1 files changed, 24 insertions, 48 deletions

diff --git a/fourier_integrator.cpp b/fourier_integrator.cpp
index 9ce47b6..96fbffd 100644
--- a/fourier_integrator.cpp
+++ b/fourier_integrator.cpp
@@ -120,13 +120,10 @@ int main(int argc, char* argv[]) {
   std::vector<Real> Rb = R;
   std::vector<Complex> Ctb = Ct;
   std::vector<Complex> Rtb = Rt;
+  Real μb = μ;
 
   Real fac = 1;
-  while (β += Δβ, β <= βₘₐₓ) {
-    Real Δμ = 1e-2;
-    Real μ₁ = 0;
-    Real μ₂ = 0;
-    while (true) {
+  while (β <= βₘₐₓ) {
     Real ΔC = 1;
     Real ΔCprev = 1000;
     unsigned it = 0;
@@ -150,6 +147,8 @@ int main(int argc, char* argv[]) {
         Rnew[i] = 0;
       }
 
+      μ *= pow(tanh(C[0]-1)+1, 0.05);
+
       ΔC = 0;
       for (unsigned i = 0; i < Cnew.size() / 2; i++) {
         ΔC += pow(Cnew[i] - C[i], 2);
@@ -164,7 +163,6 @@ int main(int argc, char* argv[]) {
         R[i] += γ * (Rnew[i] - R[i]);
       }
 
-      /*
       if (it % maxIterations == 0) {
         if (ΔC < ΔCprev) {
           γ = std::min(1.1 * γ, 1.0);
@@ -174,52 +172,23 @@ int main(int argc, char* argv[]) {
 
         ΔCprev = ΔC;
       }
-      */
 
+      std::cerr << "\x1b[2K" << "\r";
       std::cerr << μ << " " << p << " " << s << " " << τ₀ << " " << β << " " << sqrt(2 * ΔC / C.size()) << " " << γ << " " << C[0];
-      std::cerr << "\r";
-
-    }
-      if (std::isnan(C[0])) {
-        C = Cb;
-        R = Rb;
-        Ct = Ctb;
-        Rt = Rtb;
-//        μ /= sqrt(sqrt(fac*std::tanh(Cb[0]-1)+1));
-        μ *= 1.1;
-        fac /= 2;
-        μ₁ = 0;
-        μ₂ = 0;
-      } else {
-        Cb = C;
-        Rb = R;
-        Ctb = Ct;
-        Rtb = Rt;
-      if (pow(C[0] - 1, 2) < ε) {
-        break;
-      }
-      if (μ₁ == 0 || μ₂ == 0) {
-        if (C[0] > 1 && μ₁ == 0) {
-          /* We found a lower bound */
-          μ₁ = μ;
-        }
-        if (C[0] < 1 && μ₂ == 0) {
-          /* We found an upper bound */
-          μ₂ = μ;
-        }
-        μ *= sqrt(sqrt(fac*std::tanh(C[0]-1)+1));
-      } else {
-        /* Once the bounds are set, we can use bisection */
-        if (C[0] > 1) {
-          μ₁ = μ;
-        } else {
-          μ₂ = μ;
-        }
-        μ = (μ₁ + μ₂) / 2;
-      }
-      }
     }
 
+    if (std::isnan(C[0])) {
+      Δβ /= 2;
+      β -= Δβ;
+      C = Cb;
+      R = Rb;
+      Ct = Ctb;
+      Rt = Rtb;
+      μ = μb;
+    } else {
+
+    std::cerr << "\x1b[2K" << "\r";
+
     Real e = energy(C, R, p, s, λ, β, Δτ);
 
     std::cerr << "y " << β << " " << e << " " << μ << std::endl;
@@ -231,6 +200,13 @@ int main(int argc, char* argv[]) {
     std::ofstream outfileR(fourierFile("R", p, s, λ, τ₀, β, log2n, τₘₐₓ), std::ios::out | std::ios::binary);
     outfileR.write((const char*)(R.data()), (R.size() / 2) * sizeof(Real));
     outfileR.close();
+    β += Δβ;
+    Cb = C;
+    Rb = R;
+    Rtb = Rt;
+    Ctb = Ct;
+    μb = μ;
+    }
   }
 
   return 0;