diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py
index 3fef05a59b312239234f49f074096253e73a4b5d..7866f04d3e290f71b5457effbb2eaf69c3f7362e 100644
--- a/demos/critical_filtering.py
+++ b/demos/critical_filtering.py
@@ -53,7 +53,7 @@ if __name__ == "__main__":
     d_data = d.val
     ift.plot(d, name="data.png")
 
-    IC1 = ift.GradientNormController(verbose=True, iteration_limit=100,
+    IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
                                      tol_abs_gradnorm=0.1)
     minimizer = ift.RelaxedNewton(IC1)
 
diff --git a/demos/log_normal_wiener_filter.py b/demos/log_normal_wiener_filter.py
index 87563debbfae8a130727a7f2aaa9673ee30fb759..5ab513539d01bb135d5bd95d8c959245cfe73317 100644
--- a/demos/log_normal_wiener_filter.py
+++ b/demos/log_normal_wiener_filter.py
@@ -47,9 +47,8 @@ if __name__ == "__main__":
 
     # Wiener filter
     m0 = ift.Field.zeros(harmonic_space)
-    ctrl = ift.GradientNormController(verbose=False, tol_abs_gradnorm=0.0001)
-    ctrl2 = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1,
-                                       name="outer")
+    ctrl = ift.GradientNormController(tol_abs_gradnorm=0.0001)
+    ctrl2 = ift.GradientNormController(tol_abs_gradnorm=0.1, name="outer")
     inverter = ift.ConjugateGradient(controller=ctrl)
     energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic,
                                                      N, S, inverter=inverter)
diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py
index 7feef7420d90a66f0c4903e5aea20e3fb0e65031..9a67ca09ab4ff326f14e3407e63933bb8da889da 100644
--- a/demos/nonlinear_critical_filter.py
+++ b/demos/nonlinear_critical_filter.py
@@ -72,13 +72,12 @@ if __name__ == "__main__":
     t0 = ift.Field(p_space, val=-4.)
     power0 = Projection.adjoint_times(ift.exp(0.5 * t0))
 
-    IC1 = ift.GradientNormController(verbose=True, iteration_limit=100,
+    IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
                                      tol_abs_gradnorm=1e-3)
     LS = ift.LineSearchStrongWolfe(c2=0.02)
     minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
 
-    ICI = ift.GradientNormController(verbose=False, name="ICI",
-                                     iteration_limit=500,
+    ICI = ift.GradientNormController(iteration_limit=500,
                                      tol_abs_gradnorm=1e-3)
     inverter = ift.ConjugateGradient(controller=ICI)
 
diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py
index 4fd696d093a656d10289b7347925e6c376fb49c5..9b4665424be2334a6fa34b622085da6f4ed34000 100644
--- a/demos/paper_demos/cartesian_wiener_filter.py
+++ b/demos/paper_demos/cartesian_wiener_filter.py
@@ -86,7 +86,7 @@ if __name__ == "__main__":
 
     # Wiener filter
     j = R_harmonic.adjoint_times(N.inverse_times(data))
-    ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
+    ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
     inverter = ift.ConjugateGradient(controller=ctrl)
     wiener_curvature = ift.library.WienerFilterCurvature(
         S=S, N=N, R=R_harmonic, inverter=inverter)
diff --git a/demos/paper_demos/wiener_filter.py b/demos/paper_demos/wiener_filter.py
index ddf5a63a74468b782716bb5b21fc60d327046309..b0d682fcd671bb58b34c5b164eab7625661a5166 100644
--- a/demos/paper_demos/wiener_filter.py
+++ b/demos/paper_demos/wiener_filter.py
@@ -19,7 +19,7 @@ if __name__ == "__main__":
 
     signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
     harmonic_space = signal_space.get_default_codomain()
-    fft = ift.FFTOperator(harmonic_space, target=signal_space)
+    fft = ift.HarmonicTransformOperator(harmonic_space, target=signal_space)
     power_space = ift.PowerSpace(
         harmonic_space, binbounds=ift.PowerSpace.useful_binbounds(
             harmonic_space, logarithmic=True))
@@ -28,40 +28,47 @@ if __name__ == "__main__":
     S = ift.create_power_operator(harmonic_space,
                                   power_spectrum=power_spectrum)
     mock_power = ift.PS_field(power_space, power_spectrum)
-    mock_signal = fft(ift.power_synthesize(mock_power, real_signal=True))
+    mock_signal = ift.power_synthesize(mock_power, real_signal=True)
 
     # Setting up an exemplary response
     mask = np.ones(signal_space.shape)
     N10 = int(N_pixels/10)
     mask[N10*5:N10*9, N10*5:N10*9] = 0.
     mask = ift.Field(signal_space, ift.dobj.from_global_data(mask))
-    R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
-                             exposure=(mask,))
+    R = ift.GeometryRemover(signal_space)
+    R = R*ift.DiagonalOperator(mask)
+    R = R*fft
+    R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
     data_domain = R.target[0]
-    R_harmonic = R * fft
 
     # Setting up the noise covariance and drawing a random noise realization
-    ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
-    N = ift.DiagonalOperator(ndiag.weight(1))
+    ndiag = 1e-8*ift.Field.full(data_domain, fft(mock_signal).var()/signal_to_noise)
+    N = ift.DiagonalOperator(ndiag)
     noise = ift.Field.from_random(
         domain=data_domain, random_type='normal',
         std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
     data = R(mock_signal) + noise
 
     # Wiener filter
-    j = R_harmonic.adjoint_times(N.inverse_times(data))
-    ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
+    j = R.adjoint_times(N.inverse_times(data))
+    ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-6)
     inverter = ift.ConjugateGradient(controller=ctrl)
     wiener_curvature = ift.library.WienerFilterCurvature(
-        S=S, N=N, R=R_harmonic, inverter=inverter)
+        S=S, N=N, R=R, inverter=inverter)
     m_k = wiener_curvature.inverse_times(j)
     m = fft(m_k)
 
+    plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
+                "colormap": "Planck-like"}
+    ift.plot(mock_signal, name="mock_signal.png", **plotdict)
+    ift.plot(ift.Field(signal_space, val=data.val),
+             name="data.png", **plotdict)
+    ift.plot(m, name="map.png", **plotdict)
     # Probing the uncertainty
     class Proby(ift.DiagonalProberMixin, ift.Prober):
         pass
-    proby = Proby(signal_space, probe_count=1, ncpu=1)
-    proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
+    proby = Proby(harmonic_space, probe_count=1, ncpu=1)
+    proby(lambda z: wiener_curvature.inverse_times(z))
 
     sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03)
     variance = ift.sqrt(sm(proby.diagonal.weight(-1)))
diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py
index 5884b1451991bf4a8671aeda4833e9a3b6711074..09f946326056ebf2cb3d3ebee96c52b91f8c83bc 100644
--- a/demos/wiener_filter_easy.py
+++ b/demos/wiener_filter_easy.py
@@ -54,7 +54,7 @@ if __name__ == "__main__":
     # Wiener filter
 
     j = R.adjoint_times(N.inverse_times(d))
-    IC = ift.GradientNormController(verbose=True, iteration_limit=500,
+    IC = ift.GradientNormController(name="inverter", iteration_limit=500,
                                     tol_abs_gradnorm=0.1)
     inverter = ift.ConjugateGradient(controller=IC)
     S_inv = fft.adjoint*Sh.inverse*fft
diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py
index 16ee1acb47cada4981a88f52a64a312d007e6eb9..0e3309b748fe63c700f598af1755cab1d45e4270 100644
--- a/demos/wiener_filter_via_curvature.py
+++ b/demos/wiener_filter_via_curvature.py
@@ -69,7 +69,7 @@ if __name__ == "__main__":
 
     j = R.adjoint_times(N.inverse_times(data))
     ctrl = ift.GradientNormController(
-        verbose=True, tol_abs_gradnorm=1e-5/(nu.K*(nu.m**dimensionality)))
+        name="inverter", tol_abs_gradnorm=1e-5/(nu.K*(nu.m**dimensionality)))
     inverter = ift.ConjugateGradient(controller=ctrl)
     wiener_curvature = ift.library.WienerFilterCurvature(
         S=S, N=N, R=R, inverter=inverter)
diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py
index eeeff78153da2272308e16ecf28389aecc954e3f..74f5ae90c082dfd8039077e2c19456e54798806b 100644
--- a/demos/wiener_filter_via_hamiltonian.py
+++ b/demos/wiener_filter_via_hamiltonian.py
@@ -50,9 +50,9 @@ if __name__ == "__main__":
 
     # Choosing the minimization strategy
 
-    ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
+    ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
     inverter = ift.ConjugateGradient(controller=ctrl)
-    controller = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1)
+    controller = ift.GradientNormController(name="min", tol_abs_gradnorm=0.1)
     minimizer = ift.RelaxedNewton(controller=controller)
     m0 = ift.Field.zeros(h_space)
     # Initializing the Wiener Filter energy
diff --git a/nifty4/minimization/gradient_norm_controller.py b/nifty4/minimization/gradient_norm_controller.py
index 2c9c9a4f227b618161487b9e9f73b541332d6858..c06a9988dfd568479178ae8f6fa88e57a469f0a4 100644
--- a/nifty4/minimization/gradient_norm_controller.py
+++ b/nifty4/minimization/gradient_norm_controller.py
@@ -22,15 +22,13 @@ from .. import dobj
 
 class GradientNormController(IterationController):
     def __init__(self, tol_abs_gradnorm=None, tol_rel_gradnorm=None,
-                 convergence_level=1, iteration_limit=None,
-                 name=None, verbose=None):
+                 convergence_level=1, iteration_limit=None, name=None):
         super(GradientNormController, self).__init__()
         self._tol_abs_gradnorm = tol_abs_gradnorm
         self._tol_rel_gradnorm = tol_rel_gradnorm
         self._convergence_level = convergence_level
         self._iteration_limit = iteration_limit
         self._name = name
-        self._verbose = verbose
 
     def start(self, energy):
         self._itcount = -1
@@ -56,10 +54,8 @@ class GradientNormController(IterationController):
             self._ccount = max(0, self._ccount-1)
 
         # report
-        if self._verbose:
-            msg = ""
-            if self._name is not None:
-                msg += self._name+":"
+        if self._name is not None:
+            msg = self._name+":"
             msg += " Iteration #" + str(self._itcount)
             msg += " energy=" + str(energy.value)
             msg += " gradnorm=" + str(energy.gradient_norm)
diff --git a/nifty4/minimization/line_search_strong_wolfe.py b/nifty4/minimization/line_search_strong_wolfe.py
index 1a139892bc802c61f6893babc244edb28fd009e9..7987a6271a09f3e06aae74cc5d515e0cabbc56a8 100644
--- a/nifty4/minimization/line_search_strong_wolfe.py
+++ b/nifty4/minimization/line_search_strong_wolfe.py
@@ -209,8 +209,10 @@ class LineSearchStrongWolfe(LineSearch):
         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.
+        if phi_lo > phi_0 + self.c1*alpha_lo*phiprime_0:
+            raise ValueError("inconsistent data")
+        if phiprime_lo*(alpha_hi-alpha_lo) >= 0.:
+            raise ValueError("inconsistent data")
         for i in range(self.max_zoom_iterations):
             # assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0
             # assert phiprime_lo*(alpha_hi-alpha_lo)<0.
diff --git a/nifty4/operators/harmonic_transform_operator.py b/nifty4/operators/harmonic_transform_operator.py
index c7e368ccb5af9e43d378af94363658f506485f7e..d3d8cf29749f55b39c94ca51ff0afa4b25b84e41 100644
--- a/nifty4/operators/harmonic_transform_operator.py
+++ b/nifty4/operators/harmonic_transform_operator.py
@@ -149,8 +149,9 @@ class HarmonicTransformOperator(LinearOperator):
 
     def _slice_p2h(self, inp):
         rr = self.sjob.alm2map_adjoint(inp)
-        assert len(rr) == ((self.mmax+1)*(self.mmax+2))//2 + \
-                          (self.mmax+1)*(self.lmax-self.mmax)
+        if len(rr) != ((self.mmax+1)*(self.mmax+2))//2 + \
+                      (self.mmax+1)*(self.lmax-self.mmax):
+            raise ValueError("array length mismatch")
         res = np.empty(2*len(rr)-self.lmax-1, dtype=rr[0].real.dtype)
         res[0:self.lmax+1] = rr[0:self.lmax+1].real
         res[self.lmax+1::2] = np.sqrt(2)*rr[self.lmax+1:].real
@@ -159,8 +160,9 @@ class HarmonicTransformOperator(LinearOperator):
 
     def _slice_h2p(self, inp):
         res = np.empty((len(inp)+self.lmax+1)//2, dtype=(inp[0]*1j).dtype)
-        assert len(res) == ((self.mmax+1)*(self.mmax+2))//2 + \
-                           (self.mmax+1)*(self.lmax-self.mmax)
+        if len(res) != ((self.mmax+1)*(self.mmax+2))//2 + \
+                       (self.mmax+1)*(self.lmax-self.mmax):
+            raise ValueError("array length mismatch")
         res[0:self.lmax+1] = inp[0:self.lmax+1]
         res[self.lmax+1:] = np.sqrt(0.5)*(inp[self.lmax+1::2] +
                                           1j*inp[self.lmax+2::2])
diff --git a/nifty4/spaces/power_space.py b/nifty4/spaces/power_space.py
index f5b9d38a05aaa59a5f677f56890fd032264ce3c2..c4ec0b547dbce1325d0efe36d2d68f66c9e6c3ae 100644
--- a/nifty4/spaces/power_space.py
+++ b/nifty4/spaces/power_space.py
@@ -62,7 +62,8 @@ class PowerSpace(Space):
             units of the harmonic partner space.
         """
         nbin = int(nbin)
-        assert nbin >= 3, "nbin must be at least 3"
+        if nbin < 3:
+            raise ValueError("nbin must be at least 3")
         return np.linspace(float(first_bound), float(last_bound), nbin-1)
 
     @staticmethod
@@ -81,7 +82,8 @@ class PowerSpace(Space):
             units of the harmonic partner space.
         """
         nbin = int(nbin)
-        assert nbin >= 3, "nbin must be at least 3"
+        if nbin < 3:
+            raise ValueError("nbin must be at least 3")
         return np.logspace(np.log(float(first_bound)),
                            np.log(float(last_bound)),
                            nbin-1, base=np.e)
@@ -107,7 +109,8 @@ class PowerSpace(Space):
         nbin_max = int((dists[-1]-dists[0])/binsz_min)+2
         if nbin is None:
             nbin = nbin_max
-        assert nbin >= 3, "nbin must be at least 3"
+        if nbin < 3:
+            raise ValueError("nbin must be at least 3")
         if nbin > nbin_max:
             raise ValueError("nbin is too large")
         if logarithmic:
@@ -141,16 +144,19 @@ class PowerSpace(Space):
                 tbb = binbounds
             locdat = np.searchsorted(tbb, dobj.local_data(k_length_array.val))
             temp_pindex = dobj.from_local_data(
-                k_length_array.val.shape, locdat, dobj.distaxis(k_length_array.val))
+                k_length_array.val.shape, locdat,
+                dobj.distaxis(k_length_array.val))
             nbin = len(tbb)+1
             temp_rho = np.bincount(dobj.local_data(temp_pindex).ravel(),
                                    minlength=nbin)
             temp_rho = dobj.np_allreduce_sum(temp_rho)
-            assert not (temp_rho == 0).any(), "empty bins detected"
+            if (temp_rho == 0).any():
+                raise ValueError("empty bins detected")
             # The explicit conversion to float64 is necessary because bincount
             # sometimes returns its result as an integer array, even when
             # floating-point weights are present ...
-            temp_k_lengths = np.bincount(dobj.local_data(temp_pindex).ravel(),
+            temp_k_lengths = np.bincount(
+                dobj.local_data(temp_pindex).ravel(),
                 weights=dobj.local_data(k_length_array.val).ravel(),
                 minlength=nbin).astype(np.float64, copy=False)
             temp_k_lengths = dobj.np_allreduce_sum(temp_k_lengths) / temp_rho
diff --git a/test/test_minimization/test_minimizers.py b/test/test_minimization/test_minimizers.py
index 30b54e33477f6a20541d331aeb2b60f0b61cba64..31460d29270c1da5677d465250a7f331e4d72e85 100644
--- a/test/test_minimization/test_minimizers.py
+++ b/test/test_minimization/test_minimizers.py
@@ -41,7 +41,7 @@ class Test_Minimizers(unittest.TestCase):
         covariance = ift.DiagonalOperator(covariance_diagonal)
         required_result = ift.Field.ones(space, dtype=np.float64)
 
-        IC = ift.GradientNormController(verbose=True,tol_abs_gradnorm=1e-5, iteration_limit=1000)
+        IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000)
         try:
             minimizer = minimizer_class(controller=IC)
             energy = ift.QuadraticEnergy(A=covariance, b=required_result,