from ortools.linear_solver import pywraplp as OR

def find_subtour(nodes, edges):
    """Identifies the first closed loop in the list of active edges."""
    unvisited = set(nodes)
    while unvisited:
        start_node = next(iter(unvisited))
        cycle = []
        current = start_node
        while current is not None:
            cycle.append(current)
            if current in unvisited:
                unvisited.remove(current)
            # Find the next node in the path (i -> j)
            next_node = next((j for i, j in edges if i == current), None)
            if next_node == start_node:
                return cycle
            if next_node not in unvisited:
                break
            current = next_node
    return None

def solve_tsp(dist_matrix):
    n = len(dist_matrix)
    nodes = range(n)
    
    # Initialize the SCIP solver
    m = OR.Solver('TSP_Class_Demo', OR.Solver.SCIP_MIXED_INTEGER_PROGRAMMING)
    
    # 1. Variables: x[i,j] is 1 if we travel from city i to city j
    x = {}
    for i in nodes:
        for j in nodes:
            if i != j:
                x[i, j] = m.BoolVar(f'x_{i}_{j}')

    # 2. Objective: Minimize total travel distance
    m.Minimize(sum(dist_matrix[i][j] * x[i, j] for (i, j) in x))

    # 3. Degree Constraints (The Assignment Problem)
    # Every city must have exactly one outgoing and one incoming edge
    for i in nodes:
        m.Add(sum(x[i, j] for j in nodes if i != j) == 1)
        m.Add(sum(x[j, i] for j in nodes if i != j) == 1)

    iteration = 1
    while True:
        status = m.Solve()
        
        if status != OR.Solver.OPTIMAL:
            print("No feasible solution found.")
            break

        # Convert solver output to a list of chosen edges
        active_edges = [k for k, var in x.items() if var.solution_value() > 0.5]
        print(f'{active_edges=}')
        subtour = find_subtour(nodes, active_edges)

        # If the first subtour found includes all nodes, we are done!
        if len(subtour) == n:
            print(f"Optimal TSP Tour Found! Distance: {m.Objective().Value()}")
            return active_edges

        # 4. Subtour Elimination: Add a constraint to break the discovered loop
        print(f"Iteration {iteration}: Subtour {subtour} detected. Breaking it...")
        
        # Logic: The number of edges connecting nodes within this subset 
        # must be less than the number of nodes in the subset.
        m.Add(sum(x[i, j] for i in subtour for j in subtour if i != j) <= len(subtour) - 1)
        
        iteration += 1

Input:

dist_matrix = [
    [0, 10, 8, 9, 7],
    [10, 0, 10, 5, 6],
    [8, 10, 0, 8, 9],
    [9, 5, 8, 0, 6],
    [7, 6, 9, 6, 0]
]

tour = solve_tsp(dist_matrix)
print("Final Tour:", tour)

active_edges=[(0, 2), (1, 4), (2, 0), (3, 1), (4, 3)]
Iteration 1: Subtour [0, 2] detected. Breaking it...
active_edges=[(0, 2), (1, 4), (2, 3), (3, 1), (4, 0)]
Optimal TSP Tour Found! Distance: 34.0
Final Tour: [(0, 2), (1, 4), (2, 3), (3, 1), (4, 0)]