diff --git a/phi_34_gen.py b/phi_34_gen.py
index e5c2402..db5bdf3 100644
--- a/phi_34_gen.py
+++ b/phi_34_gen.py
@@ -57,22 +57,23 @@ def gen_graphs(L, m, cntd, edge2cntd, vtx2cntd, notadpoles):
                 g_bulk, frozenset(range(bulk_n)),
                 L, m, notadpoles))
 
-            # do we need frozenset here ??
-            unlabeled_graphs = frozenset(g.unlabeled_graph
-                                         for g in labeled_graphs)
+            # store graphs according to number of edges
+            stock = {k: set() for k in range(min_edges, max_edges + 1)}
 
-            for g in unlabeled_graphs:
+            for nedges, g in labeled_graphs:
                 if vtx2cntd and not g.is_vtx_2_connected:
                     continue
-                elif edge2cntd and not g.is_edge_2_connected:
+                if edge2cntd and not g.is_edge_2_connected:
                     continue
-                elif cntd and not g.is_connected:
-                    continue
-
+                # useless, as we already asked nauty for connected graphs
+                # elif cntd and not g.is_connected:
+                #    continue
                 if not vtx2cntd and notadpoles and g.is_tadpole:
                     continue
-
-                yield g
+                new_g = g.unlabeled_graph
+                if new_g not in stock[nedges]:
+                    stock[nedges].add(new_g)
+                    yield new_g
 
 
 def gen_from_bulk_g(g, vtcs_set, L, m, notadpoles):
@@ -90,7 +91,8 @@ def gen_ext_vtcs():
             yield leg_candidates
             return
 
-        if not notadpoles:  # Extra legs can still be closed with self-loops!
+        if not notadpoles:
+            # Extra legs can still be closed by adding self-loops!
             for ext_vtcs in itertools.combinations(leg_candidates, m):
                 yield frozenset(ext_vtcs)
 
@@ -101,9 +103,94 @@ def gen_ext_vtcs():
         if any(d < 0 for d in degree_defs):
             continue
 
+        # adding loops to inner vertices of valence 1 or 2 (default 3 or 2)
         selfloop_edges = [(v, v) for v, d in zip(int_vtcs, degree_defs)
                           for i in range(d // 2)]
 
         edges = g.edges + selfloop_edges
         if len(edges) - len(vtcs_set) == L - 1:
-            yield Graph(edges)
+            yield (len(edges), Graph(edges))
+
+
+# sage versions
+
+
+def gen_graphs_sage(L, m, cntd, edge2cntd, vtx2cntd, notadpoles):
+    """Generate phi^3 + phi^4 graphs with the desired parameters and
+        properties.
+        L: Loop number
+        m: Ext. legs"""
+
+    cntd = cntd | edge2cntd | vtx2cntd
+    edge2cntd = edge2cntd | vtx2cntd
+    notadpoles = notadpoles | vtx2cntd
+
+    min_n, max_n, min_edges, max_edges = calc_gen_params(L, m,
+                                                         cntd, notadpoles)
+
+    if max_n <= 0:
+        return
+
+    for bulk_n in range(min_n, max_n + 1):
+        for g_bulk in nauty_ctrl.gen_nauty_graphs_sage(
+                bulk_n, cntd, 4, min_edges, max_edges):
+            labeled_graphs = (g for g in gen_from_bulk_g_sage(
+                g_bulk, frozenset(range(bulk_n)),
+                L, m, notadpoles))
+
+            # store graphs according to number of edges
+            stock = {k: set() for k in range(min_edges, max_edges + 1)}
+
+            for nedges, g in labeled_graphs:
+                if vtx2cntd and not g.is_biconnected():
+                    continue
+
+                # TODO
+                # if edge2cntd and g.edge_connectivity(value_only=True) < 2:
+                #    continue
+
+                # useless, as we already asked nauty for connected graphs
+                # elif cntd and not g.is_connected:
+                #    continue
+
+                # TODO:
+                # if not vtx2cntd and notadpoles and g.is_tadpole:
+                #    continue
+
+                new_g = g.canonical_label().copy(immutable=True)
+                if new_g not in stock[nedges]:
+                    stock[nedges].add(new_g)
+                    yield new_g
+
+
+def gen_from_bulk_g_sage(g, vtcs_set, L, m, notadpoles):
+    leg_candidates = frozenset(v for v in g if g.degree(v) == 1)
+
+    def gen_ext_vtcs():
+        if len(leg_candidates) < m:
+            return
+
+        if len(leg_candidates) == m:
+            yield leg_candidates
+            return
+
+        if not notadpoles:
+            # Extra legs can still be closed by adding self-loops!
+            for ext_vtcs in itertools.combinations(leg_candidates, m):
+                yield frozenset(ext_vtcs)
+
+    for ext_vtcs in gen_ext_vtcs():
+
+        degree_defs = [(v, 4 - g.degree(v)) for v in g if v not in ext_vtcs]
+        if any(d < 0 for d, _ in degree_defs):
+            continue
+
+        # adding loops to inner vertices of valence 1 or 2 (default 3 or 2)
+        selfloop_edges = [(v, v) for v, d in degree_defs
+                          for i in range(d // 2)]
+
+        num_edges = g.num_edges() + len(selfloop_edges)
+        if num_edges - g.num_verts() == L - 1:
+            new_g = g.copy()
+            new_g.add_edges(selfloop_edges)
+            yield (num_edges, new_g)