add @per_step deterministic "rate"

Author: thevolatilebit

docs/src/index.md

@@ -26,6 +26,10 @@ We list common species attributes:
 | `specInitUncertainty` | uncertainty about variable's initial state (modelled as Gaussian standard deviation) |
 | `specInitVal` | initial value of a variable |
 
+Moreover, it is possible to specify the semantics of the "rate" term. By default, at each time step `n ~ Poisson(rate * dt)` instances of a given transition will be spawned. If you want to specify the rate in terms of a cycle time, you may want to use `@ct(cycle_time)`, e.g., `@ct(ex), A --> B, ...`. This is a shorthand for `1/ex, A --> B, ...`.
+
+For deterministic "rates", use `@per_step(ex)`. Here, `ex` evaluates to a deterministic number (ceiled to the nearest integer) of a transition's instances to spawn per a single integrator's step. However, note that in this case, the number doesn't scale with the step length! Moreover
+
 ```@docs
 @add_species
 @aka

src/interface/create.jl

@@ -115,16 +115,22 @@ function recursively_expand_actions!(evs, condex, event)
     else push!(evs, Event(condex, event)) end
 end
 
-function prune_rate(rate)
-    if (isexpr(rate, :macrocall) && (macroname(rate) ∈ prettynames[:transCycleTime]))
-        :(1/$(rate.args[3]))
-    else rate end
+function expand_rate(rate)
+    rate = if !(isexpr(rate, :macrocall) && (macroname(rate) == :per_step))
+        :(rand(Poisson(max($rate, 0))))
+    else rate.args[3] end
+
+    postwalk(rate) do ex
+        if (isexpr(ex, :macrocall) && (macroname(ex) ∈ prettynames[:transCycleTime]))
+            :(1/$(ex.args[3]))
+        else ex end
+    end
 end
 
 function get_transitions!(trans, reactants, pcs, exs)
     args = empty(defargs[:T])
     (rate, r_line) = exs[1:2]
-    rxs = prune_reaction_line!(pcs, reactants, r_line); rate = prune_rate(rate)
+    rxs = prune_reaction_line!(pcs, reactants, r_line); rate = expand_rate(rate)
     rxs = rxs isa Tuple ? tuple.(rate.args, rxs) : ((rate, rxs),)
 
     exs = exs[3:end]; empty!(args)

src/solvers.jl

@@ -99,8 +99,9 @@ function evolve!(u, state)
     reqs = zeros(nparts(state, :S), nparts(state, :T)); qs = zeros(nparts(state, :T))
 
     foreach(i -> qs[i] = state[i, :transRate] * state[i, :transMultiplier], 1:nparts(state, :T))
+    qs .= ceil.(Ref(Int), qs)
     for i in 1:nparts(state, :T)
-        new_instances = rand(Poisson(max(state.solverargs[:tstep] * qs[i], 0))) + state[i, :transToSpawn]
+        new_instances = state.solverargs[:tstep]*qs[i] + state[i, :transToSpawn]
         capacity = state[i, :transCapacity] - count(t -> t[:transHash] == state[i, :transHash], state.ongoing_transitions)
         (capacity < new_instances) && add_to_spawn!(state, state[i, :transHash], new_instances - capacity)
         qs[i] = min(capacity, new_instances)

test/tutorial_tests.jl

@@ -4,4 +4,5 @@ using ReactiveDynamics
 @safeinclude "joins" "../tutorial/joins/joins.jl"
 @safeinclude "loadsave" "../tutorial/loadsave/loadsave.jl"
 @safeinclude "optimize" "../tutorial/optimize/optimize.jl"
-@safeinclude "solution wrap" "../tutorial/optimize/optimize_custom.jl"
+@safeinclude "solution wrap" "../tutorial/optimize/optimize_custom.jl"
+@safeinclude "toy pharma model" "../tutorial/toy_pharma_model.jl"

tutorial/example.jl

@@ -21,13 +21,13 @@ end
 @valuation sir_acs I=.1 # for testing purpose only
 
 # atomic data: initial values, parameter values
-@prob_init sir_acs S=999 I=10 R=0
+@prob_init sir_acs S=999 I=100 R=100
 u0 = [999, 10, 0] # alternative specification
 @prob_init sir_acs u0
 @prob_params sir_acs β=0.0001 ν=0.01 γ=5
 @prob_meta sir_acs tspan=100
 
-prob = @problematize sir_acs
+#prob = @problematize sir_acs
 prob = @problematize sir_acs tspan=200
 
 sol = @solve prob trajectories=20

tutorial/toy_pharma_model.jl

@@ -0,0 +1,65 @@
+using ReactiveDynamics
+
+# model dynamics
+toy_pharma_model = @ReactionNetwork begin
+        α(candidate_compound, marketed_drug, κ), 3*@conserved(scientist) + @rate(budget) --> candidate_compound, name=>discovery, probability=>.3, cycletime=>10., priority=>.5
+        β(candidate_compound, marketed_drug), candidate_compound + 5*@conserved(scientist) + 2*@rate(budget) --> marketed_drug + 5*budget, name=>dx2market, probability=>.5+.001*@t(), cycletime=>4
+        γ*marketed_drug, marketed_drug --> ∅, name=>drug_killed
+end
+
+@periodic toy_pharma_model 1. budget += 11*marketed_drug
+
+@register α(number_candidate_compounds, number_marketed_drugs, κ) = κ + exp(-number_candidate_compounds) + exp(-number_marketed_drugs)
+@register β(number_candidate_compounds, number_marketed_drugs) = number_candidate_compounds + exp(-number_marketed_drugs)
+
+# simulation parameters
+## initial values
+@prob_init toy_pharma_model candidate_compound=5 marketed_drug=6 scientist=20 budget=100
+## parameters
+@prob_params toy_pharma_model κ=4 γ=.1
+## other arguments passed to the solver
+@prob_meta toy_pharma_model tspan=250 dt=.1
+
+
+prob = @problematize toy_pharma_model
+
+sol = @solve prob trajectories=20
+
+using Plots; plotly()
+
+@plot sol plot_type=summary
+
+@plot sol plot_type=summary show=:marketed_drug
+
+## for deterministic rates 
+
+# model dynamics
+toy_pharma_model = @ReactionNetwork begin
+        @per_step(α(candidate_compound, marketed_drug, κ)), 3*@conserved(scientist) + @rate(budget) --> candidate_compound, name=>discovery, probability=>.3, cycletime=>10., priority=>.5
+        @per_step(β(candidate_compound, marketed_drug)), candidate_compound + 5*@conserved(scientist) + 2*@rate(budget) --> marketed_drug + 5*budget, name=>dx2market, probability=>.5+.001*@t(), cycletime=>4
+        @per_step(γ*marketed_drug), marketed_drug --> ∅, name=>drug_killed
+end
+
+@periodic toy_pharma_model 0. budget += 11*marketed_drug
+
+@register α(number_candidate_compounds, number_marketed_drugs, κ) = κ + exp(-number_candidate_compounds) + exp(-number_marketed_drugs)
+@register β(number_candidate_compounds, number_marketed_drugs) = number_candidate_compounds + exp(-number_marketed_drugs)
+
+# simulation parameters
+## initial values
+@prob_init toy_pharma_model candidate_compound=5 marketed_drug=6 scientist=20 budget=100
+## parameters
+@prob_params toy_pharma_model κ=4 γ=.1
+## other arguments passed to the solver
+@prob_meta toy_pharma_model tspan=250
+
+
+prob = @problematize toy_pharma_model
+
+sol = @solve prob trajectories=20
+
+using Plots; plotly()
+
+@plot sol plot_type=summary
+
+@plot sol plot_type=summary show=:marketed_drug