There is a generic tool for solving such problems without having to think too much: Integer Linear Programming.

For each edge $uv$ of your graph, you define two binary variables (with values in {0, 1} ): one named $c(uv)$ telling that a cycle from you spanning configuration that is passing through the edge $uv$, and one names $e(uv)$ telling whether the edge the edge $uv$ is selected as an edge of the spanning configuration.

Then, you add the following constraint: for each vertex $v$, there is either exactly one edge adjacent to it that is selected to be an edge of the spanning configuration, or there are exactly two edges adjacent to it that are selected to be part of a cycle of the spanning configuration.

This constraints translates into the following linear equalities: $\forall v \in V(G), \sum_{uv \in E(G)} 2*e(uv) + c(uv) =1$ (E)

Hence, the set of solutions you want can be seen as the set of integer points of the polytope defined by the equation (E) together with the equations $0\leq e(uv)\leq1$ and $0\leq c(uv)\leq1$ for all edge $uv \in E(G)$.

The way to define this polytope in Sage is done by using the `polyhedron`

method of the `MixedIntegerLinearProgram`

class.

Here is the code of the function (some explanations will follow):

```
def spanning_elementary_subgraphs(G):
p = MixedIntegerLinearProgram(solver='ppl')
e = p.new_variable(binary=True)
c = p.new_variable(binary=True)
# those are useless constraints, only to control the order in which the variables are created, see discussion below
for edge in G.edges(labels=False):
p.add_constraint(e[edge] >= 0)
for edge in G.edges(labels=False):
p.add_constraint(c[edge] >= 0)
for v in G.vertices():
p.add_constraint(sum(2*e[edge] + c[edge] for edge in G.edges_incident(v, labels=False)) == 2)
P = p.polyhedron(backend='normaliz')
for i,L in enumerate(P.integral_points()):
edges = []
cycle_edges = []
for j,edge in enumerate(G.edges(labels=False)):
if L[j] == 1:
edges.append(edge)
if L[j+G.num_edges()] == 1:
cycle_edges.append(edge)
H = Graph(cycle_edges, format='list_of_edges')
yield {'edges': edges, 'cycles': H.connected_components()}
```

This function returns an iterator over all solutions. At each iteration, you get a dictionary

Here is it working on your example:

```
sage: G = graphs.EmptyGraph()
sage: G.add_edges([(1,2),(2,3),(3,4),(4,5),(5,1),(6,5),(6,8),(8,9),(7,9),(7,6),(7,10),(10,11),(10,12),(11,12)])
sage: s = spanning_elementary_subgraphs(G)
sage: list(s)
[{'cycles': [[1, 2, 3, 4, 5], [6, 7, 8, 9], [10, 11, 12]], 'edges': []},
{'cycles': [[1, 2, 3, 4, 5], [10, 11, 12]], 'edges': [(6, 8), (7, 9)]},
{'cycles': [[1, 2, 3, 4, 5], [10, 11, 12]], 'edges': [(6, 7), (8, 9)]},
{'cycles': [], 'edges': [(1, 2), (3, 4), (5, 6), (7, 10), (8, 9), (11, 12)]}]
```

So, you can see that there are exactly 4 possible spanning configurations.

Now, some more explanations about the code.

First, when we construct the polyhedron, we have to use the `normaliz`

backend or the computation will be too slow. To be able to use it, you have to install the `pynormaliz`

optional package, by typing from a terminal:

```
sage -i pynormaliz
```

Then, there is a total of $2*|E(G)|$ variables: for each edge $uv$, there is $e(uv)$ and $c(uv)$. Each variable will be an element of the basis of the space where the polyhedron lives.

The problem is that, in the current Sage user interface, there is no easy way to see order of the variables that are used to define that basis.

The variables are created on the fly when used for the first time. Since there is a loop over the vertices $v$ and then a loop over the edges $uv$ adjacent to the current vertex $v$ and then a constraint involving both $e(uv)$ and $c(uv)$, it is hard to guess the order of the basis.

Hence, we added some useless constraints in two flat loops in whch the edges appear in a simple-to-handle order.

Note that it trick has nothing to do with interesting algorithmics or mathematics, it is only an annoyance of the current Sage user interface and should be solved by Sage developers (a ticket will follow).