Merge remote-tracking branch 'origin/master' into linkageclasses

SciML · Jul 12, 2024 · c207d98 · c207d98
2 parents 6627760 + 68c0664
commit c207d98
Show file tree

Hide file tree

Showing 60 changed files with 4,214 additions and 2,045 deletions.
diff --git a/HISTORY.md b/HISTORY.md
@@ -3,103 +3,190 @@
 ## Catalyst unreleased (master branch)
 
 ## Catalyst 14.0
-- The `reactionparams`, `numreactionparams`, and `reactionparamsmap`  functions have been 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.
-- 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).
-- Adds a tutorial to illustrate the use of the extension.
-- Enable adding metadata to individual reactions, e.g:
-```julia
-rn = @reaction_network begin
-    @parameters η
-    k, 2X --> X2, [noise_scaling=η]
-end
-getnoisescaling(rn)
-```
-- `SDEProblem` no longer takes the `noise_scaling` argument (see above for new approach to handle noise scaling).
-- Changed fields of internal `Reaction` structure. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions.
-- Simulation of spatial ODEs now supported. For full details, please see https://github.com/SciML/Catalyst.jl/pull/644 and upcoming documentation. Note that these methods are currently considered alpha, with the interface and approach changing even in non-breaking Catalyst releases.
-- LatticeReactionSystem structure represents a spatial reaction network:
+
+#### 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.
+  ```julia
+  # do
+  t = default_t()
+  @species A(t)
+
+  # avoid
+  @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)
+- 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
-  tr = @transport_reaction D X
-  lattice = Graphs.grid([5, 5])
-  lrs = LatticeReactionSystem(rn, [tr], lattice)
-```
-- Here, if a `u0` or `p` vector is given with scalar values:
+  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 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
-  u0 = [:X => 1.0]
-  p = [:p => 1.0, :d => 0.5, :D => 0.1]
+  rn = @reaction_network begin
+      @equations begin
+          D(V) ~ 1 - V
+      end
+      (p/V,d/V), 0 <--> X
+  end
   ```
-  this value will be used across the entire system. If their values are instead vectors, different values are used across the spatial system. Here
+  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.
+
+#### Structural identifiability extension
+- Added CatalystStructuralIdentifiabilityExtension, which permits
+  StructuralIdentifiability.jl to be applied directly to Catalyst systems. E.g.
+  use
   ```julia
-  X0 = zeros(25)
-  X0[1] = 1.0
-  u0 = [:X => X0]
+  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])
   ```
-  X's value will be `1.0` in the first vertex, but `0.0` in the remaining one (the system have 25 vertexes in total). SInce th parameters `p` and `d` are part of the non-spatial reaction network, their values are tied to vertexes. However, if the `D` parameter (which governs diffusion between vertexes) is given several values, these will instead correspond to the specific edges (and transportation along those edges.)
+  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/).
 
-- 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).
-- 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
+#### 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", yguide="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.
-- 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]
-
-# Finds (the trivial) steady state, and computes stability.
-steady_state = [2.0]
-steady_state_stability(steady_state, rn, p)
-```
-Here, `steady_state_stability` take an optional argument `tol = 10*sqrt(eps())`, which is used to determine whether a eigenvalue real part is reliably less that 0.
 
 ## Catalyst 13.5
 - Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example:
@@ -658,7 +745,7 @@ hc_steady_states(wilhelm_2009_model, ps)
   field has been changed (only when created through the `@reaction_network`
   macro). Previously they were ordered according to the order with which they
   appeared in the macro. Now they are ordered according the to order with which
-  they appeard after the `end` part. E.g. in
+  they appeared after the `end` part. E.g. in
   ```julia
   rn = @reaction_network begin
     (p,d), 0 <--> X
@@ -763,7 +850,7 @@ which gives
 ![rn_complexes](https://user-images.githubusercontent.com/9385167/130252763-4418ba5a-164f-47f7-b512-a768e4f73834.png)
 
 *2.* Support for units via ModelingToolkit and
-[Uniftul.jl](https://github.com/PainterQubits/Unitful.jl) in directly constructed
+[Unitful.jl](https://github.com/PainterQubits/Unitful.jl) in directly constructed
 `ReactionSystem`s:
 ```julia
 # ]add Unitful

diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "14.0.0-DEV"
+version = "14.0.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -76,6 +76,7 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLNLSolve = "e9a6253c-8580-4d32-9898-8661bb511710"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
@@ -84,4 +85,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["BifurcationKit", "DiffEqCallbacks", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "Logging", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"]
+test = ["BifurcationKit", "DiffEqCallbacks", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "Logging", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "StaticArrays", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"]