search

SciML / Catalyst.jl

Chemical reaction network and systems biology interface for scientific machine learning (SciML). High performance, GPU-parallelized, and O(1) solvers in open source software.
https://docs.sciml.ai/Catalyst/stable/
Other
442 stars 71 forks source link

Variable rate constants in programatically constructed ReactionSystems #586

Closed joegilkes closed 1 year ago

joegilkes commented 1 year ago

I'm trying to implement temperature-dependent rate constants over a large, programatically constructed reaction network using the ReactionSystem API. One way I thought this might be possible was by implementing the change in temperature, and therefore the temperature dependence of the rates, as a constraint ODESystem coupled to the network. I wanted to avoid defining a function for each rate constant because the networks I'm working on typically have upwards of 3000 reactions, but if there's a way to make this work with a single function let me know. The programmatic construction made it a bit of a challenge, but I believe I'm 90% of the way there with the following MWE:

using Catalyst
using OrdinaryDiffEq
using IfElse

@variables t
@variables T(t) [isbcspecies=true]
@parameters hr cr
@variables (species(t))[1:2]
@variables (k(t))[1:2] [isbcspecies=true, irreducible=true]

D = Differential(t)
function Tgrad(t)
    return IfElse.ifelse(t > 15.0, cr, 
        IfElse.ifelse(t < 5.0, hr, 0.0)
    )
end
function calc_rates(T)
    base_rates = [0.5, 0.5]
    # Only the reverse reaction is temperature dependent for the sake of this example
    return [base_rates[1], base_rates[2] * T]
end

@named rate_model = ODESystem([
    D(T) ~ Tgrad(t);
    k .~ calc_rates(T)
], t)
rate_model_simplified = structural_simplify(rate_model)

rxs = [
    Reaction(k[1], [species[1]], [species[2]], [1], [1]),
    Reaction(k[2], [species[2]], [species[1]], [1], [1])
]
@named rs = ReactionSystem(rxs, t, collect(species), collect(k); checks=false, constraints=rate_model_simplified)

T0 = 0.0
u0 = [10.0, 0.0]
u0map = Pair.(collect(species), u0)
kmap = Pair.(collect(k), calc_rates(T0))
u0map = vcat(u0map, kmap)
push!(u0map, Pair(T, T0))

p = [1.0, -1.0]    # heating rate `hr` and cooling rate `cr` respectively
pmap = Pair.([hr, cr], p)

tspan = (0.0, 20.0)

oprob = ODEProblem(rs, u0map, tspan, pmap)
sol = solve(oprob, QNDF())

The above defines a simple 2-species system with a single reversible reaction connecting them. The intended result is that k[1] starts at 0.5 while k[2] starts at 0 (so species[1] gets converted to species[2]), and then as the simulation progresses the temperature rises (with associated heating rate hr = 1.0), causing k[2] to increase (so species[1] dominates]. At 5 seconds in, the temperature is supposed to plateau, and then at 15 seconds it is supposed to decrease (with associated cooling rate cr = -1.0 (so species[1] gets converted back to species[2]).

Separating out rate_model_simplified and solving it by itself works fine (note that irreducible=true has to be included for k(t) because otherwise they are marked as algebraic variables and are eliminated from the ODESystem upon structural simplification). However, running the whole thing causes the following error:

ERROR: LinearAlgebra.SingularException(4)
Stacktrace:

  [1] checknonsingular
    @ ~/.local/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:19 [inlined]
  [2] generic_lufact!(A::Matrix{Float64}, pivot::LinearAlgebra.RowMaximum; check::Bool)
    @ LinearAlgebra ~/.local/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/lu.jl:186
  [3] generic_lufact!
    @ ~/.local/julia-1.8.2/share/julia/stdlib/v1.8/LinearAlgebra/src/lu.jl:137 [inlined]
  [4] do_factorization
    @ ~/.julia/packages/LinearSolve/dxfUd/src/factorization.jl:65 [inlined]
  [5] #solve#6
    @ ~/.julia/packages/LinearSolve/dxfUd/src/factorization.jl:11 [inlined]
  [6] #solve#5
    @ ~/.julia/packages/LinearSolve/dxfUd/src/common.jl:161 [inlined]
  [7] #dolinsolve#3
    @ ~/.julia/packages/OrdinaryDiffEq/pIBDs/src/misc_utils.jl:109 [inlined]
  [8] compute_step!(nlsolver::OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, integrator::OrdinaryDiffEq.ODEIntegrator{QNDF{5, 5, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, QNDF{5, 5, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.DummyController, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.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{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/pIBDs/src/nlsolve/newton.jl:302
  [9] nlsolve!(nlsolver::OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, integrator::OrdinaryDiffEq.ODEIntegrator{QNDF{5, 5, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, QNDF{5, 5, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.DummyController, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.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{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/pIBDs/src/nlsolve/nlsolve.jl:48
 [10] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{QNDF{5, 5, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, QNDF{5, 5, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.DummyController, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.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{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/pIBDs/src/perform_step/bdf_perform_step.jl:904
 [11] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEq/pIBDs/src/perform_step/bdf_perform_step.jl:841 [inlined]
 [12] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{QNDF{5, 5, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, QNDF{5, 5, true, GenericLUFactorization{LinearAlgebra.RowMaximum}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 5}}, Vector{Float64}, Vector{Vector{NTuple{5, Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, GenericLUFactorization{LinearAlgebra.RowMaximum}, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.DummyController, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.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{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/pIBDs/src/solve.jl:520
 [13] #__solve#621
    @ ~/.julia/packages/OrdinaryDiffEq/pIBDs/src/solve.jl:6 [inlined]
 [14] #solve_call#22
    @ ~/.julia/packages/DiffEqBase/190F1/src/solve.jl:494 [inlined]
 [15] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::QNDF{5, 0, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}; kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:progress, :abstol, :reltol), Tuple{Bool, Float64, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/190F1/src/solve.jl:915
 [16] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#485"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x8a869998, 0xa29fdcc5, 0x196b4fe1, 0x66fa7b8e, 0x532f38e1)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xdd906b58, 0x3e364ce2, 0x6d731476, 0xf76c2a57, 0x58e62691)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#502#generated_observed#493"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::QNDF{5, 0, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:progress, :abstol, :reltol), Tuple{Bool, Float64, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/190F1/src/solve.jl:825

I'm at a loss with what could be causing this since it's quite deep into the solve. I've tried a variety of different solvers and they seem to have varying issues with the way I'm defining the temperature ODE, but QNDF should work fine. Is there anything I'm doing wrong in the problem definition as a whole? I appreciate this isn't exactly charted territory in terms of the documentation, just wondering if anyone has done something similar before and got it to work.

Alternatively, would it be better to just define the ReactionSystem as normal and introduce a callback on every timestep to remake() the ODEProblem with new, temperature dependent rate constants? I initially avoided this way since calling remake so many times on such a large network seemed like it would be slow.

TorkelE commented 1 year ago

I am shooting entirely from the hip here, but won't constrained systems (e.g. X <--> Y) result in singular Jacobians, which can result in the type of error produced here?

Avoiding the callback seems sensible, but this kind of MTK application is not my speciality.

isaacsas commented 1 year ago

I'd suggest waiting to run structural_simplify until you've formed your full, composed model. i.e. don't run it on rate_model, but instead do something like

@named rs = ReactionSystem(rxs, t, collect(species), collect(k); checks=false, constraints=rate_model)
osys = convert(ODESystem, rs)
osys_simplified = structural_simplify(osys)

That said, they seem equivalent in this case. However, if we instead use ODAEProblem the code seems to run

@named rs = ReactionSystem(rxs, t, collect(species), collect(k); checks=false, constraints=rate_model)
osys = structural_simplify(convert(ODESystem, rs))
oprob = ODAEProblem(osys, u0map, tspan, pmap)
sol = solve(oprob, Rodas5())   # this runs as does QNDF

Why are you marking the k(t)s irreducible? In this case it seems like it would be perfectly reasonable to substitute them out, and indeed if I let structural_simplify eliminate them we get a system of only ODEs which can be solved just fine with ODEProblem. With them being irreducible you are getting a system of DAEs, and so need to keep that in mind when thinking about solving the system.

If you know the explicit times where you will switch the value of Tgrad, and it is constant between switches, using a DiscreteCallback and adding tstops to solve should be a good approach (just make Tgrad a parameter and add the callback to the constraint system).

isaacsas commented 1 year ago

Sorry, just to clarify; I meant that you could make Tgrad a parameter and change it via a callback...

isaacsas commented 1 year ago

@TorkelE just the normal X <--> Y reaction shouldn't cause such a warning (the matrix the warning is referring to is not the same as the reaction system's Jacobian). It would be really bad if such a simple reaction system couldn't be solved by standard ODE solvers.

joegilkes commented 1 year ago

won't constrained systems (e.g. X <--> Y) result in singular Jacobians, which can result in the type of error produced here?

Very possible, I'm not very familiar with how DAEs mix with ODEs so this could very well be the problem.

if we instead use ODAEProblem the code seems to run

This combined with running structural_simplify after composing the full system also works for me, thanks @isaacsas! I'd seen ODAEProblem appear in a couple of other places but couldn't find any documentation so I'd given it a miss, but using it seems to work well.

Why are you marking the k(t)s irreducible?

This was due to the ordering of the structural_simplify. Because I was simplifying before composing the full system, it was eliminating the k(t)s so that the ReactionSystem was only considering them as parameters, not variables, and the rate constants were stuck at their initial values. With the simplification in the right place it's now solvable as a straight ODEProblem again yes, that's great thanks.

using a DiscreteCallback and adding tstops to solve should be a good approach

I agree, however I'm trying to keep it general for any temperature profile (looking at doing periodic temperature ramps etc. later on) so I think coupling the temperature ODE is still preferable.

joegilkes commented 1 year ago

Thanks both for the help, I'll try slotting this into my main network code tomorrow and see how it goes, but in theory it's just a scaled up version of the MWE so should be fine. I'll reopen if there are any other issues with this.

TorkelE commented 1 year ago

@TorkelE just the normal X <--> Y reaction shouldn't cause such a warning (the matrix the warning is referring to is not the same as the reaction system's Jacobian). It would be really bad if such a simple reaction system couldn't be solved by standard ODE solvers.

Agreed, but occasionally I run into these troubles when I use niche methods. Thinking closer though I think BifurcationKit might be the only culprit (which I probably should try to find a long-term solution to at some point...).

isaacsas commented 1 year ago

Cool. Let us know if you have any other issues! This is a nice example showing off some stuff we haven't really gotten into the tutorials yet, and which has had limited usage (so isn't super documented for users yet).

Catalyst 13 will hopefully make having / adding constraints cleaner, though it will break the interface a bit (your non-Reaction ODEs or algebraic equations should just be addable as normal equations when creating the ReactionSystem). Once 13 is ready we should consider a variant on this example as a tutorial / application in the docs.

TorkelE commented 1 year ago

Yes, having settled on constraints handling would make it possible to move forward with that.

isaacsas commented 1 year ago

BifurcationKit would give such warnings; it uses factorizations of the actual Jacobian. That is different than the ODE solvers. For your example(s) where you had issues did you try eliminating the conservation laws in Catalyst before converting to ODEs? That should have gotten rid of the singularity.

TorkelE commented 1 year ago

In the end I just added reactions like

    (xTot,1), 0 <--> X

where I knew xTot was the conserved concentration or something. Not pretty but was the fastest. But I agree that going forward we should make a tutorial how to do this according to how you suggest

isaacsas commented 1 year ago

@TorkelE, see https://docs.sciml.ai/Catalyst/stable/tutorials/homotopy_continuation_tutorial/#Systems-with-Conservation-Laws. Should be the same as for using homotopy continuation. Just pass the relevant keyword when calling convert to remove conserved quantities and get a non-singular Jacobian.

joegilkes commented 1 year ago

Unfortunately when scaling this up to one of the networks I'm working on (2492 reactions, 415 species) the ODEProblem seems to become unsolvable. I usually formulate it as a sparse problem to improve performance but the sparse linear solvers seem to have issues with the matrices being non-square (even QR solvers seem to have issues, and as far as I can tell they should work for non-square). Leaving sparsity out of it, the problem takes about 4x longer to be created and get to solving, and then the solver seems to really struggle to make progress. I haven't looked into where it's having issues yet, but after being left for an hour there was no appreciable progress.

I'm guessing there isn't an easy way to solve the non-square matrix issue? Failing that, I might try the full callback approach since I've got it working on the MWE and it seems to be fast still since it doesn't actually involve a full remake, and all the rate constants can be recalculated simultaneously at each timestep and reapplied to integrator.p, but I'm expecting a significant slowdown.

joegilkes commented 1 year ago

It turns out there was a unit conversion issue on my end that was causing numerical instability, this does actually solve when not a sparse problem! It is still very expensive to solve compared to an equivalent constant temperature run and requires considerable amounts of memory, so I'd prefer to get sparse solving working if at all possible.

Creating the ODEProblem with sparse=true and ensuring that I'm using a non-square matrix compatible linear solver (QRFactorization, which for a sparse problem should mean using SuiteSparse's SPQR) results in the following error:

ERROR: DimensionMismatch: array could not be broadcast to match destination
Stacktrace:
  [1] check_broadcast_shape
    @ ./broadcast.jl:540 [inlined]
  [2] check_broadcast_axes
    @ ./broadcast.jl:543 [inlined]
  [3] check_broadcast_axes
    @ ./broadcast.jl:546 [inlined]
  [4] instantiate
    @ ./broadcast.jl:284 [inlined]
  [5] materialize!
    @ ./broadcast.jl:871 [inlined]
  [6] materialize!
    @ ./broadcast.jl:868 [inlined]
  [7] forwarddiff_color_jacobian!(J::SparseArrays.SparseMatrixCSC{Float64, Int64}, f::SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, x::Vector{Float64}, jac_cache::SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing})
    @ SparseDiffTools ~/.julia/packages/SparseDiffTools/zGdIo/src/differentiation/compute_jacobian_ad.jl:437
  [8] jacobian!(J::SparseArrays.SparseMatrixCSC{Float64, Int64}, f::SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, x::Vector{Float64}, fx::Vector{Float64}, integrator::OrdinaryDiffEq.ODEIntegrator{QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.DummyController, 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{}, Vector{Float64}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, jac_config::SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/derivative_wrappers.jl:228
  [9] calc_J!
    @ ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/derivative_utils.jl:169 [inlined]
 [10] calc_W!(W::SparseArrays.SparseMatrixCSC{Float64, Int64}, integrator::OrdinaryDiffEq.ODEIntegrator{QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.DummyController, 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{}, Vector{Float64}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, nlsolver::OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, cache::OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, dtgamma::Float64, repeat_step::Bool, W_transform::Bool, newJW::Nothing)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/derivative_utils.jl:716
 [11] update_W!
    @ ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/derivative_utils.jl:824 [inlined]
 [12] update_W!
    @ ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/derivative_utils.jl:823 [inlined]
 [13] nlsolve!(nlsolver::OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, integrator::OrdinaryDiffEq.ODEIntegrator{QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.DummyController, 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{}, Vector{Float64}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/nlsolve/nlsolve.jl:25
 [14] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.DummyController, 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{}, Vector{Float64}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/perform_step/bdf_perform_step.jl:904
 [15] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/perform_step/bdf_perform_step.jl:841 [inlined]
 [16] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Float64, Vector{Float64}}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Vector{Tuple{Float64}}}, UnitRange{Int64}, Nothing}, LinearSolve.LinearCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, QRFactorization{NoPivot}, QR{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, Nothing}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.DummyController, 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{}, Vector{Float64}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/solve.jl:514
 [17] __solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::QNDF{5, 1, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}; kwargs::Base.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/tAXVW/src/solve.jl:6
 [18] #solve_call#21
    @ ~/.julia/packages/DiffEqBase/ENGnO/src/solve.jl:494 [inlined]
 [19] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf82da3ee, 0x0e1ee9ed, 0x3de93b5c, 0x27232f01, 0xb8548b14)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x4bcc00cb, 0xbdc85050, 0x9e96a1be, 0x9c3cc8e4, 0x1d00b6d4)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::QNDF{5, 0, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}; kwargs::Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/ENGnO/src/solve.jl:856
 [20] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf82da3ee, 0x0e1ee9ed, 0x3de93b5c, 0x27232f01, 0xb8548b14)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x4bcc00cb, 0xbdc85050, 0x9e96a1be, 0x9c3cc8e4, 0x1d00b6d4)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::QNDF{5, 0, true, QRFactorization{NoPivot}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/ENGnO/src/solve.jl:823

I'm not sure how to dig into the actual matrices that the linear solver is getting here, so I don't know how to diagnose the cause any further.

Alternatively, I've noticed that much of the time spent solving is spent in compilation, since running the same solve again after it's been compiled once is ~10x faster. I realise it's probably more of a question for the DiffEq crowd, but is there an easy way of doing the compilation on a smaller system and then scaling up the compiled routines for a larger system, or does it not work like that?

joegilkes commented 1 year ago

As for the callback method, it unfortunately (but expectedly) ruins performance to the point of not running, as it has to stop at every timestep to recalculate the temperature and rates, ruining the adaptive timestepping of the ODE solver and causing it to take the smallest timesteps possible for the whole calculation.

isaacsas commented 1 year ago

What is the non-square matrix you are referring to?

For your full model, does it simplify down after structural_simplify to just a set of ODEs, or are there still algebraic equations?

For a couple larger systems biology models we found that the best performance was when using preconditioned GMRES with sparse=true and jac=true in the ODEProblem, and then making sure to set autodiff = false in the solver, see

https://docs.sciml.ai/SciMLBenchmarksOutput/dev/Bio/BCR/

https://docs.sciml.ai/SciMLBenchmarksOutput/dev/Bio/fceri_gamma2/

isaacsas commented 1 year ago

Have you tried solving a case without any Jacobian or sparsity set, but also disabling autodiff as I mentioned above? Using autodiff can cause memory use issues for these larger models.

joegilkes commented 1 year ago

Hmm, so I'd been leaving jac=false while testing since in my experience, creating an analytic Jacobian roughly doubled compilation time while only reducing solve time by a small amount. Running with sparse=true, jac=true, linsolve=QRFactorization() and autodiff=false actually compiles in a similar amount of time to the jac=false, sparse=false case, and the solve does run but it slows down to a crawl (extrapolating it out, probably ~2 hours instead of ~2 minutes). Memory usage is much lower too.

For the jac=false, sparse=false, autodiff=false case, compilation is about the same, solving is about twice as fast, and it's still quite memory hungry.

I haven't fully got my head around setting up the preconditioner for preconditioned GMRES yet, but I'll definitely keep that in mind.

For your full model, does it simplify down after structural_simplify to just a set of ODEs, or are there still algebraic equations?

What would be the best way to examine this? I can do @show equations(osys) but the output is so vast that it's unclear.

isaacsas commented 1 year ago

Why are you using QR? You mentioned a non-square matrix before, are you saying you have a non-square Jacobian? Have you actually checked this? I don't think that should be happening. If it is square and sparse usually the recommendation is to use KLUFactorization.

For checking if you have diff eqs maybe you can do something like

all(ModelingToolkit.isdiffeq, equations(osys))

I haven't used that function before, but if your equations are in canonical form, which I think they should be after calling structural_simplify, I think it should be able to tell which are ODEs.

joegilkes commented 1 year ago

Also, is there a way to store the 'observable' rate constants which have been eliminated from the ODESystem? Currently, trying to retrieve them by doing sol[k] only returns the rate constants at what I assume is the final timestep since they get cached, but attempting to access any other timesteps (e.g. sol[k, 1] or sol[k, :]) results in a never ending computation as it tries to lazily reconstruct the requested values and has significant difficulty doing so. Oddly, the same behaviour is not true of the MWE system, as running sol[k] there retrieves both rate constants at all timesteps.

isaacsas commented 1 year ago

You could use a custom SavingCallback and just save them yourself. I think ModelingToolkit has been doing some work on that observable issue you mention, but I'm not sure where it is at (if you open an issue there they might be able to give you some guidance on how to do that).

isaacsas commented 1 year ago

What is the size of oprob.f.jac_prototype?

What is oprob.f.mass_matrix, is it just the identity or an actual matrix? If a matrix, what size? Is it sparse?

isaacsas commented 1 year ago

Also, would you be willing to share your model code with myself and @ChrisRackauckas for testing purposes (via a private repo would be fine). That would make it easier for us to play around with this, and if I conclude everything looks ok I can turn it over to Chris (who is much more of an ODE/DAE performance optimization expert than myself).

joegilkes commented 1 year ago

For your full model, does it simplify down after structural_simplify to just a set of ODEs, or are there still algebraic equations?

It's ODEs only.

Why are you using QR?

I'm still using QR because something is still registering as non-square.

What is the size of oprob.f.jac_prototype?

416x416 with jac=true and sparse=true, but 415x416 with jac=false and sparse=true, so that's the non-square matrix in question.

What is oprob.f.mass_matrix

Doing typeof(oprob.f.mass_matrix) returns UniformScaling{Bool}, so I guess it's just the identity matrix.

joegilkes commented 1 year ago

Re. sharing the model code, I'm more than happy to but it's pretty heavily tied up in a codebase which I'm unfortunately not at liberty to share right now. I can generate a loadable BSON of the ReactionSystem along with the activation energies for all the reactions, plus a script to assemble them all in a way that mimics the code, but it'll take a bit for me to get that ready.

isaacsas commented 1 year ago

Well if the jac=true and sparse=true case gives square try KLUFactorization there.

@ChrisRackauckas sounds like there might be a bug in the sparsity detection given that with jac=true the matrix is square but jac=false is giving non-square?

joegilkes commented 1 year ago

Sorry if I was unclear, I'd already switched back to KLU for the above tests, which was the condition that originally gave me the non-square matrix error. It works properly for the jac=true case.

isaacsas commented 1 year ago

So KLU with jac=true and sparse=true and autodiff=false is ok memory and compilation wise, but the solve is very slow?

joegilkes commented 1 year ago

I reran every combination since I changed up the temperature gradient to actually be continuous and I'd honestly lost track of what works and what doesn't. I've summed up my results so far below. I know it's not the most quantitative of benchmarks, I can work on proper timings and memory usage at a later date if there's interest. All runs were done with autodiff=false.

sparse jac linsolve result
false false default fast ODEProblem compilation, fast solver compilation, runs okay, slow solve, high memory usage (~ 15 GB)
true false default (KLUFactorization) fast ODEProblem compilation, fast solver compilation, errors at first timestep (LoadError: ArgumentError: KLU only accepts square matrices.)
false true default much slower ODEProblem compilation, much slower solver compilation, runs okay, fast solve, lower memory usage (~ 10 GB)
true true default (KLUFactorization) fast ODEProblem compilation, much slower solver compilation, runs okay, very fast solve, lower memory usage (~ 9 GB)
true false QRFactorization fast ODEProblem compilation, medium length solver compilation, errors at first timestep (LoadError: DimensionMismatch: array could not be broadcast to match destination in finite_difference_jacobian!()), very high memory usage (~ 12.5 GB)

I think I'm fine to stick with the jac=true, sparse=true case even though it takes much longer to compile, since I need to do a lot of these back to back and the lower memory usage will be more beneficial. I assume the DimensionMismatch in the final case is also a symptom of the non-square Jacobian, potentially where something is expecting it to be square but whatever is being passed in is rectangular?

This seems to be solved for me for now, happy to close the issue again or we can leave it open if you'd like to keep the discussion with Chris on here?

ChrisRackauckas commented 1 year ago

I mean some of the "issues" pointed out are just clear. QRFactorization is for dense matrices, KLUFactorization is for sparse matrices. So mixing them wit hthe wrong matrix will throw an error. The only one that is odd is sparse=true, jac = false, KLUFactorization should like that one, and we have many examples of that working. Can you show a reproducible MWE for that?

joegilkes commented 1 year ago

I thought using QRFactorization with sparse=true should automatically defer to SuiteSparse's SPQR solver, so it should still work?

In theory, if it errors on my scaled up code it should also error in the same way with the corrected MWE that we were using above. I have to go for a couple of hours, but I'll give it a go when I'm back.

joegilkes commented 1 year ago 
using Catalyst
using OrdinaryDiffEq
using IfElse

@variables t
@variables T(t) [isbcspecies=true]
@variables (species(t))[1:2]
@variables (k(t))[1:2] [isbcspecies=true, irreducible=false]

heatrate = 1.0
coolrate = -1.0

D = Differential(t)
function Tgrad(t)
    return IfElse.ifelse(t > 15.0, coolrate, 
        IfElse.ifelse(t < 5.0, heatrate, 0.0)
    )
end
function calc_rates(T)
    base_rates = [0.5, 0.5]
    # Only the reverse reaction is temperature dependent for the sake of this example
    return [base_rates[1], base_rates[2] * T]
end

@named rate_model = ODESystem([
    D(T) ~ Tgrad(t);
    k .~ calc_rates(T)
], t)

rxs = [
    Reaction(k[1], [species[1]], [species[2]], [1], [1]),
    Reaction(k[2], [species[2]], [species[1]], [1], [1])
]

T0 = 0.0
u0 = [10.0, 0.0]
u0map = Pair.(collect(species), u0)
kmap = Pair.(collect(k), calc_rates(T0))
u0map = vcat(u0map, kmap)
push!(u0map, Pair(T, T0))

tspan = (0.0, 20.0)

@named rs = ReactionSystem(rxs, t, collect(species), collect(k); checks=false, constraints=rate_model)
osys = structural_simplify(convert(ODESystem, rs))
oprob = ODEProblem(osys, u0map, tspan; sparse=true, jac=false)
sol = solve(oprob, QNDF())

The above MWE reproduces the bug @ChrisRackauckas is interested in.

joegilkes commented 1 year ago

Thanks @isaacsas, I'll close this here.

isaacsas commented 1 year ago

@joegilkes thanks for taking the time to work up that MWE.

Regarding compilation time, yes it gets big (or unusable) with Jacobians of larger systems. That is a known issue with large generated functions and will hopefully get addressed in the near future. Going to Krylov solvers can get around this, though I guess the current best-performing preconditioner in our benchmarking of chemical systems was ilu which did need an explicit matrix... That said, it should be possible to go fully matrix-free (I just don't know if this requires implementing your own preconditioner at the moment).

joegilkes commented 1 year ago

Thanks for the advice, It just so happens that in my workflow I end up compiling on a smaller system first and then scaling up and running again a lot of the time, which does massively reduce compilation time overall. However, I do have some larger networks (~ 10,000 reactions) which could pose a problem still, so I'll keep the matrix free approach in mind when I get to those.