add method to compute the longest (induced) cycle in a (di)graph

[dcoudert]

This PR adds a method to compute the longest (induced) cycle in a (di)graph. The method can also consider weighted cases.

This answers a request from https://ask.sagemath.org/question/75124/how-to-find-a-longest-cycle-and-a-longest-induced-cycle-in-a-graph/

📝 Checklist

• The title is concise, informative, and self-explanatory.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation accordingly.

⌛ Dependencies

[tscrim]

Some doc nitpicks.

I am also guessing you cannot say anything about how the solution improves for setting h in line 8141. In particular, it makes sense to build a new graph each iteration rather than manipulate a single (mutable) graph?

[dcoudert]

I am also guessing you cannot say anything about how the solution improves for setting h in line 8141. In particular, it makes sense to build a new graph each iteration rather than manipulate a single (mutable) graph?

Indeed, the set of selected edges can be very different from one iteration to the next. Furthermore, the time needed to build the graph is expected to be small compared to the time needed to solve the ILP.

[tscrim]

Thanks. That is what I was thinking. Positive review.

[dcoudert]

Thank you for the review.

[github-actions]

Documentation preview for this PR (built with commit 8b71477; changes) is ready! 🎉

[lichengzhang1]

Hi, can you try this example?

Z=Graph({0:[4,6,10,11],1:[5,7,10,11],2:[8,9,10,11],3:[8,9,11],4:[6,10,11],5:[7,10,11],6:[10,11],
         7:[10],8:[10],9:[11]})

Currently, my SageMath version is 10.2, and I excerpted your code as follows, but it seems that the longest cycle and the longest induced cycle are both having issues for the graph.

def longest_cycle(self, induced=False, use_edge_labels=False,
                      solver=None, verbose=0, *, integrality_tolerance=0.001):
        self._scream_if_not_simple()
        G = self
        st = f" from {G.name()}" if G.name() else ""
        name = f"longest{' induced' if induced else ''} cycle{st}"
        if use_edge_labels:
            def weight(e):
                return 1 if (len(e) < 3 or e[2] is None) else e[2]

            def total_weight(gg):
                return sum(weight(e) for e in gg.edge_iterator())
        else:
            def weight(e):
                return 1

            def total_weight(gg):
                return gg.order()

        directed = G.is_directed()
        immutable = G.is_immutable()
        if directed:
            from sage.graphs.digraph import DiGraph as MyGraph
            blocks = G.strongly_connected_components()
        else:
            from sage.graphs.graph import Graph as MyGraph
            blocks = G.blocks_and_cut_vertices()[0]
        if len(blocks) > 1:
            best = MyGraph(name=name, immutable=immutable)
            best_w = 0
            for block in blocks:
                if induced and len(block) < 4:
                    continue
                h = G.subgraph(vertices=block)
                C = h.longest_cycle(induced=induced,
                                    use_edge_labels=use_edge_labels,
                                    solver=solver, verbose=verbose,
                                    integrality_tolerance=integrality_tolerance)
                if total_weight(C) > best_w:
                    best = C
                    best_w = total_weight(C)
            return (best_w, best) if use_edge_labels else best       
        if ((induced and G.order() < 4) or
            (not induced and ((directed and G.order() < 2) or
                              (not directed and G.order() < 3)))):
            if use_edge_labels:
                return 0, MyGraph(name=name, immutable=immutable)
            return MyGraph(name=name, immutable=immutable)
        if (not induced and ((directed and G.order() == 2) or
                             (not directed and G.order() == 3))):
            answer = G.copy()
            answer.name(name)
            if use_edge_labels:
                return total_weight(answer), answer
            return answer
        if directed:
            def F(e):
                return e[:2]
        else:
            def F(e):
                return frozenset(e[:2])

        from sage.numerical.mip import MixedIntegerLinearProgram
        from sage.numerical.mip import MIPSolverException

        p = MixedIntegerLinearProgram(maximization=True,
                                      solver=solver,
                                      constraint_generation=True)
        vertex = p.new_variable(binary=True)
        edge = p.new_variable(binary=True)
        p.set_objective(p.sum(weight(e) * edge[F(e)] for e in G.edge_iterator()))
        p.add_constraint(p.sum(edge[F(e)] for e in G.edge_iterator())
                         == p.sum(vertex[u] for u in G))
        if directed:
            for u in G:
                p.add_constraint(p.sum(edge[F(e)] for e in G.outgoing_edge_iterator(u))
                                 <= vertex[u])
                p.add_constraint(p.sum(edge[F(e)] for e in G.incoming_edge_iterator(u))
                                 <= vertex[u])
        else:
            for u in G:
                p.add_constraint(p.sum(edge[F(e)] for e in G.edge_iterator(u))
                                 <= 2 * vertex[u])
        if induced:
            for e in G.edge_iterator():
                f = F(e)
                u, v = f
                p.add_constraint(edge[f] <= vertex[u])
                p.add_constraint(edge[f] <= vertex[v])
                p.add_constraint(vertex[u] + vertex[v] <= edge[f] + 1)
            p.add_constraint(p.sum(vertex[u] for u in G), min=4)

        best = MyGraph(name=name, immutable=immutable)
        best_w = 0
        while True:
            try:
                p.solve(log=verbose)
            except MIPSolverException:
                break
            b_val = p.get_values(edge, convert=bool, tolerance=integrality_tolerance)
            edges = (e for e in G.edge_iterator() if b_val[F(e)])
            h = MyGraph(edges, format='list_of_edges', name=name, immutable=immutable)
            if not h:
                break
            if directed:
                cc = h.strongly_connected_components()
            else:
                cc = h.connected_components(sort=False)
            if len(cc) == 1:
                if total_weight(h) > best_w:
                    best = h
                    best_w = total_weight(best)
                break
            for c in cc:
                if not (induced and len(c) < 4):
                    hh = h.subgraph(vertices=c)
                    if total_weight(hh) > best_w:
                        best = hh
                        best.name(name)
                        best_w = total_weight(best)
                if directed:
                    p.add_constraint(p.sum(edge[F(e)] for e in G.edge_boundary(c)), min=1)
                    c = set(c)
                    cbar = (v for v in G if v not in c)
                    p.add_constraint(p.sum(edge[F(e)] for e in G.edge_boundary(cbar, c)), min=1)
                else:
                    p.add_constraint(p.sum(edge[F(e)] for e in G.edge_boundary(c)), min=2)
                if induced:
                    p.add_constraint(p.sum(vertex[u] for u in c) <= len(c) - 1)
        if G.get_pos():
            best.set_pos({u: pp for u, pp in G.get_pos().items() if u in best})
        return (best_w, best) if use_edge_labels else best

It shows a longest cycle of Z is of length 5 and a longest induced cycle is of length 4. But we can find a cycle of Z of length 9 and a longest induced cycle of length 5. Did I miss something in my excerpt?

[dcoudert]

Right, the subtour elimination constraints are not correct for the longest cycle. I proposed another formulation in #37181.