From 8652f3d7940a740b021e2a295da111f1c10da954 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=ADma=2C=20Jan?= Date: Thu, 22 Feb 2024 19:48:46 +0100 Subject: [PATCH] Add structured agents and related examples --- Project.toml | 2 +- src/ReactiveDynamics.jl | 7 + src/compilers.jl | 33 +- src/interface/agents.jl | 52 +++ src/interface/create.jl | 2 + src/interface/reaction_parser.jl | 4 +- src/operators/equalize.jl | 5 +- src/operators/joins.jl | 28 +- src/solvers.jl | 197 +++++++++--- src/state.jl | 39 ++- tutorial/agents-integration/Project.toml | 9 + tutorial/agents-integration/agents.jl | 101 ++++++ .../agents-integration/agents_integration.jl | 302 ++++++++++++++++++ 13 files changed, 697 insertions(+), 84 deletions(-) create mode 100644 src/interface/agents.jl create mode 100644 tutorial/agents-integration/Project.toml create mode 100644 tutorial/agents-integration/agents.jl create mode 100644 tutorial/agents-integration/agents_integration.jl diff --git a/Project.toml b/Project.toml index 87db175..0ad533d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ReactiveDynamics" uuid = "c7456e7d-545a-4b79-91ea-6e93d96dd4d4" -version = "0.2.7" +version = "0.2.8" [deps] ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8" diff --git a/src/ReactiveDynamics.jl b/src/ReactiveDynamics.jl index 8a52ed1..a7192bf 100644 --- a/src/ReactiveDynamics.jl +++ b/src/ReactiveDynamics.jl @@ -37,6 +37,7 @@ TheoryReactionNetwork = BasicSchema( :ModalityAttributeT, :PcsOptT, :PrmAttributeT, + :BoolAttributeT, ], # AttrTypes [ # species @@ -47,6 +48,7 @@ TheoryReactionNetwork = BasicSchema( (:specCost, :S, :SampleableAttributeT), (:specReward, :S, :SampleableAttributeT), (:specValuation, :S, :SampleableAttributeT), + (:specStructured, :S, :BoolAttributeT), # transitions (:trans, :T, :SampleableAttributeT), (:transPriority, :T, :SampleableAttributeT), @@ -55,6 +57,7 @@ TheoryReactionNetwork = BasicSchema( (:transProbOfSuccess, :T, :SampleableAttributeT), (:transCapacity, :T, :SampleableAttributeT), (:transMaxLifeTime, :T, :SampleableAttributeT), + (:transPreAction, :T, :SampleableAttributeT), (:transPostAction, :T, :SampleableAttributeT), (:transMultiplier, :T, :SampleableAttributeT), (:transName, :T, :DescriptiveAttributeT), @@ -81,6 +84,7 @@ const ReactionNetworkSchema = FoldedReactionNetworkType{ Set{Symbol}, FoldedObservable, Any, + Bool } Base.convert(::Type{Symbol}, ex::String) = Symbol(ex) @@ -100,6 +104,7 @@ Base.convert(::Type{FoldedObservable}, ex::String) = eval(Meta.parse(ex)) prettynames = Dict( :transRate => [:rate], :specInitUncertainty => [:uncertainty, :stoch, :stochasticity], + :transPreAction => [:preAction, :pre], :transPostAction => [:postAction, :post], :transName => [:name, :interpretation], :transPriority => [:priority], @@ -117,6 +122,7 @@ defargs = Dict( :transCycleTime => 0.0, :transMaxLifeTime => Inf, :transMultiplier => 1, + :transPreAction => :(), :transPostAction => :(), :transName => missing, ), @@ -126,6 +132,7 @@ defargs = Dict( :specCost => 0.0, :specReward => 0.0, :specValuation => 0.0, + :specStructured => false, ), :P => Dict{Symbol,Any}(:prmVal => missing), :M => Dict{Symbol,Any}(:metaVal => missing), diff --git a/src/compilers.jl b/src/compilers.jl index 46cb842..0c741d8 100644 --- a/src/compilers.jl +++ b/src/compilers.jl @@ -26,16 +26,17 @@ end Recursively substitute model variables. Subsitution pairs are specified in `varmap`. """ function recursively_substitute_vars!(varmap, ex) - ex isa Symbol && return (haskey(varmap, ex) ? varmap[ex] : ex) - ex isa Expr && for i = 1:length(ex.args) - if ex.args[i] isa Expr - recursively_substitute_vars!(varmap, ex.args[i]) - else - ( - ex.args[i] isa Symbol && - haskey(varmap, ex.args[i]) && - (ex.args[i] = varmap[ex.args[i]]) - ) + if ex isa Symbol + return haskey(varmap, ex) ? varmap[ex] : ex + elseif ex isa Expr + for i = 1:length(ex.args) + if ex.args[i] isa Expr + ex.args[i] = recursively_substitute_vars!(varmap, ex.args[i]) + else + if ex.args[i] isa Symbol && haskey(varmap, ex.args[i]) + ex.args[i] = varmap[ex.args[i]] + end + end end end @@ -114,7 +115,12 @@ function wrap_expr(fex, species_names, prm_names, varmap) # the function shall be a function of the dynamic ReactionNetworkSchema structure: letex -> :(state -> $letex) # eval the expression to a Julia function, save that function into the "compiled" acset - return eval(:(state -> $letex)) + + return eval(quote + function (state, transition) + $letex + end + end) end function get_wrap_fun(acs::ReactionNetworkSchema) @@ -133,8 +139,9 @@ function skip_compile(attr) (string(attr) == "trans") end -function compile_attrs(acs::ReactionNetworkSchema) - species_names = collect(acs[:, :specName]) +function compile_attrs(acs::ReactionNetworkSchema, structured_species) + species_names = setdiff(collect(acs[:, :specName]), structured_species) + prm_names = collect(acs[:, :prmName]) varmap = Dict([name => :(state.u[$i]) for (i, name) in enumerate(species_names)]) for name in prm_names diff --git a/src/interface/agents.jl b/src/interface/agents.jl new file mode 100644 index 0000000..6233f63 --- /dev/null +++ b/src/interface/agents.jl @@ -0,0 +1,52 @@ +export AbstractStructuredSpecies, BaseStructuredSpecies +export @structured +export add_structured_species! + +# Abstract supertype of all structured species. +abstract type AbstractStructuredSpecies <: AbstractAlgebraicAgent end + +# It comes handy to keep track of the transition the entity is assigned to (if). +# In general, we will probably assume that each "structured agent" type implements this field. +# Otherwise, it would be possible to implement getter and setter interface and use it from within ReaDyn. +@aagent FreeAgent struct BaseStructuredSpecies + bound_transition::Union{Nothing, ReactiveDynamics.Transition} +end + +# We use this to let the network know that the type is structured. +function register_structured_species!(reaction_network, type) + if !(type ∈ reaction_network[:, :specName]) + add_part!(reaction_network, :S; specName = type) + end + + i = first(incident(reaction_network, type, :specName)) + reaction_network[i, :specStructured] = true + + return nothing +end + +# Convenience macro to define structured species. +macro structured(network, type) + name = Docs.namify(type.args[2]) + + quote + $(AlgebraicAgents.aagent(BaseStructuredSpecies, AbstractStructuredSpecies, type, ReactiveDynamics)) + register_structured_species!($(esc(network)), $(QuoteNode(name))) + end +end + +# Add a structured agent instance to an instance of a reaction network. +function add_structured_species!(problem::ReactionNetworkProblem, agent) + entangle!(getagent(problem, "structured/$(nameof(typeof(agent)))"), agent) +end + +import AlgebraicAgents + +# By default, structured agents have no evolutionary rule. +AlgebraicAgents._projected_to(::AbstractStructuredSpecies) = nothing +AlgebraicAgents._step!(::AbstractStructuredSpecies) = nothing + +# Tell if an agent is assigned to a transition, as a resource. +isblocked(a) = !isnothing(a.bound_transition) + +# Priority with which an unbound agent will be assigned to a transition. +priority(a, transition) = 0.0 \ No newline at end of file diff --git a/src/interface/create.jl b/src/interface/create.jl index ec49f85..099ca1e 100644 --- a/src/interface/create.jl +++ b/src/interface/create.jl @@ -295,6 +295,8 @@ function recursively_find_reactants!(reactants, pcs, ex::SampleableValues) end elseif isexpr(ex, :macrocall) recursively_find_reactants!(reactants, pcs, ex.args[3]) + elseif isexpr(ex, :call) + push!(reactants, ex.args[1]) else push!(reactants, underscorize(ex)) end diff --git a/src/interface/reaction_parser.jl b/src/interface/reaction_parser.jl index 2c566bb..a2e9ea2 100644 --- a/src/interface/reaction_parser.jl +++ b/src/interface/reaction_parser.jl @@ -3,7 +3,7 @@ using MacroTools: postwalk struct FoldedReactant - species::Symbol + species::Union{Expr, Symbol} stoich::SampleableValues modality::Set{Symbol} end @@ -71,6 +71,8 @@ function recursive_find_reactants!( 4:length(ex.args), ) recursive_find_reactants!(ex.args[3], mult, mods, reactants) + elseif isexpr(ex, :call) + push!(reactants, FoldedReactant(ex, mult, mods)) else @error("malformed reaction") end diff --git a/src/operators/equalize.jl b/src/operators/equalize.jl index c8e5724..f52672f 100644 --- a/src/operators/equalize.jl +++ b/src/operators/equalize.jl @@ -1,4 +1,4 @@ -export @equalize +export equalize!, @equalize expand_name_ff(ex) = if ex isa Expr && isexpr(ex, :macrocall) @@ -53,9 +53,10 @@ function equalize!(acs::ReactionNetworkSchema, eqs = []) for attr in propertynames(acs.subparts) attr == :specName && continue attr_ = acs[:, attr] - for i = 1:length(attr_) + for i in eachindex(attr_) attr_[i] = escape_ref(attr_[i], collect(keys(specmap))) attr_[i] = recursively_substitute_vars!(specmap, attr_[i]) + acs[i, attr] = attr_[i] end end diff --git a/src/operators/joins.jl b/src/operators/joins.jl index 5ac9808..c0183cb 100644 --- a/src/operators/joins.jl +++ b/src/operators/joins.jl @@ -1,24 +1,24 @@ # model joins -export @join +export union_acs!, @join using MacroTools using MacroTools: prewalk """ -Merge `acs2` onto `acs1`, the attributes in `acs2` taking precedence. Identify respective species given `eqs`, renaming species in `acs2`. +Merge `acs2` into `acs1`, the attributes in `acs2` taking precedence. Identify respective species given `eqs`, renaming species in `acs2`. """ -function union_acs!(acs1, acs2, name = gensym("acs_"), eqs = []) +function union_acs!(acs1, acs2, name = gensym("acs"), eqs = []) acs2 = deepcopy(acs2) prepend!(acs2, name, eqs) + for i in parts(acs2, :S) inc = incident(acs1, acs2[i, :specName], :specName) - isempty(inc) && (inc = add_part!(acs1, :S; specName = acs2[i, :specName]); - assign_defaults!(acs1)) - return (acs1, acs2) - println(first(inc)) - println(acs1[first(inc), :specModality]) - println() - println(acs2[:, :specModality]) + + if isempty(inc) + inc = add_part!(acs1, :S; specName = acs2[i, :specName]) + assign_defaults!(acs1) + end + union!(acs1[first(inc), :specModality], acs2[i, :specModality]) for attr in propertynames(acs1.subparts) @@ -30,7 +30,9 @@ function union_acs!(acs1, acs2, name = gensym("acs_"), eqs = []) new_trans_ix = add_parts!(acs1, :T, nparts(acs2, :T)) for attr in propertynames(acs2.subparts) !occursin("trans", string(attr)) && continue - acs1[new_trans_ix, attr] .= acs2[:, attr] + for (ix1, ix2) in enumerate(new_trans_ix) + acs1[ix2, attr] = acs2[ix1, attr] + end end foreach( @@ -69,10 +71,10 @@ function prepend!(acs::ReactionNetworkSchema, name = gensym("acs"), eqs = []) for attr in propertynames(acs.subparts) attr == :specName && continue attr_ = acs[:, attr] - for i = 1:length(attr_) + for i in eachindex(attr_) attr_[i] = escape_ref(attr_[i], collect(keys(specmap))) attr_[i] = recursively_substitute_vars!(specmap, attr_[i]) - attr_[i] isa Expr && (attr_[i] = prepend_obs(attr_[i], name)) + acs[i, attr] = attr_[i] end end diff --git a/src/solvers.jl b/src/solvers.jl index 270e8c7..889ddd3 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -15,7 +15,7 @@ Compute resource requirements given transition quantities. """ function get_reqs_init!(reqs, qs, state) reqs .= 0.0 - for i = 1:size(reqs, 2) + for i in axes(reqs, 2) for tok in state[i, :transLHS] !any(m -> m in tok.modality, [:rate, :nonblock]) && (reqs[tok.index, i] += qs[i] * tok.stoich) @@ -30,11 +30,14 @@ Compute resource requirements given transition quantities. """ function get_reqs_ongoing!(reqs, qs, state) reqs .= 0.0 - for i = 1:length(state.ongoing_transitions) + for i in eachindex(state.ongoing_transitions) 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.dt) + if in(:rate, tok.modality) && in(tok.species, state.structured_species) + error("Modality `:rate` is not supported for structured species in transition $(trans[:transName]).") + end in(:nonblock, tok.modality) && (reqs[tok.index, i] += qs[i] * tok.stoich) end end @@ -55,7 +58,7 @@ end function alloc_weighted!(reqs, u, priorities, state) allocs = zero(reqs) - for i = 1:size(reqs, 1) + for i in axes(reqs, 1) s = sum(reqs[i, :]) u[i] >= s && (allocs[i, :] .= reqs[i, :]; continue) foreach(j -> allocs[i, j] = reqs[i, j] * priorities[j], 1:size(reqs, 2)) @@ -69,7 +72,7 @@ end function alloc_greedy!(reqs, u, priorities, state) allocs = zero(reqs) sorted_trans = sort(1:size(reqs, 2); by = i -> -priorities[i]) - for i = 1:size(reqs, 1) + for i in axes(reqs, 1) s = sum(reqs[i, :]) u[i] >= s && (allocs[i, :] .= reqs[i, :]; continue) a = u[i] @@ -97,12 +100,14 @@ function get_frac_satisfied(allocs, reqs, state) return qs end +isinteger(x::Number) = x == trunc(x) + """ Given available allocations and qties of transitions requested to spawn, return number of spawned transitions. Update `alloc` to match actual allocation. """ function get_init_satisfied(allocs, qs, state) reqs = zero(allocs) - for i = 1:size(allocs, 2) + for i in axes(allocs, 2) all(allocs[:, i] .>= 0) || (allocs[:, i] .= 0.0; qs[i] = 0) for tok in state[i, :transLHS] !any(m -> m in tok.modality, [:rate, :nonblock]) && @@ -120,6 +125,7 @@ function get_init_satisfied(allocs, qs, state) return qs end + """ Evolve transitions, spawn new transitions. """ @@ -165,16 +171,43 @@ function evolve!(state) # add spawned transitions to the heap for i in parts(state, :T) - qs[i] != 0 && push!( - state.ongoing_transitions, - Transition( + if qs[i] != 0 + transition = Transition( string(state[i, :transName]) * "_@$(state.t)", + i, get_sampled_transition(state, i), + AbstractAlgebraicAgent[], + AbstractAlgebraicAgent[], state.t, qs[i], 0.0, - ), - ) + ) + push!(state.ongoing_transitions, transition) + + bound = transition.bound_structured_agents + + for (j, type) in enumerate(state.acs[:, :specName]) + if type ∈ state.structured_species + if !isinteger(allocs[j, i]) + error("For structured species, stoichiometry coefficient must be integer in transition $i.") + end + + all_of_type = collect(values(inners(getagent(state, "structured/$(string(type))")))) + filter!(!isblocked, all_of_type) + sort!(all_of_type, by=a->priority(a, state.acs[i, :transName]), rev=true) + + ix = 1 + while allocs[j, i] > 0 && ix <= length(all_of_type) + all_of_type[ix].bound_transition = transition + push!(bound, all_of_type[ix]) + allocs[j, i] -= 1 + ix += 1 + end + end + end + + context_eval(state, transition, state.wrap_fun(state.acs[i, :transPreAction])) + end end ## evolve ongoing transitions @@ -204,10 +237,33 @@ function evolve!(state) state.u .-= sum(allocs; dims = 2) actual_allocs .+= sum(allocs; dims = 2) - foreach( - i -> state.ongoing_transitions[i].state += qs[i] * state.dt, - eachindex(state.ongoing_transitions), - ) + for i in eachindex(state.ongoing_transitions) + transition = state.ongoing_transitions[i] + if qs[i] != 0 + transition.state += qs[i] * state.dt + bound = transition.nonblock_structured_agents + + for (j, type) in enumerate(state.acs[:, :specName]) + if type ∈ state.structured_species + if !isinteger(allocs[j, i]) + error("For structured species, stoichiometry coefficient must be integer in transition $i.") + end + + all_of_type = collect(values(inners(getagent(state, "structured/$(string(type))")))) + filter!(!isblocked, all_of_type) + sort!(all_of_type, by=a -> priority(a, state.acs[i, :transName]), rev=true) + + ix = 1 + while allocs[j, i] > 0 && ix <= length(all_of_type) + all_of_type[ix].bound_transition = transition + push!(bound, all_of_type[ix]) + allocs[j, i] -= 1 + ix += 1 + end + end + end + end + end push!(state.log, (:allocation, state.t, actual_allocs)) return push!( @@ -243,45 +299,74 @@ function finish!(state) while ix <= length(state.ongoing_transitions) trans_ = state.ongoing_transitions[ix] ((state.t - trans_.t) < trans_.trans[:transMaxLifeTime]) && - trans_.state < trans_[:transCycleTime] && + (trans_.state < trans_[:transCycleTime]) && (ix += 1; continue) - toks_rhs = [] + + q = if trans_.state >= trans_[:transCycleTime] + rand(Distributions.Binomial(Int(trans_.q), trans_[:transProbOfSuccess])) + else + 0 + end + for r in extract_reactants(trans_[:transRHS], state) - i = find_index(r.species, state) - push!( - toks_rhs, - UnfoldedReactant( - i, - r.species, - context_eval(state, state.wrap_fun(r.stoich)), - r.modality ∪ state[i, :specModality], - ), - ) + if r.species isa Expr + if !(r.species.args[1] ∈ state.structured_species) + error("Unresolved structured species species $(r.species.args[1]).") + end + + i = find_index(r.species.args[1], state) + stoich = context_eval(state, trans_, state.wrap_fun(r.stoich)) + + state.u[i] += q * stoich + val_reward += state[i, :specReward] * q * stoich + + for _ in 1:q + a = context_eval(state, trans_, state.wrap_fun(r.species)) + entangle!(getagent(state, "structured/$(r.species.args[1])"), a) + end + else + i = find_index(r.species, state) + stoich = context_eval(state, trans_, state.wrap_fun(r.stoich)) + + state.u[i] += q * stoich + val_reward += state[i, :specReward] * q * stoich + end end + for tok in trans_[:transLHS] - in(:conserved, tok.modality) && ( + if in(:conserved, tok.modality) state.u[tok.index] += trans_.q * tok.stoich * (in(:rate, tok.modality) ? trans_[:transCycleTime] : 1) - ) - end - - q = if trans_.state >= trans_[:transCycleTime] - rand(Distributions.Binomial(Int(trans_.q), trans_[:transProbOfSuccess])) - else - 0 + if tok.species ∈ state.structured_species + for _ in 1:(trans_.q * tok.stoich) + trans_.bound_structured_agents[begin].bound_transition = nothing + deleteat!(trans_.bound_structured_agents, 1) + end + end + end + + if in(:nonblock, tok.modality) + if in(:conserved, tok.modality) + error("Modalities `:conserved` and `:nonblock` cannot be specified at the same time.") + end + + state.u[tok.index] += trans_.q * tok.stoich + if tok.species ∈ state.structured_species + for _ in 1:(trans_.q * tok.stoich) + trans_.nonblock_structured_agents[begin].bound_transition = nothing + deleteat!(trans_.nonblock_structured_agents, 1) + end + end + end end - foreach( - tok -> (state.u[tok.index] += q * tok.stoich; - val_reward += state[tok.index, :specReward] * q * tok.stoich), - toks_rhs, - ) + context_eval(state, trans_, state.wrap_fun(state.acs[trans_.i, :transPostAction])) - context_eval(state, trans_.trans[:transPostAction]) terminated_all[Symbol(trans_[:transHash])] = get(terminated_all, Symbol(trans_[:transHash]), 0) + trans_.q + terminated_success[Symbol(trans_[:transHash])] = get(terminated_success, Symbol(trans_[:transHash]), 0) + q @@ -301,6 +386,14 @@ function free_blocked_species!(state) for trans in state.ongoing_transitions, tok in trans[:transLHS] in(:nonblock, tok.modality) && (state.u[tok.index] += q * tok.stoich) end + + for trans in state.ongoing_transitions + for a in trans.nonblock_structured_agents + a.bound_transition = nothing + end + + empty!(trans.nonblock_structured_agents) + end end ## resolve tspan, tstep @@ -319,7 +412,7 @@ function ReactionNetworkProblem( acs::ReactionNetworkSchema, u0 = Dict(), p = Dict(); - name = "reactive_network", + name = "reaction_network", kwargs..., ) assign_defaults!(acs) @@ -333,7 +426,10 @@ function ReactionNetworkProblem( keywords[:tspan], keywords[:tstep] = get_tcontrol(keywords[:tspan], keywords) acs = remove_choose(acs) - attrs, transitions, wrap_fun = compile_attrs(acs) + + structured_species_names = acs[filter(i -> acs[i, :specStructured], 1:nparts(acs, :S)), :specName] + + attrs, transitions, wrap_fun = compile_attrs(acs, structured_species_names) transition_recipes = transitions u0_init = zeros(nparts(acs, :S)) @@ -374,6 +470,7 @@ function ReactionNetworkProblem( u0_init, merge(p, Dict(:strategy => get(keywords, :alloc_strategy, :weighted))), keywords[:tspan][1], + Symbol[], keywords[:tspan], get(keywords, :tstep, 1), transitions, @@ -384,6 +481,14 @@ function ReactionNetworkProblem( sol, ) + entangle!(network, FreeAgent("structured")) + + structured_species = filter(i -> acs[i, :specStructured], 1:nparts(acs, :S)) + for i in structured_species + push!(network.structured_species, acs[i, :specName]) + entangle!(getagent(network, "structured"), FreeAgent(string(acs[i, :specName]))) + end + save!(network) return network @@ -400,12 +505,20 @@ function AlgebraicAgents._reinit!(state::ReactionNetworkProblem) return state end +function update_u_structured!(state) + for (i, type) in enumerate(state.acs[:, :specName]) + structured_agents_type = values(inners(getagent(state, "structured/$type"))) + state.u[i] = count(!isblocked, structured_agents_type) + end +end + function AlgebraicAgents._step!(state::ReactionNetworkProblem) free_blocked_species!(state) update_observables(state) sample_transitions!(state) evolve!(state) finish!(state) + event_action!(state) push!( diff --git a/src/state.jl b/src/state.jl index 97e85db..70ed436 100644 --- a/src/state.jl +++ b/src/state.jl @@ -12,14 +12,20 @@ end Ongoing transition auxiliary structure. """ @aagent struct Transition + i::Int + trans::Dict{Symbol,Any} + bound_structured_agents::Vector{AbstractAlgebraicAgent} + nonblock_structured_agents::Vector{AbstractAlgebraicAgent} + t::Float64 q::Float64 state::Float64 end Base.getindex(state::Transition, key) = state.trans[key] +Base.setindex!(state::Transition, val, key) = state.trans[key] = val @aagent struct Observable last::Float64 # last sampling time @@ -40,6 +46,8 @@ end p::Any t::Float64 + structured_species::Vector{Symbol} + tspan::Tuple{Float64,Float64} dt::Float64 @@ -55,17 +63,21 @@ end # get value of a numeric expression # evaluate compiled numeric expression in context of (u, p, t) -function context_eval(state::ReactionNetworkProblem, o) - o = o isa Function ? Base.invokelatest(o, state) : o - +function context_eval(state::ReactionNetworkProblem, transition, o) + o = o isa Function ? Base.invokelatest(o, state, transition) : o + return o isa Sampleable ? rand(o) : o end function Base.getindex(state::ReactionNetworkProblem, keys...) - return context_eval( - state, - (contains(string(keys[2]), "trans") ? state.transitions : state.attrs)[keys[2]][keys[1]], - ) + if any(occursin.(["transPreAction", "transPostAction"], Ref(string(keys[2])))) + return state.acs[keys[1], keys[2]] + else + return context_eval( + state, nothing, + (contains(string(keys[2]), "trans") ? state.transitions : state.attrs)[keys[2]][keys[1]], + ) + end end function init_u!(state::ReactionNetworkProblem) @@ -120,7 +132,7 @@ function resample!(state::ReactionNetworkProblem, o::Observable) o.last = state.t isempty(o.range) && (return o.val = missing) - return o.sampled = context_eval(state, sample_range(o.range, state)) + return o.sampled = context_eval(state, nothing, sample_range(o.range, state)) end resample(state::ReactionNetworkProblem, o::Symbol) = resample!(state, state.observables[o]) @@ -164,10 +176,10 @@ function sample_transitions!(state::ReactionNetworkProblem) l_line, r_line = prune_r_line(state.transition_recipes[:trans][i]) for attr in keys(state.transition_recipes) - attr ∈ [:trans, :transPostAction, :transActivated, :transHash] && continue + (attr ∈ [:trans, :transPreAction, :transPostAction, :transActivated, :transHash]) && continue push!( state.transitions[attr], - context_eval(state, state.transition_recipes[attr][i]), + context_eval(state, nothing, state.transition_recipes[attr][i]), ) end @@ -179,17 +191,20 @@ function sample_transitions!(state::ReactionNetworkProblem) UnfoldedReactant( j, r.species, - context_eval(state, state.wrap_fun(r.stoich)), + context_eval(state, nothing, state.wrap_fun(r.stoich)), r.modality ∪ state[j, :specModality], ), ) end + push!(state.transitions[:transLHS], reactants) push!(state.transitions[:transRHS], r_line) + foreach( k -> push!(state.transitions[k], state.transition_recipes[k][i]), - [:transPostAction, :transToSpawn, :transHash], + [:transPreAction, :transPostAction, :transToSpawn, :transHash], ) + state.transition_recipes[:transToSpawn] .= 0 end end diff --git a/tutorial/agents-integration/Project.toml b/tutorial/agents-integration/Project.toml new file mode 100644 index 0000000..ef794fc --- /dev/null +++ b/tutorial/agents-integration/Project.toml @@ -0,0 +1,9 @@ +[deps] +AlgebraicAgents = "f6eb0ae3-10fa-40e6-88dd-9006ba45093a" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" +MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +PlotGraphviz = "78a92bc3-407c-4e2f-aae5-75bb47a6fe36" +ReactiveDynamics = "c7456e7d-545a-4b79-91ea-6e93d96dd4d4" +SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" diff --git a/tutorial/agents-integration/agents.jl b/tutorial/agents-integration/agents.jl new file mode 100644 index 0000000..c4304b9 --- /dev/null +++ b/tutorial/agents-integration/agents.jl @@ -0,0 +1,101 @@ +# -------------------------------------------------------------------------------- +# Structured (agent-based) species + +#= +import Pkg +Pkg.activate(".") +Pkg.dev("../..") +Pkg.add(["AlgebraicAgents"]) +=# + +using ReactiveDynamics +using AlgebraicAgents + +# Define the "symbolic" reaction network. +network = @ReactionNetworkSchema + +# Below, we combine +# - "classical species" (continuous or discrete; considered as pure quantities); +# - "structured agents" (possibly with custom evolutionary function; these can appear both on LHS and RHS). + +@push network begin + # With specified intensities, generate experimental resources. + ρ1, ∅ --> R1 + ρ2, ∅ --> R2 + + # Generate "Molecule 1" (where the integer corresponds to a "state" of, e.g., experimental triage). + ρ3, ∅ --> M1(@t(), rand(4)) + + # Based on properties of particular "structured agent" assigned to the transition, + # we can update the attributes of the instance of a transition (such as probability of success). + + # Transition "Molecule 1" into "Molecule 2." + # Update transition probability based on properties of "M1," + # which was assigned as a "resource" to the transition. + ρ4, R1 + M1 --> M2(@t(), rand(4)), preAction => update_prob_transition(state, transition) +end + +@prob_init network R1 = 10 R2 = 15 + +# As for structured agents, we will need to instantiate the instances +# and add them to the instance of a network. But first, we still need to define these types. +@prob_init network M1 = 2 M2 = 0 + +@prob_params network ρ1 = 2 ρ2 = 1 ρ3 = 3 ρ4 = 4 + +@prob_meta network tspan = 100 dt = 1.0 + +# We use `@structured` macro, which is a convenience wrapper around `@aagent`), +# defined in ReactiveDynamics.jl +@structured network struct M1 + descriptor + time_created +end + +using Random + +# Type `M1` lives in the scope of ReactiveDynamics. +# Accordingly, we have to explicitly declare the scope. +using ReactiveDynamics: M1 + +ReactiveDynamics.M1(time, descriptor) = M1("M1" * randstring(4), nothing, descriptor, time) + +# We define the function which updates the transition probability. +# This has to be accessible from within the name scope of ReactiveDynamics. +@register begin + update_prob_transition = function (state, transition) + if !isnothing(transition) && !isempty(transition.bound_structured_agents) + bound_agent = first(transition.bound_structured_agents) + + transition[:transProbOfSuccess] = min(1.0, sum(bound_agent.descriptor)) + end + end +end + +# Alternatively, we can define a structured agent type using +# the usual `@aagent` macro. This must be evaluated inside the scope +# of ReactiveDynamics. +@register begin + @aagent BaseStructuredSpecies AbstractStructuredSpecies struct M2 + descriptor + time_created + end + + using Random + M2(time, descriptor) = M2("M2" * randstring(4), nothing, descriptor, time) +end + +# Let the network know that the species is structured. +ReactiveDynamics.register_structured_species!(network, :M2) + +# -------------------------------------------------------------------------------- +# Instantiate the network. +network_instance = ReactionNetworkProblem(network) + +for i in 1:2 + add_structured_species!(network_instance, ReactiveDynamics.M1(0.0, rand(4))) +end + +# -------------------------------------------------------------------------------- +# Simulate the network. +simulate(network_instance, 10) \ No newline at end of file diff --git a/tutorial/agents-integration/agents_integration.jl b/tutorial/agents-integration/agents_integration.jl new file mode 100644 index 0000000..25fb6fd --- /dev/null +++ b/tutorial/agents-integration/agents_integration.jl @@ -0,0 +1,302 @@ +import Pkg + +# Manually add dependencies +#= +Pkg.activate(".") +Pkg.develop(path = "../..") +Pkg.add(["MacroTools", "AlgebraicAgents", "JSON3", "Distributions"]) +=# + +# -------------------------------------------------------------------------------- +# Parametrize the reaction network + +# We will have two classes of entities, each of which can be in a number of different states. +# States "S" and "F" will internally denote some terminal, "success" and "failure" states, respectively. + +terminal_states = ["S", "F"] +M1_states = "M1_" .* ["A", "B", "C", "D", terminal_states...] +M2_states = "M2_" .* ["A", "B", "C", "D", "E", terminal_states...] + +experimental_resources = ["ER1", "ER2", "ER3"] + +M1_transition_probs = round.(rand(length(M1_states), length(M1_states)); digits=2) +M2_transition_probs = round.(rand(length(M2_states), length(M2_states)); digits=2) + +# Ensure forward flow between states and set zero transition prob to terminal states +for probs in [M1_transition_probs, M2_transition_probs] + for i in axes(probs, 1) + probs[i, i:end] .= 0 + end + + probs[:, end-1] .= 0 +end + +M1_priorities = round.(rand(length(M1_states)); digits=2) +M2_priorities = 2 * round.(rand(length(M2_states)); digits=2) + +M1_resources = [rand(1:5, length(experimental_resources)) for _ in M1_states] +M2_resources = [2 * rand(1:5, length(experimental_resources)) for _ in M2_states] + +M1_durations = rand(2:5, length(M1_states)) +M2_durations = rand(2:5, length(M2_states)) + +# -------------------------------------------------------------------------------- +# Add initial quantities + +M1_initial = [rand(1:10, length(M1_states)-2)..., 0, 0] +M2_initial = [rand(1:10, length(M2_states)-2)..., 0, 0] + +experimental_resources_initial = rand(1e2:3e2, size(experimental_resources)) + +# -------------------------------------------------------------------------------- +# Export to JSON + +# First, we build a dictionary. +data = Dict("species" => [], "transitions" => []) + +species, transitions = data["species"], data["transitions"] + +# -------------------------------------------------------------------------------- +# Add species with initial quantities. +# Later, we may add a couple more attributes (conserved resources, etc.). + +for (state, q) in zip(M1_states, M1_initial) + push!(species, Dict("name" => "$state", "initial" => q)) +end + +for (state, q) in zip(M2_states, M2_initial) + push!(species, Dict("name" => "$state", "initial" => q)) +end + +for (res, q) in zip(experimental_resources, experimental_resources_initial) + push!(species, Dict("name" => "$res", "initial" => q)) +end + +# -------------------------------------------------------------------------------- +# Add transitions for entity class M1. + +for (i, state) in enumerate(M1_states[begin:end-2]) + t = Dict{String, Any}("priority" => M1_priorities[i], "duration" => M1_durations[i], "rate" => state, "name" => "M1") + + # from (left-hand side) + push!(t, "from" => [Dict("name" => res, "q" => q) for (res, q) in zip(experimental_resources, M1_resources[i])]) + push!(t["from"], Dict("name" => state, "q" => 1)) + + # to (right-hand side) + push!(t, "to" => []) + to = t["to"] + + for (j, state2) in enumerate(M1_states[i:end]) + if M1_transition_probs[j+i-1, i] > 0 + push!(to, Dict("probability" => M1_transition_probs[j+i-1, i], "reactants" => [Dict("name" => state2, "q" => 1)])) + end + end + + push!(transitions, t) +end + +# Add transitions for entity class M2. + +for (i, state) in enumerate(M2_states[begin:end-2]) + t = Dict{String, Any}("priority" => M2_priorities[i], "duration" => M2_durations[i], "rate" => state, "name" => "M2") + + # from (left-hand side) + push!(t, "from" => [Dict("name" => res, "q" => q) for (res, q) in zip(experimental_resources, M2_resources[i])]) + push!(t["from"], Dict("name" => state, "q" => 1)) + + # to (right-hand side) + push!(t, "to" => []) + to = t["to"] + + for (j, state2) in enumerate(M2_states[i:end]) + if M2_transition_probs[j+i-1, i] > 0 + push!(to, Dict("probability" => M2_transition_probs[j+i-1, i], "reactants" => [Dict("name" => state2, "q" => 1)])) + end + end + + push!(transitions, t) +end + +using JSON3 +open("reaction_network.json", "w") do io + #JSON3.pretty(io, data) +end + +# -------------------------------------------------------------------------------- +# Import from JSON +# Eventually move into the module (refactor current CSV interface). + +# Species and initial values. +name = "reaction_network" +str_init = "@prob_init $name" +for s in data["species"] + str_init *= " $(s["name"])=$(s["initial"])" +end + +function get_subline(d::Vector) + if isempty(d) + return "∅" + elseif haskey(first(d), "probability") + sublines = ["($(sd["probability"]), " * join(["$(s["q"]) * $(s["name"])" for s in sd["reactants"]], " + ") * ")" for sd in d] + return "@choose(" * join(sublines, ", ") * ")" + else + return join(["$(s["q"]) * $(s["name"])" for s in d], " + ") + end +end + +# Transitions. +str_transitions = [] +for t in data["transitions"] + line = t["rate"] * ", " * get_subline(t["from"]) * " --> " * get_subline(t["to"]) + line *= ", name => $(t["name"]), priority => $(t["priority"]), cycletime => $(t["duration"])" + push!(str_transitions, line) +end + +push!(str_transitions, "i1, ∅ --> M1_A, name => M1_creation") +push!(str_transitions, "i2, ∅ --> M2_A, name => M2_creation") + +str_params = "@prob_params $name i1 = .3 i2 = .2" + +str_network_def = """ + begin + $name = @ReactionNetworkSchema + @push $name begin + $(join(str_transitions, '\n')) + end + $str_init + $str_params + end +""" + +using MacroTools: striplines +expr_network_def = striplines(Meta.parseall(str_network_def)) + +using ReactiveDynamics + +eval(expr_network_def) + +@isdefined reaction_network + +# -------------------------------------------------------------------------------- +# Solve problem + +@prob_meta reaction_network tspan = 100 dt = 1.0 + +# Convert network into an AlgAgents hierarchy. +problem = ReactionNetworkProblem(reaction_network; name = "network") + +# AlgAgents: "periodic" callback. +using AlgebraicAgents + +@aagent struct Controller + M1_λ::Float64 + M2_λ::Float64 + + log::Vector{String} + + current_time::Float64 + time_step::Float64 + + initial_time::Float64 +end + +function Controller(name::String, M1_λ::Float64, M2_λ::Float64, time::T, time_step::T) where {T<:Real} + return Controller(name, M1_λ, M2_λ, String[], time, time_step, time) +end + +using Distributions: Poisson + +function AlgebraicAgents._step!(c::Controller) + n_removed_M1 = rand(Poisson(c.time_step * c.M1_λ)) + n_removed_M2 = rand(Poisson(c.time_step * (c.M2_λ + c.M1_λ))) + + transitions = getagent(c, "../network").ongoing_transitions + M1_transitions = filter(x -> x[:transName] === :M1, transitions) + M2_transitions = filter(x -> x[:transName] === :M2, transitions) + + M1_transitions_delete = isempty(M1_transitions) ? [] : unique(rand(M1_transitions, n_removed_M1)) + for trans in M1_transitions_delete + trans.state = trans[:transCycleTime] + trans.trans[:transProbOfSuccess] = 0 + end + + M2_transitions_delete = isempty(M2_transitions) ? [] : unique(rand(M2_transitions, n_removed_M2)) + for trans in M2_transitions_delete + trans.state = trans[:transCycleTime] + trans.trans[:transProbOfSuccess] = 0 + end + + push!(c.log, "t = $(c.current_time) removed compounds: " * join(getname.(union(M1_transitions_delete, M2_transitions_delete)), ", ")) + + return c.current_time += c.time_step +end + +AlgebraicAgents._projected_to(c::Controller) = c.current_time + +c = Controller("controller", 1e-1, 2e-1, 0., 1.); + +compound_problem = ⊕(problem, c; name = "compound problem"); + +# Simulate +simulate(compound_problem, 100); + +# Access solution +compound_problem.inners["network"].sol +compound_problem.inners["controller"].log + +# Plot solution +draw(compound_problem.inners["network"]) + +# -------------------------------------------------------------------------------- +# Add resource making part to the reaction network + +eval(expr_network_def) + +n_primary_resources = 5 +primary_resources = ["R$i" for i in 1:n_primary_resources] + +str_resource_making_transitions = [ + "p_primary_$i, ∅ --> $res, name => primary_resource_maker_$i" for (i, res) in enumerate(primary_resources) +] + +for (i, res) in enumerate(experimental_resources) + stoich = rand(1:5, n_primary_resources) + rate = "p_$i * (" * join(primary_resources, " + ") * ")" + transition_lhs = join(["$(stoich[i])" * "$primary_res" for (i, primary_res) in enumerate(primary_resources)], " + ") + + push!(str_resource_making_transitions, "$rate, $transition_lhs --> $res, name => resource_maker_$i") +end + +str_resource_making_params = + "@prob_params resource_making_network " * + join(["p_$i = $(rand(1:3))" for i in 1:length(experimental_resources)], " ") * " " * + join(["p_primary_$i = $(rand(1:5))" for i in 1:length(primary_resources)], " ") + +str_network_def = """ + begin + resource_making_network = @ReactionNetworkSchema + @push resource_making_network begin + $(join(str_resource_making_transitions, '\n')) + end + $str_resource_making_params + end +""" + +expr_network_resource_making_def = striplines(Meta.parseall(str_network_def)) + +eval(expr_network_resource_making_def) + +extended_network = union_acs!(reaction_network, resource_making_network, "resource_making_network") + +#equalize!(extended_network, [Meta.parseall("$res=resource_making_network.$res") for res in experimental_resources]) +str_equalize = "@equalize extended_network " * join(["$res=resource_making_network.$res" for res in experimental_resources], " ") + +equalize_expr = striplines(Meta.parseall(str_equalize)) +eval(equalize_expr) + +@prob_meta extended_network tspan = 100 dt = 1.0 + +# Convert network into an AlgAgents hierarchy. +extended_problem = ReactionNetworkProblem(extended_network; name = "extended_network") + +simulate(extended_problem, 100); \ No newline at end of file