Skip to content

Commit 806a9c2

Browse files
TorkelETorkelE
committed
init
1 parent 9d271bf commit 806a9c2

File tree

1 file changed

+142
-94
lines changed

1 file changed

+142
-94
lines changed

‎docs/src/model_simulation/simulation_structure_interfacing.md

Lines changed: 142 additions & 94 deletions
Original file line numberDiff line numberDiff line change
@@ -1,141 +1,189 @@
11
# [Interfacing problems, integrators, and solutions](@id simulation_structure_interfacing)
2-
When simulating a model, one begins with creating a [problem](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/). Next, a simulation is performed on a problem, during which the state of the simulation is recorded through an [integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/). Finally, the simulation output is returned as a [solution](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/). This tutorial describes how to access, or modify the state, or parameter, values of problems, integrators, and solutions structures.
2+
When simulating a model, one begins with creating a [problem](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/). Next, a simulation is performed on a problem, during which the state of the simulation is recorded through an [integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/). Finally, the simulation output is returned as a [solution](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/). This tutorial describes how to access (or modify) the state (or parameter) values of problem, integrator, and solution structures.
33

4-
Generally, when we have a structure `simulation_struct` and want to interface with the unknown (or parameter) `G`, we use `simulation_struct[:G]` to access the value, and `simulation_struct[:G] = 5.0` to set it to a new value. However, see the following examples for full details.
4+
Generally, when we have a structure `simulation_struct` and want to interface with the unknown (or parameter) `x`, we use `simulation_struct[:x]` to access the value, and `simulation_struct[:x] = 5.0` to set it to a new value. For situations where a value is accessed (or changed) a large number of times, it can improve performance to first create a [specialised getter/setter function](@ref ref).
55

6-
!!! note
7-
The following tutorial will describe how to interface with problems, integrators, and solutions using `[]` notation. An alternative is to use the `ModelingToolkit.getu`, `ModelingToolkit.getp`, `ModelingToolkit.setu`, and `ModelingToolkit.setp` functions. These requires an additional step to use, however, they can also improve performance when a very large number interfaces are carried out.
8-
9-
## Interfacing problem objects
6+
## [Interfacing problem objects](@id simulation_structure_interfacing_problems)
107

11-
We begin by demonstrating how we can interface with problem objects. We will demonstrate using a `ODEProblem`, however, it works similarly for other problem types.
12-
```@example ex1
8+
We begin by demonstrating how we can interface with problem objects. First, we create an `ODEProblem` representation of a [chemical cross-coupling model](@ref ref) (where a catalyst, $C$, couples two substrates, $S₁$ and $S₂$, to form a product, $P$).
9+
```@example structure_indexing
1310
using Catalyst
14-
rn = @reaction_network begin
15-
(k1,k2), X1 <--> X2
11+
cc_system = @reaction_network begin
12+
k₁, S₁ + C --> S₁C
13+
k₂, S₁C + S₂ --> CP
14+
k₃, CP --> C + P
1615
end
1716
18-
u0 = [:X1 => 1.0, :X2 => 5.0]
19-
p = [:k1 => 5.0, :k2 => 2.0]
20-
oprob = ODEProblem(rn, u0, (0.0,10.0), p)
17+
u0 = [:S₁ => 1.0, :C => 0.05, :S₂ => 1.2, :S₁C => 0.0, :CP => 0.0, :P => 0.0]
18+
tspan = (0., 10.0)
19+
ps = [:k₁ => 5.0, :k₂ => 5.0, :k₃ => 100.0]
20+
oprob = ODEProblem(cc_system, u0, tspan, ps)
2121
nothing # hide
2222
```
2323

24-
We can find the value of a state simply by interfacing with the corresponding symbol:
25-
```@example ex1
26-
oprob[:X1]
24+
We can find a species (or [variables](@ref ref)) initial condition value by simply indexing with the species of interest as input. Here we check the initial condition value of $C$:
25+
```@example structure_indexing
26+
oprob[:C]
2727
```
28+
An almost identical notation can be used for parameters, however, here we use `oprob.ps` (rather than `oprob`):
2829
with the notation being identical for parameters:
29-
```@example ex1
30-
oprob.ps[:k1]
30+
```@example structure_indexing
31+
oprob.ps[:k₁]
3132
```
32-
33-
If we want to change an unknown's initial condition value, we use the following notation
34-
```@example ex1
35-
oprob[:X1] = 10.0
33+
To retrieve several species initial condition (or parameter) values, simply give a vector input. Here we check the values of the two substrates ($S₁$ and $S₂$):
34+
```@example structure_indexing
35+
oprob[[:S₁, :S₂]]
3636
```
37-
with parameters using the same notation.
3837

39-
!!! note
40-
When interfacing with a parameter, `.ps` must be appended to the structure uses (e.g. `oprob`). This is not done when species are interfaced with.
41-
42-
#### [Remaking problems using the `remake` function](@id simulation_structure_interfacing_remake)
43-
Typically, when modifying problems, it is recommended to use the `remake` function. Unlike when we do `oprob[:X1] = 10.0` (which modifies the problem in question), `remake` creates a new problem object. The `remake` function takes a problem as input, and any fields you wish to modify (and their new values) as optional inputs. Thus, we can do:
44-
```@example ex1
45-
using DifferentialEquations
46-
@unpack X1, X2, k1, k2 = rn
47-
oprob1 = ODEProblem(rn, u0, (0.0,10.0), p)
48-
oprob2 = remake(oprob1; u0=[X1 => 10.0, X2 => 50.0], tspan=(0.0,100.0), p=[k1 => 50.0,k2 => 20.0])
49-
nothing # hide
50-
```
51-
and we can now check the fields of `oprob2`
52-
```@example ex1
53-
oprob2.u0
38+
We can change a species's initial condition value using a similar notation. Here we increase the initial concentration of $C$ (and also confirm that the new value is stored in an updated `oprob`):
39+
```@example structure_indexing
40+
oprob[:C] = 0.1
41+
oprob[:C]
5442
```
55-
```@example ex1
56-
oprob2.tspan
43+
Again, parameter values can be changed using a similar notation, however, again requiring `oprob.ps` notation:
44+
```@example structure_indexing
45+
oprob.ps[:k₁] = 10.0
46+
oprob.ps[:k₁]
5747
```
58-
```@example ex1
59-
oprob2.p
48+
Finally, vectors can be used to update multiple quantities simultaneously
49+
```@example structure_indexing
50+
oprob[[:S₁, :S₂]] = [0.5, 0.3]
51+
oprob[[:S₁, :S₂]]
6052
```
61-
Please note that, currently, `remake` does not work while giving `Symbol`s as input (e.g `[:X1 => 10.0, :X2 => 50.0]`), but we need to unpack the symbolic variables and use them instead (please see the end of this tutorial for more information on using symbolic variables rather than `Symbol`s).
53+
Generally, when updating problems, it is often better to use the [`remake` function](@ref simulation_structure_interfacing_problems_remake) (especially when several values are updated).
6254

63-
When using `remake`, we only have to provide the fields that we actually wish to change, e.g.
64-
```@example ex1
65-
oprob3 = remake(oprob1; u0=[X1 => 10.0, X2 => 50.0])
66-
nothing # hide
55+
A problem's time span can be accessed through the `tspan` field:
56+
```@example structure_indexing
57+
oprob.tspan
6758
```
68-
will only update the initial conditions.
6959

60+
!!! note
61+
Here we have used an `ODEProblem`to demonstrate all interfacing functionality. However, identical workflows work for the other problem types.
62+
63+
### [Remaking problems using the `remake` function](@id simulation_structure_interfacing_problems_remake)
64+
The `remake` function offers an (to indexing) alternative approach for updating problems. Unlike indexing, `remake` creates a new problem (rather than updating the old one). Furthermore, it permits the updating of several values simultaneously. The `remake` function takes the following inputs:
65+
- The problem that is remade.
66+
- (optionally) `u0`: A vector with initial conditions that should be updated. The vector takes the same form as normal initial condition vectors, but does not need to be complete (in which case only a subset of the initial conditions are updated).
67+
- (optionally) `tspan`: An updated time span (using the same format as time spans normally are given in).
68+
- (optionally) `p`: A vector with parameters that should be updated. The vector takes the same form as normal parameter vectors, but does not need to be complete (in which case only a subset of the parameters are updated).
69+
70+
Here we modify our problem to increase the initial condition concentrations of the two substrates ($S₁$ and $S₂$), and also confirm that the new problem is different from the old (unchanged) one:
71+
```@example structure_indexing
72+
using OrdinaryDiffEq
73+
oprob_new = remake(oprob; u0 = [:S₁ => 5.0, :S₂ => 2.5])
74+
oprob_new == oprob
75+
```
76+
Here, we instead use `remake` to simultaneously update a large number of fields:
77+
```@example structure_indexing
78+
oprob_new_2 = remake(oprob; u0 = [:C => 0.2], tspan = (0.0, 20.0), p = [:k₁ => 2.0, :k₂ => 2.0])
79+
nothing # hide
80+
```
7081

71-
## Interfacing integrator objects
82+
## [Interfacing integrator objects](@id simulation_structure_interfacing_integrators)
7283

73-
During a simulation, the solution is stored in an integrator object, we will here describe how to interface with these. The primary circumstance under which a user may wish to do so is when using [callbacks](@ref advanced_simulations_callbacks). We can create an integrator by calling `init` on our problem ([while circumstances where the user might want to use `init` function exist](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Initialization-and-Stepping), since integrators are automatically created during simulations, these are rare):
74-
```@example ex1
84+
During a simulation, the solution is stored in an integrator object. Here, we will describe how to interface with these. The almost exclusive circumstance when integrator-interfacing is relevant is when simulation events are [implemented through callbacks](@ref ref). However, to demonstrate integrator indexing in this tutorial, we will create one through the `init` function (while circumstances where one might [want to use `init` function exist](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Initialization-and-Stepping), since integrators are automatically created during simulations, these are rare).
85+
```@example structure_indexing
7586
integrator = init(oprob)
87+
nothing # hide
7688
```
77-
Using a similar syntax to problems, we can get the current values of an unknown within the integrator:
78-
```@example ex1
79-
integrator[:X1]
89+
We can interface with our integrator using an identical syntax as [was used for problems](@ref simulation_structure_interfacing_problems) (with the exception that `remake` is not available). Here we update, and then check the values of, first the species $C$ and then the parameter $k₁$:
90+
```@example structure_indexing
91+
integrator[:C] = 0.0
92+
integrator[:C]
8093
```
8194
or a parameter:
82-
```@example ex1
83-
integrator.ps[:k1]
84-
```
85-
Similarly, we can update their values using:
86-
```@example ex1
87-
integrator[:X1] = 10.0
95+
```@example structure_indexing
96+
integrator.ps[:k₂] = 1.0
97+
integrator.ps[:k₂]
8898
```
89-
Please read [this](@ref advanced_simulations_ssa_callbacks) with regards to updating integrators of `JumpProblem`s.
99+
Note that here, species-interfacing yields (or changes) a simulation's current value for a species, not its initial condition.
90100

101+
If you are interfacing with jump simulation integrators, please read [this, highly relevant, section](@ref ref).
91102

92103
## [Interfacing solution objects](@id simulation_structure_interfacing_solutions)
93104

94105
Finally, we consider solution objects. First, we simulate our problem:
95-
```@example ex1
106+
```@example structure_indexing
96107
sol = solve(oprob)
108+
nothing # hide
97109
```
98-
For solutions, when we access an unknown, we get its whole simulation vector:
99-
```@example ex1
100-
sol[:X1]
110+
Next, we can access the simulation's values using the same notation as previously. When we check a species's value, its full across the full simulation is returned as a vector:
111+
```@example structure_indexing
112+
sol[:P]
101113
```
102-
while when we access a parameter we only get a single value:
103-
```@example ex1
104-
sol.ps[:k1]
114+
Parameter values can also be accessed (however, here we only get a single value):
115+
```@example structure_indexing
116+
sol.ps[:k₃]
105117
```
106-
Finally, we note that we cannot change the values of solution unknowns or parameters (i.e. both `sol[:X1] = 0.0` and `sol[:k1] = 0.0` generate errors).
118+
Unlike for problems and integrators, species or parameter values of solutions cannot be changed.
107119

108-
## [Interfacing using symbolic representation](@id simulation_structure_interfacing_symbolic_representation)
120+
A vector with the time values for all simulation time steps can be retrieved using
121+
```@example structure_indexing
122+
sol.t
123+
```
109124

110-
Catalyst is built on an *intermediary representation* implemented by (ModelingToolkit.jl)[https://github.com/SciML/ModelingToolkit.jl]. ModelingToolkit is a modelling framework where one first declares a set of symbolic variables and parameters using e.g.
111-
```@example ex2
112-
using Catalyst # hide
113-
using ModelingToolkit
114-
t = default_t()
115-
@parameters σ ρ β
116-
@variables x(t) y(t) z(t)
117-
nothing # hide
125+
## [Interfacing using specialised getter/setter functions](@id simulation_structure_interfacing_functions)
126+
Internally, species and parameter values are stored in vectors. Whenever e.g. `oprob[:C]` is called, Julia must first find which index in the storage vector $C$ is stored in. Next, its value can be retrieved. If `oprob[:C]` is called a large number of times, this index must be found in each call. If a large number of such accesses are carried out, and performance is essential, it can be worthwhile to pre-compute a function to carry this out.
127+
128+
There exist four different functions, each returning a function for performing a specific type of interfacing:
129+
- `ModelingToolkit.getu`: For accessing species values.
130+
- `ModelingToolkit.getp`: For accessing parameter values.
131+
- `ModelingToolkit.setu`: For changing species values.
132+
- `ModelingToolkit.setp`: For changing parameter values.
133+
134+
For each species (or parameter) we wish to interface with, a new interfacing function must be created. Here we first creates a function for retrieving the value of $C$, and then use it for this purpose:
135+
```@example structure_indexing
136+
get_C = ModelingToolkit.getu(oprob, :C)
137+
get_C(oprob)
138+
```
139+
Here, `getu` (as well as `getp`, `setu`, and `setp`) first takes the structure we wish to interface with, and then the target quantity. When using `setu` and `setp`, in the second step, we must also provide the update value:
140+
```@example structure_indexing
141+
set_C = ModelingToolkit.setu(oprob, :C)
142+
set_C(oprob, 0.2)
143+
get_C(oprob)
118144
```
119-
and then uses these to build systems of equations. Here, these symbolic variables (`x`, `y`, and `z`) and parameters (`σ`, `ρ`, and `β`) can be used to interface a `problem`, `integrator`, and `solution` object (like we did previously, but using Symbols, e.g. `:X`). Since Catalyst models are built on ModelingToolkit, these models also contain similar symbolic variables and parameters.
120-
```@example ex2
145+
146+
Like when indexing-based interfacing is used, these functions also work with vectors:
147+
```@example structure_indexing
148+
get_S = ModelingToolkit.getu(oprob, [:S₁, :S₂])
149+
get_S(oprob)
150+
```
151+
152+
## [Interfacing using symbolic representations](@id simulation_structure_interfacing_symbolic_representation)
153+
As [previously described](@ref ref), when e.g. [programmatic modelling is used], species and parameters can be represented as *symbolic variables*. These can be used to index a problem, just like `Symbol '-based representations can. Here we create a simple [two-state model](@ref ref) programmatically, and use its symbolic variables to check, and update, an initial condition:
154+
```@example structure_indexing_symbolic_variables
121155
using Catalyst
122-
rn = @reaction_network begin
123-
(k1,k2), X1 <--> X2
124-
end
156+
t = default_t()
157+
@species X1(t) X2(t)
158+
@parameters k1 k2
159+
rxs = [
160+
Reaction(k1, [X1], [X2]),
161+
Reaction(k2, [X2], [X1])
162+
]
163+
@named two_state_model = ReactionSystem(rxs, t)
164+
two_state_model = complete(two_state_model)
165+
166+
u0 = [X1 => 2.0, X2 => 0.0]
167+
tspan = (0.0, 1.0)
168+
ps = [k1 => 1.0, k2 => 2.0]
169+
oprob = ODEProblem(two_state_model, u0, tspan, ps)
125170
126-
@unpack k1,k2,X1,X2 = rn
171+
oprob[X1] = 5.0
172+
oprob[X1]
127173
```
128-
Here, we first list the parameters and variables (for reaction systems the latter are typically species) we wish to import (in this case we select all, but we could select only a subset), next we denote from which model (here `rn`) from which we wish to import from. Next, these values can be used directly to interface with e.g. an `ODEProblem`:
129-
```@example ex2
130-
u0 = [X1 => 1.0, X2 => 5.0]
131-
p = [:k1 => 5.0, :k2 => 2.0]
132-
oprob = ODEProblem(rn, u0, (0.0,10.0), p)
174+
Symbolic variables can be used to access or update species or parameters for all the cases when `Symbol`s can (including when using `remake` or e.g. `getu`).
133175

134-
oprob.ps[k1]
176+
An advantage when quantities are represented as symbolic variables is that [symbolic expressions] can be formed and used to index a structure. E.g. here we check the combined initial concentration of $X$ ($X1 + X2$) in our two-state problem:
177+
```@example structure_indexing_symbolic_variables
178+
oprob[X1 + X2]
135179
```
136-
To interface with integrators and solutions we use a similar syntax.
137180

138-
Finally, instead of using `@unpack` to access a symbolic variable or parameter, we can access it directly using `rn.X1`, and thus access an unknown of our `ODEProblem` using
139-
```@example ex2
140-
oprob[rn.X1]
181+
Just like symbolic variables can be used to directly interface with a structure, symbolic variables stored in `ReactionSystem` models can be used:
182+
```@example structure_indexing_symbolic_variables
183+
oprob[two_state_model.X1 + two_state_model.X2]
141184
```
185+
This can be used to form symbolic expressions using model quantities when a model has been created using the DSL (as an alternative to [@unpack]
186+
(@ref ref)). Alternatively, [creating an observable](@ref ref), and then interface using its `Symbol` representation, is also possible.
187+
188+
!!! warn
189+
With interfacing with a simulating structure using symbolic variables stored in a `ReactionSystem` model, ensure that the [model is complete](@ref ref).

0 commit comments

Comments
 (0)