From d27bca16e0c3d6a6fba795d91f7084caeffd0f05 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 14 Jul 2017 16:03:27 +0200
Subject: [PATCH] new stopping criterion

---
 demos/critical_filtering.py                   | 29 ++++++++++---------
 nifty/energies/line_energy.py                 |  5 ++--
 nifty/minimization/descent_minimizer.py       |  5 +++-
 .../test_descent_minimizers.py                |  3 +-
 4 files changed, 24 insertions(+), 18 deletions(-)

diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py
index 5674ec83..11356f17 100644
--- a/demos/critical_filtering.py
+++ b/demos/critical_filtering.py
@@ -76,8 +76,8 @@ if __name__ == "__main__":
 
 
     # Choosing the measurement instrument
-    Instrument = SmoothingOperator(s_space, sigma=0.01)
-    # Instrument = DiagonalOperator(s_space, diagonal=1.)
+    # Instrument = SmoothingOperator(s_space, sigma=0.01)
+    Instrument = DiagonalOperator(s_space, diagonal=1.)
     # Instrument._diagonal.val[200:400, 200:400] = 0
     #Instrument._diagonal.val[64:512-64, 64:512-64] = 0
 
@@ -100,8 +100,8 @@ if __name__ == "__main__":
     realized_power = log(sh.power_analyze(binbounds=p_space.binbounds))
     data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds))
     d_data = d.val.get_full_data().real
-    #if rank == 0:
-    #    pl.plot([go.Heatmap(z=d_data)], filename='data.html')
+    if rank == 0:
+        pl.plot([go.Heatmap(z=d_data)], filename='data.html')
 
     #  minimization strategy
 
@@ -110,17 +110,18 @@ if __name__ == "__main__":
         print (x, iteration)
 
 
-    minimizer1 = RelaxedNewton(convergence_tolerance=1e-2,
-                              convergence_level=2,
+    minimizer1 = RelaxedNewton(convergence_tolerance=1e-8,
+                              convergence_level=1,
                               iteration_limit=5,
                               callback=convergence_measure)
 
-    minimizer2 = VL_BFGS(convergence_tolerance=1e-3,
-                       iteration_limit=20,
+    minimizer2 = VL_BFGS(convergence_tolerance=1e-8,
+                       convergence_level=1,
+                       iteration_limit=1000,
                        callback=convergence_measure,
-                       max_history_length=10)
-    minimizer3 = SteepestDescent(convergence_tolerance=1e-3,
-                       iteration_limit=70,
+                       max_history_length=20)
+    minimizer3 = SteepestDescent(convergence_tolerance=1e-8,
+                       iteration_limit=500,
                        callback=convergence_measure)
 
     # Setting starting position
@@ -143,7 +144,7 @@ if __name__ == "__main__":
         # Initializing the power energy with updated parameters
         power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, smoothness_prior=10., samples=3)
 
-        (power_energy, convergence) = minimizer3(power_energy)
+        (power_energy, convergence) = minimizer2(power_energy)
 
 
         # Setting new power spectrum
@@ -151,7 +152,7 @@ if __name__ == "__main__":
 
         # Plotting current estimate
         print i
-        #if i%50 == 0:
-        #    plot_parameters(m0,t0,log(sp), data_power)
+        if i%50 == 0:
+            plot_parameters(m0,t0,log(sp), data_power)
 
 
diff --git a/nifty/energies/line_energy.py b/nifty/energies/line_energy.py
index 3f9c0f75..e816203e 100644
--- a/nifty/energies/line_energy.py
+++ b/nifty/energies/line_energy.py
@@ -108,6 +108,7 @@ class LineEnergy:
     @property
     def dd(self):
         res = self.energy.gradient.vdot(self.linedir)
-        assert abs(res-res.real)<1e-12, \
-               "directional derivative has non-negligible imaginary part"
+        if abs(res.imag)/max(abs(res.real),1.)>1e-12:
+               print "directional derivative has non-negligible " \
+                     "imaginary part:", res
         return res.real
diff --git a/nifty/minimization/descent_minimizer.py b/nifty/minimization/descent_minimizer.py
index be3abb53..10436faf 100644
--- a/nifty/minimization/descent_minimizer.py
+++ b/nifty/minimization/descent_minimizer.py
@@ -154,6 +154,10 @@ class DescentMinimizer(Loggable, object):
                                                energy=energy,
                                                pk=descent_direction,
                                                f_k_minus_1=f_k_minus_1)
+            if f_k_minus_1 is None:
+                delta=1e30
+            else:
+                delta = abs(f_k -f_k_minus_1)/max(abs(f_k),abs(f_k_minus_1),1.)
             f_k_minus_1 = energy.value
             tx1=energy.position-new_energy.position
             # check if new energy value is bigger than old energy value
@@ -165,7 +169,6 @@ class DescentMinimizer(Loggable, object):
 
             energy = new_energy
             # check convergence
-            delta = abs(gradient).max() * step_length / gradient_norm
             self.logger.debug("Iteration:%08u step_length=%3.1E "
                               "delta=%3.1E energy=%3.1E" %
                               (iteration_number, step_length, delta,
diff --git a/test/test_minimization/test_descent_minimizers.py b/test/test_minimization/test_descent_minimizers.py
index fcffa73d..81daf7eb 100644
--- a/test/test_minimization/test_descent_minimizers.py
+++ b/test/test_minimization/test_descent_minimizers.py
@@ -43,7 +43,8 @@ class Test_DescentMinimizers(unittest.TestCase):
         covariance = DiagonalOperator(space, diagonal=covariance_diagonal)
         energy = QuadraticPotential(position=starting_point,
                                     eigenvalues=covariance)
-        minimizer = minimizer_class(iteration_limit=30)
+        minimizer = minimizer_class(iteration_limit=30,
+                                    convergence_tolerance=1e-10)
 
         (energy, convergence) = minimizer(energy)
 
-- 
GitLab