From 136fcddcd38d0b8f3b40faf7c1cb7365d9b2a753 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Fri, 15 Jan 2021 14:45:19 +0100
Subject: Converted more of library to templates accepting generic Scalar types
 and p

---
 dynamics.hpp | 117 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 117 insertions(+)
 create mode 100644 dynamics.hpp

(limited to 'dynamics.hpp')

diff --git a/dynamics.hpp b/dynamics.hpp
new file mode 100644
index 0000000..85eef71
--- /dev/null
+++ b/dynamics.hpp
@@ -0,0 +1,117 @@
+#pragma once
+
+#include <exception>
+#include <eigen3/Eigen/LU>
+#include <random>
+
+#include "p-spin.hpp"
+#include "stereographic.hpp"
+
+class gradientDescentStallException: public std::exception {
+  virtual const char* what() const throw() {
+    return "Gradient descent stalled.";
+  }
+} gradientDescentStall;
+
+template <class Scalar, int p>
+std::tuple<double, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double ε, double γ0 = 1, double δγ = 2) {
+  Vector<Scalar> z = z0;
+  double γ = γ0;
+
+  auto [W, dW] = WdW(J, z);
+
+  while (W > ε) {
+    Vector<Scalar> zNewTmp = z - γ * dW.conjugate();
+    Vector<Scalar> zNew = normalize(zNewTmp);
+
+    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 gradientDescentStall;
+    }
+  }
+
+  return {W, z};
+}
+
+template <class Scalar, int p>
+Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double ε, double δW = 2, double γ0 = 1, double δγ = 2) {
+  Vector<Scalar> z = z0;
+  Vector<Scalar> ζ = euclideanToStereographic(z);
+
+  double W;
+  std::tie(W, std::ignore) = WdW(J, z);
+
+  Vector<Scalar> dH;
+  Matrix<Scalar> ddH;
+  std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z);
+
+  while (W > ε) {
+    // ddH is complex symmetric, which is (almost always) invertible, so a
+    // partial pivot LU decomposition can be used.
+    Vector<Scalar> dζ = ddH.partialPivLu().solve(dH);
+    Vector<Scalar> ζNew = ζ - dζ;
+    Vector<Scalar> zNew = stereographicToEuclidean(ζNew);
+
+    double WNew;
+    std::tie(WNew, std::ignore) = WdW(J, zNew);
+
+    if (WNew < W) { // If Newton's step lowered the objective, accept it!
+      ζ = ζNew;
+      z = zNew;
+      W = WNew;
+    } else { // Otherwise, do gradient descent until W is a factor δW smaller.
+      std::tie(W, z) = gradientDescent(J, z, W / δW, γ0, δγ);
+      ζ = euclideanToStereographic(z);
+    }
+
+    std::tie(std::ignore, dH, ddH) = stereographicHamGradHess(J, ζ, z);
+  }
+
+  return z;
+}
+
+template <class Scalar, class Distribution, class Generator>
+Vector<Scalar> randomVector(unsigned N, Distribution d, Generator& r) {
+  Vector<Scalar> z(N);
+
+  for (unsigned i = 0; i < N; i++) {
+    z(i) = d(r);
+  }
+
+  return z;
+}
+
+template <class Scalar, int p, class Distribution, class Generator>
+std::tuple<double, Vector<Scalar>> langevin(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, double T, double γ, unsigned N, Distribution d, Generator& r) {
+  Vector<Scalar> z = z0;
+
+  double W;
+  std::tie(W, std::ignore) = WdW(J, z);
+
+  std::uniform_real_distribution<double> D(0, 1);
+
+  for (unsigned i = 0; i < N; i++) {
+    Vector<Scalar> zNewTmp = z + randomVector<Scalar>(z.size(), d, r);
+    Vector<Scalar> zNew = normalize(zNewTmp);
+
+    double WNew;
+    std::tie(WNew, std::ignore) = WdW(J, zNew);
+
+    if (exp((W - WNew) / T) > D(r)) {
+      z = zNew;
+      W = WNew;
+    }
+  }
+
+  return {W, z};
+}
-- 
cgit v1.2.3-54-g00ecf