diff --git a/.gitignore b/.gitignore
index 628b6dddd1aa886fb49b022acba96dc8d1425957..211f64ebfe248b222c45d92b9f71fd5e9327d5a9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,6 +10,7 @@ setup.cfg
 .document
 .svn/
 *.csv
+.pytest_cache/
 
 # from https://github.com/github/gitignore/blob/master/Python.gitignore
 
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index d385dd1b36e7c536e31aa56e2e001be54f235047..415bd85c53cfba7f021240c1e7734aa0a9a47cce 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -34,18 +34,22 @@ build_docker_from_cache:
     - docker build -t $CONTAINER_TEST_IMAGE .
     - docker push $CONTAINER_TEST_IMAGE
 
-test:
+test_serial:
   stage: test
-  variables:
-    OMPI_MCA_btl_vader_single_copy_mechanism: none
   script:
-    - mpiexec -n 2 --bind-to none pytest-3 -q test
     - pytest-3 -q --cov=nifty5 test
     - >
       python3 -m coverage report --omit "*plotting*,*distributed_do*"
     - >
       python3 -m coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }'
 
+test_mpi:
+  stage: test
+  variables:
+    OMPI_MCA_btl_vader_single_copy_mechanism: none
+  script:
+    - mpiexec -n 2 --bind-to none pytest-3 -q test
+
 pages:
   stage: release
   before_script:
diff --git a/Dockerfile b/Dockerfile
index 93c6276b9ab30b38ede21d6f72bef8c9c540e008..aea6272a0769e62457ae3daf7d1b96e65c59b9fe 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -7,9 +7,9 @@ RUN apt-get update && apt-get install -y \
     libfftw3-dev \
     python3 python3-pip python3-dev python3-future python3-scipy cython3 \
     # Documentation build dependencies
-    python3-sphinx python3-sphinx-rtd-theme python3-numpydoc \
+    python3-sphinx python3-sphinx-rtd-theme \
     # Testing dependencies
-    python3-coverage python3-parameterized python3-pytest python3-pytest-cov \
+    python3-coverage python3-pytest python3-pytest-cov \
     # Optional NIFTy dependencies
     openmpi-bin libopenmpi-dev python3-mpi4py \
     # Packages needed for NIFTy
diff --git a/docs/generate.sh b/docs/generate.sh
index 4d9ec2932e651854b10c2e018f06c75ddd2aa9e4..2a77e15f649c8c598f27f3173c18be81750510c8 100755
--- a/docs/generate.sh
+++ b/docs/generate.sh
@@ -1,3 +1,2 @@
-rm -rf docs/build docs/source/mod
 sphinx-apidoc -l -e -d 2 -o docs/source/mod nifty5
 sphinx-build -b html docs/source/ docs/build/
diff --git a/docs/source/code.rst b/docs/source/code.rst
index ae1900660e4d54be1cf6728621a300868b8c7d16..88be215ff7e64fd74aabace144ce6a17bed4cbf2 100644
--- a/docs/source/code.rst
+++ b/docs/source/code.rst
@@ -1,4 +1,3 @@
-.. currentmodule:: nifty5
 
 =============
 Code Overview
@@ -37,9 +36,12 @@ Domains
 Abstract base class
 -------------------
 
+.. currentmodule:: nifty5.domains.domain
+
 One of the fundamental building blocks of the NIFTy5 framework is the *domain*.
-Its required capabilities are expressed by the abstract :class:`Domain` class.
+Its required capabilities are expressed by the abstract :py:class:`Domain` class.
 A domain must be able to answer the following queries:
+m
 
 - its total number of data entries (pixels), which is accessible via the
   :attr:`~Domain.size` property
@@ -51,6 +53,8 @@ A domain must be able to answer the following queries:
 Unstructured domains
 --------------------
 
+.. currentmodule:: nifty5.domains.unstructured_domain
+
 Domains can be either *structured* (i.e. there is geometrical information
 associated with them, like position in space and volume factors),
 or *unstructured* (meaning that the data points have no associated manifold).
@@ -62,6 +66,8 @@ Unstructured domains can be described by instances of NIFTy's
 Structured domains
 ------------------
 
+.. currentmodule:: nifty5.domains.structured_domain
+
 In contrast to unstructured domains, these domains have an assigned geometry.
 NIFTy requires them to provide the volume elements of their grid cells.
 The additional methods are specified in the abstract class
@@ -81,15 +87,17 @@ The additional methods are specified in the abstract class
 
 NIFTy comes with several concrete subclasses of :class:`StructuredDomain`:
 
-- :class:`RGSpace` represents a regular Cartesian grid with an arbitrary
+.. currentmodule:: nifty5.domains
+
+- :class:`rg_space.RGSpace` represents a regular Cartesian grid with an arbitrary
   number of dimensions, which is supposed to be periodic in each dimension.
-- :class:`HPSpace` and :class:`GLSpace` describe pixelisations of the
-  2-sphere; their counterpart in harmonic space is :class:`LMSpace`, which
+- :class:`hp_space.HPSpace` and :class:`gl_space.GLSpace` describe pixelisations of the
+  2-sphere; their counterpart in harmonic space is :class:`lm_space.LMSpace`, which
   contains spherical harmonic coefficients.
-- :class:`PowerSpace` is used to describe one-dimensional power spectra.
+- :class:`power_space.PowerSpace` is used to describe one-dimensional power spectra.
 
-Among these, :class:`RGSpace` can be harmonic or not (depending on constructor arguments), :class:`GLSpace`, :class:`HPSpace`, and :class:`PowerSpace` are
-pure position domains (i.e. nonharmonic), and :class:`LMSpace` is always
+Among these, :class:`rg_space.RGSpace` can be harmonic or not (depending on constructor arguments), :class:`gl_space.GLSpace`, :class:`hp_space.HPSpace`, and :class:`power_space.PowerSpace` are
+pure position domains (i.e. nonharmonic), and :class:`lm_space.LMSpace` is always
 harmonic.
 
 
@@ -158,7 +166,7 @@ be extracted first, then changed, and a new field has to be created from the
 result.
 
 Fields defined on a MultiDomain
-------------------------------
+-------------------------------
 
 The :class:`MultiField` class can be seen as a dictionary of individual
 :class:`Field` s, each identified by a name, which is defined on a
@@ -300,7 +308,7 @@ As an example one may consider the following combination of ``x``, which is an o
 
 
 Basic operators
-------------
+---------------
 # FIXME All this is outdated!
 
 Basic operator classes provided by NIFTy are
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 8174948041b6dd38fb841f0842d2264c1ed95bbb..1fac727c236f387f6ec43dd8f1f1e03d2d31158a 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -1,24 +1,16 @@
 import nifty5
-import sphinx_rtd_theme
-
-napoleon_google_docstring = False
-napoleon_numpy_docstring = True
-napoleon_use_ivar = True
-napoleon_use_param = False
-napoleon_use_keyword = False
-
-autodoc_member_order = 'groupwise'
-numpydoc_show_inherited_class_members = False
-numpydoc_class_members_toctree = False
 
 extensions = [
-    'sphinx.ext.autodoc', 'numpydoc', 'sphinx.ext.autosummary',
-    'sphinx.ext.napoleon', 'sphinx.ext.imgmath', 'sphinx.ext.viewcode'
+    'sphinx.ext.napoleon',  # Support for NumPy and Google style docstrings
+    'sphinx.ext.imgmath',  # Render math as images
+    'sphinx.ext.viewcode'  # Add links to highlighted source code
 ]
-templates_path = ['_templates']
-source_suffix = '.rst'
 master_doc = 'index'
 
+napoleon_google_docstring = False
+napoleon_numpy_docstring = True
+napoleon_use_ivar = True
+
 project = u'NIFTy5'
 copyright = u'2013-2019, Max-Planck-Society'
 author = u'Martin Reinecke'
@@ -29,19 +21,6 @@ version = release[:-2]
 language = None
 exclude_patterns = []
 add_module_names = False
-pygments_style = 'sphinx'
-todo_include_todos = True
 
 html_theme = "sphinx_rtd_theme"
-html_theme_options = {
-    'collapse_navigation': False,
-    'display_version': False,
-}
-html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 html_logo = 'nifty_logo_black.png'
-html_static_path = []
-html_last_updated_fmt = '%b %d, %Y'
-html_domain_indices = False
-html_use_index = False
-html_show_sourcelink = False
-htmlhelp_basename = 'NIFTydoc'
diff --git a/docs/source/ift.rst b/docs/source/ift.rst
index bab6307688c24c454eaf7381ce1215ef1bd3fa7e..883ff2d61a4149ac8dfff7f2a7e908cace8a276f 100644
--- a/docs/source/ift.rst
+++ b/docs/source/ift.rst
@@ -104,6 +104,7 @@ with :math:`{R}` the measurement response, which maps the continous signal field
 
 This is called a free theory, as the information Hamiltonian
 associate professor
+
 .. math::
 
     \mathcal{H}(d,s)= -\log \mathcal{P}(d,s)= \frac{1}{2} s^\dagger S^{-1} s + \frac{1}{2} (d-R\,s)^\dagger N^{-1} (d-R\,s) + \mathrm{const}
@@ -179,23 +180,22 @@ NIFTy takes advantage of this formulation in several ways:
 
 The reconstruction of a non-Gaussian signal with unknown covarinance from a non-trivial (tomographic) response is demonstrated in demos/getting_started_3.py. Here, the uncertainty of the field and the power spectrum of its generating process are probed via posterior samples provided by the MGVI algorithm.
 
-+-------------------------------------------------+
-| .. image:: images/getting_started_3_setup.png   |
-|     :width:  30 %                               |
-+-------------------------------------------------+
-| .. image:: images/getting_started_3_results.png |
-|     :width:  30 %                               |
-+-------------------------------------------------+
-| Output of tomography demo getting_started_3.py. |
-| **Top row:** Non-Gaussian signal field,         |
-| data backprojected into the image domain, power |
-| spectrum of underlying Gausssian process.       |
-| **Bottom row:** Posterior mean field signal     |
-| reconstruction, its uncertainty, and the power  |
-| spectrum of the process for different posterior |
-| samples in comparison to the correct one (thick |
-| orange line).                                   |
-+-------------------------------------------------+
-
-
-
++----------------------------------------------------+
+| **Output of tomography demo getting_started_3.py** |
++----------------------------------------------------+
+| .. image:: images/getting_started_3_setup.png      |
+|                                                    |
++----------------------------------------------------+
+| Non-Gaussian signal field,                         |
+| data backprojected into the image domain, power    |
+| spectrum of underlying Gausssian process.          |
++----------------------------------------------------+
+| .. image:: images/getting_started_3_results.png    |
+|                                                    |
++----------------------------------------------------+
+| Posterior mean field signal                        |
+| reconstruction, its uncertainty, and the power     |
+| spectrum of the process for different posterior    |
+| samples in comparison to the correct one (thick    |
+| orange line).                                      |
++----------------------------------------------------+
diff --git a/nifty5/data_objects/distributed_do.py b/nifty5/data_objects/distributed_do.py
index bbd56e3cc561b0d76983c69f71aa65bfd3505b8c..4f93fb09751d443e060bff2d4a5caf9c7b5d1c2a 100644
--- a/nifty5/data_objects/distributed_do.py
+++ b/nifty5/data_objects/distributed_do.py
@@ -17,6 +17,7 @@
 
 import sys
 from functools import reduce
+
 import numpy as np
 from mpi4py import MPI
 
diff --git a/nifty5/data_objects/numpy_do.py b/nifty5/data_objects/numpy_do.py
index 4f30a27ddc23a26ea847432069cdbc782ebe8634..8160f3b34163283a146c977cdd14d4aa36f0d370 100644
--- a/nifty5/data_objects/numpy_do.py
+++ b/nifty5/data_objects/numpy_do.py
@@ -18,11 +18,10 @@
 # Data object module for NIFTy that uses simple numpy ndarrays.
 
 import numpy as np
-from numpy import empty, empty_like, exp, full, log
+from numpy import absolute, clip, cos, cosh, empty, empty_like, exp, full, log
 from numpy import ndarray as data_object
-from numpy import ones, sqrt, tanh, vdot, zeros
-from numpy import sin, cos, tan, sinh, cosh, sinc
-from numpy import absolute, sign, clip
+from numpy import ones, sign, sin, sinc, sinh, sqrt, tan, tanh, vdot, zeros
+
 from .random import Random
 
 __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full",
diff --git a/nifty5/library/amplitude_operator.py b/nifty5/library/amplitude_operator.py
index dcca12300b5936ec8dc1fa770858fd51fcc50be6..a3b2fbfab5b3646a6bca6fb9c4326fb6dad1b836 100644
--- a/nifty5/library/amplitude_operator.py
+++ b/nifty5/library/amplitude_operator.py
@@ -111,23 +111,27 @@ def _SlopePowerSpectrum(logk_space, sm, sv, im, iv):
 
 def AmplitudeOperator(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
                       keys=['tau', 'phi'], zero_mode=True):
-    '''
+    ''' Operator for parametrizing smooth power spectra.
+
     Computes a smooth power spectrum.
     Output is defined on a PowerSpace.
 
     Parameters
     ----------
-
-    Npixdof : #pix in dof_space
-
-    ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
-                        eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
-                        a = ceps_a,  k0 = ceps_k0
-
-    sm, sv : slope_mean = expected exponent of power law (e.g. -4),
-                slope_variance (default=1)
-
-    im, iv : y-intercept_mean, y-intercept_variance  of power_slope
+    Npixdof : int
+        #pix in dof_space
+    ceps_a : float
+        Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a,  k0 = ceps_k0
+    ceps_k0 : float
+        Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a,  k0 = ceps_k0
+    sm : float
+        slope_mean = expected exponent of power law (e.g. -4)
+    sv : float
+        slope_variance (default=1)
+    im : float
+        y-intercept_mean
+    iv : float
+        y-intercept_variance  of power_slope
     '''
 
     from ..operators.exp_transform import ExpTransform
diff --git a/nifty5/library/los_response.py b/nifty5/library/los_response.py
index 9c7ab404f2bd3be2ffdadb586be33e121020192f..400eadcec66347d38d447235b65cc056e643819b 100644
--- a/nifty5/library/los_response.py
+++ b/nifty5/library/los_response.py
@@ -28,14 +28,14 @@ from ..field import Field
 from ..operators.linear_operator import LinearOperator
 
 
-def _gaussian_error_function(x):
-    return 0.5/erfc(x*np.sqrt(2.))
+def _gaussian_sf(x):
+    return 0.5*erfc(x/np.sqrt(2.))
 
 
-def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
+def _comp_traverse(start, end, shp, dist, lo, mid, hi, sig, erf):
     ndim = start.shape[0]
     nlos = start.shape[1]
-    inc = np.full(len(shp), 1)
+    inc = np.full(len(shp), 1, dtype=np.int64)
     for i in range(-2, -len(shp)-1, -1):
         inc[i] = inc[i+1]*shp[i+1]
 
@@ -59,7 +59,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
         dmin += 1e-7
         dmax -= 1e-7
         if dmin >= dmax:  # no intersection
-            out[i] = (np.full(0, 0), np.full(0, 0.))
+            out[i] = (np.full(0, 0, dtype=np.int64), np.full(0, 0.))
             continue
         # determine coordinates of first cell crossing
         c_first = np.ceil(start[:, i]+dir*dmin)
@@ -75,7 +75,7 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
                 tmp = np.arange(start=c_first[j], stop=dmax,
                                 step=abs(1./dir[j]))
                 cdist = np.append(cdist, tmp)
-                add = np.append(add, np.full(len(tmp), step))
+                add = np.append(add, np.full(len(tmp), step, dtype=np.int64))
         idx = np.argsort(cdist)
         cdist = cdist[idx]
         add = add[idx]
@@ -85,21 +85,19 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
         cdist *= corfac
         wgt = np.diff(cdist)
         mdist = 0.5*(cdist[:-1]+cdist[1:])
-        wgt = apply_erf(wgt, mdist, lo[i], mid[i], hi[i], erf)
+        wgt = apply_erf(wgt, mdist, lo[i], mid[i], hi[i], sig[i], erf)
         add = np.append(pos1, add)
         add = np.cumsum(add)
         out[i] = (add, wgt)
     return out
 
 
-def apply_erf(wgt, dist, lo, mid, hi, erf):
+def apply_erf(wgt, dist, lo, mid, hi, sig, erf):
     wgt = wgt.copy()
     mask = dist > hi
     wgt[mask] = 0.
-    mask = (dist > mid) & (dist <= hi)
-    wgt[mask] *= erf((dist[mask]-mid)/(hi-mid))
-    mask = (dist <= mid) & (dist > lo)
-    wgt[mask] *= erf((dist[mask]-mid)/(mid-lo))
+    mask = (dist > lo) & (dist <= hi)
+    wgt[mask] *= erf((-1/dist[mask]+1/mid)/sig)
     return wgt
 
 
@@ -119,17 +117,32 @@ class LOSResponse(LinearOperator):
         of sight. The first dimension must have as many entries as `domain`
         has dimensions. The second dimensions must be identical for both arrays
         and indicated the total number of lines of sight.
-    sigmas_low, sigmas_up : numpy.ndarray(float) (optional)
-        For expert use. If unsure, leave blank.
+    sigmas: numpy.ndarray(float) (optional)
+        If this is not None, the inverse of the lengths of the LOSs are assumed
+        to be Gaussian distributed with these sigmas. The start point will
+        remain the same, but the endpoint is assumed to be unknown.
+        This is a typical statistical model for astrophysical parallaxes.
+        The LOS response then returns the expected integral
+        over the input given that the length of the LOS is unknown and
+        therefore the result is averaged over different endpoints.
+        default: None
+    truncation: float (optional)
+        Use only if the sigmas keyword argument is used!
+        This truncates the probability of the endpoint lying more sigmas away
+        than the truncation. Used to speed up computation and to avoid negative
+        distances. It should hold that `1./(1./length-sigma*truncation)>0`
+        for all lengths of the LOSs and all corresponding sigma of sigmas.
+        If unsure, leave blank.
+        default: 3.
 
     Notes
     -----
-    `starts, `ends`, `sigmas_low`, and `sigmas_up` have to be identical on
+    `starts, `ends`, `sigmas`, and `truncation` have to be identical on
     every calling MPI task (i.e. the full LOS information has to be provided on
     every task).
     """
 
-    def __init__(self, domain, starts, ends, sigmas_low=None, sigmas_up=None):
+    def __init__(self, domain, starts, ends, sigmas=None, truncation=3.):
         self._domain = DomainTuple.make(domain)
         self._capability = self.TIMES | self.ADJOINT_TIMES
 
@@ -141,20 +154,15 @@ class LOSResponse(LinearOperator):
         starts = np.array(starts)
         nlos = starts.shape[1]
         ends = np.array(ends)
-        if sigmas_low is None:
-            sigmas_low = np.zeros(nlos, dtype=np.float32)
-        if sigmas_up is None:
-            sigmas_up = np.zeros(nlos, dtype=np.float32)
-        sigmas_low = np.array(sigmas_low)
-        sigmas_up = np.array(sigmas_up)
+        if sigmas is None:
+            sigmas = np.zeros(nlos, dtype=np.float32)
+        sigmas = np.array(sigmas)
         if starts.shape[0] != ndim:
             raise TypeError("dimension mismatch")
-        if nlos != sigmas_low.shape[0]:
+        if nlos != sigmas.shape[0]:
             raise TypeError("dimension mismatch")
         if starts.shape != ends.shape:
             raise TypeError("dimension mismatch")
-        if sigmas_low.shape != sigmas_up.shape:
-            raise TypeError("dimension mismatch")
 
         self._local_shape = dobj.local_shape(self.domain[0].shape)
         local_zero_point = (np.array(
@@ -164,7 +172,11 @@ class LOSResponse(LinearOperator):
         diffs = ends-starts
         difflen = np.linalg.norm(diffs, axis=0)
         diffs /= difflen
-        real_ends = ends + sigmas_up*diffs
+        real_distances = 1./(1./difflen - truncation*sigmas)
+        if np.any(real_distances < 0):
+            raise ValueError("parallax error truncation to high: "
+                             "getting negative distances")
+        real_ends = starts + diffs*real_distances
         lzp = local_zero_point.reshape((-1, 1))
         dist = np.array(self.domain[0].distances).reshape((-1, 1))
         localized_pixel_starts = (starts-lzp)/dist + 0.5
@@ -175,8 +187,11 @@ class LOSResponse(LinearOperator):
                              localized_pixel_ends,
                              self._local_shape,
                              np.array(self.domain[0].distances),
-                             difflen-sigmas_low, difflen, difflen+sigmas_up,
-                             _gaussian_error_function)
+                             1./(1./difflen+truncation*sigmas),
+                             difflen,
+                             1./(1./difflen-truncation*sigmas),
+                             sigmas,
+                             _gaussian_sf)
 
         boxsz = 16
         nlos = len(w_i)
@@ -202,7 +217,7 @@ class LOSResponse(LinearOperator):
                 tmp += (fullidx[j]//boxsz)*fct
                 fct *= self._local_shape[j]
             tmp += cnt/float(nlos)
-            tmp += iarr[ofs:ofs+nval]/float(nlos*npix)
+            tmp += iarr[ofs:ofs+nval]/(float(nlos)*float(npix))
             pri[ofs:ofs+nval] = tmp
             ofs += nval
             cnt += 1
diff --git a/nifty5/linearization.py b/nifty5/linearization.py
index 3b7c5b486db5f7572d07f6512fd907e326433fb4..4d75ec6b6199f4ca96168a5490a5bde1002c4617 100644
--- a/nifty5/linearization.py
+++ b/nifty5/linearization.py
@@ -110,11 +110,11 @@ class Linearization(object):
 
     def __truediv__(self, other):
         if isinstance(other, Linearization):
-            return self.__mul__(other.inverse())
+            return self.__mul__(other.one_over())
         return self.__mul__(1./other)
 
     def __rtruediv__(self, other):
-        return self.inverse().__mul__(other)
+        return self.one_over().__mul__(other)
 
     def __pow__(self, power):
         if not np.isscalar(power):
@@ -122,9 +122,6 @@ class Linearization(object):
         return self.new(self._val**power,
                         makeOp(self._val**(power-1)).scale(power)(self._jac))
 
-    def inverse(self):
-        return self.new(1./self._val, makeOp(-1./(self._val**2))(self._jac))
-
     def __mul__(self, other):
         from .sugar import makeOp
         if isinstance(other, Linearization):
diff --git a/nifty5/operators/domain_tuple_field_inserter.py b/nifty5/operators/domain_tuple_field_inserter.py
index c190cfb20a4696cd66046cb17d874ce321fc5bc6..74e5d555e0b49f4082647ee6ae6c192e9aa736d3 100644
--- a/nifty5/operators/domain_tuple_field_inserter.py
+++ b/nifty5/operators/domain_tuple_field_inserter.py
@@ -23,18 +23,18 @@ from .linear_operator import LinearOperator
 
 
 class DomainTupleFieldInserter(LinearOperator):
-    def __init__(self, domain, new_space, index, position):
-        '''Writes the content of a field into one slice of a DomainTuple.
+    '''Writes the content of a field into one slice of a DomainTuple.
 
-        Parameters
-        ----------
-        domain : Domain, tuple of Domain or DomainTuple
-        new_space : Domain, tuple of Domain or DomainTuple
-        index : Integer
-            Index at which new_space shall be added to domain.
-        position : tuple
-            Slice in new_space in which the input field shall be written into.
-        '''
+    Parameters
+    ----------
+    domain : Domain, tuple of Domain or DomainTuple
+    new_space : Domain, tuple of Domain or DomainTuple
+    index : Integer
+        Index at which new_space shall be added to domain.
+    position : tuple
+        Slice in new_space in which the input field shall be written into.
+    '''
+    def __init__(self, domain, new_space, index, position):
         self._domain = DomainTuple.make(domain)
         tgt = list(self.domain)
         tgt.insert(index, new_space)
diff --git a/test/common.py b/test/common.py
index e3bc49c645f5ab805c01117866d7f09bdf4557be..9b95150d7ad0ebf1ea0922348fae3fba3a791fd1 100644
--- a/test/common.py
+++ b/test/common.py
@@ -15,20 +15,12 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-from builtins import str
+import pytest
 
-import numpy as np
-from parameterized import parameterized
 
-np.seterr(all='raise', under='ignore')
+def list2fixture(lst):
+    @pytest.fixture(params=lst)
+    def myfixture(request):
+        return request.param
 
-
-def _custom_name_func(testcase_func, param_num, param):
-    return "{}_{}".format(
-        testcase_func.__name__,
-        parameterized.to_safe_name("_".join(str(x) for x in param.args)),
-    )
-
-
-def expand(*args, **kwargs):
-    return parameterized.expand(*args, func=_custom_name_func, **kwargs)
+    return myfixture
diff --git a/test/test_energies/test_consistency.py b/test/test_energies/test_consistency.py
deleted file mode 100644
index 906e61ce5aaa37d32828fab7adaa7b152e049568..0000000000000000000000000000000000000000
--- a/test/test_energies/test_consistency.py
+++ /dev/null
@@ -1,118 +0,0 @@
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-# Copyright(C) 2013-2019 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
-import numpy as np
-
-
-class Energy_Tests(unittest.TestCase):
-    def make_operator(self, **kwargs):
-        np.random.seed(kwargs['seed'])
-        S = ift.ScalingOperator(1., kwargs['space'])
-        s = S.draw_sample()
-        return ift.MultiField.from_dict({kwargs['space_key']: s})
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testGaussian(self, space, seed):
-        op = self.make_operator(
-            space_key='s1', space=space, seed=seed)['s1']
-        energy = ift.GaussianEnergy(domain=space)
-        ift.extra.check_value_gradient_consistency(energy, op)
-
-#     @expand(product(
-#         [ift.GLSpace(15),
-#          ift.RGSpace(64, distances=.789),
-#          ift.RGSpace([32, 32], distances=.789)],
-#         [4, 78, 23]
-#         ))
-#     def testQuadratic(self, type1, space, seed):
-#         np.random.seed(seed)
-#         S = ift.ScalingOperator(1., space)
-#         s = [S.draw_sample() for _ in range(3)]
-#         energy = ift.QuadraticEnergy(s[0], ift.makeOp(s[1]), s[2])
-#         ift.extra.check_value_gradient_consistency(energy)
-
-    @expand(
-        product([
-            ift.GLSpace(15),
-            ift.RGSpace(64, distances=.789),
-            ift.RGSpace([32, 32], distances=.789)
-        ], [4, 78, 23]))
-    def testInverseGammaLikelihood(self, space, seed):
-        op = self.make_operator(space_key='s1', space=space, seed=seed)['s1']
-        op = op.exp()
-        d = np.random.normal(10, size=space.shape)**2
-        d = ift.Field.from_global_data(space, d)
-        energy = ift.InverseGammaLikelihood(d)
-        ift.extra.check_value_gradient_consistency(energy, op, tol=1e-7)
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testPoissonian(self, space, seed):
-        op = self.make_operator(
-            space_key='s1', space=space, seed=seed)['s1']
-        op = op.exp()
-        d = np.random.poisson(120, size=space.shape)
-        d = ift.Field.from_global_data(space, d)
-        energy = ift.PoissonianEnergy(d)
-        ift.extra.check_value_gradient_consistency(energy, op, tol=1e-7)
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testHamiltonian_and_KL(self, space, seed):
-        op = self.make_operator(
-            space_key='s1', space=space, seed=seed)['s1']
-        op = op.exp()
-        lh = ift.GaussianEnergy(domain=space)
-        hamiltonian = ift.Hamiltonian(lh)
-        ift.extra.check_value_gradient_consistency(hamiltonian, op)
-        S = ift.ScalingOperator(1., space)
-        samps = [S.draw_sample() for i in range(3)]
-        kl = ift.SampledKullbachLeiblerDivergence(hamiltonian, samps)
-        ift.extra.check_value_gradient_consistency(kl, op)
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testBernoulli(self, space, seed):
-        op = self.make_operator(
-            space_key='s1', space=space, seed=seed)['s1']
-        op = op.sigmoid()
-        d = np.random.binomial(1, 0.1, size=space.shape)
-        d = ift.Field.from_global_data(space, d)
-        energy = ift.BernoulliEnergy(d)
-        ift.extra.check_value_gradient_consistency(energy, op, tol=1e-6)
diff --git a/test/test_energies/test_map.py b/test/test_energies/test_map.py
deleted file mode 100644
index 7953b07405abe09eec57088325ca3d59a79453b0..0000000000000000000000000000000000000000
--- a/test/test_energies/test_map.py
+++ /dev/null
@@ -1,75 +0,0 @@
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-# Copyright(C) 2013-2019 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
-import numpy as np
-
-
-def _flat_PS(k):
-    return np.ones_like(k)
-
-
-class Energy_Tests(unittest.TestCase):
-    @expand(product([ift.GLSpace(15),
-                     ift.RGSpace(64, distances=.789),
-                     ift.RGSpace([32, 32], distances=.789)],
-                    ["tanh", "exp", ""],
-                    [1, 1e-2, 1e2],
-                    [4, 78, 23]))
-    def testGaussianEnergy(self, space, nonlinearity, noise, seed):
-        np.random.seed(seed)
-        dim = len(space.shape)
-        hspace = space.get_default_codomain()
-        ht = ift.HarmonicTransformOperator(hspace, target=space)
-        binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
-        pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
-        xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
-        xi0_var = ift.Linearization.make_var(xi0)
-
-        def pspec(k): return 1 / (1 + k**2)**dim
-        pspec = ift.PS_field(pspace, pspec)
-        A = Dist(ift.sqrt(pspec))
-        N = ift.ScalingOperator(noise, space)
-        n = N.draw_sample()
-        s = ht(ift.makeOp(A)(xi0_var))
-        R = ift.ScalingOperator(10., space)
-
-        def d_model():
-            if nonlinearity == "":
-                return R(ht(ift.makeOp(A)))
-            else:
-                tmp = ht(ift.makeOp(A))
-                nonlin = getattr(tmp, nonlinearity)()
-                return R(nonlin)
-
-        d = d_model()(xi0) + n
-
-        if noise == 1:
-            N = None
-
-        energy = ift.GaussianEnergy(d, N)(d_model())
-        if nonlinearity == "":
-            ift.extra.check_value_gradient_metric_consistency(
-                energy, xi0, ntries=10)
-        else:
-            ift.extra.check_value_gradient_consistency(
-                energy, xi0, ntries=10, tol=5e-8)
diff --git a/test/test_energy_gradients.py b/test/test_energy_gradients.py
new file mode 100644
index 0000000000000000000000000000000000000000..9cbe7105769c9e8f39f26eb3a54b8c20fb674ac0
--- /dev/null
+++ b/test/test_energy_gradients.py
@@ -0,0 +1,86 @@
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+import pytest
+
+import nifty5 as ift
+from itertools import product
+
+# Currently it is not possible to parametrize fixtures. But this will
+# hopefully be fixed in the future.
+# https://docs.pytest.org/en/latest/proposals/parametrize_with_fixtures.html
+
+SPACES = [
+    ift.GLSpace(15),
+    ift.RGSpace(64, distances=.789),
+    ift.RGSpace([32, 32], distances=.789)
+]
+SEEDS = [4, 78, 23]
+PARAMS = product(SEEDS, SPACES)
+
+
+@pytest.fixture(params=PARAMS)
+def field(request):
+    np.random.seed(request.param[0])
+    S = ift.ScalingOperator(1., request.param[1])
+    s = S.draw_sample()
+    return ift.MultiField.from_dict({'s1': s})['s1']
+
+
+def test_gaussian(field):
+    energy = ift.GaussianEnergy(domain=field.domain)
+    ift.extra.check_value_gradient_consistency(energy, field)
+
+
+def test_inverse_gamma(field):
+    field = field.exp()
+    space = field.domain
+    d = np.random.normal(10, size=space.shape)**2
+    d = ift.Field.from_global_data(space, d)
+    energy = ift.InverseGammaLikelihood(d)
+    ift.extra.check_value_gradient_consistency(energy, field, tol=1e-7)
+
+
+def testPoissonian(field):
+    field = field.exp()
+    space = field.domain
+    d = np.random.poisson(120, size=space.shape)
+    d = ift.Field.from_global_data(space, d)
+    energy = ift.PoissonianEnergy(d)
+    ift.extra.check_value_gradient_consistency(energy, field, tol=1e-7)
+
+
+def test_hamiltonian_and_KL(field):
+    field = field.exp()
+    space = field.domain
+    lh = ift.GaussianEnergy(domain=space)
+    hamiltonian = ift.Hamiltonian(lh)
+    ift.extra.check_value_gradient_consistency(hamiltonian, field)
+    S = ift.ScalingOperator(1., space)
+    samps = [S.draw_sample() for i in range(3)]
+    kl = ift.SampledKullbachLeiblerDivergence(hamiltonian, samps)
+    ift.extra.check_value_gradient_consistency(kl, field)
+
+
+def test_bernoulli(field):
+    field = field.sigmoid()
+    space = field.domain
+    d = np.random.binomial(1, 0.1, size=space.shape)
+    d = ift.Field.from_global_data(space, d)
+    energy = ift.BernoulliEnergy(d)
+    ift.extra.check_value_gradient_consistency(energy, field, tol=1e-6)
diff --git a/test/test_field.py b/test/test_field.py
index 59d380db66c21c5edaf0b87dd70a4b0c436072f3..285bd07d33c591cbdff0c56e020c53729ea10e56 100644
--- a/test/test_field.py
+++ b/test/test_field.py
@@ -15,290 +15,311 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
 import numpy as np
+import pytest
 from numpy.testing import assert_allclose, assert_equal, assert_raises
 
+import nifty5 as ift
+
+pmp = pytest.mark.parametrize
 SPACES = [ift.RGSpace((4,)), ift.RGSpace((5))]
 SPACE_COMBINATIONS = [(), SPACES[0], SPACES[1], SPACES]
 
 
-class Test_Interface(unittest.TestCase):
-    @expand(product(SPACE_COMBINATIONS,
-                    [['domain', ift.DomainTuple],
-                     ['val', ift.dobj.data_object],
-                     ['shape', tuple],
-                     ['size', (np.int, np.int64)]]))
-    def test_return_types(self, domain, attribute_desired_type):
-        attribute = attribute_desired_type[0]
-        desired_type = attribute_desired_type[1]
-        f = ift.Field.full(domain, 1.)
-        assert_equal(isinstance(getattr(f, attribute), desired_type), True)
+@pmp('domain', SPACE_COMBINATIONS)
+@pmp('attribute_desired_type',
+     [['domain', ift.DomainTuple], ['val', ift.dobj.data_object],
+      ['shape', tuple], ['size', (np.int, np.int64)]])
+def test_return_types(domain, attribute_desired_type):
+    attribute = attribute_desired_type[0]
+    desired_type = attribute_desired_type[1]
+    f = ift.Field.full(domain, 1.)
+    assert_equal(isinstance(getattr(f, attribute), desired_type), True)
 
 
 def _spec1(k):
-    return 42/(1.+k)**2
+    return 42/(1. + k)**2
 
 
 def _spec2(k):
-    return 42/(1.+k)**3
-
-
-class Test_Functionality(unittest.TestCase):
-    @expand(product([ift.RGSpace((8,), harmonic=True),
-                     ift.RGSpace((8, 8), harmonic=True, distances=0.123)],
-                    [ift.RGSpace((8,), harmonic=True),
-                     ift.LMSpace(12)]))
-    def test_power_synthesize_analyze(self, space1, space2):
-        np.random.seed(11)
-
-        p1 = ift.PowerSpace(space1)
-        fp1 = ift.PS_field(p1, _spec1)
-        p2 = ift.PowerSpace(space2)
-        fp2 = ift.PS_field(p2, _spec2)
-        outer = np.outer(fp1.to_global_data(), fp2.to_global_data())
-        fp = ift.Field.from_global_data((p1, p2), outer)
-
-        op1 = ift.create_power_operator((space1, space2), _spec1, 0)
-        op2 = ift.create_power_operator((space1, space2), _spec2, 1)
-        opfull = op2(op1)
-
-        samples = 500
-        sc1 = ift.StatCalculator()
-        sc2 = ift.StatCalculator()
-        for ii in range(samples):
-            sk = opfull.draw_sample()
-
-            sp = ift.power_analyze(sk, spaces=(0, 1),
-                                   keep_phase_information=False)
-            sc1.add(sp.sum(spaces=1)/fp2.sum())
-            sc2.add(sp.sum(spaces=0)/fp1.sum())
-
-        assert_allclose(sc1.mean.local_data, fp1.local_data, rtol=0.2)
-        assert_allclose(sc2.mean.local_data, fp2.local_data, rtol=0.2)
-
-    @expand(product([ift.RGSpace((8,), harmonic=True),
-                     ift.RGSpace((8, 8), harmonic=True, distances=0.123)],
-                    [ift.RGSpace((8,), harmonic=True),
-                     ift.LMSpace(12)]))
-    def test_DiagonalOperator_power_analyze2(self, space1, space2):
-        np.random.seed(11)
-
-        fp1 = ift.PS_field(ift.PowerSpace(space1), _spec1)
-        fp2 = ift.PS_field(ift.PowerSpace(space2), _spec2)
-
-        S_1 = ift.create_power_operator((space1, space2), _spec1, 0)
-        S_2 = ift.create_power_operator((space1, space2), _spec2, 1)
-        S_full = S_2(S_1)
-
-        samples = 500
-        sc1 = ift.StatCalculator()
-        sc2 = ift.StatCalculator()
-
-        for ii in range(samples):
-            sk = S_full.draw_sample()
-            sp = ift.power_analyze(sk, spaces=(0, 1),
-                                   keep_phase_information=False)
-            sc1.add(sp.sum(spaces=1)/fp2.sum())
-            sc2.add(sp.sum(spaces=0)/fp1.sum())
-
-        assert_allclose(sc1.mean.local_data, fp1.local_data, rtol=0.2)
-        assert_allclose(sc2.mean.local_data, fp2.local_data, rtol=0.2)
-
-    @expand(product([ift.RGSpace((8,), harmonic=True), (),
-                     ift.RGSpace((8, 8), harmonic=True, distances=0.123),
-                     ift.RGSpace((2, 3, 7))]))
-    def test_norm(self, space):
-        f = ift.Field.from_random("normal", domain=space, dtype=np.complex128)
-        gd = f.to_global_data().reshape(-1)
-        assert_allclose(f.norm(), np.linalg.norm(gd))
-        assert_allclose(f.norm(1), np.linalg.norm(gd, ord=1))
-        assert_allclose(f.norm(2), np.linalg.norm(gd, ord=2))
-        assert_allclose(f.norm(3), np.linalg.norm(gd, ord=3))
-        assert_allclose(f.norm(np.inf), np.linalg.norm(gd, ord=np.inf))
-
-    def test_vdot(self):
-        s = ift.RGSpace((10,))
-        f1 = ift.Field.from_random("normal", domain=s, dtype=np.complex128)
-        f2 = ift.Field.from_random("normal", domain=s, dtype=np.complex128)
-        assert_allclose(f1.vdot(f2), f1.vdot(f2, spaces=0))
-        assert_allclose(f1.vdot(f2), np.conj(f2.vdot(f1)))
-
-    def test_vdot2(self):
-        x1 = ift.RGSpace((200,))
-        x2 = ift.RGSpace((150,))
-        m = ift.Field.full((x1, x2), .5)
-        res = m.vdot(m, spaces=1)
-        assert_allclose(res.local_data, 37.5)
-
-    def test_outer(self):
-        x1 = ift.RGSpace((9,))
-        x2 = ift.RGSpace((3,))
-        m1 = ift.Field.full(x1, .5)
-        m2 = ift.Field.full(x2, 3.)
-        res = m1.outer(m2)
-        assert_allclose(res.to_global_data(), np.full((9, 3,), 1.5))
-
-    def test_sum(self):
-        x1 = ift.RGSpace((9,), distances=2.)
-        x2 = ift.RGSpace((2, 12,), distances=(0.3,))
-        m1 = ift.Field.from_global_data(ift.makeDomain(x1), np.arange(9))
-        m2 = ift.Field.full(ift.makeDomain((x1, x2)), 0.45)
-        res1 = m1.sum()
-        res2 = m2.sum(spaces=1)
-        assert_allclose(res1, 36)
-        assert_allclose(res2.to_global_data(), np.full(9, 2*12*0.45))
-
-    def test_integrate(self):
-        x1 = ift.RGSpace((9,), distances=2.)
-        x2 = ift.RGSpace((2, 12,), distances=(0.3,))
-        m1 = ift.Field.from_global_data(ift.makeDomain(x1), np.arange(9))
-        m2 = ift.Field.full(ift.makeDomain((x1, x2)), 0.45)
-        res1 = m1.integrate()
-        res2 = m2.integrate(spaces=1)
-        assert_allclose(res1, 36*2)
-        assert_allclose(res2.to_global_data(), np.full(9, 2*12*0.45*0.3**2))
-
-    def test_dataconv(self):
-        s1 = ift.RGSpace((10,))
-        ld = np.arange(ift.dobj.local_shape(s1.shape)[0])
-        gd = np.arange(s1.shape[0])
-        assert_equal(ld, ift.from_local_data(s1, ld).local_data)
-        assert_equal(gd, ift.from_global_data(s1, gd).to_global_data())
-
-    def test_cast_domain(self):
-        s1 = ift.RGSpace((10,))
-        s2 = ift.RGSpace((10,), distances=20.)
-        d = np.arange(s1.shape[0])
-        d2 = ift.from_global_data(s1, d).cast_domain(s2).to_global_data()
-        assert_equal(d, d2)
-
-    def test_empty_domain(self):
-        f = ift.Field.full((), 5)
-        assert_equal(f.local_data, 5)
-        f = ift.Field.full(None, 5)
-        assert_equal(f.local_data, 5)
-
-    def test_trivialities(self):
-        s1 = ift.RGSpace((10,))
-        f1 = ift.Field.full(s1, 27)
-        assert_equal(f1.local_data, f1.real.local_data)
-        f1 = ift.Field.full(s1, 27.+3j)
-        assert_equal(f1.real.local_data, 27.)
-        assert_equal(f1.imag.local_data, 3.)
-        assert_equal(f1.local_data, +f1.local_data)
-        assert_equal(f1.sum(), f1.sum(0))
-        f1 = ift.from_global_data(s1, np.arange(10))
-#        assert_equal(f1.min(), 0)
-#        assert_equal(f1.max(), 9)
-        assert_equal(f1.prod(), 0)
-
-    def test_weight(self):
-        s1 = ift.RGSpace((10,))
-        f = ift.Field.full(s1, 10.)
-        f2 = f.weight(1)
-        assert_equal(f.weight(1).local_data, f2.local_data)
-        assert_equal(f.total_volume(), 1)
-        assert_equal(f.total_volume(0), 1)
-        assert_equal(f.total_volume((0,)), 1)
-        assert_equal(f.scalar_weight(), 0.1)
-        assert_equal(f.scalar_weight(0), 0.1)
-        assert_equal(f.scalar_weight((0,)), 0.1)
-        s1 = ift.GLSpace(10)
-        f = ift.Field.full(s1, 10.)
-        assert_equal(f.scalar_weight(), None)
-        assert_equal(f.scalar_weight(0), None)
-        assert_equal(f.scalar_weight((0,)), None)
-
-    @expand(product([ift.RGSpace(10), ift.GLSpace(10)],
-                    [np.float64, np.complex128]))
-    def test_reduction(self, dom, dt):
-        s1 = ift.Field.full(dom, dt(1.))
-        assert_allclose(s1.mean(), 1.)
-        assert_allclose(s1.mean(0), 1.)
-        assert_allclose(s1.var(), 0., atol=1e-14)
-        assert_allclose(s1.var(0), 0., atol=1e-14)
-        assert_allclose(s1.std(), 0., atol=1e-14)
-        assert_allclose(s1.std(0), 0., atol=1e-14)
-
-    def test_err(self):
-        s1 = ift.RGSpace((10,))
-        s2 = ift.RGSpace((11,))
-        f1 = ift.Field.full(s1, 27)
-        with assert_raises(ValueError):
-            f2 = ift.Field(ift.DomainTuple.make(s2), f1.val)
-        with assert_raises(TypeError):
-            f2 = ift.Field.full(s2, "xyz")
-        with assert_raises(TypeError):
-            if f1:
-                pass
-        with assert_raises(TypeError):
-            f1.full((2, 4, 6))
-        with assert_raises(TypeError):
-            f2 = ift.Field(None, None)
-        with assert_raises(TypeError):
-            f2 = ift.Field(s1, None)
-        with assert_raises(ValueError):
-            f1.imag
-        with assert_raises(TypeError):
-            f1.vdot(42)
-        with assert_raises(ValueError):
-            f1.vdot(ift.Field.full(s2, 1.))
-        with assert_raises(TypeError):
-            ift.full(s1, [2, 3])
-
-    def test_stdfunc(self):
-        s = ift.RGSpace((200,))
-        f = ift.Field.full(s, 27)
-        assert_equal(f.local_data, 27)
-        assert_equal(f.shape, (200,))
-        assert_equal(f.dtype, np.int)
-        fx = ift.full(f.domain, 0)
-        assert_equal(f.dtype, fx.dtype)
-        assert_equal(f.shape, fx.shape)
-        assert_equal(fx.local_data, 0)
-        fx = ift.full(f.domain, 1)
-        assert_equal(f.dtype, fx.dtype)
-        assert_equal(f.shape, fx.shape)
-        assert_equal(fx.local_data, 1)
-        fx = ift.full(f.domain, 67.)
-        assert_equal(f.shape, fx.shape)
-        assert_equal(fx.local_data, 67.)
-        f = ift.Field.from_random("normal", s)
-        f2 = ift.Field.from_random("normal", s)
-        assert_equal((f > f2).local_data, f.local_data > f2.local_data)
-        assert_equal((f >= f2).local_data, f.local_data >= f2.local_data)
-        assert_equal((f < f2).local_data, f.local_data < f2.local_data)
-        assert_equal((f <= f2).local_data, f.local_data <= f2.local_data)
-        assert_equal((f != f2).local_data, f.local_data != f2.local_data)
-        assert_equal((f == f2).local_data, f.local_data == f2.local_data)
-        assert_equal((f+f2).local_data, f.local_data+f2.local_data)
-        assert_equal((f-f2).local_data, f.local_data-f2.local_data)
-        assert_equal((f*f2).local_data, f.local_data*f2.local_data)
-        assert_equal((f/f2).local_data, f.local_data/f2.local_data)
-        assert_equal((-f).local_data, -(f.local_data))
-        assert_equal(abs(f).local_data, abs(f.local_data))
-
-    def test_emptydomain(self):
-        f = ift.Field.full((), 3.)
-        assert_equal(f.sum(), 3.)
-        assert_equal(f.prod(), 3.)
-        assert_equal(f.local_data, 3.)
-        assert_equal(f.local_data.shape, ())
-        assert_equal(f.local_data.size, 1)
-        assert_equal(f.vdot(f), 9.)
-
-    @expand(product([float(5), 5.],
-                    [ift.RGSpace((8,), harmonic=True), ()],
-                    ["exp", "log", "sin", "cos", "tan", "sinh", "cosh", "sinc",
-                     "absolute", "sign"]))
-    def test_funcs(self, num, dom, func):
-        num = 5
-        f = ift.Field.full(dom, num)
-        res = getattr(f, func)()
-        res2 = getattr(np, func)(num)
-        assert_allclose(res.local_data, res2)
+    return 42/(1. + k)**3
+
+
+@pmp('space1', [
+    ift.RGSpace((8,), harmonic=True),
+    ift.RGSpace((8, 8), harmonic=True, distances=0.123)
+])
+@pmp('space2', [ift.RGSpace((8,), harmonic=True), ift.LMSpace(12)])
+def test_power_synthesize_analyze(space1, space2):
+    np.random.seed(11)
+
+    p1 = ift.PowerSpace(space1)
+    fp1 = ift.PS_field(p1, _spec1)
+    p2 = ift.PowerSpace(space2)
+    fp2 = ift.PS_field(p2, _spec2)
+    outer = np.outer(fp1.to_global_data(), fp2.to_global_data())
+    fp = ift.Field.from_global_data((p1, p2), outer)
+
+    op1 = ift.create_power_operator((space1, space2), _spec1, 0)
+    op2 = ift.create_power_operator((space1, space2), _spec2, 1)
+    opfull = op2(op1)
+
+    samples = 500
+    sc1 = ift.StatCalculator()
+    sc2 = ift.StatCalculator()
+    for ii in range(samples):
+        sk = opfull.draw_sample()
+
+        sp = ift.power_analyze(sk, spaces=(0, 1), keep_phase_information=False)
+        sc1.add(sp.sum(spaces=1)/fp2.sum())
+        sc2.add(sp.sum(spaces=0)/fp1.sum())
+
+    assert_allclose(sc1.mean.local_data, fp1.local_data, rtol=0.2)
+    assert_allclose(sc2.mean.local_data, fp2.local_data, rtol=0.2)
+
+
+@pmp('space1', [
+    ift.RGSpace((8,), harmonic=True),
+    ift.RGSpace((8, 8), harmonic=True, distances=0.123)
+])
+@pmp('space2', [ift.RGSpace((8,), harmonic=True), ift.LMSpace(12)])
+def test_DiagonalOperator_power_analyze2(space1, space2):
+    np.random.seed(11)
+
+    fp1 = ift.PS_field(ift.PowerSpace(space1), _spec1)
+    fp2 = ift.PS_field(ift.PowerSpace(space2), _spec2)
+
+    S_1 = ift.create_power_operator((space1, space2), _spec1, 0)
+    S_2 = ift.create_power_operator((space1, space2), _spec2, 1)
+    S_full = S_2(S_1)
+
+    samples = 500
+    sc1 = ift.StatCalculator()
+    sc2 = ift.StatCalculator()
+
+    for ii in range(samples):
+        sk = S_full.draw_sample()
+        sp = ift.power_analyze(sk, spaces=(0, 1), keep_phase_information=False)
+        sc1.add(sp.sum(spaces=1)/fp2.sum())
+        sc2.add(sp.sum(spaces=0)/fp1.sum())
+
+    assert_allclose(sc1.mean.local_data, fp1.local_data, rtol=0.2)
+    assert_allclose(sc2.mean.local_data, fp2.local_data, rtol=0.2)
+
+
+@pmp('space', [
+    ift.RGSpace((8,), harmonic=True), (),
+    ift.RGSpace((8, 8), harmonic=True, distances=0.123),
+    ift.RGSpace((2, 3, 7))
+])
+def test_norm(space):
+    f = ift.Field.from_random("normal", domain=space, dtype=np.complex128)
+    gd = f.to_global_data().reshape(-1)
+    assert_allclose(f.norm(), np.linalg.norm(gd))
+    assert_allclose(f.norm(1), np.linalg.norm(gd, ord=1))
+    assert_allclose(f.norm(2), np.linalg.norm(gd, ord=2))
+    assert_allclose(f.norm(3), np.linalg.norm(gd, ord=3))
+    assert_allclose(f.norm(np.inf), np.linalg.norm(gd, ord=np.inf))
+
+
+def test_vdot():
+    s = ift.RGSpace((10,))
+    f1 = ift.Field.from_random("normal", domain=s, dtype=np.complex128)
+    f2 = ift.Field.from_random("normal", domain=s, dtype=np.complex128)
+    assert_allclose(f1.vdot(f2), f1.vdot(f2, spaces=0))
+    assert_allclose(f1.vdot(f2), np.conj(f2.vdot(f1)))
+
+
+def test_vdot2():
+    x1 = ift.RGSpace((200,))
+    x2 = ift.RGSpace((150,))
+    m = ift.Field.full((x1, x2), .5)
+    res = m.vdot(m, spaces=1)
+    assert_allclose(res.local_data, 37.5)
+
+
+def test_outer():
+    x1 = ift.RGSpace((9,))
+    x2 = ift.RGSpace((3,))
+    m1 = ift.Field.full(x1, .5)
+    m2 = ift.Field.full(x2, 3.)
+    res = m1.outer(m2)
+    assert_allclose(res.to_global_data(), np.full((9, 3), 1.5))
+
+
+def test_sum():
+    x1 = ift.RGSpace((9,), distances=2.)
+    x2 = ift.RGSpace(
+        (
+            2,
+            12,
+        ), distances=(0.3,))
+    m1 = ift.Field.from_global_data(ift.makeDomain(x1), np.arange(9))
+    m2 = ift.Field.full(ift.makeDomain((x1, x2)), 0.45)
+    res1 = m1.sum()
+    res2 = m2.sum(spaces=1)
+    assert_allclose(res1, 36)
+    assert_allclose(res2.to_global_data(), np.full(9, 2*12*0.45))
+
+
+def test_integrate():
+    x1 = ift.RGSpace((9,), distances=2.)
+    x2 = ift.RGSpace((2, 12), distances=(0.3,))
+    m1 = ift.Field.from_global_data(ift.makeDomain(x1), np.arange(9))
+    m2 = ift.Field.full(ift.makeDomain((x1, x2)), 0.45)
+    res1 = m1.integrate()
+    res2 = m2.integrate(spaces=1)
+    assert_allclose(res1, 36*2)
+    assert_allclose(res2.to_global_data(), np.full(9, 2*12*0.45*0.3**2))
+
+
+def test_dataconv():
+    s1 = ift.RGSpace((10,))
+    ld = np.arange(ift.dobj.local_shape(s1.shape)[0])
+    gd = np.arange(s1.shape[0])
+    assert_equal(ld, ift.from_local_data(s1, ld).local_data)
+    assert_equal(gd, ift.from_global_data(s1, gd).to_global_data())
+
+
+def test_cast_domain():
+    s1 = ift.RGSpace((10,))
+    s2 = ift.RGSpace((10,), distances=20.)
+    d = np.arange(s1.shape[0])
+    d2 = ift.from_global_data(s1, d).cast_domain(s2).to_global_data()
+    assert_equal(d, d2)
+
+
+def test_empty_domain():
+    f = ift.Field.full((), 5)
+    assert_equal(f.local_data, 5)
+    f = ift.Field.full(None, 5)
+    assert_equal(f.local_data, 5)
+
+
+def test_trivialities():
+    s1 = ift.RGSpace((10,))
+    f1 = ift.Field.full(s1, 27)
+    assert_equal(f1.local_data, f1.real.local_data)
+    f1 = ift.Field.full(s1, 27. + 3j)
+    assert_equal(f1.real.local_data, 27.)
+    assert_equal(f1.imag.local_data, 3.)
+    assert_equal(f1.local_data, +f1.local_data)
+    assert_equal(f1.sum(), f1.sum(0))
+    f1 = ift.from_global_data(s1, np.arange(10))
+    # assert_equal(f1.min(), 0)
+    # assert_equal(f1.max(), 9)
+    assert_equal(f1.prod(), 0)
+
+
+def test_weight():
+    s1 = ift.RGSpace((10,))
+    f = ift.Field.full(s1, 10.)
+    f2 = f.weight(1)
+    assert_equal(f.weight(1).local_data, f2.local_data)
+    assert_equal(f.total_volume(), 1)
+    assert_equal(f.total_volume(0), 1)
+    assert_equal(f.total_volume((0,)), 1)
+    assert_equal(f.scalar_weight(), 0.1)
+    assert_equal(f.scalar_weight(0), 0.1)
+    assert_equal(f.scalar_weight((0,)), 0.1)
+    s1 = ift.GLSpace(10)
+    f = ift.Field.full(s1, 10.)
+    assert_equal(f.scalar_weight(), None)
+    assert_equal(f.scalar_weight(0), None)
+    assert_equal(f.scalar_weight((0,)), None)
+
+
+@pmp('dom', [ift.RGSpace(10), ift.GLSpace(10)])
+@pmp('dt', [np.float64, np.complex128])
+def test_reduction(dom, dt):
+    s1 = ift.Field.full(dom, dt(1.))
+    assert_allclose(s1.mean(), 1.)
+    assert_allclose(s1.mean(0), 1.)
+    assert_allclose(s1.var(), 0., atol=1e-14)
+    assert_allclose(s1.var(0), 0., atol=1e-14)
+    assert_allclose(s1.std(), 0., atol=1e-14)
+    assert_allclose(s1.std(0), 0., atol=1e-14)
+
+
+def test_err():
+    s1 = ift.RGSpace((10,))
+    s2 = ift.RGSpace((11,))
+    f1 = ift.Field.full(s1, 27)
+    with assert_raises(ValueError):
+        f2 = ift.Field(ift.DomainTuple.make(s2), f1.val)
+    with assert_raises(TypeError):
+        f2 = ift.Field.full(s2, "xyz")
+    with assert_raises(TypeError):
+        if f1:
+            pass
+    with assert_raises(TypeError):
+        f1.full((2, 4, 6))
+    with assert_raises(TypeError):
+        f2 = ift.Field(None, None)
+    with assert_raises(TypeError):
+        f2 = ift.Field(s1, None)
+    with assert_raises(ValueError):
+        f1.imag
+    with assert_raises(TypeError):
+        f1.vdot(42)
+    with assert_raises(ValueError):
+        f1.vdot(ift.Field.full(s2, 1.))
+    with assert_raises(TypeError):
+        ift.full(s1, [2, 3])
+
+
+def test_stdfunc():
+    s = ift.RGSpace((200,))
+    f = ift.Field.full(s, 27)
+    assert_equal(f.local_data, 27)
+    assert_equal(f.shape, (200,))
+    assert_equal(f.dtype, np.int)
+    fx = ift.full(f.domain, 0)
+    assert_equal(f.dtype, fx.dtype)
+    assert_equal(f.shape, fx.shape)
+    assert_equal(fx.local_data, 0)
+    fx = ift.full(f.domain, 1)
+    assert_equal(f.dtype, fx.dtype)
+    assert_equal(f.shape, fx.shape)
+    assert_equal(fx.local_data, 1)
+    fx = ift.full(f.domain, 67.)
+    assert_equal(f.shape, fx.shape)
+    assert_equal(fx.local_data, 67.)
+    f = ift.Field.from_random("normal", s)
+    f2 = ift.Field.from_random("normal", s)
+    assert_equal((f > f2).local_data, f.local_data > f2.local_data)
+    assert_equal((f >= f2).local_data, f.local_data >= f2.local_data)
+    assert_equal((f < f2).local_data, f.local_data < f2.local_data)
+    assert_equal((f <= f2).local_data, f.local_data <= f2.local_data)
+    assert_equal((f != f2).local_data, f.local_data != f2.local_data)
+    assert_equal((f == f2).local_data, f.local_data == f2.local_data)
+    assert_equal((f + f2).local_data, f.local_data + f2.local_data)
+    assert_equal((f - f2).local_data, f.local_data - f2.local_data)
+    assert_equal((f*f2).local_data, f.local_data*f2.local_data)
+    assert_equal((f/f2).local_data, f.local_data/f2.local_data)
+    assert_equal((-f).local_data, -(f.local_data))
+    assert_equal(abs(f).local_data, abs(f.local_data))
+
+
+def test_emptydomain():
+    f = ift.Field.full((), 3.)
+    assert_equal(f.sum(), 3.)
+    assert_equal(f.prod(), 3.)
+    assert_equal(f.local_data, 3.)
+    assert_equal(f.local_data.shape, ())
+    assert_equal(f.local_data.size, 1)
+    assert_equal(f.vdot(f), 9.)
+
+
+@pmp('num', [float(5), 5.])
+@pmp('dom', [ift.RGSpace((8,), harmonic=True), ()])
+@pmp('func', [
+    "exp", "log", "sin", "cos", "tan", "sinh", "cosh", "sinc", "absolute",
+    "sign"
+])
+def test_funcs(num, dom, func):
+    num = 5
+    f = ift.Field.full(dom, num)
+    res = getattr(f, func)()
+    res2 = getattr(np, func)(num)
+    assert_allclose(res.local_data, res2)
diff --git a/test/test_gaussian_metric.py b/test/test_gaussian_metric.py
new file mode 100644
index 0000000000000000000000000000000000000000..aa9f4beb76ef8ee90fb8dde44daacf268202cdb5
--- /dev/null
+++ b/test/test_gaussian_metric.py
@@ -0,0 +1,81 @@
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+import pytest
+
+import nifty5 as ift
+
+
+def _flat_PS(k):
+    return np.ones_like(k)
+
+
+pmp = pytest.mark.parametrize
+
+
+@pmp('space', [
+    ift.GLSpace(15),
+    ift.RGSpace(64, distances=.789),
+    ift.RGSpace([32, 32], distances=.789)
+])
+@pmp('nonlinearity', ["tanh", "exp", ""])
+@pmp('noise', [1, 1e-2, 1e2])
+@pmp('seed', [4, 78, 23])
+def test_gaussian_energy(space, nonlinearity, noise, seed):
+    np.random.seed(seed)
+    dim = len(space.shape)
+    hspace = space.get_default_codomain()
+    ht = ift.HarmonicTransformOperator(hspace, target=space)
+    binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
+    pspace = ift.PowerSpace(hspace, binbounds=binbounds)
+    Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
+    xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
+    # FIXME Needed?
+    xi0_var = ift.Linearization.make_var(xi0)
+
+    def pspec(k):
+        return 1/(1 + k**2)**dim
+
+    pspec = ift.PS_field(pspace, pspec)
+    A = Dist(ift.sqrt(pspec))
+    N = ift.ScalingOperator(noise, space)
+    n = N.draw_sample()
+    # FIXME Needed?
+    s = ht(ift.makeOp(A)(xi0_var))
+    R = ift.ScalingOperator(10., space)
+
+    def d_model():
+        if nonlinearity == "":
+            return R(ht(ift.makeOp(A)))
+        else:
+            tmp = ht(ift.makeOp(A))
+            nonlin = getattr(tmp, nonlinearity)()
+            return R(nonlin)
+
+    d = d_model()(xi0) + n
+
+    if noise == 1:
+        N = None
+
+    energy = ift.GaussianEnergy(d, N)(d_model())
+    if nonlinearity == "":
+        ift.extra.check_value_gradient_metric_consistency(
+            energy, xi0, ntries=10)
+    else:
+        ift.extra.check_value_gradient_consistency(
+            energy, xi0, ntries=10, tol=5e-8)
diff --git a/test/test_minimization/test_minimizers.py b/test/test_minimization/test_minimizers.py
deleted file mode 100644
index 364f11ed3bf1e3df280f05af0f4e41e2cd8e851f..0000000000000000000000000000000000000000
--- a/test/test_minimization/test_minimizers.py
+++ /dev/null
@@ -1,212 +0,0 @@
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-# Copyright(C) 2013-2019 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
-import numpy as np
-from unittest import SkipTest
-from numpy.testing import assert_allclose, assert_equal
-
-IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000)
-
-spaces = [ift.RGSpace([1024], distances=0.123), ift.HPSpace(32)]
-
-minimizers = ['ift.VL_BFGS(IC)',
-              'ift.NonlinearCG(IC, "Polak-Ribiere")',
-              # 'ift.NonlinearCG(IC, "Hestenes-Stiefel"),
-              'ift.NonlinearCG(IC, "Fletcher-Reeves")',
-              'ift.NonlinearCG(IC, "5.49")',
-              'ift.L_BFGS_B(ftol=1e-10, gtol=1e-5, maxiter=1000)',
-              'ift.L_BFGS(IC)',
-              'ift.NewtonCG(IC)']
-
-newton_minimizers = ['ift.RelaxedNewton(IC)']
-quadratic_only_minimizers = ['ift.ConjugateGradient(IC)',
-                             'ift.ScipyCG(tol=1e-5, maxiter=300)']
-slow_minimizers = ['ift.SteepestDescent(IC)']
-
-
-class Test_Minimizers(unittest.TestCase):
-
-    @expand(product(minimizers + newton_minimizers +
-                    quadratic_only_minimizers + slow_minimizers, spaces))
-    def test_quadratic_minimization(self, minimizer, space):
-        np.random.seed(42)
-        starting_point = ift.Field.from_random('normal', domain=space)*10
-        covariance_diagonal = ift.Field.from_random(
-                                  'uniform', domain=space) + 0.5
-        covariance = ift.DiagonalOperator(covariance_diagonal)
-        required_result = ift.full(space, 1.)
-
-        try:
-            minimizer = eval(minimizer)
-            energy = ift.QuadraticEnergy(A=covariance, b=required_result,
-                                         position=starting_point)
-
-            (energy, convergence) = minimizer(energy)
-        except NotImplementedError:
-            raise SkipTest
-
-        assert_equal(convergence, IC.CONVERGED)
-        assert_allclose(energy.position.local_data,
-                        1./covariance_diagonal.local_data,
-                        rtol=1e-3, atol=1e-3)
-
-    @expand(product(minimizers+newton_minimizers))
-    def test_rosenbrock(self, minimizer):
-        try:
-            from scipy.optimize import rosen, rosen_der, rosen_hess_prod
-        except ImportError:
-            raise SkipTest
-        np.random.seed(42)
-        space = ift.DomainTuple.make(ift.UnstructuredDomain((2,)))
-        starting_point = ift.Field.from_random('normal', domain=space)*10
-
-        class RBEnergy(ift.Energy):
-            def __init__(self, position):
-                super(RBEnergy, self).__init__(position)
-
-            @property
-            def value(self):
-                return rosen(self._position.to_global_data_rw())
-
-            @property
-            def gradient(self):
-                inp = self._position.to_global_data_rw()
-                out = ift.Field.from_global_data(space, rosen_der(inp))
-                return out
-
-            @property
-            def metric(self):
-                class RBCurv(ift.EndomorphicOperator):
-                    def __init__(self, loc):
-                        self._loc = loc.to_global_data_rw()
-                        self._capability = self.TIMES
-                        self._domain = space
-
-                    def apply(self, x, mode):
-                        self._check_input(x, mode)
-                        inp = x.to_global_data_rw()
-                        out = ift.Field.from_global_data(
-                            space, rosen_hess_prod(self._loc.copy(), inp))
-                        return out
-
-                t1 = ift.GradientNormController(tol_abs_gradnorm=1e-5,
-                                                iteration_limit=1000)
-                return ift.InversionEnabler(RBCurv(self._position), t1)
-
-            def apply_metric(self, x):
-                inp = x.to_global_data_rw()
-                pos = self._position.to_global_data_rw()
-                return ift.Field.from_global_data(
-                    space, rosen_hess_prod(pos, inp))
-
-        try:
-            minimizer = eval(minimizer)
-            energy = RBEnergy(position=starting_point)
-
-            (energy, convergence) = minimizer(energy)
-        except NotImplementedError:
-            raise SkipTest
-
-        assert_equal(convergence, IC.CONVERGED)
-        assert_allclose(energy.position.local_data, 1.,
-                        rtol=1e-3, atol=1e-3)
-
-    @expand(product(minimizers+slow_minimizers))
-    def test_gauss(self, minimizer):
-        space = ift.UnstructuredDomain((1,))
-        starting_point = ift.Field.full(space, 3.)
-
-        class ExpEnergy(ift.Energy):
-            def __init__(self, position):
-                super(ExpEnergy, self).__init__(position)
-
-            @property
-            def value(self):
-                x = self.position.to_global_data()[0]
-                return -np.exp(-(x**2))
-
-            @property
-            def gradient(self):
-                x = self.position.to_global_data()[0]
-                return ift.Field.full(self.position.domain,
-                                      2*x*np.exp(-(x**2)))
-
-            def apply_metric(self, x):
-                p = self.position.to_global_data()[0]
-                v = (2 - 4*p*p)*np.exp(-p**2)
-                return ift.DiagonalOperator(
-                    ift.Field.full(self.position.domain, v))(x)
-
-        try:
-            minimizer = eval(minimizer)
-            energy = ExpEnergy(position=starting_point)
-
-            (energy, convergence) = minimizer(energy)
-        except NotImplementedError:
-            raise SkipTest
-
-        assert_equal(convergence, IC.CONVERGED)
-        assert_allclose(energy.position.local_data, 0.,
-                        atol=1e-3)
-
-    @expand(product(minimizers+newton_minimizers+slow_minimizers))
-    def test_cosh(self, minimizer):
-        space = ift.UnstructuredDomain((1,))
-        starting_point = ift.Field.full(space, 3.)
-
-        class CoshEnergy(ift.Energy):
-            def __init__(self, position):
-                super(CoshEnergy, self).__init__(position)
-
-            @property
-            def value(self):
-                x = self.position.to_global_data()[0]
-                return np.cosh(x)
-
-            @property
-            def gradient(self):
-                x = self.position.to_global_data()[0]
-                return ift.Field.full(self.position.domain, np.sinh(x))
-
-            @property
-            def metric(self):
-                x = self.position.to_global_data()[0]
-                v = np.cosh(x)
-                return ift.DiagonalOperator(
-                    ift.Field.full(self.position.domain, v))
-
-            def apply_metric(self, x):
-                p = self.position.to_global_data()[0]
-                v = np.cosh(p)
-                return ift.DiagonalOperator(
-                    ift.Field.full(self.position.domain, v))(x)
-
-        try:
-            minimizer = eval(minimizer)
-            energy = CoshEnergy(position=starting_point)
-
-            (energy, convergence) = minimizer(energy)
-        except NotImplementedError:
-            raise SkipTest
-
-        assert_equal(convergence, IC.CONVERGED)
-        assert_allclose(energy.position.local_data, 0., atol=1e-3)
diff --git a/test/test_minimizers.py b/test/test_minimizers.py
new file mode 100644
index 0000000000000000000000000000000000000000..9ab55cb9ae4eb7ed23ec3bdf3f306e0797e87662
--- /dev/null
+++ b/test/test_minimizers.py
@@ -0,0 +1,214 @@
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+from unittest import SkipTest
+
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose, assert_equal
+
+import nifty5 as ift
+
+pmp = pytest.mark.parametrize
+IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000)
+
+spaces = [ift.RGSpace([1024], distances=0.123), ift.HPSpace(32)]
+
+minimizers = [
+    'ift.VL_BFGS(IC)',
+    'ift.NonlinearCG(IC, "Polak-Ribiere")',
+    # 'ift.NonlinearCG(IC, "Hestenes-Stiefel"),
+    'ift.NonlinearCG(IC, "Fletcher-Reeves")',
+    'ift.NonlinearCG(IC, "5.49")',
+    'ift.L_BFGS_B(ftol=1e-10, gtol=1e-5, maxiter=1000)',
+    'ift.L_BFGS(IC)',
+    'ift.NewtonCG(IC)'
+]
+
+newton_minimizers = ['ift.RelaxedNewton(IC)']
+quadratic_only_minimizers = [
+    'ift.ConjugateGradient(IC)', 'ift.ScipyCG(tol=1e-5, maxiter=300)'
+]
+slow_minimizers = ['ift.SteepestDescent(IC)']
+
+
+@pmp('minimizer', minimizers + newton_minimizers + quadratic_only_minimizers +
+     slow_minimizers)
+@pmp('space', spaces)
+def test_quadratic_minimization(minimizer, space):
+    np.random.seed(42)
+    starting_point = ift.Field.from_random('normal', domain=space)*10
+    covariance_diagonal = ift.Field.from_random('uniform', domain=space) + 0.5
+    covariance = ift.DiagonalOperator(covariance_diagonal)
+    required_result = ift.full(space, 1.)
+
+    try:
+        minimizer = eval(minimizer)
+        energy = ift.QuadraticEnergy(
+            A=covariance, b=required_result, position=starting_point)
+
+        (energy, convergence) = minimizer(energy)
+    except NotImplementedError:
+        raise SkipTest
+
+    assert_equal(convergence, IC.CONVERGED)
+    assert_allclose(
+        energy.position.local_data,
+        1./covariance_diagonal.local_data,
+        rtol=1e-3,
+        atol=1e-3)
+
+
+@pmp('minimizer', minimizers + newton_minimizers)
+def test_rosenbrock(minimizer):
+    try:
+        from scipy.optimize import rosen, rosen_der, rosen_hess_prod
+    except ImportError:
+        raise SkipTest
+    np.random.seed(42)
+    space = ift.DomainTuple.make(ift.UnstructuredDomain((2,)))
+    starting_point = ift.Field.from_random('normal', domain=space)*10
+
+    class RBEnergy(ift.Energy):
+        def __init__(self, position):
+            super(RBEnergy, self).__init__(position)
+
+        @property
+        def value(self):
+            return rosen(self._position.to_global_data_rw())
+
+        @property
+        def gradient(self):
+            inp = self._position.to_global_data_rw()
+            out = ift.Field.from_global_data(space, rosen_der(inp))
+            return out
+
+        @property
+        def metric(self):
+            class RBCurv(ift.EndomorphicOperator):
+                def __init__(self, loc):
+                    self._loc = loc.to_global_data_rw()
+                    self._capability = self.TIMES
+                    self._domain = space
+
+                def apply(self, x, mode):
+                    self._check_input(x, mode)
+                    inp = x.to_global_data_rw()
+                    out = ift.Field.from_global_data(
+                        space, rosen_hess_prod(self._loc.copy(), inp))
+                    return out
+
+            t1 = ift.GradientNormController(
+                tol_abs_gradnorm=1e-5, iteration_limit=1000)
+            return ift.InversionEnabler(RBCurv(self._position), t1)
+
+        def apply_metric(self, x):
+            inp = x.to_global_data_rw()
+            pos = self._position.to_global_data_rw()
+            return ift.Field.from_global_data(space, rosen_hess_prod(pos, inp))
+
+    try:
+        minimizer = eval(minimizer)
+        energy = RBEnergy(position=starting_point)
+
+        (energy, convergence) = minimizer(energy)
+    except NotImplementedError:
+        raise SkipTest
+
+    assert_equal(convergence, IC.CONVERGED)
+    assert_allclose(energy.position.local_data, 1., rtol=1e-3, atol=1e-3)
+
+
+@pmp('minimizer', minimizers + slow_minimizers)
+def test_gauss(minimizer):
+    space = ift.UnstructuredDomain((1,))
+    starting_point = ift.Field.full(space, 3.)
+
+    class ExpEnergy(ift.Energy):
+        def __init__(self, position):
+            super(ExpEnergy, self).__init__(position)
+
+        @property
+        def value(self):
+            x = self.position.to_global_data()[0]
+            return -np.exp(-(x**2))
+
+        @property
+        def gradient(self):
+            x = self.position.to_global_data()[0]
+            return ift.Field.full(self.position.domain, 2*x*np.exp(-(x**2)))
+
+        def apply_metric(self, x):
+            p = self.position.to_global_data()[0]
+            v = (2 - 4*p*p)*np.exp(-p**2)
+            return ift.DiagonalOperator(
+                ift.Field.full(self.position.domain, v))(x)
+
+    try:
+        minimizer = eval(minimizer)
+        energy = ExpEnergy(position=starting_point)
+
+        (energy, convergence) = minimizer(energy)
+    except NotImplementedError:
+        raise SkipTest
+
+    assert_equal(convergence, IC.CONVERGED)
+    assert_allclose(energy.position.local_data, 0., atol=1e-3)
+
+
+@pmp('minimizer', minimizers + newton_minimizers + slow_minimizers)
+def test_cosh(minimizer):
+    space = ift.UnstructuredDomain((1,))
+    starting_point = ift.Field.full(space, 3.)
+
+    class CoshEnergy(ift.Energy):
+        def __init__(self, position):
+            super(CoshEnergy, self).__init__(position)
+
+        @property
+        def value(self):
+            x = self.position.to_global_data()[0]
+            return np.cosh(x)
+
+        @property
+        def gradient(self):
+            x = self.position.to_global_data()[0]
+            return ift.Field.full(self.position.domain, np.sinh(x))
+
+        @property
+        def metric(self):
+            x = self.position.to_global_data()[0]
+            v = np.cosh(x)
+            return ift.DiagonalOperator(
+                ift.Field.full(self.position.domain, v))
+
+        def apply_metric(self, x):
+            p = self.position.to_global_data()[0]
+            v = np.cosh(p)
+            return ift.DiagonalOperator(
+                ift.Field.full(self.position.domain, v))(x)
+
+    try:
+        minimizer = eval(minimizer)
+        energy = CoshEnergy(position=starting_point)
+
+        (energy, convergence) = minimizer(energy)
+    except NotImplementedError:
+        raise SkipTest
+
+    assert_equal(convergence, IC.CONVERGED)
+    assert_allclose(energy.position.local_data, 0., atol=1e-3)
diff --git a/test/test_model_gradients.py b/test/test_model_gradients.py
new file mode 100644
index 0000000000000000000000000000000000000000..376d607605e78c47d27d67cb160e15a3c3147072
--- /dev/null
+++ b/test/test_model_gradients.py
@@ -0,0 +1,113 @@
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+import pytest
+
+import nifty5 as ift
+
+from .common import list2fixture
+
+pmp = pytest.mark.parametrize
+space = list2fixture([
+    ift.GLSpace(15),
+    ift.RGSpace(64, distances=.789),
+    ift.RGSpace([32, 32], distances=.789)
+])
+space1 = space
+seed = list2fixture([4, 78, 23])
+
+
+def _make_linearization(type, space, seed):
+    np.random.seed(seed)
+    S = ift.ScalingOperator(1., space)
+    s = S.draw_sample()
+    if type == "Constant":
+        return ift.Linearization.make_const(s)
+    elif type == "Variable":
+        return ift.Linearization.make_var(s)
+    raise ValueError('unknown type passed')
+
+
+def testBasics(space, seed):
+    var = _make_linearization("Variable", space, seed)
+    model = ift.ScalingOperator(6., var.target)
+    ift.extra.check_value_gradient_consistency(model, var.val)
+
+
+@pmp('type1', ['Variable', 'Constant'])
+@pmp('type2', ['Variable'])
+def testBinary(type1, type2, space, seed):
+    dom1 = ift.MultiDomain.make({'s1': space})
+    # FIXME Remove?
+    lin1 = _make_linearization(type1, dom1, seed)
+    dom2 = ift.MultiDomain.make({'s2': space})
+    # FIXME Remove?
+    lin2 = _make_linearization(type2, dom2, seed)
+
+    dom = ift.MultiDomain.union((dom1, dom2))
+    select_s1 = ift.ducktape(None, dom, "s1")
+    select_s2 = ift.ducktape(None, dom, "s2")
+    model = select_s1*select_s2
+    pos = ift.from_random("normal", dom)
+    ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
+    model = select_s1 + select_s2
+    pos = ift.from_random("normal", dom)
+    ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
+    model = select_s1.scale(3.)
+    pos = ift.from_random("normal", dom1)
+    ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
+    model = ift.ScalingOperator(2.456, space)(select_s1*select_s2)
+    pos = ift.from_random("normal", dom)
+    ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
+    model = ift.sigmoid(
+        ift.ScalingOperator(2.456, space)(select_s1*select_s2))
+    pos = ift.from_random("normal", dom)
+    ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
+    pos = ift.from_random("normal", dom)
+    model = ift.OuterProduct(pos['s1'], ift.makeDomain(space))
+    ift.extra.check_value_gradient_consistency(model, pos['s2'], ntries=20)
+    if isinstance(space, ift.RGSpace):
+        model = ift.FFTOperator(space)(select_s1*select_s2)
+        pos = ift.from_random("normal", dom)
+        ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
+
+
+def testModelLibrary(space, seed):
+    # Tests amplitude model and coorelated field model
+    Npixdof, ceps_a, ceps_k, sm, sv, im, iv = 4, 0.5, 2., 3., 1.5, 1.75, 1.3
+    np.random.seed(seed)
+    model = ift.AmplitudeOperator(space, Npixdof, ceps_a, ceps_k, sm, sv, im,
+                                  iv)
+    S = ift.ScalingOperator(1., model.domain)
+    pos = S.draw_sample()
+    ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
+
+    model2 = ift.CorrelatedField(space, model)
+    S = ift.ScalingOperator(1., model2.domain)
+    pos = S.draw_sample()
+    ift.extra.check_value_gradient_consistency(model2, pos, ntries=20)
+
+
+def testPointModel(space, seed):
+    S = ift.ScalingOperator(1., space)
+    pos = S.draw_sample()
+    alpha = 1.5
+    q = 0.73
+    model = ift.InverseGammaOperator(space, alpha, q)
+    # FIXME All those cdfs and ppfs are not very accurate
+    ift.extra.check_value_gradient_consistency(model, pos, tol=1e-2, ntries=20)
diff --git a/test/test_multi_field.py b/test/test_multi_field.py
index 1a8200ca14c35468f79db3e6dd022b72da60626c..a316b3604421558ffb1c53a3c59b3fae9e7b2fcf 100644
--- a/test/test_multi_field.py
+++ b/test/test_multi_field.py
@@ -15,43 +15,43 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-
-import nifty5 as ift
 import numpy as np
 from numpy.testing import assert_allclose, assert_equal
 
+import nifty5 as ift
+
 dom = ift.makeDomain({"d1": ift.RGSpace(10)})
 
 
-class Test_Functionality(unittest.TestCase):
-    def test_vdot(self):
-        f1 = ift.from_random("normal", domain=dom, dtype=np.complex128)
-        f2 = ift.from_random("normal", domain=dom, dtype=np.complex128)
-        assert_allclose(f1.vdot(f2), np.conj(f2.vdot(f1)))
-
-    def test_func(self):
-        f1 = ift.from_random("normal", domain=dom, dtype=np.complex128)
-        assert_allclose(ift.log(ift.exp((f1)))["d1"].local_data,
-                        f1["d1"].local_data)
-
-    def test_dataconv(self):
-        f1 = ift.full(dom, 27)
-        f2 = ift.from_global_data(dom, f1.to_global_data())
-        for key, val in f1.items():
-            assert_equal(val.local_data, f2[key].local_data)
-
-    def test_blockdiagonal(self):
-        op = ift.BlockDiagonalOperator(
-            dom, (ift.ScalingOperator(20., dom["d1"]),))
-        op2 = op(op)
-        ift.extra.consistency_check(op2)
-        assert_equal(type(op2), ift.BlockDiagonalOperator)
-        f1 = op2(ift.full(dom, 1))
-        for val in f1.values():
-            assert_equal((val == 400).all(), True)
-        op2 = op+op
-        assert_equal(type(op2), ift.BlockDiagonalOperator)
-        f1 = op2(ift.full(dom, 1))
-        for val in f1.values():
-            assert_equal((val == 40).all(), True)
+def test_vdot():
+    f1 = ift.from_random("normal", domain=dom, dtype=np.complex128)
+    f2 = ift.from_random("normal", domain=dom, dtype=np.complex128)
+    assert_allclose(f1.vdot(f2), np.conj(f2.vdot(f1)))
+
+
+def test_func():
+    f1 = ift.from_random("normal", domain=dom, dtype=np.complex128)
+    assert_allclose(
+        ift.log(ift.exp((f1)))["d1"].local_data, f1["d1"].local_data)
+
+
+def test_dataconv():
+    f1 = ift.full(dom, 27)
+    f2 = ift.from_global_data(dom, f1.to_global_data())
+    for key, val in f1.items():
+        assert_equal(val.local_data, f2[key].local_data)
+
+
+def test_blockdiagonal():
+    op = ift.BlockDiagonalOperator(dom, (ift.ScalingOperator(20., dom["d1"]),))
+    op2 = op(op)
+    ift.extra.consistency_check(op2)
+    assert_equal(type(op2), ift.BlockDiagonalOperator)
+    f1 = op2(ift.full(dom, 1))
+    for val in f1.values():
+        assert_equal((val == 400).all(), True)
+    op2 = op + op
+    assert_equal(type(op2), ift.BlockDiagonalOperator)
+    f1 = op2(ift.full(dom, 1))
+    for val in f1.values():
+        assert_equal((val == 40).all(), True)
diff --git a/test/test_operators/__init__.py b/test/test_operators/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py
index 908089eb42e89dd295c47a394680c27fe8e49b4a..a3d6b5cf54d660f250ea0969ba9e9194bedc1afa 100644
--- a/test/test_operators/test_adjoint.py
+++ b/test/test_operators/test_adjoint.py
@@ -15,248 +15,262 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
+import numpy as np
+import pytest
 
 import nifty5 as ift
-import numpy as np
 
-_h_RG_spaces = [ift.RGSpace(7, distances=0.2, harmonic=True),
-                ift.RGSpace((12, 46), distances=(.2, .3), harmonic=True)]
-_h_spaces = _h_RG_spaces + [ift.LMSpace(17)]
+from ..common import list2fixture
 
-_p_RG_spaces = [ift.RGSpace(19, distances=0.7),
-                ift.RGSpace((1, 2, 3, 6), distances=(0.2, 0.25, 0.34, .8))]
+_h_RG_spaces = [
+    ift.RGSpace(7, distances=0.2, harmonic=True),
+    ift.RGSpace((12, 46), distances=(.2, .3), harmonic=True)
+]
+_h_spaces = _h_RG_spaces + [ift.LMSpace(17)]
+_p_RG_spaces = [
+    ift.RGSpace(19, distances=0.7),
+    ift.RGSpace((1, 2, 3, 6), distances=(0.2, 0.25, 0.34, .8))
+]
 _p_spaces = _p_RG_spaces + [ift.HPSpace(17), ift.GLSpace(8, 13)]
-
 _pow_spaces = [ift.PowerSpace(ift.RGSpace((17, 38), harmonic=True))]
 
+pmp = pytest.mark.parametrize
+dtype = list2fixture([np.float64, np.complex128])
+
+
+@pmp('sp', _p_RG_spaces)
+def testLOSResponse(sp, dtype):
+    starts = np.random.randn(len(sp.shape), 10)
+    ends = np.random.randn(len(sp.shape), 10)
+    sigma_low = 1e-4*np.random.randn(10)
+    sigma_ups = 1e-5*np.random.randn(10)
+    op = ift.LOSResponse(sp, starts, ends, sigma_low, sigma_ups)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
+def testOperatorCombinations(sp, dtype):
+    a = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
+    b = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
+    op = ift.SandwichOperator.make(a, b)
+    ift.extra.consistency_check(op, dtype, dtype)
+    op = a(b)
+    ift.extra.consistency_check(op, dtype, dtype)
+    op = a + b
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+def testLinearInterpolator():
+    sp = ift.RGSpace((10, 8), distances=(0.1, 3.5))
+    pos = np.random.rand(2, 23)
+    pos[0, :] *= 0.9
+    pos[1, :] *= 7*3.5
+    op = ift.LinearInterpolator(sp, pos)
+    ift.extra.consistency_check(op)
+
+
+@pmp('args', [(ift.RGSpace(10, harmonic=True), 4, 0), (ift.RGSpace(
+    (24, 31), distances=(0.4, 2.34), harmonic=True), 3, 0),
+              (ift.LMSpace(4), 10, 0)])
+def testSlopeOperator(args, dtype):
+    tmp = ift.ExpTransform(ift.PowerSpace(args[0]), args[1], args[2])
+    tgt = tmp.domain[0]
+    op = ift.SlopeOperator(tgt)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
+def testOperatorAdaptor(sp, dtype):
+    op = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
+    ift.extra.consistency_check(op.adjoint, dtype, dtype)
+    ift.extra.consistency_check(op.inverse, dtype, dtype)
+    ift.extra.consistency_check(op.inverse.adjoint, dtype, dtype)
+    ift.extra.consistency_check(op.adjoint.inverse, dtype, dtype)
+
+
+@pmp('sp1', _h_spaces + _p_spaces + _pow_spaces)
+@pmp('sp2', _h_spaces + _p_spaces + _pow_spaces)
+def testNullOperator(sp1, sp2, dtype):
+    op = ift.NullOperator(sp1, sp2)
+    ift.extra.consistency_check(op, dtype, dtype)
+    mdom1 = ift.MultiDomain.make({'a': sp1})
+    mdom2 = ift.MultiDomain.make({'b': sp2})
+    op = ift.NullOperator(mdom1, mdom2)
+    ift.extra.consistency_check(op, dtype, dtype)
+    op = ift.NullOperator(sp1, mdom2)
+    ift.extra.consistency_check(op, dtype, dtype)
+    op = ift.NullOperator(mdom1, sp2)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _p_RG_spaces)
+def testHarmonicSmoothingOperator(sp, dtype):
+    op = ift.HarmonicSmoothingOperator(sp, 0.1)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
+def testDOFDistributor(sp, dtype):
+    # TODO: Test for DomainTuple
+    if sp.size < 4:
+        return
+    dofdex = np.arange(sp.size).reshape(sp.shape) % 3
+    dofdex = ift.Field.from_global_data(sp, dofdex)
+    op = ift.DOFDistributor(dofdex)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _h_spaces)
+def testPPO(sp, dtype):
+    op = ift.PowerDistributor(target=sp)
+    ift.extra.consistency_check(op, dtype, dtype)
+    ps = ift.PowerSpace(
+        sp, ift.PowerSpace.useful_binbounds(sp, logarithmic=False, nbin=3))
+    op = ift.PowerDistributor(target=sp, power_space=ps)
+    ift.extra.consistency_check(op, dtype, dtype)
+    ps = ift.PowerSpace(
+        sp, ift.PowerSpace.useful_binbounds(sp, logarithmic=True, nbin=3))
+    op = ift.PowerDistributor(target=sp, power_space=ps)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _h_RG_spaces + _p_RG_spaces)
+def testFFT(sp, dtype):
+    op = ift.FFTOperator(sp)
+    ift.extra.consistency_check(op, dtype, dtype)
+    op = ift.FFTOperator(sp.get_default_codomain())
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _h_RG_spaces + _p_RG_spaces)
+def testHartley(sp, dtype):
+    op = ift.HartleyOperator(sp)
+    ift.extra.consistency_check(op, dtype, dtype)
+    op = ift.HartleyOperator(sp.get_default_codomain())
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _h_spaces)
+def testHarmonic(sp, dtype):
+    op = ift.HarmonicTransformOperator(sp)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _p_spaces)
+def testMask(sp, dtype):
+    # Create mask
+    f = ift.from_random('normal', sp).to_global_data()
+    mask = np.zeros_like(f)
+    mask[f > 0] = 1
+    mask = ift.Field.from_global_data(sp, mask)
+    # Test MaskOperator
+    op = ift.MaskOperator(mask)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _h_spaces + _p_spaces)
+def testDiagonal(sp, dtype):
+    op = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
+def testGeometryRemover(sp, dtype):
+    op = ift.GeometryRemover(sp)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('spaces', [0, 1, 2, 3, (0, 1), (0, 2), (0, 1, 2), (0, 2, 3), (1, 3)])
+@pmp('wgt', [0, 1, 2, -1])
+def testContractionOperator(spaces, wgt, dtype):
+    dom = (ift.RGSpace(10), ift.RGSpace(13), ift.GLSpace(5), ift.HPSpace(4))
+    op = ift.ContractionOperator(dom, spaces, wgt)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+def testDomainTupleFieldInserter():
+    domain = ift.DomainTuple.make((ift.UnstructuredDomain(12),
+                                   ift.RGSpace([4, 22])))
+    new_space = ift.UnstructuredDomain(7)
+    pos = (5,)
+    op = ift.DomainTupleFieldInserter(domain, new_space, 0, pos)
+    ift.extra.consistency_check(op)
+
+
+@pmp('space', [0, 2])
+def testSymmetrizingOperator(space, dtype):
+    dom = (ift.LogRGSpace(10, [2.], [1.]), ift.UnstructuredDomain(13),
+           ift.LogRGSpace((5, 27), [1., 2.7], [0., 4.]), ift.HPSpace(4))
+    op = ift.SymmetrizingOperator(dom, space)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('space', [0, 2])
+@pmp('factor', [1, 2, 2.7])
+@pmp('central', [False, True])
+def testZeroPadder(space, factor, dtype, central):
+    dom = (ift.RGSpace(10), ift.UnstructuredDomain(13), ift.RGSpace(7, 12),
+           ift.HPSpace(4))
+    newshape = [int(factor*l) for l in dom[space].shape]
+    op = ift.FieldZeroPadder(dom, newshape, space, central)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('space', [0, 2])
+@pmp('factor', [2, 2.7])
+def testZeroPadder2(space, factor, dtype):
+    dom = (ift.RGSpace(10), ift.UnstructuredDomain(13), ift.RGSpace(7, 12),
+           ift.HPSpace(4))
+    newshape = [int(factor*l) for l in dom[space].shape]
+    op = ift.CentralZeroPadder(dom, newshape, space)
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('args',
+     [(ift.RGSpace(10, harmonic=True), 4, 0), (ift.RGSpace(
+         (24, 31), distances=(0.4, 2.34), harmonic=True), (4, 3), 0),
+      ((ift.HPSpace(4), ift.RGSpace(27, distances=0.3, harmonic=True)),
+       (10,), 1),
+      (ift.PowerSpace(ift.RGSpace(10, distances=0.3, harmonic=True)), 6, 0)])
+def testExpTransform(args, dtype):
+    op = ift.ExpTransform(args[0], args[1], args[2])
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('args',
+     [(ift.LogRGSpace([10, 17], [2., 3.], [1., 0.]), 0),
+      ((ift.LogRGSpace(10, [2.], [1.]), ift.UnstructuredDomain(13)), 0),
+      ((ift.UnstructuredDomain(13), ift.LogRGSpace(17, [3.], [.7])), 1)])
+def testQHTOperator(args):
+    dtype = np.float64
+    tgt = ift.DomainTuple.make(args[0])
+    op = ift.QHTOperator(tgt, args[1])
+    ift.extra.consistency_check(op, dtype, dtype)
+
+
+@pmp('args', [[ift.RGSpace(
+    (13, 52, 40)), (4, 6, 25), None], [ift.RGSpace(
+        (128, 128)), (45, 48), 0], [ift.RGSpace(13), (7,), None], [
+            (ift.HPSpace(3), ift.RGSpace((12, 24), distances=0.3)), (12, 12), 1
+        ]])
+def testRegridding(args):
+    op = ift.RegriddingOperator(*args)
+    ift.extra.consistency_check(op)
+
 
-class Consistency_Tests(unittest.TestCase):
-    @expand(product(_p_RG_spaces, [np.float64, np.complex128]))
-    def testLOSResponse(self, sp, dtype):
-        starts = np.random.randn(len(sp.shape), 10)
-        ends = np.random.randn(len(sp.shape), 10)
-        sigma_low = 1e-4*np.random.randn(10)
-        sigma_ups = 1e-5*np.random.randn(10)
-        op = ift.LOSResponse(sp, starts, ends, sigma_low, sigma_ups)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_h_spaces + _p_spaces + _pow_spaces,
-                    [np.float64, np.complex128]))
-    def testOperatorCombinations(self, sp, dtype):
-        a = ift.DiagonalOperator(ift.Field.from_random("normal", sp,
-                                                       dtype=dtype))
-        b = ift.DiagonalOperator(ift.Field.from_random("normal", sp,
-                                                       dtype=dtype))
-        op = ift.SandwichOperator.make(a, b)
-        ift.extra.consistency_check(op, dtype, dtype)
-        op = a(b)
-        ift.extra.consistency_check(op, dtype, dtype)
-        op = a+b
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    def testLinearInterpolator(self):
-        sp = ift.RGSpace((10, 8), distances=(0.1, 3.5))
-        pos = np.random.rand(2, 23)
-        pos[0, :] *= 0.9
-        pos[1, :] *= 7*3.5
-        op = ift.LinearInterpolator(sp, pos)
-        ift.extra.consistency_check(op)
-
-    @expand(product([(ift.RGSpace(10, harmonic=True), 4, 0),
-                     (ift.RGSpace((24, 31), distances=(0.4, 2.34),
-                                  harmonic=True), 3, 0),
-                     (ift.LMSpace(4), 10, 0)],
-                    [np.float64, np.complex128]))
-    def testSlopeOperator(self, args, dtype):
-        tmp = ift.ExpTransform(ift.PowerSpace(args[0]), args[1], args[2])
-        tgt = tmp.domain[0]
-        op = ift.SlopeOperator(tgt)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_h_spaces + _p_spaces + _pow_spaces,
-                    [np.float64, np.complex128]))
-    def testOperatorAdaptor(self, sp, dtype):
-        op = ift.DiagonalOperator(ift.Field.from_random("normal", sp,
-                                                        dtype=dtype))
-        ift.extra.consistency_check(op.adjoint, dtype, dtype)
-        ift.extra.consistency_check(op.inverse, dtype, dtype)
-        ift.extra.consistency_check(op.inverse.adjoint, dtype, dtype)
-        ift.extra.consistency_check(op.adjoint.inverse, dtype, dtype)
-
-    @expand(product(_h_spaces + _p_spaces + _pow_spaces,
-                    _h_spaces + _p_spaces + _pow_spaces,
-                    [np.float64, np.complex128]))
-    def testNullOperator(self, sp1, sp2, dtype):
-        op = ift.NullOperator(sp1, sp2)
-        ift.extra.consistency_check(op, dtype, dtype)
-        mdom1 = ift.MultiDomain.make({'a': sp1})
-        mdom2 = ift.MultiDomain.make({'b': sp2})
-        op = ift.NullOperator(mdom1, mdom2)
-        ift.extra.consistency_check(op, dtype, dtype)
-        op = ift.NullOperator(sp1, mdom2)
-        ift.extra.consistency_check(op, dtype, dtype)
-        op = ift.NullOperator(mdom1, sp2)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_p_RG_spaces,
-                    [np.float64, np.complex128]))
-    def testHarmonicSmoothingOperator(self, sp, dtype):
-        op = ift.HarmonicSmoothingOperator(sp, 0.1)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_h_spaces + _p_spaces + _pow_spaces,
-                    [np.float64, np.complex128]))
-    def testDOFDistributor(self, sp, dtype):
-        # TODO: Test for DomainTuple
-        if sp.size < 4:
-            return
-        dofdex = np.arange(sp.size).reshape(sp.shape) % 3
-        dofdex = ift.Field.from_global_data(sp, dofdex)
-        op = ift.DOFDistributor(dofdex)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_h_spaces, [np.float64, np.complex128]))
-    def testPPO(self, sp, dtype):
-        op = ift.PowerDistributor(target=sp)
-        ift.extra.consistency_check(op, dtype, dtype)
-        ps = ift.PowerSpace(
-            sp, ift.PowerSpace.useful_binbounds(sp, logarithmic=False, nbin=3))
-        op = ift.PowerDistributor(target=sp, power_space=ps)
-        ift.extra.consistency_check(op, dtype, dtype)
-        ps = ift.PowerSpace(
-            sp, ift.PowerSpace.useful_binbounds(sp, logarithmic=True, nbin=3))
-        op = ift.PowerDistributor(target=sp, power_space=ps)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_h_RG_spaces+_p_RG_spaces,
-                    [np.float64, np.complex128]))
-    def testFFT(self, sp, dtype):
-        op = ift.FFTOperator(sp)
-        ift.extra.consistency_check(op, dtype, dtype)
-        op = ift.FFTOperator(sp.get_default_codomain())
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_h_RG_spaces+_p_RG_spaces,
-                    [np.float64, np.complex128]))
-    def testHartley(self, sp, dtype):
-        op = ift.HartleyOperator(sp)
-        ift.extra.consistency_check(op, dtype, dtype)
-        op = ift.HartleyOperator(sp.get_default_codomain())
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_h_spaces, [np.float64, np.complex128]))
-    def testHarmonic(self, sp, dtype):
-        op = ift.HarmonicTransformOperator(sp)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_p_spaces, [np.float64, np.complex128]))
-    def testMask(self, sp, dtype):
-        # Create mask
-        f = ift.from_random('normal', sp).to_global_data()
-        mask = np.zeros_like(f)
-        mask[f > 0] = 1
-        mask = ift.Field.from_global_data(sp, mask)
-        # Test MaskOperator
-        op = ift.MaskOperator(mask)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_h_spaces+_p_spaces, [np.float64, np.complex128]))
-    def testDiagonal(self, sp, dtype):
-        op = ift.DiagonalOperator(ift.Field.from_random("normal", sp,
-                                                        dtype=dtype))
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product(_h_spaces+_p_spaces+_pow_spaces,
-                    [np.float64, np.complex128]))
-    def testGeometryRemover(self, sp, dtype):
-        op = ift.GeometryRemover(sp)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product([0, 1, 2, 3, (0, 1), (0, 2), (0, 1, 2), (0, 2, 3), (1, 3)],
-                    [0, 1, 2, -1], [np.float64, np.complex128]))
-    def testContractionOperator(self, spaces, wgt, dtype):
-        dom = (ift.RGSpace(10), ift.RGSpace(13), ift.GLSpace(5),
-               ift.HPSpace(4))
-        op = ift.ContractionOperator(dom, spaces, wgt)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    def testDomainTupleFieldInserter(self):
-        domain = ift.DomainTuple.make((ift.UnstructuredDomain(12),
-                                       ift.RGSpace([4, 22])))
-        new_space = ift.UnstructuredDomain(7)
-        pos = (5,)
-        op = ift.DomainTupleFieldInserter(domain, new_space, 0, pos)
-        ift.extra.consistency_check(op)
-
-    @expand(product([0, 2], [np.float64, np.complex128]))
-    def testSymmetrizingOperator(self, space, dtype):
-        dom = (ift.LogRGSpace(10, [2.], [1.]), ift.UnstructuredDomain(13),
-               ift.LogRGSpace((5, 27), [1., 2.7], [0., 4.]), ift.HPSpace(4))
-        op = ift.SymmetrizingOperator(dom, space)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product([0, 2], [1, 2, 2.7], [np.float64, np.complex128],
-                    [False, True]))
-    def testZeroPadder(self, space, factor, dtype, central):
-        dom = (ift.RGSpace(10), ift.UnstructuredDomain(13), ift.RGSpace(7, 12),
-               ift.HPSpace(4))
-        newshape = [int(factor*l) for l in dom[space].shape]
-        op = ift.FieldZeroPadder(dom, newshape, space, central)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product([0, 2], [2, 2.7], [np.float64, np.complex128]))
-    def testZeroPadder2(self, space, factor, dtype):
-        dom = (ift.RGSpace(10), ift.UnstructuredDomain(13), ift.RGSpace(7, 12),
-               ift.HPSpace(4))
-        newshape = [int(factor*l) for l in dom[space].shape]
-        op = ift.CentralZeroPadder(dom, newshape, space)
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product([(ift.RGSpace(10, harmonic=True), 4, 0),
-                     (ift.RGSpace((24, 31), distances=(0.4, 2.34),
-                                  harmonic=True), (4, 3), 0),
-                     ((ift.HPSpace(4), ift.RGSpace(27, distances=0.3,
-                                                   harmonic=True)), (10,), 1),
-                     (ift.PowerSpace(ift.RGSpace(10, distances=0.3,
-                                     harmonic=True)), 6, 0)],
-                    [np.float64, np.complex128]))
-    def testExpTransform(self, args, dtype):
-        op = ift.ExpTransform(args[0], args[1], args[2])
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand(product([(ift.LogRGSpace([10, 17], [2., 3.], [1., 0.]), 0),
-                     ((ift.LogRGSpace(10, [2.], [1.]),
-                       ift.UnstructuredDomain(13)), 0),
-                     ((ift.UnstructuredDomain(13),
-                       ift.LogRGSpace(17, [3.], [.7])), 1)],
-                    [np.float64]))
-    def testQHTOperator(self, args, dtype):
-        tgt = ift.DomainTuple.make(args[0])
-        op = ift.QHTOperator(tgt, args[1])
-        ift.extra.consistency_check(op, dtype, dtype)
-
-    @expand([[ift.RGSpace((13, 52, 40)), (4, 6, 25), None],
-             [ift.RGSpace((128, 128)), (45, 48), 0],
-             [ift.RGSpace(13), (7,), None],
-             [(ift.HPSpace(3), ift.RGSpace((12, 24), distances=0.3)),
-              (12, 12), 1]])
-    def testRegridding(self, domain, shape, space):
-        op = ift.RegriddingOperator(domain, shape, space)
-        ift.extra.consistency_check(op)
-
-    @expand(product([ift.DomainTuple.make((ift.RGSpace((3, 5, 4)),
-                                           ift.RGSpace((16,),
-                                           distances=(7.,))),),
-                     ift.DomainTuple.make(ift.HPSpace(12),)],
-                    [ift.DomainTuple.make((ift.RGSpace((2,)),
-                                           ift.GLSpace(10)),),
-                     ift.DomainTuple.make(ift.RGSpace((10, 12),
-                                          distances=(0.1, 1.)),)]
-                    ))
-    def testOuter(self, fdomain, domain):
-        f = ift.from_random('normal', fdomain)
-        op = ift.OuterProduct(f, domain)
-        ift.extra.consistency_check(op)
+@pmp(
+    'fdomain',
+    [
+        ift.DomainTuple.make((ift.RGSpace(
+            (3, 5, 4)), ift.RGSpace((16,), distances=(7.,))),),
+        ift.DomainTuple.make(ift.HPSpace(12),)
+    ],
+)
+@pmp('domain', [
+    ift.DomainTuple.make((ift.RGSpace((2,)), ift.GLSpace(10)),),
+    ift.DomainTuple.make(ift.RGSpace((10, 12), distances=(0.1, 1.)),)
+])
+def testOuter(fdomain, domain):
+    f = ift.from_random('normal', fdomain)
+    op = ift.OuterProduct(f, domain)
+    ift.extra.consistency_check(op)
diff --git a/test/test_operators/test_composed_operator.py b/test/test_operators/test_composed_operator.py
index 63cd7e4f0e2aaae2237fb35513b46cc0cf722d24..c89ca7bfb953c38e6b408e5eb9688ac73cdae21c 100644
--- a/test/test_operators/test_composed_operator.py
+++ b/test/test_operators/test_composed_operator.py
@@ -15,77 +15,79 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
+from numpy.testing import assert_allclose, assert_equal
 
 import nifty5 as ift
-from numpy.testing import assert_allclose, assert_equal
+
+from ..common import list2fixture
+
+space1 = list2fixture([
+    ift.RGSpace(4),
+    ift.PowerSpace(ift.RGSpace((4, 4), harmonic=True)),
+    ift.LMSpace(5),
+    ift.HPSpace(4),
+    ift.GLSpace(4)
+])
+space2 = space1
+
+
+def test_times_adjoint_times(space1, space2):
+    cspace = (space1, space2)
+    diag1 = ift.Field.from_random('normal', domain=space1)
+    diag2 = ift.Field.from_random('normal', domain=space2)
+    op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
+    op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
+
+    op = op2(op1)
+
+    rand1 = ift.Field.from_random('normal', domain=(space1, space2))
+    rand2 = ift.Field.from_random('normal', domain=(space1, space2))
+
+    tt1 = rand2.vdot(op.times(rand1))
+    tt2 = rand1.vdot(op.adjoint_times(rand2))
+    assert_allclose(tt1, tt2)
+
+
+def test_times_inverse_times(space1, space2):
+    cspace = (space1, space2)
+    diag1 = ift.Field.from_random('normal', domain=space1)
+    diag2 = ift.Field.from_random('normal', domain=space2)
+    op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
+    op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
+
+    op = op2(op1)
+
+    rand1 = ift.Field.from_random('normal', domain=(space1, space2))
+    tt1 = op.inverse_times(op.times(rand1))
+
+    assert_allclose(tt1.local_data, rand1.local_data)
+
+
+def test_sum(space1):
+    op1 = ift.makeOp(ift.Field.full(space1, 2.))
+    op2 = ift.ScalingOperator(3., space1)
+    full_op = op1 + op2 - (op2 - op1) + op1 + op1 + op2
+    x = ift.Field.full(space1, 1.)
+    res = full_op(x)
+    assert_equal(isinstance(full_op, ift.DiagonalOperator), True)
+    assert_allclose(res.local_data, 11.)
+
+
+def test_chain(space1):
+    op1 = ift.makeOp(ift.Field.full(space1, 2.))
+    op2 = ift.ScalingOperator(3., space1)
+    full_op = op1(op2)(op2)(op1)(op1)(op1)(op2)
+    x = ift.Field.full(space1, 1.)
+    res = full_op(x)
+    assert_equal(isinstance(full_op, ift.DiagonalOperator), True)
+    assert_allclose(res.local_data, 432.)
 
 
-class ComposedOperator_Tests(unittest.TestCase):
-    spaces = [ift.RGSpace(4),
-              ift.PowerSpace(ift.RGSpace((4, 4), harmonic=True)),
-              ift.LMSpace(5), ift.HPSpace(4), ift.GLSpace(4)]
-
-    @expand(product(spaces, spaces))
-    def test_times_adjoint_times(self, space1, space2):
-        cspace = (space1, space2)
-        diag1 = ift.Field.from_random('normal', domain=space1)
-        diag2 = ift.Field.from_random('normal', domain=space2)
-        op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
-        op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
-
-        op = op2(op1)
-
-        rand1 = ift.Field.from_random('normal', domain=(space1, space2))
-        rand2 = ift.Field.from_random('normal', domain=(space1, space2))
-
-        tt1 = rand2.vdot(op.times(rand1))
-        tt2 = rand1.vdot(op.adjoint_times(rand2))
-        assert_allclose(tt1, tt2)
-
-    @expand(product(spaces, spaces))
-    def test_times_inverse_times(self, space1, space2):
-        cspace = (space1, space2)
-        diag1 = ift.Field.from_random('normal', domain=space1)
-        diag2 = ift.Field.from_random('normal', domain=space2)
-        op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
-        op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
-
-        op = op2(op1)
-
-        rand1 = ift.Field.from_random('normal', domain=(space1, space2))
-        tt1 = op.inverse_times(op.times(rand1))
-
-        assert_allclose(tt1.local_data, rand1.local_data)
-
-    @expand(product(spaces))
-    def test_sum(self, space):
-        op1 = ift.makeOp(ift.Field.full(space, 2.))
-        op2 = ift.ScalingOperator(3., space)
-        full_op = op1 + op2 - (op2 - op1) + op1 + op1 + op2
-        x = ift.Field.full(space, 1.)
-        res = full_op(x)
-        assert_equal(isinstance(full_op, ift.DiagonalOperator), True)
-        assert_allclose(res.local_data, 11.)
-
-    @expand(product(spaces))
-    def test_chain(self, space):
-        op1 = ift.makeOp(ift.Field.full(space, 2.))
-        op2 = ift.ScalingOperator(3., space)
-        full_op = op1(op2)(op2)(op1)(op1)(op1)(op2)
-        x = ift.Field.full(space, 1.)
-        res = full_op(x)
-        assert_equal(isinstance(full_op, ift.DiagonalOperator), True)
-        assert_allclose(res.local_data, 432.)
-
-    @expand(product(spaces))
-    def test_mix(self, space):
-        op1 = ift.makeOp(ift.Field.full(space, 2.))
-        op2 = ift.ScalingOperator(3., space)
-        full_op = op1(op2+op2)(op1)(op1) - op1(op2)
-        x = ift.Field.full(space, 1.)
-        res = full_op(x)
-        assert_equal(isinstance(full_op, ift.DiagonalOperator), True)
-        assert_allclose(res.local_data, 42.)
+def test_mix(space1):
+    op1 = ift.makeOp(ift.Field.full(space1, 2.))
+    op2 = ift.ScalingOperator(3., space1)
+    full_op = op1(op2 + op2)(op1)(op1) - op1(op2)
+    x = ift.Field.full(space1, 1.)
+    res = full_op(x)
+    assert_equal(isinstance(full_op, ift.DiagonalOperator), True)
+    assert_allclose(res.local_data, 42.)
diff --git a/test/test_operators/test_diagonal_operator.py b/test/test_operators/test_diagonal_operator.py
index 6e6114abd95a7981f5c0b0be4b8406e165025522..52f104f85691aecdbdfab786ee9cf4d111a26a44 100644
--- a/test/test_operators/test_diagonal_operator.py
+++ b/test/test_operators/test_diagonal_operator.py
@@ -15,79 +15,80 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
+from numpy.testing import assert_allclose, assert_equal
 
 import nifty5 as ift
-from numpy.testing import assert_allclose, assert_equal
+
+from ..common import list2fixture
+
+space = list2fixture([
+    ift.RGSpace(4),
+    ift.PowerSpace(ift.RGSpace((4, 4), harmonic=True)),
+    ift.LMSpace(5),
+    ift.HPSpace(4),
+    ift.GLSpace(4)
+])
+
+
+def test_property(space):
+    diag = ift.Field.from_random('normal', domain=space)
+    D = ift.DiagonalOperator(diag)
+    if D.domain[0] != space:
+        raise TypeError
+
+
+def test_times_adjoint(space):
+    rand1 = ift.Field.from_random('normal', domain=space)
+    rand2 = ift.Field.from_random('normal', domain=space)
+    diag = ift.Field.from_random('normal', domain=space)
+    D = ift.DiagonalOperator(diag)
+    tt1 = rand1.vdot(D.times(rand2))
+    tt2 = rand2.vdot(D.times(rand1))
+    assert_allclose(tt1, tt2)
+
+
+def test_times_inverse(space):
+    rand1 = ift.Field.from_random('normal', domain=space)
+    diag = ift.Field.from_random('normal', domain=space)
+    D = ift.DiagonalOperator(diag)
+    tt1 = D.times(D.inverse_times(rand1))
+    assert_allclose(rand1.local_data, tt1.local_data)
+
+
+def test_times(space):
+    rand1 = ift.Field.from_random('normal', domain=space)
+    diag = ift.Field.from_random('normal', domain=space)
+    D = ift.DiagonalOperator(diag)
+    tt = D.times(rand1)
+    assert_equal(tt.domain[0], space)
+
+
+def test_adjoint_times(space):
+    rand1 = ift.Field.from_random('normal', domain=space)
+    diag = ift.Field.from_random('normal', domain=space)
+    D = ift.DiagonalOperator(diag)
+    tt = D.adjoint_times(rand1)
+    assert_equal(tt.domain[0], space)
+
+
+def test_inverse_times(space):
+    rand1 = ift.Field.from_random('normal', domain=space)
+    diag = ift.Field.from_random('normal', domain=space)
+    D = ift.DiagonalOperator(diag)
+    tt = D.inverse_times(rand1)
+    assert_equal(tt.domain[0], space)
+
+
+def test_adjoint_inverse_times(space):
+    rand1 = ift.Field.from_random('normal', domain=space)
+    diag = ift.Field.from_random('normal', domain=space)
+    D = ift.DiagonalOperator(diag)
+    tt = D.adjoint_inverse_times(rand1)
+    assert_equal(tt.domain[0], space)
 
 
-class DiagonalOperator_Tests(unittest.TestCase):
-    spaces = [ift.RGSpace(4),
-              ift.PowerSpace(ift.RGSpace((4, 4), harmonic=True)),
-              ift.LMSpace(5), ift.HPSpace(4), ift.GLSpace(4)]
-
-    @expand(product(spaces))
-    def test_property(self, space):
-        diag = ift.Field.from_random('normal', domain=space)
-        D = ift.DiagonalOperator(diag)
-        if D.domain[0] != space:
-            raise TypeError
-
-    @expand(product(spaces))
-    def test_times_adjoint(self, space):
-        rand1 = ift.Field.from_random('normal', domain=space)
-        rand2 = ift.Field.from_random('normal', domain=space)
-        diag = ift.Field.from_random('normal', domain=space)
-        D = ift.DiagonalOperator(diag)
-        tt1 = rand1.vdot(D.times(rand2))
-        tt2 = rand2.vdot(D.times(rand1))
-        assert_allclose(tt1, tt2)
-
-    @expand(product(spaces))
-    def test_times_inverse(self, space):
-        rand1 = ift.Field.from_random('normal', domain=space)
-        diag = ift.Field.from_random('normal', domain=space)
-        D = ift.DiagonalOperator(diag)
-        tt1 = D.times(D.inverse_times(rand1))
-        assert_allclose(rand1.local_data, tt1.local_data)
-
-    @expand(product(spaces))
-    def test_times(self, space):
-        rand1 = ift.Field.from_random('normal', domain=space)
-        diag = ift.Field.from_random('normal', domain=space)
-        D = ift.DiagonalOperator(diag)
-        tt = D.times(rand1)
-        assert_equal(tt.domain[0], space)
-
-    @expand(product(spaces))
-    def test_adjoint_times(self, space):
-        rand1 = ift.Field.from_random('normal', domain=space)
-        diag = ift.Field.from_random('normal', domain=space)
-        D = ift.DiagonalOperator(diag)
-        tt = D.adjoint_times(rand1)
-        assert_equal(tt.domain[0], space)
-
-    @expand(product(spaces))
-    def test_inverse_times(self, space):
-        rand1 = ift.Field.from_random('normal', domain=space)
-        diag = ift.Field.from_random('normal', domain=space)
-        D = ift.DiagonalOperator(diag)
-        tt = D.inverse_times(rand1)
-        assert_equal(tt.domain[0], space)
-
-    @expand(product(spaces))
-    def test_adjoint_inverse_times(self, space):
-        rand1 = ift.Field.from_random('normal', domain=space)
-        diag = ift.Field.from_random('normal', domain=space)
-        D = ift.DiagonalOperator(diag)
-        tt = D.adjoint_inverse_times(rand1)
-        assert_equal(tt.domain[0], space)
-
-    @expand(product(spaces))
-    def test_diagonal(self, space):
-        diag = ift.Field.from_random('normal', domain=space)
-        D = ift.DiagonalOperator(diag)
-        diag_op = D(ift.Field.full(space, 1.))
-        assert_allclose(diag.local_data, diag_op.local_data)
+def test_diagonal(space):
+    diag = ift.Field.from_random('normal', domain=space)
+    D = ift.DiagonalOperator(diag)
+    diag_op = D(ift.Field.full(space, 1.))
+    assert_allclose(diag.local_data, diag_op.local_data)
diff --git a/test/test_operators/test_fft_operator.py b/test/test_operators/test_fft_operator.py
index 2de51d3ab46a7a2d09d7c5116b96769d3e21235c..171f26bfab0ce31629c5a3dcd9560052c4715b16 100644
--- a/test/test_operators/test_fft_operator.py
+++ b/test/test_operators/test_fft_operator.py
@@ -15,14 +15,14 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
 import numpy as np
+import pytest
 from numpy.testing import assert_allclose
 
+import nifty5 as ift
+
+from ..common import list2fixture
+
 
 def _get_rtol(tp):
     if (tp == np.float64) or (tp == np.complex128):
@@ -31,83 +31,89 @@ def _get_rtol(tp):
         return 1e-5
 
 
-class FFTOperatorTests(unittest.TestCase):
-    @expand(product([16, ], [0.1, 1, 3.7],
-                    [np.float64, np.float32, np.complex64, np.complex128],
-                    [ift.HartleyOperator, ift.FFTOperator]))
-    def test_fft1D(self, dim1, d, itp, op):
-        tol = _get_rtol(itp)
-        a = ift.RGSpace(dim1, distances=d)
-        b = ift.RGSpace(dim1, distances=1./(dim1*d), harmonic=True)
-        np.random.seed(16)
-
-        fft = op(domain=a, target=b)
-        inp = ift.Field.from_random(domain=a, random_type='normal',
-                                    std=7, mean=3, dtype=itp)
-        out = fft.inverse_times(fft.times(inp))
-        assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
-
-        a, b = b, a
-
-        fft = ift.FFTOperator(domain=a, target=b)
-        inp = ift.Field.from_random(domain=a, random_type='normal',
-                                    std=7, mean=3, dtype=itp)
-        out = fft.inverse_times(fft.times(inp))
-        assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
-
-    @expand(product([12, 15], [9, 12], [0.1, 1, 3.7],
-                    [0.4, 1, 2.7],
-                    [np.float64, np.float32, np.complex64, np.complex128],
-                    [ift.HartleyOperator, ift.FFTOperator]))
-    def test_fft2D(self, dim1, dim2, d1, d2, itp, op):
-        tol = _get_rtol(itp)
-        a = ift.RGSpace([dim1, dim2], distances=[d1, d2])
-        b = ift.RGSpace([dim1, dim2],
-                        distances=[1./(dim1*d1), 1./(dim2*d2)], harmonic=True)
-
-        fft = op(domain=a, target=b)
-        inp = ift.Field.from_random(domain=a, random_type='normal',
-                                    std=7, mean=3, dtype=itp)
-        out = fft.inverse_times(fft.times(inp))
-        assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
-
-        a, b = b, a
-
-        fft = ift.FFTOperator(domain=a, target=b)
-        inp = ift.Field.from_random(domain=a, random_type='normal',
-                                    std=7, mean=3, dtype=itp)
-        out = fft.inverse_times(fft.times(inp))
-        assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
-
-    @expand(product([0, 1, 2],
-                    [np.float64, np.float32, np.complex64, np.complex128],
-                    [ift.HartleyOperator, ift.FFTOperator]))
-    def test_composed_fft(self, index, dtype, op):
-        tol = _get_rtol(dtype)
-        a = [a1, a2, a3] = [ift.RGSpace((32,)), ift.RGSpace((4, 4)),
-                            ift.RGSpace((5, 6))]
-        fft = op(domain=a, space=index)
-
-        inp = ift.Field.from_random(domain=(a1, a2, a3), random_type='normal',
-                                    std=7, mean=3, dtype=dtype)
-        out = fft.inverse_times(fft.times(inp))
-        assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
-
-    @expand(product([ift.RGSpace(128, distances=3.76, harmonic=True),
-                     ift.RGSpace((15, 27), distances=(.7, .33), harmonic=True),
-                     ift.RGSpace(73, distances=0.5643)],
-                    [np.float64, np.float32, np.complex64, np.complex128],
-                    [ift.HartleyOperator, ift.FFTOperator]))
-    def test_normalisation(self, space, tp, op):
-        tol = 10 * _get_rtol(tp)
-        cospace = space.get_default_codomain()
-        fft = op(space, cospace)
-        inp = ift.Field.from_random(domain=space, random_type='normal',
-                                    std=1, mean=2, dtype=tp)
-        out = fft.times(inp)
-        fft2 = op(cospace, space)
-        out2 = fft2.inverse_times(inp)
-        zero_idx = tuple([0]*len(space.shape))
-        assert_allclose(inp.to_global_data()[zero_idx], out.integrate(),
-                        rtol=tol, atol=tol)
-        assert_allclose(out.local_data, out2.local_data, rtol=tol, atol=tol)
+pmp = pytest.mark.parametrize
+dtype = list2fixture([np.float64, np.float32, np.complex64, np.complex128])
+op = list2fixture([ift.HartleyOperator, ift.FFTOperator])
+
+
+@pmp('d', [0.1, 1, 3.7])
+def test_fft1D(d, dtype, op):
+    dim1 = 16
+    tol = _get_rtol(dtype)
+    a = ift.RGSpace(dim1, distances=d)
+    b = ift.RGSpace(dim1, distances=1./(dim1*d), harmonic=True)
+    np.random.seed(16)
+
+    fft = op(domain=a, target=b)
+    inp = ift.Field.from_random(
+        domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
+    out = fft.inverse_times(fft.times(inp))
+    assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
+
+    a, b = b, a
+
+    fft = ift.FFTOperator(domain=a, target=b)
+    inp = ift.Field.from_random(
+        domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
+    out = fft.inverse_times(fft.times(inp))
+    assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
+
+
+@pmp('dim1', [12, 15])
+@pmp('dim2', [9, 12])
+@pmp('d1', [0.1, 1, 3.7])
+@pmp('d2', [0.4, 1, 2.7])
+def test_fft2D(dim1, dim2, d1, d2, dtype, op):
+    tol = _get_rtol(dtype)
+    a = ift.RGSpace([dim1, dim2], distances=[d1, d2])
+    b = ift.RGSpace(
+        [dim1, dim2], distances=[1./(dim1*d1), 1./(dim2*d2)], harmonic=True)
+
+    fft = op(domain=a, target=b)
+    inp = ift.Field.from_random(
+        domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
+    out = fft.inverse_times(fft.times(inp))
+    assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
+
+    a, b = b, a
+
+    fft = ift.FFTOperator(domain=a, target=b)
+    inp = ift.Field.from_random(
+        domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
+    out = fft.inverse_times(fft.times(inp))
+    assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
+
+
+@pmp('index', [0, 1, 2])
+def test_composed_fft(index, dtype, op):
+    tol = _get_rtol(dtype)
+    a = [a1, a2,
+         a3] = [ift.RGSpace((32,)),
+                ift.RGSpace((4, 4)),
+                ift.RGSpace((5, 6))]
+    fft = op(domain=a, space=index)
+
+    inp = ift.Field.from_random(
+        domain=(a1, a2, a3), random_type='normal', std=7, mean=3, dtype=dtype)
+    out = fft.inverse_times(fft.times(inp))
+    assert_allclose(inp.local_data, out.local_data, rtol=tol, atol=tol)
+
+
+@pmp('space', [
+    ift.RGSpace(128, distances=3.76, harmonic=True),
+    ift.RGSpace((15, 27), distances=(.7, .33), harmonic=True),
+    ift.RGSpace(73, distances=0.5643)
+])
+def test_normalisation(space, dtype, op):
+    tol = 10*_get_rtol(dtype)
+    cospace = space.get_default_codomain()
+    fft = op(space, cospace)
+    inp = ift.Field.from_random(
+        domain=space, random_type='normal', std=1, mean=2, dtype=dtype)
+    out = fft.times(inp)
+    fft2 = op(cospace, space)
+    out2 = fft2.inverse_times(inp)
+    zero_idx = tuple([0]*len(space.shape))
+    assert_allclose(
+        inp.to_global_data()[zero_idx], out.integrate(), rtol=tol, atol=tol)
+    assert_allclose(out.local_data, out2.local_data, rtol=tol, atol=tol)
diff --git a/test/test_operators/test_harmonic_transform_operator.py b/test/test_operators/test_harmonic_transform_operator.py
index e78475a270c96bf0af47a9f823733154d0184915..94aab8f537dad31fc7caf046c9245b269c74ac82 100644
--- a/test/test_operators/test_harmonic_transform_operator.py
+++ b/test/test_operators/test_harmonic_transform_operator.py
@@ -15,14 +15,14 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
 import numpy as np
+import pytest
 from numpy.testing import assert_allclose
 
+import nifty5 as ift
+
+from ..common import list2fixture
+
 
 def _get_rtol(tp):
     if (tp == np.float64) or (tp == np.complex128):
@@ -31,44 +31,45 @@ def _get_rtol(tp):
         return 1e-5
 
 
-class HarmonicTransformOperatorTests(unittest.TestCase):
-    @expand(product([128, 256],
-                    [np.float64, np.float32, np.complex64, np.complex128]))
-    def test_dotsht(self, lm, tp):
-        tol = 10 * _get_rtol(tp)
-        a = ift.LMSpace(lmax=lm)
-        b = ift.GLSpace(nlat=lm+1)
-        fft = ift.HarmonicTransformOperator(domain=a, target=b)
-        inp = ift.Field.from_random(domain=a, random_type='normal',
-                                    std=1, mean=0, dtype=tp)
-        out = fft.times(inp)
-        v1 = np.sqrt(out.vdot(out))
-        v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
-        assert_allclose(v1, v2, rtol=tol, atol=tol)
+tp = list2fixture([np.float64, np.float32, np.complex64, np.complex128])
+lm = list2fixture([128, 256])
+pmp = pytest.mark.parametrize
+
+
+def test_dotsht(lm, tp):
+    tol = 10*_get_rtol(tp)
+    a = ift.LMSpace(lmax=lm)
+    b = ift.GLSpace(nlat=lm + 1)
+    fft = ift.HarmonicTransformOperator(domain=a, target=b)
+    inp = ift.Field.from_random(
+        domain=a, random_type='normal', std=1, mean=0, dtype=tp)
+    out = fft.times(inp)
+    v1 = np.sqrt(out.vdot(out))
+    v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
+    assert_allclose(v1, v2, rtol=tol, atol=tol)
+
+
+def test_dotsht2(lm, tp):
+    tol = 10*_get_rtol(tp)
+    a = ift.LMSpace(lmax=lm)
+    b = ift.HPSpace(nside=lm//2)
+    fft = ift.HarmonicTransformOperator(domain=a, target=b)
+    inp = ift.Field.from_random(
+        domain=a, random_type='normal', std=1, mean=0, dtype=tp)
+    out = fft.times(inp)
+    v1 = np.sqrt(out.vdot(out))
+    v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
+    assert_allclose(v1, v2, rtol=tol, atol=tol)
 
-    @expand(product([128, 256],
-                    [np.float64, np.float32, np.complex64, np.complex128]))
-    def test_dotsht2(self, lm, tp):
-        tol = 10 * _get_rtol(tp)
-        a = ift.LMSpace(lmax=lm)
-        b = ift.HPSpace(nside=lm//2)
-        fft = ift.HarmonicTransformOperator(domain=a, target=b)
-        inp = ift.Field.from_random(domain=a, random_type='normal',
-                                    std=1, mean=0, dtype=tp)
-        out = fft.times(inp)
-        v1 = np.sqrt(out.vdot(out))
-        v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
-        assert_allclose(v1, v2, rtol=tol, atol=tol)
 
-    @expand(product([ift.LMSpace(lmax=30, mmax=25)],
-                    [np.float64, np.float32, np.complex64, np.complex128]))
-    def test_normalisation(self, space, tp):
-        tol = 10 * _get_rtol(tp)
-        cospace = space.get_default_codomain()
-        fft = ift.HarmonicTransformOperator(space, cospace)
-        inp = ift.Field.from_random(domain=space, random_type='normal',
-                                    std=1, mean=2, dtype=tp)
-        out = fft.times(inp)
-        zero_idx = tuple([0]*len(space.shape))
-        assert_allclose(inp.to_global_data()[zero_idx], out.integrate(),
-                        rtol=tol, atol=tol)
+@pmp('space', [ift.LMSpace(lmax=30, mmax=25)])
+def test_normalisation(space, tp):
+    tol = 10*_get_rtol(tp)
+    cospace = space.get_default_codomain()
+    fft = ift.HarmonicTransformOperator(space, cospace)
+    inp = ift.Field.from_random(
+        domain=space, random_type='normal', std=1, mean=2, dtype=tp)
+    out = fft.times(inp)
+    zero_idx = tuple([0]*len(space.shape))
+    assert_allclose(
+        inp.to_global_data()[zero_idx], out.integrate(), rtol=tol, atol=tol)
diff --git a/test/test_operators/test_operator_gradients.py b/test/test_operators/test_operator_gradients.py
deleted file mode 100644
index 6bf047a4c67d7283176b1bb5417990515ea5b25f..0000000000000000000000000000000000000000
--- a/test/test_operators/test_operator_gradients.py
+++ /dev/null
@@ -1,131 +0,0 @@
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-# Copyright(C) 2013-2019 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
-import numpy as np
-
-
-class OperatorTests(unittest.TestCase):
-    @staticmethod
-    def make_linearization(type, space, seed):
-        np.random.seed(seed)
-        S = ift.ScalingOperator(1., space)
-        s = S.draw_sample()
-        if type == "Constant":
-            return ift.Linearization.make_const(s)
-        elif type == "Variable":
-            return ift.Linearization.make_var(s)
-        raise ValueError('unknown type passed')
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testBasics(self, space, seed):
-        var = self.make_linearization("Variable", space, seed)
-        op = ift.ScalingOperator(6., var.target)
-        ift.extra.check_value_gradient_consistency(op, var.val)
-
-    @expand(product(
-        ['Variable', 'Constant'],
-        ['Variable'],
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]
-        ))
-    def testBinary(self, type1, type2, space, seed):
-        dom1 = ift.MultiDomain.make({'s1': space})
-        lin1 = self.make_linearization(type1, dom1, seed)
-        dom2 = ift.MultiDomain.make({'s2': space})
-        lin2 = self.make_linearization(type2, dom2, seed)
-
-        dom = ift.MultiDomain.union((dom1, dom2))
-        select_s1 = ift.ducktape(None, dom, "s1")
-        select_s2 = ift.ducktape(None, dom, "s2")
-        op = select_s1*select_s2
-        pos = ift.from_random("normal", dom)
-        ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
-        op = select_s1+select_s2
-        pos = ift.from_random("normal", dom)
-        ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
-        op = select_s1.scale(3.)
-        pos = ift.from_random("normal", dom1)
-        ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
-        op = ift.ScalingOperator(2.456, space)(select_s1*select_s2)
-        pos = ift.from_random("normal", dom)
-        ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
-        op = ift.sigmoid(ift.ScalingOperator(2.456, space)(
-            select_s1*select_s2))
-        pos = ift.from_random("normal", dom)
-        ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
-        pos = ift.from_random("normal", dom)
-        op = ift.OuterProduct(pos['s1'], ift.makeDomain(space))
-        ift.extra.check_value_gradient_consistency(op, pos['s2'], ntries=20)
-        if isinstance(space, ift.RGSpace):
-            op = ift.FFTOperator(space)(select_s1*select_s2)
-            pos = ift.from_random("normal", dom)
-            ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4],
-        [0.5],
-        [2.],
-        [3.],
-        [1.5],
-        [1.75],
-        [1.3],
-        [4, 78, 23],
-        ))
-    def testOperatorLibrary(self, space, Npixdof, ceps_a,
-                            ceps_k, sm, sv, im, iv, seed):
-        # tests amplitude operator and coorelated field operator
-        np.random.seed(seed)
-        op = ift.AmplitudeOperator(space, Npixdof, ceps_a, ceps_k, sm,
-                                   sv, im, iv)
-        S = ift.ScalingOperator(1., op.domain)
-        pos = S.draw_sample()
-        ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
-
-        op2 = ift.CorrelatedField(space, op)
-        S = ift.ScalingOperator(1., op2.domain)
-        pos = S.draw_sample()
-        ift.extra.check_value_gradient_consistency(op2, pos, ntries=20)
-
-    @expand(product(
-        [ift.GLSpace(15),
-         ift.RGSpace(64, distances=.789),
-         ift.RGSpace([32, 32], distances=.789)],
-        [4, 78, 23]))
-    def testInvGammaOperator(self, space, seed):
-        S = ift.ScalingOperator(1., space)
-        pos = S.draw_sample()
-        alpha = 1.5
-        q = 0.73
-        op = ift.InverseGammaOperator(space, alpha, q)
-        # FIXME All those cdfs and ppfs are not very accurate
-        ift.extra.check_value_gradient_consistency(op, pos, tol=1e-2,
-                                                   ntries=20)
diff --git a/test/test_operators/test_regridding.py b/test/test_operators/test_regridding.py
index 318486a93693824e380496987bf26e40c88217b1..9f804a6221b77e78bd9d21bf9776158f8275e598 100644
--- a/test/test_operators/test_regridding.py
+++ b/test/test_operators/test_regridding.py
@@ -15,23 +15,20 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
-
 from numpy.testing import assert_allclose
 
 import nifty5 as ift
 
+from ..common import list2fixture
+
+s = list2fixture([
+    ift.RGSpace(8, distances=12.9),
+    ift.RGSpace(59, distances=.24, harmonic=True),
+    ift.RGSpace([12, 3])
+])
+
 
-class Regridding_Tests(unittest.TestCase):
-    @expand(
-        product([
-            ift.RGSpace(8, distances=12.9),
-            ift.RGSpace(59, distances=.24, harmonic=True),
-            ift.RGSpace([12, 3])
-        ]))
-    def test_value(self, s):
-        Regrid = ift.RegriddingOperator(s, s.shape)
-        f = ift.from_random('normal', Regrid.domain)
-        assert_allclose(f.to_global_data(), Regrid(f).to_global_data())
+def test_value(s):
+    Regrid = ift.RegriddingOperator(s, s.shape)
+    f = ift.from_random('normal', Regrid.domain)
+    assert_allclose(f.to_global_data(), Regrid(f).to_global_data())
diff --git a/test/test_operators/test_smoothing_operator.py b/test/test_operators/test_smoothing_operator.py
index 4b160a34bed609f04a515e4909ea984ee964007a..73c7adfdbe7514d9cb5a3d3583fd55e0ff0daf6d 100644
--- a/test/test_operators/test_smoothing_operator.py
+++ b/test/test_operators/test_smoothing_operator.py
@@ -15,14 +15,14 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
-
-import nifty5 as ift
 import numpy as np
+import pytest
 from numpy.testing import assert_allclose
 
+import nifty5 as ift
+
+from ..common import list2fixture
+
 
 def _get_rtol(tp):
     if (tp == np.float64) or (tp == np.complex128):
@@ -31,51 +31,57 @@ def _get_rtol(tp):
         return 1e-5
 
 
-class SmoothingOperator_Tests(unittest.TestCase):
-    spaces = [ift.RGSpace(128)]
-
-    @expand(product(spaces, [0., .5, 5.]))
-    def test_property(self, space, sigma):
-        op = ift.HarmonicSmoothingOperator(space, sigma=sigma)
-        if op.domain[0] != space:
-            raise TypeError
-
-    @expand(product(spaces, [0., .5, 5.]))
-    def test_adjoint_times(self, space, sigma):
-        op = ift.HarmonicSmoothingOperator(space, sigma=sigma)
-        rand1 = ift.Field.from_random('normal', domain=space)
-        rand2 = ift.Field.from_random('normal', domain=space)
-        tt1 = rand1.vdot(op.times(rand2))
-        tt2 = rand2.vdot(op.adjoint_times(rand1))
-        assert_allclose(tt1, tt2)
-
-    @expand(product(spaces, [0., .5, 5.]))
-    def test_times(self, space, sigma):
-        op = ift.HarmonicSmoothingOperator(space, sigma=sigma)
-        fld = np.zeros(space.shape, dtype=np.float64)
-        fld[0] = 1.
-        rand1 = ift.Field.from_global_data(space, fld)
-        tt1 = op.times(rand1)
-        assert_allclose(1, tt1.sum())
-
-    @expand(product([128, 256], [1, 0.4], [0., 1.,  3.7],
-                    [np.float64, np.complex128]))
-    def test_smooth_regular1(self, sz, d, sigma, tp):
-        tol = _get_rtol(tp)
-        sp = ift.RGSpace(sz, distances=d)
-        smo = ift.HarmonicSmoothingOperator(sp, sigma=sigma)
-        inp = ift.Field.from_random(domain=sp, random_type='normal', std=1,
-                                    mean=4, dtype=tp)
-        out = smo(inp)
-        assert_allclose(inp.sum(), out.sum(), rtol=tol, atol=tol)
-
-    @expand(product([10, 15], [7, 10], [1, 0.4], [2, 0.3], [0., 1.,  3.7],
-                    [np.float64, np.complex128]))
-    def test_smooth_regular2(self, sz1, sz2, d1, d2, sigma, tp):
-        tol = _get_rtol(tp)
-        sp = ift.RGSpace([sz1, sz2], distances=[d1, d2])
-        smo = ift.HarmonicSmoothingOperator(sp, sigma=sigma)
-        inp = ift.Field.from_random(domain=sp, random_type='normal', std=1,
-                                    mean=4, dtype=tp)
-        out = smo(inp)
-        assert_allclose(inp.sum(), out.sum(), rtol=tol, atol=tol)
+space = list2fixture([ift.RGSpace(128)])
+sigma = list2fixture([0., .5, 5.])
+tp = list2fixture([np.float64, np.complex128])
+pmp = pytest.mark.parametrize
+
+
+def test_property(space, sigma):
+    op = ift.HarmonicSmoothingOperator(space, sigma=sigma)
+    if op.domain[0] != space:
+        raise TypeError
+
+
+def test_adjoint_times(space, sigma):
+    op = ift.HarmonicSmoothingOperator(space, sigma=sigma)
+    rand1 = ift.Field.from_random('normal', domain=space)
+    rand2 = ift.Field.from_random('normal', domain=space)
+    tt1 = rand1.vdot(op.times(rand2))
+    tt2 = rand2.vdot(op.adjoint_times(rand1))
+    assert_allclose(tt1, tt2)
+
+
+def test_times(space, sigma):
+    op = ift.HarmonicSmoothingOperator(space, sigma=sigma)
+    fld = np.zeros(space.shape, dtype=np.float64)
+    fld[0] = 1.
+    rand1 = ift.Field.from_global_data(space, fld)
+    tt1 = op.times(rand1)
+    assert_allclose(1, tt1.sum())
+
+
+@pmp('sz', [128, 256])
+@pmp('d', [1, 0.4])
+def test_smooth_regular1(sz, d, sigma, tp):
+    tol = _get_rtol(tp)
+    sp = ift.RGSpace(sz, distances=d)
+    smo = ift.HarmonicSmoothingOperator(sp, sigma=sigma)
+    inp = ift.Field.from_random(
+        domain=sp, random_type='normal', std=1, mean=4, dtype=tp)
+    out = smo(inp)
+    assert_allclose(inp.sum(), out.sum(), rtol=tol, atol=tol)
+
+
+@pmp('sz1', [10, 15])
+@pmp('sz2', [7, 10])
+@pmp('d1', [1, 0.4])
+@pmp('d2', [2, 0.3])
+def test_smooth_regular2(sz1, sz2, d1, d2, sigma, tp):
+    tol = _get_rtol(tp)
+    sp = ift.RGSpace([sz1, sz2], distances=[d1, d2])
+    smo = ift.HarmonicSmoothingOperator(sp, sigma=sigma)
+    inp = ift.Field.from_random(
+        domain=sp, random_type='normal', std=1, mean=4, dtype=tp)
+    out = smo(inp)
+    assert_allclose(inp.sum(), out.sum(), rtol=tol, atol=tol)
diff --git a/test/test_spaces/__init__.py b/test/test_spaces/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/test/test_spaces/test_gl_space.py b/test/test_spaces/test_gl_space.py
index 03ef8f3204232ebd231a53757e8a442858c8bacc..ab35c5d4fe3bd8bbf25487fa7e4985c913cd7614 100644
--- a/test/test_spaces/test_gl_space.py
+++ b/test/test_spaces/test_gl_space.py
@@ -16,61 +16,62 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import itertools
-import unittest
-from test.common import expand
 
 import numpy as np
-from nifty5 import GLSpace
+import pytest
 from numpy.testing import (assert_, assert_almost_equal, assert_equal,
                            assert_raises)
 
+from nifty5 import GLSpace
+
+pmp = pytest.mark.parametrize
+
 # [nlat, nlon, expected]
-CONSTRUCTOR_CONFIGS = [
-        [2, None, {
-            'nlat': 2,
-            'nlon': 3,
-            'harmonic': False,
-            'shape': (6,),
-            'size': 6,
-            }],
-        [0, None, {
-            'error': ValueError
-            }]
-    ]
+CONSTRUCTOR_CONFIGS = [[
+    2, None, {
+        'nlat': 2,
+        'nlon': 3,
+        'harmonic': False,
+        'shape': (6,),
+        'size': 6,
+    }
+], [0, None, {
+    'error': ValueError
+}]]
 
 
 def get_dvol_configs():
     np.random.seed(42)
-    wgt = [2.0943951,  2.0943951]
+    wgt = [2.0943951, 2.0943951]
     # for GLSpace(nlat=2, nlon=3)
-    dvol_0 = np.array(list(itertools.chain.from_iterable(
-        itertools.repeat(x, 3) for x in wgt)))
+    dvol_0 = np.array(
+        list(
+            itertools.chain.from_iterable(
+                itertools.repeat(x, 3) for x in wgt)))
     return [
         [1, dvol_0],
-        ]
+    ]
+
 
+@pmp('attribute', ['nlat', 'nlon'])
+def test_property_ret_type(attribute):
+    g = GLSpace(2)
+    assert_(isinstance(getattr(g, attribute), int))
 
-class GLSpaceInterfaceTests(unittest.TestCase):
-    @expand([['nlat', int],
-            ['nlon', int]])
-    def test_property_ret_type(self, attribute, expected_type):
-        g = GLSpace(2)
-        assert_(isinstance(getattr(g, attribute), expected_type))
 
+@pmp('nlat, nlon, expected', CONSTRUCTOR_CONFIGS)
+def test_constructor(nlat, nlon, expected):
+    g = GLSpace(4)
 
-class GLSpaceFunctionalityTests(unittest.TestCase):
-    @expand(CONSTRUCTOR_CONFIGS)
-    def test_constructor(self, nlat, nlon, expected):
-        g = GLSpace(4)
+    if 'error' in expected:
+        with assert_raises(expected['error']):
+            GLSpace(nlat, nlon)
+    else:
+        g = GLSpace(nlat, nlon)
+        for key, value in expected.items():
+            assert_equal(getattr(g, key), value)
 
-        if 'error' in expected:
-            with assert_raises(expected['error']):
-                GLSpace(nlat, nlon)
-        else:
-            g = GLSpace(nlat, nlon)
-            for key, value in expected.items():
-                assert_equal(getattr(g, key), value)
 
-    @expand(get_dvol_configs())
-    def test_dvol(self, power, expected):
-        assert_almost_equal(GLSpace(2).dvol, expected)
+@pmp('power, expected', get_dvol_configs())
+def test_dvol(power, expected):
+    assert_almost_equal(GLSpace(2).dvol, expected)
diff --git a/test/test_spaces/test_hp_space.py b/test/test_spaces/test_hp_space.py
index ce8412ea571a93ecb059d238e24b25c0874e6faf..9da92be51cd4f070c1cc862176355e2b065f4700 100644
--- a/test/test_spaces/test_hp_space.py
+++ b/test/test_spaces/test_hp_space.py
@@ -15,57 +15,52 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from test.common import expand
-
 import numpy as np
-from nifty5 import HPSpace
+import pytest
 from numpy.testing import (assert_, assert_almost_equal, assert_equal,
                            assert_raises)
 
+from nifty5 import HPSpace
+
+pmp = pytest.mark.parametrize
 # [nside, expected]
-CONSTRUCTOR_CONFIGS = [
-        [2, {
-            'nside': 2,
-            'harmonic': False,
-            'shape': (48,),
-            'size': 48,
-            }],
-        [5, {
-            'nside': 5,
-            'harmonic': False,
-            'shape': (300,),
-            'size': 300,
-            }],
-        [1, {
-            'nside': 1,
-            'harmonic': False,
-            'shape': (12,),
-            'size': 12,
-            }],
-        [0, {
-            'error': ValueError
-            }]
-    ]
+CONSTRUCTOR_CONFIGS = [[
+    2, {
+        'nside': 2,
+        'harmonic': False,
+        'shape': (48,),
+        'size': 48,
+    }
+], [5, {
+    'nside': 5,
+    'harmonic': False,
+    'shape': (300,),
+    'size': 300,
+}], [1, {
+    'nside': 1,
+    'harmonic': False,
+    'shape': (12,),
+    'size': 12,
+}], [0, {
+    'error': ValueError
+}]]
+
 
+def test_property_ret_type():
+    x = HPSpace(2)
+    assert_(isinstance(getattr(x, 'nside'), int))
 
-class HPSpaceInterfaceTests(unittest.TestCase):
-    @expand([['nside', int]])
-    def test_property_ret_type(self, attribute, expected_type):
-        x = HPSpace(2)
-        assert_(isinstance(getattr(x, attribute), expected_type))
 
+@pmp('nside, expected', CONSTRUCTOR_CONFIGS)
+def test_constructor(nside, expected):
+    if 'error' in expected:
+        with assert_raises(expected['error']):
+            HPSpace(nside)
+    else:
+        h = HPSpace(nside)
+        for key, value in expected.items():
+            assert_equal(getattr(h, key), value)
 
-class HPSpaceFunctionalityTests(unittest.TestCase):
-    @expand(CONSTRUCTOR_CONFIGS)
-    def test_constructor(self, nside, expected):
-        if 'error' in expected:
-            with assert_raises(expected['error']):
-                HPSpace(nside)
-        else:
-            h = HPSpace(nside)
-            for key, value in expected.items():
-                assert_equal(getattr(h, key), value)
 
-    def test_dvol(self):
-        assert_almost_equal(HPSpace(2).dvol, np.pi/12)
+def test_dvol():
+    assert_almost_equal(HPSpace(2).dvol, np.pi/12)
diff --git a/test/test_spaces/test_interface.py b/test/test_spaces/test_interface.py
index fda9848f3a57d7e733c1d8f2286c6b87942f7a73..9ce9129e70e39803f733af3541cb4ead99afd165 100644
--- a/test/test_spaces/test_interface.py
+++ b/test/test_spaces/test_interface.py
@@ -15,33 +15,37 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from itertools import product
-from test.common import expand
 from types import LambdaType
 
-import nifty5 as ift
+import pytest
 from numpy.testing import assert_
 
-spaces = [ift.RGSpace(4),
-          ift.PowerSpace(ift.RGSpace((4, 4), harmonic=True)),
-          ift.LMSpace(5), ift.HPSpace(4), ift.GLSpace(4)]
-
-
-class SpaceInterfaceTests(unittest.TestCase):
-    @expand(product(spaces, [
-                    ['harmonic', bool],
-                    ['shape', tuple],
-                    ['size', int]]))
-    def test_property_ret_type(self, space, attr_expected_type):
-        assert_(isinstance(getattr(space, attr_expected_type[0]),
-                           attr_expected_type[1]))
-
-    @expand(product([ift.RGSpace(4, harmonic=True), ift.LMSpace(5)], [
-        ['get_k_length_array', ift.Field],
-        ['get_fft_smoothing_kernel_function', 2.0, LambdaType],
-        ]))
-    def test_method_ret_type(self, space, method_expected_type):
-        assert_(type(getattr(space, method_expected_type[0])(
-                          *method_expected_type[1:-1])) is
-                method_expected_type[-1])
+import nifty5 as ift
+
+pmp = pytest.mark.parametrize
+
+
+@pmp('attr_expected_type',
+     [['harmonic', bool], ['shape', tuple], ['size', int]])
+@pmp('space', [
+    ift.RGSpace(4),
+    ift.PowerSpace(ift.RGSpace((4, 4), harmonic=True)),
+    ift.LMSpace(5),
+    ift.HPSpace(4),
+    ift.GLSpace(4)
+])
+def test_property_ret_type(space, attr_expected_type):
+    assert_(
+        isinstance(
+            getattr(space, attr_expected_type[0]), attr_expected_type[1]))
+
+
+@pmp('method_expected_type',
+     [['get_k_length_array', ift.Field],
+      ['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_(
+        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 bb9b217c1f1e200d3c4d28ca4f00d1e75b66c4ba..9568fbfa15704085e3d7294cd3c0130b7d97b119 100644
--- a/test/test_spaces/test_lm_space.py
+++ b/test/test_spaces/test_lm_space.py
@@ -15,33 +15,33 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from test.common import expand
-
-import nifty5 as ift
 import numpy as np
+import pytest
 from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
 
+import nifty5 as ift
+
+pmp = pytest.mark.parametrize
 # [lmax, expected]
-CONSTRUCTOR_CONFIGS = [
-        [5, None, {
-            'lmax': 5,
-            'mmax': 5,
-            'shape': (36,),
-            'harmonic': True,
-            'size': 36,
-            }],
-        [7, 4, {
-            'lmax': 7,
-            'mmax': 4,
-            'shape': (52,),
-            'harmonic': True,
-            'size': 52,
-            }],
-        [-1, 28, {
-            'error': ValueError
-            }]
-    ]
+CONSTRUCTOR_CONFIGS = [[
+    5, None, {
+        'lmax': 5,
+        'mmax': 5,
+        'shape': (36,),
+        'harmonic': True,
+        'size': 36,
+    }
+], [
+    7, 4, {
+        'lmax': 7,
+        'mmax': 4,
+        'shape': (52,),
+        'harmonic': True,
+        'size': 52,
+    }
+], [-1, 28, {
+    'error': ValueError
+}]]
 
 
 def _k_length_array_helper(index_arr, lmax):
@@ -53,8 +53,8 @@ def _k_length_array_helper(index_arr, lmax):
         else:
             index_half = (index_arr + lmax + 1)//2
 
-    m = np.ceil(((2*lmax + 1) - np.sqrt((2*lmax + 1)**2 -
-                 8*(index_half - lmax)))/2).astype(int)
+    m = np.ceil(((2*lmax + 1) - np.sqrt((2*lmax + 1)**2 - 8 *
+                                        (index_half - lmax)))/2).astype(int)
 
     return index_half - m*(2*lmax + 1 - m)//2
 
@@ -64,30 +64,28 @@ def get_k_length_array_configs():
     return [[5, da_0]]
 
 
-class LMSpaceInterfaceTests(unittest.TestCase):
-    @expand([['lmax', int],
-            ['mmax', int],
-            ['size', int]])
-    def test_property_ret_type(self, attribute, expected_type):
-        l = ift.LMSpace(7, 5)
-        assert_(isinstance(getattr(l, attribute), expected_type))
+@pmp('attribute', ['lmax', 'mmax', 'size'])
+def test_property_ret_type(attribute):
+    l = ift.LMSpace(7, 5)
+    assert_(isinstance(getattr(l, attribute), int))
 
 
-class LMSpaceFunctionalityTests(unittest.TestCase):
-    @expand(CONSTRUCTOR_CONFIGS)
-    def test_constructor(self, lmax, mmax, expected):
-        if 'error' in expected:
-            with assert_raises(expected['error']):
-                ift.LMSpace(lmax, mmax)
-        else:
-            l = ift.LMSpace(lmax, mmax)
-            for key, value in expected.items():
-                assert_equal(getattr(l, key), value)
+@pmp('lmax, mmax, expected', CONSTRUCTOR_CONFIGS)
+def test_constructor(lmax, mmax, expected):
+    if 'error' in expected:
+        with assert_raises(expected['error']):
+            ift.LMSpace(lmax, mmax)
+    else:
+        l = ift.LMSpace(lmax, mmax)
+        for key, value in expected.items():
+            assert_equal(getattr(l, key), value)
+
+
+def test_dvol():
+    assert_allclose(ift.LMSpace(5).dvol, 1.)
 
-    def test_dvol(self):
-        assert_allclose(ift.LMSpace(5).dvol, 1.)
 
-    @expand(get_k_length_array_configs())
-    def test_k_length_array(self, lmax, expected):
-        l = ift.LMSpace(lmax)
-        assert_allclose(l.get_k_length_array().to_global_data(), expected)
+@pmp('lmax, expected', get_k_length_array_configs())
+def test_k_length_array(lmax, expected):
+    l = ift.LMSpace(lmax)
+    assert_allclose(l.get_k_length_array().to_global_data(), expected)
diff --git a/test/test_spaces/test_power_space.py b/test/test_spaces/test_power_space.py
index db5d3f6d11e63c28212e9c5d4f378477b92f2eb9..55e33ad13bf2dbb0a8f042ad3e80e9ae4e76364a 100644
--- a/test/test_spaces/test_power_space.py
+++ b/test/test_spaces/test_power_space.py
@@ -15,119 +15,137 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
 from itertools import chain, product
-from test.common import expand
 
+import pytest
 import nifty5 as ift
 import numpy as np
 from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
 
-HARMONIC_SPACES = [ift.RGSpace((8,), harmonic=True),
-                   ift.RGSpace((7, 8), harmonic=True),
-                   ift.RGSpace((6, 6), harmonic=True),
-                   ift.RGSpace((5, 5), harmonic=True),
-                   ift.RGSpace((4, 5, 7), harmonic=True),
-                   ift.LMSpace(6),
-                   ift.LMSpace(9)]
+pmp = pytest.mark.parametrize
 
+HARMONIC_SPACES = [
+    ift.RGSpace((8,), harmonic=True),
+    ift.RGSpace((7, 8), harmonic=True),
+    ift.RGSpace((6, 6), harmonic=True),
+    ift.RGSpace((5, 5), harmonic=True),
+    ift.RGSpace((4, 5, 7), harmonic=True),
+    ift.LMSpace(6),
+    ift.LMSpace(9)
+]
 
 # Try all sensible kinds of combinations of spaces and binning parameters
-CONSISTENCY_CONFIGS_IMPLICIT = product(HARMONIC_SPACES,
-                                       [None], [None, 3, 4], [True, False])
-CONSISTENCY_CONFIGS_EXPLICIT = product(HARMONIC_SPACES,
-                                       [[0., 1.3]], [None], [None])
+CONSISTENCY_CONFIGS_IMPLICIT = product(HARMONIC_SPACES, [None], [None, 3, 4],
+                                       [True, False])
+CONSISTENCY_CONFIGS_EXPLICIT = product(HARMONIC_SPACES, [[0., 1.3]], [None],
+                                       [None])
 CONSISTENCY_CONFIGS = chain(CONSISTENCY_CONFIGS_IMPLICIT,
                             CONSISTENCY_CONFIGS_EXPLICIT)
 
 # [harmonic_partner, logarithmic, nbin, binbounds, expected]
 CONSTRUCTOR_CONFIGS = [
-    [1, False, None, None, {'error': (ValueError, NotImplementedError)}],
-    [ift.RGSpace((8,)), False, None, None, {'error': ValueError}],
-    [ift.RGSpace((8,), harmonic=True), None, None, None, {
-        'harmonic': False,
-        'shape': (5,),
-        'size': 5,
-        'harmonic_partner': ift.RGSpace((8,), harmonic=True),
-        'binbounds': None,
-        'pindex': ift.dobj.from_global_data(
-            np.array([0, 1, 2, 3, 4, 3, 2, 1])),
-        'k_lengths': np.array([0., 1., 2., 3., 4.]),
-        }],
-    [ift.RGSpace((8,), harmonic=True), True, None, None, {
-        'harmonic': False,
-        'shape': (4,),
-        'size': 4,
-        'harmonic_partner': ift.RGSpace((8,), harmonic=True),
-        'binbounds': (0.5, 1.3228756555322954, 3.5),
-        'pindex': ift.dobj.from_global_data(
-            np.array([0, 1, 2, 2, 3, 2, 2, 1])),
-        'k_lengths': np.array([0., 1., 2.5, 4.]),
-        }],
-    ]
+    [1, False, None, None, {
+        'error': (ValueError, NotImplementedError)
+    }],
+    [ift.RGSpace((8,)), False, None, None, {
+        'error': ValueError
+    }],
+    [
+        ift.RGSpace((8,), harmonic=True), None, None, None, {
+            'harmonic':
+            False,
+            'shape': (5,),
+            'size':
+            5,
+            'harmonic_partner':
+            ift.RGSpace((8,), harmonic=True),
+            'binbounds':
+            None,
+            'pindex':
+            ift.dobj.from_global_data(np.array([0, 1, 2, 3, 4, 3, 2, 1])),
+            'k_lengths':
+            np.array([0., 1., 2., 3., 4.]),
+        }
+    ],
+    [
+        ift.RGSpace((8,), harmonic=True), True, None, None, {
+            'harmonic':
+            False,
+            'shape': (4,),
+            'size':
+            4,
+            'harmonic_partner':
+            ift.RGSpace((8,), harmonic=True),
+            'binbounds': (0.5, 1.3228756555322954, 3.5),
+            'pindex':
+            ift.dobj.from_global_data(np.array([0, 1, 2, 2, 3, 2, 2, 1])),
+            'k_lengths':
+            np.array([0., 1., 2.5, 4.]),
+        }
+    ],
+]
 
 
 def k_lengths_configs():
     da_0 = np.array([0, 1.0, 1.41421356, 2., 2.23606798, 2.82842712])
     return [
-        [ift.RGSpace((4, 4), harmonic=True),  da_0],
-        ]
-
-
-class PowerSpaceInterfaceTest(unittest.TestCase):
-    @expand([
-        ['harmonic_partner', ift.StructuredDomain],
-        ['binbounds', type(None)],
-        ['pindex', ift.dobj.data_object],
-        ['k_lengths', np.ndarray],
-        ])
-    def test_property_ret_type(self, attribute, expected_type):
-        r = ift.RGSpace((4, 4), harmonic=True)
-        p = ift.PowerSpace(r)
-        assert_(isinstance(getattr(p, attribute), expected_type))
-
-
-class PowerSpaceConsistencyCheck(unittest.TestCase):
-    @expand(CONSISTENCY_CONFIGS)
-    def test_rhopindexConsistency(self, harmonic_partner,
-                                  binbounds, nbin, logarithmic):
+        [ift.RGSpace((4, 4), harmonic=True), da_0],
+    ]
+
+
+@pmp('attribute, expected_type', [
+    ['harmonic_partner', ift.StructuredDomain],
+    ['binbounds', type(None)],
+    ['pindex', ift.dobj.data_object],
+    ['k_lengths', np.ndarray],
+])
+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))
+
+
+@pmp('harmonic_partner, binbounds, nbin, logarithmic', CONSISTENCY_CONFIGS)
+def test_rhopindexConsistency(harmonic_partner, binbounds, nbin, logarithmic):
+    bb = ift.PowerSpace.useful_binbounds(harmonic_partner, logarithmic, nbin)
+    p = ift.PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb)
+
+    assert_equal(
+        np.bincount(ift.dobj.to_global_data(p.pindex).ravel()),
+        p.dvol,
+        err_msg='rho is not equal to pindex degeneracy')
+
+
+@pmp('harmonic_partner, logarithmic, nbin, binbounds, expected',
+     CONSTRUCTOR_CONFIGS)
+def test_constructor(harmonic_partner, logarithmic, nbin, binbounds, expected):
+    if 'error' in expected:
+        with assert_raises(expected['error']):
+            bb = ift.PowerSpace.useful_binbounds(harmonic_partner, logarithmic,
+                                                 nbin)
+            ift.PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb)
+    else:
         bb = ift.PowerSpace.useful_binbounds(harmonic_partner, logarithmic,
                                              nbin)
         p = ift.PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb)
-
-        assert_equal(np.bincount(ift.dobj.to_global_data(p.pindex).ravel()),
-                     p.dvol, err_msg='rho is not equal to pindex degeneracy')
-
-
-class PowerSpaceFunctionalityTest(unittest.TestCase):
-    @expand(CONSTRUCTOR_CONFIGS)
-    def test_constructor(self, harmonic_partner,
-                         logarithmic, nbin, binbounds, expected):
-        if 'error' in expected:
-            with assert_raises(expected['error']):
-                bb = ift.PowerSpace.useful_binbounds(harmonic_partner,
-                                                     logarithmic, nbin)
-                ift.PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb)
-        else:
-            bb = ift.PowerSpace.useful_binbounds(harmonic_partner,
-                                                 logarithmic, nbin)
-            p = ift.PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb)
-            for key, value in expected.items():
-                if isinstance(value, np.ndarray):
-                    assert_allclose(getattr(p, key), value)
-                else:
-                    assert_equal(getattr(p, key), value)
-
-    @expand(k_lengths_configs())
-    def test_k_lengths(self, harmonic_partner, expected):
-        p = ift.PowerSpace(harmonic_partner=harmonic_partner)
-        assert_allclose(p.k_lengths, expected)
-
-    def test_dvol(self):
-        hp = ift.RGSpace(10, harmonic=True)
-        p = ift.PowerSpace(harmonic_partner=hp)
-        v1 = hp.dvol
-        v1 = hp.size*v1 if np.isscalar(v1) else np.sum(v1)
-        v2 = p.dvol
-        v2 = p.size*v2 if np.isscalar(v2) else np.sum(v2)
-        assert_allclose(v1, v2)
+        for key, value in expected.items():
+            if isinstance(value, np.ndarray):
+                assert_allclose(getattr(p, key), value)
+            else:
+                assert_equal(getattr(p, key), value)
+
+
+@pmp('harmonic_partner, expected', k_lengths_configs())
+def test_k_lengths(harmonic_partner, expected):
+    p = ift.PowerSpace(harmonic_partner=harmonic_partner)
+    assert_allclose(p.k_lengths, expected)
+
+
+def test_dvol():
+    hp = ift.RGSpace(10, harmonic=True)
+    p = ift.PowerSpace(harmonic_partner=hp)
+    v1 = hp.dvol
+    v1 = hp.size*v1 if np.isscalar(v1) else np.sum(v1)
+    v2 = p.dvol
+    v2 = p.size*v2 if np.isscalar(v2) else np.sum(v2)
+    assert_allclose(v1, v2)
diff --git a/test/test_spaces/test_rg_space.py b/test/test_spaces/test_rg_space.py
index 3475d5902e748b40d0f67a0a6d23287724bcc6c6..5a73485fd6a926fad0502a3f41c6d9ee07b31ea3 100644
--- a/test/test_spaces/test_rg_space.py
+++ b/test/test_spaces/test_rg_space.py
@@ -15,103 +15,85 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-import unittest
-from test.common import expand
-
-import nifty5 as ift
 import numpy as np
+import pytest
 from numpy.testing import assert_, assert_allclose, assert_equal
 
-# [shape, distances, harmonic, expected]
-CONSTRUCTOR_CONFIGS = [
-        [(8,), None, False,
-            {
-                'shape': (8,),
-                'distances': (0.125,),
-                'harmonic': False,
-                'size': 8,
-            }],
-        [(8,), None, True,
-            {
-                'shape': (8,),
-                'distances': (1.0,),
-                'harmonic': True,
-                'size': 8,
-            }],
-        [(8,), (12,), True,
-            {
-                'shape': (8,),
-                'distances': (12.0,),
-                'harmonic': True,
-                'size': 8,
-            }],
-        [(11, 11), None, False,
-            {
-                'shape': (11, 11),
-                'distances': (1/11, 1/11),
-                'harmonic': False,
-                'size': 121,
-            }],
-        [(12, 12), (1.3, 1.3), True,
-            {
-                'shape': (12, 12),
-                'distances': (1.3, 1.3),
-                'harmonic': True,
-                'size': 144,
-            }]
+import nifty5 as ift
 
-    ]
+pmp = pytest.mark.parametrize
+# [shape, distances, harmonic, expected]
+CONSTRUCTOR_CONFIGS = [[(8,), None, False, {
+    'shape': (8,),
+    'distances': (0.125,),
+    'harmonic': False,
+    'size': 8,
+}], [(8,), None, True, {
+    'shape': (8,),
+    'distances': (1.0,),
+    'harmonic': True,
+    'size': 8,
+}], [(8,), (12,), True, {
+    'shape': (8,),
+    'distances': (12.0,),
+    'harmonic': True,
+    'size': 8,
+}], [(11, 11), None, False, {
+    'shape': (11, 11),
+    'distances': (1/11, 1/11),
+    'harmonic': False,
+    'size': 121,
+}], [(12, 12), (1.3, 1.3), True, {
+    'shape': (12, 12),
+    'distances': (1.3, 1.3),
+    'harmonic': True,
+    'size': 144,
+}]]
 
 
 def get_k_length_array_configs():
     # for RGSpace(shape=(4, 4), distances=(0.25,0.25))
     cords_0 = np.ogrid[0:4, 0:4]
-    da_0 = ((cords_0[0] - 4 // 2) * 0.25)**2
+    da_0 = ((cords_0[0] - 4//2)*0.25)**2
     da_0 = np.fft.ifftshift(da_0)
-    temp = ((cords_0[1] - 4 // 2) * 0.25)**2
+    temp = ((cords_0[1] - 4//2)*0.25)**2
     temp = np.fft.ifftshift(temp)
     da_0 = da_0 + temp
     da_0 = np.sqrt(da_0)
     return [
         [(4, 4), (0.25, 0.25), da_0],
-        ]
+    ]
 
 
 def get_dvol_configs():
     np.random.seed(42)
-    return [
-        [(11, 11), None, False, 1],
-        [(11, 11), None, False, 1],
-        [(11, 11), (1.3, 1.3), False, 1],
-        [(11, 11), (1.3, 1.3), False, 1],
-        [(11, 11), None, True, 1],
-        [(11, 11), None, True, 1],
-        [(11, 11), (1.3, 1.3), True, 1],
-        [(11, 11), (1.3, 1.3), True, 1]
-        ]
+    return [[(11, 11), None, False, 1], [(11, 11), None, False, 1],
+            [(11, 11), (1.3, 1.3), False, 1], [(11, 11), (1.3, 1.3), False,
+                                               1], [(11, 11), None, True, 1],
+            [(11, 11), None, True, 1], [(11, 11), (1.3, 1.3), True,
+                                        1], [(11, 11), (1.3, 1.3), True, 1]]
+
+
+@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))
 
 
-class RGSpaceInterfaceTests(unittest.TestCase):
-    @expand([['distances', tuple]])
-    def test_property_ret_type(self, attribute, expected_type):
-        x = ift.RGSpace(1)
-        assert_(isinstance(getattr(x, attribute), expected_type))
+@pmp('shape, distances, harmonic, expected', CONSTRUCTOR_CONFIGS)
+def test_constructor(shape, distances, harmonic, expected):
+    x = ift.RGSpace(shape, distances, harmonic)
+    for key, value in expected.items():
+        assert_equal(getattr(x, key), value)
 
 
-class RGSpaceFunctionalityTests(unittest.TestCase):
-    @expand(CONSTRUCTOR_CONFIGS)
-    def test_constructor(self, shape, distances,
-                         harmonic, expected):
-        x = ift.RGSpace(shape, distances, harmonic)
-        for key, value in expected.items():
-            assert_equal(getattr(x, key), value)
+@pmp('shape, distances, expected', get_k_length_array_configs())
+def test_k_length_array(shape, distances, expected):
+    r = ift.RGSpace(shape=shape, distances=distances, harmonic=True)
+    assert_allclose(r.get_k_length_array().to_global_data(), expected)
 
-    @expand(get_k_length_array_configs())
-    def test_k_length_array(self, shape, distances, expected):
-        r = ift.RGSpace(shape=shape, distances=distances, harmonic=True)
-        assert_allclose(r.get_k_length_array().to_global_data(), expected)
 
-    @expand(get_dvol_configs())
-    def test_dvol(self, shape, distances, harmonic, power):
-        r = ift.RGSpace(shape=shape, distances=distances, harmonic=harmonic)
-        assert_allclose(r.dvol, np.prod(r.distances)**power)
+@pmp('shape, distances, harmonic, power', get_dvol_configs())
+def test_dvol(shape, distances, harmonic, power):
+    r = ift.RGSpace(shape=shape, distances=distances, harmonic=harmonic)
+    assert_allclose(r.dvol, np.prod(r.distances)**power)