Lots of work, refactoring.

1 files changed, 65 insertions, 0 deletions

diff --git a/tensor.hpp b/tensor.hpp
new file mode 100644
index 0000000..8345d83
--- /dev/null
+++ b/tensor.hpp
@@ -0,0 +1,65 @@
+#pragma once
+
+#include <algorithm>
+#include <array>
+#include <utility>
+#include <eigen3/unsupported/Eigen/CXX11/Tensor>
+
+#include "factorial.hpp"
+
+template <class Scalar, unsigned p, std::size_t... Indices>
+void setJ(Scalar z, Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& is, std::index_sequence<Indices...>) {
+  J(std::get<Indices>(is)...) = z;
+}
+
+template <class Scalar, unsigned p, std::size_t... Indices>
+Eigen::Tensor<Scalar, p> initializeJ(unsigned N, std::index_sequence<Indices...>) {
+  std::array<unsigned, p> Ns;
+  std::fill_n(Ns.begin(), p, N);
+  return Eigen::Tensor<Scalar, p>(std::get<Indices>(Ns)...);
+}
+
+template <class Scalar, unsigned p, class Distribution, class Generator>
+void populateCouplings(Eigen::Tensor<Scalar, p>& J, unsigned N, unsigned l, std::array<unsigned, p> is, Distribution d, Generator& r) {
+  if (l == 0) {
+    Scalar z = d(r);
+    do {
+      setJ<Scalar, p>(z, J, is, std::make_index_sequence<p>());
+    } while (std::next_permutation(is.begin(), is.end()));
+  } else {
+    unsigned iMin;
+    if (l == p) {
+      iMin = 0;
+    } else {
+      iMin = is[p - l - 1];
+    }
+    for (unsigned i = iMin; i < N; i++) {
+      std::array<unsigned, p> js = is;
+      js[p - l] = i;
+      populateCouplings<Scalar, p, Distribution, Generator>(J, N, l - 1, js, d, r);
+    }
+  }
+}
+
+template <class Scalar, unsigned p, class Distribution, class Generator>
+Eigen::Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
+  Eigen::Tensor<Scalar, p> J = initializeJ<Scalar, p>(N, std::make_index_sequence<p>());
+
+  std::array<unsigned, p> is;
+  populateCouplings<Scalar, p, Distribution, Generator>(J, N, p, is, d, r);
+
+  return J;
+}
+
+template <class Scalar>
+Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> contractDown(const Eigen::Tensor<Scalar, 2>& J, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& z) {
+  return Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>(J.data(), z.size(), z.size());
+}
+
+template <class Scalar, int r>
+Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> contractDown(const Eigen::Tensor<Scalar, r>& J, const Eigen::Matrix<Scalar, Eigen::Dynamic, 1>& z) {
+  Eigen::Tensor<Scalar, 1> zT = Eigen::TensorMap<Eigen::Tensor<const Scalar, 1>>(z.data(), {z.size()});
+  std::array<Eigen::IndexPair<int>, 1> ip = {Eigen::IndexPair<int>(0,0)};
+  Eigen::Tensor<Scalar, r - 1> Jz = J.contract(zT, ip);
+  return contractDown(Jz, z);
+}