diff options
Diffstat (limited to 'langevin.cpp')
| -rw-r--r-- | langevin.cpp | 129 |
1 files changed, 41 insertions, 88 deletions
diff --git a/langevin.cpp b/langevin.cpp index dc7c5cf..fe26332 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -1,8 +1,10 @@ #include <getopt.h> +#include <limits> #include <stereographic.hpp> #include <unordered_map> #include <list> +#include "Eigen/src/Eigenvalues/ComplexEigenSolver.h" #include "complex_normal.hpp" #include "p-spin.hpp" #include "dynamics.hpp" @@ -84,11 +86,10 @@ int main(int argc, char* argv[]) { Real T = 1; // temperature Real Rκ = 0; // real part of distribution parameter Real Iκ = 0; // imaginary part of distribution parameter - // simulation parameters Real ε = 1e-15; Real εJ = 1e-5; - Real δ = 1e-2; // threshold for determining saddle + Real δ = 1; // threshold for determining saddle Real Δ = 1e-3; Real γ = 1e-1; // step size unsigned t = 1000; // number of Langevin steps @@ -140,7 +141,7 @@ int main(int argc, char* argv[]) { complex_normal_distribution<Real> d(0, 1, 0); - ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine()); + ComplexTensor J = generateCouplings<Complex, p>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine()); std::function<Real(const ComplexTensor&, const ComplexVector&)> energyNormGrad = [] (const ComplexTensor& J, const ComplexVector& z) { @@ -149,106 +150,58 @@ int main(int argc, char* argv[]) { return W; }; - - std::function<Real(const ComplexTensor&, const ComplexVector&)> energyThres = [N, κ] - (const ComplexTensor& J, const ComplexVector& z) { - Real a = z.squaredNorm() / (Real)N; - /* - Complex ε = getHamiltonian(J, z) / (Real)N; - Real t = norm(ε) / getThresholdEnergyDensity(p, κ, ε, a); - */ - - Real amin = 1.1547; - Real amax = 1.46806; - - /* - Real tmin = 1; - Real tmax = 1.17514; - */ - - Real E = 0; - - if (a > amax) { - E += pow(a - amax, 2); - } else if (a < amin) { - E += pow(a - amin, 2); - } - - /* - if (t > tmax) { - E += pow(t - tmax, 2); - } else if (t < tmin) { - E += pow(t - tmin, 2); - } - */ - - return E; - }; - - Real largestThreshold = 0; - ComplexVector saddlePastThreshold; - bool foundSaddle = false; - -#pragma omp parallel default(none) shared(largestThreshold, saddlePastThreshold, foundSaddle, M, J, energyThres, d, N, γ, ε, κ) - { - Rng rPrivate; - ComplexVector z = normalize(randomVector<Complex>(N, d, rPrivate.engine())); - - while (largestThreshold < 1) { // Until we find a saddle past the threshold... - std::tie(std::ignore, z) = metropolis(J, z, energyThres, (Real)0.1, γ, M, d, rPrivate.engine()); - try { - ComplexVector latestSaddle = findSaddle(J, z, ε); - Real threshold = getProportionOfThreshold(κ, J, latestSaddle); - -#pragma omp critical - { - if (threshold > largestThreshold + 1e-6) { - largestThreshold = threshold; - saddlePastThreshold = latestSaddle; - std::cerr << "Found saddle with threshold porportion " << largestThreshold << std::endl;; - } - } - } catch (std::exception& e) {} - } - } - - std::cerr << "Found saddle with energy beyond threshold." << std::endl; + ComplexVector zSaddle = randomSaddle(J, d, r, ε); + std::cerr << "Found saddle." << std::endl; ComplexVector zSaddleNext; - Real closestSaddle = std::numeric_limits<Real>::infinity(); - ComplexVector z = saddlePastThreshold; - - while (closestSaddle > δ) { // Until we find two saddles sufficiently close... - std::tie(std::ignore, z) = metropolis(J, saddlePastThreshold, energyNormGrad, T, γ, 1e4, d, r.engine()); + bool foundSaddle = false; + while (!foundSaddle) { + ComplexVector z0 = normalize(zSaddle + δ * randomVector<Complex>(N, d, r.engine())); try { - zSaddleNext = findSaddle(J, z, ε); - Real saddleDistance = (zSaddleNext - saddlePastThreshold).norm(); - if (saddleDistance < closestSaddle && saddleDistance > 1e-2) { - closestSaddle = saddleDistance; - std::cerr << "Nearby saddle: found saddle at distance " << saddleDistance << std::endl; + zSaddleNext = findSaddle(J, z0, ε); + Real saddleDistance = (zSaddleNext - zSaddle).norm(); + if (saddleDistance / N > 1e-2) { + foundSaddle = true; } } catch (std::exception& e) {} } - std::cerr << "Found sufficiently nearby saddles, perturbing J to equalize Im H." << std::endl; + auto [H1, dH1, ddH1] = hamGradHess(J, zSaddle); + auto [H2, dH2, ddH2] = hamGradHess(J, zSaddleNext); - ComplexVector z1 = saddlePastThreshold; - ComplexVector z2 = zSaddleNext; + Eigen::ComplexEigenSolver<ComplexMatrix> ces; + ces.compute(ddH1); - std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r); + Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real()); + std::cerr << (zSaddle - zSaddleNext).norm() / (Real)N << " " << φ << " " << H1 * exp(Complex(0, φ)) << " " << H2 * exp(Complex(0, φ)) << std::endl; + Real smallestNorm = std::numeric_limits<Real>::infinity(); + for (Complex z : ces.eigenvalues()) { + if (norm(z) < smallestNorm) { + smallestNorm = norm(z); + } + } + std::cerr << smallestNorm << std::endl; - std::cerr << "Im H is now sufficently close, starting to relax rope." << std::endl; + J = exp(Complex(0, φ)) * J; - std::cerr << "Threshold proportions of saddles are " << getProportionOfThreshold(κ, J, z1) << " and " << getProportionOfThreshold(κ, J, z2) << std::endl; + /* + if (stokesLineTest<Real>(J, zSaddle, zSaddleNext, 10, 4)) { + std::cerr << "Found a Stokes line" << std::endl; +// stokesLineTestNew<Real>(J, zSaddle, zSaddleNext, 10, 3); + } else { + std::cerr << "Didn't find a Stokes line" << std::endl; + } + */ - Rope stokes(10, z1, z2, J); - for (unsigned i = 0; i < 9; i++) { - stokes.relaxDiscreteGradient(J, 1e6, 0.1, 0); + Cord test(J, zSaddle, zSaddleNext, 5); + test.relax(J, 20, 1, 1e5); - std::cout << stokes.n() << " " << stokes.cost(J) << " " << stokes.length() << std::endl; + std::cout << test.z0.transpose() << std::endl; + std::cout << test.z1.transpose() << std::endl; - stokes = stokes.interpolate(); + for (Vector<Complex>& g : test.gs) { + std::cout << g.transpose() << std::endl; } return 0; |