#define WOLFF_USE_FINITE_STATES

#include <cstring>
#include <fstream>
#include <iostream>
#include <sstream>

#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"
#include "wolff/lib/wolff_models/ising.hpp"

using namespace wolff;

template <typename T> std::string to_string_with_precision(const T a_value, const int n = 12) {
  std::ostringstream out;
  out.precision(n);
  out << std::fixed << a_value;
  return out.str();
}

void update_distribution_file(std::string id, const std::vector<uint64_t>& data,
                              std::string model_string) {
  std::string filename = model_string + id + ".dat";
  std::ifstream file(filename);

  std::vector<uint64_t> data_old(data.size(), 0);

  if (file.is_open()) {
    for (unsigned i = 0; i < data.size(); i++) {
      uint64_t num;
      file >> num;
      data_old[i] = num;
    }

    file.close();
  }

  std::ofstream file_out(filename);

  for (unsigned i = 0; i < data.size(); i++) {
    file_out << std::fixed << data_old[i] + data[i] << " ";
  }

  file_out.close();
}

void update_distribution_file(std::string id, const std::vector<std::vector<uint64_t>>& data,
                              std::string model_string) {
  std::string filename = model_string + id + ".dat";
  std::ifstream file(filename);

  std::vector<std::vector<uint64_t>> data_old(data.size());

  for (unsigned i = 0; i < data.size(); i++) {
    data_old[i].resize(data[0].size());
  }

  if (file.is_open()) {
    for (unsigned i = 0; i < data.size(); i++) {
      for (unsigned j = 0; j < data[0].size(); j++) {
        uint64_t num;
        file >> num;
        data_old[i][j] = num;
      }
    }

    file.close();
  }

  std::ofstream file_out(filename);

  for (unsigned i = 0; i < data.size(); i++) {
    for (unsigned j = 0; j < data[0].size(); j++) {
      file_out << std::fixed << data_old[i][j] + data[i][j] << " ";
    }
    file_out << "\n";
  }

  file_out.close();
}

class meas : public measurement<ising_t, ising_t, graph<>> {
private:
  std::vector<uint64_t> mag_dist;
  std::vector<std::vector<uint64_t>> energy_mag_dist;

  signed E;
  signed M;

  std::string model_string;

public:
  meas(const wolff::system<ising_t, ising_t, graph<>>& S, double H)
      : mag_dist(S.nv + 1, 0), energy_mag_dist(S.nv + 1) {
    M = S.nv * S.s[0];
    E = S.ne;
    for (std::vector<uint64_t>& d : energy_mag_dist) {
      d.resize(S.ne + 1, 0);
    }
    model_string = "metastable_" + std::to_string(S.G.D) + "_" + std::to_string(S.G.L) + "_" +
                   to_string_with_precision(S.T) + "_" + to_string_with_precision(H) + "_";
  }

  ~meas() {
    update_distribution_file("m", mag_dist, model_string);
    update_distribution_file("e", energy_mag_dist, model_string);
  }

  void plain_bond_visited(const wolff::system<ising_t, ising_t, graph<>>&,
                          const typename graph<>::halfedge&, const ising_t&, double dE) override {
    if (dE > 0) {
      E -= 2;
    } else {
      E += 2;
    }
  }

  void ghost_bond_visited(const wolff::system<ising_t, ising_t, graph<>>&, const typename graph<>::vertex&,
                          const ising_t& s_old, const ising_t& s_new, double dE) override {
    M += s_new - s_old;
  }

  void post_cluster(unsigned, unsigned, const wolff::system<ising_t, ising_t, graph<>>& S) override {
    mag_dist[(S.nv + M) / 2]++;
    energy_mag_dist[(S.nv + M) / 2][(E + S.ne) / 2]++;
  }
};

int main(int argc, char* argv[]) {
  unsigned N = (unsigned)1e4;
  unsigned D = 2;
  unsigned L = 128;
  double T = 2.26918531421;
  double H = 0.0;

  int opt;

  while ((opt = getopt(argc, argv, "N:D:L:T:H:")) != -1) {
    switch (opt) {
    case 'N': // number of steps
      N = (unsigned)atof(optarg);
      break;
    case 'D': // dimension
      D = atoi(optarg);
      break;
    case 'L': // linear size
      L = atoi(optarg);
      break;
    case 'T': // temperature
      T = atof(optarg);
      break;
    case 'H': // external field
      H = atof(optarg);
      break;
    default:
      exit(EXIT_FAILURE);
    }
  }

  // define the spin-spin coupling
  std::function<double(const ising_t&, const ising_t&)> Z =
      [](const ising_t& s1, const ising_t& s2) -> double { return (double)(s1 * s2); };

  // define the spin-field coupling
  std::function<double(const ising_t&)> B = [=](const ising_t& s) -> double { return H * s; };

  // initialize the lattice
  graph<> G(D, L);

  // initialize the system
  wolff::system<ising_t, ising_t, graph<>> S(G, T, Z, B);

  randutils::auto_seed_128 seeds;
  std::mt19937 rng(seeds);

  meas A(S, H);

  // run wolff N times
  S.run_wolff(N, gen_ising<graph<>>, A, rng);

  return 0;
}