text: # -*- coding: utf-8 -*-
  ###############################################################################
  """ This file is a part of the VeRyPy classical vehicle routing problem
  heuristic library and implements inter route (i.e. between two or more routes)
  local search improvement heuristics such as one and two point move etc. All 
  operators assume we start from a feasible solution. Also, all functions 
  implementing the operations have the following signature:
  
  do_Y_move(route1_data, route2_data, ... , D,C,d,L, strategy, best_delta), where
  
  * routeZ_data is a RouteData object, the data members are
          1. route as a list of node indices. Must begin and end to depot (0).
          2. current route cost
          3. current route demand, can be None is C = None
  * D is the numpy-compatibe distance matrix as with intra route heuristics.
  * C is the optional capacity constraint (can be None). If set, also give d.
  * d is a list of customer demands, and of the depot it is d[0]=0,
      The parameter must be set if  C is given.
  * L is the maximun route length/duration/cost constraint
  * strategy is either FIRST_ACCEPT (default) or BEST_ACCEPT. 
      First accept returns a modified route as soon as the first improvement is  
      encoutered. Best accept tries all possible combinations and returns the
      best one.
  * best_delta is the required level of change in the in route cost. It can be
     used to  set the upper bound (requirement) for the improvement. Usually it 
     is None, which sets the level to 0.0. In that case only improving deltas are
     accepted. If set to (large) positive value, best worsening move can also be
     returned.
     
  All inter route improvement operators return the new improved routes and the
  improvement (delta) as (n_r+1)-tuple, where n_r is the number of routes given 
  to the operator or (None,...,None) if no improvements were found."""
  ###############################################################################
  
  #TODO the efficiency of these implementations could be greatly improved using
  # numba, doubly linked lists for routes, and in place editing of routes.
  #  - use doubly linked list (dllist) as the data structure for the routes
  #  - do not copy routes, make changes in place
  #  - use numba (and numba compatible dllist) or even cython 
  
  # Written in Python 2.7, but try to maintain Python 3+ compatibility
  from __future__ import print_function
  from __future__ import division
  
  from itertools import product, permutations
  
  from routedata import RouteData
  from local_search import LSOPT, ROUTE_ORDER_SENSITIVE_OPERATORS
  from config import COST_EPSILON as S_EPS
  from config import CAPACITY_EPSILON as C_EPS
  
  __author__ = "Jussi Rasku"
  __copyright__ = "Copyright 2018, Jussi Rasku"
  __credits__ = ["Jussi Rasku"]
  __license__ = "MIT"
  __maintainer__ = "Jussi Rasku"
  __email__ = "jussi.rasku@jyu.fi"
  __status__ = "Development"
  
  # Create a decorator to mark routes sensitive to the route order
  def routeordersensitive(f):
      ROUTE_ORDER_SENSITIVE_OPERATORS.add(f)
      return f
      
  def do_2optstar_move(route1_data, route2_data, D, d=None, 
                        C=None, L=None, # constraints
                        strategy=LSOPT.FIRST_ACCEPT,
                        best_delta = None):
      
      """ 2-opt* inter-route local search operation for the symmetric distances D
      Remove 2 edges from different routes and check if swapping the edge halves
      (in two different ways) would yield an improvement, while making sure the
      move does not violate constraints.
      """
  
      # use 2-opt
      if route1_data==route2_data:
          raise ValueError("Use do_2opt_move to find intra route moves")
  
      # make sure we have the aux data for constant time feasibility checks   
      if not isinstance(route1_data, RouteData):
          route1_data = RouteData(*route1_data)
      if not isinstance(route2_data, RouteData):
          route2_data = RouteData(*route2_data)
      if not route1_data.aux_data_updated:
          route1_data.update_auxiliary_data(D,d)
      if not route2_data.aux_data_updated:
          route2_data.update_auxiliary_data(D,d)
  
      if not best_delta:
          best_delta = 0
      best_move = None
      accept_move = False
      
      #print("REMOVEME: 2opt* on %s %s"%(list(route1_data.route), list(route2_data.route)) )
      
      
      for i in xrange(0,len(route1_data.route)-1):
          for j in xrange(0,len(route2_data.route)-1):
              a = route1_data.route[i]
              b = route1_data.route[i+1]
              c = route2_data.route[j]
              d = route2_data.route[j+1]
              
              #print("REMOVEME: attempt remove %d-%d and %d-%d"%(a,b,c,d) )
              
              # a->c b->d
              #       __________
              #      /          \
              # 0->-a   b-<-0-<-c   d->-0
              #         \___________/  
              #
              delta = D[a,c] + D[b,d] \
                       -D[a,b]-D[c,d]
              if delta+S_EPS<best_delta:
                  # is an improving move, check feasibility 
                  constraint_violated = False
                  r1_new_demand = None
                  r2_new_demand = None
                  if C:
                      r1_new_demand = route1_data.fwd_d[i]+route2_data.fwd_d[j]
                      r2_new_demand = route1_data.rwd_d[i+1]+route2_data.rwd_d[j+1]
                      if r1_new_demand-C_EPS>C or r2_new_demand-C_EPS>C:
                          constraint_violated = True
                  if not constraint_violated and L and (    
                    route1_data.fwd_l[i]+D[a,c]+route2_data.fwd_l[j]-S_EPS>L or
                    route1_data.rwd_l[i+1]+D[b,d]+route2_data.rwd_l[j+1]-S_EPS>L):
                      constraint_violated = True
                  
                  if not constraint_violated:
                      # store segments, deltas, and demands
                      best_move = (((None,i+1,1), (j,None,-1),
                                    D[a,c]-D[a,b], r1_new_demand),
                                   ((None,i,-1), (j+1,None,1),
                                    D[b,d]-D[c,d], r2_new_demand))
                      best_delta = delta
                      
                      if strategy==LSOPT.FIRST_ACCEPT:
                          accept_move = True
                          break # j loop
              
              # a->d c->b
              #       ______________
              #      /              \
              # 0->-a   b-<-0-<-c   d->-0
              #         \______/  
              #
              delta = D[a,d] + D[b,c] \
                       -D[a,b]-D[c,d]
              if delta+S_EPS<best_delta:
                  # is an improving move, check feasibility 
                  constraint_violated = False
                  r1_new_demand = None
                  r2_new_demand = None
                  if C:
                      r1_new_demand = route1_data.fwd_d[i]+route2_data.rwd_d[j+1]
                      r2_new_demand = route1_data.rwd_d[i+1]+route2_data.fwd_d[j]
                      if r1_new_demand-C_EPS>C or r2_new_demand-C_EPS>C:
                          constraint_violated = True
                  if not constraint_violated and L and (
                    route1_data.fwd_l[i]+D[a,d]+route2_data.rwd_l[j+1]-S_EPS>L or
                    route1_data.rwd_l[i+1]+D[b,c]+route2_data.fwd_l[j]-S_EPS>L):
                      constraint_violated = True
                      
                  if not constraint_violated:
                      # store segments, deltas, and demands
                      best_move = (((None,i+1,1), (j+1,None,1),
                                    D[a,d]-D[a,b], r1_new_demand),
                                   ((None,i,-1), (j,None,-1),
                                    D[b,c]-D[c,d], r2_new_demand))
                      best_delta = delta
                      
                      if strategy==LSOPT.FIRST_ACCEPT:
                          accept_move = True
                          break # j loop
                       
          if accept_move:
              break # i loop
                  
      if best_move:
          # unpack the move
          ((r1_sgm1, r2_sgm1, r1_delta, r1_new_demand),
           (r1_sgm2, r2_sgm2, r2_delta, r2_new_demand)) = best_move
           
          return (
              # route 1
              RouteData(
                  route1_data.route[r1_sgm1[0]:r1_sgm1[1]:r1_sgm1[2]]+\
                  route2_data.route[r2_sgm1[0]:r2_sgm1[1]:r2_sgm1[2]],
                  route1_data.cost+r1_delta,
                  r1_new_demand),
              # route 2
              RouteData(
                  route1_data.route[r1_sgm2[0]:r1_sgm2[1]:r1_sgm2[2]]+\
                  route2_data.route[r2_sgm2[0]:r2_sgm2[1]:r2_sgm2[2]],
                  route1_data.cost+r1_delta,
                  r1_new_demand),
              # delta
              best_delta)
   
      return None, None, None
  
  # These define all possible 3opt moves. Describes which segment ends to connect
  #  This makes implementing the do_3optstar_move (below) operator simpler.
  
  #     
  #    __0  1__      route 1
  #   /        \        
  # n0---2  3--n0    route 2
  #   \__    __/      
  #      4  5        route 3  
  #
  MOVES_3OPTSTAR_3ROUTES = (
      #((0, 1), (2, 3), (4, 5)), 
      #2-opt
      ((0, 1), (2, 4), (3, 5)), 
      ((0, 1), (2, 5), (4, 3)), 
      ((0, 4), (2, 3), (1, 5)),
      ((0, 5), (2, 3), (4, 1)),
      ((0, 2), (1, 3), (4, 5)), 
      ((0, 3), (2, 1), (4, 5)), 
      #3-opt
      ((0, 2), (5, 3), (4, 1)),
      ((0, 2), (4, 3), (1, 5)),
      ((0, 3), (2, 5), (4, 1)),
      ((0, 3), (2, 4), (1, 5)),
      ((0, 4), (2, 1), (3, 5)),
      ((0, 4), (1, 3), (2, 5)),
      ((0, 5), (2, 1), (4, 3)),
      ((0, 5), (1, 3), (4, 2)),
  )
  
  #
  #    _0  1_2  3_     route 1
  #   /           \  
  # n0            n0  
  #   \____    ___/  
  #        4  5        route 2
  #
  # All moves have 3 segments on the second route
  #  (connected by the two first edges)
  MOVES_3OPTSTAR_2ROUTES = (
      #((0, 1), (2, 3), (4, 5)), 
      
      #2-opt
      ((3, 5), (0, 1), (2, 4)), # note, route 2 first
      ((4, 3), (0, 1), (2, 5)), # note, route 2 first
      ((0, 4), (3, 2), (1, 5)), 
      ((0, 5), (4, 1), (2, 3)), 
      ((4, 5), (0, 2), (1, 3)), # note, route 2 first
      #((0, 3), (2, 1), (4, 5)), would form a loop
      #3-opt
      ((3, 5), (0, 2), (1, 4)), # note, route 2 first
      ((4, 3), (0, 2), (1, 5)), # note, route 2 first
      ((0, 3), (4, 1), (2, 5)),
      ((0, 3), (4, 2), (1, 5)),
      #((0, 4), (2, 1), (3, 5)), would form a loop
      ((0, 4), (3, 1), (2, 5)),
      #((0, 5), (2, 1), (4, 3)), would form a loop
      ((0, 5), (4, 2), (1, 3)),
  )
  
  @routeordersensitive
  def do_3optstar_2route_move(route1_data, route2_data, D, demands=None, 
                        C=None, L=None, # constraints
                        strategy=LSOPT.FIRST_ACCEPT,
                        best_delta = None):
  
      if route1_data==route2_data:
          raise ValueError("Use do_3opt_move to find intra route moves")
          
      # make sure we have the aux data for constant time feasibility checks   
      if not isinstance(route1_data, RouteData):
          route1_data = RouteData(*route1_data)
      if not isinstance(route2_data, RouteData):
          route2_data = RouteData(*route2_data)
      if not route1_data.aux_data_updated:
          route1_data.update_auxiliary_data(D,demands)
      if not route2_data.aux_data_updated:
          route2_data.update_auxiliary_data(D,demands)
          
      if not best_delta:
          best_delta = 0
      best_move = None
      accept_move = False
      
      # segment end nodes and cumulative d and l at those nodes is specified 
      # by the indices i,j,k. In the two route case this means:
      #
      #    _0/i  1_2/j  3_     route 1
      #   /               \  
      # n0                n0  
      #   \____      _____/  
      #        4/k  5        route 2
      #
      end_n = [0]*6
      cum_d = [0]*6
      cum_l = [0]*6
      
      for i in xrange(0,len(route1_data.route)-1):
          end_n[0] = route1_data.route[i]
          end_n[1] = route1_data.route[i+1]
  
          # make it so that i<j
          for j in xrange(i+1 ,len(route1_data.route)-1):
              
              # the edge endpoints
              end_n[2] = route2_data.route[j]
              end_n[3] = route2_data.route[j+1]
              if C:
                  cum_d[0] = route1_data.fwd_d[i]
                  cum_d[1] = route1_data.rwd_d[i+1]\
                            -route1_data.rwd_d[j+1]
                  cum_d[2] = route1_data.fwd_d[j]\
                             -route1_data.fwd_d[i]
                  cum_d[3] = route1_data.rwd_d[j+1]
              if L:
                  cum_l[0] = route1_data.fwd_l[i]
                  cum_l[1] = route1_data.rwd_l[i+1]\
                             -route1_data.rwd_l[j]
                  cum_l[2] = route1_data.fwd_l[j]\
                             -route1_data.fwd_l[i+1]
                  cum_l[3] = route1_data.rwd_l[j+1]
                  
              
              for k in xrange(0,len(route2_data.route)-1):
                  
                  # the edge endpoints
                  end_n[4] = route2_data.route[k]
                  end_n[5] = route2_data.route[k+1]
                  if C:
                      cum_d[2] = route2_data.fwd_d[k]
                      cum_d[3] = route2_data.rwd_d[k+1]
                  if L:
                      cum_l[2] = route2_data.fwd_l[k]
                      cum_l[3] = route2_data.rwd_l[k+1]
                  
                  removed_weights = D[end_n[0],end_n[1]]+\
                                    D[end_n[2],end_n[3]]+\
                                    D[end_n[4],end_n[5]]
        
                  for e1,e2,e3 in MOVES_3OPTSTAR_2ROUTES:
                      e1_wt = D[end_n[e1[0]],end_n[e1[1]]]
                      e2_wt = D[end_n[e2[0]],end_n[e2[1]]]
                      e3_wt = D[end_n[e3[0]],end_n[e3[1]]]
                  
                      delta = e1_wt+e2_wt+e3_wt-removed_weights
                      
                      if ((delta+S_EPS<best_delta) and
                          (not C or (cum_d[e1[0]]+cum_d[e1[1]]-C_EPS<C and
                           cum_d[e2[0]]+cum_d[e2[1]]+cum_d[e3[1]]-C_EPS<C)) and
                          (not L or ( cum_l[e1[0]]+cum_l[e1[1]]+e1_wt-S_EPS<L and
                            cum_l[e2[0]]+cum_l[e2[1]]+cum_l[e3[1]]+e2_wt+e3_wt-S_EPS<L))):
                          best_move = ((i,j,k),(e1,e2,e3),delta)
                          best_delta = delta  
                          
                          if strategy==LSOPT.FIRST_ACCEPT:
                              accept_move = True
                              break # move loop
                  if accept_move:
                      break # k loop
              if accept_move:
                  break # j loop
          if accept_move:
              break # i loop
      
      
      if best_move:
          raise NotImplementedError("Not yet completely implemented")
          
          # TODO: 
          # 1. unpack the move
          # 2. construct the modified routes
          # 3. return the modified routes and delta
            
          # unpack the move
          #(best_ijk,best_edges,best_delta) = best_move
          #best_e1, best_e2, best_e3 = best_edges
      
  #         new_route1 = 
  #        
  #            if edge[0]%2 and edge[1]%2:
  #                new_route = (
  #                        routes[r1_idx].route[None:r1_i :-1]+
  #                        routes[r2_idx].route[r2_i+1:None:1])
  #                new_cost = ( D[b, d]+
  #                        routes[r1_idx].rwd_l[r1_i+1]+
  #                        routes[r2_idx].rwd_l[r2_i+1])
  #                new_demand = (
  #                        routes[r1_idx].rwd_d[r1_i+1]+
  #                        routes[r2_idx].rwd_d[r2_i+1])
  #            # combine the head of route1 with the tail of route2
  #            elif not edge[0]%2 and edge[1]%2:
  #                new_route = (
  #                        routes[r1_idx].route[None:r1_i+1:1]+
  #                        routes[r2_idx].route[r2_i+1:None:1])
  #                new_cost = (D[a, d],
  #                        routes[r1_idx].fwd_l[r1_i]+
  #                        routes[r2_idx].rwd_l[r2_i+1])
  #                new_demand = (
  #                        routes[r1_idx].fwd_d[r1_i]+
  #                        routes[r2_idx].rwd_d[r2_i+1])
          
  def do_3optstar_3route_move(route1_data, route2_data, route3_data, D, demands=None, 
                        C=None, L=None, # constraints
                        strategy=LSOPT.FIRST_ACCEPT,
                        best_delta = None):
      
      """ 3-opt* inter-route local search operation for the symmetric distances D
      Remove 3 edges from different routes and check if reconnecting them in any
      configuration would improve the solution, while also making sure the
      move does not violate constraints.
      
      For two route checks, the first and second route can be same i.e. 
      route1_data==route2_data. However, if route2_data==route3_data or 
      route1_data==route3_data a ValueError is raised.
      """
  
      if route1_data==route2_data or\
         route2_data==route3_data or\
         route1_data==route3_data:
          raise ValueError("Use do_3opt_move to find intra route moves")
     
     # make sure we have the aux data for constant time feasibility checks   
      if not isinstance(route1_data, RouteData):
          route1_data = RouteData(*route1_data)
      if not isinstance(route2_data, RouteData):
          route2_data = RouteData(*route2_data)
      if not isinstance(route3_data, RouteData):
          route3_data = RouteData(*route3_data)
      if not route1_data.aux_data_updated:
          route1_data.update_auxiliary_data(D,demands)
      if not route2_data.aux_data_updated:
          route2_data.update_auxiliary_data(D,demands)
      if not route3_data.aux_data_updated:
          route3_data.update_auxiliary_data(D,demands)
  
      if not best_delta:
          best_delta = 0
      best_move = None
      accept_move = False
      
      # segment end nodes and cumulative d and l at those nodes
      # Note the indexing of segment end nodes here:
      #     
      #    __0/i  1__      route 1
      #   /          \        
      # n0---2/j  3--n0    route 2
      #   \__      __/      
      #      4/k  5        route 3  
      #
      end_n = [0]*6
      cum_d = [0]*6
      cum_l = [0]*6
      
      for i in xrange(0,len(route1_data.route)-1):
          end_n[0] = route1_data.route[i]
          end_n[1] = route1_data.route[i+1]
  
          # make it so that i<j if route1==route2
          for j in xrange(0 ,len(route2_data.route)-1):
              
              # the edge endpoints
              end_n[2] = route2_data.route[j]
              end_n[3] = route2_data.route[j+1]
              if C:
                  cum_d[0] = route1_data.fwd_d[i]
                  cum_d[1] = route1_data.rwd_d[i+1]
                  cum_d[2] = route2_data.fwd_d[j]
                  cum_d[3] = route2_data.rwd_d[j+1]
              if L:
                  cum_l[0] = route1_data.fwd_l[i]
                  cum_l[1] = route1_data.rwd_l[i+1]
                  cum_l[2] = route2_data.fwd_l[j]
                  cum_l[3] = route2_data.rwd_l[j+1]
              
              for k in xrange(0,len(route3_data.route)-1):
                  
                  # the edge endpoints
                  end_n[4] = route3_data.route[k]
                  end_n[5] = route3_data.route[k+1]
                  if C:
                      cum_d[4] = route3_data.fwd_d[k]
                      cum_d[5] = route3_data.rwd_d[k+1]
                  if L:
                      cum_l[4] = route3_data.fwd_l[k]
                      cum_l[5] = route3_data.rwd_l[k+1]
                  
                  removed_weights = D[end_n[0],end_n[1]]+\
                                    D[end_n[2],end_n[3]]+\
                                    D[end_n[4],end_n[5]]
        
                  for e1,e2,e3 in MOVES_3OPTSTAR_3ROUTES:
                      e1_wt = D[end_n[e1[0]],end_n[e1[1]]]
                      e2_wt = D[end_n[e2[0]],end_n[e2[1]]]
                      e3_wt = D[end_n[e3[0]],end_n[e3[1]]]
                      
                      delta = e1_wt+e2_wt+e3_wt-removed_weights
  
                      if ((delta+S_EPS<best_delta) and
                         (not C or ( cum_d[e1[0]]+cum_d[e1[1]]-C_EPS<C and
                                     cum_d[e2[0]]+cum_d[e2[1]]-C_EPS<C and
                                     cum_d[e3[0]]+cum_d[e3[1]]-C_EPS<C ))
                         and
                         (not L or ( cum_l[e1[0]]+cum_l[e1[1]]+e1_wt-S_EPS<L and
                                     cum_l[e2[0]]+cum_l[e2[1]]+e2_wt-S_EPS<L and
                                     cum_l[e3[0]]+cum_l[e3[1]]+e3_wt-S_EPS<L ))):
                          best_move = ((i,j,k),(e1,e2,e3),delta)
                          best_delta = delta  
                          
                          
                          if strategy==LSOPT.FIRST_ACCEPT:
                              accept_move = True
                              break # move loop
                  if accept_move:
                      break # k loop
              if accept_move:
                  break # j loop
          if accept_move:
              break # i loop
                  
      if best_move:
          # unpack the move
          (best_ijk,best_edges,best_delta) = best_move
      
          routes = [route1_data,route2_data,route3_data]
          ret = []
          for edge in best_edges:
              
              # Unfortunately, the splicing is a bit complex, but basic idea 
              #  is to use the move and its edges in best_edges to get the
              #  route to splice from and if the splice should be inverted.
              #
              # edge[i]//2 gets the route index of the move endpoint
              # edge[i]%2 tells if it is the head segment or tail segment
              #
              #    ___a c___    r1
              #   /         \
              # n0           n1
              #   \___   ___/
              #       b d       r2
              #                 
              
              r1_idx = edge[0]//2
              r2_idx = edge[1]//2
              r1_i = best_ijk[r1_idx]
              r2_i = best_ijk[r1_idx]
              a = routes[r1_idx].route[r1_i]
              b = routes[r1_idx].route[r1_i+1]
              c = routes[r2_idx].route[r2_i]
              d = routes[r2_idx].route[r2_i+1]
              
              # combine the tails of route1 and route2
              if edge[0]%2 and edge[1]%2:
                  new_route = (
                          routes[r1_idx].route[None:r1_i :-1]+
                          routes[r2_idx].route[r2_i+1:None:1])
                  new_cost = ( D[b, d]+
                          routes[r1_idx].rwd_l[r1_i+1]+
                          routes[r2_idx].rwd_l[r2_i+1])
                  new_demand = (
                          routes[r1_idx].rwd_d[r1_i+1]+
                          routes[r2_idx].rwd_d[r2_i+1])
              # combine the head of route1 with the tail of route2
              elif not edge[0]%2 and edge[1]%2:
                  new_route = (
                          routes[r1_idx].route[None:r1_i+1:1]+
                          routes[r2_idx].route[r2_i+1:None:1])
                  new_cost = (D[a, d],
                          routes[r1_idx].fwd_l[r1_i]+
                          routes[r2_idx].rwd_l[r2_i+1])
                  new_demand = (
                          routes[r1_idx].fwd_d[r1_i]+
                          routes[r2_idx].rwd_d[r2_i+1])
              # combine the heads of route1 and route2
              elif not edge[0]%2 and not edge[1]%2:
                  new_route = (
                          routes[r1_idx].route[None:r1_i+1:1]+
                          routes[r2_idx].route[r2_i:None:-1])
                  new_cost = (D[a, c],
                          routes[r1_idx].fwd_l[r1_i]+
                          routes[r2_idx].fwd_l[r2_i])
                  new_demand = (
                          routes[r1_idx].fwd_d[r1_i]+
                          routes[r2_idx].fwd_d[r2_i])
              else:
                  assert False, "there is no move to combine the tail of"+\
                                "route1 with the (reversed) head of route2"
                              
              ret.append(RouteData(new_route, new_cost, new_demand))
          ret.append(best_delta)
          return tuple(ret)
       
      return None, None, None, None
  
  
  @routeordersensitive
  def do_1point_move(route1_data, route2_data, D, d=None,
                        C=None, L=None, # constraints
                        strategy=LSOPT.FIRST_ACCEPT,
                        best_delta = None):
      """ Move one point from route1 to route2. Tries all possible combinations
      of moving a node from route 1 to any valid position in route 2. Sometimes
      called "relocate" (e.g., Bräysy & M. Gendreau 2005, Savelsbergh 1992), but
      here we use the name "one point move" (Groër et al 2010) to differentiate
      from the intra-route version.
      
      If an improving move was found and made, operation returns new routes in
      same format as route_data inputs, but if C and d are not given the 3.
      field of the tuple is None. If there is no improving move, returns None.
  
      Groër, C., Golden, B. and Wasil, E., 2010. A library of local search
       heuristics for the vehicle routing problem. Mathematical Programming
       Computation, 2(2), pp.79-101.
      Bräysy, O. & Gendreau, M. 2005. Vehicle Routing Problem, Part I: Route 
       Construction and Local Search Algorithms Transportation Science 39(1),
       pp. 104–118
      Savelsbergh, M. W. P. 1992. The vehicle routing problem with time windows:
       Minimizing route duration. J.Comput. 4 146-154
      """
  
      if route1_data==route2_data:
          return None,None,None
      
      # unpack route, current cost, and current demand
      route1, r1_l, r1_d, _ = route1_data
      route2, r2_l, r2_d, _ = route2_data
      
      if not best_delta:
          best_delta = 0
      best_move = None
      accept_move = False
      
      for i in xrange(1,len(route1)-1):
          remove_after = route1[i-1]
          to_move = route1[i]
          remove_before = route1[i+1]
          
          remove_delta = D[remove_after,remove_before]\
                        -D[remove_after,to_move]\
                        -D[to_move,remove_before]
          
          # capacity constraint feasibility check 
          if C and r2_d+d[to_move]-C_EPS>C:
              continue
  
          for j in xrange(1,len(route2)):
              insert_after = route2[j-1]
              insert_before = route2[j]
              
              
              insert_delta = D[insert_after,to_move]+D[to_move,insert_before]\
                            -D[insert_after,insert_before]
              
              # route cost constraint feasibility check 
              if L and r2_l+insert_delta-S_EPS>L:
                  continue
              
              delta = remove_delta+insert_delta
              if delta+S_EPS<best_delta:                    
                  best_delta = delta
                  best_move = (i, j, remove_delta, insert_delta)
                  if strategy==LSOPT.FIRST_ACCEPT:
                      accept_move=True
                      # break j-loop
                      break
          # break i-loop
          if accept_move:
              break
          
      if best_move:
          # unpack best move
          i, j, remove_delta, insert_delta = best_move
          to_move = route1[i] 
          return (RouteData(route1[:i]+route1[i+1:], r1_l+remove_delta,
                   None if not C else r1_d-d[to_move]),
                  RouteData(route2[:j]+[to_move]+route2[j:], r2_l+insert_delta,
                   None if not C else r2_d+d[to_move]),
                  remove_delta+insert_delta)
                  
      return None,None,None
  
  
  def do_2point_move(route1_data, route2_data, D, d=None,
                     C=None,L=None, # constraints
                     strategy=LSOPT.FIRST_ACCEPT,
                     best_delta = None):                 
      """ Swap one point from route1 with one point on route2 if it improves the
      solution. This operation is sometimes referred to as "exchange" (e.g. in 
      Bräysy & M. Gendreau 2005, Savelsbergh 1992), but we use the name
      "two point move" (Groër et al 2010) to differentiate it from the inter-
      route one.
      
      If an improving move was found and made, operation returns new routes in
       same format as route_data inputs, but if C and d are not given the 3.
       field of the tuple is None. If there is no improving move, returns None.
  
      Groër, C., Golden, B. and Wasil, E., 2010. A library of local search
       heuristics for the vehicle routing problem. Mathematical Programming
       Computation, 2(2), pp.79-101.
      Bräysy, O. & Gendreau, M. 2005. Vehicle Routing Problem, Part I: Route 
       Construction and Local Search Algorithms Transportation Science 39(1),
       pp. 104–118
      Savelsbergh, M. W. P. 1992. The vehicle routing problem with time windows:
       Minimizing route duration. J.Comput. 4 146-154
      """
      
      # unpack route, current cost, and current demand
      route1, r1_l, r1_d, _ = route1_data
      route2, r2_l, r2_d, _ = route2_data
      
      if not best_delta:
          best_delta = 0
      best_move = None
      accept_move = False
      
      for i in xrange(1,len(route1)-1):
          to_swap1 = route1[i]
  
          swap1_after = route1[i-1]
          swap1_before = route1[i+1]     
  
  
          for j in xrange(1,len(route2)-1):
              to_swap2 = route2[j]
              
              # capacity constraint feasibility check 
              if C and (r1_d-d[to_swap1]+d[to_swap2]-C_EPS>C or
                        r2_d-d[to_swap2]+d[to_swap1]-C_EPS>C):
                  continue
              
              swap2_after = route2[j-1]
              swap2_before = route2[j+1]    
              
              route1_delta = -D[swap1_after,to_swap1]\
                             -D[to_swap1,swap1_before]\
                             +D[swap1_after,to_swap2]\
                             +D[to_swap2,swap1_before]
              route2_delta = -D[swap2_after,to_swap2]\
                             -D[to_swap2,swap2_before]\
                             +D[swap2_after,to_swap1]\
                             +D[to_swap1,swap2_before]
                             
              # route cost constraint feasibility check 
              if L and (r1_l+route1_delta>L or r2_l+route2_delta>L):
                  continue
                          
              delta = route1_delta+route2_delta
              if delta+S_EPS<best_delta:                    
                  best_delta = delta
                  best_move = (i, j, route1_delta, route2_delta)
                  if strategy==LSOPT.FIRST_ACCEPT:
                      accept_move=True
                      break # j loop
          if accept_move:
              break # i loop
          
      if best_move:
          # unpack best move
          i, j, route1_delta, route2_delta = best_move
          to_swap1 = route1[i] 
          to_swap2 = route2[j] 
          return (RouteData(route1[:i]+[to_swap2]+route1[i+1:],
                   r1_l+route1_delta,
                   None if not C else r1_d-d[to_swap1]+d[to_swap2]),
                  RouteData(route2[:j]+[to_swap1]+route2[j+1:],
                   r2_l+route2_delta,
                   None if not C else r2_d-d[to_swap2]+d[to_swap1]),
                            
                  route1_delta+route2_delta)
                  
      return None,None,None
  
  @routeordersensitive
  def do_insert_move(unrouted, recieving_route_data, D,d=None,C=None,L=None,
                     strategy=LSOPT.FIRST_ACCEPT,
                     best_delta = None):
      """ Try to insert (unrouted) node(s) on the existing recieving route. The
      operation succeeds only if all customers can be inserted on the recieving
      route. The difference to the one point move is that there is no delta for
      removing the unrouted node from its RouteData / list.
      
      Returns a 3-tuple: An empty RouteData, A copy of the updated RouteData and
       the route cost (usually route length) delta; or (None,None,None) if the
       insertion fails.
      
      * unrouted
         Can be a single customer (int), a list of those, or a route (RouteData).
         If RouteData is given its .demand field must be set. 
      * recieving_route_data
         is the route (RouteData) that the customers to insert are inserted to.
         Its .demand field must be set. 
      * D,d,C,L
         the problem definition and constraints
      * strategy
         is the applied loca search strategy to use. Please note that if there is
         more than one customer to insert, the maximum route cost constraint L
         is set, and strategy is set to first accept, the best insertion position
         for the customers are still searched for to leave the most amount of
         room for subsequent insertions.
      """
  
      if type(unrouted) is int:
          unrouted_demand = d[unrouted] if C else 0
          unrouted = [unrouted]
      if isinstance(unrouted, RouteData):
          unrouted_demand = unrouted.demand if C else 0
          unrouted = unrouted.route[1:-1]
      elif C:
          unrouted_demand = sum(d[n] for n in unrouted)
      
      # make an ansatz for fitting the inserted customers into
      ansatz_route = list(recieving_route_data.route)
      ansatz_d = recieving_route_data.demand
      ansatz_l = recieving_route_data.cost
      total_delta = 0
      
      # is not possible because of constraint C
      if C and ansatz_d+unrouted_demand-C_EPS>C:
          return None, None, None
          
      for ni, node in enumerate(unrouted):
          best_insert_pos = None
          if strategy==LSOPT.BEST_ACCEPT or L:
              # need to find a place where it can be inserted (if at all)
              best_insert_delta = None
              for i in xrange(1, len(ansatz_route)):
                  insert_after = ansatz_route[i-1]
                  insert_before = ansatz_route[i]
                  
                  insert_delta = +D[insert_after,node]\
                                 +D[node,insert_before]\
                                 -D[insert_after,insert_before]
                          
                  if not L or ansatz_l+insert_delta-S_EPS<=L:
                      if (best_insert_delta is None) or\
                         (insert_delta<best_insert_delta):
                          best_insert_pos = i
                          best_insert_delta = insert_delta 
                      # to avoid complexity, only allow early termination with 
                      #  FIRST_ACCEPT when customers are inserted one
                      #  by one
                      if strategy==LSOPT.FIRST_ACCEPT and ni==len(unrouted)-1:
                          break # i loop
          else:
              # we know it fits, just append it 
              best_insert_pos = -1
              insert_after = ansatz_route[-2] # the last non-depot node
              best_insert_delta = +D[insert_after,node]+D[node,0]-D[insert_after,0] 
          
          # no L valid routing
          if best_insert_pos is None:
              return None, None, None
          else:
              if C: ansatz_d+=d[node]
              ansatz_l+=best_insert_delta
              # we found a valid insertion location accept it
              ansatz_route = ansatz_route[:best_insert_pos]\
                             +[node]\
                             +ansatz_route[best_insert_pos:]
              total_delta+=best_insert_delta
              
              # must be better than the preset level (if applicable)
              if (best_delta is not None) and (total_delta+S_EPS>best_delta):
                  return None, None, None
              
      return RouteData(), RouteData(ansatz_route, ansatz_l, ansatz_d), total_delta
  
  @routeordersensitive
  def do_redistribute_move(redisributed_route_data,
                           receiving_route_or_routes_data,
                           D,d=None, C=None, L=None,
                           strategy=LSOPT.FIRST_ACCEPT,
                           best_delta = None,
                           recombination_level=0):
      """
      Try to insert the nodes of the first route on the other routes if possible.
      Note that second argument receiving_route_data can also be a list of routes.
      
      The insertion order matters, so it is possible to try all permutations
      of nodes to be inserted and routes to insert to (and return the best).
      However, this will be very computationally intesive as the number of
      combinations explode. Therefore, there are different options available:
          
      * recombination_level=0, insert the nodes in the order they are on the
         route that will be destroyed and consider the recieving routes in the
         order of the "routes" parameter.
      * recombination_level=1, the nodes are inserted in all possible 
         order. The recieving routes are checked in the order of the "routes"
         parameter.
      * recombination_level=2, insert the nodes in the order they are on the
         route that will be destroyed. However, every possible order of the 
         recieving routes is checked.
      * recombination_level=3, all possible orderings of the nodes to be
         inserted AND the recieving routes is checked.       
         
      There could be yet another level (try_all_insertions_level=4), where for
      each insertion of each insertion node ordering, but it has not been
      implemented.
      """
   
      # allow only single redistributing to a single other route e.g. when using
      #  (local_search API)
      if isinstance(receiving_route_or_routes_data, RouteData):
          receiving_route_or_routes_data = [receiving_route_or_routes_data]
      
      # before iterating the permuations, check that it is at all possible to 
      #  fit all from the route to be destroyed to the other routes.
      if C:
          total_d_slack = sum(C-rrd.demand for rrd in receiving_route_or_routes_data)
          if redisributed_route_data.demand-C_EPS > total_d_slack:
              return [None]*(len(receiving_route_or_routes_data)+2)
          
      if not best_delta:
          # allows move that makes the overall solution cost worse
          best_delta = float("inf")
      redistribution_succesfull = False    
      best_insertion_routes = None
   
      # unpack route, current cost, and current demand
      route1, r1l, r1d, _ = redisributed_route_data    
      if recombination_level==0:
          insertion_node_orderings = [(route1[1:-1], receiving_route_or_routes_data)]
      elif recombination_level==1:
          insertion_node_orderings = product(
              permutations(route1[1:-1]), 
              [receiving_route_or_routes_data])
      elif recombination_level==2:
          insertion_node_orderings = product(
              [route1[1:-1]], 
              permutations(receiving_route_or_routes_data))
      elif recombination_level==3:
          insertion_node_orderings = product(
              list(permutations(route1[1:-1])), 
              list(permutations(receiving_route_or_routes_data)))
      
      #print(list(po for po in insertion_node_orderings))
      for route_customers, candidate_routes in insertion_node_orderings:
          total_delta = 0
          all_succesfully_inserted = True        
          ansatz_route_data = list(candidate_routes)
          
          for to_insert in route_customers:            
              node_succesfully_inserted = False
              for rd2_ar_idx, rd2 in enumerate(ansatz_route_data):
                  _, new_rd2, delta = do_insert_move(to_insert, rd2, D,d,C,L,
                                                     strategy=strategy)
                  if delta is None:
                      continue
                  
                  #print("REMOVEME:", "succesfully inserted %d to create %s"%(to_insert, new_rd2.route))
                  ansatz_route_data[rd2_ar_idx] = new_rd2
                  total_delta+=delta
                  if total_delta+S_EPS<best_delta:
                      node_succesfully_inserted = True
                      break
                  else:
                      continue
                  
              if not node_succesfully_inserted:
                  all_succesfully_inserted = False
                  break
              
          if all_succesfully_inserted:
              #print("REMOVEME:", "all inserted to create %s (delta %.2f)"%([rd.route for rd in ansatz_route_data], total_delta))   
              if total_delta<best_delta:
                  best_delta = total_delta
                  best_insertion_routes = ansatz_route_data
                  redistribution_succesfull = True
                      
      if redistribution_succesfull:
          # return an empty route and n routes the customers were distributed to
          return tuple( [RouteData()]+best_insertion_routes+[best_delta] )
      else:
          return [None]*(len(receiving_route_or_routes_data)+2)
          
  @routeordersensitive
  def do_chain_move(route1_data, route2_data, route3_data, D, d=None,
                        C=None, L=None, # constraints
                        strategy=LSOPT.FIRST_ACCEPT,
                        best_delta = None):  
      """ This is the "pair" operation described in Wren and Holliday (1972). 
      It involves moving a repacing a node on route 2 with a node on route 1.
      The replaced node is then inserted on route 3 (if able). Route 1!=2!=3 
      
      This is a very expensive operation, corresponding to 5-opt with chain
      length of 1, use with care.
      
      Wren, A. and Holliday, A., 1972. Computer scheduling of vehicles from one
      or more depots to a number of delivery points. Journal of the Operational
      Research Society, 23(3), pp.333-344.
      """
  
      if not best_delta:
          best_delta = 0
      best_move = None
      accept_move = False
      
      route1, r1_l, r1_d, _ = route1_data    
      route2, r2_l, r2_d, _ = route2_data
      route3, r3_l, r3_d, _ = route3_data
      
      for i in xrange(1,len(route1)-1):
          
          remove_after = route1[i-1]
          to_move = route1[i]
          remove_before = route1[i+1]
          
          remove_delta =  +D[remove_after,remove_before]\
                          -D[remove_after,to_move]\
                          -D[to_move,remove_before]
          
          for j in xrange(1,len(route2)-1):
              
              replace_after = route2[j-1]
              to_replace = route2[j]
              replace_before = route2[j+1]
              
              # can do capacity constraint feasibility check here
              if C and (r2_d-d[to_replace]+d[to_move]>C or
                        r3_d+d[to_replace]>C): 
                  continue
                              
              replace_delta = +D[replace_after,to_move]\
                              +D[to_move,replace_before]\
                              -D[replace_after,to_replace]\
                              -D[to_replace,replace_before]
                              
              if L and r2_l+replace_delta>L:
                  continue
          
              for k in xrange(1,len(route3)):
                  
                  insert_after = route3[k-1]
                  insert_before = route3[k]
                  
                  insert_delta =  +D[insert_after, to_replace]\
                                  +D[to_replace, insert_before]\
                                  -D[insert_after,insert_before]
                                  
                  if L and (r3_l+insert_delta>L):
                      continue
  
                  # check if this is best so far            
                  delta = remove_delta+replace_delta+insert_delta
                  if delta+S_EPS<best_delta:                    
                      best_delta = delta
                      best_move = (i, j, k, remove_delta, replace_delta, insert_delta)
                      if strategy==LSOPT.FIRST_ACCEPT:
                          accept_move=True
                          break # k loop
              if accept_move:
                  break # j loop
          if accept_move:
              break # i loop
              
      if best_move:
          # unpack best move
          i, j, k, remove_delta, replace_delta, insert_delta = best_move
          to_move = route1[i] 
          to_replace = route2[j] 
          return (RouteData(route1[:i]+route1[i+1:],
                      r1_l+remove_delta,
                      None if not C else r1_d-d[to_move]),
                  RouteData(route2[:j]+[to_move]+route2[j+1:],
                      r2_l+replace_delta,
                      None if not C else r2_d-d[to_replace]+d[to_move]),
                  RouteData(route3[:k]+[to_replace]+route3[k:],
                      r3_l+insert_delta,
                      None if not C else r3_d+d[to_replace]),
                  best_delta)
      return None,None,None,None             
  