Finset algs

Project.toml

@@ -6,9 +6,16 @@ version = "0.0.1"
 
 [deps]
 AlgebraicDynamics = "5fd6ff03-a254-427e-8840-ba658f502e32"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-
-[compat]
-Catlab = "^0.15"
-julia = "1.9"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/AlgebraicOptimization.jl

@@ -2,6 +2,11 @@
 """
 module AlgebraicOptimization
 
-include("OpenProblems.jl")
+include("FinSetAlgebras.jl")
+include("Optimizers.jl")
+include("Objectives.jl")
+include("OpenFlowGraphs.jl")
+#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/FinSetAlgebras.jl

@@ -0,0 +1,214 @@
+module FinSetAlgebras
+
+export FinSetAlgebra, CospanAlgebra, Open, hom_map, laxator, data, portmap
+
+using LinearAlgebra, SparseArrays
+using Catlab
+import Catlab: oapply, dom, Cospan
+
+#abstract type AlgebraObject end
+
+abstract type FinSetAlgebra{T} end # A finset algebra with object type T
+
+#=function dom(X::T)::FinSet where T <: AlgebraObject
+    error("Domain not specified.")
+end=#
+
+#=function ob_map(::FinSetAlgebra{T}, N::FinSet)::T where T <: AlgebraObject
+    error("Object map not implemented.")
+end=#
+
+function hom_map(::FinSetAlgebra{T}, ϕ::FinFunction, X::T)::T where T
+    error("Morphism map not implemented.")
+end
+
+function laxator(::FinSetAlgebra{T}, Xs::Vector{T})::T where T
+    error("Laxator not implemented.")
+end
+
+function oapply(A::FinSetAlgebra{T}, ϕ::FinFunction, Xs::Vector{T})::T where T
+    return hom_map(A, ϕ, laxator(A, Xs))
+end
+
+# Example
+function pullback_matrix(f::FinFunction)
+    n = length(dom(f))
+    sparse(1:n, f.(dom(f)), ones(Int,n), dom(f).n, codom(f).n)
+end
+
+pushforward_matrix(f::FinFunction) = pullback_matrix(f)'
+
+#=struct FreeVector <: AlgebraObject
+    dim::FinSet
+    v::Vector{Float64}
+end=#
+
+#dom(v::FreeVector) = v.dim
+
+struct Pushforward <: FinSetAlgebra{Vector{Float64}} end
+
+dom(v::Vector{Float64}) = FinSet(length(v))
+
+hom_map(::Pushforward, ϕ::FinFunction, v::Vector{Float64}) = pushforward_matrix(ϕ)*v
+
+
+laxator(::Pushforward, Xs::Vector{Vector{Float64}}) = vcat(Xs...)
+
+# UWD algebras from finset algebras
+abstract type CospanAlgebra{T} end
+#abstract type SimpleCospanAlgebra{T, A<:FinSetAlgebra{T}} <: CospanAlgebra{T} end
+
+function hom_map(::CospanAlgebra{T}, ϕ::Cospan, X::T)::T where T
+    error("Morphism map not implemented.")
+end
+
+function laxator(::CospanAlgebra{T}, Xs::Vector{T})::T where T
+    error("Laxator not implemented.")
+end
+
+
+function oapply(A::CospanAlgebra{T}, ϕ::Cospan, Xs::Vector{T})::T where T
+    return hom_map(A, ϕ, laxator(A, Xs))
+end
+
+
+struct Open{T}
+    S::FinSet
+    o::T
+    m::FinFunction    
+    Open{T}(S, o, m) where T = 
+        S != codom(m) || dom(o) != S ? error("Invalid portmap.") : new(S, o, m)
+end
+
+data(obj::Open{T}) where T = obj.o
+portmap(obj::Open{T}) where T = obj.m
+
+
+function Open{T}(o::T) where T
+    Open{T}(domain(o), o, id(domain(o)))
+end
+
+dom(obj::Open{T}) where T = dom(obj.m)
+
+function hom_map(::CospanAlgebra{Open{T}}, A::FinSetAlgebra{T}, ϕ::Cospan, X::Open{T})::Open{T} where T
+    l = left(ϕ)
+    r = right(ϕ)
+    p = pushout(X.m, l)
+    pL = legs(p)[1]
+    pR = legs(p)[2]
+    return Open{T}(apex(p), hom_map(A, pL, X.o), compose(r,pR))
+end
+
+function laxator(::CospanAlgebra{Open{T}}, A::FinSetAlgebra{T}, Xs::Vector{Open{T}})::Open{T} where T
+    S = coproduct([X.S for X in Xs])
+    inclusions(i::Int) = legs(S)[i]
+    m = copair([compose(Xs[i].m, inclusions(i)) for i in 1:length(Xs)])
+    o = laxator(A, [X.o for X in Xs])
+    return Open{T}(apex(S), o, m)
+end
+
+function oapply(CA::CospanAlgebra{Open{T}}, FA::FinSetAlgebra{T}, ϕ::Cospan, Xs::Vector{Open{T}})::Open{T} where T
+    return hom_map(CA, FA, ϕ, laxator(CA, FA, Xs))
+end
+
+function uwd_to_cospan(d::AbstractUWD)
+    # Build the left leg
+    left_dom = vcat([length(ports(d, i)) for i in boxes(d)])
+    left_codom = njunctions(d)
+
+    #println(cp_dom)
+    ports_to_junctions = FinFunction[]
+    total_portmap = subpart(d, :junction)
+
+    for box in ports.([d], boxes(d))
+        push!(ports_to_junctions, FinFunction([total_portmap[p] for p in box], length(box), left_codom))
+    end
+    #println(ports_to_junctions)
+    #cp = CompositionPattern(cp_dom, cp_codom, ports_to_junctions)
+
+    left = copair(ports_to_junctions)
+    right = FinFunction(subpart(d, :outer_junction), left_codom)
+
+    return Cospan(left, right)  
+end
+
+function oapply(CA::CospanAlgebra{Open{T}}, FA::FinSetAlgebra{T}, d::AbstractUWD, Xs::Vector{Open{T}})::Open{T} where T
+    return oapply(CA, FA, uwd_to_cospan(d), Xs)
+end
+
+
+end
+#=
+# Test example
+struct UWDPushforward <: CospanAlgebra{Open{Vector{Float64}}} end
+
+const OpenVector = Open{Vector{Float64}}
+
+# UWD Interop
+function uwd_to_cospan(d::AbstractUWD)
+    # Build the left leg
+    left_dom = vcat([length(ports(d, i)) for i in boxes(d)])
+    left_codom = njunctions(d)
+
+    #println(cp_dom)
+    ports_to_junctions = FinFunction[]
+    total_portmap = subpart(d, :junction)
+
+    for box in ports.([d], boxes(d))
+        push!(ports_to_junctions, FinFunction([total_portmap[p] for p in box], length(box), left_codom))
+    end
+    #println(ports_to_junctions)
+    #cp = CompositionPattern(cp_dom, cp_codom, ports_to_junctions)
+
+    left = copair(ports_to_junctions)
+    right = FinFunction(subpart(d, :outer_junction), left_codom)
+
+    return Cospan(left, right)  
+end
+
+function oapply(CA::CospanAlgebra{Open{T}}, FA::FinSetAlgebra{T}, d::AbstractUWD, Xs::Vector{Open{T}})::Open{T} where T
+    return oapply(CA, FA, uwd_to_cospan(d), Xs)
+end
+
+# AlgebraicDynamics
+struct System
+    state_space::FinSet
+    dynamics::Function # R^ss -> R^ss
+end
+(s::System)(x::Vector) = s.dynamics(x)
+dom(s::System) = s.state_space
+
+struct Dynam <: FinSetAlgebra{System} end
+
+hom_map(::Dynam, ϕ::FinFunction, s::System) = 
+    System(codom(ϕ), x->pushforward_matrix(ϕ)*s(pullback_matrix(ϕ)*x))
+
+function laxator(::Dynam, Xs::Vector{System})
+    c = coproduct([dom(X) for X in Xs])
+    subsystems = [x -> X(pullback_matrix(l)*x) for (X,l) in zip(Xs, legs(c))]
+    function parallel_dynamics(x)
+        res = Vector{Vector}(undef, length(Xs)) # Initialize storage for results
+        Threads.@threads for i = 1:length(Xs)
+            res[i] = subsystems[i](x)
+        end
+        return vcat(res...)
+    end
+    return System(apex(c), parallel_dynamics)
+end
+
+struct UWDDynam <: CospanAlgebra{Open{System}} end # Turn Dynam into a UWD algebra in one line!
+
+A = rand(-1.0:.01:1.0, 5,5)
+B = rand(-1.0:.01:1.0, 3,3)
+C = rand(-1.0:.01:1.0, 4,4)
+
+γ = 0.1
+s1 = System(FinSet(5), x->x +γ*A*x)
+s2 = System(FinSet(3), x->x + γ*B*x)
+s3 = System(FinSet(4), x->x + γ*C*x)
+
+ϕ = FinFunction([1,2,3,4,5,2,3,6,3,6,7,8])
+
+s = oapply(Dynam(), ϕ, [s1,s2,s3])=#
+
+#end