Repository: PermutaTriangle/comb_spec_searcher (language: python)

Issue:
"not really isomorphism but something"



Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index ad2d2964c..d5b1fa042 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,9 +4,29 @@ 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).
 
-## [2.4.0]
+## [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.
+- 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.
+
+### Changed
+- Specification are now built using a set of rules
+- The json format of a spec is based storing the json format of its rule
+- 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.
+
+### Fixed
+- `forward_map` of `EquivalencePathRule`
+- `all_specifications` method of `RuleDB`
+
+## [2.4.0] - 2020-11-11
 ### Changed
-- Computation of terms in the constructor is now on a per size basis. The value for each possible paratemeters is computed at once in the new `get_terms` function of constructor. If you do not implement your own constructor this should have no effect on your code. The old `generate_object_of_size` method of the constructor is also replaced by a `get_objects` method that returns the objects for each possible combination of parameters.
+- Computation of terms in the constructor is now on a per size basis. The value for each possible parameters is computed at once in the new `get_terms` function of constructor. If you do not implement your own constructor this should have no effect on your code. The old `generate_object_of_size` method of the constructor is also replaced by a `get_objects` method that returns the objects for each possible combination of parameters.
 
 ## [2.3.0] - 2020-10-28
 ### Added
diff --git a/README.rst b/README.rst
index f118e8d22..ff634fddb 100644
--- a/README.rst
+++ b/README.rst
@@ -132,12 +132,13 @@ if the same class appears multiple times.
    ...                     hash(self.alphabet) + hash(self.just_prefix))
    ...
    ...     def __str__(self):
+   ...         prefix = self.prefix if self.prefix else '""'
    ...         if self.just_prefix:
-   ...             return "The word {}".format(self.prefix)
+   ...             return "The word {}".format(prefix)
    ...         return ("Words over {{{}}} avoiding {{{}}} with prefix {}"
    ...                 "".format(", ".join(l for l in self.alphabet),
    ...                           ", ".join(p for p in self.patterns),
-   ...                           self.prefix if self.prefix else '""'))
+   ...                           prefix))
    ...
    ...     def __repr__(self):
    ...         return "AvoidingWithPrefix({}, {}, {}".format(repr(self.prefix),
diff --git a/comb_spec_searcher/comb_spec_searcher.py b/comb_spec_searcher/comb_spec_searcher.py
index f6ba9bcbd..835d31b1d 100644
--- a/comb_spec_searcher/comb_spec_searcher.py
+++ b/comb_spec_searcher/comb_spec_searcher.py
@@ -25,12 +25,7 @@
 from logzero import logger
 from sympy import Eq, Function, var
 
-from comb_spec_searcher.typing import (
-    CombinatorialClassType,
-    CSSstrategy,
-    SpecificationClassesAndStrats,
-    WorkPacket,
-)
+from comb_spec_searcher.typing import CombinatorialClassType, CSSstrategy, WorkPacket
 
 from .class_db import ClassDB
 from .class_queue import DefaultQueue
@@ -43,6 +38,7 @@
 from .rule_db import RuleDB, RuleDBForgetStrategy
 from .rule_db.base import RuleDBBase
 from .specification import CombinatorialSpecification
+from .specification_extrator import SpecificationRuleExtractor
 from .strategies import (
     AbstractStrategy,
     StrategyFactory,
@@ -50,6 +46,7 @@
     VerificationRule,
 )
 from .strategies.rule import AbstractRule
+from .tree_searcher import Node
 from .utils import (
     cssiteratortimer,
     cssmethodtimer,
@@ -275,19 +272,6 @@ def _expand_class_with_strategy(
                         "NotImplementedError: %s",
                         e,
                     )
-            if any(
-                self.ruledb.are_equivalent(start_label, elabel) for elabel in end_labels
-            ):
-                # This says comb_class = comb_class, so we skip it, but mark
-                # every other class as empty.
-                for elabel in end_labels:
-                    if not self.ruledb.are_equivalent(start_label, elabel):
-                        self._add_empty_rule(elabel)
-                if self.debug:
-                    for elabel, child in zip(end_labels, children):
-                        if not self.ruledb.are_equivalent(start_label, elabel):
-                            assert child.is_empty()
-
             yield start_label, tuple(end_labels), rule
 
     def _add_rule(
@@ -633,7 +617,7 @@ def auto_search(self, **kwargs) -> CombinatorialSpecification:
                 round(max_expansion_time, 2),
             )
 
-    def _auto_search_rules(self) -> SpecificationClassesAndStrats:
+    def _auto_search_rules(self) -> Iterator[AbstractRule]:
         """
         A basic auto search for returning equivalence paths and rules.
 
@@ -706,44 +690,41 @@ def get_specification(
         The function will return None if no such CombinatorialSpecification
         exists in the universe.
         """
-        spec = self._get_specification_rules(minimization_time_limit, smallest)
-        if spec is None:
+        rules = self._get_specification_rules(minimization_time_limit, smallest)
+        if rules is None:
             return None
         start_class = self.classdb.get_class(self.start_label)
-        strategies, comb_class_eqv_paths = spec
         logger.info("Creating a specification.")
-        return CombinatorialSpecification(start_class, strategies, comb_class_eqv_paths)
+        return CombinatorialSpecification(start_class, rules)
 
-    @cssmethodtimer("get specification")
     def _get_specification_rules(
         self, minimization_time_limit: float = 10, smallest: bool = False
-    ) -> Optional[SpecificationClassesAndStrats]:
-        """
-        Returns the equivalence paths needed to create a
-        CombinatorialSpecification, if one exists.
-
-        The minimization_time_limit only applies when smallest is 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()
 
-        The function will return None if no such CombinatorialSpecification
-        exists in the universe.
+    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")
-                rules, eqv_paths = self.ruledb.get_smallest_specification(
-                    self.start_label
-                )
+                node = self.ruledb.get_smallest_specification(self.start_label)
             else:
-                rules, eqv_paths = self.ruledb.get_specification_rules(
+                node = self.ruledb.find_specification(
                     self.start_label,
                     minimization_time_limit=minimization_time_limit,
                     iterative=self.iterative,
                 )
         except SpecificationNotFound:
             return None
-        comb_class_eqv_paths = [
-            list(map(self.classdb.get_class, path)) for path in eqv_paths
-        ]
-        strategies = [(self.classdb.get_class(label), rule) for label, rule in rules]
-        return strategies, comb_class_eqv_paths
+        return node
diff --git a/comb_spec_searcher/equiv_db.py b/comb_spec_searcher/equiv_db.py
index 2429972be..1a2c25ae7 100644
--- a/comb_spec_searcher/equiv_db.py
+++ b/comb_spec_searcher/equiv_db.py
@@ -10,8 +10,10 @@
 integer, that the classdb gives.
 """
 
-from collections import deque
-from typing import Deque, Dict, FrozenSet, List, Set, Union
+from collections import defaultdict, deque
+from typing import Deque, Dict, Iterator, List, Set, Tuple
+
+from .utils import cssmethodtimer
 
 
 class EquivalenceDB:
@@ -33,7 +35,10 @@ def __init__(self) -> None:
         self.parents: Dict[int, int] = {}
         self.weights: Dict[int, int] = {}
         self.verified_roots: Set[int] = set()
-        self.vertices: Set[FrozenSet[int]] = set()
+        self.vertices: Dict[int, Set[int]] = defaultdict(set)
+        self._one_way_vertices: Dict[int, Set[int]] = defaultdict(set)
+        self.func_times: Dict[str, float] = defaultdict(float)
+        self.func_calls: Dict[str, int] = defaultdict(int)
 
     def __eq__(self, other: object) -> bool:
         """Check if all information stored is the same."""
@@ -66,24 +71,100 @@ def __getitem__(self, comb_class: int) -> int:
             self.parents[ancestor] = root
         return root
 
-    def union(self, t1: int, t2: int) -> None:
-        """Find sets containing t1 and t2 and merge them."""
-        if t1 == t2:
+    def get_one_way_vertices(self) -> Dict[int, Set[int]]:
+        """
+        Return the adjacency table where vertices are strong components,
+        and edges are those that connect some vertex in a strong component
+        to the other.
+
+        Note: this might not be accurate, you should use the connect_cycle
+        method first if needed.
+        """
+        res: Dict[int, Set[int]] = defaultdict(set)
+        for start, ends in self._one_way_vertices.items():
+            start = self[start]
+            for end in ends:
+                end = self[end]
+                if start != end:
+                    res[start].add(end)
+        self._one_way_vertices = res
+        return res
+
+    def add_one_way_edge(self, label: int, other_label: int) -> None:
+        """Add an edge from label to other label in the directed graph."""
+        self._add_edge(label, other_label)
+        self._one_way_vertices[self[label]].add(self[other_label])
+
+    @cssmethodtimer("find_paths")
+    def connect_cycles(self):
+        """Look for cycles using one way edges that have been added."""
+        one_way_vertices = self.get_one_way_vertices()
+        stack: List[Tuple[int, ...]] = list()
+        stack.extend([(label,) for label in one_way_vertices.keys()])
+        visited = set()
+        while stack:
+            path = stack.pop()
+            end = path[-1]
+            if end in visited:
+                continue
+            visited.add(end)
+            for new_end in one_way_vertices[end]:
+                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
+                if new_end not in path:
+                    # new end is not in the path, so we haven't tried looking for
+                    # cycles from here.
+                    stack.append(
+                        path + (new_end,)
+                    )  # appending left makes this a depth first search
+
+    def add_two_way_edge(self, label: int, other_label: int) -> None:
+        """Add a two way edge to the directed graph."""
+        self._add_edge(label, other_label)
+        self._add_edge(other_label, label)
+        self._set_equivalent(label, other_label)
+
+    def _add_edge(self, label: int, other_label: int) -> None:
+        """Add an edge from label to other_label."""
+        if label == other_label:
             return
-        self.vertices.add(frozenset((t1, t2)))
-        verified = self.is_verified(t1) or self.is_verified(t2)
-        roots = [self[t1], self[t2]]
+        self.vertices[label].add(other_label)
+
+    def find_paths(self, label: int, other_label: int) -> Iterator[Tuple[int, ...]]:
+        """Return an iterator of paths starting at other_label and ending at label."""
+        dequeue: Deque[Tuple[int, ...]] = deque()
+        dequeue.append((other_label,))
+        while dequeue:
+            path = dequeue.popleft()
+            end = path[-1]
+            if self.equivalent(end, label):
+                yield path
+                continue
+            for new_end in self.vertices[end]:
+                if new_end in path:
+                    continue
+                dequeue.append(path + (new_end,))
+
+    def _set_equivalent(self, label: int, other_label: int) -> None:
+        """Set the two labels as equivalent in UnionFind"""
+        verified = self.is_verified(label) or self.is_verified(other_label)
+        roots = [self[label], self[other_label]]
         heaviest = max([(self.weights[r], r) for r in roots])[1]
         for r in roots:
             if r != heaviest:
                 self.weights[heaviest] += self.weights[r]
                 self.parents[r] = heaviest
         if verified:
-            self.set_verified(t1)
+            self.set_verified(label)
 
-    def equivalent(self, t1: int, t2: int) -> bool:
-        """Return True if t1 and t2 are equivalent, False otherwise."""
-        return self[t1] == self[t2]
+    def equivalent(self, label: int, other_label: int) -> bool:
+        """Return True if label and other_label are equivalent, False otherwise."""
+        return self[label] == self[other_label]
 
     def set_verified(self, comb_class: int) -> None:
         """Update database that combinatorial classes equivalent to comb_class
@@ -103,65 +184,28 @@ def equivalent_set(self, comb_class: int) -> Set[int]:
                 equivalent_classes.add(t)
         return equivalent_classes
 
-    def find_path(self, comb_class: int, other_comb_class: int) -> List[int]:
+    def find_path(self, comb_class: int, other_comb_class: int) -> Tuple[int, ...]:
         """
-        BFS for shortest path.
+        BFS for shortest path from comb_class to other_comb_class.
 
         Used to find shortest explanation of why things are equivalent.
         """
         if not self.equivalent(comb_class, other_comb_class):
             raise KeyError("The classes given are not equivalent.")
-        if comb_class == other_comb_class:
-            return [comb_class]
-        equivalent_comb_classes: Dict[int, int] = {}
-        reverse_map: Dict[int, int] = {}
-
-        for x in self.parents:
-            n = len(equivalent_comb_classes)
-            if self.equivalent(comb_class, x):
-                equivalent_comb_classes[x] = n
-                reverse_map[n] = x
-
-        adjacency_list: List[List[int]] = [
-            [] for i in range(len(equivalent_comb_classes))
-        ]
-        for (t1, t2) in self.vertices:
-            if self.equivalent(t1, comb_class):
-                u = equivalent_comb_classes[t1]
-                v = equivalent_comb_classes[t2]
-                adjacency_list[u].append(v)
-                adjacency_list[v].append(u)
-
-        dequeue: Deque[int] = deque()
-
-        start = equivalent_comb_classes[comb_class]
-        end = equivalent_comb_classes[other_comb_class]
-
-        n = len(equivalent_comb_classes)
-
-        dequeue.append(start)
-        visited = [False for i in range(n)]
-        neighbour: List[Union[int, None]] = [None for i in range(n)]
+
+        dequeue: Deque[Tuple[int, ...]] = deque()
+        dequeue.append((comb_class,))
+        visited: Set[int] = set()
         while dequeue:
-            u = dequeue.popleft()
-            if u == end:
+            path = dequeue.popleft()
+            end = path[-1]
+            if end == other_comb_class:
                 break
-            if visited[u]:
+            if end in visited:
                 continue
-            visited[u] = True
-
-            for v in adjacency_list[u]:
-                if not visited[v]:
-                    dequeue.append(v)
-                    neighbour[v] = u
-
-        path = [reverse_map[end]]
-        while neighbour[end] is not None:
-            assert isinstance(end, int), "something went wrong"
-            a = neighbour[end]
-            assert isinstance(a, int), "something went wrong"
-            t = reverse_map[a]
-            path.append(t)
-            end = a
-
-        return path[::-1]
+            for new_end in self.vertices[end]:
+                if new_end in path:
+                    continue
+                dequeue.append(path + (new_end,))
+            visited.add(path[-1])
+        return path
diff --git a/comb_spec_searcher/isomorphism.py b/comb_spec_searcher/isomorphism.py
new file mode 100644
index 000000000..7a39f786f
--- /dev/null
+++ b/comb_spec_searcher/isomorphism.py
@@ -0,0 +1,273 @@
+from typing import TYPE_CHECKING, Callable, Dict, List, Optional, Set, Tuple
+
+from .combinatorial_class import CombinatorialClass, CombinatorialObject
+from .strategies.rule import AbstractRule, Rule
+
+if TYPE_CHECKING:
+    from .specification import CombinatorialSpecification
+
+
+AtomEquals = Callable[["CombinatorialClass", "CombinatorialClass"], bool]
+
+
+def _default_equals(c1: CombinatorialClass, c2: CombinatorialClass) -> bool:
+    return c1 == c2
+
+
+class Isomorphism:
+    """Isomorphism checker."""
+
+    _INVALID, _UNKNOWN, _VALID = range(-1, 2)
+
+    @classmethod
+    def check(
+        cls,
+        spec1: "CombinatorialSpecification",
+        spec2: "CombinatorialSpecification",
+        eq: AtomEquals = _default_equals,
+    ) -> bool:
+        """Check if two specs are isomorphic."""
+        return cls(spec1, spec2, eq).are_isomorphic()
+
+    def __init__(
+        self,
+        spec1: "CombinatorialSpecification",
+        spec2: "CombinatorialSpecification",
+        eq: AtomEquals = _default_equals,
+    ) -> 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)
+
+    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]:
+        """Get order map of corresponding nodes."""
+        return self._order_map[(c1, c2)]
+
+    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
+        )
+        rule1, rule2 = self._rules1[node1], self._rules2[node2]
+
+        is_base_case = self._base_cases(node1, rule1, node2, rule2)
+        if is_base_case:
+            return bool(is_base_case + 1)
+
+        self._index += 1
+        self._ancestors1[node1], self._ancestors2[node2] = (self._index, self._index)
+
+        # Check all matches of children, if any are valid then trees are isomorphic
+        n = len(rule1.children)
+        child_order: List[int] = [0] * n
+        stack = [(0, i, {i}) for i in range(n - 1, -1, -1)]
+        while stack:
+            i1, i2, in_use = stack.pop()
+            if not self._are_isomorphic(rule1.children[i1], rule2.children[i2]):
+                continue
+            child_order[i1] = i2
+            if i1 == n - 1:
+                self._order_map[(node1, node2)] = child_order
+                self._cleanup(node1, node2)
+                return True
+            Isomorphism._extend_stack(i1, n, in_use, stack)
+
+        self._cleanup(node1, node2)
+        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]
+
+        if rule1.is_equivalence():
+            node1 = rule1.children[0]
+
+        return node1, node2
+
+    def _base_cases(
+        self,
+        node1: CombinatorialClass,
+        rule1: AbstractRule,
+        node2: CombinatorialClass,
+        rule2: AbstractRule,
+    ) -> int:
+        # If different number of children
+        if len(rule1.children) != len(rule2.children):
+            return Isomorphism._INVALID
+
+        # If leaf that has no equal-ancestors
+        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
+            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:
+            return Isomorphism._INVALID
+
+        # If they have occured before and that occurrence matches in both trees
+        if ind1 == ind2 != -1:
+            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]
+
+    @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})))
+
+
+class Node:
+    """Parse tree node."""
+
+    def __init__(
+        self,
+        rule: AbstractRule,
+        obj: CombinatorialObject,
+        get_rule: Callable[[CombinatorialClass], AbstractRule],
+    ):
+        self.rule = rule
+        self.obj = obj
+        if not rule.children:
+            self.children: Tuple[Optional["Node"], ...] = 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()
+            )
+
+    def build_obj(
+        self,
+        rule: AbstractRule,
+        get_order: Callable[[CombinatorialClass, CombinatorialClass], List[int]],
+        get_rule: Callable[[CombinatorialClass], AbstractRule],
+    ) -> CombinatorialObject:
+        """Parse tree's recursive build object function."""
+        if not self.children:
+            # TODO: stop special casing verification rules!
+            return self.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)
+        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)
+        )
+        assert isinstance(rule, Rule)
+        val: CombinatorialObject = next(rule.backward_map(child_objs))
+        return val
+
+
+class ParseTree:
+    """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 build_obj(
+        self,
+        root: AbstractRule,
+        get_order: Callable[[CombinatorialClass, CombinatorialClass], List[int]],
+        get_rule: Callable[[CombinatorialClass], AbstractRule],
+    ) -> CombinatorialObject:
+        """Build object from the other specification."""
+        return self.root.build_obj(root, get_order, get_rule)
+
+
+class Bijection:
+    """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:
+            return None
+
+    def __init__(
+        self,
+        spec: "CombinatorialSpecification",
+        other: "CombinatorialSpecification",
+        eq: AtomEquals = _default_equals,
+    ):
+        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)
+        return parse_tree.build_obj(
+            self.other.root_rule, self.iso.get_order, self.other.get_rule
+        )
diff --git a/comb_spec_searcher/rule_db/base.py b/comb_spec_searcher/rule_db/base.py
index 65ccf8fc4..b7b146ca2 100644
--- a/comb_spec_searcher/rule_db/base.py
+++ b/comb_spec_searcher/rule_db/base.py
@@ -21,7 +21,7 @@
     random_proof_tree,
     smallish_random_proof_tree,
 )
-from comb_spec_searcher.typing import RuleKey, SpecificationLabelsAndStrats
+from comb_spec_searcher.typing import RuleKey
 
 __all__ = ["RuleDBBase", "RuleDB"]
 
@@ -37,13 +37,21 @@ def __init__(self) -> None:
     def rule_to_strategy(self) -> MutableMapping[RuleKey, AbstractStrategy]:
         pass
 
+    @property
+    @abc.abstractmethod
+    def eqv_rule_to_strategy(self) -> MutableMapping[RuleKey, AbstractStrategy]:
+        pass
+
     def __eq__(self, other: object) -> bool:
         """Check if all stored information is the same."""
         if not isinstance(other, self.__class__):
             return NotImplemented
-        return bool(self.rule_to_strategy == other.rule_to_strategy)
+        return (
+            self.rule_to_strategy == other.rule_to_strategy
+            and self.eqv_rule_to_strategy == other.eqv_rule_to_strategy
+        )
 
-    def add(self, start: int, ends: Iterable[int], rule: AbstractRule) -> None:
+    def add(self, start: int, ends: Tuple[int, ...], rule: AbstractRule) -> None:
         """
         Add a rule to the database.
 
@@ -54,13 +62,16 @@ def add(self, start: int, ends: Iterable[int], rule: AbstractRule) -> None:
         ends = tuple(sorted(ends))
         if isinstance(rule, VerificationRule):
             self.set_verified(start)
-        is_equiv = len(ends) == 1 and rule.strategy.can_be_equivalent()
-        if is_equiv:
-            self.set_equivalent(start, ends[0])
-        if len(ends) != 1 or is_equiv or not self.are_equivalent(start, ends[0]):
-            # to avoid overwriting an equivalence rule with a non-equivalence
-            # rule, we only save if an equivalence rule, or does not have the
-            # same start -> ends as some equivalence rule.
+        if len(ends) == 1:
+            if rule.is_two_way():
+                self.equivdb.add_two_way_edge(start, ends[0])
+                self.eqv_rule_to_strategy[(start, ends)] = rule.strategy
+                self.rule_to_strategy.pop((start, ends), None)
+                self.rule_to_strategy.pop((ends[0], (start,)), None)
+            else:
+                self.equivdb.add_one_way_edge(start, ends[0])
+                self.rule_to_strategy[(start, ends)] = rule.strategy
+        else:
             self.rule_to_strategy[(start, ends)] = rule.strategy
 
     def is_verified(self, label: int) -> bool:
@@ -77,12 +88,15 @@ def are_equivalent(self, label: int, other: int) -> bool:
 
     def set_equivalent(self, label: int, other: int) -> None:
         """Mark label and other as equivalent."""
-        self.equivdb.union(label, other)
+        self.equivdb.add_two_way_edge(label, other)
 
     def rules_up_to_equivalence(self) -> Dict[int, Set[Tuple[int, ...]]]:
         """Return a defaultdict containing all rules up to the equivalence."""
+        self.equivdb.connect_cycles()
         rules_dict: Dict[int, Set[Tuple[int, ...]]] = defaultdict(set)
         for start, ends in self:
+            if len(ends) == 1 and self.are_equivalent(start, ends[0]):
+                continue
             rules_dict[self.equivdb[start]].add(
                 tuple(sorted(self.equivdb[e] for e in ends))
             )
@@ -96,19 +110,22 @@ def all_rules(self) -> Iterator[Tuple[int, Tuple[int, ...], AbstractStrategy]]:
     def __iter__(self) -> Iterator[Tuple[int, Tuple[int, ...]]]:
         """Iterate through rules as the pairs (start, end)."""
         for start, ends in self.rule_to_strategy:
-            if len(ends) != 1 or not self.are_equivalent(start, ends[0]):
-                yield start, ends
+            yield start, ends
 
     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
+        return (start, ends) in self.rule_to_strategy or (
+            start,
+            ends,
+        ) in self.eqv_rule_to_strategy
 
     def status(self, elaborate: bool) -> str:
         """Return a string describing the status of the rule database."""
         status = "RuleDB status:\n"
         table: List[Tuple[str, str]] = []
         table.append(("Combinatorial rules", f"{len(self.rule_to_strategy):,d}"))
+        table.append(("Equivalence rules", f"{len(self.eqv_rule_to_strategy):,d}"))
         if elaborate:
             eqv_labels = set()
             strategy_verified_labels = set()
@@ -149,6 +166,12 @@ def status(self, elaborate: bool) -> str:
         status += tabulate.tabulate(
             table, headers=("", "Total number"), colalign=("left", "right")
         ).replace("\n", "\n    ")
+        time_taken = round(self.equivdb.func_times["find_paths"], 2)
+        status += (
+            f"\n\tCalled find equiv path {self.equivdb.func_calls['find_paths']}"
+            + f" times, for total time of {time_taken}"
+            + " seconds.\n"
+        )
         return status
 
     ################################################################
@@ -228,62 +251,7 @@ def find_specification(
             raise SpecificationNotFound("No specification for label {}".format(label))
         return specification
 
-    def get_specification_rules(
-        self, label: int, minimization_time_limit: float, iterative: bool = False
-    ) -> SpecificationLabelsAndStrats:
-        """
-        Return a list of pairs (label, rule) which form a specification.
-        The specification returned is random, so two calls to the function
-        may result in a a different output.
-        """
-        proof_tree_node = self.find_specification(
-            label=label,
-            iterative=iterative,
-            minimization_time_limit=minimization_time_limit,
-        )
-        return self._get_specification_rules(label, proof_tree_node)
-
-    def _get_specification_rules(
-        self, label: int, proof_tree_node: Node
-    ) -> SpecificationLabelsAndStrats:
-        children: Dict[int, Tuple[int, ...]] = dict()
-        internal_nodes = set([label])
-        logger.info("Computing rule <-> equivalence rule mapping.")
-        eqv_rules = [
-            (node.label, tuple(sorted(child.label for child in node.children)))
-            for node in proof_tree_node.nodes()
-        ]
-        rule_from_equivalence_rules = self.rule_from_equivalence_rule_dict(eqv_rules)
-        logger.info("Finding actual rules.")
-        for eqv_start, eqv_ends in eqv_rules:
-            rule = rule_from_equivalence_rules.get((eqv_start, eqv_ends))
-            if rule is not None:
-                start, ends = rule
-                children[start] = ends
-                internal_nodes.update(ends)
-        res = []
-        eqv_paths = []
-        logger.info("Computing strategies used to create rules.")
-        for start, ends in children.items():
-            for eqv_label in internal_nodes:
-                if self.are_equivalent(start, eqv_label):
-                    path = self.equivdb.find_path(eqv_label, start)
-                    for a, b in zip(path[:-1], path[1:]):
-                        try:
-                            strategy = self.rule_to_strategy[(a, (b,))]
-                            res.append((a, strategy))
-                        except KeyError:
-                            strategy = self.rule_to_strategy[(b, (a,))]
-                            res.append((b, strategy))
-                    if len(path) > 1:
-                        eqv_paths.append(path)
-            strategy = self.rule_to_strategy[(start, ends)]
-            res.append((start, strategy))
-        return res, eqv_paths
-
-    def all_specifications(
-        self, label: int, iterative: bool = False
-    ) -> Iterator[SpecificationLabelsAndStrats]:
+    def all_specifications(self, label: int, iterative: bool = False) -> Iterator[Node]:
         """
         A generator that yields all specifications in the universe for
         the given label.
@@ -297,13 +265,9 @@ def all_specifications(
         prune(rules_dict)
 
         if self.equivdb[label] in rules_dict:
-            proof_trees = proof_tree_generator_dfs(rules_dict, root=self.equivdb[label])
-        for proof_tree_node in proof_trees:
-            yield self._get_specification_rules(label, proof_tree_node)
+            yield from proof_tree_generator_dfs(rules_dict, root=self.equivdb[label])
 
-    def get_smallest_specification(
-        self, label: int, iterative: bool = False
-    ) -> SpecificationLabelsAndStrats:
+    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.
@@ -338,14 +302,23 @@ def get_smallest_specification(
             except StopIteration:
                 minimum = middle + 1
         logger.info("The smallest specification is of size %s.", len(tree))
-        return self._get_specification_rules(label, tree)
+        return tree
 
 
 class RuleDB(RuleDBBase):
     def __init__(self) -> None:
         super().__init__()
         self._rule_to_strategy: Dict[Tuple[int, Tuple[int, ...]], AbstractStrategy] = {}
+        self._eqv_rule_to_strategy: Dict[
+            Tuple[int, Tuple[int, ...]], AbstractStrategy
+        ] = {}
 
     @property
     def rule_to_strategy(self) -> Dict[Tuple[int, Tuple[int, ...]], AbstractStrategy]:
         return self._rule_to_strategy
+
+    @property
+    def eqv_rule_to_strategy(
+        self,
+    ) -> Dict[Tuple[int, Tuple[int, ...]], AbstractStrategy]:
+        return self._eqv_rule_to_strategy
diff --git a/comb_spec_searcher/rule_db/forget.py b/comb_spec_searcher/rule_db/forget.py
index 8076f4bf9..e452f65c3 100644
--- a/comb_spec_searcher/rule_db/forget.py
+++ b/comb_spec_searcher/rule_db/forget.py
@@ -22,17 +22,24 @@ class RecomputingDict(MutableMapping[RuleKey, AbstractStrategy]):
     A mapping from rules to strategies that recompute the strategy every time it's
     needed instead of storing it in order to save memory.
 
+    If only equiv is set to true the returned strategy will only create two way rules.
+
     Also in order to save memory we store flat version of the rules, i.e. for a rule
     (a, (b, c,...)) we store (a, b, c,...)
     """
 
     def __init__(
-        self, classdb: ClassDB, strat_pack: StrategyPack, equivdb: EquivalenceDB
+        self,
+        classdb: ClassDB,
+        strat_pack: StrategyPack,
+        equivdb: EquivalenceDB,
+        only_equiv: bool,
     ) -> None:
         self.classdb = classdb
         self.equivdb = equivdb
         self.pack = strat_pack
         self.rules: Set[Tuple[int, ...]] = set()
+        self.only_equiv: bool = only_equiv
 
     @staticmethod
     def _flatten(tuple_: RuleKey) -> Tuple[int, ...]:
@@ -42,18 +49,9 @@ def _flatten(tuple_: RuleKey) -> Tuple[int, ...]:
     def _unflatten(tuple_: Tuple[int, ...]) -> RuleKey:
         return (tuple_[0], tuple_[1:])
 
-    def _is_equiv(self, key: RuleKey) -> bool:
-        """
-        Returns True if the key should be the key of an equivalence rule.
-        """
-        start = key[0]
-        ends = key[1]
-        return len(ends) == 1 and self.equivdb.equivalent(start, ends[0])
-
     def __getitem__(self, key: RuleKey) -> AbstractStrategy:
         if self._flatten(key) not in self.rules:
             raise KeyError(key)
-        expect_equiv = self._is_equiv(key)
         possible_labels = (key[0],) + key[1]
         for label, strat in itertools.product(possible_labels, self.pack):
             comb_class = self.classdb.get_class(label)
@@ -80,7 +78,7 @@ def __getitem__(self, key: RuleKey) -> AbstractStrategy:
                         sorted(map(self.classdb.get_label, nonempty_children))
                     )
                     if (start_label, end_labels) == key:
-                        if expect_equiv and not rule.strategy.can_be_equivalent():
+                        if self.only_equiv and not rule.is_two_way():
                             continue
                         return rule.strategy
                 except StrategyDoesNotApply:
@@ -95,6 +93,7 @@ def __getitem__(self, key: RuleKey) -> AbstractStrategy:
         raise RuntimeError(err_message)
 
     def __setitem__(self, key: RuleKey, value: AbstractStrategy) -> None:
+        assert not self.only_equiv or len(key[1]) == 1
         self.rules.add(self._flatten(key))
 
     def __delitem__(self, key: RuleKey) -> None:
@@ -114,8 +113,17 @@ def __contains__(self, key: object) -> bool:
 class RuleDBForgetStrategy(RuleDBBase):
     def __init__(self, classdb: ClassDB, strat_pack: StrategyPack) -> None:
         super().__init__()
-        self.rules = RecomputingDict(classdb, strat_pack, self.equivdb)
+        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
+        )
+
+    @property
+    def rule_to_strategy(self) -> RecomputingDict:
+        return self._rule_to_strategy
 
     @property
-    def rule_to_strategy(self):
-        return self.rules
+    def eqv_rule_to_strategy(self) -> RecomputingDict:
+        return self._eqv_rule_to_strategy
diff --git a/comb_spec_searcher/specification.py b/comb_spec_searcher/specification.py
index cfdb72f28..ea6a0b6dd 100644
--- a/comb_spec_searcher/specification.py
+++ b/comb_spec_searcher/specification.py
@@ -5,7 +5,7 @@
 from copy import copy
 from functools import reduce
 from operator import mul
-from typing import Dict, Generic, Iterable, Iterator, List, Sequence, Set, Tuple, Union
+from typing import Dict, Generic, Iterable, Iterator, List, Optional, Set, Union
 
 import sympy
 from logzero import logger
@@ -24,14 +24,12 @@
     InvalidOperationError,
     TaylorExpansionError,
 )
+from .isomorphism import AtomEquals, Bijection, Isomorphism
 from .specification_drawer import SpecificationDrawer
 from .strategies import (
-    AbstractStrategy,
     EmptyStrategy,
     EquivalencePathRule,
-    ReverseRule,
     Rule,
-    Strategy,
     StrategyPack,
     VerificationRule,
     VerificationStrategy,
@@ -59,20 +57,13 @@ class CombinatorialSpecification(
     def __init__(
         self,
         root: CombinatorialClassType,
-        strategies: Iterable[
-            Tuple[
-                CombinatorialClassType,
-                AbstractStrategy[CombinatorialClassType, CombinatorialObjectType],
-            ]
-        ],
-        equivalence_paths: Iterable[Sequence[CombinatorialClassType]],
+        rules: Iterable[AbstractRule[CombinatorialClassType, CombinatorialObjectType]],
     ):
         self.root = root
-        self.rules_dict: Dict[CombinatorialClassType, AbstractRule] = {}
-        self._populate_rules_dict(strategies, equivalence_paths)
-        self._remove_redundant_rules()
+        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()
         self._set_subrules()
 
     def _set_subrules(self) -> None:
@@ -82,35 +73,57 @@ def _set_subrules(self) -> None:
         ):  # list as we lazily assign empty rules
             rule.set_subrecs(self.get_rule)
 
-    def _populate_rules_dict(
-        self,
-        strategies: Iterable[
-            Tuple[
-                CombinatorialClassType,
-                AbstractStrategy[CombinatorialClassType, CombinatorialObjectType],
-            ]
-        ],
-        equivalence_paths: Iterable[Sequence[CombinatorialClassType]],
-    ) -> None:
-        logger.info("Creating rules.")
-        equivalence_rules: Dict[
-            Tuple[CombinatorialClassType, CombinatorialClassType], Rule
-        ] = {}
-        for comb_class, strategy in strategies:
-            if isinstance(strategy, AlreadyVerified):
-                continue
-            rule = strategy(comb_class)
-            non_empty_children = rule.non_empty_children()
+    def _group_equiv_in_path(self) -> None:
+        """
+        Group chain of equivalence rules of the specification into equivalence paths.
+        """
+        self._ungroup_equiv_path()
+        not_hidden_classes = {self.root}
+        for rule in self.rules_dict.values():
             if rule.is_equivalence():
+                continue
+            not_hidden_classes.add(rule.comb_class)
+            not_hidden_classes.update(rule.children)
+        # Creating paths
+        comb_class_stack = [self.root]
+        visited: Set[CombinatorialClassType] = set()
+        path_rules: List[Rule] = []
+        eqv_path_rules: Dict[CombinatorialClassType, EquivalencePathRule] = {}
+        while comb_class_stack:
+            if path_rules and path_rules[-1].children[0] in not_hidden_classes:
+                new_rule: EquivalencePathRule = EquivalencePathRule(path_rules)
+                eqv_path_rules[new_rule.comb_class] = new_rule
+                path_rules.clear()
+            comb_class = comb_class_stack.pop()
+            if comb_class in not_hidden_classes and comb_class in visited:
+                continue
+            visited.add(comb_class)
+            rule = self.get_rule(comb_class)
+            comb_class_stack.extend(rule.children)
+            if rule.is_equivalence() and not isinstance(rule, EquivalencePathRule):
                 assert isinstance(rule, Rule)
-                equivalence_rules[(comb_class, non_empty_children[0])] = (
-                    rule if len(rule.children) == 1 else rule.to_equivalence_rule()
-                )
-            elif non_empty_children or isinstance(rule, VerificationRule):
-                self.rules_dict[comb_class] = rule
-            else:
-                raise ValueError("Non verification rule has no children.")
-        self._add_equivalence_path_rules(equivalence_paths, equivalence_rules)
+                assert not path_rules or path_rules[-1].children[0] == rule.comb_class
+                path_rules.append(rule)
+        # Clearing the equivalences and replacing by paths
+        self.rules_dict = {
+            comb_class: rule
+            for comb_class, rule in self.rules_dict.items()
+            if comb_class in not_hidden_classes
+        }
+        self.rules_dict.update(eqv_path_rules)
+        assert self._is_valid_spec()
+
+    def _ungroup_equiv_path(self) -> None:
+        """
+        Remove the equiv path and replacing only with normal rules.
+        """
+        new_rules: Dict[CombinatorialClassType, Rule] = {}
+        for rule in self.rules_dict.values():
+            if not isinstance(rule, EquivalencePathRule):
+                continue
+            for r in rule.rules:
+                new_rules[r.comb_class] = r
+        self.rules_dict.update(new_rules)
 
     def expand_verified(self) -> None:
         """
@@ -125,95 +138,44 @@ def expand_verified(self) -> None:
                 except InvalidOperationError:
                     logger.info("Can't expand the rule:\n%s", rule)
         if verification_packs:
-            for comb_class in verification_packs:
-                self.rules_dict.pop(comb_class)
-            self._expand_verified_comb_classes(verification_packs)
+            for comb_class, pack in verification_packs.items():
+                self.expand_comb_class(comb_class, pack)
             self.expand_verified()
+            self._group_equiv_in_path()
             self._set_subrules()
 
-    def expand_comb_class(self, comb_class: Union[int, CombinatorialClassType]) -> None:
+    def expand_comb_class(
+        self, comb_class: Union[int, CombinatorialClassType], pack: StrategyPack
+    ) -> None:
         """
         Will try to expand a particular class with respect to the strategy pack
         that the VerificationStrategy has.
         """
-        rule = self.get_rule(comb_class)
-        if isinstance(rule, VerificationRule):
-            try:
-                pack = rule.pack()
-            except InvalidOperationError:
-                logger.info("Can't expand the rule:\n%s", rule)
-                return
-            self.rules_dict.pop(rule.comb_class)
-            self._expand_verified_comb_classes({rule.comb_class: pack})
-        else:
-            logger.info("Can't expand the rule:\n%s", rule)
-
-    def _add_equivalence_path_rules(
-        self,
-        equivalence_paths: Iterable[Sequence[CombinatorialClassType]],
-        equivalence_rules: Dict[
-            Tuple[CombinatorialClassType, CombinatorialClassType], Rule
-        ],
-    ) -> None:
-        logger.info("Creating equivalence path rules.")
-        for eqv_path in equivalence_paths:
-            if len(eqv_path) > 1:
-                start = eqv_path[0]
-                rules = []
-                for a, b in zip(eqv_path[:-1], eqv_path[1:]):
-                    try:
-                        rule = equivalence_rules[(a, b)]
-                    except KeyError:
-                        rule = equivalence_rules[(b, a)].to_reverse_rule()
-                    rules.append(rule)
-                self.rules_dict[start] = EquivalencePathRule(rules)
-
-    def _remove_redundant_rules(self) -> None:
-        """
-        Any redundant rules passed to the __init__ method are removed using
-        this method.
-        """
-        rules_dict = copy(self.rules_dict)
-
-        def prune(comb_class: CombinatorialClassType) -> None:
-            try:
-                rule = rules_dict.pop(comb_class)
-            except KeyError:
-                assert comb_class in self.rules_dict or comb_class.is_empty()
-                return
-            if isinstance(rule, EquivalencePathRule):
-                try:
-                    rule = rules_dict.pop(rule.children[0])
-                except KeyError:
-                    assert comb_class in self.rules_dict
-                    return
-            for child in rule.children:
-                prune(child)
+        # pylint: disable=import-outside-toplevel
+        from .comb_spec_searcher import CombinatorialSpecificationSearcher
 
-        prune(self.root)
-
-        logger.info("Removed %s redundant rules.", len(rules_dict.values()))
-        for rule in rules_dict.values():
-            self.rules_dict.pop(rule.comb_class)
+        rule = self.get_rule(comb_class)
+        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
 
-    def _expand_verified_comb_classes(
-        self, verification_packs: Dict[CombinatorialClassType, StrategyPack]
-    ) -> None:
-        for comb_class, pack in verification_packs.items():
-            # pylint: disable=import-outside-toplevel
-            from .comb_spec_searcher import CombinatorialSpecificationSearcher
-
-            css = CombinatorialSpecificationSearcher(
-                comb_class,
-                pack.add_verification(
-                    AlreadyVerified(self.rules_dict), apply_first=True
-                ),
-            )
-            logger.info(css.run_information())
-            # pylint: disable=protected-access
-            comb_class_rules, comb_class_eqv_paths = css._auto_search_rules()
-            self._populate_rules_dict(comb_class_rules, comb_class_eqv_paths)
-        self._remove_redundant_rules()
+    def _is_valid_spec(self) -> bool:
+        """Checks that each class is on a left hand side."""
+        comb_classes = set()
+        for rule in self.rules_dict.values():
+            comb_classes.add(rule.comb_class)
+            comb_classes.update(rule.children)
+        return self.root in comb_classes and all(
+            c in self.rules_dict or c.is_empty() for c in comb_classes
+        )
 
     def get_rule(
         self, comb_class: Union[int, CombinatorialClassType]
@@ -229,7 +191,9 @@ def get_rule(
                 ) from e
         if comb_class not in self.rules_dict:
             assert comb_class.is_empty(), "rule not in the spec and not empty"
-            empty_strat = EmptyStrategy()
+            empty_strat = EmptyStrategy[
+                CombinatorialClassType, CombinatorialObjectType
+            ]()
             self.rules_dict[comb_class] = empty_strat(comb_class)
         return self.rules_dict[comb_class]
 
@@ -434,11 +398,27 @@ def sanity_check(self, length: int = 5) -> bool:
 
         Raise an SanityCheckFailure error if it fails.
         """
-        return all(
+        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)
         )
 
+    def get_bijection_to(
+        self, other: "CombinatorialSpecification", eq: Optional[AtomEquals] = None
+    ) -> Optional[Bijection]:
+        """Get bijection from self to other."""
+        if eq is None:
+            return Bijection.construct(self, other)
+        return Bijection.construct(self, other, eq)
+
+    def are_isomorphic(
+        self, other: "CombinatorialSpecification", eq: Optional[AtomEquals] = None
+    ) -> bool:
+        """Check if self is isomorphic to other."""
+        if eq is None:
+            return Isomorphism.check(self, other)
+        return Isomorphism.check(self, other, eq)
+
     def show(
         self, levels_shown: int = 0, levels_expand: int = 0, verbose: bool = False
     ) -> None:
@@ -543,30 +523,9 @@ def equations_string(self) -> str:
 
     def to_jsonable(self) -> dict:
         """Return a JSON serializable dictionary for the specification."""
-        rules = [
-            rule
-            for rule in self.rules_dict.values()
-            if not isinstance(rule, EquivalencePathRule)
-        ]
-        eqv_paths = []
-        for rule in self.rules_dict.values():
-            if isinstance(rule, EquivalencePathRule):
-                eqv_path = []
-                for c, r in rule.eqv_path_rules():
-                    rules.append(r)
-                    eqv_path.append(c.to_jsonable())
-                eqv_path.append(rule.children[0].to_jsonable())
-                eqv_paths.append(eqv_path)
-        strategies = [
-            (rule.children[0].to_jsonable(), rule.strategy.to_jsonable())
-            if isinstance(rule, ReverseRule)
-            else (rule.comb_class.to_jsonable(), rule.strategy.to_jsonable())
-            for rule in rules
-        ]
         return {
             "root": self.root.to_jsonable(),
-            "strategies": strategies,
-            "eqv_paths": eqv_paths,
+            "rules": [rule.to_jsonable() for rule in self.rules_dict.values()],
         }
 
     @classmethod
@@ -575,20 +534,9 @@ def from_dict(cls, d):
         Return the specification with the dictionary outputter by the
         'to_jsonable' method
         """
-        strategies = []
-        for comb_class_dict, strategy_dict in d["strategies"]:
-            comb_class = CombinatorialClass.from_dict(comb_class_dict)
-            strategy = Strategy.from_dict(strategy_dict)
-            assert isinstance(strategy, (Strategy, VerificationStrategy))
-            strategies.append((comb_class, strategy))
-        return cls(
-            root=CombinatorialClass.from_dict(d["root"]),
-            strategies=strategies,
-            equivalence_paths=[
-                [CombinatorialClass.from_dict(class_dict) for class_dict in eqv_path]
-                for eqv_path in d["eqv_paths"]
-            ],
-        )
+        root = CombinatorialClass.from_dict(d.pop("root"))
+        rules = [AbstractRule.from_dict(rule_dict) for rule_dict in d.pop("rules")]
+        return CombinatorialSpecification(root, rules)
 
 
 class AlreadyVerified(VerificationStrategy[CombinatorialClass, CombinatorialObject]):
diff --git a/comb_spec_searcher/specification_extrator.py b/comb_spec_searcher/specification_extrator.py
new file mode 100644
index 000000000..bb28bdbed
--- /dev/null
+++ b/comb_spec_searcher/specification_extrator.py
@@ -0,0 +1,117 @@
+import itertools
+from typing import Dict, Iterator, List, 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
+
+
+class SpecificationRuleExtractor:
+    def __init__(
+        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.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] = {}
+        self._populate_decompositions()
+        self._populate_equivalences()
+        self._check()
+
+    def _populate_decompositions(self) -> None:
+        """
+        Populate the rules dict with the labels of decomposition rules.
+        """
+        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
+            self.rules_dict[parent] = children
+            self.eqvparent_to_parent[eqvrule[0]] = parent
+
+    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:
+            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])
+                self.rules_dict[parent] = (child,)
+
+    def _check(self):
+        """
+        Perform a check to make sure that the spec make sense.
+        """
+        all_rhs = set(
+            itertools.chain.from_iterable(
+                children for children in self.rules_dict.values()
+            )
+        )
+        assert all_rhs.issubset(self.rules_dict), "Something is not on the lhs"
+
+    def _find_rule(self, parent: int, children: Tuple[int, ...]) -> AbstractRule:
+        try:
+            strategy = self.ruledb.rule_to_strategy[(parent, children)]
+        except KeyError as e:
+            if len(children) != 1:
+                raise self._missing_rule_error(parent, children) from e
+        else:
+            return strategy(self.classdb.get_class(parent))
+        # From now on we are looking for a two way strategy
+        assert len(children) == 1
+        try:
+            strategy = self.ruledb.eqv_rule_to_strategy[(parent, children)]
+        except KeyError:
+            pass
+        else:
+            rule = strategy(self.classdb.get_class(parent))
+            assert isinstance(rule, Rule)
+            return rule if len(rule.children) == 1 else rule.to_equivalence_rule()
+        # Now trying a reverse rule
+        try:
+            strategy = self.ruledb.eqv_rule_to_strategy[(children[0], (parent,))]
+        except KeyError:
+            pass
+        else:
+            rule = strategy(self.classdb.get_class(children[0]))
+            assert isinstance(rule, Rule)
+            rule = rule if len(rule.children) == 1 else rule.to_equivalence_rule()
+            return rule.to_reverse_rule(0)
+        raise self._missing_rule_error(parent, children)
+
+    def _missing_rule_error(self, parent: int, children: Tuple[int, ...]) -> ValueError:
+        """
+        Return a ValueError saying the rule was not found.
+        """
+        msg = f"Unable to retrieve rule for {(parent, children)}\n"
+        msg += f"Class {parent} is:\n{self.classdb.get_class(parent)}\n"
+        for child in children:
+            msg += f"Class {child} is:\n{self.classdb.get_class(child)}\n"
+        return ValueError(msg)
+
+    def rules(self) -> Iterator[AbstractRule]:
+        for parent, children in self.rules_dict.items():
+            yield self._find_rule(parent, children)
diff --git a/comb_spec_searcher/strategies/__init__.py b/comb_spec_searcher/strategies/__init__.py
index 965b4e182..2bba8ad17 100644
--- a/comb_spec_searcher/strategies/__init__.py
+++ b/comb_spec_searcher/strategies/__init__.py
@@ -1,4 +1,10 @@
-from .constructor import CartesianProduct, Constructor, DisjointUnion
+from .constructor import (
+    CartesianProduct,
+    Complement,
+    Constructor,
+    DisjointUnion,
+    Quotient,
+)
 from .rule import EquivalencePathRule, ReverseRule, Rule, VerificationRule
 from .strategy import (
     AbstractStrategy,
@@ -15,8 +21,10 @@
 
 __all__ = [
     "CartesianProduct",
+    "Quotient",
     "Constructor",
     "DisjointUnion",
+    "Complement",
     "EquivalencePathRule",
     "ReverseRule",
     "Rule",
diff --git a/comb_spec_searcher/strategies/constructor/__init__.py b/comb_spec_searcher/strategies/constructor/__init__.py
new file mode 100644
index 000000000..1152fbf5e
--- /dev/null
+++ b/comb_spec_searcher/strategies/constructor/__init__.py
@@ -0,0 +1,11 @@
+from .base import Constructor
+from .cartesian import CartesianProduct, Quotient
+from .disjoint import Complement, DisjointUnion
+
+__all__ = [
+    "Constructor",
+    "DisjointUnion",
+    "Complement",
+    "CartesianProduct",
+    "Quotient",
+]
diff --git a/comb_spec_searcher/strategies/constructor/base.py b/comb_spec_searcher/strategies/constructor/base.py
new file mode 100644
index 000000000..2469ce6a0
--- /dev/null
+++ b/comb_spec_searcher/strategies/constructor/base.py
@@ -0,0 +1,89 @@
+"""
+The constructor class contains all the method, and logic, needed to get the
+enumeration, generate objects, and sample objects.
+"""
+import abc
+from typing import Callable, Generic, Iterator, List, Optional, Tuple
+
+import sympy
+
+from comb_spec_searcher.typing import (
+    CombinatorialClassType,
+    CombinatorialObjectType,
+    Parameters,
+    ParametersMap,
+    RelianceProfile,
+    SubObjects,
+    SubRecs,
+    SubSamplers,
+    SubTerms,
+    Terms,
+)
+
+
+class Constructor(abc.ABC, Generic[CombinatorialClassType, CombinatorialObjectType]):
+    """The constructor is akin to the 'counting function' in the comb exp paper."""
+
+    @abc.abstractmethod
+    def get_equation(
+        self, lhs_func: sympy.Function, rhs_funcs: Tuple[sympy.Function, ...]
+    ) -> sympy.Eq:
+        """
+        Return the sympy.Eq in the form lhs_func = f(rhs_funcs).
+        """
+
+    @abc.abstractmethod
+    def reliance_profile(self, n: int, **parameters: int) -> RelianceProfile:
+        """
+        Return the reliance profile. That is for the parameters given,
+        which parameters of each individual subclass are required.
+        """
+
+    @abc.abstractmethod
+    def get_terms(
+        self, parent_terms: Callable[[int], Terms], subterms: SubTerms, n: int
+    ) -> Terms:
+        """
+        Return the terms for n given the subterms of the children.
+        """
+
+    @abc.abstractmethod
+    def get_sub_objects(
+        self, subobjs: SubObjects, n: int
+    ) -> Iterator[
+        Tuple[Parameters, Tuple[List[Optional[CombinatorialObjectType]], ...]]
+    ]:
+        """
+        Iterate over all the possibe subobjs/image of the bijection implied by
+        the constructor together with the parameters.
+        """
+
+    @abc.abstractmethod
+    def random_sample_sub_objects(
+        self,
+        parent_count: int,
+        subsamplers: SubSamplers,
+        subrecs: SubRecs,
+        n: int,
+        **parameters: int,
+    ) -> Tuple[Optional[CombinatorialObjectType], ...]:
+        """Return a randomly sampled subobjs/image of the bijection implied
+        by the constructor."""
+
+    @staticmethod
+    def build_param_map(
+        child_pos_to_parent_pos: Tuple[Tuple[int, ...], ...], num_parent_params: int
+    ) -> ParametersMap:
+        """
+        Return a parameters map according to the given pos map.
+        """
+
+        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)
+
+        return param_map
diff --git a/comb_spec_searcher/strategies/constructor.py b/comb_spec_searcher/strategies/constructor/cartesian.py
similarity index 58%
rename from comb_spec_searcher/strategies/constructor.py
rename to comb_spec_searcher/strategies/constructor/cartesian.py
index da7216308..105634964 100644
--- a/comb_spec_searcher/strategies/constructor.py
+++ b/comb_spec_searcher/strategies/constructor/cartesian.py
@@ -1,34 +1,21 @@
-"""
-The constructor class contains all the method, and logic, needed to get the
-enumeration, generate objects, and sample objects.
-
-The default constructors implemented are:
-- CartesianProduct
-- DisjointUnion
-- Empty
-
-Currently the constructors are implemented in one variable, namely 'n' which is
-used throughout to denote size.
-"""
-import abc
+import random
 from collections import defaultdict
 from itertools import product
-from random import randint
 from typing import (
     Callable,
     Counter,
     Dict,
-    Generic,
     Iterable,
     Iterator,
     List,
     Optional,
     Tuple,
     TypeVar,
+    Union,
     cast,
 )
 
-from sympy import Eq, Function
+import sympy
 
 from comb_spec_searcher import utils
 from comb_spec_searcher.typing import (
@@ -44,57 +31,11 @@
     Terms,
 )
 
-__all__ = ("Constructor", "CartesianProduct", "DisjointUnion")
+from .base import Constructor
 
 T = TypeVar("T")
 
 
-class Constructor(abc.ABC, Generic[CombinatorialClassType, CombinatorialObjectType]):
-    """The constructor is akin to the 'counting function' in the comb exp paper."""
-
-    @abc.abstractmethod
-    def get_equation(self, lhs_func: Function, rhs_funcs: Tuple[Function, ...]) -> Eq:
-        """
-        Return the sympy.Eq in the form lhs_func = f(rhs_funcs).
-        """
-
-    @abc.abstractmethod
-    def reliance_profile(self, n: int, **parameters: int) -> RelianceProfile:
-        """
-        Return the reliance profile. That is for the parameters given,
-        which parameters of each individual subclass are required.
-        """
-
-    @abc.abstractmethod
-    def get_terms(self, subterms: SubTerms, n: int) -> Terms:
-        """
-        Return the terms for n given the subterms of the children.
-        """
-
-    @abc.abstractmethod
-    def get_sub_objects(
-        self, subobjs: SubObjects, n: int
-    ) -> Iterator[
-        Tuple[Parameters, Tuple[List[Optional[CombinatorialObjectType]], ...]]
-    ]:
-        """
-        Iterate over all the possibe subobjs/image of the bijection implied by
-        the constructor together with the parameters.
-        """
-
-    @abc.abstractmethod
-    def random_sample_sub_objects(
-        self,
-        parent_count: int,
-        subsamplers: SubSamplers,
-        subrecs: SubRecs,
-        n: int,
-        **parameters: int
-    ):
-        """Return a randomly sampled subobjs/image of the bijection implied
-        by the constructor."""
-
-
 class CartesianProduct(Constructor[CombinatorialClassType, CombinatorialObjectType]):
     """
     The CartesianProduct is initialised with the children of the rule that is
@@ -124,7 +65,9 @@ def __init__(
             self.extra_parameters = tuple(extra_parameters)
         else:
             self.extra_parameters = tuple(dict() for _ in children)
-        self._children_param_maps = self._build_children_param_map(parent, children)
+        self._children_param_maps = self._build_children_param_map(
+            parent, children, self.extra_parameters
+        )
 
         self.minimum_sizes = {"n": parent.minimum_size_of_object()}
         for k in parent.extra_parameters:
@@ -167,17 +110,18 @@ def max_sizes(self) -> Tuple[Optional[int], ...]:
         """
         return tuple(d.get("n", None) for d in self.max_child_sizes)
 
+    @staticmethod
     def _build_children_param_map(
-        self,
         parent: CombinatorialClassType,
         children: Tuple[CombinatorialClassType, ...],
+        extra_parameters: Tuple[Dict[str, str], ...],
     ) -> Tuple[ParametersMap, ...]:
         map_list: List[ParametersMap] = []
         num_parent_params = len(parent.extra_parameters)
         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 child, extra_param in zip(children, extra_parameters):
             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)
@@ -191,40 +135,22 @@ def _build_children_param_map(
                 for child_param in child.extra_parameters
             )
             map_list.append(
-                self._build_param_map(child_pos_to_parent_pos, num_parent_params)
+                Constructor.build_param_map(child_pos_to_parent_pos, num_parent_params)
             )
         return tuple(map_list)
 
-    @staticmethod
-    def _build_param_map(
-        child_pos_to_parent_pos: Tuple[Tuple[int, ...], ...], num_parent_params: int
-    ) -> ParametersMap:
-        """
-        Build the ParametersMap that will map according to the given child pos to parent
-        pos map given.
-        """
-
-        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)
-
-        return param_map
-
-    def get_equation(self, lhs_func: Function, rhs_funcs: Tuple[Function, ...]) -> Eq:
+    def get_equation(
+        self, lhs_func: sympy.Function, rhs_funcs: Tuple[sympy.Function, ...]
+    ) -> sympy.Eq:
         res = 1
         for extra_parameters, rhs_func in zip(self.extra_parameters, rhs_funcs):
             res *= rhs_func.subs(
                 {child: parent for parent, child in extra_parameters.items()},
                 simultaneous=True,
             )
-        return Eq(lhs_func, res)
+        return sympy.Eq(lhs_func, res)
 
     def reliance_profile(self, n: int, **parameters: int) -> RelianceProfile:
-        #  TODO: consider when parameters are subsets of each other etc
         assert all(
             set(["n", *parameters]) == set(min_child_sizes)
             for min_child_sizes in self.min_child_sizes
@@ -332,13 +258,15 @@ def get_extra_parameters(
             res.append(extra_params)
         return res
 
-    def get_terms(self, subterms: SubTerms, n: int) -> Terms:
+    def get_terms(
+        self, parent_terms: Callable[[int], Terms], subterms: SubTerms, n: int
+    ) -> Terms:
         new_terms: Terms = Counter()
         size_compositions = utils.compositions(
             n, len(subterms), self.min_sizes, self.max_sizes
         )
         for sizes in size_compositions:
-            for param_value_pairs in self._params_value_pairs_combinations(
+            for param_value_pairs in self.params_value_pairs_combinations(
                 sizes, subterms
             ):
                 new_param = self._new_param(*(p for p, _ in param_value_pairs))
@@ -356,11 +284,11 @@ def _new_param(self, *children_params: Parameters) -> Parameters:
         return tuple(sum(vals) for vals in zip(*mapped_params))
 
     @staticmethod
-    def _params_value_pairs_combinations(
+    def params_value_pairs_combinations(
         sizes: Tuple[int, ...],
         sub_getters: Tuple[Callable[[int], Dict[Parameters, T]], ...],
     ) -> Iterator[Tuple[Tuple[Parameters, T], ...]]:
-        """"""
+        assert len(sizes) == len(sub_getters), sub_getters
         children_values = (sg(s) for sg, s in zip(sub_getters, sizes))
         return product(*(c.items() for c in children_values))
 
@@ -373,7 +301,7 @@ def get_sub_objects(
             n, len(subobjs), self.min_sizes, self.max_sizes
         )
         for sizes in size_compositions:
-            for param_objs_pairs in self._params_value_pairs_combinations(
+            for param_objs_pairs in self.params_value_pairs_combinations(
                 sizes, subobjs
             ):
                 new_param = self._new_param(*(p for p, _ in param_objs_pairs))
@@ -389,9 +317,9 @@ def random_sample_sub_objects(
         subsamplers: SubSamplers,
         subrecs: SubRecs,
         n: int,
-        **parameters: int
-    ):
-        random_choice = randint(1, parent_count)
+        **parameters: int,
+    ) -> Tuple[Optional[CombinatorialObjectType], ...]:
+        random_choice = random.randint(1, parent_count)
         total = 0
         for child_parameters in self._valid_compositions(n, **parameters):
             tmp = 1
@@ -423,164 +351,241 @@ def __str__(self) -> str:
         return "Cartesian product"
 
 
-class DisjointUnion(Constructor[CombinatorialClassType, CombinatorialObjectType]):
+class Quotient(Constructor[CombinatorialClassType, CombinatorialObjectType]):
     """
-    The DisjointUnion constructor takes as input the children. Each constructor
-    is unique up to the length of the children being used to count.
+    This constructor implements the quotient operation. A rule A -> (B, C, D)
+    should be read as A = B / (C * D).
 
-    Extra parameters are passed on using the parameters dictionaries. Each
-    dictionary's keys should be the extra variable of the child pointing the
-    variable on the parent it came from.
-    If a parents variable does not map to a child, then this variable must be 0
-    as the child contains no occurences.
+    This comes from a cartesian product rule B = A * C * D
 
-    The fixed value dictionaries passed will be used ensure that the parameter
-    of a child must take on the given value.
+    The initialiser should be made using the original cartesian product rule.
+
+    The parent is B, idx is the index of the child the rule is being rearranged
+    for and the extra parameters maps from parent to children.
+
+    To get terms we compute a polynomial equation for a fixed size a_n(x,...) =
+    b_n(x,...) * c(x, ...) where b_n is the terms for the flipped child of length n and
+    then make the division.
     """
 
     def __init__(
         self,
         parent: CombinatorialClassType,
         children: Tuple[CombinatorialClassType, ...],
+        idx: int,
         extra_parameters: Optional[Tuple[Dict[str, str], ...]] = None,
-        fixed_values: Optional[Tuple[Dict[str, int], ...]] = None,
     ):
         self.number_of_children = len(children)
         if extra_parameters is not None:
-            self.extra_parameters = extra_parameters
-            assert len(extra_parameters) == len(children)
+            self.extra_parameters = tuple(extra_parameters)
         else:
-            assert not parent.extra_parameters
-            self.extra_parameters = tuple(
-                dict() for _ in range(self.number_of_children)
-            )
-        self._children_param_maps = self._build_children_param_maps(parent, children)
-        self.zeroes = tuple(
-            frozenset(parent.extra_parameters) - frozenset(parameter.keys())
-            for parameter in self.extra_parameters
+            self.extra_parameters = tuple(dict() for _ in children)
+        self.idx = idx
+        self._children_param_maps = CartesianProduct._build_children_param_map(
+            parent, children, self.extra_parameters
         )
-        if fixed_values is not None:
-            self.fixed_values = fixed_values
-            assert len(fixed_values) == len(children)
-        else:
-            self.fixed_values = tuple({} for _ in children)
-
-    def get_equation(self, lhs_func: Function, rhs_funcs: Tuple[Function, ...]) -> Eq:
-        res = 0
-        for rhs_func, extra_parameters in zip(rhs_funcs, self.extra_parameters):
-            res += rhs_func.subs(
-                {child: parent for parent, child in extra_parameters.items()},
-                simultaneous=True,
-            )
-        return Eq(lhs_func, res)
-
-    def reliance_profile(self, n: int, **parameters: int) -> RelianceProfile:
-        # TODO: implement in multiple variables and use in get_recurrence
-        assert not parameters, "only implemented in one variable, namely 'n'"
-        return tuple({"n": (n,)} for _ in range(self.number_of_children))
-
-    def get_extra_parameters(
-        self, n: int, **parameters: int
-    ) -> List[Optional[Dict[str, int]]]:
-        """
-        Will return the extra parameters dictionary based on the parent's
-        parameters. If there is a contradiction, that is some child parameter
-        is given two or more different values that do not match, then None
-        will be returned for that child, indicating that 0 objects on the
-        child match the parents parameters.
-        """
-        res: List[Optional[Dict[str, int]]] = []
-        for i, extra_parameters in enumerate(self.extra_parameters):
-            update_params: Dict[str, int] = {**self.fixed_values[i]}
-            for parent_var, child_var in extra_parameters.items():
-                updated_value = parameters[parent_var]
-                if child_var not in update_params:
-                    update_params[child_var] = updated_value
-                elif update_params[child_var] != updated_value:
-                    break
-            else:
-                res.append(update_params)
-                continue
-            res.append(None)
-        return res
-
-    def get_terms(self, subterms: SubTerms, n: int) -> Terms:
-        new_terms: Terms = Counter()
-        for child_terms, param_map in zip(subterms, self._children_param_maps):
-            for param, value in child_terms(n).items():
-                new_terms[param_map(param)] += value
-        return new_terms
+        self._min_sizes = tuple(child.minimum_size_of_object() for child in children)
+        self._max_sizes = tuple(
+            child.minimum_size_of_object() if child.is_atom() else None
+            for child in children
+        )
+        self._parent_param_map = self._build_parent_param_map(parent, children)
+        self._parent_shift = sum(self._min_sizes) - self._min_sizes[self.idx]
+        self._num_parent_params = len(parent.extra_parameters)
 
     @staticmethod
-    def _build_param_map(
-        child_pos_to_parent_pos: Tuple[Optional[Tuple[int, ...]], ...],
-        num_parent_params: int,
+    def build_param_map(
+        child_pos_to_parent_pos: Tuple[Tuple[int, ...], ...], num_parent_params: int
     ) -> ParametersMap:
-        """
-        Build the ParametersMap that will map according to the given child pos to parent
-        pos map given.
-
-        If a child pos maps to None, then the value of the parameter must be 0 at all
-        times.
-        """
-
         def param_map(param: Parameters) -> Parameters:
-            new_params = [0 for _ in range(num_parent_params)]
+            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]
-                if parent_pos is None:
-                    assert value == 0
-                    continue
                 for p in parent_pos:
-                    assert new_params[p] == 0
+                    assert new_params[p] is None or new_params[p] == value
                     new_params[p] = value
-            return tuple(new_params)
+            assert all(p is not None for p in new_params)
+            return cast(Parameters, tuple(new_params))
 
         return param_map
 
-    def _build_children_param_maps(
+    def _build_parent_param_map(
         self,
         parent: CombinatorialClassType,
         children: Tuple[CombinatorialClassType, ...],
-    ) -> Tuple[ParametersMap, ...]:
-        map_list: List[ParametersMap] = []
-        num_parent_params = len(parent.extra_parameters)
-        parent_param_to_pos = {
-            param: pos for pos, param in enumerate(parent.extra_parameters)
+    ) -> ParametersMap:
+        extra_param = self.extra_parameters[self.idx]
+        child = children[self.idx]
+        child_param_to_pos = {
+            param: pos for pos, param in enumerate(child.extra_parameters)
         }
-        for child, extra_param in zip(children, self.extra_parameters):
-            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)
-            child_pos_to_parent_pos: Tuple[Optional[Tuple[int, ...]], ...] = tuple(
-                tuple(
-                    map(
-                        parent_param_to_pos.__getitem__,
-                        reversed_extra_param[child_param],
-                    )
-                )
-                if child_param in reversed_extra_param
-                else None
-                for child_param in child.extra_parameters
-            )
-            map_list.append(
-                self._build_param_map(child_pos_to_parent_pos, num_parent_params)
+        parent_pos_to_child_pos: Tuple[Tuple[int, ...], ...] = tuple(
+            (child_param_to_pos[extra_param[parent_var]],)
+            if parent_var in extra_param
+            else tuple()
+            for parent_var in parent.extra_parameters
+        )
+        return self.build_param_map(
+            parent_pos_to_child_pos, len(child.extra_parameters)
+        )
+
+    def _new_param(self, *children_params: Parameters) -> Parameters:
+        """
+        Computes the parameter values on the parent that the given children parameters
+        contribute to.
+        """
+        mapped_params = (
+            pmap(p) for pmap, p in zip(self._children_param_maps, children_params)
+        )
+        return tuple(sum(vals) for vals in zip(*mapped_params))
+
+    def _other_new_param(self, *children_params: Parameters) -> Parameters:
+        """
+        Computes the parameter values on the parent that the given children parameters
+        contribute to.
+        """
+        mapped_params = (
+            pmap(p)
+            for pmap, p in zip(
+                (m for i, m in enumerate(self._children_param_maps) if i != self.idx),
+                children_params,
             )
-        return tuple(map_list)
+        )
+        return tuple(sum(vals) for vals in zip(*mapped_params))
+
+    def _a(
+        self,
+        n: int,
+        parent_subterm: Callable[[int], Terms],
+        children_subterms: SubTerms,
+    ) -> Terms:
+        """
+        A terms dict representing the polynomial for the rule parent terms minus the
+        contribution of the childs when the size of the flipped child in the composition
+        is smaller than n.
+        """
+        max_sizes = (
+            self._max_sizes[: self.idx] + (n - 1,) + self._max_sizes[self.idx + 1 :]
+        )
+        possible_sizes = utils.compositions(
+            n + self._parent_shift, self.number_of_children, self._min_sizes, max_sizes
+        )
+        res: Terms = Counter(parent_subterm(n + self._parent_shift))
+        for sizes in possible_sizes:
+            for param_value_pairs in CartesianProduct.params_value_pairs_combinations(
+                sizes, children_subterms
+            ):
+                new_param = self._new_param(*(p for p, _ in param_value_pairs))
+                res[new_param] -= utils.prod((v for _, v in param_value_pairs))
+                assert res[new_param] >= 0
+        return res
+
+    def _c(
+        self,
+        children_subterms: SubTerms,
+    ) -> Terms:
+        """
+        Dict representing the polynomial that multiply the polynomial of the
+        flipped child for length n.
+        """
+        min_sizes = self._min_sizes[: self.idx] + self._min_sizes[self.idx + 1 :]
+        max_sizes = self._max_sizes[: self.idx] + self._max_sizes[self.idx + 1 :]
+        possible_sizes = utils.compositions(
+            self._parent_shift, self.number_of_children - 1, min_sizes, max_sizes
+        )
+        other_children_subterms = (
+            children_subterms[: self.idx] + children_subterms[self.idx + 1 :]
+        )
+        res: Terms = Counter()
+        for sizes in possible_sizes:
+            for param_value_pairs in CartesianProduct.params_value_pairs_combinations(
+                sizes, other_children_subterms
+            ):
+                new_param = self._other_new_param(*(p for p, _ in param_value_pairs))
+                res[new_param] += utils.prod((v for _, v in param_value_pairs))
+                assert res[new_param] >= 0
+        return res
+
+    def _terms_to_poly(self, terms: Terms) -> Union[sympy.Poly, int]:
+        variables = tuple(sympy.var(f"k_{i}") for i in range(self._num_parent_params))
+        if variables:
+            res = sympy.Poly(sympy.sympify("0"), *variables)
+        else:
+            res = 0
+        for param, value in terms.items():
+            if not variables:
+                res += value
+            else:
+                res += value * utils.prod(
+                    var ** power for var, power in zip(variables, param)
+                )
+        return res
+
+    @staticmethod
+    def _poly_to_terms(poly: Union[sympy.Poly, int]) -> Terms:
+        if isinstance(poly, int):
+            if poly == 0:
+                return Counter()
+            return Counter({tuple(): poly})
+        return Counter(poly.as_dict())
+
+    def _b(
+        self,
+        n: int,
+        parent_subterm: Callable[[int], Terms],
+        children_subterms: SubTerms,
+    ) -> Terms:
+        """
+        The terms of length n for the flipped child.
+        """
+        a = self._a(n, parent_subterm, children_subterms)
+        c = self._c(children_subterms)
+        a_poly = self._terms_to_poly(a)
+        c_poly = self._terms_to_poly(c)
+        if self._num_parent_params > 0:
+            b_poly, remainder = sympy.div(a_poly, c_poly, domain="ZZ")
+            assert remainder == 0
+        else:
+            b_poly = a_poly // c_poly
+            assert a_poly % c_poly == 0
+        return self._poly_to_terms(b_poly)
+
+    def get_terms(
+        self, parent_terms: Callable[[int], Terms], subterms: SubTerms, n: int
+    ) -> Terms:
+        new_terms: Terms = Counter()
+        children_subterms = (
+            subterms[1 : self.idx + 1] + (parent_terms,) + subterms[self.idx + 1 :]
+        )
+        for param, value in self._b(n, subterms[0], children_subterms).items():
+            mapped_params = self._parent_param_map(param)
+            assert mapped_params not in new_terms or new_terms[mapped_params] == value
+            new_terms[mapped_params] = value
+        return new_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)
+
+    def reliance_profile(self, n: int, **parameters: int) -> RelianceProfile:
+        raise NotImplementedError
 
     def get_sub_objects(
         self, subobjs: SubObjects, n: int
     ) -> Iterator[
         Tuple[Parameters, Tuple[List[Optional[CombinatorialObjectType]], ...]]
     ]:
-        res: List[List[Optional[CombinatorialObjectType]]]
-        res = [[None] for _ in range(self.number_of_children)]
-        for i, subobj in enumerate(subobjs):
-            param_map = self._children_param_maps[i]
-            for param, comb_objs in subobj(n).items():
-                res[i] = comb_objs
-                yield (param_map(param), tuple(res))
-            res[i] = [None]
+        raise NotImplementedError
 
     def random_sample_sub_objects(
         self,
@@ -588,35 +593,6 @@ def random_sample_sub_objects(
         subsamplers: SubSamplers,
         subrecs: SubRecs,
         n: int,
-        **parameters: int
-    ):
-        random_choice = randint(1, parent_count)
-        total = 0
-        for (idx, rec), subsampler, extra_params in zip(
-            enumerate(subrecs), subsamplers, self.get_extra_parameters(n, **parameters)
-        ):
-            # if a parent parameter is not mapped to by some child parameter
-            # then it is assumed that the value of the parent parameter must be 0
-            if extra_params is None or any(
-                val != 0 and k in self.zeroes[idx] for k, val in parameters.items()
-            ):
-                continue
-            total += rec(n=n, **extra_params)
-            if random_choice <= total:
-                obj = subsampler(n=n, **extra_params)
-                return (
-                    tuple(None for _ in range(idx))
-                    + (obj,)
-                    + tuple(None for _ in range(len(subrecs) - idx - 1))
-                )
-
-    @staticmethod
-    def get_eq_symbol() -> str:
-        return "="
-
-    @staticmethod
-    def get_op_symbol() -> str:
-        return "+"
-
-    def __str__(self):
-        return "disjoint union"
+        **parameters: int,
+    ) -> Tuple[Optional[CombinatorialObjectType], ...]:
+        raise NotImplementedError
diff --git a/comb_spec_searcher/strategies/constructor/disjoint.py b/comb_spec_searcher/strategies/constructor/disjoint.py
new file mode 100644
index 000000000..321d793fb
--- /dev/null
+++ b/comb_spec_searcher/strategies/constructor/disjoint.py
@@ -0,0 +1,337 @@
+from collections import defaultdict
+from random import randint
+from typing import Callable, Counter, Dict, Iterator, List, Optional, Tuple
+
+import sympy
+
+from comb_spec_searcher.typing import (
+    CombinatorialClassType,
+    CombinatorialObjectType,
+    Parameters,
+    ParametersMap,
+    RelianceProfile,
+    SubObjects,
+    SubRecs,
+    SubSamplers,
+    SubTerms,
+    Terms,
+)
+
+from .base import Constructor
+
+
+class DisjointUnion(Constructor[CombinatorialClassType, CombinatorialObjectType]):
+    """
+    The DisjointUnion constructor takes as input the children. Each constructor
+    is unique up to the length of the children being used to count.
+
+    Extra parameters are passed on using the parameters dictionaries. Each
+    dictionary's keys should be the extra variable of the child pointing the
+    variable on the parent it came from.
+    If a parents variable does not map to a child, then this variable must be 0
+    as the child contains no occurences.
+
+    The fixed value dictionaries passed will be used ensure that the parameter
+    of a child must take on the given value.
+    """
+
+    def __init__(
+        self,
+        parent: CombinatorialClassType,
+        children: Tuple[CombinatorialClassType, ...],
+        extra_parameters: Optional[Tuple[Dict[str, str], ...]] = None,
+        fixed_values: Optional[Tuple[Dict[str, int], ...]] = None,
+    ):
+        self.number_of_children = len(children)
+        if extra_parameters is not None:
+            self.extra_parameters = extra_parameters
+            assert len(extra_parameters) == len(children)
+        else:
+            assert not parent.extra_parameters
+            self.extra_parameters = tuple(
+                dict() for _ in range(self.number_of_children)
+            )
+        self._children_param_maps = self._build_children_param_maps(parent, children)
+        self.zeroes = tuple(
+            frozenset(parent.extra_parameters) - frozenset(parameter.keys())
+            for parameter in self.extra_parameters
+        )
+        if fixed_values is not None:
+            self.fixed_values = fixed_values
+            assert len(fixed_values) == len(children)
+        else:
+            self.fixed_values = tuple({} for _ in children)
+
+    def get_equation(
+        self, lhs_func: sympy.Function, rhs_funcs: Tuple[sympy.Function, ...]
+    ) -> sympy.Eq:
+        res = 0
+        for rhs_func, extra_parameters in zip(rhs_funcs, self.extra_parameters):
+            res += rhs_func.subs(
+                {child: parent for parent, child in extra_parameters.items()},
+                simultaneous=True,
+            )
+        return sympy.Eq(lhs_func, res)
+
+    def reliance_profile(self, n: int, **parameters: int) -> RelianceProfile:
+        assert not parameters, "only implemented in one variable, namely 'n'"
+        return tuple({"n": (n,)} for _ in range(self.number_of_children))
+
+    def get_extra_parameters(
+        self, n: int, **parameters: int
+    ) -> List[Optional[Dict[str, int]]]:
+        """
+        Will return the extra parameters dictionary based on the parent's
+        parameters. If there is a contradiction, that is some child parameter
+        is given two or more different values that do not match, then None
+        will be returned for that child, indicating that 0 objects on the
+        child match the parents parameters.
+        """
+        res: List[Optional[Dict[str, int]]] = []
+        for i, extra_parameters in enumerate(self.extra_parameters):
+            update_params: Dict[str, int] = {**self.fixed_values[i]}
+            for parent_var, child_var in extra_parameters.items():
+                updated_value = parameters[parent_var]
+                if child_var not in update_params:
+                    update_params[child_var] = updated_value
+                elif update_params[child_var] != updated_value:
+                    break
+            else:
+                res.append(update_params)
+                continue
+            res.append(None)
+        return res
+
+    def get_terms(
+        self, parent_terms: Callable[[int], Terms], subterms: SubTerms, n: int
+    ) -> Terms:
+        new_terms: Terms = Counter()
+        for child_terms, param_map in zip(subterms, self._children_param_maps):
+            for param, value in child_terms(n).items():
+                new_terms[param_map(param)] += value
+        return new_terms
+
+    def _build_children_param_maps(
+        self,
+        parent: CombinatorialClassType,
+        children: Tuple[CombinatorialClassType, ...],
+    ) -> Tuple[ParametersMap, ...]:
+        map_list: List[ParametersMap] = []
+        num_parent_params = len(parent.extra_parameters)
+        parent_param_to_pos = {
+            param: pos for pos, param in enumerate(parent.extra_parameters)
+        }
+        for child, extra_param in zip(children, self.extra_parameters):
+            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)
+            child_pos_to_parent_pos: Tuple[Tuple[int, ...], ...] = tuple(
+                tuple(
+                    map(
+                        parent_param_to_pos.__getitem__,
+                        reversed_extra_param[child_param],
+                    )
+                )
+                if child_param in reversed_extra_param
+                else tuple()
+                for child_param in child.extra_parameters
+            )
+            map_list.append(
+                self.build_param_map(child_pos_to_parent_pos, num_parent_params)
+            )
+        return tuple(map_list)
+
+    def get_sub_objects(
+        self, subobjs: SubObjects, n: int
+    ) -> Iterator[
+        Tuple[Parameters, Tuple[List[Optional[CombinatorialObjectType]], ...]]
+    ]:
+        res: List[List[Optional[CombinatorialObjectType]]]
+        res = [[None] for _ in range(self.number_of_children)]
+        for i, subobj in enumerate(subobjs):
+            param_map = self._children_param_maps[i]
+            for param, comb_objs in subobj(n).items():
+                res[i] = comb_objs
+                yield (param_map(param), tuple(res))
+            res[i] = [None]
+
+    def random_sample_sub_objects(
+        self,
+        parent_count: int,
+        subsamplers: SubSamplers,
+        subrecs: SubRecs,
+        n: int,
+        **parameters: int,
+    ) -> Tuple[Optional[CombinatorialObjectType], ...]:
+        random_choice = randint(1, parent_count)
+        total = 0
+        for (idx, rec), subsampler, extra_params in zip(
+            enumerate(subrecs), subsamplers, self.get_extra_parameters(n, **parameters)
+        ):
+            # if a parent parameter is not mapped to by some child parameter
+            # then it is assumed that the value of the parent parameter must be 0
+            if extra_params is None or any(
+                val != 0 and k in self.zeroes[idx] for k, val in parameters.items()
+            ):
+                continue
+            total += rec(n=n, **extra_params)
+            if random_choice <= total:
+                obj = subsampler(n=n, **extra_params)
+                return (
+                    tuple(None for _ in range(idx))
+                    + (obj,)
+                    + tuple(None for _ in range(len(subrecs) - idx - 1))
+                )
+
+    @staticmethod
+    def get_eq_symbol() -> str:
+        return "="
+
+    @staticmethod
+    def get_op_symbol() -> str:
+        return "+"
+
+    def __str__(self):
+        return "disjoint union"
+
+
+class Complement(Constructor[CombinatorialClassType, CombinatorialObjectType]):
+    """
+    This constructor implements the complement operation. A rule A -> (B, C, D)
+    should be read as A = B - C - D.
+
+    This comes from a disjoint union rule B = A + C + D
+
+    The initialiser should be made using the original disjoint union rule.
+
+    The parent is B, idx is the index of the child the rule is being rearranged
+    for and the extra parameters maps from parent to children.
+    """
+
+    def __init__(
+        self,
+        parent: CombinatorialClassType,
+        children: Tuple[CombinatorialClassType, ...],
+        idx: int,
+        extra_parameters: Optional[Tuple[Dict[str, str], ...]] = None,
+    ):
+        self.number_of_children = len(children)
+        if extra_parameters is not None:
+            self.extra_parameters = extra_parameters
+            assert len(extra_parameters) == len(children)
+        else:
+            assert not parent.extra_parameters
+            self.extra_parameters = tuple(
+                dict() for _ in range(self.number_of_children)
+            )
+        self.idx = idx
+        self._children_param_maps = self._build_children_param_maps(parent, children)
+        self._parent_param_map = self._build_parent_param_map(parent, children, idx)
+
+    def _build_children_param_maps(
+        self,
+        parent: CombinatorialClassType,
+        children: Tuple[CombinatorialClassType, ...],
+    ) -> Tuple[ParametersMap, ...]:
+        map_list: List[ParametersMap] = []
+        num_parent_params = len(parent.extra_parameters)
+        parent_param_to_pos = {
+            param: pos for pos, param in enumerate(parent.extra_parameters)
+        }
+        for child, extra_param in zip(children, self.extra_parameters):
+            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)
+            child_pos_to_parent_pos: Tuple[Tuple[int, ...], ...] = tuple(
+                tuple(
+                    map(
+                        parent_param_to_pos.__getitem__,
+                        reversed_extra_param[child_param],
+                    )
+                )
+                if child_param in reversed_extra_param
+                else tuple()
+                for child_param in child.extra_parameters
+            )
+            map_list.append(
+                self.build_param_map(child_pos_to_parent_pos, num_parent_params)
+            )
+        return tuple(map_list)
+
+    def _build_parent_param_map(
+        self,
+        parent: CombinatorialClassType,
+        children: Tuple[CombinatorialClassType, ...],
+        idx: int,
+    ):
+        extra_param = self.extra_parameters[idx]
+        child = children[self.idx]
+        child_param_to_pos = {
+            param: pos for pos, param in enumerate(child.extra_parameters)
+        }
+        parent_pos_to_child_pos: Tuple[Tuple[int, ...], ...] = tuple(
+            (child_param_to_pos[extra_param[parent_var]],)
+            if parent_var in extra_param
+            else tuple()
+            for parent_var in parent.extra_parameters
+        )
+        return self.build_param_map(
+            parent_pos_to_child_pos, len(child.extra_parameters)
+        )
+
+    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)
+
+    def reliance_profile(self, n: int, **parameters: int) -> RelianceProfile:
+        raise NotImplementedError
+
+    def get_terms(
+        self, parent_terms: Callable[[int], Terms], subterms: SubTerms, n: int
+    ) -> Terms:
+        parent_terms_mapped: Terms = Counter()
+        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)
+
+        return parent_terms_mapped
+
+    def get_sub_objects(
+        self, subobjs: SubObjects, n: int
+    ) -> Iterator[
+        Tuple[Parameters, Tuple[List[Optional[CombinatorialObjectType]], ...]]
+    ]:
+        raise NotImplementedError
+
+    def random_sample_sub_objects(
+        self,
+        parent_count: int,
+        subsamplers: SubSamplers,
+        subrecs: SubRecs,
+        n: int,
+        **parameters: int,
+    ) -> Tuple[Optional[CombinatorialObjectType], ...]:
+        raise NotImplementedError
diff --git a/comb_spec_searcher/strategies/rule.py b/comb_spec_searcher/strategies/rule.py
index 7390d3664..70b52c9c1 100644
--- a/comb_spec_searcher/strategies/rule.py
+++ b/comb_spec_searcher/strategies/rule.py
@@ -6,9 +6,10 @@
 A CombinatorialSpecification is (more or less) a set of Rule.
 """
 import abc
+import random
 from collections import defaultdict
+from importlib import import_module
 from itertools import chain, product
-from random import randint
 from typing import (
     TYPE_CHECKING,
     Callable,
@@ -19,11 +20,14 @@
     Optional,
     Sequence,
     Tuple,
+    Type,
+    Union,
     cast,
 )
 
 from sympy import Eq, Function
 
+from comb_spec_searcher.combinatorial_class import CombinatorialClass
 from comb_spec_searcher.typing import (
     Objects,
     ObjectsCache,
@@ -35,7 +39,7 @@
 
 from ..combinatorial_class import CombinatorialClassType, CombinatorialObjectType
 from ..exception import SanityCheckFailure, StrategyDoesNotApply
-from .constructor import Constructor, DisjointUnion
+from .constructor import Complement, Constructor, DisjointUnion
 
 if TYPE_CHECKING:
     from .strategy import AbstractStrategy, Strategy, VerificationStrategy
@@ -63,12 +67,30 @@ def __init__(
         self.subgenerators: Optional[
             Tuple[Callable[..., Iterator[CombinatorialObjectType]], ...]
         ] = None
-        self.subsamplers: Optional[Tuple[Callable[..., CombinatorialObjectType], ...]]
+        self.subsamplers: Optional[
+            Tuple[Callable[..., CombinatorialObjectType], ...]
+        ] = None
         self.subterms: Optional[SubTerms] = None
         self.subobjects: Optional[SubObjects] = None
         self._children = children
         self._non_empty_children: Optional[Tuple[CombinatorialClassType, ...]] = None
 
+    def to_jsonable(self) -> dict:
+        d = {
+            "class_module": self.__class__.__module__,
+            "rule_class": self.__class__.__name__,
+        }
+        return d
+
+    @classmethod
+    @abc.abstractmethod
+    def from_dict(cls, d: dict) -> "AbstractRule":
+        module = import_module(d.pop("class_module"))
+        RuleClass: Type["AbstractRule"] = getattr(module, d.pop("rule_class"))
+        if not issubclass(RuleClass, AbstractRule):
+            raise ValueError("Not a valid rule class")
+        return RuleClass.from_dict(d)
+
     @property
     def strategy(
         self,
@@ -152,6 +174,12 @@ def is_equivalence(self) -> bool:
         Returns True if the rule is an equivalence.
         """
 
+    def is_two_way(self) -> bool:
+        """
+        Returns True if it is a two way rule.
+        """
+        return self.strategy.is_two_way(self.comb_class)
+
     def non_empty_children(self) -> Tuple[CombinatorialClassType, ...]:
         """
         Return the tuple of non-empty combinatorial classes that are found
@@ -241,6 +269,12 @@ def __eq__(self, other: object) -> bool:
             and self._strategy == other._strategy
         )
 
+    def get_eq_symbol(self):
+        return self.strategy.get_eq_symbol()
+
+    def get_op_symbol(self):
+        return self.strategy.get_op_symbol()
+
     def __str__(self) -> str:
         def frontpad(res, height):
             n = max(len(s) for s in res)
@@ -266,7 +300,7 @@ def join(res, new):
             symbol_height = min(1, len(res) - 1)
             eq_symbol = (
                 ["     " for i in range(symbol_height)]
-                + ["  {}  ".format(self.strategy.get_eq_symbol())]
+                + ["  {}  ".format(self.get_eq_symbol())]
                 + ["     " for i in range(symbol_height)]
             )
             join(res, eq_symbol)
@@ -274,7 +308,7 @@ def join(res, new):
             if len(children) > 1:
                 op_symbol = (
                     ["     " for i in range(symbol_height)]
-                    + ["  {}  ".format(self.strategy.get_op_symbol())]
+                    + ["  {}  ".format(self.get_op_symbol())]
                     + ["     " for i in range(symbol_height)]
                 )
                 for child in children[1:]:
@@ -293,6 +327,26 @@ def __init__(
         super().__init__(strategy, comb_class, children)
         self._constructor: Optional[Constructor] = None
 
+    def to_jsonable(self) -> dict:
+        d = super().to_jsonable()
+        d["comb_class"] = self.comb_class.to_jsonable()
+        d["children"] = [child.to_jsonable() for child in self.children]
+        d["strategy"] = self.strategy.to_jsonable()
+        return d
+
+    @classmethod
+    def from_dict(cls, d: dict) -> "Rule":
+        # pylint: disable=import-outside-toplevel
+        from comb_spec_searcher.strategies.strategy import AbstractStrategy, Strategy
+
+        strategy = AbstractStrategy.from_dict(d.pop("strategy"))
+        assert isinstance(strategy, Strategy)
+        comb_class = CombinatorialClass.from_dict(d.pop("comb_class"))
+        comb_class = cast(CombinatorialClassType, comb_class)
+        d.pop("children")
+        assert not d
+        return cls(strategy, comb_class)
+
     @property
     def strategy(self) -> "Strategy[CombinatorialClassType, CombinatorialObjectType]":
         return cast(
@@ -317,7 +371,10 @@ def constructor(
         return self._constructor
 
     def is_equivalence(self):
-        return self.strategy.can_be_equivalent() and len(self.non_empty_children()) == 1
+        return (
+            isinstance(self.constructor, (DisjointUnion, Complement))
+            and len(self.non_empty_children()) == 1
+        )
 
     def backward_map(
         self, objs: Tuple[Optional[CombinatorialObjectType], ...]
@@ -347,21 +404,20 @@ def to_equivalence_rule(self) -> "EquivalenceRule":
         ), f"EquivalenceRule can only be created for equivalence rules\n{self}"
         return EquivalenceRule(self)
 
-    def to_reverse_rule(self) -> "ReverseRule":
+    def to_reverse_rule(self, idx: int) -> "Rule":
         """
         Return the reverse rule. At this stage, reverse rules can only be
         created for equivalence rules.
         """
-        assert (
-            self.is_equivalence()
-        ), "reverse rule can only be created for equivalence rules"
-        return ReverseRule(self)
+        return ReverseRule(self, idx)
 
     def _ensure_level(self, n: int) -> None:
         if self.subterms is None:
             raise RuntimeError("set_subrecs must be set first")
         while n >= len(self.terms_cache):
-            terms = self.constructor.get_terms(self.subterms, len(self.terms_cache))
+            terms = self.constructor.get_terms(
+                self.get_terms, self.subterms, len(self.terms_cache)
+            )
             self.terms_cache.append(terms)
 
     def get_equation(
@@ -398,12 +454,15 @@ def random_sample_object_of_size(
             total_count, self.subsamplers, self.subrecs, n, **parameters
         )
         objs = tuple(self.backward_map(subobjs))
-        idx = randint(0, len(objs) - 1)
-        return objs[idx]
+        return random.choice(objs)
 
     def sanity_check(self, n: int) -> bool:
         try:
-            return self._sanity_check_count(n) and self._sanity_check_objects(n)
+            return (
+                self._sanity_check_count(n)
+                and self._sanity_check_objects(n)
+                and self._sanity_check_random_sample(n)
+            )
         except SanityCheckFailure as e:
             raise e
 
@@ -456,6 +515,75 @@ def _sanity_check_objects(self, n: int) -> bool:
             )
         return True
 
+    def _sanity_check_random_sample(self, n: int) -> bool:
+        # pylint: disable=access-member-before-definition
+        # pylint: disable=attribute-defined-outside-init
+        actual_objects = self.comb_class.get_objects(n)
+        possible_parameters = [
+            (dict(zip(self.comb_class.extra_parameters, param)), param)
+            for param in actual_objects
+        ]
+        if not possible_parameters:
+            return True
+
+        def fake_subsampler(
+            comb_class: CombinatorialClassType,
+        ) -> Callable[..., CombinatorialObjectType]:
+            objs_cache: Dict[int, Objects] = {}
+
+            def sampler(n: int, **parameters: int) -> CombinatorialObjectType:
+                if n not in objs_cache:
+                    objs_cache[n] = comb_class.get_objects(n)
+                objs = objs_cache[n]
+                param_tuple = tuple(parameters[p] for p in comb_class.extra_parameters)
+                return cast(CombinatorialObjectType, random.choice(objs[param_tuple]))
+
+            return sampler
+
+        def fake_subrec(comb_class: CombinatorialClassType) -> Callable[..., int]:
+            terms_cache: Dict[int, Terms] = {}
+
+            def subrec(n: int, **parameters: int) -> int:
+                if n not in terms_cache:
+                    terms_cache[n] = comb_class.get_terms(n)
+                terms = terms_cache[n]
+                param_tuple = tuple(parameters[p] for p in comb_class.extra_parameters)
+                return terms[param_tuple]
+
+            return subrec
+
+        tmpsamplers = self.subsamplers
+        tmpsubrec = self.subrecs
+        tmpsubterms = self.subterms
+        self.subsamplers = tuple(fake_subsampler(child) for child in self.children)
+        self.subrecs = tuple(fake_subrec(child) for child in self.children)
+        self.subterms = tuple(child.get_terms for child in self.children)
+        try:
+            self.random_sample_object_of_size(n, **possible_parameters[0][0])
+        except NotImplementedError:
+            # Skipping testing rules that have not implemented random_sampling.
+            self.subsamplers = tmpsamplers
+            self.subrecs = tmpsubrec
+            self.subterms = tmpsubterms
+            return True
+
+        for paramd, paramt in possible_parameters:
+            obj = self.random_sample_object_of_size(n, **paramd)
+            if obj not in actual_objects[paramt]:
+                self.subsamplers = tmpsamplers
+                self.subrecs = tmpsubrec
+                self.subterms = tmpsubterms
+                raise SanityCheckFailure(
+                    f"The following rule failed sanity check:\n"
+                    f"{self}\n"
+                    f"Failed for size {n}:\n"
+                    f"Sampled the object {obj} for the parameters {paramd}\n"
+                )
+        self.subsamplers = tmpsamplers
+        self.subrecs = tmpsubrec
+        self.subterms = tmpsubterms
+        return True
+
 
 class EquivalenceRule(Rule[CombinatorialClassType, CombinatorialObjectType]):
     def __init__(self, rule: Rule):
@@ -465,31 +593,62 @@ def __init__(self, rule: Rule):
         super().__init__(rule.strategy, rule.comb_class, (child,))
         self.child_idx = rule.children.index(child)
         self.actual_children = rule.children
-        self._constructor: Optional[DisjointUnion] = None
+        self.original_rule = rule
+        self._constructor: Optional[Union[DisjointUnion, Complement]] = None
+
+    def to_jsonable(self) -> dict:
+        d = super().to_jsonable()
+        d.pop("comb_class")
+        d.pop("strategy")
+        d.pop("children")
+        d["original_rule"] = self.original_rule.to_jsonable()
+        return d
+
+    @classmethod
+    def from_dict(cls, d: dict) -> "EquivalenceRule":
+        rule = AbstractRule.from_dict(d.pop("original_rule"))
+        assert isinstance(rule, Rule)
+        assert not d, d.keys()
+        return cls(rule)
 
     @property
-    def constructor(self) -> DisjointUnion:
+    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 = cast(
-                DisjointUnion,
-                self.strategy.constructor(self.comb_class, self.actual_children),
-            )
-            assert isinstance(original_constructor, DisjointUnion)
-            self._constructor = DisjointUnion(
-                self.comb_class,
-                self.children,
-                (original_constructor.extra_parameters[self.child_idx],),
-            )
+            original_constructor = self.original_rule.constructor
+            if isinstance(original_constructor, DisjointUnion):
+                self._constructor = DisjointUnion(
+                    self.comb_class,
+                    self.children,
+                    (original_constructor.extra_parameters[self.child_idx],),
+                )
+            elif isinstance(original_constructor, Complement):
+                assert isinstance(self.original_rule, ReverseRule)
+                original_original_constructor = (
+                    self.original_rule.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,
+                )
+            else:
+                raise NotImplementedError
         return self._constructor
 
     @staticmethod
     def is_equivalence() -> bool:
         return True
 
+    def to_reverse_rule(self, idx: int) -> "EquivalenceRule":
+        assert idx == 0
+        return self.original_rule.to_reverse_rule(self.child_idx).to_equivalence_rule()
+
     @property
     def formal_step(self) -> str:
         return "{} but only the child at index {} is non-empty".format(
@@ -519,15 +678,34 @@ def forward_map(
 
 class EquivalencePathRule(Rule[CombinatorialClassType, CombinatorialObjectType]):
     """
-    A class for shortening a chain of equivalence rules into a single Rule.
+    A class for shortening a chain of disjoint union  or complement unary rules into a
+    single Rule.
     """
 
     def __init__(self, rules: Sequence[Rule]):
         assert all(rule.is_equivalence() for rule in rules)
         super().__init__(rules[0].strategy, rules[0].comb_class, rules[-1].children)
-        self.rules = rules
+        self.rules = tuple(rules)
         self._constructor: Optional[DisjointUnion] = None
 
+    def to_jsonable(self) -> dict:
+        d = super().to_jsonable()
+        d.pop("comb_class")
+        d.pop("strategy")
+        d.pop("children")
+        d["rules"] = [r.to_jsonable() for r in self.rules]
+        return d
+
+    @classmethod
+    def from_dict(cls, d: dict) -> "EquivalencePathRule":
+        rules: List[Rule] = []
+        for rule_dict in d.pop("rules"):
+            rule = AbstractRule.from_dict(rule_dict)
+            assert isinstance(rule, Rule)
+            rules.append(rule)
+        assert not d
+        return cls(rules)
+
     @property
     def constructor(self) -> DisjointUnion:
         if self._constructor is None:
@@ -536,7 +714,7 @@ def constructor(self) -> DisjointUnion:
             }
             for rule in self.rules:
                 original_constructor = rule.constructor
-                assert isinstance(original_constructor, DisjointUnion)
+                assert isinstance(original_constructor, (DisjointUnion, Complement))
                 rules_parameters = original_constructor.extra_parameters[0]
                 extra_parameters = {
                     parent_var: rules_parameters[child_var]
@@ -588,7 +766,7 @@ def forward_map(
         self, obj: CombinatorialObjectType
     ) -> Tuple[Optional[CombinatorialObjectType], ...]:
         res = cast(CombinatorialObjectType, obj)
-        for rule in reversed(self.rules):
+        for rule in self.rules:
             res = cast(CombinatorialObjectType, rule.forward_map(res)[0])
         return (res,)
 
@@ -632,55 +810,83 @@ def join(res, new):
 
 class ReverseRule(Rule[CombinatorialClassType, CombinatorialObjectType]):
     """
-    A class for creating the reverse of an equivalence rule.
+    A class for creating a reverse equivalence rule.
     """
 
-    def __init__(self, rule: Rule):
-        assert (
-            rule.is_equivalence()
-        ), "reversing a rule only works for equivalence rules"
-        super().__init__(rule.strategy, rule.children[0], (rule.comb_class,))
+    def __init__(self, rule: Rule, idx: int):
+        assert rule.is_two_way()
+        super().__init__(
+            rule.strategy,
+            rule.children[idx],
+            (rule.comb_class, *rule.children[:idx], *rule.children[idx + 1 :]),
+        )
         self.original_rule = rule
-        self._constructor: Optional[DisjointUnion] = None
+        self.idx = idx  # the idx of the child to be counted.
+        self._constructor: Optional[Constructor] = None
+
+    def to_jsonable(self) -> dict:
+        d = super().to_jsonable()
+        d.pop("comb_class")
+        d.pop("strategy")
+        d.pop("children")
+        d["original_rule"] = self.original_rule.to_jsonable()
+        d["idx"] = self.idx
+        return d
+
+    @classmethod
+    def from_dict(cls, d: dict) -> "ReverseRule":
+        rule = AbstractRule.from_dict(d.pop("original_rule"))
+        assert isinstance(rule, Rule)
+        idx = d.pop("idx")
+        assert not d
+        return cls(rule, idx)
+
+    def to_reverse_rule(self, idx: int) -> "Rule":
+        raise NotImplementedError("You don't want to do that! I promise")
 
     @property
-    def constructor(self) -> DisjointUnion:
+    def constructor(
+        self,
+    ) -> Constructor:
         if self._constructor is None:
-            original_constructor = self.original_rule.constructor
-            assert isinstance(original_constructor, DisjointUnion), (
-                "reverse rule coming from non disjoint union strategy - "
-                "you'll need to update the ReverseRule constructor!"
+            self._constructor = self.strategy.reverse_constructor(
+                self.idx,
+                self.original_rule.comb_class,
+                self.original_rule.children,
             )
-            flipped_extra_params = {
-                b: a for a, b in original_constructor.extra_parameters[0].items()
-            }
-            fixed_values = {
-                k: 0
-                for k in self.children[0].extra_parameters
-                if k not in flipped_extra_params.values()
-            }
-            self._constructor = DisjointUnion(
-                self.comb_class, self.children, (flipped_extra_params,), (fixed_values,)
-            )
-        return self._constructor
 
-    @staticmethod
-    def is_equivalence() -> bool:
-        return True
+        return self._constructor
 
     @property
     def formal_step(self) -> str:
-        return "reverse of '{}'".format(self.strategy.formal_step())
+        if len(self.children) == 1:
+            return "reverse of '{}'".format(self.strategy.formal_step())
+        return f"reverse of '{self.strategy.formal_step()}', counting child {self.idx}"
 
     def backward_map(
         self, objs: Tuple[Optional[CombinatorialObjectType], ...]
     ) -> Iterator[CombinatorialObjectType]:
-        yield cast(CombinatorialObjectType, self.original_rule.forward_map(objs[0])[0])
+        # 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]
+        )
 
     def forward_map(
         self, obj: CombinatorialObjectType
     ) -> Tuple[Optional[CombinatorialObjectType], ...]:
-        return (next(self.original_rule.backward_map((obj,))),)
+        # 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)
+        )
+
+    def get_op_symbol(self):
+        return self.strategy.get_reverse_op_symbol()
 
 
 class VerificationRule(AbstractRule[CombinatorialClassType, CombinatorialObjectType]):
@@ -698,6 +904,27 @@ def __init__(
         assert not children
         super().__init__(strat, comb_class, children)
 
+    def to_jsonable(self) -> dict:
+        d = super().to_jsonable()
+        d["comb_class"] = self.comb_class.to_jsonable()
+        d["strategy"] = self.strategy.to_jsonable()
+        return d
+
+    @classmethod
+    def from_dict(cls, d: dict) -> "VerificationRule":
+        # pylint: disable=import-outside-toplevel
+        from comb_spec_searcher.strategies.strategy import (
+            AbstractStrategy,
+            VerificationStrategy,
+        )
+
+        strategy = AbstractStrategy.from_dict(d.pop("strategy"))
+        assert isinstance(strategy, VerificationStrategy)
+        comb_class = CombinatorialClass.from_dict(d.pop("comb_class"))
+        comb_class = cast(CombinatorialClassType, comb_class)
+        assert not d
+        return cls(strategy, comb_class)
+
     @property
     def strategy(
         self,
diff --git a/comb_spec_searcher/strategies/strategy.py b/comb_spec_searcher/strategies/strategy.py
index deafc3b93..0fa89504c 100644
--- a/comb_spec_searcher/strategies/strategy.py
+++ b/comb_spec_searcher/strategies/strategy.py
@@ -75,7 +75,13 @@
     CombinatorialObjectType,
 )
 from ..exception import InvalidOperationError, ObjectMappingError, StrategyDoesNotApply
-from .constructor import CartesianProduct, Constructor, DisjointUnion
+from .constructor import (
+    CartesianProduct,
+    Complement,
+    Constructor,
+    DisjointUnion,
+    Quotient,
+)
 from .rule import AbstractRule, Rule, VerificationRule
 
 if TYPE_CHECKING:
@@ -166,13 +172,20 @@ 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
         equivalence rule.
         """
 
+    @abc.abstractmethod
+    def is_two_way(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 decomposition_function(
         self, comb_class: CombinatorialClassType
@@ -204,6 +217,10 @@ def get_op_symbol() -> str:
         """
         return "?"
 
+    @staticmethod
+    def get_reverse_op_symbol() -> str:
+        return "?"
+
     def __eq__(self, other: object) -> bool:
         if not isinstance(other, AbstractStrategy):
             return NotImplemented
@@ -310,6 +327,18 @@ def constructor(
         if children is None:
             children = self.decomposition_function(comb_class)
 
+    @abc.abstractmethod
+    def reverse_constructor(
+        self,
+        idx: int,
+        comb_class: CombinatorialClassType,
+        children: Optional[Tuple[CombinatorialClassType, ...]] = None,
+    ) -> Constructor:
+        """
+        This is where the details of the 'reliance profile' and 'counting'
+        functions are hidden.
+        """
+
     @abc.abstractmethod
     def backward_map(
         self,
@@ -393,6 +422,10 @@ def __init__(
     def can_be_equivalent() -> bool:
         return True
 
+    @staticmethod
+    def is_two_way(comb_class: CombinatorialClassType) -> bool:
+        return True
+
     def constructor(
         self,
         comb_class: CombinatorialClassType,
@@ -408,6 +441,20 @@ def constructor(
             extra_parameters=self.extra_parameters(comb_class, children),
         )
 
+    def reverse_constructor(
+        self,
+        idx: int,
+        comb_class: CombinatorialClassType,
+        children: Optional[Tuple[CombinatorialClassType, ...]] = None,
+    ):
+        if children is None:
+            children = self.decomposition_function(comb_class)
+            if children is None:
+                raise StrategyDoesNotApply("Strategy does not apply")
+        return Quotient(
+            comb_class, children, idx, self.extra_parameters(comb_class, children)
+        )
+
     @staticmethod
     def get_op_symbol() -> str:
         """
@@ -416,6 +463,10 @@ def get_op_symbol() -> str:
         """
         return "x"
 
+    @staticmethod
+    def get_reverse_op_symbol() -> str:
+        return "/"
+
 
 class DisjointUnionStrategy(Strategy[CombinatorialClassType, CombinatorialObjectType]):
     """
@@ -444,6 +495,10 @@ def __init__(
     def can_be_equivalent() -> bool:
         return True
 
+    @staticmethod
+    def is_two_way(comb_class: CombinatorialClassType) -> bool:
+        return True
+
     def constructor(
         self,
         comb_class: CombinatorialClassType,
@@ -459,6 +514,20 @@ def constructor(
             extra_parameters=self.extra_parameters(comb_class, children),
         )
 
+    def reverse_constructor(
+        self,
+        idx: int,
+        comb_class: CombinatorialClassType,
+        children: Optional[Tuple[CombinatorialClassType, ...]] = None,
+    ):
+        if children is None:
+            children = self.decomposition_function(comb_class)
+            if children is None:
+                raise StrategyDoesNotApply("Strategy does not apply")
+        return Complement(
+            comb_class, children, idx, self.extra_parameters(comb_class, children)
+        )
+
     @staticmethod
     def backward_map_index(objs: Tuple[Optional[CombinatorialObjectType], ...]) -> int:
         """
@@ -496,6 +565,10 @@ def get_op_symbol() -> str:
         """
         return "+"
 
+    @staticmethod
+    def get_reverse_op_symbol() -> str:
+        return "-"
+
 
 class SymmetryStrategy(
     DisjointUnionStrategy[CombinatorialClassType, CombinatorialObjectType]
@@ -563,6 +636,10 @@ def __call__(
     def can_be_equivalent() -> bool:
         return False
 
+    @staticmethod
+    def is_two_way(comb_class: CombinatorialClassType) -> bool:
+        return False
+
     def pack(self, comb_class: CombinatorialClassType) -> "StrategyPack":
         """
         Returns a StrategyPack that finds a proof tree for the comb_class in
@@ -740,7 +817,9 @@ def __str__(self) -> str:
         return "verify atoms"
 
 
-class EmptyStrategy(VerificationStrategy[CombinatorialClass, CombinatorialObject]):
+class EmptyStrategy(
+    VerificationStrategy[CombinatorialClassType, CombinatorialObjectType]
+):
     """
     A subclass for when a combinatorial class is equal to the empty set.
     """
@@ -758,8 +837,8 @@ def get_objects(comb_class: CombinatorialClass, n: int) -> Objects:
 
     def get_genf(
         self,
-        comb_class: CombinatorialClass,
-        funcs: Optional[Dict[CombinatorialClass, Function]] = None,
+        comb_class: CombinatorialClassType,
+        funcs: Optional[Dict[CombinatorialClassType, Function]] = None,
     ) -> Integer:
         if not self.verified(comb_class):
             raise StrategyDoesNotApply(
@@ -769,8 +848,8 @@ def get_genf(
 
     @staticmethod
     def random_sample_object_of_size(
-        comb_class: CombinatorialClass, n: int, **parameters: int
-    ) -> CombinatorialObject:
+        comb_class: CombinatorialClassType, n: int, **parameters: int
+    ) -> CombinatorialObjectType:
         raise StrategyDoesNotApply("Can't sample from empty set.")
 
     @staticmethod
diff --git a/comb_spec_searcher/tree_searcher.py b/comb_spec_searcher/tree_searcher.py
index 19565de8b..b73447827 100644
--- a/comb_spec_searcher/tree_searcher.py
+++ b/comb_spec_searcher/tree_searcher.py
@@ -232,7 +232,7 @@ def _dfs_forest(
                 new_maximum = maximum - length if maximum is not None else None
                 for seen2, trees in _dfs_forest(roots, seen1, new_maximum):
                     actual_length = length + sum(len(t) for t in trees)
-                    if maximum is not None and actual_length < maximum:
+                    if maximum is None or actual_length < maximum:
                         yield seen1.union(seen2), [tree] + trees
 
     sorted_rules_dict = {
diff --git a/comb_spec_searcher/typing.py b/comb_spec_searcher/typing.py
index 091fe586c..2e89aaff6 100644
--- a/comb_spec_searcher/typing.py
+++ b/comb_spec_searcher/typing.py
@@ -20,7 +20,23 @@
 __all__ = [
     "CombinatorialObjectType",
     "CombinatorialClassType",
-    "SpecificationClassesAndStrats",
+    "ClassKey",
+    "Key",
+    "CSSstrategy",
+    "RulesDict",
+    "WorkPacket",
+    "RuleKey",
+    "Parameters",
+    "ParametersMap",
+    "RelianceProfile",
+    "Objects",
+    "ObjectsCache",
+    "Terms",
+    "TermsCache",
+    "SubObjects",
+    "SubRecs",
+    "SubSamplers",
+    "SubTerms",
 ]
 
 CombinatorialObjectType = TypeVar(
@@ -30,12 +46,8 @@
 CombinatorialClassType = TypeVar("CombinatorialClassType", bound="CombinatorialClass")
 
 
-SpecificationClassesAndStrats = Tuple[
-    List[Tuple[CombinatorialClassType, "AbstractStrategy"]],
-    List[List[CombinatorialClassType]],
-]
 ClassKey = Union[bytes, CombinatorialClassType]
-Key = Union["CombinatorialClassType", int]
+Key = Union[CombinatorialClassType, int]
 
 CSSstrategy = Union["AbstractStrategy", "StrategyFactory"]
 RulesDict = Dict[int, Set[Tuple[int, ...]]]
@@ -47,9 +59,6 @@ class WorkPacket(NamedTuple):
     inferral: bool
 
 
-SpecificationLabelsAndStrats = Tuple[
-    List[Tuple[int, "AbstractStrategy"]], List[List[int]]
-]
 RuleKey = Tuple[int, Tuple[int, ...]]
 
 # From constructor
diff --git a/example.py b/example.py
index 981877368..43b44dc7d 100644
--- a/example.py
+++ b/example.py
@@ -123,12 +123,13 @@ def __hash__(self) -> int:
         )
 
     def __str__(self) -> str:
+        prefix = self.prefix if self.prefix else '""'
         if self.just_prefix:
-            return "The word {}".format(self.prefix)
+            return "The word {}".format(prefix)
         return "Words over {{{}}} avoiding {{{}}} with prefix {}" "".format(
             ", ".join(self.alphabet),
             ", ".join(self.patterns),
-            self.prefix if self.prefix else '""',
+            prefix,
         )
 
     def __repr__(self) -> str:
diff --git a/setup.py b/setup.py
index f35b9c2dc..860b36867 100755
--- a/setup.py
+++ b/setup.py
@@ -10,7 +10,7 @@ def read(fname):
 
 setup(
     name="comb_spec_searcher",
-    version="2.4.0",
+    version="3.0.0",
     author="Permuta Triangle",
     author_email="permutatriangle@gmail.com",
     description="A library for performing combinatorial exploration.",
@@ -41,11 +41,11 @@ def read(fname):
         "Topic :: Scientific/Engineering :: Mathematics",
     ],
     install_requires=[
-        "logzero==1.6.2",
-        "sympy==1.6.2",
-        "psutil==5.7.3",
+        "logzero==1.6.3",
+        "sympy==1.7.1",
+        "psutil==5.8.0",
         "pympler==0.9",
-        "requests==2.24.0",
+        "requests==2.25.1",
         "typing-extensions==3.7.4.3",
         "tabulate==0.8.7",
     ],
diff --git a/tox.ini b/tox.ini
index ebabfe8eb..e9740e9ac 100644
--- a/tox.ini
+++ b/tox.ini
@@ -21,10 +21,10 @@ basepython =
     py39: python3.9
     pypy36: pypy3
 deps =
-    pytest==6.1.2
+    pytest==6.2.1
     pytest-repeat==0.9.1
     docutils==0.16
-    Pygments==2.7.2
+    Pygments==2.7.3
 commands = pytest
 
 [pytest]


Test patch:
diff --git a/tests/test_rule.py b/tests/test_rule.py
new file mode 100644
index 000000000..2dd6b6a17
--- /dev/null
+++ b/tests/test_rule.py
@@ -0,0 +1,15 @@
+from example import AvoidingWithPrefix, ExpansionStrategy
+
+
+def test_reverse_equivalence():
+    comb_class = AvoidingWithPrefix("aaa", ["aaaa"], "a", False)
+    rule = ExpansionStrategy()(comb_class)
+    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]
+    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]