diff --git a/ChangeLog.md b/ChangeLog.md
index 158e409142612895353cdab926512beddced38da..c572010ec7da43b24b7c39fa74f28127faa41150 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -1,6 +1,29 @@
 Changes since NIFTy 6
 =====================
 
+CorrelatedFieldMaker interface change
+-------------------------------------
+
+The interface of `ift.CorrelatedFieldMaker.make` and
+`ift.CorrelatedFieldMaker.add_fluctuations` changed; it now expects the mean
+and the standard deviation of their various parameters not as separate
+arguments but as a tuple.
+
+Furthermore, it is now possible to disable the asperity and the flexibility
+together with the asperity in the correlated field model. Note that disabling
+only the flexibility is not possible.
+
+SimpleCorrelatedField
+---------------------
+
+A simplified version of the correlated field model was introduced which does not
+allow for multiple power spectra, the presence of a degree of freedom parameter
+`dofdex`, or `total_N` larger than zero. Except for the above mentioned
+limitations, it is equivalent to `ift.CorrelatedFieldMaker`. Hence, if one
+wants to understand the implementation idea behind the model, it is easier to
+grasp from reading `ift.SimpleCorrelatedField` than from going through
+`ift.CorrelatedFieldMaker`.
+
 Change in external dependencies
 -------------------------------
 
@@ -20,7 +43,7 @@ The implementation tests for nonlinear operators are now available in
 MetricGaussianKL interface
 --------------------------
 
-Users do not instanciate `MetricGaussianKL` by its constructor anymore. Rather
+Users do not instantiate `MetricGaussianKL` by its constructor anymore. Rather
 `MetricGaussianKL.make()` shall be used.
 
 
@@ -91,10 +114,10 @@ New approach for sampling complex numbers
 When calling draw_sample_with_dtype with a complex dtype,
 the variance is now used for the imaginary part and real part separately.
 This is done in order to be consistent with the Hamiltonian.
-Note that by this, 
+Note that by this,
 ```
 np.std(ift.from_random(domain, 'normal', dtype=np.complex128).val)
-```` 
+````
 does not give 1, but sqrt(2) as a result.
 
 
diff --git a/demos/getting_started_4_CorrelatedFields.ipynb b/demos/getting_started_4_CorrelatedFields.ipynb
index 225ad195cb073c0feb68c67c81a4bf2df84355ce..072f401a1872f01b27a4d691be5ada1d41effa52 100644
--- a/demos/getting_started_4_CorrelatedFields.ipynb
+++ b/demos/getting_started_4_CorrelatedFields.ipynb
@@ -135,7 +135,7 @@
     "## Before the Action: The Moment-Matched Log-Normal Distribution\n",
     "Many properties of the correlated field are modelled as being lognormally distributed.\n",
     "\n",
-    "The distribution models are parametrized via their means `_mean` and standard-deviations `_stddev`.\n",
+    "The distribution models are parametrized via their means and standard-deviations (first and second position in tuple).\n",
     "\n",
     "To get a feeling of how the ratio of the `mean` and `stddev` parameters influences the distribution shape,\n",
     "here are a few example histograms: (observe the x-axis!)"
@@ -189,21 +189,16 @@
     "# Neutral model parameters yielding a quasi-constant field\n",
     "cf_make_pars = {\n",
     "    'offset_mean': 0.,\n",
-    "    'offset_std_mean': 1e-3,\n",
-    "    'offset_std_std': 1e-16,\n",
+    "    'offset_std': (1e-3, 1e-16),\n",
     "    'prefix': ''\n",
     "}\n",
     "\n",
     "cf_x_fluct_pars = {\n",
     "    'target_subdomain': x_space,\n",
-    "    'fluctuations_mean': 1e-3,\n",
-    "    'fluctuations_stddev': 1e-16,\n",
-    "    'flexibility_mean': 1e-3,\n",
-    "    'flexibility_stddev': 1e-16,\n",
-    "    'asperity_mean': 1e-3,\n",
-    "    'asperity_stddev': 1e-16,\n",
-    "    'loglogavgslope_mean': 0.,\n",
-    "    'loglogavgslope_stddev': 1e-16\n",
+    "    'fluctuations': (1e-3, 1e-16),\n",
+    "    'flexibility': (1e-3, 1e-16),\n",
+    "    'asperity': (1e-3, 1e-16),\n",
+    "    'loglogavgslope': (0., 1e-16)\n",
     "}\n",
     "\n",
     "init_model(cf_make_pars, cf_x_fluct_pars)"
@@ -225,18 +220,18 @@
    "source": [
     "# Parameter Showcases\n",
     "\n",
-    "## The `fluctuations_` parameters of `add_fluctuations()`\n",
+    "## The `fluctuations` parameters of `add_fluctuations()`\n",
     "\n",
     "determine the **amplitude of variations along the field dimension**\n",
     "for which `add_fluctuations` is called.\n",
     "\n",
-    "`fluctuations_mean` set the _average_ amplitude of the fields fluctuations along the given dimension,\\\n",
-    "`fluctuations_stddev` sets the width and shape of the amplitude distribution.\n",
+    "`fluctuations[0]` set the _average_ amplitude of the fields fluctuations along the given dimension,\\\n",
+    "`fluctuations[1]` sets the width and shape of the amplitude distribution.\n",
     "\n",
     "The amplitude is modelled as being log-normally distributed,\n",
     "see `The Moment-Matched Log-Normal Distribution` above for details.\n",
     "\n",
-    "#### `fluctuations_mean`:"
+    "#### `fluctuations` mean:"
    ]
   },
   {
@@ -245,14 +240,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "vary_parameter('fluctuations_mean', [0.05, 0.5, 2.], samples_vary_in='xi')"
+    "vary_parameter('fluctuations', [(0.05, 1e-16), (0.5, 1e-16), (2., 1e-16)], samples_vary_in='xi')"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#### `fluctuations_stddev`:"
+    "#### `fluctuations` std:"
    ]
   },
   {
@@ -261,21 +256,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cf_x_fluct_pars['fluctuations_mean'] = 1.0\n",
-    "vary_parameter('fluctuations_stddev', [0.01, 0.1, 1.0], samples_vary_in='fluctuations')"
+    "vary_parameter('fluctuations', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='fluctuations')\n",
+    "cf_x_fluct_pars['fluctuations'] = (1., 1e-16)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## The `loglogavgslope_` parameters of `add_fluctuations()`\n",
+    "## The `loglogavgslope` parameters of `add_fluctuations()`\n",
     "\n",
     "determine **the slope of the loglog-linear (power law) component of the power spectrum**.\n",
     "\n",
     "The slope is modelled to be normally distributed.\n",
     "\n",
-    "#### `loglogavgslope_mean`:"
+    "#### `loglogavgslope` mean:"
    ]
   },
   {
@@ -284,14 +279,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "vary_parameter('loglogavgslope_mean', [-3., -1., 1.], samples_vary_in='xi')"
+    "vary_parameter('loglogavgslope', [(-3., 1e-16), (-1., 1e-16), (1., 1e-16)], samples_vary_in='xi')"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#### `loglogavgslope_stddev`:"
+    "#### `loglogavgslope` std:"
    ]
   },
   {
@@ -300,25 +295,25 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cf_x_fluct_pars['loglogavgslope_mean'] = -1.0\n",
-    "vary_parameter('loglogavgslope_stddev', [0.01, 0.1, 1.0], samples_vary_in='loglogavgslope')"
+    "vary_parameter('loglogavgslope', [(-1., 0.01), (-1., 0.1), (-1., 1.0)], samples_vary_in='loglogavgslope')\n",
+    "cf_x_fluct_pars['loglogavgslope'] = (-1., 1e-16)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## The `flexibility_` parameters of `add_fluctuations()`\n",
+    "## The `flexibility` parameters of `add_fluctuations()`\n",
     "\n",
     "determine **the amplitude of the integrated Wiener process component of the power spectrum**\n",
     "(how strong the power spectrum varies besides the power-law).\n",
     "\n",
-    "`flexibility_mean` sets the _average_ amplitude of the i.g.p. component,\n",
-    "`flexibility_stddev` sets how much the amplitude can vary.\\\n",
+    "`flexibility[0]` sets the _average_ amplitude of the i.g.p. component,\\\n",
+    "`flexibility[1]` sets how much the amplitude can vary.\\\n",
     "These two parameters feed into a moment-matched log-normal distribution model,\n",
     "see above for a demo of its behavior.\n",
     "\n",
-    "#### `flexibility_mean`:"
+    "#### `flexibility` mean:"
    ]
   },
   {
@@ -327,14 +322,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "vary_parameter('flexibility_mean', [0.1, 1.0, 3.0], samples_vary_in='spectrum')"
+    "vary_parameter('flexibility', [(0.1, 1e-16), (1.0, 1e-16), (3.0, 1e-16)], samples_vary_in='spectrum')"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#### `flexibility_stddev`:"
+    "#### `flexibility` std:"
    ]
   },
   {
@@ -343,23 +338,23 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cf_x_fluct_pars['flexibility_mean'] = 2.0\n",
-    "vary_parameter('flexibility_stddev', [0.01, 0.1, 1.0], samples_vary_in='flexibility')"
+    "vary_parameter('flexibility', [(2., 0.01), (2., 0.1), (2., 1.)], samples_vary_in='flexibility')\n",
+    "cf_x_fluct_pars['flexibility'] = (2., 1e-16)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## The `asperity_` parameters of `add_fluctuations()`\n",
+    "## The `asperity` parameters of `add_fluctuations()`\n",
     "\n",
-    "`asperity_` determines **how rough the integrated Wiener process component of the power spectrum is**.\n",
+    "`asperity` determines **how rough the integrated Wiener process component of the power spectrum is**.\n",
     "\n",
-    "`asperity_mean` sets the average roughness, `asperity_stddev` sets how much the roughness can vary.\\\n",
+    "`asperity[0]` sets the average roughness, `asperity[1]` sets how much the roughness can vary.\\\n",
     "These two parameters feed into a moment-matched log-normal distribution model,\n",
     "see above for a demo of its behavior.\n",
     "\n",
-    "#### `asperity_mean`:"
+    "#### `asperity` mean:"
    ]
   },
   {
@@ -368,14 +363,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "vary_parameter('asperity_mean', [0.001, 1.0, 5.0], samples_vary_in='spectrum')"
+    "vary_parameter('asperity', [(0.001, 1e-16), (1.0, 1e-16), (5., 1e-16)], samples_vary_in='spectrum')"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#### `asperity_stddev`:"
+    "#### `asperity` std:"
    ]
   },
   {
@@ -384,8 +379,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cf_x_fluct_pars['asperity_mean'] = 1.0\n",
-    "vary_parameter('asperity_stddev', [0.01, 0.1, 1.0], samples_vary_in='asperity')"
+    "vary_parameter('asperity', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='asperity')\n",
+    "cf_x_fluct_pars['asperity'] = (1., 1e-16)"
    ]
   },
   {
@@ -406,10 +401,10 @@
    "outputs": [],
    "source": [
     "# Reset model to neutral\n",
-    "cf_x_fluct_pars['fluctuations_mean'] = 1e-3\n",
-    "cf_x_fluct_pars['flexibility_mean'] = 1e-3\n",
-    "cf_x_fluct_pars['asperity_mean'] = 1e-3\n",
-    "cf_x_fluct_pars['loglogavgslope_mean'] = 1e-3"
+    "cf_x_fluct_pars['fluctuations'] = (1e-3, 1e-16)\n",
+    "cf_x_fluct_pars['flexibility'] = (1e-3, 1e-16)\n",
+    "cf_x_fluct_pars['asperity'] = (1e-3, 1e-16)\n",
+    "cf_x_fluct_pars['loglogavgslope'] = (1e-3, 1e-16)"
    ]
   },
   {
@@ -425,15 +420,15 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## The `offset_std_` parameters of `CorrelatedFieldMaker.make()`\n",
+    "## The `offset_std` parameters of `CorrelatedFieldMaker.make()`\n",
     "\n",
     "Variation of the global offset of the field are modelled as being log-normally distributed.\n",
     "See `The Moment-Matched Log-Normal Distribution` above for details.\n",
     "\n",
-    "The `offset_std_mean` parameter sets how much NIFTy will vary the offset *on average*.\\\n",
-    "The `offset_std_std` parameters defines the with and shape of the offset variation distribution.\n",
+    "The `offset_std[0]` parameter sets how much NIFTy will vary the offset *on average*.\\\n",
+    "The `offset_std[1]` parameters defines the with and shape of the offset variation distribution.\n",
     "\n",
-    "#### `offset_std_mean`:"
+    "#### `offset_std` mean:"
    ]
   },
   {
@@ -442,14 +437,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "vary_parameter('offset_std_mean', [1e-16, 0.5, 2.], samples_vary_in='xi')"
+    "vary_parameter('offset_std', [(1e-16, 1e-16), (0.5, 1e-16), (2., 1e-16)], samples_vary_in='xi')"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "#### `offset_std_std`:"
+    "#### `offset_std` std:"
    ]
   },
   {
@@ -458,8 +453,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cf_make_pars['offset_std_mean'] = 1.0\n",
-    "vary_parameter('offset_std_std', [0.01, 0.1, 1.0], samples_vary_in='zeromode')"
+    "vary_parameter('offset_std', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='zeromode')"
    ]
   }
  ],
@@ -479,7 +473,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.2"
+   "version": "3.8.3"
   }
  },
  "nbformat": 4,
diff --git a/demos/getting_started_5_mf.py b/demos/getting_started_5_mf.py
index 84d9b77cd4e6b6347c00ca67a0ea5e1d22c688ee..d4225094c281c80ac7b34e2ea1b4e57c3f452e6f 100644
--- a/demos/getting_started_5_mf.py
+++ b/demos/getting_started_5_mf.py
@@ -73,10 +73,10 @@ def main():
     sp2 = ift.RGSpace(npix2)
 
     # Set up signal model
-    cfmaker = ift.CorrelatedFieldMaker.make(0., 1e-2, 1e-6, '')
-    cfmaker.add_fluctuations(sp1, 0.1, 1e-2, 1, .1, .01, .5, -2, 1., 'amp1')
-    cfmaker.add_fluctuations(sp2, 0.1, 1e-2, 1, .1, .01, .5,
-                             -1.5, .5, 'amp2')
+    cfmaker = ift.CorrelatedFieldMaker.make(0., (1e-2, 1e-6), '')
+    cfmaker.add_fluctuations(sp1, (0.1, 1e-2), (1, .1), (.01, .5), (-2, 1.), 'amp1')
+    cfmaker.add_fluctuations(sp2, (0.1, 1e-2), (1, .1), (.01, .5),
+                             (-1.5, .5), 'amp2')
     correlated_field = cfmaker.finalize()
 
     A1 = cfmaker.normalized_amplitudes[0]
diff --git a/src/__init__.py b/src/__init__.py
index 1f255f8487e36021a99988e1f5362a7128b1a6b4..d5e0b9b9ec70e0b4714ad2ad85b930e48c3b677b 100644
--- a/src/__init__.py
+++ b/src/__init__.py
@@ -92,7 +92,7 @@ from .library.correlated_fields_simple import SimpleCorrelatedField
 
 from . import extra
 
-from .utilities import memo, frozendict
+from .utilities import memo, frozendict, myassert
 
 from .logger import logger
 
diff --git a/src/extra.py b/src/extra.py
index 69f1ecae0c0ad45d21398e0584dc5536b17e045d..488d3b6240839766b393493426b60359badb5b3f 100644
--- a/src/extra.py
+++ b/src/extra.py
@@ -18,7 +18,6 @@
 from itertools import combinations
 
 import numpy as np
-from numpy.testing import assert_
 
 from .domain_tuple import DomainTuple
 from .field import Field
@@ -29,6 +28,7 @@ from .operators.energy_operators import EnergyOperator
 from .operators.linear_operator import LinearOperator
 from .operators.operator import Operator
 from .sugar import from_random
+from .utilities import myassert
 
 __all__ = ["check_linear_operator", "check_operator",
            "assert_allclose"]
@@ -137,20 +137,6 @@ def assert_equal(f1, f2):
         assert_equal(val, f2[key])
 
 
-def _nozero(fld):
-    if isinstance(fld, Field):
-        return np.testing.assert_((fld != 0).s_all())
-    for val in fld.values():
-        _nozero(val)
-
-
-def _allzero(fld):
-    if isinstance(fld, Field):
-        return np.testing.assert_((fld == 0.).s_all())
-    for val in fld.values():
-        _allzero(val)
-
-
 def _adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol,
                             only_r_linear):
     needed_cap = op.TIMES | op.ADJOINT_TIMES
@@ -206,32 +192,32 @@ def _domain_check_linear(op, domain_dtype=None, inp=None):
         inp = from_random(op.domain, "normal", dtype=domain_dtype)
     elif inp is None:
         raise ValueError('Need to specify either dtype or inp')
-    assert_(inp.domain is op.domain)
-    assert_(op(inp).domain is op.target)
+    myassert(inp.domain is op.domain)
+    myassert(op(inp).domain is op.target)
 
 
 def _domain_check_nonlinear(op, loc):
     _domain_check(op)
-    assert_(isinstance(loc, (Field, MultiField)))
-    assert_(loc.domain is op.domain)
+    myassert(isinstance(loc, (Field, MultiField)))
+    myassert(loc.domain is op.domain)
     for wm in [False, True]:
         lin = Linearization.make_var(loc, wm)
         reslin = op(lin)
-        assert_(lin.domain is op.domain)
-        assert_(lin.target is op.domain)
-        assert_(lin.val.domain is lin.domain)
-        assert_(reslin.domain is op.domain)
-        assert_(reslin.target is op.target)
-        assert_(reslin.val.domain is reslin.target)
-        assert_(reslin.target is op.target)
-        assert_(reslin.jac.domain is reslin.domain)
-        assert_(reslin.jac.target is reslin.target)
-        assert_(lin.want_metric == reslin.want_metric)
+        myassert(lin.domain is op.domain)
+        myassert(lin.target is op.domain)
+        myassert(lin.val.domain is lin.domain)
+        myassert(reslin.domain is op.domain)
+        myassert(reslin.target is op.target)
+        myassert(reslin.val.domain is reslin.target)
+        myassert(reslin.target is op.target)
+        myassert(reslin.jac.domain is reslin.domain)
+        myassert(reslin.jac.target is reslin.target)
+        myassert(lin.want_metric == reslin.want_metric)
         _domain_check_linear(reslin.jac, inp=loc)
         _domain_check_linear(reslin.jac.adjoint, inp=reslin.jac(loc))
         if reslin.metric is not None:
-            assert_(reslin.metric.domain is reslin.metric.target)
-            assert_(reslin.metric.domain is op.domain)
+            myassert(reslin.metric.domain is reslin.metric.target)
+            myassert(reslin.metric.domain is op.domain)
 
 
 def _domain_check(op):
diff --git a/src/library/correlated_fields.py b/src/library/correlated_fields.py
index 886484066b9a574bb2b69c1a51f464c3012c88c0..5a7b1c8446f7e8899df65823a36e1544e23c5ec6 100644
--- a/src/library/correlated_fields.py
+++ b/src/library/correlated_fields.py
@@ -209,13 +209,13 @@ class _Amplitude(Operator):
                  loglogavgslope, azm, totvol, key, dofdex):
         """
         fluctuations > 0
-        flexibility > 0
-        asperity > 0
+        flexibility > 0 or None
+        asperity > 0 or None
         loglogavgslope probably negative
         """
         assert isinstance(fluctuations, Operator)
-        assert isinstance(flexibility, Operator)
-        assert isinstance(asperity, Operator)
+        assert isinstance(flexibility, Operator) or flexibility is None
+        assert isinstance(asperity, Operator) or asperity is None
         assert isinstance(loglogavgslope, Operator)
 
         if len(dofdex) > 0:
@@ -240,43 +240,53 @@ class _Amplitude(Operator):
         ps_expander = ContractionOperator(twolog.target, spaces=space).adjoint
 
         # Prepare constant fields
-        foo = np.zeros(shp)
-        foo[0] = foo[1] = np.sqrt(_log_vol(target[space]))
-        vflex = DiagonalOperator(makeField(dom[space], foo), dom, space)
+        vflex = np.zeros(shp)
+        vflex[0] = vflex[1] = np.sqrt(_log_vol(target[space]))
+        vflex = DiagonalOperator(makeField(dom[space], vflex), dom, space)
 
-        foo = np.zeros(shp, dtype=np.float64)
-        foo[0] += 1
-        vasp = DiagonalOperator(makeField(dom[space], foo), dom, space)
+        vasp = np.zeros(shp, dtype=np.float64)
+        vasp[0] += 1
+        vasp = DiagonalOperator(makeField(dom[space], vasp), dom, space)
 
-        foo = np.ones(shp)
-        foo[0] = _log_vol(target[space])**2/12.
-        shift = DiagonalOperator(makeField(dom[space], foo), dom, space)
+        shift = np.ones(shp)
+        shift[0] = _log_vol(target[space])**2 / 12.
+        shift = DiagonalOperator(makeField(dom[space], shift), dom, space)
+        shift = shift(full(shift.domain, 1))
 
         vslope = DiagonalOperator(
             makeField(target[space], _relative_log_k_lengths(target[space])),
             target, space)
 
-        foo, bar = [np.zeros(target[space].shape) for _ in range(2)]
-        bar[1:] = foo[0] = totvol
+        vol0, vol1 = [np.zeros(target[space].shape) for _ in range(2)]
+        vol1[1:] = vol0[0] = totvol
         vol0, vol1 = [
             DiagonalOperator(makeField(target[space], aa), target, space)
-            for aa in (foo, bar)
+            for aa in (vol0, vol1)
         ]
-
-        # Prepare fields for Adder
-        shift, vol0 = [op(full(op.domain, 1)) for op in (shift, vol0)]
+        vol0 = vol0(full(vol0.domain, 1))
         # End prepare constant fields
 
         slope = vslope @ ps_expander @ loglogavgslope
-        sig_flex = vflex @ expander @ flexibility
-        sig_asp = vasp @ expander @ asperity
+        sig_flex = vflex @ expander @ flexibility if flexibility is not None else None
+        sig_asp = vasp @ expander @ asperity if asperity is not None else None
         sig_fluc = vol1 @ ps_expander @ fluctuations
         sig_fluc = vol1 @ ps_expander @ fluctuations
 
-        xi = ducktape(dom, None, key)
-        sigma = sig_flex*(Adder(shift) @ sig_asp).ptw("sqrt")
-        smooth = _SlopeRemover(target, space) @ twolog @ (sigma*xi)
-        op = _Normalization(target, space) @ (slope + smooth)
+        if sig_asp is None and sig_flex is None:
+            op = _Normalization(target, space) @ slope
+        elif sig_asp is None:
+            xi = ducktape(dom, None, key)
+            sigma = DiagonalOperator(shift.ptw("sqrt"), dom, space) @ sig_flex
+            smooth = _SlopeRemover(target, space) @ twolog @ (sigma * xi)
+            op = _Normalization(target, space) @ (slope + smooth)
+        elif sig_flex is None:
+            raise ValueError("flexibility may not be disabled on its own")
+        else:
+            xi = ducktape(dom, None, key)
+            sigma = sig_flex * (Adder(shift) @ sig_asp).ptw("sqrt")
+            smooth = _SlopeRemover(target, space) @ twolog @ (sigma * xi)
+            op = _Normalization(target, space) @ (slope + smooth)
+
         if N_copies > 0:
             op = Distributor @ op
             sig_fluc = Distributor @ sig_fluc
@@ -301,7 +311,7 @@ class _Amplitude(Operator):
 
 
 class CorrelatedFieldMaker:
-    """Constrution helper for hirarchical correlated field models.
+    """Construction helper for hierarchical correlated field models.
 
     The correlated field models are parametrized by creating
     power spectrum operators ("amplitudes") via calls to
@@ -347,18 +357,16 @@ class CorrelatedFieldMaker:
         self._total_N = total_N
 
     @staticmethod
-    def make(offset_mean, offset_std_mean, offset_std_std, prefix, total_N=0,
-             dofdex=None):
+    def make(offset_mean, offset_std, prefix, total_N=0, dofdex=None):
         """Returns a CorrelatedFieldMaker object.
 
         Parameters
         ----------
         offset_mean : float
             Mean offset from zero of the correlated field to be made.
-        offset_std_mean : float
-            Mean standard deviation of the offset.
-        offset_std_std : float
-            Standard deviation of the stddev of the offset.
+        offset_std : tuple of float
+            Mean standard deviation and standard deviation of the standard
+            deviation of the offset. No, this is not a word duplication.
         prefix : string
             Prefix to the names of the domains of the cf operator to be made.
             This determines the names of the operator domain.
@@ -371,7 +379,7 @@ class CorrelatedFieldMaker:
             An integer array specifying the zero mode models used if
             total_N > 1. It needs to have length of total_N. If total_N=3 and
             dofdex=[0,0,1], that means that two models for the zero mode are
-            instanciated; the first one is used for the first and second
+            instantiated; the first one is used for the first and second
             field model and the second is used for the third field model.
             *If not specified*, use the same zero mode model for all
             constructed field models.
@@ -381,22 +389,19 @@ class CorrelatedFieldMaker:
         elif len(dofdex) != total_N:
             raise ValueError("length of dofdex needs to match total_N")
         N = max(dofdex) + 1 if total_N > 0 else 0
-        zm = LognormalTransform(offset_std_mean, offset_std_std,
-                                prefix + 'zeromode', N)
+        if len(offset_std) != 2:
+            raise TypeError
+        zm = LognormalTransform(*offset_std, prefix + 'zeromode', N)
         if total_N > 0:
             zm = _Distributor(dofdex, zm.target, UnstructuredDomain(total_N)) @ zm
         return CorrelatedFieldMaker(offset_mean, zm, prefix, total_N)
 
     def add_fluctuations(self,
                          target_subdomain,
-                         fluctuations_mean,
-                         fluctuations_stddev,
-                         flexibility_mean,
-                         flexibility_stddev,
-                         asperity_mean,
-                         asperity_stddev,
-                         loglogavgslope_mean,
-                         loglogavgslope_stddev,
+                         fluctuations,
+                         flexibility,
+                         asperity,
+                         loglogavgslope,
                          prefix='',
                          index=None,
                          dofdex=None,
@@ -411,8 +416,8 @@ class CorrelatedFieldMaker:
         used on the target field subdomain `target_subdomain`.
         It is assembled as the sum of a power law component
         (linear slope in log-log power-frequency-space),
-        a smooth varying component (integrated wiener process) and
-        a ragged componenent (unintegrated wiener process).
+        a smooth varying component (integrated Wiener process) and
+        a ragged component (un-integrated Wiener process).
 
         Multiple calls to `add_fluctuations` are possible, in which case
         the constructed field will have the outer product of the individual
@@ -424,14 +429,14 @@ class CorrelatedFieldMaker:
                            :class:`~nifty7.domain_tuple.DomainTuple`
             Target subdomain on which the correlation structure defined
             in this call should hold.
-        fluctuations_{mean,stddev} : float
+        fluctuations : tuple of float
             Total spectral energy -> Amplitude of the fluctuations
-        flexibility_{mean,stddev} : float
+        flexibility : tuple of float or None
             Amplitude of the non-power-law power spectrum component
-        asperity_{mean,stddev} : float
+        asperity : tuple of float or None
             Roughness of the non-power-law power spectrum component
-            Used to accomodate single frequency peaks
-        loglogavgslope_{mean,stddev} : float
+            Used to accommodate single frequency peaks
+        loglogavgslope : tuple of float
             Power law component exponent
         prefix : string
             prefix of the power spectrum parameter domain names
@@ -442,7 +447,7 @@ class CorrelatedFieldMaker:
             An integer array specifying the power spectrum models used if
             total_N > 1. It needs to have length of total_N. If total_N=3 and
             dofdex=[0,0,1], that means that two power spectrum models are
-            instanciated; the first one is used for the first and second
+            instantiated; the first one is used for the first and second
             field model and the second one is used for the third field model.
             *If not given*, use the same power spectrum model for all
             constructed field models.
@@ -467,21 +472,37 @@ class CorrelatedFieldMaker:
         else:
             N = 0
             target_subdomain = makeDomain(target_subdomain)
-        prefix = str(prefix)
-        # assert isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace)
-
-        fluct = LognormalTransform(fluctuations_mean, fluctuations_stddev,
-                                   self._prefix + prefix + 'fluctuations', N)
-        flex = LognormalTransform(flexibility_mean, flexibility_stddev,
-                                  self._prefix + prefix + 'flexibility', N)
-        asp = LognormalTransform(asperity_mean, asperity_stddev,
-                                 self._prefix + prefix + 'asperity', N)
-        avgsl = NormalTransform(loglogavgslope_mean, loglogavgslope_stddev,
-                                self._prefix + prefix + 'loglogavgslope', N)
+        # assert isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace))
+
+        for arg in [fluctuations, loglogavgslope]:
+            if len(arg) != 2:
+                raise TypeError
+        for kw, arg in [("flexibility", flexibility), ("asperity", asperity)]:
+            if arg is None:
+                continue
+            if len(arg) != 2:
+                raise TypeError
+            if len(arg) == 2 and (arg[0] <= 0. or arg[1] <= 0.):
+                ve = "{0} must be strictly positive (or None)"
+                raise ValueError(ve.format(kw))
+        if flexibility is None and asperity is not None:
+            raise ValueError("flexibility may not be disabled on its own")
+
+        pre = self._prefix + str(prefix)
+        fluct = LognormalTransform(*fluctuations, pre + 'fluctuations', N)
+        if flexibility is not None:
+            flex = LognormalTransform(*flexibility, pre + 'flexibility', N)
+        else:
+            flex = None
+        if asperity is not None:
+            asp = LognormalTransform(*asperity, pre + 'asperity', N)
+        else:
+            asp = None
+        avgsl = NormalTransform(*loglogavgslope, pre + 'loglogavgslope', N)
 
         amp = _Amplitude(PowerSpace(harmonic_partner), fluct, flex, asp, avgsl,
                          self._azm, target_subdomain[-1].total_volume,
-                         self._prefix + prefix + 'spectrum', dofdex)
+                         pre + 'spectrum', dofdex)
 
         if index is not None:
             self._a.insert(index, amp)
diff --git a/src/library/correlated_fields_simple.py b/src/library/correlated_fields_simple.py
index 3209b8c8e19969059ffec3d124cb765743b67c36..6b96d2259ea489db58ffb4a3f70b200e74f8dfbc 100644
--- a/src/library/correlated_fields_simple.py
+++ b/src/library/correlated_fields_simple.py
@@ -78,7 +78,7 @@ class SimpleCorrelatedField(Operator):
             xi = ducktape(dom, None, prefix + 'spectrum')
 
             shift = np.ones(dom.shape)
-            shift[0] = _log_vol(pspace)**2/12.
+            shift[0] = _log_vol(pspace)**2 / 12.
             shift = makeField(dom, shift)
             if asperity is None:
                 asp = makeOp(shift.ptw("sqrt")) @ (xi*sig_flex)
diff --git a/src/utilities.py b/src/utilities.py
index 67f4948aeb14e3e07539e2ea08ee1b17b0f03134..7f843fce03f9655ccf45dfddd2fc90a4a272497c 100644
--- a/src/utilities.py
+++ b/src/utilities.py
@@ -405,3 +405,10 @@ def lognormal_moments(mean, sigma, N=0):
     logsigma = np.sqrt(np.log1p((sigma / mean)**2))
     logmean = np.log(mean) - logsigma**2 / 2
     return logmean, logsigma
+
+
+def myassert(val):
+    """Safe alternative to python's assert statement which is active even if
+    `__debug__` is False."""
+    if not val:
+        raise AssertionError
diff --git a/test/test_kl.py b/test/test_kl.py
index 35d426c736df06bfa2dc38849a9c0907f628d41f..887942709f47d003da051837a66e180bfb6b65d5 100644
--- a/test/test_kl.py
+++ b/test/test_kl.py
@@ -17,9 +17,10 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_allclose, assert_raises
+from numpy.testing import assert_allclose, assert_raises
 
 import nifty7 as ift
+from nifty7 import myassert
 
 from .common import setup_function, teardown_function
 
@@ -55,13 +56,13 @@ def test_kl(constants, point_estimates, mirror_samples, mf):
             ift.MetricGaussianKL.make(**args)
         return
     kl = ift.MetricGaussianKL.make(**args)
-    assert_(len(ic.history) > 0)
-    assert_(len(ic.history) == len(ic.history.time_stamps))
-    assert_(len(ic.history) == len(ic.history.energy_values))
+    myassert(len(ic.history) > 0)
+    myassert(len(ic.history) == len(ic.history.time_stamps))
+    myassert(len(ic.history) == len(ic.history.energy_values))
     ic.history.reset()
-    assert_(len(ic.history) == 0)
-    assert_(len(ic.history) == len(ic.history.time_stamps))
-    assert_(len(ic.history) == len(ic.history.energy_values))
+    myassert(len(ic.history) == 0)
+    myassert(len(ic.history) == len(ic.history.time_stamps))
+    myassert(len(ic.history) == len(ic.history.energy_values))
 
     locsamp = kl._local_samples
     if isinstance(mean0, ift.MultiField):
@@ -74,7 +75,7 @@ def test_kl(constants, point_estimates, mirror_samples, mf):
 
     # Test number of samples
     expected_nsamps = 2*nsamps if mirror_samples else nsamps
-    assert_(len(tuple(kl.samples)) == expected_nsamps)
+    myassert(len(tuple(kl.samples)) == expected_nsamps)
 
     # Test value
     assert_allclose(kl.value, klpure.value)
diff --git a/test/test_linearization.py b/test/test_linearization.py
index 5a3cf0fe7d69f3372739600f36625fd8a8a4e0a7..1bc564ddfcf0b9744ee35f614bfeb4dac13d1079 100644
--- a/test/test_linearization.py
+++ b/test/test_linearization.py
@@ -17,7 +17,7 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_allclose
+from numpy.testing import assert_allclose
 
 import nifty7 as ift
 from .common import setup_function, teardown_function
@@ -44,7 +44,7 @@ def test_special_gradients():
 
     assert_allclose(
         _lin2grad(ift.Linearization.make_var(0*f).ptw("sinc")), np.zeros(s.shape))
-    assert_(np.isnan(_lin2grad(ift.Linearization.make_var(0*f).ptw("abs"))))
+    ift.myassert(np.isnan(_lin2grad(ift.Linearization.make_var(0*f).ptw("abs"))))
     assert_allclose(
         _lin2grad(ift.Linearization.make_var(0*f + 10).ptw("abs")),
         np.ones(s.shape))
diff --git a/test/test_mpi/test_kl.py b/test/test_mpi/test_kl.py
index 557acd4a248e65dbe453cc282822ea698a1f51d2..44c1da13b3acea87e4dac24fe650437388885a85 100644
--- a/test/test_mpi/test_kl.py
+++ b/test/test_mpi/test_kl.py
@@ -18,7 +18,7 @@
 import numpy as np
 import pytest
 from mpi4py import MPI
-from numpy.testing import assert_, assert_equal, assert_raises
+from numpy.testing import assert_equal, assert_raises
 
 import nifty7 as ift
 
@@ -84,8 +84,8 @@ def test_kl(constants, point_estimates, mirror_samples, mode, mf):
 
     # Test number of samples
     expected_nsamps = 2*nsamps if mirror_samples else nsamps
-    assert_(len(tuple(kl0.samples)) == expected_nsamps)
-    assert_(len(tuple(kl1.samples)) == expected_nsamps)
+    ift.myassert(len(tuple(kl0.samples)) == expected_nsamps)
+    ift.myassert(len(tuple(kl1.samples)) == expected_nsamps)
 
     # Test value
     assert_equal(kl0.value, kl1.value)
diff --git a/test/test_operator_tree_optimiser.py b/test/test_operator_tree_optimiser.py
index 7ad35e094d148b06f46321a843cd43b344e2f201..beb34776ecfe184d814510beab62bdd574df53f4 100644
--- a/test/test_operator_tree_optimiser.py
+++ b/test/test_operator_tree_optimiser.py
@@ -16,14 +16,16 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-from numpy.testing import assert_, assert_allclose
-import numpy as np
 from copy import deepcopy
+
+import numpy as np
+from numpy.testing import assert_allclose
+
 import nifty7 as ift
 
 
 class CountingOp(ift.Operator):
-    #FIXME: Not a LinearOperator since ChainOps not supported yet
+    # FIXME: Not a LinearOperator since ChainOps not supported yet
     def __init__(self, domain):
         self._domain = self._target = ift.sugar.makeDomain(domain)
         self._count = 0
@@ -56,6 +58,6 @@ def test_operator_tree_optimiser():
     op_orig = deepcopy(op)
     op = ift.operator_tree_optimiser._optimise_operator(op)
     assert_allclose(op(fld).val, op_orig(fld).val, rtol=np.finfo(np.float64).eps)
-    assert_(1 == ( (cop4.count-1) * cop3.count * cop2.count * cop1.count))
+    ift.myassert(1 == ((cop4.count-1) * cop3.count * cop2.count * cop1.count))
     # test testing
     ift.optimise_operator(op_orig)
diff --git a/test/test_operators/test_correlated_fields.py b/test/test_operators/test_correlated_fields.py
index 3eb26b72198f2ddd3f863b2b94d1fdfd6e7de240..5f486469134cdc7978119e17a7279440732ffa1c 100644
--- a/test/test_operators/test_correlated_fields.py
+++ b/test/test_operators/test_correlated_fields.py
@@ -17,7 +17,7 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_allclose
+from numpy.testing import assert_allclose
 
 import nifty7 as ift
 
@@ -71,11 +71,11 @@ def testAmplitudesInvariants(sspace, N):
 
     astds = 0.2, 1.2
     offset_std_mean = 1.3
-    fa = ift.CorrelatedFieldMaker.make(1.2, offset_std_mean, 1e-2, '', N,
+    fa = ift.CorrelatedFieldMaker.make(1.2, (offset_std_mean, 1e-2), '', N,
                                        dofdex1)
-    fa.add_fluctuations(sspace, astds[0], 1e-2, 1.1, 2., 2.1, .5, -2, 1.,
+    fa.add_fluctuations(sspace, (astds[0], 1e-2), (1.1, 2.), (2.1, .5), (-2, 1.),
                         'spatial', dofdex=dofdex2)
-    fa.add_fluctuations(fsspace, astds[1], 1e-2, 3.1, 1., .5, .1, -4, 1.,
+    fa.add_fluctuations(fsspace, (astds[1], 1e-2), (3.1, 1.), (.5, .1), (-4, 1.),
                         'freq', dofdex=dofdex3)
     op = fa.finalize()
 
@@ -103,19 +103,19 @@ def testAmplitudesInvariants(sspace, N):
     assert_allclose(slice_fluct_std0, sl_fluct_space, rtol=0.5)
     assert_allclose(slice_fluct_std1, sl_fluct_freq, rtol=0.5)
 
-    fa = ift.CorrelatedFieldMaker.make(0., offset_std_mean, .1, '', N, dofdex1)
-    fa.add_fluctuations(fsspace, astds[1], 1., 3.1, 1., .5, .1, -4, 1., 'freq',
+    fa = ift.CorrelatedFieldMaker.make(0., (offset_std_mean, .1), '', N, dofdex1)
+    fa.add_fluctuations(fsspace, (astds[1], 1.), (3.1, 1.), (.5, .1), (-4, 1.), 'freq',
                         dofdex=dofdex3)
     m = 3.
     x = fa.moment_slice_to_average(m)
-    fa.add_fluctuations(sspace, x, 1.5, 1.1, 2., 2.1, .5, -2, 1., 'spatial', 0,
+    fa.add_fluctuations(sspace, (x, 1.5), (1.1, 2.), (2.1, .5), (-2, 1.), 'spatial', 0,
                         dofdex=dofdex2)
     op = fa.finalize()
     em, estd = _stats(fa.slice_fluctuation(0), samples)
 
     assert_allclose(m, em, rtol=0.5)
-    assert_(op.target[-2] == sspace)
-    assert_(op.target[-1] == fsspace)
+    assert op.target[-2] == sspace
+    assert op.target[-1] == fsspace
 
     for ampl in fa.normalized_amplitudes:
         ift.extra.check_operator(ampl, 0.1*ift.from_random(ampl.domain), ntries=10)
@@ -124,71 +124,59 @@ def testAmplitudesInvariants(sspace, N):
 
 @pmp('seed', [42, 31])
 @pmp('domain', spaces)
-def test_complicated_vs_simple(seed, domain):
+@pmp('without', (('asperity', ), ('flexibility', ), ('flexibility', 'asperity')))
+def test_complicated_vs_simple(seed, domain, without):
     with ift.random.Context(seed):
         offset_mean = _rand()
-        offset_std_mean = _posrand()
-        offset_std_std = _posrand()
-        fluctuations_mean = _posrand()
-        fluctuations_stddev = _posrand()
-        flexibility_mean = _posrand()
-        flexibility_stddev = _posrand()
-        asperity_mean = _posrand()
-        asperity_stddev = _posrand()
-        loglogavgslope_mean = _posrand()
-        loglogavgslope_stddev = _posrand()
+        offset_std = _posrand(), _posrand()
+        fluctuations = _posrand(), _posrand()
+        if "flexibility" in without:
+            flexibility = None
+        else:
+            flexibility = _posrand(), _posrand()
+        if "asperity" in without:
+            asperity = None
+        else:
+            asperity = _posrand(), _posrand()
+        loglogavgslope = _posrand(), _posrand()
         prefix = 'foobar'
         hspace = domain.get_default_codomain()
-        scf = ift.SimpleCorrelatedField(
+        scf_args = (
             domain,
-            offset_mean, (offset_std_mean, offset_std_std),
-            (fluctuations_mean, fluctuations_stddev),
-            (flexibility_mean, flexibility_stddev),
-            (asperity_mean, asperity_stddev),
-            (loglogavgslope_mean, loglogavgslope_stddev),
-            prefix=prefix,
-            harmonic_partner=hspace)
-        cfm = ift.CorrelatedFieldMaker.make(offset_mean, offset_std_mean,
-                                            offset_std_std, prefix)
-        cfm.add_fluctuations(domain,
-                             fluctuations_mean,
-                             fluctuations_stddev,
-                             flexibility_mean,
-                             flexibility_stddev,
-                             asperity_mean,
-                             asperity_stddev,
-                             loglogavgslope_mean,
-                             loglogavgslope_stddev,
-                             prefix='',
+            offset_mean,
+            offset_std,
+            fluctuations,
+            flexibility,
+            asperity,
+            loglogavgslope
+        )
+        add_fluct_args = (
+            domain,
+            fluctuations,
+            flexibility,
+            asperity,
+            loglogavgslope
+        )
+        cfm = ift.CorrelatedFieldMaker.make(offset_mean, offset_std, prefix)
+        if asperity is not None and flexibility is None:
+            with pytest.raises(ValueError):
+                scf = ift.SimpleCorrelatedField(*scf_args, prefix=prefix,
+                                                harmonic_partner=hspace)
+            with pytest.raises(ValueError):
+                cfm.add_fluctuations(*add_fluct_args, prefix='',
+                                     harmonic_partner=hspace)
+            return
+        scf = ift.SimpleCorrelatedField(*scf_args, prefix=prefix,
+                                        harmonic_partner=hspace)
+        cfm.add_fluctuations(*add_fluct_args, prefix='',
                              harmonic_partner=hspace)
         inp = ift.from_random(scf.domain)
         op1 = cfm.finalize()
-        assert_(scf.domain is op1.domain)
+        assert scf.domain is op1.domain
         ift.extra.assert_allclose(scf(inp), op1(inp))
         ift.extra.check_operator(scf, inp, ntries=10)
 
         op1 = cfm.amplitude
         op0 = scf.amplitude
-        assert_(op0.domain is op1.domain)
+        assert op0.domain is op1.domain
         ift.extra.assert_allclose(op0.force(inp), op1.force(inp))
-
-
-@pmp('asperity', [None, (_posrand(), _posrand())])
-@pmp('flexibility', [None, (_posrand(), _posrand())])
-def test_simple_without_asp_fluct(asperity, flexibility):
-    domain = ift.RGSpace((4, 4), (0.123, 0.4))
-    offset_mean = _rand()
-    offset_std = _posrand(), _posrand()
-    fluctuations = _posrand(), _posrand()
-    loglogavgslope = _posrand(), _posrand()
-    prefix = 'foobar'
-    hspace = domain.get_default_codomain()
-    args = (domain, offset_mean, offset_std, fluctuations, flexibility,
-            asperity, loglogavgslope, prefix, hspace)
-    if asperity is not None and flexibility is None:
-        with pytest.raises(ValueError):
-            scf = ift.SimpleCorrelatedField(*args)
-    else:
-        scf = ift.SimpleCorrelatedField(*args)
-        inp = ift.from_random(scf.domain)
-        ift.extra.check_operator(scf, inp, ntries=10)
diff --git a/test/test_operators/test_einsum.py b/test/test_operators/test_einsum.py
index c426d5226c4ea30e2c6eefcbed03047015ee7992..d1360358a64993d40d9b3b803ba38254b444943a 100644
--- a/test/test_operators/test_einsum.py
+++ b/test/test_operators/test_einsum.py
@@ -17,10 +17,10 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import numpy as np
-from numpy.testing import assert_, assert_allclose
+from numpy.testing import assert_allclose
 
 import nifty7 as ift
-from nifty7.extra import check_operator, check_linear_operator
+from nifty7.extra import check_linear_operator, check_operator
 
 from ..common import list2fixture, setup_function, teardown_function
 
@@ -42,7 +42,7 @@ def test_linear_einsum_outer(space1, space2, dtype, n_invocations=10):
     ss = "i,ij,j->ij"
     key_order = ("dom01", "dom02")
     le = ift.LinearEinsum(space2, mf, ss, key_order=key_order)
-    assert_(check_linear_operator(le, domain_dtype=dtype, target_dtype=dtype) is None)
+    ift.myassert(check_linear_operator(le, domain_dtype=dtype, target_dtype=dtype) is None)
 
     le_ift = ift.DiagonalOperator(mf["dom01"], domain=mf_dom["dom02"], spaces=0) @ ift.DiagonalOperator(mf["dom02"])
     le_ift = le_ift @ ift.OuterProduct(ift.DomainTuple.make(mf_dom["dom02"][1]),
@@ -63,7 +63,7 @@ def test_linear_einsum_contraction(space1, space2, dtype, n_invocations=10):
     ss = "i,ij,j->i"
     key_order = ("dom01", "dom02")
     le = ift.LinearEinsum(space2, mf, ss, key_order=key_order)
-    assert_(check_linear_operator(le, domain_dtype=dtype, target_dtype=dtype) is None)
+    ift.myassert(check_linear_operator(le, domain_dtype=dtype, target_dtype=dtype) is None)
 
     le_ift = ift.ContractionOperator(mf_dom["dom02"], 1)
     le_ift = le_ift @ ift.DiagonalOperator(mf["dom01"], domain=mf_dom["dom02"], spaces=0)
@@ -116,7 +116,7 @@ def test_linear_einsum_transpose(space1, space2, dtype, n_invocations=10):
     mf = ift.MultiField.from_dict({})
     ss = "ij->ji"
     le = ift.LinearEinsum(dom, mf, ss)
-    assert_(check_linear_operator(le, domain_dtype=dtype, target_dtype=dtype) is None)
+    ift.myassert(check_linear_operator(le, domain_dtype=dtype, target_dtype=dtype) is None)
 
     # SwitchSpacesOperator is equivalent to LinearEinsum with "ij->ji"
     le_ift = _SwitchSpacesOperator(dom, 1)
diff --git a/test/test_operators/test_fft_operator.py b/test/test_operators/test_fft_operator.py
index bdfccc43bf3fbbf8f2c4fd9a08a203f12cad0a23..fffb7a4caf50c80a0a1ddbe2a395a9c78545328d 100644
--- a/test/test_operators/test_fft_operator.py
+++ b/test/test_operators/test_fft_operator.py
@@ -17,7 +17,7 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_allclose
+from numpy.testing import assert_allclose
 
 import nifty7 as ift
 
@@ -63,7 +63,7 @@ def test_fft1D(d, dtype, op):
 @pmp('nthreads', [-1, 1, 2, 3, 4])
 def test_fft2D(dim1, dim2, d1, d2, dtype, op, nthreads):
     ift.fft.set_nthreads(nthreads)
-    assert_(ift.fft.nthreads() == nthreads)
+    ift.myassert(ift.fft.nthreads() == nthreads)
     tol = _get_rtol(dtype)
     a = ift.RGSpace([dim1, dim2], distances=[d1, d2])
     b = ift.RGSpace(
diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py
index 9b166bb574e143dc5060fdb4d0d6f3812287f04b..f64562c5f105450154e7bc9abedcee8e6c29f983 100644
--- a/test/test_operators/test_nft.py
+++ b/test/test_operators/test_nft.py
@@ -17,9 +17,9 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_
 
 import nifty7 as ift
+
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
@@ -57,7 +57,7 @@ def test_gridding(nu, nv, N, eps):
     for i in range(N):
         dft += (
             vis[i]*np.exp(2j*np.pi*(x*uv[i, 0]*dstx + y*uv[i, 1]*dsty))).real
-    assert_(_l2error(dft, pynu) < eps)
+    ift.myassert(_l2error(dft, pynu) < eps)
 
 
 def test_cartesian():
diff --git a/test/test_operators/test_partial_multifield_insert.py b/test/test_operators/test_partial_multifield_insert.py
index e38dbbc7ccd980d95f8d0841172a9160fb5bc9c7..69d809245e39507e8febd8e2f15581f5f3b5e1f7 100644
--- a/test/test_operators/test_partial_multifield_insert.py
+++ b/test/test_operators/test_partial_multifield_insert.py
@@ -17,7 +17,7 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_allclose
+from numpy.testing import assert_allclose
 
 import nifty7 as ift
 
@@ -40,9 +40,9 @@ def test_part_mf_insert():
     op = a.partial_insert(b)
     fld = ift.from_random(op.domain, 'normal')
     ift.extra.check_operator(op, fld, ntries=ntries)
-    assert_(op.domain is ift.MultiDomain.union(
+    ift.myassert(op.domain is ift.MultiDomain.union(
         [op1.domain, op2.domain, op4.domain, op5.domain]))
-    assert_(op.target is ift.MultiDomain.union(
+    ift.myassert(op.target is ift.MultiDomain.union(
         [op1.target, op2.target, op3.target, op5.target]))
     x, y = fld.val, op(fld).val
     assert_allclose(y['a1'], x['a']*1.32)
diff --git a/test/test_operators/test_simplification.py b/test/test_operators/test_simplification.py
index f1fa7420c48718d7fe074e9344d4d2332cc67ed6..7f0c202f0e3a4c34f398dab1e5f344a76111f644 100644
--- a/test/test_operators/test_simplification.py
+++ b/test/test_operators/test_simplification.py
@@ -15,7 +15,7 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-from numpy.testing import assert_, assert_allclose, assert_raises
+from numpy.testing import assert_allclose, assert_raises
 
 import nifty7 as ift
 from nifty7.operators.simplify_for_const import ConstantOperator
@@ -26,7 +26,7 @@ def test_simplification():
     f1 = ift.full(dom, 2.)
     op = ift.FFTOperator(f1.domain["a"]).ducktape("a")
     _, op2 = op.simplify_for_constant_input(f1)
-    assert_(isinstance(op2, ConstantOperator))
+    ift.myassert(isinstance(op2, ConstantOperator))
     assert_allclose(op(f1).val, op2.force(f1).val)
 
     dom = {"a": ift.RGSpace(10), "b": ift.RGSpace(5)}
diff --git a/test/test_operators/test_value_inserter.py b/test/test_operators/test_value_inserter.py
index 1bfb45bd68d018d80013f25f6aa97d5cf3db4389..1e44d69837d4aacf20ba069ea3ba5f0946cfff67 100644
--- a/test/test_operators/test_value_inserter.py
+++ b/test/test_operators/test_value_inserter.py
@@ -17,9 +17,9 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_
 
 import nifty7 as ift
+
 from ..common import setup_function, teardown_function
 
 
@@ -38,5 +38,5 @@ def test_value_inserter(sp, seed):
         f = ift.from_random(op.domain, 'normal')
     inp = f.val
     ret = op(f).val
-    assert_(ret[ind] == inp)
-    assert_(np.sum(ret) == inp)
+    ift.myassert(ret[ind] == inp)
+    ift.myassert(np.sum(ret) == inp)
diff --git a/test/test_spaces/test_gl_space.py b/test/test_spaces/test_gl_space.py
index 7536e2f490802c83447f9db318d1e8c7cf4a3498..026b7b53f5a388d387945695c6f4638357346f26 100644
--- a/test/test_spaces/test_gl_space.py
+++ b/test/test_spaces/test_gl_space.py
@@ -19,10 +19,10 @@ import itertools
 
 import numpy as np
 import pytest
-from numpy.testing import (assert_, assert_almost_equal, assert_equal,
-                           assert_raises)
+from numpy.testing import assert_almost_equal, assert_equal, assert_raises
+
+from nifty7 import GLSpace, myassert
 
-from nifty7 import GLSpace
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
@@ -56,7 +56,7 @@ def get_dvol_configs():
 @pmp('attribute', ['nlat', 'nlon'])
 def test_property_ret_type(attribute):
     g = GLSpace(2)
-    assert_(isinstance(getattr(g, attribute), int))
+    myassert(isinstance(getattr(g, attribute), int))
 
 
 @pmp('nlat, nlon, expected', CONSTRUCTOR_CONFIGS)
diff --git a/test/test_spaces/test_hp_space.py b/test/test_spaces/test_hp_space.py
index 24dee06f1a6cad72799d84f0260c3b7e71b4d1fd..3eb17b27269c5c31c96c63e0bd066f71920dba78 100644
--- a/test/test_spaces/test_hp_space.py
+++ b/test/test_spaces/test_hp_space.py
@@ -17,10 +17,10 @@
 
 import numpy as np
 import pytest
-from numpy.testing import (assert_, assert_almost_equal, assert_equal,
-                           assert_raises)
+from numpy.testing import assert_almost_equal, assert_equal, assert_raises
+
+from nifty7 import HPSpace, myassert
 
-from nifty7 import HPSpace
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
@@ -49,7 +49,7 @@ CONSTRUCTOR_CONFIGS = [[
 
 def test_property_ret_type():
     x = HPSpace(2)
-    assert_(isinstance(getattr(x, 'nside'), int))
+    myassert(isinstance(getattr(x, 'nside'), int))
 
 
 @pmp('nside, expected', CONSTRUCTOR_CONFIGS)
diff --git a/test/test_spaces/test_interface.py b/test/test_spaces/test_interface.py
index 815c8b0d05f2f79e69d67f5ce1ab93dfe6889b5b..46ab83c67f90981b8c32e1baf925ef3189179d54 100644
--- a/test/test_spaces/test_interface.py
+++ b/test/test_spaces/test_interface.py
@@ -18,9 +18,9 @@
 from types import LambdaType
 
 import pytest
-from numpy.testing import assert_
 
 import nifty7 as ift
+
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
@@ -36,7 +36,7 @@ pmp = pytest.mark.parametrize
     ift.GLSpace(4)
 ])
 def test_property_ret_type(space, attr_expected_type):
-    assert_(
+    ift.myassert(
         isinstance(
             getattr(space, attr_expected_type[0]), attr_expected_type[1]))
 
@@ -46,7 +46,7 @@ def test_property_ret_type(space, attr_expected_type):
       ['get_fft_smoothing_kernel_function', 2.0, LambdaType]])
 @pmp('space', [ift.RGSpace(4, harmonic=True), ift.LMSpace(5)])
 def test_method_ret_type(space, method_expected_type):
-    assert_(
+    ift.myassert(
         type(
             getattr(space, method_expected_type[0])
             (*method_expected_type[1:-1])) is method_expected_type[-1])
diff --git a/test/test_spaces/test_lm_space.py b/test/test_spaces/test_lm_space.py
index e8ddc217679a168911cc3ef94ede46a7f4e289bb..12dbb3efc3ac2ad1fd13e8ac54243539a0fc3abd 100644
--- a/test/test_spaces/test_lm_space.py
+++ b/test/test_spaces/test_lm_space.py
@@ -17,9 +17,10 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
+from numpy.testing import assert_allclose, assert_equal, assert_raises
 
 import nifty7 as ift
+
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
@@ -67,7 +68,7 @@ def get_k_length_array_configs():
 
 @pmp('attribute', ['lmax', 'mmax', 'size'])
 def test_property_ret_type(attribute):
-    assert_(isinstance(getattr(ift.LMSpace(7, 5), attribute), int))
+    ift.myassert(isinstance(getattr(ift.LMSpace(7, 5), attribute), int))
 
 
 @pmp('lmax, mmax, expected', CONSTRUCTOR_CONFIGS)
diff --git a/test/test_spaces/test_power_space.py b/test/test_spaces/test_power_space.py
index a6eaf2007c526361fd7186b571f4aab0dc9d4a1e..a062eba6bc044ef24050aa0fcddfd32402170de2 100644
--- a/test/test_spaces/test_power_space.py
+++ b/test/test_spaces/test_power_space.py
@@ -19,9 +19,10 @@ from itertools import chain, product
 
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
+from numpy.testing import assert_allclose, assert_equal, assert_raises
 
 import nifty7 as ift
+
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
@@ -104,7 +105,7 @@ def k_lengths_configs():
 def test_property_ret_type(attribute, expected_type):
     r = ift.RGSpace((4, 4), harmonic=True)
     p = ift.PowerSpace(r)
-    assert_(isinstance(getattr(p, attribute), expected_type))
+    ift.myassert(isinstance(getattr(p, attribute), expected_type))
 
 
 @pmp('harmonic_partner, binbounds, nbin, logarithmic', CONSISTENCY_CONFIGS)
diff --git a/test/test_spaces/test_rg_space.py b/test/test_spaces/test_rg_space.py
index 788970a9528f56008b9bbf2257fe486c66612ec5..8343a170ee21048510574b17e88115aa18933bd1 100644
--- a/test/test_spaces/test_rg_space.py
+++ b/test/test_spaces/test_rg_space.py
@@ -17,9 +17,10 @@
 
 import numpy as np
 import pytest
-from numpy.testing import assert_, assert_allclose, assert_equal
+from numpy.testing import assert_allclose, assert_equal
 
 import nifty7 as ift
+
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
@@ -77,7 +78,7 @@ def get_dvol_configs():
 @pmp('attribute, expected_type', [['distances', tuple]])
 def test_property_ret_type(attribute, expected_type):
     x = ift.RGSpace(1)
-    assert_(isinstance(getattr(x, attribute), expected_type))
+    ift.myassert(isinstance(getattr(x, attribute), expected_type))
 
 
 @pmp('shape, distances, harmonic, expected', CONSTRUCTOR_CONFIGS)