#include <list>
#include <queue>

#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/src/Core/CwiseNullaryOp.h>

#include "pcg-cpp/include/pcg_random.hpp"
#include "randutils/randutils.hpp"

template <int D> using Vector = Eigen::Matrix<int, D, 1>;
template <int D> using Matrix = Eigen::Matrix<int, D, D>;

using Rng = randutils::random_generator<pcg32>;

template <unsigned D> Vector<D> randomVector(unsigned L, Rng& r) {
  Vector<D> x;
  for (unsigned i = 0; i < D; i++) {
    x[i] = r.uniform((unsigned)0, L - 1);
  }
  return x;
}

int iPow(int x, unsigned p) {
  if (p == 0)
    return 1;
  if (p == 1)
    return x;

  int tmp = iPow(x, p / 2);
  if (p % 2 == 0) {
    return tmp * tmp;
  } else {
    return x * tmp * tmp;
  }
}

unsigned mod(signed a, unsigned b) { return ((a < 0) ? (a + (1 - a / (signed)b) * b) : a) % b; }

template <int D, typename Derived> Vector<D> mod(const Eigen::MatrixBase<Derived>& v, unsigned b) {
  Vector<D> u;
  for (unsigned i = 0; i < D; i++) {
    u(i) = mod(v(i), b);
  }
  return u;
}

template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
  if (level > 0) {
    unsigned new_level = level - 1;
    unsigned old_length = sequences.size();
    for (std::array<double, D>& sequence : sequences) {
      std::array<double, D> new_sequence = sequence;
      new_sequence[new_level] = -1;
      sequences.push_front(new_sequence);
    }
    one_sequences<D>(sequences, new_level);
  }
}

template <unsigned D> std::vector<Matrix<D>> generateTorusMatrices() {
  std::vector<Matrix<D>> mats;

  std::array<double, D> ini_sequence;
  ini_sequence.fill(1);
  std::list<std::array<double, D>> sequences;
  sequences.push_back(ini_sequence);

  one_sequences<D>(sequences, D);

  sequences.pop_back(); // don't want the identity matrix!

  for (std::array<double, D> sequence : sequences) {
    Matrix<D> m;
    for (unsigned i = 0; i < D; i++) {
      for (unsigned j = 0; j < D; j++) {
        if (i == j) {
          m(i, j) = sequence[i];
        } else {
          m(i, j) = 0;
        }
      }
    }

    mats.push_back(m);
  }

  for (unsigned i = 0; i < D; i++) {
    for (unsigned j = 0; j < D; j++) {
      if (i != j) {
        Matrix<D> m;
        for (unsigned k = 0; k < D; k++) {
          for (unsigned l = 0; l < D; l++) {
            if ((k == i && l == j) || (k == j && l == i)) {
              if (i < j) {
                m(k, l) = 1;
              } else {
                m(k, l) = -1;
              }
            } else if (k == l && (k != i && k != j)) {
              m(k, l) = 1;
            } else {
              m(k, l) = 0;
            }
          }
        }
        mats.push_back(m);
      }
    }
  }

  return mats;
}

template <unsigned D> class CiamarraState;

template <unsigned D> class Transformation {
public:
  unsigned L;
  Matrix<D> m;
  Vector<D> v;

  Transformation(unsigned L) : L(L) {
    m.setIdentity();
    v.setZero();
  }

  Transformation(unsigned L, const Matrix<D>& m, const Vector<D>& v) : L(L), m(m), v(v) {}

  Transformation(unsigned L, const std::vector<Matrix<D>>& ms, Rng& r) : m(r.pick(ms)), L(L) {
    for (unsigned i = 0; i < D; i++) {
      v[i] = r.uniform((unsigned)0, L - 1);
    }

    v = v - m * v;
  }

  Vector<D> apply(const Vector<D>& x) const { return mod<D>(v + m * x, L); }

  Transformation<D> apply(const Transformation<D>& t) const {
    Transformation<D> tNew(L);

    tNew.m = m * t.m;
    tNew.v = apply(t.v);

    return tNew;
  }

  Transformation<D> inverse() const {
    return Transformation<D>(L, m.transpose(), -m.transpose() * v);
  }

  CiamarraState<D> apply(const CiamarraState<D>& s) const {
    CiamarraState<D> s2(m * s);
    return s2;
  }

  template <class A> A apply(const A& s) const {
    return s;
  };
};

template <typename T> concept State = requires(T v, const T& v2) {
  { v.empty() } -> std::same_as<bool>;
  {v.remove()};
  { v.operator==(v2) } -> std::same_as<bool>;
};

template <unsigned D, State S> class HalfEdge;

template <unsigned D, State S> class Vertex {
public:
  Vector<D> position;
  S state;
  std::vector<HalfEdge<D, S>> adjacentEdges;
  bool marked;

  bool empty() const { return state.empty(); }
};

template <unsigned D, State S> class HalfEdge {
public:
  Vertex<D, S>& neighbor;
  Vector<D> Δx;

  HalfEdge(Vertex<D, S>& n, const Vector<D>& d) : neighbor(n), Δx(d) {}
};

template <unsigned D, State S> class System {
public:
  const unsigned L;
  unsigned N;
  std::vector<Vertex<D, S>> vertices;
  Transformation<D> orientation;

  unsigned vectorToIndex(const Vector<D>& x) const {
    unsigned i = 0;
    for (unsigned d = 0; d < D; d++) {
      i += x[d] * iPow(L, d);
    }
    return i;
  }

  Vector<D> indexToVector(unsigned i) const {
    Vector<D> x;
    for (unsigned d = 0; d < D; d++) {
      x[d] = (i / iPow(L, d)) % L;
    }
    return x;
  }

  System(unsigned L) : L(L), N(0), vertices(iPow(L, D)), orientation(L) {
    for (unsigned i = 0; i < iPow(L, D); i++) {
      vertices[i].position = indexToVector(i);
      vertices[i].adjacentEdges.reserve(2 * D);
      vertices[i].marked = false;
    }

    for (unsigned d = 0; d < D; d++) {
      Vector<D> Δx = Vector<D>::Zero();
      Δx[d] = 1;
      for (signed i = 0; i < iPow(L, D); i++) {
        unsigned j = iPow(L, d + 1) * (i / iPow(L, d + 1)) + mod(i + iPow(L, d), pow(L, d + 1));
        vertices[i].adjacentEdges.push_back(HalfEdge<D, S>(vertices[j], Δx));
        vertices[j].adjacentEdges.push_back(HalfEdge<D, S>(vertices[i], -Δx));
      }
    }
  }

  unsigned size() const { return vertices.size(); }

  double density() const { return N / pow(L, D); }

  unsigned maxOccupation() const {
    return iPow(L, D);
  }

  int overlap(const System<D, S>& s) const {
    unsigned o = 0;
    Transformation<D> oRel = orientation.apply(s.orientation.inverse());

    for (unsigned i = 0; i < vertices.size(); i++) {
      if (!vertices[i].empty()) {
        S s1 = vertices[i].state; 
        S s2 = oRel.apply(s.vertices[vectorToIndex(oRel.inverse().apply(indexToVector(i)))].state);

        if (s1 == s2) {
          o += 1;
        }
      }
    }

    return o;
  }
};

template <unsigned D, State S> class SoftSystem : public System<D, S> {
public:
  using System<D, S>::System;

  virtual double pairEnergy(const S&, const S&) const { return 0; }

  double siteEnergy(const Vertex<D, S>& v) const {
    double E = 0;

    for (const HalfEdge<D, S>& e : v.adjacentEdges) {
      E += pairEnergy(v.state, e.neighbor.state);
    }

    return E;
  }

  double energy() const {
    double E = 0;

    for (const Vertex<D, S>& v : this->vertices) {
      E += siteEnergy(v);
    }

    return E / 2;
  }

  bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& v2, double T, Rng& r) {
    double E0 = siteEnergy(v1) + siteEnergy(v2);
    std::swap(v1.state, v2.state);
    double E1 = siteEnergy(v1) + siteEnergy(v2);
    double ΔE = E1 - E0;

    if (exp(-ΔE / T) > r.uniform(0.0, 1.0)) {
      /* Accept the swap. */
      return true;
    } else {
      /* Revert the swap. */
      std::swap(v1.state, v2.state);
      return false;
    }
  }

  bool tryRandomMove(double T, Rng& r) {
    Vertex<D, S>& v1 = r.pick(this->vertices);

    if (v1.empty()) {
      return false;
    }

    Vertex<D, S>& v2 = (r.pick(v1.adjacentEdges)).neighbor;

    if (!v2.empty()) {
      return false;
    }

    return trySwap(v1, v2, T, r);
  }

  bool tryRandomSwap(double T, Rng& r) {
    Vertex<D, S>& v1 = r.pick(this->vertices);
    Vertex<D, S>& v2 = r.pick(this->vertices);

    if (v1.state != v2.state) {
      return trySwap(v1, v2, T, r);
    } else {
      return false;
    }
  }

  void sweepLocal(Rng& r) {
    for (unsigned i = 0; i < iPow(this->L, D); i++) {
      tryRandomMove(r);
    }
  }

  void sweepSwap(Rng& r) {
    for (unsigned i = 0; i < iPow(this->L, D); i++) {
      tryRandomSwap(r);
    }
  }

  unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, double T, Rng& r,
                       bool dry = false) {
    std::queue<std::array<std::reference_wrapper<Vertex<D, S>>, 2>> q;
    Vector<D> x0New = R.apply(v0.position);
    Vertex<D, S>& v0New = SoftSystem::vertices[SoftSystem::vectorToIndex(x0New)];
    q.push({v0, v0New});

    unsigned n = 0;

    while (!q.empty()) {
      auto [vR, vNewR] = q.front();
      q.pop();

      Vertex<D, S>& v = vR;
      Vertex<D, S>& vNew = vNewR;

      if (!v.marked && !vNew.marked) {
        v.marked = true;
        vNew.marked = true;

        for (HalfEdge<D, S>& e : v.adjacentEdges) {
          Vertex<D, S>& vn = e.neighbor;

          Vector<D> xnNew = R.apply(vn.position);
          Vertex<D, S>& vnNew = SoftSystem::vertices[SoftSystem::vectorToIndex(xnNew)];

          if (!vn.marked && !vnNew.marked) {
            double E0 = pairEnergy(v.state, vn.state) + pairEnergy(vNew.state, vnNew.state);
            double E1 = pairEnergy(vNew.state, vn.state) + pairEnergy(v.state, vnNew.state);
            double ΔE = E1 - E0;

            if (exp(-ΔE / T) < r.uniform(0.0, 1.0)) {
              q.push({vn, vnNew});
            }
          }
        }

        if (!dry) {
          std::swap(v.state, vNew.state);
        }

        if (!v.empty()) {
          n += 1;
        }

        if (!vNew.empty()) {
          n += 1;
        }
      }
    }

    return n;
  }

  unsigned wolff(const Transformation<D>& R, double T, Rng& r) {
    unsigned n = flipCluster(R, r.pick(this->vertices), T, r);
    if (n > SoftSystem::N / 2) {
      SoftSystem::orientation = R.apply(SoftSystem::orientation);
    }
    for (Vertex<D, S>& v : this->vertices) {
      v.marked = false;
    }
    return n;
  }

  void swendsenWang(const Transformation<D>& R, double T, Rng& r) {
    for (Vertex<D, S>& v : this->vertices) {
      if (!v.marked) {
        bool dry = 0.5 < r.uniform(0.0, 1.0);
        flipCluster(R, v, T, r, dry);
      }
    }

    for (Vertex<D, S>& v : this->vertices) {
      v.marked = false;
    }
  }
};

template <unsigned D, State S> class HardSystem : public System<D, S> {
public:

  HardSystem(unsigned L) : System<D, S>(L) {}

  virtual std::list<std::reference_wrapper<Vertex<D, S>>> overlaps(Vertex<D, S>&, const S&,
                                                        bool = false) { return {}; }

  bool insert(Vertex<D, S>& v, const S& s) {
    if (overlaps(v, s).empty()) {
      v.state = s;
      HardSystem::N++;
      return true;
    } else {
      return false;
    }
  }

  bool tryDeletion(Vertex<D, S>& v) {
    if (v.empty()) {
      return false;
    } else {
      v.state.remove();
      HardSystem::N--;
      return true;
    }
  }

  bool tryRandomMove(Rng& r) {
    Vertex<D, S>& v = r.pick(HardSystem::vertices);
    S oldState = v.state;

    if (!tryDeletion(v)) {
      return false;
    }

    if (1.0 / (2.0 * D) > r.uniform(0.0, 1.0)) {
      for (HalfEdge<D, S>& e : v.adjacentEdges) {
        if (1 == e.Δx.dot(oldState)) {
          if (insert(e.neighbor, oldState.flip())) {
            return true;
          }
          break;
        }
      }
    } else {
      S newState(r);
      while (newState == oldState) {
        newState = S(r);
      }
      if (insert(v, newState)) {
        return true;
      }
    }
    v.state = oldState;
    HardSystem::N++;
    return false;
  }

  bool tryRandomNonlocalMove(Rng& r) {
    Vertex<D, S>& v = r.pick(HardSystem::vertices);
    S oldState = v.state;

    if (!tryDeletion(v)) {
      return false;
    }

    Vertex<D, S>& vNew = r.pick(HardSystem::vertices);
    S newState(r);

    if (insert(vNew, newState)) {
      return true;
    }

    v.state = oldState;
    HardSystem::N++;
    return false;
  }

  bool trySwap(Vertex<D, S>& v1, Vertex<D, S>& v2) {
    if (v1.state == v2.state) {
      return false;
    }
    if (overlaps(v1, v2.state, true).size() == 0 && overlaps(v2, v1.state, true).size() == 0) {
      std::swap(v1.state, v2.state);
      return true;
    } else {
      return false;
    }
  }

  bool tryRandomSwap(Rng& r) {
    Vertex<D, S>& v1 = r.pick(HardSystem::vertices);
    Vertex<D, S>& v2 = r.pick(HardSystem::vertices);

    return trySwap(v1, v2);
  }


  bool compatible() {
    for (Vertex<D, S>& v : HardSystem::vertices) {
      if (overlaps(v, v.state, true).size() > 0) {
        return false;
      }
    }

    return true;
  }

  void sweepGrandCanonical(double z, Rng& r) {
    for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) {
      if (0.5 < r.uniform(0.0, 1.0)) {
        double pIns = HardSystem::maxOccupation() * z / (HardSystem::N + 1);

        if (pIns > r.uniform(0.0, 1.0)) {
          while (true) {
            Vertex<D, S>& v = r.pick(HardSystem::vertices);
            if (v.empty()) {
              insert(v, S(r));
              break;
            }
          }
        }
      } else {

        double pDel = HardSystem::N / (z * HardSystem::maxOccupation());

        if (pDel > r.uniform(0.0, 1.0)) {
          tryDeletion(r.pick(HardSystem::vertices));
        }
      }

      tryRandomMove(r);
    }
  }

  void sweepLocal(Rng& r) {
    for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) {
      tryRandomMove(r);
    }
  }

  void sweepSwap(Rng& r) {
    for (unsigned i = 0; i < iPow(HardSystem::L, D); i++) {
      tryRandomSwap(r);
    }
  }

  unsigned flipCluster(const Transformation<D>& R, Vertex<D, S>& v0, bool dry = false) {
    std::queue<std::reference_wrapper<Vertex<D, S>>> q;
    q.push(v0);

    unsigned n = 0;

    while (!q.empty()) {
      Vertex<D, S>& v = q.front();
      q.pop();

      if (!v.marked) {
        Vector<D> xNew = R.apply(v.position);
        Vertex<D, S>& vNew = HardSystem::vertices[HardSystem::vectorToIndex(xNew)];

        S s = R.apply(v.state);
        S sNew = R.apply(vNew.state);

        if (s != vNew.state) {
          v.marked = true;
          vNew.marked = true;

          std::list<std::reference_wrapper<Vertex<D, S>>> overlaps1 = overlaps(vNew, s, true);
          std::list<std::reference_wrapper<Vertex<D, S>>> overlaps2 = overlaps(v, sNew, true);
          overlaps1.splice(overlaps1.begin(), overlaps2);

          for (Vertex<D, S>& vn : overlaps1) {
            if (!vn.marked) {
              q.push(vn);
            }
          }

          if (!dry) {
            v.state = sNew;
            vNew.state = s;
          }

          if (!v.empty()) {
            n += 1;
          }

          if (!vNew.empty()) {
            n += 1;
          }
        }
      }
    }

    return n;
  }

  unsigned wolff(const Transformation<D>& R, Vertex<D, S>& v0) {
    unsigned n = flipCluster(R, v0);
    if (n > HardSystem::N / 2) {
      HardSystem::orientation = R.apply(HardSystem::orientation);
    }
    for (Vertex<D, S>& v : this->vertices) {
      v.marked = false;
    }
    return n;
  }

  unsigned wolff(const Transformation<D>& R, Rng& r) {
    unsigned n = flipCluster(R, r.pick(this->vertices));
    if (n > HardSystem::N / 2) {
      HardSystem::orientation = R.apply(HardSystem::orientation);
    }
    for (Vertex<D, S>& v : this->vertices) {
      v.marked = false;
    }
    return n;
  }

  void swendsenWang(const Transformation<D>& R, Rng& r) {
    for (Vertex<D, S>& v : HardSystem::vertices) {
      if (!v.marked) {
        bool dry = 0.5 < r.uniform(0.0, 1.0);
        unsigned n = flipCluster(R, v, dry);
        if (n > pow(HardSystem::L, D) / 4 && !dry) {
          HardSystem::orientation = R.apply(HardSystem::orientation);
        }
      }
    }

    for (Vertex<D, S>& v : HardSystem::vertices) {
      v.marked = false;
    }
  }

};