summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2022-01-04 17:42:30 +0100
committerJaron Kent-Dobias <jaron@kent-dobias.com>2022-01-04 17:42:30 +0100
commit94b6ed7d9c591ae2df0957a25e73bd2e5c69f6e1 (patch)
treedda5b4b9f7a8a08896260315f4f5bda1889ececc
parent874770326dc2063b56b0058f2f86ff47eccc6406 (diff)
downloadcode-94b6ed7d9c591ae2df0957a25e73bd2e5c69f6e1.tar.gz
code-94b6ed7d9c591ae2df0957a25e73bd2e5c69f6e1.tar.bz2
code-94b6ed7d9c591ae2df0957a25e73bd2e5c69f6e1.zip
Added option for spiking the interaction tensor to study minima below the threshold.HEADmaster
-rw-r--r--collectStokesData.hpp38
-rw-r--r--mixedStokes.cpp2
-rw-r--r--pureStokes.cpp8
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;