Commit dd8e0f42 by Martin Reinecke

### Merge branch 'move_getting_started_to_docs' into 'NIFTy_8'

Add getting started to docs

See merge request !666
parents 88258c6e b530bdb1
Pipeline #106640 passed with stages
in 34 minutes and 30 seconds
 # never store the git version file git_version.py docs/source/user/getting_started_0.rst # custom *.txt setup.cfg ... ...
 ... ... @@ -8,7 +8,7 @@ } }, "source": [ "# A NIFTy demonstration" "# Code example: Wiener filter" ] }, { ... ... @@ -19,7 +19,7 @@ } }, "source": [ "## IFT: Big Picture\n", "## Introduction\n", "IFT starting point:\n", "\n", "$$d = Rs+n$$\n", ... ... @@ -28,9 +28,6 @@ "\n", "IFT aims at **inverting** the above uninvertible problem in the **best possible way** using Bayesian statistics.\n", "\n", "\n", "## NIFTy\n", "\n", "NIFTy (Numerical Information Field Theory) is a Python framework in which IFT problems can be tackled easily.\n", "\n", "Main Interfaces:\n", ... ... @@ -48,7 +45,7 @@ } }, "source": [ "## Wiener Filter: Formulae\n", "## Wiener filter on one-dimensional fields\n", "\n", "### Assumptions\n", "\n", ... ... @@ -61,11 +58,9 @@ "$$\\mathcal P (s|d) \\propto P(s,d) = \\mathcal G(d-Rs,N) \\,\\mathcal G(s,S) \\propto \\mathcal G (s-m,D)$$\n", "\n", "where\n", "\\begin{align}\n", "m &= Dj \\\\\n", "D^{-1}&= (S^{-1} +R^\\dagger N^{-1} R )\\\\\n", "j &= R^\\dagger N^{-1} d\n", "\\end{align}\n", "$$m = Dj$$\n", "with\n", "$$D = (S^{-1} +R^\\dagger N^{-1} R)^{-1} , \\quad j = R^\\dagger N^{-1} d.$$\n", "\n", "Let us implement this in NIFTy!" ] ... ... @@ -78,7 +73,7 @@ } }, "source": [ "## Wiener Filter: Example\n", "### In NIFTy\n", "\n", "- We assume statistical homogeneity and isotropy. Therefore the signal covariance $S$ is diagonal in harmonic space, and is described by a one-dimensional power spectrum, assumed here as $$P(k) = P_0\\,\\left(1+\\left(\\frac{k}{k_0}\\right)^2\\right)^{-\\gamma /2},$$\n", "with $P_0 = 0.2, k_0 = 5, \\gamma = 4$.\n", ... ... @@ -114,7 +109,7 @@ } }, "source": [ "## Wiener Filter: Implementation" "### Implementation" ] }, { ... ... @@ -125,7 +120,7 @@ } }, "source": [ "### Import Modules" "#### Import Modules" ] }, { ... ... @@ -142,7 +137,8 @@ "import numpy as np\n", "import nifty8 as ift\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" "plt.rcParams['figure.dpi'] = 100\n", "plt.style.use(\"seaborn-notebook\")" ] }, { ... ... @@ -153,7 +149,7 @@ } }, "source": [ "### Implement Propagator" "#### Implement Propagator" ] }, { ... ... @@ -182,7 +178,7 @@ } }, "source": [ "### Conjugate Gradient Preconditioning\n", "#### Conjugate Gradient Preconditioning\n", "\n", "- $D$ is defined via:\n", "$$D^{-1} = \\mathcal S_h^{-1} + R^\\dagger N^{-1} R.$$\n", ... ... @@ -213,7 +209,7 @@ } }, "source": [ "### Generate Mock data\n", "#### Generate Mock data\n", "\n", "- Generate a field $s$ and $n$ with given covariances.\n", "- Calculate $d$." ... ... @@ -257,7 +253,7 @@ } }, "source": [ "### Run Wiener Filter" "#### Run Wiener Filter" ] }, { ... ... @@ -281,7 +277,7 @@ } }, "source": [ "### Signal Reconstruction" "#### Results" ] }, { ... ... @@ -299,10 +295,9 @@ "m_data = HT(m).val\n", "d_data = d.val\n", "\n", "plt.figure(figsize=(15,10))\n", "plt.plot(s_data, 'r', label=\"Signal\", linewidth=3)\n", "plt.plot(s_data, 'r', label=\"Signal\", linewidth=2)\n", "plt.plot(d_data, 'k.', label=\"Data\")\n", "plt.plot(m_data, 'k', label=\"Reconstruction\",linewidth=3)\n", "plt.plot(m_data, 'k', label=\"Reconstruction\",linewidth=2)\n", "plt.title(\"Reconstruction\")\n", "plt.legend()\n", "plt.show()" ... ... @@ -318,10 +313,9 @@ }, "outputs": [], "source": [ "plt.figure(figsize=(15,10))\n", "plt.plot(s_data - s_data, 'r', label=\"Signal\", linewidth=3)\n", "plt.plot(s_data - s_data, 'r', label=\"Signal\", linewidth=2)\n", "plt.plot(d_data - s_data, 'k.', label=\"Data\")\n", "plt.plot(m_data - s_data, 'k', label=\"Reconstruction\",linewidth=3)\n", "plt.plot(m_data - s_data, 'k', label=\"Reconstruction\",linewidth=2)\n", "plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5)\n", "plt.title(\"Residuals\")\n", "plt.legend()\n", ... ... @@ -336,7 +330,7 @@ } }, "source": [ "### Power Spectrum" "#### Power Spectrum" ] }, { ... ... @@ -351,7 +345,6 @@ "source": [ "s_power_data = ift.power_analyze(sh).val\n", "m_power_data = ift.power_analyze(m).val\n", "plt.figure(figsize=(15,10))\n", "plt.loglog()\n", "plt.xlim(1, int(N_pixels/2))\n", "ymin = min(m_power_data)\n", ... ... @@ -516,12 +509,11 @@ }, "outputs": [], "source": [ "fig = plt.figure(figsize=(15,10))\n", "plt.axvspan(l, h, facecolor='0.8',alpha=0.5)\n", "plt.fill_between(range(N_pixels), m_data - uncertainty, m_data + uncertainty, facecolor='0.5', alpha=0.5)\n", "plt.plot(s_data, 'r', label=\"Signal\", alpha=1, linewidth=3)\n", "plt.plot(s_data, 'r', label=\"Signal\", alpha=1, linewidth=2)\n", "plt.plot(d_data, 'k.', label=\"Data\")\n", "plt.plot(m_data, 'k', label=\"Reconstruction\", linewidth=3)\n", "plt.plot(m_data, 'k', label=\"Reconstruction\", linewidth=2)\n", "plt.title(\"Reconstruction of incomplete data\")\n", "plt.legend()" ] ... ... @@ -534,7 +526,7 @@ } }, "source": [ "# 2d Example" "## Wiener filter on two-dimensional field" ] }, { ... ... @@ -618,18 +610,18 @@ }, "outputs": [], "source": [ "cm = ['magma', 'inferno', 'plasma', 'viridis'][1]\n", "cmap = ['magma', 'inferno', 'plasma', 'viridis'][1]\n", "\n", "mi = np.min(s_data)\n", "ma = np.max(s_data)\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(15, 7))\n", "fig, axes = plt.subplots(1, 2)\n", "\n", "data = [s_data, d_data]\n", "caption = [\"Signal\", \"Data\"]\n", "\n", "for ax in axes.flat:\n", " im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi,\n", " im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cmap, vmin=mi,\n", " vmax=ma)\n", " ax.set_title(caption.pop(0))\n", "\n", ... ... @@ -651,7 +643,7 @@ "mi = np.min(s_data)\n", "ma = np.max(s_data)\n", "\n", "fig, axes = plt.subplots(3, 2, figsize=(15, 22.5))\n", "fig, axes = plt.subplots(3, 2, figsize=(10, 15))\n", "sample = HT(curv.draw_sample(from_inverse=True)+m).val\n", "post_mean = (m_mean + HT(m)).val\n", "\n", ... ... @@ -659,7 +651,7 @@ "caption = [\"Signal\", \"Reconstruction\", \"Posterior mean\", \"Sample\", \"Residuals\", \"Uncertainty Map\"]\n", "\n", "for ax in axes.flat:\n", " im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi, vmax=ma)\n", " im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cmap, vmin=mi, vmax=ma)\n", " ax.set_title(caption.pop(0))\n", "\n", "fig.subplots_adjust(right=0.8)\n", ... ... @@ -691,31 +683,9 @@ "precise = (np.abs(s_data-m_data) < uncertainty)\n", "print(\"Error within uncertainty map bounds: \" + str(np.sum(precise) * 100 / N_pixels**2) + \"%\")\n", "\n", "plt.figure(figsize=(15,10))\n", "plt.imshow(precise.astype(float), cmap=\"brg\")\n", "plt.colorbar()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Start Coding\n", "## NIFTy Repository + Installation guide\n", "\n", "https://gitlab.mpcdf.mpg.de/ift/NIFTy\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { ... ... @@ -735,7 +705,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.2" "version": "3.9.2" } }, "nbformat": 4, ... ...
 %% Cell type:markdown id: tags: # A NIFTy demonstration # Code example: Wiener filter %% Cell type:markdown id: tags: ## IFT: Big Picture ## Introduction IFT starting point: $$d = Rs+n$$ Typically, $s$ is a continuous field, $d$ a discrete data vector. Particularly, $R$ is not invertible. IFT aims at **inverting** the above uninvertible problem in the **best possible way** using Bayesian statistics. ## NIFTy NIFTy (Numerical Information Field Theory) is a Python framework in which IFT problems can be tackled easily. Main Interfaces: - **Spaces**: Cartesian, 2-Spheres (Healpix, Gauss-Legendre) and their respective harmonic spaces. - **Fields**: Defined on spaces. - **Operators**: Acting on fields. %% Cell type:markdown id: tags: ## Wiener Filter: Formulae ## Wiener filter on one-dimensional fields ### Assumptions - $d=Rs+n$, $R$ linear operator. - $\mathcal P (s) = \mathcal G (s,S)$, $\mathcal P (n) = \mathcal G (n,N)$ where $S, N$ are positive definite matrices. ### Posterior The Posterior is given by: $$\mathcal P (s|d) \propto P(s,d) = \mathcal G(d-Rs,N) \,\mathcal G(s,S) \propto \mathcal G (s-m,D)$$ where \begin{align} m &= Dj \\ D^{-1}&= (S^{-1} +R^\dagger N^{-1} R )\\ j &= R^\dagger N^{-1} d \end{align} $$m = Dj$$ with $$D = (S^{-1} +R^\dagger N^{-1} R)^{-1} , \quad j = R^\dagger N^{-1} d.$$ Let us implement this in NIFTy! %% Cell type:markdown id: tags: ## Wiener Filter: Example ### In NIFTy - We assume statistical homogeneity and isotropy. Therefore the signal covariance $S$ is diagonal in harmonic space, and is described by a one-dimensional power spectrum, assumed here as $$P(k) = P_0\,\left(1+\left(\frac{k}{k_0}\right)^2\right)^{-\gamma /2},$$ with $P_0 = 0.2, k_0 = 5, \gamma = 4$. - $N = 0.2 \cdot \mathbb{1}$. - Number of data points $N_{pix} = 512$. - reconstruction in harmonic space. - Response operator: $$R = FFT_{\text{harmonic} \rightarrow \text{position}}$$ %% Cell type:code id: tags:  python N_pixels = 512 # Number of pixels def pow_spec(k): P0, k0, gamma = [.2, 5, 4] return P0 / ((1. + (k/k0)**2)**(gamma / 2))  %% Cell type:markdown id: tags: ## Wiener Filter: Implementation ### Implementation %% Cell type:markdown id: tags: ### Import Modules #### Import Modules %% Cell type:code id: tags:  python import numpy as np import nifty8 as ift import matplotlib.pyplot as plt %matplotlib inline plt.rcParams['figure.dpi'] = 100 plt.style.use("seaborn-notebook")  %% Cell type:markdown id: tags: ### Implement Propagator #### Implement Propagator %% Cell type:code id: tags:  python def Curvature(R, N, Sh): IC = ift.GradientNormController(iteration_limit=50000, tol_abs_gradnorm=0.1) # WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy # helper methods. return ift.WienerFilterCurvature(R,N,Sh,iteration_controller=IC,iteration_controller_sampling=IC)  %% Cell type:markdown id: tags: ### Conjugate Gradient Preconditioning #### Conjugate Gradient Preconditioning - $D$ is defined via: $$D^{-1} = \mathcal S_h^{-1} + R^\dagger N^{-1} R.$$ In the end, we want to apply $D$ to $j$, i.e. we need the inverse action of $D^{-1}$. This is done numerically (algorithm: *Conjugate Gradient*). %% Cell type:markdown id: tags: ### Generate Mock data #### Generate Mock data - Generate a field $s$ and $n$ with given covariances. - Calculate $d$. %% Cell type:code id: tags:  python s_space = ift.RGSpace(N_pixels) h_space = s_space.get_default_codomain() HT = ift.HarmonicTransformOperator(h_space, target=s_space) # Operators Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec) R = HT # @ ift.create_harmonic_smoothing_operator((h_space,), 0, 0.02) # Fields and data sh = Sh.draw_sample_with_dtype(dtype=np.float64) noiseless_data=R(sh) noise_amplitude = np.sqrt(0.2) N = ift.ScalingOperator(s_space, noise_amplitude**2) n = ift.Field.from_random(domain=s_space, random_type='normal', std=noise_amplitude, mean=0) d = noiseless_data + n j = R.adjoint_times(N.inverse_times(d)) curv = Curvature(R=R, N=N, Sh=Sh) D = curv.inverse  %% Cell type:markdown id: tags: ### Run Wiener Filter #### Run Wiener Filter %% Cell type:code id: tags:  python m = D(j)  %% Cell type:markdown id: tags: ### Signal Reconstruction #### Results %% Cell type:code id: tags:  python # Get signal data and reconstruction data s_data = HT(sh).val m_data = HT(m).val d_data = d.val plt.figure(figsize=(15,10)) plt.plot(s_data, 'r', label="Signal", linewidth=3) plt.plot(s_data, 'r', label="Signal", linewidth=2) plt.plot(d_data, 'k.', label="Data") plt.plot(m_data, 'k', label="Reconstruction",linewidth=3) plt.plot(m_data, 'k', label="Reconstruction",linewidth=2) plt.title("Reconstruction") plt.legend() plt.show()  %% Cell type:code id: tags:  python plt.figure(figsize=(15,10)) plt.plot(s_data - s_data, 'r', label="Signal", linewidth=3) plt.plot(s_data - s_data, 'r', label="Signal", linewidth=2) plt.plot(d_data - s_data, 'k.', label="Data") plt.plot(m_data - s_data, 'k', label="Reconstruction",linewidth=3) plt.plot(m_data - s_data, 'k', label="Reconstruction",linewidth=2) plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5) plt.title("Residuals") plt.legend() plt.show()  %% Cell type:markdown id: tags: ### Power Spectrum #### Power Spectrum %% Cell type:code id: tags:  python s_power_data = ift.power_analyze(sh).val m_power_data = ift.power_analyze(m).val plt.figure(figsize=(15,10)) plt.loglog() plt.xlim(1, int(N_pixels/2)) ymin = min(m_power_data) plt.ylim(ymin, 1) xs = np.arange(1,int(N_pixels/2),.1) plt.plot(xs, pow_spec(xs), label="True Power Spectrum", color='k',alpha=0.5) plt.plot(s_power_data, 'r', label="Signal") plt.plot(m_power_data, 'k', label="Reconstruction") plt.axhline(noise_amplitude**2 / N_pixels, color="k", linestyle='--', label="Noise level", alpha=.5) plt.axhspan(noise_amplitude**2 / N_pixels, ymin, facecolor='0.9', alpha=.5) plt.title("Power Spectrum") plt.legend() plt.show()  %% Cell type:markdown id: tags: ## Wiener Filter on Incomplete Data %% Cell type:code id: tags:  python # Operators Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec) N = ift.ScalingOperator(s_space, noise_amplitude**2) # R is defined below # Fields sh = Sh.draw_sample_with_dtype(dtype=np.float64) s = HT(sh) n = ift.Field.from_random(domain=s_space, random_type='normal', std=noise_amplitude, mean=0)  %% Cell type:markdown id: tags: ### Partially Lose Data %% Cell type:code id: tags:  python l = int(N_pixels * 0.2) h = int(N_pixels * 0.2 * 2) mask = np.full(s_space.shape, 1.) mask[l:h] = 0 mask = ift.Field.from_raw(s_space, mask) R = ift.DiagonalOperator(mask)(HT) n = n.val_rw() n[l:h] = 0 n = ift.Field.from_raw(s_space, n) d = R(sh) + n  %% Cell type:code id: tags:  python curv = Curvature(R=R, N=N, Sh=Sh) D = curv.inverse j = R.adjoint_times(N.inverse_times(d)) m = D(j)  %% Cell type:markdown id: tags: ### Compute Uncertainty %% Cell type:code id: tags:  python m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 200, np.float64)  %% Cell type:markdown id: tags: ### Get data %% Cell type:code id: tags:  python # Get signal data and reconstruction data s_data = s.val m_data = HT(m).val m_var_data = m_var.val uncertainty = np.sqrt(m_var_data) d_data = d.val_rw() # Set lost data to NaN for proper plotting d_data[d_data == 0] = np.nan  %% Cell type:code id: tags:  python fig = plt.figure(figsize=(15,10)) plt.axvspan(l, h, facecolor='0.8',alpha=0.5) plt.fill_between(range(N_pixels), m_data - uncertainty, m_data + uncertainty, facecolor='0.5', alpha=0.5) plt.plot(s_data, 'r', label="Signal", alpha=1, linewidth=3) plt.plot(s_data, 'r', label="Signal", alpha=1, linewidth=2) plt.plot(d_data, 'k.', label="Data") plt.plot(m_data, 'k', label="Reconstruction", linewidth=3) plt.plot(m_data, 'k', label="Reconstruction", linewidth=2) plt.title("Reconstruction of incomplete data") plt.legend()  %% Cell type:markdown id: tags: # 2d Example ## Wiener filter on two-dimensional field %% Cell type:code id: tags:  python N_pixels = 256 # Number of pixels sigma2 = 2. # Noise variance def pow_spec(k): P0, k0, gamma = [.2, 2, 4] return P0 * (1. + (k/k0)**2)**(-gamma/2) s_space = ift.RGSpace([N_pixels, N_pixels])  %% Cell type:code id: tags:  python h_space = s_space.get_default_codomain() HT = ift.HarmonicTransformOperator(h_space,s_space) # Operators Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec) N = ift.ScalingOperator(s_space, sigma2) # Fields and data sh = Sh.draw_sample_with_dtype(dtype=np.float64) n = ift.Field.from_random(domain=s_space, random_type='normal', std=np.sqrt(sigma2), mean=0) # Lose some data l = int(N_pixels * 0.33) h = int(N_pixels * 0.33 * 2) mask = np.full(s_space.shape, 1.) mask[l:h,l:h] = 0. mask = ift.Field.from_raw(s_space, mask) R = ift.DiagonalOperator(mask)(HT) n = n.val_rw() n[l:h, l:h] = 0 n = ift.Field.from_raw(s_space, n) curv = Curvature(R=R, N=N, Sh=Sh) D = curv.inverse d = R(sh) + n j = R.adjoint_times(N.inverse_times(d)) # Run Wiener filter m = D(j) # Uncertainty m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 20, np.float64) # Get data s_data = HT(sh).val m_data = HT(m).val m_var_data = m_var.val d_data = d.val uncertainty = np.sqrt(np.abs(m_var_data))  %% Cell type:code id: tags:  python cm = ['magma', 'inferno', 'plasma', 'viridis'][1] cmap = ['magma', 'inferno', 'plasma', 'viridis'][1] mi = np.min(s_data) ma = np.max(s_data) fig, axes = plt.subplots(1, 2, figsize=(15, 7)) fig, axes = plt.subplots(1, 2) data = [s_data, d_data] caption = ["Signal", "Data"] for ax in axes.flat: im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi, im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cmap, vmin=mi, vmax=ma) ax.set_title(caption.pop(0)) fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax)  %% Cell type:code id: tags:  python mi = np.min(s_data) ma = np.max(s_data) fig, axes = plt.subplots(3, 2, figsize=(15, 22.5)) fig, axes = plt.subplots(3, 2, figsize=(10, 15)) sample = HT(curv.draw_sample(from_inverse=True)+m).val post_mean = (m_mean + HT(m)).val data = [s_data, m_data, post_mean, sample, s_data - m_data, uncertainty] caption = ["Signal", "Reconstruction", "Posterior mean", "Sample", "Residuals", "Uncertainty Map"] for ax in axes.flat: im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi, vmax=ma) im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cmap, vmin=mi, vmax=ma) ax.set_title(caption.pop(0)) fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([.85, 0.15, 0.05, 0.7]) fig.colorbar(im, cax=cbar_ax)  %% Cell type:markdown id: tags: ### Is the uncertainty map reliable? %% Cell type:code id: tags:  python precise = (np.abs(s_data-m_data) < uncertainty) print("Error within uncertainty map bounds: " + str(np.sum(precise) * 100 / N_pixels**2) + "%") plt.figure(figsize=(15,10)) plt.imshow(precise.astype(float), cmap="brg") plt.colorbar()  %% Cell type:markdown id: tags: # Start Coding ## NIFTy Repository + Installation guide https://gitlab.mpcdf.mpg.de/ift/NIFTy %% Cell type:code id: tags:  python  ... ...
 ... ... @@ -31,6 +31,8 @@ "source": [ "import nifty8 as ift\n", "import matplotlib.pyplot as plt\n", "plt.rcParams['figure.dpi'] = 100\n", "plt.style.use(\"seaborn-notebook\")\n", "import numpy as np\n", "ift.random.push_sseq_from_seed(43)\n", "\n", ... ... @@ -49,7 +51,7 @@ " # Plotting preparation is normally handled by nifty8.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", " fig = plt.figure(tight_layout=True, figsize=(10, 3))\n", " if title is not None:\n", " fig.suptitle(title, fontsize=14)\n", " \n", ... ...
</
 %% Cell type:markdown id: tags: # Notebook showcasing the NIFTy 6 Correlated Field model **Skip to Parameter Showcases for the meat/veggies ;)** The field model roughly works like this: f = HT( A * zero_mode * xi ) + offset A is a spectral power field which is constructed from power spectra that hold on subdomains of the target domain. It is scaled by a zero mode operator and then pointwise multiplied by a gaussian excitation field, yielding a representation of the field in harmonic space. It is then transformed into the target real space and a offset added. The power spectra A is constructed of are in turn constructed as the sum of a power law component and an integrated Wiener process whose amplitude and roughness can be set. ## Setup code %% Cell type:code id: tags:  python import nifty8 as ift import matplotlib.pyplot as plt plt.rcParams['figure.dpi'] = 100 plt.style.use("seaborn-notebook") import numpy as np ift.random.push_sseq_from_seed(43) n_pix = 256 x_space = ift.RGSpace(n_pix)  %% Cell type:code id: tags:  python # Plotting routine def plot(fields, spectra, title=None): # Plotting preparation is normally handled by nifty8.Plot # It is done manually here to be able to tweak details # Fields are assumed to have identical domains fig = plt.figure(tight_layout=True, figsize=(12, 3.5)) fig = plt.figure(tight_layout=True, figsize=(10, 3)) if title is not None: fig.suptitle(title, fontsize=14) # Field ax1 = fig.add_subplot(1, 2, 1) ax1.axhline(y=0., color='k', linestyle='--', alpha=0.25) for field in fields: dom = field.domain[0] xcoord = np.arange(dom.shape[0]) * dom.distances[0] ax1.plot(xcoord, field.val) ax1.set_xlim(xcoord[0], xcoord[-1]) ax1.set_ylim(-5., 5.) ax1.set_xlabel('x') ax1.set_ylabel('f(x)') ax1.set_title('Field realizations') # Spectrum ax2 = fig.add_subplot(1, 2, 2) for spectrum in spectra: xcoord = spectrum.domain[0].k_lengths ycoord = spectrum.val_rw() ycoord[0] = ycoord[1] ax2.plot(xcoord, ycoord) ax2.set_ylim(1e-6, 10.) ax2.set_xscale('log') ax2.set_yscale('log') ax2.set_xlabel('k') ax2.set_ylabel('p(k)') ax2.set_title('Power Spectrum') fig.align_labels() plt.show() # Helper: draw main sample main_sample = None def init_model(m_pars, fl_pars): global main_sample cf = ift.CorrelatedFieldMaker(m_pars["prefix"]) cf.set_amplitude_total_offset(m_pars["offset_mean"], m_pars["offset_std"]) cf.add_fluctuations(**fl_pars) field = cf.finalize(prior_info=0) main_sample = ift.from_random(field.domain) print("model domain keys:", field.domain.keys()) # Helper: field and spectrum from parameter dictionaries + plotting def eval_model(m_pars, fl_pars, title=None, samples=None): cf = ift.CorrelatedFieldMaker(m_pars["prefix"]) cf.set_amplitude_total_offset(m_pars["offset_mean"], m_pars["offset_std"]) cf.add_fluctuations(**fl_pars) field = cf.finalize(prior_info=0) spectrum = cf.amplitude if samples is None: samples = [main_sample] field_realizations = [field(s) for s in samples] spectrum_realizations = [spectrum.force(s) for s in samples] plot(field_realizations, spectrum_realizations, title) def gen_samples(key_to_vary): if key_to_vary is None: return [main_sample] dct = main_sample.to_dict() subdom_to_vary = dct.pop(key_to_vary).domain samples = [] for i in range(8): d = dct.copy() d[key_to_vary] = ift.from_random(subdom_to_vary) samples.append(ift.MultiField.from_dict(d)) return samples def vary_parameter(parameter_key, values, samples_vary_in=None): s = gen_samples(samples_vary_in) for v in values: if parameter_key in cf_make_pars.keys(): m_pars = {**cf_make_pars, parameter_key: v} eval_model(m_pars, cf_x_fluct_pars, f"{parameter_key} = {v}", s) else: fl_pars = {**cf_x_fluct_pars, parameter_key: v} eval_model(cf_make_pars, fl_pars, f"{parameter_key} = {v}", s)  %% Cell type:markdown id: tags: ## Before the Action: The Moment-Matched Log-Normal Distribution Many properties of the correlated field are modelled as being lognormally distributed. The distribution models are parametrized via their means and standard-deviations (first and second position in tuple). To get a feeling of how the ratio of the mean and stddev parameters influences the distribution shape, here are a few example histograms: (observe the x-axis!) %% Cell type:code id: tags:  python fig = plt.figure(figsize=(13, 3.5)) mean = 1.0 sigmas = [1.0, 0.5, 0.1] for i in range(3): op = ift.LognormalTransform(mean=mean, sigma=sigmas[i], key='foo', N_copies=0) op_samples = np.array( [op(s).val for s in [ift.from_random(op.domain) for i in range(10000)]]) ax = fig.add_subplot(1, 3, i + 1) ax.hist(op_samples, bins=50) ax.set_title(f"mean = {mean}, sigma = {sigmas[i]}") ax.set_xlabel('x') del op_samples plt.show()  %% Cell type:markdown id: tags: ## The Neutral Field To demonstrate the effect of all parameters, first a 'neutral' set of parameters is defined which then are varied one by one, showing the effect of the variation on the generated field realizations and the underlying power spectrum from which they were drawn. As a neutral field, a model with a white power spectrum and vanishing spectral power was chosen. %% Cell type:code id: tags:  python # Neutral model parameters yielding a quasi-constant field cf_make_pars = { 'offset_mean': 0., 'offset_std': (1e-3, 1e-16), 'prefix': '' } cf_x_fluct_pars = { 'target_subdomain': x_space, 'fluctuations': (1e-3, 1e-16), 'flexibility': (1e-3, 1e-16), 'asperity': (1e-3, 1e-16), 'loglogavgslope': (0., 1e-16) } init_model(cf_make_pars, cf_x_fluct_pars)  %% Cell type:code id: tags:  python # Show neutral field eval_model(cf_make_pars, cf_x_fluct_pars, "Neutral Field")  %% Cell type:markdown id: tags: # Parameter Showcases ## The fluctuations parameters of add_fluctuations() determine the **amplitude of variations along the field dimension** for which add_fluctuations is called. fluctuations[0] set the _average_ amplitude of the fields fluctuations along the given dimension,\ fluctuations[1] sets the width and shape of the amplitude distribution. The amplitude is modelled as being log-normally distributed, see The Moment-Matched Log-Normal Distribution above for details. #### fluctuations mean: %% Cell type:code id: tags:  python vary_parameter('fluctuations', [(0.05, 1e-16), (0.5, 1e-16), (2., 1e-16)], samples_vary_in='xi')  %% Cell type:markdown id: tags: #### fluctuations std: %% Cell type:code id: tags:  python vary_parameter('fluctuations', [(1., 0.01), (1., 0.1), (1., 1.)], samples_vary_in='fluctuations') cf_x_fluct_pars['fluctuations'] = (1., 1e-16)  %% Cell type:markdown id: tags: ## The loglogavgslope parameters of add_fluctuations() determine **the slope of the loglog-linear (power law) component of the power spectrum**. The slope is modelled to be normally distributed. #### loglogavgslope mean: %% Cell type:code id: tags:  python vary_parameter('loglogavgslope', [(-6., 1e-16), (-2., 1e-16), (2., 1e-16)], samples_vary_in='xi')  %% Cell type:markdown id: tags: #### loglogavgslope std: %% Cell type:code id: tags:  python vary_parameter('loglogavgslope', [(-2., 0.02), (-2., 0.2), (-2., 2.0)], samples_vary_in='loglogavgslope') cf_x_fluct_pars['loglogavgslope'] = (-2., 1e-16)  %% Cell type:markdown id: tags: ## The flexibility parameters of add_fluctuations() determine **the amplitude of the integrated Wiener process component of the power spectrum** (how strong the power spectrum varies besides the power-law). flexibility[0] sets the _average_ amplitude of the i.g.p. component,\ flexibility[1] sets how much the amplitude can vary.\ These two parameters feed into a moment-matched log-normal distribution model, see above for a demo of its behavior. #### flexibility mean: %% Cell type:code id: tags:  python vary_parameter('flexibility', [(0.4, 1e-16), (4.0, 1e-16), (12.0, 1e-16)], samples_vary_in='spectrum')  %% Cell type:markdown id: tags: #### flexibility std: %% Cell type:code id: tags:  python vary_parameter('flexibility', [(4., 0.02), (4., 0.2), (4., 2.)], samples_vary_in='flexibility') cf_x_fluct_pars['flexibility'] = (4., 1e-16)  %% Cell type:markdown id: tags: ## The asperity parameters of add_fluctuations() asperity determines **how rough the integrated Wiener process component of the power spectrum is**. asperity[0] sets the average roughness, asperity[1] sets how much the roughness can vary.\ These two parameters feed into a moment-matched log-normal distribution model, see above for a demo of its behavior. #### asperity mean: %% Cell type:code id: tags:  python vary_parameter('asperity', [(0.001, 1e-16), (1.0, 1e-16), (5., 1e-16)], samples_vary_in='spectrum')  %% Cell type:markdown id: tags: