solve

Project.toml

@@ -9,18 +9,12 @@ Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 Crayons = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-DocumenterMarkdown = "997ab1e6-3595-5248-9280-8efb232c3433"
 GeneratedExpressions = "84d730a5-1eb9-4187-a799-27dd07f33a14"
 IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781"
 PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
@@ -32,26 +26,19 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 CSV = "0.10"
 Catlab = "0.14"
 ComponentArrays = "0.14"
 Crayons = "4.1"
 DataFrames = "1.6"
-DiffEqBase = "6.128"
-DifferentialEquations = "7.9"
 Distributions = "0.25"
-Documenter = "0.27"
-DocumenterMarkdown = "0.2"
 GeneratedExpressions = "0.1"
 IJulia = "1.24"
 JLD2 = "0.4"
 JSON = "0.21"
 MacroTools = "0.5"
-NLopt = "1.0"
-OrdinaryDiffEq = "6.55"
 Plots = "1.39"
 Pluto = "0.19"
 PlutoUI = "0.7"

src/ReactiveDynamics.jl

@@ -3,7 +3,6 @@ module ReactiveDynamics
 using Catlab, Catlab.CategoricalAlgebra, Catlab.Present
 using Reexport
 using MacroTools
-using NLopt
 using ComponentArrays
 
 @reexport using GeneratedExpressions

src/solvers.jl

@@ -1,8 +1,3 @@
-# assortment of SciML-compatible problem solvers
-
-export DiscreteProblem
-
-using DiffEqBase, DifferentialEquations
 using Distributions
 using Random
 
@@ -37,7 +32,7 @@ function get_reqs_ongoing!(reqs, qs, state)
         for tok in state.ongoing_transitions[i][:transLHS]
             in(:rate, tok.modality) &&
                 (state.ongoing_transitions[i][:transCycleTime] > 0) &&
-                (reqs[tok.index, i] += qs[i] * tok.stoich * state.solverargs[:tstep])
+                (reqs[tok.index, i] += qs[i] * tok.stoich * state.dt)
             in(:nonblock, tok.modality) && (reqs[tok.index, i] += qs[i] * tok.stoich)
         end
     end
@@ -124,9 +119,8 @@ end
 """
 Evolve transitions, spawn new transitions.
 """
-function evolve!(u, state)
-    update_u!(state, u)
-    actual_allocs = zero(u)
+function evolve!(state)
+    actual_allocs = zero(state.u)
 
     ## schedule new transitions
     reqs = zeros(nparts(state, :S), nparts(state, :T))
@@ -137,8 +131,9 @@ function evolve!(u, state)
         1:nparts(state, :T),
     )
     qs .= ceil.(Ref(Int), qs)
+    @show qs
     for i = 1:nparts(state, :T)
-        new_instances = state.solverargs[:tstep] * qs[i] + state[i, :transToSpawn]
+        new_instances = state.dt * qs[i] + state[i, :transToSpawn]
         capacity =
             state[i, :transCapacity] -
             count(t -> t[:transHash] == state[i, :transHash], state.ongoing_transitions)
@@ -148,10 +143,13 @@ function evolve!(u, state)
     end
 
     reqs = get_reqs_init!(reqs, qs, state)
+    @show reqs
     allocs =
-        get_allocs!(reqs, u, state, state[:, :transPriority], state.solverargs[:strategy])
+        get_allocs!(reqs, state.u, state, state[:, :transPriority], state.p[:strategy])
+    @show allocs
     qs .= get_init_satisfied(allocs, qs, state)
-
+    @show qs
+    println("====")
     push!(
         state.log,
         (
@@ -160,7 +158,7 @@ function evolve!(u, state)
             [(hash, q) for (hash, q) in zip(state[:, :transHash], qs)]...,
         ),
     )
-    u .-= sum(allocs; dims = 2)
+    state.u .-= sum(allocs; dims = 2)
     actual_allocs .+= sum(allocs; dims = 2)
 
     # add spawned transitions to the heap
@@ -171,18 +169,17 @@ function evolve!(u, state)
         )
     end
 
-    update_u!(state, u)
     ## evolve ongoing transitions 
     reqs = zeros(nparts(state, :S), length(state.ongoing_transitions))
     qs = map(t -> t.q, state.ongoing_transitions)
 
     get_reqs_ongoing!(reqs, qs, state)
     allocs = get_allocs!(
         reqs,
-        u,
+        state.u,
         state,
         map(t -> t[:transPriority], state.ongoing_transitions),
-        state.solverargs[:strategy],
+        state.p[:strategy],
     )
     qs .= get_frac_satisfied(allocs, reqs, state)
     push!(
@@ -196,11 +193,11 @@ function evolve!(u, state)
             ]...,
         ),
     )
-    u .-= sum(allocs; dims = 2)
+    state.u .-= sum(allocs; dims = 2)
     actual_allocs .+= sum(allocs; dims = 2)
 
     foreach(
-        i -> state.ongoing_transitions[i].state += qs[i] * state.solverargs[:tstep],
+        i -> state.ongoing_transitions[i].state += qs[i] * state.dt,
         eachindex(state.ongoing_transitions),
     )
 
@@ -229,8 +226,7 @@ function event_action!(state)
 end
 
 # collect terminated transitions
-function finish!(u, state)
-    update_u!(state, u)
+function finish!(state)
     val_reward = 0
     terminated_all = Dict{Symbol,Float64}()
     terminated_success = Dict{Symbol,Float64}()
@@ -256,7 +252,7 @@ function finish!(u, state)
         end
         for tok in trans_[:transLHS]
             in(:conserved, tok.modality) && (
-                u[tok.index] +=
+                state.u[tok.index] +=
                     trans_.q *
                     tok.stoich *
                     (in(:rate, tok.modality) ? trans_[:transCycleTime] : 1)
@@ -268,12 +264,11 @@ function finish!(u, state)
             0
         end
         foreach(
-            tok -> (u[tok.index] += q * tok.stoich;
+            tok -> (state.u[tok.index] += q * tok.stoich;
             val_reward += state[tok.index, :specReward] * q * tok.stoich),
             toks_rhs,
         )
 
-        update_u!(state, u)
         context_eval(state, trans_.trans[:transPostAction])
         terminated_all[trans_[:transHash]] =
             get(terminated_all, trans_[:transHash], 0) + trans_.q
@@ -289,7 +284,7 @@ function finish!(u, state)
     push!(state.log, (:terminated_success, state.t, terminated_success...))
     push!(state.log, (:valuation_reward, state.t, val_reward))
 
-    return u
+    return state.u
 end
 
 function free_blocked_species!(state)
@@ -306,15 +301,15 @@ function get_tcontrol(tspan, args)
     tunit = get(args, :tunit, oneunit(tspan))
     tspan = tspan / tunit
 
-    tstep = get(args, :tstep, haskey(args, :tstops) ? tspan / args[:tstops] : tunit) / tunit
+    dt = get(args, :dt, haskey(args, :tstops) ? tspan / args[:tstops] : tunit) / tunit
 
-    return ((0.0, tspan), tstep)
+    return ((0.0, tspan), dt)
 end
 
 function ReactiveNetwork(
     acs::ReactionNetwork,
     u0 = Dict(),
-    p = DiffEqBase.NullParameters();
+    p = Dict();
     name = "reactive_network",
     kwargs...,
 )
@@ -325,39 +320,21 @@ function ReactiveNetwork(
     ])
     merge!(keywords, Dict(collect(kwargs)))
     merge!(keywords, Dict(:strategy => get(keywords, :alloc_strategy, :weighted)))
+
     keywords[:tspan], keywords[:tstep] = get_tcontrol(keywords[:tspan], keywords)
 
     acs = remove_choose(acs)
     attrs, transitions, wrap_fun = compile_attrs(acs)
-
-    init_u!(state)
-    save!(state)
+    transition_recipes = transitions
+    u0_init = zeros(nparts(acs, :S))
 
-    u0_init = zeros(nparts(state, :S))
-
-    u0 isa Dict && foreach(
-        i ->
-            if !isnothing(acs[i, :specName]) && haskey(u0, acs[i, :specName])
-                u0_init[i] = u0[acs[i, :specName]]
-            end,
-        1:nparts(state, :S),
-    )
-
-    p_ = p == DiffEqBase.NullParameters() ? Dict() : Dict(k => v for (k, v) in p)
-    prob = remake(
-        prob;
-        u0 = prob.u0,
-        tspan = keywords[:tspan],
-        dt = get(keywords, :tstep, 1),
-        p = merge(
-            prob.p,
-            p_,
-            Dict(
-                :tstep => get(keywords, :tstep, 1),
-                :strategy => get(keywords, :alloc_strategy, :weighted),
-            ),
-        ),
-    )
+    for i in 1:nparts(acs, :S)
+        if !isnothing(acs[i, :specName]) && haskey(u0, acs[i, :specName])
+            u0_init[i] = u0[acs[i, :specName]]
+        else
+            u0_init[i] = acs[i, :specInitVal]
+        end
+    end
 
     ongoing_transitions = Transition[]
     log = NamedTuple[]
@@ -369,21 +346,19 @@ function ReactiveNetwork(
         ) ∪ [:transLHS, :transRHS, :transToSpawn, :transHash]
     transitions = Dict{Symbol,Vector}(a => [] for a in transitions_attrs)
 
-    return ReactiveNetwork(
+    network =  ReactiveNetwork(
         name,
         acs,
         attrs,
         transition_recipes,
         u0_init, 
         merge(
-            prob.p,
-            p_,
+            p,
             Dict(
                 :tstep => get(keywords, :tstep, 1),
                 :strategy => get(keywords, :alloc_strategy, :weighted),
             ),
         ),
-        t, 
         keywords[:tspan][1],
         keywords[:tspan],
         get(keywords, :tstep, 1),
@@ -396,31 +371,32 @@ function ReactiveNetwork(
         Vector{Float64}[],
         Float64[],
     )
+
+    save!(network)
+
+    return network
 end
 
 function AlgebraicAgents.step!(state::ReactiveNetwork)
-    du = copy(state.u)
+    #du = copy(state.u)
     free_blocked_species!(state)
-    du .= state.u
+    #du .= state.u
     update_observables(state)
     sample_transitions!(state)
-    evolve!(du, state)
-    finish!(du, state)
-    update_u!(state, du)
+    evolve!(state)
+    finish!(state)
+    #update_u!(state, u)
     event_action!(state)
 
-    du .= state.u
     push!(
         state.log,
-        (:valuation, t, du' * [state[i, :specValuation] for i = 1:nparts(state, :S)]),
+        (:valuation, state.t, state.u' * [state[i, :specValuation] for i = 1:nparts(state, :S)]),
     )
 
-    t = (state.t += state.solverargs[:tstep])
-    update_u!(state, du)
+    #update_u!(state, du)
     save!(state)
-    sync_p!(p, state)
 
-    state.u .= du
+    #state.u .= du
     state.t += state.dt
 end
 

src/state.jl

@@ -1,5 +1,3 @@
-using DiffEqBase: NullParameters
-
 using AlgebraicAgents
 
 struct UnfoldedReactant
@@ -77,7 +75,7 @@ function init_u!(state::ReactiveNetwork)
     state.u = u)
 end
 function save!(state::ReactiveNetwork)
-    return (push!(state.history_u, state.u); push!(state.history_t, state.t))
+    return (push!(state.history_u, copy(state.u)); push!(state.history_t, state.t))
 end
 
 function compile_observables(acs::ReactionNetwork)

tutorial/basics.jl

@@ -1,23 +1,32 @@
-using ReactiveDynamics
+using ReactiveDynamics, AlgebraicAgents
 
 # acs = @ReactionNetwork begin
 #     1.0, X ⟺ Y
 # end
 
 acs = @ReactionNetwork begin
-    1.0, X ⟺ Y, name => "transition1"
+    1.0, X --> Y, name => "transition1"
 end
 
 @prob_init acs X = 10 Y = 20
 @prob_params acs
 @prob_meta acs tspan = 250 dt = 0.1
 
-prob = @problematize acs
 
 # sol = ReactiveDynamics.solve(prob)
 
-sol = @solve prob
+#sol = @solve prob
+
+prob = @agentize acs
 
+for i in 1:30
+    AlgebraicAgents.step!(prob)
+end
+prob.history_u
 using Plots
 
 @plot sol plot_type = summary
+∪
+
+prob = @problematize acs
+@solve prob