#include <iostream>
#include <fstream>

#include "ising.hpp"
#include "quantity.hpp"

class Autocorrelation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
public:
  Quantity E;
  Autocorrelation(unsigned lag) : E(lag) {}

  void post_cluster(const isingModel& m) override {
    double E_tmp = 0;

    TorusGroup<signed, D> s0Inv = m.s0.inverse();

    for (const isingSpin* si : m.s) {
      E_tmp -= m.B(s0Inv.act(*si));

      for (const isingSpin* sj : m.dict.neighbors(si->x)) {
        if (si != sj) {
          E_tmp -= m.Z(*si, *sj) / 2;
        }
      }
    }

    E.add(E_tmp);
  }
};

int main(int argc, char* argv[]) {
  unsigned L = 32;
  unsigned N = 1000;
  unsigned mod = 0;
  unsigned multi = 1e4;
  double mag = 0.5;
  double pop = 1.0;
  double T = 2.0 / log(1.0 + sqrt(2.0));
  double H = 1.0;
  double ε = 0.1;
  unsigned lag = 1e2;

  int opt;

  while ((opt = getopt(argc, argv, "N:L:T:H:e:m:M:r:p:l:")) != -1) {
    switch (opt) {
    case 'N':
      N = (unsigned)atof(optarg);
      break;
    case 'L':
      L = atoi(optarg);
      break;
    case 'T':
      T = atof(optarg);
      break;
    case 'H':
      H = atof(optarg);
      break;
    case 'e':
      ε = atof(optarg);
      break;
    case 'm':
      mod = atoi(optarg);
      break;
    case 'M':
      multi = atoi(optarg);
      break;
    case 'r':
      mag = atof(optarg);
      break;
    case 'p':
      pop = atof(optarg);
      break;
    case 'l':
      lag = (unsigned)atof(optarg);
      break;
    default:
      exit(1);
    }
  }

  isingModel ising(L, isingZ(L), isingBMod(L, mod, H));
  isingPopulate(ising, L, pop, mag);

  auto g = isingGen(L);

  measurement<signed, D, TorusGroup<signed, D>, signed> A_empty;

  ising.wolff(T, {g}, A_empty, N);

  Autocorrelation A(lag);

  while (true) {
    ising.wolff(T, {g}, A, N);
    std::array<double, 2> τ = A.E.τ();
    std::cout << A.E.num_added() << " " << A.E.avg() << " " << τ[0] << " " << τ[1] << " " << τ[1] / τ[0] << "\n";
    if (τ[1] / τ[0] < ε && τ[0] * multi < A.E.num_added()) {
      break;
    }
  }

  std::ofstream outfile;
  outfile.open("out.dat", std::ios::app);

  std::array<double, 2> act = A.E.τ();
  std::vector<double> ρ = A.E.ρ();

  outfile << L << " " << T << " " << mod << " " << H << " " << A.E.num_added() << " " <<
  A.E.avg() << " " << A.E.serr() << " " << act[0] << " " << act[1]; for (double ρi : ρ) {
    outfile << " " << ρi;
  }
  outfile << "\n";

  return 0;
}