feat: add BCR jacobian benchmark

Author: AayushSabharwal

Diff for: benchmarks/Symbolics/BCR.jmd

@@ -0,0 +1,60 @@
+---
+title: BCR Symbolic Jacobian
+author: Aayush Sabharwal, Bowen Zhu, Chris Rackauckas
+---
+
+The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff
+chemical reaction network modeling the BCR signaling network from [Barua et
+al.](https://doi.org/10.4049/jimmunol.1102003). 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.
+
+The resultant large model is used to benchmark the time taken to compute a symbolic
+jacobian, generate a function to calculate it and call the function.
+
+
+```julia
+using OrdinaryDiffEq, Catalyst, ReactionNetworkImporters,
+    Plots, TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools,
+    LinearSolve, Symbolics, SymbolicUtils.Code, SparseArrays
+
+datadir  = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr")
+const to = TimerOutput()
+tf       = 100000.0
+
+# generate ModelingToolkit ODEs
+prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net"))
+show(to)
+rn    = complete(prnbng.rn; split = false)
+obs = [eq.lhs for eq in observed(rn)]
+osys = convert(ODESystem, rn)
+
+rhs = [eq.rhs for eq in full_equations(osys)]
+vars = unknowns(osys)
+pars = parameters(osys)
+@timeit to "Calculate jacobian" jac = Symbolics.jacobian(rhs, vars);
+jac = sparse(jac)
+
+args = (vars, pars, ModelingToolkit.get_iv(osys))
+# out of place versions run into an error saying the expression is too large
+# due to the `SymbolicUtils.Code.create_array` call.
+@timeit to "Build jacobian - no CSE" _, jac_nocse_iip = build_function(jac, args...; cse = false, expression = Val{false});
+@timeit to "Build jacobian - CSE" _, jac_cse_iip = build_function(jac, args...; cse = true, expression = Val{false});
+
+defs = defaults(osys)
+u = Float64[Symbolics.fixpoint_sub(var, defs) for var in vars]
+buffer = SparseMatrixCSC{Float64, Int}(undef, length(u), length(u))
+p = Float64[Symbolics.fixpoint_sub(par, defs) for par in pars]
+tt = 0.0
+
+@timeit to "Compile jacobian - CSE" jac_cse_iip(buffer, u, p, tt)
+@timeit to "Compute jacobian - CSE" jac_cse_iip(buffer, u, p, t)
+
+@timeit to "Compile jacobian - no CSE" jac_nocse_iip(buffer, u, p, tt)
+@timeit to "Compute jacobian - no CSE" jac_nocse_iip(buffer, u, p, t)
+
+show(to)
+```

Diff for: benchmarks/Symbolics/Project.toml

@@ -0,0 +1,15 @@
+[deps]
+Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+
+[compat]
+SymbolicUtils = "3.15"
+Symbolics = "6.29.1"
+ModelingToolkit = "9"