diff --git a/nifty_interfaces.py b/nifty_interfaces.py
new file mode 100644
index 0000000000000000000000000000000000000000..463a7f4eb3d629298973626b82990f35e92d8ac3
--- /dev/null
+++ b/nifty_interfaces.py
@@ -0,0 +1,68 @@
+import nifty8 as ift
+import numpy as np
+
+# Domains / spaces
+sp = ift.RGSpace(shape=(2, 3),
+                 distances=(1., 1.5))
+print(sp.shape)
+print(sp.distances)
+print(sp.total_volume)
+print(sp.scalar_dvol)
+print()
+
+sp1 = ift.UnstructuredDomain(4)
+print(sp1)
+print()
+
+# Fields
+d_arr = np.array([12.1, 4.3, 21., 1110.])
+d = ift.makeField(sp1, d_arr)
+print(d)
+s = ift.from_random(sp)  # draws std normal samples
+print(s)
+print(s.integrate())
+print(s.sum())
+print()
+print(s.domain)
+print(s.val)
+#print(s.val_rw())
+
+
+# Product domains
+dom = ift.makeDomain((sp, sp1))
+print(dom)
+print(dom.shape)
+print(dom.size)
+
+# MultiDomain: dictionaries of Domains
+dom = ift.makeDomain({"key0": sp, "key1": sp1})
+print(dom)
+print(dom.size)
+
+mfld = ift.from_random(dom)
+mfld = mfld.sin()
+print(mfld.domain)
+print(mfld.val)
+
+
+
+print()
+print()
+print()
+print()
+# Operators
+sp = ift.RGSpace((2, 2), (1, 1))
+mask = ift.makeField(sp,
+         np.array([[False, False],
+                   [False, True]])
+         )
+op = ift.MaskOperator(mask)
+print(op)
+print(op.domain)
+print(op.target)
+
+inp = ift.from_random(op.domain)
+print("Input")
+print(inp)
+print("Output")
+print(op(inp))
diff --git a/nifty_intro.py b/nifty_intro.py
new file mode 100644
index 0000000000000000000000000000000000000000..50e1130edb491ec0e76744e02a4e11a8da9d686b
--- /dev/null
+++ b/nifty_intro.py
@@ -0,0 +1,75 @@
+import nifty8 as ift
+from matplotlib import pyplot as plt
+import numpy as np
+
+
+plt.ion()
+position_space = ift.RGSpace([128,128])
+harmonic_space = position_space.get_default_codomain()
+
+HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
+
+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_field = pd(a)
+A = ift.makeOp(A_field)
+xi = ift.FieldAdapter(harmonic_space,'xi')
+
+args = {
+    'offset_mean': 0,
+    'offset_std': (1e-3, 1e-6),
+
+    # Amplitude of field fluctuations
+    'fluctuations': (1., 0.8),  # 1.0, 1e-2
+
+    # Exponent of power law power spectrum component
+    'loglogavgslope': (-3., 1),  # -6.0, 1
+
+    # Amplitude of integrated Wiener process power spectrum component
+    'flexibility': (2, 1.),  # 1.0, 0.5
+
+    # How ragged the integrated Wiener process component is
+    'asperity': (0.5, 0.4)  # 0.1, 0.5
+}
+
+correlated_field = ift.SimpleCorrelatedField(position_space, **args)
+
+# GP = HT @ A @ xi
+
+xi_true = ift.from_random(correlated_field.domain)
+
+s = correlated_field(xi_true)
+
+sky = ift.exp(correlated_field)
+
+sky_true = sky(xi_true)
+
+data_raw = np.random.poisson(sky_true.val)
+
+data = ift.makeField(sky.target, data_raw)
+
+log_likelihood = ift.PoissonianEnergy(data) @ sky
+
+ic_sampling = ift.AbsDeltaEnergyController(name="Sampling (linear)",
+                                               deltaE=0.05, iteration_limit=100)
+
+ic_sampling_nl = ift.AbsDeltaEnergyController(name="Sampling (non-linear)",
+                                               deltaE=0.05, iteration_limit=10)
+
+
+ic_newton = ift.AbsDeltaEnergyController(name='Newton', deltaE=0.5,
+                                             convergence_level=2, iteration_limit=5)
+
+minimizer = ift.NewtonCG(ic_newton)
+minimizer_sampling = (lambda iiter: None if iiter < 3 else ift.NewtonCG(ic_sampling_nl))
+
+n_iterations = 5
+n_samples = 5
+
+samples = ift.optimize_kl(log_likelihood, n_iterations, n_samples, minimizer,
+                          ic_sampling, minimizer_sampling)