Merge pull request #784 from SciML/reactionsystem_completeness

Handle test errors from MTKv9 update (including `ReactionSystem` completeness)

Author: TorkelE

Project.toml

@@ -44,7 +44,7 @@ HomotopyContinuation = "2.9"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"
-ModelingToolkit = "9.5"
+ModelingToolkit = "9.7.1"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"

docs/Project.toml

@@ -17,7 +17,6 @@ OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -47,7 +46,6 @@ OptimizationNLopt = "0.1.8"
 OptimizationOptimJL = "0.1.14"
 OptimizationOptimisers = "0.1.1"
 OrdinaryDiffEq = "6"
-PEtab = "2"
 Plots = "1.36"
 SciMLBase = "2.13"
 SciMLSensitivity = "7.19"

docs/pages.jl

@@ -1,6 +1,6 @@
 pages = Any["Home" => "index.md",
             "Introduction to Catalyst" => Any["introduction_to_catalyst/catalyst_for_new_julia_users.md",
-                                              "introduction_to_catalyst/introduction_to_catalyst.md" ],
+                                              "introduction_to_catalyst/introduction_to_catalyst.md"],
             "Catalyst Functionality" => Any["catalyst_functionality/dsl_description.md",
                                             "catalyst_functionality/programmatic_CRN_construction.md",
                                             "catalyst_functionality/compositional_modeling.md",
@@ -17,7 +17,7 @@ pages = Any["Home" => "index.md",
                                            "catalyst_applications/nonlinear_solve.md",
                                            "catalyst_applications/bifurcation_diagrams.md"],
             "Inverse Problems" => Any["inverse_problems/optimization_ode_param_fitting.md",
-                                        "inverse_problems/petab_ode_param_fitting.md",
+                                        #"inverse_problems/petab_ode_param_fitting.md",
                                         "inverse_problems/structural_identifiability.md",
                                         "Inverse problem examples" => Any["inverse_problems/examples/ode_fitting_oscillation.md"]],
             "FAQs" => "faqs.md",

docs/src/catalyst_functionality/compositional_modeling.md

@@ -5,16 +5,44 @@ can construct the earlier repressilator model by composing together three
 identically repressed genes, and how to use compositional modeling to create
 compartments.
 
+## [A note on *completeness*](@id completeness_note)
+Catalyst `ReactionSystem` can either be *complete* or *incomplete*. When created using the `@reaction_network` DSL they are *created as complete*. Here, only complete `ReactionSystem`s can be used to create the various problem types (e.g. `ODEProblem`). However, only incomplete `ReactionSystem`s can be composed using the features described below. Hence, for compositional modeling, `ReactionSystem` must be created as incomplete, and later set to complete before simulation.
+
+To create a `ReactionSystem`s for use in compositional modeling via the DSL, simply use the `@network_component` macro instead of `@reaction_network`:
+```@example ex0
+using Catalyst
+degradation_component = @network_component begin
+  d, X --> 0
+end
+```
+Alternatively one can just build the `ReactionSystem` via the symbolic interface.
+```@example ex0
+@parameters d
+@variable t
+@species X(t)
+rx = Reaction(d, [X], nothing)
+@named degradation_component = ReactionSystem([rs], t)
+```
+We can test whether a system is complete using the `ModelingToolkit.iscomplete` function:
+```@example ex0
+ModelingToolkit.iscomplete(degradation_component)
+```
+To mark a system as complete, after which is should be considered as representing a finalized model, use the `complete` function
+```@example ex0
+degradation_component_complete = complete(degradation_component)
+ModelingToolkit.iscomplete(degradation_component_complete)
+```
+
 ## Compositional modeling tooling
 Catalyst supports two ModelingToolkit interfaces for composing multiple
 [`ReactionSystem`](@ref)s together into a full model. The first mechanism for
 extending a system is the `extend` command
 ```@example ex1
 using Catalyst
-basern = @reaction_network rn1 begin
+basern = @network_component rn1 begin
   k, A + B --> C
 end
-newrn = @reaction_network rn2 begin
+newrn = @network_component rn2 begin
   r, C --> A + B
 end
 @named rn = extend(newrn, basern)
@@ -27,9 +55,9 @@ The second main compositional modeling tool is the use of subsystems. Suppose we
 now add to `basern` two subsystems, `newrn` and `newestrn`, we get a
 different result:
 ```@example ex1
-newestrn = @reaction_network rn3 begin
-            v, A + D --> 2D
-           end
+newestrn = @network_component rn3 begin
+  v, A + D --> 2D
+end
 @named rn = compose(basern, [newrn, newestrn])
 ```
 Here we have created a new `ReactionSystem` that adds `newrn` and `newestrn` as
@@ -99,7 +127,7 @@ modular fashion. We start by defining a function that creates a negatively
 repressed gene, taking the repressor as input
 ```@example ex1
 function repressed_gene(; R, name)
-  @reaction_network $name begin
+  @network_component $name begin
     hillr($R,α,K,n), ∅ --> m
     (δ,γ), m <--> ∅
     β, m --> m + P
@@ -156,13 +184,13 @@ In our model we'll therefore add the conversions of the last column to properly
 account for compartment volumes:
 ```@example ex1
 # transcription and regulation
-nuc = @reaction_network nuc begin
+nuc = @network_component nuc begin
   α, G --> G + M
   (κ₊/V,κ₋), D + G <--> DG
 end
 
 # translation and dimerization
-cyto = @reaction_network cyto begin
+cyto = @network_component cyto begin
   β, M --> M + P
   (k₊/V,k₋), 2P <--> D
   σ, P --> 0
@@ -171,7 +199,7 @@ end
 
 # export reactions,
 # γ,δ=probability per time to be exported/imported
-model = @reaction_network model begin
+model = @network_component model begin
   γ, $(nuc.M) --> $(cyto.M)
   δ, $(cyto.D) --> $(nuc.D)
 end

docs/src/catalyst_functionality/programmatic_CRN_construction.md

@@ -61,6 +61,9 @@ system to be the same as the name of the variable storing the system.
 Alternatively, one can use the `name = :repressilator` keyword argument to the
 `ReactionSystem` constructor.
 
+!!! warn
+       All `ReactionSystem`s created via the symbolic interface (i.e. by calling `ReactionSystem` with some input, rather than using `@reaction_network`) are not marked as complete. To simulate them, they must first be marked as *complete*, indicating to Catalyst and ModelingToolkit that they represent finalized models. This can be done using the `complete` function, i.e. by calling `repressilator = complete(repressilator)`. An expanded description on *completeness* can be found [here](@ref completeness_note).
+
 We can check that this is the same model as the one we defined via the DSL as
 follows (this requires that we use the same names for rates, species and the
 system)

docs/src/introduction_to_catalyst/introduction_to_catalyst.md

@@ -96,6 +96,7 @@ Let's now use our `ReactionSystem` to generate and solve a corresponding mass
 action ODE model. We first convert the system to a `ModelingToolkit.ODESystem`
 by
 ```@example tut1
+repressilator = complete(repressilator)
 odesys = convert(ODESystem, repressilator)
 ```
 (Here Latexify is used automatically to display `odesys` in Latex within Markdown
@@ -138,6 +139,7 @@ By passing `repressilator` directly to the `ODEProblem`, Catalyst has to
 (internally) call `convert(ODESystem, repressilator)` again to generate the
 symbolic ODEs. We could instead pass `odesys` directly like
 ```@example tut1
+odesys = complete(odesys)
 oprob2 = ODEProblem(odesys, u₀symmap, tspan, psymmap)
 nothing   # hide
 ```
@@ -152,6 +154,10 @@ underlying problem.
     always be converted to `symbolic` mappings using [`symmap_to_varmap`](@ref),
     see the [Basic Syntax](@ref basic_examples) section for more details.
 
+
+!!! note
+    Above we have used `repressilator = complete(repressilator)` and `odesys = complete(odesys)` to mark these systems as *complete*, indicating to Catalyst and ModelingToolkit that these models are finalized. This must be done before any system is given as input to a `convert` call or some problem type. `ReactionSystem` models created through the @reaction_network` DSL (which is introduced elsewhere, and primarily used throughout these documentation) are always marked as complete when generated. Hence `complete` does not need to be called on them. Symbolically generated `ReactionSystem`s, `ReactionSystem`s generated via the `@network_component` macro, and any ModelingToolkit system generated by `convert` always needs to be manually marked as `complete` as we do for `odesys` above. An expanded description on *completeness* can be found [here](@ref completeness_note).
+
 At this point we are all set to solve the ODEs. We can now use any ODE solver
 from within the
 [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -13,7 +13,7 @@ function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
 
     # Creates NonlinearSystem.
     Catalyst.conservationlaw_errorcheck(rs, vcat(ps, u0))
-    nsys = convert(NonlinearSystem, rs; remove_conserved=true, defaults=Dict(u0))
+    nsys = complete(convert(NonlinearSystem, rs; remove_conserved=true, defaults=Dict(u0)))
 
     # Makes BifurcationProblem (this call goes through the ModelingToolkit-based BifurcationKit extension).
     return BK.BifurcationProblem(nsys, u0_bif, ps, bif_par, args...; plot_var=plot_var,

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -45,7 +45,7 @@ end
 # For a given reaction system, parameter values, and initial conditions, find the polynomial that HC solves to find steady states.
 function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     rs = Catalyst.expand_registered_functions(rs)
-    ns = convert(NonlinearSystem, rs; remove_conserved = true)
+    ns = complete(convert(NonlinearSystem, rs; remove_conserved = true))
     pre_varmap = [symmap_to_varmap(rs,u0)..., symmap_to_varmap(rs,ps)...]
     Catalyst.conservationlaw_errorcheck(rs, pre_varmap)
     p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, parameters(ns); defaults = ModelingToolkit.defaults(ns))

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

@@ -157,8 +157,9 @@ end
 function make_osys(rs::ReactionSystem; remove_conserved=true)
     # Creates the ODESystem corresponding to the ReactionSystem (expanding functions and flattening it).
     # Creates a list of the systems all symbols (unknowns and parameters).
-    rs = Catalyst.expand_registered_functions(flatten(rs))
-    osys = convert(ODESystem, rs; remove_conserved)
+    ModelingToolkit.iscomplete(rs) || error("Identifiability should only be computed for complete systems. A ReactionSystem can be marked as complete using the `complete` function.")
+    rs = complete(Catalyst.expand_registered_functions(flatten(rs)))
+    osys = complete(convert(ODESystem, rs; remove_conserved))
     vars = [unknowns(rs); parameters(rs)]
 
     # Computes equations for system conservation laws.

src/Catalyst.jl

@@ -33,7 +33,7 @@ import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_v
 # internal but needed ModelingToolkit functions
 import ModelingToolkit: check_variables,
                         check_parameters, _iszero, _merge, check_units,
-                        get_unit, check_equations
+                        get_unit, check_equations, iscomplete
 
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
 import MacroTools, Graphs
@@ -48,8 +48,8 @@ end
 function default_t()
     return ModelingToolkit.t_nounits
 end
-const DEFAULT_t = default_t()
-const DEFAULT_IV_SYM = Symbol(DEFAULT_t)
+const DEFAULT_IV = default_t()
+const DEFAULT_IV_SYM = Symbol(DEFAULT_IV)
 export default_t, default_time_deriv
 
 # as used in Catlab
@@ -77,7 +77,7 @@ export ODEProblem,
 const ExprValues = Union{Expr, Symbol, Float64, Int, Bool}
 include("expression_utils.jl")
 include("reaction_network.jl")
-export @reaction_network, @add_reactions, @reaction, @species
+export @reaction_network, @network_component, @add_reactions, @reaction, @species
 
 # registers CRN specific functions using Symbolics.jl
 include("registered_functions.jl")