mereg

Merge branch 'master' of github.com:SciML/Catalyst.jl into spatial_SSA_support

.github/workflows/CI.yml

@@ -17,7 +17,7 @@ jobs:
           - '1'
           - '1.6'
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}

.github/workflows/Documentation.yml

@@ -11,7 +11,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'

.github/workflows/FormatCheck.yml

@@ -21,7 +21,7 @@ jobs:
         with:
           version: ${{ matrix.julia-version }}
 
-      - uses: actions/checkout@v3
+      - 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:
         #

.github/workflows/Invalidations.yml

@@ -19,12 +19,12 @@ jobs:
     - uses: julia-actions/setup-julia@v1
       with:
         version: '1'
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - uses: julia-actions/julia-buildpkg@v1
     - uses: julia-actions/julia-invalidations@v1
       id: invs_pr
 
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         ref: ${{ github.event.repository.default_branch }}
     - uses: julia-actions/julia-buildpkg@v1

HISTORY.md

@@ -14,6 +14,55 @@
 
 
 
+## Catalyst 13.4
+- Added the ability to create species that represent chemical compounds and know
+  their constituents. For example, water can be created and queried as
+  ```julia
+  @variables t
+  @species H(t) O(t)
+  @compound H2O(t) 2*H O
+  iscompound(H2O) == true
+  isspecies(H2O) == true   # compounds are also species, so can be used in reactions
+  isequal(components(H2O), [H, O])
+  coefficients(H2O) == [2, 1]
+  ```
+- Added reaction balancing via the `balance_reaction` command, which returns a
+  vector of balanced reaction versions, i.e.
+  ```julia
+  @variables t
+  @species H(t) O(t) C(t)
+  @compound CH4(t) C 4H
+  @compound O2(t) 2O
+  @compound CO2(t) C 2O
+  @compound H2O(t) 2H O
+
+  # unbalanced reaction to balance
+  rx = Reaction(1.0, [CH4, O2], [CO2, H2O])
+
+  # calculate a balanced version, this returns a vector
+  # storing a single balanced version of the reaction in this case
+  brxs = balance_reaction(rx)
+
+  # what one would calculate by hand
+  balanced_rx = Reaction(1.0, [CH4, O2], [CO2, H2O], [1, 2], [1, 2])
+
+  # testing equality
+  @test isequal(balanced_rx, first(brxs))
+  ```
+- Note that balancing works via calculating the nullspace of an associated
+  integer matrix that stores in entry `(i,j)` a signed integer representing the
+  number of times the `i`'th atom appears within the `j`th compound. The entry
+  is positive for a substrate and negative for a product. One cannot balance a
+  reaction involving compounds of compounds currently. A non-empty solution
+  vector is returned if the reaction can be balanced in exactly one way with
+  minimal coefficients while preserving the set of substrates and products, i.e.
+  if the dimension of the nullspace is one. If the dimension is greater than one
+  we return a `Reaction` for each nullspace basis vector, but note that they may
+  currently interchange substrates and products (i.e. we do not solve for if
+  there is a linear combination of them that preserves the set of substrates and
+  products). An empty `Reaction` vector indicates it is not possible to balance
+  the reaction.
+
 ## Catalyst 13.2
 - Array parameters, species, and variables can be use in the DSL if explicitly
   declared with `@parameters`, `@species`, or `@variables` respectively, i.e.

Project.toml

@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "13.3.0"
+version = "13.4.1"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -30,7 +30,7 @@ JumpProcesses = "9.3.2"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"
-ModelingToolkit = "8.47"
+ModelingToolkit = "8.66"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"

docs/Project.toml

@@ -28,7 +28,7 @@ HomotopyContinuation = "2.6"
 Latexify = "0.15, 0.16"
 ModelingToolkit = "8.47"
 NonlinearSolve = "1.6.0"
-Optimization = "3.11"
+Optimization = "~3.15"
 OptimizationOptimisers = "0.1.1"
 OrdinaryDiffEq = "6"
 Plots = "1.36"

docs/src/catalyst_applications/simulation_structure_interfacing.md

@@ -1,12 +1,12 @@
 # [Interfacing problems, integrators, and solutions](@id simulation_structure_interfacing)
-When simulating a model, one begins with creating a [problem](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/). Next, a simulation is performed on a problem, during which the state of the simulation is recorded through an [integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/). Finally, the simulation output is returned as a [solution](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/). This tutorial describes how to access, or modify the state, or parameter, values of problems, integrators, and solutions structures. 
+When simulating a model, one begins with creating a [problem](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/). Next, a simulation is performed on a problem, during which the state of the simulation is recorded through an [integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/). Finally, the simulation output is returned as a [solution](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/). This tutorial describes how to access, or modify the state, or parameter, values of problems, integrators, and solutions structures.
 
 Generally, when we have a structure `simulation_struct` and want to interface with the state (or parameter) `G`, we use `simulation_struct[:G]` to access the value, and `simulation_struct[:G] = 5.0` to set it to a new value. However, see the following examples for full details.
 
 ## Interfacing problem objects
 
 We begin by demonstrating how we can interface with problem objects. We will demonstrate using a `ODEProblem`, however, it works similarily for other problem types.
-```example ex1
+```@example ex1
 using Catalyst
 rn = @reaction_network begin
     (k1,k2), X1 <--> X2
@@ -19,43 +19,43 @@ nothing    # hide
 ```
 
 We can find the value of a state simply by interfacing with the corresponding symbol:
-```example ex1
+```@example ex1
 oprob[:X1]
 ```
 with the notation being identical for parameters:
-```example ex1
+```@example ex1
 oprob[:k1]
 ```
 
 If we want to change a state's initial condition value, we use the following notation
-```example ex1
+```@example ex1
 oprob[:X1] = 10.0
 ```
 with parameters using the same notation.
 
 #### Remaking problems using the `remake` function
 Typically, when modifying problems, it is recommended to use the `remake` function. Unlike when we do `oprob[:X1] = 10.0` (which modifies the problem in question), `remake` creates a new problem object. The `remake` function takes a problem as input, and any fields you wish to modify (and their new values) as optional inputs. Thus, we can do:
-```example ex1
+```@example ex1
 using DifferentialEquations
 @unpack X1, X2, k1, k2 = rn
 oprob1 = ODEProblem(rn, u0, (0.0,10.0), p)
 oprob2 = remake(oprob1; u0=[X1 => 10.0, X2 => 50.0], tspan=(0.0,100.0), p=[k1 => 50.0,k2 => 20.0])
 nothing    # hide
 ```
 and we can now check the fields of `oprob2`
-```example ex1
+```@example ex1
 oprob2.u0
 ```
-```example ex1
+```@example ex1
 oprob2.tspan
 ```
-```example ex1
+```@example ex1
 oprob2.p
 ```
 Please note that, currently, `remake` does not work while giving `Symbol`s as input (e.g `[:X1 => 10.0, :X2 => 50.0]`), but we need to unpack the symbolic variables and use them instead (please see the end of this tutorial for more information on using symbolic variables rather than `Symbol`s).
 
 When using `remake`, we only have to provide the fields that we actually wish to change, e.g.
-```example ex1
+```@example ex1
 oprob3 = remake(oprob1; u0=[X1 => 10.0, X2 => 50.0])
 nothing    # hide
 ```
@@ -65,19 +65,19 @@ will only update the initial conditions.
 ## Interfacing integrator objects
 
 During a simulation, the solution is stored in an integrator object, we will here describe how to interface with these. The primary circumstance under which a user may wish to do so is when using [callbacks](@ref advanced_simulations_callbacks). We can create an integrator by calling `init` on our problem ([while circumstances where the user might want to use `init` function exist](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Initialization-and-Stepping), since integrators are automatically created during simulations, these are rare):
-```example ex1
+```@example ex1
 integrator = init(oprob)
 ```
 Using a similar syntax to problems, we can get the current values of a state within the integrator:
-```example ex1
+```@example ex1
 integrator[:X1]
 ```
 or a parameter:
-```example ex1
+```@example ex1
 integrator[:k1]
 ```
 Similarly, we can update their values using:
-```example ex1
+```@example ex1
 integrator[:X1] = 10.0
 ```
 Please read [this](@ref advanced_simulations_ssa_callbacks) with regards to updating integrators of `JumpProblem`s.
@@ -86,30 +86,30 @@ Please read [this](@ref advanced_simulations_ssa_callbacks) with regards to upda
 ## Interfacing solution objects
 
 Finally, we consider solution objects. First, we simulate our problem:
-```example ex1
+```@example ex1
 sol = solve(oprob)
 ```
 For solutions, when we access a state, we get its whole simulation vector:
-```example ex1
+```@example ex1
 sol[:X1]
 ```
 while when we access a parameter we only get a single value:
-```example ex1
+```@example ex1
 sol[:k1]
 ```
 Finally, we note that we cannot change the values of solution states or parameters (i.e. both `sol[:X1] = 0.0` and `sol[:k1] = 0.0` generate errors).
 
 ## Interfacing using symbolic representation
 
 Catalyst is built on an *intermediary representation* implemented by (ModelingToolkit.jl)[https://github.com/SciML/ModelingToolkit.jl]. ModelingToolkit is a modelling framework where one first declares a set of symbolic variables and parameters using e.g.
-```example ex2
+```@example ex2
 using ModelingToolkit
 @parameters σ ρ β
 @variables t x(t) y(t) z(t)
 nothing # hide
 ```
-and then uses these to build systems of equations. Here, these symbolic variables (`x`, `y`, and `z`) and parameters (`σ`, `ρ`, and `β`) can be used to interface a `problem`, `integrator`, and `solution` object (like we did previously, but using Symbols, e.g. `:X`). Since Catalyst models are built on ModelingToolkit, these models also contain similar symbolic variables and parameters. 
-```example ex2
+and then uses these to build systems of equations. Here, these symbolic variables (`x`, `y`, and `z`) and parameters (`σ`, `ρ`, and `β`) can be used to interface a `problem`, `integrator`, and `solution` object (like we did previously, but using Symbols, e.g. `:X`). Since Catalyst models are built on ModelingToolkit, these models also contain similar symbolic variables and parameters.
+```@example ex2
 using Catalyst
 rn = @reaction_network begin
     (k1,k2), X1 <--> X2
@@ -118,7 +118,7 @@ end
 @unpack k1,k2,X1,X2 = rn
 ```
 Here, we first list the parameters and variables (for reaction systems the latter are typically species) we wish to import (in this case we select all, but we could select only a subset), next we denote from which model (here `rn`) from which we wish to import from. Next, these values can be used directly to interface with e.g. an `ODEProblem`:
-```example ex2
+```@example ex2
 u0 = [X1 => 1.0, X2 => 5.0]
 p = [:k1 => 5.0, :k2 => 2.0]
 oprob = ODEProblem(rn, u0, (0.0,10.0), p)
@@ -128,6 +128,6 @@ oprob[k1]
 To interface with integrators and solutions we use a similar syntax.
 
 Finally, instead of using `@unpack` to access a symbolic variable or parameter, we can access it directly using `rn.X1`, and thus access a state of our `ODEProblem` using
-```example ex2
+```@example ex2
 oprob[rn.X1]
 ```

docs/src/catalyst_functionality/dsl_description.md

@@ -412,6 +412,26 @@ sol = solve(oprob)
 plot(sol)
 ```
 
+## Setting initial conditions that depend on parameters
+It is possible to set the initial condition of one (or several) species so that they depend on some system parameter. This is done in a similar way as default initial conditions, but giving the parameter instead of a value. When doing this, we also need to ensure that the initial condition parameter is a variable of the system:
+```@example tut2
+rn = @reaction_network begin
+  @parameters X0
+  @species X(t)=X0
+  p, 0 --> X
+  d, X --> ∅
+end
+```
+We can now simulate the network without providing any initial conditions:
+```@example tut2
+u0 = []
+tspan = (0.0, 10.0)
+p = [:p => 2.0, :d => .1, :X0 => 1.0]  
+oprob = ODEProblem(rn, u0, tspan, p)
+sol = solve(oprob)
+plot(sol)
+```
+
 ## Naming the generated `ReactionSystem`
 ModelingToolkit uses system names to allow for compositional and hierarchical
 models. To specify a name for the generated `ReactionSystem` via the

docs/src/example_networks/hodgkin_huxley_equation.md

@@ -14,7 +14,7 @@ We begin by importing some necessary packages.
 ```@example hh1
 using ModelingToolkit, Catalyst, NonlinearSolve
 using DifferentialEquations, Symbolics
-using Plots, IfElse
+using Plots
 ```
 
 We'll build a simple Hodgkin-Huxley model for a single neuron, with the voltage,
@@ -32,7 +32,7 @@ The transition rate functions, which depend on the voltage, $V(t)$, are given by
 ```@example hh1
 function αₘ(V)
     theta = (V + 45) / 10
-    IfElse.ifelse(theta == 0.0, 1.0, theta/(1 - exp(-theta)))
+    ifelse(theta == 0.0, 1.0, theta/(1 - exp(-theta)))
 end
 βₘ(V) = 4*exp(-(V + 70)/18)
 
@@ -41,14 +41,14 @@ end
 
 function αₙ(V)
     theta = (V + 60) / 10
-    IfElse.ifelse(theta == 0.0, .1, .1*theta / (1 - exp(-theta)))
+    ifelse(theta == 0.0, .1, .1*theta / (1 - exp(-theta)))
 end
 βₙ(V) = .125 * exp(-(V + 70)/80)
 nothing # hide
 ```
 
- We now declare the symbolic variable, `V(t)`, that will represent the
- transmembrane potential
+We now declare the symbolic variable, `V(t)`, that will represent the
+transmembrane potential
 
 ```@example hh1
 @variables t V(t)

docs/src/index.md

@@ -4,7 +4,7 @@ Catalyst.jl is a symbolic modeling package for analysis and high performance
 simulation of chemical reaction networks. Catalyst defines symbolic
 [`ReactionSystem`](@ref)s, which can be created programmatically or easily
 specified using Catalyst's domain specific language (DSL). Leveraging
-[ModelingToolkit](https://docs.sciml.ai/ModelingToolkit/stable/) and
+[ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) and
 [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/), Catalyst enables
 large-scale simulations through auto-vectorization and parallelism. Symbolic
 `ReactionSystem`s can be used to generate ModelingToolkit-based models, allowing

docs/src/introduction_to_catalyst/introduction_to_catalyst.md

@@ -194,10 +194,10 @@ We see that oscillations remain, but become much noisier. Note, in constructing
 the `JumpProblem` we could have used any of the SSAs that are part of JumpProcesses
 instead of the `Direct` method, see the list of SSAs (i.e., constant rate jump
 aggregators) in the
-[documentation](https://docs.sciml.ai/stable/modules/JumpProcesses/jump_types/#Constant-Rate-Jump-Aggregators).
+[documentation](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-for-Exact-Simulation).
 
 Common questions that arise in using the JumpProcesses SSAs (i.e. Gillespie methods)
-are collated in the [JumpProcesses FAQ](https://docs.sciml.ai/stable/modules/JumpProcesses/faq/).
+are collated in the [JumpProcesses FAQ](https://docs.sciml.ai/JumpProcesses/stable/faq/).
 
 ---
 ## Chemical Langevin equation (CLE) stochastic differential equation (SDE) models

src/Catalyst.jl

@@ -6,7 +6,8 @@ module Catalyst
 using DocStringExtensions
 using SparseArrays, DiffEqBase, Reexport
 using LaTeXStrings, Latexify, Requires
-using JumpProcesses: JumpProcesses, JumpProblem, MassActionJump, ConstantRateJump,
+using JumpProcesses: JumpProcesses,
+                     JumpProblem, MassActionJump, ConstantRateJump,
                      VariableRateJump
 
 # ModelingToolkit imports and convenience functions we use
@@ -26,7 +27,8 @@ import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_v
                         namespace_parameters, rename, renamespace, getname, flatten
 
 # internal but needed ModelingToolkit functions
-import ModelingToolkit: check_variables, check_parameters, _iszero, _merge, check_units,
+import ModelingToolkit: check_variables,
+                        check_parameters, _iszero, _merge, check_units,
                         get_unit, check_equations
 
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
@@ -56,7 +58,8 @@ end
 include("reactionsystem.jl")
 export isspecies
 export Reaction, ReactionSystem, ismassaction, oderatelaw, jumpratelaw, isspatial
-export ODEProblem, SDEProblem, JumpProblem, NonlinearProblem, DiscreteProblem,
+export ODEProblem,
+       SDEProblem, JumpProblem, NonlinearProblem, DiscreteProblem,
        SteadyStateProblem
 
 # reaction_network macro
@@ -94,9 +97,10 @@ include("graphs.jl")
 export Graph, savegraph, complexgraph
 
 # for creating compounds
-include("compound.jl")
+include("chemistry_functionality.jl")
 export @compound
 export components, iscompound, coefficients
+export balance_reaction
 
 ### Spatial Reaction Networks ###