Module tsp.tsp_cpx

"""
    TSP - Solving the symmetric TSP by means of CPLEX's MIP solver and docplex
"""
import networkx as nx
from tsp.tsp_plot import plotEdgeList

from docplex.mp.model import Model
from docplex.mp.callbacks.cb_mixin import ConstraintCallbackMixin
from docplex.mp.solution import SolveSolution
from cplex.callbacks import UserCutCallback
from cplex.callbacks import LazyConstraintCallback

from itertools import combinations as combinat

cpx_wait = None

#---------------------------------------------------------------------

class __FindSCEcuts( ConstraintCallbackMixin, UserCutCallback ):
    """
    Finds sub-cycle elimination (SCE) constraints violated by a 
    fractional solution to the current node LP. Limit maybe to first or 
    first three nodes in search tree!!!
    """
    
    def __init__(self, env):
        UserCutCallback.__init__(self, env)
        ConstraintCallbackMixin.__init__(self)
        self.eps = 1e-5
        self.nb_cuts = 0
    
    def __call__(self):
        
        x = self.x         # docplex edge decision variables 
        G = self.G         # graph obtained from problem.get_graph()
        n = len(G)         # number of nodes in the problem/graph
        model = self.model # handle to the model
        
        # Get list of edges with positive value in the current LP solution
        elist = list(e for e in x.keys() if self.get_values(x[e]._index) > self.eps)
        
        # If with animation, show the (fractional) solution
        if self.anim:
            objv = self.get_objective_value()
            efrac = [ e for e in elist if self.get_values(x[e]._index) < 0.99 ]
            plotEdgeList(self.problem,elist,specialEdges=efrac,wait_func=cpx_wait,\
                         title="Fractional LP solution of objval="+str(objv) )
                       
        # Build the solution's support graph, i.e. the sub-graph induced by
        # the edges in the list "elist" found above.
        supG = nx.edge_subgraph(G, elist)
        
        # If the support graph G is not connected, build the SCE constraint
        # based on the graph's smallest component. Note that this cannot
        # happen if the 2-commodity flow formulation is used
        S = None
        if self.modtype != '2cf':
            S = min( nx.connected_components(supG), key=len )
            weight = 1.0
        
        # If the graph is connected, obtain the minimum-weight cut set
        # using the algorithm of Stoer and Wagner
        if (S is None) or (len(S)==n):
            # Set the edge weights equal to the variable's solution values
            for e in elist:
                supG[e[0]][e[1]]['weight'] = self.get_values( x[e]._index )
            weight, part = nx.stoer_wagner(supG)    
            S = part[0] if len(part[0]) < len(part[1]) else part[1]            
        
        # If the cutset constraint is violated, include the cut in 
        # form of the equivalent subcycle elimination constraint
        if ( len(S) > 2 ) and ( weight < 1.98):
            E = list(combinat(sorted(S),2))
            sce = self.linear_ct_to_cplex( model.sum(x[e] for e in E) <= len(S)-1 )
            self.add( *sce )      
            self.nb_cuts += 1
            
#---------------------------------------------------------------------

class __FindLazySCEs( ConstraintCallbackMixin, LazyConstraintCallback ):
    """
    Finds sub-cycle elimination (SCE) constraints violated by an
    integer solution to the current node LP. Note that this can only
    be the case if the 2-commodity flow formulation is not used
    and instead the fractional 2-matching problem is the initial
    relaxation.
    """
    
    def __init__(self, env):
        LazyConstraintCallback.__init__(self, env)
        ConstraintCallbackMixin.__init__(self)
        self.one = 0.99999
        self.nb_lazy = 0
        
    def __call__(self):     
        """ 
        Solution is integer. Check for sub-cycles and if there is one
        store the smallest and include the corresponding subcycle
        elimination constraint.
        """
        x = self.x         # docplex edge decision variables 
        G = self.G         # graph obtained from problem.get_graph()
        n = len(G)         # number of nodes in the problem/graph
        model = self.model # handle to the docplex model
        
        # Get list of edges selected in the integer solution
        elist = list(e for e in x.keys() if self.get_values(x[e]._index) > self.one)
        
        # If with animation, show the (fractional) solution
        if self.anim:
            objv = self.get_objective_value()
            __title = "Checking lazy constraints. Integer solution objval="+str(objv)
            plotEdgeList(self.problem, elist, title=__title, wait_func=cpx_wait )
        
        # Find smallest component of the support graph
        S = min( nx.connected_components(nx.edge_subgraph(G,elist)), key=len )
        
        # Add subcycle elimination constraint if not connected
        if len(S) < n:
            E = list(combinat(sorted(S),2))
            sce = self.linear_ct_to_cplex( model.sum(x[e] for e in E) <= len(S)-1 )
            self.add( *sce )      
            self.nb_lazy += 1
        
#---------------------------------------------------------------------

def __makeRoute( edges ):
    """
    Transforms solution given as edge list to a solution in permutation 
    form.
    
    Parameters:
        edges : list of int
            The list of edges selected in the solution.
            
    Returns:
        route : list of int
            The route as a list of nodes.             
    """
    e = edges.pop(0)
    route = [ e[0], e[1] ]
    while len(edges) > 0:
        last = route[-1]
        edge = next( e for e in edges if last in e )
        succ = edge[1] if last==edge[0] else edge[0]
        route.append(succ)
        edges.remove(edge)        
    return( route )

#---------------------------------------------------------------------

def __build_model( problem, use2CNF ):
    """
    Set up the mathematical program of the TSP using either just the
    2-matching problem or the 2-commodity network flow formulation.
    
    Parameters :
    
        problem : class
            TSP problem object has returned by function load_problem
            of package tsplib95
        use2CNF : boolean
            If true, the 2-commodity flow formulations is used.
            Otherwise 2-matching+subcycle elimination
    
    Returns :
        model : class
            handle to the model
        x     : list of binary variables
    """
    
    # number of nodes 
    n = problem.dimension 
    
    # Create Cplex mathematical programming model
    if use2CNF:
        model = Model("2-commodity flow TSP formulation")
    else:
        model = Model("2-matching + subcycle elimination")    
    
    # List of nodes and edges
    nodes = list(problem.get_nodes())
    edges = list( combinat( nodes, 2 ) )
    
    # Function to access edges incident to a node
    get_edges = lambda node : list( (i,node) for i in nodes if i < node ) \
                             +list( (node,i) for i in nodes if i > node )

    # list of binary edge variables
    x = model.binary_var_dict( edges, name = 'x', key_format = '_%s') 
      
    # Add the objective function to the model
    model.minimize( model.dotf( x, lambda e: problem.wfunc(e[0],e[1]) ) )
    
    # Add degree constraint for each node
    model.add_constraints_( (model.sum(x[e] for e in get_edges(i)) == 2) for i in nodes )
    
    if use2CNF:
        # Matrix of flow variables f(i,j)
        f = model.continuous_var_matrix(nodes, nodes, name = 'f', key_format='_%s')
  
        # Add flow constraints for each node except the first
        model.add_constraints_( (model.sum((f[i,j]-f[j,i]) for i in nodes if i!=j) == 2)\
                                 for j in nodes[1:] )
        
        # Add constraints fij + fji = n*x(ij) for each edge (ij)
        model.add_constraints_( (f[e[0],e[1]]+f[e[1],e[0]] - n*x[e] == 0) for e in edges )
         
    return model, x

#---------------------------------------------------------------------

def __setStartSol ( G, model, x, route ):
    """
    Pass an initial solution given by "route" to the cplex solver.
    
    Parameters:
        G     : networkx graph
            The graph representing the problem instance
        model : class
            The model as returned by docplex.mp.model.Model()
        x     : dict of variables
            The edge variables
        route : list of int
            The initial solution give as sequence of nodes
    """
    routeLen = sum([G[route[i-1]][route[i]]['weight'] for i in range(1,len(route))]) 
    routeLen += G[route[0]][route[-1]]['weight']
    
    var_value = dict( zip(x.values(),[0]*G.number_of_edges()) )
    
    edge = lambda u,v: (u,v) if u < v else (v,u)
    for u,v in zip( route[:-1],route[1:] ):
        var_value[x[edge(u,v)]] = 1
    var_value[x[edge(route[-1],route[0])]] = 1
                      
    sol = SolveSolution(model, var_value_map=var_value )
    model.add_mip_start(sol)
    
    return routeLen
    
#---------------------------------------------------------------------

def CPXtspSolve( problem, G = None, startRoute=None, formulation='SCE',\
                 nodeLim=None, anim=False, wait_func=None, log_output=False ):
    """
    Uses CPLEX's MIP solver via docplex for solving the symmetric TSP. 
    Subcycle elimination constraints are included via cut-callback functions.
      
    Parameters:
        **problem** : class  
            TSP problem object has returned by function load_problem
            of package tsplib95.  
        **G** : object  
            If not None, the networkx graph of the problem.
            Default = None.  
        **startRoute** : list of int, optional  
            If not None it specifies a feasible route given as the 
            sequence in which the nodes are visited. Default=None  
        **formulation** : string, optional  
            If equal to 'SCE', the initial problem formulation
            is the 2-matching problem. Subcycle elimination constraints
            are then generated on the fly as user cuts and lazy constraints.
            If equal to '2CF', the 2-commodity flow formulation is applied
            and no subcycle elimination constraints are generated.
            If equal to '2CF-SCE', again the 2-commodity flow formulation
            is used but subcycle elimination constraints are generated as cuts.  
        **nodeLim** : int, optional  
            Limit on the nodes to be enumerated. Default=None                 
        **anim** : bool, optional  
            If true, the (fractional) LP solution obtained before 
            calling a cut callback is displayed (attention: can be 
            time consuming). Default=False.        
        **wait_func** : function, optional  
            Function that is called to wait for any user action if run
            in animation mode (anim=True), so that execution first continues
            on return from the wait_func. Default=None      
        **log_output** : bool, optional  
            If true, the Cplex solver sends information about the solution progress
            to the screen. Otherwise, it is just silent.                  

    Returns:
        **lowbnd** : float  
            A lower bound or the optimal solution value.  
        **routeLen** : int  
            Length of the computed route.  
        **route** : list of int  
            The route computed.
    """
    
    # Waiting function if used in animated way
    global cpx_wait
    cpx_wait = wait_func
    
    # Get ther right edge weight function from tsplib95
    try:
        _ = problem.wfunc # tsplib95 version before 0.7.1
    except:
        problem.wfunc = problem._wfunc # tsplib95, version 0.7.1

    # Creat the problem graph if not passed as argument
    G_local = ( G is None )
    if G_local:
        G = problem.get_graph()
        G.remove_edges_from(nx.selfloop_edges(G))
                        
    # Build the model
    use2CF = '2CF' in formulation.upper()
    addSCE = (not use2CF) or ('SCE' in formulation.upper())
    model, x  = __build_model( problem, use2CF )
    
    #Initial route provided?
    routeLen = None
    route = startRoute
    if (not route is None): routeLen = __setStartSol(G, model, x, route)
   
    # Set the node limit if given
    if not nodeLim is None: model.parameters.mip.limits.nodes.set(nodeLim)

    # Register the lazy constraint callback
    if not use2CF:
        model.register_callback( __FindLazySCEs )
        __FindLazySCEs.G = G
        __FindLazySCEs.x = x
        __FindLazySCEs.anim = anim
        if anim: __FindLazySCEs.problem = problem
    
    # Register the user cut callback
    if addSCE:
        model.register_callback( __FindSCEcuts )
        __FindSCEcuts.G = G
        __FindSCEcuts.x = x 
        __FindSCEcuts.anim = anim
        __FindSCEcuts.modtype = '2cf' if use2CF else 'SCE'
        if anim: __FindSCEcuts.problem = problem
    # Solve problem with lazy constraints
    sol = model.solve(log_output=log_output)
    
    # Obtain status of the solution
    best_bnd = model.solve_details.best_bound
    if not sol is None:
        objval = model.objective_value
        if (route is None) or (objval < routeLen ):
            edges = list ( (sol.get_value_dict(x,keep_zeros=False)).keys() )
            route = __makeRoute( edges )
            routeLen = objval

    model.end()
    cpx_wait = None

    # If graph is not a local object, remove the edge attribute "var" again. 
    # Restore also the original weights
    if not G_local:
        for e in G.edges: G[e[0]][e[1]]['weight'] = problem.wfunc(e[0],e[1])

    return best_bnd, routeLen, route