From e6b4d83097f4442edf2290236e1b092724d8fd74 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 12 Jan 2021 17:39:27 +0100 Subject: Changed code to now find nearby saddles and perturb J. --- langevin.cpp | 49 ++++++++++++++++++++++++++++++++++------------ tensor.hpp | 64 +++++++++++++++++++++++++++++++++++++++++++++++------------- 2 files changed, 86 insertions(+), 27 deletions(-) diff --git a/langevin.cpp b/langevin.cpp index cf61b85..acac6c7 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -1,6 +1,7 @@ #include #include +#include #include "complex_normal.hpp" #include "p-spin.hpp" @@ -8,6 +9,7 @@ #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" +#include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h" using Rng = randutils::random_generator; @@ -129,7 +131,9 @@ int main(int argc, char* argv[]) { // simulation parameters double ε = 1e-4; + double εJ = 5e-2; double δ = 1e-2; // threshold for determining saddle + double Δ = 1e-3; double γ = 1e-2; // step size unsigned t = 1000; // number of Langevin steps unsigned M = 100; @@ -182,24 +186,43 @@ int main(int argc, char* argv[]) { Vector z0 = normalize(randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine())); Vector zSaddle = findSaddle(J, z0, ε); - - double W; + Vector zSaddlePrev = Vector::Zero(N); Vector z = zSaddle; - for (unsigned i = 0; i < n; i++) { - std::tie(W, z) = langevin(J, z, T, γ, M, r); + + while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close... + std::tie(std::ignore, z) = langevin(J, z, T, γ, M, r); try { - Vector zNewSaddle = findSaddle(J, z, ε); - Scalar H; - Matrix ddH; - std::tie(H, std::ignore, ddH) = hamGradHess(J, zNewSaddle); - Eigen::SelfAdjointEigenSolver es(ddH.adjoint() * ddH); - std::cout << N << "\t" << M * (i+1) << "\t" << H << "\t" - << zNewSaddle.transpose() << "\t" << es.eigenvalues().transpose() - << std::endl; - std::cerr << M * (i+1) << " steps taken to move " << (zNewSaddle - zSaddle).norm() << ", saddle information saved." << std::endl; + Vector zSaddleNext = findSaddle(J, z, ε); + if (Δ < (zSaddleNext - zSaddle).norm()) { // Ensure we are finding distinct saddles. + zSaddlePrev = zSaddle; + zSaddle = zSaddleNext; + } } catch (std::exception& e) { std::cerr << "Could not find a saddle: " << e.what() << std::endl; } + + std::cerr << "Current saddles are " << (zSaddle - zSaddlePrev).norm() << " apart." << std::endl; + } + + std::cerr << "Found sufficiently nearby saddles, perturbing J." << std::endl; + + complex_normal_distribution<> dJ(0, εJ * σ, 0); + + std::function)> perturbJ = + [&dJ, &r] (Tensor& JJ, std::array ind) { + Scalar Ji = getJ(JJ, ind); + setJ(JJ, ind, Ji + dJ(r.engine())); + }; + + for (unsigned i = 0; i < n; i++) { + Tensor Jp = J; + + iterateOver(Jp, perturbJ); + + Vector zSaddleNew = findSaddle(Jp, zSaddle, ε); + Vector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε); + + std::cout << (zSaddleNew - zSaddlePrevNew).norm() << std::endl; } return 0; diff --git a/tensor.hpp b/tensor.hpp index 528a51a..05fec36 100644 --- a/tensor.hpp +++ b/tensor.hpp @@ -1,27 +1,53 @@ #pragma once #include +#include #include +#include #include "factorial.hpp" template -Eigen::Tensor initializeJ(unsigned N, std::index_sequence) { +Eigen::Tensor initializeJHelper(unsigned N, std::index_sequence) { std::array Ns; std::fill_n(Ns.begin(), p, N); return Eigen::Tensor(std::get(Ns)...); } -template -void populateCouplings(Eigen::Tensor& J, unsigned N, unsigned l, - std::array is, Distribution d, Generator& r, - std::index_sequence ii) { +template +Eigen::Tensor initializeJ(unsigned N) { + return initializeJHelper(N, std::make_index_sequence

()); +} + +template +void setJHelper(Eigen::Tensor& J, const std::array& ind, Scalar val, std::index_sequence) { + J(std::get(ind)...) = val; +} + +template +void setJ(Eigen::Tensor& J, std::array ind, Scalar val) { + do { + setJHelper(J, ind, val, std::make_index_sequence

()); + } while (std::next_permutation(ind.begin(), ind.end())); +} + +template +Scalar getJHelper(const Eigen::Tensor& J, const std::array& ind, std::index_sequence) { + return J(std::get(ind)...); +} + +template +Scalar getJ(const Eigen::Tensor& J, const std::array& ind) { + return getJHelper(J, ind, std::make_index_sequence

()); +} + +template +void iterateOverHelper(Eigen::Tensor& J, + std::function&, std::array)>& f, + unsigned l, std::array is) { if (l == 0) { - Scalar z = d(r); - do { - J(std::get(is)...) = z; - } while (std::next_permutation(is.begin(), is.end())); + f(J, is); } else { unsigned iMin; if (l == p) { @@ -29,20 +55,30 @@ void populateCouplings(Eigen::Tensor& J, unsigned N, unsigned l, } else { iMin = is[p - l - 1]; } - for (unsigned i = iMin; i < N; i++) { + for (unsigned i = iMin; i < J.dimension(p - 1); i++) { std::array js = is; js[p - l] = i; - populateCouplings(J, N, l - 1, js, d, r, ii); + iterateOverHelper(J, f, l - 1, js); } } } +template +void iterateOver(Eigen::Tensor& J, std::function&, std::array)>& f) { + std::array is; + iterateOverHelper(J, f, p, is); +} + template Eigen::Tensor generateCouplings(unsigned N, Distribution d, Generator& r) { - Eigen::Tensor J = initializeJ(N, std::make_index_sequence

()); + Eigen::Tensor J = initializeJ(N); - std::array is; - populateCouplings(J, N, p, is, d, r, std::make_index_sequence

()); + std::function&, std::array)> setRandom = + [&d, &r] (Eigen::Tensor& JJ, std::array ind) { + setJ(JJ, ind, d(r)); + }; + + iterateOver(J, setRandom); return J; } -- cgit v1.2.3-70-g09d2