Merge remote-tracking branch 'origin/master' into BK

docs/Project.toml

@@ -26,7 +26,12 @@ OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
 OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
+OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -64,7 +69,12 @@ OptimizationBBO = "0.4"
 OptimizationNLopt = "0.3"
 OptimizationOptimJL = "0.4"
 OptimizationOptimisers = "0.3"
-OrdinaryDiffEq = "6.80.1"
+OrdinaryDiffEqBDF = "1"
+OrdinaryDiffEqDefault = "1"
+OrdinaryDiffEqRosenbrock = "1"
+OrdinaryDiffEqSDIRK = "1"
+OrdinaryDiffEqTsit5 = "1"
+OrdinaryDiffEqVerner = "1"
 Plots = "1.40"
 QuasiMonteCarlo = "0.3"
 SciMLBase = "2.46"

docs/src/api.md

@@ -35,7 +35,7 @@ corresponding chemical reaction ODE models, chemical Langevin equation SDE
 models, and stochastic chemical kinetics jump process models.
 
 ```@example ex1
-using Catalyst, OrdinaryDiffEq, StochasticDiffEq, JumpProcesses, Plots
+using Catalyst, OrdinaryDiffEqTsit5, StochasticDiffEq, JumpProcesses, Plots
 t = default_t()
 @parameters β γ
 @species S(t) I(t) R(t)
@@ -393,4 +393,4 @@ Finally, we provide the following helper functions to plot and animate spatial l
 lattice_plot
 lattice_animation
 lattice_kymograph
-```
+```

docs/src/faqs.md

@@ -5,7 +5,7 @@ One can directly use symbolic variables to index into SciML solution objects.
 Moreover, observables can also be evaluated in this way. For example,
 consider the system
 ```@example faq1
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqTsit5, Plots
 rn = @reaction_network ABtoC begin
   (k₊,k₋), A + B <--> C
 end
@@ -132,7 +132,7 @@ When directly constructing a `ReactionSystem`, we can set the symbolic values to
 have the desired default values, and this will automatically be propagated
 through to the equation solvers:
 ```@example faq3
-using Catalyst, Plots, OrdinaryDiffEq
+using Catalyst, Plots, OrdinaryDiffEqTsit5
 t = default_t()
 @parameters β=1e-4 ν=.01
 @species S(t)=999.0 I(t)=1.0 R(t)=0.0
@@ -175,7 +175,7 @@ Julia `Symbol`s corresponding to each variable/parameter to their values, or
 from ModelingToolkit symbolic variables/parameters to their values. Using
 `Symbol`s we have
 ```@example faq4
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqTsit5
 rn = @reaction_network begin
     α, S + I --> 2I
     β, I --> R

docs/src/index.md

@@ -32,7 +32,7 @@ etc).
 - [Steady states](@ref homotopy_continuation) (and their [stabilities](@ref steady_state_stability)) can be computed for model ODE representations.
 
 #### [Features of Catalyst composing with other packages](@id doc_index_features_composed)
-- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) Can be used to numerically solver generated reaction rate equation ODE models.
+- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) Can be used to numerically solve generated reaction rate equation ODE models.
 - [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) can be used to numerically solve generated Chemical Langevin Equation SDE models.
 - [JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl) can be used to numerically sample generated Stochastic Chemical Kinetics Jump Process models.
 - Support for [parallelization of all simulations](@ref ode_simulation_performance_parallelisation), including parallelization of [ODE simulations on GPUs](@ref ode_simulation_performance_parallelisation_GPU) using [DiffEqGPU.jl](https://github.com/SciML/DiffEqGPU.jl).
@@ -92,7 +92,7 @@ Pkg.add("Catalyst")
 
 Many Catalyst features require the installation of additional packages. E.g. for ODE-solving and simulation plotting
 ```julia
-Pkg.add("OrdinaryDiffEq")
+Pkg.add("OrdinaryDiffEqDefault")
 Pkg.add("Plots")
 ```
 is also needed.
@@ -124,7 +124,7 @@ an ordinary differential equation.
 
 ```@example home_simple_example
 # Fetch required packages.
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqDefault, Plots
 
 # Create model.
 model = @reaction_network begin

docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md

@@ -55,15 +55,15 @@ To import a Julia package into a session, you can use the `using PackageName` co
 using Pkg
 Pkg.add("Catalyst")
 ```
-Here, the Julia package manager package (`Pkg`) is by default installed on your computer when Julia is installed, and can be activated directly. Next, we also wish to install the `OrdinaryDiffEq` and `Plots` packages (for numeric simulation of models, and plotting, respectively).
+Here, the Julia package manager package (`Pkg`) is by default installed on your computer when Julia is installed, and can be activated directly. Next, we install an ODE solver from a sub-library of the larger `OrdinaryDiffEq` package, and install the `Plots` package for making graphs. We will import the recommended default solver from the `OrdinaryDiffEqDefault` sub-library. A full list of `OrdinaryDiffEq` solver sublibraries can be found on the sidebar of [this page](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 ```julia
-Pkg.add("OrdinaryDiffEq")
+Pkg.add("OrdinaryDiffEqDefault")
 Pkg.add("Plots")
 ```
 Once a package has been installed through the `Pkg.add` command, this command does not have to be repeated if we restart our Julia session. We can now import all three packages into our current session with:
 ```@example ex2
 using Catalyst
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 using Plots
 ```
 Here, if we restart Julia, these `using` commands *must be rerun*.
@@ -253,4 +253,4 @@ If you are a new Julia user who has used this tutorial, and there was something
 ---
 ## References
 [^1]: [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)
-[^2]: [Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah, *Julia: A Fresh Approach to Numerical Computing*, SIAM Review (2017).](https://epubs.siam.org/doi/abs/10.1137/141000671)
+[^2]: [Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah, *Julia: A Fresh Approach to Numerical Computing*, SIAM Review (2017).](https://epubs.siam.org/doi/abs/10.1137/141000671)

docs/src/introduction_to_catalyst/introduction_to_catalyst.md

@@ -20,7 +20,7 @@ Pkg.activate("catalyst_introduction")
 
 # packages we will use in this tutorial
 Pkg.add("Catalyst")
-Pkg.add("OrdinaryDiffEq")
+Pkg.add("OrdinaryDiffEqTsit5")
 Pkg.add("Plots")
 Pkg.add("Latexify")
 Pkg.add("JumpProcesses")
@@ -29,7 +29,7 @@ Pkg.add("StochasticDiffEq")
 
 We next load the basic packages we'll need for our first example:
 ```@example tut1
-using Catalyst, OrdinaryDiffEq, Plots, Latexify
+using Catalyst, OrdinaryDiffEqTsit5, Plots, Latexify
 ```
 
 Let's start by using the Catalyst [`@reaction_network`](@ref) macro to specify a
@@ -386,4 +386,4 @@ A more detailed summary of the precise mathematical equations Catalyst can gener
 
 ---
 ## References
-1. [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)
+1. [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)

docs/src/inverse_problems/behaviour_optimisation.md

@@ -17,7 +17,7 @@ end
 ```
 To demonstrate this pulsing behaviour we will simulate the system for an example parameter set. We select an initial condition (`u0`) so the system begins in a steady state.
 ```@example behaviour_optimization
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqTsit5, Plots
 example_p = [:pX => 0.1, :pY => 1.0, :pZ => 1.0]
 tspan = (0.0, 50.0)
 example_u0 = [:X => 0.1, :Y => 0.1, :Z => 1.0]
@@ -107,4 +107,4 @@ If you use this functionality in your research, please cite the following paper
 ---
 ## References
 [^1]: [Mykel J. Kochenderfer, Tim A. Wheeler *Algorithms for Optimization*, The MIT Press (2019).](https://algorithmsbook.com/optimization/files/optimization.pdf)
-[^2]: [Lea Goentoro, Oren Shoval, Marc W Kirschner, Uri Alon *The incoherent feedforward loop can provide fold-change detection in gene regulation*, Molecular Cell (2009).](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2896310/)
+[^2]: [Lea Goentoro, Oren Shoval, Marc W Kirschner, Uri Alon *The incoherent feedforward loop can provide fold-change detection in gene regulation*, Molecular Cell (2009).](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2896310/)

docs/src/inverse_problems/examples/ode_fitting_oscillation.md

@@ -4,7 +4,7 @@ In this example we will use [Optimization.jl](https://github.com/SciML/Optimizat
 First, we fetch the required packages.
 ```@example pe_osc_example
 using Catalyst
-using OrdinaryDiffEq
+using OrdinaryDiffEqRosenbrock
 using Optimization
 using OptimizationOptimisers # Required for the ADAM optimizer.
 using SciMLSensitivity # Required for `Optimization.AutoZygote()` automatic differentiation option.

docs/src/inverse_problems/global_sensitivity_analysis.md

@@ -22,7 +22,7 @@ end
 ```
 We will study the peak number of infected cases's ($max(I(t))$) sensitivity to the system's three parameters. We create a function which simulates the system from a given initial condition and measures this property:
 ```@example gsa_1
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 
 u0 = [:S => 999.0, :I => 1.0, :E => 0.0, :R => 0.0]
 p_dummy = [:β => 0.0, :a => 0.0, :γ => 0.0]
@@ -157,4 +157,4 @@ If you use this functionality in your research, [in addition to Catalyst](@ref d
 
 ---
 ## References
-[^1]: [Saltelli, A et al. *Global Sensitivity Analysis. The Primer*, Wiley (2008).](http://www.andreasaltelli.eu/file/repository/A_Saltelli_Marco_Ratto_Terry_Andres_Francesca_Campolongo_Jessica_Cariboni_Debora_Gatelli_Michaela_Saisana_Stefano_Tarantola_Global_Sensitivity_Analysis_The_Primer_Wiley_Interscience_2008_.pdf)
+[^1]: [Saltelli, A et al. *Global Sensitivity Analysis. The Primer*, Wiley (2008).](http://www.andreasaltelli.eu/file/repository/A_Saltelli_Marco_Ratto_Terry_Andres_Francesca_Campolongo_Jessica_Cariboni_Debora_Gatelli_Michaela_Saisana_Stefano_Tarantola_Global_Sensitivity_Analysis_The_Primer_Wiley_Interscience_2008_.pdf)

docs/src/inverse_problems/optimization_ode_param_fitting.md

@@ -21,7 +21,7 @@ u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0]
 ps_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5]
 
 # Generate synthetic data.
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 oprob_true = ODEProblem(rn, u0, (0.0, 10.0), ps_true)
 true_sol = solve(oprob_true)
 data_sol = solve(oprob_true; saveat=1.0)
@@ -39,7 +39,7 @@ Catalyst.PNG(plot(plt; fmt = :png, dpi = 200)) # hide
 
 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
+using DiffEqParamEstim, Optimization, OrdinaryDiffEqTsit5
 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(); 

docs/src/model_creation/conservation_laws.md

@@ -10,7 +10,7 @@ end
 ```
 If we simulate it, we note that while the concentrations of $X₁$ and $X₂$ change throughout the simulation, the total concentration of $X$ ($= X₁ + X₂$) is constant:
 ```@example conservation_laws
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X₁ => 80.0, :X₂ => 20.0]
 ps = [:k₁ => 10.0, :k₂ => 2.0]
 oprob = ODEProblem(rs, u0, (0.0, 1.0), ps)
@@ -47,7 +47,7 @@ Here, Catalyst encodes all conserved quantities in a single, [vector-valued](@re
 
 Practically, the `remove_conserved = true` argument can be provided when a `ReactionSystem` is converted to an `ODEProblem`:
 ```@example conservation_laws
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X₁ => 80.0, :X₂ => 20.0]
 ps = [:k₁ => 10.0, :k₂ => 2.0]
 oprob = ODEProblem(rs, u0, (0.0, 1.0), ps; remove_conserved = true)
@@ -88,4 +88,4 @@ Finally, the `conservationlaws` function yields a $m \times n$ matrix, where $n$
 ```@example conservation_laws
 conservationlaws(rs)
 ```
-I.e. in this case we have a single conserved quantity, which contains a single copy each of the system's two species.
+I.e. in this case we have a single conserved quantity, which contains a single copy each of the system's two species.

docs/src/model_creation/constraint_equations.md

@@ -28,7 +28,7 @@ creating these two systems.
 Here, to create differentials with respect to time (for our differential equations), we must import the time differential operator from Catalyst. We do this through `D = default_time_deriv()`. Here, `D(V)` denotes the differential of the variable `V` with respect to time.
 
 ```@example ceq1
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqTsit5, Plots
 t = default_t()
 D = default_time_deriv()
 
@@ -73,7 +73,7 @@ plot(sol)
 As an alternative to the previous approach, we could have constructed our
 `ReactionSystem` all at once by directly using the symbolic interface:
 ```@example ceq2
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqTsit5, Plots
 t = default_t()
 D = default_time_deriv()
 
@@ -106,7 +106,7 @@ callback interface is illustrated [here](https://docs.sciml.ai/DiffEqDocs/stable
 
 Let's first create our equations and unknowns/species again
 ```@example ceq3
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqTsit5, Plots
 t = default_t()
 D = default_time_deriv()
 
@@ -172,4 +172,4 @@ plot(sol)
 For a detailed discussion on how to directly use the lower-level but more
 flexible DifferentialEquations.jl event/callback interface, see the
 [tutorial](https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Event-handling-using-callbacks)
-on event handling using callbacks.
+on event handling using callbacks.

docs/src/model_creation/dsl_advanced.md

@@ -100,7 +100,7 @@ end
 ```
 Next, if we simulate the model, we do not need to provide values for species or parameters that have default values. In this case all have default values, so both `u0` and `ps` can be empty vectors:
 ```@example dsl_advanced_defaults
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = []
 tspan = (0.0, 10.0)
 p = []
@@ -142,6 +142,7 @@ end
 ```
 Please note that as the parameter `X₀` does not occur as part of any reactions, Catalyst's DSL cannot infer whether it is a species or a parameter. This must hence be explicitly declared. We can now simulate our model while providing `X`'s value through the `X₀` parameter:
 ```@example dsl_advanced_defaults
+using OrdinaryDiffEqTsit5
 u0 = []
 p = [:X₀ => 1.0, :p => 1.0, :d => 0.5]
 oprob = ODEProblem(rn, u0, tspan, p)
@@ -262,7 +263,7 @@ end
 ```
 Now, we can also declare our initial conditions and parameter values as vectors as well:
 ```@example dsl_advanced_vector_variables
-using OrdinaryDiffEq, Plots # hide
+using OrdinaryDiffEqDefault, Plots # hide
 u0 = [:X => [0.0, 2.0]]
 tspan = (0.0, 1.0)
 ps = [:k => [1.0, 2.0]]
@@ -332,7 +333,7 @@ end
 ```
 We can now simulate our model using normal syntax (initial condition values for observables should not, and can not, be provided):
 ```@example dsl_advanced_observables
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 u0 = [:X => 1.0, :Y => 2.0, :XY => 0.0]
 tspan = (0.0, 10.0)
 ps = [:kB => 1.0, :kD => 1.5]
@@ -489,7 +490,7 @@ X
 ```
 Next, we can now use these to e.g. designate initial conditions and parameter values for model simulations:
 ```@example dsl_advanced_programmatic_unpack
-using OrdinaryDiffEq, Plots # hide
+using OrdinaryDiffEqDefault, Plots # hide
 u0 = [X => 0.1]
 tspan = (0.0, 10.0)
 ps = [p => 1.0, d => 0.2]
@@ -538,4 +539,4 @@ nothing # hide
 ```
 
 !!! note
-    When using interpolation, expressions like `2$spec` won't work; the multiplication symbol must be explicitly included like `2*$spec`.
+    When using interpolation, expressions like `2$spec` won't work; the multiplication symbol must be explicitly included like `2*$spec`.