diff options
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); |