From f51fad699957867c6eff1c4d59390513b594cb8c Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 11 Nov 2021 17:10:44 +0100 Subject: Cleaned up code and algorithmic improvements for data collection. --- collectStokesData.hpp | 76 +++++++++++++++++++------ dynamics.hpp | 2 +- langevin.cpp | 142 ----------------------------------------------- pureStokesFromMinima.cpp | 54 ++++++++++++++++++ stokes.hpp | 56 +++++++++---------- 5 files changed, 141 insertions(+), 189 deletions(-) delete mode 100644 langevin.cpp create mode 100644 pureStokesFromMinima.cpp diff --git a/collectStokesData.hpp b/collectStokesData.hpp index 50aad24..433f5ad 100644 --- a/collectStokesData.hpp +++ b/collectStokesData.hpp @@ -1,21 +1,34 @@ #include "stokes.hpp" -#include +#include "complex_normal.hpp" +#include using Complex = std::complex; template -void collectStokesData(unsigned N, Generator& r, double ε, T... μs) { +void collectStokesData(std::ofstream& file, unsigned N, Generator& r, double ε, double δz, bool minimum, T... μs) { + unsigned nGs = 8; + unsigned nTs = 20; + Real newSaddleThres = 1e-3; + pSpinModel ReM(N, r, μs...); std::normal_distribution Red(0, 1); - Vector xMin = randomMinimum(ReM, Red, r, ε); + Vector xMin(N); + + if (minimum) { + xMin = randomMinimum(ReM, Red, r, ε); + } else { + xMin = randomSaddle(ReM, Red, r, ε); + } Real Hx; Vector dHx; Matrix ddHx; std::tie(Hx, dHx, ddHx, std::ignore) = ReM.hamGradHess(xMin); - Eigen::EigenSolver> eigenS(ddHx - xMin.dot(dHx) * Matrix::Identity(xMin.size(), xMin.size()) / Real(N)); + Matrix xHess = ddHx - xMin.dot(dHx) * Matrix::Identity(xMin.size(), xMin.size()) / (Real)N; + xHess -= (xHess * xMin) * xMin.transpose() / (Real)N; + Eigen::EigenSolver> eigenS(xHess); complex_normal_distribution d(0, 1, 0); @@ -24,29 +37,56 @@ void collectStokesData(unsigned N, Generator& r, double ε, T... μs) { Vector zSaddle = zMin; - while ((zSaddle - zMin).norm() < 1e-3 * N || abs(imag(M.getHamiltonian(zSaddle))) < 1e-10) { - Vector z0 = normalize(zSaddle + 0.5 * randomVector(N, d, r)); + while ((zSaddle - zMin).norm() < newSaddleThres * N || abs(imag(M.getHamiltonian(zSaddle))) < 1e-10) { + Vector z0 = normalize(zSaddle + δz * randomVector(N, d, r)); zSaddle = findSaddle(M, z0, ε); } - Complex H1 = M.getHamiltonian(zMin); - Complex H2 = M.getHamiltonian(zSaddle); + Complex Hz; + Vector dHz; + Matrix ddHz; + std::tie(Hz, dHz, ddHz, std::ignore) = M.hamGradHess(zSaddle); + Matrix zHess = ddHz - (zSaddle.transpose() * dHz) * Matrix::Identity(N, N) / (Real)N; + zHess -= (zHess * zSaddle) * zSaddle.adjoint() / zSaddle.squaredNorm(); + Eigen::ComplexEigenSolver> eigenSz(zHess); - Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real()); + Real φ = atan2(Hz.imag(), Hx - Hz.real()); M *= exp(Complex(0, φ)); - Cord c(M, zMin, zSaddle, 3); - c.relaxNewton(10, 1, 1e4); + file.precision(15); - ((std::cout << ps), ...) << std::endl;; - ((std::cout << μs), ...) << std::endl;; + file << N << std::endl; + ((file << ps), ...) << std::endl;; + ((file << μs), ...) << std::endl;; - std::apply([](const Tensor&... Js) -> void { - ((std::cout << Js << std::endl), ...); + std::apply([&file](const Tensor&... Js) -> void { + ((file << Js << std::endl), ...); } , ReM.Js); - std::cout << xMin.transpose() << std::endl; - std::cout << Hx << std::endl; - std::cout << eigenS.eigenvalues().real().transpose() << std::endl; + file << xMin.transpose() << std::endl; + file << Hx << std::endl; + file << eigenS.eigenvalues().real().transpose() << std::endl; + file << zSaddle.transpose() << std::endl; + file << Hz << std::endl; + file << eigenSz.eigenvalues().transpose() << std::endl; + file << φ << std::endl; + + Cord c(M, zMin, zSaddle, nGs); + c.relaxNewton(nTs, 1, 1e-10, 1e3); + + Complex constraintError = 0; + Real imEnergyError = 0; + + for (unsigned i = 0; i < 100 * nTs; i++) { + Vector zi = c.f((i + 1.0) / (100.0 * nTs + 1.0)); + imEnergyError += pow(std::imag(M.getHamiltonian(zi) - M.getHamiltonian(zMin)), 2); + constraintError += pow(((Complex)(zi.transpose() * zi) - (Complex)N), 2); + } + + file << nGs << " " << nTs << std::endl;; + file << c.totalCost(100 * nTs) / (100 * nTs) << " " << sqrt(imEnergyError) / (100 * nTs) << " " << sqrt(constraintError) / (100.0 * nTs) << std::endl; + for (const Vector& gi : c.gs) { + file << gi.transpose() << std::endl; + } } diff --git a/dynamics.hpp b/dynamics.hpp index 4f33c1a..e578cdb 100644 --- a/dynamics.hpp +++ b/dynamics.hpp @@ -108,7 +108,7 @@ Vector randomSaddle(const pSpinModel& M, Distribution d, Ge bool foundSaddle = false; while (!foundSaddle) { - Vector z0 = normalize(randomVector(M.dimension(), d, r.engine())); + Vector z0 = normalize(randomVector(M.dimension(), d, r)); try { zSaddle = findSaddle(M, z0, ε); diff --git a/langevin.cpp b/langevin.cpp deleted file mode 100644 index ea28d65..0000000 --- a/langevin.cpp +++ /dev/null @@ -1,142 +0,0 @@ -#include -#include -#include -#include - -#include "Eigen/Dense" -#include "Eigen/src/Eigenvalues/ComplexEigenSolver.h" -#include "Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h" -#include "complex_normal.hpp" -#include "p-spin.hpp" -#include "dynamics.hpp" -#include "stokes.hpp" - -#include "collectStokesData.hpp" - -#include "pcg-cpp/include/pcg_random.hpp" -#include "randutils/randutils.hpp" -#include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h" - -#define PSPIN_P 3 -const int p = PSPIN_P; // polynomial degree of Hamiltonian -using Complex = std::complex; -using RealVector = Vector; -using RealMatrix = Matrix; -using RealTensor = Tensor; -using ComplexVector = Vector; -using ComplexMatrix = Matrix; -using ComplexTensor = Tensor; - -using Rng = randutils::random_generator; - -int main(int argc, char* argv[]) { - // model parameters - unsigned N = 10; // number of spins - 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 δ = 1; // threshold for determining saddle - Real Δ = 1e-3; - Real γ = 1e-1; // step size - unsigned t = 1000; // number of Langevin steps - unsigned M = 100; - unsigned n = 100; - - int opt; - - while ((opt = getopt(argc, argv, "N:M:n:T:e:r:i:g:t:E:")) != -1) { - switch (opt) { - case 'N': - N = (unsigned)atof(optarg); - break; - case 't': - t = (unsigned)atof(optarg); - break; - case 'T': - T = atof(optarg); - break; - case 'e': - δ = atof(optarg); - break; - case 'E': - ε = atof(optarg); - break; - case 'g': - γ = atof(optarg); - case 'r': - Rκ = atof(optarg); - break; - case 'i': - Iκ = atof(optarg); - break; - case 'n': - n = (unsigned)atof(optarg); - break; - case 'M': - M = (unsigned)atof(optarg); - break; - default: - exit(1); - } - } - - Rng r; - - collectStokesData<2, 4>(N, r.engine(), 1e-15, 1.0, 0.01); - - pSpinModelpSpin(N, r.engine(), 1, 0.01); - - std::normal_distribution Red(0, 1); - - RealVector zMin = randomMinimum(pSpin, Red, r.engine(), ε); - Real Hr; - RealVector dHr; - RealMatrix ddHr; - std::tie(Hr, dHr, ddHr, std::ignore) = pSpin.hamGradHess(zMin); - RealMatrix M1 = ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size(); - RealMatrix M2 = ddHr - zMin.dot(dHr) * Matrix::Identity(zMin.size(), zMin.size()) / Real(N); - Eigen::EigenSolver> eigenS(M2); - for (unsigned i = 0; i < N; i++) { - RealVector zNew = normalize(zMin + 1e-5 * eigenS.eigenvectors().col(i).real()); - std::cout << eigenS.eigenvalues()(i) << " " << 2.0 * (pSpin.getHamiltonian(zNew) - Hr) / 1e-10 << " " << 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();; - ComplexVector zSaddle = zMin.cast(); - - ComplexVector zSaddleNext; - bool foundSaddle = false; - while (!foundSaddle) { - ComplexVector z0 = normalize(zSaddle + δ * randomVector(N, d, r.engine())); - try { - zSaddleNext = findSaddle(complexPSpin, z0, ε); - Real saddleDistance = (zSaddleNext - zSaddle).norm(); - if (saddleDistance / N > 1e-2) { - std::cout << saddleDistance << std::endl; - foundSaddle = true; - } - } catch (std::exception& e) {} - } - - 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; - - complexPSpin *= exp(Complex(0, φ)); - - Cord test(complexPSpin, zSaddle, zSaddleNext, 3); - test.relaxNewton(10, 1, 1e4); - - std::cout << test.totalCost(10) << " " << test.totalCost(20) << std::endl; - - return 0; -} diff --git a/pureStokesFromMinima.cpp b/pureStokesFromMinima.cpp new file mode 100644 index 0000000..8f815cf --- /dev/null +++ b/pureStokesFromMinima.cpp @@ -0,0 +1,54 @@ +#include +#include +#include + +#include "collectStokesData.hpp" + +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" +#include "unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h" + +#define PSPIN_P 3 +const int p = PSPIN_P; // polynomial degree of Hamiltonian + +using Rng = randutils::random_generator; + +int main(int argc, char* argv[]) { + // model parameters + unsigned N = 10; // number of spins + // simulation parameters + Real ε = 1e-15; + Real δ = 1; + unsigned n = 10; + + int opt; + + while ((opt = getopt(argc, argv, "N:e:d:n:")) != -1) { + switch (opt) { + case 'N': + N = (unsigned)atof(optarg); + break; + case 'e': + ε = atof(optarg); + break; + case 'd': + δ = atof(optarg); + break; + case 'n': + n = atof(optarg); + break; + default: + exit(1); + } + } + + Rng r; + + for (unsigned i = 0; i < n; i++) { + auto tag = std::chrono::high_resolution_clock::now(); + std::ofstream output("stokes_" + std::to_string(tag.time_since_epoch().count()) + ".dat"); + collectStokesData<3>(output, N, r.engine(), ε, δ, true, 1.0); + } + + return 0; +} diff --git a/stokes.hpp b/stokes.hpp index 9fd6f01..8c0c1fb 100644 --- a/stokes.hpp +++ b/stokes.hpp @@ -1,9 +1,10 @@ #pragma once #include "p-spin.hpp" -#include "complex_normal.hpp" #include "dynamics.hpp" +#include + template class Cord { public: @@ -17,7 +18,7 @@ public: Scalar H2 = M.getHamiltonian(z2); Scalar H3 = M.getHamiltonian(z3); - if (real(H2) > real(H3)) { + if (std::real(H2) > std::real(H3)) { z0 = z2; z1 = z3; } else { @@ -62,7 +63,7 @@ public: Vector ż = zDot(z, dH); Vector dz = df(t); - return 1 - real(ż.dot(dz)) / ż.norm() / dz.norm(); + return 1 - std::real(ż.dot(dz)) / ż.norm() / dz.norm(); } Real totalCost(unsigned nt) const { @@ -117,24 +118,22 @@ public: Real dfdgn = dgCoeff(i, t); Vector dC = - 0.5 / ż.norm() / dz.norm() * ( dfdgn * ż.conjugate() + fdgn * dżc * dz + fdgn * dż * dz.conjugate() - - real(dz.dot(ż)) * ( + - std::real(dz.dot(ż)) * ( dfdgn * dz.conjugate() / dz.squaredNorm() + fdgn * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm() ) ); - for (unsigned j = 0; j < dC.size(); j++) { - x(z.size() * i + j) = dC(j); - x(z.size() * gs.size() + z.size() * i + j) = conj(dC(j)); - } + x.segment(z.size() * i, z.size()) = dC; + x.segment(z.size() * gs.size() + z.size() * i, z.size()) = dC.conjugate(); for (unsigned j = 0; j < gs.size(); j++) { Real fdgm = gCoeff(j, t); Real dfdgm = dgCoeff(j, t); - Scalar CC = - real(dz.dot(ż)) / ż.norm() / dz.norm(); + Scalar CC = - std::real(dz.dot(ż)) / ż.norm() / dz.norm(); Vector dCm = - 0.5 / ż.norm() / dz.norm() * ( dfdgm * ż.conjugate() + fdgm * dżc * dz + fdgm * dż * dz.conjugate() - - real(dz.dot(ż)) * ( + - std::real(dz.dot(ż)) * ( dfdgm * dz.conjugate() / dz.squaredNorm() + fdgm * (dżc * ż + dż * ż.conjugate()) / ż.squaredNorm() ) @@ -195,21 +194,17 @@ public: ) ; - for (unsigned k = 0; k < z.size(); k++) { - for (unsigned l = 0; l < z.size(); l++) { - M(z.size() * i + k, z.size() * j + l) = dcdC(k, l); - M(z.size() * gs.size() + z.size() * i + k, z.size() * j + l) = conj(ddC(k, l)); - M(z.size() * i + k, z.size() * gs.size() + z.size() * j + l) = ddC(k, l); - M(z.size() * gs.size() + z.size() * i + k, z.size() * gs.size() + z.size() * j + l) = conj(dcdC(k, l)); - } - } + M.block(z.size() * i, z.size() * j, z.size(), z.size()) = dcdC; + M.block(z.size() * gs.size() + z.size() * i, z.size() * j, z.size(), z.size()) = ddC.conjugate(); + M.block(z.size() * i, z.size() * gs.size() + z.size() * j, z.size(), z.size()) = ddC; + M.block(z.size() * gs.size() + z.size() * i, z.size() * gs.size() + z.size() * j, z.size(), z.size()) = dcdC.conjugate(); } } return {x, M}; } - std::pair relaxStepNewton(unsigned nt, Real δ₀) { + std::tuple 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()); @@ -231,32 +226,37 @@ public: Matrix I = abs(HessTot.diagonal().array()).matrix().asDiagonal(); + Vector δg(gs.size() * z0.size()); + while (newCost > oldCost) { - Vector δg = (HessTot + δ * I).partialPivLu().solve(dgsTot); + δg = (HessTot + δ * I).partialPivLu().solve(dgsTot).tail(gs.size() * z0.size()); for (unsigned i = 0; i < gs.size(); i++) { - for (unsigned j = 0; j < z0.size(); j++) { - cNew.gs[i](j) = gs[i](j) - δg[z0.size() * gs.size() + z0.size() * i + j]; - } + cNew.gs[i] = gs[i] - δg.segment(z0.size() * i, z0.size()); } newCost = cNew.totalCost(nt); - δ *= 2; + δ *= 1.5; } + std::cerr << newCost << " " << stepSize << " " << δ << std::endl; + gs = cNew.gs; - return {δ / 2, stepSize}; + return {δ / 1.5, stepSize}; } - void relaxNewton(unsigned nt, Real δ₀, unsigned maxSteps) { + void relaxNewton(unsigned nt, Real δ₀, Real ε, unsigned maxSteps) { Real δ = δ₀; Real stepSize = std::numeric_limits::infinity(); unsigned steps = 0; - while (stepSize > 1e-7 && steps < maxSteps) { + while (stepSize > ε && steps < maxSteps) { std::tie(δ, stepSize) = relaxStepNewton(nt, δ); - δ /= 3; + if (δ > 100) { + nt++; + } + δ /= 2; steps++; } } -- cgit v1.2.3-54-g00ecf