From 7009b98eacce42b9f013ddfea592e9bec242f93c Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Thu, 11 Jul 2024 11:18:34 -0400 Subject: [PATCH 1/2] doc updates --- HISTORY.md | 271 +++++++++++++++++++------------- docs/src/api.md | 1 + docs/src/v14_migration_guide.md | 41 +++-- src/dsl.jl | 3 +- 4 files changed, 198 insertions(+), 118 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index 753403e7af..ff92e9b5f5 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -4,18 +4,44 @@ ## Catalyst 14.0 -#### Breaking changes -Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which introduced several breaking changes to Catalyst. A summary of these (and how to handle them) can be found [here](https://docs.sciml.ai/Catalyst/stable/v14_migration_guide/). These are briefly summarised in the following bullet points: -- `ReactionSystem`s must now be marked *complete* before they are exposed to most forms of simulation and analysis. With the exception of `ReactionSystem`s created through the `@reaction_network` macro, all `ReactionSystem`s are *not* marked complete upon construction. The `complete` function can be used to mark `ReactionSystem`s as complete. To construct a `ReactionSystem` that is not marked complete via the DSL the new `@network_component` macro can be used. -- The `states` function has been replaced with `unknowns`. The `get_states` function has been replaced with `get_unknowns`. -- Support for most units (with the exception of `s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`) has currently been dropped by ModelingToolkit, and hence they are unavailable via Catalyst too. Its is expected that eventually support for relevant chemical units such as molar will return to ModelingToolkit (and should then immediately work in Catalyst too). -- Problem parameter values are now accessed through `prob.ps[p]` (rather than `prob[p]`). -- A significant bug prevents the safe application of the `remake` function on problems for which `remove_conserved = true` was used when updating the values of initial conditions. Instead, the values of each conserved constant must be directly specified. -- The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been deprecated and removed. -- To be more consistent with ModelingToolkit's immutability requirement for systems, we have removed API functions that mutate `ReactionSystem`s such as `addparam!`, `addreaction!`, `addspecies`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystem`s from multiple component systems. +#### Breaking changes +Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which +introduced several breaking changes to Catalyst. A summary of these (and how to +handle them) can be found +[here](https://docs.sciml.ai/Catalyst/stable/v14_migration_guide/). These are +briefly summarised in the following bullet points: +- `ReactionSystem`s must now be marked *complete* before they are exposed to + most forms of simulation and analysis. With the exception of `ReactionSystem`s + created through the `@reaction_network` macro, all `ReactionSystem`s are *not* + marked complete upon construction. The `complete` function can be used to mark + `ReactionSystem`s as complete. To construct a `ReactionSystem` that is not + marked complete via the DSL the new `@network_component` macro can be used. +- The `states` function has been replaced with `unknowns`. The `get_states` + function has been replaced with `get_unknowns`. +- Support for most units (with the exception of `s`, `m`, `kg`, `A`, `K`, `mol`, + and `cd`) has currently been dropped by ModelingToolkit, and hence they are + unavailable via Catalyst too. Its is expected that eventually support for + relevant chemical units such as molar will return to ModelingToolkit (and + should then immediately work in Catalyst too). +- Problem parameter values are now accessed through `prob.ps[p]` (rather than + `prob[p]`). +- ModelingToolkit currently does not support the safe application of the + `remake` function, or safe direct mutation, for problems for which + `remove_conserved = true` was used when updating the values of initial + conditions. Instead, the values of each conserved constant must be directly + specified. +- The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions + have been deprecated and removed. +- To be more consistent with ModelingToolkit's immutability requirement for + systems, we have removed API functions that mutate `ReactionSystem`s such as + `addparam!`, `addreaction!`, `addspecies`, `@add_reactions`, and `merge!`. + Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate + new merged and/or composed `ReactionSystem`s from multiple component systems. #### General changes -- The `default_t()` and `default_time_deriv()` functions are now the preferred approaches for creating the default time independent variable and its differential. i.e. +- The `default_t()` and `default_time_deriv()` functions are now the preferred + approaches for creating the default time independent variable and its + differential. i.e. ```julia # do t = default_t() @@ -25,110 +51,141 @@ Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which @variables t @species A(t) - It is now possible to add metadata to individual reactions, e.g. using: -```julia -rn = @reaction_network begin - @parameters η - k, 2X --> X2, [description="Dimerisation"] -end -getdescription(rn) -``` -a more detailed description can be found [here](https://docs.sciml.ai/Catalyst/dev/model_creation/dsl_advanced/#dsl_advanced_options_reaction_metadata). -- `SDEProblem` no longer takes the `noise_scaling` argument. Noise scaling is now handled through the `noise_scaling` metadata (described in more detail [here](https://docs.sciml.ai/Catalyst/stable/model_simulation/simulation_introduction/#simulation_intro_SDEs_noise_saling)) -- Fields of the internal `Reaction` structure have been changed. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions. -- A new function, `save_reactionsystem`, which permits the writing of `ReactionSystem` models to files, has been created. A thorough description of this function can be found [here](https://docs.sciml.ai/Catalyst/stable/model_creation/model_file_loading_and_export/#Saving-Catalyst-models-to,-and-loading-them-from,-Julia-files) -- Update how compounds are created. E.g. use -```julia -@variables t C(t) O(t) -@compound CO2 ~ C + 2O -``` -to create a compound species `CO2` that consists of `C` and 2 `O`. -- Added documentation for chemistry-related functionality (compound creation and reaction balancing). + ```julia + rn = @reaction_network begin + @parameters η + k, 2X --> X2, [description="Dimerisation"] + end + getdescription(rn) + ``` + a more detailed description can be found [here](https://docs.sciml.ai/Catalyst/dev/model_creation/dsl_advanced/#dsl_advanced_options_reaction_metadata). +- `SDEProblem` no longer takes the `noise_scaling` argument. Noise scaling is + now handled through the `noise_scaling` metadata (described in more detail + [here](https://docs.sciml.ai/Catalyst/stable/model_simulation/simulation_introduction/#simulation_intro_SDEs_noise_saling)) +- Fields of the internal `Reaction` structure have been changed. + `ReactionSystems`s saved using `serialize` on previous Catalyst versions + cannot be loaded using this (or later) versions. +- A new function, `save_reactionsystem`, which permits the writing of + `ReactionSystem` models to files, has been created. A thorough description of + this function can be found + [here](https://docs.sciml.ai/Catalyst/stable/model_creation/model_file_loading_and_export/#Saving-Catalyst-models-to,-and-loading-them-from,-Julia-files) +- Updated how compounds are created. E.g. use + ```julia + @variables t C(t) O(t) + @compound CO2 ~ C + 2O + ``` + to create a compound species `CO2` that consists of `C` and two `O`. +- Added documentation for chemistry-related functionality (compound creation and + reaction balancing). - Added function `isautonomous` to check if a `ReactionSystem` is autonomous. -- Added function `steady_state_stability` to compute stability for steady states. Example: -```julia -# Creates model. -rn = @reaction_network begin - (p,d), 0 <--> X -end -p = [:p => 1.0, :d => 0.5] +- Added function `steady_state_stability` to compute stability for steady + states. Example: + ```julia + # Creates model. + rn = @reaction_network begin + (p,d), 0 <--> X + end + p = [:p => 1.0, :d => 0.5] -# Finds (the trivial) steady state, and computes stability. -steady_state = [2.0] -steady_state_stability(steady_state, rn, p) -``` -Here, `steady_state_stability` takes an optional argument `tol = 10*sqrt(eps())`, which is used to check that the real part of all eigenvalues are at least `tol` away from zero. Eigenvalues within `tol` of zero indicate that stability may not be reliably calculated. -- Added a DSL option, `@combinatoric_ratelaws`, which can be used to toggle whether to use combinatorial rate laws within the DSL (this feature was already supported for programmatic modelling). Example: -```julia -# Creates model. -rn = @reaction_network begin - @combinatoric_ratelaws false - (kB,kD), 2X <--> X2 -end -``` -- Added a DSL option, `@observables` for [creating observables](https://docs.sciml.ai/Catalyst/stable/model_creation/dsl_advanced/#dsl_advanced_options_observables) (this feature was already supported for programmatic modelling). -- Added DSL options `@continuous_events` and `@discrete_events` to add events to a model as part of its creation (this feature was already supported for programmatic modelling). Example: -```julia -rn = @reaction_network begin - @continuous_events begin - [X ~ 1.0] => [X ~ X + 1.0] - end - d, X --> 0 -end -``` -- Added DSL option `@equations` to add (algebraic or differential) equations to a model as part of its creation (this feature was already supported for programmatic modelling). Example: -```julia -rn = @reaction_network begin - @equations begin - D(V) ~ 1 - V - end - (p/V,d/V), 0 <--> X -end -``` -- Coupled reaction network + differential equation (or algebraic differential equation) systems can now be converted to `SDESystem`s and `NonlinearSystem`s. - -#### Structural identifiability extension -- Added CatalystStructuralIdentifiabilityExtension, which permits StructuralIdentifiability.jl function to be applied directly to Catalyst systems. E.g. use -```julia -using Catalyst, StructuralIdentifiability -goodwind_oscillator = @reaction_network begin - (mmr(P,pₘ,1), dₘ), 0 <--> M - (pₑ*M,dₑ), 0 <--> E - (pₚ*E,dₚ), 0 <--> P -end -assess_identifiability(goodwind_oscillator; measured_quantities=[:M]) -``` -to assess (global) structural identifiability for all parameters and variables of the `goodwind_oscillator` model (under the presumption that we can measure `M` only). -- Automatically handles conservation laws for structural identifiability problems (eliminates these internally to speed up computations). -- A more detailed of how this extension works can be found [here](https://docs.sciml.ai/Catalyst/stable/inverse_problems/structural_identifiability/). + # Finds (the trivial) steady state, and computes stability. + steady_state = [2.0] + steady_state_stability(steady_state, rn, p) + ``` + Here, `steady_state_stability` takes an optional keyword argument `tol = + 10*sqrt(eps())`, which is used to check that the real part of all eigenvalues + are at least `tol` away from zero. Eigenvalues within `tol` of zero indicate + that stability may not be reliably calculated. +- Added a DSL option, `@combinatoric_ratelaws`, which can be used to toggle + whether to use combinatorial rate laws within the DSL (this feature was + already supported for programmatic modelling). Example: + ```julia + # Creates model. + rn = @reaction_network begin + @combinatoric_ratelaws false + (kB,kD), 2X <--> X2 + end + ``` +- Added a DSL option, `@observables` for [creating + observables](https://docs.sciml.ai/Catalyst/stable/model_creation/dsl_advanced/#dsl_advanced_options_observables) + (this feature was already supported for programmatic modelling). +- Added DSL options `@continuous_events` and `@discrete_events` to add events to + a model as part of its creation (this feature was already supported for + programmatic modelling). Example: + ```julia + rn = @reaction_network begin + @continuous_events begin + [X ~ 1.0] => [X ~ X + 1.0] + end + d, X --> 0 + end + ``` +- Added DSL option `@equations` to add (algebraic or differential) equations to + a model as part of its creation (this feature was already supported for + programmatic modelling). Example: + ```julia + rn = @reaction_network begin + @equations begin + D(V) ~ 1 - V + end + (p/V,d/V), 0 <--> X + end + ``` + couples the ODE $dV/dt = 1 - V$ to the reaction system. +- Coupled reaction networks and differential equation (or algebraic differential + equation) systems can now be converted to `SDESystem`s and `NonlinearSystem`s. -#### Bifurcation analysis extension -- Add a CatalystBifurcationKitExtension, permitting BifurcationKit's `BifurcationProblem`s to be created from Catalyst reaction networks. Example usage: -```julia -using Catalyst -wilhelm_2009_model = @reaction_network begin - k1, Y --> 2X - k2, 2X --> X + Y - k3, X + Y --> Y - k4, X --> 0 - k5, 0 --> X -end +#### Structural identifiability extension +- Added CatalystStructuralIdentifiabilityExtension, which permits + StructuralIdentifiability.jl to be applied directly to Catalyst systems. E.g. + use + ```julia + using Catalyst, StructuralIdentifiability + goodwind_oscillator = @reaction_network begin + (mmr(P,pₘ,1), dₘ), 0 <--> M + (pₑ*M,dₑ), 0 <--> E + (pₚ*E,dₚ), 0 <--> P + end + assess_identifiability(goodwind_oscillator; measured_quantities=[:M]) + ``` + to assess (global) structural identifiability for all parameters and variables + of the `goodwind_oscillator` model (under the presumption that we can measure + `M` only). +- Automatically handles conservation laws for structural identifiability + problems (eliminates these internally to speed up computations). +- A more detailed of how this extension works can be found + [here](https://docs.sciml.ai/Catalyst/stable/inverse_problems/structural_identifiability/). + +#### Bifurcation analysis extension +- Add a CatalystBifurcationKitExtension, permitting BifurcationKit's + `BifurcationProblem`s to be created from Catalyst reaction networks. Example + usage: + ```julia + using Catalyst + wilhelm_2009_model = @reaction_network begin + k1, Y --> 2X + k2, 2X --> X + Y + k3, X + Y --> Y + k4, X --> 0 + k5, 0 --> X + end -using BifurcationKit -bif_par = :k1 -u_guess = [:X => 5.0, :Y => 2.0] -p_start = [:k1 => 4.0, :k2 => 1.0, :k3 => 1.0, :k4 => 1.5, :k5 => 1.25] -plot_var = :X -bprob = BifurcationProblem(wilhelm_2009_model, u_guess, p_start, bif_par; plot_var = plot_var) + using BifurcationKit + bif_par = :k1 + u_guess = [:X => 5.0, :Y => 2.0] + p_start = [:k1 => 4.0, :k2 => 1.0, :k3 => 1.0, :k4 => 1.5, :k5 => 1.25] + plot_var = :X + bprob = BifurcationProblem(wilhelm_2009_model, u_guess, p_start, bif_par; plot_var = plot_var) -p_span = (2.0, 20.0) -opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000) + p_span = (2.0, 20.0) + opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000) -bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true) + bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true) -using Plots -plot(bif_dia; xguide = "k1", guide = "X") -``` -- Automatically handles elimination of conservation laws for computing bifurcation diagrams. + using Plots + plot(bif_dia; xguide = "k1", guide = "X") + ``` +- Automatically handles elimination of conservation laws for computing + bifurcation diagrams. - Updated Bifurcation documentation with respect to this new feature. ## Catalyst 13.5 diff --git a/docs/src/api.md b/docs/src/api.md index 8a5e543973..459dbd1c1b 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -77,6 +77,7 @@ plot(p1, p2, p3; layout = (3,1)) ```@docs @reaction_network +@network_component make_empty_network @reaction Reaction diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md index ed47a2a439..a1e41b277a 100644 --- a/docs/src/v14_migration_guide.md +++ b/docs/src/v14_migration_guide.md @@ -16,7 +16,7 @@ A model's completeness depends on how it was created: - To *use the DSL to create models that are not marked as complete*, use the `@network_component` macro (which in all other aspects is identical to `@reaction_network`). - Models generated through the `compose` and `extend` functions are *not marked as complete*. -Furthermore, any systems generated through e.g. `convert(ODESystem, rs)` are *not marked as complete*. +Furthermore, any systems generated through e.g. `convert(ODESystem, rs)` are *not marked as complete*. Complete models can be generated from incomplete models through the `complete` function. Here is a workflow where we take completeness into account in the simulation of a simple birth-death process. ```@example v14_migration_1 @@ -64,8 +64,30 @@ sol = solve(oprob) plot(sol) ``` +Note, that if we had instead used the [`@reaction_network`](@ref) DSL macro to build +our model, i.e. +```@example v14_migration_1 +rs2 = @reaction_network rs begin + p, ∅ --> X + d, X --> ∅ +end +``` +then the model is automatically marked as complete +```@example v14_migration_1 +Catalyst.iscomplete(rs2) +``` +In contrast, if we used the [`@network_component`](@ref) DSL macro to build our +model it is not marked as complete, and is equivalent to our original definition of `rs` +```@example v14_migration_1 +rs3 = @network_component rs begin + p, ∅ --> X + d, X --> ∅ +end +Catalyst.iscomplete(rs3) +``` + ## Unknowns instead of states -Previously, "states" was used as a term for system variables (both species and non-species variables). MTKv9 has switched to using the term "unknowns" instead. This means that there have been a number of changes to function names (e.g. `states` => `unknowns` and `get_states` => `get_unknowns`). +Previously, "states" was used as a term for system variables (both species and non-species variables). MTKv9 has switched to using the term "unknowns" instead. This means that there have been a number of changes to function names (e.g. `states` => `unknowns` and `get_states` => `get_unknowns`). E.g. here we declare a `ReactionSystem` model containing both species and non-species unknowns: ```@example v14_migration_2 @@ -100,10 +122,10 @@ As part of its v9 update, ModelingToolkit changed how units were handled. This i While this should lead to long-term improvements, unfortunately, as part of the process support for most units was removed. Currently, only the main SI units are supported (`s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`). Composite units (e.g. `N = kg/(m^2)`) are no longer supported. Furthermore, prefix units (e.g. `mm = m/1000`) are not supported either. This means that most units relevant to Catalyst (such as `µM`) cannot be used directly. While composite units can still be written out in full and used (e.g. `kg/(m^2)`) this is hardly user-friendly. -The maintainers of ModelingTOolkit have been notified of this issue. We are unsure when this will be fixed, however, we do not think it will be a permanent change. +The maintainers of ModelingToolkit have been notified of this issue. We are unsure when this will be fixed, however, we do not think it will be a permanent change. ## Removed support for system-mutating functions -According to the ModelingToolkit system API, systems should not be mutable. In accordance with this, the following functions have been deprecated: `addparam!`, `addreaction!`, `addspecies!`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystems` from multiple component systems. +According to the ModelingToolkit system API, systems should not be mutable. In accordance with this, the following functions have been deprecated and removed: `addparam!`, `addreaction!`, `addspecies!`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystems` from multiple component systems. It is still possible to add default values to a created `ReactionSystem`, i.e. the `setdefaults!` function is still supported. @@ -132,9 +154,9 @@ nothing # hide ``` !!! note - If you look at ModelingToolkit documentation, these defaults are instead retrieved using `using ModelingToolkit: t_nounits as t, D_nounits as D`. This will also work, however, in Catalyst we have opted to instead use `default_t` and `default_time_deriv` as our main approach. + If you look at ModelingToolkit documentation, these defaults are instead retrieved using `using ModelingToolkit: t_nounits as t, D_nounits as D`. This will also work, however, in Catalyst we have opted to instead use the functions `default_t()` and `default_time_deriv()` as our main approach. -## New interface for accessing problem/integrator/solution parameter (and species) value +## New interface for accessing problem/integrator/solution parameter (and species) values Previously, it was possible to directly index problems to query them for their parameter values. e.g. ```@example v14_migration_4 using Catalyst @@ -154,12 +176,11 @@ This is *no longer supported*. When you wish to query a problem (or integrator o oprob.ps[:p] ``` -Furthermore, a few new functions (`getp`, `getu`, `setp`, `setu`) have been introduced. These can improve performance when querying a structure for a value multiple times (especially for very large models). These are described in more detail [here](@ref simulation_structure_interfacing_functions). +Furthermore, a few new functions (`getp`, `getu`, `setp`, `setu`) have been introduced from [SymbolicIndexingInterface](https://github.com/SciML/SymbolicIndexingInterface.jl) to support efficient and systematic querying and/or updating of symbolic unknown/parameter values. Using these can *significantly* improve performance when querying or updating a value multiple times, for example within a callback. These are described in more detail [here](@ref simulation_structure_interfacing_functions). For more details on how to query various structures for parameter and species values, please read [this documentation page](@ref simulation_structure_interfacing). -## Other notes -Finally, here are some additional, minor, notes regarding the new update. +## Other changes #### Modification of problems with conservation laws broken While it is possible to update e.g. `ODEProblem`s using the [`remake`](@ref simulation_structure_interfacing_problems_remake) function, this is currently not possible if the `remove_conserved = true` option was used. E.g. while @@ -174,7 +195,7 @@ oprob = ODEProblem(rn, u0, (0.0, 10.0), ps; remove_conserved = true) solve(oprob) # hide ``` -is perfectly fine, attempting to then modify any initial conditions or the value of the conservation constant in `oprob` will silently fail: +is perfectly fine, attempting to then modify any initial conditions or the value of the conservation constant in `oprob` will likely silently fail: ```@example v14_migration_5 oprob_remade = remake(oprob; u0 = [:X1 => 5.0]) # NEVER do this. solve(oprob_remade) diff --git a/src/dsl.jl b/src/dsl.jl index 633ea7f64b..6148eb4814 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -175,7 +175,8 @@ end """ @network_component -As @reaction_network, but the output system is not complete. +Equivalent to `@reaction_network` except the generated `ReactionSystem` is not marked as +complete. """ macro network_component(name::Symbol, ex::Expr) make_reaction_system(MacroTools.striplines(ex); name = :($(QuoteNode(name)))) From 91e2bc22f6b96bd76437dce45402ce4c7e92bfab Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Thu, 11 Jul 2024 11:23:20 -0400 Subject: [PATCH 2/2] Update docs/src/v14_migration_guide.md --- docs/src/v14_migration_guide.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md index a1e41b277a..ed5074839b 100644 --- a/docs/src/v14_migration_guide.md +++ b/docs/src/v14_migration_guide.md @@ -64,7 +64,7 @@ sol = solve(oprob) plot(sol) ``` -Note, that if we had instead used the [`@reaction_network`](@ref) DSL macro to build +Note, if we had instead used the [`@reaction_network`](@ref) DSL macro to build our model, i.e. ```@example v14_migration_1 rs2 = @reaction_network rs begin