Merge 591ff99 into 151d00f

docs/src/api.md

@@ -272,6 +272,7 @@ hillar
 ```@docs
 Base.convert
 ModelingToolkit.structural_simplify
+set_default_noise_scaling
 ```
 
 ## Chemistry-related functionalities

docs/src/faqs.md

@@ -59,7 +59,7 @@ plot(sol; idxs = [A, B])
 ```
 
 ## How to disable rescaling of reaction rates in rate laws?
-As explained in the [Reaction rate laws used in simulations](@ref) section, for
+As explained in the [Reaction rate laws used in simulations](@ref introduction_to_catalyst_ratelaws) section, for
 a reaction such as `k, 2X --> 0`, the generated rate law will rescale the rate
 constant, giving `k*X^2/2` instead of `k*X^2` for ODEs and `k*X*(X-1)/2` instead
 of `k*X*(X-1)` for jumps. This can be disabled when directly `convert`ing a
@@ -105,7 +105,7 @@ parametric_stoichiometry) section.
 
 ## How to set default values for initial conditions and parameters?
 How to set defaults when using the `@reaction_network` macro is described in
-more detail [here](@ref dsl_description_defaults). There are several ways to do
+more detail [here](@ref dsl_advanced_options_default_vals). There are several ways to do
 this. Using the DSL, one can use the `@species` and `@parameters` options:
 ```@example faq3
 using Catalyst

docs/src/introduction_to_catalyst/introduction_to_catalyst.md

@@ -1,7 +1,7 @@
 # [Introduction to Catalyst](@id introduction_to_catalyst)
 In this tutorial we provide an introduction to using Catalyst to specify
 chemical reaction networks, and then to solve ODE, jump, and SDE models
-generated from them. At the end we show what mathematical rate laws and
+generated from them[1]. At the end we show what mathematical rate laws and
 transition rate functions (i.e. intensities or propensities) are generated by
 Catalyst for ODE, SDE and jump process models.
 
@@ -151,8 +151,7 @@ underlying problem.
     variable-based parameter mappings, `u₀symmap` and `psymmap`, while when
     directly passing `repressilator` we could use either those or the
     `Symbol`-based mappings, `u₀map` and `pmap`. `Symbol`-based mappings can
-    always be converted to `symbolic` mappings using [`symmap_to_varmap`](@ref),
-    see the [Basic Syntax](@ref basic_examples) section for more details.
+    always be converted to `symbolic` mappings using [`symmap_to_varmap`](@ref).
 
 
 !!! note

docs/src/inverse_problems/behaviour_optimisation.md

@@ -1,5 +1,5 @@
 # [Optimization for non-data fitting purposes](@id behaviour_optimisation)
-In previous tutorials we have described how to use [PEtab.jl](@ref petab_parameter_fitting) and [Optimization.jl](@ref optimization_parameter_fitting) for parameter fitting. This involves solving an optimisation problem (to find the parameter set yielding the best model-to-data fit). There are, however, other situations that require solving optimisation problems. Typically, these involve the creation of a custom cost function, which optimum can then be found using Optimization.jl. In this tutorial we will describe this process, demonstrating how parameter space can be searched to find values that achieve a desired system behaviour. A more throughout description on how to solve these problems is provided by [Optimization.jl's documentation](https://docs.sciml.ai/Optimization/stable/) and the literature[^1]. 
+In previous tutorials we have described how to use PEtab.jl and [Optimization.jl](@ref optimization_parameter_fitting) for parameter fitting. This involves solving an optimisation problem (to find the parameter set yielding the best model-to-data fit). There are, however, other situations that require solving optimisation problems. Typically, these involve the creation of a custom cost function, which optimum can then be found using Optimization.jl. In this tutorial we will describe this process, demonstrating how parameter space can be searched to find values that achieve a desired system behaviour. A more throughout description on how to solve these problems is provided by [Optimization.jl's documentation](https://docs.sciml.ai/Optimization/stable/) and the literature[^1]. 
 
 ## [Maximising the pulse amplitude of an incoherent feed forward loop](@id behaviour_optimisation_IFFL_example)
 Incoherent feedforward loops (network motifs where a single component both activates and deactivates a downstream component) are able to generate pulses in response to step inputs[^2]. In this tutorial we will consider such an incoherent feedforward loop, attempting to generate a system with as prominent a response pulse as possible.

docs/src/inverse_problems/global_sensitivity_analysis.md

@@ -1,6 +1,6 @@
 # [Global Sensitivity Analysis](@id global_sensitivity_analysis)
 *Global sensitivity analysis* (GSA) is used to study the sensitivity of a function's outputs with respect to its input[^1]. Within the context of chemical reaction network modelling it is primarily used for two purposes:
-- [When fitting a model's parameters to data](@ref petab_parameter_fitting), it can be applied to the cost function of the optimisation problem. Here, GSA helps determine which parameters do, and do not, affect the model's fit to the data. This can be used to identify parameters that are less relevant to the observed data.
+- When fitting a model's parameters to data, it can be applied to the cost function of the optimisation problem. Here, GSA helps determine which parameters do, and do not, affect the model's fit to the data. This can be used to identify parameters that are less relevant to the observed data.
 - [When measuring some system behaviour or property](@ref behaviour_optimisation), it can help determine which parameters influence that property. E.g. for a model of a biofuel-producing circuit in a synthetic organism, GSA could determine which system parameters have the largest impact on the total rate of biofuel production.
 
 GSA can be carried out using the [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) package. This tutorial contains a brief introduction of how to use it for GSA on Catalyst models, with [GlobalSensitivity providing a more complete documentation](https://docs.sciml.ai/GlobalSensitivity/stable/).
@@ -52,7 +52,7 @@ on the domain $10^β ∈ (-3.0,-1.0)$, $10^a ∈ (-2.0,0.0)$, $10^γ ∈ (-2.0,0
 
 !!! note
     We should make a couple of notes about the example above:
-    - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref optimization_parameter_fitting_logarithmic_scale), this is advantageous in the context of inverse problems such as this one.
+    - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref optimization_parameter_fitting_log_scale), this is advantageous in the context of inverse problems such as this one.
     - For GSA, where a function is evaluated a large number of times, it is ideal to write it as performant as possible. Hence, we initially create a base `ODEProblem`, and then apply the [`remake`](@ref simulation_structure_interfacing_problems_remake) function to it in each evaluation of `peak_cases` to generate a problem which is solved for that specific parameter set.
     - Again, as [previously described in other inverse problem tutorials](@ref optimization_parameter_fitting_basics), when exploring a function over large parameter spaces, we will likely simulate our model for unsuitable parameter sets. To reduce time spent on these, and to avoid excessive warning messages, we provide the `maxiters = 100000` and `verbose = false` arguments to `solve`.
     - As we have encountered in [a few other cases](@ref optimization_parameter_fitting_basics), the `gsa` function is not able to take parameter inputs of the map form usually used for Catalyst. Hence, as a first step in `peak_cases` we convert the parameter vector to this form. Next, we remember that the order of the parameters when we e.g. evaluate the GSA output, or set the parameter bounds, corresponds to the order used in `ps = [:β => p[1], :a => p[2], :γ => p[3]]`.

docs/src/inverse_problems/optimization_ode_param_fitting.md

@@ -129,13 +129,13 @@ If we from previous knowledge know that $kD = 0.1$, and only want to fit the val
 fixed_p_prob_generator(prob, p) = remake(prob; p = vcat(p[1], 0.1, p[2]))
 nothing # hide
 ```
-Here, it takes the `ODEProblem` (`prob`) we simulate, and the parameter set used, during the optimisation process (`p`), and creates a modified `ODEProblem` (by setting a customised parameter vector [using `remake`](@ref simulation_structure_interfacing_remake)). Now we create our modified loss function:
+Here, it takes the `ODEProblem` (`prob`) we simulate, and the parameter set used, during the optimisation process (`p`), and creates a modified `ODEProblem` (by setting a customised parameter vector [using `remake`](@ref simulation_structure_interfacing_problems_remake)). Now we create our modified loss function:
 ```@example diffeq_param_estim_1 
 loss_function_fixed_kD = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, data_vals), Optimization.AutoForwardDiff(); prob_generator = fixed_p_prob_generator, maxiters=10000, verbose=false, save_idxs=4)
 nothing # hide
 ```
 
-We can create an `OptimizationProblem` from this one like previously, but keep in mind that it (and its output results) only contains two parameter values ($k$* and $kP$):
+We can create an `OptimizationProblem` from this one like previously, but keep in mind that it (and its output results) only contains two parameter values ($k$ and $kP$):
 ```@example diffeq_param_estim_1 
 optprob_fixed_kD = OptimizationProblem(loss_function_fixed_kD, [1.0, 1.0])
 optsol_fixed_kD = solve(optprob_fixed_kD, Optim.NelderMead())
@@ -160,7 +160,7 @@ nothing # hide
 corresponds to the same true parameter values as used previously (`[:kB => 1.0, :kD => 0.1, :kP => 0.5]`).
 
 ## [Parameter fitting to multiple experiments](@id optimization_parameter_fitting_multiple_experiments)
-Say that we had measured our model for several different initial conditions, and would like to fit our model to all these measurements simultaneously. This can be done by first creating a [corresponding `EnsembleProblem`](@ref advanced_simulations_ensemble_problems). How to then create loss functions for these are described in more detail [here](https://docs.sciml.ai/DiffEqParamEstim/stable/tutorials/ensemble/).
+Say that we had measured our model for several different initial conditions, and would like to fit our model to all these measurements simultaneously. This can be done by first creating a [corresponding `EnsembleProblem`](@ref ensemble_simulations). How to then create loss functions for these are described in more detail [here](https://docs.sciml.ai/DiffEqParamEstim/stable/tutorials/ensemble/).
 
 ## [Optimisation solver options](@id optimization_parameter_fitting_solver_options)
 Optimization.jl supports various [optimisation solver options](https://docs.sciml.ai/Optimization/stable/API/solve/) that can be supplied to the `solve` command. For example, to set a maximum number of seconds (after which the optimisation process is terminated), you can use the `maxtime` argument:

docs/src/model_creation/compositional_modeling.md

@@ -139,7 +139,7 @@ nothing # hide
 Here we assume the user will pass in the repressor species as a ModelingToolkit
 variable, and specify a name for the network. We use Catalyst's interpolation
 ability to substitute the value of these variables into the DSL (see
-[Interpolation of Julia Variables](@ref dsl_description_interpolation_of_variables)). To make the repressilator we now make
+[Interpolation of Julia Variables](@ref dsl_advanced_options_symbolics_and_DSL_interpolation)). To make the repressilator we now make
 three genes, and then compose them together
 ```@example ex1
 t = default_t()

docs/src/model_creation/constraint_equations.md

@@ -52,7 +52,7 @@ end
 Notice, here we interpolated the variable `V` with `$V` to ensure we use the
 same symbolic unknown variable in the `rn` as we used in building `osys`. See the
 doc section on [interpolation of variables](@ref
-dsl_description_interpolation_of_variables) for more information.
+dsl_advanced_options_symbolics_and_DSL_interpolation) for more information.
 
 We can now merge the two systems into one complete `ReactionSystem` model using
 [`ModelingToolkit.extend`](@ref):

docs/src/model_creation/dsl_advanced.md

@@ -463,3 +463,77 @@ A reaction's metadata can be accessed using specific functions, e.g. `Catalyst.h
 rx = @reaction p, 0 --> X, [description="A production reaction"]
 Catalyst.getdescription(rx)
 ```
+
+## [Working with symbolic variables and the DSL](@id dsl_advanced_options_symbolics_and_DSL)
+We have previously described how Catalyst represents its models symbolically (enabling e.g. symbolic differentiation of expressions stored in models). While Catalyst utilises this for many internal operation, these symbolic representations can also be accessed and harnessed by the user. Primarily, doing so is much easier during programmatic (as opposed to DSL-based) modelling. Indeed, the section on [programmatic modelling](@ref programmatic_CRN_construction) goes into more details about symbolic representation in models, and how these can be used. It is, however, also ways to utilise these methods during DSL-based modelling. Below we briefly describe two methods for doing so.
+
+### [Using `@unpack` to extract symbolic variables from `ReactionSystem`s](@id dsl_advanced_options_symbolics_and_DSL_unpack)
+Let us consider a simple [birth-death process](@ref basic_CRN_library_bd) created using the DSL:
+```@example dsl_advanced_programmatic_unpack
+using Catalyst # hide
+bd_model = @reaction_network begin
+    (p,d), 0 <--> X
+end
+nothing # hide
+```
+Since we have not explicitly declared `p`, `d`, and `X` using `@parameters` and `@species`, we cannot represent these symbolically (only using `Symbol`s). If we wish to do so, however, we can fetch these into our current scope using the `@unpack` macro:
+```@example dsl_advanced_programmatic_unpack
+@unpack p, d, X = bd_model
+nothing # hide
+```
+This lists first the quantities we wish to fetch (does not need to be the model's full set of parameters and species), then `=`, followed by the model variable. `p`, `d` and `X` are now symbolic variables in the current scope, just as if they had been declared using `@parameters` or `@species`. We can confirm this:
+```@example dsl_advanced_programmatic_unpack
+X
+```
+Next, we can now use these to e.g. designate initial conditions and parameter values for model simulations:
+```@example dsl_advanced_programmatic_unpack
+using OrdinaryDiffEq, Plots # hide
+u0 = [X => 0.1]
+tspan = (0.0, 10.0)
+ps = [p => 1.0, d => 0.2]
+oprob = ODEProblem(bd_model, u0, tspan, ps)
+sol = solve(oprob)
+plot(sol)
+```
+
+!!! warn
+    Just like when using `@parameters` and `@species`, `@unpack` will overwrite any variables in the current scope which share name with the imported quantities.
+
+### [Interpolating variables into the DSL](@id dsl_advanced_options_symbolics_and_DSL_interpolation)
+Catalyst's DSL allows Julia variables to be interpolated for the network name, within rate constant expressions, or for species/stoichiometries within reactions. Using the lower-level symbolic interface we can then define symbolic variables and parameters outside of `@reaction_network`, which can then be used within expressions in the DSL. 
+
+Interpolation is carried out by pre-appending the interpolating variable with a `$`. For example, here we declare the parameters and species of a birth-death model, and interpolate these into the model:
+```@example dsl_advanced_programmatic_interpolation
+using Catalyst # hide
+t = default_t()
+@species X(t)
+@parameters p d
+bd_model = @reaction_network begin
+    ($p, $d), 0 <--> $X
+end
+```
+Additional information (such as default values or metadata) supplied to `p`, `d`, and `X` is carried through to the DSL. However, interpolation for this purpose is of limited value, as such information [can be declared within the DSL](@ref dsl_advanced_options_declaring_species_and_parameters). However, it is possible to interpolate larger algebraic expressions into the DSL, e.g. here
+```@example dsl_advanced_programmatic_interpolation
+@species X1(t) X2(t) X3(t) E(t)
+@parameters d
+d_rate = d/(1 + E)
+degradation_model = @reaction_network begin
+    $d_rate, X1 --> 0
+    $d_rate, X2 --> 0
+    $d_rate, X3 --> 0
+end
+```
+we declare an expression `d_rate`, which then can be inserted into the DSL via interpolation.
+
+It is also possible to use interpolation in combination with the `@reaction` macro. E.g. the reactions of the above network can be declared individually using
+```@example dsl_advanced_programmatic_interpolation
+rxs = [
+    @reaction $d_rate, $X1 --> 0
+    @reaction $d_rate, $X2 --> 0
+    @reaction $d_rate, $X3 --> 0
+]
+nothing # hide 
+```
+
+!!! note
+    When using interpolation, expressions like `2$spec` won't work; the multiplication symbol must be explicitly included like `2*$spec`.

docs/src/model_simulation/ode_simulation_performance.md

@@ -8,7 +8,7 @@ Generally, this short checklist provides a quick guide for dealing with ODE perf
 1. If performance is not critical, use [the default solver choice](@ref ode_simulation_performance_solvers) and do not worry further about the issue.
 2. If improved performance would be useful, read about solver selection (both in [this tutorial](@ref ode_simulation_performance_solvers) and [OrdinaryDiffEq's documentation](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)) and then try a few different solvers to find one with good performance.
 3. If you have a large ODE (approximately 100 variables or more), try the [various options for efficient Jacobian computation](@ref ode_simulation_performance_jacobian) (noting that some are non-trivial to use, and should only be investigated if truly required).
-4. If you plan to simulate your ODE many times, try [parallelise it on CPUs or GPUs](@ref investigating) (with preference for the former, which is easier to use).
+4. If you plan to simulate your ODE many times, try [parallelise it on CPUs or GPUs](@ref ode_simulation_performance_parallelisation) (with preference for the former, which is easier to use).
 
 ## [Regarding stiff and non-stiff problems and solvers](@id ode_simulation_performance_stiffness)
 Generally, ODE problems can be categorised into [*stiff ODEs* and *non-stiff ODEs*](https://en.wikipedia.org/wiki/Stiff_equation). This categorisation is important due to stiff ODEs requiring specialised solvers. A common cause of failure to simulate an ODE is the use of a non-stiff solver for a stiff problem. There is no exact way to determine whether a given ODE is stiff or not, however, systems with several different time scales (e.g. a CRN with both slow and fast reactions) typically generate stiff ODEs.

docs/src/model_simulation/simulation_introduction.md

@@ -275,7 +275,7 @@ nothing # hide
 ```
 If the `@default_noise_scaling` option is used, that term is only applied to reactions *without* `noise_scaling` metadata.
 
-While the `@default_noise_scaling` option is unavailable for [programmatically created models](@ref programmatic_CRN_construction), the [`remake_reactionsystem`](@ref) function can be used to achieve a similar effect.
+While the `@default_noise_scaling` option is unavailable for [programmatically created models](@ref programmatic_CRN_construction), the `set_default_noise_scaling` function can be used to achieve a similar effect.
 
 ## [Performing jump simulations using stochastic chemical kinetics](@id simulation_intro_jumps)