Many fixes.

demos/log_normal_wiener_filter.py

@@ -2,14 +2,15 @@
 
 from nifty import *
 
+
 if __name__ == "__main__":
     nifty_configuration['default_distribution_strategy'] = 'fftw'
 
     # Setting up parameters    |\label{code:wf_parameters}|
     correlation_length = 1.     # Typical distance over which the field is correlated
     field_variance = 2.         # Variance of field in position space
-    response_sigma = 0.05       # Smoothing length of response (in same unit as L)
-    signal_to_noise = 1.5       # The signal to noise ratio
+    response_sigma = 0.02       # Smoothing length of response (in same unit as L)
+    signal_to_noise = 100       # The signal to noise ratio
     np.random.seed(43)          # Fixing the random seed
     def power_spectrum(k):      # Defining the power spectrum
         a = 4 * correlation_length * field_variance**2
@@ -17,8 +18,9 @@ if __name__ == "__main__":
 
     # Setting up the geometry |\label{code:wf_geometry}|
     L = 2. # Total side-length of the domain
-    N_pixels = 512 # Grid resolution (pixels per axis)
+    N_pixels = 128 # Grid resolution (pixels per axis)
     signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
+    #signal_space = HPSpace(16)
     harmonic_space = FFTOperator.get_default_codomain(signal_space)
     fft = FFTOperator(harmonic_space, target=signal_space, target_dtype=np.float)
     power_space = PowerSpace(harmonic_space)
@@ -31,7 +33,7 @@ if __name__ == "__main__":
     # Setting up an exemplary response
     mask = Field(signal_space, val=1.)
     N10 = int(N_pixels/10)
-    mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
+    #mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
     R = ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
     data_domain = R.target[0]
     R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
@@ -47,15 +49,27 @@ if __name__ == "__main__":
     energy = library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S)
 
 
-    minimizer = VL_BFGS(convergence_tolerance=0,
-                        iteration_limit=50,
-                        #callback=convergence_measure,
-                        max_history_length=3)
+    minimizer1 = VL_BFGS(convergence_tolerance=1e-4,
+                         iteration_limit=1000,
+                         #callback=convergence_measure,
+                         max_history_length=3)
+
+    minimizer2 = RelaxedNewton(convergence_tolerance=1e-4,
+                               iteration_limit=1,
+                               #callback=convergence_measure
+                               )
+    minimizer3 = SteepestDescent(convergence_tolerance=1e-4, iteration_limit=1000)
+
+
+    me1 = minimizer1(energy)
+    me2 = minimizer2(energy)
+    me3 = minimizer3(energy)
+
+    m1 = fft(me1[0].position)
+    m2 = fft(me2[0].position)
+    m3 = fft(me3[0].position)
+
 
-    minimizer = RelaxedNewton(convergence_tolerance=0,
-                              iteration_limit=1,
-                              #callback=convergence_measure
-                              )
 
 #    # Probing the variance
 #    class Proby(DiagonalProberMixin, Prober): pass
@@ -67,18 +81,19 @@ if __name__ == "__main__":
 
     #Plotting #|\label{code:wf_plotting}|
     plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
+    #plotter = plotting.HealpixPlotter(color_map=plotting.colormaps.PlankCmap())
+
     plotter.figure.xaxis = plotting.Axis(label='Pixel Index')
     plotter.figure.yaxis = plotting.Axis(label='Pixel Index')
 
-    #plotter.plot.zmax = exp(mock_signal.max()); plotter.plot.zmin = 0
+    plotter.plot.zmax = 5; plotter.plot.zmin = -5
 #    plotter(variance, path = 'variance.html')
     #plotter.plot.zmin = exp(mock_signal.min());
-    plotter(mock_signal, path='mock_signal.html')
-    plotter(Field(signal_space, val=np.log(data.val.get_full_data()).reshape(signal_space.shape)),
-            path = 'data.html')
-#    plotter(m, path = 'map.html')
-
-
-
+    plotter(mock_signal.real, path='mock_signal.html')
+    plotter(Field(signal_space, val=np.log(data.val.get_full_data().real).reshape(signal_space.shape)),
+            path = 'log_of_data.html')
 
+    plotter(m1.real, path='m_LBFGS.html')
+    plotter(m2.real, path='m_Newton.html')
+    plotter(m3.real, path='m_SteepestDescent.html')
 

nifty/basic_arithmetics.py

@@ -23,7 +23,7 @@ from nifty.field import Field
 
 __all__ = ['cos', 'sin', 'cosh', 'sinh', 'tan', 'tanh', 'arccos', 'arcsin',
            'arccosh', 'arcsinh', 'arctan', 'arctanh', 'sqrt', 'exp', 'log',
-           'conjugate']
+           'conjugate', 'clipped_exp']
 
 
 def _math_helper(x, function):
@@ -95,6 +95,10 @@ def exp(x):
     return _math_helper(x, np.exp)
 
 
+def clipped_exp(x):
+    return _math_helper(x, lambda z: np.exp(np.minimum(200, z)))
+
+
 def log(x, base=None):
     result = _math_helper(x, np.log)
     if base is not None:

nifty/field.py

@@ -831,6 +831,24 @@ class Field(Loggable, Versionable, object):
         except TypeError:
             return 0.
 
+    @property
+    def real(self):
+        """ The real part of the field (data is not copied).
+        """
+        real_part = self.val.real
+        result = self.copy_empty(dtype=real_part.dtype)
+        result.set_val(new_val=real_part, copy=False)
+        return result
+
+    @property
+    def imag(self):
+        """ The imaginary part of the field (data is not copied).
+        """
+        real_part = self.val.imag
+        result = self.copy_empty(dtype=real_part.dtype)
+        result.set_val(new_val=real_part, copy=False)
+        return result
+
     # ---Special unary/binary operations---
 
     def cast(self, x=None, dtype=None):

nifty/library/log_normal_wiener_filter/log_normal_wiener_filter_curvature.py

 from nifty.operators import EndomorphicOperator,\
                             InvertibleOperatorMixin
 from nifty.energies.memoization import memo
-from nifty.basic_arithmetics import exp
+from nifty.basic_arithmetics import clipped_exp
 
 
 class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
@@ -56,15 +56,15 @@ class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
 
     def _times(self, x, spaces):
         part1 = self.S.inverse_times(x)
-        part2 = self._exppRNRexppd * x
+        # part2 = self._exppRNRexppd * x
         part3 = self._expp * self.R.adjoint_times(
                                 self.N.inverse_times(self.R(self._expp * x)))
-        return part1 + part2 + part3
+        return part1 + part3  # + part2
 
     @property
     @memo
     def _expp(self):
-        return exp(self.position)
+        return clipped_exp(self.position)
 
     @property
     @memo

nifty/library/log_normal_wiener_filter/log_normal_wiener_filter_energy.py

@@ -38,7 +38,7 @@ class LogNormalWienerFilterEnergy(Energy):
     @property
     @memo
     def value(self):
-        return 0.5*(self.position.vdot(self._Sp) -
+        return 0.5*(self.position.vdot(self._Sp) +
                     self._Rexppd.vdot(self._NRexppd))
 
     @property

nifty/minimization/descent_minimizer.py

@@ -149,17 +149,23 @@ class DescentMinimizer(Loggable, object):
             descent_direction = self.get_descent_direction(energy)
             # compute the step length, which minimizes energy.value along the
             # search direction
-            step_length, f_k, new_energy = \
-                self.line_searcher.perform_line_search(
-                                               energy=energy,
-                                               pk=descent_direction,
-                                               f_k_minus_1=f_k_minus_1)
+            try:
+                step_length, f_k, new_energy = \
+                    self.line_searcher.perform_line_search(
+                                                   energy=energy,
+                                                   pk=descent_direction,
+                                                   f_k_minus_1=f_k_minus_1)
+            except RuntimeError:
+                self.logger.warn(
+                        "Stopping because of RuntimeError in line-search")
+                break
+
             if f_k_minus_1 is None:
-                delta=1e30
+                delta = 1e30
             else:
-                delta = abs(f_k -f_k_minus_1)/max(abs(f_k),abs(f_k_minus_1),1.)
+                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
             if (new_energy.value - energy.value) > 0:
                 self.logger.info("Line search algorithm returned a new energy "

nifty/minimization/line_searching/line_search_strong_wolfe.py

@@ -111,7 +111,9 @@ class LineSearchStrongWolfe(LineSearch):
         le_0 = self.line_energy.at(0)
         phi_0 = le_0.value
         phiprime_0 = le_0.directional_derivative
-        assert phiprime_0<0, "input direction must be a descent direction"
+        if phiprime_0 >= 0:
+            self.logger.error("Input direction must be a descent direction")
+            raise RuntimeError
 
         # set alphas
         alpha0 = 0.
@@ -234,12 +236,14 @@ class LineSearchStrongWolfe(LineSearch):
         cubic_delta = 0.2  # cubic
         quad_delta = 0.1  # quadratic
         phiprime_alphaj = 0.
+        alpha_recent = None
+        phi_recent = None
 
         assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0
-        assert phiprime_lo*(alpha_hi-alpha_lo)<0.
+        assert phiprime_lo*(alpha_hi-alpha_lo) < 0.
         for i in xrange(max_iterations):
-            #assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0
-            #assert phiprime_lo*(alpha_hi-alpha_lo)<0.
+            # assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0
+            # assert phiprime_lo*(alpha_hi-alpha_lo)<0.
             delta_alpha = alpha_hi - alpha_lo
             if delta_alpha < 0:
                 a, b = alpha_hi, alpha_lo

nifty/plotting/plots/heatmaps/glmollweide.py

@@ -32,7 +32,9 @@ class GLMollweide(Heatmap):
                            xsize=self.xsize,
                            color_map=self.color_map,
                            webgl=self.webgl,
-                           smoothing=self.smoothing)
+                           smoothing=self.smoothing,
+                           zmin=self.zmin,
+                           zmax=self.zmax)
 
     @staticmethod
     def _find_closest(A, target):

nifty/plotting/plots/heatmaps/hpmollweide.py

@@ -30,7 +30,9 @@ class HPMollweide(Heatmap):
                            xsize=self.xsize,
                            color_map=self.color_map,
                            webgl=self.webgl,
-                           smoothing=self.smoothing)
+                           smoothing=self.smoothing,
+                           zmin=self.zmin,
+                           zmax=self.zmax)
 
     def _mollview(self, x):
         xsize = self.xsize

[Martin Reinecke]

I have a really bad feeling about this ...

The fact that such a function shows up is a strong indication that we are doing something incorrect somewhere else; I'd rather find the original error than use stop-gap measures like this.

This hack introduces non-differentiability and potentially vanishing gradients into the energy, which give the minimizers a hard time; strictly speaking, this breaks function properties required by our minimizers.

It also destroys the identity log(exp(a))==a.