#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 {
private:
  unsigned index;

public:
  std::vector<HalfEdge*> neighbors;

  void setIndex(unsigned i) {index = i;}
  unsigned getIndex() const {return index;}
  void addEdge(HalfEdge& e) {
    neighbors.push_back(&e);
  }
};

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

public:
  double oldX;
  double X;
  double w;
  void setIndex(unsigned i) {index = i;}
  unsigned getIndex() const {return index;}
  void setVertices(Vertex& vt, Vertex& vh) {
    tail = &vt;
    head = &vh;
  }
  Vertex* getHead() {
    return head;
  }
  Vertex* getTail() {
    return tail;
  }
  const Vertex* getHead() const {
    return head;
  }
  const Vertex* getTail() const {
    return tail;
  }
};

class BipartiteGraph {
public:
  std::vector<Vertex> redVertices;
  std::vector<Vertex> blueVertices;
  std::vector<HalfEdge> redEdges;
  std::vector<HalfEdge> blueEdges;

  BipartiteGraph(unsigned n, Rng& r) : redVertices(n * (n + 1)), blueVertices(n * (n + 1)), redEdges(pow(2 * n, 2)), blueEdges(pow(2 * n, 2)) {
    for (unsigned i = 0; i < redVertices.size(); i++) {
      redVertices[i].setIndex(i);
      blueVertices[i].setIndex(i);
    }
    for (unsigned i = 0; i < redEdges.size(); i++) {
      redEdges[i].setIndex(i);
      blueEdges[i].setIndex(i);
      Vertex& blueVertex = blueVertices[(i % (2 * n)) / 2 + n * (((i + 2*n) / 4) / n)];
      Vertex& redVertex = redVertices[(1+(i % (2 * n))) / 2 + (n + 1) * ((i / 4) / (n))];
      redEdges[i].setVertices(redVertex, blueVertex);
      blueEdges[i].setVertices(blueVertex, redVertex);
      redVertex.addEdge(redEdges[i]);
      blueVertex.addEdge(blueEdges[i]);
      redEdges[i].X = 0;
      blueEdges[i].X = 0;
      redEdges[i].oldX = 0;
      blueEdges[i].oldX = 0;
      double w = r.variate<double, std::exponential_distribution>(1);
      redEdges[i].w = w;
      blueEdges[i].w = w;
    }
  }

  void propagateBeliefs() {
    for (HalfEdge& e : redEdges) {
      double Xt = std::numeric_limits<double>::infinity();
      const Vertex* head = e.getHead();
      for (const HalfEdge* en : head->neighbors) {
        if (e.getTail()->getIndex() != en->getHead()->getIndex()) {
          Xt = std::min(en->w - en->oldX, Xt);
        }
      }
      e.X = Xt;
    }
    for (HalfEdge& e : blueEdges) {
      double Xt = std::numeric_limits<double>::infinity();
      const Vertex* head = e.getHead();
      for (const HalfEdge* en : head->neighbors) {
        if (e.getTail()->getIndex() != en->getHead()->getIndex()) {
          Xt = std::min(en->w - en->oldX, Xt);
        }
      }
      e.X = Xt;
    }
    for (HalfEdge& e : redEdges) {
      e.oldX = e.X;
    }
    for (HalfEdge& e : blueEdges) {
      e.oldX = e.X;
    }
  }
};

int main() {
  Rng r;
  BipartiteGraph G(10, r);

  for (unsigned i = 0; i < 1000; i++) {
    G.propagateBeliefs();
//    std::cout << G.redEdges[40].X << " " << G.blueEdges[40].X << std::endl;
  }

  for (unsigned i = 0; i < G.blueEdges.size(); i++) {
    std::cout << G.blueEdges[i].getTail()->getIndex() << " " << G.blueEdges[i].getHead()->getIndex() + G.blueVertices.size() << std::endl;
  }
  for (unsigned i = 0; i < G.blueEdges.size(); i++) {
    if (G.blueEdges[i].X + G.redEdges[i].X >= 0) {
      std::cout << G.blueEdges[i].getTail()->getIndex() << " " << G.blueEdges[i].getHead()->getIndex() + G.blueVertices.size() << std::endl;
    }
  }

}