Llm/catlab0.15

Project.toml

@@ -4,11 +4,13 @@ authors = ["James Fairbanks", "Andrew Baas", "Evan Patterson", "Luke Morris", "G
 version = "0.2.1"
 
 [deps]
+ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
 AlgebraicRewriting = "725a01d3-f174-5bbd-84e1-b9417bad95d9"
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
@@ -21,12 +23,13 @@ MultiScaleArrays = "f9640e96-87f6-5992-9c3b-0743c6a49ffa"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Unicode = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
 [compat]
-AlgebraicRewriting = "0.1"
-Catlab = "0.14"
-CombinatorialSpaces = "0.4"
+AlgebraicRewriting = "0.2"
+Catlab = "0.15"
+CombinatorialSpaces = "0.5"
 DataStructures = "0.18.13"
 Distributions = "0.25"
 FileIO = "1.16"
@@ -39,7 +42,7 @@ MultiScaleArrays = "1.10"
 OrdinaryDiffEq = "6.47"
 PreallocationTools = "0.4"
 Requires = "1.3"
-julia = "1.8"
+julia = "1.9"
 
 [extras]
 AlgebraicPetri = "4f99eebe-17bf-4e98-b6a1-2c4f205a959b"

examples/climate/budyko_sellers.jl

@@ -0,0 +1,189 @@
+# AlgebraicJulia Dependencies
+using Catlab
+using Catlab.Graphics
+using CombinatorialSpaces
+using Decapodes
+
+# External Dependencies
+using MLStyle
+using MultiScaleArrays
+using LinearAlgebra
+using OrdinaryDiffEq
+using JLD2
+using GLMakie
+using GeometryBasics: Point2
+Point2D = Point2{Float64}
+
+####################
+# Define the model #
+####################
+
+# ϕ := Latitude
+# Tₛ(ϕ,t) := Surface temperature
+# Q(ϕ,t) := Insolation
+# C(ϕ) := Effective heat capacity
+# α(ϕ) := Albedo
+# A := Longwave emissions at 0°C
+# B := Increase in emissions per degree
+# D := Horizontal diffusivity
+budyko_sellers = @decapode begin
+    (Q,Tₛ)::Form0
+    (α,A,B,C,D,cosϕᵖ,cosϕᵈ)::Constant
+
+    Tₛ̇ == ∂ₜ(Tₛ)
+    ASR == (1 .- α) .* Q
+    OLR == A .+ (B .* Tₛ)
+    HT == (D ./ cosϕᵖ) .* ⋆(d(cosϕᵈ .* ⋆(d(Tₛ))))
+
+    Tₛ̇ == (ASR - OLR + HT) ./ C
+end
+
+# Infer the forms of dependent variables, and resolve which versions of DEC
+# operators to use.
+infer_types!(budyko_sellers, op1_inf_rules_1D, op2_inf_rules_1D)
+resolve_overloads!(budyko_sellers, op1_res_rules_1D, op2_res_rules_1D)
+
+# Visualize.
+to_graphviz(budyko_sellers)
+
+# Demonstrate storing as JSON.
+write_json_acset(budyko_sellers, "budyko_sellers.json")
+# When reading back in, we specify that all attributes are "Symbol"s.
+budyko_sellers2 = read_json_acset(SummationDecapode{Symbol,Symbol,Symbol}, "budyko_sellers.json")
+# Or, you could choose to interpret the data as "String"s.
+budyko_sellers3 = read_json_acset(SummationDecapode{String,String,String}, "budyko_sellers.json")
+
+###################
+# Define the mesh #
+###################
+
+s′ = EmbeddedDeltaSet1D{Bool, Point2D}()
+add_vertices!(s′, 30, point=Point2D.(range(-π/2 + π/32, π/2 - π/32, length=30), 0))
+add_edges!(s′, 1:nv(s′)-1, 2:nv(s′))
+orient!(s′)
+s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s′)
+subdivide_duals!(s, Circumcenter())
+
+########################################################
+# Define constants, parameters, and initial conditions #
+########################################################
+
+# This is a primal 0-form, with values at vertices.
+cosϕᵖ = map(x -> cos(x[1]), point(s′))
+# This is a dual 0-form, with values at edge centers.
+cosϕᵈ = map(edges(s′)) do e
+    (cos(point(s′, src(s′, e))[1]) + cos(point(s′, tgt(s′, e))[1])) / 2
+end
+
+α₀ = 0.354
+α₂ = 0.25
+α = map(point(s′)) do ϕ
+    α₀ + α₂*((1/2)*(3*ϕ[1]^2 - 1))
+end
+A = 210
+B = 2
+f = 0.70
+ρ = 1025
+cw = 4186
+H = 70
+C = map(point(s′)) do ϕ
+    f * ρ * cw * H
+end
+D = 0.6
+Q = map(point(s′)) do ϕ
+    450*cos(ϕ[1])
+end
+
+#Tₛ₀ = map(point(s′)) do ϕ
+#    12 .- 40*((1/2)*(3*(sin(ϕ[1]))^2 - 1))
+#end
+# Isothermal:
+Tₛ₀ = map(point(s′)) do ϕ
+    15
+end
+
+# Store these values to be passed to the solver.
+u₀ = construct(PhysicsState, [VectorForm(Q), VectorForm(Tₛ₀)], Float64[], [:Q, :Tₛ])
+constants_and_parameters = (
+    α = α,
+    A = A,
+    B = B,
+    C = C,
+    D = D,
+    cosϕᵖ = cosϕᵖ,
+    cosϕᵈ = cosϕᵈ)
+
+#############################################
+# Define how symbols map to Julia functions #
+#############################################
+
+hodge = GeometricHodge()
+function generate(sd, my_symbol; hodge=GeometricHodge())
+  op = @match my_symbol begin
+    :d₀ => x -> begin
+      d₀ = d(s,0)
+      d₀ * x
+    end
+    :dual_d₀ => x -> begin
+      dual_d₀ = dual_derivative(s,0)
+      dual_d₀ * x
+    end
+    :⋆₁ => x -> begin
+      ⋆₁ = ⋆(s,1)
+      ⋆₁ * x
+    end
+    :⋆₀⁻¹ => x -> begin
+      ⋆₀⁻¹ = inv_hodge_star(s,0)
+      ⋆₀⁻¹ * x
+    end
+    x => error("Unmatched operator $my_symbol")
+  end
+  return (args...) -> op(args...)
+end
+
+#######################
+# Generate simulation #
+#######################
+
+sim = eval(gensim(budyko_sellers, dimension=1))
+fₘ = sim(s, generate)
+
+##################
+# Run simulation #
+##################
+
+tₑ = 1e6
+
+# Julia will pre-compile the generated simulation the first time it is run.
+@info("Precompiling Solver")
+prob = ODEProblem(fₘ, u₀, (0, 1e-4), constants_and_parameters)
+soln = solve(prob, Tsit5())
+soln.retcode != :Unstable || error("Solver was not stable")
+
+# This next run should be fast.
+@info("Solving")
+prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters)
+soln = solve(prob, Tsit5())
+@info("Done")
+
+@save "budyko_sellers.jld2" soln
+
+#############
+# Visualize #
+#############
+
+lines(map(x -> x[1], point(s′)), findnode(soln(0.0), :Tₛ))
+lines(map(x -> x[1], point(s′)), findnode(soln(tₑ), :Tₛ))
+
+# Initial frame
+fig = Figure(resolution = (800, 800))
+ax1 = Axis(fig[1,1])
+xlims!(ax1, extrema(map(x -> x[1], point(s′))))
+ylims!(ax1, extrema(findnode(soln(tₑ), :Tₛ)))
+Label(fig[1,1,Top()], "Tₛ")
+
+# Animation
+frames = 100
+record(fig, "budyko_sellers.gif", range(0.0, tₑ; length=frames); framerate = 15) do t
+  lines!(fig[1,1], map(x -> x[1], point(s′)), findnode(soln(t), :Tₛ))
+end

src/Decapodes.jl

@@ -3,7 +3,6 @@ module Decapodes
 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

src/composition.jl

@@ -196,10 +196,10 @@ end
 # Infinite loop:
 oapply(r::RelationDiagram, podes::Vector{D}) where {D<:OpenSummationDecapode} =
   oapply_rename(r, podes)
-# oapply(r::RelationDiagram, podes::Vector{D}) where {D<:OpenSummationDecapode} =
-  # invoke(oapply,
-    # Tuple{UndirectedWiringDiagram, Vector{<:StructuredMulticospan{L}} where L},
-    # r, oapply_rename(r, podes))
+#oapply(r::RelationDiagram, podes::Vector{D}) where {D<:OpenSummationDecapode} =
+#  invoke(oapply,
+#    Tuple{UndirectedWiringDiagram, Vector{<:StructuredMulticospan{L}} where L},
+#    r, oapply_rename(r, podes))
 
 #oapply(r::RelationDiagram, pode::OpenSummationDecapode) = oapply(r, [pode])
 # Luke changed the above line to the below line, for e.g. the case: (Note H should be renamed to N.)

src/decapodeacset.jl

@@ -1,3 +1,6 @@
+using Catlab.Theories: FreeSchema
+using ACSets
+using ACSets.DenseACSets
 using DataStructures
 
 @present SchDecapode(FreeSchema) begin
@@ -290,9 +293,17 @@ op2_inf_rules_1D = [
   (proj1_type = :Literal, proj2_type = :Form1, res_type = :Form1, op_names = [:/, :./, :*, :.*]),
 
   (proj1_type = :Form0, proj2_type = :Literal, res_type = :Form0, op_names = [:/, :./, :*, :.*]),
-  (proj1_type = :Form1, proj2_type = :Literal, res_type = :Form1, op_names = [:/, :./, :*, :.*])]
-
+  (proj1_type = :Form1, proj2_type = :Literal, res_type = :Form1, op_names = [:/, :./, :*, :.*]),
+
+  (proj1_type = :Constant, proj2_type = :Form0, res_type = :Form0, op_names = [:/, :./, :*, :.*]),
+  (proj1_type = :Constant, proj2_type = :Form1, res_type = :Form1, op_names = [:/, :./, :*, :.*]),
+  (proj1_type = :Form0, proj2_type = :Constant, res_type = :Form0, op_names = [:/, :./, :*, :.*]),
+  (proj1_type = :Form1, proj2_type = :Constant, res_type = :Form1, op_names = [:/, :./, :*, :.*]),
 
+  (proj1_type = :Constant, proj2_type = :DualForm0, res_type = :DualForm0, op_names = [:/, :./, :*, :.*]),
+  (proj1_type = :Constant, proj2_type = :DualForm1, res_type = :DualForm1, op_names = [:/, :./, :*, :.*]),
+  (proj1_type = :DualForm0, proj2_type = :Constant, res_type = :DualForm0, op_names = [:/, :./, :*, :.*]),
+  (proj1_type = :DualForm1, proj2_type = :Constant, res_type = :DualForm1, op_names = [:/, :./, :*, :.*])]
 
 """
 These are the default rules used to do type inference in the 2D exterior calculus.
@@ -411,7 +422,7 @@ function infer_summands_and_summations!(d::SummationDecapode)
     all(t == :infer for t in types) && continue # We can  not infer
     inferred_type = types[findfirst(!=(:infer), types)]
     to_infer_idxs = filter(i -> d[:type][i] == :infer, idxs)
-    d[:type][to_infer_idxs] .= inferred_type
+    d[to_infer_idxs, :type] = inferred_type
     applied = true
   end
   return applied
@@ -563,7 +574,7 @@ function resolve_overloads!(d::SummationDecapode, op1_rules::Vector{NamedTuple{(
     src_type = d[:type][src]; tgt_type = d[:type][tgt]
     for rule in op1_rules
       if op1 == rule[:op] && src_type == rule[:src_type] && tgt_type == rule[:tgt_type]
-        d[:op1][op1_idx] = rule[:resolved_name]
+        d[op1_idx, :op1] = rule[:resolved_name]
         break
       end
     end
@@ -574,7 +585,7 @@ function resolve_overloads!(d::SummationDecapode, op1_rules::Vector{NamedTuple{(
     proj1_type = d[:type][proj1]; proj2_type = d[:type][proj2]; res_type = d[:type][res]
     for rule in op2_rules
       if op2 == rule[:op] && proj1_type == rule[:proj1_type] && proj2_type == rule[:proj2_type] && res_type == rule[:res_type]
-        d[:op2][op2_idx] = rule[:resolved_name]
+        d[op2_idx, :op2] = rule[:resolved_name]
         break
       end
     end

src/language.jl

@@ -232,7 +232,7 @@ function SummationDecapode(e::DecaExpr)
     recognize_types(d)
 
     fill_names!(d)
-    d[:name] .= normalize_unicode.(d[:name])
+    d[:name] = normalize_unicode.(d[:name])
     make_sum_mult_unique!(d)
     return d
 end

src/meshes.jl

@@ -1,6 +1,5 @@
 using Artifacts
 using Catlab
-using Catlab.CategoricalAlgebra
 using CombinatorialSpaces
 using FileIO
 using JSON
@@ -21,7 +20,7 @@ Icosphere(n) = Icosphere(n, 1.0)
 function loadmesh(s::Icosphere)
   1 <= s.n <= 5 || error("The only icosphere divisions supported are 1-5")
   m = loadmesh_helper("UnitIcosphere$(s.n).obj")
-  m[:point] .*= s.r
+  m[:point] = m[:point]*s.r
   return m
 end
 

src/rewrite.jl

@@ -1,4 +1,4 @@
-using Catlab, Catlab.Graphs, Catlab.CategoricalAlgebra
+using Catlab
 using AlgebraicRewriting
 using AlgebraicRewriting: rewrite