diff --git a/.gitignore b/.gitignore index 3be7528..d15634d 100644 --- a/.gitignore +++ b/.gitignore @@ -27,3 +27,7 @@ Manifest.toml .DS_Store .vscode .ipynb_checkpoints + +cache/ + +figures/SIRSDiscretetrajectory.png diff --git a/Project.toml b/Project.toml index e1ba5c7..3b36cc3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,19 @@ -name = "AlgebraicABMS" -uuid = "1301f41b-2c8a-4430-9d85-ab97d23a62cd" +name = "AlgebraicABMs" +uuid = "5a5e3447-9604-46e6-8d91-cb86f5f51721" license = "MIT" -authors = ["AlgebraicJulia Developer "] +authors = ["Kris Brown "] version = "0.0.1" [deps] +AlgebraicPetri = "4f99eebe-17bf-4e98-b6a1-2c4f205a959b" +AlgebraicRewriting = "725a01d3-f174-5bbd-84e1-b9417bad95d9" Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Fleck = "5bb9b785-358c-4fee-af0f-b94a146244a8" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] Catlab = "^0.16" diff --git a/README.md b/README.md index 775c574..e79d900 100644 --- a/README.md +++ b/README.md @@ -1,42 +1,37 @@ -# AlgebraicABMS.jl +# AlgebraicABMs.jl -[![Stable Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://AlgebraicJulia.github.io/AlgebraicABMS.jl/stable) -[![Development Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://AlgebraicJulia.github.io/AlgebraicABMS.jl/dev) -[![Code Coverage](https://codecov.io/gh/AlgebraicJulia/AlgebraicABMS.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/AlgebraicJulia/AlgebraicABMS.jl) -[![CI/CD](https://github.com/AlgebraicJulia/AlgebraicABMS.jl/actions/workflows/julia_ci.yml/badge.svg)](https://github.com/AlgebraicJulia/AlgebraicABMS.jl/actions/workflows/julia_ci.yml) +[![Stable Documentation](https://img.shields.io/badge/docs-stable-blue.svg)](https://AlgebraicJulia.github.io/AlgebraicABMs.jl/stable) +[![Development Documentation](https://img.shields.io/badge/docs-dev-blue.svg)](https://AlgebraicJulia.github.io/AlgebraicABMs.jl/dev) +[![Code Coverage](https://codecov.io/gh/AlgebraicJulia/AlgebraicABMs.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/AlgebraicJulia/AlgebraicABMs.jl) +[![CI/CD](https://github.com/AlgebraicJulia/AlgebraicABMs.jl/actions/workflows/julia_ci.yml/badge.svg)](https://github.com/AlgebraicJulia/AlgebraicABMs.jl/actions/workflows/julia_ci.yml) -A template repository for making a new AlgebraicJulia package. +Important example files: `docs/literate/petri_example.jl` and `test/ABMs.jl` -## 🛠️ Usage +## Caveats + +[Fleck.jl](https://github.com/adolgert/Fleck.jl) is not officially released, so you must clone that repo and `dev` it. -1. Use the "Use this template" dropdown to select "Create a new repository" -2. In the new page select "AlgebraicJulia" as the owner, give the repository a name, and create a new repository from the template -3. Set up Codecov credentials for code coverage (If you have trouble, reach out to an AlgebraicJulia organization owner to help with this) +This will need to be addressed before the tests can be run / documentation can be built automatically by Github actions. + +## 🛠️ Usage - 1. Log into [Codecov](https://codecov.io) with your GitHub account (this requires that you are a member of the AlgebraicJulia organization) - 2. Navigate to the [AlgebraicJulia organization](https://app.codecov.io/gh/AlgebraicJulia) - 3. Select your new repository from the list (e.x. "AlgebraicX") - 4. Note down the `CODECOV_TOKEN` value (It may be in the "Settings" tab if it doesn't show up immediately) - 5. Navigate back to your new GitHub repository and go to the Settings tab - 6. Go to "Security", "Secrets and variables", and "Actions" and click the "New repository secret" button - 7. Give the secret name `CODECOV_TOKEN` and the Secret value is the value you noted from the Codecov settings - 8. Click "Add secret" +To locally build the documentation and the literate code examples, run the following in the command line: +``` +julia --project=docs -e "using AlgebraicABMs, LiveServer; servedocs(literate_dir=\"docs/literate\",skip_dir=\"docs/src/generated\")" +``` -4. Clone the new repository, for example in the terminal: - ```sh - git clone https://github.com/AlgebraicJulia/AlgebraicX.jl.git - cd AlgebraicX.jl - ``` -5. Inspect for yourself and run `init.sh` with the new repository name and (optional) UUID are parameters. This script will substitute all instances of `AlgebraicX` with your new repository name and the default UUID with a new one or, if available, the UUID provided. -6. Go back to your repository and wait until the tests have passed, you can check the status by going to the "Actions" tab in the repository +To locally run the test suite, run the following command +``` +julia --project=test test/runtests.jl +``` -### Buildkite +### To-do: Buildkite AlgebraicJulia uses [Buildkite](https://buildkite.com/) to submit resource-intensive processes such as building documentation and executing tests to the [HiPerGator](https://www.rc.ufl.edu/about/hipergator/) computing cluster. While this template comes with a preconfigured `.buildkite/pipeline.yml` file, this repository is not integrated with Buildkite by default. If you would like your repository to use Buildkite to run processes on HiPerGator, tag an issue with @AlgebraicJulia/SysAdmins. -### 📔 Set Up GitHub Pages (Public Repos Only) +### 📔 To-do: Set Up GitHub Pages (Public Repos Only) 1. Follow the Usage steps above to set up a new template, make sure all initial GitHub Actions have passed 2. Navigate to the repository settings and go to "Code and automation", "Pages" @@ -45,21 +40,3 @@ While this template comes with a preconfigured `.buildkite/pipeline.yml` file, t 5. Go back to the main page of your repository and click the gear to the right of the "About" section in the right side column 6. Under "Website" check the checkbox that says "Use your GitHub Pages website" and click "Save changes" 7. You will now see a URL in the "About" section that will link to your package's documentation - -### 🛡️ Set Up Branch Protection (Public Repos Only) - -1. Follow the Usage steps above to set up a new template, make sure all initial GitHub Actions have passed -2. Navigate to the repository settings and go to "Code and automation", "Branches" -3. Click "Add branch protection rule" to start adding branch protection -4. Under "Branch name pattern" put `main`, this will add protection to the main branch -5. Make sure to set the following options: - - Check the "Require a pull request before merging" - - Check the "Request status checks to pass before merging" and make sure the following status checks are added to the required list: - - CI / Documentation - - CI / Julia 1 - ubuntu-latest - x64 - push - - CI / Julia 1 - ubuntu-latest - x86 - push - - CI / Julia 1 - windows-latest - x64 - push - - CI / Julia 1 - windows-latest - x86 - push - - CI / Julia 1 - macOS-latest - x64 - push - - Check the "Restrict who can push to matching branches" and add `algebraicjuliabot` to the list of people with push access -6. Click "Save changes" to enable the branch protection diff --git a/docs/Project.toml b/docs/Project.toml index 3e617d2..18eb94c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,11 @@ [deps] -AlgebraicABMS = "1301f41b-2c8a-4430-9d85-ab97d23a62cd" +AlgebraicABMs = "5a5e3447-9604-46e6-8d91-cb86f5f51721" +AlgebraicPetri = "4f99eebe-17bf-4e98-b6a1-2c4f205a959b" +AlgebraicRewriting = "725a01d3-f174-5bbd-84e1-b9417bad95d9" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" diff --git a/docs/literate/game_of_life.jl b/docs/literate/game_of_life.jl new file mode 100644 index 0000000..3ed5492 --- /dev/null +++ b/docs/literate/game_of_life.jl @@ -0,0 +1,41 @@ +# # Petri Net rewriting +# +# First we want to load our package with `using` + +using AlgebraicABMs +using Catlab, AlgebraicPetri # for declaring model building blocks +using Distributions # for defining hazard rates + +# ## Schema +# +# We create a regular grid + +@present SchLife(FreeSchema) begin + Cell::Ob + (N,E,W,S)::Hom(Cell,Cell) + Life::Ob + live::Hom(Life,Cell) +end + +# We define a network of cells + +@present SchLifeGraph <: SchSymmetricGraph begin + Life::Ob + live::Hom(Life,V) +end + +@acset_type Life(SchLifeGraph) <: AbstractSymmetricGraph + +function living_neighbors(n::Int; alive=false) + X = Life(1) + alive && add_part!(X, :Live, live=1) + for _ in 1:n + v = add_part!(X, :V) + add_part!(X, :Life, live=v) + add_edge!(X, v, 1) + end + X +end + +living_neighbors(2) + diff --git a/docs/literate/literate_example.jl b/docs/literate/literate_example.jl deleted file mode 100644 index 7f5496f..0000000 --- a/docs/literate/literate_example.jl +++ /dev/null @@ -1,13 +0,0 @@ -# # Code Example -# -# This is an example of adding a code example compiled with Literate.jl in the docs. -# -# First we want to load our package with `using` - -using AlgebraicABMS - -# ## Using `hello()` -# -# We provide the `hello(string)` method which prints "Hello, `string`!" - -hello("World") diff --git a/docs/literate/petri_example.jl b/docs/literate/petri_example.jl new file mode 100644 index 0000000..f50139c --- /dev/null +++ b/docs/literate/petri_example.jl @@ -0,0 +1,68 @@ +# # Petri Net rewriting +# +# First we want to load our package with `using` + +using AlgebraicABMs +using Catlab, AlgebraicPetri # for declaring model building blocks +using Distributions # for defining hazard rates +using Makie, CairoMakie # visualization + +# ## Petri-net based model +# +# We define an SIRS model with birth and death + +sir_pn= @acset LabelledPetriNet begin + S=3; sname=[:S,:I,:R] + T=7; tname=[:inf,:rec,:birth,:deathS,:deathI,:deathR,:wane] + I=7; it=[1,1,2,4,5,6,7]; is=[1,2,2,1,2,3,3] + O=5; ot=[1,1,2,3,7]; os=[2,2,3,1,1] +end + +to_graphviz(sir_pn) + +# A "state of the world" in this model is just a finite set of Susceptible, Infected, and Recovered people. Thus we can specify a model with 3 integers. + +init = PetriNetCSet(sir_pn, S=100, I=5); # Initial state + + +# We declare parameters to specify the random waiting times +pS, pI, pR = 95,5,0 +pop = 100 +lifespan = 65*365 +μ = 1/lifespan +β = 0.001 +wane = 60; + +# We associate with each transition a function which takes a time point and return a distribution of waiting times +clockdists = Dict{Symbol,Function}(); + +## the Exponential clocks (Markov) +clockdists[:inf] = (t) -> Exponential(1 / β) +clockdists[:birth] = (t) -> Exponential(1 / (μ*pop)) +clockdists[:deathS] = (t) -> Exponential(1 / μ) +clockdists[:deathI] = (t) -> Exponential(1 / μ) +clockdists[:deathR] = (t) -> Exponential(1 / μ) +clockdists[:wane] = (t) -> Exponential(wane) + +## the Weibull clock (non-Markov) +α, θ = weibullpar(30, 5) +clockdists[:rec] = (t) -> Weibull(α, θ); + +# We make a reporting function that extracts information from each time step. +count(acs::ACSet) = + NamedTuple(Dict([o => nparts(acs, o) for o in ob(acset_schema(acs))])); + + +# We simulate a model +sirout = run!(sir_pn, clockdists, init; save=count, maxevent=2000) +X = first.(sirout) +SIR = [getindex.(last.(sirout), x) for x in 1:3] +f = Figure(); +Legend(f[1, 2], lines!.(Ref(Axis(f[1,1])), Ref(X), SIR), ["S", "I","R"]) + +## We can also save the figure to the filesystem. +isdir("figures") || mkdir("figures"); +Makie.save("figures/SIRSDiscretetrajectory.png", f, px_per_unit=1, size=(800,600)); + +f + diff --git a/docs/make.jl b/docs/make.jl index 501aeb5..0018e6a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,11 +1,11 @@ using Documenter using Literate -const literate_dir = joinpath(@__DIR__, "literate") -const generated_dir = joinpath(@__DIR__, "src", "generated") +const LITERATE_INPUT = joinpath(@__DIR__, "literate") +const LITERATE_OUTPUT = joinpath(@__DIR__, "src", "generated") -@info "Loading AlgebraicABMS" -using AlgebraicABMS +@info "Loading AlgebraicABMs" +using AlgebraicABMs const no_literate = "--no-literate" in ARGS if !no_literate @@ -14,35 +14,47 @@ if !no_literate # Set Literate.jl config if not being compiled on recognized service. config = Dict{String,String}() if !(haskey(ENV, "GITHUB_ACTIONS") || haskey(ENV, "GITLAB_CI")) - config["nbviewer_root_url"] = "https://nbviewer.jupyter.org/github/AlgebraicJulia/AlgebraicABMS.jl/blob/gh-pages/dev" - config["repo_root_url"] = "https://github.com/AlgebraicJulia/AlgebraicABMS.jl/blob/main/docs" + config["nbviewer_root_url"] = "https://nbviewer.jupyter.org/github/AlgebraicJulia/AlgebraicABMs.jl/blob/gh-pages/dev" + config["repo_root_url"] = "https://github.com/AlgebraicJulia/AlgebraicABMs.jl/blob/main/docs" end - for (root, dirs, files) in walkdir(literate_dir) - out_dir = joinpath(generated_dir, relpath(root, literate_dir)) - for file in files - f, l = splitext(file) - if l == ".jl" && !startswith(f, "_") - Literate.markdown(joinpath(root, file), out_dir; - config=config, documenter=true, credit=false) - Literate.notebook(joinpath(root, file), out_dir; - execute=true, documenter=true, credit=false) - end - end + # for (root, dirs, files) in walkdir(literate_dir) + # out_dir = joinpath(generated_dir, relpath(root, literate_dir)) + # for file in files + # f, l = splitext(file) + # if l == ".jl" && !startswith(f, "_") + # Literate.markdown(joinpath(root, file), out_dir; + # config=config, documenter=true, credit=false) + # Literate.notebook(joinpath(root, file), out_dir; + # execute=true, documenter=true, credit=false) + # end + # end + # end + + for (root, _, files) ∈ walkdir(LITERATE_INPUT), file ∈ files + # ignore non julia files + splitext(file)[2] == ".jl" || continue + # full path to a literate script + ipath = joinpath(root, file) + # generated output path + opath = splitdir(replace(ipath, LITERATE_INPUT=>LITERATE_OUTPUT))[1] + # generate the markdown file calling Literate + Literate.markdown(ipath, opath) end + end @info "Building Documenter.jl docs" makedocs( - modules=[AlgebraicABMS], + modules=[AlgebraicABMs], format=Documenter.HTML(), - sitename="AlgebraicABMS.jl", + sitename="AlgebraicABMs.jl", doctest=false, checkdocs=:none, pages=Any[ - "AlgebraicABMS.jl"=>"index.md", + "AlgebraicABMs.jl"=>"index.md", "Examples"=>Any[ - "generated/literate_example.md", + "generated/petri_example.md", ], "Library Reference"=>"api.md", ] @@ -51,6 +63,6 @@ makedocs( @info "Deploying docs" deploydocs( target="build", - repo="github.com/AlgebraicJulia/AlgebraicABMS.jl.git", + repo="github.com/AlgebraicJulia/AlgebraicABMs.jl.git", branch="gh-pages" ) diff --git a/docs/src/api.md b/docs/src/api.md index f17ac0e..283c4ce 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -1,5 +1,9 @@ # Library Reference ```@autodocs -Modules = [AlgebraicABMS] +Modules = [ + AlgebraicABMs.Distributions, + AlgebraicABMs.PetriInterface, + AlgebraicABMs.RewriteSemiMarkov +] ``` diff --git a/docs/src/index.md b/docs/src/index.md index 7c1c1f2..6432026 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,7 +1,10 @@ -# AlgebraicABMS.jl +# AlgebraicABMs.jl ```@meta -CurrentModule = AlgebraicABMS +CurrentModule = AlgebraicABMs ``` -`AlgebraicABMS.jl` is a Julia library for... +`AlgebraicABMs.jl` is a Julia library for creating agent-based models. We ultimate want to provide capabilities on par with software like AnyLogic, NetLogo, and Agents.jl - all while offering a mostly declarative interface such that: +- The model's logic is transparent (rather than hidden away in complicated interactions of code blocks) +- Models can be built compositionally (we can naturally 'glue' models together at a high level, without worrying about implementation details and edge-cases) +- Models can be _migrated_ at a high level, whether interpersonally (collaboration with others who have a different ontology/vocabulary) or intrapersonally (one updates one's own model of the of world and wishes to reuse one's old model under new assumptions, without having to manually refactor code and dig into implementation details). diff --git a/init.sh b/init.sh index db99e1a..f5ead58 100755 --- a/init.sh +++ b/init.sh @@ -1,7 +1,7 @@ #!/bin/bash -DEFAULT_REPO='AlgebraicABMS' -DEFAULT_UUID='1301f41b-2c8a-4430-9d85-ab97d23a62cd' +DEFAULT_REPO='AlgebraicABMs' +DEFAULT_UUID='5A5E3447-9604-46E6-8D91-CB86F5F51721' usage="This script is for initializing the template with the new repository name and UUID. Please provide the new repository name and UUID in that order. The repository name cannot be 'Test.'\n Example:\n diff --git a/src/ABMs.jl b/src/ABMs.jl new file mode 100644 index 0000000..fac5eda --- /dev/null +++ b/src/ABMs.jl @@ -0,0 +1,384 @@ + +module ABMs + +export ABM + +using Distributions, Fleck, Random +using DataStructures: DefaultDict + +using Catlab, AlgebraicRewriting +using AlgebraicPetri: AbstractReactionNet +using AlgebraicRewriting.Incremental: connected_acset_components, key_dict +using AlgebraicRewriting.Rewrite.Migration: pres_hash +import Catlab: acset_schema, right, is_isomorphic +import ..PetriInterface: run! +import AlgebraicRewriting: get_match +import AlgebraicRewriting.Incremental: addition!, deletion! + +# Possibly upstream +################### +is_isomorphic(f::FinFunction) = is_monic(f) && is_epic(f) + +pattern(r::Rule) = codom(left(r)) + +acset_schema(r::Rule) = acset_schema(pattern(r)) + +Base.pairs(h::IncHomSet) = [k => h[k] for k in keys(key_dict(h))] + +""" +Extract data in representable cache as a dictionary. +This is an intermediate step in yoneda_cache that should be factored out. +""" +function repr_dict(T::Type, S=nothing; cache="cache")::Dict{Symbol, Tuple{ACSet, Int}} + S = Presentation(isnothing(S) ? T : S) + yoneda_cache(T, S; cache) + cache_dir = joinpath(cache, "$(nameof(T))_$(pres_hash(S))") + Dict(map(nameof.(generators(S, :Ob))) do name + path, ipath = joinpath.(cache_dir, ["$name.json", "_id_$name.json"]) + name => (read_json_acset(T, path), parse(Int,open(io->read(io, String), ipath))) + end) +end + + +# Timers +######## + +""" +Something that can produce a ACSetTransformation × clocktime → hazard_rate +""" +abstract type AbsTimer end +abstract type StateDependentTimer <: AbsTimer end +state_dep(t::AbsTimer) = t isa StateDependentTimer + +""" +A closure which accepts a ACSetTransformation and returns a function of type +clocktime → hazard_rate +""" +struct FullClosure <: StateDependentTimer + val::Function # ACSetTransformation → clocktime → hazard_rate +end +(c::FullClosure)(m::ACSetTransformation, t::Float64) = c.val(m,t) + +""" +A closure which accepts a clocktime and returns a hazard_rate. This is a timer +which cannot depend on the match data nor ACSet state. +""" +struct ClosureTime <: AbsTimer + val::Function # clocktime → hazard_rate +end +(c::ClosureTime)(t::Float64) = c.val(t) + +""" +A closure which accepts a match morphism and returns a hazard_rate. This is a +timer which cannot depend on the absolute clock time. +""" +struct ClosureState <: StateDependentTimer + val::Function # clocktime → hazard_rate +end +(c::ClosureState)(m::ACSetTransformation) = c.val(m) + +abstract type AbsHazard <: AbsTimer end + +struct DiscreteHazard <: AbsHazard + val::Distribution{Univariate, Discrete} +end + +DiscreteHazard(t::Number) = DiscreteHazard(Dirac(t)) + +struct ContinuousHazard <: AbsHazard + val::Distribution{Univariate, Continuous} +end + +ContinuousHazard(p::Number) = ContinuousHazard(Exponential(p)) + +get_hazard(m::ACSetTransformation, t::Float64, h::FullClosure) = h(m,t) +get_hazard(::ACSetTransformation, t::Float64, h::ClosureTime) = h(t) +get_hazard(m::ACSetTransformation, ::Float64, h::ClosureState) = h(m) +get_hazard(::ACSetTransformation, ::Float64, h::AbsHazard) = h.val + +# Rules +####### +abstract type PatternType end + +"""Empty patterns have (one) trivial pattern match""" +struct EmptyP <: PatternType end + +""" +Default case, where pattern matches should be found via (incremental) +homomorphism search and represented explicitly, each with own events getting +scheduled. +""" +struct RegularP <: PatternType end + +# """ +# Special case of homsearch where no backtracking is needed. The only nonempty +# sets in L are those for objects with no outgoing homs. There may be attributes, +# however, so at runtime we must filter the sets before picking random elements. +# E.g. for labeled set L = {:a, :a, AttrVar(1)} we randomly pick two elements with +# label :a and one arbitrary element. + +# WARNING: this is only viable if the timer associated with the rewrite rule is +# symmteric with respect to the discrete parts. +# """ +# struct DiscreteP <: PatternType +# parts::Dict{Symbol, Int} +# end + +""" +A pattern match from a coproduct of representables is just a choice of parts +in the codomain. E.g. matching L = •→• • • is just a random choice of edge and +two random vertices. + +The vector of ints refers to parts of L which are the counits of the left kan +extensions that define the representables (usually this is just wherever the +colimit leg sends 1, as there is often just one X part in the representable X). + +WARNING: this is only viable if the timer associated with the rewrite rule is +symmteric with respect to the disjoint representables. +""" +struct RepresentableP <: PatternType + parts::Dict{Symbol, Vector{Int}} +end +Base.keys(p::RepresentableP) = keys(p.parts) + +""" +Analyze a pattern to find the most efficient pattern type for it. +""" +function pattern_type(r::Rule) + p = pattern(r) + S = acset_schema(p) + + # Check empty case + isempty(p) && return EmptyP() + + # Determine if pattern is a coproduct of representables + repr_loc = DefaultDict{Symbol, Vector{Int}}(() -> Int[]) + reprs = repr_dict(typeof(p), S) + ccs, iso′ = connected_acset_components(p) + iso = invert_iso(iso′) + for cc_leg in legs(ccs) + found = false + for (o, (repr, i)) in pairs(reprs) + α = isomorphism(repr, dom(cc_leg)) + if !isnothing(α) + push!(repr_loc[o], iso[o](cc_leg[o](α[o](i)))) + found = true + break + end + end + found || break + end + length(ccs) == sum(length.(values(repr_loc))) && return RepresentableP(repr_loc) + + # Determine if pattern is discrete + # all(ob(S)) do o + # nparts(p, o) == 0 || isempty(homs(S, from=o)) + # end && return DiscreteP(Dict(o => nparts(p, o) for o in ob(S))) + + return RegularP() # no special case found +end + +""" +A stochastic rewrite rule with a dependent hazard rate +""" +struct ABMRule + rule::Rule + timer::AbsTimer + pattern_type::PatternType + ABMRule(r::Rule, t::AbsTimer) = new(r, t, pattern_type(r)) +end +getrule(r::ABMRule) = r.rule +pattern_type(r::ABMRule) = r.pattern_type +pattern(r::ABMRule) = pattern(getrule(r)) +right(r::ABMRule) = right(getrule(r)) + +abstract type AbsDynamics end + +"""Use a petri net with rates""" +struct PetriDynamics <: AbsDynamics + val::AbstractReactionNet +end + +"""Continuous dynamics""" +struct ABMFlow + pat::ACSet + dyn::AbsDynamics + mapping::Vector{Pair{Symbol, Int}} # pair pat's variables w/ dyn quantities +end + +""" +An agent-based model. +""" +struct ABM + rules::Vector{ABMRule} + dyn::Vector{ABMFlow} + ABM(rules, dyn=[]) = new(rules, dyn) +end +additions(abm::ABM) = right.(abm.rules) + +"""A collection of timers associated at runtime w/ an ABMRule""" +abstract type AbsHomSet end +struct EmptyHomSet <: AbsHomSet end +struct DiscreteHomSet <: AbsHomSet end +struct ExplicitHomSet <: AbsHomSet val::IncHomSet end +Base.keys(h::ExplicitHomSet) = keys(h.val) +Base.pairs(h::ExplicitHomSet) = pairs(h.val) +Base.getindex(h::ExplicitHomSet, i) = h.val[i] +deletion!(h::ExplicitHomSet, m) = deletion!(h.val, m) +addition!(h::ExplicitHomSet, k, r, u) = addition!(h.val, k, r, u) + +"""Initialize runtime hom-set given the rule and the initial state""" +function init_homset(rule::ABMRule, state::ACSet, additions::Vector{<:ACSetTransformation}) + p, sd = pattern_type(getrule(rule)), state_dep(rule.timer) + p == EmptyP() && return EmptyHomSet() + (sd || p == RegularP()) && return ExplicitHomSet(IncHomSet(pattern(rule), additions, state)) + @assert p isa RepresentableP "$(typeof(p))" + return DiscreteHomSet() +end + +const KeyType = Union{Pair{Int,Int}} # connected component homset + Tuple{Int,Vector{Pair{Int,Int}}} # multi-component homset +const default_sampler = FirstToFire{ + Union{Pair{Int, Nothing}, # non-explicit homset + Pair{Int, KeyType}}, # explicit homset + Float64} + +""" +Data structure for maintaining simulation information while running an ABM +""" +mutable struct RuntimeABM + state::ACSet + const clocks::Vector{AbsHomSet} + tnow::Float64 + nevent::Int + const sampler::SSA # stochastic simulation algorithm + const rng::Distributions.AbstractRNG + function RuntimeABM(abm::ABM, init::T; sampler=default_sampler) where T<:ACSet + rt = new(init, init_homset.(abm.rules, Ref(init), Ref(additions(abm))), + 0., 0, sampler(), Random.RandomDevice()) + for (i, homset) in enumerate(rt.clocks) + kv = homset isa ExplicitHomSet ? pairs(homset) : [nothing => create(init)] + for (key, val) in kv + haz = get_hazard(val, 0., abm.rules[i].timer) + enable!(rt.sampler, i => key, haz, 0., 0., rt.rng) + end + end + rt + end +end + +state(r::RuntimeABM) = r.state +Base.haskey(rt::RuntimeABM, k) = haskey(rt.sampler.transition_entry, k) + +"""Pick the next random event, advance the clock""" +function Fleck.next(rt::RuntimeABM) + rt.nevent += 1 + (rt.tnow, which) = next(rt.sampler, rt.tnow, rt.rng) + which +end + + +""" +Get match returns a randomly chosen morphism for the aggregate rule + +TODO incorporate the number of possibilities as a multiplier for the rate +""" +get_match(::EmptyP, L::ACSet, G::ACSet, ::EmptyHomSet, ::Nothing) = create(G) +function get_match(P::RepresentableP, L::T, G::ACSet, ::DiscreteHomSet, ::Nothing) where T<:ACSet + initial = Dict(map(collect(pairs(P.parts))) do (o, idxs) + o => Dict(idx => rand(parts(G, o)) for idx in idxs) + end) + homomorphism(L, G; initial) +end +get_match(::RegularP, ::ACSet, ::ACSet, hs::ExplicitHomSet, key::KeyType) = hs[key] + +""" +A trajectory of an ABM: each event time and result of `save`. +""" +struct Traj + init::ACSet + events::Vector{Tuple{Float64, Int, Any}} + hist::Vector{Span{<:ACSet}} +end +Traj(x::ACSet) = Traj(x, Pair{Float64, Any}[], Span{ACSet}[]) + +function Base.push!(t::Traj, τ::Float64,rule::Int, v::Any, sp::Span{<:ACSet}) + push!(t.events, (τ, rule, v)) + isempty(t.hist) || codom(left(sp)) == codom(right(last(t.hist))) || error( + "Bad history \n$(codom(left(sp))) \n!= \n$(right(last(t.hist)))" + ) + push!(t.hist, sp) +end + +Base.isempty(t::Traj) = isempty(t.events) + +Base.length(t::Traj) = length(t.events) + +const MAXEVENT = 10 + +""" +Run an ABM, creating a fresh trajectory. +""" +function run!(abm::ABM, init::T; save=deepcopy, maxevent=MAXEVENT, maxtime=Inf, + kw...) where T<:ACSet + run!(abm::ABM, RuntimeABM(abm, init; kw...); save, maxevent) +end + +function run!(abm::ABM, rt::RuntimeABM, output::Union{Traj,Nothing}=nothing; + save=deepcopy, maxevent=MAXEVENT, maxtime=Inf) + output = isnothing(output) ? Traj(rt.state) : output + + log!(rule::Int, sp::Span) = push!(output, rt.tnow, rule, save(rt.state), sp) + + disable!′(which) = disable!(rt.sampler, which, rt.tnow) + + function enable!′(m::ACSetTransformation, rule::Int, key=nothing) + haz = get_hazard(m, rt.tnow, abm.rules[rule].timer) + enable!(rt.sampler, rule => key, haz, rt.tnow, rt.tnow, rt.rng) + end + + while rt.nevent < maxevent && rt.tnow < maxtime + which = next(rt) # get next event, update time + isnothing(which) && return output # end b/c no more events! + + event, key = which + rule, clock = abm.rules[event], rt.clocks[event] + @debug "$(length(output)): event $event fired @ $(rt.tnow)" + + m = get_match(pattern_type(rule), pattern(rule), rt.state, clock, key) + update_maps = rewrite_match_maps(getrule(abm.rules[event]), m) + rh, kh, kg = update_maps[:rh], update_maps[:kh], update_maps[:kg] + rt.state = codom(rh) # update runtime state + log!(event, Span(kg, kh)) # record state after event + + if pattern_type(rule) == EmptyP() # "always enabled" need special treatment + disable!′(which) # their hom-set won't change, so clocks + enable!′(create(rt.state), event) # won't be reset, so do it manually + end + + # update matches for all events + for (t, (ruleₜ, clockₜ)) in enumerate(zip(abm.rules, rt.clocks)) + pt = pattern_type(ruleₜ) + if pt == RegularP() # update explicit hom-set + homs = clockₜ.val + for d in Incremental.deletion!(clockₜ, kg) + disable!′(t => d) # disable clocks which are invalidated + end + for a in Incremental.addition!(clockₜ, event, rh, kh) + enable!′(clockₜ[a], t, a) + end + elseif pt isa RepresentableP + if !all(o->all(is_isomorphic, [kg[o], kh[o]]), keys(pt)) + haskey(rt, t => nothing) && disable!′(t => nothing) + if all(>(0), nparts.(Ref(rt.state), collect(keys(pt)))) + enable!′(create(rt.state), t) + end + end + end + end + end + return output +end + + +end # module \ No newline at end of file diff --git a/src/AlgebraicABMS.jl b/src/AlgebraicABMS.jl deleted file mode 100644 index 5f66e5e..0000000 --- a/src/AlgebraicABMS.jl +++ /dev/null @@ -1,15 +0,0 @@ -""" Some description of ths package -""" -module AlgebraicABMS - -export hello - -using Catlab - -""" hello(name::String) - -Returns the string "Hello, !" where `` is replaced with the provided parameter -""" -hello(name::String) = string("Hello, ", name, "!") - -end diff --git a/src/AlgebraicABMs.jl b/src/AlgebraicABMs.jl new file mode 100644 index 0000000..176d263 --- /dev/null +++ b/src/AlgebraicABMs.jl @@ -0,0 +1,17 @@ +""" Some description of ths package +""" +module AlgebraicABMs + + +using Reexport + +include("Distributions.jl") +include("RewriteSemiMarkov.jl") +include("PetriInterface.jl") +include("ABMs.jl") + +@reexport using .Distributions +@reexport using .PetriInterface +@reexport using .RewriteSemiMarkov +@reexport using .ABMs +end diff --git a/src/Distributions.jl b/src/Distributions.jl new file mode 100644 index 0000000..8436294 --- /dev/null +++ b/src/Distributions.jl @@ -0,0 +1,31 @@ +module Distributions +using Distributions, SpecialFunctions +export weibullpar +""" +Get the shape and scale parameters for a Weibull distribution with a given mean and variance. +This code is a translation of the R "mixdist" package's +[function](https://github.com/cran/mixdist/blob/master/R/weibullpar.R) +""" +function weibullpar(mu, sigma) + cv = sigma / mu + if cv < 1e-06 + nu = cv / (sqrt(trigamma(1)) - cv * digamma(1)) + shape = 1 / nu + scale = mu / (1 + nu * digamma(1)) + else + aa = log(cv^2 + 1) + nu = 2 * cv / (1 + cv) + while true + gb = (lgamma(1 + 2 * nu) - 2 * lgamma(1 + nu) - aa) / (2 * (digamma(1 + 2 * nu) - digamma(1 + nu))) + nu -= gb + if abs(gb) < 1e-12 + break + end + end + shape = 1 / nu + scale = exp(log(mu) - lgamma(1 + nu)) + end + return shape, scale +end + +end # module diff --git a/src/PetriInterface.jl b/src/PetriInterface.jl new file mode 100644 index 0000000..d30bf49 --- /dev/null +++ b/src/PetriInterface.jl @@ -0,0 +1,114 @@ +""" +Petri nets as a presentation of certain kinds of rewrite rules. +""" +module PetriInterface +export PetriNetCSet + +using AlgebraicRewriting, Catlab, AlgebraicPetri +using Fleck +using Random +import AlgebraicRewriting.Incremental: IncSumHomSet + +using ..RewriteSemiMarkov: ClockSystem, ClockKeyType, sampl, rng +import ..RewriteSemiMarkov: run! + +"""Give a name for a Petri Net transition (name = label)""" +ob_name(pn::LabelledPetriNet, s::Int)::Symbol = pn[s, :sname] + +"""Give a name for a Petri Net transition (e.g. name = "S3" for species #3)""" +ob_name(::PetriNet, s::Int)::Symbol = Symbol("S$s") + +""" +Creates a discrete C-Set from a Petri net with one object for each species in +the Petri net. By default, this creates an *empty* C-Set instance, but there +are two ways one may also wish to specify how many tokens are in each species. +One can give a vector, where the indices correspond to the indices of the S +table of the petri net. Alternatively, one can give keyword arguments where +the keys are the names of the species (as determined by `ob_name`). + +For example, `PetriNetCSet(sir_labeled_pn, S=20, I=1)` would create an *instance* +of a C-Set that has three tables ("S","I","R"), no morphisms nor attributes, +and that instance would have 20 rows in the "S" table and 1 row in the "I" +table. In general, instances on this schema are effectively named tuples +(S::Int,I::Int,R::Int). +""" +function PetriNetCSet(pn::AbstractPetriNet, args=[]; kw...) + res = AnonACSet(BasicSchema(ob_name.(Ref(pn), parts(pn, :S)),[])) + for (arg, s) in zip(args, parts(pn, :S)) + add_parts!(res, ob_name(pn, s), arg) # Add tokens from Vector{Int} + end + for (s, arg) in pairs(kw) + add_parts!(res, s, arg) # Add tokens by name + end + res +end + +"""Assumes that tokens are deleted and recreated, rather than preserved""" +function make_rule(pn::Union{PetriNet, LabelledPetriNet}, t::Int) + L, R = LR = [PetriNetCSet(pn) for _ in 1:2] + add_part!.(Ref(L), ob_name.(Ref(pn), pn[incident(pn, t, :it), :is])) + add_part!.(Ref(R), ob_name.(Ref(pn), pn[incident(pn, t, :ot), :os])) + Rule(create.(LR)...) +end + +make_rules(pn) = make_rule.(Ref(pn), parts(pn, :T)) + +""" +Convert a Petri net into a `ClockSystem`, given a dictionary of timers +(corresponding to the transitions) as well as an initial state (such that the +incremental hom-sets can be initialized). The schema for `init` is expected to +be the `MarkedLabelledPetriNet` schema if `spn::T` is itself a +`MarkedLabelledPetriNet`, but if `spn` is an unmarked Petri Net then it is +assumed `init` will be the schema generated by `PetriNetCSet` above. + +The result can be used with `run!(::ClockSystem, ::ACSet)` where the second +parameter should be the same `init` ACSet used to generate the `ClockSystem`. +""" +function to_clocksys(spn::AbstractPetriNet, init::ACSet, clockdists) + clock = @acset ClockSystem begin + Global=1 + rng=[Random.RandomDevice()] + sampler=[FirstToFire{ClockKeyType, Float64}()] # Fleck sampler + Event=nt(spn) + name=spn[:,:tname] + end + + # add rules, distributions + for t in parts(clock,:Event) + clock[t,:rule] = make_rule(spn, t) + clock[t,:dist] = clockdists[clock[t,:name]] + end + + # make incremental homsets after all rules are made + for t in parts(clock, :Event) + clock[t,:match] = IncSumHomSet(IncHomSet(codom(left(clock[t,:rule])), + right.(clock[:,:rule]), init)) + end + + # which rules are always enabled? + aa = [t for t in parts(spn, :T) if isempty(incident(spn, t, :it))] + add_parts!(clock, :AlwaysEnabled, length(aa); always_enabled=aa) + + for t in parts(clock,:Event) + newkeys = [(t,k) for k in keys(clock[t,:match])] + add_parts!( + clock, :Clock, length(newkeys); + key = newkeys, event = fill(t, length(newkeys)) + ) + for c in newkeys + enable!(sampl(clock), c, clock[t, :dist](0.), 0., 0., rng(clock)) + end + end + + return clock +end + +""" +Convert a Petri net into a system of rewrite rules +""" +function run!(mlpn::T, clockdists, init::ACSet; kw...) where {T<:AbstractPetriNet} + run!(to_clocksys(mlpn, init, clockdists), init; kw...) +end + + +end # module \ No newline at end of file diff --git a/src/RewriteSemiMarkov.jl b/src/RewriteSemiMarkov.jl new file mode 100644 index 0000000..0373bf9 --- /dev/null +++ b/src/RewriteSemiMarkov.jl @@ -0,0 +1,118 @@ +module RewriteSemiMarkov + +export run! + +using Catlab, AlgebraicRewriting, AlgebraicPetri +using Random +using Fleck +# -------------------------------------------------------------------------------- +# we want something to store the rules, clocks associated to each, and their type + +""" + A presentation of a clock system, which is responsible for maintaining all the homsets for each +event possible in the system, and the set of enabled clocks for homset. It also stores the functions that +return a distribution over possible waiting times given a current simulation time. +""" +@present SchClockSystem(FreeSchema) begin + Clock :: Ob # a single ID + Event :: Ob # an entire class of clocks (e.g. "typed" clocks) + event::Hom(Clock,Event) + (NameType, RuleType, DistType, MatchType, KeyType, RngType, SampType)::AttrType + name::Attr(Event,NameType) # each event has a name + rule::Attr(Event,RuleType) # each event has a rewrite rule + dist::Attr(Event,DistType) # each event has a function that returns a firing time when it's enabled + match::Attr(Event,MatchType) # each event has a homset + key::Attr(Clock,KeyType) + + AlwaysEnabled :: Ob # a subset of Events + always_enabled::Hom(AlwaysEnabled, Event) + + Global :: Ob # a set with one element, for ACSet-level parameters + rng::Attr(Global, RngType) # random number generator + sampler::Attr(Global, SampType) +end + + +# const ClockKeyType = Tuple{Int,Int,Int} # for single connected component homsets +const ClockKeyType = Tuple{Int,Vector{Pair{Int,Int}}} + +""" +The acset type that stores instances of `SchClockSystem` +""" +@acset_type AbsClockSystem(SchClockSystem, index=[:event], unique_index=[:name,:key]) +const ClockSystem=AbsClockSystem{Symbol,Rule,Function,Incremental.IncSumHomSet,ClockKeyType, AbstractRNG, SSA} + +rng(c::ClockSystem) = only(c[:rng]) +sampl(c::ClockSystem) = only(c[:sampler]) + +""" +Given a `spn<:AbstractLabelledPetriNet` and a dictionary mapping each transition name to +a function that takes in parameter `t` and returns a distribution object, sample a single trajectory +of the stochastic dynamics on the petri net until `maxevent` events occur. Print stuff (or not) for debugging +with `verbose`. +""" +function run!(clocksys::ClockSystem, state::T; save=deepcopy, + maxevent=1000, verbose=false) where {T<:ACSet} + nevent, tnow = 0, 0.0 + sample, RNG = sampl(clocksys), rng(clocksys) + + output = [0. => save(state)] # record initial state + + # when and what will happen next? + (tnow, which) = next(sample, tnow, RNG) + nevent += 1 + + while nevent ≤ maxevent + !verbose || println("event $(first(which)) fired at $tnow, total number of events: $(length(output))") + event = first(which) + update_maps = rewrite_match_maps( + clocksys[event, :rule], + clocksys[event, :match][last(which)] + ) + state = codom(update_maps[:rh]) + + # record state after event + push!(output, tnow => save(state)) + + # "always enabled" transitions need special treatment when they fire + # because their hom-set won't change, the clocks won't be reset, so we do it manually + if !isempty(incident(clocksys, event, :always_enabled)) + disable!(sample, which, tnow) + enable!(sample, which, clocksys[event, :dist](tnow), tnow, tnow, RNG) + end + + # update matches for all events + for t in parts(clocksys, :Event) + # update the hom-set + del = Incremental.deletion!(clocksys[t,:match], update_maps[:kg]) + add = Incremental.addition!(clocksys[t,:match], event, update_maps[:rh], update_maps[:kh]) + del = map(del) do k + (t,k) + end + add = map(add) do k + (t,k) + end + # clocks that are disabled in the new marking (state): + # 1: disable in the sampler + # 2: remove from the clock system + for c in del + disable!(sample, c, tnow) + end + rem_parts!(clocksys, :Clock, sort(vcat(incident(clocksys, del, :key)...))) + # clocks that are newly enabled in the new marking (state): + # 1: enable in the sampler + # 2: add to the clock system + for c in add + enable!(sample, c, clocksys[t, :dist](tnow), tnow, tnow, RNG) + end + add_parts!(clocksys, :Clock, length(add), key = add, event = t) + + end + (tnow, which) = next(sample, tnow, RNG) + nevent += 1 + end + + return output +end + +end # module diff --git a/test/ABMS.jl b/test/ABMS.jl new file mode 100644 index 0000000..d31947d --- /dev/null +++ b/test/ABMS.jl @@ -0,0 +1,35 @@ +module TestABMs + +# ENV["JULIA_DEBUG"] = "AlgebraicABMs" + +using Test +using AlgebraicABMs + +using Catlab, AlgebraicRewriting + +using AlgebraicABMs.ABMs: ABMRule, DiscreteHazard, ContinuousHazard, RegularP, + EmptyP, RepresentableP, RuntimeABM + +create_vertex = ABMRule(Rule(id(Graph()), create(ob(terminal(Graph)))), DiscreteHazard(1.)) +@test create_vertex.pattern_type == EmptyP() + +add_loop = ABMRule(Rule(id(Graph(1)), delete(Graph(1))), DiscreteHazard(1.5)) +@test add_loop.pattern_type isa RepresentableP + +rem_loop = ABMRule(Rule(delete(Graph(1)), id(Graph(1))), DiscreteHazard(2)) +@test rem_loop.pattern_type == RegularP() + +rem_edge = ABMRule(Rule(homomorphism(Graph(2), path_graph(Graph, 2); + monic=true), id(Graph(2))), ContinuousHazard(1)) +@test rem_edge.pattern_type isa RepresentableP + + +G = @acset Graph begin V=3; E=3; src=[1,1,1]; tgt=[1,1,2] end +abm = ABM([create_vertex, add_loop, rem_loop, rem_edge]) + +@test only(RuntimeABM(abm, G).clocks[3].val.match_vect)[1][:E](1) == 1 # One cached hom + + +traj = run!(abm, G); + +end # module diff --git a/test/PetriInterface.jl b/test/PetriInterface.jl new file mode 100644 index 0000000..432d309 --- /dev/null +++ b/test/PetriInterface.jl @@ -0,0 +1,16 @@ +module TestPetriInterface + +using Test +using AlgebraicABMs, Catlab, AlgebraicPetri + +sir_pn= @acset LabelledPetriNet begin + S=3; sname=[:S,:I,:R] + T=7; tname=[:inf,:rec,:birth,:deathS,:deathI,:deathR,:wane] + I=7; it=[1,1,2,4,5,6,7]; is=[1,2,2,1,2,3,3] + O=5; ot=[1,1,2,3,7]; os=[2,2,3,1,1] +end +t1, t2, t3, t4, t5, t6, t7 = AlgebraicABMs.PetriInterface.make_rules(sir_pn) +@test nparts.(Ref(codom(left(t1))), [:S,:I,:R]) == [1,1,0] # S, I + + +end # module diff --git a/test/Project.toml b/test/Project.toml index 2975b4c..aa9f405 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,12 @@ [deps] -AlgebraicABMS = "1301f41b-2c8a-4430-9d85-ab97d23a62cd" +AlgebraicABMs = "5a5e3447-9604-46e6-8d91-cb86f5f51721" +AlgebraicPetri = "4f99eebe-17bf-4e98-b6a1-2c4f205a959b" +AlgebraicRewriting = "725a01d3-f174-5bbd-84e1-b9417bad95d9" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Fleck = "5bb9b785-358c-4fee-af0f-b94a146244a8" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/aqua.jl b/test/aqua.jl index e290f4a..31ab5a9 100644 --- a/test/aqua.jl +++ b/test/aqua.jl @@ -1,3 +1,3 @@ -using Aqua, AlgebraicABMS +using Aqua, AlgebraicABMs -Aqua.test_all(AlgebraicABMS, ambiguities=false,) +Aqua.test_all(AlgebraicABMs, ambiguities=false,) diff --git a/test/core.jl b/test/core.jl deleted file mode 100644 index d9bf789..0000000 --- a/test/core.jl +++ /dev/null @@ -1 +0,0 @@ -@test hello("World") == "Hello, World!" diff --git a/test/runtests.jl b/test/runtests.jl index 16e23d8..82102fb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,7 @@ using Test @testset "Code Quality (Aqua.jl)" begin - include("aqua.jl") + # include("aqua.jl") end @testset "Core" begin