#include "hadamard_mcmc.hpp"
#include "matrices.hpp"

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

std::string getFilename(unsigned n, double β1, double β2, double θ₀, unsigned m, unsigned s) {
  auto tag = std::chrono::high_resolution_clock::now();
  return "aging_" + std::to_string(n) + "_" + std::to_string(β1) + "_" + std::to_string(β2) + "_" +
    std::to_string(θ₀) + "_" + std::to_string(m) + "_" + std::to_string(s) + "_" + std::to_string(tag.time_since_epoch().count()) + ".dat";
}

class Correlation {
  public:
    Orthogonal M0;
    std::list<double> C;

    Correlation(const Orthogonal& M) : M0(M) {}

    void add(const Orthogonal& M) {
      C.push_back(M0 * M);
    }
};

class MeasureCorrelation : public Measurement {
public:
  unsigned s;
  unsigned n;
  unsigned N;
  std::list<Correlation> Cs;

  MeasureCorrelation(unsigned m, unsigned skip) : N(m), s(skip) {
    n = 0;
  }

  void after_sweep(double, double, const Orthogonal& M) override {
    if (n % N == 0) {
      Cs.push_back(Correlation(M));
    }

    if (n % s == 0) {
      for (Correlation& C : Cs) {
        C.add(M);
      }
    }

    n++;
  }
};

int main(int argc, char* argv[]) {
  // model parameters
  unsigned n = 2; // matrix size over four

  // simulation settings
  double β1 = 5; // temperature
  double β2 = 15; // temperature
  unsigned M = 1e4; // number of relaxation sweeps
  unsigned N = 1e4; // number of measurement sweeps
  unsigned m = 1e3;
  unsigned s = 10;
  double θ₀ = 0.05; // standard deviation of step size

  bool loadDataFromFile = false;

  int opt;

  while ((opt = getopt(argc, argv, "n:b:c:M:N:m:t:s:")) != -1) {
    switch (opt) {
    case 'n':
      n = atoi(optarg);
      break;
    case 'b':
      β1 = atof(optarg);
      break;
    case 'c':
      β2 = atof(optarg);
      break;
    case 'M':
      M = (unsigned)atof(optarg);
      break;
    case 'N':
      N = (unsigned)atof(optarg);
      break;
    case 'm':
      m = (unsigned)atof(optarg);
      break;
    case 't':
      θ₀ = atof(optarg);
      break;
    case 's':
      s = (unsigned)atof(optarg);
      break;
    default:
      exit(1);
    }
  }

  MeasureCorrelation Q(m, s);
  MCMC simulation(n, β1, Q, θ₀);

  simulation.run(M, true);
  std::cout << "Finished initial relaxation at " << β1 << ", quenching to " << β2 << "." << std::endl;

  simulation.β = β2;
  simulation.run(N * m);

  std::cout << "Finished simulation." << std::endl;

  std::string filename = getFilename(n, β1, β2, θ₀, m, s);
  std::ofstream file(filename);

  for (const Correlation& C : Q.Cs) {
    for (double x : C.C) {
      file << x << " ";
    }
    file << std::endl;
  }

  return 0;
}