/**
 * @file EulerHamiltonianLoops.cpp
 * @brief finding Euler and Hamiltonian loops in undirected graphs
 * @author Dmitri Kuvshinov
 */

#include <stack>
#include <vector>
#include <limits>
#include <algorithm>

/*
Получение некоторого эйлерова цикла.
*/

/// algorithm finds some Euler loop in a connected undirected graph
/**
 * @param graph adjacencies vector ("AdjVec") of the input graph, all degrees must be even
 * @param out output iterator for writing out found Euler loop
 */
template<class AdjVec, class OutIt>
void findEulerLoop(const AdjVec &graph;, OutIt out)
{
  // будем удалять рёбра, поэтому создадим локальную копию
  AdjVec G = graph;
  // вспомогательный стек
  std::stack<unsigned> S;
  // начинает с 0-й вершины (можно с любой)
  S.push(0); 
  while (!S.empty())
  {
    // смотрим вершину на стеке
    const unsigned v = S.top();
    if (!G[v].empty()) // у неё есть соседи?
    {
      // выберем первую же вершину из списка смежности v
      const unsigned u = *(G[v].begin()); 
      S.push(u); // положим её в стек
      // удалим ребро u-v
      G[v].erase(u);
      G[u].erase(v);
    }
    else // соседей (уже) нет
    { 
      // удалим из стека и запишем в результат
      S.pop();
      *out++ = v;
    }
  }
}


/*
Инструмент, позволяющий перебирать все гамильтоновы циклы.
Вызов next() генерирует следующий цикл. Если такой был найден,
  то индексы вершин, составляющих цикл, задаются диапазоном
  итераторов begin(), end().
AdjVec -- граф, представленный как вектор наборов соседей.
Является основой для MinimalHamiltonianLoop.
*/

/// the object provides facility to enumerate all Hamiltonian loops in a graph
template<class AdjVec>
class HamiltonianLoopEnumerator
{
  const AdjVec &G; // ссылка на исходный граф
  unsigned N;      // количество вершин
  
  AdjVec spare;    // неиспользованные соседи (ещё не включенные в цепь)
  std::vector<unsigned> chain; // текущая цепь вершин
  std::vector<bool> status;    // статус "вершина использована"

  typedef typename AdjVec::value_type Adjacency;    // набор соседей
  typedef typename Adjacency::iterator AdjacencyIt; // итератор набора соседей
  std::vector<AdjacencyIt> cur; // текущие позиции в spare для устранения рекурсии
  unsigned k;      // количество вершин в цепи

  // откат к k - 1, "убирает" 1 или 2 вершины
  void rollback()
  {
    status[chain[k]] = false;
    ++cur[--k];
    status[chain[k]] = false;
  }

  // подготовка структур данных для k + 1, заодно увеличивает k
  // кусок кода вынесен для упрощения восприятия кода next()
  void nextK()
  {
    const unsigned y = *(cur[k]); // выберем следующего соседа
    status[y] = true;  // пометим как использованного
    chain[k] = y;  // поместим его в цепь

    // построим набор доступных вершин для выбора следующей вершины в цепи
    Adjacency &nextSpare; = spare[++k]; // ! увеличили k на 1
    nextSpare.clear();
    for (typename AdjVec::value_type::const_iterator
          q = G[y].begin(), qend = G[y].end(); q != qend; ++q)
    {
      if (!status[*q]) // выбираем только среди ещё неиспользованных
        nextSpare.insert(*q);
    }

    cur[k] = spare[k].begin();
  }

public:

  // подготовить объект к поиску гамильтоновых циклов
  /// prepare internals according to the graph given
  /**
   * graph must not be destroyed before any next() calls
   */
  HamiltonianLoopEnumerator(const AdjVec &graph;)
    : G(graph)
    , N(graph.size())
    , spare(N)
    , chain(N + 1, 0)
    , status(N, false)
    , cur(N)
    , k(1)
  {
    status[0] = true; // 0-я вершина всегда включена и считается начальной
    // всегда chain[0] == chain[N] == 0 -- в цикле N+1 вершина
    spare[1] = graph[0];
    cur[1] = spare[1].begin();
  }

  // начало построенного с помощью next() гамильтонова цикла
  /// get the beginning of the current Hamiltonian loop built by next()
  std::vector<unsigned>::const_iterator
    begin() const
  {
    return chain.begin();
  }

  // конец построенного с помощью next() гамильтонова цикла
  /// get the end of the current Hamiltonian loop built by next()
  std::vector<unsigned>::const_iterator
    end() const
  {
    return chain.end();
  }

  // собственно алгоритм (рекурсивный перебор с возвратом)
  /// build next Hamiltonian loop
  /**
   * @return success indicator, false means no more loops are available
   */
  bool next()
  {
    if (k < N - 1) // достроим цепь
    {
      for (;;)
      {
        while (cur[k] == spare[k].end() && k > 1) // перебрали всех соседей
          rollback(); // откат

        if (cur[k] == spare[k].end() && k == 1)  // откатились "до упора"?
          return false;

        nextK(); // перейдём к следующему k

        // повторим для следующей вершины в цепи
        const bool found = next();

        if (found) // есть цикл -- вернём true
          return true; // нет -- следующая итерация for

        status[chain[k]] = false; // откат, k уменьшил неудачный вызов next()
        // следующая итерация/вызов начнёт со следующей вершины из spare[k] (или откатится, если вершины кончились)
        ++cur[k];  
      }
    }
    else // k == N - 1, попытаемся замкнуть цепь в цикл, осталось не более одной вершины
    {
      const AdjacencyIt cend = spare[k].end();
      AdjacencyIt &cur;_k = cur[k];
      for (;; ++cur_k)
      {
        if (cur_k == cend)
        {
          --k;
          return false; // замкнуть цепь не удалось
        }

        // если можно перейти из *(cur[k]) в 0, остановимся на ней
        if (G[*cur_k].find(0) != G[*cur_k].end()) // 0-я есть среди соседей
        {
          chain[k] = *cur_k;
          rollback();
          return true; // замкнули цепь
        }
      }
    }
  }
};



/*
Поиск минимального гамильтонова цикла перебором (O(N!) операций).
Надстройка над HamiltonianLoopEnumerator.
По сути нужен только для проверки правильности реализации продвинутых методов.
*/

/// minimal Hamiltonian loop finder, use run() then begin(), end() to get the loop
template<class AdjVec>
class MinimalHamiltonianLoop
{
  HamiltonianLoopEnumerator<AdjVec> hle; // средство перебора циклов
  double curw; // current weight, вес текущего цикла
  std::vector<unsigned> loop; // текущий цикл

  // сравнить loop с текущим циклом hle, заменить, если он легче
  template<class WeightMatrix>
  void checkMinimal(const WeightMatrix &W)
  {
    double w = 0.0;
    for (std::vector<unsigned>::const_iterator
      p = hle.begin(), pend = hle.end() - 1; p != pend; ++p)
    {
      w += W(*p, *(p + 1));
    }

    if (w < curw)
    {
      curw = w;
      std::copy(hle.begin(), hle.end(), loop.begin());
    }
  }

public:

  /// graph must survive until run() call
  MinimalHamiltonianLoop(const AdjVec &graph;)
    : hle(graph)
    , curw(std::numeric_limits<double>::infinity())
    , loop(graph.size() + 1)
  {
  }

  /// get the beginning of the minimal loop after run() call
  std::vector<unsigned>::const_iterator
    begin() const
  {
    return loop.begin();
  }

  /// get the end of the minimal loop after run() call
  std::vector<unsigned>::const_iterator
    end() const
  {
    return loop.end();
  }

  /// weight of the current Hamiltonian loop given by begin(), end()
  double weight() const
  {
    return curw;
  }

  // перебрать все циклы, сохранив цикл минимального веса
  /// find out the minimal Hamiltonian loop according to weight matrix W
  template<class WeightMatrix>
  void run(const WeightMatrix &W)
  {
    while (hle.next())
      checkMinimal(W);
  }

  // перебрать не более maxloops циклов, сохранив цикл минимального среди них веса
  /// find out the minimal Hamiltonian loop amidst some no more than maxloops loops according to weight matrix W
  template<class WeightMatrix>
  void run(const WeightMatrix &W, std::size_t maxloops)
  {
    while (maxloops-- != 0 && hle.next())
      checkMinimal(W);
  }
};

///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <iterator>
#include <set>
#include <ctime>

#include "Graph2.h"


// тесты
class TestCase3
{
  static const unsigned N = 9; // количество вершин
  NaiveMatrix<> weights;

  // представление графов, подходящее в качестве AdjVec
  typedef std::vector< std::set<unsigned> > MyGraph;
  MyGraph graph;

public:

  // создадим случайную матрицу весов
  TestCase3()
    : weights(N, N, std::numeric_limits<double>::infinity())
    , graph(N)
  {
    // случайная матрица для полного графа
    for (unsigned i = 0; i < N; ++i)
    {
      for (unsigned j = i + 1; j < N; ++j)
      {
        weights(i, j) = std::rand();
        weights(j, i) = weights(i, j);

        graph[i].insert(j);
        graph[j].insert(i);
      }
    }

    if (N <= 10)
      weights.print();
    std::cout << "\n";
  }


  // тест 1 -- эйлеров цикл
  void testEulerLoop()
  {
    findEulerLoop(graph, std::ostream_iterator<unsigned>(std::cout, ", "));
    std::cout << "\n";
  }

  // тест 2 -- гамильтоновы циклы
  void testHLE()
  {
    HamiltonianLoopEnumerator<MyGraph> hle(graph);
    for (int i = 0; i < 10; ++i)
    {
      hle.next();
      std::copy(hle.begin(), hle.end(),
        std::ostream_iterator<unsigned>(std::cout, ", "));
      std::cout << "\n";
    }
  }

  // тест 3 -- минимальный гамильтонов цикл
  void testMHL()
  {
    MinimalHamiltonianLoop<MyGraph> mhl(graph);
    mhl.run(weights);
    std::copy(mhl.begin(), mhl.end(),
      std::ostream_iterator<unsigned>(std::cout, ", "));
    std::cout << "\t" << mhl.weight() << "\n";
  }

  // тест 4 -- заданный минимальный гамильтонов цикл
  void testMHL2()
  {
    unsigned loop[N + 1] =
    { 0, 2, 1, 4, 3, 6, 5, 8, 7, 0 };

    for (unsigned i = 0; i < N; ++i)
    {
      weights(loop[i], loop[i + 1]) = 0;
      weights(loop[i + 1], loop[i]) = 0;
    }

    std::cout << "\n";
    std::copy(loop, loop + N + 1,
      std::ostream_iterator<unsigned>(std::cout, ", "));
    std::cout << "\n";
    testMHL();
  }
};


int main()
{
  using namespace std;
  srand(time(NULL));
  TestCase3 tc;
  tc.testEulerLoop();
  tc.testHLE();
  tc.testMHL();
  tc.testMHL2();

  cin.get();
}