diff --git a/src/Debug.jl b/src/Debug.jl deleted file mode 100644 index 6c8299f..0000000 --- a/src/Debug.jl +++ /dev/null @@ -1,127 +0,0 @@ -module Debug - using CombinatorialSpaces - using CombinatorialSpaces: volume - using Catlab.Graphics - using Catlab.WiringDiagrams - using Catlab.CategoricalAlgebra - using Catlab.Theories - using GeometryBasics - using ...CairoMakie - export sim_key, get_wire, draw_wire - - function sim_key(dwd; kw...) - n_dwd = deepcopy(dwd) - outer_ports = nparts(n_dwd.diagram, :OuterInPort) - n_dwd.diagram[:outer_in_port_type] .= "inp" - n_dwd.diagram[:out_port_type] .= (1:nparts(n_dwd.diagram, :OutPort)) - to_graphviz(n_dwd; labels=true, kw...) - end - - function get_wire(dwd, vf, u0, wire; p=[], t=0) - vf(copy(u0), u0, p, t) - return vf.mem.contents[wire] - end - - function draw_wire(s, sd, dwd, vf, u0, wire; p=[], t=0, colorrange=nothing, xlim=nothing, ylim=nothing, axisaspect=nothing, kw...) - w_data = get_wire(dwd, vf, u0, wire; p=p, t=t) - w_type = supertype(typeof(dwd.diagram[wire, :out_port_type][:type])) - plot(w_type, s, sd, w_data, colorrange, xlim, ylim, axisaspect; kw...) - end - - function plot(::Type{ObExpr{:Form0}}, s, sd, vals, colorrange, xlim, ylim, axisaspect; kw...) - fig = Figure() - ax, ob = mesh(fig[1,1], s, color=vals) - if !isnothing(colorrange) - ob.colorrange = colorrange - end - if !isnothing(xlim) - xlims!(ax, xlim) - end - if !isnothing(ylim) - ylims!(ax, ylim) - end - if !isnothing(axisaspect) - ax.aspect = AxisAspect(axisaspect) - end - Colorbar(fig[1, 2], ob) - fig, ax, ob - end - - function plot(::Type{ObExpr{:DualForm2}}, s, sd, vals, colorrange, xlim, ylim, axisaspect; kw...) - vals ./= [sum(dual_volume(Val{2}, sd, elementary_duals(Val{0},sd,v))) for v in 1:nv(s)] - plot(ObExpr{:Form0}, s, sd, vals, colorrange, xlim, ylim, axisaspect; kw...) - end - - function plot(::Type{ObExpr{:DualForm1}}, s, sd, vals, colorrange, xlim, ylim, axisaspect; kw...) - vals .*= [volume(Val{1},sd,e) / sum(dual_volume(Val{1}, sd, elementary_duals(Val{1},sd,e))) - for e in 1:ne(s)] - plot(ObExpr{:Form1}, s, sd, vals, colorrange, xlim, ylim, axisaspect; kw...) - end - - function plot(::Type{ObExpr{:DualForm0}}, s, sd, vals, colorrange, xlim, ylim, axisaspect; kw...) - vals .*= [volume(Val{2},sd,t) for t in 1:ntriangles(s)] - plot(ObExpr{:Form2}, s, sd, vals, colorrange, xlim, ylim, axisaspect; kw...) - end - - function plot(::Type{ObExpr{:Form1}}, s, sd, vals, colorrange, xlim, ylim, axisaspect; use_arrows=false, n_arrows = 100, kw...) - vals ./= [volume(Val{1},sd,e) for e in 1:ne(s)] - if isnothing(colorrange) - colorrange = (0.0, maximum(abs.(vals))) - end - fig = Figure() - ax, ob = if use_arrows - orient_vals = vals .* [eo ? 1 : -1 for eo in s[:edge_orientation]] - signs = sign.(orient_vals) - locs = s[s[:∂v1], :point] .* (0.5 .- (signs * 0.35)) .+ s[s[:∂v0], :point] .* (0.5 .+ (signs * 0.35)) - mag = ((s[s[:∂v1], :point] .- s[s[:∂v0], :point]) * 0.5) .* signs - inds = collect(1:ne(s)) - if !isnothing(xlim) - filter!(i -> xlim[1] <= s[s[i, :∂v0], :point][1] <= xlim[2], inds) - end - if !isnothing(ylim) - filter!(i -> ylim[1] <= s[s[i, :∂v0], :point][2] <= ylim[2], inds) - end - inds = inds[1:min(n_arrows, length(inds))] - - arrows(fig[1,1], locs[inds], mag[inds], color=abs.(vals[inds]), colorrange=colorrange) - else - linesegments(fig[1,1], s, color=abs.(vals), colorrange=colorrange) - end - if !isnothing(xlim) - xlims!(ax, xlim) - end - if !isnothing(ylim) - ylims!(ax, ylim) - end - if !isnothing(axisaspect) - ax.aspect = AxisAspect(axisaspect) - end - Colorbar(fig[1, 2], limits=colorrange) - fig, ax, ob - end - - function plot(::Type{ObExpr{:Form2}}, s, sd, vals, colorrange, xlim, ylim, axisaspect; kw...) - # Split mesh into component triangles - m = GeometryBasics.Mesh(s) - x = faces(m) - m_points = m.position[vcat([[t[1],t[2],t[3]] for t in x]...)] - m_faces = TriangleFace{Int}[[((t-1) * 3) .+ (1,2,3) for t in 1:length(x)]...] - new_m = GeometryBasics.Mesh(Point{3, Float64}[m_points...], m_faces) - fig = Figure() - ax, ob = mesh(fig[1,1], new_m, color=vcat([[v,v,v] for v in vals]...)) - if !isnothing(colorrange) - ob.colorrange = colorrange - end - if !isnothing(xlim) - xlims!(ax, xlim) - end - if !isnothing(ylim) - ylims!(ax, ylim) - end - if !isnothing(axisaspect) - ax.aspect = AxisAspect(axisaspect) - end - Colorbar(fig[1, 2], ob) - fig, ax, ob - end -end diff --git a/src/Decapodes.jl b/src/Decapodes.jl index 18cd770..684e562 100644 --- a/src/Decapodes.jl +++ b/src/Decapodes.jl @@ -1,16 +1,41 @@ module Decapodes -using Requires +using Catlab +using Catlab.Theories +import Catlab.Theories: otimes, oplus, compose, ⊗, ⊕, ⋅, associate, associate_unit, Ob, Hom, dom, codom +using Catlab.Present +using Catlab.Programs +using Catlab.CategoricalAlgebra +using Catlab.WiringDiagrams +using Catlab.WiringDiagrams.DirectedWiringDiagrams +using LinearAlgebra +using MLStyle +using Base.Iterators -include("Decapodes2/decapodes.jl") -include("Diagrams.jl") -include("OpenDiagrams.jl") -include("Schedules.jl") -include("Simulations.jl") -include("Examples.jl") +import Unicode + +export normalize_unicode, DerivOp, append_dot, + SchDecapode, SchNamedDecapode, AbstractDecapode, AbstractNamedDecapode, Decapode, NamedDecapode, SummationDecapode, fill_names!, expand_operators, add_constant!, add_parameter, infer_types!, resolve_overloads!, op1_inf_rules_1D, op2_inf_rules_1D, op1_inf_rules_2D, op2_inf_rules_2D, op1_res_rules_1D, op2_res_rules_1D, op1_res_rules_2D, op2_res_rules_2D, + Term, Var, Judge, Eq, AppCirc1, AppCirc2, App1, App2, Plus, Tan, term, parse_decapode, + VectorForm, PhysicsState, findname, findnode, + compile, compile_env, gensim, closest_point, flat_op, + AbstractMeshKey, loadmesh, Icosphere, Rectangle_30x10, Torus_30x10, Point_Map, + Open, OpenSummationDecapodeOb, OpenSummationDecapode, unique_by, unique_by!, oapply, + CartesianPoint, SpherePoint, r, theta, phi, TangentBasis, θhat, ϕhat, + average_rewrite + +normalize_unicode(s::String) = Unicode.normalize(s, compose=true, stable=true, chartransform=Unicode.julia_chartransform) +normalize_unicode(s::Symbol) = Symbol(normalize_unicode(String(s))) +DerivOp = Symbol("∂ₜ") +append_dot(s::Symbol) = Symbol(string(s)*'\U0307') + +include("decapodeacset.jl") +include("language.jl") +include("composition.jl") +include("coordinates.jl") +include("visualization.jl") +include("simulation.jl") +include("meshes.jl") +include("rewrite.jl") -function __init__() - @require CairoMakie="13f3f980-e62b-5c42-98c6-ff1f3baf88f0" include("Debug.jl") - @require AlgebraicPetri="4f99eebe-17bf-4e98-b6a1-2c4f205a959b" include("PetriNets.jl") -end end diff --git a/src/Decapodes2/decapodes.jl b/src/Decapodes2/decapodes.jl deleted file mode 100644 index 43e15ca..0000000 --- a/src/Decapodes2/decapodes.jl +++ /dev/null @@ -1,38 +0,0 @@ -using Catlab -using Catlab.Theories -import Catlab.Theories: otimes, oplus, compose, ⊗, ⊕, ⋅, associate, associate_unit, Ob, Hom, dom, codom -using Catlab.Present -using Catlab.Programs -using Catlab.CategoricalAlgebra -using Catlab.WiringDiagrams -using Catlab.WiringDiagrams.DirectedWiringDiagrams -using LinearAlgebra -using MLStyle -using Base.Iterators - -import Unicode - -export normalize_unicode, DerivOp, append_dot, - SchDecapode, SchNamedDecapode, AbstractDecapode, AbstractNamedDecapode, Decapode, NamedDecapode, SummationDecapode, fill_names!, expand_operators, add_constant!, add_parameter, infer_types!, resolve_overloads!, op1_inf_rules_1D, op2_inf_rules_1D, op1_inf_rules_2D, op2_inf_rules_2D, op1_res_rules_1D, op2_res_rules_1D, op1_res_rules_2D, op2_res_rules_2D, - Term, Var, Judge, Eq, AppCirc1, AppCirc2, App1, App2, Plus, Tan, term, parse_decapode, - VectorForm, PhysicsState, findname, findnode, - compile, compile_env, gensim, closest_point, flat_op, - AbstractMeshKey, loadmesh, Icosphere, Rectangle_30x10, Torus_30x10, Point_Map, - Open, OpenSummationDecapodeOb, OpenSummationDecapode, unique_by, unique_by!, oapply, - CartesianPoint, SpherePoint, r, theta, phi, TangentBasis, θhat, ϕhat, - average_rewrite - -normalize_unicode(s::String) = Unicode.normalize(s, compose=true, stable=true, chartransform=Unicode.julia_chartransform) -normalize_unicode(s::Symbol) = Symbol(normalize_unicode(String(s))) -DerivOp = Symbol("∂ₜ") -append_dot(s::Symbol) = Symbol(string(s)*'\U0307') - -include("decapodeacset.jl") -include("language.jl") -include("composition.jl") -include("coordinates.jl") -include("visualization.jl") -include("simulation.jl") -include("meshes.jl") -include("rewrite.jl") - diff --git a/src/Diagrams.jl b/src/Diagrams.jl deleted file mode 100644 index 98124d1..0000000 --- a/src/Diagrams.jl +++ /dev/null @@ -1,329 +0,0 @@ -""" Provides the tooling for defining Decapodes conveniently with symbolic -expressions -""" -module Diagrams - -using CombinatorialSpaces.ExteriorCalculus -using Catlab -using Catlab.CategoricalAlgebra -using Catlab.Theories -using Catlab.Programs.DiagrammaticPrograms.AST -using Catlab.Programs.DiagrammaticPrograms: parse_diagram, parse_gat_expr, parse_diagram_ast -using Unicode - -export @decapode, Decapodes2D - -""" decapode(cat, body) - -This macro is called as: -``` -@decapode presentation begin - expressions -end -``` -where `presentation` is a Presentation containing the syntax for the Decapode -(operators and form types), and `expressions` is symbolic expressions of the -relations within the Decapode. - -""" -macro decapode(cat, body) - :(parse_exp_diagram($(esc(cat)), $(Meta.quot(body)))) -end - -function parse_exp_diagram(cat, body; kw...) - objs = filter(exp -> !(exp isa LineNumberNode) && (exp.head == :(::)), body.args) - obj_map = Dict{Symbol, Expr}() - for obj in objs - names = obj.args[1] - type = obj.args[2] - # Expand if object is a tuple or single object - if names isa Symbol - obj_map[names] = type - else - for name in names.args - obj_map[name] = type - end - end - end - new_body = quote end - for exp in body.args - if !(exp isa LineNumberNode) && exp.head == :call && exp.args[1] == :(==) - exprs = [] - if exp.args[2] isa Symbol - res, _ = expand_expr!(exprs, exp.args[3], cat, obj_map; make_var = false) - push!(exprs, :($(exp.args[2]) == $(res))) - elseif exp.args[3] isa Symbol - res, _ = expand_expr!(exprs, exp.args[2], cat, obj_map; make_var = false) - push!(exprs, :($(exp.args[3]) == $(res))) - else - res1, _ = expand_expr!(exprs, exp.args[2], cat, obj_map; make_var = true) - res2, _ = expand_expr!(exprs, exp.args[3], cat, obj_map; make_var = false) - push!(exprs, :($(res1) == $(res2))) - end - append!(new_body.args, exprs) - else - push!(new_body.args, exp) - end - end - free_new_body = AST.Diagram(parse_diagram_ast(new_body, free=true)) - parse_diagram(cat, free_new_body) -end - -i2sub = Dict( - '0'=>'₀', '1'=>'₁', '2'=>'₂', '3'=>'₃', '4'=>'₄', '5'=>'₅', - '6'=>'₆', '7'=>'₇', '8'=>'₈', '9'=>'₉', '-'=>'₋' -) -i2sup = Dict( - '0'=>'⁰', '1'=>'¹', '2'=>'²', '3'=>'³', '4'=>'⁴', '5'=>'⁵', - '6'=>'⁶', '7'=>'⁷', '8'=>'⁸', '9'=>'⁹', '-'=>'⁻' -) - -form_sup(type::Expr) = form_sup(type.args[1]) -form_sub(type::Expr) = form_sub(type.args[1]) - -form_sup(type::Symbol) = startswith("$type", "Dual") ? - i2sup[last("$type")] * '̃' : i2sup[last("$type")] -form_sub(type::Symbol) = startswith("$type", "Dual") ? i2sub[last("$type")] * '̃' : i2sub[last("$type")] - - -# Currently every projection needs its own index -function make_proj(ind, types) - op = "proj" * i2sub["$ind"[1]] * "_" * join(form_sup.(types)) * form_sub(types[ind]) - Symbol(Unicode.normalize(op)) -end - -function parse_type(type) - if head(type) == :otimes - Meta.parse("otimes{$(join(["$(head(t)){$(t.args[1])}" for t in type.args], ","))}") - else - Meta.parse("$(head(type)){$(type.args[1])}") - end -end - - -get_name(expr) = expr isa Symbol ? expr : get_name(expr.args[1]) - -function expand_expr!(expr_arr, expr, cat, obj_map; make_var = true) - if expr isa Symbol - expr, obj_map[expr] - elseif length(expr.args) > 2 - if expr.args[1] == :+ - # Only process 2 of these at a time (sum first, evaluate second as sum) - res1 = expand_expr!(expr_arr, expr.args[2], cat, obj_map; make_var = true) - - res2 = length(expr.args) > 3 ? - expand_expr!(expr_arr, :(+($(expr.args[3:end]...))), cat, obj_map; make_var = true) : - expand_expr!(expr_arr, expr.args[3], cat, obj_map; make_var = true) - type = res1[2] - - sum_name = Symbol(Unicode.normalize("sum" * form_sub(type))) - new_var = gensym() - nv_type = :(otimes{$(res1[2]), $(res2[2])}) - push!(expr_arr, :($(new_var)::$(nv_type))) - push!(expr_arr, :($(res1[1]) == $(make_proj(1, [res1[2], res2[2]]))($(new_var)))) - push!(expr_arr, :($(res2[1]) == $(make_proj(2, [res1[2], res2[2]]))($(new_var)))) - if make_var - name = gensym() - push!(expr_arr, :($(name)::$(type))) - push!(expr_arr, :($(name) == $(sum_name)($(new_var)))) - name, type - else - :($(sum_name)($(new_var))), type - end - else - res = map(expr.args[2:end]) do ex - expand_expr!(expr_arr, ex, cat, obj_map; make_var = true) - end - new_var = gensym() - nv_type = Meta.parse("otimes{$(join(last.(res), ","))}") - push!(expr_arr, :($(new_var)::$(nv_type))) - append!(expr_arr, map(enumerate(res)) do (i, obj) - :($(obj[1]) == $(make_proj(i, last.(res)))($(new_var))) - end) - type = parse_type(codom(parse_gat_expr(FinCat(cat), expr.args[1]))) - if make_var - name = gensym() - push!(expr_arr, :($(name)::$(type))) - push!(expr_arr, :($(name) == $(expr.args[1])($(new_var)))) - name, type - else - :($(expr.args[1])($(new_var))), type - end - end - else - type = parse_type(codom(parse_gat_expr(FinCat(cat), expr.args[1]))) - arg_var, arg_type = expand_expr!(expr_arr, expr.args[2], cat, obj_map; make_var = false) - expr.args[2] = arg_var - cur_var = expr - if make_var - name = gensym() - push!(expr_arr, :($(name)::$(type))) - push!(expr_arr, :($(name) == $(expr))) - cur_var = name - end - cur_var, type - end -end - -""" Decapodes2D -A schema which includes any homomorphisms that may be added by the @decapode -macro. - -TODO: This should be chipped away at as more of this tooling takes advantage -of the Catlab GAT system -""" -@present Decapodes2D(FreeExtCalc2D) begin - X::Space - proj₁_⁰⁰₀::Hom(Form0(X)⊗Form0(X),Form0(X)) - proj₂_⁰⁰₀::Hom(Form0(X)⊗Form0(X),Form0(X)) - proj₁_⁰⁰₀⁺::Hom(Form0(X)⊕Form0(X),Form0(X)) - proj₂_⁰⁰₀⁺::Hom(Form0(X)⊕Form0(X),Form0(X)) - proj₁_⁰¹₀::Hom(Form0(X)⊗Form1(X),Form0(X)) - proj₂_⁰¹₁::Hom(Form0(X)⊗Form1(X),Form1(X)) - proj₁_⁰¹₀⁺::Hom(Form0(X)⊕Form1(X),Form0(X)) - proj₂_⁰¹₁⁺::Hom(Form0(X)⊕Form1(X),Form1(X)) - proj₁_⁰²₀::Hom(Form0(X)⊗Form2(X),Form0(X)) - proj₂_⁰²₂::Hom(Form0(X)⊗Form2(X),Form2(X)) - proj₁_⁰²₀⁺::Hom(Form0(X)⊕Form2(X),Form0(X)) - proj₂_⁰²₂⁺::Hom(Form0(X)⊕Form2(X),Form2(X)) - proj₁_⁰⁰̃₀::Hom(Form0(X)⊗DualForm0(X),Form0(X)) - proj₂_⁰⁰̃₀̃::Hom(Form0(X)⊗DualForm0(X),DualForm0(X)) - proj₁_⁰⁰̃₀⁺::Hom(Form0(X)⊕DualForm0(X),Form0(X)) - proj₂_⁰⁰̃₀̃⁺::Hom(Form0(X)⊕DualForm0(X),DualForm0(X)) - proj₁_⁰¹̃₀::Hom(Form0(X)⊗DualForm1(X),Form0(X)) - proj₂_⁰¹̃₁̃::Hom(Form0(X)⊗DualForm1(X),DualForm1(X)) - proj₁_⁰¹̃₀⁺::Hom(Form0(X)⊕DualForm1(X),Form0(X)) - proj₂_⁰¹̃₁̃⁺::Hom(Form0(X)⊕DualForm1(X),DualForm1(X)) - proj₁_⁰²̃₀::Hom(Form0(X)⊗DualForm2(X),Form0(X)) - proj₂_⁰²̃₂̃::Hom(Form0(X)⊗DualForm2(X),DualForm2(X)) - proj₁_⁰²̃₀⁺::Hom(Form0(X)⊕DualForm2(X),Form0(X)) - proj₂_⁰²̃₂̃⁺::Hom(Form0(X)⊕DualForm2(X),DualForm2(X)) - proj₁_¹⁰₁::Hom(Form1(X)⊗Form0(X),Form1(X)) - proj₂_¹⁰₀::Hom(Form1(X)⊗Form0(X),Form0(X)) - proj₁_¹⁰₁⁺::Hom(Form1(X)⊕Form0(X),Form1(X)) - proj₂_¹⁰₀⁺::Hom(Form1(X)⊕Form0(X),Form0(X)) - proj₁_¹¹₁::Hom(Form1(X)⊗Form1(X),Form1(X)) - proj₂_¹¹₁::Hom(Form1(X)⊗Form1(X),Form1(X)) - proj₁_¹¹₁⁺::Hom(Form1(X)⊕Form1(X),Form1(X)) - proj₂_¹¹₁⁺::Hom(Form1(X)⊕Form1(X),Form1(X)) - proj₁_¹²₁::Hom(Form1(X)⊗Form2(X),Form1(X)) - proj₂_¹²₂::Hom(Form1(X)⊗Form2(X),Form2(X)) - proj₁_¹²₁⁺::Hom(Form1(X)⊕Form2(X),Form1(X)) - proj₂_¹²₂⁺::Hom(Form1(X)⊕Form2(X),Form2(X)) - proj₁_¹⁰̃₁::Hom(Form1(X)⊗DualForm0(X),Form1(X)) - proj₂_¹⁰̃₀̃::Hom(Form1(X)⊗DualForm0(X),DualForm0(X)) - proj₁_¹⁰̃₁⁺::Hom(Form1(X)⊕DualForm0(X),Form1(X)) - proj₂_¹⁰̃₀̃⁺::Hom(Form1(X)⊕DualForm0(X),DualForm0(X)) - proj₁_¹¹̃₁::Hom(Form1(X)⊗DualForm1(X),Form1(X)) - proj₂_¹¹̃₁̃::Hom(Form1(X)⊗DualForm1(X),DualForm1(X)) - proj₁_¹¹̃₁⁺::Hom(Form1(X)⊕DualForm1(X),Form1(X)) - proj₂_¹¹̃₁̃⁺::Hom(Form1(X)⊕DualForm1(X),DualForm1(X)) - proj₁_¹²̃₁::Hom(Form1(X)⊗DualForm2(X),Form1(X)) - proj₂_¹²̃₂̃::Hom(Form1(X)⊗DualForm2(X),DualForm2(X)) - proj₁_¹²̃₁⁺::Hom(Form1(X)⊕DualForm2(X),Form1(X)) - proj₂_¹²̃₂̃⁺::Hom(Form1(X)⊕DualForm2(X),DualForm2(X)) - proj₁_²⁰₂::Hom(Form2(X)⊗Form0(X),Form2(X)) - proj₂_²⁰₀::Hom(Form2(X)⊗Form0(X),Form0(X)) - proj₁_²⁰₂⁺::Hom(Form2(X)⊕Form0(X),Form2(X)) - proj₂_²⁰₀⁺::Hom(Form2(X)⊕Form0(X),Form0(X)) - proj₁_²¹₂::Hom(Form2(X)⊗Form1(X),Form2(X)) - proj₂_²¹₁::Hom(Form2(X)⊗Form1(X),Form1(X)) - proj₁_²¹₂⁺::Hom(Form2(X)⊕Form1(X),Form2(X)) - proj₂_²¹₁⁺::Hom(Form2(X)⊕Form1(X),Form1(X)) - proj₁_²²₂::Hom(Form2(X)⊗Form2(X),Form2(X)) - proj₂_²²₂::Hom(Form2(X)⊗Form2(X),Form2(X)) - proj₁_²²₂⁺::Hom(Form2(X)⊕Form2(X),Form2(X)) - proj₂_²²₂⁺::Hom(Form2(X)⊕Form2(X),Form2(X)) - proj₁_²⁰̃₂::Hom(Form2(X)⊗DualForm0(X),Form2(X)) - proj₂_²⁰̃₀̃::Hom(Form2(X)⊗DualForm0(X),DualForm0(X)) - proj₁_²⁰̃₂⁺::Hom(Form2(X)⊕DualForm0(X),Form2(X)) - proj₂_²⁰̃₀̃⁺::Hom(Form2(X)⊕DualForm0(X),DualForm0(X)) - proj₁_²¹̃₂::Hom(Form2(X)⊗DualForm1(X),Form2(X)) - proj₂_²¹̃₁̃::Hom(Form2(X)⊗DualForm1(X),DualForm1(X)) - proj₁_²¹̃₂⁺::Hom(Form2(X)⊕DualForm1(X),Form2(X)) - proj₂_²¹̃₁̃⁺::Hom(Form2(X)⊕DualForm1(X),DualForm1(X)) - proj₁_²²̃₂::Hom(Form2(X)⊗DualForm2(X),Form2(X)) - proj₂_²²̃₂̃::Hom(Form2(X)⊗DualForm2(X),DualForm2(X)) - proj₁_²²̃₂⁺::Hom(Form2(X)⊕DualForm2(X),Form2(X)) - proj₂_²²̃₂̃⁺::Hom(Form2(X)⊕DualForm2(X),DualForm2(X)) - proj₁_⁰̃⁰₀̃::Hom(DualForm0(X)⊗Form0(X),DualForm0(X)) - proj₂_⁰̃⁰₀::Hom(DualForm0(X)⊗Form0(X),Form0(X)) - proj₁_⁰̃⁰₀̃⁺::Hom(DualForm0(X)⊕Form0(X),DualForm0(X)) - proj₂_⁰̃⁰₀⁺::Hom(DualForm0(X)⊕Form0(X),Form0(X)) - proj₁_⁰̃¹₀̃::Hom(DualForm0(X)⊗Form1(X),DualForm0(X)) - proj₂_⁰̃¹₁::Hom(DualForm0(X)⊗Form1(X),Form1(X)) - proj₁_⁰̃¹₀̃⁺::Hom(DualForm0(X)⊕Form1(X),DualForm0(X)) - proj₂_⁰̃¹₁⁺::Hom(DualForm0(X)⊕Form1(X),Form1(X)) - proj₁_⁰̃²₀̃::Hom(DualForm0(X)⊗Form2(X),DualForm0(X)) - proj₂_⁰̃²₂::Hom(DualForm0(X)⊗Form2(X),Form2(X)) - proj₁_⁰̃²₀̃⁺::Hom(DualForm0(X)⊕Form2(X),DualForm0(X)) - proj₂_⁰̃²₂⁺::Hom(DualForm0(X)⊕Form2(X),Form2(X)) - proj₁_⁰̃⁰̃₀̃::Hom(DualForm0(X)⊗DualForm0(X),DualForm0(X)) - proj₂_⁰̃⁰̃₀̃::Hom(DualForm0(X)⊗DualForm0(X),DualForm0(X)) - proj₁_⁰̃⁰̃₀̃⁺::Hom(DualForm0(X)⊕DualForm0(X),DualForm0(X)) - proj₂_⁰̃⁰̃₀̃⁺::Hom(DualForm0(X)⊕DualForm0(X),DualForm0(X)) - proj₁_⁰̃¹̃₀̃::Hom(DualForm0(X)⊗DualForm1(X),DualForm0(X)) - proj₂_⁰̃¹̃₁̃::Hom(DualForm0(X)⊗DualForm1(X),DualForm1(X)) - proj₁_⁰̃¹̃₀̃⁺::Hom(DualForm0(X)⊕DualForm1(X),DualForm0(X)) - proj₂_⁰̃¹̃₁̃⁺::Hom(DualForm0(X)⊕DualForm1(X),DualForm1(X)) - proj₁_⁰̃²̃₀̃::Hom(DualForm0(X)⊗DualForm2(X),DualForm0(X)) - proj₂_⁰̃²̃₂̃::Hom(DualForm0(X)⊗DualForm2(X),DualForm2(X)) - proj₁_⁰̃²̃₀̃⁺::Hom(DualForm0(X)⊕DualForm2(X),DualForm0(X)) - proj₂_⁰̃²̃₂̃⁺::Hom(DualForm0(X)⊕DualForm2(X),DualForm2(X)) - proj₁_¹̃⁰₁̃::Hom(DualForm1(X)⊗Form0(X),DualForm1(X)) - proj₂_¹̃⁰₀::Hom(DualForm1(X)⊗Form0(X),Form0(X)) - proj₁_¹̃⁰₁̃⁺::Hom(DualForm1(X)⊕Form0(X),DualForm1(X)) - proj₂_¹̃⁰₀⁺::Hom(DualForm1(X)⊕Form0(X),Form0(X)) - proj₁_¹̃¹₁̃::Hom(DualForm1(X)⊗Form1(X),DualForm1(X)) - proj₂_¹̃¹₁::Hom(DualForm1(X)⊗Form1(X),Form1(X)) - proj₁_¹̃¹₁̃⁺::Hom(DualForm1(X)⊕Form1(X),DualForm1(X)) - proj₂_¹̃¹₁⁺::Hom(DualForm1(X)⊕Form1(X),Form1(X)) - proj₁_¹̃²₁̃::Hom(DualForm1(X)⊗Form2(X),DualForm1(X)) - proj₂_¹̃²₂::Hom(DualForm1(X)⊗Form2(X),Form2(X)) - proj₁_¹̃²₁̃⁺::Hom(DualForm1(X)⊕Form2(X),DualForm1(X)) - proj₂_¹̃²₂⁺::Hom(DualForm1(X)⊕Form2(X),Form2(X)) - proj₁_¹̃⁰̃₁̃::Hom(DualForm1(X)⊗DualForm0(X),DualForm1(X)) - proj₂_¹̃⁰̃₀̃::Hom(DualForm1(X)⊗DualForm0(X),DualForm0(X)) - proj₁_¹̃⁰̃₁̃⁺::Hom(DualForm1(X)⊕DualForm0(X),DualForm1(X)) - proj₂_¹̃⁰̃₀̃⁺::Hom(DualForm1(X)⊕DualForm0(X),DualForm0(X)) - proj₁_¹̃¹̃₁̃::Hom(DualForm1(X)⊗DualForm1(X),DualForm1(X)) - proj₂_¹̃¹̃₁̃::Hom(DualForm1(X)⊗DualForm1(X),DualForm1(X)) - proj₁_¹̃¹̃₁̃⁺::Hom(DualForm1(X)⊕DualForm1(X),DualForm1(X)) - proj₂_¹̃¹̃₁̃⁺::Hom(DualForm1(X)⊕DualForm1(X),DualForm1(X)) - proj₁_¹̃²̃₁̃::Hom(DualForm1(X)⊗DualForm2(X),DualForm1(X)) - proj₂_¹̃²̃₂̃::Hom(DualForm1(X)⊗DualForm2(X),DualForm2(X)) - proj₁_¹̃²̃₁̃⁺::Hom(DualForm1(X)⊕DualForm2(X),DualForm1(X)) - proj₂_¹̃²̃₂̃⁺::Hom(DualForm1(X)⊕DualForm2(X),DualForm2(X)) - proj₁_²̃⁰₂̃::Hom(DualForm2(X)⊗Form0(X),DualForm2(X)) - proj₂_²̃⁰₀::Hom(DualForm2(X)⊗Form0(X),Form0(X)) - proj₁_²̃⁰₂̃⁺::Hom(DualForm2(X)⊕Form0(X),DualForm2(X)) - proj₂_²̃⁰₀⁺::Hom(DualForm2(X)⊕Form0(X),Form0(X)) - proj₁_²̃¹₂̃::Hom(DualForm2(X)⊗Form1(X),DualForm2(X)) - proj₂_²̃¹₁::Hom(DualForm2(X)⊗Form1(X),Form1(X)) - proj₁_²̃¹₂̃⁺::Hom(DualForm2(X)⊕Form1(X),DualForm2(X)) - proj₂_²̃¹₁⁺::Hom(DualForm2(X)⊕Form1(X),Form1(X)) - proj₁_²̃²₂̃::Hom(DualForm2(X)⊗Form2(X),DualForm2(X)) - proj₂_²̃²₂::Hom(DualForm2(X)⊗Form2(X),Form2(X)) - proj₁_²̃²₂̃⁺::Hom(DualForm2(X)⊕Form2(X),DualForm2(X)) - proj₂_²̃²₂⁺::Hom(DualForm2(X)⊕Form2(X),Form2(X)) - proj₁_²̃⁰̃₂̃::Hom(DualForm2(X)⊗DualForm0(X),DualForm2(X)) - proj₂_²̃⁰̃₀̃::Hom(DualForm2(X)⊗DualForm0(X),DualForm0(X)) - proj₁_²̃⁰̃₂̃⁺::Hom(DualForm2(X)⊕DualForm0(X),DualForm2(X)) - proj₂_²̃⁰̃₀̃⁺::Hom(DualForm2(X)⊕DualForm0(X),DualForm0(X)) - proj₁_²̃¹̃₂̃::Hom(DualForm2(X)⊗DualForm1(X),DualForm2(X)) - proj₂_²̃¹̃₁̃::Hom(DualForm2(X)⊗DualForm1(X),DualForm1(X)) - proj₁_²̃¹̃₂̃⁺::Hom(DualForm2(X)⊕DualForm1(X),DualForm2(X)) - proj₂_²̃¹̃₁̃⁺::Hom(DualForm2(X)⊕DualForm1(X),DualForm1(X)) - proj₁_²̃²̃₂̃::Hom(DualForm2(X)⊗DualForm2(X),DualForm2(X)) - proj₂_²̃²̃₂̃::Hom(DualForm2(X)⊗DualForm2(X),DualForm2(X)) - proj₁_²̃²̃₂̃⁺::Hom(DualForm2(X)⊕DualForm2(X),DualForm2(X)) - proj₂_²̃²̃₂̃⁺::Hom(DualForm2(X)⊕DualForm2(X),DualForm2(X)) - sum₀::Hom(Form0(X)⊗Form0(X),Form0(X)) - sum₁::Hom(Form1(X)⊗Form1(X),Form1(X)) - sum₂::Hom(Form2(X)⊗Form2(X),Form2(X)) - sum₀̃::Hom(DualForm0(X)⊗DualForm0(X),DualForm0(X)) - sum₁̃::Hom(DualForm1(X)⊗DualForm1(X),DualForm1(X)) - sum₂̃::Hom(DualForm2(X)⊗DualForm2(X),DualForm2(X)) -end - -end diff --git a/src/Examples.jl b/src/Examples.jl deleted file mode 100644 index 3f18980..0000000 --- a/src/Examples.jl +++ /dev/null @@ -1,330 +0,0 @@ -""" Provides helper functions for generating and optimizing Decapode simulations - - -*TODO*: This module was originally intended to be "an example of open diagrams -applied to DEC", but this should be given a new name as its methods have become -more critical -""" - -module Examples - -using CombinatorialSpaces -using CombinatorialSpaces.ExteriorCalculus -using Catlab -using Catlab.Present -using Catlab.Programs -using Catlab.WiringDiagrams -using Catlab.Graphics -using Catlab.CategoricalAlgebra -using Decapodes.Simulations -using Decapodes.Diagrams -using Decapodes.Schedules -using CombinatorialSpaces: ∧ - -using LinearAlgebra - -import Catlab.Graphics: to_graphviz - -export dual, sym2func, to_graphviz, gen_dec_rules, contract_matrices, expand_dwd, zip_dwd - -"""dual(s::EmbeddedDeltaSet2D{O, P}) - -Generates a dual mesh for the provided delta set `s`. - -*TODO*: This tooling should be upstreamed to CombinatorialSpaces -""" -function dual(s::EmbeddedDeltaSet2D{O, P}) where {O, P} - sd = EmbeddedDeltaDualComplex2D{O, eltype(P), P}(s) - subdivide_duals!(sd, Barycenter()) - sd -end - -sym2func(sd) = begin - s2f = Dict(:d₀=>Dict(:operator => d(Val{0}, sd), :type => MatrixFunc()), - :d₁=>Dict(:operator => d(Val{1}, sd), :type => MatrixFunc()), - :dual_d₀=>Dict(:operator => dual_derivative(Val{0}, sd), - :type => MatrixFunc()), - :dual_d₁=>Dict(:operator => dual_derivative(Val{1}, sd), - :type => MatrixFunc()), - :⋆₀=>Dict(:operator => hodge_star(Val{0}, sd), :type => MatrixFunc()), - :⋆₁=>Dict(:operator => hodge_star(Val{1}, sd), :type => MatrixFunc()), - :⋆₂=>Dict(:operator => hodge_star(Val{2}, sd), :type => MatrixFunc()), - :⋆₀⁻¹=>Dict(:operator => inv_hodge_star(Val{0}, sd), :type => MatrixFunc()), - :⋆₁⁻¹=>Dict(:operator => inv_hodge_star(Val{1}, sd; hodge=DiagonalHodge()), - :type => MatrixFunc()), - :⋆₂⁻¹=>Dict(:operator => inv_hodge_star(Val{2}, sd), :type => MatrixFunc()), - :∧₁₀=>Dict(:operator => (γ, α,β)->(γ .= ∧(Tuple{1,0}, sd, α, β)), - :type=>InPlaceFunc()), - :∧₁₁=>Dict(:operator => (γ, α,β)->(γ .= ∧(Tuple{1,1}, sd, α, β)), - :type=>InPlaceFunc()), - :plus => Dict(:operator => (x,y)->(x′ .= x .+ y), :type => InPlaceFunc()), - :plus => Dict(:operator => (x,y)->(x+y), :type => ElementwiseFunc()), - :L₀ => Dict(:operator => (v,u)->(lie_derivative_flat(Val{2}, sd, v, u)), - :type => ArbitraryFunc())) - for s in [:sum₀, :sum₁, :sum₂, :sum₀̃, :sum₁̃, :sum₂̃] - s2f[s] = Dict(:operator => (x′,x,y)->(x′ .= x .+ y), :type => InPlaceFunc()) - end - s2f -end - - -""" expand_dwd(dwd_orig::WiringDiagram, patterns::Dict{Symbol, WiringDiagram}) - -Replaces any boxes in `dwd_orig` which map to a key in `patterns` with their -corresponding diagram in `patterns`. This allows for replacement of complex -operations with their definition. -""" -function expand_dwd(dwd_orig, patterns) - dwd = deepcopy(dwd_orig) - has_updated = true - while(has_updated) - updates = map(1:nboxes(dwd)) do b - v = box(dwd, b).value - if v ∈ keys(patterns) - dwd.diagram[b, :value] = patterns[v] - true - else - false - end - end - dwd = substitute(dwd) - has_updated = any(updates) - end - dwd -end - -""" contract_matrices(dwd::WiringDiagram, s2f::Dict) - -Pre-computes matrix products when two matrix multiplication operations -immediately follow each other. -""" - -function contract_matrices!(dwd, s2f) - has_updated = true - is_matrix(b) = begin - b >= 1 || return false - v = Symbol(first(split("$(box(dwd, b).value)", "⋅"))) - v != :remove && (s2f[v][:type] isa MatrixFunc) && - length(input_ports(dwd, b)) == 1 && - length(output_ports(dwd, b)) == 1 - end - while(has_updated) - has_updated = false - # Join parallel matrix boxes and mark for removal - for b in 1:nboxes(dwd) - if is_matrix(b) - ows = out_wires(dwd, b) - if length(ows) == 1 - ow = ows[1] - if is_matrix(ow.target.box) - has_updated = true - tb = ow.target.box - dwd.diagram[b, :value] = Symbol(dwd.diagram[tb, :value], :⋅, dwd.diagram[b, :value]) - dwd.diagram[tb, :value] = :remove - add_wires!(dwd, map(w -> Wire(w.value, Port(b, OutputPort, 1), w.target), - out_wires(dwd, tb))) - dwd.diagram[incident(dwd.diagram, b, :out_port_box)[1], :out_port_type] = output_ports(dwd, tb)[1] - end - end - end - end - rem_boxes_boot!(dwd, findall(b -> dwd.diagram[b, :value] == :remove, 1:nboxes(dwd))) - end -end - -""" rem_boxes_boot! - -This function removes boxes from a DWD in a way that preserves port ordering. -This is only a temporary fix, and is dependent on [this -issue](https://github.com/AlgebraicJulia/Catlab.jl/issues/530) in the upstream -dependency Catlab -""" -function rem_boxes_boot!(dwd, box_to_rem) - if isempty(box_to_rem) - return - end - boxes = collect(parts(dwd.diagram, :Box)) - filter!(b -> !(b ∈ box_to_rem), boxes) - append!(boxes, box_to_rem) - in_ports = incident(dwd.diagram, boxes, :in_port_box) - out_ports = incident(dwd.diagram, boxes, :out_port_box) - - in_port_trans = fill(1, nparts(dwd.diagram, :InPort)) - out_port_trans = fill(1, nparts(dwd.diagram, :OutPort)) - in_port_trans[vcat(in_ports...)] .= parts(dwd.diagram, :InPort) - out_port_trans[vcat(out_ports...)] .= parts(dwd.diagram, :OutPort) - - set_subparts!(dwd.diagram, parts(dwd.diagram, :InPort), - in_port_box = copy(dwd.diagram[vcat(in_ports...), :in_port_box]), - in_port_type= copy(dwd.diagram[vcat(in_ports...), :in_port_type])) - set_subparts!(dwd.diagram, parts(dwd.diagram, :OutPort), - out_port_box = copy(dwd.diagram[vcat(out_ports...), :out_port_box]), - out_port_type= copy(dwd.diagram[vcat(out_ports...), :out_port_type])) - - set_subparts!(dwd.diagram, parts(dwd.diagram, :Wire); - src = out_port_trans[dwd.diagram[:src]], - tgt = in_port_trans[dwd.diagram[:tgt]]) - set_subparts!(dwd.diagram, parts(dwd.diagram, :InWire), - in_tgt = in_port_trans[dwd.diagram[:in_tgt]]) - set_subparts!(dwd.diagram, parts(dwd.diagram, :OutWire), - out_src = out_port_trans[dwd.diagram[:out_src]]) - - rem_boxes!(dwd, box_to_rem) -end - -""" zip_dwd!(dwd::WiringDiagram) - -This function optimizes the wiring diagram by removing redundant operations. -Takes advantage of the fact that in a deterministic wiring diagram, the -identical operation being performed on the same value will result in the same -value. This means that parallel computations of the same operation can be -"zipped" into a single execution of the operation. -""" - -function zip_dwd!(dwd) - updated = true - while updated - updated = false - box_inputs = map(parts(dwd.diagram, :Box)) do b - [dwd.diagram[b, :value], [(w.source.box, w.source.port) for w in in_wires(dwd, b)]] - end - for b1 in 1:length(box_inputs) - for b2 in (b1+1):length(box_inputs) - if box_inputs[b1] == box_inputs[b2] - b1_oports = incident(dwd.diagram, b1, :out_port_box) - b2_oports = incident(dwd.diagram, b2, :out_port_box) - for p in 1:length(b1_oports) - b2wires = incident(dwd.diagram, b2_oports[p], :src) - set_subparts!(dwd.diagram, b2wires, src=b1_oports[p]) - b2outwires = incident(dwd.diagram, b2_oports[p], :out_src) - set_subparts!(dwd.diagram, b2outwires, out_src=b1_oports[p]) - end - rem_boxes_boot!(dwd, [b2]) - updated = true - break - end - end - if updated - break - end - end - end -end - - -function boundary_inds(::Type{Val{1}}, s) - collect(findall(x -> x != 0, boundary(Val{2},s) * fill(1,ntriangles(s)))) -end - -function boundary_inds(::Type{Val{0}}, s) - ∂1_inds = boundary_inds(Val{1}, s) - unique(vcat(s[∂1_inds,:∂v0],s[∂1_inds,:∂v1])) -end - -function boundary_inds(::Type{Val{2}}, s) - ∂1_inds = boundary_inds(Val{1}, s) - inds = map([:∂e0, :∂e1, :∂e2]) do esym - vcat(incident(s, ∂1_inds, esym)...) - end - unique(vcat(inds...)) -end - - -function bound_edges(s, ∂₀) - te = vcat(incident(s, ∂₀, :∂v1)...) - se = vcat(incident(s, ∂₀, :∂v0)...) - intersect(te, se) -end - -function adj_edges(s, ∂₀) - te = vcat(incident(s, ∂₀, :∂v1)...) - se = vcat(incident(s, ∂₀, :∂v0)...) - unique(vcat(te, se)) -end - -function gen_dec_rules() - @present ExtendedOperators(FreeExtCalc2D) begin - X::Space - neg₁̃::Hom(DualForm1(X), DualForm1(X)) # negative - half₁̃::Hom(DualForm1(X), DualForm1(X)) # half - proj₁_¹⁰₁::Hom(Form1(X) ⊗ Form0(X), Form1(X)) - proj₂_¹⁰₀::Hom(Form1(X) ⊗ Form0(X), Form0(X)) - proj₁_¹²₁::Hom(Form1(X) ⊗ Form2(X), Form1(X)) - proj₂_¹²₂::Hom(Form1(X) ⊗ Form2(X), Form2(X)) - proj₁_¹²̃₁::Hom(Form1(X) ⊗ DualForm2(X), Form1(X)) - proj₂_¹²̃₂̃::Hom(Form1(X) ⊗ DualForm2(X), DualForm2(X)) - proj₁_¹¹̃₁::Hom(Form1(X) ⊗ DualForm1(X), Form1(X)) - proj₂_¹¹̃₁̃::Hom(Form1(X) ⊗ DualForm1(X), DualForm1(X)) - proj₁_¹̃¹̃₁̃::Hom(DualForm1(X) ⊗ DualForm1(X), DualForm1(X)) - proj₂_¹̃¹̃₁̃::Hom(DualForm1(X) ⊗ DualForm1(X), DualForm1(X)) - proj₁_¹¹₁::Hom(Form1(X)⊗Form1(X), Form1(X)) - proj₂_¹¹₁::Hom(Form1(X)⊗Form1(X), Form1(X)) - sum₁̃::Hom(DualForm1(X)⊗DualForm1(X), DualForm1(X)) - sum₁::Hom(Form1(X)⊗Form1(X), Form1(X)) - L₀::Hom(Form1(X)⊗DualForm2(X), DualForm2(X)) - L₁::Hom(Form1(X)⊗DualForm1(X), DualForm1(X)) - i₀::Hom(Form1(X)⊗DualForm2(X), DualForm1(X)) - i₁::Hom(Form1(X)⊗DualForm1(X), DualForm0(X)) - end - - Lie0Imp = @decapode ExtendedOperators begin - (dḞ2, dF2)::DualForm2{X} - F1::Form1{X} - - dḞ2 == dual_d₁{X}(i₀(F1, dF2)) - end - lie0_imp = diag2dwd(Lie0Imp, in_vars = [:F1, :dF2], out_vars = [:dḞ2]) - - Lie1Imp = @decapode ExtendedOperators begin - (dḞ1, dF1)::DualForm1{X} - F1::Form1{X} - dḞ1 == i₀(F1, dual_d₁{X}(dF1)) + dual_d₀{X}(i₁(F1, dF1)) - end - lie1_imp = diag2dwd(Lie1Imp, in_vars = [:F1, :dF1], out_vars = [:dḞ1]) - - I0Imp = @decapode ExtendedOperators begin - F1::Form1{X} - dF1::DualForm1{X} - dF2::DualForm2{X} - dF1 == ⋆₁{X}(∧₁₀{X}(F1, ⋆₀⁻¹{X}(dF2))) - end - i0_imp = diag2dwd(I0Imp, in_vars = [:F1, :dF2], out_vars = [:dF1]) - - I1Imp = @decapode ExtendedOperators begin - F1::Form1{X} - dF1::DualForm1{X} - dF0::DualForm0{X} - dF0 == ⋆₂{X}(∧₁₁{X}(F1, ⋆₁⁻¹{X}(dF1))) - end - i1_imp = diag2dwd(I1Imp, in_vars = [:F1, :dF1], out_vars = [:dF0]) - - δ₁Imp = @decapode ExtendedOperators begin - F1::Form1{X} - F0::Form0{X} - F0 == ⋆₀⁻¹{X}(dual_d₁{X}(⋆₁{X}(F1))) - end - δ₁_imp = diag2dwd(δ₁Imp, in_vars = [:F1], out_vars = [:F0]) - - δ₂Imp = @decapode ExtendedOperators begin - F1::Form1{X} - F2::Form2{X} - F1 == ⋆₁⁻¹{X}(dual_d₀{X}(⋆₂{X}(F2))) - end - δ₂_imp = diag2dwd(δ₂Imp, in_vars = [:F2], out_vars = [:F1]) - - Δ0Imp = @decapode ExtendedOperators begin - (F0, Ḟ0)::Form0{X} - Ḟ0 == δ₁{X}(d₀{X}(F0)) - end - Δ0_imp = diag2dwd(Δ0Imp, in_vars = [:F0], out_vars = [:Ḟ0]) - - Δ1Imp = @decapode ExtendedOperators begin - (F1, Ḟ1)::Form1{X} - Ḟ1 == δ₂{X}(d₁{X}(F1)) + d₀{X}(δ₁{X}(F1)) - end - Δ1_imp = diag2dwd(Δ1Imp, in_vars = [:F1], out_vars = [:Ḟ1]) - - Dict(:L₀ => lie0_imp, :i₀ => i0_imp, :L₁ => lie1_imp, :i₁ => i1_imp, - :δ₁ => δ₁_imp, :δ₂ => δ₂_imp, :Δ₀ => Δ0_imp, :Δ₁ => Δ1_imp) -end -end diff --git a/src/OpenDiagrams.jl b/src/OpenDiagrams.jl deleted file mode 100644 index e355553..0000000 --- a/src/OpenDiagrams.jl +++ /dev/null @@ -1,166 +0,0 @@ -module OpenDiagrams -export OpenDiagram, oapply, draw_diagram - -using Catlab, Catlab.CategoricalAlgebra, Catlab.WiringDiagrams, Catlab.Programs -import Catlab.CategoricalAlgebra: apex, legs, feet -import Catlab.WiringDiagrams: oapply -using Catlab.Programs.DiagrammaticPrograms: NamedGraph # FIXME: Should export? -using Catlab.Graphics - -using CombinatorialSpaces.ExteriorCalculus: ThExtCalc2D -using CombinatorialSpaces -using Catlab.CategoricalAlgebra.Categories: Functor, Cat -using Catlab.CategoricalAlgebra.FinCats: FinCatSize, FinCatPresentation - -const Diagram2D = Functor{Dom, Codom} where - {Ob, Hom, Dom<:(Cat{Ob, Hom, FinCatSize} where - {Ob, Hom}), Codom<:FinCatPresentation{ThExtCalc2D, Ob, Hom}} - -""" Open diagram as a structured multicospan in R-form. - -An open diagram is a diagram, represented as a `FinDomFunctor`, together with a -cospan of finite sets whose apex is the object set of the diagram's indexing -category. -""" -struct OpenDiagram{F<:FinDomFunctor,Cosp<:Multicospan} - functor::F - cospan::Cosp -end - -apex(diagram::OpenDiagram) = diagram.functor -legs(diagram::OpenDiagram) = legs(diagram.cospan) -feet(diagram::OpenDiagram) = feet(diagram.cospan) - -function OpenDiagram(F::FinDomFunctor, names::AbstractVector{Symbol}) - g = graph(dom(F)) - legs = map(name -> FinFunction([incident(g, name, :vname)], nv(g)), names) - OpenDiagram(F, Multicospan(FinSet(1), legs)) -end - -function OpenDiagram(d::Diagram, args...) - OpenDiagram(diagram(d), args...) -end - -const OpenFreeDiagramOb, OpenFreeDiagram = OpenACSetTypes(FreeDiagram, :V) - -OpenFreeDiagram(d::FreeDiagram{Ob,Hom}, args...) where {Ob,Hom} = - OpenFreeDiagram{Ob,Hom}(d, args...) - -function oapply(uwd::RelationDiagram, diagrams::AbstractVector{<:OpenDiagram}) - open_free_diagram = oapply(uwd, map(diagrams) do d - OpenFreeDiagram(FreeDiagram(apex(d)), d.cospan) - end) - free_diagram = apex(open_free_diagram) - leg_functions = map(leg -> leg[:V], legs(open_free_diagram)) - - g = NamedGraph{Symbol,Union{Symbol,Nothing}}() - copy_parts!(g, free_diagram) - g[:vname] = collect(apex(colimit_vertex_names(uwd, diagrams))) - g[:ename] = nothing # FIXME: Coproduct edge names, allowing for nulls. - - diagram = FinFunctor(ob(free_diagram), hom(free_diagram), - FinCat(g), codom(apex(first(diagrams)))) - OpenDiagram(diagram, Multicospan(leg_functions)) -end - -function colimit_vertex_names(uwd::RelationDiagram, - diagrams::AbstractVector{<:OpenDiagram}) - @assert nboxes(uwd) == length(diagrams) - bfd = BipartiteFreeDiagram{FinSet{<:Any,Symbol},FinFunction}() - add_vertices₁!(bfd, njunctions(uwd), ob₁=map(uwd[:variable]) do var - FinSet([var]) - end) - add_vertices₂!(bfd, nboxes(uwd), ob₂=map(diagrams) do diagram - FinSet(graph(dom(diagram.functor))[:vname]) - end) - for (b, diagram) in zip(boxes(uwd), diagrams) - g = graph(dom(diagram.functor)) - for (p, leg) in zip(ports(uwd, b), legs(diagram.cospan)) - j, v = junction(uwd, p), only(collect(leg)) - f = FinFunction(Dict(uwd[j, :variable] => g[v, :vname]), ob₂(bfd, b)) - add_edge!(bfd, j, b, hom=f) - end - end - colimit(bfd) -end - -#################### -# Graphics Tooling # -#################### - -draw_diagram(d::Diagram2D) = to_graphviz(d, node_labels=true, - prog="neato", - node_attrs=Dict(:shape=>"oval"), - graph_attrs=Dict(:nodesep=>"4.0")) -using Catlab.Graphs -using Catlab.Graphs.BasicGraphs -function to_graph(J::FinCat) - g = BasicGraphs.Graph() - copy_parts!(g, graph(J)) - return g -end - -Graphics.to_graphviz(F::Diagram2D; kw...) = -to_graphviz(GraphvizGraphs.to_graphviz_property_graph(F; kw...)) - -function GraphvizGraphs.to_graphviz_property_graph(F::Diagram2D; kw...) - simplify_vname(s) = begin - table = Dict( - "Form0(X)" => "Ω₀", - "Form1(X)" => "Ω₁", - "Form2(X)" => "Ω₂", - "DualForm0(X)" => "Ω̃₀", - "DualForm1(X)" => "Ω̃₁", - "DualForm2(X)" => "Ω̃₂" - ) - st = string(s) - m = match(r"otimes\((.*)\)", st) - if isnothing(m) - return table[st] - else - forms = only(m.captures) - join(map(f -> table[f], split(forms, ",")), "×") - end - end - - simplify_ename(s) = begin - b = IOBuffer() - show_unicode(b, s) - return replace(String(take!(b)), r"{X}"=>"") - end - - J = dom(F) - G = graph(J) - pg = GraphvizGraphs.to_graphviz_property_graph(to_graph(J); kw...) - for v in vertices(G) - lᵥ = G[v, :vname] - tᵥ = simplify_vname(F.ob_map[v]) - set_vprops!(pg, v, Dict(:label => "$(lᵥ):$tᵥ")) - end - for e in edges(G) - tₑ = F.hom_map[e] - set_eprops!(pg, e, Dict(:label => "$(simplify_ename(tₑ))")) - end - pg -end - - -Graphics.to_graphviz(cosp::OpenDiagram; kw...) = -to_graphviz(GraphvizGraphs.to_graphviz_property_graph(cosp; kw...)) - - -function GraphvizGraphs.to_graphviz_property_graph(cosp::OpenDiagram; kw...) - pg = GraphvizGraphs.to_graphviz_property_graph(cosp.functor; kw...) - label(I, l, i) = length(dom(l)) > 1 ? "$(i):$I" : "$I" - for (I,l) in enumerate(legs(cosp.cospan)) - for i in dom(l) - v = add_vertex!(pg) - set_vprops!(pg, v, Dict(:label=>label(I, l,i), :shape=>"circle", :color=>"blue")) - e = add_edge!(pg, v, l(i)) - set_eprops!(pg, e, Dict(:style=>"dashed")) - end - end - return pg -end - -end diff --git a/src/PetriNets.jl b/src/PetriNets.jl deleted file mode 100644 index 8f8c44e..0000000 --- a/src/PetriNets.jl +++ /dev/null @@ -1,177 +0,0 @@ -""" Generating Decapodes from Petri nets - -This module provides support for generating a simulatable Decapode from Petri -nets generated in the AlgebraicPetri library. -""" - -module PetriNets - -using Catlab, Catlab.CategoricalAlgebra, Catlab.Graphics, Catlab.Programs -using CombinatorialSpaces.ExteriorCalculus -using Catlab.Programs.DiagrammaticPrograms: NamedGraph -using ..Simulations -using ..AlgebraicPetri -using LinearAlgebra -using Unicode - -export pn2dec, expand_pres!, gen_functions - -i2sub = Dict( - '0'=>'₀', '1'=>'₁', '2'=>'₂', '3'=>'₃', '4'=>'₄', '5'=>'₅', - '6'=>'₆', '7'=>'₇', '8'=>'₈', '9'=>'₉', '-'=>'₋' -) -i2sup = Dict( - '0'=>'⁰', '1'=>'¹', '2'=>'²', '3'=>'³', '4'=>'⁴', '5'=>'⁵', - '6'=>'⁶', '7'=>'⁷', '8'=>'⁸', '9'=>'⁹', '-'=>'⁻' -) - -function neg_term(graph, obs, homs, prs, val::Int64) - res = add_part!(graph, :V, vname=Symbol("$(gensym())")) - push!(obs, prs.syntax.Form0(prs[:X])) - add_part!(graph, :E, ename=nothing, src=val, tgt=res) - push!(homs, prs[Symbol("neg₀")]) - return res -end - -function scale_term(graph, obs, homs, prs, val::Int64, sc::Int64) - res = add_part!(graph, :V, vname="$(gensym())") - push!(obs, prs.syntax.Form0(prs[:X])) - add_part!(graph, :E, ename=nothing, src=val, tgt=res) - push!(homs, prs[Symbol("mult_$(sc)")]) - return res -end - -function exp_term(graph, obs, homs, prs, val::Int64, sc::Int64) - res = add_part!(graph, :V, vname=Symbol("$(gensym())")) - push!(obs, prs.syntax.Form0(prs[:X])) - add_part!(graph, :E, ename=nothing, src=val, tgt=res) - push!(homs, prs[Symbol("^$(sc)")]) - return res -end - -function rate_term(graph, obs, homs, prs, val::Int64, tr_name::Symbol, tgt::Int64) - add_part!(graph, :E, ename=nothing, src=val, tgt=tgt) - push!(homs, prs[Symbol("k_", tr_name)]) - return tgt -end - -function rate_term(graph, obs, homs, prs, val::Int64, tr_name::Symbol) - res = add_part!(graph, :V, vname=Symbol("$(gensym())")) - push!(obs, prs.syntax.Form0(prs[:X])) - rate_term(graph, obs, homs, prs, val, tr_name, res) -end - -function sum_dec!(graph, obs, homs, prs, vals::Int64...) - bin_exp!(graph, obs, homs, prs, :sum₀, vals...) -end - -function prod_dec!(graph, obs, homs, prs, vals::Int64...) - bin_exp!(graph, obs, homs, prs, :prod₀, vals...) -end - -function bin_exp!(graph, obs, homs, prs, exp, vals::Int64...) - if length(vals) == 1 - first(vals) - else - form0 = prs.syntax.Form0(prs[:X]) - form0x0 = prs.syntax.otimes(form0, form0) - cur_val = first(vals) - bin_int = add_parts!(graph, :V, length(vals)-1, - vname=[Symbol("$(gensym())") for i in 1:(length(vals)-1)]) - bin_res = add_parts!(graph, :V, length(vals)-1, - vname=[Symbol("$(gensym())") for i in 1:(length(vals)-1)]) - append!(obs, fill(form0x0, length(vals)-1)) - append!(obs, fill(form0, length(vals)-1)) - for i in 1:(length(vals)-1) - add_parts!(graph, :E, 3, - ename =[nothing, nothing, nothing], - tgt=[vals[i+1], cur_val, bin_res[i]], - src=[bin_int[i], bin_int[i], bin_int[i]]) - append!(homs, [prs[:proj₁_⁰⁰₀], prs[:proj₂_⁰⁰₀], prs[exp]]) - cur_val = bin_res[i] - end - cur_val - end -end - -""" pn2dec(prs::Presentation, pn::LabelledReactionNet) - -Generates a Decapode diagram which represents the law of mass action applied to -the petri net `pn` with the syntax from the presentation `prs`. -""" -function pn2dec(prs, pn::LabelledReactionNet) - synt = prs.syntax - - # Initial state variables - graph = NamedGraph{Symbol,Union{Symbol,Nothing}}() - svars = add_parts!(graph, :V, ns(pn), vname=snames(pn)) - tvars = add_parts!(graph, :V, nt(pn), vname=tnames(pn)) - obs = Vector{Any}(fill(synt.Form0(prs[:X]), ns(pn) + nt(pn))) #"($(join(vcat(snames(pn), tnames(pn)), ", ")))::Form0{X}" - homs = Vector{Any}() - map(1:nt(pn)) do t - inputs = pn[incident(pn, t, :it), :is] - res = prod_dec!(graph, obs, homs, prs, svars[inputs]...) - rate_term(graph, obs, homs, prs, res, tname(pn, t), tvars[t]) - end - - res_s = map(1:ns(pn)) do s - inps = tvars[pn[incident(pn, s, :os), :ot]] - otps = tvars[pn[incident(pn, s, :is), :it]] - otp_scaled = [neg_term(graph, obs, homs, prs, o) for o in otps] - if isempty(vcat(inps, otp_scaled)) - res = add_part!(graph, :V, vname=Symbol("$(gensym())")) - push!(obs, prs.syntax.Form0(prs[:X])) - res - else - sum_dec!(graph, obs, homs, prs, vcat(inps, otp_scaled)...) - end - end - - set_subparts!(graph, res_s, - vname=[length("$name") > 1 ? - Symbol("∂ₜ", "$name") : Symbol(Unicode.normalize("$name"* '̇')) - for name in snames(pn)]) - FinFunctor(obs, homs, FinCat(graph), FinCat(prs)) -end - -""" expand_pres!(pres::Presentation, pn::LabelledReactionNet) - -Expands the syntax of the presentation `pres` to include the necessary operators for expressing mass action on the Petri net `pn`. -""" -function expand_pres!(pres, pn::LabelledReactionNet) - synt = pres.syntax - form0 = synt.Form0(pres[:X]) - form0x0 = synt.otimes(form0, form0) - gens = Dict(:neg₀ => synt.Hom(:neg₀, form0, form0), - :sum₀ => synt.Hom(:sum₀, form0x0, form0), - :prod₀ => synt.Hom(:prod₀, form0x0, form0), - :proj₁_⁰⁰₀ => synt.Hom(:proj₁_⁰⁰₀, form0x0, form0), - :proj₂_⁰⁰₀ => synt.Hom(:proj₂_⁰⁰₀, form0x0, form0)) - for t in tnames(pn) - gens[t] = synt.Hom(Symbol("k_$t"), form0, form0) - end - - for (k,v) in gens - if !has_generator(pres, k) - add_generator!(pres, v) - end - end -end - -""" gen_functions(pn::LabelledReactionNet, s::EmbeddedDeltaDualComplex2D) - -Generates the functions necessary for executing the Petri net after it has been -converted to a Decapode. -""" -function gen_functions(pn::LabelledReactionNet, s) - funcs = Dict{Symbol, Any}() - funcs[:sum₀] = Dict(:operator => (s,a,b) -> (s .= a .+ b), :type => InPlaceFunc()) - funcs[:prod₀] = Dict(:operator => (s,a,b) -> (s .= a .* b), :type => InPlaceFunc()) - funcs[:neg₀] = Dict(:operator => (s,a) -> (s .= -1 * a), :type => InPlaceFunc()) - for t in 1:nt(pn) - funcs[Symbol("k_", tname(pn, t))] = Dict(:operator => rate(pn, t) * I(nv(s)), :type => MatrixFunc()) - end - funcs -end - -end diff --git a/src/Schedules.jl b/src/Schedules.jl deleted file mode 100644 index 7038716..0000000 --- a/src/Schedules.jl +++ /dev/null @@ -1,245 +0,0 @@ -""" Schedules - -This module contains the tooling which converts diagrams to DWDs. - -*TODO* This module will also include optimization tooling relating to DWD -schedules (compressing matrix operations, pre-computing constant values, etc). -""" -module Schedules -using Catlab.Syntax -using Catlab.WiringDiagrams.DirectedWiringDiagrams -using Catlab.Theories -using Catlab.CategoricalAlgebra -using Catlab.Programs.DiagrammaticPrograms: NamedGraph - -using ..Diagrams - -export diag2dwd - -sp_otimes(expr) = expr isa ObExpr{:otimes} ? expr.args : [expr] -n_args(expr) = length(sp_otimes(expr)) - -isa_proj(hom) = startswith(string(hom), "proj") -proj_ind(hom) = unsub(split(string(hom), "_")[1][5:end]) - -super2int = Dict( - '⁰'=>'0', '¹'=>'1', '²'=>'2', '³'=>'3', '⁴'=>'4', '⁵'=>'5', - '⁶'=>'6', '⁷'=>'7', '⁸'=>'8', '⁹'=>'9' -) -sub2int = Dict( - '₀'=>'0', '₁'=>'1', '₂'=>'2', '₃'=>'3', '₄'=>'4', '₅'=>'5', - '₆'=>'6', '₇'=>'7', '₈'=>'8', '₉'=>'9' -) - -unsub(str) = begin - parse(Int64, join([sub2int[a] for a in str])) -end - -function eval_deps!(dwd, graph, el, w2b, el2p, obs; in_els = Dict{Int, Int}(), boundaries = []) - if el in keys(el2p) - return el2p[el] - end - - inputs = incident(graph, el, :tgt) - - labels = graph[:ename] - el_types = sp_otimes(obs[el]) - outp = if length(inputs) > 1 - sort!(inputs, by = i->proj_ind(labels[i])) - vcat(map(w -> eval_deps!(dwd, graph, graph[w, :src], w2b, el2p, obs; in_els = in_els), inputs)...) - elseif length(inputs) == 0 - el ∈ keys(in_els) || error("Element $el has no dependencies, but is not defined in `in_els`") - [(input_id(dwd), in_els[el])] - else - in_ps = eval_deps!(dwd, graph, graph[inputs[1], :src], w2b, el2p, obs; in_els = in_els) - if isa_proj(graph[inputs[1], :ename]) - [in_ps[proj_ind(graph[inputs[1], :ename])]] - else - wires = map(enumerate(in_ps)) do (ip, op) - Wire(port_value(dwd, Port(op[1], OutputPort, op[2])), op, (w2b[inputs[1]], ip)) - end - add_wires!(dwd, wires) - map(1:length(output_ports(dwd, w2b[inputs[1]]))) do p - (w2b[inputs[1]], p) - end - end - end - - if el in keys(el2p) - return el2p[el] - end - - # Evaluate Boundary Conditions - bcs = incident(graph, el, :src) - filter!(a->"$(Symbol(graph[a, :ename]))"[1] == '∂', bcs) - for bc in bcs - wires = map(enumerate(outp)) do (ip, op) - Wire(port_value(dwd, Port(op[1], OutputPort, op[2])), op, (w2b[bc], ip)) - end - add_wires!(dwd, wires) - outp .= map(1:length(output_ports(dwd, w2b[bc]))) do p - (w2b[bc], p) - end - end - el2p[el] = outp - outp -end - -name(a::HomExpr) = head(a) == :generator ? args(a)[1] : head(a) - -""" diag2dwd(diagram; clean = false, calc_states = [], out_vars=[], in_vars=[]) - -Generates a directed wiring diagram (DWD) from a Decapode `diagram`. This -method generates a DWD which represents a single explicit time step of the -system, computed by the dependency graph of each time derivative within the -system. This does not necessarily work for all Decapodes, and only generates -an explicit time-step solution. - -*TODO*: Functions like this can be created to generate other solution styles, -like implicit solutions. -""" -function diag2dwd(diagram; clean = false, calc_states = [], out_vars=[], in_vars=[]) - homs = Vector{Any}(copy(diagram.hom_map)) - obs = Vector{Any}(copy(diagram.ob_map)) - graph = NamedGraph{Any, Any}() - copy_parts!(graph, dom(diagram).graph) - graph[:ename] .= homs - - - flipped = fill(false, ne(graph)) - # Flip projections for graph eval - for (i, h) in enumerate(homs) - if isa_proj(h) && !flipped[i] - iw = copy(incident(graph, graph[i, :src], :tgt)) - ow = copy(incident(graph, graph[i, :src], :src)) - if isempty(iw) - for w in ow - if isa_proj(homs[w]) && !flipped[w] - wsrc = graph[w, :src] - graph[w, :src] = graph[w, :tgt] - graph[w, :tgt] = wsrc - flipped[w] = true - end - end - end - end - end - - # Expand homs which are composition - comp_to_remove = [] - for h in parts(graph, :E) - h_name = graph[h, :ename] - elsrc, eltgt = graph[h, :src], graph[h, :tgt] - if h_name isa HomExpr{:compose} - args = h_name.args - push!(comp_to_remove, h) - verts = add_parts!(graph, :V, length(args) - 1, vname = [Symbol("anon_", gensym()) for i in 1:(length(args)-1)]) - append!(obs, codom.(args)[1:(end-1)]) - add_parts!(graph, :E, length(args), ename = args, - src = vcat([elsrc], verts), - tgt = vcat(verts, [eltgt])) - end - end - rem_parts!(graph, :E, comp_to_remove) - - - pres = presentation(codom(diagram)) - - # Decide on State and Parameter variables - time_arrs = findall(h-> h isa HomExpr{:∂ₜ}, graph[:ename]) - params = findall(e -> isempty(incident(graph, e, :tgt)), parts(graph, :V)) - state_vals = graph[time_arrs, :src] - - # Mix State and Parameters into single in_els (need to ensure alignment between input and outputs - in_els = unique(vcat(state_vals, params)) - out_els = copy(graph[time_arrs, :tgt]) - out_names = [Symbol(:∂ₜ, graph[v, :vname]) for v in state_vals] - - - # TODO: Establish some kind of consistency expectations here - if !isempty(in_vars) && !isempty(out_vars) - in_els = [incident(graph, a, :vname) for a in in_vars] - out_els = [incident(graph, a, :vname) for a in out_vars] - out_names = graph[out_els, :vname] - elseif !isempty(in_vars) - # This ensures that the output values are in-order with input values - new_inds = [findfirst(v -> graph[v, :vname] == a, in_els) for a in in_vars] - non_param_inds = filter(i -> in_els[i] in state_vals , new_inds) - in_els .= in_els[new_inds] - out_els .= out_els[non_param_inds] - out_names .= out_names[non_param_inds] - elseif !isempty(out_vars) - out_els = [incident(graph, a, :vname) for a in out_vars] - out_names = graph[out_els, out_els] - end - if !isempty(calc_states) - calc_inds = findall(v->graph[v, :vname] ∈ calc_states, state_vals) - out_els = out_els[calc_inds] - out_names = out_names[calc_inds] - end - - in_types = [Dict(:name => graph[v, :vname], - :type => obs[v]) - for v in in_els] - out_types = [Dict(:name => out_names[v], - :type => obs[out_els[v]]) - for v in 1:length(out_els)] - - # FIXME: Hacky solution to ensuring time_arrs aren't included later - # need more thoughtful approach to sorting through edges here - rem_parts!(graph, :E, time_arrs) - - dwd = WiringDiagram(collect(in_types), collect(out_types)) - - # Add all necessary boxes to the DWD based on arrows in the diagram - w2b = map(parts(graph, :E)) do a - w_type = graph[a, :ename] - # Fix this if-case. This was meant for times when multiple arguments would - # go into a single box - if isa_proj(w_type) - nothing - elseif w_type isa HomExpr - src_el = graph[a, :src] - tgt_el = graph[a, :tgt] - in_ports = sp_otimes(obs[src_el]) - out_ports = sp_otimes(obs[tgt_el]) - add_box!(dwd, Box(name(w_type), [Dict(:type=>ip) for ip in in_ports], - [Dict(:type=>op) for op in out_ports])) - end - end - - # Mapping from elements (vertices in `graph`) to ports in the DWD - # This ensures that elements which have already been computed are - # not recomputed - el2p = Dict{Int, Vector}() - in_ps = vcat(map(out_els) do el - eval_deps!(dwd, graph, el, w2b, el2p, obs; in_els=Dict(in_els[i] => i for i in 1:length(in_els))) - end...) - - wires = map(enumerate(in_ps)) do (ip, op) - Wire(out_types[ip], op, (output_id(dwd), ip)) - end - add_wires!(dwd, wires) - - # Remove any boxes without any connecting wires - if clean - to_rem = Vector{Int64}() - for b in 1:nparts(dwd.diagram, :Box) - if isempty(vcat(incident(dwd.diagram, b, [:tgt, :in_port_box]), - incident(dwd.diagram, b, [:in_tgt, :in_port_box]))) && - isempty(vcat(incident(dwd.diagram, b, [:src, :out_port_box]), - incident(dwd.diagram, b, [:out_src, :out_port_box]))) - push!(to_rem, b) - end - end - rem_boxes!(dwd, to_rem) - end - - dwd -end - -diag2dwd(simple_diagram::T; args...) where {T <: SimpleDiagram} = - diag2dwd(diagram(simple_diagram); args...) - - -end diff --git a/src/Simulations.jl b/src/Simulations.jl deleted file mode 100644 index 9e80193..0000000 --- a/src/Simulations.jl +++ /dev/null @@ -1,289 +0,0 @@ -""" Generating simulation functions from directed wiring diagrams - -This module currently provides support for generating a simple explicit -time-stepping function based on a directed wiring diagram. The generated -function is compatible with the DifferentialEquations.jl interface, which is -used in demo execution. -""" - -module Simulations - -using CombinatorialSpaces -using Catlab -using Catlab.WiringDiagrams -using Catlab.CategoricalAlgebra -using Catlab.Theories -using PreallocationTools -using LinearAlgebra - -export gen_sim, - BoxFunc, MatrixFunc, ElementwiseFunc, ArbitraryFunc, InPlaceFunc, - ConstantFunc, TDInPlaceFunc - -abstract type BoxFunc end - -""" MatrixFunc - -Linear operator applied through matrix multiplication. -""" -struct MatrixFunc <: BoxFunc end - -""" ElementwiseFunc - -Function applied to each element of one or more vectors. -""" -struct ElementwiseFunc <: BoxFunc end - -""" ArbitraryFunc - -Out of place function which applies to any number of arguments. -""" -struct ArbitraryFunc <: BoxFunc end - -""" ConstantFunc - -Function which always returns the same value. -""" -struct ConstantFunc <: BoxFunc end - -""" InPlaceFunc - -Function which operates on arguments in-place, not returning a value. -""" -struct InPlaceFunc <: BoxFunc end - -""" TDInPlaceFunc - -Function which accepts the current time as the last argument. Used to set up -time-dependent physics or boundary conditions. -""" -struct TDInPlaceFunc <: BoxFunc end - -form2dim_def = Dict(:Scalar => x->1, - :Form0 => nv, - :Form1 => ne, - :Form2 => ntriangles, - :DualForm2 => nv, - :DualForm1 => ne, - :DualForm0 => ntriangles) -dims(x, form2dim) = begin - k = filter(i -> x[:type] isa ObExpr{i}, keys(form2dim)) - length(k) == 1 || error("Object $x has multiple keys in `form2dim`") - form2dim[first(k)] -end - -""" gen_sim(dwd::WiringDiagram, name2func::Dict, - s::EmbeddedDeltaDualComplex2D; - form2dim=form2dim, params=[], autodiff=false) - -This function generates a function which evaluates the acylcic directed wiring -diagram `dwd`, using the functions in `name2func` for operators to execute -corresonding to each box and information from the mesh `s` to pre-allocate -necessary memory. This operator can generate a function which is compatible -with the autodifferentiation solvers in DifferentialEquations.jl. -""" -function gen_sim(dwd::WiringDiagram, name2func, s; form2dim=form2dim_def, params=[], autodiff=false) - check_consistency(dwd) - d = dwd.diagram - - # Generate cached memory. These variables are indexed by their respective - # output ports - mem = [zeros(Float64, dims(f, form2dim)(s)) for f in d[:out_port_type]] - if autodiff - mem = [dualcache(zeros(Float64, dims(f, form2dim)(s))) for f in d[:out_port_type]] - end - tgt2src = Dict{Int64, Int64}() - for w in 1:nparts(d, :Wire) - d[w, :tgt] ∈ keys(tgt2src) && error("Two wires input to port $(d[w, :src])") - tgt2src[d[w, :tgt]] = d[w, :src] - end - - # Get the indices in input/output arrays for specific input/output forms - input_inds = Vector{Tuple{Int64, Int64, Bool}}() - output_inds = Vector{Tuple{Int64, Int64}}() - - cur_ind = 0 - cur_param = 0 - for ft in d[:outer_in_port_type] - mag = dims(ft, form2dim)(s) - if ft[:name] ∈ params - push!(input_inds, (cur_param+1, mag + cur_param, true)) - cur_param += mag - else - push!(input_inds, (cur_ind+1, mag + cur_ind, false)) - cur_ind += mag - end - end - cur_ind = 0 - for ft in d[:outer_out_port_type] - mag = dims(ft, form2dim)(s) - push!(output_inds, (cur_ind+1, mag + cur_ind)) - cur_ind += mag - end - - # Make local variables for any MatrixFunc/ArbitraryFunc functions to be - # included in the function closure - n2f = deepcopy(name2func) - matrices = Vector{Any}() - funcs = Vector{Function}() - - # First develop a list of needed functions/matrices - used_funcs = unique(d[:value]) - - for k in used_funcs - ops = split("$k", "⋅") - if length(ops) > 1 - mat_val = foldl(*, map(o->name2func[Symbol(o)][:operator], ops)) - push!(matrices, mat_val) - n2f[k] = Dict(:operator => :(matrices[$(length(matrices))]), :type => MatrixFunc()) - else - v = name2func[k] - if(v[:type] isa MatrixFunc) - push!(matrices, v[:operator]) - n2f[k] = Dict(:operator => :(matrices[$(length(matrices))]), :type => v[:type]) - elseif(v[:type] isa Union{ElementwiseFunc, ArbitraryFunc, InPlaceFunc, TDInPlaceFunc}) - push!(funcs, v[:operator]) - n2f[k] = Dict(:operator => :(funcs[$(length(funcs))]),:type => v[:type]) - end - end - end - - # Compute composed matrices - for v in d[:value] - - end - - # Topological sort of DWD for scheduling - execution_order = topological_sort(internal_graph(dwd)) - - # Generate function expression - body = quote - end - - # Assign variables for each dual memory location - # Allows for get_tmp to be called for each cached memory location - if autodiff - dual_init = map(1:length(mem)) do i - :($(Symbol("d_$i")) = get_tmp(mem[$i], u)) - end - append!(body.args, dual_init) - end - - exec_dwd = map(execution_order) do b - iw = in_wires(dwd, b) - - input_args = map(iw) do w - p = w.source - if p.box == -2 - if input_inds[p.port][3] - :(view(p, $(input_inds[p.port][1]):$(input_inds[p.port][2]))) - else - :(view(u, $(input_inds[p.port][1]):$(input_inds[p.port][2]))) - end - else - if autodiff - :($(Symbol("d_$(incident(d, p.box, :out_port_box)[p.port])"))) - else - :(mem[$(incident(d, p.box, :out_port_box)[p.port])]) - end - end - end - input_args[[w.target.port for w in iw]] .= input_args - - output_args = map(incident(d, b, :out_port_box)) do p - if autodiff - :($(Symbol("d_$p"))) - else - :(mem[$p]) - end - end - - gen_func(input_args, output_args, n2f[d[b, :value]]) - end - append!(body.args, exec_dwd) - - # set du values - du_set = map(1:nparts(d, :OuterOutPort)) do i - iw = incident(d, i, :out_tgt) - length(iw) == 1 || error(Output) - w = first(iw) - pt = d[w, :out_tgt] - ps = d[w, :out_src] - if autodiff - :(setindex!(du, Symbol("d_$ps"), $(collect(output_inds[pt][1]:output_inds[pt][2])))) - else - :(setindex!(du, mem[$ps], $(collect(output_inds[pt][1]:output_inds[pt][2])))) - end - end - append!(body.args, du_set) - - # Currently this generates a function at the module (can be accessed from the module) - # Is this alright, or is this going to just lead to memory problems? - fname = gensym() - res = :($fname(du, u, p, t, mem, matrices, funcs) = $body) - eval(res) - f_local = deepcopy(eval(fname)) - function sim_func(du, u, p, t) - f_local(du, u, p, t, mem, matrices, funcs) - end, res -end - -function gen_func(input_args, output_args, func_info) - _gen_func(func_info[:type], input_args, output_args, func_info[:operator]) -end - -function _gen_func(::MatrixFunc, input_args, output_args, func) - length(input_args) == 1 || error("Length of input for Matrix function is $(length(input_args)) ≠ 1") - length(output_args) == 1 || error("Length of output for Matrix function is $(length(output_args)) ≠ 1") - :(mul!($(output_args[1]), $(func), $(input_args[1]))) -end - -function _gen_func(::ArbitraryFunc, input_args, output_args, func) - - length(output_args) == 0 && error("Length of output for Arbitrary function is $(length(output_args))") - :($(Meta.parse(join(output_args, ","))) = $(Expr(:call, func, input_args...))) -end - -function _gen_func(::InPlaceFunc, input_args, output_args, func) - length(output_args) == 0 && error("Length of output for Arbitrary function is $(length(output_args))") - :($(Expr(:call, func, vcat(output_args, input_args)...))) -end - -function _gen_func(::TDInPlaceFunc, input_args, output_args, func) - _gen_func(InPlaceFunc(), vcat(input_args, [:t]), output_args, func) -end - -function _gen_func(::ConstantFunc, input_args, output_args, func) - length(output_args) == 0 && error("Length of output for Constant function is $(length(output_args))") - :($(Meta.parse(join(output_args, ","))) .= $(func)) -end - - -function _gen_func(::ElementwiseFunc, input_args, output_args, func) - - length(output_args) == 0 && error("Length of output for Arbitrary function is $(length(output_args))") - body = :($(Meta.parse(join([Expr(:ref, o, :i) for o in output_args], ","))) = $(Expr(:call, func, [Expr(:ref, n, :i) for n in input_args]...))) - quote - for i in eachindex($(input_args[1])) - $body - end - end -end - -function check_consistency(dwd::WiringDiagram) - d = dwd.diagram - for w in 1:nparts(d, :Wire) - sp = d[w, :src] - tp = d[w, :tgt] - stype = d[sp, :out_port_type] - ttype = d[tp, :in_port_type] - sb = d[sp, :out_port_box] - tb = d[tp, :in_port_box] - - stype == ttype || - error("""Wire $w maps a port of type $stype to a port of type $ttype - This wire is between box $(sb)($(d[sb, :value])) and box $(tb)($(d[tb, :value]))""") - end -end - -end diff --git a/src/Decapodes2/composition.jl b/src/composition.jl similarity index 100% rename from src/Decapodes2/composition.jl rename to src/composition.jl diff --git a/src/Decapodes2/coordinates.jl b/src/coordinates.jl similarity index 100% rename from src/Decapodes2/coordinates.jl rename to src/coordinates.jl diff --git a/src/Decapodes2/decapodeacset.jl b/src/decapodeacset.jl similarity index 100% rename from src/Decapodes2/decapodeacset.jl rename to src/decapodeacset.jl diff --git a/src/Decapodes2/language.jl b/src/language.jl similarity index 100% rename from src/Decapodes2/language.jl rename to src/language.jl diff --git a/src/Decapodes2/meshes.jl b/src/meshes.jl similarity index 100% rename from src/Decapodes2/meshes.jl rename to src/meshes.jl diff --git a/src/Decapodes2/rewrite.jl b/src/rewrite.jl similarity index 100% rename from src/Decapodes2/rewrite.jl rename to src/rewrite.jl diff --git a/src/Decapodes2/simulation.jl b/src/simulation.jl similarity index 100% rename from src/Decapodes2/simulation.jl rename to src/simulation.jl diff --git a/src/Decapodes2/visualization.jl b/src/visualization.jl similarity index 100% rename from src/Decapodes2/visualization.jl rename to src/visualization.jl diff --git a/test/Decapodes.jl b/test/Decapodes.jl deleted file mode 100644 index be75d0f..0000000 --- a/test/Decapodes.jl +++ /dev/null @@ -1,82 +0,0 @@ -module DecapodesTest - - using Catlab - using Catlab.Present - using Catlab.Programs - using CombinatorialSpaces - using CombinatorialSpaces.ExteriorCalculus - using GeometryBasics - using LinearAlgebra - using Test - - using Decapodes.Simulations - using Decapodes.Diagrams - using Decapodes.Schedules - using Decapodes.OpenDiagrams - using Decapodes.Examples - using Random - - @present DiffusionSpace2D(FreeExtCalc2D) begin - X::Space - k::Hom(Form1(X), Form1(X)) # diffusivity of space, usually constant (scalar multiplication) - proj¹_⁰⁰₀::Hom(Form0(X) ⊗ Form0(X), Form0(X)) - proj²_⁰⁰₀::Hom(Form0(X) ⊗ Form0(X), Form0(X)) - sum₀::Hom(Form0(X) ⊗ Form0(X), Form0(X)) - prod₀::Hom(Form0(X) ⊗ Form0(X), Form0(X)) - end - - Diffusion = @free_diagram DiffusionSpace2D begin - (C, Ċ)::Form0{X} - - Ċ == ⋆₀⁻¹{X}(dual_d₁{X}(⋆₁{X}(k(d₀{X}(C))))) - end - - Dynamics = @free_diagram DiffusionSpace2D begin - (Ċ, C)::Form0{X} - - Ċ == ∂ₜ{Form0{X}}(C) - end; - - compose_diff = @relation (C,) begin - diff(C, Ċ) - dynamics(C, Ċ) - end - - composed_diff = oapply(compose_diff, [OpenDiagram(Diffusion, [:C, :Ċ]), OpenDiagram(Dynamics, [:C, :Ċ])]); - - res = diag2dwd(composed_diff.functor, in_vars = [:C], out_vars = [:Ċ]) - - s = EmbeddedDeltaSet2D{Bool, Point{3,Float64}}() - points = [(0,0,0),(0,0,1),(0,1,0),(0,1,1)] - - add_vertices!(s, 4, point=points) - glue_sorted_triangle!(s, 1,2,3) - glue_sorted_triangle!(s, 2,3,4) - s[:edge_orientation] = false - s[:tri_orientation] = false - orient!(s) - sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3{Float64}}(s) - subdivide_duals!(sd, Barycenter()) - #sd = dual(s); - - funcs = sym2func(sd) - funcs[:k] = Dict(:operator => 1.0 * I(ne(sd)), :type => MatrixFunc()) - Examples.contract_matrices!(res, funcs) - - - sim, _ = gen_sim(res, funcs, sd); - - Random.seed!(42) - u = [rand() for i in 1:nv(s)] - du = zeros(Float64,nv(s)) - - tempu = copy(u) - dt = 0.01 - for i in 1:1000 - sim(du, tempu, [],0) - tempu .+= du * dt - # Check that mass is conserved - end - @test all(tempu .!= u) - @test all(sum(⋆(Val{0}, sd)*tempu)-sum(⋆(Val{0}, sd)*u) < 1e-6) -end diff --git a/test/Decapodes2/runtests.jl b/test/Decapodes2/runtests.jl deleted file mode 100644 index 70e3da3..0000000 --- a/test/Decapodes2/runtests.jl +++ /dev/null @@ -1,37 +0,0 @@ -using Test - -@testset "Composition" begin - include(joinpath(@__DIR__, "composition.jl")) -end - -@testset "Coordinates" begin - include(joinpath(@__DIR__, "coordinates.jl")) -end - -@testset "Diag2DWD" begin - include(joinpath(@__DIR__, "diag2dwd.jl")) -end - -@testset "Meshes" begin - include(joinpath(@__DIR__, "meshes.jl")) -end - -@testset "MultiScaleArrays" begin - include(joinpath(@__DIR__, "multiscalearrays.jl")) -end - -@testset "Summation" begin - include(joinpath(@__DIR__, "summation.jl")) -end - -@testset "Simulation" begin - include(joinpath(@__DIR__, "simulation.jl")) -end - -@testset "Visualization" begin - include(joinpath(@__DIR__, "visualization.jl")) -end - -@testset "Rewrite" begin - include(joinpath(@__DIR__, "rewrite.jl")) -end \ No newline at end of file diff --git a/test/Examples.jl b/test/Examples.jl deleted file mode 100644 index d283fb6..0000000 --- a/test/Examples.jl +++ /dev/null @@ -1,18 +0,0 @@ -module ExamplesTest - - using Catlab - using Catlab.Present - using Catlab.Programs - using CombinatorialSpaces - using Catlab.WiringDiagrams - using CombinatorialSpaces.ExteriorCalculus - using Test - - using Decapodes.Examples - - rules = gen_dec_rules() - for v in values(rules) - @test v isa WiringDiagram - end - -end diff --git a/test/PetriNets.jl b/test/PetriNets.jl deleted file mode 100644 index ed95ef7..0000000 --- a/test/PetriNets.jl +++ /dev/null @@ -1,110 +0,0 @@ -module PetriNetsTest - - using Catlab - using Catlab.Present - using Catlab.Programs - using CombinatorialSpaces - using CombinatorialSpaces.ExteriorCalculus - using GeometryBasics - using LinearAlgebra - using AlgebraicPetri - using Test - - using Decapodes.Simulations - using Decapodes.Diagrams - using Decapodes.Schedules - using Decapodes.Examples - using Decapodes.OpenDiagrams - using Decapodes.PetriNets - using Random - @present DiffusionSpace2D(FreeExtCalc2D) begin - X::Space - k::Hom(Form1(X), Form1(X)) # diffusivity of space, usually constant (scalar multiplication) - proj₁_⁰⁰₀::Hom(Form0(X) ⊗ Form0(X), Form0(X)) - proj₂_⁰⁰₀::Hom(Form0(X) ⊗ Form0(X), Form0(X)) - sum₀::Hom(Form0(X) ⊗ Form0(X), Form0(X)) - prod₀::Hom(Form0(X) ⊗ Form0(X), Form0(X)) - end - - SIRD = LabelledReactionNet{Float64, Float64}([:S=>0.0, :I=>0.0, :R=>0.0, :D=>0.0], - (:inf=>0.01)=>((:S,:I)=>(:I,:I)), - (:rec=>0.001)=>(:I=>:R), - (:death=>0.01)=>(:I=>:D)) - - expand_pres!(DiffusionSpace2D, SIRD) - - Diffusion = @decapode DiffusionSpace2D begin - (C, Ċ)::Form0{X} - - Ċ == ⋆₀⁻¹{X}(dual_d₁{X}(⋆₁{X}(k(d₀{X}(C))))) - end - - Superposition = @decapode DiffusionSpace2D begin - (C, Ċ₁, Ċ₂)::Form0{X} - - ∂ₜ{Form0{X}}(C) == Ċ₁ + Ċ₂ - end - - - sird_dec = pn2dec(DiffusionSpace2D, SIRD) - - compose_epi = @relation (S,I,R) begin - epi(S,I,R, D, Ṡ₁,İ₁, Ṙ₁, Ḋ₁) - diff(S, Ṡ₂) - diff(I, İ₂) - diff(R, Ṙ₂) - diff(D, Ḋ₂) - superpos(Ṡ₁,Ṡ₂, S) - superpos(İ₁,İ₂, I) - superpos(Ṙ₁,Ṙ₂, R) - superpos(Ḋ₁,Ḋ₂, D) - end - - - composed_epi = oapply(compose_epi, vcat( - [OpenDiagram(sird_dec, [:S, :I, :R, :D, :Ṡ, :İ, :Ṙ, :Ḋ])], - [OpenDiagram(Diffusion, [:C, :Ċ]) for i in 1:ns(SIRD)], - [OpenDiagram(Superposition, [:Ċ₁, :Ċ₂, :C]) for i in 1:ns(SIRD)] - )); - - res = diag2dwd(composed_epi.functor, in_vars = [:S, :I, :R, :D]) - - s = EmbeddedDeltaSet2D{Bool, Point{3,Float64}}() - points = [(0,0,0),(0,0,1),(0,1,0),(0,1,1)] - - add_vertices!(s, 4, point=points) - glue_sorted_triangle!(s, 1,2,3) - glue_sorted_triangle!(s, 2,3,4) - s[:edge_orientation] = false - s[:tri_orientation] = false - orient!(s) - sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3{Float64}}(s) - subdivide_duals!(sd, Barycenter()) - - #sd = dual(s); - - funcs = sym2func(sd) - funcs[:k] = Dict(:operator => 1.0 * I(ne(sd)), :type => MatrixFunc()) - merge!(funcs, gen_functions(SIRD, s)) - - Examples.contract_matrices!(res, funcs) - - sim, _ = gen_sim(res, funcs, sd); - - Random.seed!(42) - u = rand(nv(s)*4) - du = zeros(Float64,nv(s)*4) - - tempu = copy(u) - dt = 0.01 - for i in 1:1000 - sim(du, tempu, [],0) - tempu .+= du * dt - # Check that population is conserved - end - @test all(tempu .!= u) - - mass_diff = sum(vcat([⋆(Val{0}, sd)*(tempu[(1:nv(s)) .+ (i-1) * nv(s)]) for i in 1:4]...)) - - sum(vcat([⋆(Val{0}, sd)*(u[(1:nv(s)) .+ (i-1) * nv(s)]) for i in 1:4]...)) - @test mass_diff < 1e-6 -end diff --git a/test/Decapodes2/composition.jl b/test/composition.jl similarity index 100% rename from test/Decapodes2/composition.jl rename to test/composition.jl diff --git a/test/Decapodes2/coordinates.jl b/test/coordinates.jl similarity index 100% rename from test/Decapodes2/coordinates.jl rename to test/coordinates.jl diff --git a/test/Decapodes2/diag2dwd.jl b/test/diag2dwd.jl similarity index 100% rename from test/Decapodes2/diag2dwd.jl rename to test/diag2dwd.jl diff --git a/test/Decapodes2/meshes.jl b/test/meshes.jl similarity index 100% rename from test/Decapodes2/meshes.jl rename to test/meshes.jl diff --git a/test/Decapodes2/multiscalearrays.jl b/test/multiscalearrays.jl similarity index 100% rename from test/Decapodes2/multiscalearrays.jl rename to test/multiscalearrays.jl diff --git a/test/Decapodes2/rewrite.jl b/test/rewrite.jl similarity index 100% rename from test/Decapodes2/rewrite.jl rename to test/rewrite.jl diff --git a/test/runtests.jl b/test/runtests.jl index 78da514..ae6ff98 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,17 +1,38 @@ using Test -@testset "Decapodes" begin - include("Decapodes.jl") +@testset "Construction" begin + include("diag2dwd.jl") end -@testset "PetriNets" begin - include("PetriNets.jl") +@testset "SummationDecapode Construction" begin + include("summation.jl") end -@testset "Examples" begin - include("Examples.jl") +@testset "Composition" begin + include("composition.jl") end -@testset verbose = true "Decapodes2" begin - include(joinpath(@__DIR__, "Decapodes2", "runtests.jl")) +@testset "Coordinates" begin + include("coordinates.jl") end + +@testset "Mesh Loading" begin + include("meshes.jl") +end + +@testset "MultiScaleArrays.jl Integration" begin + include("multiscalearrays.jl") +end + +@testset "Average Rewriting" begin + include("rewrite.jl") +end + +@testset "Simulation" begin + include("simulation.jl") +end + +@testset "Visualization" begin + include("visualization.jl") +end + diff --git a/test/Decapodes2/simulation.jl b/test/simulation.jl similarity index 100% rename from test/Decapodes2/simulation.jl rename to test/simulation.jl diff --git a/test/Decapodes2/summation.jl b/test/summation.jl similarity index 100% rename from test/Decapodes2/summation.jl rename to test/summation.jl diff --git a/test/Decapodes2/visualization.jl b/test/visualization.jl similarity index 100% rename from test/Decapodes2/visualization.jl rename to test/visualization.jl