From 4791a81b95997ee88a8cfc9142103d065f775e6e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=ADma=2C=20Jan?= Date: Tue, 12 Mar 2024 14:45:31 +0100 Subject: [PATCH] Structured species - An agent now implements a property pointing at the species. - Using `@structured` macro, it is possible to set call a custom agent constructor and assign the species. - Using `@move` macro, it is possible to reuse LHS agents on the RHS and modify their species. --- src/compilers.jl | 14 ++- src/interface/agents.jl | 48 +++++--- src/interface/reaction_parser.jl | 5 +- src/solvers.jl | 170 +++++++++++++++++++------- src/state.jl | 3 +- tutorial/agents-integration/agents.jl | 29 +++-- 6 files changed, 195 insertions(+), 74 deletions(-) diff --git a/src/compilers.jl b/src/compilers.jl index b4faba8..eb651ba 100644 --- a/src/compilers.jl +++ b/src/compilers.jl @@ -66,7 +66,6 @@ end reserved_names = [:t, :state, :obs, :resample, :solverarg, :take, :log, :periodic, :set_params] -push!(reserved_names, :state) function escape_ref(ex, species) return if ex isa Symbol @@ -103,6 +102,15 @@ function wrap_expr(fex, species_names, prm_names, varmap) end end + fex = prewalk(fex) do x + # here we convert the query metalanguage: @t() -> time(state) etc. + if isexpr(x, :macrocall) && (macroname(x) == :transition) + :transition + else + x + end + end + # substitute the species names with "pointers" into the state space: S -> state.u[1] fex = recursively_substitute_vars!(varmap, fex) # substitute the params names with "pointers" into the parameter space: β -> state.p[:β] @@ -139,8 +147,8 @@ function skip_compile(attr) (string(attr) == "trans") end -function compile_attrs(acs::ReactionNetworkSchema, structured_species) - species_names = setdiff(collect(acs[:, :specName]), structured_species) +function compile_attrs(acs::ReactionNetworkSchema, structured_token) + species_names = setdiff(collect(acs[:, :specName]), structured_token) prm_names = collect(acs[:, :prmName]) varmap = Dict([name => :(state.u[$i]) for (i, name) in enumerate(species_names)]) diff --git a/src/interface/agents.jl b/src/interface/agents.jl index 842a7f8..fcd5558 100644 --- a/src/interface/agents.jl +++ b/src/interface/agents.jl @@ -1,15 +1,17 @@ -export AbstractStructuredSpecies, BaseStructuredSpecies -export @structured -export add_structured_species! +export AbstractStructuredToken, BaseStructuredToken +export @structured_token +export add_structured_token! # Abstract supertype of all structured species. -abstract type AbstractStructuredSpecies <: AbstractAlgebraicAgent end +abstract type AbstractStructuredToken <: 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 +@aagent FreeAgent struct BaseStructuredToken + species::Union{Nothing,Symbol} bound_transition::Union{Nothing,ReactiveDynamics.Transition} + past_bonds::Vector{Tuple{Symbol,Float64,Transition}} end # We use this to let the network know that the type is structured. @@ -25,33 +27,45 @@ function register_structured_species!(reaction_network, type) end # Convenience macro to define structured species. -macro structured(network, type) - name = Docs.namify(type.args[2]) - +macro structured_token(network, type) quote $(AlgebraicAgents.aagent( - BaseStructuredSpecies, - AbstractStructuredSpecies, + BaseStructuredToken, + AbstractStructuredToken, 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) - return entangle!(getagent(problem, "structured/$(nameof(typeof(agent)))"), agent) +function add_structured_token!(problem::ReactionNetworkProblem, agent) + return entangle!(getagent(problem, "structured"), agent) end import AlgebraicAgents # By default, structured agents have no evolutionary rule. -AlgebraicAgents._projected_to(::AbstractStructuredSpecies) = nothing -AlgebraicAgents._step!(::AbstractStructuredSpecies) = nothing +AlgebraicAgents._projected_to(::AbstractStructuredToken) = nothing +AlgebraicAgents._step!(::AbstractStructuredToken) = nothing # Tell if an agent is assigned to a transition, as a resource. -isblocked(a) = !isnothing(a.bound_transition) +isblocked(a::AbstractStructuredToken) = !isnothing(get_bound_transition(a)) + +# Add a record that an agent was used as "species" in a "transition". +function add_to_log!(a::AbstractStructuredToken, species::Symbol, t, transition::Transition) + return push!(a.past_bonds, (species, Float64(t), transition)) +end + +# Set the transition a token is bound to. +get_bound_transition(a::AbstractStructuredToken) = a.bound_transition +function set_bound_transition!(a::AbstractStructuredToken, t::Union{Nothing,Transition}) + return a.bound_transition = t +end # Priority with which an unbound agent will be assigned to a transition. -priority(a, transition) = 0.0 +priority(a::AbstractStructuredToken, transition) = 0.0 + +# What species (place) is an agent currently assigned to. +get_species(a::AbstractStructuredToken) = a.species +set_species!(a::AbstractStructuredToken, species::Symbol) = a.species = species diff --git a/src/interface/reaction_parser.jl b/src/interface/reaction_parser.jl index 8ad0651..1f83113 100644 --- a/src/interface/reaction_parser.jl +++ b/src/interface/reaction_parser.jl @@ -63,6 +63,9 @@ function recursive_find_reactants!( for i = 2:length(ex.args) recursive_find_reactants!(ex.args[i], mult, mods, reactants) end + elseif isexpr(ex, :call) || + (ex.head == :macrocall && macroname(ex) ∈ [:structured, :move]) + push!(reactants, FoldedReactant(ex, mult, mods)) elseif ex.head == :macrocall mods = copy(mods) macroname(ex) in species_modalities && push!(mods, macroname(ex)) @@ -71,8 +74,6 @@ 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/solvers.jl b/src/solvers.jl index cb86797..9bf982c 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -35,7 +35,7 @@ function get_reqs_ongoing!(reqs, qs, state) 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) + if in(:rate, tok.modality) && in(tok.species, state.structured_token) error( "Modality `:rate` is not supported for structured species in transition $(trans[:transName]).", ) @@ -170,6 +170,8 @@ function evolve!(state) state.u .-= sum(allocs; dims = 2) actual_allocs .+= sum(allocs; dims = 2) + structured_token = collect(values(inners(getagent(state, "structured")))) + # add spawned transitions to the heap for i in parts(state, :T) if qs[i] != 0 @@ -179,6 +181,7 @@ function evolve!(state) get_sampled_transition(state, i), AbstractAlgebraicAgent[], AbstractAlgebraicAgent[], + [], state.t, qs[i], 0.0, @@ -186,29 +189,35 @@ function evolve!(state) push!(state.ongoing_transitions, transition) bound = transition.bound_structured_agents + structured_to_agents = transition.structured_to_agents for (j, type) in enumerate(state.acs[:, :specName]) - if type ∈ state.structured_species + if type ∈ state.structured_token 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))"))), + available_species = filter( + a -> get_species(a) == type && !isblocked(a), + structured_token, ) - filter!(!isblocked, all_of_type) + sort!( - all_of_type; + available_species; 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]) + while allocs[j, i] > 0 && ix <= length(available_species) + set_bound_transition!(available_species[ix], transition) + + push!(bound, available_species[ix]) + push!(structured_to_agents, type => available_species[ix]) + add_to_log!(available_species[ix], type, state.t, transition) + allocs[j, i] -= 1 ix += 1 end @@ -250,30 +259,40 @@ function evolve!(state) transition = state.ongoing_transitions[i] if qs[i] != 0 transition.state += qs[i] * state.dt + bound = transition.nonblock_structured_agents + structured_to_agents = transition.structured_to_agents for (j, type) in enumerate(state.acs[:, :specName]) - if type ∈ state.structured_species + if type ∈ state.structured_token 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))"))), + available_species = filter( + a -> get_species(a) == type && !isblocked(a), + structured_token, ) - filter!(!isblocked, all_of_type) + sort!( - all_of_type; + available_species; 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]) + while allocs[j, i] > 0 && ix <= length(available_species) + set_bound_transition!( + available_species[ix].bound_transition, + transition, + ) + + push!(bound, available_species[ix]) + push!(structured_to_agents, type => available_species[ix]) + add_to_log!(available_species[ix], type, state.t, transition) + allocs[j, i] -= 1 ix += 1 end @@ -306,6 +325,64 @@ function event_action!(state) end end +function allocate_for_move(t::Transition, s::Symbol) + return t.bound_structured_agents ∩ + map(x -> x[2], filter(x -> x[1] == s, t.structured_to_agents)) +end + +function structured_rhs(expr::Expr, state, transition) + if isexpr(expr, :macrocall) && macroname(expr) == :structured + expr = quote + token = $(expr.args[end-1]) + species = $(expr.args[end]) + + return token, species + end + + token, species = context_eval(state, transition, state.wrap_fun(expr)) + set_species!(token, Symbol(species)) + + entangle!(getagent(state, "structured"), token) + + return token, get_species(token) + elseif isexpr(expr, :macrocall) && macroname(expr) == :move + expr = quote + species_from = $(expr.args[end-1]) + species_to = $(expr.args[end]) + + return species_from, species_to + end + + species_from, species_to = + Symbol.(context_eval(state, transition, state.wrap_fun(expr))) + + tokens = + filter(x -> get_species(x) == species_from, transition.bound_structured_agents) + if !isempty(tokens) + token = first(tokens) + entangle!(getagent(state, "structured"), token) + + set_species!(token, species_to) + ix = findfirst( + i -> transition.bound_structured_agents[i] == token, + eachindex(transition.bound_structured_agents), + ) + deleteat!(transition.bound_structured_agents, ix) + set_bound_transition!(token, nothing) + + return token, species_to + else + @error "Not enough tokens to allocate for a move." + end + + else + token = context_eval(state, transition, state.wrap_fun(expr)) + entangle!(getagent(state, "structured"), token) + + return token, get_species(token) + end +end + # collect terminated transitions function finish!(state) val_reward = 0 @@ -327,19 +404,13 @@ function finish!(state) for r in extract_reactants(trans_[:transRHS], state) 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 _ = 1:q - a = context_eval(state, trans_, state.wrap_fun(r.species)) - entangle!(getagent(state, "structured/$(r.species.args[1])"), a) + for _ = 1:(q*stoich) + token, species = structured_rhs(r.species, state, trans_) + i = find_index(species, state) + state.u[i] += 1 + val_reward += state[i, :specReward] end else i = find_index(r.species, state) @@ -356,9 +427,13 @@ function finish!(state) trans_.q * tok.stoich * (in(:rate, tok.modality) ? trans_[:transCycleTime] : 1) - if tok.species ∈ state.structured_species + if tok.species ∈ state.structured_token for _ = 1:(trans_.q*tok.stoich) - trans_.bound_structured_agents[begin].bound_transition = nothing + isempty(trans_.bound_structured_agents) && break + set_bound_transition!( + trans_.bound_structured_agents[begin].bound_transition, + nothing, + ) deleteat!(trans_.bound_structured_agents, 1) end end @@ -372,9 +447,12 @@ function finish!(state) end state.u[tok.index] += trans_.q * tok.stoich - if tok.species ∈ state.structured_species + if tok.species ∈ state.structured_token for _ = 1:(trans_.q*tok.stoich) - trans_.nonblock_structured_agents[begin].bound_transition = nothing + set_bound_transition!( + trans_.nonblock_structured_agents[begin].bound_transition, + nothing, + ) deleteat!(trans_.nonblock_structured_agents, 1) end end @@ -439,6 +517,7 @@ function ReactionNetworkProblem( acs[i, :metaKeyword] => acs[i, :metaVal] for i in parts(acs, :M) if !isnothing(acs[i, :metaKeyword]) && !isnothing(acs[i, :metaVal]) ]) + merge!(keywords, Dict(collect(kwargs))) merge!(keywords, Dict(:strategy => get(keywords, :alloc_strategy, :weighted))) @@ -446,10 +525,10 @@ function ReactionNetworkProblem( acs = remove_choose(acs) - structured_species_names = + structured_token_names = acs[filter(i -> acs[i, :specStructured], 1:nparts(acs, :S)), :specName] - attrs, transitions, wrap_fun = compile_attrs(acs, structured_species_names) + attrs, transitions, wrap_fun = compile_attrs(acs, structured_token_names) transition_recipes = transitions u0_init = zeros(nparts(acs, :S)) @@ -482,6 +561,7 @@ function ReactionNetworkProblem( "t" => Float64[], (string(name) => Float64[] for name in acs[:, :specName])..., ) + network = ReactionNetworkProblem( name, acs, @@ -490,7 +570,7 @@ function ReactionNetworkProblem( u0_init, merge(p, Dict(:strategy => get(keywords, :alloc_strategy, :weighted))), keywords[:tspan][1], - Symbol[], + structured_token_names, keywords[:tspan], get(keywords, :tstep, 1), transitions, @@ -503,12 +583,6 @@ function ReactionNetworkProblem( 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 @@ -526,18 +600,26 @@ function AlgebraicAgents._reinit!(state::ReactionNetworkProblem) 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) + structured_tokens = collect(values(inners(getagent(state, "structured")))) + for (i, species) in enumerate(state.acs[:, :specName]) + if state.acs[i, :specStructured] + state.u[i] = + count(a -> get_species(a) == species && !isblocked(a), structured_tokens) + end end + + return state.u end function AlgebraicAgents._step!(state::ReactionNetworkProblem) free_blocked_species!(state) + update_u_structured!(state) update_observables(state) sample_transitions!(state) evolve!(state) + update_u_structured!(state) finish!(state) + update_u_structured!(state) event_action!(state) diff --git a/src/state.jl b/src/state.jl index ed7a7f9..d0eef41 100644 --- a/src/state.jl +++ b/src/state.jl @@ -18,6 +18,7 @@ Ongoing transition auxiliary structure. bound_structured_agents::Vector{AbstractAlgebraicAgent} nonblock_structured_agents::Vector{AbstractAlgebraicAgent} + structured_to_agents::Vector t::Float64 q::Float64 @@ -46,7 +47,7 @@ end p::Any t::Float64 - structured_species::Vector{Symbol} + structured_token::Vector{Symbol} tspan::Tuple{Float64,Float64} dt::Float64 diff --git a/tutorial/agents-integration/agents.jl b/tutorial/agents-integration/agents.jl index 5300f0d..b021720 100644 --- a/tutorial/agents-integration/agents.jl +++ b/tutorial/agents-integration/agents.jl @@ -35,6 +35,13 @@ network = @ReactionNetworkSchema ρ4, R1 + M1 --> M2(@t(), rand(4)), preAction => update_prob_transition(state, transition) + + # R2 + M1 --> M2(@t(), rand(4)) + 5.0, R2 + M1 --> @structured(M2(@t(), rand(4)), :A) + 1.0, R2 + A --> @structured(M2(@t(), rand(4)), f_species(@transition)) + + 1.0, R2 + M1 --> @move(M1, :M2) + 1.0, R2 + M1 --> @move(M1, :C) end @prob_init network R1 = 10 R2 = 15 @@ -49,7 +56,7 @@ end # We use `@structured` macro, which is a convenience wrapper around `@aagent`), # defined in ReactiveDynamics.jl -@structured network struct M1 +@structured_token network struct M1 descriptor::Any time_created::Any end @@ -60,7 +67,9 @@ using Random # Accordingly, we have to explicitly declare the scope. using ReactiveDynamics: M1 -ReactiveDynamics.M1(time, descriptor) = M1("M1" * randstring(4), nothing, descriptor, time) +function ReactiveDynamics.M1(time, descriptor) + return M1("M1" * randstring(4), :M1, nothing, [], descriptor, time) +end # We define the function which updates the transition probability. # This has to be accessible from within the name scope of ReactiveDynamics. @@ -69,35 +78,41 @@ ReactiveDynamics.M1(time, descriptor) = M1("M1" * randstring(4), nothing, descri if !isnothing(transition) && !isempty(transition.bound_structured_agents) bound_agent = first(transition.bound_structured_agents) - transition[:transProbOfSuccess] = min(1.0, sum(bound_agent.descriptor)) + transition[:transProbOfSuccess] = 1.0#min(1.0, sum(bound_agent.descriptor)) end end + + f_species(transition) = :B 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 + @aagent BaseStructuredToken AbstractStructuredToken struct M2 descriptor::Any time_created::Any end using Random - M2(time, descriptor) = M2("M2" * randstring(4), nothing, descriptor, time) + M2(time, descriptor) = M2("M2" * randstring(4), :M2, nothing, [], descriptor, time) end # Let the network know that the species is structured. -ReactiveDynamics.register_structured_species!(network, :M2) +for species in [:M1, :M2, :A, :B, :C] + ReactiveDynamics.register_structured_species!(network, species) +end # -------------------------------------------------------------------------------- # Instantiate the network. network_instance = ReactionNetworkProblem(network) for i = 1:2 - add_structured_species!(network_instance, ReactiveDynamics.M1(0.0, rand(4))) + add_structured_token!(network_instance, ReactiveDynamics.M1(0.0, rand(4))) end # -------------------------------------------------------------------------------- # Simulate the network. simulate(network_instance, 10) + +# tokens = collect(values(inners(getagent(network_instance, "structured"))))