#pragma once

#include <array>
#include <list>
#include <vector>

#include "matrix.hpp"
#include "vector.hpp"
#include "euclidean.hpp"
#include "transformation.hpp"

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 <class U, int D> std::vector<Vector<U, D>> torus_vecs(U L) {
  std::vector<Vector<U, D>> vecs;
  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_front(); // don't want the identity matrix!

  for (std::array<double, D> sequence : sequences) {
    Vector<U, D> v;
    for (unsigned i = 0; i < D; i++) {
      if (sequence[i] == 1) {
        v(i) = 0;
      } else {
        v(i) = L / 2;
      }
    }

    vecs.push_back(v);
  }

  return vecs;
}

template <class U, int D> std::vector<Matrix<U, D>> torus_mats() {
  std::vector<Matrix<U, 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<U, 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<U, 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 {
              m(k, l) = 0;
            }
          }
        }
        mats.push_back(m);
      }
    }
  }

  return mats;
}

template <int D, class S>
Gen<double, D, TorusGroup<double, D>, S> uniformGenTorus(double L) {
  std::vector<Vector<double, D>> torusVectors = torus_vecs<double, D>(L);
  std::vector<Matrix<double, D>> torusMatrices = torus_mats<double, D>();
  return [L, torusVectors,
          torusMatrices](Model<double, D, TorusGroup<double, D>, S>& M,
                         Rng& r) -> Transformation<double, D, TorusGroup<double, D>, S>* {
    Matrix<double, D> m = r.pick(torusMatrices);

    Vector<double, D> t;
    for (unsigned i = 0; i < D; i++) {
      t(i) = r.uniform<double>(0, L);
    }

    TorusGroup<double, D> g = TorusGroup<double, D>({(double)L, t - m * t, m});

    Spin<double, D, S>* ss = r.pick(M.s);

    return new SpinFlip<double, D, TorusGroup<double, D>, S>(M, g, ss);
  };
}