Removed external μ control

2 files changed, 80 insertions, 139 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;
diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index b15808b..531565c 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -93,108 +93,73 @@ int main(int argc, char* argv[]) {
   Real fac = 1;
   Real β = 0;
   while (β < βₘₐₓ) {
-    Real μ₁ = 0;
-    Real μ₂ = 0;
-    while (true) {
-      Real ΔC = 100;
-      Real ΔC₀ = 100;
-      unsigned it = 0;
-      while (ΔC > ε) {
-        std::vector<Real> dfC(N);
-        std::vector<Real> RddfC(N);
-        for (unsigned n = 0; n < N; n++) {
-          RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]);
-          dfC[n] = df(λ, p, s, Cₜ[n]);
-        }
-        std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
-        std::vector<Complex> dfCt = fft.fourier(dfC, true);
+    Real ΔC = 100;
+    Real ΔC₀ = 100;
+    unsigned it = 0;
+    while (ΔC > ε) {
+      std::vector<Real> dfC(N);
+      std::vector<Real> RddfC(N);
+      for (unsigned n = 0; n < N; n++) {
+        RddfC[n] = Rₜ[n] * ddf(λ, p, s, Cₜ[n]);
+        dfC[n] = df(λ, p, s, Cₜ[n]);
+      }
+      std::vector<Complex> RddfCt = fft.fourier(RddfC, false);
+      std::vector<Complex> dfCt = fft.fourier(dfC, true);
 
-        std::vector<Complex> Ȓₜ₊₁(N);
-        std::vector<Complex> Ĉₜ₊₁(N);
+      std::vector<Complex> Ȓₜ₊₁(N);
+      std::vector<Complex> Ĉₜ₊₁(N);
 
-        for (unsigned n = 0; n < N; n++) {
-          Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n));
+      for (unsigned n = 0; n < N; n++) {
+        Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μ + 1i * fft.ν(n));
 //          Ĉₜ₊₁[n] = - 2 * Γ₀ * Ȓₜ₊₁[n].imag() / (1 + pow(τ₀ * fft.ν(n), 2)) / fft.ν(n);
-          Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
-        }
-
-        std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
-
-        for (unsigned n = 0; n < N; n++) {
-          RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]);
-        }
-        RddfCt = fft.fourier(RddfC, false);
-        for (unsigned n = 0; n < N; n++) {
-          Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
-        }
-        std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
-
-        ΔC = 0;
-        for (unsigned i = 0; i < N; i++) {
-          ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]);
-          ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]);
-        }
-        ΔC = sqrt(ΔC) / (2*N);
+        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
+      }
 
-        for (unsigned i = 0; i < N; i++) {
-          Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
-          Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
-          Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
-          Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
-        }
+      std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
 
-        /*
-        if (ΔC < ΔC₀) {
-          ΔC₀ = ΔC;
-          it = 0;
-          γ = std::min(1.001 * γ, 1.0);
-        } else {
-          it++;
-        }
+      for (unsigned n = 0; n < N; n++) {
+        RddfC[n] = Rₜ₊₁[n] * ddf(λ, p, s, Cₜ[n]);
+      }
+      RddfCt = fft.fourier(RddfC, false);
+      for (unsigned n = 0; n < N; n++) {
+        Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ₊₁[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ₊₁[n]))) / (μ + 1i * fft.ν(n));
+      }
+      std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
 
-        if (it > 100) {
-          γ = std::max(0.5 * γ, 1e-3);
-          it = 0;
-          ΔC₀ = ΔC;
-        }
-        */
+      μ *= pow(tanh(Cₜ₊₁[0]-1)+1, 0.05);
 
-        std::cerr << β << " " << μ << " " << ΔC << " " << γ << " " << Cₜ[0];
-        std::cerr << "\r";
+      ΔC = 0;
+      for (unsigned i = 0; i < N; i++) {
+        ΔC += std::norm(Ĉₜ[i] - Ĉₜ₊₁[i]);
+        ΔC += std::norm(Ȓₜ[i] - Ȓₜ₊₁[i]);
       }
-      if (std::isnan(Cₜ[0])) {
-        Cₜ = Cₜ₋₁;
-        Rₜ = Rₜ₋₁;
-        Ĉₜ = Ĉₜ₋₁;
-        Ȓₜ = Ȓₜ₋₁;
-        μ *= 1.1;
-        fac /= 2;
-        μ₁ = 0;
-        μ₂ = 0;
-      } else {
-      if (pow(Cₜ[0] - 1, 2) < ε) {
-        break;
+      ΔC = sqrt(ΔC) / (2*N);
+
+      for (unsigned i = 0; i < N; i++) {
+        Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
+        Rₜ[i] += γ * (Rₜ₊₁[i] - Rₜ[i]);
+        Ĉₜ[i] += γ * (Ĉₜ₊₁[i] - Ĉₜ[i]);
+        Ȓₜ[i] += γ * (Ȓₜ₊₁[i] - Ȓₜ[i]);
       }
-      if (μ₁ == 0 || μ₂ == 0) {
-        if (Cₜ[0] > 1 && μ₁ == 0) {
-          /* We found a lower bound */
-          μ₁ = μ;
-        }
-        if (Cₜ[0] < 1 && μ₂ == 0) {
-          /* We found an upper bound */
-          μ₂ = μ;
-        }
-        μ *= fac*std::tanh(Cₜ[0]-1)+1;
+
+      /*
+      if (ΔC < ΔC₀) {
+        ΔC₀ = ΔC;
+        it = 0;
+        γ = std::min(1.001 * γ, 1.0);
       } else {
-        /* Once the bounds are set, we can use bisection */
-        if (Cₜ[0] > 1) {
-          μ₁ = μ;
-        } else {
-          μ₂ = μ;
-        }
-        μ = (μ₁ + μ₂) / 2;
+        it++;
       }
+
+      if (it > 100) {
+        γ = std::max(0.5 * γ, 1e-3);
+        it = 0;
+        ΔC₀ = ΔC;
       }
+      */
+
+      std::cerr << "\x1b[2K" << "\r";
+      std::cerr << β << " " << μ << " " << ΔC << " " << γ << " " << Cₜ[0];
     }
 
     /* Integrate the energy using Simpson's rule */
@@ -213,9 +178,9 @@ int main(int argc, char* argv[]) {
     }
     E *= β;
 
+      std::cerr << "\x1b[2K" << "\r";
     std::cerr << β << " " << μ << " " << Ĉₜ[0].real() << " " << E << " " << γ << std::endl;
     β += Δβ;
-    μ *= 2;
     Cₜ₋₁ = Cₜ;
     Rₜ₋₁ = Rₜ;
     Ĉₜ₋₁ = Ĉₜ;