diff --git a/demos/bernoulli_demo.py b/demos/bernoulli_demo.py
index 76fec1a116c4f7d804ff108e6935645d6e4e06ef..bd727db9eca92b2589a7c207e3b5783a6e69a84a 100644
--- a/demos/bernoulli_demo.py
+++ b/demos/bernoulli_demo.py
@@ -43,7 +43,7 @@ if __name__ == '__main__':
     harmonic_space = position_space.get_default_codomain()
     HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
 
-    position = ift.from_random('normal', harmonic_space)
+    position = ift.from_random(harmonic_space, 'normal')
 
     # Define power spectrum and amplitudes
     def sqrtpspec(k):
@@ -58,13 +58,13 @@ if __name__ == '__main__':
 
     # Generate mock data
     p = R(sky)
-    mock_position = ift.from_random('normal', harmonic_space)
+    mock_position = ift.from_random(harmonic_space, 'normal')
     tmp = p(mock_position).val.astype(np.float64)
     data = ift.random.current_rng().binomial(1, tmp)
     data = ift.Field.from_raw(R.target, data)
 
     # Compute likelihood and Hamiltonian
-    position = ift.from_random('normal', harmonic_space)
+    position = ift.from_random(harmonic_space, 'normal')
     likelihood = ift.BernoulliEnergy(data) @ p
     ic_newton = ift.DeltaEnergyController(
         name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py
index c00bbe7436f5ad2ddcd14041db142536e5ec4e14..74b2a1b7c3308ae28c3af92300ebb31bf1a1d952 100644
--- a/demos/getting_started_1.py
+++ b/demos/getting_started_1.py
@@ -40,7 +40,7 @@ def make_checkerboard_mask(position_space):
 
 def make_random_mask():
     # Random mask for spherical mode
-    mask = ift.from_random('pm1', position_space)
+    mask = ift.from_random(position_space, 'pm1')
     mask = (mask + 1)/2
     return mask.val
 
diff --git a/demos/getting_started_2.py b/demos/getting_started_2.py
index 7856049d9efac5bb6ca51ce3745be05f0b5a446a..8b7c0cb6c3f781e54728e433df768ade93c42fcd 100644
--- a/demos/getting_started_2.py
+++ b/demos/getting_started_2.py
@@ -90,7 +90,7 @@ if __name__ == '__main__':
     # Generate mock data and define likelihood operator
     d_space = R.target[0]
     lamb = R(sky)
-    mock_position = ift.from_random('normal', domain)
+    mock_position = ift.from_random(domain, 'normal')
     data = lamb(mock_position)
     data = ift.random.current_rng().poisson(data.val.astype(np.float64))
     data = ift.Field.from_raw(d_space, data)
@@ -103,7 +103,7 @@ if __name__ == '__main__':
 
     # Compute MAP solution by minimizing the information Hamiltonian
     H = ift.StandardHamiltonian(likelihood)
-    initial_position = ift.from_random('normal', domain)
+    initial_position = ift.from_random(domain, 'normal')
     H = ift.EnergyAdapter(initial_position, H, want_metric=True)
     H, convergence = minimizer(H)
 
diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py
index 1d5ec727fce130eaa4e13bfbf5bbfba0926425c1..2dcffa17be798937832eee63edf9b2639d230fca 100644
--- a/demos/getting_started_3.py
+++ b/demos/getting_started_3.py
@@ -98,7 +98,7 @@ if __name__ == '__main__':
     N = ift.ScalingOperator(data_space, noise)
 
     # Generate mock signal and data
-    mock_position = ift.from_random('normal', signal_response.domain)
+    mock_position = ift.from_random(signal_response.domain, 'normal')
     data = signal_response(mock_position) + N.draw_sample_with_dtype(dtype=np.float64)
 
     # Minimization parameters
diff --git a/demos/getting_started_mf.py b/demos/getting_started_mf.py
index c66467cbd8bce9eb490fb533ce9226a228bb2106..217528dfbac99c9fed4d84dcc94de523df5404bf 100644
--- a/demos/getting_started_mf.py
+++ b/demos/getting_started_mf.py
@@ -97,7 +97,7 @@ if __name__ == '__main__':
     N = ift.ScalingOperator(data_space, noise)
 
     # Generate mock signal and data
-    mock_position = ift.from_random('normal', signal_response.domain)
+    mock_position = ift.from_random(signal_response.domain, 'normal')
     data = signal_response(mock_position) + N.draw_sample_with_dtype(dtype=np.float64)
 
     plot = ift.Plot()
diff --git a/nifty6/extra.py b/nifty6/extra.py
index e5a303fbdeb3aafad1cb2d3755170bead8d1e1c0..75f0849e86a1988db09b37905dcae15086f03cd4 100644
--- a/nifty6/extra.py
+++ b/nifty6/extra.py
@@ -42,8 +42,8 @@ def _adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol,
     needed_cap = op.TIMES | op.ADJOINT_TIMES
     if (op.capability & needed_cap) != needed_cap:
         return
-    f1 = from_random("normal", op.domain, dtype=domain_dtype)
-    f2 = from_random("normal", op.target, dtype=target_dtype)
+    f1 = from_random(op.domain, "normal", dtype=domain_dtype)
+    f2 = from_random(op.target, "normal", dtype=target_dtype)
     res1 = f1.s_vdot(op.adjoint_times(f2))
     res2 = op.times(f1).s_vdot(f2)
     if only_r_linear:
@@ -55,11 +55,11 @@ def _inverse_implementation(op, domain_dtype, target_dtype, atol, rtol):
     needed_cap = op.TIMES | op.INVERSE_TIMES
     if (op.capability & needed_cap) != needed_cap:
         return
-    foo = from_random("normal", op.target, dtype=target_dtype)
+    foo = from_random(op.target, "normal", dtype=target_dtype)
     res = op(op.inverse_times(foo))
     assert_allclose(res, foo, atol=atol, rtol=rtol)
 
-    foo = from_random("normal", op.domain, dtype=domain_dtype)
+    foo = from_random(op.domain, "normal", dtype=domain_dtype)
     res = op.inverse_times(op(foo))
     assert_allclose(res, foo, atol=atol, rtol=rtol)
 
@@ -75,8 +75,8 @@ def _check_linearity(op, domain_dtype, atol, rtol):
     needed_cap = op.TIMES
     if (op.capability & needed_cap) != needed_cap:
         return
-    fld1 = from_random("normal", op.domain, dtype=domain_dtype)
-    fld2 = from_random("normal", op.domain, dtype=domain_dtype)
+    fld1 = from_random(op.domain, "normal", dtype=domain_dtype)
+    fld2 = from_random(op.domain, "normal", dtype=domain_dtype)
     alpha = np.random.random()  # FIXME: this can break badly with MPI!
     val1 = op(alpha*fld1+fld2)
     val2 = alpha*op(fld1)+op(fld2)
@@ -88,7 +88,7 @@ def _actual_domain_check_linear(op, domain_dtype=None, inp=None):
     if (op.capability & needed_cap) != needed_cap:
         return
     if domain_dtype is not None:
-        inp = from_random("normal", op.domain, dtype=domain_dtype)
+        inp = from_random(op.domain, "normal", dtype=domain_dtype)
     elif inp is None:
         raise ValueError('Need to specify either dtype or inp')
     assert_(inp.domain is op.domain)
@@ -219,7 +219,7 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
 def _get_acceptable_location(op, loc, lin):
     if not np.isfinite(lin.val.s_sum()):
         raise ValueError('Initial value must be finite')
-    dir = from_random("normal", loc.domain)
+    dir = from_random(loc.domain, "normal")
     dirder = lin.jac(dir)
     if dirder.norm() == 0:
         dir = dir * (lin.val.norm()*1e-5)
diff --git a/nifty6/field.py b/nifty6/field.py
index f04a7a05ecae3e6772ec1d50962f931874a7d07b..e1cde6caf3c6143c6a9f911f23f94e5e48c4642d 100644
--- a/nifty6/field.py
+++ b/nifty6/field.py
@@ -124,7 +124,7 @@ class Field(Operator):
         return Field(DomainTuple.make(new_domain), self._val)
 
     @staticmethod
-    def from_random(random_type, domain, dtype=np.float64, **kwargs):
+    def from_random(domain, random_type='normal', dtype=np.float64, **kwargs):
         """Draws a random field with the given parameters.
 
         Parameters
diff --git a/nifty6/library/correlated_fields.py b/nifty6/library/correlated_fields.py
index 830fe8acfe51e72002330c6c6304ed6a74b6c224..fd0bee7b3de33c31c1192e3214879c71ee813b3d 100644
--- a/nifty6/library/correlated_fields.py
+++ b/nifty6/library/correlated_fields.py
@@ -524,7 +524,7 @@ class CorrelatedFieldMaker:
         for kk, op in lst:
             sc = StatCalculator()
             for _ in range(prior_info):
-                sc.add(op(from_random('normal', op.domain)))
+                sc.add(op(from_random(op.domain, 'normal')))
             mean = sc.mean.val
             stddev = sc.var.ptw("sqrt").val
             for m, s in zip(mean.flatten(), stddev.flatten()):
@@ -539,7 +539,7 @@ class CorrelatedFieldMaker:
         scm = 1.
         for a in self._a:
             op = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
-            res = np.array([op(from_random('normal', op.domain)).val
+            res = np.array([op(from_random(op.domain, 'normal')).val
                             for _ in range(nsamples)])
             scm *= res**2 + 1.
         return fluctuations_slice_mean/np.mean(np.sqrt(scm))
diff --git a/nifty6/multi_field.py b/nifty6/multi_field.py
index 6826607bd0b68efbe06ba42077966172d54e4050..50bcc8caab0da4e038e9aec42c36ca6a4faf7617 100644
--- a/nifty6/multi_field.py
+++ b/nifty6/multi_field.py
@@ -101,7 +101,7 @@ class MultiField(Operator):
         return self._transform(lambda x: x.imag)
 
     @staticmethod
-    def from_random(random_type, domain, dtype=np.float64, **kwargs):
+    def from_random(domain, random_type='normal', dtype=np.float64, **kwargs):
         """Draws a random multi-field with the given parameters.
 
         Parameters
@@ -131,7 +131,7 @@ class MultiField(Operator):
         else:
             dtype = np.dtype(dtype)
             dtype = {kk: dtype for kk in domain.keys()}
-        dct = {kk: Field.from_random(random_type, domain[kk], dtype[kk], **kwargs)
+        dct = {kk: Field.from_random(domain[kk], random_type, dtype[kk], **kwargs)
                for kk in domain.keys()}
         return MultiField.from_dict(dct)
 
diff --git a/nifty6/operators/block_diagonal_operator.py b/nifty6/operators/block_diagonal_operator.py
index 37023b97be9ae43c920216ff2a35f4d741c84cdf..34f6ffca36fd8ea3e0099d8cedb8d5171cb3a7dc 100644
--- a/nifty6/operators/block_diagonal_operator.py
+++ b/nifty6/operators/block_diagonal_operator.py
@@ -61,7 +61,7 @@ class BlockDiagonalOperator(EndomorphicOperator):
         val = tuple(
             op.draw_sample_with_dtype(dtype[key], from_inverse)
             if op is not None
-            else from_random('normal', self._domain[key], dtype=dtype)
+            else from_random(self._domain[key], 'normal', dtype=dtype)
             for op, key in zip(self._ops, self._domain.keys()))
         return MultiField(self._domain, val)
 
diff --git a/nifty6/operators/diagonal_operator.py b/nifty6/operators/diagonal_operator.py
index 21055550d2180e9f6d8b46380068019b6b44b691..7cb6e0c22f17b0cbe15f4477a710d41b53f60f1f 100644
--- a/nifty6/operators/diagonal_operator.py
+++ b/nifty6/operators/diagonal_operator.py
@@ -161,8 +161,7 @@ class DiagonalOperator(EndomorphicOperator):
         return Field(self._domain, res)
 
     def draw_sample_with_dtype(self, dtype, from_inverse=False):
-        res = Field.from_random(random_type="normal", domain=self._domain,
-                                dtype=dtype)
+        res = Field.from_random(domain=self._domain, random_type="normal", dtype=dtype)
         return self.process_sample(res, from_inverse)
 
     def __repr__(self):
diff --git a/nifty6/operators/domain_tuple_field_inserter.py b/nifty6/operators/domain_tuple_field_inserter.py
index 49909b2bff1304fd7f2daf498e53f1b27357f727..0809488545a1a470b6398966877e104a32222b37 100644
--- a/nifty6/operators/domain_tuple_field_inserter.py
+++ b/nifty6/operators/domain_tuple_field_inserter.py
@@ -35,7 +35,7 @@ class DomainTupleFieldInserter(LinearOperator):
         Slice in new sub-domain in which the input field shall be written into.
     """
 
-    def __init__(self, target, space, pos):
+    def __init__(self, target, space, index):
         if not space <= len(target) or space < 0:
             raise ValueError
         self._target = DomainTuple.make(target)
@@ -48,13 +48,13 @@ class DomainTupleFieldInserter(LinearOperator):
         nshp = new_space.shape
         fst_dims = sum(len(dd.shape) for dd in self.target[:space])
 
-        if len(pos) != len(nshp):
+        if len(index) != len(nshp):
             raise ValueError("shape mismatch between new_space and position")
-        for s, p in zip(nshp, pos):
+        for s, p in zip(nshp, index):
             if p < 0 or p >= s:
                 raise ValueError("bad position value")
 
-        self._slc = (slice(None),)*fst_dims + pos
+        self._slc = (slice(None),)*fst_dims + index
 
     def apply(self, x, mode):
         self._check_input(x, mode)
diff --git a/nifty6/operators/scaling_operator.py b/nifty6/operators/scaling_operator.py
index 29573ace69b435c80a60e825039fcc6519838381..978ba1e271974849d6979d0bdc61969948a7aacb 100644
--- a/nifty6/operators/scaling_operator.py
+++ b/nifty6/operators/scaling_operator.py
@@ -93,8 +93,7 @@ class ScalingOperator(EndomorphicOperator):
 
     def draw_sample_with_dtype(self, dtype, from_inverse=False):
         from ..sugar import from_random
-        return from_random(random_type="normal", domain=self._domain,
-                           std=self._get_fct(from_inverse), dtype=dtype)
+        return from_random(domain=self._domain, random_type="normal", dtype=dtype, std=self._get_fct(from_inverse))
 
     def __call__(self, other):
         res = EndomorphicOperator.__call__(self, other)
diff --git a/nifty6/probing.py b/nifty6/probing.py
index 7aa43402ca2f6647cc4da5d53c1a3cb3a583c131..46968ae8c2d26d8953947e338bac032dc22576b0 100644
--- a/nifty6/probing.py
+++ b/nifty6/probing.py
@@ -134,7 +134,7 @@ def probe_diagonal(op, nprobes, random_type="pm1"):
     '''
     sc = StatCalculator()
     for i in range(nprobes):
-        x = from_random(random_type, op.domain)
+        x = from_random(op.domain, random_type)
         sc.add(op(x).conjugate()*x)
     return sc.mean
 
diff --git a/nifty6/sugar.py b/nifty6/sugar.py
index 65dacdcb1c02e8618bf2f1828c6e944146eef5da..719cccd51d397cf48f8dac52c85c5dde2ae1d79b 100644
--- a/nifty6/sugar.py
+++ b/nifty6/sugar.py
@@ -256,7 +256,7 @@ def full(domain, val):
     return Field.full(domain, val)
 
 
-def from_random(random_type, domain, dtype=np.float64, **kwargs):
+def from_random(domain, random_type='normal', dtype=np.float64, **kwargs):
     """Convenience function creating Fields/MultiFields with random values.
 
     Parameters
@@ -283,8 +283,8 @@ def from_random(random_type, domain, dtype=np.float64, **kwargs):
     random numbers, even for the same initial RNG state.
     """
     if isinstance(domain, (dict, MultiDomain)):
-        return MultiField.from_random(random_type, domain, dtype, **kwargs)
-    return Field.from_random(random_type, domain, dtype, **kwargs)
+        return MultiField.from_random(domain, random_type, dtype, **kwargs)
+    return Field.from_random(domain, random_type, dtype, **kwargs)
 
 
 def makeField(domain, arr):
@@ -466,7 +466,7 @@ def exec_time(obj, want_metric=True):
         logger.info('Energy.metric(position): {}'.format(time() - t0))
     elif isinstance(obj, Operator):
         want_metric = bool(want_metric)
-        pos = from_random('normal', obj.domain)
+        pos = from_random(obj.domain, 'normal')
         t0 = time()
         obj(pos)
         logger.info('Operator call with field: {}'.format(time() - t0))
@@ -509,7 +509,7 @@ def calculate_position(operator, output):
     lh = GaussianEnergy(d, invcov) @ operator
     H = StandardHamiltonian(
         lh, ic_samp=GradientNormController(iteration_limit=200))
-    pos = 0.1*from_random('normal', operator.domain)
+    pos = 0.1 * from_random(operator.domain, 'normal')
     minimizer = NewtonCG(GradientNormController(iteration_limit=10))
     for ii in range(3):
         kl = MetricGaussianKL(pos, H, 3, mirror_samples=True)
diff --git a/test/test_energy_gradients.py b/test/test_energy_gradients.py
index 32ce00e24601d2f663a2c0a56ab89bd8d71e1294..41e582a24c8f90507913b09fc826087e475a23f3 100644
--- a/test/test_energy_gradients.py
+++ b/test/test_energy_gradients.py
@@ -26,7 +26,7 @@ spaces = [ift.GLSpace(5),
           ift.MultiDomain.make({'': ift.RGSpace(5, distances=.789)}),
           (ift.RGSpace(3, distances=.789), ift.UnstructuredDomain(2))]
 pmp = pytest.mark.parametrize
-field = list2fixture([ift.from_random('normal', sp) for sp in spaces])
+field = list2fixture([ift.from_random(sp, 'normal') for sp in spaces])
 ntries = 10
 
 
@@ -62,7 +62,7 @@ def test_studentt(field):
         return
     energy = ift.StudentTEnergy(domain=field.domain, theta=.5)
     ift.extra.check_jacobian_consistency(energy, field, tol=1e-6)
-    theta = ift.from_random('normal',field.domain).exp()
+    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)
 
@@ -73,7 +73,7 @@ def test_hamiltonian_and_KL(field):
     lh = ift.GaussianEnergy(domain=space)
     hamiltonian = ift.StandardHamiltonian(lh)
     ift.extra.check_jacobian_consistency(hamiltonian, field, ntries=ntries)
-    samps = [ift.from_random('normal', space) for i in range(2)]
+    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)
 
diff --git a/test/test_field.py b/test/test_field.py
index 1394cbee46e6706481ffe2e267197574b1313e3e..6006f74f5424ab74dfe09f9cc7d0e53af9707c49 100644
--- a/test/test_field.py
+++ b/test/test_field.py
@@ -104,7 +104,7 @@ def test_DiagonalOperator_power_analyze2(space1, space2):
     ift.RGSpace((2, 3, 7))
 ])
 def test_norm(space):
-    f = ift.Field.from_random("normal", domain=space, dtype=np.complex128)
+    f = ift.Field.from_random(domain=space, random_type="normal", dtype=np.complex128)
     gd = f.val.reshape(-1)
     assert_allclose(f.norm(), np.linalg.norm(gd))
     assert_allclose(f.norm(1), np.linalg.norm(gd, ord=1))
@@ -115,8 +115,8 @@ def test_norm(space):
 
 def test_vdot():
     s = ift.RGSpace((10,))
-    f1 = ift.Field.from_random("normal", domain=s, dtype=np.complex128)
-    f2 = ift.Field.from_random("normal", domain=s, dtype=np.complex128)
+    f1 = ift.Field.from_random(domain=s, random_type="normal", dtype=np.complex128)
+    f2 = ift.Field.from_random(domain=s, random_type="normal", dtype=np.complex128)
     assert_allclose(f1.s_vdot(f2), f1.vdot(f2, spaces=0).val)
     assert_allclose(f1.s_vdot(f2), np.conj(f2.s_vdot(f1)))
 
@@ -297,8 +297,8 @@ def test_stdfunc():
     fx = ift.full(f.domain, 67.)
     assert_equal(f.shape, fx.shape)
     assert_equal(fx.val, 67.)
-    f = ift.Field.from_random("normal", s)
-    f2 = ift.Field.from_random("normal", s)
+    f = ift.Field.from_random(s, "normal")
+    f2 = ift.Field.from_random(s, "normal")
     assert_equal((f > f2).val, f.val > f2.val)
     assert_equal((f >= f2).val, f.val >= f2.val)
     assert_equal((f < f2).val, f.val < f2.val)
@@ -341,7 +341,7 @@ def test_funcs(num, dom, func):
 @pmp('dtype', [np.float64, np.complex128])
 def test_from_random(rtype, dtype):
     sp = ift.RGSpace(3)
-    ift.Field.from_random(rtype, sp, dtype=dtype)
+    ift.Field.from_random(sp, rtype, dtype=dtype)
 
 
 def test_field_of_objects():
diff --git a/test/test_kl.py b/test/test_kl.py
index d0c656a8aeeacd25690e6056b43f4510d91ffab9..75cb0276eb5ac37073f941eb550ae5c93a25ce64 100644
--- a/test/test_kl.py
+++ b/test/test_kl.py
@@ -38,7 +38,7 @@ def test_kl(constants, point_estimates, mirror_samples, mf):
     lh = ift.GaussianEnergy(domain=op.target, sampling_dtype=np.float64) @ op
     ic = ift.GradientNormController(iteration_limit=5)
     h = ift.StandardHamiltonian(lh, ic_samp=ic)
-    mean0 = ift.from_random('normal', h.domain)
+    mean0 = ift.from_random(h.domain, 'normal')
 
     nsamps = 2
     kl = ift.MetricGaussianKL(mean0,
diff --git a/test/test_minimizers.py b/test/test_minimizers.py
index df537c723651427cf4f708d926c37457d84a1d40..3c847fd5add15881ca03e8fd3f12b70010133d55 100644
--- a/test/test_minimizers.py
+++ b/test/test_minimizers.py
@@ -52,8 +52,8 @@ slow_minimizers = ['ift.SteepestDescent(IC)']
      slow_minimizers)
 @pmp('space', spaces)
 def test_quadratic_minimization(minimizer, space):
-    starting_point = ift.Field.from_random('normal', domain=space)*10
-    covariance_diagonal = ift.Field.from_random('uniform', domain=space) + 0.5
+    starting_point = ift.Field.from_random(domain=space, random_type='normal') * 10
+    covariance_diagonal = ift.Field.from_random(domain=space, random_type='uniform') + 0.5
     covariance = ift.DiagonalOperator(covariance_diagonal)
     required_result = ift.full(space, 1.)
 
@@ -78,11 +78,11 @@ def test_quadratic_minimization(minimizer, space):
 def test_WF_curvature(space):
     required_result = ift.full(space, 1.)
 
-    s = ift.Field.from_random('uniform', domain=space) + 0.5
+    s = ift.Field.from_random(domain=space, random_type='uniform') + 0.5
     S = ift.DiagonalOperator(s)
-    r = ift.Field.from_random('uniform', domain=space)
+    r = ift.Field.from_random(domain=space, random_type='uniform')
     R = ift.DiagonalOperator(r)
-    n = ift.Field.from_random('uniform', domain=space) + 0.5
+    n = ift.Field.from_random(domain=space, random_type='uniform') + 0.5
     N = ift.DiagonalOperator(n)
     all_diag = 1./s + r**2/n
     curv = ift.WienerFilterCurvature(R, N, S, iteration_controller=IC,
@@ -100,7 +100,7 @@ def test_WF_curvature(space):
 
     if len(space.shape) == 1:
         R = ift.ValueInserter(space, [0])
-        n = ift.from_random('uniform', R.domain) + 0.5
+        n = ift.from_random(R.domain, 'uniform') + 0.5
         N = ift.DiagonalOperator(n)
         all_diag = 1./s + R(1/n)
         curv = ift.WienerFilterCurvature(R.adjoint, N, S,
@@ -125,7 +125,7 @@ def test_rosenbrock(minimizer):
     except ImportError:
         raise SkipTest
     space = ift.DomainTuple.make(ift.UnstructuredDomain((2,)))
-    starting_point = ift.Field.from_random('normal', domain=space)*10
+    starting_point = ift.Field.from_random(domain=space, random_type='normal') * 10
 
     class RBEnergy(ift.Energy):
         def __init__(self, position):
diff --git a/test/test_mpi/test_kl.py b/test/test_mpi/test_kl.py
index ccd713641424425a06b232169984c42eece89a01..36a2a68a5e5ba9039fdcbb522dd0e5aaa4a13402 100644
--- a/test/test_mpi/test_kl.py
+++ b/test/test_mpi/test_kl.py
@@ -50,7 +50,7 @@ def test_kl(constants, point_estimates, mirror_samples, mode, mf):
     lh = ift.GaussianEnergy(domain=op.target, sampling_dtype=np.float64) @ op
     ic = ift.GradientNormController(iteration_limit=5)
     h = ift.StandardHamiltonian(lh, ic_samp=ic)
-    mean0 = ift.from_random('normal', h.domain)
+    mean0 = ift.from_random(h.domain, 'normal')
     nsamps = 2
     args = {'constants': constants,
             'point_estimates': point_estimates,
diff --git a/test/test_multi_field.py b/test/test_multi_field.py
index 9469596c9cc155c81060c189de4309d0121425e1..16f26e2e7ec61baad5f9a36ebdfd9157960fb3c7 100644
--- a/test/test_multi_field.py
+++ b/test/test_multi_field.py
@@ -25,13 +25,13 @@ dom = ift.makeDomain({"d1": ift.RGSpace(10)})
 
 
 def test_vdot():
-    f1 = ift.from_random("normal", domain=dom, dtype=np.complex128)
-    f2 = ift.from_random("normal", domain=dom, dtype=np.complex128)
+    f1 = ift.from_random(domain=dom, random_type="normal", dtype=np.complex128)
+    f2 = ift.from_random(domain=dom, random_type="normal", dtype=np.complex128)
     assert_allclose(f1.s_vdot(f2), np.conj(f2.s_vdot(f1)))
 
 
 def test_func():
-    f1 = ift.from_random("normal", domain=dom, dtype=np.complex128)
+    f1 = ift.from_random(domain=dom, random_type="normal", dtype=np.complex128)
     assert_allclose(
         f1.ptw("exp").ptw("log")["d1"].val, f1["d1"].val)
 
diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py
index f605df6aba2481389df1e021997a6d1f67c7a088..c410fc71627792f79ef616ed77ac836ce3dd1531 100644
--- a/test/test_operators/test_adjoint.py
+++ b/test/test_operators/test_adjoint.py
@@ -54,8 +54,8 @@ def testLOSResponse(sp, dtype):
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 def testOperatorCombinations(sp, dtype):
-    a = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
-    b = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
+    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)
     op = a(b)
@@ -91,7 +91,7 @@ def testConjugationOperator(sp):
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 def testOperatorAdaptor(sp, dtype):
-    op = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=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)
@@ -169,7 +169,7 @@ def testHarmonic(sp, dtype):
 @pmp('sp', _p_spaces)
 def testMask(sp, dtype):
     # Create mask
-    f = ift.from_random('normal', sp).val
+    f = ift.from_random(sp, 'normal').val
     mask = np.zeros_like(f)
     mask[f > 0] = 1
     mask = ift.Field.from_raw(sp, mask)
@@ -180,7 +180,7 @@ def testMask(sp, dtype):
 
 @pmp('sp', _h_spaces + _p_spaces)
 def testDiagonal(sp, dtype):
-    op = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
+    op = ift.DiagonalOperator(ift.Field.from_random(sp, "normal", dtype=dtype))
     ift.extra.consistency_check(op, dtype, dtype)
 
 
@@ -240,8 +240,8 @@ def testRegridding(args):
     ift.DomainTuple.make(ift.RGSpace((10, 12), distances=(0.1, 1.)),)
 ])
 def testOuter(fdomain, domain):
-    f = ift.from_random('normal', fdomain)
-    op = ift.OuterProduct(f, domain)
+    f = ift.from_random(fdomain, 'normal')
+    op = ift.OuterProduct(domain, f)
     ift.extra.consistency_check(op)
 
 
diff --git a/test/test_operators/test_composed_operator.py b/test/test_operators/test_composed_operator.py
index 63169c27954049f6b325c9bc8e751a5c9b310768..55b517a7648160875591ff9e1ad5aab7cf03442c 100644
--- a/test/test_operators/test_composed_operator.py
+++ b/test/test_operators/test_composed_operator.py
@@ -33,15 +33,15 @@ space2 = space1
 
 def test_times_adjoint_times(space1, space2):
     cspace = (space1, space2)
-    diag1 = ift.Field.from_random('normal', domain=space1)
-    diag2 = ift.Field.from_random('normal', domain=space2)
+    diag1 = ift.Field.from_random(domain=space1, random_type='normal')
+    diag2 = ift.Field.from_random(domain=space2, random_type='normal')
     op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
     op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
 
     op = op2(op1)
 
-    rand1 = ift.Field.from_random('normal', domain=(space1, space2))
-    rand2 = ift.Field.from_random('normal', domain=(space1, space2))
+    rand1 = ift.Field.from_random(domain=(space1, space2), random_type='normal')
+    rand2 = ift.Field.from_random(domain=(space1, space2), random_type='normal')
 
     tt1 = rand2.s_vdot(op.times(rand1))
     tt2 = rand1.s_vdot(op.adjoint_times(rand2))
@@ -50,14 +50,14 @@ def test_times_adjoint_times(space1, space2):
 
 def test_times_inverse_times(space1, space2):
     cspace = (space1, space2)
-    diag1 = ift.Field.from_random('normal', domain=space1)
-    diag2 = ift.Field.from_random('normal', domain=space2)
+    diag1 = ift.Field.from_random(domain=space1, random_type='normal')
+    diag2 = ift.Field.from_random(domain=space2, random_type='normal')
     op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
     op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
 
     op = op2(op1)
 
-    rand1 = ift.Field.from_random('normal', domain=(space1, space2))
+    rand1 = ift.Field.from_random(domain=(space1, space2), random_type='normal')
     tt1 = op.inverse_times(op.times(rand1))
 
     assert_allclose(tt1.val, rand1.val)
diff --git a/test/test_operators/test_convolution_operators.py b/test/test_operators/test_convolution_operators.py
index 52497023c3b003d50b1b464e9db5af9bb62ee415..6f4f6495c044862241fe6883b6fba76ca5be3e8b 100644
--- a/test/test_operators/test_convolution_operators.py
+++ b/test/test_operators/test_convolution_operators.py
@@ -30,7 +30,7 @@ space = list2fixture([
 
 
 def test_const_func(space):
-    sig = ift.Field.from_random('normal', domain=space)
+    sig = ift.Field.from_random(domain=space, random_type='normal')
     fco_op = ift.FuncConvolutionOperator(space, lambda x: np.ones(x.shape))
     vals = fco_op(sig).val
     vals = np.round(vals, decimals=5)
@@ -46,7 +46,7 @@ def test_gaussian_smoothing():
     N = 128
     sigma = N / 10**4
     dom = ift.RGSpace(N)
-    sig = ift.Field.from_random('normal', dom).ptw("exp")
+    sig = ift.Field.from_random(dom, 'normal').ptw("exp")
     fco_op = ift.FuncConvolutionOperator(dom, lambda x: gauss(x, sigma))
     sm_op = ift.HarmonicSmoothingOperator(dom, sigma)
     assert_allclose(fco_op(sig).val,
diff --git a/test/test_operators/test_correlated_fields.py b/test/test_operators/test_correlated_fields.py
index f9985133770b2840d9d5f9b3d21f02711c8c218e..12f3008cc9c82b44d2e4ed4cd2b0ea8984b0b3de 100644
--- a/test/test_operators/test_correlated_fields.py
+++ b/test/test_operators/test_correlated_fields.py
@@ -52,7 +52,7 @@ def testAmplitudesInvariants(sspace, N):
                         'freq', dofdex=dofdex3)
     op = fa.finalize()
 
-    samples = [ift.from_random('normal', op.domain) for _ in range(100)]
+    samples = [ift.from_random(op.domain, 'normal') for _ in range(100)]
     tot_flm, _ = _stats(fa.total_fluctuation, samples)
     offset_amp_std, _ = _stats(fa.amplitude_total_offset, samples)
     intergated_fluct_std0, _ = _stats(fa.average_fluctuation(0), samples)
diff --git a/test/test_operators/test_diagonal_operator.py b/test/test_operators/test_diagonal_operator.py
index 82a74de38b84ea3f2ea376629f485801f0f7e188..f165a22731e6804a12d928f9ad842554eddd4eb4 100644
--- a/test/test_operators/test_diagonal_operator.py
+++ b/test/test_operators/test_diagonal_operator.py
@@ -31,16 +31,16 @@ space = list2fixture([
 
 
 def test_property(space):
-    diag = ift.Field.from_random('normal', domain=space)
+    diag = ift.Field.from_random(domain=space, random_type='normal')
     D = ift.DiagonalOperator(diag)
     if D.domain[0] != space:
         raise TypeError
 
 
 def test_times_adjoint(space):
-    rand1 = ift.Field.from_random('normal', domain=space)
-    rand2 = ift.Field.from_random('normal', domain=space)
-    diag = ift.Field.from_random('normal', domain=space)
+    rand1 = ift.Field.from_random(domain=space, random_type='normal')
+    rand2 = ift.Field.from_random(domain=space, random_type='normal')
+    diag = ift.Field.from_random(domain=space, random_type='normal')
     D = ift.DiagonalOperator(diag)
     tt1 = rand1.s_vdot(D.times(rand2))
     tt2 = rand2.s_vdot(D.times(rand1))
@@ -48,47 +48,47 @@ def test_times_adjoint(space):
 
 
 def test_times_inverse(space):
-    rand1 = ift.Field.from_random('normal', domain=space)
-    diag = ift.Field.from_random('normal', domain=space)
+    rand1 = ift.Field.from_random(domain=space, random_type='normal')
+    diag = ift.Field.from_random(domain=space, random_type='normal')
     D = ift.DiagonalOperator(diag)
     tt1 = D.times(D.inverse_times(rand1))
     assert_allclose(rand1.val, tt1.val)
 
 
 def test_times(space):
-    rand1 = ift.Field.from_random('normal', domain=space)
-    diag = ift.Field.from_random('normal', domain=space)
+    rand1 = ift.Field.from_random(domain=space, random_type='normal')
+    diag = ift.Field.from_random(domain=space, random_type='normal')
     D = ift.DiagonalOperator(diag)
     tt = D.times(rand1)
     assert_equal(tt.domain[0], space)
 
 
 def test_adjoint_times(space):
-    rand1 = ift.Field.from_random('normal', domain=space)
-    diag = ift.Field.from_random('normal', domain=space)
+    rand1 = ift.Field.from_random(domain=space, random_type='normal')
+    diag = ift.Field.from_random(domain=space, random_type='normal')
     D = ift.DiagonalOperator(diag)
     tt = D.adjoint_times(rand1)
     assert_equal(tt.domain[0], space)
 
 
 def test_inverse_times(space):
-    rand1 = ift.Field.from_random('normal', domain=space)
-    diag = ift.Field.from_random('normal', domain=space)
+    rand1 = ift.Field.from_random(domain=space, random_type='normal')
+    diag = ift.Field.from_random(domain=space, random_type='normal')
     D = ift.DiagonalOperator(diag)
     tt = D.inverse_times(rand1)
     assert_equal(tt.domain[0], space)
 
 
 def test_adjoint_inverse_times(space):
-    rand1 = ift.Field.from_random('normal', domain=space)
-    diag = ift.Field.from_random('normal', domain=space)
+    rand1 = ift.Field.from_random(domain=space, random_type='normal')
+    diag = ift.Field.from_random(domain=space, random_type='normal')
     D = ift.DiagonalOperator(diag)
     tt = D.adjoint_inverse_times(rand1)
     assert_equal(tt.domain[0], space)
 
 
 def test_diagonal(space):
-    diag = ift.Field.from_random('normal', domain=space)
+    diag = ift.Field.from_random(domain=space, random_type='normal')
     D = ift.DiagonalOperator(diag)
     diag_op = D(ift.Field.full(space, 1.))
     assert_allclose(diag.val, diag_op.val)
diff --git a/test/test_operators/test_einsum.py b/test/test_operators/test_einsum.py
index 0f732131ae9aa95f14e6e4ada2a8fa225386d2e1..e2b67a6afa7df75cacbcb01b44d3ff8ef1e226fc 100644
--- a/test/test_operators/test_einsum.py
+++ b/test/test_operators/test_einsum.py
@@ -38,20 +38,20 @@ def test_linear_einsum_outer(space1, space2, dtype, n_invocations=10):
     mf_dom = ift.MultiDomain.make({
         "dom01": space1,
         "dom02": ift.DomainTuple.make((space1, space2))})
-    mf = ift.from_random("normal", mf_dom, dtype=dtype)
+    mf = ift.from_random(mf_dom, "normal", dtype=dtype)
     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)
 
     le_ift = ift.DiagonalOperator(mf["dom01"], domain=mf_dom["dom02"], spaces=0) @ ift.DiagonalOperator(mf["dom02"])
-    le_ift = le_ift @ ift.OuterProduct(ift.full(mf_dom["dom01"], 1.),
-                                       ift.DomainTuple.make(mf_dom["dom02"][1]))
+    le_ift = le_ift @ ift.OuterProduct(ift.DomainTuple.make(mf_dom["dom02"][1]),
+                                       ift.full(mf_dom["dom01"], 1.))
 
     for _ in range(n_invocations):
-        r = ift.from_random("normal", le.domain, dtype=dtype)
+        r = ift.from_random(le.domain, "normal", dtype=dtype)
         assert_allclose(le(r).val, le_ift(r).val)
-        r_adj = ift.from_random("normal", le.target, dtype=dtype)
+        r_adj = ift.from_random(le.target, "normal", dtype=dtype)
         assert_allclose(le.adjoint(r_adj).val, le_ift.adjoint(r_adj).val)
 
 
@@ -59,7 +59,7 @@ def test_linear_einsum_contraction(space1, space2, dtype, n_invocations=10):
     mf_dom = ift.MultiDomain.make({
         "dom01": space1,
         "dom02": ift.DomainTuple.make((space1, space2))})
-    mf = ift.from_random("normal", mf_dom, dtype=dtype)
+    mf = ift.from_random(mf_dom, "normal", dtype=dtype)
     ss = "i,ij,j->i"
     key_order = ("dom01", "dom02")
     le = ift.LinearEinsum(space2, mf, ss, key_order=key_order)
@@ -68,13 +68,13 @@ def test_linear_einsum_contraction(space1, space2, dtype, n_invocations=10):
     le_ift = ift.ContractionOperator(mf_dom["dom02"], 1)
     le_ift = le_ift @ ift.DiagonalOperator(mf["dom01"], domain=mf_dom["dom02"], spaces=0)
     le_ift = le_ift @ ift.DiagonalOperator(mf["dom02"])
-    le_ift = le_ift @ ift.OuterProduct(ift.full(mf_dom["dom01"], 1.),
-                                       ift.DomainTuple.make(mf_dom["dom02"][1]))
+    le_ift = le_ift @ ift.OuterProduct(ift.DomainTuple.make(mf_dom["dom02"][1]),
+                                       ift.full(mf_dom["dom01"], 1.),)
 
     for _ in range(n_invocations):
-        r = ift.from_random("normal", le.domain, dtype=dtype)
+        r = ift.from_random(le.domain, "normal", dtype=dtype)
         assert_allclose(le(r).val, le_ift(r).val)
-        r_adj = ift.from_random("normal", le.target, dtype=dtype)
+        r_adj = ift.from_random(le.target, "normal", dtype=dtype)
         assert_allclose(le.adjoint(r_adj).val, le_ift.adjoint(r_adj).val)
 
 
@@ -121,13 +121,13 @@ 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("normal", mle.domain, dtype=dtype), ntries=ntries)
+    check_jacobian_consistency(mle, ift.from_random(mle.domain, "normal", dtype=dtype), ntries=ntries)
 
     outer_i = ift.OuterProduct(
-        ift.full(mf_dom["dom03"], 1.), ift.DomainTuple.make(mf_dom["dom02"][0])
+        ift.DomainTuple.make(mf_dom["dom02"][0]), ift.full(mf_dom["dom03"], 1.)
     )
     outer_j = ift.OuterProduct(
-        ift.full(mf_dom["dom01"], 1.), ift.DomainTuple.make(mf_dom["dom02"][1])
+        ift.DomainTuple.make(mf_dom["dom02"][1]), ift.full(mf_dom["dom01"], 1.)
     )
     # SwitchSpacesOperator is equivalent to LinearEinsum with "ij->ji"
     mle_ift = _SwitchSpacesOperator(outer_i.target, 1) @ outer_i @ \
@@ -136,12 +136,12 @@ def test_multi_linear_einsum_outer(space1, space2, dtype):
         (outer_j @ ift.FieldAdapter(mf_dom["dom03"], "dom03"))
 
     for _ in range(n_invocations):
-        rl = ift.Linearization.make_var(ift.from_random("normal", mle.domain, dtype=dtype))
+        rl = ift.Linearization.make_var(ift.from_random(mle.domain, "normal", dtype=dtype))
         mle_rl, mle_ift_rl = mle(rl), mle_ift(rl)
         assert_allclose(mle_rl.val.val, mle_ift_rl.val.val)
         assert_allclose(mle_rl.jac(rl.val).val, mle_ift_rl.jac(rl.val).val)
 
-        rj_adj = ift.from_random("normal", mle_rl.jac.target, dtype=dtype)
+        rj_adj = ift.from_random(mle_rl.jac.target, "normal", dtype=dtype)
         mle_j_val = mle_rl.jac.adjoint(rj_adj).val
         mle_ift_j_val = mle_ift_rl.jac.adjoint(rj_adj).val
         for k in mle_ift.domain.keys():
diff --git a/test/test_operators/test_fft_operator.py b/test/test_operators/test_fft_operator.py
index 9e7263719096fcfd828dd42d6560cc20d1538d1b..54c0a8de4fe7711c923d1411fa01489a08eaa161 100644
--- a/test/test_operators/test_fft_operator.py
+++ b/test/test_operators/test_fft_operator.py
@@ -44,16 +44,14 @@ def test_fft1D(d, dtype, op):
     b = ift.RGSpace(dim1, distances=1./(dim1*d), harmonic=True)
 
     fft = op(domain=a, target=b)
-    inp = ift.Field.from_random(
-        domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
+    inp = ift.Field.from_random(domain=a, random_type='normal', dtype=dtype, std=7, mean=3)
     out = fft.inverse_times(fft.times(inp))
     assert_allclose(inp.val, out.val, rtol=tol, atol=tol)
 
     a, b = b, a
 
     fft = ift.FFTOperator(domain=a, target=b)
-    inp = ift.Field.from_random(
-        domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
+    inp = ift.Field.from_random(domain=a, random_type='normal', dtype=dtype, std=7, mean=3)
     out = fft.inverse_times(fft.times(inp))
     assert_allclose(inp.val, out.val, rtol=tol, atol=tol)
 
@@ -72,16 +70,14 @@ def test_fft2D(dim1, dim2, d1, d2, dtype, op, nthreads):
         [dim1, dim2], distances=[1./(dim1*d1), 1./(dim2*d2)], harmonic=True)
 
     fft = op(domain=a, target=b)
-    inp = ift.Field.from_random(
-        domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
+    inp = ift.Field.from_random(domain=a, random_type='normal', dtype=dtype, std=7, mean=3)
     out = fft.inverse_times(fft.times(inp))
     assert_allclose(inp.val, out.val, rtol=tol, atol=tol)
 
     a, b = b, a
 
     fft = ift.FFTOperator(domain=a, target=b)
-    inp = ift.Field.from_random(
-        domain=a, random_type='normal', std=7, mean=3, dtype=dtype)
+    inp = ift.Field.from_random(domain=a, random_type='normal', dtype=dtype, std=7, mean=3)
     out = fft.inverse_times(fft.times(inp))
     assert_allclose(inp.val, out.val, rtol=tol, atol=tol)
     ift.fft.set_nthreads(1)
@@ -96,8 +92,7 @@ def test_composed_fft(index, dtype, op):
                 ift.RGSpace((5, 6))]
     fft = op(domain=a, space=index)
 
-    inp = ift.Field.from_random(
-        domain=(a1, a2, a3), random_type='normal', std=7, mean=3, dtype=dtype)
+    inp = ift.Field.from_random(domain=(a1, a2, a3), random_type='normal', dtype=dtype, std=7, mean=3)
     out = fft.inverse_times(fft.times(inp))
     assert_allclose(inp.val, out.val, rtol=tol, atol=tol)
 
@@ -111,8 +106,7 @@ def test_normalisation(space, dtype, op):
     tol = 10*_get_rtol(dtype)
     cospace = space.get_default_codomain()
     fft = op(space, cospace)
-    inp = ift.Field.from_random(
-        domain=space, random_type='normal', std=1, mean=2, dtype=dtype)
+    inp = ift.Field.from_random(domain=space, random_type='normal', dtype=dtype, std=1, mean=2)
     out = fft.times(inp)
     fft2 = op(cospace, space)
     out2 = fft2.inverse_times(inp)
diff --git a/test/test_operators/test_harmonic_transform_operator.py b/test/test_operators/test_harmonic_transform_operator.py
index d9b7a44660b28be66e6ecaa40d7605d0b9d4edbc..71a597d9a99680c3ee238ed62110460c74a5eca4 100644
--- a/test/test_operators/test_harmonic_transform_operator.py
+++ b/test/test_operators/test_harmonic_transform_operator.py
@@ -41,8 +41,7 @@ def test_dotsht(lm, tp):
     a = ift.LMSpace(lmax=lm)
     b = ift.GLSpace(nlat=lm + 1)
     fft = ift.HarmonicTransformOperator(domain=a, target=b)
-    inp = ift.Field.from_random(
-        domain=a, random_type='normal', std=1, mean=0, dtype=tp)
+    inp = ift.Field.from_random(domain=a, random_type='normal', dtype=tp, std=1, mean=0)
     out = fft.times(inp)
     v1 = np.sqrt(out.s_vdot(out))
     v2 = np.sqrt(inp.s_vdot(fft.adjoint_times(out)))
@@ -54,8 +53,7 @@ def test_dotsht2(lm, tp):
     a = ift.LMSpace(lmax=lm)
     b = ift.HPSpace(nside=lm//2)
     fft = ift.HarmonicTransformOperator(domain=a, target=b)
-    inp = ift.Field.from_random(
-        domain=a, random_type='normal', std=1, mean=0, dtype=tp)
+    inp = ift.Field.from_random(domain=a, random_type='normal', dtype=tp, std=1, mean=0)
     out = fft.times(inp)
     v1 = np.sqrt(out.s_vdot(out))
     v2 = np.sqrt(inp.s_vdot(fft.adjoint_times(out)))
@@ -67,8 +65,7 @@ def test_normalisation(space, tp):
     tol = 10*_get_rtol(tp)
     cospace = space.get_default_codomain()
     fft = ift.HarmonicTransformOperator(space, cospace)
-    inp = ift.Field.from_random(
-        domain=space, random_type='normal', std=1, mean=2, dtype=tp)
+    inp = ift.Field.from_random(domain=space, random_type='normal', dtype=tp, std=1, mean=2)
     out = fft.times(inp)
     zero_idx = tuple([0]*len(space.shape))
     assert_allclose(
diff --git a/test/test_operators/test_interpolated.py b/test/test_operators/test_interpolated.py
index 985f69c8ef15ea2310fb860a7a827b285dbbcc0f..b7d3d47e6c9357075ac160cac4ef322aa300dca9 100644
--- a/test/test_operators/test_interpolated.py
+++ b/test/test_operators/test_interpolated.py
@@ -33,7 +33,7 @@ seed = list2fixture([4, 78, 23])
 
 
 def testInterpolationAccuracy(space, seed):
-    pos = ift.from_random('normal', space)
+    pos = ift.from_random(space, 'normal')
     alpha = 1.5
     qs = [0.73, pos.ptw("exp").val]
     for q in qs:
diff --git a/test/test_operators/test_jacobian.py b/test/test_operators/test_jacobian.py
index c6faa07113ce7f1792b0cd05665a31dda5bbf43d..4369560bb21a01869daed1120199937fd3fc092f 100644
--- a/test/test_operators/test_jacobian.py
+++ b/test/test_operators/test_jacobian.py
@@ -39,7 +39,7 @@ ntries = 10
 
 def testBasics(space, seed):
     with ift.random.Context(seed):
-        s = ift.from_random('normal', space)
+        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)
@@ -55,42 +55,42 @@ def testBinary(type1, type2, space, seed):
         select_s1 = ift.ducktape(None, dom1, "s1")
         select_s2 = ift.ducktape(None, dom2, "s2")
         model = select_s1*select_s2
-        pos = ift.from_random("normal", dom)
+        pos = ift.from_random(dom, "normal")
         ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
         model = select_s1 + select_s2
-        pos = ift.from_random("normal", dom)
+        pos = ift.from_random(dom, "normal")
         ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
         model = select_s1.scale(3.)
-        pos = ift.from_random("normal", dom1)
+        pos = ift.from_random(dom1, "normal")
         ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
         model = ift.ScalingOperator(space, 2.456)(select_s1*select_s2)
-        pos = ift.from_random("normal", dom)
+        pos = ift.from_random(dom, "normal")
         ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
         model = (2.456*(select_s1*select_s2)).ptw("sigmoid")
-        pos = ift.from_random("normal", dom)
+        pos = ift.from_random(dom, "normal")
         ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
-        pos = ift.from_random("normal", dom)
-        model = ift.OuterProduct(pos['s1'], ift.makeDomain(space))
+        pos = ift.from_random(dom, "normal")
+        model = ift.OuterProduct(ift.makeDomain(space), pos['s1'])
         ift.extra.check_jacobian_consistency(model, pos['s2'], ntries=ntries)
         model = select_s1**2
-        pos = ift.from_random("normal", dom1)
+        pos = ift.from_random(dom1, "normal")
         ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
         model = select_s1.clip(-1, 1)
-        pos = ift.from_random("normal", dom1)
+        pos = ift.from_random(dom1, "normal")
         ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
-        f = ift.from_random("normal", space)
+        f = ift.from_random(space, "normal")
         model = select_s1.clip(f-0.1, f+1.)
-        pos = ift.from_random("normal", dom1)
+        pos = ift.from_random(dom1, "normal")
         ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
         if isinstance(space, ift.RGSpace):
             model = ift.FFTOperator(space)(select_s1*select_s2)
-            pos = ift.from_random("normal", dom)
+            pos = ift.from_random(dom, "normal")
             ift.extra.check_jacobian_consistency(model, pos, ntries=ntries)
 
 
 def testSpecialDistributionOps(space, seed):
     with ift.random.Context(seed):
-        pos = ift.from_random('normal', space)
+        pos = ift.from_random(space, 'normal')
         alpha = 1.5
         q = 0.73
         model = ift.InverseGammaOperator(space, alpha, q)
@@ -102,8 +102,8 @@ def testSpecialDistributionOps(space, seed):
 @pmp('neg', [True, False])
 def testAdder(space, seed, neg):
     with ift.random.Context(seed):
-        f = ift.from_random('normal', space)
-        f1 = ift.from_random('normal', space)
+        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)
         op = ift.Adder(f1.val.ravel()[0], neg=neg, domain=space)
@@ -124,7 +124,7 @@ def testDynamicModel(target, causal, minimum_phase, seed):
                'causal': causal,
                'minimum_phase': minimum_phase}
         model, _ = ift.dynamic_operator(**dct)
-        pos = ift.from_random('normal', model.domain)
+        pos = ift.from_random(model.domain, 'normal')
         ift.extra.check_jacobian_consistency(model, pos, tol=1e-7, ntries=ntries)
         if len(target.shape) > 1:
             dct = {
@@ -143,7 +143,7 @@ def testDynamicModel(target, causal, minimum_phase, seed):
             dct['sigc'] = 1.
             dct['quant'] = 5
             model, _ = ift.dynamic_lightcone_operator(**dct)
-            pos = ift.from_random('normal', model.domain)
+            pos = ift.from_random(model.domain, 'normal')
             ift.extra.check_jacobian_consistency(model, pos, tol=1e-7, ntries=ntries)
 
 
@@ -159,5 +159,5 @@ def testNormalization(h_space, specialbinbounds, logarithmic, nbin):
         binbounds = ift.PowerSpace.useful_binbounds(h_space, logarithmic, nbin)
     dom = ift.PowerSpace(h_space, binbounds)
     op = ift.library.correlated_fields._Normalization(dom)
-    pos = 0.1*ift.from_random('normal', op.domain)
+    pos = 0.1 * ift.from_random(op.domain, 'normal')
     ift.extra.check_jacobian_consistency(op, pos, ntries=10)
diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py
index 70a22845d2648c924b0b1f48dfa0d1f16ff38afb..c1e0511ca15fedae41dbff3e5ef829bfcae63236 100644
--- a/test/test_operators/test_nft.py
+++ b/test/test_operators/test_nft.py
@@ -72,7 +72,7 @@ def test_cartesian():
     GM = ift.GridderMaker(dom, uv=uv)
     op = GM.getFull().adjoint
 
-    fld = ift.from_random('normal', dom)
+    fld = ift.from_random(dom, 'normal')
     arr = fld.val
 
     fld2 = ift.makeField(dom, np.roll(arr, (nx//2, ny//2), axis=(0, 1)))
diff --git a/test/test_operators/test_partial_multifield_insert.py b/test/test_operators/test_partial_multifield_insert.py
index 2fc37d46350850008a414104c9183f7e0bcf77c5..dd07b085fa5edd73d58ebb856c9c6d48375ef1bd 100644
--- a/test/test_operators/test_partial_multifield_insert.py
+++ b/test/test_operators/test_partial_multifield_insert.py
@@ -38,7 +38,7 @@ def test_part_mf_insert():
     a = op1 + op2 + op3
     b = op4 + op5
     op = a.partial_insert(b)
-    fld = ift.from_random('normal', op.domain)
+    fld = ift.from_random(op.domain, 'normal')
     ift.extra.check_jacobian_consistency(op, fld, ntries=ntries)
     assert_(op.domain is ift.MultiDomain.union(
         [op1.domain, op2.domain, op4.domain, op5.domain]))
diff --git a/test/test_operators/test_regridding.py b/test/test_operators/test_regridding.py
index 444efdf7092c3ef957de46cc14edf619e94b66a3..1ad8903b6e106254487d82803dad2c25fb7ad815 100644
--- a/test/test_operators/test_regridding.py
+++ b/test/test_operators/test_regridding.py
@@ -30,5 +30,5 @@ s = list2fixture([
 
 def test_value(s):
     Regrid = ift.RegriddingOperator(s, s.shape)
-    f = ift.from_random('normal', Regrid.domain)
+    f = ift.from_random(Regrid.domain, 'normal')
     assert_allclose(f.val, Regrid(f).val)
diff --git a/test/test_operators/test_representation.py b/test/test_operators/test_representation.py
index e791ebb5eec2004ad874fc7602b5340204418403..45c402126be077a88bbc2656a5d79a9138fd94ba 100644
--- a/test/test_operators/test_representation.py
+++ b/test/test_operators/test_representation.py
@@ -53,8 +53,8 @@ def testLOSResponse(sp, dtype):
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 def testOperatorCombinations(sp, dtype):
-    a = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
-    b = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
+    a = ift.DiagonalOperator(ift.Field.from_random(sp, "normal", dtype=dtype))
+    b = ift.DiagonalOperator(ift.Field.from_random(sp, "normal", dtype=dtype))
     _check_repr(ift.SandwichOperator.make(a, b))
     _check_repr(a(b))
     _check_repr(a+b)
@@ -73,7 +73,7 @@ def testLinearInterpolator():
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
 def testOperatorAdaptor(sp, dtype):
-    op = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
+    op = ift.DiagonalOperator(ift.Field.from_random(sp, "normal", dtype=dtype))
     _check_repr(op.adjoint)
     _check_repr(op.inverse)
     _check_repr(op.inverse.adjoint)
@@ -137,7 +137,7 @@ def testHarmonic(sp, dtype):
 @pmp('sp', _p_spaces)
 def testMask(sp, dtype):
     # Create mask
-    f = ift.from_random('normal', sp).val
+    f = ift.from_random(sp, 'normal').val
     mask = np.zeros_like(f)
     mask[f > 0] = 1
     mask = ift.Field.from_raw(sp, mask)
@@ -147,7 +147,7 @@ def testMask(sp, dtype):
 
 @pmp('sp', _h_spaces + _p_spaces)
 def testDiagonal(sp, dtype):
-    op = ift.DiagonalOperator(ift.Field.from_random("normal", sp, dtype=dtype))
+    op = ift.DiagonalOperator(ift.Field.from_random(sp, "normal", dtype=dtype))
     _check_repr(op)
 
 
@@ -199,8 +199,8 @@ def testRegridding(args):
     ift.DomainTuple.make(ift.RGSpace((10, 12), distances=(0.1, 1.)),)
 ])
 def testOuter(fdomain, domain):
-    f = ift.from_random('normal', fdomain)
-    _check_repr(ift.OuterProduct(f, domain))
+    f = ift.from_random(fdomain, 'normal')
+    _check_repr(ift.OuterProduct(domain, f))
 
 
 @pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
diff --git a/test/test_operators/test_selection_operators.py b/test/test_operators/test_selection_operators.py
index eda7e8f75cea2b65a045c2fed02d116598e2cdc3..0e6fa555e3f730033beb74e27d2e1f344a0f2af1 100644
--- a/test/test_operators/test_selection_operators.py
+++ b/test/test_operators/test_selection_operators.py
@@ -52,7 +52,7 @@ def test_split_operator_first_axes_without_intersections(
     )
     assert consistency_check(split) is None
 
-    r = ift.from_random("normal", dom)
+    r = ift.from_random(dom, "normal")
     split_r = split(r)
     # This relies on the keys of the target domain either being in the order of
     # insertion or being alphabetically sorted
@@ -81,7 +81,7 @@ def test_split_operator_first_axes_with_intersections(
     print(split_idx)
     assert consistency_check(split) is None
 
-    r = ift.from_random("normal", dom)
+    r = ift.from_random(dom, "normal")
     split_r = split(r)
     # This relies on the keys of the target domain either being in the order of
     # insertion or being alphabetically sorted
diff --git a/test/test_operators/test_smoothing_operator.py b/test/test_operators/test_smoothing_operator.py
index e40a77ef83272bef35ade38e356ee435876149c5..2ce2e6df30bdecc5270562d59d22aa3152e88ca0 100644
--- a/test/test_operators/test_smoothing_operator.py
+++ b/test/test_operators/test_smoothing_operator.py
@@ -45,8 +45,8 @@ def test_property(space, sigma):
 
 def test_adjoint_times(space, sigma):
     op = ift.HarmonicSmoothingOperator(space, sigma=sigma)
-    rand1 = ift.Field.from_random('normal', domain=space)
-    rand2 = ift.Field.from_random('normal', domain=space)
+    rand1 = ift.Field.from_random(domain=space, random_type='normal')
+    rand2 = ift.Field.from_random(domain=space, random_type='normal')
     tt1 = rand1.s_vdot(op.times(rand2))
     tt2 = rand2.s_vdot(op.adjoint_times(rand1))
     assert_allclose(tt1, tt2)
@@ -67,8 +67,7 @@ def test_smooth_regular1(sz, d, sigma, tp):
     tol = _get_rtol(tp)
     sp = ift.RGSpace(sz, distances=d)
     smo = ift.HarmonicSmoothingOperator(sp, sigma=sigma)
-    inp = ift.Field.from_random(
-        domain=sp, random_type='normal', std=1, mean=4, dtype=tp)
+    inp = ift.Field.from_random(domain=sp, random_type='normal', dtype=tp, std=1, mean=4)
     out = smo(inp)
     assert_allclose(inp.s_sum(), out.s_sum(), rtol=tol, atol=tol)
 
@@ -81,7 +80,6 @@ def test_smooth_regular2(sz1, sz2, d1, d2, sigma, tp):
     tol = _get_rtol(tp)
     sp = ift.RGSpace([sz1, sz2], distances=[d1, d2])
     smo = ift.HarmonicSmoothingOperator(sp, sigma=sigma)
-    inp = ift.Field.from_random(
-        domain=sp, random_type='normal', std=1, mean=4, dtype=tp)
+    inp = ift.Field.from_random(domain=sp, random_type='normal', dtype=tp, std=1, mean=4)
     out = smo(inp)
     assert_allclose(inp.s_sum(), out.s_sum(), rtol=tol, atol=tol)
diff --git a/test/test_operators/test_value_inserter.py b/test/test_operators/test_value_inserter.py
index 11247989abb32db8e0e7b678134d75c911ea84e4..5d15e27f1f2c72c7ca183378dac3d554a0fd31af 100644
--- a/test/test_operators/test_value_inserter.py
+++ b/test/test_operators/test_value_inserter.py
@@ -35,7 +35,7 @@ def test_value_inserter(sp, seed):
     with ift.random.Context(seed):
         ind = tuple([int(ift.random.current_rng().integers(0, ss - 1)) for ss in sp.shape])
         op = ift.ValueInserter(sp, ind)
-        f = ift.from_random('normal', op.domain)
+        f = ift.from_random(op.domain, 'normal')
     inp = f.val
     ret = op(f).val
     assert_(ret[ind] == inp)
diff --git a/test/test_plot.py b/test/test_plot.py
index f85b70aaae3e3d79506b88c2b68b195e0eef873a..d554f8850c6658b65661909165452a5cdc932b4d 100644
--- a/test/test_plot.py
+++ b/test/test_plot.py
@@ -35,11 +35,11 @@ def test_plots():
 
     fft = ift.FFTOperator(rg_space2)
 
-    field_rg1_1 = ift.from_random('normal', rg_space1)
-    field_rg1_2 = ift.from_random('normal', rg_space1)
-    field_rg2 = ift.from_random('normal', rg_space2)
-    field_hp = ift.from_random('normal', hp_space)
-    field_gl = ift.from_random('normal', gl_space)
+    field_rg1_1 = ift.from_random(rg_space1, 'normal')
+    field_rg1_2 = ift.from_random(rg_space1, 'normal')
+    field_rg2 = ift.from_random(rg_space2, 'normal')
+    field_hp = ift.from_random(hp_space, 'normal')
+    field_gl = ift.from_random(gl_space, 'normal')
     field_ps = ift.power_analyze(fft.times(field_rg2))
 
     plot = ift.Plot()
@@ -69,7 +69,7 @@ def test_mf_plot():
     d1 = ift.DomainTuple.make([x_space, f_space])
     d2 = ift.DomainTuple.make([f_space, x_space])
 
-    f1 = ift.from_random('normal', d1)
+    f1 = ift.from_random(d1, 'normal')
     f2 = ift.makeField(d2, np.moveaxis(f1.val, -1, 0))
 
     plot = ift.Plot()
diff --git a/test/test_sugar.py b/test/test_sugar.py
index 174a3f0d07f3cf32f1cedcf5fc168ff2c8b99b12..850268f6f9c98fe450b0453ac8009e3390cd60bb 100644
--- a/test/test_sugar.py
+++ b/test/test_sugar.py
@@ -55,6 +55,6 @@ def test_exec_time():
 def test_calc_pos():
     dom = ift.RGSpace(12, harmonic=True)
     op = ift.HarmonicTransformOperator(dom).ptw("exp")
-    fld = op(0.1*ift.from_random('normal', op.domain))
+    fld = op(0.1 * ift.from_random(op.domain, 'normal'))
     pos = ift.calculate_position(op, fld)
     ift.extra.assert_allclose(op(pos), fld, 1e-1, 1e-1)