up

Author: TorkelE

docs/src/model_creation/dsl_basics.md

@@ -264,6 +264,17 @@ Catalyst comes with the following predefined functions:
 - The repressive Hill function: $hillr(X,v,K,n) = v * (K^n)/(X^n + K^n)$.
 - The activating/repressive Hill function: $hillar(X,Y,v,K,n) = v * (X^n)/(X^n + Y^n + K^n)$.
 
+### [Registration of non-algebraic functions](@id dsl_description_nonconstant_rates_function_registration)
+Previously we showed how user-defined functions [can be used in rates directly](@ref dsl_description_nonconstant_rates_available_functions). For functions containing more complicated syntax (e.g. `for` loops or `if` statements), we must add an additional step: registering it using the `@register_symbolic` macro. Below we define a function which output depends on whether `X` is smaller or larger than a threshold value. Next, we register it using `@register_symbolic`, after which we can use it within the DSL.
+```@example dsl_basics
+threshold_func(X) = (X < 10) ? X : 10.0
+@register_symbolic threshold_func(X)
+rn = @reaction_network begin
+    threshold_func(X), 0 --> X
+    d, X --> 0
+end
+```
+
 ### [Time-dependant rates](@id dsl_description_nonconstant_rates_time)
 Previously we have assumed that the rates are independent of the time variable, $t$. However, time-dependent reactions are also possible. Here, simply use `t` to represent the time variable. E.g., to create a production/degradation model where the production rate decays as time progresses, we can use:
 ```@example dsl_basics

docs/src/model_creation/functional_parameters.md

@@ -1,16 +1,18 @@
 # [Inputs and time-dependent (or functional) parameters](@id time_dependent_parameters)
-Catalyst supports the usage of "functional parameters". In practice these are parameters that are given by time-dependent functions, representing a way to inject custom functions into models. They can be used when rates depend on real data, or to represent complicated functions (which use e.g. `for` loops or random number generation). Here, the function's values are declared as a data interpolation, which is then used as the functional parameter's value in the simulation. On this page, we first show how to create time-dependent functional parameters, and then give an example where the functional parameter depends on a species value.
+Catalyst supports the usage of "functional parameters". In practice, these are parameters that are given by (typically) time-dependent functions (they can also depend on e.g. species values, as discussed [here](@ref functional_parameters_sir)). They are a way to inject custom functions into models. Functional parameters can be used when rates depend on real data, or to represent complicated functions (which use e.g. `for` loops or random number generation). Here, the function's values are declared as a data interpolation (which interpolates discrete samples to a continuous function). This is then used as the functional parameter's value in the simulation. This tutorial first shows how to create time-dependent functional parameters, and then gives an example where the functional parameter depends on a species value.
+
+An alternative approach for representing complicated functions is by [using `@register_symbolic`](@ref dsl_description_nonconstant_rates_function_registration).
 
 ## [Basic example](@id functional_parameters_basic_example)
-Let us first consider an easy, quick-start example. We will consider a simple [birth-death model](@ref basic_CRN_library_bd), but where the birth rate is determined by an input parameter (for which the value depends on time). First, we [define the input parameter programmatically](@ref programmatic_CRN_construction), and its values across all time points using the [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) package. In this example we will use the input function $pIn(t) = (2 + t)/(1 + t)$.
+Let us first consider an easy, quick-start example (the next section will discuss what is going on in more detail). We will consider a simple [birth-death model](@ref basic_CRN_library_bd), but where the birth rate is determined by an input parameter (for which the value depends on time). First, we [define the input parameter programmatically](@ref programmatic_CRN_construction), and its values across all time points using the [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) package. In this example we will use the input function $pIn(t) = (2 + t)/(1 + t)$. Finally, we plot the input function, demonstrating how while it is defined at discrete points, DataInterpolations.jl generalises this to a continuous function.
 ```@example functional_parameters_basic_example
-using Catalyst, DataInterpolations
+using Catalyst, DataInterpolations, Plots
 t = default_t()
 tend = 10.0
 ts = collect(0.0:0.01:tend)
 spline = LinearInterpolation((2 .+ ts) ./ (1 .+ ts), ts)
 @parameters (pIn::typeof(spline))(..)
-nothing # hide
+plot(spline)
 ```
 Next, we create our model, [interpolating](@ref dsl_advanced_options_symbolics_and_DSL_interpolation) the input parameter into it (making it a function of `t`).
 ```@example functional_parameters_basic_example
@@ -20,9 +22,9 @@ bd_model = @reaction_network begin
     d, X --> 0
 end
 ```
-Finally, we can simulate our model as normal, but setting the value of the `pIn` parameter to our interpolated data.
+Finally, we can simulate our model as normal (but where we set the value of the `pIn` parameter to our interpolated data).
 ```@example functional_parameters_basic_example
-using OrdinaryDiffEqDefault, Plots
+using OrdinaryDiffEqDefault
 u0 = [:X => 0.5]
 ps = [:d => 2.0, :pIn => spline]
 oprob = ODEProblem(bd_model, u0, tend, ps)
@@ -48,7 +50,7 @@ using DataInterpolations
 interpolated_light = LinearInterpolation(light, ts)
 plot(interpolated_light)
 ```
-We are now ready to declare our model. We will consider a protein with an active and an inactive form ($Pₐ$ and $Pᵢ$) where the activation is driven by the presence of sunlight. In this example we we create our model using the [programmatic approach](@ref programmatic_CRN_construction). Do note the special syntax we use to declare our input parameter, where we both designate it as a generic function and its type as the type of our interpolated input. Also note that, within the model, we mark the input parameter (`light_in`) a function of `t`.
+We are now ready to declare our model. We will consider a protein with an active and an inactive form ($Pₐ$ and $Pᵢ$) where the activation is driven by the presence of sunlight. In this example we we create our model using the [programmatic approach](@ref programmatic_CRN_construction). Do note the special syntax we use to declare our input parameter, where we both designate it as a generic function and its type as the type of our interpolated input. Also note that, within the model, we mark the input parameter (`light_in`) as a function of `t`.
 ```@example functional_parameters_circ_rhythm
 using Catalyst
 t = default_t()
@@ -73,24 +75,24 @@ plot(sol)
 ```
 
 ### [Interpolating the input into the DSL](@id functional_parameters_circ_rhythm_dsl)
-It is possible to use time-dependent inputs when creating models [through the DSL](@ref dsl_description) as well. However, it can still be convenient to declare the input parameter programmatically as above. Next, we form an expression of it as a function of time, and then [interpolate](@ref dsl_advanced_options_symbolics_and_DSL_interpolation) it into our DSL-declaration:
+It is possible to use time-dependent inputs when creating models [through the DSL](@ref dsl_description) as well. However, it can still be convenient to declare the input parameter programmatically as above. Using it, we form an expression of it as a function of time, and then [interpolate](@ref dsl_advanced_options_symbolics_and_DSL_interpolation) it into our DSL-declaration:
 ```@example functional_parameters_circ_rhythm
 input = light_in(t)
 rs_dsl = @reaction_network rs begin
     (kA*$input, kD), Pᵢ <--> Pₐ
 end
 ```
-We can confirm that this model is identical to our programmatic one (and should we wish to, we can simulate it using identical syntax syntax).
+We can confirm that this model is identical to our programmatic one (and should we wish to, we can simulate it using identical syntax).
 ```@example functional_parameters_circ_rhythm
 rs == rs_dsl
 ```
 
 ## [Non-time functional parameters](@id functional_parameters_sir)
-Previously we have demonstrated functional parameters that are a function of time. However, functional parameters can be functions of any variable (however, currently, more than one argument is not supported). Here we will demonstrate this using a [SIR model](@ref basic_CRN_library_sir), but instead of having the infection rate scale linearly with the number of infected individuals, we instead assume we have measured data of the infection rate (as dependent on the number of infected individuals) and wish to use this instead. Normally we use the following infection reaction in the SIR model:
+Previously we have demonstrated functional parameters that are functions of time. However, functional parameters can be functions of any variable (however, currently, more than one argument is not supported). Here we will demonstrate this using a [SIR model](@ref basic_CRN_library_sir), but instead of having the infection rate scale linearly with the number of infected individuals, we instead assume we have measured data of the infection rate (as dependent on the number of infected individuals) and wish to use this instead. Normally we use the following infection reaction in the SIR model:
 ```julia
 @reaction k1, S + I --> 2I
 ```
-In practise, this is identical to
+For ODE models, this would give the same equations as
 ```julia
 @reaction k1*I, S --> I
 ```
@@ -104,27 +106,28 @@ We start by declaring the functional parameter that describes how the infection
 ```@example functional_parameters_sir
 using DataInterpolations, Plots
 I_grid = collect(0.0:5.0:100.0)
-I_meassured = 300.0 *(0.8*rand(length(I_grid)) .+ 0.6) .* I_grid ./ (300 .+ I_grid)
-I_rate = LinearInterpolation(I_meassured, I_grid)
-plot(I_rate; label = "Meassured infection rate")
+I_measured = 300.0 *(0.8*rand(length(I_grid)) .+ 0.6) .* I_grid ./ (300 .+ I_grid)
+I_rate = LinearInterpolation(I_measured, I_grid)
+plot(I_rate; label = "Measured infection rate")
 plot!(I_grid, I_grid; label = "Normal SIR infection rate")
 ```
-Next, we create our model (using the DSL-approach). As `I_rate` will be a function of $I$ we will need to declare this species first as well.
+Next, we create our model (using the DSL approach). As `I_rate` will be a function of $I$ we will need to declare this species first as well.
 ```@example functional_parameters_sir
 using Catalyst
+@parameters (inf_rate::typeof(I_rate))(..)
 @species I(default_t())
-inf_rate = I_rate(I)
+inf_rate_in = inf_rate(I)
 sir = @reaction_network rs begin
-    k1*$inf_rate, S --> I
+    k1*$inf_rate_in, S --> I
     k2, I --> R
 end
 nothing # hide
 ```
-Finally, we can simualte our model.
+Finally, we can simulate our model.
 ```@example functional_parameters_sir
 using OrdinaryDiffEqDefault
 u0 = [:S => 99.0, :I => 1.0, :R => 0.0]
-ps = [:k1 => 0.002, :k2 => 0.01, inf_rate => I_rate]
+ps = [:k1 => 0.002, :k2 => 0.01, :inf_rate => I_rate]
 oprob = ODEProblem(sir, u0, 250.0, ps)
 sol = solve(oprob)
 plot(sol)
@@ -145,4 +148,4 @@ plot(spline_linear)
 plot!(spline_cubuc)
 plot!(spline_const)
 ```
-Finally, DataInterpolation.jl also allows various [extrapolation methods](https://docs.sciml.ai/DataInterpolations/stable/extrapolation_methods/).
+Finally, DataInterpolations.jl also allows various [extrapolation methods](https://docs.sciml.ai/DataInterpolations/stable/extrapolation_methods/).

test/reactionsystem_core/custom_crn_functions.jl

@@ -49,7 +49,7 @@ let
         ps = rnd_ps(custom_function_network_1, rng; factor)
         t = rand(rng)
         @test f_eval(custom_function_network_1, u0, ps, t) ≈ f_eval(custom_function_network_2, u0, ps, t)
-        @test_broken jac_eval(custom_function_network_1, u0, ps, t) ≈ jac_eval(custom_function_network_2, u0, ps, t) # Weird error dur to some SciML update. Reported in https://github.com/SciML/ModelingToolkit.jl/issues/3447.
+        @test jac_eval(custom_function_network_1, u0, ps, t) ≈ jac_eval(custom_function_network_2, u0, ps, t)
         @test g_eval(custom_function_network_1, u0, ps, t) ≈ g_eval(custom_function_network_2, u0, ps, t)
     end
 end

test/runtests.jl

@@ -21,10 +21,10 @@ end
         @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end
         @time @safetestset "Symbolic Stoichiometry" begin include("reactionsystem_core/symbolic_stoichiometry.jl") end
         @time @safetestset "Parameter Type Designation" begin include("reactionsystem_core/parameter_type_designation.jl") end
-        @time @safetestset "Functional Parameters" begin include("reactionsystem_core/functional_parameters.jl") end
         @time @safetestset "Custom CRN Functions" begin include("reactionsystem_core/custom_crn_functions.jl") end
         @time @safetestset "Coupled CRN/Equation Systems" begin include("reactionsystem_core/coupled_equation_crn_systems.jl") end
         @time @safetestset "Events" begin include("reactionsystem_core/events.jl") end
+        @time @safetestset "Functional Parameters" begin include("reactionsystem_core/functional_parameters.jl") end
 
         # Tests model creation via the @reaction_network DSL.
         @time @safetestset "DSL Basic Model Construction" begin include("dsl/dsl_basic_model_construction.jl") end