Open djinnome opened 11 months ago
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)
I would try out Barycenter() instead of Circumcenter() as a quick check
Barycenter gives the same error. :(
@ChrisRackauckas any progress to report on this issue?
Any updates on this issue?
I'm finding someone to investigate this.
I don't suppose there's a more minimal reproducible example?
Here's a slightly modified reproducer that throws an Internal error: encountered unexpected error in runtime
and segfaults.
Invocation:
julia --project --startup-file=no -e 'using Revise; includet("segfault_reproducer.jl")'
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
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.
@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?
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 https://github.com/SciML/PreallocationTools.jl/pull/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.
With https://github.com/SciML/PreallocationTools.jl/pull/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 Vector
s as output when given CategoricalVector
s 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.
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?
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)
That was after moving the function construction into the loss?
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))))))
So it looks like we can apply that patch to main Decapodes and this will work?
When I applied that patch locally this example hung for me.
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.
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.
@LilithHafner These ought to be exported when using using DiagrammaticEquations
and using DiagrammaticEquations.Deca
. Does this work on your environment?
Yes. The OP with
using DiagrammaticEquations
using DiagrammaticEquations.Deca
added also reproduces for me.
Here are some solutions I've tried which did not work:
FixedSizeDiffCache gives the error in the OP
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.
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.
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]
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.dynamicsh 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.
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?
Do you think that if we removed caching entirely we would see the increased calls to the dynamics?
Im trying to understand if the increased number of calls to the dynamics is because of caching or because of the backwards pass autodiff inherently.
We could remove the caching and just have a less memory efficient but simpler runtime if that helps.
@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.
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
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?
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.
Which vjp method is chosen for the reverse mode?
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
Oh if it's ForwardDiff then it's just forward mode and going to be slow. That should be changed to a Zygote call.
The time for forward mode is O((np)), so the scaling will be extremely bad with that ( O((np)^3) with an implicit method)
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.
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?
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.
Change
soln = solve(prob, Tsit5())
@info("Done")
return soln(tₑ)
to
soln = solve(prob, Tsit5(), saveat = tₑ)
return Array(soln)
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.
what's the function you're generating now?
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
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.
Mark all of the setup as non-differentiable
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
No, it should be via ChainRules https://juliadiff.org/ChainRulesCore.jl/stable/rule_author/tips_for_packages.html
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"
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 https://github.com/SciML/SciMLSensitivity.jl/pull/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
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.
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.
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)
I've handed this off to @ChrisRackauckas. I'm not sure what the current status is.
Yeah oops I accidentally included some random changes in there. This one should just run.
<@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