Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Jump to bottom

ODE-jump models with implicit solvers not working #459

Closed
isaacsas opened this issue Oct 28, 2024 · 17 comments · Fixed by SciML/LinearSolve.jl#555
Closed

ODE-jump models with implicit solvers not working #459

isaacsas opened this issue Oct 28, 2024 · 17 comments · Fixed by SciML/LinearSolve.jl#555
Labels
bug

Comments

@isaacsas
Copy link
Member

Motivated by https://discourse.julialang.org/t/jumpprocesses-jl-solvers/121883 which seems to get an error instead, I am finding

using JumpProcesses, OrdinaryDiffEq
p = [2.0, 1.0]
u0 = [1.0]
tspan = (0.0, 4.0)

function ode_fxn(du, u, p, t)
    du .= 0
    nothing
end

b_rate(u, p, t) = (u[1] * p[1])
function birth!(integrator)
    integrator.u[1] += 1
    nothing
end
b_jump = VariableRateJump(b_rate, birth!)

d_rate(u, p, t) = (u[1] * p[2])
function death!(integrator)
    integrator.u[1] -= 1
    nothing
end
d_jump = VariableRateJump(d_rate, death!)

ode_prob = ODEProblem(ode_fxn, u0, tspan, p)
sjm_prob = JumpProblem(ode_prob, b_jump, d_jump)
sol = solve(sjm_prob, ImplicitEuler())

gives

┌ Warning: At t=0.0, dt was forced below floating point epsilon 5.0e-324, and step error estimate = 1.0. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase ~/.julia/packages/SciMLBase/fqdZV/src/integrator_interface.jl:623
retcode: Unstable
Interpolation: 3rd order Hermite
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}:
 [1.0]

while

julia> sol = solve(sjm_prob, Rodas5())
retcode: Success
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 724-element Vector{Float64}:
 0.0
 0.004127525454249106
 0.0059311904982310764
 0.0077675951067161225
 0.009603999715201168
 0.011440404323686213
 0.013276808932171258
 0.015113213540656303
 0.016949618149141348
 0.018786022757626395
 0.020622427366111442
 0.02245883197459649
 0.024295236583081536
 0.026131641191566583
 0.02796804580005163
 0.029804450408536676
 0.03164085501702172
 0.03347725962550677
 0.035313664233991814
 0.03715006884247686
 0.03898647345096191
 0.040822878059446954
 0.042659282667932
 0.04449568727641705
 0.046332091884902095
 0.04816849649338714
 0.05000490110187219
 0.051841305710357236
 0.05367771031884228
 0.05551411492732733
 0.057350519535812376
 0.05918692414429742
 0.06102332875278247
 ⋮
 3.535887868639177
 3.550138043902331
 3.5644426672952085
 3.5788019468572734
 3.593216091422883
 3.6076853116218945
 3.622209817895197
 3.6367898203665967
 3.6514255310796524
 3.666117162888113
 3.680864929459014
 3.6956690452757828
 3.7105297256413605
 3.7254471866813312
 3.7404216453470664
 3.7554533194188795
 3.7705424275091937
 3.785689187900973
 3.8008938220407074
 3.816156549882591
 3.831477593390414
 3.8468571765586956
 3.8622955218724875
 3.877792853849562
 3.8933493990617847
 3.9089653825541872
 3.9246410314279094
 3.94037657486182
 3.9561722404903854
 3.9720282580278385
 3.987944858066102
 4.0
u: 724-element Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}:
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 ⋮
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]

But this works fine with Tsit5().
@isaacsas isaacsas added the bug label Oct 28, 2024
@isaacsas
Copy link
Member Author

For reference, the mean solution here is exp(t) on [0,4] which Tsit5() captures.

@hhindley
Copy link

hhindley commented Oct 29, 2024 •
edited
Loading

With this code:

using DifferentialEquations, JumpProcesses

function ode!(du, u, p, t)
    du[1] = p[1] * u[1]
end

u0 = [1.0, 10.0, 5.0, 3.0, 7.0, 2.0, 6.0, 4.0, 8.0, 1.0]

tspan = (0.0, 100.0)

p = [0.05, 0.1, 0.2, 0.05, 0.03] 

ode_prob = ODEProblem(ode!, u0, tspan, p)

S1, S2, S3, S4, S5, S6, S7, S8, S9 = 1, 2, 3, 4, 5, 6, 7, 8, 9

k_degrade = 0.1  
k_produce = 0.2  
k_massaction1 = 0.05
k_massaction2 = 0.03

# variable rate jumps 
rate1(u,p,t) = k_degrade * u[S1]
rate2(u,p,t) = k_degrade * u[S2]
rate3(u,p,t) = k_degrade * u[S3]
rate4(u,p,t) = k_degrade * u[S4]
rate5(u,p,t) = k_degrade * u[S5]
rate6(u,p,t) = k_degrade * u[S6]
function affect1!(integrator)
    integrator.u[S1] -= 1
end
function affect2!(integrator)
    integrator.u[S2] -= 1
end
function affect3!(integrator)
    integrator.u[S3] -= 1
end
function affect4!(integrator)
    integrator.u[S4] -= 1
end
function affect5!(integrator)
    integrator.u[S5] -= 1
end
function affect6!(integrator)
    integrator.u[S6] -= 1
end

r1 = VariableRateJump(rate1, affect1!)
r2 = VariableRateJump(rate2, affect2!)
r3 = VariableRateJump(rate3, affect3!)
r4 = VariableRateJump(rate4, affect4!)
r5 = VariableRateJump(rate5, affect5!)
r6 = VariableRateJump(rate6, affect6!)

# constant rate jumps 
rate7(u,p,t) = k_produce
rate8(u,p,t) = k_produce
function affect7!(integrator)
    integrator.u[S7] += 1
end
function affect8!(integrator)
    integrator.u[S8] += 1
end
# Constant production of S7 and S8
r7 = ConstantRateJump(rate7, affect7!)
r8 = ConstantRateJump(rate8, affect8!)

# mass action jumps
reactant_stoich = [
    [S1 => 1, S2 => 1],   # S1 + S2
    [S3 => 1, S4 => 1],   # S3 + S4
    [S5 => 1, S6 => 1],   # S5 + S6
    [S1 => 1, S7 => 1],   # S1 + S7
    [S2 => 1, S8 => 1]    # S2 + S8
]

net_stoich = [
    [S1 => -1, S2 => -1, S9 => 1],   # S1 + S2 -> S9
    [S3 => -1, S4 => -1, S7 => 1],   # S3 + S4 -> S7
    [S5 => -1, S6 => -1, S8 => 1],   # S5 + S6 -> S8
    [S1 => -1, S7 => -1, S9 => 1],   # S1 + S7 -> S9
    [S2 => -1, S8 => -1, S6 => 1]    # S2 + S8 -> S6
]


r9 = MassActionJump(reactant_stoich, net_stoich, param_idxs = [1])

jumps = JumpSet(r1, r2, r3, r4, r5, r6, r7, r8, r9)

jump_prob = JumpProblem(ode_prob, Direct(), jumps)

sol = solve(jump_prob, Rodas4())

I get the error:

ERROR: MethodError: no method matching strides(::ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}})

Closest candidates are:
  strides(::StaticArrayInterface.AbstractArray2)
   @ StaticArrayInterface ~/.julia/packages/StaticArrayInterface/lkDPR/src/StaticArrayInterface.jl:285
  strides(::StaticArrayInterface.AbstractArray2, ::Any)
   @ StaticArrayInterface ~/.julia/packages/StaticArrayInterface/lkDPR/src/StaticArrayInterface.jl:289
  strides(::PermutedDimsArray{T, N, perm}) where {T, N, perm}
   @ Base permuteddimsarray.jl:64
  ...

Stacktrace:
  [1] stride(A::ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}, k::Int64)
    @ Base ./abstractarray.jl:595
  [2] stride1(x::ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/LinearAlgebra.jl:215
  [3] _chkstride1 (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/LinearAlgebra.jl:221 [inlined]
  [4] chkstride1
    @ ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/LinearAlgebra.jl:219 [inlined]
  [5] aa_getrs!(trans::Char, A::Matrix{…}, ipiv::Vector{…}, B::ExtendedJumpArray{…}; info::Base.RefValue{…})
    @ LinearSolve ~/.julia/packages/LinearSolve/B82H3/src/appleaccelerate.jl:170
  [6] solve!(cache::LinearSolve.LinearCache{…}, alg::AppleAccelerateLUFactorization; kwargs::@Kwargs{…})
    @ LinearSolve ~/.julia/packages/LinearSolve/B82H3/src/appleaccelerate.jl:258
  [7] macro expansion
    @ ~/.julia/packages/LinearSolve/B82H3/src/default.jl:363 [inlined]
  [8] solve!(::LinearSolve.LinearCache{…}, ::LinearSolve.DefaultLinearSolver; assump::OperatorAssumptions{…}, kwargs::@Kwargs{…})
    @ LinearSolve ~/.julia/packages/LinearSolve/B82H3/src/default.jl:356
  [9] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqRosenbrock.RosenbrockCache{…}, repeat_step::Bool)
    @ OrdinaryDiffEqRosenbrock ~/.julia/packages/LinearSolve/B82H3/src/default.jl:0
 [10] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEqRosenbrock/8sxst/src/rosenbrock_perform_step.jl:1300 [inlined]
 [11] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/hCoet/src/solve.jl:551
 [12] __solve(jump_prob::JumpProblem{…}, alg::Rodas4{…}, timeseries::Vector{…}, ts::Vector{…}, ks::Vector{…}, recompile::Type{…}; kwargs::@Kwargs{})
    @ JumpProcesses ~/.julia/packages/JumpProcesses/MBJ3Z/src/solve.jl:6
 [13] __solve (repeats 5 times)
    @ ~/.julia/packages/JumpProcesses/MBJ3Z/src/solve.jl:1 [inlined]
 [14] #solve#58
    @ ~/.julia/packages/DiffEqBase/slKc5/src/solve.jl:1109 [inlined]
 [15] solve(prob::JumpProblem{…}, args::Rodas4{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/slKc5/src/solve.jl:1108
 [16] top-level scope
    @ ~/phd/gillespe_stochastic/examples/test.jl:92
Some type information was truncated. Use `show(err)` to see complete types.

However using solve(jump_prob, Tsit5()) the model is solved successfully. I could not copy my exact code here due to confidentiality reasons, but I tried to reproduce the error. It seems to only appear with an increased number of species and therefore reactions, with a reduced model of 6 species I did not get this error.

@isaacsas
Copy link
Member Author

Very strange. This is a really busy week for me, so I might not get to help investigate this till Friday or the weekend. If you wanted to play around you could try changing the linear solver to a matrix-free Krylov method and see if that helps: https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Using-Jacobian-Free-Newton-Krylov

@ChrisRackauckas
Copy link
Member

For both of you, are you on a Mac or Intel machine?

The following should work for you for now:

sol = solve(jump_prob, Rodas4(linsolve = LUFactorization()))
sol = solve(jump_prob, Rodas4(linsolve = RFLUFactorization()))

Since I just moved I don't have all of my machines around, but I can see on my Mac that the issue is that we too readily default to AppleAccelerate overloads on a non-contiguous b vector.

@ChrisRackauckas
Copy link
Member

The other case that might be problematic which I cannot test right now is MKL, which is why I asked for if anyone has an Intel. Even if it's just an X86 though you can at least try and check:

sol = solve(jump_prob, Rodas4(linsolve = MKLLUFactorization()))

My guess is that it will likely have the same behavior as the Apple one.

I'm trying to think though if there's a better solution than just not defaulting to the tuned BLAS in these cases. That might be the easiest robust solution though.

@isaacsas
Copy link
Member Author

I'm on a mac.

@isaacsas
Copy link
Member Author

@ChrisRackauckas do you have any idea about the issue with the simple birth-death example I posted? Is it degenerate for some reason due to du .= 0? (I don't see why that should cause problems though -- the solutions show lots of time points but no jumps occurring.)

@isaacsas
Copy link
Member Author

(We actually use that one in the tests, so once it is working I will update the tests to check correctness against one or more implicit methods too.)

@hhindley
Copy link

I am also on a mac, using those above solvers does get rid of the error but I still hit MaxIters when run for a longer timespan.

@ChrisRackauckas
Copy link
Member

It's degenerate because the starting rate value ends up negative:

julia> sol[1]
3-element ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}:
  1.0
 -1.2064836293851149
 -0.8345225253385252

so it never rootfinds. So @isaacsas your case is different from @hhindley 's. Why that starts negative is a bit of a mystery though

@isaacsas
Copy link
Member Author

The method relies on them being negative initially though? (i.e. we use the initial values -randexp(), and then the idea is it integrates the intensity until the time where the integral plus these values are zero.)

@ChrisRackauckas
Copy link
Member

Then the rate should be positive, but we somehow have a negative rate here.

@isaacsas
Copy link
Member Author

isaacsas commented Oct 30, 2024 •
edited
Loading

Could the nonlinear and/or linear solve be introducing small negative values even though f(du,u,p,t) just sets du .= 0? (Perhaps from the coupling to the integrated intensity values, which are not zero.)

@isaacsas
Copy link
Member Author

I am also on a mac, using those above solvers does get rid of the error but I still hit MaxIters when run for a longer timespan.

But is this now consistent for both implicit and explicit solvers? If so, is it due to the frequency of jump firing? You can't "skip over" jumps with an implicit solver, so you'll still need to resolve each one. If you don't need / want to resolve each exact jump you'd need to use a tau-leaping approach instead.

@ChrisRackauckas
Copy link
Member

This is a bug. I haven't taken the time to figure out why it's going in the wrong direction, but it is.

@ChrisRackauckas
Copy link
Member

Could the nonlinear and/or linear solve be introducing small negative values even though f(du,u,p,t) just sets du .= 0? (Perhaps from the coupling to the integrated intensity values, which are not zero.)

Yes that's possible due to LU factorization error growth. SciML/OrdinaryDiffEq.jl#1651 details exactly how that can happen.

@ChrisRackauckas
Copy link
Member 

using JumpProcesses, OrdinaryDiffEq
p = [2.0, 1.0]
u0 = [1.0]
tspan = (0.0, 4.0)

function ode_fxn(du, u, p, t)
    du .= 0
    nothing
end

b_rate(u, p, t) = (u[1] * p[1])
function birth!(integrator)
    integrator.u[1] += 1
    nothing
end
b_jump = VariableRateJump(b_rate, birth!)

d_rate(u, p, t) = (u[1] * p[2])
function death!(integrator)
    integrator.u[1] -= 1
    nothing
end
d_jump = VariableRateJump(d_rate, death!)

ode_prob = ODEProblem(ode_fxn, u0, tspan, p)
sjm_prob = JumpProblem(ode_prob, b_jump, d_jump)

using LinearSolve
sol = solve(sjm_prob, Rodas5P(linsolve = QRFactorization()))
retcode: Success
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 827-element Vector{Float64}:
 0.0
 0.033531107529250564
 0.033531107529250564
 ⋮
 3.9993533353919233
 3.9993533353919233
 4.0
u: 827-element Vector{ExtendedJumpArray{Float64, 1, Vector{Float64}, Vector{Float64}}}:
 [1.0]
 [1.0]
 [2.0]
 ⋮
 [115.99999999999999]
 [114.99999999999999]
 [114.99999999999999]

Correcting for the numerical error via QRFactorization fixes @isaacsas 's case as expected due to the issue linked above.

@ChrisRackauckas ChrisRackauckas mentioned this issue Oct 30, 2024
Merged
@isaacsas isaacsas mentioned this issue Nov 11, 2024
Merged
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
bug
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Guard BLAS defaults more, allowing only Array and BLASFloats SciML/LinearSolve.jl
3 participants
@ChrisRackauckas @isaacsas @hhindley