#include <stack>
#include <list>
#include <cstdlib>
#include <getopt.h>

#include "eigen/Eigen/Dense"
#include "randutils/randutils.hpp"
#include "pcg-cpp/include/pcg_random.hpp"

using Rng = randutils::random_generator<pcg32>;

using Real = double;
using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>;
using Matrix = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;

class Vertex;
class Edge;

class HalfEdge {
public:
  Edge& parent;
  Vertex& neighbor;
  HalfEdge(Edge& e, Vertex& v) : parent(e), neighbor(v) {};
};

class Vertex {
public:
  std::vector<HalfEdge> bonds;
  unsigned index;
  unsigned cluster;
};

class Edge {
public:
  bool active;
};

class Cluster : public std::list<std::reference_wrapper<Vertex>> {};
class Clusters {
public:
  std::vector<Cluster> lists;
  std::vector<unsigned> indices;
};

class Graph {
public:
  unsigned D;
  unsigned L;
  std::vector<Vertex> vertices;
  std::vector<Edge> edges;
  std::vector<Cluster> clusters;

  unsigned Nv() const {
    return vertices.size();
  }

  unsigned Ne() const {
    return edges.size();
  }

  unsigned squaredDistance(unsigned vi, unsigned vj) const {
    unsigned sd = 0;
    for (unsigned i = 0; i < D; i++) {
      int x1 = (unsigned)fmod((unsigned)(vi / (unsigned)pow(L, i)), pow(L, i + 1));
      int x2 = (unsigned)fmod((unsigned)(vj / (unsigned)pow(L, i)), pow(L, i + 1));
      unsigned Δx = abs(x1 - x2);
      if (Δx > L / 2) {
        Δx = L - Δx;
      }
      sd += pow(Δx, 2);
    }
    return sd;
  }

  /* Initialize a square lattice */
  Graph(unsigned D, unsigned L) : D(D), L(L), vertices(pow(L, D)), edges(D * pow(L, D)) {
    for (unsigned i = 0; i < Nv(); i++) {
      vertices[i].index = i;
      vertices[i].bonds.reserve(2 * D);
      for (unsigned j = 0; j < D; j++) {
        unsigned n1 = (unsigned)pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) +
                                          (unsigned)fmod(i + pow(L, j), pow(L, j + 1));
        unsigned n2 = (unsigned)pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) +
                          (unsigned)fmod(pow(L, j + 1) + i - pow(L, j), pow(L, j + 1));

        unsigned e1 = j * Nv() + i;
        unsigned e2 = j * Nv() + n2;

        HalfEdge he1(edges[e1], vertices[n1]);
        HalfEdge he2(edges[e2], vertices[n2]);

        vertices[i].bonds.push_back(he1);
        vertices[i].bonds.push_back(he2);
      }
    }
  }

  void dilute(Real p, Rng& r) {
    for (Edge& e : edges) {
      e.active = p > r.variate<Real, std::uniform_real_distribution>(0.0, 1.0);
    }
  }

  void markClusters() {
    clusters = {};
    for (Vertex& v : vertices) {
      v.cluster = 0;
    };
    unsigned nClusters = 0;
    for (Vertex& v : vertices) {
      if (v.cluster == 0) {
        nClusters++;
        clusters.push_back({});
        std::stack<std::reference_wrapper<Vertex>> inCluster;
        inCluster.push(v);
        while (!inCluster.empty()) {
          Vertex& vn = inCluster.top();
          inCluster.pop();
          if (vn.cluster == 0) {
            vn.cluster = nClusters;
            clusters[nClusters - 1].push_back(vn);
          }
          for (HalfEdge& e : vn.bonds) {
            if (e.parent.active && e.neighbor.cluster == 0) {
              inCluster.push(e.neighbor);
            }
          }
        }
      }
    }
  }

  Matrix laplacian() const {
    Matrix L = Matrix::Zero(Nv(), Nv());

    for (const Vertex& v : vertices) {
      for (const HalfEdge& e : v.bonds) {
        if (e.parent.active) {
          L(v.index, v.index) += 1;
          L(v.index, e.neighbor.index) -= 1;
        }
      }
    }

    return L;
  }
};

int main(int argc, char* argv[]) {
  Rng r;

  unsigned D = 2;
  unsigned L = 8;
  double p = 0.5;
  unsigned n = 10;

  int opt;

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

  std::vector<unsigned> multiplicities(D * L * L / 4);

  Graph G(D, L);

  for (const Vertex& v1 : G.vertices) {
    for (const Vertex& v2 : G.vertices) {
      unsigned dist = G.squaredDistance(v1.index, v2.index);
      if (dist > 0) {
        multiplicities[dist - 1]++;
      }
    }
  }

  for (unsigned trials = 0; trials < n; trials++) {
    G.dilute(p, r);
    G.markClusters();
    Matrix laplacian = G.laplacian();

    std::vector<Real> conductivities(L * L * D / 4);

    for (const Cluster& c : G.clusters) {
      std::vector<unsigned> inds;
      inds.reserve(c.size());
      for (const Vertex& v : c){
        inds.push_back(v.index);
      }

      Matrix subLaplacian = laplacian(inds, inds);
      subLaplacian(0,0)++; /* Set voltage of first node to zero */
      Matrix subLaplacianInv = subLaplacian.inverse();

      for (unsigned i = 0; i < c.size(); i++) {
        for (unsigned j = 0; j < i; j++) {
          Vector input = Vector::Zero(inds.size());
          input(i) = 1;
          input(j) = -1;

          Vector output = subLaplacianInv * input;
          double ΔV = std::abs(output(i) - output(j));

          conductivities[G.squaredDistance(inds[i], inds[j]) - 1] += 1 / ΔV;
        }
      }
    }

    for (unsigned i = 0; i < conductivities.size(); i++) {
      if (multiplicities[i] != 0) {
        std::cout << conductivities[i] / multiplicities[i] << " ";
      } else {
        std::cout << 0 << " ";
      }
    }
    std::cout << std::endl;
  }

  return 0;
}