Netflow

Project.toml

@@ -6,8 +6,14 @@ version = "0.0.1"
 
 [deps]
 AlgebraicDynamics = "5fd6ff03-a254-427e-8840-ba658f502e32"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 
 [compat]
 Catlab = "^0.15"

src/AlgebraicOptimization.jl

@@ -3,5 +3,6 @@
 module AlgebraicOptimization
 
 include("OpenProblems.jl")
+include("FlowGraphs.jl")
 
 end

src/Experiments.jl

@@ -0,0 +1,144 @@
+using NLsolve
+using LinearAlgebra
+using ForwardDiff
+using Catlab
+
+function node_incidence_matrix(g::Graph)
+    V = nv(g)
+    E = ne(g)
+    A = zeros(V,E)
+    for (v,e) in Iterators.product(1:V, 1:E)
+        if src(g, e) == tgt(g, e) && tgt(g, e) == v
+            continue
+        elseif src(g,e) == v
+            A[v,e] = 1
+        elseif tgt(g,e) == v
+            A[v,e] = -1
+        end
+    end
+    return A
+end
+
+N_subprob = 10
+A1 = node_incidence_matrix(wheel_graph(Graph, N_subprob))
+A2 = A1
+A3 = A1
+
+
+
+
+
+
+
+function draw2(n::Int)
+    a = rand(1:n)
+    b = rand(1:n-1)
+    b += (b >= a)
+    return a, b
+end
+
+function random_column(length::Int)
+    a,b = draw2(length)
+    res = zeros(length)
+    res[a] = 1
+    res[b] = -1
+    return res
+end
+
+E = 50
+V = 30
+num_sources = 5
+num_sinks = 5
+
+A = hcat([random_column(V) for i in 1:E]...)
+
+f(x) = x^2
+
+b = zeros(V)
+
+for i in 1:num_sources
+    src_vertex = rand(1:V)
+    val = rand(1:0.01:10)
+    b[src_vertex] = val
+end
+
+for i in 1:num_sinks
+    sink_vertex = rand(1:V)
+    val = rand(1:0.01:10)
+    b[sink_vertex] = -val
+end
+λ = randn(V)
+total_L(x) = sum([f(x_i) for x_i in x]) + λ'*(A*x-b) 
+
+L(i) = x -> f(x) + λ'*(A[:,i]*x - b)
+
+function grad_flow_L(x::Vector)
+    ForwardDiff.gradient(total_L, x)
+end
+
+function grad_flow_L(i::Int)
+    x -> ForwardDiff.derivative(L(i), x)
+end
+
+function grad_descent(f, x0, γ, max_iters, ϵ)
+    x_prev = x0
+    x_cur = x0
+    for i in 1:max_iters
+        x_cur = x_prev - γ*ForwardDiff.gradient(f, x_prev)
+        if norm(f(x_cur) - f(x_prev)) < ϵ
+            #println("Terminated in $i iterations.")
+            return x_cur
+        end
+        x_prev = x_cur
+    end
+    println("Did not converge.")
+    return x_cur
+end
+
+function grad_descent_1D(f, x0, γ, max_iters, ϵ)
+    x_prev = x0
+    x_cur = x0
+    for i in 1:max_iters
+        x_cur = x_prev - γ*ForwardDiff.derivative(f, x_prev)
+        if norm(f(x_cur) - f(x_prev)) < ϵ
+            #println("Terminated in $i iterations.")
+            return x_cur
+        end
+        x_prev = x_cur
+    end
+    println("Did not converge.")
+    return x_cur
+end
+
+sol_nl_total = nlsolve(grad_flow_L, repeat([10.0], E), iterations=1000000, xtol=0.01).zero
+@time nlsolve(grad_flow_L, repeat([10.0], E), iterations=1000000, xtol=0.01)
+
+function sol_nl(i::Int) 
+    f = grad_flow_L(i)
+    sol = nlsolve(n_ary(f), [10.0], xtol=0.01)
+    return sol.zero[1]
+end
+
+sol_nl_distributed = zeros(E)
+for i in 1:E
+    sol_nl_distributed[i]=sol_nl(i)
+end
+
+@time for i in 1:E
+    sol_nl(i)
+end
+
+sol_total = grad_descent(total_L, repeat([10.0], E), 0.1,100000, 0.0001)
+@time grad_descent(total_L, repeat([10.0], E), 0.1,100000, 0.01)
+sol(i) = grad_descent_1D(L(i), 10.0, 0.1, 100000, 0.01)
+
+sol_distributed = zeros(E)
+
+for i in 1:E
+    sol_distributed[i] = sol(i)
+end
+
+@time for i in 1:E
+    sol(i)
+end
+

src/FlowGraphs.jl

@@ -0,0 +1,107 @@
+""" Define the UWD algebra of open flow graphs
+"""
+module FlowGraphs
+
+export FlowGraph, to_problem, node_incidence_matrix
+
+using Catlab
+import Catlab: src, tgt
+
+using ..OpenProblems
+
+struct FlowGraph
+    exposed_nodes::Int
+    nodes::Int
+    edges::Int
+    src::FinFunction # E -> V
+    tgt::FinFunction # E -> V
+    edge_costs::Vector{Function}
+    flows::Vector{Float64}
+    p::FinFunction
+    FlowGraph(ens, ns, es, src, tgt, e_cs, fs, p) = 
+        dom(p) != FinSet(ens) || codom(p) != FinSet(ns) ?
+            error("Invalid portmap") : new(ens, ns, es, src, tgt, e_cs, fs, p)
+end
+
+# Construct a flow graph from a Catlab graph
+function FlowGraph(g::Graph, costs, flows, portmap)
+    V = nv(g)
+    exposed_V = length(dom(portmap))
+    return FlowGraph(exposed_V, V, ne(g), FinFunction(src(g)), FinFunction(tgt(g)), costs, flows, portmap)
+end
+
+# Convenience constructor when `nodes`==`exposed_nodes` and `p` is the identity function
+FlowGraph(nodes, src, tgt, edge_costs, flows) =
+    FlowGraph(nodes, nodes, src, tgt, edge_costs, flows, FinFunction(1:nodes))
+
+nvertices(g::FlowGraph) = g.nodes
+n_exposed_vertices(g::FlowGraph) = g.exposed_nodes
+n_edges(g::FlowGraph) = g.edges
+portmap(g::FlowGraph) = g.p
+src(g::FlowGraph, e::Int) = g.src(e)
+tgt(g::FlowGraph, e::Int) = g.tgt(e)
+
+function node_incidence_matrix(g::FlowGraph)
+    V = nvertices(g)
+    E = n_edges(g)
+    A = zeros(V,E)
+    for (v,e) in Iterators.product(1:V, 1:E)
+        if src(g, e) == tgt(g, e) && tgt(g, e) == v
+            continue
+        elseif src(g,e) == v
+            A[v,e] = 1
+        elseif tgt(g,e) == v
+            A[v,e] = -1
+        end
+    end
+    return A
+end
+
+function to_problem(g::FlowGraph)
+    A = node_incidence_matrix(g)
+    obj = (x,λ) -> sum([g.edge_costs[i](x[i]) for i in 1:n_edges(g)]) + λ'*(A*x - g.flows)
+
+    return DualMaxProblem(
+        n_exposed_vertices(g),
+        nvertices(g),
+        n_edges(g),
+        obj,
+        portmap(g)
+    )
+end
+
+fills(d::AbstractUWD, b::Int, g::FlowGraph) = 
+    b <= nparts(d, :Box) ? length(incident(d, b, :box)) == n_exposed_vertices(g) : 
+        error("Trying to fill box $b when $d has fewer than $b boxes")
+
+### oapply helper functions
+
+induced_ports(d::AbstractUWD) = nparts(d, :OuterPort)
+induced_ports(d::RelationDiagram) = subpart(d, [:outer_junction, :variable])
+
+# Returns the pushout which induces the new set of vertices
+function induced_vertices(d::AbstractUWD, gs::Vector{FlowGraph}, inclusions::Function)
+    for b in parts(d, :Box)
+        fills(d, b, gs[b]) || error("$(gs[b]) does not fill box $b")
+    end
+
+    total_portmap = copair([compose(portmap(gs[i]), inclusions(i)) for i in 1:length(gs)])
+
+    #return pushout(FinFunction(subpart(d, :junction), nparts(d, :Junction)), total_portmap)
+    return pushout(total_portmap, FinFunction(subpart(d, :junction), nparts(d, :Junction)))
+end
+
+function induced_graph(d::AbstractUWD, gs::Vector{FlowGraph}, var_map::FinFunction, inclusions::Function)
+    #proj_mats = Matrix[]
+    for b in parts(d, :Box)
+        inc = compose(inclusions(b), var_map)
+        push!(proj_mats, induced_matrix(inc))
+    end
+    return z::Vector -> sum([ps[b](proj_mats[b]*z) for b in 1:length(ps)])
+end
+
+
+
+
+
+end # module

src/FreeVectorSpaces.jl

@@ -0,0 +1,7 @@
+module FreeVectorSpaces
+
+export pushforward, pullback
+
+
+
+end # module