Solving highly stiff chemical kinetics IVP using DiffEqBiological

[bjodah]

Hi,

I’ve been looking into using DifferentialEquations.jl to solve some IVPs (chemical kinetics problems more precisely, something I have extensive experience working with using the scientific python stack).

I was delighted to see there’s a DSL for chemical kinetics: DiffEqBiological

So I went ahead trying to model a problem of mine which has been demanding for Boost’s odeint, GNU Scientific Library, as well as Sundials’ CVode (unless tweaking and tuning tolerances & constraints extensively).

I find that the Rosenbrock23 algorithm as implemented in DifferentialEquations fared best. Here’s my problem definition:

using DifferentialEquations, DiffEqBiological, Plots

rn = @reaction_network begin
    k1*A^2, 2B + 2A ⇒ C + 2D
    k2, 2E --> C
    k3, 2F --> G
    k4*E*A, E + B + A ⇒ C + D
    k5, A + F --> D
    k6, E + F --> B
    k7, A + G --> F + D
    k8, A + H --> I
    k9*I*A, 2B + I + A ⇒ G + 2D
    k10, A + J --> K
    k11, E + G --> F + B
    k12, E + H --> J
    k13, E + J --> 2F
    k14, E + I --> K
    k15, F + G --> J + B
    k16, F + I --> H + D
    k17, F + J --> H + B
    k18, 2J --> G + H
    k19*I*J, B + J + I ⇒ G + H + D
    k20*G, 2G ⇒ H + 2B
    k21, F + K --> B + I
    k22, M + G --> D + J
    k23, M + K --> D + I
    k24, M + C --> D + E
    k25, B --> L + D
    k26, L + D --> B
    k27, G --> L + K
    k28, L + K --> G
    k29, G + D --> K + B
    k30, K + B --> G + D
    k31, F --> L + M
    k32, L + M --> F
    k33, F + D --> M + B
    k34, M + B --> F + D
    k35, J --> L + I
    k36, L + I --> J
    k37, J + D --> I + B
    k38, I + B --> J + D
    k39, E --> L + A
    k40, L + A --> E
    k41, E + D --> A + B
    k42, A + B --> E + D
    k43, E + B --> C + F
    k44, C + F --> E + B
    k45, M + H --> N
    k46, N --> M + H
    k47, S + D --> B + R
    k48, B + R --> S + D
    k49, R + F --> B + Q
    k50, R + M --> Q + D
    k51, E + R --> S + A
    k52, S + A --> E + R
    k53, G + Q --> J + R
    k54, K + Q --> R + I
    k55, J + Q --> R + H
    k56, E + Q --> R
    k57*H*T^2, 2B + 4T + H ⇒ 4L + 4V
    k58*E*T, E + 2T ⇒ P + F
    k59, U + F --> L + W
    k60, U + I --> V + H
    k61*U^2, B + 2U ⇒ 2L + V + W
    k62, V + F --> U + D
    k63, W + A --> X
    k64*X, B + X ⇒ U + 2D
    k65, S + X --> R + U + D
    k66, L + X --> U + D
    k67, X + H --> W + I
    k68, P + A --> O + M
    k69*A*V, B + V + A ⇒ T + 2D
    k70, Q + H --> B + T
    k71, U + M --> W
    k72, T + I --> W
    k73*T*U, B + T + U ⇒ 2L + 2V
    k74, 0 ⇒ A
    k75, 0 ⇒ G
    k76, 0 ⇒ C
    k77, 0 ⇒ F
    k78, 0 ⇒ E
    k79, 0 ⇒ L
    k80, B ⇒ 0
    k81, 0 ⇒ H
end k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20 k21 k22 k23 k24 k25 k26 k27 k28 k29 k30 k31 k32 k33 k34 k35 k36 k37 k38 k39 k40 k41 k42 k43 k44 k45 k46 k47 k48 k49 k50 k51 k52 k53 k54 k55 k56 k57 k58 k59 k60 k61 k62 k63 k64 k65 k66 k67 k68 k69 k70 k71 k72 k73 k74 k75 k76 k77 k78 k79 k80 k81
p = Dict([(:k1, 7.226e+09), (:k2, 5.126e+09), (:k3, 4.806e+09), (:k4, 2.754e+10), (:k5, 3.546e+10), (:k6, 1.092e+10), (:k7, 1.357e+10), (:k8, 2.285e+10), (:k9, 1.295e+10), (:k10, 1.295e+10), (:k11, 3.635e+07), (:k12, 1.304e+10), (:k13, 1.135e+10), (:k14, 1.135e+10), (:k15, 2.911e+07), (:k16, 1.095e+10), (:k17, 8.828e+09), (:k18, 8.365e+05), (:k19, 1.01e+08), (:k20, 1.281e-07), (:k21, 4.057e+09), (:k22, 4.057e+09), (:k23, 7.825e+08), (:k24, 1.276e+08), (:k25, 2.086e-05), (:k26, 1.174e+11), (:k27, 0.09327), (:k28, 5.007e+10), (:k29, 1.326e+10), (:k30, 1.264e+06), (:k31, 0.09327), (:k32, 5.007e+10), (:k33, 1.326e+10), (:k34, 1.264e+06), (:k35, 7.692e+05), (:k36, 5.007e+10), (:k37, 1.326e+10), (:k38, 0.1533), (:k39, 5.752), (:k40, 2.092e+10), (:k41, 2.421e+07), (:k42, 15.64), (:k43, 4.481e-05), (:k44, 3.937e+07), (:k45, 3.739e+09), (:k46, 2593), (:k47, 1.326e+10), (:k48, 2.407e+05), (:k49, 7.982e+07), (:k50, 7.982e+06), (:k51, 9.935e+05), (:k52, 6.659e+05), (:k53, 8.974e+07), (:k54, 8.974e+07), (:k55, 8.974e+07), (:k56, 1.994e+10), (:k57, 5.983e+06), (:k58, 1.097e+10), (:k59, 4.487e+09), (:k60, 9.972e+07), (:k61, 9.972e+07), (:k62, 9.672e+08), (:k63, 9.6e+09), (:k64, 9.972e+04), (:k65, 1.994e+08), (:k66, 1.994e+10), (:k67, 1.994e+08), (:k68, 9.6e+09), (:k69, 3.576e+09), (:k70, 1.097e+09), (:k71, 3.49e+09), (:k72, 6.681e+09), (:k73, 7.24e+06), (:k74, 4.266e-08), (:k75, 1.104e-08), (:k76, 6.793e-09), (:k77, 4.353e-08), (:k78, 9.373e-09), (:k79, 4.266e-08), (:k80, 6.561e-08), (:k81, 4.066e-13)])
ics = Dict([(:B, 55.4), (:R, 0.1), (:L, 1e-11), (:D, 0.001), (:A, 1e-24), (:C, 1e-24), (:E, 1e-24), (:F, 1e-24), (:G, 1e-24), (:H, 1e-24), (:I, 1e-24), (:J, 1e-24), (:K, 1e-24), (:M, 1e-24), (:N, 1e-24), (:O, 1e-24), (:P, 1e-24), (:Q, 1e-24), (:S, 1e-24), (:T, 1e-24), (:U, 1e-24), (:V, 1e-24), (:W, 1e-24), (:X, 1e-24)])

Solving this succeeds (albeit at somewhat high computational cost, and I’d actually like tspan to go up to 1e13 – something which currently fails):

tspan = (0., 3.6e7)
u0 = Array([get(ics, k, 0.0) for k=keys(speciesmap(rn))])
parr = Array([p[k] for k=keys(paramsmap(rn))])
oprob = ODEProblem(rn, u0, tspan, parr)

sol = solve(oprob, reltol=1e-9, abstol=1e-10,
    Rosenbrock23()
)
sol.retcode

Looking at sol.destats I see ~100k rhs evaluations and ~30k jacobian creations (to be compared with ~30k rhs evals & ~3k jacobian evaluations when solving using my hand-crafted CVode solver).

Plotting the result gives a hint why the computational cost is high (high frequency oscillatory behavior in a region which should be smooth):

gr()
lgt = log10.(clamp.(sol.t[2:end], 0, Inf))
lgu = log10.(clamp.(hcat(sol.u[2:end]...)', 0, Inf))
plot(lgt, lgu, fmt = :png, legend=:topleft, legendfontsize=8)

lowering the tolerances seem to make the solution deteriorate.

So my question is: what other options do DifferentialEquations/DiffEqBiological offer for tuning the solution behavior (e.g. are positivity constraints implied from using DiffEqBiological, or does one need to pass that infomation on when solving)?

[ChrisRackauckas]

DiffEqBiological master now automatically will form sparse Jacobian information and functions. Let me get that tagged. @YingboMa can take a look at this. I would think that something like Rodas5 should end up faster than Rosenbrock23, and sparse Jacobians should really improve it.

@isaacsas if you have time too.

@bjodah

Looking at sol.destats I see ~100k rhs evaluations and ~30k jacobian creations (to be compared with ~30k rhs evals & ~3k jacobian evaluations when solving using my hand-crafted CVode solver).

Rosenbrock23 doesn’t take that large of steps, while Rodas5 can. Changing the algorithm would probably drop the evaluations a lot. A Rosenbrock method will always have more Jacobian creations though since it cannot do Jacobian reuse (though a W-method can). However, at small enough systems like this, making a new/more accurate Jacobian and re-factorizing can be faster (while at larger sizes, minimizing the number of factorizations really is key).

[bjodah]

Thanks!

Rodas5 was indeed somewhat more efficient than Rosenbrock23 when I re-tried it (I did find Rosenbrock23 to be more robust though during some initial experimenting).

For reference, I can solve for tspan=(0, 3.6e13) using CVode:

I would of course be delighted if I could match this using Julia’s stack which you have been developing.

[isaacsas]

@bjodah, when you say you are solving with CVODE do you mean in your custom C++ code, or using the CVODE_BDF2 solver within Julia? If the former, could you try solving with the Julia interface to CVODE_BDF2 and see what that gives?

BTW, with the latest DiffEqBio release I’m seeing that SymEngine is somehow changing rate laws like k4*E to k4*exp(1). I’m not sure if this is a problem for the older version you are using, but this is definitely a bug we’ll have to fix or you won’t get the right reaction network.

[bjodah]

@isaacsas, I was referring to my own custom C++ code. Sure, I’ll try CVODE_BDF2 and report back.

Ah, that’s unfortunate, symengine-0.4.0 includes a new parser (https://github.com/symengine/symengine/pull/1211) which could be to blame for that.

[isaacsas]

As long as you don’t update to DiffEqBio v3.8 I think you should be ok (the old version converted to internal symbol names and hence avoided this issue I think – I’ll try to update the current code to do something similar in the next few days).

[isaacsas]

One other gotcha to look out for – as the size of chemical systems increase, the time to compile the underlying Jacobian evaluation function (whether sparse or dense) can grow substantially. If you manually call the ODE f function and Jacobian function generated by DiffEqBio one time before calling solve, you can factor out their contribution to the running time (as this removes the compilation overhead from the solve call). This is useful in understanding if you are seeing compilation delays, or performance issues with the ODE solvers.

[j_verzani]

E is exported here https://github.com/symengine/SymEngine.jl/blob/master/src/mathops.jl#L45

You may need to rename the variable internally in the code.

[bjodah]

@isaacsas, so I’ve prepended all substance keys with a “s” and tried with CVODE_BDF() (using Sundials doesn’t seem to export any CVODE_BDF2 ?).

It’s actually considerably more efficient than either Rodas5 or Rosenbrock23 for the tolerances I tried. Here’s what I ended up with:
http://hera.physchem.kth.se/~bjorn/dahlgren_ode_2019-Copy1.htm

It still fails at integrating up to 1e13 though…

[isaacsas]

@bjodah Are you finding that the Julia solutions have negative values in them?

In your C++ code do you have constraints to enforce positivity? DEBio does not add positivity-preserving callbacks.

Is your C++ code also using CVODE with an explicit Jacobian function and the same absolute/relative error tolerances?

Some things you could check to try to isolate the problem:

Check that the generated ODE rhs function in Julia gives similar derivative values for a test vector as your C++ code. If it isn’t this would indicate a bug in your network specification, or in DEBio.

Check that the generated ODE Jacobian function in Julia gives similar values for a test vector as your C++ code.

Finally, you might try using a min_reaction_network and then using addodes! with the build_jac=false kwarg to construct the ODEProblem without a Jacobian. i.e. something like:

rn = @min_reaction_network begin
...
end
# setup initial condition and parameters as you currently do

# before calling ODEProblem call addodes!
addodes!(rn; build_jac=false)
oprob = ODEProblem(rn, u0, tspan, parr)
...

This would check if the problem is with the generated Jacobian from DEBio.

If you post a text version of your updated network I can take a look at it, for some reason I can’t select text in the html version you linked above.

[bjodah]

@isaacsas, thanks a lot for the suggestions. But for now (due to time constraints) I simply note that sundials seems to be more efficient than the Julia exclusive algortihms (at least for this single problem I have at hand). And since I’ve already jumped through hoops to get sundials to integrate the system from my C++ code, I’m not that tempted by re-implementing those tricks on the Julia side (on that note, I do set per-component absolute tolerances in my other code, use an analytic Jacobian, and enforce positivity by using sundials built-in support for constraints as well as clipping in the rhs callback).

And thanks for the offer to look closer at this, but for now I think we can leave it as is. In case you’re interested in the network as a benchmark/test you’re more than welcome to use/include it. Here’s a text version of it:

using DifferentialEquations, DiffEqBiological, Plots
rn = @reaction_network begin
    k1*sA^2, 2sB + 2sA ⇒ sC + 2sD
    k2, 2sE --> sC
    k3, 2sF --> sG
    k4*sE*sA, sE + sB + sA ⇒ sC + sD
    k5, sA + sF --> sD
    k6, sE + sF --> sB
    k7, sA + sG --> sF + sD
    k8, sA + sH --> sI
    k9*sI*sA, 2sB + sI + sA ⇒ sG + 2sD
    k10, sA + sJ --> sK
    k11, sE + sG --> sF + sB
    k12, sE + sH --> sJ
    k13, sE + sJ --> 2sF
    k14, sE + sI --> sK
    k15, sF + sG --> sJ + sB
    k16, sF + sI --> sH + sD
    k17, sF + sJ --> sH + sB
    k18, 2sJ --> sG + sH
    k19*sJ*sI, sB + sJ + sI ⇒ sG + sH + sD
    k20*sG, 2sG ⇒ sH + 2sB
    k21, sF + sK --> sB + sI
    k22, sM + sG --> sD + sJ
    k23, sM + sK --> sD + sI
    k24, sM + sC --> sD + sE
    k25, sB --> sL + sD
    k26, sL + sD --> sB
    k27, sG --> sL + sK
    k28, sL + sK --> sG
    k29, sG + sD --> sK + sB
    k30, sK + sB --> sG + sD
    k31, sF --> sL + sM
    k32, sL + sM --> sF
    k33, sF + sD --> sM + sB
    k34, sM + sB --> sF + sD
    k35, sJ --> sL + sI
    k36, sL + sI --> sJ
    k37, sJ + sD --> sI + sB
    k38, sI + sB --> sJ + sD
    k39, sE --> sL + sA
    k40, sL + sA --> sE
    k41, sE + sD --> sA + sB
    k42, sA + sB --> sE + sD
    k43, sE + sB --> sC + sF
    k44, sC + sF --> sE + sB
    k45, sM + sH --> sN
    k46, sN --> sM + sH
    k47, sS + sD --> sB + sR
    k48, sB + sR --> sS + sD
    k49, sR + sF --> sB + sQ
    k50, sR + sM --> sQ + sD
    k51, sE + sR --> sS + sA
    k52, sS + sA --> sE + sR
    k53, sG + sQ --> sJ + sR
    k54, sK + sQ --> sR + sI
    k55, sJ + sQ --> sR + sH
    k56, sE + sQ --> sR
    k57*sT^2*sH, 2sB + 4sT + sH ⇒ 4sL + 4sV
    k58*sE*sT, sE + 2sT ⇒ sP + sF
    k59, sU + sF --> sL + sW
    k60, sU + sI --> sV + sH
    k61*sU^2, sB + 2sU ⇒ 2sL + sV + sW
    k62, sV + sF --> sU + sD
    k63, sW + sA --> sX
    k64*sX, sB + sX ⇒ sU + 2sD
    k65, sS + sX --> sR + sU + sD
    k66, sL + sX --> sU + sD
    k67, sX + sH --> sW + sI
    k68, sP + sA --> sO + sM
    k69*sV*sA, sB + sV + sA ⇒ sT + 2sD
    k70, sQ + sH --> sB + sT
    k71, sU + sM --> sW
    k72, sT + sI --> sW
    k73*sT*sU, sB + sT + sU ⇒ 2sL + 2sV
    k74, 0 ⇒ sA
    k75, 0 ⇒ sG
    k76, 0 ⇒ sC
    k77, 0 ⇒ sF
    k78, 0 ⇒ sE
    k79, 0 ⇒ sL
    k80, sB ⇒ 0
    k81, 0 ⇒ sH
end k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20 k21 k22 k23 k24 k25 k26 k27 k28 k29 k30 k31 k32 k33 k34 k35 k36 k37 k38 k39 k40 k41 k42 k43 k44 k45 k46 k47 k48 k49 k50 k51 k52 k53 k54 k55 k56 k57 k58 k59 k60 k61 k62 k63 k64 k65 k66 k67 k68 k69 k70 k71 k72 k73 k74 k75 k76 k77 k78 k79 k80 k81
p = Dict([(:k1, 7.226e+09), (:k2, 5.126e+09), (:k3, 4.806e+09), (:k4, 2.754e+10), (:k5, 3.546e+10), (:k6, 1.092e+10), (:k7, 1.357e+10), (:k8, 2.285e+10), (:k9, 1.295e+10), (:k10, 1.295e+10), (:k11, 3.635e+07), (:k12, 1.304e+10), (:k13, 1.135e+10), (:k14, 1.135e+10), (:k15, 2.911e+07), (:k16, 1.095e+10), (:k17, 8.828e+09), (:k18, 8.365e+05), (:k19, 1.01e+08), (:k20, 1.281e-07), (:k21, 4.057e+09), (:k22, 4.057e+09), (:k23, 7.825e+08), (:k24, 1.276e+08), (:k25, 2.086e-05), (:k26, 1.174e+11), (:k27, 0.09327), (:k28, 5.007e+10), (:k29, 1.326e+10), (:k30, 1.264e+06), (:k31, 0.09327), (:k32, 5.007e+10), (:k33, 1.326e+10), (:k34, 1.264e+06), (:k35, 7.692e+05), (:k36, 5.007e+10), (:k37, 1.326e+10), (:k38, 0.1533), (:k39, 5.752), (:k40, 2.092e+10), (:k41, 2.421e+07), (:k42, 15.64), (:k43, 4.481e-05), (:k44, 3.937e+07), (:k45, 3.739e+09), (:k46, 2593), (:k47, 1.326e+10), (:k48, 2.407e+05), (:k49, 7.982e+07), (:k50, 7.982e+06), (:k51, 9.935e+05), (:k52, 6.659e+05), (:k53, 8.974e+07), (:k54, 8.974e+07), (:k55, 8.974e+07),(:k56, 1.994e+10), (:k57, 5.983e+06), (:k58, 1.097e+10), (:k59, 4.487e+09), (:k60, 9.972e+07), (:k61, 9.972e+07), (:k62, 9.672e+08), (:k63, 9.6e+09), (:k64, 9.972e+04), (:k65, 1.994e+08), (:k66, 1.994e+10), (:k67, 1.994e+08), (:k68, 9.6e+09), (:k69, 3.576e+09), (:k70, 1.097e+09), (:k71, 3.49e+09), (:k72, 6.681e+09), (:k73, 7.24e+06), (:k74, 4.266e-08), (:k75, 1.104e-08), (:k76, 6.793e-09), (:k77, 4.353e-08), (:k78, 9.373e-09), (:k79, 4.266e-08), (:k80, 6.561e-08), (:k81, 4.066e-13)])
ics = Dict([(:sB, 55.4), (:sR, 0.1), (:sL, 1e-11), (:sD, 0.001)])
tspan = (0., 3.6e7)
u0 = Array([get(ics, k, 1e-28) for k=keys(speciesmap(rn))])
parr = Array([p[k] for k=keys(paramsmap(rn))])
oprob = ODEProblem(rn, u0, tspan, parr)
sol1 = solve(oprob, reltol=1e-6, abstol=1e-12,
    Rosenbrock23()
)
sol1.retcode, sol1.destats
sol2 = solve(oprob, reltol=1e-7, abstol=1e-12,
    Rodas5()
)
sol2.retcode, sol2.destats
pyplot()
function plot_sol(sol)
    lgt = log10.(clamp.(sol.t[2:end], 0, Inf))
    lgu = log10.(clamp.(hcat(sol.u[2:end]...)', 0, Inf))
    plot(lgt, lgu, fmt = :png, legend=:topleft, legendfontsize=8)
    xlims!((-10, 13))
    ylims!((-20, 2))
end

plot_sol (generic function with 1 method)
plot_sol(sol1)
plot_sol(sol2)
using Sundials
sol3 = solve(oprob, reltol=1e-6, abstol=1e-18,
    CVODE_BDF()
)
sol3.retcode, sol3.destats
plot_sol(sol3)

#oprob4 = ODEProblem(rn, u0, (0, 1e13), parr)
#sol4 = solve(oprob4, reltol=1e-6, abstol=1e-18,
#    CVODE_BDF()
#)
#sol4.retcode, sol4.destats
#plot_sol(sol4)

If I return to trying Julia for solving these kinds of problems (which I am most definitely keen on doing at some point soon), I’ll try the things you’ve suggested here.

All the best, and thanks all for your help and suggestions,
Björn

[ChrisRackauckas]

@bjodah

on that note, I do set per-component absolute tolerances in my other code, use an analytic Jacobian, and enforce positivity by using sundials built-in support for constraints as well as clipping in the rhs callback)

No worries at all. Though I do want to point out that none of these things were done in the Julia code, but all exist in the library. In the documentation, you can find that abstol is allowed to be a vector to do per-component absolute tolerances (this likely would be a massive speed change, since abstol=1e-18 is very small for float64). Analytical sparse Jacobians are created by the reaction network DSL as shown here:

SciML/Catalyst.jl/blob/d496aea0f6367dd0d6150bafd240b2f62d096350/test/misc_tests.jl#L74-L85

Positivity can be enforced by the isoutofdomain argument:

http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html#Miscellaneous-1

or the PositiveDomain callback:

http://docs.juliadiffeq.org/latest/features/callback_library.html#PositiveDomain-1

Indeed, if you already have things working in another language, that’s fine! But yes if you want to see similar speeds, you’ll need to do these kinds of things in Julia as well since it’s giving the solver dramatically more information to specialize on.

[bjodah]

Interesting, curiosity got the best of me so I gave it another shot, I get a stacktrace though (not sure if I’m doing something wrong):

rn = @min_reaction_network begin
    p0*y0^2, 2y1 + 2y0 ⇒ y2 + 2y3
    p1, 2y4 --> y2
    p2, 2y5 --> y6
    p3*y0*y4, y4 + y1 + y0 ⇒ y2 + y3
    p4, y0 + y5 --> y3
    p5, y4 + y5 --> y1
    p6, y0 + y6 --> y5 + y3
    p7, y0 + y7 --> y8
    p8*y0*y8, 2y1 + y8 + y0 ⇒ y6 + 2y3
    p9, y0 + y9 --> y10
    p10, y4 + y6 --> y5 + y1
    p11, y4 + y7 --> y9
    p12, y4 + y9 --> 2y5
    p13, y4 + y8 --> y10
    p14, y5 + y6 --> y9 + y1
    p15, y5 + y8 --> y7 + y3
    p16, y5 + y9 --> y7 + y1
    p17, 2y9 --> y6 + y7
    p18*y9*y8, y1 + y9 + y8 ⇒ y6 + y7 + y3
    p19*y6, 2y6 ⇒ y7 + 2y1
    p20, y5 + y10 --> y1 + y8
    p21, y12 + y6 --> y3 + y9
    p22, y12 + y10 --> y3 + y8
    p23, y12 + y2 --> y3 + y4
    p24, y1 --> y11 + y3
    p25, y11 + y3 --> y1
    p26, y6 --> y11 + y10
    p27, y11 + y10 --> y6
    p28, y6 + y3 --> y10 + y1
    p29, y10 + y1 --> y6 + y3
    p30, y5 --> y11 + y12
    p31, y11 + y12 --> y5
    p32, y5 + y3 --> y12 + y1
    p33, y12 + y1 --> y5 + y3
    p34, y9 --> y11 + y8
    p35, y11 + y8 --> y9
    p36, y9 + y3 --> y8 + y1
    p37, y8 + y1 --> y9 + y3
    p38, y4 --> y11 + y0
    p39, y11 + y0 --> y4
    p40, y4 + y3 --> y0 + y1
    p41, y0 + y1 --> y4 + y3
    p42, y4 + y1 --> y2 + y5
    p43, y2 + y5 --> y4 + y1
    p44, y12 + y7 --> y13
    p45, y13 --> y12 + y7
    p46, y18 + y3 --> y1 + y17
    p47, y1 + y17 --> y18 + y3
    p48, y17 + y5 --> y1 + y16
    p49, y17 + y12 --> y16 + y3
    p50, y4 + y17 --> y18 + y0
    p51, y18 + y0 --> y4 + y17
    p52, y6 + y16 --> y9 + y17
    p53, y10 + y16 --> y17 + y8
    p54, y9 + y16 --> y17 + y7
    p55, y4 + y16 --> y17
    p56*y19^2*y7, 2y1 + 4y19 + y7 ⇒ 4y11 + 4y21
    p57*y19*y4, y4 + 2y19 ⇒ y15 + y5
    p58, y20 + y5 --> y11 + y22
    p59, y20 + y8 --> y21 + y7
    p60*y20^2, y1 + 2y20 ⇒ 2y11 + y21 + y22
    p61, y21 + y5 --> y20 + y3
    p62, y22 + y0 --> y23
    p63*y23, y1 + y23 ⇒ y20 + 2y3
    p64, y18 + y23 --> y17 + y20 + y3
    p65, y11 + y23 --> y20 + y3
    p66, y23 + y7 --> y22 + y8
    p67, y15 + y0 --> y14 + y12
    p68*y0*y21, y1 + y21 + y0 ⇒ y19 + 2y3
    p69, y16 + y7 --> y1 + y19
    p70, y20 + y12 --> y22
    p71, y19 + y8 --> y22
    p72*y20*y19, y1 + y19 + y20 ⇒ 2y11 + 2y21
    p73, 0 ⇒ y0
    p74, 0 ⇒ y6
    p75, 0 ⇒ y2
    p76, 0 ⇒ y5
    p77, 0 ⇒ y4
    p78, 0 ⇒ y11
    p79, y1 ⇒ 0
    p80, 0 ⇒ y7
end p0 p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 p13 p14 p15 p16 p17 p18 p19 p20 p21 p22 p23 p24 p25 p26 p27 p28 p29 p30 p31 p32 p33 p34 p35 p36 p37 p38 p39 p40 p41 p42 p43 p44 p45 p46 p47 p48 p49 p50 p51 p52 p53 p54 p55 p56 p57 p58 p59 p60 p61 p62 p63 p64 p65 p66 p67 p68 p69 p70 p71 p72 p73 p74 p75 p76 p77 p78 p79 p80
addodes!(rn, sparse_jac=true)

It says it’s a boundserror:

BoundsError: attempt to access 24-element Array{Int64,1} at index [25]

Stacktrace:
 [1] getindex at ./array.jl:729 [inlined]
 [2] nzrange at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/SparseArrays/src/sparsematrix.jl:149 [inlined]
 [3] calculate_sparse_jac(::Array{DiffEqBiological.ReactionStruct,1}, ::OrderedCollections.OrderedDict{Symbol,Int64}, ::OrderedCollections.OrderedDict{Symbol,Int64}) at /home/bjorn/.julia/packages/DiffEqBiological/9LwJI/src/reaction_network.jl:651
 [4] #get_jacs#18(::Bool, ::Bool, ::Bool, ::Bool, ::Function, ::Array{Union{Float64, Int64, Expr, Symbol},1}, ::Array{Symbol,1}, ::Array{DiffEqBiological.ReactionStruct,1}, ::OrderedCollections.OrderedDict{Symbol,Int64}, ::OrderedCollections.OrderedDict{Symbol,Int64}) at /home/bjorn/.julia/packages/DiffEqBiological/9LwJI/src/reaction_network.jl:452
 [5] #get_jacs at /home/bjorn/.julia/packages/DiffEqBiological/9LwJI/src/reaction_network.jl:0 [inlined]
 [6] #genode_exprs#11(::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Function, ::Array{DiffEqBiological.ReactionStruct,1}, ::OrderedCollections.OrderedDict{Symbol,Int64}, ::OrderedCollections.OrderedDict{Symbol,Int64}, ::Array{Symbol,1}) at /home/bjorn/.julia/packages/DiffEqBiological/9LwJI/src/reaction_network.jl:191
 [7] #genode_exprs at ./none:0 [inlined]
 [8] #addodes!#30(::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:sparse_jac,),Tuple{Bool}}}, ::Function, ::min_reaction_network) at /home/bjorn/.julia/packages/DiffEqBiological/9LwJI/src/maketype.jl:343
 [9] (::getfield(DiffEqBiological, Symbol("#kw##addodes!")))(::NamedTuple{(:sparse_jac,),Tuple{Bool}}, ::typeof(addodes!), ::min_reaction_network) at ./none:0
 [10] top-level scope at In[2]:84

Using @reaction_network does work.

I’m experimenting with either passing:
callback=PositiveDomain(ones(length(u0)), abstol=1e19.*abstol)
or
isoutofdomain=(u,p,t) -> any(x -> x < 0, u)
to solve(), but so far they haven’t proven to be the silver bullet for this problem.

[isaacsas]

@bjodah

rn = @min_reaction_network begin p0y0^2, 2y1 + 2y0 ⇒ y2 + 2y3 p1, 2y4 → y2 p2, 2y5 → y6 p3y0y4, y4 + y1 + y0 ⇒ y2 + y3 p4, y0 + y5 → y3 p5, y4 + y5 → y1 p6, y0 + y6 → y5 + y3 p7, y0 + y7 → y8 p8y0y8, 2y1 + y8 + y0 ⇒ y6 + 2y3 p9, y0 + y9 → y10 p10, y4 + y6 → y5 + y1 p11, y4 + y7 → y9 p12, y4 + y9 → 2y5 p13, y4 + y8 → y10 p14, y5 + y6 → y9 + y1 p15, y5 + y8 → y7 + y3 p16, y5 + y9 → y7 + y1 p17, 2y9 → y6 + y7 p18y9y8, y1 + y9 + y8 ⇒ y6 + y7 + y3 p19y6, 2y6 ⇒ y7 + 2y1 p20, y5 + y10 → y1 + y8 p21, y12 + y6 → y3 + y9 p22, y12 + y10 → y3 + y8 p23, y12 + y2 → y3 + y4 p24, y1 → y11 + y3 p25, y11 + y3 → y1 p26, y6 → y11 + y10 p27, y11 + y10 → y6 p28, y6 + y3 → y10 + y1 p29, y10 + y1 → y6 + y3 p30, y5 → y11 + y12 p31, y11 + y12 → y5 p32, y5 + y3 → y12 + y1 p33, y12 + y1 → y5 + y3 p34, y9 → y11 + y8 p35, y11 + y8 → y9 p36, y9 + y3 → y8 + y1 p37, y8 + y1 → y9 + y3 p38, y4 → y11 + y0 p39, y11 + y0 → y4 p40, y4 + y3 → y0 + y1 p41, y0 + y1 → y4 + y3 p42, y4 + y1 → y2 + y5 p43, y2 + y5 → y4 + y1 p44, y12 + y7 → y13 p45, y13 → y12 + y7 p46, y18 + y3 → y1 + y17 p47, y1 + y17 → y18 + y3 p48, y17 + y5 → y1 + y16 p49, y17 + y12 → y16 + y3 p50, y4 + y17 → y18 + y0 p51, y18 + y0 → y4 + y17 p52, y6 + y16 → y9 + y17 p53, y10 + y16 → y17 + y8 p54, y9 + y16 → y17 + y7 p55, y4 + y16 → y17 p56y19^2y7, 2y1 + 4y19 + y7 ⇒ 4y11 + 4y21 p57y19y4, y4 + 2y19 ⇒ y15 + y5 p58, y20 + y5 → y11 + y22 p59, y20 + y8 → y21 + y7 p60y20^2, y1 + 2y20 ⇒ 2y11 + y21 + y22 p61, y21 + y5 → y20 + y3 p62, y22 + y0 → y23 p63y23, y1 + y23 ⇒ y20 + 2y3 p64, y18 + y23 → y17 + y20 + y3 p65, y11 + y23 → y20 + y3 p66, y23 + y7 → y22 + y8 p67, y15 + y0 → y14 + y12 p68y0y21, y1 + y21 + y0 ⇒ y19 + 2y3 p69, y16 + y7 → y1 + y19 p70, y20 + y12 → y22 p71, y19 + y8 → y22 p72y20y19, y1 + y19 + y20 ⇒ 2y11 + 2y21 p73, 0 ⇒ y0 p74, 0 ⇒ y6 p75, 0 ⇒ y2 p76, 0 ⇒ y5 p77, 0 ⇒ y4 p78, 0 ⇒ y11 p79, y1 ⇒ 0 p80, 0 ⇒ y7 end p0 p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 p13 p14 p15 p16 p17 p18 p19 p20 p21 p22 p23 p24 p25 p26 p27 p28 p29 p30 p31 p32 p33 p34 p35 p36 p37 p38 p39 p40 p41 p42 p43 p44 p45 p46 p47 p48 p49 p50 p51 p52 p53 p54 p55 p56 p57 p58 p59 p60 p61 p62 p63 p64 p65 p66 p67 p68 p69 p70 p71 p72 p73 p74 p75 p76 p77 p78 p79 p80 addodes!(rn, sparse_jac=true)

I’ll look into that. The sparse Jacobian stuff has so far only been used by me I think (it only went into master a few months ago, and was only released recently).

I’m not sure if it would really be helpful for you anyways. If you only have 24 species your Jacobian should 24x24, which would probably be better as a dense matrix…

[bjodah]

Indeed 24**3 floating point operations is nothing, so a dense LU is the way to go, and as you point out: it’s not even particularly sparse.

[isaacsas]

In the playing around I did with a B cell signaling network with ~1000 species, sparse was still noticeably slower than dense.

[ChrisRackauckas]

In those circumstances, we probably want to use a different factorization than SuperLUMT.

[isaacsas]

@bjodah Ok, I think sparse isn’t working for you because the last column is all zeros, and the current code base didn’t explicitly set the sparse matrix size (i.e. it implicitly assumed there was at least one nonzero value in each column/row). I’ll put in a fix for that.

More generally, should your network have a singular Jacobian? (I’m seeing this in the dense Jacobian with your earlier network you specified too.) I guess this is arising from the sO species, which is not involved in any rate law.

[isaacsas]

What’s the recommended one to use in Julia then? I thought the default in SparseArrays is a SuiteSparse solver?

[ChrisRackauckas]

In most cases, SuperLUMT, and that is the default. However, what I am saying is that for networks of that size, there has been studies that show that KLU works better on DAEs derived from circuits. So we might possibly want to specialize better. Sundials has an option to build with KLU for this reason.