From 350d0f6fb2a98e9997bee17699cb497d5eb5de7e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Radek=20Hu=C5=A1ek?= <husek@iuuk.mff.cuni.cz>
Date: Wed, 3 Jun 2020 09:58:27 +0200
Subject: [PATCH] refactor a lot

---
 analyze-cut.py             |   2 +-
 analyze-cycle.py           |   5 +-
 count_cdc.py               |   2 +-
 flower_snarks.py           | 485 -------------------------------------
 graph_parameters.py        | 347 --------------------------
 graph_tools/__init__.py    |   1 +
 graph_tools/__main__.py    |  27 +++
 graph_tools/all.py         |   5 +
 graph_tools/definitions.py | 243 +++++++++++++++++++
 graph_tools/misc.py        | 114 +++++++++
 graph_tools/parameters.py  | 312 ++++++++++++++++++++++++
 graph_tools/sequences.py   | 189 +++++++++++++++
 graph_tools/tests.py       |  21 ++
 graph_tools/utils.py       | 146 +++++++++++
 test_cycle.py              |   7 +-
 utils.py                   | 146 -----------
 16 files changed, 1069 insertions(+), 983 deletions(-)
 delete mode 100755 flower_snarks.py
 delete mode 100644 graph_parameters.py
 create mode 100644 graph_tools/__init__.py
 create mode 100644 graph_tools/__main__.py
 create mode 100644 graph_tools/all.py
 create mode 100644 graph_tools/definitions.py
 create mode 100644 graph_tools/misc.py
 create mode 100644 graph_tools/parameters.py
 create mode 100644 graph_tools/sequences.py
 create mode 100644 graph_tools/tests.py
 create mode 100644 graph_tools/utils.py
 delete mode 100644 utils.py

diff --git a/analyze-cut.py b/analyze-cut.py
index e8bc688..7ce1d1c 100644
--- a/analyze-cut.py
+++ b/analyze-cut.py
@@ -1,4 +1,4 @@
-from flower_snarks import *
+from graph_tools.parameters import *
 from sage.all import *
 from collections import namedtuple
 from itertools import *
diff --git a/analyze-cycle.py b/analyze-cycle.py
index b27feaf..6e91900 100644
--- a/analyze-cycle.py
+++ b/analyze-cycle.py
@@ -1,4 +1,5 @@
-from flower_snarks import *
+from graph_tools.definitions import *
+from graph_tools.parameters import *
 from sage.all import *
 from collections import namedtuple
 from itertools import *
@@ -29,6 +30,8 @@ def calc(s):
   ret.map = lambda x: _bmap[x]
 
   return ret
+
+
 def C_k(k):  
   return Gadget.join(
     [CUBIC_VERTEX]*k,
diff --git a/count_cdc.py b/count_cdc.py
index 3c53e37..8ed6121 100755
--- a/count_cdc.py
+++ b/count_cdc.py
@@ -3,7 +3,7 @@
 import sys
 import os
 import json
-from flower_snarks import *
+from graph_tools.misc import count_cdc_naive
 from sage.all import Graph
 from parmap import parmap
 
diff --git a/flower_snarks.py b/flower_snarks.py
deleted file mode 100755
index 2b4023d..0000000
--- a/flower_snarks.py
+++ /dev/null
@@ -1,485 +0,0 @@
-#!/usr/bin/python3
-
-from fractions import Fraction
-
-from utils import *
-from graph_parameters import *
-
-class Gadget:
-  def __init__(self, size):
-    self._size = size
-
-  def __hash__(self): return id(self)
-  def __eq__(self, b): return self is b
-
-  def size(self): return self._size
-
-  def is_graph(self):
-    return self._size == 0
-
-  def __str__(self):
-    return "Gadget{%i}" % (self._size,)
-
-  def __repr__(self):
-    return self.__str__()
-
-  _gadget_cache = DynamicLRU()
-
-  @staticmethod
-  def join(gadgets, joins, outs):
-    assert all( isinstance(g, Gadget) for g in gadgets ), \
-      "Not all gadgets of the same size: %s" % gadgets
-
-    info = (tuple(gadgets), tuple(joins), tuple(outs))
-    return Gadget._gadget_cache.get(info, lambda: JoinGadget(*info))
-
-class FakeGadget(Gadget):
-  def __init__(self, size, value):
-    super().__init__(size)
-    self.value = value
-
-  def eval(self, _):
-    return self.value
-
-CUBIC_VERTEX = None
-FREE_EDGE = None
-
-class BaseGadget(Gadget):
-  def eval(self, parameter):
-    try:
-      if self is CUBIC_VERTEX:
-        return parameter.CUBIC_VERTEX
-      elif self is FREE_EDGE:
-        return parameter.FREE_EDGE
-    except: pass
-    raise RuntimeError("Unsupported base gadget")
-
-CUBIC_VERTEX = BaseGadget(3)
-FREE_EDGE = BaseGadget(2)
-
-
-class JoinGadget(Gadget):
-  _cache = DynamicLRU()
-
-  def __init__(self, gadgets, joins, outs):
-    super().__init__(len(outs))
-
-    self.offsets = [0]
-    for g in gadgets:
-      self.offsets.append(self.offsets[-1] + g.size())
-
-    joins = tuple(
-      (self.offsets[g1-1] + i1 - 1, self.offsets[g2-1] + i2 - 1)
-      for (g1, i1), (g2, i2) in joins
-    )
-    outs = tuple( self.offsets[g-1] + i - 1 for g, i in outs )
-
-    self._info = (tuple(gadgets), joins, outs)
-
-  def eval(self, parameter):
-    gadgets, joins, outs = self._info
-
-    def do_eval():
-      ret = parameter.eval_join(
-        [ g.eval(parameter) for g in gadgets ],
-        gadgets, joins, outs, self.offsets
-      )
-      if self.is_graph():
-        return parameter.finalize(ret)
-      return ret
-
-    return self._cache.get((parameter, self._info), do_eval)
-
-
-
-
-def ParametrizedGraphSequence(cls):
-  class GraphSequence(cls, GraphSequenceBase):
-    def __init__(self, *args, **kwargs):
-      GraphSequenceBase.__init__(self)
-      cls.__init__(self, *args, **kwargs)
-  return GraphSequence
-
-def GraphSequence(cls):
-  return ParametrizedGraphSequence(cls)()
-
-class GraphSequenceBase:
-  sequence_start = 0
-
-  def __init__(self):
-    self._gadget_cache = DynamicLRU()
-    self._graph_cache = DynamicLRU()
-
-  def _join_gadgets(self, a, b):
-    return Gadget.join([a, b], self.step_join, self.step_out)
-
-  def _next_gadget(self, gadget):
-    return self._join_gadgets(gadget , self.step_gadget)
-
-  def _make_graph(self, gadget):
-    return Gadget.join([ gadget ], self.final_join, [])
-
-  def gadget(self, n):
-    assert type(n) is int and n >= self.sequence_start, "n must be non-negative integer"
-
-    def make_gadget():
-      if n == self.sequence_start: return self.base_gadget
-      return self._next_gadget(self.gadget(n-1))
-
-    return self._gadget_cache.get(n, make_gadget)
-
-  def graph(self, n):
-    return self._graph_cache.get(n, lambda: self._make_graph(self.gadget(n)))
-
-  def stabilize(self, parameter):
-    from sage.all import matrix, QQ
-    boundaries = set()
-    new_boundaries = set( b.boundary for b in self.base_gadget.eval(parameter) )
-
-    i = 0
-    while True:
-      i += 1
-      new_boundaries.update(boundaries)
-      if new_boundaries == boundaries: break
-      boundaries = new_boundaries
-      gadget = FakeGadget(
-        self.base_gadget.size(),
-        [ Boundary(k, 1) for k in boundaries ]
-      )
-      out = self._next_gadget(gadget).eval(parameter)
-      new_boundaries = set( b.boundary for b in out )
-      print("Loop %i done - %i boundaries (%i new)" % \
-        (i, len(new_boundaries), len(new_boundaries.difference(boundaries))))
-
-    sorted_boundaries = sorted(boundaries)
-    rename = { b: i for i, b in enumerate(sorted_boundaries) }
-
-    size = len(sorted_boundaries)
-
-    m = matrix(QQ, size, size, {
-      (rename[row], rename[col]): value
-      for col in sorted_boundaries
-      for row, value in self._next_gadget(FakeGadget(self.base_gadget.size(), [ Boundary(col, 1) ])).eval(parameter)
-    })
-
-    iv = matrix(QQ, size, 1, {
-      (rename[row], 0): v for row, v in self.base_gadget.eval(parameter)
-    })
-
-    fin = matrix(QQ, 1, size, {
-      (0, rename[col]): self._make_graph(FakeGadget(self.base_gadget.size(), [ Boundary(col, 1) ])).eval(parameter)
-      for col in sorted_boundaries
-    })
-
-    return {
-      'variables': sorted_boundaries,
-      'step_matrix': m,
-      'initial_vector': iv,
-      'finalize': fin
-    }
-
-
-
-@GraphSequence
-class FlowerSnark:
-  sequence_start = 2
-
-  step_join = [
-    ((1, 2), (2, 1)),
-    ((1, 4), (2, 3)),
-    ((1, 6), (2, 5))
-  ]
-  step_out = [
-    (1, 1), (2, 2),
-    (1, 3), (2, 4),
-    (1, 5), (2, 6)
-  ]
-  final_join = [
-    ((1, 2), (1, 1)),
-    ((1, 4), (1, 3)),
-    ((1, 6), (1, 5)),
-  ]
-
-  BASIC_GADGET = Gadget.join(
-    [ CUBIC_VERTEX ] * 4,
-    [
-      ((1, 3), (2, 1)),
-      ((2, 2), (3, 3)),
-      ((2, 3), (4, 1))
-    ],
-    [
-      (1, 1), (1, 2),
-      (3, 1), (3, 2),
-      (4, 2), (4, 3)
-    ]
-  )
-  DOUBLE_GADGET = Gadget.join([ BASIC_GADGET ]*2, step_join, step_out)
-  SPECIAL_GADGET = Gadget.join(
-    [ BASIC_GADGET ] * 2,
-    [
-      ((1, 2), (2, 3)),
-      ((1, 4), (2, 1)),
-      ((1, 6), (2, 5))
-    ],
-    step_out
-  )
-  base_gadget = SPECIAL_GADGET
-  step_gadget = BASIC_GADGET
-
-@GraphSequence
-class FlowerSnarkAlt:
-  sequence_start = 1
-
-  step_join = [
-    ((1, 2), (2, 1)),
-    ((1, 4), (2, 3)),
-    ((1, 6), (2, 5))
-  ]
-  step_out = [
-    (1, 1), (2, 2),
-    (1, 3), (2, 4),
-    (1, 5), (2, 6)
-  ]
-  final_join = [
-    ((1, 2), (1, 1)),
-    ((1, 4), (1, 3)),
-    ((1, 6), (1, 5)),
-  ]
-
-  BASIC_GADGET = Gadget.join(
-    [ CUBIC_VERTEX ] * 4,
-    [
-      ((1, 3), (2, 1)),
-      ((2, 2), (3, 3)),
-      ((2, 3), (4, 1))
-    ],
-    [
-      (3, 1), (1, 2),
-      (1, 1), (3, 2),
-      (4, 2), (4, 3)
-    ]
-  )
-  base_gadget = BASIC_GADGET
-  step_gadget = BASIC_GADGET
-
-@GraphSequence
-class Necklace:
-  sequence_start = 1
-
-  step_join = [ ((1,2), (2, 1)) ]
-  step_out = [ (1,1), (2,2) ]
-  final_join = [ ((1,1), (1,2)) ]
-
-  D = Gadget.join(
-    [ CUBIC_VERTEX ] * 4,
-    [
-      ((1, 1), (2, 1)),
-      ((1, 2), (3, 1)),
-      ((1, 3), (4, 1)),
-      ((2, 2), (3, 2)),
-      ((2, 3), (4, 2))
-    ],
-    [(3,3), (4,3)]
-  )
-  base_gadget = D
-  step_gadget = D
-
-
-@GraphSequence
-class Flower:
-  sequence_start = 2
-
-  SCV = FakeGadget(3, [ Boundary(((1,), (2,), (1, 2), None), 1) ])
-
-  D = Gadget.join(
-    [ SCV, CUBIC_VERTEX ],
-    [ ((1,3), (2,3)) ],
-    [ (1,1), (2,1), (2,2), (1,2) ]
-  )
-
-  step_join = [ ((1,3), (2,2)), ((1,4), (2,1)) ]
-  step_out = [ (1,1), (1,2), (2,3), (2,4) ]
-  final_join = [ ((1,1), (1,4)), ((1,2), (1,3)) ]
-
-  base_gadget = Gadget.join([D, D], step_join, step_out)
-  step_gadget = D
-
-
-@ParametrizedGraphSequence
-class GeneralizedPetersen:
-  sequence_start = 1
-
-  def __init__(self, k):
-    assert k >= 1
-
-    self.k = k
-    self.GADGET = Gadget.join(
-      [ CUBIC_VERTEX ]*2 + [ FREE_EDGE ]*(k-1),
-      [ ((1,1), (2,1)) ],
-      [ (2,2) ] +
-      [ (e + 3, i + 1) for e in range(k-1) for i in range(2) ] +
-      [ (2,3), (1,2), (1,3) ]
-    )
-    self.step_gadget = self.GADGET
-    self.base_gadget = self.GADGET
-
-    self.step_join = [ ((1,i+1), (2,i)) for i in range(1, 2*k+3, 2) ]
-    self.step_out = [ (2 - (i % 2), i) for i in range(1, 2*k+3) ]
-    self.final_join = [ ((1,i+1), (1,i)) for i in range(1, 2*k+3, 2) ]
-
-  def graph(self, n):
-    assert n > self.k
-    return super().graph(n)
-
-
-Petersen = Gadget.join(
-  [ CUBIC_VERTEX ] * 10,
-  [
-    ((1,1), (5,3)),
-    ((1,2), (6,1)),
-    ((1,3), (2,1)),
-    ((2,2), (7,1)),
-    ((2,3), (3,1)),
-    ((3,2), (8,1)),
-    ((3,3), (4,1)),
-    ((4,2), (9,1)),
-    ((4,3), (5,1)),
-    ((5,2), (10,1)),
-    ((6,2), (8,3)),
-    ((6,3), (9,2)),
-    ((7,2), (10,2)),
-    ((7,3), (9,3)),
-    ((8,2), (10,3))
-  ],
-  []
-)
-
-
-def sun(k):
-  return Gadget.join(
-    [ CUBIC_VERTEX ] * k,
-    [ ((i + 1, 2), ((i+1) % k + 1, 1)) for i in range(k) ],
-    [ (i + 1, 3) for i in range(k) ]
-  )
-
-
-def count_cdc_naive(G):
-  """Count circuit double covers of G.
-
-  Requires G to be cubic.
-
-  EXAMPLE:
-  >>> from sage.all import *
-  >>> count_cdc_naive(graphs.PetersenGraph())
-  (52, 16)
-  """
-
-  assert(G.degree() == [3]*G.num_verts())
-
-  vert_map = { v: i for i, v in enumerate(G.vertices()) }
-  vert_used = [ 1 ] * len(vert_map)
-
-  joins = []
-  for u, v, _ in G.edges():
-    u = vert_map[u]
-    v = vert_map[v]
-    joins.append(((u+1, vert_used[u]), (v+1, vert_used[v])))
-    vert_used[u] += 1
-    vert_used[v] += 1
-
-  b = Gadget.join([ CUBIC_VERTEX ] * G.num_verts(), joins, [])
-  assert b.is_graph()
-  return (b.eval(CircuitDoubleCover), 2**(G.num_verts() // 2) // 2)
-
-
-def graph_to_gadget(G, decomposition = None):
-  """Transform graph into a Gadget.
-
-  EXAMPLES:
-  >>> from sage.all import *
-  >>> P = graphs.PetersenGraph()
-  >>> graph_to_gadget(P).eval(CircuitDoubleCover)
-  52
-  >>> d = [ list(range(5)), list(range(5, 10)) ]
-  >>> P.is_isomorphic(graph_to_gadget(P, d).eval(UnderlayingGraph))
-  True
-  >>> graph_to_gadget(P, d).eval(CircuitDoubleCover)
-  52
-  """
-  assert(G.degree() == [3]*G.num_verts())
-
-  verts = { v: i for i, v in enumerate(G.vertices()) }
-  edges = [ (verts[u], verts[v]) for u, v, _ in G.edges() ]
-
-  if decomposition is None:
-    decomposition = G.vertices()
-
-  vert_map = { v: [0, 1, 2] for v in range(len(verts)) }
-
-  def process(dec):
-    nonlocal edges
-
-    if not isinstance(dec, list): return (CUBIC_VERTEX, set([ verts[dec] ]))
-
-    gadgets = [ process(x) for x in dec ]
-    vert_to_gadget = dict()
-    gadget_verts = set()
-
-    for i, (_, vs) in enumerate(gadgets):
-      gadget_verts.update(vs)
-      for v in vs:
-        vert_to_gadget[v] = i
-
-    joins = []
-
-    new_edges = []
-    for u, v in edges:
-      if u not in gadget_verts or v not in gadget_verts:
-        new_edges.append((u, v))
-        continue
-
-      joins.append((
-        (vert_to_gadget[u]+1, vert_map[u].pop()+1),
-        (vert_to_gadget[v]+1, vert_map[v].pop()+1)
-      ))
-
-    edges = new_edges
-
-    outs = []
-    for v in gadget_verts:
-      old = vert_map[v]
-      vert_map[v] = list(range(len(outs), len(outs) + len(old)))
-      for i in old:
-        outs.append((vert_to_gadget[v]+1, i+1))
-
-    return (Gadget.join([ g for g, _ in gadgets ], joins, outs), gadget_verts)
-
-  b, _ = process(decomposition)
-  assert b.is_graph()
-  return b
-
-
-def _tests_():
-  """Helper to run tests.
-
-  Doctest does not find tests in docstrings of decorated classes
-  so we put them here instead.
-
-  >>> all( Necklace.graph(i).eval(CircuitDoubleCover) == 2 * 4**(i-1) for i in range(1, 10) )
-  True
-  >>> GeneralizedPetersen(2).graph(5).eval(CircuitDoubleCover) == 52
-  True
-  >>> FlowerSnark.graph(3).eval(CircuitDoubleCover)
-  104
-  """
-  pass
-
-
-if __name__ == "__main__":
-  import doctest
-  (f, t) = doctest.testmod()
-  print("%s: %i tests of %i failed." % (__file__, f, t))
-  if f > 0:
-    from sys import exit
-    exit(1)
-
diff --git a/graph_parameters.py b/graph_parameters.py
deleted file mode 100644
index d96fb30..0000000
--- a/graph_parameters.py
+++ /dev/null
@@ -1,347 +0,0 @@
-from utils import *
-from collections import namedtuple
-from fractions import Fraction
-from itertools import combinations
-
-Boundary = namedtuple("Boundary", ['boundary', 'value'])
-
-
-class GraphParameterBase:
-  def __hash__(self): return id(self)
-  def __eq__(self, b): return self is b
-
-  def eval_join(self, gadget_boundaries, gadgets, joins, outs, offsets):
-    it = chain.from_iterable(
-      self.join_boundaries(x, offsets) for x in product(*gadget_boundaries)
-    )
-
-    def join(gen, e1, e2):
-      for b in gen:
-        yield from self.join_edges(b, e1, e2)
-
-    for e1, e2 in joins: it = join(it, e1, e2)
-
-    it = ( self.project_and_canonize(outs, b) for b in it )
-    
-    ret = dict()
-    for b in it:
-      if b.boundary in ret:
-        ret[b.boundary] = self.merge_values(ret[b.boundary], b.value)
-      else:
-        ret[b.boundary] = b.value
-
-    return [ Boundary(k, v) for k, v in ret.items() ]
-
-  def merge_values(self, a, b):
-    return a + b
-
-
-class SimpleParameterBase(GraphParameterBase):
-  def join_boundaries(self, x, _):
-    yield Boundary(
-      tuple(chain.from_iterable( b.boundary for b in x )),
-      reduce(lambda acc, b: acc*b.value, x, 1)
-    )
-
-  def join_edges(self, b, e1, e2):
-    if self.is_compatible(b.boundary, e1, e2):
-      yield b
-
-  def project_and_canonize(self, sel, b):
-    return Boundary(select(sel, b.boundary), b.value)
-
-  def finalize(self, b):
-    b = list(b)
-
-    assert len(b) <= 1
-    if len(b) == 0: return 0
-    
-    assert b[0].boundary == tuple()
-    return b[0].value
-
-
-@OneOfKind
-class GroupFlow(SimpleParameterBase):
-  def __init__(self, group, nowherezero):
-    self.group = group
-    self.nowherezero = nowherezero
-    self.FREE_EDGE = [ Boundary((i, i), 1) for i in group if not nowherezero or i != 0 ]
-    self.CUBIC_VERTEX = [
-      Boundary((i, j, l), 1) for i in group for j in group
-      for l in group if i + j + l == 0 and
-      (not nowherezero or i != 0 and j != 0 and l != 0)
-    ]
-
-  def is_compatible(self, boundary, e_a, e_b):
-    return boundary[e_a] == boundary[e_b]
-
-@OneOfKind
-class CycleDoubleCover(SimpleParameterBase):
-  def __init__(self, k):
-    self.k = k
-    self.FREE_EDGE = [ Boundary(((i, j), (i, j)), 1) for i in range(k) for j in range(i+1, k) ]
-    self.CUBIC_VERTEX = [
-      Boundary((tuple(sorted((i, j))), tuple(sorted((i, l))), tuple(sorted((j, l)))), 1)
-      for i in range(k) for j in range(k) for l in range(k)
-      if i != j and j != l and i != l
-    ]
-
-  def is_compatible(self, boundary, e_a, e_b):
-    return boundary[e_a] == boundary[e_b]
-    
-
-@OneOfKind
-class EdgeColoring(SimpleParameterBase):
-  def __init__(self, k):
-    self.k = k
-    self.FREE_EDGE = [ Boundary((i, i), 1) for i in range(k) ]
-    self.CUBIC_VERTEX = [
-      Boundary((i, j, l), 1) for i in range(k) for j in range(k)
-      for l in range(k) if i != j and j != l and i != l
-    ]
-
-  def is_compatible(self, boundary, e_a, e_b):
-    return boundary[e_a] == boundary[e_b]
-    
-
-@OneOfKind
-class VertexColoring(SimpleParameterBase):
-  def __init__(self, k):
-    self.k = k
-    self.CUBIC_VERTEX = [ Boundary((i, i, i), 1) for i in range(k) ]
-
-  def is_compatible(self, boundary, e_a, e_b):
-    return boundary[e_a] != boundary[e_b]
-
-
-@Singleton
-class UnderlayingGraph(GraphParameterBase):
-  from sage.all import Graph
-
-  CUBIC_VERTEX = [
-    Boundary((0, 0, 0), Graph([[0], []], multiedges=True, loops=True, immutable=True))
-  ]
-
-  def join_boundaries(self, x, _):
-    offsets = [0]
-    for b in x:
-      offsets.append(offsets[-1] + b.value.num_verts())
-
-    gg = (
-      b.value.relabel(lambda v: v + offsets[i], inplace=False)
-      for i, b in enumerate(x)
-    )
-
-    yield Boundary(
-      tuple(chain.from_iterable(
-        [ v + offsets[i] for v in b.boundary ] for i, b in enumerate(x)
-      )),
-      reduce(lambda acc, b: acc.union(b), gg)
-    )
-
-  def join_edges(self, b, e1, e2):
-    g = b.value.copy(immutable=False)
-    g.add_edge(b.boundary[e1], b.boundary[e2])
-    yield Boundary(b.boundary, g)
-
-  def project_and_canonize(self, sel, b):
-    return Boundary(select(sel, b.boundary), b.value)
-
-  def merge_values(self, a, b): assert False
-
-  def finalize(self, b):
-    assert len(b) == 1
-    assert b[0].boundary == tuple()
-    return b[0].value
-
-
-@Singleton
-class VertexCount(GraphParameterBase):
-  CUBIC_VERTEX = [ Boundary(None, 1) ]
-  FREE_EDGE = [ Boundary(None, 0) ]
-
-  def join_boundaries(self, x, _):
-    yield Boundary(None, sum( b.value for b in x ))
-
-  def join_edges(self, b, _, __): yield b
-  def project_and_canonize(self, _, b): return b
-  def finalize(self, b): return b[0].value
-
-
-@Singleton
-class HamiltonianCycle(GraphParameterBase):
-  FREE_EDGE = [ Boundary((1, 1), 1), Boundary((None, None), 1) ]
-  CUBIC_VERTEX = [
-    Boundary((1, 1, None), 1),
-    Boundary((1, None, 1), 1),
-    Boundary((None, 1, 1), 1)
-  ]
-
-  def join_boundaries(self, x, _):
-    offsets = [ 0 ]
-    try:
-      for b in x:
-        offsets.append(offsets[-1] + max_(b.boundary))
-    except ValueError:
-      return
-
-    def shift(k, off):
-      if k is None: return None
-      return k + off
-
-    ret_b = tuple(chain(*[
-      [ shift(x, offsets[i]) for x in b.boundary ] for i, b in enumerate(x)
-    ]))
-
-    yield Boundary(ret_b, prod( b.value for b in x ))
-
-  def join_edges(self, b, e1, e2):
-    # print("%s %s %s" % (b, e1, e2))
-    c1 = b.boundary[e1]
-    c2 = b.boundary[e2]
-
-    if c1 is None and c2 is None: yield b
-    if c1 is None or c2 is None: return
-
-    if c1 != c2:
-      if c1 > c2: c1, c2 = c2, c1
-      boundary = tuple( v if v != c2 else c1 for v in b.boundary )
-      yield Boundary(boundary, b.value)
-    elif max_(b.boundary) == 1:
-      yield b
-
-  def project_and_canonize(self, p, b):
-    boundary = select(p, b.boundary)
-    rename = { None: None }
-    i = 1
-    for x in boundary:
-      if x not in rename:
-        rename[x] = i
-        i += 1
-    boundary = tuple( rename[x] for x in boundary )
-    return Boundary(boundary, b.value)
-
-  def finalize(self, b):
-    assert len(b) <= 1
-    if len(b) == 0: return 0
-    assert b[0].boundary == tuple()
-    return b[0].value
-
-
-
-@Singleton
-class CircuitDoubleCover(GraphParameterBase):
-  FREE_EDGE = [ Boundary(((1,2), (1,2), None), Fraction(1, 2)) ]
-  CUBIC_VERTEX = [ Boundary(((1,2), (2,3), (1,3), None), 1) ]
-
-  def join_boundaries(self, x, offsets):
-    def shift(k, off):
-      if k is None: return None
-      if len(k) == 1: return (k[0] + off,)
-      return (k[0] + off, k[1] + off)
-
-    tmp = tuple(zip(*[
-      split_on_none([ shift(k, offsets[i]) for k in b.boundary ])
-      for i, b in enumerate(x)
-    ]))
-    ret_b = tuple(chain(*tmp[0], (None,), *tmp[1]))
-
-    yield Boundary(ret_b, prod( b.value for b in x ))
-
-  def join_edges(self, b, e1, e2):
-    assert len(b.boundary[e1]) == len(b.boundary[e2])
-
-    def do_join(uf):
-      ret_b = []
-      for i, x in enumerate(b.boundary):
-        if x is None:
-          ret_b.append(None)
-        else:
-          t = tuple( uf.find(i) for i in x )
-          if len(t) == 2 and t[0] == t[1]: break
-          ret_b.append(t)
-      else:
-        yield Boundary(tuple(ret_b), b.value)
-
-    if len(b.boundary[e1]) == 1:
-      uf = UnionFind(len(b.boundary))
-      uf.union(b.boundary[e1][0], b.boundary[e2][0])
-      yield from do_join(uf)
-      return
-
-    c1, c2 = b.boundary[e1]
-    c3, c4 = b.boundary[e2]
-
-    for _ in range(2):
-      uf = UnionFind(len(b.boundary))
-      uf.union(c1, c3)
-      uf.union(c2, c4)
-      yield from do_join(uf)
-      c3, c4 = c4, c3
-
-  def project_and_canonize(self, p, bo):
-    x = bo.boundary
-    l, r = split_on_none(x)
-    r = set(r).union(l)
-    l = select(p, l)
-    r = r.difference(l)
-    r = { x for x in r if len(x) == 2 }
-
-    rename = defaultdict(list)
-
-    for i, x in enumerate(l):
-      for a in x: rename[a].append(i)
-
-    rename = sorted( (v, k) for k, v in rename.items() )
-    rename = { x[1]: i + 1 for i, x in enumerate(rename) }
-
-    l = [ tuple(sorted([ rename[a] for a in x ])) for x in l ]
-    r = { tuple(sorted([rename[a], rename[b]])) for a, b in r \
-      if a in rename and b in rename }
-
-    r.difference_update(set(l))
-
-    return Boundary(tuple(l + [None] + sorted(r)), bo.value)
-
-  def finalize(self, b):
-    assert len(b) <= 1
-    if len(b) == 0: return 0
-    assert b[0].boundary == (None,)
-    return b[0].value
-
-  def enumerate_boundaries(self, n):
-    if n == 0: yield (None,)
-    if n <= 1: return
-
-    used = set([1, 2])
-    first_unused = 3
-    rng = list(range(n))
-    
-    def add_internal(b):
-      X = set(combinations(range(1, n + 1), 2)).difference(set(b))
-      b = tuple(b + [None])
-      for x in powerset(X):
-        ret = b + tuple(sorted(x))
-        if self.project_and_canonize(rng, Boundary(ret, 0)).boundary == ret:
-          yield ret
-
-    def add_edge(b, used, fu):
-      if len(b) == n:
-        if len(used) == 0:
-          yield from add_internal(b)
-        return
-
-      s = set(used)
-      if fu <= n: s.add(fu)
-      if fu < n: s.add(fu + 1)
-      for x, y in combinations(sorted(s), 2):
-        if fu + 1 in [x, y] and fu not in [x, y]: continue
-        b_ = b + [(x, y)]
-        s_ = used.symmetric_difference(set([x, y]))
-        fu_ = fu
-        if fu + 1 in [x, y]: fu_ += 2
-        elif fu in [x, y]: fu_ += 1
-        yield from add_edge(b_, s_, fu_)
-
-    yield from add_edge([(1,2)], used, first_unused)
-
diff --git a/graph_tools/__init__.py b/graph_tools/__init__.py
new file mode 100644
index 0000000..8b13789
--- /dev/null
+++ b/graph_tools/__init__.py
@@ -0,0 +1 @@
+
diff --git a/graph_tools/__main__.py b/graph_tools/__main__.py
new file mode 100644
index 0000000..7bc5bb1
--- /dev/null
+++ b/graph_tools/__main__.py
@@ -0,0 +1,27 @@
+import doctest
+from importlib import import_module
+from sys import argv
+
+to_test = [
+  "utils",
+  "definitions",
+  "sequences",
+  "parameters",
+  "misc",
+  "tests"
+]
+
+for mod in to_test:
+  m = import_module("." + mod, __package__)
+  (f, t) = doctest.testmod(m)
+  print("%s: %i tests of %i failed." % (mod, f, t))
+  if f > 0:
+    if argv[-1:] == [ "-d" ]:
+      try:
+        doctest.testmod(m, raise_on_error=True)
+      except doctest.DocTestFailure as e:
+        exec("import pdb; pdb.set_trace()", e.test.globs)
+
+    from sys import exit
+    exit(1)
+
diff --git a/graph_tools/all.py b/graph_tools/all.py
new file mode 100644
index 0000000..321f392
--- /dev/null
+++ b/graph_tools/all.py
@@ -0,0 +1,5 @@
+from .definitions import *
+from .sequences import *
+from .parameters import *
+from .misc import *
+
diff --git a/graph_tools/definitions.py b/graph_tools/definitions.py
new file mode 100644
index 0000000..0d34745
--- /dev/null
+++ b/graph_tools/definitions.py
@@ -0,0 +1,243 @@
+def _init_():
+  global \
+    Boundary, \
+    GraphParameterBase, SimpleParameterBase, \
+    Gadget, FakeGadget, CUBIC_VERTEX, FREE_EDGE, \
+    ParametrizedGraphSequence, GraphSequence
+
+  from . import utils
+  from collections import namedtuple
+  from fractions import Fraction
+  from itertools import combinations, chain, product
+  from functools import reduce
+
+  Boundary = namedtuple("Boundary", ['boundary', 'value'])
+
+
+  class GraphParameterBase:
+    def __hash__(self): return id(self)
+    def __eq__(self, b): return self is b
+
+    def eval_join(self, gadget_boundaries, gadgets, joins, outs, offsets):
+      it = chain.from_iterable(
+        self.join_boundaries(x, offsets) for x in product(*gadget_boundaries)
+      )
+
+      def join(gen, e1, e2):
+        for b in gen:
+          yield from self.join_edges(b, e1, e2)
+
+      for e1, e2 in joins: it = join(it, e1, e2)
+
+      it = ( self.project_and_canonize(outs, b) for b in it )
+
+      ret = dict()
+      for b in it:
+        if b.boundary in ret:
+          ret[b.boundary] = self.merge_values(ret[b.boundary], b.value)
+        else:
+          ret[b.boundary] = b.value
+
+      return [ Boundary(k, v) for k, v in ret.items() ]
+
+    def merge_values(self, a, b):
+      return a + b
+
+
+  class SimpleParameterBase(GraphParameterBase):
+    def join_boundaries(self, x, _):
+      yield Boundary(
+        tuple(chain.from_iterable( b.boundary for b in x )),
+        reduce(lambda acc, b: acc*b.value, x, 1)
+      )
+
+    def join_edges(self, b, e1, e2):
+      if self.is_compatible(b.boundary, e1, e2):
+        yield b
+
+    def project_and_canonize(self, sel, b):
+      return Boundary(utils.select(sel, b.boundary), b.value)
+
+    def finalize(self, b):
+      b = list(b)
+
+      assert len(b) <= 1
+      if len(b) == 0: return 0
+
+      assert b[0].boundary == tuple()
+      return b[0].value
+
+  class Gadget:
+    def __init__(self, size):
+      self._size = size
+
+    def __hash__(self): return id(self)
+    def __eq__(self, b): return self is b
+
+    def size(self): return self._size
+
+    def is_graph(self):
+      return self._size == 0
+
+    def __str__(self):
+      return "Gadget{%i}" % (self._size,)
+
+    def __repr__(self):
+      return self.__str__()
+
+    _gadget_cache = utils.DynamicLRU()
+
+    @staticmethod
+    def join(gadgets, joins, outs):
+      assert all( isinstance(g, Gadget) for g in gadgets ), \
+        "Not all gadgets of the same size: %s" % gadgets
+
+      info = (tuple(gadgets), tuple(joins), tuple(outs))
+      return Gadget._gadget_cache.get(info, lambda: JoinGadget(*info))
+
+  class FakeGadget(Gadget):
+    def __init__(self, size, value):
+      super().__init__(size)
+      self.value = value
+
+    def eval(self, _):
+      return self.value
+
+
+  class BaseGadget(Gadget):
+    def eval(self, parameter):
+      try:
+        if self is CUBIC_VERTEX:
+          return parameter.CUBIC_VERTEX
+        elif self is FREE_EDGE:
+          return parameter.FREE_EDGE
+      except: pass
+      raise RuntimeError("Unsupported base gadget")
+
+  CUBIC_VERTEX = BaseGadget(3)
+  FREE_EDGE = BaseGadget(2)
+
+
+  class JoinGadget(Gadget):
+    _cache = utils.DynamicLRU()
+
+    def __init__(self, gadgets, joins, outs):
+      super().__init__(len(outs))
+
+      self.offsets = [0]
+      for g in gadgets:
+        self.offsets.append(self.offsets[-1] + g.size())
+
+      joins = tuple(
+        (self.offsets[g1-1] + i1 - 1, self.offsets[g2-1] + i2 - 1)
+        for (g1, i1), (g2, i2) in joins
+      )
+      outs = tuple( self.offsets[g-1] + i - 1 for g, i in outs )
+
+      self._info = (tuple(gadgets), joins, outs)
+
+    def eval(self, parameter):
+      gadgets, joins, outs = self._info
+
+      def do_eval():
+        ret = parameter.eval_join(
+          [ g.eval(parameter) for g in gadgets ],
+          gadgets, joins, outs, self.offsets
+        )
+        if self.is_graph():
+          return parameter.finalize(ret)
+        return ret
+
+      return self._cache.get((parameter, self._info), do_eval)
+
+
+
+
+  def ParametrizedGraphSequence(cls):
+    class GraphSequence(cls, GraphSequenceBase):
+      def __init__(self, *args, **kwargs):
+        GraphSequenceBase.__init__(self)
+        cls.__init__(self, *args, **kwargs)
+    return GraphSequence
+
+  def GraphSequence(cls):
+    return ParametrizedGraphSequence(cls)()
+
+  class GraphSequenceBase:
+    sequence_start = 0
+
+    def __init__(self):
+      self._gadget_cache = utils.DynamicLRU()
+      self._graph_cache = utils.DynamicLRU()
+
+    def _join_gadgets(self, a, b):
+      return Gadget.join([a, b], self.step_join, self.step_out)
+
+    def _next_gadget(self, gadget):
+      return self._join_gadgets(gadget, self.step_gadget)
+
+    def _make_graph(self, gadget):
+      return Gadget.join([ gadget ], self.final_join, [])
+
+    def gadget(self, n):
+      assert type(n) is int and n >= self.sequence_start, "n must be non-negative integer"
+
+      def make_gadget():
+        if n == self.sequence_start: return self.base_gadget
+        return self._next_gadget(self.gadget(n-1))
+
+      return self._gadget_cache.get(n, make_gadget)
+
+    def graph(self, n):
+      return self._graph_cache.get(n, lambda: self._make_graph(self.gadget(n)))
+
+    def stabilize(self, parameter):
+      from sage.all import matrix, QQ
+      boundaries = set()
+      new_boundaries = set( b.boundary for b in self.base_gadget.eval(parameter) )
+
+      i = 0
+      while True:
+        i += 1
+        new_boundaries.update(boundaries)
+        if new_boundaries == boundaries: break
+        boundaries = new_boundaries
+        gadget = FakeGadget(
+          self.base_gadget.size(),
+          [ Boundary(k, 1) for k in boundaries ]
+        )
+        out = self._next_gadget(gadget).eval(parameter)
+        new_boundaries = set( b.boundary for b in out )
+        print("Loop %i done - %i boundaries (%i new)" % \
+          (i, len(new_boundaries), len(new_boundaries.difference(boundaries))))
+
+      sorted_boundaries = sorted(boundaries)
+      rename = { b: i for i, b in enumerate(sorted_boundaries) }
+
+      size = len(sorted_boundaries)
+
+      m = matrix(QQ, size, size, {
+        (rename[row], rename[col]): value
+        for col in sorted_boundaries
+        for row, value in self._next_gadget(FakeGadget(self.base_gadget.size(), [ Boundary(col, 1) ])).eval(parameter)
+      })
+
+      iv = matrix(QQ, size, 1, {
+        (rename[row], 0): v for row, v in self.base_gadget.eval(parameter)
+      })
+
+      fin = matrix(QQ, 1, size, {
+        (0, rename[col]): self._make_graph(FakeGadget(self.base_gadget.size(), [ Boundary(col, 1) ])).eval(parameter)
+        for col in sorted_boundaries
+      })
+
+      return {
+        'variables': sorted_boundaries,
+        'step_matrix': m,
+        'initial_vector': iv,
+        'finalize': fin
+      }
+
+
+_init_()
+
diff --git a/graph_tools/misc.py b/graph_tools/misc.py
new file mode 100644
index 0000000..4631d81
--- /dev/null
+++ b/graph_tools/misc.py
@@ -0,0 +1,114 @@
+def _init_():
+  global sun, count_cdc_naive, graph_to_gadget
+
+  from .definitions import Gadget, CUBIC_VERTEX
+  from .parameters import CircuitDoubleCover
+
+
+  def sun(k):
+    return Gadget.join(
+      [ CUBIC_VERTEX ] * k,
+      [ ((i + 1, 2), ((i+1) % k + 1, 1)) for i in range(k) ],
+      [ (i + 1, 3) for i in range(k) ]
+    )
+
+
+  def count_cdc_naive(G):
+    """Count circuit double covers of G.
+
+    Requires G to be cubic.
+
+    EXAMPLE:
+    >>> from sage.all import *
+    >>> count_cdc_naive(graphs.PetersenGraph())
+    (52, 16)
+    """
+
+    assert(G.degree() == [3]*G.num_verts())
+
+    vert_map = { v: i for i, v in enumerate(G.vertices()) }
+    vert_used = [ 1 ] * len(vert_map)
+
+    joins = []
+    for u, v, _ in G.edges():
+      u = vert_map[u]
+      v = vert_map[v]
+      joins.append(((u+1, vert_used[u]), (v+1, vert_used[v])))
+      vert_used[u] += 1
+      vert_used[v] += 1
+
+    b = Gadget.join([ CUBIC_VERTEX ] * G.num_verts(), joins, [])
+    assert b.is_graph()
+    return (b.eval(CircuitDoubleCover), 2**(G.num_verts() // 2) // 2)
+
+
+  def graph_to_gadget(G, decomposition = None):
+    """Transform graph into a Gadget.
+
+    EXAMPLES:
+    >>> from sage.all import *
+    >>> from .parameters import CircuitDoubleCover, UnderlayingGraph
+    >>> P = graphs.PetersenGraph()
+    >>> graph_to_gadget(P).eval(CircuitDoubleCover)
+    52
+    >>> d = [ list(range(5)), list(range(5, 10)) ]
+    >>> P.is_isomorphic(graph_to_gadget(P, d).eval(UnderlayingGraph))
+    True
+    >>> graph_to_gadget(P, d).eval(CircuitDoubleCover)
+    52
+    """
+    assert(G.degree() == [3]*G.num_verts())
+
+    verts = { v: i for i, v in enumerate(G.vertices()) }
+    edges = [ (verts[u], verts[v]) for u, v, _ in G.edges() ]
+
+    if decomposition is None:
+      decomposition = G.vertices()
+
+    vert_map = { v: [0, 1, 2] for v in range(len(verts)) }
+
+    def process(dec):
+      nonlocal edges
+
+      if not isinstance(dec, list): return (CUBIC_VERTEX, set([ verts[dec] ]))
+
+      gadgets = [ process(x) for x in dec ]
+      vert_to_gadget = dict()
+      gadget_verts = set()
+
+      for i, (_, vs) in enumerate(gadgets):
+        gadget_verts.update(vs)
+        for v in vs:
+          vert_to_gadget[v] = i
+
+      joins = []
+
+      new_edges = []
+      for u, v in edges:
+        if u not in gadget_verts or v not in gadget_verts:
+          new_edges.append((u, v))
+          continue
+
+        joins.append((
+          (vert_to_gadget[u]+1, vert_map[u].pop()+1),
+          (vert_to_gadget[v]+1, vert_map[v].pop()+1)
+        ))
+
+      edges = new_edges
+
+      outs = []
+      for v in gadget_verts:
+        old = vert_map[v]
+        vert_map[v] = list(range(len(outs), len(outs) + len(old)))
+        for i in old:
+          outs.append((vert_to_gadget[v]+1, i+1))
+
+      return (Gadget.join([ g for g, _ in gadgets ], joins, outs), gadget_verts)
+
+    b, _ = process(decomposition)
+    assert b.is_graph()
+    return b
+
+
+_init_()
+
diff --git a/graph_tools/parameters.py b/graph_tools/parameters.py
new file mode 100644
index 0000000..01bfe07
--- /dev/null
+++ b/graph_tools/parameters.py
@@ -0,0 +1,312 @@
+def _init_():
+  global \
+    GroupFlow, \
+    EdgeColoring, VertexColoring, \
+    UnderlayingGraph, VertexCount, \
+    HamiltonianCycle, \
+    CycleDoubleCover, CircuitDoubleCover
+
+
+  from .utils import \
+    OneOfKind, Singleton, \
+    split_on_none, prod, select, powerset, \
+    UnionFind
+
+  from .definitions import \
+    Boundary, \
+    SimpleParameterBase, GraphParameterBase
+
+  from fractions import Fraction
+  from itertools import chain, combinations
+  from functools import reduce
+  from collections import defaultdict
+
+
+  @OneOfKind
+  class GroupFlow(SimpleParameterBase):
+    def __init__(self, group, nowherezero):
+      self.group = group
+      self.nowherezero = nowherezero
+      self.FREE_EDGE = [ Boundary((i, i), 1) for i in group if not nowherezero or i != 0 ]
+      self.CUBIC_VERTEX = [
+        Boundary((i, j, l), 1) for i in group for j in group
+        for l in group if i + j + l == 0 and
+        (not nowherezero or i != 0 and j != 0 and l != 0)
+      ]
+
+    def is_compatible(self, boundary, e_a, e_b):
+      return boundary[e_a] == boundary[e_b]
+
+  @OneOfKind
+  class CycleDoubleCover(SimpleParameterBase):
+    def __init__(self, k):
+      self.k = k
+      self.FREE_EDGE = [ Boundary(((i, j), (i, j)), 1) for i in range(k) for j in range(i+1, k) ]
+      self.CUBIC_VERTEX = [
+        Boundary((tuple(sorted((i, j))), tuple(sorted((i, l))), tuple(sorted((j, l)))), 1)
+        for i in range(k) for j in range(k) for l in range(k)
+        if i != j and j != l and i != l
+      ]
+
+    def is_compatible(self, boundary, e_a, e_b):
+      return boundary[e_a] == boundary[e_b]
+
+
+  @OneOfKind
+  class EdgeColoring(SimpleParameterBase):
+    def __init__(self, k):
+      self.k = k
+      self.FREE_EDGE = [ Boundary((i, i), 1) for i in range(k) ]
+      self.CUBIC_VERTEX = [
+        Boundary((i, j, l), 1) for i in range(k) for j in range(k)
+        for l in range(k) if i != j and j != l and i != l
+      ]
+
+    def is_compatible(self, boundary, e_a, e_b):
+      return boundary[e_a] == boundary[e_b]
+
+
+  @OneOfKind
+  class VertexColoring(SimpleParameterBase):
+    def __init__(self, k):
+      self.k = k
+      self.CUBIC_VERTEX = [ Boundary((i, i, i), 1) for i in range(k) ]
+
+    def is_compatible(self, boundary, e_a, e_b):
+      return boundary[e_a] != boundary[e_b]
+
+
+  @Singleton
+  class UnderlayingGraph(GraphParameterBase):
+    from sage.all import Graph
+
+    CUBIC_VERTEX = [
+      Boundary((0, 0, 0), Graph([[0], []], multiedges=True, loops=True, immutable=True))
+    ]
+
+    def join_boundaries(self, x, _):
+      offsets = [0]
+      for b in x:
+        offsets.append(offsets[-1] + b.value.num_verts())
+
+      gg = (
+        b.value.relabel(lambda v: v + offsets[i], inplace=False)
+        for i, b in enumerate(x)
+      )
+
+      yield Boundary(
+        tuple(chain.from_iterable(
+          [ v + offsets[i] for v in b.boundary ] for i, b in enumerate(x)
+        )),
+        reduce(lambda acc, b: acc.union(b), gg)
+      )
+
+    def join_edges(self, b, e1, e2):
+      g = b.value.copy(immutable=False)
+      g.add_edge(b.boundary[e1], b.boundary[e2])
+      yield Boundary(b.boundary, g)
+
+    def project_and_canonize(self, sel, b):
+      return Boundary(select(sel, b.boundary), b.value)
+
+    def merge_values(self, a, b): assert False
+
+    def finalize(self, b):
+      assert len(b) == 1
+      assert b[0].boundary == tuple()
+      return b[0].value
+
+
+  @Singleton
+  class VertexCount(GraphParameterBase):
+    CUBIC_VERTEX = [ Boundary(None, 1) ]
+    FREE_EDGE = [ Boundary(None, 0) ]
+
+    def join_boundaries(self, x, _):
+      yield Boundary(None, sum( b.value for b in x ))
+
+    def join_edges(self, b, _, __): yield b
+    def project_and_canonize(self, _, b): return b
+    def finalize(self, b): return b[0].value
+
+
+  @Singleton
+  class HamiltonianCycle(GraphParameterBase):
+    FREE_EDGE = [ Boundary((1, 1), 1), Boundary((None, None), 1) ]
+    CUBIC_VERTEX = [
+      Boundary((1, 1, None), 1),
+      Boundary((1, None, 1), 1),
+      Boundary((None, 1, 1), 1)
+    ]
+
+    def join_boundaries(self, x, _):
+      offsets = [ 0 ]
+      try:
+        for b in x:
+          offsets.append(offsets[-1] + max_(b.boundary))
+      except ValueError:
+        return
+
+      def shift(k, off):
+        if k is None: return None
+        return k + off
+
+      ret_b = tuple(chain(*[
+        [ shift(x, offsets[i]) for x in b.boundary ] for i, b in enumerate(x)
+      ]))
+
+      yield Boundary(ret_b, prod( b.value for b in x ))
+
+    def join_edges(self, b, e1, e2):
+      # print("%s %s %s" % (b, e1, e2))
+      c1 = b.boundary[e1]
+      c2 = b.boundary[e2]
+
+      if c1 is None and c2 is None: yield b
+      if c1 is None or c2 is None: return
+
+      if c1 != c2:
+        if c1 > c2: c1, c2 = c2, c1
+        boundary = tuple( v if v != c2 else c1 for v in b.boundary )
+        yield Boundary(boundary, b.value)
+      elif max_(b.boundary) == 1:
+        yield b
+
+    def project_and_canonize(self, p, b):
+      boundary = select(p, b.boundary)
+      rename = { None: None }
+      i = 1
+      for x in boundary:
+        if x not in rename:
+          rename[x] = i
+          i += 1
+      boundary = tuple( rename[x] for x in boundary )
+      return Boundary(boundary, b.value)
+
+    def finalize(self, b):
+      assert len(b) <= 1
+      if len(b) == 0: return 0
+      assert b[0].boundary == tuple()
+      return b[0].value
+
+
+
+  @Singleton
+  class CircuitDoubleCover(GraphParameterBase):
+    FREE_EDGE = [ Boundary(((1,2), (1,2), None), Fraction(1, 2)) ]
+    CUBIC_VERTEX = [ Boundary(((1,2), (2,3), (1,3), None), 1) ]
+
+    def join_boundaries(self, x, offsets):
+      def shift(k, off):
+        if k is None: return None
+        if len(k) == 1: return (k[0] + off,)
+        return (k[0] + off, k[1] + off)
+
+      tmp = tuple(zip(*[
+        split_on_none([ shift(k, offsets[i]) for k in b.boundary ])
+        for i, b in enumerate(x)
+      ]))
+      ret_b = tuple(chain(*tmp[0], (None,), *tmp[1]))
+
+      yield Boundary(ret_b, prod( b.value for b in x ))
+
+    def join_edges(self, b, e1, e2):
+      assert len(b.boundary[e1]) == len(b.boundary[e2])
+
+      def do_join(uf):
+        ret_b = []
+        for i, x in enumerate(b.boundary):
+          if x is None:
+            ret_b.append(None)
+          else:
+            t = tuple( uf.find(i) for i in x )
+            if len(t) == 2 and t[0] == t[1]: break
+            ret_b.append(t)
+        else:
+          yield Boundary(tuple(ret_b), b.value)
+
+      if len(b.boundary[e1]) == 1:
+        uf = UnionFind(len(b.boundary))
+        uf.union(b.boundary[e1][0], b.boundary[e2][0])
+        yield from do_join(uf)
+        return
+
+      c1, c2 = b.boundary[e1]
+      c3, c4 = b.boundary[e2]
+
+      for _ in range(2):
+        uf = UnionFind(len(b.boundary))
+        uf.union(c1, c3)
+        uf.union(c2, c4)
+        yield from do_join(uf)
+        c3, c4 = c4, c3
+
+    def project_and_canonize(self, p, bo):
+      x = bo.boundary
+      l, r = split_on_none(x)
+      r = set(r).union(l)
+      l = select(p, l)
+      r = r.difference(l)
+      r = { x for x in r if len(x) == 2 }
+
+      rename = defaultdict(list)
+
+      for i, x in enumerate(l):
+        for a in x: rename[a].append(i)
+
+      rename = sorted( (v, k) for k, v in rename.items() )
+      rename = { x[1]: i + 1 for i, x in enumerate(rename) }
+
+      l = [ tuple(sorted([ rename[a] for a in x ])) for x in l ]
+      r = { tuple(sorted([rename[a], rename[b]])) for a, b in r \
+        if a in rename and b in rename }
+
+      r.difference_update(set(l))
+
+      return Boundary(tuple(l + [None] + sorted(r)), bo.value)
+
+    def finalize(self, b):
+      assert len(b) <= 1
+      if len(b) == 0: return 0
+      assert b[0].boundary == (None,)
+      return b[0].value
+
+    def enumerate_boundaries(self, n):
+      if n == 0: yield (None,)
+      if n <= 1: return
+
+      used = set([1, 2])
+      first_unused = 3
+      rng = list(range(n))
+
+      def add_internal(b):
+        X = set(combinations(range(1, n + 1), 2)).difference(set(b))
+        b = tuple(b + [None])
+        for x in powerset(X):
+          ret = b + tuple(sorted(x))
+          if self.project_and_canonize(rng, Boundary(ret, 0)).boundary == ret:
+            yield ret
+
+      def add_edge(b, used, fu):
+        if len(b) == n:
+          if len(used) == 0:
+            yield from add_internal(b)
+          return
+
+        s = set(used)
+        if fu <= n: s.add(fu)
+        if fu < n: s.add(fu + 1)
+        for x, y in combinations(sorted(s), 2):
+          if fu + 1 in [x, y] and fu not in [x, y]: continue
+          b_ = b + [(x, y)]
+          s_ = used.symmetric_difference(set([x, y]))
+          fu_ = fu
+          if fu + 1 in [x, y]: fu_ += 2
+          elif fu in [x, y]: fu_ += 1
+          yield from add_edge(b_, s_, fu_)
+
+      yield from add_edge([(1,2)], used, first_unused)
+
+
+_init_()
+
diff --git a/graph_tools/sequences.py b/graph_tools/sequences.py
new file mode 100644
index 0000000..e63a15f
--- /dev/null
+++ b/graph_tools/sequences.py
@@ -0,0 +1,189 @@
+def _init_():
+  global \
+    FlowerSnark, FlowerSnarkAlt, \
+    Necklace, Flower, \
+    GeneralizedPetersen, Petersen
+
+  from .definitions import \
+    Boundary, \
+    Gadget, FakeGadget, CUBIC_VERTEX, FREE_EDGE, \
+    GraphSequence, ParametrizedGraphSequence
+
+
+  @GraphSequence
+  class FlowerSnark:
+    sequence_start = 2
+
+    step_join = [
+      ((1, 2), (2, 1)),
+      ((1, 4), (2, 3)),
+      ((1, 6), (2, 5))
+    ]
+    step_out = [
+      (1, 1), (2, 2),
+      (1, 3), (2, 4),
+      (1, 5), (2, 6)
+    ]
+    final_join = [
+      ((1, 2), (1, 1)),
+      ((1, 4), (1, 3)),
+      ((1, 6), (1, 5)),
+    ]
+
+    BASIC_GADGET = Gadget.join(
+      [ CUBIC_VERTEX ] * 4,
+      [
+        ((1, 3), (2, 1)),
+        ((2, 2), (3, 3)),
+        ((2, 3), (4, 1))
+      ],
+      [
+        (1, 1), (1, 2),
+        (3, 1), (3, 2),
+        (4, 2), (4, 3)
+      ]
+    )
+    DOUBLE_GADGET = Gadget.join([ BASIC_GADGET ]*2, step_join, step_out)
+    SPECIAL_GADGET = Gadget.join(
+      [ BASIC_GADGET ] * 2,
+      [
+        ((1, 2), (2, 3)),
+        ((1, 4), (2, 1)),
+        ((1, 6), (2, 5))
+      ],
+      step_out
+    )
+    base_gadget = SPECIAL_GADGET
+    step_gadget = BASIC_GADGET
+
+  @GraphSequence
+  class FlowerSnarkAlt:
+    sequence_start = 1
+
+    step_join = [
+      ((1, 2), (2, 1)),
+      ((1, 4), (2, 3)),
+      ((1, 6), (2, 5))
+    ]
+    step_out = [
+      (1, 1), (2, 2),
+      (1, 3), (2, 4),
+      (1, 5), (2, 6)
+    ]
+    final_join = [
+      ((1, 2), (1, 1)),
+      ((1, 4), (1, 3)),
+      ((1, 6), (1, 5)),
+    ]
+
+    BASIC_GADGET = Gadget.join(
+      [ CUBIC_VERTEX ] * 4,
+      [
+        ((1, 3), (2, 1)),
+        ((2, 2), (3, 3)),
+        ((2, 3), (4, 1))
+      ],
+      [
+        (3, 1), (1, 2),
+        (1, 1), (3, 2),
+        (4, 2), (4, 3)
+      ]
+    )
+    base_gadget = BASIC_GADGET
+    step_gadget = BASIC_GADGET
+
+  @GraphSequence
+  class Necklace:
+    sequence_start = 1
+
+    step_join = [ ((1,2), (2, 1)) ]
+    step_out = [ (1,1), (2,2) ]
+    final_join = [ ((1,1), (1,2)) ]
+
+    D = Gadget.join(
+      [ CUBIC_VERTEX ] * 4,
+      [
+        ((1, 1), (2, 1)),
+        ((1, 2), (3, 1)),
+        ((1, 3), (4, 1)),
+        ((2, 2), (3, 2)),
+        ((2, 3), (4, 2))
+      ],
+      [(3,3), (4,3)]
+    )
+    base_gadget = D
+    step_gadget = D
+
+
+  @GraphSequence
+  class Flower:
+    sequence_start = 2
+
+    SCV = FakeGadget(3, [ Boundary(((1,), (2,), (1, 2), None), 1) ])
+
+    D = Gadget.join(
+      [ SCV, CUBIC_VERTEX ],
+      [ ((1,3), (2,3)) ],
+      [ (1,1), (2,1), (2,2), (1,2) ]
+    )
+
+    step_join = [ ((1,3), (2,2)), ((1,4), (2,1)) ]
+    step_out = [ (1,1), (1,2), (2,3), (2,4) ]
+    final_join = [ ((1,1), (1,4)), ((1,2), (1,3)) ]
+
+    base_gadget = Gadget.join([D, D], step_join, step_out)
+    step_gadget = D
+
+
+  @ParametrizedGraphSequence
+  class GeneralizedPetersen:
+    sequence_start = 1
+
+    def __init__(self, k):
+      assert k >= 1
+
+      self.k = k
+      self.GADGET = Gadget.join(
+        [ CUBIC_VERTEX ]*2 + [ FREE_EDGE ]*(k-1),
+        [ ((1,1), (2,1)) ],
+        [ (2,2) ] +
+        [ (e + 3, i + 1) for e in range(k-1) for i in range(2) ] +
+        [ (2,3), (1,2), (1,3) ]
+      )
+      self.step_gadget = self.GADGET
+      self.base_gadget = self.GADGET
+
+      self.step_join = [ ((1,i+1), (2,i)) for i in range(1, 2*k+3, 2) ]
+      self.step_out = [ (2 - (i % 2), i) for i in range(1, 2*k+3) ]
+      self.final_join = [ ((1,i+1), (1,i)) for i in range(1, 2*k+3, 2) ]
+
+    def graph(self, n):
+      assert n > self.k
+      return super().graph(n)
+
+
+  Petersen = Gadget.join(
+    [ CUBIC_VERTEX ] * 10,
+    [
+      ((1,1), (5,3)),
+      ((1,2), (6,1)),
+      ((1,3), (2,1)),
+      ((2,2), (7,1)),
+      ((2,3), (3,1)),
+      ((3,2), (8,1)),
+      ((3,3), (4,1)),
+      ((4,2), (9,1)),
+      ((4,3), (5,1)),
+      ((5,2), (10,1)),
+      ((6,2), (8,3)),
+      ((6,3), (9,2)),
+      ((7,2), (10,2)),
+      ((7,3), (9,3)),
+      ((8,2), (10,3))
+    ],
+    []
+  )
+
+
+_init_()
+
diff --git a/graph_tools/tests.py b/graph_tools/tests.py
new file mode 100644
index 0000000..0a75806
--- /dev/null
+++ b/graph_tools/tests.py
@@ -0,0 +1,21 @@
+"""Helper to run tests.
+
+Doctest does not find tests in docstrings of decorated classes
+so we put them here instead.
+
+
+>>> from .all import *
+>>> from sage.all import *
+
+>>> graphs.CompleteGraph(4).is_isomorphic(Necklace.graph(1).eval(UnderlayingGraph))
+
+>>> all( Necklace.graph(i).eval(CircuitDoubleCover) == 2 * 4**(i-1) for i in range(1, 10) )
+True
+>>> GeneralizedPetersen(2).graph(5).eval(CircuitDoubleCover) == 52
+True
+>>> FlowerSnark.graph(3).eval(CircuitDoubleCover)
+104
+
+
+"""
+
diff --git a/graph_tools/utils.py b/graph_tools/utils.py
new file mode 100644
index 0000000..16361cf
--- /dev/null
+++ b/graph_tools/utils.py
@@ -0,0 +1,146 @@
+def _init_():
+  global \
+    Singleton, OneOfKind, \
+    UnionFind, DynamicLRU, \
+    max_, select, prod, split_on_none, powerset
+
+  from itertools import chain, combinations
+  from functools import reduce
+  from collections import defaultdict
+  import operator
+  import weakref
+
+  def Singleton(cls):
+    """Decorator which replaces the class with object of it."""
+    return cls()
+
+  def OneOfKind(cls):
+    """Decorator which ensures that there exisits at most one
+    object of given class for every tuple of constuctor
+    arguments."""
+
+    cache = DynamicLRU()
+
+    def call(*args, **kwargs):
+      key = (tuple(args), tuple(sorted(kwargs.items())))
+      return cache.get(key, lambda: cls(*args, **kwargs))
+
+    return call
+
+  class UnionFind:
+    """Simple union-find datastructure.
+
+    It is implemented with pointers. No path compression
+    or ranks are implemented.
+
+    EXAMPLES:
+    >>> uf = UnionFind(5)
+    >>> any( uf.same(i, j) for i in range(5) for j in range(i) )
+    False
+    >>> uf.union(3, 2)
+    >>> any( uf.same(i, j) for i in range(5) for j in range(i) )
+    True
+    >>> uf.same(2, 3)
+    True
+    >>> uf.same(1, 2)
+    False
+    >>> uf.union(2, 0)
+    >>> uf.same(0, 3)
+    True
+    """
+
+    def __init__(self, size):
+      self.x = [ None ] * size
+
+    def union(self, a, b):
+      a = self.find(a)
+      b = self.find(b)
+      if a != b: self.x[b] = a
+
+    def find(self, a):
+      while self.x[a] is not None:
+        a = self.x[a]
+
+      return a
+
+    def same(self, a, b):
+      return self.find(a) == self.find(b)
+
+
+  def max_(l):
+    """Version of `max` which ignores `None` in the list of arguments.
+
+    EXPAMPLES:
+    >>> max_([1, 5, 2, 3])
+    5
+    >>> max_([0, -14, 2, 7, None, 9])
+    9
+    """
+    return max( x for x in l if x is not None )
+
+  def select(pat, l):
+    """Shuffle `l` according to pattern `pat`.
+
+    Returns list-like object of the same type as `pat`
+    and of the same length as `pat` where on `i`-th
+    position is element `l[pat[i]]`.
+
+    EXAMPLES:
+    >>> select([5, 3, 1], [1, 2, 3, 4, 5, 6])
+    [6, 4, 2]
+    >>> select((5, 3, 1), [1, 2, 3, 4, 5, 6])
+    (6, 4, 2)
+    """
+
+    return type(pat)( l[i] for i in pat )
+
+
+  def prod(iterable):
+    """Calculate product of elements of `iterable`."""
+    return reduce(operator.mul, iterable, 1)
+
+
+  class DynamicLRU:
+    def __init__(self, capacity = 100):
+      self._capacity = capacity
+      self._data = dict()
+      self._time = 0
+
+    def get(self, key, make_value):
+      self._time += 1
+
+      if key not in self._data:
+        self._data[key] = (make_value(), self._time)
+
+      return self._data[key][0]
+
+
+  def split_on_none(t):
+    """Split list (or tuple) at first `None`.
+
+    EXAMPLE:
+    >>> split_on_none([1, 2, 3, None, 4, 5])
+    ([1, 2, 3], [4, 5])
+    >>> split_on_none([1, 2])
+    ([1, 2], [])
+    """
+
+    i = 0
+    while i < len(t) and t[i] is not None: i += 1
+    return (t[:i], t[i+1:])
+
+
+  def powerset(iterable):
+      """Enumerate all subsets of `iterable`.
+
+      EXAMPLE:
+      >>> list(powerset([1,2,3]))
+      [(), (1,), (2,), (3,), (1, 2), (1, 3), (2, 3), (1, 2, 3)]
+      """
+
+      s = list(iterable)
+      return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))
+
+
+_init_()
+
diff --git a/test_cycle.py b/test_cycle.py
index 9bcd135..6362071 100644
--- a/test_cycle.py
+++ b/test_cycle.py
@@ -1,5 +1,8 @@
-from flower_snarks import *
-from sage.all import *
+from graph_tools.definitions import *
+from graph_tools.parameters import CircuitDoubleCover, VertexCount
+from graph_tools.misc import sun
+
+from sage.all import MixedIntegerLinearProgram
 from itertools import permutations
 
 def solve_gadget(g, alt_gadgets):
diff --git a/utils.py b/utils.py
deleted file mode 100644
index ea9a370..0000000
--- a/utils.py
+++ /dev/null
@@ -1,146 +0,0 @@
-from itertools import *
-from functools import reduce
-import operator
-from collections import defaultdict
-import weakref
-
-def Singleton(cls):
-  """Decorator which replaces the class with object of it."""
-  return cls()
-
-def OneOfKind(cls):
-  """Decorator which ensures that there exisits at most one
-  object of given class for every tuple of constuctor
-  arguments."""
-
-  cache = DynamicLRU()
-
-  def call(*args, **kwargs):
-    key = (tuple(args), tuple(sorted(kwargs.items())))
-    return cache.get(key, lambda: cls(*args, **kwargs))
-
-  return call
-
-class UnionFind:
-  """Simple union-find datastructure.
-
-  It is implemented with pointers. No path compression
-  or ranks are implemented.
-
-  EXAMPLES:
-  >>> uf = UnionFind(5)
-  >>> any( uf.same(i, j) for i in range(5) for j in range(i) )
-  False
-  >>> uf.union(3, 2)
-  >>> any( uf.same(i, j) for i in range(5) for j in range(i) )
-  True
-  >>> uf.same(2, 3)
-  True
-  >>> uf.same(1, 2)
-  False
-  >>> uf.union(2, 0)
-  >>> uf.same(0, 3)
-  True
-  """
-
-  def __init__(self, size):
-    self.x = [ None ] * size
-
-  def union(self, a, b):
-    a = self.find(a)
-    b = self.find(b)
-    if a != b: self.x[b] = a
-
-  def find(self, a):
-    while self.x[a] is not None:
-      a = self.x[a]
-
-    return a
-
-  def same(self, a, b):
-    return self.find(a) == self.find(b)
-
-
-def max_(l):
-  """Version of `max` which ignores `None` in the list of arguments.
-
-  EXPAMPLES:
-  >>> max_([1, 5, 2, 3])
-  5
-  >>> max_([0, -14, 2, 7, None, 9])
-  9
-  """
-  return max( x for x in l if x is not None )
-
-def select(pat, l):
-  """Shuffle `l` according to pattern `pat`.
-
-  Returns list-like object of the same type as `pat`
-  and of the same length as `pat` where on `i`-th
-  position is element `l[pat[i]]`.
-
-  EXAMPLES:
-  >>> select([5, 3, 1], [1, 2, 3, 4, 5, 6])
-  [6, 4, 2]
-  >>> select((5, 3, 1), [1, 2, 3, 4, 5, 6])
-  (6, 4, 2)
-  """
-
-  return type(pat)( l[i] for i in pat )
-
-
-def prod(iterable):
-  """Calculate product of elements of `iterable`."""
-  return reduce(operator.mul, iterable, 1)
-
-
-class DynamicLRU:
-  def __init__(self, capacity = 100):
-    self._capacity = capacity
-    self._data = dict()
-    self._time = 0
-
-  def get(self, key, make_value):
-    self._time += 1
-
-    if key not in self._data:
-      self._data[key] = (make_value(), self._time)
-
-    return self._data[key][0]
-
-
-def split_on_none(t):
-  """Split list (or tuple) at first `None`.
-
-  EXAMPLE:
-  >>> split_on_none([1, 2, 3, None, 4, 5])
-  ([1, 2, 3], [4, 5])
-  >>> split_on_none([1, 2])
-  ([1, 2], [])
-  """
-
-  i = 0
-  while i < len(t) and t[i] is not None: i += 1
-  return (t[:i], t[i+1:])
-
-
-def powerset(iterable):
-    """Enumerate all subsets of `iterable`.
-
-    EXAMPLE:
-    >>> list(powerset([1,2,3]))
-    [(), (1,), (2,), (3,), (1, 2), (1, 3), (2, 3), (1, 2, 3)]
-    """
-
-    s = list(iterable)
-    return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))
-
-
-if __name__ == "__main__":
-  import doctest
-  (f, t) = doctest.testmod()
-  print("%s: %i tests of %i failed." % (__file__, f, t))
-  if f > 0:
-    from sys import exit
-    exit(1)
-
-- 
GitLab