#include <fstream>
#include <iostream>
#include <chrono>

#include "space_wolff.hpp"
#include "spheres.hpp"
#include "animation.hpp"

const unsigned D = 2;
typedef Model<double, D, Euclidean<double, D>, double> model;

class SaveFlip : public measurement<double, D, Euclidean<double, D>, Radius> {
  std::ofstream snapfile;
  unsigned n;

  public:
  SaveFlip() {}

  void pre_cluster(const Model<double, D, Euclidean<double, D>, Radius>& m, unsigned,
                   const Transformation<double, D, Euclidean<double, D>, Radius>* t) override {
    snapfile.open("sphere_flip.dat");
    n = 0;
    for (const Sphere<D>* s : m.s) {
      snapfile << s << " " << s->s << " " << s->x.transpose() << " ";
    }
    snapfile << "\n";
  }

  void plain_site_transformed(const Model<double, D, Euclidean<double, D>, Radius>& m,
                              const Transformation<double, D, Euclidean<double, D>, Radius>& t) override {
    for (const Sphere<D>* s : t.current()) {
      snapfile << s << " ";
    }
    snapfile << "\n";
    n++;
  }

  void post_cluster(const Model<double, D, Euclidean<double, D>, Radius>& m) override {
    for (const Sphere<D>* s : m.s) {
      snapfile << s << " " << s->s << " " << s->x.transpose() << " ";
    }
    snapfile << "\n";
    snapfile.close();
    std::cout << n << "\n";
    if (2 < n && n < 20) {
    getchar();
    }
  }

};

int main(int argc, char* argv[]) {
  const unsigned D = 2;

  double L = 32;
  unsigned N = 1000;
  double T = 2.0 / log(1.0 + sqrt(2.0));
  double H = 1.0;
  unsigned n = 25;
  unsigned wait = 1000;

  double k = 1e2;
  double a = 0.1;

  int opt;

  while ((opt = getopt(argc, argv, "n:N:L:T:H:a:k:w:")) != -1) {
    switch (opt) {
    case 'n':
      n = (unsigned)atof(optarg);
      break;
    case 'N':
      N = (unsigned)atof(optarg);
      break;
    case 'L':
      L = atof(optarg);
      break;
    case 'T':
      T = atof(optarg);
      break;
    case 'H':
      H = atof(optarg);
      break;
    case 'a':
      a = atof(optarg);
      break;
    case 'k':
      k = atof(optarg);
      break;
    case 'w':
      wait = atoi(optarg);
      break;
    default:
      exit(1);
    }
  }

  std::function<double(Spin<double, D, double>)> B_hard = [L, H](Spin<double, D, double> s) -> double {
    if (fabs(s.x(0)) < 3 * L / 4 && fabs(s.x(1)) < 3 * L / 4) {
      return 0;
    } else {
      return std::numeric_limits<double>::infinity();
    }
  };

  auto g1 = nudgeGen<D, Radius>(1);
  auto g2 = swapGen<D, Radius>(0.01);
  auto g3 = accrossGen<D, Radius>(0.1);
  auto g4 = centerGen<D, Radius>(0);

  auto tag = std::chrono::high_resolution_clock::now();

  std::string filename = "flips_" + std::to_string(n) + "_" + std::to_string(T) + "_" + std::to_string(H) + "_" + std::to_string(a) + "_" + std::to_string(k) + "_" +
                         std::to_string(tag.time_since_epoch().count()) + ".dat";

  std::ofstream file;
  file.open(filename);

  Animation<double, D, Euclidean<double, D>, Radius> A(L, 750, argc, argv, wait);
  model sphere(1.0, zSpheres<D>(a, k), bCenter<D, Radius>(H));

  Rng rng;

  sphere.s.resize(n);

  unsigned nx = floor(sqrt(n));
  for (unsigned i = 0; i < sphere.s.size(); i++) {
    Spin<double, 2, double>* ss = new Spin<double, 2, double>();
    ss->x = {(i / nx) * L / nx - L / 2, (i % nx) * L / nx - L / 2};
    ss->s = rng.pick({0.45, 0.45});
    sphere.s[i] = ss;
    sphere.dict.insert(ss);
  }

  measurement<double, D, Euclidean<double, D>, Radius> A_tmp;

  sphere.wolff(T, {g1, g2, g3, g4}, A_tmp, N);

  SaveFlip A_new;
  sphere.wolff(T, {g2}, A_new, 10000);


  return 0;
}