A lot of work, and fixed a huge bug regarding the meaning of .dot in the Eigen library for complex vectors.

1 files changed, 66 insertions, 16 deletions

diff --git a/langevin.cpp b/langevin.cpp
index 9da4ef6..3bc12ed 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -9,6 +9,7 @@
 
 #include "complex_normal.hpp"
 #include "p-spin.hpp"
+#include "stereographic.hpp"
 
 using Rng = randutils::random_generator<pcg32>;
 
@@ -23,6 +24,57 @@ Vector randomVector(unsigned N, Distribution d, Generator& r) {
   return z;
 }
 
+Vector findSaddle(const Tensor& J, const Vector& z0, double γ0, double δ, double ε, Rng& r) {
+  Vector z = z0;
+
+  double W;
+  Vector dW;
+  std::tie(W, dW) = WdW(J, z);
+
+  while (W > δ) {
+    double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2);
+
+    Vector zNew = z - γ * dW.conjugate();
+    zNew *= sqrt(z.size()) / sqrt((Scalar)(zNew.transpose() * zNew));
+
+    double WNew;
+    Vector dWNew;
+    std::tie(WNew, dWNew) = WdW(J, zNew);
+
+    if (WNew < W) {
+      z = zNew;
+      W = WNew;
+      dW = dWNew;
+      std::cout << W << std::endl;
+    } 
+  }
+
+  Vector ζ = euclideanToStereographic(z);
+  Vector dH;
+  Matrix ddH;
+  std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ);
+
+  while (W > ε) {
+    double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2);
+
+    Vector ζNew = ζ - ddH.ldlt().solve(dH);
+
+    double WNew;
+    Vector dWNew;
+    std::tie(WNew, dWNew) = WdW(J, stereographicToEuclidean(ζNew));
+
+    if (WNew < W) {
+      ζ = ζNew;
+      W = WNew;
+      dW = dWNew;
+      std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ);
+      std::cout << WNew << std::endl;
+    } 
+  }
+
+  return stereographicToEuclidean(ζ);
+}
+
 Vector langevin(const Tensor& J, const Vector& z0, double T, double γ0,
                 std::function<bool(double, unsigned)> quit, Rng& r) {
   Vector z = z0;
@@ -32,24 +84,16 @@ Vector langevin(const Tensor& J, const Vector& z0, double T, double γ0,
   std::tie(W, dW) = WdW(J, z);
 
   unsigned steps = 0;
-  complex_normal_distribution<> d(0, T, 0);
+  complex_normal_distribution<> d(0, sqrt(T), 0);
 
   while (!quit(W, steps)) {
     double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2);
     Vector η = randomVector(z.size(), d, r.engine());
 
-    Vector zNext = z - γ * dW + η;
-    zNext *= sqrt(zNext.size()) / sqrt(zNext.dot(zNext));
+    z -= γ * (dW.conjugate() / pow((double)z.size(), 2) + η);
+    z *= sqrt(z.size()) / sqrt((Scalar)(z.transpose() * z));
 
-    double WNext;
-    Vector dWNext;
-    std::tie(WNext, dWNext) = WdW(J, zNext);
-
-    if (exp((W - WNext) / T) > r.uniform(0.0, 1.0)) {
-      z = zNext;
-      W = WNext;
-      dW = dWNext;
-    }
+    std::tie(W, dW) = WdW(J, z);
   }
 
   return z;
@@ -63,13 +107,14 @@ int main(int argc, char* argv[]) {
   double Iκ = 0;   // imaginary part of distribution parameter
 
   // simulation parameters
+  double ε = 1e-4;
   double δ = 1e-2;   // threshold for determining saddle
   double γ = 1e-2;   // step size
   unsigned t = 1000; // number of Langevin steps
 
   int opt;
 
-  while ((opt = getopt(argc, argv, "N:T:e:r:i:g:t:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:T:e:r:i:g:t:E:")) != -1) {
     switch (opt) {
     case 'N':
       N = (unsigned)atof(optarg);
@@ -83,6 +128,9 @@ int main(int argc, char* argv[]) {
     case 'e':
       δ = atof(optarg);
       break;
+    case 'E':
+      ε = atof(optarg);
+      break;
     case 'g':
       γ = atof(optarg);
     case 'r':
@@ -103,14 +151,15 @@ int main(int argc, char* argv[]) {
 
   Tensor J = generateCouplings<Scalar, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine());
   Vector z = randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine());
-  z *= sqrt(N) / sqrt(z.dot(z)); // Normalize.
+  z *= sqrt(N) / sqrt((Scalar)(z.transpose() * z)); // Normalize.
 
-  std::function<bool(double, unsigned)> findSaddle = [δ](double W, unsigned) {
+  std::function<bool(double, unsigned)> f = [δ](double W, unsigned) {
     std::cout << W << std::endl;
     return W < δ;
   };
 
-  Vector zm = langevin(J, z, T, γ, findSaddle, r);
+  //Vector zm = langevin(J, z, T, γ, f, r);
+  Vector zm = findSaddle(J, z, γ, δ, ε, r);
 
   Scalar H;
   Vector dH;
@@ -119,5 +168,6 @@ int main(int argc, char* argv[]) {
 
   Vector constraint = dH - ((double)p * H / (double)N) * zm;
 
+  std::cout << constraint << std::endl;
   std::cout << H / (double)N << std::endl;
 }