1 files changed, 16 insertions, 79 deletions

diff --git a/langevin.cpp b/langevin.cpp
index 86aeb9e..dc3d07d 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -27,63 +27,6 @@ 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, Real δ) {
-  std::list<std::array<ComplexVector, 2>> pairs;
-
-  for (auto it1 = saddles.begin(); it1 != saddles.end(); it1++) {
-    for (auto it2 = std::next(it1); it2 != saddles.end(); it2++) {
-      if ((it1->second - it2->second).norm() < δ) {
-        pairs.push_back({it1->second, it2->second});
-      }
-    }
-  }
-
-  return pairs;
-}
-
-template <class Generator>
-std::tuple<ComplexTensor, ComplexVector, ComplexVector> matchImaginaryEnergies(const ComplexTensor& J0, const ComplexVector& z10, const ComplexVector& z20, Real ε, Real Δ, Generator& r) {
-  Real σ = sqrt(factorial(p) / (Real(2) * pow(z10.size(), p - 1)));
-  complex_normal_distribution<Real> dJ(0, σ, 0);
-
-  ComplexTensor J = J0;
-  ComplexVector z1, z2;
-  Real ΔH = abs(imag(getHamiltonian(J, z10) - getHamiltonian(J, z20))) / z10.size();
-  Real γ = ΔH;
-
-  std::function<void(ComplexTensor&, std::array<unsigned, p>)> perturbJ =
-    [&γ, &dJ, &r] (ComplexTensor& JJ, std::array<unsigned, p> ind) {
-      Complex Ji = getJ<Complex, p>(JJ, ind);
-      setJ<Complex, p>(JJ, ind, Ji + γ * dJ(r.engine()));
-    };
-
-  while (ΔH > ε) {
-    ComplexTensor Jp = J;
-    iterateOver<Complex, p>(Jp, perturbJ);
-
-    try {
-      z1 = findSaddle(Jp, z10, Δ);
-      z2 = findSaddle(Jp, z20, Δ);
-
-      Real Δz = (z1 - z2).norm();
-
-      if (Δz > 1e-2) {
-        Real ΔHNew = abs(imag(getHamiltonian(Jp, z1) - getHamiltonian(Jp, z2))) / z1.size();
-
-        if (ΔHNew < ΔH) {
-          J = Jp;
-          ΔH = ΔHNew;
-          γ = ΔH;
-
-          std::cerr << "Match imaginary energies: Found couplings with ΔH = " << ΔH << std::endl;
-        }
-      }
-    } catch (std::exception& e) {}
-  }
-
-  return {J, z1, z2};
-}
-
 int main(int argc, char* argv[]) {
   // model parameters
   unsigned N = 10; // number of spins
@@ -143,24 +86,29 @@ int main(int argc, char* argv[]) {
 
   Rng r;
 
-  Tensor<Real, p> ReJ = generateCouplings<Real, p>(N, std::normal_distribution<Real>(0, σ), r.engine());
+  pSpinModel<Real, p> pSpin;
+
+  std::get<0>(pSpin.Js) = generateCouplings<Real, p>(N, std::normal_distribution<Real>(0, σ), r.engine());
 
   std::normal_distribution<Real> Red(0, 1);
 
-  RealVector zMin = randomMinimum(ReJ, Red, r, ε);
-  auto [Hr, dHr, ddHr] = hamGradHess(ReJ, zMin);
+  RealVector zMin = randomMinimum(pSpin, Red, r, ε);
+  Real Hr;
+  RealVector dHr;
+  RealMatrix ddHr;
+  std::tie(Hr, dHr, ddHr, std::ignore) = pSpin.hamGradHess(zMin);
   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++) {
     RealVector 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 << pSpin.getHamiltonian(zNew) - Hr << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << std::endl;
   }
   std::cout << std::endl;
   getchar();
 
   complex_normal_distribution<Real> d(0, 1, 0);
 
-  ComplexTensor J = ReJ.cast<Complex>();
+  pSpinModel<Complex, p> complexPSpin = pSpin.cast<Complex>();;
   ComplexVector zSaddle = zMin.cast<Complex>();
 
   ComplexVector zSaddleNext;
@@ -168,7 +116,7 @@ int main(int argc, char* argv[]) {
   while (!foundSaddle) {
     ComplexVector z0 = normalize(zSaddle + δ * randomVector<Complex>(N, d, r.engine()));
     try {
-      zSaddleNext = findSaddle(J, z0, ε);
+      zSaddleNext = findSaddle(complexPSpin, z0, ε);
       Real saddleDistance = (zSaddleNext - zSaddle).norm();
       if (saddleDistance / N > 1e-2) {
         std::cout << saddleDistance << std::endl;
@@ -177,27 +125,16 @@ int main(int argc, char* argv[]) {
     } catch (std::exception& e) {}
   }
 
-  auto [H1, dH1, ddH1] = hamGradHess(J, zSaddle);
-  auto [H2, dH2, ddH2] = hamGradHess(J, zSaddleNext);
-
-  Eigen::ComplexEigenSolver<ComplexMatrix> ces;
-  ces.compute(ddH1);
+  Complex H1 = complexPSpin.getHamiltonian(zSaddle);
+  Complex H2 = complexPSpin.getHamiltonian(zSaddleNext);
 
   Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());
   std::cerr << (zSaddle - zSaddleNext).norm() / (Real)N << " " << φ << " " << H1 * exp(Complex(0, φ)) << " " << H2 * exp(Complex(0, φ)) << std::endl;
-  Real smallestNorm = std::numeric_limits<Real>::infinity();
-  for (Complex z : ces.eigenvalues()) {
-    if (norm(z) < smallestNorm) {
-      smallestNorm = norm(z);
-    }
-  }
-  std::cerr << smallestNorm << std::endl;
-  getchar();
 
-  J = exp(Complex(0, φ)) * J;
+  std::get<0>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<0>(complexPSpin.Js);
 
-  Cord test(J, zSaddle, zSaddleNext, 3);
-  test.relaxNewton(J, 10, 1, 1e4);
+  Cord test(complexPSpin, zSaddle, zSaddleNext, 3);
+  test.relaxNewton(10, 1, 1e4);
 
   std::cout << test.z0.transpose() << std::endl;
   std::cout << test.z1.transpose() << std::endl;