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

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

std::string getFilename(unsigned n, double β, double θ₀, unsigned m, unsigned lag, unsigned skip) {
  return "correlation_" + std::to_string(n) + "_" + std::to_string(β) + "_" +
    std::to_string(θ₀) + "_" + std::to_string(m) + "_" + std::to_string(lag) + "_" + std::to_string(skip)
    + ".dat";
}

class MeasureCorrelation : public Measurement {
public:
  Quantity<Orthogonal> Q;

  MeasureCorrelation(unsigned lag, unsigned skip) : Q(lag, skip) {}

  void after_sweep(double, double E, const Orthogonal& M) override {
    Q.add(M);
  }
};

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

  // simulation settings
  double β = 1; // temperature
  unsigned m = 1e4; // number of relaxation sweeps
  unsigned N = 1e4; // number of measurement sweeps
  double θ₀ = 0.05; // standard deviation of step size

  // measurement settings
  unsigned l = 1e3; // lag in autocorrelation function
  unsigned s = 1; // skip in autocorrelation function

  bool loadDataFromFile = false;

  int opt;

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

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

  if (loadDataFromFile) {
    Q.Q.read(getFilename(n, β, θ₀, m, l, s));
    std::cout << "Imported data from file, number of steps " << Q.Q.num_added() << "." << std::endl;
  }

  simulation.run(m, true);
  std::cout << "Finished initial relaxation at " << β << "." << std::endl;
  simulation.run(N);

  std::cout << "Finished simulation of " << β << ", writing output file." << std::endl;
  std::string filename = getFilename(n, β, θ₀, m, l, s);
  Q.Q.write(filename);

  return 0;
}