Roe refactor

roe/Project.toml

@@ -0,0 +1,46 @@
+name = "DEC"
+uuid = "e670b126-1168-4583-bd0f-29deacb39f6a"
+authors = ["Owen Lynch <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+AssociatedTests = "e00e7eca-deca-4415-9dc4-9b3efe792c16"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
+ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
+Crayons = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
+Decapodes = "679ab3ea-c928-4fe6-8d59-fd451142d391"
+Dtries = "fb203528-72e9-47b6-b8ee-6a26b9f77273"
+GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
+Metatheory = "e9d8d322-4543-424a-9be4-0cc815abe26c"
+Moshi = "2e0e35c7-a2e4-4343-998d-7ef72827ed2d"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+StructEquality = "6ec83bb0-ed9f-11e9-3b4c-2b04cb4e219c"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
+TermInterface = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+AbstractTrees = "0.4.5"
+AssociatedTests = "0.1.0"
+CairoMakie = "0.12.5"
+Colors = "0.12.11"
+CombinatorialSpaces = "0.6.7"
+Crayons = "4.1.1"
+Decapodes = "0.5.5"
+GeometryBasics = "0.4.11"
+MLStyle = "0.4.17"
+Moshi = "0.3.1"
+OrderedCollections = "1.6.3"
+OrdinaryDiffEq = "6.86.0"
+Reexport = "1.2.2"
+StructEquality = "2.1.0"
+SymbolicUtils = "1.5.1"
+Test = "1.11.0"

roe/README.md

@@ -0,0 +1,133 @@
+# Roe
+
+This is a refactor of the core of diagrammatic equations, attempting to achieve a "more Julionic" approach to the problem of typed computer algebra via direct use of multiple dispatch.
+
+## Signature
+
+The "signature" of the DEC is encoded in a module `ThDEC`, in the following way.
+
+First, we make a type `Sort`, elements of which represent types in the discrete exterior calculus, for instance scalars, dual/primal forms of any degree, and vector fields.
+
+Then, we make a Julia function for each sort in the DEC. We define these Julia functions to *act on sorts*. So for instance, for the wedge product, we write a function definition like
+
+```julia
+function ∧(s1::Sort, s2::Sort)
+  @match (s1, s2) begin
+    (Form(i, isdual), Scalar()) || (Scalar(), Form(i, isdual)) => Form(i, isdual)
+    (Form(i1, isdual), Form(i2, isdual)) =>
+      if i1 + i2 <= 2
+        Form(i1 + i2, isdual)
+      else
+        throw(SortError("Can only take a wedge product when the dimensions of the forms add to less than 2: tried to wedge product $i1 and $i2"))
+      end
+    _ => throw(SortError("Can only take a wedge product of two forms of the same duality"))
+  end
+end
+```
+
+The advantage of encoding the signature in this way is twofold.
+
+1. We can give high-quality, context-specific errors when types fail to match.
+2. It doesn't depend on any external libraries (except for MLStyle for convenience); it is just Plain Old Julia.
+
+## Using the signature to wrap symbolic algebra
+
+We can then wrap various symbolic frameworks by including the sorts as type parameters/metadata. For instance, for Metatheory, we create a wrapper struct `Var{s::Sort}` which wraps a `Metatheory.Id` (note that `s::Sort` rather than `s<:Sort`) We then define our methods on this struct as
+
+```julia
+unop_dec = [:∂ₜ, :d, :★, :-, :♯, :♭]
+for unop in unop_dec
+  @eval begin
+    @nospecialize
+    function $unop(v::Var{s}) where s
+      s′ = $unop(s)
+      Var{s′}(roe(v), addcall!(graph(v), $unop, (id(v),)))
+    end
+
+    export $unop
+  end
+end
+```
+
+SymbolicUtils.jl is totally analogous:
+
+```julia
+unop_dec = [:∂ₜ, :d, :★, :♯, :♭, :-]
+for unop in unop_dec
+  @eval begin
+    @nospecialize
+    function $unop(v::SymbolicUtils.BasicSymbolic{DECVar{s}}) where s
+      s′ = $unop(s)
+      SymbolicUtils.Term{DECVar{s′}}($unop, [v])
+    end
+
+    export $unop
+  end
+end
+```
+
+SymbolicUtils gets confused when the type parameter to `BasicSymbolic` is not a type: we work around this by passing in `DECVar{s}` (name subject to change), which is a zero-field struct that wraps a `Sort` as a type parameter.
+
+## Models and namespacing
+
+Models are then plain old Julia functions that accept as their first argument a "roe." The point of the "Roe" is to record variables that are created and equations that are asserted by the function. It looks something like:
+
+```julia
+struct Roe{T}
+  vars::Dtry{T}
+  eqs::Vector{Tuple{T, T}}
+end
+```
+
+The type parameter `T` could be instantiated with `BasicSymbolic{<:DECVar}` or `Var`, depending on whether we are working with SymbolicUtils or Metatheory.
+
+So, for instance, the Klausmeier model might look like
+
+```julia
+function klausmeier(roe::Roe)
+  @vars roe n::DualForm0 w::DualForm0 dX::Form1 a::Constant{DualForm0} ν::Constant{DualForm0}
+  @vars roe m::Number
+  # The equations for our model
+  @eq roe (∂ₜ(w) == a + w + (w * (n^2)) + ν * L(dX,w))
+  @eq roe (∂ₜ(n) == w * n^2 - m*n + Δ(n))
+end
+```
+
+Namespacing is achieved by moving the roe into a namespace before passing it into submodels. So, for instance, to make a model with two Klausmeier submodels that share the same `m`, we could do:
+
+```julia
+function double_klausmeier(roe::Roe)
+  klausmeier(namespaced(roe, :k1))
+  klausmeier(namespaced(roe, :k2))
+
+  @eq roe (roe.k1.m == roe.k2.m)
+end
+```
+
+The implementation of `namespace` would be something like
+
+```julia
+function namespace(roe::Roe{T}, name::Symbol) where {T}
+  Roe(get(roe.vars, name, Dtry{T}()), roe.eqs)
+end
+```
+
+Here, `get` either gets a pre-existing subnamespace at `name`, or creates a new subnamespace and inserts it at `name`. An alternative implementation would have a `namespace::Vec{Symbol}` field on `Roe` which is used to prefix every newly created variable.
+
+An alternative to adding the equation `roe.k1.m == roe.k2.m` would be to have `klausmeier` take `m` as a parameter, which would look like
+
+```julia
+function klausmeier(roe::Roe, m)
+  @vars roe n::DualForm0 w::DualForm0 dX::Form1 a::Constant{DualForm0} ν::Constant{DualForm0}
+  # The equations for our model
+  @eq roe (∂ₜ(w) == a + w + (w * (n^2)) + ν * L(dX,w))
+  @eq roe (∂ₜ(n) == w * n^2 - m*n + Δ(n))
+end
+
+function double_klausemeier(roe::Roe)
+  @vars roe m::Number
+
+  klausmeier(namespaced(roe, :k1), m)
+  klausmeier(namespaced(roe, :k2), m)
+end
+```

roe/docs/literate/egraphs.jl

@@ -0,0 +1,22 @@
+# This lesson covers the internals of Metatheory.jl-style E-Graphs. Let's reuse the heat_equation model on a new roe.
+roe = Roe(DEC.ThDEC.Sort);
+function heat_equation(roe)
+  u = fresh!(roe, PrimalForm(0), :u)
+
+  ∂ₜ(u) ≐ Δ(u)
+
+  ([u], [])
+end
+
+# We apply the model to the roe and collect its state variables.
+(state, _) = heat_equation(roe)
+
+# Recall from the Introduction that an E-Graph is a bipartite graph of ENodes and EClasses. Let's look at the EClasses: 
+classes = roe.graph.classes
+# The keys are Metatheory Id types which store an Id. The values are EClasses, which are implementations of equivalence classes. Nodes which share the same EClass are considered equivalent.
+
+
+
+# The constants in Roe are a dictionary of hashes of functions and constants. Let's extract just the values again: 
+vals = collect(values(e.graph.constants))
+# The `u` is ::RootVar{Form} and ∂ₜ, ★, d are all functions defined in ThDEC/signature.jl file. 

roe/docs/literate/heatequation.jl

@@ -0,0 +1,74 @@
+# Load AlgebraicJulia dependencies
+using DEC
+import DEC.ThDEC: Δ # conflicts with CombinatorialSpaces
+
+# load other dependencies
+using ComponentArrays
+using CombinatorialSpaces 
+using GeometryBasics
+using OrdinaryDiffEq
+Point2D = Point2{Float64}
+Point3D = Point3{Float64}
+using CairoMakie
+
+## Here we define the 1D heat equation model with one state variable and no parameters. That is, given an e-graph "roe," we define `u` to be a primal 0-form. The root variable carries a reference to the e-graph which it resides in. We then assert that the time derivative of the state is just its Laplacian. We return the state variable. 
+function heat_equation(roe)
+    u = fresh!(roe, PrimalForm(0), :u)
+
+    ∂ₜ(u) ≐ Δ(u)
+
+    ([u], [])
+end
+
+# Since this is a model in the DEC, we need to initialize the primal and dual meshes.
+rect = triangulated_grid(100, 100, 1, 1, Point3D);
+d_rect = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(rect);
+subdivide_duals!(d_rect, Circumcenter());
+
+# Now that we have a dual mesh, we can associate operators in our theory with precomputed matrices from Decapodes.jl.
+op_lookup = ThDEC.precompute_matrices(d_rect, DiagonalHodge())
+
+# Now we produce a "vector field" function which, given a model and operators in a theory, returns a function to be passed to the ODESolver. In stages, this function 
+#
+#   1) extracts the Root Variables (state or parameter term) and runs the extractor along the e-graph, 
+#   2) extracts the derivative terms from the model into an SSA
+#   3) yields a function accepting derivative terms, state terms, and parameter terms, whose body is both the lines, and derivatives. 
+vf = vfield(heat_equation, op_lookup)
+
+# Let's initialize the 
+U = first.(d_rect[:point]);
+
+# TODO component arrays
+constants_and_parameters = ()
+
+# We will run this for 500 timesteps.
+t0 = 500.0
+
+@info("Precompiling Solver")
+prob = ODEProblem(vf, U, (0, t0), constants_and_parameters);
+soln = solve(prob, Tsit5());
+
+## 1-D HEAT EQUATION WITH DIFFUSIVITY
+
+function heat_equation_with_constants(roe)
+  u = fresh!(roe, PrimalForm(0), :u)
+  k = fresh!(roe, Scalar(), :k)
+  ℓ = fresh!(roe, Scalar(), :ℓ)
+
+  ∂ₜ(u) ≐ k * Δ(u)
+
+  ([u], [k])
+end
+
+# we can reuse the mesh and operator lookup
+vf = vfield(heat_equation_with_constants, operator_lookup)
+
+# we can reuse the initial condition U but are specifying diffusivity constants. 
+constants_and_parameters = ComponentArray(k=0.25,);
+t0 = 500
+
+@info("Precompiling solver")
+prob = ODEProblem(vf, U, (0, t0), constants_and_parameters);
+soln = solve(prob, Tsit5());
+
+

roe/docs/literate/tutorial.jl

@@ -0,0 +1,52 @@
+# This tutorial is a slower-paced introduction into the design. Here, we will construct a simple exponential model.
+using DEC
+using Test
+using Metatheory.EGraphs
+using ComponentArrays
+using GeometryBasics
+using OrdinaryDiffEq
+Point2D = Point2{Float64}
+Point3D = Point3{Float64}
+
+using CairoMakie
+
+# We define our model of exponential growth. This model is a function which accepts a Roe and returns a tuple of State and Parameter variables. Let's break it down:
+#
+# 1. Function adds root variables (::RootVar) to the Roe. The root variables have no child nodes.
+# 2. Our model makes claims about what terms equal one another. The "≐" operator is an infix of "equate!" which claims unites the ids of the left and right VecExprs.
+# 3. The State and Parameter variables are returned. Each variable points to the same parent Roe.
+#
+#
+# Each variable points to the same Roe. 
+function exp_growth(roe)
+  u = fresh!(roe, PrimalForm(0), :u)
+  k = fresh!(roe, Scalar(), :k)
+
+  ∂ₜ(u) ≐ k * u
+
+  ([u], [k])
+end
+
+# We now need to initialize the primal and dual meshes we'll need to compute with.
+rect = triangulated_grid(100, 100, 1, 1, Point3D);
+d_rect = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(rect);
+subdivide_duals!(d_rect, Circumcenter());
+
+# For the theory of the DEC, we will need to associate each operator to the precomputed matrix specific to our dual mesh.
+operator_lookup = ThDEC.create_dynamic_model(d_rect, DiagonalHodge())
+
+# We now need to convert our model to an ODEProblem. In our case, ``vfield`` produces 
+vf = vfield(exp_growth, operator_lookup)
+
+U = first.(d_rect[:point]);
+
+constants_and_parameters = ComponentArray(k=-0.5,)
+
+t0 = 50.0
+
+@info("Precompiling Solver")
+prob = ODEProblem(vf, U, (0, t0), constants_and_parameters);
+soln = solve(prob, Tsit5());
+
+save_dynamics(soln, "decay.gif")
+

roe/docs/make.jl

@@ -0,0 +1,57 @@
+using Documenter
+using Literate
+using Distributed
+
+using DEC
+
+using CairoMakie
+
+# Set Literate.jl config if not being compiled on recognized service.
+# config = Dict{String,String}()
+# if !(haskey(ENV, "GITHUB_ACTIONS") || haskey(ENV, "GITLAB_CI"))
+#   config["nbviewer_root_url"] = "https://nbviewer.jupyter.org/github/AlgebraicJulia/DEC.jl/blob/gh-pages/dev"
+#   config["repo_root_url"] = "https://github.com/AlgebraicJulia/Decapodes.jl/blob/main/docs"
+# end
+
+const literate_dir = joinpath(@__DIR__, "..", "examples")
+const generated_dir = joinpath(@__DIR__, "src", "examples")
+
+@info "Building literate files"
+for (root, dirs, files) in walkdir(literate_dir)
+  out_dir = joinpath(generated_dir, relpath(root, literate_dir))
+  pmap(files) do file
+    f,l = splitext(file)
+    if l == ".jl" && !startswith(f, "_")
+      Literate.markdown(joinpath(root, file), out_dir;
+        config=config, documenter=true, credit=false)
+      Literate.notebook(joinpath(root, file), out_dir;
+        execute=true, documenter=true, credit=false)
+    end
+  end
+end
+@info "Completed literate"
+
+pages = Any[]
+push!(pages, "DEC.jl"      => "index.md")
+push!(pages, "Library Reference" => "api.md")
+
+@info "Building Documenter.jl docs"
+makedocs(
+  modules   = [DEC],
+  format    = Documenter.HTML(
+    assets = ["assets/analytics.js"],
+  ),
+  remotes   = nothing,
+  sitename  = "DEC.jl",
+  doctest   = false,
+  checkdocs = :none,
+  pages     = pages)
+
+
+@info "Deploying docs"
+deploydocs(
+  target = "build",
+  repo   = "github.com/AlgebraicJulia/DEC.jl.git",
+  branch = "gh-pages",
+  devbranch = "main"
+)