diff --git a/docs/src/inverse_problems/optimization_ode_param_fitting.md b/docs/src/inverse_problems/optimization_ode_param_fitting.md index 8f0d530e0c..c579a68aa7 100644 --- a/docs/src/inverse_problems/optimization_ode_param_fitting.md +++ b/docs/src/inverse_problems/optimization_ode_param_fitting.md @@ -1,11 +1,11 @@ # [Parameter Fitting for ODEs using SciML/Optimization.jl and DiffEqParamEstim.jl](@id optimization_parameter_fitting) Fitting parameters to data involves solving an optimisation problem (that is, finding the parameter set that optimally fits your model to your data, typically by minimising a cost function). The SciML ecosystem's primary package for solving optimisation problems is [Optimization.jl](https://github.com/SciML/Optimization.jl). It provides access to a variety of solvers via a single common interface by wrapping a large number of optimisation libraries that have been implemented in Julia. -This tutorial demonstrates both how to create parameter fitting cost functions using the [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl) package, and how to use Optimization.jl to minimise these. Optimization.jl can also be used in other contexts, such as finding parameter sets that maximise the magnitude of some system behaviour. More details on how to use these packages can be found in their [respective](https://docs.sciml.ai/Optimization/stable/) [documentations](https://docs.sciml.ai/DiffEqParamEstim/stable/). +This tutorial demonstrates both how to create parameter fitting cost functions using the [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl) package, and how to use Optimization.jl to minimise these. Optimization.jl can also be used in other contexts, such as [finding parameter sets that maximise the magnitude of some system behaviour](@ref ref). More details on how to use these packages can be found in their [respective](https://docs.sciml.ai/Optimization/stable/) [documentations](https://docs.sciml.ai/DiffEqParamEstim/stable/). ## [Basic example](@id optimization_parameter_fitting_basics) -Let us consider a simple catalysis network, where an enzyme ($E$) turns a substrate ($S$) into a product ($P$): +Let us consider a [Michaelis-Menten enzyme kinetics model](@ref ref), where an enzyme ($E$) converts a substrate ($S$) into a product ($P$): ```@example diffeq_param_estim_1 using Catalyst rn = @reaction_network begin @@ -18,17 +18,17 @@ From some known initial condition, and a true parameter set (which we later want ```@example diffeq_param_estim_1 # Define initial conditions and parameters. u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0] -p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] +ps_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5] # Generate synthetic data. using OrdinaryDiffEq -oprob_true = ODEProblem(rn, u0, (0.0, 10.0), p_true) -true_sol = solve(oprob_true, Tsit5()) -data_sol = solve(oprob_true, Tsit5(); saveat=1.0) +oprob_true = ODEProblem(rn, u0, (0.0, 10.0), ps_true) +true_sol = solve(oprob_true) +data_sol = solve(oprob_true; saveat=1.0) data_ts = data_sol.t[2:end] data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end] -# Plots the true solutions and the (synthetic data) measurements. +# Plots the true solutions and the (synthetic) data measurements. using Plots plot(true_sol; idxs=:P, label="True solution", lw=8) plot!(data_ts, data_vals; label="Measurements", seriestype=:scatter, ms=6, color=:blue) @@ -37,22 +37,22 @@ plot!(data_ts, data_vals; label="Measurements", seriestype=:scatter, ms=6, color Next, we will use DiffEqParamEstim to build a loss function to measure how well our model's solutions fit the data. ```@example diffeq_param_estim_1 using DiffEqParamEstim, Optimization -p_dummy = [:kB => 0.0, :kD => 0.0, :kP => 0.0] -oprob = ODEProblem(rn, u0, (0.0, 10.0), p_dummy) +ps_dummy = [:kB => 0.0, :kD => 0.0, :kP => 0.0] +oprob = ODEProblem(rn, u0, (0.0, 10.0), ps_dummy) loss_function = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, data_vals), Optimization.AutoForwardDiff(); - maxiters=10000, verbose=false, save_idxs=4) + maxiters = 10000, verbose = false, save_idxs = 4) nothing # hide ``` To `build_loss_objective` we provide the following arguments: - `oprob`: The `ODEProblem` with which we simulate our model (using some dummy parameter values, since we do not know these). -- `Tsit5()`: The numeric integrator we wish to simulate our model with. +- `Tsit5()`: The [numeric integrator](@ref ref) we wish to simulate our model with. - `L2Loss(data_ts, data_vals)`: Defines the loss function. While [other alternatives](https://docs.sciml.ai/DiffEqParamEstim/stable/getting_started/#Alternative-Cost-Functions-for-Increased-Robustness) are available, `L2Loss` is the simplest one (measuring the sum of squared distances between model simulations and data measurements). Its first argument is the time points at which the data is collected, and the second is the data's values. - `Optimization.AutoForwardDiff()`: Our choice of [automatic differentiation](https://en.wikipedia.org/wiki/Automatic_differentiation) framework. Furthermore, we can pass any number of additional optional arguments, these are then passed to the internal `solve()` function (which is used to solve our ODE). Here we provide the following additional arguments: -- `maxiters=10000`: If the ODE integrator takes a very large number of steps, that can be a sign of a very poor fit (or stiffness in the ODEs, but that is not a concern for our current example). Reducing the `maxiters` threshold reduces the time we waste on evaluating such models. -- `verbose=false`: The simulation of models with highly unsuitable parameter sets typically generate various warnings (such as premature simulation termination due to reaching `maxiters` timesteps). To avoid an overflow of such (here unnecessary) warnings, as we evaluate a large number of parameter sets, we turn warnings off. -- `save_idxs=4`: The measured species ($P$) is the 4th species in our species vector (`species(rn)`). Since we only assume data is available for $P(t)$ there is no reason to save any other species. +- `maxiters = 10000`: If the ODE integrator takes a very large number of steps, that can be a sign of a very poor fit (or stiffness in the ODEs, but that is not a concern for our current example). Reducing the `maxiters` threshold reduces the time we waste on evaluating such models. +- `verbose = false`: The simulation of models with highly unsuitable parameter sets typically generate various warnings (such as premature simulation termination due to reaching `maxiters` time steps). To avoid an overflow of such (here unnecessary) warnings, as we evaluate a large number of parameter sets, we turn warnings off. +- `save_idxs = 4`: The measured species ($P$) is the 4th species in our species vector (`species(rn)`). Since data is available for $P(t)$, we will only save the value of this species. Now we can create an `OptimizationProblem` using our `loss_function` and some initial guess of parameter values from which the optimiser will start: ```@example diffeq_param_estim_1 @@ -63,7 +63,7 @@ nothing # hide !!! note `OptimizationProblem` cannot currently accept parameter values in the form of a map (e.g. `[:kB => 1.0, :kD => 1.0, :kP => 1.0]`). These must be provided as individual values (using the same order as the parameters occur in in the `parameters(rs)` vector). Similarly, `build_loss_objective`'s `save_idxs` uses the species' indexes, rather than the species directly. These inconsistencies should be remedied in future DiffEqParamEstim releases. -Finally, we can optimise `optprob` to find the parameter set that best fits our data. Optimization.jl only provides a few optimisation methods natively. However, for each supported optimisation package, it provides a corresponding wrapper-package to import that optimisation package for use with Optimization.jl. E.g., if we wish to use [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl)'s [Nelder-Mead](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method) method, we must install and import the OptimizationOptimJL package. A summary of all, by Optimization.jl, supported optimisation packages can be found [here](https://docs.sciml.ai/Optimization/stable/#Overview-of-the-Optimizers). Here, we import the Optim.jl package and uses it to minimise our cost function (thus finding a parameter set that fits the data): +Finally, we can optimise `optprob` to find the parameter set that best fits our data. Optimization.jl only provides a few optimisation methods natively. However, for each supported optimisation package, it provides a corresponding wrapper-package to import that optimisation package for use with Optimization.jl. E.g., if we wish to use [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl)'s [Nelder-Mead](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method) method, we must install and import the OptimizationOptimJL package. A summary of all, by Optimization.jl supported, optimisation packages can be found [here](https://docs.sciml.ai/Optimization/stable/#Overview-of-the-Optimizers). Here, we import the Optim.jl package and uses it to minimise our cost function (thus finding a parameter set that fits the data): ```@example diffeq_param_estim_1 using OptimizationOptimJL optsol = solve(optprob, Optim.NelderMead()) @@ -87,7 +87,7 @@ sol = solve(optprob, NLopt.LD_LBFGS()) nothing # hide ``` -## Optimisation problems with data for multiple species +## [Optimisation problems with data for multiple species](@id optimization_parameter_fitting_multiple_species) Imagine that, in our previous example, we had measurements of the concentration of both *S* and *P*: ```@example diffeq_param_estim_1 data_vals_S = (0.8 .+ 0.4*rand(10)) .* data_sol[:S][2:end] @@ -98,32 +98,32 @@ plot!(data_ts, data_vals_S; label="Measured S", seriestype=:scatter, ms=6, color plot!(data_ts, data_vals_P; label="Measured P", seriestype=:scatter, ms=6, color=:red) ``` -In this case we would have to use the `L2Loss(data_ts, hcat(data_vals_S, data_vals_P))` and `save_idxs=[1,4]` arguments in `loss_function`: +In this case we would have to use the `L2Loss(data_ts, hcat(data_vals_S, data_vals_P))` and `save_idxs=[1, 4]` arguments in `loss_function`: ```@example diffeq_param_estim_1 -loss_function_S_P = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, Array(hcat(data_vals_S, data_vals_P)')), Optimization.AutoForwardDiff(); maxiters=10000, verbose=false, save_idxs=[1,4]) +loss_function_S_P = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, Array(hcat(data_vals_S, data_vals_P)')), Optimization.AutoForwardDiff(); maxiters=10000, verbose=false, save_idxs=[1, 4]) nothing # hide ``` Here, `Array(hcat(data_vals_S, data_vals_P)')` is required to put the data in the right form (in this case, a 2x10 matrix). We can now fit our model to data and plot the results: ```@example diffeq_param_estim_1 -optprob_S_P = OptimizationProblem(loss_function_S_P, [1.0,1.0, 1.0]) +optprob_S_P = OptimizationProblem(loss_function_S_P, [1.0, 1.0, 1.0]) optsol_S_P = solve(optprob_S_P, Optim.NelderMead()) -oprob_fitted_S_P = remake(oprob; p=optsol_S_P.u) -fitted_sol_S_P = solve(oprob_fitted_S_P, Tsit5()) -plot!(fitted_sol_S_P; idxs=[:S, :P], label="Fitted solution", linestyle=:dash, lw=6, color=[:lightblue :pink]) +oprob_fitted_S_P = remake(oprob; p = optsol_S_P.u) +fitted_sol_S_P = solve(oprob_fitted_S_P) +plot!(fitted_sol_S_P; idxs=[:S, :P], label="Fitted solution", linestyle = :dash, lw = 6, color = [:lightblue :pink]) ``` -## Setting parameter constraints and boundaries -Sometimes, it is desired to set boundaries on parameter values. Indeed, this can speed up the optimisation process (by preventing searching through unfeasible parts of parameter space), and can also be a requirement for some optimisation methods. This can be done by passing the `lb` (lower bounds) and `up` (upper bounds) arguments to `OptimizationProblem`. These are vectors (of the same length as the number of parameters), with each argument corresponding to the boundary value of the parameter with the same index (as used in the `parameters(rn)` vector). If we wish to constrain each parameter to the interval $(0.1, 10.0)$ this can be done through: +## [Setting parameter constraints and boundaries](@id optimization_parameter_fitting_constraints) +Sometimes, it is desirable to set boundaries on parameter values. Indeed, this can speed up the optimisation process (by preventing searching through unfeasible parts of parameter space), and can also be a requirement for some optimisation methods. This can be done by passing the `lb` (lower bounds) and `up` (upper bounds) arguments to `OptimizationProblem`. These are vectors (of the same length as the number of parameters), with each argument corresponding to the boundary value of the parameter with the same index (as used in the `parameters(rn)` vector). If we wish to constrain each parameter to the interval $(0.1, 10.0)$ this can be done through: ```@example diffeq_param_estim_1 -optprob = OptimizationProblem(loss_function, [1.0, 1.0, 1.0]; lb = [0.1, 0.1, 0.1], ub = [10.0, 10.0, 10.0]) +optprob = OptimizationProblem(loss_function, [1.0, 1.0, 1.0]; lb = [1e-1, 1e-1, 1e-1], ub = [1e1, 1e1, 1e1]) nothing # hide ``` In addition to boundaries, Optimization.jl also supports setting [linear and non-linear constraints](https://docs.sciml.ai/Optimization/stable/tutorials/constraints/#constraints) on its output solution for some optimizers. -## Parameter fitting with known parameters +## [Parameter fitting with known parameters](@id optimization_parameter_fitting_known_parameters) If we from previous knowledge know that $kD = 0.1$, and only want to fit the values of $kB$ and $kP$, this can be achieved through `build_loss_objective`'s `prob_generator` argument. First, we create a function (`fixed_p_prob_generator`) that modifies our `ODEProblem` to incorporate this knowledge: ```@example diffeq_param_estim_1 fixed_p_prob_generator(prob, p) = remake(prob; p = vcat(p[1], 0.1, p[2])) @@ -142,8 +142,8 @@ optsol_fixed_kD = solve(optprob_fixed_kD, Optim.NelderMead()) nothing # hide ``` -## [Fitting parameters on the logarithmic scale](@id optimization_parameter_fitting_logarithmic_scale) -Often it can be advantageous to fit parameters on a [logarithmic, rather than linear, scale](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008646). The best way to proceed is to simply replace each parameter in the model definition by its logarithmic version: +## [Fitting parameters on the logarithmic scale](@id optimization_parameter_fitting_log_scale) +Often it can be advantageous to fit parameters on a [logarithmic scale (rather than on a linear scale)](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008646). The best way to do this is to simply replace each parameter in the model definition by its logarithmic version: ```@example diffeq_param_estim_2 using Catalyst rn = @reaction_network begin @@ -152,23 +152,39 @@ rn = @reaction_network begin 10^kP, SE --> P + E end ``` -And then proceeding, by keeping in mind that parameter values are logarithmic. Here, setting +And then going forward, by keeping in mind that parameter values are logarithmic. Here, setting ```@example diffeq_param_estim_2 p_true = [:kB => 0.0, :kD => -1.0, :kP => 10^(0.5)] nothing # hide ``` corresponds to the same true parameter values as used previously (`[:kB => 1.0, :kD => 0.1, :kP => 0.5]`). -## Parameter fitting to multiple experiments +## [Parameter fitting to multiple experiments](@id optimization_parameter_fitting_multiple_experiments) Say that we had measured our model for several different initial conditions, and would like to fit our model to all these measurements simultaneously. This can be done by first creating a [corresponding `EnsembleProblem`](@ref advanced_simulations_ensemble_problems). How to then create loss functions for these are described in more detail [here](https://docs.sciml.ai/DiffEqParamEstim/stable/tutorials/ensemble/). -## Optimisation solver options +## [Optimisation solver options](@id optimization_parameter_fitting_solver_options) Optimization.jl supports various [optimisation solver options](https://docs.sciml.ai/Optimization/stable/API/solve/) that can be supplied to the `solve` command. For example, to set a maximum number of seconds (after which the optimisation process is terminated), you can use the `maxtime` argument: ```@example diffeq_param_estim_1 -optsol_fixed_kD = solve(optprob, Optim.NelderMead(); maxtime=100) +optsol_fixed_kD = solve(optprob, Optim.NelderMead(); maxtime = 100) nothing # hide ``` +--- +## [Citation](@id structural_identifiability_citation) +If you use this functionality in your research, please cite the following paper to support the authors of the Optimization.jl package: +``` +@software{vaibhav_kumar_dixit_2023_7738525, + author = {Vaibhav Kumar Dixit and Christopher Rackauckas}, + month = mar, + publisher = {Zenodo}, + title = {Optimization.jl: A Unified Optimization Package}, + version = {v3.12.1}, + doi = {10.5281/zenodo.7738525}, + url = {https://doi.org/10.5281/zenodo.7738525}, + year = 2023 +} +``` + --- ## References [^1]: [Alejandro F. Villaverde, Dilan Pathirana, Fabian Fröhlich, Jan Hasenauer, Julio R. Banga, *A protocol for dynamic model calibration*, Briefings in Bioinformatics (2023).](https://academic.oup.com/bib/article/23/1/bbab387/6383562?login=false) \ No newline at end of file diff --git a/test/reactionsystem_core/events.jl b/test/reactionsystem_core/events.jl index 5a5c4b509f..ab15dd68ae 100644 --- a/test/reactionsystem_core/events.jl +++ b/test/reactionsystem_core/events.jl @@ -243,10 +243,6 @@ let sol_dsl = solve(ODEProblem(rn_dsl, u0, tspan, ps), Tsit5()) sol_prog = solve(ODEProblem(rn_prog, u0, tspan, ps), Tsit5()) @test sol_dsl == sol_prog - - sol_dsl = solve(SDEProblem(rn_dsl, u0, tspan, ps), ImplicitEM(); seed = 1234) - sol_prog = solve(SDEProblem(rn_prog, u0, tspan, ps), ImplicitEM(); seed = 1234) - @test sol_dsl == sol_prog end # Checks that misformatted events yields errors in the DSL. @@ -313,6 +309,65 @@ end ### Additional Correctness Tests ### +# Tests that events are properly triggered for SDEs. +# Tests for continuous events, and all three types of discrete events. +let + # Creates model with all types of events. The `e` parameters track whether events are triggered. + rn = @reaction_network begin + @parameters e1=0 e2=0 e3=0 e4=0 + @continuous_events begin + [X ~ 1000.0] => [e1 ~ 1] + end + @discrete_events begin + [1.0] => [e2 ~ 1] + 1.0 => [e3 ~ 1] + (Y > 1000.0) & (e4==0) => [e4 ~ 1] + end + (p,d), 0 <--> X + (p,d), 0 <--> Y + end + + # Simulates the model for conditions where it *definitely* will cross `X = 1000.0` + u0 = [:X => 999.9, :Y => 999.9] + ps = [:p => 10.0, :d => 0.001] + sprob = SDEProblem(rn, u0, (0.0, 2.0), ps) + sol = solve(sprob, ImplicitEM(); seed) + + # Checks that all `e` parameters have been updated properly. + @test sol.ps[:e1] == 1 + @test sol.ps[:e2] == 1 + @test sol.ps[:e3] == 1 + @test sol.ps[:e4] == 1 +end + +# Tests that events are properly triggered for Jump simulations. +# Tests for all three types of discrete events. +let + # Creates model with all types of events. The `e` parameters track whether events are triggered. + rn = @reaction_network begin + @parameters e1=0 e2=0 e3=0 + @discrete_events begin + [1.0] => [e1 ~ 1] + # 1.0 => [e2 ~ 1] + (X > 1000.0) & (e3==0) => [e3 ~ 1] + end + (p,d), 0 <--> X + end + + # Simulates the model for conditions where it *definitely* will cross `X = 1000.0` + u0 = [:X => 999] + ps = [:p => 10.0, :d => 0.001] + dprob = DiscreteProblem(rn, u0, (0.0, 2.0), ps) + jprob = JumpProblem(rn, dprob, Direct(); rng) + sol = solve(jprob, SSAStepper(); seed) + + # Checks that all `e` parameters have been updated properly. + # Note that periodic discrete events are currently broken for jump processes. + @test sol.ps[:e1] == 1 + @test_broken sol.ps[:e2] == 1 + @test sol.ps[:e3] == 1 +end + # Compares simulations using MTK type events with those generated through callbacks. # Jump simulations must be handles differently (since these only accepts discrete callbacks). # Checks for all types of discrete callbacks, and for continuous callbacks. diff --git a/test/reactionsystem_core/higher_order_reactions.jl b/test/reactionsystem_core/higher_order_reactions.jl index f575cd5f55..7922905792 100644 --- a/test/reactionsystem_core/higher_order_reactions.jl +++ b/test/reactionsystem_core/higher_order_reactions.jl @@ -103,7 +103,7 @@ let dprob_alt2 = DiscreteProblem(u0_alt2, (0.0, 100.0), ps_alt2) jprob_alt2 = JumpProblem(dprob_alt2, Direct(), higher_order_network_alt2...; rng = StableRNG(1234)) - # Simualtes the models. + # Simulates the models. sol_base = solve(jprob_base, SSAStepper(); seed, saveat = 1.0) sol_alt1 = solve(jprob_alt1, SSAStepper(); seed, saveat = 1.0) sol_alt2 = solve(jprob_alt2, SSAStepper(); seed, saveat = 1.0) diff --git a/test/reactionsystem_core/symbolic_stoichiometry.jl b/test/reactionsystem_core/symbolic_stoichiometry.jl index 745a9c50af..797a78f239 100644 --- a/test/reactionsystem_core/symbolic_stoichiometry.jl +++ b/test/reactionsystem_core/symbolic_stoichiometry.jl @@ -192,64 +192,68 @@ end # Tests symbolic stoichiometries in simulations. # Tests for decimal numbered symbolic stoichiometries. let + # Declares models. The references models have the `n` parameters so they can use the same + # parameter vectors as the non-reference ones. rs_int = @reaction_network begin - @parameters n1::Int64 n2::Int64 - p, 0 -->X - k, n1*X --> n2*Y - d, Y --> 0 + @parameters n::Int64 + (k1, k2), n*X1 <--> X2 end rs_dec = @reaction_network begin - @parameters n1::Float64 n2::Float64 - p, 0 -->X - k, n1*X --> n2*Y - d, Y --> 0 + @parameters n::Float64 + (k1, k2), n*X1 <--> X2 end rs_ref_int = @reaction_network begin - @parameters n1::Int64 n2::Int64 - p, 0 -->X - k, 3*X --> 4*Y - d, Y --> 0 + @parameters n::Int64 + (k1, k2), 3*X1 <--> X2 end rs_ref_dec = @reaction_network begin - @parameters n1::Float64 n2::Float64 - p, 0 -->X - k, 0.5*X --> 2.5*Y - d, Y --> 0 + @parameters n::Float64 + (k1, k2), 2.5*X1 <--> X2 end - u0 = [:X => 8, :Y => 8] - ps_int = (:p => 2.0, :k => 0.01, :n1 => 3, :n2 => 4, :d => 0.2) - ps_dec = [:p => 2.0, :k => 0.01, :n1 => 0.5, :n2 => 2.5, :d => 0.2] - tspan = (0.0, 10.0) + # Set simulation settings. Initial conditions are design to start, more or less, at + # steady state concentrations. + # Values are selected so that stochastic tests should always pass within the bounds (independent + # of seed). + u0_int = [:X1 => 150, :X2 => 600] + u0_dec = [:X1 => 100, :X2 => 600] + tspan_det = (0.0, 1.0) + tspan_stoch = (0.0, 10000.0) + ps_int = (:k1 => 0.00001, :k2 => 0.01, :n => 3) + ps_dec = (:k1 => 0.00001, :k2 => 0.01, :n => 2.5) # Test ODE simulations with integer coefficients. - oprob_int = ODEProblem(rs_int, u0, tspan, ps_int) - oprob_int_ref = ODEProblem(rs_ref_int, u0, tspan, ps_int) - @test solve(oprob_int, Tsit5()) ≈ solve(oprob_int_ref, Tsit5()) + oprob_int = ODEProblem(rs_int, u0_int, tspan_det, ps_int) + oprob_int_ref = ODEProblem(rs_ref_int, u0_int, tspan_det, ps_int) + @test solve(oprob_int, Tsit5())[:X1][end] ≈ solve(oprob_int_ref, Tsit5())[:X1][end] # Test ODE simulations with decimal coefficients. - oprob_dec = ODEProblem(rs_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false) - oprob_dec_ref = ODEProblem(rs_ref_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false) - @test solve(oprob_dec, Tsit5()) ≈ solve(oprob_dec_ref, Tsit5()) + oprob_dec = ODEProblem(rs_dec, u0_dec, tspan_det, ps_dec; combinatoric_ratelaws = false) + oprob_dec_ref = ODEProblem(rs_ref_dec, u0_dec, tspan_det, ps_dec; combinatoric_ratelaws = false) + @test solve(oprob_dec, Tsit5())[:X1][end] ≈ solve(oprob_dec_ref, Tsit5())[:X1][end] # Test SDE simulations with integer coefficients. - sprob_int = SDEProblem(rs_int, u0, tspan, ps_int) - sprob_int_ref = SDEProblem(rs_ref_int, u0, tspan, ps_int) - @test solve(sprob_int, ImplicitEM(); seed) ≈ solve(sprob_int_ref, ImplicitEM(); seed) + sprob_int = SDEProblem(rs_int, u0_int, tspan_stoch, ps_int) + sprob_int_ref = SDEProblem(rs_ref_int, u0_dec, tspan_stoch, ps_int) + ssol_int = solve(sprob_int, ImplicitEM(); seed) + ssol_int_ref = solve(sprob_int_ref, ImplicitEM(); seed) + @test mean(ssol_int[:X1]) ≈ mean(ssol_int_ref[:X1]) atol = 2*1e0 # Test SDE simulations with decimal coefficients. - sprob_dec = SDEProblem(rs_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false) - sprob_dec_ref = SDEProblem(rs_ref_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false) - @test solve(sprob_dec, ImplicitEM(); seed) ≈ solve(sprob_dec_ref, ImplicitEM(); seed) - - # Tests jump simulations with integer coefficients. - dprob_int = DiscreteProblem(rs_int, u0, (0.0, 100000.0), ps_int) - dprob_int_ref = DiscreteProblem(rs_ref_int, u0, (0.0, 100000.0), ps_int) - jprob_int = JumpProblem(rs_int, dprob_int, Direct(); rng) - jprob_int_ref = JumpProblem(rs_ref_int, dprob_int_ref, Direct(); rng) - sol_int = solve(jprob_int, SSAStepper(); seed) - sol_int_ref = solve(jprob_int_ref, SSAStepper(); seed) - @test mean(sol_int[:Y]) ≈ mean(sol_int_ref[:Y]) atol = 1e-2 rtol = 1e-2 + sprob_dec = SDEProblem(rs_dec, u0_dec, tspan_stoch, ps_dec; combinatoric_ratelaws = false) + sprob_dec_ref = SDEProblem(rs_ref_dec, u0_dec, tspan_stoch, ps_dec; combinatoric_ratelaws = false) + ssol_dec = solve(sprob_dec, ImplicitEM(); seed) + ssol_dec_ref = solve(sprob_dec_ref, ImplicitEM(); seed) + @test mean(ssol_dec[:X1]) ≈ mean(ssol_dec_ref[:X1]) atol = 2*1e0 + + # Test Jump simulations with integer coefficients. + dprob_int = DiscreteProblem(rs_int, u0_int, tspan_stoch, ps_int) + dprob_int_ref = DiscreteProblem(rs_ref_int, u0_int, tspan_stoch, ps_int) + jprob_int = JumpProblem(rs_int, dprob_int, Direct(); rng, save_positions = (false, false)) + jprob_int_ref = JumpProblem(rs_ref_int, dprob_int_ref, Direct(); rng, save_positions = (false, false)) + jsol_int = solve(jprob_int, SSAStepper(); seed, saveat = 1.0) + jsol_int_ref = solve(jprob_int_ref, SSAStepper(); seed, saveat = 1.0) + @test mean(jsol_int[:X1]) ≈ mean(jsol_int_ref[:X1]) atol = 1e-2 rtol = 1e-2 end # Check that jump simulations (implemented with and without symbolic stoichiometries) yield simulations diff --git a/test/simulation_and_solving/simulate_SDEs.jl b/test/simulation_and_solving/simulate_SDEs.jl index 4b1c935435..5684db4132 100644 --- a/test/simulation_and_solving/simulate_SDEs.jl +++ b/test/simulation_and_solving/simulate_SDEs.jl @@ -128,17 +128,6 @@ let duW = zeros(size(rn_manual.nrp)) rn_manual.g(duW, u0_2, ps_2, 0.0) @test duW ≈ g_eval(rn_catalyst, u0_1, ps_1, 0.0) - - # Compares simulation with identical seed. - # Cannot test the third network (requires very fiddly parameters for CLE to be defined). - if nameof(rn_catalyst) != :rnc9 - sprob_1 = SDEProblem(rn_catalyst, u0_1, (0.0, 1.0), ps_1) - sprob_2 = SDEProblem(rn_manual.f, rn_manual.g, u0_2, (0.0, 1.0), ps_2; noise_rate_prototype = rn_manual.nrp) - seed = seed = rand(rng, 1:100) - sol1 = solve(sprob_1, ImplicitEM(); seed) - sol2 = solve(sprob_2, ImplicitEM(); seed) - @test sol1[u0_sym] ≈ sol2.u - end end end end diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl index 6425758b89..9c8c02db8a 100644 --- a/test/simulation_and_solving/simulate_jumps.jl +++ b/test/simulation_and_solving/simulate_jumps.jl @@ -6,6 +6,7 @@ using Catalyst, JumpProcesses, Statistics, Test # Sets stable rng number. using StableRNGs rng = StableRNG(12345) +seed = rand(rng, 1:100) # Fetch test functions and networks. include("../test_functions.jl") @@ -16,12 +17,17 @@ include("../test_networks.jl") # Compares jump simulations generated through hCatalyst, and manually created systems. let - # Manually declares jumps to compare Catalyst-generated jump simulations to. + # Declares vectors to store information about each network in. Initial conditions are approximately + # at the steady state. catalyst_networks = [] manual_networks = [] u0_syms = [] ps_syms = [] + u0s = [] + ps = [] + sps = [] + # Manual declaration of the first network (and its simulation settings). rate_1_1(u, p, t) = p[1] rate_1_2(u, p, t) = p[2] * u[1] rate_1_3(u, p, t) = p[3] * u[2] @@ -47,12 +53,16 @@ let jump_1_7 = ConstantRateJump(rate_1_7, affect_1_7!) jump_1_8 = ConstantRateJump(rate_1_8, affect_1_8!) jumps_1 = (jump_1_1, jump_1_2, jump_1_3, jump_1_4, jump_1_5, jump_1_6, jump_1_7, - jump_1_8) + jump_1_8) push!(catalyst_networks, reaction_networks_standard[5]) push!(manual_networks, jumps_1) push!(u0_syms, [:X1, :X2, :X3, :X4]) push!(ps_syms, [:p, :k1, :k2, :k3, :k4, :k5, :k6, :d]) + push!(u0s, [:X1 => 40, :X2 => 30, :X3 => 20, :X4 => 10]) + push!(ps, [:p => 10.0, :k1 => 1.0, :k2 => 1.0, :k3 => 1.0, :k4 => 1.0, :k5 => 1.0, :k6 => 1.0, :d => 1.0]) + push!(sps, :X4) + # Manual declaration of the second network (and its simulation settings). rate_2_1(u, p, t) = p[1] / 10 + p[1] * (u[1]^p[3]) / (u[1]^p[3] + p[2]^p[3]) rate_2_2(u, p, t) = p[4] * u[1] * u[2] rate_2_3(u, p, t) = p[5] * u[3] @@ -83,7 +93,11 @@ let push!(manual_networks, jumps_2) push!(u0_syms, [:X1, :X2, :X3]) push!(ps_syms, [:v, :K, :n, :k1, :k2, :k3, :d]) + push!(u0s, [:X1 => 10, :X2 => 10, :X3 => 10]) + push!(ps, [:v => 20.0, :K => 2.0, :n => 2, :k1 => 0.1, :k2 => 0.1, :k3 => 0.1, :d => 1.0]) + push!(sps, :X3) + # Manual declaration of the third network (and its simulation settings). rate_3_1(u, p, t) = p[1] * binomial(u[1], 1) rate_3_2(u, p, t) = p[2] * binomial(u[2], 2) rate_3_3(u, p, t) = p[3] * binomial(u[2], 2) @@ -107,33 +121,29 @@ let push!(manual_networks, jumps_3) push!(u0_syms, [:X1, :X2, :X3, :X4]) push!(ps_syms, [:k1, :k2, :k3, :k4, :k5, :k6]) + push!(u0s, [:X1 => 15, :X2 => 5, :X3 => 5, :X4 => 5]) + push!(ps, [:k1 => 0.1, :k2 => 0.1, :k3 => 0.1, :k4 => 0.1, :k5 => 0.1, :k6 => 0.1]) + push!(sps, :X4) # Loops through all cases, checks that identical simulations are generated with/without Catalyst. - for (rn_catalyst, rn_manual, u0_sym, ps_sym) in zip(catalyst_networks, manual_networks, u0_syms, ps_syms) - for factor in [5, 50] - seed = rand(rng, 1:100) - u0_1 = rnd_u0_Int64(rn_catalyst, rng; n = factor) - ps_1 = rnd_ps(rn_catalyst, rng; factor = factor/100.0) - dprob_1 = DiscreteProblem(rn_catalyst, u0_1, (0.0, 100.0), ps_1) - jprob_1 = JumpProblem(rn_catalyst, dprob_1, Direct(); rng) - sol1 = solve(jprob_1, SSAStepper(); seed, saveat = 1.0) - - u0_2 = map_to_vec(u0_1, u0_sym) - ps_2 = map_to_vec(ps_1, ps_sym) - dprob_2 = DiscreteProblem(u0_2, (0.0, 100.0), ps_2) - jprob_2 = JumpProblem(dprob_2, Direct(), rn_manual...; rng) - sol2 = solve(jprob_2, SSAStepper(); seed, saveat = 1.0) + for (rn_catalyst, rn_manual, u0_sym, ps_sym, u0_1, ps_1, sp) in + zip(catalyst_networks, manual_networks, u0_syms, ps_syms, u0s, ps, sps) - if nameof(rn_catalyst) == :rnh7 - # Have spent a few hours figuring this one out. For certain seeds it actually works, - # but others not. This feels weird, and I didn't get any longer. I tried using non-random - # parameters/initial conditions, and removing the non-hill function reactions. Problem - # still persists. - @test_broken sol1[u0_sym] == sol2.u - else - @test sol1[u0_sym] == sol2.u - end - end + # Simulates the Catalyst-created model. + dprob_1 = DiscreteProblem(rn_catalyst, u0_1, (0.0, 10000.0), ps_1) + jprob_1 = JumpProblem(rn_catalyst, dprob_1, Direct(); rng) + sol1 = solve(jprob_1, SSAStepper(); seed, saveat = 1.0) + + # Simulates the manually written model + u0_2 = map_to_vec(u0_1, u0_sym) + ps_2 = map_to_vec(ps_1, ps_sym) + dprob_2 = DiscreteProblem(u0_2, (0.0, 10000.0), ps_2) + jprob_2 = JumpProblem(dprob_2, Direct(), rn_manual...; rng) + sol2 = solve(jprob_2, SSAStepper(); seed, saveat = 1.0) + + # Checks that the means are similar (the test have been check that it holds across a large + # number of simulates, even without seed). + @test mean(sol1[sp]) ≈ mean(sol2[findfirst(u0_sym .== sp),:]) rtol = 1e-1 end end