diff options
Diffstat (limited to 'langevin.cpp')
| -rw-r--r-- | langevin.cpp | 156 |
1 files changed, 26 insertions, 130 deletions
diff --git a/langevin.cpp b/langevin.cpp index a4c9cf3..870879b 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -1,127 +1,21 @@ -#include <exception> #include <getopt.h> -#include <eigen3/Eigen/LU> -#include <utility> - #include "complex_normal.hpp" #include "p-spin.hpp" -#include "stereographic.hpp" +#include "dynamics.hpp" #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" #include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h" -using Rng = randutils::random_generator<pcg32>; - -Vector normalize(const Vector& z) { - return z * sqrt((double)z.size() / (Scalar)(z.transpose() * z)); -} - -template <class Distribution, class Generator> -Vector randomVector(unsigned N, Distribution d, Generator& r) { - Vector z(N); - - for (unsigned i = 0; i < N; i++) { - z(i) = d(r); - } - - return z; -} - -class gradientDescentStallException: public std::exception { - virtual const char* what() const throw() { - return "Gradient descent stalled."; - } -} gradientDescentStall; - -std::tuple<double, Vector> gradientDescent(const Tensor& J, const Vector& z0, double ε, double γ0 = 1, double δγ = 2) { - Vector z = z0; - double γ = γ0; - - auto [W, dW] = WdW(J, z); - - while (W > ε) { - Vector zNew = normalize(z - γ * dW.conjugate()); - - auto [WNew, dWNew] = WdW(J, zNew); - - if (WNew < W) { // If the step lowered the objective, accept it! - z = zNew; - W = WNew; - dW = dWNew; - γ = γ0; - } else { // Otherwise, shrink the step and try again. - γ /= δγ; - } - - if (γ < 1e-50) { - throw gradientDescentStall; - } - } +#define PSPIN_P 3 +const int p = PSPIN_P; // polynomial degree of Hamiltonian +using Complex = std::complex<double>; +using ComplexVector = Vector<Complex>; +using ComplexMatrix = Matrix<Complex>; +using ComplexTensor = Tensor<Complex, p>; - return {W, z}; -} - -Vector findSaddle(const Tensor& J, const Vector& z0, double ε, double δW = 2, double γ0 = 1, double δγ = 2) { - Vector z = z0; - Vector ζ = euclideanToStereographic(z); - - double W; - std::tie(W, std::ignore) = WdW(J, z); - - Vector dH; - Matrix ddH; - std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z); - - while (W > ε) { - // ddH is complex symmetric, which is (almost always) invertible, so a - // partial pivot LU decomposition can be used. - Vector dζ = ddH.partialPivLu().solve(dH); - Vector ζNew = ζ - dζ; - Vector zNew = stereographicToEuclidean(ζNew); - - double WNew; - std::tie(WNew, std::ignore) = WdW(J, zNew); - - if (WNew < W) { // If Newton's step lowered the objective, accept it! - ζ = ζNew; - z = zNew; - W = WNew; - } else { // Otherwise, do gradient descent until W is a factor δW smaller. - std::tie(W, z) = gradientDescent(J, z, W / δW, γ0, δγ); - ζ = euclideanToStereographic(z); - } - - std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z); - } - - return z; -} - -std::tuple<double, Vector> langevin(const Tensor& J, const Vector& z0, double T, double γ, unsigned N, Rng& r) { - Vector z = z0; - - double W; - std::tie(W, std::ignore) = WdW(J, z); - - complex_normal_distribution<> d(0, γ, 0); - - for (unsigned i = 0; i < N; i++) { - Vector dz = randomVector(z.size(), d, r.engine()); - Vector zNew = normalize(z + dz); - - double WNew; - std::tie(WNew, std::ignore) = WdW(J, zNew); - - if (exp((W - WNew) / T) > r.uniform(0.0, 1.0)) { - z = zNew; - W = WNew; - } - } - - return {W, z}; -} +using Rng = randutils::random_generator<pcg32>; int main(int argc, char* argv[]) { // model parameters @@ -178,22 +72,24 @@ int main(int argc, char* argv[]) { } } - Scalar κ(Rκ, Iκ); + Complex κ(Rκ, Iκ); double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1))); Rng r; - Tensor J = generateCouplings<Scalar, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine()); - Vector z0 = normalize(randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine())); + 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())); - Vector zSaddle = findSaddle(J, z0, ε); - Vector zSaddlePrev = Vector::Zero(N); - Vector z = zSaddle; + ComplexVector zSaddle = findSaddle(J, z0, ε); + ComplexVector zSaddlePrev = ComplexVector::Zero(N); + ComplexVector z = zSaddle; while (δ < (zSaddle - zSaddlePrev).norm()) { // Until we find two saddles sufficiently close... - std::tie(std::ignore, z) = langevin(J, z, T, γ, M, r); + std::tie(std::ignore, z) = langevin(J, z, T, γ, M, d, r.engine()); try { - Vector zSaddleNext = findSaddle(J, z, ε); + ComplexVector zSaddleNext = findSaddle(J, z, ε); if (Δ < (zSaddleNext - zSaddle).norm()) { // Ensure we are finding distinct saddles. zSaddlePrev = zSaddle; zSaddle = zSaddleNext; @@ -209,20 +105,20 @@ int main(int argc, char* argv[]) { complex_normal_distribution<> dJ(0, εJ * σ, 0); - std::function<void(Tensor&, std::array<unsigned, p>)> perturbJ = - [&dJ, &r] (Tensor& JJ, std::array<unsigned, p> ind) { - Scalar Ji = getJ<Scalar, p>(JJ, ind); - setJ<Scalar, p>(JJ, ind, Ji + dJ(r.engine())); + std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ = + [&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())); }; for (unsigned i = 0; i < n; i++) { - Tensor Jp = J; + ComplexTensor Jp = J; - iterateOver<Scalar, p>(Jp, perturbJ); + iterateOver<Complex, p>(Jp, perturbJ); try { - Vector zSaddleNew = findSaddle(Jp, zSaddle, ε); - Vector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε); + ComplexVector zSaddleNew = findSaddle(Jp, zSaddle, ε); + ComplexVector zSaddlePrevNew = findSaddle(Jp, zSaddlePrev, ε); std::cout << zSaddleNew.transpose() << " " << zSaddlePrevNew.transpose() << std::endl; } catch (std::exception& e) { |