From 9c3ac9c97abeca3ebcb11a1a2c8a6e3cb3791735 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 16 Feb 2021 16:42:53 +0100 Subject: Progress towards a working Stokes line finder. --- langevin.cpp | 92 ++++++++++++++++++++++++++++-------------- stokes.hpp | 129 ++++++++++++++++++++++++++++++++++++++++++----------------- 2 files changed, 154 insertions(+), 67 deletions(-) diff --git a/langevin.cpp b/langevin.cpp index 5ca2d72..5d9acc9 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -1,8 +1,10 @@ #include +#include #include "complex_normal.hpp" #include "p-spin.hpp" #include "dynamics.hpp" +#include "stokes.hpp" #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" @@ -26,7 +28,7 @@ int main(int argc, char* argv[]) { // simulation parameters double ε = 1e-4; - double εJ = 1e-2; + double εJ = 1e-5; double δ = 1e-2; // threshold for determining saddle double Δ = 1e-3; double γ = 1e-2; // step size @@ -92,28 +94,6 @@ int main(int argc, char* argv[]) { 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... @@ -131,7 +111,11 @@ int main(int argc, char* argv[]) { std::cerr << "Current saddles are " << (zSaddle - zSaddlePrev).norm() << " apart." << std::endl; } - std::cerr << "Found sufficiently nearby saddles, perturbing J." << std::endl; + Rope 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); @@ -141,20 +125,68 @@ int main(int argc, char* argv[]) { setJ(JJ, ind, Ji + dJ(r.engine())); }; - for (unsigned i = 0; i < n; i++) { - ComplexTensor Jp = J; + ComplexTensor Jp = J; + ComplexVector z1, z2; + Complex H1, H2; + double prevdist = 100; - iterateOver(Jp, perturbJ); + while (true) { + ComplexTensor Jpp = Jp; + iterateOver(Jpp, perturbJ); try { - ComplexVector zSaddleNew = findSaddle(Jp, zSaddle, ε); - ComplexVector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε); + z1 = findSaddle(Jpp, zSaddle, ε); + z2 = findSaddle(Jpp, zSaddlePrev, ε); + + double dist = (z1 - z2).norm(); + + std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jpp, z1); + std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jpp, z2); + + if (abs(imag(H1 - H2)) < prevdist && dist > 1e-2) { + Jp = Jpp; + prevdist = abs(imag(H1 - H2)); + + 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; + break; + } + } + + - std::cout << zSaddleNew.transpose() << " " << zSaddlePrevNew.transpose() << std::endl; } catch (std::exception& e) { std::cerr << "Couldn't find a saddle with new couplings, skipping." << std::endl; } } + Rope stokes(20, z1, z2); + + std::cout << stokes.cost(Jp) << std::endl; + + stokes.relax(Jp, 10000, 0.0001); + + std::cout << stokes.cost(Jp) << std::endl; + + Rope stokes2 = stokes.interpolate(); + + stokes2.relax(Jp, 10000, 0.001); + + std::cout << stokes2.cost(Jp) << std::endl; + + stokes2 = stokes2.interpolate(); + + stokes2.relax(Jp, 10000, 0.001); + + std::cout << stokes2.cost(Jp) << std::endl; + + stokes2 = stokes2.interpolate(); + + stokes2.relax(Jp, 10000, 0.001); + + std::cout << stokes2.cost(Jp) << std::endl; + + return 0; } diff --git a/stokes.hpp b/stokes.hpp index b60ae0e..43a69e3 100644 --- a/stokes.hpp +++ b/stokes.hpp @@ -1,47 +1,96 @@ -#include "stereographic.hpp" +#include "p-spin.hpp" +#include +#include "dynamics.hpp" template -class Rope { - private: - std::vector> z; +Vector zDot(const Vector& z, const Vector& dH) { + return -dH.conjugate() + (dH.dot(z) / z.squaredNorm()) * z.conjugate(); +} - template - Vector δz(const Tensor& J) const { - Vector dz(z.size()); - - for (unsigned i = 1; i < z.size() - 1; i++) { - Vector z12 = (z[i] + z[i - 1]) / 2.0; - Vector g; - std::tie(std::ignore, g, std::ignore) = stereographicHamGradHessHess(J, z12); +template +double segmentCost(const Vector& z, const Vector& dz, const Vector& dH) { + Vector zD = zDot(z, dH); + return 1.0 - pow(real(zD.dot(dz)), 2) / zD.squaredNorm() / dz.squaredNorm(); +} - Vector Δ = z[i] - z[i - 1]; +template +Vector variation(const Vector& z, const Vector& z´, const Vector& z´´, const Vector& dH, const Matrix& ddH) { + double z² = z.squaredNorm(); + double z´² = z´.squaredNorm(); + + Vector ż = zDot(z, dH); + double ż² = ż.squaredNorm(); + + double Reż·z´² = real(pow(ż.dot(z´), 2)); + + Matrix dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²; + Matrix dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z² + + (z.dot(dH) / z²) * (Matrix::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z²); + + Vector dLdz = - 0.5 * ( + ż.dot(z´) * (dżc * z´) + z´.dot(ż) * (dż * z´.conjugate()) + - (dż * ż.conjugate() + dżc * ż) * Reż·z´² / ż² + ) / (ż² * z´²); + + Vector ż´ = -(ddH * z´).conjugate() + ((ddH * z´).dot(z) / z²) * z.conjugate() + ( + dH.dot(z) * z´.conjugate() + dH.dot(z´) * z.conjugate() - (dH.dot(z) * (z´.dot(z) + z.dot(z´)) / z²) * z.conjugate() + ) / z²; + + Vector ddtdLdz´ = -0.5 * ( + ż´.dot(z´) * ż.conjugate() + ż.dot(z´´) * ż.conjugate() + ż.dot(z´) * ż´.conjugate() + - ( + Reż·z´² * z´´.conjugate() + real(2.0 * ż.dot(z´) * (ż.dot(z´´) + ż´.dot(z´))) * z´.conjugate() + - Reż·z´² * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() / z´² + ) / z´² + - ( + (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´² + ) * ( + ż.dot(z´) * ż.conjugate() - Reż·z´² / z´² * z´.conjugate() + ) + ) / (ż² * z´²); + + return dLdz - ddtdLdz´; +} - g.normalize(); - Δ.normalize(); +template +class Rope { + public: + std::vector> z; - if (abs(arg(g.dot(Δ))) < M_PI / 2) { - dz[i] = g - Δ; - } else { - dz[i] = -g - Δ; - } + Rope(unsigned N, const Vector& z1, const Vector& z2) : z(N + 2) { + for (unsigned i = 0; i < N + 2; i++) { + z[i] = z1 + (z2 - z1) * ((double)i / (N + 1.0)); + z[i] = normalize(z[i]); } + } + + Vector dz(unsigned i) const { + return z[i + 1] - z[i - 1]; + } - return dz; + Vector ddz(unsigned i) const { + return 4.0 * (z[i + 1] + z[i - 1] - 2.0 * z[i]); } template void perturb(const Tensor& J, double δ) { - z += δ * this->δz(J); + for (unsigned i = 1; i < z.size() - 1; i++) { + Vector dH; + Matrix ddH; + std::tie(std::ignore, dH, ddH) = hamGradHess(J, z[i]); + + Vector δz = variation(z[i], this->dz(i), this->ddz(i), dH, ddH); + + z[i] -= δ * δz.conjugate(); + z[i] = normalize(z[i]); + } } void spread() { - std::vector d(z.size() - 1); double l = 0; for (unsigned i= 0; i < z.size() - 1; i++) { - double Δ = (z[i + 1] - z[i]).norm(); - d[i] = Δ; - l += Δ; + l += (z[i + 1] - z[i]).norm(); } double a = 0; @@ -49,19 +98,16 @@ class Rope { for (unsigned i = 1; i < z.size() - 1; i++) { double b = i * l / (z.size() - 1); - while (b >= a) { - a += d[pos]; + + while (b > a) { pos++; + a += (z[pos] - z[pos - 1]).norm(); } - z[i] = z[pos] - (a - b) * (z[pos] - z[pos - 1]).normalized(); - } - } + Vector δz = z[pos] - z[pos - 1]; - public: - Rope(unsigned N, const Vector& z1, const Vector& z2) : z(N + 2) { - for (unsigned i = 0; i < N + 2; i++) { - z[i] = z1 + (z2 - z1) * ((double)i / (N + 1.0)); + z[i] = z[pos] - (a - b) / δz.norm() * δz; + z[i] = normalize(z[i]); } } @@ -73,7 +119,16 @@ class Rope { template double cost(const Tensor& J) const { - return this->δz(J).norm(); + double c = 0; + + for (unsigned i = 1; i < z.size() - 1; i++) { + Vector dH; + std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z[i]); + + c += segmentCost(z[i], this->dz(i), dH); + } + + return c; } template @@ -84,7 +139,7 @@ class Rope { } Rope interpolate() const { - Rope r(2 * z.size() - 1); + Rope r(2 * z.size() - 1, z[0], z[z.size() - 1]); for (unsigned i = 0; i < z.size(); i++) { r.z[2 * i] = z[i]; -- cgit v1.2.3-54-g00ecf