From bcc5c327016a3d188a9bdcb0069cc63b0bb8163f Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 9 Nov 2021 15:41:56 +0100 Subject: Moved some repeated code into the pSpinModel class. --- langevin.cpp | 23 ++++--------- p-spin.hpp | 109 +++++++++++++++++++++++++++++++++++------------------------ tensor.hpp | 24 ++++--------- 3 files changed, 77 insertions(+), 79 deletions(-) diff --git a/langevin.cpp b/langevin.cpp index 8ecd69c..0e5a616 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -81,17 +81,10 @@ int main(int argc, char* argv[]) { } } - Complex κ(Rκ, Iκ); - Real σp = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1))); - Real σ2 = sqrt(factorial(2) / ((Real)2 * pow(N, 2 - 1))); - Real σ4 = sqrt(factorial(4) / ((Real)2 * pow(N, 4 - 1))); - Rng r; - pSpinModel pSpin; - - std::get<0>(pSpin.Js) = generateCouplings(N, std::normal_distribution(0, σ2), r.engine()); - std::get<1>(pSpin.Js) = generateCouplings(N, std::normal_distribution(0, σ4 / 10), r.engine()); +// pSpinModel pSpin({J2, J4}); + pSpinModelpSpin(N, r.engine(), 1, 0.1); std::normal_distribution Red(0, 1); @@ -100,18 +93,17 @@ int main(int argc, char* argv[]) { 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; + Eigen::EigenSolver> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + zMin * zMin.transpose() / (Real)N); for (unsigned i = 0; i < N; i++) { - RealVector zNew = normalize(zMin + 0.01 * eigenS.eigenvectors().col(i).real()); - std::cout << pSpin.getHamiltonian(zNew) - Hr << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << std::endl; + RealVector zNew = normalize(zMin + 1e-3 * eigenS.eigenvectors().col(i).real()); + std::cout << eigenS.eigenvalues()(i) << " " << 2.0 * (pSpin.getHamiltonian(zNew) - Hr) / 1e-6 << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << " " << eigenS.eigenvectors().col(0).dot(eigenS.eigenvectors().col(i)) << std::endl; } std::cout << std::endl; getchar(); complex_normal_distribution d(0, 1, 0); - pSpinModel complexPSpin = pSpin.cast();; + pSpinModel complexPSpin = pSpin.cast();; ComplexVector zSaddle = zMin.cast(); ComplexVector zSaddleNext; @@ -134,8 +126,7 @@ int main(int argc, char* argv[]) { 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; - std::get<0>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<0>(complexPSpin.Js); - std::get<1>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<1>(complexPSpin.Js); + complexPSpin *= exp(Complex(0, φ)); Cord test(complexPSpin, zSaddle, zSaddleNext, 3); test.relaxNewton(10, 1, 1e4); diff --git a/p-spin.hpp b/p-spin.hpp index 7a49bca..544abf1 100644 --- a/p-spin.hpp +++ b/p-spin.hpp @@ -12,73 +12,92 @@ Vector normalize(const Eigen::MatrixBase& z) return z * sqrt((Real)z.size() / (typename Derived::Scalar)(z.transpose() * 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 +class pSpinModel { +private: + std::tuple, Tensor> hamGradTensorHelper(const Vector& z, const Tensor& J) const { + Tensor J3(z.size(), z.size(), z.size());; + J3.setZero(); + Matrix Jz = Eigen::Map>(J.data(), z.size(), z.size()); -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}; + } - return {Jz, J3}; -} + template + std::tuple, Tensor> hamGradTensorHelper(const Vector& z, const Tensor& J) const { + 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()); -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, dddHs] = hamGradHessHelper(z, Js...); - H += Hs; - dH += dHs; - ddH += ddHs; - dddH += dddHs; + return {Jz, J3}; } - return {H, dH, ddH, dddH}; -} + template + std::tuple, Matrix, Tensor> hamGradHessHelper(const Vector& z, const Tensor& J, const Tensor& ...Js) const { + 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, dddHs] = hamGradHessHelper(z, Js...); + H += Hs; + dH += dHs; + ddH += ddHs; + dddH += dddHs; + } + + return {H, dH, ddH, dddH}; + } -template -class pSpinModel { public: std::tuple...> Js; + pSpinModel() {} + + pSpinModel(const std::tuple...>& Js) : Js(Js) {} + + template + pSpinModel(unsigned N, Generator& r, T... μs) { + Js = std::make_tuple(μs * generateRealPSpinCouplings(N, r)...); + } + unsigned dimension() const { return std::get<0>(Js).dimension(0); } template pSpinModel cast() const { - pSpinModel M; - M.Js = std::apply( + return std::apply( [] (const Tensor& ...Ks) -> std::tuple...> { return std::make_tuple(Ks.template cast()...); }, Js - ); + ); + } + + template + pSpinModel& operator*=(T x) { + std::tuple...> newJs = std::apply( + [x] (const Tensor& ...Ks) -> std::tuple...> { + return std::make_tuple((x * Ks)...); + }, Js + ); + + Js = newJs; - return M; + return *this; } 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); + return std::apply([&z, this](const Tensor& ...Ks) -> std::tuple, Matrix, Tensor> { return hamGradHessHelper(z, Ks...); }, Js); } Scalar getHamiltonian(const Vector& z) const { diff --git a/tensor.hpp b/tensor.hpp index fc99042..1212bda 100644 --- a/tensor.hpp +++ b/tensor.hpp @@ -5,6 +5,8 @@ #include +#include "factorial.hpp" + template using Vector = Eigen::Matrix; @@ -89,25 +91,11 @@ Tensor generateCouplings(unsigned N, Distribution d, Generator& r) { return J; } -template -Tensor plantState(const Tensor& J, const Vector& z, double β) { - Tensor JPlanted = J; - - std::function&, std::array)> plant = - [&z, β] (Tensor& JJ, std::array ind) { - Scalar Ji = getJ(JJ, ind); - Scalar zzz = 1; - - for (unsigned i : ind) { - zzz *= z(i); - } - - setJ(JJ, ind, Ji - β * zzz / pow(zzz.size(), p)); - }; - - iterateOver(JPlanted, plant); +template +Tensor generateRealPSpinCouplings(unsigned N, Generator& r) { + Real σp = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1))); - return JPlanted; + return generateCouplings(N, std::normal_distribution(0, σp), r); } template -- cgit v1.2.3-70-g09d2