mincost_flow (#2)

* added mincost func

* added default solver

* added tests for mincost, export function

* default sink and source, circulation tests

* changed LP solver, added doc, conservation tests

Author: matbesancon

REQUIRE

@@ -1,3 +1,7 @@
 julia 0.6
 LightGraphs
 SimpleTraits
+JuMP
+MathProgBase
+GLPK
+GLPKMathProgInterface

src/LightGraphsFlows.jl

@@ -5,6 +5,10 @@ const lg = LightGraphs
 using SimpleTraits: @traitfn, @traitimpl
 import SimpleTraits
 
+import JuMP
+using MathProgBase.SolverInterface: AbstractMathProgSolver
+using GLPKMathProgInterface: GLPKSolverLP
+
 import Base: getindex, size, transpose, ctranspose
 
 include("maximum_flow.jl")
@@ -15,10 +19,10 @@ include("push_relabel.jl")
 include("multiroute_flow.jl")
 include("kishimoto.jl")
 include("ext_multiroute_flow.jl")
+include("mincost.jl")
 
-# flow
 export
 maximum_flow, EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm,
-multiroute_flow, KishimotoAlgorithm, ExtendedMultirouteFlowAlgorithm
+multiroute_flow, KishimotoAlgorithm, ExtendedMultirouteFlowAlgorithm, mincost_flow
 
 end

src/mincost.jl

@@ -0,0 +1,91 @@
+"""
+    mincost_flow(graph, capacity, demand, cost[, source][, sink][, solver])
+
+Find a flow satisfying the `demand` and `capacity` constraints for each edge
+while minimizing the `sum(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.
+
+- The problem can be seen as a linear programming problem and uses a LP 
+solver under the hood. The default solver is GLPKSolverLP.
+
+Returns a flow matrix, flow[i,j] corresponds to the flow on the (i,j) arc.
+
+### Usage Example:
+
+```jldoctest
+julia> g = lg.DiGraph(6) # Create a flow-graph
+julia> lg.add_edge!(g, 5, 1)
+julia> lg.add_edge!(g, 5, 2)
+julia> lg.add_edge!(g, 3, 6)
+julia> lg.add_edge!(g, 4, 6)
+julia> lg.add_edge!(g, 1, 3)
+julia> lg.add_edge!(g, 1, 4)
+julia> lg.add_edge!(g, 2, 3)
+julia> lg.add_edge!(g, 2, 4)
+julia> w = zeros(6,6)
+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> capacity = ones(6,6)
+julia> flow = mincost_flow(g, capacity, demand, w, 5, 6)
+```
+"""
+function mincost_flow(g::lg.DiGraph, 
+		capacity::AbstractMatrix,
+		demand::AbstractMatrix,
+		cost::AbstractMatrix,
+		source::Int = -1, # if source and/or sink omitted or not in nodes, circulation problem 
+		sink::Int = -1,
+		solver::AbstractMathProgSolver = GLPKSolverLP())
+	m = JuMP.Model(solver = solver)
+	flow = JuMP.@variable(m, x[1:lg.nv(g),1:lg.nv(g)] >= 0.)
+	JuMP.@constraint(m,
+		[i=1:lg.nv(g),j=1:lg.nv(g)],
+		flow[i,j] - demand[i,j] >= 0.
+	)
+	JuMP.@constraint(m,
+		[i=1:lg.nv(g),j=1:lg.nv(g)],
+		capacity[i,j] - flow[i,j] >= 0.
+	)
+	for n1 in 1:lg.nv(g)
+		for n2 in 1:lg.nv(g)
+			if !lg.has_edge(g,n1,n2)
+				JuMP.@constraint(m, flow[n1,n2] <= 0.)
+			end
+		end
+	end
+	# flow conservation constraints
+	for node in 1:lg.nv(g)
+		if node != source && node != sink
+			JuMP.@constraint(m,
+				sum(flow[node,j] for j in 1:lg.nv(g)) -
+				sum(flow[i,node] for i in 1:lg.nv(g)) == 0
+			)
+		end
+	end
+	# source is a net flow producer
+	if source ∈ 1:lg.nv(g)
+		JuMP.@constraint(m,
+			sum(flow[source,j] for j in 1:lg.nv(g)) -
+			sum(flow[i,source] for i in 1:lg.nv(g)) >= 0
+		)
+	end
+	# source is a net flow consumer
+	if sink ∈ 1:lg.nv(g)
+		JuMP.@constraint(m,
+			sum(flow[sink,j] for j in 1:lg.nv(g)) -
+			sum(flow[i,sink] for i in 1:lg.nv(g)) <= 0
+		)
+	end
+	JuMP.@objective(m, :Min, sum(flow[i,j]*cost[i,j] for i=1:lg.nv(g),j=1:lg.nv(g)))
+	solution = JuMP.solve(m)
+	return JuMP.getvalue(flow)
+end

test/mincost.jl

@@ -0,0 +1,95 @@
+@testset "Minimum-cost flow" 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, 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, 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, -1, -1)
+    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, -1, -1)
+    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

@@ -15,7 +15,8 @@ for t in [
         "boykov_kolmogorov",
         "push_relabel",
         "maximum_flow",
-        "multiroute_flow"]
+        "multiroute_flow",
+        "mincost"]
     tp = joinpath(testdir, "$(t).jl")
     include(tp)
 end