#include <iostream>
#include <fstream>

#include "glass.hpp"

void print(const BiroliSystem<2>& s) {
  for (const Vertex<2, BiroliState>& v : s.vertices) {
      std::cerr << v.state.type;

    if (v.position(0) == s.L - 1) {
      std::cerr << std::endl;
    }
  }
}

int main() {
  const unsigned D = 3;
  unsigned L = 30;
  unsigned Nmin = 2e2;
  unsigned Nmax = 2e5;
  double Tmin = 0.04;
  double Tmax = 0.2;
  double δT = 0.02;

  BiroliSystem<D> s(L);

  Rng r;

  double z = exp(1 / 0.2);

  if (!s.compatible()) {
    std::cerr << "Storted incompatible!" << std::endl;
    return 1;
  } 

  while (s.density() < 0.57) {
    s.sweepGrandCanonical(z, r);
  }

  if (!s.compatible()) {
    std::cerr << "Not compatible!" << std::endl;
    return 1;
  } 

  std::cerr << "Found state with appropriate density." << std::endl;

  BiroliSystem<D> s0 = s;

  std::vector<Matrix<D>> ms = generateTorusMatrices<D>();

  std::vector<unsigned> clusterDist(s.size() + 1);

  unsigned n = 1;
  unsigned i = 0;
  double nC = 0;
  while (nC / s.size() < 1e5) {
    if (n < 20 * log(i + 1)) {
      n++;
      std::cout << nC / s.size() << " " << (double)s.overlap(s0) / s.size() << std::endl;
    }
    unsigned nn = s.flipCluster(Transformation<D>(L, ms, r), r.pick(s.vertices));
    nC += nn;
    clusterDist[nn]++;
//    s.sweepLocal(r);
//    nC += s.size();
//    s.sweepSwap(r);
//    s.swendsenWang(Transformation<D>(L, ms, r), r);
    i++;
  }

  if (!s.compatible()) {
    std::cerr << "Not compatible!" << std::endl;
    return 1;
  } 

  std::ofstream file("dist.dat");
  for (unsigned i : clusterDist) {
    file << i << " ";
  }
  file.close();

  return 0;
}