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

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

#include "blossom5-v2.05.src/PerfectMatching.h"

using Rng = randutils::random_generator<pcg32>;

class Edge;
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 HalfEdge {
private:
  Vertex* tail;
  Vertex* head;

public:
  const Edge& edge;

  HalfEdge(const Edge& e) : edge(e) {}
  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]);
  }
};

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);
    }
  }
};

int main(int argc, char* argv[]) {
  unsigned n = 100;

  int opt;

  while ((opt = getopt(argc, argv, "n:")) != -1) {
    switch (opt) {
    case 'n':
      n = atoi(optarg);
      break;
    default:
      exit(1);
    }
  }

  Rng r;
  Graph G(n, r);

  PerfectMatching pm(G.vertices.size(), G.edges.size());

  for (const Edge& e : G.edges) {
    pm.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, e.weight);
  }

  pm.options.verbose = false;
  pm.Solve();

  std::reference_wrapper<Edge> eFlip = r.pick(G.edges);

  while (pm.GetMatch(eFlip.get().halfedges[0].getTail().index) != eFlip.get().halfedges[0].getHead().index) {
    eFlip = r.pick(G.edges);
  }

  eFlip.get().weight = 1e10;

  PerfectMatching pm2(G.vertices.size(), G.edges.size());

  for (const Edge& e : G.edges) {
    pm2.AddEdge(e.halfedges[0].getHead().index, e.halfedges[0].getTail().index, e.weight);
  }

  pm2.options.verbose = false;
  pm2.Solve();

  std::cout << n << std::endl;

  for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
    unsigned j1 = pm.GetMatch(i);

    std::cout
      << G.vertices[i].coordinate[0] << " "
      << G.vertices[i].coordinate[1] << " "
      << G.vertices[j1].coordinate[0] << " "
      << G.vertices[j1].coordinate[1] << std::endl;
  }

  unsigned m = 0;
  for (unsigned i = 0; i < G.vertices.size() / 2; i++) {
    unsigned j1 = pm.GetMatch(i);
    unsigned j2 = pm2.GetMatch(i);

    if (j1 != j2) {
      m++;
      std::cout
        << G.vertices[i].coordinate[0] << " "
        << G.vertices[i].coordinate[1] << " "
        << G.vertices[j2].coordinate[0] << " "
        << G.vertices[j2].coordinate[1] << std::endl;
    }
  }

  std::cerr << "Size of excitation: " << m << std::endl;

  return 0;
}