#include <iostream>

#include "glass.hpp"
#include "quantity.hpp"

class DistinguishableState {
public:
  unsigned type;

  DistinguishableState() { type = 0; }

  DistinguishableState(unsigned t) { type = t; }

  bool empty() const { return type == 0; }
  bool operator==(const DistinguishableState& s) const { return type == s.type; }
  void remove() { type = 0; }
};

template <unsigned D>
class DistinguishableSystem : public SoftSystem<D, DistinguishableState> {
public:
  std::vector<double> interaction;

  DistinguishableSystem(unsigned L, unsigned N, Rng& r)
      : SoftSystem<D, DistinguishableState>(L), interaction(N * N) {
    DistinguishableSystem::N = N;

    for (double& V : interaction) {
      V = r.uniform(-0.5, 0.5);
    }

    for (unsigned i = 0; i < this->N; i++) {
      DistinguishableSystem::vertices[i].state.type = i + 1;
    }
  }

  double pairEnergy(const DistinguishableState& s1, const DistinguishableState& s2) const override {
    if (s1.empty() || s2.empty()) {
      return 0;
    } else {
      if (s1.type < s2.type) {
        return interaction[(s1.type - 1) * this->N + (s2.type - 1)];
      } else {
        return interaction[(s2.type - 1) * this->N + (s1.type - 1)];
      }
    }
  }
};

int main() {
  const unsigned D = 2;
  unsigned L = 40;
  unsigned N = 1560;
  double Tmin = 0;
  double Tmax = 0.4;
  double δT = 0.05;

  Rng r;

  DistinguishableSystem<D> s(L, N, r);

  std::vector<Matrix<D>> ms = generateTorusMatrices<D>();

  for (unsigned j = 0; j < 10 * s.vertices.size(); j++) {
    // unsigned nC = s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
    s.tryRandomSwap(1000, r);
  }

  for (unsigned i = 0; i < 1e5; i++) {
    for (unsigned j = 0; j < s.vertices.size(); j++) {
      s.tryRandomMove(Tmax, r);
    }
  }


  unsigned n = 1;
  for (double T = Tmax; T > Tmin; T -= δT) {

    /*
    for (unsigned i = 0; i < 1e5; i++) {
      Matrix<D> m = r.pick(ms);
      Vertex<D, DistinguishableState>& v = r.pick(s.vertices);
      s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
    }
    */
    /*
     */
    std::cout << T << " ";
    for (unsigned i = 0; i < 1e3; i++) {
      for (unsigned j = 0; j < s.vertices.size(); j++) {
        s.tryRandomSwap(T, r);
      }
    }
    s.setInitialPosition();
    DistinguishableSystem<D> s0 = s;
    auto start = std::chrono::high_resolution_clock::now();
//    Quantity<double> energy(1e3);
    for (unsigned i = 0; i < 1e5; i++) {
      /*
      for (unsigned j = 0; j < s.vertices.size(); j++) {
        s.tryRandomSwap(T, r);
      }
      */
      Matrix<D> m = r.pick(ms);
      Vertex<D, DistinguishableState>& v = r.pick(s.vertices);
      s.wolff(Transformation<D>(L, m, v.position - m * v.position), T, r);
      if (i % 10 == 0) {
//        energy.add(s.energy());
        auto stop = std::chrono::high_resolution_clock::now();
        auto duration = duration_cast<std::chrono::microseconds>(stop - start);
        std::cout << duration.count() << " " << s.overlap(s0) << " ";
      }
    }
    /*
    for (unsigned i = 0; i < 1e5; i++) {
      for (unsigned j = 0; j < s.vertices.size(); j++) {
        s.tryRandomMove(T, r);
      }
      std::cout << s.selfIntermediateScattering() << " ";
    }
    */
    std::cout << std::endl;
    /*
    std::vector<double> rho = energy.ρ();
    for (const double& x : rho) {
      std::cout << x << " ";
    }
    std::cout << std::endl;
    std::cerr << T << " " << s.energy() / N << std::endl;
    */
    //    s.sweepLocal(r);
    //    s.sweepSwap(r);
    //    s.swendsenWang(Transformation<D>(L, ms, r), r);
  }

  return 0;
}