diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 5961be29e5643cff0750b4459dd02a3cb048c79b..bb716306045552759da1cd5b9862ab206d58bcc3 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -37,7 +37,7 @@ build_docker_from_cache:
 test_serial:
   stage: test
   script:
-    - pytest-3 -q --cov=nifty6 test
+    - pytest-3 -q --cov=nifty7 test
     - >
       python3 -m coverage report --omit "*plot*" | tee coverage.txt
     - >
diff --git a/ChangeLog b/ChangeLog.md
similarity index 88%
rename from ChangeLog
rename to ChangeLog.md
index 06ae02cfecdfe5f2451619755a7c5e5be79e7412..f95bc720c70aa816b04b2a80f59c1a165adb7f5d 100644
--- a/ChangeLog
+++ b/ChangeLog.md
@@ -1,11 +1,29 @@
-Changes since NIFTy 5:
+Changes since NIFTy 6
+=====================
+
+Naming of operator tests
+------------------------
+
+The implementation tests for nonlinear operators are now available in
+`ift.extra.check_operator()` and for linear operators
+`ift.extra.check_linear_operator()`.
+
+MetricGaussianKL interface
+--------------------------
+
+Users do not instanciate `MetricGaussianKL` by its constructor anymore. Rather
+`MetricGaussianKL.make()` shall be used.
+
+
+Changes since NIFTy 5
+=====================
 
 Minimum Python version increased to 3.6
-=======================================
+---------------------------------------
 
 
 New operators
-=============
+-------------
 
 In addition to the below changes, the following operators were introduced:
 
@@ -20,14 +38,14 @@ In addition to the below changes, the following operators were introduced:
 * IntegrationOperator: Integrates over subspaces of fields
 
 FFT convention adjusted
-=======================
+-----------------------
 
 When going to harmonic space, NIFTy's FFT operator now uses a minus sign in the
 exponent (and, consequently, a plus sign on the adjoint transform). This
 convention is consistent with almost all other numerical FFT libraries.
 
 Interface change in EndomorphicOperator.draw_sample()
-=====================================================
+-----------------------------------------------------
 
 Both complex-valued and real-valued Gaussian probability distributions have
 Hermitian and positive endomorphisms as covariance. Just by looking at an
@@ -43,7 +61,7 @@ In order to dive into those subtleties I suggest running the following code and
 playing around with the dtypes.
 
 ```
-import nifty6 as ift
+import nifty7 as ift
 import numpy as np
 
 dom = ift.UnstructuredDomain(5)
@@ -72,30 +90,30 @@ does not give 1, but sqrt(2) as a result.
 
 
 MPI parallelisation over samples in MetricGaussianKL
-====================================================
+----------------------------------------------------
 
 The classes `MetricGaussianKL` and `MetricGaussianKL_MPI` have been unified
 into one `MetricGaussianKL` class which has MPI support built in.
 
 New approach for random number generation
-=========================================
+-----------------------------------------
 
 The code now uses `numpy`'s new `SeedSequence` and `Generator` classes for the
 production of random numbers (introduced in numpy 1.17. This greatly simplifies
 the generation of reproducible random numbers in the presence of MPI parallelism
 and leads to cleaner code overall. Please see the documentation of
-`nifty6.random` for details.
+`nifty7.random` for details.
 
 
 Interface Change for from_random and OuterProduct
-=================================================
+-------------------------------------------------
 
 The sugar.from_random, Field.from_random, MultiField.from_random now take domain
 as the first argument and default to 'normal' for the second argument.
 Likewise OuterProduct takes domain as the first argument and a field as the second.
 
 Interface Change for non-linear Operators
-=========================================
+-----------------------------------------
 
 The method `Operator.apply()` takes a `Linearization` or a `Field` or a
 `MultiField` as input. This has not changed. However, now each non-linear
@@ -112,7 +130,7 @@ behaviour since both `Operator._check_input()` and
 fulfilled.
 
 Special functions for complete Field reduction operations
-=========================================================
+---------------------------------------------------------
 
 So far, reduction operations called on Fields (like `vdot`, `sum`, `integrate`,
 `mean`, `var`, `std`, `prod` etc.) returned a scalar when the reduction was
@@ -124,7 +142,7 @@ operate over all subdomains and therefore don't take a `spaces` argument; they
 are named `s_vdot`, `s_sum` etc. and always return a scalar.
 
 Updates regarding correlated fields
-===================================
+-----------------------------------
 
 The most commonly used model for homogeneous and isotropic correlated fields in
 nifty5 has been `SLAmplitude` combined with `CorrelatedField`. This model
@@ -141,7 +159,7 @@ via `napprox` breaks the inference scheme with the new model so `napprox` may no
 be used here.
 
 Removal of the standard MPI parallelization scheme:
-===================================================
+---------------------------------------------------
 
 When several MPI tasks are present, NIFTy5 distributes every Field over these
 tasks by splitting it along the first axis. This approach to parallelism is not
diff --git a/README.md b/README.md
index 20ac02efa9b5abd45fbe83bc75b8f6ae941b2454..a2e88606c3a88bdcc4a8e880c168a08186483cc0 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
 NIFTy - Numerical Information Field Theory
 ==========================================
-[![pipeline status](https://gitlab.mpcdf.mpg.de/ift/nifty/badges/NIFTy_6/pipeline.svg)](https://gitlab.mpcdf.mpg.de/ift/nifty/-/commits/NIFTy_6)
-[![coverage report](https://gitlab.mpcdf.mpg.de/ift/nifty/badges/NIFTy_6/coverage.svg)](https://gitlab.mpcdf.mpg.de/ift/nifty/-/commits/NIFTy_6)
+[![pipeline status](https://gitlab.mpcdf.mpg.de/ift/nifty/badges/NIFTy_7/pipeline.svg)](https://gitlab.mpcdf.mpg.de/ift/nifty/-/commits/NIFTy_7)
+[![coverage report](https://gitlab.mpcdf.mpg.de/ift/nifty/badges/NIFTy_7/coverage.svg)](https://gitlab.mpcdf.mpg.de/ift/nifty/-/commits/NIFTy_7)
 
 **NIFTy** project homepage:
 [http://ift.pages.mpcdf.de/nifty](http://ift.pages.mpcdf.de/nifty)
@@ -59,8 +59,8 @@ Optional dependencies:
 
 ### Sources
 
-The current version of NIFTy6 can be obtained by cloning the repository and
-switching to the NIFTy_6 branch:
+The current version of NIFTy7 can be obtained by cloning the repository and
+switching to the NIFTy_7 branch:
 
     git clone https://gitlab.mpcdf.mpg.de/ift/nifty.git
 
@@ -69,10 +69,10 @@ switching to the NIFTy_6 branch:
 In the following, we assume a Debian-based distribution. For other
 distributions, the "apt" lines will need slight changes.
 
-NIFTy6 and its mandatory dependencies can be installed via:
+NIFTy7 and its mandatory dependencies can be installed via:
 
     sudo apt-get install git python3 python3-pip python3-dev
-    pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_6
+    pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_7
 
 Plotting support is added via:
 
@@ -108,7 +108,7 @@ To run the tests, additional packages are required:
 Afterwards the tests (including a coverage report) can be run using the
 following command in the repository root:
 
-    pytest-3 --cov=nifty6 test
+    pytest-3 --cov=nifty7 test
 
 ### First Steps
 
diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py
index b36c0d500cc68638e29a71a1e5650f8a77f7d72d..be6a57332ae77780b05bcd73c9c0be1e9340d32c 100644
--- a/demos/bench_gridder.py
+++ b/demos/bench_gridder.py
@@ -11,7 +11,7 @@
 # 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) 2019 Max-Planck-Society
+# Copyright(C) 2019-2020 Max-Planck-Society
 # Author: Martin Reinecke
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
@@ -21,70 +21,75 @@ from time import time
 import matplotlib.pyplot as plt
 import numpy as np
 
-import nifty6 as ift
-
-N0s, a0s, b0s, c0s = [], [], [], []
-
-for ii in range(10, 26):
-    nu = 1024
-    nv = 1024
-    N = int(2**ii)
-    print('N = {}'.format(N))
-
-    rng = ift.random.current_rng()
-    uv = rng.uniform(-.5, .5, (N,2))
-    vis = rng.normal(0., 1., N) + 1j*rng.normal(0., 1., N)
-
-    uvspace = ift.RGSpace((nu, nv))
-
-    visspace = ift.UnstructuredDomain(N)
-
-    img = rng.standard_normal((nu, nv))
-    img = ift.makeField(uvspace, img)
-
-    t0 = time()
-    GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
-    vis = ift.makeField(visspace, vis)
-    op = GM.getFull().adjoint
-    t1 = time()
-    op(img).val
-    t2 = time()
-    op.adjoint(vis).val
-    t3 = time()
-    print(t2-t1, t3-t2)
-    N0s.append(N)
-    a0s.append(t1 - t0)
-    b0s.append(t2 - t1)
-    c0s.append(t3 - t2)
-
-print('Measure rest operator')
-sc = ift.StatCalculator()
-op = GM.getRest().adjoint
-for _ in range(10):
-    t0 = time()
-    res = op(img)
-    sc.add(time() - t0)
-t_fft = sc.mean
-print('FFT shape', res.shape)
-
-plt.scatter(N0s, a0s, label='Gridder mr')
-plt.legend()
-# no idea why this is necessary, but if it is omitted, the range is wrong
-plt.ylim(min(a0s), max(a0s))
-plt.ylabel('time [s]')
-plt.title('Initialization')
-plt.loglog()
-plt.savefig('bench0.png')
-plt.close()
-
-plt.scatter(N0s, b0s, color='k', marker='^', label='Gridder mr times')
-plt.scatter(N0s, c0s, color='k', label='Gridder mr adjoint times')
-plt.axhline(sc.mean, label='FFT')
-plt.axhline(sc.mean + np.sqrt(sc.var))
-plt.axhline(sc.mean - np.sqrt(sc.var))
-plt.legend()
-plt.ylabel('time [s]')
-plt.title('Apply')
-plt.loglog()
-plt.savefig('bench1.png')
-plt.close()
+import nifty7 as ift
+
+
+def main():
+    N0s, a0s, b0s, c0s = [], [], [], []
+
+    for ii in range(10, 26):
+        nu = 1024
+        nv = 1024
+        N = int(2**ii)
+        print('N = {}'.format(N))
+
+        rng = ift.random.current_rng()
+        uv = rng.uniform(-.5, .5, (N, 2))
+        vis = rng.normal(0., 1., N) + 1j*rng.normal(0., 1., N)
+
+        uvspace = ift.RGSpace((nu, nv))
+
+        visspace = ift.UnstructuredDomain(N)
+
+        img = rng.standard_normal((nu, nv))
+        img = ift.makeField(uvspace, img)
+
+        t0 = time()
+        GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
+        vis = ift.makeField(visspace, vis)
+        op = GM.getFull().adjoint
+        t1 = time()
+        op(img).val
+        t2 = time()
+        op.adjoint(vis).val
+        t3 = time()
+        print(t2-t1, t3-t2)
+        N0s.append(N)
+        a0s.append(t1 - t0)
+        b0s.append(t2 - t1)
+        c0s.append(t3 - t2)
+
+    print('Measure rest operator')
+    sc = ift.StatCalculator()
+    op = GM.getRest().adjoint
+    for _ in range(10):
+        t0 = time()
+        res = op(img)
+        sc.add(time() - t0)
+    print('FFT shape', res.shape)
+
+    plt.scatter(N0s, a0s, label='Gridder mr')
+    plt.legend()
+    # no idea why this is necessary, but if it is omitted, the range is wrong
+    plt.ylim(min(a0s), max(a0s))
+    plt.ylabel('time [s]')
+    plt.title('Initialization')
+    plt.loglog()
+    plt.savefig('bench0.png')
+    plt.close()
+
+    plt.scatter(N0s, b0s, color='k', marker='^', label='Gridder mr times')
+    plt.scatter(N0s, c0s, color='k', label='Gridder mr adjoint times')
+    plt.axhline(sc.mean, label='FFT')
+    plt.axhline(sc.mean + np.sqrt(sc.var))
+    plt.axhline(sc.mean - np.sqrt(sc.var))
+    plt.legend()
+    plt.ylabel('time [s]')
+    plt.title('Apply')
+    plt.loglog()
+    plt.savefig('bench1.png')
+    plt.close()
+
+
+if __name__ == '__main__':
+    main()
diff --git a/demos/bernoulli_demo.py b/demos/bernoulli_demo.py
index bd727db9eca92b2589a7c207e3b5783a6e69a84a..9aaca61bf6b3eee62af7c0c28465c747d03471a1 100644
--- a/demos/bernoulli_demo.py
+++ b/demos/bernoulli_demo.py
@@ -11,7 +11,7 @@
 # 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
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -24,9 +24,10 @@
 
 import numpy as np
 
-import nifty6 as ift
+import nifty7 as ift
 
-if __name__ == '__main__':
+
+def main():
     # Set up the position space of the signal
     mode = 2
     if mode == 0:
@@ -84,3 +85,7 @@ if __name__ == '__main__':
     plot.add(GR.adjoint_times(data), title='data')
     plot.add(sky(mock_position), title='truth')
     plot.output(nx=3, xsize=16, ysize=9, title="results", name="bernoulli.png")
+
+
+if __name__ == '__main__':
+    main()
diff --git a/demos/getting_started_0.ipynb b/demos/getting_started_0.ipynb
index 63ce572518e499b2ebf1a85477f00ec1d292fd5b..85678373c3ca1dc3d55f1a2bdb08a19aef862b8b 100644
--- a/demos/getting_started_0.ipynb
+++ b/demos/getting_started_0.ipynb
@@ -140,7 +140,7 @@
    "outputs": [],
    "source": [
     "import numpy as np\n",
-    "import nifty6 as ift\n",
+    "import nifty7 as ift\n",
     "import matplotlib.pyplot as plt\n",
     "%matplotlib inline"
    ]
diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py
index 74b2a1b7c3308ae28c3af92300ebb31bf1a1d952..f3d44b26a602d96c58ef7ec8701556d24b6054c8 100644
--- a/demos/getting_started_1.py
+++ b/demos/getting_started_1.py
@@ -11,7 +11,7 @@
 # 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
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -25,7 +25,7 @@ import sys
 
 import numpy as np
 
-import nifty6 as ift
+import nifty7 as ift
 
 
 def make_checkerboard_mask(position_space):
@@ -38,14 +38,14 @@ def make_checkerboard_mask(position_space):
     return mask
 
 
-def make_random_mask():
+def make_random_mask(domain):
     # Random mask for spherical mode
-    mask = ift.from_random(position_space, 'pm1')
+    mask = ift.from_random(domain, 'pm1')
     mask = (mask + 1)/2
     return mask.val
 
 
-if __name__ == '__main__':
+def main():
     # Choose space on which the signal field is defined
     if len(sys.argv) == 2:
         mode = int(sys.argv[1])
@@ -63,7 +63,7 @@ if __name__ == '__main__':
     else:
         # Sphere with half of its pixels randomly masked
         position_space = ift.HPSpace(128)
-        mask = make_random_mask()
+        mask = make_random_mask(position_space)
 
     # Specify harmonic space corresponding to signal
     harmonic_space = position_space.get_default_codomain()
@@ -148,3 +148,7 @@ if __name__ == '__main__':
                  title='Residuals')
         plot.output(nx=2, ny=2, xsize=10, ysize=10, name=filename)
     print("Saved results as '{}'.".format(filename))
+
+
+if __name__ == '__main__':
+    main()
diff --git a/demos/getting_started_2.py b/demos/getting_started_2.py
index 8b7c0cb6c3f781e54728e433df768ade93c42fcd..c792a14803665b5ab81df8a3f40f6d84ef9f9d6b 100644
--- a/demos/getting_started_2.py
+++ b/demos/getting_started_2.py
@@ -11,7 +11,7 @@
 # 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
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -25,25 +25,23 @@ import sys
 
 import numpy as np
 
-import nifty6 as ift
+import nifty7 as ift
 
 
-def exposure_2d():
+def exposure_2d(domain):
     # Structured exposure for 2D mode
-    x_shape, y_shape = position_space.shape
-
-    exposure = np.ones(position_space.shape)
+    x_shape, y_shape = domain.shape
+    exposure = np.ones(domain.shape)
     exposure[x_shape//3:x_shape//2, :] *= 2.
     exposure[x_shape*4//5:x_shape, :] *= .1
     exposure[x_shape//2:x_shape*3//2, :] *= 3.
     exposure[:, x_shape//3:x_shape//2] *= 2.
     exposure[:, x_shape*4//5:x_shape] *= .1
     exposure[:, x_shape//2:x_shape*3//2] *= 3.
-
-    return ift.Field.from_raw(position_space, exposure)
+    return ift.Field.from_raw(domain, exposure)
 
 
-if __name__ == '__main__':
+def main():
     # Choose space on which the signal field is defined
     if len(sys.argv) == 2:
         mode = int(sys.argv[1])
@@ -57,7 +55,7 @@ if __name__ == '__main__':
     elif mode == 1:
         # Two-dimensional regular grid with inhomogeneous exposure
         position_space = ift.RGSpace([512, 512])
-        exposure = exposure_2d()
+        exposure = exposure_2d(position_space)
     else:
         # Sphere with uniform exposure of 100
         position_space = ift.HPSpace(128)
@@ -118,3 +116,7 @@ if __name__ == '__main__':
     plot.add(reconst - signal, title='Residuals')
     plot.output(xsize=12, ysize=10, name=filename)
     print("Saved results as '{}'.".format(filename))
+
+
+if __name__ == '__main__':
+    main()
diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py
index 74d96e516c66d03f81f11509ac37aa98e5edfdc6..21f90e563ab9bcbdd8852496843f69a4bd8a8df9 100644
--- a/demos/getting_started_3.py
+++ b/demos/getting_started_3.py
@@ -11,7 +11,7 @@
 # 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
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -29,7 +29,7 @@ import sys
 
 import numpy as np
 
-import nifty6 as ift
+import nifty7 as ift
 
 
 def random_los(n_los):
@@ -44,7 +44,7 @@ def radial_los(n_los):
     return starts, ends
 
 
-if __name__ == '__main__':
+def main():
     # Choose between random line-of-sight response (mode=0) and radial lines
     # of sight (mode=1)
     if len(sys.argv) == 2:
@@ -60,10 +60,10 @@ if __name__ == '__main__':
     #  see 'getting_started_4_CorrelatedFields.ipynb'.
 
     cfmaker = ift.CorrelatedFieldMaker.make(
-        offset_mean =      0.0,  # 0.
-        offset_std_mean = 1e-3,  # 1e-3
-        offset_std_std =  1e-6,  # 1e-6
-        prefix = '')
+        offset_mean=      0.0,
+        offset_std_mean= 1e-3,
+        offset_std_std=  1e-6,
+        prefix='')
 
     fluctuations_dict = {
         # Amplitude of field fluctuations
@@ -72,7 +72,7 @@ if __name__ == '__main__':
 
         # Exponent of power law power spectrum component
         'loglogavgslope_mean': -2.0,  # -3.0
-        'loglogavgslope_stddev': 0.5,  #  0.5
+        'loglogavgslope_stddev': 0.5,  # 0.5
 
         # Amplitude of integrated Wiener process power spectrum component
         'flexibility_mean':   2.5,  # 1.0
@@ -131,7 +131,7 @@ if __name__ == '__main__':
     # Draw new samples to approximate the KL five times
     for i in range(5):
         # Draw new samples and minimize KL
-        KL = ift.MetricGaussianKL(mean, H, N_samples)
+        KL = ift.MetricGaussianKL.make(mean, H, N_samples)
         KL, convergence = minimizer(KL)
         mean = KL.position
 
@@ -144,7 +144,7 @@ if __name__ == '__main__':
                     name=filename.format("loop_{:02d}".format(i)))
 
     # Draw posterior samples
-    KL = ift.MetricGaussianKL(mean, H, N_samples)
+    KL = ift.MetricGaussianKL.make(mean, H, N_samples)
     sc = ift.StatCalculator()
     for sample in KL.samples:
         sc.add(signal(sample + KL.position))
@@ -163,3 +163,7 @@ if __name__ == '__main__':
         linewidth=[1.]*len(powers) + [3., 3.])
     plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res)
     print("Saved results as '{}'.".format(filename_res))
+
+
+if __name__ == '__main__':
+    main()
diff --git a/demos/getting_started_4_CorrelatedFields.ipynb b/demos/getting_started_4_CorrelatedFields.ipynb
index f8476df63503c4e6f0f4ad8121a454f2d407b94f..225ad195cb073c0feb68c67c81a4bf2df84355ce 100644
--- a/demos/getting_started_4_CorrelatedFields.ipynb
+++ b/demos/getting_started_4_CorrelatedFields.ipynb
@@ -29,7 +29,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import nifty6 as ift\n",
+    "import nifty7 as ift\n",
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "ift.random.push_sseq_from_seed(43)\n",
@@ -46,7 +46,7 @@
    "source": [
     "# Plotting routine\n",
     "def plot(fields, spectra, title=None):\n",
-    "    # Plotting preparation is normally handled by nifty6.Plot\n",
+    "    # Plotting preparation is normally handled by nifty7.Plot\n",
     "    # It is done manually here to be able to tweak details\n",
     "    # Fields are assumed to have identical domains\n",
     "    fig = plt.figure(tight_layout=True, figsize=(12, 3.5))\n",
diff --git a/demos/getting_started_5_mf.py b/demos/getting_started_5_mf.py
index 5a00b38cfd52c1c396e93eb9a02e3a5124999971..84d9b77cd4e6b6347c00ca67a0ea5e1d22c688ee 100644
--- a/demos/getting_started_5_mf.py
+++ b/demos/getting_started_5_mf.py
@@ -11,7 +11,7 @@
 # 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
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -29,7 +29,7 @@ import sys
 
 import numpy as np
 
-import nifty6 as ift
+import nifty7 as ift
 
 
 class SingleDomain(ift.LinearOperator):
@@ -55,7 +55,7 @@ def radial_los(n_los):
     return starts, ends
 
 
-if __name__ == '__main__':
+def main():
     # Choose between random line-of-sight response (mode=0) and radial lines
     # of sight (mode=1)
     if len(sys.argv) == 2:
@@ -83,7 +83,7 @@ if __name__ == '__main__':
     A2 = cfmaker.normalized_amplitudes[1]
     DC = SingleDomain(correlated_field.target, position_space)
 
-    ## Apply a nonlinearity
+    # Apply a nonlinearity
     signal = DC @ ift.sigmoid(correlated_field)
 
     # Build the line-of-sight response and define signal response
@@ -118,7 +118,7 @@ if __name__ == '__main__':
     ic_newton.enable_logging()
     minimizer = ift.NewtonCG(ic_newton, enable_logging=True)
 
-    ## number of samples used to estimate the KL
+    # number of samples used to estimate the KL
     N_samples = 20
 
     # Set up likelihood and information Hamiltonian
@@ -131,7 +131,7 @@ if __name__ == '__main__':
 
     for i in range(10):
         # Draw new samples and minimize KL
-        KL = ift.MetricGaussianKL(mean, H, N_samples)
+        KL = ift.MetricGaussianKL.make(mean, H, N_samples)
         KL, convergence = minimizer(KL)
         mean = KL.position
 
@@ -157,8 +157,7 @@ if __name__ == '__main__':
                     name=filename.format("loop_{:02d}".format(i)))
 
     # Done, draw posterior samples
-    Nsamples = 20
-    KL = ift.MetricGaussianKL(mean, H, N_samples)
+    KL = ift.MetricGaussianKL.make(mean, H, N_samples)
     sc = ift.StatCalculator()
     scA1 = ift.StatCalculator()
     scA2 = ift.StatCalculator()
@@ -189,3 +188,7 @@ if __name__ == '__main__':
              linewidth=[1.]*len(powers2) + [3., 3.])
     plot.output(ny=2, nx=2, xsize=15, ysize=15, name=filename_res)
     print("Saved results as '{}'.".format(filename_res))
+
+
+if __name__ == '__main__':
+    main()
diff --git a/demos/mgvi_visualized.py b/demos/mgvi_visualized.py
index bd67bc48536367e870e4dfbbadda9c01f12f80bb..69e0373ede3aed5895ee51a87132e9fb9536dca0 100644
--- a/demos/mgvi_visualized.py
+++ b/demos/mgvi_visualized.py
@@ -32,9 +32,10 @@ import numpy as np
 import pylab as plt
 from matplotlib.colors import LogNorm
 
-import nifty6 as ift
+import nifty7 as ift
 
-if __name__ == '__main__':
+
+def main():
     dom = ift.UnstructuredDomain(1)
     scale = 10
 
@@ -90,7 +91,7 @@ if __name__ == '__main__':
     plt.figure(figsize=[12, 8])
     for ii in range(15):
         if ii % 3 == 0:
-            mgkl = ift.MetricGaussianKL(pos, ham, 40)
+            mgkl = ift.MetricGaussianKL.make(pos, ham, 40)
 
         plt.cla()
         plt.imshow(z.T, origin='lower', norm=LogNorm(), vmin=1e-3,
@@ -119,3 +120,7 @@ if __name__ == '__main__':
     ift.logger.info('Finished')
     # Uncomment the following line in order to leave the plots open
     # plt.show()
+
+
+if __name__ == '__main__':
+    main()
diff --git a/demos/misc/convolution.py b/demos/misc/convolution.py
index 268435e1fc9a544649916f7748dc843b3a4fcb10..7a83d264d9e51093179f0f8f9b1f9184e2fde955 100644
--- a/demos/misc/convolution.py
+++ b/demos/misc/convolution.py
@@ -11,12 +11,12 @@
 # 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) 2019 Max-Planck-Society
+# Copyright(C) 2019-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import numpy as np
-import nifty6 as ift
+import nifty7 as ift
 
 
 def convtest(test_signal, delta, func):
@@ -37,55 +37,52 @@ def convtest(test_signal, delta, func):
 
     # Plot results
     plot = ift.Plot()
-    plot.add(signal, title='Signal')
+    plot.add(test_signal, title='Signal')
     plot.add(conv_signal, title='Signal Convolved')
     plot.add(cac_signal, title='Signal, Conv, Adj-Conv')
     plot.add(conv_delta, title='Kernel')
     plot.output()
 
 
-# Healpix test
-nside = 64
-npix = 12 * nside * nside
+def main():
+    # Healpix test
+    nside = 64
+    npix = 12 * nside * nside
 
-domain = ift.HPSpace(nside)
+    domain = ift.HPSpace(nside)
 
-# Define test signal (some point sources)
-signal_vals = np.zeros(npix, dtype=np.float64)
-for i in range(0, npix, npix//12 + 27):
-    signal_vals[i] = 500.
+    # Define test signal (some point sources)
+    signal_vals = np.zeros(npix, dtype=np.float64)
+    for i in range(0, npix, npix//12 + 27):
+        signal_vals[i] = 500.
 
-signal = ift.makeField(domain, signal_vals)
+    signal = ift.makeField(domain, signal_vals)
 
-delta_vals = np.zeros(npix, dtype=np.float64)
-delta_vals[0] = 1.0
-delta = ift.makeField(domain, delta_vals)
+    delta_vals = np.zeros(npix, dtype=np.float64)
+    delta_vals[0] = 1.0
+    delta = ift.makeField(domain, delta_vals)
 
+    # Define kernel function
+    def func(theta):
+        ct = np.cos(theta)
+        return 1. * np.logical_and(ct > 0.7, ct <= 0.8)
 
-# Define kernel function
-def func(theta):
-    ct = np.cos(theta)
-    return 1. * np.logical_and(ct > 0.7, ct <= 0.8)
+    convtest(signal, delta, func)
 
+    domain = ift.RGSpace((100, 100))
+    # Define test signal (some point sources)
+    signal_vals = np.zeros(domain.shape, dtype=np.float64)
+    signal_vals[35, 70] = 5000.
+    signal_vals[45, 8] = 5000.
+    signal = ift.makeField(domain, signal_vals)
 
-convtest(signal, delta, func)
+    # Define delta signal, generate kernel image
+    delta_vals = np.zeros(domain.shape, dtype=np.float64)
+    delta_vals[0, 0] = 1.0
+    delta = ift.makeField(domain, delta_vals)
 
-domain = ift.RGSpace((100, 100))
-# Define test signal (some point sources)
-signal_vals = np.zeros(domain.shape, dtype=np.float64)
-signal_vals[35, 70] = 5000.
-signal_vals[45, 8] = 5000.
-signal = ift.makeField(domain, signal_vals)
+    convtest(signal, delta, lambda d: 1. * np.logical_and(d > 0.1, d <= 0.2))
 
-# Define delta signal, generate kernel image
-delta_vals = np.zeros(domain.shape, dtype=np.float64)
-delta_vals[0, 0] = 1.0
-delta = ift.makeField(domain, delta_vals)
 
-
-# Define kernel function
-def func(dist):
-    return 1. * np.logical_and(dist > 0.1, dist <= 0.2)
-
-
-convtest(signal, delta, func)
+if __name__ == '__main__':
+    main()
diff --git a/demos/polynomial_fit.py b/demos/polynomial_fit.py
index 2db06700ca9f6067726ed4e8383a95e04c579acf..8b1e860744bf285e95156064cc206fa767d0f930 100644
--- a/demos/polynomial_fit.py
+++ b/demos/polynomial_fit.py
@@ -11,14 +11,15 @@
 # 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
+# Copyright(C) 2013-2020 Max-Planck-Society
+# Author: Philipp Arras
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import matplotlib.pyplot as plt
 import numpy as np
 
-import nifty6 as ift
+import nifty7 as ift
 
 
 def polynomial(coefficients, sampling_points):
@@ -56,7 +57,7 @@ class PolynomialResponse(ift.LinearOperator):
 
     def __init__(self, domain, sampling_points):
         if not (isinstance(domain, ift.UnstructuredDomain)
-                and isinstance(x, np.ndarray)):
+                and isinstance(sampling_points, np.ndarray)):
             raise TypeError
         self._domain = ift.DomainTuple.make(domain)
         tgt = ift.UnstructuredDomain(sampling_points.shape)
@@ -80,7 +81,7 @@ class PolynomialResponse(ift.LinearOperator):
         return ift.makeField(self._tgt(mode), out)
 
 
-if __name__ == '__main__':
+def main():
     # Generate some mock data
     N_params = 10
     N_samples = 100
@@ -96,7 +97,7 @@ if __name__ == '__main__':
     p_space = ift.UnstructuredDomain(N_params)
     params = ift.full(p_space, 0.)
     R = PolynomialResponse(p_space, x)
-    ift.extra.consistency_check(R)
+    ift.extra.check_linear_operator(R)
 
     d_space = R.target
     d = ift.makeField(d_space, y)
@@ -140,3 +141,7 @@ if __name__ == '__main__':
     sigma = np.sqrt(sc.var.val)
     for ii in range(len(mean)):
         print('Coefficient x**{}: {:.2E} +/- {:.2E}'.format(ii, mean[ii], sigma[ii]))
+
+
+if __name__ == '__main__':
+    main()
diff --git a/docs/generate.sh b/docs/generate.sh
index 75e9f57855eb437a51af27d6b2e196256c963559..d9f07447976f6542bbc46267721f8ab7bd65e63c 100755
--- a/docs/generate.sh
+++ b/docs/generate.sh
@@ -1,4 +1,4 @@
 rm -rf docs/build docs/source/mod
-EXCLUDE="nifty6/logger.py nifty6/git_version.py"
-sphinx-apidoc -e -o docs/source/mod nifty6 ${EXCLUDE}
+EXCLUDE="nifty7/logger.py nifty7/git_version.py"
+sphinx-apidoc -e -o docs/source/mod nifty7 ${EXCLUDE}
 sphinx-build -b html docs/source/ docs/build/
diff --git a/docs/source/citations.rst b/docs/source/citations.rst
index 802329a504ee7181086622e3f145719b893e8797..b67d33cfba75f68ae3e042b2632abeaa3a36c65c 100644
--- a/docs/source/citations.rst
+++ b/docs/source/citations.rst
@@ -13,7 +13,7 @@ NIFTy-related publications
       author = {{Martin Reinecke, Theo Steininger, Marco Selig}},
       title = {NIFTy -- Numerical Information Field TheorY},
       url = {https://gitlab.mpcdf.mpg.de/ift/NIFTy},
-      version = {nifty6},
+      version = {nifty7},
       date = {2018-04-05},
       }
 
diff --git a/docs/source/code.rst b/docs/source/code.rst
index ea08950906a500bfee710b0c01ba33a6548d7989..554c7beec5319178f710d84eb6c618da3cd6b8fd 100644
--- a/docs/source/code.rst
+++ b/docs/source/code.rst
@@ -36,9 +36,9 @@ Domains
 Abstract base class
 -------------------
 
-.. currentmodule:: nifty6.domains.domain
+.. currentmodule:: nifty7.domains.domain
 
-One of the fundamental building blocks of the NIFTy6 framework is the *domain*.
+One of the fundamental building blocks of the NIFTy7 framework is the *domain*.
 Its required capabilities are expressed by the abstract :py:class:`Domain` class.
 A domain must be able to answer the following queries:
 
@@ -52,7 +52,7 @@ A domain must be able to answer the following queries:
 Unstructured domains
 --------------------
 
-.. currentmodule:: nifty6.domains.unstructured_domain
+.. currentmodule:: nifty7.domains.unstructured_domain
 
 Domains can be either *structured* (i.e. there is geometrical information
 associated with them, like position in space and volume factors),
@@ -65,7 +65,7 @@ Unstructured domains can be described by instances of NIFTy's
 Structured domains
 ------------------
 
-.. currentmodule:: nifty6.domains.structured_domain
+.. currentmodule:: nifty7.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.
@@ -86,7 +86,7 @@ The additional methods are specified in the abstract class
 
 NIFTy comes with several concrete subclasses of :class:`StructuredDomain`:
 
-.. currentmodule:: nifty6.domains
+.. currentmodule:: nifty7.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.
@@ -123,7 +123,7 @@ Some examples are:
   (on a harmonic domain) and a few inferred model parameters (e.g. on an
   unstructured grid).
 
-.. currentmodule:: nifty6
+.. currentmodule:: nifty7
 
 Consequently, NIFTy defines a class called :class:`~domain_tuple.DomainTuple`
 holding a sequence of :class:`~domains.domain.Domain` objects. The full domain is
@@ -228,7 +228,7 @@ Advanced operators
 NIFTy provides a library of commonly employed operators which can be used for
 specific inference problems. Currently these are:
 
-.. currentmodule:: nifty6.library
+.. currentmodule:: nifty7.library
 
 - :class:`~smooth_linear_amplitude.SLAmplitude`, which returns a smooth power spectrum.
 - :class:`~inverse_gamma_operator.InverseGammaOperator`, which models point sources which are
@@ -240,9 +240,9 @@ specific inference problems. Currently these are:
 Linear Operators
 ================
 
-.. currentmodule:: nifty6.operators
+.. currentmodule:: nifty7.operators
 
-A linear operator (represented by NIFTy6's abstract
+A linear operator (represented by NIFTy7's abstract
 :class:`~linear_operator.LinearOperator` class) is derived from
 :class:`~operator.Operator` and can be interpreted as an (implicitly defined)
 matrix. Since its operation is linear, it can provide some additional
@@ -282,7 +282,7 @@ This functionality is provided by NIFTy's
 :class:`~inversion_enabler.InversionEnabler` class, which is itself a linear
 operator.
 
-.. currentmodule:: nifty6.operators.operator
+.. currentmodule:: nifty7.operators.operator
 
 Direct multiplication and adjoint inverse multiplication transform a field
 defined on the operator's :attr:`~Operator.domain` to one defined on the
@@ -290,7 +290,7 @@ operator's :attr:`~Operator.target`, whereas adjoint multiplication and inverse
 multiplication transform from :attr:`~Operator.target` to
 :attr:`~Operator.domain`.
 
-.. currentmodule:: nifty6.operators
+.. currentmodule:: nifty7.operators
 
 Operators with identical domain and target can be derived from
 :class:`~endomorphic_operator.EndomorphicOperator`. Typical examples for this
@@ -299,7 +299,7 @@ multiplies its input by a scalar value, and
 :class:`~diagonal_operator.DiagonalOperator`, which multiplies every value of
 its input field with potentially different values.
 
-.. currentmodule:: nifty6
+.. currentmodule:: nifty7
 
 Further operator classes provided by NIFTy are
 
@@ -325,7 +325,7 @@ and ``f1`` and ``f2`` are of type :class:`~field.Field`, writing::
     X = A(B.inverse(A.adjoint)) + C
     f2 = X(f1)
 
-.. currentmodule:: nifty6.operators.linear_operator
+.. currentmodule:: nifty7.operators.linear_operator
 
 will perform the operation suggested intuitively by the notation, checking
 domain compatibility while building the composed operator.
@@ -349,12 +349,12 @@ Minimization
 Most problems in IFT are solved by (possibly nested) minimizations of
 high-dimensional functions, which are often nonlinear.
 
-.. currentmodule:: nifty6.minimization
+.. currentmodule:: nifty7.minimization
 
 Energy functionals
 ------------------
 
-In NIFTy6 such functions are represented by objects of type
+In NIFTy7 such functions are represented by objects of type
 :class:`~energy.Energy`. These hold the prescription how to calculate the
 function's :attr:`~energy.Energy.value`, :attr:`~energy.Energy.gradient` and
 (optionally) :attr:`~energy.Energy.metric` at any given
@@ -362,11 +362,11 @@ function's :attr:`~energy.Energy.value`, :attr:`~energy.Energy.gradient` and
 floating-point scalars, gradients have the form of fields defined on the energy's
 position domain, and metrics are represented by linear operator objects.
 
-.. currentmodule:: nifty6
+.. currentmodule:: nifty7
 
 Energies are classes that typically have to be provided by the user when
 tackling new IFT problems. An example of concrete energy classes delivered with
-NIFTy6 is :class:`~minimization.quadratic_energy.QuadraticEnergy` (with
+NIFTy7 is :class:`~minimization.quadratic_energy.QuadraticEnergy` (with
 position-independent metric, mainly used with conjugate gradient minimization).
 
 For MGVI, NIFTy provides the :class:`~energy.Energy` subclass
@@ -394,14 +394,14 @@ functional for MGVI and minimize it.
 Iteration control
 -----------------
 
-.. currentmodule:: nifty6.minimization.iteration_controllers
+.. currentmodule:: nifty7.minimization.iteration_controllers
 
 Iterative minimization of an energy requires some means of
 checking the quality of the current solution estimate and stopping once
 it is sufficiently accurate. In case of numerical problems, the iteration needs
 to be terminated as well, returning a suitable error description.
 
-In NIFTy6, this functionality is encapsulated in the abstract
+In NIFTy7, this functionality is encapsulated in the abstract
 :class:`IterationController` class, which is provided with the initial energy
 object before starting the minimization, and is updated with the improved
 energy after every iteration. Based on this information, it can either continue
@@ -418,7 +418,7 @@ in NIFTy can be found below, but users have complete freedom to implement custom
 Minimization algorithms
 -----------------------
 
-.. currentmodule:: nifty6.minimization
+.. currentmodule:: nifty7.minimization
 
 All minimization algorithms in NIFTy inherit from the abstract
 :class:`~minimizer.Minimizer` class, which presents a minimalistic interface
@@ -459,7 +459,7 @@ available under the names :class:`~scipy_minimizer.ScipyCG` and
 Application to operator inversion
 ---------------------------------
 
-.. currentmodule:: nifty6
+.. currentmodule:: nifty7
 
 The machinery presented here cannot only be used for minimizing functionals
 derived from IFT, but also for the numerical inversion of linear operators, if
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 1a0af1502216dac01ade1468a3c984b8a0b66763..d9cd9beb803d2dc3195b9a13173f6a84a17558e1 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -1,4 +1,4 @@
-import nifty6
+import nifty7
 
 extensions = [
     'sphinx.ext.napoleon',  # Support for NumPy and Google style docstrings
@@ -15,11 +15,11 @@ napoleon_use_admonition_for_examples = True
 napoleon_use_admonition_for_references = True
 napoleon_include_special_with_doc = True
 
-project = u'NIFTy6'
+project = u'NIFTy7'
 copyright = u'2013-2019, Max-Planck-Society'
 author = u'Martin Reinecke'
 
-release = nifty6.version.__version__
+release = nifty7.version.__version__
 version = release[:-2]
 
 language = None
@@ -30,5 +30,5 @@ html_theme = "sphinx_rtd_theme"
 html_logo = 'nifty_logo_black.png'
 
 exclude_patterns = [
-    'mod/modules.rst', 'mod/nifty6.git_version.rst', 'mod/nifty6.logger.rst'
+    'mod/modules.rst', 'mod/nifty7.git_version.rst', 'mod/nifty7.logger.rst'
 ]
diff --git a/docs/source/ift.rst b/docs/source/ift.rst
index d3dbd6649247a26ff4982e3f0861a130f389cf57..7d151e53abcebb508323c01755372096113312fa 100644
--- a/docs/source/ift.rst
+++ b/docs/source/ift.rst
@@ -174,7 +174,7 @@ It only requires minimizing the information Hamiltonian, e.g. by a gradient desc
 
     \frac{\partial \mathcal{H}(d,\xi)}{\partial \xi} = 0.
 
-NIFTy6 automatically calculates the necessary gradient from a generative model of the signal and the data and uses this to minimize the Hamiltonian.
+NIFTy7 automatically calculates the necessary gradient from a generative model of the signal and the data and uses this to minimize the Hamiltonian.
 
 However, MAP often provides unsatisfactory results in cases of deep hirachical Bayesian networks.
 The reason for this is that MAP ignores the volume factors in parameter space, which are not to be neglected in deciding whether a solution is reasonable or not.
@@ -224,7 +224,7 @@ Thus, only the gradient of the KL is needed with respect to this, which can be e
 
 We stochastically estimate the KL-divergence and gradients with a set of samples drawn from the approximate posterior distribution.
 The particular structure of the covariance allows us to draw independent samples solving a certain system of equations.
-This KL-divergence for MGVI is implemented in the class :class:`~minimization.metric_gaussian_kl.MetricGaussianKL` within NIFTy6.
+This KL-divergence for MGVI is implemented in the class :class:`~minimization.metric_gaussian_kl.MetricGaussianKL` within NIFTy7.
 
 
 The demo `getting_started_3.py` for example not only infers a field this way, but also the power spectrum of the process that has generated the field.
diff --git a/docs/source/index.rst b/docs/source/index.rst
index bf6ae63256c1b1191fd4db8da4e4770365c0a774..6f157b23bbddb176894ea6fbb87b1e3b9e224318 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -31,4 +31,4 @@ Contents
    installation
    code
    citations
-   Package Documentation <mod/nifty6>
+   Package Documentation <mod/nifty7>
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
index 0313bd7a669bbcbc8bdb073593de5c2a59894633..a2c61c921b2529adf76d834c9019deac3f036ec1 100644
--- a/docs/source/installation.rst
+++ b/docs/source/installation.rst
@@ -5,10 +5,10 @@ Installation
 In the following, we assume a Debian-based Linux distribution. For other
 distributions, the "apt" lines will need slight changes.
 
-NIFTy6 and its mandatory dependencies can be installed via::
+NIFTy7 and its mandatory dependencies can be installed via::
 
     sudo apt-get install git python3 python3-pip python3-dev
-    pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_6
+    pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_7
 
 Plotting support is added via::
 
diff --git a/docs/source/volume.rst b/docs/source/volume.rst
index 2b1318f09b9a41e1dcd51648b258ec7c476b758f..6cdd2e8a73e639b9cdef32385a57bb24d967aed1 100644
--- a/docs/source/volume.rst
+++ b/docs/source/volume.rst
@@ -148,7 +148,7 @@ A visualisation of this can be seen in figure 2, which displays the MAP inferenc
 Implementation in NIFTy
 .......................
 
-.. currentmodule:: nifty6
+.. currentmodule:: nifty7
 
 Most codes in NIFTy will contain the description of a measurement process or,
 more generally, a log-likelihood.
diff --git a/nifty7 b/nifty7
new file mode 120000
index 0000000000000000000000000000000000000000..aa8e45f12bce4c9ccf69652a97bd54ffc43c7123
--- /dev/null
+++ b/nifty7
@@ -0,0 +1 @@
+src/
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 0e405659bc875c400116180fb7d36db3f25e531b..84fc04372952b6cb7dca1d64eed0281f10e78b43 100644
--- a/setup.py
+++ b/setup.py
@@ -23,20 +23,20 @@ def write_version():
     p = subprocess.Popen(["git", "describe", "--dirty", "--tags", "--always"],
                          stdout=subprocess.PIPE)
     res = p.communicate()[0].strip().decode('utf-8')
-    with open("nifty6/git_version.py", "w") as file:
+    with open("nifty7/git_version.py", "w") as file:
         file.write('gitversion = "{}"\n'.format(res))
 
 
 write_version()
-exec(open('nifty6/version.py').read())
+exec(open('nifty7/version.py').read())
 
-setup(name="nifty6",
+setup(name="nifty7",
       version=__version__,
       author="Theo Steininger, Martin Reinecke",
       author_email="martin@mpa-garching.mpg.de",
       description="Numerical Information Field Theory",
       url="http://www.mpa-garching.mpg.de/ift/nifty/",
-      packages=find_packages(include=["nifty6", "nifty6.*"]),
+      packages=find_packages(include=["nifty7", "nifty7.*"]),
       zip_safe=True,
       license="GPLv3",
       setup_requires=['scipy>=1.4.1', 'numpy>=1.17'],
diff --git a/nifty6/__init__.py b/src/__init__.py
similarity index 98%
rename from nifty6/__init__.py
rename to src/__init__.py
index 36075af38f9e6be492c80ddb6b3945634bcc27a6..7614f1e6670c91d72d4e251915cb0e808e93f316 100644
--- a/nifty6/__init__.py
+++ b/src/__init__.py
@@ -102,4 +102,4 @@ from .operator_spectrum import operator_spectrum
 from .operator_tree_optimiser import optimise_operator
 
 # We deliberately don't set __all__ here, because we don't want people to do a
-# "from nifty6 import *"; that would swamp the global namespace.
+# "from nifty7 import *"; that would swamp the global namespace.
diff --git a/nifty6/domain_tuple.py b/src/domain_tuple.py
similarity index 99%
rename from nifty6/domain_tuple.py
rename to src/domain_tuple.py
index 84df832ca7a6a7c38a808f7048355fad9d113181..456c5d3924a3f367826959454e43b05f33a31ba0 100644
--- a/nifty6/domain_tuple.py
+++ b/src/domain_tuple.py
@@ -16,11 +16,12 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 from functools import reduce
-from . import utilities
-from .domains.domain import Domain
 
 import numpy as np
 
+from . import utilities
+from .domains.domain import Domain
+
 
 class DomainTuple(object):
     """Ordered sequence of Domain objects.
diff --git a/nifty6/domains/__init__.py b/src/domains/__init__.py
similarity index 100%
rename from nifty6/domains/__init__.py
rename to src/domains/__init__.py
diff --git a/nifty6/domains/dof_space.py b/src/domains/dof_space.py
similarity index 100%
rename from nifty6/domains/dof_space.py
rename to src/domains/dof_space.py
diff --git a/nifty6/domains/domain.py b/src/domains/domain.py
similarity index 99%
rename from nifty6/domains/domain.py
rename to src/domains/domain.py
index a5bed83066c2ad43797d58a95b2d12391aed9179..92ea0201766c07144b7989071be2af8c318dbb9c 100644
--- a/nifty6/domains/domain.py
+++ b/src/domains/domain.py
@@ -16,6 +16,7 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 from functools import reduce
+
 from ..utilities import NiftyMeta
 
 
diff --git a/nifty6/domains/gl_space.py b/src/domains/gl_space.py
similarity index 98%
rename from nifty6/domains/gl_space.py
rename to src/domains/gl_space.py
index 8172af35135cf4a3fdd90a6b2e626b8fc0c9acf9..2239c3602f2385e061ddc6dab2d394c4548ab3e8 100644
--- a/nifty6/domains/gl_space.py
+++ b/src/domains/gl_space.py
@@ -24,7 +24,7 @@ class GLSpace(StructuredDomain):
     """Represents a 2-sphere with Gauss-Legendre pixelization.
 
     Its harmonic partner domain is the
-    :class:`~nifty6.domains.lm_space.LMSpace`.
+    :class:`~nifty7.domains.lm_space.LMSpace`.
 
     Parameters
     ----------
diff --git a/nifty6/domains/hp_space.py b/src/domains/hp_space.py
similarity index 98%
rename from nifty6/domains/hp_space.py
rename to src/domains/hp_space.py
index a07a4e5c8034409bd35ef1921018cc7b027f019e..7e35755cbc84c53b5a000f3c90da6bfae0826856 100644
--- a/nifty6/domains/hp_space.py
+++ b/src/domains/hp_space.py
@@ -24,7 +24,7 @@ class HPSpace(StructuredDomain):
     """Represents 2-sphere with HEALPix discretization.
 
     Its harmonic partner domain is the
-    :class:`~nifty6.domains.lm_space.LMSpace`.
+    :class:`~nifty7.domains.lm_space.LMSpace`.
 
     Parameters
     ----------
diff --git a/nifty6/domains/lm_space.py b/src/domains/lm_space.py
similarity index 97%
rename from nifty6/domains/lm_space.py
rename to src/domains/lm_space.py
index 3cad266e3c8c7a0815c6040a2ba3a14a7515ebd7..6a7446f887c3bbdeb5e9290851e2efa9da1803b9 100644
--- a/nifty6/domains/lm_space.py
+++ b/src/domains/lm_space.py
@@ -24,8 +24,8 @@ from .structured_domain import StructuredDomain
 class LMSpace(StructuredDomain):
     """Represents a set of spherical harmonic coefficients.
 
-    Its harmonic partner spaces are :class:`~nifty6.domains.hp_space.HPSpace`
-    and :class:`~nifty6.domains.gl_space.GLSpace`.
+    Its harmonic partner spaces are :class:`~nifty7.domains.hp_space.HPSpace`
+    and :class:`~nifty7.domains.gl_space.GLSpace`.
 
     Parameters
     ----------
@@ -152,7 +152,7 @@ class LMSpace(StructuredDomain):
         return self._mmax
 
     def get_default_codomain(self):
-        """Returns a :class:`~nifty6.domains.gl_space.GLSpace` object, which is
+        """Returns a :class:`~nifty7.domains.gl_space.GLSpace` object, which is
         capable of storing an accurate representation of data residing on
         `self`.
 
diff --git a/nifty6/domains/power_space.py b/src/domains/power_space.py
similarity index 100%
rename from nifty6/domains/power_space.py
rename to src/domains/power_space.py
diff --git a/nifty6/domains/rg_space.py b/src/domains/rg_space.py
similarity index 99%
rename from nifty6/domains/rg_space.py
rename to src/domains/rg_space.py
index 4adc200c63b7137a36d3c452dbf162486abd6d18..20f2adc7dbb9eda1d3040e582e032ed5d4378dd7 100644
--- a/nifty6/domains/rg_space.py
+++ b/src/domains/rg_space.py
@@ -16,6 +16,7 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 from functools import reduce
+
 import numpy as np
 
 from ..field import Field
diff --git a/nifty6/domains/structured_domain.py b/src/domains/structured_domain.py
similarity index 100%
rename from nifty6/domains/structured_domain.py
rename to src/domains/structured_domain.py
diff --git a/nifty6/domains/unstructured_domain.py b/src/domains/unstructured_domain.py
similarity index 95%
rename from nifty6/domains/unstructured_domain.py
rename to src/domains/unstructured_domain.py
index 60029938f24c0716c2872a87d7d20f06cbaf91ba..218dfe940d287371b71deb9930b1927ef85b0000 100644
--- a/nifty6/domains/unstructured_domain.py
+++ b/src/domains/unstructured_domain.py
@@ -16,11 +16,12 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 from functools import reduce
+
 from .domain import Domain
 
 
 class UnstructuredDomain(Domain):
-    """A :class:`~nifty6.domains.domain.Domain` subclass for spaces with no
+    """A :class:`~nifty7.domains.domain.Domain` subclass for spaces with no
     associated geometry.
 
     Typically used for data spaces.
diff --git a/nifty6/extra.py b/src/extra.py
similarity index 72%
rename from nifty6/extra.py
rename to src/extra.py
index ae8da3e7f73b7a9f487b71d53cf1adcb73ba114f..052f7a37a2c7cb6c2a7f0a498738b48b2d864815 100644
--- a/nifty6/extra.py
+++ b/src/extra.py
@@ -15,6 +15,8 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
+from itertools import combinations
+
 import numpy as np
 from numpy.testing import assert_
 
@@ -23,13 +25,103 @@ from .field import Field
 from .linearization import Linearization
 from .multi_domain import MultiDomain
 from .multi_field import MultiField
+from .operators.energy_operators import EnergyOperator
 from .operators.linear_operator import LinearOperator
-from .sugar import from_random
+from .operators.operator import Operator
+from .sugar import from_random, makeDomain
 
-__all__ = ["consistency_check", "check_jacobian_consistency",
+__all__ = ["check_linear_operator", "check_operator",
            "assert_allclose"]
 
 
+def check_linear_operator(op, domain_dtype=np.float64, target_dtype=np.float64,
+                          atol=1e-12, rtol=1e-12, only_r_linear=False):
+    """
+    Checks an operator for algebraic consistency of its capabilities.
+
+    Checks whether times(), adjoint_times(), inverse_times() and
+    adjoint_inverse_times() (if in capability list) is implemented
+    consistently. Additionally, it checks whether the operator is linear.
+
+    Parameters
+    ----------
+    op : LinearOperator
+        Operator which shall be checked.
+    domain_dtype : dtype
+        The data type of the random vectors in the operator's domain. Default
+        is `np.float64`.
+    target_dtype : dtype
+        The data type of the random vectors in the operator's target. Default
+        is `np.float64`.
+    atol : float
+        Absolute tolerance for the check. If rtol is specified,
+        then satisfying any tolerance will let the check pass.
+        Default: 0.
+    rtol : float
+        Relative tolerance for the check. If atol is specified,
+        then satisfying any tolerance will let the check pass.
+        Default: 0.
+    only_r_linear: bool
+        set to True if the operator is only R-linear, not C-linear.
+        This will relax the adjointness test accordingly.
+    """
+    if not isinstance(op, LinearOperator):
+        raise TypeError('This test tests only linear operators.')
+    _domain_check_linear(op, domain_dtype)
+    _domain_check_linear(op.adjoint, target_dtype)
+    _domain_check_linear(op.inverse, target_dtype)
+    _domain_check_linear(op.adjoint.inverse, domain_dtype)
+    _check_linearity(op, domain_dtype, atol, rtol)
+    _check_linearity(op.adjoint, target_dtype, atol, rtol)
+    _check_linearity(op.inverse, target_dtype, atol, rtol)
+    _check_linearity(op.adjoint.inverse, domain_dtype, atol, rtol)
+    _full_implementation(op, domain_dtype, target_dtype, atol, rtol,
+                         only_r_linear)
+    _full_implementation(op.adjoint, target_dtype, domain_dtype, atol, rtol,
+                         only_r_linear)
+    _full_implementation(op.inverse, target_dtype, domain_dtype, atol, rtol,
+                         only_r_linear)
+    _full_implementation(op.adjoint.inverse, domain_dtype, target_dtype, atol,
+                         rtol, only_r_linear)
+
+
+def check_operator(op, loc, tol=1e-12, ntries=100, perf_check=True,
+                   only_r_differentiable=True, metric_sampling=True):
+    """
+    Performs various checks of the implementation of linear and nonlinear
+    operators.
+
+    Computes the Jacobian with finite differences and compares it to the
+    implemented Jacobian.
+
+    Parameters
+    ----------
+    op : Operator
+        Operator which shall be checked.
+    loc : Field or MultiField
+        An Field or MultiField instance which has the same domain
+        as op. The location at which the gradient is checked
+    tol : float
+        Tolerance for the check.
+    perf_check : Boolean
+        Do performance check. May be disabled for very unimportant operators.
+    only_r_differentiable : Boolean
+        Jacobians of C-differentiable operators need to be C-linear.
+        Default: True
+    metric_sampling: Boolean
+        If op is an EnergyOperator, metric_sampling determines whether the
+        test shall try to sample from the metric or not.
+    """
+    if not isinstance(op, Operator):
+        raise TypeError('This test tests only linear operators.')
+    _domain_check_nonlinear(op, loc)
+    _performance_check(op, loc, bool(perf_check))
+    _linearization_value_consistency(op, loc)
+    _jac_vs_finite_differences(op, loc, np.sqrt(tol), ntries, only_r_differentiable)
+    _check_nontrivial_constant(op, loc, tol, ntries, only_r_differentiable,
+                               metric_sampling)
+
+
 def assert_allclose(f1, f2, atol, rtol):
     if isinstance(f1, Field):
         return np.testing.assert_allclose(f1.val, f2.val, atol=atol, rtol=rtol)
@@ -37,6 +129,27 @@ def assert_allclose(f1, f2, atol, rtol):
         assert_allclose(val, f2[key], atol=atol, rtol=rtol)
 
 
+def assert_equal(f1, f2):
+    if isinstance(f1, Field):
+        return np.testing.assert_equal(f1.val, f2.val)
+    for key, val in f1.items():
+        assert_equal(val, f2[key])
+
+
+def _nozero(fld):
+    if isinstance(fld, Field):
+        return np.testing.assert_((fld != 0).s_all())
+    for val in fld.values():
+        _nozero(val)
+
+
+def _allzero(fld):
+    if isinstance(fld, Field):
+        return np.testing.assert_((fld == 0.).s_all())
+    for val in fld.values():
+        _allzero(val)
+
+
 def _adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol,
                             only_r_linear):
     needed_cap = op.TIMES | op.ADJOINT_TIMES
@@ -83,7 +196,8 @@ def _check_linearity(op, domain_dtype, atol, rtol):
     assert_allclose(val1, val2, atol=atol, rtol=rtol)
 
 
-def _actual_domain_check_linear(op, domain_dtype=None, inp=None):
+def _domain_check_linear(op, domain_dtype=None, inp=None):
+    _domain_check(op)
     needed_cap = op.TIMES
     if (op.capability & needed_cap) != needed_cap:
         return
@@ -95,8 +209,9 @@ def _actual_domain_check_linear(op, domain_dtype=None, inp=None):
     assert_(op(inp).domain is op.target)
 
 
-def _actual_domain_check_nonlinear(op, loc):
-    assert isinstance(loc, (Field, MultiField))
+def _domain_check_nonlinear(op, loc):
+    _domain_check(op)
+    assert_(isinstance(loc, (Field, MultiField)))
     assert_(loc.domain is op.domain)
     for wm in [False, True]:
         lin = Linearization.make_var(loc, wm)
@@ -111,8 +226,8 @@ def _actual_domain_check_nonlinear(op, loc):
         assert_(reslin.jac.domain is reslin.domain)
         assert_(reslin.jac.target is reslin.target)
         assert_(lin.want_metric == reslin.want_metric)
-        _actual_domain_check_linear(reslin.jac, inp=loc)
-        _actual_domain_check_linear(reslin.jac.adjoint, inp=reslin.jac(loc))
+        _domain_check_linear(reslin.jac, inp=loc)
+        _domain_check_linear(reslin.jac.adjoint, inp=reslin.jac(loc))
         if reslin.metric is not None:
             assert_(reslin.metric.domain is reslin.metric.target)
             assert_(reslin.metric.domain is op.domain)
@@ -164,58 +279,6 @@ def _performance_check(op, pos, raise_on_fail):
                 raise RuntimeError(s)
 
 
-def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
-                      atol=0, rtol=1e-7, only_r_linear=False):
-    """
-    Checks an operator for algebraic consistency of its capabilities.
-
-    Checks whether times(), adjoint_times(), inverse_times() and
-    adjoint_inverse_times() (if in capability list) is implemented
-    consistently. Additionally, it checks whether the operator is linear.
-
-    Parameters
-    ----------
-    op : LinearOperator
-        Operator which shall be checked.
-    domain_dtype : dtype
-        The data type of the random vectors in the operator's domain. Default
-        is `np.float64`.
-    target_dtype : dtype
-        The data type of the random vectors in the operator's target. Default
-        is `np.float64`.
-    atol : float
-        Absolute tolerance for the check. If rtol is specified,
-        then satisfying any tolerance will let the check pass.
-        Default: 0.
-    rtol : float
-        Relative tolerance for the check. If atol is specified,
-        then satisfying any tolerance will let the check pass.
-        Default: 0.
-    only_r_linear: bool
-        set to True if the operator is only R-linear, not C-linear.
-        This will relax the adjointness test accordingly.
-    """
-    if not isinstance(op, LinearOperator):
-        raise TypeError('This test tests only linear operators.')
-    _domain_check(op)
-    _actual_domain_check_linear(op, domain_dtype)
-    _actual_domain_check_linear(op.adjoint, target_dtype)
-    _actual_domain_check_linear(op.inverse, target_dtype)
-    _actual_domain_check_linear(op.adjoint.inverse, domain_dtype)
-    _check_linearity(op, domain_dtype, atol, rtol)
-    _check_linearity(op.adjoint, target_dtype, atol, rtol)
-    _check_linearity(op.inverse, target_dtype, atol, rtol)
-    _check_linearity(op.adjoint.inverse, domain_dtype, atol, rtol)
-    _full_implementation(op, domain_dtype, target_dtype, atol, rtol,
-                         only_r_linear)
-    _full_implementation(op.adjoint, target_dtype, domain_dtype, atol, rtol,
-                         only_r_linear)
-    _full_implementation(op.inverse, target_dtype, domain_dtype, atol, rtol,
-                         only_r_linear)
-    _full_implementation(op.adjoint.inverse, domain_dtype, target_dtype, atol,
-                         rtol, only_r_linear)
-
-
 def _get_acceptable_location(op, loc, lin):
     if not np.isfinite(lin.val.s_sum()):
         raise ValueError('Initial value must be finite')
@@ -248,34 +311,48 @@ def _linearization_value_consistency(op, loc):
         assert_allclose(fld0, fld1, 0, 1e-7)
 
 
-def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100, perf_check=True,
-        only_r_differentiable=True):
-    """
-    Checks the Jacobian of an operator against its finite difference
-    approximation.
-
-    Computes the Jacobian with finite differences and compares it to the
-    implemented Jacobian.
-
-    Parameters
-    ----------
-    op : Operator
-        Operator which shall be checked.
-    loc : Field or MultiField
-        An Field or MultiField instance which has the same domain
-        as op. The location at which the gradient is checked
-    tol : float
-        Tolerance for the check.
-    perf_check : Boolean
-        Do performance check. May be disabled for very unimportant operators.
-    only_r_differentiable : Boolean
-        Jacobians of C-differentiable operators need to be C-linear. 
-        Default: True
-    """
-    _domain_check(op)
-    _actual_domain_check_nonlinear(op, loc)
-    _performance_check(op, loc, bool(perf_check))
-    _linearization_value_consistency(op, loc)
+def _check_nontrivial_constant(op, loc, tol, ntries, only_r_differentiable,
+                               metric_sampling):
+    return  # FIXME
+    # Assumes that the operator is not constant
+    if isinstance(op.domain, DomainTuple):
+        return
+    keys = op.domain.keys()
+    for ll in range(0, len(keys)):
+        for cstkeys in combinations(keys, ll):
+            cstdom, vardom = {}, {}
+            for kk, dd in op.domain.items():
+                if kk in cstkeys:
+                    cstdom[kk] = dd
+                else:
+                    vardom[kk] = dd
+            cstdom, vardom = makeDomain(cstdom), makeDomain(vardom)
+            cstloc = loc.extract(cstdom)
+
+            val0 = op(loc)
+            _, op0 = op.simplify_for_constant_input(cstloc)
+            val1 = op0(loc)
+            # MR FIXME: This tests something we don't promise!
+#            val2 = op0(loc.unite(cstloc))
+#            assert_equal(val1, val2)
+            assert_equal(val0, val1)
+
+            lin = Linearization.make_var(loc, want_metric=True)
+            oplin = op0(lin)
+            if isinstance(op, EnergyOperator):
+                _allzero(oplin.gradient.extract(cstdom))
+            # MR FIXME: This tests something we don't promise!
+#            _allzero(oplin.jac(from_random(cstdom).unite(full(vardom, 0))))
+
+            if isinstance(op, EnergyOperator) and metric_sampling:
+                samp0 = oplin.metric.draw_sample()
+                _allzero(samp0.extract(cstdom))
+                _nozero(samp0.extract(vardom))
+
+            _jac_vs_finite_differences(op0, loc, np.sqrt(tol), ntries, only_r_differentiable)
+
+
+def _jac_vs_finite_differences(op, loc, tol, ntries, only_r_differentiable):
     for _ in range(ntries):
         lin = op(Linearization.make_var(loc))
         loc2, lin2 = _get_acceptable_location(op, loc, lin)
@@ -300,8 +377,6 @@ def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100, perf_check=True,
             print(hist)
             raise ValueError("gradient and value seem inconsistent")
         loc = locnext
-
-        ddtype = loc.values()[0].dtype if isinstance(loc, MultiField) else loc.dtype
-        tdtype = dirder.values()[0].dtype if isinstance(dirder, MultiField) else dirder.dtype
-        consistency_check(linmid.jac, domain_dtype=ddtype, target_dtype=tdtype,
-                          only_r_linear=only_r_differentiable)
+        check_linear_operator(linmid.jac, domain_dtype=loc.dtype,
+                              target_dtype=dirder.dtype,
+                              only_r_linear=only_r_differentiable)
diff --git a/nifty6/fft.py b/src/fft.py
similarity index 100%
rename from nifty6/fft.py
rename to src/fft.py
diff --git a/nifty6/field.py b/src/field.py
similarity index 99%
rename from nifty6/field.py
rename to src/field.py
index 1cd9253841ba80126bd49e2c0e956316e1e6d123..749a5b2e858429a69562b62ba1859ab6b5436a02 100644
--- a/nifty6/field.py
+++ b/src/field.py
@@ -16,6 +16,7 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 from functools import reduce
+
 import numpy as np
 
 from . import utilities
@@ -656,10 +657,10 @@ class Field(Operator):
         return np.sqrt(self.s_var())
 
     def __repr__(self):
-        return "<nifty6.Field>"
+        return "<nifty7.Field>"
 
     def __str__(self):
-        return "nifty6.Field instance\n- domain      = " + \
+        return "nifty7.Field instance\n- domain      = " + \
                self._domain.__str__() + \
                "\n- val         = " + repr(self._val)
 
diff --git a/nifty6/library/__init__.py b/src/library/__init__.py
similarity index 100%
rename from nifty6/library/__init__.py
rename to src/library/__init__.py
diff --git a/nifty6/library/adjust_variances.py b/src/library/adjust_variances.py
similarity index 100%
rename from nifty6/library/adjust_variances.py
rename to src/library/adjust_variances.py
diff --git a/nifty6/library/correlated_fields.py b/src/library/correlated_fields.py
similarity index 98%
rename from nifty6/library/correlated_fields.py
rename to src/library/correlated_fields.py
index 5b0c22f492523acc347b725b9fa2b92fec86ea13..886484066b9a574bb2b69c1a51f464c3012c88c0 100644
--- a/nifty6/library/correlated_fields.py
+++ b/src/library/correlated_fields.py
@@ -21,6 +21,7 @@ from operator import mul
 
 import numpy as np
 
+from .. import utilities
 from ..domain_tuple import DomainTuple
 from ..domains.power_space import PowerSpace
 from ..domains.unstructured_domain import UnstructuredDomain
@@ -39,7 +40,6 @@ from ..operators.simple_linear_operators import ducktape
 from ..operators.normal_operators import NormalTransform, LognormalTransform
 from ..probing import StatCalculator
 from ..sugar import full, makeDomain, makeField, makeOp
-from .. import utilities
 
 
 def _log_k_lengths(pspace):
@@ -312,11 +312,11 @@ class CorrelatedFieldMaker:
     around this offset is parametrized.
 
     The resulting correlated field model operator has a
-    :class:`~nifty6.multi_domain.MultiDomain` as its domain and
+    :class:`~nifty7.multi_domain.MultiDomain` as its domain and
     expects its input values to be univariately gaussian.
 
     The target of the constructed operator will be a
-    :class:`~nifty6.domain_tuple.DomainTuple` containing the
+    :class:`~nifty7.domain_tuple.DomainTuple` containing the
     `target_subdomains` of the added fluctuations in the order of
     the `add_fluctuations` calls.
 
@@ -420,8 +420,8 @@ class CorrelatedFieldMaker:
 
         Parameters
         ----------
-        target_subdomain : :class:`~nifty6.domain.Domain`, \
-                           :class:`~nifty6.domain_tuple.DomainTuple`
+        target_subdomain : :class:`~nifty7.domain.Domain`, \
+                           :class:`~nifty7.domain_tuple.DomainTuple`
             Target subdomain on which the correlation structure defined
             in this call should hold.
         fluctuations_{mean,stddev} : float
@@ -446,8 +446,8 @@ class CorrelatedFieldMaker:
             field model and the second one is used for the third field model.
             *If not given*, use the same power spectrum model for all
             constructed field models.
-        harmonic_partner : :class:`~nifty6.domain.Domain`, \
-                           :class:`~nifty6.domain_tuple.DomainTuple`
+        harmonic_partner : :class:`~nifty7.domain.Domain`, \
+                           :class:`~nifty7.domain_tuple.DomainTuple`
             In which harmonic space to define the power spectrum
         """
         if harmonic_partner is None:
diff --git a/nifty6/library/dynamic_operator.py b/src/library/dynamic_operator.py
similarity index 100%
rename from nifty6/library/dynamic_operator.py
rename to src/library/dynamic_operator.py
diff --git a/nifty6/library/gridder.py b/src/library/gridder.py
similarity index 99%
rename from nifty6/library/gridder.py
rename to src/library/gridder.py
index 1af2d314c2186a8870254b732b2bf95f6ea3f13b..595ba25d07f859c4fab39d3fe26b61785f604fb2 100644
--- a/nifty6/library/gridder.py
+++ b/src/library/gridder.py
@@ -21,7 +21,7 @@ from ..domain_tuple import DomainTuple
 from ..domains.rg_space import RGSpace
 from ..domains.unstructured_domain import UnstructuredDomain
 from ..operators.linear_operator import LinearOperator
-from ..sugar import makeField, makeDomain
+from ..sugar import makeDomain, makeField
 
 
 class GridderMaker(object):
diff --git a/nifty6/library/light_cone_operator.py b/src/library/light_cone_operator.py
similarity index 100%
rename from nifty6/library/light_cone_operator.py
rename to src/library/light_cone_operator.py
diff --git a/nifty6/library/los_response.py b/src/library/los_response.py
similarity index 100%
rename from nifty6/library/los_response.py
rename to src/library/los_response.py
diff --git a/nifty6/library/special_distributions.py b/src/library/special_distributions.py
similarity index 100%
rename from nifty6/library/special_distributions.py
rename to src/library/special_distributions.py
index c06b15beca3076792e6b31d4d31dee368a8959f2..49b37ca61b94b3d50b5c0ceaa0d401e41c7c1150 100644
--- a/nifty6/library/special_distributions.py
+++ b/src/library/special_distributions.py
@@ -16,15 +16,15 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import numpy as np
-from scipy.stats import invgamma, norm
 from scipy.interpolate import CubicSpline
+from scipy.stats import invgamma, norm
 
+from .. import random
 from ..domain_tuple import DomainTuple
 from ..domains.unstructured_domain import UnstructuredDomain
 from ..field import Field
 from ..operators.operator import Operator
 from ..sugar import makeOp
-from .. import random
 
 
 def _f_on_np(f, arr):
diff --git a/nifty6/library/wiener_filter_curvature.py b/src/library/wiener_filter_curvature.py
similarity index 100%
rename from nifty6/library/wiener_filter_curvature.py
rename to src/library/wiener_filter_curvature.py
diff --git a/nifty6/linearization.py b/src/linearization.py
similarity index 100%
rename from nifty6/linearization.py
rename to src/linearization.py
index 281342afc06811b2992443d8cbcb122c53d527f1..f0bffe2588cf8260f24888d6c148dc8d74917680 100644
--- a/nifty6/linearization.py
+++ b/src/linearization.py
@@ -17,8 +17,8 @@
 
 import numpy as np
 
-from .sugar import makeOp
 from .operators.operator import Operator
+from .sugar import makeOp
 
 
 class Linearization(Operator):
diff --git a/nifty6/logger.py b/src/logger.py
similarity index 96%
rename from nifty6/logger.py
rename to src/logger.py
index 1754d08c2efd2f2210b3964d70337f9c28d8c74b..27c59c36f366e40c331b253b0322d0ffb71d393c 100644
--- a/nifty6/logger.py
+++ b/src/logger.py
@@ -18,7 +18,7 @@
 
 def _logger_init():
     import logging
-    res = logging.getLogger('NIFTy6')
+    res = logging.getLogger('NIFTy7')
     res.setLevel(logging.DEBUG)
     res.propagate = False
     ch = logging.StreamHandler()
diff --git a/nifty6/minimization/__init__.py b/src/minimization/__init__.py
similarity index 100%
rename from nifty6/minimization/__init__.py
rename to src/minimization/__init__.py
diff --git a/nifty6/minimization/conjugate_gradient.py b/src/minimization/conjugate_gradient.py
similarity index 98%
rename from nifty6/minimization/conjugate_gradient.py
rename to src/minimization/conjugate_gradient.py
index 6b67d14432aefdd44a804054b596830318c9f9c6..57d2171bbe4fe7cc90e6b93d8feb76d8cb0b55d5 100644
--- a/nifty6/minimization/conjugate_gradient.py
+++ b/src/minimization/conjugate_gradient.py
@@ -27,7 +27,7 @@ class ConjugateGradient(Minimizer):
 
     Parameters
     ----------
-    controller : :py:class:`nifty6.IterationController`
+    controller : :py:class:`nifty7.IterationController`
         Object that decides when to terminate the minimization.
     nreset : int
         every `nreset` CG steps the residual will be recomputed accurately
diff --git a/nifty6/minimization/descent_minimizers.py b/src/minimization/descent_minimizers.py
similarity index 100%
rename from nifty6/minimization/descent_minimizers.py
rename to src/minimization/descent_minimizers.py
diff --git a/nifty6/minimization/energy.py b/src/minimization/energy.py
similarity index 100%
rename from nifty6/minimization/energy.py
rename to src/minimization/energy.py
diff --git a/nifty6/minimization/energy_adapter.py b/src/minimization/energy_adapter.py
similarity index 69%
rename from nifty6/minimization/energy_adapter.py
rename to src/minimization/energy_adapter.py
index 2624151f7b5421fffdcad205e16bb40ac988842f..284d9fab08d25be2294dd45d26ad5cbc93cafc76 100644
--- a/nifty6/minimization/energy_adapter.py
+++ b/src/minimization/energy_adapter.py
@@ -11,12 +11,15 @@
 # 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
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
+import numpy as np
+
 from ..linearization import Linearization
 from ..minimization.energy import Energy
+from ..sugar import makeDomain
 
 
 class EnergyAdapter(Energy):
@@ -38,22 +41,34 @@ class EnergyAdapter(Energy):
         If True, the class will provide a `metric` property. This should only
         be enabled if it is required, because it will most likely consume
         additional resources. Default: False.
+    nanisinf : bool
+        If true, nan energies which can happen due to overflows in the forward
+        model are interpreted as inf. Thereby, the code does not crash on
+        these occaisions but rather the minimizer is told that the position it
+        has tried is not sensible.
     """
 
-    def __init__(self, position, op, constants=[], want_metric=False):
+    def __init__(self, position, op, constants=[], want_metric=False,
+                 nanisinf=False):
         super(EnergyAdapter, self).__init__(position)
         self._op = op
-        self._constants = constants
+        if len(constants) > 0:
+            dom = makeDomain({kk: vv for kk, vv in position.domain.items()
+                              if kk in constants})
+            _, self._op = op.simplify_for_constant_input(position.extract(dom))
         self._want_metric = want_metric
-        lin = Linearization.make_partial_var(position, constants, want_metric)
+        lin = Linearization.make_var(position, want_metric)
         tmp = self._op(lin)
         self._val = tmp.val.val[()]
         self._grad = tmp.gradient
         self._metric = tmp._metric
+        self._nanisinf = bool(nanisinf)
+        if self._nanisinf and np.isnan(self._val):
+            self._val = np.inf
 
     def at(self, position):
-        return EnergyAdapter(position, self._op, self._constants,
-                             self._want_metric)
+        return EnergyAdapter(position, self._op, want_metric=self._want_metric,
+                             nanisinf=self._nanisinf)
 
     @property
     def value(self):
diff --git a/nifty6/minimization/iteration_controllers.py b/src/minimization/iteration_controllers.py
similarity index 100%
rename from nifty6/minimization/iteration_controllers.py
rename to src/minimization/iteration_controllers.py
diff --git a/nifty6/minimization/line_search.py b/src/minimization/line_search.py
similarity index 99%
rename from nifty6/minimization/line_search.py
rename to src/minimization/line_search.py
index 7aa4af3276ca5c362249bb467423cd9ce1dce544..590cdb86252a1ba89ffb8b8d5c9227e9233699ef 100644
--- a/nifty6/minimization/line_search.py
+++ b/src/minimization/line_search.py
@@ -15,10 +15,11 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-from ..utilities import NiftyMeta
-from ..logger import logger
 import numpy as np
 
+from ..logger import logger
+from ..utilities import NiftyMeta
+
 
 class LineEnergy(object):
     """Evaluates an underlying Energy along a certain line direction.
diff --git a/nifty6/minimization/metric_gaussian_kl.py b/src/minimization/metric_gaussian_kl.py
similarity index 50%
rename from nifty6/minimization/metric_gaussian_kl.py
rename to src/minimization/metric_gaussian_kl.py
index baeb53539afee850e229662695c0a6d57416ec3c..46ca5ea260eb41e8616fd0702db92bf87328cbeb 100644
--- a/nifty6/minimization/metric_gaussian_kl.py
+++ b/src/minimization/metric_gaussian_kl.py
@@ -18,13 +18,13 @@
 import numpy as np
 
 from .. import random, utilities
-from ..field import Field
 from ..linearization import Linearization
+from ..logger import logger
 from ..multi_field import MultiField
 from ..operators.endomorphic_operator import EndomorphicOperator
 from ..operators.energy_operators import StandardHamiltonian
 from ..probing import approximation2endo
-from ..sugar import full, makeOp
+from ..sugar import makeOp
 from .energy import Energy
 
 
@@ -42,6 +42,11 @@ class _KLMetric(EndomorphicOperator):
         return self._KL._metric_sample(from_inverse)
 
 
+def _get_lo_hi(comm, n_samples):
+    ntask, rank, _ = utilities.get_MPI_params_from_comm(comm)
+    return utilities.shareRange(n_samples, ntask, rank)
+
+
 class MetricGaussianKL(Energy):
     """Provides the sampled Kullback-Leibler divergence between a distribution
     and a Metric Gaussian.
@@ -58,120 +63,129 @@ class MetricGaussianKL(Energy):
     true probability distribution the standard parametrization is assumed.
     The samples of this class can be distributed among MPI tasks.
 
-    Parameters
-    ----------
-    mean : Field
-        Mean of the Gaussian probability distribution.
-    hamiltonian : StandardHamiltonian
-        Hamiltonian of the approximated probability distribution.
-    n_samples : integer
-        Number of samples used to stochastically estimate the KL.
-    constants : list
-        List of parameter keys that are kept constant during optimization.
-        Default is no constants.
-    point_estimates : list
-        List of parameter keys for which no samples are drawn, but that are
-        (possibly) optimized for, corresponding to point estimates of these.
-        Default is to draw samples for the complete domain.
-    mirror_samples : boolean
-        Whether the negative of the drawn samples are also used,
-        as they are equally legitimate samples. If true, the number of used
-        samples doubles. Mirroring samples stabilizes the KL estimate as
-        extreme sample variation is counterbalanced. Default is False.
-    napprox : int
-        Number of samples for computing preconditioner for sampling. No
-        preconditioning is done by default.
-    comm : MPI communicator or None
-        If not None, samples will be distributed as evenly as possible
-        across this communicator. If `mirror_samples` is set, then a sample and
-        its mirror image will always reside on the same task.
-    nanisinf : bool
-        If true, nan energies which can happen due to overflows in the forward
-        model are interpreted as inf. Thereby, the code does not crash on
-        these occaisions but rather the minimizer is told that the position it
-        has tried is not sensible.
-    _local_samples : None
-        Only a parameter for internal uses. Typically not to be set by users.
-
-    Note
-    ----
-    The two lists `constants` and `point_estimates` are independent from each
-    other. It is possible to sample along domains which are kept constant
-    during minimization and vice versa.
+    Notes
+    -----
 
+    DomainTuples should never be created using the constructor, but rather
+    via the factory function :attr:`make`!
     See also
     --------
     `Metric Gaussian Variational Inference`, Jakob Knollmüller,
     Torsten A. Enßlin, `<https://arxiv.org/abs/1901.11033>`_
     """
-
-    def __init__(self, mean, hamiltonian, n_samples, constants=[],
-                 point_estimates=[], mirror_samples=False,
-                 napprox=0, comm=None, _local_samples=None,
-                 nanisinf=False):
+    def __init__(self, mean, hamiltonian, n_samples, mirror_samples, comm,
+                 local_samples, nanisinf, _callingfrommake=False):
+        if not _callingfrommake:
+            raise NotImplementedError
         super(MetricGaussianKL, self).__init__(mean)
-
-        if not isinstance(hamiltonian, StandardHamiltonian):
-            raise TypeError
-        if hamiltonian.domain is not mean.domain:
-            raise ValueError
-        if not isinstance(n_samples, int):
-            raise TypeError
-        self._constants = tuple(constants)
-        self._point_estimates = tuple(point_estimates)
-        self._mitigate_nans = nanisinf
-        if not isinstance(mirror_samples, bool):
-            raise TypeError
-
         self._hamiltonian = hamiltonian
-
         self._n_samples = int(n_samples)
+        self._mirror_samples = bool(mirror_samples)
         self._comm = comm
-        ntask, rank, _ = utilities.get_MPI_params_from_comm(self._comm)
-        self._lo, self._hi = utilities.shareRange(self._n_samples, ntask, rank)
+        self._local_samples = local_samples
+        self._nanisinf = bool(nanisinf)
 
-        self._mirror_samples = bool(mirror_samples)
-        self._n_eff_samples = self._n_samples
-        if self._mirror_samples:
-            self._n_eff_samples *= 2
-
-        if _local_samples is None:
-            met = hamiltonian(Linearization.make_partial_var(
-                mean, self._point_estimates, True)).metric
-            if napprox >= 1:
-                met._approximation = makeOp(approximation2endo(met, napprox))
-            _local_samples = []
-            sseq = random.spawn_sseq(self._n_samples)
-            for i in range(self._lo, self._hi):
-                with random.Context(sseq[i]):
-                    _local_samples.append(met.draw_sample(from_inverse=True))
-            _local_samples = tuple(_local_samples)
-        else:
-            if len(_local_samples) != self._hi-self._lo:
-                raise ValueError("# of samples mismatch")
-        self._local_samples = _local_samples
-        self._lin = Linearization.make_partial_var(mean, self._constants)
+        lin = Linearization.make_var(mean)
         v, g = [], []
         for s in self._local_samples:
-            tmp = self._hamiltonian(self._lin+s)
+            tmp = hamiltonian(lin+s)
             tv = tmp.val.val
             tg = tmp.gradient
-            if self._mirror_samples:
-                tmp = self._hamiltonian(self._lin-s)
+            if mirror_samples:
+                tmp = hamiltonian(lin-s)
                 tv = tv + tmp.val.val
                 tg = tg + tmp.gradient
             v.append(tv)
             g.append(tg)
-        self._val = utilities.allreduce_sum(v, self._comm)[()]/self._n_eff_samples
-        if np.isnan(self._val) and self._mitigate_nans:
+        self._val = utilities.allreduce_sum(v, self._comm)[()]/self.n_eff_samples
+        if np.isnan(self._val) and self._nanisinf:
             self._val = np.inf
-        self._grad = utilities.allreduce_sum(g, self._comm)/self._n_eff_samples
+        self._grad = utilities.allreduce_sum(g, self._comm)/self.n_eff_samples
+
+    @staticmethod
+    def make(mean, hamiltonian, n_samples, constants=[], point_estimates=[],
+             mirror_samples=False, napprox=0, comm=None, nanisinf=False):
+        """Return instance of :class:`MetricGaussianKL`.
+
+        Parameters
+        ----------
+        mean : Field
+            Mean of the Gaussian probability distribution.
+        hamiltonian : StandardHamiltonian
+            Hamiltonian of the approximated probability distribution.
+        n_samples : integer
+            Number of samples used to stochastically estimate the KL.
+        constants : list
+            List of parameter keys that are kept constant during optimization.
+            Default is no constants.
+        point_estimates : list
+            List of parameter keys for which no samples are drawn, but that are
+            (possibly) optimized for, corresponding to point estimates of these.
+            Default is to draw samples for the complete domain.
+        mirror_samples : boolean
+            Whether the negative of the drawn samples are also used,
+            as they are equally legitimate samples. If true, the number of used
+            samples doubles. Mirroring samples stabilizes the KL estimate as
+            extreme sample variation is counterbalanced. Default is False.
+        napprox : int
+            Number of samples for computing preconditioner for sampling. No
+            preconditioning is done by default.
+        comm : MPI communicator or None
+            If not None, samples will be distributed as evenly as possible
+            across this communicator. If `mirror_samples` is set, then a sample and
+            its mirror image will always reside on the same task.
+        nanisinf : bool
+            If true, nan energies which can happen due to overflows in the forward
+            model are interpreted as inf. Thereby, the code does not crash on
+            these occaisions but rather the minimizer is told that the position it
+            has tried is not sensible.
+
+        Note
+        ----
+        The two lists `constants` and `point_estimates` are independent from each
+        other. It is possible to sample along domains which are kept constant
+        during minimization and vice versa.
+        """
+
+        if not isinstance(hamiltonian, StandardHamiltonian):
+            raise TypeError
+        if hamiltonian.domain is not mean.domain:
+            raise ValueError
+        if not isinstance(n_samples, int):
+            raise TypeError
+        if not isinstance(mirror_samples, bool):
+            raise TypeError
+        if isinstance(mean, MultiField) and set(point_estimates) == set(mean.keys()):
+            raise RuntimeError(
+                'Point estimates for whole domain. Use EnergyAdapter instead.')
+        n_samples = int(n_samples)
+        mirror_samples = bool(mirror_samples)
+
+        if isinstance(mean, MultiField):
+            cstpos = mean.extract_by_keys(point_estimates)
+            _, ham_sampling = hamiltonian.simplify_for_constant_input(cstpos)
+        else:
+            ham_sampling = hamiltonian
+        met = ham_sampling(Linearization.make_var(mean, True)).metric
+        if napprox >= 1:
+            met._approximation = makeOp(approximation2endo(met, napprox))
+        local_samples = []
+        sseq = random.spawn_sseq(n_samples)
+        for i in range(*_get_lo_hi(comm, n_samples)):
+            with random.Context(sseq[i]):
+                local_samples.append(met.draw_sample(from_inverse=True))
+        local_samples = tuple(local_samples)
+
+        if isinstance(mean, MultiField):
+            _, hamiltonian = hamiltonian.simplify_for_constant_input(mean.extract_by_keys(constants))
+        return MetricGaussianKL(
+            mean, hamiltonian, n_samples, mirror_samples, comm, local_samples,
+            nanisinf, _callingfrommake=True)
 
     def at(self, position):
         return MetricGaussianKL(
-            position, self._hamiltonian, self._n_samples, self._constants,
-            self._point_estimates, self._mirror_samples, comm=self._comm,
-            _local_samples=self._local_samples, nanisinf=self._mitigate_nans)
+            position, self._hamiltonian, self._n_samples, self._mirror_samples,
+            self._comm, self._local_samples, self._nanisinf, True)
 
     @property
     def value(self):
@@ -182,14 +196,20 @@ class MetricGaussianKL(Energy):
         return self._grad
 
     def apply_metric(self, x):
-        lin = self._lin.with_want_metric()
+        lin = Linearization.make_var(self.position, want_metric=True)
         res = []
         for s in self._local_samples:
             tmp = self._hamiltonian(lin+s).metric(x)
             if self._mirror_samples:
                 tmp = tmp + self._hamiltonian(lin-s).metric(x)
             res.append(tmp)
-        return utilities.allreduce_sum(res, self._comm)/self._n_eff_samples
+        return utilities.allreduce_sum(res, self._comm)/self.n_eff_samples
+
+    @property
+    def n_eff_samples(self):
+        if self._mirror_samples:
+            return 2*self._n_samples
+        return self._n_samples
 
     @property
     def metric(self):
@@ -205,9 +225,10 @@ class MetricGaussianKL(Energy):
                     yield -s
         else:
             rank_lo_hi = [utilities.shareRange(self._n_samples, ntask, i) for i in range(ntask)]
+            lo, _ = _get_lo_hi(self._comm, self._n_samples)
             for itask, (l, h) in enumerate(rank_lo_hi):
                 for i in range(l, h):
-                    data = self._local_samples[i-self._lo] if rank == itask else None
+                    data = self._local_samples[i-lo] if rank == itask else None
                     s = self._comm.bcast(data, root=itask)
                     yield s
                     if self._mirror_samples:
@@ -216,7 +237,11 @@ class MetricGaussianKL(Energy):
     def _metric_sample(self, from_inverse=False):
         if from_inverse:
             raise NotImplementedError()
-        lin = self._lin.with_want_metric()
+        s = ('This draws from the Hamiltonian used for evaluation and does '
+             ' not take point_estimates into accout. Make sure that this '
+             'is your intended use.')
+        logger.warning(s)
+        lin = Linearization.make_var(self.position, True)
         samp = []
         sseq = random.spawn_sseq(self._n_samples)
         for i, v in enumerate(self._local_samples):
@@ -225,4 +250,4 @@ class MetricGaussianKL(Energy):
                 if self._mirror_samples:
                     tmp = tmp + self._hamiltonian(lin-v).metric.draw_sample(from_inverse=False)
                 samp.append(tmp)
-        return utilities.allreduce_sum(samp, self._comm)/self._n_eff_samples
+        return utilities.allreduce_sum(samp, self._comm)/self.n_eff_samples
diff --git a/nifty6/minimization/minimizer.py b/src/minimization/minimizer.py
similarity index 100%
rename from nifty6/minimization/minimizer.py
rename to src/minimization/minimizer.py
diff --git a/nifty6/minimization/nonlinear_cg.py b/src/minimization/nonlinear_cg.py
similarity index 100%
rename from nifty6/minimization/nonlinear_cg.py
rename to src/minimization/nonlinear_cg.py
diff --git a/nifty6/minimization/quadratic_energy.py b/src/minimization/quadratic_energy.py
similarity index 100%
rename from nifty6/minimization/quadratic_energy.py
rename to src/minimization/quadratic_energy.py
diff --git a/nifty6/minimization/scipy_minimizer.py b/src/minimization/scipy_minimizer.py
similarity index 100%
rename from nifty6/minimization/scipy_minimizer.py
rename to src/minimization/scipy_minimizer.py
diff --git a/nifty6/multi_domain.py b/src/multi_domain.py
similarity index 100%
rename from nifty6/multi_domain.py
rename to src/multi_domain.py
diff --git a/nifty6/multi_field.py b/src/multi_field.py
similarity index 97%
rename from nifty6/multi_field.py
rename to src/multi_field.py
index e880fa9f84945befe66743fb6124e106e711c008..582936aea70a9b7760a1fb4cef6c2e8edd094c71 100644
--- a/nifty6/multi_field.py
+++ b/src/multi_field.py
@@ -18,9 +18,9 @@
 import numpy as np
 
 from . import utilities
+from .domain_tuple import DomainTuple
 from .field import Field
 from .multi_domain import MultiDomain
-from .domain_tuple import DomainTuple
 from .operators.operator import Operator
 
 
@@ -250,11 +250,17 @@ class MultiField(Operator):
         return MultiField(subset,
                           tuple(self[key] for key in subset.keys()))
 
+    def extract_by_keys(self, keys):
+        dom = MultiDomain.make({kk: vv for kk, vv in self.domain.items() if kk in keys})
+        return self.extract(dom)
+
     def extract_part(self, subset):
         if subset is self._domain:
             return self
-        return MultiField.from_dict({key: self[key] for key in subset.keys()
-                                     if key in self})
+        dct = {key: self[key] for key in subset.keys() if key in self}
+        if len(dct) == 0:
+            return None
+        return MultiField.from_dict(dct)
 
     def unite(self, other):
         """Merges two MultiFields on potentially different MultiDomains.
diff --git a/nifty6/operator_spectrum.py b/src/operator_spectrum.py
similarity index 99%
rename from nifty6/operator_spectrum.py
rename to src/operator_spectrum.py
index 4c3ff20457877472a90f338f03716fe5761990d7..7954ca2a93977617dd64f1b339093efa9c86f074 100644
--- a/nifty6/operator_spectrum.py
+++ b/src/operator_spectrum.py
@@ -23,7 +23,7 @@ from .multi_domain import MultiDomain
 from .multi_field import MultiField
 from .operators.linear_operator import LinearOperator
 from .operators.sandwich_operator import SandwichOperator
-from .sugar import makeField, makeDomain
+from .sugar import makeDomain, makeField
 
 
 class _DomRemover(LinearOperator):
diff --git a/nifty6/operator_tree_optimiser.py b/src/operator_tree_optimiser.py
similarity index 98%
rename from nifty6/operator_tree_optimiser.py
rename to src/operator_tree_optimiser.py
index 61a5fef6bc7bc05eea5a8b5c9d623f85c13c0a82..a571054b71a0887c5e8f03809bb99b06ef59dabe 100644
--- a/nifty6/operator_tree_optimiser.py
+++ b/src/operator_tree_optimiser.py
@@ -15,9 +15,15 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-from .operators.operator import _OpChain, _OpSum, _OpProd
-from .sugar import domain_union
+from copy import deepcopy
+
+from numpy import allclose
+
+from .multi_field import MultiField
+from .operators.operator import _OpChain, _OpProd, _OpSum
 from .operators.simple_linear_operators import FieldAdapter
+from .sugar import domain_union, from_random
+
 
 def _optimise_operator(op):
     """
@@ -260,10 +266,6 @@ def _optimise_operator(op):
     return op
 
 
-from copy import deepcopy
-from .sugar import from_random
-from .multi_field import MultiField
-from numpy import allclose
 
 
 
@@ -290,7 +292,7 @@ def optimise_operator(op):
     -----
     Operators are compared only by id, so best results are achieved when the following code
 
-    >>> from nifty6 import UniformOperator, DomainTuple
+    >>> from nifty7 import UniformOperator, DomainTuple
     >>> uni1 = UniformOperator(DomainTuple.scalar_domain()
     >>> uni2 = UniformOperator(DomainTuple.scalar_domain()
     >>> op = (uni1 + uni2)*(uni1 + uni2)
diff --git a/nifty6/operators/__init__.py b/src/operators/__init__.py
similarity index 100%
rename from nifty6/operators/__init__.py
rename to src/operators/__init__.py
diff --git a/nifty6/operators/adder.py b/src/operators/adder.py
similarity index 100%
rename from nifty6/operators/adder.py
rename to src/operators/adder.py
diff --git a/nifty6/operators/block_diagonal_operator.py b/src/operators/block_diagonal_operator.py
similarity index 100%
rename from nifty6/operators/block_diagonal_operator.py
rename to src/operators/block_diagonal_operator.py
index 34f6ffca36fd8ea3e0099d8cedb8d5171cb3a7dc..a204a9be730a9bf4f149a3db7daa9a3074a4521b 100644
--- a/nifty6/operators/block_diagonal_operator.py
+++ b/src/operators/block_diagonal_operator.py
@@ -20,6 +20,7 @@ from ..multi_field import MultiField
 from .endomorphic_operator import EndomorphicOperator
 from .linear_operator import LinearOperator
 
+
 class BlockDiagonalOperator(EndomorphicOperator):
     """
     Parameters
@@ -45,7 +46,6 @@ class BlockDiagonalOperator(EndomorphicOperator):
                 else:
                     raise TypeError("LinearOperator expected")
 
-
     def apply(self, x, mode):
         self._check_input(x, mode)
         val = tuple(op.apply(v, mode=mode) if op is not None else v
diff --git a/nifty6/operators/chain_operator.py b/src/operators/chain_operator.py
similarity index 100%
rename from nifty6/operators/chain_operator.py
rename to src/operators/chain_operator.py
diff --git a/nifty6/operators/contraction_operator.py b/src/operators/contraction_operator.py
similarity index 100%
rename from nifty6/operators/contraction_operator.py
rename to src/operators/contraction_operator.py
diff --git a/nifty6/operators/convolution_operators.py b/src/operators/convolution_operators.py
similarity index 100%
rename from nifty6/operators/convolution_operators.py
rename to src/operators/convolution_operators.py
diff --git a/nifty6/operators/diagonal_operator.py b/src/operators/diagonal_operator.py
similarity index 100%
rename from nifty6/operators/diagonal_operator.py
rename to src/operators/diagonal_operator.py
diff --git a/nifty6/operators/distributors.py b/src/operators/distributors.py
similarity index 100%
rename from nifty6/operators/distributors.py
rename to src/operators/distributors.py
diff --git a/nifty6/operators/domain_tuple_field_inserter.py b/src/operators/domain_tuple_field_inserter.py
similarity index 100%
rename from nifty6/operators/domain_tuple_field_inserter.py
rename to src/operators/domain_tuple_field_inserter.py
diff --git a/nifty6/operators/einsum.py b/src/operators/einsum.py
similarity index 100%
rename from nifty6/operators/einsum.py
rename to src/operators/einsum.py
diff --git a/nifty6/operators/endomorphic_operator.py b/src/operators/endomorphic_operator.py
similarity index 100%
rename from nifty6/operators/endomorphic_operator.py
rename to src/operators/endomorphic_operator.py
diff --git a/nifty6/operators/energy_operators.py b/src/operators/energy_operators.py
similarity index 85%
rename from nifty6/operators/energy_operators.py
rename to src/operators/energy_operators.py
index 21ece3648ad781d20b154722d6906baf92e0068f..7137aa47a6cac1ea6dcf1326327f9ea467608bb4 100644
--- a/nifty6/operators/energy_operators.py
+++ b/src/operators/energy_operators.py
@@ -48,6 +48,10 @@ def _check_sampling_dtype(domain, dtypes):
     raise TypeError
 
 
+def _iscomplex(dtype):
+    return np.issubdtype(dtype, np.complexfloating)
+
+
 def _field_to_dtype(field):
     if isinstance(field, Field):
         dt = field.dtype
@@ -127,10 +131,10 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
     The covariance is assumed to be diagonal.
 
     .. math ::
-        E(s,D) = - \\log G(s, D) = 0.5 (s)^\\dagger D^{-1} (s) + 0.5 tr log(D),
+        E(s,D) = - \\log G(s, C) = 0.5 (s)^\\dagger C (s) - 0.5 tr log(C),
 
     an information energy for a Gaussian distribution with residual s and
-    diagonal covariance D.
+    inverse diagonal covariance C.
     The domain of this energy will be a MultiDomain with two keys,
     the target will be the scalar domain.
 
@@ -139,10 +143,10 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
     domain : Domain, DomainTuple, tuple of Domain
         domain of the residual and domain of the covariance diagonal.
 
-    residual : key
+    residual_key : key
         Residual key of the Gaussian.
 
-    inverse_covariance : key
+    inverse_covariance_key : key
         Inverse covariance diagonal key of the Gaussian.
 
     sampling_dtype : np.dtype
@@ -150,21 +154,67 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
     """
 
     def __init__(self, domain, residual_key, inverse_covariance_key, sampling_dtype):
-        self._r = str(residual_key)
-        self._icov = str(inverse_covariance_key)
+        self._kr = str(residual_key)
+        self._ki = str(inverse_covariance_key)
         dom = DomainTuple.make(domain)
-        self._domain = MultiDomain.make({self._r: dom, self._icov: dom})
-        self._sampling_dtype = sampling_dtype
-        _check_sampling_dtype(self._domain, sampling_dtype)
+        self._domain = MultiDomain.make({self._kr: dom, self._ki: dom})
+        self._dt = {self._kr: sampling_dtype, self._ki: np.float64}
+        _check_sampling_dtype(self._domain, self._dt)
+        self._cplx = _iscomplex(sampling_dtype)
+
+    def apply(self, x):
+        self._check_input(x)
+        r, i = x[self._kr], x[self._ki]
+        if self._cplx:
+            res = 0.5*r.vdot(r*i.real).real - i.ptw("log").sum()
+        else:
+            res = 0.5*(r.vdot(r*i) - i.ptw("log").sum())
+        if not x.want_metric:
+            return res
+        met = i.val if self._cplx else 0.5*i.val
+        met = MultiField.from_dict({self._kr: i.val, self._ki: met**(-2)})
+        return res.add_metric(SamplingDtypeSetter(makeOp(met), self._dt))
+
+    def _simplify_for_constant_input_nontrivial(self, c_inp):
+        from .simplify_for_const import ConstantEnergyOperator
+        assert len(c_inp.keys()) == 1
+        key = c_inp.keys()[0]
+        assert key in self._domain.keys()
+        cst = c_inp[key]
+        if key == self._kr:
+            res = _SpecialGammaEnergy(cst).ducktape(self._ki)
+        else:
+            dt = self._dt[self._kr]
+            res = GaussianEnergy(inverse_covariance=makeOp(cst),
+                                 sampling_dtype=dt).ducktape(self._kr)
+            trlog = cst.log().sum().val_rw()
+            if not _iscomplex(dt):
+                trlog /= 2
+            res = res + ConstantEnergyOperator(res.domain, -trlog)
+        res = res + ConstantEnergyOperator(self._domain, 0.)
+        assert res.domain is self.domain
+        assert res.target is self.target
+        return None, res
+
+
+class _SpecialGammaEnergy(EnergyOperator):
+    def __init__(self, residual):
+        self._domain = DomainTuple.make(residual.domain)
+        self._resi = residual
+        self._cplx = _iscomplex(self._resi.dtype)
+        self._scale = ScalingOperator(self._domain, 1 if self._cplx else .5)
 
     def apply(self, x):
         self._check_input(x)
-        res = 0.5*(x[self._r].vdot(x[self._r]*x[self._icov]).real - x[self._icov].ptw("log").sum())
+        r = self._resi
+        if self._cplx:
+            res = 0.5*(r*x.real).vdot(r).real - x.ptw("log").sum()
+        else:
+            res = 0.5*((r*x).vdot(r) - x.ptw("log").sum())
         if not x.want_metric:
             return res
-        mf = {self._r: x.val[self._icov], self._icov: .5*x.val[self._icov]**(-2)}
-        met = makeOp(MultiField.from_dict(mf))
-        return res.add_metric(SamplingDtypeSetter(met, self._sampling_dtype))
+        met = makeOp((self._scale(x.val))**(-2))
+        return res.add_metric(SamplingDtypeSetter(met, self._resi.dtype))
 
 
 class GaussianEnergy(EnergyOperator):
@@ -227,6 +277,7 @@ class GaussianEnergy(EnergyOperator):
                 if sampling_dtype != _field_to_dtype(self._mean):
                     raise ValueError("Sampling dtype and mean not compatible")
 
+        self._icov = inverse_covariance
         if inverse_covariance is None:
             self._op = Squared2NormOperator(self._domain).scale(0.5)
             self._met = ScalingOperator(self._domain, 1)
@@ -252,6 +303,10 @@ class GaussianEnergy(EnergyOperator):
             return res.add_metric(self._met)
         return res
 
+    def __repr__(self):
+        dom = '()' if isinstance(self.domain, DomainTuple) else self.domain.keys()
+        return f'GaussianEnergy {dom}'
+
 
 class PoissonianEnergy(EnergyOperator):
     """Computes likelihood Hamiltonians of expected count field constrained by
diff --git a/nifty6/operators/field_zero_padder.py b/src/operators/field_zero_padder.py
similarity index 100%
rename from nifty6/operators/field_zero_padder.py
rename to src/operators/field_zero_padder.py
diff --git a/nifty6/operators/harmonic_operators.py b/src/operators/harmonic_operators.py
similarity index 100%
rename from nifty6/operators/harmonic_operators.py
rename to src/operators/harmonic_operators.py
diff --git a/nifty6/operators/inversion_enabler.py b/src/operators/inversion_enabler.py
similarity index 100%
rename from nifty6/operators/inversion_enabler.py
rename to src/operators/inversion_enabler.py
diff --git a/nifty6/operators/linear_interpolation.py b/src/operators/linear_interpolation.py
similarity index 100%
rename from nifty6/operators/linear_interpolation.py
rename to src/operators/linear_interpolation.py
diff --git a/nifty6/operators/linear_operator.py b/src/operators/linear_operator.py
similarity index 100%
rename from nifty6/operators/linear_operator.py
rename to src/operators/linear_operator.py
diff --git a/nifty6/operators/mask_operator.py b/src/operators/mask_operator.py
similarity index 100%
rename from nifty6/operators/mask_operator.py
rename to src/operators/mask_operator.py
diff --git a/nifty6/operators/matrix_product_operator.py b/src/operators/matrix_product_operator.py
similarity index 99%
rename from nifty6/operators/matrix_product_operator.py
rename to src/operators/matrix_product_operator.py
index 66357de1b0248cc4b1fecf0b29dbfc5dcb93cbe0..b9fa34300eade7044939339e5b80c8678fe90398 100644
--- a/nifty6/operators/matrix_product_operator.py
+++ b/src/operators/matrix_product_operator.py
@@ -15,11 +15,12 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
+import numpy as np
+
+from .. import utilities
 from ..domain_tuple import DomainTuple
 from ..field import Field
 from .endomorphic_operator import EndomorphicOperator
-from .. import utilities
-import numpy as np
 
 
 class MatrixProductOperator(EndomorphicOperator):
diff --git a/nifty6/operators/normal_operators.py b/src/operators/normal_operators.py
similarity index 95%
rename from nifty6/operators/normal_operators.py
rename to src/operators/normal_operators.py
index 11263b47c356b4e014d5a2fa3e461b07afa71d9d..37e5078d73b886b70222f8a42490191356a41497 100644
--- a/nifty6/operators/normal_operators.py
+++ b/src/operators/normal_operators.py
@@ -41,7 +41,7 @@ def NormalTransform(mean, sigma, key, N_copies=0):
     N_copies : integer
         If == 0, target will be a scalar field.
         If >= 1, target will be an
-        :class:`~nifty6.unstructured_domain.UnstructuredDomain`.
+        :class:`~nifty7.unstructured_domain.UnstructuredDomain`.
     """
     if N_copies == 0:
         domain = DomainTuple.scalar_domain()
@@ -71,7 +71,7 @@ def LognormalTransform(mean, sigma, key, N_copies):
     N_copies : integer
         If == 0, target will be a scalar field.
         If >= 1, target will be an
-        :class:`~nifty6.unstructured_domain.UnstructuredDomain`.
+        :class:`~nifty7.unstructured_domain.UnstructuredDomain`.
     """
     logmean, logsigma = lognormal_moments(mean, sigma, N_copies)
     return NormalTransform(logmean, logsigma, key, N_copies).ptw("exp")
diff --git a/nifty6/operators/operator.py b/src/operators/operator.py
similarity index 86%
rename from nifty6/operators/operator.py
rename to src/operators/operator.py
index 034e940dff4f7d81f572f3c97de2b6d8f181c4b5..d4417a82146411533595f38dfd5a45c5816aedfa 100644
--- a/nifty6/operators/operator.py
+++ b/src/operators/operator.py
@@ -11,14 +11,16 @@
 # 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
+# Copyright(C) 2013-2020 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import numpy as np
 
-from ..utilities import NiftyMeta, indent
 from .. import pointwise
+from ..logger import logger
+from ..multi_domain import MultiDomain
+from ..utilities import NiftyMeta, indent
 
 
 class Operator(metaclass=NiftyMeta):
@@ -269,15 +271,39 @@ class Operator(metaclass=NiftyMeta):
         return self.__class__.__name__
 
     def simplify_for_constant_input(self, c_inp):
+        from .energy_operators import EnergyOperator
+        from .simplify_for_const import ConstantEnergyOperator, ConstantOperator
         if c_inp is None:
             return None, self
-        if c_inp.domain == self.domain:
-            op = _ConstantOperator(self.domain, self(c_inp))
+        dom = c_inp.domain
+        if isinstance(dom, MultiDomain) and len(dom) == 0:
+            return None, self
+
+        # Convention: If c_inp is MultiField, it needs to be defined on a
+        # subdomain of self._domain
+        if isinstance(self.domain, MultiDomain):
+            assert isinstance(dom, MultiDomain)
+            if set(c_inp.keys()) > set(self.domain.keys()):
+                raise ValueError
+
+        if dom is self.domain:
+            if isinstance(self, EnergyOperator):
+                op = ConstantEnergyOperator(self.domain, self(c_inp))
+            else:
+                op = ConstantOperator(self.domain, self(c_inp))
             return op(c_inp), op
+        if not isinstance(dom, MultiDomain):
+            raise RuntimeError
         return self._simplify_for_constant_input_nontrivial(c_inp)
 
     def _simplify_for_constant_input_nontrivial(self, c_inp):
-        return None, self
+        from .simplify_for_const import SlowPartialConstantOperator
+        s = ('SlowPartialConstantOperator used. You might want to consider'
+             ' implementing `_simplify_for_constant_input_nontrivial()` for'
+             ' this operator:')
+        logger.warning(s)
+        logger.warning(self.__repr__())
+        return None, self @ SlowPartialConstantOperator(self.domain, c_inp.keys())
 
     def ptw(self, op, *args, **kwargs):
         return _OpChain.make((_FunctionApplier(self.target, op, *args, **kwargs), self))
@@ -291,65 +317,6 @@ for f in pointwise.ptw_dict.keys():
     setattr(Operator, f, func(f))
 
 
-class _ConstCollector(object):
-    def __init__(self):
-        self._const = None
-        self._nc = set()
-
-    def mult(self, const, fulldom):
-        if const is None:
-            self._nc |= set(fulldom)
-        else:
-            self._nc |= set(fulldom) - set(const)
-            if self._const is None:
-                from ..multi_field import MultiField
-                self._const = MultiField.from_dict(
-                    {key: const[key] for key in const if key not in self._nc})
-            else:
-                from ..multi_field import MultiField
-                self._const = MultiField.from_dict(
-                    {key: self._const[key]*const[key]
-                     for key in const if key not in self._nc})
-
-    def add(self, const, fulldom):
-        if const is None:
-            self._nc |= set(fulldom.keys())
-        else:
-            from ..multi_field import MultiField
-            self._nc |= set(fulldom.keys()) - set(const.keys())
-            if self._const is None:
-                self._const = MultiField.from_dict(
-                    {key: const[key]
-                     for key in const.keys() if key not in self._nc})
-            else:
-                self._const = self._const.unite(const)
-                self._const = MultiField.from_dict(
-                    {key: self._const[key]
-                     for key in self._const if key not in self._nc})
-
-    @property
-    def constfield(self):
-        return self._const
-
-
-class _ConstantOperator(Operator):
-    def __init__(self, dom, output):
-        from ..sugar import makeDomain
-        self._domain = makeDomain(dom)
-        self._target = output.domain
-        self._output = output
-
-    def apply(self, x):
-        from .simple_linear_operators import NullOperator
-        self._check_input(x)
-        if x.jac is not None:
-            return x.new(self._output, NullOperator(self._domain, self._target))
-        return self._output
-
-    def __repr__(self):
-        return 'ConstantOperator <- {}'.format(self.domain.keys())
-
-
 class _FunctionApplier(Operator):
     def __init__(self, domain, funcname, *args, **kwargs):
         from ..sugar import makeDomain
@@ -447,16 +414,16 @@ class _OpProd(Operator):
         return lin1.new(lin1._val*lin2._val, jac)
 
     def _simplify_for_constant_input_nontrivial(self, c_inp):
+        from ..multi_domain import MultiDomain
+        from .simplify_for_const import ConstCollector
+
         f1, o1 = self._op1.simplify_for_constant_input(
             c_inp.extract_part(self._op1.domain))
         f2, o2 = self._op2.simplify_for_constant_input(
             c_inp.extract_part(self._op2.domain))
-
-        from ..multi_domain import MultiDomain
         if not isinstance(self._target, MultiDomain):
             return None, _OpProd(o1, o2)
-
-        cc = _ConstCollector()
+        cc = ConstCollector()
         cc.mult(f1, o1.target)
         cc.mult(f2, o2.target)
         return cc.constfield, _OpProd(o1, o2)
@@ -493,16 +460,16 @@ class _OpSum(Operator):
         return res
 
     def _simplify_for_constant_input_nontrivial(self, c_inp):
+        from ..multi_domain import MultiDomain
+        from .simplify_for_const import ConstCollector
+
         f1, o1 = self._op1.simplify_for_constant_input(
             c_inp.extract_part(self._op1.domain))
         f2, o2 = self._op2.simplify_for_constant_input(
             c_inp.extract_part(self._op2.domain))
-
-        from ..multi_domain import MultiDomain
         if not isinstance(self._target, MultiDomain):
             return None, _OpSum(o1, o2)
-
-        cc = _ConstCollector()
+        cc = ConstCollector()
         cc.add(f1, o1.target)
         cc.add(f2, o2.target)
         return cc.constfield, _OpSum(o1, o2)
diff --git a/nifty6/operators/operator_adapter.py b/src/operators/operator_adapter.py
similarity index 100%
rename from nifty6/operators/operator_adapter.py
rename to src/operators/operator_adapter.py
diff --git a/nifty6/operators/outer_product_operator.py b/src/operators/outer_product_operator.py
similarity index 100%
rename from nifty6/operators/outer_product_operator.py
rename to src/operators/outer_product_operator.py
diff --git a/nifty6/operators/partial_conjugate.py b/src/operators/partial_conjugate.py
similarity index 99%
rename from nifty6/operators/partial_conjugate.py
rename to src/operators/partial_conjugate.py
index 9cdf9aa9733ceaeb08ec94b91ba7ddd532efad5f..48e6ad0ccba31a68962268ee612a4217f93cdce0 100644
--- a/nifty6/operators/partial_conjugate.py
+++ b/src/operators/partial_conjugate.py
@@ -16,9 +16,10 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-from .endomorphic_operator import EndomorphicOperator
 from ..multi_domain import MultiDomain
 from ..multi_field import MultiField
+from .endomorphic_operator import EndomorphicOperator
+
 
 class PartialConjugate(EndomorphicOperator):
     """Perform partial conjugation of a :class:`MultiField`
diff --git a/nifty6/operators/regridding_operator.py b/src/operators/regridding_operator.py
similarity index 100%
rename from nifty6/operators/regridding_operator.py
rename to src/operators/regridding_operator.py
diff --git a/nifty6/operators/sampling_enabler.py b/src/operators/sampling_enabler.py
similarity index 100%
rename from nifty6/operators/sampling_enabler.py
rename to src/operators/sampling_enabler.py
diff --git a/nifty6/operators/sandwich_operator.py b/src/operators/sandwich_operator.py
similarity index 100%
rename from nifty6/operators/sandwich_operator.py
rename to src/operators/sandwich_operator.py
diff --git a/nifty6/operators/scaling_operator.py b/src/operators/scaling_operator.py
similarity index 100%
rename from nifty6/operators/scaling_operator.py
rename to src/operators/scaling_operator.py
diff --git a/nifty6/operators/selection_operators.py b/src/operators/selection_operators.py
similarity index 98%
rename from nifty6/operators/selection_operators.py
rename to src/operators/selection_operators.py
index 6c8714b1c3d91d9649f860c0af5c808f87f994ad..ccc92e5c06772c31c21754134ffe73557539f347 100644
--- a/nifty6/operators/selection_operators.py
+++ b/src/operators/selection_operators.py
@@ -17,6 +17,7 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import numpy as np
+
 from ..domain_tuple import DomainTuple
 from ..domains.rg_space import RGSpace
 from ..domains.unstructured_domain import UnstructuredDomain
@@ -101,10 +102,8 @@ class SliceOperator(LinearOperator):
         return Field.from_raw(self.domain, res)
 
     def __str__(self):
-        ss = (
-            f"{self.__class__.__name__}"
-            f"({self.domain.shape} -> {self.target.shape})"
-        )
+        ss = (f"{self.__class__.__name__}"
+              f"({self.domain.shape} -> {self.target.shape})")
         return ss
 
 
diff --git a/nifty6/operators/simple_linear_operators.py b/src/operators/simple_linear_operators.py
similarity index 94%
rename from nifty6/operators/simple_linear_operators.py
rename to src/operators/simple_linear_operators.py
index 4bea5d9725ee6f8b32fca75fe2725c9d9e079234..dfe5b9911c077edc079674101c49a7c899d91d4e 100644
--- a/nifty6/operators/simple_linear_operators.py
+++ b/src/operators/simple_linear_operators.py
@@ -173,16 +173,9 @@ class FieldAdapter(LinearOperator):
             return MultiField(self._tgt(mode), (x,))
 
     def __repr__(self):
-        s = 'FieldAdapter'
-        dom = isinstance(self._domain, MultiDomain)
-        tgt = isinstance(self._target, MultiDomain)
-        if dom and tgt:
-            s += ' {} <- {}'.format(self._target.keys(), self._domain.keys())
-        elif dom:
-            s += ' <- {}'.format(self._domain.keys())
-        elif tgt:
-            s += ' {} <-'.format(self._target.keys())
-        return s
+        dom = self.domain.keys() if isinstance(self.domain, MultiDomain) else '()'
+        tgt = self.target.keys() if isinstance(self.target, MultiDomain) else '()'
+        return f'{tgt} <- {dom}'
 
 
 class _SlowFieldAdapter(LinearOperator):
@@ -349,6 +342,17 @@ class NullOperator(LinearOperator):
         self._check_input(x, mode)
         return self._nullfield(self._tgt(mode))
 
+    def __repr__(self):
+        dom = self.domain.keys() if isinstance(self.domain, MultiDomain) else '()'
+        tgt = self.target.keys() if isinstance(self.target, MultiDomain) else '()'
+        return f'{tgt} <- NullOperator <- {dom}'
+
+    def draw_sample(self, from_inverse=False):
+        if self._domain is not self._target:
+            raise RuntimeError
+        from ..sugar import full
+        return full(self._domain, 0.)
+
 
 class PartialExtractor(LinearOperator):
     def __init__(self, domain, target):
@@ -373,3 +377,6 @@ class PartialExtractor(LinearOperator):
         res0 = MultiField.from_dict({key: x[key] for key in x.domain.keys()})
         res1 = MultiField.full(self._compldomain, 0.)
         return res0.unite(res1)
+
+    def __repr__(self):
+        return f'{self.target.keys()} <- {self.domain.keys()}'
diff --git a/src/operators/simplify_for_const.py b/src/operators/simplify_for_const.py
new file mode 100644
index 0000000000000000000000000000000000000000..99f4f830b7f870326898acb5486cb5e08e39fbf2
--- /dev/null
+++ b/src/operators/simplify_for_const.py
@@ -0,0 +1,125 @@
+# 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-2020 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+from ..multi_domain import MultiDomain
+from .block_diagonal_operator import BlockDiagonalOperator
+from .energy_operators import EnergyOperator
+from .operator import Operator
+from .scaling_operator import ScalingOperator
+from .simple_linear_operators import NullOperator
+
+
+class ConstCollector(object):
+    def __init__(self):
+        self._const = None  # MultiField on the part of the MultiDomain that could be constant
+        self._nc = set()  # NoConstant - set of keys that we know cannot be constant
+
+    def mult(self, const, fulldom):
+        if const is None:
+            self._nc |= set(fulldom.keys())
+        else:
+            from ..multi_field import MultiField
+            self._nc |= set(fulldom.keys()) - set(const.keys())
+            if self._const is None:
+                self._const = MultiField.from_dict(
+                    {key: const[key]
+                     for key in const.keys() if key not in self._nc})
+            else:  # we know that the domains are identical for products
+                self._const = MultiField.from_dict(
+                    {key: self._const[key]*const[key]
+                     for key in const.keys() if key not in self._nc})
+
+    def add(self, const, fulldom):
+        if const is None:
+            self._nc |= set(fulldom.keys())
+        else:
+            from ..multi_field import MultiField
+            self._nc |= set(fulldom.keys()) - set(const.keys())
+            self._const = const if self._const is None else self._const.unite(const)
+            self._const = MultiField.from_dict(
+                {key: const[key]
+                 for key in const.keys() if key not in self._nc})
+
+    @property
+    def constfield(self):
+        return self._const
+
+
+class ConstantOperator(Operator):
+    def __init__(self, dom, output):
+        from ..sugar import makeDomain
+        self._domain = makeDomain(dom)
+        self._target = output.domain
+        self._output = output
+
+    def apply(self, x):
+        from .simple_linear_operators import NullOperator
+        self._check_input(x)
+        if x.jac is not None:
+            return x.new(self._output, NullOperator(self._domain, self._target))
+        return self._output
+
+    def __repr__(self):
+        dom = self.domain.keys() if isinstance(self.domain, MultiDomain) else '()'
+        tgt = self.target.keys() if isinstance(self.target, MultiDomain) else '()'
+        return f'{tgt} <- ConstantOperator <- {dom}'
+
+
+class SlowPartialConstantOperator(Operator):
+    def __init__(self, domain, constant_keys):
+        from ..sugar import makeDomain
+        if not isinstance(domain, MultiDomain):
+            raise TypeError
+        if set(constant_keys) > set(domain.keys()) or len(constant_keys) == 0:
+            raise ValueError
+        self._keys = set(constant_keys) & set(domain.keys())
+        self._domain = self._target = makeDomain(domain)
+
+    def apply(self, x):
+        self._check_input(x)
+        if x.jac is None:
+            return x
+        jac = {kk: ScalingOperator(dd, 0 if kk in self._keys else 1)
+               for kk, dd in self._domain.items()}
+        return x.prepend_jac(BlockDiagonalOperator(x.jac.domain, jac))
+
+    def __repr__(self):
+        return f'SlowPartialConstantOperator ({self._keys})'
+
+
+class ConstantEnergyOperator(EnergyOperator):
+    def __init__(self, dom, output):
+        from ..sugar import makeDomain
+        from ..field import Field
+        self._domain = makeDomain(dom)
+        if not isinstance(output, Field):
+            output = Field.scalar(float(output))
+        if self.target is not output.domain:
+            raise TypeError
+        self._output = output
+
+    def apply(self, x):
+        self._check_input(x)
+        if x.jac is not None:
+            val = self._output
+            jac = NullOperator(self._domain, self._target)
+            met = NullOperator(self._domain, self._domain) if x.want_metric else None
+            return x.new(val, jac, met)
+        return self._output
+
+    def __repr__(self):
+        return 'ConstantEnergyOperator <- {}'.format(self.domain.keys())
diff --git a/nifty6/operators/sum_operator.py b/src/operators/sum_operator.py
similarity index 98%
rename from nifty6/operators/sum_operator.py
rename to src/operators/sum_operator.py
index e758d57d6b8a9f0be11401e56141fd3de0813af6..2a7fb7fb8e2a5a580370a7dd5570bcf4d13f84f3 100644
--- a/nifty6/operators/sum_operator.py
+++ b/src/operators/sum_operator.py
@@ -224,8 +224,8 @@ class SumOperator(LinearOperator):
                 fullop = op if fullop is None else fullop + op
             return None, fullop
 
-        from .operator import _ConstCollector
-        cc = _ConstCollector()
+        from .simplify_for_const import ConstCollector
+        cc = ConstCollector()
         fullop = None
         for tf, to, n in zip(f, o, self._neg):
             cc.add(tf, to.target)
diff --git a/nifty6/operators/value_inserter.py b/src/operators/value_inserter.py
similarity index 100%
rename from nifty6/operators/value_inserter.py
rename to src/operators/value_inserter.py
diff --git a/nifty6/plot.py b/src/plot.py
similarity index 100%
rename from nifty6/plot.py
rename to src/plot.py
diff --git a/nifty6/pointwise.py b/src/pointwise.py
similarity index 100%
rename from nifty6/pointwise.py
rename to src/pointwise.py
diff --git a/nifty6/probing.py b/src/probing.py
similarity index 99%
rename from nifty6/probing.py
rename to src/probing.py
index 46968ae8c2d26d8953947e338bac032dc22576b0..4f224cc651e4db4520308624bf6e5cd24cd6103f 100644
--- a/nifty6/probing.py
+++ b/src/probing.py
@@ -18,7 +18,7 @@
 from .multi_field import MultiField
 from .operators.endomorphic_operator import EndomorphicOperator
 from .operators.operator import Operator
-from .sugar import makeField, from_random
+from .sugar import from_random, makeField
 
 
 class StatCalculator(object):
diff --git a/nifty6/random.py b/src/random.py
similarity index 98%
rename from nifty6/random.py
rename to src/random.py
index 82b9d718f42c774be73c37da0d53eed3a19ab470..99f117c9229fc7003f5dcfd66b616cc200c77fc4 100644
--- a/nifty6/random.py
+++ b/src/random.py
@@ -21,7 +21,7 @@ Some remarks on NIFTy's treatment of random numbers
 NIFTy makes use of the `Generator` and `SeedSequence` classes introduced to
 `numpy.random` in numpy 1.17.
 
-On first load of the `nifty6.random` module, it creates a stack of
+On first load of the `nifty7.random` module, it creates a stack of
 `SeedSequence` objects which contains a single `SeedSequence` with a fixed seed,
 and also a stack of `Generator` objects, which contains a single generator
 derived from the above seed sequence. Without user intervention, this generator
@@ -87,7 +87,7 @@ def _fix_seed(seed):
 
 # Stack of SeedSequence objects. Will always start out with a well-defined
 # default. Users can change the "random seed" used by a calculation by pushing
-# a different SeedSequence before invoking any other nifty6.random calls
+# a different SeedSequence before invoking any other nifty7.random calls
 _sseq = [np.random.SeedSequence(_fix_seed(42))]
 # Stack of random number generators associated with _sseq.
 _rng = [np.random.default_rng(_sseq[-1])]
diff --git a/nifty6/sugar.py b/src/sugar.py
similarity index 99%
rename from nifty6/sugar.py
rename to src/sugar.py
index d384af8598828dd7fc27132a9f9e8bb2d4c418bc..b5cecfa8d3802e2a92d54972300cee1a3891f03f 100644
--- a/nifty6/sugar.py
+++ b/src/sugar.py
@@ -522,7 +522,7 @@ def calculate_position(operator, output):
     minimizer = NewtonCG(GradientNormController(iteration_limit=10, name='findpos'))
     for ii in range(3):
         logger.info(f'Start iteration {ii+1}/3')
-        kl = MetricGaussianKL(pos, H, 3, mirror_samples=True)
+        kl = MetricGaussianKL.make(pos, H, 3, mirror_samples=True)
         kl, _ = minimizer(kl)
         pos = kl.position
     return pos
diff --git a/nifty6/utilities.py b/src/utilities.py
similarity index 99%
rename from nifty6/utilities.py
rename to src/utilities.py
index ab9bac60b2cb1aab4ec475c9e3ab025c572ba50d..67f4948aeb14e3e07539e2ea08ee1b17b0f03134 100644
--- a/nifty6/utilities.py
+++ b/src/utilities.py
@@ -16,8 +16,8 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import collections
-from itertools import product
 from functools import reduce
+from itertools import product
 
 import numpy as np
 
@@ -395,7 +395,7 @@ def lognormal_moments(mean, sigma, N=0):
     """Calculates the parameters for a normal distribution `n(x)`
     such that `exp(n)(x)` has the mean and standard deviation given.
 
-    Used in :func:`~nifty6.normal_operators.LognormalTransform`."""
+    Used in :func:`~nifty7.normal_operators.LognormalTransform`."""
     mean, sigma = (value_reshaper(param, N) for param in (mean, sigma))
     if not np.all(mean > 0):
         raise ValueError("mean must be greater 0; got {!r}".format(mean))
diff --git a/nifty6/version.py b/src/version.py
similarity index 100%
rename from nifty6/version.py
rename to src/version.py
diff --git a/test/common.py b/test/common.py
index ba5224fe61a222cd9ba87f9ac5d6c6f46d78116d..0985bf92711ea122087484f5e5b7dd948d6c916f 100644
--- a/test/common.py
+++ b/test/common.py
@@ -27,10 +27,10 @@ def list2fixture(lst):
 
 
 def setup_function():
-    import nifty6 as ift
+    import nifty7 as ift
     ift.random.push_sseq_from_seed(42)
 
 
 def teardown_function():
-    import nifty6 as ift
+    import nifty7 as ift
     ift.random.pop_sseq()
diff --git a/test/test_energy_gradients.py b/test/test_energy_gradients.py
index 41e582a24c8f90507913b09fc826087e475a23f3..0dc0d0b0b69d6eddcf9334bda1f1ede194420d81 100644
--- a/test/test_energy_gradients.py
+++ b/test/test_energy_gradients.py
@@ -18,7 +18,7 @@
 import numpy as np
 import pytest
 
-import nifty6 as ift
+import nifty7 as ift
 
 from .common import list2fixture, setup_function, teardown_function
 
@@ -32,13 +32,13 @@ ntries = 10
 
 def test_gaussian(field):
     energy = ift.GaussianEnergy(domain=field.domain)
-    ift.extra.check_jacobian_consistency(energy, field)
+    ift.extra.check_operator(energy, field)
 
 
 def test_ScaledEnergy(field):
     icov = ift.ScalingOperator(field.domain, 1.2)
     energy = ift.GaussianEnergy(inverse_covariance=icov, sampling_dtype=np.float64)
-    ift.extra.check_jacobian_consistency(energy.scale(0.3), field)
+    ift.extra.check_operator(energy.scale(0.3), field)
 
     lin = ift.Linearization.make_var(field, want_metric=True)
     met1 = energy(lin).metric
@@ -54,17 +54,17 @@ def test_QuadraticFormOperator(field):
     op = ift.ScalingOperator(field.domain, 1.2)
     endo = ift.makeOp(op.draw_sample_with_dtype(dtype=np.float64))
     energy = ift.QuadraticFormOperator(endo)
-    ift.extra.check_jacobian_consistency(energy, field)
+    ift.extra.check_operator(energy, field)
 
 
 def test_studentt(field):
     if isinstance(field.domain, ift.MultiDomain):
         return
     energy = ift.StudentTEnergy(domain=field.domain, theta=.5)
-    ift.extra.check_jacobian_consistency(energy, field, tol=1e-6)
+    ift.extra.check_operator(energy, field)
     theta = ift.from_random(field.domain, 'normal').exp()
     energy = ift.StudentTEnergy(domain=field.domain, theta=theta)
-    ift.extra.check_jacobian_consistency(energy, field, tol=1e-6, ntries=ntries)
+    ift.extra.check_operator(energy, field, ntries=ntries)
 
 
 def test_hamiltonian_and_KL(field):
@@ -72,10 +72,10 @@ def test_hamiltonian_and_KL(field):
     space = field.domain
     lh = ift.GaussianEnergy(domain=space)
     hamiltonian = ift.StandardHamiltonian(lh)
-    ift.extra.check_jacobian_consistency(hamiltonian, field, ntries=ntries)
+    ift.extra.check_operator(hamiltonian, field, ntries=ntries)
     samps = [ift.from_random(space, 'normal') for i in range(2)]
     kl = ift.AveragedEnergy(hamiltonian, samps)
-    ift.extra.check_jacobian_consistency(kl, field, ntries=ntries)
+    ift.extra.check_operator(kl, field, ntries=ntries)
 
 
 def test_variablecovariancegaussian(field):
@@ -84,10 +84,19 @@ def test_variablecovariancegaussian(field):
     dc = {'a': field, 'b': field.ptw("exp")}
     mf = ift.MultiField.from_dict(dc)
     energy = ift.VariableCovarianceGaussianEnergy(field.domain, 'a', 'b', np.float64)
-    ift.extra.check_jacobian_consistency(energy, mf, tol=1e-6, ntries=ntries)
+    ift.extra.check_operator(energy, mf, ntries=ntries)
     energy(ift.Linearization.make_var(mf, want_metric=True)).metric.draw_sample()
 
 
+def test_specialgamma(field):
+    if isinstance(field.domain, ift.MultiDomain):
+        return
+    energy = ift.operators.energy_operators._SpecialGammaEnergy(field)
+    loc = ift.from_random(energy.domain).exp()
+    ift.extra.check_operator(energy, loc, ntries=ntries)
+    energy(ift.Linearization.make_var(loc, want_metric=True)).metric.draw_sample()
+
+
 def test_inverse_gamma(field):
     if isinstance(field.domain, ift.MultiDomain):
         return
@@ -96,7 +105,7 @@ def test_inverse_gamma(field):
     d = ift.random.current_rng().normal(10, size=space.shape)**2
     d = ift.Field(space, d)
     energy = ift.InverseGammaLikelihood(d)
-    ift.extra.check_jacobian_consistency(energy, field, tol=1e-5)
+    ift.extra.check_operator(energy, field, tol=1e-10)
 
 
 def testPoissonian(field):
@@ -107,7 +116,7 @@ def testPoissonian(field):
     d = ift.random.current_rng().poisson(120, size=space.shape)
     d = ift.Field(space, d)
     energy = ift.PoissonianEnergy(d)
-    ift.extra.check_jacobian_consistency(energy, field, tol=1e-6)
+    ift.extra.check_operator(energy, field)
 
 
 def test_bernoulli(field):
@@ -118,4 +127,4 @@ def test_bernoulli(field):
     d = ift.random.current_rng().binomial(1, 0.1, size=space.shape)
     d = ift.Field(space, d)
     energy = ift.BernoulliEnergy(d)
-    ift.extra.check_jacobian_consistency(energy, field, tol=1e-5)
+    ift.extra.check_operator(energy, field, tol=1e-10)
diff --git a/test/test_field.py b/test/test_field.py
index 7096cc5319f74eee471e1887a961ee38cabec8ce..f9c6e811421e9e2e63165c5528a4abe4dbcb0d99 100644
--- a/test/test_field.py
+++ b/test/test_field.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_allclose, assert_equal, assert_raises
 
-import nifty6 as ift
+import nifty7 as ift
 
 from .common import setup_function, teardown_function
 
diff --git a/test/test_gaussian_energy.py b/test/test_gaussian_energy.py
index 83894b3f27dfa75820a70aed85def412baaaa81d..a473dc302680972426b811850b76695d69871ad6 100644
--- a/test/test_gaussian_energy.py
+++ b/test/test_gaussian_energy.py
@@ -18,7 +18,7 @@
 import numpy as np
 import pytest
 
-import nifty6 as ift
+import nifty7 as ift
 
 from .common import setup_function, teardown_function
 
@@ -28,6 +28,7 @@ def _flat_PS(k):
 
 
 pmp = pytest.mark.parametrize
+ntries = 10
 
 
 @pmp('space', [ift.GLSpace(5),
@@ -69,5 +70,56 @@ def test_gaussian_energy(space, nonlinearity, noise, seed):
             N = None
 
         energy = ift.GaussianEnergy(d, N) @ d_model()
-        ift.extra.check_jacobian_consistency(
-            energy, xi0, ntries=10, tol=1e-6)
+        ift.extra.check_operator(
+            energy, xi0, ntries=ntries, tol=1e-6)
+
+
+@pmp('cplx', [True, False])
+def testgaussianenergy_compatibility(cplx):
+    dt = np.complex128 if cplx else np.float64
+    dom = ift.UnstructuredDomain(3)
+    e = ift.VariableCovarianceGaussianEnergy(dom, 'resi', 'icov', dt)
+    resi = ift.from_random(dom)
+    if cplx:
+        resi = resi + 1j*ift.from_random(dom)
+    loc0 = ift.MultiField.from_dict({'resi': resi})
+    loc1 = ift.MultiField.from_dict({'icov': ift.from_random(dom).exp()})
+    loc = loc0.unite(loc1)
+    val0 = e(loc).val
+
+    _, e0 = e.simplify_for_constant_input(loc0)
+    val1 = e0(loc).val
+    val2 = e0(loc.unite(loc0)).val
+    np.testing.assert_equal(val1, val2)
+    np.testing.assert_equal(val0, val1)
+
+    _, e1 = e.simplify_for_constant_input(loc1)
+    val1 = e1(loc).val
+    val2 = e1(loc.unite(loc1)).val
+    np.testing.assert_equal(val0, val1)
+    np.testing.assert_equal(val1, val2)
+
+    ift.extra.check_operator(e, loc, ntries=ntries)
+    ift.extra.check_operator(e0, loc, ntries=ntries, tol=1e-7)
+    ift.extra.check_operator(e1, loc, ntries=ntries)
+
+    # Test jacobian is zero
+    lin = ift.Linearization.make_var(loc, want_metric=True)
+    grad = e(lin).gradient.val
+    grad0 = e0(lin).gradient.val
+    grad1 = e1(lin).gradient.val
+    samp = e(lin).metric.draw_sample().val
+    samp0 = e0(lin).metric.draw_sample().val
+    samp1 = e1(lin).metric.draw_sample().val
+    np.testing.assert_equal(samp0['resi'], 0.)
+    np.testing.assert_equal(samp1['icov'], 0.)
+    np.testing.assert_equal(grad0['resi'], 0.)
+    np.testing.assert_equal(grad1['icov'], 0.)
+    np.testing.assert_(all(samp['resi'] != 0))
+    np.testing.assert_(all(samp['icov'] != 0))
+    np.testing.assert_(all(samp0['icov'] != 0))
+    np.testing.assert_(all(samp1['resi'] != 0))
+    np.testing.assert_(all(grad['resi'] != 0))
+    np.testing.assert_(all(grad['icov'] != 0))
+    np.testing.assert_(all(grad0['icov'] != 0))
+    np.testing.assert_(all(grad1['resi'] != 0))
diff --git a/test/test_kl.py b/test/test_kl.py
index 9bbebf51cc8ca308530a5facb8ed59b450e3bc19..b8f07a3f54101299adda6d2b8faf123f74a6dfa3 100644
--- a/test/test_kl.py
+++ b/test/test_kl.py
@@ -16,9 +16,9 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import pytest
-from numpy.testing import assert_, assert_allclose
+from numpy.testing import assert_, assert_allclose, assert_raises
 
-import nifty6 as ift
+import nifty7 as ift
 
 from .common import setup_function, teardown_function
 
@@ -44,13 +44,17 @@ def test_kl(constants, point_estimates, mirror_samples, mf):
     mean0 = ift.from_random(h.domain, 'normal')
 
     nsamps = 2
-    kl = ift.MetricGaussianKL(mean0,
-                              h,
-                              nsamps,
-                              constants=constants,
-                              point_estimates=point_estimates,
-                              mirror_samples=mirror_samples,
-                              napprox=0)
+    args = {'constants': constants,
+            'point_estimates': point_estimates,
+            'mirror_samples': mirror_samples,
+            'n_samples': nsamps,
+            'mean': mean0,
+            'hamiltonian': h}
+    if isinstance(mean0, ift.MultiField) and set(point_estimates) == set(mean0.keys()):
+        with assert_raises(RuntimeError):
+            ift.MetricGaussianKL.make(**args)
+        return
+    kl = ift.MetricGaussianKL.make(**args)
     assert_(len(ic.history) > 0)
     assert_(len(ic.history) == len(ic.history.time_stamps))
     assert_(len(ic.history) == len(ic.history.energy_values))
@@ -60,12 +64,11 @@ def test_kl(constants, point_estimates, mirror_samples, mf):
     assert_(len(ic.history) == len(ic.history.energy_values))
 
     locsamp = kl._local_samples
-    klpure = ift.MetricGaussianKL(mean0,
-                                  h,
-                                  nsamps,
-                                  mirror_samples=mirror_samples,
-                                  napprox=0,
-                                  _local_samples=locsamp)
+    if isinstance(mean0, ift.MultiField):
+        _, tmph = h.simplify_for_constant_input(mean0.extract_by_keys(constants))
+    else:
+        tmph = h
+    klpure = ift.MetricGaussianKL(mean0, tmph, nsamps, mirror_samples, None, locsamp, False, True)
 
     # Test number of samples
     expected_nsamps = 2*nsamps if mirror_samples else nsamps
diff --git a/test/test_linearization.py b/test/test_linearization.py
index 151557f5cc2750b69c39e4e2273d86a547cc9e2b..5a3cf0fe7d69f3372739600f36625fd8a8a4e0a7 100644
--- a/test/test_linearization.py
+++ b/test/test_linearization.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_, assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 from .common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
diff --git a/test/test_minimizers.py b/test/test_minimizers.py
index 3c847fd5add15881ca03e8fd3f12b70010133d55..a5da0a9aea9251476af7cf3242b999014e1fd537 100644
--- a/test/test_minimizers.py
+++ b/test/test_minimizers.py
@@ -21,7 +21,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_allclose, assert_equal
 
-import nifty6 as ift
+import nifty7 as ift
 from .common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
diff --git a/test/test_mpi/test_kl.py b/test/test_mpi/test_kl.py
index 7e3b4e934dd51c89ac7221085134fe19e08eacaa..f0eef5655cd42106c05c5ca79da554da35a2c58f 100644
--- a/test/test_mpi/test_kl.py
+++ b/test/test_mpi/test_kl.py
@@ -18,9 +18,9 @@
 import numpy as np
 import pytest
 from mpi4py import MPI
-from numpy.testing import assert_, assert_equal
+from numpy.testing import assert_, assert_equal, assert_raises
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import setup_function, teardown_function
 
@@ -58,17 +58,29 @@ def test_kl(constants, point_estimates, mirror_samples, mode, mf):
             'n_samples': 2,
             'mean': mean0,
             'hamiltonian': h}
+    if isinstance(mean0, ift.MultiField) and set(point_estimates) == set(mean0.keys()):
+        with assert_raises(RuntimeError):
+            ift.MetricGaussianKL.make(**args, comm=comm)
+        return
     if mode == 0:
-        kl0 = ift.MetricGaussianKL(**args, comm=comm)
+        kl0 = ift.MetricGaussianKL.make(**args, comm=comm)
         locsamp = kl0._local_samples
-        kl1 = ift.MetricGaussianKL(**args, comm=comm, _local_samples=locsamp)
+        if isinstance(mean0, ift.MultiField):
+            _, tmph = h.simplify_for_constant_input(mean0.extract_by_keys(constants))
+        else:
+            tmph = h
+        kl1 = ift.MetricGaussianKL(mean0, tmph, 2, mirror_samples, comm, locsamp, False, True)
     elif mode == 1:
-        kl0 = ift.MetricGaussianKL(**args)
+        kl0 = ift.MetricGaussianKL.make(**args)
         samples = kl0._local_samples
         ii = len(samples)//2
         slc = slice(None, ii) if rank == 0 else slice(ii, None)
         locsamp = samples[slc]
-        kl1 = ift.MetricGaussianKL(**args, comm=comm, _local_samples=locsamp)
+        if isinstance(mean0, ift.MultiField):
+            _, tmph = h.simplify_for_constant_input(mean0.extract_by_keys(constants))
+        else:
+            tmph = h
+        kl1 = ift.MetricGaussianKL(mean0, tmph, 2, mirror_samples, comm, locsamp, False, True)
 
     # Test number of samples
     expected_nsamps = 2*nsamps if mirror_samples else nsamps
diff --git a/test/test_multi_field.py b/test/test_multi_field.py
index 16f26e2e7ec61baad5f9a36ebdfd9157960fb3c7..5923d455efee4a9a8605dcc3affc752a7a6986b1 100644
--- a/test/test_multi_field.py
+++ b/test/test_multi_field.py
@@ -18,7 +18,7 @@
 import numpy as np
 from numpy.testing import assert_allclose, assert_equal
 
-import nifty6 as ift
+import nifty7 as ift
 from .common import setup_function, teardown_function
 
 dom = ift.makeDomain({"d1": ift.RGSpace(10)})
@@ -61,7 +61,7 @@ def test_blockdiagonal():
     op = ift.BlockDiagonalOperator(
         dom, {"d1": ift.ScalingOperator(dom["d1"], 20.)})
     op2 = op(op)
-    ift.extra.consistency_check(op2)
+    ift.extra.check_linear_operator(op2)
     assert_equal(type(op2), ift.BlockDiagonalOperator)
     f1 = op2(ift.full(dom, 1))
     for val in f1.values():
diff --git a/test/test_operator_tree_optimiser.py b/test/test_operator_tree_optimiser.py
index dab6cd61fd8fd9961cba8b1c47999efd70f5c153..7ad35e094d148b06f46321a843cd43b344e2f201 100644
--- a/test/test_operator_tree_optimiser.py
+++ b/test/test_operator_tree_optimiser.py
@@ -19,7 +19,7 @@
 from numpy.testing import assert_, assert_allclose
 import numpy as np
 from copy import deepcopy
-import nifty6 as ift
+import nifty7 as ift
 
 
 class CountingOp(ift.Operator):
diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py
index d47255a69f372ecdf9753368210913f6b3769569..66422a84caff64fb793391b5cdc892c86adb4082 100644
--- a/test/test_operators/test_adjoint.py
+++ b/test/test_operators/test_adjoint.py
@@ -18,7 +18,7 @@
 import numpy as np
 import pytest
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
@@ -42,7 +42,7 @@ def testLOSResponse(sp, dtype):
     sigma_low = 1e-4*ift.random.current_rng().standard_normal(10)
     sigma_ups = 1e-5*ift.random.current_rng().standard_normal(10)
     op = ift.LOSResponse(sp, starts, ends, sigma_low, sigma_ups)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
@@ -50,13 +50,13 @@ def testOperatorCombinations(sp, dtype):
     a = ift.DiagonalOperator(ift.Field.from_random(sp, "normal", dtype=dtype))
     b = ift.DiagonalOperator(ift.Field.from_random(sp, "normal", dtype=dtype))
     op = ift.SandwichOperator.make(a, b)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
     op = a(b)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
     op = a + b
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
     op = a - b
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 def testLinearInterpolator():
@@ -65,60 +65,60 @@ def testLinearInterpolator():
     pos[0, :] *= 0.9
     pos[1, :] *= 7*3.5
     op = ift.LinearInterpolator(sp, pos)
-    ift.extra.consistency_check(op)
+    ift.extra.check_linear_operator(op)
 
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 def testRealizer(sp):
     op = ift.Realizer(sp)
-    ift.extra.consistency_check(op, np.complex128, np.float64,
+    ift.extra.check_linear_operator(op, np.complex128, np.float64,
                                 only_r_linear=True)
 
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 def testImaginizer(sp):
     op = ift.Imaginizer(sp)
-    ift.extra.consistency_check(op, np.complex128, np.float64,
+    ift.extra.check_linear_operator(op, np.complex128, np.float64,
                                 only_r_linear=True)
     loc = ift.from_random(op.domain, dtype=np.complex128)
-    ift.extra.check_jacobian_consistency(op, loc)
+    ift.extra.check_operator(op, loc)
 
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 def testConjugationOperator(sp):
     op = ift.ConjugationOperator(sp)
-    ift.extra.consistency_check(op, np.complex128, np.complex128,
+    ift.extra.check_linear_operator(op, np.complex128, np.complex128,
                                 only_r_linear=True)
 
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 def testOperatorAdaptor(sp, dtype):
     op = ift.DiagonalOperator(ift.Field.from_random(sp, "normal", 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)
+    ift.extra.check_linear_operator(op.adjoint, dtype, dtype)
+    ift.extra.check_linear_operator(op.inverse, dtype, dtype)
+    ift.extra.check_linear_operator(op.inverse.adjoint, dtype, dtype)
+    ift.extra.check_linear_operator(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)
+    ift.extra.check_linear_operator(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)
+    ift.extra.check_linear_operator(op, dtype, dtype)
     op = ift.NullOperator(sp1, mdom2)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
     op = ift.NullOperator(mdom1, sp2)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(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)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
@@ -129,43 +129,43 @@ def testDOFDistributor(sp, dtype):
     dofdex = np.arange(sp.size).reshape(sp.shape) % 3
     dofdex = ift.Field.from_raw(sp, dofdex)
     op = ift.DOFDistributor(dofdex)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 @pmp('sp', _h_spaces)
 def testPPO(sp, dtype):
     op = ift.PowerDistributor(target=sp)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(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)
+    ift.extra.check_linear_operator(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)
+    ift.extra.check_linear_operator(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)
+    ift.extra.check_linear_operator(op, dtype, dtype)
     op = ift.FFTOperator(sp.get_default_codomain())
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(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)
+    ift.extra.check_linear_operator(op, dtype, dtype)
     op = ift.HartleyOperator(sp.get_default_codomain())
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 @pmp('sp', _h_spaces)
 def testHarmonic(sp, dtype):
     op = ift.HarmonicTransformOperator(sp)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 @pmp('sp', _p_spaces)
@@ -175,19 +175,19 @@ def testMask(sp, dtype):
     mask[f > 0] = 1
     mask = ift.Field.from_raw(sp, mask)
     op = ift.MaskOperator(mask)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 @pmp('sp', _h_spaces + _p_spaces)
 def testDiagonal(sp, dtype):
     op = ift.DiagonalOperator(ift.Field.from_random(sp, dtype=dtype))
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(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)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 @pmp('spaces', [0, 1, 2, 3, (0, 1), (0, 2), (0, 1, 2), (0, 2, 3), (1, 3)])
@@ -195,7 +195,7 @@ def testGeometryRemover(sp, dtype):
 def testContractionOperator(spaces, wgt, dtype):
     dom = (ift.RGSpace(1), ift.RGSpace(2), ift.GLSpace(3), ift.HPSpace(2))
     op = ift.ContractionOperator(dom, spaces, wgt)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 def testDomainTupleFieldInserter():
@@ -203,7 +203,7 @@ def testDomainTupleFieldInserter():
                                    ift.UnstructuredDomain(7),
                                    ift.RGSpace([4, 22])))
     op = ift.DomainTupleFieldInserter(target, 1, (5,))
-    ift.extra.consistency_check(op)
+    ift.extra.check_linear_operator(op)
 
 
 @pmp('space', [0, 2])
@@ -214,7 +214,7 @@ def testZeroPadder(space, factor, dtype, central):
            ift.HPSpace(2))
     newshape = [int(factor*ll) for ll in dom[space].shape]
     op = ift.FieldZeroPadder(dom, newshape, space, central)
-    ift.extra.consistency_check(op, dtype, dtype)
+    ift.extra.check_linear_operator(op, dtype, dtype)
 
 
 @pmp('args', [[ift.RGSpace((13, 52, 40)), (4, 6, 25), None],
@@ -223,7 +223,7 @@ def testZeroPadder(space, factor, dtype, central):
               [(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)
+    ift.extra.check_linear_operator(op)
 
 
 @pmp('fdomain', [(ift.RGSpace((2, 3, 2)), ift.RGSpace((4,), distances=(7.,))),
@@ -233,7 +233,7 @@ def testRegridding(args):
 def testOuter(fdomain, domain):
     f = ift.from_random(ift.makeDomain(fdomain), 'normal')
     op = ift.OuterProduct(domain, f)
-    ift.extra.consistency_check(op)
+    ift.extra.check_linear_operator(op)
 
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
@@ -247,25 +247,25 @@ def testValueInserter(sp, seed):
             else:
                 ind.append(int(ift.random.current_rng().integers(0, ss-1)))
         op = ift.ValueInserter(sp, ind)
-        ift.extra.consistency_check(op)
+        ift.extra.check_linear_operator(op)
 
 
 @pmp('sp', _pow_spaces)
 def testSlopeRemover(sp):
     op = ift.library.correlated_fields._SlopeRemover(sp)
-    ift.extra.consistency_check(op)
+    ift.extra.check_linear_operator(op)
 
 
 @pmp('sp', _pow_spaces)
 def testTwoLogIntegrations(sp):
     op = ift.library.correlated_fields._TwoLogIntegrations(sp)
-    ift.extra.consistency_check(op)
+    ift.extra.check_linear_operator(op)
 
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 def testSpecialSum(sp):
     op = ift.library.correlated_fields._SpecialSum(sp)
-    ift.extra.consistency_check(op)
+    ift.extra.check_linear_operator(op)
 
 
 @pmp('dofdex', [(0,), (1,), (0, 1), (1, 0)])
@@ -273,17 +273,17 @@ def testCorFldDistr(dofdex):
     tgt = ift.UnstructuredDomain(len(dofdex))
     dom = ift.UnstructuredDomain(2)
     op = ift.library.correlated_fields._Distributor(dofdex, dom, tgt)
-    ift.extra.consistency_check(op)
+    ift.extra.check_linear_operator(op)
 
 
 def metatestMatrixProductOperator(sp, mat_shape, seed, **kwargs):
     with ift.random.Context(seed):
         mat = ift.random.current_rng().standard_normal(mat_shape)
         op = ift.MatrixProductOperator(sp, mat, **kwargs)
-        ift.extra.consistency_check(op)
+        ift.extra.check_linear_operator(op)
         mat = mat + 1j*ift.random.current_rng().standard_normal(mat_shape)
         op = ift.MatrixProductOperator(sp, mat, **kwargs)
-        ift.extra.consistency_check(op)
+        ift.extra.check_linear_operator(op)
 
 
 @pmp('sp', [ift.RGSpace(10)])
@@ -323,11 +323,11 @@ def testPartialExtractor(seed):
         dom = ift.MultiDomain.make(dom)
         tgt = ift.MultiDomain.make(tgt)
         op = ift.PartialExtractor(dom, tgt)
-        ift.extra.consistency_check(op)
+        ift.extra.check_linear_operator(op)
 
 
 @pmp('seed', [12, 3])
 def testSlowFieldAdapter(seed):
     dom = {'a': ift.RGSpace(1), 'b': ift.RGSpace(2)}
     op = ift.operators.simple_linear_operators._SlowFieldAdapter(dom, 'a')
-    ift.extra.consistency_check(op)
+    ift.extra.check_linear_operator(op)
diff --git a/test/test_operators/test_composed_operator.py b/test/test_operators/test_composed_operator.py
index 55b517a7648160875591ff9e1ad5aab7cf03442c..874b23867f3a2d797d660df5cff9bd8fbb049db1 100644
--- a/test/test_operators/test_composed_operator.py
+++ b/test/test_operators/test_composed_operator.py
@@ -17,7 +17,7 @@
 
 from numpy.testing import assert_allclose, assert_equal
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
diff --git a/test/test_operators/test_conjugation_operator.py b/test/test_operators/test_conjugation_operator.py
index 6955c447e7bf040b26586ff26f3ae34f999b40dc..a0aea1dd28351bbc02aaff90f59b674e0fba5861 100644
--- a/test/test_operators/test_conjugation_operator.py
+++ b/test/test_operators/test_conjugation_operator.py
@@ -18,7 +18,7 @@
 import numpy as np
 from numpy.testing import assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import setup_function, teardown_function
 
@@ -34,7 +34,7 @@ def test_conjugation_operator():
     res3 = arr.conjugate()
     assert_allclose(res1.val, res2.val)
     assert_allclose(res1.val, res3)
-    ift.extra.consistency_check(op, domain_dtype=np.float64,
-                                target_dtype=np.float64)
-    ift.extra.consistency_check(op, domain_dtype=np.complex128,
-                                target_dtype=np.complex128, only_r_linear=True)
+    ift.extra.check_linear_operator(op, domain_dtype=np.float64,
+                                    target_dtype=np.float64)
+    ift.extra.check_linear_operator(op, domain_dtype=np.complex128,
+                                    target_dtype=np.complex128, only_r_linear=True)
diff --git a/test/test_operators/test_contraction_operator.py b/test/test_operators/test_contraction_operator.py
index 9a01be466dfef7391d0cc241348d41b82b451f1c..b1fab0d3896a5db85d260f62ad73c1cf5933a115 100644
--- a/test/test_operators/test_contraction_operator.py
+++ b/test/test_operators/test_contraction_operator.py
@@ -18,7 +18,7 @@
 import numpy as np
 from numpy.testing import assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import setup_function, teardown_function
 
@@ -51,7 +51,7 @@ def test_operator_sum():
     assert_allclose(res7.val, res8.val)
     assert_allclose(res7.val, res9)
     for op in [op1, op2, op3]:
-        ift.extra.consistency_check(op, domain_dtype=np.float64,
-                                    target_dtype=np.float64)
-        ift.extra.consistency_check(op, domain_dtype=np.complex128,
-                                    target_dtype=np.complex128)
+        ift.extra.check_linear_operator(op, domain_dtype=np.float64,
+                                        target_dtype=np.float64)
+        ift.extra.check_linear_operator(op, domain_dtype=np.complex128,
+                                        target_dtype=np.complex128)
diff --git a/test/test_operators/test_convolution_operators.py b/test/test_operators/test_convolution_operators.py
index 6f4f6495c044862241fe6883b6fba76ca5be3e8b..8e66f6bc0a32bb338cfce1759a0a4566b1a9accd 100644
--- a/test/test_operators/test_convolution_operators.py
+++ b/test/test_operators/test_convolution_operators.py
@@ -17,7 +17,7 @@
 
 from numpy.testing import assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 import numpy as np
 
 from ..common import list2fixture, setup_function, teardown_function
diff --git a/test/test_operators/test_correlated_fields.py b/test/test_operators/test_correlated_fields.py
index 20c7d12362036cb1e4a47a5e055c942650d48eaf..8d1ef09d6f1ca6483c8c93183c0d1811b189c463 100644
--- a/test/test_operators/test_correlated_fields.py
+++ b/test/test_operators/test_correlated_fields.py
@@ -18,7 +18,7 @@
 import pytest
 from numpy.testing import assert_, assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import setup_function, teardown_function
 
@@ -43,7 +43,7 @@ def testDistributor(dofdex, seed):
         target = ift.makeDomain((ift.UnstructuredDomain(N_copies), dom))
         op = ift.library.correlated_fields._Distributor(
             dofdex, target, distributed_target)
-        ift.extra.consistency_check(op)
+        ift.extra.check_linear_operator(op)
 
 
 @pmp('sspace', [
@@ -108,9 +108,5 @@ def testAmplitudesInvariants(sspace, N):
     assert_(op.target[-1] == fsspace)
 
     for ampl in fa.normalized_amplitudes:
-        ift.extra.check_jacobian_consistency(ampl,
-                                             ift.from_random(ampl.domain),
-                                             ntries=10)
-    ift.extra.check_jacobian_consistency(op,
-                                         ift.from_random(op.domain),
-                                         ntries=10)
+        ift.extra.check_operator(ampl, ift.from_random(ampl.domain), ntries=10)
+    ift.extra.check_operator(op, ift.from_random(op.domain), ntries=10)
diff --git a/test/test_operators/test_diagonal_operator.py b/test/test_operators/test_diagonal_operator.py
index f165a22731e6804a12d928f9ad842554eddd4eb4..0b79975ca681d78849a6e55d50ed11da256a02a0 100644
--- a/test/test_operators/test_diagonal_operator.py
+++ b/test/test_operators/test_diagonal_operator.py
@@ -17,7 +17,7 @@
 
 from numpy.testing import assert_allclose, assert_equal
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
diff --git a/test/test_operators/test_einsum.py b/test/test_operators/test_einsum.py
index 41df830aab224fa9d1de026053225052ca9a385f..c426d5226c4ea30e2c6eefcbed03047015ee7992 100644
--- a/test/test_operators/test_einsum.py
+++ b/test/test_operators/test_einsum.py
@@ -19,8 +19,8 @@
 import numpy as np
 from numpy.testing import assert_, assert_allclose
 
-import nifty6 as ift
-from nifty6.extra import check_jacobian_consistency, consistency_check
+import nifty7 as ift
+from nifty7.extra import check_operator, check_linear_operator
 
 from ..common import list2fixture, setup_function, teardown_function
 
@@ -42,7 +42,7 @@ def test_linear_einsum_outer(space1, space2, dtype, n_invocations=10):
     ss = "i,ij,j->ij"
     key_order = ("dom01", "dom02")
     le = ift.LinearEinsum(space2, mf, ss, key_order=key_order)
-    assert_(consistency_check(le, domain_dtype=dtype, target_dtype=dtype) is None)
+    assert_(check_linear_operator(le, domain_dtype=dtype, target_dtype=dtype) is None)
 
     le_ift = ift.DiagonalOperator(mf["dom01"], domain=mf_dom["dom02"], spaces=0) @ ift.DiagonalOperator(mf["dom02"])
     le_ift = le_ift @ ift.OuterProduct(ift.DomainTuple.make(mf_dom["dom02"][1]),
@@ -63,7 +63,7 @@ def test_linear_einsum_contraction(space1, space2, dtype, n_invocations=10):
     ss = "i,ij,j->i"
     key_order = ("dom01", "dom02")
     le = ift.LinearEinsum(space2, mf, ss, key_order=key_order)
-    assert_(consistency_check(le, domain_dtype=dtype, target_dtype=dtype) is None)
+    assert_(check_linear_operator(le, domain_dtype=dtype, target_dtype=dtype) is None)
 
     le_ift = ift.ContractionOperator(mf_dom["dom02"], 1)
     le_ift = le_ift @ ift.DiagonalOperator(mf["dom01"], domain=mf_dom["dom02"], spaces=0)
@@ -116,7 +116,7 @@ def test_linear_einsum_transpose(space1, space2, dtype, n_invocations=10):
     mf = ift.MultiField.from_dict({})
     ss = "ij->ji"
     le = ift.LinearEinsum(dom, mf, ss)
-    assert_(consistency_check(le, domain_dtype=dtype, target_dtype=dtype) is None)
+    assert_(check_linear_operator(le, domain_dtype=dtype, target_dtype=dtype) is None)
 
     # SwitchSpacesOperator is equivalent to LinearEinsum with "ij->ji"
     le_ift = _SwitchSpacesOperator(dom, 1)
@@ -138,7 +138,7 @@ def test_multi_linear_einsum_outer(space1, space2, dtype):
     ss = "i,ij,j->ij"
     key_order = ("dom01", "dom02", "dom03")
     mle = ift.MultiLinearEinsum(mf_dom, ss, key_order=key_order)
-    check_jacobian_consistency(mle, ift.from_random(mle.domain, "normal", dtype=dtype), ntries=ntries)
+    check_operator(mle, ift.from_random(mle.domain, "normal", dtype=dtype), ntries=ntries)
 
     outer_i = ift.OuterProduct(
         ift.DomainTuple.make(mf_dom["dom02"][0]), ift.full(mf_dom["dom03"], 1.)
diff --git a/test/test_operators/test_fft_operator.py b/test/test_operators/test_fft_operator.py
index 54c0a8de4fe7711c923d1411fa01489a08eaa161..bdfccc43bf3fbbf8f2c4fd9a08a203f12cad0a23 100644
--- a/test/test_operators/test_fft_operator.py
+++ b/test/test_operators/test_fft_operator.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_, assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
diff --git a/test/test_operators/test_harmonic_transform_operator.py b/test/test_operators/test_harmonic_transform_operator.py
index 71a597d9a99680c3ee238ed62110460c74a5eca4..1cc13aca1be0ec6b65f91af6c05135578c4ceef5 100644
--- a/test/test_operators/test_harmonic_transform_operator.py
+++ b/test/test_operators/test_harmonic_transform_operator.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
diff --git a/test/test_operators/test_integration.py b/test/test_operators/test_integration.py
index 6749bebf38eafd28a8a79114678af35871a030f2..5d8dfc148d8b7cb95909bea959295ea9f9993d35 100644
--- a/test/test_operators/test_integration.py
+++ b/test/test_operators/test_integration.py
@@ -18,7 +18,7 @@
 import numpy as np
 from numpy.testing import assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import setup_function, teardown_function
 
@@ -49,7 +49,7 @@ def test_integration_operator():
     assert_allclose(res7.val, res8.val)
     assert_allclose(res7.val, res9)
     for op in [op1, op2, op3]:
-        ift.extra.consistency_check(op, domain_dtype=np.float64,
-                                    target_dtype=np.float64)
-        ift.extra.consistency_check(op, domain_dtype=np.complex128,
-                                    target_dtype=np.complex128)
+        ift.extra.check_linear_operator(op, domain_dtype=np.float64,
+                                        target_dtype=np.float64)
+        ift.extra.check_linear_operator(op, domain_dtype=np.complex128,
+                                        target_dtype=np.complex128)
diff --git a/test/test_operators/test_interpolated.py b/test/test_operators/test_interpolated.py
index b7d3d47e6c9357075ac160cac4ef322aa300dca9..f15c26b7b92ef2a4ffe923119c7283c782d47130 100644
--- a/test/test_operators/test_interpolated.py
+++ b/test/test_operators/test_interpolated.py
@@ -20,7 +20,7 @@ import pytest
 from numpy.testing import assert_allclose
 from scipy.stats import invgamma, norm
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
diff --git a/test/test_operators/test_jacobian.py b/test/test_operators/test_jacobian.py
index c5715029b6b2762f19f3b0a49ae89b369af57e4b..f262d78b6c8033820132ad99051d9dec206dcd1e 100644
--- a/test/test_operators/test_jacobian.py
+++ b/test/test_operators/test_jacobian.py
@@ -17,7 +17,7 @@
 
 import pytest
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
@@ -42,7 +42,7 @@ def testBasics(space, seed):
         s = ift.from_random(space, 'normal')
         var = ift.Linearization.make_var(s)
         model = ift.ScalingOperator(var.target, 6.)
-        ift.extra.check_jacobian_consistency(model, var.val, ntries=ntries)
+        ift.extra.check_operator(model, var.val, ntries=ntries)
 
 
 @pmp('type1', ['Variable', 'Constant'])
@@ -56,36 +56,36 @@ def testBinary(type1, type2, space, seed):
         select_s2 = ift.ducktape(None, dom2, "s2")
         model = select_s1*select_s2
         pos = ift.from_random(dom, "normal")
-        ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
+        ift.extra.check_operator(model, pos, ntries=ntries)
         model = select_s1 + select_s2
         pos = ift.from_random(dom, "normal")
-        ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
+        ift.extra.check_operator(model, pos, ntries=ntries)
         model = select_s1.scale(3.)
         pos = ift.from_random(dom1, "normal")
-        ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
+        ift.extra.check_operator(model, pos, ntries=ntries)
         model = ift.ScalingOperator(space, 2.456)(select_s1*select_s2)
         pos = ift.from_random(dom, "normal")
-        ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
+        ift.extra.check_operator(model, pos, ntries=ntries)
         model = (2.456*(select_s1*select_s2)).ptw("sigmoid")
         pos = ift.from_random(dom, "normal")
-        ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
+        ift.extra.check_operator(model, pos, ntries=ntries)
         pos = ift.from_random(dom, "normal")
         model = ift.OuterProduct(ift.makeDomain(space), pos['s1'])
-        ift.extra.check_jacobian_consistency(model, pos['s2'], ntries=ntries)
+        ift.extra.check_operator(model, pos['s2'], ntries=ntries)
         model = select_s1**2
         pos = ift.from_random(dom1, "normal")
-        ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
+        ift.extra.check_operator(model, pos, ntries=ntries)
         model = select_s1.clip(-1, 1)
         pos = ift.from_random(dom1, "normal")
-        ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
+        ift.extra.check_operator(model, pos, ntries=ntries)
         f = ift.from_random(space, "normal")
         model = select_s1.clip(f-0.1, f+1.)
         pos = ift.from_random(dom1, "normal")
-        ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
+        ift.extra.check_operator(model, pos, ntries=ntries)
         if isinstance(space, ift.RGSpace):
             model = ift.FFTOperator(space)(select_s1*select_s2)
             pos = ift.from_random(dom, "normal")
-            ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
+            ift.extra.check_operator(model, pos, ntries=ntries)
 
 
 def testSpecialDistributionOps(space, seed):
@@ -94,9 +94,9 @@ def testSpecialDistributionOps(space, seed):
         alpha = 1.5
         q = 0.73
         model = ift.InverseGammaOperator(space, alpha, q)
-        ift.extra.check_jacobian_consistency(model, pos, ntries=20)
+        ift.extra.check_operator(model, pos, ntries=20)
         model = ift.UniformOperator(space, alpha, q)
-        ift.extra.check_jacobian_consistency(model, pos, ntries=20)
+        ift.extra.check_operator(model, pos, ntries=20)
 
 
 @pmp('neg', [True, False])
@@ -105,9 +105,9 @@ def testAdder(space, seed, neg):
         f = ift.from_random(space, 'normal')
         f1 = ift.from_random(space, 'normal')
         op = ift.Adder(f1, neg)
-        ift.extra.check_jacobian_consistency(op, f, ntries=ntries)
+        ift.extra.check_operator(op, f, ntries=ntries)
         op = ift.Adder(f1.val.ravel()[0], neg=neg, domain=space)
-        ift.extra.check_jacobian_consistency(op, f, ntries=ntries)
+        ift.extra.check_operator(op, f, ntries=ntries)
 
 
 @pmp('target', [ift.RGSpace(64, distances=.789, harmonic=True),
@@ -125,7 +125,7 @@ def testDynamicModel(target, causal, minimum_phase, seed):
                'minimum_phase': minimum_phase}
         model, _ = ift.dynamic_operator(**dct)
         pos = ift.from_random(model.domain, 'normal')
-        ift.extra.check_jacobian_consistency(model, pos, tol=1e-7, ntries=ntries)
+        ift.extra.check_operator(model, pos, ntries=ntries)
         if len(target.shape) > 1:
             dct = {
                 'target': target,
@@ -144,7 +144,7 @@ def testDynamicModel(target, causal, minimum_phase, seed):
             dct['quant'] = 5
             model, _ = ift.dynamic_lightcone_operator(**dct)
             pos = ift.from_random(model.domain, 'normal')
-            ift.extra.check_jacobian_consistency(model, pos, tol=1e-7, ntries=ntries)
+            ift.extra.check_operator(model, pos, ntries=ntries)
 
 
 @pmp('h_space', _h_spaces)
@@ -160,11 +160,11 @@ def testNormalization(h_space, specialbinbounds, logarithmic, nbin):
     dom = ift.PowerSpace(h_space, binbounds)
     op = ift.library.correlated_fields._Normalization(dom)
     pos = 0.1 * ift.from_random(op.domain, 'normal')
-    ift.extra.check_jacobian_consistency(op, pos, ntries=10)
+    ift.extra.check_operator(op, pos, ntries=10)
 
 
 @pmp('N', [1, 20])
 def testLognormalTransform(N):
     op = ift.LognormalTransform(1, 0.2, '', N)
     loc = ift.from_random(op.domain)
-    ift.extra.check_jacobian_consistency(op, loc, ntries=10)
+    ift.extra.check_operator(op, loc, ntries=10)
diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py
index 4ea0d26fd0435b40b656c763ad4396819a749c38..d636ba67f8204f3bbd17819eee480e3b3b31d6ed 100644
--- a/test/test_operators/test_nft.py
+++ b/test/test_operators/test_nft.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_
 
-import nifty6 as ift
+import nifty7 as ift
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
@@ -104,7 +104,7 @@ def test_build(nu, nv, N, eps):
     # Consistency checks
     flt = np.float64
     cmplx = np.complex128
-    ift.extra.consistency_check(R0, cmplx, flt, only_r_linear=True)
-    ift.extra.consistency_check(R1, flt, flt)
-    ift.extra.consistency_check(R, cmplx, flt, only_r_linear=True)
-    ift.extra.consistency_check(RF, cmplx, flt, only_r_linear=True)
+    ift.extra.check_linear_operator(R0, cmplx, flt, only_r_linear=True)
+    ift.extra.check_linear_operator(R1, flt, flt)
+    ift.extra.check_linear_operator(R, cmplx, flt, only_r_linear=True)
+    ift.extra.check_linear_operator(RF, cmplx, flt, only_r_linear=True)
diff --git a/test/test_operators/test_normal_operators.py b/test/test_operators/test_normal_operators.py
index 6895a8e928813a17cbfc81768aee60145939fd00..227700546dcacf6cef2220848e6bb12c5acbae35 100644
--- a/test/test_operators/test_normal_operators.py
+++ b/test/test_operators/test_normal_operators.py
@@ -18,12 +18,13 @@
 import pytest
 from numpy.testing import assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
 
+
 @pmp('mean', [3., -2])
 @pmp('std', [1.5, 0.1])
 @pmp('seed', [7, 21])
@@ -39,7 +40,8 @@ def test_normal_transform(mean, std, seed):
         assert_allclose(res.val.std(), std, rtol=0.1)
 
         loc = ift.from_random(op.domain)
-        ift.extra.check_jacobian_consistency(op, loc)
+        ift.extra.check_operator(op, loc)
+
 
 @pmp('mean', [0.01, 10.])
 @pmp('std_fct', [0.01, 0.1])
@@ -54,4 +56,4 @@ def test_lognormal_transform(mean, std_fct, seed):
         assert_allclose(res.val.std(), std, rtol=0.1)
 
         loc = ift.from_random(op.domain)
-        ift.extra.check_jacobian_consistency(op, loc)
+        ift.extra.check_operator(op, loc)
diff --git a/test/test_operators/test_partial_multifield_insert.py b/test/test_operators/test_partial_multifield_insert.py
index dd07b085fa5edd73d58ebb856c9c6d48375ef1bd..e38dbbc7ccd980d95f8d0841172a9160fb5bc9c7 100644
--- a/test/test_operators/test_partial_multifield_insert.py
+++ b/test/test_operators/test_partial_multifield_insert.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_, assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
@@ -39,7 +39,7 @@ def test_part_mf_insert():
     b = op4 + op5
     op = a.partial_insert(b)
     fld = ift.from_random(op.domain, 'normal')
-    ift.extra.check_jacobian_consistency(op, fld, ntries=ntries)
+    ift.extra.check_operator(op, fld, ntries=ntries)
     assert_(op.domain is ift.MultiDomain.union(
         [op1.domain, op2.domain, op4.domain, op5.domain]))
     assert_(op.target is ift.MultiDomain.union(
diff --git a/test/test_operators/test_regridding.py b/test/test_operators/test_regridding.py
index 1ad8903b6e106254487d82803dad2c25fb7ad815..def6637b9d8f8335ad95968c76e8a9e50eb2e6be 100644
--- a/test/test_operators/test_regridding.py
+++ b/test/test_operators/test_regridding.py
@@ -17,7 +17,7 @@
 
 from numpy.testing import assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
diff --git a/test/test_operators/test_representation.py b/test/test_operators/test_representation.py
index 45c402126be077a88bbc2656a5d79a9138fd94ba..aa139085b8e161028e914e2f077b90ac8b425b68 100644
--- a/test/test_operators/test_representation.py
+++ b/test/test_operators/test_representation.py
@@ -18,7 +18,7 @@
 import numpy as np
 import pytest
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
@@ -83,7 +83,7 @@ def testOperatorAdaptor(sp, dtype):
 @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)
+    ift.extra.check_linear_operator(op, dtype, dtype)
     mdom1 = ift.MultiDomain.make({'a': sp1})
     mdom2 = ift.MultiDomain.make({'b': sp2})
     _check_repr(ift.NullOperator(mdom1, mdom2))
diff --git a/test/test_operators/test_selection_operators.py b/test/test_operators/test_selection_operators.py
index 0e6fa555e3f730033beb74e27d2e1f344a0f2af1..5554a574eee4bd6f7fdc818887ed077a65ba78aa 100644
--- a/test/test_operators/test_selection_operators.py
+++ b/test/test_operators/test_selection_operators.py
@@ -16,12 +16,13 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
+import numpy as np
 import pytest
 from numpy.testing import assert_allclose, assert_array_equal
-from nifty6.extra import consistency_check
 
-import numpy as np
-import nifty6 as ift
+import nifty7 as ift
+from nifty7.extra import check_linear_operator
+
 from ..common import list2fixture, setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
@@ -50,7 +51,7 @@ def test_split_operator_first_axes_without_intersections(
         dom, {f"{i:06d}": (si, )
               for i, si in enumerate(split_idx)}
     )
-    assert consistency_check(split) is None
+    assert check_linear_operator(split) is None
 
     r = ift.from_random(dom, "normal")
     split_r = split(r)
@@ -79,7 +80,7 @@ def test_split_operator_first_axes_with_intersections(
               for i, si in enumerate(split_idx)}
     )
     print(split_idx)
-    assert consistency_check(split) is None
+    assert check_linear_operator(split) is None
 
     r = ift.from_random(dom, "normal")
     split_r = split(r)
diff --git a/test/test_operators/test_simplification.py b/test/test_operators/test_simplification.py
index bd3d8c4d6139a96ee06360d6083d33bd168fb40f..e8122e00f74c4c89bba8a7abd3e7608893e24d40 100644
--- a/test/test_operators/test_simplification.py
+++ b/test/test_operators/test_simplification.py
@@ -17,23 +17,23 @@
 
 from numpy.testing import assert_allclose, assert_equal
 
-import nifty6 as ift
+import nifty7 as ift
 from ..common import setup_function, teardown_function
 
 
 def test_simplification():
-    from nifty6.operators.operator import _ConstantOperator
+    from nifty7.operators.simplify_for_const import ConstantOperator
     f1 = ift.Field.full(ift.RGSpace(10), 2.)
     op = ift.FFTOperator(f1.domain)
     _, op2 = op.simplify_for_constant_input(f1)
-    assert_equal(isinstance(op2, _ConstantOperator), True)
+    assert_equal(isinstance(op2, ConstantOperator), True)
     assert_allclose(op(f1).val, op2(f1).val)
 
     dom = {"a": ift.RGSpace(10)}
     f1 = ift.full(dom, 2.)
     op = ift.FFTOperator(f1.domain["a"]).ducktape("a")
     _, op2 = op.simplify_for_constant_input(f1)
-    assert_equal(isinstance(op2, _ConstantOperator), True)
+    assert_equal(isinstance(op2, ConstantOperator), True)
     assert_allclose(op(f1).val, op2(f1).val)
 
     dom = {"a": ift.RGSpace(10), "b": ift.RGSpace(5)}
@@ -45,7 +45,7 @@ def test_simplification():
     op = (o1.ducktape("a").ducktape_left("a") +
           o2.ducktape("b").ducktape_left("b"))
     _, op2 = op.simplify_for_constant_input(f2)
-    assert_equal(isinstance(op2._op1, _ConstantOperator), True)
+    assert_equal(isinstance(op2._op1, ConstantOperator), True)
     assert_allclose(op(f1)["a"].val, op2(f1)["a"].val)
     assert_allclose(op(f1)["b"].val, op2(f1)["b"].val)
     lin = ift.Linearization.make_var(ift.MultiField.full(op2.domain, 2.), True)
diff --git a/test/test_operators/test_smoothing_operator.py b/test/test_operators/test_smoothing_operator.py
index 2ce2e6df30bdecc5270562d59d22aa3152e88ca0..b53dac7ed94770966bffad22adb2c7b74e93280c 100644
--- a/test/test_operators/test_smoothing_operator.py
+++ b/test/test_operators/test_smoothing_operator.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import list2fixture, setup_function, teardown_function
 
diff --git a/test/test_operators/test_value_inserter.py b/test/test_operators/test_value_inserter.py
index 5d15e27f1f2c72c7ca183378dac3d554a0fd31af..1bfb45bd68d018d80013f25f6aa97d5cf3db4389 100644
--- a/test/test_operators/test_value_inserter.py
+++ b/test/test_operators/test_value_inserter.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_
 
-import nifty6 as ift
+import nifty7 as ift
 from ..common import setup_function, teardown_function
 
 
diff --git a/test/test_operators/test_vdot_operator.py b/test/test_operators/test_vdot_operator.py
index c5394f517a0a69831eb9d4244d74c3ec7a23af78..140b43e856886866c30a186c5f193495340a9806 100644
--- a/test/test_operators/test_vdot_operator.py
+++ b/test/test_operators/test_vdot_operator.py
@@ -18,7 +18,7 @@
 import numpy as np
 from numpy.testing import assert_allclose
 
-import nifty6 as ift
+import nifty7 as ift
 
 from ..common import setup_function, teardown_function
 
@@ -36,4 +36,4 @@ def test_vdot_operator():
     res3 = np.vdot(arr1, arr2)
     assert_allclose(res1.val, res2.val)
     assert_allclose(res1.val, res3)
-    ift.extra.check_jacobian_consistency(op1, f)
+    ift.extra.check_operator(op1, f)
diff --git a/test/test_plot.py b/test/test_plot.py
index 239f7a25bd0a574c424255a06500c9c11fd00e1b..f2fcbfeba846e044d75dd5434761aad54480b772 100644
--- a/test/test_plot.py
+++ b/test/test_plot.py
@@ -19,7 +19,7 @@ from itertools import count
 
 import numpy as np
 
-import nifty6 as ift
+import nifty7 as ift
 
 from .common import setup_function, teardown_function
 
diff --git a/test/test_random.py b/test/test_random.py
index 62d48b1984c1cc45902281dc554148e06aca4980..21d576f631a9b713e8df9c87d53eaf9eaf2be654 100644
--- a/test/test_random.py
+++ b/test/test_random.py
@@ -17,7 +17,7 @@
 
 import numpy as np
 
-import nifty6 as ift
+import nifty7 as ift
 
 
 def check_state_back_to_orig():
diff --git a/test/test_spaces/test_gl_space.py b/test/test_spaces/test_gl_space.py
index b5ac1eec0559b6860a90d50c3b3a80746438c64c..7536e2f490802c83447f9db318d1e8c7cf4a3498 100644
--- a/test/test_spaces/test_gl_space.py
+++ b/test/test_spaces/test_gl_space.py
@@ -22,7 +22,7 @@ import pytest
 from numpy.testing import (assert_, assert_almost_equal, assert_equal,
                            assert_raises)
 
-from nifty6 import GLSpace
+from nifty7 import GLSpace
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
diff --git a/test/test_spaces/test_hp_space.py b/test/test_spaces/test_hp_space.py
index 363d815e7d168691c6ae2a37d3a7c446533d6745..24dee06f1a6cad72799d84f0260c3b7e71b4d1fd 100644
--- a/test/test_spaces/test_hp_space.py
+++ b/test/test_spaces/test_hp_space.py
@@ -20,7 +20,7 @@ import pytest
 from numpy.testing import (assert_, assert_almost_equal, assert_equal,
                            assert_raises)
 
-from nifty6 import HPSpace
+from nifty7 import HPSpace
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
diff --git a/test/test_spaces/test_interface.py b/test/test_spaces/test_interface.py
index 780e6cca0da2dc16601876519f1f7c44c4e974f3..815c8b0d05f2f79e69d67f5ce1ab93dfe6889b5b 100644
--- a/test/test_spaces/test_interface.py
+++ b/test/test_spaces/test_interface.py
@@ -20,7 +20,7 @@ from types import LambdaType
 import pytest
 from numpy.testing import assert_
 
-import nifty6 as ift
+import nifty7 as ift
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
diff --git a/test/test_spaces/test_lm_space.py b/test/test_spaces/test_lm_space.py
index 1c04a27b47dbec1a346733be81e5a379a2529296..e8ddc217679a168911cc3ef94ede46a7f4e289bb 100644
--- a/test/test_spaces/test_lm_space.py
+++ b/test/test_spaces/test_lm_space.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
 
-import nifty6 as ift
+import nifty7 as ift
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
diff --git a/test/test_spaces/test_power_space.py b/test/test_spaces/test_power_space.py
index a8b9979fa920950ab1808aed6605b9430e1a9872..a6eaf2007c526361fd7186b571f4aab0dc9d4a1e 100644
--- a/test/test_spaces/test_power_space.py
+++ b/test/test_spaces/test_power_space.py
@@ -21,7 +21,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
 
-import nifty6 as ift
+import nifty7 as ift
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
diff --git a/test/test_spaces/test_rg_space.py b/test/test_spaces/test_rg_space.py
index e95c5d8b9df473a3592d3bda311f77e2a75b8f92..788970a9528f56008b9bbf2257fe486c66612ec5 100644
--- a/test/test_spaces/test_rg_space.py
+++ b/test/test_spaces/test_rg_space.py
@@ -19,7 +19,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_, assert_allclose, assert_equal
 
-import nifty6 as ift
+import nifty7 as ift
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
diff --git a/test/test_sugar.py b/test/test_sugar.py
index 151b75572738ba95d4be1ad6fd128e8401f615c7..61b59a47e2cca99b0cf3abd80461bdc278a24276 100644
--- a/test/test_sugar.py
+++ b/test/test_sugar.py
@@ -16,11 +16,15 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 import numpy as np
+import pytest
 from numpy.testing import assert_equal
 
-import nifty6 as ift
+import nifty7 as ift
+
 from .common import setup_function, teardown_function
 
+pmp = pytest.mark.parametrize
+
 
 def test_get_signal_variance():
     space = ift.RGSpace(3)
@@ -45,15 +49,13 @@ def test_exec_time():
     lh = ift.GaussianEnergy(domain=op.target, sampling_dtype=np.float64) @ op1
     ic = ift.GradientNormController(iteration_limit=2)
     ham = ift.StandardHamiltonian(lh, ic_samp=ic)
-    kl = ift.MetricGaussianKL(ift.full(ham.domain, 0.), ham, 1)
+    kl = ift.MetricGaussianKL.make(ift.full(ham.domain, 0.), ham, 1)
     ops = [op, op1, lh, ham, kl]
     for oo in ops:
         for wm in [True, False]:
             ift.exec_time(oo, wm)
 
 
-import pytest
-pmp = pytest.mark.parametrize
 @pmp('mf', [False, True])
 @pmp('cplx', [False, True])
 def test_calc_pos(mf, cplx):