Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Decapodes conversion to SciML fails for Halfar #141

Open
djinnome opened this issue Nov 28, 2023 — with Slack · 50 comments
Open

Decapodes conversion to SciML fails for Halfar #141

djinnome opened this issue Nov 28, 2023 — with Slack · 50 comments
Assignees

Comments

Copy link

<@U03T6HR1328> <@U03JKH9FNBE> <@U03TU2CFAQY> am now running into:
ERROR: ArgumentError: cannot reinterpret an `ForwardDiff.Dual{nothing, Float64, 11}` array to `ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}` whose first dimension has size `99`. The resulting array would have non-integral first dimension.
More details in thread. Am still investigating.

Slack Message

@djinnome
Copy link
Author

minimal" reproducing example (not minimal because I don't have a sufficient sense of decapode goings-on to easily reduce further).

import Pkg

Pkg.add(name="Decapodes", rev="main")
Pkg.add(name="ACSets", rev="main")

Pkg.add([
  "Catlab",
  "CombinatorialSpaces",
  "MLStyle",
  "MultiScaleArrays",
  "LinearAlgebra",
  "OrdinaryDiffEq",
  "JLD2",
  "SparseArrays",
  "Statistics",
  "GLMakie",
  "GeometryBasics",
  "ComponentArrays"])

# AlgebraicJulia Dependencies
using Catlab
using Catlab.Graphics
using CombinatorialSpaces
using Decapodes
using ComponentArrays

# External Dependencies
using MLStyle
using MultiScaleArrays
using LinearAlgebra
using OrdinaryDiffEq
using JLD2
using SparseArrays
using Statistics
using GLMakie # Just for visualization
using GeometryBasics: Point2, Point3
Point2D = Point2{Float64};
Point3D = Point3{Float64};

@info("Packages Loaded")

halfar_eq2 = @decapode begin
  h::Form0
  Γ::Form1
  n::Constant== ∂ₜ(h)
  ḣ == (, d, )(Γ * d(h) * avg₀₁(mag((d(h)))^(n-1)) * avg₀₁(h^(n+2)))
end
glens_law = @decapode begin
  Γ::Form1
  (A,ρ,g,n)::Constant

  Γ == (2/(n+2))*A**g)^n
end

@info("Decapodes Defined")

ice_dynamics_composition_diagram = @relation () begin
  dynamics(Γ,n)
  stress(Γ,n)
end
ice_dynamics_cospan = oapply(ice_dynamics_composition_diagram,
  [Open(halfar_eq2, [,:n]),
  Open(glens_law, [,:n])])
ice_dynamics = apex(ice_dynamics_cospan)
ice_dynamics1D = expand_operators(ice_dynamics)
infer_types!(ice_dynamics1D, op1_inf_rules_1D, op2_inf_rules_1D)
resolve_overloads!(ice_dynamics1D, op1_res_rules_1D, op2_res_rules_1D)

s_prime = EmbeddedDeltaSet1D{Bool, Point2D}()
add_vertices!(s_prime, 100, point=Point2D.(range(-2, 2, length=100), 0))
add_edges!(s_prime, 1:nv(s_prime)-1, 2:nv(s_prime))
orient!(s_prime)
s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s_prime)
subdivide_duals!(s, Circumcenter())

@info("Spaces Defined")

function generate(sd, my_symbol; hodge=GeometricHodge())
  e_vecs = map(edges(sd)) do e
    point(sd, sd[e, :∂v0]) - point(sd, sd[e, :∂v1])
  end
  neighbors = map(vertices(sd)) do v
    union(incident(sd, v, :∂v0), incident(sd, v, :∂v1))
  end
  n_vecs = map(neighbors) do es
    [e_vecs[e] for e in es]
  end
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for e in 1:ne(s)
      append!(J, [s[e,:∂v0],s[e,:∂v1]])
      append!(I, [e,e])
      append!(V, [0.5, 0.5])
  end
  avg_mat = sparse(I,J,V)
  op = @match my_symbol begin
    :♯ => x -> begin
      map(neighbors, n_vecs) do es, nvs
        sum([nv*norm(nv)*x[e] for (e,nv) in zip(es,nvs)]) / sum(norm.(nvs))
      end
    end
    :mag => x -> begin
      norm.(x)
    end
    :avg₀₁ => x -> begin
      avg_mat * x
    end
    :^ => (x,y) -> x .^ y
    :* => (x,y) -> x .* y
    :show => x -> begin
      @show x
      x
    end
    x => error("Unmatched operator $my_symbol")
  end
  return (args...) -> op(args...)
end

sim = eval(gensim(ice_dynamics1D, dimension=1))
fₘ = sim(s, generate)

function f(flow_rate, ice_density, u_init_arr)
  n = 3
  ρ = ice_density
  g = 9.8101
  A = fill(flow_rate, ne(s))
  tₑ = 5e13

  u₀ = ComponentArray(dynamics_h = u_init_arr)
  
  constants_and_parameters = (
    n = n,
    stress_ρ = ρ,
    stress_g = g,
    stress_A = A)
  prob = ODEProblem{true, SciMLBase.FullSpecialize}(fₘ, u₀, (0, tₑ), constants_and_parameters)
  @info("Solving")
  soln = solve(prob, Tsit5())
  @info("Done")

  return soln(tₑ)
end


h₀ = map(x -> try sqrt(1. - x[1]^2) catch DomainError return 0.0 end, point(s_prime))

y = f(1e-16, 910., h₀)

Pkg.add("ForwardDiff")
using ForwardDiff: gradient;

gradient(x -> f(1e-16, 910., x), h₀)

results in this stacktrace:

ERROR: ArgumentError: cannot reinterpret an `ForwardDiff.Dual{nothing, Float64, 11}` array to `ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}` whose first dimension has size `99`.
The resulting array would have non-integral first dimension.
Stacktrace:
  [1] (::Base.var"#thrownonint#326")(S::Type, T::Type, dim::Int64)
    @ Base ./reinterpretarray.jl:53
  [2] reinterpret
    @ ./reinterpretarray.jl:71 [inlined]
  [3] get_tmp(dc::PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, u::ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}})
    @ PreallocationTools ~/.julia/packages/PreallocationTools/nhCNl/src/PreallocationTools.jl:51
  [4] (::var"#f#34"{PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, var"#20#33"{var"#17#30"}, var"#20#33"{var"#16#29"{SparseMatrixCSC{Float64, Int64}}}, var"#20#33"{var"#15#28"}, var"#20#33"{var"#12#25"{Vector{Vector{Point2{Float64}}}, Vector{Vector{Int64}}}}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}})(du::ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, u::ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, p::NamedTuple{(:n, :stress_ρ, :stress_g, :stress_A), Tuple{Int64, Float64, Float64, Vector{Float64}}}, t::Float64)
    @ Main ~/.julia/packages/Decapodes/AQT8L/src/simulation.jl:432
  [5] ODEFunction
    @ ~/.julia/packages/SciMLBase/eK30d/src/scimlfunctions.jl:2407 [inlined]
  [6] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, Nothing, Float64, NamedTuple{(:n, :stress_ρ, :stress_g, :stress_A), Tuple{Int64, Float64, Float64, Vector{Float64}}}, Float64, Float64, Float64, Float64, Vector{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}}, ODESolution{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, 2, Vector{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}}}, ODEProblem{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, Tuple{Float64, Float64}, true, NamedTuple{(:n, :stress_ρ, :stress_g, :stress_A), Tuple{Int64, Float64, Float64, Vector{Float64}}}, ODEFunction{true, SciMLBase.FullSpecialize, var"#f#34"{PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, var"#20#33"{var"#17#30"}, var"#20#33"{var"#16#29"{SparseMatrixCSC{Float64, Int64}}}, var"#20#33"{var"#15#28"}, var"#20#33"{var"#12#25"{Vector{Vector{Point2{Float64}}}, Vector{Vector{Int64}}}}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, var"#f#34"{PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, var"#20#33"{var"#17#30"}, var"#20#33"{var"#16#29"{SparseMatrixCSC{Float64, Int64}}}, var"#20#33"{var"#15#28"}, var"#20#33"{var"#12#25"{Vector{Vector{Point2{Float64}}}, Vector{Vector{Int64}}}}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}}, Vector{Float64}, Vector{Vector{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}}}, OrdinaryDiffEq.Tsit5Cache{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, SciMLBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, var"#f#34"{PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, var"#20#33"{var"#17#30"}, var"#20#33"{var"#16#29"{SparseMatrixCSC{Float64, Int64}}}, var"#20#33"{var"#15#28"}, var"#20#33"{var"#12#25"{Vector{Vector{Point2{Float64}}}, Vector{Vector{Int64}}}}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/JJd6g/src/perform_step/low_order_rk_perform_step.jl:792
  [7] __init(prob::ODEProblem{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, Tuple{Float64, Float64}, true, NamedTuple{(:n, :stress_ρ, :stress_g, :stress_A), Tuple{Int64, Float64, Float64, Vector{Float64}}}, ODEFunction{true, SciMLBase.FullSpecialize, var"#f#34"{PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, var"#20#33"{var"#17#30"}, var"#20#33"{var"#16#29"{SparseMatrixCSC{Float64, Int64}}}, var"#20#33"{var"#15#28"}, var"#20#33"{var"#12#25"{Vector{Vector{Point2{Float64}}}, Vector{Vector{Int64}}}}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/JJd6g/src/solve.jl:502
  [8] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEq/JJd6g/src/solve.jl:10 [inlined]
  [9] __solve(::ODEProblem{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, Tuple{Float64, Float64}, true, NamedTuple{(:n, :stress_ρ, :stress_g, :stress_A), Tuple{Int64, Float64, Float64, Vector{Float64}}}, ODEFunction{true, SciMLBase.FullSpecialize, var"#f#34"{PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, var"#20#33"{var"#17#30"}, var"#20#33"{var"#16#29"{SparseMatrixCSC{Float64, Int64}}}, var"#20#33"{var"#15#28"}, var"#20#33"{var"#12#25"{Vector{Vector{Point2{Float64}}}, Vector{Vector{Int64}}}}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/JJd6g/src/solve.jl:5
 [10] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/JJd6g/src/solve.jl:1 [inlined]
 [11] #solve_call#34
    @ ~/.julia/packages/DiffEqBase/NpZ7U/src/solve.jl:557 [inlined]
 [12] solve_call
    @ ~/.julia/packages/DiffEqBase/NpZ7U/src/solve.jl:523 [inlined]
 [13] #solve_up#42
    @ ~/.julia/packages/DiffEqBase/NpZ7U/src/solve.jl:1006 [inlined]
 [14] solve_up
    @ ~/.julia/packages/DiffEqBase/NpZ7U/src/solve.jl:992 [inlined]
 [15] #solve#40
    @ ~/.julia/packages/DiffEqBase/NpZ7U/src/solve.jl:929 [inlined]
 [16] solve(prob::ODEProblem{ComponentVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, Tuple{Float64, Float64}, true, NamedTuple{(:n, :stress_ρ, :stress_g, :stress_A), Tuple{Int64, Float64, Float64, Vector{Float64}}}, ODEFunction{true, SciMLBase.FullSpecialize, var"#f#34"{PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 12}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, PreallocationTools.FixedSizeDiffCache{Vector{Float64}, Vector{ForwardDiff.Dual{nothing, Float64, 11}}}, var"#20#33"{var"#17#30"}, var"#20#33"{var"#16#29"{SparseMatrixCSC{Float64, Int64}}}, var"#20#33"{var"#15#28"}, var"#20#33"{var"#12#25"{Vector{Vector{Point2{Float64}}}, Vector{Vector{Int64}}}}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}, Diagonal{Float64, Vector{Float64}}, SparseMatrixCSC{Int64, Int64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/NpZ7U/src/solve.jl:919
 [17] f(flow_rate::Float64, ice_density::Float64, u_init_arr::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}})
    @ Main ~/GitRepo/causal_pyro/docs/source/pde_scratch/src/halfar_ice/ice_hackathon.jl:141
 [18] #39
    @ ~/GitRepo/causal_pyro/docs/source/pde_scratch/src/halfar_ice/ice_hackathon.jl:156 [inlined]
 [19] chunk_mode_gradient(f::var"#39#40", x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:123
 [20] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:0
 [21] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#39#40", Float64}, Float64, 12}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:17
 [22] gradient(f::Function, x::Vector{Float64})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:17
 [23] top-level scope
    @ ~/GitRepo/causal_pyro/docs/source/pde_scratch/src/halfar_ice/ice_hackathon.jl:156
Collapse
 This snippet was truncated for display; [see it in full](https://askemgroup.slack.com/files/U05A31F749Z/F066TS847NE/stacktrace.txt?origin_team=T03JAGZ2ZHU&origin_channel=C03MVAWRMPU)

@lukem12345
Copy link

I would try out Barycenter() instead of Circumcenter() as a quick check

@azane
Copy link

azane commented Nov 29, 2023

Barycenter gives the same error. :(

@djinnome
Copy link
Author

djinnome commented Dec 6, 2023

@ChrisRackauckas any progress to report on this issue?

@hagabbar
Copy link

Any updates on this issue?

@joshday
Copy link
Collaborator

joshday commented Dec 21, 2023

I'm finding someone to investigate this.

I don't suppose there's a more minimal reproducible example?

@LilithHafner
Copy link

Here's a slightly modified reproducer that throws an Internal error: encountered unexpected error in runtime and segfaults.

segfault_reproducer.jl
# AlgebraicJulia Dependencies
using Catlab
using CombinatorialSpaces
using Decapodes
using ComponentArrays

# External Dependencies
using MLStyle
using MultiScaleArrays
using LinearAlgebra
using OrdinaryDiffEq
using JLD2
using SparseArrays
using Statistics
using GeometryBasics: Point2, Point3
const Point2D = Point2{Float64};

@info("Packages Loaded")

halfar_eq2 = @decapode begin
  h::Form0
  Γ::Form1
  n::Constant== ∂ₜ(h)
  ḣ == (, d, )(Γ * d(h) * avg₀₁(mag((d(h)))^(n-1)) * avg₀₁(h^(n+2)))
end
glens_law = @decapode begin
  Γ::Form1
  (A,ρ,g,n)::Constant

  Γ == (2/(n+2))*A**g)^n
end

@info("Decapodes Defined")

ice_dynamics_composition_diagram = @relation () begin
  dynamics(Γ,n)
  stress(Γ,n)
end
ice_dynamics_cospan = oapply(ice_dynamics_composition_diagram,
  [Open(halfar_eq2, [,:n]),
  Open(glens_law, [,:n])])
ice_dynamics = apex(ice_dynamics_cospan)
ice_dynamics1D = expand_operators(ice_dynamics)
infer_types!(ice_dynamics1D, op1_inf_rules_1D, op2_inf_rules_1D)
resolve_overloads!(ice_dynamics1D, op1_res_rules_1D, op2_res_rules_1D)

s_prime = EmbeddedDeltaSet1D{Bool, Point2D}()
add_vertices!(s_prime, 100, point=Point2D.(range(-2, 2, length=100), 0))
add_edges!(s_prime, 1:nv(s_prime)-1, 2:nv(s_prime))
orient!(s_prime)
s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s_prime)
subdivide_duals!(s, Circumcenter())

@info("Spaces Defined")

function generate(sd, my_symbol; hodge=GeometricHodge())
  e_vecs = map(edges(sd)) do e
    point(sd, sd[e, :∂v0]) - point(sd, sd[e, :∂v1])
  end
  neighbors = map(vertices(sd)) do v
    union(incident(sd, v, :∂v0), incident(sd, v, :∂v1))
  end
  n_vecs = map(neighbors) do es
    [e_vecs[e] for e in es]
  end
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for e in 1:ne(s)
      append!(J, [s[e,:∂v0],s[e,:∂v1]])
      append!(I, [e,e])
      append!(V, [0.5, 0.5])
  end
  avg_mat = sparse(I,J,V)
  op = @match my_symbol begin
    :♯ => x -> begin
      map(neighbors, n_vecs) do es, nvs
        sum([nv*norm(nv)*x[e] for (e,nv) in zip(es,nvs)]) / sum(norm.(nvs))
      end
    end
    :mag => x -> begin
      norm.(x)
    end
    :avg₀₁ => x -> begin
      avg_mat * x
    end
    :^ => (x,y) -> x .^ y
    :* => (x,y) -> x .* y
    :show => x -> begin
      @show x
      x
    end
    x => error("Unmatched operator $my_symbol")
  end
  return (args...) -> op(args...)
end

sim = eval(gensim(ice_dynamics1D, dimension=1))
fₘ = sim(s, generate)

function f(u_init_arr)
  n = 3
  ρ = 1e-16
  g = 9.8101
  A = fill(910, ne(s))
  tₑ = 5e10

  u₀ = ComponentArray(dynamics_h = u_init_arr)

  constants_and_parameters = (
    n = n,
    stress_ρ = ρ,
    stress_g = g,
    stress_A = A)
  prob = ODEProblem{true, SciMLBase.FullSpecialize}(fₘ, u₀, (0, tₑ), constants_and_parameters)
  @info("Solving")
  soln = solve(prob, Tsit5())
  @info("Done")

  return soln(tₑ)
end


h₀ = map(x -> sqrt(max(0.0, 1.0 - x[1]^2)), point(s_prime))

y = f(h₀)

using ForwardDiff: gradient;

Invocation:

julia --project --startup-file=no -e 'using Revise; includet("segfault_reproducer.jl")'
Project.toml
[deps]
ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
Decapodes = "679ab3ea-c928-4fe6-8d59-fd451142d391"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
MultiScaleArrays = "f9640e96-87f6-5992-9c3b-0743c6a49ffa"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
Manifest.toml
# This file is machine-generated - editing it directly is not advised

julia_version = "1.9.4"
manifest_format = "2.0"
project_hash = "11c93f18f6ec7fc7b3bead5c17bad6976308b446"

[[deps.ACSets]]
deps = ["AlgebraicInterfaces", "Base64", "CompTime", "DataStructures", "JSON3", "MLStyle", "OrderedCollections", "Permutations", "Pkg", "PrettyTables", "Random", "Reexport", "SHA", "StaticArrays", "StructEquality", "Tables"]
git-tree-sha1 = "f6d27610156a7496c237c9375d663ad1453fecb3"
repo-rev = "main"
repo-url = "https://github.com/AlgebraicJulia/ACSets.jl.git"
uuid = "227ef7b5-1206-438b-ac65-934d6da304b8"
version = "0.2.11"

    [deps.ACSets.extensions]
    NautyACSetsExt = "nauty_jll"
    XLSXACSetsExt = "XLSX"

    [deps.ACSets.weakdeps]
    XLSX = "fdbf4ff8-1666-58a4-91e7-1b58723a45e0"
    nauty_jll = "55c6dc9b-343a-50ca-8ff2-b71adb3733d5"

[[deps.ADTypes]]
git-tree-sha1 = "332e5d7baeff8497b923b730b994fa480601efc7"
uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
version = "0.2.5"

[[deps.Adapt]]
deps = ["LinearAlgebra", "Requires"]
git-tree-sha1 = "cde29ddf7e5726c9fb511f340244ea3481267608"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.7.2"
weakdeps = ["StaticArrays"]

    [deps.Adapt.extensions]
    AdaptStaticArraysExt = "StaticArrays"

[[deps.AlgebraicInterfaces]]
git-tree-sha1 = "a81b76ea8d1801494562dd057315e4b7b25b8de1"
uuid = "23cfdc9f-0504-424a-be1f-4892b28e2f0c"
version = "0.1.1"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "bbec08a37f8722786d87bedf84eae19c020c4efa"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "7.7.0"

    [deps.ArrayInterface.extensions]
    ArrayInterfaceBandedMatricesExt = "BandedMatrices"
    ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
    ArrayInterfaceCUDAExt = "CUDA"
    ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
    ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
    ArrayInterfaceTrackerExt = "Tracker"

    [deps.ArrayInterface.weakdeps]
    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
    GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"

[[deps.ArrayLayouts]]
deps = ["FillArrays", "LinearAlgebra"]
git-tree-sha1 = "b08a4043e1c14096ef8efe4dd97e07de5cacf240"
uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
version = "1.4.5"
weakdeps = ["SparseArrays"]

    [deps.ArrayLayouts.extensions]
    ArrayLayoutsSparseArraysExt = "SparseArrays"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.BitTwiddlingConvenienceFunctions]]
deps = ["Static"]
git-tree-sha1 = "0c5f81f47bbbcf4aea7b2959135713459170798b"
uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b"
version = "0.1.5"

[[deps.CPUSummary]]
deps = ["CpuId", "IfElse", "PrecompileTools", "Static"]
git-tree-sha1 = "601f7e7b3d36f18790e2caf83a882d88e9b71ff1"
uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
version = "0.2.4"

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[deps.Catlab]]
deps = ["ACSets", "AlgebraicInterfaces", "Colors", "CompTime", "Compose", "DataStructures", "GATlab", "GeneralizedGenerated", "JSON3", "LightXML", "LinearAlgebra", "Logging", "MLStyle", "PrettyTables", "Random", "Reexport", "SparseArrays", "StaticArrays", "Statistics", "StructEquality", "Tables"]
git-tree-sha1 = "3286deb17c5799b1956f815b11487fd1aed97d3e"
uuid = "134e5e36-593f-5add-ad60-77f754baafbe"
version = "0.16.4"

    [deps.Catlab.extensions]
    CatlabConvexExt = "Convex"
    CatlabDataFramesExt = "DataFrames"
    CatlabGraphsExt = "Graphs"
    CatlabGraphvizExt = "Graphviz_jll"
    CatlabMetaGraphsExt = "MetaGraphs"
    CatlabSCSExt = "SCS"
    CatlabTikzPicturesExt = "TikzPictures"

    [deps.Catlab.weakdeps]
    Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
    DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
    Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
    Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85"
    MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
    SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
    TikzPictures = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra"]
git-tree-sha1 = "2118cb2765f8197b08e5958cdd17c165427425ee"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.19.0"
weakdeps = ["SparseArrays"]

    [deps.ChainRulesCore.extensions]
    ChainRulesCoreSparseArraysExt = "SparseArrays"

[[deps.CloseOpenIntervals]]
deps = ["Static", "StaticArrayInterface"]
git-tree-sha1 = "70232f82ffaab9dc52585e0dd043b5e0c6b714f1"
uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9"
version = "0.1.12"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.4"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.10"

[[deps.CombinatorialSpaces]]
deps = ["ACSets", "Catlab", "FileIO", "GeometryBasics", "LazyArrays", "LinearAlgebra", "MeshIO", "Reexport", "Requires", "SparseArrays", "StaticArrays"]
git-tree-sha1 = "0cfaa8852c89ea5a5b9f771f2be22ff06256e54f"
uuid = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
version = "0.6.0"

[[deps.Combinatorics]]
git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
version = "1.0.2"

[[deps.CommonSolve]]
git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c"
uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
version = "0.2.4"

[[deps.CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[deps.CompTime]]
deps = ["MLStyle", "MacroTools"]
git-tree-sha1 = "8c05059bc293a17f71cae4cd58b1fc18d4ede271"
uuid = "0fb5dd42-039a-4ca4-a1d7-89a96eae6d39"
version = "0.1.2"

[[deps.Compat]]
deps = ["UUIDs"]
git-tree-sha1 = "886826d76ea9e72b35fcd000e535588f7b60f21d"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "4.10.1"
weakdeps = ["Dates", "LinearAlgebra"]

    [deps.Compat.extensions]
    CompatLinearAlgebraExt = "LinearAlgebra"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "1.0.5+0"

[[deps.ComponentArrays]]
deps = ["ArrayInterface", "ChainRulesCore", "ForwardDiff", "Functors", "LinearAlgebra", "PackageExtensionCompat", "StaticArrayInterface", "StaticArraysCore"]
git-tree-sha1 = "6b4fd3b1c46e87c60f538fee25e8da1716870195"
uuid = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
version = "0.15.6"

    [deps.ComponentArrays.extensions]
    ComponentArraysAdaptExt = "Adapt"
    ComponentArraysConstructionBaseExt = "ConstructionBase"
    ComponentArraysGPUArraysExt = "GPUArrays"
    ComponentArraysRecursiveArrayToolsExt = "RecursiveArrayTools"
    ComponentArraysReverseDiffExt = "ReverseDiff"
    ComponentArraysSciMLBaseExt = "SciMLBase"
    ComponentArraysTrackerExt = "Tracker"
    ComponentArraysZygoteExt = "Zygote"

    [deps.ComponentArrays.weakdeps]
    Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
    ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
    GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
    RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
    SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
    Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[[deps.Compose]]
deps = ["Base64", "Colors", "DataStructures", "Dates", "IterTools", "JSON", "LinearAlgebra", "Measures", "Printf", "Random", "Requires", "Statistics", "UUIDs"]
git-tree-sha1 = "bf6570a34c850f99407b494757f5d7ad233a7257"
uuid = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
version = "0.9.5"

[[deps.ConcreteStructs]]
git-tree-sha1 = "f749037478283d372048690eb3b5f92a79432b34"
uuid = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
version = "0.2.3"

[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c53fc348ca4d40d7b371e71fd52251839080cbc9"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.5.4"

    [deps.ConstructionBase.extensions]
    ConstructionBaseIntervalSetsExt = "IntervalSets"
    ConstructionBaseStaticArraysExt = "StaticArrays"

    [deps.ConstructionBase.weakdeps]
    IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

[[deps.CpuId]]
deps = ["Markdown"]
git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406"
uuid = "adafc99b-e345-5852-983c-f28acb93d879"
version = "0.3.1"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.15.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.15"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Decapodes]]
deps = ["ACSets", "Artifacts", "Catlab", "CombinatorialSpaces", "ComponentArrays", "DataStructures", "FileIO", "GeometryBasics", "JLD2", "JSON", "LinearAlgebra", "MLStyle", "OrdinaryDiffEq", "PreallocationTools", "Requires", "Unicode"]
path = "/Users/x/.julia/dev/Decapodes"
uuid = "679ab3ea-c928-4fe6-8d59-fd451142d391"
version = "0.5.0"

[[deps.DiffEqBase]]
deps = ["ArrayInterface", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces"]
git-tree-sha1 = "83114bb5158ca74ac9dee3d01edf87cc0363b5db"
uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
version = "6.144.1"

    [deps.DiffEqBase.extensions]
    DiffEqBaseChainRulesCoreExt = "ChainRulesCore"
    DiffEqBaseDistributionsExt = "Distributions"
    DiffEqBaseEnzymeExt = ["ChainRulesCore", "Enzyme"]
    DiffEqBaseGeneralizedGeneratedExt = "GeneralizedGenerated"
    DiffEqBaseMPIExt = "MPI"
    DiffEqBaseMeasurementsExt = "Measurements"
    DiffEqBaseMonteCarloMeasurementsExt = "MonteCarloMeasurements"
    DiffEqBaseReverseDiffExt = "ReverseDiff"
    DiffEqBaseTrackerExt = "Tracker"
    DiffEqBaseUnitfulExt = "Unitful"

    [deps.DiffEqBase.weakdeps]
    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
    Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
    Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
    GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
    MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
    Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
    MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
    Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

[[deps.DiffEqNoiseProcess]]
deps = ["DiffEqBase", "Distributions", "GPUArraysCore", "LinearAlgebra", "Markdown", "Optim", "PoissonRandom", "QuadGK", "Random", "Random123", "RandomNumbers", "RecipesBase", "RecursiveArrayTools", "Requires", "ResettableStacks", "SciMLBase", "StaticArraysCore", "Statistics"]
git-tree-sha1 = "319377c927a4aa1f491228b2ac23f3554a3497c6"
uuid = "77a26b50-5914-5dd7-bc55-306e6241c503"
version = "5.20.0"

    [deps.DiffEqNoiseProcess.extensions]
    DiffEqNoiseProcessReverseDiffExt = "ReverseDiff"

    [deps.DiffEqNoiseProcess.weakdeps]
    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"

[[deps.DiffResults]]
deps = ["StaticArraysCore"]
git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.1.0"

[[deps.DiffRules]]
deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.15.1"

[[deps.Distances]]
deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
git-tree-sha1 = "66c4c81f259586e8f002eacebc177e1fb06363b0"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.11"
weakdeps = ["ChainRulesCore", "SparseArrays"]

    [deps.Distances.extensions]
    DistancesChainRulesCoreExt = "ChainRulesCore"
    DistancesSparseArraysExt = "SparseArrays"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
git-tree-sha1 = "9242eec9b7e2e14f9952e8ea1c7e31a50501d587"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.104"

    [deps.Distributions.extensions]
    DistributionsChainRulesCoreExt = "ChainRulesCore"
    DistributionsDensityInterfaceExt = "DensityInterface"
    DistributionsTestExt = "Test"

    [deps.Distributions.weakdeps]
    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
    DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
    Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.9.3"

[[deps.Downloads]]
deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
version = "1.6.0"

[[deps.DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.8"

[[deps.EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e3290f2d49e661fbd94046d7e3726ffcb2d41053"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.4+0"

[[deps.EnumX]]
git-tree-sha1 = "bdb1942cd4c45e3c678fd11569d5cccd80976237"
uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
version = "1.0.4"

[[deps.EnzymeCore]]
deps = ["Adapt"]
git-tree-sha1 = "2efe862de93cd87f620ad6ac9c9e3f83f1b2841b"
uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
version = "0.6.4"

[[deps.ExponentialUtilities]]
deps = ["Adapt", "ArrayInterface", "GPUArraysCore", "GenericSchur", "LinearAlgebra", "PrecompileTools", "Printf", "SparseArrays", "libblastrampoline_jll"]
git-tree-sha1 = "602e4585bcbd5a25bc06f514724593d13ff9e862"
uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18"
version = "1.25.0"

[[deps.ExprTools]]
git-tree-sha1 = "27415f162e6028e81c72b82ef756bf321213b6ec"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.10"

[[deps.Extents]]
git-tree-sha1 = "2140cd04483da90b2da7f99b2add0750504fc39c"
uuid = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
version = "0.1.2"

[[deps.FastBroadcast]]
deps = ["ArrayInterface", "LinearAlgebra", "Polyester", "Static", "StaticArrayInterface", "StrideArraysCore"]
git-tree-sha1 = "a6e756a880fc419c8b41592010aebe6a5ce09136"
uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
version = "0.2.8"

[[deps.FastClosures]]
git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef"
uuid = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
version = "0.3.2"

[[deps.FastLapackInterface]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "b12f05108e405dadcc2aff0008db7f831374e051"
uuid = "29a986be-02c6-4525-aec4-84b980013641"
version = "2.0.0"

[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "299dc33549f68299137e51e6d49a13b5b1da9673"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.16.1"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random"]
git-tree-sha1 = "5b93957f6dcd33fc343044af3d48c215be2562f1"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "1.9.3"
weakdeps = ["PDMats", "SparseArrays", "Statistics"]

    [deps.FillArrays.extensions]
    FillArraysPDMatsExt = "PDMats"
    FillArraysSparseArraysExt = "SparseArrays"
    FillArraysStatisticsExt = "Statistics"

[[deps.FiniteDiff]]
deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays"]
git-tree-sha1 = "c6e4a1fbe73b31a3dea94b1da449503b8830c306"
uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
version = "2.21.1"

    [deps.FiniteDiff.extensions]
    FiniteDiffBandedMatricesExt = "BandedMatrices"
    FiniteDiffBlockBandedMatricesExt = "BlockBandedMatrices"
    FiniteDiffStaticArraysExt = "StaticArrays"

    [deps.FiniteDiff.weakdeps]
    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
git-tree-sha1 = "cf0fe81336da9fb90944683b8c41984b08793dad"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.36"
weakdeps = ["StaticArrays"]

    [deps.ForwardDiff.extensions]
    ForwardDiffStaticArraysExt = "StaticArrays"

[[deps.FunctionWrappers]]
git-tree-sha1 = "d62485945ce5ae9c0c48f124a84998d755bae00e"
uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
version = "1.1.3"

[[deps.FunctionWrappersWrappers]]
deps = ["FunctionWrappers"]
git-tree-sha1 = "b104d487b34566608f8b4e1c39fb0b10aa279ff8"
uuid = "77dc65aa-8811-40c2-897b-53d922fa7daf"
version = "0.1.3"

[[deps.Functors]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "9a68d75d466ccc1218d0552a8e1631151c569545"
uuid = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
version = "0.4.5"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.GATlab]]
deps = ["AlgebraicInterfaces", "Colors", "Crayons", "DataStructures", "JSON", "MLStyle", "Random", "Reexport", "StructEquality", "UUIDs"]
git-tree-sha1 = "2ee8cc36b69ab432d42eec913fcdff60d1228dc0"
uuid = "f0ffcf3b-d13a-433e-917c-cc44ccf5ead2"
version = "0.1.0"

[[deps.GPUArraysCore]]
deps = ["Adapt"]
git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0"
uuid = "46192b85-c4d5-4398-a991-12ede77f4527"
version = "0.1.5"

[[deps.GeneralizedGenerated]]
deps = ["DataStructures", "JuliaVariables", "MLStyle", "Serialization"]
git-tree-sha1 = "60f1fa1696129205873c41763e7d0920ac7d6f1f"
uuid = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
version = "0.3.3"

[[deps.GenericSchur]]
deps = ["LinearAlgebra", "Printf"]
git-tree-sha1 = "fb69b2a645fa69ba5f474af09221b9308b160ce6"
uuid = "c145ed77-6b09-5dd9-b285-bf645a82121e"
version = "0.5.3"

[[deps.GeoInterface]]
deps = ["Extents"]
git-tree-sha1 = "d53480c0793b13341c40199190f92c611aa2e93c"
uuid = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
version = "1.3.2"

[[deps.GeometryBasics]]
deps = ["EarCut_jll", "Extents", "GeoInterface", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "424a5a6ce7c5d97cca7bcc4eac551b97294c54af"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.9"

[[deps.Graphs]]
deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "899050ace26649433ef1af25bc17a815b3db52b7"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.9.0"

[[deps.HostCPUFeatures]]
deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
git-tree-sha1 = "eb8fed28f4994600e29beef49744639d985a04b2"
uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0"
version = "0.1.16"

[[deps.HypergeometricFunctions]]
deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"]
git-tree-sha1 = "f218fe3736ddf977e0e772bc9a586b2383da2685"
uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
version = "0.3.23"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.Inflate]]
git-tree-sha1 = "ea8031dea4aff6bd41f1df8f2fdfb25b33626381"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.4"

[[deps.IntelOpenMP_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl"]
git-tree-sha1 = "31d6adb719886d4e32e38197aae466e98881320b"
uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
version = "2024.0.0+0"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.IrrationalConstants]]
git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.2.2"

[[deps.IterTools]]
git-tree-sha1 = "274c38bd733f9d29036d0a73658fff1dc1d3a065"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.9.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLD2]]
deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "PrecompileTools", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"]
git-tree-sha1 = "c2d0f45afcb5f6209155670bffd100c3b4937ea3"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.4.40"

[[deps.JLLWrappers]]
deps = ["Artifacts", "Preferences"]
git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.5.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.4"

[[deps.JSON3]]
deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"]
git-tree-sha1 = "eb3edce0ed4fa32f75a0a11217433c31d56bd48b"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.14.0"

    [deps.JSON3.extensions]
    JSON3ArrowExt = ["ArrowTypes"]

    [deps.JSON3.weakdeps]
    ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd"

[[deps.JuliaVariables]]
deps = ["MLStyle", "NameResolution"]
git-tree-sha1 = "49fb3cb53362ddadb4415e9b73926d6b40709e70"
uuid = "b14d175d-62b4-44ba-8fb7-3064adc8c3ec"
version = "0.2.4"

[[deps.JumpProcesses]]
deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "FunctionWrappers", "Graphs", "LinearAlgebra", "Markdown", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "StaticArrays", "UnPack"]
git-tree-sha1 = "c451feb97251965a9fe40bacd62551a72cc5902c"
uuid = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
version = "9.10.1"
weakdeps = ["FastBroadcast"]

    [deps.JumpProcesses.extensions]
    JumpProcessFastBroadcastExt = "FastBroadcast"

[[deps.KLU]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse_jll"]
git-tree-sha1 = "884c2968c2e8e7e6bf5956af88cb46aa745c854b"
uuid = "ef3ab10e-7fda-4108-b977-705223b18434"
version = "0.4.1"

[[deps.Krylov]]
deps = ["LinearAlgebra", "Printf", "SparseArrays"]
git-tree-sha1 = "8a6837ec02fe5fb3def1abc907bb802ef11a0729"
uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
version = "0.9.5"

[[deps.LaTeXStrings]]
git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.1"

[[deps.LayoutPointers]]
deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface"]
git-tree-sha1 = "62edfee3211981241b57ff1cedf4d74d79519277"
uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
version = "0.1.15"

[[deps.Lazy]]
deps = ["MacroTools"]
git-tree-sha1 = "1370f8202dac30758f3c345f9909b97f53d87d3f"
uuid = "50d2b5c4-7a5e-59d5-8109-a42b560f39c0"
version = "0.15.1"

[[deps.LazyArrays]]
deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "MatrixFactorizations", "SparseArrays"]
git-tree-sha1 = "9cfca23ab83b0dfac93cb1a1ef3331ab9fe596a5"
uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02"
version = "1.8.3"
weakdeps = ["StaticArrays"]

    [deps.LazyArrays.extensions]
    LazyArraysStaticArraysExt = "StaticArrays"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[deps.LevyArea]]
deps = ["LinearAlgebra", "Random", "SpecialFunctions"]
git-tree-sha1 = "56513a09b8e0ae6485f34401ea9e2f31357958ec"
uuid = "2d8b4e74-eb68-11e8-0fb9-d5eb67b50637"
version = "1.0.0"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
version = "0.6.4"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
version = "8.4.0+0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
version = "1.11.0+1"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl"]
git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.17.0+0"

[[deps.LightXML]]
deps = ["Libdl", "XML2_jll"]
git-tree-sha1 = "3a994404d3f6709610701c7dabfc03fed87a81f8"
uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
version = "0.9.1"

[[deps.LineSearches]]
deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
git-tree-sha1 = "7bbea35cec17305fc70a0e5b4641477dc0789d9d"
uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
version = "7.2.0"

[[deps.LinearAlgebra]]
deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LinearSolve]]
deps = ["ArrayInterface", "ConcreteStructs", "DocStringExtensions", "EnumX", "FastLapackInterface", "GPUArraysCore", "InteractiveUtils", "KLU", "Krylov", "Libdl", "LinearAlgebra", "MKL_jll", "PrecompileTools", "Preferences", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Sparspak", "StaticArraysCore", "UnPack"]
git-tree-sha1 = "b1148a6596c5cd9b2d848c26b500c79d102ffc5d"
uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
version = "2.21.1"

    [deps.LinearSolve.extensions]
    LinearSolveBandedMatricesExt = "BandedMatrices"
    LinearSolveBlockDiagonalsExt = "BlockDiagonals"
    LinearSolveCUDAExt = "CUDA"
    LinearSolveEnzymeExt = ["Enzyme", "EnzymeCore"]
    LinearSolveFastAlmostBandedMatricesExt = ["FastAlmostBandedMatrices"]
    LinearSolveHYPREExt = "HYPRE"
    LinearSolveIterativeSolversExt = "IterativeSolvers"
    LinearSolveKernelAbstractionsExt = "KernelAbstractions"
    LinearSolveKrylovKitExt = "KrylovKit"
    LinearSolveMetalExt = "Metal"
    LinearSolvePardisoExt = "Pardiso"
    LinearSolveRecursiveArrayToolsExt = "RecursiveArrayTools"

    [deps.LinearSolve.weakdeps]
    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
    BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
    Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
    EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
    FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
    HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771"
    IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
    KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
    KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
    Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
    Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
    RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"

[[deps.LogExpFunctions]]
deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.26"

    [deps.LogExpFunctions.extensions]
    LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
    LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables"
    LogExpFunctionsInverseFunctionsExt = "InverseFunctions"

    [deps.LogExpFunctions.weakdeps]
    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
    ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
    InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.LoopVectorization]]
deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"]
git-tree-sha1 = "0f5648fbae0d015e3abe5867bca2b362f67a5894"
uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
version = "0.12.166"
weakdeps = ["ChainRulesCore", "ForwardDiff", "SpecialFunctions"]

    [deps.LoopVectorization.extensions]
    ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"]
    SpecialFunctionsExt = "SpecialFunctions"

[[deps.MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl"]
git-tree-sha1 = "72dc3cf284559eb8f53aa593fe62cb33f83ed0c0"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2024.0.0+0"

[[deps.MLStyle]]
git-tree-sha1 = "bc38dff0548128765760c79eb7388a4b37fae2c8"
uuid = "d8e11817-5142-5d16-987a-aa16d5891078"
version = "0.4.17"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "b211c553c199c111d998ecdaf7623d1b89b69f93"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.12"

[[deps.ManualMemory]]
git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd"
uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667"
version = "0.1.8"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MatrixFactorizations]]
deps = ["ArrayLayouts", "LinearAlgebra", "Printf", "Random"]
git-tree-sha1 = "78f6e33434939b0ac9ba1df81e6d005ee85a7396"
uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
version = "2.1.0"

[[deps.MaybeInplace]]
deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "SparseArrays"]
git-tree-sha1 = "a85c6a98c9e5a2a7046bc1bb89f28a3241e1de4d"
uuid = "bb5d69b7-63fc-4a16-80bd-7e42200c7bdb"
version = "0.1.1"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
version = "2.28.2+0"

[[deps.Measures]]
git-tree-sha1 = "c13304c81eec1ed3af7fc20e75fb6b26092a1102"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.2"

[[deps.MeshIO]]
deps = ["ColorTypes", "FileIO", "GeometryBasics", "Printf"]
git-tree-sha1 = "8be09d84a2d597c7c0c34d7d604c039c9763e48c"
uuid = "7269a6da-0436-5bbc-96c2-40638cbb6118"
version = "0.4.10"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.1.0"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
version = "2022.10.11"

[[deps.MuladdMacro]]
git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab"
uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
version = "0.2.4"

[[deps.MultiScaleArrays]]
deps = ["DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "OrdinaryDiffEq", "Random", "RecursiveArrayTools", "SparseDiffTools", "Statistics", "StochasticDiffEq"]
git-tree-sha1 = "2c23fdd74aae9a7b5c4132e380e1ea01bf174473"
uuid = "f9640e96-87f6-5992-9c3b-0743c6a49ffa"
version = "1.12.0"

[[deps.NLSolversBase]]
deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
git-tree-sha1 = "a0b464d183da839699f4c79e7606d9d186ec172c"
uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
version = "7.8.3"

[[deps.NLsolve]]
deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"]
git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1"
uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
version = "4.5.1"

[[deps.NaNMath]]
deps = ["OpenLibm_jll"]
git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "1.0.2"

[[deps.NameResolution]]
deps = ["PrettyPrint"]
git-tree-sha1 = "1a0fa0e9613f46c9b8c11eee38ebb4f590013c5e"
uuid = "71a1bf82-56d0-4bbc-8a3c-48b961074391"
version = "0.1.5"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
version = "1.2.0"

[[deps.NonlinearSolve]]
deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "EnumX", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LazyArrays", "LineSearches", "LinearAlgebra", "LinearSolve", "MaybeInplace", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SimpleNonlinearSolve", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
git-tree-sha1 = "2ed2c035a8186344663d98337ea1b5c06a697591"
uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
version = "3.1.1"

    [deps.NonlinearSolve.extensions]
    NonlinearSolveBandedMatricesExt = "BandedMatrices"
    NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
    NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim"
    NonlinearSolveMINPACKExt = "MINPACK"
    NonlinearSolveNLsolveExt = "NLsolve"
    NonlinearSolveSymbolicsExt = "Symbolics"
    NonlinearSolveZygoteExt = "Zygote"

    [deps.NonlinearSolve.weakdeps]
    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
    FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce"
    LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
    MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
    NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
    Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
    Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[[deps.OffsetArrays]]
git-tree-sha1 = "6a731f2b5c03157418a20c12195eb4b74c8f8621"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.13.0"
weakdeps = ["Adapt"]

    [deps.OffsetArrays.extensions]
    OffsetArraysAdaptExt = "Adapt"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
version = "0.3.21+4"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
version = "0.8.1+0"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.Optim]]
deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
git-tree-sha1 = "01f85d9269b13fedc61e63cc72ee2213565f7a72"
uuid = "429524aa-4258-5aef-a3af-852621145aeb"
version = "1.7.8"

[[deps.OrderedCollections]]
git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.6.3"

[[deps.OrdinaryDiffEq]]
deps = ["ADTypes", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastBroadcast", "FastClosures", "FillArrays", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "IfElse", "InteractiveUtils", "LineSearches", "LinearAlgebra", "LinearSolve", "Logging", "LoopVectorization", "MacroTools", "MuladdMacro", "NonlinearSolve", "Polyester", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SimpleNonlinearSolve", "SimpleUnPack", "SparseArrays", "SparseDiffTools", "StaticArrayInterface", "StaticArrays", "TruncatedStacktraces"]
git-tree-sha1 = "a5a029dffca35d4afc28c75e8a074d49a7bc993e"
uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
version = "6.63.0"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "949347156c25054de2db3b166c52ac4728cbad65"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.31"

[[deps.PackageExtensionCompat]]
git-tree-sha1 = "fb28e33b8a95c4cee25ce296c817d89cc2e53518"
uuid = "65ce6f38-6b18-4e1d-a461-8949797d7930"
version = "1.0.2"
weakdeps = ["Requires", "TOML"]

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[deps.Parsers]]
deps = ["Dates", "PrecompileTools", "UUIDs"]
git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.8.1"

[[deps.Permutations]]
deps = ["Combinatorics", "LinearAlgebra", "Random"]
git-tree-sha1 = "c7745750b8a829bc6039b7f1f0981bcda526a946"
uuid = "2ae35dd2-176d-5d53-8349-f30d82d94d4f"
version = "0.4.19"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
version = "1.9.2"

[[deps.PoissonRandom]]
deps = ["Random"]
git-tree-sha1 = "a0f1159c33f846aa77c3f30ebbc69795e5327152"
uuid = "e409e4f3-bfea-5376-8464-e040bb5c01ab"
version = "0.4.4"

[[deps.Polyester]]
deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Requires", "Static", "StaticArrayInterface", "StrideArraysCore", "ThreadingUtilities"]
git-tree-sha1 = "fca25670784a1ae44546bcb17288218310af2778"
uuid = "f517fe37-dbe3-4b94-8317-1923a5111588"
version = "0.7.9"

[[deps.PolyesterWeave]]
deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"]
git-tree-sha1 = "240d7170f5ffdb285f9427b92333c3463bf65bf6"
uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad"
version = "0.2.1"

[[deps.PositiveFactorizations]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20"
uuid = "85a6dd25-e78a-55b7-8502-1745935b8125"
version = "0.2.4"

[[deps.PreallocationTools]]
deps = ["Adapt", "ArrayInterface", "ForwardDiff", "Requires"]
git-tree-sha1 = "01ac95fca7daabe77a9cb705862bd87016af9ddb"
uuid = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
version = "0.4.13"

    [deps.PreallocationTools.extensions]
    PreallocationToolsReverseDiffExt = "ReverseDiff"

    [deps.PreallocationTools.weakdeps]
    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"

[[deps.PrecompileTools]]
deps = ["Preferences"]
git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
version = "1.2.0"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.4.1"

[[deps.PrettyPrint]]
git-tree-sha1 = "632eb4abab3449ab30c5e1afaa874f0b98b586e4"
uuid = "8162dcfd-2161-5ef2-ae6c-7681170c5f98"
version = "0.2.0"

[[deps.PrettyTables]]
deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "Reexport", "StringManipulation", "Tables"]
git-tree-sha1 = "88b895d13d53b5577fd53379d913b9ab9ac82660"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "2.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "9ebcd48c498668c7fa0e97a9cae873fbee7bfee1"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.9.1"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Random123]]
deps = ["Random", "RandomNumbers"]
git-tree-sha1 = "552f30e847641591ba3f39fd1bed559b9deb0ef3"
uuid = "74087812-796a-5b5d-8853-05524746bad3"
version = "1.6.1"

[[deps.RandomNumbers]]
deps = ["Random", "Requires"]
git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111"
uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
version = "1.5.3"

[[deps.RecipesBase]]
deps = ["PrecompileTools"]
git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.3.4"

[[deps.RecursiveArrayTools]]
deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "Requires", "SparseArrays", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"]
git-tree-sha1 = "1afbe5e94fe6cb61500c5255db018ce010ef6d9f"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "3.2.1"

    [deps.RecursiveArrayTools.extensions]
    RecursiveArrayToolsMeasurementsExt = "Measurements"
    RecursiveArrayToolsMonteCarloMeasurementsExt = "MonteCarloMeasurements"
    RecursiveArrayToolsTrackerExt = "Tracker"
    RecursiveArrayToolsZygoteExt = "Zygote"

    [deps.RecursiveArrayTools.weakdeps]
    Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
    MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
    Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[[deps.RecursiveFactorization]]
deps = ["LinearAlgebra", "LoopVectorization", "Polyester", "PrecompileTools", "StrideArraysCore", "TriangularSolve"]
git-tree-sha1 = "8bc86c78c7d8e2a5fe559e3721c0f9c9e303b2ed"
uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
version = "0.2.21"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.ResettableStacks]]
deps = ["StaticArrays"]
git-tree-sha1 = "256eeeec186fa7f26f2801732774ccf277f05db9"
uuid = "ae5879a3-cd67-5da8-be7f-38c6eb64a37b"
version = "1.1.1"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "f65dcb5fa46aee0cf9ed6274ccbd597adc49aa7b"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.1"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6ed52fdd3382cf21947b15e8870ac0ddbff736da"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.4.0+0"

[[deps.RuntimeGeneratedFunctions]]
deps = ["ExprTools", "SHA", "Serialization"]
git-tree-sha1 = "6aacc5eefe8415f47b3e34214c1d79d2674a0ba2"
uuid = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
version = "0.5.12"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
version = "0.7.0"

[[deps.SIMDTypes]]
git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c"
uuid = "94e857df-77ce-4151-89e5-788b33177be4"
version = "0.1.0"

[[deps.SLEEFPirates]]
deps = ["IfElse", "Static", "VectorizationBase"]
git-tree-sha1 = "3aac6d68c5e57449f5b9b865c9ba50ac2970c4cf"
uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
version = "0.6.42"

[[deps.SciMLBase]]
deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FillArrays", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables", "TruncatedStacktraces"]
git-tree-sha1 = "9e3e8a2ee407795dd93a8c452080822c6451521f"
uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
version = "2.11.4"

    [deps.SciMLBase.extensions]
    SciMLBaseChainRulesCoreExt = "ChainRulesCore"
    SciMLBasePartialFunctionsExt = "PartialFunctions"
    SciMLBasePyCallExt = "PyCall"
    SciMLBasePythonCallExt = "PythonCall"
    SciMLBaseRCallExt = "RCall"
    SciMLBaseZygoteExt = "Zygote"

    [deps.SciMLBase.weakdeps]
    ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
    PartialFunctions = "570af359-4316-4cb7-8c74-252c00c2016b"
    PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
    PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
    RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
    Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[[deps.SciMLOperators]]
deps = ["ArrayInterface", "DocStringExtensions", "Lazy", "LinearAlgebra", "Setfield", "SparseArrays", "StaticArraysCore", "Tricks"]
git-tree-sha1 = "51ae235ff058a64815e0a2c34b1db7578a06813d"
uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
version = "0.3.7"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"]
git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "1.1.1"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleNonlinearSolve]]
deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "DiffEqBase", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "Reexport", "SciMLBase", "StaticArraysCore"]
git-tree-sha1 = "d3009620f6aa42016c06bc2aaa766ca8bc4ae5db"
uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
version = "1.0.3"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.SimpleUnPack]]
git-tree-sha1 = "58e6353e72cde29b90a69527e56df1b5c3d8c437"
uuid = "ce78b400-467f-4804-87d8-8f486da07d0a"
version = "1.1.0"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "5165dfb9fd131cf0c6957a3a7605dede376e7b63"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.2.0"

[[deps.SparseArrays]]
deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SparseDiffTools]]
deps = ["ADTypes", "Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "Graphs", "LinearAlgebra", "PackageExtensionCompat", "Random", "Reexport", "SciMLOperators", "Setfield", "SparseArrays", "StaticArrayInterface", "StaticArrays", "Tricks", "UnPack", "VertexSafeGraphs"]
git-tree-sha1 = "c281e11db4eacb36a292a054bac83c5a0aca2a26"
uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
version = "2.15.0"

    [deps.SparseDiffTools.extensions]
    SparseDiffToolsEnzymeExt = "Enzyme"
    SparseDiffToolsSymbolicsExt = "Symbolics"
    SparseDiffToolsZygoteExt = "Zygote"

    [deps.SparseDiffTools.weakdeps]
    Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
    Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
    Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[[deps.Sparspak]]
deps = ["Libdl", "LinearAlgebra", "Logging", "OffsetArrays", "Printf", "SparseArrays", "Test"]
git-tree-sha1 = "342cf4b449c299d8d1ceaf00b7a49f4fbc7940e7"
uuid = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
version = "0.3.9"

[[deps.SpecialFunctions]]
deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.3.1"
weakdeps = ["ChainRulesCore"]

    [deps.SpecialFunctions.extensions]
    SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.8.8"

[[deps.StaticArrayInterface]]
deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"]
git-tree-sha1 = "5d66818a39bb04bf328e92bc933ec5b4ee88e436"
uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718"
version = "1.5.0"
weakdeps = ["OffsetArrays", "StaticArrays"]

    [deps.StaticArrayInterface.extensions]
    StaticArrayInterfaceOffsetArraysExt = "OffsetArrays"
    StaticArrayInterfaceStaticArraysExt = "StaticArrays"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
git-tree-sha1 = "fba11dbe2562eecdfcac49a05246af09ee64d055"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.8.1"
weakdeps = ["ChainRulesCore", "Statistics"]

    [deps.StaticArrays.extensions]
    StaticArraysChainRulesCoreExt = "ChainRulesCore"
    StaticArraysStatisticsExt = "Statistics"

[[deps.StaticArraysCore]]
git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
version = "1.4.2"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
version = "1.9.0"

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "1ff449ad350c9c4cbc756624d6f8a8c3ef56d3ed"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.7.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "1d77abd07f617c4868c33d4f5b9e1dbb2643c9cf"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.34.2"

[[deps.StatsFuns]]
deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "f625d686d5a88bcd2b15cd81f18f98186fdc0c9a"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "1.3.0"

    [deps.StatsFuns.extensions]
    StatsFunsChainRulesCoreExt = "ChainRulesCore"
    StatsFunsInverseFunctionsExt = "InverseFunctions"

    [deps.StatsFuns.weakdeps]
    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
    InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"

[[deps.StochasticDiffEq]]
deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqNoiseProcess", "DocStringExtensions", "FiniteDiff", "ForwardDiff", "JumpProcesses", "LevyArea", "LinearAlgebra", "Logging", "MuladdMacro", "NLsolve", "OrdinaryDiffEq", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
git-tree-sha1 = "753219de57ac7aab0feb88871d3c51e0eb5e3b03"
uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
version = "6.64.0"

[[deps.StrideArraysCore]]
deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface", "ThreadingUtilities"]
git-tree-sha1 = "d6415f66f3d89c615929af907fdc6a3e17af0d8c"
uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da"
version = "0.5.2"

[[deps.StringManipulation]]
deps = ["PrecompileTools"]
git-tree-sha1 = "a04cabe79c5f01f4d723cc6704070ada0b9d46d5"
uuid = "892a3eda-7b42-436c-8928-eab12a02cf0e"
version = "0.3.4"

[[deps.StructArrays]]
deps = ["Adapt", "ConstructionBase", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"]
git-tree-sha1 = "0a3db38e4cce3c54fe7a71f831cd7b6194a54213"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.16"

[[deps.StructEquality]]
deps = ["Compat"]
git-tree-sha1 = "192a9f1de3cfef80ab1a4ba7b150bb0e11ceedcf"
uuid = "6ec83bb0-ed9f-11e9-3b4c-2b04cb4e219c"
version = "2.1.0"

[[deps.StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "ca4bccb03acf9faaf4137a9abc1881ed1841aa70"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.10.0"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.SuiteSparse_jll]]
deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
version = "5.10.1+6"

[[deps.SymbolicIndexingInterface]]
git-tree-sha1 = "65f4ed0f9e3125e0836df12c231cea3dd98bb165"
uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
version = "0.3.0"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
version = "1.0.3"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits"]
git-tree-sha1 = "cb76cf677714c095e535e3501ac7954732aeea2d"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.11.1"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
version = "1.10.0"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.ThreadingUtilities]]
deps = ["ManualMemory"]
git-tree-sha1 = "eda08f7e9818eb53661b3deb74e3159460dfbc27"
uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
version = "0.5.2"

[[deps.TranscodingStreams]]
git-tree-sha1 = "1fbeaaca45801b4ba17c251dd8603ef24801dd84"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.10.2"
weakdeps = ["Random", "Test"]

    [deps.TranscodingStreams.extensions]
    TestExt = ["Test", "Random"]

[[deps.TriangularSolve]]
deps = ["CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "LoopVectorization", "Polyester", "Static", "VectorizationBase"]
git-tree-sha1 = "fadebab77bf3ae041f77346dd1c290173da5a443"
uuid = "d5829a12-d9aa-46ab-831f-fb7c9ab06edf"
version = "0.1.20"

[[deps.Tricks]]
git-tree-sha1 = "eae1bb484cd63b36999ee58be2de6c178105112f"
uuid = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775"
version = "0.1.8"

[[deps.TruncatedStacktraces]]
deps = ["InteractiveUtils", "MacroTools", "Preferences"]
git-tree-sha1 = "ea3e54c2bdde39062abf5a9758a23735558705e1"
uuid = "781d530d-4396-4725-bb49-402e4bee1e77"
version = "1.4.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.VectorizationBase]]
deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static", "StaticArrayInterface"]
git-tree-sha1 = "7209df901e6ed7489fe9b7aa3e46fb788e15db85"
uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
version = "0.21.65"

[[deps.VertexSafeGraphs]]
deps = ["Graphs"]
git-tree-sha1 = "8351f8d73d7e880bfc042a8b6922684ebeafb35c"
uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f"
version = "0.2.0"

[[deps.XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
git-tree-sha1 = "801cbe47eae69adc50f36c3caec4758d2650741b"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.12.2+0"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
version = "1.2.13+0"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
version = "5.8.0+0"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
version = "1.52.0+1"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
version = "17.4.0+0"

versioninfo()

Julia Version 1.9.4
Commit 8e5136fa297 (2023-11-14 08:46 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 8 × Apple M2
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, apple-m1)
  Threads: 1 on 4 virtual cores

@ChrisRackauckas
Copy link
Contributor

This has another piece beyond the v6.65.

        begin
            #= /Users/chrisrackauckas/.julia/packages/Decapodes/AQT8L/src/simulation.jl:205 =#
            var"__dynamics_•1" = Decapodes.FixedSizeDiffCache(Vector{Float64}(undef, nparts(mesh, :E)))
            var"__dynamics_•6" = Decapodes.FixedSizeDiffCache(Vector{Float64}(undef, nparts(mesh, :E)))
            var"__•_8_1" = Decapodes.FixedSizeDiffCache(Vector{Float64}(undef, nparts(mesh, :E)))
            var"__•_8_2" = Decapodes.FixedSizeDiffCache(Vector{Float64}(undef, nparts(mesh, :V)))
            __dynamics_ḣ = Decapodes.FixedSizeDiffCache(Vector{Float64}(undef, nparts(mesh, :V)))
        end

those are inappropriately eltyped.

@LilithHafner
Copy link

@ChrisRackauckas, I'm looking into this issue but having trouble making sense of what is going on or how to get started. Do you have any advice on how to get up to speed and contribute effectively to the resolution of this issue?

@ChrisRackauckas
Copy link
Contributor

Take a look at:

sim = eval(gensim(ice_dynamics1D, dimension=1))
fₘ = sim(s, generate)

function f(flow_rate, ice_density, u_init_arr)
  n = 3
  ρ = ice_density
  g = 9.8101
  A = fill(flow_rate, ne(s))
  tₑ = 5e13

  u₀ = ComponentArray(dynamics_h = u_init_arr)
  
  constants_and_parameters = (
    n = n,
    stress_ρ = ρ,
    stress_g = g,
    stress_A = A)
  prob = ODEProblem{true, SciMLBase.FullSpecialize}(fₘ, u₀, (0, tₑ), constants_and_parameters)
  @info("Solving")
  soln = solve(prob, Tsit5())
  @info("Done")

  return soln(tₑ)
end

Some choices have to be made since there's no way that can work. Effectively, sim = eval(gensim(ice_dynamics1D, dimension=1)) needs to make choices about cache generation, the choice it makes is https://github.com/AlgebraicJulia/Decapodes.jl/blob/main/src/simulation.jl#L92 where the element type is derived from https://github.com/AlgebraicJulia/Decapodes.jl/blob/main/src/simulation.jl#L398. However, when inside of a ForwardDiff context it needs to know that it's going to nest duals similar to https://github.com/SciML/PreallocationTools.jl/blob/master/test/core_nesteddual.jl#L54-L62. But as you can see, this version that is fully preallocated needs to know before hand whether it's being differentiated. That of course is simply a requirement because when you're doing AD you're going to need a different sized cache, and so if you're generating all caches before knowing how the function is going to be used that is not possible to do.

Now there's a few ways to do this. One way to do this which I investigated was to move the gensim like:

fₘ = sim(s, generate)

function f(flow_rate, ice_density, u_init_arr)
  n = 3
  ρ = ice_density
  g = 9.8101
  A = fill(flow_rate, ne(s))
  tₑ = 5e13

  u₀ = ComponentArray(dynamics_h = u_init_arr)
  sim = eval(gensim(ice_dynamics1D, dimension=1, stateeltype = eltype(u₀)))
  constants_and_parameters = (
    n = n,
    stress_ρ = ρ,
    stress_g = g,
    stress_A = A)
  prob = ODEProblem{true, SciMLBase.FullSpecialize}(fₘ, u₀, (0, tₑ), constants_and_parameters)
  @info("Solving")
  soln = solve(prob, Tsit5())
  @info("Done")

  return soln(tₑ)
end

The downside to this is that the user has to be aware of what's going on here, and it's not all that fast because the function generation and recompilation isn't fast. Also, you cannot eval that and have to use RuntimeGeneratedFunctions, so it's a bit complex on the user.

Also, the biggest issue, is that https://github.com/SciML/PreallocationTools.jl/blob/master/src/PreallocationTools.jl#L74 the restructure operation is effectively just a reshape in this case, but reshape into building a new ComponentArray seems to have issues because of the ReinterpretedArray. We'd need to find some nice way in order to unsafe wrap that back into a normal array in order to get the right dispatching to fix the u.y syntax, since otherwise it cannot index like a ComponentArray which then is where the errors come from. Doing this naively will just allocate new arrays which then defeats the purpose of all this. This is something I think we should fix reguardless of whether it's used in the solution here.

Now the other approach is to develop caches that are not context-dependent. The simplest thing is just to switch from a FixedSizeDiffCache to a DiffCache. This will resize when it's hit with an element type that is sized to be larger bits than expected. It will throw a warning about the resize the first time but other than that it's good to go. This requires the latest version here SciML/PreallocationTools.jl#90 which was the purpose of it, since you'd then need to change https://github.com/AlgebraicJulia/Decapodes.jl/blob/main/src/simulation.jl#L378 to use the element types, i.e. promote_type(eltype(u),eltype(p)) instead of u. However, it still runs into the issue of ReinterpretedArray on ComponentArray.

So the simplest thing here might be to just swap it to a LazyBufferCache https://github.com/SciML/PreallocationTools.jl?tab=readme-ov-file#lazybuffercache, since that doesn't require knowing anything other than the size and a construction mechanism upfront. Since that uses similar https://github.com/SciML/PreallocationTools.jl/blob/master/src/PreallocationTools.jl#L216 it should be fine. That does mean that the first time it's hit with new cache types it will allocate them on demand and we will be introducing a ~100ns overhead on cache retrieval, but I think for PDEs that should be sufficient. So I would try swapping in a LazyBufferCache and seeing if that's sufficient.

If you really want to get crazy you could introduce a bump allocator https://github.com/MasonProtter/Bumper.jl and that would be a general solution here. But I would recommend trying the LBC first and then benchmark it. If the overhead is low enough in comparison to other operations it's a good stable and robust way of handling the AD nesting so we can leave it at that, though the annoyances of DiffCache have nerd sniped me into wanting to make sure DiffCache+ComponentArrays works better.

@LilithHafner
Copy link

With SciML/PreallocationTools.jl#93 IIUC the only way to get substantial (more than a few dictionary lookups per solve) savings compared to LazyBufferCache is to re-use the same memory for collections of different size and/or type or across different problems, none of which which seem likely to be worth the effort and complexity.

Currently, Decapodes uses get_tmp to get a temporary array from PreallocationTools objects and LazyBufferCache uses getindex. That API difference makes it a slightly non-trivial drop-in replacement. Additionally, FixedSizeDiffCache and DiffCache produce Vectors as output when given CategoricalVectors as input while LazyBufferCache producees a result similar to its input.

IIUC, the disadvantage of using mis-sized DiffCache instead of a correctly sizes FixedSizeDiffCache is mostly a one-off cost of resizing the cache, but in this workflow the cache is accessed many many times, so that doesn't really matter.

TL;DR, switching from FixedSizeDiffCache to DiffCache at https://github.com/AlgebraicJulia/Decapodes.jl/blob/main/src/simulation.jl#L92 "just works" while using LazyBufferCache requires a bit more work.

When I switched FixedSizeDiffCache to DiffCache there the forward pass continued to run successfully in about 12s and the gradient appears to run successfully, but either hangs or takes more than 10 minutes to complete which is likely slower than finite differencing would be. What sort of runtimes do the domain experts expect for computing this gradient?

I have yet to get the forward pass to succeed with LazyBufferCache, though it should be possible, and maybe more performant given the problems with _restructure that @ChrisRackauckas pointed out.

Logs
julia> @time y = f(1e-16, 910., h₀)
[ Info: Solving
[ Info: Done
 12.422065 seconds (4.07 M allocations: 391.523 MiB, 41.71% gc time, 87.68% compilation time)
ComponentVector{Float64}(dynamics_h = [0.31331064770525385, 0.3149084648262946, 0.31753532021923586, 0.32045398600516295, 0.32365951931173176, 0.3270459856949837, 0.3305699140970918, 0.33419138122981706, 0.3378821098663329, 0.3416193021978081  …  0.3416193021978081, 0.3378821098663329, 0.33419138122981706, 0.3305699140970918, 0.3270459856949837, 0.32365951931173176, 0.32045398600516295, 0.31753532021923586, 0.3149084648262946, 0.31331064770525385])

julia> println(time()); @time gradient(x -> f(1e-16, 910., x), h₀)
1.70431743882972e9
[ Info: Solving
┌ Warning: The supplied DiffCache was too small and was enlarged. This incurs allocations
│     on the first call to `get_tmp`. If few calls to `get_tmp` occur and optimal performance is essential,
│     consider changing 'N'/chunk size of this DiffCache to 12.
└ @ PreallocationTools ~/.julia/packages/PreallocationTools/pKtIj/src/PreallocationTools.jl:191
┌ Warning: The supplied DiffCache was too small and was enlarged. This incurs allocations
│     on the first call to `get_tmp`. If few calls to `get_tmp` occur and optimal performance is essential,
│     consider changing 'N'/chunk size of this DiffCache to 12.
└ @ PreallocationTools ~/.julia/packages/PreallocationTools/pKtIj/src/PreallocationTools.jl:191
┌ Warning: The supplied DiffCache was too small and was enlarged. This incurs allocations
│     on the first call to `get_tmp`. If few calls to `get_tmp` occur and optimal performance is essential,
│     consider changing 'N'/chunk size of this DiffCache to 12.
└ @ PreallocationTools ~/.julia/packages/PreallocationTools/pKtIj/src/PreallocationTools.jl:191
^CERROR: InterruptException:                                                  
Stacktrace:
  [1] Array
    @ ./boot.jl:477 [inlined]
  [2] Array
    @ ./boot.jl:486 [inlined]
  [3] Array
    @ ./boot.jl:494 [inlined]
  [4] similar
    @ ./abstractarray.jl:874 [inlined]
  [5] similar
    @ ./abstractarray.jl:873 [inlined]
  [6] similar
    @ ./broadcast.jl:224 [inlined]
  [7] similar
    @ ./broadcast.jl:223 [inlined]
  [8] copy
    @ ./broadcast.jl:928 [inlined]
  [9] materialize
    @ ./broadcast.jl:903 [inlined]
 [10] (::var"#422#435")(x::SubArray{ForwardDiff.Dual{…}, 1, Vector{…}, Tuple{…}, true}, y::Float64)
    @ Main ./REPL[517]:32
 [11] #425
    @ ./REPL[517]:40 [inlined]
 [12] (::var"#f#439"{…})(du::ComponentVector{…}, u::ComponentVector{…}, p::@NamedTuple{…}, t::Float64)
    @ Main ~/.julia/dev/Decapodes/src/simulation.jl:432
 [13] (::ODEFunction{…})(::ComponentVector{…}, ::Vararg{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/2HZ5m/src/scimlfunctions.jl:2355 [inlined]
 [14] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.Tsit5Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/kiJn6/src/perform_step/low_order_rk_perform_step.jl:806
 [15] perform_step!(integrator::Any, cache::OrdinaryDiffEq.Tsit5Cache, repeat_step::Any)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/kiJn6/src/perform_step/low_order_rk_perform_step.jl:798 [inlined]
 [16] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/kiJn6/src/solve.jl:544
 [17] __solve(prob::Union{…}, alg::Union{…}, args::Vararg{…}; kwargs...)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/kiJn6/src/solve.jl:6 [inlined]
 [18] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/kiJn6/src/solve.jl:1 [inlined]
 [19] solve_call(_prob::NonlinearProblem{…}, args::NonlinearSolve.Broyden{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/3QlWZ/src/solve.jl:608 [inlined]
 [20] solve_call
    @ ~/.julia/packages/DiffEqBase/3QlWZ/src/solve.jl:566 [inlined]
 [21] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::ComponentVector{…}, p::@NamedTuple{…}, args::Tsit5{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/3QlWZ/src/solve.jl:1057
 [22] solve_up
    @ DiffEqBase ~/.julia/packages/DiffEqBase/3QlWZ/src/solve.jl:1043 [inlined]
 [23] #solve#40
    @ DiffEqBase ~/.julia/packages/DiffEqBase/3QlWZ/src/solve.jl:980 [inlined]
 [24] solve(prob::ODEProblem{…}, args::Tsit5{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/3QlWZ/src/solve.jl:970
 [25] f(flow_rate::Float64, ice_density::Float64, u_init_arr::Vector{ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 12}})
    @ Main ./REPL[520]:17
 [26] #442
    @ ./REPL[523]:1 [inlined]
 [27] chunk_mode_gradient(f::var"#442#443", x::Vector{…}, cfg::ForwardDiff.GradientConfig{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:123
 [28] gradient(f::Function, x::Vector{…}, cfg::ForwardDiff.GradientConfig{…}, ::Val{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:0
 [29] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{…}, Float64, 12, Vector{…}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:17
 [30] macro expansion
    @ ./timing.jl:279 [inlined]
 [31] top-level scope
    @ ./REPL[523]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> time()
1.704318269410412e9

julia> 1.704318269410412e9-1.70431743882972e9
830.5806920528412

@ChrisRackauckas
Copy link
Contributor

TL;DR, switching from FixedSizeDiffCache to DiffCache at https://github.com/AlgebraicJulia/Decapodes.jl/blob/main/src/simulation.jl#L92 "just works" while using LazyBufferCache requires a bit more work.

Interesting. You didn't hit any reshape issues on ComponentArrays?

When I switched FixedSizeDiffCache to DiffCache there the forward pass continued to run successfully in about 12s and the gradient appears to run successfully, but either hangs or takes more than 10 minutes to complete which is likely slower than finite differencing would be. What sort of runtimes do the domain experts expect for computing this gradient?

It should be under a minute. How are you calculating the gradient? This is the length 30 example?

@LilithHafner
Copy link

No reshape issues, but when I interrupted it threw a stack trace that was allocating which may indicate something isn't working with the pre-allocations.

How are you calculating the gradient? This is the length 30 example?

I pasted the code in the OP into a REPL (excluding the Pkg operations and using a dev'ed version of Decapodes that uses DiffCache instead of FixedSizeDiffCache)

@ChrisRackauckas
Copy link
Contributor

That was after moving the function construction into the loss?

@LilithHafner
Copy link

Before moving. I ran this:

# AlgebraicJulia Dependencies
using Catlab
using Catlab.Graphics
using CombinatorialSpaces
using Decapodes
using ComponentArrays

# External Dependencies
using MLStyle
using MultiScaleArrays
using LinearAlgebra
using OrdinaryDiffEq
using JLD2
using SparseArrays
using Statistics
using GLMakie # Just for visualization
using GeometryBasics: Point2, Point3
Point2D = Point2{Float64};
Point3D = Point3{Float64};

@info("Packages Loaded")

halfar_eq2 = @decapode begin
  h::Form0
  Γ::Form1
  n::Constant== ∂ₜ(h)
  ḣ == (, d, )(Γ * d(h) * avg₀₁(mag((d(h)))^(n-1)) * avg₀₁(h^(n+2)))
end
glens_law = @decapode begin
  Γ::Form1
  (A,ρ,g,n)::Constant

  Γ == (2/(n+2))*A**g)^n
end

@info("Decapodes Defined")

ice_dynamics_composition_diagram = @relation () begin
  dynamics(Γ,n)
  stress(Γ,n)
end
ice_dynamics_cospan = oapply(ice_dynamics_composition_diagram,
  [Open(halfar_eq2, [,:n]),
  Open(glens_law, [,:n])])
ice_dynamics = apex(ice_dynamics_cospan)
ice_dynamics1D = expand_operators(ice_dynamics)
infer_types!(ice_dynamics1D, op1_inf_rules_1D, op2_inf_rules_1D)
resolve_overloads!(ice_dynamics1D, op1_res_rules_1D, op2_res_rules_1D)

s_prime = EmbeddedDeltaSet1D{Bool, Point2D}()
add_vertices!(s_prime, 100, point=Point2D.(range(-2, 2, length=100), 0))
add_edges!(s_prime, 1:nv(s_prime)-1, 2:nv(s_prime))
orient!(s_prime)
s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s_prime)
subdivide_duals!(s, Circumcenter())

@info("Spaces Defined")

function generate(sd, my_symbol; hodge=GeometricHodge())
  e_vecs = map(edges(sd)) do e
    point(sd, sd[e, :∂v0]) - point(sd, sd[e, :∂v1])
  end
  neighbors = map(vertices(sd)) do v
    union(incident(sd, v, :∂v0), incident(sd, v, :∂v1))
  end
  n_vecs = map(neighbors) do es
    [e_vecs[e] for e in es]
  end
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for e in 1:ne(s)
      append!(J, [s[e,:∂v0],s[e,:∂v1]])
      append!(I, [e,e])
      append!(V, [0.5, 0.5])
  end
  avg_mat = sparse(I,J,V)
  op = @match my_symbol begin
    :♯ => x -> begin
      map(neighbors, n_vecs) do es, nvs
        sum([nv*norm(nv)*x[e] for (e,nv) in zip(es,nvs)]) / sum(norm.(nvs))
      end
    end
    :mag => x -> begin
      norm.(x)
    end
    :avg₀₁ => x -> begin
      avg_mat * x
    end
    :^ => (x,y) -> x .^ y
    :* => (x,y) -> x .* y
    :show => x -> begin
      @show x
      x
    end
    x => error("Unmatched operator $my_symbol")
  end
  return (args...) -> op(args...)
end

sim = eval(gensim(ice_dynamics1D, dimension=1))
fₘ = sim(s, generate)

function f(flow_rate, ice_density, u_init_arr)
  n = 3
  ρ = ice_density
  g = 9.8101
  A = fill(flow_rate, ne(s))
  tₑ = 5e13

  u₀ = ComponentArray(dynamics_h = u_init_arr)
  
  constants_and_parameters = (
    n = n,
    stress_ρ = ρ,
    stress_g = g,
    stress_A = A)
  prob = ODEProblem{true, SciMLBase.FullSpecialize}(fₘ, u₀, (0, tₑ), constants_and_parameters)
  @info("Solving")
  soln = solve(prob, Tsit5())
  @info("Done")

  return soln(tₑ)
end


h₀ = map(x -> try sqrt(1. - x[1]^2) catch DomainError return 0.0 end, point(s_prime))

y = f(1e-16, 910., h₀)

using ForwardDiff: gradient;

gradient(x -> f(1e-16, 910., x), h₀)

In this environment

(@DARPA-ASKEM-141) pkg> st
Status `~/.julia/environments/DARPA-ASKEM-141/Project.toml`
  [227ef7b5] ACSets v0.2.12 `https://github.com/AlgebraicJulia/ACSets.jl.git#main`
  [134e5e36] Catlab v0.16.4
  [b1c52339] CombinatorialSpaces v0.6.0
  [b0b7db55] ComponentArrays v0.15.7
  [679ab3ea] Decapodes v0.5.0 `../../dev/Decapodes`
  [f6369f11] ForwardDiff v0.10.36
  [e9467ef8] GLMakie v0.9.4
  [5c1252a2] GeometryBasics v0.4.9
  [033835bb] JLD2 v0.4.41
  [d8e11817] MLStyle v0.4.17
  [f9640e96] MultiScaleArrays v1.12.0
  [1dea7af3] OrdinaryDiffEq v6.67.0
  [37e2e46d] LinearAlgebra
  [2f01184e] SparseArrays v1.10.0
  [10745b16] Statistics v1.10.0

julia> versioninfo()
Julia Version 1.10.0
Commit 3120989f39b (2023-12-25 18:01 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (aarch64-linux-gnu)
  CPU: 8 × unknown
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, generic)
  Threads: 1 on 8 virtual cores

Where the diff from Decapodes main is

-    :($(Symbol(:__,c.name)) = Decapodes.FixedSizeDiffCache(Vector{$(c.T)}(undef, nparts(mesh, $(QuoteNode(resolved_form))))))
+    :($(Symbol(:__,c.name)) = Decapodes.DiffCache(Vector{$(c.T)}(undef, nparts(mesh, $(QuoteNode(resolved_form))))))

@jpfairbanks
Copy link

So it looks like we can apply that patch to main Decapodes and this will work?

@LilithHafner
Copy link

When I applied that patch locally this example hung for me.

@LilithHafner
Copy link

I can no longer reproduce the OP because the latest version of Decapodes is no longer compatible with CombinatorialSpaces v0.6.0. I am investigating that compatibility regression.

methodology
(Decapodes) pkg> activate --temp
  Activating new project at `/tmp/jl_VnYzgL`

(jl_VnYzgL) pkg> add https://github.com/AlgebraicJulia/Decapodes.jl#main
     Cloning git-repo `https://github.com/AlgebraicJulia/Decapodes.jl`
    Updating git-repo `https://github.com/AlgebraicJulia/Decapodes.jl`
   Resolving package versions...
    Updating `/tmp/jl_VnYzgL/Project.toml`
  [679ab3ea] + Decapodes v0.5.0 `https://github.com/AlgebraicJulia/Decapodes.jl#main`
    Updating `/tmp/jl_VnYzgL/Manifest.toml`
  [227ef7b5] + ACSets v0.2.13
  [47edcb42] + ADTypes v0.2.6
⌅ [79e6a3ab] + Adapt v3.7.2
  [23cfdc9f] + AlgebraicInterfaces v0.1.1
⌅ [725a01d3] + AlgebraicRewriting v0.2.1
  [ec485272] + ArnoldiMethod v0.2.0
  [4fba245c] + ArrayInterface v7.7.0
  [4c555306] + ArrayLayouts v1.5.3
  [62783981] + BitTwiddlingConvenienceFunctions v0.1.5
  [2a0fbf3d] + CPUSummary v0.2.4
  [49dc2e85] + Calculus v0.5.1
⌅ [134e5e36] + Catlab v0.15.5
  [d360d2e6] + ChainRulesCore v1.20.1
  [fb6a15b2] + CloseOpenIntervals v0.1.12
  [3da002f7] + ColorTypes v0.11.4
  [5ae59095] + Colors v0.12.10
⌃ [b1c52339] + CombinatorialSpaces v0.5.0
  [861a8166] + Combinatorics v1.0.2
  [38540f10] + CommonSolve v0.2.4
  [bbf7d656] + CommonSubexpressions v0.3.0
  [0fb5dd42] + CompTime v0.1.2
  [34da2185] + Compat v4.12.0
  [b0b7db55] + ComponentArrays v0.15.7
  [a81c6b42] + Compose v0.9.5
  [2569d6c7] + ConcreteStructs v0.2.3
  [187b0558] + ConstructionBase v1.5.4
  [adafc99b] + CpuId v0.3.1
  [a8cc5b0e] + Crayons v4.1.1
  [9a962f9c] + DataAPI v1.16.0
  [864edb3b] + DataStructures v0.18.16
  [e2d170a0] + DataValueInterfaces v1.0.0
  [679ab3ea] + Decapodes v0.5.0 `https://github.com/AlgebraicJulia/Decapodes.jl#main`
  [6f00c28b] + DiagrammaticEquations v0.1.0
  [2b5f629d] + DiffEqBase v6.146.0
  [163ba53b] + DiffResults v1.1.0
  [b552c78f] + DiffRules v1.15.1
  [31c24e10] + Distributions v0.25.107
  [ffbed154] + DocStringExtensions v0.9.3
  [fa6b7ba4] + DualNumbers v0.6.8
  [4e289a0a] + EnumX v1.0.4
  [f151be2c] + EnzymeCore v0.6.5
  [d4d017d3] + ExponentialUtilities v1.25.0
  [e2ba6199] + ExprTools v0.1.10
  [411431e0] + Extents v0.1.2
  [7034ab61] + FastBroadcast v0.2.8
  [9aa1b823] + FastClosures v0.3.2
  [29a986be] + FastLapackInterface v2.0.1
  [5789e2e9] + FileIO v1.16.2
  [1a297f60] + FillArrays v1.9.3
  [6a86dc24] + FiniteDiff v2.22.0
  [53c48c17] + FixedPointNumbers v0.8.4
  [f6369f11] + ForwardDiff v0.10.36
  [069b7b12] + FunctionWrappers v1.1.3
  [77dc65aa] + FunctionWrappersWrappers v0.1.3
  [d9f16b24] + Functors v0.4.5
⌃ [46192b85] + GPUArraysCore v0.1.5
  [6b9d7cbe] + GeneralizedGenerated v0.3.3
  [c145ed77] + GenericSchur v0.5.3
  [cf35fbd7] + GeoInterface v1.3.3
  [5c1252a2] + GeometryBasics v0.4.10
  [86223c79] + Graphs v1.9.0
  [3e5b6fbb] + HostCPUFeatures v0.1.16
  [34004b35] + HypergeometricFunctions v0.3.23
  [615f187c] + IfElse v0.1.1
  [d25df0c9] + Inflate v0.1.4
  [92d709cd] + IrrationalConstants v0.2.2
  [c8e1da08] + IterTools v1.10.0
  [82899510] + IteratorInterfaceExtensions v1.0.0
  [033835bb] + JLD2 v0.4.45
  [692b3bcd] + JLLWrappers v1.5.0
  [682c06a0] + JSON v0.21.4
  [0f8b85d8] + JSON3 v1.14.0
  [b14d175d] + JuliaVariables v0.2.4
  [ef3ab10e] + KLU v0.4.1
  [ba0b0d4f] + Krylov v0.9.5
  [b964fa9f] + LaTeXStrings v1.3.1
  [10f19ff3] + LayoutPointers v0.1.15
  [50d2b5c4] + Lazy v0.15.1
  [5078a376] + LazyArrays v1.8.3
  [9c8b4983] + LightXML v0.9.1
  [d3d80556] + LineSearches v7.2.0
  [7ed4a6bd] + LinearSolve v2.22.1
  [2ab3a3ac] + LogExpFunctions v0.3.26
  [bdcacae8] + LoopVectorization v0.12.166
  [d8e11817] + MLStyle v0.4.17
  [1914dd2f] + MacroTools v0.5.13
  [d125e4d3] + ManualMemory v0.1.8
  [a3b82374] + MatrixFactorizations v2.1.0
  [bb5d69b7] + MaybeInplace v0.1.1
  [442fdcdd] + Measures v0.3.2
  [7269a6da] + MeshIO v0.4.10
  [e1d29d7a] + Missings v1.1.0
  [46d2c3a1] + MuladdMacro v0.2.4
  [d41bc354] + NLSolversBase v7.8.3
  [77ba4419] + NaNMath v1.0.2
  [71a1bf82] + NameResolution v0.1.5
  [8913a72c] + NonlinearSolve v3.5.3
  [6fe1bfb0] + OffsetArrays v1.13.0
  [bac558e1] + OrderedCollections v1.6.3
  [1dea7af3] + OrdinaryDiffEq v6.70.0
  [90014a1f] + PDMats v0.11.31
  [65ce6f38] + PackageExtensionCompat v1.0.2
  [d96e819e] + Parameters v0.12.3
  [69de0a69] + Parsers v2.8.1
  [2ae35dd2] + Permutations v0.4.20
  [f517fe37] + Polyester v0.7.9
  [1d0040c9] + PolyesterWeave v0.2.1
  [d236fae5] + PreallocationTools v0.4.17
  [aea7be01] + PrecompileTools v1.2.0
  [21216c6a] + Preferences v1.4.1
  [8162dcfd] + PrettyPrint v0.2.0
  [08abe8d2] + PrettyTables v2.3.1
  [1fd47b50] + QuadGK v2.9.4
  [3cdcf5f2] + RecipesBase v1.3.4
  [731186ca] + RecursiveArrayTools v3.7.0
  [f2c3362d] + RecursiveFactorization v0.2.21
  [189a3867] + Reexport v1.2.2
  [ae029012] + Requires v1.3.0
  [79098fc4] + Rmath v0.7.1
  [7e49a35a] + RuntimeGeneratedFunctions v0.5.12
  [94e857df] + SIMDTypes v0.1.0
  [476501e8] + SLEEFPirates v0.6.42
  [0bca4576] + SciMLBase v2.23.0
  [c0aeaf25] + SciMLOperators v0.3.7
  [efcf1570] + Setfield v1.1.1
  [727e6d20] + SimpleNonlinearSolve v1.3.1
  [699a6c99] + SimpleTraits v0.9.4
  [ce78b400] + SimpleUnPack v1.1.0
  [a2af1166] + SortingAlgorithms v1.2.1
  [47a9eef4] + SparseDiffTools v2.16.0
  [e56a9233] + Sparspak v0.3.9
  [276daf66] + SpecialFunctions v2.3.1
  [aedffcd0] + Static v0.8.8
  [0d7ed370] + StaticArrayInterface v1.5.0
  [90137ffa] + StaticArrays v1.9.2
  [1e83bf80] + StaticArraysCore v1.4.2
  [82ae8749] + StatsAPI v1.7.0
  [2913bbd2] + StatsBase v0.34.2
  [4c63d2b9] + StatsFuns v1.3.0
  [7792a7ef] + StrideArraysCore v0.5.2
  [892a3eda] + StringManipulation v0.3.4
  [09ab397b] + StructArrays v0.6.17
  [6ec83bb0] + StructEquality v2.1.0
  [856f2bd8] + StructTypes v1.10.0
  [2efcf032] + SymbolicIndexingInterface v0.3.5
  [3783bdb8] + TableTraits v1.0.1
  [bd369af6] + Tables v1.11.1
  [8290d209] + ThreadingUtilities v0.5.2
  [a759f4b9] + TimerOutputs v0.5.23
  [3bb67fe8] + TranscodingStreams v0.10.2
  [d5829a12] + TriangularSolve v0.1.20
  [410a4b4d] + Tricks v0.1.8
  [781d530d] + TruncatedStacktraces v1.4.0
  [3a884ed6] + UnPack v1.0.2
  [3d5dd08c] + VectorizationBase v0.21.65
  [19fa3120] + VertexSafeGraphs v0.2.0
  [5ae413db] + EarCut_jll v2.2.4+0
  [1d5cc7b8] + IntelOpenMP_jll v2024.0.2+0
  [94ce4f54] + Libiconv_jll v1.17.0+0
  [856f044c] + MKL_jll v2024.0.0+0
  [efe28fd5] + OpenSpecFun_jll v0.5.5+0
  [f50d1b31] + Rmath_jll v0.4.0+0
  [02c8fc9c] + XML2_jll v2.12.2+0
  [0dad84c5] + ArgTools v1.1.1
  [56f22d72] + Artifacts
  [2a0f44e3] + Base64
  [ade2ca70] + Dates
  [8ba89e20] + Distributed
  [f43a241f] + Downloads v1.6.0
  [7b1f6079] + FileWatching
  [9fa8497b] + Future
  [b77e0a4c] + InteractiveUtils
  [4af54fe1] + LazyArtifacts
  [b27032c2] + LibCURL v0.6.4
  [76f85450] + LibGit2
  [8f399da3] + Libdl
  [37e2e46d] + LinearAlgebra
  [56ddb016] + Logging
  [d6f4376e] + Markdown
  [a63ad114] + Mmap
  [ca575930] + NetworkOptions v1.2.0
  [44cfe95a] + Pkg v1.10.0
  [de0858da] + Printf
  [3fa0cd96] + REPL
  [9a3f8284] + Random
  [ea8e919c] + SHA v0.7.0
  [9e88b42a] + Serialization
  [1a1011a3] + SharedArrays
  [6462fe0b] + Sockets
  [2f01184e] + SparseArrays v1.10.0
  [10745b16] + Statistics v1.10.0
  [4607b0f0] + SuiteSparse
  [fa267f1f] + TOML v1.0.3
  [a4e569a6] + Tar v1.10.0
  [8dfed614] + Test
  [cf7118a7] + UUIDs
  [4ec0a83e] + Unicode
  [e66e0078] + CompilerSupportLibraries_jll v1.0.5+1
  [deac9b47] + LibCURL_jll v8.4.0+0
  [e37daf67] + LibGit2_jll v1.6.4+0
  [29816b5a] + LibSSH2_jll v1.11.0+1
  [c8ffd9c3] + MbedTLS_jll v2.28.2+1
  [14a3606d] + MozillaCACerts_jll v2023.1.10
  [4536629a] + OpenBLAS_jll v0.3.23+2
  [05823500] + OpenLibm_jll v0.8.1+2
  [bea87d4a] + SuiteSparse_jll v7.2.1+1
  [83775a58] + Zlib_jll v1.2.13+1
  [8e850b90] + libblastrampoline_jll v5.8.0+1
  [8e850ede] + nghttp2_jll v1.52.0+1
  [3f19e933] + p7zip_jll v17.4.0+2
        Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`
┌ Warning: attempting to remove probably stale pidfile
│   path = "/home/x/.julia/compiled/v1.10/Catlab/fBQ1G_WmKiG.ji.pidfile"
└ @ FileWatching.Pidfile ~/.julia/juliaup/julia-1.10.0+0.aarch64.linux.gnu/share/julia/stdlib/v1.10/FileWatching/src/pidfile.jl:244
Precompiling project...
  6 dependencies successfully precompiled in 25 seconds. 207 already precompiled.
  1 dependency had output during precompilation:
┌ DiagrammaticEquations
│  WARNING: using Deca.infer_states in module DiagrammaticEquations conflicts with an existing identifier.
└  

(jl_VnYzgL) pkg> st -m
Status `/tmp/jl_VnYzgL/Manifest.toml`
  [227ef7b5] ACSets v0.2.13
  [47edcb42] ADTypes v0.2.6
⌅ [79e6a3ab] Adapt v3.7.2
  [23cfdc9f] AlgebraicInterfaces v0.1.1
⌅ [725a01d3] AlgebraicRewriting v0.2.1
  [ec485272] ArnoldiMethod v0.2.0
  [4fba245c] ArrayInterface v7.7.0
  [4c555306] ArrayLayouts v1.5.3
  [62783981] BitTwiddlingConvenienceFunctions v0.1.5
  [2a0fbf3d] CPUSummary v0.2.4
  [49dc2e85] Calculus v0.5.1
⌅ [134e5e36] Catlab v0.15.5
  [d360d2e6] ChainRulesCore v1.20.1
  [fb6a15b2] CloseOpenIntervals v0.1.12
  [3da002f7] ColorTypes v0.11.4
  [5ae59095] Colors v0.12.10
⌃ [b1c52339] CombinatorialSpaces v0.5.0
  [861a8166] Combinatorics v1.0.2
  [38540f10] CommonSolve v0.2.4
  [bbf7d656] CommonSubexpressions v0.3.0
  [0fb5dd42] CompTime v0.1.2
  [34da2185] Compat v4.12.0
  [b0b7db55] ComponentArrays v0.15.7
  [a81c6b42] Compose v0.9.5
  [2569d6c7] ConcreteStructs v0.2.3
  [187b0558] ConstructionBase v1.5.4
  [adafc99b] CpuId v0.3.1
  [a8cc5b0e] Crayons v4.1.1
  [9a962f9c] DataAPI v1.16.0
  [864edb3b] DataStructures v0.18.16
  [e2d170a0] DataValueInterfaces v1.0.0
  [679ab3ea] Decapodes v0.5.0 `https://github.com/AlgebraicJulia/Decapodes.jl#main`
  [6f00c28b] DiagrammaticEquations v0.1.0
  [2b5f629d] DiffEqBase v6.146.0
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [31c24e10] Distributions v0.25.107
  [ffbed154] DocStringExtensions v0.9.3
  [fa6b7ba4] DualNumbers v0.6.8
  [4e289a0a] EnumX v1.0.4
  [f151be2c] EnzymeCore v0.6.5
  [d4d017d3] ExponentialUtilities v1.25.0
  [e2ba6199] ExprTools v0.1.10
  [411431e0] Extents v0.1.2
  [7034ab61] FastBroadcast v0.2.8
  [9aa1b823] FastClosures v0.3.2
  [29a986be] FastLapackInterface v2.0.1
  [5789e2e9] FileIO v1.16.2
  [1a297f60] FillArrays v1.9.3
  [6a86dc24] FiniteDiff v2.22.0
  [53c48c17] FixedPointNumbers v0.8.4
  [f6369f11] ForwardDiff v0.10.36
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
  [d9f16b24] Functors v0.4.5
⌃ [46192b85] GPUArraysCore v0.1.5
  [6b9d7cbe] GeneralizedGenerated v0.3.3
  [c145ed77] GenericSchur v0.5.3
  [cf35fbd7] GeoInterface v1.3.3
  [5c1252a2] GeometryBasics v0.4.10
  [86223c79] Graphs v1.9.0
  [3e5b6fbb] HostCPUFeatures v0.1.16
  [34004b35] HypergeometricFunctions v0.3.23
  [615f187c] IfElse v0.1.1
  [d25df0c9] Inflate v0.1.4
  [92d709cd] IrrationalConstants v0.2.2
  [c8e1da08] IterTools v1.10.0
  [82899510] IteratorInterfaceExtensions v1.0.0
  [033835bb] JLD2 v0.4.45
  [692b3bcd] JLLWrappers v1.5.0
  [682c06a0] JSON v0.21.4
  [0f8b85d8] JSON3 v1.14.0
  [b14d175d] JuliaVariables v0.2.4
  [ef3ab10e] KLU v0.4.1
  [ba0b0d4f] Krylov v0.9.5
  [b964fa9f] LaTeXStrings v1.3.1
  [10f19ff3] LayoutPointers v0.1.15
  [50d2b5c4] Lazy v0.15.1
  [5078a376] LazyArrays v1.8.3
  [9c8b4983] LightXML v0.9.1
  [d3d80556] LineSearches v7.2.0
  [7ed4a6bd] LinearSolve v2.22.1
  [2ab3a3ac] LogExpFunctions v0.3.26
  [bdcacae8] LoopVectorization v0.12.166
  [d8e11817] MLStyle v0.4.17
  [1914dd2f] MacroTools v0.5.13
  [d125e4d3] ManualMemory v0.1.8
  [a3b82374] MatrixFactorizations v2.1.0
  [bb5d69b7] MaybeInplace v0.1.1
  [442fdcdd] Measures v0.3.2
  [7269a6da] MeshIO v0.4.10
  [e1d29d7a] Missings v1.1.0
  [46d2c3a1] MuladdMacro v0.2.4
  [d41bc354] NLSolversBase v7.8.3
  [77ba4419] NaNMath v1.0.2
  [71a1bf82] NameResolution v0.1.5
  [8913a72c] NonlinearSolve v3.5.3
  [6fe1bfb0] OffsetArrays v1.13.0
  [bac558e1] OrderedCollections v1.6.3
  [1dea7af3] OrdinaryDiffEq v6.70.0
  [90014a1f] PDMats v0.11.31
  [65ce6f38] PackageExtensionCompat v1.0.2
  [d96e819e] Parameters v0.12.3
  [69de0a69] Parsers v2.8.1
  [2ae35dd2] Permutations v0.4.20
  [f517fe37] Polyester v0.7.9
  [1d0040c9] PolyesterWeave v0.2.1
  [d236fae5] PreallocationTools v0.4.17
  [aea7be01] PrecompileTools v1.2.0
  [21216c6a] Preferences v1.4.1
  [8162dcfd] PrettyPrint v0.2.0
  [08abe8d2] PrettyTables v2.3.1
  [1fd47b50] QuadGK v2.9.4
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.7.0
  [f2c3362d] RecursiveFactorization v0.2.21
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.0
  [79098fc4] Rmath v0.7.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.12
  [94e857df] SIMDTypes v0.1.0
  [476501e8] SLEEFPirates v0.6.42
  [0bca4576] SciMLBase v2.23.0
  [c0aeaf25] SciMLOperators v0.3.7
  [efcf1570] Setfield v1.1.1
  [727e6d20] SimpleNonlinearSolve v1.3.1
  [699a6c99] SimpleTraits v0.9.4
  [ce78b400] SimpleUnPack v1.1.0
  [a2af1166] SortingAlgorithms v1.2.1
  [47a9eef4] SparseDiffTools v2.16.0
  [e56a9233] Sparspak v0.3.9
  [276daf66] SpecialFunctions v2.3.1
  [aedffcd0] Static v0.8.8
  [0d7ed370] StaticArrayInterface v1.5.0
  [90137ffa] StaticArrays v1.9.2
  [1e83bf80] StaticArraysCore v1.4.2
  [82ae8749] StatsAPI v1.7.0
  [2913bbd2] StatsBase v0.34.2
  [4c63d2b9] StatsFuns v1.3.0
  [7792a7ef] StrideArraysCore v0.5.2
  [892a3eda] StringManipulation v0.3.4
  [09ab397b] StructArrays v0.6.17
  [6ec83bb0] StructEquality v2.1.0
  [856f2bd8] StructTypes v1.10.0
  [2efcf032] SymbolicIndexingInterface v0.3.5
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.11.1
  [8290d209] ThreadingUtilities v0.5.2
  [a759f4b9] TimerOutputs v0.5.23
  [3bb67fe8] TranscodingStreams v0.10.2
  [d5829a12] TriangularSolve v0.1.20
  [410a4b4d] Tricks v0.1.8
  [781d530d] TruncatedStacktraces v1.4.0
  [3a884ed6] UnPack v1.0.2
  [3d5dd08c] VectorizationBase v0.21.65
  [19fa3120] VertexSafeGraphs v0.2.0
  [5ae413db] EarCut_jll v2.2.4+0
  [1d5cc7b8] IntelOpenMP_jll v2024.0.2+0
  [94ce4f54] Libiconv_jll v1.17.0+0
  [856f044c] MKL_jll v2024.0.0+0
  [efe28fd5] OpenSpecFun_jll v0.5.5+0
  [f50d1b31] Rmath_jll v0.4.0+0
  [02c8fc9c] XML2_jll v2.12.2+0
  [0dad84c5] ArgTools v1.1.1
  [56f22d72] Artifacts
  [2a0f44e3] Base64
  [ade2ca70] Dates
  [8ba89e20] Distributed
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching
  [9fa8497b] Future
  [b77e0a4c] InteractiveUtils
  [4af54fe1] LazyArtifacts
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2
  [8f399da3] Libdl
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [d6f4376e] Markdown
  [a63ad114] Mmap
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.10.0
  [de0858da] Printf
  [3fa0cd96] REPL
  [9a3f8284] Random
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization
  [1a1011a3] SharedArrays
  [6462fe0b] Sockets
  [2f01184e] SparseArrays v1.10.0
  [10745b16] Statistics v1.10.0
  [4607b0f0] SuiteSparse
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test
  [cf7118a7] UUIDs
  [4ec0a83e] Unicode
  [e66e0078] CompilerSupportLibraries_jll v1.0.5+1
  [deac9b47] LibCURL_jll v8.4.0+0
  [e37daf67] LibGit2_jll v1.6.4+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.2+1
  [14a3606d] MozillaCACerts_jll v2023.1.10
  [4536629a] OpenBLAS_jll v0.3.23+2
  [05823500] OpenLibm_jll v0.8.1+2
  [bea87d4a] SuiteSparse_jll v7.2.1+1
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.8.0+1
  [8e850ede] nghttp2_jll v1.52.0+1
  [3f19e933] p7zip_jll v17.4.0+2
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`

@LilithHafner
Copy link

LilithHafner commented Feb 5, 2024

I can once more reproduce the OP, though I now need to add the statement

using Decapodes: @decapode, Open, expand_operators, infer_types!, resolve_overloads!,
  op1_inf_rules_1D, op2_inf_rules_1D, op1_res_rules_1D, op2_res_rules_1D

as these names appear to no longer be exported.

@lukem12345
Copy link

@LilithHafner These ought to be exported when using using DiagrammaticEquations and using DiagrammaticEquations.Deca. Does this work on your environment?

@LilithHafner
Copy link

Yes. The OP with

using DiagrammaticEquations
using DiagrammaticEquations.Deca

added also reproduces for me.

@LilithHafner
Copy link

Here are some solutions I've tried which did not work:

  1. FixedSizeDiffCache gives the error in the OP

  2. DiffCache succeeds in the forward computation but says "Warning: The supplied DiffCache was too small and was enlarged." and hangs computing the gradient. The hang cannot be easily interrupted.

  3. LazyBufferCache does not support get_tmp, but if we change get_tmp to LazyBufferCache's getindex API we get "expected ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(dynamics_h = 1:100,)}}}, got a value of type Vector{Float64}" because LazyBufferCache looks at eltype and size but not at the full type of it's input.

  4. Using this custom cache:

struct MyBufferCache
    cache::Dict{Any, Any}
end
MyBufferCache() = MyBufferCache(Dict())
function get_tmp(cache::MyBufferCache, x::AbstractArray)
    get(cache.cache, (axes(x), typeof(x))) do
        copy(x)
    end
end

Gives DimensionMismatch on mul!(var"dynamics_•1", var"GenSim-M_d₀", dynamics_h) because of [^1]

  1. Using this custom cache
struct MyFixedSizeBufferCache{N}
    cache::Dict{Any, Any}
    size::NTuple{N, Int}
end
MyFixedSizeBufferCache(size...) = MyFixedSizeBufferCache(Dict(), size)
function PreallocationTools.get_tmp(cache::MyFixedSizeBufferCache, x)
    get(cache.cache, typeof(x)) do
        similar(x, cache.size...)
    end
end

Succeeds in the forward solve and hangs computing the gradient, though keyboard interruption works to stop the hang


I plan to investigate why the error using FixedSizeDiffCache is thrown and why the warning using DiffCache is printed. Happy to hear ideas from other folks who might know what's going on or be able to interpret these results.


[^1] A simplified excerpt of the generated dynamics is

dynamics_h = u.dynamics_h
...
var"dynamics_•1" = Decapodes.get_tmp(var"__dynamics_•1", u)
...
# var"GenSim-M_d₀" is a 99x100 matrix
mul!(var"dynamics_•1", var"GenSim-M_d₀", dynamics_h) # <=== This line throws
...
dynamics_mult_3 .= Γ .* var"dynamics_•1"
dynamics_mult_1 = dynamics_mult_3 .* var"dynamics_•2"
dynamics_mult_2 = dynamics_mult_1 .* var"dynamics_•8"
# var"GenSim-M_GenSim-ConMat_1" is a 100x99 matrix
mul!(dynamics_ḣ, var"GenSim-M_GenSim-ConMat_1", dynamics_mult_2)

I think the issue is that all the temporary variables are constructed in terms of u
which does not provide enough information to determine what the size of the scratch
spaces should be. For example, in this case, u.dynamics_h has length 100, but the
intermediary var"dynamics_•1" needs length 99. We could do some really hacky things at the
"get_tmp" level to make this sort of work, but an efficient robust solution should track the
requisite size of all intermediary variables and pre-allocate the correct size.

The existing implementation does that by storing dims in the cache object and MyBufferCache
does not track this size.

@LilithHafner
Copy link

When I change the caching system to avoid the errors in the gradient computation, I'm seeing the dynamics run substantially more times computing the gradient than computing the forward pass. Specifically, thy run between 13k and 27k times on the forward pass and more than 16 million times on the reverse pass. From @ChrisRackauckas "It should be under a minute. " I'm guessing this 500x increase is not expected. Do folks have any Idea why the dynamics would run so many times?

@ChrisRackauckas
Copy link
Contributor

@LilithHafner Is it in mutating form? Are you getting a fallback to ReverseDiff non-compiled mode for the vjp calculation? It should be throwing warnings about this.

@ChrisRackauckas
Copy link
Contributor

We could remove the caching and just have a less memory efficient but simpler runtime if that helps.

Having two dispatches could be helpful for the reverse mode. Reverse mode in some sense "necessarily" allocates in the vjp, though with enough sophistication we can work around that, though my gut says there's a warning being thrown in the hasbranching check

@LilithHafner
Copy link

Do you think that if we removed caching entirely we would see the increased calls to the dynamics?

Removing caching should have not semantic impact whatsoever (unless there is a caching bug like simultaneous uses of a temp from the same cache) I'll check though.

Is it in mutating form?

Yes. The dynamic function is defined as

#= /home/x/.julia/dev/Decapodes/src/simulation.jl:597 =#
f(du, u, p, t) = begin
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:597 =#
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:598 =#
        [...]
        mul!(dynamics_ḣ, var"GenSim-M_GenSim-ConMat_1", dynamics_mult_2)
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:600 =#
        getproperty(du, :dynamics_h) .= dynamics_ḣ
    end

Are you getting a fallback to ReverseDiff non-compiled mode for the vjp calculation? It should be throwing warnings about this.

It's not throwing warnings, so probably not?

@LilithHafner
Copy link

Replacing all uses of

a = get_temp(...)
mul!(a, b, c)

with

a = b * c

and all uses of

a = get_temp(...)
a .= ...

with

a = ...

Has no impact.


Also switching f from the 4-arg mutating mode to a three-arg version that returns a freshly allocated gradient also does not help with computing the gradient.


Removing all this caching increases the forward pass runtime by 1.4x from 92ms to 130ms.

I think that (with some upstream fixes, too) switching from a FixedSizeDiffCache to a LazyBufferCache will fix the caching issue, but I think there is also a non-caching-related issue which causes gradient computation to either hang or run for a long time.

@ChrisRackauckas
Copy link
Contributor

Which vjp method is chosen for the reverse mode?

@LilithHafner
Copy link

How can I find out which vjp method is chosen for the reverse mode? (this is also ForwardDiff.gradient, so maybe there is no reverse mode?) Here's a stack trace. It includes a frame chunk_mode_gradient if that's what you're looking for.

Stacktrace:
  [1] (::var"#f#373"{…})(du::ComponentVector{…}, u::ComponentVector{…}, p::@NamedTuple{…}, t::Float64)
    @ Main ~/.julia/dev/Decapodes/src/simulation.jl:579
  [2] ODEFunction
    @ ~/.julia/packages/SciMLBase/Avnpi/src/scimlfunctions.jl:2180 [inlined]
  [3] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.Tsit5Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/pWCEq/src/perform_step/low_order_rk_perform_step.jl:821
  [4] perform_step!(integrator::Any, cache::OrdinaryDiffEq.Tsit5Cache, repeat_step::Any)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/pWCEq/src/perform_step/low_order_rk_perform_step.jl:798 [inlined]
  [5] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/pWCEq/src/solve.jl:538
  [6] __solve(prob::Union{…}, alg::Union{…}, args::Vararg{…}; kwargs...)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/pWCEq/src/solve.jl:6 [inlined]
  [7] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/pWCEq/src/solve.jl:1 [inlined]
  [8] solve_call(_prob::Any, args::Vararg{Any}; merge_callbacks::Any, kwargshandle::Any, kwargs...)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:609 [inlined]
 [9] solve_call
    @ ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:567 [inlined]
 [10] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::ComponentVector{…}, p::@NamedTuple{…}, args::Tsit5{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:1058
 [11] solve_up
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:1044 [inlined]
 [12] #solve#40
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:981 [inlined]
 [13] solve(prob::ODEProblem{…}, args::Tsit5{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:971
 [14] f(flow_rate::Float64, ice_density::Float64, u_init_arr::Vector{ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 12}})
    @ Main ~/.julia/dev/decapodes-issue/example.jl:131
 [15] #376
    @ ~/.julia/dev/decapodes-issue/example.jl:145 [inlined]
 [16] chunk_mode_gradient(f::var"#376#377", x::Vector{…}, cfg::ForwardDiff.GradientConfig{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:123
 [17] gradient(f::Function, x::Vector{…}, cfg::ForwardDiff.GradientConfig{…}, ::Val{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:0
 [18] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{…}, Float64, 12, Vector{…}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:17
 [19] top-level scope
    @ ~/.julia/dev/decapodes-issue/example.jl:145

@ChrisRackauckas
Copy link
Contributor

Oh if it's ForwardDiff then it's just forward mode and going to be slow. That should be changed to a Zygote call.

@ChrisRackauckas
Copy link
Contributor

The time for forward mode is O((np)), so the scaling will be extremely bad with that ( O((np)^3) with an implicit method)

@LilithHafner
Copy link

LilithHafner commented Feb 10, 2024

That's fantastic news! I'll verify that this works well with Zygote and we should be good to go (with one or two more forthcoming upstream PRs)

EDIT: I'm seeing a new batch of issues with Zygote. I'll investigate more on Monday.

@LilithHafner
Copy link

Zygote can't AD this because of mutation. It seems like folks have made a fairly effective fallback system for when AD fails, but that system is still slow and has some inherent limitations.

I agree with @ChrisRackauckas's comment that

Having two dispatches could be helpful for the reverse mode.

Specifically, I think we should add a non-mutating option.

If anyone is curious, the latest error I'm dealing with is

┌ Warning: Potential performance improvement omitted. EnzymeVJP tried and failed in the automated AD choice algorithm. To show the stack trace, set SciMLSensitivity.STACKTRACE_WITH_VJPWARN[] = true. To turn off this printing, add `verbose = false` to the `solve` call.
└ @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/concrete_solve.jl:24

┌ Warning: Potential performance improvement omitted. ReverseDiffVJP tried and failed in the automated AD choice algorithm. To show the stack trace, set SciMLSensitivity.STACKTRACE_WITH_VJPWARN[] = true. To turn off this printing, add `verbose = false` to the `solve` call.
└ @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/concrete_solve.jl:65

┌ Warning: Reverse-Mode AD VJP choices all failed. Falling back to numerical VJPs
└ @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/concrete_solve.jl:167
[ Info: Done
ERROR: LoadError: BoundsError: attempt to access 1321-element Vector{Float64} at index [0]
Stacktrace:
  [1] getindex
    @ ./essentials.jl:13 [inlined]
  [2] (::SciMLSensitivity.ReverseLossCallback{…})(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/adjoint_common.jl:514
  [3] #111
    @ ~/.julia/packages/DiffEqCallbacks/uVI0B/src/preset_time.jl:32 [inlined]
  [4] apply_discrete_callback!
    @ ~/.julia/packages/DiffEqBase/UoaYd/src/callbacks.jl:605 [inlined]
  [5] handle_callbacks!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/integrators/integrator_utils.jl:346
  [6] _loopfooter!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/integrators/integrator_utils.jl:253
  [7] loopfooter!
    @ ~/.julia/dev/OrdinaryDiffEq/src/integrators/integrator_utils.jl:206 [inlined]
  [8] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:539
  [9] __solve(prob::Union{…}, alg::Union{…}, args::Vararg{…}; kwargs...)
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:6 [inlined]
 [10] __solve
    @ ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:1 [inlined]
 [11] solve_call(_prob::ODEProblem{…}, args::Tsit5{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:609
 [12] solve_call
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:567 [inlined]
 [13] #solve_up#42
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:1058 [inlined]
 [14] solve_up
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:1044 [inlined]
 [15] #solve#40
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:981 [inlined]
 [16] _adjoint_sensitivities(sol::ODESolution{…}, sensealg::InterpolatingAdjoint{…}, alg::Tsit5{…}; t::Vector{…}, dgdu_discrete::Function, dgdp_discrete::Nothing, dgdu_continuous::Nothing, dgdp_continuous::Nothing, g::Nothing, abstol::Float64, reltol::Float64, checkpoints::Vector{…}, corfunc_analytical::Nothing, callback::Nothing, kwargs::@Kwargs{…})
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/sensitivity_interface.jl:432
 [17] _adjoint_sensitivities
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/sensitivity_interface.jl:390 [inlined]
 [18] #adjoint_sensitivities#63
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/sensitivity_interface.jl:386 [inlined]
 [19] (::SciMLSensitivity.var"#adjoint_sensitivity_backpass#307"{…})(Δ::ODESolution{…})
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/concrete_solve.jl:536
 [20] ZBack
    @ ~/.julia/packages/Zygote/jxHJc/src/compiler/chainrules.jl:211 [inlined]
 [21] (::Zygote.var"#291#292"{…})(Δ::ODESolution{…})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/lib.jl:206
 [22] (::Zygote.var"#2169#back#293"{…})(Δ::ODESolution{…})
    @ Zygote ~/.julia/packages/ZygoteRules/M4xmc/src/adjoint.jl:72
 [23] #solve#40
    @ ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:981 [inlined]
 [24] (::Zygote.Pullback{…})(Δ::ODESolution{…})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [25] (::Zygote.var"#291#292")(Δ::Any)
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/lib.jl:206 [inlined]
 [26] #2169#back
    @ ~/.julia/packages/ZygoteRules/M4xmc/src/adjoint.jl:72 [inlined]
 [27] solve
    @ ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:971 [inlined]
 [28] (::Zygote.Pullback{…})(Δ::ODESolution{…})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [29] f
    @ ~/.julia/dev/decapodes-issue/example.jl:132 [inlined]
 [30] (::Zygote.Pullback{Tuple{typeof(f), Float64, Float64, Vector{Float64}}, Any})(Δ::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [31] #210
    @ Main ~/.julia/dev/decapodes-issue/example.jl:148 [inlined]
 [32] #291
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/lib.jl:206 [inlined]
 [33] (::Zygote.var"#2169#back#293"{Zygote.var"#291#292"{Tuple{…}, Zygote.Pullback{…}}})(Δ::Vector{Float64})
    @ Zygote ~/.julia/packages/ZygoteRules/M4xmc/src/adjoint.jl:72
 [34] call_composed
    @ ./operators.jl:1045 [inlined]
 [35] (::Zygote.Pullback{Tuple{typeof(Base.call_composed), Tuple{…}, Tuple{…}, @Kwargs{}}, Any})(Δ::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [36] call_composed
    @ ./operators.jl:1044 [inlined]
 [37] #_#103
    @ ./operators.jl:1041 [inlined]
 [38] (::Zygote.Pullback{Tuple{…}, Tuple{…}})(Δ::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [39] (::Zygote.var"#291#292")(Δ::Any)
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/lib.jl:206 [inlined]
 [40] #2169#back
    @ ~/.julia/packages/ZygoteRules/M4xmc/src/adjoint.jl:72 [inlined]
 [41] ComposedFunction
    @ ./operators.jl:1041 [inlined]
 [42] (::Zygote.Pullback{Tuple{…}, Tuple{…}})(Δ::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [43] (::Zygote.var"#75#76"{Zygote.Pullback{Tuple{…}, Tuple{…}}})(Δ::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface.jl:91
 [44] withjacobian(f::Function, args::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/grad.jl:150
 [45] jacobian(f::Function, args::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/grad.jl:128
 [46] top-level scope
    @ ~/.julia/dev/decapodes-issue/example.jl:148
 [47] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [48] top-level scope
    @ REPL[12]:1
in expression starting at /home/x/.julia/dev/decapodes-issue/example.jl:148
Some type information was truncated. Use `show(err)` to see complete types.

Using the code

# AlgebraicJulia Dependencies
using Catlab
using Catlab.Graphics
using CombinatorialSpaces
using Decapodes
using ComponentArrays

# External Dependencies
using MLStyle
using MultiScaleArrays
using LinearAlgebra
using OrdinaryDiffEq
using JLD2
using SparseArrays
using Statistics
using GLMakie # Just for visualization
using GeometryBasics: Point2, Point3
Point2D = Point2{Float64};
Point3D = Point3{Float64};

using DiagrammaticEquations
using DiagrammaticEquations.Deca

@info("Packages Loaded")

halfar_eq2 = @decapode begin
  h::Form0
  Γ::Form1
  n::Constant== ∂ₜ(h)
  ḣ == (, d, )(Γ * d(h) * avg₀₁(mag((d(h)))^(n-1)) * avg₀₁(h^(n+2)))
end
glens_law = @decapode begin
  Γ::Form1
  (A,ρ,g,n)::Constant

  Γ == (2/(n+2))*A**g)^n
end

@info("Decapodes Defined")

ice_dynamics_composition_diagram = @relation () begin
  dynamics(Γ,n)
  stress(Γ,n)
end
ice_dynamics_cospan = oapply(ice_dynamics_composition_diagram,
  [Open(halfar_eq2, [,:n]),
  Open(glens_law, [,:n])])
ice_dynamics = apex(ice_dynamics_cospan)
ice_dynamics1D = expand_operators(ice_dynamics)
infer_types!(ice_dynamics1D, op1_inf_rules_1D, op2_inf_rules_1D)
resolve_overloads!(ice_dynamics1D, op1_res_rules_1D, op2_res_rules_1D)

s_prime = EmbeddedDeltaSet1D{Bool, Point2D}()
add_vertices!(s_prime, 100, point=Point2D.(range(-2, 2, length=100), 0))
add_edges!(s_prime, 1:nv(s_prime)-1, 2:nv(s_prime))
orient!(s_prime)
s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s_prime)
subdivide_duals!(s, Circumcenter())

@info("Spaces Defined")

function generate(sd, my_symbol; hodge=GeometricHodge())
  e_vecs = map(edges(sd)) do e
    point(sd, sd[e, :∂v0]) - point(sd, sd[e, :∂v1])
  end
  neighbors = map(vertices(sd)) do v
    union(incident(sd, v, :∂v0), incident(sd, v, :∂v1))
  end
  n_vecs = map(neighbors) do es
    [e_vecs[e] for e in es]
  end
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for e in 1:ne(s)
      append!(J, [s[e,:∂v0],s[e,:∂v1]])
      append!(I, [e,e])
      append!(V, [0.5, 0.5])
  end
  avg_mat = sparse(I,J,V)
  op = @match my_symbol begin
    :♯ => x -> begin
      map(neighbors, n_vecs) do es, nvs
        sum([nv*norm(nv)*x[e] for (e,nv) in zip(es,nvs)]) / sum(norm.(nvs))
      end
    end
    :mag => x -> begin
      norm.(x)
    end
    :avg₀₁ => x -> begin
      avg_mat * x
    end
    :^ => (x,y) -> max.(x, 0) .^ y # x is sometimes a very small negative number due to numerical error
    :* => (x,y) -> x .* y
    :show => x -> begin
      @show x
      x
    end
    x => error("Unmatched operator $my_symbol")
  end
  return (args...) -> op(args...)
end

sim = eval(gensim(ice_dynamics1D, dimension=1))

fₘ = sim(s, generate)

function f(flow_rate, ice_density, u_init_arr)
  n = 3
  ρ = ice_density
  g = 9.8101
  A = fill(flow_rate, ne(s))
  tₑ = 5e13

  u₀ = ComponentArray(dynamics_h = u_init_arr)

  # Note that this must be a ComponentArray to differentiate
  constants_and_parameters = ComponentArray(
    n = n,
    stress_ρ = ρ,
    stress_g = g,
    stress_A = A)
  prob = ODEProblem{true, SciMLBase.FullSpecialize}(fₘ, u₀, (0, tₑ), constants_and_parameters)
  @info("Solving")
  soln = solve(prob, Tsit5())
  @info("Done")

  # return soln(tₑ)
  last(soln) # last, not soln(tₑ) because to avoid interpolation fails when AD fails.
end


h₀ = map(x -> try sqrt(1. - x[1]^2) catch DomainError return 0.0 end, point(s_prime))

y = f(1e-16, 910., h₀)

# Must use Zygote—ForwardDiff is asymtotically slow
using SciMLSensitivity
using Zygote

Zygote.jacobian(x -> f(1e-16, 910., x), h₀)

With Decapodes configured to use LazyBufferCache and the patch from SciML/PreallocationTools.jl#102.

@ChrisRackauckas, would you like me to look into adding an option to Decapodes to generate mutation-free code?

@jpfairbanks
Copy link

I think that would be really helpful for benchmarking too. I definitely expected that our preallocations were saving more than 1.4X in runtime but we didn't implement a nonpreallocating version because we just assumed that memory traffic would be a significant performance problem that needed a solution. Removing the preallocation would make the code a lot more flexible and if it only costs 2X in runtime it is worth supporting that version too.

We probably want to use a compiler configuration abstract type and dispatch rather than a bunch of kwargs and if/else to future proof as we might need more compiler options.

@ChrisRackauckas
Copy link
Contributor

Change

  soln = solve(prob, Tsit5())
  @info("Done")

  return soln(tₑ)

to

  soln = solve(prob, Tsit5(), saveat = tₑ)
  return Array(soln)

@LilithHafner
Copy link

we just assumed that memory traffic would be a significant performance problem that needed a solution

Shout out to the folks who make the Julia GC

Change...

After that change I get

[ Info: Packages Loaded
[ Info: Decapodes Defined
[ Info: Spaces Defined
[ Info: Solving
[ Info: Solving
┌ Warning: Potential performance improvement omitted. EnzymeVJP tried and failed in the automated AD choice algorithm. To show the stack trace, set SciMLSensitivity.STACKTRACE_WITH_VJPWARN[] = true. To turn off this printing, add `verbose = false` to the `solve` call.
└ @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/concrete_solve.jl:24

┌ Warning: Potential performance improvement omitted. ReverseDiffVJP tried and failed in the automated AD choice algorithm. To show the stack trace, set SciMLSensitivity.STACKTRACE_WITH_VJPWARN[] = true. To turn off this printing, add `verbose = false` to the `solve` call.
└ @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/concrete_solve.jl:65

┌ Warning: Reverse-Mode AD VJP choices all failed. Falling back to numerical VJPs
└ @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/concrete_solve.jl:167
ERROR: LoadError: BoundsError: attempt to access 2-element StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64} at index [0]
Stacktrace:
  [1] throw_boundserror(A::StepRangeLen{Float64, Base.TwicePrecision{…}, Base.TwicePrecision{…}, Int64}, I::Tuple{Int64})
    @ Base ./abstractarray.jl:734
  [2] checkbounds
    @ ./abstractarray.jl:699 [inlined]
  [3] getindex
    @ ./range.jl:955 [inlined]
  [4] (::SciMLSensitivity.ReverseLossCallback{…})(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/adjoint_common.jl:514
  [5] #111
    @ ~/.julia/packages/DiffEqCallbacks/uVI0B/src/preset_time.jl:32 [inlined]
  [6] apply_discrete_callback!
    @ ~/.julia/packages/DiffEqBase/UoaYd/src/callbacks.jl:605 [inlined]
  [7] handle_callbacks!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/integrators/integrator_utils.jl:346
  [8] _loopfooter!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/integrators/integrator_utils.jl:253
  [9] loopfooter!
    @ ~/.julia/dev/OrdinaryDiffEq/src/integrators/integrator_utils.jl:206 [inlined]
 [10] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:539
 [11] __solve(prob::Union{…}, alg::Union{…}, args::Vararg{…}; kwargs...)
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:6 [inlined]
 [12] __solve
    @ ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:1 [inlined]
 [13] solve_call(_prob::ODEProblem{…}, args::Tsit5{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:609
 [14] solve_call
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:567 [inlined]
 [15] #solve_up#42
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:1058 [inlined]
 [16] solve_up
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:1044 [inlined]
 [17] #solve#40
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:981 [inlined]
 [18] _adjoint_sensitivities(sol::ODESolution{…}, sensealg::InterpolatingAdjoint{…}, alg::Tsit5{…}; t::StepRangeLen{…}, dgdu_discrete::Function, dgdp_discrete::Nothing, dgdu_continuous::Nothing, dgdp_continuous::Nothing, g::Nothing, abstol::Float64, reltol::Float64, checkpoints::Vector{…}, corfunc_analytical::Nothing, callback::Nothing, kwargs::@Kwargs{…})
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/sensitivity_interface.jl:432
 [19] _adjoint_sensitivities
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/sensitivity_interface.jl:390 [inlined]
 [20] #adjoint_sensitivities#63
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/sensitivity_interface.jl:386 [inlined]
 [21] (::SciMLSensitivity.var"#adjoint_sensitivity_backpass#307"{…})(Δ::ODESolution{…})
    @ SciMLSensitivity ~/.julia/packages/SciMLSensitivity/450Qi/src/concrete_solve.jl:536
 [22] ZBack
    @ ~/.julia/packages/Zygote/jxHJc/src/compiler/chainrules.jl:211 [inlined]
 [23] (::Zygote.var"#kw_zpullback#53"{…})(dy::ODESolution{…})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/chainrules.jl:237
 [24] #291
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/lib.jl:206 [inlined]
 [25] (::Zygote.var"#2169#back#293"{…})(Δ::ODESolution{…})
    @ Zygote ~/.julia/packages/ZygoteRules/M4xmc/src/adjoint.jl:72
 [26] #solve#40
    @ ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:981 [inlined]
 [27] (::Zygote.Pullback{…})(Δ::ODESolution{…})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [28] #291
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/lib.jl:206 [inlined]
 [29] (::Zygote.var"#2169#back#293"{…})(Δ::ODESolution{…})
    @ Zygote ~/.julia/packages/ZygoteRules/M4xmc/src/adjoint.jl:72
 [30] solve
    @ ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:971 [inlined]
 [31] (::Zygote.Pullback{…})(Δ::ODESolution{…})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [32] f
    @ ~/.julia/dev/decapodes-issue/example.jl:127 [inlined]
 [33] (::Zygote.Pullback{Tuple{typeof(f), Float64, Float64, Vector{Float64}}, Any})(Δ::Matrix{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [34] #42
    @ Main ~/.julia/dev/decapodes-issue/example.jl:144 [inlined]
 [35] #291
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/lib.jl:206 [inlined]
 [36] (::Zygote.var"#2169#back#293"{Zygote.var"#291#292"{Tuple{…}, Zygote.Pullback{…}}})(Δ::Matrix{Float64})
    @ Zygote ~/.julia/packages/ZygoteRules/M4xmc/src/adjoint.jl:72
 [37] call_composed
    @ ./operators.jl:1045 [inlined]
 [38] (::Zygote.Pullback{Tuple{typeof(Base.call_composed), Tuple{…}, Tuple{…}, @Kwargs{}}, Any})(Δ::Matrix{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [39] call_composed
    @ ./operators.jl:1044 [inlined]
 [40] #_#103
    @ ./operators.jl:1041 [inlined]
 [41] (::Zygote.Pullback{Tuple{…}, Tuple{…}})(Δ::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [42] (::Zygote.var"#291#292")(Δ::Any)
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/lib.jl:206 [inlined]
 [43] #2169#back
    @ ~/.julia/packages/ZygoteRules/M4xmc/src/adjoint.jl:72 [inlined]
 [44] ComposedFunction
    @ ./operators.jl:1041 [inlined]
 [45] (::Zygote.Pullback{Tuple{…}, Tuple{…}})(Δ::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface2.jl:0
 [46] (::Zygote.var"#75#76"{Zygote.Pullback{Tuple{…}, Tuple{…}}})(Δ::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/compiler/interface.jl:91
 [47] withjacobian(f::Function, args::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/grad.jl:150
 [48] jacobian(f::Function, args::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/jxHJc/src/lib/grad.jl:128
 [49] top-level scope
    @ ~/.julia/dev/decapodes-issue/example.jl:144
 [50] include(fname::String)
    @ Base.MainInclude ./client.jl:489
in expression starting at /home/x/.julia/dev/decapodes-issue/example.jl:144
Some type information was truncated. Use `show(err)` to see complete types.

@ChrisRackauckas
Copy link
Contributor

what's the function you're generating now?

@LilithHafner
Copy link

quote
    #= /home/x/.julia/dev/Decapodes/src/simulation.jl:594 =#
    function simulate(mesh, operators, hodge = GeometricHodge())
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:594 =#
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:595 =#
        begin
            #= /home/x/.julia/dev/Decapodes/src/simulation.jl:185 =#
            (var"GenSim-M_d₀", d₀) = default_dec_matrix_generate(mesh, :d₀, hodge)
            (var"GenSim-M_⋆₁", ₁) = default_dec_matrix_generate(mesh, :₁, hodge)
            (var"GenSim-M_dual_d₀", dual_d₀) = default_dec_matrix_generate(mesh, :dual_d₀, hodge)
            (var"GenSim-M_⋆₀⁻¹", ₀⁻¹) = default_dec_matrix_generate(mesh, :₀⁻¹, hodge)
            ♯ = operators(mesh, :♯)
            mag = operators(mesh, :mag)
            avg₀₁ = operators(mesh, :avg₀₁)
            (^) = operators(mesh, :^)
        end
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:596 =#
        begin
            #= /home/x/.julia/dev/Decapodes/src/simulation.jl:477 =#
            var"GenSim-M_GenSim-ConMat_1" = var"GenSim-M_⋆₀⁻¹" * var"GenSim-M_dual_d₀" * var"GenSim-M_⋆₁"
            var"GenSim-ConMat_1" = (x->var"GenSim-M_GenSim-ConMat_1" * x)
        end
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:597 =#
        begin
            #= /home/x/.julia/dev/Decapodes/src/simulation.jl:231 =#
            var"__dynamics_•1" = Decapodes.LazyBufferCache(Returns(nparts(mesh, :E)))
            var"__dynamics_•6" = Decapodes.LazyBufferCache(Returns(nparts(mesh, :E)))
            __dynamics_sum_1 = Decapodes.LazyBufferCache(Returns(nparts(mesh, :V)))
            __dynamics_mult_3 = Decapodes.LazyBufferCache(Returns(nparts(mesh, :E)))
            __dynamics_ḣ = Decapodes.LazyBufferCache(Returns(nparts(mesh, :V)))
        end
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:598 =#
        f(du, u, p, t) = begin
                #= /home/x/.julia/dev/Decapodes/src/simulation.jl:598 =#
                #= /home/x/.julia/dev/Decapodes/src/simulation.jl:599 =#
                begin
                    #= /home/x/.julia/dev/Decapodes/src/simulation.jl:256 =#
                    dynamics_h = u.dynamics_h
                    n = p.n
                    stress_A = p.stress_A
                    stress_ρ = p.stress_ρ
                    stress_g = p.stress_g
                    var"1" = 1.0
                    var"2" = 2.0
                    var"2" = 2.0
                end
                #= /home/x/.julia/dev/Decapodes/src/simulation.jl:600 =#
                var"dynamics_•1" = Decapodes.get_tmp(var"__dynamics_•1", u)
                var"dynamics_•6" = Decapodes.get_tmp(var"__dynamics_•6", u)
                dynamics_sum_1 = Decapodes.get_tmp(__dynamics_sum_1, u)
                dynamics_mult_3 = Decapodes.get_tmp(__dynamics_mult_3, u)
                dynamics_ḣ = Decapodes.get_tmp(__dynamics_ḣ, u)
                mul!(var"dynamics_•1", var"GenSim-M_d₀", dynamics_h)
                mul!(var"dynamics_•6", var"GenSim-M_d₀", dynamics_h)
                var"dynamics_•5" = (var"dynamics_•6")
                var"dynamics_•4" = mag(var"dynamics_•5")
                var"dynamics_•7" = n .- var"1"
                var"dynamics_•3" = var"dynamics_•4" ^ var"dynamics_•7"
                var"stress_•3" = stress_ρ .* stress_g
                var"stress_•2" = var"stress_•3" ^ n
                dynamics_sum_1 .= (.+)(n, var"2")
                stress_sum_1 = (.+)(n, var"2")
                var"dynamics_•2" = avg₀₁(var"dynamics_•3")
                var"dynamics_•9" = dynamics_h ^ dynamics_sum_1
                var"stress_•1" = var"2" / stress_sum_1
                stress_mult_1 = var"stress_•1" .* stress_A
                Γ = stress_mult_1 .* var"stress_•2"
                var"dynamics_•8" = avg₀₁(var"dynamics_•9")
                dynamics_mult_3 .= Γ .* var"dynamics_•1"
                dynamics_mult_1 = dynamics_mult_3 .* var"dynamics_•2"
                dynamics_mult_2 = dynamics_mult_1 .* var"dynamics_•8"
                mul!(dynamics_ḣ, var"GenSim-M_GenSim-ConMat_1", dynamics_mult_2)
                #= /home/x/.julia/dev/Decapodes/src/simulation.jl:601 =#
                getproperty(du, :dynamics_h) .= dynamics_ḣ
            end
    end
end

@LilithHafner
Copy link

LilithHafner commented Feb 11, 2024

And changing that generated function to

quote
    #= /home/x/.julia/dev/Decapodes/src/simulation.jl:594 =#
    function simulate(mesh, operators, hodge = GeometricHodge())
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:594 =#
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:595 =#
        begin
            #= /home/x/.julia/dev/Decapodes/src/simulation.jl:185 =#
            (var"GenSim-M_d₀", d₀) = default_dec_matrix_generate(mesh, :d₀, hodge)
            (var"GenSim-M_⋆₁", ₁) = default_dec_matrix_generate(mesh, :₁, hodge)
            (var"GenSim-M_dual_d₀", dual_d₀) = default_dec_matrix_generate(mesh, :dual_d₀, hodge)
            (var"GenSim-M_⋆₀⁻¹", ₀⁻¹) = default_dec_matrix_generate(mesh, :₀⁻¹, hodge)
            ♯ = operators(mesh, :♯)
            mag = operators(mesh, :mag)
            avg₀₁ = operators(mesh, :avg₀₁)
            (^) = operators(mesh, :^)
        end
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:596 =#
        begin
            #= /home/x/.julia/dev/Decapodes/src/simulation.jl:477 =#
            var"GenSim-M_GenSim-ConMat_1" = var"GenSim-M_⋆₀⁻¹" * var"GenSim-M_dual_d₀" * var"GenSim-M_⋆₁"
            var"GenSim-ConMat_1" = (x->var"GenSim-M_GenSim-ConMat_1" * x)
        end
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:597 =#
        begin
            #= /home/x/.julia/dev/Decapodes/src/simulation.jl:231 =#
            var"__dynamics_•1" = Decapodes.LazyBufferCache(Returns(nparts(mesh, :E)))
            var"__dynamics_•6" = Decapodes.LazyBufferCache(Returns(nparts(mesh, :E)))
            __dynamics_sum_1 = Decapodes.LazyBufferCache(Returns(nparts(mesh, :V)))
            __dynamics_mult_3 = Decapodes.LazyBufferCache(Returns(nparts(mesh, :E)))
            __dynamics_ḣ = Decapodes.LazyBufferCache(Returns(nparts(mesh, :V)))
        end
        #= /home/x/.julia/dev/Decapodes/src/simulation.jl:598 =#
        f(u, p, t) = begin
                #= /home/x/.julia/dev/Decapodes/src/simulation.jl:598 =#
                #= /home/x/.julia/dev/Decapodes/src/simulation.jl:599 =#
                begin
                    #= /home/x/.julia/dev/Decapodes/src/simulation.jl:256 =#
                    dynamics_h = u.dynamics_h
                    n = p.n
                    stress_A = p.stress_A
                    stress_ρ = p.stress_ρ
                    stress_g = p.stress_g
                    var"1" = 1.0
                    var"2" = 2.0
                    var"2" = 2.0
                end
                #= /home/x/.julia/dev/Decapodes/src/simulation.jl:600 =#
                # var"dynamics_•1" = Decapodes.get_tmp(var"__dynamics_•1", u)
                # var"dynamics_•6" = Decapodes.get_tmp(var"__dynamics_•6", u)
                # dynamics_sum_1 = Decapodes.get_tmp(__dynamics_sum_1, u)
                # dynamics_mult_3 = Decapodes.get_tmp(__dynamics_mult_3, u)
                dynamics_ḣ = Decapodes.get_tmp(__dynamics_ḣ, u)
                var"dynamics_•1" = var"GenSim-M_d₀" * dynamics_h
                var"dynamics_•6" = var"GenSim-M_d₀" * dynamics_h
                var"dynamics_•5" = (var"dynamics_•6")
                var"dynamics_•4" = mag(var"dynamics_•5")
                var"dynamics_•7" = n .- var"1"
                var"dynamics_•3" = var"dynamics_•4" ^ var"dynamics_•7"
                var"stress_•3" = stress_ρ .* stress_g
                var"stress_•2" = var"stress_•3" ^ n
                dynamics_sum_1 = (.+)(n, var"2")
                stress_sum_1 = (.+)(n, var"2")
                var"dynamics_•2" = avg₀₁(var"dynamics_•3")
                var"dynamics_•9" = dynamics_h ^ dynamics_sum_1
                var"stress_•1" = var"2" / stress_sum_1
                stress_mult_1 = var"stress_•1" .* stress_A
                Γ = stress_mult_1 .* var"stress_•2"
                var"dynamics_•8" = avg₀₁(var"dynamics_•9")
                dynamics_mult_3 = Γ .* var"dynamics_•1"
                dynamics_mult_1 = dynamics_mult_3 .* var"dynamics_•2"
                dynamics_mult_2 = dynamics_mult_1 .* var"dynamics_•8"
                dynamics_ḣ = var"GenSim-M_GenSim-ConMat_1" * dynamics_mult_2
                #= /home/x/.julia/dev/Decapodes/src/simulation.jl:601 =#
                ComponentArray(dynamics_h = dynamics_ḣ)
            end
    end
end

Results in the same errors and warnings, so I expect there is also mutation elsewhere.

@ChrisRackauckas
Copy link
Contributor

Mark all of the setup as non-differentiable

@LilithHafner
Copy link

I can't find the documentation of Zygote.@Nograd, if that is what you are referring to.

When I apply Zygote.@nograd to the simmulate function, I get the error on the forward pass, which is likely because I'm misusing this macro.

ERROR: UndefVarError: `GenSim-M_d₀` not defined
Stacktrace:
  [1] (::var"#f#41"{…})(u::ComponentVector{…}, p::ComponentVector{…}, t::Float64)
    @ Main ./REPL[5]:55
  [2] (::ODEFunction{…})(::ComponentVector{…}, ::Vararg{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/Avnpi/src/scimlfunctions.jl:2180
  [3] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.Tsit5ConstantCache)
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/perform_step/low_order_rk_perform_step.jl:726
  [4] __init(prob::ODEProblem{…}, alg::Tsit5{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:512
  [5] __init (repeats 5 times)
    @ ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:10 [inlined]
  [6] __solve(::ODEProblem{…}, ::Tsit5{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:5
  [7] __solve
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:1 [inlined]
  [8] #solve_call#34
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:609 [inlined]
  [9] solve_call
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:567 [inlined]
 [10] #solve_up#42
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:1058 [inlined]
 [11] solve_up
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:1044 [inlined]
 [12] #solve#40
    @ DiffEqBase ~/.julia/packages/DiffEqBase/UoaYd/src/solve.jl:981 [inlined]
 [13] f(flow_rate::Float64, ice_density::Float64, u_init_arr::Vector{Float64})
    @ Main ./REPL[8]:18
 [14] top-level scope
    @ REPL[10]:1

@ChrisRackauckas
Copy link
Contributor

No, it should be via ChainRules https://juliadiff.org/ChainRulesCore.jl/stable/rule_author/tips_for_packages.html

@LilithHafner
Copy link

Thanks! That makes more sense. Unfortunately, it seems to have no effect. Using this as the result of gensim gives the same warnings and errors as before:

sim = eval(quote
           #= /home/x/.julia/dev/Decapodes/src/simulation.jl:593 =#
           function simulate(mesh, operators, hodge = GeometricHodge())
               #= /home/x/.julia/dev/Decapodes/src/simulation.jl:593 =#
               #= /home/x/.julia/dev/Decapodes/src/simulation.jl:594 =#
               var"GenSim-M_d₀", var"GenSim-M_GenSim-ConMat_1", ♯, mag, avg₀₁, ^ = ChainRulesCore.ignore_derivatives() do
                   #= /home/x/.julia/dev/Decapodes/src/simulation.jl:185 =#
                   (var"GenSim-M_d₀", d₀) = default_dec_matrix_generate(mesh, :d₀, hodge)
                   (var"GenSim-M_⋆₁", ₁) = default_dec_matrix_generate(mesh, :₁, hodge)
                   (var"GenSim-M_dual_d₀", dual_d₀) = default_dec_matrix_generate(mesh, :dual_d₀, hodge)
                   (var"GenSim-M_⋆₀⁻¹", ₀⁻¹) = default_dec_matrix_generate(mesh, :₀⁻¹, hodge)
                   ♯ = operators(mesh, :♯)
                   mag = operators(mesh, :mag)
                   avg₀₁ = operators(mesh, :avg₀₁)
                   (^) = operators(mesh, :^)

               #= /home/x/.julia/dev/Decapodes/src/simulation.jl:595 =#

                   #= /home/x/.julia/dev/Decapodes/src/simulation.jl:477 =#
                   var"GenSim-M_GenSim-ConMat_1" = var"GenSim-M_⋆₀⁻¹" * var"GenSim-M_dual_d₀" * var"GenSim-M_⋆₁"
                   var"GenSim-ConMat_1" = (x->var"GenSim-M_GenSim-ConMat_1" * x)

               #= /home/x/.julia/dev/Decapodes/src/simulation.jl:596 =#

                   #= /home/x/.julia/dev/Decapodes/src/simulation.jl:231 =#
                   var"__dynamics_•1" = Decapodes.LazyBufferCache(Returns(nparts(mesh, :E)))
                   var"__dynamics_•6" = Decapodes.LazyBufferCache(Returns(nparts(mesh, :E)))
                   __dynamics_sum_1 = Decapodes.LazyBufferCache(Returns(nparts(mesh, :V)))
                   __dynamics_mult_3 = Decapodes.LazyBufferCache(Returns(nparts(mesh, :E)))
                   __dynamics_ḣ = Decapodes.LazyBufferCache(Returns(nparts(mesh, :V)))

                   var"GenSim-M_d₀", var"GenSim-M_GenSim-ConMat_1", ♯, mag, avg₀₁, ^
               end
               #= /home/x/.julia/dev/Decapodes/src/simulation.jl:597 =#
               f(u, p, t) = begin
                       #= /home/x/.julia/dev/Decapodes/src/simulation.jl:597 =#
                       #= /home/x/.julia/dev/Decapodes/src/simulation.jl:598 =#
                       begin
                           #= /home/x/.julia/dev/Decapodes/src/simulation.jl:256 =#
                           dynamics_h = u.dynamics_h
                           n = p.n
                           stress_A = p.stress_A
                           stress_ρ = p.stress_ρ
                           stress_g = p.stress_g
                           var"1" = 1.0
                           var"2" = 2.0
                           var"2" = 2.0
                       end
                       #= /home/x/.julia/dev/Decapodes/src/simulation.jl:599 =#
                       # var"dynamics_•1" = Decapodes.get_tmp(var"__dynamics_•1", u)
                       # var"dynamics_•6" = Decapodes.get_tmp(var"__dynamics_•6", u)
                       # dynamics_sum_1 = Decapodes.get_tmp(__dynamics_sum_1, u)
                       # dynamics_mult_3 = Decapodes.get_tmp(__dynamics_mult_3, u)
                       # dynamics_ḣ = Decapodes.get_tmp(__dynamics_ḣ, u)

                       # mul!(var"dynamics_•1", var"GenSim-M_d₀", dynamics_h)
                       var"dynamics_•1" = var"GenSim-M_d₀" * dynamics_h
                       # mul!(var"dynamics_•6", var"GenSim-M_d₀", dynamics_h)
                       var"dynamics_•6" = var"GenSim-M_d₀" * dynamics_h

                       var"dynamics_•5" = (var"dynamics_•6")
                       var"dynamics_•4" = mag(var"dynamics_•5")
                       var"dynamics_•7" = n .- var"1"
                       var"dynamics_•3" = var"dynamics_•4" ^ var"dynamics_•7"
                       var"stress_•3" = stress_ρ .* stress_g
                       var"stress_•2" = var"stress_•3" ^ n
                       # dynamics_sum_1 .= (.+)(n, var"2")
                       dynamics_sum_1 = (.+)(n, var"2")
                       stress_sum_1 = (.+)(n, var"2")
                       var"dynamics_•2" = avg₀₁(var"dynamics_•3")
                       var"dynamics_•9" = dynamics_h ^ dynamics_sum_1
                       var"stress_•1" = var"2" / stress_sum_1
                       stress_mult_1 = var"stress_•1" .* stress_A
                       Γ = stress_mult_1 .* var"stress_•2"
                       var"dynamics_•8" = avg₀₁(var"dynamics_•9")
                       # dynamics_mult_3 .= Γ .* var"dynamics_•1"
                       dynamics_mult_3 = Γ .* var"dynamics_•1"
                       dynamics_mult_1 = dynamics_mult_3 .* var"dynamics_•2"
                       dynamics_mult_2 = dynamics_mult_1 .* var"dynamics_•8"
                       # mul!(dynamics_ḣ, var"GenSim-M_GenSim-ConMat_1", dynamics_mult_2)
                       dynamics_ḣ = var"GenSim-M_GenSim-ConMat_1" * dynamics_mult_2
                       #= /home/x/.julia/dev/Decapodes/src/simulation.jl:600 =#
                       # getproperty(du, :dynamics_h) .= dynamics_ḣ
                       ComponentArray(dynamics_h=dynamics_ḣ)
                   end
           end
           simulate
       end)

Notably including "Warning: Reverse-Mode AD VJP choices all failed. Falling back to numerical VJPs"

@ChrisRackauckas
Copy link
Contributor

To work around the mutation issues, we're going to get enough Enzyme rules to Enzyme everywhere. The ODE solver mostly has this done, and we're working on doing a bunch of unit tests SciML/SciMLSensitivity.jl#905 and fixing what we find there before doing it on this example. I assume we can get something of that merging either tonight or tomorrow, then try changing this to Enzyme.

CC @ArnoStrouwen

@ChrisRackauckas
Copy link
Contributor

We have a gradient!

cd(@__DIR__)
using Pkg
Pkg.activate(".")
Pkg.instantiate()

# AlgebraicJulia Dependencies
using Catlab
using Catlab.Graphics
using CombinatorialSpaces
using Decapodes
using ComponentArrays

# External Dependencies
using MLStyle
using MultiScaleArrays
using LinearAlgebra
using OrdinaryDiffEq
using JLD2
using SparseArrays
using Statistics
using GLMakie # Just for visualization
using GeometryBasics: Point2, Point3
Point2D = Point2{Float64};
Point3D = Point3{Float64};

using DiagrammaticEquations
using DiagrammaticEquations.Deca

@info("Packages Loaded")

halfar_eq2 = @decapode begin
  h::Form0
  Γ::Form1
  n::Constant== ∂ₜ(h)
  ḣ == (, d, )(Γ * d(h) * avg₀₁(mag((d(h)))^(n-1)) * avg₀₁(h^(n+2)))
end
glens_law = @decapode begin
  Γ::Form1
  (A,ρ,g,n)::Constant

  Γ == (2/(n+2))*A**g)^n
end

@info("Decapodes Defined")

ice_dynamics_composition_diagram = @relation () begin
  dynamics(Γ,n)
  stress(Γ,n)
end
ice_dynamics_cospan = oapply(ice_dynamics_composition_diagram,
  [Open(halfar_eq2, [,:n]),
  Open(glens_law, [,:n])])
ice_dynamics = apex(ice_dynamics_cospan)
ice_dynamics1D = expand_operators(ice_dynamics)
infer_types!(ice_dynamics1D, op1_inf_rules_1D, op2_inf_rules_1D)
resolve_overloads!(ice_dynamics1D, op1_res_rules_1D, op2_res_rules_1D)

s_prime = EmbeddedDeltaSet1D{Bool, Point2D}()
add_vertices!(s_prime, 100, point=Point2D.(range(-2, 2, length=100), 0))
add_edges!(s_prime, 1:nv(s_prime)-1, 2:nv(s_prime))
orient!(s_prime)
s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s_prime)
subdivide_duals!(s, Circumcenter())

@info("Spaces Defined")

function generate(sd, my_symbol; hodge=GeometricHodge())
  e_vecs = map(edges(sd)) do e
    point(sd, sd[e, :∂v0]) - point(sd, sd[e, :∂v1])
  end
  neighbors = map(vertices(sd)) do v
    union(incident(sd, v, :∂v0), incident(sd, v, :∂v1))
  end
  n_vecs = map(neighbors) do es
    [e_vecs[e] for e in es]
  end
  I = Vector{Int64}()
  J = Vector{Int64}()
  V = Vector{Float64}()
  for e in 1:ne(s)
      append!(J, [s[e,:∂v0],s[e,:∂v1]])
      append!(I, [e,e])
      append!(V, [0.5, 0.5])
  end
  avg_mat = sparse(I,J,V)
  op = @match my_symbol begin
    :♯ => x -> begin
      map(neighbors, n_vecs) do es, nvs
        sum([nv*norm(nv)*x[e] for (e,nv) in zip(es,nvs)]) / sum(norm.(nvs))
      end
    end
    :mag => x -> begin
      norm.(x)
    end
    :avg₀₁ => x -> begin
      avg_mat * x
    end
    :^ => (x,y) -> max.(x, 0) .^ y # x is sometimes a very small negative number due to numerical error
    :* => (x,y) -> x .* y
    :show => x -> begin
      @show x
      x
    end
    x => error("Unmatched operator $my_symbol")
  end
  return (args...) -> op(args...)
end

sim = eval(gensim(ice_dynamics1D, dimension=1))

fₘ = sim(s, generate)

function f(constants_and_parameters)
  prob = ODEProblem{true, SciMLBase.FullSpecialize}(fₘ, u₀, (0, tₑ), constants_and_parameters)
  @info("Solving")
  soln = solve(prob, Tsit5())
  @info("Done")

  # return soln(tₑ)
  sum(last(soln)) # last, not soln(tₑ) because to avoid interpolation fails when AD fails.
end

h₀ = map(x -> try sqrt(1. - x[1]^2) catch DomainError return 0.0 end, point(s_prime))

flow_rate, ice_density, u_init_arr = 1e-16, 910., h₀
n = 3
ρ = ice_density
g = 9.8101
A = fill(flow_rate, ne(s))
tₑ = 5e13

u₀ = ComponentArray(dynamics_h = u_init_arr)

# Note that this must be a ComponentArray to differentiate
constants_and_parameters = ComponentArray(
  n = n,
  stress_ρ = ρ,
  stress_g = g,
  stress_A = A)

y = f(constants_and_parameters)

# Must use Zygote—ForwardDiff is asymtotically slow
using SciMLSensitivity
using Zygote

Zygote.gradient(f, constants_and_parameters)

Still more to do, it's using some fallbacks, but this is at least a starting point.

@ChrisRackauckas
Copy link
Contributor

DecapodesDiff.zip

Hand modified decapode:

decapode_f = begin
    #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:531 =#
    function simulate(mesh, operators, hodge = GeometricHodge())
        #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:531 =#
        #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:532 =#
        begin
            #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:157 =#
            (var"GenSim-M_d₀", d₀) = default_dec_matrix_generate(mesh, :d₀, hodge)
            (var"GenSim-M_⋆₁", ₁) = default_dec_matrix_generate(mesh, :₁, hodge)
            (var"GenSim-M_dual_d₀", dual_d₀) = default_dec_matrix_generate(mesh, :dual_d₀, hodge)
            (var"GenSim-M_⋆₀⁻¹", ₀⁻¹) = default_dec_matrix_generate(mesh, :₀⁻¹, hodge)
            ♯ = operators(mesh, :♯)
            mag = operators(mesh, :mag)
            avg₀₁ = operators(mesh, :avg₀₁)
            (^) = operators(mesh, :^)
        end
        #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:533 =#
        begin
            #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:449 =#
            var"GenSim-M_GenSim-ConMat_1" = var"GenSim-M_⋆₀⁻¹" * var"GenSim-M_dual_d₀" * var"GenSim-M_⋆₁"
            var"GenSim-ConMat_1" = (x->var"GenSim-M_GenSim-ConMat_1" * x)
        end
        #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:534 =#
        begin
            #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:203 =#
            var"__dynamics_•1" = Decapodes.LazyBufferCache(x->nparts(mesh, :E))
            var"__dynamics_•6" = Decapodes.LazyBufferCache(x->nparts(mesh, :E))
            __dynamics_sum_1 = Decapodes.LazyBufferCache(x->nparts(mesh, :V))
            __dynamics_mult_3 = Decapodes.LazyBufferCache(x->nparts(mesh, :E))
            __dynamics_ḣ = Decapodes.LazyBufferCache(x->nparts(mesh, :V))
        end
        #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:535 =#
        f(du, u, p, t) = begin
                #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:535 =#
                #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:536 =#
                begin
                    #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:228 =#
                    dynamics_h = u.dynamics_h
                    n = p.n
                    stress_A = p.stress_A
                    stress_ρ = p.stress_ρ
                    stress_g = p.stress_g
                    var"1" = 1.0
                    var"2" = 2.0
                    var"2" = 2.0
                end
                #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:537 =#
                var"dynamics_•1" = Decapodes.get_tmp(var"__dynamics_•1", u)
                var"dynamics_•6" = Decapodes.get_tmp(var"__dynamics_•6", u)
                dynamics_sum_1 = Decapodes.get_tmp(__dynamics_sum_1, u)
                dynamics_mult_3 = Decapodes.get_tmp(__dynamics_mult_3, u)
                dynamics_ḣ = Decapodes.get_tmp(__dynamics_ḣ, u)
                mul!(var"dynamics_•1", var"GenSim-M_d₀", dynamics_h)
                mul!(var"dynamics_•6", var"GenSim-M_d₀", dynamics_h)
                var"dynamics_•5" = (var"dynamics_•6")
                var"dynamics_•4" = mag(var"dynamics_•5")
                var"dynamics_•7" = n .- var"1"
                var"dynamics_•3" = var"dynamics_•4" ^ var"dynamics_•7"
                var"stress_•3" = stress_ρ .* stress_g
                var"stress_•2" = var"stress_•3" ^ n
                dynamics_sum_1 .= (.+)(n, var"2")
                stress_sum_1 = (.+)(n, var"2")
                var"dynamics_•2" = avg₀₁(var"dynamics_•3")
                @show typeof(dynamics_h)
                @show typeof(dynamics_sum_1)
                var"dynamics_•9" = dynamics_h ^ dynamics_sum_1
                var"stress_•1" = var"2" / stress_sum_1
                stress_mult_1 = var"stress_•1" .* stress_A
                Γ = stress_mult_1 .* var"stress_•2"
                var"dynamics_•8" = avg₀₁(var"dynamics_•9")
                dynamics_mult_3 .= Γ .* var"dynamics_•1"
                dynamics_mult_1 = dynamics_mult_3 .* var"dynamics_•2"
                dynamics_mult_2 = dynamics_mult_1 .* var"dynamics_•8"
                mul!(dynamics_ḣ, var"GenSim-M_GenSim-ConMat_1", dynamics_mult_2)
                @show size(du)
                @show size(dynamics_ḣ)
                #= /Users/chrisrackauckas/.julia/packages/Decapodes/GugTA/src/simulation.jl:538 =#
                du .= dynamics_ḣ
            end
    end
end

This is starting to do stuff but I think the AD there is doing something odd, so I'm debugging the lazy buffer cache handling.

@djinnome
Copy link
Author

djinnome commented Mar 4, 2024

Hi @ChrisRackauckas @LilithHafner ,

What is the current status of this issue?

I tried running the test.jl code by unzipping DecapodesDiff.zip above, but I got this error message:

julia> include("test.jl")
  Activating project at `~/Projects/Proposals/ASKEM/DecapodesDiff`
[ Info: Packages Loaded
[ Info: Decapodes Defined
[ Info: Spaces Defined
ERROR: LoadError: UndefVarError: `sim` not defined
Stacktrace:
 [1] top-level scope
   @ ~/Projects/Proposals/ASKEM/DecapodesDiff/test.jl:117
 [2] include(fname::String)
   @ Base.MainInclude ./client.jl:489
 [3] top-level scope
   @ REPL[13]:1
in expression starting at /Users/zuck016/Projects/Proposals/ASKEM/DecapodesDiff/test.jl:117

Presumably, this is due to the fact that the call to sim on line 117 is not defined. Is the code missing an appropriate using or include statement? If not, which package should have contained sim?

fₘ = sim(s, generate)

@LilithHafner
Copy link

I've handed this off to @ChrisRackauckas. I'm not sure what the current status is.

@ChrisRackauckas
Copy link
Contributor

Yeah oops I accidentally included some random changes in there. This one should just run.

DecapodesDiff.zip

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
Status: No status
Development

No branches or pull requests

8 participants