Merge pull request #771 from SciML/enable_dsl_observed

Enable dsl observed

Author: TorkelE

Project.toml

@@ -64,6 +64,7 @@ HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 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"
@@ -77,4 +78,4 @@ 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", "StructuralIdentifiability", "Test", "Unitful"]
+test = ["BifurcationKit", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"]

docs/src/catalyst_functionality/dsl_description.md

@@ -567,3 +567,69 @@ species(rn)
 !!! note
     When using interpolation, expressions like `2$spec` won't work; the
     multiplication symbol must be explicitly included like `2*$spec`.
+
+## Including observables
+Sometimes, one might want to include observable variables. These are variables that can be computed directly from the other system variables (rather than having their values implicitly given through some differential equation). These can be introduced through the `@observables` option.
+
+Let us consider a simple example where two species ($X$ and $Y$) are produced and degraded at constant rates. They can also bind, forming a complex ($XY$). If we want to access the total amount of $X$ in the system we can create an observable that denotes this quantity ($Xtot = X + XY$). Here, we create observables for the total amount of $X$ and $Y$:
+```@example obs1
+using Catalyst # hide
+rn = @reaction_network begin
+  @observables begin
+    Xtot ~ X + XY
+    Ytot ~ Y + XY
+  end
+  (pX,dX), 0 <--> X
+  (pY,dY), 0 <--> Y
+  (kB,kD), X + Y <--> XY
+end
+```
+The `@observables` option is followed by one line for each observable formula (enclosed by a `begin ... end` block). The left-hand sides indicate the observables' names, and the right-hand sides how their values are computed. The two sides are separated by a `~`. 
+
+If we now simulate our model:
+```@example obs1
+using DifferentialEquations # hide
+u0 = [:X => 0.0, :Y => 0.0, :XY => 0.0]
+tspan = (0.0, 10.0)
+ps = [:pX => 1.0, :dX => 0.2, :pY => 1.0, :dY => 0.5, :kB => 1.0, :kD => 0.2]
+oprob = ODEProblem(rn, u0, tspan, ps)
+sol = solve(oprob)
+nothing # hide
+```
+we can index the solution using our observables (just like for [other variables](@ref simulation_structure_interfacing_solutions)). E.g. we can receive a vector with all $Xtot$ values using
+```@example obs1
+sol[:Xtot]
+```
+similarly, we can plot the values of $Xtot$ and $Ytot$ using
+```@example obs1
+plot(sol; idxs=[:Xtot, :Ytot], label=["Total X" "Total Y"])
+```
+
+If we only wish to provide a single observable, the `begin ... end` block is note required. E.g., to record only the total amount of $X$ we can use:
+```@example obs1
+using Catalyst # hide
+rn = @reaction_network begin
+  @observables Xtot ~ X + XY
+  (pX,dX), 0 <--> X
+  (pY,dY), 0 <--> Y
+  (kB,kD), X + Y <--> XY
+end
+```
+
+Finally, some general rules for creating observables:
+- Observables can depend on any species, parameters, or variables, but not on other observables.
+- All observables components appearing on the right side of the `~` must be declared somewhere (i.e., they cannot only appear as a part of the observables formula).
+- Only a single `@observables` option block can be used in each `@reaction_network` call.
+- The left-hand side of the observables expression must be a single symbol, indicating the observable's name.
+- Metadata can, however, be provided, e.g through `@observables (Xtot, [description="Total amount of X"]) ~ X + XY`.
+- The right-hand side of the observables expression can be any valid algebraic expression.
+- Observables are (by default, but this can be changed) considered `variables` (and not `species`). This can be changed by e.g. pre-declaring them using the `@species` option:
+```@example obs2
+using Catalyst # hide
+rn = @reaction_network begin
+  @species Xtot(t)
+  @observables Xtot ~ X1 + X2
+  (k1,k2), X1 <--> X2
+end
+nothing # hide
+```

src/expression_utils.jl

@@ -24,6 +24,9 @@ end
 # The independent variables are given as a vector (empty if none given).
 # Does not support e.g. "X [metadata=true]" (when metadata does not have a comma before).
 function find_varinfo_in_declaration(expr)
+    # Handles the $(Expr(:escape, :Y)) case:
+    is_escaped_expr(expr) && (return find_varinfo_in_declaration(expr.args[1]))
+
     # Case: X
     (expr isa Symbol) && (return expr, [], nothing, nothing)                          
     # Case: X(t)         
@@ -84,6 +87,10 @@ function insert_independent_variable(expr_in, iv_expr)
     return expr
 end
 
+# Checks if an expression is an escaped expression (e.g. on the form `$(Expr(:escape, :Y))`)
+function is_escaped_expr(expr)
+    return (expr isa Expr) && (expr.head == :escape) && (length(expr.args) == 1)
+end
 
 ### Old Stuff ###
 

src/reaction_network.jl

@@ -123,7 +123,7 @@ const forbidden_variables_error = let
 end
 
 # Declares the keys used for various options.
-const option_keys = (:species, :parameters, :variables, :ivs, :compounds)
+const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables)
 
 ### The @species macro, basically a copy of the @variables macro. ###
 macro species(ex...)
@@ -355,10 +355,11 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     option_lines = Expr[x for x in ex.args if x.head == :macrocall]
 
     # Get macro options.
+    length(unique(arg.args[1] for arg in option_lines)) < length(option_lines) && error("Some options where given multiple times.")
     options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg,
                        option_lines))
 
-    # Reads compounds options.
+    # Reads options.
     compound_expr, compound_species = read_compound_options(options)
 
     # Parses reactions, species, and parameters.
@@ -378,7 +379,11 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     end
     tiv = ivs[1]
     sivs = (length(ivs) > 1) ? Expr(:vect, ivs[2:end]...) : nothing
+    all_ivs = (isnothing(sivs) ? [tiv] : [tiv; sivs.args])
 
+    # Reads more options.
+    observed_vars, observed_eqs, obs_syms = read_observed_options(options, [species_declared; variables], all_ivs)
+    
     declared_syms = Set(Iterators.flatten((parameters_declared, species_declared,
                                            variables)))
     species_extracted, parameters_extracted = extract_species_and_parameters!(reactions,
@@ -408,17 +413,17 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         push!(rxexprs.args, get_rxexprs(reaction))
     end
 
-    # Returns the rephrased expression.
     quote
         $ps
         $ivexpr
         $vars
         $sps
+        $observed_vars
         $comps
 
-        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, union($spssym, $varssym, $compssym),
-                                              $pssym; name = $name,
-                                              spatial_ivs = $sivs)
+        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
+                                              $pssym; name = $name, spatial_ivs = $sivs,
+                                              observed = $observed_eqs)
     end
 end
 
@@ -682,6 +687,88 @@ function recursive_find_reactants!(ex::ExprValues, mult::ExprValues,
     reactants
 end
 
+### DSL Options Handling ###
+# Most options handled in previous sections, when code re-organised, these should ideally be moved to the same place.
+
+# Reads the observables options. Outputs an expression ofr creating the obervable variables, and a vector of observable equations.
+function read_observed_options(options, species_n_vars_declared, ivs_sorted)
+    if haskey(options, :observables)
+        # Gets list of observable equations and prepares variable declaration expression.
+        # (`options[:observables]` inlucdes `@observables`, `.args[3]` removes this part)
+        observed_eqs = make_observed_eqs(options[:observables].args[3])
+        observed_vars = Expr(:block, :(@variables))
+        obs_syms = :([])
+        
+        for (idx, obs_eq) in enumerate(observed_eqs.args)
+            # Extract the observable, checks errors, and continues the loop if the observable has been declared.        
+            obs_name, ivs, defaults, metadata = find_varinfo_in_declaration(obs_eq.args[2])
+            isempty(ivs) || error("An observable ($obs_name) was given independent variable(s). These should not be given, as they are inferred automatically.")
+            isnothing(defaults) || error("An observable ($obs_name) was given a default value. This is forbidden.")
+            in(obs_name, forbidden_symbols_error) && error("A forbidden symbol ($(obs_eq.args[2])) was used as an observable name.")
+            
+            # Error checks.
+            if (obs_name in species_n_vars_declared) && is_escaped_expr(obs_eq.args[2])
+                error("An interpoalted observable have been used, which has also been explicitly delcared within the system using eitehr @species or @variables. This is not permited.")
+            end
+            if ((obs_name in species_n_vars_declared) || is_escaped_expr(obs_eq.args[2])) && !isnothing(metadata)
+                error("Metadata was provided to observable $obs_name in the `@observables` macro. However, the obervable was also declared separately (using either @species or @variables). When this is done, metadata should instead be provided within the original @species or @variable declaration.")
+            end
+
+            # This bits adds the observables to the @variables vector which is given as output.
+            # For Observables that have already been declared using @species/@variables,
+            # or are interpolated, this parts should not be carried out.
+            if !((obs_name in species_n_vars_declared) || is_escaped_expr(obs_eq.args[2]))
+                # Appends (..) to the observable (which is later replaced with the extracted ivs).
+                # Adds the observable to the first line of the output expression (starting with `@variables`).
+                    obs_expr = insert_independent_variable(obs_eq.args[2], :(..))
+                    push!(observed_vars.args[1].args, obs_expr)
+                
+                # Adds a line to the `observed_vars` expression, setting the ivs for this observable.
+                # Cannot extract directly using e.g. "getfield.(dependants_structs, :reactant)" because 
+                # then we get something like :([:X1, :X2]), rather than :([X1, X2]).
+                dep_var_expr = :(filter(!MT.isparameter, Symbolics.get_variables($(obs_eq.args[3]))))
+                ivs_get_expr = :(unique(reduce(vcat,[arguments(MT.unwrap(dep)) for dep in $dep_var_expr])))    
+                ivs_get_expr_sorted = :(sort($(ivs_get_expr); by = iv -> findfirst(MT.getname(iv) == ivs for ivs in $ivs_sorted)))
+                push!(observed_vars.args, :($obs_name = $(obs_name)($(ivs_get_expr_sorted)...)))
+            end
+
+            # In case metadata was given, this must be cleared from `observed_eqs`.
+            # For interpolated observables (I.e. $X ~ ...) this should and cannot be done.
+            is_escaped_expr(obs_eq.args[2]) || (observed_eqs.args[idx].args[2] = obs_name)
+
+            # Adds the observable to the list of observable names.
+            # This is required for filtering away so these are not added to the ReactionSystem's species list.
+            # Again, avoid this check if we have interpoalted teh variable.
+            is_escaped_expr(obs_eq.args[2]) || push!(obs_syms.args, obs_name)
+        end
+
+        # If nothing was added to `observed_vars`, it has to be modified not to throw an error.
+        (length(observed_vars.args) == 1) && (observed_vars = :())
+    else  
+        # If option is not used, return empty expression and vector.
+        observed_vars = :()
+        observed_eqs = :([])
+        obs_syms = :([])
+    end
+    return observed_vars, observed_eqs, obs_syms
+end
+
+# From the input to the @observables options, creates a vector containing one equation for each observable.
+# Checks separate cases for "@obervables O ~ ..." and "@obervables begin ... end". Other cases errors.
+function make_observed_eqs(observables_expr)
+    if observables_expr.head == :call
+        return :([$(observables_expr)])
+    elseif observables_expr.head == :block
+        observed_eqs = :([])
+        for arg in observables_expr.args
+            push!(observed_eqs.args, arg)
+        end
+        return observed_eqs
+    else
+        error("Malformed observables option usage: $(observables_expr).")
+    end 
+end
+
 ### Functionality for expanding function call to custom and specific functions ###
 
 #Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.

src/reactionsystem.jl

@@ -516,7 +516,7 @@ struct ReactionSystem{V <: NetworkProperties} <:
     function ReactionSystem(eqs, rxs, iv, sivs, states, spcs, ps, var_to_name, observed,
                             name, systems, defaults, connection_type, nps, cls, cevs, devs,
                             complete::Bool = false; checks::Bool = true)
-
+                            
         # unit checks are for ODEs and Reactions only currently
         nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
         if checks && isempty(sivs)