#include "title.hpp"
#include "smiley.hpp"
#include "space_wolff.hpp"
#include "torus_symmetries.hpp"

const unsigned D = 2;
typedef Model<signed, D, TorusGroup<signed, D>, signed> isingModel;
typedef Spin<signed, D, signed> isingSpin;

std::function<double(const isingSpin&, const isingSpin&)> isingZ(unsigned L) {
  return [L](const isingSpin& s1, const isingSpin& s2) -> double {
    bool old_one_one = false;
    bool old_many_ones = false;
    bool old_any_two = false;

    Vector<signed, D> old_diff = diff<signed, D>(L, s1.x, s2.x);

    for (unsigned i = 0; i < D; i++) {
      if (old_diff(i) == 1 && !old_one_one) {
        old_one_one = true;
      } else if (old_diff(i) == 1 && old_one_one) {
        old_many_ones = true;
      } else if (old_diff(i) > 1) {
        old_any_two = true;
      }
    }

    bool were_on_someone = !old_one_one && !old_any_two;
    bool were_nearest_neighbors = old_one_one && !old_many_ones && !old_any_two;

    if (were_on_someone) {
      return -std::numeric_limits<double>::infinity();
      ;
    } else if (were_nearest_neighbors) {
      return s1.s * s2.s;
    } else {
      return 0.0;
    }
  };
}

std::function<double(isingSpin)> isingBFace(unsigned L, double H) {
  return [L, H](const isingSpin& s) -> double {
    return H * s.s * phenomena[511 - s.x(1) * 512 / L][s.x(0) * 512 / L];
  };
}

std::function<double(isingSpin)> isingBMod(unsigned L, unsigned mod, double H) {
  return [mod, L, H](const isingSpin& s) -> double {
    return H * s.s * sin(s.x(0) * mod * 2 * M_PI / L);
  };
}

Gen<signed, D, TorusGroup<signed, D>, signed> isingGen(unsigned L) {
  std::vector<Vector<signed, 2>> torusVectors = torus_vecs<signed, 2>(L);
  std::vector<Matrix<signed, 2>> torusMatrices = torus_mats<signed, 2>();
  return [L, torusVectors,
          torusMatrices](Model<signed, D, TorusGroup<signed, D>, signed>& M,
                         Rng& r) -> Transformation<signed, D, TorusGroup<signed, D>, signed>* {
    Matrix<signed, 2> m;
    Vector<signed, 2> t;

    m = r.pick(torusMatrices);
    t(0) = r.uniform<signed>(0, L);
    t(1) = r.uniform<signed>(0, L);
    t = t - m * t;

    TorusGroup<signed, 2> g = TorusGroup<signed, 2>({(signed)L, t, m});

    Spin<signed, 2, signed>* ss = r.pick(M.s);
    Spin<signed, 2, signed> s2 = g.act(*ss);

    while (ss->x(0) == s2.x(0) && ss->x(1) == s2.x(1)) {
      ss = r.pick(M.s);
      s2 = g.act(*ss);
    }

    std::set<Spin<signed, 2, signed>*> s2s = M.dict.at(s2.x);

    if (s2s.empty()) {
      return new SpinFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss);
    } else {
      return new PairFlip<signed, 2, TorusGroup<signed, 2>, signed>(M, g, ss, *s2s.begin());
    }
  };
}

void isingPopulate(isingModel& m, unsigned L, double pop, double mag) {
  Rng rng;

  for (unsigned i = 0; i < L; i++) {
    for (unsigned j = 0; j < L; j++) {
      if (rng.uniform<double>(0, 1) < pop) {
        Spin<signed, D, signed>* ss = new Spin<signed, D, signed>();
        ss->x = {i, j};
        if (rng.uniform<double>(0, 1) < mag) {
          ss->s = 1;
        } else {
          ss->s = -1;
        }
        m.s.push_back(ss);
        m.dict.insert(ss);
      }
    }
  }
}