Number of disjoint hamiltonian cycles in a graph

edit

This method will return the maximum number of edge disjoint Hamiltonian cycles. This is certainly not the best possible formulation.

def disjoint_hamiltonian_cycles(g, solver=None, verbose=0,
                                *, integrality_tolerance=1e-3):
    """
    """
    from sage.numerical.mip import MixedIntegerLinearProgram
    from sage.numerical.mip import MIPSolverException
    from sage.graphs.graph import Graph

    UB = min(g.degree()) // 2
    if not UB:
        raise EmptySetError("the given graph is not Hamiltonian")

    eps = 1 / (2*Integer(g.order()))
    x = next(g.vertex_iterator())

    while True:
        # Search for UB edge disjoint Hamiltonian cycles
        p = MixedIntegerLinearProgram(solver=solver)

        f = p.new_variable(binary=True)
        r = p.new_variable(nonnegative=True)

        # An edge belongs to a unique cycle
        for e in g.edge_iterator(labels=False):
            p.add_constraint(p.sum(f[frozenset(e), i] for i in range(UB)) <= 1)

        # All the vertices of cycle i have degree 2
        for i in range(UB):
            for v in g:
                p.add_constraint(p.sum(f[frozenset((u, v)), i] for u in g.neighbor_iterator(v)), min=2, max=2)

        # Connectivity constraints based on the maximum average degree
        for i in range(UB):
            # r is greater than f
            for u, v in g.edge_iterator(labels=None):
                p.add_constraint(r[u, v, i] + r[v, u, i] - f[frozenset((u, v)), i], min=0)

            # no cycle which does not contain x
            for v in g:
                if v != x:
                    p.add_constraint(p.sum(r[u, v, i] for u in g.neighbor_iterator(v)), max=1 - eps)

        try:
            p.solve(log=verbose)
            f_val = p.get_values(f, convert=bool, tolerance=integrality_tolerance)
            tsp = [Graph() for i in range(UB)]
            for i in range(UB):
                tsp[i].add_vertices(g.vertex_iterator())
                tsp[i].set_pos(g.get_pos())
            for (e, i), b in f_val.items():
                if b:
                    tsp[i].add_edge(e)
            return tsp

        except MIPSolverException:
            UB -= 1
            if not UB:
                raise EmptySetError("the given graph is not Hamiltonian")

You get:

sage: G = graphs.CompleteGraph(5)
sage: T = disjoint_hamiltonian_cycles(G)
sage: T
[Graph on 5 vertices, Graph on 5 vertices]
sage: [t.edges(labels=False, sort=True) for t in T]
[[(0, 1), (0, 3), (1, 4), (2, 3), (2, 4)],
 [(0, 2), (0, 4), (1, 2), (1, 3), (3, 4)]]
sage: G = graphs.CompleteGraph(7)
sage: T = disjoint_hamiltonian_cycles(G)
sage: [t.edges(labels=False, sort=True) for t in T]
[[(0, 1), (0, 4), (1, 6), (2, 3), (2, 5), (3, 6), (4, 5)],
 [(0, 2), (0, 5), (1, 3), (1, 4), (2, 6), (3, 5), (4, 6)],
 [(0, 3), (0, 6), (1, 2), (1, 5), (2, 4), (3, 4), (5, 6)]]