Centralizing Operators (#291)

.gitignore

@@ -1,5 +1,7 @@
 .DS_Store
 .ipynb_checkpoints
+*.ipynb
+*.mp4
 /Manifest.toml
 /dev/
 *.mem

docs/make.jl

@@ -15,7 +15,6 @@ if !(haskey(ENV, "GITHUB_ACTIONS") || haskey(ENV, "GITLAB_CI"))
   config["repo_root_url"] = "https://github.com/AlgebraicJulia/Decapodes.jl/blob/main/docs"
 end
 
-# const literate_dir = joinpath(@__DIR__, "..", "examples")
 const literate_dir = joinpath(@__DIR__, "literate")
 const generated_dir = joinpath(@__DIR__, "src", "examples")
 
@@ -36,8 +35,7 @@ end
 
 bib = CitationBibliography(
     joinpath(@__DIR__, "src", "decapodes_documenter.bib");
-    style=:numeric
-)
+    style=:numeric)
 
 @info "Building Documenter.jl docs"
 makedocs(
@@ -89,13 +87,11 @@ makedocs(
     "Canonical Models" => "canon.md",
     "Library Reference" => "api.md"
   ],
-  plugins=[bib]
-)
+  plugins=[bib])
 
 @info "Deploying docs"
 deploydocs(
   target = "build",
   repo   = "github.com/AlgebraicJulia/Decapodes.jl.git",
   branch = "gh-pages",
-  devbranch = "main"
-)
+  devbranch = "main")

docs/src/bsh/budyko_sellers_halfar.md

@@ -297,8 +297,7 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
         sum([nv*norm(nv)*x[e] for (e,nv) in zip(es,nvs)]) / sum(norm.(nvs))
       end
     end
-    :mag => x -> norm.(x)
-    x => error("Unmatched operator $my_symbol")
+    x => default_dec_generate(sd, my_symbol, hodge)
   end
   return (args...) -> op(args...)
 end

docs/src/cism/cism.md

@@ -57,7 +57,7 @@ halfar_eq2 = @decapode begin
   n::Constant
 
   ḣ == ∂ₜ(h)
-  ḣ == Γ * ∘(⋆, d, ⋆)(d(h) ∧₁₀ ((mag(♯(d(h)))^(n-1)) ∧₀₀ h^(n+2)))
+  ḣ == Γ * ∘(⋆, d, ⋆)(d(h) ∧₁₀ ((mag(♯ᵖᵖ(d(h)))^(n-1)) ∧₀₀ h^(n+2)))
 end
 
 to_graphviz(halfar_eq2)
@@ -192,29 +192,11 @@ constants_and_parameters = (
   stress_A = A)
 ```
 
-We provide here the mapping from symbols to differential operators. As more of the differential operators of the DEC are implemented, they are upstreamed to the Decapodes and CombinatorialSpaces libraries. Of course, users can also provide their own implementations of these operators and others as they see fit.
-
-## Setting up the Simulation
-```@example DEC
-function generate(sd, my_symbol; hodge=GeometricHodge())
-  # We pre-allocate matrices that encode differential operators.
-  op = @match my_symbol begin
-    :mag => x -> norm.(x)
-    :♯  => begin
-      sharp_mat = ♯_mat(sd, AltPPSharp())
-      x -> sharp_mat * x
-    end
-    x => error("Unmatched operator $my_symbol")
-  end
-  return (args...) -> op(args...)
-end
-```
-
 The `gensim` function takes our high-level representation of the physics equations and produces compiled simulation code. It performs optimizations such as allocating memory for intermediate variables, and so on.
 
 ```@example DEC
 sim = eval(gensim(ice_dynamics2D))
-fₘ = sim(sd, generate)
+fₘ = sim(sd, nothing)
 ```
 
 Julia is a "Just-In-Time" compiled language. That means that functions are compiled the first time they are called, and later calls to those functions skip this step. To get a feel for just how fast this simulation is, we will run the dynamics twice, once for a very short timespan to trigger pre-compilation, and then again for the actual dynamics.

docs/src/ebm_melt/ebm_melt.md

@@ -93,7 +93,7 @@ halfar_eq2 = @decapode begin
   Γ::Form1
   n::Constant
 
-  ∂ₜ(h)  == ∘(⋆, d, ⋆)(Γ * d(h) * avg₀₁(mag(♯(d(h)))^(n-1)) * avg₀₁(h^(n+2))) - melt
+  ∂ₜ(h)  == ∘(⋆, d, ⋆)(Γ * d(h) * avg₀₁(mag(♯ᵖᵖ(d(h)))^(n-1)) * avg₀₁(h^(n+2))) - melt
 end
 
 glens_law = @decapode begin
@@ -274,21 +274,6 @@ constants_and_parameters = (
   halfar_stress_ρ = halfar_ρ,
   halfar_stress_g = g,
   melting_Dₕ₂ₒ = Dₕ₂ₒ)
-
-# Define how symbols map to Julia functions
-function generate(sd, my_symbol; hodge=GeometricHodge())
-  op = @match my_symbol begin
-    :♯ => begin
-      sharp_mat = ♯_mat(sd, AltPPSharp())
-      x -> sharp_mat * x
-    end
-    :mag => x -> begin
-      norm.(x)
-    end
-    x => error("Unmatched operator $my_symbol")
-  end
-  return (args...) -> op(args...)
-end
 ```
 
 ## Generate and run simulation 
@@ -299,7 +284,7 @@ We will save the final ice thickness data in a .jld2 file, an HDF5-compatible fi
 
 ``` @example DEC
 sim = eval(gensim(budyko_sellers_halfar_water))
-fₘ = sim(s, generate)
+fₘ = sim(s, nothing, DiagonalHodge())
 
 tₑ = 100.0
 

docs/src/grigoriev/grigoriev.md

@@ -145,23 +145,11 @@ ice_dynamics = apex(ice_dynamics_cospan)
 to_graphviz(ice_dynamics)
 ```
 
-## Define our functions
-
-``` @example DEC
-function generate(sd, my_symbol; hodge=GeometricHodge())
-  op = @match my_symbol begin
-    :mag => x -> norm.(x)
-    x => error("Unmatched operator $my_symbol")
-  end
-  return op
-end
-```
-
 ## Generate simulation
 
 ``` @example DEC
 sim = eval(gensim(ice_dynamics, dimension=2))
-fₘ = sim(sd, generate)
+fₘ = sim(sd, nothing)
 ```
 
 ## Run

docs/src/halmo/halmo.md

@@ -145,7 +145,6 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
   op = @match my_symbol begin
     # This is a new function.
     :σ => sigmoid
-    :mag => x -> norm.(x)
     # Remaining operations (such as our differential operators) are built-in.
     _ => error("Unmatched operator $my_symbol")
   end

docs/src/ice_dynamics/ice_dynamics.md

@@ -53,7 +53,7 @@ halfar_eq2 = @decapode begin
   n::Constant
 
   ḣ == ∂ₜ(h)
-  ḣ == Γ * ∘(⋆, d, ⋆)(d(h) ∧₁₀ ((mag(♯(d(h)))^(n-1)) ∧₀₀ h^(n+2)))
+  ḣ == Γ * ∘(⋆, d, ⋆)(d(h) ∧₁₀ ((mag(♯ᵖᵖ(d(h)))^(n-1)) ∧₀₀ h^(n+2)))
 end
 
 to_graphviz(halfar_eq2)
@@ -165,7 +165,7 @@ In order to solve our equations, we will need numerical linear operators that gi
 ``` @example DEC
 function generate(sd, my_symbol; hodge=GeometricHodge())
   op = @match my_symbol begin
-    :♯ => x -> begin
+    :♯ᵖᵖ => x -> begin
       # This is an implementation of the "sharp" operator from the exterior
       # calculus, which takes co-vector fields to vector fields.
       # This could be up-streamed to the CombinatorialSpaces.jl library. (i.e.
@@ -183,8 +183,7 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
         sum([nv*norm(nv)*x[e] for (e,nv) in zip(es,nvs)]) / sum(norm.(nvs))
       end
     end
-    :mag => x -> norm.(x)
-    x => error("Unmatched operator $my_symbol")
+    x => default_dec_generate(sd, my_symbol, hodge)
   end
   return (args...) -> op(args...)
 end
@@ -331,27 +330,11 @@ constants_and_parameters = (
   stress_A = A)
 ```
 
-## Define our functions
-
-``` @example DEC
-function generate(sd, my_symbol; hodge=GeometricHodge())
-  op = @match my_symbol begin
-    :♯ => begin
-      sharp_mat = ♯_mat(sd, AltPPSharp())
-      x -> sharp_mat * x
-    end
-    :mag => x -> norm.(x)
-    x => error("Unmatched operator $my_symbol")
-  end
-  return (args...) -> op(args...)
-end
-```
-
 ## Generate simulation
 
 ``` @example DEC
 sim = eval(gensim(ice_dynamics2D, dimension=2))
-fₘ = sim(sd, generate)
+fₘ = sim(sd, nothing)
 ```
 
 ## Pre-compile and run 2D
@@ -447,7 +430,7 @@ constants_and_parameters = (
 
 ``` @example DEC
 sim = eval(gensim(ice_dynamics2D, dimension=2))
-fₘ = sim(sd, generate)
+fₘ = sim(sd, nothing)
 ```
 
 For brevity's sake, we'll skip the pre-compilation cell.

src/operators.jl

@@ -19,9 +19,11 @@ function default_dec_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::Discret
 
   op = @match my_symbol begin
 
-    # :plus => (+)
+    # Misc.
+    :plus => (+)
     :(-) || :neg => x -> -1 .* x
     :ln => (x -> log.(x))
+    :mag => x -> norm.(x)
 
     # Musical Isomorphisms
     :♯ᵖᵖ => dec_♯_p(sd)
@@ -99,6 +101,8 @@ function default_dec_matrix_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::
     :∧₂₀ => dec_pair_wedge_product(Tuple{2,0}, sd)
     :∧₁₁ => dec_pair_wedge_product(Tuple{1,1}, sd)
 
+    :♭♯ => ♭♯_mat(sd)
+
     # Averaging Operator
     :avg₀₁ => dec_avg₀₁(sd)
 

src/simulation.jl

@@ -601,7 +601,7 @@ const NONMATRIX_OPTIMIZABLE_DEC_OPERATORS = Set([:⋆₁⁻¹, :∧₀₁, :∧
 
 const NON_OPTIMIZABLE_CPU_OPERATORS = Set([:♯ᵖᵖ, :♯ᵈᵈ, :♭ᵈᵖ,
                                            :∧ᵖᵈ₁₁, :∧ᵖᵈ₀₁, :∧ᵈᵖ₁₁, :∧ᵈᵖ₁₀, :∧ᵈᵈ₁₁, :∧ᵈᵈ₁₀, :∧ᵈᵈ₀₁,
-                                           :ι₁₁, :ι₁₂, :ℒ₁, :Δᵈ₀ , :Δᵈ₁, :Δ₀⁻¹, :neg])
+                                           :ι₁₁, :ι₁₂, :ℒ₁, :Δᵈ₀ , :Δᵈ₁, :Δ₀⁻¹, :neg, :mag])
 const NON_OPTIMIZABLE_CUDA_OPERATORS = Set{Symbol}()