1 files changed, 31 insertions, 27 deletions

diff --git a/topology.cpp b/topology.cpp
index 6650ea6..1d46649 100644
--- a/topology.cpp
+++ b/topology.cpp
@@ -103,15 +103,40 @@ public:
     return {∂L, ∂∂L};
   }
 
-  Vector newtonMethod(const Vector& v0, Real γ = 1, Real ε = 1e-12) const {
+  Vector newtonMethod(const Vector& v0, Real γ = 1, Real ε = 1e-12, unsigned maxSteps = 1e2) const {
     Vector v = v0;
+    unsigned steps = 0;
 
     Vector ∂L;
     Matrix ∂∂L;
 
     while (std::tie(∂L, ∂∂L) = ∂L_∂∂L(v), ∂L.squaredNorm() > ε) {
+      if (v.tail(M + 1).squaredNorm() > 1 / ε || steps > N * maxSteps) {
+        return Vector::Zero(N + M + 1);
+      }
+
       v -= γ * ∂∂L.partialPivLu().solve(∂L);
-      v.head(N) = normalize(v.head(N)); // might as well stay on the sphere
+
+      v.head(N) = normalize(v.head(N)); // Stay on the sphere
+      steps++;
+    }
+
+    return v;
+  }
+
+  Vector randomStationaryPoint(Rng& r, Real γ = 1, Real ε = 1e-12, unsigned maxTries = 1e1) const {
+    Vector v = Vector::Zero(N + M + 1);
+    unsigned tries = 0;
+
+    while (v.squaredNorm() < ε && tries < maxTries) {
+      Vector v0(N + M + 1);
+      for (Real& v0ᵢ : v0) {
+        v0ᵢ = r.variate<Real, std::normal_distribution>();
+      }
+      v0.head(N) = normalize(v0.head(N));
+
+      v = newtonMethod(v0, γ, ε);
+      tries++;
     }
 
     return v;
@@ -235,30 +260,15 @@ int main(int argc, char* argv[]) {
   std::cout << std::setprecision(15);
 
   if (quadratic) {
-    Vector x = Vector::Zero(N);
-    x(0) = sqrt(N);
-
     for (unsigned sample = 0; sample < samples; sample++) {
-      Vector v0(N + M + 1);
-      v0.head(N) = x;
-      for (Real& v0ᵢ : v0.tail(M + 1)) {
-        v0ᵢ = r.variate<Real, std::normal_distribution>();
-      }
       ExtensiveQuadratic S(N, M, E, r);
-      Vector v = S.newtonMethod(v0, γ);
-      std::cout << S.overlap(v) << " " << S.cost(v) << std::endl;
+      Vector v = S.randomStationaryPoint(r, γ);
+      std::cout << S.overlap(v) << std::endl;
     }
   } else if (!count) {
-    Vector x = Vector::Zero(N);
-    x(0) = sqrt(N);
-
     for (unsigned sample = 0; sample < samples; sample++) {
-      Vector v0(N + 2);
-      v0 << x,
-        r.variate<Real, std::normal_distribution>(),
-        r.variate<Real, std::normal_distribution>();
       Spherical3Spin S(N, E, r);
-      Vector v = S.newtonMethod(v0, γ);
+      Vector v = S.randomStationaryPoint(r, γ);
       std::cout << S.overlap(v) << std::endl;
     }
   } else {
@@ -266,13 +276,7 @@ int main(int argc, char* argv[]) {
     std::list<Vector> points;
 
     for (unsigned sample = 0; sample < samples; sample++) {
-      Vector v0(N + 2);
-      for (Real& v0i : v0) {
-        v0i = r.variate<Real, std::normal_distribution>();
-      }
-      v0.head(N) = normalize(v0.head(N));
-
-      Vector v = S.newtonMethod(v0, γ);
+      Vector v = S.randomStationaryPoint(r, γ);
       if (v(N) == v(N)) {
         bool inList = false;
         for (const Vector& u : points) {