diff --git a/src/ReactiveDynamics.jl b/src/ReactiveDynamics.jl index a7192bf..c10b7a0 100644 --- a/src/ReactiveDynamics.jl +++ b/src/ReactiveDynamics.jl @@ -84,7 +84,7 @@ const ReactionNetworkSchema = FoldedReactionNetworkType{ Set{Symbol}, FoldedObservable, Any, - Bool + Bool, } Base.convert(::Type{Symbol}, ex::String) = Symbol(ex) diff --git a/src/compilers.jl b/src/compilers.jl index 0c741d8..b4faba8 100644 --- a/src/compilers.jl +++ b/src/compilers.jl @@ -118,7 +118,7 @@ function wrap_expr(fex, species_names, prm_names, varmap) return eval(quote function (state, transition) - $letex + return $letex end end) end diff --git a/src/interface/agents.jl b/src/interface/agents.jl index 6233f63..842a7f8 100644 --- a/src/interface/agents.jl +++ b/src/interface/agents.jl @@ -9,7 +9,7 @@ abstract type AbstractStructuredSpecies <: AbstractAlgebraicAgent end # 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} + bound_transition::Union{Nothing,ReactiveDynamics.Transition} end # We use this to let the network know that the type is structured. @@ -27,16 +27,21 @@ end # Convenience macro to define structured species. macro structured(network, type) name = Docs.namify(type.args[2]) - - quote - $(AlgebraicAgents.aagent(BaseStructuredSpecies, AbstractStructuredSpecies, type, ReactiveDynamics)) + + 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) + return entangle!(getagent(problem, "structured/$(nameof(typeof(agent)))"), agent) end import AlgebraicAgents @@ -49,4 +54,4 @@ AlgebraicAgents._step!(::AbstractStructuredSpecies) = nothing 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 +priority(a, transition) = 0.0 diff --git a/src/interface/reaction_parser.jl b/src/interface/reaction_parser.jl index a2e9ea2..8ad0651 100644 --- a/src/interface/reaction_parser.jl +++ b/src/interface/reaction_parser.jl @@ -3,7 +3,7 @@ using MacroTools: postwalk struct FoldedReactant - species::Union{Expr, Symbol} + species::Union{Expr,Symbol} stoich::SampleableValues modality::Set{Symbol} end diff --git a/src/solvers.jl b/src/solvers.jl index 889ddd3..cb86797 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -36,7 +36,9 @@ function get_reqs_ongoing!(reqs, qs, state) (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]).") + 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 @@ -125,7 +127,6 @@ function get_init_satisfied(allocs, qs, state) return qs end - """ Evolve transitions, spawn new transitions. """ @@ -189,12 +190,20 @@ function evolve!(state) 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.") + error( + "For structured species, stoichiometry coefficient must be integer in transition $i.", + ) end - all_of_type = collect(values(inners(getagent(state, "structured/$(string(type))")))) + 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) + 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) @@ -246,12 +255,20 @@ function evolve!(state) 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.") + error( + "For structured species, stoichiometry coefficient must be integer in transition $i.", + ) end - all_of_type = collect(values(inners(getagent(state, "structured/$(string(type))")))) + 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) + 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) @@ -301,7 +318,7 @@ function finish!(state) ((state.t - trans_.t) < trans_.trans[:transMaxLifeTime]) && (trans_.state < trans_[:transCycleTime]) && (ix += 1; continue) - + q = if trans_.state >= trans_[:transCycleTime] rand(Distributions.Binomial(Int(trans_.q), trans_[:transProbOfSuccess])) else @@ -320,7 +337,7 @@ function finish!(state) state.u[i] += q * stoich val_reward += state[i, :specReward] * q * stoich - for _ in 1:q + for _ = 1:q a = context_eval(state, trans_, state.wrap_fun(r.species)) entangle!(getagent(state, "structured/$(r.species.args[1])"), a) end @@ -340,7 +357,7 @@ function finish!(state) tok.stoich * (in(:rate, tok.modality) ? trans_[:transCycleTime] : 1) if tok.species ∈ state.structured_species - for _ in 1:(trans_.q * tok.stoich) + for _ = 1:(trans_.q*tok.stoich) trans_.bound_structured_agents[begin].bound_transition = nothing deleteat!(trans_.bound_structured_agents, 1) end @@ -349,12 +366,14 @@ function finish!(state) if in(:nonblock, tok.modality) if in(:conserved, tok.modality) - error("Modalities `:conserved` and `:nonblock` cannot be specified at the same time.") + error( + "Modalities `:conserved` and `:nonblock` cannot be specified at the same time.", + ) end - state.u[tok.index] += trans_.q * tok.stoich + state.u[tok.index] += trans_.q * tok.stoich if tok.species ∈ state.structured_species - for _ in 1:(trans_.q * tok.stoich) + for _ = 1:(trans_.q*tok.stoich) trans_.nonblock_structured_agents[begin].bound_transition = nothing deleteat!(trans_.nonblock_structured_agents, 1) end @@ -366,7 +385,7 @@ function finish!(state) 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 @@ -427,7 +446,8 @@ function ReactionNetworkProblem( acs = remove_choose(acs) - structured_species_names = acs[filter(i -> acs[i, :specStructured], 1:nparts(acs, :S)), :specName] + 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 diff --git a/src/state.jl b/src/state.jl index 70ed436..ed7a7f9 100644 --- a/src/state.jl +++ b/src/state.jl @@ -65,7 +65,7 @@ end # evaluate compiled numeric expression in context of (u, p, t) 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 @@ -74,7 +74,8 @@ function Base.getindex(state::ReactionNetworkProblem, keys...) return state.acs[keys[1], keys[2]] else return context_eval( - state, nothing, + state, + nothing, (contains(string(keys[2]), "trans") ? state.transitions : state.attrs)[keys[2]][keys[1]], ) end @@ -176,7 +177,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, :transPreAction, :transPostAction, :transActivated, :transHash]) && continue + ( + attr ∈ + [:trans, :transPreAction, :transPostAction, :transActivated, :transHash] + ) && continue push!( state.transitions[attr], context_eval(state, nothing, state.transition_recipes[attr][i]), @@ -196,10 +200,10 @@ function sample_transitions!(state::ReactionNetworkProblem) ), ) end - + push!(state.transitions[:transLHS], reactants) push!(state.transitions[:transRHS], r_line) - + foreach( k -> push!(state.transitions[k], state.transition_recipes[k][i]), [:transPreAction, :transPostAction, :transToSpawn, :transHash], diff --git a/tutorial/agents-integration/agents.jl b/tutorial/agents-integration/agents.jl index c4304b9..5300f0d 100644 --- a/tutorial/agents-integration/agents.jl +++ b/tutorial/agents-integration/agents.jl @@ -22,17 +22,19 @@ network = @ReactionNetworkSchema # 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) + ρ4, + R1 + M1 --> M2(@t(), rand(4)), + preAction => update_prob_transition(state, transition) end @prob_init network R1 = 10 R2 = 15 @@ -48,8 +50,8 @@ end # We use `@structured` macro, which is a convenience wrapper around `@aagent`), # defined in ReactiveDynamics.jl @structured network struct M1 - descriptor - time_created + descriptor::Any + time_created::Any end using Random @@ -66,7 +68,7 @@ ReactiveDynamics.M1(time, descriptor) = M1("M1" * randstring(4), nothing, descri 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 @@ -77,8 +79,8 @@ end # of ReactiveDynamics. @register begin @aagent BaseStructuredSpecies AbstractStructuredSpecies struct M2 - descriptor - time_created + descriptor::Any + time_created::Any end using Random @@ -92,10 +94,10 @@ ReactiveDynamics.register_structured_species!(network, :M2) # Instantiate the network. network_instance = ReactionNetworkProblem(network) -for i in 1:2 +for i = 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 +simulate(network_instance, 10) diff --git a/tutorial/agents-integration/agents_integration.jl b/tutorial/agents-integration/agents_integration.jl index 25fb6fd..3279902 100644 --- a/tutorial/agents-integration/agents_integration.jl +++ b/tutorial/agents-integration/agents_integration.jl @@ -19,8 +19,8 @@ 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) +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] @@ -31,8 +31,8 @@ for probs in [M1_transition_probs, M2_transition_probs] 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_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] @@ -43,8 +43,8 @@ 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] +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)) @@ -76,10 +76,21 @@ 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") + 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" => 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) @@ -88,7 +99,13 @@ for (i, state) in enumerate(M1_states[begin:end-2]) 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)])) + push!( + to, + Dict( + "probability" => M1_transition_probs[j+i-1, i], + "reactants" => [Dict("name" => state2, "q" => 1)], + ), + ) end end @@ -98,10 +115,21 @@ 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") + 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" => 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) @@ -110,7 +138,13 @@ for (i, state) in enumerate(M2_states[begin:end-2]) 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)])) + push!( + to, + Dict( + "probability" => M2_transition_probs[j+i-1, i], + "reactants" => [Dict("name" => state2, "q" => 1)], + ), + ) end end @@ -137,7 +171,11 @@ 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] + 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], " + ") @@ -200,7 +238,13 @@ using AlgebraicAgents initial_time::Float64 end -function Controller(name::String, M1_λ::Float64, M2_λ::Float64, time::T, time_step::T) where {T<:Real} +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 @@ -214,26 +258,32 @@ function AlgebraicAgents._step!(c::Controller) 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)) + 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)) + 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)), ", ")) + 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.); +c = Controller("controller", 1e-1, 2e-1, 0.0, 1.0); compound_problem = ⊕(problem, c; name = "compound problem"); @@ -253,24 +303,35 @@ draw(compound_problem.inners["network"]) eval(expr_network_def) n_primary_resources = 5 -primary_resources = ["R$i" for i in 1:n_primary_resources] +primary_resources = ["R$i" for i = 1:n_primary_resources] str_resource_making_transitions = [ - "p_primary_$i, ∅ --> $res, name => primary_resource_maker_$i" for (i, res) in enumerate(primary_resources) -] + "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") + 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 = +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)], " ") + join(["p_$i = $(rand(1:3))" for i = 1:length(experimental_resources)], " ") * + " " * + join(["p_primary_$i = $(rand(1:5))" for i = 1:length(primary_resources)], " ") str_network_def = """ begin @@ -286,10 +347,13 @@ 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") +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], " ") +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) @@ -299,4 +363,4 @@ eval(equalize_expr) # Convert network into an AlgAgents hierarchy. extended_problem = ReactionNetworkProblem(extended_network; name = "extended_network") -simulate(extended_problem, 100); \ No newline at end of file +simulate(extended_problem, 100);