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. --- langevin.cpp | 142 ----------------------------------------------------------- 1 file changed, 142 deletions(-) delete mode 100644 langevin.cpp (limited to 'langevin.cpp') 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; -} -- cgit v1.2.3-54-g00ecf