diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index 1cb1c01..eb6ff10 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -12,6 +12,7 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v3 + - uses: ts-graphviz/setup-graphviz@v2 - uses: julia-actions/setup-julia@latest with: version: '1.10' diff --git a/Project.toml b/Project.toml index 2ec8279..91ef736 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "AlgebraicAgents" uuid = "f6eb0ae3-10fa-40e6-88dd-9006ba45093a" -version = "0.3.23" +version = "0.3.24" [deps] Crayons = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" @@ -12,8 +12,8 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [compat] +julia = "1.9" Crayons = "4.1" Glob = "1.3" MacroTools = "0.5" Requires = "1.3" -julia = "1.8" diff --git a/docs/Project.toml b/docs/Project.toml index e863d3a..cce1231 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,6 +9,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterMarkdown = "997ab1e6-3595-5248-9280-8efb232c3433" DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1" +Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" diff --git a/docs/make.jl b/docs/make.jl index 8de81b8..4f2e6a2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -33,7 +33,7 @@ end # Literate for tutorials const literate_dir = joinpath(@__DIR__, "..", "tutorials") const generated_dir = joinpath(@__DIR__, "src", "sketches") -const skip_dirs = ["traces"] +const skip_dirs = ["traces", "wires"] for (root, dirs, files) in walkdir(literate_dir) if any(occursin.(skip_dirs, root)) || startswith(root, "_") diff --git a/docs/sketches/agents/agents.md b/docs/sketches/agents/agents.md new file mode 100644 index 0000000..4b96324 --- /dev/null +++ b/docs/sketches/agents/agents.md @@ -0,0 +1,223 @@ +```@meta +EditURL = "../../../../tutorials/agents/agents.jl" +``` + +# Agents.jl Integration + +We instantiate an agent-based SIR model based on [Agents.jl: SIR model for the spread of COVID-19](https://juliadynamics.github.io/Agents.jl/stable/examples/sir/) make use of a SIR model constructor from an Agents.jl' [SIR model for the spread of COVID-19](https://juliadynamics.github.io/Agents.jl/stable/examples/sir/), and then we simulate the model using AlgebraicAgents.jl + +## SIR Model in Agents.jl + +To begin with, we define the Agents.jl model: + +````@example agents +# SIR model for the spread of COVID-19 +# taken from https://juliadynamics.github.io/Agents.jl/stable/examples/sir/ +using AlgebraicAgents +using Agents, Random +using Agents.DataFrames, Agents.Graphs +using Distributions: Poisson, DiscreteNonParametric +using DrWatson: @dict +using Plots + +@agent PoorSoul GraphAgent begin + days_infected::Int ## number of days since is infected + status::Symbol ## 1: S, 2: I, 3:R +end +```` + +Let's provide the constructors: + +````@example agents +function model_initiation(; + Ns, + migration_rates, + β_und, + β_det, + infection_period = 30, + reinfection_probability = 0.05, + detection_time = 14, + death_rate = 0.02, + Is = [zeros(Int, length(Ns) - 1)..., 1], + seed = 0) + rng = MersenneTwister(seed) + @assert length(Ns)== + length(Is)== + length(β_und)== + length(β_det)== + size(migration_rates, 1) "length of Ns, Is, and B, and number of rows/columns in migration_rates should be the same " + @assert size(migration_rates, 1)==size(migration_rates, 2) "migration_rates rates should be a square matrix" + + C = length(Ns) + # normalize migration_rates + migration_rates_sum = sum(migration_rates, dims = 2) + for c in 1:C + migration_rates[c, :] ./= migration_rates_sum[c] + end + + properties = @dict(Ns, + Is, + β_und, + β_det, + β_det, + migration_rates, + infection_period, + infection_period, + reinfection_probability, + detection_time, + C, + death_rate) + space = GraphSpace(complete_digraph(C)) + model = ABM(PoorSoul, space; properties, rng) + + # Add initial individuals + for city in 1:C, n in 1:Ns[city] + ind = add_agent!(city, model, 0, :S) ## Susceptible + end + # add infected individuals + for city in 1:C + inds = ids_in_position(city, model) + for n in 1:Is[city] + agent = model[inds[n]] + agent.status = :I ## Infected + agent.days_infected = 1 + end + end + return model +end + +using LinearAlgebra: diagind + +function create_params(; + C, + max_travel_rate, + infection_period = 30, + reinfection_probability = 0.05, + detection_time = 14, + death_rate = 0.02, + Is = [zeros(Int, C - 1)..., 1], + seed = 19) + Random.seed!(seed) + Ns = rand(50:5000, C) + β_und = rand(0.3:0.02:0.6, C) + β_det = β_und ./ 10 + + Random.seed!(seed) + migration_rates = zeros(C, C) + for c in 1:C + for c2 in 1:C + migration_rates[c, c2] = (Ns[c] + Ns[c2]) / Ns[c] + end + end + maxM = maximum(migration_rates) + migration_rates = (migration_rates .* max_travel_rate) ./ maxM + migration_rates[diagind(migration_rates)] .= 1.0 + + params = @dict(Ns, + β_und, + β_det, + migration_rates, + infection_period, + reinfection_probability, + detection_time, + death_rate, + Is) + + return params +end +```` + +It remains to provide the SIR stepping functions: + +````@example agents +function agent_step!(agent, model) + @get_model model + extract_agent(model, agent) + migrate!(agent, model) + transmit!(agent, model) + update!(agent, model) + recover_or_die!(agent, model) +end + +function migrate!(agent, model) + pid = agent.pos + d = DiscreteNonParametric(1:(model.C), model.migration_rates[pid, :]) + m = rand(model.rng, d) + if m ≠ pid + move_agent!(agent, m, model) + end +end + +function transmit!(agent, model) + agent.status == :S && return + rate = if agent.days_infected < model.detection_time + model.β_und[agent.pos] + else + model.β_det[agent.pos] + end + + d = Poisson(rate) + n = rand(model.rng, d) + n == 0 && return + + for contactID in ids_in_position(agent, model) + contact = model[contactID] + if contact.status == :S || + (contact.status == :R && rand(model.rng) ≤ model.reinfection_probability) + contact.status = :I + n -= 1 + n == 0 && return + end + end +end + +update!(agent, model) = agent.status == :I && (agent.days_infected += 1) + +function recover_or_die!(agent, model) + if agent.days_infected ≥ model.infection_period + if rand(model.rng) ≤ model.death_rate + @a kill_agent!(agent, model) + else + agent.status = :R + agent.days_infected = 0 + end + end +end +```` + +That's it! + +## Simulation Using AlgebraicAgents.jl + +We instantiate a sample `ABM` model: + +````@example agents +# create a sample agent based model +params = create_params(C = 8, max_travel_rate = 0.01) +abm = model_initiation(; params...) +```` + +Let's specify what data to collect: + +````@example agents +infected(x) = count(i == :I for i in x) +recovered(x) = count(i == :R for i in x) +to_collect = [(:status, f) for f in (infected, recovered, length)] +```` + +We wrap the model as an agent: + +````@example agents +m = ABMAgent("sir_model", abm; agent_step!, tspan=(0., 100.), adata=to_collect) +```` + +And we simulate the dynamics: + +````@example agents +simulate(m) +```` + +````@example agents +draw(m) +```` + diff --git a/docs/sketches/algebraicdynamics/algebraicdynamics.md b/docs/sketches/algebraicdynamics/algebraicdynamics.md new file mode 100644 index 0000000..28feeec --- /dev/null +++ b/docs/sketches/algebraicdynamics/algebraicdynamics.md @@ -0,0 +1,180 @@ +```@meta +EditURL = "../../../../tutorials/algebraicdynamics/algebraicdynamics.jl" +``` + +# Lotka-Voltera Two Ways + +We demonstrate an integration of [AlgebraicDynamics.jl](https://github.com/AlgebraicJulia/AlgebraicDynamics.jl). + +The tutorial is based on [AlgebraicDynamics.jl: Lotka-Volterra Three Ways](https://algebraicjulia.github.io/AlgebraicDynamics.jl/dev/examples/Lotka-Volterra/). + +## Undirected Composition + +````@example algebraicdynamics +using AlgebraicDynamics +using AlgebraicDynamics.UWDDynam +using Catlab.WiringDiagrams, Catlab.Programs +using LabelledArrays +using Plots + +const UWD = UndirectedWiringDiagram +```` + +````@example algebraicdynamics +using AlgebraicAgents +```` + +Define the primitive systems + +````@example algebraicdynamics +dotr(u,p,t) = p.α*u +dotrf(u,p,t) = [-p.β*u[1]*u[2], p.γ*u[1]*u[2]] +dotf(u,p,t) = -p.δ*u + +rabbit_growth = wrap_system("rabbit_growth", ContinuousResourceSharer{Float64}(1, dotr)) +rabbitfox_predation = wrap_system("rabbitfox_predation", ContinuousResourceSharer{Float64}(2, dotrf)) +fox_decline = wrap_system("fox_decline", ContinuousResourceSharer{Float64}(1, dotf)) +```` + +Define the composition pattern + +````@example algebraicdynamics +rf = @relation (rabbits,foxes) begin + growth(rabbits) + predation(rabbits,foxes) + decline(foxes) +end +```` + +Compose + +````@example algebraicdynamics +rabbitfox_system = ⊕(rabbit_growth, rabbitfox_predation, fox_decline, diagram=rf, name="rabbitfox_system") +```` + +Solve and plot + +````@example algebraicdynamics +u0 = [10.0, 100.0] +params = LVector(α=.3, β=0.015, γ=0.015, δ=0.7) +tspan = (0.0, 100.0) +```` + +````@example algebraicdynamics +import DifferentialEquations +prob = DiffEqAgent(rabbitfox_system, u0, tspan, params) +```` + +````@example algebraicdynamics +sol = simulate(prob) +```` + +````@example algebraicdynamics +draw(sol; label=["rabbits" "foxes"]) +```` + +## Directed Composition + +````@example algebraicdynamics +using AlgebraicDynamics, AlgebraicDynamics.DWDDynam +using Catlab.WiringDiagrams, Catlab.Programs +using LabelledArrays +using Plots +```` + +````@example algebraicdynamics +using AlgebraicAgents, DifferentialEquations +```` + +Define the primitive systems + +````@example algebraicdynamics +dotr(u, x, p, t) = [p.α*u[1] - p.β*u[1]*x[1]] +dotf(u, x, p, t) = [p.γ*u[1]*x[1] - p.δ*u[1]] + +rabbit = wrap_system("rabbit", ContinuousMachine{Float64}(1,1,1, dotr, (u, p, t) -> u)) +fox = wrap_system("fox", ContinuousMachine{Float64}(1,1,1, dotf, (u, p, t) -> u)) +```` + +Define the composition pattern + +````@example algebraicdynamics +rabbitfox_pattern = WiringDiagram([], [:rabbits, :foxes]) +rabbit_box = add_box!(rabbitfox_pattern, Box(:rabbit, [:pop], [:pop])) +fox_box = add_box!(rabbitfox_pattern, Box(:fox, [:pop], [:pop])) + +add_wires!(rabbitfox_pattern, Pair[ + (rabbit_box, 1) => (fox_box, 1), + (fox_box, 1) => (rabbit_box, 1), + (rabbit_box, 1) => (output_id(rabbitfox_pattern), 1), + (fox_box, 1) => (output_id(rabbitfox_pattern), 2) +]) +```` + +Compose + +````@example algebraicdynamics +rabbitfox_system = ⊕(rabbit, fox; diagram=rabbitfox_pattern, name="rabbitfox_system") +```` + +Solve and plot + +````@example algebraicdynamics +u0 = [10.0, 100.0] +params = LVector(α=.3, β=0.015, γ=0.015, δ=0.7) +tspan = (0.0, 100.0) +```` + +````@example algebraicdynamics +# convert the system to a problem +prob = DiffEqAgent(rabbitfox_system, u0, tspan, params) +```` + +````@example algebraicdynamics +# solve the problem +simulate(prob) +```` + +````@example algebraicdynamics +# plot +draw(prob; label=["rabbits" "foxes"]) +```` + +## Open CPG + +````@example algebraicdynamics +using AlgebraicDynamics.CPortGraphDynam +using AlgebraicDynamics.CPortGraphDynam: barbell +```` + +Define the composition pattern + +````@example algebraicdynamics +rabbitfox_pattern = barbell(1) + +rabbitfox_system = ⊕(rabbit, fox; diagram=rabbitfox_pattern, name="rabbitfox_system") +```` + +Solve and plot + +````@example algebraicdynamics +u0 = [10.0, 100.0] +params = LVector(α=.3, β=0.015, γ=0.015, δ=0.7) +tspan = (0.0, 100.0) +```` + +````@example algebraicdynamics +# convert the system to a problem +prob = DiffEqAgent(rabbitfox_system, u0, tspan, params) +```` + +````@example algebraicdynamics +# solve the problem +simulate(prob) +```` + +````@example algebraicdynamics +# plot +draw(prob; label=["rabbits" "foxes"]) +```` + diff --git a/docs/sketches/molecules/molecules.md b/docs/sketches/molecules/molecules.md new file mode 100644 index 0000000..ab140ca --- /dev/null +++ b/docs/sketches/molecules/molecules.md @@ -0,0 +1,335 @@ +```@meta +EditURL = "../../../../tutorials/molecules/molecules.jl" +``` + +# A Toy Pharma Model + +We implement a toy pharma model. To that end, we have the following type hierarchy: + + - overarching **pharma model** (represented by a `FreeAgent` span type), + - **therapeutic area** (represented by a `FreeAgent`), + - **molecules** (small, large - to demonstrate dynamic dispatch; alternatively, marketed drugs; a drug may drop out from the system), + - **discovery unit** (per therapeutic area); these generate new molecules according to a Poisson counting process, + - **market demand**; this will be represented by a stochastic differential equation implemented in [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl). + +## Agent Types + +We define the type system and supply the stepping functions. + +````@example molecules +using AlgebraicAgents + +using DataFrames +using Distributions, Random +using Plots + +# type hierarchy +"Drug entity, lives in a therapeutic area." +abstract type Molecule <: AbstractAlgebraicAgent end + +# molecule params granularity +const N = 3; + +# drug entity, lives in a therapeutic area +"Parametrized drug entity, lives in a therapeutic area." +@aagent FreeAgent Molecule struct SmallMolecule + age::Float64 + birth_time::Float64 + time_removed::Float64 + + mol::AbstractString + profile::NTuple{N, Float64} + + sales::Float64 + df_sales::DataFrame +end +```` + +Note the use of a conveniency macro `@aagent` which appends additional fields expected (not required, though) by default interface methods. + +We proceed with other agent types. + +````@example molecules +# drug entity, lives in a therapeutic area +@doc (@doc SmallMolecule) +@aagent FreeAgent Molecule struct LargeMolecule + age::Float64 + birth_time::Float64 + time_removed::Float64 + + mol::AbstractString + profile::NTuple{N, Float64} + + sales::Float64 + df_sales::DataFrame +end + +# toy discovery unit - outputs small/large molecules to a given therapeutic area +"Toy discovery unit; outputs small and large molecules." +@aagent struct Discovery + rate_small::Float64 + rate_large::Float64 + discovery_intensity::Float64 + + t::Float64 + dt::Float64 + + t0::Float64 + + removed_mols::Vector{Tuple{String, Float64}} + + df_output::DataFrame +end + +# `Discovery` constructor +"Initialize a discovery unit, parametrized by small/large molecules production rate." +function Discovery(name, rate_small, rate_large, t = 0.0; dt = 2.0) + df_output = DataFrame(time = Float64[], small = Int[], large = Int[], removed = Int[]) + + Discovery(name, rate_small, rate_large, 0.0, t, dt, t, Tuple{String, Float64}[], + df_output) +end +```` + +## Stepping Functions + +Next we provide an evolutionary law for the agent types. This is done by extending the interface function [`AlgebraicAgents._step!`](@ref). + +````@example molecules +# Return initial sales volume of a molecule. +function sales0_from_params end + +const sales0_factor_small = rand(N) +const sales0_factor_large = rand(N) + +# dispatch on molecule type +sales0_from_params(profile) = 10^3 * (1 + collect(profile)' * sales0_factor_small) +sales0_from_params(profile) = 10^5 * (1 + collect(profile)' * sales0_factor_large) + +const sales_decay_small = 0.9 +const sales_decay_large = 0.7 + +# implement evolution +function AlgebraicAgents._step!(mol::SmallMolecule) + t = projected_to(mol) + push!(mol.df_sales, (t, mol.sales)) + mol.age += 1 + mol.sales *= sales_decay_small + + if (mol.sales <= 10) || (rand() >= exp(-0.2 * mol.age)) + mol.time_removed = t + push!(getagent(mol, "../dx").removed_mols, (mol.mol, t)) + + # move to removed candidates + rm_mols = getagent(mol, "../removed-molecules") + disentangle!(mol) + entangle!(rm_mols, mol) + end +end + +# implement common interface" +# molecules" +function AlgebraicAgents._step!(mol::LargeMolecule) + t = projected_to(mol) + push!(mol.df_sales, (t, mol.sales)) + + mol.age += 1 + mol.sales *= sales_decay_large + + if (mol.sales <= 10) || (rand() >= exp(-0.3 * mol.age)) + mol.time_removed = t + push!(getagent(mol, "../dx").removed_mols, (mol.mol, t)) + + # move to removed candidates + rm_mols = getagent(mol, "../removed-molecules") + disentangle!(mol) + entangle!(rm_mols, mol) + end +end + +# discovery +function AlgebraicAgents._step!(dx::Discovery) + t = projected_to(dx) + # sync with market demand + dx.discovery_intensity = getobservable(getagent(dx, "../demand"), "demand") + ν = dx.discovery_intensity * dx.dt + small, large = rand(Poisson(ν * dx.rate_small)), rand(Poisson(ν * dx.rate_large)) + removeed = 0 + ix = 1 + while ix <= length(dx.removed_mols) + if (dx.removed_mols[ix][2] < t) + removeed += 1 + deleteat!(dx.removed_mols, ix) + else + ix += 1 + end + end + push!(dx.df_output, (t, small, large, removeed)) + + for _ in 1:small + mol = release_molecule(randstring(5), Tuple(rand(N)), t, SmallMolecule) + entangle!(getparent(dx), mol) + end + + for _ in 1:large + mol = release_molecule(randstring(5), Tuple(rand(N)), t, LargeMolecule) + entangle!(getparent(dx), mol) + end + + dx.t += dx.dt +end + +"Initialize a new molecule." +function release_molecule(mol, profile, t, ::Type{T}) where {T <: Molecule} + T(mol, 0.0, t, Inf, mol, profile, sales0_from_params(profile), + DataFrame(time = Float64[], sales = Float64[])) +end +```` + +We provide additional interface methods, such as [`AlgebraicAgents._reinit!`](@ref) and [`AlgebraicAgents._projected_to`](@ref). + +````@example molecules +AlgebraicAgents._reinit!(mol::Molecule) = disentangle!(mol) + +function AlgebraicAgents._reinit!(dx::Discovery) + dx.t = dx.t0 + dx.discovery_intensity = 0.0 + empty!(dx.df_output) + + dx +end + +function AlgebraicAgents._projected_to(mol::Molecule) + if mol.time_removed < Inf + # candidate removed, do not step further + true + else + mol.age + mol.birth_time + end +end + +AlgebraicAgents._projected_to(dx::Discovery) = dx.t +```` + +We may also provide a custom plotting recipe. As the internal log is modeled as a DataFrame, we want to use [`AlgebraicAgents.@draw_df`](@ref) convenience macro. + +````@example molecules +# implement plots +AlgebraicAgents.@draw_df Discovery df_output +```` + +## Model Instantiation + +Next step is to instantiate a dynamical system. + +````@example molecules +# define therapeutic areas +therapeutic_area1 = FreeAgent("therapeutic_area1") +therapeutic_area2 = FreeAgent("therapeutic_area2") + +# join therapeutic models into a pharma model +pharma_model = ⊕(therapeutic_area1, therapeutic_area2; name="pharma_model") + +# initialize and push discovery units to therapeutic areas +# discovery units evolve at different pace +entangle!(therapeutic_area1, Discovery("dx", 5.2, 10.; dt=3.)) +entangle!(therapeutic_area2, Discovery("dx", 6., 8.; dt=5.)) +# log removed candidates +entangle!(therapeutic_area1, FreeAgent("removed-molecules")) +entangle!(therapeutic_area2, FreeAgent("removed-molecules")) +```` + +### Integration with SciML + +Let's define toy market demand model and represent this as a stochastic differential equation defined in `DifferentialEquations.jl` + +````@example molecules +# add SDE models for drug demand in respective areas +using DifferentialEquations + +dt = 1//2^(4); tspan = (0.0,100.) +f(u,p,t) = p[1]*u; g(u,p,t) = p[2]*u + +prob_1 = SDEProblem(f,g,.9,tspan,[.01, .01]) +prob_2 = SDEProblem(f,g,1.2,tspan,[.01, .01]) +```` + +Internally, a discovery unit will adjust the candidate generating process intensity according to the observed market demand: + +````@example molecules +# add SDE models for drug demand in respective areas +demand_model_1 = DiffEqAgent("demand", prob_1, EM(); observables=Dict("demand" => 1), dt) +demand_model_2 = DiffEqAgent("demand", prob_2, EM(); observables=Dict("demand" => 1), dt) + +# push market demand units to therapeutic areas +entangle!(therapeutic_area1, demand_model_1) +entangle!(therapeutic_area2, demand_model_2) + +# sync with market demand +getobservable(getagent(pharma_model, "therapeutic_area1/demand"), "demand") +```` + +Let's inspect the composite model: + +````@example molecules +# show the model +pharma_model +```` + +````@example molecules +getagent(pharma_model, glob"therapeutic_area?/") +```` + +## Simulating the System + +Let's next evolve the composite model over a hundred time units. The last argument is optional here; see `?simulate` for the details. + +````@example molecules +# let the problem evolve +simulate(pharma_model, 100) + +getagent(pharma_model, "therapeutic_area1/dx") + +getagent(pharma_model, "therapeutic_area1/demand") +```` + +## Plotting + +We draw the statistics of a Discovery unit in Therapeutic Area 1: + +````@example molecules +draw(getagent(pharma_model, "therapeutic_area1/dx")) +```` + +## Queries + +Let's now query the simulated system. + +To find out which molecules were discovered after time `t=10` and removed from the track before time `t=30`, write + +````@example molecules +pharma_model |> @filter(_.birth_time > 10 && _.time_removed < 30) +```` + +Equivalently, we could make use of f(ilter)-strings, see [`@f_str`](@ref), and write + +````@example molecules +from = 10; to = 30 +pharma_model |> @filter(f"_.birth_time > $from && _.time_removed < $to"); +nothing #hide +```` + +Let's investigate the average life time: + +````@example molecules +# get molecules already removed from the system +removed_molecules = pharma_model |> @filter(f"_.time_removed < Inf") +# calculate `time_removed - birth_time` +# we iterate over `removed_molecules`, and apply the (named) transform function on each agent +# a given agent is referenced to as `_` +life_times = removed_molecules |> @transform(area = getname(getagent(_, "../..")), time=(_.time_removed - _.birth_time)) + +using Statistics: mean +avg_life_time = mean(x -> x.time, life_times) +```` + diff --git a/docs/sketches/sciml/sciml.md b/docs/sketches/sciml/sciml.md new file mode 100644 index 0000000..49e8791 --- /dev/null +++ b/docs/sketches/sciml/sciml.md @@ -0,0 +1,93 @@ +```@meta +EditURL = "../../../../tutorials/sciml/sciml.jl" +``` + +# SciML Integration + +````@example sciml +using AlgebraicAgents +```` + +````@example sciml +# declare problems (models in AA's type system) +using DifferentialEquations + +# vanilla function +f(u,p,t) = 1.01*u +u0 = 1/2 +tspan = (0.0,10.0) +prob = ODEProblem(f,u0,tspan) +```` + +## Atomic Models + +````@example sciml +m1 = DiffEqAgent("model1", prob) +m2 = DiffEqAgent("model2", prob) +m3 = DiffEqAgent("model3", prob) + +# declare observables (out ports) for a model +# it will be possible to reference m3's first variable as both `o1`, `o2` +push!(observables(m3), "o1" => 1, "o2" => 1) + +# simple function, calls to which will be scheduled during the model integration +custom_function(agent, t) = println("inside $agent at time $t") + +# a bit more intricate logic - +function f_(u,p,t) + # access the wrapping agent (hierarchy bond) + agent = extract_agent(p) + + # access observables + o1 = getobservable(getagent(agent, "../model3"), "o1") + o2 = getobservable(getagent(agent, "../model3"), "o2") + # fetch observable's value at **a given time point in the past** + o3 = gettimeobservable(getagent(agent, "../model3"), t/2, 1) + + # schedule interaction + # first, schedule a call to `_interact!(agent)` with priority 0 + # this is the default behavior + poke(agent) + # alternatively, provide a function call f(args...) + # this will be expanded to a call f(agent, args...) + @call agent custom_function(agent, t) + + min(2., 1.01*u + o1 + o2 + o3) +end +```` + +## Another Atomic Model + +````@example sciml +m4 = DiffEqAgent("model4", ODEProblem(f_,u0,tspan)) +```` + +## Hierarchical Sum of Atomic Models + +````@example sciml +m = ⊕(m1, m2; name="diagram1") ⊕ ⊕(m3, m4; name="diagram2") +```` + +````@example sciml +# explore path-like structure of agents + +# index by unix-like path +getagent(m, "diagram1/model1/") +getagent(m, "diagram1/model1") +getagent(m1, "../model2") +getagent(m1, "../../diagram2/model3") + +# index by regex expression +getagent(m, r"model.*") + +# index by glob expression +getagent(m, glob"**/model?/") +getagent(m, glob"**/model?"s) +```` + +## Solving + +````@example sciml +sol = simulate(m) +```` + diff --git a/docs/sketches/stochastic_simulation/anderson.md b/docs/sketches/stochastic_simulation/anderson.md new file mode 100644 index 0000000..8a5a658 --- /dev/null +++ b/docs/sketches/stochastic_simulation/anderson.md @@ -0,0 +1,389 @@ +```@meta +EditURL = "../../../../tutorials/stochastic_simulation/anderson.jl" +``` + +# Simulating Stochastic Reaction Systems + +To demonstrate the generality of AlgebraicAgents.jl, we demonstrate here how +to use the package to set up a type system capable of simulating continuous +time discrete state stochastic processes using the method described by +[Anderson (2007)](https://people.math.wisc.edu/~dfanderson/papers/AndNRM.pdf). + +We begin by importing packages we will use. + +````@example anderson +using AlgebraicAgents, Distributions, DataFrames, Plots +```` + +## Reaction System +We use the `@aagent` struct to define a new type which is a concrete subtype of `AbstractAlgebraicAgent` +called `ReactionSystem`. It contains data members: + + - `t`: current simulation time + - `Δ`: the interarrival time between events + - `X`: current system state + - `X0`: initial system state + - `df_output`: a `DataFrame` contining the sampled trajectory + +````@example anderson +@aagent struct ReactionSystem{T,S} + t::T + Δ::T + X::Vector{S} + X0::Vector{S} + df_output::DataFrame +end +```` + +We define a method `make_reactionsystem` which constructs a concrete instantiation +of the `ReactionSystem` type, with a given `name` and initial state `X0`. +Note that we use `entangle!` to add an instantiation of the `FreeAgent` type +exported from AlgebraicAgents to the object. This agent is called "clocks" and +will contain agents which jointly make up the stochastic dynamics of the system. + +````@example anderson +function make_reactionsystem(name::T, X0::Vector{S}) where {T,S} + df_output = DataFrame(time=Float64[],clock=String[]) + for i in eachindex(X0) + insertcols!(df_output, Symbol("X"*string(i))=>S[]) + end + rs = ReactionSystem{Float64,S}(name, 0.0, 0.0, X0, X0, df_output) + entangle!(rs, FreeAgent("clocks")) + return rs +end +```` + +Because the `ReactionSystem` itself has no dynamics (it represents the "world state"), +its implementation of `AlgebraicAgents._step!` does nothing. + +````@example anderson +AlgebraicAgents._step!(a::ReactionSystem) = nothing +```` + +We also need to implement `AlgebraicAgents._projected_to` for `ReactionSystem`. In +this case because the times to which the system is solved are determined by the individual +stochastic processes which make up the system (defined later), +we can set it to the trivial implementation which does nothing. + +````@example anderson +AlgebraicAgents._projected_to(a::ReactionSystem) = nothing +```` + +## Clock Process + +The key concept in [Anderson (2007)](https://people.math.wisc.edu/~dfanderson/papers/AndNRM.pdf) +is that the stochastic system is defined by a set of clocks, each of which fires +at the points of an inhomogeneous Poisson process. Strictly speaking, each clock +process $k$ also has an associated marking $\nu_{k}$, which updates the state of +the system. Let $M$ be the total number of clock processes. + +Consider the state $X(t)$ at a time $t$ to be a vector of integers. Let each clock process +have an associated counting process $R_{k}(t)$ which tells us the number of times +it has fired up to time $t$. Then we can write the current model state as a function +of the set of all counting processes, their markings, and the initial condition as: + +```math +X(t) = X(0) + \sum_{k=1}^{M} R_{k}(t) \nu_{k} +``` + +We can write each counting process as arising from an inhomogenous Poisson process +with intensity function $a_{k}(X)$. Specific forms of $a_{k}$ and $\nu_{k}$ will make meaningful +models for chemical reactions, ecological systems, epidemiological processes, sociology, +or other domains. + +```math +R_{k}(t) = Y_{k}\left(\int_{0}^{t}a_{k}(X(s))ds\right) +``` +The above is therefore an expression of the counting process in terms of a unit rate +Poisson process $Y_{k}$. Via the random time change theorem by transforming time according +to the integrated intensity $T_{k}(t) = \int_{0}^{t}a_{k}(X(s))ds$ we get the proper +inhomogeneous Poisson process, such that when the intensity is high, more events occur +(i.e. the inter-arrival time is smaller), and vice-versa for low intensity. + +To apply the method of [Anderson (2007)](https://people.math.wisc.edu/~dfanderson/papers/AndNRM.pdf) +we need one more definition, $P_{k}$, which is the _next_ firing time of $T_{k}$, applying +the random time change such that it advances at unit exponentially-distributed increments. + +```math + P_{k} = \{s > T_{k} : Y_{k}(s) > Y_{k}(T_{k}) \} +``` + +Now let us define $R_{k}$ using the `@aagent` macro from AlgebraicAgents. +It contains data members: + + - `P`: the next internal firing time of the homogeneous Poisson process + - `T`: the internal time of the homogeneous Poisson process + - `Δt`: absolute (wall-clock) time time before next putative firing + - `τ`: next asbolute putative time to fire + - `a`: current value of $a_{k}$ + - `intensity`: the intensity function $a_{k}$, which accepts as input $X$ and returns floating point value + - `ν`: the marking which updates state (a vector) + +````@example anderson +@aagent struct Clock{N<:Real,Fn<:Function,S<:Number} + P::N + T::N + Δt::N + τ::N + a::N + intensity::Fn + ν::Vector{S} +end +```` + +We must now define a method which adds a clock process to an object of +type `ReactionSystem`. This "initialization" method implements steps 1-5 +of Algorithm 3 in Anderson's paper. Readers should note that we include +step 5 as initialization and then compute it again at the _end_ of the +loop via using a control interaction because of how AlgebraicAgents structures +updates. + +````@example anderson +function add_clock!(rs::ReactionSystem, name::T, intensity::U, ν::Vector{S}) where {T,U,S} + c = Clock{Float64,U,S}(name, 0.0, 0.0, 0.0, 0.0, 0.0, intensity, ν) + + c.a = c.intensity(rs.X) + c.P = rand(Exponential()) + c.Δt = (c.P - c.T) / c.a + c.τ += c.Δt + + entangle!(inners(rs)["clocks"], c) + + add_control!(rs, () -> control_clock(c), "control " * name) +end +```` + +We must implement `AlgebraicAgents._projected_to` for the type `Clock`. +Here it will return the putative times to fire. This is because the simulation method +(and many other continuous time discrete event simulation algorithms) updates +via a race condition, where the next event to cause a state change is the one +that "actually happens". AlgebraicAgents selects the agent(s) to `_step!` as the +ones whose projected time is the minimum of all projected times, so that the +clock that fires first will be the one whose dynamics occur for this next iteration +of the simulation loop. Because the next times are sampled from the points of Poisson processes +they almost surely occur at unique times, and therefore conflicts cannot occur. + +It is interesting to note here that AlgebraicAgents can implement any algorithm +that depends on such a race condition. + +````@example anderson +function AlgebraicAgents._projected_to(c::Clock) + c.τ +end +```` + +Now we implement `AlgebraicAgents._step!` for type `Clock`. In this method, steps 6,7, and 9 +of Algorithm 3 from the paper are implemented (specific order of 8 and 9 is not important). +Basically, we update the global time to the time this clock fired, update the state $X$, +and draw the new next firing time $P_{k}$. We also push output to the top level agent, which +is of type `ReactionSystem`. Each update log will have the time of the event, the name of +the clock that caused it, and the new system state. + +````@example anderson +function AlgebraicAgents._step!(c::Clock) + + if isinf(c.τ) + return nothing + end + + topmost(c).Δ = c.τ - topmost(c).t + topmost(c).t = c.τ + topmost(c).X += c.ν + c.P += rand(Exponential()) + + push!(topmost(c).df_output, [topmost(c).t, getname(c), topmost(c).X...]) +end +```` + +Finally we must implement the control interaction which is applied to each clock at the end +of an iteration in the loop. This implements steps 8,9, and 5 of Algorithm 3 (note that we +are allowed to move step 5 to the end because we also included it in the "initialization" phase +earlier). It also updates the putative next firing time. + +````@example anderson +function control_clock(c::Clock) + c.T += c.a * topmost(c).Δ + c.a = c.intensity(topmost(c).X) + c.Δt = (c.P - c.T) / c.a + c.τ = topmost(c).t + c.Δt +end +```` + +## Simulation + +We will simulate a continuous time stochastic SIR (Susceptible-Infectious-Recovered) +model. For mathematical details on this model, please consult [Allen (2017)](https://www.sciencedirect.com/science/article/pii/S2468042716300495). +Let us use parameters $\beta$ to represent the effective contact rate, and +$\gamma$ to represent the recovery rate. + +````@example anderson +β = 0.05*10.0/1000 +γ = 0.25 +```` + +Now we make a `ReactionSystem` object, and initialize it to a population with +990 susceptible persons, 10 infectious persons, and 0 recovered persons. +The two events in the system are infection and recovery, which fire according +to the rates given by the anonymous functions passed to `add_clock!`. + +````@example anderson +rs = make_reactionsystem("SIR", [990, 10, 0]) +add_clock!(rs, "infection", (x) -> β*x[2]*x[1], [-1,1,0]) +add_clock!(rs, "recovery", (x) -> γ*x[2], [0,-1,1]) +```` + +Now we call `simulate` on the constructed system. Because a clock will return a next +event time of `Inf` when its rate is 0 (meaning it will never fire again), when all clocks +return `Inf` it means the simulation is over, because nothing else can happen. Therefore +we pass as the second argument to `simulate` the largest representable floating point +value. When all clocks return `Inf`, the minimum will be larger than this value and +the simulation loop will end. + +````@example anderson +simulate(rs, floatmax(Float64)) +```` + +After simulation is complete, we can extract the simulated trajectory. + +````@example anderson +df_out = select(rs.df_output, Not(:clock)); +plot(df_out[!,:time], Matrix(df_out[:,[:X1,:X2,:X3]]), label = ["S" "I" "R"]) +```` + +## Stochastic Petri Net + +Stochastic Petri nets (SPN) are a mathematical language to describe distributed systems which evolve according +to a stochastic trajectory. There are many ways to define them, and for a comprehensive overview of their modeling +power, we reccomend [Haas (2002)](https://link.springer.com/book/10.1007/b97265). We will implement +a very simple SPN to set up a state transition system. Our SPN is nearly identical to the category +of Petri net proposed by [Kock (2023)](https://arxiv.org/abs/2005.05108), with the addition of a rate +parameter associated with each transition. When we assume that overall transition rates occur according to the +mass action law multiplied by the rate constant associated with that transition, we will be able to +produce a `ReactionSystem` that can be simulated using the code above. + +````@example anderson +@aagent struct StochasticPetriNet + P::Vector{Symbol} + T::Vector{Symbol} + I::Int + O::Int + + ip::Vector{Symbol} + it::Vector{Symbol} + op::Vector{Symbol} + ot::Vector{Symbol} + + rate::Dict{Symbol,Float64} +end +```` + +The `StochasticPetriNet` has objects corresponding to (Sets) of places, transitions, input and output arcs. There +are mappings (Functions) which indicate which place or transition each input (output) arc is connected to. For example +`ip` is of length `I`, such that each input arc identifies which place is is connected to (likewise for `it`, but for +transitions). Instead of arc multiplicites, we duplicate arcs, which has the same effect, and simplifies the code. + +We write a helper function to construct SPNs. It is only responsible for checking our input makes sense. + +````@example anderson +function make_stochasticpetrinet(name, P, T, I, O, ip, it, op, ot, rate) + @assert length(T) == length(rate) + @assert all([p ∈ P for p in ip]) + @assert all([t ∈ T for t in it]) + @assert all([p ∈ P for p in op]) + @assert all([t ∈ T for t in ot]) + @assert I == length(ip) + @assert I == length(it) + @assert O == length(op) + @assert O == length(ot) + StochasticPetriNet(name, P, T, I, O, ip, it, op, ot, rate) +end +```` + +The structural components of the SIR model are all in the SPN generated below. Note that there are two output arcs +from the "infection" transition, to the "I" compartment. This is the same as having a single arc of multiplicity 2, +we model arcs "individually" here only to make the code cleaner and more readable. + +````@example anderson +sir_spn = make_stochasticpetrinet( + "SIR", [:S,:I,:R], [:inf,:rec], + 3, 3, + [:S,:I,:I], [:inf,:inf,:rec], + [:I,:I,:R], [:inf,:inf,:rec], + Dict((:inf => β), (:rec => γ)) +) +```` + +Now we can write a function which generates a `ReactionSystem` from our SPN, assuming the law of mass action. +The argument `X0` is the initial marking. + +Note that in this simple example, we do not check the logical "enabling rules" for each transition, we directly +compute the current rate/intensity. Because the net assumes the law of mass action, the computed rate will +equal zero when the transition is not enabled, but this is not true of more general SPNs. A complete implementation +would compute enabling rules from input arcs, and require the user to specify the rate as a `Function` that computed +the intensity of that transition if the enabling rule for that transition evaluated to `true`. We would also +want to apply the "consumption" of input tokens and the "production" of output tokens seperately, rather +than compute the difference of consumption and production as the overall difference, as done here. + +````@example anderson +function generate_reaction_system(spn::StochasticPetriNet, X0) + + mass_action_rs = make_reactionsystem(getname(spn), X0) + + # for each transition, we must make a stochastic clock in the reaction system + for t in spn.T + # get the vector of preconditions (number of times each place is an input for this transition) + precond = zeros(Int, length(spn.P)) + # get a vector of input indices + precond_ix = Int[] + for i in eachindex(spn.it) + if spn.it[i] != t + continue + else + push!(precond_ix, findfirst(isequal(spn.ip[i]), spn.P)) + precond[precond_ix[end]] += 1 + end + end + # get the vector of postconditions (number of times each places is an output for this transition) + postcond = zeros(Int, length(spn.P)) + for i in eachindex(spn.ot) + if spn.ot[i] != t + continue + else + postcond[findfirst(isequal(spn.op[i]), spn.P)] += 1 + end + end + # total change to the marking as a result of transition t + change = postcond - precond + # add a stochastic clock to the reaction system for transition t + add_clock!( + mass_action_rs, String(t), (x) -> prod(x[precond_ix])*spn.rate[t], change + ) + end + + return mass_action_rs +end +```` + +Now we can generate the reaction system which implements the stochastic dynamics of the SIR model from +the Petri net representing the structural constraints of the SIR model. In this way, we seperate specification +of structure from specification of dynamics. We use the same initial condition as before. + +````@example anderson +x0 = [990, 10, 0] +sir_rs = generate_reaction_system(sir_spn, x0) +```` + +We now run another simulation. + +````@example anderson +simulate(sir_rs, floatmax(Float64)) +```` + +We can make another plot. Although the parameters are the same, the stochastic trajectory should look a little different, +due to the randomness in the two driving Poisson processes. + +````@example anderson +df_out = select(sir_rs.df_output, Not(:clock)); +plot(df_out[!,:time], Matrix(df_out[:,[:X1,:X2,:X3]]), label = ["S" "I" "R"]) +```` + diff --git a/docs/src/design.md b/docs/src/design.md index 1346d98..6f134d7 100644 --- a/docs/src/design.md +++ b/docs/src/design.md @@ -157,4 +157,111 @@ poke(agent, 1.) # call `_interact!(agent)`; this call is added to the instantiou # `@call` bob_agent = only(getagent(agent, r"bob")) @call agent wake_up(bob_agent) # translates into `() -> wake_up(bob_agent)` with priority 0 -``` \ No newline at end of file +``` + +### Wires + +It is possible to explicitly establish oriented "wires" along which information flows between different agents in a hierarchy. Note that it is possible to retrieve and modify the state of any other agent from within any agent, in any way. However, in some cases, it may be desirable to explicitly specify that certain agents observe a particular state variable of another agent. + +Consider the following example. First, we set up the hierarchy. + +```@setup wires +using AlgebraicAgents +@aagent struct MyAgentType end +``` + +```@example wires +alice = MyAgentType("alice") +alice1 = MyAgentType("alice1") +entangle!(alice, alice1) + +bob = MyAgentType("bob") +bob1 = MyAgentType("bob1") +entangle!(bob, bob1) + +joint_system = ⊕(alice, bob, name = "joint system") +``` + +We then add the wires. Note that the agents connected by a wire can be specified using the respective agent objects, relative paths, or their UUIDs. + +Additionally, you can assign names to the edges of a wire (which are `nothing` by default). These names can subsequently be used to fetch the information incoming along an edge, a process that we will describe below. + +```@example wires +add_wire!(joint_system; from=alice, to=bob, from_var_name="alice_x", to_var_name="bob_x") +add_wire!(joint_system; from=bob, to=alice, from_var_name="bob_y", to_var_name="alice_y") + +add_wire!(joint_system; from=alice, to=alice1, from_var_name="alice_x", to_var_name="alice1_x") +add_wire!(joint_system; from=bob, to=bob1, from_var_name="bob_x", to_var_name="bob1_x") +``` + +We list the wires going from and to `alice` and `alice1`, respectively. + +```@example wires +get_wires_from(alice) +``` + +```@example wires +get_wires_to(alice1) +``` + +All the wires within an hierarchy can be retrieved as follows: + +```@example wires +getopera(joint_system).wires +``` +Given an agent, if [`getobservable`](@ref) is implemented for all agents that feed information into the specific agent, we can fetch the values incoming to it. + +```@example wires +AlgebraicAgents.getobservable(a::MyAgentType, args...) = getname(a) + +retrieve_input_vars(alice1) +``` +Additionally, we can plot the agents in the hierarchy, displaying the links between parents and children, along with the wires. The output can be adjusted as needed. Note that it is also possible to export an agent hierarchy as a Mermaid diagram. See [`typetree_mmd`](@ref) and [`agent_hierarchy_mmd`](@ref). + +In what follows, [`wiring_diagram`](@ref) generates a visually appealing Graphviz diagram. + +```@example wires +graph1 = wiring_diagram(joint_system) + +run_graphviz("gv1.svg", graph1) +``` + +![](gv1.svg) + +```@example wires +# Do not show edges between parents and children. +graph2 = wiring_diagram(joint_system; parentship_edges=false) + +run_graphviz("gv2.svg", graph2) +``` + +![](gv2.svg) + + +```@example wires +# Only show listed agents. +graph3 = wiring_diagram([alice, alice1, bob, bob1]) + +run_graphviz("gv3.svg", graph3) +``` + +![](gv3.svg) + +```@example wires +# Group agents into two clusters. +graph4 = wiring_diagram([[alice, alice1], [bob, bob1]]) + +run_graphviz("gv4.svg", graph4) +``` + +![](gv4.svg) + +```@example wires +# Provide labels for clusters. +graph5 = wiring_diagram([[alice, alice1], [bob, bob1]]; group_labels=["alice", "bob"], parentship_edges=false) + +run_graphviz("gv5.svg", graph5) +``` + +![](gv5.svg) + diff --git a/docs/src/design_mmd.md b/docs/src/design_mmd.md new file mode 100644 index 0000000..cd19c64 --- /dev/null +++ b/docs/src/design_mmd.md @@ -0,0 +1,273 @@ +```@raw html + +``` +```@raw html + +``` +```@raw html + +``` +# Framework design + +Here we describe the design principles of the AlgebraicAgents. It should be of most use to advanced users and persons interested in contributing to the software. New users are encouraged to start by reading one of the tutorials ("sketches"). + +## Simulation loop + +We describe here the main simulation loop which steps models built in AlgebraicAgents forward in time. + +AlgebraicAgents keeps agents synchronized by ensuring that the model will only simulate the agent(s) whose projected time (e.g. the maximum time for which that agent's trajectory has been solved) is the minimum of all agents. For example, if there are 3 agents, whose internal step sizes are of 1, 1.5, and 3 time units, respectively, then at the end of the first step, their projected time will be 1, 1.5, and 3 (assuming they all start at time 0). The simulation will find the minimum of those times and only simulate the agent(s) whose projected time(s) are equal to the minimum. In this case, it is the first agent, who is now projected to time 2 on the second step (other agents are not simulated). Now, the new minimum time is 1.5, and the second agent is projected to time 3 on the third step, etc. The simulation continues until all agents have been projected to their final time point, or the minimum of projected times reaches a maximum time horizon. If all agents take time steps of the same size, then they will all be updated each global step. + +There are several functions in the interface for an [`AbstractAlgebraicAgent`](@ref) which implement these dynamics. When defining new agent types, one should implement the [`AlgebraicAgents._step!`](@ref) method, which will step that agent forward if its projected time is equal to the least projected time, among all agents in the hierarchy. Agent types also need to implement [`AlgebraicAgents._projected_to`](@ref), which is crucial to keeping the simulation synchronized. It will return: + + * `nothing` if the agent does not implement its own `_step!` rule (e.g. [`FreeAgent`](@ref) which is a container of agents) + * `true` if the agent has been projected to the final time point (`_step!` will not be called again) + * a value of `Number`, giving the time point to which the agent has been projected + +These are collected into `ret`, which is an object that will be `true` if and only if all agents have returned `true`, and is otherwise the minimum of the numeric values (projection times) returned from each inner agent's step. + +```@raw html +
+ +flowchart TD + + Start((Enter Program))-->Project[Set t equal to minimum \n projected time]:::GreenNode + + Project-->RootDecision1{is root?}:::YellowNode + + RootDecision1 -->|yes| PreWalk[Prestep inner agents]:::GreenNode + + RootDecision1 -->|no| Step[Step inner agents]:::GreenNode + + PreWalk -.->|_prestep!| Inners([<:AbstractAlgebraicAgent]):::RedNode + + PreWalk --> Step + + Step -.->|step!| Inners + + subgraph inners + Inners + end + + Ret([ret]):::RedNode + + Inners -.->|_projected_to| Ret + + Step --> LocalDecision{local projected time == t\n equals the min projected time}:::YellowNode + + LocalDecision -->|yes| LocalStep[Local step]:::GreenNode + LocalDecision -->|no| RootDecision2{is root?}:::YellowNode + + LocalStep -.->|_projected_to| Ret + + LocalStep --> RootDecision2 + + subgraph Opera + + RootDecision2 -->|yes| InstantOpera[Execute instantaneous interactions]:::GreenNode + InstantOpera --> FutureOpera[Execute delayed interactions]:::GreenNode + FutureOpera --> ControlOpera[Execute control interactions]:::GreenNode + end + + Opera -.->|_projected_to| Ret + + RootDecision2 -->|no| Stop + + ControlOpera --> Stop((Exit program and\n return ret)) + + classDef GreenNode fill:#D5E8D4,stroke:#82B366; + classDef RedNode fill:#F8CECC,stroke:#B85450; + classDef YellowNode fill:#FFE6CC,stroke:#D79B00; + +
+``` + +Above we show a caricature of the main simulation loop. "Enter program" corresponds to the call to `simulate`, the value of `ret` is (typically) initialized to `0.0`. The simulation continues to step while `ret` is not `true` (meaning the maximum time horizon has been reached by the slowest agent), or has not exceeded some maximum. + +The inner area enclosed by a dashed border represents where program control is given to the `step!` method. The root agent applies `_prestep!` recurvisely to all of its inner (enclosed) agents. After this, `step!` is then applied to all inner agents, and `ret` is updated by each of them. Then the agent applies its own local update `_step!` if its own projected time is equal to the minimum of all inner agent projected times (not shown). Then the Opera module for additional interactions is called for the root agent. + +## Opera + +The Opera system allows interactions between agents to be scheduled. By default, AlgebraicAgents.jl provides support for three types of interactions: + + * **futures (delayed interactions)** + * **system controls** + * **instantious interactions** + +For more details, see the API documentation of [`Opera`](@ref) and our tests. + +### Futures +You may schedule function calls, to be executed at predetermined points of time. +An action is modeled as a tuple `(id, call, time)`, where `id` is an optional textual identifier of the action and `call` is a (parameterless) anonymous function, which will be called at the given `time`. +Once the action is executed, the return value with corresponding action id and execution time is added to `futures_log` field of `Opera` instance. + +See [`add_future!`](@ref) and [`@future`](@ref). + +#### Example + +```julia +alice = MyAgentType("alice") +interact = agent -> wake_up!(agent) +@future alice 5.0 interact(alice) "alice_schedule" +``` + +The solver will stop at `t=5` and call the function `() -> interact(alice)` (a closure is taken at the time when `@future` is invoked). This interaction is identified as `"alice_schedule"`. + +### Control Interactions +You may schedule control function calls, to be executed at every step of the model. +An action is modeled as a tuple `(id, call)`, where `id` is an optional textual identifier of the action, and `call` is a (parameterless) anonymous function. +Once the action is executed, the return value with corresponding action id and execution time is added to `controls_log` field of `Opera` instance. + +See [`add_control!`](@ref) and [`@control`](@ref). + +#### Example + +```julia +system = MyAgentType("system") +control = agent -> agent.temp > 100 && cool!(agent) +@control system control(system) "temperature control" +``` + +At each step, the solver will call the function `() -> control(system)` (a closure is taken at the time when `@future` is invoked). + +### Instantious Interactions +You may schedule additional interactions which exist within a single step of the model; +such actions are modeled as named tuples `(id, priority=0., call)`. Here, `call` is a (parameterless) anonymous function. + +They exist within a single step of the model and are executed after the calls +to `_prestep!` and `_step!` finish, in the order of the assigned priorities. + +In particular, you may schedule interactions of two kinds: + + - `poke(agent, priority)`, which will translate into a call `() -> _interact!(agent)`, with the specified priority, + - `@call opera expresion priority`, which will translate into a call `() -> expression`, with the specified priority. + +See [`poke`](@ref) and [`@call`](@ref). + +#### Examples + +```julia +# `poke` +poke(agent, 1.) # call `_interact!(agent)`; this call is added to the instantious priority queue with priority 1 +``` + +```julia +# `@call` +bob_agent = only(getagent(agent, r"bob")) +@call agent wake_up(bob_agent) # translates into `() -> wake_up(bob_agent)` with priority 0 +``` + +### Wires + +It is possible to explicitly establish oriented "wires" along which information flows between different agents in a hierarchy. Note that it is possible to retrieve and modify the state of any other agent from within any agent, in any way. However, in some cases, it may be desirable to explicitly specify that certain agents observe a particular state variable of another agent. + +Consider the following example. First, we set up the hierarchy. + +```@setup wires +using AlgebraicAgents +@aagent struct MyAgentType end +``` + +```@example wires +alice = MyAgentType("alice") +alice1 = MyAgentType("alice1") +entangle!(alice, alice1) + +bob = MyAgentType("bob") +bob1 = MyAgentType("bob1") +entangle!(bob, bob1) + +joint_system = ⊕(alice, bob, name = "joint system") +``` + +We then add the wires. Note that the agents connected by a wire can be specified using the respective agent objects, relative paths, or their UUIDs. + +Additionally, you can assign names to the edges of a wire (which are `nothing` by default). These names can subsequently be used to fetch the information incoming along an edge, a process that we will describe below. + +```@example wires +add_wire!(joint_system; from=alice, to=bob, from_var_name="alice_x", to_var_name="bob_x") +add_wire!(joint_system; from=bob, to=alice, from_var_name="bob_y", to_var_name="alice_y") + +add_wire!(joint_system; from=alice, to=alice1, from_var_name="alice_x", to_var_name="alice1_x") +add_wire!(joint_system; from=bob, to=bob1, from_var_name="bob_x", to_var_name="bob1_x") +``` + +We list the wires going from and to `alice` and `alice1`, respectively. + +```@example wires +get_wires_from(alice) +``` + +```@example wires +get_wires_to(alice1) +``` + +All the wires within an hierarchy can be retrieved as follows: + +```@example wires +getopera(joint_system).wires +``` +Given an agent, if [`getobservable`](@ref) is implemented for all agents that feed information into the specific agent, we can fetch the values incoming to it. + +```@example wires +AlgebraicAgents.getobservable(a::MyAgentType, args...) = getname(a) + +retrieve_input_vars(alice1) +``` +Additionally, we can plot the agents in the hierarchy, displaying the links between parents and children, along with the wires. The output can be adjusted as needed. Note that it is also possible to export an agent hierarchy as a Mermaid diagram. See [`typetree_mmd`](@ref) and [`agent_hierarchy_mmd`](@ref). + +In what follows, [`wiring_diagram`](@ref) generates a visually appealing Graphviz diagram. + +```@example wires +graph1 = wiring_diagram(joint_system) + +run_graphviz("gv1.svg", graph1) +``` + +![](gv1.svg) + +```@example wires +# Do not show edges between parents and children. +graph2 = wiring_diagram(joint_system; parentship_edges=false) + +run_graphviz("gv2.svg", graph2) +``` + +![](gv2.svg) + + +```@example wires +# Only show listed agents. +graph3 = wiring_diagram([alice, alice1, bob, bob1]) + +run_graphviz("gv3.svg", graph3) +``` + +![](gv3.svg) + +```@example wires +# Group agents into two clusters. +graph4 = wiring_diagram([[alice, alice1], [bob, bob1]]) + +run_graphviz("gv4.svg", graph4) +``` + +![](gv4.svg) + +```@example wires +# Provide labels for clusters. +graph5 = wiring_diagram([[alice, alice1], [bob, bob1]]; group_labels=["alice", "bob"], parentship_edges=false) + +run_graphviz("gv5.svg", graph5) +``` + +![](gv5.svg) + diff --git a/docs/src/index.md b/docs/src/index.md index 4826060..cf2c60d 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -177,6 +177,17 @@ flatten @draw_df ``` +### Wires and wiring diagrams + +```@docs +get_wires_from +get_wires_to +add_wire! +delete_wires! +retrieve_input_vars +wiring_diagram +``` + ### Helper functions for Mermaid diagrams ```@docs diff --git a/docs/src/sketches/agents/agents.md b/docs/src/sketches/agents/agents.md index 4b96324..ea54737 100644 --- a/docs/src/sketches/agents/agents.md +++ b/docs/src/sketches/agents/agents.md @@ -20,7 +20,7 @@ using Distributions: Poisson, DiscreteNonParametric using DrWatson: @dict using Plots -@agent PoorSoul GraphAgent begin +@agent struct PoorSoul(GraphAgent) days_infected::Int ## number of days since is infected status::Symbol ## 1: S, 2: I, 3:R end @@ -68,7 +68,7 @@ function model_initiation(; C, death_rate) space = GraphSpace(complete_digraph(C)) - model = ABM(PoorSoul, space; properties, rng) + model = ABM(PoorSoul, space; properties, rng, model_step! = identity) # Add initial individuals for city in 1:C, n in 1:Ns[city] diff --git a/src/AlgebraicAgents.jl b/src/AlgebraicAgents.jl index 4a0604e..7de1681 100644 --- a/src/AlgebraicAgents.jl +++ b/src/AlgebraicAgents.jl @@ -58,6 +58,14 @@ export draw ## save and load export save, load +# Annotate and display the "input" and "output" of each agent in a hierarchy, +# similar to how "wires" connect the agents. +include("wires.jl") +export get_wires_from, get_wires_to +export add_wire!, delete_wires! +export retrieve_input_vars +export wiring_diagram, run_graphviz + # convenient agent subtyping include("agents.jl") export @aagent @@ -78,6 +86,6 @@ export @f_str, @filter, filter export TransformQuery export @transform, transform -include("integrations/requires.jl") +include("requires.jl") end diff --git a/src/integrations/AgentsIntegration/core.jl b/src/integrations/AgentsIntegration/core.jl index 17f86f3..3023bbc 100644 --- a/src/integrations/AgentsIntegration/core.jl +++ b/src/integrations/AgentsIntegration/core.jl @@ -13,8 +13,7 @@ Initialize `ABMAgent`, incl. hierarchy of ABM's agents. Configure the evolutionary step, logging, and step size by keyword arguments below. # Arguments - - `agent_step!`, `model_step!`: same meaning as in `Agents.step!` - - in general, any kwarg accepted by `Agents.run!`, incl. `adata`, `mdata` + - any kwarg accepted by `Agents.run!`, incl. `adata`, `mdata` - `when`, `when_model`: when to collect agents data, model data true by default, and performs data collection at every step if an `AbstractVector`, checks if `t ∈ when`; otherwise a function (model, t) -> ::Bool @@ -34,8 +33,6 @@ mutable struct ABMAgent <: AbstractAlgebraicAgent abm::Agents.AgentBasedModel - agent_step!::Any - model_step!::Any # evolutionary functions kwargs::Any # kwargs propagated to `run!` (incl. `adata`, `mdata`) when::Any when_model::Any # when to collect agents data, model data @@ -54,7 +51,6 @@ mutable struct ABMAgent <: AbstractAlgebraicAgent ## implement constructor function ABMAgent(name::AbstractString, abm::Agents.AgentBasedModel; - agent_step! = Agents.dummystep, model_step! = Agents.dummystep, when = true, when_model = when, step_size = 1.0, tspan::NTuple{2, Float64} = (0.0, Inf), kwargs...) @@ -63,8 +59,6 @@ mutable struct ABMAgent <: AbstractAlgebraicAgent setup_agent!(i, name) i.abm = abm - i.agent_step! = agent_step! - i.model_step! = model_step! i.kwargs = kwargs i.when = when i.when_model = when_model @@ -75,12 +69,12 @@ mutable struct ABMAgent <: AbstractAlgebraicAgent i.df_agents = DataFrames.DataFrame() i.df_model = DataFrames.DataFrame() - i.abm.properties[:__aagent__] = i + Agents.abmproperties(i.abm)[:__aagent__] = i i.abm0 = deepcopy(i.abm) i.t0 = i.t # initialize contained agents - for (id, _) in abm.agents + for id in Agents.allids(abm) entangle!(i, AAgent(string(id))) end @@ -102,34 +96,24 @@ function _step!(a::ABMAgent) collect_model = a.when_model isa AbstractVector ? (t ∈ a.when_model) : a.when isa Bool ? a.when : a.when_model(a.abm, t) - df_agents, df_model = Agents.run!(a.abm, a.agent_step!, a.model_step!, 1; - a.kwargs...) + df_agents, df_model = Agents.run!(a.abm, 1.0; a.kwargs...) + # append collected data ## df_agents - if collect_agents && ("step" ∈ names(df_agents)) + if collect_agents && ("time" ∈ names(df_agents)) if a.t == a.tspan[1] - df_agents_0 = df_agents[df_agents.step .== 0.0, :] - df_agents_0[!, :step] = convert.(Float64, df_agents_0[!, :step]) - df_agents_0[!, :step] .+= a.t - append!(a.df_agents, df_agents_0) + append!(a.df_agents, df_agents) + else + push!(a.df_agents, df_agents[end, :]) end - df_agents = df_agents[df_agents.step .== 1.0, :] - append!(a.df_agents, df_agents) - a.df_agents[(end - DataFrames.nrow(df_agents) + 1):end, :step] .+= a.t + - step_size - 1 end ## df_model - if collect_model && ("step" ∈ names(df_model)) + if collect_model && ("time" ∈ names(df_model)) if a.t == a.tspan[1] - df_model_0 = df_model[df_model.step .== 0.0, :] - df_model_0[!, :step] = convert.(Float64, df_model_0[!, :step]) - df_model_0[!, :step] .+= a.t - append!(a.df_model, df_model_0) + append!(a.df_model, df_model) + else + push!(a.df_model, df_model[end, :]) end - df_model = df_model[df_model.step .== 1.0, :] - append!(a.df_model, df_model) - a.df_model[(end - DataFrames.nrow(df_model) + 1):end, :step] .+= a.t + - step_size - 1 end a.t += step_size @@ -137,23 +121,23 @@ end # if step is a float, need to retype the dataframe function fix_float!(df, val) - if eltype(df[!, :step]) <: Int && !isa(val, Int) - df[!, :step] = convert.(Float64, df[!, :step]) + if eltype(df[!, :time]) <: Int && !isa(val, Int) + df[!, :time] = convert.(Float64, df[!, :time]) end end _projected_to(a::ABMAgent) = a.tspan[2] <= a.t ? true : a.t function getobservable(a::ABMAgent, obs) - getproperty(a.abm.properties, Symbol(obs)) + getproperty(abmproperties(a.abm), Symbol(obs)) end function gettimeobservable(a::ABMAgent, t::Float64, obs) df = a.df_model - @assert ("step" ∈ names(df)) && (string(obs) ∈ names(df)) + @assert ("time" ∈ names(df)) && (string(obs) ∈ names(df)) # query dataframe - df[df.step .== Int(t), obs] |> first + df[df.time .== Int(t), obs] |> first end function _reinit!(a::ABMAgent) @@ -174,7 +158,7 @@ Base.propertynames(::AAgent) = fieldnames(AAgent) ∪ [:agent] function Base.getproperty(a::AAgent, prop::Symbol) if prop == :agent - getparent(a).abm.agents[parse(Int, getname(a))] + (getparent(a).abm)[parse(Int, getname(a))] else getfield(a, prop) end @@ -190,10 +174,10 @@ end function gettimeobservable(a::AAgent, t::Float64, obs) df = getparent(a).df_agents - @assert ("step" ∈ names(df)) && (string(obs) ∈ names(df)) + @assert ("time" ∈ names(df)) && (string(obs) ∈ names(df)) # query df - df[(df.step .== Int(t)) .& (df.id .== a.agent.id), obs] |> first + df[(df.time .== Int(t)) .& (df.id .== a.agent.id), obs] |> first end function print_custom(io::IO, mime::MIME"text/plain", a::ABMAgent) @@ -220,7 +204,7 @@ end # retrieve algebraic agent as a property of the core dynamical system function extract_agent(model::Agents.ABM, agent::Agents.AbstractAgent) - model.properties[:__aagent__].inners[string(agent.id)] + abmproperties(model)[:__aagent__].inners[string(agent.id)] end """ @@ -233,7 +217,7 @@ algebraic_model = @get_model abm_model ``` """ macro get_model(model) - :($(esc(model)).properties[:__aagent__]) + :(abmproperties($(esc(model)))[:__aagent__]) end # macros to add, kill agents @@ -263,7 +247,7 @@ macro a(call) if call.args[1] == :add_agent! quote model = $(esc(model_call)) - omodel = model isa ABMAgent ? model : model.properties[:__aagent__] + omodel = model isa ABMAgent ? model : abmproperties(model)[:__aagent__] agent = $(esc(call)) entangle!(omodel, ABAModel(string(agent.id), a)) @@ -272,7 +256,7 @@ macro a(call) agent = model_call.args[2] quote model = $(esc(model_call)) - omodel = model isa ABMAgent ? model : model.properties[:__aagent__] + omodel = model isa ABMAgent ? model : abmproperties(model)[:__aagent__] agent = $(esc(agent)) agent = agent isa Number ? string(agent) : string(agent.id) diff --git a/src/integrations/SciMLIntegration/core.jl b/src/integrations/SciMLIntegration/core.jl index f39f8b4..cd1c96f 100644 --- a/src/integrations/SciMLIntegration/core.jl +++ b/src/integrations/SciMLIntegration/core.jl @@ -136,7 +136,7 @@ function print_custom(io::IO, mime::MIME"text/plain", a::DiffEqAgent) print_observables(IOContext(io, :indent => get(io, :indent, 0) + 3), mime, a) end -"Print in/out observables of a `DiffEqAgent`." +"""Print observables (positional indices and pretty names of "exported variables") of a `DiffEqAgent`.""" function print_observables(io::IO, ::MIME"text/plain", a::DiffEqAgent) indent = get(io, :indent, 0) diff --git a/src/opera.jl b/src/opera.jl index e1ca684..23d6043 100644 --- a/src/opera.jl +++ b/src/opera.jl @@ -12,6 +12,10 @@ const FutureLog = NamedTuple{(:id, :time, :retval), const Control = NamedTuple{(:id, :call), <:Tuple{AbstractString, Function}} const ControlLog = NamedTuple{(:id, :time, :retval), <:Tuple{AbstractString, Any, Any}} +## wire type +const Wire = NamedTuple{(:from, :from_var_name, :to, :to_var_name), + <:Tuple{AbstractAlgebraicAgent, Any, AbstractAlgebraicAgent, Any}} + """ Opera(uuid2agent_pairs...) A dynamic structure that @@ -96,6 +100,8 @@ mutable struct Opera controls::Vector{Control} controls_log::Vector{ControlLog} n_controls::Ref{UInt} + # wires + wires::Vector{Wire} function Opera(uuid2agent_pairs...) new(Dict{UUID, AbstractAlgebraicAgent}(uuid2agent_pairs...), @@ -107,7 +113,8 @@ mutable struct Opera 0, Vector{Control}(undef, 0), Vector{ControlLog}(undef, 0), - 0) + 0, + Vector{Wire}(undef, 0)) end end @@ -389,7 +396,7 @@ function execute_controls!(opera::Opera, time) end # if `expr` is a string, parse it as an expression -function get_expr(expr) +function get_expr(expr; eval_scope = @__MODULE__) if expr isa AbstractString Base.eval(eval_scope, Meta.parseall(expr)) else @@ -423,7 +430,7 @@ function load_opera!(opera::Opera, dump::AbstractDict; eval_scope = @__MODULE__) # instantious interactions for interaction in get(dump, "instantious", []) add_instantious!(opera, - get_expr(interaction["call"]), + get_expr(interaction["call"]; eval_scope), get(interaction, "priority", 0), get(interaction, "id", @@ -434,14 +441,14 @@ function load_opera!(opera::Opera, dump::AbstractDict; eval_scope = @__MODULE__) for interaction in get(dump, "futures", []) add_future!(opera, interaction["time"], - get_expr(interaction["call"]), + get_expr(interaction["call"]; eval_scope), get(interaction, "id", "future_" * get_count(opera, :n_futures))) end # controls for interaction in get(dump, "controls", []) add_control!(opera, - get_expr(interaction["call"]), + get_expr(interaction["call"]; eval_scope), get(interaction, "id", "control_" * get_count(opera, :n_controls))) end diff --git a/src/integrations/requires.jl b/src/requires.jl similarity index 63% rename from src/integrations/requires.jl rename to src/requires.jl index c22baba..06247a9 100644 --- a/src/integrations/requires.jl +++ b/src/requires.jl @@ -1,30 +1,35 @@ function __init__() # SciMLIntegration - @require DifferentialEquations="0c46a032-eb83-5123-abaf-570d42b7fbaa" include("SciMLIntegration/core.jl") + @require DifferentialEquations="0c46a032-eb83-5123-abaf-570d42b7fbaa" include("integrations/SciMLIntegration/core.jl") ## plotting @require DifferentialEquations="0c46a032-eb83-5123-abaf-570d42b7fbaa" begin - @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("SciMLIntegration/plots.jl") + @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("integrations/SciMLIntegration/plots.jl") end # AgentsIntegration - @require Agents="46ada45e-f475-11e8-01d0-f70cc89e6671" include("AgentsIntegration/core.jl") + @require Agents="46ada45e-f475-11e8-01d0-f70cc89e6671" include("integrations/AgentsIntegration/core.jl") ## plotting @require Agents="46ada45e-f475-11e8-01d0-f70cc89e6671" begin - @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("AgentsIntegration/plots.jl") + @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("integrations/AgentsIntegration/plots.jl") end # AlgebraicDynamicsIntegration - @require AlgebraicDynamics="5fd6ff03-a254-427e-8840-ba658f502e32" include("AlgebraicDynamicsIntegration/core.jl") + @require AlgebraicDynamics="5fd6ff03-a254-427e-8840-ba658f502e32" include("integrations/AlgebraicDynamicsIntegration/core.jl") ## SciML transforms @require AlgebraicDynamics="5fd6ff03-a254-427e-8840-ba658f502e32" begin - @require DifferentialEquations="0c46a032-eb83-5123-abaf-570d42b7fbaa" include("AlgebraicDynamicsIntegration/sciml_transform.jl") + @require DifferentialEquations="0c46a032-eb83-5123-abaf-570d42b7fbaa" include("integrations/AlgebraicDynamicsIntegration/sciml_transform.jl") end # DataFrame log out-of-the-box plots @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin - @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("utils_plots.jl") + @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("integrations/utils_plots.jl") + end + + # Use Graphviz_jll in `run_graphviz` + @require Graphviz_jll="3c863552-8265-54e4-a6dc-903eb78fde85" begin + include("wires_GraphvizExt.jl") end end diff --git a/src/wires.jl b/src/wires.jl new file mode 100644 index 0000000..a7e7056 --- /dev/null +++ b/src/wires.jl @@ -0,0 +1,262 @@ +""" + get_wires_to(a) +Get wires going into an agent. +""" +function get_wires_to(a::AbstractAlgebraicAgent) + return filter(w -> w.to == a, getopera(a).wires) +end + +""" + get_wires_from(a) +Get wires originating from an agent. +""" +function get_wires_from(a::AbstractAlgebraicAgent) + return filter(w -> w.from == a, getopera(a).wires) +end + +""" + retrieve_input_vars(a) +Return a dictionary with values along wires going into `a`, +specificed as `target => value` pairs. + +```julia +retrieve_input_vars(alice) # Dict("alice1_x" => "alice") +``` +""" +function retrieve_input_vars(a::AbstractAlgebraicAgent) + wires = filter(w -> (w.to == a) && !isnothing(w.from_var_name), getopera(a).wires) + + observables = map(wires) do w + w.to_var_name => getobservable(w.from, w.from_var_name) + end + + return Dict(observables...) +end + +# If `a` is a string or an UUID, retrieve the agent. +function to_agent(a, relative) + if !isa(relative, AbstractAlgebraicAgent) + return getagent(a, relative) + else + return relative + end +end + +""" + add_wire!(a; from, to, from_var_name, to_var_name) +Add a wire connecting two agents. + +# Examples +```julia +add_wire!(joint_system; from=alice, to=bob, from_var_name="alice_x", to_var_name="bob_x") +``` +""" +function add_wire!(a::AbstractAlgebraicAgent; + from::T, + to::T, + from_var_name = nothing, + to_var_name = nothing) where {T <: + Union{AbstractAlgebraicAgent, AbstractString, UUID}} + from, to = to_agent(a, from), to_agent(a, to) + from_var_name = something(from_var_name, to_var_name) + to_var_name = something(to_var_name, from_var_name) + + wire = (; from, from_var_name, to, to_var_name) + + return add_wire!(a, wire) +end + +add_wire!(a::AbstractAlgebraicAgent, w::Wire) = push!(getopera(a).wires, w) + +""" + delete_wires!(a; from, to, from_var_name, to_var_name) +Delete wires connecting two agents. Optionally, specify source and target variables. + +# Examples +```julia +delete_wires!(joint_system; from=alice, to=bob) +delete_wires!(joint_system; from=alice, to=bob, from_var_name="alice_x", to_var_name="bob_x") +``` +""" +function delete_wires!(a::AbstractAlgebraicAgent; + from::T, + to::T, + from_var_name = nothing, + to_var_name = nothing) where {T <: + Union{AbstractAlgebraicAgent, AbstractString, UUID}} + from, to = to_agent(a, from), to_agent(a, to) + + ixs = findall(x -> x.from == from && x.to == to, getopera(a).wires) + + if !isnothing(from_var_name) + ixs = ixs ∩ findall(x -> x.from_var_name == from_var_name, getopera(a).wires) + end + + if !isnothing(to_var_name) + ixs = ixs ∩ findall(x -> x.to_var_name == to_var_name, getopera(a).wires) + end + + foreach(ix -> deleteat!(getopera(a).wires, ix), ixs) + + return getopera(a).wires +end + +""" + wiring_diagram(agent; parentship_edges=true, wires=true) + wiring_diagram(agents; parentship_edges=true, wires=true) + wiring_diagram(groups; group_labels=nothing, parentship_edges=true, wires=true) + +Render a Graphviz graph of agents in an hierarchy. +By default, the graph shows edges between parent and child agents, +and annotated wires. + +Also see [`agent_hierarchy_mmd`](@ref). + +# Examples +```julia +# Build a compound problem. +joint_system = ⊕(alice, bob, name = "joint system") + +wiring_diagram(joint_system) + +# Do not show edges between parents and children. +wiring_diagram(joint_system; parentship_edges=false) + +# Only show listed agents. +wiring_diagram([alice, alice1, bob, bob1]) + +# Group agents into two clusters. +wiring_diagram([[alice, alice1], [bob, bob1]]) +# Provide labels for clusters. +wiring_diagram([[alice, alice1], [bob, bob1]]; group_labels=["alice", "bob"], parentship_edges=false) +``` +""" +function wiring_diagram end + +function wiring_diagram(agent::AbstractAlgebraicAgent; kwargs...) + all_agents = collect(values(getdirectory(agent))) + + return wiring_diagram(all_agents; kwargs...) +end + +function wiring_diagram(agents::Vector{T}; + parentship_edges = true, + wires = true) where {T <: AbstractAlgebraicAgent} + wiring_diagram([agents]; parentship_edges, wires) +end + +function wiring_diagram(groups; + group_labels = nothing, + parentship_edges = true, + wires = true) + nodes = build_nodes(groups; group_labels) + edges_parentship, edges_wires = build_edges(vcat(groups...)) + + edges = "" + if parentship_edges + edges *= "\n" * edges_parentship + end + + if wires + edges *= "\n" * edges_wires + end + + return """ + digraph "algagents" { + compound=true; + node[color=Teal, fontsize=7.0, width=0.5, height=0.5, shape=circle]; + $nodes + $edges + }""" +end + +# see https://github.com/AlgebraicJulia/Catlab.jl/blob/8d35ef724f0f0864ecdccef173eec958329f43e5/src/graphics/Graphviz.jl + +gv_backend(backend::Symbol, prog) = gv_backend(Val{backend}, prog) +gv_backend(::Type{<:Val}, prog) = prog + +""" + run_graphviz(io::IO, graph::AbstractString; prog::Symbol=:dot, format::String="svg") + run_graphviz(path::AbstractString, graph::AbstractString; prog::Symbol=:dot, format::String="svg") +Run the Graphviz program to render the graph and stream the results into `io`. + +This requires either `prog` (e.g., `dot`) to be available in your path (see https://graphviz.org) +or for the `Graphviz_jll` package to be installed and loaded before calling this function. + +See [`wiring_diagram`](@ref) to obtain the Graphviz wiring diagram for an agent hierarchy. +""" +function run_graphviz( + io::IO, graph::AbstractString; prog::Symbol = :dot, format::String = "svg") + prog = gv_backend(:graphviz_jll, prog) + + return open(gv -> print(gv, graph), `$prog -T$format`, io, write = true) +end + +function run_graphviz(path::AbstractString, graph::AbstractString; + prog::Symbol = :dot, format::String = "svg") + open(path, "w+") do io + run_graphviz(io, graph; prog, format) + end +end + +# Build code for nodes or subgraphs, in case multiple groups are specified. +function build_nodes(groups; group_labels = nothing) + if length(groups) == 1 + nodes = ["""$i [label="$(getname(a))"]""" for (i, a) in enumerate(only(groups))] + + return join(nodes, "\n") + else + subgraphs = [] + j = 0 + for (i_group, group) in enumerate(groups) + nodes = ["""$(j+i) [label="$(getname(a))"]""" for (i, a) in enumerate(group)] + + push!(subgraphs, + """ + subgraph cluster_$i_group {\n + """ * + (!isnothing(group_labels) ? """label="$(group_labels[i_group])" \n""" : + "") * + join(nodes, "\n") * + "\n}") + j += length(group) + end + + return join(subgraphs, "\n") + end +end + +# Build parentship edges and wires. +function build_edges(all_agents) + parents = map(all_agents) do a + parent = getparent(a) + while !isnothing(parent) && !in(parent, all_agents) + parent = getparent(parent) + end + + parent + end + + edges_parentship = [] + for (i, a) in enumerate(all_agents) + if !isnothing(parents[i]) + p = parents[i] + ix1, ix2 = findfirst(==(p), all_agents), findfirst(==(a), all_agents) + push!(edges_parentship, + "$ix1 -> $ix2 [len=1, penwidth=0.5, arrowsize=0.4, arrowtype=normal, style=dashed, fontsize=5.0, color=grey]") + end + end + + edges_wires = [] + for a in all_agents, b in all_agents + ix1, ix2 = findfirst(==(a), all_agents), findfirst(==(b), all_agents) + + oriented_wires_between = get_wires_from(a) ∩ get_wires_to(b) + for wire in oriented_wires_between + push!(edges_wires, + "$ix1 -> $ix2 [len=1, headlabel=$(wire.from_var_name), taillabel=$(wire.to_var_name), arrowsize=0.3, arrow=normal, fontsize=7.0]") + end + end + + return join(edges_parentship, "\n"), join(edges_wires, "\n") +end diff --git a/src/wires_GraphvizExt.jl b/src/wires_GraphvizExt.jl new file mode 100644 index 0000000..b743cf1 --- /dev/null +++ b/src/wires_GraphvizExt.jl @@ -0,0 +1,15 @@ +# see https://github.com/AlgebraicJulia/Catlab.jl/blob/8d35ef724f0f0864ecdccef173eec958329f43e5/ext/CatlabGraphvizExt.jl + +using .Graphviz_jll + +println("running") + +function gv_backend(::Type{Val{:graphviz_jll}}, prog) + getfield(Graphviz_jll, Symbol(prog))(identity) +end + +let cfg = joinpath(Graphviz_jll.artifact_dir, "lib", "graphviz", "config6") + if !isfile(cfg) + Graphviz_jll.dot(path -> run(`$path -c`)) + end +end diff --git a/test/integrations/agents_sir.jl b/test/integrations/agents_sir.jl new file mode 100644 index 0000000..14d8253 --- /dev/null +++ b/test/integrations/agents_sir.jl @@ -0,0 +1,148 @@ +# Taken from https://github.com/JuliaDynamics/Agents.jl/blob/dc2ce2c8b9e805e7c0f6b2ead4d120f0b1590ef9/src/models/sir.jl + +using LinearAlgebra +using StatsBase +using Random: Xoshiro + +@agent struct PoorSoul(GraphAgent) + days_infected::Int # number of days since is infected + status::Symbol # 1: S, 2: I, 3:R +end + +""" +```julia +sir(; + C = 8, + max_travel_rate = 0.01, + Ns = rand(50:5000, C), + β_und = rand(0.3:0.02:0.6, C), + β_det = β_und ./ 10, + infection_period = 30, + reinfection_probability = 0.05, + detection_time = 14, + death_rate = 0.02, + Is = [zeros(Int, length(Ns) - 1)..., 1], + seed = 19, +) +``` +Same as in [SIR model for the spread of COVID-19](@ref). +""" +function sir(; + C = 8, + max_travel_rate = 0.01, + Ns = rand(50:5000, C), + β_und = rand(0.3:0.02:0.6, C), + β_det = β_und ./ 10, + infection_period = 30, + reinfection_probability = 0.05, + detection_time = 14, + death_rate = 0.02, + Is = [zeros(Int, length(Ns) - 1)..., 1], + seed = 19 +) + rng = Xoshiro(seed) + migration_rates = zeros(C, C) + @assert length(Ns)== + length(Is)== + length(β_und)== + length(β_det)== + size(migration_rates, 1) "length of Ns, Is, and B, and number of rows/columns in migration_rates should be the same " + @assert size(migration_rates, 1)==size(migration_rates, 2) "migration_rates rates should be a square matrix" + + for c in 1:C + for c2 in 1:C + migration_rates[c, c2] = (Ns[c] + Ns[c2]) / Ns[c] + end + end + maxM = maximum(migration_rates) + migration_rates = (migration_rates .* max_travel_rate) ./ maxM + migration_rates[diagind(migration_rates)] .= 1.0 + + ## normalize migration_rates + migration_rates_sum = sum(migration_rates, dims = 2) + for c in 1:C + migration_rates[c, :] ./= migration_rates_sum[c] + end + + properties = Dict( + :Ns => Ns, + :Is => Is, + :β_und => β_und, + :β_det => β_det, + :migration_rates => migration_rates, + :infection_period => infection_period, + :infection_period => infection_period, + :reinfection_probability => reinfection_probability, + :detection_time => detection_time, + :C => C, + :death_rate => death_rate + ) + + space = GraphSpace(Agents.Graphs.complete_graph(C)) + model = ABM(PoorSoul, space; agent_step! = sir_agent_step!, properties, rng) + + ## Add initial individuals + for city in 1:C, n in 1:Ns[city] + ind = add_agent!(city, model, 0, :S) # Susceptible + end + ## add infected individuals + for city in 1:C + inds = ids_in_position(city, model) + for n in 1:Is[city] + agent = model[inds[n]] + agent.status = :I # Infected + agent.days_infected = 1 + end + end + return model, sir_agent_step!, dummystep +end + +function sir_agent_step!(agent, model) + sir_migrate!(agent, model) + sir_transmit!(agent, model) + sir_update!(agent, model) + sir_recover_or_die!(agent, model) +end + +function sir_migrate!(agent, model) + pid = agent.pos + m = sample(abmrng(model), 1:(model.C), Weights(model.migration_rates[pid, :])) + if m ≠ pid + move_agent!(agent, m, model) + end +end + +function sir_transmit!(agent, model) + agent.status == :S && return + rate = if agent.days_infected < model.detection_time + model.β_und[agent.pos] + else + model.β_det[agent.pos] + end + + n = rate * abs(randn(abmrng(model))) + n <= 0 && return + + for contactID in ids_in_position(agent, model) + contact = model[contactID] + if contact.status == :S || + (contact.status == :R && rand(abmrng(model)) ≤ model.reinfection_probability) + contact.status = :I + n -= 1 + n <= 0 && return + end + end +end + +sir_update!(agent, model) = agent.status == :I && (agent.days_infected += 1) + +function sir_recover_or_die!(agent, model) + if agent.days_infected ≥ model.infection_period + if rand(abmrng(model)) ≤ model.death_rate + remove_agent!(agent, model) + else + agent.status = :R + agent.days_infected = 0 + end + end +end diff --git a/test/integrations/agents_test.jl b/test/integrations/agents_test.jl index 169c376..312f8fb 100644 --- a/test/integrations/agents_test.jl +++ b/test/integrations/agents_test.jl @@ -3,13 +3,15 @@ using Plots import Random import StatsBase: sample, Weights +include("agents_sir.jl") + # test pure Agents.jl solution vs AlgebraicAgents.jl wrap # Agents.jl Random.seed!(2023) # use Agents.jl predefined model, in https://juliadynamics.github.io/Agents.jl/stable/models/#Predefined-Models-1 -abm_agents, agent_step, _ = Agents.Models.sir() +abm_agents, agent_step, _ = sir() # data to collect infected(x) = count(i == :I for i in x) @@ -20,70 +22,27 @@ to_collect = [(:status, f) for f in (infected, recovered, length)] Random.seed!(2023) # use Agents.jl predefined model, in https://juliadynamics.github.io/Agents.jl/stable/models/#Predefined-Models-1 -abm_algebraic, _, _ = Agents.Models.sir() +abm_algebraic, _, _ = sir() # modify stepping functions function agent_step!(agent, model) @get_model model extract_agent(model, agent) - migrate!(agent, model) - transmit!(agent, model) - update!(agent, model) - recover_or_die!(agent, model) -end - -function migrate!(agent, model) - pid = agent.pos - m = sample(model.rng, 1:(model.C), Weights(model.migration_rates[pid, :])) - if m ≠ pid - move_agent!(agent, m, model) - end -end - -function transmit!(agent, model) - agent.status == :S && return - rate = if agent.days_infected < model.detection_time - model.β_und[agent.pos] - else - model.β_det[agent.pos] - end - - n = rate * abs(randn(model.rng)) - n <= 0 && return - - for contactID in ids_in_position(agent, model) - contact = model[contactID] - if contact.status == :S || - (contact.status == :R && rand(model.rng) ≤ model.reinfection_probability) - contact.status = :I - n -= 1 - n <= 0 && return - end - end -end - -update!(agent, model) = agent.status == :I && (agent.days_infected += 1) - -function recover_or_die!(agent, model) - if agent.days_infected ≥ model.infection_period - if rand(model.rng) ≤ model.death_rate - @a kill_agent!(agent, model) - else - agent.status = :R - agent.days_infected = 0 - end - end + sir_migrate!(agent, model) + sir_transmit!(agent, model) + sir_update!(agent, model) + sir_recover_or_die!(agent, model) end @testset "Agents.jl and AlgebraicAgents.jl solution are equal" begin Random.seed!(1) - abm_algebraic_wrap = ABMAgent("sir_model", abm_algebraic; agent_step!, + abm_algebraic_wrap = ABMAgent("sir_model", abm_algebraic; tspan = (0.0, 10.0), adata = to_collect) simulate(abm_algebraic_wrap) data_algebraic = abm_algebraic_wrap.df_agents Random.seed!(1) - data_agent, _ = run!(abm_agents, agent_step, 10; adata = to_collect) + data_agent, _ = run!(abm_agents, 10; adata = to_collect) @test abm_algebraic_wrap.t == 10.0 @test data_algebraic == data_agent @@ -92,7 +51,7 @@ end end @testset "plotting for ABM wraps" begin - abm_algebraic_wrap = ABMAgent("sir_model", abm_algebraic; agent_step!, + abm_algebraic_wrap = ABMAgent("sir_model", abm_algebraic; tspan = (0.0, 10.0), adata = to_collect) simulate(abm_algebraic_wrap) diff --git a/test/opera.jl b/test/opera.jl index d571c12..3ac5994 100644 --- a/test/opera.jl +++ b/test/opera.jl @@ -188,3 +188,52 @@ end @test opera.controls_log[1].retval == "alice" @test opera.controls_log[2].retval == 1.0 end + +@testset "wires" begin + @aagent struct MyAgent4 end + + alice = MyAgent4("alice") + alice1 = MyAgent4("alice1") + entangle!(alice, alice1) + + bob = MyAgent4("bob") + bob1 = MyAgent4("bob1") + entangle!(bob, bob1) + + joint_system = ⊕(alice, bob, name = "joint system") + + # Add wires. + add_wire!(joint_system; + from = alice, + to = bob, + from_var_name = "alice_x", + to_var_name = "bob_x") + add_wire!(joint_system; + from = bob, + to = alice, + from_var_name = "bob_y", + to_var_name = "alice_y") + + add_wire!(joint_system; + from = alice, + to = alice1, + from_var_name = "alice_x", + to_var_name = "alice1_x") + add_wire!(joint_system; + from = bob, + to = bob1, + from_var_name = "bob_x", + to_var_name = "bob1_x") + + # Show wires. + @test length(get_wires_from(alice)) == 2 + @test length(get_wires_to(alice1)) == 1 + + # Retrieve variables along input wires. + AlgebraicAgents.getobservable(a::MyAgent4, args...) = getname(a) + + @test retrieve_input_vars(alice1) == Dict("alice1_x" => "alice") + + # Delete wires. + @test length(delete_wires!(joint_system; from = alice, to = alice1)) == 3 +end diff --git a/tutorials/agents/agents.jl b/tutorials/agents/agents.jl index 27d5d53..c73eeb7 100644 --- a/tutorials/agents/agents.jl +++ b/tutorials/agents/agents.jl @@ -15,7 +15,7 @@ using Distributions: Poisson, DiscreteNonParametric using DrWatson: @dict using Plots -@agent PoorSoul GraphAgent begin +@agent struct PoorSoul(GraphAgent) days_infected::Int ## number of days since is infected status::Symbol ## 1: S, 2: I, 3:R end @@ -61,7 +61,7 @@ function model_initiation(; C, death_rate) space = GraphSpace(complete_digraph(C)) - model = ABM(PoorSoul, space; properties, rng) + model = ABM(PoorSoul, space; properties, rng, model_step! = identity) ## Add initial individuals for city in 1:C, n in 1:Ns[city] @@ -193,7 +193,7 @@ to_collect = [(:status, f) for f in (infected, recovered, length)] # We wrap the model as an agent: -m = ABMAgent("sir_model", abm; agent_step!, tspan=(0., 100.), adata=to_collect) +m = ABMAgent("sir_model", abm; tspan=(0., 100.), adata=to_collect) # And we simulate the dynamics: diff --git a/tutorials/wires/Project.toml b/tutorials/wires/Project.toml new file mode 100644 index 0000000..b0eeac4 --- /dev/null +++ b/tutorials/wires/Project.toml @@ -0,0 +1,3 @@ +[deps] +AlgebraicAgents = "f6eb0ae3-10fa-40e6-88dd-9006ba45093a" +Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85" diff --git a/tutorials/wires/wires.jl b/tutorials/wires/wires.jl new file mode 100644 index 0000000..905f0de --- /dev/null +++ b/tutorials/wires/wires.jl @@ -0,0 +1,49 @@ +using AlgebraicAgents + +# Define a simple hierarchy. + +@aagent struct MyAgent end + +alice = MyAgent("alice") +alice1 = MyAgent("alice1") +entangle!(alice, alice1) + +bob = MyAgent("bob") +bob1 = MyAgent("bob1") +entangle!(bob, bob1) + +joint_system = ⊕(alice, bob, name = "joint system") + +# Add wires. +add_wire!(joint_system; from=alice, to=bob, from_var_name="alice_x", to_var_name="bob_x") +add_wire!(joint_system; from=bob, to=alice, from_var_name="bob_y", to_var_name="alice_y") + +add_wire!(joint_system; from=alice, to=alice1, from_var_name="alice_x", to_var_name="alice1_x") +add_wire!(joint_system; from=bob, to=bob1, from_var_name="bob_x", to_var_name="bob1_x") + +# Show wires. +get_wires_from(alice) +get_wires_to(alice1) + +# Retrieve variables along input wires. +AlgebraicAgents.getobservable(a::MyAgent, args...) = getname(a) + +retrieve_input_vars(alice1) + +# Plot wires. +graph1 = wiring_diagram(joint_system) + +graph2 = wiring_diagram(joint_system; parentship_edges=false) + +graph3 = wiring_diagram([alice, alice1, bob, bob1]) + +graph4 = wiring_diagram([[alice, alice1], [bob, bob1]]) + +graph5 = wiring_diagram([[alice, alice1], [bob, bob1]]; group_labels=["alice", "bob"], parentship_edges=false) + +# using Graphviz_jll + +run_graphviz("graph5.svg", graph5) + +# Delete wires. +delete_wires!(joint_system; from=alice, to=alice1) \ No newline at end of file