From 94b6ed7d9c591ae2df0957a25e73bd2e5c69f6e1 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 4 Jan 2022 17:42:30 +0100
Subject: Added option for spiking the interaction tensor to study minima below
 the threshold.

---
 collectStokesData.hpp | 38 +++++++++++++++++++++++++++++++++++++-
 mixedStokes.cpp       |  2 +-
 pureStokes.cpp        |  8 ++++++--
 3 files changed, 44 insertions(+), 4 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);
diff --git a/mixedStokes.cpp b/mixedStokes.cpp
index a437ff3..6a48f42 100644
--- a/mixedStokes.cpp
+++ b/mixedStokes.cpp
@@ -45,7 +45,7 @@ int main(int argc, char* argv[]) {
 
   for (unsigned i = 0; i < n; i++) {
     auto tag = std::chrono::high_resolution_clock::now();
-    collectStokesData<2, 4>(std::to_string(tag.time_since_epoch().count()), N, r.engine(), ε, δz₀, useMinima, 1.0, 0.01);
+    collectStokesData<2, 4>(std::to_string(tag.time_since_epoch().count()), N, r.engine(), ε, δz₀, useMinima, false, 1.0, 0.01);
   }
 
   return 0;
diff --git a/pureStokes.cpp b/pureStokes.cpp
index 448d5ef..b1ed1dc 100644
--- a/pureStokes.cpp
+++ b/pureStokes.cpp
@@ -16,10 +16,11 @@ int main(int argc, char* argv[]) {
   Real δz₀ = 1e-3;
   unsigned n = 10;
   bool useMinima = false;
+  bool useSpike = false;
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:e:d:n:m")) != -1) {
+  while ((opt = getopt(argc, argv, "N:e:d:n:ms")) != -1) {
     switch (opt) {
     case 'N':
       N = (unsigned)atof(optarg);
@@ -36,6 +37,9 @@ int main(int argc, char* argv[]) {
     case 'm':
       useMinima = true;
       break;
+    case 's':
+      useSpike = true;
+      break;
     default:
       exit(1);
     }
@@ -45,7 +49,7 @@ int main(int argc, char* argv[]) {
 
   for (unsigned i = 0; i < n; i++) {
     auto tag = std::chrono::high_resolution_clock::now();
-    collectStokesData<3>(std::to_string(tag.time_since_epoch().count()), N, r.engine(), ε, δz₀, useMinima, 1.0);
+    collectStokesData<3>(std::to_string(tag.time_since_epoch().count()), N, r.engine(), ε, δz₀, useMinima, useSpike, 1.0);
   }
 
   return 0;
-- 
cgit v1.2.3-54-g00ecf