diff --git a/benchmark/Graphs.jl b/benchmark/Graphs.jl
index 1ddedc166..826d926d5 100644
--- a/benchmark/Graphs.jl
+++ b/benchmark/Graphs.jl
@@ -11,6 +11,33 @@ using Catlab.Graphs.BasicGraphs: TheoryGraph
 using Catlab.WiringDiagrams: query
 using Catlab.Programs: @relation
 
+testdatadir = joinpath(dirname(@__FILE__), "..", "test", "testdata")
+
+# Example Graphs
+#
+################
+
+# Stolen from the Lightgraphs benchmark suite
+
+dg1fn = joinpath(testdatadir, "graph-50-500.jgz")
+
+LG_GRAPHS = Dict{String,LG.DiGraph}(
+    "complete100"   => LG.complete_digraph(100),
+    # "5000-50000"    => LG.loadgraph(dg1fn)["graph-5000-50000"],
+    "path500"       => LG.path_digraph(500)
+)
+
+GRAPHS = Dict(k => from_lightgraph(g) for (k,g) in LG_GRAPHS)
+
+LG_SYMGRAPHS = Dict{String,LG.Graph}(
+    "complete100"   => LG.complete_graph(100),
+    "tutte"         => LG.smallgraph(:tutte),
+    "path500"       => LG.path_graph(500),
+    # "5000-49947"    => LG.SimpleGraph(DIGRAPHS["5000-50000"])
+)
+
+SYMGRAPHS = Dict(k => from_lightgraph(g) for (k,g) in LG_SYMGRAPHS)
+
 # Helpers
 #########
 
@@ -49,7 +76,7 @@ end
 function bench_iter_neighbors(g)
   count = 0
   for v in vertices(g)
-    count += length(neighbors(g, v))
+    count += length(neighbors(g,v))
   end
   count
 end
@@ -65,6 +92,8 @@ end
 function Graphs.connected_component_projection(g::LG.AbstractGraph)
   label = Vector{Int}(undef, LG.nv(g))
   LG.connected_components!(label, g)
+  normalized = searchsortedfirst.(Ref(unique!(sort(label))), label)
+  FinFunction(normalized)
 end
 
 abstract type FindTrianglesAlgorithm end
@@ -116,29 +145,43 @@ lgbench["has-edge"] = @benchmarkable bench_has_edge($lg)
 clbench["iter-neighbors"] = @benchmarkable bench_iter_neighbors($g)
 lgbench["iter-neighbors"] = @benchmarkable bench_iter_neighbors($lg)
 
+
+bench = SUITE["GraphConnComponents"] = BenchmarkGroup()
+clbench = bench["Catlab"] = BenchmarkGroup()
+lgbench = bench["LightGraphs"] = BenchmarkGroup()
+
 n₀ = 2000
 g₀ = path_graph(Graph, n₀)
 g = ob(coproduct(fill(g₀, 5)))
 lg = LG.DiGraph(g)
-clbench["path-graph-components"] =
-  @benchmarkable connected_component_projection($g)
-lgbench["path-graph-components"] =
+clbench["path-graph"] =
+  @benchmarkable connected_component_projection_bfs($g)
+lgbench["path-graph"] =
   @benchmarkable connected_component_projection($lg)
 
 g₀ = star_graph(Graph, n₀)
 g = ob(coproduct(fill(g₀, 5)))
 lg = LG.DiGraph(g)
-clbench["star-graph-components"] =
-  @benchmarkable connected_component_projection($g)
-lgbench["star-graph-components"] =
+clbench["star-graph"] =
+  @benchmarkable connected_component_projection_bfs($g)
+lgbench["star-graph"] =
   @benchmarkable connected_component_projection($lg)
 
+for gn in keys(GRAPHS)
+  clbench[gn] = @benchmarkable connected_component_projection_bfs($(GRAPHS[gn]))
+  lgbench[gn] = @benchmarkable connected_component_projection($(LG_GRAPHS[gn]))
+end
+
+bench = SUITE["GraphTriangles"] = BenchmarkGroup()
+clbench = bench["Catlab"] = BenchmarkGroup()
+lgbench = bench["LightGraphs"] = BenchmarkGroup()
+
 n = 100
 g = wheel_graph(Graph, n)
 lg = LG.DiGraph(g)
 clbench["wheel-graph-triangles-hom"] =
   @benchmarkable ntriangles($g, TriangleBacktrackingSearch())
-clbench["wheel-graph-triangles-query"] =
+clbench["wheel-graph-triangles"] =
   @benchmarkable ntriangles($g, TriangleQuery())
 
 # Symmetric graphs
@@ -170,6 +213,10 @@ lgbench["has-edge"] = @benchmarkable bench_has_edge($lg)
 clbench["iter-neighbors"] = @benchmarkable bench_iter_neighbors($g)
 lgbench["iter-neighbors"] = @benchmarkable bench_iter_neighbors($lg)
 
+bench = SUITE["SymmetricGraphConnComponent"] = BenchmarkGroup()
+clbench = bench["Catlab"] = BenchmarkGroup()
+lgbench = bench["LightGraphs"] = BenchmarkGroup()
+
 n₀ = 2000
 g₀ = path_graph(SymmetricGraph, n₀)
 g = ob(coproduct(fill(g₀, 5)))
@@ -187,12 +234,22 @@ clbench["star-graph-components"] =
 lgbench["star-graph-components"] =
   @benchmarkable connected_component_projection($lg)
 
+for gn in keys(SYMGRAPHS)
+  clbench[gn] = @benchmarkable connected_component_projection_bfs($(SYMGRAPHS[gn]))
+  lgbench[gn] = @benchmarkable connected_component_projection($(LG_SYMGRAPHS[gn]))
+end
+
+bench = SUITE["SymmetricGraphTriangles"] = BenchmarkGroup()
+clbench = bench["Catlab"] = BenchmarkGroup()
+lgbench = bench["LightGraphs"] = BenchmarkGroup()
+
 n = 100
 g = wheel_graph(SymmetricGraph, n)
 lg = LG.Graph(g)
 clbench["wheel-graph-triangles-hom"] =
   @benchmarkable ntriangles($g, TriangleBacktrackingSearch())
-clbench["wheel-graph-triangles-query"] =
+# clbench["wheel-graph-triangles-query"] =
+clbench["wheel-graph-triangles"] =
   @benchmarkable ntriangles($g, TriangleQuery())
 lgbench["wheel-graph-triangles"] = @benchmarkable sum(LG.triangles($lg))
 
@@ -202,7 +259,7 @@ lgbench["wheel-graph-triangles"] = @benchmarkable sum(LG.triangles($lg))
 bench = SUITE["WeightedGraph"] = BenchmarkGroup()
 clbench = bench["Catlab"] = BenchmarkGroup()
 clvecbench = bench["Catlab-vectorized"] = BenchmarkGroup()
-lgbench = bench["LightGraphs"] = BenchmarkGroup()
+lgbench = bench["MetaGraphs"] = BenchmarkGroup()
 
 n = 10000
 g = path_graph(WeightedGraph{Float64}, n; E=(weight=range(0,1,length=n-1),))
@@ -305,3 +362,51 @@ lgbench["indexed-lookup"] = @benchmarkable begin
     @assert $mg["v$i", :label] == i
   end
 end
+
+# Random Graphs
+###############
+
+bench = SUITE["RandomGraph"] = BenchmarkGroup()
+clbench = bench["Catlab"] = BenchmarkGroup()
+lgbench = bench["LightGraphs"] = BenchmarkGroup()
+
+sizes = [10000]
+ps = [0.001]
+for size in sizes, p in ps
+  clbench["erdos_renyi-$size-$p"] =
+    @benchmarkable erdos_renyi($Graph, $size, $(p/2))
+  lgbench["erdos_renyi-$size-$p"] =
+    @benchmarkable LightGraphs.erdos_renyi($size, $p)
+end
+
+ks = [10]
+
+for size in sizes, k in ks
+  clbench["expected_degree_graph-$size-$k"] =
+    @benchmarkable expected_degree_graph($Graph, $([min(k,size-1) for _ in 1:size]))
+  lgbench["expected_degree_graph-$size-$k"] =
+    @benchmarkable LightGraphs.expected_degree_graph($([min(k,size-1) for _ in 1:size]))
+end
+
+for size in sizes, k in ks
+  clbench["watts_strogatz-$size-$k"] =
+    @benchmarkable watts_strogatz($Graph, $size, $(min(k,size-1)), 0.5)
+  lgbench["watts_strogatz-$size-$k"] =
+    @benchmarkable LightGraphs.watts_strogatz($size, $(min(k,size-1)), 0.5)
+end
+
+# Searching
+###########
+
+bench = SUITE["Searching"] = BenchmarkGroup()
+clbench = bench["Catlab"] = BenchmarkGroup()
+lgbench = bench["LightGraphs"] = BenchmarkGroup()
+
+for size in sizes, p in ps
+  local g = erdos_renyi(Graph, size, p)
+  local lg = LightGraphs.SimpleDiGraph(g)
+  clbench["bfs_erdos_renyi-$size-$p"] = @benchmarkable bfs_parents($g,1)
+  lgbench["bfs_erdos_renyi-$size-$p"] = @benchmarkable LightGraphs.bfs_parents($lg,1)
+  clbench["dfs_erdos_renyi-$size-$p"] = @benchmarkable dfs_parents($g,1)
+  lgbench["dfs_erdos_renyi-$size-$p"] = @benchmarkable LightGraphs.dfs_parents($lg,1)
+end
diff --git a/benchmark/make_plots.jl b/benchmark/make_plots.jl
new file mode 100644
index 000000000..e09878e80
--- /dev/null
+++ b/benchmark/make_plots.jl
@@ -0,0 +1,7 @@
+include("benchmarks.jl")
+
+include("plots.jl")
+
+results = run(SUITE)
+data = graphbench_data(results)
+plot_all_subcats(data)
diff --git a/benchmark/plots.jl b/benchmark/plots.jl
new file mode 100644
index 000000000..15fbf8fa2
--- /dev/null
+++ b/benchmark/plots.jl
@@ -0,0 +1,57 @@
+using Plots, StatsPlots
+using Plots.PlotMeasures
+using DataFrames, Query
+using BenchmarkTools
+
+function graphbench_data(suite)
+  data = DataFrame(subcat=String[],bench=String[],platform=String[],
+                   mt_normalized=Float64[],mediantime=Float64[])
+  graphbenches = suite["Graphs"]
+  noncatlab_times = Dict{Tuple{String,String},Float64}()
+  for (subcat,subsuite) in graphbenches
+    for (platform,results) in subsuite
+      for (bench,result) in results
+        if platform ∈ ["LightGraphs", "MetaGraphs"]
+            noncatlab_times[(subcat,bench)] = median(result).time
+        end
+        new_row = (subcat=subcat,
+                    bench=bench,
+                    platform=platform,
+                    mt_normalized=0.,
+                    mediantime=median(result).time)
+        push!(data, new_row)
+      end
+    end
+  end
+  for i in 1:length(data.subcat)
+    key = (data[i,:subcat], data[i,:bench])
+    if key ∈ keys(noncatlab_times)
+      data[i,:mt_normalized] = data[i,:mediantime] / noncatlab_times[key]
+    end
+  end
+  data
+end
+
+function subcat_data(dat,subcat)
+  dat |>
+    @filter(_.subcat==subcat) |>
+    @filter(_.platform == ["Catlab"]) |>
+    @orderby((_.bench,_.platform)) |>
+    @select(-:subcat) |>
+    DataFrame
+end
+
+function plot_subcat(dat,subcat,yscale=:linear)
+  subcat_data(dat,subcat) |>
+    @df groupedbar(:bench,:mt_normalized,group=:platform,
+                   xrotation=45,legend=:outerright,bar_width=0.5,
+                   yscale=yscale, yguide="Rel. time", bottom_margin=50px)
+end
+
+function plot_all_subcats(dat)
+  for subcat in unique(dat[!,:subcat])
+    yscale = subcat ∈ ["WeightedGraph", "LabeledGraph"] ? :log : :linear
+    fig = plot_subcat(dat,subcat,yscale)
+    savefig(fig, string("figures/",subcat,".pdf"))
+  end
+end
diff --git a/src/categorical_algebra/ACSetInterface.jl b/src/categorical_algebra/ACSetInterface.jl
index f1b279df8..480ff4dd2 100644
--- a/src/categorical_algebra/ACSetInterface.jl
+++ b/src/categorical_algebra/ACSetInterface.jl
@@ -62,12 +62,12 @@ of `[:src,:vattr]` as two independent columns does not even make sense, since
 they have different domains (belong to different tables).
 """
 function  subpart end
-@inline subpart(acs, part, name) = view_or_slice(subpart(acs, name), part)
+@inline Base.@propagate_inbounds subpart(acs, part, name) = view_or_slice(subpart(acs, name), part)
 
 function view_or_slice end
 @inline view_or_slice(x::AbstractVector, i::Union{Integer,StaticArray}) = x[i]
 @inline view_or_slice(x::AbstractVector, ::Colon) = x
-@inline view_or_slice(x::AbstractVector, i) = @view x[i]
+@inline Base.@propagate_inbounds view_or_slice(x::AbstractVector, i) = @view x[i]
 
 @inline subpart(acs, expr::GATExpr{:generator}) = subpart(acs, first(expr))
 @inline subpart(acs, expr::GATExpr{:id}) = parts(acs, first(dom(expr)))
@@ -120,16 +120,18 @@ end
 incident(acs, part, expr::GATExpr; kw...) =
   incident(acs, part, subpart_names(expr); kw...)
 
+@inline add_part!(acs, type; kw...) = add_part!(acs, type, (;kw...))
+
 """ Add part of given type to acset, optionally setting its subparts.
 
 Returns the ID of the added part.
 
 See also: [`add_parts!`](@ref).
 """
-@inline function add_part!(acs, type::Symbol, args...; kw...)
+@inline function add_part!(acs, type::Symbol, kw)
   part = only(add_parts!(acs,type,1))
   try
-    set_subparts!(acs, part, args...; kw...)
+    set_subparts!(acs, part, kw)
   catch e
     rem_part!(acs, type, part)
     rethrow(e)
@@ -145,10 +147,12 @@ See also: [`add_part!`](@ref).
 """
 function add_parts! end
 
-@inline function add_parts!(acs, type::Symbol, n::Int, args...; kw...)
+@inline add_parts!(acs, type::Symbol, n::Int; kw...) = add_parts!(acs, type, n, (;kw...))
+
+@inline function add_parts!(acs, type::Symbol, n::Int, kw)
   parts = add_parts!(acs, type, n)
   try
-    set_subparts!(acs, parts, args...; kw...)
+    set_subparts!(acs, parts, kw)
   catch e
     rem_parts!(acs, type, parts)
     rethrow(e)
@@ -183,10 +187,8 @@ Both single and vectorized assignment are supported.
 
 See also: [`set_subpart!`](@ref).
 """
-@inline function set_subparts!(acs, part, kw::NamedTuple)
-  for name in keys(kw)
-    set_subpart!(acs, part, name, kw[name])
-  end
+@inline @generated function set_subparts!(acs, part, kw::NamedTuple{keys}) where {keys}
+  Expr(:block,[:(set_subpart!(acs, part, $(Expr(:quote, name)), kw.$name)) for name in keys]...)
 end
 
 @inline set_subparts!(acs, part; kw...) = set_subparts!(acs, part, (;kw...))
diff --git a/src/categorical_algebra/CSetDataStructures.jl b/src/categorical_algebra/CSetDataStructures.jl
index 8b3b71471..87a80c291 100644
--- a/src/categorical_algebra/CSetDataStructures.jl
+++ b/src/categorical_algebra/CSetDataStructures.jl
@@ -245,7 +245,7 @@ end
 @inline Base.getindex(acs::StructACSet, args...) = ACSetInterface.subpart(acs, args...)
 
 @inline ACSetInterface.incident(acs::StructACSet, part, f::Symbol; copy::Bool=false) =
-  _incident(acs, part, Val{f}; copy=copy)
+  _incident(acs, part, Val{f}, Val{copy})
 
 broadcast_findall(xs, array::AbstractArray) =
   broadcast(x -> findall(y -> x == y, array), xs)
@@ -262,17 +262,18 @@ We keep the main body of the code generating out of the @generated function
 so that the code-generating function only needs to be compiled once.
 """
 function incident_body(s::SchemaDesc, idxed::AbstractDict{Symbol,Bool},
-                       unique_idxed::AbstractDict{Symbol,Bool}, f::Symbol)
+                       unique_idxed::AbstractDict{Symbol,Bool}, f::Symbol,
+                       copy::Bool)
   if f ∈ s.homs
     if idxed[f]
       quote
         indices = $(GlobalRef(ACSetInterface,:view_or_slice))(acs.hom_indices.$f, part)
-        copy ? Base.copy.(indices) : indices
+        $(copy ? :(Base.copy.(indices)) : :(indices))
       end
     elseif unique_idxed[f]
       quote
         indices = $(GlobalRef(ACSetInterface,:view_or_slice))(acs.hom_unique_indices.$f, part)
-        copy ? Base.copy.(indices) : indices
+        $(copy ? :(Base.copy.(indices)) : :(indices))
       end
     else
       :(broadcast_findall(part, acs.homs.$f))
@@ -281,7 +282,7 @@ function incident_body(s::SchemaDesc, idxed::AbstractDict{Symbol,Bool},
     if idxed[f]
       quote
         indices = get_attr_index(acs.attr_indices.$f, part)
-        copy ? Base.copy.(indices) : indices
+        $(copy ? :(Base.copy.(indices)) : :(indices))
       end
     elseif unique_idxed[f]
       quote
@@ -296,18 +297,41 @@ function incident_body(s::SchemaDesc, idxed::AbstractDict{Symbol,Bool},
 end
 
 @generated function _incident(acs::StructACSet{S,Ts,Idxed,UniqueIdxed},
-                              part, ::Type{Val{f}}; copy::Bool=false) where
-  {S,Ts,Idxed,UniqueIdxed,f}
-  incident_body(SchemaDesc(S),pairs(Idxed),pairs(UniqueIdxed),f)
+                              part, ::Type{Val{f}}, ::Type{Val{copy}}) where
+  {S,Ts,Idxed,UniqueIdxed,f,copy}
+  incident_body(SchemaDesc(S),pairs(Idxed),pairs(UniqueIdxed),f,copy)
 end
 
 # Mutators
 ##########
 
-@inline ACSetInterface.add_parts!(acs::StructACSet, ob::Symbol, n::Int) = _add_parts!(acs, Val{ob}, n)
+"""
+This is a specialized function to add parts to an ACSet and preallocate the indices of
+morphisms leading from those parts. This is useful if you want to reduce your
+total number of allocations when allocating an acset if you already know ahead of
+time a reasonable bound on the size of the preimages of the morphisms that you are using.
+
+For instance, if you are making a cyclic graph, then you know that the preimages of
+src and tgt will all be of size 1, and hence you can avoid allocating a zero-size
+array, and then again allocating a 1-size array and instead just allocate a
+1-size array off the bat.
+
+This function is currently exposed, but is not well-integrated with wrappers
+around add_parts; only use this if you really need it for performance and
+understand what you are doing. Additionally, the only guarantee w.r.t. to this
+is that it works the same semantically as `add_parts!`; it might make your code
+faster, but it also might not. Only use this if you have the benchmarks to back
+it up.
+"""
+@inline add_parts_with_indices!(acs::StructACSet, ob::Symbol, n::Int, index_sizes::NamedTuple) =
+  _add_parts!(acs, Val{ob}, n, index_sizes)
+
+@inline ACSetInterface.add_parts!(acs::StructACSet, ob::Symbol, n::Int) =
+  _add_parts!(acs, Val{ob}, n, (;))
 
 function add_parts_body(s::SchemaDesc, idxed::AbstractDict,
-                        unique_idxed::AbstractDict, ob::Symbol)
+                        unique_idxed::AbstractDict, ob::Symbol,
+                        index_sized_homs::Vector)
   code = quote
     m = acs.obs[$(ob_num(s, ob))]
     nparts = m + n
@@ -322,10 +346,16 @@ function add_parts_body(s::SchemaDesc, idxed::AbstractDict,
             end)
     end
     if s.codoms[f] == ob && idxed[f]
+      size = if f ∈ index_sized_homs
+        :(index_sizes[$(Expr(:quote, f))])
+      else
+        0
+      end
       push!(code.args, quote
             resize!(acs.hom_indices.$f, nparts)
             for i in newparts
-              acs.hom_indices.$f[i] = Int[]
+              acs.hom_indices.$f[i] = Array{Int}(undef, $size)
+              empty!(acs.hom_indices.$f[i])
             end
             end)
     elseif s.codoms[f] == ob && unique_idxed[f]
@@ -343,16 +373,17 @@ function add_parts_body(s::SchemaDesc, idxed::AbstractDict,
       push!(code.args,:(resize!(acs.attrs.$a, nparts)))
     end
   end
-  push!(code.args, :(newparts))
+  push!(code.args, :(return newparts))
   code
 end
 
 """ This generates the _add_parts! methods for a specific object of a `StructACSet`.
 """
 @generated function _add_parts!(acs::StructACSet{S,Ts,Idxed,UniqueIdxed},
-                                ::Type{Val{ob}}, n::Int) where
-  {S, Ts, Idxed, UniqueIdxed, ob}
-  add_parts_body(SchemaDesc(S),pairs(Idxed),pairs(UniqueIdxed),ob)
+                                ::Type{Val{ob}}, n::Int,
+                                index_sizes::NamedTuple{index_sized_homs}) where
+  {S, Ts, Idxed, UniqueIdxed, ob, index_sized_homs}
+  add_parts_body(SchemaDesc(S),pairs(Idxed),pairs(UniqueIdxed),ob,[index_sized_homs...])
 end
 
 @inline ACSetInterface.set_subpart!(acs::StructACSet, part::Int, f::Symbol, subpart) =
@@ -365,14 +396,15 @@ function set_subpart_body(s::SchemaDesc, idxed::AbstractDict{Symbol,Bool},
     if idxed[f]
       quote
         @assert 0 <= subpart <= acs.obs[$(ob_num(s, s.codoms[f]))]
-        old = acs.homs.$f[part]
-        acs.homs.$f[part] = subpart
+        @inbounds old = acs.homs.$f[part]
+        @inbounds acs.homs.$f[part] = subpart
         if old > 0
           @assert deletesorted!(acs.hom_indices.$f[old], part)
         end
         if subpart > 0
           insertsorted!(acs.hom_indices.$f[subpart], part)
         end
+        subpart
       end
     elseif unique_idxed[f]
       quote
@@ -384,11 +416,13 @@ function set_subpart_body(s::SchemaDesc, idxed::AbstractDict{Symbol,Bool},
         end
         acs.homs.$f[part] = subpart
         acs.hom_unique_indices.$f[subpart] = part
+        subpart
       end
     else
       quote
         @assert 0 <= subpart <= acs.obs[$(ob_num(s, s.codoms[f]))]
         acs.homs.$f[part] = subpart
+        subpart
       end
     end
   elseif f ∈ s.attrs
@@ -400,16 +434,18 @@ function set_subpart_body(s::SchemaDesc, idxed::AbstractDict{Symbol,Bool},
         end
         acs.attrs.$f[part] = subpart
         set_attr_index!(acs.attr_indices.$f, subpart, part)
+        subpart
       end
     elseif unique_idxed[f]
       quote
-        @assert subpart ∉ keys(acs.attr_unique_indices.$f) "subpart not unique"
+        @boundscheck @assert subpart ∉ keys(acs.attr_unique_indices.$f) "subpart not unique"
         if isassigned(acs.attrs.$f, part)
           old = acs.attrs.$f[part]
           delete!(acs.attr_unique_indices.$f, old)
         end
         acs.attrs.$f[part] = subpart
         acs.attr_unique_indices.$f[subpart] = part
+        subpart
       end
     else
       :(acs.attrs.$f[part] = subpart)
@@ -433,11 +469,6 @@ end
 @inline ACSetInterface.rem_part!(acs::StructACSet, type::Symbol, part::Int) =
   _rem_part!(acs, Val{type}, part)
 
-function getassigned(acs::StructACSet, arrows, i)
-  assigned_subparts = filter(f -> isassigned(subpart(acs,f),i), arrows)
-  Dict(f => subpart(acs,i,f) for f in assigned_subparts)
-end
-
 function rem_part_body(s::SchemaDesc, idxed, ob::Symbol)
   in_homs = filter(hom -> s.codoms[hom] == ob, s.homs)
   out_homs = filter(f -> s.doms[f] == ob, s.homs)
@@ -445,7 +476,7 @@ function rem_part_body(s::SchemaDesc, idxed, ob::Symbol)
   indexed_out_homs = filter(hom -> s.doms[hom] == ob && idxed[hom], s.homs)
   indexed_attrs = filter(attr -> s.doms[attr] == ob && idxed[attr], s.attrs)
   quote
-    last_part = acs.obs[$(ob_num(s, ob))]
+    last_part = @inbounds acs.obs[$(ob_num(s, ob))]
     @assert 1 <= part <= last_part
     # Unassign superparts of the part to be removed and also reassign superparts
     # of the last part to this part.
@@ -453,15 +484,27 @@ function rem_part_body(s::SchemaDesc, idxed, ob::Symbol)
       set_subpart!(acs, incident(acs, part, hom, copy=true), hom, 0)
       set_subpart!(acs, incident(acs, last_part, hom, copy=true), hom, part)
     end
-    last_row = getassigned(acs, $([out_homs;out_attrs]), last_part)
+
+    # This is a hack to avoid allocating a named tuple, because these parts
+    # are a union type, so there would be dynamic dispatch
+    $(Expr(:block, (map([out_homs; out_attrs]) do f
+         :($(Symbol("last_row_" * string(f))) =
+           if isassigned(subpart(acs, $(Expr(:quote,f))), last_part)
+             subpart(acs, last_part, $(Expr(:quote,f)))
+           else
+             nothing
+           end)
+       end)...))
 
     # Clear any morphism and data attribute indices for last part.
-    for hom in $(Tuple(indexed_out_homs))
-      set_subpart!(acs, last_part, hom, 0)
-    end
+    $(Expr(:block,
+           (map(indexed_out_homs) do hom
+              :(set_subpart!(acs, last_part, $(Expr(:quote, hom)), 0))
+            end)...))
+
     for attr in $(Tuple(indexed_attrs))
-      if haskey(last_row, attr)
-        unset_attr_index!(acs.attr_indices[attr], last_row[attr], last_part)
+      if isassigned(subpart(acs, attr), last_part)
+        unset_attr_index!(acs.attr_indices[attr], subpart(acs, last_part, attr), last_part)
       end
     end
 
@@ -472,9 +515,17 @@ function rem_part_body(s::SchemaDesc, idxed, ob::Symbol)
     for a in $(Tuple(out_attrs))
       resize!(acs.attrs[a], last_part - 1)
     end
-    acs.obs[$(ob_num(s, ob))] -= 1
+    @inbounds acs.obs[$(ob_num(s, ob))] -= 1
     if part < last_part
-      set_subparts!(acs, part, (;last_row...))
+      $(Expr(:block,
+             (map([out_homs; out_attrs]) do f
+                quote
+                  x = $(Symbol("last_row_" * string(f)))
+                  if !isnothing(x)
+                    set_subpart!(acs, part, $(Expr(:quote, f)), x)
+                  end
+                end
+             end)...))
     end
   end
 end
diff --git a/src/categorical_algebra/GraphCategories.jl b/src/categorical_algebra/GraphCategories.jl
index 936097af6..b6f04002b 100644
--- a/src/categorical_algebra/GraphCategories.jl
+++ b/src/categorical_algebra/GraphCategories.jl
@@ -2,13 +2,42 @@
 """
 module GraphCategories
 
+using DataStructures
+
 using ..FinSets, ...ACSetInterface, ..Limits
 using ...Graphs.BasicGraphs
-import ...Graphs.GraphAlgorithms: connected_component_projection
+import ...Graphs.GraphAlgorithms: connected_component_projection,
+  connected_component_projection_bfs
 
 function connected_component_projection(g::ACSet)::FinFunction
   proj(coequalizer(FinFunction(src(g), nv(g)),
                    FinFunction(tgt(g), nv(g))))
 end
 
+# This algorithm is linear in the number of vertices of g, so it should be
+# significantly faster than the previous one in some cases.
+function connected_component_projection_bfs(g::ACSet)
+  label = zeros(Int, nv(g))
+
+  q = Queue{Int}()
+  for v in 1:nv(g)
+    label[v] != 0 && continue
+    label[v] = v
+    empty!(q)
+    enqueue!(q, v)
+    while !isempty(q)
+      src = dequeue!(q)
+      for vertex in neighbors(g, src)
+        if label[vertex] == 0
+          enqueue!(q,vertex)
+          label[vertex] = v
+        end
+      end
+    end
+  end
+
+  normalized = searchsortedfirst.(Ref(unique!(sort(label))), label)
+  FinFunction(normalized)
+end
+
 end
diff --git a/src/graphs/BasicGraphs.jl b/src/graphs/BasicGraphs.jl
index b4285c98b..2719bdd11 100644
--- a/src/graphs/BasicGraphs.jl
+++ b/src/graphs/BasicGraphs.jl
@@ -8,17 +8,18 @@ departures due to differences between the data structures.
 """
 module BasicGraphs
 export HasVertices, HasGraph,
-  AbstractGraph, Graph, nv, ne, src, tgt, edges, vertices,
-  has_edge, has_vertex, add_edge!, add_edges!, add_vertex!, add_vertices!,
+  AbstractGraph, Graph, nv, ne, src, tgt, edges, inedges, outedges, vertices,
+  has_edge, has_vertex, add_edge!, add_edges!, add_vertex!, add_vertices!, add_vertices_with_indices!,
   rem_edge!, rem_edges!, rem_vertex!, rem_vertices!,
-  neighbors, inneighbors, outneighbors, all_neighbors, induced_subgraph,
+  neighbors, inneighbors, outneighbors, all_neighbors, degree, induced_subgraph,
   AbstractSymmetricGraph, SymmetricGraph, inv,
   AbstractReflexiveGraph, ReflexiveGraph, refl,
   AbstractSymmetricReflexiveGraph, SymmetricReflexiveGraph,
   AbstractHalfEdgeGraph, HalfEdgeGraph, vertex, half_edges,
   add_dangling_edge!, add_dangling_edges!,
   AbstractWeightedGraph, WeightedGraph, weight,
-  AbstractSymmetricWeightedGraph, SymmetricWeightedGraph
+  AbstractSymmetricWeightedGraph, SymmetricWeightedGraph,
+  from_lightgraph
 
 import Base: inv
 using Requires
@@ -101,6 +102,14 @@ edges(g::HasGraph) = parts(g, :E)
 edges(g::HasGraph, src::Int, tgt::Int) =
   (e for e in incident(g, src, :src) if subpart(g, e, :tgt) == tgt)
 
+""" Edges coming out of a vertex
+"""
+outedges(g::HasGraph, v) = incident(g, v, :src)
+
+""" Edges coming into a vertex
+"""
+inedges(g::HasGraph, v) = incident(g, v, :tgt)
+
 """ Whether the graph has the given vertex.
 """
 has_vertex(g::HasVertices, v) = has_part(g, :V, v)
@@ -108,8 +117,13 @@ has_vertex(g::HasVertices, v) = has_part(g, :V, v)
 """ Whether the graph has the given edge, or an edge between two vertices.
 """
 has_edge(g::HasGraph, e) = has_part(g, :E, e)
-has_edge(g::HasGraph, src::Int, tgt::Int) =
-  has_vertex(g, src) && tgt ∈ outneighbors(g, src)
+function has_edge(g::HasGraph, s::Int, t::Int)
+  (1 <= s <= nv(g)) || return false
+  for e in outedges(g,s)
+    (tgt(g,e) == t) && return true
+  end
+  false
+end
 
 """ Add a vertex to a graph.
 """
@@ -119,17 +133,24 @@ add_vertex!(g::HasVertices; kw...) = add_part!(g, :V; kw...)
 """
 add_vertices!(g::HasVertices, n::Int; kw...) = add_parts!(g, :V, n; kw...)
 
+""" Add vertices with preallocated src/tgt indexes
+"""
+function add_vertices_with_indices!(g::HasVertices, n::Int, k::Int; kw...)
+  CSetDataStructures.add_parts_with_indices!(g, :V, n, (src=k,tgt=k))
+  set_subparts!(g, :V; kw...)
+end
+
 """ Add an edge to a graph.
 """
 add_edge!(g::HasGraph, src::Int, tgt::Int; kw...) =
-  add_part!(g, :E; src=src, tgt=tgt, kw...)
+  add_part!(g, :E, (src=src, tgt=tgt, kw...))
 
 """ Add multiple edges to a graph.
 """
 function add_edges!(g::HasGraph, srcs::AbstractVector{Int},
                     tgts::AbstractVector{Int}; kw...)
   @assert (n = length(srcs)) == length(tgts)
-  add_parts!(g, :E, n; src=srcs, tgt=tgts, kw...)
+  add_parts!(g, :E, n, (src=srcs, tgt=tgts, kw...))
 end
 
 """ Remove a vertex from a graph.
@@ -173,21 +194,27 @@ distinction is moot.
 In the presence of multiple edges, neighboring vertices are given *with
 multiplicity*. To get the unique neighbors, call `unique(neighbors(g))`.
 """
-neighbors(g::AbstractGraph, v::Int) = outneighbors(g, v)
+@inline neighbors(g::AbstractGraph, v::Int) = outneighbors(g, v)
 
 """ In-neighbors of vertex in a graph.
 """
-inneighbors(g::AbstractGraph, v::Int) = subpart(g, incident(g, v, :tgt), :src)
+@inline inneighbors(g::AbstractGraph, v::Int) = @inbounds subpart(g, incident(g, v, :tgt), :src)
 
 """ Out-neighbors of vertex in a graph.
 """
-outneighbors(g::AbstractGraph, v::Int) = subpart(g, incident(g, v, :src), :tgt)
+@inline outneighbors(g::AbstractGraph, v::Int) = @inbounds subpart(g, incident(g, v, :src), :tgt)
 
 """ Union of in-neighbors and out-neighbors in a graph.
 """
 all_neighbors(g::AbstractGraph, v::Int) =
   Iterators.flatten((inneighbors(g, v), outneighbors(g, v)))
 
+""" Total degree of a vertex
+
+Equivalent to length(all_neighbors(g,v)) but faster
+"""
+degree(g,v) = length(incident(g,v,:tgt)) + length(incident(g,v,:src))
+
 """ Subgraph induced by a set of a vertices.
 
 The [induced subgraph](https://en.wikipedia.org/wiki/Induced_subgraph) consists
@@ -524,6 +551,23 @@ function __init__()
       lg
     end
 
+    function from_lightgraph(lg::SimpleDiGraph)
+      g = Graph(LightGraphs.nv(lg))
+      for e in LightGraphs.edges(lg)
+        add_edge!(g,LightGraphs.src(e),LightGraphs.dst(e))
+      end
+      g
+    end
+
+    function from_lightgraph(lg::SimpleGraph)
+      g = SymmetricGraph(LightGraphs.nv(lg))
+      for e in LightGraphs.edges(lg)
+        add_edge!(g,LightGraphs.src(e),LightGraphs.dst(e))
+      end
+      g
+    end
+
+
     function SimpleGraph(g::AbstractHalfEdgeGraph)
       lg = SimpleGraph(nv(g))
       for e in half_edges(g)
diff --git a/src/graphs/GraphAlgorithms.jl b/src/graphs/GraphAlgorithms.jl
index 75d66bbeb..98c4fe98c 100644
--- a/src/graphs/GraphAlgorithms.jl
+++ b/src/graphs/GraphAlgorithms.jl
@@ -1,8 +1,8 @@
 """ Algorithms on graphs based on C-sets.
 """
 module GraphAlgorithms
-export connected_components, connected_component_projection, topological_sort,
-  transitive_reduction!, enumerate_paths
+export connected_components, connected_component_projection, connected_component_projection_bfs,
+  topological_sort, transitive_reduction!, enumerate_paths
 
 using DataStructures: Stack, DefaultDict
 
@@ -31,6 +31,8 @@ end
 Returns a function in FinSet{Int} from the vertex set to the set of components.
 """
 function connected_component_projection end
+
+function connected_component_projection_bfs end
 # Implemented elsewhere, where coequalizers are available.
 
 # DAGs
diff --git a/src/graphs/GraphGenerators.jl b/src/graphs/GraphGenerators.jl
index 375ecf898..3b1462b84 100644
--- a/src/graphs/GraphGenerators.jl
+++ b/src/graphs/GraphGenerators.jl
@@ -1,9 +1,11 @@
 module GraphGenerators
 export path_graph, cycle_graph, complete_graph, star_graph, wheel_graph,
-  parallel_arrows
+  parallel_arrows, erdos_renyi, expected_degree_graph, watts_strogatz
 
 using ...CSetDataStructures, ..BasicGraphs
 using ...CSetDataStructures: hom
+using Random
+using Random: GLOBAL_RNG
 
 """ Path graph on ``n`` vertices.
 """
@@ -72,4 +74,246 @@ end
 # Should this be exported from `BasicGraphs`?
 @generated is_directed(::Type{T}) where {S, T<:StructACSet{S}} = :inv ∉ hom(S)
 
+getRNG(seed::Integer, rng::AbstractRNG) = seed >= 0 ? MersenneTwister(seed) : rng
+
+"""
+    randbn(n, p, seed=-1)
+
+Return a binomally-distribted random number with parameters `n` and `p` and optional `seed`.
+
+### Optional Arguments
+- `seed=-1`: set the RNG seed.
+- `rng`: set the RNG directly
+
+### References
+- "Non-Uniform Random Variate Generation," Luc Devroye, p. 522. Retrieved via http://www.eirene.de/Devroye.pdf.
+- http://stackoverflow.com/questions/23561551/a-efficient-binomial-random-number-generator-code-in-java
+- https://github.com/JuliaGraphs/LightGraphs.jl/blob/2a644c2b15b444e7f32f73021ec276aa9fc8ba30/src/SimpleGraphs/generators/randgraphs.jl#L90
+"""
+function randbn(n::Integer, p::Real; rng::AbstractRNG=GLOBAL_RNG)
+  log_q = log(1.0 - p)
+  x = 0
+  sum = 0.0
+  while true
+    sum += log(rand(rng)) / (n - x)
+    sum < log_q && break
+    x += 1
+  end
+  return x
+end
+
+"""
+    erdos_renyi(GraphType, n, p)
+
+Create an [Erdős–Rényi](http://en.wikipedia.org/wiki/Erdős–Rényi_model)
+random graph with `n` vertices. Edges are added between pairs of vertices with
+probability `p`.
+
+### Optional Arguments
+- `seed=-1`: set the RNG seed.
+- `rng`: set the RNG directly
+
+### References
+- https://github.com/JuliaGraphs/LightGraphs.jl/blob/2a644c2b15b444e7f32f73021ec276aa9fc8ba30/src/SimpleGraphs/generators/randgraphs.jl
+"""
+function erdos_renyi(::Type{T}, n::Int, p::Real; V=(;),
+                     seed::Int=-1, rng::AbstractRNG=GLOBAL_RNG) where T <: ACSet
+  rng = getRNG(seed,rng)
+  _erdos_renyi(T,n,p,V,rng)
+end
+
+function _erdos_renyi(::Type{T}, n::Int, p::Real, V, rng::AbstractRNG) where T <: ACSet
+  p >= 1 && return complete_graph(T, n)
+  maxe = n * (n-1)
+  m = randbn(maxe,p;rng=rng)
+  return _erdos_renyi(T,n,m,V,rng)
+end
+
+
+"""
+    erdos_renyi(GraphType, n, ne)
+
+Create an [Erdős–Rényi](http://en.wikipedia.org/wiki/Erdős–Rényi_model) random
+graph with `n` vertices and `ne` edges.
+
+### References
+- https://github.com/JuliaGraphs/LightGraphs.jl/blob/2a644c2b15b444e7f32f73021ec276aa9fc8ba30/src/SimpleGraphs/generators/randgraphs.jl
+"""
+function erdos_renyi(::Type{T}, n::Int, m::Int; V=(;),
+                     seed::Int=-1, rng::AbstractRNG=GLOBAL_RNG) where T <: ACSet
+  rng = getRNG(seed, rng)
+  _erdos_renyi(T,n,m,V,rng)
+end
+
+function _erdos_renyi(::Type{T}, n::Int, m::Int, V, rng::AbstractRNG) where T <: ACSet
+  maxe = n * (n-1)
+  maxe == m && return complete_graph(T, n)
+  @assert(m <= maxe, "Maximum number of edges for this graph is $maxe")
+  # In the case of a symmetric graph, the edges are double-counted
+  totale = is_directed(T) ? m : 2*m
+  k = Int(ceil(m/n))
+
+  g = T()
+  add_vertices_with_indices!(g, n, k)
+  set_subparts!(g,:V,V)
+
+  while ne(g) < totale
+    src = rand(rng, 1:n)
+    tgt = rand(rng, 1:n)
+    src != tgt && !has_edge(g,src,tgt) && add_edge!(g,src,tgt)
+  end
+
+  return g
+end
+
+
+"""
+    expected_degree_graph(GraphType, ω)
+
+Given a vector of expected degrees `ω` indexed by vertex, create a random undirected graph in which vertices `i` and `j` are
+connected with probability `ω[i]*ω[j]/sum(ω)`.
+
+### Optional Arguments
+- `seed=-1`: set the RNG seed.
+- `rng`: set the RNG directly
+
+### Implementation Notes
+The algorithm should work well for `maximum(ω) << sum(ω)`. As `maximum(ω)` approaches `sum(ω)`, some deviations
+from the expected values are likely.
+
+### References
+- Connected Components in Random Graphs with Given Expected Degree Sequences, Linyuan Lu and Fan Chung. [https://link.springer.com/article/10.1007%2FPL00012580](https://link.springer.com/article/10.1007%2FPL00012580)
+- Efficient Generation of Networks with Given Expected Degrees, Joel C. Miller and Aric Hagberg. [https://doi.org/10.1007/978-3-642-21286-4_10](https://doi.org/10.1007/978-3-642-21286-4_10)
+- https://github.com/JuliaGraphs/LightGraphs.jl/blob/2a644c2b15b444e7f32f73021ec276aa9fc8ba30/src/SimpleGraphs/generators/randgraphs.jl#L187
+"""
+function expected_degree_graph(::Type{T},ω::Vector{<:Real}; V=(;),
+                               seed::Int=-1, rng::AbstractRNG=GLOBAL_RNG) where T <: ACSet
+  rng = getRNG(seed, rng)
+  g = T()
+  add_vertices!(g,length(ω);V...)
+  expected_degree_graph!(g, ω, rng=rng)
+end
+
+function expected_degree_graph!(g::T, ω::Vector{U};
+                                seed::Int=-1, rng::AbstractRNG=GLOBAL_RNG) where {T <: ACSet, U <: Real}
+  rng = getRNG(seed, rng)
+  n = length(ω)
+  @assert all(zero(U) .<= ω .<= n - one(U)) "Elements of ω need to be at least 0 and at most n-1"
+
+  π = sortperm(ω, rev=true)
+
+  S = sum(ω)
+
+  for u = 1:(n - 1)
+    v = u + 1
+    p = min(ω[π[u]] * ω[π[v]] / S, one(U))
+    while v <= n && p > zero(p)
+      if p != one(U)
+        v += floor(Int, log(rand(rng)) / log(one(U) - p))
+      end
+      if v <= n
+        q = min(ω[π[u]] * ω[π[v]] / S, one(U))
+        if rand(rng) < q / p
+          add_edge!(g, π[u], π[v])
+        end
+        p = q
+        v += 1
+      end
+    end
+  end
+  return g
+end
+
+"""
+    watts_strogatz(n, k, β)
+
+Return a [Watts-Strogatz](https://en.wikipedia.org/wiki/Watts_and_Strogatz_model)
+small world random graph with `n` vertices, each with expected degree `k` (or `k
+- 1` if `k` is odd). Edges are randomized per the model based on probability `β`.
+
+The algorithm proceeds as follows. First, a perfect 1-lattice is constructed,
+where each vertex has exacly `div(k, 2)` neighbors on each side (i.e., `k` or
+`k - 1` in total). Then the following steps are repeated for a hop length `i` of
+`1` through `div(k, 2)`.
+
+1. Consider each vertex `s` in turn, along with the edge to its `i`th nearest
+   neighbor `t`, in a clockwise sense.
+
+2. Generate a uniformly random number `r`. If `r ≥ β`, then the edge `(s, t)` is
+   left unaltered. Otherwise, the edge is deleted and *rewired* so that `s` is
+   connected to some vertex `d`, chosen uniformly at random from the entire
+   graph, excluding `s` and its neighbors. (Note that `t` is a valid candidate.)
+
+For `β = 0`, the graph will remain a 1-lattice, and for `β = 1`, all edges will
+be rewired randomly.
+
+### Optional Arguments
+- `is_directed=false`: if true, return a directed graph.
+- `seed=-1`: set the RNG seed.
+
+### References
+- Collective dynamics of ‘small-world’ networks, Duncan J. Watts, Steven H. Strogatz. [https://doi.org/10.1038/30918](https://doi.org/10.1038/30918)
+- Small Worlds, Duncan J. watts. [https://en.wikipedia.org/wiki/Special:BookSources?isbn=978-0691005416](https://en.wikipedia.org/wiki/Special:BookSources?isbn=978-0691005416)
+- https://github.com/JuliaGraphs/LightGraphs.jl/blob/2a644c2b15b444e7f32f73021ec276aa9fc8ba30/src/SimpleGraphs/generators/randgraphs.jl#L187
+"""
+function watts_strogatz(::Type{T}, n::Integer, k::Integer, β::Real;
+                        seed::Int=-1, rng::AbstractRNG=GLOBAL_RNG) where T <: ACSet
+  rng = getRNG(seed, rng)
+  _watts_strogatz(T,n,k,β,rng)
+end
+
+function _watts_strogatz(::Type{T}, n::Integer, k::Integer, β::Real,
+                         rng::AbstractRNG) where T <: ACSet
+  @assert k < n
+
+  # If we have n - 1 neighbors (exactly k/2 on each side), then the graph is
+  # necessarily complete. No need to run the Watts-Strogatz procedure:
+  if k == n - 1 && iseven(k)
+    return complete_graph(T,n)
+  end
+
+  g = T()
+  add_vertices_with_indices!(g, n, k)
+
+  # The ith next vertex, in clockwise order.
+  # (Reduce to zero-based indexing, so the modulo works, by subtracting 1
+  # before and adding 1 after.)
+  @inline target(s, i) = ((s + i - 1) % n) + 1
+
+  # Phase 1: For each step size i, add an edge from each vertex s to the ith
+  # next vertex, in clockwise order.
+
+  for i = 1:div(k, 2), s = 1:n
+    add_edge!(g, s, target(s, i))
+  end
+
+  # Phase 2: For each step size i and each vertex s, consider the edge to the
+  # ith next vertex, in clockwise order. With probability β, delete the edge
+  # and rewire it to any (valid) target, chosen uniformly at random.
+
+  for i = 1:div(k, 2), s = 1:n
+
+    # We only rewire with a probability β, and we only worry about rewiring
+    # if there is some vertex not connected to s; otherwise, the only valid
+    # rewiring is to reconnect to the ith next vertex, and there is no work
+    # to do.
+    (rand(rng) < β && degree(g, s) < n - 1) || continue
+
+    t = target(s, i)
+
+    while true
+      d = rand(1:n)               # Tentative new target
+      d == s && continue          # Self-loops prohibited
+      d == t && break             # Rewired to original target
+      if !(has_edge(g,s,d))       # Was this valid (i.e., unconnected)?
+        add_edge!(g, s, d)        # True rewiring: Add new edge
+        rem_edge!(g, s, t)        # True rewiring: Delete original edge
+        break                     # We found a valid target
+      end
+    end
+
+  end
+  return g
+end
+
 end
diff --git a/src/graphs/Graphs.jl b/src/graphs/Graphs.jl
index 8d44d87ad..300056816 100644
--- a/src/graphs/Graphs.jl
+++ b/src/graphs/Graphs.jl
@@ -7,11 +7,13 @@ include("BipartiteGraphs.jl")
 include("PropertyGraphs.jl")
 include("GraphAlgorithms.jl")
 include("GraphGenerators.jl")
+include("Searching.jl")
 
 @reexport using .BasicGraphs
 @reexport using .BipartiteGraphs
 @reexport using .PropertyGraphs
 @reexport using .GraphAlgorithms
-@reexport using. GraphGenerators
+@reexport using .GraphGenerators
+@reexport using .Searching
 
 end
diff --git a/src/graphs/Searching.jl b/src/graphs/Searching.jl
new file mode 100644
index 000000000..461d1d6b2
--- /dev/null
+++ b/src/graphs/Searching.jl
@@ -0,0 +1,125 @@
+module Searching
+export bfs_parents, bfs_tree, dfs_parents, dfs_tree
+
+using ...CSetDataStructures, ..BasicGraphs
+
+"""
+    tree(parents)
+Convert a parents array into a directed graph.
+"""
+function tree(parents::AbstractVector{Int})
+    n = T(length(parents))
+    t = Graph(n)
+    for (v, u) in enumerate(parents)
+        if u > 0 && u != v
+            add_edge!(t, u, v)
+        end
+    end
+    return t
+end
+
+"""
+    bfs_parents(g, s[; dir=:out])
+
+Perform a breadth-first search of graph `g` starting from vertex `s`.
+Return a vector of parent vertices indexed by vertex. If `dir` is specified,
+use the corresponding edge direction (`:in` and `:out` are acceptable values).
+
+
+
+### Performance
+This implementation is designed to perform well on large graphs. There are
+implementations which are marginally faster in practice for smaller graphs,
+but the performance improvements using this implementation on large graphs
+can be significant.
+"""
+bfs_parents(g::ACSet, s::Int; dir = :out) =
+    (dir == :out) ? _bfs_parents(g, s, outneighbors) : _bfs_parents(g, s, inneighbors)
+
+function _bfs_parents(g::ACSet, source, neighborfn::Function)
+    n = nv(g)
+    visited = falses(n)
+    parents = zeros(Int, nv(g))
+    cur_level = Int[]
+    sizehint!(cur_level, n)
+    next_level = Int[]
+    sizehint!(next_level, n)
+    @inbounds for s in source
+        visited[s] = true
+        push!(cur_level, s)
+        parents[s] = s
+    end
+    while !isempty(cur_level)
+        @inbounds for v in cur_level
+            @inbounds @simd for i in neighborfn(g, v)
+                if !visited[i]
+                    push!(next_level, i)
+                    parents[i] = v
+                    visited[i] = true
+                end
+            end
+        end
+        empty!(cur_level)
+        cur_level, next_level = next_level, cur_level
+        sort!(cur_level)
+    end
+    return parents
+end
+
+"""
+    bfs_tree(g, s[; dir=:out])
+Provide a breadth-first traversal of the graph `g` starting with source vertex `s`,
+and return a directed acyclic graph of vertices in the order they were discovered.
+If `dir` is specified, use the corresponding edge direction (`:in` and `:out` are
+acceptable values).
+"""
+bfs_tree(g::ACSet, s::Integer; dir = :out) = tree(bfs_parents(g, s; dir = dir))
+
+"""
+    dfs_parents(g, s[; dir=:out])
+
+Perform a depth-first search of graph `g` starting from vertex `s`.
+Return a vector of parent vertices indexed by vertex. If `dir` is specified,
+use the corresponding edge direction (`:in` and `:out` are acceptable values).
+
+### Implementation Notes
+This version of DFS is iterative.
+"""
+dfs_parents(g::ACSet, s::Integer; dir=:out) =
+    (dir == :out) ? _dfs_parents(g, s, outneighbors) : _dfs_parents(g, s, inneighbors)
+
+function _dfs_parents(g::ACSet, s::Int, neighborfn::Function)
+    parents = zeros(Int, nv(g))
+    seen = zeros(Bool, nv(g))
+    S = [s]
+    seen[s] = true
+    parents[s] = s
+    while !isempty(S)
+        v = S[end]
+        u = 0
+        for n in neighborfn(g, v)
+            if !seen[n]
+                u = n
+                break
+            end
+        end
+        if u == 0
+            pop!(S)
+        else
+            seen[u] = true
+            push!(S, u)
+            parents[u] = v
+        end
+    end
+    return parents
+end
+
+"""
+    dfs_tree(g, s)
+
+Return a directed acyclic graph based on
+depth-first traversal of the graph `g` starting with source vertex `s`.
+"""
+dfs_tree(g::AbstractGraph, s::Integer; dir=:out) = tree(dfs_parents(g, s; dir=dir))
+
+end
diff --git a/test/graphs/BasicGraphs.jl b/test/graphs/BasicGraphs.jl
index c1be06f8d..3b21e6e83 100644
--- a/test/graphs/BasicGraphs.jl
+++ b/test/graphs/BasicGraphs.jl
@@ -23,6 +23,7 @@ add_edge!(g, 2, 3)
 @test !has_edge(g, 1, 3)
 @test outneighbors(g, 2) == [3]
 @test inneighbors(g, 2) == [1]
+@test degree(g, 2) == 2
 @test collect(all_neighbors(g, 2)) == [1,3]
 
 add_edge!(g, 1, 2)
@@ -31,6 +32,7 @@ add_edge!(g, 1, 2)
 @test collect(edges(g, 1, 2)) == [1,3]
 @test outneighbors(g, 1) == [2,2]
 @test inneighbors(g, 1) == []
+@test degree(g, 1) == 2
 @test LG.DiGraph(g) == LG.path_digraph(3)
 
 g = Graph(4)
diff --git a/test/graphs/GraphGenerators.jl b/test/graphs/GraphGenerators.jl
index 6ea638ee0..dc20d8216 100644
--- a/test/graphs/GraphGenerators.jl
+++ b/test/graphs/GraphGenerators.jl
@@ -64,4 +64,28 @@ g = parallel_arrows(Graph, n)
 @test all(==(1), src(g))
 @test all(==(2), tgt(g))
 
+# Erdos-Renyi graphs
+#-------------------
+
+g = erdos_renyi(Graph, 100, 0.0)
+@test (nv(g), ne(g)) == (100, 0)
+g = erdos_renyi(Graph, 100, 0)
+@test (nv(g), ne(g)) == (100, 0)
+g = erdos_renyi(Graph, 100, 5)
+@test (nv(g), ne(g)) == (100, 5)
+g = erdos_renyi(Graph, 100, 1.0)
+@test g == complete_graph(Graph, 100)
+
+# Expected Degree Graph
+#----------------------
+
+g = expected_degree_graph(Graph, 5 .* ones(100))
+@test nv(g) == 100
+
+# Watts Strogatz Graph
+#---------------------
+
+g = watts_strogatz(Graph, 100, 4, 0.2)
+@test (nv(g), ne(g)) == (100, 200)
+
 end
diff --git a/test/graphs/Searching.jl b/test/graphs/Searching.jl
new file mode 100644
index 000000000..dd809bd87
--- /dev/null
+++ b/test/graphs/Searching.jl
@@ -0,0 +1,41 @@
+module TestSearching
+using Test
+
+using Catlab.Graphs.BasicGraphs, Catlab.Graphs.Searching
+import Catlab.Graphs.Searching: tree
+
+# Tests stolen from LightGraphs
+
+# BFS
+#----
+
+g = Graph(4)
+add_edges!(g, [1,2,1,3], [2,3,3,4])
+z = @inferred(bfs_tree(g, 1))
+t = bfs_parents(g, 1)
+@test t == [1,1,1,3]
+@test nv(z) == 4 && ne(z) == 3 && !has_edge(z,2,3)
+
+g = Graph(5) # house graph
+add_edges!(g, [1,1,2,3,3,4], [2,3,4,4,5,5])
+n = nv(g)
+parents = bfs_parents(g, 1)
+@test length(parents) == n
+t1 = @inferred(bfs_tree(g, 1))
+t2 = tree(parents)
+@test t1 == t2
+@test is_directed(t2)
+@test typeof(t2) <: AbstractGraph
+@test ne(t2) < nv(t2)
+
+# DFS
+#----
+
+g = Graph(4)
+add_edges!(g, [1,2,1,3], [2,3,3,4])
+z = @inferred(dfs_tree(g, 1))
+@test ne(z) == 3 && nv(z) == 4
+@test !has_edge(z, 1, 3)
+@test !is_cyclic(g)
+
+end