Llm/mem mgmt

Project.toml

@@ -14,6 +14,8 @@ 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"
@@ -26,6 +28,8 @@ 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"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Unicode = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
 [compat]

examples/Testing Range/AdvDiff.jl

@@ -0,0 +1,134 @@
+using Catlab
+using CombinatorialSpaces
+using CombinatorialSpaces.ExteriorCalculus
+using Decapodes
+using MultiScaleArrays
+using OrdinaryDiffEq
+using MLStyle
+using Distributions
+using LinearAlgebra
+using GeometryBasics: Point3
+
+Point3D = Point3{Float64}
+
+flatten(vfield::Function, mesh) =  ♭(mesh, DualVectorField(vfield.(mesh[triangle_center(mesh),:dual_point])))
+
+function generate(sd, my_symbol; hodge=GeometricHodge())
+    op = @match my_symbol begin
+      :k => x->2000x
+      _ => default_dec_generate(sd, my_symbol, hodge)
+    end
+
+    return (args...) ->  op(args...)
+  end
+
+begin
+    RADIUS = 6371+90
+    primal_earth = loadmesh(Icosphere(4, RADIUS))
+    # primal_earth = EmbeddedDeltaSet2D("Icosphere8.obj")
+    nploc = argmax(x -> x[3], primal_earth[:point])
+    primal_earth[:edge_orientation] .= false
+    orient!(primal_earth)
+    earth = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(primal_earth)
+    subdivide_duals!(earth, Circumcenter())
+end
+
+begin
+    AdvDiff = quote
+        C::Form0{X}
+        Ċ::Form0{X}
+        V::Form1{X}
+        ϕ::Form1{X}
+        ϕ₁::Form1{X}
+        ϕ₂::Form1{X}
+        # starC::DualForm2{X}
+        # lvc::Form1{X}
+        # Fick's first law
+        ϕ₁ ==  ∘(d₀,k,⋆₁)(C)
+        # Advective Flux
+        ϕ₂ == -(L₀(V, C))
+        # Superposition Principle
+        ϕ == plus(ϕ₁ , ϕ₂)
+        # Conservation of Mass
+        Ċ == ∘(dual_d₁,⋆₀⁻¹)(ϕ)
+        ∂ₜ(C) == Ċ
+    end
+
+    advdiff = parse_decapode(AdvDiff)
+    advdiffdp = SummationDecapode(advdiff)
+end
+
+#= function man_simulate(mesh, operators)
+    begin
+        #= d₀ = generate(mesh, :d₀)
+        k = generate(mesh, :k)
+        (⋆₁) = generate(mesh, :⋆₁)
+        (-) = generate(mesh, :-)
+        dual_d₁ = generate(mesh, :dual_d₁)
+        (⋆₀⁻¹) = generate(mesh, :⋆₀⁻¹)=#
+        L₀ = operators(mesh, :L₀)
+
+        hd1 = ⋆(1,mesh,hodge=GeometricHodge())
+        d0 = d(0,mesh)
+
+        dd1 = dual_derivative(1,mesh)
+        invhd0 = inv_hodge_star(0,mesh,hodge=GeometricHodge())
+
+        ϕ₁ = Vector{Float64}(undef, ne(mesh))
+        var"2" = Vector{Float64}(undef, ne(mesh))
+        var"3" = Vector{Float64}(undef, ne(mesh))
+        var"4" = Vector{Float64}(undef, nv(mesh))
+        Ċ = Vector{Float64}(undef, nv(mesh))
+        var"•1" = Vector{Float64}(undef, ne(mesh))
+        ϕ₂ = Vector{Float64}(undef, ne(mesh))
+        ϕ = Vector{Float64}(undef, ne(mesh))
+    end
+    return begin
+            f(du, u, p, t) = begin
+                    begin
+                        C = (findnode(u, :C)).values
+                        V = (findnode(u, :V)).values
+                    end
+                    println("--------------------")
+                    # ϕ₁ = (∘(⋆₁, k, d₀))(C)
+                    mul!(var"2", d0, C)
+                    var"3" .= 2000 .* var"2"
+                    mul!(ϕ₁, hd1, var"3")
+                    var"•1" .= L₀(V, C)
+                    ϕ₂ .= -var"•1"
+                    ϕ .= ϕ₁ .+ ϕ₂
+                    # Ċ = ((⋆₀⁻¹) ∘ dual_d₁)(ϕ)
+                    mul!(var"4", dd1, ϕ)
+                    mul!(Ċ, invhd0, var"4")
+                    # mul!(Ċ, mat2, ϕ)
+                    du .= 0.0
+                    begin
+                        (findnode(du, :C)).values .= Ċ
+                    end
+                end
+        end
+end =#
+
+sim = eval(gensim(advdiffdp, [:C, :V]))
+fₘ = sim(earth, generate)
+# fₘ = man_simulate(earth, generate)
+
+begin
+    c_dist = MvNormal(nploc[[1,2]], 100[1, 1])
+    c = [pdf(c_dist, [p[1], p[2]]./√RADIUS) for p in earth[:point]]
+
+    vmag = 500
+    velocity(p) = TangentBasis(CartesianPoint(p))((vmag/4, vmag/4))
+    v = flatten(velocity, earth)
+
+    u₀ = construct(PhysicsState, [VectorForm(c), VectorForm(collect(v))],Float64[], [:C, :V])
+    tₑ = 10
+    prob = ODEProblem(fₘ,u₀,(0, tₑ))
+    soln = solve(prob, Tsit5())
+end
+
+fig, ax, ob = GLMakie.mesh(primal_earth, color = findnode(soln(0), :C))
+for t in range(0.0, tₑ; length=300)
+    sleep(0.01)
+    ob.color = findnode(soln(t), :C)
+end

examples/Testing Range/Brusselator.jl

@@ -0,0 +1,209 @@
+using Catlab
+using CombinatorialSpaces
+using CombinatorialSpaces.ExteriorCalculus
+using Decapodes
+using MultiScaleArrays
+using OrdinaryDiffEq
+using MLStyle
+using Distributions
+using LinearAlgebra
+using GeometryBasics: Point3
+using TerminalLoggers
+using Logging
+using BenchmarkTools
+using GLMakie
+
+global_logger(TerminalLogger())
+
+Point3D = Point3{Float64}
+
+flatten(vfield::Function, mesh) =  ♭(mesh, DualVectorField(vfield.(mesh[triangle_center(mesh),:dual_point])))
+
+function generate(sd, my_symbol; hodge=GeometricHodge())
+    op = @match my_symbol begin
+        # :my_id => x -> id(x)
+        _ => default_dec_generate(sd, my_symbol, hodge)
+    end
+
+    return (args...) ->  op(args...)
+end
+
+begin
+    primal_earth = loadmesh(Icosphere(3))
+    # primal_earth = EmbeddedDeltaSet2D("Icosphere8.obj")
+    nploc = argmax(x -> x[3], primal_earth[:point])
+    primal_earth[:edge_orientation] .= false
+    orient!(primal_earth)
+    earth = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(primal_earth)
+    subdivide_duals!(earth, Circumcenter());
+end
+
+begin
+    Brusselator = SummationDecapode(parse_decapode(quote
+        (U, V)::Form0{X} 
+        (U2V, aTU)::Form0{X}
+        (U̇, V̇)::Form0{X}
+
+        (α, One)::Constant{X}
+        (F)::Parameter{X}
+
+        U2V == (U .* U) .* V
+        aTU == α * Δ(U)
+
+        U̇ == One + U2V - (4.4 * U) + aTU + F
+        V̇ == (3.4 * U) - U2V + aTU
+
+        ∂ₜ(U) == U̇
+        ∂ₜ(V) == V̇
+    end))
+
+    infer_types!(Brusselator)
+
+    resolve_overloads!(Brusselator)
+end
+
+#= function man_simulate(mesh, operators)
+    begin
+        # Δ₀ = generate(mesh, :Δ₀)
+        lpdr0   = Δ(0, mesh)
+        var"•2" = Vector{Float64}(undef, nv(mesh))
+        U2V     = Vector{Float64}(undef, nv(mesh))
+        var"•1" = Vector{Float64}(undef, nv(mesh))
+        aTU     = Vector{Float64}(undef, nv(mesh))
+        var"•4" = Vector{Float64}(undef, nv(mesh))
+        var"•6" = Vector{Float64}(undef, nv(mesh))
+        var"•5" = Vector{Float64}(undef, nv(mesh))
+        sum_1   = Vector{Float64}(undef, nv(mesh))
+        var"•3" = Vector{Float64}(undef, nv(mesh))
+        #U̇ = Vector{Float64}(undef, nv(mesh))
+        #V̇ = Vector{Float64}(undef, nv(mesh))
+    end
+    return begin
+        f(du, u, p, t) = begin
+            begin
+                U = (findnode(u, :U)).values
+                V = (findnode(u, :V)).values
+                One = p.One
+                var"4.4" = 4.4
+                var"3.4" = 3.4
+                α = p.α
+                F = p.F(t)
+            end
+            # println("--------------------")
+            # var"•2" = Δ₀(U)
+            mul!(var"•2", lpdr0, U)
+            #mul!(b, d₀, a)
+            var"•1" .= U .* U
+            U2V .= var"•1" .* V
+            aTU .= α .* var"•2"
+            var"•4" .= var"4.4" .* U
+            var"•6" .= var"3.4" .* U
+            var"•5" .= var"•6" .- U2V
+            sum_1 .= One .+ U2V
+            (findnode(du, :V)).values .= var"•5" .+ aTU
+            var"•3" .= sum_1 .- var"•4"
+            (findnode(du, :U)).values .= var"•3" .+ aTU .+ F
+        end
+    end
+end =#
+
+sim = eval(gensim(Brusselator))
+#= sim = eval_gensim(Brusselator, CUDA=true, no_infers=true,
+    verbose_errors=true, precision=Float64,
+    primals=[:J, :B], duals=[:E]) =#
+fₘ = sim(earth, generate)
+# fₘ = sim(earth)
+# fₘ = simulate(earth, generate)
+
+begin
+    U = map(earth[:point]) do (_,y,_)
+      abs(y)
+    end
+
+    V = map(earth[:point]) do (x,_,_)
+      abs(x)
+    end
+
+    One = ones(nv(earth))
+
+    F₁ = map(earth[:point]) do (_,_,z)
+      z ≥ 0.8 ? 5.0 : 0.0
+    end
+
+    F₂ = zeros(nv(earth))
+
+    constants_and_parameters = (
+        α = 0.001,
+        F = t -> t ≥ 1.1 ? F₁ : F₂,
+        One = One)
+
+    u₀ = construct(PhysicsState, [VectorForm(U), VectorForm(V)],Float64[], [:U, :V])
+    tₑ = 11.5
+    prob = ODEProblem(fₘ,u₀,(0, tₑ), constants_and_parameters)
+    soln = solve(prob, Tsit5())
+end
+
+fig, ax, ob = GLMakie.mesh(primal_earth, color = findnode(soln(0), :U))
+for t in range(0.0, tₑ; length=300)
+    sleep(0.01)
+    ob.color = findnode(soln(t), :U)
+end
+
+
+#UnicodePlots.scatterplot(soln.t[begin+1:end] - soln.t[begin:end-1])
+     ┌─────────────────────────────────────────────┐ 
+0.07 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡏⠂⠀⠀⠀⠀⠀⠀⠀⠀⢀⡄⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⢡⣄⣠⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⡤⠤⠖⠁⣦⣀⣀⣀⣀⣀⣀⣀⣀⡼⣁⣀⣀⡀⠀⠀⠀⠀⠀│ 
+     │⠚⠀⢉⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠚⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠀⠀│ 
+   0 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
+     └─────────────────────────────────────────────┘ 
+     ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀600⠀ 
+
+
+# GPU on Icosphere 5 Brusselator for 15.0 sim time.
+#     BenchmarkTools.Trial: 8 samples with 1 evaluation.
+# Range (min … max):  622.786 ms … 668.150 ms  ┊ GC (min … max): 0.00% … 0.00%
+# Time  (median):     633.515 ms               ┊ GC (median):    0.00%
+# Time  (mean ± σ):   639.551 ms ±  15.785 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%
+#
+#  █         ██ █ █   █                             █          █
+#  █▁▁▁▁▁▁▁▁▁██▁█▁█▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁█ ▁
+#  623 ms           Histogram: frequency by time          668 ms <
+#
+# Memory estimate: 22.76 MiB, allocs estimate: 391054.
+
+# GPU on Icosphere 5 Brusselator for 15.0 sim time, w/ @inbounds
+#BenchmarkTools.Trial: 8 samples with 1 evaluation.
+# Range (min … max):  603.889 ms … 655.395 ms  ┊ GC (min … max): 0.00% 
+#… 0.00%
+# Time  (median):     633.734 ms               ┊ GC (median):    0.00% 
+# Time  (mean ± σ):   633.116 ms ±  15.405 ms  ┊ GC (mean ± σ):  0.00% 
+#± 0.00%
+#
+#  █                     █       █  █  █       █  █            █       
+#  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█▁▁█▁▁█▁▁▁▁▁▁▁█▁▁█▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
+#  604 ms           Histogram: frequency by time          655 ms <     
+#
+# Memory estimate: 22.76 MiB, allocs estimate: 391115.
+
+# CPU on Icosphere 5 Brusselator for 15.0 sim time.
+#BenchmarkTools.Trial: 45 samples with 1 evaluation.
+# Range (min … max):   65.203 ms … 195.998 ms  ┊ GC (min … max):  0.00% … 47.45%
+# Time  (median):     108.710 ms               ┊ GC (median):     0.00% Time  (mean ± σ):   111.422 ms ±  30.407 ms  ┊ GC (mean ± σ):  12.80% ± 18.05%
+#
+#  ▄                  ▄▄█ 
+#  █▆▁▄▁▁▁▁▁▁▁▁▁▄▁▁▆█▆███▄▁▁▁▄▁▄▁▁▁▁▁▄▄▄▁▁▁▁▄▆▁▁▁▁▁▄▁▁▁▁▄▁▁▁▁▁▁▄ ▁     
+#  65.2 ms          Histogram: frequency by time          196 ms <     
+#
+# Memory estimate: 70.99 MiB, allocs estimate: 12102.