diff --git a/Project.toml b/Project.toml
index 707269f10..e0939e186 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,7 +28,6 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
 Colors = "0.12"
-CompTime = "0.1"
 Compose = "0.7, 0.8, 0.9"
 DataStructures = "0.17, 0.18"
 GeneralizedGenerated = "0.2, 0.3"
diff --git a/src/categorical_algebra/CSets.jl b/src/categorical_algebra/CSets.jl
index a80b780bd..d2ed7f387 100644
--- a/src/categorical_algebra/CSets.jl
+++ b/src/categorical_algebra/CSets.jl
@@ -8,11 +8,12 @@ export ACSetTransformation, CSetTransformation,
   isomorphism, isomorphisms, is_isomorphic,
   generate_json_acset, parse_json_acset, read_json_acset, write_json_acset,
   generate_json_acset_schema, parse_json_acset_schema,
-  read_json_acset_schema, write_json_acset_schema, acset_schema_json_schema
+  read_json_acset_schema, write_json_acset_schema, acset_schema_json_schema,
+  uncurry, curry, ACSetCat
 
 using Base.Iterators: flatten
 using Base.Meta: quot
-using DataStructures: OrderedDict
+using DataStructures: OrderedDict, DefaultDict
 using StructEquality
 import JSON
 using Reexport
@@ -30,8 +31,11 @@ using ..FreeDiagrams, ..Limits, ..Subobjects, ..FinSets, ..FinCats
 import ..Limits: limit, colimit, universal
 import ..Subobjects: Subobject, implies, ⟹, subtract, \, negate, ¬, non, ~
 import ..Sets: SetOb, SetFunction, TypeSet
+import ..FinCats: FinDomFunctor, components, is_natural, FinCatGraph, FinCatPresentation, homomorphisms
+using ...Graphs
+using ..FinCats: normalize, Path
 import ..FinSets: FinSet, FinFunction, FinDomFunction, force, predicate, is_monic, is_epic
-import ..FinCats: FinDomFunctor, components, is_natural
+import ..FinCats: FinDomFunctor, components, is_natural, FinTransformationMap, FinDomFunctorMap
 
 # Sets interop
 ##############
@@ -134,17 +138,24 @@ end
 const ACSetDomCat = FinCats.FinCatPresentation{
   Symbol, Union{FreeSchema.Ob,FreeSchema.AttrType},
   Union{FreeSchema.Hom,FreeSchema.Attr,FreeSchema.AttrType}}
+const FinSetCat = TypeCat{SetOb,FinDomFunction{Int}}
 
-""" Wrapper type to interpret attributed C-set as a functor.
-"""
 @struct_hash_equal struct ACSetFunctor{ACS<:ACSet} <:
-    Functor{ACSetDomCat,TypeCat{SetOb,FinDomFunction{Int}}}
+    Functor{ACSetDomCat,FinSetCat}
   acset::ACS
+  eqs::Vector{Pair}
 end
-FinDomFunctor(X::ACSet) = ACSetFunctor(X)
+FinDomFunctor(X::ACSet; eqs=Pair[]) = ACSetFunctor(X, eqs)
 
-dom(F::ACSetFunctor) = FinCat(Presentation(F.acset))
-codom(F::ACSetFunctor) = TypeCat{SetOb,FinDomFunction{Int}}()
+function dom(F::ACSetFunctor)
+  p = Presentation(F.acset)
+  for (l,r) in F.eqs
+    add_equation!(p, l, r)
+  end
+  FinCat(p)
+end
+
+codom(F::ACSetFunctor) = FinSetCat()
 
 Categories.do_ob_map(F::ACSetFunctor, x) = SetOb(F.acset, functor_key(x))
 Categories.do_hom_map(F::ACSetFunctor, f) = SetFunction(F.acset, functor_key(f))
@@ -165,13 +176,19 @@ function (::Type{ACS})(F::FinDomFunctor) where ACS <: ACSet
   return X
 end
 
+function (C::Type{ACS})(F::FinTransformationMap) where ACS <: ACSet
+  Cd, CCd = C(dom(F)), C(codom(F))
+  return CSetTransformation(Cd, CCd; components(F)...)
+end
+
+
 """ Copy parts from a set-valued `FinDomFunctor` to an `ACSet`.
 """
 function ACSetInterface.copy_parts!(X::ACSet, F::FinDomFunctor)
   pres = presentation(dom(F))
   added = Dict(Iterators.map(generators(pres, :Ob)) do c
     c = nameof(c)
-    c => add_parts!(X, c, length(ob_map(F, c)::FinSet{Int}))
+    c => add_parts!(X, c, length(ob_map(F, Symbol(c))::FinSet{Int}))
   end)
   for f in generators(pres, :Hom)
     dom_parts, codom_parts = added[nameof(dom(f))], added[nameof(codom(f))]
@@ -234,6 +251,11 @@ end
 components(α::ACSetTransformation) = α.components
 force(α::ACSetTransformation) = map_components(force, α)
 
+FinTransformationMap(f::ACSetTransformation; eqs=Pair[]) =
+  FinTransformationMap(components(f),
+                    FinDomFunctor(dom(f); eqs=eqs),
+                    FinDomFunctor(codom(f); eqs=eqs))
+
 """ Transformation between C-sets.
 
 Recall that a C-set homomorphism is a natural transformation: a transformation
@@ -738,6 +760,94 @@ partial_assignments(x::AbstractVector) =
 in_hom(S, c) = [dom(S,f) => f for f in hom(S) if codom(S,f) == c]
 out_hom(S, c) = [f => codom(S,f) for f in hom(S) if dom(S,f) == c]
 
+
+
+"""
+Convert a FinFunctor between graph FinCats and add labels from
+FinCatPresentations
+"""
+function grph_fun_to_pres_fun(F, X,Y)
+  hs = map(F.hom_map) do p
+    if isempty(p.edges)
+      return id(ob_generators(Y)[p.src])
+    else
+      compose(hom_generators(Y)[p.edges])
+    end
+  end
+  FinFunctor(Dict(zip(ob_generators(X), ob_generators(Y)[F.ob_map])),
+             Dict(zip(hom_generators(X),hs)), X, Y)
+end
+
+# Convert presentation inputs into graph inputs, opposite for outputs
+homomorphisms(X::FinCatPresentation, Y::FinCatPresentation; kwargs...) =
+  [grph_fun_to_pres_fun(F,X,Y) for F in
+   homomorphisms(FinCatGraph(X),FinCatGraph(Y); kwargs...)]
+
+"""
+Search for FinFunctors between FinCats. I.e. map objects onto objects and
+generating morphisms onto paths. Note that there can be infinite paths in a
+FinCat with a cyclic underlying graph (although equations make make it finite).
+`nmax` restricts the paths found to those which pass over the same object at
+most `nmax` times.
+
+Because morphism equality modulo codomain equations is not implemented, there
+may be "duplicate" morphisms in the results. To be safe, we error if the
+`monic`/`iso` constraints are used when the codomain has equations.
+ob and hom maps can be initialized. Hom maps can be restricted to be a
+particular length (e.g. 1 means a generator must map onto another generator, not
+a composite).
+"""
+function homomorphisms(gX::FinCatGraph, gY::FinCatGraph; n_max::Int=3,
+                       monic_obs=false, epic_obs=false, init_obs=nothing,
+                       init_homs=nothing, hom_lens=nothing)
+
+  y_pths = map(collect(enumerate_paths_cyclic(gY.graph; n_max=n_max))) do (k,v)
+    k => unique(h->normalize(gY, Path(h, k[1], k[2])), v)
+  end |> Dict
+
+  yGrph = Graph(nv(gY.graph))
+
+  for (i,j) in Iterators.product(vertices(gY.graph), vertices(gY.graph))
+    if !isempty(y_pths[(i, j)])
+      add_edge!(yGrph, i, j)
+    end
+  end
+  no,nh = [length(f(gX)) for f in [ob_generators,hom_generators]]
+  init_obs = isnothing(init_obs) ? fill(nothing,no) : init_obs
+  init_homs = isnothing(init_homs) ? fill(nothing,nh) : init_homs
+  hom_lens = isnothing(hom_lens) ? fill(nothing,nh) : hom_lens
+  res = []
+  kwargs = Dict{Symbol,Any}(:initial=>(V=Dict([
+    i=>v for (i, v) in enumerate(init_obs) if !isnothing(v)]),))
+  if monic_obs kwargs[:monic] = [:V] end
+  if epic_obs error("depends on is_surjective PR") end
+  for h in homomorphisms(gX.graph, yGrph; kwargs...)
+    om = collect(h[:V])
+    pths = map(zip(collect(h[:E]), init_homs, hom_lens)) do (e, ih, hl)
+      p = y_pths[(src(yGrph,e),tgt(yGrph,e))]
+      # Apply init_homs and hom_lens constraints to possible paths
+      if !isnothing(hl)
+        p = filter(z->length(z)==hl, p)
+      end
+      if !isnothing(ih)
+        p = filter(x->is_hom_equal(gY,x,ih), p)
+      end
+      return p
+    end
+    # This should be done in a backtracking style
+    for combo in Iterators.product(pths...) # pick a path for each edge in X
+      hm = map(enumerate(combo)) do (ei,z)
+        isempty(z) ? id(gY, om[src(gX.graph,ei)]) : z
+      end
+      F = FinFunctor((V=om, E=hm), gX, gY)
+      if is_functorial(F; check_equations=true)
+        push!(res, F)
+      end
+    end
+  end
+  return res
+end
+
 # Limits and colimits
 #####################
 
@@ -1093,6 +1203,123 @@ end
     end...))
 end
 
+
+# Tensor-hom adjunction (currying of diagrams in C-Set)
+#######################################################
+const ACSetCat{S} = TypeCat{S, ACSetTransformation}
+
+"""    uncurry(d::FinFunctor{D, ACSetCat{S}}) where {D,S}
+Undoing currying on objects of a functor category. C->D->Set ==> (CxD)->Set
+"""
+function uncurry(d::FinFunctor{D, ACSetCat{S}}) where {D,S}
+  shapelim = product(FinCatPresentation[dom(d), FinCat(Presentation(S))])
+  shape_ind, part_ind = legs(shapelim)
+  asl = apex(shapelim)
+  omap = Dict(map(ob_generators(asl)) do o
+    x = ob_map(shape_ind, o)
+    y = ob_map(part_ind, o)
+    o => FinSet(ob_map(d, x), Symbol(y))
+  end)
+
+  hmap = Dict(map(hom_generators(asl)) do o
+    x = hom_map(shape_ind, o)
+    y = hom_map(part_ind, o)
+    if first(typeof(x).parameters) == :id
+      o => FinFunction(ob_map(d, only(x.args)), Symbol(y))
+    elseif first(typeof(y).parameters) == :id
+      o => hom_map(d, x)[Symbol(only(y.args))]
+    else
+      error("x $x $(typeof(x)) y $y $(typeof(y))")
+    end
+  end)
+
+  FinDomFunctor(omap,hmap,asl,FinSetCat())
+end
+
+"""
+Currying a FinFunctor into Set: (CxD)->Set ==> C->D->Set
+
+An example FinDomFunctor (in the original curried format) is required.
+"""
+function curry(d::FinDomFunctor{D1, FinSetCat},
+                 old_d::FinDomFunctor{D2, ACSetCat{S}}) where {D1,D2,S}
+  # Recover schema for d as a product, not just the apex
+  shapelim = product([dom(old_d), FinCat(Presentation(S))])
+  asl = apex(shapelim)
+  shape_ind, part_ind = legs(shapelim)
+
+  cset_type = typeof(first(old_d.ob_map)[2])
+  omap = Dict(map(ob_generators(dom(old_d))) do o
+    x = Base.invokelatest(cset_type)
+    for o_ in ob_generators(asl)
+      if ob_map(shape_ind, o_) == o
+        add_parts!(x, Symbol(ob_map(part_ind, o_)), length(ob_map(d, o_)))
+      end
+    end
+    for h in hom_generators(asl)
+      h_ = hom_map(shape_ind, h)
+      if h_ == id(o)
+        set_subpart!(x, Symbol(hom_map(part_ind, h)), collect(hom_map(d, h)))
+      end
+    end
+    o => x
+  end)
+  hmap = Dict(map(hom_generators(dom(old_d))) do h
+    comps = Dict()
+    for h_ in hom_generators(asl)
+      if hom_map(shape_ind, h_) == h
+        comps[Symbol(only(hom_map(part_ind, h_).args))] = hom_map(d, h_)
+      end
+    end
+    dom_, codom_ = [omap[get(h)] for get in [dom, codom]]
+    h => ACSetTransformation(dom_,codom_; comps...)
+  end)
+  FinDomFunctor(omap,hmap,dom(old_d),ACSetCat{S}())
+end
+
+"""    uncurry(d::FinFunctor{D, ACSetCat{S}}) where {D,S}
+Uncurrying on morphisms of a functor category with an ACSetCat as codom
+"""
+function uncurry(ϕ::FinTransformationMap{D, ACSetCat{S}}) where {D,S}
+  cur_d, cur_cd = uncurry.([dom(ϕ), codom(ϕ)])
+  shapelim = product([dom(dom(ϕ)), FinCat(Presentation(S))])
+  shape_ind, part_ind = legs(shapelim)
+  comps = Dict(map(ob_generators(apex(shapelim))) do o
+    oshape, opart = Symbol(shape_ind(o)), Symbol(part_ind(o))
+    Symbol(o) => components(ϕ)[oshape][opart]
+  end)
+  FinTransformationMap(comps,cur_d,cur_cd)
+end
+
+"""    curry(d::FinTransformationMap, old_d::FinTransformationMap{D, ACSetCat{S}}) where {D, S}
+Currying on morphisms of a functor category with an ACSetCat as codom
+An example FinDomFunctor (in the original curried format) is required.
+"""
+function curry(d::FinTransformationMap,
+                 old_d::FinDomFunctor{D, ACSetCat{S}}
+                 ) where {D, S}
+  # Recover schema for d as a product, not just the apex
+  shapelim = product([dom(old_d), FinCat(Presentation(S))])
+  shape_ind, part_ind = legs(shapelim)
+
+  αcomps = Dict(o => DefaultDict{Symbol,Vector{Int}}(()->Int[])
+                for o in Symbol.(ob_generators(dom(old_d))))
+  for o in (ob_generators(apex(shapelim)))
+    dic = αcomps[Symbol(ob_map(shape_ind, o))]
+    dic[Symbol(ob_map(part_ind, o))] = collect(components(d)[Symbol(o)])
+  end
+
+  uc_d, uc_cd = [curry(get(d), old_d) for get in [dom, codom]]
+
+  α = Dict(map(collect(αcomps)) do (o, comps)
+    o => ACSetTransformation(ob_map(uc_d,   o),
+                             ob_map(uc_cd, o); comps...)
+  end)
+
+
+  FinTransformationMap(α, uc_d, uc_cd)
+end
+
 # ACSet serialization
 #####################
 
diff --git a/src/categorical_algebra/Categories.jl b/src/categorical_algebra/Categories.jl
index ea6586502..82908fd01 100644
--- a/src/categorical_algebra/Categories.jl
+++ b/src/categorical_algebra/Categories.jl
@@ -100,6 +100,7 @@ categories, checking equality of morphisms may involve nontrivial reasoning.
 is_hom_equal(C::Cat, f, g) = is_hom_equal(f, g)
 is_hom_equal(f, g) = f == g
 
+
 # Instances
 #----------
 
diff --git a/src/categorical_algebra/Diagrams.jl b/src/categorical_algebra/Diagrams.jl
index 99ddf5da5..d12efb2b0 100644
--- a/src/categorical_algebra/Diagrams.jl
+++ b/src/categorical_algebra/Diagrams.jl
@@ -10,10 +10,12 @@ import ...Theories: dom, codom, id, compose, ⋅, ∘, munit
 using ...Theories: ThCategory, composeH
 import ..Categories: ob_map, hom_map, op, co
 using ..FinCats, ..FreeDiagrams
-using ..FinCats: mapvals
+using ..FinCats: mapvals, FinDomFunctorMap, FinTransformationMap
 import ..FinCats: force, collect_ob, collect_hom
 import ..Limits: limit, colimit, universal
-
+using ..CSets: curry, uncurry, ACSetTransformation
+import ..CSets: homomorphisms, is_natural
+using ...CSetDataStructures: AnonACSetType
 # Data types
 ############
 
@@ -159,6 +161,11 @@ function Base.show(io::IO, f::DiagramHom{T}) where T
   print(io, ")")
 end
 
+function is_natural(f::DiagramHom)
+  all(is_functorial,[shape_map(f), diagram.([dom(f), codom(f)])...]
+     ) && is_natural(diagram_map(f))
+end
+
 # Categories of diagrams
 ########################
 
@@ -225,6 +232,58 @@ function compose(f::DiagramHom{T}, F::Functor; kw...) where T
                 compose(f.precomposed_diagram, F; kw...))
 end
 
+const ACSetCat{S} = TypeCat{S, ACSetTransformation}
+
+"""
+Search for natural transformations between Functors F,G: C->(D-Set) by
+translating the problem into natural transformations in CxD->Set.
+"""
+function homomorphisms(x::FinDomFunctorMap{D,ACSetCat{S}},
+                       y::FinDomFunctorMap{D,ACSetCat{S}}; kwargs...) where {D,S}
+    cx, cy = uncurry.([x,y])
+    dt = AnonACSetType(presentation(dom(cx)))
+    [curry(FinTransformationMap(h), x) for h in homomorphisms(dt(cx), dt(cy); kwargs...)]
+end
+
+function homomorphism(x::FinDomFunctorMap{D,ACSetCat{S}},
+                      y::FinDomFunctorMap{D,ACSetCat{S}}; kwargs...) where {D,S}
+  cx, cy = uncurry.([x,y])
+  dt = AnonACSetType(presentation(dom(cx)))
+  h = homomorphism(dt(cx), dt(cy); kwargs...)
+  return isnothing(h) ? nothing : curry(FinTransformationMap(h), x)
+end
+
+
+"""
+Diagram morphisms by first finding shape maps and then, for each,
+finding a diagram map.
+"""
+function homomorphisms(X::Diagram{T}, Y::Diagram{T}; n_max::Int=3,
+                       monic_obs=false, epic_obs=false,
+                       init_obs=nothing,init_homs=nothing,
+                       hom_lens=nothing, diag_kws=(;)) where T
+  fs = homomorphisms(shape(X),shape(Y); init_obs=init_obs, monic_obs=monic_obs,
+                     epic_obs=epic_obs, init_homs=init_homs, n_max=n_max,
+                     hom_lens=hom_lens)
+
+  res = []
+  for shapemap in fs
+    if T == id
+      hargs = [diagram(X), shapemap ⋅ diagram(Y)]
+    elseif T == op
+      hargs = [shapemap ⋅ diagram(Y), diagram(X)]
+    else
+      error("$T not supported")
+    end
+    tt = typeof(hargs[1]).parameters[2]
+    natural_transformations = homomorphisms(hargs...; diag_kws...)
+    for nt in natural_transformations
+      push!(res, DiagramHom{T}(shapemap, nt.components, X, Y))
+    end
+  end
+  return res
+end
+
 # Limits and colimits
 #####################
 
diff --git a/src/categorical_algebra/FinCats.jl b/src/categorical_algebra/FinCats.jl
index 441b486c1..b6ab3b134 100644
--- a/src/categorical_algebra/FinCats.jl
+++ b/src/categorical_algebra/FinCats.jl
@@ -14,7 +14,7 @@ export FinCat, FinCatGraph, Path, ob_generator, hom_generator,
   ob_generator_name, hom_generator_name, ob_generators, hom_generators,
   equations, is_discrete, is_free, graph, edges, src, tgt, presentation,
   FinFunctor, FinDomFunctor, is_functorial, collect_ob, collect_hom, force,
-  FinTransformation, components, is_natural, is_initial
+  FinTransformation, components, is_natural, is_initial, homomorphisms
 
 using StructEquality
 using Reexport
@@ -23,12 +23,12 @@ using DataStructures: IntDisjointSets, in_same_set, num_groups
 
 @reexport using ..Categories
 using ...GAT, ...Present, ...Syntax
-import ...Present: equations
+import ...Present: equations, generators
 using ...Theories: ThCategory, ThSchema, ObExpr, HomExpr, AttrExpr, AttrTypeExpr
 import ...Theories: dom, codom, id, compose, ⋅, ∘
 using ...CSetDataStructures, ...Graphs
 import ...Graphs: edges, src, tgt, enumerate_paths
-import ..Categories: CatSize, ob, hom, ob_map, hom_map, component, op
+import ..Categories: CatSize, ob, hom, ob_map, hom_map, component, op, is_hom_equal
 
 # Categories
 ############
@@ -259,6 +259,7 @@ See (Spivak, 2014, *Category theory for the sciences*, §4.5).
 end
 
 equations(C::FinCatGraphEq) = C.equations
+equations(C::FinCatGraph) = Pair[]
 
 function FinCatGraph(g::HasGraph, eqs::AbstractVector)
   eqs = map(eqs) do eq
@@ -270,6 +271,10 @@ function FinCatGraph(g::HasGraph, eqs::AbstractVector)
   FinCatGraphEq(g, eqs)
 end
 
+
+FinCatGraph(p::FinCatGraph) = p
+
+
 # Symbolic categories
 #####################
 
@@ -345,6 +350,31 @@ function Base.show(io::IO, C::FinCatPresentation)
   print(io, ")")
 end
 
+function FinCatGraph(p::FinCatPresentation)
+  pres = p.presentation
+  eqs = map(equations(pres)) do e12
+    e1,e2 = map(e12) do e
+      s, t = [findfirst(==(st), ob_generators(p)) for st in e.type_args]
+      et = only(typeof(e).parameters)
+      if et == :compose
+        return Path(Int[findfirst(==(x), hom_generators(p)) for x in e.args], s, t)
+      elseif et == :id
+        i = findfirst(==(e), ob_generators(p))
+        return Path(Int[], s, t)
+      else
+        Path(Int[findfirst(==(et), hom_generators(p))], s, t)
+      end
+    end
+    return e1 => e2
+  end
+
+  g = Graph(length(generators(pres, :Ob)))
+  hs = generators(pres, :Hom)
+  add_edges!(g, map(f -> generator_index(pres, first(gat_type_args(f))), hs),
+                map(f -> generator_index(pres, last(gat_type_args(f))), hs))
+  return isempty(eqs) ? FinCatGraph(g) : FinCatGraphEq(g,eqs)
+end
+
 # Functors
 ##########
 
@@ -639,6 +669,175 @@ function Base.show(io::IO, α::FinTransformationMap)
   print(io, ")")
 end
 
+
+
+"""
+Internal state for backtracking search for FinTransformations between
+FinFunctors between FinCatGraphs.
+"""
+struct BacktrackingState
+  """ The current assignment, a partially-defined homomorphism of ACSets. """
+  assignment::Vector{Union{Nothing,Path}}
+  """ Depth in search tree at which assignments were made. """
+  assignment_depth::Vector{Int}
+  """ Inverse assignment for monic components or if finding a monomorphism. """
+  inv_assignment::Union{Nothing,Vector{Int}}
+  """ Domain FinCat: the "variables" in the CSP. """
+  dom::FinFunctor
+  """ Codomain FinCat: the "values" in the CSP. """
+  codom::FinFunctor
+  """ Domain FinCat: the "variables" in the CSP. """
+  src::FinCatGraph
+  """ Codomain FinCat: the "values" in the CSP. """
+  tgt::FinCatGraph
+  cd_pths::Dict{Tuple{Int64, Int64}, Vector{Path}}
+end
+
+"""
+Search for FinTransformations between FinFunctors F,G: C->D. This requires
+sending each object A ∈ Ob(C) to a path of morphism generators in the codomain
+category, αₐ, such that the naturality square below commutes for all morphisms
+f ∈ Hom(A,B):
+          αₐ
+    F(A) --> G(A)
+ F(f)|        | G(f)
+     v        v
+    F(β) --> G(β)
+          αᵦ
+
+We first compute Hom(A,B) for all pairs of objects in D. This could potentially
+be infinite, so we restrict the max # of cycles possible in any morphism in this
+set with keyword `n_max`. This Hom-Set is quotiented by the equations of D,
+using a naive algorithm by default.
+
+We simply assign elements to objects of C in order (e.g. look in Hom(F(A),G(A))
+for each A). Naively we only check the naturality condition at the very end,
+though this ought to be done whenever a new assignment is made (a set of
+morphisms in C whose naturality has yet to be checked should be maintained).
+
+Therefore this code currently loops through all roughly |Path(D)|^Ob(C)
+possibilities and ought to be fixed up before applied to any problem of
+nontrivial scale.
+"""
+function homomorphisms(X::FinFunctor, Y::FinFunctor; kw...)
+  results = FinTransformation[]
+  backtracking_search(X, Y; kw...) do α
+    if is_natural(α)
+      push!(results, deepcopy(α));
+    end
+    return false
+  end
+  results
+end
+
+function backtracking_search(f, X::FinFunctor, Y::FinFunctor;
+                             n_max=2, initial::AbstractDict=Dict())
+  D, CD = dom(X), codom(X) # assume Y is also D->CD
+  Y_pths = Dict(map(collect(enumerate_paths_cyclic(CD.graph; n_max=n_max))) do (k,v)
+    v_ = [Path(x, k[1], k[2]) for x in v]
+    k => unique(h->normalize(CD, h), v_)
+  end)
+
+  # Initialize state variables for search.
+  assignment = fill(nothing, length(ob_generators(D)))
+  assignment_depth = zeros(Int, length(ob_generators(D)))
+  inv_assignment = nothing # currently not used
+  state = BacktrackingState(assignment, assignment_depth, inv_assignment,
+                            X, Y, D, CD, Y_pths)
+
+  # Start the main recursion for backtracking search.
+  for (k,v) in collect(initial)
+    assign_elem!(state, 0, k, v) || return false
+  end
+
+  backtracking_search(f, state, 1)
+end
+
+"""
+To reduce the branching factor, the element to branch on ought be the one with
+the smallest hom-set: Hom(F(x),G(x)). TODO.
+"""
+function backtracking_search(f, state::BacktrackingState, depth::Int) where {S}
+  # Choose the next unassigned element.
+  x = findfirst(isnothing, state.assignment)
+  if isnothing(x)
+    # No unassigned elements remain, so we have a complete assignment.
+    return f(FinTransformation(state.assignment, state.dom, state.codom))
+  end
+
+  # Attempt all assignments of the chosen element.
+  Y = state.codom
+  for y in state.cd_pths[(ob_map(state.dom,x), ob_map(state.codom,x))]
+    assign_elem!(state, depth, x, y) &&
+      backtracking_search(f, state, depth + 1) &&
+      return true
+    unassign_elem!(state, depth, x)
+  end
+  return false
+end
+
+
+""" Check whether element x can be assigned to y in current assignment.
+"""
+function can_assign_elem(state::BacktrackingState, depth, x, y)
+  ok = assign_elem!(state, depth, x, y)
+  unassign_elem!(state, depth, x)
+  return ok
+end
+
+""" Attempt to assign element (c,x) to (c,y) in the current assignment.
+
+Returns whether the assignment succeeded. Note that the backtracking state can
+be mutated even when the assignment fails.
+"""
+function assign_elem!(state::BacktrackingState, depth, x, y)
+  y′ = state.assignment[x]
+  y′ == y && return true  # If x is already assigned to y, return immediately.
+  isnothing(y′) || return false # Otherwise, x must be unassigned.
+
+  if !isnothing(state.inv_assignment) && state.inv_assignment[y] != 0
+    # Also, y must unassigned in the inverse assignment.
+    return false
+  end
+
+  # Make the assignment and recursively assign subparts.
+  state.assignment[x] = y
+  state.assignment_depth[x] = depth
+  if !isnothing(state.inv_assignment)
+    state.inv_assignment[y] = x
+  end
+  return true
+end
+
+""" Unassign the element (c,x) in the current assignment.
+"""
+function unassign_elem!(state::BacktrackingState, depth, x)
+  state.assignment[x] == 0 && return
+  assign_depth = state.assignment_depth[x]
+  @assert assign_depth <= depth
+  if assign_depth == depth
+    X = state.dom
+    if !isnothing(state.inv_assignment)
+      y = state.assignment[x]
+      state.inv_assignment[y] = 0
+    end
+    state.assignment[x] = nothing
+    state.assignment_depth[x] = 0
+  end
+end
+
+""" Get assignment pairs from partially specified component of C-set morphism.
+"""
+partial_assignments(x::AbstractDict) = pairs(x)
+partial_assignments(x::AbstractVector) =
+  ((i,y) for (i,y) in enumerate(x) if !isnothing(y) && y > 0)
+
+# FIXME: Should these accessors go elsewhere?
+in_hom(S, c) = [dom(S,f) => f for f in hom_generators(S) if codom(S,f) == c]
+out_hom(S, c) = [f => codom(S,f) for f in hom_generators(S) if dom(S,f) == c]
+
+
+
 # Initial functors
 ##################
 
@@ -746,4 +945,86 @@ function make_map(f, xs, ::Type{T}) where T
   end
 end
 
+# Word problems
+
+"""
+String rewriting (term rewriting with only unary operators) is an undecidable
+problem for the general case of arbitrary generators and equations.
+A MetaTheory.jl interop could give a more satisfying approximate solution.
+
+For now, we repeatedly apply equations as rewrite rules that make a
+term smaller or equal size.
+"""
+function is_hom_equal(C::FinCat, f, g)
+  is_hom_equal(normalize(C,f),  normalize(C,g))
+end
+
+
+term_size(t::Path) = length(edges(t))
+
+"""Used to orient equations into rewrite rules"""
+function term_size(t)
+  typ = only(typeof(t).parameters)
+  if typ == :id
+    0
+  elseif typ == :generator
+    1
+  elseif typ == :compose
+    length(t.args)
+  else error("is_hom_equal does not support terms of this type: $t")
+  end
+end
+"""Apply substitution to a list"""
+function sublist(xs, pat, rep)
+  pat != rep || error("Pattern and replacement must differ")
+  for i in 1:(length(xs)-(length(pat)-1))
+    if xs[i:i+length(pat)-1] == pat
+      return vcat(xs[1:i-1], rep, xs[i+length(pat):end]) => true
+    end
+  end
+  xs => false
+end
+
+function sub(t::Path,pat,rep)
+  new_edges, changed = sublist(edges(t), edges(pat), edges(rep))
+  return changed ? Path(new_edges, t.src, t.tgt) => true : t => false
+end
+
+"""Substitute pattern in a FreeSchema term"""
+function sub(t,pat,rep)
+  ttyp,pattyp, reptyp = [only(typeof(x).parameters) for x in [t,pat,rep]]
+  ttyp == :compose || return t => false
+  pattyp == :compose || error("pattern must be compose $pat")
+  if reptyp == :id
+    r = []
+  elseif reptyp == :generator
+    r = [rep]
+  else
+    r = rep.args
+  end
+  new_args, changed = sublist(t.args, pat.args, r)
+  if !changed
+    return t => false
+  elseif isempty(new_args)
+    return id(t.type_args[1]) => true
+  elseif length(new_args) == 1
+    return only(new_args) => true
+  else
+    return typeof(t)(new_args, t.type_args) => true
+  end
+end
+"""Use path equalities to make term as small as possible"""
+function normalize(C::FinCat, m)
+  changed = true
+  while changed
+    changed = false
+    for eqterms in equations(C)
+      r, l = sort(collect(eqterms); by=term_size)
+      m, was_changed = sub(m, l, r)
+      changed |= was_changed
+    end
+  end
+  m
+end
+
 end
diff --git a/src/categorical_algebra/Limits.jl b/src/categorical_algebra/Limits.jl
index 0b7cdc499..86c52e9b1 100644
--- a/src/categorical_algebra/Limits.jl
+++ b/src/categorical_algebra/Limits.jl
@@ -18,13 +18,13 @@ export AbstractLimit, AbstractColimit, Limit, Colimit,
 using StructEquality
 using StaticArrays: StaticVector, SVector
 
-using ...GAT, ...Theories
+using ...GAT, ...Theories, ...Present
 import ...Theories: ob, terminal, product, proj1, proj2, equalizer, incl,
   initial, coproduct, coproj1, coproj2, coequalizer, proj,
   delete, create, pair, copair, factorize, universal
 using ...CSetDataStructures, ..FinCats, ..FreeDiagrams
 import ..FreeDiagrams: apex, legs
-
+using ..FinCats: FinCatPresentation
 # Data types for limits
 #######################
 
@@ -606,4 +606,134 @@ function universal(colim::BipartiteColimit, cocone::Multicospan)
   universal(colim, cocone)
 end
 
+# Limits of FinCatPresentations
+###############################
+"""
+Decompose a product morphism into a composition of canonical generators, i.e.
+morphisms that have exactly one non-ID component
+Assumes we have a dictionary that takes us, e.g., from a tuple (f: A->B, idₓ) to
+the product morphism (f,idₓ): A×X⟶B×X
+"""
+function product_decompose(f::Vector, hdict::Dict)
+  curr = id.(dom.(f))
+  res, n = [], length(f)
+  for (i,m) in enumerate(f)
+    if (m isa HomExpr{:compose}) margs = m.args
+    elseif (m isa HomExpr{:generator}) margs = [m]
+    else margs = []
+    end
+    for marg in margs
+      args = tuple([j == i ? marg : curr[j] for j in 1:n]...)
+      push!(res, hdict[args])
+    end
+    curr[i] = id(codom(m))
+  end
+  return compose(res)
+end
+
+"""
+Product of finitely-presented categories has the cartesian product of ob
+generators as its objects. Its morphisms are generated by, for each morphism in
+the underlying categories, taking the product of that morphism with the
+identity morphism of all objects of all other categories. For example:
+ h                  f   g
+X->Y multiplied by A->B<-C is:
+                 (f,idY) (g, idY)
+                 YA -> YB  <-  YC
+         (h,idA) |     |(h,idB) | (h,idC)
+                 XA -> XB  <-  XC
+                 (f,idX) (g, idX)
+For any pair (e.g. f:A->B,h:X->Y), we get a naturality square
+                (f,idX)
+          A x X  ---->  B x X
+            |             |
+  (id(A),h) |             | (id(B),h)
+          A x Y  --->  B x Y
+              (f,id(Y))
+For any triple, we get a naturality cube (six naturality squares), and so on.
+TODO: figure out whether or not we also require 8 path equations to go from one
+corner of the cube to the other, or if that is derivable from the equalities on
+the faces of the cube. My intuition is: yes we need these additional equations.
+Currently we have no tests that fail due to this.
+"""
+function product(Xs::AbstractVector{<: FinCatPresentation}; kw...)
+  # Get cartesian product of obs and hosm
+  obs = collect(Iterators.product([ob_generators(x)  for x in Xs]...))[:]
+
+  homs = vcat(map(enumerate(Xs)) do (i,X)
+    vcat(map(hom_generators(X)) do h
+      p = Iterators.product([id.(ob_generators(Y)) for (j,Y) in enumerate(Xs) if j!=i]...)
+      map(collect.(collect(p)[:])) do hgens
+        tuple(insert!(Vector{Any}(hgens), i, h)...)
+      end
+    end...)
+  end...)
+
+  obdict = Dict([v=>k for (k,v) in enumerate(obs)])
+
+  # Create new presentation with tuple-looking names
+  p = Presentation(FreeSchema)
+
+  ogens = [Ob(FreeSchema, Symbol(o)) for o in obs]
+  map(ogens) do g add_generator!(p, g) end
+  hgens = Dict{Any,Any}(map(homs) do hs
+    src, tgt = map([dom, codom]) do get
+      ogens[obdict[tuple([get(X, h) for (h, X) in zip(hs,Xs)]...)]]
+    end
+    hs => add_generator!(p, Hom(Symbol(hs), src, tgt))
+  end)
+  for (k,v) in zip(obs, ogens)
+    hgens[tuple(id.(k)...)] = id(v)
+  end
+
+  # Add naturality squares
+  for i in 1:(length(Xs)-1)
+    hgi = hom_generators(Xs[i])
+    for j in i+1:length(Xs)
+      hgj = hom_generators(Xs[j])
+      fun(iarg, jarg, args) = map(1:length(Xs)) do k
+        if k==i iarg elseif k==j jarg else args[k] end end
+      for hij in Iterators.product(hgi, hgj)
+        hi, hj = hij
+        args = fun([nothing], [nothing], [id.(ob_generators(x))  for x in Xs])
+        for hs in Iterators.product(args...)
+          (di, cdi), (dj, cdj) = [id.([dom(x), codom(x)]) for x in hij]
+          xs = [hi=>dj, cdi=>hj, di=>hj, hi=>cdj] # naturality square
+          a1,a2,b1,b2 = [product_decompose(fun(x1,x2,hs), hgens) for (x1,x2) in xs]
+          add_equation!(p, compose(a1,a2), compose(b1,b2))
+        end
+      end
+    end
+  end
+
+  # Add equations from base categories
+  for (i, X) in enumerate(Xs)
+    for lr in equations(X)
+      is = [j == i ? [nothing] : id.(ob_generators(x))
+            for (j, x) in enumerate(Xs)]
+      for bkgrnd in Iterators.product(is...)
+        l_, r_ = map(lr) do t
+          comps = (t isa HomExpr{:id} || t isa HomExpr{:generator}) ? [t] : t.args
+          comps_ = map(comps) do comp
+            hgens[tuple([k==i ? comp : c for (k,c) in enumerate(bkgrnd)]...)]
+          end
+          compose(comps_)
+        end
+        add_equation!(p, l_, r_)
+      end
+    end
+  end
+
+  # Create projection maps
+  apx = FinCat(p)
+  ls = map(enumerate(Xs)) do (i, x)
+    os, hs = map([obs, homs]) do oldgens
+      Dict([Symbol(o) => o[i] for o in oldgens])
+    end
+    FinDomFunctor(os, hs, apx, x)
+  end
+  Limit(DiscreteDiagram(Xs), Multispan(ls))
+end
+
+
 end
diff --git a/src/graphs/GraphAlgorithms.jl b/src/graphs/GraphAlgorithms.jl
index d4dd59542..b3ffcc3ad 100644
--- a/src/graphs/GraphAlgorithms.jl
+++ b/src/graphs/GraphAlgorithms.jl
@@ -2,7 +2,7 @@
 """
 module GraphAlgorithms
 export connected_components, connected_component_projection, connected_component_projection_bfs,
-  topological_sort, transitive_reduction!, enumerate_paths
+  topological_sort, transitive_reduction!, enumerate_paths, enumerate_paths_cyclic
 
 using DataStructures: Stack, DefaultDict
 
@@ -147,4 +147,43 @@ function enumerate_paths(G::Graph;
   return res
 end
 
+"""
+Idea: have a DAG <: Graph type that stores the graph as well as
+an ordering. That way dispatch could use the correct function.
+
+Because cyclic graphs have an infinite number of paths, a cap on the the number of loops is required
+"""
+const Pth = Vector{Int}
+function enumerate_paths_cyclic(G::Graph; n_max=2)
+  ijs = collect(Iterators.product(vertices(G),vertices(G)))
+  es = Dict([(i,j)=>i==j ? [Int[]] : Pth[] for (i,j) in ijs])
+  n,done = 0,false
+
+  """False iff any vertex is visited more than n_max times"""
+  function count_cycles(p::Vector{Int})
+    cnt = zeros(Int, nv(G))
+    for e in p
+      cnt[src(G,e)] += 1
+      if cnt[src(G,e)] > n_max return false end
+    end
+    return true
+  end
+
+  while !done
+    done = true
+    n += 1 # we now add paths of length n
+    new_unseen = Dict([ij=>Int[] for ij in ijs])
+    for e in edges(G) # try to postcompose this edge w/ len n-1 paths
+      s, t = src(G,e), tgt(G,e)
+      for src_v in vertices(G)
+        for u in filter(u->length(u)==n-1 && count_cycles(u),es[(src_v, s)])
+          push!(es[(src_v, t)], [u; [e]])
+          done = false
+        end
+      end
+    end
+  end
+  return es
 end
+
+end # module
diff --git a/test/categorical_algebra/CSets.jl b/test/categorical_algebra/CSets.jl
index 51086a61d..43641d717 100644
--- a/test/categorical_algebra/CSets.jl
+++ b/test/categorical_algebra/CSets.jl
@@ -5,6 +5,7 @@ using JSON
 import JSONSchema
 
 using Catlab, Catlab.Theories, Catlab.Graphs, Catlab.CategoricalAlgebra
+using Catlab.CategoricalAlgebra.CSetDataStructures: AnonACSetType
 
 @present SchDDS(FreeSchema) begin
   X::Ob
@@ -489,6 +490,38 @@ A = Subobject(S, X=[3,4,5])
 @test ¬Subobject(ι₂) |> force == Subobject(ι₁)
 @test ~A |> force == ⊤(S) |> force
 
+# Currying
+##########
+
+const Grph = ACSetCat{Graph}
+g1 = Graph(1)
+ar = @acset Graph begin V=2; E=2; src=[1,2]; tgt=[2,2] end
+t1 = apex(terminal(Graph))
+t1_ar = homomorphism(t1, ar)
+_, g1_arr2 = homomorphisms(g1, ar)
+
+@present CSpanPres_(FreeSchema) begin
+  (C1, C2, C3)::Ob; c1::Hom(C1, C2); c2::Hom(C3,C2)
+end
+CSpan = FinCat(CSpanPres_)
+
+# Example FinFunctor into Grph
+CG_t1ar = FinDomFunctor(Dict(:C1=>t1,:C2=>ar,:C3=>g1),
+                        Dict(:c1=>t1_ar,:c2=>g1_arr2),
+                        CSpan, Grph());
+
+# FinFunctor to Grph -> FinFunctor to Set
+cspan_graph_ex = uncurry(CG_t1ar);
+# (reversible)
+@test Catlab.CategoricalAlgebra.curry(cspan_graph_ex, CG_t1ar) == CG_t1ar
+
+# FinFunctor to Set -> C-Set
+cg = AnonACSetType(presentation(dom(cspan_graph_ex)))
+cg_cset = cg(cspan_graph_ex)
+
+# C-Set -> FinFunctor to Set -> FinFunctor to Grph
+@test Catlab.CategoricalAlgebra.curry(FinDomFunctor(cg_cset), CG_t1ar) == CG_t1ar
+
 # Acset serialization
 #####################
 
diff --git a/test/categorical_algebra/Diagrams.jl b/test/categorical_algebra/Diagrams.jl
index 2af578c94..b674eb13e 100644
--- a/test/categorical_algebra/Diagrams.jl
+++ b/test/categorical_algebra/Diagrams.jl
@@ -1,7 +1,8 @@
 module TestDiagrams
+
 using Test
 
-using Catlab.Theories, Catlab.Graphs, Catlab.CategoricalAlgebra
+using Catlab.Theories, Catlab.Graphs, Catlab.CategoricalAlgebra, Catlab.Present
 
 const SchSGraph = SchSymmetricGraph
 
@@ -9,13 +10,10 @@ const SchSGraph = SchSymmetricGraph
 ##########
 
 # Diagram for paths of length 2.
-C = FinCat(@acset Graph begin
-  V = 3
-  E = 2
-  src = [1,2]
-  tgt = [3,3]
-end)
+C = FinCat(@acset Graph begin V = 3; E = 2; src = [1,2]; tgt = [3,3] end)
 D = FinDomFunctor([:E,:E,:V], [:tgt,:src], C, FinCat(SchSGraph))
+
+
 d = Diagram(D)
 @test shape(d) == C
 @test ob_map(d, 3) == SchSGraph[:V]
@@ -66,6 +64,93 @@ d = dom(f)
 @test codom(op(f)) == op(dom(f))
 @test op(g)⋅op(f) == op(f⋅g)
 
+# Morphism search
+##################
+
+# Diagrams in a FinCat
+
+"""
+Shapes:       •         •
+              ↓    ==>  ↓
+              •         •
+
+living        1 --> 2
+over          |     |
+square:       v     v
+              3 --> 4
+A diagram morphism in the forward direction is equivalent to a
+natural transformation. When the shape map is identity, this is
+the calculation of FinTransformations into Squarish, as found in the FinCats test suite. However, there are also shape maps which
+send  both objects to the first or the second object. Only in the
+second case does there exist a diagram morphism.
+
+"""
+Sq = @acset(Graph, begin V=4;E=4;src=[1,1,2,3]; tgt=[2,3,4,4] end)
+Squarish = FinCatGraph(Sq, [[[1,3],[2,4]]])
+
+Arr = FinCat(path_graph(Graph, 2))
+XF = FinFunctor((V=[1,3], E=[[2]]), Arr, Squarish)
+XG = FinFunctor((V=[2,4], E=[[3]]), Arr, Squarish)
+@test all(is_functorial,[XF,XG])
+F, G = Diagram.([XF,XG])
+hs = homomorphisms(F,G)
+@test length(hs) == 2
+@test length(homomorphisms(F,G; monic_obs=true))==1
+
+# There are no morphisms from RHS to LHS, so there are no op functors
+@test isempty(homomorphisms((Diagram{op}.([XF,XG]))...))
+
+# But if we look for op morphisms from the RHS *to* the LHS, we get 2
+# one of which has id shape map and the other which send both objs to 1.
+hs = homomorphisms((Diagram{op}.([XG,XF]))...)
+@test length(hs) == 2
+@test length(homomorphisms((Diagram{op}.([XG,XF]))...; monic_obs=true))==1
+
+# Diagrams in C-Set
+
+# dots shapes are just C-Sets
+@present SchDot(FreeSchema) begin X::Ob end
+toDot(x::StructACSet) = FinDomFunctor(Dict(:X=>x), nothing, FinCat(SchDot), ACSetCat{typeof(x)}())
+ar = path_graph(Graph, 2)
+g1 = toDot(ar)
+g2 = toDot(Sq)
+@test length(homomorphisms(g1,g2)) == length(homomorphisms(ar,Sq))
+
+# more complicated shapes
+
+const Grph = ACSetCat{Graph}
+@present SchArr(FreeSchema) begin (A,B)::Ob; f::Hom(A,B) end
+Arr = FinCat(SchArr)
+A,B = SchArr.generators[:Ob]
+
+toArr(f::ACSetTransformation) = FinDomFunctor(
+  Dict(:A=>dom(f), :B=>codom(f)),
+  Dict(:f=>f),
+  Arr, ACSetCat{typeof(dom(f))}())
+f1 = toArr(homomorphism(Graph(1), Graph(2)))
+f2 = toArr(homomorphism(Graph(2), Graph(1)))
+
+res = homomorphisms(f1,f2)
+@test length(res) == 2
+
+trm = terminal(Graph) |> apex
+f1 = toArr(homomorphism(Graph(2), ar; monic=true))
+f2 = toArr(homomorphism(Graph(2), trm ⊕ trm; monic=true))
+res = homomorphisms(f1,f2)
+@test length(res) == 2
+
+twoC, fourC =  [cycle_graph(Graph, i) for i in [2,4]]
+otherG = @acset Graph begin V=2; E=3; src=[1,1,2]; tgt=[1,2,1] end
+xF = FinDomFunctor(Dict(:A=>ar,:B=>fourC), Dict(:f=>homomorphism(ar,fourC)),
+                   Arr, Grph())
+xG = FinDomFunctor(Dict(:A=>twoC,:B=>otherG),
+                   Dict(:f=>homomorphism(twoC,otherG;monic=true)),
+                   Arr, Grph())
+@test all(is_functorial, [xF,xG])
+F,G=Diagram.([xF,xG]);
+res = homomorphisms(F,G; diag_kws=(monic=Symbol.(["(A, V)"]),));
+
+
 # Monads of diagrams
 ####################
 
diff --git a/test/categorical_algebra/FinCats.jl b/test/categorical_algebra/FinCats.jl
index fe44e1e66..d1b25f5b6 100644
--- a/test/categorical_algebra/FinCats.jl
+++ b/test/categorical_algebra/FinCats.jl
@@ -1,6 +1,6 @@
 module TestFinCats
-using Test
 
+using Test
 using Catlab, Catlab.Theories, Catlab.CategoricalAlgebra, Catlab.Graphs
 
 # Categories on graphs
@@ -18,6 +18,12 @@ C = FinCat(g)
 @test hom_generators(C) == 1:3
 @test startswith(sprint(show, C), "FinCat($(Graph)")
 
+# all gens to id (src or tgt) or 3*3*3 possibilities for V map = [1,2]
+@test length(homomorphisms(C,C)) == 29
+@test length(homomorphisms(C,C; init_obs=[1,2])) == 27
+@test length(homomorphisms(C,C; init_homs=[[1],[1],nothing])) == 3
+@test length(homomorphisms(C,C; hom_lens=[0, nothing, nothing])) == 2
+
 C_op = op(C)
 @test C_op isa FinCat
 @test (ob(C_op, 1), ob_generator(C_op, 1)) == (1, 1)
@@ -36,6 +42,9 @@ f = id(D, 2)
 g = compose(D, 1, 3)
 @test edges(g) == [1,3]
 
+homs = homomorphisms(D,D)
+@test all(is_functorial, homs)
+
 D_op = op(D)
 @test dom(D_op, 1) == 2
 @test codom(D_op, 1) == 1
@@ -45,12 +54,16 @@ D_op = op(D)
 # Functors between free categories.
 C = FinCat(parallel_arrows(Graph, 2))
 F = FinFunctor((V=[1,4], E=[[1,3], [2,4]]), C, D)
+@test length(homomorphisms(F,F))==1
 @test dom(F) == C
 @test codom(F) == D
 @test is_functorial(F)
 @test Ob(F) == FinFunction([1,4], FinSet(4))
 @test startswith(sprint(show, F), "FinFunctor($([1,4]),")
 
+homs = homomorphisms(C,D; monic_obs=true);
+@test F ∈ homs
+
 @test ob_map(F, 2) == 4
 @test hom_map(F, 1) == Path(h, [1,3])
 @test collect_ob(F) == [1,4]
@@ -136,6 +149,17 @@ end
 @test !is_free(Δ¹)
 @test startswith(sprint(show, Δ¹), "FinCat(")
 
+SG = FinCat(SchGraph)
+# Because of the equations, we must map V to V and E to E. There are four ways
+# to map the δ's onto each other.
+@test length(homomorphisms(Δ¹,Δ¹; n_max=0, monic_obs=true))==4
+
+
+SG = FinCat(SchGraph)
+# four ways to map src+tgt onto themselves
+@test length(homomorphisms(SG,SG; n_max=2, monic_obs=true))==4
+
+
 # Graph as set-valued functor on a free category.
 F = FinDomFunctor(path_graph(Graph, 3))
 C = dom(F)
@@ -204,39 +228,56 @@ G = FinDomFunctor(g)
 @test ob_map(G, :Weight) == TypeSet(Float64)
 @test hom_map(G, :weight) == FinDomFunction([0.5, 1.5])
 
+# FinTransformation search
+##########################
+"""
+Only one natural transformation between the functor which sends the arrow to
+the left hand side and another which sends it to the right hand side of a
+commutative square. If we add another edge on top yet don't make the composite
+path (from 1->4) commute with the lower route, then the naturality condition
+will forbid returning a result that uses this other arrow.
+
+         α_src
+        1  ⇉  2
+F(arr)  |     | G(arr)
+        v     v
+        3 --> 4
+         α_tgt
+"""
+Squarish = FinCatGraph(@acset(Graph, begin
+  V=4; E=5; src=[1,1,1,2,3]; tgt=[2,2,3,4,4]
+end), [[[1,4],[3,5]]])
+Arr = FinCat(path_graph(Graph, 2))
+F = FinFunctor((V=[1,3], E=[[3]]), Arr, Squarish)
+G = FinFunctor((V=[2,4], E=[[4]]), Arr, Squarish)
+FGs = homomorphisms(F,G)
+@test only(FGs) == FinTransformation([[1], [5]], F, G)
+@test FGs == homomorphisms(F,G, initial=Dict(2=>Path([5],3,4)))
+
+# try to force it to use the non-commuting arrow
+@test isempty(homomorphisms(F,G; initial=Dict(1=>Path([2],1,2))))
+
 # Initial functors
 ##################
 
 # Commutative square diagram: with 1→2→4 and 1→3→4
 S = FinCat(@acset Graph begin
-  V = 4
-  E = 4
-  src = [1,1,2,3]
-  tgt = [2,3,4,4]
+  V = 4; E = 4; src = [1,1,2,3]; tgt = [2,3,4,4]
 end)
 
 # Equalizer diagram: 1→2⇉3
 T = FinCat(@acset Graph begin
-  V = 3
-  E = 3
-  src = [1,2,2]
-  tgt = [2,3,3]
+  V = 3; E = 3; src = [1,2,2]; tgt = [2,3,3]
 end)
 
 # Extra bit added to beginning equalizer diagram: 4→1→2⇉3
 T2 = FinCat(@acset Graph begin
-  V = 4
-  E = 4
-  src = [1,2,2,4]
-  tgt = [2,3,3,1]
+  V = 4; E = 4; src = [1,2,2,4]; tgt = [2,3,3,1]
 end)
 
 # Extra bit added to end of equalizer diagram: 1→2⇉3→4
 T3 = FinCat(@acset Graph begin
-  V = 4
-  E = 4
-  src = [1,2,2,3]
-  tgt = [2,3,3,4]
+  V = 4; E = 4; src = [1,2,2,3]; tgt = [2,3,3,4]
 end)
 
 # Opposite square corners folded on top of each other
@@ -254,11 +295,11 @@ F4 = FinFunctor([1,2,3], [1,2,3], T, T2)
 # Same as F1, but there is an additional piece of data in codomain, ignored
 F5 = FinFunctor([1,2,3], [1,2,3], T, T2)
 
-@test all(is_functorial.([F1,F2,F3,F4]))
+@test all(is_functorial.([F1,F2,F3,F4,F5]))
 @test is_initial(F1)
 @test !is_initial(F2)
 @test !is_initial(F3)
 @test !is_initial(F4)
 @test !is_initial(F5)
 
-end
+end # module