diff options
Diffstat (limited to 'langevin.cpp')
| -rw-r--r-- | langevin.cpp | 77 |
1 files changed, 48 insertions, 29 deletions
diff --git a/langevin.cpp b/langevin.cpp index 5d9acc9..e1464ec 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -1,5 +1,7 @@ #include <getopt.h> #include <stereographic.hpp> +#include <unordered_map> +#include <list> #include "complex_normal.hpp" #include "p-spin.hpp" @@ -19,6 +21,20 @@ using ComplexTensor = Tensor<Complex, p>; using Rng = randutils::random_generator<pcg32>; +std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_map<uint64_t, ComplexVector>& saddles, double δ) { + std::list<std::array<ComplexVector, 2>> pairs; + + for (auto it1 = saddles.begin(); it1 != saddles.end(); it1++) { + for (auto it2 = std::next(it1); it2 != saddles.end(); it2++) { + if ((it1->second - it2->second).norm() < δ) { + pairs.push_back({it1->second, it2->second}); + } + } + } + + return pairs; +} + int main(int argc, char* argv[]) { // model parameters unsigned N = 10; // number of spins @@ -82,10 +98,7 @@ int main(int argc, char* argv[]) { complex_normal_distribution<> d(0, 1, 0); ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine()); - ComplexVector z0 = normalize(randomVector<Complex>(N, d, r.engine())); - - ComplexVector zSaddle = findSaddle(J, z0, ε); - ComplexVector z = zSaddle; + ComplexVector z = normalize(randomVector<Complex>(N, d, r.engine())); std::function<double(const ComplexTensor&, const ComplexVector&)> energyNormGrad = [] (const ComplexTensor& J, const ComplexVector& z) { @@ -94,49 +107,55 @@ int main(int argc, char* argv[]) { return W; }; - ComplexVector zSaddlePrev = ComplexVector::Zero(N); + std::unordered_map<uint64_t, ComplexVector> saddles; + std::list<std::array<ComplexVector, 2>> nearbySaddles; - while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close... + while (true) { // Until we find two saddles sufficiently close... 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. - zSaddlePrev = zSaddle; - zSaddle = zSaddleNext; + uint64_t saddleKey = round(1e2 * real(zSaddleNext(0))); + auto got = saddles.find(saddleKey); + + if (got == saddles.end()) { + saddles[saddleKey] = zSaddleNext; + nearbySaddles = saddlesCloserThan(saddles, δ); + if (nearbySaddles.size() > 0) { + break; + } + std::cerr << "Found " << saddles.size() << " distinct saddles." << std::endl; } } catch (std::exception& e) { - std::cerr << "Could not find a saddle: " << e.what() << std::endl; +// std::cerr << "Could not find a saddle: " << e.what() << std::endl; } - std::cerr << "Current saddles are " << (zSaddle - zSaddlePrev).norm() << " apart." << std::endl; } - Rope<Complex> stokest(20, zSaddle, zSaddlePrev); - - stokest.relax(J, 10000, 0.0001); - std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl; complex_normal_distribution<> dJ(0, εJ * σ, 0); std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ = - [&dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) { + [&εJ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) { Complex Ji = getJ<Complex, p>(JJ, ind); - setJ<Complex, p>(JJ, ind, Ji + dJ(r.engine())); + setJ<Complex, p>(JJ, ind, Ji + εJ * dJ(r.engine())); }; ComplexTensor Jp = J; - ComplexVector z1, z2; Complex H1, H2; - double prevdist = 100; + ComplexVector z1, z2; + std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, nearbySaddles.front()[0]); + std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, nearbySaddles.front()[1]); + double prevdist = abs(imag(H1-H2)); + εJ = 1e4 * prevdist; while (true) { ComplexTensor Jpp = Jp; iterateOver<Complex, p>(Jpp, perturbJ); try { - z1 = findSaddle(Jpp, zSaddle, ε); - z2 = findSaddle(Jpp, zSaddlePrev, ε); + z1 = findSaddle(Jpp, nearbySaddles.front()[0], ε); + z2 = findSaddle(Jpp, nearbySaddles.front()[1], ε); double dist = (z1 - z2).norm(); @@ -146,10 +165,10 @@ int main(int argc, char* argv[]) { if (abs(imag(H1 - H2)) < prevdist && dist > 1e-2) { Jp = Jpp; prevdist = abs(imag(H1 - H2)); + εJ = 1e4 * prevdist; - std::cout << abs(imag(H1 - H2)) << " " << dist << std::endl; - if (abs(imag(H1 - H2)) < 1e-8 && dist > 1e-2) { - std::cout << "Found distinct critical points with sufficiently similar Im H." << std::endl; + if (abs(imag(H1 - H2)) < 1e-10 && dist > 1e-2) { + std::cerr << "Found distinct critical points with sufficiently similar Im H." << std::endl; break; } } @@ -161,29 +180,29 @@ int main(int argc, char* argv[]) { } } - Rope<Complex> stokes(20, z1, z2); + Rope<Complex> stokes(10, z1, z2); std::cout << stokes.cost(Jp) << std::endl; - stokes.relax(Jp, 10000, 0.0001); + stokes.relax(Jp, 10000, 1e-4); std::cout << stokes.cost(Jp) << std::endl; Rope<Complex> stokes2 = stokes.interpolate(); - stokes2.relax(Jp, 10000, 0.001); + stokes2.relax(Jp, 10000, 1e-4); std::cout << stokes2.cost(Jp) << std::endl; stokes2 = stokes2.interpolate(); - stokes2.relax(Jp, 10000, 0.001); + stokes2.relax(Jp, 10000, 1e-4); std::cout << stokes2.cost(Jp) << std::endl; stokes2 = stokes2.interpolate(); - stokes2.relax(Jp, 10000, 0.001); + stokes2.relax(Jp, 10000, 1e-4); std::cout << stokes2.cost(Jp) << std::endl; |