From 7bf6e952b53699977f5091a78f0f9f48f7b359c5 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 9 Nov 2021 10:13:59 +0100 Subject: Generalized code to easily allow mixed p-spins. --- dynamics.hpp | 50 ++++++++++++++++------------- langevin.cpp | 95 ++++++++++--------------------------------------------- p-spin.hpp | 101 +++++++++++++++++++++++++++++++++++++---------------------- stokes.hpp | 70 +++++++++++------------------------------ tensor.hpp | 18 ++--------- 5 files changed, 129 insertions(+), 205 deletions(-) diff --git a/dynamics.hpp b/dynamics.hpp index 7d04b1a..36b2dfa 100644 --- a/dynamics.hpp +++ b/dynamics.hpp @@ -6,22 +6,28 @@ #include "p-spin.hpp" -template -Vector findMinimum(const Tensor& J, const Vector& z0, Real ε) { +template +Vector findMinimum(const pSpinModel& M, const Vector& z0, Real ε) { Vector z = z0; Real λ = 10; - auto [H, dH, ddH] = hamGradHess(J, z); + Real H; + Vector dH; + Matrix ddH; + std::tie(H, dH, ddH, std::ignore) = M.hamGradHess(z); Vector g = dH - z.dot(dH) * z / (Real)z.size(); - Matrix M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose(); + Matrix m = ddH - (dH * z.transpose() + z.dot(dH) * Matrix::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose(); while (g.norm() / z.size() > ε && λ < 1e8) { - Vector dz = (M + λ * (Matrix)abs(M.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g); + Vector dz = (m + λ * (Matrix)abs(m.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g); dz -= z.dot(dz) * z / (Real)z.size(); Vector zNew = normalize(z - dz); - auto [HNew, dHNew, ddHNew] = hamGradHess(J, zNew); + Real HNew; + Vector dHNew; + Matrix ddHNew; + std::tie(HNew, dHNew, ddHNew, std::ignore) = M.hamGradHess(zNew); if (HNew * 1.0001 <= H) { z = zNew; @@ -30,7 +36,7 @@ Vector findMinimum(const Tensor& J, const Vector& z0, Real ddH = ddHNew; g = dH - z.dot(dH) * z / (Real)z.size(); - M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose(); + m = ddH - (dH * z.transpose() + z.dot(dH) * Matrix::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose(); λ /= 2; } else { @@ -41,28 +47,28 @@ Vector findMinimum(const Tensor& J, const Vector& z0, Real return z; } -template -Vector findSaddle(const Tensor& J, const Vector& z0, Real ε) { +template +Vector findSaddle(const pSpinModel& M, const Vector& z0, Real ε) { Vector z = z0; Vector dH; Matrix ddH; - std::tie(std::ignore, dH, ddH) = hamGradHess(J, z); + std::tie(std::ignore, dH, ddH, std::ignore) = M.hamGradHess(z); Scalar zz = z.transpose() * z; Vector g = dH - (z.transpose() * dH) * z / (Real)z.size() + z * (zz - (Real)z.size()); - Matrix M = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + Matrix::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose(); + Matrix m = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + Matrix::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose(); while (g.norm() / z.size() > ε) { - Vector dz = M.partialPivLu().solve(g); + Vector dz = m.partialPivLu().solve(g); dz -= z.conjugate().dot(dz) / z.squaredNorm() * z.conjugate(); z = normalize(z - dz); - std::tie(std::ignore, dH, ddH) = hamGradHess(J, z); + std::tie(std::ignore, dH, ddH, std::ignore) = M.hamGradHess(z); zz = z.transpose() * z; g = dH - (z.transpose() * dH) * z / (Real)z.size() + z * (zz - (Real)z.size()); - M = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + Matrix::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose(); + m = ddH - (dH * z.transpose() + (z.transpose() * dH) * Matrix::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + Matrix::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose(); } return z; @@ -79,16 +85,16 @@ Vector randomVector(unsigned N, Distribution d, Generator& r) { return z; } -template -Vector randomMinimum(const Tensor& J, Distribution d, Generator& r, Real ε) { +template +Vector randomMinimum(const pSpinModel& M, Distribution d, Generator& r, Real ε) { Vector zSaddle; bool foundSaddle = false; while (!foundSaddle) { - Vector z0 = normalize(randomVector(J.dimension(0), d, r.engine())); + Vector z0 = normalize(randomVector(M.dimension(), d, r.engine())); try { - zSaddle = findMinimum(J, z0, ε); + zSaddle = findMinimum(M, z0, ε); foundSaddle = true; } catch (std::exception& e) {} } @@ -96,16 +102,16 @@ Vector randomMinimum(const Tensor& J, Distribution d, Generator& return zSaddle; } -template -Vector randomSaddle(const Tensor& J, Distribution d, Generator& r, Real ε) { +template +Vector randomSaddle(const pSpinModel& M, Distribution d, Generator& r, Real ε) { Vector zSaddle; bool foundSaddle = false; while (!foundSaddle) { - Vector z0 = normalize(randomVector(J.dimension(0), d, r.engine())); + Vector z0 = normalize(randomVector(M.dimension(), d, r.engine())); try { - zSaddle = findSaddle(J, z0, ε); + zSaddle = findSaddle(M, z0, ε); foundSaddle = true; } catch (std::exception& e) {} } diff --git a/langevin.cpp b/langevin.cpp index 86aeb9e..dc3d07d 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -27,63 +27,6 @@ using ComplexTensor = Tensor; using Rng = randutils::random_generator; -std::list> saddlesCloserThan(const std::unordered_map& saddles, Real δ) { - std::list> 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; -} - -template -std::tuple matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, Real ε, Real Δ, Generator& r) { - Real σ = sqrt(factorial(p) / (Real(2) * pow(z10.size(), p - 1))); - complex_normal_distribution dJ(0, σ, 0); - - ComplexTensor J = J0; - ComplexVector z1, z2; - Real ΔH = abs(imag(getHamiltonian(J, z10) - getHamiltonian(J, z20))) / z10.size(); - Real γ = ΔH; - - std::function)> perturbJ = - [&γ, &dJ, &r] (ComplexTensor& JJ, std::array ind) { - Complex Ji = getJ(JJ, ind); - setJ(JJ, ind, Ji + γ * dJ(r.engine())); - }; - - while (ΔH > ε) { - ComplexTensor Jp = J; - iterateOver(Jp, perturbJ); - - try { - z1 = findSaddle(Jp, z10, Δ); - z2 = findSaddle(Jp, z20, Δ); - - Real Δz = (z1 - z2).norm(); - - if (Δz > 1e-2) { - Real ΔHNew = abs(imag(getHamiltonian(Jp, z1) - getHamiltonian(Jp, z2))) / z1.size(); - - if (ΔHNew < ΔH) { - J = Jp; - ΔH = ΔHNew; - γ = ΔH; - - std::cerr << "Match imaginary energies: Found couplings with ΔH = " << ΔH << std::endl; - } - } - } catch (std::exception& e) {} - } - - return {J, z1, z2}; -} - int main(int argc, char* argv[]) { // model parameters unsigned N = 10; // number of spins @@ -143,24 +86,29 @@ int main(int argc, char* argv[]) { Rng r; - Tensor ReJ = generateCouplings(N, std::normal_distribution(0, σ), r.engine()); + pSpinModel pSpin; + + std::get<0>(pSpin.Js) = generateCouplings(N, std::normal_distribution(0, σ), r.engine()); std::normal_distribution Red(0, 1); - RealVector zMin = randomMinimum(ReJ, Red, r, ε); - auto [Hr, dHr, ddHr] = hamGradHess(ReJ, zMin); + RealVector zMin = randomMinimum(pSpin, Red, r, ε); + Real Hr; + RealVector dHr; + RealMatrix ddHr; + std::tie(Hr, dHr, ddHr, std::ignore) = pSpin.hamGradHess(zMin); Eigen::EigenSolver> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + 2.0 * zMin * zMin.transpose()); std::cout << eigenS.eigenvalues().transpose() << std::endl; for (unsigned i = 0; i < N; i++) { RealVector zNew = normalize(zMin + 0.01 * eigenS.eigenvectors().col(i).real()); - std::cout << getHamiltonian(ReJ, zNew) - Hr << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << std::endl; + std::cout << pSpin.getHamiltonian(zNew) - Hr << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << std::endl; } std::cout << std::endl; getchar(); complex_normal_distribution d(0, 1, 0); - ComplexTensor J = ReJ.cast(); + pSpinModel complexPSpin = pSpin.cast();; ComplexVector zSaddle = zMin.cast(); ComplexVector zSaddleNext; @@ -168,7 +116,7 @@ int main(int argc, char* argv[]) { while (!foundSaddle) { ComplexVector z0 = normalize(zSaddle + δ * randomVector(N, d, r.engine())); try { - zSaddleNext = findSaddle(J, z0, ε); + zSaddleNext = findSaddle(complexPSpin, z0, ε); Real saddleDistance = (zSaddleNext - zSaddle).norm(); if (saddleDistance / N > 1e-2) { std::cout << saddleDistance << std::endl; @@ -177,27 +125,16 @@ int main(int argc, char* argv[]) { } catch (std::exception& e) {} } - auto [H1, dH1, ddH1] = hamGradHess(J, zSaddle); - auto [H2, dH2, ddH2] = hamGradHess(J, zSaddleNext); - - Eigen::ComplexEigenSolver ces; - ces.compute(ddH1); + Complex H1 = complexPSpin.getHamiltonian(zSaddle); + Complex H2 = complexPSpin.getHamiltonian(zSaddleNext); 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::infinity(); - for (Complex z : ces.eigenvalues()) { - if (norm(z) < smallestNorm) { - smallestNorm = norm(z); - } - } - std::cerr << smallestNorm << std::endl; - getchar(); - J = exp(Complex(0, φ)) * J; + std::get<0>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<0>(complexPSpin.Js); - Cord test(J, zSaddle, zSaddleNext, 3); - test.relaxNewton(J, 10, 1, 1e4); + Cord test(complexPSpin, zSaddle, zSaddleNext, 3); + test.relaxNewton(10, 1, 1e4); std::cout << test.z0.transpose() << std::endl; std::cout << test.z1.transpose() << std::endl; diff --git a/p-spin.hpp b/p-spin.hpp index 6111b75..a5f6c9a 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -12,66 +12,93 @@ Vector normalize(const Eigen::MatrixBase& z) return z * sqrt((Real)z.size() / (typename Derived::Scalar)(z.transpose() * z)); } -template -std::tuple, Matrix> hamGradHess(const Tensor& J, const Tensor& ... Js, const Vector& z) { - Matrix Jz = contractDown(J, z); // Contracts J into p - 2 copies of z. +template +std::tuple, Tensor> hamGradTensorHelper(const Vector& z, const Tensor& J) { + Tensor J3(z.size(), z.size(), z.size());; + J3.setZero(); + Matrix Jz = Eigen::Map>(J.data(), z.size(), z.size()); + + return {Jz, J3}; +} + +template +std::tuple, Tensor> hamGradTensorHelper(const Vector& z, const Tensor& J) { + Tensor J3 = contractDown(J, z); + Tensor zT = Eigen::TensorMap>(z.data(), {z.size()}); + Tensor J3zT = J3.contract(zT, ip00); + Matrix Jz = Eigen::Map>(J3zT.data(), z.size(), z.size()); + + return {Jz, J3}; +} + +template +std::tuple, Matrix, Tensor> hamGradHessHelper(const Vector& z, const Tensor& J, const Tensor& ...Js) { + auto [Jz, J3] = hamGradTensorHelper(z, J); + Vector Jzz = Jz * z; Scalar Jzzz = Jzz.transpose() * z; Real pBang = factorial(p); + Tensor dddH = ((p - 2) * (p - 1) * p / pBang) * J3; Matrix ddH = ((p - 1) * p / pBang) * Jz; Vector dH = (p / pBang) * Jzz; Scalar H = Jzzz / pBang; if constexpr (sizeof...(Js) > 0) { - auto [Hs, dHs, ddHs] = hamGradHess(Js..., z); + auto [Hs, dHs, ddHs, dddHs] = hamGradHess(z, Js...); H += Hs; dH += dHs; ddH += ddHs; + dddH += dddHs; } - return {H, dH, ddH}; + return {H, dH, ddH, dddH}; } -template -Scalar getHamiltonian(const Tensor& J, const Vector& z) { - Scalar H; - std::tie(H, std::ignore, std::ignore) = hamGradHess(J, z); - return H; -} +template +class pSpinModel { +public: + std::tuple...> Js; -template -Vector getGradient(const Tensor& J, const Vector& z) { - Vector dH; - std::tie(std::ignore, dH, std::ignore) = hamGradHess(J, z); - return dH; -} + unsigned dimension() const { + return std::get<0>(Js).dimension(0); + } -template -Matrix getHessian(const Tensor& J, const Vector& z) { - Matrix ddH; - std::tie(std::ignore, std::ignore, ddH) = hamGradHess(J, z); - return ddH; -} + template + pSpinModel cast() const { + pSpinModel M; + M.Js = std::apply( + [] (const Tensor& ...Ks) -> std::tuple...> { + return std::make_tuple(Ks.template cast()...); + }, Js + ); -template -Real getThresholdEnergyDensity(unsigned p, Scalar κ, Scalar ε, Real a) { - Real apm2 = pow(a, p - 2); - Scalar δ = κ / apm2; - Real θ = arg(κ) + 2 * arg(ε); + return M; + } - return (p - 1) * apm2 / (2 * p) * pow(1 - norm(δ), 2) / (1 + norm(δ) - 2 * abs(δ) * cos(θ)); -} + std::tuple, Matrix, Tensor> hamGradHess(const Vector& z) const { + return std::apply([&z](const Tensor& ...Ks) -> std::tuple, Matrix, Tensor> { return hamGradHessHelper(z, Ks...); }, Js); + } -template -Real getProportionOfThreshold(Scalar κ, const Tensor& J, const Vector& z) { - Real N = z.size(); - Scalar ε = getHamiltonian(J, z) / N; - Real a = z.squaredNorm() / N; + Scalar getHamiltonian(const Vector& z) const { + Scalar H; + std::tie(H, std::ignore, std::ignore, std::ignore) = hamGradHess(z); + return H; + } - return norm(ε) / getThresholdEnergyDensity(p, κ, ε, a); -} + Vector getGradient(const Vector& z) const { + Vector dH; + std::tie(std::ignore, dH, std::ignore, std::ignore) = hamGradHess(z); + return dH; + } + + Matrix getHessian(const Vector& z) const { + Matrix ddH; + std::tie(std::ignore, std::ignore, ddH, std::ignore) = hamGradHess(z); + return ddH; + } +}; template Vector zDot(const Vector& z, const Vector& dH) { diff --git a/stokes.hpp b/stokes.hpp index b3ae90b..c2b0ec5 100644 --- a/stokes.hpp +++ b/stokes.hpp @@ -4,17 +4,18 @@ #include "complex_normal.hpp" #include "dynamics.hpp" -template +template class Cord { public: + const pSpinModel& M; std::vector> gs; Vector z0; Vector z1; + double φ; - template - Cord(const Tensor& J, const Vector& z2, const Vector& z3, unsigned ng) : gs(ng, Vector::Zero(z2.size())) { - Scalar H2 = getHamiltonian(J, z2); - Scalar H3 = getHamiltonian(J, z3); + Cord(const pSpinModel& M, const Vector& z2, const Vector& z3, unsigned ng) : M(M), gs(ng, Vector::Zero(z2.size())) { + Scalar H2 = M.getHamiltonian(z2); + Scalar H3 = M.getHamiltonian(z3); if (real(H2) > real(H3)) { z0 = z2; @@ -53,64 +54,31 @@ public: return z; } - template - Real cost(const Tensor& J, Real t) const { + Real cost(Real t) const { Vector z = f(t); Scalar H; Vector dH; - std::tie(H, dH, std::ignore) = hamGradHess(J, z); + std::tie(H, dH, std::ignore, std::ignore) = M.hamGradHess(z); Vector ż = zDot(z, dH); Vector dz = df(t); return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm(); } - template - Real totalCost(const Tensor& J, unsigned nt) const { + Real totalCost(unsigned nt) const { Real tc = 0; for (unsigned i = 0; i < nt; i++) { Real t = (i + 1.0) / (nt + 1.0); - tc += cost(J, t); + tc += cost(t); } return tc; } - template - std::vector> dgs(const Tensor& J, Real t) const { + std::pair, Matrix> gsGradHess(Real t) const { Vector z = f(t); - auto [H, dH, ddH] = hamGradHess(J, z); - Vector ż = zDot(z, dH); - Vector dz = df(t); - Matrix dż = dzDot(z, dH); - Matrix dżc = dzDotConjugate(z, dH, ddH); - - std::vector> x; - x.reserve(gs.size()); - - for (unsigned i = 0; i < gs.size(); i++) { - Real fdg = gCoeff(i, t); - Real dfdg = dgCoeff(i, t); - Vector dC = - 0.5 / ż.norm() / dz.norm() * ( - dfdg * ż.conjugate() + fdg * dżc * dz + fdg * dż * dz.conjugate() - - real(dz.dot(ż)) * ( - dfdg * dz.conjugate() / dz.squaredNorm() + - fdg * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm() - ) - ); - - x.push_back(dC.conjugate()); - } - - return x; - } - - template - std::pair, Matrix> gsGradHess(const Tensor& J, Real t) const { - Vector z = f(t); - auto [H, dH, ddH] = hamGradHess(J, z); - Tensor dddH = ((p - 2) * (p - 1) * p / (Real)factorial(p)) * J; // BUG IF P != 3 + auto [H, dH, ddH, dddH] = M.hamGradHess(z); Vector ż = zDot(z, dH); Tensor żT = Eigen::TensorMap>(ż.data(), {z.size()}); Vector dz = df(t); @@ -241,14 +209,13 @@ public: return {x, M}; } - template - std::pair relaxStepNewton(const Tensor& J, unsigned nt, Real δ₀) { + std::pair relaxStepNewton(unsigned nt, Real δ₀) { Vector dgsTot = Vector::Zero(2 * z0.size() * gs.size()); Matrix HessTot = Matrix::Zero(2 * z0.size() * gs.size(), 2 * z0.size() * gs.size()); for (unsigned i = 0; i < nt; i++) { Real t = (i + 1.0) / (nt + 1.0); - auto [dgsi, Hessi] = gsGradHess(J, t); + auto [dgsi, Hessi] = gsGradHess(t); dgsTot += dgsi / nt; HessTot += Hessi / nt; @@ -259,7 +226,7 @@ public: Cord cNew(*this); Real δ = δ₀; - Real oldCost = totalCost(J, nt); + Real oldCost = totalCost(nt); Real newCost = std::numeric_limits::infinity(); Matrix I = abs(HessTot.diagonal().array()).matrix().asDiagonal(); @@ -273,7 +240,7 @@ public: } } - newCost = cNew.totalCost(J, nt); + newCost = cNew.totalCost(nt); δ *= 2; } @@ -283,13 +250,12 @@ public: return {δ / 2, stepSize}; } - template - void relaxNewton(const Tensor& J, unsigned nt, Real δ₀, unsigned maxSteps) { + void relaxNewton(unsigned nt, Real δ₀, unsigned maxSteps) { Real δ = δ₀; Real stepSize = std::numeric_limits::infinity(); unsigned steps = 0; while (stepSize > 1e-7 && steps < maxSteps) { - std::tie(δ, stepSize) = relaxStepNewton(J, nt, δ); + std::tie(δ, stepSize) = relaxStepNewton(nt, δ); δ /= 3; steps++; } diff --git a/tensor.hpp b/tensor.hpp index aa33069..fc99042 100644 --- a/tensor.hpp +++ b/tensor.hpp @@ -111,27 +111,15 @@ Tensor plantState(const Tensor& J, const Vector& z } template -Matrix contractDown(const Tensor& J, const Vector& z) { - return Eigen::Map>(J.data(), z.size(), z.size()); +Tensor contractDown(const Tensor& J, const Vector& z) { + return J; } const std::array, 1> ip00 = {Eigen::IndexPair(0, 0)}; template -Matrix contractDown(const Tensor& J, const Vector& z) { +Tensor 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); } - -template -Tensor contractDownTo(const Tensor& J, const Vector& z) { - return J; -} - -template -Tensor contractDownTo(const Tensor& J, const Vector& z) { - Tensor zT = Eigen::TensorMap>(z.data(), {z.size()}); - Tensor Jz = J.contract(zT, ip00); - return contractDownTo(Jz, z); -} -- cgit v1.2.3-70-g09d2