Merge branch 'master' into spatial___documentation

.github/workflows/CI.yml

@@ -17,10 +17,10 @@ jobs:
           - '1'
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
-      - uses: actions/cache@v3
+      - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:
@@ -33,9 +33,11 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
           file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: false
       - uses: coverallsapp/github-action@master
         with:
           github-token: ${{ secrets.GITHUB_TOKEN }}

.github/workflows/CompatHelper.yml

@@ -15,7 +15,7 @@ jobs:
         run: which julia
         continue-on-error: true
       - name: Install Julia, but only if it is not already available in the PATH
-        uses: julia-actions/setup-julia@v1
+        uses: julia-actions/setup-julia@v2
         with:
           version: '1'
           arch: ${{ runner.arch }}

.github/workflows/Documentation.yml

@@ -24,7 +24,8 @@ jobs:
           GKSwstype: "100" # https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988
         run: julia --project=docs/ --code-coverage=user docs/make.jl
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
           file: lcov.info
-
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: false

.github/workflows/FormatCheck.yml

@@ -1,4 +1,4 @@
-name: format-check
+name: "Format Check"
 
 on:
   push:
@@ -9,34 +9,6 @@ on:
   pull_request:
 
 jobs:
-  build:
-    runs-on: ${{ matrix.os }}
-    strategy:
-      matrix:
-        julia-version: [1]
-        julia-arch: [x86]
-        os: [ubuntu-latest]
-    steps:
-      - uses: julia-actions/setup-julia@latest
-        with:
-          version: ${{ matrix.julia-version }}
-
-      - uses: actions/checkout@v4
-      - name: Install JuliaFormatter and format
-        # This will use the latest version by default but you can set the version like so:
-        #
-        # julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))'
-        run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="1.0.32"))'
-          julia  -e 'using JuliaFormatter; format(".", verbose=true)'
-      - name: Format check
-        run: |
-          julia -e '
-          out = Cmd(`git diff --name-only`) |> read |> String
-          if out == ""
-              exit(0)
-          else
-              @error "Some files have not been formatted !!!"
-              write(stdout, out)
-              exit(1)
-          end'
+  format-check:
+    name: "Format Check"
+    uses: "SciML/.github/.github/workflows/format-check.yml@v1"

.github/workflows/Invalidations.yml

@@ -4,37 +4,10 @@ on:
   pull_request:
 
 concurrency:
-  # Skip intermediate builds: always.
-  # Cancel intermediate builds: always.
   group: ${{ github.workflow }}-${{ github.ref }}
   cancel-in-progress: true
 
 jobs:
-  evaluate:
-    # Only run on PRs to the default branch.
-    # In the PR trigger above branches can be specified only explicitly whereas this check should work for master, main, or any other default branch
-    if: github.base_ref == github.event.repository.default_branch
-    runs-on: ubuntu-latest
-    steps:
-    - uses: julia-actions/setup-julia@v1
-      with:
-        version: '1'
-    - uses: actions/checkout@v4
-    - uses: julia-actions/julia-buildpkg@v1
-    - uses: julia-actions/julia-invalidations@v1
-      id: invs_pr
-
-    - uses: actions/checkout@v4
-      with:
-        ref: ${{ github.event.repository.default_branch }}
-    - uses: julia-actions/julia-buildpkg@v1
-    - uses: julia-actions/julia-invalidations@v1
-      id: invs_default
-
-    - name: Report invalidation counts
-      run: |
-        echo "Invalidations on default branch: ${{ steps.invs_default.outputs.total }} (${{ steps.invs_default.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
-        echo "This branch: ${{ steps.invs_pr.outputs.total }} (${{ steps.invs_pr.outputs.deps }} via deps)" >> $GITHUB_STEP_SUMMARY
-    - name: Check if the PR does increase number of invalidations
-      if: steps.invs_pr.outputs.total > steps.invs_default.outputs.total
-      run: exit 1
+  evaluate-invalidations:
+    name: "Evaluate Invalidations"
+    uses: "SciML/.github/.github/workflows/invalidations.yml@v1"

HISTORY.md

@@ -1,29 +1,107 @@
 # Breaking updates and feature summaries across releases
 
 ## Catalyst unreleased (master branch)
-- 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:
+
+## 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.
+
+#### 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
-  rn = @reaction_network begin
-      (p,d), 0 <--> X
-  end
-  tr = @transport_reaction D X
-  lattice = Graphs.grid([5, 5])
-  lrs = LatticeReactionSystem(rn, [tr], lattice)
+  # 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)
 ```
-- Here, if a `u0` or `p` vector is given with scalar values:
-  ```julia
-  u0 = [:X => 1.0]
-  p = [:p => 1.0, :d => 0.5, :D => 0.1]
-  ```
-  this value will be used across the entire system. If their values are instead vectors, different values are used across the spatial system. Here
-  ```julia
-  X0 = zeros(25)
-  X0[1] = 1.0
-  u0 = [:X => X0]
-  ```
-  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.)
+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).
+- 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` 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/).
+
+#### Bifurcation analysis extension 
 - Add a CatalystBifurcationKitExtension, permitting BifurcationKit's `BifurcationProblem`s to be created from Catalyst reaction networks. Example usage:
 ```julia
 using Catalyst
@@ -35,21 +113,20 @@ wilhelm_2009_model = @reaction_network begin
     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)
+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)
+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")
+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.
@@ -611,7 +688,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
@@ -716,7 +793,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

Project.toml

@@ -1,11 +1,14 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "13.5.1"
+version = "14.0.0-DEV"
 
 [deps]
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
@@ -19,58 +22,67 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 
 [extensions]
 CatalystBifurcationKitExtension = "BifurcationKit"
 CatalystHomotopyContinuationExtension = "HomotopyContinuation"
+CatalystStructuralIdentifiabilityExtension = "StructuralIdentifiability"
 
 [compat]
 BifurcationKit = "0.3"
+Combinatorics = "1.0.2"
 DataStructures = "0.18"
 DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
+DynamicPolynomials = "0.5"
+DynamicQuantities = "0.13.2"
 Graphs = "1.4"
-JumpProcesses = "9.3.2"
 HomotopyContinuation = "2.9"
+JumpProcesses = "9.3.2"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"
-ModelingToolkit = "8.73"
+ModelingToolkit = "9.16.0"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
 Setfield = "1"
-SymbolicUtils = "1.0.3"
-Symbolics = "5.0.3"
+StructuralIdentifiability = "0.5.8"
+Symbolics = "5.30.1"
 Unitful = "1.12.4"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 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"
+StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["BifurcationKit", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "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"]