Travelling Salesman in Qt

Question

I am writing a recursive function based on the Travelling Salesman Problem.

Is this a correct way of doing things, or am I creating memory leaks?

Prerequisites are:

• a two dimensional matrix distance_matrix[][] to store the distances between the cities
• a global int called current_minimum_cost which will be MAX_INT for starters (or any high 'impossible' number), it keeps track of the cost of the current 'optimal' path
• remaining_cities: a QList of ints which at starting point contains all the cities to visit (except the HOME_CITY) (every city is represented by an integer, the same integer is used to look up the distance between two cities by using distance_matrix[cityA][cityB].
• path: a QList of ints which stores the current path followed (so 0 , 7, 2, 5) means the salesman goes from his HOME_CITY to city #7 to city #2 to city #5 . At starting point it only contains the HOME_CITY (which is 0)
• after doing his tour, the salesman has to return to his HOME_CITY (0)

void tsp(QList<int>* path, int total_cost, QList<int>* remaining_cities)
{
    int last_visited_city = path->last();

    if(total_cost > current_minimum_cost)
    {
        //cost higher than current_minimum_cost => ABORT
        return;
    }

    if(remaining_cities->empty())
    {
       //all cities visited

       //only add the distance back to the home city (salesman returns home)
       //path += HOME_CITY en total_cost += terug naar HOME_CITY";
       path->push_back(HOME_CITY);

       total_cost += distance_matrix[last_visited_city][HOME_CITY];

       if (total_cost <= current_minimum_cost)
       {
           // we have a new winner!
           current_minimum_cost = total_cost;

           // forget about the previous 'winner' and save the visited path
           tsp_solution->clear();
           foreach(int element , *path)
               tsp_solution->push_back(element);
       }
       // if not, than it means that "returning home" makes the trip more costly than the previous 'winner' (we always 'return home', so this is also taken into account with the previous 'winner')
       return;
    }

    // calculate the new possibilities to travel to from the current city:    
    foreach(int next_city_to_visit , *remaining_cities)
    {
      QList<int>*newPATH = new QList<int>();
      newPATH->append(*path);
      newPATH->push_back(next_city_to_visit);

      int total_cost_for_this_new_path = total_cost + distance_matrix[last_visited_city][next_city_to_visit];

      // make a list of 'remaining cities' to visti:
      QList<int>* newREMAINING_CITIES = new QList<int>();
      newREMAINING_CITIES->append(*remaining_cities);

      // of course we do not have to visit 'next_city_to_visit' twice
      newREMAINING_CITIES->removeOne(next_city_to_visit);

      // recursion:
      tsp ( newPATH , total_cost_for_this_new_path , newREMAINING_CITIES );
      delete newPATH;
      delete newREMAINING_CITIES;
    }
}

My questions:

1. Does this function make sense? It seems to work for all my test cases.
2. Once I go up to 20 cities, it takes a long while to calculate, over an hour.
3. On memory management: is issuing delete newPATH enough to clear the memory? Or should I issue a delete for all members of the list newPATH? (I guess not as the members are all ints)

Answer 1

1. As a brute force solution to the Traveling Salesman it seems to make sense. The algorithm seems to test every single combination of cities.

   The method used for listing remaining cities seems to be a simple one. You could avoid a lot of allocations and adding/removing to lists by using another method to track visited cities.

   For example, instead of creating a new list of remaining cities before the next recursion, just mark the city as visited using another table (make a city be a Class containing it's id in the table, and a bool visited).

   Then before the recursive tsp call you mark the current city as "visited", as the city will already have been visited for every successive recursive call under it. When that tsp call returns you simply mark the city as "not visited" again.

   Thus as the recursion goes deeper more and more cities are marked as "visited" until they all are, and you perform the last HOME_CITY check. And as they unwind (move back up) they all mark themselves as "not visited" again.

   This way you avoid not only a very expensive "new" and "delete" for EVERY recursion (see point 2.), but also the removeOne call, which actually has to iterate the whole list looking for the city to remove.

2. That a list of 20 cities takes over an hour might not be so strange. 20 cities might not sound like a lot, but this is a brute force algorithm. That means that it will check EVERY possible combination of traveling between the cities.

   How many combinations are that you ask?

   Well, for 19 cities(We don't count the HOME_CITY), we have to call a recursive tsp on all of them. Each of these 19 recursive tsp calls then have to call a recursive tsp call on 18 cities.

   Thus we get a Factorial: \$19! = 19*18*17...*3*2*1 = 121645100408832000\$.

   Yes, the tsp function is called 121645100408832000 times! That is 243290200817664000 calls to new and delete!

   121645100408832000 times filling the newREMAINING_CITIES with remaining_cities. 121645100408832000 times iterating the newREMAINING_CITIES list to remove the current city.

   Do you see how this, when running on only one cpu, can tend to use over an hour?

3. Your memory management looks ok. The rule is for every call to new, there should be a call to delete.

   You don't need to delete the entries in newPATH, as they have already deleted themself before you delete your newPATH.

Answer 2

I'd have serious doubts about its being the fastest possible code, but just for fun let's consider the simplest possible code to generate correct answers.

The standard library has next_permutation, which can generate all the possible permutations of a collection. So, let's consider how we do this. We start with a 2D array of distances between cities. We then have a collection to represent a route through the cities. Since each city must appear once in the collection, we can start with it just containing numbers from 0 to N in order. We then permute those, find the distance for each, keep the shortest, and we're (eventually) done.

#include <algorithm>
#include <vector>
#include <iostream>
#include <numeric>
#include <limits>

static const int max = 13;

void setup(int distances[max][max]) {
    for (int i = 0; i < max; i++) {
        distances[i][i] = 0;
        for (int j = i + 1; j < max; j++) {
            int d = rand() % 50;

            distances[i][j] = d;
            distances[j][i] = d;
        }
    }
}

int main() {
    int distances[max][max];

    setup(distances);

    std::vector<int> route(max);

    std::iota(route.begin(), route.end(), 0);

    int best_distance = std::numeric_limits<int>::max();
    do {
        int distance = 0;
        for (int i = 1; i < route.size(); i++) {
            distance += distances[route[i - 1]][route[i]];
            if (distance > best_distance)
                break;
        }
        distance += distances[route[max - 1]][0];

        if (distance < best_distance) best_distance = distance;
    } while (std::next_permutation(route.begin() + 1, route.end()));

    std::cout << "best distance: " << best_distance << "\n";
}

Now some might wonder why I'd question this being the fastest way to do the job. Generic algorithms are, after all, intended to provide nearly the speed of custom-written code--their runtime overhead is usually quite minimal.

While that's true, in this case an important optimization is missed. If we view the routes as a tree or a graph, we can stop traversing a branch as soon as we reach a node where the current distance is greater than the best distance so far. This can avoid visiting a large number of paths based on a single test.

Using next_permutation, however we don't get that optimization. Since we don't know the relationships between paths, information we obtain from traversing one path can't be used to eliminate other possible paths.

While probably pretty fast for small sets of data, this is likely to be substantially slower than many others when N gets large at all.

Note: this does not depend upon the triangle inequality (and the test data it generates almost certainly doesn't follow the triangle inequality either).

Answer 3

A quick note: I haven't tested any of the code in this answer, so there may be bugs hiding in there.

OK, first things first. Let's figure out what we're looking at:

You want to write a function that solves a given instance of the traveling salesman problem. The algorithm you've decided to use is to enumerate the permutations of cities, terminating consideration of a permutation prefix if its partial cost exceeds the cost of the best full permutation so far.

I'll touch more on the choice of algorithms below, but let's start by writing a class to represent an instance of the problem.

/**
 * Represents an instance of the traveling salesman problem.
 *
 * Any solution must start and finish at city 0.
 */
class TravelingSalesmanProblem {
  public:
    template <std::size_t N>
    TravelingSalesmanProblem(const int (&distance_matrix)[N][N]);

    int distance(std::size_t city1, std::size_t city2) const;
    std::size_t num_cities() const { return num_cities_; }

  private:
    std::vector<std::vector<int>> distances_;
    std::size_t num_cities_;
};

One reason to write a class instead of passing around handles to the distance matrix directly is that I can imagine two reasonable implementations of this class, and I'd need a profiler to choose one over the other.

Version 1:

template <std::size_t N>
TravelingSalesmanProblem::TravelingSalesmanProblem(const int (&distance_matrix)[N][N]) 
    : distances_(N, std::vector<int>(N)), num_cities_(N) {
  for (std::size_t i = 0; i < N; ++i) {
    for (std::size_t j = i; j < N; ++j) {
      assert(distance_matrix[i][j] == distance_matrix[j][i]);  // You might rewrite this using exceptions instead.
      distances_[i][j] = distances_[j][i] = distance_matrix[i][j];
    }
  }
}

int TravelingSalesmanProblem::distance(std::size_t city1, std::size_t city2) const {
  return distances_.at(city1).at(city2);
}

Version 2 requires changing the declaration of TravelingSalesmanProblem::distance_ to

std::vector<int> distance_;

Then the implementation of TravelingSalesmanProblem is:

template <std::size_t N>
TravelingSalesmanProblem::TravelingSalesmanProblem(const int (&distance_matrix)[N][N])
    : distances_(N * (N + 1) / 2), num_cities_(N) {
  for (std::size_t i = 0; i < N; ++i) {
    for (std::size_t j = 0; j <= i; ++j) {
      assert(distance_matrix[i][j] == distance_matrix[j][i]);
      distances_[i * (i + 1) / 2 + j] = distance_matrix[i][j];
    }
  }
}

int TravelingSalesmanProblem::distance(std::size_t city1, std::size_t city2) const {
  std::size_t i, j;
  std::tie(j, i) = std::minmax(city1, city2);
  return distances_.at(i * (i + 1) / 2 + j);
}

Version 2 requires less indirection and less memory but requires a little extra computation at each call to distance. You might also make TravelingSalesmanProblem a template class parameterized on N so that you can make distance_ an int[N][N], but this will require that the problem size is known at compile time (which is fine for this example but makes the code less reusable).

Now let's write solve_tsp (I'm writing it here as a standalone function, but you could reasonably make it a member of TravelingSalesmanProblem).

The first question is which data structures to use for various purposes. The first is the path concept. Here, I suggest std::vector<int>. You used QList<int>, but that includes memory indirection and likely non-local allocations of elements in memory, which I think may be a significant factor in running time. std::vector allocates a single block of contiguous memory where its members are stored; this means that access to elements requires only one memory indirection (as opposed to two for QList) and that contiguous elements are contiguous in memory as well, which may improve memory caching.

You also used QList<int> for the list of remaining cities. QList doesn't really buy you anything over std::vector here, since you never prepend; see also my arguments about path. Therefore, I also choose std::vector<int> for this purpose.

solve_tsp is just a sort of wrapper function that initializes the data structures that we'll need and calls out to a helper, which uses recursion to iterate over all permutations of the cities.

std::vector<int> solve_tsp(const TravelingSalesmanProblem& problem) {
  std::vector<int> path;
  path.reserve(problem.num_cities() + 1);
  path.push_back(0);

  std::vector<int> remaining_cities(problem.num_cities() - 1);
  std::iota(remaining_cities.begin(), remaining_cities.end(), 1);

  solve_tsp_helper(path, remaining_cities, 0, std::numeric_limits<int>::max());

  return path;
}

// path is an in-out parameter; if there is any solution with path as a prefix whose
// cost is less than best_previous_path, on return path will contain the best such
// solution, and the return value will be the cost of that solution. If no such
// solution exists, path will be unmodified after return and the return value will
// be at least as large as best_previous_path.
int solve_tsp_helper(const TravelingSalesmanProblem& problem,
                     std::vector<int>& path,
                     const std::vector<int>& remaining_cities,
                     const int cost_so_far, int best_previous_path_cost) {
  const int last_city = path.back();
  const auto current_path_length = path.size();
  const int full_path_cost = cost_so_far + problem.distance(last_city, 0);
  if (remaining_cities.empty()) {
    path.push_back(0);
    return full_path_cost;
  }
  // Note: The validity of this early exit relies on problem adhering to the triangle inequality.
  if (full_path_cost >= best_previous_path_cost) {
    return full_path_cost;
  }
  auto best_path = path;
  // This next bit of code is a bit more complicated than the naive
  // "for (int next_city : remaining_cities)" to prevent copying
  // remaining_cities more than once.
  std::vector<int> other_remaining_cities(std::next(remaining_cities.begin()),
                                          remaining_cities.end());
  for (std::size_t i = 0; i < remaining_cities.size(); ++i) {
    // This loop maintains the invariant that either best_previous_path_cost has its 
    // original value and best_path has path's original value, or 
    // best_previous_path_cost has reduced in value and best_path is a full path
    // with cost best_previous_path_cost.
    const int next_city = remaining_cities[i];

    if (i != 0) {
      other_remaining_cities[i - 1] = remaining_cities[i - 1];
    }

    path.resize(current_path_length);
    path.push_back(next_city);

    const auto cost = solve_tsp_helper(
        problem, path, other_remaining_cities,
        cost_so_far + problem.distance(last_city, next_city), best_previous_path_cost);
    if (cost < best_previous_path) {
      best_previous_path_cost = cost;
      swap(path, best_path);
    }
  }
  // Because of the loop invariant, either best_previous_path_cost has not changed
  // and best_path is the original value of path, or best_path is a solution to problem
  // with cost equal to best_previous_path_cost. In either case, we fulfill the
  // contract of the function by swapping best_path into path and returning
  // best_previous_path_cost.
  swap(path, best_path);
  return best_previous_path_cost;
}

Some notes on this implementation. First of all, it avoids any calls to new or delete. It also copies remaining_cities and path once per call to solve_tsp_helper instead of once per element of remaining_cities. So even though the algorithm hasn't changed (much -- see next paragraph), the code should be faster. With some effort, you could probably eliminate any copies of path except when a better path is found.

I've also slightly changed your early exit condition, from cost_so_far > best_previous_path_cost to cost_so_far + distance(last_city, 0) >= best_previous_path_cost. This is valid whenever the triangle inequality holds for your data (that is, for any three cities A, B, C, dist(A, C) <= dist(A, B) + dist(B, C); this condition usually holds).

As for the algorithm, this SO question discusses faster exact algorithms. They are considerably more complicated than this algorithm but can yield much better performance. You can also consider heuristic answers, which can produce reasonably good solutions though are not guaranteed to provide the exact optimal solution.