Changed code to now find nearby saddles and perturb J.

1 files changed, 50 insertions, 14 deletions

diff --git a/tensor.hpp b/tensor.hpp
index 528a51a..05fec36 100644
--- a/tensor.hpp
+++ b/tensor.hpp
@@ -1,27 +1,53 @@
 #pragma once
 
 #include <array>
+#include <functional>
 
 #include <eigen3/unsupported/Eigen/CXX11/Tensor>
+#include <utility>
 
 #include "factorial.hpp"
 
 template <class Scalar, unsigned p, std::size_t... Indices>
-Eigen::Tensor<Scalar, p> initializeJ(unsigned N, std::index_sequence<Indices...>) {
+Eigen::Tensor<Scalar, p> initializeJHelper(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, std::size_t... Indices>
-void populateCouplings(Eigen::Tensor<Scalar, p>& J, unsigned N, unsigned l,
-                       std::array<unsigned, p> is, Distribution d, Generator& r,
-                       std::index_sequence<Indices...> ii) {
+template <class Scalar, unsigned p>
+Eigen::Tensor<Scalar, p> initializeJ(unsigned N) {
+  return initializeJHelper<Scalar, p>(N, std::make_index_sequence<p>());
+}
+
+template <class Scalar, unsigned p, std::size_t... Indices>
+void setJHelper(Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, Scalar val, std::index_sequence<Indices...>) {
+  J(std::get<Indices>(ind)...) = val;
+}
+
+template <class Scalar, unsigned p>
+void setJ(Eigen::Tensor<Scalar, p>& J, std::array<unsigned, p> ind, Scalar val) {
+  do {
+    setJHelper<Scalar, p>(J, ind, val, std::make_index_sequence<p>());
+  } while (std::next_permutation(ind.begin(), ind.end()));
+}
+
+template <class Scalar, unsigned p, std::size_t... Indices>
+Scalar getJHelper(const Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, std::index_sequence<Indices...>) {
+  return J(std::get<Indices>(ind)...);
+}
+
+template <class Scalar, unsigned p>
+Scalar getJ(const Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind) {
+  return getJHelper<Scalar, p>(J, ind, std::make_index_sequence<p>());
+}
+
+template <class Scalar, unsigned p>
+void iterateOverHelper(Eigen::Tensor<Scalar, p>& J,
+    std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)>& f,
+    unsigned l, std::array<unsigned, p> is) {
   if (l == 0) {
-    Scalar z = d(r);
-    do {
-      J(std::get<Indices>(is)...) = z;
-    } while (std::next_permutation(is.begin(), is.end()));
+    f(J, is);
   } else {
     unsigned iMin;
     if (l == p) {
@@ -29,20 +55,30 @@ void populateCouplings(Eigen::Tensor<Scalar, p>& J, unsigned N, unsigned l,
     } else {
       iMin = is[p - l - 1];
     }
-    for (unsigned i = iMin; i < N; i++) {
+    for (unsigned i = iMin; i < J.dimension(p - 1); i++) {
       std::array<unsigned, p> js = is;
       js[p - l] = i;
-      populateCouplings<Scalar, p>(J, N, l - 1, js, d, r, ii);
+      iterateOverHelper<Scalar, p>(J, f, l - 1, js);
     }
   }
 }
 
+template <class Scalar, unsigned p>
+void iterateOver(Eigen::Tensor<Scalar, p>& J, std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)>& f) {
+  std::array<unsigned, p> is;
+  iterateOverHelper<Scalar, p>(J, f, p, is);
+}
+
 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>());
+  Eigen::Tensor<Scalar, p> J = initializeJ<Scalar, p>(N);
 
-  std::array<unsigned, p> is;
-  populateCouplings<Scalar, p>(J, N, p, is, d, r, std::make_index_sequence<p>());
+  std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)> setRandom =
+    [&d, &r] (Eigen::Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
+      setJ<Scalar, p>(JJ, ind, d(r));
+    };
+
+  iterateOver<Scalar, p>(J, setRandom);
 
   return J;
 }