1 files changed, 37 insertions, 1 deletions
diff --git a/collectStokesData.hpp b/collectStokesData.hpp
index 54062f5..ca52386 100644
--- a/collectStokesData.hpp
+++ b/collectStokesData.hpp
@@ -4,8 +4,32 @@
 
 using Complex = std::complex<Real>;
 
+template<int p>
+Tensor<Real, p> makeSpikeHelper(const Vector<Real>& x) {
+  Eigen::array<Eigen::IndexPair<int>, 0> ei = {};
+  Tensor<Real, 1> xT = Eigen::TensorMap<Tensor<const Real, 1>>(x.data(), {x.size()});
+
+  return xT.contract(makeSpikeHelper<p - 1>(x), ei);
+}
+
+template <>
+Tensor<Real, 1> makeSpikeHelper(const Vector<Real>& x) {
+  Tensor<Real, 1> xT = Eigen::TensorMap<Tensor<const Real, 1>>(x.data(), {x.size()});
+  return xT;
+}
+
+template<int p, int ...ps>
+Tensor<Real, p> makeSpike(const Vector<Real>& x) {
+  return makeSpikeHelper<p>(x);
+}
+
+template <int p, int ...ps>
+int getFirstP() {
+  return p;
+}
+
 template<int ...ps, class Generator, typename... T>
-void collectStokesData(std::string tag, unsigned N, Generator& r, Real ε, Real δz₀, bool minimum, T... μs) {
+void collectStokesData(std::string tag, unsigned N, Generator& r, Real ε, Real δz₀, bool minimum, bool spike, T... μs) {
   unsigned nIts = 5;
   unsigned nGs = 4;
   unsigned coDim = 16;
@@ -24,6 +48,18 @@ void collectStokesData(std::string tag, unsigned N, Generator& r, Real ε, Real
   }
 
   Real Hx;
+  std::tie(Hx, std::ignore, std::ignore, std::ignore) = ReM.hamGradHess(xMin);
+
+  if (spike) {
+    const unsigned p = getFirstP<ps...>();
+    Real εgs = sqrt((p - 1) * log(p - 1) / (p - 2));
+    Real εth = sqrt((p - 1) * 2.0 / p);
+
+    Real εTarget = std::uniform_real_distribution<Real>(εth, εgs)(r);
+
+    std::get<0>(ReM.Js) -= factorial(p) * (εTarget * N + Hx) * makeSpike<ps...>(xMin) / pow(N, p);
+  }
+
   Vector<Real> dHx;
   Matrix<Real> ddHx;
   std::tie(Hx, dHx, ddHx, std::ignore) = ReM.hamGradHess(xMin);