diff options
| -rw-r--r-- | langevin.cpp | 49 | ||||
| -rw-r--r-- | 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 <getopt.h> #include <eigen3/Eigen/LU> +#include <utility> #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<pcg32>; @@ -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<Matrix> 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<void(Tensor&, std::array<unsigned, p>)> perturbJ = + [&dJ, &r] (Tensor& JJ, std::array<unsigned, p> ind) { + Scalar Ji = getJ<Scalar, p>(JJ, ind); + setJ<Scalar, p>(JJ, ind, Ji + dJ(r.engine())); + }; + + for (unsigned i = 0; i < n; i++) { + Tensor Jp = J; + + iterateOver<Scalar, p>(Jp, perturbJ); + + Vector zSaddleNew = findSaddle(Jp, zSaddle, ε); + Vector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε); + + std::cout << (zSaddleNew - zSaddlePrevNew).norm() << std::endl; } return 0; @@ -1,27 +1,53 @@ #pragma once #include <array> +#include <functional> #include <eigen3/unsupported/Eigen/CXX11/Tensor> +#include <utility> #include "factorial.hpp" template <class Scalar, unsigned p, std::size_t... Indices> -Eigen::Tensor<Scalar, p> initializeJ(unsigned N, std::index_sequence<Indices...>) { +Eigen::Tensor<Scalar, p> initializeJHelper(unsigned N, std::index_sequence<Indices...>) { std::array<unsigned, p> Ns; std::fill_n(Ns.begin(), p, N); return Eigen::Tensor<Scalar, p>(std::get<Indices>(Ns)...); } -template <class Scalar, unsigned p, class Distribution, class Generator, std::size_t... Indices> -void populateCouplings(Eigen::Tensor<Scalar, p>& J, unsigned N, unsigned l, - std::array<unsigned, p> is, Distribution d, Generator& r, - std::index_sequence<Indices...> ii) { +template <class Scalar, unsigned p> +Eigen::Tensor<Scalar, p> initializeJ(unsigned N) { + return initializeJHelper<Scalar, p>(N, std::make_index_sequence<p>()); +} + +template <class Scalar, unsigned p, std::size_t... Indices> +void setJHelper(Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, Scalar val, std::index_sequence<Indices...>) { + J(std::get<Indices>(ind)...) = val; +} + +template <class Scalar, unsigned p> +void setJ(Eigen::Tensor<Scalar, p>& J, std::array<unsigned, p> ind, Scalar val) { + do { + setJHelper<Scalar, p>(J, ind, val, std::make_index_sequence<p>()); + } while (std::next_permutation(ind.begin(), ind.end())); +} + +template <class Scalar, unsigned p, std::size_t... Indices> +Scalar getJHelper(const Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind, std::index_sequence<Indices...>) { + return J(std::get<Indices>(ind)...); +} + +template <class Scalar, unsigned p> +Scalar getJ(const Eigen::Tensor<Scalar, p>& J, const std::array<unsigned, p>& ind) { + return getJHelper<Scalar, p>(J, ind, std::make_index_sequence<p>()); +} + +template <class Scalar, unsigned p> +void iterateOverHelper(Eigen::Tensor<Scalar, p>& J, + std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)>& f, + unsigned l, std::array<unsigned, p> is) { if (l == 0) { - Scalar z = d(r); - do { - J(std::get<Indices>(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<Scalar, p>& 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<unsigned, p> js = is; js[p - l] = i; - populateCouplings<Scalar, p>(J, N, l - 1, js, d, r, ii); + iterateOverHelper<Scalar, p>(J, f, l - 1, js); } } } +template <class Scalar, unsigned p> +void iterateOver(Eigen::Tensor<Scalar, p>& J, std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)>& f) { + std::array<unsigned, p> is; + iterateOverHelper<Scalar, p>(J, f, p, is); +} + template <class Scalar, unsigned p, class Distribution, class Generator> Eigen::Tensor<Scalar, p> generateCouplings(unsigned N, Distribution d, Generator& r) { - Eigen::Tensor<Scalar, p> J = initializeJ<Scalar, p>(N, std::make_index_sequence<p>()); + Eigen::Tensor<Scalar, p> J = initializeJ<Scalar, p>(N); - std::array<unsigned, p> is; - populateCouplings<Scalar, p>(J, N, p, is, d, r, std::make_index_sequence<p>()); + std::function<void(Eigen::Tensor<Scalar, p>&, std::array<unsigned, p>)> setRandom = + [&d, &r] (Eigen::Tensor<Scalar, p>& JJ, std::array<unsigned, p> ind) { + setJ<Scalar, p>(JJ, ind, d(r)); + }; + + iterateOver<Scalar, p>(J, setRandom); return J; } |