@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