# HG changeset patch
# User Peter Kovacs <kpeter@inf.elte.hu>
# Date 1294582015 -3600
# Node ID dff32ce3db712f2d6e262e75fae4d9404798afa8
# Parent  07682e24c4e8062ec5e7b937652ce9bcf82cbf7f
Make InsertionTsp much faster and improve docs (#386)

diff -r 07682e24c4e8 -r dff32ce3db71 doc/groups.dox
--- a/doc/groups.dox	Sun Jan 09 00:57:12 2011 +0100
+++ b/doc/groups.dox	Sun Jan 09 15:06:55 2011 +0100
@@ -572,6 +572,16 @@
  - \ref ChristofidesTsp Christofides algorithm
  - \ref Opt2Tsp 2-opt algorithm
 
+\ref NearestNeighborTsp, \ref GreedyTsp, and \ref InsertionTsp are the fastest
+solution methods. Furthermore, \ref InsertionTsp is usually quite effective.
+
+\ref ChristofidesTsp is somewhat slower, but it has the best guaranteed
+approximation factor: 3/2.
+
+\ref Opt2Tsp usually provides the best results in practice, but
+it is the slowest method. It can also be used to improve given tours,
+for example, the results of other algorithms.
+
 \image html tsp.png
 \image latex tsp.eps "Traveling salesman problem" width=\textwidth
 */
diff -r 07682e24c4e8 -r dff32ce3db71 lemon/christofides_tsp.h
--- a/lemon/christofides_tsp.h	Sun Jan 09 00:57:12 2011 +0100
+++ b/lemon/christofides_tsp.h	Sun Jan 09 15:06:55 2011 +0100
@@ -40,8 +40,9 @@
   ///
   /// This a well-known approximation method for the TSP problem with
   /// metric cost function.
-  /// It yields a tour whose total cost is at most 3/2 of the optimum,
-  /// but it is usually much better.
+  /// It has a guaranteed approximation factor of 3/2 (i.e. it finds a tour
+  /// whose total cost is at most 3/2 of the optimum), but it usually
+  /// provides better solutions in practice.
   /// This implementation runs in O(n<sup>3</sup>log(n)) time.
   ///
   /// The algorithm starts with a \ref spantree "minimum cost spanning tree" and
diff -r 07682e24c4e8 -r dff32ce3db71 lemon/greedy_tsp.h
--- a/lemon/greedy_tsp.h	Sun Jan 09 00:57:12 2011 +0100
+++ b/lemon/greedy_tsp.h	Sun Jan 09 15:06:55 2011 +0100
@@ -43,9 +43,10 @@
   /// as long as it does not create a cycle of less than n edges and it does
   /// not increase the degree of any node above two.
   ///
-  /// This method runs in O(n<sup>2</sup>log(n)) time.
-  /// It quickly finds a short tour for most TSP instances, but in special
-  /// cases, it could yield a really bad (or even the worst) solution.
+  /// This method runs in O(n<sup>2</sup>) time.
+  /// It quickly finds a relatively short tour for most TSP instances,
+  /// but it could also yield a really bad (or even the worst) solution
+  /// in special cases.
   ///
   /// \tparam CM Type of the cost map.
   template <typename CM>
diff -r 07682e24c4e8 -r dff32ce3db71 lemon/insertion_tsp.h
--- a/lemon/insertion_tsp.h	Sun Jan 09 00:57:12 2011 +0100
+++ b/lemon/insertion_tsp.h	Sun Jan 09 15:06:55 2011 +0100
@@ -24,6 +24,7 @@
 /// \brief Insertion algorithm for symmetric TSP
 
 #include <vector>
+#include <functional>
 #include <lemon/full_graph.h>
 #include <lemon/maps.h>
 #include <lemon/random.h>
@@ -37,13 +38,20 @@
   /// InsertionTsp implements the insertion heuristic for solving
   /// symmetric \ref tsp "TSP".
   ///
-  /// This is a basic TSP heuristic that has many variants.
+  /// This is a fast and effective tour construction method that has
+  /// many variants.
   /// It starts with a subtour containing a few nodes of the graph and it
   /// iteratively inserts the other nodes into this subtour according to a
   /// certain node selection rule.
   ///
-  /// This implementation provides four different node selection rules,
-  /// from which the most powerful one is used by default.
+  /// This method is among the fastest TSP algorithms, and it typically
+  /// provides quite good solutions (usually much better than
+  /// \ref NearestNeighborTsp and \ref GreedyTsp).
+  ///
+  /// InsertionTsp implements four different node selection rules,
+  /// from which the most effective one (\e farthest \e node \e selection)
+  /// is used by default.
+  /// With this choice, the algorithm runs in O(n<sup>2</sup>) time.
   /// For more information, see \ref SelectionRule.
   ///
   /// \tparam CM Type of the cost map.
@@ -64,7 +72,7 @@
       const FullGraph &_gr;
       const CostMap &_cost;
       std::vector<Node> _notused;
-      std::vector<Node> _path;
+      std::vector<Node> _tour;
       Cost _sum;
 
     public:
@@ -76,9 +84,10 @@
       ///
       /// During the algorithm, nodes are selected for addition to the current
       /// subtour according to the applied rule.
-      /// In general, the FARTHEST method yields the best tours, thus it is the
-      /// default option. The RANDOM rule usually gives somewhat worse results,
-      /// but it is much faster than the others and it is the most robust.
+      /// The FARTHEST method is one of the fastest selection rules, and
+      /// it is typically the most effective, thus it is the default
+      /// option. The RANDOM rule usually gives slightly worse results,
+      /// but it is more robust.
       ///
       /// The desired selection rule can be specified as a parameter of the
       /// \ref run() function.
@@ -86,20 +95,21 @@
 
         /// An unvisited node having minimum distance from the current
         /// subtour is selected at each step.
-        /// The algorithm runs in O(n<sup>3</sup>) time using this
+        /// The algorithm runs in O(n<sup>2</sup>) time using this
         /// selection rule.
         NEAREST,
 
         /// An unvisited node having maximum distance from the current
         /// subtour is selected at each step.
-        /// The algorithm runs in O(n<sup>3</sup>) time using this
+        /// The algorithm runs in O(n<sup>2</sup>) time using this
         /// selection rule.
         FARTHEST,
 
         /// An unvisited node whose insertion results in the least
         /// increase of the subtour's total cost is selected at each step.
         /// The algorithm runs in O(n<sup>3</sup>) time using this
-        /// selection rule.
+        /// selection rule, but in most cases, it is almost as fast as
+        /// with other rules.
         CHEAPEST,
 
         /// An unvisited node is selected randomly without any evaluation
@@ -134,22 +144,24 @@
       ///
       /// \return The total cost of the found tour.
       Cost run(SelectionRule rule = FARTHEST) {
-        _path.clear();
+        _tour.clear();
 
         if (_gr.nodeNum() == 0) return _sum = 0;
         else if (_gr.nodeNum() == 1) {
-          _path.push_back(_gr(0));
+          _tour.push_back(_gr(0));
           return _sum = 0;
         }
 
         switch (rule) {
           case NEAREST:
             init(true);
-            start<NearestSelection, DefaultInsertion>();
+            start<ComparingSelection<std::less<Cost> >,
+                  DefaultInsertion>();
             break;
           case FARTHEST:
             init(false);
-            start<FarthestSelection, DefaultInsertion>();
+            start<ComparingSelection<std::greater<Cost> >,
+                  DefaultInsertion>();
             break;
           case CHEAPEST:
             init(true);
@@ -185,7 +197,7 @@
       ///
       /// \pre run() must be called before using this function.
       const std::vector<Node>& tourNodes() const {
-        return _path;
+        return _tour;
       }
 
       /// \brief Gives back the node sequence of the found tour.
@@ -196,7 +208,7 @@
       /// \pre run() must be called before using this function.
       template <typename Container>
       void tourNodes(Container &container) const {
-        container.assign(_path.begin(), _path.end());
+        container.assign(_tour.begin(), _tour.end());
       }
 
       /// \brief Gives back the found tour as a path.
@@ -208,11 +220,11 @@
       template <typename Path>
       void tour(Path &path) const {
         path.clear();
-        for (int i = 0; i < int(_path.size()) - 1; ++i) {
-          path.addBack(_gr.arc(_path[i], _path[i+1]));
+        for (int i = 0; i < int(_tour.size()) - 1; ++i) {
+          path.addBack(_gr.arc(_tour[i], _tour[i+1]));
         }
-        if (int(_path.size()) >= 2) {
-          path.addBack(_gr.arc(_path.back(), _path.front()));
+        if (int(_tour.size()) >= 2) {
+          path.addBack(_gr.arc(_tour.back(), _tour.front()));
         }
       }
 
@@ -224,9 +236,9 @@
       void init(bool min) {
         Edge min_edge = min ? mapMin(_gr, _cost) : mapMax(_gr, _cost);
 
-        _path.clear();
-        _path.push_back(_gr.u(min_edge));
-        _path.push_back(_gr.v(min_edge));
+        _tour.clear();
+        _tour.push_back(_gr.u(min_edge));
+        _tour.push_back(_gr.v(min_edge));
 
         _notused.clear();
         for (NodeIt n(_gr); n!=INVALID; ++n) {
@@ -241,106 +253,82 @@
       // Executes the algorithm
       template <class SelectionFunctor, class InsertionFunctor>
       void start() {
-        SelectionFunctor selectNode(_gr, _cost, _path, _notused);
-        InsertionFunctor insertNode(_gr, _cost, _path, _sum);
+        SelectionFunctor selectNode(_gr, _cost, _tour, _notused);
+        InsertionFunctor insertNode(_gr, _cost, _tour, _sum);
 
         for (int i=0; i<_gr.nodeNum()-2; ++i) {
           insertNode.insert(selectNode.select());
         }
 
-        _sum = _cost[_gr.edge(_path.back(), _path.front())];
-        for (int i = 0; i < int(_path.size())-1; ++i) {
-          _sum += _cost[_gr.edge(_path[i], _path[i+1])];
+        _sum = _cost[_gr.edge(_tour.back(), _tour.front())];
+        for (int i = 0; i < int(_tour.size())-1; ++i) {
+          _sum += _cost[_gr.edge(_tour[i], _tour[i+1])];
         }
       }
 
 
-      // Implementation of the nearest selection rule
-      class NearestSelection {
+      // Implementation of the nearest and farthest selection rule
+      template <typename Comparator>
+      class ComparingSelection {
         public:
-          NearestSelection(const FullGraph &gr, const CostMap &cost,
-                           std::vector<Node> &path, std::vector<Node> &notused)
-            : _gr(gr), _cost(cost), _path(path), _notused(notused) {}
-
-          Node select() const {
-            Cost insert_val = 0;
-            int insert_node = -1;
-
+          ComparingSelection(const FullGraph &gr, const CostMap &cost,
+                  std::vector<Node> &tour, std::vector<Node> &notused)
+            : _gr(gr), _cost(cost), _tour(tour), _notused(notused),
+              _dist(gr, 0), _compare()
+          {
+            // Compute initial distances for the unused nodes
             for (unsigned int i=0; i<_notused.size(); ++i) {
-              Cost min_val = _cost[_gr.edge(_notused[i], _path[0])];
-              int min_node = 0;
-
-              for (unsigned int j=1; j<_path.size(); ++j) {
-                Cost curr = _cost[_gr.edge(_notused[i], _path[j])];
-                if (min_val > curr) {
-                    min_val = curr;
-                    min_node = j;
+              Node u = _notused[i];
+              Cost min_dist = _cost[_gr.edge(u, _tour[0])];
+              for (unsigned int j=1; j<_tour.size(); ++j) {
+                Cost curr = _cost[_gr.edge(u, _tour[j])];
+                if (curr < min_dist) {
+                  min_dist = curr;
                 }
               }
+              _dist[u] = min_dist;
+            }
+          }
 
-              if (insert_val > min_val || insert_node == -1) {
-                insert_val = min_val;
-                insert_node = i;
+          Node select() {
+
+            // Select an used node with minimum distance
+            Cost ins_dist = 0;
+            int ins_node = -1;
+            for (unsigned int i=0; i<_notused.size(); ++i) {
+              Cost curr = _dist[_notused[i]];
+              if (_compare(curr, ins_dist) || ins_node == -1) {
+                ins_dist = curr;
+                ins_node = i;
               }
             }
 
-            Node n = _notused[insert_node];
-            _notused.erase(_notused.begin()+insert_node);
+            // Remove the selected node from the unused vector
+            Node sn = _notused[ins_node];
+            _notused[ins_node] = _notused.back();
+            _notused.pop_back();
 
-            return n;
+            // Update the distances of the remaining nodes
+            for (unsigned int i=0; i<_notused.size(); ++i) {
+              Node u = _notused[i];
+              Cost nc = _cost[_gr.edge(sn, u)];
+              if (nc < _dist[u]) {
+                _dist[u] = nc;
+              }
+            }
+
+            return sn;
           }
 
         private:
           const FullGraph &_gr;
           const CostMap &_cost;
-          std::vector<Node> &_path;
+          std::vector<Node> &_tour;
           std::vector<Node> &_notused;
+          FullGraph::NodeMap<Cost> _dist;
+          Comparator _compare;
       };
 
-
-      // Implementation of the farthest selection rule
-      class FarthestSelection {
-        public:
-          FarthestSelection(const FullGraph &gr, const CostMap &cost,
-                            std::vector<Node> &path, std::vector<Node> &notused)
-            : _gr(gr), _cost(cost), _path(path), _notused(notused) {}
-
-          Node select() const {
-            Cost insert_val = 0;
-            int insert_node = -1;
-
-            for (unsigned int i=0; i<_notused.size(); ++i) {
-              Cost min_val = _cost[_gr.edge(_notused[i], _path[0])];
-              int min_node = 0;
-
-              for (unsigned int j=1; j<_path.size(); ++j) {
-                Cost curr = _cost[_gr.edge(_notused[i], _path[j])];
-                if (min_val > curr) {
-                  min_val = curr;
-                  min_node = j;
-                }
-              }
-
-              if (insert_val < min_val || insert_node == -1) {
-                insert_val = min_val;
-                insert_node = i;
-              }
-            }
-
-            Node n = _notused[insert_node];
-            _notused.erase(_notused.begin()+insert_node);
-
-            return n;
-          }
-
-        private:
-          const FullGraph &_gr;
-          const CostMap &_cost;
-          std::vector<Node> &_path;
-          std::vector<Node> &_notused;
-      };
-
-
       // Implementation of the cheapest selection rule
       class CheapestSelection {
         private:
@@ -353,50 +341,102 @@
 
         public:
           CheapestSelection(const FullGraph &gr, const CostMap &cost,
-                            std::vector<Node> &path, std::vector<Node> &notused)
-            : _gr(gr), _cost(cost), _path(path), _notused(notused) {}
-
-          Cost select() const {
-            int insert_notused = -1;
-            int best_insert_index = -1;
-            Cost insert_val = 0;
-
+                            std::vector<Node> &tour, std::vector<Node> &notused)
+            : _gr(gr), _cost(cost), _tour(tour), _notused(notused),
+              _ins_cost(gr, 0), _ins_pos(gr, -1)
+          {
+            // Compute insertion cost and position for the unused nodes
             for (unsigned int i=0; i<_notused.size(); ++i) {
-              int min = i;
-              int best_insert_tmp = 0;
-              Cost min_val =
-                costDiff(_path.front(), _path.back(), _notused[i]);
-
-              for (unsigned int j=1; j<_path.size(); ++j) {
-                Cost tmp =
-                  costDiff(_path[j-1], _path[j], _notused[i]);
-
-                if (min_val > tmp) {
-                  min = i;
-                  min_val = tmp;
-                  best_insert_tmp = j;
+              Node u = _notused[i];
+              Cost min_cost = costDiff(_tour.back(), _tour.front(), u);
+              int min_pos = 0;              
+              for (unsigned int j=1; j<_tour.size(); ++j) {
+                Cost curr_cost = costDiff(_tour[j-1], _tour[j], u);
+                if (curr_cost < min_cost) {
+                  min_cost = curr_cost;
+                  min_pos = j;
                 }
               }
+              _ins_cost[u] = min_cost;
+              _ins_pos[u] = min_pos;
+            }
+          }
 
-              if (insert_val > min_val || insert_notused == -1) {
-                insert_notused = min;
-                insert_val = min_val;
-                best_insert_index = best_insert_tmp;
+          Cost select() {
+
+            // Select an used node with minimum insertion cost
+            Cost min_cost = 0;
+            int min_node = -1;
+            for (unsigned int i=0; i<_notused.size(); ++i) {
+              Cost curr_cost = _ins_cost[_notused[i]];
+              if (curr_cost < min_cost || min_node == -1) {
+                min_cost = curr_cost;
+                min_node = i;
               }
             }
 
-            _path.insert(_path.begin()+best_insert_index,
-                         _notused[insert_notused]);
-            _notused.erase(_notused.begin()+insert_notused);
+            // Remove the selected node from the unused vector
+            Node sn = _notused[min_node];
+            _notused[min_node] = _notused.back();
+            _notused.pop_back();
+            
+            // Insert the selected node into the tour
+            const int ipos = _ins_pos[sn];
+            _tour.insert(_tour.begin() + ipos, sn);
 
-            return insert_val;
+            // Update the insertion cost and position of the remaining nodes
+            for (unsigned int i=0; i<_notused.size(); ++i) {
+              Node u = _notused[i];
+              Cost curr_cost = _ins_cost[u];
+              int curr_pos = _ins_pos[u];
+
+              int ipos_prev = ipos == 0 ? _tour.size()-1 : ipos-1;
+              int ipos_next = ipos == int(_tour.size())-1 ? 0 : ipos+1;
+              Cost nc1 = costDiff(_tour[ipos_prev], _tour[ipos], u);
+              Cost nc2 = costDiff(_tour[ipos], _tour[ipos_next], u);
+              
+              if (nc1 <= curr_cost || nc2 <= curr_cost) {
+                // A new position is better than the old one
+                if (nc1 <= nc2) {
+                  curr_cost = nc1;
+                  curr_pos = ipos;
+                } else {
+                  curr_cost = nc2;
+                  curr_pos = ipos_next;
+                }
+              }
+              else {
+                if (curr_pos == ipos) {
+                  // The minimum should be found again
+                  curr_cost = costDiff(_tour.back(), _tour.front(), u);
+                  curr_pos = 0;              
+                  for (unsigned int j=1; j<_tour.size(); ++j) {
+                    Cost tmp_cost = costDiff(_tour[j-1], _tour[j], u);
+                    if (tmp_cost < curr_cost) {
+                      curr_cost = tmp_cost;
+                      curr_pos = j;
+                    }
+                  }
+                }
+                else if (curr_pos > ipos) {
+                  ++curr_pos;
+                }
+              }
+              
+              _ins_cost[u] = curr_cost;
+              _ins_pos[u] = curr_pos;
+            }
+
+            return min_cost;
           }
 
         private:
           const FullGraph &_gr;
           const CostMap &_cost;
-          std::vector<Node> &_path;
+          std::vector<Node> &_tour;
           std::vector<Node> &_notused;
+          FullGraph::NodeMap<Cost> _ins_cost;
+          FullGraph::NodeMap<int> _ins_pos;
       };
 
       // Implementation of the random selection rule
@@ -409,9 +449,11 @@
           Node select() const {
             const int index = rnd[_notused.size()];
             Node n = _notused[index];
-            _notused.erase(_notused.begin()+index);
+            _notused[index] = _notused.back();
+            _notused.pop_back();
             return n;
           }
+
         private:
           std::vector<Node> &_notused;
       };
@@ -429,30 +471,30 @@
 
         public:
           DefaultInsertion(const FullGraph &gr, const CostMap &cost,
-                           std::vector<Node> &path, Cost &total_cost) :
-            _gr(gr), _cost(cost), _path(path), _total(total_cost) {}
+                           std::vector<Node> &tour, Cost &total_cost) :
+            _gr(gr), _cost(cost), _tour(tour), _total(total_cost) {}
 
           void insert(Node n) const {
             int min = 0;
             Cost min_val =
-              costDiff(_path.front(), _path.back(), n);
+              costDiff(_tour.front(), _tour.back(), n);
 
-            for (unsigned int i=1; i<_path.size(); ++i) {
-              Cost tmp = costDiff(_path[i-1], _path[i], n);
+            for (unsigned int i=1; i<_tour.size(); ++i) {
+              Cost tmp = costDiff(_tour[i-1], _tour[i], n);
               if (tmp < min_val) {
                 min = i;
                 min_val = tmp;
               }
             }
 
-            _path.insert(_path.begin()+min, n);
+            _tour.insert(_tour.begin()+min, n);
             _total += min_val;
           }
 
         private:
           const FullGraph &_gr;
           const CostMap &_cost;
-          std::vector<Node> &_path;
+          std::vector<Node> &_tour;
           Cost &_total;
       };
 
diff -r 07682e24c4e8 -r dff32ce3db71 lemon/nearest_neighbor_tsp.h
--- a/lemon/nearest_neighbor_tsp.h	Sun Jan 09 00:57:12 2011 +0100
+++ b/lemon/nearest_neighbor_tsp.h	Sun Jan 09 15:06:55 2011 +0100
@@ -44,8 +44,9 @@
   /// Finally, it connects the two end points of the path to form a tour.
   ///
   /// This method runs in O(n<sup>2</sup>) time.
-  /// It quickly finds a short tour for most TSP instances, but in special
-  /// cases, it could yield a really bad (or even the worst) solution.
+  /// It quickly finds a relatively short tour for most TSP instances,
+  /// but it could also yield a really bad (or even the worst) solution
+  /// in special cases.
   ///
   /// \tparam CM Type of the cost map.
   template <typename CM>
diff -r 07682e24c4e8 -r dff32ce3db71 lemon/opt2_tsp.h
--- a/lemon/opt2_tsp.h	Sun Jan 09 00:57:12 2011 +0100
+++ b/lemon/opt2_tsp.h	Sun Jan 09 15:06:55 2011 +0100
@@ -45,9 +45,9 @@
   /// algorithm uses the node sequence determined by the node IDs.
   /// Oherwise, it starts with the given tour.
   ///
-  /// This is a relatively slow but powerful method. 
-  /// A typical usage of it is the improvement of a solution that is resulted
-  /// by a fast tour construction heuristic (e.g. the InsertionTsp algorithm).
+  /// This is a rather slow but effective method.
+  /// Its typical usage is the improvement of the result of a fast tour
+  /// construction heuristic (e.g. the InsertionTsp algorithm).
   ///
   /// \tparam CM Type of the cost map.
   template <typename CM>