Implemented planting.

1 files changed, 21 insertions, 0 deletions

diff --git a/tensor.hpp b/tensor.hpp
index 6a5b2c2..1d90c78 100644
--- a/tensor.hpp
+++ b/tensor.hpp
@@ -89,6 +89,27 @@ Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) {
   return J;
 }
 
+template <class Scalar, int p, class Distribution, class Generator>
+Tensor<Scalar, p> plantState(const Tensor<Scalar, p>& J, const Vector<Scalar>& z, double β) {
+  Tensor<Scalar, p> JPlanted = J;
+
+  std::function<void(Tensor<Scalar, p>&, std::array<unsigned, p>)> plant =
+    [&z, β] (Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) {
+      Scalar Ji = getJ<Scalar, p>(JJ, ind);
+      Scalar zzz = 1;
+
+      for (unsigned i : ind) {
+        zzz *= z(i);
+      }
+
+      setJ<Scalar, p>(JJ, ind, Ji - β * zzz / pow(zzz.size(), p));
+    };
+
+  iterateOver<Scalar, p>(JPlanted, plant);
+
+  return JPlanted;
+}
+
 template <class Scalar>
 Matrix<Scalar> contractDown(const Tensor<Scalar, 2>& J, const Vector<Scalar>& z) {
   return Eigen::Map<const Matrix<Scalar>>(J.data(), z.size(), z.size());