Mathprog (#5)

* 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

* added mincost docs

* replace JuMP with basic MathProg

* removed JuMP from require

Author: matbesancon

REQUIRE

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

docs/make.jl

@@ -11,5 +11,6 @@ makedocs(
 		"Getting started"    => "index.md",
         "Maxflow algorithms" => "maxflow.md",
         "Multiroute flows"   => "multiroute.md",
+        "Min-cost flows"   => "mincost.md",
     ]
 )

docs/src/mincost.md

@@ -0,0 +1,7 @@
+## Min-cost flow
+
+```@autodocs
+Modules = [LightGraphsFlows]
+Pages = ["mincost.jl"]
+Order = [:function, :type]
+```

src/LightGraphsFlows.jl

@@ -5,7 +5,7 @@ const lg = LightGraphs
 using SimpleTraits: @traitfn, @traitimpl
 import SimpleTraits
 
-import JuMP
+using MathProgBase.HighLevelInterface: linprog
 using MathProgBase.SolverInterface: AbstractMathProgSolver
 using GLPKMathProgInterface: GLPKSolverLP
 

src/mincost.jl

@@ -15,7 +15,7 @@ Returns a flow matrix, flow[i,j] corresponds to the flow on the (i,j) arc.
 
 ### Usage Example:
 
-```jldoctest
+```julia
 julia> g = lg.DiGraph(6) # Create a flow-graph
 julia> lg.add_edge!(g, 5, 1)
 julia> lg.add_edge!(g, 5, 2)
@@ -45,47 +45,37 @@ function mincost_flow(g::lg.DiGraph,
 		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
+	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
-	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
-			)
+		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
 		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
-		)
+	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
 	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)
+	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]
 end

test/boykov_kolmogorov.jl

@@ -10,8 +10,8 @@
 
     for g in testdigraphs(gg)
     # default capacity
-      capacity_matrix = lg.DefaultCapacity(g)
-      residual_graph = @inferred(lg.residual(g))
+      capacity_matrix = LightGraphsFlows.DefaultCapacity(g)
+      residual_graph = @inferred(LightGraphsFlows.residual(g))
       T = eltype(g)
       flow_matrix = zeros(T, 3, 3)
       TREE = zeros(T, 3)
@@ -21,7 +21,7 @@
       A = [T(source), T(target)]
 # see https://github.com/JuliaLang/julia/issues/21077
 # @show("testing $g with eltype $T, residual_graph type is $(eltype(residual_graph)), flow_matrix type is $(eltype(flow_matrix)), capacity_matrix type is $(eltype(capacity_matrix))")
-      path = lg.find_path!(
+      path = LightGraphsFlows.find_path!(
         residual_graph, source, target, flow_matrix,
         capacity_matrix, PARENT, TREE, A)
 

test/dinic.jl

@@ -19,7 +19,7 @@
 
     # Construct the residual graph
     for fg in (flow_graph, lg.DiGraph{UInt8}(flow_graph), lg.DiGraph{Int16}(flow_graph))
-      residual_graph = @inferred(lg.residual(fg))
+      residual_graph = @inferred(LightGraphsFlows.residual(fg))
 
       # Test with default distances
       @test @inferred(LightGraphsFlows.dinic_impl(residual_graph, 1, 8, LightGraphsFlows.DefaultCapacity(residual_graph)))[1] == 3