/*Найдите приближенное решение метрической неориентированной задачи коммивояжера
в полном графе (на плоскости) с помощью минимального остовного дерева,
построенного в первой задаче. Оцените качество приближения на случайном наборе
точек, нормально распределенном на плоскости с дисперсией 1. Нормально
распределенный набор точек получайте с помощью std::normal_distribution. При
фиксированном N, количестве вершин графа, несколько раз запустите оценку
качества приближения. Вычислите среднее значение и среднеквадратичное отклонение
качества приближения для данного N. Запустите данный эксперимент для всех N в
некотором диапазоне, например, [2, 10]. Автоматизируйте запуск экспериментов. В
решении требуется разумно разделить код на файлы. Каждому классу - свой
заголовочный файл и файл с реализацией.
*/

#include "Graph.h"
#include <math.h>
#include <iostream>
#include <limits>
#include <memory>
#include <optional>
#include <random>
#include <set>
#include <stack>
#include <vector>

using IntPair = std::pair<int, int>;
using IntVector = std::vector<int>;

bool operator<(Edge edge1, Edge edge2) {
  auto edge1_pair = std::make_pair(edge1.vertex1_, edge1.vertex2_);
  auto edge2_pair = std::make_pair(edge2.vertex1_, edge2.vertex2_);
  return edge1_pair < edge2_pair;
}

Graph::Graph(int vertex_count)
    : adjacency_matrix_(vertex_count, std::vector<double>(vertex_count, 0)) {}

void Graph::AddEdge(int vertex1, int vertex2, double weight) {
  adjacency_matrix_[vertex1][vertex2] = weight;
  adjacency_matrix_[vertex2][vertex1] = weight;
  edges_.emplace(vertex1, vertex2, weight);
  edges_.emplace(vertex2, vertex1, weight);
}

int Graph::GetVertexCount() const { return adjacency_matrix_.size(); }

void Graph::FindComponentVertices(std::vector<int>& component_indexes,
                                  int component_representative,
                                  std::vector<bool>& is_visited,
                                  int component_index) const {
  component_indexes[component_representative] = component_index;
  is_visited[component_representative] = true;
  for (int vertex = 0; vertex < GetVertexCount(); ++vertex) {
    if (!is_visited[vertex] &&
        adjacency_matrix_[component_representative][vertex] > 0) {
      FindComponentVertices(component_indexes, vertex, is_visited,
                            component_index);
    }
  }
}

Graph Graph::GetForest() const {
  auto forest = std::make_shared<Graph>(GetVertexCount());
  return *forest;
}

inline bool Graph::IsUpdateNeeded(
    int vertex1_component_index,
    const std::vector<std::optional<IntPair>>& new_edges,
    double new_weight) const {
  return !new_edges[vertex1_component_index].has_value() ||
         adjacency_matrix_[new_edges[vertex1_component_index].value().first]
                          [new_edges[vertex1_component_index].value().second] >
             new_weight;
}

inline bool Graph::IsAddedEdge(const std::set<IntPair>& added_edges,
                               const std::optional<IntPair> new_edge) const {
  return added_edges.find(new_edge.value()) == added_edges.end();
}

void Graph::AddNewEdgesToForest(
    Graph& forest, int& added_edge_count, std::set<IntPair>& added_edges,
    const std::vector<std::optional<IntPair>>& new_edges,
    double& minimal_spanning_tree_weight) const {
  for (auto&& new_edge : new_edges) {
    if (added_edge_count == GetVertexCount() - 1) {
      break;
    }
    if (IsAddedEdge(added_edges, new_edge)) {
      minimal_spanning_tree_weight +=
          adjacency_matrix_[new_edge.value().first][new_edge.value().second];
      added_edge_count++;
      forest.AddEdge(
          new_edge.value().first, new_edge.value().second,
          adjacency_matrix_[new_edge.value().first][new_edge.value().second]);
      added_edges.emplace(new_edge.value());
      added_edges.emplace(new_edge.value().second, new_edge.value().first);
    }
  }
}

Graph Graph::GetMinimalSpanningTree() {
  std::set<IntPair> added_edges;
  double minimal_spanning_tree_weight = 0;
  auto forest = GetForest();
  int added_edge_count = 0;

  while (added_edge_count < GetVertexCount() - 1) {
    std::vector<int> component_indexes(GetVertexCount());
    std::vector<std::optional<IntPair>> new_edges(GetVertexCount() -
                                                  added_edge_count);
    int component_index = 0;
    std::vector<bool> is_visited(GetVertexCount(), false);

    for (int vertex = 0; vertex < GetVertexCount(); ++vertex) {
      if (!is_visited[vertex]) {
        forest.FindComponentVertices(component_indexes, vertex, is_visited,
                                     component_index);
        component_index++;
      }
    }
    for (int vertex1 = 0; vertex1 < GetVertexCount(); ++vertex1) {
      for (int vertex2 = 0; vertex2 < adjacency_matrix_[vertex1].size();
           ++vertex2) {
        if (component_indexes[vertex1] != component_indexes[vertex2]) {
          int vertex1_component_index = component_indexes[vertex1];
          if (IsUpdateNeeded(vertex1_component_index, new_edges,
                             adjacency_matrix_[vertex1][vertex2])) {
            new_edges[vertex1_component_index].emplace(vertex1, vertex2);
          }
        }
      }
    }
    AddNewEdgesToForest(forest, added_edge_count, added_edges, new_edges,
                        minimal_spanning_tree_weight);
  }
  return forest;
}

std::vector<int> Graph::GetEylerianCycle() {
  std::set<Edge> erased_edges;
  std::vector<int> eylerian_cycle;
  std::stack<int> added_vertices;
  added_vertices.emplace(0);

  while (!added_vertices.empty()) {
    auto current_vertex = added_vertices.top();
    for (auto&& edge : edges_) {
      if (edge.vertex1_ == current_vertex) {
        added_vertices.emplace(edge.vertex2_);
        if (erased_edges.find(edge) == erased_edges.end()) {
          erased_edges.emplace(edge);
          erased_edges.emplace(edge.vertex2_, edge.vertex1_, edge.weight_);
          break;
        }
        edges_.erase(edge);
        auto edge_copy =
            std::make_shared<Edge>(edge.vertex2_, edge.vertex1_, edge.weight_);
        edges_.erase(*edge_copy);
        break;
      }
    }
    if (current_vertex == added_vertices.top()) {
      added_vertices.pop();
      eylerian_cycle.emplace_back(current_vertex);
    }
  }
  return eylerian_cycle;
}

double Graph::GetApproximateCycleWeight(IntVector eylerian_cycle) const {
  double tour_cost = 0;
  int prev_vertex = eylerian_cycle[0];
  std::vector<bool> is_visited(GetVertexCount(), false);
  is_visited[eylerian_cycle[0]] = true;
  for (auto&& vertex : eylerian_cycle) {
    if (!is_visited[vertex]) {
      is_visited[vertex] = true;
      tour_cost += adjacency_matrix_[prev_vertex][vertex];
      prev_vertex = vertex;
    }
  }
  tour_cost += adjacency_matrix_[prev_vertex][eylerian_cycle[0]];
  return tour_cost;
}

int GetTwoPower(int power) {
  int two_in_power = 1;
  for (int i = 0; i < power; ++i) {
    two_in_power *= 2;
  }
  return two_in_power;
}

int GetBit(int bit_index, int mask) {
  for (int i = 0; i < bit_index; ++i) {
    mask /= 2;
  }
  return mask % 2;
}

double Graph::GetGamiltonicPathWeight(
    int mask, int to,
    std::vector<std::vector<std::optional<double>>>& gamiltonic_path_weights)
    const {
  double gamilton_path_weight = std::numeric_limits<double>::max();
  if (mask == 1 && to == 0) {
    return 0;
  }
  if (GetBit(0, mask) == 1 && GetBit(to, mask) == 1 && to > 0) {
    for (int vertex = 0; vertex < GetVertexCount(); ++vertex) {
      if (adjacency_matrix_[to][vertex] > 0) {
        if (GetBit(vertex, mask) == 1) {
          if (!gamiltonic_path_weights[mask ^ GetTwoPower(to)][vertex]
                   .has_value()) {
            gamiltonic_path_weights[mask ^ GetTwoPower(to)][vertex].emplace(
                GetGamiltonicPathWeight(mask ^ GetTwoPower(to), vertex,
                                        gamiltonic_path_weights));
          }
          if (gamiltonic_path_weights[mask ^ GetTwoPower(to)][vertex].value() +
                  adjacency_matrix_[to][vertex] <
              gamilton_path_weight) {
            gamilton_path_weight =
                gamiltonic_path_weights[mask ^ GetTwoPower(to)][vertex]
                    .value() +
                adjacency_matrix_[to][vertex];
          }
        }
      }
    }
    return gamilton_path_weight;
  }
  return std::numeric_limits<double>::max();
}

double Graph::GetGamiltonicCycleWeight() const {
  std::vector<std::vector<std::optional<double>>> gamiltonic_way_weights(
      GetTwoPower(GetVertexCount()),
      std::vector<std::optional<double>>(GetVertexCount()));
  double gamiltonic_cycle_weight = std::numeric_limits<double>::max();

  for (int vertex = 1; vertex < GetVertexCount(); ++vertex) {
    if (adjacency_matrix_[0][vertex] > 0) {
      if (!gamiltonic_way_weights[GetTwoPower(GetVertexCount()) - 1][vertex]
               .has_value()) {
        gamiltonic_way_weights[GetTwoPower(GetVertexCount()) - 1][vertex]
            .emplace(GetGamiltonicPathWeight(GetTwoPower(GetVertexCount()) - 1,
                                             vertex, gamiltonic_way_weights));
      }
      if (gamiltonic_way_weights[GetTwoPower(GetVertexCount()) - 1][vertex]
                  .value() +
              adjacency_matrix_[0][vertex] <
          gamiltonic_cycle_weight) {
        gamiltonic_cycle_weight =
            gamiltonic_way_weights[GetTwoPower(GetVertexCount()) - 1][vertex]
                .value() +
            adjacency_matrix_[0][vertex];
      }
    }
  }
  return gamiltonic_cycle_weight;
}

double GetStandardDeviation(double average_value,
                            const std::vector<double>& values) {
  double standard_deviation = 0;
  for (auto&& value : values) {
    standard_deviation +=
        (average_value - value) * (average_value - value) / (values.size());
  }
  return sqrt(standard_deviation);
}

void CarryOutExperiment() {
  const int experiment_count = 10;
  const int maximal_point_count = 10;
  std::default_random_engine generator;
  std::normal_distribution<double> distributor{0, 1};
  for (int i = 2; i < maximal_point_count + 1; ++i) {
    double coefficient_sum = 0;
    std::vector<double> experiment_results(experiment_count);
    for (int j = 0; j < experiment_count; ++j) {
      Graph graph = Graph(i);
      std::vector<std::pair<double, double>> points;
      points.reserve(i);
      for (int point_index = 0; point_index < i; ++point_index) {
        double x = distributor(generator);
        double y = distributor(generator);
        points.emplace_back(x, y);
        for (int k = 0; k < points.size(); ++k) {
          double distance =
              sqrt((x - points[k].first) * (x - points[k].first) +
                   (y - points[k].second) * (y - points[k].second));
          graph.AddEdge(k, point_index, distance);
        }
      }
      Graph minimal_spanning_tree = graph.GetMinimalSpanningTree();
      IntVector eylerian_cycle = minimal_spanning_tree.GetEylerianCycle();
      double approximate_cycle_weight =
          graph.GetApproximateCycleWeight(eylerian_cycle);
      double gamiltonic_cycle_weight = graph.GetGamiltonicCycleWeight();
      coefficient_sum += approximate_cycle_weight / gamiltonic_cycle_weight;
      experiment_results[j] =
          approximate_cycle_weight / gamiltonic_cycle_weight;
    }
    std::cout << "Average approximation coefficient is "
              << coefficient_sum / experiment_count << std::endl;
    std::cout << "Standard deviation is "
              << GetStandardDeviation(coefficient_sum / experiment_count,
                                      experiment_results)
              << std::endl;
  }
}

int main() {
  CarryOutExperiment();
  return 0;
}