diff --git a/Project.toml b/Project.toml
index d8bba62..dc31d84 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,7 @@ name = "AlgebraicRewriting"
 uuid = "725a01d3-f174-5bbd-84e1-b9417bad95d9"
 license = "MIT"
 authors = ["Kris Brown <kris@topos.institute>"]
-version = "0.3.5"
+version = "0.3.6"
 
 [deps]
 ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
@@ -25,7 +25,7 @@ AlgebraicRewritingLuxorExt = "Luxor"
 
 [compat]
 ACSets = "0.2.20"
-Catlab = "0.16.11"
+Catlab = "0.16.16"
 CompTime = "0.1"
 DataMigrations = "0.0.3,0.1"
 DataStructures = "0.17, 0.18"
diff --git a/README.md b/README.md
index 6737cb7..26a98f9 100644
--- a/README.md
+++ b/README.md
@@ -9,6 +9,10 @@ This package defines DPO, SPO, and SqPO for C-Sets, as described in
 [[Brown 2022](https://arxiv.org/abs/2111.03784)]. See the [documentation](https://algebraicjulia.github.io/AlgebraicRewriting.jl/dev/) for more details.
 
 
+To locally build the [documentation](https://algebraicjulia.github.io/AlgebraicRewriting.jl/dev) and the literate code examples, run the following in the command line:
+```
+julia --project=docs -e "using AlgebraicRewriting, LiveServer; servedocs(literate_dir=\"docs/literate\",skip_dir=\"docs/src/generated\")"
+```
 
 ## NOTE
 This library is currently under active development, and so is not yet at a
diff --git a/docs/literate/full_demo.jl b/docs/literate/full_demo.jl
index d025e4f..ff172d6 100644
--- a/docs/literate/full_demo.jl
+++ b/docs/literate/full_demo.jl
@@ -1,6 +1,6 @@
 # # Full Demo
 
-using AlgebraicRewriting, Catlab, AlgebraicPetri, DataMigrations
+using AlgebraicRewriting, Catlab, DataMigrations, AlgebraicPetri
 using Test
 
 
@@ -101,7 +101,7 @@ rule2 = Rule(@migration(SchRulel, SchGraph, begin
     l => begin
       v => src(e)
     end
-  end), yG)
+  end), yG);
 
 # We can also rewrite without a match (and let it pick an arbitrary match).
 
@@ -114,7 +114,8 @@ rule2 = Rule(@migration(SchRulel, SchGraph, begin
 rule_spo = Rule{:SPO}(l, r)  # (same data as before)
 
 @test length(get_matches(rule_spo, G)) == 4 # there are now four matches
-res = rewrite(rule_spo, G)
+m = get_matches(rule_spo, G)[1]
+res = rewrite_match(rule_spo, m)
 to_graphviz(res)
 @test is_isomorphic(res, path_graph(Graph, 3) ⊕ R)
 
@@ -127,19 +128,23 @@ to_graphviz(res)
 
 L = Graph(1)
 I = Graph(2)
-R = path_graph(Graph, 2)
+R = path_graph(Graph, 2);
 
 #=
-We can use automated homomorphism search to reduce the tedium of specifying data manually. In this case, there is a unique option.
+We can use automated homomorphism search to reduce the tedium of specifying data 
+manually. In this case, there is a unique option. In general, `homomorphism`
+will throw an error if there is *more* than one homomorphism.
 =#
 
-l = homomorphism(I, L)
+l = homomorphism(I, L);
 
 #=
-There are many constraints we can put on the search, such as being monic.
+There are many constraints we can put on the search, such as being monic. Here
+there are two monic homomorphisms (sending vertices 1 and 2 to (1,2) and (2,1)),
+so we add the keyword `any=true` to avoid throwing an error.
 =#
 
-r = homomorphism(I, R; monic=true)
+r = homomorphism(I, R; monic=true, any=true)
 
 rule_sqpo = Rule{:SqPO}(l, r) # same data as before)
 
@@ -154,56 +159,64 @@ to_graphviz(res; prog="neato")
 # # 4. PBPO+
 
 #=
-PBPO+ requires not merely a span but also additional data for L and K which can be thought of as type graphs. The graph G that we rewrite will be typed over the L' type graph to determine how it is rewritten. 
+PBPO+ requires a span just like the other kinds of rewriting.
 =#
 
 L = Graph(1)
-K = Graph(2)
-l = homomorphism(K, L)
-r = id(K)
-
-# We allow edges into and out of the matched vertex as well as edges 
-# between the vertices incident to the matched vertex 
-L′ = @acset Graph begin
-  V = 3
-  E = 6
-  src = [1, 1, 1, 2, 3, 3]
-  tgt = [1, 2, 3, 3, 3, 1]
-end
-tl = ACSetTransformation(L, L′; V=[2]) # 2 is the matched vertex
-to_graphviz(L′; node_labels=true)
+K = R = Graph(2)
+l, r = homomorphism(K,L), id(K);
 
-# The outneighbors of the matched vertex are duplicated (an edge connects the 
-# old ones to the new ones) and the matched vertex is duplicated. The new copy 
-# of the matched vertex points at the new ones. It does not have any inneighbors.
+# However, it also requires more data. The graph G that we rewrite will be typed *over* the L' type graph which controls how various parts of the context (not merely the matched pattern) are rewritten. 
+
+# 1 = root of the deep copy, 2 = children of #1, 3 = everything else 
+L′ = @acset Graph begin V=3;E=5;src=[1,2,3,3,3];tgt=[2,2,1,2,3] end
+to_graphviz(L′; node_labels=true)
 
-K′ = @acset Graph begin
-  V = 5
-  E = 9
-  src = [1, 1, 1, 2, 3, 3, 3, 4, 5]
-  tgt = [1, 2, 3, 3, 3, 1, 5, 5, 5]
+tl = ACSetTransformation(L,L′;V=[1])
+K′ = @acset Graph begin V=5;E=9 ;
+  src=[1,2,3,3,4,3,5,3,3];tgt=[2,2,3,2,5,5,5,1,4] 
+end
+tk = ACSetTransformation(K,K′;V=[1,4])
+l′ = homomorphism(K′,L′; initial=(V=[1,2,3,1,2],));
+
+"""Given a match L → G, compute what the typing map G → L' should be"""
+function get_adherence(m::ACSetTransformation) 
+  root, G, descendents  = only(collect(m[:V])), codom(m), Set()
+  queue = [root]
+  while !isempty(queue)
+    nxt = pop!(queue)
+    union!(descendents, outneighbors(G,nxt))
+    union!(queue, outneighbors(G,nxt))
+  end
+  return (V = map(parts(codom(m),:V)) do v_G 
+    if     v_G == root       return 1
+    elseif v_G ∈ descendents return 2
+    else                     return 3
+    end 
+  end,)
 end
-tk = ACSetTransformation(K, K′; V=[2, 4])
-to_graphviz(K′; node_labels=true)
 
-l′ = homomorphism(K′, L′; initial=(V=[1, 2, 3, 2, 3],))
+rule = PBPORule(l, r, tl, tk, l′; adherence=get_adherence);
 
-prule = PBPORule(l, r, tl, tk, l′)
+# Think of the following graph as a file system 
 
-# Apply to an example vertex (#3) with two inneighbors and one outneighbor.
-G = @acset Graph begin
-  V = 4
-  E = 5
-  src = [1, 1, 2, 3, 4]
-  tgt = [2, 3, 3, 4, 4]
+G = @acset Graph begin V=8; E=8; 
+  src=[1,1,2,2,3,4,4,5]; tgt=[2,3,4,5,6,5,7,8] 
 end
+
 to_graphviz(G; node_labels=true)
 
-m = get_match(prule, G; initial=(V=[3],) => Dict())
+# Executing this rule (forcing the pattern to match at vertex 2) performs a 
+# "deepcopy" operation, copying vertex 2 and everything underneath it.
+
+expected = @acset Graph begin V=13;E=14;
+  src=[7,7,7,1,1,3,3,4,2,2,10,10,11,8]; tgt=[1,2,8,3,4,5,4,6,10,11,12,11,13,9] 
+end
+
+@test is_isomorphic(expected, rewrite(rule, G; initial=(V=[2],)))
+
+to_graphviz(expected; node_labels=true)
 
-res = rewrite_match(prule, m)
-# V1 is copied to V2. Outneighbor V5 (w/ loop) is copied to V6, creating an edge
-to_graphviz(res; node_labels=true)
 
 # # 5. Generalizing Graphs
 
@@ -213,7 +226,6 @@ Any data structure which implements the required functions we need can, in princ
 Here we'll do rewriting in graphs sliced over •⇆•, which is isomorphic to the category of (whole-grain) Petri nets, with States and Transitions.
 =#
 
-
 function graph_slice(s::Slice)
   h = s.slice
   V, E = collect.([h[:V], h[:E]])
@@ -222,10 +234,7 @@ function graph_slice(s::Slice)
   nS, nT, nI, nO = length.([S, T, I, O])
   findS, findT = [x -> findfirst(==(x), X) for X in [S, T]]
   to_graphviz(@acset AlgebraicPetri.PetriNet begin
-    S = nS
-    T = nT
-    I = nI
-    O = nO
+    S = nS; T = nT; I = nI; O = nO
     is = findS.(g[I, :src])
     it = findT.(g[I, :tgt])
     ot = findT.(g[O, :src])
@@ -236,18 +245,18 @@ end;
 # This is the graph we are slicing over.
 
 two = @acset Graph begin
-  V = 2
-  E = 2
-  src = [1, 2]
-  tgt = [2, 1]
+  V = 2; E = 2; src = [1, 2]; tgt = [2, 1]
 end
 
-# Define a rule which deletes a [T] -> S edge
+to_graphviz(two)
+
+# Define a rule which deletes a [T] -> S edge. Start with the pattern, L.
 
 L_ = path_graph(Graph, 2)
 L = Slice(ACSetTransformation(L_, two, V=[2, 1], E=[2])) # [T] ⟶ (S)
 graph_slice(L)
 
+# Then define I and R 
 I_ = Graph(1)
 I = Slice(ACSetTransformation(I_, two, V=[2])) # [T]
 R_ = Graph(2)
@@ -274,11 +283,13 @@ While the vast majority of functionality is focused on ACSets at the present mom
 #=
 We can construct commutative diagrams with certain edges left unspecified or marked with ∀ or ∃. If only one edge is left free, we can treat the diagram as a boolean function which tests whether the morphism makes the specified paths commute (or not commute). This generalizes positive/negative application conditions and lifting conditions, but because those are most common there are constructors AppCond and LiftCond to make these directly.
 
+```
          ∀ 
   [↻•]   →  ?
     ↓    ↗ ∃ ↓ 
   [↻•⟶•]  → [↻•⟶•⟵•↺]
-
+```
+`
 Every vertex with a loop also has a map to the vertex marked by the bottom map.
 =#
 
@@ -291,25 +302,19 @@ end
 
 v = homomorphism(t, looparr)
 loop_csp = @acset Graph begin
-  V = 3
-  E = 4
-  src = [1, 3, 1, 3]
-  tgt = [1, 3, 2, 2]
+  V = 3; E = 4; src = [1, 3, 1, 3]; tgt = [1, 3, 2, 2]
 end
-b = homomorphism(looparr, loop_csp; monic=true)
+b = homomorphism(looparr, loop_csp; initial=(V=[2,1],))
 constr = LiftCond(v, b)
 
-@test !apply_constraint(constr, homomorphism(t, loop_csp))
+@test !apply_constraint(constr, homomorphism(t, loop_csp; initial=(V=[1],)))
 @test apply_constraint(constr, b)
 
 # We can combining constraints with logical combinators.
 
 # match vertex iff it has 2 or 3 self loops
 one, two, three, four, five = [@acset(Graph, begin
-  V = 1
-  E = n
-  src = 1
-  tgt = 1
+  V = 1; E = n; src = 1; tgt = 1 
 end) for n in 1:5]
 
 c2 = AppCond(homomorphism(Graph(1), two); monic=true)         # PAC
@@ -343,16 +348,12 @@ L = @acset_colim yWG begin
 end
 I = WeightedGraph{Int}(2)
 R = @acset WeightedGraph{Int} begin
-  V = 2
-  E = 1
-  Weight = 1
-  src = 1
-  tgt = 2
-  weight = [AttrVar(1)]
+  V = 2; E = 1; Weight = 1
+  src = [1]; tgt = [2]; weight = [AttrVar(1)]
 end
 
-l = homomorphism(I, L; monic=true)
-r = homomorphism(I, R; monic=true)
+l = homomorphism(I, L; initial=(V=1:2,))
+r = homomorphism(I, R; initial=(V=1:2,))
 rule = Rule(l, r; monic=[:E], expr=Dict(:Weight => [xs -> xs[1] + xs[2]]))
 
 G = @acset WeightedGraph{Int} begin
@@ -363,12 +364,10 @@ G = @acset WeightedGraph{Int} begin
   weight = [10, 20, 100]
 end
 
-@test rewrite(rule, G) == @acset WeightedGraph{Int} begin
-  V = 1
-  E = 2
-  src = 1
-  tgt = 1
-  weight = [30, 100]
+m = get_matches(rule,G)[1]
+@test rewrite_match(rule, m) == @acset WeightedGraph{Int} begin
+  V = 1; E = 2
+  src = [1]; tgt = [1]; weight = [30, 100]
 end
 
 # # 8. Graph processes
@@ -381,7 +380,7 @@ via analyzing the colimit of all the partial maps induced by the rewrites.
 
 using AlgebraicRewriting.Processes: RWStep, find_deps
 
-G0, G1, G2, G3 = Graph.([0, 1, 2, 3])
+G0, G1, G2, G3 = Graph.([0, 1, 2, 3]);
 # Delete a node
 Rule1 = Span(create(G1), id(G0));
 # Merge two nodes
@@ -389,7 +388,7 @@ Rule2 = Span(id(G2), homomorphism(G2, G1));
 # Add a node
 Rule3 = Span(id(G0), create(G1))
 
-R1, R2, R3 = [Rule(l, r) for (l, r) in [Rule1, Rule2, Rule3]]
+R1, R2, R3 = [Rule(l, r) for (l, r) in [Rule1, Rule2, Rule3]];
 
 # # 9. Trajectory
 
@@ -397,32 +396,28 @@ R1, R2, R3 = [Rule(l, r) for (l, r) in [Rule1, Rule2, Rule3]]
 M1 = create(G2)
 CM1 = ACSetTransformation(G1, G3; V=[3])
 Pmap1 = Span(id(G2), ACSetTransformation(G2, G3; V=[1, 2]))
-RS1 = RWStep(Rule3, Pmap1, M1, CM1)
+RS1 = RWStep(Rule3, Pmap1, M1, CM1);
 
 # Step 2: merge node 2 and 3 to yield a G2
 M2 = ACSetTransformation(G2, G3; V=[2, 3])
 CM2 = ACSetTransformation(G1, G2; V=[2])
 Pmap2 = Span(id(G3), ACSetTransformation(G3, G2; V=[1, 2, 2]))
-RS2 = RWStep(Rule2, Pmap2, M2, CM2)
+RS2 = RWStep(Rule2, Pmap2, M2, CM2);
 
 # Step 3: delete vertex 1 
 M3 = ACSetTransformation(G1, G2; V=[1])
 CM3 = create(G1)
 Pmap3 = Span(ACSetTransformation(G1, G2; V=[2]), id(G1))
-RS3 = RWStep(Rule1, Pmap3, M3, CM3)
+RS3 = RWStep(Rule1, Pmap3, M3, CM3);
 
 
-steps = [RS1, RS2, RS3]
+steps = [RS1, RS2, RS3];
 
 g = find_deps(steps)
 to_graphviz(g; node_labels=true)
 
-expected = @acset Graph begin
-  V = 3
-  E = 1
-  src = 1
-  tgt = 2
-end
+# Confirm this what we expect
+expected = @acset Graph begin V = 3; E = 1; src = 1; tgt = 2 end
 @test expected == g
 
 # Interface that just uses rules and match morphisms:
@@ -432,4 +427,4 @@ end
 σ₃ = ACSetTransformation(G3, G3; V=[3, 1, 2])
 
 g′ = find_deps([R3 => M1, R2 => M2 ⋅ σ₃, R1 => M3 ⋅ σ₂])
-@test g′ == g
\ No newline at end of file
+@test g′ == g
diff --git a/docs/literate/game_of_life.jl b/docs/literate/game_of_life.jl
index cc3098c..309f53b 100644
--- a/docs/literate/game_of_life.jl
+++ b/docs/literate/game_of_life.jl
@@ -6,17 +6,21 @@ This is a demonstration of the game of life as an agent-based model.
 We start with importing some libraries.
 =#
 
-using AlgebraicRewriting
-using Catlab, Catlab.Graphs, Catlab.CategoricalAlgebra, Catlab.Theories
+using AlgebraicRewriting, Catlab
 import Catlab.Graphics: to_graphviz
 using Catlab.Graphics.Graphviz: Attributes, Statement, Node, Edge, Digraph
-using PrettyTables
-using Luxor
+using PrettyTables, Luxor
 
 #=
-The game of life has two rules: one which turns living things dead, and one that brings dead things to life. We model the terrain as a symmetric graph: cells are vertices. Neighboring cells have edges between them.
-
-Implementation wise, if we are going to update cells one at a time, we must keep track of two bits of information (the cell's living status for the *current* timestep and whether it will be alive in the *next* timestep). Thus we need helper rule to overwrite the "current" life status with the "next" life status at the end of each timestep.
+The game of life has two rules: one which turns living things dead, and one that 
+brings dead things to life. We model the terrain as a symmetric graph: cells are 
+vertices. Neighboring cells have edges between them.
+
+Implementation wise, if we are going to update cells one at a time, we must keep 
+track of two bits of information (the cell's living status for the *current* 
+timestep and whether it will be alive in the *next* timestep). Thus we need 
+helper rule to overwrite the "current" life status with the "next" life status 
+at the end of each timestep.
 =#
 
 # # Ontology
@@ -161,7 +165,7 @@ data migration, `F`.
 =#
 
 A = Life(1)
-view_life(homomorphism(F(A), init)) |> println
+view_life(homomorphism(F(A), init; any=true)) |> println
 
 #=
 We must also work with miniature game states that are *not* grids in order for 
@@ -244,7 +248,7 @@ rewrite rules.
 BirthP1 = living_neighbors(3) # must have 3 neighbors
 BirthN1 = living_neighbors(4) # forbid the cell to have 4 neighbors
 BirthN2 = Curr() # forbid the cell to be alive (i.e. it's currently dead)
-BP1, BN1, BN2 = homomorphism.(Ref(Life(1)), [BirthP1, BirthN1, BirthN2])
+BP1, BN1, BN2 = homomorphism.(Ref(Life(1)), [BirthP1, BirthN1, BirthN2]; initial=(V=[1],))
 bac = [AppCond(BP1; monic=true), AppCond.([BN1, BN2], false; monic=true)...]
 Birth = Rule(id(Life(1)), to_next(); ac=bac);
 
@@ -254,7 +258,7 @@ PersistR = @acset Life begin
 end
 PersistP1 = living_neighbors(2; alive=true)
 PersistN1 = living_neighbors(4; alive=true)
-DR, DP1, DN1 = homomorphism.(Ref(Curr()), [PersistR, PersistP1, PersistN1])
+DR, DP1, DN1 = homomorphism.(Ref(Curr()), [PersistR, PersistP1, PersistN1]; initial=(V=[1],))
 pac = [AppCond(DP1; monic=true), AppCond(DN1, false; monic=true)]
 Persist = Rule(id(Curr()), DR; ac=pac);
 
diff --git a/docs/literate/lotka_volterra.jl b/docs/literate/lotka_volterra.jl
index ebd63ba..7ad1596 100644
--- a/docs/literate/lotka_volterra.jl
+++ b/docs/literate/lotka_volterra.jl
@@ -491,7 +491,7 @@ sheep_reprod_rule = Rule(
 );
 
 sheep_reprod = RuleApp(:reproduce, sheep_reprod_rule,
-  id(S), hom(S, s_reprod_r)) |> tryrule;
+  id(S), hom(S, s_reprod_r; any=true)) |> tryrule;
 
 # #### Reproduction test
 
@@ -541,8 +541,9 @@ begin
     sheep_loc = 2; sheep_eng = [3]; sheep_dir = [:N]
     countdown = [0, 10, 2]; dir = fill(:N, 2)
   end
-  @test is_isomorphic(rewrite(g_inc_rule, ex), expected)
-  rewrite!(g_inc_rule, ex)
+  m = get_matches(g_inc_rule, ex)[1]
+  @test is_isomorphic(rewrite_match(g_inc_rule, m), expected)
+  rewrite_match!(g_inc_rule, m)
   @test is_isomorphic(ex, expected)
 end;
 
diff --git a/src/categorical_algebra/CSets.jl b/src/categorical_algebra/CSets.jl
index 7789ff4..6689dd8 100644
--- a/src/categorical_algebra/CSets.jl
+++ b/src/categorical_algebra/CSets.jl
@@ -194,7 +194,7 @@ function pushout_complement(pair::ComposablePair{<:ACSet, <:TightACSetTransforma
         G[g_components[o](val), f]
       catch e 
         vᵢ = findfirst(==(AttrVar(v)), k_components[at])
-        m[at](AttrVar(l[at](vᵢ)))
+        m[at](l[at](AttrVar(vᵢ)))
       end 
     end)
   end
diff --git a/src/incremental/IncrementalConstraints.jl b/src/incremental/IncrementalConstraints.jl
index 82133c4..902663a 100644
--- a/src/incremental/IncrementalConstraints.jl
+++ b/src/incremental/IncrementalConstraints.jl
@@ -150,13 +150,13 @@ function can_match(constr::IncConstraints, m::ACSetTransformation;
 
   if pac 
     for ac in constr.pac
-      isnothing(extend_morphism(m, ac.m;  monic=ac.monic)) && return false
+      isnothing(extend_morphism(m, ac.m; any=true, monic=ac.monic)) && return false
     end
   end
   
   if nac 
     for ac in constr.nac
-      isnothing(extend_morphism(m, ac.m;  monic=ac.monic)) || return false
+      isnothing(extend_morphism(m, ac.m; any=true, monic=ac.monic)) || return false
     end
   end
   
diff --git a/src/rewrite/Inplace.jl b/src/rewrite/Inplace.jl
index 6e27e9a..8296610 100644
--- a/src/rewrite/Inplace.jl
+++ b/src/rewrite/Inplace.jl
@@ -288,7 +288,7 @@ end
 
 rewrite_match!(r::Rule{:DPO}, ::Nothing; kw...) = nothing
 
-rewrite!(r::Rule{:DPO}, G; initial=nothing, random=false, kw...) =
+rewrite!(r::Rule{:DPO}, G; initial=(;), random=false, kw...) =
   rewrite_match!(r, get_match(r, G; initial, random); kw...)
 
 end
diff --git a/src/rewrite/PBPO.jl b/src/rewrite/PBPO.jl
index 7355eec..1f0dd5e 100644
--- a/src/rewrite/PBPO.jl
+++ b/src/rewrite/PBPO.jl
@@ -116,30 +116,15 @@ The "strong match" condition we enforce is that: tl⁻¹(α(A)) = a⁻¹(A). Thi
 we can deduce precisely what m is by looking at α.
 
 """
-function get_matches(rule::PBPORule, G::ACSet;  initial=nothing, 
-                     α_unique=true, random=false, n=-1, kw...)
-  S = acset_schema(G)
+function get_matches(rule::PBPORule, G::ACSet;  initial=(;), 
+                     α_unique=true, random=false, take=nothing)
   res = [] # Quadruples of of (m, Labs, abs, α)
   L = codom(left(rule))
 
-  # Process the initial constraints for match morphism and typing morphism
-  if isnothing(initial)
-    matchinit, typinit = (;), (;)
-  elseif initial isa Union{NamedTuple,AbstractDict}
-    matchinit, typinit = Dict(pairs(initial)), (;)
-  elseif length(initial)==2
-    matchinit, typinit = [Dict(pairs(x)) for x in initial]
-  else 
-    error("Unexpected type for `initial` keyword: $initial")
-  end  
-
   # Search for each match morphism
-  backtracking_search(L, G; monic=rule.monic, initial=NamedTuple(matchinit),
-                      random) do m
-    m_seen = false # keeps track if α_unique is violated for each new m
-    if all(ac->apply_constraint(ac, m), rule.acs)
-      @debug "m:  $([k=>collect(v) for (k,v) in pairs(components(m))])"
-      
+  backtracking_search(L, G; monic=rule.monic, initial=initial, random) do ms
+    for m in ms       
+      all(ac -> apply_constraint(ac, m), rule.acs) || continue
       # Construct partially-abtract version of G. Labs: L->A and abs: A->G 
       Labs, abs = partial_abstract(m)
       A = codom(Labs)
@@ -150,46 +135,49 @@ function get_matches(rule::PBPORule, G::ACSet;  initial=nothing,
         # Return nothing if failure
         if !isnothing(init)
           αs = homomorphisms(A, codom(rule.tl); initial=init)
-          # Also return nothing if the result is not unique
-          if length(αs) ==1 
-            push!(res, deepcopy((m, Labs, abs, only(αs))))
-          end 
+          # Also return nothing if the result isn't unique
+          length(αs) == 1 && push!(res, deepcopy((m, Labs, abs, only(αs))))
         end
       else 
         # Search for adherence morphisms: A -> L′
         init = extend_morphism_constraints(rule.tl, Labs)
         isnothing(init) && return false
-        backtracking_search(A, codom(rule.tl); initial=init, kw...) do α
-          @debug "\tα: ", [k=>collect(v) for (k,v) in pairs(components(α))] 
-
-          # Check strong match condition
-          strong_match = all(types(S)) do o
-            prt = o ∈ ob(S) ? identity : AttrVar
-            all(prt.(parts(A,o))) do i 
-              p1 = preimage(rule.tl[o],α[o](i))
-              p2 = preimage(Labs[o], i)
-              p1 == p2
-            end
-          end
-          if strong_match && all(lc -> apply_constraint(lc, α), rule.lcs)
-            all(is_natural, [m, Labs, abs, α]) || error("Unnatural match")
-            if m_seen  error("Multiple α for a single match $m") end 
-            @debug "\tSUCCESS"
-            push!(res, deepcopy((m, Labs, abs, α)))
-            m_seen |= α_unique
-            return length(res) == n
-          else
-            @debug "\tFAILURE (strong $strong_match)"
-            return false
-          end
+        kwargs = if α_unique 
+          (max = 1,) 
+        else 
+          (take = isnothing(take) ? nothing : take-length(res),)
+        end
+        for α in homomorphisms(A, codom(rule.tl); initial=init, kwargs...,
+                               filter=test_adherence(rule, Labs)) 
+          push!(res, deepcopy((m, Labs, abs, α)))
+          length(res) == take && return true
         end
       end
     end
-    return length(res) == n
+    return false
   end 
   return res
 end
 
+""" See `get_matches(::PBPORule, ::ACSet)` """
+function test_adherence(rule::PBPORule, Labs::ACSetTransformation)
+  A = codom(Labs)
+  S = acset_schema(A)
+  """ For use in the homomorphism search for adherence morphisms """
+  function filter(α::ACSetTransformation)::Bool
+    # Check strong match condition
+    strong_match = all(types(S)) do o
+      prt = o ∈ ob(S) ? identity : AttrVar
+      all(prt.(parts(A, o))) do i 
+        p1 = preimage(rule.tl[o],α[o](i))
+        p2 = preimage(Labs[o], i)
+        p1 == p2
+      end
+    end
+    return strong_match && all(lc -> apply_constraint(lc, α), rule.lcs)
+  end
+end
+
 
 """
 This construction addresses the following problem: ideally when we 'abstract' 
diff --git a/src/rewrite/Utils.jl b/src/rewrite/Utils.jl
index 8ec765d..5a36a4f 100644
--- a/src/rewrite/Utils.jl
+++ b/src/rewrite/Utils.jl
@@ -144,23 +144,27 @@ has_comp(monic::Vector{Symbol}, c::Symbol) = c ∈ monic
 """
 Returns nothing if the match is acceptable for rewriting according to the
 rule, otherwise returns the reason why it should be rejected
+
+homsearch = if we know ahead of time that m was obtained m via automatic hom 
+            search, then we do not need to make certain checks
 """
-function can_match(r::Rule{T}, m; initial=Dict()) where T
+function can_match(r::Rule{T}, m; homsearch=false, initial=Dict()) where T
   S = acset_schema(dom(m))
-  for k in ob(S)
-    if has_comp(r.monic, k) && !is_monic(m[k])
-      return ("Match is not injective", k, m[k])
+  if !homsearch
+    for k in ob(S)
+      if has_comp(r.monic, k) && !is_monic(m[k])
+        return ("Match is not injective", k, m[k])
+      end
     end
-  end
-  for (k, vs) in collect(initial)
-    errs = check_initial(vs, collect(m[k]))
-    if !isempty(errs)
-      return ("Initial condition violated", k, errs)
+    for (k, vs) in collect(initial)
+      errs = check_initial(vs, collect(m[k]))
+      if !isempty(errs)
+        return ("Initial condition violated", k, errs)
+      end
     end
+    is_natural(m) || return ("Match is not natural", m)
   end
 
-  is_natural(m) || return ("Match is not natural", m)
-
   if T == :DPO
     gc = gluing_conditions(ComposablePair(r.L, m))
     if !isempty(gc)
@@ -183,7 +187,7 @@ function can_match(r::Rule{T}, m; initial=Dict()) where T
 end
 
 """Get one match (if any exist) otherwise return """
-get_match(args...; kw...) = let x = get_matches(args...; n=1, kw...);
+get_match(args...; kw...) = let x = get_matches(args...; take=1, kw...);
   isempty(x) ? nothing : only(x) end 
 
 """
@@ -193,32 +197,16 @@ This function has the same behavior as the generic `get_matches`, but it is
 more performant because we do not have to query all homomorphisms before finding 
 a valid match, in case n=1. 
 """
-function get_matches(r::Rule{T}, G::ACSet; initial=nothing,
-                     random=false, n=-1) where T
-  initial = isnothing(initial) ? Dict() : initial
-
-  hs = []
-  backtracking_search(codom(r.L), G; monic=r.monic, initial=NamedTuple(initial), 
-                      random) do h 
-    cm = can_match(r, h)
-    if isnothing(cm)
-      push!(hs, h)
-      return length(hs) == n # we stop the search Hom(L,G) when this holds
-    else
-      @debug "$([k => collect(v) for (k,v) in pairs(components(h))]): $cm"
-      return false
-    end
-  end 
-  return hs
-end
+get_matches(r::Rule, G::ACSet; kw...) =
+  homomorphisms(codom(r.L), G; kw..., monic=r.monic, 
+                filter= m -> isnothing(can_match(r, m; homsearch=true)))
 
 """If not rewriting ACSets, we have to compute entire Hom(L,G)."""
-function get_matches(r::Rule{T}, G; initial=nothing, random=false, n=-1) where T 
-  initial = NamedTuple(isnothing(initial) ? Dict() : initial)
-  ms = homomorphisms(codom(left(r)), G; monic=r.monic, initial, random)
+function get_matches(r::Rule, G; take=nothing, kw...)
+  ms = homomorphisms(codom(left(r)), G; kw..., monic=r.monic)
   res = []
   for m in ms 
-    if (n < 0 || length(res) < n) && isnothing(can_match(r, m))
+    if (isnothing(take) || length(res) < take) && isnothing(can_match(r, m))
       push!(res, m)
     end
   end
@@ -276,7 +264,7 @@ function rewrite_match_maps end  # to be implemented for each T
 """    rewrite(r::Rule, G; kw...)
 Perform a rewrite (automatically finding an arbitrary match) and return result.
 """
-function rewrite(r::AbsRule, G; initial=nothing, random=false, kw...)
+function rewrite(r::AbsRule, G; initial=(;), random=false, kw...)
   m = get_match(r, G; initial, random)
   isnothing(m) ? nothing : rewrite_match(r, m; kw...)
 end
diff --git a/test/categorical_algebra/CSets.jl b/test/categorical_algebra/CSets.jl
index 43cd14c..0a8f8fc 100644
--- a/test/categorical_algebra/CSets.jl
+++ b/test/categorical_algebra/CSets.jl
@@ -31,7 +31,7 @@ const LSet = LFinSetType{Symbol}
 
 I = @acset LSet begin X=1; D=1; f=[AttrVar(1)] end
 G = @acset LSet begin X=2; f=[:x,:y] end
-f = homomorphism(I,G)
+f = homomorphism(I,G; initial=(X=[1],))
 
 kg = last(pushout_complement(id(I),f))
 @test dom(kg) == @acset LSet begin X=2; D=1; f=[AttrVar(1),:y] end
diff --git a/test/incremental/Benchmark.jl b/test/incremental/Benchmark.jl
index 74e3a4d..68afe13 100644
--- a/test/incremental/Benchmark.jl
+++ b/test/incremental/Benchmark.jl
@@ -21,10 +21,10 @@ while true
   I = rand(path_graph.(Graph, 2:4))
   NV=200
   start = erdos_renyi(Graph, NV, 2*NV)
-  l = homomorphism(I, L; monic=true)
-  r = homomorphism(I, R; monic=true)
+  l = homomorphism(I, L; monic=true, any=true)
+  r = homomorphism(I, R; monic=true, any=true)
   isnothing(r) && continue
-  m = homomorphism(I, start)
+  m = homomorphism(I, start; any=true)
   isnothing(m) && continue
   res = rewrite_match_maps(Rule(id(I), r), m);
   (pl, pr), rmap = get_pmap(:DPO, res), get_rmap(:DPO, res);
@@ -50,14 +50,14 @@ DDS(i::Int) = @acset DDS begin X=i; Φ=[rand(1:i) for _ in 1:i] end
 
 while true
   L, R, I, A, B = DDS.([5, 5, 5, 3, 3])
-  l1,l2,r1,r2 = hs = [homomorphism(x...; monic=true) for x in 
+  l1,l2,r1,r2 = hs = [homomorphism(x...; monic=true, any=true) for x in 
                       [(A,L),(A,I),(B,I),(B,R)]]
   all(!isnothing, hs) || continue
   (_, l), (r, _) = pushout(l1,l2), pushout(r1,r2)
   rand_rule = Rule(l, r)
 
   start, pattern = DDS(2000), DDS(5)
-  m = homomorphism(codom(l), start)
+  m = homomorphism(codom(l), start; any=true)
   (!isnothing(m) && isnothing(can_match(rand_rule, m))) || continue 
 
   res = rewrite_match_maps(rand_rule, m)
diff --git a/test/incremental/IncrementalCC.jl b/test/incremental/IncrementalCC.jl
index f0a779f..c27649a 100644
--- a/test/incremental/IncrementalCC.jl
+++ b/test/incremental/IncrementalCC.jl
@@ -13,7 +13,7 @@ G2 = Graph(2)
                                                                   #  • ⇉ •
 e, ee = path_graph.(Graph, 2:3)                                   #   ↘ ↙
 A = @acset Graph begin V=3; E=4; src=[1,1,1,2]; tgt=[2,2,3,3] end #    •
-A_rule = Rule(id(e), homomorphism(e, A));
+A_rule = Rule(id(e), homomorphism(e, A; initial=(E=[1],)));
 
 # Empty edge case
 #----------------
@@ -35,8 +35,8 @@ del, add = rewrite!(hset, A_rule, m)
 @test length.(match_vect(hset)) == [3,0,6]
 @test validate(hset)
 
-m = homomorphism(e, state(hset); monic=true)
-rewrite!(hset, A_rule)
+m = homomorphism(e, state(hset); initial=(V=1:2,E=[1]))
+rewrite!(hset, A_rule, m)
 @test validate(hset)
 @test length.(match_vect(hset)) == [3, 0, 6, 0, 8]
 @test !haskey(hset, 3=>7)
@@ -53,8 +53,8 @@ roundtrip = IncCCHomSet(IncSumHomSet(hset));
 #----------------
 tri = @acset Graph begin V=3;E=3;src=[1,1,2];tgt=[3,2,3]end
 X = @acset Graph begin V=2; E=2; src=[1,2]; tgt=[2,2] end
-omap = homomorphism(e, X)
-r = homomorphism(e, tri)
+omap = homomorphism(e, X; initial=(V=1:2,))
+r = homomorphism(e, tri; initial=(V=1:2,))
 hset = IncHomSet(ee, [r], X);
 addition!(hset, r, omap)
 @test validate(hset)
@@ -75,7 +75,8 @@ del = delete(Graph(1))
 mset = IncHomSet(Graph(1), [del], G2⊕T; nac=[del]);
 @test length(keys(mset)) == 2
 M_rule = Rule(id(Graph(1)), delete(Graph(1)); ac=[AppCond(del, false)])
-rewrite!(mset, M_rule)
+m = ACSetTransformation(Graph(1), G2⊕T; V=[1])
+rewrite!(mset, M_rule, m)
 @test length(keys(mset)) == 1
 rewrite!(mset, M_rule)
 @test length(keys(mset)) == 0
@@ -84,15 +85,15 @@ rewrite!(mset, M_rule)
 # Application conditions: NAC adding morphisms during deletion!
 #--------------------------------------------------------------
 del = homomorphism(⊕(Graph[fill(T, 2); Graph(1)]), 
-                   state(mset); monic=true) # delete one loop
+                   state(mset); initial=(V=1:3,)) # delete one loop
 deletion!(mset, del)
 @test length(keys(mset)) == 1
 del = homomorphism(⊕(Graph[T; G2]), 
-                   state(mset); monic=true) # delete another loop
+                   state(mset); initial=(V=1:3,)) # delete another loop
 deletion!(mset, del)
 @test length(keys(mset)) == 2
 
-del = homomorphism(Graph(3), state(mset); monic=true) # delete another loop
+del = homomorphism(Graph(3), state(mset); initial=(V=1:3,)) # delete another loop
 deletion!(mset, del)
 @test length(keys(mset)) == 3
 
@@ -101,7 +102,7 @@ deletion!(mset, del)
 edge_loop = @acset Graph begin V=2; E=2; src=[1,1]; tgt=[1,2] end
 to_edge_loop = homomorphism(e, edge_loop; monic=true)
 # rem edge, not if src has loop
-r = Rule(homomorphism(G2, e; monic=true), id(G2);
+r = Rule(homomorphism(G2, e; initial=(V=1:2,)), id(G2);
          ac=[AppCond(to_edge_loop, false; monic=true)]);
 
 mset = IncHomSet(r, edge_loop);
@@ -114,7 +115,7 @@ rewrite!(mset, r)
 # Application conditions: PAC removing morphisms during deletion!
 #----------------------------------------------------------------
 # Remove edge, only if src has loop (no monic constraint on PAC)
-r = Rule(homomorphism(G2, e; monic=true), id(G2);
+r = Rule(homomorphism(G2, e; initial=(V=1:2,)), id(G2);
          ac=[AppCond(to_edge_loop)]);
 mset = IncHomSet(r, edge_loop ⊕ e);
 m1, m2 = get_matches(r, state(mset)) # first one removes the loop
diff --git a/test/incremental/IncrementalSum.jl b/test/incremental/IncrementalSum.jl
index 5cd04d3..b706a86 100644
--- a/test/incremental/IncrementalSum.jl
+++ b/test/incremental/IncrementalSum.jl
@@ -9,7 +9,7 @@ using AlgebraicRewriting.Incremental.IncrementalHom: runtime, state
                                                                   #  • ⇉ •
 e, ee = path_graph.(Graph, 2:3)                                   #   ↘ ↙
 A = @acset Graph begin V=3; E=4; src=[1,1,1,2]; tgt=[2,2,3,3] end #    •
-A_rule = Rule(id(e), homomorphism(e, A));
+A_rule = Rule(id(e), homomorphism(e, A; initial=(E=[1],)));
 start = @acset Graph begin V=3; E=3; src=[1,2,3]; tgt=[2,3,3] end
 
 hset = IncHomSet(ee ⊕ e, [A_rule.R], start);
@@ -18,7 +18,7 @@ hset = IncHomSet(ee ⊕ e, [A_rule.R], start);
 @test !haskey(hset, [2=>2, 1=>2])
 @test length(keys(hset)) == 9
 @test hset[[1=>3,1=>3]] == hset[9]
-del, add = rewrite!(hset, A_rule, homomorphisms(e, start)[2]);
+del, add = rewrite!(hset, A_rule, homomorphism(e, start; initial=(V=2:3,)));
 
 @test isempty(del)
 
@@ -28,7 +28,8 @@ del, add = rewrite!(hset, A_rule, homomorphisms(e, start)[2]);
 @test validate(hset)
 
 @test Set(matches(hset)) == Set(homomorphisms(ee ⊕ e, state(hset)))
-rewrite!(hset, A_rule);
+m = ACSetTransformation(e, state(hset); V=1:2, E=[1])
+rewrite!(hset, A_rule, m);
 @test validate(hset)
 @test Set(matches(hset)) == Set(homomorphisms(ee ⊕ e, state(hset)))
 
@@ -42,7 +43,8 @@ roundtrip = IncSumHomSet(IncCCHomSet(hset));
 hset = IncHomSet(Graph(1) ⊕ e, [A_rule.R], start);
 rewrite!(hset, A_rule, homomorphisms(e, start)[2]);
 @test validate(hset)
-rewrite!(hset, A_rule)
+m = ACSetTransformation(e, state(hset); V=1:2, E=[1])
+rewrite!(hset, A_rule, m)
 @test validate(hset)
 @test length(keys(hset)) == 45
 
@@ -56,7 +58,7 @@ DDS(phi::Vector{Int}) = @acset DDS begin X=(length(phi)); Φ=phi end
 
 p2 = DDS([2,2])
 p22 = p2 ⊕ p2
-r = homomorphism(p22, DDS([2,2,4,4,4]); monic=true)
+r = homomorphism(p22, DDS([2,2,4,4,4]); initial=(X=1:4,))
 hset = IncHomSet(p22, [r], p22);
 rewrite!(hset, Rule(id(p22), r), id(p22))
 @test validate(hset)
diff --git a/test/rewrite/Constraints.jl b/test/rewrite/Constraints.jl
index 2e13671..1397204 100644
--- a/test/rewrite/Constraints.jl
+++ b/test/rewrite/Constraints.jl
@@ -50,10 +50,10 @@ Every vertex with a loop also has a map to the vertex marked by the bottom map.
 t = terminal(Graph)|>apex
 v = homomorphism(t, looparr)
 loop_csp = @acset Graph begin V=3;E=4; src=[1,3,1,3]; tgt=[1,3,2,2] end 
-b = homomorphism(looparr, loop_csp; monic=true)
+b = homomorphism(looparr, loop_csp; initial=(V=1:2,))
 constr = LiftCond(v, b)
 
-@test !apply_constraint(constr,homomorphism(t, loop_csp))
+@test !apply_constraint(constr,homomorphism(t, loop_csp; initial=(V=[1],)))
 @test apply_constraint(constr,b)
 
 
@@ -68,8 +68,8 @@ Every vertex mapped into the leftmost vertex must also have an outgoing edge to
 a vertex that is is mapped into the middle vertex.
 """
 
-constr = LiftCond(homomorphism(Graph(1),p2), 
-                  homomorphism(p2, loop_csp; monic=true))
+constr = LiftCond(ACSetTransformation(Graph(1),p2; V=[1]), 
+                  homomorphism(p2, loop_csp; initial=(V=[1,2],)))
 
 G = @acset Graph begin V=3; E=3; src=[1,1,3]; tgt=[1,2,3] end
 h1,h2,h3,h4 = homomorphisms(G, loop_csp; initial=(V=Dict(1=>1),))
diff --git a/test/rewrite/DPO.jl b/test/rewrite/DPO.jl
index 30314d8..b436a5b 100644
--- a/test/rewrite/DPO.jl
+++ b/test/rewrite/DPO.jl
@@ -1,8 +1,6 @@
 module TestDPO 
 
-using Test
-using Catlab, Catlab.CategoricalAlgebra, Catlab.Programs, Catlab.Graphs
-using AlgebraicRewriting
+using Test, Catlab, AlgebraicRewriting
 
 # Graphs
 ########
@@ -243,13 +241,15 @@ R = @acset SSet begin
   tgt=[2,3,4,4,3]
 end
 edge = @acset SSet begin E=1; V=2; src=[1]; tgt=[2] end
-edge_left = homomorphism(edge, L; initial=Dict([:V=>[1,3]]))
-edge_left_R = homomorphism(edge, R; initial=Dict([:V=>[1,3]]))
-edge_right = homomorphism(edge, L; initial=Dict([:V=>[2,4]]))
+edge_left = only(homomorphisms(edge, L; initial=Dict([:V=>[1,3]])))
+edge_left_R = only(homomorphisms(edge, R; initial=Dict([:V=>[1,3]])))
+edge_right = only(homomorphisms(edge, L; initial=Dict([:V=>[2,4]])))
 G = apex(pushout(edge_left, edge_right))
-r = Rule(homomorphism(I, L; monic=true), homomorphism(I, R; monic=true);
+r = Rule(homomorphism(I, L; initial=(V=1:4,)), 
+         homomorphism(I, R; initial=(V=1:4,));
          monic=true)
-res =  rewrite(r, G)
+m = get_matches(r, G)[1]
+res =  rewrite_match(r, m)
 expected = apex(pushout(edge_left_R, edge_right))
 @test !is_isomorphic(res, G) # it changed
 @test is_isomorphic(res, expected)
@@ -287,19 +287,20 @@ I = WeightedGraph{Int}(2)
 R = @acset WeightedGraph{Int} begin V=2; E=1; Weight=1; src=1; tgt=2; 
                                     weight=[AttrVar(1)] end
 
-l = homomorphism(I,L; monic=true)
-r = homomorphism(I,R; monic=true)
+l = homomorphism(I,L; initial=(V=1:2,))
+r = homomorphism(I,R; initial=(V=1:2,))
 rule = Rule(l, r; monic=[:E], expr=Dict(:Weight=>[xs->xs[1]+xs[2]]))
 
 G = @acset WeightedGraph{Int} begin V=1; E=3; src=1; tgt=1; 
                                     weight=[10,20,100] end
-
-@test rewrite(rule, G) == @acset WeightedGraph{Int} begin 
+m = get_matches(rule, G)[1]
+@test rewrite_match(rule, m) == @acset WeightedGraph{Int} begin 
   V=1; E=2; src=1; tgt=1; weight=[30, 100] end
 
 # or, introduce free variables 
 rule = Rule(l, r; monic=[:E], freevar=true)
-@test rewrite(rule, G) == @acset WeightedGraph{Int} begin 
+m = get_matches(rule, G)[1]
+@test rewrite_match(rule, m) == @acset WeightedGraph{Int} begin 
   V=1; E=2; Weight=1; src=1; tgt=1; weight=[AttrVar(1), 100] end
 
 # Rewriting with induced equations between AttrVars in the rule
@@ -312,13 +313,14 @@ const Foo = AbsFoo{Bool}
 L = @acset Foo begin X=2; f=[false, false] end 
 I = @acset Foo begin X=2;D=2; f=AttrVar.(1:2) end 
 R = @acset Foo begin X=2; f=[false, true]end 
-rule = Rule(homomorphism(I, L; monic=[:X]), homomorphism(I, R; monic=[:X]))
+rule = Rule(homomorphism(I, L; initial=(X=1:2,)), 
+            homomorphism(I, R; initial=(X=1:2,)))
 
 # we cannot match both X of L to the same part in G because this would yield 
 # an inconsistent result
 @test length(get_matches(rule, L)) == 2
-
-res = rewrite(rule, L)
+m = get_matches(rule, L)[1]
+res = rewrite_match(rule, m)
 @test is_isomorphic(res, R)
 
 # Now with arbitrary functions
diff --git a/test/rewrite/Inplace.jl b/test/rewrite/Inplace.jl
index 8579412..fefba9b 100644
--- a/test/rewrite/Inplace.jl
+++ b/test/rewrite/Inplace.jl
@@ -37,11 +37,6 @@ prog = compile_rewrite(toggle)
 @test F == rewrite(toggle, T)
 
 m= homomorphism(Pat, T)
-is_natural(m)
-can_match(toggle, m)
-
-get_match(toggle, T)
-
 rewrite_match!(toggle, m; prog);
 @test T == F
 
@@ -66,7 +61,7 @@ G = @acset MADIntGraph begin V=1; E=3; src=1; tgt=1; weight=[10,20,100] end
 # to_graphviz(G; edge_labels=:weight)
 
 # Rule
-l, r = homomorphism.(Ref(MADIntGraph(2)), [L, R]; monic=true)
+l, r = homomorphism.(Ref(MADIntGraph(2)), [L, R]; initial=(V=1:2,))
 
 plus(xs) = xs[1] + xs[2]
 
@@ -75,7 +70,8 @@ rule = Rule(l, r; monic=[:E], expr=Dict(:Weight=>[plus]))
 # Apply rewrite
 prog = compile_rewrite(rule)
 
-f = rewrite!(rule, G)
+m = get_matches(rule, G)[1]
+f = rewrite_match!(rule, m)
 
 @test f.components[:Weight](AttrVar(1)) == 30
 
diff --git a/test/rewrite/PBPO.jl b/test/rewrite/PBPO.jl
index 962db4f..7e6f929 100644
--- a/test/rewrite/PBPO.jl
+++ b/test/rewrite/PBPO.jl
@@ -1,8 +1,6 @@
 module TestPBPO 
 
-using Test
-using AlgebraicRewriting
-using Catlab
+using Test, AlgebraicRewriting, Catlab
 
 using AlgebraicRewriting.Rewrite.PBPO: partial_abstract
 
@@ -50,7 +48,7 @@ G = @acset Graph begin V=5;E=11;
   src=[1,1,1,2,2,2,3,4,5,5,5]; tgt=[2,4,5,1,3,3,2,5,3,3,5] 
 end
 
-res = rewrite(rule,G; initial=(V=[1,2,3],)=>(V=[1,2,3,4,4],))
+res = rewrite(rule, G; initial=(V=[1,2,3],))
 expected = @acset Graph begin V=6;E=10;
   src=[1,3,3,4,5,6,6,6,6,6];tgt=[3,4,5,5,5,3,3,4,5,6] 
 end
@@ -155,7 +153,7 @@ G = @acset Graph begin V=8; E=8;
   src=[1,1,2,2,3,4,4,5]; tgt=[2,3,4,5,6,5,7,8] 
 end
 
-init = (initial=Dict(:V=>[2])=>Dict(),)
+init = (initial=(V=[2],),)
 
 @test length(get_matches(rule_no_condition, G; α_unique=false, init...)) > 1
 @test_throws ErrorException get_matches(rule_no_condition, G; init...) 
@@ -172,12 +170,11 @@ end
 #############################################################
 L = apex(terminal(Graph))
 K = R = Graph(1)
-l = homomorphism(K,L)
-r = homomorphism(K,R)
+l, r = homomorphism.(Ref(K), [L, R])
 L′ = L ⊕ L 
 K′ = K ⊕ L 
-tl = homomorphism(L,L′)
-tk = homomorphism(K,K′)
+tl = homomorphism(L, L′; initial=(V=[1],))
+tk = ACSetTransformation(K,K′; V=[1])
 l′ = homomorphism(K′,L′; initial=(V=[1,2],)) 
 rule = PBPORule(l,r,tl,tk,l′)
 
@@ -197,7 +194,7 @@ l = homomorphism(K,L)
 r = homomorphism(K,R)
 L′ = @acset Graph begin V=2; E=3; src=[1,2,2]; tgt=[1,2,1] end 
 K′ = @acset Graph begin V=2; E=2; src=[2,2]; tgt=[2,1] end 
-tl = homomorphism(L,L′)
+tl = homomorphism(L,L′; initial=(V=[1],))
 tk = ACSetTransformation(K,K′; V=[1])
 l′ = homomorphism(K′,L′; initial=(V=[1,2],))
 rule = PBPORule(l,r,tl,tk,l′)
@@ -213,12 +210,12 @@ expected = @acset Graph begin V=3; E=4; src=[1,2,1,2]; tgt=[1,2,2,3] end
 L = @acset WG begin V=2; E=1; Weight=1; src=1; tgt=2; weight=[AttrVar(1)] end
 K = WG(2)
 R = WG(1)
-l = homomorphism(K, L; monic=true)
+l = homomorphism(K, L; initial=(V=[1,2],))
 r = homomorphism(K, R)
 L′ = @acset WG begin V=4; E=10; Weight=10;
   src=[1,1,2,3,3,3,3,4,4,4]; tgt=[2,4,4,1,4,3,2,2,3,4]; weight=AttrVar.(1:10)
 end
-tl = homomorphism(L, L′)
+tl = homomorphism(L, L′; initial=(V=[1,2],))
 K′ = @acset WG begin V=4; E=9; Weight=9;
   src=[1,2,3,3,3,3,4,4,4]; tgt=[4,4,1,4,3,2,2,3,4]; weight=AttrVar.(1:9)
 end
@@ -237,7 +234,7 @@ ac = AppCond(homomorphism(L,loop), false) # cannot bind pattern to loop
            [   •    ]
 
 """
-lc = LiftCond(homomorphism(R, L), # vertical
+lc = LiftCond(homomorphism(R, L; initial=(V=[1],)), # vertical
               homomorphism(L, L′; initial=(E=[4],)))
 
 kx = Any[fill(nothing, 9)...]
@@ -251,13 +248,13 @@ expected = @acset WG begin V=4; E=6; src=[2,2,2,3,3,4]; tgt=[1,3,4,1,1,1];
   weight=[6., 5., 6., 8., 7., 9.]
 end
 
-init = Dict(:V => [1, 2]) => Dict()
+initial = Dict(:V => [1, 2])
 
-@test length(get_matches(rule, G; initial=init))==1
+@test length(get_matches(rule, G; initial))==1
 @test isempty(homomorphisms(codom(left(rule)), G; monic=true))
 
-@test only(get_matches(rule, G; initial=init)) == get_match(rule, G;initial=init)
-@test is_isomorphic(expected, rewrite(rule, G; initial=init))
+@test only(get_matches(rule, G; initial)) == get_match(rule, G;initial)
+@test is_isomorphic(expected, rewrite(rule, G; initial))
 
 # Test canonization: TODO
 
diff --git a/test/rewrite/Representable.jl b/test/rewrite/Representable.jl
index ca5f72a..0567c9f 100644
--- a/test/rewrite/Representable.jl
+++ b/test/rewrite/Representable.jl
@@ -1,8 +1,6 @@
 module TestRepresentable 
 
-using Test
-using AlgebraicRewriting, DataMigrations
-using Catlab.Graphs, Catlab.CategoricalAlgebra, Catlab.Programs
+using Test, AlgebraicRewriting, DataMigrations, Catlab
 
 yWG = yoneda_cache(WeightedGraph{Float64}; clear=true);
 
@@ -42,10 +40,10 @@ r = Rule(d, yWG; semantics=:SPO, expr=(Weight=[vs->0.],), monic=true)
 G = path_graph(WeightedGraph{Float64}, 2)
 add_edges!(G, [2,2], [2,2])
 G[:weight] = [-1.,10,10]
-m = get_match(r, G)
+m = get_matches(r, G)[1]
 expected = @acset WeightedGraph{Float64} begin 
   V=3; E=2; src=[1,2]; tgt=[2,3]; weight=[-1.,0.] 
 end
-@test is_isomorphic(rewrite(r, G), expected)
+@test is_isomorphic(rewrite_match(r, m), expected)
 
 end # module 
diff --git a/test/rewrite/SPO.jl b/test/rewrite/SPO.jl
index 37d02fa..4a20876 100644
--- a/test/rewrite/SPO.jl
+++ b/test/rewrite/SPO.jl
@@ -1,18 +1,17 @@
 module TestSPO
 
-using Test
-using Catlab, Catlab.Graphs, Catlab.CategoricalAlgebra
-using AlgebraicRewriting
+using Test, Catlab, AlgebraicRewriting
 
 
 # Removing edges 
 #---------------
 p2, g2 = path_graph(Graph, 2), Graph(2)
-f = homomorphism(g2, p2; monic=true)
+f = homomorphism(g2, p2; initial=(V=1:2,))
 r = Rule{:SPO}(f, id(g2))
 r2 = Rule{:SPO}(create(Graph(1)), id(Graph()))
 @test rewrite(r, p2) == Graph(2)
-@test rewrite(r2, p2) == Graph(1)
+m = get_match(r2, p2)
+@test rewrite_match(r2, m) == Graph(1)
 
 
 
@@ -41,12 +40,12 @@ spr = rewrite_match(Rule{:SPO}(ka,kb), ac)
 # Semisimplicial sets
 #####################
 @present ThSemisimplicialSet(FreeSchema) begin
-(V,E,T) :: Ob
-(d1,d2,d3)::Hom(T,E)
-(src,tgt) :: Hom(E,V)
-compose(d1, src) == compose(d2, src)
-compose(d1, tgt) == compose(d3, tgt)
-compose(d2, tgt) == compose(d3, src)
+  (V,E,T) :: Ob
+  (d1,d2,d3)::Hom(T,E)
+  (src,tgt) :: Hom(E,V)
+  compose(d1, src) == compose(d2, src)
+  compose(d1, tgt) == compose(d3, tgt)
+  compose(d2, tgt) == compose(d3, src)
 end
 @acset_type SSet(ThSemisimplicialSet)
 
@@ -61,37 +60,37 @@ end
 
 L = quadrangle  # We defined quadrilateral above.
 I = @acset SSet begin
-E=4; V=4
-src=[1,1,2,3]
-tgt=[2,3,4,4]
+  E=4; V=4
+  src=[1,1,2,3]
+  tgt=[2,3,4,4]
 end
 R = @acset SSet begin
-T=2; E=5; V=4
-d1=[2,3]
-d2=[1,5]
-d3=[5,4]
-src=[1,1,2,3,2]
-tgt=[2,3,4,4,3]
+  T=2; E=5; V=4
+  d1=[2,3]
+  d2=[1,5]
+  d3=[5,4]
+  src=[1,1,2,3,2]
+  tgt=[2,3,4,4,3]
 end
 edge = @acset SSet begin E=1; V=2; src=[1]; tgt=[2] end
 
 Tri = @acset SSet begin
-T=1; E=3; V=3;
-d1=[1]; d2=[2]; d3=[3];
-src=[1,1,2]; tgt=[3,2,3]
+  T=1; E=3; V=3;
+  d1=[1]; d2=[2]; d3=[3];
+  src=[1,1,2]; tgt=[3,2,3]
 end
 
 
 r = Rule{:SPO}(homomorphisms(edge, Tri)[2], id(edge))
 r_dpo = Rule(r.L, r.R)
 
-m = homomorphism(Tri, quadrangle)
+m = homomorphism(Tri, quadrangle; initial=(V=[1,2,4],))
 
 # This does not make sense for DPO
 @test !can_pushout_complement(ComposablePair(r.L, m))
 @test_throws ErrorException rewrite_match_maps(r_dpo, m; check=true)
-@test is_isomorphic(rewrite_match(r,m),
-                  @acset SSet begin E=2; V=3; src=1; tgt=[2,3] end)
+@test is_isomorphic(rewrite_match(r, m),
+                    @acset SSet begin E=2; V=3; src=1; tgt=[2,3] end)
 
 # Attributed rewrite
 #-------------------
diff --git a/test/rewrite/SqPO.jl b/test/rewrite/SqPO.jl
index 47d7c60..4ca1b32 100644
--- a/test/rewrite/SqPO.jl
+++ b/test/rewrite/SqPO.jl
@@ -1,9 +1,7 @@
 
 module TestSqPO 
 
-using Test
-using Catlab, Catlab.Graphs, Catlab.CategoricalAlgebra
-using AlgebraicRewriting 
+using Test, Catlab, AlgebraicRewriting 
 using AlgebraicRewriting.Rewrite.SqPO: final_pullback_complement
 
 # Sesqui Pushout Tests
diff --git a/test/schedules/Eval.jl b/test/schedules/Eval.jl
index 80ab366..e398600 100644
--- a/test/schedules/Eval.jl
+++ b/test/schedules/Eval.jl
@@ -42,7 +42,7 @@ N["•→•"] = ar
 Dot, A = Symbol.([N[g1],N[ar]]) 
 
 av = RuleApp(:add_vertex, Rule(id(z), create(g1)))
-g2ar = homomorphism(g2, ar; monic=true)
+g2ar = homomorphism(g2, ar; initial=(V=1:2,))
 de = loop_rule(RuleApp(:del_edge, Rule(g2ar, id(g2))))
 coin = uniform(2, z)
 sched = coin ⋅ (tryrule(av) ⊗ id([z])) ⋅ merge_wires(z) ⋅ de
@@ -194,7 +194,8 @@ inc_E_ = @acset LG begin Cell=2; V=6; E=7; Life=1; Eng=2
   live=[true,AttrVar(1)]; eng=AttrVar.(1:2)
 end
 inc_E = Rule(id(inc_E_), id(inc_E_); expr=(Eng=[es->es[1],es->es[2]+1],))
-inc_E_rule = RuleApp(:incE, inc_E, homomorphism(Cell, inc_E_)) |> tryrule
+inc_E_rule = RuleApp(
+  :incE, inc_E, homomorphism(Cell, inc_E_; initial=(V=[1,2,4,5],))) |> tryrule
 
 # Assemble a schedule 
 sched = agent(inc_E_rule, Cell, ret=Cell)