Merge pull request #678 from SciML/improve_noise_scaling_implementation

[v14 ready] Improve noise scaling implementation

Author: TorkelE

HISTORY.md

@@ -22,8 +22,9 @@ rn = @reaction_network begin
     @parameters η
     k, 2X --> X2, [noise_scaling=η]
 end
-getnoisescaling(rn)
+get_noise_scaling(rn)
 ```
+- `SDEProblem` no longer takes the `noise_scaling` argument (see above for new approach to handle noise scaling).
 - Changed fields of internal `Reaction` structure. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions.
 - Simulation of spatial ODEs now supported. For full details, please see https://github.com/SciML/Catalyst.jl/pull/644 and upcoming documentation. Note that these methods are currently considered alpha, with the interface and approach changing even in non-breaking Catalyst releases.
 - LatticeReactionSystem structure represents a spatial reaction network:

docs/src/catalyst_applications/advanced_simulations.md

@@ -292,8 +292,9 @@ plot(sol)
 
 ## Scaling the noise magnitude in the chemical Langevin equations
 When using the CLE to generate SDEs from a CRN, it can sometimes be desirable to
-scale the magnitude of the noise terms. This can be done by introducing a *noise
-scaling parameter*. First, we simulate a simple two-state CRN model using the
+scale the magnitude of the noise terms. Here, each reaction of the system generates a separate noise term in the CLE. If you require identical scaling for all reactions, the `@default_noise_scaling` option can be used. Else, you can supply a `noise_scaling` metadata for each individual reaction, describing how to scale the noise for that reaction.
+
+We begin with considering the first approach. First, we simulate a simple two-state CRN model using the
 CLE:
 ```@example ex3
 using Catalyst, StochasticDiffEq, Plots
@@ -306,58 +307,61 @@ tspan = (0.0, 10.0)
 p_1 = [:k1 => 1.0, :k2 => 1.0]
 
 sprob_1 = SDEProblem(rn_1, u0, tspan, p_1)
-sol_1 = solve(sprob_1)
+sol_1 = solve(sprob_1, ImplicitEM())
 plot(sol_1; idxs = :X1, ylimit = (0.0, 20.0))
 ```
-Here we can see that the `X` concentration fluctuates around a steady state of $X≈10.0$.
-
-Next, we wish to introduce a noise scaling parameter ,`η`. This will scale the
-noise magnitude so that for $η=0.0$ the system lacks noise (and its SDE
-simulations are identical to its ODE simulations) and for $η=1.0$ noise is not
-scaled (and SDE simulations are identical to as if no noise scaling was used).
-Setting $η<1.0$ will reduce noise and $η>1.0$ will increase noise. The syntax
-for setting a noise scaling parameter `η` is
+Here we can see that the $X$ concentration fluctuates around a steady state of $X≈10.0$.
+
+Next, we wish increase the amount of noise by a factor 2. To do so, we use the `@default_noise_scaling` option, to which we provide the desired scaling 
 ```@example ex3
 rn_2 = @reaction_network begin
-    @parameters η
+    @default_noise_scaling 2
     (k1,k2), X1 <--> X2
 end
-u0 = [:X1 => 10.0, :X2 => 10.0]
-tspan = (0.0, 10.0)
-p_2 = [:k1 => 1.0, :k2 => 1.0, :η => 0.1]
-
-sprob_2 = SDEProblem(rn_2, u0, tspan, p_2; noise_scaling = (@parameters η)[1])
 ```
-Here, we first need to add `η` as a parameter to the system using the
-`@parameters η` option. Next, we pass the `noise_scaling = (@parameters η)[1]`
-argument to the `SDEProblem`. We can now simulate our system and confirm that
-noise is reduced:
+If we re-simulate the system we see that the amount of noise have increased:
 ```@example ex3
-sol_2 = solve(sprob_2)
-plot(sol_2; idxs = :X1, ylimit = (0.0, 20.0))
+sprob_1 = SDEProblem(rn_2, u0, tspan, p_1)
+sol_1 = solve(sprob_1, ImplicitEM())
+plot(sol_1; idxs = :X1, ylimit = (0.0, 20.0))
 ```
 
-Finally, it is possible to set individual noise scaling parameters for each
-reaction of the system. Our model has two reactions (`X1 --> X2` and `X2 -->
-X1`) so we will use two noise scaling parameters, `η1` and `η2`. We use the
-following syntax:
+It is possible to scale the amount of noise using any expression. A common use of this is to set a parameter which determines the amount of noise. Here we create a parameter $η$, and uses its value to scale the noise.
 ```@example ex3
+using Catalyst, StochasticDiffEq, Plots
+
 rn_3 = @reaction_network begin
-    @parameters η1 η2
+    @parameters η
+    @default_noise_scaling η
     (k1,k2), X1 <--> X2
 end
 u0 = [:X1 => 10.0, :X2 => 10.0]
 tspan = (0.0, 10.0)
-p_3 = [:k1 => 1.0, :k2 => 1.0, :η1 => 0.1, :η2 => 1.0]
+p_3 = [:k1 => 1.0, :k2 => 1.0, :η => 0.2]
 
-sprob_3 = SDEProblem(rn_3, u0, tspan, p_3; noise_scaling = @parameters η1 η2)
+sprob_3 = SDEProblem(rn_3, u0, tspan, p_3)
+sol_3 = solve(sprob_3, ImplicitEM())
+plot(sol_3; idxs = :X1, ylimit = (0.0, 20.0))
 ```
-plotting the results, we see that we have less fluctuation than for the first
-simulation, but more as compared to the second one (which is as expected):
+Here we saw how, by setting a small $η$ value, the amount of noise was reduced.
+
+It is possible to use a different noise scaling expression for each reaction. Here, each reaction's noise scaling expression is provided using the `noise_scaling` metadata. In the following example, we use this to tune the noise of for both reactions involving the species $Y$.
+
 ```@example ex3
-sol_3 = solve(sprob_3)
-plot(sol_3; idxs = :X1, ylimit = (0.0, 20.0))
+rn_4 = @reaction_network begin
+    (p, d), 0 <--> X
+    (p, d), 0 <--> Y, ([noise_scaling=0.0, noise_scaling=0.0])
+end
+
+u0_4 = [:X => 10.0, :Y => 10.0]
+tspan = (0.0, 10.0)
+p_4 = [p => 10.0, d => 1.]
+
+sprob_4 = SDEProblem(rn_4, u0_4, tspan, p_4)
+sol_4 = solve(sprob_4, ImplicitEM())
+plot(sol_4; ylimit = (0.0, 20.0))
 ```
+Here, we not that there is n fluctuation in the value of $Y$. If the `@default_noise_scaling` option is used, its value is used for all reactions for which the `noise_scaling` metadata is unused. If `@default_noise_scaling` is not used, the default noise scaling value is `1.0` (i.e. no scaling).
 
 ## Useful plotting options
 Catalyst, just like DifferentialEquations, uses the Plots package for all

src/Catalyst.jl

@@ -91,6 +91,7 @@ export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsm
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
 export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants
 export isequivalent
+export set_default_noise_scaling, get_noise_scaling, has_noise_scaling
 
 # depreciated functions to remove in future releases
 export params, numparams

src/expression_utils.jl

@@ -93,6 +93,19 @@ function is_escaped_expr(expr)
     return (expr isa Expr) && (expr.head == :escape) && (length(expr.args) == 1)
 end
 
+
+### Generic Expression Handling ###
+
+# Convert an expression that is a vector with symbols that have values assigned using `=` 
+# (e.g. :([a=1.0, b=2.0])) to a vector where the assignment instead uses symbols and pairs 
+# (e.g. :([a=>1.0, b=>2.0])). Used to e.g. convert default reaction metadata to a form that can be 
+# evaluated as actual code.
+function expr_equal_vector_to_pairs(expr_vec)
+    pair_vector = :([])
+    foreach(arg -> push!(pair_vector.args, arg.args[1] => arg.args[2]), expr_vec.args) 
+    return pair_vector
+end
+
 ### Old Stuff ###
 
 #This will be called whenever a function stored in funcdict is called.

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, :observables)
+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...)
@@ -361,7 +361,10 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
                        option_lines))
 
     # Reads options.
+    default_reaction_metadata = :([])
     compound_expr, compound_species = read_compound_options(options)
+    check_default_noise_scaling!(default_reaction_metadata, options)
+    default_reaction_metadata = expr_equal_vector_to_pairs(default_reaction_metadata)
 
     # Parses reactions, species, and parameters.
     reactions = get_reactions(reaction_lines)
@@ -422,9 +425,12 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         $observed_vars
         $comps
 
-        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
-                                              $pssym; name = $name, spatial_ivs = $sivs,
-                                              observed = $observed_eqs)
+        Catalyst.remake_ReactionSystem_internal(
+            Catalyst.make_ReactionSystem_internal(
+                $rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
+                $pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs);
+            default_reaction_metadata = $default_reaction_metadata
+        )
     end
 end
 
@@ -637,7 +643,6 @@ end
 
 # Takes a reaction line and creates reaction(s) from it and pushes those to the reaction array.
 # Used to create multiple reactions from, for instance, `k, (X,Y) --> 0`.
-# Handles metadata, e.g. `1.0, Z --> 0, [noisescaling=η]`.
 function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues, 
                          rate::ExprValues, metadata::ExprValues, arrow::Symbol)  
     # The rates, substrates, products, and metadata may be in a tupple form (e.g. `k, (X,Y) --> 0`).
@@ -651,7 +656,7 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues
     for i in 1:maximum(lengs)                       
         # If the `only_use_rate` metadata was not provided, this has to be infered from the arrow used.
         metadata_i = get_tup_arg(metadata, i)
-        if all(arg.args[1] != :only_use_rate for arg in metadata_i.args)
+        if all(arg -> arg.args[1] != :only_use_rate, metadata_i.args)
             push!(metadata_i.args, :(only_use_rate = $(in(arrow, pure_rate_arrows))))
         end
 
@@ -798,6 +803,17 @@ function make_observed_eqs(observables_expr)
     end 
 end
 
+# Checks if the `@default_noise_scaling` option is used. If so, read its input and adds it as a
+# default metadata value to the `default_reaction_metadata` vector.
+function check_default_noise_scaling!(default_reaction_metadata, options)
+    if haskey(options, :default_noise_scaling)
+        if (length(options[:default_noise_scaling].args) != 3) # Becasue of how expressions are, 3 is used.
+            error("@default_noise_scaling should only have a single input, this appears not to be the case: \"$(options[:default_noise_scaling])\"")
+        end
+        push!(default_reaction_metadata.args, :(noise_scaling=$(options[:default_noise_scaling].args[3])))
+    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.