Playing with new step adjuster

1 files changed, 20 insertions, 1 deletions

diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp
index 0e05366..a061e1f 100644
--- a/log-fourier_integrator.cpp
+++ b/log-fourier_integrator.cpp
@@ -133,17 +133,26 @@ int main(int argc, char* argv[]) {
     Real ΔCmin = 1000;
     Real ΔCₜ = 100;
     unsigned stepsUp = 0;
+    Real cost;
+    Real oldCost = 1000;
     while (ΔCₜ > ε) {
       auto [RddfCt, dfCt] = RddfCtdfCt(fft, Cₜ, Rₜ, p, s, λ);
 
       std::vector<Complex> Ĉₜ₊₁(N);
       std::vector<Complex> Ȓₜ₊₁(N);
+      cost = 0;
       for (unsigned n = 0; n < N; n++) {
+        cost += std::norm((μₜ + II * fft.ν(n)) * Ȓₜ[n] - ((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]));
+        cost += std::norm((μₜ + II * fft.ν(n)) * Ĉₜ[n] - ((2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n])))));
         Ȓₜ₊₁[n] = ((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n));
         Ĉₜ₊₁[n] = ((2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n))).real();
       }
+      cost = sqrt(cost);
+      cost /= 2*N;
+
       std::vector<Real> Rₜ₊₁ = fft.inverse(Ȓₜ₊₁);
       std::vector<Real> Cₜ₊₁ = fft.inverse(Ĉₜ₊₁);
+      cost += std::abs(Cₜ₊₁[0] - 1);
 
       Real C₀ = Cₜ₊₁[0];
 
@@ -202,6 +211,15 @@ int main(int argc, char* argv[]) {
       }
       ΔCₜ = sqrt(ΔCₜ) / (2*N);
 
+      if (cost < oldCost) {
+        γ = std::min(1.01 * γ, γ₀);
+      } else {
+        γ = std::max(γ / 2, (Real)1e-2);;
+      }
+
+      oldCost = cost;
+
+      /*
       if (ΔCₜ < 0.9 * ΔCmin) {
         ΔCmin = ΔCₜ;
         stepsUp = 0;
@@ -214,6 +232,7 @@ int main(int argc, char* argv[]) {
         stepsUp = 0;
         ΔCmin = ΔCₜ;
       }
+      */
 
       for (unsigned i = 0; i < N; i++) {
         Cₜ[i] += γ * (Cₜ₊₁[i] - Cₜ[i]);
@@ -223,7 +242,7 @@ int main(int argc, char* argv[]) {
       }
 
       std::cerr << "\x1b[2K" << "\r";
-      std::cerr << β << " " << μₜ << " " << ΔCₜ << " " << γ << " " << C₀;
+      std::cerr << β << " " << μₜ << " " << ΔCₜ << " " << γ << " " << C₀ << " " << cost;
     }
 
     if (std::isnan(Cₜ[0])) {