up

Author: TorkelE

docs/src/catalyst_functionality/dsl_description.md

@@ -623,4 +623,13 @@ Finally, some general rules for creating observables:
 - 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.
+- 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

@@ -382,8 +382,8 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     all_ivs = (isnothing(sivs) ? [tiv] : [tiv; sivs.args])
 
     # Reads more options.
-    observed_vars, observed_eqs = read_observed_options(options, [species_declared; variables], all_ivs)
-
+    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,
@@ -421,7 +421,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         $observed_vars
         $comps
 
-        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, union($spssym, $varssym, $compssym),
+        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
                                               $pssym; name = $name, spatial_ivs = $sivs,
                                               observed = $observed_eqs)
     end
@@ -691,49 +691,60 @@ end
 # 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_declared, ivs_sorted)
+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.")
-            if (obs_name in species_declared) 
+            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.")
-                continue
             end
 
-            # 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)...)))
-
-            # In case metadata was given, this must be cleared from `observed_eqs`
-            observed_eqs.args[idx].args[2] = obs_name
+            # 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.
+            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
+    return observed_vars, observed_eqs, obs_syms
 end
 
 # From the input to the @observables options, creates a vector containing one equation for each observable.

src/reactionsystem.jl

@@ -516,11 +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)
-        # Filters away any potential obervables from `states` and `spcs`.
-        obs_vars = [obs_eq.lhs for obs_eq in observed]
-        states = filter(state -> !any(isequal(state, obs_var) for obs_var in obs_vars), states)
-        spcs = filter(spc -> !any(isequal(spc, obs_var) for obs_var in obs_vars), spcs)
-
+                            
         # 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)

test/dsl/dsl_options.jl

@@ -510,6 +510,57 @@ let
     @test isequal(states(rn3), states(rn4))
 end
 
+# Tests for interpolation into the observables option.
+let
+    # Interpolation into lhs.
+    @species X [description="An observable"]
+    rn1 = @reaction_network begin
+        @observables $X ~ X1 + X2
+        (k1, k2), X1 <--> X2
+    end
+    @test isequal(observed(rn1)[1].lhs, X)
+    @test getdescription(rn1.X) == "An observable"
+    @test isspecies(rn1.X)
+    @test length(states(rn1)) == 2
+
+    # Interpolation into rhs.
+    @parameters n [description="A parameter"]
+    @variables t
+    @species S(t)
+    rn2 = @reaction_network begin
+        @observables Stot ~ $S + $n*Sn
+        (kB, kD), $n*S <--> Sn
+    end
+    @unpack Stot, Sn, kD, kB = rn2
+
+    u0 = Dict([S => 5.0, Sn => 1.0])
+    ps = Dict([n => 2, kB => 1.0, kD => 1.0])
+    oprob = ODEProblem(rn2, u0, (0.0, 1.0), ps)
+
+    @test issetequal(Symbolics.get_variables(observed(rn2)[1].rhs), [S, n, Sn])
+    @test oprob[Stot] == u0[S] + ps[n]*u0[Sn]
+    @test length(states(rn2)) == 2
+end
+
+# Tests specific declaration of Observables as species/variables
+let
+    rn = @reaction_network begin
+        @species X(t)
+        @variables Y(t)
+        @observables begin
+            X ~ X + 2X2
+            Y ~ Y1 + Y2
+            Z ~ X + Y
+        end
+        (kB,kD), 2X <--> X2
+        (k1,k2), Y1 <--> Y2
+    end
+    
+    @test isspecies(rn.X)
+    @test !isspecies(rn.Y)
+    @test !isspecies(rn.Z)
+end
+
 # Tests various erroneous declarations throw errors.
 let 
     # Independent variable in @observables.