/**
 * @file WeightMatrix.cpp
 * @brief minimization algorithms on graphs represented by weight matrices
 * @author Dmitri Kuvshinov
 * @version rev.2
 */

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

#include <ctime>
#include <cstdlib>
#include <iostream>

/*
Построение остова минимального веса алгоритмом Ярника-Прима-Дейкстры.
Вход: N -- количество вершин в графе, W -- матрица весов рёбер.
Выход: Pr -- массив предшественников: остов содержит рёбра вида (Pr[i], i), D -- хранит веса этих рёбер.
Непосредственно функция возвращает вес построенного остова (сумма элементов D).
*/

/// Jarnik-Prim-Dijkstra
template<class Previous, class WeightVector, class WeightMatrix>
double minimalSpanningTree
  (Previous &Pr;, WeightVector &D, const WeightMatrix &W, unsigned N)
{
  if (N == 0) return 0.0; // нет вершин -- нет остова

  // признак "неиспользованности" вершины
  std::vector<bool> spare(N, true);

  // начинаем с первой вершины (индекс 0)
  spare[0] = false; 
  for (unsigned v = 0; v < N; ++v)
  {
    Pr[v] = 0;
    D[v] = W(v, 0);
  }

  // N - 1 раунд минимизации, каждый раунд из 
  // неиспользованных вершин выбирается одна ближайшая
  for (unsigned round = 1; round < N; ++round)
  {
    // найдём ближайшую простым перебором
    unsigned v = 0;
    double minD = std::numeric_limits<double>::infinity();
    for (unsigned u = 0; u < N; ++u)
      if (spare[u] && D[u] < minD)
      {
        v = u;
        minD = D[u];
      }

    // пометим найденную
    // ребро (Pr[v], v) входит в остов
    spare[v] = false;

    // обновим D и Pr для неиспользованных вершин
    // относительно найденной вершины v
    for (unsigned u = 0; u < N; ++u)
    {
      if (spare[u] && D[u] > W(u, v))
      {
        Pr[u] = v;
        D[u] = W(u, v);
      }
    }
  }

  // посчитаем сумму весов рёбер остова
  double weight = 0.0;
  for (unsigned i = 0; i < N; ++i)
    weight += D[i];
  return weight;
}


/*
Алгоритм Дейкстры поиска кратчайших путей их вершины S до всех остальных вершин графа.
Требуются неотрицательные веса рёбер матрицы W.
Вход: N -- количество вершин в графе, S -- начальная вершина, W -- матрица весов рёбер.
Выход: Pr -- массив предшественников: кратчайший путь из S в i приходит из вершины Pr[i],
  D -- массив длин путей, D[i] -- длина пути из S в i.
*/

/// Dijkstra's algorithm: distances from s to all other
template<class Previous, class WeightVector, class WeightMatrix>
void dijkstra(
    Previous &Pr;, 
    WeightVector &D,
    const WeightMatrix &W,
    unsigned S,
    unsigned N)
{
  // признак "непройденности" вершин
  std::vector<bool> spare(N, true);
  spare[S] = false;

  // подготовка
  for (unsigned i = 0; i < N; ++i)
  {
    D[i] = W(S, i);
    Pr[i] = S;
  }
 
  D[S] = 0;
  

  // раунды минимизации, каждый раунд
  // выбирается вершина, ближайшая из непройденных
  for (unsigned round = 1; round < N; ++round)
  {
    // найдём ближайшую простым перебором
    unsigned v = 0;
    double minD = std::numeric_limits<double>::infinity();
    for (unsigned u = 0; u < N; ++u)
      if (spare[u] && D[u] < minD)
      {
        v = u;
        minD = D[u];
      }

    // пометим её как пройденную
    spare[v] = false;
    
    // обновим D и Pr для непройденных вершин
    // относительно найденной вершины v
    for (unsigned u = 0; u < N; ++u)
      if (spare[u])
      {
        const double weight = D[v] + W(v, u);
        if (weight < D[u])
        {
          D[u] = weight;
          Pr[u] = v;
        }
      }
  }
}


/*
Алгоритм Флойда-Уоршелла поиска кратчайших путей между всеми парами вершин.
Вход: N -- количество вершин в графе, W -- матрица весов рёбер.
Выход: Pr -- матрица предшественников: на пути из i в j вершина Pr(i, j) предшествует j,
  D -- матрица длин путей: D(i, j) -- длина кратчайшего пути из i в j.
*/

/// Floyd–Warshall algorithm: distances between all pairs of vertices
template<class PreviousMatrix, class WeightMatrix>
void floyd(
    PreviousMatrix &Pr;, 
    WeightMatrix &D, 
    const WeightMatrix &W,
    unsigned N)
{
  // подготовка
  D = W;
  for (unsigned i = 0; i < N; ++i)
    for (unsigned j = 0; j < N; ++j)
      Pr(i, j) = i;

  // минимизация
  for (unsigned k = 0; k < N; ++k)
    for (unsigned i = 0; i < N; ++i)
    {
      const double D_ik = D(i, k);
      for (unsigned j = 0; j < N; ++j)
      {
        const double weight = D_ik + D(k, j);
        if (weight < D(i, j))
        {
          D(i, j) = weight;
          Pr(i, j) = Pr(k, j);
        }
      }
    }
}


/*
Простейшая "наивная" реализация матрицы как вектора векторов.
*/

/// simplest matrix implemenation over std::vector
template<class T = double>
class NaiveMatrix
{
  std::vector< std::vector<T> > data;

public:
  // создать пустую матрицу
  NaiveMatrix() {}

  // создать матрицу N на M, элементы которой будут равны value
  NaiveMatrix(unsigned N, unsigned M, T value = T(0))
    : data(N)
  {
    for (unsigned i = 0; i < M; ++i)
      data[i].resize(M, value);
  }

  // доступ к элементу i, j
  T& operator()(unsigned i, unsigned j)
  {
    return data[i][j];
  }

  // константный доступ к элементу i, j
  T operator()(unsigned i, unsigned j) const
  {
    return data[i][j];
  }

  // размер матрицы по первой координате
  unsigned size1() const
  {
    return data.size();
  }

  // размер матрицы по второй координате
  unsigned size2() const
  {
    return data.front().size();
  }

  // вывод матрицы в cout
  void print() const
  {
    for (unsigned i = 0; i < data.size(); ++i)
    {
      const std::vector<T> &row; = data[i];
      for (unsigned j = 0; j < row.size(); ++j)
        std::cout << row[j] << '\t';
      std::cout << '\n';
    }
  }
};



/*
Случайное заполнение матрицы.
Вход: N -- количество вершин, k -- количество исходящих рёбер из каждой вершины.
Выход: W -- матрица весов (надо предварительно заполнить "бесконечностями").
*/

/// make random k-halfdegree graph, N <= RAND_MAX, k < N
template<class WeightMatrix>
void randomMatrix(WeightMatrix &W, unsigned k, unsigned N)
{
  std::vector<bool> pos(N - 1, false);
  std::fill(pos.begin(), pos.begin() + k, true);

  for (unsigned i = 0; i < N; ++i)
  {
    std::random_shuffle(pos.begin(), pos.end());
    for (unsigned j = 0; j < i; ++j)
      if (pos[j])
        W(i, j) = std::rand();

    W(i, i) = 0;

    for (unsigned j = i + 1; j < N; ++j)
      if (pos[j - 1])
        W(i, j) = std::rand();
  }
}


/*
Простая реализация таймера для замера промежутков времени.
*/

/// simple timer implementation
class Timer
{
  clock_t t;
public:
  Timer() : t(clock()) {}
  clock_t operator()()
  {
    clock_t d = clock() - t;
    t += d;
    return d;
  }
};


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

public:

  // создадим случайную матрицу весов
  TestCase()
    : weights(N, N, std::numeric_limits<double>::infinity())
  {
    randomMatrix(weights, N / 4, N);
    if (N <= 10)
      weights.print();
  }


  // тест 1
  void testSpanningTree()
  {
    using namespace std;
    // минимальный остов
    vector<unsigned> prev(N, 0);
    vector<double> dist(N, numeric_limits<double>::infinity());
    cout << minimalSpanningTree(prev, dist, weights, N) << endl;

    // замерить скорость алгоритма построения минимального остова
    Timer timer;
    for (int i = N / 6; i; --i) // N/6 повторов
    {
      // очистить prev и dist
      fill(prev.begin(), prev.end(), 0);
      fill(dist.begin(), dist.end(), numeric_limits<double>::infinity());

      minimalSpanningTree(prev, dist, weights, N);
    }

    cout << "spanning tree time = " << timer() << endl;
  }

  
  // тест 2
  void testDijkstra()
  {
    using namespace std;
    // найти кратчайшие пути из вершины 0
    vector<unsigned> prev(N, 0);
    vector<double> dist(N, std::numeric_limits<double>::infinity());
    dijkstra(prev, dist, weights, 0, N);
    cout << dist[N / 2] << endl;

    // замерить скорость алгоритма Дейкстры
    Timer timer;
    for (int i = N / 2; i; --i)  // N/2 повторов
    {
      // очистить prev и dist
      fill(prev.begin(), prev.end(), 0);
      fill(dist.begin(), dist.end(), numeric_limits<double>::infinity());

      dijkstra(prev, dist, weights, i - 1, N);
    }

    cout << "dijkstra time = " << timer() << endl;
  }


  // тест 3
  void testFloyd()
  {
    using namespace std;
    // найти кратчайшие пути между всеми парами вершин
    NaiveMatrix<unsigned> prevmat(N, N);
    NaiveMatrix<> distmat(N, N);
    // одновременно замерим скорость одного прогона алгоритма Флойда-Уоршелла
    Timer timer;

    floyd(prevmat, distmat, weights, N);
    cout << distmat(0, N / 2);
    cout << "\nfloyd time = " << timer() << endl;
  }
};


int main()
{
  using namespace std;
  
  srand(time(NULL));
  TestCase tc;
  tc.testSpanningTree();
  tc.testDijkstra();
  tc.testFloyd();

  cin.get();
  return 0;
}