Thursday, March 30, 2017

2-opt with Neighbor Lists and Don't Look Bits

proc One_City_2_Opt_ND(tour: var Tour_Array;
                       basePos: Tour_Index;
                       neighbor: Neighbor_Matrix;
                       position: var City_Position_In_Tour;
                       numberOfNeigbors: int;
                       DontLook: var DLB_Array): bool =
  var
    improved: bool
    pos_Y2, i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number

  improved = false
  block two_loops:
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
      else:
        i  = (N + basePos - 1) mod N
        X2 = tour[i]
        X1 = tour[(i+1) mod N]  # == tour[basePos]

      for neighbor_number in 0 .. numberOfNeigbors-1:
        Y1 = neighbor[X1][neighbor_number]
        if direction == forward:
          j = position[Y1]
          Y2 = tour[(j+1) mod N]
        else:
          j  = (N + position[Y1] - 1) mod N  # pos[Y1] == j+1
          Y2 = tour[j]

        if (X2 == Y1) or (Y2 == X1):
          continue

        if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
          Set_DLB_off(DontLook, [X1, X2, Y1, Y2])
          Make_2_Opt_Move(tour, i, j, position)
          improved = true
          break two_loops
  result = improved
proc LS_2_Opt_ND(tour: var Tour_Array;
                 neighbor: Neighbor_Matrix;
                 neighborListLen: int = N-1) =
  ## Optimizes the given tour using 2-opt with neighbor lists
  ## and dDon't Look Bits (DLB)
  var
    locallyOptimal: bool = false
    improved: bool
    DontLook: DLB_Array
    baseCity: City_Number
    position: City_Position_In_Tour

  let numberOfNeigbors = min(neighborListLen, len(neighbor[0]))
  position = Create_Position_In_Tour(tour)

  Set_DLB_off(DontLook, tour)
  
  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      baseCity = tour[basePos]
      if isDLB_on(DontLook, baseCity):
        continue
      improved = One_City_2_Opt_ND(tour, basePos, neighbor, position,
                                   numberOfNeigbors, DontLook)
      if not improved:
        Set_DLB_on(DontLook, baseCity)
      else:
        locallyOptimal = false

Wednesday, March 29, 2017

2-opt with Neighbor Lists and Fixed Radius

proc One_City_2_Opt_NR(tour: var Tour_Array;
                       basePos: Tour_Index;
                       neighbor: Neighbor_Matrix;
                       position: var City_Position_In_Tour;
                       numberOfNeigbors: int): bool =
  var
    improved: bool
    pos_Y2, i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number
    radius: Length

  improved = false
  block two_loops:
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
      else:
        i  = (N + basePos - 1) mod N
        X2 = tour[i]
        X1 = tour[(i+1) mod N]  # == tour[basePos]
      radius = distance(X1, X2)

      for neighbor_number in 0 .. numberOfNeigbors-1:
        Y1 = neighbor[X1][neighbor_number]
        if direction == forward:
          j = position[Y1]
          Y2 = tour[(j+1) mod N]
        else:
          j  = (N + position[Y1] - 1) mod N  # pos[Y1] == j+1
          Y2 = tour[j]

        if (X2 == Y1) or (Y2 == X1):
          continue
          
        if distance(X1, Y1) > radius:
          break

        if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
          Make_2_Opt_Move(tour, i, j, position)
          improved = true
          break two_loops
  result = improved
proc LS_2_Opt_NR(tour: var Tour_Array;
                 neighbor: Neighbor_Matrix;
                 neighborListLen: int = N-1) =
  ## Optimizes tour using 2-opt with neighbor lists and fixed radius
  var
    locallyOptimal: bool = false
    improved: bool
    position: City_Position_In_Tour

  let numberOfNeigbors = min(neighborListLen, len(neighbor[0]))
  position = Create_Position_In_Tour(tour)

  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      improved = One_City_2_Opt_NR(tour, basePos, neighbor,
                                   position, numberOfNeigbors)
      if improved:
        locallyOptimal = false

Monday, March 27, 2017

2-opt with Neighbor Lists

proc One_City_2_Opt_N(tour: var Tour_Array;
                      basePos: Tour_Index;
                      neighbor: Neighbor_Matrix;
                      position: var City_Position_In_Tour;
                      numberOfNeigbors: int): bool =
  ## Optimizes the given tour using 2-opt with neighours lists
  # For each edge looks for improvement considering neigbours
  # in ascending order of their distance
  var
    improved: bool
    pos_Y2, i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number

  improved = false
  block two_loops:
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
      else:
        i  = (N + basePos - 1) mod N
        X2 = tour[i]
        X1 = tour[(i+1) mod N]  # == tour[basePos]

      for neighbor_number in 0 .. numberOfNeigbors-1:
        # after 2-opt move new tour would contain direct link X1-Y1,
        # so we look for city Y1 among cities that are close to city X1
        Y1 = neighbor[X1][neighbor_number]
        if direction == forward:
          j = position[Y1]
          Y2 = tour[(j+1) mod N]
        else:
          j  = (N + position[Y1] - 1) mod N  # pos[Y1] == j+1
          Y2 = tour[j]

        if (X2 == Y1) or (Y2 == X1):
          # 2-opt is for non-adjacent links only
          # by constuction X1 != Y1
          continue

        if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
          Make_2_Opt_Move(tour, i, j, position)
          improved = true
          break two_loops
  result = improved
proc LS_2_Opt_N(tour: var Tour_Array;
                       neighbor: Neighbor_Matrix;
                       neighborListLen: int = N-1) =
  ## Optimizes the given tour using 2-opt with neighbor lists
  var
    locallyOptimal: bool = false
    improved: bool
    position: City_Position_In_Tour

  let numberOfNeigbors = min(neighborListLen, len(neighbor[0]))
  position = Create_Position_In_Tour(tour)

  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      improved = One_City_2_Opt_N(tour, basePos, neighbor,
                                  position, numberOfNeigbors)
      if improved:
        locallyOptimal = false

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by 2-opt and 2-opt with neighbor lists (number of neighbors in parentheses). Random planar problems, N=100.
Problem # Method Time Best Avg Worst
1NN100100.00113.18126.03
NN + 2opt125392.3095.93101.34
NN + 2opt-N(99)93992.9895.97101.72
NN + 2opt-N(48)52092.8995.89101.72
NN + 2opt-N(24)31392.8995.56101.72
2NN100100.00106.57115.38
NN + 2opt87586.0290.9896.83
NN + 2opt-N(99)107486.1891.2396.04
NN + 2opt-N(48)58786.1891.2494.55
NN + 2opt-N(24)29386.9991.2395.95
3NN100100.00106.76115.62
NN + 2opt156682.1084.8792.41
NN + 2opt-N(99)93382.8085.4190.72
NN + 2opt-N(48)52683.0585.7091.73
NN + 2opt-N(24)30682.1085.6690.70
4NN100100.00108.63115.55
NN + 2opt166688.4891.5798.65
NN + 2opt-N(99)94088.6392.4697.13
NN + 2opt-N(48)52088.6392.6998.92
NN + 2opt-N(24)20688.3792.6697.79
5NN100100.00107.56119.38
NN + 2opt104691.2794.1798.35
NN + 2opt-N(99)93390.9994.80100.14
NN + 2opt-N(48)41991.8594.7398.68
NN + 2opt-N(24)20690.7094.69100.00
6NN100100.00110.60121.69
NN + 2opt156687.9690.8694.66
NN + 2opt-N(99)93387.7992.3396.87
NN + 2opt-N(48)42088.0992.2896.82
NN + 2opt-N(24)31387.7391.8496.76
7NN100100.00111.52123.61
NN + 2opt114690.0694.4597.38
NN + 2opt-N(99)83391.0094.7298.71
NN + 2opt-N(48)41990.6994.6797.86
NN + 2opt-N(24)31390.6494.5297.86
8NN100100.00111.82121.50
NN + 2opt136086.5790.2894.59
NN + 2opt-N(99)104086.7790.2895.18
NN + 2opt-N(48)51986.6490.0994.81
NN + 2opt-N(24)31386.7789.8695.55
9NN100100.00116.59127.95
NN + 2opt146090.1995.30104.15
NN + 2opt-N(99)104090.8296.62105.77
NN + 2opt-N(48)62690.8296.74105.77
NN + 2opt-N(24)31390.8296.60105.77
10NN100100.00108.10118.75
NN + 2opt88189.6894.3598.28
NN + 2opt-N(99)87489.9694.6497.56
NN + 2opt-N(48)48791.8494.7399.41
NN + 2opt-N(24)20091.2394.2399.17

Position of given city in tour

Using neighbor lists forces us to use some convenient and quick method to determine for given city (a neighbor) what is its position in current tour. We use an array indexed by city numbers and storing corresponding positions.

type
  City_Position_In_Tour = array[City_Number, Tour_Index]
proc Create_Position_In_Tour(tour: Tour_Array): City_Position_In_Tour =
  ## Creates array with position of each city in tour
  var
    position: City_Position_In_Tour

  for pos in 0 .. N-1:
    position[tour[pos]] = pos
  result = position

Note that array of positions should be updated during any action that changes the tour, be it reversing a segment:

proc Reverse_Segment(tour: var Tour_Array;
                     startIndex, endIndex: Tour_Index;
                     position: var City_Position_In_Tour) =
  ## Reverses order of elements in segment [startIndex, endIndex]
  ## of tour, and updates positon[] array used with neighbour lists
  var
    left, right: Tour_Index

  let inversionSize = ((endIndex - startIndex + 1 + N) mod N) div 2
 
  left  = startIndex
  right = endIndex

  for counter in 1 .. inversionSize:
    swap(tour[left], tour[right])

    position[tour[left]]  = left
    position[tour[right]] = right

    left  = (left + 1) mod N
    right = (N + right - 1) mod N

...or shifting a segment:

proc Shift_Segment(tour: var Tour_Array; i, j, k: Tour_Index;
                   position: var City_Position_In_Tour) =
  ## Shifts the segment of tour:
  ## cities from t[i+1] to t[j] from their current position to position
  ## after current city t[k], that is between cities t[k] and t[k+1].
  ## Accordingly updates positon[] array.
  # Assumes:  k, k+1 are not within the segment [i+1..j]
  let
    segmentSize = (j - i + N) mod N
    shiftSize = ((k - i + N)  - segmentSize + N) mod N
    offset = i + 1 + shiftSize
  var
    pos, posNew: Tour_Index
    segment: seq[City_Number] = newSeq[City_Number](segmentSize)

  # make a copy of the segment before shift
  for pos in 0 .. segmentSize-1:
    segment[pos] = tour[(pos + i+1) mod N]

  # shift to the left by segmentSize all cities between old position
  # of right end of the segment and new position of its left end
  pos  = (i + 1) mod N
  for counter in 1 .. shiftSize:
    tour[pos] = tour[(pos + segmentSize) mod N]
    position[tour[pos]]  = pos
    pos  = (pos + 1) mod N
  #end_for counter loop

  # put the copy of the segment into its new place in the tour
  for pos in 0 .. segmentSize-1:
    posNew = (pos + offset) mod N 
    tour[posNew] = segment[pos]    
    position[tour[posNew]] = posNew

Sunday, March 26, 2017

Nearest neighbor lists

Note that neighbor lists are used to speed up the process of optimization, and using them may cost slightly worse quality of resulting tours, even if it not always occurs.

type
  Neighbor_Matrix = array [City_Number, seq[City_Number]]
proc Build_Neighbors_Matrix(No_Of_Neigbors: int = N-1): Neighbor_Matrix =
  ## Builds array of No_Of_Neigbors nearest neighbors of each city
  # neighbor[i][j] is j-th nearest neighbour of city number i.
  var
    currentCity, otherCity: City_Number
    idx: int
    neighbors_of_currentCity: seq[City_Number]
    neighbor: Neighbor_Matrix

  newSeq(neighbors_of_currentCity, N-1)
  for currentCity in 0 .. N-1:
    # create sequence of all neighbors of current city:
    idx = 0
    for otherCity in 0 .. N-1:
      if otherCity != currentCity:
        neighbors_of_currentCity[idx] = otherCity
        idx = idx + 1

    # Sort elements in neighbors_of_currentCity by ascending distance
    # to currentCity (use any convenient method of sorting), then take
    # first No_Of_Neigbors cities from the result.
    sort(neighbors_of_currentCity) do (x,y: City_Number) -> int:
      result = cmp(distance(currentCity, x), distance(currentCity, y))

    neighbor[currentCity] = @[]
    setlen(neighbor[currentCity], No_Of_Neigbors)
    for idx in 0 .. No_Of_Neigbors-1:
      neighbor[currentCity][idx] = neighbors_of_currentCity[idx]
  result = neighbor

2-opt with DLB and Fixed Radius - alternative

While the inequity which a 2-opt move must meet to improve a tour::

holds if at least one of the two conditions is met:

or

The inequity (a) means that for given (X1, X2) we can search for Y2 among those cities that are closer to X2 than is X1.

But alternatively the inequity between the two sums may be considered as a pair of inequities:

or

The inequity (c) means that for given (X1, X2) we can search for Y1 among those cities that are closer to X1 than is X2.

In fixed radius search we can use either (a) or (c), the radius is the same. When fixed radius searching is used with DLB, then probably searching around X1 for Y1 would be better – because we would rather prefer to not miss a good move for base city (X1) before activating its don’t look bit. When we want to continue searching and moving in given direction, then probably searching for Y2 around X2 could be better idea.

proc One_City_2_Opt_DLBR2(tour: var Tour_Array,
                          basePos: Tour_Index,
                          DontLook: var DLB_Array): bool =
  var
    improved: bool
    i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number
    radius: Length

  improved = false
  block two_loops:
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
      else:
        i  = (N + basePos - 1) mod N
        X2 = tour[i]
        X1 = tour[basePos]  # == tour[(i+1) mod N]
      radius = distance(X1, X2)

      for counter_2 in 0 .. N-1:
        j = counter_2
        Y1 = tour[j]
        Y2 = tour[(j+1) mod N]
        if direction == backward:
          swap(Y1, Y2)
        if (X1 == Y1) or (X2 == Y1) or (Y2 == X1):
          continue

        if distance(X1, Y1) > radius:  # NOTE
          continue

        if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
          Set_DLB_off(DontLook, [X1, X2, Y1, Y2])
          Make_2_Opt_Move(tour, i, j)
          improved = true
          break two_loops
  result = improved

2-opt with DLB and Fixed Radius

Note

This example of code does not contain any new ideas and can be skipped by a reader.

proc One_City_2_Opt_DLBR(tour: var Tour_Array,
                         basePos: Tour_Index,
                         DontLook: var DLB_Array): bool =
  var
    improved: bool
    i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number
    radius: Length

  improved = false
  block two_loops:
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
      else:
        i  = (N + basePos - 1) mod N
        X2 = tour[i]
        X1 = tour[basePos]  # == tour[(i+1) mod N]
      radius = distance(X1, X2)

      for counter_2 in 0 .. N-1:
        j = counter_2
        Y1 = tour[j]
        Y2 = tour[(j+1) mod N]
        if direction == backward:
          swap(Y1, Y2)
        if (X1 == Y1) or (X2 == Y1) or (Y2 == X1):
          continue

        if distance(X2, Y2) > radius:
          continue

        if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
          Set_DLB_off(DontLook, [X1, X2, Y1, Y2])
          Make_2_Opt_Move(tour, i, j)
          improved = true
          break two_loops
  result = improved
proc LS_2_Opt_DLBR(tour: var Tour_Array) =
  ## Optimizes the given tour using 2-opt with Don't Look Bits (DLB)
  ## and fixed radius search
  var
    locallyOptimal: bool = false
    improved: bool
    DontLook: DLB_Array
    baseCity: City_Number

  Set_DLB_off(DontLook, tour)

  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      baseCity = tour[basePos]
      if isDLB_on(DontLook, baseCity):
        continue
      improved = One_City_2_Opt_DLBR(tour, basePos, DontLook)
      if not improved:
        Set_DLB_on(DontLook, baseCity)
      else:
        locallyOptimal = false

Friday, March 24, 2017

2-opt with Fixed Radius Search

Fixed radius search uses the property of the inequity which a 2-opt move must meet to improve a tour. We can benefit from this even more, but only when we consider links in the two orientations and the two symmetric forms.

forward (ABCD)backward (DCBA)
A B C D A B C D
implementation
look forwardlook backward
base base + 1 j j + 1 j j + 1 base - 1 base
X1 X2 Y1 Y2 Y2 Y1 X2 X1

Using this scheme we cannot miss an improving move even if we restrict possible choices to one inequity.

Example implementation:

type
  Orientation = enum forward, backward
proc One_City_2_Opt_R(tour: var Tour_Array,
                         basePos: Tour_Index): bool =
  var
    improved: bool
    i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number
    radius: Length

  improved = false
  block two_loops:
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
      else:
        i  = (N + basePos - 1) mod N
        X2 = tour[i]
        X1 = tour[basePos]  # == tour[(i+1) mod N]
      radius = distance(X1, X2)

      for counter_2 in 0 .. N-1:
        j = counter_2
        Y1 = tour[j]
        Y2 = tour[(j+1) mod N]
        if direction == backward:
          swap(Y1, Y2)
        if (X1 == Y1) or (X2 == Y1) or (Y2 == X1):
          continue

        if distance(X2, Y2) > radius:
          continue

        if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
          Make_2_Opt_Move(tour, i, j)
          improved = true
          break two_loops
  result = improved
proc LS_2_Opt_R(tour: var Tour_Array) =
  ## Optimizes the given tour using 2-opt with fixed radius search
  var
    locallyOptimal: bool = false
    improved: bool

  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      improved = One_City_2_Opt_R(tour, basePos)
      if improved:
        locallyOptimal = false

Thursday, March 23, 2017

Fixed radius search – general idea

The condition

Suppose that for some X1, X2, Y1, Y2 occurs:

Hence:

Therefore the inequity which a 2-opt move must meet to improve a tour::

holds if at least one of the two conditions is met:

or

Simple radius search with 2-opt

Note that for simplicity of the code below we do not change orientation and instead use both conditions.

proc LS_2_Opt_simpleRadius(tour: var Tour_Array) =
  var
    locallyOptimal: bool = false
    i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number
    radius: Length

  while not locallyOptimal:
    locallyOptimal = true
    block two_loops:
      for counter_1 in 0 .. N-1:
        i  = counter_1
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
        radius = distance(X1, X2)

        for counter_2 in 0 .. N-1:
          j = counter_2
          Y1 = tour[j]
          Y2 = tour[(j+1) mod N]
          if (X1 == Y1) or (X2 == Y1) or (Y2 == X1):
            continue

          if distance(X2, Y2) > radius:
            if distance(X1, Y1) > radius:
              continue

          if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
            Make_2_Opt_Move(tour, i, j)
            locallyOptimal = false
            break two_loops

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by 2-opt and 2-opt with simple radius search algorithms. Random planar problems, N=100.
Problem # Method Time Best Avg Worst
1NN100100.00107.63117.95
NN + 2opt94086.9589.1894.09
NN + 2opt-simpleR72685.3788.8693.34
2NN100100.00107.02114.02
NN + 2opt78187.5290.9298.22
NN + 2opt-simpleR58787.9090.6495.50
3NN100100.00108.09121.54
NN + 2opt97487.5390.6796.83
NN + 2opt-simpleR58787.5391.0597.84
4NN100100.00109.31118.36
NN + 2opt87486.4791.6096.67
NN + 2opt-simpleR49386.8590.8897.71
5NN100100.00108.73120.07
NN + 2opt83385.7789.4496.25
NN + 2opt-simpleR62686.1289.3493.99
6NN100100.00106.37114.14
NN + 2opt83388.8391.2795.10
NN + 2opt-simpleR52687.2890.6795.23
7NN100100.00108.36116.70
NN + 2opt107482.5186.7391.88
NN + 2opt-simpleR68182.8186.5391.64
8NN100100.00109.35120.17
NN + 2opt68195.2097.46103.04
NN + 2opt-simpleR48794.3297.70102.40
9NN100100.00112.93125.15
NN + 2opt97590.5194.84100.73
NN + 2opt-simpleR58790.9795.1499.15
10NN100100.00110.36125.80
NN + 2opt83387.3892.3397.47
NN + 2opt-simpleR62687.8092.0297.50

Tuesday, March 21, 2017

3-opt with DLB

Procedures to calculate the gain from 3-opt move and to make 3-opt move assume that their parameters (cities and their indices in tour) are given in the order of ascending indices of tour array (modulo N). Therefore we need a function to recognize whether three indices are in this order, or: whether in this orientation the middle index is between the two boundary indices.

proc Between(a, x, b: Tour_Index): bool =
  ## Returns true if `x` is between `a` and `b` in cyclic sequence
  ## with established orientation of nondescending indices
  # That is: when one begins a forward traversal of the tour
  # at position `a', then position `x` is reached before position `b'.
  # Returns true if and only if:
  #   a < x < b  or  b < a < x  or  x < b < a
  if (b > a):
    result = (x > a) and (x < b)
  elif (b < a):
    result = (x > a) or (x < b)
  else:
    result = false

Note the additional speed-up: the link between base city X1 and its tour neighbor X2 is not considered when at least one of these cities has its don’t-look bit set to on.

proc One_City_3_Opt_DLB(tour: var Tour_Array,
                        basePos: Tour_Index,
                        DontLook: var DLB_Array): bool =
  type
    Move3 = object
      i, j, k: Tour_Index
      optCase: Reconnection_3_optCase
  var
    improved: bool
    i, j, k: Tour_Index
    X1, X2, Y1, Y2, Z1, Z2: City_Number
    gainExpected: Length_Gain
    optCase: Reconnection_3_optCase
    goodMove: Move3

  improved = false
  block three_loops:
    for direction in [forward, backward]:
      if direction == forward:
        i  = basePos
      else:
        i  = (N + basePos - 1) mod N
      X1 = tour[i]
      X2 = tour[(i+1) mod N]

      if isDLB_on(DontLook,X1) or isDLB_on(DontLook,X2):  # NOTE
        continue

      for counter_2 in 0 .. N-1:
        j = counter_2 # second cut after j
        Y1 = tour[j]
        Y2 = tour[(j+1) mod N]

        if (X1 == Y1) or (X2 == Y1) or (Y2 == X1):
          continue

        # examine 2-opt move (see: opt3_case_1 .. opt3_case_3)
        gainExpected = Gain_From_2_Opt(X1, X2, Y1, Y2)
        if gainExpected > 0:
          goodMove = Move3(i: i,  j: j,  k: j,
                           optCase: opt3_case_1)
          improved = true
          break three_loops

        for counter_3 in 0 .. N-1:
          k = counter_3  # third cut after k
          Z1 = tour[k]
          Z2 = tour[(k+1) mod N]
          if (X1 == Z1) or (Y1 == Z1):
            continue

          if not Between(i, j, k):
            continue

          # examine pure 3-opt cases
          for optCase in [opt3_case_6, opt3_case_7]:
            gainExpected = Gain_From_3_Opt(X1, X2, Y1, Y2, Z1, Z2,
                                           optCase)
            if gainExpected > 0:
              goodMove = Move3(i: i,  j: j,  k: k,
                               optCase: optCase)
              improved = true
              break three_loops
      # end_block three_loops
  if improved:
    Set_DLB_off(DontLook, [X1, X2, Y1, Y2])
    if goodMove.optCase > opt3_case_3:
      Set_DLB_off(DontLook, [Z1, Z2])
    Make_3_Opt_Move(tour,
                    goodMove.i,  goodMove.j,  goodMove.k,
                    goodMove.optCase)
  result = improved
proc LS_3_Opt_DLB(tour: var Tour_Array) =
  ## Optimizes the given tour using 3-opt with Don't Look Bits (DLB)
  var
    locallyOptimal: bool = false
    improved: bool
    DontLook: DLB_Array
    baseCity: City_Number

  Set_DLB_off(DontLook, tour)

  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      baseCity = tour[basePos]
      if isDLB_on(DontLook, baseCity):
        continue
      improved = One_City_3_Opt_DLB(tour, basePos, DontLook)
      if not improved:
        Set_DLB_on(DontLook, baseCity)
      else:
        locallyOptimal = false

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by 3-opt and 3-opt with DLB algorithms. Random planar problems, N=100.
Problem # Method Time Best Avg Worst
1NN100100.00107.94117.16
NN + 2opt126890.0692.88101.40
NN + 3opt7509388.6189.9291.77
NN + 3opt-DLB2812588.6190.0792.72
2NN100100.00109.67122.62
NN + 2opt114688.9194.14101.71
NN + 3opt7364688.8990.4193.56
NN + 3opt-DLB2791988.9990.8594.98
3NN100100.00107.04118.01
NN + 2opt94091.1194.7498.25
NN + 3opt7292088.7791.1993.11
NN + 3opt-DLB2593388.8691.5393.85
4NN100100.00107.40118.32
NN + 2opt97589.0392.5697.91
NN + 3opt7285687.9989.5692.33
NN + 3opt-DLB2568187.9989.8093.41
5NN100100.00104.79113.38
NN + 2opt146081.1984.6089.37
NN + 3opt7270679.8281.7684.73
NN + 3opt-DLB2916679.8281.9885.00
6NN100100.00104.66115.85
NN + 2opt136082.9886.0390.66
NN + 3opt6853982.3983.5985.12
NN + 3opt-DLB2697982.3983.5985.56
7NN100100.00109.45118.42
NN + 2opt136088.1092.8897.84
NN + 3opt7863986.6688.8391.93
NN + 3opt-DLB2802686.6689.1492.78
8NN100100.00104.37112.02
NN + 2opt104688.5091.5997.77
NN + 3opt7531386.4088.3491.11
NN + 3opt-DLB2604086.7588.9295.75
9NN100100.00108.32118.64
NN + 2opt156691.3194.95100.77
NN + 3opt6646090.4692.0595.71
NN + 3opt-DLB2572690.5992.1295.83
10NN100100.00107.87122.61
NN + 2opt114689.1293.2199.11
NN + 3opt6906687.9889.4393.58
NN + 3opt-DLB2666687.9889.6792.78

Monday, March 20, 2017

2.5-opt with DLB

proc One_City_25_Opt_DLB(tour: var Tour_Array,
                         basePos: Tour_Index,
                         DontLook: var DLB_Array): bool =
  var
    improved: bool
    i, j: Tour_Index
    X1, X2, Y1, Y2, V0: City_Number

  improved = false
  block two_loops:
    for direction in [forward, backward]:  # NOTE: both tour neighbors
      # consider link between t[i] and t[i+1],
      # then between t[i-1] and t[i]
      if direction == forward:
        i  = basePos
      else:
        i  = (N + basePos - 1) mod N
      X1 = tour[i]
      X2 = tour[(i+1) mod N]

      for counter_2 in 0 .. N-1:
        j = counter_2
        Y1 = tour[j]
        Y2 = tour[(j+1) mod N]
        if (X1 == Y1) or (X2 == Y1) or (Y2 == X1):
          continue

        if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
          Make_2_Opt_Move(tour, i, j)
          improved = true
        else:
          V0 = tour[(i+2) mod N]
          if V0 != Y1:
            if Gain_From_Node_Shift(X1, X2, V0, Y1, Y2) > 0:
              Set_DLB_off(DontLook, V0)
              Make_Node_Shift_Move(tour, (i+1) mod N, j)
              improved = true
          else:
            V0 = tour[(N+j-1) mod N]
            if V0 != X2:
              if Gain_From_Node_Shift(V0, Y1, Y2, X1, X2) > 0:
                Set_DLB_off(DontLook, V0)
                Make_Node_Shift_Move(tour, j, i)
                improved = true
        #end_if Gain_From_2_Opt...
        if improved:
          Set_DLB_off(DontLook, [X1, X2, Y1, Y2])
          break two_loops
  result = improved
proc LS_25_Opt_DLB(tour: var Tour_Array) =
  ## Optimizes the given tour using 2.5-opt with Don't Look Bits (DLB)
  var
    locallyOptimal: bool = false
    improved: bool
    DontLook: DLB_Array
    baseCity: City_Number

  Set_DLB_off(DontLook, tour)

  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      baseCity = tour[basePos]
      if isDLB_on(DontLook, baseCity):
        continue
      improved = One_City_25_Opt_DLB(tour, basePos, DontLook)
      if not improved:
        Set_DLB_on(DontLook, baseCity)
      else:
        locallyOptimal = false

Saturday, March 18, 2017

2-opt with DLB

type
  DLB_Array = array [City_Number, bool]
  Orientation = enum forward, backward
proc isDLB_on(DontLook: DLB_Array,
              city: City_Number): bool =
  return DontLook[city]

proc Set_DLB_on(DontLook: var DLB_Array,
                city: City_Number) =
  DontLook[city] = true

proc Set_DLB_off(DontLook: var DLB_Array,
                 city: City_Number) =
  DontLook[city] = false

proc Set_DLB_off(DontLook: var DLB_Array,
                 cities: openArray[City_Number]) =
  # turns off DLB for given cities
  for city in 0 .. len(cities)-1:
    DontLook[cities[city]] = false
proc One_City_2_Opt_DLB(tour: var Tour_Array,
                        basePos: Tour_Index,
                        DontLook: var DLB_Array): bool =
  var
    improved: bool
    i, j: Tour_Index
    X1, X2, Y1, Y2: City_Number

  improved = false
  block two_loops:
    for direction in [forward, backward]:  # NOTE: both tour neighbors
      # consider link between t[i] and t[i+1],
      # then between t[i-1] and t[i]
      if direction == forward:
        i  = basePos
      else:
        i  = (N + basePos - 1) mod N
      X1 = tour[i]
      X2 = tour[(i+1) mod N]

      for counter_2 in 0 .. N-1:
        j = counter_2
        Y1 = tour[j]
        Y2 = tour[(j+1) mod N]

        if (X1 == Y1) or (X2 == Y1) or (Y2 == X1):
          continue

        if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
          Set_DLB_off(DontLook, [X1, X2, Y1, Y2])
          Make_2_Opt_Move(tour, i, j)
          improved = true
          break two_loops
  result = improved
proc LS_2_Opt_DLB(tour: var Tour_Array) =
  ## Optimizes the given tour using 2-opt with Don't Look Bits (DLB)
  var
    locallyOptimal: bool = false
    improved: bool
    DontLook: DLB_Array
    baseCity: City_Number

  Set_DLB_off(DontLook, tour)

  while not locallyOptimal:
    locallyOptimal = true
    for basePos in 0 .. N-1:
      baseCity = tour[basePos]
      if isDLB_on(DontLook, baseCity):
        continue
      improved = One_City_2_Opt_DLB(tour, basePos, DontLook)
      if not improved:
        Set_DLB_on(DontLook, baseCity)
      else:
        locallyOptimal = false

Friday, March 17, 2017

Don't Look Bits (DLB) – general idea

The idea

A local improvement algorithm can be made faster when we notice that the running time includes the time to needed perform the move, but also the time needed to examine possible good moves and find an improving move (if there is one).

If for a given city X1 we previously did not find any improving move and it still has the same neighbors (none of the links between this city and its neighbors has been changed), then chances that we will find an improving move for this city in current iteration are small (although non-zero, see below).

The concept based on this observation is known as don’t look bits (DLB). We use special flag for each of the cities. Initially all the flags are turned off, which means that we allow searching for improving moves starting from given city. When a search for improving move starting from city C fails, then the bit for this city is turned on. When a move is performed in which this city is involved as one of the endpoints of the exchanged links, then the bit for this city is turned off. When we consider candidates for X1, the first city of an improving move, we skip all cities which have their don’t-look bits set to on.

Comments

Consider also what happens with DLB when we examine moves “in one direction” only (see 2-opt basic algorithn). In such case we consider only one of the two links of X1, that is the link to its successor, the link between tour[i] and tour[i+1]. If we failed to found any improving move, then we turn the don’t-look bit for X1 on. If then in some subsequent iteration a move is performed in which X1 is one of the endpoints of the exchanged links, then we turn the don’t-look bit off. It works as expected. But what if a move changes the orientation of the segment in which X1 is inside, but is not at one of its ends? The don’t-look bit remains unchanged, although now the previous predecessor of X1 is its current successor, and we examined X1 only with its succcessor to verify that there are no improvements for X1. So now there could be some improving move for X1 in “forward” orientation, but it will not be examined, because the don’t-look bit for this city is still on.

That is why in the code below we consider both tour neighbors of a given city.

For the same reason of not missing improving moves we do not limit searching to the cases where j>i, but cosider all possible pairs of (i, j).

Aside from the above, it may seem impossible that new, improving move for given city could appear when both neighbors of this city remain unchanged, but let us take a closer look:

X0 Y1 Z1 X0_pred X0_succ Y2 Z2 no gain found from X0
X0 Y1 Z1 X0_pred X0_succ Y2 Z2 move with gain, not from X0
X0 Y1 Z1 X0_pred X0_succ Y2 Z2 after the move: possible gain from X0 with new links

However we should notice that in this case cities Y1 and Z1 have their don’t-look bits set to off after making the move. Therefore the new, possibly improving move: from X0 and involving the link between Y1 and Z1, will be examined in one of the subsequent iterations, that is as a move from Y1 or as a move from Z1.

Using Dont' Looks Bits with a simple 2-opt

type
  DLB_Array = array[City_Number, bool]
  Orientation = enum forward, backward
proc LS_2_Opt_DLB_General_Idea(tour: var Tour_Array) =
  ## Optimizes the given tour using 2-opt with Don't Look Bits (DLB) 
  var
    locallyOptimal: bool = false
    i, j: Tour_Index
    pos_2_Limit: Tour_Index
    X1, X2, Y1, Y2: City_Number
    DontLook: DLB_Array

  for city in 0 .. N-1:
    DontLook[city] = false

  while not locallyOptimal:
    locallyOptimal = true
    block three_loops:
      for counter_1 in 0 .. N-1:

        if DontLook[tour[counter_1]]:
          continue

        for direction in [forward, backward]:  # NOTE: both tour neighbors
          # consider link between t[i] and t[i+1],
          # then between t[i-1] and t[i]
          if direction == forward:
            i  = counter_1
          else:
            i  = (N + counter_1 - 1) mod N
          X1 = tour[i]
          X2 = tour[(i+1) mod N]

          for counter_2 in 0 .. N-1:
            j = counter_2  # 2nd cut between t[j] and t[j+1]
            Y1 = tour[j]
            Y2 = tour[(j+1) mod N]
            if (X1 == Y1) or (X2 == Y1) or (Y2 == X1):
              # the same city or adjacent cities
              continue

            if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
              DontLook[X1] = false
              DontLook[X2] = false
              DontLook[Y1] = false
              DontLook[Y2] = false
              Make_2_Opt_Move(tour, i, j)
              locallyOptimal = false
              break three_loops
        # end_for backward loop
        DontLook[tour[counter_1]] = true
      # end_for counter_1 loop

It is expected that the resulting tours will be different from the ones from 2-opt First Take procedure.

See also an alternative

Instead of don't look bits speedup we can use an array indexed by city numbers, in which each element contains length of tour after improving from given city. This would make implementation simpler, but would require more memory, because of storing numbers (lengths), which are longer than booleans.

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by 2-opt and 2-opt with DLBs algorithms. Random planar problems, N=100.
Problem # Method Time Best Avg Worst
1NN100100.00112.31122.36
NN + 2opt114693.0997.73104.89
NN + 2opt-DLB31394.4299.86104.98
2NN100100.00108.83124.17
NN + 2opt146286.1690.2097.26
NN + 2opt-DLB29387.8291.2696.06
3NN100100.00106.21115.53
NN + 2opt156682.1985.7789.56
NN + 2opt-DLB31382.6687.3992.24
4NN100100.00105.83112.36
NN + 2opt114690.0294.7898.19
NN + 2opt-DLB21391.6494.9199.78
5NN100100.00111.19119.81
NN + 2opt78191.9195.4399.81
NN + 2opt-DLB29391.8196.31102.66
6NN100100.00112.89126.03
NN + 2opt136094.3599.05104.02
NN + 2opt-DLB20696.0299.56105.29
7NN100100.00109.35117.75
NN + 2opt114690.4793.2899.21
NN + 2opt-DLB31390.5593.8099.71
8NN100100.00109.02120.22
NN + 2opt125386.5091.8898.31
NN + 2opt-DLB31388.4592.6495.84
9NN100100.00108.79116.77
NN + 2opt94095.4697.58101.63
NN + 2opt-DLB42095.5098.14102.77
10NN100100.00105.40113.11
NN + 2opt126882.1285.6689.87
NN + 2opt-DLB29381.7886.4090.22

Thursday, March 16, 2017

2.5-opt

2.5-opt algorithm, also known as 2h-opt, is an extension of 2-opt, that takes into account two simple alternative moves of 3-opt type. In its most internal loop it examines not only a 2-opt move, but also two Node Shifts:

before
X1 X2 Y1 Y2
after
variant Avariant Bvariant C
X1 X2 Y1 Y2 X1 X2 Y1 Y2 X1 X2 Y1 Y2

2.5-optimization

proc LS_25_Opt_Take_First(tour: var Tour_Array) =
  ## Optimizes the given tour using 2.5-opt
  # Shortens tour by repeating 2.5-opt moves until no improvement
  # can by done: in every iteration immediatelly makes permanent
  # the first move found that gives any length gain.
  var
    locallyOptimal: bool = false
    i, j: Tour_Index
    X1, X2, Y1, Y2, V0: City_Number

  while not locallyOptimal:
    locallyOptimal = true
    block two_loops:
      for counter_1 in 0 .. N-3:
        i = counter_1
        X1 = tour[i]
        X2 = tour[(i+1) mod N]
        for counter_2 in (counter_1+2) .. N-1:
          j = counter_2
          Y1 = tour[j]
          Y2 = tour[(j+1) mod N]

          if Gain_From_2_Opt(X1, X2, Y1, Y2) > 0:
            Make_2_Opt_Move(tour, i, j)
            locallyOptimal = false
            break two_loops
          else:
            V0 = tour[(i+2) mod N]
            if V0 != Y1:
              if Gain_From_Node_Shift(X1, X2, V0, Y1, Y2) > 0:
                Make_Node_Shift_Move(tour, (i+1) mod N, j)
                locallyOptimal = false
                break two_loops
            else:
              V0 = tour[(N+j-1) mod N]
              if V0 != X2:
                if Gain_From_Node_Shift(V0, Y1, Y2, X1, X2) > 0:
                  Make_Node_Shift_Move(tour, j, i)
                  locallyOptimal = false
                  break two_loops

Note that this code could be improved to be more efficient:

var
    del_2_Length: Length
    dX1Y1, dX2Y2, dX2Y1: Length
...
        for counter_2 in (counter_1+2) .. N-1:
...
          del_2_Length = distance(X1, X2) + distance(Y1, Y2)
          dX1Y1 = distance(X1, Y1)
          dX2Y2 = distance(X2, Y2)
          if del_2_Length - (dX1Y1 + dX2Y2) > 0:
            Make_2_Opt_Move(tour, i, j)
            locallyOptimal = false
            break two_loops
          else:
            dX2Y1 = distance(X2, Y1)
            V0 = tour[(i+2) mod N]
            if V0 != Y1:
              if (del_2_Length + distance(X2, V0)) - 
                 (dX2Y2 + dX2Y1 + distance(X1, V0)) > 0:
                Make_Node_Shift_Move(tour, (i+1) mod N, j)
                locallyOptimal = false
                break two_loops
            else:
              V0 = tour[(N+j-1) mod N]
              if V0 != X2:
                if (del_2_Length + distance(Y1, V0)) - 
                   (dX1Y1 + dX2Y1 + distance(Y2, V0)) > 0:
                  Make_Node_Shift_Move(tour, j, i)
                  locallyOptimal = false
                  break two_loops

This reduces running times by about 25%. But since for other algorithms we used statistics for unimproved versions, the data below are also taken from runs of the basic version.

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by 2-opt and 2.5-opt algorithms. Random planar problems, N=100.
Problem # Method Time Best Avg Worst
1NN100100.00108.79119.13
NN + 2opt136082.4186.2891.53
NN + 2.5opt197382.2486.1090.26
2NN100100.00112.19126.04
NN + 2opt107590.7495.0798.61
NN + 2.5opt156290.6794.3798.81
3NN100100.00106.74115.58
NN + 2opt125389.0892.9698.44
NN + 2.5opt281388.8792.0394.75
4NN100100.00106.17113.25
NN + 2opt198087.0291.5095.35
NN + 2.5opt177387.0090.5993.88
5NN100100.00106.40116.05
NN + 2opt104687.4992.3097.11
NN + 2.5opt155987.0391.5496.61
6NN100100.00108.92117.85
NN + 2opt136090.2895.28100.53
NN + 2.5opt156091.0194.90101.26
7NN100100.00105.35112.60
NN + 2opt97588.5592.3596.38
NN + 2.5opt234388.7791.7197.12
8NN100100.00114.47122.94
NN + 2opt104691.3394.8099.07
NN + 2.5opt187391.1094.69101.03
9NN100100.00106.94119.14
NN + 2opt146089.7093.7999.39
NN + 2.5opt187388.8192.8097.89
10NN100100.00108.41117.03
NN + 2opt146088.8693.7497.66
NN + 2.5opt166688.6492.8398.34

Tuesday, March 14, 2017

Or-opt

Or-opt idea

Or-opt is an algorithm proposed by Ilhan Or in 1976: a segment consisting of 3, 2 or 1 consecutive cities is shifted to other position in tour. This means that Or-opt is a restricted version of 3-opt and therefore it should be somewhat faster, but produce worse tours than full 3-optimization.

Or-opt move

Each move used by Or-opt is just Segment Shift:

proc Make_Segment_Shift_Move(tour: var Tour_Array;
                             i, j, k: Tour_Index) =
  ## Performs given Segment Shift move
  Shift_Segment(tour, i, j, k)

Gain from Or-opt move

proc Gain_From_Segment_Shift(X1, X2, Y1, Y2, Z1, Z2: City_Number):
                             Length_Gain =
  ## Gain of tour length which can be obtained by performing Segment
  ## Shift
  # Cities from X2 to Y1 would be moved from its current position,
  # between X1 and Y2, to position between cities Z1 and Z2.
  # Assumes: X1!=Z1
  #          X2==successor(X1); Y2==successor(Y1); Z2==successor(Z1)
  let del_Length = distance(X1, X2) + distance(Y1, Y2) + distance(Z1, Z2)
  let add_Length = distance(X1, Y2) + distance(Z1, X2) + distance(Y1, Z2)
  result = del_Length - add_Length

Or-opt using the first improving move

proc LS_Or_opt_Take_First(tour: var Tour_Array) =
  ## Optimizes the given tour using Or-opt
  # Shortens the tour by repeating Segment Shift moves for segment
  # length equal 3, 2, 1 until no improvement can by done: in every
  # iteration immediatelly makes permanent the first move found that
  # gives any length gain.
  var
    locallyOptimal: bool = false
    i, j, k: Tour_Index
    X1, X2, Y1, Y2, Z1, Z2: City_Number

  while not locallyOptimal:
    locallyOptimal = true

    for segmentLen in countdown(3, 1):
    
      block two_loops:
        for pos in 0 .. N-1:
          i = pos
          X1 = tour[i]
          X2 = tour[(i + 1) mod N]
          
          j = (i + segmentLen) mod N
          Y1 = tour[j]
          Y2 = tour[(j + 1) mod N]

          for shift in segmentLen+1 .. N-1:
            k  = (i + shift) mod N
            Z1 = tour[k]
            Z2 = tour[(k + 1) mod N]

            if Gain_From_Segment_Shift(X1, X2, Y1, Y2, Z1, Z2) > 0:
              Make_Segment_Shift_Move(tour, i, j, k)
              locallyOptimal = false
              break two_loops

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by 2-opt and Or-opt algorithms. Random planar problems, N=100.
Problem # Method Time Best Avg Worst
1NN100100.00111.23122.80
NN + 2opt146093.9696.65100.64
NN + Or-opt646094.6497.26102.54
NN + 2opt + Or-opt312692.2994.5997.28
2NN100100.00110.64119.25
NN + 2opt126886.9390.1995.17
NN + Or-opt810689.0093.09101.05
NN + 2opt + Or-opt321885.3088.3391.93
3NN100100.00109.22115.52
NN + 2opt107585.4889.2093.58
NN + Or-opt654386.9192.3598.14
NN + 2opt + Or-opt312585.0187.1990.52
4NN100100.00111.96119.75
NN + 2opt136087.3791.0296.94
NN + Or-opt707987.9295.32104.80
NN + 2opt + Or-opt291985.8389.1994.70
5NN100100.00109.49117.46
NN + 2opt107593.4396.66103.79
NN + Or-opt468793.1297.78103.31
NN + 2opt + Or-opt283192.2894.29100.36
6NN100100.00105.35112.56
NN + 2opt125387.0290.9798.85
NN + Or-opt614687.4392.79100.48
NN + 2opt + Or-opt291386.8489.1894.75
7NN100100.00107.32117.78
NN + 2opt125390.0694.1497.89
NN + Or-opt635392.1895.90100.65
NN + 2opt + Or-opt343988.4491.2895.44
8NN100100.00115.11124.44
NN + 2opt166688.1091.0897.24
NN + Or-opt687990.6596.09102.98
NN + 2opt + Or-opt374686.1688.8492.96
9NN100100.00108.46113.81
NN + 2opt114692.3597.09101.10
NN + Or-opt542094.4099.33105.41
NN + 2opt + Or-opt291391.6794.7798.97
10NN100100.00108.04117.50
NN + 2opt114688.2093.3098.54
NN + Or-opt583388.5895.29102.31
NN + 2opt + Or-opt281388.1291.3494.55

Monday, March 13, 2017

Shifting a segment

Shifting a segment is a generalization of the concept of shifting a node: instead of shifting a single city, we shift a whole segment (chain of cities) to other position in tour.

The idea:

A B
A B

Diagrams:

X1 X2 Y1 Y2 Z1 Z2
X1 X2 Y1 Y2 Z1 Z2

Segment Shift move is a the same as the symmetric pure 3-opt move (case 7).

Node Shift is a special case of Segment Shift, when a segment has length of 1 (j==i+1).

proc Shift_Segment(tour: var Tour_Array; i, j, k: Tour_Index) =
  ## Shifts the segment of tour:
  # cities from t[i+1] to t[j] from their current position to position
  # after current city t[k], that is between cities t[k] and t[k+1].
  # Assumes:  k, k+1 are not within the segment [i+1..j]
  let
    segmentSize = (j - i + N) mod N
    shiftSize = ((k - i + N)  - segmentSize + N) mod N
    offset = i + 1 + shiftSize
  var
    pos: Tour_Index
    segment: seq[City_Number] = newSeq[City_Number](segmentSize)

  # make a copy of the segment before shift
  for counter in 0 .. segmentSize-1:
    segment[pos] = tour[(pos + i+1) mod N]

  # shift to the left by segmentSize all cities between old position
  # of right end of the segment and new position of its left end
  pos  = (i + 1) mod N
  for counter in 1 .. shiftSize:
    tour[pos] = tour[(pos + segmentSize) mod N]
    pos  = (pos + 1) mod N

  # put the copy of the segment into its new place in the tour
  for pos in 0 .. segmentSize-1:
    tour[(pos + offset) mod N] = segment[pos]

Saturday, March 11, 2017

Node Shift

Node Shift move

Node Shift is a simple tour optimization heuristics: it consists in shifting a city to other position in tour.

The idea:

A
A

Diagrams:

X0 X0_pred X0_succ Y1 Y2
X0 X0_pred X0_succ Y1 Y2

Node Shift move is obtained by removing three links and adding three new links, and therefore it is a special case of 3-opt move:

Gain from Node Shift move

proc Gain_From_Node_Shift(X0_pred, X0, X0_succ,
                          Y1, Y2: City_Number): Length_Gain =
  ## Gain of tour length which can be obtained by performing Node Shift
  # City X0 would be moved from its current position, between X0_pred
  # and X0_succ, to position between cities Y1 and Y2.
  # Assumes: X0_pred!=Y1 (the same as: X0!=Y2);
  #   X0==successor(X0_pred); X0_succ==successor(X0); Y2==successor(Y1)
  let del_Length = distance(X0_pred, X0) + distance(X0, X0_succ) +
                   distance(Y1, Y2)
  let add_Length = distance(X0_pred, X0_succ) +
                   distance(Y1, X0) + distance(X0, Y2)
  result = del_Length - add_Length

Node Shift move implementation

proc Make_Node_Shift_Move(tour: var Tour_Array;
                          i, j: Tour_Index) =
  ## Performs given Node Shift move
  # Shifts the city t[i] from its current position to position
  # after current city t[j], i.e. between cities t[j] and t[j+1].
  # This is a special, simple case of segment shift, thus a special
  # case of 3-opt move.
  var
    left, right: Tour_Index
    X0: City_Number

  let shiftSize = (j - i + 1 + N) mod N

  X0 = tour[i]
  left = i
  for counter in 1 .. shiftSize:
    right = (left + 1) mod N
    tour[left] = tour[right]
    left = right
  # end_for counter loop
  tour[j] = X0

Iterative Node Shift optimization

proc LS_Node_Shift(tour: var Tour_Array) =
  ## Iteratively optimizes given tour using Node_Shift moves.
  # Shortens tour by repeating Node Shift moves until no improvement
  # can by done: in every iteration immediatelly makes permanent
  # the first move found that gives any length gain.
  var
    locallyOptimal: bool = false
    i, j: Tour_Index
    X0_pred, X0, X0_succ, Y1, Y2: City_Number

  while not locallyOptimal:
    locallyOptimal = true

    for counter_1 in 0 .. N-1:
      i = counter_1  # current position of city to shift
      X0_pred = tour[(i + N - 1) mod N]
      X0      = tour[i]
      X0_succ = tour[(i + 1) mod N]

      for counter_2 in 1 .. N-2:
        j  = (i + counter_2) mod N  # new position of the city
        Y1 = tour[j]
        Y2 = tour[(j + 1) mod N]

        if Gain_From_Node_Shift(X0_pred, X0, X0_succ, Y1, Y2) > 0:
          Make_Node_Shift_Move(tour, i, j)
          locallyOptimal = false
          break

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by 2-opt and 2-opt plus Node Shift algorithms. Random planar problems, N=100.
Problem # Method Time Best Avg Worst
1NN100100.00106.01114.87
NN + 2opt125384.3889.8396.82
NN + 2opt + NodeShift145983.9487.9091.71
2NN100100.00107.45116.93
NN + 2opt156288.0592.4096.60
NN + 2opt + NodeShift156287.1690.3494.04
3NN100100.00111.21121.46
NN + 2opt156687.4191.5295.76
NN + 2opt + NodeShift166684.9988.7992.92
4NN100100.00106.14116.83
NN + 2opt114690.5293.4296.88
NN + 2opt + NodeShift135390.1891.9095.45
5NN100100.00106.30111.75
NN + 2opt83389.6393.3598.55
NN + 2opt + NodeShift125388.5891.6895.80
6NN100100.00112.73120.18
NN + 2opt125389.4693.6098.49
NN + 2opt + NodeShift145988.2792.0096.30
7NN100100.00109.50124.34
NN + 2opt146088.6691.0596.24
NN + 2opt + NodeShift156686.4789.2792.40
8NN100100.00108.63117.53
NN + 2opt114692.4796.80102.24
NN + 2opt + NodeShift156691.8695.44101.14
9NN100100.00105.75112.56
NN + 2opt40088.7091.6896.63
NN + 2opt + NodeShift46587.8590.0394.57
10NN100100.00115.43123.97
NN + 2opt146091.7795.11100.35
NN + 2opt + NodeShift145990.8893.8798.72

Node Swap

Node Swap move

Node Swap is a simple tour optimization heuristics: it consists in swapping two cities in tour.

The idea:

A B
B A

Diagrams:

X0 Y0 X0_pred X0_succ Y0_pred Y0_succ
X0 Y0 X0_pred X0_succ Y0_pred Y0_succ

Node Swap move is a special case of two subsequent 2-opt moves: the first including both cities and the second without them. It involves removing four links and adding four new links, therefore is a specific type of 4-opt move.

Gain from Node Swap move

proc Gain_From_Node_Swap(X0_pred, X0, X0_succ,
                         Y0_pred, Y0, Y0_succ: City_Number): Length_Gain =
  ## Gain of tour length which can be obtained by performing Node Swap
  ## City X0 would be moved from its current position, between X0_pred
  ## and X0_succ, to the position of city Y0 and vice versa
  # Assumes: X0!=Y0
  #          X0=successor(X0_pred); X0_succ==successor(X0);
  #          Y0=successor(Y0_pred); Y0_succ==successor(Y0);
  var
    del_Length, add_Length: Length_Gain
  if (Y0==X0_succ) or (X0_succ==Y0_pred):
    del_Length = distance(X0_pred, X0) + distance(Y0, Y0_succ)
    add_Length = distance(X0_pred, Y0) + distance(X0, Y0_succ)
  elif (X0==Y0_succ) or (Y0_succ==X0_pred):
    del_Length = distance(Y0_pred, Y0) + distance(X0, X0_succ)
    add_Length = distance(Y0_pred, X0) + distance(Y0, X0_succ)
  else:
    del_Length = distance(X0_pred, X0) + distance(X0, X0_succ) +
                 distance(Y0_pred, Y0) + distance(Y0, Y0_succ)
    add_Length = distance(X0_pred, Y0) + distance(Y0, X0_succ) +
                 distance(Y0_pred, X0) + distance(X0, Y0_succ)
  #end_if
  result = del_Length - add_Length

Node Swap move implementation

proc Make_Node_Swap_Move(tour: var Tour_Array;
                         i, j: Tour_Index) =
  ## Performs given Node Swap move: swaps city t[i] with city city t[j]
  # This is a special, simple case of two subsequent 2-opt moves:
  # the first including both cities and the second without them
  #   -a-t[i]-b-t[j]-a-
  #   -a-(t[i]-b-t[j])'-a-  == -a-t[j]-b'-t[i]-a-
  #   -a-t[j]-(b')'-t[i]-a- == -a-t[j]-b-t[i]-a-
  swap(tour[i], tour[j])

Iterative Node Swap optimization

proc LS_Node_Swap(tour: var Tour_Array) =
  ## Optimizes the given tour using Node_Swap.
  # Shortens tour by repeating Node Swap moves until no improvement
  # can by done: in every iteration immediatelly makes permanent
  # the first move found that gives any length gain.
  var
    locallyOptimal: bool = false
    i, j: Tour_Index
    X0_pred, X0, X0_succ, Y0_pred, Y0, Y0_succ: City_Number

  while not locallyOptimal:
    locallyOptimal = true

    for counter_1 in 0 .. N-1:
      i = counter_1  # position of the first city
      X0_pred = tour[(N + i - 1) mod N]
      X0      = tour[i]
      X0_succ = tour[(i + 1) mod N]

      for counter_2 in counter_1+1 .. N-1:
        j  = counter_2  # position of the second city
        Y0_pred = tour[(N + j - 1) mod N]
        Y0 = tour[j]
        Y0_succ = tour[(j + 1) mod N]

        if Gain_From_Node_Swap(X0_pred, X0, X0_succ,
                               Y0_pred, Y0, Y0_succ) > 0:
          Make_Node_Swap_Move(tour, i, j)
          locallyOptimal = false
          break

Some simple statistics

Relative time, best length, average length and worst length for all tours created by NN heuristic and improved by 2-opt and Node Swap algorithms. Random planar problems, N=100.
Problem # Method Time Best Avg Worst
1NN100100.00107.17119.19
NN + 2opt97589.5392.7799.31
NN + NodeSwap29394.61103.80116.17
2NN100100.00107.49118.89
NN + 2opt114690.7894.5198.18
NN + NodeSwap41995.23101.53110.37
3NN100100.00107.68114.62
NN + 2opt117589.9793.0897.98
NN + NodeSwap29394.55102.50111.75
4NN100100.00109.62122.88
NN + 2opt107590.7493.5598.18
NN + NodeSwap29397.54106.47118.40
5NN100100.00106.67113.50
NN + 2opt117590.9093.9197.43
NN + NodeSwap29394.35101.54107.83
6NN100100.00109.94120.93
NN + 2opt125387.1490.8497.60
NN + NodeSwap41993.50101.63110.62
7NN100100.00109.60121.26
NN + 2opt146086.2290.3294.58
NN + NodeSwap42095.34103.06114.27
8NN100100.00110.14121.53
NN + 2opt97586.0590.9298.23
NN + NodeSwap48791.82105.68115.13
9NN100100.00110.14121.53
NN + 2opt107586.0590.9298.23
NN + NodeSwap39391.82105.68115.13
10NN100100.00111.65120.86
NN + 2opt83395.0997.70102.49
NN + NodeSwap41997.68106.43117.24