Removed some unused files.

5 files changed, 9 insertions, 210 deletions

diff --git a/dynamics.hpp b/dynamics.hpp
index 21afd30..7d04b1a 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -4,45 +4,7 @@
 #include <eigen3/Eigen/LU>
 #include <random>
 
-#include <iostream>
-
 #include "p-spin.hpp"
-#include "stereographic.hpp"
-
-class gradientDescentStallException: public std::exception {
-  virtual const char* what() const throw() {
-    return "Gradient descent stalled.";
-  }
-};
-
-template <class Real, class Scalar, int p>
-std::tuple<Real, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, Real ε, Real γ0 = 1, Real δγ = 2) {
-  Vector<Scalar> z = z0;
-  Real γ = γ0;
-
-  auto [W, dW] = WdW(J, z);
-
-  while (W > ε) {
-    Vector<Scalar> 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 <class Real, int p>
 Vector<Real> findMinimum(const Tensor<Real, p>& J, const Vector<Real>& z0, Real ε) {
@@ -98,8 +60,6 @@ Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& 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<Scalar>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size()  + Matrix<Scalar>::Identity(z.size(), z.size()) * (zz - (Real)z.size()) + 2.0 * z * z.transpose();
@@ -119,30 +79,6 @@ Vector<Scalar> randomVector(unsigned N, Distribution d, Generator& r) {
   return z;
 }
 
-template <class Real, class Scalar, int p, class Distribution, class Generator>
-std::tuple<Real, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0,
-    std::function<Real(const Tensor<Scalar, p>&, const Vector<Scalar>&)>& energy,
-    Real T, Real γ, unsigned N, Distribution d, Generator& r) {
-  Vector<Scalar> z = z0;
-
-  Real E = energy(J, z);
-
-  std::uniform_real_distribution<Real> D(0, 1);
-
-  for (unsigned i = 0; i < N; i++) {
-    Vector<Scalar> zNew = normalize(z + γ * randomVector<Scalar>(z.size(), d, r));
-
-    Real ENew = energy(J, zNew);
-
-    if (E - ENew > T * log(D(r))) {
-      z = zNew;
-      E = ENew;
-    }
-  }
-
-  return {E, z};
-}
-
 template <class Real, int p, class Distribution, class Generator>
 Vector<Real> randomMinimum(const Tensor<Real, p>& J, Distribution d, Generator& r, Real ε) {
   Vector<Real> 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<Real>;
+using RealVector = Vector<Real>;
+using RealMatrix = Matrix<Real>;
+using RealTensor = Tensor<Real, p>;
 using ComplexVector = Vector<Complex>;
 using ComplexMatrix = Matrix<Complex>;
 using ComplexTensor = Tensor<Complex, p>;
@@ -144,12 +147,12 @@ int main(int argc, char* argv[]) {
 
   std::normal_distribution<Real> Red(0, 1);
 
-  Vector<Real> zMin = randomMinimum(ReJ, Red, r, ε);
+  RealVector zMin = randomMinimum(ReJ, Red, r, ε);
   auto [Hr, dHr, ddHr] = hamGradHess(ReJ, zMin);
   Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::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<Real> 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<Real> d(0, 1, 0);
 
-  ComplexTensor J = generateCouplings<Complex, p>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine());
-
-  std::function<Real(const ComplexTensor&, const ComplexVector&)> 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<Complex>();
+  ComplexVector zSaddle = zMin.cast<Complex>();
 
   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<Real>(J, zSaddle, zSaddleNext, 10, 4)) {
-    std::cerr << "Found a Stokes line" << std::endl;
-//    stokesLineTestNew<Real>(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 <eigen3/Eigen/Cholesky>
-
-#include "p-spin.hpp"
-
-template <class Scalar>
-Vector<Scalar> stereographicToEuclidean(const Vector<Scalar>& ζ) {
-  unsigned N = ζ.size() + 1;
-  Vector<Scalar> 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 <class Scalar>
-Vector<Scalar> euclideanToStereographic(const Vector<Scalar>& z) {
-  unsigned N = z.size();
-  Vector<Scalar> ζ(N - 1);
-
-  for (unsigned i = 0; i < N - 1; i++) {
-    ζ(i) = z(i);
-  }
-
-  return ζ / (sqrt((Real)N) - z(N - 1));
-}
-
-template <class Scalar>
-Matrix<Scalar> stereographicJacobian(const Vector<Scalar>& ζ) {
-  unsigned N = ζ.size();
-  Matrix<Scalar> 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 <class Scalar, int p>
-std::tuple<Scalar, Vector<Scalar>, Matrix<Scalar>> stereographicHamGradHess(const Tensor<Scalar, p>& J, const Vector<Scalar>& ζ, const Vector<Scalar>& z) {
-  auto [hamiltonian, gradZ, hessZ] = hamGradHess(J, z);
-  Matrix<Scalar> jacobian = stereographicJacobian(ζ);
-
-  Matrix<Scalar> metric = jacobian * jacobian.adjoint();
-
-  // The metric is Hermitian and positive definite, so a Cholesky decomposition can be used.
-  Vector<Scalar> grad = metric.llt().solve(jacobian) * gradZ;
-  Matrix<Scalar> 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<Vector<Scalar>, Matrix<Scalar>> gsGradHess(const Tensor<Scalar, p>& J, Real t) const {
     Vector<Scalar> z = f(t);
     auto [H, dH, ddH] = hamGradHess(J, z);
-    Tensor<Scalar, 3> dddH = ((p - 2) * (p - 1) * p / (Real)factorial(p)) * J;
+    Tensor<Scalar, 3> dddH = ((p - 2) * (p - 1) * p / (Real)factorial(p)) * J; // BUG IF P != 3
     Vector<Scalar> ż = zDot(z, dH);
     Tensor<Scalar, 1> żT = Eigen::TensorMap<Tensor<const Scalar, 1>>(ż.data(), {z.size()});
     Vector<Scalar> 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<std::complex<Real>> z(2), dz(2), ddz(2), dH(2);
-  Matrix<std::complex<Real>> 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;
-}