iterative version
diff --git a/src/isomorphism.rs b/src/isomorphism.rs
index 7eecd85..0d307b9 100644
--- a/src/isomorphism.rs
+++ b/src/isomorphism.rs
@@ -17,7 +17,7 @@
 struct Vf2State<Ty, Ix> {
     /// The current mapping M(s) of nodes from G0 → G1 and G1 → G0,
     /// NodeIndex::end() for no mapping.
-    mapping: Vec<NodeIndex<Ix>>,
+    pub mapping: Vec<NodeIndex<Ix>>,
     /// out[i] is non-zero if i is in either M_0(s) or Tout_0(s)
     /// These are all the next vertices that are not mapped yet, but
     /// have an outgoing edge from the mapping.
@@ -173,7 +173,12 @@
     }
 
     let mut st = [Vf2State::new(g0), Vf2State::new(g1)];
-    try_match(&mut st, g0, g1, &mut NoSemanticMatch, &mut NoSemanticMatch).unwrap_or(false)
+    //let y = try_match(&mut st, g0, g1, &mut NoSemanticMatch, &mut NoSemanticMatch).unwrap_or(false);
+    //st = [Vf2State::new(g0), Vf2State::new(g1)];
+    let x = try_match_iter(&mut st, g0, g1, &mut NoSemanticMatch, &mut NoSemanticMatch).unwrap_or(false);
+    //println!("{:?}", (&x, &y));
+    //debug_assert!(x == y);
+    return x;
 }
 
 /// [Graph] Return `true` if the graphs `g0` and `g1` are isomorphic.
@@ -221,6 +226,279 @@
 }
 
 /// Return Some(bool) if isomorphism is decided, else None.
+fn try_match_iter<N, E, Ty, Ix, F, G>(mut st: &mut [Vf2State<Ty, Ix>; 2],
+                                 g0: &Graph<N, E, Ty, Ix>,
+                                 g1: &Graph<N, E, Ty, Ix>,
+                                 node_match: &mut F,
+                                 edge_match: &mut G)
+    -> Option<bool>
+    where Ty: EdgeType,
+          Ix: IndexType,
+          F: SemanticMatcher<N>,
+          G: SemanticMatcher<E>,
+{
+    if st[0].is_complete() {
+        return Some(true);
+    }
+    let g = [g0, g1];
+    let graph_indices = 0..2;
+    let end = NodeIndex::end();
+
+    // A "depth first" search of a valid mapping from graph 1 to graph 2
+
+    // F(s, n, m) -- evaluate state s and add mapping n <-> m
+
+    // Find least T1out node (in st.out[1] but not in M[1])
+    #[derive(Copy, Clone, PartialEq, Debug)]
+    enum OpenList {
+        Out,
+        In,
+        Other,
+    }
+
+    #[derive(Clone, PartialEq, Debug)]
+    enum Frame<N> {
+        Outer,
+        Inner{ nodes: [N; 2], open_list: OpenList },
+        Unwind{ nodes: [N; 2], open_list: OpenList },
+    }
+
+    let next_candidate = |st: &mut[Vf2State<Ty, Ix>; 2]|
+        -> Option<(NodeIndex<Ix>, NodeIndex<Ix>, OpenList)> {
+        let mut to_index;
+        let mut from_index = None;
+        let mut open_list = OpenList::Out;
+        // Try the out list
+        to_index = st[1].next_out_index(0);
+
+        if to_index.is_some() {
+            from_index = st[0].next_out_index(0);
+            open_list = OpenList::Out;
+        }
+        // Try the in list
+        if to_index.is_none() || from_index.is_none() {
+            to_index = st[1].next_in_index(0);
+
+            if to_index.is_some() {
+                from_index = st[0].next_in_index(0);
+                open_list = OpenList::In;
+            }
+        }
+        // Try the other list -- disconnected graph
+        if to_index.is_none() || from_index.is_none() {
+            to_index = st[1].next_rest_index(0);
+            if to_index.is_some() {
+                from_index = st[0].next_rest_index(0);
+                open_list = OpenList::Other;
+            }
+        }
+        match (from_index, to_index) {
+            (Some(n), Some(m)) => Some((NodeIndex::new(n), NodeIndex::new(m), open_list)),
+            // No more candidates
+            _ => None,
+        }
+    };
+    let next_from_ix = |st: &mut[Vf2State<Ty, Ix>; 2], nx: NodeIndex<Ix>, open_list: OpenList| -> Option<NodeIndex<Ix>> {
+        // Find the next node index to try on the `from` side of the mapping
+        let start = nx.index() + 1;
+        let cand0 = match open_list {
+            OpenList::Out => st[0].next_out_index(start),
+            OpenList::In => st[0].next_in_index(start),
+            OpenList::Other => st[0].next_rest_index(start),
+        }.map(|c| c + start); // compensate for start offset.
+        match cand0 {
+            None => None, // no more candidates
+            Some(ix) => {
+                debug_assert!(ix >= start);
+                Some(NodeIndex::new(ix))
+            },
+        }
+    };
+    //fn pop_state(nodes: [NodeIndex<Ix>; 2]) {
+    let pop_state = |st: &mut[Vf2State<Ty, Ix>; 2], nodes: [NodeIndex<Ix>; 2]| {
+        // Restore state.
+        for j in graph_indices.clone() {
+            st[j].pop_mapping(nodes[j], g[j]);
+        }
+    };
+    //fn push_state(nodes: [NodeIndex<Ix>; 2]) {
+    let push_state = |st: &mut[Vf2State<Ty, Ix>; 2],nodes: [NodeIndex<Ix>; 2]| {
+        // Add mapping nx <-> mx to the state
+        for j in graph_indices.clone() {
+            st[j].push_mapping(nodes[j], nodes[1-j], g[j]);
+        }
+    };
+    //fn is_feasible(nodes: [NodeIndex<Ix>; 2]) -> bool {
+    let mut is_feasible = |st: &mut[Vf2State<Ty, Ix>; 2], nodes: [NodeIndex<Ix>; 2]| -> bool {
+        // Check syntactic feasibility of mapping by ensuring adjacencies
+        // of nx map to adjacencies of mx.
+        //
+        // nx == map to => mx
+        //
+        // R_succ
+        //
+        // Check that every neighbor of nx is mapped to a neighbor of mx,
+        // then check the reverse, from mx to nx. Check that they have the same
+        // count of edges.
+        //
+        // Note: We want to check the lookahead measures here if we can,
+        // R_out: Equal for G0, G1: Card(Succ(G, n) ^ Tout); for both Succ and Pred
+        // R_in: Same with Tin
+        // R_new: Equal for G0, G1: Ñ n Pred(G, n); both Succ and Pred,
+        //      Ñ is G0 - M - Tin - Tout
+        // last attempt to add these did not speed up any of the testcases
+        let mut succ_count = [0, 0];
+        for j in graph_indices.clone() {
+            for n_neigh in g[j].neighbors(nodes[j]) {
+                succ_count[j] += 1;
+                // handle the self loop case; it's not in the mapping (yet)
+                let m_neigh = if nodes[j] != n_neigh {
+                    st[j].mapping[n_neigh.index()]
+                } else {
+                    nodes[1 - j]
+                };
+                if m_neigh == end {
+                    continue;
+                }
+                let has_edge = g[1-j].is_adjacent(&st[1-j].adjacency_matrix, nodes[1-j], m_neigh);
+                if !has_edge {
+                    return false;
+                }
+            }
+        }
+        if succ_count[0] != succ_count[1] {
+            return false;
+        }
+        // R_pred
+        if g[0].is_directed() {
+            let mut pred_count = [0, 0];
+            for j in graph_indices.clone() {
+                for n_neigh in g[j].neighbors_directed(nodes[j], Incoming) {
+                    pred_count[j] += 1;
+                    // the self loop case is handled in outgoing
+                    let m_neigh = st[j].mapping[n_neigh.index()];
+                    if m_neigh == end {
+                        continue;
+                    }
+                    let has_edge = g[1-j].is_adjacent(&st[1-j].adjacency_matrix, m_neigh, nodes[1-j]);
+                    if !has_edge {
+                        return false;
+                    }
+                }
+            }
+            if pred_count[0] != pred_count[1] {
+                return false;
+            }
+        }
+        // semantic feasibility: compare associated data for nodes
+        if F::enabled() && !node_match.eq(&g[0][nodes[0]], &g[1][nodes[1]]) {
+            return false;
+        }
+        // semantic feasibility: compare associated data for edges
+        if G::enabled() {
+            // outgoing edges
+            for j in graph_indices.clone() {
+                let mut edges = g[j].neighbors(nodes[j]).detach();
+                while let Some((n_edge, n_neigh)) = edges.next(g[j]) {
+                    // handle the self loop case; it's not in the mapping (yet)
+                    let m_neigh = if nodes[j] != n_neigh {
+                        st[j].mapping[n_neigh.index()]
+                    } else {
+                        nodes[1 - j]
+                    };
+                    if m_neigh == end {
+                        continue;
+                    }
+                    match g[1-j].find_edge(nodes[1 - j], m_neigh) {
+                        Some(m_edge) => {
+                            if !edge_match.eq(&g[j][n_edge], &g[1-j][m_edge]) {
+                                return false;
+                            }
+                        }
+                        None => unreachable!() // covered by syntactic check
+                    }
+                }
+            }
+            // incoming edges
+            if g[0].is_directed() {
+                for j in graph_indices.clone() {
+                    let mut edges = g[j].neighbors_directed(nodes[j], Incoming).detach();
+                    while let Some((n_edge, n_neigh)) = edges.next(g[j]) {
+                        // the self loop case is handled in outgoing
+                        let m_neigh = st[j].mapping[n_neigh.index()];
+                        if m_neigh == end {
+                            continue;
+                        }
+                        match g[1-j].find_edge(m_neigh, nodes[1-j]) {
+                            Some(m_edge) => {
+                                if !edge_match.eq(&g[j][n_edge], &g[1-j][m_edge]) {
+                                    return false;
+                                }
+                            }
+                            None => unreachable!() // covered by syntactic check
+                        }
+                    }
+                }
+            }
+        }
+        true
+    };
+    let mut stack = vec![Frame::Outer];
+    println!("readysetgo! {:?}", &stack);
+
+    while let Some(frame) = stack.pop() {
+        println!("\t{:?}", &frame);
+        match frame {
+            Frame::Unwind{nodes, open_list: ol} => {
+                pop_state(&mut st, nodes);
+
+                match next_from_ix(&mut st, nodes[0], ol) {
+                    None => continue,
+                    Some(nx) => {
+                        let f = Frame::Inner{nodes: [nx, nodes[1]], open_list: ol};
+                        stack.push(f);
+                    }
+                }
+            }
+            Frame::Outer => {
+                match next_candidate(&mut st) {
+                    None => continue,
+                    Some((nx, mx, ol)) => {
+                        let f = Frame::Inner{nodes: [nx, mx], open_list: ol};
+                        stack.push(f); 
+                    }
+                }
+            }
+            Frame::Inner{nodes, open_list: ol} => {
+                if is_feasible(&mut st, nodes) {
+                    push_state(&mut st, nodes);
+                    if st[0].is_complete() {
+                        return Some(true);
+                    }
+                    // Check cardinalities of Tin, Tout sets
+                    if st[0].out_size == st[1].out_size &&
+                       st[0].ins_size == st[1].ins_size {
+                           let f0 = Frame::Unwind{nodes: nodes, open_list: ol};
+                           stack.push(f0);
+                           stack.push(Frame::Outer);
+                           continue;
+                    }
+                    pop_state(&mut st, nodes);
+                }
+                match next_from_ix(&mut st, nodes[0], ol) {
+                    None => continue,
+                    Some(nx) => {
+                        let f = Frame::Inner{nodes: [nx, nodes[1]], open_list: ol};
+                        stack.push(f);
+                    }
+                }
+            }
+        }
+    }
+    None
+}
+
+/// Return Some(bool) if isomorphism is decided, else None.
 fn try_match<N, E, Ty, Ix, F, G>(st: &mut [Vf2State<Ty, Ix>; 2],
                                  g0: &Graph<N, E, Ty, Ix>,
                                  g1: &Graph<N, E, Ty, Ix>,
@@ -291,6 +569,7 @@
 
     let mut nx = NodeIndex::new(cand0);
     let mx = NodeIndex::new(cand1);
+    println!("{:?}\n{:?}\n{:?}\n", (cand0, cand1), &st[0].mapping, &st[1].mapping);
 
     let mut first = true;
 
@@ -303,6 +582,7 @@
                 OpenList::In => st[0].next_in_index(start),
                 OpenList::Other => st[0].next_rest_index(start),
             }.map(|c| c + start); // compensate for start offset.
+            println!("  {:?}", (cand0, cand1));
             nx = match cand0 {
                 None => break, // no more candidates
                 Some(ix) => NodeIndex::new(ix),
@@ -428,6 +708,7 @@
                 }
             }
         }
+        println!("\tfeasible");
 
         // Add mapping nx <-> mx to the state
         for j in graph_indices.clone() {
@@ -438,13 +719,14 @@
         if st[0].out_size == st[1].out_size &&
            st[0].ins_size == st[1].ins_size
         {
-
+            println!("\trecurse");
             // Recurse
             match try_match(st, g0, g1, node_match, edge_match) {
                 None => {}
                 result => return result,
             }
         }
+        println!("\tpop");
 
         // Restore state.
         for j in graph_indices.clone() {