#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 15 11:48:39 2025

@author: jf
"""
import numpy as np
from matplotlib import pyplot as plt
import scipy.sparse as sp
import networkx as nx
import gurobipy as gp
from edges import * 

#------------------------------------------------------------------------------
# TSP symétrique 
# Formulation DFJ symétrique avec élimination des sous-circuits
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
# Données distances sur les arêtes
#------------------------------------------------------------------------------
n=9
m=n*(n-1)//2  
distance=np.array([2, 5, 3, 3, 1, 2, 1, 5, 1, 4, 3, 1, 1, 2, 2, 1, 4, 1, 4, 1, 5, 5,
        1, 2, 4, 5, 4, 3, 2, 2, 3, 2, 5, 3, 5, 5])
#------------------------------------------------------------------------------
# Option d'affichage (graphe)
#------------------------------------------------------------------------------
disp_graph=True
#------------------------------------------------------------------------------
# Graphe complet, matrice d'incidence (sommet-aretes) complète
#------------------------------------------------------------------------------
G = nx.complete_graph(n)
id_edges = np.array(G.edges)

#------------------------------------------------------------------------------
# Construction de la matrice E (sommets/sommets/aretes) qui permet de retrouver 
# l'arete k associée aux sommets (i,j) : E[i,j]=k
#------------------------------------------------------------------------------
E=nodes_edges(n,id_edges)
#------------------------------------------------------------------------------
# Définition du modèle de base (avec sous-circuits) et resolution
#------------------------------------------------------------------------------
model = gp.Model("TSP_SYM")
x = model.addMVar(m, vtype=gp.GRB.BINARY)        # variables binaires
model.setObjective(distance @ x, gp.GRB.MINIMIZE)
model.addConstr( A@x == 2)

def subtourelim(model, where):

    if where == gp.GRB.Callback.MIPSOL:
        # Construction de la liste d'arêtes associée à la solution
        edgesol = model.cbGetSolution(x)
        
        # Construction du graphe associé à la solution
        esol=np.where(np.array(edgesol)==1)[0]
        G = nx.Graph()
        G.add_edges_from(id_edges[esol,:])
        
        # ncc : nombre de composantes connexes du graphe G
        ncc=nx.number_connected_components(G)
        if ncc >1:
            
            B = sp.lil_array((ncc, m), dtype=int)
            b = np.zeros(ncc, dtype=int)
    
            . . . A  COMPLETER
	    
            model.cbLazy( B@x <= b )    # ajout de la contrainte 

            
# Optimisation par Gurobi avec des lazy constraints
model.Params.lazyConstraints = 1
model.optimize(subtourelim)   # La fonction subtourelim est passée en argument

# Recuperation de la solution
sol_edges = model.X
dmin = model.ObjVal
print("    distance minimale =", dmin)

#------------------------------------------------------------------------------
# Construction du graphe associé à la solution
#------------------------------------------------------------------------------
esol=np.where(np.array(sol_edges)==1)[0]
G = nx.Graph()
G.add_edges_from(id_edges[esol,:])

#------------------------------------------------------------------------------
# Affichage du graphe
#------------------------------------------------------------------------------
if disp_graph:    
    plt.close('all')
    
    plt.figure()
    plt.title('Chemin/graphe optimal')
    pos = nx.shell_layout(G)
    nx.draw(G, pos=pos, with_labels=True, font_weight='bold')
    plt.axis('square')
    plt.pause(0.1)
    plt.show()

print('graphe connexe :',nx.is_connected(G))
ncc=nx.number_connected_components(G)
print('number of cc :', ncc)