Merge pull request #37 from AlgebraicJulia/subtimer

"Basis" for rewrite rules

docs/literate/showcase.jl

@@ -0,0 +1,47 @@
+# # Showcasing AlgebraicABMs features
+# 
+# This file is an incomplete showcase of some less obvious features in 
+# AlgebraicABMs. It is not an intro-tutorial (the other demos are more geared 
+# towards that). To begin, we want to load our packages with `using
+
+using AlgebraicABMs, Catlab, AlgebraicRewriting
+ENV["JULIA_DEBUG"] = "AlgebraicABMs"; # turn off @debug messages for this package
+
+# ## Basis 
+
+# Rewrote rules have associated timers which are scheduled to fire in the future.
+# For normal rules, at any point in time, there is exactly one such timer per 
+# match morphism from the pattern $L$ of the rule into the present state of the 
+# world, $X$. This is often what we want, though in some circumstances, we want the 
+# 'per-what?' of the rule to be something different from the pattern. Thus we 
+# have the ability to explicitly give a *basis* for the rule, given as a 
+# map $B \rightarrowtail L$ into the rule's pattern. The rule will have timers 
+# per every 'match' $B \rightarrow X$ (which is just a *partial* match 
+# $L \nrightarrow X$). Once it's time to fire, the rest of the match is 
+# completed at random (if possible, otherwise nothing happens).
+
+# For example, we may want to associate the timers of an infection process with 
+# simply the infected people, even though the pattern of an infection rule is 
+# a suceptible + an infected person. E.g. "every infected person tries to infect
+# *someone* once per day", which would not be affected by how many susceptible 
+# people there are. Let's demonstrate this:
+
+@present SchSI(FreeSchema) begin
+  (S,I)::Ob
+end
+
+@acset_type SI(SchSI)
+
+si   = @acset SI begin S=1; I=1 end 
+i    = @acset SI begin I=1 end 
+ii   = @acset SI begin I=2 end 
+init = @acset SI begin S=10000;I=1 end 
+inf_rule = Rule(homomorphism(i, si), homomorphism(i, ii; any=true))
+basis_hom = homomorphism(i, si) # our actual basis is just a single infected
+
+abm_rule = ABMRule(inf_rule, DiscreteHazard(1); basis = basis_hom)
+abm = ABM([abm_rule])
+run!(abm, init; maxevent=3);
+
+# We can see that after the first timestep there is 1 infection, then 2 on the 
+# next timestep, then four on the next one.

docs/make.jl

@@ -54,6 +54,7 @@ makedocs(
   pages=Any[
     "AlgebraicABMs.jl"=>"index.md",
     "Examples"=>Any[
+     "generated/showcase.md",
       "generated/sir_petri.md",
       "generated/game_of_life.md",
       "generated/lotka_volterra.md",

src/ABMs.jl

@@ -12,11 +12,12 @@ using StructEquality
 using Catlab, AlgebraicRewriting
 using AlgebraicRewriting.Incremental.Algorithms: connected_acset_components, pull_back
 using AlgebraicRewriting.Rewrite.Migration: repr_dict
+using Catlab.CategoricalAlgebra.Chase: extend_morphism_constraints
 using AlgebraicRewriting.Rewrite.Utils: get_pmap, get_rmap, get_expr_binding_map
-import Catlab: right
+import Catlab: left, right
 import AlgebraicRewriting: get_match, ruletype, addition!, deletion!, get_matches
 
-import ..Upstream: pattern, pops!
+import ..Upstream: pattern, pops!, IncHomSet_basis
 
 # Timers
 ########
@@ -199,14 +200,19 @@ const Maybe{T} = Union{Nothing, T}
 
 """
 A stochastic rewrite rule with a dependent hazard rate
+
+A basis is a subobject of the pattern of the rule for which we want a timer 
+per match. By default, the basis ↣ pattern map is just id(pattern).
+
 """
 @struct_hash_equal struct ABMRule
   rule::Rule
-  timer::AbsTimer 
+  timer::AbsTimer
+  basis::Maybe{ACSetTransformation}
   name::Maybe{Symbol}
   pattern_type::PatternType
-  ABMRule(r::Rule, t::AbsTimer; name=nothing) = 
-    new(r, t, name, pattern_type(r, is_exp(t)))
+  ABMRule(r::Rule, t::AbsTimer; basis=nothing, name=nothing) = 
+    new(r, t, basis, name, pattern_type(r, is_exp(t)))
 end
 
 # Give name as first arg rather than as kwarg
@@ -221,15 +227,20 @@ pattern_type(r::ABMRule) = r.pattern_type
 
 pattern(r::ABMRule) = pattern(getrule(r))
 
+left(r::ABMRule) = left(getrule(r))
 right(r::ABMRule) = right(getrule(r))
 
 ruletype(r::ABMRule) = ruletype(getrule(r))
 
+basis(r::ABMRule) = r.basis
+
+basis_pattern(r::ABMRule) = isnothing(r.basis) ? codom(left(r)) : dom(basis(r))
+
 get_matches(r::ABMRule, args...; kw...) = 
   get_matches(getrule(r), args...; kw...)
 
 (F::Migrate)(r::ABMRule) = 
-  ABMRule(F(r.rule), r.timer; name=r.name)
+  ABMRule(F(r.rule), r.timer; basis=F(r.basis), name=r.name)
 
 """
 A type which implements AbsDynamics must be able to compiled to an ODE for some 
@@ -252,8 +263,8 @@ end
 end 
 
 # Accessing an IncHomSet
-const KeyType = Union{Pair{Int, Int}}       # connected comp. homset
-                      Vector{Pair{Int,Int}} # multi-component homset
+const KeyType = Union{Pair{Int, Int},        # connected comp. homset
+                      Vector{Pair{Int,Int}}} # multi-component homset
 
 """
 An agent-based model.
@@ -306,7 +317,8 @@ function init_homset(rule::ABMRule, state::ACSet,
   p, sd = pattern_type(rule), state_dep(rule.timer)
   p == EmptyP() && return EmptyHomSet()
   (sd || p == RegularP()  
-   ) && return ExplicitHomSet(IncHomSet(getrule(rule), state,  additions))
+   ) && return ExplicitHomSet(IncHomSet_basis(getrule(rule), state,  additions; 
+                                        basis=basis_pattern(rule)))
   @assert p isa RepresentableP  "$(typeof(p))"
   return RepresentableHomSet()
 end 
@@ -318,7 +330,7 @@ const default_sampler = FirstToFire{
   Float64}
 
 """
-Data @struct_hash_equal structure for maintaining simulation information while running an ABM
+Data structure for maintaining simulation information while running an ABM
 """
 mutable struct RuntimeABM
   state::ACSet
@@ -395,6 +407,14 @@ function pops!(rt::RuntimeABM)::Vector{Pair{Int, Maybe{KeyType}}}
 end
 
 
+function get_match(pat::PatternType, L::ACSet, G::ACSet, timer::AbsHomSet, key; 
+                   basis::Maybe{ACSetTransformation}) 
+  isnothing(basis) && return get_match(pat, L, G, timer, key)
+  m = get_match(pat, dom(basis), G, timer, key)
+  initial = extend_morphism_constraints(m, basis)
+  rand(homomorphisms(L, G; initial))
+end
+
 """
 Get match returns a randomly chosen morphism for the aggregate rule
 """
@@ -458,7 +478,7 @@ function run!(abm::ABM, rt::RuntimeABM, output::Traj;
     push!(output, (rt.tnow, rule, getname(rule), save(rt.state), sp))
   disable!′(key::Pair) = disable!(rt.sampler, key, rt.tnow)
   disable!′(i::Int) = disable!′(i => nothing)
-  function enable!′(m::ACSetTransformation, rule_id::Int, key=nothing) 
+  function enable!′(m::ACSetTransformation, rule_id::Int, key::Maybe{KeyType}=nothing) 
     rule = abm.rules[rule_id]
     haz = get_hazard(pattern_type(rule), m, rt.tnow, rule.timer)
     enable!(rt.sampler, rule_id => key, haz, rt.tnow, rt.tnow, rt.rng)
@@ -494,7 +514,8 @@ function run!(abm::ABM, rt::RuntimeABM, output::Traj;
         rule::ABMRule, clocks::AbsHomSet = abm.rules[event], rt.clocks[event]
         rule′::Rule, rule_type::Symbol = getrule(rule), ruletype(rule)
         # If RegularPattern, we have an explicit match, otherwise randomly pick one
-        m = get_match(pattern_type(rule), pattern(rule), rt.state, clocks, key)
+        m = get_match(pattern_type(rule), pattern(rule), rt.state, clocks, key; 
+                      basis=basis(rule))
         # bring the match 'up to speed' given the previous (simultanous) updates
         for (l, r) in first.(update_data)
           m = pull_back(l, m) ⋅ r

src/Upstream.jl

@@ -4,6 +4,8 @@ module Upstream
 using Catlab, AlgebraicRewriting
 import Catlab: is_isomorphic, Presentation
 using AlgebraicRewriting.Rewrite.Migration: pres_hash
+import AlgebraicRewriting: IncHomSet
+using AlgebraicRewriting.Incremental.IncrementalConstraints: AC, PAC, NAC
 using CompetingClocks: FirstToFire, disable!, next
 using Distributions: AbstractRNG
 
@@ -28,6 +30,20 @@ end
 
 # Upstream to AlgRewriting
 ##########################
+"""Optionally use a different pattern than the L of the rule"""
+function IncHomSet_basis(rule::Rule{T}, state::ACSet, additions=ACSetTransformation[]; 
+                         basis=nothing) where T
+  pac, nac = [], []
+  dpo = (T == :DPO) ? [left(rule)] : ACSetTransformation[]
+  right(rule) ∈ additions || push!(additions, right(rule))
+  for c in AC.(rule.conditions, Ref(additions), Ref(dpo))
+    c isa PAC && push!(pac, c)
+    c isa NAC && push!(nac, c)
+  end
+  pat = isnothing(basis) ? codom(left(rule)) : basis
+  IncHomSet(pat, additions, state; monic=rule.monic, pac, nac)
+end
+
 
 # CompetingClocks
 #################

src/Visualization.jl

@@ -1,7 +1,9 @@
 module Visualization
 
+export graphviz_write
+
 import AlgebraicRewriting: view_traj 
-using Catlab: codom, right
+using Catlab: codom, right, to_graphviz
 
 using ..ABMs
 using ..ABMs: Traj
@@ -29,4 +31,12 @@ function Base.view(t::Traj, viewer; dirname="default")
   end
 end 
 
+function graphviz_write(x, dirname="default")
+  G = to_graphviz(x)
+  open(dirname, "w") do io
+    show(io, "image/svg+xml", G)
+  end
+  G
+end
+
 end # module

test/ABMs.jl

@@ -62,8 +62,28 @@ create_vert = ABMRule(Rule(id(Graph()),  create(Graph(1))), DiscreteHazard(1.));
 abm = ABM([create_loop, create_vert]);
 traj = run!(abm, Graph(); maxtime=5);
 
-# ODEs
-######
+
+# Basis
+#######
+using AlgebraicABMs, Catlab, AlgebraicRewriting
+
+# Rule which fires once per vertex (it tries to add an edge to an 
+# arbitrary other vertex) - rate is based on vertex pattern, not pattern for 
+# the rule itself, which has two 
+add_edge = ABMRule(Rule(id(Graph(2)), 
+                        homomorphism(Graph(2), path_graph(Graph, 2); 
+                                     any=true, monic=true)),
+                   DiscreteHazard(1.),
+                   basis=homomorphism(Graph(1), Graph(2); any=true))
+abm = ABM([add_edge])
+rt = RuntimeABM(abm, Graph(3))
+traj = run!(abm, Graph(3); maxtime=3);
+
+view(traj, graphviz_write)
+
+
+# ODEs (IN PROGRESS)
+#####################
 using AlgebraicABMs, AlgebraicRewriting, Catlab
 
 # State of world: a set of free-floating Float64s