#pragma once

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

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

template <class U, int D, class R, class S> class Model {
public:
  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)
      : 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>*> cluster = t0->current();

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

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

      A.plain_site_transformed(*this, *t);
      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);

      A.pre_cluster(*this, i, t);

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

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