From 3276bdd1e9796fec71e169e6c41d77da72b3a4fb Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 25 Feb 2021 15:28:11 +0100
Subject: Many changes.

---
 dynamics.hpp | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

(limited to 'dynamics.hpp')

diff --git a/dynamics.hpp b/dynamics.hpp
index d421d13..10b1be2 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -13,10 +13,10 @@ class gradientDescentStallException: public std::exception {
   }
 };
 
-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) {
+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;
-  double γ = γ0;
+  Real γ = γ0;
 
   auto [W, dW] = WdW(J, z);
 
@@ -42,12 +42,12 @@ std::tuple<double, Vector<Scalar>> gradientDescent(const Tensor<Scalar, p>& J, c
   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) {
+template <class Real, class Scalar, int p>
+Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0, Real ε, Real δW = 2, Real γ0 = 1, Real δγ = 2) {
   Vector<Scalar> z = z0;
   Vector<Scalar> ζ = euclideanToStereographic(z);
 
-  double W;
+  Real W;
   std::tie(W, std::ignore) = WdW(J, z);
 
   Vector<Scalar> dH;
@@ -61,7 +61,7 @@ Vector<Scalar> findSaddle(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0,
     Vector<Scalar> ζNew = ζ - dζ;
     Vector<Scalar> zNew = stereographicToEuclidean(ζNew);
 
-    double WNew;
+    Real WNew;
     std::tie(WNew, std::ignore) = WdW(J, zNew);
 
     if (WNew < W) { // If Newton's step lowered the objective, accept it!
@@ -90,20 +90,20 @@ Vector<Scalar> randomVector(unsigned N, Distribution d, Generator& r) {
   return z;
 }
 
-template <class Scalar, int p, class Distribution, class Generator>
-std::tuple<double, Vector<Scalar>> metropolis(const Tensor<Scalar, p>& J, const Vector<Scalar>& z0,
-    std::function<double(const Tensor<Scalar, p>&, const Vector<Scalar>&)>& energy,
-    double T, double γ, unsigned N, Distribution d, Generator& r) {
+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;
 
-  double E = energy(J, z);
+  Real E = energy(J, z);
 
-  std::uniform_real_distribution<double> D(0, 1);
+  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));
 
-    double ENew = energy(J, zNew);
+    Real ENew = energy(J, zNew);
 
     if (E - ENew > T * log(D(r))) {
       z = zNew;
-- 
cgit v1.2.3-54-g00ecf