"""
tsp.py:  solve the traveling salesman problem 

minimize the travel cost for visiting n customers exactly once
approach:
    - start with assignment model
    - add cuts until there are no sub-cycles

Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
"""
import math
import time
import random
import networkx
from gurobipy import *

def solve_tsp(V,c):
    """solve_tsp -- solve the traveling salesman problem 
       - start with assignment model
       - add cuts until there are no sub-cycles
    Parameters:
        - V: set/list of nodes in the graph
        - c[i,j]: cost for traversing edge (i,j)
    Returns the optimum objective value and the list of edges used.
    """
    
    def addcut(cut_edges):
        G = networkx.Graph()
        G.add_edges_from(cut_edges)
        Components = networkx.connected_components(G)

        if len(Components) == 1:
            return False
        for S in Components:
            model.addConstr(quicksum(x[i,j] for i in S for j in S if j>i) <= len(S)-1)
            print "cut: len(%s) <= %s" % (S,len(S)-1)
        return True


    def addcut2(cut_edges):
        G = networkx.Graph()
        G.add_edges_from(cut_edges)
        Components = networkx.connected_components(G)

        if len(Components) == 1:
            return False
        for S in Components:
            T = set(V) - set(S)
            print "S:",S
            print "T:",T
            model.addConstr(quicksum(x[i,j] for i in S for j in T if j>i) >= 2)
            print "cut: %s <--> %s >= 2" % (S,T), [(i,j) for i in S for j in T if j>i]
        return True


    # main part of the solution process:
    model = Model("tsp")

    # model.Params.OutputFlag = 0 # silent/verbose mode
    x = {}
    for i in V:
        for j in V:
            if j > i:
                x[i,j] = model.addVar(ub=1, name="x(%s,%s)"%(i,j))
    model.update()
    
    for i in V:
        model.addConstr(quicksum(x[j,i] for j in V if j < i) + \
                        quicksum(x[i,j] for j in V if j > i) == 2, "Degree(%s)"%i)

    model.setObjective(quicksum(c[i,j]*x[i,j] for i in V for j in V if j > i), GRB.MINIMIZE)

    EPS = 1.e-6
    while True:
        model.optimize()
        edges = []
        for (i,j) in x:
            if x[i,j].X > EPS:
                edges.append( (i,j) )

        if addcut(edges) == False:
            if model.IsMIP:     # integer variables, components connected: solution found
                break
            for (i,j) in x:     # all components connected, switch to integer model
                x[i,j].VType = "B"
            model.update()

    return model.ObjVal,edges


def distance(x1,y1,x2,y2):
    """distance: euclidean distance between (x1,y1) and (x2,y2)"""
    return math.sqrt((x2-x1)**2 + (y2-y1)**2)

def make_data(n):
    """make_data: compute matrix distance based on euclidean distance"""
    V = range(1,n+1)
    x = dict([(i,random.random()) for i in V])
    y = dict([(i,random.random()) for i in V])
    c = {}
    for i in V:
        for j in V:
            if j > i:
                c[i,j] = distance(x[i],y[i],x[j],y[j])
    return V,c
    
            
if __name__ == "__main__":
    import sys

    # Parse argument
    if len(sys.argv) < 2:
        print "Usage: %s instance" % sys.argv[0]
        exit(1)
        # n = 200
        # seed = 1
        # random.seed(seed)
        # V,c = make_data(n)

    from read_tsplib import read_tsplib
    try:
        V,c,x,y = read_tsplib(sys.argv[1])
    except:
        print "Cannot read TSPLIB file",sys.argv[1]
        exit(1)

    obj,edges = solve_tsp(V,c)

    print
    print "Optimal tour:",edges
    print "Optimal cost:",obj
    print