#pragma once

#include "hadamard_mcmc.hpp"
#include <array>

Orthogonal walsh(unsigned k) {
  if (k == 0) {
    return Orthogonal(1);
  } else {
    Orthogonal s = walsh(k - 1);
    Orthogonal t = Orthogonal(2 * s.size());

    for (unsigned i = 0; i < s.size(); i++) {
      for (unsigned j = 0; j < s.size(); j++) {
        double sij = s(i, j);

        t(i, j) = sij;
        t(s.size() + i, j) = sij;
        t(i, s.size() + j) = sij;
        t(s.size() + i, s.size() + j) = -sij;
      }
    }

    return t;
  }
}

Orthogonal m4 = walsh(2);
Orthogonal m8 = walsh(3);
Orthogonal m16 = walsh(4);
Orthogonal m32 = walsh(5);

Orthogonal m12(
    {
       1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
       1, 1,-1, 1,-1,-1,-1, 1, 1, 1,-1, 1,
       1, 1, 1,-1, 1,-1,-1,-1, 1, 1, 1,-1,
       1,-1, 1, 1,-1, 1,-1,-1,-1, 1, 1, 1,
       1, 1,-1, 1, 1,-1, 1,-1,-1,-1, 1, 1,
       1, 1, 1,-1, 1, 1,-1, 1,-1,-1,-1, 1,
       1, 1, 1, 1,-1, 1, 1,-1, 1,-1,-1,-1,
       1,-1, 1, 1, 1,-1, 1, 1,-1, 1,-1,-1,
       1,-1,-1, 1, 1, 1,-1, 1, 1,-1, 1,-1,
       1,-1,-1,-1, 1, 1, 1,-1, 1, 1,-1, 1,
       1, 1,-1,-1,-1, 1, 1, 1,-1, 1, 1,-1,
       1,-1, 1,-1,-1,-1, 1, 1, 1,-1, 1, 1
    }
    );

Orthogonal m20(
    {
       1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
       1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1,
       1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1,
       1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,
       1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,
       1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1,
       1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1,
       1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1,
       1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1,
       1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,
       1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1,
       1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1, 1,
       1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,-1,
       1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,-1,
       1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1,
       1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,
       1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1, 1,
       1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1,
       1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1,
       1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1, 1
    }
    );

Orthogonal m24(
    {
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
       1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,
       1, 1,-1,-1,-1,-1, 1, 1, 1, 1,-1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1,-1,-1,-1,
       1, 1, 1, 1,-1,-1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1,
       1, 1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1, 1, 1,-1,-1,
       1, 1,-1,-1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1,
       1, 1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1,-1,-1, 1, 1,
       1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1,
       1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1,-1,-1, 1, 1, 1, 1,
       1, 1,-1,-1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1, 1, 1,-1,-1,-1,-1, 1, 1,
       1, 1,-1,-1, 1, 1,-1,-1,-1,-1, 1, 1, 1, 1,-1,-1,-1,-1, 1, 1,-1,-1, 1, 1,
       1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,
       1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,
       1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1,
       1,-1,-1, 1,-1, 1, 1,-1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1, 1,-1,-1, 1,-1, 1,
       1,-1, 1,-1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1,
       1,-1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1,
       1,-1,-1, 1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, 1,-1,-1, 1,
       1,-1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1,-1, 1, 1,-1,
       1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1,
       1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1, 1,-1,
       1,-1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1,-1, 1, 1,-1,
       1,-1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1
    }
    );

Orthogonal m28(
    {
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1,-1,-1, 1, 1, 1,-1,-1,-1, 1, 1, 1,-1,-1,-1, 1, 1, 1,-1,-1,-1, 1, 1, 1, 1,-1,-1,-1,
       1,-1,-1,-1,-1,-1, 1, 1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1, 1,-1, 1, 1, 1, 1,-1,-1,-1,
       1,-1,-1,-1,-1,-1, 1, 1, 1,-1, 1, 1,-1,-1, 1,-1, 1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,
       1,-1,-1, 1, 1,-1, 1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1,-1,-1,-1, 1, 1,-1, 1,-1, 1, 1,-1,
       1,-1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1, 1,-1,-1,-1, 1,-1, 1, 1,-1, 1,-1, 1, 1, 1,-1,
       1,-1,-1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1, 1,-1,-1,-1, 1, 1, 1,-1,-1, 1, 1, 1, 1,-1,
      -1, 1,-1, 1, 1,-1,-1,-1, 1,-1, 1,-1, 1,-1, 1,-1,-1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,
      -1, 1,-1,-1,-1, 1,-1, 1, 1,-1, 1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1, 1,-1, 1,-1, 1,-1,-1,
      -1, 1,-1,-1,-1, 1, 1,-1, 1, 1,-1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1,-1,
      -1, 1,-1, 1,-1, 1, 1,-1,-1,-1,-1, 1,-1, 1, 1, 1,-1, 1, 1, 1,-1, 1,-1,-1, 1,-1, 1,-1,
      -1, 1,-1,-1, 1,-1, 1, 1,-1, 1, 1,-1,-1, 1,-1,-1, 1, 1, 1,-1, 1,-1,-1, 1, 1,-1, 1,-1,
      -1, 1,-1, 1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1, 1, 1,-1,-1, 1,-1,-1, 1, 1,-1,
      -1,-1, 1,-1, 1, 1,-1, 1,-1,-1,-1, 1, 1, 1, 1,-1, 1,-1, 1,-1, 1, 1, 1,-1, 1,-1,-1,-1,
      -1,-1, 1, 1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1, 1,-1,-1,
      -1,-1, 1,-1, 1,-1, 1, 1,-1,-1, 1, 1, 1,-1, 1, 1,-1, 1,-1, 1,-1,-1,-1, 1, 1, 1,-1,-1,
      -1,-1, 1, 1,-1,-1,-1, 1, 1, 1,-1, 1,-1, 1, 1, 1,-1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1,-1,
      -1,-1, 1, 1,-1,-1, 1,-1, 1, 1, 1, 1, 1,-1,-1,-1, 1,-1,-1, 1, 1, 1,-1,-1, 1,-1, 1,-1,
      -1,-1, 1,-1, 1, 1,-1,-1, 1, 1, 1,-1, 1, 1,-1, 1,-1, 1, 1,-1,-1, 1,-1,-1,-1, 1, 1,-1,
       1, 1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1,-1,-1, 1,-1, 1, 1, 1, 1,-1,-1, 1,-1,-1,
       1, 1, 1, 1,-1,-1,-1, 1,-1, 1, 1,-1,-1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1,-1, 1, 1,-1,-1,
       1, 1, 1, 1,-1,-1,-1,-1, 1,-1,-1, 1, 1, 1,-1,-1, 1, 1, 1,-1,-1,-1,-1, 1, 1, 1,-1,-1,
       1, 1, 1,-1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1,
       1, 1, 1,-1,-1, 1, 1,-1,-1,-1, 1, 1, 1,-1,-1, 1,-1,-1, 1,-1, 1,-1, 1, 1,-1,-1, 1,-1,
       1, 1, 1,-1, 1,-1,-1,-1, 1,-1, 1,-1,-1, 1, 1, 1, 1,-1,-1, 1,-1,-1, 1,-1, 1,-1, 1,-1,
       1, 1,-1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
       1,-1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,
      -1, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1,-1
    }
    );

std::vector<Orthogonal> hadamards = {m4, m8, m12, m16, m20, m24, m28, m32};