From 199b129c08242be6a2726aae3c9918ca2f2484f7 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 15 Jan 2021 16:29:59 +0100 Subject: More changes. --- dynamics.hpp | 18 +++++++++--------- langevin.cpp | 34 ++++++++++++++++++++++++++++++++-- p-spin.hpp | 9 --------- tensor.hpp | 49 ++++++++++++++++++++++++++++--------------------- 4 files changed, 69 insertions(+), 41 deletions(-) diff --git a/dynamics.hpp b/dynamics.hpp index 85eef71..b373659 100644 --- a/dynamics.hpp +++ b/dynamics.hpp @@ -92,26 +92,26 @@ Vector randomVector(unsigned N, Distribution d, Generator& r) { } template -std::tuple> langevin(const Tensor& J, const Vector& z0, double T, double γ, unsigned N, Distribution d, Generator& r) { +std::tuple> metropolis(const Tensor& J, const Vector& z0, + std::function&, const Vector&)>& energy, + double T, double γ, unsigned N, Distribution d, Generator& r) { Vector z = z0; - double W; - std::tie(W, std::ignore) = WdW(J, z); + double E = energy(J, z); std::uniform_real_distribution D(0, 1); for (unsigned i = 0; i < N; i++) { - Vector zNewTmp = z + randomVector(z.size(), d, r); + Vector zNewTmp = z + γ * randomVector(z.size(), d, r); Vector zNew = normalize(zNewTmp); - double WNew; - std::tie(WNew, std::ignore) = WdW(J, zNew); + double ENew = energy(J, zNew); - if (exp((W - WNew) / T) > D(r)) { + if (E - ENew > T * log(D(r))) { z = zNew; - W = WNew; + E = ENew; } } - return {W, z}; + return {E, z}; } diff --git a/langevin.cpp b/langevin.cpp index 870879b..5ca2d72 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -83,11 +83,41 @@ int main(int argc, char* argv[]) { ComplexVector z0 = normalize(randomVector(N, d, r.engine())); ComplexVector zSaddle = findSaddle(J, z0, ε); - ComplexVector zSaddlePrev = ComplexVector::Zero(N); ComplexVector z = zSaddle; + std::function energyNormGrad = [] + (const ComplexTensor& J, const ComplexVector& z) { + double W; + std::tie(W, std::ignore) = WdW(J, z); + return W; + }; + + double aGoal = 1e3; + + std::function energyInvA = [aGoal] + (const ComplexTensor& J, const ComplexVector& z) { + double a = z.squaredNorm(); + if (a > aGoal) { + return -aGoal; + } else { + return -z.squaredNorm(); + } + }; + + while (zSaddle.squaredNorm() < aGoal) { + std::tie(std::ignore, z) = metropolis(J, z, energyInvA, T, γ, 100, d, r.engine()); + try { + std::cerr << "Starting descent from " << z.squaredNorm() << "." << std::endl; + zSaddle = findSaddle(J, z, ε); + } catch (std::exception& e) { + } + std::cerr << "Current saddle is of size " << zSaddle.squaredNorm() << "." << std::endl; + } + + ComplexVector zSaddlePrev = ComplexVector::Zero(N); + while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close... - std::tie(std::ignore, z) = langevin(J, z, T, γ, M, d, r.engine()); + std::tie(std::ignore, z) = metropolis(J, z, energyNormGrad, T, γ, M, d, r.engine()); try { ComplexVector zSaddleNext = findSaddle(J, z, ε); if (Δ < (zSaddleNext - zSaddle).norm()) { // Ensure we are finding distinct saddles. diff --git a/p-spin.hpp b/p-spin.hpp index 4621db6..e5d4f94 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -5,15 +5,6 @@ #include "tensor.hpp" #include "factorial.hpp" -template -using Vector = Eigen::Matrix; - -template -using Matrix = Eigen::Matrix; - -template -using Tensor = Eigen::Tensor; - template std::tuple, Matrix> hamGradHess(const Tensor& J, const Vector& z) { Matrix Jz = contractDown(J, z); // Contracts J into p - 2 copies of z. diff --git a/tensor.hpp b/tensor.hpp index 21f2a89..6a5b2c2 100644 --- a/tensor.hpp +++ b/tensor.hpp @@ -5,43 +5,52 @@ #include +template +using Vector = Eigen::Matrix; + +template +using Matrix = Eigen::Matrix; + +template +using Tensor = Eigen::Tensor; + template -Eigen::Tensor initializeJHelper(unsigned N, std::index_sequence) { +Tensor initializeJHelper(unsigned N, std::index_sequence) { std::array Ns; std::fill_n(Ns.begin(), p, N); - return Eigen::Tensor(std::get(Ns)...); + return Tensor(std::get(Ns)...); } template -Eigen::Tensor initializeJ(unsigned N) { +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) { +void setJHelper(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) { +void setJ(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) { +Scalar getJHelper(const 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) { +Scalar getJ(const 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, +void iterateOverHelper(Tensor& J, + std::function&, std::array)>& f, unsigned l, std::array is) { if (l == 0) { f(J, is); @@ -61,17 +70,17 @@ void iterateOverHelper(Eigen::Tensor& J, } template -void iterateOver(Eigen::Tensor& J, std::function&, std::array)>& f) { +void iterateOver(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); +Tensor generateCouplings(unsigned N, Distribution d, Generator& r) { + Tensor J = initializeJ(N); - std::function&, std::array)> setRandom = - [&d, &r] (Eigen::Tensor& JJ, std::array ind) { + std::function&, std::array)> setRandom = + [&d, &r] (Tensor& JJ, std::array ind) { setJ(JJ, ind, d(r)); }; @@ -81,17 +90,15 @@ Eigen::Tensor generateCouplings(unsigned N, Distribution d, Generator } template -Eigen::Matrix -contractDown(const Eigen::Tensor& J, const Eigen::Matrix& z) { - return Eigen::Map>(J.data(), z.size(), z.size()); +Matrix contractDown(const Tensor& J, const Vector& z) { + return Eigen::Map>(J.data(), z.size(), z.size()); } const std::array, 1> ip00 = {Eigen::IndexPair(0, 0)}; template -Eigen::Matrix -contractDown(const Eigen::Tensor& J, const Eigen::Matrix& z) { - Eigen::Tensor zT = Eigen::TensorMap>(z.data(), {z.size()}); - Eigen::Tensor Jz = J.contract(zT, ip00); +Matrix contractDown(const Tensor& J, const Vector& z) { + Tensor zT = Eigen::TensorMap>(z.data(), {z.size()}); + Tensor Jz = J.contract(zT, ip00); return contractDown(Jz, z); } -- cgit v1.2.3-70-g09d2