From 97c7c745c54b309f45a4c299b500699157b43dd5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?B=C3=ADma=2C=20Jan?= <bima@merck.com>
Date: Sat, 16 Sep 2023 03:49:46 +0200
Subject: [PATCH] fixes initialization, stepping function

- fix support for initial values specified in acs
- fix support for custom tsteps
- fix req/alloc
---
 src/interface/create.jl |  2 +-
 src/interface/update.jl |  2 +-
 src/solvers.jl          | 16 +++++++---------
 src/state.jl            |  2 +-
 tutorial/basics.jl      |  5 ++++-
 5 files changed, 14 insertions(+), 13 deletions(-)

diff --git a/src/interface/create.jl b/src/interface/create.jl
index 1a82374..fe23701 100644
--- a/src/interface/create.jl
+++ b/src/interface/create.jl
@@ -147,7 +147,7 @@ end
 
 function expand_rate(rate)
     rate = if !(isexpr(rate, :macrocall) && (macroname(rate) == :per_step))
-        :(rand(Poisson(max($rate, 0))))
+        :(rand(Poisson(max(state.dt * $rate, 0))))
     else
         rate.args[3]
     end
diff --git a/src/interface/update.jl b/src/interface/update.jl
index 326c3f6..ce4e811 100644
--- a/src/interface/update.jl
+++ b/src/interface/update.jl
@@ -489,7 +489,7 @@ Add a jump process (with specified Poisson intensity per unit time step) to a mo
 macro jump(acsex, inex, acex)
     return push_to_acs!(
         acsex,
-        Expr(:&&, Expr(:call, :rand, :(Poisson(max(@solverarg(:tstep) * $inex, 0)))), acex),
+        Expr(:&&, Expr(:call, :rand, :(Poisson(max(state.dt * $inex, 0)))), acex),
     )
 end
 
diff --git a/src/solvers.jl b/src/solvers.jl
index 601530a..27cefcf 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -107,9 +107,11 @@ function get_init_satisfied(allocs, qs, state)
                 (reqs[tok.index, i] += tok.stoich)
         end
     end
+    @show 2 reqs
     for i in eachindex(allocs)
         allocs[i] = reqs[i] == 0.0 ? Inf : floor(allocs[i] / reqs[i])
     end
+    @show allocs
     foreach(i -> qs[i] = min(qs[i], minimum(allocs[:, i])), 1:size(reqs, 2))
     foreach(i -> allocs[:, i] .= reqs[:, i] * qs[i], 1:size(reqs, 2))
 
@@ -133,7 +135,7 @@ function evolve!(state)
     qs .= ceil.(Ref(Int), qs)
     @show qs
     for i = 1:nparts(state, :T)
-        new_instances = state.dt * qs[i] + state[i, :transToSpawn]
+        new_instances = qs[i] + state[i, :transToSpawn]
         capacity =
             state[i, :transCapacity] -
             count(t -> t[:transHash] == state[i, :transHash], state.ongoing_transitions)
@@ -165,7 +167,7 @@ function evolve!(state)
     for i = 1:nparts(state, :T)
         qs[i] != 0 && push!(
             state.ongoing_transitions,
-            Transition(state[i, :transName] * "_@$(state.t)", get_sampled_transition(state, i), state.t, qs[i], 0.0),
+            Transition(string(state[i, :transName]) * "_@$(state.t)", get_sampled_transition(state, i), state.t, qs[i], 0.0),
         )
     end
 
@@ -258,11 +260,13 @@ function finish!(state)
                     (in(:rate, tok.modality) ? trans_[:transCycleTime] : 1)
             )
         end
+        
         q = if trans_.state >= trans_[:transCycleTime]
             rand(Distributions.Binomial(Int(trans_.q), trans_[:transProbOfSuccess]))
         else
             0
         end
+
         foreach(
             tok -> (state.u[tok.index] += q * tok.stoich;
             val_reward += state[tok.index, :specReward] * q * tok.stoich),
@@ -378,14 +382,11 @@ function ReactiveNetwork(
 end
 
 function AlgebraicAgents.step!(state::ReactiveNetwork)
-    #du = copy(state.u)
     free_blocked_species!(state)
-    #du .= state.u
     update_observables(state)
     sample_transitions!(state)
     evolve!(state)
     finish!(state)
-    #update_u!(state, u)
     event_action!(state)
 
     push!(
@@ -393,11 +394,8 @@ function AlgebraicAgents.step!(state::ReactiveNetwork)
         (:valuation, state.t, state.u' * [state[i, :specValuation] for i = 1:nparts(state, :S)]),
     )
 
-    #update_u!(state, du)
     save!(state)
-
-    #state.u .= du
-    state.t += state.dt
+    return state.t += state.dt
 end
 
 function AlgebraicAgents._projected_to(state::ReactiveNetwork)
diff --git a/src/state.jl b/src/state.jl
index 2236f51..d70d3bf 100644
--- a/src/state.jl
+++ b/src/state.jl
@@ -217,7 +217,7 @@ end
 ## query the state
 
 t(state::ReactiveNetwork) = state.t
-solverarg(state::ReactiveNetwork, arg) = state.solverargs[arg]
+solverarg(state::ReactiveNetwork, arg) = state.p[arg]
 take(state::ReactiveNetwork, pcs::Symbol) = state.observables[pcs].sampled
 log(state::ReactiveNetwork, msg) = (println(msg); push!(state.log, (:log, msg)))
 state(state::ReactiveNetwork) = state
diff --git a/tutorial/basics.jl b/tutorial/basics.jl
index 1fda02f..e05f7f7 100644
--- a/tutorial/basics.jl
+++ b/tutorial/basics.jl
@@ -10,7 +10,7 @@ end
 
 @prob_init acs X = 10 Y = 20
 @prob_params acs
-@prob_meta acs tspan = 250 dt = 0.1
+@prob_meta acs tspan = 250 dt = 0.11
 
 
 # sol = ReactiveDynamics.solve(prob)
@@ -23,6 +23,9 @@ for i in 1:30
     AlgebraicAgents.step!(prob)
 end
 prob.history_u
+
+
+
 using Plots
 
 @plot sol plot_type = summary