I have some trouble with my c++ implementation of a solver for TSP. I basically have two problems.
The first one regards the comparison in the twoOptSearch() between min_change and change. I need to use that small offset in order to avoid an infinite loop. I don't know why I have to use it and I don't know why there is no need in the corresponding comparison in threeOptSearch.
The second one regards the threeOptSearch. From the theory, I know that a tour that is 3-opt optimal is also 2-opt optimal. But if I call twoOptSearch after threeOptSearch, it is able to reduce the length of the tour.
Please, if you are able to help me with these troubles or simply you have some suggestion to give me, feel free to answer. I would appreciate it a lot.
These are the files that are involved.
tsp.h
#ifndef _TSP
#define _TSP
#include <vector> // for vector
#include <math.h> // for sqrt.
class City{
public:
int id;
double x_coord;
double y_coord;
};
class Tour{
int size; // number of cities in tour, that is the total number of different cities plus one (the tour ends with the same city from which it started)
std::vector< std::vector<double> > distances; // distance between cities
double length; // length of the tour (it's updated with each tour change)
void twoOptMove(int, int);
void threeOptMove(char, int, int, int);
public:
std::vector<City> cities;
Tour(int);
int getSize(){return size;}; // return size
void computeDistances(void); // compute distances matrix, that has dimensions (size-1)x(size-1)
double getDistance(int i, int j){return distances[i-1][j-1];}; // given two ids, return distance between corresponding cities
void computeLength(void); // compute length
double getLength(void){return length;}; // return length
void nearestNeighbourSearch(void); // do the nearest neighbour search
void twoOptSearch(void); // do the 2-opt search
void threeOptSearch(void); // do the 3-opt search
~Tour();
};
#endif
tsp.cpp
#include "tsp.h"
Tour::Tour(int tour_size) : size(tour_size){
cities.resize(size);
distances.resize(size - 1);
for (int i = 0; i < size - 1; ++i)
distances[i].resize(size - 1);
}
void Tour::computeDistances(){
int h_x, h_y;
for(int i = 0; i < size - 1; i++){
for(int j = i; j < size - 1; j++){
h_x = cities[i].x_coord - cities[j].x_coord;
h_y = cities[i].y_coord - cities[j].y_coord;
h_x *= h_x;
h_y *= h_y;
distances[j][i] = distances[i][j] = sqrt(h_x + h_y);
}
}
}
void Tour::computeLength(){
for(int i = length = 0; i < size - 1; i++)
length += getDistance(cities[i].id, cities[i+1].id);
}
void Tour::nearestNeighbourSearch(){
// Do nearest neighbor search
int curr = 0;
int succ;
int flag;
std::vector<City> help_cities;
help_cities.resize(size);
// set start city
succ = curr;
help_cities[0] = cities[curr];
for(int i = 1; i < size - 1; i++){
for(int j = 0; j < size - 1; j++){
flag = true;
// verify that it is not already present in tour
for(int k = 0; flag && k < i; k++)
if(cities[j].id == help_cities[k].id) flag = false;
// choose the first unvisited city if a candidate city has not been selected yet
if(flag && curr == succ)
succ = j;
// if a candidate city has been selected, verify which one is the best
if(flag && getDistance(cities[curr].id, cities[succ].id) > getDistance(cities[curr].id, cities[j].id))
succ = j;
}
// chose the best city as the current one and save its position
curr = succ;
help_cities[i] = cities[curr];
}
// set end city equal to the start one
help_cities[size-1] = help_cities[0];
// update tour
for(int i = 0; i < size; i++){
cities[i] = help_cities[i];
}
computeLength();
}
void Tour::twoOptSearch(){
/* In each iteration, apply a 2-opt move. That is, find pair of edges
* (a,a+1) and (b,b+1) such that replacing them with (a,b) and (a+1,b+1)
* minimizes tour length.
*
* a b a b
* ────→● ●←──── ────→●───●────→
* ╲ ╱
* ╳ ==>
* ╱ ╲
* ←────● ●────→ ←────●───●←────
* b+1 a+1 b+1 a+1
*/
double min_change, change;
int a_min, b_min;
double od1, od2;
double nd1, nd2;
// repeat until there are not new improvement
while(min_change){
min_change = 0;
for(int a = 0; a < size - 2; a++){
for(int b = a + 1; b < size - 1; b++){
od1 = getDistance(cities[a].id, cities[a+1].id);
od2 = getDistance(cities[b].id, cities[b+1].id);
nd1 = getDistance(cities[a].id, cities[b].id);
nd2 = getDistance(cities[a+1].id, cities[b+1].id);
change = nd1 + nd2 - od1 - od2;
// ERROR: without this small offset there is an infinite loop.
if(min_change - 0.0000000000001 > change){
min_change = change;
a_min = a;
b_min = b;
}
}
}
if(min_change) twoOptMove(a_min, b_min);
}
// update length
computeLength();
}
void Tour::twoOptMove(int a, int b){
/* - take tour[0] to tour[a] and add them in order to new tour
* - take tour[a+1] to tour[b] and add them in reverse order to new tour
* - take tour[b+1] to end and add them in order to new tour
*
* a b a b
* ────→● ●←──── ────→●───●────→
* ╲ ╱
* ╳ ==>
* ╱ ╲
* ←────● ●────→ ←────●───●←────
* b+1 a+1 b+1 a+1
*/
std::vector<City> help_cities;
help_cities.resize(size);
// compute new tour
for(int k = 0; k <= a; k++){
help_cities[k] = cities[k];
}
for(int k = 0; k < b - a; k++){
help_cities[a+k+1] = cities[b-k];
}
for(int k = b + 1; k < size; k++){
help_cities[k] = cities[k];
}
// update tour
for(int i = 0; i < size; i++){
cities[i] = help_cities[i];
}
}
void Tour::threeOptSearch(){
/* In each iteration, apply a 3-opt move. That is, find three edges
* (a,a+1), (b,b+1) and (c,c+1) such that replacing them with three
* other edges minimizes tour length. Each 3-opt move can be performed
* with one, two or three 2-opt valid move.
* REMARK: a move is called valid if it reduces the tour length.
*/
double min_change, change;
double change_a, change_b, change_c, change_d, change_e, change_f, change_g;
char move_type_min, move_type;
int a_min, b_min, c_min;
double od1, od2, od3;
double nd1, nd2, nd3;
// repeat until there are not new improvement
while(min_change){
min_change = 0;
for(int a = 0; a < size - 3; a++){
for(int b = a + 1; b < size - 2; b++){
for(int c = b + 1; c < size - 1; c++){
od1 = getDistance(cities[a].id, cities[a+1].id);
od2 = getDistance(cities[b].id, cities[b+1].id);
od2 = getDistance(cities[c].id, cities[c+1].id);
/* a)
* ↑ ↓ ↑ ↓
* b+1 ● ● c b+1 ● ● c
* ╲ ╱ ╱ ╲
* a ╳ a+1 3 * 2-opt moves ╱ ╲
* ──→●───╱─╲───●──→ ====> ──→● a a+1 ●──→
* ╱ ╲
* ←────● ●←──── ←────●─────●←────
* c+1 b c+1 b
*/
nd1 = getDistance(cities[a].id, cities[b+1].id);
nd2 = getDistance(cities[c].id, cities[a+1].id);
nd3 = getDistance(cities[b].id, cities[c+1].id);
change_a = nd1 + nd2 + nd3 - od1 - od2 - od3;
/* b)
* ↑ ↓ ↑ ↓
* b+1 ● ● c b+1 ● ● c
* │ │ ╱ ╲
* a └───┼──┐ 2 * 2-opt moves ╱ ╲
* ──→●┐ ┌──┘ ●←── ====> ──→● a b ●──→
* └──┼───┐ b
* ←────●─┘ ●────→ ←────●─────●←────
* c+1 a+1 c+1 a+1
*/
nd1 = getDistance(cities[a].id, cities[b+1].id);
nd2 = getDistance(cities[c].id, cities[b].id);
nd3 = getDistance(cities[a+1].id, cities[c+1].id);
change_b = nd1 + nd2 + nd3 - od1 - od2 - od3;
/* c)
* ↓ ↑ ↑ ↓
* b ● ● a+1 b ● ● a+1
* │ │ ╱ ╲
* ┌──┼───┘ c 2 * 2-opt moves ╱ ╲
* ──→● └──┐ ┌●←── ====> ──→● a c ●──→
* a ┌───┼──┘
* ←────● └─●────→ ←────●─────●←────
* c+1 b+1 c+1 b+1
*/
nd1 = getDistance(cities[a].id, cities[b].id);
nd2 = getDistance(cities[a+1].id, cities[c].id);
nd3 = getDistance(cities[b+1].id, cities[c+1].id);
change_c = nd1 + nd2 + nd3 - od1 - od2 - od3;
/* d)
* ↓ ↑ ↑ ↓
* c ● ● b+1 c ● ● b+1
* │ │ ╱ ╲
* a │ │ a+1 2 * 2-opt moves ╱ ╲
* ──→●──┼───┼──●──→ ====> ──→● a a+1 ●──→
* │ │
* ←────●┘ └●←──── ←────●─────●←────
* c+1 b c+1 b
*/
nd1 = getDistance(cities[a].id, cities[c].id);
nd2 = getDistance(cities[b+1].id, cities[a+1].id);
nd3 = getDistance(cities[b].id, cities[c+1].id);
change_d = nd1 + nd2 + nd3 - od1 - od2 - od3;
/* e)
* ↓ ↑ ↑ ↓
* c ● ● b+1 c ● ● b+1
* │ ╲ ╱ ╲
* a │ ╲ 1 * 2-opt move ╱ ╲
* ──→●──┼────┐ ●←── ====> ──→● a b ●──→
* │ │ b
* ←────●┘ ●────→ ←────●─────●←────
* c+1 a+1 c+1 a+1
*/
nd1 = getDistance(cities[a].id, cities[c].id);
nd2 = getDistance(cities[b+1].id, cities[b].id);
nd3 = getDistance(cities[a+1].id, cities[c+1].id);
change_e = nd1 + nd2 + nd3 - od1 - od2 - od3;
/* f)
* ↑ ↓ ↑ ↓
* a+1 ● ● b a+1 ● ● b
* ╱ │ ╱ ╲
* ╱ │ 1 * 2-opt move ╱ ╲
* ──→● ┌───┼──●←── ====> ──→● a c ●──→
* a │ │ c
* ←────●┘ └●────→ ←────●─────●←────
* c+1 b+1 c+1 b+1
*/
nd1 = getDistance(cities[a].id, cities[a+1].id);
nd2 = getDistance(cities[b].id, cities[c].id);
nd3 = getDistance(cities[b+1].id, cities[c+1].id);
change_f = nd1 + nd2 + nd3 - od1 - od2 - od3;
/* g)
* ↓ ↑ ↑ ↓
* b ● ● a+1 b ● ● a+1
* │ │ ╱ ╲
* └───┼──┐ 1 * 2-opt move ╱ ╲
* ──→●──────┘ ●──→ ====> ──→● a b+1 ●──→
* a b+1
* ←────●─────●←──── ←────●─────●←────
* c+1 c c+1 c
*/
nd1 = getDistance(cities[a].id, cities[b].id);
nd2 = getDistance(cities[a+1].id, cities[b+1].id);
nd3 = getDistance(cities[c].id, cities[c+1].id);
change_g = nd1 + nd2 + nd3 - od1 - od2 - od3;
// select the best change for the current a, b, c.
change = 0;
if(change > change_a){
change = change_a;
move_type = 'a';
}
if(change > change_b){
change = change_b;
move_type = 'b';
}
if(change > change_c){
change = change_c;
move_type = 'c';
}
if(change > change_d){
change = change_d;
move_type = 'd';
}
if(change > change_e){
change = change_e;
move_type = 'e';
}
if(change > change_f){
change = change_f;
move_type = 'f';
}
if(change > change_g){
change = change_g;
move_type = 'g';
}
// WARNING: there is no need for the small offset in this case. I don't know why.
if(min_change > change){
min_change = change;
a_min = a;
b_min = b;
c_min = c;
move_type_min = move_type;
}
}
}
}
if(min_change) threeOptMove(move_type_min, a_min, b_min, c_min);
}
// update length
computeLength();
}
void Tour::threeOptMove(char move_type, int a, int b, int c){
std::vector<City> help_cities;
help_cities.resize(size);
switch(move_type){
case 'a':
/* a)
* ↑ ↓ ↑ ↓
* b+1 ● ● c b+1 ● ● c
* ╲ ╱ ╱ ╲
* a ╳ a+1 3 * 2-opt moves ╱ ╲
* ──→●───╱─╲───●──→ ====> ──→● a a+1 ●──→
* ╱ ╲
* ←────● ●←──── ←────●─────●←────
* c+1 b c+1 b
*
* - take tour[0] to tour[a] and add them in order to new tour
* - take tour[b+1] to tour[c] and add them in order to new tour
* - take tour[a+1] to tour[b] and add them in order to new tour
* - take tour[c+1] to end and add them in order to new tour
*/
// compute new_tour
for(int k = 0; k <= a; k++){
help_cities[k] = cities[k];
}
for(int k = 0; k < c - b; k++){
help_cities[a+k+1] = cities[b+k+1];
}
for(int k = 0; k < b - a; k++){
help_cities[a+c-b+k+1] = cities[a+k+1];
}
for(int k = c + 1; k < size; k++){
help_cities[k] = cities[k];
}
break;
case 'b':
/* b)
* ↑ ↓ ↑ ↓
* b+1 ● ● c b+1 ● ● c
* │ │ ╱ ╲
* a └───┼──┐ 2 * 2-opt moves ╱ ╲
* ──→●┐ ┌──┘ ●←── ====> ──→● a b ●──→
* └──┼───┐ b
* ←────●─┘ ●────→ ←────●─────●←────
* c+1 a+1 c+1 a+1
*
* - take tour[0] to tour[a] and add them in order to new tour
* - take tour[b+1] to tour[c] and add them in order to new tour
* - take tour[a+1] to tour[b] and add them in reverse order to new tour
* - take tour[c+1] to end and add them in order to new tour
*/
// compute new_tour
for(int k = 0; k <= a; k++){
help_cities[k] = cities[k];
}
for(int k = 0; k < c - b; k++){
help_cities[a+k+1] = cities[b+k+1];
}
for(int k = 0; k < b - a; k++){
help_cities[a+c-b+k+1] = cities[b-k];
}
for(int k = c + 1; k < size; k++){
help_cities[k] = cities[k];
}
break;
case 'c':
/* c)
* ↓ ↑ ↑ ↓
* b ● ● a+1 b ● ● a+1
* │ │ ╱ ╲
* ┌──┼───┘ c 2 * 2-opt moves ╱ ╲
* ──→● └──┐ ┌●←── ====> ──→● a c ●──→
* a ┌───┼──┘
* ←────● └─●────→ ←────●─────●←────
* c+1 b+1 c+1 b+1
*
* - take tour[0] to tour[a] and add them in order to new tour
* - take tour[a+1] to tour[b] and add them in reverse order to new tour
* - take tour[b+1] to tour[c] and add them in reverse order to new tour
* - take tour[c+1] to end and add them in order to new tour
*/
// compute new_tour
for(int k = 0; k <= a; k++){
help_cities[k] = cities[k];
}
for(int k = 0; k < b - a; k++){
help_cities[a+k+1] = cities[b-k];
}
for(int k = 0; k < c - b; k++){
help_cities[b+k+1] = cities[c-k];
}
for(int k = c + 1; k < size; k++){
help_cities[k] = cities[k];
}
break;
case 'd':
/* d)
* ↓ ↑ ↑ ↓
* c ● ● b+1 c ● ● b+1
* │ │ ╱ ╲
* a │ │ a+1 2 * 2-opt moves ╱ ╲
* ──→●──┼───┼──●──→ ====> ──→● a a+1 ●──→
* │ │
* ←────●┘ └●←──── ←────●─────●←────
* c+1 b c+1 b
*
* - take tour[0] to tour[a] and add them in order to new tour
* - take tour[b+1] to tour[c] and add them in reverse order to new tour
* - take tour[a+1] to tour[b] and add them in order to new tour
* - take tour[c+1] to end and add them in order to new tour
*/
// compute new_tour
for(int k = 0; k <= a; k++){
help_cities[k] = cities[k];
}
for(int k = 0; k < c - b; k++){
help_cities[a+k+1] = cities[c-k];
}
for(int k = 0; k < b - a; k++){
help_cities[a+c-b+k+1] = cities[a+k+1];
}
for(int k = c + 1; k < size; k++){
help_cities[k] = cities[k];
}
break;
case 'e':
/* e)
* ↓ ↑ ↑ ↓
* c ● ● b+1 c ● ● b+1
* │ ╲ ╱ ╲
* a │ ╲ 1 * 2-opt move ╱ ╲
* ──→●──┼────┐ ●←── ====> ──→● a b ●──→
* │ │ b
* ←────●┘ ●────→ ←────●─────●←────
* c+1 a+1 c+1 a+1
*/
// compute new_tour
for(int k = 0; k <= a; k++){
help_cities[k] = cities[k];
}
for(int k = 0; k < c - a; k++){
help_cities[a+k+1] = cities[c-k];
}
for(int k = c + 1; k < size; k++){
help_cities[k] = cities[k];
}
break;
case 'f':
/* f)
* ↑ ↓ ↑ ↓
* a+1 ● ● b a+1 ● ● b
* ╱ │ ╱ ╲
* ╱ │ 1 * 2-opt move ╱ ╲
* ──→● ┌───┼──●←── ====> ──→● a c ●──→
* a │ │ c
* ←────●┘ └●────→ ←────●─────●←────
* c+1 b+1 c+1 b+1
*/
// compute new_tour
for(int k = 0; k <= b; k++){
help_cities[k] = cities[k];
}
for(int k = 0; k < c - b; k++){
help_cities[b+k+1] = cities[c-k];
}
for(int k = c + 1; k < size; k++){
help_cities[k] = cities[k];
}
break;
case 'g':
/* g)
* ↓ ↑ ↑ ↓
* b ● ● a+1 b ● ● a+1
* │ │ ╱ ╲
* └───┼──┐ 1 * 2-opt move ╱ ╲
* ──→●──────┘ ●──→ ====> ──→● a b+1 ●──→
* a b+1
* ←────●─────●←──── ←────●─────●←────
* c+1 c c+1 c
*/
// compute new_tour
for(int k = 0; k <= a; k++){
help_cities[k] = cities[k];
}
for(int k = 0; k < b - a; k++){
help_cities[a+k+1] = cities[b-k];
}
for(int k = b + 1; k < size; k++){
help_cities[k] = cities[k];
}
break;
}
// update tour
for(int i = 0; i < size; i++){
cities[i] = help_cities[i];
}
}
Tour::~Tour(){
// do nothing
}
P.S. the tour is filled with cities with incremental index starting from 1.
P.P.S. I tried to use 3-opt and 2-opt as tags, but I don't have enough reputation
