Significant improvement to mincost_flow (#23)

* Fixed mincost_flow to have much better performace..  around 10000x faster forcases with 500 nodes.  Changed mincost_flow arguments to be more consistent with traditional problems.

* Made requested fixes to documentation, code spacing, and changed variable names to increase clarity.

* Changed argument names to add clarity

* Update testing to new interface, maintained testing for depreciated interface in seperate file.

* Fix missing test

Author: ndinsmore

src/LightGraphsFlows.jl

@@ -8,7 +8,7 @@ import SimpleTraits
 using MathProgBase.HighLevelInterface: linprog
 using MathProgBase.SolverInterface: AbstractMathProgSolver
 
-using SparseArrays: spzeros
+using SparseArrays: spzeros, sparse, sparsevec
 using Markdown: @doc_str
 
 import Base: getindex, size, transpose, adjoint
@@ -26,6 +26,7 @@ include("mincut.jl")
 
 export
 maximum_flow, EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm,
-multiroute_flow, KishimotoAlgorithm, ExtendedMultirouteFlowAlgorithm, mincost_flow, mincut
+multiroute_flow, KishimotoAlgorithm, ExtendedMultirouteFlowAlgorithm, mincost_flow
+
 
 end

src/mincost.jl

@@ -1,17 +1,22 @@
 """
-    mincost_flow(graph, capacity, demand, cost, solver, [, source][, sink])
+    mincost_flow(graph, node_demand, edge_capacity, edge_cost, solver,<keyword arguments>)
 
-Find a flow satisfying the `demand` and `capacity` constraints for each edge
-while minimizing the `sum(cost.*flow)`.
+Find a flow satisfying the `node_demand` at each node and `edge_capacity` constraints for each edge
+while minimizing the `sum(edge_cost.*flow)`.
 
-- If `source` and `sink` are specified, they are allowed a net flow production,
-consumption respectively. All other nodes must respect the flow conservation
-property.
+-`node_demand` is a vector of nodal values, which should be positive for sources and
+ negative at sink nodes, any node with zero demand is a transport node
 
-- The problem can be seen as a linear programming problem and uses a LP 
+- The problem can be seen as a linear programming problem and uses a LP
 solver under the hood. We use Clp in the examples and tests.
 
 Returns a flow matrix, flow[i,j] corresponds to the flow on the (i,j) arc.
+# Arguments
+- `edge_demand::AbstractMatrix`: demand a minimum flow for a given edge.
+- `edge_demand_exact=true`: changes the capacity of a non zero edge to the demanded value
+- `source_nodes::AbstractVector` Sources at which the nodal netflow is allowed to be greater than nodal demand **
+- `sink_nodes::AbstractVector` Sinks at which the nodal netflow is allowed to be less than nodal demand **
+   ** source_nodes & sink_nodes are only needed when nodal flow are not explictly set in node_demand
 
 ### Usage Example:
 
@@ -31,52 +36,87 @@ julia> w[1,3] = 10.
 julia> w[1,4] = 5.
 julia> w[2,3] = 2.
 julia> w[2,4] = 2.
-julia> # v2 -> sink have demand of one
-julia> demand = spzeros(6,6)
-julia> demand[3,6] = 1
-julia> demand[4,6] = 1
+julia> demand = spzeros(6)
+julia> demand[5] = 2
+julia> demand[6] = -2
 julia> capacity = ones(6,6)
-julia> flow = mincost_flow(g, capacity, demand, w, ClpSolver(), 5, 6)
+julia> flow = mincost_flow(g, demand, capacity, w, ClpSolver())
 ```
 """
+
+
 function mincost_flow(g::lg.DiGraph,
-		capacity::AbstractMatrix,
-		demand::AbstractMatrix,
-		cost::AbstractMatrix,
-		solver::AbstractMathProgSolver,
-		source::Int = -1, # if source and/or sink omitted or not in nodes, circulation problem
-		sink::Int = -1)
-	flat_cap = collect(Iterators.flatten(capacity))
-	flat_dem = collect(Iterators.flatten(demand))
-	flat_cost = collect(Iterators.flatten(cost))
-	n = lg.nv(g)
-	b = zeros(n+n*n)
-	A = spzeros(n+n*n,n*n)
-	sense = ['=' for _ in 1:n]
-	append!(sense, ['>' for _ in 1:n*n])
-	for node in 1:n
-		if sink == node
-			sense[sink] = '>'
-		elseif source == node
-			sense[source] = '<'
-		end
-		col_idx = (node-1)*n+1:node*n
-		line_idx = node:n:n*n
-		for jdx in col_idx
-			A[node,jdx] = A[node,jdx]+1.0
-		end
-		for idx in line_idx
-			A[node,idx] = A[node,idx]-1.0
+		node_demand::AbstractVector,
+		edge_capacity::AbstractMatrix,
+		edge_cost::AbstractMatrix,
+		solver::AbstractMathProgSolver;
+		edge_demand::Union{Nothing,AbstractMatrix} = nothing,
+		edge_demand_exact::Bool = false, #If true changes capacity of that node to the value in edge demand
+		source_nodes::AbstractVector{<:Integer} = Vector{Integer}(), #Source nodes at which to allow a netflow greater than nodal demand
+		sink_nodes::AbstractVector{<:Integer} = Vector{Integer}()	 #Sink nodes at which to allow a netflow less than nodal demand
+		)
+	T = eltype(g)
+	VT = promote_type(eltype(edge_capacity),eltype(node_demand),eltype(edge_cost))
+
+	nv = lg.nv(g)
+	ne = lg.ne(g)
+	b = Vector{VT}(undef,nv)
+	AI = Vector{T}(undef,2*ne)
+	AJ = Vector{T}(undef,2*ne)
+	AV = Vector{VT}(undef,2*ne)
+	edge_costs = Vector{VT}(undef,ne)
+	upperbounds = Vector{VT}(undef,ne)
+	lowerbounds = Vector{VT}(undef,ne)
+	sense = ['=' for _ in 1:nv]
+
+	has_edge_demand= edge_demand !== nothing
+
+	inds = [-1,0]
+	for (n,e) in enumerate(lg.edges(g))
+		inds .+= 2
+		s = lg.src(e)
+		t = lg.dst(e)
+		AI[inds] = [s,t] #the node index
+		AJ[inds] = [n,n]  #the edge index
+		AV[inds] = [1,-1]
+		upperbounds[n] = edge_capacity[s,t]
+		lowerbounds[n] = has_edge_demand ? edge_demand[s,t] : zero(VT)
+		if edge_demand_exact && has_edge_demand && lowerbounds[n] != 0
+			upperbounds[n] = lowerbounds[n]
 		end
+
+		edge_costs[n] = edge_cost[s,t]
 	end
-	for src in 1:n
-		for dst in 1:n
-			if lg.Edge(src,dst) ∉ lg.edges(g)
-				A[n+src+n*(dst-1),src+n*(dst-1)] = 1
-				sense[n+src+n*(dst-1)] = '<'
-			end
-		end
+	for n = 1:nv
+		b[n] = node_demand[n]
 	end
-	sol = linprog(flat_cost, A, sense, b, flat_dem, flat_cap, solver)
-	[sol.sol[idx + n*(jdx-1)] for idx=1:n,jdx=1:n]
+
+	for n in source_nodes
+		sense[n] = '>'
+	end
+	for n in sink_nodes
+		sense[n] = '<'
+	end
+
+
+	A = sparse(AI,AJ,AV,nv,ne)
+	# return (edge_costs, A, sense, b, lowerbounds, upperbounds, solver)
+	sol = linprog(edge_costs, A, sense, b, lowerbounds, upperbounds, solver)
+	sol_sparse = sparse(view(AI,1:2:2*ne),view(AI,2:2:2*ne),sol.sol,nv,nv)
+	return sol_sparse
 end
+
+@deprecate mincost_flow(g::lg.DiGraph,
+		capacity::AbstractMatrix,
+		demand::AbstractMatrix,
+		cost::AbstractMatrix,
+		solver::AbstractMathProgSolver,
+		source::Int, # if source and/or sink omitted or not in nodes, circulation problem
+		sink::Int)  mincost_flow(g,spzeros(lg.nv(g)),capacity,cost,solver,edge_demand=demand,source_nodes=[source],sink_nodes=[sink])
+
+@deprecate mincost_flow(g::lg.DiGraph,
+		capacity::AbstractMatrix,
+		demand::AbstractMatrix,
+		cost::AbstractMatrix,
+		solver::AbstractMathProgSolver
+		)  mincost_flow(g,spzeros(lg.nv(g)),capacity,cost,solver,edge_demand=demand)

test/mincost.jl

@@ -25,7 +25,7 @@ using SparseArrays: spzeros
     demand[4,6] = 1
     capacity = ones(6,6)
 
-    flow = mincost_flow(g, capacity, demand, w, ClpSolver(), 5, 6)
+    flow = mincost_flow(g, spzeros(lg.nv(g)),capacity, w, ClpSolver(), edge_demand=demand,source_nodes=[5],sink_nodes= [6])
     @test flow[5,1] == 1
     @test flow[5,2] == 1
     @test flow[3,6] == 1
@@ -41,9 +41,66 @@ using SparseArrays: spzeros
         @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
     end
 
+    #same test with exact flows
+    flow = mincost_flow(g, spzeros(lg.nv(g)),capacity, w, ClpSolver(),edge_demand_exact=true, edge_demand=demand,source_nodes=[5],sink_nodes= [6])
+    @test flow[5,1] == 1
+    @test flow[5,2] == 1
+    @test flow[3,6] == 1
+    @test flow[4,6] == 1
+    @test flow[1,4] == 1
+    @test flow[2,4] == 0
+    @test flow[2,3] == 1
+    @test flow[1,3] == 0
+    @test sum(diag(flow)) == 0
+
+    # flow conservation property
+    for n1 = 1:4
+        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
+    end
+
+    #--- Same test without edge_demands & positive cost should result in null flow
+    flow = mincost_flow(g, spzeros(lg.nv(g)),capacity, w, ClpSolver(),source_nodes=[5],sink_nodes= [6])
+    @test flow[5,1] == 0
+    @test flow[5,2] == 0
+    @test flow[3,6] == 0
+    @test flow[4,6] == 0
+    @test flow[1,4] == 0
+    @test flow[2,4] == 0
+    @test flow[2,3] == 0
+    @test flow[1,3] == 0
+    @test sum(diag(flow)) == 0
+
+    # flow conservation property
+    for n1 = 1:4
+        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
+    end
+
+
+    #run sam test again with nodal demands
+    node_demand=spzeros(6)
+    node_demand[5]=2
+    node_demand[6]=-2
+    #--- Same test without edge_demands
+    flow = mincost_flow(g, node_demand,capacity, w, ClpSolver())
+    @test flow[5,1] == 1
+    @test flow[5,2] == 1
+    @test flow[3,6] == 1
+    @test flow[4,6] == 1
+    @test flow[1,4] == 1
+    @test flow[2,4] == 0
+    @test flow[2,3] == 1
+    @test flow[1,3] == 0
+    @test sum(diag(flow)) == 0
+
+    # flow conservation property
+    for n1 = 1:4
+        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
+    end
+
+
     # no demand => null flow
     d2 = spzeros(6,6)
-    flow = mincost_flow(g, capacity, d2, w, ClpSolver(), 5, 6)
+    flow = mincost_flow(g, spzeros(lg.nv(g)),capacity, w, ClpSolver(), edge_demand=d2,source_nodes=[5],sink_nodes= [6])
     for idx in 1:6
         for jdx in 1:6
             @test flow[idx,jdx] ≈ 0.0
@@ -63,7 +120,7 @@ using SparseArrays: spzeros
     demand = spzeros(6,6)
     demand[1,2] = 1
     costs = ones(6,6)
-    flow = mincost_flow(g, capacity, demand, costs, ClpSolver())
+    flow = mincost_flow(g,spzeros(lg.nv(g)), capacity, costs, ClpSolver(), edge_demand=demand)
     active_flows = [(1,2), (2,5), (5,6),(6,1)]
     for s in 1:6
         for t in 1:6
@@ -80,7 +137,7 @@ using SparseArrays: spzeros
     end
     # higher short-circuit cost
     costs[2,5] = 10.
-    flow = mincost_flow(g, capacity, demand, costs, ClpSolver())
+    flow = mincost_flow(g,spzeros(lg.nv(g)), capacity, costs, ClpSolver(), edge_demand=demand)
     active_flows = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,1)]
     for s in 1:6
         for t in 1:6

test/mincost_dep.jl

@@ -0,0 +1,98 @@
+using Clp: ClpSolver
+using SparseArrays: spzeros
+#This is for testing of the depreciated interface.   Please remove once that interface is removed
+@testset "Minimum-cost flow Depreciated" begin
+
+    # bipartite oriented + source & sink
+    # source 5, sink 6, v1 {1, 2} v2 {3, 4}
+    g = lg.DiGraph(6)
+    lg.add_edge!(g, 5, 1)
+    lg.add_edge!(g, 5, 2)
+    lg.add_edge!(g, 3, 6)
+    lg.add_edge!(g, 4, 6)
+    lg.add_edge!(g, 1, 3)
+    lg.add_edge!(g, 1, 4)
+    lg.add_edge!(g, 2, 3)
+    lg.add_edge!(g, 2, 4)
+    w = zeros(6,6)
+    w[1,3] = 10.
+    w[1,4] = 5.
+    w[2,3] = 2.
+    w[2,4] = 2.
+    # v2 -> sink have demand of one
+    demand = spzeros(6,6)
+    demand[3,6] = 1
+    demand[4,6] = 1
+    capacity = ones(6,6)
+
+    flow = mincost_flow(g, capacity, demand, w, ClpSolver(), 5, 6)
+    @test flow[5,1] == 1
+    @test flow[5,2] == 1
+    @test flow[3,6] == 1
+    @test flow[4,6] == 1
+    @test flow[1,4] == 1
+    @test flow[2,4] == 0
+    @test flow[2,3] == 1
+    @test flow[1,3] == 0
+    @test sum(diag(flow)) == 0
+
+    # flow conservation property
+    for n1 = 1:4
+        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
+    end
+
+    # no demand => null flow
+    d2 = spzeros(6,6)
+    flow = mincost_flow(g, capacity, d2, w, ClpSolver(), 5, 6)
+    for idx in 1:6
+        for jdx in 1:6
+            @test flow[idx,jdx] ≈ 0.0
+        end
+    end
+
+    # graph without sink or source
+    g = lg.DiGraph(6)
+    # create circuit
+    for s in 1:5
+        lg.add_edge!(g, s, s+1)
+    end
+    lg.add_edge!(g, 6, 1)
+    lg.add_edge!(g, 2, 5) # shortcut
+
+    capacity = ones(6,6)
+    demand = spzeros(6,6)
+    demand[1,2] = 1
+    costs = ones(6,6)
+    flow = mincost_flow(g, capacity, demand, costs, ClpSolver())
+    active_flows = [(1,2), (2,5), (5,6),(6,1)]
+    for s in 1:6
+        for t in 1:6
+            if (s,t) in active_flows
+                @test flow[s,t] ≈ 1.
+            else
+                @test flow[s,t] ≈ 0.
+            end
+        end
+    end
+    # flow conservation property
+    for n1 = 1:6
+        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
+    end
+    # higher short-circuit cost
+    costs[2,5] = 10.
+    flow = mincost_flow(g, capacity, demand, costs, ClpSolver())
+    active_flows = [(1,2),(2,3),(3,4),(4,5),(5,6),(6,1)]
+    for s in 1:6
+        for t in 1:6
+            if (s,t) in active_flows
+                @test flow[s,t] ≈ 1.
+            else
+                @test flow[s,t] ≈ 0.
+            end
+        end
+    end
+        # flow conservation property
+    for n1 = 1:6
+        @test sum(flow[n1,:]) ≈ sum(flow[:,n1])
+    end
+end

test/runtests.jl

@@ -19,6 +19,7 @@ for t in [
         "maximum_flow",
         "multiroute_flow",
         "mincost",
+        "mincost_dep",  #Please remove this once depreciated mincost_flow iterface is removed
         "mincut",
         ]
     tp = joinpath(testdir, "$(t).jl")