update with new completeness handling

Author: TorkelE

docs/src/catalyst_functionality/compositional_modeling.md

@@ -6,13 +6,12 @@ identically repressed genes, and how to use compositional modeling to create
 compartments.
 
 ## A note on *completeness*
-Catalyst `ReactionSystem` can either be *complete* or *incomplete*. *By default 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.
+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 set a `ReactionSystem` created via the DSL as complete, use the `@incomplete` option:
+To create incomplete `ReactionSystem`s using the DSL as complete, use the `@network_component` instead of `@reaction_network`:
 ```@example ex0
 using Catalyst
-degradation_component = @reaction_network begin
-  @incomplete
+degradation_component = @network_component begin
   d, X --> 0
 end
 ```
@@ -40,12 +39,10 @@ Catalyst supports two ModelingToolkit interfaces for composing multiple
 extending a system is the `extend` command
 ```@example ex1
 using Catalyst
-basern = @reaction_network rn1 begin
-  @incomplete
+basern = @network_component rn1 begin
   k, A + B --> C
 end
-newrn = @reaction_network rn2 begin
-  @incomplete
+newrn = @network_component rn2 begin
   r, C --> A + B
 end
 @named rn = extend(newrn, basern)
@@ -58,8 +55,7 @@ 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
-  @incomplete
+newestrn = @network_component rn3 begin
   v, A + D --> 2D
 end
 @named rn = compose(basern, [newrn, newestrn])
@@ -131,8 +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
-    @incomplete
+  @network_component $name begin
     hillr($R,α,K,n), ∅ --> m
     (δ,γ), m <--> ∅
     β, m --> m + P
@@ -189,15 +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
-  @incomplete
+nuc = @network_component nuc begin
   α, G --> G + M
   (κ₊/V,κ₋), D + G <--> DG
 end
 
 # translation and dimerization
-cyto = @reaction_network cyto begin
-  @incomplete
+cyto = @network_component cyto begin
   β, M --> M + P
   (k₊/V,k₋), 2P <--> D
   σ, P --> 0
@@ -206,8 +199,7 @@ end
 
 # export reactions,
 # γ,δ=probability per time to be exported/imported
-model = @reaction_network model begin
-  @incomplete
+model = @network_component model begin
   γ, $(nuc.M) --> $(cyto.M)
   δ, $(cyto.D) --> $(nuc.D)
 end

src/Catalyst.jl

@@ -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")

src/reaction_network.jl

@@ -1,45 +1,4 @@
 ### Temporary deprecation warning - Eventually to be removed. ###
-deprication_message = """
-@reaction_network notation where parameters are declared after "end", e.g. like:
-
-```julia
-@reaction_network begin
-    p, 0 --> X
-    d, X --> 0
-end p d
-```
-
-has been deprecated in favor of a notation where the parameters are inferred, e.g:
-
-```julia
-@reaction_network begin
-    p, 0 --> X
-    d, X --> 0
-end
-```
-
-Parameters and species can be explicitly indicated using the @parameters and @species
-macros, e.g:
-
-```julia
-@reaction_network begin
-    @parameters p d
-    @species X(t)
-    p, 0 --> X
-    d, X --> 0
-end
-```
-"""
-macro reaction_network(name::Symbol, ex::Expr, parameters...)
-    error(deprication_message)
-end
-macro reaction_network(name::Expr, ex::Expr, parameters...)
-    error(deprication_message)
-end
-macro reaction_network(ex::Expr, parameters...)
-    error(deprication_message)
-end
-
 """
 Macro that inputs an expression corresponding to a reaction network and outputs
 a `ReactionNetwork` that can be used as input to generation of ODE, SDE, and
@@ -123,7 +82,7 @@ const forbidden_variables_error = let
 end
 
 # Declares the keys used for various options.
-const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling, :incomplete)
+const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling)
 
 ### The @species macro, basically a copy of the @variables macro. ###
 macro species(ex...)
@@ -186,19 +145,55 @@ emptyrn = @reaction_network empty
 # an empty network with random generated name
 emptyrn = @reaction_network
 ```
+
+ReactionSystems generated through `@reaction_network` are compelte.
 """
 macro reaction_network(name::Symbol, ex::Expr)
-    make_reaction_system(MacroTools.striplines(ex); name = :($(QuoteNode(name))))
+    :(complete($(make_reaction_system(MacroTools.striplines(ex); name = :($(QuoteNode(name)))))))
 end
 
 # allows @reaction_network $name begin ... to interpolate variables storing a name
 macro reaction_network(name::Expr, ex::Expr)
-    make_reaction_system(MacroTools.striplines(ex); name = :($(esc(name.args[1]))))
+    :(complete($(make_reaction_system(MacroTools.striplines(ex); name = :($(esc(name.args[1])))))))
 end
 
 macro reaction_network(ex::Expr)
     ex = MacroTools.striplines(ex)
 
+    # no name but equations: @reaction_network begin ... end ...
+    if ex.head == :block
+        :(complete($(make_reaction_system(ex))))
+    else  # empty but has interpolated name: @reaction_network $name
+        networkname = :($(esc(ex.args[1])))
+        return Expr(:block, :(@parameters t),
+                    :(complete(ReactionSystem(Reaction[], t, [], []; name = $networkname))))
+    end
+end
+
+#Returns a empty network (with, or without, a declared name)
+macro reaction_network(name::Symbol = gensym(:ReactionSystem))
+    return Expr(:block, :(@parameters t),
+                :(complete(ReactionSystem(Reaction[], t, [], []; name = $(QuoteNode(name))))))
+end
+
+
+"""
+    @network_component
+
+As @reaction_network, but the output system is not complete.
+"""
+macro network_component(name::Symbol, ex::Expr)
+    make_reaction_system(MacroTools.striplines(ex); name = :($(QuoteNode(name))))
+end
+
+# allows @reaction_network $name begin ... to interpolate variables storing a name
+macro network_component(name::Expr, ex::Expr)
+    make_reaction_system(MacroTools.striplines(ex); name = :($(esc(name.args[1]))))
+end
+
+macro network_component(ex::Expr)
+    ex = MacroTools.striplines(ex)
+
     # no name but equations: @reaction_network begin ... end ...
     if ex.head == :block
         make_reaction_system(ex)
@@ -210,12 +205,12 @@ macro reaction_network(ex::Expr)
 end
 
 #Returns a empty network (with, or without, a declared name)
-# @reaction_network name
-macro reaction_network(name::Symbol = gensym(:ReactionSystem))
+macro network_component(name::Symbol = gensym(:ReactionSystem))
     return Expr(:block, :(@parameters t),
                 :(ReactionSystem(Reaction[], t, [], []; name = $(QuoteNode(name)))))
 end
 
+
 ### Macros used for manipulating, and successively builing up, reaction systems. ###
 @doc raw"""
     @reaction

src/reactionsystem.jl

@@ -544,7 +544,7 @@ struct ReactionSystem{V <: NetworkProperties} <:
     # inner constructor is considered private and may change between non-breaking releases.
     function ReactionSystem(eqs, rxs, iv, sivs, unknowns, spcs, ps, var_to_name, observed,
                             name, systems, defaults, connection_type, nps, cls, cevs, devs,
-                            metadata = nothing, complete::Bool = false; checks::Bool = true)
+                            metadata, complete; checks::Bool = true)
 
         # unit checks are for ODEs and Reactions only currently
         nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
@@ -610,8 +610,7 @@ function ReactionSystem(eqs, iv, unknowns, ps;
                         spatial_ivs = nothing,
                         continuous_events = nothing,
                         discrete_events = nothing,
-                        metadata = nothing,
-                        complete = true)
+                        metadata = nothing)
 
     name === nothing &&
         throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
@@ -681,7 +680,7 @@ function ReactionSystem(eqs, iv, unknowns, ps;
 
     ReactionSystem(eqs′, rxs, iv′, sivs′, unknowns′, spcs, ps′, var_to_name, observed, name,
                    systems, defaults, connection_type, nps, combinatoric_ratelaws,
-                   ccallbacks, dcallbacks, metadata, complete; checks = checks)
+                   ccallbacks, dcallbacks, metadata, false; checks = checks)
 end
 
 function ReactionSystem(rxs::Vector, iv = Catalyst.DEFAULT_IV; kwargs...)
@@ -1472,6 +1471,8 @@ function spatial_convert_err(rs::ReactionSystem, systype)
     isspatial(rs) && error("Conversion to $systype is not supported for spatial networks.")
 end
 
+COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be converted to other system types. A ReactionSystem can be marked as complete using the `compelte` function."
+
 """
 ```julia
 Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
@@ -1493,8 +1494,9 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
                       include_zero_odes = true, remove_conserved = false, checks = false,
                       default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
                       kwargs...)
+    iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, ODESystem)
-    fullrs = Catalyst.flatten(rs; complete = true)
+    fullrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(fullrs)
     ists, ispcs = get_indep_sts(fullrs, remove_conserved)
     eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
@@ -1509,7 +1511,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
               continuous_events = MT.get_continuous_events(fullrs),
               discrete_events = MT.get_discrete_events(fullrs),
               kwargs...)
-    return iscomplete(rs) ? complete(osys) : osys
+    return osys
 end
 
 """
@@ -1534,6 +1536,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
                       include_zero_odes = true, remove_conserved = false, checks = false,
                       default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
                       kwargs...)
+    iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, NonlinearSystem)
     fullrs = Catalyst.flatten(rs; complete = true)
     remove_conserved && conservationlaws(fullrs)
@@ -1549,7 +1552,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
                     defaults = _merge(defaults,defs),
                     checks,
                     kwargs...)
-    return iscomplete(rs) ? complete(nlsys) : nlsys
+    return nlsys
 end
 
 """
@@ -1576,6 +1579,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
                       include_zero_odes = true, checks = false, remove_conserved = false,
                       default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
                       kwargs...)
+    iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, SDESystem)
 
     flatrs = Catalyst.flatten(rs; complete = true)
@@ -1600,7 +1604,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
               continuous_events = MT.get_continuous_events(flatrs),
               discrete_events = MT.get_discrete_events(flatrs),
               kwargs...)
-    return iscomplete(rs) ? complete(ssys) : ssys
+    return ssys
 end
 
 """
@@ -1626,6 +1630,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
                       remove_conserved = nothing, checks = false,
                       default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
                       kwargs...)
+                      iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, JumpSystem)
 
     (remove_conserved !== nothing) &&
@@ -1651,7 +1656,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
                checks,
                discrete_events = MT.discrete_events(flatrs),
                kwargs...)
-    return iscomplete(rs) ? complete(jsys) : jsys
+    return jsys
 end
 
 ### Converts a reaction system to ODE or SDE problems ###
@@ -1667,6 +1672,7 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
     pmap = symmap_to_varmap(rs, p)
     osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
                    remove_conserved)
+    osys = complete(osys)
     return ODEProblem(osys, u0map, tspan, pmap, args...; check_length, kwargs...)
 end
 
@@ -1681,6 +1687,7 @@ function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0,
     pmap = symmap_to_varmap(rs, p)
     nlsys = convert(NonlinearSystem, rs; name, combinatoric_ratelaws, include_zero_odes,
                     checks, remove_conserved)
+    nlsys = complete(nlsys)
     return NonlinearProblem(nlsys, u0map, pmap, args...; check_length, kwargs...)
 end
 
@@ -1695,6 +1702,7 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan,
     pmap = symmap_to_varmap(rs, p)
     sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws,
                       include_zero_odes, checks, remove_conserved)
+    sde_sys = complete(sde_sys)
     p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(rs))
     return SDEProblem(sde_sys, u0map, tspan, pmap, args...; check_length,
                       noise_rate_prototype = p_matrix, kwargs...)
@@ -1709,6 +1717,7 @@ function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0, tspan::Tuple,
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, checks)
+    jsys = complete(jsys)
     return DiscreteProblem(jsys, u0map, tspan, pmap, args...; kwargs...)
 end
 
@@ -1718,6 +1727,7 @@ function JumpProcesses.JumpProblem(rs::ReactionSystem, prob, aggregator, args...
                                    combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
                                    checks = false, kwargs...)
     jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, checks)
+    jsys = complete(jsys)
     return JumpProblem(jsys, prob, aggregator, args...; kwargs...)
 end
 
@@ -1732,6 +1742,7 @@ function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0,
     pmap = symmap_to_varmap(rs, p)
     osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
                    remove_conserved)
+    osys = complete(osys)
     return SteadyStateProblem(osys, u0map, pmap, args...; check_length, kwargs...)
 end
 
@@ -1792,7 +1803,7 @@ Notes:
 - The default value of `combinatoric_ratelaws` will be the logical or of all
   `ReactionSystem`s.
 """
-function MT.flatten(rs::ReactionSystem; name = nameof(rs), complete = false)
+function MT.flatten(rs::ReactionSystem; name = nameof(rs))
     isempty(get_systems(rs)) && return rs
 
     # right now only NonlinearSystems and ODESystems can be handled as subsystems
@@ -1810,8 +1821,7 @@ function MT.flatten(rs::ReactionSystem; name = nameof(rs), complete = false)
                    balanced_bc_check = false,
                    spatial_ivs = get_sivs(rs),
                    continuous_events = MT.continuous_events(rs),
-                   discrete_events = MT.discrete_events(rs),
-                   complete = complete)
+                   discrete_events = MT.discre)
 end
 
 """
@@ -1868,6 +1878,5 @@ function ModelingToolkit.extend(sys::MT.AbstractSystem, rs::ReactionSystem;
                    balanced_bc_check = false,
                    spatial_ivs = sivs,
                    continuous_events,
-                   discrete_events,
-                   complete = false)
+                   discrete_events)
 end