#include <vector>
#include <list>
#include <set>
#include <iostream>
#include <fstream>
#include <functional>
#include <random>
#include <queue>
#include <eigen3/Eigen/Dense>
#include "randutils/randutils.hpp"

const std::array<std::array<unsigned, 16>, 16> smiley = 
  {{
    {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}},
    {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
    {{0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0}},
    {{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}},
    {{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1}},
    {{1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1}},
    {{1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1}},
    {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
    {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}},
    {{1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1}},
    {{1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1}},
    {{1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1}},
    {{0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0}},
    {{0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0}},
    {{0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0}},
    {{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0}}
  }};

template <class T, int D>
using vector = Eigen::Matrix<T, D, 1>;

template <class T, int D>
using matrix = Eigen::Matrix<T, D, D>;

template <class T, int D, class state>
class spin {
  public:
    vector<T, D> x;
    state s;
};

template <class T, int D>
class euclidean {
  private:
    unsigned L;

  public:
    vector<T, D> t;
    matrix<T, D> r;
    euclidean(unsigned L) : L(L) {
      for (unsigned i = 0; i < D; i++) {
        t(i) = 0;
        r(i, i) = 1;
        for (unsigned j = 1; j < D; j++) {
          r(i, (i + j) % D) = 0;
        }
      }
    }

    euclidean(unsigned L, vector<T, D> t0, matrix<T, D> r0) : L(L) {
      t = t0;
      r = r0;
    }

    template <class state>
    spin<T, D, state> act(const spin<T, D, state>& s) const {
      spin<T, D, state> s_new;

      s_new.x = t + r * s.x;
      s_new.s = s.s;

      for (unsigned i = 0; i < D; i++) {
        s_new.x(i) = fmod(L + s_new.x(i), L);
      }

      return s_new;
    }

    euclidean act(const euclidean& x) const {
      vector<T, D> tnew = r * x.t + t;
      matrix<T, D> rnew = r * x.r;

      euclidean pnew(this->L, tnew, rnew);

      return pnew;
    }

    euclidean inverse() const {
      vector<T, D> tnew = - r.transpose() * t;
      matrix<T, D> rnew = r.transpose();

      euclidean pnew(this->L, tnew, rnew);

      return pnew;
    }
};

// to-do: clean up nnn... code
template <int D>
class dictionary {
  private:
    unsigned L;
    std::vector<std::set<unsigned>> d;

  public:
    dictionary(unsigned Li) : L(Li), d(pow(Li, D)) {};

    template <class T>
    unsigned dictionary_index(vector<T, D> x) const {
      unsigned pos_ind = 0;

      for (unsigned i = 0; i < D; i++) {
        pos_ind += pow(L, i) * (unsigned)x(i);
      };

      return pos_ind;
    }

    template <class T>
    void record(vector<T, D> x, unsigned ind) {
      d[dictionary_index<T>(x)].insert(ind);
    };

    template <class T>
    void remove(vector<T, D> x, unsigned ind) {
      d[dictionary_index<T>(x)].erase(ind);
    };

    template <class T>
    std::set<unsigned> neighbors(vector<T, D> x, unsigned depth) const {
      return nearest_neighbors_of(dictionary_index<T>(x), depth, {});
    };

    std::set<unsigned> nearest_neighbors_of(unsigned ind, unsigned depth, std::list<unsigned> ignores) const {
      std::set<unsigned> ns = d[ind];

      if (depth > 0) {
        for (unsigned i = 0; i < D; i++) {
          if (std::none_of(ignores.begin(), ignores.end(), [i](unsigned k){return i==k;})) {
            std::list<unsigned> ignores_new = ignores;
            ignores_new.push_back(i);
            for (signed j : {-1, 1}) {
              unsigned ni = pow(L, i + 1) * (ind / ((unsigned)pow(L, i + 1))) + fmod(pow(L, i + 1) + ind + j * pow(L, i), pow(L, i + 1));

              std::set<unsigned> nns = nearest_neighbors_of(ni, depth - 1, ignores_new);

              for (unsigned guy : nns) {
                ns.insert(guy);
              }
            }
          }
        }
      }

      return ns;
    };
};

class quantity {
  private:
    double total;
    double total2;
    std::vector<double> C;
    unsigned wait;

  public:
    unsigned n;
    std::list<double> hist;

    quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
      n = 0;
      total = 0;
      total2 = 0;
    }

    void add(double x) {
      if (n > wait) {
        hist.push_front(x);
        if (hist.size() > C.size()) {
          hist.pop_back();
          unsigned t = 0;
          for (double a : hist) {
            C[t] += a * x;
            t++;
          }
          total += x;
          total2 += pow(x, 2);
        }
      }
      n++;
    }

    double avg() const {
      return total / (n - wait);
    }

    double avg2() const {
      return total2 / (n - wait);
    }

    std::vector<double> ρ() const {
      double C0 = C.front() / (n - wait);
      double avg2 = pow(total / (n - wait), 2);

      std::vector<double> ρtmp;

      for (double Ct : C) {
        ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
      }

      return ρtmp;
    }

    std::array<double, 2> τ() const {
      double τtmp = 0.5;
      unsigned M = 1;
      double c = 8.0;

      std::vector<double> ρ_tmp = this->ρ();

      while (c * τtmp > M && M < C.size()) {
        τtmp += ρ_tmp[M];
        M++;
      }

      return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
    }

    double σ() const {
      return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
    }

    double serr() const {
      return sqrt(this->σ());
    }
};

template <class U, int D, class state>
class model {
  public:
    unsigned L;
    euclidean<U, D> s0;
    std::vector<spin<U, D, state>> s;
    dictionary<D> dict;
    std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors;
    std::function<double(spin<U, D, state>, spin<U, D, state>)> Z;
    std::function<double(spin<U, D, state>)> B;
    long double E;
    quantity Eq;
    quantity Cq;

    model(unsigned L, std::function<double(spin<U, D, state>, spin<U, D, state>)> Z,
        std::function<double(spin<U, D, state>)> B,
        std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns) : 
      L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) {
        }

    void update_energy() {
      E = 0;
      for (unsigned i = 0; i < s.size(); i++) {
        for (unsigned j = 0; j < i; j++) {
          E -= Z(s[i], s[j]);
        }

        E -= B(s0.inverse().act(s[i]));
      }
    }

    void step(double T, unsigned ind, euclidean<U, D> r, std::mt19937& rng) {
      unsigned cluster_size = 0;
      std::uniform_real_distribution<double> dist(0.0, 1.0);

      std::queue<unsigned> queue;
      queue.push(ind);

      std::vector<bool> visited(s.size() + 1, false);

      while (!queue.empty()) {
        unsigned i = queue.front();
        queue.pop();

        if (!visited[i]) {
          visited[i] = true;

          bool we_are_ghost = i == s.size();

          spin<U, D, state> si_new;
          euclidean<U, D> s0_new(L);

          if (we_are_ghost) {
            s0_new = r.act(s0);
          } else {
            si_new = r.act(s[i]);
          }

          for (unsigned j : neighbors(*this, i, si_new)) {
            if (j != i) {
              double dE;
              bool neighbor_is_ghost = j == s.size();

              if (we_are_ghost || neighbor_is_ghost) {
                spin<U, D, state> s0s_old, s0s_new;
                unsigned non_ghost;

                if (neighbor_is_ghost) {
                  non_ghost = i;
                  s0s_old = s0.inverse().act(s[i]);
                  s0s_new = s0.inverse().act(si_new);
                } else {
                  non_ghost = j;
                  s0s_old = s0.inverse().act(s[j]);
                  s0s_new = s0_new.inverse().act(s[j]);
                }

                dE = B(s0s_old) - B(s0s_new);
              } else {
                dE = Z(s[i], s[j]) - Z(si_new, s[j]);
              }

              double p = 1.0 - exp(-dE / T);

              if (dist(rng) < p) {
                queue.push(j);
              }
            }
          }

          if (we_are_ghost) {
            s0 = s0_new;
          } else {
            dict.remove(s[i].x, i);
            s[i] = si_new;
            dict.record(s[i].x, i);
            cluster_size++;
          }
        }
      }

      Cq.add(cluster_size);
    }

    void wolff(double T, unsigned N, std::mt19937& rng) {
      std::uniform_real_distribution<double> t_dist(0, L);
      std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
      std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);

      for (unsigned i = 0; i < N; i++) {
        vector<U, D> t;
        matrix<U, D> m;
        for (unsigned j = 0; j < D; j++) {
          t(j) = (U)t_dist(rng);
        }

        unsigned flip_D1 = r_dist(rng);
        unsigned flip_D2 = r_dist(rng);

        for (unsigned j = 0; j < D; j++) {
          for (unsigned k = 0; k < D; k++) {
            if ((j == flip_D1 && k == flip_D2) || (j == flip_D2 && k == flip_D1)) {
              if (flip_D1 <= flip_D2) {
                m(j, k) = -1;
              } else {
                m(j, k) = 1;
              }
            } else if ((j == k && j != flip_D1) && j != flip_D2) {
              m(j, k) = 1;
            } else {
              m(j, k) = 0;
            }
          }
        }

        euclidean<U, D> g(L, t, m);

        this->step(T, ind_dist(rng), g, rng);

        this->update_energy();
        Eq.add(E);
      }
    }
};