Reconfiguring Decapodes Analysis for new CCL Interop

packages/algjulia-interop/Project.toml

@@ -5,6 +5,15 @@ authors = ["CatColab team"]
 version = "0.1.1"
 
 [deps]
+HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
+JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
+Oxygen = "df9a0d86-3283-4920-82dc-4555fc0d1d8b"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
+TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[weakdeps]
 ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
@@ -14,40 +23,21 @@ Decapodes = "679ab3ea-c928-4fe6-8d59-fd451142d391"
 DiagrammaticEquations = "6f00c28b-6bed-4403-80fa-30e0dc12f317"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
-JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-Preferences = "21216c6a-2e73-6563-6e65-726566657250"
-REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
-[weakdeps]
-PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
-
 [extensions]
-SysImageExt = "PackageCompiler"
+CatlabExt = ["Catlab", "ACSets"]
+DecapodesExt = ["DiagrammaticEquations", "Decapodes", "ACSets", "CombinatorialSpaces", "JSON3", "StaticArrays", "LinearAlgebra", "ComponentArrays", "Distributions", "CoordRefSystems", "GeometryBasics", "OrdinaryDiffEq"]
 
 [compat]
-ACSets = "0.2.21"
-Catlab = "0.16.20"
-CombinatorialSpaces = "0.7.4"
-ComponentArrays = "0.15"
-CoordRefSystems = "0.18.9"
-Decapodes = "0.6"
-DiagrammaticEquations = "0.2"
-Distributions = "0.25"
-GeometryBasics = "0.5.7"
-IJulia = "1.26.0"
-JSON3 = "1"
-LinearAlgebra = "1"
-MLStyle = "0.4"
-OrdinaryDiffEq = "6.101.0"
-PackageCompiler = "2.2.1"
-Preferences = "1.5.0"
-REPL = "1.11.0"
+Catlab = "0.17.3"
+HTTP = "1.10.19"
+JSON3 = "1.14.3"
+MLStyle = "0.4.17"
+Oxygen = "1.7.5"
 Reexport = "1.2.2"
-StaticArrays = "1"
+StructTypes = "1.11.0"
+TOML = "1.0.3"
 julia = "1.11"

packages/algjulia-interop/ext/CatlabExt.jl

@@ -0,0 +1,86 @@
+module CatlabExt
+
+using ACSets
+using Catlab: Presentation, FreeSchema
+import Catlab: id, dom
+using Catlab.CategoricalAlgebra.Pointwise.FunctorialDataMigrations.Yoneda: 
+  yoneda, colimit_representables, DiagramData
+using CatColabInterop, Oxygen, HTTP
+import CatColabInterop: endpoint
+
+# workaround, could be upstreamed
+""" 
+If we call `id(X::FreeSchema.Ob{:generator})` we get a FreeSchema.Hom. However,
+there is no analogue for attribute types. For `list_to_hom` (a function within)
+`colimit_representables` to be concisely defined, we need such an analogue. 
+`IdAttr` is workaround type intended to fill this absence.
+"""
+struct IdAttr{T} x::FreeSchema.AttrType{:generator} end
+dom(i::IdAttr) = i.x
+id(x::FreeSchema.AttrType{:generator}) = IdAttr{:id}(x)
+
+
+function model_to_schema(m::Model)::Tuple{Schema, Dict{String,Symbol}}
+  obs, homs, attrtypes, attrs = Symbol[],[],[],[]
+  names = Dict{String, Symbol}()
+  for stmt in m.obGenerators
+    names[stmt.id] = Symbol(only(stmt.label))
+    if stmt.obType.content == "Entity"
+        push!(obs, names[stmt.id])
+    elseif stmt.obType.content == "AttrType"
+        push!(attrtypes, names[stmt.id])
+    else 
+      error(stmt.obType)
+    end
+  end
+  for stmt in m.morGenerators
+    h = (Symbol(only(stmt.label)), names[stmt.dom.content], names[stmt.cod.content])
+    names[stmt.id] = h[1]
+    if stmt.morType.content == "Attr"
+      push!(attrs, h)
+    else 
+      push!(homs, h)
+    end
+  end
+  (Schema(Presentation(BasicSchema{Symbol}(obs, homs, attrtypes, attrs, []))), names)
+end
+
+function diagram_to_data(d::Types.Diagram, names::Dict{String,Symbol})::DiagramData
+  data = DiagramData()
+  for o in d.obGenerators
+    names[o.id] = Symbol(only(o.label))
+    push!(data.reprs[names[o.over.content]], names[o.id])
+  end
+  for m in d.morGenerators
+    p1 = names[m.cod.content] => Symbol[]
+    p2 = names[m.dom.content] => [names[m.over.content]]
+    push!(data.eqs, p1 => p2)
+  end
+  data
+end
+
+""" Receiver of the data already knows the schema """
+function acset_to_json(X::ACSet, S::Schema)::AbstractDict
+  Dict{Symbol, Union{Int,Vector{Int}}}(
+    [t => nparts(X, t) for t in types(S)] 
+    ∪ [f => X[f] for f in homs(S; just_names=true)] 
+    ∪ [f => getvalue.(X[f]) for f in attrs(S; just_names=true)] )
+end
+
+""" 
+Compute ACSet colimit of a diagrammatic instance. Return a tabular representation 
+"""
+function endpoint(::Val{:ACSetColim})
+  @post "/acsetcolim" function(req::HTTP.Request)
+    payload = json(req, ModelDiagram)
+    schema, names = model_to_schema(payload.model)
+    data = diagram_to_data(payload.diagram, names)
+    acset_type = AnonACSet(
+      schema; type_assignment=Dict(a=>Nothing for a in schema.attrtypes))
+    y = yoneda(constructor(acset_type))
+    res = colimit_representables(data, y)
+    acset_to_json(res, schema)
+  end
+end
+
+end # module

packages/algjulia-interop/ext/DecapodesExt/DecapodesExt.jl

@@ -0,0 +1,65 @@
+module DecapodesExt 
+
+import Base: run
+
+# algebraicjulia dependencies
+using ACSets
+using DiagrammaticEquations
+using Decapodes
+using CombinatorialSpaces
+
+# dependencies
+import JSON3
+using StaticArrays
+using MLStyle
+using LinearAlgebra
+using ComponentArrays
+using Distributions # for initial conditions
+
+# meshing
+using CoordRefSystems
+using GeometryBasics: Point2, Point3
+Point3D = Point3{Float64}
+
+# simulation
+using OrdinaryDiffEq
+
+
+using CatColabInterop, Oxygen, HTTP
+import CatColabInterop: Model, ModelDiagram, ObGenerator, MorGenerator, DiagramObGenerator, DiagramMorGenerator
+import CatColabInterop: endpoint
+
+# necessary to export
+export infer_types!, evalsim, default_dec_generate, default_dec_matrix_generate, DiagonalHodge, ComponentArray
+
+include("interop.jl")
+
+# functions for geometry and initial conditions
+include("geometry.jl")
+
+# helper functions for Navier-Stokes
+include("ns_helper.jl")
+
+# constructing initial conditions
+include("initial_conditions.jl")
+
+# constructs a simulation from the payload
+include("simulation.jl")
+
+# executes the analysis 
+include("execute.jl")
+
+"""
+"""
+function endpoint(::Val{:DecapodeSimulation})
+    @post "/decapodes-simulation" function(req::HTTP.Request)
+        payload = json(req, ModelDiagrma)
+        simulation = DecapodesSimulation(payload)
+        sim = evalsim(simulation.pode)
+        f = sim(simulation.geometry.dualmesh, simulation.generate, DiagonalHodge())
+        res = run(f, simulation, ComponentArray(k=0.5,))
+        sim_to_json(res)
+    end
+end
+
+end

packages/algjulia-interop/ext/DecapodesExt/execute.jl

@@ -0,0 +1,73 @@
+""" """
+function Base.run(fm, sim::DecapodeSimulation, constparam)
+    prob = ODEProblem(fm, sim.init, sim.duration, constparam)
+    soln = solve(prob, Tsit5(), saveat=0.01)
+    SimResult(soln, sim)
+end
+
+abstract type AbstractResult end
+
+struct SimResult <: AbstractResult
+    retcode::Enum{Int32} # ReturnCode
+    time::Vector{Float64}
+    state::Dict{String, Vector{AbstractArray{SVector{3, Float64}}}}
+    x::Vector{Float64} # axis
+    y::Vector{Float64}
+end
+export SimResult
+
+function SimResult(soln::ODESolution, system::DecapodeSimulation)
+    idx_bounds = indexing_bounds(system)
+    state_val_dict = variables_state(soln, system) # Dict("UUID1" => VectorMatrixSVectr...)
+    SimResult(soln.retcode, soln.t, state_val_dict, 0:idx_bounds.x, 0:idx_bounds.y)
+end
+# TODO generalize to HasDeltaSet
+
+points(system::DecapodeSimulation) = collect(values(system.geometry.dualmesh.subparts.point.m))
+indexing_bounds(system::DecapodeSimulation) = indexing_bounds(system.geometry.domain)
+
+""" for the variables in a system, associate them to their state values over the duration of the simulation """
+function variables_state(soln::ODESolution, system::DecapodeSimulation)
+    plottedVars = [ k for (k, v) in system.plotVariables if v == true ]
+    uuid2symb = Dict([ v => k for (k, v) in system.uuiddict]) # TODO why reverse again?
+    Dict([ String(uuid2symb[var]) => state_entire_sim(soln, system, uuid2symb[var]) for var ∈ plottedVars ])
+end
+
+""" given a simulation, a domain, and a variable, gets the state values over the duration of a simulation. 
+Called by `variables_state`[@ref] """
+function state_entire_sim(soln::ODESolution, system::DecapodeSimulation, var::Symbol)
+    map(1:length(soln.t)) do i
+        state_at_time(soln, system, var, i)
+    end
+end
+
+# TODO type `points`
+function state_at_time(soln::ODESolution, system::DecapodeSimulation, plotvar::Symbol, t::Int)
+    @match system.geometry.domain begin
+        # TODO check time indexing here
+        domain::Rectangle => state_at_time(soln, domain, plotvar, t) 
+        domain::Sphere => state_at_time(soln, domain, plotvar, t, points(system)) 
+        _ => throw(ImplError("state_at_time function for domain $domain"))
+    end
+end
+
+function state_at_time(soln::ODESolution, domain::Rectangle, var::Symbol, t::Int)
+    (x, y) = indexing_bounds(domain)
+    [SVector(i, j, getproperty(soln.u[t], var)[(x+1)*(i-1) + j]) for i in 1:x+1, j in 1:y+1]
+end
+
+# TODO just separated this from the SimResult function and added type parameters, but need to generalize
+function grid(pt3::Point3, grid_size::Vector{Int})
+    pt2 = [(pt3[1]+1)/2, (pt3[2]+1)/2]
+    [round(Int, pt2[1]*grid_size[1]), round(Int, pt2[2]*grid_size[2])]
+end
+
+function state_at_time(soln::ODESolution, domain::Sphere, var::Symbol, t::Int, points)
+    l , _ = indexing_bounds(domain) # TODO this is hardcoded to return 100, 100
+    northern_indices = filter(i -> points[i][3] > 0, keys(points)) 
+    map(northern_indices) do n
+        i, j = grid(points[n], [l, l]) # TODO
+        SVector(i, j, getproperty(soln.u[t], var)[n])
+    end
+end
+

packages/algjulia-interop/src/decapodes-service/geometry.jl → packages/algjulia-interop/ext/DecapodesExt/geometry.jl

@@ -83,6 +83,16 @@ function Geometry(json_object::AbstractDict)
     Geometry(domain)
 end
 
+function Geometry(analysis::Analysis)
+    domain = PREDEFINED_MESHES[Symbol(analysis.analysis[:mesh])]
+    Geometry(domain)
+end
+
+# function Geometry(payload::DiagramPayload{ThDecapode})
+#     domain = PREDEFINED_MESHES[Symbol(payload.data[:mesh])]
+#     Geometry(domain)
+# end
+
 # function Geometry(d::Domain, args...)
 #     throw(ImplError("The mesh ($(d)) is"))
 # end

packages/algjulia-interop/src/decapodes-service/analysis/initial_conditions.jl → packages/algjulia-interop/ext/DecapodesExt/initial_conditions.jl

@@ -80,7 +80,7 @@ end
 
 function initial_conditions(ics::GaussianIC, geometry::Geometry)
     c_dist = MvNormal(ics.ξ)
-    c = [pdf(c_dist, [p[1], p[2]]) for p ∈ geometry.dualmesh[:point]]
+    c = [Distributions.pdf(c_dist, [p[1], p[2]]) for p ∈ geometry.dualmesh[:point]]
     return c
 end