diff options
| -rw-r--r-- | collectStokesData.hpp | 38 | ||||
| -rw-r--r-- | mixedStokes.cpp | 2 | ||||
| -rw-r--r-- | 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; |