Ability to move past unstable points in parameter space

xztt · June 5, 2022, 1:57pm

Hi,

This is a fairly general question but I can provide some more concrete code if needed.

Say I have 100 samples in a n dimensional parameter space between upper and lower bounds these bounds say 30% away from a “truth” parameter point, I am looking to solve an ODE system at each one of these samples in parameter space.

The issue I am encountering is that naturally not all of these parameter points leads to a stable solution, I have identified the parameter point which causes the error and simulated it separately and it returns the same error.

As I am iterating through these parameter samples is there any way for me to be able to mark that parameter sample as unstable and move onto the next one so the whole thing does not crash?

Cheers,
Harry

ChrisRackauckas · June 5, 2022, 3:37pm

Here’s the docs page for it: Handling Divergent and Unstable Trajectories · SciMLSensitivity.jl

xztt · June 5, 2022, 4:29pm

Looks great thank you however I am still encountering an issue the code I am using is

function circ_local(x)
    tmp_prob = remake(prob,u0=x[1:3],p=x[4:end])
    print(_prob.p) 
    print(_prob.u0)
    tmp_sol = solve(tmp_prob, Tsit5(), saveat=savetime)
    if tmp_sol.retcode == :Success
        return [(maximum(tmp_sol[LV.V])-minimum(tmp_sol[LV.V])),mean(tmp_sol[LV.V]),(maximum(tmp_sol[LV.p])-minimum(tmp_sol[LV.p])),(maximum(tmp_sol[C1.in.p])-minimum(tmp_sol[C1.in.p])),mean(tmp_sol[C1.in.p]), maximum(tmp_sol[C1.q]), maximum(tmp_sol[LV.p])/mean(tmp_sol[LV.p]),maximum(tmp_sol[C1.in.p])/mean(tmp_sol[C1.in.p]), mean(tmp_sol[LV.p])]
    else
    else
        return Inf
    end
end

S = zeros(9,12,100)
for i in 1:100
    S[:,:,i] = ForwardDiff.jacobian(circ_local,[su[:,i];sp[:,i]])
end

su and sp are 3x100 and 9x100 respectively each col containing the respective sample at that point. I should add I get the same error when I use

S[:,:,i] = Calculus.finite_difference_jacobian(circ_local,[u0;p])

I am still having the error returned to me.

ERROR: The nonlinear solver failed with the return code MAXITERS_EXCEED.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] nlsolve_failure(rc::Symbol)
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/vZUs5/src/structural_transformation/utils.jl:331
  [3] numerical_nlsolve(f::Function, u0::Float64, p::Tuple{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/vZUs5/src/structural_transformation/utils.jl:338
  [4] macro expansion
    @ ~/.julia/packages/ModelingToolkit/vZUs5/src/structural_transformation/codegen.jl:202 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:351 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
  [7] macro expansion
    @ ./none:0 [inlined]
  [8] generated_callfunc
    @ ./none:0 [inlined]
  [9] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#623"), Symbol("##arg#11434496959051576183"), Symbol("##arg#10207861574440612533"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x21a10f99, 0xa5ac6f87, 0xd3994102, 0x3f768d41, 0x75df6f19)})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:117
 [10] ODEFunction
    @ ~/.julia/packages/SciMLBase/chsnh/src/scimlfunctions.jl:1593 [inlined]
 [11] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Nothing, Float64, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Float64, Float64, Float64, Float64, Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}, ODESolution{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, 2, Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}}, ODEProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Tuple{Float64, Float64}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#623"), Symbol("##arg#11434496959051576183"), Symbol("##arg#10207861574440612533"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x21a10f99, 0xa5ac6f87, 0xd3994102, 0x3f768d41, 0x75df6f19)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#226"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#623"), Symbol("##arg#11434496959051576183"), Symbol("##arg#10207861574440612533"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x21a10f99, 0xa5ac6f87, 0xd3994102, 0x3f768d41, 0x75df6f19)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#226"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}, Vector{Float64}, Vector{Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#623"), Symbol("##arg#11434496959051576183"), Symbol("##arg#10207861574440612533"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x21a10f99, 0xa5ac6f87, 0xd3994102, 0x3f768d41, 0x75df6f19)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#226"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, LinRange{Float64, Int64}, Tuple{}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/irVAX/src/perform_step/low_order_rk_perform_step.jl:689
 [12] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEq/irVAX/src/perform_step/low_order_rk_perform_step.jl:675 [inlined]
 [13] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Nothing, Float64, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Float64, Float64, Float64, Float64, Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}, ODESolution{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, 2, Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}}, ODEProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Tuple{Float64, Float64}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#623"), Symbol("##arg#11434496959051576183"), Symbol("##arg#10207861574440612533"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x21a10f99, 0xa5ac6f87, 0xd3994102, 0x3f768d41, 0x75df6f19)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#226"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#623"), Symbol("##arg#11434496959051576183"), Symbol("##arg#10207861574440612533"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x21a10f99, 0xa5ac6f87, 0xd3994102, 0x3f768d41, 0x75df6f19)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#226"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}, Vector{Float64}, Vector{Vector{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}}, OrdinaryDiffEq.Tsit5Cache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#623"), Symbol("##arg#11434496959051576183"), Symbol("##arg#10207861574440612533"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x21a10f99, 0xa5ac6f87, 0xd3994102, 0x3f768d41, 0x75df6f19)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#226"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, LinRange{Float64, Int64}, Tuple{}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/irVAX/src/solve.jl:477
 [14] __solve(::ODEProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, Tuple{Float64, Float64}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#623"), Symbol("##arg#11434496959051576183"), Symbol("##arg#10207861574440612533"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x21a10f99, 0xa5ac6f87, 0xd3994102, 0x3f768d41, 0x75df6f19)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#226"{Bool, ModelingToolkit.BipartiteGraphs.Matching{ModelingToolkit.BipartiteGraphs.Unassigned, Vector{Union{ModelingToolkit.BipartiteGraphs.Unassigned, Int64}}}, TearingState{ODESystem}, Dict{Any, Any}, BitVector, Vector{SymbolicUtils.Code.Assignment}, Tuple{Vector{Vector{Int64}}, Vector{BitSet}}, SymbolicUtils.Code.NameState, Dict{Any, Int64}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, LinRange{Float64, Int64}, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{LinRange{Float64, Int64}}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/irVAX/src/solve.jl:5
 [15] #solve_call#28
    @ ~/.julia/packages/DiffEqBase/LsFow/src/solve.jl:388 [inlined]
 [16] #solve_up#30
    @ ~/.julia/packages/DiffEqBase/LsFow/src/solve.jl:686 [inlined]
 [17] #solve#29
    @ ~/.julia/packages/DiffEqBase/LsFow/src/solve.jl:670 [inlined]
 [18] circ_local(x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}})
    @ Main ~/Desktop/PhD/Year 1/Sem 1/Valve and Ventricle analysis/AcausalModel_refactored.jl:1205
 [19] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/wAaVJ/src/apiutils.jl:37 [inlined]
 [20] vector_mode_jacobian(f::typeof(circ_local), x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/wAaVJ/src/jacobian.jl:148
 [21] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/wAaVJ/src/jacobian.jl:21
 [22] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(circ_local), Float64}, Float64, 12}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/wAaVJ/src/jacobian.jl:19
 [23] top-level scope
    @ ~/Desktop/PhD/Year 1/Sem 1/Valve and Ventricle analysis/AcausalModel_refactored.jl:1222

I could be missing something but unsure

xztt · June 7, 2022, 10:19am

After fixing some quite embarrassing errors in the code I am still facing the same problem however I have narrowed it down to where the problem may be

I tried the below

function circ_local(x)
    tmp_prob = remake(prob,u0=x[1:3],p=x[4:end])
    tmp_sol = solve(tmp_prob, Rodas4P(),  saveat=savetime)
    if tmp_sol.retcode == :Success
        return [(maximum(tmp_sol[LV.V])-minimum(tmp_sol[LV.V])),mean(tmp_sol[LV.V]),(maximum(tmp_sol[LV.p])-minimum(tmp_sol[LV.p])),(maximum(tmp_sol[C1.in.p])-minimum(tmp_sol[C1.in.p])),mean(tmp_sol[C1.in.p]), maximum(tmp_sol[C1.q]), maximum(tmp_sol[LV.p])/mean(tmp_sol[LV.p]),maximum(tmp_sol[C1.in.p])/mean(tmp_sol[C1.in.p]), mean(tmp_sol[LV.p])]
    else 
        return [0,0,0,0,0,0,0,0,0]
    end
end 

S = zeros(9,12,100)

for i in 1:100
    x = [su[:,i];sp[:,i]]
    S[:,:,i] = ForwardDiff.jacobian(circ_local,x)
    if S[:,:,i].retcode == :success
        return ForwardDiff.jacobian(circ_local,x)
    else 
        return zeros(9,12)
    end 
    #S[:,:,i] = Calculus.finite_difference_jacobian(circ_local,x)
end

However Array has no code ret code, I feel like the function is ignoring how I circle through the different parameter samples so maybe it is a problem as I am trying to incorporate arrays? Just thought I would add these here as any others looking it may lend a hand

ChrisRackauckas · June 7, 2022, 10:25am

It looks like you generated it via ODAEProblem. What happens if you generate the code from MTK using ODEProblem?

xztt · June 7, 2022, 11:24am

Using ODE problem I am actually have a problem with the solution been returned as unstable however returns fine when using ODAE problem.

The warning is

Warning: dt(2.220446049250313e-16) <= dtmin(2.220446049250313e-16) at t=0.3051060816585229. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/chsnh/src/integrator_interface.jl:366
  0.002520 seconds (3.80 k allocations: 139.750 KiB)
retcode: DtLessThanMin
Interpolation: specialized 3rd order "free" stiffness-aware interpolation

I have tried with a variety of initial conditions and solvers model equations after structural simplify are

 Differential(t)(LV₊p(t)) ~ ((LV₊Eₘᵢₙ + (LV₊Eₘₐₓ - LV₊Eₘᵢₙ)*((1//2)*(1 - cos((3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ)) / LV₊τₑₛ))*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₛ) + (1//2)*(1 + cos((3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ) - 3.141592653589793LV₊τₑₛ) / (LV₊τₑₚ - LV₊τₑₛ)))*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₚ)*(rem(t + 0.999999LV₊τ, LV₊τ) > LV₊τₑₛ)))*(Rd₊q(t) - AV₊q(t))) / (1.0 + 1.0e-6LV₊Eₘᵢₙ + 1.0e-6(LV₊Eₘₐₓ - LV₊Eₘᵢₙ)*((1//2)*(1 - cos((3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ)) / LV₊τₑₛ))*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₛ) + (1//2)*(1 + cos((3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ) - 3.141592653589793LV₊τₑₛ) / (LV₊τₑₚ - LV₊τₑₛ)))*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₚ)*(rem(t + 0.999999LV₊τ, LV₊τ) > LV₊τₑₛ))) + ((LV₊Eₘₐₓ - LV₊Eₘᵢₙ)*((1.5707963267948966sin((3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ)) / LV₊τₑₛ)*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₛ)) / LV₊τₑₛ + (1.5707963267948966sin((3.141592653589793LV₊τₑₛ - 3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ)) / (LV₊τₑₚ - LV₊τₑₛ))*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₚ)*(rem(t + 0.999999LV₊τ, LV₊τ) > LV₊τₑₛ)) / (LV₊τₑₚ - LV₊τₑₛ))*LV₊p(t)) / ((LV₊Eₘᵢₙ + (LV₊Eₘₐₓ - LV₊Eₘᵢₙ)*((1//2)*(1 - cos((3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ)) / LV₊τₑₛ))*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₛ) + (1//2)*(1 + cos((3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ) - 3.141592653589793LV₊τₑₛ) / (LV₊τₑₚ - LV₊τₑₛ)))*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₚ)*(rem(t + 0.999999LV₊τ, LV₊τ) > LV₊τₑₛ)))*(1.0 + 1.0e-6LV₊Eₘᵢₙ + 1.0e-6(LV₊Eₘₐₓ - LV₊Eₘᵢₙ)*((1//2)*(1 - cos((3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ)) / LV₊τₑₛ))*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₛ) + (1//2)*(1 + cos((3.141592653589793rem(t + 0.999999LV₊τ, LV₊τ) - 3.141592653589793LV₊τₑₛ) / (LV₊τₑₚ - LV₊τₑₛ)))*(rem(t + 0.999999LV₊τ, LV₊τ) <= LV₊τₑₚ)*(rem(t + 0.999999LV₊τ, LV₊τ) > LV₊τₑₛ))))
 0 ~ -Rd₊Δp(t) - Rd₊R*Rd₊q(t)
 Differential(t)(L1₊q(t)) ~ (-L1₊Δp(t)) / L1₊L
 Differential(t)(C1₊Δp(t)) ~ (-C1₊q(t)) / C1₊C
 0 ~ AV₊q(t) - C1₊q(t) - Rd₊q(t)

I will try with a toy problem though to see where the issue lays.

xztt · June 7, 2022, 2:38pm

Managed to get my model to run just by messing with some with some parameter values however for this problem that is not really relevant.

When the problem is generated using ODAE the error of MAXITERS appears. However, when my problem is generated using ODE problem, using my above code, the parameter samples for which the solution diverges the output is a matrix of zeros as specified by my code.

So I think the issue is when I generate the problem using ODAE. Do you know of a work round for this?

ChrisRackauckas · June 7, 2022, 9:40pm

No, can you open an issue? Indeed, when you generate an ODAE problem function, it embeds nonlinear solvers inside of the RHS. You seem to have a case where one of those internal nlsolvers diverges, and it seems it throws an error instead of a warning (to be safe, that’s good at least). The way normal nonlinear iterations work in an ODE solver though, a nonlinear solver divergence would trigger a step decrease, but this cannot do it because it’s baked into the f. We need to do a few things here, also memory for the initial conditions, so if you can open an issue this would serve as a nice MWE for handling some of this stuff.

xztt · June 8, 2022, 8:54am

Of course, issue can be found here.

https://github.com/SciML/ModelingToolkit.jl/issues/1633