GoL name changes, auto-resample matches that just fired if valid

docs/literate/game_of_life.jl

@@ -14,7 +14,7 @@ using AlgebraicABMs, Catlab, AlgebraicRewriting
   live::Hom(Life,V)
 end
 
-@acset_type Life(SchLifeGraph) <: AbstractSymmetricGraph;
+@acset_type LifeState(SchLifeGraph) <: AbstractSymmetricGraph;
 
 to_graphviz(SchLifeGraph)
 
@@ -27,15 +27,15 @@ to_graphviz(SchLifeGraph)
   coords::Attr(V, Coords)
 end
 
-@acset_type LifeCoords(SchLifeCoords){Tuple{Int,Int}} <: AbstractSymmetricGraph;
+@acset_type LifeStateCoords(SchLifeCoords){Tuple{Int,Int}} <: AbstractSymmetricGraph;
 
 to_graphviz(SchLifeCoords)
 
 # Create a regular grid with periodic boundary conditions with an initial state
 function make_grid(curr::AbstractMatrix)
   n, m = size(curr)
   n == m || error("Must be square")
-  X, coords = LifeCoords(), Dict()
+  X, coords = LifeStateCoords(), Dict()
   for i in 1:n, j in 1:n
     coords[i=>j] = add_vertex!(X; coords=(i, j))
     Bool(curr[i, j]) && add_part!(X, :Life, live=coords[i=>j])
@@ -54,14 +54,14 @@ make_grid(n::Int, random=true) =
   make_grid((random ? rand : zeros)(Bool, (n, n)));
 
 # Visualize a game state (using coordinates if available)
-function view_life(X::Union{Life, LifeCoords}, pth=tempname())
+function view_life(X::Union{LifeState, LifeStateCoords}, pth=tempname())
   pg = PropertyGraph{Any}(; prog="neato", graph=Dict(),
     node=Dict(:shape => "circle", :style => "filled", :margin => "0"),
     edge=Dict(:dir => "none", :minlen => "1"))
   add_vertices!(pg, nparts(X, :V))
   for v in vertices(X)
     set_vprop!(pg, v, :fillcolor, isempty(incident(X, v, :live)) ? "red" : "green")
-    if X isa LifeCoords 
+    if X isa LifeStateCoords 
       x, y = X[v, :coords]
       set_vprop!(pg, v, :pos, "$x,$(y)!")
     end
@@ -87,23 +87,23 @@ is to assign "variables" for the values of the coordinates).
 
 idₒ = Dict(x => x for x in Symbol.(generators(SchLifeGraph, :Ob)))
 idₘ = Dict(x => x for x in Symbol.(generators(SchLifeGraph, :Hom)))
-AddCoords = Migrate′(idₒ, idₘ, SchLifeGraph, Life, SchLifeCoords, LifeCoords; delta=false);
+AddCoords = Migrate′(idₒ, idₘ, SchLifeGraph, LifeState, SchLifeCoords, LifeStateCoords; delta=false);
 RemCoords = DeltaMigration(FinFunctor(idₒ, idₘ, SchLifeGraph, SchLifeCoords))
 # ## Helper constants and functions 
-const Dead = Life(1) # a single dead cell
-const Live = @acset Life begin V=1; Life=1; live=1 end # a single living cell
-const to_life = homomorphism(Dead, Live)  # the unique map Dead → Live
+const DeadCell = LifeState(1) # a single dead cell
+const LiveCell = @acset LifeState begin V=1; Life=1; live=1 end # a single living cell
+const to_life = homomorphism(DeadCell, LiveCell)  # the unique map Dead → Live
 
 """Create a context of n living neighbors for either a dead or alive cell"""
 function living_neighbors(n::Int; alive=true)::ACSetTransformation
-  X = Life(1)
+  X = LifeState(1)
   alive && add_part!(X, :Life, live=1)
   for _ in 1:n
     v = add_part!(X, :V)
     add_part!(X, :Life, live=v)
     add_edge!(X, v, 1)
   end
-  homomorphism(alive ? Live : Dead, X; initial=(V=[1],))
+  homomorphism(alive ? LiveCell : DeadCell, X; initial=(V=[1],))
 end;
 
 PAC(m) = AppCond(m; monic=true) # Positive Application condition
@@ -114,13 +114,13 @@ TickRule(args...; kw...) = # Rule which fires on 1.0, 2.0, ...
 # ## Create model by defining update rules
 
 # A cell is born iff it has three living neighbors
-birth = TickRule(id(Dead), to_life; 
+birth = TickRule(id(DeadCell), to_life; 
                  ac=[PAC(living_neighbors(3; alive=false)),
                      NAC(living_neighbors(4; alive=false)),
                      NAC(to_life)]);
 
 # A cell is born iff it has ≥ 2 living neighbors but < 4 living neighbors
-death = TickRule(to_life, id(Dead); 
+death = TickRule(to_life, id(DeadCell); 
                  ac=[PAC(living_neighbors(2)), 
                      NAC(living_neighbors(4))]);
 
@@ -136,6 +136,4 @@ G = make_grid([1 0 1 0 1;
 view_life(G);
 G = make_grid(ones(1,1))
 # Run the model
-migrate(Life, G, RemCoords)
-get_matches(birth.rule, migrate(Life, G, RemCoords))
 res = run!(AddCoords(GoL), G; maxevent=2);

src/ABMs.jl

@@ -142,7 +142,8 @@ schema as an argument in order to compute representables that would otherwise
 be infinite.
 
 Even if the pattern is a coproduct of representables, we cannot use the 
-efficient encoding unless the distribution is either an exponential.
+efficient encoding unless the distribution is either an exponential 
+(or a single dirac delta - not yet supported).
 """
 function pattern_type(r::Rule, is_exp::Bool; schema=nothing)
   p = pattern(r)
@@ -227,6 +228,11 @@ end
   mapping::Vector{Pair{Symbol, Int}} # pair pat's variables w/ dyn quantities
 end 
 
+const Maybe{T} = Union{Nothing, T}
+
+const KeyType = Union{Pair{Int, Int}}                  # connected comp. homset
+                      Tuple{Int,Vector{Pair{Int,Int}}} # multi-component homset
+
 """
 An agent-based model.
 """
@@ -251,6 +257,8 @@ abstract type AbsHomSet end
 
 Base.keys(h::ExplicitHomSet) = keys(h.val)
 
+Base.haskey(h::ExplicitHomSet, k::KeyType) = haskey(h.val, k)
+
 Base.pairs(h::ExplicitHomSet) = pairs(h.val)
 
 Base.getindex(h::ExplicitHomSet, i) = h.val[i]
@@ -268,11 +276,6 @@ function init_homset(rule::ABMRule, state::ACSet, additions::Vector{<:ACSetTrans
   return DiscreteHomSet()
 end 
 
-const Maybe{T} = Union{Nothing, T}
-
-const KeyType = Union{Pair{Int, Int}}                  # connected comp. 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
@@ -309,10 +312,10 @@ Base.haskey(rt::RuntimeABM, k::Pair) = haskey(rt.sampler.transition_entry, k)
 Base.haskey(rt::RuntimeABM, k::Int) = 
   haskey(rt.sampler.transition_entry, k => nothing)
 
-"""Pick the next random event, advance the clock"""
+"""Pop 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)
+  (rt.tnow, which) = pop!(rt.sampler, rt.rng, rt.tnow)
   return which
 end
 
@@ -364,10 +367,10 @@ Run an ABM, creating a fresh runtime + trajectory.
 """
 function run!(abm::ABM, init::T; save=deepcopy, maxevent=MAXEVENT, maxtime=Inf, 
               kw...) where T<:ACSet 
-  run!(abm::ABM, RuntimeABM(abm, init; kw...), Traj(init); save, maxevent)
+  run!(abm::ABM, RuntimeABM(abm, init; kw...), Traj(init); save, maxtime, maxevent)
 end
 
-function run!(abm::ABM, rt::RuntimeABM, output::Traj; 
+function run!(abm::ABM, rt::RuntimeABM, output::Traj;
               save=deepcopy, maxevent=MAXEVENT, maxtime=Inf)
 
   # Helper functions that automatically incorporate the runtime `rt`
@@ -381,55 +384,55 @@ function run!(abm::ABM, rt::RuntimeABM, output::Traj;
 
   # Main loop
   while rt.nevent < maxevent && rt.tnow < maxtime
-    # get next event + update clock time
-    which = next(rt)
-
-    if isnothing(which)
+    if length(rt.sampler) == 0 
       @info "Stochastic scheduling algorithm ran out of events"
       return output
     end
 
-    # Unpack data associated with the current event
-    event::Int, key::Maybe{KeyType} = which
+    # Get next event + unpack data
+    event::Int, key::Maybe{KeyType} = next(rt) # updates the clock time
     rule::ABMRule, clocks::AbsHomSet = abm.rules[event], rt.clocks[event]
     rule_type::Symbol = ruletype(rule) # DPO, SPO, etc.
 
     @debug "$(length(output)): event $event fired @ $(rt.tnow)"
-
+    show(stdout,"text/plain", rt.state)
     # If RegularPattern, we have an explicit match, otherwise randomly pick one
     m = get_match(pattern_type(rule), pattern(rule), rt.state, clocks, key)
 
     # Excute rewrite rule and unpack results
     rw_result = (rule_type, rewrite_match_maps(getrule(rule), m))
     rmap_ = get_rmap(rw_result...)
-    xmap_ = get_expr_binding_map(getrule(rule), m, rw_result[2])
+    xmap = get_expr_binding_map(getrule(rule), m, rw_result[2])
     (lft, rght_) = get_pmap(rw_result...)
-    rmap, rght = compose.([rmap_,rght_],Ref(xmap_))
+    rmap, rght = compose.([rmap_,rght_], Ref(xmap))
     pmap = Span(lft, rght)
     rt.state = codom(rmap) # update runtime state
     log!(event, pmap)      # record event result
 
-    # update matches for all events 
-    #------------------------------
-    # The only time EmptyPattern rules update is when they are fired
-    if pattern_type(rule) == EmptyP()   
-      disable!′(which)
-      enable!′(create(rt.state), event)
-    end
     # All other rules can potentially update in response to the current event
     for (i, (ruleᵢ, clocksᵢ)) in enumerate(zip(abm.rules, rt.clocks))
       pt = pattern_type(ruleᵢ)
-      if pt == RegularP() # update explicit hom-set w/r/t span Xₙ ↩ • -> Xₙ₊₁
+      if pt == EmptyP()
+        i == event && enable!′(create(rt.state), i)
+      elseif pt == RegularP() # update explicit hom-set w/r/t span Xₙ ↩ • -> Xₙ₊₁
         for d in deletion!(clocksᵢ, lft)
-          disable!′(i => d) # disable clocks which are invalidated
+          # disable clocks which are invalidated
+          # note that (event,key) is already diabled
+          (i,d)==(event,key) || disable!′(i => d) 
         end
         for a in addition!(clocksᵢ, event, rmap, rght) # rght: R → Xₙ₊₁
           enable!′(clocksᵢ[a], i, a)
         end
+        # If the match that just fired is still preserved
+        if i == event && haskey(clocksᵢ, key)
+          enable!′(clocksᵢ[key], i, key)
+        end
       elseif pt isa RepresentableP
         relevant_obs = keys(pt)
         # we need to update current timer if # of parts has changed
-        if !all(ob -> allequal(nparts.(codom.(pmap), Ref(ob))), relevant_obs)
+        if i == event && all(>(0), nparts.(Ref(rt.state), relevant_obs))
+          enable!′(create(rt.state), i)
+        elseif !all(ob -> allequal(nparts.(codom.(pmap), Ref(ob))), relevant_obs)
           currently_enabled = haskey(rt, i)
           currently_enabled && disable!′(i) # Disable if active
           # enable new timer if possible to apply rule

src/Upstream.jl

@@ -6,6 +6,8 @@ using Catlab, AlgebraicRewriting
 import Catlab: acset_schema, is_isomorphic, Presentation
 using AlgebraicRewriting.Rewrite.Migration: pres_hash
 using AlgebraicRewriting.Incremental: key_dict
+using Fleck: SSA, disable!, next
+using Distributions: AbstractRNG
 
 # Upstream to Catlab
 ####################
@@ -68,5 +70,14 @@ struct Migrate′
     new(Migrate(o, h, t1, t2; delta), s1, isnothing(s2) ? s1 : s2)
 end
 
+# Fleck
+#######
+"""Get the next event and disable it"""
+function Base.pop!(sampler::SSA{K,T}, rng::AbstractRNG, 
+                   tnow::T)::Tuple{T,K} where {K,T}
+  (new_time, which) = next(sampler, tnow, rng)
+  disable!(sampler, which, new_time)
+  return (new_time, which)
+end 
 
 end # module

test/ABMS.jl

@@ -45,4 +45,12 @@ abm = ABM([create_loop, add_loop, rem_loop, rem_edge])
 traj = run!(abm, G; maxevent=10);
 @test length(traj) == 10
 
+traj = run!(ABM([rem_edge]), G);
+@test length(traj) == 3
+
+
+traj = run!(ABM([add_loop]), G);
+@test length(traj) > 3 # after we add a loop, the match persists and is resampled
+
+
 end # module