egfr_net.jmd

---
title: Egfr_net Work-Precision Diagrams
author: Torkel Loman
---

The following benchmark is of 356 ODEs with 3749 terms that describe a
chemical reaction network. This egfr_net model was used as a benchmark model in [Gupta et
al.](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6013266/). It describes the epidermal growth factor receptor signalling system [Blinov et
al.](https://pubmed.ncbi.nlm.nih.gov/16233948/). We use
[`ReactionNetworkImporters`](https://github.com/isaacsas/ReactionNetworkImporters.jl)
to load the BioNetGen model files as a
[Catalyst](https://github.com/SciML/Catalyst.jl) model, and then use
[ModelingToolkit](https://github.com/SciML/ModelingToolkit.jl) to convert the
Catalyst network model to ODEs.


```julia
using DiffEqBase, OrdinaryDiffEq, Catalyst, ReactionNetworkImporters,
      Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq,
      LSODA, TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools,
      LinearSolve, RecursiveFactorization

gr()
const to = TimerOutput()
tf       = 10.0

# generate ModelingToolkit ODEs
@timeit to "Parse Network" prnbng = loadrxnetwork(BNGNetwork(), joinpath(@__DIR__, "Models/egfr_net.net"))
show(to)
rn    = complete(prnbng.rn)
obs = [eq.lhs for eq in observed(rn)]

@timeit to "Create ODESys" osys = complete(convert(ODESystem, rn))
show(to)

tspan = (0.,tf)
@timeit to "ODEProb No Jac" oprob = ODEProblem{true, SciMLBase.FullSpecialize}(osys, Float64[], tspan, Float64[])
show(to);
```

```julia
@timeit to "ODEProb SparseJac" sparsejacprob = ODEProblem{true, SciMLBase.FullSpecialize}(osys, Float64[], tspan, Float64[], jac=true, sparse=true)
show(to)
```


```julia
@show numspecies(rn) # Number of ODEs
@show numreactions(rn) # Apprx. number of terms in the ODE
@show length(parameters(rn)); # Number of Parameters
```

## Time ODE derivative function compilation
As compiling the ODE derivative functions has in the past taken longer than
running a simulation, we first force compilation by evaluating these functions
one time.
```julia
u  = oprob.u0
du = copy(u)
p  = oprob.p
@timeit to "ODE rhs Eval1" oprob.f(du,u,p,0.)
@timeit to "ODE rhs Eval2" oprob.f(du,u,p,0.)
sparsejacprob.f(du,u,p,0.)
```

We also time the ODE rhs function with BenchmarkTools as it is more accurate
given how fast evaluating `f` is:
```julia
@btime oprob.f($du,$u,$p,0.)
```

## Picture of the solution

```julia
sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5)
plot(sol; idxs=obs, legend=false, fmt=:png)
```

For these benchmarks we will be using the time-series error with these saving
points.

## Generate Test Solution

```julia
@time sol = solve(oprob, CVODE_BDF(), abstol=1/10^14, reltol=1/10^14)
test_sol  = TestSolution(sol);
```
## Setups

#### Sets plotting defaults

```julia
default(legendfontsize=7,framestyle=:box,gridalpha=0.3,gridlinewidth=2.5)
```

#### Sets tolerances

```julia
abstols = 1.0 ./ 10.0 .^ (6:10)
reltols = 1.0 ./ 10.0 .^ (6:10);
```

## Implicit Work-Precision Diagrams

Benchmarks for implicit solvers.

#### Declare solvers (using default linear solver)

We designate the solvers we wish to use.
```julia
setups = [
          Dict(:alg=>lsoda()),
          Dict(:alg=>CVODE_BDF()),
          Dict(:alg=>CVODE_BDF(linear_solver=:LapackDense)),
          Dict(:alg=>CVODE_Adams()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>QNDF()),
          Dict(:alg=>FBDF()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>Rodas5P())
          ];
```


#### Plot Work-Precision Diagram (using default linear solver)

Finally, we generate a work-precision diagram for the selection of solvers.
```julia
wp = WorkPrecisionSet(oprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e9),numruns=100)

names = ["lsoda" "CVODE_BDF" "CVODE_BDF (Lapack Dense)" "CVODE_Adams" "TRBDF2" "QNDF" "FBDF" "KenCarp4" "Rosenbrock23" "Rodas4" "Rodas5P"]
plot(wp;label=names)
```

#### Declare solvers (using GMRES linear solver)

We designate the solvers we wish to use.
```julia
setups = [
          Dict(:alg=>lsoda()),
          Dict(:alg=>CVODE_BDF(linear_solver=:GMRES)),
          Dict(:alg=>TRBDF2(linsolve=KrylovJL_GMRES())),
          Dict(:alg=>QNDF(linsolve=KrylovJL_GMRES())),
          Dict(:alg=>FBDF(linsolve=KrylovJL_GMRES())),
          Dict(:alg=>KenCarp4(linsolve=KrylovJL_GMRES())),
          Dict(:alg=>Rosenbrock23(linsolve=KrylovJL_GMRES())),
          Dict(:alg=>Rodas4(linsolve=KrylovJL_GMRES())),
          Dict(:alg=>Rodas5P(linsolve=KrylovJL_GMRES()))
          ];
```


#### Plot Work-Precision Diagram (using GMRES linear solver)

Finally, we generate a work-precision diagram for the selection of solvers.
```julia
wp = WorkPrecisionSet(oprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e9),numruns=100)

names = ["lsoda" "CVODE_BDF (GMRES)" "TRBDF2 (GMRES)" "QNDF (GMRES)" "FBDF (GMRES)" "KenCarp4 (GMRES)" "Rosenbrock23 (GMRES)" "Rodas4 (GMRES)" "Rodas5P (GMRES)"]
plot(wp;label=names)
```

#### Declare solvers (using sparse jacobian)

We designate the solvers we wish to use.
```julia
setups = [
          Dict(:alg=>CVODE_BDF(linear_solver=:KLU)),
          Dict(:alg=>TRBDF2(linsolve=KLUFactorization())),
          Dict(:alg=>QNDF(linsolve=KLUFactorization())),
          Dict(:alg=>FBDF(linsolve=KLUFactorization())),
          Dict(:alg=>KenCarp4(linsolve=KLUFactorization())),
          Dict(:alg=>Rosenbrock23(linsolve=KLUFactorization())),
          Dict(:alg=>Rodas4(linsolve=KLUFactorization())),
          Dict(:alg=>Rodas5P(linsolve=KLUFactorization()))
          ];
```


#### Plot Work-Precision Diagram (using sparse jacobian)

Finally, we generate a work-precision diagram for the selection of solvers.
```julia
wp = WorkPrecisionSet(sparsejacprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e9),numruns=100)

names = ["CVODE_BDF (KLU, sparse jac)" "TRBDF2 (KLU, sparse jac)" "QNDF (KLU, sparse jac)" "FBDF (KLU, sparse jac)" "KenCarp4 (KLU, sparse jac)" "Rosenbrock23 (KLU, sparse jac)" "Rodas4 (KLU, sparse jac)" "Rodas5P (KLU, sparse jac)"]
plot(wp;label=names)
```

## Explicit Work-Precision Diagram

Benchmarks for explicit solvers.

#### Declare solvers

We designate the solvers we wish to use, this also includes lsoda and CVODE.
```julia
setups = [
          Dict(:alg=>lsoda()),
          Dict(:alg=>CVODE_Adams()),
          Dict(:alg=>Tsit5()),
          Dict(:alg=>BS5()),
          Dict(:alg=>VCABM()),
          Dict(:alg=>Vern6()),
          Dict(:alg=>Vern7()),
          Dict(:alg=>Vern8()),
          Dict(:alg=>Vern9()),
          Dict(:alg=>ROCK4())
          ];
```

#### Plot Work-Precision Diagram

```julia
wp = WorkPrecisionSet(oprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e9),numruns=200)

names = ["lsoda" "CVODE_Adams" "Tsit5" "BS5" "VCABM" "Vern6" "Vern7" "Vern8" "Vern9" "ROCK4"]
plot(wp;label=names)
```

#### Additional explicit solvers
One additional explicit solver, `ROCK2`, performs noticeably worse as compared to the other ones.
```julia
setups = [Dict(:alg=>ROCK2())];
wp = WorkPrecisionSet(oprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e9),numruns=200)
names = ["ROCK2"]
plot(wp;label=names)
```

## Summary of results
Finally, we compute a single diagram comparing the various solvers used.

#### Declare solvers
We designate the solvers we wish to compare.
```julia
setups = [
          Dict(:alg=>lsoda(), :prob_choice => 1),
          Dict(:alg=>CVODE_BDF(linear_solver=:GMRES), :prob_choice => 1),
          Dict(:alg=>QNDF(linsolve=KrylovJL_GMRES()), :prob_choice => 2),
          Dict(:alg=>QNDF(linsolve=KLUFactorization()), :prob_choice => 2),
          Dict(:alg=>BS5(), :prob_choice => 1),
          Dict(:alg=>Vern6(), :prob_choice => 1),
          Dict(:alg=>ROCK4(), :prob_choice => 1)
          ];
```

#### Plot Work-Precision Diagram

For these, we generate a work-precision diagram for the selection of solvers.
```julia
wp = WorkPrecisionSet([oprob,sparsejacprob],abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=[test_sol,test_sol],maxiters=Int(1e9),numruns=200)

names = ["lsoda" "CVODE_BDF (GMRES)" "QNDF (GMRES)" "QNDF (KLU)" "BS5" "Vern6" "ROCK4"]
colors = [:seagreen1 :darkgreen :deepskyblue1 :deepskyblue4 :thistle2 :lightslateblue :purple4]
markershapes = [:star4 :rect :hexagon :octagon :star8 :rtriangle :square]
plot(wp;label=names,left_margin=10Plots.mm,right_margin=10Plots.mm,xticks=[1e-9,1e-8,1e-7,1e-6,1e-5,1e-4,1e-3,1e-2],yticks=[1e-2,1e-1],color=colors,markershape=markershapes,legendfontsize=15,tickfontsize=15,guidefontsize=15, legend=:topright, lw=20, la=0.8, markersize=20,markerstrokealpha=1.0, markerstrokewidth=1.5, gridalpha=0.3, gridlinewidth=7.5,size=(1100,1000))
```


```julia, echo = false
using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
```