Merge pull request #31 from AlgebraicJulia/random_gol

Stochastic game of life example + adjusting to new hom search interface

Project.toml

@@ -26,8 +26,8 @@ PetriInterface = "AlgebraicPetri"
 
 [compat]
 AlgebraicPetri = "0.9"
-AlgebraicRewriting = "^0.3.5"
-Catlab = "^0.16.11"
+AlgebraicRewriting = "^0.4.0"
+Catlab = "^0.16.16"
 CompetingClocks = "0.1"
 DataStructures = "0.18"
 DifferentialEquations = "7"

benchmark/ABMs.jl

@@ -0,0 +1 @@
+# TODO

benchmark/README.md

@@ -0,0 +1,18 @@
+# AlgebraicABMs.jl benchmarks
+
+This directory contains benchmarks for different parts of AlgebraicABMs. To run all the
+benchmarks, launch `julia --project=benchmark` and enter:
+
+``` julia
+using PkgBenchmark
+import AlgebraicABMs
+
+benchmarkpkg(AlgebraicABMs)
+```
+
+To run a specific set of benchmarks, use the `script` keyword argument, for
+example:
+
+``` julia
+benchmarkpkg(AlgebraicABMs; script="benchmark/ABMs.jl")
+```

benchmark/benchmarks.jl

@@ -0,0 +1,7 @@
+using BenchmarkTools
+
+const SUITE = BenchmarkGroup()
+
+module BenchmarkABMs
+  include("ABMs.jl")
+end

docs/literate/game_of_life.jl

@@ -4,8 +4,9 @@
 # The first step of running a Julia program is to load the external libraries 
 # one will be using. We do this with a `using` statement.
 
-using AlgebraicABMs, Catlab, AlgebraicRewriting
+using AlgebraicABMs, Catlab, AlgebraicRewriting, Random, Test
 ENV["JULIA_DEBUG"] = "AlgebraicABMs"; # hide
+Random.seed!(100)
 
 #=
 ## Schema 
@@ -415,8 +416,8 @@ match_coords(birth, G)
 # This is easy! Pass in our model and our initial state. We optionally could 
 # limit the run via some maximum number of steps or time, but this one will 
 # achieve steady state within 3 time steps.
-
 res = run!(GoL_coords, G);
+@test length(res) == 13
 
 # ## View results
 
@@ -478,3 +479,35 @@ imgs[13]
 # 
 
 imgs[14]
+
+# ## Bonus: stochastic game of life
+
+#=
+Rather than having all cells update in lockstep, we could change the probability
+of firing from a Dirac delta distribution at t=1 to an exponential distribution,
+where the expected value is firing at t=1. This means cells will update one at 
+a time, as it is almost impossible for two events to occur at the same time.
+The change involved for this is simply replacing the `DiscreteHazard` of 
+`TickRule` with a `ContinuousHazard`.
+=#
+
+continuous_abm = ABM([ABMRule(r.name, r.rule, ContinuousHazard(1)) for r in GoL_coords.rules]) 
+
+res = run!(continuous_abm, G; maxevent=5)
+imgs = view(res, view_life);
+
+# Here is our starting point.
+
+imgs[1]
+
+# Let's look at the next few steps.
+
+imgs[2]
+
+# Then
+
+imgs[3]
+
+# Then
+
+imgs[4]

docs/literate/lotka_volterra.jl

@@ -67,13 +67,13 @@ well as be oriented in a particular direction.
 end;
 
 #=
-The vertices (patches of grass) have a "countdown" attribute,
-which measures how long it takes until the grass is ready to eat at a location. 
+A set of time counters are associated with vertices via "countdown"
+which tracks how long it takes until the grass is ready to eat at a location. 
 =#
 
 @present SchWSG <: SchWS begin
-  Countdown::AttrType
-  countdown::Attr(V, Countdown)
+  Time::Ob
+  countdown::Hom(Time, V)
 end;
 
 #=
@@ -87,7 +87,7 @@ eating food.
   sheep_eng::Attr(Sheep, Eng)
 end;
 
-@acset_type LV(SchLV){Int, Int} <: AbstractSymmetricGraph # Create LV datatype
+@acset_type LV(SchLV){Int} <: AbstractSymmetricGraph # Create LV datatype
 
 to_graphviz(SchLV; prog="dot")
 
@@ -108,7 +108,7 @@ model with the LV_viz schema.
   dirname::Attr(Direction, Name)
 end
 
-@acset_type LV_Viz(SchLV_Viz){Int, Int, Tuple{Int,Int}, String} <: AbstractSymmetricGraph
+@acset_type LV_Viz(SchLV_Viz){Int, Tuple{Int,Int}, String} <: AbstractSymmetricGraph
 
 const LV′ = Union{LV, LV_Viz};
 
@@ -124,7 +124,8 @@ function create_grid(n::Int)::LV_Viz
   N, W, S, E = add_parts!(lv, :Direction, 4; left=[2,3,4,1], right=[4,1,2,3], # hide
                           dirname=["N","W","S","E"]) # hide
   for i in 0:n-1, j in 0:n-1 # hide
-    coords[i=>j] = add_vertex!(lv; countdown=max(0, rand(-30:30)), coord=(i,j)) # hide
+    coords[i=>j] = add_vertex!(lv; coord=(i,j)) # hide
+    add_parts!(lv, :Time, max(0, rand(-30:30)); countdown=coords[i=>j]) # hide
   end # hide
   for i in 0:n-1, j in 0:n-1 # hide
     _, e = add_edge!(lv, coords[i=>j], coords[mod(i + 1, n)=>j]; dir=E) # hide
@@ -166,7 +167,7 @@ function view_LV(p::LV′, pth=tempname(); name="G", title="")
   pstr = ["$(i),$(j)!" for (i, j) in p[:coord]] # hide
   stmts = Statement[] # hide
   for s in 1:nv(p) # hide
-    gv = p[s, :countdown] # hide
+    gv = length(incident(p, s, :countdown)) # hide
     col = gv == 0 ? "lightgreen" : "tan" # hide
     push!(stmts, Node("v$s", Attributes( # hide
       :label => gv == 0 ? "" : string(gv), :shape => "circle", # hide
@@ -275,7 +276,7 @@ sheep_rr = Rule(DVS_N, DVS_W);
 
 # #### Test rotation
 ex = @acset_colim yLV begin
-  e::E; s::Sheep; countdown(src(e)) == 0; countdown(tgt(e)) == 0
+  e::E; s::Sheep;
   sheep_loc(s)==src(e); sheep_eng(s)==100; dir(e)==left(sheep_dir(s))
 end;
 
@@ -305,7 +306,6 @@ sheep_fwd_rule = Rule(first(homs(E, s_fwd_l; monic=true)),
 ex = @acset_colim yLV begin
   (e1,e2)::E; (s::Sheep)
   src(e2)==tgt(e1); sheep_loc(s)==src(e1)
-  countdown(src(e1))==0; countdown(tgt(e1))==0; countdown(tgt(e2))==0
   sheep_eng(s)==10
   sheep_dir(s)==dir(e1); dir(e1)==left(dir(e2))
 end
@@ -317,38 +317,45 @@ expected[1, :sheep_eng] = 9
 
 
 # ### Sheep eat grass
-
-s_eat_L = @acset_colim yLV begin s::Sheep; countdown(sheep_loc(s)) == 0 end;
-s_eat_R = @acset_colim yLV begin s::Sheep; countdown(sheep_loc(s)) == 30 end;
-se_rule = Rule(hom(S, s_eat_L), hom(S, s_eat_R); expr=(Eng=[((vₛ,),) -> vₛ + 4],));
+s_eat_N = @acset_colim yLV begin 
+  s::Sheep; t::Time; countdown(t) == sheep_loc(s) 
+end
+s_eat_L = @acset_colim yLV begin s::Sheep;  end;
+s_eat_R = deepcopy(s_eat_L)
+add_parts!(s_eat_R, :Time, 30; countdown=1)
+se_left = hom(G⊕D, s_eat_L; initial=(Direction=1:4,))
+se_right = hom(G⊕D, s_eat_R;  initial=(Direction=1:4,))
+se_nac = hom(s_eat_L, s_eat_N)
+se_rule = Rule(se_left, se_right; 
+               ac=[AppCond(se_nac, false)], expr=(Eng=[((vₛ,),) -> vₛ + 4],));
 
 # #### Sheep eating test
 ex = @acset_colim yLV begin
-  e::E; s::Sheep
-  dir(e)==sheep_dir(s); countdown(src(e))==10; countdown(tgt(e)) == 0
+  e::E; s::Sheep; t::Time
+  dir(e)==sheep_dir(s); countdown(t)==src(e);
   sheep_loc(s)==tgt(e); sheep_eng(s) == 3
 end
 
 expected = copy(ex)
-expected[2,:countdown] = 30
+add_parts!(expected, :Time, 30; countdown=2)
 expected[1,:sheep_eng] = 7
-
 @test is_isomorphic(expected, rewrite(se_rule, ex))
-
 
 # ### Wolves eat sheep
 
 w_eat_l = @acset_colim yLV begin
   s::Sheep; w::Wolf; sheep_loc(s) == wolf_loc(w)
 end;
 
-we_rule = Rule(hom(W⊕D, w_eat_l), id(W⊕D); 
+we_left = hom(G⊕D⊕D, w_eat_l; initial=(Direction=Dict(1=>1, 5=>5),))
+we_right = hom(G⊕D⊕D, D⊕W; initial=(Direction=Dict(1=>1, 5=>5),))
+we_rule = Rule(we_left, we_right; 
                expr=(Eng=[((vₛ, vᵩ),) -> vᵩ + 20],));
 
 # #### Wolf eating test
 ex = @acset_colim yLV begin 
-  (s::Sheep); (w::Wolf); (e::E); 
-  countdown(src(e))==9; countdown(tgt(e))==10
+  (s::Sheep); (w::Wolf); (e::E); (t1,t2,t3)::Time
+  countdown(t1)==src(e); countdown(t2)==src(e); countdown(t3)==tgt(e)
   sheep_dir(s)==left(wolf_dir(w))
   sheep_dir(s)==right(dir(e))
   sheep_eng(s)==3; wolf_eng(w)==16
@@ -363,12 +370,12 @@ rem_part!(expected, :Sheep, 1)
 
 # ### Sheep starvation
 s_die_l = @acset_colim yLV begin s::Sheep; sheep_eng(s) == 0 end;
-sheep_die_rule = Rule(hom(G⊕D, s_die_l), id(G⊕D));
+sheep_die_rule = Rule(hom(G⊕D, s_die_l; any=true), id(G⊕D));
 
 # #### Sheep starvation test
 ex = @acset_colim yLV begin 
-  s::Sheep; w::Wolf
-  countdown(sheep_loc(s))==20; countdown(wolf_loc(w))==10
+  s::Sheep; w::Wolf; (t,t2)::Time
+  countdown(t)==sheep_loc(s); countdown(t2)==wolf_loc(w)
   sheep_eng(s)==0; wolf_eng(w)==10; sheep_dir(s) == right(wolf_dir(w))
 end
 expected = copy(ex)
@@ -383,17 +390,17 @@ s_reprod_r = @acset_colim yLV begin
 end;
 
 sheep_reprod_rule = Rule(
-  hom(G⊕D, S),
-  hom(G⊕D, s_reprod_r);
+  hom(G⊕D, S; any=true),
+  hom(G⊕D, s_reprod_r; any=true);
   expr=(Dir=fill(vs->only(vs) ,2), 
         Eng=fill(vs -> round(Int, vs[1] / 2, RoundUp), 2),)
 );
 
 # #### Reproduction test
 
 ex = @acset_colim yLV begin 
-  s::Sheep; w::Wolf
-  countdown(sheep_loc(s))==20; countdown(wolf_loc(w))==10
+  s::Sheep; w::Wolf; t::Time
+  countdown(t)==sheep_loc(s);
   sheep_eng(s)==10; wolf_eng(w)==20; sheep_dir(s) == right(wolf_dir(w))
 end
 
@@ -410,21 +417,16 @@ can_match(sheep_reprod_rule, m)
 
 # ### Grass increments
 
-g_inc_n = @acset LV begin V=1; countdown=0 end
-
-g_inc_rule = Rule(G; ac=[AppCond(hom(G, g_inc_n), false)],
-                  expr=(Countdown=[((vᵥ,),) -> vᵥ - 1],));
+g_inc_L = @acset_colim yLV begin t::Time end
+rem_time = hom(G, g_inc_L)
+g_inc_rule = Rule(rem_time, id(G));
 
 # #### Grass incrementing test
 ex = @acset_colim yLV begin
-  e::E; countdown(src(e)) == 0; countdown(tgt(e)) == 2
-end
-
-expected = @acset_colim yLV begin
-  e::E; countdown(src(e)) == 0; countdown(tgt(e)) == 1
+  e::E; t::Time; countdown(t) == tgt(e)
 end
 
-@test is_isomorphic(rewrite(g_inc_rule, ex), expected)
+@test is_isomorphic(rewrite(g_inc_rule, ex), E)
 
 # # Adding timers to the rules and making the model
 

src/ABMs.jl

@@ -55,7 +55,7 @@ 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
+  val::Function # ACSetTransformation → hazard_rate
 end
 
 (c::ClosureState)(m::ACSetTransformation) = c.val(m)
@@ -73,9 +73,8 @@ DiscreteHazard(t::Number) = DiscreteHazard(Dirac(t))
 end
 
 """Check if a hazard rate is a simple exponential"""
-is_exp(::Distributions.Exponential) = true
 is_exp(h::ContinuousHazard) = h.val isa Distributions.Exponential
-is_exp(h::AbsHazard) = false
+is_exp(h::AbsTimer) = false
 
 ContinuousHazard(p::Number) = ContinuousHazard(Exponential(p))
 
@@ -193,7 +192,7 @@ function get_hazard(r::RepresentableP, f::ACSetTransformation, ::Float64,
                     h::ContinuousHazard) 
    err = "Representable patterns must have simple exponential rules"
    X = codom(f)
-   is_exp(h.val) ? Exponential(h.val.θ/multiplier(r,X)) : error(err)
+   is_exp(h) ? Exponential(h.val.θ/multiplier(r,X)) : error(err)
 end
 
 const Maybe{T} = Union{Nothing, T}
@@ -476,7 +475,7 @@ function run!(abm::ABM, rt::RuntimeABM, output::Traj;
     new_time = first(next(rt.sampler, rt.tnow, rt.rng))
     if !isempty(abm.dyn) && dt < new_time 
       error("HERE")
-    else  
+    else
       # Get next event + unpack data
       events::Vector{Pair{Int,Maybe{KeyType}}} = pops!(rt) # updates the clock time
       N = length(rt.sampler)

test/ABMs.jl

@@ -30,7 +30,7 @@ rem_loop = ABMRule(:RemLoop, Rule(delete(Graph(1)), id(Graph(1))), DiscreteHazar
 @test rem_loop.pattern_type == RegularP()
 
 # •→• ⇽ • → •
-rem_edge = ABMRule(:RemEdge, Rule(homomorphism(Graph(2), path_graph(Graph, 2); monic=true), 
+rem_edge = ABMRule(:RemEdge, Rule(homomorphism(Graph(2), path_graph(Graph, 2); initial=(V=1:2,)), 
                         id(Graph(2))), 
                    ContinuousHazard(1))
 @test rem_edge.pattern_type == RepresentableP(Dict(:E=>[1]))
@@ -73,7 +73,7 @@ using AlgebraicABMs, AlgebraicRewriting, Catlab
 # Rule: copy a variable
 v = @acset LSet begin X=1; D=1; f=[AttrVar(1)] end
 v2 = @acset LSet begin X=2; D=1; f=[AttrVar(1), AttrVar(1)] end
-dup_vertex = ABMRule(Rule(id(v), homomorphism(v, v2)), DiscreteHazard(1.))
+dup_vertex = ABMRule(Rule(id(v), homomorphism(v, v2; initial=(X=[1],))), DiscreteHazard(1.))
 
 # Dynamics: for an individual variable, it grows linearly w/ time
 flow = ABMFlow(v, RawODE([_ -> 1.0]), :Grow, [], [(:D => 1)])