|
| 1 | +--- |
| 2 | +title: BCR Symbolic Jacobian |
| 3 | +author: Aayush Sabharwal, Bowen Zhu, Chris Rackauckas |
| 4 | +--- |
| 5 | + |
| 6 | +The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff |
| 7 | +chemical reaction network modeling the BCR signaling network from [Barua et |
| 8 | +al.](https://doi.org/10.4049/jimmunol.1102003). We use |
| 9 | +[`ReactionNetworkImporters`](https://github.com/isaacsas/ReactionNetworkImporters.jl) |
| 10 | +to load the BioNetGen model files as a |
| 11 | +[Catalyst](https://github.com/SciML/Catalyst.jl) model, and then use |
| 12 | +[ModelingToolkit](https://github.com/SciML/ModelingToolkit.jl) to convert the |
| 13 | +Catalyst network model to ODEs. |
| 14 | + |
| 15 | +The resultant large model is used to benchmark the time taken to compute a symbolic |
| 16 | +jacobian, generate a function to calculate it and call the function. |
| 17 | + |
| 18 | + |
| 19 | +```julia |
| 20 | +using OrdinaryDiffEq, Catalyst, ReactionNetworkImporters, |
| 21 | + Plots, TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools, |
| 22 | + LinearSolve, Symbolics, SymbolicUtils.Code, SparseArrays |
| 23 | + |
| 24 | +datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") |
| 25 | +const to = TimerOutput() |
| 26 | +tf = 100000.0 |
| 27 | + |
| 28 | +# generate ModelingToolkit ODEs |
| 29 | +prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net")) |
| 30 | +show(to) |
| 31 | +rn = complete(prnbng.rn; split = false) |
| 32 | +obs = [eq.lhs for eq in observed(rn)] |
| 33 | +osys = convert(ODESystem, rn) |
| 34 | + |
| 35 | +rhs = [eq.rhs for eq in full_equations(osys)] |
| 36 | +vars = unknowns(osys) |
| 37 | +pars = parameters(osys) |
| 38 | +@timeit to "Calculate jacobian" jac = Symbolics.jacobian(rhs, vars); |
| 39 | +jac = sparse(jac) |
| 40 | + |
| 41 | +args = (vars, pars, ModelingToolkit.get_iv(osys)) |
| 42 | +# out of place versions run into an error saying the expression is too large |
| 43 | +# due to the `SymbolicUtils.Code.create_array` call. |
| 44 | +@timeit to "Build jacobian - no CSE" _, jac_nocse_iip = build_function(jac, args...; cse = false, expression = Val{false}); |
| 45 | +@timeit to "Build jacobian - CSE" _, jac_cse_iip = build_function(jac, args...; cse = true, expression = Val{false}); |
| 46 | + |
| 47 | +defs = defaults(osys) |
| 48 | +u = Float64[Symbolics.fixpoint_sub(var, defs) for var in vars] |
| 49 | +buffer = SparseMatrixCSC{Float64, Int}(undef, length(u), length(u)) |
| 50 | +p = Float64[Symbolics.fixpoint_sub(par, defs) for par in pars] |
| 51 | +tt = 0.0 |
| 52 | + |
| 53 | +@timeit to "Compile jacobian - CSE" jac_cse_iip(buffer, u, p, tt) |
| 54 | +@timeit to "Compute jacobian - CSE" jac_cse_iip(buffer, u, p, t) |
| 55 | + |
| 56 | +@timeit to "Compile jacobian - no CSE" jac_nocse_iip(buffer, u, p, tt) |
| 57 | +@timeit to "Compute jacobian - no CSE" jac_nocse_iip(buffer, u, p, t) |
| 58 | + |
| 59 | +show(to) |
| 60 | +``` |
0 commit comments