diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index da37b41c424dc58259bc22de8476ee6ce631696d..ca674370334fc6678dd2562ccc41df3516f41a17 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -5,11 +5,20 @@ variables:
   OMP_NUM_THREADS: 1
 
 stages:
+  - static_checks
   - build_docker
   - test
   - release
   - demo_runs
 
+check_no_asserts:
+  image: debian:testing-slim
+  stage: static_checks
+  before_script:
+    - ls
+  script:
+    - if [ `grep -r "^[[:space:]]*assert[ (]" src demos | wc -l` -ne 0 ]; then echo "Have found assert statements. Don't use them! Use \`utilities.myassert\` instead." && exit 1; fi
+
 build_docker_from_scratch:
   only:
     - schedules
diff --git a/src/extra.py b/src/extra.py
index 2a0512d9fbe2e21fa808cf7f3b16de76aa180327..85918eda948e6e8515006e845e41d87e0a92569c 100644
--- a/src/extra.py
+++ b/src/extra.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-2020 Max-Planck-Society
+# Copyright(C) 2013-2021 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -341,7 +341,7 @@ def _check_nontrivial_constant(op, loc, tol, ntries, only_r_differentiable,
 
         val0 = op(loc)
         _, op0 = op.simplify_for_constant_input(cstloc)
-        assert op0.domain is varloc.domain
+        myassert(op0.domain is varloc.domain)
         val1 = op0(varloc)
         assert_equal(val0, val1)
 
@@ -350,7 +350,7 @@ def _check_nontrivial_constant(op, loc, tol, ntries, only_r_differentiable,
         oplin0 = op0(lin0)
         oplin = op(lin)
 
-        assert oplin.jac.target is oplin0.jac.target
+        myassert(oplin.jac.target is oplin0.jac.target)
         rndinp = from_random(oplin.jac.target)
         assert_allclose(oplin.jac.adjoint(rndinp).extract(varloc.domain),
                         oplin0.jac.adjoint(rndinp), 1e-13, 1e-13)
diff --git a/src/library/correlated_fields.py b/src/library/correlated_fields.py
index b09fa2353c26650a8f7f7b2641d19e52d99eb5a3..a920020915564a5a28b4cc26cab4caea4dcbc19e 100644
--- a/src/library/correlated_fields.py
+++ b/src/library/correlated_fields.py
@@ -43,6 +43,7 @@ from ..operators.operator import Operator
 from ..operators.simple_linear_operators import VdotOperator, ducktape
 from ..probing import StatCalculator
 from ..sugar import full, makeDomain, makeField, makeOp
+from ..utilities import myassert
 
 
 def _log_k_lengths(pspace):
@@ -54,17 +55,17 @@ def _relative_log_k_lengths(power_space):
     """Log-distance to first bin
     logkl.shape==power_space.shape, logkl[0]=logkl[1]=0"""
     power_space = DomainTuple.make(power_space)
-    assert isinstance(power_space[0], PowerSpace)
-    assert len(power_space) == 1
+    myassert(isinstance(power_space[0], PowerSpace))
+    myassert(len(power_space) == 1)
     logkl = _log_k_lengths(power_space[0])
-    assert logkl.shape[0] == power_space[0].shape[0] - 1
+    myassert(logkl.shape[0] == power_space[0].shape[0] - 1)
     logkl -= logkl[0]
     return np.insert(logkl, 0, 0)
 
 
 def _log_vol(power_space):
     power_space = makeDomain(power_space)
-    assert isinstance(power_space[0], PowerSpace)
+    myassert(isinstance(power_space[0], PowerSpace))
     logk_lengths = _log_k_lengths(power_space[0])
     return logk_lengths[1:] - logk_lengths[:-1]
 
@@ -89,7 +90,7 @@ def _total_fluctuation_realized(samples):
 class _SlopeRemover(EndomorphicOperator):
     def __init__(self, domain, space=0):
         self._domain = makeDomain(domain)
-        assert isinstance(self._domain[space], PowerSpace)
+        myassert(isinstance(self._domain[space], PowerSpace))
         logkl = _relative_log_k_lengths(self._domain[space])
         sc = logkl/float(logkl[-1])
 
@@ -114,7 +115,7 @@ class _SlopeRemover(EndomorphicOperator):
 class _TwoLogIntegrations(LinearOperator):
     def __init__(self, target, space=0):
         self._target = makeDomain(target)
-        assert isinstance(self.target[space], PowerSpace)
+        myassert(isinstance(self.target[space], PowerSpace))
         dom = list(self._target)
         dom[space] = UnstructuredDomain((2, self.target[space].shape[0]-2))
         self._domain = makeDomain(dom)
@@ -173,7 +174,7 @@ class _Normalization(Operator):
     """
     def __init__(self, domain, space=0):
         self._domain = self._target = DomainTuple.make(domain)
-        assert isinstance(self._domain[space], PowerSpace)
+        myassert(isinstance(self._domain[space], PowerSpace))
         hspace = list(self._domain)
         hspace[space] = hspace[space].harmonic_partner
         hspace = makeDomain(hspace)
@@ -280,10 +281,10 @@ class _Amplitude(Operator):
         asperity > 0 or None
         loglogavgslope probably negative
         """
-        assert isinstance(fluctuations, Operator)
-        assert isinstance(flexibility, Operator) or flexibility is None
-        assert isinstance(asperity, Operator) or asperity is None
-        assert isinstance(loglogavgslope, Operator)
+        myassert(isinstance(fluctuations, Operator))
+        myassert(isinstance(flexibility, Operator) or flexibility is None)
+        myassert(isinstance(asperity, Operator) or asperity is None)
+        myassert(isinstance(loglogavgslope, Operator))
 
         if len(dofdex) > 0:
             N_copies = max(dofdex) + 1
@@ -296,7 +297,7 @@ class _Amplitude(Operator):
             N_copies = 0
             space = 0
             distributed_tgt = target = makeDomain(target)
-        assert isinstance(target[space], PowerSpace)
+        myassert(isinstance(target[space], PowerSpace))
 
         twolog = _TwoLogIntegrations(target, space)
         dom = twolog.domain
@@ -514,7 +515,7 @@ class CorrelatedFieldMaker:
         else:
             N = 0
             target_subdomain = makeDomain(target_subdomain)
-        # assert isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace))
+        # myassert(isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace)))
 
         for arg in [fluctuations, loglogavgslope]:
             if len(arg) != 2:
@@ -803,7 +804,7 @@ class CorrelatedFieldMaker:
             a_target = amp.target
             a_space = 0 if not hasattr(amp, "_space") else amp._space
             a_pp = amp.target[a_space]
-            assert isinstance(a_pp, PowerSpace)
+            myassert(isinstance(a_pp, PowerSpace))
 
             azm_expander = ContractionOperator(
                 a_target, spaces=a_space
diff --git a/src/minimization/line_search.py b/src/minimization/line_search.py
index 7eb93f89f2051a7adde223453021573f886aedce..76b6882c312f1d0f2a90665ca7736cf02a6f073d 100644
--- a/src/minimization/line_search.py
+++ b/src/minimization/line_search.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-2019, 2021 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -294,8 +294,8 @@ class LineSearch(metaclass=NiftyMeta):
         if phiprime_lo*(alpha_hi-alpha_lo) >= 0.:
             raise ValueError("inconsistent data")
         for i in range(self.max_zoom_iterations):
-            # assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0
-            # assert phiprime_lo*(alpha_hi-alpha_lo)<0.
+            # myassert(phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0)
+            # myassert(phiprime_lo*(alpha_hi-alpha_lo)<0.)
             delta_alpha = alpha_hi - alpha_lo
             a, b = min(alpha_lo, alpha_hi), max(alpha_lo, alpha_hi)
 
diff --git a/src/minimization/metric_gaussian_kl.py b/src/minimization/metric_gaussian_kl.py
index 6f6866425a37099ea853e205a413cd745b4a8bb4..dd8efc1b2b0bb34ae2daa0b2e46c5acb082fb3a8 100644
--- a/src/minimization/metric_gaussian_kl.py
+++ b/src/minimization/metric_gaussian_kl.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-2020 Max-Planck-Society
+# Copyright(C) 2013-2021 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -25,6 +25,7 @@ from ..operators.endomorphic_operator import EndomorphicOperator
 from ..operators.energy_operators import StandardHamiltonian
 from ..probing import approximation2endo
 from ..sugar import makeOp
+from ..utilities import myassert
 from .energy import Energy
 
 
@@ -47,9 +48,9 @@ def _get_lo_hi(comm, n_samples):
 def _modify_sample_domain(sample, domain):
     """Takes only keys from sample which are also in domain and inserts zeros
     for keys which are not in sample.domain."""
-    from ..multi_domain import MultiDomain
-    from ..field import Field
     from ..domain_tuple import DomainTuple
+    from ..field import Field
+    from ..multi_domain import MultiDomain
     from ..sugar import makeDomain
     domain = makeDomain(domain)
     if isinstance(domain, DomainTuple) and isinstance(sample, Field):
@@ -96,7 +97,7 @@ class MetricGaussianKL(Energy):
         if not _callingfrommake:
             raise NotImplementedError
         super(MetricGaussianKL, self).__init__(mean)
-        assert mean.domain is hamiltonian.domain
+        myassert(mean.domain is hamiltonian.domain)
         self._hamiltonian = hamiltonian
         self._n_samples = int(n_samples)
         self._mirror_samples = bool(mirror_samples)
diff --git a/src/operator_tree_optimiser.py b/src/operator_tree_optimiser.py
index a571054b71a0887c5e8f03809bb99b06ef59dabe..d986174003115f0d87a15b11cef7a15897ebfc50 100644
--- a/src/operator_tree_optimiser.py
+++ b/src/operator_tree_optimiser.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-2020 Max-Planck-Society
+# Copyright(C) 2013-2021 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -23,6 +23,7 @@ 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
+from .utilities import myassert
 
 
 def _optimise_operator(op):
@@ -312,7 +313,7 @@ def optimise_operator(op):
     test_field = from_random(op.domain)
     if isinstance(op(test_field), MultiField):
         for key in op(test_field).keys():
-            assert allclose(op(test_field).val[key], op_optimised(test_field).val[key], 1e-10)
+            myassert(allclose(op(test_field).val[key], op_optimised(test_field).val[key], 1e-10))
     else:
-        assert allclose(op(test_field).val, op_optimised(test_field).val, 1e-10)
+        myassert(allclose(op(test_field).val, op_optimised(test_field).val, 1e-10))
     return op_optimised
diff --git a/src/operators/einsum.py b/src/operators/einsum.py
index 230f7cea5c83c89dde5b7404f84c0b72ac86bf74..03aa39989fb47ca3175491333690eba594b280be 100644
--- a/src/operators/einsum.py
+++ b/src/operators/einsum.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-2020 Max-Planck-Society
+# Copyright(C) 2013-2021 Max-Planck-Society
 # Authors: Gordian Edenhofer, Philipp Frank
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
@@ -25,6 +25,7 @@ from ..field import Field
 from ..linearization import Linearization
 from ..multi_domain import MultiDomain
 from ..multi_field import MultiField
+from ..utilities import myassert
 from .linear_operator import LinearOperator
 from .operator import Operator
 
@@ -248,7 +249,7 @@ class LinearEinsum(LinearOperator):
                 if k_hit in _key_order:
                     tgt += [self._mf.domain[k_hit][dom_k_idx]]
                 else:
-                    assert k_hit == id(self)
+                    myassert(k_hit == id(self))
                     tgt += [self._domain[dom_k_idx]]
                 numpy_subscripts += "".join(subscriptmap[o])
             _target = DomainTuple.make(tgt)
diff --git a/src/operators/energy_operators.py b/src/operators/energy_operators.py
index 92b9bdd33e0fdfafd39cfbe5364223f561d7e320..2859e20120f74019b92d0c41f87d860976c253b7 100644
--- a/src/operators/energy_operators.py
+++ b/src/operators/energy_operators.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-2020 Max-Planck-Society
+# Copyright(C) 2013-2021 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -23,6 +23,7 @@ from ..field import Field
 from ..multi_domain import MultiDomain
 from ..multi_field import MultiField
 from ..sugar import makeDomain, makeOp
+from ..utilities import myassert
 from .linear_operator import LinearOperator
 from .operator import Operator
 from .sampling_enabler import SamplingDtypeSetter, SamplingEnabler
@@ -178,9 +179,9 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
 
     def _simplify_for_constant_input_nontrivial(self, c_inp):
         from .simplify_for_const import ConstantEnergyOperator
-        assert len(c_inp.keys()) == 1
+        myassert(len(c_inp.keys()) == 1)
         key = c_inp.keys()[0]
-        assert key in self._domain.keys()
+        myassert(key in self._domain.keys())
         cst = c_inp[key]
         if key == self._kr:
             res = _SpecialGammaEnergy(cst).ducktape(self._ki)
@@ -193,7 +194,7 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
                 trlog /= 2
             res = res + ConstantEnergyOperator(-trlog)
         res = res + ConstantEnergyOperator(0.)
-        assert res.target is self.target
+        myassert(res.target is self.target)
         return None, res
 
 
diff --git a/src/operators/operator.py b/src/operators/operator.py
index 965c33badff421d5f5bff9b646568550c26e66a5..dc9100075280b0a74e32d21d96e5ad038f5a66f5 100644
--- a/src/operators/operator.py
+++ b/src/operators/operator.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-2020 Max-Planck-Society
+# Copyright(C) 2013-2021 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -20,7 +20,7 @@ import numpy as np
 from .. import pointwise
 from ..logger import logger
 from ..multi_domain import MultiDomain
-from ..utilities import NiftyMeta, indent
+from ..utilities import NiftyMeta, indent, myassert
 
 
 class Operator(metaclass=NiftyMeta):
@@ -270,8 +270,8 @@ class Operator(metaclass=NiftyMeta):
         return self @ ducktape(self, None, name)
 
     def ducktape_left(self, name):
+        from ..sugar import is_fieldlike, is_linearization, is_operator
         from .simple_linear_operators import ducktape
-        from ..sugar import is_operator, is_fieldlike, is_linearization
         if is_operator(self):
             return ducktape(None, self, name) @ self
         if is_fieldlike(self) or is_linearization(self):
@@ -281,11 +281,12 @@ 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
-        from ..multi_field import MultiField
         from ..domain_tuple import DomainTuple
+        from ..multi_field import MultiField
         from ..sugar import makeDomain
+        from .energy_operators import EnergyOperator
+        from .simplify_for_const import (ConstantEnergyOperator,
+                                         ConstantOperator)
         if c_inp is None or (isinstance(c_inp, MultiField) and len(c_inp.keys()) == 0):
             return None, self
         dom = c_inp.domain
@@ -295,7 +296,7 @@ class Operator(metaclass=NiftyMeta):
         # 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)
+            myassert(isinstance(dom, MultiDomain))
             if not set(c_inp.keys()) <= set(self.domain.keys()):
                 raise ValueError
 
@@ -312,13 +313,13 @@ class Operator(metaclass=NiftyMeta):
         c_out, op = self._simplify_for_constant_input_nontrivial(c_inp)
         vardom = makeDomain({kk: vv for kk, vv in self.domain.items()
                              if kk not in c_inp.keys()})
-        assert op.domain is vardom
-        assert op.target is self.target
-        assert isinstance(op, Operator)
+        myassert(op.domain is vardom)
+        myassert(op.target is self.target)
+        myassert(isinstance(op, Operator))
         if c_out is not None:
-            assert isinstance(c_out, MultiField)
-            assert len(set(c_out.keys()) & self.domain.keys()) == 0
-            assert set(c_out.keys()) <= set(c_inp.keys())
+            myassert(isinstance(c_out, MultiField))
+            myassert(len(set(c_out.keys()) & self.domain.keys()) == 0)
+            myassert(set(c_out.keys()) <= set(c_inp.keys()))
         return c_out, op
 
     def _simplify_for_constant_input_nontrivial(self, c_inp):
diff --git a/src/operators/simplify_for_const.py b/src/operators/simplify_for_const.py
index 1f75371f8316969df692add024549aec92e65c14..1e68d6be5d9f79af94bb069bfa690c92624dd38a 100644
--- a/src/operators/simplify_for_const.py
+++ b/src/operators/simplify_for_const.py
@@ -11,11 +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) 2013-2020 Max-Planck-Society
+# Copyright(C) 2013-201 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 from ..multi_domain import MultiDomain
+from ..utilities import myassert
 from .block_diagonal_operator import BlockDiagonalOperator
 from .energy_operators import EnergyOperator
 from .operator import Operator
@@ -80,8 +81,8 @@ class ConstantOperator(Operator):
 
 class ConstantEnergyOperator(EnergyOperator):
     def __init__(self, output):
-        from ..sugar import makeDomain
         from ..field import Field
+        from ..sugar import makeDomain
         self._domain = makeDomain({})
         if not isinstance(output, Field):
             output = Field.scalar(float(output))
@@ -116,7 +117,7 @@ class InsertionOperator(Operator):
         self._jac = BlockDiagonalOperator(self._domain, jac) + NullOperator(makeDomain({}), cstdom)
 
     def apply(self, x):
-        assert len(set(self._cst.keys()) & set(x.domain.keys())) == 0
+        myassert(len(set(self._cst.keys()) & set(x.domain.keys())) == 0)
         val = x if x.jac is None else x.val
         val = val.unite(self._cst)
         if x.jac is None: