#pragma once

#include <array>
#include <functional>
#include <list>
#include <queue>
#include <vector>

#include "euclidean.hpp"
#include "octree.hpp"
#include "quantity.hpp"
#include "spin.hpp"
#include "random.hpp"
#include "transformation.hpp"

template <class U, int D, class R, class S> class Model;

template <class U, int D, class R, class S> class measurement {
public:
  virtual void pre_cluster(const Model<U, D, R, S>&, unsigned, const R&){};
  virtual void plain_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>*,
                                  const Spin<U, D, S>*, const Spin<U, D, S>&, double){};
  virtual void plain_site_transformed(const Model<U, D, R, S>&, const Spin<U, D, S>*,
                                      const Spin<U, D, S>&){};

  virtual void ghost_bond_visited(const Model<U, D, R, S>&, const Spin<U, D, S>&,
                                  const Spin<U, D, S>&, double){};
  virtual void ghost_site_transformed(const Model<U, D, R, S>&, const R&){};

  virtual void post_cluster(const Model<U, D, R, S>&){};
};

template <class U, int D, class R, class S> class Model {
public:
  U L;
  R s0;
  std::vector<Spin<U, D, S>*> s;
  Octree<U, D, S> dict;
  std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z;
  std::function<double(const Spin<U, D, S>&)> B;
  Rng rng;

  Model(U L, std::function<double(const Spin<U, D, S>&, const Spin<U, D, S>&)> Z,
        std::function<double(const Spin<U, D, S>&)> B)
      : L(L), s0(L), Z(Z), B(B), rng(), dict(L, std::floor(L)) {}

  void step(double T, Transformation<U, D, R, S>* t0, measurement<U, D, R, S>& A) {
    std::queue<Transformation<U, D, R, S>*> queue;
    queue.push(t0);

    std::set<Spin<U, D, S>*> visited;
    std::set<Spin<U, D, S>*> fused;

    while (!queue.empty()) {
      Transformation<U, D, R, S>* t = queue.front();
      queue.pop();

      std::set<Spin<U, D, S>*> c = t->current();
      if (!std::any_of(c.begin(), c.end(), [&visited](Spin<U, D, S>* s) { return visited.contains(s); }) && !(c.size() == 1 && std::any_of(c.begin(), c.end(), [&fused](Spin<U, D, S>* s) { return fused.contains(s); }))) {
        visited.insert(c.begin(), c.end());
        fused.insert(c.begin(), c.end());

        for (Spin<U, D, S>* s : t->toConsider()) {
          if (!fused.contains(s)) {
          double ΔE = t->ΔE(s);
          if (ΔE > 0 && 1.0 - exp(-ΔE / T) > rng.uniform<double>(0, 1)) {
            queue.push(t->createNew(T, fused, s, rng));
          }
          }
        }

        t->apply();
      }
      delete t;
    }
  }

  void resize(double T, double P) {}

  void wolff(double T, std::vector<Gen<U, D, R, S>> gs,
             measurement<U, D, R, S>& A, unsigned N) {
    for (unsigned i = 0; i < N; i++) {
      Gen<U, D, R, S>& g = rng.pick(gs);
      Transformation<U, D, R, S> *t = g(*this, rng);

      this->step(T, t, A);

      A.post_cluster(*this);
    }
  }
};