diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2022-01-04 17:42:30 +0100 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2022-01-04 17:42:30 +0100 |
commit | 94b6ed7d9c591ae2df0957a25e73bd2e5c69f6e1 (patch) | |
tree | dda5b4b9f7a8a08896260315f4f5bda1889ececc /collectStokesData.hpp | |
parent | 874770326dc2063b56b0058f2f86ff47eccc6406 (diff) | |
download | code-94b6ed7d9c591ae2df0957a25e73bd2e5c69f6e1.tar.gz code-94b6ed7d9c591ae2df0957a25e73bd2e5c69f6e1.tar.bz2 code-94b6ed7d9c591ae2df0957a25e73bd2e5c69f6e1.zip |
Added option for spiking the interaction tensor to study minima below the threshold.
Diffstat (limited to 'collectStokesData.hpp')
-rw-r--r-- | collectStokesData.hpp | 38 |
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); |