From 31c4470039d0610028eb2311e06cd69a562597cb Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 11 Jul 2024 16:07:36 -0400 Subject: [PATCH] finish interfacing, ode, and jump doc sections. --- .../lattice_reaction_systems.md | 130 ++++-------------- ...ttice_simulation_structure_ interaction.md | 115 +++++++--------- .../spatial_jump_simulations.md | 25 +++- .../spatial_ode_simulations.md | 36 ++++- .../lattice_sim_struct_interfacing.jl | 45 +++--- 5 files changed, 158 insertions(+), 193 deletions(-) diff --git a/docs/src/spatial_modelling/lattice_reaction_systems.md b/docs/src/spatial_modelling/lattice_reaction_systems.md index eef69e2b37..afce905a3b 100644 --- a/docs/src/spatial_modelling/lattice_reaction_systems.md +++ b/docs/src/spatial_modelling/lattice_reaction_systems.md @@ -4,12 +4,12 @@ Catalyst supports the expansion of non-spatial `ReactionSystem` (created using e - Discrete spatial domains. - Constant-rate transportation reactions (a species moving spatially at constant rates). -Features which support is planned in future updates include: +Features for which support is planned in future updates include: - Models on continuous domains with automatic discretisation (these models can already be simulated if the user provides a discretisation). - SDE simulations. - Transportation reactions with non-constant rates and general spatial reactions. -This tutorial introduces spatial modelling on discrete domains, here called lattices. It describes the basics of creating and simulating such models. To do so, it uses ODE simulations as examples. Additional tutorials provide further details how to interact with [spatial simulation structures](@ref lattice_simulation_plotting) and [plot spatial simulations](@ref ref), and also provide further details on [ODE](@ref spatial_lattice_ode_simulations) and [jump](@ref spatial_lattice_jump_simulations) simulations, respectively. +This tutorial introduces spatial modelling on discrete domains, here called lattices. It describes the basics of creating and simulating such models. To do so, it uses ODE simulations as examples. Additional tutorials provide further details on how to interact with [spatial simulation structures](@ref lattice_simulation_plotting) and [plot spatial simulations](@ref lattice_simulation_plotting), and also provide further details on [ODE](@ref spatial_lattice_ode_simulations) and [jump](@ref spatial_lattice_jump_simulations) simulations, respectively. ## [Basic example of spatial simulations on a discrete domain](@id spatial_lattice_modelling_intro_example) @@ -27,12 +27,12 @@ brusselator = @reaction_network begin B, X --> Y 1, X --> ∅ end -diffusion_reaction = @transport_reaction D X +diffusion_rx = @transport_reaction D X lattice = CartesianGrid((20,20)) -lrs = LatticeReactionSystem(brusselator, [diffusion_reaction], lattice) +lrs = LatticeReactionSystem(brusselator, [diffusion_rx], lattice) nothing # hide ``` -Here we define +This model contains: - A single spatial reaction, a transport reaction where $X$ moves at constant rate $D$ between adjacent compartments. - A 2d Cartesian grid of $20x20$ compartments to simulate our model on. @@ -46,14 +46,15 @@ ps = [:A => 1.0, :B => 4.0, :D => 0.2] oprob = ODEProblem(lrs, u0, tspan, ps) nothing # hide ``` -In this example we used non-uniform values for $X$'s initial condition, but uniform values for the remaining initial condition and parameter values. More details of uniform and non-uniform initial conditions and parameter values are provided [here](@ref spatial_lattice_modelling_intro_simulation_inputs). We also note that the diffusion reaction introduces a parameter $D$ (determining $X$'s diffusion rate) which value must designate in the parameter vector. +In this example we used non-uniform values for $X$'s initial condition, but uniform values for the remaining initial condition and parameter values. More details of uniform and non-uniform initial conditions and parameter values are provided [here](@ref spatial_lattice_modelling_intro_simulation_inputs). We also note that the diffusion reaction introduces a new parameter $D$ (determining $X$'s diffusion rate) whose value must be designated in the parameter vector. We can now simulate our model: ```@example spatial_intro_1 -sol = solve(oprob, QNDF()) +using OrdinaryDiffEq +sol = solve(oprob, FBDF()) nothing # hide ``` -We note that simulations of spatial models are often computationally expensive. Here we use the `QNDF` solver, which typically performs well for large systems. More advice on the performance of spatial simulations is provided [here](@ref ref). +We note that simulations of spatial models are often computationally expensive. Here we use the `FBDF` solver, which typically performs well for large systems. More advice on the performance of spatial ODE simulations is provided [here](@ref ref). Finally, we can access "$X$'s value across the simulation using ```@example spatial_intro_1 @@ -71,10 +72,10 @@ Spatial reactions describe reaction events which involve species across two conn - A rate at which it occurs. As for non-spatial reactions, this can be any expression. However, currently, it may only consist of parameters and other constants. - A single species which is transported from one compartment to an adjacent one. -At the occurrence of a transport reaction, the specific species moves to the adjacent compartment. Most common spatial models can be represented using transport reactions only. These can model phenomena such as diffusion or constant flux. A transportation reaction can be created using the `@transportation_reaction` macro. E.g. above we used +At the occurrence of a transport reaction, the specific species moves to the adjacent compartment. Many common spatial models can be represented using transport reactions only. These can model phenomena such as diffusion or constant flux. A transportation reaction can be created using the `@transportation_reaction` macro. E.g. above we used ```@example spatial_intro_3 using Catalyst # hide -diffusion_reaction = @transport_reaction D X +diffusion_rx = @transport_reaction D X nothing # hide ``` to create a reaction where species $X$ moves at a constant rate $D$ between adjacent compartments (this rate also scales linearly with the concentration of $X$ in the source compartment). Transport reactions may have rates depending on several parameters. E.g. to model a Brusselator where both species ($X$ and $Y$) are transported at a rate which depends both on the species, but also on some non-uniform parameter which is unique to each connection (e.g. representing the area connecting two cells in a tissue) we could do: @@ -89,22 +90,23 @@ If models are created [programmatically](@ref programmatic_CRN_construction) it ```@example spatial_intro_3 @variables t @species X(t) Y(t) -@parameters A B D +@parameters A B D [edgeparameter=true] tr_X = TransportationReaction(D, X) nothing # hide ``` +Note that in this example, we specifically designate $D$ an an [edge parameter](@ref spatial_lattice_modelling_intro_simulation_edge_parameters). ## [Defining discrete spatial domains (lattices)](@id spatial_lattice_modelling_intro_lattices) Discrete spatial domains can represent: -1. Systems which are composed of a (finite number of) compartments, where each compartment can be considered well-mixed (e.g. can be modelled non-spatially) and (potentially) species can move between different compartments. Tissues, where each compartment corresponds to a biological cell, are examples of such systems. -2. Systems that are continuous in nature, but have been approximated as a discrete domain. Future updates will include the ability for definition, and automatic discretisation, of continuous domains. Currently, however, the user has to perform this discretisation themselves. +1. Systems which are composed of a (finite number of) compartments, where each compartment can be considered well-mixed (e.g. can be modelled non-spatially) and where (potentially) species can move between adjacent compartments. Tissues, where each compartment corresponds to a biological cell, are examples of such systems. +2. Systems that are continuous in nature, but have been approximated as a discrete domain. Future updates will include the ability for the definition, and automatic discretisation, of continuous domains. Currently, however, the user has to perform this discretisation themselves. Catalyst supports three distinct types of lattices: - Cartesian lattices. These are grids where each grid point corresponds to a compartment. Spatial transportation is permitted between adjacent compartments. - Masked lattices. In these grids, only a subset of the grid point actually corresponds to compartments. Spatial transportation is permitted between adjacent compartments. - Unstructured (or graph) lattices. These are defined by graphs, where vertices correspond to compartments and edges connect adjacent compartments. -Here, Cartesian lattices are a subset of the masked lattices, which are a subset of the unstructured lattices. If possible, it is advantageous to use as narrow a lattice definition as possible (this may both improve simulation performance and simplifies syntax). Cartesian and masked lattices can be defined as one, two, and three dimensional. By default, these lattices assume that diagonally neighbouring compartments are non-adjacent (do not permit direct movement of species in between themselves). To change this, provide the `diagonally_adjacent = true` argument to your `LatticeReactionSystem` when it is created. +Here, Cartesian lattices are a subset of the masked lattices, which are a subset of the unstructured lattices. If possible, it is advantageous to use as narrow a lattice definition as possible (this may both improve simulation performance and simplify syntax). Cartesian and masked lattices can be defined as one, two, and three-dimensional. By default, these lattices assume that diagonally neighbouring compartments are non-adjacent (do not permit direct movement of species in between themselves). To change this, provide the `diagonally_adjacent = true` argument to your `LatticeReactionSystem` when it is created. ### [Defining Cartesian lattices](@id spatial_lattice_modelling_intro_lattices_cartesian) A Cartesian lattice is defined using the `CartesianGrid` function, which takes a single argument. For a 1d grid, simply provide the length of the grid as a single argument: @@ -113,7 +115,7 @@ using Catalyst # hide cgrid_1d = CartesianGrid(5) nothing # hide ``` -For 2d and 3d grids, instead provide a Tuple with the length of the grid in each dimension: +For 2d and 3d grids, we instead provide a Tuple with the length of the grid in each dimension: ```@example spatial_intro_3 cgrid_2d = CartesianGrid((3, 9)) cgrid_3d = CartesianGrid((2, 4, 8)) @@ -126,7 +128,7 @@ Masked lattices are defined through 1d, 2d, or 3d Boolean arrays. Each position rgrid_1d = [true, true, true, false, true, true, true] nothing # hide ``` -To define a 2d grid corresponding to the shape of an "8", we can use: +To define a 2d grid corresponding to the shape of an (laying) "8", we can use: ```@example spatial_intro_3 rgrid_2d = [ true true true true true; @@ -142,13 +144,13 @@ nothing # hide ``` ### [Defining unstructured lattices](@id spatial_lattice_modelling_intro_lattices_graph) -To define unstructured lattices, we must first import the [Graphs.jl](https://github.com/JuliaGraphs/Graphs.jl) package. Next, we can either use some [pre-defined formula for building graphs](https://juliagraphs.org/Graphs.jl/dev/core_functions/simplegraphs_generators/#Generators-for-common-graphs), or [build a graph from scratch](https://juliagraphs.org/Graphs.jl/dev/first_steps/construction/). Here we create a cyclic graph (where each compartment is connected to exactly two other compartments): +To define unstructured lattices, we must first import the [Graphs.jl](https://github.com/JuliaGraphs/Graphs.jl) package. Next, we can either use some [pre-defined formula for building graphs](https://juliagraphs.org/Graphs.jl/stable/core_functions/simplegraphs_generators/#Generators-for-common-graphs), or [build a graph from scratch](https://juliagraphs.org/Graphs.jl/stable/first_steps/construction/). Here we create a cyclic graph (where each compartment is connected to exactly two other compartments): ```@example spatial_intro_3 cycle_graph(7) nothing # hide ``` -Since graphs can represent any network of connected compartments, they do not have dimensions (like Cartesian or masked lattices). Another feature of graph lattices is that they can have non-symmetric connections, i.e. pairs of compartments where spatial movement of species is only permitted in one direction (in practise, this can be done for Cartesian and masked lattices as well, by [defining non-uniform spatial rates](@ref spatial_lattice_modelling_intro_simulation_inputs) and setting them to zero in one direction). This can be done by using a [*directed graph*](https://juliagraphs.org/Graphs.jl/dev/algorithms/digraph/) as input. E.g. here we define a directed cyclic graph, where movement is only allowed in one direction of the cycle: +Since graphs can represent any network of connected compartments, they do not have dimensions (like Cartesian or masked lattices). Another feature of graph lattices is that they can have non-symmetric connections, i.e. pairs of compartments where spatial movement of species is only permitted in one direction (in practice, this can be done for Cartesian and masked lattices as well, by [defining non-uniform spatial rates](@ref spatial_lattice_modelling_intro_simulation_inputs) and setting them to zero in one direction). This can be done by using a [*directed graph*](https://juliagraphs.org/Graphs.jl/dev/algorithms/digraph/) as input. E.g. here we define a directed cyclic graph, where movement is only allowed in one direction of the cycle: ```@example spatial_intro_3 cycle_digraph(7) nothing # hide @@ -166,7 +168,7 @@ The initial condition will be $1.0$ for $X$ across compartments, and $2.0$ for $ Below we describe how to set non-uniform values in the various cases. ### [Non-uniform compartment values for Cartesian lattices](@id spatial_lattice_modelling_intro_simulation_inputs_cartesian) -To provide non-uniform values across a Cartesian lattices, simply provide the values in an array of the same dimension and size as the Cartesian lattice. E.g. for a $5x10$ Cartesian lattice: +To provide non-uniform values across a Cartesian lattice, simply provide the values in an array of the same dimension and size as the Cartesian lattice. E.g. for a $5x10$ Cartesian lattice: ```@example spatial_intro_4 ccart_lattices = CartesianGrid((5, 10)) nothing # hide @@ -179,12 +181,12 @@ nothing # hide Non-uniform values for parameters (which values are tied to compartments) are provided similarly. ### [Non-uniform compartment values for masked lattices](@id spatial_lattice_modelling_intro_simulation_inputs_masked) -Non-uniform values for masked lattices are provided in the same manner as for Cartesian lattices (however, values at coordinates that do not hold compartments are ignored). E.g. To provide random values for a masked lattices contained within a $5x10$ Cartesian lattices we can again set: +Non-uniform values for masked lattices are provided in the same manner as for Cartesian lattices (however, values at coordinates that do not hold compartments are ignored). E.g. To provide random values for a masked lattice contained within a $5x10$ Cartesian lattices we can again set: ```@example spatial_intro_4 [:X => rand(5, 10), :Y => 10.0] nothing # hide ``` -If we want, it is also possible to provide the values as a [*sparse array*](https://github.com/JuliaSparse/SparseArrays.jl) with values only in the coordinates corresponding to compartments. +If we want, it is also possible to provide the values as a [*sparse array*](https://github.com/JuliaSparse/SparseArrays.jl) with values only in the coordinates that corresponds to compartments. ### [Non-uniform compartment values for unstructured lattices](@id spatial_lattice_modelling_intro_simulation_inputs_graphs) In graphs (which are used to represent unstructured lattices) each vertex (i.e. compartment) has a specific index. To set non-uniform values for unstructured lattices, provide a vector where the i'th value corresponds to the value in the compartment with index i in the graph. E.g. for a graph with 5 vertexes, where we want $X$ to be zero in all compartments bar one (where it is $1.0$) we use: @@ -196,7 +198,7 @@ nothing # hide ### [Non-uniform values for edge-parameters](@id spatial_lattice_modelling_intro_simulation_inputs_edge_parameters) Adjacent compartments are connected by edges (with which compartments are connected by edges being defined by the lattice). For unstructured lattices, it is possible (if a directed graph was used) to have edges from one compartment to another, but not in the opposite direction. For a lattice with $N$ compartments, edge values are set by a $NxN$ matrix, where value $(i,j)$ corresponds to the parameter's values in the edge going *from* compartment i *to* compartment j. This matrix can be either [sparse or non-sparse](https://docs.julialang.org/en/v1/stdlib/SparseArrays/). In the latter cases, values corresponding to non-existing edges are ignored. -E.g. let's consider a 1d Cartesian lattices with 4 compartments. Here, an edge parameter's values are provided in a $4x4$ matrix. For [the Brusselator model described previously](@ref spatial_lattice_modelling_intro_example), $D$'s value was tied to edges. If we wish to set the value of $D$ to various values between $0.1$ and $0.4$ we can do: +Let's consider a 1d Cartesian lattice with 4 compartments. Here, an edge parameter's values are provided in a $4x4$ matrix. For [the Brusselator model described previously](@ref spatial_lattice_modelling_intro_example), $D$'s value was tied to edges. If we wish to set the value of $D$ to various values between $0.1$ and $0.4$ we can do: ```@example spatial_intro_4 ps = [:A => 1.0, :B => 4.0, :D => [ @@ -207,7 +209,7 @@ ps = [:A => 1.0, :B => 4.0, ] nothing # hide ``` -Here, the value at index $i,j$ correspond to $D$'s value in the edge from compartment $i$ to compartment $j$. `0.0` is used for elements that do not correspond to an edge. The [`make_edge_p_values`](@ref) and [`make_directed_edge_values`](@ref) provides convenient interfaces for generating non-uniform edge parameter values. +Here, the value at index $i,j$ corresponds to $D$'s value in the edge from compartment $i$ to compartment $j$. `0.0` is used for elements that do not correspond to an edge. The [`make_edge_p_values`](@ref) and [`make_directed_edge_values`](@ref) provide convenient interfaces for generating non-uniform edge parameter values. ## [Edge parameters and compartment parameters](@id spatial_lattice_modelling_intro_simulation_edge_parameters) Parameters can be divided into *edge parameters* and *compartment parameters* (initial condition values are always tied to compartments). Here, edge parameters have their values tied to edges, while compartment parameters have their values tied to compartments. All parameters that are part of the rates (or stoichiometries) of non-spatial reactions must be compartment parameters. Parameters that are part of spatial reactions can be either compartment parameters or edge parameters. When a spatial reaction's rate is computed, edge parameters fetch their values for from the edge of the transition, and compartment parameters from the compartment from which *the edge originates*. @@ -232,87 +234,13 @@ brusselator = @reaction_network begin B, X --> Y 1, X --> ∅ end -diffusion_reaction = @transport_reaction D X +diffusion_rx = @transport_reaction D X lattice = CartesianGrid((20,20)) -lrs = LatticeReactionSystem(brusselator, [diffusion_reaction], lattice) +lrs = LatticeReactionSystem(brusselator, [diffusion_rx], lattice) edge_parameters(lrs) ``` ## [Spatial modelling limitations](@id spatial_lattice_modelling_intro_limitations) -Many features which are supported for non-spatial `ReactionSystem`s are currently unsupported for `LatticeReactionSystem`s. This includes [observables](@ref dsl_advanced_options_observables), [algebraic and differential equations](@ref constraint_equations), [hierarchical models](@ref compositional_modeling), and [events](@ref constraint_equations_events). It is possible that these features will be supported in the future. +Many features which are supported for non-spatial `ReactionSystem`s are currently unsupported for `LatticeReactionSystem`s. This includes [observables](@ref dsl_advanced_options_observables), [algebraic and differential equations](@ref constraint_equations), [hierarchical models](@ref compositional_modeling), and [events](@ref constraint_equations_events). It is possible that these features will be supported in the future. Furthermore, [removal of conserved quantities](@ref network_analysis_deficiency) is not supported when creating spatial `ODEProblem`s. -If you are using Catalyst's features for spatial modelling, please give us feedback on how we can improve these features. Additionally, just letting us know that you use these features is useful, as it helps inform us whether continued development of spatial modelling features will be worthwhile. - - - - - - -## Interacting with spatial problems and integrators -We have previously described how to [check, and update, values stored in problems, integrators, and solution objects](@ref ref). Due to internal representations, these operations are more difficult for spatial model. Generally, the following limitations apply to spatial systems: -- Parameter values stored in `ODEProblem`s and corresponding integrators can be accessed as normal. -- While parameter values stored in `ODEProblem`s and corresponding integrators can be updated, this must typically be followed by calling the `rebuild_lat_internals!` function for the update to take effect. -- Parameter values stored in `DiscreteProblem`s, `JumpProblem`s, and corresponding integrators cannot be accessed nor updated. -- Species initial conditions stored in problems cannot be accessed nor updated. -- Species values stored in integrators cannot be updated, and can only be access through the specialised `lat_getu` function (described in more detail [here](@ref ref)). -- Species values stored in solution objects can only be accessed through the specialised `lat_getu` function. - -Below we describe these cases - -## Accessing and updating parameter values stored in problems and integrators -Let us consider a simple spatial [birth-death] process where $p$'ve value is non-uniform (while $d$'s value is uniform). -```@example spatial_intro_5 -using Catalyst -bd_model = @reaction_network begin - (p,d), 0 <--> X -end -tr = @transport_reaction D X -lattice = CartesianGrid((2,2)) -lrs = LatticeReactionSystem(rs, [tr], lattice) - -u0 = [:X => 0.1] -tspan = (0.0, 200.0) -ps = [:p => [1.1 1.2; 1.3 1.4], :d => 0.5, :D => 0.2] -oprob = ODEProblem(lrs, u0, tspan, ps) -nothing # hide -``` -Here, we can check the parameter values by indexing with the corresponding symbol: -```@example spatial_intro_5 -oprob.ps[:p] -``` -!!! note - For non-spatial models, a parameter's [symbolic representation](@ref ref) can be used for indexing. However, for spatial ones, the symbol form must be used. You can convert from the former to the latter through `ModelingToolkit.getname(p)`. - -We can update a parameter value using similar syntax: -```@example spatial_intro_5 -oprob.ps[:p] = [2.1 2.2; 2.3 2.4] -nothing # hide -``` -When assigning a uniform value to a parameter, *that value must be enclosed in a vector*. I.e. to update the value of `d` to `1.0` we must use -```@example spatial_intro_5 -oprob.ps[:d] = [1.0] -nothing # hide -``` -I.e. this -```@example spatial_intro_5 -oprob.ps[:d] = 1.0 -nothing # hide -``` -does not work. - -When updating parameter values in either of these two manners: -- Setting non-uniform values for a previously uniform parameter. -- Setting uniform values for a previously non-uniform parameter. -- Setting values for an [edge parameter](@ref spatial_lattice_modelling_intro_simulation_edge_parameters). - -You must rebuild the problem for the changes to have any effect. This is done through the `rebuild_lat_internals!` function. I.e. here we gives $d$ non-uniform values, and then rebuilds the `oprob`: -```@example spatial_intro_5 -oprob.ps[:d] = [0.4 0.5; 0.6 0.7] -rebuild_lat_internals!(oprob) -nothing # hide -``` - -Parameter values stored in [integrators](@ref ref) can be accessed and updated in the same manner as those stored in problems. This can primarily be used to update parameter values during a simulation as result of the [triggering of a callback](@ref ref). - -!!! warning - Unlike for `ODEProblem`s and their integrator, the `lat_setp!` and `rebuild_lat_internals!` functions are currently not supported for `JumpProblem`s and their integrators. \ No newline at end of file +If you are using Catalyst's features for spatial modelling, please give us feedback on how we can improve these features. Additionally, just letting us know that you use these features is useful, as it helps inform us whether continued development of spatial modelling features is worthwhile. diff --git a/docs/src/spatial_modelling/lattice_simulation_structure_ interaction.md b/docs/src/spatial_modelling/lattice_simulation_structure_ interaction.md index 9ab2e46e5f..a9ad9c5b4e 100644 --- a/docs/src/spatial_modelling/lattice_simulation_structure_ interaction.md +++ b/docs/src/spatial_modelling/lattice_simulation_structure_ interaction.md @@ -9,13 +9,13 @@ there are no equally straightforward interfaces for spatial simulations. Typical ```julia lat_getu(sol, :X, lrs) ``` -furthermore, there are some cases of interfacing which are currently not supported (e.g. updating parameter values in `JumpProblem`s). It is likely that these interfaces will be improved in the future (i.e. using a similar syntax as non-spatial ones). Finally, we note that many of the functions presented below have not been as extensively optimised for performance as other parts of the Catalyst package. Hence, you should take care when designing workflows which requires using them a large number of times. +Furthermore, there are some cases of interfacing which are currently not supported (e.g. updating parameter values in `JumpProblem`s). It is likely that these interfaces will be improved in the future (i.e. by introducing a similar syntax to the current non-spatial one). Finally, we note that many of the functions presented below have not been as extensively optimised for performance as other parts of the Catalyst package. Hence, you should take care when designing workflows which requires using them a large number of times. !!! note - Below we will describe various features using ODE simulations as examples. However, all interfaces (unless where else is stated) works identically for jump simulations. + Below we will describe various features using ODE simulations as examples. However, all interfaces (unless where else is stated) work identically for jump simulations. ## [Retrieving values from lattice simulations](@ref lattice_simulation_structure_interaction_simulation_species) -Let us consider a simulation from a `LatticeReactionSystem` +Let us consider a simulation of a `LatticeReactionSystem` ```@example lattice_struct_interaction_sims using Catalyst, OrdinaryDiffEq two_state_model = @reaction_network begin @@ -32,7 +32,7 @@ oprob = ODEProblem(lrs, u0, tspan, ps) sol = solve(oprob, Tsit5()) nothing # hide ``` -to retrieve the values of $X1$ across the simulation we use the `lat_getu` function. It takes three arguments: +To retrieve the values of $X1$ across the simulation we use the `lat_getu` function. It takes three arguments: - The solution objects from which we wish to retrieve values. - The species which values we wish to retrieve. - The `LatticeReactionSystem` which was simulated. @@ -40,19 +40,22 @@ to retrieve the values of $X1$ across the simulation we use the `lat_getu` funct ```@example lattice_struct_interaction_sims lat_getu(sol, :X1, lrs) ``` -Here, the output is a vector with $X$'s value at each simulation time step. How the specie's value is represented in each time step depends on the lattice which was used to create the `LatticeReactionSystem` originally: +Here, the output is a vector with $X$'s value at each simulation time step. How the specie's value is represented in each time step depends on the lattice which was originally used to create the `LatticeReactionSystem`: - For Cartesian lattices, an array of the same size as the Cartesian lattice is used. Each array element corresponds to the species's value in the corresponding compartment. -- For masked lattices, a [sparse array]https://docs.julialang.org/en/v1/stdlib/SparseArrays/ of the same size as the masked lattice is used. Each filled array element corresponds to the species's value in the corresponding compartment. Unfilled array elements corresponds to positions without compartments. -- For unstructured (graph) lattices, vectors are used. The i'th element in the vectors correspond to the species's value in the i'th compartment. +- For masked lattices, a [sparse array]https://docs.julialang.org/en/v1/stdlib/SparseArrays/ of the same size as the masked lattice is used. Each filled array element corresponds to the species's value in the corresponding compartment. Unfilled array elements correspond to positions without compartments. +- For unstructured (graph) lattices, vectors are used. The i'th element in the vectors corresponds to the species's value in the i'th compartment. Unlike for non-spatial simulations, `lat_getu` does not take vector (e.g. `lat_getu(sol, [:X1, :X2], lrs)`) or symbolic expression (e.g. `lat_getu(sol, [X1 + X2], lrs)`) inputs. However, it is possible to use symbolic variables as input (e.g. `lat_getu(sol, two_state_model.X1, lrs)`). ### [Retrieving lattice simulations values at specific time points](@ref lattice_simulation_structure_interaction_simulation_species_ts) +Just like for non-spatial solutions, it is possible to access the simulation's values at designated time points. This is possible even if the simulation did not stop at those specific values (in which case an interpolated value is returned). To do this, the desired time points to sample are provided as a vector to `lat_getu` using the optional argument `t`. E.g. here we retrieve the simulation's (interpolated) values at time points `0.5` and `0.75`: - +```@example lattice_struct_interaction_sims +lat_getu(sol, :X1, lrs; t = [0.5, 0.75]) +``` ## [Retrieving and updating species values in problems and integrators](@ref lattice_simulation_structure_interaction_prob_int_species) -Let us consider a spatial `ODEProblem` and its corresponding [integrator](@ref simulation_structure_interfacing_integrators). +Let us consider a spatial `ODEProblem` ```@example lattice_struct_interaction_prob_ints using Catalyst, OrdinaryDiffEq two_state_model = @reaction_network begin @@ -66,14 +69,13 @@ u0 = [:X1 => [0.0 0.0 0.0; 2.0 2.0 2.0], :X2 => 0.0] tspan = (0.0, 1.0) ps = [:k1 => 2.0, :k2 => 1.0, :D => 0.1] oprob = ODEProblem(lrs, u0, tspan, ps) -oint = init(oprob, Tsit5()) nothing # hide ``` -We can retrieve the species values stored in either `oprob` or `oint` using the `lat_getu` function. It uses [identical syntax as for simulations](@ref lattice_simulation_structure_interaction_simulation_species) (except that you cannot specify a time point). However, it returns a single set of species values (while for simulations it returns a vector across all time steps): +We can retrieve the species values stored in `oprob` using the `lat_getu` function. It uses [identical syntax as for simulations](@ref lattice_simulation_structure_interaction_simulation_species) (except that you cannot specify a time point). However, it returns a single set of species values (while for simulations it returns a vector across all time steps): ```@example lattice_struct_interaction_prob_ints lat_getu(oprob, :X1, lrs) ``` -again, the format used correspond to the lattice used to create the original `LatticeReactionSystem`. Here, even if a species have homogeneous values, the full format is used. +again, the format used corresponds to the lattice used to create the original `LatticeReactionSystem`. Here, even if a species has homogeneous values, the full format is used. ```@example lattice_struct_interaction_prob_ints lat_getu(oprob, :X2, lrs) ``` @@ -82,75 +84,54 @@ For both problems and integrators, species values can be updated using the `lat_ ```@example lattice_struct_interaction_prob_ints lat_setu!(oprob, :X1, lrs, [1.0 2.0 3.0; 4.0 5.0 6.0]) ``` -Here, the same format (which depends on the used lattice) is used for the species new values, as which is used when initially designating their initial conditions. I.e. to make $X1$'s initial condition values uniform we can call +Here, the same format (which depends on the used lattice) is used for the species's new values, as which is used when initially designating their initial conditions. I.e. to make $X1$'s initial condition values uniform we can call ```@example lattice_struct_interaction_prob_ints lat_setu!(oprob, :X1, lrs, 1.0) ``` !!! note - It is currently not possible to change a species value at a single compartment only. To do so, you'd first have to retrieve its values across all compartments using `lat_getu`, the modify this at the desired compartment, and then used this modified array as an input to `lat_setu!`. + It is currently not possible to change a species value at a single compartment only. To do so, you must first retrieve its values across all compartments using `lat_getu`, then modify this at the desired compartment, and then use the modified version as input to `lat_setu!`. + +Species values in [integrators](@ref simulation_structure_interfacing_integrators) can be interfaced with using identical syntax as for problems. ## [Retrieving and updating parameter values in problems and integrators](@ref lattice_simulation_structure_interaction_prob_int_parameters) -The retrieval and updating of parameter values works a bit differently as for species: -- It uses an interface much similar as that for non-spatial models. -- It is only supported for `ODEProblem`s and their integrators, not for jump-based ones. -- Uniform parameter values are returned enclosed in a vector (e.g. `[1.0]`), and must also be provided in this format. -- After a problem or integrator has been updated, the `rebuild_lat_internals!` function must be called upon it. -- Indexing can only be carried out using symbols (e.g. `:X`). - -Let us consider a simple spatial model where $p$'s value is non-uniform and $d$'s is uniform. -```@example spatial_intro_5 -using Catalyst -bd_model = @reaction_network begin - (p,d), 0 <--> X -end -tr = @transport_reaction D X -lattice = CartesianGrid((2,2)) -lrs = LatticeReactionSystem(rs, [tr], lattice) +Retrieval and updating of parameter values for problems and integrators works similarly as for species, but with the following differences: +- The `lat_getp` and `lat_setp!` functions are used. +- It is currently not possible to interface with parameter values of `JumpProblem`s and their integrators. +- After parameter values are modified, the `rebuild_lat_internals!` function must be applied before the problem/integrators can be used for further analysis. +- Updating of [edge parameters](@ref spatial_lattice_modelling_intro_simulation_edge_parameters) is limited and uses a different interface. -u0 = [:X => 0.1] -tspan = (0.0, 200.0) -ps = [:p => [1.1 1.2; 1.3 1.4], :d => 0.5, :D => 0.2] -oprob = ODEProblem(lrs, u0, tspan, ps) -nothing # hide +Let us consider the spatial `ODEProblem` we previously declared. We can check the value of $k1$ by using `lat_getp` +```@example lattice_struct_interaction_prob_ints +lat_getp(oprob, :k1, lrs) ``` -Here, we can check the parameter values by indexing with the corresponding symbol: -```@example spatial_intro_5 -oprob.ps[:p] +Next, we can update it using `lat_setp!` (here we also confirm that it now has the updated values): +```@example lattice_struct_interaction_prob_ints +lat_setp!(oprob, :k1, lrs, [1.0 2.0 3.0; 4.0 5.0 6.0]) +lat_getp(oprob, :k1, lrs) ``` -!!! note - For non-spatial models, a parameter's [symbolic representation](@ref programmatic_CRN_construction) can be used for indexing. However, for spatial ones, the symbol form must be used. You can convert from the former to the latter through `ModelingToolkit.getname(p)`. -We can update a parameter value using similar syntax: -```@example spatial_intro_5 -oprob.ps[:p] = [2.1 2.2; 2.3 2.4] -nothing # hide -``` -When assigning a uniform value to a parameter, *that value must be enclosed in a vector*. I.e. to update the value of `d` to `1.0` we must use -```@example spatial_intro_5 -oprob.ps[:d] = [1.0] -nothing # hide -``` -I.e. this -```@example spatial_intro_5 -oprob.ps[:d] = 1.0 -nothing # hide +If we now were to simulate `oprob`, the simulation would not take the updated value of $k1$ into account. For our changes to take effect we might first need to call `rebuild_lat_internals!` with `oprob` as an input +```@example lattice_struct_interaction_prob_ints +rebuild_lat_internals!(oprob) ``` -does not work. +There are two different circumstances when `rebuild_lat_internals!` must be called: +- When modifying the value of an [edge parameter](@ref spatial_lattice_modelling_intro_simulation_edge_parameters). +- When changing a parameter from having spatially uniform values to spatially non-uniform values, or the other way around. -When updating parameter values in either of these two manners: -- Setting non-uniform values for a previously uniform parameter. -- Setting uniform values for a previously non-uniform parameter. -- Setting values for an [edge parameter](@ref spatial_lattice_modelling_intro_simulation_edge_parameters). +Parameter values of integrators can be interfaced with just like for problems (this is primarily relevant when using [*callbacks*](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/)). Again, after doing so, the `rebuild_lat_internals!` function might need to be applied to the integrator. -You must rebuild the problem for the changes to have any effect. This is done through the `rebuild_lat_internals!` function. I.e. here we gives $d$ non-uniform values, and then rebuilds the `oprob`: -```@example spatial_intro_5 -oprob.ps[:d] = [0.4 0.5; 0.6 0.7] -rebuild_lat_internals!(oprob) +The `lat_getp` and `lat_setp!` functions cannot currently be applied to [edge parameters](@ref spatial_lattice_modelling_intro_simulation_edge_parameters). Instead, to access the value of an edge parameter, use +```@example lattice_struct_interaction_prob_ints +oprob.ps[:D] +``` +To update an edge parameter's value, use +```@example lattice_struct_interaction_prob_ints +oprob.ps[:D] = [0.2] nothing # hide ``` - -Parameter values stored in [integrators](@ref ref) can be accessed and updated in the same manner as those stored in problems. This can primarily be used to update parameter values during a simulation as result of the [triggering of a callback](@ref ref). - - +This interface is somewhat limited, and the following aspects should be noted: +- Edge parameter values can only be interfaced with if the edge parameter's value is spatially uniform. +- When accessing an (spatially uniform) edge parameter's value, its single value will be encapsulated in a vector. +- When setting an (spatially uniform) edge parameter's value, you must encapsulate the new value in a vector. diff --git a/docs/src/spatial_modelling/spatial_jump_simulations.md b/docs/src/spatial_modelling/spatial_jump_simulations.md index daa65f5905..6aa048cd4b 100644 --- a/docs/src/spatial_modelling/spatial_jump_simulations.md +++ b/docs/src/spatial_modelling/spatial_jump_simulations.md @@ -1 +1,24 @@ -# [Spatial jump simulations](@id spatial_lattice_jump_simulations) \ No newline at end of file +# [Spatial jump simulations](@id spatial_lattice_jump_simulations) +Our [introduction to spatial lattice simulations](@ref spatial_lattice_modelling_intro) has already described how to simulate `LatticeReactionSystem`s using ODEs. Jump simulations of `LatticeReactionSystem` are carried out using an almost identical approach. However, just like for non-spatial models, we must first create a `DiscreteProblem`, which is then used as input to our `JumpProblem`. Furthermore, a spatial [jump aggregator](@ref simulation_intro_jumps_solver_designation) (like `NSM`) can be used. Spatial jump simulations in Catalyst are built on top of JumpProcesses.jl's spatial jump simulators, more details on which can be found [here](https://docs.sciml.ai/JumpProcesses/stable/tutorials/spatial/). + +Below we perform a spatial jump simulation of a simple [birth-death process](@ref basic_CRN_library_bd). Note that we use our `LatticeReactionSystem` as input to both our `DiscreteProblem` and `JumpProblem`, and that we [provide] the spatial `NSM` jump aggregator to `JumpProblem`. +```@example spatial_jump +using Catalyst, JumpProcesses +bd_model = @reaction_network begin + (p,d), 0 <--> X +end +diffusion_rx = @transport_reaction D X +lattice = CartesianGrid((10,10)) +lrs = LatticeReactionSystem(bd_model, [diffusion_rx], lattice) + +u0 = [:X => 0] +tspan = (0.0, 200.0) +ps = [:p => 10.0, :d => 1.0, :D => 0.1] +dprob = DiscreteProblem(lrs, u0, tspan, ps) +jprob = JumpProblem(lrs, dprob, NSM()) +sol = solve(jprob, SSAStepper()) +nothing # hide +``` +We can now access the values of `sol` using the interfaces described [here](@ref lattice_simulation_structure_interaction_simulation_species), and plot it using the functions described [here](@ref lattice_simulation_plotting). + +Currently, the only available spatial jump aggregators are `NSM` and `DirectCRDirect`, with `DirectCRDirect` expected to perform better for large networks. It is also possible to use non-spatial solvers (like `Direct` and `RSSA`) when creating spatial `JumpProblem`s. \ No newline at end of file diff --git a/docs/src/spatial_modelling/spatial_ode_simulations.md b/docs/src/spatial_modelling/spatial_ode_simulations.md index 87b9d61caf..f3e0c992c5 100644 --- a/docs/src/spatial_modelling/spatial_ode_simulations.md +++ b/docs/src/spatial_modelling/spatial_ode_simulations.md @@ -1 +1,35 @@ -# [Spatial ODE simulations](@id spatial_lattice_ode_simulations) \ No newline at end of file +# [Spatial ODE simulations](@id spatial_lattice_ode_simulations) +Our [introduction to spatial lattice simulations](@ref spatial_lattice_modelling_intro) has already provided an extensive description of how to simulate `LatticeReactionSystem`s using ODEs. Further tutorials have also shown how to [retrieve values from simulations](@ref lattice_simulation_structure_interaction_simulation_species), or how to [plot them](@ref lattice_simulation_plotting). Here we will build on this, primarily discussing strategies for increasing ODE simulation performance. This is especially important for spatial simulations, as these typically are more computationally demanding as compared to non-spatial ones. While focusing on non-spatial simulations, this [ODE performance tutorial](@ref ode_simulation_performance) can also be useful to read. + +## [Solver selection for spatial ODE simulations](@id spatial_lattice_ode_simulations_solvers) +Previously we have described [how to select ODE solvers, and how this can impact simulation performance](@ref ode_simulation_performance_solvers). This is especially relevant for spatial simulations. For stiff problems, `FBDF` is a good first solver to try. + +## [Jacobian options for spatial ODE simulations](@id spatial_lattice_ode_simulations_jacobians) +We have previously described how, when [implicit solvers are used to solve stiff ODEs](@ref ode_simulation_performance_stiffness), the [strategy for computing the system Jacobian](@ref ode_simulation_performance_jacobian) is important. This is especially the case for spatial simulations, where the Jacobian often is large and highly sparse. Catalyst implements special methods for spatial Jacobians. To utilise these, provide the `jac = true` argument to your `ODEProblem` when it is created (if `jac = false`, which is the default, [*automatic differentiation*](https://en.wikipedia.org/wiki/Automatic_differentiation) will be used for Jacobian computation). Here we simulate a [Brusselator](@ref basic_CRN_library_brusselator) while designating to use Catalyst's computed Jacobian: +```@example spatial_ode +using Catalyst, OrdinaryDiffEq +brusselator = @reaction_network begin + A, ∅ --> X + 1, 2X + Y --> 3X + B, X --> Y + 1, X --> ∅ +end +diffusion_rx = @transport_reaction D X +lattice = CartesianGrid((20,20)) +lrs = LatticeReactionSystem(brusselator, [diffusion_rx], lattice) + +u0 = [:X => rand(20, 20), :Y => 10.0] +tspan = (0.0, 200.0) +ps = [:A => 1.0, :B => 4.0, :D => 0.2] +oprob = ODEProblem(lrs, u0, tspan, ps; jac = true) +sol = solve(oprob, FBDF()) +nothing # hide +``` +For large systems, building a dense Jacobian can be problematic, in which case a [*sparse*](@ref ode_simulation_performance_sparse_jacobian) Jacobian can be designated using `sparse = true`: +```@example spatial_ode +oprob = ODEProblem(lrs, u0, tspan, ps; jac = true, sparse = true) +sol = solve(oprob, FBDF()) +nothing # hide +``` + +It is possible to use `sparse = true` while `jac = false`, in which case a sparse Jacobian is computed using automatic differentiation. \ No newline at end of file diff --git a/src/spatial_reaction_systems/lattice_sim_struct_interfacing.jl b/src/spatial_reaction_systems/lattice_sim_struct_interfacing.jl index 7f3cac9b11..9b0b4b9023 100644 --- a/src/spatial_reaction_systems/lattice_sim_struct_interfacing.jl +++ b/src/spatial_reaction_systems/lattice_sim_struct_interfacing.jl @@ -1,18 +1,17 @@ ### Lattice Simulation Structure Species Getters/Setters ### - """ lat_setp!(sim_struct, p, lrs::LatticeReactionSystem, p_val) -For a problem or integrators, updates its `p` vector with the input `p_val`. For non-lattice models, +For a problem or integrators, update its `p` vector with the input `p_val`. For non-lattice models, this is can be done through direct interfacing (e.g. `prob[p] = 1.0`). However, for `LatticeReactionSystem`-based problems and integrators, this function must be used instead. Arguments: - `sim_struct`: The simulation structure which `u` value we wish to update. Can be either a `ODEProblem`, `JumpProblem`, or an integrator derived from either of these. -- `p`: The species which value we wish to update. Can be provided either in its symbolic form, or -as a symbol. +- `p`: The species which value we wish to update. Can be provided either in its symbolic form or +a symbol. - `lrs`: The `LatticeReactionSystem` which was used to generate the structure we wish to modify. - `p_val`: The parameter's new values. Must be given in a form which is also a valid initial input to the `ODEProblem`/`JumpProblem`. @@ -43,7 +42,7 @@ function lat_setp!(sim_struct, p, lrs::LatticeReactionSystem, p_vals) (p_vals isa Number) || check_lattice_format(extract_lattice(lrs), p_vals) edge_param_check(p, lrs) - # Converts symbol parameter to symbolic and find correct species index and numbers. + # Converts symbol parameter to symbolic and find the correct species index and numbers. (p isa Symbol) && (p = _symbol_to_var(lrs, p)) (p isa Num) && (p = Symbolics.unwrap(p)) p_idx, p_tot = get_p_idxs(p, lrs) @@ -54,7 +53,7 @@ function lat_setp!(sim_struct, p, lrs::LatticeReactionSystem, p_vals) end # Note: currently, `lat_setp!(oprob::ODEProblem, ...`) and `lat_setp!(SciMLBase.AbstractODEIntegrator, ...`) -# are identical and could be merged to a singe function. +# are identical and could be merged to a single function. function lat_setp!(oprob::ODEProblem, p_idx::Int64, p_vals, num_verts) if length(p_vals) == 1 foreach(idx -> (oprob.p[p_idx][idx] = p_vals[1]), 1:num_verts) @@ -89,14 +88,14 @@ For a problem or integrators, retrieves its `p` values. For non-lattice models, this is can be done through direct interfacing (e.g. `prob[p]`). However, for `LatticeReactionSystem`-based problems and integrators, this function must be used instead. The output format depends on the lattice (a dense array for cartesian grid lattices, a sparse array for -masked grid lattices, and a vector for graph lattices). This format is similar to which is used to +masked grid lattices, and a vector for graph lattices). This format is similar to what is used to designate parameter initial values. Arguments: - `sim_struct`: The simulation structure which `p` value we wish to retrieve. Can be either a `ODEProblem`, `JumpProblem`, or an integrator derived from either of these. -- `p`: The parameter which value we wish to update. Can be provided either in its symbolic form, or -as a symbol. +- `p`: The parameter which value we wish to update. Can be provided either in its symbolic form or +a symbol. - `lrs`: The `LatticeReactionSystem` which was used to generate the structure we wish to modify. Notes: @@ -150,15 +149,15 @@ end """ lat_setu!(sim_struct, sp, lrs::LatticeReactionSystem, u) -For a problem or integrators, updates its `u` vector with the input `u`. For non-lattice models, +For a problem or integrators, update its `u` vector with the input `u`. For non-lattice models, this is can be done through direct interfacing (e.g. `prob[X] = 1.0`). However, for `LatticeReactionSystem`-based problems and integrators, this function must be used instead. Arguments: - `sim_struct`: The simulation structure which `u` value we wish to update. Can be either a `ODEProblem`, `JumpProblem`, or an integrator derived from either of these. -- `sp`: The species which value we wish to update. Can be provided either in its symbolic form, or -as a symbol. +- `sp`: The species which value we wish to update. Can be provided either in its symbolic form or +a symbol. - `lrs`: The `LatticeReactionSystem` which was used to generate the structure we wish to modify. - `u`: The species's new values. Must be given in a form which is also a valid initial input to the `ODEProblem`/`JumpProblem`. @@ -188,7 +187,7 @@ function lat_setu!(sim_struct, sp, lrs::LatticeReactionSystem, u) # Checks that if u is non-uniform, it has the correct format for the system's lattice. (u isa Number) || check_lattice_format(extract_lattice(lrs), u) - # Converts symbol species to symbolic and find correct species index and numbers. + # Converts symbol species to symbolic and finds correct species index and numbers. (sp isa Symbol) && (sp = _symbol_to_var(lrs, sp)) (sp isa Num) && (sp = Symbolics.unwrap(sp)) sp_idx, sp_tot = get_sp_idxs(sp, lrs) @@ -242,8 +241,8 @@ designate species initial conditions. Arguments: - `sim_struct`: The simulation structure which `u` value we wish to retrieve. Can be either a `ODEProblem`, `JumpProblem`, or an integrator derived from either of these. -- `sp`: The species which value we wish to update. Can be provided either in its symbolic form, or -as a symbol. +- `sp`: The species which value we wish to update. Can be provided either in its symbolic form or +a symbol. - `lrs`: The `LatticeReactionSystem` which was used to generate the structure we wish to modify. Notes: @@ -312,8 +311,8 @@ Arguments: - `sp`: The species which values we wish to retrieve. Can be either a symbol (e.g. `:X`) or a symbolic variable (e.g. `X`). - `lrs`: The `LatticeReactionSystem` which was simulated to generate the solution. -- `t = nothing`: If `nothing`, we simply return the solution across all saved time steps. If `t` -instead is a vector (or range of values), returns the solutions interpolated at these time points. +- `t = nothing`: If `nothing`, we simply return the solution across all saved time steps (default). +If `t` instead is a vector (or range of values), returns the solution interpolated at these time points. Notes: - The `lat_getu` is not optimised for performance. However, it should still be quite performant, @@ -351,7 +350,7 @@ end # across all sample points). function lat_getu(sol::ODESolution, lattice, t::Nothing, sp_idx::Int64, sp_tot::Int64) # ODE simulations contain, in each data point, all values in a single vector. Jump simulations - # instead in a matrix (NxM, where N is the number of species and M the number of vertices). We + # instead in a matrix (NxM, where N is the number of species and M is the number of vertices). We # must consider each case separately. if sol.prob isa ODEProblem return [reshape_vals(vals[sp_idx:sp_tot:end], lattice) for vals in sol.u] @@ -366,13 +365,13 @@ end # value at all designated time points. function lat_getu(sol::ODESolution, lattice, t::AbstractVector{T}, sp_idx::Int64, sp_tot::Int64) where {T <: Number} - # Checks that an appropriate `t` is provided (however, DiffEq does permit out of range `t`s). + # Checks that an appropriate `t` is provided (however, DiffEq does permit out-of-range `t`s). if (minimum(t) < sol.t[1]) || (maximum(t) > sol.t[end]) error("The range of the t values provided for sampling, ($(minimum(t)),$(maximum(t))) is not fully within the range of the simulation time span ($(sol.t[1]),$(sol.t[end])).") end # ODE simulations contain, in each data point, all values in a single vector. Jump simulations - # instead in a matrix (NxM, where N is the number of species and M the number of vertices). We + # instead in a matrix (NxM, where N is the number of species and M is the number of vertices). We # must consider each case separately. if sol.prob isa ODEProblem return [reshape_vals(sol(ti)[sp_idx:sp_tot:end], lattice) for ti in t] @@ -388,7 +387,7 @@ end """ rebuild_lat_internals!(sciml_struct) -Rebuilds the internal functions for simulating a LatticeReactionSystem. Wenever a problem or +Rebuilds the internal functions for simulating a LatticeReactionSystem. Whenever a problem or integrator has had its parameter values updated, this function should be called for the update to be taken into account. For ODE simulations, `rebuild_lat_internals!` needs only to be called when - An edge parameter has been updated. @@ -573,7 +572,7 @@ end # Throws an error when interfacing with an edge parameter. function edge_param_check(p, lrs) (p isa Symbol) && (p = _symbol_to_var(lrs, p)) - if isedgeparameter(p) + if isedgeparameter(p) throw(ArgumentError("The `lat_getp` and `lat_setp!` functions currently does not support edge parameter updating. If you require this functionality, please raise an issue on the Catalyst GitHub page and we can add this feature.")) end -end \ No newline at end of file +end