#include <iostream>
#include <cmath>
#include <functional>
#include <random>
#include <vector>
#include <limits>

#include "randutils/randutils.hpp"
#include "pcg-cpp/include/pcg_random.hpp"

using Rng = randutils::random_generator<pcg32>;

class HalfEdge;

class Vertex {
public:
  unsigned index;
  std::vector<std::reference_wrapper<HalfEdge>> neighbors;
  std::array<unsigned, 2> coordinate;

  void addEdge(HalfEdge& e) {
    neighbors.push_back(e);
  }
};

class Edge;

class HalfEdge {
private:
  Vertex* tail;
  Vertex* head;

public:
  const Edge& edge;
  double oldX;
  double X;

  HalfEdge(const Edge& e) : edge(e) {
    X = 0;
    oldX = X;
  }
  void setVertices(Vertex& vt, Vertex& vh) {
    tail = &vt;
    head = &vh;
  }
  const Vertex& getHead() const {
    return *head;
  }
  const Vertex& getTail() const {
    return *tail;
  }
};

class Edge {
public:
  std::array<HalfEdge, 2> halfedges;
  double weight;

  Edge() : halfedges{*this, *this} {};
  void setVertices(Vertex& red, Vertex& blue) {
    halfedges[0].setVertices(red, blue);
    halfedges[1].setVertices(blue, red);
    red.addEdge(halfedges[0]);
    blue.addEdge(halfedges[1]);
  }
  double belief() const {
    return halfedges[0].X + halfedges[1].X;
  }
  bool active() const {
    return belief() >= 0;
  }
};

class Graph {
public:
  std::vector<Vertex> vertices;
  std::vector<Edge> edges;

  Graph(unsigned n, Rng& r) : vertices(2 * n * (n + 1)), edges(pow(2 * n, 2)) {
    unsigned M = vertices.size() / 2;
    for (unsigned i = 0; i < M; i++) {
      vertices[i].index = i;
      vertices[i].coordinate = {2 * (i % (n + 1)), 2 * (i / (n + 1)) + 1};
      vertices[M + i].index = M + i;
      vertices[M + i].coordinate = {2 * (i % n) + 1, 2 * (i / n)};
    }
    for (unsigned i = 0; i < edges.size(); i++) {
      Vertex& redVertex = vertices[(1 + (i % (2 * n))) / 2 + (n + 1) * ((i / 4) / n)];
      Vertex& blueVertex = vertices[M + (i % (2 * n)) / 2 + n * (((i + 2 * n) / 4) / n)];
      edges[i].setVertices(redVertex, blueVertex);
      edges[i].weight = r.variate<double, std::exponential_distribution>(1);
    }
  }

  void propagateBeliefs() {
    for (Edge& e : edges) {
      for (HalfEdge& h : e.halfedges) {
        h.X = std::numeric_limits<double>::infinity();
        for (const HalfEdge& hn : h.getHead().neighbors) {
          if (h.getTail().index != hn.getHead().index) {
            h.X = std::min(hn.edge.weight - hn.oldX, h.X);
          }
        }
      }
    }

    for (Edge& e : edges) {
      for (HalfEdge& he : e.halfedges) {
        he.oldX = he.X;
      }
    }
  }
};

int main() {
  unsigned n = 100;
  unsigned steps = 1e5;
  Rng r;
  Graph G(n, r);

  for (unsigned i = 0; i < steps; i++) {
    G.propagateBeliefs();
  }

  for (unsigned i = 0; i < G.edges.size(); i++) {
    if (G.edges[i].active()) {
      std::cout
        << G.edges[i].halfedges[0].getTail().coordinate[0] << " "
        << G.edges[i].halfedges[0].getTail().coordinate[1] << " "
        << G.edges[i].halfedges[0].getHead().coordinate[0] << " "
        << G.edges[i].halfedges[0].getHead().coordinate[1] << std::endl;
    }
  }

  return 0;
}