#pragma once
#include "hadamard_mcmc.hpp"
#include <execution>
#include <list>

void swap(MCMC& s1, MCMC& s2) {
  std::swap(s1.M, s2.M);
  std::swap(s1.E, s2.E);
}

class ParallelMeasurement {
public:
  virtual void after_step(bool, unsigned, double, double, const MCMC&, const MCMC&){};
};

class PT {
private:
  Rng rng;

public:
  std::vector<MCMC> Ms;
  ParallelMeasurement& B;
  std::vector<Measurement*>& As;

  PT(double β₀, double β₁, unsigned N, unsigned n, ParallelMeasurement& B,
     std::vector<Measurement*>& As)
      : B(B), As(As) {
    Ms.reserve(N);
    for (unsigned i = 0; i < N; i++) {
      double β = β₀ + i * (β₁ - β₀) / (N - 1);
      Ms.push_back(MCMC(n, β, *As[i]));
    }
  }

  double T(unsigned i) { return 1 / Ms[Ms.size() - i - 1].β; }

  bool step(unsigned i, unsigned j, bool dry = false) {
    double Δβ = Ms[i].β - Ms[j].β;
    double ΔE = Ms[i].E - Ms[j].E;

    bool accepted = Δβ * ΔE > 0 || exp(Δβ * ΔE) > rng.uniform((double)0.0, 1.0);

    if (accepted)
      swap(Ms[i], Ms[j]);

    if (!dry) {
      B.after_step(accepted, i, Δβ, ΔE, Ms[i], Ms[j]);
    }
    return accepted;
  }

  std::vector<double> tune(unsigned n0, unsigned m, double ε, double ε2, double y = 0.2) {
    unsigned n = n0;

    while (true) {
      std::vector<color> colors(Ms.size(), none);
      colors.front() = down;
      colors.back() = up;

      std::vector<unsigned> nu(Ms.size(), 0);
      std::vector<unsigned> nd(Ms.size(), 0);

      for (unsigned i = 0; i < n; i++) {
        std::for_each(std::execution::par_unseq, Ms.begin(), Ms.end(),
                      [m, ε](MCMC& M) { M.tune(m, ε); });

        for (unsigned k = 0; k < m * Ms.size(); k++) {
          unsigned j = rng.uniform((unsigned)0, (unsigned)(Ms.size() - 2));

          if (this->step(j, j + 1, true)) {
            std::swap(colors[j], colors[j + 1]);
            colors.front() = down;
            colors.back() = up;
          }
        }

        for (unsigned j = 0; j < Ms.size(); j++) {
          if (colors[j] == up) {
            nu[j]++;
          } else if (colors[j] == down) {
            nd[j]++;
          }
        }
      }

      std::vector<double> f(Ms.size());
      std::vector<double> δf(Ms.size());

      std::vector<unsigned> f_keep;
      double f_last = 0;
      for (unsigned i = 0; i < Ms.size(); i++) {
        f[i] = nu[i] / (double)(nu[i] + nd[i]);
        δf[i] = f[i] * sqrt(1.0/nu[i] + 1.0/(nu[i]+nd[i]));
        if (f[i] >= f_last) {
          f_keep.push_back(i);
          f_last = f[i];
        }
      }

      for (signed i = 0; i < f_keep.size() - 1; i++) {
        for (unsigned j = f_keep[i]; j < f_keep[i + 1]; j++) {
          f[j] = f[f_keep[i]] + (f[f_keep[i + 1]] - f[f_keep[i]]) /
                                    (Ms[f_keep[i + 1]].β - Ms[f_keep[i]].β) *
                                    (Ms[j].β - Ms[f_keep[i]].β);
          δf[j] = sqrt(pow(δf[f_keep[i]], 2) + pow(sqrt(pow(δf[f_keep[i + 1]], 2) + pow(δf[f_keep[i]], 2)) /
                                    (Ms[f_keep[i + 1]].β - Ms[f_keep[i]].β) *
                                    (Ms[j].β - Ms[f_keep[i]].β), 2));
        }
      }

      // measure the difference between f and a straight line
      double Δf² = 0;
      double δΔf⁴ = 0;

      for (unsigned j = 1; j < f.size(); j++) {
        Δf² += pow(f[j] - (j + 1.0) / (Ms.size() + 1.0), 2);
        δΔf⁴ += pow(2 * δf[j] * f[j], 2);
      }

      std::vector<double> η(Ms.size() - 1);

      for (unsigned i = 0; i < Ms.size() - 1; i++) {
        η[η.size() - i - 1] = sqrt((f[i + 1] - f[i])) / (1 / Ms[i].β - 1 / Ms[i + 1].β);
      }

      double C = 0;
      for (unsigned i = 0; i < η.size(); i++) {
        C += η[i] * (T(i + 1) - T(i));
      }

      std::vector<double> T1(Ms.size() - 2);

      double x = 0;
      unsigned j = 1;

      for (unsigned i = 0; i < η.size(); i++) {
        double xnew = x + η[i] * (T(i + 1) - T(i)) / C;
        while (j < xnew * η.size() && j < η.size()) {
          T1[j - 1] = T(i) + (j / (double)η.size() - x) / η[i] * C;
          j++;
        }
        x = xnew;
      }

      double err = 0;
      for (unsigned i = 0; i < T1.size(); i++) {
        err += fabs(T1[i] - 1 / Ms[Ms.size() - i - 2].β);
        Ms[Ms.size() - i - 2].β = (1 - y) * Ms[Ms.size() - i - 2].β + y / T1[i];
      }

      double Δf = sqrt(Δf²);
      double δΔf = sqrt(δΔf⁴) / 2 / Δf;
      double relErr = Δf / Ms.size();
      // double relErr = err / T1.size() * Ms.size() / (1 / Ms.front().β - 1 / Ms.back().β);

      std::cout << "RMS difference from ideal transit flow is " << relErr << " ± " << δΔf / Ms.size() << ".\n";
      if (relErr < ε2 && !(δΔf != δΔf)) {
        return f;
      } else {
        if (5 * δΔf > Δf || δΔf != δΔf) {
          n *= 2;
          std::cout << "Error in RMS difference too close to difference, increasing tuning to " << n << ".\n";
        }
      }
    }
  }

  void run(unsigned n, unsigned m, bool dry = false) {
    for (unsigned i = 0; i < n; i++) {
      std::for_each(std::execution::par_unseq, Ms.begin(), Ms.end(),
                    [m, dry](MCMC& M) { M.run(m, dry); });
      for (unsigned j = 0; j < Ms.size(); j++) {
        unsigned k = rng.uniform((unsigned)0, (unsigned)(Ms.size() - 2));
        this->step(k, k + 1, dry);
      }
    }
  }
};