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.

---
 langevin.cpp | 56 +++++++++++++++++++++++++++-----------------------------
 1 file changed, 27 insertions(+), 29 deletions(-)

(limited to 'langevin.cpp')

diff --git a/langevin.cpp b/langevin.cpp
index 51cb2a0..8c191f3 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -14,14 +14,14 @@
 
 #define PSPIN_P 3
 const int p = PSPIN_P; // polynomial degree of Hamiltonian
-using Complex = std::complex<double>;
+using Complex = std::complex<Real>;
 using ComplexVector = Vector<Complex>;
 using ComplexMatrix = Matrix<Complex>;
 using ComplexTensor = Tensor<Complex, p>;
 
 using Rng = randutils::random_generator<pcg32>;
 
-std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_map<uint64_t, ComplexVector>& saddles, double δ) {
+std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_map<uint64_t, ComplexVector>& saddles, Real δ) {
   std::list<std::array<ComplexVector, 2>> pairs;
 
   for (auto it1 = saddles.begin(); it1 != saddles.end(); it1++) {
@@ -36,17 +36,17 @@ std::list<std::array<ComplexVector, 2>> saddlesCloserThan(const std::unordered_m
 }
 
 template <class Generator>
-std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, double ε, double Δ, Generator& r) {
-  double σ = sqrt(factorial(p) / (2.0 * pow(z10.size(), p - 1)));
-  complex_normal_distribution<> dJ(0, σ, 0);
+std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, Real ε, Real Δ, Generator& r) {
+  Real σ = sqrt(factorial(p) / (2.0 * pow(z10.size(), p - 1)));
+  complex_normal_distribution<Real> dJ(0, σ, 0);
 
   ComplexTensor J = J0;
   Complex H1, H2;
   ComplexVector z1, z2;
   std::tie(H1, std::ignore, std::ignore) = hamGradHess(J, z10);
   std::tie(H2, std::ignore, std::ignore) = hamGradHess(J, z20);
-  double prevdist = abs(imag(H1-H2));
-  double γ = 0.1 * prevdist;
+  Real prevdist = abs(imag(H1-H2));
+  Real γ = 0.1 * prevdist;
 
   std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ =
     [&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) {
@@ -62,7 +62,7 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
       z1 = findSaddle(Jp, z10, Δ);
       z2 = findSaddle(Jp, z20, Δ);
 
-      double dist = (z1 - z2).norm();
+      Real dist = (z1 - z2).norm();
 
       std::tie(H1, std::ignore, std::ignore) = hamGradHess(Jp, z1);
       std::tie(H2, std::ignore, std::ignore) = hamGradHess(Jp, z2);
@@ -87,16 +87,16 @@ std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(c
 int main(int argc, char* argv[]) {
   // model parameters
   unsigned N = 10; // number of spins
-  double T = 1;    // temperature
-  double Rκ = 1;   // real part of distribution parameter
-  double Iκ = 0;   // imaginary part of distribution parameter
+  Real T = 1;    // temperature
+  Real Rκ = 1;   // real part of distribution parameter
+  Real Iκ = 0;   // imaginary part of distribution parameter
 
   // simulation parameters
-  double ε = 1e-12;
-  double εJ = 1e-5;
-  double δ = 1e-2;   // threshold for determining saddle
-  double Δ = 1e-3;
-  double γ = 1e-2;   // step size
+  Real ε = 1e-12;
+  Real εJ = 1e-5;
+  Real δ = 1e-2;   // threshold for determining saddle
+  Real Δ = 1e-3;
+  Real γ = 1e-2;   // step size
   unsigned t = 1000; // number of Langevin steps
   unsigned M = 100;
   unsigned n = 100;
@@ -140,18 +140,18 @@ int main(int argc, char* argv[]) {
   }
 
   Complex κ(Rκ, Iκ);
-  double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1)));
+  Real σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1)));
 
   Rng r;
 
-  complex_normal_distribution<> d(0, 1, 0);
+  complex_normal_distribution<Real> d(0, 1, 0);
 
-  ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine());
+  ComplexTensor J = generateCouplings<Complex, PSPIN_P>(N, complex_normal_distribution<Real>(0, σ, κ), r.engine());
   ComplexVector z = normalize(randomVector<Complex>(N, d, r.engine()));
 
-  std::function<double(const ComplexTensor&, const ComplexVector&)> energyNormGrad = []
+  std::function<Real(const ComplexTensor&, const ComplexVector&)> energyNormGrad = []
     (const ComplexTensor& J, const ComplexVector& z) {
-      double W;
+      Real W;
       std::tie(W, std::ignore) = WdW(J, z);
       return W;
     };
@@ -185,18 +185,16 @@ int main(int argc, char* argv[]) {
   ComplexVector z1 = nearbySaddles.front()[0];
   ComplexVector z2 = nearbySaddles.front()[1];
 
-  for (unsigned j = 2; j < 8; j++) {
-    std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, pow(10, -2.0 * j), ε, r);
+  std::tie(J, z1, z2) = matchImaginaryEnergies(J, z1, z2, 1e-14, ε, r);
 
-    Rope<Complex> stokes(10, z1, z2, J);
+  Rope stokes(10, z1, z2, J);
 
-    for (unsigned i = 0; i < 8; i++) {
-      stokes.relax(J, 1e5, 0.1);
+  for (unsigned i = 0; i < 9; i++) {
+    stokes.relaxDiscreteGradient(J, 1e6, 0.1, 0);
 
-      std::cout << stokes.n() << " " << stokes.cost(J) << " " << stokes.length() << " " << stokes.error(J) << std::endl;
+    std::cout << stokes.n() << " " << stokes.cost(J) << std::endl;
 
-      stokes = stokes.interpolate();
-    }
+    stokes = stokes.interpolate();
   }
 
   return 0;
-- 
cgit v1.2.3-70-g09d2