Ticket #386: 386-tsp-7ed08f6da055.patch

lemon/Makefile.am

@@ -64 +64 @@
         lemon/bucket_heap.h \
         lemon/capacity_scaling.h \
         lemon/cbc.h \
+        lemon/christofides_tsp.h \
         lemon/circulation.h \
         lemon/clp.h \
         lemon/color.h \
@@ -89 +90 @@
         lemon/glpk.h \
         lemon/gomory_hu.h \
         lemon/graph_to_eps.h \
+        lemon/greedy_tsp.h \
         lemon/grid_graph.h \
         lemon/hartmann_orlin_mmc.h \
         lemon/howard_mmc.h \
         lemon/hypercube_graph.h \
+        lemon/insertion_tsp.h \
         lemon/karp_mmc.h \
         lemon/kruskal.h \
         lemon/hao_orlin.h \
@@ -107 +110 @@
         lemon/math.h \
         lemon/min_cost_arborescence.h \
         lemon/nauty_reader.h \
+        lemon/nearest_neighbor_tsp.h \
         lemon/network_simplex.h \
+        lemon/opt2_tsp.h \
         lemon/pairing_heap.h \
         lemon/path.h \
         lemon/planarity.h \

new file lemon/christofides_tsp.h

@@ +1 @@
+#ifndef LEMON_CHRISTOFIDES_TSP
+#define LEMON_CHRISTOFIDES_TSP
+
+#include <lemon/full_graph.h>
+#include <lemon/smart_graph.h>
+#include <lemon/path.h>
+#include <lemon/kruskal.h>
+#include <lemon/matching.h>
+#include <lemon/adaptors.h>
+#include <lemon/maps.h>
+#include <lemon/euler.h>
+#include <lemon/tsp_utils.h>
+
+namespace lemon {
+  
+  namespace christofides_helper {
+    template <class L>
+    L vectorConvert(const std::vector<FullGraph::Node> &_path) {
+      return L(_path.begin(), _path.end());
+    }
+  
+    template <>
+    std::vector<FullGraph::Node> vectorConvert(const std::vector<FullGraph::Node> &_path) {
+      return _path;
+    }
+  }
+  
+  template <typename CM>
+  class ChristofidesTsp {
+    private:
+      GRAPH_TYPEDEFS(SmartGraph);
+
+    public:
+      typedef typename CM::Value Cost;
+      typedef SmartGraph::EdgeMap<Cost> CostMap;
+    
+      ChristofidesTsp(const FullGraph &gr, const CM &cost) : _cost(_gr), fullg(gr), fullcost(cost), nr(_gr) {
+        graphCopy(gr, _gr).nodeCrossRef(nr).edgeMap(cost, _cost).run();
+      }
+
+      Cost run() {
+        _path.clear();
+        
+        SmartGraph::EdgeMap<bool> tree(_gr);
+        kruskal(_gr, _cost, tree);
+        
+        FilterEdges<SmartGraph> treefiltered(_gr, tree);
+        InDegMap<FilterEdges<SmartGraph> > deg(treefiltered);
+        
+        SmartGraph::NodeMap<bool> oddfilter(_gr, false);
+        FilterNodes<SmartGraph> oddNodes(_gr, oddfilter);
+        
+        for (NodeIt n(_gr); n!=INVALID; ++n) {
+          if (deg[n]%2 == 1) {
+            oddNodes.enable(n);
+          }
+        }
+        
+        NegMap<CostMap> negmap(_cost);
+        MaxWeightedPerfectMatching<FilterNodes<SmartGraph>,
+                  NegMap<CostMap> > perfmatch(oddNodes, negmap);
+        perfmatch.run();
+        
+        for (FilterNodes<SmartGraph>::EdgeIt e(oddNodes); e!=INVALID; ++e) {
+          if (perfmatch.matching(e)) {
+            treefiltered.enable(_gr.addEdge(_gr.u(e), _gr.v(e)));
+          }
+        }
+        
+        FilterEdges<SmartGraph>::NodeMap<bool> seen(treefiltered, false);
+        for (EulerIt<FilterEdges<SmartGraph> > e(treefiltered); e!=INVALID; ++e) {
+          if (seen[treefiltered.target(e)]==false) {
+            _path.push_back(nr[treefiltered.target(e)]);
+            seen[treefiltered.target(e)] = true;
+          }
+        }
+
+        _sum = fullcost[ fullg.edge(_path.back(), _path.front()) ];
+        for (unsigned int i=0; i<_path.size()-1; ++i)
+          _sum += fullcost[ fullg.edge(_path[i], _path[i+1]) ];
+
+        return _sum;
+      }
+
+      template <typename L>
+      void tourNodes(L &container) {
+        container(christofides_helper::vectorConvert<L>(_path));
+      }
+      
+      template <template <typename> class L>
+      L<FullGraph::Node> tourNodes() {
+        return christofides_helper::vectorConvert<L<FullGraph::Node> >(_path);
+      }
+
+      const std::vector<Node>& tourNodes() {
+        return _path;
+      }
+      
+      Path<FullGraph> tour() {
+        Path<FullGraph> p;
+        if (_path.size()<2)
+          return p;
+
+        for (unsigned int i=0; i<_path.size()-1; ++i) {
+          p.addBack(fullg.arc(_path[i], _path[i+1]));
+        }
+        p.addBack(fullg.arc(_path.back(), _path.front()));
+        
+        return p;
+      }
+      
+      Cost tourCost() {
+        return _sum;
+      }
+      
+
+  private:
+    SmartGraph _gr;
+    CostMap _cost;
+    Cost _sum;
+    const FullGraph &fullg;
+    const CM &fullcost;
+    std::vector<FullGraph::Node> _path;
+    SmartGraph::NodeMap<FullGraph::Node> nr;
+  };
+
+
+}; // namespace lemon
+
+#endif

new file lemon/greedy_tsp.h

@@ +1 @@
+#ifndef LEMON_GREEDY_TSP
+#define LEMON_GREEDY_TSP
+
+#include <lemon/path.h>
+#include <lemon/full_graph.h>
+#include <lemon/unionfind.h>
+#include <lemon/tsp_utils.h>
+#include <algorithm>
+
+namespace lemon {
+
+  namespace greedy_tsp_helper {
+
+    template <typename CostMap>
+    class KeyComp {
+      typedef typename CostMap::Key Key;
+      const CostMap &cost;
+      
+      public:
+        KeyComp(const CostMap &_cost) : cost(_cost) {}
+      
+        bool operator() (const Key &a, const Key &b) const {
+          return cost[a] < cost[b];
+        }
+    };
+
+  
+    template <class L>
+    L vectorConvert(const std::vector<FullGraph::Node> &_path) {
+      return L(_path.begin(), _path.end());
+    }
+  
+    template <>
+    std::vector<FullGraph::Node> vectorConvert(const std::vector<FullGraph::Node> &_path) {
+      return _path;
+    }
+    
+  }
+
+
+  template <typename CM>
+  class GreedyTsp {
+    private:
+      GRAPH_TYPEDEFS(FullGraph);
+    
+    public:
+      typedef CM CostMap;
+      typedef typename CM::Value Cost;
+      
+      GreedyTsp(const FullGraph &gr, const CostMap &cost) : _gr(gr), _cost(cost) {}  
+
+      Cost run() {
+        typedef UnionFind<FullGraph::NodeMap<int> > Union;
+        _nodes.clear();
+        
+        std::vector<int> path;
+        path.resize(_gr.nodeNum()*2, -1);
+        
+        std::vector<typename CostMap::Key> sorted_edges;
+        sorted_edges.reserve(_gr.edgeNum());
+        for (EdgeIt n(_gr); n != INVALID; ++n)
+          sorted_edges.push_back(n);
+
+        std::sort(sorted_edges.begin(), sorted_edges.end(), greedy_tsp_helper::KeyComp<CostMap>(_cost));
+
+        FullGraph::NodeMap<int> nodemap(_gr);
+        Union unionfind(nodemap);
+
+        for (NodeIt n(_gr); n != INVALID; ++n)
+          unionfind.insert(n);
+        
+        FullGraph::NodeMap<int> degree(_gr, 0);
+
+        int nodesNum = 0, i = 0;
+
+        while ( nodesNum != _gr.nodeNum()-1 ) {
+          const Edge &e = sorted_edges[i];
+          
+          const Node u = _gr.u(e),
+                     v = _gr.v(e);
+          
+          if (degree[u]<=1 && degree[v]<=1) {
+            if (unionfind.join(u, v)) {
+              ++degree[u];
+              ++degree[v];
+              ++nodesNum;
+              
+              const int uid = _gr.id(u),
+                        vid = _gr.id(v);
+              
+              
+              path[uid*2 + (path[uid*2]==-1 ? 0 : 1)] = vid;
+              path[vid*2 + (path[vid*2]==-1 ? 0 : 1)] = uid;
+            }
+          }
+
+          ++i;
+        }
+
+
+        for (int i=0, n=-1; i<_gr.nodeNum()*2; ++i) {
+          if (path[i] == -1) {
+            if (n==-1) {
+              n = i;
+            } else {
+              path[n] = i/2;
+              path[i] = n/2;
+              break;
+            }
+          }
+        }
+
+
+        for (int i=0, j=0, last=-1; i!=_gr.nodeNum(); ++i) {
+          _nodes.push_back(_gr.nodeFromId(j));
+          
+          if (path[2*j] != last) {
+            last = j;
+            j = path[2*j];
+          } else {
+            last = j;
+            j = path[2*j+1];
+          }
+        }
+
+        _sum = _cost[_gr.edge(_nodes.back(), _nodes.front())];
+        for (unsigned int i=0; i<_nodes.size()-1; ++i)
+          _sum += _cost[_gr.edge(_nodes[i], _nodes[i+1])];
+
+        return _sum;
+      }
+
+
+
+      template <typename L>
+      void tourNodes(L &container) {
+        container(greedy_tsp_helper::vectorConvert<L>(_nodes));
+      }
+      
+      template <template <typename> class L>
+      L<Node> tourNodes() {
+        return greedy_tsp_helper::vectorConvert<L<Node> >(_nodes);
+      }
+      
+      const std::vector<Node>& tourNodes() {
+        return _nodes;
+      }
+      
+      Path<FullGraph> tour() {
+        Path<FullGraph> p;
+        if (_nodes.size()<2)
+          return p;
+        
+        for (unsigned int i=0; i<_nodes.size()-1; ++i) {
+          p.addBack(_gr.arc(_nodes[i], _nodes[i+1]));
+        }
+        
+        p.addBack(_gr.arc(_nodes.back(), _nodes.front()));
+        
+        return p;
+      }
+      
+      Cost tourCost() {
+        return _sum;
+      }
+      
+
+    private:
+      const FullGraph &_gr;
+      const CostMap &_cost;
+      Cost _sum;
+      std::vector<Node> _nodes;
+  };
+
+}; // namespace lemon
+
+#endif

new file lemon/insertion_tsp.h

@@ +1 @@
+#ifndef LEMON_INSERTION_TSP
+#define LEMON_INSERTION_TSP
+
+#include <lemon/full_graph.h>
+#include <lemon/path.h>
+#include <lemon/maps.h>
+#include <lemon/random.h>
+#include <lemon/tsp_utils.h>
+#include <vector>
+
+namespace lemon {
+
+  namespace insertion_tsp_helper {
+  
+    template <class L>
+    L vectorConvert(const std::vector<FullGraph::Node> &_path) {
+      return L(_path.begin(), _path.end());
+    };
+  
+    template <>
+    std::vector<FullGraph::Node> vectorConvert(
+        const std::vector<FullGraph::Node> &_path) {
+      return _path;
+    };
+    
+  };
+
+  template <typename CM>
+  class InsertionTsp {
+    private:
+      GRAPH_TYPEDEFS(FullGraph);
+
+    public:      
+      typedef CM CostMap;
+      typedef typename CM::Value Cost;
+      
+      InsertionTsp(const FullGraph &gr, const CostMap &cost) : 
+            _gr(gr), _cost(cost) {}
+      
+      enum InsertionMethod {
+        INSERT_NEAREST,
+        INSERT_FARTHEST,
+        INSERT_CHEAPEST,
+        INSERT_RANDOM
+      };     
+      
+      Cost run(InsertionMethod method = INSERT_FARTHEST) {
+        switch (method) {
+          case INSERT_NEAREST:
+            start<InitTour<true>, NearestSelection, DefaultInsert>();
+            break;
+          case INSERT_FARTHEST:
+            start<InitTour<false>, FarthestSelection, DefaultInsert>();
+            break;
+          case INSERT_CHEAPEST:
+            start<InitTour<true>, CheapestSelection, CheapestInsert>();
+            break;
+          case INSERT_RANDOM:
+            start<InitTour<true>, RandomSelection, DefaultInsert>();
+            break;
+        }
+        return sum;
+      }
+
+      template <typename L>
+      void tourNodes(L &container) {
+        container(insertion_tsp_helper::vectorConvert<L>(nodesPath));
+      }
+      
+      template <template <typename> class L>
+      L<Node> tourNodes() {
+        return insertion_tsp_helper::vectorConvert<L<Node> >(nodesPath);
+      }
+      
+      const std::vector<Node>& tourNodes() {
+        return nodesPath;
+      }
+      
+      Path<FullGraph> tour() {
+        Path<FullGraph> p;
+        if (nodesPath.size()<2)
+          return p;
+        
+        for (unsigned int i=0; i<nodesPath.size()-1; ++i)
+          p.addBack(_gr.arc(nodesPath[i], nodesPath[i+1]));
+        
+        p.addBack(_gr.arc(nodesPath.back(), nodesPath.front()));
+        return p;
+      }
+      
+      Cost tourCost() {
+        return sum;
+      }
+
+
+    private:
+
+      template <bool MIN>
+      class InitTour {
+        public:
+          InitTour(const FullGraph &gr, const CostMap &cost,
+                   std::vector<Node> &used, std::vector<Node> &notused,
+                   Cost &fullcost) :
+              _gr(gr), _cost(cost), _used(used), _notused(notused),
+              _fullcost(fullcost) {}
+
+          std::vector<Node> init() const {
+            Edge min = (MIN) ? mapMin(_gr, _cost) : mapMax(_gr, _cost);
+            std::vector<Node> path;
+            path.push_back(_gr.u(min));
+            path.push_back(_gr.v(min));
+            
+            _used.clear();
+            _notused.clear();
+            for (NodeIt n(_gr); n!=INVALID; ++n) {
+              if (n != _gr.u(min) && n != _gr.v(min)) {
+                _notused.push_back(n);
+              }
+            }
+            _used.push_back(_gr.u(min));
+            _used.push_back(_gr.v(min));
+            
+            _fullcost = _cost[min]*2;
+            return path;
+          }
+
+        private:
+          const FullGraph &_gr;
+          const CostMap &_cost;
+          std::vector<Node> &_used;
+          std::vector<Node> &_notused;
+          Cost &_fullcost;
+      };
+
+      class NearestSelection {
+        public:
+          NearestSelection(const FullGraph &gr, const CostMap &cost,
+                           std::vector<Node> &used, std::vector<Node> &notused) : 
+              _gr(gr), _cost(cost), _used(used), _notused(notused) {}
+      
+          Node select() const {
+
+            Cost insert_val =
+              std::numeric_limits<Cost>::max();
+            int insert_node = -1;
+            
+            for (unsigned int i=0; i<_notused.size(); ++i) {
+              Cost min_val = _cost[_gr.edge(_notused[i], _used[0])];
+              int min_node = 0;
+            
+              for (unsigned int j=1; j<_used.size(); ++j) {
+                if (min_val > _cost[_gr.edge(_notused[i], _used[j])]) {
+                    min_val = _cost[_gr.edge(_notused[i], _used[j])];
+                    min_node = j;
+                }
+              }
+              
+              if (insert_val > min_val) {
+                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> &_used;
+          std::vector<Node> &_notused;
+      };
+
+      class FarthestSelection {
+        public:
+          FarthestSelection(const FullGraph &gr, const CostMap &cost,
+                            std::vector<Node> &used, std::vector<Node> &notused) : 
+              _gr(gr), _cost(cost), _used(used), _notused(notused) {}
+      
+          Node select() const {
+
+            Cost insert_val =
+              std::numeric_limits<Cost>::min();
+            int insert_node = -1;
+            
+            for (unsigned int i=0; i<_notused.size(); ++i) {
+              Cost min_val = _cost[_gr.edge(_notused[i], _used[0])];
+              int min_node = 0;
+              
+              for (unsigned int j=1; j<_used.size(); ++j) {
+                if (min_val > _cost[_gr.edge(_notused[i], _used[j])]) {
+                  min_val = _cost[_gr.edge(_notused[i], _used[j])];
+                  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> &_used;
+          std::vector<Node> &_notused;
+      };
+
+
+      class CheapestSelection {
+        private:
+          Cost costDiff(Node u, Node v, Node w) const {
+            return 
+              _cost[_gr.edge(u, w)] +
+              _cost[_gr.edge(v, w)] -
+              _cost[_gr.edge(u, v)];
+          }
+
+        public:
+          CheapestSelection(const FullGraph &gr, const CostMap &cost,
+                            std::vector<Node> &used, std::vector<Node> &notused) : 
+              _gr(gr), _cost(cost), _used(used), _notused(notused) {}
+        
+          Cost select() const {
+            int insert_notused = -1;
+            int best_insert_index = -1;
+            Cost insert_val =
+              std::numeric_limits<Cost>::max();
+            
+            for (unsigned int i=0; i<_notused.size(); ++i) {
+              int min = i;
+              int best_insert_tmp = 0;
+              Cost min_val =
+                costDiff(_used.front(), _used.back(), _notused[i]);
+              
+              for (unsigned int j=1; j<_used.size(); ++j) {
+                Cost tmp =
+                  costDiff(_used[j-1], _used[j], _notused[i]);
+
+                if (min_val > tmp) {
+                  min = i;
+                  min_val = tmp;
+                  best_insert_tmp = j;
+                }
+              }
+
+              if (insert_val > min_val) {
+                insert_notused = min;
+                insert_val = min_val;
+                best_insert_index = best_insert_tmp;
+              }
+            }
+
+            _used.insert(_used.begin()+best_insert_index, _notused[insert_notused]);
+            _notused.erase(_notused.begin()+insert_notused);
+
+            return insert_val;
+          }
+          
+        private:
+          const FullGraph &_gr;
+          const CostMap &_cost;
+          std::vector<Node> &_used;
+          std::vector<Node> &_notused;
+      };
+
+      class RandomSelection {
+        public:
+          RandomSelection(const FullGraph &, const CostMap &,
+                            std::vector<Node> &, std::vector<Node> &notused) : 
+                           _notused(notused) {
+            rnd.seedFromTime();
+          }
+          
+          Node select() const {
+            const int index = rnd[_notused.size()];
+            Node n = _notused[index];
+            _notused.erase(_notused.begin()+index);
+            return n;
+          }
+        private:
+          std::vector<Node> &_notused;
+      };
+
+
+      class DefaultInsert {
+        private:
+          Cost costDiff(Node u, Node v, Node w) const {
+            return 
+              _cost[_gr.edge(u, w)] +
+              _cost[_gr.edge(v, w)] -
+              _cost[_gr.edge(u, v)];
+          }
+  
+        public:
+          DefaultInsert(const FullGraph &gr, const CostMap &cost,
+                        std::vector<Node> &path, Cost &fullcost) : 
+            _gr(gr), _cost(cost), _path(path), _fullcost(fullcost) {}
+          
+          void insert(Node n) const {
+            int min = 0;
+            Cost min_val =
+              costDiff(_path.front(), _path.back(), n);
+            
+            for (unsigned int i=1; i<_path.size(); ++i) {
+              Cost tmp = costDiff(_path[i-1], _path[i], n);
+              if (tmp < min_val) {
+                min = i;
+                min_val = tmp;
+              }
+            }
+            
+            _path.insert(_path.begin()+min, n);
+            _fullcost += min_val;
+          }
+        
+        private:
+          const FullGraph &_gr;
+          const CostMap &_cost;
+          std::vector<Node> &_path;
+          Cost &_fullcost;
+      };
+  
+      class CheapestInsert {
+        TEMPLATE_GRAPH_TYPEDEFS(FullGraph);
+        public:
+          CheapestInsert(const FullGraph &, const CostMap &,
+                         std::vector<Node> &, Cost &fullcost) : 
+            _fullcost(fullcost) {}
+          
+          void insert(Cost diff) const {
+            _fullcost += diff;
+          }
+
+        private:
+          Cost &_fullcost;
+      };  
+      
+      
+      template <class InitFunctor, class SelectionFunctor, class InsertFunctor>
+      void start() {
+        InitFunctor init(_gr, _cost, nodesPath, notused, sum);
+        SelectionFunctor selectNode(_gr, _cost, nodesPath, notused);
+        InsertFunctor insertNode(_gr, _cost, nodesPath, sum);
+        
+        nodesPath = init.init();
+        
+        for (int i=0; i<_gr.nodeNum()-2; ++i) {
+          insertNode.insert(selectNode.select());
+        }
+        
+        sum = _cost[ _gr.edge(nodesPath.front(), nodesPath.back()) ];
+        for (unsigned int i=0; i<nodesPath.size()-1; ++i)
+          sum += _cost[ _gr.edge(nodesPath[i], nodesPath[i+1]) ];
+      }
+    
+      const FullGraph &_gr;
+      const CostMap &_cost;
+      std::vector<Node> notused;
+      std::vector<Node> nodesPath;
+      Cost sum;
+  };
+  
+}; // namespace lemon
+
+#endif

new file lemon/nearest_neighbor_tsp.h

@@ +1 @@
+#ifndef LEMON_NEAREST_NEIGHBOUR_TSP
+#define LEMON_NEAREST_NEIGHBOUR_TSP
+
+#include <deque>
+#include <lemon/full_graph.h>
+#include <lemon/path.h>
+#include <lemon/maps.h>
+#include <lemon/tsp_utils.h>
+
+namespace lemon {
+  
+  namespace nn_helper {
+    template <class L>
+    L dequeConvert(const std::deque<FullGraph::Node> &_path) {
+      return L(_path.begin(), _path.end());
+    }
+
+    template <>
+    std::deque<FullGraph::Node> dequeConvert(const std::deque<FullGraph::Node> &_path) {
+      return _path;
+    }
+  }
+  
+  template <typename CM>
+  class NearestNeighborTsp {
+    private:
+      GRAPH_TYPEDEFS(FullGraph);
+
+    public:
+      typedef CM CostMap;
+      typedef typename CM::Value Cost;
+    
+      NearestNeighborTsp(const FullGraph &gr, const CostMap &cost) : _gr(gr), _cost(cost) {}
+
+      Cost run() {
+        _path.clear();
+
+        Edge min_edge1 = INVALID,
+             min_edge2 = INVALID;
+        
+        min_edge1 = mapMin(_gr, _cost);
+
+        FullGraph::NodeMap<bool> used(_gr, false);
+
+        Node n1 = _gr.u(min_edge1), 
+             n2 = _gr.v(min_edge1);
+        
+        _path.push_back(n1);
+        _path.push_back(n2);
+
+        used[n1] = true;
+        used[n2] = true;
+
+        min_edge1 = INVALID;
+
+        while (int(_path.size()) != _gr.nodeNum()) {
+          if (min_edge1 == INVALID) {
+            for (IncEdgeIt e(_gr, n1); e!=INVALID; ++e) {
+              if (!used[_gr.runningNode(e)]) {
+                if (min_edge1==INVALID || _cost[min_edge1] > _cost[e])
+                  min_edge1 = e;
+              }
+            }
+          }
+
+          if (min_edge2 == INVALID) {
+            for (IncEdgeIt e(_gr, n2); e!=INVALID; ++e) {
+              if (!used[_gr.runningNode(e)]) {
+                if (min_edge2==INVALID || _cost[min_edge2] > _cost[e])
+                  min_edge2 = e;
+              }
+            }
+          }
+
+          if ( _cost[min_edge1] < _cost[min_edge2] ) {
+            n1 = (_gr.u(min_edge1) == n1) ? _gr.v(min_edge1) : _gr.u(min_edge1);
+            _path.push_front(n1);
+
+            used[n1] = true;
+            min_edge1 = INVALID;
+
+            if (_gr.u(min_edge2)==n1 || _gr.v(min_edge2)==n1)
+              min_edge2 = INVALID;
+          } else {
+            n2 = (_gr.u(min_edge2) == n2) ? _gr.v(min_edge2) : _gr.u(min_edge2);
+            _path.push_back(n2);
+
+            used[n2] = true;
+            min_edge2 = INVALID;
+
+            if (_gr.u(min_edge1)==n2 || _gr.v(min_edge1)==n2)
+              min_edge1 = INVALID;
+          }
+        }
+
+        _sum = _cost[ _gr.edge(_path.back(), _path.front()) ];
+        for (unsigned int i=0; i<_path.size()-1; ++i)
+          _sum += _cost[ _gr.edge(_path[i], _path[i+1]) ];
+
+        return _sum;
+      }
+
+      
+      template <typename L>
+      void tourNodes(L &container) {
+        container(nn_helper::dequeConvert<L>(_path));
+      }
+      
+      template <template <typename> class L>
+      L<Node> tourNodes() {
+        return nn_helper::dequeConvert<L<Node> >(_path);
+      }      
+
+      const std::deque<Node>& tourNodes() {
+        return _path;
+      }
+      
+      Path<FullGraph> tour() {
+        Path<FullGraph> p;
+        if (_path.size()<2)
+          return p;
+        
+        for (unsigned int i=0; i<_path.size()-1; ++i) {
+          p.addBack(_gr.arc(_path[i], _path[i+1]));
+        }
+        p.addBack(_gr.arc(_path.back(), _path.front()));
+        
+        return p;
+      }
+      
+      Cost tourCost() {
+        return _sum;
+      }
+      
+
+  private:
+    const FullGraph &_gr;
+    const CostMap &_cost;
+    Cost _sum;
+    std::deque<Node> _path;
+  };
+
+
+}; // namespace lemon
+
+#endif

new file lemon/opt2_tsp.h

@@ +1 @@
+#ifndef LEMON_OPT2_TSP
+#define LEMON_OPT2_TSP
+
+#include <vector>
+#include <lemon/full_graph.h>
+#include <lemon/path.h>
+#include <lemon/tsp_utils.h>
+
+namespace lemon {
+  
+  namespace opt2_helper {
+    template <class L>
+    L vectorConvert(const std::vector<FullGraph::Node> &_path) {
+      return L(_path.begin(), _path.end());
+    }
+  
+    template <>
+    std::vector<FullGraph::Node> vectorConvert(const std::vector<FullGraph::Node> &_path) {
+      return _path;
+    }
+  }
+  
+  template <typename CM>
+  class Opt2Tsp {
+    private:
+      GRAPH_TYPEDEFS(FullGraph);
+
+    public:
+      typedef CM CostMap;
+      typedef typename CM::Value Cost;
+      
+    
+      Opt2Tsp(const FullGraph &gr, const CostMap &cost) : _gr(gr), _cost(cost), 
+            tmppath(_gr.nodeNum()*2) {
+            
+        for (int i=1; i<_gr.nodeNum()-1; ++i) {
+          tmppath[2*i] = i-1;
+          tmppath[2*i+1] = i+1;
+        }
+        tmppath[0] = _gr.nodeNum()-1;
+        tmppath[1] = 1;
+        tmppath[2*(_gr.nodeNum()-1)] = _gr.nodeNum()-2;
+        tmppath[2*(_gr.nodeNum()-1)+1] = 0;
+      }
+      
+      Opt2Tsp(const FullGraph &gr, const CostMap &cost, const std::vector<Node> &path) : 
+              _gr(gr), _cost(cost), tmppath(_gr.nodeNum()*2) {
+
+        for (unsigned int i=1; i<path.size()-1; ++i) {
+          tmppath[2*_gr.id(path[i])] = _gr.id(path[i-1]);
+          tmppath[2*_gr.id(path[i])+1] = _gr.id(path[i+1]);
+        }
+        
+        tmppath[2*_gr.id(path[0])] = _gr.id(path.back());
+        tmppath[2*_gr.id(path[0])+1] = _gr.id(path[1]);
+        tmppath[2*_gr.id(path.back())] = _gr.id(path[path.size()-2]);
+        tmppath[2*_gr.id(path.back())+1] = _gr.id(path.front());
+      }
+
+    private:
+      Cost c(int u, int v) {
+        return _cost[_gr.edge(_gr.nodeFromId(u), _gr.nodeFromId(v))];
+      }
+      
+      class It {
+        public:
+          It(const std::vector<int> &path, int i=0) : tmppath(path), act(i), last(tmppath[2*act]) {}
+          It(const std::vector<int> &path, int i, int l) : tmppath(path), act(i), last(l) {}
+
+          int next_index() const {
+            return (tmppath[2*act]==last)? 2*act+1 : 2*act;
+          }
+          
+          int prev_index() const {
+            return (tmppath[2*act]==last)? 2*act : 2*act+1;
+          }
+          
+          int next() const {
+            return tmppath[next_index()];
+          }
+
+          int prev() const {
+            return tmppath[prev_index()];
+          }
+          
+          It& operator++() {
+            int tmp = act;
+            act = next();
+            last = tmp;
+            return *this;
+          }
+          
+          operator int() const {
+            return act;
+          }
+          
+        private:
+          const std::vector<int> &tmppath;
+          int act;
+          int last;
+      };
+
+      bool check(std::vector<int> &path, It i, It j) {
+        if (c(i, i.next()) + c(j, j.next()) > 
+            c(i, j) + c(j.next(), i.next())) {
+
+            path[ It(path, i.next(), i).prev_index() ] = j.next();
+            path[ It(path, j.next(), j).prev_index() ] = i.next();
+
+            path[i.next_index()] = j;
+            path[j.next_index()] = i;
+
+            return true;
+        }
+        return false;
+      }
+      
+    public:
+      
+      Cost run() {
+        _path.clear();
+
+        if (_gr.nodeNum()>3) {
+
+opt2_tsp_label:
+          It i(tmppath);
+          It j(tmppath, i, i.prev());
+          ++j; ++j;
+          for (; j.next()!=0; ++j) {
+            if (check(tmppath, i, j))
+              goto opt2_tsp_label;
+          }
+          
+          for (++i; i.next()!=0; ++i) {
+            It j(tmppath, i, i.prev());
+            if (++j==0)
+              break;
+            if (++j==0)
+              break;
+            
+            for (; j!=0; ++j) {
+              if (check(tmppath, i, j))
+                goto opt2_tsp_label;
+            }
+          }
+        }
+
+        It k(tmppath);
+        _path.push_back(_gr.nodeFromId(k));
+        for (++k; k!=0; ++k)
+          _path.push_back(_gr.nodeFromId(k));
+
+        
+
+        _sum = _cost[ _gr.edge(_path.back(), _path.front()) ];
+        for (unsigned int i=0; i<_path.size()-1; ++i)
+          _sum += _cost[ _gr.edge(_path[i], _path[i+1]) ];
+        return _sum;
+      }
+
+      
+      template <typename L>
+      void tourNodes(L &container) {
+        container(opt2_helper::vectorConvert<L>(_path));
+      }
+      
+      template <template <typename> class L>
+      L<Node> tourNodes() {
+        return opt2_helper::vectorConvert<L<Node> >(_path);
+      }
+
+      const std::vector<Node>& tourNodes() {
+        return _path;
+      }
+      
+      Path<FullGraph> tour() {
+        Path<FullGraph> p;
+        if (_path.size()<2)
+          return p;
+
+        for (unsigned int i=0; i<_path.size()-1; ++i) {
+          p.addBack(_gr.arc(_path[i], _path[i+1]));
+        }
+        p.addBack(_gr.arc(_path.back(), _path.front()));
+        return p;
+      }
+      
+      Cost tourCost() {
+        return _sum;
+      }
+      
+
+  private:
+    const FullGraph &_gr;
+    const CostMap &_cost;
+    Cost _sum;
+    std::vector<int> tmppath;
+    std::vector<Node> _path;
+  };
+
+
+}; // namespace lemon
+
+#endif