Skip to content

Commit

Permalink
small changes
Browse files Browse the repository at this point in the history
  • Loading branch information
TorkelE committed May 27, 2024
1 parent b9731d8 commit 4506fdf
Show file tree
Hide file tree
Showing 3 changed files with 16 additions and 10 deletions.
15 changes: 10 additions & 5 deletions docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md
Original file line number Diff line number Diff line change
Expand Up @@ -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 `DifferentialEquations` 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 also wish to install the `OrdinaryDiffEq` and `Plots` packages (for numeric simulation of models, and plotting, respectively).
```julia
Pkg.add("DifferentialEquations")
Pkg.add("OrdinaryDiffEq")
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 DifferentialEquations
using OrdinaryDiffEq
using Plots
```
Here, if we restart Julia, these `using` commands *must be rerun*.
Expand Down Expand Up @@ -130,7 +130,11 @@ For more information about the numerical simulation package, please see the [Dif
## Additional modelling example
To make this introduction more comprehensive, we here provide another example, using a more complicated model. Instead of simulating our model as concentrations evolve over time, we will now simulate the individual reaction events through the [Gillespie algorithm](https://en.wikipedia.org/wiki/Gillespie_algorithm) (a common approach for adding *noise* to models).

Remember (unless we have restarted Julia) we do not need to activate our packages (through the `using` command) again.
Remember (unless we have restarted Julia) we do not need to activate our packages (through the `using` command) again. However, we do need to install, and then import, the JumpProcesses package (just to perform Gillespie, and other jump, simulations)
```julia
Pkg.add("JumpProcesses")
using JumpProcesses
```

This time, we will declare a so-called [SIR model for an infectious disease](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model). Note that even if this model does not describe a set of chemical reactions, it can be modelled using the same framework. The model consists of 3 species:
* $S$, the amount of *susceptible* individuals.
Expand Down Expand Up @@ -164,12 +168,13 @@ nothing # hide

Previously we have bundled this information into an `ODEProblem` (denoting a deterministic *ordinary differential equation*). Now we wish to simulate our model as a jump process (where each reaction event corresponds to a single jump in the state of the system). We do this by first creating a `DiscreteProblem`, and then using this as an input to a `JumpProblem`.
```@example ex2
using JumpProcesses # hide
dprob = DiscreteProblem(sir_model, u0, tspan, params)
jprob = JumpProblem(sir_model, dprob, Direct())
```
Again, the order in which the inputs are given to the `DiscreteProblem` and the `JumpProblem` is important. The last argument to the `JumpProblem` (`Direct()`) denotes which simulation method we wish to use. For now, we recommend that users simply use the `Direct()` option, and then consider alternative ones (see the [JumpProcesses.jl docs](https://docs.sciml.ai/JumpProcesses/stable/)) when they are more familiar with modelling in Catalyst and Julia.

Finally, we can simulate our model using the `solve` function, and plot the solution using the `plot` function. Here, the `solve` function also has a second argument (`SSAStepper()`). This is a time-stepping algorithm that calls the `Direct` solver to advance a simulation. Again, we recommend at this stage you simply use this option, and then explore exactly what this means at a later stage.
Finally, we can simulate our model using the `solve` function, and plot the solution using the `plot` function. For jump simulations, the `solve` function also requires a second argument (`SSAStepper()`). This is a time-stepping algorithm that calls the `Direct` solver to advance a simulation. Again, we recommend at this stage you simply use this option, and then explore exactly what this means at a later stage.
```@example ex2
sol = solve(jprob, SSAStepper())
sol = solve(jprob, SSAStepper(); seed=1234) # hide
Expand Down
10 changes: 5 additions & 5 deletions docs/src/model_simulation/ode_simulation_performance.md
Original file line number Diff line number Diff line change
Expand Up @@ -235,17 +235,17 @@ plot(0.01:0.01:1.0, map(sol -> sol[:P][end], esol.u), xguide = "kP", yguide = "P
```

Above, we have simply used `EnsembleProblem` as a convenient interface to run a large number of similar simulations. However, these problems have the advantage that they allow the passing of an *ensemble algorithm* to the `solve` command, which describes a strategy for parallelising the simulations. By default, `EnsembleThreads` is used. This parallelises the simulations using [multithreading](https://en.wikipedia.org/wiki/Multithreading_(computer_architecture)) (parallelisation within a single process), which is typically advantageous for small problems on shared memory devices. An alternative is `EnsembleDistributed` which instead parallelises the simulations using [multiprocessing](https://en.wikipedia.org/wiki/Multiprocessing) (parallelisation across multiple processes). To do this, we simply supply this additional solver to the solve command:
```@example ode_simulation_performance_4
```julia
esol = solve(eprob, Tsit5(), EnsembleDistributed(); trajectories=100)
nothing # hide
```
To utilise multiple processes, you must first give Julia access to these. You can check how many processes are available using the `nprocs` (which requires the [Distributed.jl](https://github.com/JuliaLang/Distributed.jl) package):
```@example ode_simulation_performance_4
```julia
using Distributed
nprocs()
```
Next, more processes can be added using `addprocs`. E.g. here we add an additional 4 processes:
```@example ode_simulation_performance_4
```julia
addprocs(4)
nothing # hide
```
Expand All @@ -268,7 +268,7 @@ Furthermore (while not required) to receive good performance, we should also mak
- We should designate all our vectors (i.e. initial conditions and parameter values) as [static vectors](https://github.com/JuliaArrays/StaticArrays.jl).

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:
```@example ode_simulation_performance_5
```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/).
Expand Down Expand Up @@ -306,7 +306,7 @@ We can now simulate our model using a GPU-based ensemble algorithm. Currently, t
- While `EnsembleGPUArray` can use standard ODE solvers, `EnsembleGPUKernel` requires specialised versions (such as `GPUTsit5`). A list of available such solvers can be found [here](https://docs.sciml.ai/DiffEqGPU/dev/manual/ensemblegpukernel/#specialsolvers).

Generally, it is recommended to use `EnsembleGPUArray` for large models (that have at least $100$ variables), and `EnsembleGPUKernel` for smaller ones. Here we simulate our model using both approaches (noting that `EnsembleGPUKernel` requires `GPUTsit5`):
```@example ode_simulation_performance_5
```julia
esol1 = solve(eprob, Tsit5(), EnsembleGPUArray(CUDA.CUDABackend()); trajectories = 10000)
esol2 = solve(eprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()); trajectories = 10000)
nothing # hide
Expand Down
1 change: 1 addition & 0 deletions test/dsl/dsl_basic_model_construction.jl
Original file line number Diff line number Diff line change
Expand Up @@ -438,3 +438,4 @@ let
@test_throws LoadError @eval @reaction k, 0 --> nothing
end


0 comments on commit 4506fdf

Please sign in to comment.