#ifndef LEMON_CHRISTOFIDES_TSP_H
#define LEMON_CHRISTOFIDES_TSP_H

#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>

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