From 2ac29c8dc674976729b9496e03696bfe5bb0a4a5 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 20 Feb 2025 18:08:22 +0530 Subject: [PATCH 01/31] feat: add BCR jacobian benchmark --- benchmarks/Symbolics/BCR.jmd | 77 +++++++++++++++++++++++++++++++ benchmarks/Symbolics/Project.toml | 17 +++++++ 2 files changed, 94 insertions(+) create mode 100644 benchmarks/Symbolics/BCR.jmd create mode 100644 benchmarks/Symbolics/Project.toml diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd new file mode 100644 index 000000000..6a45625af --- /dev/null +++ b/benchmarks/Symbolics/BCR.jmd @@ -0,0 +1,77 @@ +--- +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) + +SymbolicUtils.ENABLE_HASHCONSING[] = false +@timeit to "Calculate jacobian - without hashconsing" jac_nohc = Symbolics.sparsejacobian(rhs, vars); +SymbolicUtils.ENABLE_HASHCONSING[] = true +SymbolicUtils.toggle_caching!(Symbolics.occursin_info, false) +@timeit to "Calculate jacobian - hashconsing, without caching" jac_hc_nocache = Symbolics.sparsejacobian(rhs, vars); +SymbolicUtils.toggle_caching!(Symbolics.occursin_info, true) +stats = SymbolicUtils.get_stats(Symbolics.occursin_info) +@assert stats.hits == stats.misses == 0 +@timeit to "Calculate jacobian - hashconsing and caching" jac_hc_cache = Symbolics.sparsejacobian(rhs, vars); + +@assert isequal(jac_nohc, jac_hc_nocache) +@assert isequal(jac_hc_nocache, jac_hc_cache) + +jac = jac_hc_cache +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. `iip_config` prevents it +# from trying to build the function. +kwargs = (; iip_config = (false, true), expression = Val{false}) +@timeit to "Build jacobian - no CSE" _, jac_nocse_iip = build_function(jac, args...; cse = false, kwargs...); +@timeit to "Build jacobian - CSE" _, jac_cse_iip = build_function(jac, args...; cse = true, kwargs...); + +defs = defaults(osys) +u = Float64[Symbolics.fixpoint_sub(var, defs) for var in vars] +buffer_cse = SparseMatrixCSC{Float64, Int}(undef, length(u), length(u)) +buffer_nocse = 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_cse, u, p, tt) +@timeit to "Compute jacobian - CSE" jac_cse_iip(buffer_cse, u, p, t) + +@timeit to "Compile jacobian - no CSE" jac_nocse_iip(buffer_nocse, u, p, tt) +@timeit to "Compute jacobian - no CSE" jac_nocse_iip(buffer_nocse, u, p, t) + +@assert isapprox(buffer_cse, buffer_nocse, rtol = 1e-10) + +show(to) +``` diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml new file mode 100644 index 000000000..61d84b08b --- /dev/null +++ b/benchmarks/Symbolics/Project.toml @@ -0,0 +1,17 @@ +[deps] +Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +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] +ModelingToolkit = "9" +SymbolicUtils = "3.17" +Symbolics = "6.29.1" From 0e74194c48801e08e74e6d78c88f381d533e6149 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 13 Mar 2025 23:59:10 +0530 Subject: [PATCH 02/31] fix: fix preallocation of buffers in BCR benchmark --- benchmarks/Symbolics/BCR.jmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 6a45625af..9bea261a5 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -60,8 +60,8 @@ kwargs = (; iip_config = (false, true), expression = Val{false}) defs = defaults(osys) u = Float64[Symbolics.fixpoint_sub(var, defs) for var in vars] -buffer_cse = SparseMatrixCSC{Float64, Int}(undef, length(u), length(u)) -buffer_nocse = SparseMatrixCSC{Float64, Int}(undef, length(u), length(u)) +buffer_cse = similar(Float64, jac) +buffer_nocse = similar(Float64, jac) p = Float64[Symbolics.fixpoint_sub(par, defs) for par in pars] tt = 0.0 From ac0c0a12ae2819d334fb7005b208c79ae96b7e8c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 13 Mar 2025 23:59:21 +0530 Subject: [PATCH 03/31] feat: add ThermalFluid jacobian scaling benchmark --- benchmarks/Symbolics/Project.toml | 5 + benchmarks/Symbolics/ThermalFluid.jmd | 309 ++++++++++++++++++++++++++ 2 files changed, 314 insertions(+) create mode 100644 benchmarks/Symbolics/ThermalFluid.jmd diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml index 61d84b08b..1ccb42a57 100644 --- a/benchmarks/Symbolics/Project.toml +++ b/benchmarks/Symbolics/Project.toml @@ -1,4 +1,5 @@ [deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" @@ -6,10 +7,14 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881" +SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +XSteam = "95ff35a0-be81-11e9-2ca3-5b4e338e8476" [compat] ModelingToolkit = "9" diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd new file mode 100644 index 000000000..309dd93b4 --- /dev/null +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -0,0 +1,309 @@ +--- +title: Thermal Fluid Symbolic Jacobian Scaling +author: Aayush Sabharwal, Bowen Zhu, Chris Rackauckas +--- + +This is a 1D advection-diffusion-source PDE that uses a second order upwind scheme. + +```julia +using Pkg +# Rev fixes precompilation https://github.com/hzgzh/XSteam.jl/pull/2 +Pkg.add(Pkg.PackageSpec(;name="XSteam", rev="f2a1c589054cfd6bba307985a3a534b6f5a1863b")) + +using ModelingToolkit, Symbolics, XSteam, Polynomials, CairoMakie +using SparseArrays +using ModelingToolkit: t_nounits as t, D_nounits as D +``` + +## Setup Julia Code + +```julia +# o o o o o o o < heat capacitors +# | | | | | | | < heat conductors +# o o o o o o o +# | | | | | | | +#Source -> o--o--o--o--o--o--o -> Sink +# advection diff source PDE + +m_flow_source(t) = 2.75 +T_source(t) = (t > 12 * 3600) * 56.0 + 12.0 +@register_symbolic m_flow_source(t) +@register_symbolic T_source(t) + +#build polynomial liquid-water property only dependent on Temperature +p_l = 5 #bar +T_vec = collect(1:1:150); +@generated kin_visc_T(t) = :(Base.evalpoly(t, $(fit(T_vec, my_pT.(p_l, T_vec) ./ rho_pT.(p_l, T_vec), 5).coeffs...,))) +@generated lambda_T(t) = :(Base.evalpoly(t, $(fit(T_vec, tc_pT.(p_l, T_vec), 3).coeffs...,))) +@generated Pr_T(t) = :(Base.evalpoly(t, $(fit(T_vec, 1e3 * Cp_pT.(p_l, T_vec) .* my_pT.(p_l, T_vec) ./ tc_pT.(p_l, T_vec), 5).coeffs...,))) +@generated rho_T(t) = :(Base.evalpoly(t, $(fit(T_vec, rho_pT.(p_l, T_vec), 4).coeffs...,))) +@generated rhocp_T(t) = :(Base.evalpoly(t, $(fit(T_vec, 1000 * rho_pT.(p_l, T_vec) .* Cp_pT.(p_l, T_vec), 5).coeffs...,))) +@register_symbolic kin_visc_T(t) +@register_symbolic lambda_T(t) +@register_symbolic Pr_T(t) +@register_symbolic rho_T(t) +@register_symbolic rhocp_T(t) + +@connector function FluidPort(; name, p=101325.0, m=0.0, T=0.0) + sts = @variables p(t) = p m(t) = m [connect = Flow] T(t) = T [connect = Stream] + ODESystem(Equation[], t, sts, []; name=name) +end + +@connector function VectorHeatPort(; name, N=100, T0=0.0, Q0=0.0) + sts = @variables (T(t))[1:N] = T0 (Q(t))[1:N] = Q0 [connect = Flow] + ODESystem(Equation[], t, [T; Q], []; name=name) +end + +@register_symbolic Dxx_coeff(u, d, T) +#Taylor-aris dispersion model +function Dxx_coeff(u, d, T) + Re = abs(u) * d / kin_visc_T(T) + 0.1 + if Re < 1000.0 + (d^2 / 4) * u^2 / 48 / 0.14e-6 + else + d * u * (1.17e9 * Re^(-2.5) + 0.41) + end +end + +@register_symbolic Nusselt(Re, Pr, f) +#Nusselt number model +function Nusselt(Re, Pr, f) + if Re <= 2300.0 + 3.66 + elseif Re <= 3100.0 + 3.5239 * (Re / 1000)^4 - 45.158 * (Re / 1000)^3 + 212.13 * (Re / 1000)^2 - 427.45 * (Re / 1000) + 316.08 + else + f / 8 * ((Re - 1000) * Pr) / (1 + 12.7 * (f / 8)^(1 / 2) * (Pr^(2 / 3) - 1)) + end +end + +@register_symbolic Churchill_f(Re, epsilon, d) +#Darcy weisbach friction factor +function Churchill_f(Re, epsilon, d) + theta_1 = (-2.457 * log(((7 / Re)^0.9) + (0.27 * (epsilon / d))))^16 + theta_2 = (37530 / Re)^16 + 8 * ((((8 / Re)^12) + (1 / ((theta_1 + theta_2)^1.5)))^(1 / 12)) +end + +function FluidRegion(; name, L=1.0, dn=0.05, N=100, T0=0.0, + lumped_T=50, diffusion=true, e=1e-4) + @named inlet = FluidPort() + @named outlet = FluidPort() + @named heatport = VectorHeatPort(N=N) + + dx = L / N + c = [-1 / 8, -3 / 8, -3 / 8] # advection stencil coefficients + A = pi * dn^2 / 4 + + p = @parameters C_shift = 0.0 Rw = 0.0 # stuff for latter + @variables begin + (T(t))[1:N] = fill(T0, N) + Twall(t)[1:N] = fill(T0, N) + (S(t))[1:N] = fill(T0, N) + (C(t))[1:N] = fill(1.0, N) + u(t) = 1e-6 + Re(t) = 1000.0 + Dxx(t) = 0.0 + Pr(t) = 1.0 + alpha(t) = 1.0 + f(t) = 1.0 + end + + sts = vcat(T, Twall, S, C, Num[u], Num[Re], Num[Dxx], Num[Pr], Num[alpha], Num[f]) + + eqs = Equation[ + Re ~ 0.1 + dn * abs(u) / kin_visc_T(lumped_T) + Pr ~ Pr_T(lumped_T) + f ~ Churchill_f(Re, e, dn) #Darcy-weisbach + alpha ~ Nusselt(Re, Pr, f) * lambda_T(lumped_T) / dn + Dxx ~ diffusion * Dxx_coeff(u, dn, lumped_T) + inlet.m ~ -outlet.m + inlet.p ~ outlet.p + inlet.T ~ instream(inlet.T) + outlet.T ~ T[N] + u ~ inlet.m / rho_T(inlet.T) / A + [C[i] ~ dx * A * rhocp_T(T[i]) for i in 1:N] + [S[i] ~ heatport.Q[i] for i in 1:N] + [Twall[i] ~ heatport.T[i] for i in 1:N] + + #source term + [S[i] ~ (1 / (1 / (alpha * dn * pi * dx) + abs(Rw / 1000))) * (Twall[i] - T[i]) for i in 1:N] + + #second order upwind + diffusion + source + D(T[1]) ~ u / dx * (inlet.T - T[1]) + Dxx * (T[2] - T[1]) / dx^2 + S[1] / (C[1] - C_shift) + D(T[2]) ~ u / dx * (c[1] * inlet.T - sum(c) * T[1] + c[2] * T[2] + c[3] * T[3]) + Dxx * (T[1] - 2 * T[2] + T[3]) / dx^2 + S[2] / (C[2] - C_shift) + [D(T[i]) ~ u / dx * (c[1] * T[i-2] - sum(c) * T[i-1] + c[2] * T[i] + c[3] * T[i+1]) + Dxx * (T[i-1] - 2 * T[i] + T[i+1]) / dx^2 + S[i] / (C[i] - C_shift) for i in 3:N-1] + D(T[N]) ~ u / dx * (T[N-1] - T[N]) + Dxx * (T[N-1] - T[N]) / dx^2 + S[N] / (C[N] - C_shift) + ] + + ODESystem(eqs, t, sts, p; systems=[inlet, outlet, heatport], name=name) +end + +@register_symbolic Cn_circular_wall_inner(d, D, cp, ρ) +function Cn_circular_wall_inner(d, D, cp, ρ) + C = pi / 4 * (D^2 - d^2) * cp * ρ + return C / 2 +end + +@register_symbolic Cn_circular_wall_outer(d, D, cp, ρ) +function Cn_circular_wall_outer(d, D, cp, ρ) + C = pi / 4 * (D^2 - d^2) * cp * ρ + return C / 2 +end + +@register_symbolic Ke_circular_wall(d, D, λ) +function Ke_circular_wall(d, D, λ) + 2 * pi * λ / log(D / d) +end + +function CircularWallFEM(; name, L=100, N=10, d=0.05, t_layer=[0.002], + λ=[50], cp=[500], ρ=[7850], T0=0.0) + @named inner_heatport = VectorHeatPort(N=N) + @named outer_heatport = VectorHeatPort(N=N) + dx = L / N + Ne = length(t_layer) + Nn = Ne + 1 + dn = vcat(d, d .+ 2.0 .* cumsum(t_layer)) + Cn = zeros(Nn) + Cn[1:Ne] += Cn_circular_wall_inner.(dn[1:Ne], dn[2:Nn], cp, ρ) .* dx + Cn[2:Nn] += Cn_circular_wall_outer.(dn[1:Ne], dn[2:Nn], cp, ρ) .* dx + p = @parameters C_shift = 0.0 + Ke = Ke_circular_wall.(dn[1:Ne], dn[2:Nn], λ) .* dx + @variables begin + (Tn(t))[1:N, 1:Nn] = fill(T0, N, Nn) + (Qe(t))[1:N, 1:Ne] = fill(T0, N, Ne) + end + sts = [vec(Tn); vec(Qe)] + e0 = Equation[inner_heatport.T[i] ~ Tn[i, 1] for i in 1:N] + e1 = Equation[outer_heatport.T[i] ~ Tn[i, Nn] for i in 1:N] + e2 = Equation[Qe[i, j] ~ Ke[j] * (-Tn[i, j+1] + Tn[i, j]) for i in 1:N for j in 1:Ne] + e3 = Equation[D(Tn[i, 1]) * (Cn[1] + C_shift) ~ inner_heatport.Q[i] - Qe[i, 1] for i in 1:N] + e4 = Equation[D(Tn[i, j]) * Cn[j] ~ Qe[i, j-1] - Qe[i, j] for i in 1:N for j in 2:Nn-1] + e5 = Equation[D(Tn[i, Nn]) * Cn[Nn] ~ Qe[i, Ne] + outer_heatport.Q[i] for i in 1:N] + eqs = vcat(e0, e1, e2, e3, e4, e5) + ODESystem(eqs, t, sts, p; systems=[inner_heatport, outer_heatport], name=name) +end + +function CylindricalSurfaceConvection(; name, L=100, N=100, d=1.0, α=5.0) + dx = L / N + S = pi * d * dx + @named heatport = VectorHeatPort(N=N) + sts = @variables Tenv(t) = 0.0 + eqs = [ + Tenv ~ 18.0 + [heatport.Q[i] ~ α * S * (heatport.T[i] - Tenv) for i in 1:N] + ] + + ODESystem(eqs, t, sts, []; systems=[heatport], name=name) +end + +function PreinsulatedPipe(; name, L=100.0, N=100.0, dn=0.05, T0=0.0, t_layer=[0.004, 0.013], + λ=[50, 0.04], cp=[500, 1200], ρ=[7800, 40], α=5.0, + e=1e-4, lumped_T=50, diffusion=true) + @named inlet = FluidPort() + @named outlet = FluidPort() + @named fluid_region = FluidRegion(L=L, N=N, dn=dn, e=e, lumped_T=lumped_T, diffusion=diffusion) + @named shell = CircularWallFEM(L=L, N=N, d=dn, t_layer=t_layer, λ=λ, cp=cp, ρ=ρ) + @named surfconv = CylindricalSurfaceConvection(L=L, N=N, d=dn + 2.0 * sum(t_layer), α=α) + systems = [fluid_region, shell, inlet, outlet, surfconv] + eqs = [ + connect(fluid_region.inlet, inlet) + connect(fluid_region.outlet, outlet) + connect(fluid_region.heatport, shell.inner_heatport) + connect(shell.outer_heatport, surfconv.heatport) + ] + ODESystem(eqs, t, [], []; systems=systems, name=name) +end + +function Source(; name, p_feed=100000) + @named outlet = FluidPort() + sts = @variables m_flow(t) = 1e-6 + eqs = [ + m_flow ~ m_flow_source(t) + outlet.m ~ -m_flow + outlet.p ~ p_feed + outlet.T ~ T_source(t) + ] + compose(ODESystem(eqs, t, sts, []; name=name), [outlet]) +end + +function Sink(; name) + @named inlet = FluidPort() + eqs = [ + inlet.T ~ instream(inlet.T) + ] + compose(ODESystem(eqs, t, [], []; name=name), [inlet]) +end + +function TestBenchPreinsulated(; name, L=1.0, dn=0.05, t_layer=[0.0056, 0.013], N=100, diffusion=true, lumped_T=20) + @named pipe = PreinsulatedPipe(L=L, dn=dn, N=N, diffusion=diffusion, t_layer=t_layer, lumped_T=lumped_T) + @named source = Source() + @named sink = Sink() + subs = [source, pipe, sink] + eqs = [ + connect(source.outlet, pipe.inlet) + connect(pipe.outlet, sink.inlet) + ] + compose(ODESystem(eqs, t, [], []; name=name), subs) +end + +function run_and_time!(jacobian_times, build_times, first_call_times, second_call_times, i, N) + @mtkbuild sys = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058]) + jac_result = @timed calculate_jacobian(sys; sparse = true) + jac = jac_result.value + jacobian_times[i] = jac_result.time + + dvs = unknowns(sys) + ps = parameters(sys) + defs = defaults(sys) + u0 = Float64[Symbolics.fixpoint_sub(v, defs) for v in dvs] + p = Float64[Symbolics.fixpoint_sub(v, defs) for v in ps] + t0 = 0.0 + buffer_nocse = similar(jac, Float64) + buffer_nocse.nzval .= 0.0 + buffer_cse = similar(jac, Float64) + buffer_cse.nzval .= 0.0 + + build_result_nocse = @timed build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false) + f_jac_nocse = build_result_nocse.value[2] + build_times[1][i] = build_result_nocse.time + first_call_result_nocse = @timed f_jac_nocse(buffer_nocse, u0, p, t0) + first_call_times[1][i] = first_call_result_nocse.time + second_call_result_nocse = @timed f_jac_nocse(buffer_nocse, u0, p, t0) + second_call_times[1][i] = second_call_result_nocse.time + + build_result_cse = @timed build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true) + f_jac_cse = build_result_cse.value[2] + build_times[2][i] = build_result_cse.time + first_call_result_cse = @timed f_jac_cse(buffer_cse, u0, p, t0) + first_call_times[2][i] = first_call_result_cse.time + second_call_result_cse = @timed f_jac_cse(buffer_cse, u0, p, t0) + second_call_times[2][i] = second_call_result_cse.time + + return nothing +end +``` + +```julia +N = [5, 10, 20, 40, 60, 80, 160, 320, 480]; +jacobian_times = zeros(Float64, length(N)) +# [without_cse_times, with_cse_times] +build_times = [copy(jacobian_times), copy(jacobian_times)] +first_call_times = copy.(build_times) +second_call_times = copy.(build_times) +``` + +## Julia Timings + +```julia +for (i, n) in enumerate(N) + run_and_time!(jacobian_times, build_times, first_call_times, second_call_times, i, n) +end +``` + +## Appendix + +```julia, echo = false +using SciMLBenchmarks +SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file]) +``` From 5181954530b70efa33948b1c381de68f1e1323f1 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 19 Mar 2025 15:39:45 +0530 Subject: [PATCH 04/31] feat: add plots for Symbolics/ThermalFluid benchmark --- benchmarks/Symbolics/Project.toml | 1 + benchmarks/Symbolics/ThermalFluid.jmd | 23 +++++++++++++++++++++-- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml index 1ccb42a57..cf3ca56b5 100644 --- a/benchmarks/Symbolics/Project.toml +++ b/benchmarks/Symbolics/Project.toml @@ -8,6 +8,7 @@ OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" +PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881" SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 309dd93b4..009febd1a 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -10,7 +10,7 @@ using Pkg # Rev fixes precompilation https://github.com/hzgzh/XSteam.jl/pull/2 Pkg.add(Pkg.PackageSpec(;name="XSteam", rev="f2a1c589054cfd6bba307985a3a534b6f5a1863b")) -using ModelingToolkit, Symbolics, XSteam, Polynomials, CairoMakie +using ModelingToolkit, Symbolics, XSteam, Polynomials, CairoMakie, PrettyTables using SparseArrays using ModelingToolkit: t_nounits as t, D_nounits as D ``` @@ -293,7 +293,7 @@ first_call_times = copy.(build_times) second_call_times = copy.(build_times) ``` -## Julia Timings +## Timings ```julia for (i, n) in enumerate(N) @@ -301,6 +301,25 @@ for (i, n) in enumerate(N) end ``` +## Results + +```julia +data = hcat(N, jacobian_times, build_times, first_call_times, second_call_times) +header = ["N", "Jacobian calculation time", "`build_function` time", "First call time", "Second call time"] +pretty_table(data; header) +``` + +```julia +f = Figure(size = (800, 1200)) +labels = header[2:end] +times = [jacobian_times, build_times, first_call_times, second_call_times] +for (i, (label, data)) in enumerate(zip(labels, times)) + ax = Axis(f[i, 1], xscale = log10, yscale = log10, ylabel = label) + lines!(ax, N, data) +end +f +``` + ## Appendix ```julia, echo = false From cdd4b281e4d9b1e8be39d3e5df20a7d423f7e1f8 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 19 Mar 2025 20:11:06 +0530 Subject: [PATCH 05/31] feat: track allocation statistics, measure with and without hashconsing --- benchmarks/Symbolics/ThermalFluid.jmd | 52 ++++++++++++++++++--------- 1 file changed, 36 insertions(+), 16 deletions(-) diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 009febd1a..607dd7ce2 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -247,11 +247,24 @@ function TestBenchPreinsulated(; name, L=1.0, dn=0.05, t_layer=[0.0056, 0.013], compose(ODESystem(eqs, t, [], []; name=name), subs) end -function run_and_time!(jacobian_times, build_times, first_call_times, second_call_times, i, N) +function run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, i, N) @mtkbuild sys = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058]) + SymbolicUtils.ENABLE_HASHCONSING[] = false jac_result = @timed calculate_jacobian(sys; sparse = true) - jac = jac_result.value - jacobian_times[i] = jac_result.time + jac_nocse = jac_result.value + jacobian_times[1][i] = jac_result.time + jacobian_gctimes[1][i] = jac_result.gctime + jacobian_allocs[1][i] = jac_result.bytes + SymbolicUtils.ENABLE_HASHCONSING[] = true + SymbolicUtils.clear_cache!(Symbolics.occursin_info) + jac_result = @timed calculate_jacobian(sys; sparse = true) + jac_cse = jac_result.value + jacobian_times[2][i] = jac_result.time + jacobian_gctimes[2][i] = jac_result.gctime + jacobian_allocs[2][i] = jac_result.bytes + + @assert isequal(jac_nocse, jac_cse) + jac = jac_cse dvs = unknowns(sys) ps = parameters(sys) @@ -286,36 +299,43 @@ end ```julia N = [5, 10, 20, 40, 60, 80, 160, 320, 480]; -jacobian_times = zeros(Float64, length(N)) +jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] +jacobian_gctimes = copy.(jacobian_times) +jacobian_allocs = [zeros(Int64, length(N)), zeros(Int64, length(N))] # [without_cse_times, with_cse_times] -build_times = [copy(jacobian_times), copy(jacobian_times)] -first_call_times = copy.(build_times) -second_call_times = copy.(build_times) +build_times = copy.(jacobian_times) +first_call_times = copy.(jacobian_times) +second_call_times = copy.(jacobian_times) ``` ## Timings ```julia +# compile +run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, 1, 5) for (i, n) in enumerate(N) - run_and_time!(jacobian_times, build_times, first_call_times, second_call_times, i, n) + run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, i, n) end ``` ## Results ```julia -data = hcat(N, jacobian_times, build_times, first_call_times, second_call_times) -header = ["N", "Jacobian calculation time", "`build_function` time", "First call time", "Second call time"] -pretty_table(data; header) +tabledata = hcat(N, jacobian_times..., jacobian_gctimes..., jacobian_allocs..., build_times..., first_call_times..., second_call_times...) +header = ["N", "Jacobian time (no hashconsing)", "Jacobian time (hashconsing)", "Jacobian GC time (no hashconsing)", "Jacobian GC time (hashconsing)", "Jacobian allocated memory (no hashconsing) (B)", "Jacobian allocated memory (hashconsing) (B)", "`build_function` time (no CSE)", "`build_function` time (CSE)", "First call time (no CSE)", "First call time (CSE)", "Second call time (no CSE)", "Second call time (CSE)"] +pretty_table(tabledata; header, backend = Val(:html)) ``` ```julia -f = Figure(size = (800, 1200)) -labels = header[2:end] -times = [jacobian_times, build_times, first_call_times, second_call_times] +f = Figure(size = (800, 1200)); +labels = ["Symbolic jacobian time", "Symbolic jacobian\nGC time", "Symbolic jacobian\nallocated memory (B)", "`build_function` time", "First call time", "Second call time"] +times = [jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times] for (i, (label, data)) in enumerate(zip(labels, times)) - ax = Axis(f[i, 1], xscale = log10, yscale = log10, ylabel = label) - lines!(ax, N, data) + ax = Axis(f[i, 1], xscale = log10, yscale = log10, ylabel = label, xticks = N) + l1 = lines!(ax, N, data[1]) + l2 = lines!(ax, N, data[2]) + legend_entries = startswith(label, "Symbolic") ? ["without hashconsing", "with hashconsing"] : ["without CSE", "with CSE"] + Legend(f[i, 2], [l1, l2], legend_entries) end f ``` From a15f26562f1654d6c64c2abac4575118cc92728d Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 19 Mar 2025 20:12:01 +0530 Subject: [PATCH 06/31] build: bump Symbolics, SymbolicUtils compats --- benchmarks/Symbolics/Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml index cf3ca56b5..e4cbcc1df 100644 --- a/benchmarks/Symbolics/Project.toml +++ b/benchmarks/Symbolics/Project.toml @@ -19,5 +19,5 @@ XSteam = "95ff35a0-be81-11e9-2ca3-5b4e338e8476" [compat] ModelingToolkit = "9" -SymbolicUtils = "3.17" -Symbolics = "6.29.1" +SymbolicUtils = "3.20" +Symbolics = "6.31" From b2fd9c544683cb2f4aa90ecaee876d8107293a52 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 20 Mar 2025 01:04:24 +0530 Subject: [PATCH 07/31] refactor: improve ThermalFluid benchmark --- benchmarks/Symbolics/Project.toml | 2 + benchmarks/Symbolics/ThermalFluid.jmd | 60 ++++++++++++++------------- 2 files changed, 34 insertions(+), 28 deletions(-) diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml index e4cbcc1df..58cbfb5d0 100644 --- a/benchmarks/Symbolics/Project.toml +++ b/benchmarks/Symbolics/Project.toml @@ -1,6 +1,8 @@ [deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" +Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 607dd7ce2..1b87a3d2f 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -10,8 +10,8 @@ using Pkg # Rev fixes precompilation https://github.com/hzgzh/XSteam.jl/pull/2 Pkg.add(Pkg.PackageSpec(;name="XSteam", rev="f2a1c589054cfd6bba307985a3a534b6f5a1863b")) -using ModelingToolkit, Symbolics, XSteam, Polynomials, CairoMakie, PrettyTables -using SparseArrays +using ModelingToolkit, Symbolics, SymbolicUtils, XSteam, Polynomials, CairoMakie, PrettyTables +using SparseArrays, Chairmarks using ModelingToolkit: t_nounits as t, D_nounits as D ``` @@ -247,21 +247,25 @@ function TestBenchPreinsulated(; name, L=1.0, dn=0.05, t_layer=[0.0056, 0.013], compose(ODESystem(eqs, t, [], []; name=name), subs) end +function call(fn, args...) + fn(args...) +end + function run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, i, N) @mtkbuild sys = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058]) SymbolicUtils.ENABLE_HASHCONSING[] = false - jac_result = @timed calculate_jacobian(sys; sparse = true) - jac_nocse = jac_result.value - jacobian_times[1][i] = jac_result.time - jacobian_gctimes[1][i] = jac_result.gctime - jacobian_allocs[1][i] = jac_result.bytes + jac_result = @be calculate_jacobian(sys; sparse = true) + jac_nocse = calculate_jacobian(sys; sparse = true) + jacobian_times[1][i] = minimum(x -> x.time, jac_result.samples) + jacobian_gctimes[1][i] = minimum(x -> x.time * x.gc_fraction, jac_result.samples) + jacobian_allocs[1][i] = minimum(x -> x.bytes, jac_result.samples) SymbolicUtils.ENABLE_HASHCONSING[] = true SymbolicUtils.clear_cache!(Symbolics.occursin_info) - jac_result = @timed calculate_jacobian(sys; sparse = true) - jac_cse = jac_result.value - jacobian_times[2][i] = jac_result.time - jacobian_gctimes[2][i] = jac_result.gctime - jacobian_allocs[2][i] = jac_result.bytes + jac_result = @be calculate_jacobian(sys; sparse = true) + jac_cse = calculate_jacobian(sys; sparse = true) + jacobian_times[2][i] = minimum(x -> x.time, jac_result.samples) + jacobian_gctimes[2][i] = minimum(x -> x.time * x.gc_fraction, jac_result.samples) + jacobian_allocs[2][i] = minimum(x -> x.bytes, jac_result.samples) @assert isequal(jac_nocse, jac_cse) jac = jac_cse @@ -277,21 +281,21 @@ function run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_ buffer_cse = similar(jac, Float64) buffer_cse.nzval .= 0.0 - build_result_nocse = @timed build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false) - f_jac_nocse = build_result_nocse.value[2] - build_times[1][i] = build_result_nocse.time - first_call_result_nocse = @timed f_jac_nocse(buffer_nocse, u0, p, t0) - first_call_times[1][i] = first_call_result_nocse.time - second_call_result_nocse = @timed f_jac_nocse(buffer_nocse, u0, p, t0) - second_call_times[1][i] = second_call_result_nocse.time - - build_result_cse = @timed build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true) - f_jac_cse = build_result_cse.value[2] - build_times[2][i] = build_result_cse.time - first_call_result_cse = @timed f_jac_cse(buffer_cse, u0, p, t0) - first_call_times[2][i] = first_call_result_cse.time - second_call_result_cse = @timed f_jac_cse(buffer_cse, u0, p, t0) - second_call_times[2][i] = second_call_result_cse.time + build_result_nocse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false) + f_jac_nocse = build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false)[2] + build_times[1][i] = minimum(x -> x.time, build_result_nocse.samples) + first_call_result_nocse = @be call(f_jac_nocse, buffer_nocse, u0, p, t0) + first_call_times[1][i] = minimum(x -> x.time, first_call_result_nocse.samples) + second_call_result_nocse = @be call(f_jac_nocse, buffer_nocse, u0, p, t0) + second_call_times[1][i] = minimum(x -> x.time, second_call_result_nocse.samples) + + build_result_cse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true) + f_jac_cse = build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true)[2] + build_times[2][i] = minimum(x -> x.time, build_result_cse.samples) + first_call_result_cse = @be call(f_jac_cse, buffer_cse, u0, p, t0) + first_call_times[2][i] = minimum(x -> x.time, first_call_result_cse.samples) + second_call_result_cse = @be call(f_jac_cse, buffer_cse, u0, p, t0) + second_call_times[2][i] = minimum(x -> x.time, second_call_result_cse.samples) return nothing end @@ -301,7 +305,7 @@ end N = [5, 10, 20, 40, 60, 80, 160, 320, 480]; jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] jacobian_gctimes = copy.(jacobian_times) -jacobian_allocs = [zeros(Int64, length(N)), zeros(Int64, length(N))] +jacobian_allocs = copy.(jacobian_times) # [without_cse_times, with_cse_times] build_times = copy.(jacobian_times) first_call_times = copy.(jacobian_times) From d5cd9d1a4e994467987f0a8a9dd0762d64a86080 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 20 Mar 2025 01:05:08 +0530 Subject: [PATCH 08/31] fix: fix BCR benchmark --- benchmarks/Symbolics/BCR.jmd | 96 ++++++++++++++++++++++++++++++++++-- 1 file changed, 92 insertions(+), 4 deletions(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 9bea261a5..f7de8b6a9 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -18,8 +18,9 @@ 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 + TimerOutputs, LinearAlgebra, ModelingToolkit, Chairmarks, + LinearSolve, Symbolics, SymbolicUtils, SymbolicUtils.Code, SparseArrays, CairoMakie, + PrettyTables datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() @@ -60,8 +61,8 @@ kwargs = (; iip_config = (false, true), expression = Val{false}) defs = defaults(osys) u = Float64[Symbolics.fixpoint_sub(var, defs) for var in vars] -buffer_cse = similar(Float64, jac) -buffer_nocse = similar(Float64, jac) +buffer_cse = similar(jac, Float64) +buffer_nocse = similar(jac, Float64) p = Float64[Symbolics.fixpoint_sub(par, defs) for par in pars] tt = 0.0 @@ -75,3 +76,90 @@ tt = 0.0 show(to) ``` + +We'll also measure scaling. + + +```julia +function run_and_time!(rhs, vars, pars, iv, u0, p, t0, N, i, jac_times, jac_allocs, build_times, first_call_times, second_call_times) + outputs = rhs[1:N] + SymbolicUtils.ENABLE_HASHCONSING[] = false + jac_result = @be Symbolics.sparsejacobian(outputs, vars) + jac_times[1][i] = minimum(x -> x.time, jac_result.samples) + jac_allocs[1][i] = minimum(x -> x.bytes, jac_result.samples) + jac_nohc = Symbolics.sparsejacobian(outputs, vars) + + SymbolicUtils.ENABLE_HASHCONSING[] = true + SymbolicUtils.clear_cache!(Symbolics.occursin_info) + jac_result = @be Symbolics.sparsejacobian(outputs, vars) + jac_times[2][i] = minimum(x -> x.time, jac_result.samples) + jac_allocs[2][i] = minimum(x -> x.bytes, jac_result.samples) + jac_hc = Symbolics.sparsejacobian(outputs, vars) + + @assert isequal(jac_nohc, jac_hc) + jac = jac_hc + args = (vars, pars, iv) + kwargs = (; iip_config = (false, true), expression = Val{false}) + + build_result = @be build_function(jac, args...; cse = false, kwargs...); + build_times[1][i] = minimum(x -> x.time, build_result.samples) + jacfn_nocse = build_function(jac, args...; cse = false, kwargs...)[2] + + build_result = @be build_function(jac, args...; cse = true, kwargs...); + build_times[2][i] = minimum(x -> x.time, build_result.samples) + jacfn_cse = build_function(jac, args...; cse = false, kwargs...)[2] + + buffer = similar(jac, Float64) + call_result = @be jacfn_nocse(buffer, u0, p, t0) + first_call_times[1][i] = minimum(x -> x.time, call_result.samples) + call_result = @be jacfn_cse(buffer, u0, p, t0) + first_call_times[2][i] = minimum(x -> x.time, call_result.samples) + + call_result = @be jacfn_nocse(buffer, u0, p, t0) + second_call_times[1][i] = minimum(x -> x.time, call_result.samples) + call_result = @be jacfn_cse(buffer, u0, p, t0) + second_call_times[2][i] = minimum(x -> x.time, call_result.samples) + + return nothing +end +``` + +# Run benchmark + +```julia +N = [10, 20, 40, 60, 100, 200, 300, 400] +jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] +jacobian_allocs = copy.(jacobian_times) +# [without_cse_times, with_cse_times] +build_times = copy.(jacobian_times) +first_call_times = copy.(jacobian_times) +second_call_times = copy.(jacobian_times) + +iv = ModelingToolkit.get_iv(osys) +run_and_time!(rhs, vars, pars, iv, u, p, tt, 10, 1, jacobian_times, jacobian_allocs, build_times, first_call_times, second_call_times) +for (i, n) in enumerate(N) + run_and_time!(rhs, vars, pars, iv, u, p, tt, n, i, jacobian_times, jacobian_allocs, build_times, first_call_times, second_call_times) +end +``` + +# Plot figures + +```julia +tabledata = hcat(N, jacobian_times..., jacobian_allocs..., build_times..., first_call_times..., second_call_times...) +header = ["N", "Jacobian time (no hashconsing)", "Jacobian time (hashconsing)", "Jacobian allocated memory (no hashconsing) (B)", "Jacobian allocated memory (hashconsing) (B)", "`build_function` time (no CSE)", "`build_function` time (CSE)", "First call time (no CSE)", "First call time (CSE)", "Second call time (no CSE)", "Second call time (CSE)"] +pretty_table(tabledata; header, backend = Val(:html)) +``` + +```julia +f = Figure(size = (800, 1200)); +labels = ["Symbolic jacobian time", "Symbolic jacobian\nallocated memory (B)", "`build_function` time", "First call time", "Second call time"] +times = [jacobian_times, jacobian_allocs, build_times, first_call_times, second_call_times] +for (i, (label, data)) in enumerate(zip(labels, times)) + ax = Axis(f[i, 1], xscale = log10, yscale = log10, ylabel = label, xticks = N) + l1 = lines!(ax, N, data[1]) + l2 = lines!(ax, N, data[2]) + legend_entries = startswith(label, "Symbolic") ? ["without hashconsing", "with hashconsing"] : ["without CSE", "with CSE"] + Legend(f[i, 2], [l1, l2], legend_entries) +end +f +``` From f34959356f1492914d2afe2b6a4b7cf0c358d33c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 20 Mar 2025 11:42:07 +0530 Subject: [PATCH 09/31] fix: fix BCR benchmark --- benchmarks/Symbolics/BCR.jmd | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index f7de8b6a9..06dbe4e0b 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -67,10 +67,10 @@ p = Float64[Symbolics.fixpoint_sub(par, defs) for par in pars] tt = 0.0 @timeit to "Compile jacobian - CSE" jac_cse_iip(buffer_cse, u, p, tt) -@timeit to "Compute jacobian - CSE" jac_cse_iip(buffer_cse, u, p, t) +@timeit to "Compute jacobian - CSE" jac_cse_iip(buffer_cse, u, p, tt) @timeit to "Compile jacobian - no CSE" jac_nocse_iip(buffer_nocse, u, p, tt) -@timeit to "Compute jacobian - no CSE" jac_nocse_iip(buffer_nocse, u, p, t) +@timeit to "Compute jacobian - no CSE" jac_nocse_iip(buffer_nocse, u, p, tt) @assert isapprox(buffer_cse, buffer_nocse, rtol = 1e-10) @@ -107,7 +107,7 @@ function run_and_time!(rhs, vars, pars, iv, u0, p, t0, N, i, jac_times, jac_allo build_result = @be build_function(jac, args...; cse = true, kwargs...); build_times[2][i] = minimum(x -> x.time, build_result.samples) - jacfn_cse = build_function(jac, args...; cse = false, kwargs...)[2] + jacfn_cse = build_function(jac, args...; cse = true, kwargs...)[2] buffer = similar(jac, Float64) call_result = @be jacfn_nocse(buffer, u0, p, t0) From 274b3faaffc5a6e586465dfa241bddeef6c7374d Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sun, 23 Mar 2025 08:55:38 +0530 Subject: [PATCH 10/31] build: bump Symbolics, SymbolicUtils compats --- benchmarks/Symbolics/Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml index 58cbfb5d0..ecbd75ba0 100644 --- a/benchmarks/Symbolics/Project.toml +++ b/benchmarks/Symbolics/Project.toml @@ -21,5 +21,5 @@ XSteam = "95ff35a0-be81-11e9-2ca3-5b4e338e8476" [compat] ModelingToolkit = "9" -SymbolicUtils = "3.20" -Symbolics = "6.31" +SymbolicUtils = "3.21" +Symbolics = "6.33" From ccd7a105b9360a92a38ecfb46c007d42077fea6c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sun, 23 Mar 2025 13:08:36 +0530 Subject: [PATCH 11/31] refactor: use `mean` for ThermalFluid timings --- benchmarks/Symbolics/ThermalFluid.jmd | 30 +++++++++++++-------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 1b87a3d2f..5066d3c40 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -11,7 +11,7 @@ using Pkg Pkg.add(Pkg.PackageSpec(;name="XSteam", rev="f2a1c589054cfd6bba307985a3a534b6f5a1863b")) using ModelingToolkit, Symbolics, SymbolicUtils, XSteam, Polynomials, CairoMakie, PrettyTables -using SparseArrays, Chairmarks +using SparseArrays, Chairmarks, Statistics using ModelingToolkit: t_nounits as t, D_nounits as D ``` @@ -256,16 +256,16 @@ function run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_ SymbolicUtils.ENABLE_HASHCONSING[] = false jac_result = @be calculate_jacobian(sys; sparse = true) jac_nocse = calculate_jacobian(sys; sparse = true) - jacobian_times[1][i] = minimum(x -> x.time, jac_result.samples) - jacobian_gctimes[1][i] = minimum(x -> x.time * x.gc_fraction, jac_result.samples) - jacobian_allocs[1][i] = minimum(x -> x.bytes, jac_result.samples) + jacobian_times[1][i] = mean(x -> x.time, jac_result.samples) + jacobian_gctimes[1][i] = mean(x -> x.time * x.gc_fraction, jac_result.samples) + jacobian_allocs[1][i] = mean(x -> x.bytes, jac_result.samples) SymbolicUtils.ENABLE_HASHCONSING[] = true SymbolicUtils.clear_cache!(Symbolics.occursin_info) jac_result = @be calculate_jacobian(sys; sparse = true) jac_cse = calculate_jacobian(sys; sparse = true) - jacobian_times[2][i] = minimum(x -> x.time, jac_result.samples) - jacobian_gctimes[2][i] = minimum(x -> x.time * x.gc_fraction, jac_result.samples) - jacobian_allocs[2][i] = minimum(x -> x.bytes, jac_result.samples) + jacobian_times[2][i] = mean(x -> x.time, jac_result.samples) + jacobian_gctimes[2][i] = mean(x -> x.time * x.gc_fraction, jac_result.samples) + jacobian_allocs[2][i] = mean(x -> x.bytes, jac_result.samples) @assert isequal(jac_nocse, jac_cse) jac = jac_cse @@ -281,21 +281,21 @@ function run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_ buffer_cse = similar(jac, Float64) buffer_cse.nzval .= 0.0 - build_result_nocse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false) f_jac_nocse = build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false)[2] - build_times[1][i] = minimum(x -> x.time, build_result_nocse.samples) + build_result_nocse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false) + build_times[1][i] = mean(x -> x.time, build_result_nocse.samples) first_call_result_nocse = @be call(f_jac_nocse, buffer_nocse, u0, p, t0) - first_call_times[1][i] = minimum(x -> x.time, first_call_result_nocse.samples) + first_call_times[1][i] = mean(x -> x.time, first_call_result_nocse.samples) second_call_result_nocse = @be call(f_jac_nocse, buffer_nocse, u0, p, t0) - second_call_times[1][i] = minimum(x -> x.time, second_call_result_nocse.samples) + second_call_times[1][i] = mean(x -> x.time, second_call_result_nocse.samples) - build_result_cse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true) f_jac_cse = build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true)[2] - build_times[2][i] = minimum(x -> x.time, build_result_cse.samples) + build_result_cse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true) + build_times[2][i] = mean(x -> x.time, build_result_cse.samples) first_call_result_cse = @be call(f_jac_cse, buffer_cse, u0, p, t0) - first_call_times[2][i] = minimum(x -> x.time, first_call_result_cse.samples) + first_call_times[2][i] = mean(x -> x.time, first_call_result_cse.samples) second_call_result_cse = @be call(f_jac_cse, buffer_cse, u0, p, t0) - second_call_times[2][i] = minimum(x -> x.time, second_call_result_cse.samples) + second_call_times[2][i] = mean(x -> x.time, second_call_result_cse.samples) return nothing end From d161cab4f96ed6ad7d8862128535f44e6f09fac9 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sun, 23 Mar 2025 13:08:50 +0530 Subject: [PATCH 12/31] refactor: avoid registering trivial functions in ThermalFluid --- benchmarks/Symbolics/ThermalFluid.jmd | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 5066d3c40..94af82490 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -27,8 +27,8 @@ using ModelingToolkit: t_nounits as t, D_nounits as D m_flow_source(t) = 2.75 T_source(t) = (t > 12 * 3600) * 56.0 + 12.0 -@register_symbolic m_flow_source(t) -@register_symbolic T_source(t) +# @register_symbolic m_flow_source(t) +# @register_symbolic T_source(t) #build polynomial liquid-water property only dependent on Temperature p_l = 5 #bar @@ -38,11 +38,11 @@ T_vec = collect(1:1:150); @generated Pr_T(t) = :(Base.evalpoly(t, $(fit(T_vec, 1e3 * Cp_pT.(p_l, T_vec) .* my_pT.(p_l, T_vec) ./ tc_pT.(p_l, T_vec), 5).coeffs...,))) @generated rho_T(t) = :(Base.evalpoly(t, $(fit(T_vec, rho_pT.(p_l, T_vec), 4).coeffs...,))) @generated rhocp_T(t) = :(Base.evalpoly(t, $(fit(T_vec, 1000 * rho_pT.(p_l, T_vec) .* Cp_pT.(p_l, T_vec), 5).coeffs...,))) -@register_symbolic kin_visc_T(t) -@register_symbolic lambda_T(t) -@register_symbolic Pr_T(t) -@register_symbolic rho_T(t) -@register_symbolic rhocp_T(t) +# @register_symbolic kin_visc_T(t) +# @register_symbolic lambda_T(t) +# @register_symbolic Pr_T(t) +# @register_symbolic rho_T(t) +# @register_symbolic rhocp_T(t) @connector function FluidPort(; name, p=101325.0, m=0.0, T=0.0) sts = @variables p(t) = p m(t) = m [connect = Flow] T(t) = T [connect = Stream] @@ -77,7 +77,7 @@ function Nusselt(Re, Pr, f) end end -@register_symbolic Churchill_f(Re, epsilon, d) +# @register_symbolic Churchill_f(Re, epsilon, d) #Darcy weisbach friction factor function Churchill_f(Re, epsilon, d) theta_1 = (-2.457 * log(((7 / Re)^0.9) + (0.27 * (epsilon / d))))^16 @@ -139,19 +139,19 @@ function FluidRegion(; name, L=1.0, dn=0.05, N=100, T0=0.0, ODESystem(eqs, t, sts, p; systems=[inlet, outlet, heatport], name=name) end -@register_symbolic Cn_circular_wall_inner(d, D, cp, ρ) +# @register_symbolic Cn_circular_wall_inner(d, D, cp, ρ) function Cn_circular_wall_inner(d, D, cp, ρ) C = pi / 4 * (D^2 - d^2) * cp * ρ return C / 2 end -@register_symbolic Cn_circular_wall_outer(d, D, cp, ρ) +# @register_symbolic Cn_circular_wall_outer(d, D, cp, ρ) function Cn_circular_wall_outer(d, D, cp, ρ) C = pi / 4 * (D^2 - d^2) * cp * ρ return C / 2 end -@register_symbolic Ke_circular_wall(d, D, λ) +# @register_symbolic Ke_circular_wall(d, D, λ) function Ke_circular_wall(d, D, λ) 2 * pi * λ / log(D / d) end From 00cb3084c6bec5de522991111f0128940732a384 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sun, 23 Mar 2025 16:05:13 +0530 Subject: [PATCH 13/31] build: bump Symbolics compat --- benchmarks/Symbolics/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml index ecbd75ba0..d882ab18f 100644 --- a/benchmarks/Symbolics/Project.toml +++ b/benchmarks/Symbolics/Project.toml @@ -22,4 +22,4 @@ XSteam = "95ff35a0-be81-11e9-2ca3-5b4e338e8476" [compat] ModelingToolkit = "9" SymbolicUtils = "3.21" -Symbolics = "6.33" +Symbolics = "6.33.1" From 39889ca0301ee594d8a925e29b5e0ef912f1c9cf Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 10:02:40 +0530 Subject: [PATCH 14/31] refactor: run for smaller `N`, add logging --- benchmarks/Symbolics/ThermalFluid.jmd | 31 +++++++++++++++++++++------ 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 94af82490..820156126 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -253,20 +253,26 @@ end function run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, i, N) @mtkbuild sys = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058]) + @info "Built system" N SymbolicUtils.ENABLE_HASHCONSING[] = false jac_result = @be calculate_jacobian(sys; sparse = true) + @info "No hashconsing benchmark" jac_nocse = calculate_jacobian(sys; sparse = true) + @info "No hashconsing result" jacobian_times[1][i] = mean(x -> x.time, jac_result.samples) jacobian_gctimes[1][i] = mean(x -> x.time * x.gc_fraction, jac_result.samples) jacobian_allocs[1][i] = mean(x -> x.bytes, jac_result.samples) + @info "times" jacobian_times[1][i] jacobian_gctimes[1][i] jacobian_allocs[1][i] SymbolicUtils.ENABLE_HASHCONSING[] = true SymbolicUtils.clear_cache!(Symbolics.occursin_info) jac_result = @be calculate_jacobian(sys; sparse = true) + @info "Hashconsing benchmark" jac_cse = calculate_jacobian(sys; sparse = true) + @info "Hashconsing result" jacobian_times[2][i] = mean(x -> x.time, jac_result.samples) jacobian_gctimes[2][i] = mean(x -> x.time * x.gc_fraction, jac_result.samples) jacobian_allocs[2][i] = mean(x -> x.bytes, jac_result.samples) - + @info "times" jacobian_times[2][i] jacobian_gctimes[2][i] jacobian_allocs[2][i] @assert isequal(jac_nocse, jac_cse) jac = jac_cse @@ -282,27 +288,37 @@ function run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_ buffer_cse.nzval .= 0.0 f_jac_nocse = build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false)[2] + @info "No CSE build_function result" build_result_nocse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false) + @info "No CSE build_function benchmark" build_times[1][i] = mean(x -> x.time, build_result_nocse.samples) - first_call_result_nocse = @be call(f_jac_nocse, buffer_nocse, u0, p, t0) - first_call_times[1][i] = mean(x -> x.time, first_call_result_nocse.samples) + @info "time" build_times[1][i] + first_call_result_nocse = @timed call(f_jac_nocse, buffer_nocse, u0, p, t0) + first_call_times[1][i] = first_call_result_nocse.time + @info "time" first_call_times[1][i] second_call_result_nocse = @be call(f_jac_nocse, buffer_nocse, u0, p, t0) second_call_times[1][i] = mean(x -> x.time, second_call_result_nocse.samples) + @info "time" second_call_times[1][i] f_jac_cse = build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true)[2] + @info "CSE build_function result" build_result_cse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true) + @info "CSE build_function benchmark" build_times[2][i] = mean(x -> x.time, build_result_cse.samples) - first_call_result_cse = @be call(f_jac_cse, buffer_cse, u0, p, t0) - first_call_times[2][i] = mean(x -> x.time, first_call_result_cse.samples) + @info "time" build_times[2][i] + first_call_result_cse = @timed call(f_jac_cse, buffer_cse, u0, p, t0) + first_call_times[2][i] = first_call_result_cse.time + @info "time" first_call_times[2][i] second_call_result_cse = @be call(f_jac_cse, buffer_cse, u0, p, t0) second_call_times[2][i] = mean(x -> x.time, second_call_result_cse.samples) + @info "time" second_call_times[2][i] return nothing end ``` ```julia -N = [5, 10, 20, 40, 60, 80, 160, 320, 480]; +N = [5, 10, 20, 40, 80, 160, 320]; jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] jacobian_gctimes = copy.(jacobian_times) jacobian_allocs = copy.(jacobian_times) @@ -318,7 +334,8 @@ second_call_times = copy.(jacobian_times) # compile run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, 1, 5) for (i, n) in enumerate(N) - run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, i, n) + @info i n + @time run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, i, n) end ``` From 47b41c9e6f579791c1c816ee96587062942db3f1 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 10:35:50 +0530 Subject: [PATCH 15/31] fix: fix first call time calculation --- benchmarks/Symbolics/BCR.jmd | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 06dbe4e0b..6bbb94c11 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -110,10 +110,10 @@ function run_and_time!(rhs, vars, pars, iv, u0, p, t0, N, i, jac_times, jac_allo jacfn_cse = build_function(jac, args...; cse = true, kwargs...)[2] buffer = similar(jac, Float64) - call_result = @be jacfn_nocse(buffer, u0, p, t0) - first_call_times[1][i] = minimum(x -> x.time, call_result.samples) - call_result = @be jacfn_cse(buffer, u0, p, t0) - first_call_times[2][i] = minimum(x -> x.time, call_result.samples) + call_result = @timed jacfn_nocse(buffer, u0, p, t0) + first_call_times[1][i] = call_result.time + call_result = @timed jacfn_cse(buffer, u0, p, t0) + first_call_times[2][i] = call_result.time call_result = @be jacfn_nocse(buffer, u0, p, t0) second_call_times[1][i] = minimum(x -> x.time, call_result.samples) From de7ff53f33c599d826533945d804ba3edf5c899f Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 16:35:48 +0530 Subject: [PATCH 16/31] fix: use `eval` instead of RGFs --- benchmarks/Symbolics/BCR.jmd | 59 ++++++++++++++++------ benchmarks/Symbolics/ThermalFluid.jmd | 73 ++++++++++++++++++--------- 2 files changed, 92 insertions(+), 40 deletions(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 6bbb94c11..b98ce3c87 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -33,7 +33,7 @@ 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)] +rhs = [eq.rhs for eq in full_equations(osys)][1:10] vars = unknowns(osys) pars = parameters(osys) @@ -55,10 +55,13 @@ 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. `iip_config` prevents it # from trying to build the function. -kwargs = (; iip_config = (false, true), expression = Val{false}) +kwargs = (; iip_config = (false, true), expression = Val{true}) @timeit to "Build jacobian - no CSE" _, jac_nocse_iip = build_function(jac, args...; cse = false, kwargs...); @timeit to "Build jacobian - CSE" _, jac_cse_iip = build_function(jac, args...; cse = true, kwargs...); +jac_nocse_iip = eval(jac_nocse_iip) +jac_cse_iip = eval(jac_cse_iip) + defs = defaults(osys) u = Float64[Symbolics.fixpoint_sub(var, defs) for var in vars] buffer_cse = similar(jac, Float64) @@ -81,7 +84,7 @@ We'll also measure scaling. ```julia -function run_and_time!(rhs, vars, pars, iv, u0, p, t0, N, i, jac_times, jac_allocs, build_times, first_call_times, second_call_times) +function run_and_time_construct!(rhs, vars, pars, iv, N, i, jac_times, jac_allocs, build_times, functions) outputs = rhs[1:N] SymbolicUtils.ENABLE_HASHCONSING[] = false jac_result = @be Symbolics.sparsejacobian(outputs, vars) @@ -99,46 +102,70 @@ function run_and_time!(rhs, vars, pars, iv, u0, p, t0, N, i, jac_times, jac_allo @assert isequal(jac_nohc, jac_hc) jac = jac_hc args = (vars, pars, iv) - kwargs = (; iip_config = (false, true), expression = Val{false}) + kwargs = (; iip_config = (false, true), expression = Val{true}) build_result = @be build_function(jac, args...; cse = false, kwargs...); build_times[1][i] = minimum(x -> x.time, build_result.samples) - jacfn_nocse = build_function(jac, args...; cse = false, kwargs...)[2] + jacfn_nocse = eval(build_function(jac, args...; cse = false, kwargs...)[2]) build_result = @be build_function(jac, args...; cse = true, kwargs...); build_times[2][i] = minimum(x -> x.time, build_result.samples) - jacfn_cse = build_function(jac, args...; cse = true, kwargs...)[2] + jacfn_cse = eval(build_function(jac, args...; cse = true, kwargs...)[2]) + + functions[1][i] = let buffer = similar(jac, Float64), fn = jacfn_nocse + function nocse(u, p, t) + fn(buffer, u, p, t) + buffer + end + end + functions[2][i] = let buffer = similar(jac, Float64), fn = jacfn_cse + function cse(u, p, t) + fn(buffer, u, p, t) + buffer + end + end + + return nothing +end + +function run_and_time_call!(i, u, p, tt, functions, first_call_times, second_call_times) + jacfn_nocse = functions[1][i] + jacfn_cse = functions[2][i] - buffer = similar(jac, Float64) - call_result = @timed jacfn_nocse(buffer, u0, p, t0) + call_result = @timed jacfn_nocse(u, p, tt) first_call_times[1][i] = call_result.time - call_result = @timed jacfn_cse(buffer, u0, p, t0) + call_result = @timed jacfn_cse(u, p, tt) first_call_times[2][i] = call_result.time - call_result = @be jacfn_nocse(buffer, u0, p, t0) + call_result = @be jacfn_nocse(u, p, tt) second_call_times[1][i] = minimum(x -> x.time, call_result.samples) - call_result = @be jacfn_cse(buffer, u0, p, t0) + call_result = @be jacfn_cse(u, p, tt) second_call_times[2][i] = minimum(x -> x.time, call_result.samples) - - return nothing end ``` # Run benchmark ```julia -N = [10, 20, 40, 60, 100, 200, 300, 400] +N = [10, 20, 40, 80] jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] jacobian_allocs = copy.(jacobian_times) +functions = [Vector{Any}(undef, length(N)), Vector{Any}(undef, length(N))] # [without_cse_times, with_cse_times] build_times = copy.(jacobian_times) first_call_times = copy.(jacobian_times) second_call_times = copy.(jacobian_times) iv = ModelingToolkit.get_iv(osys) -run_and_time!(rhs, vars, pars, iv, u, p, tt, 10, 1, jacobian_times, jacobian_allocs, build_times, first_call_times, second_call_times) +run_and_time_construct!(rhs, vars, pars, iv, 10, 1, jacobian_times, jacobian_allocs, build_times, functions) +run_and_time_call!(1, u, p, tt, functions, first_call_times, second_call_times) +for (i, n) in enumerate(N) + @info i n + run_and_time_construct!(rhs, vars, pars, iv, n, i, jacobian_times, jacobian_allocs, build_times, functions) +end for (i, n) in enumerate(N) - run_and_time!(rhs, vars, pars, iv, u, p, tt, n, i, jacobian_times, jacobian_allocs, build_times, first_call_times, second_call_times) + @info i n + run_and_time_call!(i, u, p, tt, functions, first_call_times, second_call_times) end ``` diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 820156126..388723ad2 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -251,9 +251,9 @@ function call(fn, args...) fn(args...) end -function run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, i, N) +function run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, functions, i, N) @mtkbuild sys = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058]) - @info "Built system" N + @info "Built system" SymbolicUtils.ENABLE_HASHCONSING[] = false jac_result = @be calculate_jacobian(sys; sparse = true) @info "No hashconsing benchmark" @@ -287,39 +287,58 @@ function run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_ buffer_cse = similar(jac, Float64) buffer_cse.nzval .= 0.0 - f_jac_nocse = build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false)[2] + f_jac_nocse = eval(build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{true}, cse = false)[2]) + functions[1][i] = let buffer_nocse = buffer_nocse, u0 = u0, p = p, t0 = t0, f_jac_nocse = f_jac_nocse + function nocse() + f_jac_nocse(buffer_nocse, u0, p, t0) + buffer_nocse + end + end @info "No CSE build_function result" - build_result_nocse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = false) + build_result_nocse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{true}, cse = false) @info "No CSE build_function benchmark" build_times[1][i] = mean(x -> x.time, build_result_nocse.samples) - @info "time" build_times[1][i] - first_call_result_nocse = @timed call(f_jac_nocse, buffer_nocse, u0, p, t0) - first_call_times[1][i] = first_call_result_nocse.time - @info "time" first_call_times[1][i] - second_call_result_nocse = @be call(f_jac_nocse, buffer_nocse, u0, p, t0) - second_call_times[1][i] = mean(x -> x.time, second_call_result_nocse.samples) - @info "time" second_call_times[1][i] - - f_jac_cse = build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true)[2] + @info "build_function time" build_times[1][i] + + f_jac_cse = eval(build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{true}, cse = true)[2]) + functions[2][i] = let buffer_cse = buffer_cse, u0 = u0, p = p, t0 = t0, f_jac_cse = f_jac_cse + function nocse() + f_jac_cse(buffer_cse, u0, p, t0) + buffer_cse + end + end @info "CSE build_function result" - build_result_cse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{false}, cse = true) + build_result_cse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{true}, cse = true) @info "CSE build_function benchmark" build_times[2][i] = mean(x -> x.time, build_result_cse.samples) - @info "time" build_times[2][i] - first_call_result_cse = @timed call(f_jac_cse, buffer_cse, u0, p, t0) - first_call_times[2][i] = first_call_result_cse.time - @info "time" first_call_times[2][i] - second_call_result_cse = @be call(f_jac_cse, buffer_cse, u0, p, t0) - second_call_times[2][i] = mean(x -> x.time, second_call_result_cse.samples) - @info "time" second_call_times[2][i] + @info "build_function time" build_times[2][i] return nothing end + +function run_and_time_call!(functions, first_call_times, second_call_times, i) + fnocse = functions[1][i] + fcse = functions[2][i] + first_call_result_nocse = @timed fnocse() + first_call_times[1][i] = first_call_result_nocse.time + @info "First call time" first_call_times[1][i] + second_call_result_nocse = @be fnocse() + second_call_times[1][i] = mean(x -> x.time, second_call_result_nocse.samples) + @info "Runtime" second_call_times[1][i] + + first_call_result_cse = @timed fcse() + first_call_times[2][i] = first_call_result_cse.time + @info "First call time" first_call_times[2][i] + second_call_result_cse = @be fcse() + second_call_times[2][i] = mean(x -> x.time, second_call_result_cse.samples) + @info "Runtime" second_call_times[2][i] +end ``` ```julia -N = [5, 10, 20, 40, 80, 160, 320]; +N = [5, 10, 20, 40]; jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] +functions = [Vector{Any}(undef, length(N)), Vector{Any}(undef, length(N))] jacobian_gctimes = copy.(jacobian_times) jacobian_allocs = copy.(jacobian_times) # [without_cse_times, with_cse_times] @@ -332,10 +351,16 @@ second_call_times = copy.(jacobian_times) ```julia # compile -run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, 1, 5) +run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, functions, 1, 5) +run_and_time_call!(functions, first_call_times, second_call_times, 1) +for (i, n) in enumerate(N) + @info i n + @time run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, functions, i, n) +end + for (i, n) in enumerate(N) @info i n - @time run_and_time!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times, i, n) + run_and_time_call!(functions, first_call_times, second_call_times, i) end ``` From 8220e4250ab12913109224a367560bc7a5f70fd9 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 16:38:39 +0530 Subject: [PATCH 17/31] fix: remove debugging --- benchmarks/Symbolics/BCR.jmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index b98ce3c87..f740fc3f2 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -33,7 +33,7 @@ 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)][1:10] +rhs = [eq.rhs for eq in full_equations(osys)] vars = unknowns(osys) pars = parameters(osys) From 8f9013dac9e3b773ce10451ff49cd99066244d04 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 16:39:19 +0530 Subject: [PATCH 18/31] fix: run for more sizes --- benchmarks/Symbolics/BCR.jmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index f740fc3f2..07e9e6a85 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -147,7 +147,7 @@ end # Run benchmark ```julia -N = [10, 20, 40, 80] +N = [10, 20, 40, 80, 160, 320] jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] jacobian_allocs = copy.(jacobian_times) functions = [Vector{Any}(undef, length(N)), Vector{Any}(undef, length(N))] From be70a601c9013af41f067142aad38196972180ad Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 16:39:42 +0530 Subject: [PATCH 19/31] fix: run for more sizes --- benchmarks/Symbolics/ThermalFluid.jmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 388723ad2..f58f6d501 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -336,7 +336,7 @@ end ``` ```julia -N = [5, 10, 20, 40]; +N = [5, 10, 20, 40, 80, 160, 320]; jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] functions = [Vector{Any}(undef, length(N)), Vector{Any}(undef, length(N))] jacobian_gctimes = copy.(jacobian_times) From ca6d6742da549db47ea54b2eebdd968565ab57db Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 25 Mar 2025 00:19:45 +0530 Subject: [PATCH 20/31] fix: clear `occursin_info` cache in `@be` macro --- benchmarks/Symbolics/BCR.jmd | 4 ++-- benchmarks/Symbolics/ThermalFluid.jmd | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 07e9e6a85..2426d0216 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -93,8 +93,7 @@ function run_and_time_construct!(rhs, vars, pars, iv, N, i, jac_times, jac_alloc jac_nohc = Symbolics.sparsejacobian(outputs, vars) SymbolicUtils.ENABLE_HASHCONSING[] = true - SymbolicUtils.clear_cache!(Symbolics.occursin_info) - jac_result = @be Symbolics.sparsejacobian(outputs, vars) + jac_result = @be (SymbolicUtils.clear_cache!(Symbolics.occursin_info); Symbolics.sparsejacobian(outputs, vars)) jac_times[2][i] = minimum(x -> x.time, jac_result.samples) jac_allocs[2][i] = minimum(x -> x.bytes, jac_result.samples) jac_hc = Symbolics.sparsejacobian(outputs, vars) @@ -147,6 +146,7 @@ end # Run benchmark ```julia +Chairmarks.DEFAULTS.seconds = 15.0 N = [10, 20, 40, 80, 160, 320] jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] jacobian_allocs = copy.(jacobian_times) diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index f58f6d501..0ad60ab5f 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -264,8 +264,7 @@ function run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_a jacobian_allocs[1][i] = mean(x -> x.bytes, jac_result.samples) @info "times" jacobian_times[1][i] jacobian_gctimes[1][i] jacobian_allocs[1][i] SymbolicUtils.ENABLE_HASHCONSING[] = true - SymbolicUtils.clear_cache!(Symbolics.occursin_info) - jac_result = @be calculate_jacobian(sys; sparse = true) + jac_result = @be (SymbolicUtils.clear_cache!(Symbolics.occursin_info); calculate_jacobian(sys; sparse = true)) @info "Hashconsing benchmark" jac_cse = calculate_jacobian(sys; sparse = true) @info "Hashconsing result" @@ -350,6 +349,7 @@ second_call_times = copy.(jacobian_times) ## Timings ```julia +Chairmarks.DEFAULTS.seconds = 15.0 # compile run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, functions, 1, 5) run_and_time_call!(functions, first_call_times, second_call_times, 1) From dcf628eee4be2848fcfb2c549c05b62f166b366b Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 26 Mar 2025 23:50:17 +0530 Subject: [PATCH 21/31] refactor: rewrite benchmarks to account for new derivative calculation --- benchmarks/Symbolics/BCR.jmd | 37 ++-- benchmarks/Symbolics/ThermalFluid.jmd | 20 +- benchmarks/Symbolics/old_sparse_jacobian.jl | 227 ++++++++++++++++++++ 3 files changed, 266 insertions(+), 18 deletions(-) create mode 100644 benchmarks/Symbolics/old_sparse_jacobian.jl diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 2426d0216..07c76ab4d 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -17,7 +17,7 @@ jacobian, generate a function to calculate it and call the function. ```julia -using OrdinaryDiffEq, Catalyst, ReactionNetworkImporters, +using Catalyst, ReactionNetworkImporters, TimerOutputs, LinearAlgebra, ModelingToolkit, Chairmarks, LinearSolve, Symbolics, SymbolicUtils, SymbolicUtils.Code, SparseArrays, CairoMakie, PrettyTables @@ -36,15 +36,29 @@ osys = convert(ODESystem, rn) rhs = [eq.rhs for eq in full_equations(osys)] vars = unknowns(osys) pars = parameters(osys) +``` + +The `sparsejacobian` function in Symbolics.jl is optimized for hashconsing and caching, and as such +performs very poorly without either of those features. We use the old implementation, optimized without +hashconsing, to benchmark performance without hashconsing and without caching to avoid biasing the results. +```julia +include("old_sparse_jacobian.jl") +``` + +```julia SymbolicUtils.ENABLE_HASHCONSING[] = false -@timeit to "Calculate jacobian - without hashconsing" jac_nohc = Symbolics.sparsejacobian(rhs, vars); +@timeit to "Calculate jacobian - without hashconsing" jac_nohc = old_sparsejacobian(rhs, vars); SymbolicUtils.ENABLE_HASHCONSING[] = true -SymbolicUtils.toggle_caching!(Symbolics.occursin_info, false) -@timeit to "Calculate jacobian - hashconsing, without caching" jac_hc_nocache = Symbolics.sparsejacobian(rhs, vars); -SymbolicUtils.toggle_caching!(Symbolics.occursin_info, true) -stats = SymbolicUtils.get_stats(Symbolics.occursin_info) -@assert stats.hits == stats.misses == 0 +Symbolics.toggle_derivative_caching!(false) +Symbolics.clear_derivative_caches!() +@timeit to "Calculate jacobian - hashconsing, without caching" jac_hc_nocache = old_sparsejacobian(rhs, vars); +Symbolics.toggle_derivative_caching!(true) +for fn in Symbolics.cached_derivative_functions() + stats = SymbolicUtils.get_stats(fn) + @assert stats.hits == stats.misses == 0 +end +Symbolics.clear_derivative_caches!() @timeit to "Calculate jacobian - hashconsing and caching" jac_hc_cache = Symbolics.sparsejacobian(rhs, vars); @assert isequal(jac_nohc, jac_hc_nocache) @@ -87,19 +101,16 @@ We'll also measure scaling. function run_and_time_construct!(rhs, vars, pars, iv, N, i, jac_times, jac_allocs, build_times, functions) outputs = rhs[1:N] SymbolicUtils.ENABLE_HASHCONSING[] = false - jac_result = @be Symbolics.sparsejacobian(outputs, vars) + jac_result = @be old_sparsejacobian(outputs, vars) jac_times[1][i] = minimum(x -> x.time, jac_result.samples) jac_allocs[1][i] = minimum(x -> x.bytes, jac_result.samples) - jac_nohc = Symbolics.sparsejacobian(outputs, vars) SymbolicUtils.ENABLE_HASHCONSING[] = true - jac_result = @be (SymbolicUtils.clear_cache!(Symbolics.occursin_info); Symbolics.sparsejacobian(outputs, vars)) + jac_result = @be (Symbolics.clear_derivative_caches!(); Symbolics.sparsejacobian(outputs, vars)) jac_times[2][i] = minimum(x -> x.time, jac_result.samples) jac_allocs[2][i] = minimum(x -> x.bytes, jac_result.samples) - jac_hc = Symbolics.sparsejacobian(outputs, vars) - @assert isequal(jac_nohc, jac_hc) - jac = jac_hc + jac = Symbolics.sparsejacobian(outputs, vars) args = (vars, pars, iv) kwargs = (; iip_config = (false, true), expression = Val{true}) diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 0ad60ab5f..aa2c94c54 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -250,23 +250,34 @@ end function call(fn, args...) fn(args...) end +``` + +The `sparsejacobian` function in Symbolics.jl is optimized for hashconsing and caching, and as such +performs very poorly without either of those features. We use the old implementation, optimized without +hashconsing, to benchmark performance without hashconsing and without caching to avoid biasing the results. + +```julia +include("old_sparse_jacobian.jl") function run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, functions, i, N) @mtkbuild sys = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058]) + rhs = [eq.rhs for eq in full_equations(sys)] + dvs = unknowns(sys) + @info "Built system" SymbolicUtils.ENABLE_HASHCONSING[] = false - jac_result = @be calculate_jacobian(sys; sparse = true) + jac_result = @be old_sparsejacobian(rhs, dvs) @info "No hashconsing benchmark" - jac_nocse = calculate_jacobian(sys; sparse = true) + jac_nocse = old_sparsejacobian(rhs, dvs) @info "No hashconsing result" jacobian_times[1][i] = mean(x -> x.time, jac_result.samples) jacobian_gctimes[1][i] = mean(x -> x.time * x.gc_fraction, jac_result.samples) jacobian_allocs[1][i] = mean(x -> x.bytes, jac_result.samples) @info "times" jacobian_times[1][i] jacobian_gctimes[1][i] jacobian_allocs[1][i] SymbolicUtils.ENABLE_HASHCONSING[] = true - jac_result = @be (SymbolicUtils.clear_cache!(Symbolics.occursin_info); calculate_jacobian(sys; sparse = true)) + jac_result = @be (Symbolics.clear_derivative_caches!(); Symbolics.sparsejacobian(rhs, dvs)) @info "Hashconsing benchmark" - jac_cse = calculate_jacobian(sys; sparse = true) + jac_cse = Symbolics.sparsejacobian(rhs, dvs) @info "Hashconsing result" jacobian_times[2][i] = mean(x -> x.time, jac_result.samples) jacobian_gctimes[2][i] = mean(x -> x.time * x.gc_fraction, jac_result.samples) @@ -275,7 +286,6 @@ function run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_a @assert isequal(jac_nocse, jac_cse) jac = jac_cse - dvs = unknowns(sys) ps = parameters(sys) defs = defaults(sys) u0 = Float64[Symbolics.fixpoint_sub(v, defs) for v in dvs] diff --git a/benchmarks/Symbolics/old_sparse_jacobian.jl b/benchmarks/Symbolics/old_sparse_jacobian.jl new file mode 100644 index 000000000..9a009473c --- /dev/null +++ b/benchmarks/Symbolics/old_sparse_jacobian.jl @@ -0,0 +1,227 @@ +function old_sparsejacobian(ops::AbstractVector, vars::AbstractVector) + sp = Symbolics.jacobian_sparsity(ops, vars) + I,J,_ = findnz(sp) + + exprs = old_sparsejacobian_vals(ops, vars, I, J) + + sparse(I, J, exprs, length(ops), length(vars)) +end + +function old_sparsejacobian_vals(ops::AbstractVector, vars::AbstractVector, I::AbstractVector, J::AbstractVector; simplify::Bool=false, kwargs...) + exprs = Num[] + sizehint!(exprs, length(I)) + + for (i,j) in zip(I, J) + push!(exprs, Num(old_expand_derivatives(Differential(vars[j])(ops[i]), simplify; kwargs...))) + end + exprs +end + + +function old_expand_derivatives(O::SymbolicUtils.Symbolic, simplify=false; throw_no_derivative=false) + if iscall(O) && isa(operation(O), Differential) + arg = only(arguments(O)) + arg = old_expand_derivatives(arg, false; throw_no_derivative) + return old_executediff(operation(O), arg, simplify; throw_no_derivative) + elseif iscall(O) && isa(operation(O), Integral) + return operation(O)(old_expand_derivatives(arguments(O)[1]; throw_no_derivative)) + elseif !Symbolics.hasderiv(O) + return O + else + args = map(a->old_expand_derivatives(a, false; throw_no_derivative), arguments(O)) + O1 = operation(O)(args...) + return simplify ? SymbolicUtils.simplify(O1) : O1 + end +end +function old_expand_derivatives(n::Num, simplify=false; kwargs...) + Symbolics.wrap(old_expand_derivatives(Symbolics.value(n), simplify; kwargs...)) +end + +function old_occursin_info(x, expr, fail = true) + if SymbolicUtils.symtype(expr) <: AbstractArray + if fail + error("Differentiation with array expressions is not yet supported") + else + return occursin(x, expr) + end + end + + # Allow scalarized expressions + function is_scalar_indexed(ex) + (iscall(ex) && operation(ex) == getindex && !(SymbolicUtils.symtype(ex) <: AbstractArray)) || + (iscall(ex) && (SymbolicUtils.issym(operation(ex)) || iscall(operation(ex))) && + is_scalar_indexed(operation(ex))) + end + + # x[1] == x[1] but not x[2] + if is_scalar_indexed(x) && is_scalar_indexed(expr) && + isequal(first(arguments(x)), first(arguments(expr))) + return isequal(operation(x), operation(expr)) && + isequal(arguments(x), arguments(expr)) + end + + if is_scalar_indexed(x) && is_scalar_indexed(expr) && + !occursin(first(arguments(x)), first(arguments(expr))) + return false + end + + if is_scalar_indexed(expr) && !is_scalar_indexed(x) && !occursin(x, expr) + return false + end + + !iscall(expr) && return isequal(x, expr) + if isequal(x, expr) + true + else + args = map(a->old_occursin_info(x, a, operation(expr) !== getindex), arguments(expr)) + if all(_isfalse, args) + return false + end + Term{Real}(true, args) + end +end + +function old_occursin_info(x, expr::Sym, fail) + if SymbolicUtils.symtype(expr) <: AbstractArray && fail + error("Differentiation of expressions involving arrays and array variables is not yet supported.") + end + isequal(x, expr) +end + +_isfalse(occ::Bool) = occ === false +_isfalse(occ::SymbolicUtils.Symbolic) = iscall(occ) && _isfalse(operation(occ)) + +_iszero(x) = false +_isone(x) = false +_iszero(x::Number) = iszero(x) +_isone(x::Number) = isone(x) +_iszero(::SymbolicUtils.Symbolic) = false +_isone(::SymbolicUtils.Symbolic) = false +_iszero(x::Num) = _iszero(value(x))::Bool +_isone(x::Num) = _isone(value(x))::Bool + + +function old_executediff(D, arg, simplify=false; occurrences=nothing, throw_no_derivative=false) + if occurrences == nothing + occurrences = old_occursin_info(D.x, arg) + end + + _isfalse(occurrences) && return 0 + occurrences isa Bool && return 1 # means it's a `true` + + if !iscall(arg) + return D(arg) # Cannot expand + elseif (op = operation(arg); SymbolicUtils.issym(op)) + inner_args = arguments(arg) + if any(isequal(D.x), inner_args) + return D(arg) # base case if any argument is directly equal to the i.v. + else + return sum(inner_args, init=0) do a + return old_executediff(Differential(a), arg; throw_no_derivative) * + old_executediff(D, a; throw_no_derivative) + end + end + elseif op === getindex + inner_args = arguments(arguments(arg)[1]) + c = 0 + for a in inner_args + if isequal(a, D.x) + return D(arg) + else + c += Differential(a)(arg) * D(a) + end + end + return old_expand_derivatives(c) + elseif op === ifelse + args = arguments(arg) + O = op(args[1], + old_executediff(D, args[2], simplify; occurrences=arguments(occurrences)[2], throw_no_derivative), + old_executediff(D, args[3], simplify; occurrences=arguments(occurrences)[3], throw_no_derivative)) + return O + elseif isa(op, Differential) + # The recursive expand_derivatives was not able to remove + # a nested Differential. We can attempt to differentiate the + # inner expression wrt to the outer iv. And leave the + # unexpandable Differential outside. + if isequal(op.x, D.x) + return D(arg) + else + inner = old_executediff(D, arguments(arg)[1], false; throw_no_derivative) + # if the inner expression is not expandable either, return + if iscall(inner) && operation(inner) isa Differential + return D(arg) + else + # otherwise give the nested Differential another try + return old_executediff(op, inner, simplify; throw_no_derivative) + end + end + elseif isa(op, Integral) + if isa(op.domain.domain, Symbolics.AbstractInterval) + domain = op.domain.domain + a, b = Symbolics.DomainSets.endpoints(domain) + c = 0 + inner_function = arguments(arg)[1] + if iscall(value(a)) + t1 = SymbolicUtils.substitute(inner_function, Dict(op.domain.variables => value(a))) + t2 = D(a) + c -= t1*t2 + end + if iscall(value(b)) + t1 = SymbolicUtils.substitute(inner_function, Dict(op.domain.variables => value(b))) + t2 = D(b) + c += t1*t2 + end + inner = old_executediff(D, arguments(arg)[1]; throw_no_derivative) + c += op(inner) + return Symbolics.value(c) + end + end + + inner_args = arguments(arg) + l = length(inner_args) + exprs = [] + c = 0 + + for i in 1:l + t2 = old_executediff(D, inner_args[i],false; occurrences=arguments(occurrences)[i], throw_no_derivative) + + x = if _iszero(t2) + t2 + elseif _isone(t2) + d = Symbolics.derivative_idx(arg, i) + if d isa Symbolics.NoDeriv + throw_no_derivative && error((arg, i)) + D(arg) + else + d + end + else + t1 = Symbolics.derivative_idx(arg, i) + t1 = if t1 isa Symbolics.NoDeriv + throw_no_derivative && error((arg, i)) + D(arg) + else + t1 + end + t1 * t2 + end + + if _iszero(x) + continue + elseif x isa SymbolicUtils.Symbolic + push!(exprs, x) + else + c += x + end + end + + if isempty(exprs) + return c + elseif length(exprs) == 1 + term = (simplify ? SymbolicUtils.simplify(exprs[1]) : exprs[1]) + return _iszero(c) ? term : c + term + else + x = +((!_iszero(c) ? vcat(c, exprs) : exprs)...) + return simplify ? SymbolicUtils.simplify(x) : x + end +end From ae3f7b63c115e402d86e834889a850ce75fc14a4 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 26 Mar 2025 23:50:33 +0530 Subject: [PATCH 22/31] build: bump Symbolics compat --- benchmarks/Symbolics/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml index d882ab18f..53acd2cd7 100644 --- a/benchmarks/Symbolics/Project.toml +++ b/benchmarks/Symbolics/Project.toml @@ -22,4 +22,4 @@ XSteam = "95ff35a0-be81-11e9-2ca3-5b4e338e8476" [compat] ModelingToolkit = "9" SymbolicUtils = "3.21" -Symbolics = "6.33.1" +Symbolics = "6.34" From 2b70c8b4df2b8abf1b0ad5ffbc261b18594d05a5 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 3 Apr 2025 00:13:33 +0530 Subject: [PATCH 23/31] TEMP run benchmark with uint-id --- benchmarks/Symbolics/BCR.jmd | 6 +++++- benchmarks/Symbolics/ThermalFluid.jmd | 1 + 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 07c76ab4d..5edffe293 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -17,10 +17,14 @@ jacobian, generate a function to calculate it and call the function. ```julia +using Pkg + +Pkg.add(Pkg.PackageSpec(;name="SymbolicUtils", rev="775efa7556f986087562f2a32873fa1f7bcf173c")) + using Catalyst, ReactionNetworkImporters, TimerOutputs, LinearAlgebra, ModelingToolkit, Chairmarks, LinearSolve, Symbolics, SymbolicUtils, SymbolicUtils.Code, SparseArrays, CairoMakie, - PrettyTables + PrettyTables, NoImport datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index aa2c94c54..86e2e2ae3 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -9,6 +9,7 @@ This is a 1D advection-diffusion-source PDE that uses a second order upwind sche using Pkg # Rev fixes precompilation https://github.com/hzgzh/XSteam.jl/pull/2 Pkg.add(Pkg.PackageSpec(;name="XSteam", rev="f2a1c589054cfd6bba307985a3a534b6f5a1863b")) +Pkg.add(Pkg.PackageSpec(;name="SymbolicUtils", rev="775efa7556f986087562f2a32873fa1f7bcf173c")) using ModelingToolkit, Symbolics, SymbolicUtils, XSteam, Polynomials, CairoMakie, PrettyTables using SparseArrays, Chairmarks, Statistics From a7b8ba913896c3f1a283cfb414d617bbdaeb7860 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 3 Apr 2025 11:35:42 +0530 Subject: [PATCH 24/31] fixup! TEMP run benchmark with uint-id --- benchmarks/Symbolics/BCR.jmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 5edffe293..56d033113 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -24,7 +24,7 @@ Pkg.add(Pkg.PackageSpec(;name="SymbolicUtils", rev="775efa7556f986087562f2a32873 using Catalyst, ReactionNetworkImporters, TimerOutputs, LinearAlgebra, ModelingToolkit, Chairmarks, LinearSolve, Symbolics, SymbolicUtils, SymbolicUtils.Code, SparseArrays, CairoMakie, - PrettyTables, NoImport + PrettyTables datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() From fb66afbc833dcaa146f08472040aa99cfba178b4 Mon Sep 17 00:00:00 2001 From: "Bowen S. Zhu" Date: Wed, 9 Apr 2025 15:23:36 -0400 Subject: [PATCH 25/31] refactor: update BCR benchmark visualization with improved axes and labels --- benchmarks/Symbolics/BCR.jmd | 54 ++++++++++++++++++++++++++++++------ 1 file changed, 46 insertions(+), 8 deletions(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 56d033113..113d41e7d 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -193,15 +193,53 @@ pretty_table(tabledata; header, backend = Val(:html)) ``` ```julia -f = Figure(size = (800, 1200)); -labels = ["Symbolic jacobian time", "Symbolic jacobian\nallocated memory (B)", "`build_function` time", "First call time", "Second call time"] +f = Figure(size = (750, 400)) +titles = [ + "Jacobian symbolic computation", "Jacobian symbolic computation", "Code generation", + "Numerical function compilation", "Numerical function evaluation"] +labels = ["Time (seconds)", "Allocated memory (bytes)", + "Time (seconds)", "Time (seconds)", "Time (seconds)"] times = [jacobian_times, jacobian_allocs, build_times, first_call_times, second_call_times] -for (i, (label, data)) in enumerate(zip(labels, times)) - ax = Axis(f[i, 1], xscale = log10, yscale = log10, ylabel = label, xticks = N) - l1 = lines!(ax, N, data[1]) - l2 = lines!(ax, N, data[2]) - legend_entries = startswith(label, "Symbolic") ? ["without hashconsing", "with hashconsing"] : ["without CSE", "with CSE"] - Legend(f[i, 2], [l1, l2], legend_entries) +axes = Axis[] +for i in 1:2 + label = labels[i] + data = times[i] + ax = Axis(f[1, i], xscale = log10, yscale = log10, xlabel = "model size", + xlabelsize = 10, ylabel = label, ylabelsize = 10, xticks = N, + title = titles[i], titlesize = 12, xticklabelsize = 10, yticklabelsize = 10) + push!(axes, ax) + l1 = scatterlines!(ax, N, data[1], label = "without hashconsing") + l2 = scatterlines!(ax, N, data[2], label = "with hashconsing") end +Legend(f[1, 3], axes[1], "Methods", tellwidth = false, labelsize = 12, titlesize = 15) +axes2 = Axis[] +# make equal y-axis unit length +mn3, mx3 = extrema(reduce(vcat, times[3])) +xn3 = log10(mx3 / mn3) +mn4, mx4 = extrema(reduce(vcat, times[4])) +xn4 = log10(mx4 / mn4) +mn5, mx5 = extrema(reduce(vcat, times[5])) +xn5 = log10(mx5 / mn5) +xn = max(xn3, xn4, xn5) +xn += 0.2 +hxn = xn / 2 +hxn3 = (log10(mx3) + log10(mn3)) / 2 +hxn4 = (log10(mx4) + log10(mn4)) / 2 +hxn5 = (log10(mx5) + log10(mn5)) / 2 +ylims = [(exp10(hxn3 - hxn), exp10(hxn3 + hxn)), (exp10(hxn4 - hxn), exp10(hxn4 + hxn)), + (exp10(hxn5 - hxn), exp10(hxn5 + hxn))] +for i in 1:3 + ir = i + 2 + label = labels[ir] + data = times[ir] + ax = Axis(f[2, i], xscale = log10, yscale = log10, xlabel = "model size", + xlabelsize = 10, ylabel = label, ylabelsize = 10, xticks = N, + title = titles[ir], titlesize = 12, xticklabelsize = 10, yticklabelsize = 10) + ylims!(ax, ylims[i]...) + push!(axes2, ax) + l1 = scatterlines!(ax, N, data[1], label = "without hashconsing") + l2 = scatterlines!(ax, N, data[2], label = "with hashconsing") +end +save("bcr.pdf", f) f ``` From c6e1468dcb5b7d7cc29b3913e669886ad48b7d50 Mon Sep 17 00:00:00 2001 From: "Bowen S. Zhu" Date: Wed, 9 Apr 2025 15:26:03 -0400 Subject: [PATCH 26/31] refactor: enhance ThermalFluid benchmark visualization with improved layout and labels --- benchmarks/Symbolics/ThermalFluid.jmd | 56 ++++++++++++++++++++++----- 1 file changed, 47 insertions(+), 9 deletions(-) diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 86e2e2ae3..3dc1f4e51 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -384,16 +384,54 @@ pretty_table(tabledata; header, backend = Val(:html)) ``` ```julia -f = Figure(size = (800, 1200)); -labels = ["Symbolic jacobian time", "Symbolic jacobian\nGC time", "Symbolic jacobian\nallocated memory (B)", "`build_function` time", "First call time", "Second call time"] -times = [jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, first_call_times, second_call_times] -for (i, (label, data)) in enumerate(zip(labels, times)) - ax = Axis(f[i, 1], xscale = log10, yscale = log10, ylabel = label, xticks = N) - l1 = lines!(ax, N, data[1]) - l2 = lines!(ax, N, data[2]) - legend_entries = startswith(label, "Symbolic") ? ["without hashconsing", "with hashconsing"] : ["without CSE", "with CSE"] - Legend(f[i, 2], [l1, l2], legend_entries) +f = Figure(size = (750, 400)) +titles = [ + "Jacobian symbolic computation", "Jacobian symbolic computation", "Code generation", + "Numerical function compilation", "Numerical function evaluation"] +labels = ["Time (seconds)", "Allocated memory (bytes)", + "Time (seconds)", "Time (seconds)", "Time (seconds)"] +times = [jacobian_times, jacobian_allocs, build_times, first_call_times, second_call_times] +axes = Axis[] +for i in 1:2 + label = labels[i] + data = times[i] + ax = Axis(f[1, i], xscale = log10, yscale = log10, xlabel = "model size", + xlabelsize = 10, ylabel = label, ylabelsize = 10, xticks = N, + title = titles[i], titlesize = 12, xticklabelsize = 10, yticklabelsize = 10) + push!(axes, ax) + l1 = scatterlines!(ax, N, data[1], label = "without hashconsing") + l2 = scatterlines!(ax, N, data[2], label = "with hashconsing") end +Legend(f[1, 3], axes[1], "Methods", tellwidth = false, labelsize = 12, titlesize = 15) +axes2 = Axis[] +# make equal y-axis unit length +mn3, mx3 = extrema(reduce(vcat, times[3])) +xn3 = log10(mx3 / mn3) +mn4, mx4 = extrema(reduce(vcat, times[4])) +xn4 = log10(mx4 / mn4) +mn5, mx5 = extrema(reduce(vcat, times[5])) +xn5 = log10(mx5 / mn5) +xn = max(xn3, xn4, xn5) +xn += 0.2 +hxn = xn / 2 +hxn3 = (log10(mx3) + log10(mn3)) / 2 +hxn4 = (log10(mx4) + log10(mn4)) / 2 +hxn5 = (log10(mx5) + log10(mn5)) / 2 +ylims = [(exp10(hxn3 - hxn), exp10(hxn3 + hxn)), (exp10(hxn4 - hxn), exp10(hxn4 + hxn)), + (exp10(hxn5 - hxn), exp10(hxn5 + hxn))] +for i in 1:3 + ir = i + 2 + label = labels[ir] + data = times[ir] + ax = Axis(f[2, i], xscale = log10, yscale = log10, xlabel = "model size", + xlabelsize = 10, ylabel = label, ylabelsize = 10, xticks = N, + title = titles[ir], titlesize = 12, xticklabelsize = 10, yticklabelsize = 10) + ylims!(ax, ylims[i]...) + push!(axes2, ax) + l1 = scatterlines!(ax, N, data[1], label = "without hashconsing") + l2 = scatterlines!(ax, N, data[2], label = "with hashconsing") +end +save("thermal_fluid.pdf", f) f ``` From 8287ac4fb8ab0d5f8442eb4cce08ff32029da68c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 14 Apr 2025 16:32:46 +0530 Subject: [PATCH 27/31] build: bump SymbolicUtils compat --- benchmarks/Symbolics/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml index 53acd2cd7..f434b6b74 100644 --- a/benchmarks/Symbolics/Project.toml +++ b/benchmarks/Symbolics/Project.toml @@ -21,5 +21,5 @@ XSteam = "95ff35a0-be81-11e9-2ca3-5b4e338e8476" [compat] ModelingToolkit = "9" -SymbolicUtils = "3.21" +SymbolicUtils = "3.26" Symbolics = "6.34" From 793306a0916f3db0d18415edb5b589c0c2ccb410 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 14 Apr 2025 16:33:42 +0530 Subject: [PATCH 28/31] fix: remove usage of non-primary SymbolicUtils branch --- benchmarks/Symbolics/BCR.jmd | 4 ---- benchmarks/Symbolics/ThermalFluid.jmd | 1 - 2 files changed, 5 deletions(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 113d41e7d..a717190ad 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -17,10 +17,6 @@ jacobian, generate a function to calculate it and call the function. ```julia -using Pkg - -Pkg.add(Pkg.PackageSpec(;name="SymbolicUtils", rev="775efa7556f986087562f2a32873fa1f7bcf173c")) - using Catalyst, ReactionNetworkImporters, TimerOutputs, LinearAlgebra, ModelingToolkit, Chairmarks, LinearSolve, Symbolics, SymbolicUtils, SymbolicUtils.Code, SparseArrays, CairoMakie, diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 3dc1f4e51..bf7031a7c 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -9,7 +9,6 @@ This is a 1D advection-diffusion-source PDE that uses a second order upwind sche using Pkg # Rev fixes precompilation https://github.com/hzgzh/XSteam.jl/pull/2 Pkg.add(Pkg.PackageSpec(;name="XSteam", rev="f2a1c589054cfd6bba307985a3a534b6f5a1863b")) -Pkg.add(Pkg.PackageSpec(;name="SymbolicUtils", rev="775efa7556f986087562f2a32873fa1f7bcf173c")) using ModelingToolkit, Symbolics, SymbolicUtils, XSteam, Polynomials, CairoMakie, PrettyTables using SparseArrays, Chairmarks, Statistics From e4bd0fa789227bc7d9f80869d0194dfd8658a2c4 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 17 Apr 2025 19:03:27 +0530 Subject: [PATCH 29/31] TEMP: use branch of SymbolicUtils --- benchmarks/Symbolics/BCR.jmd | 2 ++ benchmarks/Symbolics/ThermalFluid.jmd | 1 + 2 files changed, 3 insertions(+) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index a717190ad..82c217f6d 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -22,6 +22,8 @@ using Catalyst, ReactionNetworkImporters, LinearSolve, Symbolics, SymbolicUtils, SymbolicUtils.Code, SparseArrays, CairoMakie, PrettyTables +Pkg.add(Pkg.PackageSpec(; name = "SymbolicUtils", rev="9d52810c547276a60c0aae502ffee7f93ff99143") + datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() tf = 100000.0 diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index bf7031a7c..9dff2c52b 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -9,6 +9,7 @@ This is a 1D advection-diffusion-source PDE that uses a second order upwind sche using Pkg # Rev fixes precompilation https://github.com/hzgzh/XSteam.jl/pull/2 Pkg.add(Pkg.PackageSpec(;name="XSteam", rev="f2a1c589054cfd6bba307985a3a534b6f5a1863b")) +Pkg.add(Pkg.PackageSpec(; name = "SymbolicUtils", rev="9d52810c547276a60c0aae502ffee7f93ff99143") using ModelingToolkit, Symbolics, SymbolicUtils, XSteam, Polynomials, CairoMakie, PrettyTables using SparseArrays, Chairmarks, Statistics From ca0924ad0d90847c1096f9ce33c5077ac25224fb Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 17 Apr 2025 20:20:51 +0530 Subject: [PATCH 30/31] fixup! TEMP: use branch of SymbolicUtils --- benchmarks/Symbolics/BCR.jmd | 2 +- benchmarks/Symbolics/ThermalFluid.jmd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 82c217f6d..6e75d6130 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -22,7 +22,7 @@ using Catalyst, ReactionNetworkImporters, LinearSolve, Symbolics, SymbolicUtils, SymbolicUtils.Code, SparseArrays, CairoMakie, PrettyTables -Pkg.add(Pkg.PackageSpec(; name = "SymbolicUtils", rev="9d52810c547276a60c0aae502ffee7f93ff99143") +Pkg.add(Pkg.PackageSpec(; name = "SymbolicUtils", rev="9d52810c547276a60c0aae502ffee7f93ff99143")) datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd index 9dff2c52b..22066a25a 100644 --- a/benchmarks/Symbolics/ThermalFluid.jmd +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -9,7 +9,7 @@ This is a 1D advection-diffusion-source PDE that uses a second order upwind sche using Pkg # Rev fixes precompilation https://github.com/hzgzh/XSteam.jl/pull/2 Pkg.add(Pkg.PackageSpec(;name="XSteam", rev="f2a1c589054cfd6bba307985a3a534b6f5a1863b")) -Pkg.add(Pkg.PackageSpec(; name = "SymbolicUtils", rev="9d52810c547276a60c0aae502ffee7f93ff99143") +Pkg.add(Pkg.PackageSpec(; name = "SymbolicUtils", rev="9d52810c547276a60c0aae502ffee7f93ff99143")) using ModelingToolkit, Symbolics, SymbolicUtils, XSteam, Polynomials, CairoMakie, PrettyTables using SparseArrays, Chairmarks, Statistics From 2f4644699ba37c025f3ff00ef3a6a0752c2386e9 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 18 Apr 2025 13:21:34 +0530 Subject: [PATCH 31/31] fixup! TEMP: use branch of SymbolicUtils --- benchmarks/Symbolics/BCR.jmd | 1 + 1 file changed, 1 insertion(+) diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd index 6e75d6130..ad229da89 100644 --- a/benchmarks/Symbolics/BCR.jmd +++ b/benchmarks/Symbolics/BCR.jmd @@ -22,6 +22,7 @@ using Catalyst, ReactionNetworkImporters, LinearSolve, Symbolics, SymbolicUtils, SymbolicUtils.Code, SparseArrays, CairoMakie, PrettyTables +using Pkg Pkg.add(Pkg.PackageSpec(; name = "SymbolicUtils", rev="9d52810c547276a60c0aae502ffee7f93ff99143")) datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr")