From c80f2f852f02fd83b6846678172f5b375c95aab8 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sun, 18 Dec 2022 10:03:54 -0500
Subject: [PATCH 01/35] Implementation for planarity test First step to fixing
 #153

---
 src/planarity.jl  | 452 ++++++++++++++++++++++++++++++++++++++++++++++
 test/planarity.jl | 121 +++++++++++++
 2 files changed, 573 insertions(+)
 create mode 100644 src/planarity.jl
 create mode 100644 test/planarity.jl

diff --git a/src/planarity.jl b/src/planarity.jl
new file mode 100644
index 000000000..14ed66f3f
--- /dev/null
+++ b/src/planarity.jl
@@ -0,0 +1,452 @@
+# Planarity algorithm for Julia graphs. 
+# Algorithm from https://www.uni-konstanz.de/algo/publications/b-lrpt-sub.pdf
+# The implementation is heavily influenced by the recursive implementation in Networkx (https://networkx.org/documentation/stable/_modules/networkx/algorithms/planarity.html)
+
+using DataStructures
+
+import Base: isempty
+
+export is_planar
+
+"""
+    is_planar(g)
+
+Determines whether or not the graph `g` is [planar](https://en.wikipedia.org/wiki/Planar_graph).
+
+Uses the [left-right planarity test](https://en.wikipedia.org/wiki/Left-right_planarity_test).
+
+### References 
+- [Brandes 2009](https://www.uni-konstanz.de/algo/publications/b-lrpt-sub.pdf)
+"""
+function is_planar(g::SimpleGraph)
+    lrp = LRPlanarity(g)
+    lr_planarity!(lrp)
+end
+
+#if g is directed, analyse the undirected version 
+function is_planar(dg::SimpleDiGraph)
+    g = SimpleGraph(dg)
+    is_planar(g)
+end
+
+#aux functions 
+function add_vertices_from!(input_g, output_g)
+    for v in vertices(input_g)
+        add_vertex!(output_g)
+    end
+end
+
+function add_edges_from!(input_g, output_g)
+    for e in edges(input_g)
+        add_edge!(output_g, e)
+    end
+end
+
+#Simple structs to be used in algorithm. Keep private for now. 
+function empty_edge(T)
+    Edge{T}(0, 0)
+end
+
+function isempty(e::Edge{T}) where {T}
+    e.src == zero(T) && e.dst == zero(T)
+end
+
+mutable struct Interval{T}
+    high::Edge{T}
+    low::Edge{T}
+end
+
+function empty_interval(T)
+    Interval(empty_edge(T), empty_edge(T))
+end
+
+function isempty(interval::Interval)
+    isempty(interval.high) && isempty(interval.low)
+end
+
+function conflicting(interval::Interval{T}, b, lrp_state) where {T}
+    !isempty(interval) && (lrp_state.lowpt[interval.high] > lrp_state.lowpt[b])
+end
+
+mutable struct Pair{T}
+    L::Interval{T}
+    R::Interval{T}
+end
+
+function empty_pair(T)
+    Pair(empty_interval(T), empty_interval(T))
+end
+
+function swap!(self::Pair)
+    #Swap left and right intervals
+    temp = self.L
+    self.L = self.R
+    self.R = temp
+end
+
+function root_pair(T)
+    #returns the "root pair" of type T 
+    e = Edge{T}(0, 0)
+    Pair(Interval(e, e), Interval(e, e))
+end
+
+function isempty(p::Pair)
+    isempty(p.L) && isempty(p.R)
+end
+
+struct LRPlanarity{T}
+    #State class for the planarity test 
+    #We index by Edge structs throughout as it is easier than switching between
+    #Edges and tuples
+    G::SimpleGraph{T} #Copy of the input graph
+    roots::Vector{T} #Vector of roots for disconnected graphs. Normally size(roots, 1) == 1
+    height::DefaultDict{T,Int64} #DefaultDict of heights <: Int, indexed by node. default is -1
+    lowpt::Dict{Edge{T},Int64} #Dict of low points, indexed by Edge
+    lowpt2::Dict{Edge{T},Int64} #Dict of low points (copy), indexed by Edge 
+    nesting_depth::Dict{Edge{T},Int64} #Dict of nesting depths, indexed by Edge
+    parent_edge::DefaultDict{T,Edge{T}} #Dict of parent edges, indexed by node
+    DG::SimpleDiGraph{T} #Directed graph for the orientation phase
+    adjs::Dict{T,Vector{T}} #Dict of neighbors of nodes, indexed by node
+    ordered_adjs::Dict{T,Vector{T}} #Dict of neighbors of nodes sorted by nesting depth, indexed by node
+    ref::DefaultDict{Edge{T},Edge{T}} #DefaultDict of Edges, indexed by Edge
+    side::DefaultDict{Edge{T},Int8} #DefaultDict of +/- 1, indexed by edge
+    S::Stack{Pair{T}} #Stack of tuples of Edges
+    stack_bottom::Dict{Edge{T},Pair{T}} #Dict of Tuples of Edges, indexed by Edge
+    lowpt_edge::Dict{Edge{T},Edge{T}} #Dict of Edges, indexed by Edge 
+    left_ref::Dict{T,Edge{T}} #Dict of Edges, indexed by node
+    right_ref::Dict{T,Edge{T}} #Dict of Edges, indexed by node
+    # skip embedding for now 
+end
+
+#outer constructor for LRPlanarity
+function LRPlanarity(g)
+    #record nodetype of g
+    T = eltype(g)
+    # copy G without adding self-loops
+    output_g = SimpleGraph{T}()
+    add_vertices_from!(g, output_g)
+    for e in edges(g)
+        if e.src != e.dst
+            add_edge!(output_g, e)
+        end
+    end
+
+    N = nv(output_g)
+
+    roots = T[]
+
+    # distance from tree root
+    height = DefaultDict{T,Int64}(-1)
+
+    lowpt = Dict{Edge{T},Int64}()  # height of lowest return point of an edge
+    lowpt2 = Dict{Edge{T},Int64}()  # height of second lowest return point
+    nesting_depth = Dict{Edge{T},Int64}()  # for nesting order
+
+    # None == Edge(0, 0)  for our type-stable algo
+
+    parent_edge = DefaultDict{T,Edge{T}}(Edge(zero(T), zero(T)))
+
+    # oriented DFS graph
+    DG = SimpleDiGraph{T}()
+    add_vertices_from!(g, DG)
+
+    adjs = Dict{T,Vector{T}}()
+    ordered_adjs = Dict{T,Vector{T}}()
+
+    ref = DefaultDict{Edge{T},Edge{T}}(empty_edge(T))
+    side = DefaultDict{Edge{T},Int8}(one(Int8))
+
+    # stack of conflict pairs
+    S = Stack{Pair{T}}()
+    stack_bottom = Dict{Edge{T},Pair{T}}()
+    lowpt_edge = Dict{Edge{T},Edge{T}}()
+    left_ref = Dict{T,Edge{T}}()
+    right_ref = Dict{T,Edge{T}}()
+
+    #self.embedding = PlanarEmbedding()
+    LRPlanarity(
+        g,
+        roots,
+        height,
+        lowpt,
+        lowpt2,
+        nesting_depth,
+        parent_edge,
+        DG,
+        adjs,
+        ordered_adjs,
+        ref,
+        side,
+        S,
+        stack_bottom,
+        lowpt_edge,
+        left_ref,
+        right_ref,
+    )
+end
+
+function lowest(self::Pair, planarity_state::LRPlanarity)
+    #Returns the lowest lowpoint of a conflict pair
+    if isempty(self.L)
+        return planarity_state.lowpt[self.R.low]
+    end
+
+    if isempty(self.R)
+        return planarity_state.lowpt[self.L.low]
+    end
+
+    return min(planarity_state.lowpt[self.L.low], planarity_state.lowpt[self.R.low])
+end
+
+function lr_planarity!(self::LRPlanarity{T}) where {T}
+    V = nv(self.G)
+    E = ne(self.G)
+
+    if V > 2 && (E > (3V - 6))
+        # graph is not planar
+        return false
+    end
+
+
+    # make adjacency lists for dfs
+    for v = 1:nv(self.G) #for all vertices in G,
+        self.adjs[v] = neighbors(self.G, v) ##neighbourhood of v
+    end
+
+    # orientation of the graph by depth first search traversal
+    for v in vertices(self.G)
+        if self.height[v] == -one(T) #using -1 rather than nothing for type stability. 
+            self.height[v] = zero(T)
+            push!(self.roots, v)
+            dfs_orientation!(self, v)
+        end
+    end
+
+    #Testing stage
+    #First, sort the ordered_adjs by nesting depth
+    for v = 1:nv(self.G) #for all vertices in G,
+        #get neighboring nodes
+        neighboring_nodes = T[]
+        neighboring_nesting_depths = Int64[]
+        for (k, value) in self.nesting_depth
+            if k.src == v
+                push!(neighboring_nodes, k.dst)
+                push!(neighboring_nesting_depths, value)
+            end
+        end
+        neighboring_nodes .= neighboring_nodes[sortperm(neighboring_nesting_depths)]
+        self.ordered_adjs[v] = neighboring_nodes
+    end
+    for s in self.roots
+        if !dfs_testing!(self, s)
+            return false
+        end
+    end
+
+    #if the algorithm finishes, the graph is planar. Return true
+    true
+end
+
+
+
+function dfs_orientation!(self::LRPlanarity, v)
+    # get the parent edge of v. 
+    # if v is a root, the parent_edge dict 
+    # will return Edge(0, 0)
+    e = self.parent_edge[v] #get the parent edge of v. 
+    #orient all edges in graph recursively 
+    for w in self.adjs[v]
+        # see if vw = Edge(v, w) has already been oriented 
+        vw = Edge(v, w)
+        wv = Edge(w, v) #Need to consider approach from both direction
+        if vw in edges(self.DG) || wv in edges(self.DG)
+            continue
+        end
+
+        #otherwise, appended to DG
+        add_edge!(self.DG, vw)
+
+        #record lowpoints
+        self.lowpt[vw] = self.height[v]
+        self.lowpt2[vw] = self.height[v]
+        #if height == -1, i.e. we are at a tree edge, then 
+        # record the height accordingly
+        if self.height[w] == -1 ##tree edge
+            self.parent_edge[w] = vw
+            self.height[w] = self.height[v] + 1
+            dfs_orientation!(self, w)
+        else
+            #at a back edge - no need to 
+            #go through a DFS
+            self.lowpt[vw] = self.height[w]
+        end
+
+        #determine nesting depth with formulae from Brandes
+        #note that this will only be carried out 
+        #once per edge 
+        # if the edge is chordal, use the alternative formula
+        self.nesting_depth[vw] = 2 * self.lowpt[vw]
+        if self.lowpt2[vw] < self.height[v]
+            #chordal 
+            self.nesting_depth[vw] += 1
+        end
+
+        #update lowpoints of parent 
+        if !isempty(e) #if e != root
+            if self.lowpt[vw] < self.lowpt[e]
+                self.lowpt2[e] = min(self.lowpt[e], self.lowpt2[vw])
+                self.lowpt[e] = self.lowpt[vw]
+            elseif self.lowpt[vw] > self.lowpt[e]
+                self.lowpt2[e] = min(self.lowpt2[e], self.lowpt[vw])
+            else
+                self.lowpt2[e] = min(self.lowpt2[e], self.lowpt2[vw])
+            end
+        end
+    end
+
+end
+
+function dfs_testing!(self, v)
+    T = typeof(v)
+    e = self.parent_edge[v]
+    for w in self.ordered_adjs[v] #already ordered 
+        ei = Edge(v, w)
+        if !isempty(self.S) #stack is not empty
+            self.stack_bottom[ei] = top(self.S)
+        else #stack is empty
+            self.stack_bottom[ei] = root_pair(T)
+        end
+
+        if ei == self.parent_edge[ei.dst] #tree edge
+            if !dfs_testing!(self, ei.dst) #half if testing fails
+                return false
+            end
+        else #back edge
+            self.lowpt_edge[ei] = ei
+            push!(self.S, Pair(empty_interval(T), Interval(ei, ei)))
+        end
+
+        #integrate new return edges 
+        if self.lowpt[ei] < self.height[v] #ei has return edge 
+            e1 = Edge(v, first(self.ordered_adjs[v]))
+            if ei == e1
+                self.lowpt_edge[e] = self.lowpt_edge[ei] #in Brandes this is e <- e1. Corrected in Python source?
+            else
+                #add contraints (algo 4)
+                if !edge_constraints!(self, ei, e) #half if fails
+                    return false
+                end
+            end
+        end
+    end
+
+    #remove back edges returning to parent 
+    if !isempty(e)#v is not root
+        u = e.src
+        #trim edges ending at parent u, algo 5 
+        trim_back_edges!(self, u)
+        #side of e is side of highest returning edge 
+        if self.lowpt[e] < self.height[u] #e has return edge
+            hl = top(self.S).L.high
+            hr = top(self.S).R.high
+            if !isempty(hl) && (isempty(hr) || (self.lowpt[hl] > self.lowpt[hr]))
+                self.ref[e] = hl
+            else
+                self.ref[e] = hr
+            end
+        end
+    end
+    true
+end
+
+function edge_constraints!(self, ei, e)
+    T = eltype(ei)
+    P = empty_pair(T)
+    #merge return edges of ei into P.R
+    while top(self.S) != self.stack_bottom[ei]
+        Q = pop!(self.S)
+        if !isempty(Q.L)
+            swap!(Q)
+        end
+        if !isempty(Q.L) #not planar
+            return false
+        else
+            if self.lowpt[Q.R.low] > self.lowpt[e] #merge intervals
+                if isempty(P.R) #topmost interval
+                    P.R.high = Q.R.high
+                else
+                    self.ref[P.R.low] = Q.R.high
+                end
+                P.R.low = Q.R.low
+            else #align
+                self.ref[Q.R.low] = self.lowpt_edge[e]
+            end
+        end
+    end
+
+    #merge conflicting return edges of <e...> into P.LRPlanarity
+    while conflicting(top(self.S).L, ei, self) || conflicting(top(self.S).R, ei, self)
+        Q = pop!(self.S)
+        if conflicting(Q.R, ei, self)
+            swap!(Q)
+        end
+        if conflicting(Q.R, ei, self) #not planar
+            return false
+        else #merge interval below into P.R
+            self.ref[P.R.low] = Q.R.high
+            if !isempty(Q.R.low)
+                P.R.low = Q.R.low
+            end
+        end
+        if isempty(P.L) #topmost interval
+            P.L.high = Q.L.high
+        else
+            self.ref[P.L.low] = Q.L.high
+        end
+        P.L.low = Q.L.low
+    end
+    if !isempty(P)
+        push!(self.S, P)
+    end
+    true
+end
+
+function trim_back_edges!(self, u)
+    #trim back edges ending at u 
+    #drop entire conflict pairs 
+    while !isempty(self.S) && (lowest(top(self.S), self) == self.height[u])
+        P = pop!(self.S)
+        if !isempty(P.L.low)
+            self.side[P.L.low] = -1
+        end
+    end
+
+    if !isempty(self.S) #one more conflict pair to consider 
+        P = pop!(self.S)
+        #trim left interval 
+        while !isempty(P.L.high) && P.L.high.dst == u
+            P.L.high = self.ref[P.L.high]
+        end
+
+        if isempty(P.L.high) && !isempty(P.L.low) #just emptied
+            self.ref[P.L.low] = P.R.low
+            self.side[P.L.low] = -1
+            T = typeof(u)
+            P.L.low = empty_edge(T)
+        end
+
+        #trim right interval 
+        while !isempty(P.R.high) && P.R.high.dst == u
+            P.R.high = self.ref[P.R.high]
+        end
+
+        if isempty(P.R.high) && !isempty(P.R.low) #just emptied
+            self.ref[P.R.low] = P.L.low
+            self.side[P.R.low] = -1
+            T = typeof(u)
+            P.R.low = empty_edge(T)
+        end
+        push!(self.S, P)
+    end
+
+end
diff --git a/test/planarity.jl b/test/planarity.jl
new file mode 100644
index 000000000..b51c2a211
--- /dev/null
+++ b/test/planarity.jl
@@ -0,0 +1,121 @@
+@testset "Planarity tests" begin 
+    @testset "Aux functions tests" begin
+        g = complete_graph(10)
+        f = SimpleGraph()
+        add_vertices_from!(g, f)
+        @test nv(g) == nv(f)
+        add_edges_from!(g, f)
+        @test g == f
+    end
+    
+    @testset "LRP constructor tests" begin
+        function lrp_constructor_test(g)
+            try
+                lrp = LRPlanarity(g)
+                return true
+            catch e
+                return false
+            end
+        end
+        g = SimpleGraph(10, 10)
+        @test lrp_constructor_test(g)
+        g = SimpleGraph{Int8}(10, 10)
+        @test lrp_constructor_test(g)
+    end
+    
+    @testset "DFS orientation tests" begin
+        dfs_test_edges = [
+            (1, 2),
+            (1, 4),
+            (1, 6),
+            (2, 3),
+            (2, 4),
+            (2, 5),
+            (3, 5),
+            (3, 6),
+            (4, 5),
+            (4, 6),
+            (5, 6),
+        ]
+    
+        dfs_g = Graph(6)
+        for edge in dfs_test_edges
+            add_edge!(dfs_g, edge)
+        end
+    
+        self = LRPlanarity(dfs_g)
+        #want to test dfs orientation 
+        # make adjacency lists for dfs
+        for v = 1:nv(self.G) #for all vertices in G,
+            self.adjs[v] = neighbors(self.G, v) ##neighbourhood of v
+        end
+        T = eltype(self.G)
+    
+        # orientation of the graph by depth first search traversal
+        for v in vertices(self.G)
+            if self.height[v] == -one(T) #using -1 rather than nothing for type stability. 
+                self.height[v] = zero(T)
+                push!(self.roots, v)
+                dfs_orientation!(self, v)
+            end
+        end
+    
+        #correct data
+        parent_edges = Dict([
+            (1, Edge(0, 0)),
+            (2, Edge(1, 2)),
+            (3, Edge(2, 3)),
+            (4, Edge(5, 4)),
+            (5, Edge(3, 5)),
+            (6, Edge(4, 6)),
+        ])
+    
+        @test parent_edges == self.parent_edge
+    
+        correct_heights = Dict([(1, 0), (2, 1), (3, 2), (4, 4), (5, 3), (6, 5)])
+    
+        @test correct_heights == self.height
+    end
+    
+    @testset "Planarity results" begin
+        #Construct example planar graph
+        planar_edges = [(1, 2), (1, 3), (2, 4), (2, 7), (3, 4), (3, 5), (4, 6), (4, 7), (5, 6)]
+    
+        g = Graph(7)
+    
+        for edge in planar_edges
+            add_edge!(g, edge)
+        end
+    
+        @test is_planar(g) == true
+    
+        #another planar graph 
+        cl = circular_ladder_graph(8)
+        @test is_planar(cl) == true
+    
+        # one more planar graph 
+        w = wheel_graph(10)
+        @test is_planar(w) == true
+    
+        #Construct some non-planar graphs
+        g = complete_graph(10)
+        @test is_planar(g) == false
+    
+        petersen = smallgraph(:petersen)
+        @test is_planar(petersen) == false
+    
+        d = smallgraph(:desargues)
+        @test is_planar(d) == false
+    
+        h = smallgraph(:heawood)
+        @test is_planar(h) == false
+    
+        mb = smallgraph(:moebiuskantor)
+        @test is_planar(mb) == false
+    
+        #Directed, planar example
+        dg = SimpleDiGraphFromIterator([Edge(1, 2), Edge(2, 3), Edge(3, 1)])
+    
+        @test is_planar(dg) == true
+    end
+end
\ No newline at end of file

From 9a5179cccf2eed51db2ae8febdc1637e2ce8d5e6 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Mon, 19 Dec 2022 20:21:05 -0500
Subject: [PATCH 02/35] Tweak planarity check; add PMFG algorithm

---
 src/Graphs.jl             |  7 ++++-
 src/planarity.jl          | 64 +++++++++++++++++----------------------
 src/spanningtrees/pmfg.jl | 34 +++++++++++++++++++++
 test/planarity.jl         | 64 ++++++++++++++++++++-------------------
 4 files changed, 101 insertions(+), 68 deletions(-)
 create mode 100644 src/spanningtrees/pmfg.jl

diff --git a/src/Graphs.jl b/src/Graphs.jl
index 171204112..66b95fab3 100644
--- a/src/Graphs.jl
+++ b/src/Graphs.jl
@@ -413,8 +413,11 @@ export
     independent_set,
 
     # vertexcover
-    vertex_cover
+    vertex_cover,
 
+    # planarity
+    is_planar,
+    pmfg
 """
     Graphs
 
@@ -528,6 +531,8 @@ include("vertexcover/degree_vertex_cover.jl")
 include("vertexcover/random_vertex_cover.jl")
 include("Experimental/Experimental.jl")
 include("Parallel/Parallel.jl")
+include("planarity.jl")
+include("spanningtrees/pmfg.jl")
 
 using .LinAlg
 end # module
diff --git a/src/planarity.jl b/src/planarity.jl
index 14ed66f3f..df74384c6 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -2,12 +2,9 @@
 # Algorithm from https://www.uni-konstanz.de/algo/publications/b-lrpt-sub.pdf
 # The implementation is heavily influenced by the recursive implementation in Networkx (https://networkx.org/documentation/stable/_modules/networkx/algorithms/planarity.html)
 
-using DataStructures
-
+import DataStructures: DefaultDict, Stack, top
 import Base: isempty
 
-export is_planar
-
 """
     is_planar(g)
 
@@ -20,13 +17,13 @@ Uses the [left-right planarity test](https://en.wikipedia.org/wiki/Left-right_pl
 """
 function is_planar(g::SimpleGraph)
     lrp = LRPlanarity(g)
-    lr_planarity!(lrp)
+    return lr_planarity!(lrp)
 end
 
 #if g is directed, analyse the undirected version 
 function is_planar(dg::SimpleDiGraph)
     g = SimpleGraph(dg)
-    is_planar(g)
+    return is_planar(g)
 end
 
 #aux functions 
@@ -44,11 +41,11 @@ end
 
 #Simple structs to be used in algorithm. Keep private for now. 
 function empty_edge(T)
-    Edge{T}(0, 0)
+    return Edge{T}(0, 0)
 end
 
 function isempty(e::Edge{T}) where {T}
-    e.src == zero(T) && e.dst == zero(T)
+    return e.src == zero(T) && e.dst == zero(T)
 end
 
 mutable struct Interval{T}
@@ -57,41 +54,41 @@ mutable struct Interval{T}
 end
 
 function empty_interval(T)
-    Interval(empty_edge(T), empty_edge(T))
+    return Interval(empty_edge(T), empty_edge(T))
 end
 
 function isempty(interval::Interval)
-    isempty(interval.high) && isempty(interval.low)
+    return isempty(interval.high) && isempty(interval.low)
 end
 
 function conflicting(interval::Interval{T}, b, lrp_state) where {T}
-    !isempty(interval) && (lrp_state.lowpt[interval.high] > lrp_state.lowpt[b])
+    return !isempty(interval) && (lrp_state.lowpt[interval.high] > lrp_state.lowpt[b])
 end
 
-mutable struct Pair{T}
+mutable struct ConflictPair{T}
     L::Interval{T}
     R::Interval{T}
 end
 
 function empty_pair(T)
-    Pair(empty_interval(T), empty_interval(T))
+    return ConflictPair(empty_interval(T), empty_interval(T))
 end
 
-function swap!(self::Pair)
+function swap!(self::ConflictPair)
     #Swap left and right intervals
     temp = self.L
     self.L = self.R
-    self.R = temp
+    return self.R = temp
 end
 
 function root_pair(T)
     #returns the "root pair" of type T 
     e = Edge{T}(0, 0)
-    Pair(Interval(e, e), Interval(e, e))
+    return ConflictPair(Interval(e, e), Interval(e, e))
 end
 
-function isempty(p::Pair)
-    isempty(p.L) && isempty(p.R)
+function isempty(p::ConflictPair)
+    return isempty(p.L) && isempty(p.R)
 end
 
 struct LRPlanarity{T}
@@ -110,8 +107,8 @@ struct LRPlanarity{T}
     ordered_adjs::Dict{T,Vector{T}} #Dict of neighbors of nodes sorted by nesting depth, indexed by node
     ref::DefaultDict{Edge{T},Edge{T}} #DefaultDict of Edges, indexed by Edge
     side::DefaultDict{Edge{T},Int8} #DefaultDict of +/- 1, indexed by edge
-    S::Stack{Pair{T}} #Stack of tuples of Edges
-    stack_bottom::Dict{Edge{T},Pair{T}} #Dict of Tuples of Edges, indexed by Edge
+    S::Stack{ConflictPair{T}} #Stack of tuples of Edges
+    stack_bottom::Dict{Edge{T},ConflictPair{T}} #Dict of Tuples of Edges, indexed by Edge
     lowpt_edge::Dict{Edge{T},Edge{T}} #Dict of Edges, indexed by Edge 
     left_ref::Dict{T,Edge{T}} #Dict of Edges, indexed by node
     right_ref::Dict{T,Edge{T}} #Dict of Edges, indexed by node
@@ -142,7 +139,7 @@ function LRPlanarity(g)
     lowpt2 = Dict{Edge{T},Int64}()  # height of second lowest return point
     nesting_depth = Dict{Edge{T},Int64}()  # for nesting order
 
-    # None == Edge(0, 0)  for our type-stable algo
+    # None == Edge(-1, 0)  for our type-stable algo
 
     parent_edge = DefaultDict{T,Edge{T}}(Edge(zero(T), zero(T)))
 
@@ -157,14 +154,14 @@ function LRPlanarity(g)
     side = DefaultDict{Edge{T},Int8}(one(Int8))
 
     # stack of conflict pairs
-    S = Stack{Pair{T}}()
-    stack_bottom = Dict{Edge{T},Pair{T}}()
+    S = Stack{ConflictPair{T}}()
+    stack_bottom = Dict{Edge{T},ConflictPair{T}}()
     lowpt_edge = Dict{Edge{T},Edge{T}}()
     left_ref = Dict{T,Edge{T}}()
     right_ref = Dict{T,Edge{T}}()
 
     #self.embedding = PlanarEmbedding()
-    LRPlanarity(
+    return LRPlanarity(
         g,
         roots,
         height,
@@ -185,7 +182,7 @@ function LRPlanarity(g)
     )
 end
 
-function lowest(self::Pair, planarity_state::LRPlanarity)
+function lowest(self::ConflictPair, planarity_state::LRPlanarity)
     #Returns the lowest lowpoint of a conflict pair
     if isempty(self.L)
         return planarity_state.lowpt[self.R.low]
@@ -207,9 +204,8 @@ function lr_planarity!(self::LRPlanarity{T}) where {T}
         return false
     end
 
-
     # make adjacency lists for dfs
-    for v = 1:nv(self.G) #for all vertices in G,
+    for v in 1:nv(self.G) #for all vertices in G,
         self.adjs[v] = neighbors(self.G, v) ##neighbourhood of v
     end
 
@@ -224,7 +220,7 @@ function lr_planarity!(self::LRPlanarity{T}) where {T}
 
     #Testing stage
     #First, sort the ordered_adjs by nesting depth
-    for v = 1:nv(self.G) #for all vertices in G,
+    for v in 1:nv(self.G) #for all vertices in G,
         #get neighboring nodes
         neighboring_nodes = T[]
         neighboring_nesting_depths = Int64[]
@@ -244,11 +240,9 @@ function lr_planarity!(self::LRPlanarity{T}) where {T}
     end
 
     #if the algorithm finishes, the graph is planar. Return true
-    true
+    return true
 end
 
-
-
 function dfs_orientation!(self::LRPlanarity, v)
     # get the parent edge of v. 
     # if v is a root, the parent_edge dict 
@@ -303,7 +297,6 @@ function dfs_orientation!(self::LRPlanarity, v)
             end
         end
     end
-
 end
 
 function dfs_testing!(self, v)
@@ -323,7 +316,7 @@ function dfs_testing!(self, v)
             end
         else #back edge
             self.lowpt_edge[ei] = ei
-            push!(self.S, Pair(empty_interval(T), Interval(ei, ei)))
+            push!(self.S, ConflictPair(empty_interval(T), Interval(ei, ei)))
         end
 
         #integrate new return edges 
@@ -356,7 +349,7 @@ function dfs_testing!(self, v)
             end
         end
     end
-    true
+    return true
 end
 
 function edge_constraints!(self, ei, e)
@@ -408,7 +401,7 @@ function edge_constraints!(self, ei, e)
     if !isempty(P)
         push!(self.S, P)
     end
-    true
+    return true
 end
 
 function trim_back_edges!(self, u)
@@ -448,5 +441,4 @@ function trim_back_edges!(self, u)
         end
         push!(self.S, P)
     end
-
 end
diff --git a/src/spanningtrees/pmfg.jl b/src/spanningtrees/pmfg.jl
new file mode 100644
index 000000000..9b4517240
--- /dev/null
+++ b/src/spanningtrees/pmfg.jl
@@ -0,0 +1,34 @@
+"""
+	```pmfg(g, dmtx = weights(g)```
+
+Compete the Planar Maximally Filtered Graph (PMFG) of `g`, given an optional set of weights `dmtx` that form a distance matrix.
+
+### Examples
+```
+using Graphs, SimpleWeightedGraphs
+N = 20
+M = Symmetric(randn(N, N))
+g = SimpleWeightedGraph(M) 
+p_g = pmfg(g)
+```
+
+### References 
+- Tuminello et al. 2005
+"""
+function pmfg(g, dmtx=weights(g))
+    T = eltype(g)
+    N = nv(g)
+    out = SimpleGraph{T}(N)
+    #create list of edges w/ weights
+    edge_list = collect(edges(g))
+    sort!(edge_list; by=x -> x.weight) #leaves biggest weight last
+    while !isempty(edge_list)
+        new_edge_weighted = pop!(edge_list)
+        new_edge = Edge(new_edge_weighted.src, new_edge_weighted.dst)
+        add_edge!(out, new_edge)
+        if !is_planar(out)
+            rem_edge!(out, new_edge)
+        end
+    end
+    return out
+end
diff --git a/test/planarity.jl b/test/planarity.jl
index b51c2a211..958c3be0c 100644
--- a/test/planarity.jl
+++ b/test/planarity.jl
@@ -1,17 +1,17 @@
-@testset "Planarity tests" begin 
+@testset "Planarity tests" begin
     @testset "Aux functions tests" begin
         g = complete_graph(10)
         f = SimpleGraph()
-        add_vertices_from!(g, f)
+        Graphs.add_vertices_from!(g, f)
         @test nv(g) == nv(f)
-        add_edges_from!(g, f)
+        Graphs.add_edges_from!(g, f)
         @test g == f
     end
-    
+
     @testset "LRP constructor tests" begin
         function lrp_constructor_test(g)
             try
-                lrp = LRPlanarity(g)
+                lrp = Graphs.LRPlanarity(g)
                 return true
             catch e
                 return false
@@ -22,7 +22,7 @@
         g = SimpleGraph{Int8}(10, 10)
         @test lrp_constructor_test(g)
     end
-    
+
     @testset "DFS orientation tests" begin
         dfs_test_edges = [
             (1, 2),
@@ -37,29 +37,29 @@
             (4, 6),
             (5, 6),
         ]
-    
+
         dfs_g = Graph(6)
         for edge in dfs_test_edges
             add_edge!(dfs_g, edge)
         end
-    
-        self = LRPlanarity(dfs_g)
+
+        self = Graphs.LRPlanarity(dfs_g)
         #want to test dfs orientation 
         # make adjacency lists for dfs
-        for v = 1:nv(self.G) #for all vertices in G,
+        for v in 1:nv(self.G) #for all vertices in G,
             self.adjs[v] = neighbors(self.G, v) ##neighbourhood of v
         end
         T = eltype(self.G)
-    
+
         # orientation of the graph by depth first search traversal
         for v in vertices(self.G)
             if self.height[v] == -one(T) #using -1 rather than nothing for type stability. 
                 self.height[v] = zero(T)
                 push!(self.roots, v)
-                dfs_orientation!(self, v)
+                Graphs.dfs_orientation!(self, v)
             end
         end
-    
+
         #correct data
         parent_edges = Dict([
             (1, Edge(0, 0)),
@@ -69,53 +69,55 @@
             (5, Edge(3, 5)),
             (6, Edge(4, 6)),
         ])
-    
+
         @test parent_edges == self.parent_edge
-    
+
         correct_heights = Dict([(1, 0), (2, 1), (3, 2), (4, 4), (5, 3), (6, 5)])
-    
+
         @test correct_heights == self.height
     end
-    
+
     @testset "Planarity results" begin
         #Construct example planar graph
-        planar_edges = [(1, 2), (1, 3), (2, 4), (2, 7), (3, 4), (3, 5), (4, 6), (4, 7), (5, 6)]
-    
+        planar_edges = [
+            (1, 2), (1, 3), (2, 4), (2, 7), (3, 4), (3, 5), (4, 6), (4, 7), (5, 6)
+        ]
+
         g = Graph(7)
-    
+
         for edge in planar_edges
             add_edge!(g, edge)
         end
-    
+
         @test is_planar(g) == true
-    
+
         #another planar graph 
         cl = circular_ladder_graph(8)
         @test is_planar(cl) == true
-    
+
         # one more planar graph 
         w = wheel_graph(10)
         @test is_planar(w) == true
-    
+
         #Construct some non-planar graphs
         g = complete_graph(10)
         @test is_planar(g) == false
-    
+
         petersen = smallgraph(:petersen)
         @test is_planar(petersen) == false
-    
+
         d = smallgraph(:desargues)
         @test is_planar(d) == false
-    
+
         h = smallgraph(:heawood)
         @test is_planar(h) == false
-    
+
         mb = smallgraph(:moebiuskantor)
         @test is_planar(mb) == false
-    
+
         #Directed, planar example
         dg = SimpleDiGraphFromIterator([Edge(1, 2), Edge(2, 3), Edge(3, 1)])
-    
+
         @test is_planar(dg) == true
     end
-end
\ No newline at end of file
+end

From 2188295a3f0de54af4579389cee4407245b6b8d5 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Mon, 19 Dec 2022 20:22:58 -0500
Subject: [PATCH 03/35] Fix missing tests code

---
 test/runtests.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 786bdb949..666b5ca58 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -117,6 +117,7 @@ tests = [
     "vertexcover/degree_vertex_cover",
     "vertexcover/random_vertex_cover",
     "experimental/experimental",
+    "planarity",
 ]
 
 @testset verbose = true "Graphs" begin

From 392336c1fe2c8df6ab5c17ead7afcc11af2b4564 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Mon, 19 Dec 2022 20:29:12 -0500
Subject: [PATCH 04/35] Fix type in planarity test

---
 src/planarity.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index df74384c6..5ea4d5304 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -139,7 +139,7 @@ function LRPlanarity(g)
     lowpt2 = Dict{Edge{T},Int64}()  # height of second lowest return point
     nesting_depth = Dict{Edge{T},Int64}()  # for nesting order
 
-    # None == Edge(-1, 0)  for our type-stable algo
+    # None == Edge(0, 0)  for our type-stable algo
 
     parent_edge = DefaultDict{T,Edge{T}}(Edge(zero(T), zero(T)))
 

From 3540b780b26d041d2d3f8f5f2817c37946597518 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Mon, 19 Dec 2022 20:34:12 -0500
Subject: [PATCH 05/35] Fix DataStructures deprecations

---
 src/planarity.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index 5ea4d5304..771b6e490 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -2,7 +2,7 @@
 # Algorithm from https://www.uni-konstanz.de/algo/publications/b-lrpt-sub.pdf
 # The implementation is heavily influenced by the recursive implementation in Networkx (https://networkx.org/documentation/stable/_modules/networkx/algorithms/planarity.html)
 
-import DataStructures: DefaultDict, Stack, top
+import DataStructures: DefaultDict, Stack 
 import Base: isempty
 
 """
@@ -305,7 +305,7 @@ function dfs_testing!(self, v)
     for w in self.ordered_adjs[v] #already ordered 
         ei = Edge(v, w)
         if !isempty(self.S) #stack is not empty
-            self.stack_bottom[ei] = top(self.S)
+            self.stack_bottom[ei] = first(self.S)
         else #stack is empty
             self.stack_bottom[ei] = root_pair(T)
         end
@@ -340,8 +340,8 @@ function dfs_testing!(self, v)
         trim_back_edges!(self, u)
         #side of e is side of highest returning edge 
         if self.lowpt[e] < self.height[u] #e has return edge
-            hl = top(self.S).L.high
-            hr = top(self.S).R.high
+            hl = first(self.S).L.high
+            hr = first(self.S).R.high
             if !isempty(hl) && (isempty(hr) || (self.lowpt[hl] > self.lowpt[hr]))
                 self.ref[e] = hl
             else
@@ -356,7 +356,7 @@ function edge_constraints!(self, ei, e)
     T = eltype(ei)
     P = empty_pair(T)
     #merge return edges of ei into P.R
-    while top(self.S) != self.stack_bottom[ei]
+    while first(self.S) != self.stack_bottom[ei]
         Q = pop!(self.S)
         if !isempty(Q.L)
             swap!(Q)
@@ -378,7 +378,7 @@ function edge_constraints!(self, ei, e)
     end
 
     #merge conflicting return edges of <e...> into P.LRPlanarity
-    while conflicting(top(self.S).L, ei, self) || conflicting(top(self.S).R, ei, self)
+    while conflicting(first(self.S).L, ei, self) || conflicting(first(self.S).R, ei, self)
         Q = pop!(self.S)
         if conflicting(Q.R, ei, self)
             swap!(Q)
@@ -407,7 +407,7 @@ end
 function trim_back_edges!(self, u)
     #trim back edges ending at u 
     #drop entire conflict pairs 
-    while !isempty(self.S) && (lowest(top(self.S), self) == self.height[u])
+    while !isempty(self.S) && (lowest(first(self.S), self) == self.height[u])
         P = pop!(self.S)
         if !isempty(P.L.low)
             self.side[P.L.low] = -1

From c4c56fb0fd39ecef309800406b7a3010c1433d52 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Mon, 2 Jan 2023 21:39:53 +0000
Subject: [PATCH 06/35] Fix planarity typos; add PMFG algo

---
 src/planarity.jl          | 41 +++++++++++++++---------
 src/spanningtrees/pmfg.jl | 67 +++++++++++++++++++++++++++++++--------
 2 files changed, 79 insertions(+), 29 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index 771b6e490..f91e8956f 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -15,7 +15,7 @@ Uses the [left-right planarity test](https://en.wikipedia.org/wiki/Left-right_pl
 ### References 
 - [Brandes 2009](https://www.uni-konstanz.de/algo/publications/b-lrpt-sub.pdf)
 """
-function is_planar(g::SimpleGraph)
+function is_planar(g)
     lrp = LRPlanarity(g)
     return lr_planarity!(lrp)
 end
@@ -39,6 +39,20 @@ function add_edges_from!(input_g, output_g)
     end
 end
 
+function add_simple_edges_from!(input_g, output_g)
+    for e in edges(input_g)
+        add_edge!(output_g, e.dst, e.src)
+    end
+end
+
+#= function add__simple_edges_from!(input_g::SimpleWeightedGraph, output_g)
+    for e in edges(input_g)
+        if e.dst != e.src
+            add_edge!(output_g, e.src, e.dst)
+        end
+    end
+end =#
+
 #Simple structs to be used in algorithm. Keep private for now. 
 function empty_edge(T)
     return Edge{T}(0, 0)
@@ -91,7 +105,7 @@ function isempty(p::ConflictPair)
     return isempty(p.L) && isempty(p.R)
 end
 
-struct LRPlanarity{T}
+struct LRPlanarity{T<:Integer} 
     #State class for the planarity test 
     #We index by Edge structs throughout as it is easier than switching between
     #Edges and tuples
@@ -110,8 +124,8 @@ struct LRPlanarity{T}
     S::Stack{ConflictPair{T}} #Stack of tuples of Edges
     stack_bottom::Dict{Edge{T},ConflictPair{T}} #Dict of Tuples of Edges, indexed by Edge
     lowpt_edge::Dict{Edge{T},Edge{T}} #Dict of Edges, indexed by Edge 
-    left_ref::Dict{T,Edge{T}} #Dict of Edges, indexed by node
-    right_ref::Dict{T,Edge{T}} #Dict of Edges, indexed by node
+    #left_ref::Dict{T,Edge{T}} #Dict of Edges, indexed by node
+    #right_ref::Dict{T,Edge{T}} #Dict of Edges, indexed by node
     # skip embedding for now 
 end
 
@@ -122,11 +136,7 @@ function LRPlanarity(g)
     # copy G without adding self-loops
     output_g = SimpleGraph{T}()
     add_vertices_from!(g, output_g)
-    for e in edges(g)
-        if e.src != e.dst
-            add_edge!(output_g, e)
-        end
-    end
+    add_simple_edges_from!(g, output_g)
 
     N = nv(output_g)
 
@@ -157,8 +167,8 @@ function LRPlanarity(g)
     S = Stack{ConflictPair{T}}()
     stack_bottom = Dict{Edge{T},ConflictPair{T}}()
     lowpt_edge = Dict{Edge{T},Edge{T}}()
-    left_ref = Dict{T,Edge{T}}()
-    right_ref = Dict{T,Edge{T}}()
+    #left_ref = Dict{T,Edge{T}}()
+    #right_ref = Dict{T,Edge{T}}()
 
     #self.embedding = PlanarEmbedding()
     return LRPlanarity(
@@ -177,8 +187,8 @@ function LRPlanarity(g)
         S,
         stack_bottom,
         lowpt_edge,
-        left_ref,
-        right_ref,
+        #left_ref,
+        #right_ref,
     )
 end
 
@@ -196,8 +206,8 @@ function lowest(self::ConflictPair, planarity_state::LRPlanarity)
 end
 
 function lr_planarity!(self::LRPlanarity{T}) where {T}
-    V = nv(self.G)
-    E = ne(self.G)
+    V::Int64 = nv(self.G)
+    E::Int64 = ne(self.G)
 
     if V > 2 && (E > (3V - 6))
         # graph is not planar
@@ -208,6 +218,7 @@ function lr_planarity!(self::LRPlanarity{T}) where {T}
     for v in 1:nv(self.G) #for all vertices in G,
         self.adjs[v] = neighbors(self.G, v) ##neighbourhood of v
     end
+    # TODO this could be done during the alloc phase
 
     # orientation of the graph by depth first search traversal
     for v in vertices(self.G)
diff --git a/src/spanningtrees/pmfg.jl b/src/spanningtrees/pmfg.jl
index 9b4517240..2c57c9fa8 100644
--- a/src/spanningtrees/pmfg.jl
+++ b/src/spanningtrees/pmfg.jl
@@ -1,7 +1,7 @@
 """
-	```pmfg(g, dmtx = weights(g)```
+	pmfg(g)
 
-Compete the Planar Maximally Filtered Graph (PMFG) of `g`, given an optional set of weights `dmtx` that form a distance matrix.
+Compete the Planar Maximally Filtered Graph (PMFG) of weighted graph `g`.
 
 ### Examples
 ```
@@ -15,20 +15,59 @@ p_g = pmfg(g)
 ### References 
 - Tuminello et al. 2005
 """
-function pmfg(g, dmtx=weights(g))
-    T = eltype(g)
-    N = nv(g)
-    out = SimpleGraph{T}(N)
-    #create list of edges w/ weights
+function better_mst(
+    g::AG, distmx::AbstractMatrix{T}=weights(g); minimize=true
+) where {T<:Real,U,AG<:AbstractGraph{U}}
+    #check graph is not directed
+    if is_directed(g)
+        error("PMFG only supports non-directed graphs")
+    end
+
+    #construct a list of edge weights
+    weights = Vector{T}()
+    sizehint!(weights, ne(g))
     edge_list = collect(edges(g))
-    sort!(edge_list; by=x -> x.weight) #leaves biggest weight last
+    for e in edge_list
+        push!(weights, distmx[src(e), dst(e)])
+    end
+    #sort the set of edges by weight
+    #we try to maximally filter the graph and assume that weights
+    #represent distances. Thus we want to add edges with the 
+    #smallest distances first. Given that we pop the edge list, 
+    #we want the smallest weight edges at the end of the edge list
+    #after sorting, which means reversing the usual direction of
+    #the sort.
+    edge_list .= edge_list[sortperm(weights, rev = !minimize)]
+
+    #construct an initial graph
+    test_graph = SimpleGraph(nv(g))
+
+    #go through the edge list 
     while !isempty(edge_list)
-        new_edge_weighted = pop!(edge_list)
-        new_edge = Edge(new_edge_weighted.src, new_edge_weighted.dst)
-        add_edge!(out, new_edge)
-        if !is_planar(out)
-            rem_edge!(out, new_edge)
+        e = pop!(edge_list) #get most weighted edge
+        add_edge!(test_graph, e.src, e.dst) #add it to graph
+        if !is_planar(test_graph) #if resulting graph is not planar, remove it again
+            rem_edge!(test_graph, e.src, e.dst)
         end
+        (ne(test_graph) >= 3*nv(test_graph) - 6) && break #break if limit reached
     end
-    return out
+
+    return test_graph
 end
+
+#= 
+This could be improved a lot by not reallocating
+the LRP construct. 
+Things to reset on each planarity retest:
+heights 
+lowpts(2)
+nesting_depth 
+parent_edge
+DG (shame...)
+adjs (could just be re-edited?)
+ordered_adjs (same)
+Ref
+side 
+S 
+stack_bottom 
+lowpt_edge =#

From 2fae5322e164d90aa16252ef1c3654c7aa4b9606 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Mon, 2 Jan 2023 21:44:31 +0000
Subject: [PATCH 07/35] Typo in PMFG

---
 src/spanningtrees/pmfg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/spanningtrees/pmfg.jl b/src/spanningtrees/pmfg.jl
index 2c57c9fa8..4b6434855 100644
--- a/src/spanningtrees/pmfg.jl
+++ b/src/spanningtrees/pmfg.jl
@@ -15,7 +15,7 @@ p_g = pmfg(g)
 ### References 
 - Tuminello et al. 2005
 """
-function better_mst(
+function pmfg(
     g::AG, distmx::AbstractMatrix{T}=weights(g); minimize=true
 ) where {T<:Real,U,AG<:AbstractGraph{U}}
     #check graph is not directed

From 5b1cf8484f56b76f568501925de8a1494d287996 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Mon, 2 Jan 2023 22:49:49 +0000
Subject: [PATCH 08/35] Typos

---
 src/planarity.jl          | 4 ++--
 src/spanningtrees/pmfg.jl | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index f91e8956f..d159e5797 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -2,7 +2,7 @@
 # Algorithm from https://www.uni-konstanz.de/algo/publications/b-lrpt-sub.pdf
 # The implementation is heavily influenced by the recursive implementation in Networkx (https://networkx.org/documentation/stable/_modules/networkx/algorithms/planarity.html)
 
-import DataStructures: DefaultDict, Stack 
+import DataStructures: DefaultDict, Stack
 import Base: isempty
 
 """
@@ -105,7 +105,7 @@ function isempty(p::ConflictPair)
     return isempty(p.L) && isempty(p.R)
 end
 
-struct LRPlanarity{T<:Integer} 
+struct LRPlanarity{T<:Integer}
     #State class for the planarity test 
     #We index by Edge structs throughout as it is easier than switching between
     #Edges and tuples
diff --git a/src/spanningtrees/pmfg.jl b/src/spanningtrees/pmfg.jl
index 4b6434855..a7852c753 100644
--- a/src/spanningtrees/pmfg.jl
+++ b/src/spanningtrees/pmfg.jl
@@ -37,7 +37,7 @@ function pmfg(
     #we want the smallest weight edges at the end of the edge list
     #after sorting, which means reversing the usual direction of
     #the sort.
-    edge_list .= edge_list[sortperm(weights, rev = !minimize)]
+    edge_list .= edge_list[sortperm(weights; rev=!minimize)]
 
     #construct an initial graph
     test_graph = SimpleGraph(nv(g))
@@ -49,7 +49,7 @@ function pmfg(
         if !is_planar(test_graph) #if resulting graph is not planar, remove it again
             rem_edge!(test_graph, e.src, e.dst)
         end
-        (ne(test_graph) >= 3*nv(test_graph) - 6) && break #break if limit reached
+        (ne(test_graph) >= 3 * nv(test_graph) - 6) && break #break if limit reached
     end
 
     return test_graph

From 37547c7200b9f0b6c44a3de2e11e1f23e895bdd1 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Fri, 6 Jan 2023 10:25:22 +0000
Subject: [PATCH 09/35] Reduced memory usage of planarity test and PMFG
 alogorithm

---
 src/planarity.jl          | 37 ++++++++++++++++++++++++-------------
 src/spanningtrees/pmfg.jl |  2 +-
 2 files changed, 25 insertions(+), 14 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index d159e5797..5e3bef285 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -22,7 +22,7 @@ end
 
 #if g is directed, analyse the undirected version 
 function is_planar(dg::SimpleDiGraph)
-    g = SimpleGraph(dg)
+    g = SimpleGraph(dg) #TODO - do we need a separate method?
     return is_planar(g)
 end
 
@@ -109,7 +109,9 @@ struct LRPlanarity{T<:Integer}
     #State class for the planarity test 
     #We index by Edge structs throughout as it is easier than switching between
     #Edges and tuples
-    G::SimpleGraph{T} #Copy of the input graph
+    #G::SimpleGraph{T} #Copy of the input graph
+    V::Int64 
+    E::Int64 
     roots::Vector{T} #Vector of roots for disconnected graphs. Normally size(roots, 1) == 1
     height::DefaultDict{T,Int64} #DefaultDict of heights <: Int, indexed by node. default is -1
     lowpt::Dict{Edge{T},Int64} #Dict of low points, indexed by Edge
@@ -131,14 +133,17 @@ end
 
 #outer constructor for LRPlanarity
 function LRPlanarity(g)
+    V = Int64(nv(g)) #needs promoting
+    E = Int64(ne(g)) #JIC
     #record nodetype of g
     T = eltype(g)
+    #= 
     # copy G without adding self-loops
     output_g = SimpleGraph{T}()
     add_vertices_from!(g, output_g)
-    add_simple_edges_from!(g, output_g)
+    add_simple_edges_from!(g, output_g) =#
 
-    N = nv(output_g)
+    N = nv(g)
 
     roots = T[]
 
@@ -154,10 +159,14 @@ function LRPlanarity(g)
     parent_edge = DefaultDict{T,Edge{T}}(Edge(zero(T), zero(T)))
 
     # oriented DFS graph
-    DG = SimpleDiGraph{T}()
-    add_vertices_from!(g, DG)
+    DG = SimpleDiGraph{T}(N)
 
     adjs = Dict{T,Vector{T}}()
+    # make adjacency lists for dfs
+    for v in 1:nv(g) #for all vertices in G,
+        adjs[v] = neighbors(g, v) ##neighbourhood of v
+    end
+    
     ordered_adjs = Dict{T,Vector{T}}()
 
     ref = DefaultDict{Edge{T},Edge{T}}(empty_edge(T))
@@ -172,7 +181,9 @@ function LRPlanarity(g)
 
     #self.embedding = PlanarEmbedding()
     return LRPlanarity(
-        g,
+        #g,
+        V, 
+        E,
         roots,
         height,
         lowpt,
@@ -206,22 +217,22 @@ function lowest(self::ConflictPair, planarity_state::LRPlanarity)
 end
 
 function lr_planarity!(self::LRPlanarity{T}) where {T}
-    V::Int64 = nv(self.G)
-    E::Int64 = ne(self.G)
+    V = self.V
+    E = self.E
 
     if V > 2 && (E > (3V - 6))
         # graph is not planar
         return false
     end
 
-    # make adjacency lists for dfs
+    #= # make adjacency lists for dfs
     for v in 1:nv(self.G) #for all vertices in G,
         self.adjs[v] = neighbors(self.G, v) ##neighbourhood of v
     end
-    # TODO this could be done during the alloc phase
+    # TODO this could be done during the alloc phase =#
 
     # orientation of the graph by depth first search traversal
-    for v in vertices(self.G)
+    for v in 1:V
         if self.height[v] == -one(T) #using -1 rather than nothing for type stability. 
             self.height[v] = zero(T)
             push!(self.roots, v)
@@ -231,7 +242,7 @@ function lr_planarity!(self::LRPlanarity{T}) where {T}
 
     #Testing stage
     #First, sort the ordered_adjs by nesting depth
-    for v in 1:nv(self.G) #for all vertices in G,
+    for v in 1:V #for all vertices in G,
         #get neighboring nodes
         neighboring_nodes = T[]
         neighboring_nesting_depths = Int64[]
diff --git a/src/spanningtrees/pmfg.jl b/src/spanningtrees/pmfg.jl
index a7852c753..3b25a736d 100644
--- a/src/spanningtrees/pmfg.jl
+++ b/src/spanningtrees/pmfg.jl
@@ -37,7 +37,7 @@ function pmfg(
     #we want the smallest weight edges at the end of the edge list
     #after sorting, which means reversing the usual direction of
     #the sort.
-    edge_list .= edge_list[sortperm(weights; rev=!minimize)]
+    edge_list .= edge_list[sortperm(weights; rev=minimize)]
 
     #construct an initial graph
     test_graph = SimpleGraph(nv(g))

From 9937ef550edc461a893663e71fe6c468da1a88de Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Fri, 6 Jan 2023 10:35:03 +0000
Subject: [PATCH 10/35] Add tests for PMFG

---
 test/planarity.jl | 37 ++++++++++++++++++++++++++++++++-----
 1 file changed, 32 insertions(+), 5 deletions(-)

diff --git a/test/planarity.jl b/test/planarity.jl
index 958c3be0c..ac3baef50 100644
--- a/test/planarity.jl
+++ b/test/planarity.jl
@@ -46,13 +46,13 @@
         self = Graphs.LRPlanarity(dfs_g)
         #want to test dfs orientation 
         # make adjacency lists for dfs
-        for v in 1:nv(self.G) #for all vertices in G,
-            self.adjs[v] = neighbors(self.G, v) ##neighbourhood of v
+        for v in 1:self.V #for all vertices in G,
+            self.adjs[v] = neighbors(dfs_g, v) ##neighbourhood of v
         end
-        T = eltype(self.G)
+        T = eltype(dfs_g)
 
         # orientation of the graph by depth first search traversal
-        for v in vertices(self.G)
+        for v in 1:self.V
             if self.height[v] == -one(T) #using -1 rather than nothing for type stability. 
                 self.height[v] = zero(T)
                 push!(self.roots, v)
@@ -120,4 +120,31 @@
 
         @test is_planar(dg) == true
     end
-end
+
+    @testset "PMFG tests" begin 
+        k5 = complete_graph(5)
+        k5_am = Matrix{Float64}(adjacency_matrix(k5))
+
+        #random weights
+        for i in CartesianIndices(k5_am)
+            if k5_am[i] == 1
+                k5_am[i] = rand()
+            end
+        end
+
+        #let's make 1->5 very distant 
+        k5_am[1, 5] = 10.
+
+        k5_am .= Symmetric(k5_am)
+        k5 = SimpleWeightedGraph(k5_am)
+
+        #correct result of PMFG
+        correct_am = [0  1  1  1  0
+        1  0  1  1  1
+        1  1  0  1  1
+        1  1  1  0  1
+        0  1  1  1  0]
+
+        @test correct_am == adjacency_matrix(pmfg(k5))
+    end
+end
\ No newline at end of file

From d83dea3ecf28928543098f5aa071fac3bc8ae4e9 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Fri, 6 Jan 2023 19:24:52 +0000
Subject: [PATCH 11/35] Fixed trait function

---
 src/spanningtrees/pmfg.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/spanningtrees/pmfg.jl b/src/spanningtrees/pmfg.jl
index 3b25a736d..02c592775 100644
--- a/src/spanningtrees/pmfg.jl
+++ b/src/spanningtrees/pmfg.jl
@@ -15,8 +15,9 @@ p_g = pmfg(g)
 ### References 
 - Tuminello et al. 2005
 """
-function pmfg(
-    g::AG, distmx::AbstractMatrix{T}=weights(g); minimize=true
+function pmfg end
+@traitfn function pmfg(
+    g::AG::(!IsDirected), distmx::AbstractMatrix{T}=weights(g); minimize=true
 ) where {T<:Real,U,AG<:AbstractGraph{U}}
     #check graph is not directed
     if is_directed(g)

From 1a1316adfe29dda6d88e54bfa86860ea5092e484 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Sat, 7 Jan 2023 12:08:15 +0000
Subject: [PATCH 12/35] Fix test failures, format, and refactor tests into
 correct folder

---
 src/spanningtrees/pmfg.jl  |  2 +-
 test/planarity.jl          | 33 +++------------------------------
 test/runtests.jl           |  1 +
 test/spanningtrees/pmfg.jl | 26 ++++++++++++++++++++++++++
 4 files changed, 31 insertions(+), 31 deletions(-)
 create mode 100644 test/spanningtrees/pmfg.jl

diff --git a/src/spanningtrees/pmfg.jl b/src/spanningtrees/pmfg.jl
index 02c592775..e7b15470b 100644
--- a/src/spanningtrees/pmfg.jl
+++ b/src/spanningtrees/pmfg.jl
@@ -17,7 +17,7 @@ p_g = pmfg(g)
 """
 function pmfg end
 @traitfn function pmfg(
-    g::AG::(!IsDirected), distmx::AbstractMatrix{T}=weights(g); minimize=true
+    g::AG::(!IsDirected); distmx::AbstractMatrix{T}=weights(g), minimize=true
 ) where {T<:Real,U,AG<:AbstractGraph{U}}
     #check graph is not directed
     if is_directed(g)
diff --git a/test/planarity.jl b/test/planarity.jl
index ac3baef50..7bd0189ea 100644
--- a/test/planarity.jl
+++ b/test/planarity.jl
@@ -46,13 +46,13 @@
         self = Graphs.LRPlanarity(dfs_g)
         #want to test dfs orientation 
         # make adjacency lists for dfs
-        for v in 1:self.V #for all vertices in G,
+        for v in 1:(self.V) #for all vertices in G,
             self.adjs[v] = neighbors(dfs_g, v) ##neighbourhood of v
         end
         T = eltype(dfs_g)
 
         # orientation of the graph by depth first search traversal
-        for v in 1:self.V
+        for v in 1:(self.V)
             if self.height[v] == -one(T) #using -1 rather than nothing for type stability. 
                 self.height[v] = zero(T)
                 push!(self.roots, v)
@@ -120,31 +120,4 @@
 
         @test is_planar(dg) == true
     end
-
-    @testset "PMFG tests" begin 
-        k5 = complete_graph(5)
-        k5_am = Matrix{Float64}(adjacency_matrix(k5))
-
-        #random weights
-        for i in CartesianIndices(k5_am)
-            if k5_am[i] == 1
-                k5_am[i] = rand()
-            end
-        end
-
-        #let's make 1->5 very distant 
-        k5_am[1, 5] = 10.
-
-        k5_am .= Symmetric(k5_am)
-        k5 = SimpleWeightedGraph(k5_am)
-
-        #correct result of PMFG
-        correct_am = [0  1  1  1  0
-        1  0  1  1  1
-        1  1  0  1  1
-        1  1  1  0  1
-        0  1  1  1  0]
-
-        @test correct_am == adjacency_matrix(pmfg(k5))
-    end
-end
\ No newline at end of file
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 666b5ca58..59ab10422 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -104,6 +104,7 @@ tests = [
     "spanningtrees/boruvka",
     "spanningtrees/kruskal",
     "spanningtrees/prim",
+    "spanningtrees/pmfg.jl",
     "steinertree/steiner_tree",
     "biconnectivity/articulation",
     "biconnectivity/biconnect",
diff --git a/test/spanningtrees/pmfg.jl b/test/spanningtrees/pmfg.jl
new file mode 100644
index 000000000..8a0b1081a
--- /dev/null
+++ b/test/spanningtrees/pmfg.jl
@@ -0,0 +1,26 @@
+@testset "PMFG tests" begin
+    k5 = complete_graph(5)
+    k5_am = Matrix{Float64}(adjacency_matrix(k5))
+
+    #random weights
+    for i in CartesianIndices(k5_am)
+        if k5_am[i] == 1
+            k5_am[i] = rand()
+        end
+    end
+
+    #let's make 1->5 very distant 
+    k5_am[1, 5] = 10.0
+    k5_am[5, 1] = 10.0
+
+    #correct result of PMFG
+    correct_am = [
+        0 1 1 1 0
+        1 0 1 1 1
+        1 1 0 1 1
+        1 1 1 0 1
+        0 1 1 1 0
+    ]
+
+    @test correct_am == Matrix(adjacency_matrix(pmfg(k5; distmx=k5_am)))
+end

From eaa91802fdc3d50fc1a6fa840d2251e3682e1465 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Sat, 7 Jan 2023 12:10:51 +0000
Subject: [PATCH 13/35] Change pmfg to planar_maximally_filtered_graph

---
 src/Graphs.jl              | 4 ++--
 src/spanningtrees/pmfg.jl  | 8 ++++----
 test/runtests.jl           | 2 +-
 test/spanningtrees/pmfg.jl | 2 +-
 4 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/Graphs.jl b/src/Graphs.jl
index 66b95fab3..7068ee61a 100644
--- a/src/Graphs.jl
+++ b/src/Graphs.jl
@@ -417,7 +417,7 @@ export
 
     # planarity
     is_planar,
-    pmfg
+    planar_maximally_filtered_graph
 """
     Graphs
 
@@ -532,7 +532,7 @@ include("vertexcover/random_vertex_cover.jl")
 include("Experimental/Experimental.jl")
 include("Parallel/Parallel.jl")
 include("planarity.jl")
-include("spanningtrees/pmfg.jl")
+include("spanningtrees/planar_maximally_filtered_graph.jl")
 
 using .LinAlg
 end # module
diff --git a/src/spanningtrees/pmfg.jl b/src/spanningtrees/pmfg.jl
index e7b15470b..017d4d1e8 100644
--- a/src/spanningtrees/pmfg.jl
+++ b/src/spanningtrees/pmfg.jl
@@ -1,5 +1,5 @@
 """
-	pmfg(g)
+	planar_maximally_filtered_graph(g)
 
 Compete the Planar Maximally Filtered Graph (PMFG) of weighted graph `g`.
 
@@ -9,14 +9,14 @@ using Graphs, SimpleWeightedGraphs
 N = 20
 M = Symmetric(randn(N, N))
 g = SimpleWeightedGraph(M) 
-p_g = pmfg(g)
+p_g = planar_maximally_filtered_graph(g)
 ```
 
 ### References 
 - Tuminello et al. 2005
 """
-function pmfg end
-@traitfn function pmfg(
+function planar_maximally_filtered_graph end
+@traitfn function planar_maximally_filtered_graph(
     g::AG::(!IsDirected); distmx::AbstractMatrix{T}=weights(g), minimize=true
 ) where {T<:Real,U,AG<:AbstractGraph{U}}
     #check graph is not directed
diff --git a/test/runtests.jl b/test/runtests.jl
index 59ab10422..5713f3243 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -104,7 +104,7 @@ tests = [
     "spanningtrees/boruvka",
     "spanningtrees/kruskal",
     "spanningtrees/prim",
-    "spanningtrees/pmfg.jl",
+    "spanningtrees/planar_maximally_filtered_graph.jl",
     "steinertree/steiner_tree",
     "biconnectivity/articulation",
     "biconnectivity/biconnect",
diff --git a/test/spanningtrees/pmfg.jl b/test/spanningtrees/pmfg.jl
index 8a0b1081a..c92e66a13 100644
--- a/test/spanningtrees/pmfg.jl
+++ b/test/spanningtrees/pmfg.jl
@@ -22,5 +22,5 @@
         0 1 1 1 0
     ]
 
-    @test correct_am == Matrix(adjacency_matrix(pmfg(k5; distmx=k5_am)))
+    @test correct_am == Matrix(adjacency_matrix(planar_maximally_filtered_graph(k5; distmx=k5_am)))
 end

From dae5ea2a877f8b6e53e3040c9bd4786b338ba1eb Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Sat, 7 Jan 2023 12:12:16 +0000
Subject: [PATCH 14/35] Typo in tests

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 5713f3243..bc45b0e18 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -104,7 +104,7 @@ tests = [
     "spanningtrees/boruvka",
     "spanningtrees/kruskal",
     "spanningtrees/prim",
-    "spanningtrees/planar_maximally_filtered_graph.jl",
+    "spanningtrees/planar_maximally_filtered_graph",
     "steinertree/steiner_tree",
     "biconnectivity/articulation",
     "biconnectivity/biconnect",

From 40064846b63c02718b02d1798f48ba18fbba5f73 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Sat, 7 Jan 2023 12:14:25 +0000
Subject: [PATCH 15/35] Fix file names

---
 src/spanningtrees/{pmfg.jl => planar_maximally_filtered_graph.jl} | 0
 .../spanningtrees/{pmfg.jl => planar_maximally_filtered_graph.jl} | 0
 2 files changed, 0 insertions(+), 0 deletions(-)
 rename src/spanningtrees/{pmfg.jl => planar_maximally_filtered_graph.jl} (100%)
 rename test/spanningtrees/{pmfg.jl => planar_maximally_filtered_graph.jl} (100%)

diff --git a/src/spanningtrees/pmfg.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
similarity index 100%
rename from src/spanningtrees/pmfg.jl
rename to src/spanningtrees/planar_maximally_filtered_graph.jl
diff --git a/test/spanningtrees/pmfg.jl b/test/spanningtrees/planar_maximally_filtered_graph.jl
similarity index 100%
rename from test/spanningtrees/pmfg.jl
rename to test/spanningtrees/planar_maximally_filtered_graph.jl

From 2dcc0d79c924fb7091eb0a60202f6965f73546d2 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <76974735+josephcbradley@users.noreply.github.com>
Date: Sat, 7 Jan 2023 12:24:42 +0000
Subject: [PATCH 16/35] Fix formatting issues

---
 src/planarity.jl                                      | 8 ++++----
 src/spanningtrees/planar_maximally_filtered_graph.jl  | 2 +-
 test/spanningtrees/planar_maximally_filtered_graph.jl | 2 +-
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index 5e3bef285..b9056cea1 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -110,8 +110,8 @@ struct LRPlanarity{T<:Integer}
     #We index by Edge structs throughout as it is easier than switching between
     #Edges and tuples
     #G::SimpleGraph{T} #Copy of the input graph
-    V::Int64 
-    E::Int64 
+    V::Int64
+    E::Int64
     roots::Vector{T} #Vector of roots for disconnected graphs. Normally size(roots, 1) == 1
     height::DefaultDict{T,Int64} #DefaultDict of heights <: Int, indexed by node. default is -1
     lowpt::Dict{Edge{T},Int64} #Dict of low points, indexed by Edge
@@ -166,7 +166,7 @@ function LRPlanarity(g)
     for v in 1:nv(g) #for all vertices in G,
         adjs[v] = neighbors(g, v) ##neighbourhood of v
     end
-    
+
     ordered_adjs = Dict{T,Vector{T}}()
 
     ref = DefaultDict{Edge{T},Edge{T}}(empty_edge(T))
@@ -182,7 +182,7 @@ function LRPlanarity(g)
     #self.embedding = PlanarEmbedding()
     return LRPlanarity(
         #g,
-        V, 
+        V,
         E,
         roots,
         height,
diff --git a/src/spanningtrees/planar_maximally_filtered_graph.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
index 017d4d1e8..c72c1b2e6 100644
--- a/src/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/src/spanningtrees/planar_maximally_filtered_graph.jl
@@ -17,7 +17,7 @@ p_g = planar_maximally_filtered_graph(g)
 """
 function planar_maximally_filtered_graph end
 @traitfn function planar_maximally_filtered_graph(
-    g::AG::(!IsDirected); distmx::AbstractMatrix{T}=weights(g), minimize=true
+    g::AG::(!IsDirected), distmx::AbstractMatrix{T}=weights(g); minimize=true
 ) where {T<:Real,U,AG<:AbstractGraph{U}}
     #check graph is not directed
     if is_directed(g)
diff --git a/test/spanningtrees/planar_maximally_filtered_graph.jl b/test/spanningtrees/planar_maximally_filtered_graph.jl
index c92e66a13..cb56df172 100644
--- a/test/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/test/spanningtrees/planar_maximally_filtered_graph.jl
@@ -22,5 +22,5 @@
         0 1 1 1 0
     ]
 
-    @test correct_am == Matrix(adjacency_matrix(planar_maximally_filtered_graph(k5; distmx=k5_am)))
+    @test correct_am == Matrix(adjacency_matrix(planar_maximally_filtered_graph(k5, k5_am)))
 end

From 739374e3a10413de4362e9a7d384b34dec3fc2fd Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sat, 7 Jan 2023 15:30:53 +0000
Subject: [PATCH 17/35] Remove redundant type check

---
 src/spanningtrees/planar_maximally_filtered_graph.jl | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/spanningtrees/planar_maximally_filtered_graph.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
index c72c1b2e6..71741ff1b 100644
--- a/src/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/src/spanningtrees/planar_maximally_filtered_graph.jl
@@ -19,10 +19,6 @@ function planar_maximally_filtered_graph end
 @traitfn function planar_maximally_filtered_graph(
     g::AG::(!IsDirected), distmx::AbstractMatrix{T}=weights(g); minimize=true
 ) where {T<:Real,U,AG<:AbstractGraph{U}}
-    #check graph is not directed
-    if is_directed(g)
-        error("PMFG only supports non-directed graphs")
-    end
 
     #construct a list of edge weights
     weights = Vector{T}()

From 274619ee7dcbbf93206aa2ef8b9d4dc16337a25f Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sat, 7 Jan 2023 15:32:01 +0000
Subject: [PATCH 18/35] Add link to paper for PMFG

---
 src/spanningtrees/planar_maximally_filtered_graph.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/spanningtrees/planar_maximally_filtered_graph.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
index 71741ff1b..315abaa2e 100644
--- a/src/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/src/spanningtrees/planar_maximally_filtered_graph.jl
@@ -13,7 +13,7 @@ p_g = planar_maximally_filtered_graph(g)
 ```
 
 ### References 
-- Tuminello et al. 2005
+- [Tuminello et al. 2005](https://doi.org/10.1073/pnas.0500298102)
 """
 function planar_maximally_filtered_graph end
 @traitfn function planar_maximally_filtered_graph(

From 6ed143f8415e0196d8b9551b5d28adbf6d0f2477 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sat, 7 Jan 2023 15:50:42 +0000
Subject: [PATCH 19/35] Document and enforce return type

---
 src/spanningtrees/planar_maximally_filtered_graph.jl  | 4 ++--
 test/spanningtrees/planar_maximally_filtered_graph.jl | 7 +++++++
 2 files changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/spanningtrees/planar_maximally_filtered_graph.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
index 315abaa2e..316b8b7c4 100644
--- a/src/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/src/spanningtrees/planar_maximally_filtered_graph.jl
@@ -1,7 +1,7 @@
 """
 	planar_maximally_filtered_graph(g)
 
-Compete the Planar Maximally Filtered Graph (PMFG) of weighted graph `g`.
+Compete the Planar Maximally Filtered Graph (PMFG) of weighted graph `g`. Returns a `SimpleGraph{eltype(g)}`.
 
 ### Examples
 ```
@@ -37,7 +37,7 @@ function planar_maximally_filtered_graph end
     edge_list .= edge_list[sortperm(weights; rev=minimize)]
 
     #construct an initial graph
-    test_graph = SimpleGraph(nv(g))
+    test_graph = SimpleGraph{U}(nv(g))
 
     #go through the edge list 
     while !isempty(edge_list)
diff --git a/test/spanningtrees/planar_maximally_filtered_graph.jl b/test/spanningtrees/planar_maximally_filtered_graph.jl
index cb56df172..9a97cdd70 100644
--- a/test/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/test/spanningtrees/planar_maximally_filtered_graph.jl
@@ -23,4 +23,11 @@
     ]
 
     @test correct_am == Matrix(adjacency_matrix(planar_maximally_filtered_graph(k5, k5_am)))
+
+    #type test 
+    N = 10
+    g = SimpleGraph{Int16}(N)
+    X = rand(N, N); C = X'*X
+    p = planar_maximally_filtered_graph(g, C)
+    @test typeof(p) == SimpleGraph{eltype(g)}
 end

From 7b30b014605f194f47d038eaa6be3844126e173c Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sat, 7 Jan 2023 15:53:16 +0000
Subject: [PATCH 20/35] Formatting

---
 test/spanningtrees/planar_maximally_filtered_graph.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/spanningtrees/planar_maximally_filtered_graph.jl b/test/spanningtrees/planar_maximally_filtered_graph.jl
index 9a97cdd70..e8dce4a9f 100644
--- a/test/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/test/spanningtrees/planar_maximally_filtered_graph.jl
@@ -27,7 +27,8 @@
     #type test 
     N = 10
     g = SimpleGraph{Int16}(N)
-    X = rand(N, N); C = X'*X
+    X = rand(N, N)
+    C = X' * X
     p = planar_maximally_filtered_graph(g, C)
     @test typeof(p) == SimpleGraph{eltype(g)}
 end

From bc0d3fabab47f6a516bb35d2c08f8b558d3a1aaf Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sat, 7 Jan 2023 16:19:02 +0000
Subject: [PATCH 21/35] `is_planar` can now dispatch directly on directed
 graphs

---
 src/planarity.jl | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index b9056cea1..0fe97cafe 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -20,12 +20,6 @@ function is_planar(g)
     return lr_planarity!(lrp)
 end
 
-#if g is directed, analyse the undirected version 
-function is_planar(dg::SimpleDiGraph)
-    g = SimpleGraph(dg) #TODO - do we need a separate method?
-    return is_planar(g)
-end
-
 #aux functions 
 function add_vertices_from!(input_g, output_g)
     for v in vertices(input_g)
@@ -164,7 +158,7 @@ function LRPlanarity(g)
     adjs = Dict{T,Vector{T}}()
     # make adjacency lists for dfs
     for v in 1:nv(g) #for all vertices in G,
-        adjs[v] = neighbors(g, v) ##neighbourhood of v
+        adjs[v] = all_neighbors(g, v) ##neighbourhood of v
     end
 
     ordered_adjs = Dict{T,Vector{T}}()

From dd37b65f41dab8ff4111d8677a8bcc27974d1e32 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sat, 7 Jan 2023 16:20:17 +0000
Subject: [PATCH 22/35] Type annotation of state constructor

---
 src/planarity.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index 0fe97cafe..983166409 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -126,7 +126,7 @@ struct LRPlanarity{T<:Integer}
 end
 
 #outer constructor for LRPlanarity
-function LRPlanarity(g)
+function LRPlanarity(g::AG) where AG<:AbstractGraph
     V = Int64(nv(g)) #needs promoting
     E = Int64(ne(g)) #JIC
     #record nodetype of g

From e1b5e4434c9d0a6860603012ab9d05ad88374f2b Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sat, 7 Jan 2023 16:55:36 +0000
Subject: [PATCH 23/35] Remove trait and refactor according to suggestions

---
 src/planarity.jl                                |  2 +-
 .../planar_maximally_filtered_graph.jl          | 17 ++++++-----------
 2 files changed, 7 insertions(+), 12 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index 983166409..a6c9fd206 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -126,7 +126,7 @@ struct LRPlanarity{T<:Integer}
 end
 
 #outer constructor for LRPlanarity
-function LRPlanarity(g::AG) where AG<:AbstractGraph
+function LRPlanarity(g::AG) where {AG<:AbstractGraph}
     V = Int64(nv(g)) #needs promoting
     E = Int64(ne(g)) #JIC
     #record nodetype of g
diff --git a/src/spanningtrees/planar_maximally_filtered_graph.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
index 316b8b7c4..528d171c7 100644
--- a/src/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/src/spanningtrees/planar_maximally_filtered_graph.jl
@@ -15,18 +15,14 @@ p_g = planar_maximally_filtered_graph(g)
 ### References 
 - [Tuminello et al. 2005](https://doi.org/10.1073/pnas.0500298102)
 """
-function planar_maximally_filtered_graph end
-@traitfn function planar_maximally_filtered_graph(
-    g::AG::(!IsDirected), distmx::AbstractMatrix{T}=weights(g); minimize=true
+
+function planar_maximally_filtered_graph(
+    g::AG, distmx::AbstractMatrix{T}=weights(g); minimize=true
 ) where {T<:Real,U,AG<:AbstractGraph{U}}
 
     #construct a list of edge weights
-    weights = Vector{T}()
-    sizehint!(weights, ne(g))
     edge_list = collect(edges(g))
-    for e in edge_list
-        push!(weights, distmx[src(e), dst(e)])
-    end
+    weights = [distmx[src(e), dst(e)] for e in edge_list]
     #sort the set of edges by weight
     #we try to maximally filter the graph and assume that weights
     #represent distances. Thus we want to add edges with the 
@@ -34,14 +30,13 @@ function planar_maximally_filtered_graph end
     #we want the smallest weight edges at the end of the edge list
     #after sorting, which means reversing the usual direction of
     #the sort.
-    edge_list .= edge_list[sortperm(weights; rev=minimize)]
+    edge_list .= edge_list[sortperm(weights; rev=!minimize)]
 
     #construct an initial graph
     test_graph = SimpleGraph{U}(nv(g))
 
     #go through the edge list 
-    while !isempty(edge_list)
-        e = pop!(edge_list) #get most weighted edge
+    for e in edge_list
         add_edge!(test_graph, e.src, e.dst) #add it to graph
         if !is_planar(test_graph) #if resulting graph is not planar, remove it again
             rem_edge!(test_graph, e.src, e.dst)

From f557d0069784875ae8ce447931cfc44be1a95fd5 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sat, 7 Jan 2023 17:10:19 +0000
Subject: [PATCH 24/35] Add docs to makefile

---
 docs/src/algorithms/spanningtrees.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/src/algorithms/spanningtrees.md b/docs/src/algorithms/spanningtrees.md
index d037abd25..4a0cb54a4 100644
--- a/docs/src/algorithms/spanningtrees.md
+++ b/docs/src/algorithms/spanningtrees.md
@@ -16,6 +16,8 @@ Pages = [
     "spanningtrees/boruvka.jl",
     "spanningtrees/kruskal.jl",
     "spanningtrees/prim.jl",
+    "planarity.jl",
+    "spanningtrees/planar_maximally_filtered_graph.jl",
 ]
 
 ```

From e059c216aa46c1a5402ac5f0ce777ec5db6b7e4c Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sun, 8 Jan 2023 16:28:04 +0000
Subject: [PATCH 25/35] Typo in docs

---
 src/spanningtrees/planar_maximally_filtered_graph.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/spanningtrees/planar_maximally_filtered_graph.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
index 528d171c7..06107b82d 100644
--- a/src/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/src/spanningtrees/planar_maximally_filtered_graph.jl
@@ -1,7 +1,7 @@
 """
 	planar_maximally_filtered_graph(g)
 
-Compete the Planar Maximally Filtered Graph (PMFG) of weighted graph `g`. Returns a `SimpleGraph{eltype(g)}`.
+Compute the Planar Maximally Filtered Graph (PMFG) of weighted graph `g`. Returns a `SimpleGraph{eltype(g)}`.
 
 ### Examples
 ```

From 7a2899c611a511cba5f5598b05818690ffe05f5f Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sun, 8 Jan 2023 16:29:02 +0000
Subject: [PATCH 26/35] Trim unnecessary Matrix call

---
 test/spanningtrees/planar_maximally_filtered_graph.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/spanningtrees/planar_maximally_filtered_graph.jl b/test/spanningtrees/planar_maximally_filtered_graph.jl
index e8dce4a9f..813ec5aa9 100644
--- a/test/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/test/spanningtrees/planar_maximally_filtered_graph.jl
@@ -22,7 +22,7 @@
         0 1 1 1 0
     ]
 
-    @test correct_am == Matrix(adjacency_matrix(planar_maximally_filtered_graph(k5, k5_am)))
+    @test correct_am == adjacency_matrix(planar_maximally_filtered_graph(k5, k5_am))
 
     #type test 
     N = 10

From e2b486226fa67b0ab80c52ff2a0232d91f4b680f Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Tue, 10 Jan 2023 12:52:02 +0000
Subject: [PATCH 27/35] Add further test of PMFG

---
 .../planar_maximally_filtered_graph.jl           | 16 +++++++++++++++-
 1 file changed, 15 insertions(+), 1 deletion(-)

diff --git a/test/spanningtrees/planar_maximally_filtered_graph.jl b/test/spanningtrees/planar_maximally_filtered_graph.jl
index 813ec5aa9..a2caf387e 100644
--- a/test/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/test/spanningtrees/planar_maximally_filtered_graph.jl
@@ -31,4 +31,18 @@
     C = X' * X
     p = planar_maximally_filtered_graph(g, C)
     @test typeof(p) == SimpleGraph{eltype(g)}
-end
+
+    #Test that MST is a subset of the PMFG 
+    N = 50
+    X = rand(N, N); D = X'*X
+    c = complete_graph(N)
+    p = planar_maximally_filtered_graph(c, D)
+    mst_edges = kruskal_mst(c, D)
+    is_subgraph = true 
+    for mst_edge in mst_edges
+        if mst_edge ∉ edges(p)
+            is_subgraph = false
+        end
+    end
+    @test is_subgraph
+end
\ No newline at end of file

From 61c1a68edcc8e565fbdb37bab8493030e426552a Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Tue, 10 Jan 2023 15:45:27 +0000
Subject: [PATCH 28/35] Fix for crash on higher order graphs

---
 src/planarity.jl  | 35 +++++++++++++++++++++++++++++++++--
 test/planarity.jl |  9 ++++++++-
 2 files changed, 41 insertions(+), 3 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index a6c9fd206..367665df2 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -99,6 +99,37 @@ function isempty(p::ConflictPair)
     return isempty(p.L) && isempty(p.R)
 end
 
+# ATM (Julia 1.8.4, DataStructures v0.18.13), DefaultDict crashes
+# for large order matrices when we attempt to trim the back edges. 
+#To fix this we create a manual version of the DefaultDict that seems to be more stable
+
+struct ManualDict{A, B}
+    d::Dict{A, B}
+    default::B
+end
+
+function ManualDict(A, B, default)
+    ManualDict(Dict{A, B}(), default)
+end
+
+import Base: getindex
+
+function getindex(md::ManualDict, x)
+    d = md.d
+    if haskey(d, x)
+        d[x]
+    else
+        d[x] = md.default
+        md.default
+    end
+end
+
+function setindex!(md::ManualDict, X, key)
+    setindex!(md.d, X, key)
+end
+
+
+
 struct LRPlanarity{T<:Integer}
     #State class for the planarity test 
     #We index by Edge structs throughout as it is easier than switching between
@@ -115,7 +146,7 @@ struct LRPlanarity{T<:Integer}
     DG::SimpleDiGraph{T} #Directed graph for the orientation phase
     adjs::Dict{T,Vector{T}} #Dict of neighbors of nodes, indexed by node
     ordered_adjs::Dict{T,Vector{T}} #Dict of neighbors of nodes sorted by nesting depth, indexed by node
-    ref::DefaultDict{Edge{T},Edge{T}} #DefaultDict of Edges, indexed by Edge
+    ref::ManualDict{Edge{T},Edge{T}} #ManualDict of Edges, indexed by Edge
     side::DefaultDict{Edge{T},Int8} #DefaultDict of +/- 1, indexed by edge
     S::Stack{ConflictPair{T}} #Stack of tuples of Edges
     stack_bottom::Dict{Edge{T},ConflictPair{T}} #Dict of Tuples of Edges, indexed by Edge
@@ -163,7 +194,7 @@ function LRPlanarity(g::AG) where {AG<:AbstractGraph}
 
     ordered_adjs = Dict{T,Vector{T}}()
 
-    ref = DefaultDict{Edge{T},Edge{T}}(empty_edge(T))
+    ref = ManualDict(Edge{T},Edge{T}, empty_edge(T))
     side = DefaultDict{Edge{T},Int8}(one(Int8))
 
     # stack of conflict pairs
diff --git a/test/planarity.jl b/test/planarity.jl
index 7bd0189ea..e79c2c2df 100644
--- a/test/planarity.jl
+++ b/test/planarity.jl
@@ -120,4 +120,11 @@
 
         @test is_planar(dg) == true
     end
-end
+
+    @testset "ManualDict tests" begin 
+        md = Graphs.ManualDict(Int64, Int64, 0)
+        @test md[0] == 0
+        md[0] = 1
+        @test md[0] == 1
+    end
+end
\ No newline at end of file

From 18c9a30c73edc15b6ba3020e5b388ac533e1c415 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Tue, 10 Jan 2023 15:48:40 +0000
Subject: [PATCH 29/35] Tidy and format

---
 src/planarity.jl                              | 51 ++-----------------
 test/planarity.jl                             | 13 +----
 .../planar_maximally_filtered_graph.jl        |  7 +--
 3 files changed, 11 insertions(+), 60 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index 367665df2..ff7bd2a65 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -20,33 +20,6 @@ function is_planar(g)
     return lr_planarity!(lrp)
 end
 
-#aux functions 
-function add_vertices_from!(input_g, output_g)
-    for v in vertices(input_g)
-        add_vertex!(output_g)
-    end
-end
-
-function add_edges_from!(input_g, output_g)
-    for e in edges(input_g)
-        add_edge!(output_g, e)
-    end
-end
-
-function add_simple_edges_from!(input_g, output_g)
-    for e in edges(input_g)
-        add_edge!(output_g, e.dst, e.src)
-    end
-end
-
-#= function add__simple_edges_from!(input_g::SimpleWeightedGraph, output_g)
-    for e in edges(input_g)
-        if e.dst != e.src
-            add_edge!(output_g, e.src, e.dst)
-        end
-    end
-end =#
-
 #Simple structs to be used in algorithm. Keep private for now. 
 function empty_edge(T)
     return Edge{T}(0, 0)
@@ -103,13 +76,13 @@ end
 # for large order matrices when we attempt to trim the back edges. 
 #To fix this we create a manual version of the DefaultDict that seems to be more stable
 
-struct ManualDict{A, B}
-    d::Dict{A, B}
+struct ManualDict{A,B}
+    d::Dict{A,B}
     default::B
 end
 
 function ManualDict(A, B, default)
-    ManualDict(Dict{A, B}(), default)
+    return ManualDict(Dict{A,B}(), default)
 end
 
 import Base: getindex
@@ -125,11 +98,9 @@ function getindex(md::ManualDict, x)
 end
 
 function setindex!(md::ManualDict, X, key)
-    setindex!(md.d, X, key)
+    return setindex!(md.d, X, key)
 end
 
-
-
 struct LRPlanarity{T<:Integer}
     #State class for the planarity test 
     #We index by Edge structs throughout as it is easier than switching between
@@ -162,12 +133,6 @@ function LRPlanarity(g::AG) where {AG<:AbstractGraph}
     E = Int64(ne(g)) #JIC
     #record nodetype of g
     T = eltype(g)
-    #= 
-    # copy G without adding self-loops
-    output_g = SimpleGraph{T}()
-    add_vertices_from!(g, output_g)
-    add_simple_edges_from!(g, output_g) =#
-
     N = nv(g)
 
     roots = T[]
@@ -194,7 +159,7 @@ function LRPlanarity(g::AG) where {AG<:AbstractGraph}
 
     ordered_adjs = Dict{T,Vector{T}}()
 
-    ref = ManualDict(Edge{T},Edge{T}, empty_edge(T))
+    ref = ManualDict(Edge{T}, Edge{T}, empty_edge(T))
     side = DefaultDict{Edge{T},Int8}(one(Int8))
 
     # stack of conflict pairs
@@ -250,12 +215,6 @@ function lr_planarity!(self::LRPlanarity{T}) where {T}
         return false
     end
 
-    #= # make adjacency lists for dfs
-    for v in 1:nv(self.G) #for all vertices in G,
-        self.adjs[v] = neighbors(self.G, v) ##neighbourhood of v
-    end
-    # TODO this could be done during the alloc phase =#
-
     # orientation of the graph by depth first search traversal
     for v in 1:V
         if self.height[v] == -one(T) #using -1 rather than nothing for type stability. 
diff --git a/test/planarity.jl b/test/planarity.jl
index e79c2c2df..ab3571cc3 100644
--- a/test/planarity.jl
+++ b/test/planarity.jl
@@ -1,13 +1,4 @@
 @testset "Planarity tests" begin
-    @testset "Aux functions tests" begin
-        g = complete_graph(10)
-        f = SimpleGraph()
-        Graphs.add_vertices_from!(g, f)
-        @test nv(g) == nv(f)
-        Graphs.add_edges_from!(g, f)
-        @test g == f
-    end
-
     @testset "LRP constructor tests" begin
         function lrp_constructor_test(g)
             try
@@ -121,10 +112,10 @@
         @test is_planar(dg) == true
     end
 
-    @testset "ManualDict tests" begin 
+    @testset "ManualDict tests" begin
         md = Graphs.ManualDict(Int64, Int64, 0)
         @test md[0] == 0
         md[0] = 1
         @test md[0] == 1
     end
-end
\ No newline at end of file
+end
diff --git a/test/spanningtrees/planar_maximally_filtered_graph.jl b/test/spanningtrees/planar_maximally_filtered_graph.jl
index a2caf387e..72c9a3220 100644
--- a/test/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/test/spanningtrees/planar_maximally_filtered_graph.jl
@@ -34,15 +34,16 @@
 
     #Test that MST is a subset of the PMFG 
     N = 50
-    X = rand(N, N); D = X'*X
+    X = rand(N, N)
+    D = X' * X
     c = complete_graph(N)
     p = planar_maximally_filtered_graph(c, D)
     mst_edges = kruskal_mst(c, D)
-    is_subgraph = true 
+    is_subgraph = true
     for mst_edge in mst_edges
         if mst_edge ∉ edges(p)
             is_subgraph = false
         end
     end
     @test is_subgraph
-end
\ No newline at end of file
+end

From 2d1f0c0d9857ac3a85a04c91c106428883a2f55c Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Wed, 8 Feb 2023 16:10:49 +0000
Subject: [PATCH 30/35] Make src and dst generic

---
 src/planarity.jl                                     | 6 +++---
 src/spanningtrees/planar_maximally_filtered_graph.jl | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index ff7bd2a65..f8209602f 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -26,7 +26,7 @@ function empty_edge(T)
 end
 
 function isempty(e::Edge{T}) where {T}
-    return e.src == zero(T) && e.dst == zero(T)
+    return src(e) == zero(T) && dst(e) == zero(T)
 end
 
 mutable struct Interval{T}
@@ -341,7 +341,7 @@ function dfs_testing!(self, v)
 
     #remove back edges returning to parent 
     if !isempty(e)#v is not root
-        u = e.src
+        u = src(e)
         #trim edges ending at parent u, algo 5 
         trim_back_edges!(self, u)
         #side of e is side of highest returning edge 
@@ -410,7 +410,7 @@ function edge_constraints!(self, ei, e)
     return true
 end
 
-function trim_back_edges!(self, u)
+@noinline function trim_back_edges!(self, u)
     #trim back edges ending at u 
     #drop entire conflict pairs 
     while !isempty(self.S) && (lowest(first(self.S), self) == self.height[u])
diff --git a/src/spanningtrees/planar_maximally_filtered_graph.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
index 06107b82d..f36a1471e 100644
--- a/src/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/src/spanningtrees/planar_maximally_filtered_graph.jl
@@ -37,9 +37,9 @@ function planar_maximally_filtered_graph(
 
     #go through the edge list 
     for e in edge_list
-        add_edge!(test_graph, e.src, e.dst) #add it to graph
+        add_edge!(test_graph, src(e), dst(e)) #add it to graph
         if !is_planar(test_graph) #if resulting graph is not planar, remove it again
-            rem_edge!(test_graph, e.src, e.dst)
+            rem_edge!(test_graph, src(e), dst(e))
         end
         (ne(test_graph) >= 3 * nv(test_graph) - 6) && break #break if limit reached
     end

From 89154306fb0482c9f6e5c6278a0bf0f8be341ca6 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sun, 26 Feb 2023 14:29:54 +0000
Subject: [PATCH 31/35] Add first six edges without testing

---
 .../planar_maximally_filtered_graph.jl         | 18 ++++++++++++++++--
 1 file changed, 16 insertions(+), 2 deletions(-)

diff --git a/src/spanningtrees/planar_maximally_filtered_graph.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
index f36a1471e..dfb7f90e0 100644
--- a/src/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/src/spanningtrees/planar_maximally_filtered_graph.jl
@@ -20,6 +20,15 @@ function planar_maximally_filtered_graph(
     g::AG, distmx::AbstractMatrix{T}=weights(g); minimize=true
 ) where {T<:Real,U,AG<:AbstractGraph{U}}
 
+    #if graph has <= 6 edges, just return it
+    if ne(g) <= 6
+        test_graph = SimpleGraph{U}(nv(g))
+        for e in edges(g)
+            add_edge!(g, e)
+        end
+        return test_graph
+    end
+
     #construct a list of edge weights
     edge_list = collect(edges(g))
     weights = [distmx[src(e), dst(e)] for e in edge_list]
@@ -35,8 +44,13 @@ function planar_maximally_filtered_graph(
     #construct an initial graph
     test_graph = SimpleGraph{U}(nv(g))
 
-    #go through the edge list 
-    for e in edge_list
+    for e in edge_list[1:6]
+        #we can always add the first six edges of a graph 
+        add_edge!(test_graph, src(e), dst(e))
+    end
+
+    #go through the rest of the edge list 
+    for e in edge_list[7:end]
         add_edge!(test_graph, src(e), dst(e)) #add it to graph
         if !is_planar(test_graph) #if resulting graph is not planar, remove it again
             rem_edge!(test_graph, src(e), dst(e))

From 25de62877060776a35356c3d8cc734b30d787a40 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sun, 26 Feb 2023 15:42:37 +0000
Subject: [PATCH 32/35] Consistent syntax

---
 src/planarity.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index f8209602f..a4f7d6174 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -146,7 +146,7 @@ function LRPlanarity(g::AG) where {AG<:AbstractGraph}
 
     # None == Edge(0, 0)  for our type-stable algo
 
-    parent_edge = DefaultDict{T,Edge{T}}(Edge(zero(T), zero(T)))
+    parent_edge = DefaultDict{T,Edge{T}}(empty_edge(T))
 
     # oriented DFS graph
     DG = SimpleDiGraph{T}(N)

From 1b093234dcf1f8d42209dcd6f5496ce223a80f35 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sun, 26 Feb 2023 16:09:05 +0000
Subject: [PATCH 33/35] Revert to Default Dict

---
 src/planarity.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index a4f7d6174..10aac3eac 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -117,7 +117,7 @@ struct LRPlanarity{T<:Integer}
     DG::SimpleDiGraph{T} #Directed graph for the orientation phase
     adjs::Dict{T,Vector{T}} #Dict of neighbors of nodes, indexed by node
     ordered_adjs::Dict{T,Vector{T}} #Dict of neighbors of nodes sorted by nesting depth, indexed by node
-    ref::ManualDict{Edge{T},Edge{T}} #ManualDict of Edges, indexed by Edge
+    ref::DefaultDict{Edge{T},Edge{T}} #ManualDict of Edges, indexed by Edge
     side::DefaultDict{Edge{T},Int8} #DefaultDict of +/- 1, indexed by edge
     S::Stack{ConflictPair{T}} #Stack of tuples of Edges
     stack_bottom::Dict{Edge{T},ConflictPair{T}} #Dict of Tuples of Edges, indexed by Edge
@@ -159,7 +159,7 @@ function LRPlanarity(g::AG) where {AG<:AbstractGraph}
 
     ordered_adjs = Dict{T,Vector{T}}()
 
-    ref = ManualDict(Edge{T}, Edge{T}, empty_edge(T))
+    ref = DefaultDict{Edge{T}, Edge{T}}(empty_edge(T))
     side = DefaultDict{Edge{T},Int8}(one(Int8))
 
     # stack of conflict pairs

From 3344f3fff7299ed1965a4f13d57623632e7ac710 Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sun, 26 Feb 2023 16:09:27 +0000
Subject: [PATCH 34/35] Add reset lrp state function

---
 src/planarity.jl | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/src/planarity.jl b/src/planarity.jl
index 10aac3eac..80fa9dfd0 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -193,6 +193,38 @@ function LRPlanarity(g::AG) where {AG<:AbstractGraph}
     )
 end
 
+function lrp_type(lrp::LRPlanarity{T}) where T
+    T
+end
+
+function reset_lrp_state!(lrp_state)
+    T = lrp_type(lrp_state)
+    #resets the LRP state 
+    #reset heights
+    for k ∈ keys(lrp_state.height)
+        lrp_state.height[k] = -1
+    end
+
+    for k in keys(lrp_state.parent_edge)
+        lrp_state.parent_edge[k] = empty_edge(T)
+    end
+
+    for e in edges(lrp_state.DG)
+        rem_edge!(lrp_state.DG, e)
+    end
+
+    for k ∈ keys(lrp_state.ref)
+        lrp_state.ref[k] = empty_edge(T)
+    end
+
+    for k ∈ keys(lrp_state.side)
+        lrp_state.side[k] = one(Int8)
+    end
+
+    empty!(lrp_state.S)
+end
+
+
 function lowest(self::ConflictPair, planarity_state::LRPlanarity)
     #Returns the lowest lowpoint of a conflict pair
     if isempty(self.L)

From 9ea7ec80e240663205ab55a64d590306cc41d22a Mon Sep 17 00:00:00 2001
From: Joseph Bradley <josephbradley16@googlemail.com>
Date: Sun, 26 Feb 2023 20:17:46 +0000
Subject: [PATCH 35/35] Reset planarity struct to speed up PMFG

---
 src/planarity.jl                              | 21 ++++++++++++++++---
 .../planar_maximally_filtered_graph.jl        | 14 ++++++-------
 2 files changed, 24 insertions(+), 11 deletions(-)

diff --git a/src/planarity.jl b/src/planarity.jl
index 80fa9dfd0..0c437059d 100644
--- a/src/planarity.jl
+++ b/src/planarity.jl
@@ -101,7 +101,7 @@ function setindex!(md::ManualDict, X, key)
     return setindex!(md.d, X, key)
 end
 
-struct LRPlanarity{T<:Integer}
+mutable struct LRPlanarity{T<:Integer}
     #State class for the planarity test 
     #We index by Edge structs throughout as it is easier than switching between
     #Edges and tuples
@@ -197,9 +197,20 @@ function lrp_type(lrp::LRPlanarity{T}) where T
     T
 end
 
-function reset_lrp_state!(lrp_state)
+function reset_lrp_state!(lrp_state, g)
     T = lrp_type(lrp_state)
     #resets the LRP state 
+
+
+    #reset roots 
+    empty!(lrp_state.roots)
+
+    #reset lowpts 
+    #empty!(lrp_state.lowpt)
+    #empty!(lrp_state.lowpt2)
+    #reset nesting depth
+    empty!(lrp_state.nesting_depth)
+
     #reset heights
     for k ∈ keys(lrp_state.height)
         lrp_state.height[k] = -1
@@ -213,6 +224,10 @@ function reset_lrp_state!(lrp_state)
         rem_edge!(lrp_state.DG, e)
     end
 
+    for v in 1:nv(g) #for all vertices in G,
+        lrp_state.adjs[v] = all_neighbors(g, v) ##neighbourhood of v
+    end
+
     for k ∈ keys(lrp_state.ref)
         lrp_state.ref[k] = empty_edge(T)
     end
@@ -442,7 +457,7 @@ function edge_constraints!(self, ei, e)
     return true
 end
 
-@noinline function trim_back_edges!(self, u)
+function trim_back_edges!(self, u)
     #trim back edges ending at u 
     #drop entire conflict pairs 
     while !isempty(self.S) && (lowest(first(self.S), self) == self.height[u])
diff --git a/src/spanningtrees/planar_maximally_filtered_graph.jl b/src/spanningtrees/planar_maximally_filtered_graph.jl
index dfb7f90e0..e3beb0ee1 100644
--- a/src/spanningtrees/planar_maximally_filtered_graph.jl
+++ b/src/spanningtrees/planar_maximally_filtered_graph.jl
@@ -49,10 +49,13 @@ function planar_maximally_filtered_graph(
         add_edge!(test_graph, src(e), dst(e))
     end
 
+    #generate lrp state
+    lrp_state = LRPlanarity(test_graph)
     #go through the rest of the edge list 
     for e in edge_list[7:end]
         add_edge!(test_graph, src(e), dst(e)) #add it to graph
-        if !is_planar(test_graph) #if resulting graph is not planar, remove it again
+        reset_lrp_state!(lrp_state, test_graph)
+        if !lr_planarity!(lrp_state) #if resulting graph is not planar, remove it again
             rem_edge!(test_graph, src(e), dst(e))
         end
         (ne(test_graph) >= 3 * nv(test_graph) - 6) && break #break if limit reached
@@ -66,14 +69,9 @@ This could be improved a lot by not reallocating
 the LRP construct. 
 Things to reset on each planarity retest:
 heights 
-lowpts(2)
-nesting_depth 
 parent_edge
 DG (shame...)
-adjs (could just be re-edited?)
-ordered_adjs (same)
 Ref
 side 
-S 
-stack_bottom 
-lowpt_edge =#
+S  =#
+