Repository: PermutaTriangle/comb_spec_searcher (language: python)

Issue:
Constructor equiv refactor



Discussion / hints:

Fix patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index ee3d6a7e1..3578ab4c9 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,5 +1,5 @@
 repos:
   - repo: https://github.com/psf/black
-    rev: 20.8b1
+    rev: 21.4b0
     hooks:
       - id: black
diff --git a/.zenodo.json b/.zenodo.json
new file mode 100644
index 000000000..1e64adbb3
--- /dev/null
+++ b/.zenodo.json
@@ -0,0 +1,25 @@
+{
+    "upload_type": "software",
+    "creators": [
+        {
+            "affiliation": "Reykjavik University",
+            "name": "Émile Nadeau"
+        },
+        {
+            "affiliation": "Reykjavik University",
+            "name": "Christian Bean"
+        },
+        {
+            "affiliation": "Reykjavik University",
+            "name": "Henning Ulfarsson"
+        },
+        {
+            "affiliation": "Reykjavik University",
+            "name": "Jon Steinn Eliasson"
+        },
+        {
+            "affiliation": "Marquette University",
+            "name": "Jay Pantone"
+        }
+    ]
+}
\ No newline at end of file
diff --git a/CHANGELOG.md b/CHANGELOG.md
index d5b1fa042..a63ae1c6b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,14 +4,66 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/)
 and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+## [4.0.0] - 2021-06-14
+### Added
+- Reverse rule default back to the equation of the original rule in case
+  `NotImplementedError`
+- `find_bijection_between` tries to find a bijection between classes given
+  a `CombinatorialSpecificationSearcher` object for both.
+- Added Forest searching capability to the css. Those are specification that can
+  you reverse rule that are not equivalences.
+- Special forward and backward maps, called indexed forward and backward maps.
+  They are to be used for bijections and their purposes is to support bijections
+  for non-injective forward maps by labelling the resulting objects or map from
+  an labelled object.
+- `NonBijectiveRule` will implement a labelling system for indexed forward
+  and backward maps.
+- `EqPathParallelSpecFinder` that, on top of the base class, validates any
+  potential path contained in equivalence labels.
+- `PartialSpecificationRuleExtractor` that extracts rules from partially built
+  specifications and two subclasses with specific applications of that.
+
+### Changed
+- If a rule in a specification cannot be sanity checked (e.g., counting is not
+  implemented), a warning is printed and sanity checking continues, instead of returning
+  the exception.
+- More responsibility have been delegated to the rule database. This will allow
+  for more flexibility in the implementation of the rule database. Things like
+  computing specification rule now happens at the rule db level.
+
+### Fixed
+- Removed a debug print
+- Sharing of a specification html via gofile API.
+- Moves local `Constructor.param_map` function outward so that specifications can be
+  pickled.
+- Fix bug in complement constructor counting
+- Fixing a bug in counting for equivalence rule from a reverse rule
+- Fixing in the formal step of reverse equivalence rule. The formal step now
+  state that the rule is reversed.
+- Removed bug in equation of `DisjointUnion` that ignores multiple CVs mapping to the
+  same in a child
+
+### Deprecated
+- Python 3.6 is no longer supported
+- `LimitedStrategyRuleDB` is not longer part of css
+
+
 ## [3.0.0] - 2021-01-04
 ### Added
-- Automatic bijection between equivalent specifications through the functions `get_bijection_to` and `are_isomorphic` of the specifications class. The bijection object holds a `map` function that performs the actual mapping.
+- Automatic bijection between equivalent specifications through the functions
+  `get_bijection_to` and `are_isomorphic` of the specifications class. The bijection
+  object holds a `map` function that performs the actual mapping.
 - Sanity check for random sampling on rule
 - Strategy must not define a `is_two_way` method in order to decide if they can
   be used to find the count of a children knowing the parent's and other
   children' count. If so the constructor returned by the new method
   `reverse_constructor` is used.
+- Adds `expand_all_verified_with_pack` to `Specification` to attempt expansion of
+  verified nodes with a given pack and time limit.
+- Adds `unexpanded_verified_classes` to `Specification` to return the set of verified
+  classes.
 
 ### Changed
 - Specification are now built using a set of rules
@@ -19,6 +71,8 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
 - Streamlined the extraction of a specification from a searcher.
 - All unary strategy are now store at the level of the equivalence database in
   order to avoid some productivity issue with catalytic variables.
+- Adds a `max_expansion_time` optional parameter to
+  `comb_spec_searcher._auto_search_rules`.
 
 ### Fixed
 - `forward_map` of `EquivalencePathRule`
diff --git a/README.rst b/README.rst
index ff634fddb..1117d16fc 100644
--- a/README.rst
+++ b/README.rst
@@ -13,14 +13,14 @@ Combinatorial Specification Searcher
     :target: https://pypi.python.org/pypi/comb_spec_searcher
 .. image:: https://img.shields.io/pypi/pyversions/comb_spec_searcher.svg
     :target: https://pypi.python.org/pypi/comb_spec_searcher
-
 .. image:: http://img.shields.io/badge/readme-tested-brightgreen.svg
     :alt: Travis
     :target: https://travis-ci.org/PermutaTriangle/comb_spec_searcher
-
 .. image:: https://requires.io/github/PermutaTriangle/comb_spec_searcher/requirements.svg?branch=master
      :target: https://requires.io/github/PermutaTriangle/comb_spec_searcher/requirements/?branch=master
      :alt: Requirements Status
+.. image:: https://zenodo.org/badge/121520109.svg
+   :target: https://zenodo.org/badge/latestdoi/121520109
 
 The ``comb_spec_searcher`` package contains code for combinatorial
 exploration.
@@ -78,8 +78,14 @@ using ``in``.
 
 .. code:: python
 
-   >>> w = "acbabcabbb"
-   >>> p = "abcab"
+   >>> from comb_spec_searcher import CombinatorialObject
+   
+   >>> class Word(str, CombinatorialObject):
+   ...     def size(self):
+   ...         return str.__len__(self)
+
+   >>> w = Word("acbabcabbb")
+   >>> p = Word("abcab")
    >>> p in w
    True
 
@@ -101,11 +107,11 @@ create a new python ``class`` representing this that inherits from
    >>> from comb_spec_searcher import CombinatorialClass
 
 
-   >>> class AvoidingWithPrefix(CombinatorialClass):
+   >>> class AvoidingWithPrefix(CombinatorialClass[Word]):
    ...     def __init__(self, prefix, patterns, alphabet, just_prefix=False):
    ...         self.alphabet = frozenset(alphabet)
-   ...         self.prefix = prefix
-   ...         self.patterns = frozenset(patterns)
+   ...         self.prefix = Word(prefix)
+   ...         self.patterns = frozenset(map(Word, patterns))
    ...         self.just_prefix = just_prefix # this will be needed later
 
 Inheriting from ``CombinatorialClass`` requires you to implement a few
@@ -233,7 +239,7 @@ create multiple rules, and as such should be implemented as generators.
    >>> from comb_spec_searcher import DisjointUnionStrategy
 
 
-   >>> class ExpansionStrategy(DisjointUnionStrategy):
+   >>> class ExpansionStrategy(DisjointUnionStrategy[AvoidingWithPrefix, Word]):
    ...     def decomposition_function(self, avoiding_with_prefix):
    ...        if not avoiding_with_prefix.just_prefix:
    ...           alphabet, prefix, patterns = (
@@ -265,7 +271,7 @@ create multiple rules, and as such should be implemented as generators.
    ...           if word[: len(child.prefix)] == child.prefix:
    ...                 return (
    ...                    tuple(None for _ in range(idx + 1))
-   ...                    + (child,)
+   ...                    + (word,)
    ...                    + tuple(None for _ in range(len(children) - idx - 1))
    ...                 )
    ...
@@ -292,7 +298,7 @@ constructor is cartesian product a ``DecompositionRule``.
    >>> from comb_spec_searcher import CartesianProductStrategy
 
 
-   >>> class RemoveFrontOfPrefix(CartesianProductStrategy):
+   >>> class RemoveFrontOfPrefix(CartesianProductStrategy[AvoidingWithPrefix, Word]):
    ...     def decomposition_function(self, avoiding_with_prefix):
    ...        """If the k is the maximum size of a pattern to be avoided, then any
    ...        occurrence using indices further to the right of the prefix can use at
@@ -318,7 +324,8 @@ constructor is cartesian product a ``DecompositionRule``.
    ...        )
    ...        # safe will be the index of the prefix in which we can remove upto without
    ...        # affecting the avoidance conditions
-   ...        safe = max(0, len(prefix) - max(len(p) for p in patterns) + 1)
+   ...        m = max(len(p) for p in patterns) if patterns else 1
+   ...        safe = max(0, len(prefix) - m + 1)
    ...        for i in range(safe, len(prefix)):
    ...           end = prefix[i:]
    ...           if any(end == patt[: len(end)] for patt in patterns):
@@ -340,7 +347,7 @@ constructor is cartesian product a ``DecompositionRule``.
    ...        if children is None:
    ...           children = self.decomposition_function(avoiding_with_prefix)
    ...           assert children is not None
-   ...        return Word(words[0] + words[1])
+   ...        yield Word(words[0] + words[1])
    ...
    ...     def forward_map(self, comb_class, word, children=None):
    ...        """
@@ -399,6 +406,7 @@ will).
 
    >>> searcher = CombinatorialSpecificationSearcher(start_class, pack)
    >>> spec = searcher.auto_search()
+   >>> # spec.show() will display the specification in your browser
 
 Now that we have a ``CombinatorialSpecification``, the obvious
 thing we want to do is find the generating function for the class that
@@ -422,11 +430,11 @@ of a given size in the class.
    ...            return
    ...         for letters in product(self.alphabet,
    ...                                 repeat=size - len(self.prefix)):
-   ...             yield self.prefix + "".join(a for a in letters)
+   ...             yield Word(self.prefix + "".join(a for a in letters))
    ...
    ...     if self.just_prefix:
    ...         if size == len(self.prefix) and not self.is_empty():
-   ...             yield self.prefix
+   ...             yield Word(self.prefix)
    ...         return
    ...     for word in possible_words():
    ...         if all(patt not in word for patt in self.patterns):
@@ -453,3 +461,69 @@ Moreover, we can get directly the number of objects by size with the method
 
 You can now try this yourself using the file ``example.py``, which can
 count any set of words avoiding consecutive patterns.
+
+Now we will demonstrate how a bijection can be found between classes.
+We will first need a couple of imports.
+
+.. code:: python
+
+   >>> from comb_spec_searcher.bijection import ParallelSpecFinder
+   >>> from comb_spec_searcher.isomorphism import Bijection
+
+We start by defining our two classes that we wish to find a bijection between.
+
+.. code:: python
+
+   >>> prefix1 = ''
+   >>> patterns1 = ["00"]
+   >>> alphabet1 = ['0', '1']
+   >>> class1 = AvoidingWithPrefix(prefix1, patterns1, alphabet1)
+   >>> prefix2 = ''
+   >>> patterns2 = ["bb"]
+   >>> alphabet2 = ['a', 'b']
+   >>> class2 = AvoidingWithPrefix(prefix2, patterns2, alphabet2)
+
+To find a bijection we expand the universe given a pack for both classes
+and try to construct specifications that are parallel.
+
+.. code:: python
+
+   >>> searcher1 = CombinatorialSpecificationSearcher(class1, pack)
+   >>> searcher2 = CombinatorialSpecificationSearcher(class2, pack)
+
+We get two parallel specs if successful, ``None`` otherwise.
+
+.. code:: python
+
+   >>> specs = ParallelSpecFinder(searcher1, searcher2).find()
+
+We then construct the bijection from the parallel specifications.
+
+.. code:: python
+
+   >>> bijection = Bijection.construct(*specs)
+
+We can use the `Bijection` object to map (either way) sampled objects
+from the sepcifications.
+
+.. code:: python
+
+   >>> for i in range(10):
+   ...     for w in bijection.domain.generate_objects_of_size(i):
+   ...         assert w == bijection.inverse_map(bijection.map(w))
+   ...     for w in bijection.codomain.generate_objects_of_size(i):
+   ...         assert w == bijection.map(bijection.inverse_map(w))
+   ...
+
+Whether we find a bijection or not (when one exists) is highly 
+dependent on the packs chosen.
+
+Citing
+######
+
+If you found this library helpful with your research and would like to cite us, 
+you can use the following `BibTeX`_ or go to `Zenodo`_ for alternative formats. 
+
+.. _BibTex: https://zenodo.org/record/4944021/export/hx#.YMcpIC2l30o
+
+.. _Zenodo: https://doi.org/10.5281/zenodo.4944020
\ No newline at end of file
diff --git a/comb_spec_searcher/__init__.py b/comb_spec_searcher/__init__.py
index 9f99e5a73..7dda478e6 100644
--- a/comb_spec_searcher/__init__.py
+++ b/comb_spec_searcher/__init__.py
@@ -16,6 +16,8 @@
     VerificationStrategy,
 )
 
+__version__ = "4.0.0"
+
 __all__ = [
     "CombinatorialSpecificationSearcher",
     "CombinatorialClass",
diff --git a/comb_spec_searcher/bijection.py b/comb_spec_searcher/bijection.py
new file mode 100644
index 000000000..283a5a8c7
--- /dev/null
+++ b/comb_spec_searcher/bijection.py
@@ -0,0 +1,790 @@
+from collections import defaultdict, deque
+from typing import (
+    TYPE_CHECKING,
+    DefaultDict,
+    Dict,
+    Generic,
+    List,
+    Optional,
+    Set,
+    Tuple,
+    TypeVar,
+)
+
+from comb_spec_searcher.comb_spec_searcher import CombinatorialSpecificationSearcher
+from comb_spec_searcher.exception import NoMoreClassesToExpandError
+from comb_spec_searcher.rule_db.base import RuleDB
+from comb_spec_searcher.specification import CombinatorialSpecification
+from comb_spec_searcher.specification_extrator import (
+    EquivalenceRuleExtractor,
+    SpecificationRuleExtractor,
+)
+from comb_spec_searcher.strategies.rule import AbstractRule, Rule
+from comb_spec_searcher.tree_searcher import Node, prune
+from comb_spec_searcher.typing import (
+    CombinatorialClassType,
+    CombinatorialObjectType,
+    RuleKey,
+    Terms,
+)
+
+if TYPE_CHECKING:
+    from comb_spec_searcher.combinatorial_class import (
+        CombinatorialClass,
+        CombinatorialObject,
+    )
+
+
+ClassType1 = TypeVar("ClassType1", bound="CombinatorialClass")
+ObjType1 = TypeVar("ObjType1", bound="CombinatorialObject")
+ClassType2 = TypeVar("ClassType2", bound="CombinatorialClass")
+ObjType2 = TypeVar("ObjType2", bound="CombinatorialObject")
+NO_CONSTRUCTOR = None.__class__
+MatchingInfo = DefaultDict[
+    Tuple[int, int], Dict[Tuple[Tuple[int, ...], Tuple[int, ...]], List[int]]
+]
+MatchingInfoSingle = DefaultDict[int, DefaultDict[int, Set[Tuple[int, ...]]]]
+SpecMap = Dict[int, Tuple[int, ...]]
+RuleClassification = DefaultDict[
+    int,
+    DefaultDict[
+        int,
+        List[
+            Tuple[
+                Tuple[int, ...],
+                AbstractRule[CombinatorialClassType, CombinatorialObjectType],
+            ]
+        ],
+    ],
+]
+EqPathTracker = DefaultDict[
+    Tuple[int, int],
+    DefaultDict[
+        Tuple[int, int],
+        Dict[Tuple[Tuple[int, ...], Tuple[int, ...]], bool],
+    ],
+]
+
+
+class ParallelInfo(Generic[CombinatorialClassType, CombinatorialObjectType]):
+    """Information from searcher essential to finding bijection with another and then
+    to create specifications.
+    """
+
+    def __init__(
+        self,
+        searcher: CombinatorialSpecificationSearcher[CombinatorialClassType],
+        additional_levels: int = 0,
+    ):
+        self.searcher = searcher
+        self._expand(additional_levels)
+        if not isinstance(self.searcher.ruledb, RuleDB):
+            raise RuntimeError("Only searcher supported rule db is `RuleDB`.")
+        self.ruledb: RuleDB = self.searcher.ruledb
+
+        # Root eq label and actual class
+        self.root_eq_label: int = self.searcher.ruledb.equivdb[
+            self.searcher.start_label
+        ]
+        self.root_class = self.searcher.classdb.get_class(self.searcher.start_label)
+
+        # Atom information for comparing them
+        self.atom_map: Dict[int, Tuple[int, Terms]] = {}
+
+        # Used to filter out impossible matchings to avoid unnecessary expansions
+        self.eq_label_rules = self._construct_eq_label_rules()
+
+    def _expand(self, additional_levels: int) -> None:
+        """Expand until at least one spec exists and after that some
+        additional expansions can be done (optionally)."""
+        self._expand_until_spec()
+        self._additional_expands(additional_levels)
+
+    def _expand_until_spec(self) -> None:
+        try:
+            while not self.searcher.has_specification():
+                self.searcher.do_level()
+        except NoMoreClassesToExpandError as ex:
+            if not self.searcher.has_specification():
+                raise ValueError("No specifications were found") from ex
+
+    def _additional_expands(self, additional_levels: int) -> None:
+        try:
+            for _ in range(additional_levels):
+                self.searcher.do_level()
+        except NoMoreClassesToExpandError:
+            pass
+
+    def _construct_eq_label_rules(
+        self,
+    ) -> RuleClassification[CombinatorialClassType, CombinatorialObjectType]:
+        """
+        Creates a dictionary d such that
+            d[eq_label][child_count]
+        is a set of possible resulting children eq_labels and the corresponding
+        constructor if any.
+        """
+        lis = self._pruned_rules_up_to_eq()
+        rule_dict = self.ruledb.rule_from_equivalence_rule_dict(lis)
+        eq_label_rules: RuleClassification[
+            CombinatorialClassType, CombinatorialObjectType
+        ] = defaultdict(lambda: defaultdict(list))
+
+        for eq_par, eq_chi in lis:
+            parent, rule = self._get_class_and_rule(eq_par, eq_chi, rule_dict)
+            assert not parent.is_empty()
+            if parent.is_atom():
+                sz = next(
+                    parent.objects_of_size(parent.minimum_size_of_object())
+                ).size()
+                self.atom_map[eq_par] = (sz, rule.get_terms(sz))
+            else:
+                if not isinstance(rule, Rule):
+                    raise ValueError("Only atoms can be verified.")
+                assert len(eq_chi) > 0
+            eq_label_rules[eq_par][len(eq_chi)].append((eq_chi, rule))
+        return eq_label_rules
+
+    def _get_class_and_rule(
+        self, eq_par: int, eq_chi: Tuple[int, ...], rule_dict: Dict[RuleKey, RuleKey]
+    ) -> Tuple[
+        CombinatorialClassType,
+        AbstractRule[CombinatorialClassType, CombinatorialObjectType],
+    ]:
+        actual_par, actual_children = rule_dict[(eq_par, eq_chi)]
+        strategy = self.ruledb.rule_to_strategy[(actual_par, actual_children)]
+        parent = self.searcher.classdb.get_class(actual_par)
+        rule: AbstractRule[CombinatorialClassType, CombinatorialObjectType] = strategy(
+            parent
+        )
+        return parent, rule
+
+    def _pruned_rules_up_to_eq(
+        self,
+    ) -> List[Tuple[int, Tuple[int, ...]]]:
+        rules_up_to_eq = self.ruledb.rules_up_to_equivalence()
+        prune(rules_up_to_eq)
+        return [(k, c) for k, v in rules_up_to_eq.items() for c in v]
+
+
+class ParallelSpecFinder(Generic[ClassType1, ObjType1, ClassType2, ObjType2]):
+    """
+    Finder for paralell specs. This version assumes that any classes that share
+    equivalence labels are in fact equivalent.
+    """
+
+    _INVALID, _UNKNOWN, _VALID = range(-1, 2)
+
+    def __init__(
+        self,
+        searcher1: CombinatorialSpecificationSearcher[ClassType1],
+        searcher2: CombinatorialSpecificationSearcher[ClassType2],
+    ):
+        self._pre_expand(searcher1, searcher2)
+        self._pi1 = ParallelInfo[ClassType1, ObjType1](searcher1)
+        self._pi2 = ParallelInfo[ClassType2, ObjType2](searcher2)
+        self._ancestors: Set[Tuple[int, int]] = set()
+
+    def _pre_expand(
+        self,
+        searcher1: CombinatorialSpecificationSearcher[ClassType1],
+        searcher2: CombinatorialSpecificationSearcher[ClassType2],
+    ) -> None:
+        """When overwritten, this method can be implemented and will be called
+        prior to the searcher expanding the universe."""
+
+    # First search
+
+    def find(
+        self,
+    ) -> Optional[
+        Tuple[
+            CombinatorialSpecification[ClassType1, ObjType1],
+            CombinatorialSpecification[ClassType2, ObjType2],
+        ]
+    ]:
+        """Find bijections between the two universes."""
+        # d[(p1,p2)][(c1,c2)] = permutation of children such that they match
+        matching_info: MatchingInfo = defaultdict(dict)
+
+        # Failed matches
+        visited: Set[Tuple[int, int]] = set()
+
+        # Find matching eq lables
+        found = self._find(
+            self._pi1.root_eq_label, self._pi2.root_eq_label, matching_info, visited
+        )
+
+        # Try to construct valid specs from said eq-label matchings
+        return self._matching_info_to_specs(found, matching_info)
+
+    def _find(
+        self,
+        id1: int,
+        id2: int,
+        matching_info: MatchingInfo,
+        visited: Set[Tuple[int, int]],
+    ) -> bool:
+        """Recursion call to find matching for eq labels (id1,id2)."""
+
+        # Check for base cases (invalid = -1, valid = 1)
+        known = self._base_case(id1, id2, matching_info, visited)
+        if known:
+            return bool(known + 1)
+
+        # Store pair as ancestors ids checked with further recursions
+        self._ancestors.add((id1, id2))
+
+        for children1, children2, n in self._potential_children(id1, id2):
+            # The permutation of matched children indices
+            child_order: List[int] = [-1] * n
+
+            # Iterative stack with elements:
+            # (index1, index2, indices of the second rule's children already matched)
+            stack = [(0, i, {i}) for i in range(n - 1, -1, -1)]
+
+            # To avoid repeated index pairs, since outside choices have no effect.
+            blacklist: Set[Tuple[int, int]] = set()
+
+            # Check all pairings of children. Any pairing hat exhaust all nonempty
+            # children, such that all are valid will result in current nodes being
+            # matched. This is done with backtracking.
+            while stack:
+                i1, i2, in_use = stack.pop()
+
+                # Already failed
+                if (i1, i2) in blacklist:
+                    continue
+
+                # Recursion. Blacklist children pair if we fail.
+                if not self._find(children1[i1], children2[i2], matching_info, visited):
+                    blacklist.add((i1, i2))
+                    continue
+
+                # Update permutation
+                child_order[i2] = i1
+
+                # If we reach the last index we have succeeded. Gather this matching and
+                # keep trying other children as we want to collect all that are valid.
+                if i1 == n - 1:
+                    matching_info[(id1, id2)][(children1, children2)] = child_order
+                    break
+
+                # Extand stack if we are not done.
+                ParallelSpecFinder._extend_stack(i1, n, in_use, stack)
+
+        # Remove self as ancestors if we are done.
+        self._ancestors.remove((id1, id2))
+
+        # Store pair in memory of failures.
+        visited.add((id1, id2))
+
+        # Return true if any valid matching was found.
+        return (id1, id2) in matching_info
+
+    def _base_case(
+        self,
+        id1: int,
+        id2: int,
+        matching_info: MatchingInfo,
+        visited: Set[Tuple[int, int]],
+    ) -> int:
+        # If no children we compare atoms
+        if (
+            id1 in self._pi1.atom_map
+            and id2 in self._pi2.atom_map
+            and self._atom_match(id1, id2)
+        ):
+            matching_info[(id1, id2)] = {((), ()): []}
+            return ParallelSpecFinder._VALID
+
+        # Already matched
+        if (id1, id2) in matching_info:
+            return ParallelSpecFinder._VALID
+        # Already failed to match
+        if (id1, id2) in visited:
+            return ParallelSpecFinder._INVALID
+        # Recursive match
+        if (id1, id2) in self._ancestors:
+            return ParallelSpecFinder._VALID
+        return ParallelSpecFinder._UNKNOWN
+
+    def _atom_match(self, id1: int, id2: int) -> bool:
+        """Returns true if atoms match."""
+        sz1, terms1 = self._pi1.atom_map[id1]
+        sz2, terms2 = self._pi2.atom_map[id2]
+        return sz1 == sz2 and terms1 == terms2
+
+    def _rule_match(
+        self,
+        rule1: AbstractRule[ClassType1, ObjType1],
+        rule2: AbstractRule[ClassType2, ObjType2],
+    ) -> bool:
+        # pylint: disable=no-self-use
+        if not isinstance(rule1, Rule):
+            return not isinstance(rule2, Rule)
+        if not isinstance(rule2, Rule):
+            return False
+        return rule1.constructor.equiv(rule2.constructor)[0]
+
+    def _potential_children(
+        self, id1: int, id2: int
+    ) -> List[Tuple[Tuple[int, ...], Tuple[int, ...], int]]:
+        """Children are potential children if all of the following holds they have the
+        same number of non-empty children and the constructor's GF matches.
+        """
+        return [
+            (c1, c2, n)
+            for n, child_rule_pairs in self._pi1.eq_label_rules[id1].items()
+            for c1, r1 in child_rule_pairs
+            for c2, r2 in self._pi2.eq_label_rules[id2][n]
+            if self._rule_match(r1, r2)
+        ]
+
+    @staticmethod
+    def _extend_stack(
+        i1: int, n: int, in_use: Set[int], stack: List[Tuple[int, int, Set[int]]]
+    ) -> None:
+        for i in range(n - 1, -1, -1):
+            if i in in_use:
+                continue
+            stack.append((i1 + 1, i, in_use.union({i})))
+
+    # Second search
+
+    def _matching_info_to_specs(
+        self,
+        found: bool,
+        matching_info: MatchingInfo,
+    ) -> Optional[
+        Tuple[
+            CombinatorialSpecification[ClassType1, ObjType1],
+            CombinatorialSpecification[ClassType2, ObjType2],
+        ]
+    ]:
+        # If first search was unsuccessful
+        if not found:
+            return None
+        # Attempt to construct specifications from matches
+        spec_rules = self._search_matching_info(matching_info)
+        if spec_rules is None:
+            return None
+        return (
+            ParallelSpecFinder._create_spec(spec_rules[0], self._pi1),
+            ParallelSpecFinder._create_spec(spec_rules[1], self._pi2),
+        )
+
+    def _search_matching_info(
+        self, matching_info: MatchingInfo
+    ) -> Optional[Tuple[SpecMap, SpecMap]]:
+        # Matching info 1 and 2 are flatten dictionaries of matching_info
+        # for the perspective of each spec. The maps sp1 and sp2 are the
+        # gathered specifications (in terms of eq labels), initially empty.
+        (
+            matching_info1,
+            matching_info2,
+            sp1,
+            sp2,
+        ) = ParallelSpecFinder._search_matching_info_init(matching_info)
+
+        # Recursive helper function that tries to populate the spec maps in a way
+        # that no label appears on the LHS more than one time.
+        def _rec(
+            id1: int,
+            id2: int,
+            id_sets: Tuple[Set[int], Set[int]],
+        ) -> bool:
+            # Check base cases (invalid = -1, valid = 1)
+            bc = ParallelSpecFinder._search_matching_info_recursion_base_cases(
+                id1, id2, matching_info, matching_info1, matching_info2, sp1, sp2
+            )
+            if bc:
+                return bool(bc + 1)
+
+            # If ids have appeared before, we don't want to clean them when done
+            rec1, rec2 = id1 in sp1, id2 in sp2
+
+            # For all children that matches for given parent pair where they are either
+            # not in spec maps or they are mapped to the children we are looking at.
+            for children1, children2 in filter(
+                lambda c: (id1 not in sp1 or c[0] == sp1[id1])
+                and (id2 not in sp2 or c[1] == sp2[id2]),
+                matching_info[(id1, id2)],
+            ):
+                # Add rule for spec to (id -> children) for both specs
+                sp1[id1], sp2[id2] = children1, children2
+
+                # Construct cleaning sets that are passed down the recursion
+                to_clean: Tuple[Set[int], Set[int]] = (set(), set())
+
+                # If all the children are valid, compared one at a time in the order
+                # they were matched in the first search.
+                if all(
+                    _rec(child1, child2, to_clean)
+                    for child1, child2 in zip(
+                        (
+                            children1[i]
+                            for i in matching_info[(id1, id2)][(children1, children2)]
+                        ),
+                        children2,
+                    )
+                ):
+                    # Update the cleaning set so any failing ancestor can remove.
+                    id_sets[0].update(to_clean[0], () if rec1 else (id1,))
+                    id_sets[1].update(to_clean[1], () if rec2 else (id2,))
+                    return True
+                # If failed, remove all descendants that populated the spec maps.
+                ParallelSpecFinder._clean_descendants(
+                    *to_clean, id1, id2, sp1, sp2, rec1, rec2
+                )
+            return False
+
+        if _rec(
+            self._pi1.root_eq_label,
+            self._pi2.root_eq_label,
+            (set(), set()),
+        ):
+            return sp1, sp2
+        return None
+
+    @staticmethod
+    def _search_matching_info_init(
+        matching_info: MatchingInfo,
+    ) -> Tuple[MatchingInfoSingle, MatchingInfoSingle, SpecMap, SpecMap]:
+        matching_info_1: MatchingInfoSingle = defaultdict(lambda: defaultdict(set))
+        matching_info_2: MatchingInfoSingle = defaultdict(lambda: defaultdict(set))
+        # convert (p1,p2) -> children to p1 -> p2 -> children and p2-> p1 -> children.
+        for (p1, p2), children in matching_info.items():
+            for ch1, ch2 in children:
+                matching_info_1[p1][p2].add(ch1)
+                matching_info_2[p2][p1].add(ch2)
+        sp1: SpecMap = {}
+        sp2: SpecMap = {}
+        return matching_info_1, matching_info_2, sp1, sp2
+
+    @staticmethod
+    def _search_matching_info_recursion_base_cases(
+        id1: int,
+        id2: int,
+        matching_info: MatchingInfo,
+        matching_info1: MatchingInfoSingle,
+        matching_info2: MatchingInfoSingle,
+        sp1: SpecMap,
+        sp2: SpecMap,
+    ) -> int:
+        if ParallelSpecFinder._inconsistent_with_matching_info(
+            id1, id2, matching_info, matching_info1, matching_info2, sp1, sp2
+        ):
+            return ParallelSpecFinder._INVALID
+        # If atoms, there is only one possible outcome
+        if ((), ()) in matching_info[(id1, id2)]:
+            sp1[id1], sp2[id2] = (), ()
+            return ParallelSpecFinder._VALID
+        # If both are assigned, we are done
+        if id1 in sp1 and id2 in sp2:
+            return ParallelSpecFinder._VALID
+        return ParallelSpecFinder._UNKNOWN
+
+    @staticmethod
+    def _inconsistent_with_matching_info(
+        id1: int,
+        id2: int,
+        matching_info: MatchingInfo,
+        matching_info1: MatchingInfoSingle,
+        matching_info2: MatchingInfoSingle,
+        sp1: SpecMap,
+        sp2: SpecMap,
+    ):
+        """If ids are not in the matching order or if one is assigned but those
+        children aren't compatible with the other.
+        """
+        return (
+            (id1, id2) not in matching_info
+            or (id1 in sp1 and sp1[id1] not in matching_info1[id1][id2])
+            or (id2 in sp2 and sp2[id2] not in matching_info2[id2][id1])
+        )
+
+    @staticmethod
+    def _clean_descendants(
+        to_clean1: Set[int],
+        to_clean2: Set[int],
+        id1: int,
+        id2: int,
+        sp1: SpecMap,
+        sp2: SpecMap,
+        rec1: bool,
+        rec2: bool,
+    ):
+        # Clean descendants
+        for i in to_clean1:
+            if i in sp1:
+                del sp1[i]
+        for i in to_clean2:
+            if i in sp2:
+                del sp2[i]
+        # Clean self if no prior assignments
+        if not rec1 and id1 in sp1:
+            del sp1[id1]
+        if not rec2 and id2 in sp2:
+            del sp2[id2]
+
+    # specification construction helpers
+
+    @staticmethod
+    def _create_spec(d: SpecMap, pi: ParallelInfo) -> CombinatorialSpecification:
+        rules = SpecificationRuleExtractor(
+            pi.root_eq_label,
+            ParallelSpecFinder._create_tree(d, pi.root_eq_label),
+            pi.ruledb,
+            pi.searcher.classdb,
+        ).rules()
+        return CombinatorialSpecification(pi.root_class, rules)
+
+    @staticmethod
+    def _create_tree(d: SpecMap, root_eq_label: int) -> Node:
+        visited: Set[int] = set()
+        root_node = Node(root_eq_label)
+        queue = deque([root_node])
+        # BFS to consume spec map.
+        while queue:
+            v = queue.popleft()
+            rule = d.get(v.label, ())
+            if not (v.label in visited or rule == ()):
+                children = [Node(i) for i in rule]
+                queue.extend(children)
+                v.children = children
+            visited.add(v.label)
+        return root_node
+
+
+class EqPathParallelSpecFinder(
+    ParallelSpecFinder[ClassType1, ObjType1, ClassType2, ObjType2]
+):
+    """
+    A version of ParallelSpecFinder that supports nonequivalent
+    classes sharing equivalence labels.
+    """
+
+    def __init__(
+        self,
+        searcher1: CombinatorialSpecificationSearcher[ClassType1],
+        searcher2: CombinatorialSpecificationSearcher[ClassType2],
+    ):
+        super().__init__(searcher1, searcher2)
+
+        # Tracks the path taken in the second search
+        self._path: List[Tuple[int, int, int, int]] = [(-1, -1, -1, -1)]
+        self._path_ancestors: Set[Tuple[int, int]] = set()
+
+    def _search_matching_info(
+        self, matching_info: MatchingInfo
+    ) -> Optional[Tuple[SpecMap, SpecMap]]:
+        # This will differ (from super's) that equiv lables must
+        # be validated given the path that has been taken.
+
+        (
+            matching_info1,
+            matching_info2,
+            sp1,
+            sp2,
+        ) = EqPathParallelSpecFinder._search_matching_info_init(matching_info)
+
+        # For storing results of eq path checking.
+        eq_path_tracker: EqPathTracker = defaultdict(lambda: defaultdict(dict))
+
+        def _rec(
+            id1: int,
+            id2: int,
+            id_sets: Tuple[Set[int], Set[int]],
+        ) -> bool:
+            bc = self._search_matching_info_recursion_base_cases_eq(
+                id1,
+                id2,
+                matching_info,
+                matching_info1,
+                matching_info2,
+                sp1,
+                sp2,
+                self._path[-1],
+                eq_path_tracker,
+            )
+            if bc:
+                return bool(bc + 1)
+            rec1, rec2 = id1 in sp1, id2 in sp2
+            for children1, children2 in filter(
+                lambda c: (id1 not in sp1 or c[0] == sp1[id1])
+                and (id2 not in sp2 or c[1] == sp2[id2]),
+                matching_info[(id1, id2)],
+            ):
+                sp1[id1], sp2[id2] = children1, children2
+                to_clean: Tuple[Set[int], Set[int]] = (set(), set())
+
+                # Check if the path required for the eq labels actually matches
+                if self._eq_path_matches(
+                    id1, id2, *self._path[-1], sp1, sp2, eq_path_tracker
+                ):
+                    valid = True
+                    for j2, ((j1, child1), child2) in enumerate(
+                        zip(
+                            (
+                                (i, children1[i])
+                                for i in matching_info[(id1, id2)][
+                                    (children1, children2)
+                                ]
+                            ),
+                            children2,
+                        )
+                    ):
+                        self._path.append((id1, id2, j1, j2))
+                        self._path_ancestors.add((id1, id2))
+                        if not _rec(child1, child2, to_clean):
+                            valid = False
+                            self._path.pop()
+                            self._path_ancestors.remove((id1, id2))
+                            break
+                        self._path.pop()
+                        self._path_ancestors.remove((id1, id2))
+                    if valid:
+                        id_sets[0].update(to_clean[0], () if rec1 else (id1,))
+                        id_sets[1].update(to_clean[1], () if rec2 else (id2,))
+                        return True
+                EqPathParallelSpecFinder._clean_descendants(
+                    *to_clean, id1, id2, sp1, sp2, rec1, rec2
+                )
+            return False
+
+        if _rec(
+            self._pi1.root_eq_label,
+            self._pi2.root_eq_label,
+            (set(), set()),
+        ):
+            return sp1, sp2
+        return None
+
+    def _search_matching_info_recursion_base_cases_eq(
+        self,
+        id1: int,
+        id2: int,
+        matching_info: MatchingInfo,
+        matching_info1: MatchingInfoSingle,
+        matching_info2: MatchingInfoSingle,
+        sp1: SpecMap,
+        sp2: SpecMap,
+        relations: Tuple[int, int, int, int],
+        eq_path_tracker: EqPathTracker,
+    ) -> int:
+        pid1, pid2, idx1, idx2 = relations
+        if EqPathParallelSpecFinder._inconsistent_with_matching_info(
+            id1, id2, matching_info, matching_info1, matching_info2, sp1, sp2
+        ):
+            return EqPathParallelSpecFinder._INVALID
+        # These atoms will have matches in the first searcher but we do offer
+        # additional checks which will, by default, always return true.
+        if ((), ()) in matching_info[(id1, id2)]:
+            if not self._atom_path_match(id1, id2, sp1, sp2):
+                return EqPathParallelSpecFinder._INVALID
+            sp1[id1], sp2[id2] = (), ()
+            return EqPathParallelSpecFinder._VALID
+        # If both ids are set we still need to check if they are valid in terms
+        # of eq paths since we may have arrived from different parents.
+        if id1 in sp1 and id2 in sp2:
+            if self._eq_path_matches(
+                id1, id2, pid1, pid2, idx1, idx2, sp1, sp2, eq_path_tracker
+            ):
+                # If we have already matched id1 and id2 but in a different path, they
+                # can still be incompatible due to this current path's atom labels.
+                should_check = (id1, id2) not in self._path_ancestors
+                if should_check and not self._validate_atoms_for_existing_entries(
+                    id1, id2, sp1, sp2, matching_info, set()
+                ):
+                    return EqPathParallelSpecFinder._INVALID
+                return EqPathParallelSpecFinder._VALID
+            return EqPathParallelSpecFinder._INVALID
+        return EqPathParallelSpecFinder._UNKNOWN
+
+    def _validate_atoms_for_existing_entries(
+        self,
+        id1: int,
+        id2: int,
+        sp1: SpecMap,
+        sp2: SpecMap,
+        matching_info: MatchingInfo,
+        mem: Set[Tuple[int, int]],
+    ) -> bool:
+        """Check if atoms are valid for a given path, when initial ids have already
+        been added to the spec maps elsewhere."""
+        if (id1, id2) in mem:
+            return True
+        children1, children2 = sp1.get(id1, None), sp2.get(id2, None)
+        if children1 is None or children2 is None:
+            return True
+        if children1 == () == children2:
+            return self._atom_path_match(id1, id2, sp1, sp2)
+        mem.add((id1, id2))
+        for j2, ((j1, child1), child2) in enumerate(
+            zip(
+                (
+                    (i, children1[i])
+                    for i in matching_info[(id1, id2)][(children1, children2)]
+                ),
+                children2,
+            )
+        ):
+            self._path.append((id1, id2, j1, j2))
+            check_descendant = self._validate_atoms_for_existing_entries(
+                child1, child2, sp1, sp2, matching_info, mem
+            )
+            self._path.pop()
+            if not check_descendant:
+                return False
+        return True
+
+    def _atom_path_match(self, id1: int, id2: int, sp1: SpecMap, sp2: SpecMap) -> bool:
+        """This can be overwritten if the path can affect the validity of the pair."""
+        # pylint: disable=no-self-use
+        return True
+
+    def _eq_path_matches(
+        self,
+        id1,
+        id2,
+        pid1,
+        pid2,
+        idx1: int,
+        idx2: int,
+        sp1: SpecMap,
+        sp2: SpecMap,
+        cache: EqPathTracker,
+    ) -> bool:
+        """
+        Given (grandparent, parent, children) eq labels in both specs, this will check
+        if the actual rules required (disregarding actual equivalence rules) to traverse
+        this path match.
+        """
+        children = (sp1[id1], sp2[id2])
+        children_cache = cache[(id1, id2)][(pid1, pid2)]
+        if children not in children_cache:
+            path1 = EquivalenceRuleExtractor(
+                self._pi1.root_eq_label,
+                self._pi1.searcher.start_label,
+                EqPathParallelSpecFinder._create_tree(sp1, self._pi1.root_eq_label),
+                self._pi1.ruledb,
+                self._pi1.searcher.classdb,
+                id1,
+                pid1,
+                idx1,
+            ).nonequivalent_rules_in_equiv_path()
+            path2 = EquivalenceRuleExtractor(
+                self._pi2.root_eq_label,
+                self._pi2.searcher.start_label,
+                EqPathParallelSpecFinder._create_tree(sp2, self._pi2.root_eq_label),
+                self._pi2.ruledb,
+                self._pi2.searcher.classdb,
+                id2,
+                pid2,
+                idx2,
+            ).nonequivalent_rules_in_equiv_path()
+            children_cache[children] = len(path1) == len(path2) and all(
+                self._rule_match(r1, r2) for r1, r2 in zip(path1, path2)
+            )
+        return children_cache[children]
diff --git a/comb_spec_searcher/class_db.py b/comb_spec_searcher/class_db.py
index 554819b39..766be2edb 100644
--- a/comb_spec_searcher/class_db.py
+++ b/comb_spec_searcher/class_db.py
@@ -8,7 +8,9 @@
 if is_empty has been checked.
 """
 
+import time
 import zlib
+from datetime import timedelta
 from typing import Dict, Generic, Iterator, Optional, Type, cast
 
 from comb_spec_searcher.typing import ClassKey, CombinatorialClassType, Key
@@ -56,6 +58,8 @@ def __init__(self, combinatorial_class: Type[CombinatorialClassType]):
         self.class_to_info: Dict[ClassKey, Info] = {}
         self.label_to_info: Dict[int, Info] = {}
         self.combinatorial_class = combinatorial_class
+        self._empty_time = 0.0
+        self._empty_num_application = 0
 
     def __iter__(self) -> Iterator[int]:
         """
@@ -177,7 +181,6 @@ def is_empty(
     ) -> bool:
         """
         Return True if combinatorial class is empty set, False if not.
-        Return None if status not set.
         """
         if label is None:
             info = self._get_info(comb_class)
@@ -195,7 +198,11 @@ def is_empty(
     def _is_empty(self, comb_class: CombinatorialClassType) -> bool:
         if not isinstance(comb_class, self.combinatorial_class):
             comb_class = self.get_class(comb_class)
-        return comb_class.is_empty()
+        self._empty_num_application += 1
+        start = time.time()
+        is_empty = comb_class.is_empty()
+        self._empty_time += time.time() - start
+        return is_empty
 
     def set_empty(self, key: Key, empty: bool = True) -> None:
         """
@@ -210,5 +217,7 @@ def status(self) -> str:
         """
         status = "ClassDB status:\n"
         status += "\tTotal number of combinatorial classes found is"
-        status += f" {len(self.label_to_info):,d}"
+        status += f" {len(self.label_to_info):,d}\n"
+        status += f"\tis_empty check applied {self._empty_num_application} time. "
+        status += f"Time spent: {timedelta(seconds=int(self._empty_time))}"
         return status
diff --git a/comb_spec_searcher/comb_spec_searcher.py b/comb_spec_searcher/comb_spec_searcher.py
index 835d31b1d..ab0980086 100644
--- a/comb_spec_searcher/comb_spec_searcher.py
+++ b/comb_spec_searcher/comb_spec_searcher.py
@@ -16,16 +16,14 @@
     Sequence,
     Set,
     Tuple,
-    Union,
     cast,
 )
 
 import logzero
 import tabulate
 from logzero import logger
-from sympy import Eq, Function, var
 
-from comb_spec_searcher.typing import CombinatorialClassType, CSSstrategy, WorkPacket
+from comb_spec_searcher.typing import CombinatorialClassType, CSSstrategy
 
 from .class_db import ClassDB
 from .class_queue import DefaultQueue
@@ -35,18 +33,11 @@
     SpecificationNotFound,
     StrategyDoesNotApply,
 )
-from .rule_db import RuleDB, RuleDBForgetStrategy
-from .rule_db.base import RuleDBBase
+from .rule_db import RuleDB
+from .rule_db.base import RuleDBAbstract
 from .specification import CombinatorialSpecification
-from .specification_extrator import SpecificationRuleExtractor
-from .strategies import (
-    AbstractStrategy,
-    StrategyFactory,
-    StrategyPack,
-    VerificationRule,
-)
+from .strategies import AbstractStrategy, StrategyFactory, StrategyPack
 from .strategies.rule import AbstractRule
-from .tree_searcher import Node
 from .utils import (
     cssiteratortimer,
     cssmethodtimer,
@@ -77,92 +68,85 @@ def __init__(
         self,
         start_class: CombinatorialClassType,
         strategy_pack: StrategyPack,
-        ruledb: Optional[Union[str, RuleDB]] = None,
-        **kwargs,
+        *,
+        ruledb: Optional[RuleDBAbstract] = None,
+        expand_verified: bool = False,
+        debug: bool = False,
     ):
         """
         Initialise CombinatorialSpecificationSearcher.
 
-        OTHER INPUT:
-            - `ruledb`: a string to specify the type of ruledb to use for the
-            search. Default to `None` but can be changed to "forget" for a ruledb that
-            saves more memory.
-            - `expand_verified`: if True, every verified combinatorial class will
-              still be expanded using the strategies in strategy pack
-            - `debug`: if True every rule found will be sanity checked and logged
-              to logging.DEBUG
-            - `function_kwargs` are passed to the call method of strategies
+        INPUTS:
+          - `start_class`: the combinatorial class to search a specification for.
+          - `strategy_pack`: a set of strategy to use for the search.
+          - `ruledb`: a string to specify the type of ruledb to use for the
+          search. Default to `None` but can be changed to "forget" for a ruledb that
+          saves more memory.
+          - `expand_verified`: if True, every verified combinatorial class will
+            still be expanded using the strategies in strategy pack
+          - `debug`: if True every rule found will be sanity checked and logged
+            to logging.DEBUG
         """
         self.strategy_pack = strategy_pack
-        self.debug = kwargs.get("debug", False)
-        self.expand_verified = kwargs.get("expand_verified", False)
+        self.debug = debug
+        self.expand_verified = expand_verified
         if self.debug:
             logzero.loglevel(logging.DEBUG, True)
-        self.kwargs = kwargs.get("function_kwargs", dict())
 
         self.func_times: Dict[str, float] = defaultdict(float)
         self.func_calls: Dict[str, int] = defaultdict(int)
 
-        self.kwargs["symmetry"] = bool(strategy_pack.symmetries)
-
         self.classdb = ClassDB[CombinatorialClassType](type(start_class))
         self.classqueue = DefaultQueue(strategy_pack)
-
-        if ruledb is None:
-            self.ruledb: RuleDBBase = RuleDB()
-        elif ruledb == "forget":
-            self.ruledb = RuleDBForgetStrategy(self.classdb, self.strategy_pack)
-        elif isinstance(ruledb, RuleDBBase):
-            self.ruledb = ruledb
-        else:
-            raise ValueError(
-                "ruledb argument should be None or 'forget' or a RuleDB object"
-            )
+        self.ruledb: RuleDBAbstract = ruledb if ruledb is not None else RuleDB()
+        self.ruledb.link_searcher(self)
 
         # initialise the run with start_class
         self.start_label = self.classdb.get_label(start_class)
-        self._add_to_queue(self.start_label)
+        self.classqueue.add(self.start_label)
         self.tried_to_verify: Set[int] = set()
         self.symmetry_expanded: Set[int] = set()
+        self.inferral_expanded: Set[int] = set()
         self.try_verify(start_class, self.start_label)
         if self.symmetries:
             self._symmetry_expand(start_class, self.start_label)
 
+    def __eq__(self, other: object) -> bool:
+        if not isinstance(other, CombinatorialSpecificationSearcher):
+            return NotImplemented
+        return self.__class__ == other.__class__ and self.__dict__ == other.__dict__
+
     @property
     def verification_strategies(self) -> Sequence[CSSstrategy]:
         """The verification strategies from the strategy pack."""
         return self.strategy_pack.ver_strats
 
-    @property
-    def iterative(self) -> bool:
-        """The iterative parameter from the strategy pack."""
-        return self.strategy_pack.iterative
-
     @property
     def symmetries(self) -> Sequence[CSSstrategy]:
         """The symmetries functions for the strategy pack."""
         return self.strategy_pack.symmetries
 
+    @property
+    def start_class(self) -> CombinatorialClassType:
+        """Returns the start class of the searcher."""
+        return self.classdb.get_class(self.start_label)
+
     def try_verify(self, comb_class: CombinatorialClassType, label: int) -> None:
         """
         Try to verify the combinatorial class.
         """
-        if label not in self.tried_to_verify:
-            for strategy in self.verification_strategies:
-                if self.ruledb.is_verified(label):
-                    break
-                for start_label, end_labels, rule in self._expand_class_with_strategy(
-                    comb_class, strategy, label
-                ):
-                    self._add_rule(start_label, end_labels, rule)
-            self.tried_to_verify.add(label)
-
-    @cssmethodtimer("is empty")
-    def is_empty(self, comb_class: CombinatorialClassType, label: int) -> bool:
-        """Return True if a combinatorial class contains no objects, False
-        otherwise."""
-        empty = self.classdb.is_empty(comb_class, label)
-        return empty
+        if label in self.tried_to_verify:
+            return
+        self.tried_to_verify.add(label)
+        if self.classdb.is_empty(comb_class, label):
+            return
+        for strategy in self.verification_strategies:
+            if self.ruledb.is_verified(label):
+                return
+            for start_label, end_labels, rule in self._expand_class_with_strategy(
+                comb_class, strategy, label
+            ):
+                self.add_rule(start_label, end_labels, rule)
 
     def _expand(
         self,
@@ -182,31 +166,32 @@ def _expand(
                 for start_label, end_labels, rule in self._expand_class_with_strategy(
                     comb_class, strategy_generator, label
                 ):
-                    self._add_rule(start_label, end_labels, rule)
+                    self.add_rule(start_label, end_labels, rule)
 
+    @staticmethod
     def _rules_from_strategy(
-        self, comb_class: CombinatorialClassType, strategy: CSSstrategy
+        comb_class: CombinatorialClassType, strategy: CSSstrategy
     ) -> Iterator[AbstractRule]:
-        """Yield all the rules given by a strategy/strategy generator."""
+        """Yield all the rules given by a strategy/strategy factory."""
         if isinstance(strategy, AbstractStrategy):
             try:
-                yield strategy(comb_class, **self.kwargs)
+                yield strategy(comb_class)
             except StrategyDoesNotApply:
                 pass
         elif isinstance(strategy, StrategyFactory):
-            for strat in strategy(comb_class, **self.kwargs):
+            for strat in strategy(comb_class):
                 if isinstance(strat, AbstractRule):
                     yield strat
                 elif isinstance(strat, AbstractStrategy):
                     try:
-                        yield strat(comb_class, **self.kwargs)
+                        yield strat(comb_class)
                     except StrategyDoesNotApply:
                         continue
                 else:
                     raise InvalidOperationError(
                         "Attempting to add non Rule type. A Strategy "
-                        "Generator's __call__ method should yield Strategy or "
-                        "Strategy(comb_class, children) object."
+                        "Factory's __call__ method should yield strategy or "
+                        "a rule."
                     )
         else:
             raise InvalidOperationError(
@@ -274,54 +259,33 @@ def _expand_class_with_strategy(
                     )
             yield start_label, tuple(end_labels), rule
 
-    def _add_rule(
+    @cssmethodtimer("add rule")
+    def add_rule(
         self, start_label: int, end_labels: Tuple[int, ...], rule: AbstractRule
     ) -> None:
         """
-        Add the cleaned rules labels.
+        Add the rule to the searcher
 
         - try to verify children combinatorial classes
         - set workability of combinatorial classes
-        - remove empty combinatorial classes
         - symmetry expand combinatorial classes
         - add class to classqueue
         """
-        cleaned_end_labels = []
         for comb_class, child_label in zip(rule.children, end_labels):
             if self.symmetries and child_label not in self.symmetry_expanded:
                 self._symmetry_expand(
                     comb_class, child_label
                 )  # TODO: mark symmetries as empty where appropriate
-            # Only applying is_empty check to comb classes that are
-            # possibly empty.
-            if rule.possibly_empty and self.is_empty(comb_class, child_label):
-                logger.debug("Label %s is empty.", child_label)
-                continue
             if rule.workable:
-                self._add_to_queue(child_label)
+                self.classqueue.add(child_label)
             if not rule.inferrable:
-                self._not_inferrable(child_label)
+                self.classqueue.set_not_inferrable(child_label)
             if not rule.possibly_empty:
                 self.classdb.set_empty(child_label, empty=False)
-
             self.try_verify(comb_class, child_label)
-
-            cleaned_end_labels.append(child_label)
-
-        if cleaned_end_labels == [start_label]:
-            return
-
         if rule.ignore_parent:
-            self._stop_yielding(start_label)
-
-        if not cleaned_end_labels:
-            # this must be a verification strategy!
-            assert isinstance(rule, VerificationRule), rule.formal_step
-        self.ruledb.add(start_label, tuple(cleaned_end_labels), rule)
-
-    def _add_empty_rule(self, label: int) -> None:
-        """Mark label as empty. Treated as verified as can count empty set."""
-        self.classdb.set_empty(label, empty=True)
+            self.classqueue.set_stop_yielding(start_label)
+        self.ruledb.add(start_label, end_labels, rule)
 
     def _symmetry_expand(self, comb_class: CombinatorialClassType, label: int) -> None:
         """Add symmetries of combinatorial class to the database."""
@@ -332,7 +296,7 @@ def _symmetry_expand(self, comb_class: CombinatorialClassType, label: int) -> No
             ):
                 sym_label = end_labels[0]
                 self.ruledb.add(start_label, (sym_label,), rule)
-                self._stop_yielding(sym_label)
+                self.classqueue.set_stop_yielding(sym_label)
                 sym_labels.add(sym_label)
         self.symmetry_expanded.update(sym_labels)
 
@@ -352,6 +316,9 @@ def _inferral_expand(
         """
         if self.debug:
             assert comb_class == self.classdb.get_class(label)
+        if label in self.inferral_expanded:
+            return
+        self.inferral_expanded.add(label)
         for i, strategy_generator in enumerate(inferral_strategies):
             if strategy_generator == skip:
                 continue
@@ -360,8 +327,8 @@ def _inferral_expand(
             ):
                 inf_class = rule.children[0]
                 inf_label = end_labels[0]
-                self._add_rule(start_label, end_labels, rule)
-                self._not_inferrable(start_label)
+                self.add_rule(start_label, end_labels, rule)
+                self.classqueue.set_not_inferrable(start_label)
                 inferral_strategies = (
                     inferral_strategies[i + 1 :] + inferral_strategies[0 : i + 1]
                 )
@@ -369,35 +336,11 @@ def _inferral_expand(
                     inf_class, inf_label, inferral_strategies, skip=strategy_generator
                 )
                 break
-        self._not_inferrable(label)
-
-    def get_equations(self, **kwargs) -> Set[Eq]:
-        """
-        Returns a set of equations for all rules currently found.
-        """
-        x = var("x")
-
-        def get_function(comb_class):
-            label = self.classdb.get_label(comb_class)
-            eqv_label = self.ruledb.equivdb[label]
-            return Function("F_{}".format(eqv_label))(x)
-
-        eqs = set()
-        for start, ends, strategy in self.ruledb.all_rules():
-            parent = self.classdb.get_class(start)
-            children = tuple(map(self.classdb.get_class, ends))
-            rule = strategy(parent, children)
-            try:
-                eq = rule.get_equation(get_function)
-            except NotImplementedError:
-                logger.info(
-                    "can't find generating function for %s." " The comb class is:\n%s",
-                    get_function(rule.comb_class),
-                    rule.comb_class,
-                )
-                eq = Eq(get_function(rule.comb_class), Function("NOTIMPLEMENTED")(x))
-            eqs.add(eq)
-        return eqs
+            else:
+                continue
+            # only breaks if inner loop breaks
+            break
+        self.classqueue.set_not_inferrable(label)
 
     def do_level(self) -> None:
         """Expand combinatorial classes in current queue. Combintorial classes
@@ -406,18 +349,6 @@ def do_level(self) -> None:
             comb_class = self.classdb.get_class(label)
             self._expand(comb_class, label, strategies, inferral)
 
-    def _labels_to_expand(self) -> Iterator[WorkPacket]:
-        yield from self.classqueue
-
-    def _add_to_queue(self, label: int):
-        self.classqueue.add(label)
-
-    def _not_inferrable(self, label: int):
-        self.classqueue.set_not_inferrable(label)
-
-    def _stop_yielding(self, label: int):
-        self.classqueue.set_stop_yielding(label)
-
     @cssmethodtimer("status")
     def status(self, elaborate: bool) -> str:
         """
@@ -437,10 +368,9 @@ def status(self, elaborate: bool) -> str:
         status += self.classdb.status() + "\n"
         status += self.classqueue.status() + "\n"
         status += self.ruledb.status(elaborate) + "\n"
-        status += self.mem_status(elaborate)
+        status += self._mem_status(elaborate)
         return status
 
-    @cssmethodtimer("status")
     def _css_status(self, total: float) -> str:
         table: List[Tuple[str, str, timedelta, str]] = []
         for explanation in self.func_calls:
@@ -462,8 +392,7 @@ def _css_status(self, total: float) -> str:
             + "\n"
         )
 
-    @cssmethodtimer("status")
-    def mem_status(self, elaborate: bool) -> str:
+    def _mem_status(self, elaborate: bool) -> str:
         """Provide status information related to memory usage."""
 
         status = "Memory Status:\n"
@@ -501,17 +430,17 @@ def mem_status(self, elaborate: bool) -> str:
         return status
 
     def run_information(self) -> str:
-        """Return string detailing what CombSpecSearcher is looking for."""
+        """Return a string detailing what CombSpecSearcher is looking for."""
         start_string = (
             "Initialising CombSpecSearcher for the combinatorial"
-            " class:\n{}\n".format(self.classdb.get_class(self.start_label))
+            " class:\n{}\n".format(self.start_class)
         )
         start_string += str(self.strategy_pack)
         return start_string
 
     def _log_spec_found(
         self, specification: CombinatorialSpecification, start_time: float
-    ):
+    ) -> None:
         found_string = "Specification built {}\n".format(
             time.strftime("%a, %d %b %Y %H:%M:%S", time.gmtime())
         )
@@ -565,84 +494,69 @@ def auto_search(self, **kwargs) -> CombinatorialSpecification:
         tree that is as small as possible.
         """
         auto_search_start = time.time()
-
-        perc = kwargs.get("perc", 1)
-        if not 0 < perc <= 100:
-            logger.warning(
-                (
-                    "Percentage not between 0 and 100, so assuming 1%"
-                    " search percentage."
-                ),
-            )
-            perc = 1
-        status_update = kwargs.get("status_update", None)
-        max_time = kwargs.get("max_time", None)
-        status_start = time.time()
         start_string = "Auto search started {}\n".format(
             time.strftime("%a, %d %b %Y %H:%M:%S", time.gmtime())
         )
         start_string += self.run_information()
         logger.info(start_string)
+        spec_rules = self._auto_search_rules(**kwargs)
+        if spec_rules is not None:
+            specification = CombinatorialSpecification(self.start_class, spec_rules)
+            self._log_spec_found(specification, auto_search_start)
+            return specification
+        raise SpecificationNotFound
 
-        max_expansion_time = 0
-        expanding = True
-
-        while expanding:
-            expanding, status_start = self._expand_classes_for(
-                max_expansion_time, status_update, status_start, auto_search_start
-            )
-
-            spec_search_start = time.time()
-            logger.debug("Searching for specification.")
-            specification = self.get_specification(
-                smallest=kwargs.get("smallest", False),
-                minimization_time_limit=0.01 * (time.time() - auto_search_start),
-            )
-            if specification is not None:
-                self._log_spec_found(specification, auto_search_start)
-                return specification
-            logger.debug("No specification found.")
-            if max_time is not None:
-                if time.time() - auto_search_start > max_time:
-                    raise ExceededMaxtimeError(
-                        "Exceeded maximum time. Aborting auto search.",
-                    )
-            # worst case, search every hour
-            multiplier = 100 // perc
-            max_expansion_time = min(
-                multiplier * (time.time() - spec_search_start), 3600
-            )
-            logger.debug(
-                "Will expand for %s seconds.",
-                round(max_expansion_time, 2),
-            )
-
-    def _auto_search_rules(self) -> Iterator[AbstractRule]:
+    def _auto_search_rules(
+        self,
+        *,
+        max_expansion_time: Optional[float] = None,
+        perc: int = 1,
+        smallest: bool = False,
+        status_update: Optional[int] = None,
+    ) -> Iterator[AbstractRule]:
         """
-        A basic auto search for returning equivalence paths and rules.
+        The core functionality of the auto_search method.
 
         This method is used by CombinatorialSpecification for expanding
-        verified classes.
-
-        Will raise SpecificationNotFound error if no specification is found
-        after running out of classes to expand.
+        verified classes.  Will raise SpecificationNotFound error if no
+        specification is found after running out of classes to expand.
         """
-        status_start = time.time()
-        status_update = None  # this prevents status updates happening
+        if not 0 < perc <= 100:
+            logger.warning(
+                "Percentage not between 0 and 100, so assuming 1% search percentage."
+            )
+            perc = 1
         auto_search_start = time.time()
-        max_expansion_time: float = 0
+        expansion_time: float = 0
+        status_start = time.time()
         expanding = True
         while expanding:
             expanding, status_start = self._expand_classes_for(
-                max_expansion_time, status_update, status_start, auto_search_start
+                expansion_time, status_update, status_start, auto_search_start
             )
             spec_search_start = time.time()
-            spec = self._get_specification_rules(
-                0.01 * (time.time() - auto_search_start)
+            logger.debug("Searching for specification.")
+            if self.has_specification():
+                logger.info("Specification detected.")
+                return self.ruledb.get_specification_rules(
+                    smallest=smallest,
+                    minimization_time_limit=0.01 * (time.time() - auto_search_start),
+                )
+            logger.debug("No specification found.")
+            if (
+                max_expansion_time is not None
+                and time.time() - auto_search_start > max_expansion_time
+            ):
+                raise ExceededMaxtimeError(
+                    "Exceeded maximum time. Aborting auto search.",
+                )
+            # worst case, search every hour
+            multiplier = 100 / perc
+            expansion_time = min(multiplier * (time.time() - spec_search_start), 3600.0)
+            logger.debug(
+                "Will expand for %s seconds.",
+                round(expansion_time, 2),
             )
-            if spec is not None:
-                return spec
-            max_expansion_time = min(0.01 * (time.time() - spec_search_start), 3600)
         raise SpecificationNotFound
 
     def _expand_classes_for(
@@ -656,13 +570,13 @@ def _expand_classes_for(
         Will expand classes for `expansion_time` seconds.
 
         It will return a pair (bool, time), where the bool is True if there are
-        more classes to expand, False otherwise. The `time` is the initial time
-        for checking whether to post a status update.
+        more classes to expand, False otherwise. The `time` is the time that the
+        last status update printed.
         """
         expansion_start = time.time()
         last_label = None
         expanding = True
-        for label, strategies, inferral in self._labels_to_expand():
+        for label, strategies, inferral in self.classqueue:
             if label != last_label:
                 comb_class = self.classdb.get_class(label)
                 last_label = label
@@ -678,10 +592,14 @@ def _expand_classes_for(
             logger.info("No more classes to expand.")
         return expanding, status_start
 
+    @cssmethodtimer("has specification")
+    def has_specification(self) -> bool:
+        return self.ruledb.has_specification()
+
     @cssmethodtimer("get specification")
     def get_specification(
         self, minimization_time_limit: float = 10, smallest: bool = False
-    ) -> Optional[CombinatorialSpecification]:
+    ) -> CombinatorialSpecification:
         """
         Return a CombinatorialSpecification if the universe contains one.
 
@@ -690,41 +608,12 @@ def get_specification(
         The function will return None if no such CombinatorialSpecification
         exists in the universe.
         """
-        rules = self._get_specification_rules(minimization_time_limit, smallest)
-        if rules is None:
-            return None
-        start_class = self.classdb.get_class(self.start_label)
+        if not self.ruledb.has_specification():
+            raise SpecificationNotFound
+        kwargs = {
+            "minimization_time_limit": minimization_time_limit,
+            "smallest": smallest,
+        }
+        rules = self.ruledb.get_specification_rules(**kwargs)
         logger.info("Creating a specification.")
-        return CombinatorialSpecification(start_class, rules)
-
-    def _get_specification_rules(
-        self, minimization_time_limit: float = 10, smallest: bool = False
-    ) -> Optional[Iterator[AbstractRule]]:
-        node = self._get_specification_node(minimization_time_limit, smallest)
-        if node is None:
-            return None
-        spec_extractor = SpecificationRuleExtractor(
-            self.start_label, node, self.ruledb, self.classdb
-        )
-        return spec_extractor.rules()
-
-    def _get_specification_node(
-        self, minimization_time_limit: float = 10, smallest: bool = False
-    ) -> Optional[Node]:
-        """
-        Return the a specification node if one exists.
-        """
-        try:
-            if smallest:
-                if self.iterative:
-                    raise InvalidOperationError("can't use iterative and smallest")
-                node = self.ruledb.get_smallest_specification(self.start_label)
-            else:
-                node = self.ruledb.find_specification(
-                    self.start_label,
-                    minimization_time_limit=minimization_time_limit,
-                    iterative=self.iterative,
-                )
-        except SpecificationNotFound:
-            return None
-        return node
+        return CombinatorialSpecification(self.start_class, rules)
diff --git a/comb_spec_searcher/combinatorial_class.py b/comb_spec_searcher/combinatorial_class.py
index 8e51ee384..8e82742c3 100644
--- a/comb_spec_searcher/combinatorial_class.py
+++ b/comb_spec_searcher/combinatorial_class.py
@@ -5,9 +5,19 @@
 from functools import reduce
 from importlib import import_module
 from operator import mul
-from typing import Counter, DefaultDict, Dict, Generic, Iterator, List, Tuple, Type
+from typing import (
+    Callable,
+    Counter,
+    DefaultDict,
+    Dict,
+    Generic,
+    Iterator,
+    List,
+    Tuple,
+    Type,
+)
 
-from sympy import Expr, Number, var
+import sympy
 
 from comb_spec_searcher.typing import (
     CombinatorialClassType,
@@ -65,7 +75,18 @@ def to_jsonable(self) -> dict:
             "comb_class": c.__name__,
         }
 
-    def get_genf(self, *args, **kwargs) -> Expr:
+    def get_function(
+        self: CombinatorialClassType, get_label: Callable[[CombinatorialClassType], int]
+    ) -> sympy.Function:
+        """
+        Return a sympy function for the comb class, using the label it is
+        assigned.
+        """
+        x = sympy.var("x")
+        extra_parameters = [sympy.var(k) for k in self.extra_parameters]
+        return sympy.Function("F_{}".format(get_label(self)))(x, *extra_parameters)
+
+    def get_genf(self, *args, **kwargs) -> sympy.Expr:
         """Return the generating function for the combinatorial class"""
         raise NotImplementedError(
             (
@@ -143,20 +164,22 @@ def get_objects(self, n: int) -> Objects:
             objects[param].append(obj)  # pylint: disable=invalid-sequence-index
         return objects
 
-    def initial_conditions(self, check: int = 6) -> List[Expr]:
+    def initial_conditions(self, check: int = 6) -> List[sympy.Expr]:
         """
         Returns a list with the initial conditions to size `check` of the
         CombinatorialClass.
         """
 
-        def monomial(parameters: Dict[str, int]) -> Expr:
+        def monomial(parameters: Dict[str, int]) -> sympy.Expr:
             return reduce(
-                mul, [var(k) ** val for k, val in parameters.items()], Number(1)
+                mul,
+                [sympy.var(k) ** val for k, val in parameters.items()],
+                sympy.Number(1),
             )
 
         return [
             sum(
-                sum(Number(1) for _ in self.objects_of_size(n, **parameters))
+                sum(sympy.Number(1) for _ in self.objects_of_size(n, **parameters))
                 * monomial(parameters)
                 for parameters in self.possible_parameters(n)
             )
@@ -202,10 +225,6 @@ def from_bytes(
         """
         raise NotImplementedError
 
-    @abc.abstractmethod
-    def __init__(self):
-        return
-
     @abc.abstractmethod
     def __eq__(self, other) -> bool:
         return False
diff --git a/comb_spec_searcher/equiv_db.py b/comb_spec_searcher/equiv_db.py
index 1a2c25ae7..16c3d6907 100644
--- a/comb_spec_searcher/equiv_db.py
+++ b/comb_spec_searcher/equiv_db.py
@@ -112,7 +112,6 @@ def connect_cycles(self):
                 for i, vertex in enumerate(path[:-1]):
                     if self.equivalent(vertex, new_end):
                         # we've detected a cycle!
-                        print("CYCLE", path, new_end)
                         for eqv_vertix in path[i:]:
                             self._set_equivalent(eqv_vertix, new_end)
                         break
diff --git a/comb_spec_searcher/isomorphism.py b/comb_spec_searcher/isomorphism.py
index 7a39f786f..7a2c571f3 100644
--- a/comb_spec_searcher/isomorphism.py
+++ b/comb_spec_searcher/isomorphism.py
@@ -1,20 +1,34 @@
-from typing import TYPE_CHECKING, Callable, Dict, List, Optional, Set, Tuple
-
-from .combinatorial_class import CombinatorialClass, CombinatorialObject
+from itertools import product
+from typing import (
+    TYPE_CHECKING,
+    Any,
+    Callable,
+    Dict,
+    Generic,
+    List,
+    Optional,
+    Set,
+    Tuple,
+    TypeVar,
+)
+
+from .combinatorial_class import CombinatorialClass
 from .strategies.rule import AbstractRule, Rule
 
 if TYPE_CHECKING:
+    from .combinatorial_class import CombinatorialObject
     from .specification import CombinatorialSpecification
 
 
-AtomEquals = Callable[["CombinatorialClass", "CombinatorialClass"], bool]
-
-
-def _default_equals(c1: CombinatorialClass, c2: CombinatorialClass) -> bool:
-    return c1 == c2
+ObjType1 = TypeVar("ObjType1", bound="CombinatorialObject")
+ClassType1 = TypeVar("ClassType1", bound="CombinatorialClass")
+ClassType2 = TypeVar("ClassType2", bound="CombinatorialClass")
+ObjType2 = TypeVar("ObjType2", bound="CombinatorialObject")
+OrderMap = Dict[Tuple[ClassType1, ClassType2], List[int]]
+IndexDataMap = Optional[Dict[Tuple[ClassType1, ClassType2], object]]
 
 
-class Isomorphism:
+class Isomorphism(Generic[ClassType1, ObjType1, ClassType2, ObjType2]):
     """Isomorphism checker."""
 
     _INVALID, _UNKNOWN, _VALID = range(-1, 2)
@@ -22,133 +36,242 @@ class Isomorphism:
     @classmethod
     def check(
         cls,
-        spec1: "CombinatorialSpecification",
-        spec2: "CombinatorialSpecification",
-        eq: AtomEquals = _default_equals,
+        spec1: "CombinatorialSpecification[ClassType1, ObjType1]",
+        spec2: "CombinatorialSpecification[ClassType2, ObjType2]",
     ) -> bool:
         """Check if two specs are isomorphic."""
-        return cls(spec1, spec2, eq).are_isomorphic()
+        return cls(spec1, spec2).are_isomorphic()
 
     def __init__(
         self,
-        spec1: "CombinatorialSpecification",
-        spec2: "CombinatorialSpecification",
-        eq: AtomEquals = _default_equals,
+        spec1: "CombinatorialSpecification[ClassType1, ObjType1]",
+        spec2: "CombinatorialSpecification[ClassType2, ObjType2]",
     ) -> None:
-        self._eq = eq
-        self._index = 0
-        self._root1: CombinatorialClass = spec1.root
-        self._rules1: Dict[CombinatorialClass, AbstractRule] = spec1.rules_dict
-        self._ancestors1: Dict[CombinatorialClass, int] = {}
-        self._root2: CombinatorialClass = spec2.root
-        self._rules2: Dict[CombinatorialClass, AbstractRule] = spec2.rules_dict
-        self._ancestors2: Dict[CombinatorialClass, int] = {}
-        self._order_map: Dict[
-            Tuple[CombinatorialClass, CombinatorialClass], List[int]
-        ] = {}
-        self._isomorphic = self._are_isomorphic(self._root1, self._root2)
+        # The specs
+        self._rules1: Dict[
+            ClassType1, AbstractRule[ClassType1, ObjType1]
+        ] = spec1.rules_dict
+        self._rules2: Dict[
+            ClassType2, AbstractRule[ClassType2, ObjType2]
+        ] = spec2.rules_dict
+        self.root1 = spec1.root
+        self.root2 = spec2.root
+
+        # Tracking recursion
+        self._ancestors: Set[Tuple[ClassType1, ClassType2]] = set()
+
+        # Tracking matchin permutations
+        self._order_map: Dict[Tuple[ClassType1, ClassType2], List[int]] = {}
+
+        # Tracking failed matchings
+        self._failed: Set[Tuple[ClassType1, ClassType2]] = set()
+
+        # Tracking additional data
+        self._index_data: Dict[Tuple[ClassType1, ClassType2], object] = {}
+
+        # Tracking paths in terms of child indices
+        self._path_tracker: List[Tuple[Tuple[int, ...], Tuple[int, ...]]] = []
+
+        # Store results
+        self._isomorphic = self._are_isomorphic(spec1.root, spec2.root)
+
+    def get_order_data(
+        self,
+    ) -> Dict[Tuple[ClassType1, ClassType2], object]:
+        """Return any index data for nonbijective matches."""
+        return self._index_data
 
     def are_isomorphic(self) -> bool:
         """Check if the two specs are isomorphic."""
         return self._isomorphic
 
-    def get_order(self, c1: CombinatorialClass, c2: CombinatorialClass) -> List[int]:
+    def get_order(
+        self,
+    ) -> OrderMap[ClassType1, ClassType2]:
         """Get order map of corresponding nodes."""
-        return self._order_map[(c1, c2)]
+        return self._order_map
 
-    def _are_isomorphic(
-        self, node1: CombinatorialClass, node2: CombinatorialClass
-    ) -> bool:
-        # If there are equivilances, we use the 'latest' one
-        node1, node2 = Isomorphism._get_eq_descendant(
-            node1, self._rules1, node2, self._rules2
+    def _are_isomorphic(self, node1: ClassType1, node2: ClassType2) -> bool:
+        # Both ends of equivilance paths (len = 1 if none)
+        eq_path1, eq_path2 = self._get_eq_descendant(node1, node2)
+
+        # The ones that are not on the LHS of equivalence rules
+        curr1, curr2 = eq_path1[-1], eq_path2[-1]
+
+        # Rules from curr1 and curr2 to their respective children
+        rule1, rule2 = self._rules1[curr1], self._rules2[curr2]
+
+        # Any empty children are ignored
+        non_empty_ind1, non_empty_ind2 = (
+            tuple(i for i, c in enumerate(rule1.children) if not c.is_empty()),
+            tuple(i for i, c in enumerate(rule2.children) if not c.is_empty()),
         )
-        rule1, rule2 = self._rules1[node1], self._rules2[node2]
 
-        is_base_case = self._base_cases(node1, rule1, node2, rule2)
+        # Check for base cases (invalid = -1, valid = 1)
+        is_base_case = self._base_cases(
+            eq_path1, rule1, non_empty_ind1, eq_path2, rule2, non_empty_ind2
+        )
         if is_base_case:
             return bool(is_base_case + 1)
 
-        self._index += 1
-        self._ancestors1[node1], self._ancestors2[node2] = (self._index, self._index)
+        # Update ancestors for recursion
+        self._ancestors.update(product(eq_path1, eq_path2))
+
+        # The number of nonempty children
+        n = len(non_empty_ind1)
+
+        # The permutation of matched children indices
+        child_order: List[int] = [-1] * n
+
+        # To avoid repeated index pairs, since outside choices have no effect.
+        blacklist: Set[Tuple[int, int]] = set()
 
-        # Check all matches of children, if any are valid then trees are isomorphic
-        n = len(rule1.children)
-        child_order: List[int] = [0] * n
+        # Iterative stack with elements:
+        # (index1, index2, indices of the second rule's children already matched)
         stack = [(0, i, {i}) for i in range(n - 1, -1, -1)]
+
+        # Check all pairings of children. Any pairing hat exhaust all nonempty children,
+        # such that all are valid will result in current nodes being matched. This is
+        # done with backtracking.
         while stack:
             i1, i2, in_use = stack.pop()
-            if not self._are_isomorphic(rule1.children[i1], rule2.children[i2]):
+
+            # Already failed
+            if (i1, i2) in blacklist:
+                continue
+
+            # Track path
+            self._path_tracker_push(
+                non_empty_ind1[i1],
+                non_empty_ind2[i2],
+                len(eq_path1) > 1,
+                len(eq_path2) > 1,
+            )
+
+            # Call _are_isomorphic recursively for current pair
+            # and add to blacklist and try again if we fail.
+            if not self._are_isomorphic(
+                rule1.children[non_empty_ind1[i1]], rule2.children[non_empty_ind2[i2]]
+            ):
+                blacklist.add((i1, i2))
+                self._path_tracker_pop()
                 continue
-            child_order[i1] = i2
+            self._path_tracker_pop()
+
+            # Update permutation
+            child_order[i2] = i1
+
+            # If we reach the last index we have succeeded
             if i1 == n - 1:
-                self._order_map[(node1, node2)] = child_order
-                self._cleanup(node1, node2)
+                self._order_map[(curr1, curr2)] = child_order
+                # Since we did not conclude this match by recursion we can remove our
+                # current ids from the ancestor set since the next ones we check are
+                # not descendants of the current ids (or their equivalences).
+                self._ancestors.difference_update(product(eq_path1, eq_path2))
                 return True
+
+            # If we have indices to expand, that is done here
             Isomorphism._extend_stack(i1, n, in_use, stack)
 
-        self._cleanup(node1, node2)
+        # If stack has been exhausted we have failed to match the ids. We remove them
+        # from the ancestors set, memorize the failure and remove any data gathered.
+        self._ancestors.difference_update(product(eq_path1, eq_path2))
+        self._failed.add((curr1, curr2))
+        self._index_data.pop((curr1, curr2), None)
         return False
 
-    @staticmethod
     def _get_eq_descendant(
-        node1: CombinatorialClass,
-        rules1: Dict[CombinatorialClass, AbstractRule],
-        node2: CombinatorialClass,
-        rules2: Dict[CombinatorialClass, AbstractRule],
-    ) -> Tuple[CombinatorialClass, CombinatorialClass]:
-        rule1, rule2 = rules1[node1], rules2[node2]
-
-        if rule2.is_equivalence():
-            node2 = rule2.children[0]
-
+        self, node1: ClassType1, node2: ClassType2
+    ) -> Tuple[List[ClassType1], List[ClassType2]]:
+        """Get path to nodes that are not on the LHS of equivalence rules."""
+        rule1, rule2 = self._rules1[node1], self._rules2[node2]
+        nodes1, nodes2 = [node1], [node2]
         if rule1.is_equivalence():
-            node1 = rule1.children[0]
-
-        return node1, node2
+            nodes1.append(rule1.children[0])
+        if rule2.is_equivalence():
+            nodes2.append(rule2.children[0])
+        return nodes1, nodes2
 
     def _base_cases(
         self,
-        node1: CombinatorialClass,
-        rule1: AbstractRule,
-        node2: CombinatorialClass,
-        rule2: AbstractRule,
+        eq_nodes1: List[ClassType1],
+        rule1: AbstractRule[ClassType1, ObjType1],
+        non_empty_children1: Tuple[int, ...],
+        eq_nodes2: List[ClassType2],
+        rule2: AbstractRule[ClassType2, ObjType2],
+        non_empty_children2: Tuple[int, ...],
     ) -> int:
-        # If different number of children
-        if len(rule1.children) != len(rule2.children):
+        # The ones that are not on the LHS of equivalence rules
+        curr1, curr2 = eq_nodes1[-1], eq_nodes2[-1]
+
+        # Already matched
+        if (curr1, curr2) in self._order_map:
+            return Isomorphism._VALID
+
+        # Already failed
+        if (curr1, curr2) in self._failed:
+            return Isomorphism._INVALID
+
+        # If different number of non-empty children
+        if len(non_empty_children1) != len(non_empty_children2):
             return Isomorphism._INVALID
 
-        # If leaf that has no equal-ancestors
+        # Atoms
         if not rule1.children and not rule2.children:
-            # If one is atom, both should be and should also be equal
-            if node1.is_atom() and node2.is_atom() and self._eq(node1, node2):
-                return Isomorphism._VALID
+            if curr1.is_atom() and curr2.is_atom():
+                if self._atom_match(curr1, curr2, rule1, rule2):
+                    return Isomorphism._VALID
             return Isomorphism._INVALID
 
         # If different type of rules are applied, the trees are not isomorphic
-        # TODO: Use == when implemented
         assert isinstance(rule1, Rule) and isinstance(rule2, Rule)
-        if rule1.constructor.__class__ != rule2.constructor.__class__:
-            return Isomorphism._INVALID
-
-        # Check if we have seen this node before
-        ind1, ind2 = self._ancestors1.get(node1, -1), self._ancestors2.get(node2, -1)
-
-        # If one is found and the other not or if
-        # they occured in different places before
-        if ind1 != ind2:
+        if not self._constructor_match(rule1, rule2, curr1, curr2):
             return Isomorphism._INVALID
 
-        # If they have occured before and that occurrence matches in both trees
-        if ind1 == ind2 != -1:
+        # Check for recursive match
+        if any((n1, n2) in self._ancestors for n1, n2 in product(eq_nodes1, eq_nodes2)):
             return Isomorphism._VALID
 
         return Isomorphism._UNKNOWN
 
-    def _cleanup(self, node1: CombinatorialClass, node2: CombinatorialClass) -> None:
-        self._index -= 1
-        del self._ancestors1[node1]
-        del self._ancestors2[node2]
+    def _path_tracker_push(
+        self, c1: int, c2: int, eq_path1: bool, eq_path2: bool
+    ) -> None:
+        """Store the path being taken."""
+        t1 = (0, c1) if eq_path1 else (c1,)
+        t2 = (0, c2) if eq_path2 else (c2,)
+        self._path_tracker.append((t1, t2))
+
+    def _path_tracker_pop(self) -> None:
+        """Remove last stored step in path."""
+        self._path_tracker.pop()
+
+    def _atom_match(
+        self,
+        atom1: ClassType1,
+        atom2: ClassType2,
+        rule1: AbstractRule[ClassType1, ObjType1],
+        rule2: AbstractRule[ClassType2, ObjType2],
+    ) -> bool:
+        """Returns true if atoms match, false otherwise."""
+        # pylint: disable=no-self-use
+        sz1 = next(atom1.objects_of_size(atom1.minimum_size_of_object())).size()
+        sz2 = next(atom2.objects_of_size(atom2.minimum_size_of_object())).size()
+        return sz1 == sz2 and rule1.get_terms(sz1) == rule2.get_terms(sz2)
+
+    def _constructor_match(
+        self,
+        rule1: Rule[ClassType1, ObjType1],
+        rule2: Rule[ClassType2, ObjType2],
+        curr1: ClassType1,
+        curr2: ClassType2,
+    ) -> bool:
+        are_eq, data = rule1.constructor.equiv(rule2.constructor)
+        if not are_eq:
+            return False
+        if data is not None:
+            self._index_data[(curr1, curr2)] = data
+        return True
 
     @staticmethod
     def _extend_stack(
@@ -160,114 +283,246 @@ def _extend_stack(
             stack.append((i1 + 1, i, in_use.union({i})))
 
 
-class Node:
+class Node(Generic[ClassType1, ObjType1, ClassType2, ObjType2]):
     """Parse tree node."""
 
     def __init__(
         self,
-        rule: AbstractRule,
-        obj: CombinatorialObject,
-        get_rule: Callable[[CombinatorialClass], AbstractRule],
+        rule: AbstractRule[ClassType1, ObjType1],
+        obj: ObjType1,
+        get_rule: Callable[[ClassType1], AbstractRule[ClassType1, ObjType1]],
     ):
-        self.rule = rule
         self.obj = obj
+        self._rule = rule
         if not rule.children:
-            self.children: Tuple[Optional["Node"], ...] = tuple()
+            self._children: Tuple[
+                Optional["Node[ClassType1, ObjType1, ClassType2, ObjType2]"], ...
+            ] = tuple()
         else:
             assert isinstance(rule, Rule)
-            self.children = (
-                tuple(
-                    type(self)(get_rule(child), child_obj, get_rule)
-                    if child_obj is not None
-                    else None
-                    for child, child_obj in zip(rule.children, rule.forward_map(obj))
-                )
-                if rule.children
-                else tuple()
+            children, self.idx = rule.indexed_forward_map(obj)
+            self._children = tuple(
+                type(self)(get_rule(child), child_obj, get_rule)
+                if child_obj is not None
+                else None
+                for child, child_obj in zip(rule.children, children)
+                if not child.is_empty()
             )
 
     def build_obj(
         self,
-        rule: AbstractRule,
-        get_order: Callable[[CombinatorialClass, CombinatorialClass], List[int]],
-        get_rule: Callable[[CombinatorialClass], AbstractRule],
-    ) -> CombinatorialObject:
+        rule: AbstractRule[ClassType2, ObjType2],
+        get_order: OrderMap[ClassType1, ClassType2],
+        get_rule: Callable[[ClassType2], AbstractRule[ClassType2, ObjType2]],
+        index_data: Dict[Tuple[ClassType1, ClassType2], object],
+    ) -> ObjType2:
         """Parse tree's recursive build object function."""
-        if not self.children:
-            # TODO: stop special casing verification rules!
-            return self.obj
+
+        atom = not self._children
+        if atom:
+            obj: ObjType2 = next(
+                rule.comb_class.objects_of_size(
+                    rule.comb_class.minimum_size_of_object()
+                )
+            )
+            return obj
+
         if rule.is_equivalence():
-            if not self.rule.is_equivalence():
-                return self.build_obj(get_rule(rule.children[0]), get_order, get_rule)
-            order = [0]
-        elif self.rule.is_equivalence():
-            assert self.children[0] is not None
-            return self.children[0].build_obj(rule, get_order, get_rule)
+            if not self._rule.is_equivalence():
+                assert isinstance(rule, Rule)
+                val: ObjType2 = rule.indexed_backward_map(
+                    (
+                        self.build_obj(
+                            get_rule(rule.children[0]), get_order, get_rule, index_data
+                        ),
+                    ),
+                    self.idx,
+                )
+                return val
+            order = [0]  # If both are equivalences
+        elif self._rule.is_equivalence():
+            assert self._children[0] is not None
+            return self._children[0].build_obj(rule, get_order, get_rule, index_data)
         else:
-            order = get_order(self.rule.comb_class, rule.comb_class)
-        children = tuple(self.children[idx] for idx in order)
-        child_objs = tuple(
-            None
-            if child is None
-            else child.build_obj(get_rule(child_class), get_order, get_rule)
-            for child, child_class in zip(children, rule.children)
-        )
+            order = get_order[(self._rule.comb_class, rule.comb_class)]
+        children = (self._children[idx] for idx in order)
         assert isinstance(rule, Rule)
-        val: CombinatorialObject = next(rule.backward_map(child_objs))
+
+        # Deal with avoiding empty children and supplying ibm with data if any
+        val = rule.indexed_backward_map(
+            tuple(
+                c[0].build_obj(get_rule(c[1]), get_order, get_rule, index_data)
+                if c is not None and c[0] is not None
+                else None
+                for c in (
+                    None if child.is_empty() else (next(children), child)
+                    for child in rule.children
+                )
+            ),
+            self.idx,
+            index_data.get((self._rule.comb_class, rule.comb_class)),
+        )
         return val
 
 
-class ParseTree:
+class ParseTree(Generic[ClassType1, ObjType1, ClassType2, ObjType2]):
     """Parse tree that parses an object and can reconstruct it for the other
     specification."""
 
-    def __init__(self, obj: CombinatorialObject, spec: "CombinatorialSpecification"):
-        self.root = Node(spec.root_rule, obj, spec.get_rule)
+    def __init__(
+        self, obj: ObjType1, spec: "CombinatorialSpecification[ClassType1, ObjType1]"
+    ):
+        self._root = Node[ClassType1, ObjType1, ClassType2, ObjType2](
+            spec.root_rule, obj, spec.get_rule
+        )
 
     def build_obj(
         self,
-        root: AbstractRule,
-        get_order: Callable[[CombinatorialClass, CombinatorialClass], List[int]],
-        get_rule: Callable[[CombinatorialClass], AbstractRule],
-    ) -> CombinatorialObject:
+        root: AbstractRule[ClassType2, ObjType2],
+        get_order: OrderMap[ClassType1, ClassType2],
+        get_rule: Callable[[ClassType2], AbstractRule[ClassType2, ObjType2]],
+        index_data: Dict[Tuple[ClassType1, ClassType2], object],
+    ) -> ObjType2:
         """Build object from the other specification."""
-        return self.root.build_obj(root, get_order, get_rule)
+        return self._root.build_obj(root, get_order, get_rule, index_data)
 
 
-class Bijection:
+class Bijection(Generic[ClassType1, ObjType1, ClassType2, ObjType2]):
     """A bijection object that contains a map between two specifications."""
 
     @classmethod
     def construct(
         cls,
-        spec: "CombinatorialSpecification",
-        other: "CombinatorialSpecification",
-        eq: AtomEquals = _default_equals,
-    ) -> Optional["Bijection"]:
-        """Create a bijection object between two specifications if possible.
-
-        raises: ValueError if specifications are not isomorphic.
-        """
-        try:
-            return cls(spec, other, eq)
-        except ValueError:
+        spec: "CombinatorialSpecification[ClassType1, ObjType1]",
+        other: "CombinatorialSpecification[ClassType2, ObjType2]",
+    ) -> Optional["Bijection[ClassType1, ObjType1, ClassType2, ObjType2]"]:
+        """Create a bijection object between two specifications if possible."""
+        iso = Isomorphism[ClassType1, ObjType1, ClassType2, ObjType2](spec, other)
+        if not iso.are_isomorphic():
             return None
+        return cls(spec, other, iso.get_order(), iso.get_order_data())
 
     def __init__(
         self,
-        spec: "CombinatorialSpecification",
-        other: "CombinatorialSpecification",
-        eq: AtomEquals = _default_equals,
+        spec: "CombinatorialSpecification[ClassType1, ObjType1]",
+        other: "CombinatorialSpecification[ClassType2, ObjType2]",
+        get_order: OrderMap[ClassType1, ClassType2],
+        index_data: IndexDataMap = None,
     ):
-        self.spec = spec
-        self.other = other
-        self.iso = Isomorphism(spec, other, eq)
-        if not self.iso.are_isomorphic():
-            raise ValueError("Specifications are not isomorphic")
-
-    def map(self, obj: CombinatorialObject) -> CombinatorialObject:
-        """Map an object of the root of one specification to the root of the other."""
-        parse_tree = ParseTree(obj, self.spec)
+        self._index_data = {} if index_data is None else index_data
+        self._inv_index_data = {
+            (t2, t1): data for (t1, t2), data in self._index_data.items()
+        }
+        self._spec = spec
+        self._other = other
+        self._get_order = get_order
+        self._get_inverse_order = {
+            (c2, c1): Bijection._perm_inv(lis) for ((c1, c2), lis) in get_order.items()
+        }
+
+    @property
+    def domain(self) -> "CombinatorialSpecification[ClassType1, ObjType1]":
+        """The root class of the first specification."""
+        return self._spec
+
+    @property
+    def codomain(self) -> "CombinatorialSpecification[ClassType2, ObjType2]":
+        """The root class of the second specification."""
+        return self._other
+
+    @staticmethod
+    def _perm_inv(perm: List[int]) -> List[int]:
+        inv = [0] * len(perm)
+        for i, v in enumerate(perm):
+            inv[v] = i
+        return inv
+
+    def map(self, obj: ObjType1) -> ObjType2:
+        """Map an object of the domain to an object of the codomain."""
+        parse_tree = ParseTree[ClassType1, ObjType1, ClassType2, ObjType2](
+            obj, self._spec
+        )
+        return parse_tree.build_obj(
+            self._other.root_rule,
+            self._get_order,
+            self._other.get_rule,
+            self._index_data,
+        )
+
+    def inverse_map(self, obj: ObjType2) -> ObjType1:
+        """Map an object of the codomain to an object of the domain."""
+        # Swap types and use inverse versions of tuple maps
+        parse_tree = ParseTree[ClassType2, ObjType2, ClassType1, ObjType1](
+            obj, self._other
+        )
         return parse_tree.build_obj(
-            self.other.root_rule, self.iso.get_order, self.other.get_rule
+            self._spec.root_rule,
+            self._get_inverse_order,
+            self._spec.get_rule,
+            self._inv_index_data,
         )
+
+    def to_jsonable(self) -> dict:
+        """Return a JSON serializable dictionary for the bijection."""
+        id_map, classes = self._classes_to_array()
+        reconstructed_order: Dict[str, Dict[str, List[int]]] = {}
+        reconstructed_data: Dict[str, Dict[str, object]] = {}
+        Bijection._populate_json_map(self._get_order, reconstructed_order, id_map)
+        Bijection._populate_json_map(self._index_data, reconstructed_data, id_map)
+        return {
+            "spec": self._spec.to_jsonable(),
+            "other": self._other.to_jsonable(),
+            "order": reconstructed_order,
+            "index_data": reconstructed_data,
+            "classes": classes,
+        }
+
+    def _classes_to_array(self) -> Tuple[Dict[CombinatorialClass, str], List[dict]]:
+        """Dump classes into an array with indices as keys to simplify maps."""
+        id_map: Dict[CombinatorialClass, str] = {}
+        classes: List[dict] = []
+        for c1, c2 in self._get_order.keys():
+            if c1 not in id_map:
+                id_map[c1] = f"{len(classes)}"
+                classes.append(c1.to_jsonable())
+            if c2 not in id_map:
+                id_map[c2] = f"{len(classes)}"
+                classes.append(c2.to_jsonable())
+        return id_map, classes
+
+    @staticmethod
+    def _populate_json_map(
+        tuple_map: Dict[Tuple[ClassType1, ClassType2], Any],
+        json_map: Dict[str, Dict[str, Any]],
+        id_map: Dict[CombinatorialClass, str],
+    ) -> None:
+        """Convert tuple maps into valid json maps."""
+        for (c1, c2), value in tuple_map.items():
+            i1, i2 = id_map[c1], id_map[c2]
+            if i1 not in json_map:
+                json_map[i1] = {i2: value}
+            else:
+                json_map[i1][i2] = value
+
+    @classmethod
+    def from_dict(cls, d) -> "Bijection":
+        """Return the bijection with the dictionary outputter by
+        the 'to_jsonable' method.
+        """
+        # pylint: disable=import-outside-toplevel
+        from comb_spec_searcher import CombinatorialSpecification
+
+        spec1 = CombinatorialSpecification[ClassType1, ObjType1].from_dict(d["spec"])
+        spec2 = CombinatorialSpecification[ClassType2, ObjType2].from_dict(d["other"])
+        order = [CombinatorialClass.from_dict(c) for c in d["classes"]]
+        get_order = {
+            (order[int(idx1)], order[int(idx2)]): lis
+            for idx1, sub_map in d["order"].items()
+            for idx2, lis in sub_map.items()
+        }
+        index_data = {
+            (order[int(idx1)], order[int(idx2)]): data
+            for idx1, sub_map in d["index_data"].items()
+            for idx2, data in sub_map.items()
+        }
+        return cls(spec1, spec2, get_order, index_data)  # type: ignore
diff --git a/comb_spec_searcher/rule_db/__init__.py b/comb_spec_searcher/rule_db/__init__.py
index b4db4401a..394eb1524 100644
--- a/comb_spec_searcher/rule_db/__init__.py
+++ b/comb_spec_searcher/rule_db/__init__.py
@@ -1,5 +1,5 @@
 from .base import RuleDB
+from .forest import RuleDBForest
 from .forget import RuleDBForgetStrategy
-from .limited_strategy import LimitedStrategyRuleDB
 
-__all__ = ["RuleDB", "RuleDBForgetStrategy", "LimitedStrategyRuleDB"]
+__all__ = ["RuleDB", "RuleDBForgetStrategy", "RuleDBForest"]
diff --git a/comb_spec_searcher/rule_db/abstract.py b/comb_spec_searcher/rule_db/abstract.py
new file mode 100644
index 000000000..9388ff427
--- /dev/null
+++ b/comb_spec_searcher/rule_db/abstract.py
@@ -0,0 +1,103 @@
+"""
+Definition of the base feature of any ruledb.
+"""
+import abc
+from typing import (
+    TYPE_CHECKING,
+    Any,
+    Callable,
+    Iterator,
+    Optional,
+    Tuple,
+    TypeVar,
+    cast,
+)
+
+from comb_spec_searcher.class_db import ClassDB
+from comb_spec_searcher.exception import SpecificationNotFound
+from comb_spec_searcher.strategies.rule import AbstractRule
+from comb_spec_searcher.strategies.strategy_pack import StrategyPack
+
+if TYPE_CHECKING:
+    from comb_spec_searcher import CombinatorialSpecificationSearcher
+
+RuleDBMethod = TypeVar("RuleDBMethod", bound=Callable[..., Any])
+
+
+__all__ = ["ensure_specification", "RuleDBAbstract"]
+
+
+class RuleDBAbstract(abc.ABC):
+    NOT_LINK_ERROR = RuntimeError(
+        "RuleDB is not linked with the searcher. Call `link_searcher` before using it."
+    )
+
+    def __init__(self) -> None:
+        self._searcher: Optional["CombinatorialSpecificationSearcher"] = None
+
+    def link_searcher(self, searcher: "CombinatorialSpecificationSearcher") -> None:
+        if self._searcher is not None:
+            raise RuntimeError("Searcher is alreay linked")
+        self._searcher = searcher
+
+    @property
+    def searcher(self) -> "CombinatorialSpecificationSearcher":
+        if self._searcher is None:
+            raise RuleDBAbstract.NOT_LINK_ERROR
+        return self._searcher
+
+    @property
+    def classdb(self) -> ClassDB:
+        return self.searcher.classdb
+
+    @property
+    def strategy_pack(self) -> StrategyPack:
+        return self.searcher.strategy_pack
+
+    @property
+    def root_label(self) -> int:
+        return self.searcher.start_label
+
+    @abc.abstractmethod
+    def is_verified(self, label: int) -> bool:
+        """Return True if label has been verified."""
+
+    @abc.abstractmethod
+    def add(self, start: int, ends: Tuple[int, ...], rule: AbstractRule) -> None:
+        """
+        Add a rule to the database.
+
+        - start is a single integer.
+        - ends is a tuple of integers, representing the children.
+        - rule is a Rule that creates start -> ends.
+        """
+
+    @abc.abstractmethod
+    def status(self, elaborate: bool) -> str:
+        """Return a string describing the status of the rule database."""
+
+    @abc.abstractmethod
+    def has_specification(self) -> bool:
+        """Indicate if the rule db contains a specification for the given label."""
+
+    @abc.abstractmethod
+    def get_specification_rules(self, **kwargs) -> Iterator[AbstractRule]:
+        """
+        Return the rules that constitute a specification.
+
+        Raise a SpecificationNotFound error if the ruledb does not contain a
+        specification.
+
+        The kwargs should to used to pass some option for the specification search that
+        are specific to each type of ruledb.
+        """
+
+
+def ensure_specification(func: RuleDBMethod) -> RuleDBMethod:
+    def inner(ruledb: RuleDBAbstract, *args, **kwargs):
+        assert isinstance(ruledb, RuleDBAbstract)
+        if not ruledb.has_specification():
+            raise SpecificationNotFound
+        return func(ruledb, *args, **kwargs)
+
+    return cast(RuleDBMethod, inner)
diff --git a/comb_spec_searcher/rule_db/base.py b/comb_spec_searcher/rule_db/base.py
index b7b146ca2..befbd5776 100644
--- a/comb_spec_searcher/rule_db/base.py
+++ b/comb_spec_searcher/rule_db/base.py
@@ -1,15 +1,19 @@
 """
-A database for rules.
+A database to search for tree.
 """
 import abc
+import itertools
 from collections import defaultdict
 from typing import Dict, Iterable, Iterator, List, MutableMapping, Optional, Set, Tuple
 
+import sympy
 import tabulate
 from logzero import logger
 
 from comb_spec_searcher.equiv_db import EquivalenceDB
-from comb_spec_searcher.exception import SpecificationNotFound
+from comb_spec_searcher.exception import InvalidOperationError
+from comb_spec_searcher.rule_db.abstract import RuleDBAbstract, ensure_specification
+from comb_spec_searcher.specification_extrator import SpecificationRuleExtractor
 from comb_spec_searcher.strategies import AbstractStrategy, VerificationRule
 from comb_spec_searcher.strategies.rule import AbstractRule
 from comb_spec_searcher.tree_searcher import (
@@ -18,19 +22,50 @@
     iterative_prune,
     proof_tree_generator_dfs,
     prune,
-    random_proof_tree,
     smallish_random_proof_tree,
 )
-from comb_spec_searcher.typing import RuleKey
+from comb_spec_searcher.typing import RuleKey, RulesDict
 
 __all__ = ["RuleDBBase", "RuleDB"]
 
 
-class RuleDBBase(abc.ABC):
-    """A database for rules found."""
+class RuleDBBase(RuleDBAbstract):
+    """
+    A database to find specification using the pruning method.
+    """
 
     def __init__(self) -> None:
+        super().__init__()
         self.equivdb = EquivalenceDB()
+        # Store the pruned dict. Should be set back to None every time the ruledb is
+        # edited. Use the propertu for clean access.
+        self._pruned_dict: Optional[RulesDict] = None
+
+    @property
+    def iterative(self) -> bool:
+        """
+        Indicate whether we are searching for an iterative specification.
+        """
+        return self.strategy_pack.iterative
+
+    @property
+    def pruned_dict(self) -> RulesDict:
+        """
+        Returned the prune dict of all the rules that are in a specification.
+
+        The dict needs to be recomputed every time the database changes and it
+        is costly to so. Call with moderation.
+        """
+        if self._pruned_dict is None:
+            rules_dict = self.rules_up_to_equivalence()
+            if self.iterative:
+                rules_dict = iterative_prune(rules_dict, root=self.root_label)
+            else:
+                prune(rules_dict)
+            self._pruned_dict = rules_dict
+            for ver_label in rules_dict.keys():
+                self.equivdb.set_verified(ver_label)
+        return self._pruned_dict
 
     @property
     @abc.abstractmethod
@@ -52,16 +87,12 @@ def __eq__(self, other: object) -> bool:
         )
 
     def add(self, start: int, ends: Tuple[int, ...], rule: AbstractRule) -> None:
-        """
-        Add a rule to the database.
-
-        - start is a single integer.
-        - ends is a tuple of integers, representing the non-empty children.
-        - rule is a Rule that creates start -> ends.
-        """
-        ends = tuple(sorted(ends))
+        self._pruned_dict = None
+        ends = self._clean_labels(ends, rule)
+        if ends == [start]:
+            return
         if isinstance(rule, VerificationRule):
-            self.set_verified(start)
+            self.equivdb.set_verified(start)
         if len(ends) == 1:
             if rule.is_two_way():
                 self.equivdb.add_two_way_edge(start, ends[0])
@@ -74,21 +105,28 @@ def add(self, start: int, ends: Tuple[int, ...], rule: AbstractRule) -> None:
         else:
             self.rule_to_strategy[(start, ends)] = rule.strategy
 
+    def _clean_labels(
+        self, ends: Tuple[int, ...], rule: AbstractRule
+    ) -> Tuple[int, ...]:
+        """
+        Remove the empty label and sort the remaining one.
+        """
+        cleaned_ends = []
+        for comb_class, child_label in zip(rule.children, ends):
+            if rule.possibly_empty and self.classdb.is_empty(comb_class, child_label):
+                logger.debug("Label %s is empty.", child_label)
+                self.searcher.classqueue.set_stop_yielding(child_label)
+                continue
+            cleaned_ends.append(child_label)
+        return tuple(sorted(cleaned_ends))
+
     def is_verified(self, label: int) -> bool:
         """Return True if label has been verified."""
-        return bool(self.equivdb.is_verified(label))
-
-    def set_verified(self, label: int) -> None:
-        """Mark label as verified."""
-        self.equivdb.set_verified(label)
+        return self.equivdb.is_verified(label)
 
     def are_equivalent(self, label: int, other: int) -> bool:
         """Return true if label and other are equivalent."""
-        return bool(self.equivdb.equivalent(label, other))
-
-    def set_equivalent(self, label: int, other: int) -> None:
-        """Mark label and other as equivalent."""
-        self.equivdb.add_two_way_edge(label, other)
+        return self.equivdb.equivalent(label, other)
 
     def rules_up_to_equivalence(self) -> Dict[int, Set[Tuple[int, ...]]]:
         """Return a defaultdict containing all rules up to the equivalence."""
@@ -102,23 +140,14 @@ def rules_up_to_equivalence(self) -> Dict[int, Set[Tuple[int, ...]]]:
             )
         return rules_dict
 
-    def all_rules(self) -> Iterator[Tuple[int, Tuple[int, ...], AbstractStrategy]]:
-        """Yield all the rules found so far."""
-        for start, ends in self:
-            yield start, ends, self.rule_to_strategy[(start, ends)]
-
     def __iter__(self) -> Iterator[Tuple[int, Tuple[int, ...]]]:
-        """Iterate through rules as the pairs (start, end)."""
-        for start, ends in self.rule_to_strategy:
-            yield start, ends
+        """Iterate through rules as the pairs (start, ends)."""
+        return itertools.chain(self.rule_to_strategy, self.eqv_rule_to_strategy)
 
     def contains(self, start: int, ends: Tuple[int, ...]) -> bool:
         """Return true if the rule start -> ends is in the database."""
-        ends = tuple(sorted(ends))
-        return (start, ends) in self.rule_to_strategy or (
-            start,
-            ends,
-        ) in self.eqv_rule_to_strategy
+        key = (start, sorted(ends))
+        return key in self.rule_to_strategy or key in self.eqv_rule_to_strategy
 
     def status(self, elaborate: bool) -> str:
         """Return a string describing the status of the rule database."""
@@ -178,9 +207,9 @@ def status(self, elaborate: bool) -> str:
     # Below are methods for finding a combinatorial specification. #
     ################################################################
 
-    def has_specification(self, label: int) -> bool:
+    def has_specification(self) -> bool:
         """Return True if a specification has been found, false otherwise."""
-        return self.is_verified(label)
+        return self.equivdb[self.root_label] in self.pruned_dict
 
     def rule_from_equivalence_rule(
         self, eqv_start: int, eqv_ends: Iterable[int]
@@ -201,11 +230,12 @@ def rule_from_equivalence_rule(
         return None
 
     def rule_from_equivalence_rule_dict(
-        self, eqv_rules: List[RuleKey]
+        self, eqv_rules: Iterable[RuleKey]
     ) -> Dict[RuleKey, RuleKey]:
         """
         Return a dictionary pointing from an equivalence rule to an actual rule.
         """
+        eqv_rules = set(eqv_rules)
         res: Dict[RuleKey, RuleKey] = {}
         for start, ends in self.rule_to_strategy:
             eqv_start = self.equivdb[start]
@@ -215,94 +245,108 @@ def rule_from_equivalence_rule_dict(
                 res[eqv_key] = (start, ends)
         return res
 
-    def find_specification(
-        self, label: int, minimization_time_limit: float, iterative: bool = False
+    @ensure_specification
+    def get_specification_rules(  # pylint: disable=arguments-differ
+        self,
+        *,
+        minimization_time_limit: float = 10,
+        smallest: bool = False,
+        **kwargs,
+    ) -> Iterator[AbstractRule]:
+        node = self._get_specification_node(minimization_time_limit, smallest)
+        logger.info("Found specification with %s rules.", len(node.labels()))
+        spec_extractor = SpecificationRuleExtractor(
+            self.root_label, node, self, self.classdb
+        )
+        return spec_extractor.rules()
+
+    @ensure_specification
+    def _get_specification_node(
+        self, minimization_time_limit: float, smallest: bool
     ) -> Node:
-        """Search for a specification based on current data found."""
-        rules_dict = self.rules_up_to_equivalence()
-        # Prune all unverified labels (recursively)
-        if iterative:
-            rules_dict = iterative_prune(rules_dict, root=label)
+        """
+        Return a specification node for the given label.
+        """
+        if self.iterative:
+            if smallest:
+                raise InvalidOperationError("can't use iterative and smallest")
+            node = self._get_iterative_node()
         else:
-            prune(rules_dict)  # this function removes rules not in a specification.
-
-        # only verified labels in rules_dict, in particular, there is a
-        # specification if a label is in the rules_dict
-        for ver_label in rules_dict.keys():
-            self.set_verified(ver_label)
-
-        if self.equivdb[label] in rules_dict:
-            logger.info("Specification detected.")
-            if iterative:
-                specification = iterative_proof_tree_finder(
-                    rules_dict, root=self.equivdb[label]
-                )
+            if smallest:
+                node = self._get_smallest_node(minimization_time_limit)
             else:
-                logger.info(
-                    "Minimizing for %s seconds.", round(minimization_time_limit)
-                )
-                specification = smallish_random_proof_tree(
-                    rules_dict, self.equivdb[label], minimization_time_limit
+                node = self._get_smallish_node(
+                    minimization_time_limit=minimization_time_limit,
                 )
-            logger.info(
-                "Found specification with %s rules.", len(specification.labels())
-            )
-        else:
-            raise SpecificationNotFound("No specification for label {}".format(label))
-        return specification
+        return node
+
+    @ensure_specification
+    def _get_iterative_node(self) -> Node:
+        if not self.iterative:
+            raise InvalidOperationError("Not supported in non-iterative mode")
+        return iterative_proof_tree_finder(
+            self.pruned_dict, root=self.equivdb[self.root_label]
+        )
 
-    def all_specifications(self, label: int, iterative: bool = False) -> Iterator[Node]:
-        """
-        A generator that yields all specifications in the universe for
-        the given label.
+    @ensure_specification
+    def _get_smallish_node(self, minimization_time_limit: float) -> Node:
+        """Search for a smallish node for the given minimization time."""
+        if self.iterative:
+            raise InvalidOperationError("Not supported in iterative mode.")
+        logger.info("Minimizing for %s seconds.", round(minimization_time_limit))
+        return smallish_random_proof_tree(
+            self.pruned_dict, self.equivdb[self.root_label], minimization_time_limit
+        )
+
+    @ensure_specification
+    def _get_smallest_node(self, minimization_time_limit: float) -> Node:
         """
-        if iterative:
-            raise NotImplementedError(
-                "There is no method for yielding all iterative proof trees."
-            )
-        rules_dict = self.rules_up_to_equivalence()
-        # Prune all unverified labels (recursively)
-        prune(rules_dict)
+        Return the smallest specification node in the universe.
 
-        if self.equivdb[label] in rules_dict:
-            yield from proof_tree_generator_dfs(rules_dict, root=self.equivdb[label])
+        Will look for a smallish spec for minimization_time_limit and then use
+        binary search to find it to find the smallest one.
 
-    def get_smallest_specification(self, label: int, iterative: bool = False) -> Node:
-        """
-        Return the smallest specification in the universe for label. It uses
-        exponential search to find it.
         This doesn't consider the length of the equivalence paths.
         """
-        if iterative:
-            raise NotImplementedError(
-                "There is no method for finding smallest iterative proof trees."
-            )
-        rules_dict = self.rules_up_to_equivalence()
-        prune(rules_dict)
-
-        if not self.equivdb[label] in rules_dict:
-            raise SpecificationNotFound("No specification for label {}".format(label))
-        tree = random_proof_tree(rules_dict, root=self.equivdb[label])
+        if self.iterative:
+            raise InvalidOperationError("Not supported in iterative mode.")
+        node = self._get_smallish_node(minimization_time_limit)
         minimum = 1
-        maximum = len(tree)
+        maximum = len(node)
         logger.info(
-            "Found a specification of size %s. Looking for the smallest.", len(tree)
+            "Found a specification of size %s. Looking for the smallest.", len(node)
         )
         # Binary search to find a smallest proof tree.
         while minimum < maximum:
             middle = (minimum + maximum) // 2
             logger.info("Looking for specification of size %s", middle)
             try:
-                tree = next(
+                node = next(
                     proof_tree_generator_dfs(
-                        rules_dict, root=self.equivdb[label], maximum=middle
+                        self.pruned_dict,
+                        root=self.equivdb[self.root_label],
+                        maximum=middle,
                     )
                 )
-                maximum = min(middle, len(tree))
+                maximum = min(middle, len(node))
             except StopIteration:
                 minimum = middle + 1
-        logger.info("The smallest specification is of size %s.", len(tree))
-        return tree
+        logger.info("The smallest specification is of size %s.", len(node))
+        return node
+
+    @ensure_specification
+    def _all_nodes(self, iterative: bool = False) -> Iterator[Node]:
+        """
+        A generator that yields all specifications in the universe for
+        the given label.
+        """
+        if self.iterative:
+            raise NotImplementedError(
+                "There is no method for yielding all iterative proof trees."
+            )
+        yield from proof_tree_generator_dfs(
+            self.pruned_dict, root=self.equivdb[self.root_label]
+        )
 
 
 class RuleDB(RuleDBBase):
@@ -322,3 +366,38 @@ def eqv_rule_to_strategy(
         self,
     ) -> Dict[Tuple[int, Tuple[int, ...]], AbstractStrategy]:
         return self._eqv_rule_to_strategy
+
+    def all_rules(self) -> Iterator[AbstractRule]:
+        """Yield all the rules found so far."""
+        for rulekey in self:
+            parent = self.classdb.get_class(rulekey[0])
+            try:
+                strategy = self.rule_to_strategy[rulekey]
+            except KeyError:
+                strategy = self.eqv_rule_to_strategy[rulekey]
+            yield strategy(parent)
+
+    def all_equations(self) -> Set[sympy.Eq]:
+        """
+        Returns a set of equations for all rules currently found.
+        """
+
+        def get_function(comb_class) -> sympy.Function:
+            return comb_class.get_function(self.classdb.get_label)
+
+        eqs = set()
+        for rule in self.all_rules():
+            try:
+                eq = rule.get_equation(get_function)
+            except NotImplementedError:
+                logger.info(
+                    "can't find generating function for %s." " The comb class is:\n%s",
+                    get_function(rule.comb_class),
+                    rule.comb_class,
+                )
+                eq = sympy.Eq(
+                    get_function(rule.comb_class),
+                    sympy.Function("NOTIMPLEMENTED")(sympy.var("x")),
+                )
+            eqs.add(eq)
+        return eqs
diff --git a/comb_spec_searcher/rule_db/forest.py b/comb_spec_searcher/rule_db/forest.py
new file mode 100644
index 000000000..470cd58ec
--- /dev/null
+++ b/comb_spec_searcher/rule_db/forest.py
@@ -0,0 +1,657 @@
+import itertools
+import time
+from datetime import timedelta
+from typing import (
+    Callable,
+    Deque,
+    Dict,
+    Generic,
+    Iterable,
+    Iterator,
+    List,
+    Optional,
+    Set,
+    Tuple,
+    TypeVar,
+    Union,
+)
+
+from logzero import logger
+
+from comb_spec_searcher.class_db import ClassDB
+from comb_spec_searcher.exception import StrategyDoesNotApply
+from comb_spec_searcher.rule_db.abstract import RuleDBAbstract, ensure_specification
+from comb_spec_searcher.strategies.rule import AbstractRule, Rule
+from comb_spec_searcher.strategies.strategy import (
+    AbstractStrategy,
+    EmptyStrategy,
+    StrategyFactory,
+)
+from comb_spec_searcher.strategies.strategy_pack import StrategyPack
+from comb_spec_searcher.typing import CSSstrategy, ForestRuleKey, RuleBucket, RuleKey
+
+T = TypeVar("T")
+RuleWithShifts = Tuple[RuleKey, Tuple[int, ...]]
+SortedRWS = Dict[RuleBucket, List[ForestRuleKey]]
+empty_strategy: EmptyStrategy = EmptyStrategy()
+
+
+class DefaultList(Generic[T]):
+    """
+    A list data structure get automatically gets longer if an index not existing is
+    requested.
+
+    When getting longer the list is filled by calling the provied `default_factory`.
+    This is similar to the `collections.defaultdict` data structure.
+    """
+
+    def __init__(self, default_factory: Callable[[], T]):
+        self._default_factory = default_factory
+        self._list: List[T] = []
+
+    def _increase_list_len(self, key: int) -> None:
+        """
+        Increase the length of the list so that the given is valid.
+        """
+        num_new_entry = key - len(self._list) + 1
+        self._list.extend((self._default_factory() for _ in range(num_new_entry)))
+
+    def __getitem__(self, key: int) -> T:
+        try:
+            return self._list[key]
+        except IndexError:
+            self._increase_list_len(key)
+            return self._list[key]
+
+    def __setitem__(self, key: int, value: T) -> None:
+        self._list[key] = value
+
+    def __iter__(self) -> Iterator[T]:
+        return iter(self._list)
+
+    def __str__(self) -> str:
+        return str(self._list)
+
+
+class Function:
+    """
+    A python representation of a function.
+
+    The function maps natural number to a natural number or infinity (represented
+    by the use of None)
+
+    The default value of the function is 0.
+    """
+
+    def __init__(self) -> None:
+        self._value: List[Optional[int]] = []
+        self._preimage_count: DefaultList[int] = DefaultList(int)
+        self._infinity_count: int = 0
+
+    @property
+    def preimage_count(self) -> List[int]:
+        """
+        Return the number of classes for each value of the function.
+        """
+        count = list(self._preimage_count)
+        while count and count[-1] == 0:
+            count.pop()
+        return count
+
+    @property
+    def infinity_count(self) -> int:
+        return self._infinity_count
+
+    def __getitem__(self, key: int) -> Optional[int]:
+        """
+        Return the value of the function for the given key.
+        """
+        try:
+            return self._value[key]
+        except IndexError:
+            self._increase_list_len(key)
+            return 0
+
+    def _increase_list_len(self, key: int) -> None:
+        num_new_entry = key - len(self._value) + 1
+        self._value.extend((0 for _ in range(num_new_entry)))
+        self._preimage_count[0] += num_new_entry
+
+    def increase_value(self, key: int) -> None:
+        """
+        Increase by one the value of the function for the given key.
+        """
+        try:
+            old_value = self._value[key]
+            if old_value is None:
+                raise ValueError(f"The function is already infinity for {key}")
+        except IndexError:
+            self._increase_list_len(key)
+            old_value = 0
+        self._value[key] = old_value + 1
+        self._preimage_count[old_value] -= 1
+        self._preimage_count[old_value + 1] += 1
+
+    def set_infinite(self, key: int) -> None:
+        """
+        Set the value of function for the given key to infinity.
+        """
+        try:
+            old_value = self._value[key]
+            if old_value is None:
+                raise ValueError(f"The function is already infinity for {key}")
+        except IndexError:
+            self._increase_list_len(key)
+            old_value = 0
+        self._value[key] = None
+        self._preimage_count[old_value] -= 1
+        self._infinity_count += 1
+
+    def preimage_gap(self, length: int) -> int:
+        """
+        Return the smallest k such that the preimage of the interval
+        [k, k+length-1] is empty.
+        """
+        if length <= 0:
+            raise ValueError("length argument must be positive")
+        last_non_zero = -1
+        for i, v in enumerate(self._preimage_count):
+            if v != 0:
+                last_non_zero = i
+            elif i - last_non_zero >= length:
+                return last_non_zero + 1
+        return last_non_zero + 1
+
+    def preimage(self, value: Optional[int]) -> Iterator[int]:
+        """
+        Return the preimage of the function for the given value.
+        """
+        if value == 0:
+            raise ValueError("The preimage of 0 is infinite.")
+        return (k for k, v in enumerate(self._value) if v == value)
+
+    def to_dict(self) -> Dict[int, Optional[int]]:
+        """
+        Return a dictionary view of the function with only the non-zero value.
+        """
+        return {i: v for i, v in enumerate(self._value) if v != 0}
+
+    def __str__(self) -> str:
+        parts = (
+            f"{i} -> {v if v is not None else '∞'}" for i, v in enumerate(self._value)
+        )
+        return "\n".join(parts)
+
+
+class TableMethod:
+    def __init__(self) -> None:
+        self._rules: List[ForestRuleKey] = []
+        self._shifts: List[List[Optional[int]]] = []
+        self._function: Function = Function()
+        self._gap_size: int = 1
+        self._rules_using_class: DefaultList[List[Tuple[int, int]]] = DefaultList(list)
+        self._rules_pumping_class: DefaultList[List[int]] = DefaultList(list)
+        self._processing_queue: Deque[int] = Deque()
+        self._current_gap: Tuple[int, int] = (1, 1)
+        self._rule_holding_extra_terms: Set[int] = set()
+
+    @property
+    def function(self) -> Dict[int, Optional[int]]:
+        """
+        Return a dict representing the number of term that are computable for each
+        class.
+        Only the class where it can get at least a term are included. If a class is map
+        to None, then all the terms of the enumeration are computable.
+        """
+        return self._function.to_dict()
+
+    def add_rule_key(
+        self,
+        rule_key: ForestRuleKey,
+    ):
+        """
+        Add the rule to the database.
+
+        INPUTS:
+          - `rule_key`
+          - `shifts_for_zero`: The values of the shifts if no information was known
+          about any of the classes.
+          - `rule_bucket` the type of rule
+        """
+        self._rules.append(rule_key)
+        self._shifts.append(self._compute_shift(rule_key.key, rule_key.shifts))
+        max_gap = max((abs(s) for s in rule_key.shifts), default=0)
+        if max_gap > self._gap_size:
+            self._gap_size = max_gap
+            self._correct_gap()
+        if self._function[rule_key.parent] is not None:
+            rule_idx = len(self._rules) - 1
+            self._rules_pumping_class[rule_key.parent].append(rule_idx)
+            for child_idx, child in enumerate(rule_key.children):
+                if self._function[child] is not None:
+                    self._rules_using_class[child].append((rule_idx, child_idx))
+            self._processing_queue.append(rule_idx)
+        self._process_queue()
+
+    def is_pumping(self, label: int) -> bool:
+        """
+        Determine if the comb_class is pumping in the current universe.
+        """
+        return self._function[label] is None
+
+    def status(self) -> str:
+        s = f"\tSize of the gap: {self._gap_size}\n"
+        s += f"\tSize of the stable subset: {self._function.infinity_count}\n"
+        s += f"\tSizes of the pre-images: {self._function.preimage_count}\n"
+        return s
+
+    def stable_subset(self) -> Iterator[int]:
+        return self._function.preimage(None)
+
+    def pumping_subuniverse(
+        self,
+    ) -> Iterator[ForestRuleKey]:
+        """
+        Iterator over all the forest rule keys that contain only pumping
+        combinatorial classes.
+        """
+        stable_subset = set(self.stable_subset())
+        for forest_key in self._rules:
+            if forest_key.parent in stable_subset and stable_subset.issuperset(
+                forest_key.children
+            ):
+                yield forest_key
+
+    def _compute_shift(
+        self,
+        rule_key: RuleKey,
+        shifts_for_zero: Tuple[int, ...],
+    ) -> List[Optional[int]]:
+        """
+        Compute the initial value for the shifts a rule based on the current state of
+        the function.
+        """
+        parent_current_value = self._function[rule_key[0]]
+        if parent_current_value is None:
+            return [None for _ in shifts_for_zero]
+        chidlren_function_value = map(self._function.__getitem__, rule_key[1])
+        return [
+            fvalue + sfz - parent_current_value if fvalue is not None else None
+            for fvalue, sfz in zip(chidlren_function_value, shifts_for_zero)
+        ]
+
+    def _correct_gap(self) -> None:
+        """
+        Correct the gap and if needed queue rules for the classes that were previously
+        on the right hand  side of the gap.
+
+        This should be toggled every time the gap changes whether the size changes or
+        the some value changes of the function caused the gap to change.
+        """
+        k = self._function.preimage_gap(self._gap_size)
+        new_gap = (k, k + self._gap_size - 1)
+        if new_gap[1] > self._current_gap[1]:
+            self._processing_queue.extend(self._rule_holding_extra_terms)
+            self._rule_holding_extra_terms.clear()
+        self._current_gap = new_gap
+
+    def _process_queue(self) -> None:
+        """
+        Try to make improvement with all the class in the processing queue.
+        """
+        while self._processing_queue or self._rule_holding_extra_terms:
+            while self._processing_queue:
+                rule_idx = self._processing_queue.popleft()
+                shifts = self._shifts[rule_idx]
+                if self._can_give_terms(shifts):
+                    parent = self._rules[rule_idx].parent
+                    self._increase_value(parent, rule_idx)
+            if self._rule_holding_extra_terms:
+                rule_idx = self._rule_holding_extra_terms.pop()
+                parent = self._rules[rule_idx].parent
+                self._set_infinite(parent)
+
+    @staticmethod
+    def _can_give_terms(shifts: List[Optional[int]]) -> bool:
+        """
+        Return True if the shifts indicate that a new terms can be computed.
+        """
+        return all(s is None or s > 0 for s in shifts)
+
+    def _increase_value(self, comb_class: int, rule_idx: int) -> None:
+        """
+        Increase the value of the comb_class and put on the processing stack any rule
+        that can now give a new term.
+
+        The rule_idx must indicate the rule used to justify the increase.
+        """
+        current_value = self._function[comb_class]
+        if current_value is None:
+            return
+        if current_value > self._current_gap[1]:
+            self._rule_holding_extra_terms.add(rule_idx)
+            return
+        self._function.increase_value(comb_class)
+        # Correction of the gap
+        gap_start = self._function.preimage_gap(self._gap_size)
+        if self._current_gap[0] != gap_start:
+            self._correct_gap()
+        # Correction of the shifts for rule pumping comb_class
+        for r_idx in self._rules_pumping_class[comb_class]:
+            shifts = self._shifts[r_idx]
+            for i, v in enumerate(shifts):
+                shifts[i] = v - 1 if v is not None else None
+            if self._can_give_terms(shifts):
+                self._processing_queue.append(r_idx)
+        # Correction of the shifts for rules using comb_class to pump
+        for r_idx, class_idx in self._rules_using_class[comb_class]:
+            shifts = self._shifts[r_idx]
+            current_shift = shifts[class_idx]
+            assert current_shift is not None
+            shifts[class_idx] = current_shift + 1
+            if self._can_give_terms(shifts):
+                self._processing_queue.append(r_idx)
+
+    def _set_infinite(self, comb_class: int) -> None:
+        """
+        Set the value if the class to infinity.
+
+        This should happen when we know that we cannot pump anything on the left side
+        of the gap.
+        """
+        current_value = self._function[comb_class]
+        if current_value is None:
+            return
+        assert current_value > self._current_gap[1]
+        assert not self._processing_queue
+        self._function.set_infinite(comb_class)
+        # This class will never be increased again so we remove any occurrence
+        # of the rule of any rule for that class from _rules_using_class and
+        # _rules_pumping_class
+        for rule_idx in self._rules_pumping_class[comb_class]:
+            for child in self._rules[rule_idx].children:
+                self._rules_using_class[child] = [
+                    (ri, ci)
+                    for ri, ci in self._rules_using_class[child]
+                    if ri != rule_idx
+                ]
+        self._rules_pumping_class[comb_class].clear()
+        # Correction of the shifts for rules using comb_class to pump
+        for rule_idx, class_idx in self._rules_using_class[comb_class]:
+            shifts = self._shifts[rule_idx]
+            shifts[class_idx] = None
+            if self._can_give_terms(shifts):
+                self._processing_queue.append(rule_idx)
+        self._rules_using_class[comb_class].clear()
+
+    def rule_info(self, rule_idx: int) -> str:
+        """
+        Return a string with information about a particular rule.
+        Mostly intended for debugging.
+        """
+
+        def v_to_str(v: Optional[int]) -> str:
+            """Return a string for the integer and infinity if None"""
+            if v is None:
+                return "∞"
+            return str(v)
+
+        rule_key = self._rules[rule_idx]
+        current_value = f"{v_to_str(self._function[rule_key.parent])} -> " + ", ".join(
+            map(v_to_str, (self._function[c] for c in rule_key.children))
+        )
+        shifts = map(v_to_str, self._shifts[rule_idx])
+        child_with_shift = ", ".join(
+            f"({c}, {s})" for c, s in zip(rule_key.children, shifts)
+        )
+        return f"{rule_key.parent} -> {child_with_shift} || {current_value}"
+
+
+class ForestRuleExtractor:
+    MINIMIZE_ORDER = (
+        RuleBucket.REVERSE,
+        RuleBucket.NORMAL,
+        RuleBucket.EQUIV,
+        RuleBucket.VERIFICATION,
+    )
+
+    def __init__(
+        self,
+        root_label: int,
+        ruledb: "RuleDBForest",
+        classdb: ClassDB,
+        pack: StrategyPack,
+    ):
+        self.pack = pack
+        self.classdb = classdb
+        self.root_label = root_label
+        self.rule_by_bucket = self._sorted_stable_rules(ruledb.table_method)
+        assert set(ForestRuleExtractor.MINIMIZE_ORDER) == set(self.rule_by_bucket)
+        self.needed_rules: List[ForestRuleKey] = list()
+        self._minimize()
+
+    def check(self) -> bool:
+        """
+        Make a sanity check of the status of the extractor.
+        """
+        lhs = set(rk.parent for rk in self.needed_rules)
+        assert len(lhs) == len(self.needed_rules)
+        assert self._is_productive(self.needed_rules)
+
+    def rules(self) -> Iterator[AbstractRule]:
+        """
+        Return all the rules of the specification.
+
+        The empty rule are ignored as they be produced as needed by the specification.
+        """
+        for rk in self.needed_rules:
+            rule = self._find_rule(rk)
+            if isinstance(rule.strategy, EmptyStrategy):
+                continue
+            if rule.is_equivalence():
+                assert isinstance(rule, Rule)
+                yield rule.to_equivalence_rule()
+            else:
+                yield rule
+
+    def _minimize(self):
+        """
+        Perform the complete minimization of the forest.
+        """
+        for key in ForestRuleExtractor.MINIMIZE_ORDER:
+            self._minimize_key(key)
+
+    def _minimize_key(self, key: RuleBucket) -> None:
+        """
+        Minimize the number of rules used for the type of rule given by key.
+
+        The list of rule in `self.rule_by_bucket[key]` is cleared and a
+        minimal set from theses is added to `self.needed_rules`.
+        """
+        logger.info("Minimizing %s", key.name)
+        maybe_useful: List[ForestRuleKey] = []
+        not_minimizing: List[List[ForestRuleKey]] = [
+            self.needed_rules,
+            maybe_useful,
+        ]
+        not_minimizing.extend(
+            rules for k, rules in self.rule_by_bucket.items() if k != key
+        )
+        minimizing = self.rule_by_bucket[key]
+        while minimizing:
+            tb = TableMethod()
+            # Add the rule we are not trying to minimize
+            for rk in itertools.chain.from_iterable(not_minimizing):
+                tb.add_rule_key(rk)
+            if tb.is_pumping(self.root_label):
+                minimizing.clear()
+                break
+            # Add rule until it gets productive
+            for i, rk in enumerate(minimizing):
+                tb.add_rule_key(rk)
+                if tb.is_pumping(self.root_label):
+                    break
+            else:
+                raise RuntimeError("Not pumping after adding all rules")
+            maybe_useful.append(rk)
+            assert minimizing, "variable i won't be set"
+            # pylint: disable=undefined-loop-variable
+            for _ in range(i, len(minimizing)):
+                minimizing.pop()
+        counter = 0
+        while maybe_useful:
+            rk = maybe_useful.pop()
+            if not self._is_productive(itertools.chain.from_iterable(not_minimizing)):
+                self.needed_rules.append(rk)
+                counter += 1
+        logger.info("Using %s rule for %s", counter, key.name)
+
+    def _is_productive(self, rule_keys: Iterable[ForestRuleKey]) -> bool:
+        """
+        Check if the given set of rules is productive.
+        """
+        ruledb = TableMethod()
+        for rk in rule_keys:
+            ruledb.add_rule_key(rk)
+        return ruledb.is_pumping(self.root_label)
+
+    def _sorted_stable_rules(self, ruledb: TableMethod) -> SortedRWS:
+        """
+        Extract all the rule from the stable subuniverse and return all of them in a
+        dict sorted by type.
+        """
+        res: SortedRWS = {bucket: [] for bucket in self.MINIMIZE_ORDER}
+        for forest_key in ruledb.pumping_subuniverse():
+            try:
+                res[forest_key.bucket].append(forest_key)
+            except KeyError as e:
+                msg = (
+                    f"{forest_key.bucket} type is not currently supported "
+                    "by the extractor"
+                )
+                raise RuntimeError(msg) from e
+        return res
+
+    def _find_rule(self, rule_key: ForestRuleKey) -> AbstractRule:
+        """
+        Find a rule that have the given rule key.
+        """
+        all_classes = (rule_key.parent,) + rule_key.children
+        all_normal_rules = itertools.chain.from_iterable(
+            self._rules_for_class(c) for c in all_classes
+        )
+        for normal_rule in all_normal_rules:
+            potential_rules = [normal_rule]
+            if normal_rule.is_reversible():
+                assert isinstance(normal_rule, Rule)
+                potential_rules.extend(
+                    normal_rule.to_reverse_rule(i)
+                    for i in range(len(normal_rule.children))
+                )
+            for rule in potential_rules:
+                if rule.forest_key(self.classdb.get_label) == rule_key:
+                    return rule
+        err = f"Can't find a rule for {rule_key}\n"
+        err += f"Parent:\n{self.classdb.get_class(rule_key.parent)}\n"
+        for i, l in enumerate(rule_key.children):
+            err += f"Child {i}:\n{self.classdb.get_class(l)}\n"
+        raise RuntimeError(err)
+
+    def _rules_for_class(self, label: int) -> Iterator[AbstractRule]:
+        """
+        Return all the rule created for that class with the pack.
+        """
+        comb_class = self.classdb.get_class(label)
+        strats: Iterator[CSSstrategy] = itertools.chain([EmptyStrategy()], self.pack)
+        for strat in strats:
+            if isinstance(strat, StrategyFactory):
+                strats_or_rules: Iterable[
+                    Union[AbstractRule, AbstractStrategy]
+                ] = strat(comb_class)
+            else:
+                strats_or_rules = [strat]
+            for x in strats_or_rules:
+                if isinstance(x, AbstractStrategy):
+                    try:
+                        yield x(comb_class)
+                    except StrategyDoesNotApply:
+                        continue
+                else:
+                    yield x
+
+
+class RuleDBForest(RuleDBAbstract):
+    """
+    The rule database that provides live information on which class are pumping with the
+    current rule in the database.
+    """
+
+    def __init__(self) -> None:
+        super().__init__()
+        self._num_rules = 0
+        self._time_table_method = 0.0
+        self._time_key = 0.0
+        self.table_method = TableMethod()
+        self._already_empty: Set[int] = set()
+
+    # Implementation of RuleDBAbstract
+
+    def status(self, elaborate: bool) -> str:
+        s = "RuleDB status:\n"
+        s += f"\tAdded from {self._num_rules} normal rules\n"
+        key_time = timedelta(seconds=int(self._time_key))
+        tm_time = timedelta(seconds=int(self._time_table_method))
+        s += f"\tTime spent computing forest keys: {key_time}\n"
+        s += f"\tTime spent running the table method: {tm_time}\n"
+        s += self.table_method.status()
+        return s
+
+    def is_verified(self, label: int) -> bool:
+        """
+        Determine if the comb_class is pumping in the current universe.
+        """
+        return self.table_method.is_pumping(label)
+
+    def has_specification(self) -> bool:
+        return self.is_verified(self.root_label)
+
+    def add(self, start: int, ends: Tuple[int, ...], rule: AbstractRule) -> None:
+        self._add_empty_rule(ends, rule)
+        self._num_rules += 1
+        start_time = time.time()
+        new_rule_keys = [rule.forest_key(self.classdb.get_label)]
+        if rule.is_reversible():
+            assert isinstance(rule, Rule)
+            new_rule_keys.extend(
+                rule.to_reverse_rule(i).forest_key(self.classdb.get_label)
+                for i in range(len(rule.children))
+            )
+        self._time_key += time.time() - start_time
+        start_time = time.time()
+        for new_key in new_rule_keys:
+            self.table_method.add_rule_key(new_key)
+        self._time_table_method += time.time() - start_time
+
+    @ensure_specification
+    def get_specification_rules(self, **kwargs) -> Iterator[AbstractRule]:
+        extractor = ForestRuleExtractor(
+            self.root_label, self, self.classdb, self.strategy_pack
+        )
+        extractor.check()
+        return extractor.rules()
+
+    # Other methods
+
+    def _add_empty_rule(self, ends: Iterable[int], rule: AbstractRule) -> None:
+        """
+        Add empty rule for the children of the rule if needed.
+        """
+        if not rule.possibly_empty:
+            return
+        for label, comb_class in zip(ends, rule.children):
+            if label not in self._already_empty and self.classdb.is_empty(
+                comb_class, label
+            ):
+                rule = empty_strategy(comb_class)
+                self._already_empty.add(label)
+                self.searcher.add_rule(label, (), rule)
diff --git a/comb_spec_searcher/rule_db/forget.py b/comb_spec_searcher/rule_db/forget.py
index e452f65c3..1b51870bf 100644
--- a/comb_spec_searcher/rule_db/forget.py
+++ b/comb_spec_searcher/rule_db/forget.py
@@ -1,11 +1,22 @@
 """
-A database for rules.
+A database to search for tree.
+
+The database do not store the strategy to save memory.
 """
 import itertools
-from typing import Iterable, Iterator, MutableMapping, Set, Tuple, Union, cast
+from typing import (
+    TYPE_CHECKING,
+    Iterable,
+    Iterator,
+    MutableMapping,
+    Optional,
+    Set,
+    Tuple,
+    Union,
+    cast,
+)
 
 from comb_spec_searcher.class_db import ClassDB
-from comb_spec_searcher.equiv_db import EquivalenceDB
 from comb_spec_searcher.exception import StrategyDoesNotApply
 from comb_spec_searcher.strategies.rule import AbstractRule
 from comb_spec_searcher.strategies.strategy import AbstractStrategy, StrategyFactory
@@ -14,9 +25,13 @@
 
 from .base import RuleDBBase
 
+if TYPE_CHECKING:
+    from comb_spec_searcher import CombinatorialSpecificationSearcher
+
 __all__ = ["RuleDBForgetStrategy"]
 
 
+# pylint: disable=too-many-ancestors
 class RecomputingDict(MutableMapping[RuleKey, AbstractStrategy]):
     """
     A mapping from rules to strategies that recompute the strategy every time it's
@@ -30,17 +45,31 @@ class RecomputingDict(MutableMapping[RuleKey, AbstractStrategy]):
 
     def __init__(
         self,
-        classdb: ClassDB,
-        strat_pack: StrategyPack,
-        equivdb: EquivalenceDB,
         only_equiv: bool,
     ) -> None:
-        self.classdb = classdb
-        self.equivdb = equivdb
-        self.pack = strat_pack
+        self._classdb: Optional[ClassDB] = None
+        self._pack: Optional[StrategyPack] = None
         self.rules: Set[Tuple[int, ...]] = set()
         self.only_equiv: bool = only_equiv
 
+    def link_searcher(self, classdb: ClassDB, strat_pack: StrategyPack) -> None:
+        if self._classdb is not None or self._pack is not None:
+            raise RuntimeError("Searcher is alreay linked")
+        self._classdb = classdb
+        self._pack = strat_pack
+
+    @property
+    def classdb(self) -> ClassDB:
+        if self._classdb is None:
+            raise RuntimeError("Not linked with classdb")
+        return self._classdb
+
+    @property
+    def pack(self) -> StrategyPack:
+        if self._pack is None:
+            raise RuntimeError("Not linked with pack")
+        return self._pack
+
     @staticmethod
     def _flatten(tuple_: RuleKey) -> Tuple[int, ...]:
         return (tuple_[0],) + tuple_[1]
@@ -111,14 +140,17 @@ def __contains__(self, key: object) -> bool:
 
 
 class RuleDBForgetStrategy(RuleDBBase):
-    def __init__(self, classdb: ClassDB, strat_pack: StrategyPack) -> None:
+    def __init__(self) -> None:
         super().__init__()
-        self._rule_to_strategy = RecomputingDict(
-            classdb, strat_pack, self.equivdb, only_equiv=False
-        )
-        self._eqv_rule_to_strategy = RecomputingDict(
-            classdb, strat_pack, self.equivdb, only_equiv=True
-        )
+        self._rule_to_strategy = RecomputingDict(only_equiv=False)
+        self._eqv_rule_to_strategy = RecomputingDict(only_equiv=True)
+
+    def link_searcher(self, searcher: "CombinatorialSpecificationSearcher") -> None:
+        super().link_searcher(searcher)
+        classdb = searcher.classdb
+        strat_pack = searcher.strategy_pack
+        self.rule_to_strategy.link_searcher(classdb, strat_pack)
+        self.eqv_rule_to_strategy.link_searcher(classdb, strat_pack)
 
     @property
     def rule_to_strategy(self) -> RecomputingDict:
diff --git a/comb_spec_searcher/rule_db/limited_strategy.py b/comb_spec_searcher/rule_db/limited_strategy.py
deleted file mode 100644
index ba49b7d3b..000000000
--- a/comb_spec_searcher/rule_db/limited_strategy.py
+++ /dev/null
@@ -1,125 +0,0 @@
-"""
-    An alternative ruledb that permits restricting the number of applications of a given
-    set of strategies to a given limit.
-"""
-from copy import deepcopy
-from itertools import combinations
-from typing import Iterable, Set, Type
-
-from logzero import logger
-
-from comb_spec_searcher.exception import SpecificationNotFound
-from comb_spec_searcher.rule_db import RuleDB
-from comb_spec_searcher.strategies import AbstractStrategy
-from comb_spec_searcher.tree_searcher import (
-    Node,
-    iterative_proof_tree_finder,
-    iterative_prune,
-    prune,
-    smallish_random_proof_tree,
-)
-from comb_spec_searcher.typing import RuleKey
-
-
-class LimitedStrategyRuleDB(RuleDB):
-    """
-    An alternative ruledb that permits restricting the number of applications of a
-    given set of strategies to a given limit. The limit is a single limit for the
-    total number of strategies whose type is in the given set; it is not a per-type
-    limit.
-    """
-
-    def __init__(
-        self,
-        strategies_to_limit: Iterable[Type[AbstractStrategy]],
-        limit: int,
-        mark_verified: bool,
-    ) -> None:
-        self.strategies_to_limit = set(strategies_to_limit)
-        self.limit = limit
-        self.mark_verified = mark_verified
-        super().__init__()
-
-    def get_rules_up_to_equiv_using_strategies(self) -> Set[RuleKey]:
-        """
-        returns a set of rules, up to equivalence, that come from
-        <self.strategies_to_limit>
-        """
-        eqv_rules_using_strategies: Set[RuleKey] = set()
-        for label_rule in self:
-            rule_strat = self.rule_to_strategy[label_rule]
-            start, ends = label_rule
-            if isinstance(rule_strat, tuple(self.strategies_to_limit)):
-                temp_start = self.equivdb[start]
-                temp_ends = tuple(sorted(map(self.equivdb.__getitem__, ends)))
-                eqv_rules_using_strategies.add((temp_start, temp_ends))
-        return eqv_rules_using_strategies
-
-    def find_specification(
-        self, label: int, minimization_time_limit: float, iterative: bool = False
-    ) -> Node:
-        """Search for a specification based on current data found."""
-        rules_dict = self.rules_up_to_equivalence()
-        # Prune all unverified labels (recursively)
-        # Afterward, only verified labels remain in rules_dict. In particular, there
-        # is a specification if <label> is in the rules_dict
-        if iterative:
-            rules_dict = iterative_prune(rules_dict, root=label)
-        else:
-            prune(rules_dict)  # this function removes rules not in a specification.
-
-        # if the root is not in the pruned dict, then there's nothing to be done
-        if self.equivdb[label] not in rules_dict:
-            raise SpecificationNotFound("No specification for label {}".format(label))
-
-        # find the rules that we want to limit, and consider only those in a spec
-        rules_to_isolate = self.get_rules_up_to_equiv_using_strategies()
-        rules_to_isolate = set(
-            (lhs, rhs)
-            for (lhs, rhs) in rules_to_isolate
-            if lhs in rules_dict and rhs in rules_dict[lhs]
-        )
-
-        logger.debug("Found %s rules to isolate.", len(rules_to_isolate))
-
-        if len(rules_to_isolate) < self.limit:
-            rule_combos = [tuple(rules_to_isolate)]
-        else:
-            rule_combos = list(combinations(rules_to_isolate, self.limit))
-        num_combos = len(rule_combos)
-        logger.debug("%s combinations to check", num_combos)
-
-        # for each combination of <self.limit> rules, remove all but those and
-        # reprune to see if a spec still exists
-        for i, combo in enumerate(rule_combos):
-            logger.debug("[ %s / %s ]", i + 1, num_combos)
-            temp_rules_dict = deepcopy(rules_dict)
-            for lhs, rhs in rules_to_isolate:
-                if (lhs, rhs) in combo:
-                    continue
-                try:
-                    temp_rules_dict[lhs].remove(rhs)
-                except KeyError:
-                    pass
-                if len(temp_rules_dict[lhs]) == 0:
-                    del temp_rules_dict[lhs]
-            prune(temp_rules_dict)
-
-            if self.mark_verified:
-                for ver_label in temp_rules_dict.keys():
-                    self.set_verified(ver_label)
-
-            if self.equivdb[label] in temp_rules_dict:
-                for label_in_spec in temp_rules_dict.keys():
-                    self.set_verified(label_in_spec)
-                if iterative:
-                    specification = iterative_proof_tree_finder(
-                        temp_rules_dict, root=self.equivdb[label]
-                    )
-                else:
-                    specification = smallish_random_proof_tree(
-                        temp_rules_dict, self.equivdb[label], minimization_time_limit
-                    )
-                return specification
-
-        raise SpecificationNotFound("No specification for label {}".format(label))
diff --git a/comb_spec_searcher/specification.py b/comb_spec_searcher/specification.py
index ea6a0b6dd..d157a4680 100644
--- a/comb_spec_searcher/specification.py
+++ b/comb_spec_searcher/specification.py
@@ -11,6 +11,7 @@
 from logzero import logger
 from sympy import Eq, Expr, Function, Number, solve, var
 
+from comb_spec_searcher.exception import SpecificationNotFound
 from comb_spec_searcher.typing import (
     CombinatorialClassType,
     CombinatorialObjectType,
@@ -24,7 +25,7 @@
     InvalidOperationError,
     TaylorExpansionError,
 )
-from .isomorphism import AtomEquals, Bijection, Isomorphism
+from .isomorphism import Bijection, Isomorphism
 from .specification_drawer import SpecificationDrawer
 from .strategies import (
     EmptyStrategy,
@@ -58,12 +59,14 @@ def __init__(
         self,
         root: CombinatorialClassType,
         rules: Iterable[AbstractRule[CombinatorialClassType, CombinatorialObjectType]],
+        group_equiv: bool = True,
     ):
         self.root = root
         self.rules_dict = {rule.comb_class: rule for rule in rules}
         self.labels: Dict[CombinatorialClassType, int] = {}
         self._label_to_tiling: Dict[int, CombinatorialClassType] = {}
-        self._group_equiv_in_path()
+        if group_equiv:
+            self._group_equiv_in_path()
         self._set_subrules()
 
     def _set_subrules(self) -> None:
@@ -132,7 +135,9 @@ def expand_verified(self) -> None:
         """
         verification_packs: Dict[CombinatorialClassType, StrategyPack] = {}
         for comb_class, rule in self.rules_dict.items():
-            if isinstance(rule, VerificationRule):
+            if isinstance(rule, VerificationRule) and not isinstance(
+                rule.strategy, EmptyStrategy
+            ):
                 try:
                     verification_packs[comb_class] = rule.pack()
                 except InvalidOperationError:
@@ -144,8 +149,38 @@ def expand_verified(self) -> None:
             self._group_equiv_in_path()
             self._set_subrules()
 
+    def unexpanded_verified_classes(self) -> Set[CombinatorialClassType]:
+        """
+        Returns the set of verified classes so you can determine which classes still
+        need to be expanded.
+        """
+        return set(
+            comb_class
+            for comb_class, rule in self.rules_dict.items()
+            if isinstance(rule, VerificationRule)
+            and not isinstance(rule.strategy, EmptyStrategy)
+        )
+
+    def expand_all_verified_with_pack(
+        self, pack: StrategyPack, max_expansion_time: float = 0
+    ) -> None:
+        """
+        Will attempt to expand all verified class with the given strategy pack, limiting
+        expansion to a particular time limit if max_expansion_time is specified.
+        """
+        for comb_class, rule in list(self.rules_dict.items()):
+            if isinstance(rule, VerificationRule) and not isinstance(
+                rule.strategy, EmptyStrategy
+            ):
+                self.expand_comb_class(comb_class, pack, max_expansion_time)
+        self._group_equiv_in_path()
+        self._set_subrules()
+
     def expand_comb_class(
-        self, comb_class: Union[int, CombinatorialClassType], pack: StrategyPack
+        self,
+        comb_class: Union[int, CombinatorialClassType],
+        pack: StrategyPack,
+        max_expansion_time: Optional[float] = None,
     ) -> None:
         """
         Will try to expand a particular class with respect to the strategy pack
@@ -155,17 +190,25 @@ def expand_comb_class(
         from .comb_spec_searcher import CombinatorialSpecificationSearcher
 
         rule = self.get_rule(comb_class)
-        self.rules_dict.pop(rule.comb_class)
+
+        removed_rule = self.rules_dict.pop(rule.comb_class)
+
         pack = pack.add_verification(AlreadyVerified(self.rules_dict), apply_first=True)
         css = CombinatorialSpecificationSearcher(rule.comb_class, pack)
         logger.info(css.run_information())
-        rules = css._auto_search_rules()  # pylint: disable=protected-access
-        for rule in rules:
-            if not isinstance(rule.strategy, AlreadyVerified):
-                assert (
-                    rule.comb_class not in self.rules_dict or rule.comb_class.is_empty()
-                )
-                self.rules_dict[rule.comb_class] = rule
+        try:
+            rules = css._auto_search_rules(  # pylint: disable=protected-access
+                max_expansion_time=max_expansion_time
+            )
+            for rule in rules:
+                if not isinstance(rule.strategy, AlreadyVerified):
+                    assert (
+                        rule.comb_class not in self.rules_dict
+                        or rule.comb_class.is_empty()
+                    )
+                    self.rules_dict[rule.comb_class] = rule
+        except SpecificationNotFound:
+            self.rules_dict[rule.comb_class] = removed_rule
 
     def _is_valid_spec(self) -> bool:
         """Checks that each class is on a left hand side."""
@@ -190,7 +233,9 @@ def get_rule(
                     " in the specification."
                 ) from e
         if comb_class not in self.rules_dict:
-            assert comb_class.is_empty(), "rule not in the spec and not empty"
+            assert (
+                comb_class.is_empty()
+            ), f"rule not in the spec and not empty\n{comb_class}"
             empty_strat = EmptyStrategy[
                 CombinatorialClassType, CombinatorialObjectType
             ]()
@@ -230,12 +275,8 @@ def get_function(self, comb_class: CombinatorialClassType) -> Function:
         """
         Return a sympy function for the comb class, using the label it is
         assigned.
-
-        TODO: call comb_class for its parameters - 'x' is reserved for size.
         """
-        x = var("x")
-        extra_parameters = [var(k) for k in comb_class.extra_parameters]
-        return Function("F_{}".format(self.get_label(comb_class)))(x, *extra_parameters)
+        return comb_class.get_function(self.get_label)
 
     def get_equations(self) -> Iterator[Eq]:
         """
@@ -309,7 +350,7 @@ def get_genf(self, check: int = 6) -> Expr:
         logger.info("Solving...")
         solutions = solve(
             eqs,
-            tuple([eq.lhs for eq in eqs]),
+            tuple(eq.lhs for eq in eqs),
             dict=True,
             cubics=False,
             quartics=False,
@@ -398,26 +439,33 @@ def sanity_check(self, length: int = 5) -> bool:
 
         Raise an SanityCheckFailure error if it fails.
         """
-        return self._is_valid_spec() and all(
-            all(rule.sanity_check(n) for rule in self.rules_dict.values())
-            for n in range(length + 1)
-        )
+        if not self._is_valid_spec():
+            return False
+
+        for rule in self.rules_dict.values():
+            try:
+                for n in range(length + 1):
+                    if not rule.sanity_check(n):
+                        return False
+            except NotImplementedError:
+                logger.warning(
+                    "Can't sanity check the rule %s -> %s, which is\n" "%s",
+                    self.get_label(rule.comb_class),
+                    tuple(self.get_label(child) for child in rule.children),
+                    rule,
+                )
+
+        return True
 
     def get_bijection_to(
-        self, other: "CombinatorialSpecification", eq: Optional[AtomEquals] = None
+        self, other: "CombinatorialSpecification"
     ) -> Optional[Bijection]:
         """Get bijection from self to other."""
-        if eq is None:
-            return Bijection.construct(self, other)
-        return Bijection.construct(self, other, eq)
+        return Bijection.construct(self, other)
 
-    def are_isomorphic(
-        self, other: "CombinatorialSpecification", eq: Optional[AtomEquals] = None
-    ) -> bool:
+    def are_isomorphic(self, other: "CombinatorialSpecification") -> bool:
         """Check if self is isomorphic to other."""
-        if eq is None:
-            return Isomorphism.check(self, other)
-        return Isomorphism.check(self, other, eq)
+        return Isomorphism.check(self, other)
 
     def show(
         self, levels_shown: int = 0, levels_expand: int = 0, verbose: bool = False
@@ -461,7 +509,8 @@ def share(
             levels_expand=levels_expand,
             verbose=verbose,
         )
-        sd.share()
+        file_url = sd.share()
+        print(f"The specification is available at {file_url}.")
 
     def __eq__(self, other: object) -> bool:
         if not isinstance(other, CombinatorialSpecification):
@@ -536,7 +585,7 @@ def from_dict(cls, d):
         """
         root = CombinatorialClass.from_dict(d.pop("root"))
         rules = [AbstractRule.from_dict(rule_dict) for rule_dict in d.pop("rules")]
-        return CombinatorialSpecification(root, rules)
+        return CombinatorialSpecification(root, rules, group_equiv=False)
 
 
 class AlreadyVerified(VerificationStrategy[CombinatorialClass, CombinatorialObject]):
diff --git a/comb_spec_searcher/specification_drawer.py b/comb_spec_searcher/specification_drawer.py
index bb888c564..5c9fcdb55 100644
--- a/comb_spec_searcher/specification_drawer.py
+++ b/comb_spec_searcher/specification_drawer.py
@@ -9,14 +9,14 @@
 import uuid
 import webbrowser
 from copy import copy
-from typing import TYPE_CHECKING, ClassVar, Dict, Iterable, List, Tuple
+from typing import TYPE_CHECKING, ClassVar, Dict, Iterable, List, Tuple, cast
 
 import requests
 from logzero import logger
 from typing_extensions import TypedDict
 
 from .combinatorial_class import CombinatorialClass
-from .strategies import EquivalencePathRule, Rule, VerificationRule
+from .strategies import EquivalencePathRule, ReverseRule, Rule, VerificationRule
 
 if TYPE_CHECKING:
     from .specification import CombinatorialSpecification
@@ -78,7 +78,7 @@ def _get_new_node_id() -> str:
         return str(new_id)
 
     def _to_tree(
-        self, comb_class: CombinatorialClass, rec_depth: int = 0
+        self, comb_class: CombinatorialClass, flipped: bool = False, rec_depth: int = 0
     ) -> TreantNode:
         """Create the subtree rooted in comb_class"""
         try:
@@ -103,22 +103,44 @@ def _to_tree(
         else:
             comb_classes = [comb_class]
 
+        is_flip_rule = isinstance(rule, ReverseRule)
         if recursed:
             children = []
         else:
-            children = [self._to_tree(child, rec_depth + 1) for child in rule.children]
+            children = [
+                self._to_tree(
+                    child,
+                    is_flip_rule
+                    and child == cast(ReverseRule, rule).original_rule.comb_class,
+                    rec_depth + 1,
+                )
+                for child in rule.children
+            ]
 
         # if a node has children we create a delimiter node between them
         # in order to separate tooltip describing the steps taking in the tree
 
         if children and isinstance(rule, Rule):
+            delimiter_additional_style = ""
+            if is_flip_rule:
+                delimiter_additional_style += "background-color: #ff8080;"
+                delimiter_additional_style += "border-color: red;"
+                delimiter_additional_style += "border-width: 2px; "
             delimiter_node, node_id = self._create_delimiter_node(
-                rule, children, rec_depth
+                rule, children, rec_depth, delimiter_additional_style
             )
             self._create_delimiter_tooltip(rule, node_id)
             children = [delimiter_node]
 
-        html_node = self.comb_classes_to_html_node(comb_classes)
+        node_additional_style = ""
+        node_label_additional_style = ""
+        if flipped:
+            node_label_additional_style += "background-color: #ff8080;"
+            node_additional_style += "border-color: red;"
+            node_additional_style += "border-width: 3px; "
+        html_node = self.comb_classes_to_html_node(
+            comb_classes, node_additional_style, node_label_additional_style
+        )
         treant_node, node_id = self._create_standard_node(html_node, children)
         self._create_standard_tooltip(comb_classes, node_id)
         return treant_node
@@ -272,13 +294,14 @@ def show(self) -> None:
         fsd = ForestSpecificationDrawer([self])
         fsd.show()
 
-    def share(self) -> None:
+    def share(self) -> str:
         """
-        Upload the html of the specification on a file server and displays a link to the
+        Upload the html of the specification on a file server and return the link to the
+        file.
         file.
         """
         fsd = ForestSpecificationDrawer([self])
-        fsd.share()
+        return fsd.share()
 
     def to_html(self) -> str:
         """Returns a html string that contains the whole tree"""
@@ -450,9 +473,8 @@ def to_html(self) -> str:
     @staticmethod
     def export_html(html: str, file_name: str = "tree.html") -> None:
         """Creates a html file in current directory"""
-        text_file = open(file_name, "w", encoding="UTF8")
-        text_file.write(html)
-        text_file.close()
+        with open(file_name, "w", encoding="UTF8") as fp:
+            fp.write(html)
 
     def show(self) -> None:
         """Displays the forest in the web browser"""
@@ -460,9 +482,9 @@ def show(self) -> None:
         viewer = HTMLViewer()
         viewer.open_html(html_string)
 
-    def share(self) -> None:
+    def share(self) -> str:
         """
-        Upload the html of the forest on a file server and displays a link to the
+        Upload the html of the forest on a file server and return a link to the
         file.
         """
         html_string = self.to_html()
@@ -475,22 +497,20 @@ def share(self) -> None:
 
         # ask gofile.io which server it wants us to use
         logger.info("Sending specification to file host.")
-        req1 = requests.get("https://apiv2.gofile.io/getServer")
+        req1 = requests.get("https://api.gofile.io/getServer")
+        req1.raise_for_status()
         server = req1.json()["data"]["server"]
-        upload_url = f"https://{server}.gofile.io/upload"
+        upload_url = f"https://{server}.gofile.io/uploadFile"
 
         with open(file_path, "rb") as f:
             req2 = requests.post(upload_url, files={"filesUploaded": f})
+        req2.raise_for_status()
 
         os.remove(file_path)
         file_code = req2.json()["data"]["code"]
-        _, file_name = os.path.split(file_path)
-
-        file_enable_url = f"https://gofile.io/d/{file_code}"
-        file_url = f"https://{server}.gofile.io/download/{file_code}/{file_name}"
-
-        print(f"First click this: {file_enable_url}")
-        print(f"Then this link will work: {file_url}")
+        file_url = f"https://gofile.io/d/{file_code}"
+        logger.info("File uploaded to: %s", file_url)
+        return file_url
 
 
 class HTMLViewer:
diff --git a/comb_spec_searcher/specification_extrator.py b/comb_spec_searcher/specification_extrator.py
index bb28bdbed..b8dec2aad 100644
--- a/comb_spec_searcher/specification_extrator.py
+++ b/comb_spec_searcher/specification_extrator.py
@@ -1,24 +1,25 @@
 import itertools
-from typing import Dict, Iterator, List, Tuple
+from collections import deque
+from operator import itemgetter
+from typing import TYPE_CHECKING, Deque, Dict, Iterator, List, Set, Tuple
 
 from comb_spec_searcher.class_db import ClassDB
-from comb_spec_searcher.rule_db.base import RuleDBBase
 from comb_spec_searcher.strategies.rule import AbstractRule, Rule
 from comb_spec_searcher.tree_searcher import Node
 from comb_spec_searcher.typing import RuleKey
 
+if TYPE_CHECKING:
+    from comb_spec_searcher.rule_db.base import RuleDBBase
+
 
 class SpecificationRuleExtractor:
     def __init__(
-        self, root_label: int, root_node: Node, ruledb: RuleDBBase, classdb: ClassDB
+        self, root_label: int, root_node: Node, ruledb: "RuleDBBase", classdb: ClassDB
     ):
         self.ruledb = ruledb
         self.classdb = classdb
         self.root_label = root_label
-        self.eqv_rulekeys: List[RuleKey] = [
-            (node.label, tuple(sorted(child.label for child in node.children)))
-            for node in root_node.nodes()
-        ]
+        self.eqv_rulekeys = root_node.rule_keys()
         self.rules_dict: Dict[int, Tuple[int, ...]] = {}
         # A map from equiv label to an equivalent label actually in the tree.
         self.eqvparent_to_parent: Dict[int, int] = {}
@@ -32,10 +33,7 @@ def _populate_decompositions(self) -> None:
         """
         eqvrule_to_rule = self.ruledb.rule_from_equivalence_rule_dict(self.eqv_rulekeys)
         for eqvrule in self.eqv_rulekeys:
-            try:
-                parent, children = eqvrule_to_rule[eqvrule]
-            except KeyError:
-                continue
+            parent, children = eqvrule_to_rule[eqvrule]
             self.rules_dict[parent] = children
             self.eqvparent_to_parent[eqvrule[0]] = parent
 
@@ -44,23 +42,30 @@ def _populate_equivalences(self) -> None:
         Populate the rules dict with the labels of equivalences rules needed to make
         every class be on a right hand side.
         """
-        all_rhs = itertools.chain.from_iterable(self.rules_dict.values())
-        no_lhs_labels = set(
-            itertools.filterfalse(self.rules_dict.__contains__, all_rhs)
-        )
-        if self.root_label not in self.rules_dict:
-            no_lhs_labels.add(self.root_label)
-        for label in no_lhs_labels:
+        for label in self._no_lhs_labels():
             eqv_label = self.ruledb.equivdb[label]
             path = self.ruledb.equivdb.find_path(
                 label, self.eqvparent_to_parent[eqv_label]
             )
             for parent, child in zip(path[:-1], path[1:]):
-                assert parent not in self.rules_dict or self.rules_dict[parent] == (
-                    child,
-                ), (parent, child, self.rules_dict[parent])
+                if parent in self.rules_dict:
+                    # No need to keep adding the path since the class already connects
+                    # with something in the specification
+                    break
                 self.rules_dict[parent] = (child,)
 
+    def _no_lhs_labels(self) -> Set[int]:
+        """
+        Find labels that are not a lhs of any rule.
+        """
+        all_rhs = itertools.chain.from_iterable(self.rules_dict.values())
+        no_lhs_labels = set(
+            itertools.filterfalse(self.rules_dict.__contains__, all_rhs)
+        )
+        if self.root_label not in self.rules_dict:
+            no_lhs_labels.add(self.root_label)
+        return no_lhs_labels
+
     def _check(self):
         """
         Perform a check to make sure that the spec make sense.
@@ -71,6 +76,8 @@ def _check(self):
             )
         )
         assert all_rhs.issubset(self.rules_dict), "Something is not on the lhs"
+        extra_lhs = {c for c in self.rules_dict if c not in all_rhs}
+        assert extra_lhs.issubset({self.root_label})
 
     def _find_rule(self, parent: int, children: Tuple[int, ...]) -> AbstractRule:
         try:
@@ -115,3 +122,189 @@ def _missing_rule_error(self, parent: int, children: Tuple[int, ...]) -> ValueEr
     def rules(self) -> Iterator[AbstractRule]:
         for parent, children in self.rules_dict.items():
             yield self._find_rule(parent, children)
+
+
+class PartialSpecificationRuleExtractor(SpecificationRuleExtractor):
+    """A base class for extractors for incomplete specifications."""
+
+    def _populate_decompositions(self) -> None:
+        pass
+
+    def _populate_equivalences(self) -> None:
+        # SpecificationRuleExtractor's with a try-catch do deal with incomplete spec
+        for label in self._no_lhs_labels():
+            eqv_label = self.ruledb.equivdb[label]
+            # Might not exist because spec is incomplete.
+            try:
+                path = self.ruledb.equivdb.find_path(
+                    label, self.eqvparent_to_parent[eqv_label]
+                )
+            except KeyError:
+                continue
+            for parent, child in zip(path[:-1], path[1:]):
+                if parent in self.rules_dict:
+                    break
+                self.rules_dict[parent] = (child,)
+
+    def _check(self):
+        pass
+
+    def _ordered_eqvrule_to_rule(
+        self,
+    ) -> Tuple[Dict[RuleKey, RuleKey], Dict[RuleKey, Tuple[int, ...]]]:
+        """Altered version of rule_from_equivalence_rule_dict that deals with order."""
+        eqv_rules = set(self.eqv_rulekeys)
+        eqvrule_to_rule: Dict[RuleKey, RuleKey] = {}
+        # Maps our children index order to that of the actual rules
+        index_order_map: Dict[RuleKey, Tuple[int, ...]] = {}
+        for start, ends in self.ruledb.rule_to_strategy:
+            eqv_start = self.ruledb.equivdb[start]
+            # Sort but keep knowledge of original order
+            index_order, eqv_ends = (
+                zip(
+                    *sorted(
+                        enumerate(map(self.ruledb.equivdb.__getitem__, ends)),
+                        key=itemgetter(1),
+                    )
+                )
+                if ends
+                else ((), ())
+            )
+
+            eqv_key = (eqv_start, eqv_ends)
+            index_order_map[eqv_key] = index_order
+            if eqv_key in eqv_rules:
+                eqvrule_to_rule[eqv_key] = (start, ends)
+        return eqvrule_to_rule, index_order_map
+
+
+class EquivalenceRuleExtractor(PartialSpecificationRuleExtractor):
+    """
+    This class extracts the actual rules that would be created for an equivalence
+    label given a partial spec that is traversable from root. Suppose we want the
+    rules for the parent in `grandparent-parent-children` than what rules are used
+    will depend on the grandparent and children and we are looking for the rules
+    that have the same equivalance label on both the LHS and RHS and in the order they
+    would be traversed from the root.
+    """
+
+    def __init__(
+        self,
+        root_eq_label: int,
+        root_class_label: int,
+        root_node: Node,
+        ruledb: "RuleDBBase",
+        classdb: ClassDB,
+        target: int,
+        parent_of_target: int,
+        idx: int,
+    ):
+        # parent
+        self.target = target
+        # grandparent
+        self.parent_of_target = parent_of_target
+        # grandparent.children[idx] is parent
+        self.idx = idx
+        # class label of the LHS of the first rule for the eq path, updated later
+        self.start = root_class_label
+        # class label of the RHS of the last rule for the eq path, updated later
+        self.end = -1
+        # Children are included in the tree
+        super().__init__(root_eq_label, root_node, ruledb, classdb)
+
+    def _populate_decompositions(self) -> None:
+        eqvrule_to_rule, index_order_map = self._ordered_eqvrule_to_rule()
+        for eqvrule in self.eqv_rulekeys:
+            # Might not exist because spec is incomplete.
+            try:
+                parent, children = eqvrule_to_rule[eqvrule]
+            except KeyError:
+                continue
+            self.rules_dict[parent] = children
+            self.eqvparent_to_parent[eqvrule[0]] = parent
+            # If grandparent, then we can find our starting point
+            # If no grandparent exists, this will always be the root.
+            if eqvrule[0] == self.parent_of_target:
+                self.start = children[index_order_map[eqvrule][self.idx]]
+            # If parent, we have our end point
+            if eqvrule[0] == self.target:
+                self.end = parent
+
+    def nonequivalent_rules_in_equiv_path(self) -> List[Rule]:
+        return list(self._nonequivalent_rules_in_equiv_path())
+
+    def _nonequivalent_rules_in_equiv_path(self) -> Iterator[Rule]:
+        path = self.ruledb.equivdb.find_path(self.start, self.end)
+        for p, c in zip(path, path[1:]):
+            rule = self._find_rule(p, (c,))
+            # If rule is an actual equivalence rule, we don't care and ignore it.
+            if not rule.is_equivalence():
+                assert isinstance(rule, Rule)
+                yield rule
+
+
+class RulePathToAtomExtractor(PartialSpecificationRuleExtractor):
+    """
+    This class extracts rules from the root down to an atom in a specification (in terms
+    of eq label maps) that is not complete but contains a path from root to said atom.
+    """
+
+    def __init__(
+        self,
+        root_class_label: int,
+        root_label: int,
+        root_node: Node,
+        ruledb: "RuleDBBase",
+        classdb: ClassDB,
+        path: List[Tuple[int, int]],
+    ):
+        self.root_class_label = root_class_label
+        self.path = path
+        self.path.pop()  # remove padded -1
+        self.search_order: Dict[int, Tuple[int, ...]] = {}
+        super().__init__(root_label, root_node, ruledb, classdb)
+
+    def _populate_decompositions(self) -> None:
+        eqvrule_to_rule, index_order_map = self._ordered_eqvrule_to_rule()
+        for eqvrule in self.eqv_rulekeys:
+            try:
+                parent, children = eqvrule_to_rule[eqvrule]
+            except KeyError:
+                continue
+            self.rules_dict[parent] = children
+            self.eqvparent_to_parent[eqvrule[0]] = parent
+
+            # Order in searcher
+            order = index_order_map[eqvrule]
+            self.search_order[parent] = tuple(
+                children[order[i]] for i in range(len(children))
+            )
+
+    def rule_path(self) -> Deque[Tuple[AbstractRule, int]]:
+        curr_class = self.root_class_label
+        path: Deque[Tuple[AbstractRule, int]] = deque([])
+        while self.path:
+            children = self.rules_dict[curr_class]
+            if (
+                len(children) == 1
+                and self.ruledb.equivdb[curr_class] == self.ruledb.equivdb[children[0]]
+            ):
+                # If rule is within eq class, we add it with idx 0
+                path.append((self._find_rule(curr_class, children), 0))
+                curr_class = children[0]
+            else:
+                _, idx = self.path.pop()
+                rule = self._find_rule(curr_class, children)
+                # The idx given by the searcher might differ from that of rule_dict
+                # and we use search_order to get the corresponding index in rule_dict
+                next_child = self.search_order[curr_class][idx]
+                # Get the first child that matches the label of the child we want.
+                # If there are multiple, it does not matter which one.
+                rule_idx = next(
+                    i
+                    for i, c in enumerate(rule.children)
+                    if self.classdb.get_label(c) == next_child
+                )
+                path.append((rule, rule_idx))
+                curr_class = next_child
+        return path
diff --git a/comb_spec_searcher/strategies/__init__.py b/comb_spec_searcher/strategies/__init__.py
index 2bba8ad17..8a67a3ff2 100644
--- a/comb_spec_searcher/strategies/__init__.py
+++ b/comb_spec_searcher/strategies/__init__.py
@@ -5,7 +5,13 @@
     DisjointUnion,
     Quotient,
 )
-from .rule import EquivalencePathRule, ReverseRule, Rule, VerificationRule
+from .rule import (
+    EquivalencePathRule,
+    NonBijectiveRule,
+    ReverseRule,
+    Rule,
+    VerificationRule,
+)
 from .strategy import (
     AbstractStrategy,
     AtomStrategy,
@@ -28,6 +34,7 @@
     "EquivalencePathRule",
     "ReverseRule",
     "Rule",
+    "NonBijectiveRule",
     "VerificationRule",
     "AbstractStrategy",
     "AtomStrategy",
diff --git a/comb_spec_searcher/strategies/constructor/base.py b/comb_spec_searcher/strategies/constructor/base.py
index 2469ce6a0..7245fb35a 100644
--- a/comb_spec_searcher/strategies/constructor/base.py
+++ b/comb_spec_searcher/strategies/constructor/base.py
@@ -3,7 +3,9 @@
 enumeration, generate objects, and sample objects.
 """
 import abc
-from typing import Callable, Generic, Iterator, List, Optional, Tuple
+from collections import Counter
+from functools import partial
+from typing import Callable, Dict, Generic, Iterator, List, Optional, Set, Tuple
 
 import sympy
 
@@ -70,6 +72,19 @@ def random_sample_sub_objects(
         """Return a randomly sampled subobjs/image of the bijection implied
         by the constructor."""
 
+    @staticmethod
+    def param_map(
+        child_pos_to_parent_pos: Tuple[Tuple[int, ...], ...],
+        num_parent_params: int,
+        param: Parameters,
+    ) -> Parameters:
+        new_params = [0 for _ in range(num_parent_params)]
+        for pos, value in enumerate(param):
+            parent_pos = child_pos_to_parent_pos[pos]
+            for p in parent_pos:
+                new_params[p] += value
+        return tuple(new_params)
+
     @staticmethod
     def build_param_map(
         child_pos_to_parent_pos: Tuple[Tuple[int, ...], ...], num_parent_params: int
@@ -77,13 +92,83 @@ def build_param_map(
         """
         Return a parameters map according to the given pos map.
         """
+        return partial(
+            Constructor.param_map,
+            child_pos_to_parent_pos,
+            num_parent_params,
+        )
+
+    @abc.abstractmethod
+    def equiv(
+        self, other: "Constructor", data: Optional[object] = None
+    ) -> Tuple[bool, Optional[object]]:
+        """Two constructors are equiv if A <self> B and C <other> D are isomorphic for
+        any combinatorial classes A, B, C and D where A and C are isomorphic and B and
+        D. This must be implemented for bijections. The data argument is to optionally
+        pass additional information to the equiv function. The second returned value is
+        optional data for bijections. It is used to pass additional arguments to
+        determine index ordering. It needs to be JSON compatible for bijection's
+        `to_jsonable` and `from_dict` to work.
+        """
+
+    @staticmethod
+    def extra_params_equiv(
+        params1: Tuple[Dict[str, str], ...],
+        params2: Tuple[Dict[str, str], ...],
+    ) -> bool:
+        """A base equiv for extra params that will suffice for most constructors."""
+        non_empty1 = tuple(par for par in params1 if par)
+        non_empty2 = tuple(par for par in params2 if par)
+        n = len(non_empty1)
+        if len(non_empty2) != n:
+            return False
+        if n == 0:
+            return True
+        return Constructor._extra_params_match_bijection(n, non_empty1, non_empty2)
+
+    @staticmethod
+    def _extra_params_match_single(par1: Dict[str, str], par2: Dict[str, str]) -> bool:
+        """Two parameter dictionaries are equivalent if they have the same number
+        of keys if there is a value with k keys mapping to it, then there is a
+        corresponding value with k keys mapping to it in the other one.
+        """
+        return len(par1) == len(par2) and sorted(
+            Counter(par1.values()).values()
+        ) == sorted(Counter(par2.values()).values())
+
+    @staticmethod
+    def _extra_params_match_bijection(
+        n: int, p1: Tuple[Dict[str, str], ...], p2: Tuple[Dict[str, str], ...]
+    ) -> bool:
+        """Check if there is a bijection between the two tuples of params
+        such that they are pairwise equivalent."""
+
+        # Shared variables in recursion
+        bijection: List[int] = []
+        in_use: Set[int] = set()
+        cache: Dict[Tuple[int, int], bool] = {}
+
+        def _backtracking() -> bool:
+            for i in range(n):
+                if i in in_use:
+                    continue
+                bijection.append(i)
+                in_use.add(i)
+
+                k = len(bijection)
+                key = (k - 1, i)
+                res = cache.get(key, None)
+                if res is None:
+                    res = Constructor._extra_params_match_single(p1[key[0]], p2[key[1]])
+                    cache[key] = res
+
+                # Check if consistent and either done or recursive success.
+                res = res and (k == n or _backtracking())
 
-        def param_map(param: Parameters) -> Parameters:
-            new_params = [0 for _ in range(num_parent_params)]
-            for pos, value in enumerate(param):
-                parent_pos = child_pos_to_parent_pos[pos]
-                for p in parent_pos:
-                    new_params[p] += value
-            return tuple(new_params)
+                bijection.pop()
+                in_use.remove(i)
+                if res:
+                    return True
+            return False
 
-        return param_map
+        return _backtracking()
diff --git a/comb_spec_searcher/strategies/constructor/cartesian.py b/comb_spec_searcher/strategies/constructor/cartesian.py
index 105634964..8f85d3ba0 100644
--- a/comb_spec_searcher/strategies/constructor/cartesian.py
+++ b/comb_spec_searcher/strategies/constructor/cartesian.py
@@ -1,5 +1,6 @@
 import random
 from collections import defaultdict
+from functools import partial
 from itertools import product
 from typing import (
     Callable,
@@ -350,6 +351,17 @@ def get_op_symbol() -> str:
     def __str__(self) -> str:
         return "Cartesian product"
 
+    def equiv(
+        self, other: "Constructor", data: Optional[object] = None
+    ) -> Tuple[bool, Optional[object]]:
+        return (
+            isinstance(other, type(self))
+            and CartesianProduct.extra_params_equiv(
+                self.extra_parameters, other.extra_parameters
+            ),
+            None,
+        )
+
 
 class Quotient(Constructor[CombinatorialClassType, CombinatorialObjectType]):
     """
@@ -393,21 +405,27 @@ def __init__(
         self._parent_shift = sum(self._min_sizes) - self._min_sizes[self.idx]
         self._num_parent_params = len(parent.extra_parameters)
 
+    @staticmethod
+    def param_map(
+        child_pos_to_parent_pos: Tuple[Tuple[int, ...], ...],
+        num_parent_params: int,
+        param: Parameters,
+    ) -> Parameters:
+        new_params: List[Optional[int]] = [None for _ in range(num_parent_params)]
+        for pos, value in enumerate(param):
+            parent_pos = child_pos_to_parent_pos[pos]
+            for p in parent_pos:
+                assert new_params[p] is None or new_params[p] == value
+                new_params[p] = value
+        assert all(p is not None for p in new_params)
+        return cast(Parameters, tuple(new_params))
+
     @staticmethod
     def build_param_map(
         child_pos_to_parent_pos: Tuple[Tuple[int, ...], ...], num_parent_params: int
     ) -> ParametersMap:
-        def param_map(param: Parameters) -> Parameters:
-            new_params: List[Optional[int]] = [None for _ in range(num_parent_params)]
-            for pos, value in enumerate(param):
-                parent_pos = child_pos_to_parent_pos[pos]
-                for p in parent_pos:
-                    assert new_params[p] is None or new_params[p] == value
-                    new_params[p] = value
-            assert all(p is not None for p in new_params)
-            return cast(Parameters, tuple(new_params))
-
-        return param_map
+
+        return partial(Quotient.param_map, child_pos_to_parent_pos, num_parent_params)
 
     def _build_parent_param_map(
         self,
@@ -523,10 +541,10 @@ def _terms_to_poly(self, terms: Terms) -> Union[sympy.Poly, int]:
 
     @staticmethod
     def _poly_to_terms(poly: Union[sympy.Poly, int]) -> Terms:
-        if isinstance(poly, int):
+        if isinstance(poly, (int, sympy.core.numbers.Integer)):
             if poly == 0:
                 return Counter()
-            return Counter({tuple(): poly})
+            return Counter({tuple(): int(poly)})
         return Counter(poly.as_dict())
 
     def _b(
@@ -566,16 +584,15 @@ def get_terms(
     def get_equation(
         self, lhs_func: sympy.Function, rhs_funcs: Tuple[sympy.Function, ...]
     ) -> sympy.Eq:
-        res = lhs_func.subs(self.extra_parameters[self.idx])
-        for (idx, rhs_func), extra_parameters in zip(
-            enumerate(rhs_funcs), self.extra_parameters
-        ):
-            if self.idx != idx:
-                res /= rhs_func.subs(
-                    {child: parent for parent, child in extra_parameters.items()},
-                    simultaneous=True,
-                ).subs(self.extra_parameters[self.idx], simultaneous=True)
-        return sympy.Eq(rhs_funcs[self.idx], res)
+        if any(self.extra_parameters):
+            raise NotImplementedError(
+                "Quotient equation is not implemented with extra parameters. "
+                "You can fall back on the cartesian equation."
+            )
+        res = rhs_funcs[0]
+        for rhs_func in rhs_funcs[1:]:
+            res /= rhs_func
+        return sympy.Eq(lhs_func, res)
 
     def reliance_profile(self, n: int, **parameters: int) -> RelianceProfile:
         raise NotImplementedError
@@ -596,3 +613,17 @@ def random_sample_sub_objects(
         **parameters: int,
     ) -> Tuple[Optional[CombinatorialObjectType], ...]:
         raise NotImplementedError
+
+    def __str__(self):
+        return "quotient"
+
+    def equiv(
+        self, other: "Constructor", data: Optional[object] = None
+    ) -> Tuple[bool, Optional[object]]:
+        return (
+            isinstance(other, type(self))
+            and Quotient.extra_params_equiv(
+                self.extra_parameters, other.extra_parameters
+            ),
+            None,
+        )
diff --git a/comb_spec_searcher/strategies/constructor/disjoint.py b/comb_spec_searcher/strategies/constructor/disjoint.py
index 321d793fb..edaf41d1b 100644
--- a/comb_spec_searcher/strategies/constructor/disjoint.py
+++ b/comb_spec_searcher/strategies/constructor/disjoint.py
@@ -67,8 +67,14 @@ def get_equation(
     ) -> sympy.Eq:
         res = 0
         for rhs_func, extra_parameters in zip(rhs_funcs, self.extra_parameters):
+            subs: Dict[str, sympy.Expr] = dict()
+            for parent, child in extra_parameters.items():
+                if child in subs:
+                    subs[child] *= sympy.var(parent)
+                else:
+                    subs[child] = sympy.var(parent)
             res += rhs_func.subs(
-                {child: parent for parent, child in extra_parameters.items()},
+                subs,
                 simultaneous=True,
             )
         return sympy.Eq(lhs_func, res)
@@ -194,6 +200,17 @@ def get_op_symbol() -> str:
     def __str__(self):
         return "disjoint union"
 
+    def equiv(
+        self, other: "Constructor", data: Optional[object] = None
+    ) -> Tuple[bool, Optional[object]]:
+        return (
+            isinstance(other, type(self))
+            and DisjointUnion.extra_params_equiv(
+                self.extra_parameters, other.extra_parameters
+            ),
+            None,
+        )
+
 
 class Complement(Constructor[CombinatorialClassType, CombinatorialObjectType]):
     """
@@ -238,7 +255,11 @@ def _build_children_param_maps(
         parent_param_to_pos = {
             param: pos for pos, param in enumerate(parent.extra_parameters)
         }
-        for child, extra_param in zip(children, self.extra_parameters):
+        for (idx, child), extra_param in zip(
+            enumerate(children), self.extra_parameters
+        ):
+            if idx == self.idx:
+                continue
             reversed_extra_param: Dict[str, List[str]] = defaultdict(list)
             for parent_var, child_var in extra_param.items():
                 reversed_extra_param[child_var].append(parent_var)
@@ -282,16 +303,15 @@ def _build_parent_param_map(
     def get_equation(
         self, lhs_func: sympy.Function, rhs_funcs: Tuple[sympy.Function, ...]
     ) -> sympy.Eq:
-        res = lhs_func.subs(self.extra_parameters[self.idx])
-        for (idx, rhs_func), extra_parameters in zip(
-            enumerate(rhs_funcs), self.extra_parameters
-        ):
-            if self.idx != idx:
-                res -= rhs_func.subs(
-                    {child: parent for parent, child in extra_parameters.items()},
-                    simultaneous=True,
-                ).subs(self.extra_parameters[self.idx], simultaneous=True)
-        return sympy.Eq(rhs_funcs[self.idx], res)
+        if any(self.extra_parameters):
+            raise NotImplementedError(
+                "Complement equation is not implemented with extra parameters. "
+                "You can fall back on the union equations."
+            )
+        res = rhs_funcs[0]
+        for rhs_func in rhs_funcs[1:]:
+            res -= rhs_func
+        return sympy.Eq(lhs_func, res)
 
     def reliance_profile(self, n: int, **parameters: int) -> RelianceProfile:
         raise NotImplementedError
@@ -303,19 +323,15 @@ def get_terms(
         for param, value in subterms[0](n).items():
             if value:
                 parent_terms_mapped[self._parent_param_map(param)] += value
-
         children_terms = subterms[1:]
-        for (idx, child_terms), param_map in zip(
-            enumerate(children_terms), self._children_param_maps
-        ):
-            if self.idx != idx:
-                # we subtract from total
-                for param, value in child_terms(n).items():
-                    mapped_param = self._parent_param_map(param_map(param))
-                    parent_terms_mapped[mapped_param] -= value
-                    assert parent_terms_mapped[mapped_param] >= 0
-                    if parent_terms_mapped[mapped_param] == 0:
-                        parent_terms_mapped.pop(mapped_param)
+        for child_terms, param_map in zip(children_terms, self._children_param_maps):
+            # we subtract from total
+            for param, value in child_terms(n).items():
+                mapped_param = self._parent_param_map(param_map(param))
+                parent_terms_mapped[mapped_param] -= value
+                assert parent_terms_mapped[mapped_param] >= 0
+                if parent_terms_mapped[mapped_param] == 0:
+                    parent_terms_mapped.pop(mapped_param)
 
         return parent_terms_mapped
 
@@ -335,3 +351,17 @@ def random_sample_sub_objects(
         **parameters: int,
     ) -> Tuple[Optional[CombinatorialObjectType], ...]:
         raise NotImplementedError
+
+    def __str__(self):
+        return "complement"
+
+    def equiv(
+        self, other: "Constructor", data: Optional[object] = None
+    ) -> Tuple[bool, Optional[object]]:
+        return (
+            isinstance(other, type(self))
+            and Complement.extra_params_equiv(
+                self.extra_parameters, other.extra_parameters
+            ),
+            None,
+        )
diff --git a/comb_spec_searcher/strategies/rule.py b/comb_spec_searcher/strategies/rule.py
index 70b52c9c1..6042f565b 100644
--- a/comb_spec_searcher/strategies/rule.py
+++ b/comb_spec_searcher/strategies/rule.py
@@ -1,10 +1,11 @@
+# pylint: disable=too-many-lines
 """
 The rule class is used for a specific application of a strategy on a combinatorial
 class. This is not something the user should implement, as it is just a wrapper for
 calling the Strategy class and storing its results.
-
-A CombinatorialSpecification is (more or less) a set of Rule.
 """
+
+
 import abc
 import random
 from collections import defaultdict
@@ -29,8 +30,10 @@
 
 from comb_spec_searcher.combinatorial_class import CombinatorialClass
 from comb_spec_searcher.typing import (
+    ForestRuleKey,
     Objects,
     ObjectsCache,
+    RuleBucket,
     SubObjects,
     SubTerms,
     Terms,
@@ -50,7 +53,8 @@
 
 class AbstractRule(abc.ABC, Generic[CombinatorialClassType, CombinatorialObjectType]):
     """
-    An instance of Rule is created by the __call__ method of strategy.
+    An application of a strategy to a comb_class. If the children are not provided they
+    are computed from the strategy.
     """
 
     def __init__(
@@ -74,6 +78,7 @@ def __init__(
         self.subobjects: Optional[SubObjects] = None
         self._children = children
         self._non_empty_children: Optional[Tuple[CombinatorialClassType, ...]] = None
+        self._shifts: Optional[Tuple[int, ...]] = None
 
     def to_jsonable(self) -> dict:
         d = {
@@ -180,6 +185,24 @@ def is_two_way(self) -> bool:
         """
         return self.strategy.is_two_way(self.comb_class)
 
+    def is_reversible(self) -> bool:
+        """
+        Returns True if it is possible to compute the enumeration of a child using the
+        enumeration of the parent and the other children.
+        """
+        return self.strategy.is_reversible(self.comb_class)
+
+    def shifts(self) -> Tuple[int, ...]:
+        if self._shifts is None:
+            self._shifts = self.strategy.shifts(self.comb_class, self.children)
+        return self._shifts
+
+    @abc.abstractmethod
+    def forest_key(
+        self, get_label: Callable[[CombinatorialClassType], int]
+    ) -> ForestRuleKey:
+        """Return the forest key of the rule."""
+
     def non_empty_children(self) -> Tuple[CombinatorialClassType, ...]:
         """
         Return the tuple of non-empty combinatorial classes that are found
@@ -258,7 +281,6 @@ def sanity_check(self, n: int) -> bool:
 
         Raise a SanityCheckFailure error if the sanity_check fails.
         """
-        raise NotImplementedError
 
     def __eq__(self, other: object) -> bool:
         if not isinstance(other, AbstractRule):
@@ -295,7 +317,7 @@ def join(res, new):
 
         res = str(self.comb_class).split("\n")
         backpad(res)
-        if isinstance(self, Rule):
+        if self.children:
             children = [str(child).split("\n") for child in self.children]
             symbol_height = min(1, len(res) - 1)
             eq_symbol = (
@@ -371,9 +393,16 @@ def constructor(
         return self._constructor
 
     def is_equivalence(self):
-        return (
-            isinstance(self.constructor, (DisjointUnion, Complement))
-            and len(self.non_empty_children()) == 1
+        return self.strategy.can_be_equivalent() and len(self.non_empty_children()) == 1
+
+    def forest_key(
+        self, get_label: Callable[[CombinatorialClassType], int]
+    ) -> ForestRuleKey:
+        return ForestRuleKey(
+            get_label(self.comb_class),
+            tuple(map(get_label, self.children)),
+            self.shifts(),
+            RuleBucket.EQUIV if self.is_equivalence() else RuleBucket.NORMAL,
         )
 
     def backward_map(
@@ -394,10 +423,32 @@ def forward_map(
         """
         return self.strategy.forward_map(self.comb_class, obj, self.children)
 
+    def indexed_forward_map(
+        self,
+        obj: CombinatorialObjectType,
+        data: Optional[object] = None,
+    ) -> Tuple[Tuple[Optional[CombinatorialObjectType], ...], int]:
+        """
+        A version of forward_map that is used for bijections. The idx is used
+        create a bijection from a nonbijective rule.
+        """
+        return self.forward_map(obj), 0
+
+    def indexed_backward_map(
+        self,
+        objs: Tuple[Optional[CombinatorialObjectType], ...],
+        idx: int,
+        data: Optional[object] = None,
+    ) -> CombinatorialObjectType:
+        """
+        A version of backward_map that is used for bijections. The idx is used
+        create a bijection from a nonbijective rule.
+        """
+        return next(self.backward_map(objs))
+
     def to_equivalence_rule(self) -> "EquivalenceRule":
         """
-        Return the reverse rule. At this stage, reverse rules can only be
-        created for equivalence rules.
+        Return the version of the rule where the empty children are removed.
         """
         assert (
             self.is_equivalence()
@@ -406,9 +457,9 @@ def to_equivalence_rule(self) -> "EquivalenceRule":
 
     def to_reverse_rule(self, idx: int) -> "Rule":
         """
-        Return the reverse rule. At this stage, reverse rules can only be
-        created for equivalence rules.
+        Return the reverse rule where the child at the given index is the parent.
         """
+        assert self.is_reversible()
         return ReverseRule(self, idx)
 
     def _ensure_level(self, n: int) -> None:
@@ -585,8 +636,51 @@ def subrec(n: int, **parameters: int) -> int:
         return True
 
 
+class NonBijectiveRule(Rule[CombinatorialClassType, CombinatorialObjectType]):
+    def indexed_forward_map(
+        self,
+        obj: CombinatorialObjectType,
+        data: Optional[object] = None,
+    ) -> Tuple[Tuple[Optional[CombinatorialObjectType], ...], int]:
+        img = self.forward_map(obj)
+        return img, self._forward_order(obj, img, data)
+
+    def indexed_backward_map(
+        self,
+        objs: Tuple[Optional[CombinatorialObjectType], ...],
+        idx: int,
+        data: Optional[object] = None,
+    ) -> CombinatorialObjectType:
+        return self._backward_order_item(idx, objs, data)
+
+    @abc.abstractmethod
+    def _forward_order(
+        self,
+        obj: CombinatorialObjectType,
+        image: Tuple[Optional[CombinatorialObjectType], ...],
+        data: Optional[object] = None,
+    ) -> int:
+        """For rules that do not have an injective forward_map, we can mark the
+        resulting object with a number i (that depends on the object we map from)
+        in order to achieve a bijection. The backward map will then take said
+        number and yield the i-th element (by some predefined order). This requires
+        the backward_map's order to be deterministic.
+        """
+
+    @abc.abstractmethod
+    def _backward_order_item(
+        self,
+        idx: int,
+        objs: Tuple[Optional[CombinatorialObjectType], ...],
+        data: Optional[object] = None,
+    ) -> CombinatorialObjectType:
+        """Given an iterator of all possible domain elements for an indexed forward map
+        and an index, this method will return the element that maps to said index.
+        """
+
+
 class EquivalenceRule(Rule[CombinatorialClassType, CombinatorialObjectType]):
-    def __init__(self, rule: Rule):
+    def __init__(self, rule: Rule[CombinatorialClassType, CombinatorialObjectType]):
         non_empty_children = rule.non_empty_children()
         assert rule.is_equivalence(), "not an equivalence rule: {}".format(str(rule))
         child = non_empty_children[0]
@@ -613,10 +707,6 @@ def from_dict(cls, d: dict) -> "EquivalenceRule":
 
     @property
     def constructor(self) -> Union[DisjointUnion, Complement]:
-        """
-        Return the constructor, that contains all the information about how to
-        count/generate objects from the rule.
-        """
         if self._constructor is None:
             original_constructor = self.original_rule.constructor
             if isinstance(original_constructor, DisjointUnion):
@@ -627,15 +717,21 @@ def constructor(self) -> Union[DisjointUnion, Complement]:
                 )
             elif isinstance(original_constructor, Complement):
                 assert isinstance(self.original_rule, ReverseRule)
-                original_original_constructor = (
-                    self.original_rule.original_rule.constructor
+                original_original_rule = self.original_rule.original_rule
+                original_original_child_idx = (
+                    original_original_rule.to_equivalence_rule().child_idx
                 )
+                original_original_constructor = original_original_rule.constructor
                 assert isinstance(original_original_constructor, DisjointUnion)
                 self._constructor = Complement(
-                    self.original_rule.original_rule.comb_class,
-                    self.original_rule.original_rule.children,
-                    self.original_rule.idx,
-                    original_original_constructor.extra_parameters,
+                    self.children[0],
+                    (self.comb_class,),
+                    0,
+                    (
+                        original_original_constructor.extra_parameters[
+                            original_original_child_idx
+                        ],
+                    ),
                 )
             else:
                 raise NotImplementedError
@@ -649,11 +745,20 @@ def to_reverse_rule(self, idx: int) -> "EquivalenceRule":
         assert idx == 0
         return self.original_rule.to_reverse_rule(self.child_idx).to_equivalence_rule()
 
+    def to_equivalence_rule(self) -> "EquivalenceRule":
+        raise NotImplementedError("You don't want to do that! I promise")
+
     @property
     def formal_step(self) -> str:
-        return "{} but only the child at index {} is non-empty".format(
-            self.strategy.formal_step(), self.child_idx
-        )
+        strat = self.strategy.formal_step()
+        if isinstance(self.original_rule, ReverseRule):
+            idx = self.original_rule.idx
+        else:
+            idx = self.child_idx
+        s = f"{strat} but only child and index {idx} is non-empty"
+        if isinstance(self.original_rule, ReverseRule):
+            return f"reverse of '{s}'"
+        return s
 
     def backward_map(
         self, objs: Tuple[Optional[CombinatorialObjectType], ...]
@@ -662,18 +767,12 @@ def backward_map(
             objs[0] if i == self.child_idx else None
             for i in range(len(self.actual_children))
         )
-        yield from self.strategy.backward_map(
-            self.comb_class, actual_objs, self.actual_children
-        )
+        yield from self.original_rule.backward_map(actual_objs)
 
     def forward_map(
         self, obj: CombinatorialObjectType
     ) -> Tuple[Optional[CombinatorialObjectType], ...]:
-        return (
-            self.strategy.forward_map(self.comb_class, obj, self.actual_children)[
-                self.child_idx
-            ],
-        )
+        return (self.original_rule.forward_map(obj)[self.child_idx],)
 
 
 class EquivalencePathRule(Rule[CombinatorialClassType, CombinatorialObjectType]):
@@ -682,8 +781,11 @@ class EquivalencePathRule(Rule[CombinatorialClassType, CombinatorialObjectType])
     single Rule.
     """
 
-    def __init__(self, rules: Sequence[Rule]):
+    def __init__(
+        self, rules: Sequence[Rule[CombinatorialClassType, CombinatorialObjectType]]
+    ):
         assert all(rule.is_equivalence() for rule in rules)
+        assert all(len(rule.children) == 1 for rule in rules)
         super().__init__(rules[0].strategy, rules[0].comb_class, rules[-1].children)
         self.rules = tuple(rules)
         self._constructor: Optional[DisjointUnion] = None
@@ -751,6 +853,12 @@ def eqv_path_rules(self) -> List[Tuple[CombinatorialClassType, Rule]]:
             curr = rule.children[0]
         return eqv_path_rules
 
+    def to_equivalence_rule(self) -> "EquivalenceRule":
+        raise NotImplementedError("You don't want to do that! I promise")
+
+    def to_reverse_rule(self, idx: int) -> "Rule":
+        raise NotImplementedError("You don't want to do that! I promise")
+
     def backward_map(
         self, objs: Tuple[Optional[CombinatorialObjectType], ...]
     ) -> Iterator[CombinatorialObjectType]:
@@ -809,12 +917,10 @@ def join(res, new):
 
 
 class ReverseRule(Rule[CombinatorialClassType, CombinatorialObjectType]):
-    """
-    A class for creating a reverse equivalence rule.
-    """
-
-    def __init__(self, rule: Rule, idx: int):
-        assert rule.is_two_way()
+    def __init__(
+        self, rule: Rule[CombinatorialClassType, CombinatorialObjectType], idx: int
+    ):
+        assert rule.is_reversible()
         super().__init__(
             rule.strategy,
             rule.children[idx],
@@ -844,6 +950,32 @@ def from_dict(cls, d: dict) -> "ReverseRule":
     def to_reverse_rule(self, idx: int) -> "Rule":
         raise NotImplementedError("You don't want to do that! I promise")
 
+    def shifts(self) -> Tuple[int, ...]:
+        original_shifts = self.original_rule.shifts()
+        pshift = -original_shifts[self.idx]
+        return (pshift,) + tuple(
+            s + pshift
+            for s in (s for i, s in enumerate(original_shifts) if i != self.idx)
+        )
+
+    def forest_key(
+        self, get_label: Callable[[CombinatorialClassType], int]
+    ) -> ForestRuleKey:
+        return ForestRuleKey(
+            get_label(self.comb_class),
+            tuple(map(get_label, self.children)),
+            self.shifts(),
+            RuleBucket.EQUIV if self.is_equivalence() else RuleBucket.REVERSE,
+        )
+
+    def get_equation(
+        self, get_function: Callable[[CombinatorialClassType], Function]
+    ) -> Eq:
+        try:
+            return super().get_equation(get_function)
+        except NotImplementedError:
+            return self.original_rule.get_equation(get_function)
+
     @property
     def constructor(
         self,
@@ -854,7 +986,6 @@ def constructor(
                 self.original_rule.comb_class,
                 self.original_rule.children,
             )
-
         return self._constructor
 
     @property
@@ -866,24 +997,26 @@ def formal_step(self) -> str:
     def backward_map(
         self, objs: Tuple[Optional[CombinatorialObjectType], ...]
     ) -> Iterator[CombinatorialObjectType]:
-        # We expect the object to be in only the first and in particular when
-        # mapped should be in the idx child.
-        if not isinstance(self.constructor, DisjointUnion) and len(self.children) != 1:
-            raise NotImplementedError("Cannot map for non equivalence rule.")
-        yield cast(
-            CombinatorialObjectType, self.original_rule.forward_map(objs[0])[self.idx]
-        )
+        if not len(self.original_rule.non_empty_children()) == 1:
+            raise NotImplementedError("Cannot map forward for non equivalence rule.")
+        assert all(obj is None for obj in objs[1:])
+        assert objs[0] is not None
+        res = self.original_rule.forward_map(objs[0])[self.idx]
+        assert res is not None
+        yield res
 
     def forward_map(
         self, obj: CombinatorialObjectType
     ) -> Tuple[Optional[CombinatorialObjectType], ...]:
-        # We map forward to the first child.
-        # It is assumed that obj is ONLY on this child
-        if not isinstance(self.constructor, DisjointUnion) and len(self.children) != 1:
-            raise NotImplementedError("Cannot map for non equivalence rule.")
-        return (next(self.original_rule.backward_map((obj,))),) + tuple(
-            None for _ in range(len(self.children) - 1)
-        )
+        if not len(self.original_rule.non_empty_children()) == 1:
+            raise NotImplementedError("Cannot map forward for non equivalence rule.")
+        objs: List[Optional[CombinatorialObjectType]] = [
+            None for _ in self.original_rule.children
+        ]
+        objs[self.idx] = obj
+        orig_res = list(self.original_rule.backward_map(tuple(objs)))
+        assert len(orig_res) == 1, "More than one results from back map"
+        return tuple(orig_res) + tuple(None for _ in range(len(self.children) - 1))
 
     def get_op_symbol(self):
         return self.strategy.get_reverse_op_symbol()
@@ -895,13 +1028,14 @@ class VerificationRule(AbstractRule[CombinatorialClassType, CombinatorialObjectT
     empty tuple if it applies, else None.
     """
 
+    # The signature is refined from AbstractRule
+    # pylint: disable=useless-super-delegation
     def __init__(
         self,
         strat: "VerificationStrategy[CombinatorialClassType,CombinatorialObjectType]",
         comb_class: CombinatorialClassType,
         children: Optional[Tuple[CombinatorialClassType, ...]] = None,
     ):
-        assert not children
         super().__init__(strat, comb_class, children)
 
     def to_jsonable(self) -> dict:
@@ -953,6 +1087,17 @@ def _ensure_level_objects(self, n: int) -> None:
     def pack(self) -> "StrategyPack":
         return self.strategy.pack(self.comb_class)
 
+    def forest_key(
+        self, get_label: Callable[[CombinatorialClassType], int]
+    ) -> ForestRuleKey:
+        return ForestRuleKey(
+            get_label(self.comb_class),
+            tuple(map(get_label, self.children)),
+            self.shifts(),
+            RuleBucket.VERIFICATION,
+        )
+
+    # pylint: disable=arguments-differ
     def get_equation(
         self,
         get_function: Callable[[CombinatorialClassType], Function],
@@ -961,6 +1106,10 @@ def get_equation(
         lhs_func = get_function(self.comb_class)
         return Eq(lhs_func, self.strategy.get_genf(self.comb_class, funcs))
 
+    @staticmethod
+    def get_eq_symbol() -> str:
+        return "↝"
+
     def random_sample_object_of_size(
         self, n: int, **parameters: int
     ) -> CombinatorialObjectType:
diff --git a/comb_spec_searcher/strategies/strategy.py b/comb_spec_searcher/strategies/strategy.py
index 0fa89504c..64874d93a 100644
--- a/comb_spec_searcher/strategies/strategy.py
+++ b/comb_spec_searcher/strategies/strategy.py
@@ -1,3 +1,4 @@
+# pylint: disable=too-many-lines
 """
 The Strategy class is essentially following the mantra of 'strategy' from the
 combinatorial explanation paper.
@@ -137,7 +138,6 @@ def __call__(
         self,
         comb_class: CombinatorialClassType,
         children: Tuple[CombinatorialClassType, ...] = None,
-        **kwargs,
     ) -> AbstractRule[CombinatorialClassType, CombinatorialObjectType]:
         """
         Return the rule formed by using the strategy.
@@ -172,6 +172,7 @@ def workable(self) -> bool:
         """
         return self._workable
 
+    @abc.abstractmethod
     def can_be_equivalent(self) -> bool:
         """
         Return True if every Rule returned with one non-empty child is an
@@ -186,6 +187,27 @@ def is_two_way(self, comb_class: CombinatorialClassType) -> bool:
         the rhs.
         """
 
+    @abc.abstractmethod
+    def is_reversible(self, comb_class: CombinatorialClassType) -> bool:
+        """
+        Return True if knowing the enumeration of the lhs and all bar one
+        of the rhs implies knowing the enumeration of all of the classes on
+        the rhs.
+        """
+
+    @abc.abstractmethod
+    def shifts(
+        self,
+        comb_class: CombinatorialClassType,
+        children: Optional[Tuple[CombinatorialClassType, ...]],
+    ) -> Tuple[int, ...]:
+        """
+        Return the reliance on the children needed for counting.
+
+        Positive shift indicate a reliance on shorted length than the one on the parent
+        while negative number indicate reliance on length greater length.
+        """
+
     @abc.abstractmethod
     def decomposition_function(
         self, comb_class: CombinatorialClassType
@@ -306,7 +328,6 @@ def __call__(
         self,
         comb_class: CombinatorialClassType,
         children: Tuple[CombinatorialClassType, ...] = None,
-        **kwargs,
     ) -> Rule[CombinatorialClassType, CombinatorialObjectType]:
         if children is None:
             children = self.decomposition_function(comb_class)
@@ -426,6 +447,23 @@ def can_be_equivalent() -> bool:
     def is_two_way(comb_class: CombinatorialClassType) -> bool:
         return True
 
+    @staticmethod
+    def is_reversible(comb_class: CombinatorialClassType) -> bool:
+        return True
+
+    def shifts(
+        self,
+        comb_class: CombinatorialClassType,
+        children: Optional[Tuple[CombinatorialClassType, ...]] = None,
+    ) -> Tuple[int, ...]:
+        if children is None:
+            children = self.decomposition_function(comb_class)
+            if children is None:
+                raise StrategyDoesNotApply("Strategy does not apply")
+        min_points = tuple(c.minimum_size_of_object() for c in children)
+        point_sum = sum(min_points)
+        return tuple(point_sum - mpoint for mpoint in min_points)
+
     def constructor(
         self,
         comb_class: CombinatorialClassType,
@@ -499,6 +537,21 @@ def can_be_equivalent() -> bool:
     def is_two_way(comb_class: CombinatorialClassType) -> bool:
         return True
 
+    @staticmethod
+    def is_reversible(comb_class: CombinatorialClassType) -> bool:
+        return True
+
+    def shifts(
+        self,
+        comb_class: CombinatorialClassType,
+        children: Optional[Tuple[CombinatorialClassType, ...]] = None,
+    ) -> Tuple[int, ...]:
+        if children is None:
+            children = self.decomposition_function(comb_class)
+            if children is None:
+                raise StrategyDoesNotApply("Strategy does not apply")
+        return tuple(0 for _ in children)
+
     def constructor(
         self,
         comb_class: CombinatorialClassType,
@@ -624,7 +677,6 @@ def __call__(
         self,
         comb_class: CombinatorialClassType,
         children: Tuple[CombinatorialClassType, ...] = None,
-        **kwargs,
     ) -> VerificationRule[CombinatorialClassType, CombinatorialObjectType]:
         if children is None:
             children = self.decomposition_function(comb_class)
@@ -640,6 +692,17 @@ def can_be_equivalent() -> bool:
     def is_two_way(comb_class: CombinatorialClassType) -> bool:
         return False
 
+    @staticmethod
+    def is_reversible(comb_class: CombinatorialClassType) -> bool:
+        return False
+
+    def shifts(
+        self,
+        comb_class: CombinatorialClassType,
+        children: Optional[Tuple[CombinatorialClassType, ...]] = None,
+    ) -> Tuple[int, ...]:
+        return ()
+
     def pack(self, comb_class: CombinatorialClassType) -> "StrategyPack":
         """
         Returns a StrategyPack that finds a proof tree for the comb_class in
@@ -692,9 +755,11 @@ def decomposition_function(
         self, comb_class: CombinatorialClassType
     ) -> Union[Tuple[CombinatorialClassType, ...], None]:
         """
-        A combinatorial class C is marked as verified by returning a rule
-        C -> (). This ensures that C is in a combinatorial specification as it
-        appears exactly once on the left hand side.
+        A verification rule usually returns a rule C -> () with no children that
+        guarantees that C is in a combinatorial specification.
+
+        In special cases, the rule could have children to mark dependencies of the
+        expansion on some other combinatorial classes.
 
         The function returns None if the verification strategy doesn't apply.
         """
@@ -739,6 +804,9 @@ def to_jsonable(self) -> dict:
         d.pop("workable")
         return d
 
+    def __repr__(self) -> str:
+        return f"{self.__class__.__name__}(ignore_parent={self.ignore_parent})"
+
 
 class AtomStrategy(VerificationStrategy[CombinatorialClass, CombinatorialObject]):
     """
@@ -889,7 +957,7 @@ class StrategyFactory(abc.ABC, Generic[CombinatorialClassType]):
 
     @abc.abstractmethod
     def __call__(
-        self, comb_class: CombinatorialClassType, **kwargs
+        self, comb_class: CombinatorialClassType
     ) -> Iterator[Union[AbstractRule, AbstractStrategy]]:
         """
         Returns the results of the strategy on a comb_class.
diff --git a/comb_spec_searcher/tree_searcher.py b/comb_spec_searcher/tree_searcher.py
index b73447827..c249a20ff 100644
--- a/comb_spec_searcher/tree_searcher.py
+++ b/comb_spec_searcher/tree_searcher.py
@@ -8,7 +8,7 @@
 from random import choice, shuffle
 from typing import Dict, FrozenSet, Iterator, List, Optional, Sequence, Set, Tuple
 
-from comb_spec_searcher.typing import RulesDict
+from comb_spec_searcher.typing import RuleKey, RulesDict
 
 __all__ = ("prune", "proof_tree_generator_dfs", "proof_tree_generator_bfs")
 
@@ -34,6 +34,30 @@ def nodes(self) -> Iterator["Node"]:
         for node in self.children:
             yield from node.nodes()
 
+    def rule_keys(self) -> Set[RuleKey]:
+        """
+        Yields all the rule keys in the proof tree.
+
+        We remove rules from recurring node.
+        """
+        rulekeys = (
+            (node.label, tuple(sorted(child.label for child in node.children)))
+            for node in self.nodes()
+        )
+        res: Dict[int, Tuple[int, ...]] = {}
+        for rulekey in rulekeys:
+            parent, children = rulekey
+            if parent in res:
+                if not res[parent]:
+                    # This was probably inserted because of a recursion node.
+                    # We replace it.
+                    res[parent] = children
+                else:
+                    assert not children or children == res[parent]
+            else:
+                res[parent] = children
+        return set(res.items())
+
     def __str__(self) -> str:
         return "".join(["(", str(self.label), *map(str, self.children), ")"])
 
diff --git a/comb_spec_searcher/typing.py b/comb_spec_searcher/typing.py
index 2e89aaff6..a20fcd07d 100644
--- a/comb_spec_searcher/typing.py
+++ b/comb_spec_searcher/typing.py
@@ -1,3 +1,4 @@
+import enum
 from typing import (
     TYPE_CHECKING,
     Callable,
@@ -29,6 +30,8 @@
     "Parameters",
     "ParametersMap",
     "RelianceProfile",
+    "RuleBucket",
+    "ForestRuleKey",
     "Objects",
     "ObjectsCache",
     "Terms",
@@ -59,8 +62,29 @@ class WorkPacket(NamedTuple):
     inferral: bool
 
 
+@enum.unique
+class RuleBucket(enum.Enum):
+    UNDEFINED = enum.auto()
+    VERIFICATION = enum.auto()
+    EQUIV = enum.auto()
+    NORMAL = enum.auto()
+    REVERSE = enum.auto()
+
+
 RuleKey = Tuple[int, Tuple[int, ...]]
 
+
+class ForestRuleKey(NamedTuple):
+    parent: int
+    children: Tuple[int, ...]
+    shifts: Tuple[int, ...]
+    bucket: RuleBucket
+
+    @property
+    def key(self) -> RuleKey:
+        return (self.parent, self.children)
+
+
 # From constructor
 Parameters = Tuple[int, ...]
 ParametersMap = Callable[[Parameters], Parameters]
diff --git a/example.py b/example.py
index 43b44dc7d..f5ea46659 100644
--- a/example.py
+++ b/example.py
@@ -47,10 +47,13 @@
     DisjointUnionStrategy,
     StrategyPack,
 )
+from comb_spec_searcher.bijection import ParallelSpecFinder
+from comb_spec_searcher.isomorphism import Bijection
 
 
 class Word(str, CombinatorialObject):
-    pass
+    def size(self):
+        return str.__len__(self)
 
 
 class AvoidingWithPrefix(CombinatorialClass[Word]):
@@ -207,7 +210,7 @@ def forward_map(
             if word[: len(child.prefix)] == child.prefix:
                 return (
                     tuple(None for _ in range(idx + 1))
-                    + (child,)
+                    + (word,)
                     + tuple(None for _ in range(len(children) - idx - 1))
                 )
 
@@ -250,7 +253,8 @@ def index_safe_to_remove_up_to(self, avoiding_with_prefix: AvoidingWithPrefix):
         )
         # safe will be the index of the prefix in which we can remove upto without
         # affecting the avoidance conditions
-        safe = max(0, len(prefix) - max(len(p) for p in patterns) + 1)
+        m = max(len(p) for p in patterns) if patterns else 1
+        safe = max(0, len(prefix) - m + 1)
         for i in range(safe, len(prefix)):
             end = prefix[i:]
             if any(end == patt[: len(end)] for patt in patterns):
@@ -356,3 +360,26 @@ def __repr__(self) -> str:
         random_word = spec.random_sample_object_of_size(n)
         print(random_word)
         print("Time to sample:", round(time.time() - start_time, 2), "seconds")
+
+    input(
+        '\nBijection example between "00" and "11" avoiding binary strings '
+        "(press any key to continue)"
+    )
+
+    specs = ParallelSpecFinder[AvoidingWithPrefix, Word, AvoidingWithPrefix, Word](
+        CombinatorialSpecificationSearcher(
+            AvoidingWithPrefix(Word(), ["00"], ["0", "1"]), pack
+        ),
+        CombinatorialSpecificationSearcher(
+            AvoidingWithPrefix(Word(), ["11"], ["0", "1"]), pack
+        ),
+    ).find()
+    assert specs is not None
+    spec1, spec2 = specs
+    bijection = Bijection.construct(spec1, spec2)
+    assert bijection is not None
+    for i in range(5):
+        for word in bijection.domain.generate_objects_of_size(i):
+            mapped_to = bijection.map(word)
+            assert bijection.inverse_map(mapped_to) == word
+            print(f"{word} -> {mapped_to}")
diff --git a/setup.py b/setup.py
index 860b36867..6bbbb5263 100755
--- a/setup.py
+++ b/setup.py
@@ -8,9 +8,17 @@ def read(fname):
     return open(os.path.join(os.path.dirname(__file__), fname)).read()
 
 
+def get_version(rel_path):
+    for line in read(rel_path).splitlines():
+        if line.startswith("__version__"):
+            delim = '"' if '"' in line else "'"
+            return line.split(delim)[1]
+    raise RuntimeError("Unable to find version string.")
+
+
 setup(
     name="comb_spec_searcher",
-    version="3.0.0",
+    version=get_version("comb_spec_searcher/__init__.py"),
     author="Permuta Triangle",
     author_email="permutatriangle@gmail.com",
     description="A library for performing combinatorial exploration.",
@@ -21,10 +29,10 @@ def read(fname):
         "Source": "https://github.com/PermutaTriangle/comb_spec_searcher",
         "Tracker": ("https://github.com/PermutaTriangle/comb_spec_searcher" "/issues"),
     },
-    packages=find_packages(),
+    packages=find_packages(exclude=["*.tests", "*.tests.*", "tests.*", "tests"]),
     package_data={"comb_spec_searcher": ["py.typed"]},
     long_description=read("README.rst"),
-    python_requires=">=3.6",
+    python_requires=">=3.7",
     include_package_data=True,
     classifiers=[
         "Development Status :: 5 - Production/Stable",
@@ -32,21 +40,21 @@ def read(fname):
         "Intended Audience :: Science/Research",
         "License :: OSI Approved :: BSD License",
         "Programming Language :: Python :: 3",
-        "Programming Language :: Python :: 3.6",
         "Programming Language :: Python :: 3.7",
         "Programming Language :: Python :: 3.8",
+        "Programming Language :: Python :: 3.9",
         "Programming Language :: Python :: Implementation :: CPython",
         "Programming Language :: Python :: Implementation :: PyPy",
         "Topic :: Education",
         "Topic :: Scientific/Engineering :: Mathematics",
     ],
     install_requires=[
-        "logzero==1.6.3",
-        "sympy==1.7.1",
+        "logzero==1.7.0",
+        "sympy==1.8",
         "psutil==5.8.0",
         "pympler==0.9",
         "requests==2.25.1",
-        "typing-extensions==3.7.4.3",
-        "tabulate==0.8.7",
+        "typing-extensions==3.10.0.0",
+        "tabulate==0.8.9",
     ],
 )
diff --git a/tox.ini b/tox.ini
index e9740e9ac..34b285da0 100644
--- a/tox.ini
+++ b/tox.ini
@@ -6,8 +6,8 @@
 [tox]
 envlist =
     flake8, mypy, pylint, black, tilescope
-    py{36,37,38,39},
-    pypy36
+    py{37,38,39},
+    pypy37
 
 [default]
 basepython=python3.8
@@ -15,16 +15,15 @@ basepython=python3.8
 [testenv]
 description = run test
 basepython =
-    py36: python3.6
     py37: python3.7
     py38: python3.8
     py39: python3.9
-    pypy36: pypy3
+    pypy37: pypy3
 deps =
-    pytest==6.2.1
+    pytest==6.2.4
     pytest-repeat==0.9.1
-    docutils==0.16
-    Pygments==2.7.3
+    docutils==0.17.1
+    Pygments==2.9.0
 commands = pytest
 
 [pytest]
@@ -46,7 +45,7 @@ description = run flake8 (linter)
 basepython = {[default]basepython}
 skip_install = True
 deps =
-    flake8==3.8.4
+    flake8==3.9.2
     flake8-isort==4.0.0
 commands =
     flake8 --isort-show-traceback comb_spec_searcher tests setup.py example.py
@@ -55,14 +54,16 @@ commands =
 description = run pylint (static code analysis)
 basepython = {[default]basepython}
 deps =
-    pylint==2.6.0
+    pylint==2.8.3
 commands = pylint comb_spec_searcher
 
 [testenv:mypy]
 description = run mypy (static type checker)
 basepython = {[default]basepython}
 deps =
-    mypy==0.790
+    mypy==0.902
+    types-requests==0.1.11
+    types-tabulate==0.1.0
 commands =
     mypy
 
@@ -70,5 +71,5 @@ commands =
 description = check that comply with autoformating
 basepython = {[default]basepython}
 deps =
-    black==20.8b1
+    black==21.6b0
 commands = black --check --diff .


Test patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 2909d2b60..427b4007f 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -25,20 +25,17 @@ jobs:
             toxenv: tilescope
             os: ubuntu-latest
 
-          - python: 3.6
-            toxenv: py36
-            os: ubuntu-latest
           - python: 3.7
             toxenv: py37
             os: ubuntu-latest
           - python: 3.8
             toxenv: py38
             os: ubuntu-latest
-          - python: 3.9-dev
+          - python: 3.9
             toxenv: py39
             os: ubuntu-latest
-          - python: pypy3
-            toxenv: pypy36
+          - python: pypy-3.7
+            toxenv: pypy37
             os: ubuntu-latest
 
           - python: 3.8
diff --git a/tests/test_bijections.py b/tests/test_bijections.py
new file mode 100644
index 000000000..067df1995
--- /dev/null
+++ b/tests/test_bijections.py
@@ -0,0 +1,100 @@
+import json
+from typing import Optional
+
+from comb_spec_searcher import (
+    AtomStrategy,
+    CombinatorialSpecificationSearcher,
+    StrategyPack,
+)
+from comb_spec_searcher.bijection import ParallelSpecFinder
+from comb_spec_searcher.isomorphism import Bijection
+from example import AvoidingWithPrefix, ExpansionStrategy, RemoveFrontOfPrefix, Word
+
+
+def find_bijection_between(
+    searcher1: CombinatorialSpecificationSearcher[AvoidingWithPrefix],
+    searcher2: CombinatorialSpecificationSearcher[AvoidingWithPrefix],
+) -> Optional[Bijection]:
+    specs = ParallelSpecFinder[AvoidingWithPrefix, Word, AvoidingWithPrefix, Word](
+        searcher1, searcher2
+    ).find()
+    if specs is not None:
+        s1, s2 = specs
+        return Bijection.construct(s1, s2)
+
+
+def get_word_searcher(avoid, alphabet):
+    pack = StrategyPack(
+        initial_strats=[RemoveFrontOfPrefix()],
+        inferral_strats=[],
+        expansion_strats=[[ExpansionStrategy()]],
+        ver_strats=[AtomStrategy()],
+        name=("Finding specification for words avoiding consecutive patterns."),
+    )
+    start_class = AvoidingWithPrefix("", avoid, alphabet)
+    searcher = CombinatorialSpecificationSearcher(start_class, pack)
+    return searcher
+
+
+def get_bijection(alphabet1, avoid1, alphabet2, avoid2):
+    s1, s2 = get_word_searcher(avoid1, alphabet1), get_word_searcher(avoid2, alphabet2)
+    return find_bijection_between(s1, s2)
+
+
+def assert_mappings(bijection):
+    assert bijection is not None
+    assert all(
+        set(bijection.map(w) for w in bijection.domain.generate_objects_of_size(i))
+        == set(bijection.codomain.generate_objects_of_size(i))
+        and set(
+            bijection.inverse_map(w)
+            for w in bijection.codomain.generate_objects_of_size(i)
+        )
+        == set(bijection.domain.generate_objects_of_size(i))
+        and all(
+            w == bijection.inverse_map(bijection.map(w))
+            for w in bijection.domain.generate_objects_of_size(i)
+        )
+        for i in range(10)
+    )
+
+
+def test_finding_bijection_and_map_in_words():
+    bijection = get_bijection(["0", "1"], ["00"], ["0", "1"], ["11"])
+    assert bijection is not None
+    assert_mappings(bijection)
+
+
+def test_finding_bijection_and_map_in_words2():
+    bijection = get_bijection(["a", "b"], ["aa"], ["0", "1"], ["11"])
+    assert bijection is not None
+    assert_mappings(bijection)
+
+
+def test_finding_bijection_and_map_in_words3():
+    bijection = get_bijection(["0", "1"], ["1011", "101"], ["0", "1"], ["010", "0101"])
+    assert bijection is not None
+    assert_mappings(bijection)
+
+
+def test_non_equal_classes():
+    bijection = get_bijection(["0", "1"], ["00"], ["0", "1"], ["101"])
+    assert bijection is None
+
+
+def test_constructing_bijection_non_eq_classes():
+    assert (
+        Bijection.construct(
+            get_word_searcher(["00"], ["0", "1"]).auto_search(),
+            get_word_searcher(["101"], ["0", "1"]).auto_search(),
+        )
+        is None
+    )
+
+
+def test_bijection_jsonable():
+    bijection = get_bijection(["0", "1"], ["1011", "101"], ["0", "1"], ["010", "0101"])
+    assert bijection is not None
+    assert_mappings(
+        Bijection.from_dict(json.loads(json.dumps(bijection.to_jsonable())))
+    )
diff --git a/tests/test_forest.py b/tests/test_forest.py
new file mode 100644
index 000000000..b4a3268d5
--- /dev/null
+++ b/tests/test_forest.py
@@ -0,0 +1,342 @@
+import itertools
+from typing import Dict, List, Union
+
+import pytest
+
+from comb_spec_searcher import CombinatorialSpecificationSearcher
+from comb_spec_searcher.rule_db.forest import Function, RuleDBForest, TableMethod
+from comb_spec_searcher.strategies.strategy import EmptyStrategy
+from comb_spec_searcher.typing import ForestRuleKey, RuleBucket
+from example import AvoidingWithPrefix, pack
+
+
+def assert_function_values(
+    function: Function, values: Union[Dict[int, int], List[int]]
+):
+    """
+    Check that the function values matches the given values.
+    """
+    if isinstance(values, dict):
+        values = [values[i] if i in values else 0 for i in range(max(values) + 1)]
+    failing_index = []
+    for n, (fv, gv) in enumerate(
+        itertools.zip_longest(function._value, values, fillvalue=0)
+    ):
+        if not (fv == gv):
+            failing_index.append(n)
+
+    if failing_index:
+        m = f"Function differs on index {failing_index}\n"
+        m += f"The function values are {function._value}\n"
+        m += f"The expected values are {values}"
+        raise AssertionError(m)
+
+
+class TestFunction:
+    def test_add_value(self):
+        f = Function()
+        assert f[0] == 0
+        assert f[4] == 0
+        f.increase_value(0)
+        assert f[0] == 1
+        f.increase_value(3)
+        assert f[4] == 0
+        f.increase_value(4)
+        assert f[0] == 1
+        assert f[1] == 0
+        assert f[2] == 0
+        assert f[3] == 1
+        assert f[4] == 1
+        assert f[5] == 0
+        assert f[6] == 0
+        with pytest.raises(ValueError):
+            assert sorted(f.preimage(0)) == []
+        assert sorted(f.preimage(1)) == [0, 3, 4]
+        assert sorted(f.preimage(2)) == []
+
+    def test_infinity(self):
+        f = Function()
+        f.increase_value(0)
+        f.increase_value(3)
+        f.increase_value(4)
+        f.set_infinite(3)
+        assert f[0] == 1
+        assert f[1] == 0
+        assert f[2] == 0
+        assert f[3] is None
+        assert f[4] == 1
+        assert f.preimage_gap(100) == 2
+        assert sorted(f.preimage(None)) == [3]
+        assert sorted(f.preimage(1)) == [0, 4]
+
+    def test_find_gap(self):
+        f = Function()
+        f.increase_value(0)
+        f.increase_value(0)
+        f.increase_value(0)
+        f.increase_value(0)
+        f.increase_value(1)
+        f.increase_value(2)
+        assert f[0] == 4
+        assert f[1] == 1
+        assert f[2] == 1
+        assert f[3] == 0
+        assert f[4] == 0
+        assert f[5] == 0
+        assert f[6] == 0
+        assert f.preimage_gap(1) == 2
+        assert f.preimage_gap(2) == 2
+        assert f.preimage_gap(3) == 5
+        with pytest.raises(ValueError):
+            f.preimage_gap(0)
+        with pytest.raises(ValueError):
+            f.preimage_gap(-1)
+
+    def test_find_gap2(self):
+        f = Function()
+        f.increase_value(2)
+        f.increase_value(3)
+        f.increase_value(4)
+        f.increase_value(5)
+        f.increase_value(5)
+        f.increase_value(5)
+        print(f)
+        print(f._preimage_count)
+        assert f[0] == 0
+        assert f[1] == 0
+        assert f[2] == 1
+        assert f[3] == 1
+        assert f[4] == 1
+        assert f[5] == 3
+        assert f.preimage_gap(1) == 2
+        assert f.preimage_gap(2) == 4
+        assert f.preimage_gap(3) == 4
+
+
+# Test of the table method
+
+
+def test_132_universe_pumping():
+    """
+    The universe consist of the rule of the usual 132 tree plus a dummy rule that is
+    useless.
+    """
+    tb = TableMethod()
+    rules = [
+        ForestRuleKey(0, (1, 2), (0, 0), RuleBucket.NORMAL),
+        ForestRuleKey(1, (), (), RuleBucket.VERIFICATION),  # Empty verif
+        ForestRuleKey(2, (3,), (0,), RuleBucket.NORMAL),  # point placement
+        ForestRuleKey(3, (4,), (0,), RuleBucket.NORMAL),  # row col sep
+        ForestRuleKey(4, (5, 0, 0), (0, 1, 1), RuleBucket.NORMAL),  # factor
+        ForestRuleKey(5, (), (), RuleBucket.VERIFICATION),  # point verif
+        ForestRuleKey(2, (6,), (2,), RuleBucket.UNDEFINED),  # dumb rule
+    ]
+    for rule in rules:
+        tb.add_rule_key(rule)
+    assert tb.function == {i: None for i in range(6)}
+    assert all(tb.is_pumping(c) for c in range(6))
+    assert not tb.is_pumping(6)
+    assert sorted(fk.key for fk in tb.pumping_subuniverse()) == [
+        (0, (1, 2)),
+        (1, tuple()),
+        (2, (3,)),
+        (3, (4,)),
+        (4, (5, 0, 0)),
+        (5, tuple()),
+    ]
+
+
+def test_132_universe_pumping_progressive():
+    """
+    The universe consist of the rule of the usual 132 tree plus a dummy rule that is
+    useless.
+
+    We add rule progressively and make sure the function is always up to date.
+    """
+    tb = TableMethod()
+
+    # Point insertion
+    tb.add_rule_key(ForestRuleKey(0, (1, 2), (0, 0), RuleBucket.NORMAL))
+    assert tb.function == {}
+
+    # Empty verif
+    tb.add_rule_key(ForestRuleKey(1, (), (), RuleBucket.VERIFICATION))
+    assert tb.function == {1: None}
+
+    # Point placement
+    tb.add_rule_key(ForestRuleKey(2, (3,), (0,), RuleBucket.NORMAL))
+    assert tb.function == {1: None}
+
+    # Row col sep
+    tb.add_rule_key(ForestRuleKey(3, (4,), (0,), RuleBucket.NORMAL))
+    assert tb.function == {1: None}
+
+    # Point verif
+    tb.add_rule_key(ForestRuleKey(5, (), (), RuleBucket.VERIFICATION))
+    assert tb.function == {1: None, 5: None}
+
+    # Dumb rule
+    tb.add_rule_key(ForestRuleKey(2, (6,), (-2,), RuleBucket.UNDEFINED))
+    assert tb.function == {1: None, 5: None}
+
+    # Dumb rule. This will pump 2 and 0 a little bit
+    tb.add_rule_key(ForestRuleKey(2, (7,), (2,), RuleBucket.UNDEFINED))
+    assert tb.function == {0: 2, 1: None, 2: 2, 5: None}
+
+    # Factor
+    tb.add_rule_key(ForestRuleKey(4, (5, 0, 0), (0, 1, 1), RuleBucket.NORMAL))
+    assert tb.function == {i: None for i in range(6)}
+
+
+def test_universe_not_pumping():
+    tb = TableMethod()
+    rules = [
+        ForestRuleKey(0, (1, 2), (0, 0), RuleBucket.NORMAL),  # point insertion
+        ForestRuleKey(5, (), (), RuleBucket.VERIFICATION),  # point verif
+        ForestRuleKey(2, (3,), (0,), RuleBucket.NORMAL),  # point placement
+        ForestRuleKey(3, (4,), (0,), RuleBucket.NORMAL),  # row col sep
+        ForestRuleKey(4, (5, 0, 0), (0, 1, 1), RuleBucket.NORMAL),  # factor
+    ]
+    for rule in rules:
+        tb.add_rule_key(rule)
+    assert tb.function == {2: 1, 3: 1, 4: 1, 5: None}
+
+
+def test_segmented():
+    """The segemented forest."""
+    tb = TableMethod()
+    tb.add_rule_key(ForestRuleKey(0, (1, 2), (0, 0), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(1, (4, 14), (0, 0), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(2, tuple(), tuple(), RuleBucket.UNDEFINED))
+    assert tb.function == {2: None}
+
+    tb.add_rule_key(ForestRuleKey(3, (16, 5), (1, 0), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(4, tuple(), tuple(), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(5, tuple(), tuple(), RuleBucket.UNDEFINED))
+    assert tb.function == {2: None, 3: 1, 4: None, 5: None}
+
+    # Induced a gap size change
+    tb.add_rule_key(ForestRuleKey(6, (7, 5, 17), (2, 1, 1), RuleBucket.UNDEFINED))
+    assert tb.function == {2: None, 3: 1, 4: None, 5: None, 6: 1}
+
+    tb.add_rule_key(ForestRuleKey(16, (6,), (0,), RuleBucket.UNDEFINED))
+    assert tb.function == {2: None, 3: 2, 4: None, 5: None, 6: 1, 16: 1}
+
+    tb.add_rule_key(ForestRuleKey(7, tuple(), tuple(), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(8, (9, 5), (1, 0), RuleBucket.UNDEFINED))
+    assert tb.function == {
+        2: None,
+        3: 2,
+        4: None,
+        5: None,
+        6: 1,
+        7: None,
+        8: 1,
+        16: 1,
+    }
+
+    tb.add_rule_key(ForestRuleKey(12, (20, 5), (-1, 0), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(20, (13,), (0,), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(13, (15, 2, 5), (-1, 1, 0), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(15, (1,), (0,), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(14, (3,), (0,), RuleBucket.UNDEFINED))
+    assert tb.function == {
+        0: 2,
+        1: 2,
+        2: None,
+        3: 2,
+        4: None,
+        5: None,
+        6: 1,
+        7: None,
+        8: 1,
+        13: 1,
+        14: 2,
+        15: 2,
+        16: 1,
+        20: 1,
+    }
+
+    tb.add_rule_key(ForestRuleKey(18, (8,), (0,), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(11, (12, 18), (0, 0), RuleBucket.UNDEFINED))
+    assert tb.function == {
+        0: 2,
+        1: 2,
+        2: None,
+        3: 2,
+        4: None,
+        5: None,
+        6: 1,
+        7: None,
+        8: 1,
+        13: 1,
+        14: 2,
+        15: 2,
+        16: 1,
+        18: 1,
+        20: 1,
+    }
+
+    tb.add_rule_key(ForestRuleKey(17, (8,), (0,), RuleBucket.UNDEFINED))
+    assert tb.function == {
+        0: 3,
+        1: 3,
+        2: None,
+        3: 3,
+        4: None,
+        5: None,
+        6: 2,
+        7: None,
+        8: 1,
+        11: 1,
+        12: 1,
+        13: 2,
+        14: 3,
+        15: 3,
+        16: 2,
+        17: 1,
+        18: 1,
+        20: 2,
+    }
+
+    tb.add_rule_key(ForestRuleKey(9, (0, 19), (0, 0), RuleBucket.UNDEFINED))
+    tb.add_rule_key(ForestRuleKey(10, (5, 11), (0, 1), RuleBucket.UNDEFINED))
+    assert tb.function == {
+        0: 3,
+        1: 3,
+        2: None,
+        3: 3,
+        4: None,
+        5: None,
+        6: 2,
+        7: None,
+        8: 1,
+        10: 2,
+        11: 1,
+        12: 1,
+        13: 2,
+        14: 3,
+        15: 3,
+        16: 2,
+        17: 1,
+        18: 1,
+        20: 2,
+    }
+
+    tb.add_rule_key(ForestRuleKey(19, (10,), (0,), RuleBucket.UNDEFINED))
+    assert tb.function == {i: None for i in range(21)}
+    assert all(tb.is_pumping(c) for c in range(21))
+
+
+# Test of the extractor
+
+
+def test_emtpy_rule_first():
+    """
+    When recomputing the rules. The forest extractor should favor the empty rule.
+    """
+    empty_class = AvoidingWithPrefix("ab", ["ab"], ["a", "b"])
+    assert empty_class.is_empty()
+    css = CombinatorialSpecificationSearcher(empty_class, pack, ruledb=RuleDBForest())
+    spec = css.auto_search()
+    assert isinstance(spec.rules_dict[empty_class].strategy, EmptyStrategy)
diff --git a/tests/test_pickle.py b/tests/test_pickle.py
index 3e6e9b428..cba8ac4c7 100644
--- a/tests/test_pickle.py
+++ b/tests/test_pickle.py
@@ -30,3 +30,23 @@ def test_pickle_queue_not_exhausted():
     new_queue = pickle.loads(pickle.dumps(queue))
     assert new_queue == queue
     assert list(queue) == list(new_queue)
+
+
+def tests_pickling_css():
+    alphabet = ["a", "b"]
+    start_class = AvoidingWithPrefix("", ["ababa", "babb"], alphabet)
+    searcher = CombinatorialSpecificationSearcher[AvoidingWithPrefix](start_class, pack)
+    spec = searcher.auto_search()
+    new_searcher = pickle.loads(pickle.dumps(searcher))
+    assert searcher == new_searcher
+    assert spec == new_searcher.get_specification()
+    assert spec == new_searcher.auto_search(max_expansion_time=2)
+
+
+def test_pickle_specification():
+    alphabet = ["a", "b"]
+    start_class = AvoidingWithPrefix("", ["ababa", "babb"], alphabet)
+    searcher = CombinatorialSpecificationSearcher(start_class, pack)
+    spec = searcher.auto_search()
+    spec.count_objects_of_size(10)
+    pickle.dumps(spec)
diff --git a/tests/test_rule.py b/tests/test_rule.py
index 2dd6b6a17..d7257c2ff 100644
--- a/tests/test_rule.py
+++ b/tests/test_rule.py
@@ -1,15 +1,88 @@
-from example import AvoidingWithPrefix, ExpansionStrategy
+import sympy
+
+from comb_spec_searcher.strategies.rule import EquivalenceRule
+from example import AvoidingWithPrefix, ExpansionStrategy, RemoveFrontOfPrefix
 
 
 def test_reverse_equivalence():
     comb_class = AvoidingWithPrefix("aaa", ["aaaa"], "a", False)
     rule = ExpansionStrategy()(comb_class)
+    eqv_rule = rule.to_equivalence_rule()
     assert len(rule.non_empty_children()) == 1
     equiv_then_reverse = rule.to_equivalence_rule().to_reverse_rule(0)
     assert equiv_then_reverse.is_equivalence()
     assert equiv_then_reverse.children == (comb_class,)
     assert equiv_then_reverse.comb_class == rule.non_empty_children()[0]
+    assert isinstance(equiv_then_reverse, EquivalenceRule)
     reverse_then_equiv = rule.to_reverse_rule(0).to_equivalence_rule()
     assert reverse_then_equiv.is_equivalence()
     assert reverse_then_equiv.children == (comb_class,)
     assert reverse_then_equiv.comb_class == rule.non_empty_children()[0]
+    assert isinstance(reverse_then_equiv, EquivalenceRule)
+    assert (
+        eqv_rule.formal_step == "Either just the prefix, or append a"
+        " letter from the alphabet but only child and index 0 is non-empty"
+    )
+    assert reverse_then_equiv.formal_step == f"reverse of '{eqv_rule.formal_step}'"
+    assert equiv_then_reverse.formal_step == f"reverse of '{eqv_rule.formal_step}'"
+
+
+def test_complement_rule():
+    comb_class = AvoidingWithPrefix("aaa", ["aaaa"], "abc", False)
+    rule = ExpansionStrategy()(comb_class)
+    all_comb_classes = (rule.comb_class,) + rule.children
+    x = sympy.var("x")
+
+    def get_function(w):
+        i = all_comb_classes.index(w)
+        return sympy.Function(f"F_{i}")(x)
+
+    reverse0 = rule.to_reverse_rule(0)
+    reverse1 = rule.to_reverse_rule(1)
+    reverse2 = rule.to_reverse_rule(2)
+    reverse3 = rule.to_reverse_rule(3)
+    for i in range(4):
+        rule.sanity_check(i)
+        reverse0.sanity_check(i)
+        reverse1.sanity_check(i)
+        reverse2.sanity_check(i)
+        reverse3.sanity_check(i)
+    assert rule.get_equation(get_function) == sympy.sympify(
+        "Eq(F_0(x), F_1(x) + F_2(x) + F_3(x) + F_4(x))"
+    )
+    assert reverse0.get_equation(get_function) == sympy.sympify(
+        "Eq(F_1(x), F_0(x) - F_2(x) - F_3(x) - F_4(x))"
+    )
+    assert reverse1.get_equation(get_function) == sympy.sympify(
+        "Eq(F_2(x), F_0(x) - F_1(x) - F_3(x) - F_4(x))"
+    )
+    assert reverse2.get_equation(get_function) == sympy.sympify(
+        "Eq(F_3(x), F_0(x) - F_1(x) - F_2(x) - F_4(x))"
+    )
+    assert reverse3.get_equation(get_function) == sympy.sympify(
+        "Eq(F_4(x), F_0(x) - F_1(x) - F_2(x) - F_3(x))"
+    )
+
+
+def test_quotient_rule():
+    comb_class = AvoidingWithPrefix("aaa", ["aa"], "a", False)
+    rule = RemoveFrontOfPrefix()(comb_class)
+
+    all_comb_classes = (rule.comb_class,) + rule.children
+    x = sympy.var("x")
+
+    def get_function(w):
+        i = all_comb_classes.index(w)
+        return sympy.Function(f"F_{i}")(x)
+
+    reverse0 = rule.to_reverse_rule(0)
+    reverse1 = rule.to_reverse_rule(1)
+    assert rule.get_equation(get_function) == sympy.sympify(
+        "Eq(F_0(x), F_1(x) * F_2(x))"
+    )
+    assert reverse0.get_equation(get_function) == sympy.sympify(
+        "Eq(F_1(x), F_0(x) / F_2(x))"
+    )
+    assert reverse1.get_equation(get_function) == sympy.sympify(
+        "Eq(F_2(x), F_0(x) / F_1(x))"
+    )
diff --git a/tests/test_searcher.py b/tests/test_searcher.py
index c03160b9b..f32100dbb 100644
--- a/tests/test_searcher.py
+++ b/tests/test_searcher.py
@@ -1,7 +1,10 @@
 import pytest
 
 from comb_spec_searcher import CombinatorialSpecificationSearcher
-from comb_spec_searcher.exception import NoMoreClassesToExpandError
+from comb_spec_searcher.exception import (
+    NoMoreClassesToExpandError,
+    SpecificationNotFound,
+)
 from example import AvoidingWithPrefix, pack
 
 
@@ -14,3 +17,24 @@ def test_do_level():
         while True:
             searcher.do_level()
     assert all(searcher.classqueue.queue_sizes)
+
+
+def test_no_spec_exception():
+    alphabet = ["a", "b"]
+    start_class = AvoidingWithPrefix("", ["aabb", "bbbbab"], alphabet)
+    searcher = CombinatorialSpecificationSearcher(start_class, pack)
+    for _ in range(2):
+        searcher.do_level()
+    assert not searcher.ruledb.has_specification()
+    with pytest.raises(SpecificationNotFound):
+        searcher.get_specification()
+    with pytest.raises(SpecificationNotFound):
+        searcher.ruledb.get_specification_rules()
+
+
+def test_iterative():
+    alphabet = ["a", "b"]
+    start_class = AvoidingWithPrefix("", ["a"], alphabet)
+    it_pack = pack.make_iterative("iterative")
+    searcher = CombinatorialSpecificationSearcher(start_class, it_pack)
+    searcher.auto_search()
diff --git a/tests/test_specification.py b/tests/test_specification.py
index 59bea7eb4..44cbae018 100644
--- a/tests/test_specification.py
+++ b/tests/test_specification.py
@@ -6,6 +6,7 @@
     CombinatorialSpecification,
     CombinatorialSpecificationSearcher,
 )
+from comb_spec_searcher.rule_db import RuleDBForest, RuleDBForgetStrategy
 from comb_spec_searcher.utils import taylor_expand
 from example import AvoidingWithPrefix, Word, pack
 
@@ -88,8 +89,23 @@ def test_random_sample(specification):
 def test_forget_ruledb():
     alphabet = ["a", "b"]
     start_class = AvoidingWithPrefix("", ["ababa", "babb"], alphabet)
-    searcher = CombinatorialSpecificationSearcher(start_class, pack, ruledb="forget")
-    return searcher.auto_search()
+    ruledb = RuleDBForgetStrategy()
+    searcher = CombinatorialSpecificationSearcher(start_class, pack, ruledb=ruledb)
+    spec = searcher.auto_search()
+    expected_count = [1, 2, 4, 8, 15, 27, 48, 87, 157, 283]
+    count = [spec.count_objects_of_size(n) for n in range(10)]
+    assert count == expected_count
+
+
+def test_forest_ruledb():
+    alphabet = ["a", "b"]
+    start_class = AvoidingWithPrefix("", ["ababa", "babb"], alphabet)
+    ruledb = RuleDBForest()
+    searcher = CombinatorialSpecificationSearcher(start_class, pack, ruledb=ruledb)
+    spec = searcher.auto_search()
+    expected_count = [1, 2, 4, 8, 15, 27, 48, 87, 157, 283]
+    count = [spec.count_objects_of_size(n) for n in range(10)]
+    assert count == expected_count
 
 
 def test_sancheck(specification):