diff --git a/demos/bernoulli_demo.py b/demos/bernoulli_demo.py
index dbcd16cfd6334c4e691b9fd6cb103281ca9d4163..5ae269ead480611024001231f488d53aa1d00af3 100644
--- a/demos/bernoulli_demo.py
+++ b/demos/bernoulli_demo.py
@@ -47,25 +47,20 @@ if __name__ == '__main__':
     def sqrtpspec(k):
         return 1. / (20. + k**2)
 
-    p_space = ift.PowerSpace(harmonic_space)
-    pd = ift.PowerDistributor(harmonic_space, p_space)
-    a = ift.PS_field(p_space, sqrtpspec)
-    A = pd(a)
+    A = ift.create_power_operator(harmonic_space, sqrtpspec)
 
     # Set up a sky model
-    sky = HT(ift.makeOp(A)).positive_tanh()
+    sky = ift.positive_tanh(HT(A))
 
     GR = ift.GeometryRemover(position_space)
     # Set up instrumental response
     R = GR
 
     # Generate mock data
-    d_space = R.target[0]
     p = R(sky)
     mock_position = ift.from_random('normal', harmonic_space)
-    pp = p(mock_position)
-    data = np.random.binomial(1, pp.to_global_data().astype(np.float64))
-    data = ift.Field.from_global_data(d_space, data)
+    data = np.random.binomial(1, p(mock_position).local_data.astype(np.float64))
+    data = ift.Field.from_local_data(R.target, data)
 
     # Compute likelihood and Hamiltonian
     position = ift.from_random('normal', harmonic_space)
@@ -79,7 +74,7 @@ if __name__ == '__main__':
     # Minimize the Hamiltonian
     H = ift.Hamiltonian(likelihood, ic_sampling)
     H = ift.EnergyAdapter(position, H, ic_cg)
-    # minimizer = ift.SteepestDescent(ic_newton)
+    # minimizer = ift.L_BFGS(ic_newton)
     H, convergence = minimizer(H)
 
     reconstruction = sky(H.position)
diff --git a/test/test_models/test_model_gradients.py b/test/test_models/test_model_gradients.py
index 236ba145c8a6d3b11e945b76e41abe52885100fe..47befaa2d55060e4c0b3313dcad936610c907323 100644
--- a/test/test_models/test_model_gradients.py
+++ b/test/test_models/test_model_gradients.py
@@ -69,7 +69,7 @@ class Model_Tests(unittest.TestCase):
         pos = ift.from_random("normal", dom)
         ift.extra.check_value_gradient_consistency(model, pos)
         model = ift.FieldAdapter(dom, "s1").scale(3.)
-        pos = ift.from_random("normal", dom)
+        pos = ift.from_random("normal", dom1)
         ift.extra.check_value_gradient_consistency(model, pos)
         model = ift.ScalingOperator(2.456, space)(
             ift.FieldAdapter(dom, "s1")*ift.FieldAdapter(dom, "s2"))