Merge branch 'master' into src___refactoring___dsl_file

Author: TorkelE

HISTORY.md

@@ -7,6 +7,29 @@
   (at the time the release is made). If you need a dependency version increased,
   please open an issue and we can update it and make a new Catalyst release once
   testing against the newer dependency version is complete. 
+- New formula for inferring variables from equations (declared using the `@equations` options) in the DSL. The order of inference of species/variables/parameters is now:
+    (1) Every symbol explicitly declared using `@species`, `@variables`, and `@parameters` are assigned to the correct category.
+    (2) Every symbol used as a reaction reactant is inferred as a species.
+    (3) Every symbol not declared in (1) or (2) that occurs in an expression provided after `@equations` is inferred as a variable.
+    (4) Every symbol not declared in (1), (2), or (3) that occurs either as a reaction rate or stoichiometric coefficient is inferred to be a parameter.
+E.g. in
+```julia
+@reaction_network begin
+    @equations V1 + S ~ V2^2
+    (p + S + V1), S --> 0
+end
+```
+`S` is inferred as a species, `V1` and `V2` as variables, and `p` as a parameter. The previous special cases for the `@observables`, `@compounds`, and `@differentials` options still hold. Finally, the `@require_declaration` options (described in more detail below) can now be used to require everything to be explicitly declared.
+- New formula for determining whether the default differentials have been used within an `@equations` option. Now, if any expression `D(...)` is encountered (where `...` can be anything), this is inferred as usage of the default differential D. E.g. in the following equations `D` is inferred as a differential with respect to the default independent variable:
+```julia
+@reaction_network begin
+  @equations D(V) + V ~ 1
+end
+@reaction_network begin
+  @equations D(D(V)) ~ 1
+end
+```
+Please note that this cannot be used at the same time as `D` is used to represent a species, variable, or parameter (including is these are implicitly designated as such by e.g. appearing as a reaction reactant). 
 - Array symbolics support is more consistent with ModelingToolkit v9. Parameter
   arrays are no longer scalarized by Catalyst, while species and variables
   arrays still are (as in ModelingToolkit). As such, parameter arrays should now

docs/src/inverse_problems/examples/ode_fitting_oscillation.md

@@ -56,7 +56,7 @@ function optimise_p(pinit, tend)
         newprob = remake(prob; tspan = (0.0, tend), p = p)
         sol = Array(solve(newprob, Rosenbrock23(); saveat = newtimes))
         loss = sum(abs2, sol .- sample_vals[:, 1:size(sol,2)])
-        return loss, sol
+        return loss
     end
 
     # optimize for the parameters that minimize the loss

src/dsl.jl

@@ -354,6 +354,10 @@ function make_reaction_system(ex::Expr, name)
     default_reaction_metadata = read_default_noise_scaling_option(options)
     combinatoric_ratelaws = read_combinatoric_ratelaws_option(options)
 
+    # Reads observables.
+    observed_vars, observed_eqs, obs_syms = read_observed_options(
+        options, [species_declared; variables], all_ivs; requiredec)
+
     # Checks for input errors.
     forbidden_symbol_check(union(species, parameters))
     forbidden_variable_check(variables)
@@ -732,7 +736,7 @@ end
 # `vars_extracted`: A vector with extracted variables (lhs in pure differential equations only).
 # `dtexpr`: If a differential equation is defined, the default derivative (D ~ Differential(t)) must be defined.
 # `equations`: a vector with the equations provided.
-function read_equations_options(options, variables_declared; requiredec = false)
+function read_equations_options(options, syms_declared, parameters_extracted; requiredec = false)
     # Prepares the equations. First, extracts equations from provided option (converting to block form if required).
     # Next, uses MTK's `parse_equations!` function to split input into a vector with the equations.
     eqs_input = haskey(options, :equations) ? options[:equations].args[3] : :(begin end)
@@ -744,34 +748,40 @@ function read_equations_options(options, variables_declared; requiredec = false)
     # Loops through all equations, checks for lhs of the form `D(X) ~ ...`.
     # When this is the case, the variable X and differential D are extracted (for automatic declaration).
     # Also performs simple error checks.
-    vars_extracted = Vector{Symbol}()
+    vars_extracted = OrderedSet{Union{Symbol, Expr}}()
     add_default_diff = false
     for eq in equations
         if (eq.head != :call) || (eq.args[1] != :~)
             error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")
         end
 
-        # Checks if the equation have the format D(X) ~ ... (where X is a symbol). This means that the
-        # default differential has been used. X is added as a declared variable to the system, and
-        # we make a note that a differential D = Differential(iv) should be made as well.
-        lhs = eq.args[2]
-        # if lhs: is an expression. Is a function call. The function's name is D. Calls a single symbol.
-        if (lhs isa Expr) && (lhs.head == :call) && (lhs.args[1] == :D) &&
-           (lhs.args[2] isa Symbol)
-            diff_var = lhs.args[2]
-            if in(diff_var, forbidden_symbols_error)
-                error("A forbidden symbol ($(diff_var)) was used as an variable in this differential equation: $eq")
-            elseif (!in(diff_var, variables_declared)) && requiredec
-                throw(UndeclaredSymbolicError(
-                                              "Unrecognized symbol $(diff_var) was used as a variable in an equation: \"$eq\". Since the @require_declaration flag is set, all variables in equations must be explicitly declared via @variables, @species, or @parameters."))
-            else
-                add_default_diff = true
-                in(diff_var, variables_declared) || push!(vars_extracted, diff_var)
-            end
+        # If the default differential (`D`) is used, record that it should be decalred later on.
+        if (:D ∉ union(syms_declared, parameters_extracted)) && find_D_call(eq)
+            requiredec && throw(UndeclaredSymbolicError(
+                "Unrecognized symbol D was used as a differential in an equation: \"$eq\". Since the @require_declaration flag is set, all differentials in equations must be explicitly declared using the @differentials option."))
+            add_default_diff = true
+            push!(syms_declared, :D)
         end
+
+        # Any undecalred symbolic variables encountered should be extracted as variables.
+        add_syms_from_expr!(vars_extracted, eq, syms_declared)
+        (!isempty(vars_extracted) && requiredec) && throw(UndeclaredSymbolicError(
+            "Unrecognized symbolic variables $(join(vars_extracted, ", ")) detected in equation expression: \"$(string(eq))\". Since the flag @require_declaration is declared, all symbolic variables must be explicitly declared with the @species, @variables, and @parameters options."))
     end
 
-    return vars_extracted, add_default_diff, equations
+    return collect(vars_extracted), add_default_diff, equations
+end
+
+# Searches an expresion `expr` and returns true if it have any subexpression `D(...)` (where `...` can be anything).
+# Used to determine whether the default differential D has been used in any equation provided to `@equations`.
+function find_D_call(expr)
+    return if Base.isexpr(expr, :call) && expr.args[1] == :D
+        true
+    elseif expr isa Expr
+        any(find_D_call, expr.args)
+    else
+        false
+    end
 end
 
 # Creates an expression declaring differentials. Here, `tiv` is the time independent variables,

src/network_analysis.jl

@@ -206,9 +206,6 @@ function laplacianmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false)
     D * K
 end
 
-Base.zero(::Type{Union{R, Symbolics.BasicSymbolic{Real}}}) where R <: Real = zero(R)
-Base.one(::Type{Union{R, Symbolics.BasicSymbolic{Real}}}) where R <: Real = one(R)
-
 @doc raw"""
     fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false)
 
@@ -218,15 +215,20 @@ Base.one(::Type{Union{R, Symbolics.BasicSymbolic{Real}}}) where R <: Real = one(
     **Warning**: Unlike other Catalyst functions, the `fluxmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work.
 """
 function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false)
+    deps = Set()
+    for (i, rx) in enumerate(reactions(rn))
+        empty!(deps)
+        get_variables!(deps, rx.rate, species(rn))
+        (!isempty(deps)) && (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `fluxmat` cannot support rate constants of this form."))
+    end
+
     rates = if isempty(pmap)
         reactionrates(rn)
     else
         substitutevals(rn, pmap, parameters(rn), reactionrates(rn))
     end
 
     rcmap = reactioncomplexmap(rn)
-    nc = length(rcmap)
-    nr = length(rates)
     mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates)
     if sparse
         return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates)
@@ -295,6 +297,10 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric
     rxs = reactions(rn)
     sm = speciesmap(rn)
 
+    for rx in rxs
+        !ismassaction(rx, rn) && error("The supplied ReactionSystem has non-mass action reaction $rx. The `massactionvector` can only be constructed for mass action networks.")
+    end
+
     specs = if isempty(scmap) 
         species(rn)
     else
@@ -936,7 +942,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, re
     complexes, D = reactioncomplexes(rs)
     img = incidencematgraph(rs)
     undir_img = SimpleGraph(incidencematgraph(rs))
-    K = ratematrix(rs, pmap)
+    K = adjacencymat(rs, pmap)
 
     spanning_forest = Graphs.kruskal_mst(undir_img)
     outofforest_edges = setdiff(collect(edges(undir_img)), spanning_forest)
@@ -993,52 +999,25 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap)
 end
 
 """
-    iscomplexbalanced(rs::ReactionSystem, parametermap)
+    iscomplexbalanced(rs::ReactionSystem, parametermap; sparse = false)
 
 Constructively compute whether a network will have complex-balanced equilibrium
 solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). 
-Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
-"""
-function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict)
-    if length(parametermap) != numparams(rs)
-        error("Incorrect number of parameters specified.")
-    end
-
-    pmap = symmap_to_varmap(rs, parametermap)
-    pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap)
 
-    sm = speciesmap(rs)
-    cm = reactioncomplexmap(rs)
-    complexes, D = reactioncomplexes(rs)
+Requires the specification of values for the parameters/rate constants. Accepts a dictionary, vector, or tuple of parameter-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
+"""
+function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict; sparse = false)
     rxns = reactions(rs)
-    nc = length(complexes)
-    nr = numreactions(rs)
-    nm = numspecies(rs)
-
     if !all(r -> ismassaction(r, rs), rxns)
         error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.")
     end
 
-    rates = [substitute(rate, pmap) for rate in reactionrates(rs)]
-
-    # Construct kinetic matrix, K
-    K = zeros(nr, nc)
-    for c in 1:nc
-        complex = complexes[c]
-        for (r, dir) in cm[complex]
-            rxn = rxns[r]
-            if dir == -1
-                K[r, c] = rates[r]
-            end
-        end
-    end
-
-    L = -D * K
-    S = netstoichmat(rs)
+    D = incidencemat(rs; sparse)
+    S = netstoichmat(rs; sparse)
 
     # Compute ρ using the matrix-tree theorem
     g = incidencematgraph(rs)
-    R = ratematrix(rs, rates)
+    R = adjacencymat(rs, parametermap; sparse)
     ρ = matrixtree(g, R)
 
     # Determine if 1) ρ is positive and 2) D^T Ln ρ lies in the image of S^T
@@ -1069,50 +1048,81 @@ function iscomplexbalanced(rs::ReactionSystem, parametermap)
 end
 
 """
-    ratematrix(rs::ReactionSystem, parametermap)
+    adjacencymat(rs::ReactionSystem, pmap = Dict(); sparse = false)
 
     Given a reaction system with n complexes, outputs an n-by-n matrix where R_{ij} is the rate 
     constant of the reaction between complex i and complex j. Accepts a dictionary, vector, or tuple 
     of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
+
+    Equivalent to the adjacency matrix of the weighted graph whose weights are the 
+    reaction rates. 
+    Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. 
+
+    The difference between this matrix and [`fluxmat`](@ref) is quite subtle. The non-zero entries to both matrices are the rate constants. However:
+        - In `fluxmat`, the rows represent complexes and columns represent reactions, and an entry (c, r) is non-zero if reaction r's source complex is c. 
+        - In `adjacencymat`, the rows and columns both represent complexes, and an entry (c1, c2) is non-zero if there is a reaction c1 --> c2.
 """
-function ratematrix(rs::ReactionSystem, rates::Vector{Float64})
-    complexes, D = reactioncomplexes(rs)
-    n = length(complexes)
-    rxns = reactions(rs)
-    ratematrix = zeros(n, n)
+function adjacencymat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false)
+    deps = Set()
+    for (i, rx) in enumerate(reactions(rn))
+        empty!(deps)
+        get_variables!(deps, rx.rate, species(rn))
+        (!isempty(deps)) && (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `adjacencymat` cannot support rate constants of this form."))
+    end
 
-    for r in 1:length(rxns)
-        rxn = rxns[r]
-        s = findfirst(==(-1), @view D[:, r])
-        p = findfirst(==(1), @view D[:, r])
-        ratematrix[s, p] = rates[r]
+    rates = if isempty(pmap)
+        reactionrates(rn)
+    else
+        substitutevals(rn, pmap, parameters(rn), reactionrates(rn))
+    end
+    mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates)
+
+    if sparse
+        return adjacencymat(SparseMatrixCSC{mtype, Int}, incidencemat(rn), rates)
+    else
+        return adjacencymat(Matrix{mtype}, incidencemat(rn), rates)
     end
-    ratematrix
 end
 
-function ratematrix(rs::ReactionSystem, parametermap::Dict)
-    if length(parametermap) != numparams(rs)
-        error("Incorrect number of parameters specified.")
+function adjacencymat(::Type{SparseMatrixCSC{T, Int}}, D, rates) where T
+    Is = Int[]
+    Js = Int[]
+    Vs = T[]
+    nc = size(D, 1)
+
+    for r in 1:length(rates)
+        s = findfirst(==(-1), @view D[:, r])
+        p = findfirst(==(1), @view D[:, r])
+        push!(Is, s)
+        push!(Js, p)
+        push!(Vs, rates[r])
     end
+    A = sparse(Is, Js, Vs, nc, nc)
+end
 
-    pmap = symmap_to_varmap(rs, parametermap)
-    pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap)
+function adjacencymat(::Type{Matrix{T}}, D, rates) where T
+    nc = size(D, 1)
+    A = zeros(T, nc, nc)
 
-    rates = [substitute(rate, pmap) for rate in reactionrates(rs)]
-    ratematrix(rs, rates)
+    for r in 1:length(rates)
+        s = findfirst(==(-1), @view D[:, r])
+        p = findfirst(==(1), @view D[:, r])
+        A[s, p] = rates[r]
+    end
+    A
 end
 
-function ratematrix(rs::ReactionSystem, parametermap::Vector{<:Pair})
-    pdict = Dict(parametermap)
-    ratematrix(rs, pdict)
+function adjacencymat(rn::ReactionSystem, pmap::Vector{<:Pair}; sparse=false)
+    pdict = Dict(pmap)
+    adjacencymat(rn, pdict; sparse)
 end
 
-function ratematrix(rs::ReactionSystem, parametermap::Tuple)
-    pdict = Dict(parametermap)
-    ratematrix(rs, pdict)
+function adjacencymat(rn::ReactionSystem, pmap::Tuple; sparse=false)
+    pdict = Dict(pmap)
+    adjacencymat(rn, pdict; sparse)
 end
 
-function ratematrix(rs::ReactionSystem, parametermap)
+function adjacencymat(rn::ReactionSystem, pmap)
     error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.")
 end