Merge pull request #845 from SciML/sde_performance_advice

SDE performance tutorial

Author: TorkelE

docs/pages.jl

@@ -27,7 +27,8 @@ pages = Any[
         "model_simulation/simulation_plotting.md",
         "model_simulation/simulation_structure_interfacing.md",
         "model_simulation/ensemble_simulations.md",
-        "model_simulation/ode_simulation_performance.md"
+        "model_simulation/ode_simulation_performance.md",
+        "model_simulation/sde_simulation_performance.md"
     ],
     "Steady state analysis" => Any[
         "steady_state_functionality/homotopy_continuation.md",

docs/src/model_simulation/sde_simulation_performance.md

@@ -0,0 +1,58 @@
+# [Advice for performant SDE simulations](@id sde_simulation_performance)
+While there exist relatively straightforward approaches to manage performance for [ODE](@ref ode_simulation_performance) and jump simulations, this is generally not the case for SDE simulations. Below, we briefly describe some options. However, as one starts to investigate these, one quickly reaches what is (or could be) active areas of research.
+
+## [SDE solver selection](@id sde_simulation_performance_solvers)
+We have previously described how [ODE solver selection](@ref ode_simulation_performance_solvers) can impact simulation performance. Again, it can be worthwhile to investigate solver selection's impact on performance for SDE simulations. Throughout this documentation, we generally use the `STrapezoid` solver as the default choice. However, if the `DifferentialEquations` package is loaded
+```julia
+using DifferentialEquations
+```
+automatic SDE solver selection enabled (just like is the case for ODEs by default). Generally, the automatic SDE solver choice enabled by `DifferentialEquations` is better than just using `STrapezoid`. Next, if performance is critical, it can be worthwhile to check the [list of available SDE solvers](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/) to find one with advantageous performance for a given problem. When doing so, it is important to pick a solver compatible with *non-diagonal noise* and with [*Ito problems*](https://en.wikipedia.org/wiki/It%C3%B4_calculus). 
+
+## [Options for Jacobian computation](@id sde_simulation_performance_jacobian)
+In the section on ODE simulation performance, we describe various [options for computing the system Jacobian](@ref ode_simulation_performance_jacobian), and how these could be used to improve performance for [implicit solvers](@ref ode_simulation_performance_stiffness). These can be used in tandem with implicit SDE solvers (such as `STrapezoid`). However, due to additional considerations during SDE simulations, it is much less certain whether these will actually have any impact on performance. So while these options might be worth reading about and trialling, there is no guarantee that they will be beneficial.
+
+## [Parallelisation on CPUs and GPUs](@id sde_simulation_performance_parallelisation)
+We have previously described how simulation parallelisation can be used to [improve performance when multiple ODE simulations are carried out](@ref ode_simulation_performance_parallelisation). The same approaches can be used for SDE simulations. Indeed, it is often more relevant for SDEs, as these are often re-simulated using identical simulation conditions (to investigate their typical behaviour across many samples). CPU parallelisation of SDE simulations uses the [same approach as ODEs](@ref ode_simulation_performance_parallelisation_CPU). GPU parallelisation requires some additional considerations, which are described below.
+
+### [GPU parallelisation of SDE simulations](@id sde_simulation_performance_parallelisation_GPU)
+GPU parallelisation of SDE simulations uses a similar approach as that for [ODE simulations](@ref ode_simulation_performance_parallelisation_GPU). The main differences are that SDE parallelisation requires a GPU SDE solver (like `GPUEM`) and fixed time stepping.
+
+We will assume that we are using the CUDA GPU hardware, so we will first load the [CUDA.jl](https://github.com/JuliaGPU/CUDA.jl) backend package, as well as DiffEqGPU:
+```julia
+using CUDA, DiffEqGPU
+```
+Which backend package you should use depends on your available hardware, with the alternatives being listed [here](https://docs.sciml.ai/DiffEqGPU/stable/manual/backends/).
+
+Next, we create the `SDEProblem` which we wish to simulate. Like for ODEs, we ensure that all vectors are [static vectors](https://github.com/JuliaArrays/StaticArrays.jl) and that all values are `Float32`s. Here we prepare the parallel simulations of a simple [birth-death process](@ref basic_CRN_library_bd).
+```@example sde_simulation_performance_gpu
+using Catalyst, StochasticDiffEq, StaticArrays
+bd_model = @reaction_network begin
+    (p,d), 0 <--> X
+end
+@unpack X, p, d = bd_model
+
+u0 = @SVector [X => 20.0f0]
+tspan = (0.0f0, 10.0f0)
+ps = @SVector [p => 10.0f0, d => 1.0f0]
+sprob = SDEProblem(bd_model, u0, tspan, ps)
+nothing # hide
+```
+The `SDEProblem` is then used to [create an `EnsembleProblem`](@ref ensemble_simulations_monte_carlo).
+```@example sde_simulation_performance_gpu
+eprob = EnsembleProblem(sprob)
+nothing # hide
+```
+Finally, we can solve our `EnsembleProblem` while:
+- Using a valid GPU SDE solver (either [`GPUEM`](https://docs.sciml.ai/DiffEqGPU/stable/manual/ensemblegpukernel/#DiffEqGPU.GPUEM) or [`GPUSIEA`](https://docs.sciml.ai/DiffEqGPU/stable/manual/ensemblegpukernel/#DiffEqGPU.GPUSIEA)).
+- Designating the GPU ensemble method, `EnsembleGPUKernel` (with the correct GPU backend as input).
+- Designating the number of trajectories we wish to simulate.
+- Designating a fixed time step size.
+
+```julia
+esol = solve(eprob, GPUEM(), EnsembleGPUKernel(CUDA.CUDABackend()); trajectories = 10000, dt = 0.01)
+```
+
+Above we parallelise GPU simulations with identical initial conditions and parameter values. However, [varying these](@ref ensemble_simulations_varying_conditions) is also possible.
+
+### [Multilevel Monte Carlo](@id sde_simulation_performance_parallelisation_mlmc)
+An approach for speeding up parallel stochastic simulations is so-called [*multilevel Monte Carlo approaches*](https://en.wikipedia.org/wiki/Multilevel_Monte_Carlo_method) (MLMC). These are used when a stochastic process is simulated repeatedly using identical simulation conditions. Here, instead of performing all simulations using identical [tolerance](@ref ode_simulation_performance_error), the ensemble is simulated using a range of tolerances (primarily lower ones, which yields faster simulations). Currently, [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) do not have a native implementation for performing MLMC simulations (this will hopefully be added in the future). However, if high performance of parallel SDE simulations is required, these approaches may be worth investigating.

docs/src/model_simulation/simulation_introduction.md

@@ -1,7 +1,7 @@
 # [Model Simulation Introduction](@id simulation_intro)
 Catalyst's core functionality is the creation of *chemical reaction network* (CRN) models that can be simulated using ODE, SDE, and jump simulations. How such simulations are carried out has already been described in [Catalyst's introduction](@ref introduction_to_catalyst). This page provides a deeper introduction, giving some additional background and introducing various simulation-related options. 
 
-Here we will focus on the basics, with other sections of the simulation documentation describing various specialised features, or giving advice on performance. Anyone who plans on using Catalyst's simulation functionality extensively is recommended to also read the documentation on [solution plotting](@ref simulation_plotting), and on how to [interact with simulation problems, integrators, and solutions](@ref simulation_structure_interfacing). Anyone with an application for which performance is critical should consider reading the corresponding page on performance advice for [ODEs](@ref ode_simulation_performance).
+Here we will focus on the basics, with other sections of the simulation documentation describing various specialised features, or giving advice on performance. Anyone who plans on using Catalyst's simulation functionality extensively is recommended to also read the documentation on [solution plotting](@ref simulation_plotting), and on how to [interact with simulation problems, integrators, and solutions](@ref simulation_structure_interfacing). Anyone with an application for which performance is critical should consider reading the corresponding page on performance advice for [ODEs](@ref ode_simulation_performance) or [SDEs](@ref sde_simulation_performance).
 
 ### [Background to CRN simulations](@id simulation_intro_theory)
 This section provides some brief theory on CRN simulations. For details on how to carry out these simulations in actual code, please skip to the following sections.
@@ -169,7 +169,7 @@ nothing # hide
 ```
 
 ## [Performing SDE simulations](@id simulation_intro_SDEs)
-Catalyst uses the [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) package to perform SDE simulations. This section provides a brief introduction, with [StochasticDiffEq's documentation](https://docs.sciml.ai/StochasticDiffEq/stable/) providing a more extensive description. sBy default, Catalyst generates SDEs from CRN models using the chemical Langevin equation.
+Catalyst uses the [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) package to perform SDE simulations. This section provides a brief introduction, with [StochasticDiffEq's documentation](https://docs.sciml.ai/StochasticDiffEq/stable/) providing a more extensive description. By default, Catalyst generates SDEs from CRN models using the chemical Langevin equation. A dedicated section giving advice on how to optimise SDE simulation performance can be found [here](@ref sde_simulation_performance).
 
 SDE simulations are performed in a similar manner to ODE simulations. The only exception is that an `SDEProblem` is created (rather than an `ODEProblem`). Furthermore, the [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) package (rather than the OrdinaryDiffEq package) is required for performing simulations. Here we simulate the two-state model for the same parameter set as previously used:
 ```@example simulation_intro_sde

src/reactionsystem_conversions.jl

@@ -590,7 +590,7 @@ function nonlinear_convert_differentials_check(rs::ReactionSystem)
                     (2) The right-hand side does not contain any differentials.
                 This is generally not permitted.
 
-                If you still would like to perform this conversions, please use the `all_differentials_permitted = true` option. In this case, all differential will be set to `0`. 
+                If you still would like to perform this conversion, please use the `all_differentials_permitted = true` option. In this case, all differentials will be set to `0`. 
                 However, it is recommended to proceed with caution to ensure that the produced nonlinear equation makes sense for your intended application."
             )
         end