From 54416465850a9280a121127e0b5301d97ad13a61 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 8 Nov 2021 23:19:54 +0100 Subject: Removed some unused files. --- dynamics.hpp | 64 -------------------------------------------------- langevin.cpp | 58 +++++++-------------------------------------- stereographic.hpp | 70 ------------------------------------------------------- stokes.hpp | 2 +- stokes_test.cpp | 25 -------------------- 5 files changed, 9 insertions(+), 210 deletions(-) delete mode 100644 stereographic.hpp delete mode 100644 stokes_test.cpp diff --git a/dynamics.hpp b/dynamics.hpp index 21afd30..7d04b1a 100644 --- a/dynamics.hpp +++ b/dynamics.hpp @@ -4,45 +4,7 @@ #include #include -#include - #include "p-spin.hpp" -#include "stereographic.hpp" - -class gradientDescentStallException: public std::exception { - virtual const char* what() const throw() { - return "Gradient descent stalled."; - } -}; - -template -std::tuple> gradientDescent(const Tensor& J, const Vector& z0, Real ε, Real γ0 = 1, Real δγ = 2) { - Vector z = z0; - Real γ = γ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 gradientDescentStallException(); - } - } - - return {W, z}; -} template Vector findMinimum(const Tensor& J, const Vector& z0, Real ε) { @@ -98,8 +60,6 @@ Vector findSaddle(const Tensor& J, const Vector& z0, std::tie(std::ignore, dH, ddH) = hamGradHess(J, z); - std::cout << "error : " << g.norm() / z.size() << std::endl; - 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(); @@ -119,30 +79,6 @@ Vector randomVector(unsigned N, Distribution d, Generator& r) { return z; } -template -std::tuple> metropolis(const Tensor& J, const Vector& z0, - std::function&, const Vector&)>& energy, - Real T, Real γ, unsigned N, Distribution d, Generator& r) { - Vector z = z0; - - Real E = energy(J, z); - - std::uniform_real_distribution D(0, 1); - - for (unsigned i = 0; i < N; i++) { - Vector zNew = normalize(z + γ * randomVector(z.size(), d, r)); - - Real ENew = energy(J, zNew); - - if (E - ENew > T * log(D(r))) { - z = zNew; - E = ENew; - } - } - - return {E, z}; -} - template Vector randomMinimum(const Tensor& J, Distribution d, Generator& r, Real ε) { Vector zSaddle; diff --git a/langevin.cpp b/langevin.cpp index 1181173..9aec7c8 100644 --- a/langevin.cpp +++ b/langevin.cpp @@ -18,6 +18,9 @@ #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; @@ -144,12 +147,12 @@ int main(int argc, char* argv[]) { std::normal_distribution Red(0, 1); - Vector zMin = randomMinimum(ReJ, Red, r, ε); + RealVector zMin = randomMinimum(ReJ, Red, r, ε); auto [Hr, dHr, ddHr] = hamGradHess(ReJ, 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++) { - Vector zNew = normalize(zMin + 0.01 * eigenS.eigenvectors().col(i).real()); + 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 << std::endl; @@ -157,17 +160,8 @@ int main(int argc, char* argv[]) { complex_normal_distribution d(0, 1, 0); - ComplexTensor J = generateCouplings(N, complex_normal_distribution(0, σ, κ), r.engine()); - - std::function energyNormGrad = [] - (const ComplexTensor& J, const ComplexVector& z) { - Real W; - std::tie(W, std::ignore) = WdW(J, z); - return W; - }; - - ComplexVector zSaddle = randomSaddle(J, d, r, ε); - std::cerr << "Found saddle." << std::endl; + ComplexTensor J = ReJ.cast(); + ComplexVector zSaddle = zMin.cast(); ComplexVector zSaddleNext; bool foundSaddle = false; @@ -198,46 +192,10 @@ int main(int argc, char* argv[]) { } } std::cerr << smallestNorm << std::endl; + getchar(); J = exp(Complex(0, φ)) * J; - /* - if (stokesLineTest(J, zSaddle, zSaddleNext, 10, 4)) { - std::cerr << "Found a Stokes line" << std::endl; -// stokesLineTestNew(J, zSaddle, zSaddleNext, 10, 3); - } else { - std::cerr << "Didn't find a Stokes line" << std::endl; - } - */ - - /* - J(0,0,0) = Complex(2, 3); - J(1,1,1) = Complex(-2, 0.3); - J(0,1,1) = Complex(4, 0.2); - J(1,0,1) = Complex(4, 0.2); - J(1,1,0) = Complex(4, 0.2); - J(1,0,0) = Complex(0.7, 0.4); - J(0,1,0) = Complex(0.7, 0.4); - J(0,0,1) = Complex(0.7, 0.4); - - ComplexVector z0(2);; - z0 << Complex(0.8, 0.3), Complex(0.7, 0.2); - ComplexVector z1(2); - z1 << Complex(-0.5, 0.2), Complex(1.0, 1.0); - - Cord test(J, z0, z1, 2); - - test.gs[0](0) = Complex(0.2, 0.2); - test.gs[0](1) = Complex(0.4, 0.4); - test.gs[1](0) = Complex(0.1, 0.2); - test.gs[1](1) = Complex(0.3, 0.4); - - auto [dgs, ddgs] = test.gsGradHess(J, 0.7); - - std::cout << dgs << std::endl; - std::cout << ddgs << std::endl; - */ - Cord test(J, zSaddle, zSaddleNext, 5); test.relaxNewton(J, 25, 1, 1e4); diff --git a/stereographic.hpp b/stereographic.hpp deleted file mode 100644 index 58c6b31..0000000 --- a/stereographic.hpp +++ /dev/null @@ -1,70 +0,0 @@ -#pragma once - -#include - -#include "p-spin.hpp" - -template -Vector stereographicToEuclidean(const Vector& ζ) { - unsigned N = ζ.size() + 1; - Vector z(N); - - Scalar a = ζ.transpose() * ζ; - Scalar b = 2 * sqrt((Real)N) / ((Real)1 + a); - - for (unsigned i = 0; i < N - 1; i++) { - z(i) = ζ(i); - } - - z(N - 1) = (a - (Real)1) / (Real)2; - - return b * z; -} - -template -Vector euclideanToStereographic(const Vector& z) { - unsigned N = z.size(); - Vector ζ(N - 1); - - for (unsigned i = 0; i < N - 1; i++) { - ζ(i) = z(i); - } - - return ζ / (sqrt((Real)N) - z(N - 1)); -} - -template -Matrix stereographicJacobian(const Vector& ζ) { - unsigned N = ζ.size(); - Matrix J(N, N + 1); - - Scalar b = (Real)1 + (Scalar)(ζ.transpose() * ζ); - - for (unsigned i = 0; i < N; i++) { - for (unsigned j = 0; j < N; j++) { - J(i, j) = - ζ(i) * ζ(j); - - if (i == j) { - J(i, j) += b / (Real)2; - } - } - - J(i, N) = ζ(i); - } - - return 4 * sqrt(N + 1) * J / pow(b, 2); -} - -template -std::tuple, Matrix> stereographicHamGradHess(const Tensor& J, const Vector& ζ, const Vector& z) { - auto [hamiltonian, gradZ, hessZ] = hamGradHess(J, z); - Matrix jacobian = stereographicJacobian(ζ); - - Matrix metric = jacobian * jacobian.adjoint(); - - // The metric is Hermitian and positive definite, so a Cholesky decomposition can be used. - Vector grad = metric.llt().solve(jacobian) * gradZ; - Matrix hess = jacobian * hessZ * jacobian.transpose(); - - return {hamiltonian, grad, hess}; -} diff --git a/stokes.hpp b/stokes.hpp index e14db52..cf68d9b 100644 --- a/stokes.hpp +++ b/stokes.hpp @@ -336,7 +336,7 @@ public: 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; + Tensor dddH = ((p - 2) * (p - 1) * p / (Real)factorial(p)) * J; // BUG IF P != 3 Vector ż = zDot(z, dH); Tensor żT = Eigen::TensorMap>(ż.data(), {z.size()}); Vector dz = df(t); diff --git a/stokes_test.cpp b/stokes_test.cpp deleted file mode 100644 index a9347a7..0000000 --- a/stokes_test.cpp +++ /dev/null @@ -1,25 +0,0 @@ -#include "stokes.hpp" - -using namespace std::complex_literals; - -int main() { - - Vector> z(2), dz(2), ddz(2), dH(2); - Matrix> ddH(2, 2); - - z << -0.75067 - 0.190132 * 1i, -0.625994 + 0.665987 * 1i; - dz << -0.0636149 - 0.469166 * 1i, 0.820037 - 0.449064 * 1i; - ddz << 0.55777 + 0.730164 * 1i, 0.361959 - 0.463245 * 1i; - dH << 0.967613 - 0.907519 * 1i, 0.712336 - 0.649056 * 1i; - - ddH << -0.371925 - 0.280788 * 1i, -0.163888 - 0.141297 * 1i, - -0.163888 - 0.141297 * 1i, 0.230969 + 0.942449 * 1i; - - std::cout << zDot(z, dH) << std::endl << std::endl; - - std::cout << segmentCost(z, dz, dH) << std::endl << std::endl; - - std::cout << variation(z, dz, ddz, dH, ddH) << std::endl; - - return 0; -} -- cgit v1.2.3-70-g09d2