From 6b42ce01cd289567a6109d831fe86a0667d18f37 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Mon, 8 Nov 2021 23:05:15 +0100
Subject: Working minimization.

---
 dynamics.hpp | 8 +++-----
 langevin.cpp | 7 ++++++-
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/dynamics.hpp b/dynamics.hpp
index 639803b..21afd30 100644
--- a/dynamics.hpp
+++ b/dynamics.hpp
@@ -55,13 +55,13 @@ Vector<Real> findMinimum(const Tensor<Real, p>& J, const Vector<Real>& z0, Real
   Matrix<Real> M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose();
 
   while (g.norm() / z.size() > ε && λ < 1e8) {
-    Vector<Real> dz = (M + λ * Matrix<Real>::Identity(z.size(), z.size())).partialPivLu().solve(g);
+    Vector<Real> dz = (M + λ * (Matrix<Real>)abs(M.diagonal().array()).matrix().asDiagonal()).partialPivLu().solve(g);
     dz -= z.dot(dz) * z / (Real)z.size();
     Vector<Real> zNew = normalize(z - dz);
 
     auto [HNew, dHNew, ddHNew] = hamGradHess(J, zNew);
 
-    if (HNew <= H * 1.01) {
+    if (HNew * 1.0001 <= H) {
       z = zNew;
       H = HNew;
       dH = dHNew;
@@ -70,12 +70,10 @@ Vector<Real> findMinimum(const Tensor<Real, p>& J, const Vector<Real>& z0, Real
       g = dH - z.dot(dH) * z / (Real)z.size();
       M = ddH - (dH * z.transpose() + z.dot(dH) * Matrix<Real>::Identity(z.size(), z.size()) + (ddH * z) * z.transpose()) / (Real)z.size() + 2.0 * z * z.transpose();
 
-      λ /= 1.001;
+      λ /= 2;
     } else {
       λ *= 1.5;
     }
-
-    std::cout << "error : " << H << " " << g.norm() / z.size() << " " << λ << std::endl;
   }
 
   return z;
diff --git a/langevin.cpp b/langevin.cpp
index a775774..1181173 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -146,8 +146,13 @@ int main(int argc, char* argv[]) {
 
   Vector<Real> zMin = randomMinimum(ReJ, Red, r, ε);
   auto [Hr, dHr, ddHr] = hamGradHess(ReJ, zMin);
-  Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - ((ddHr * zMin) * zMin.transpose()) / (Real)zMin.size());
+  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());
+    std::cout << getHamiltonian(ReJ, zNew) - Hr << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << std::endl;
+  }
+  std::cout << std::endl;
   getchar();
 
   complex_normal_distribution<Real> d(0, 1, 0);
-- 
cgit v1.2.3-54-g00ecf