From 81cab5d61991f51f3de0426bd3af33b9e01dd561 Mon Sep 17 00:00:00 2001 From: Philipp Arras Date: Thu, 1 Feb 2018 12:11:43 +0100 Subject: [PATCH 1/8] Wiener filter demo Works only with python2 and on commit 1d10be4674a42945f8548f3b68688bf0f0d753fe. --- demos/Wiener Filter.ipynb | 879 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 879 insertions(+) create mode 100644 demos/Wiener Filter.ipynb diff --git a/demos/Wiener Filter.ipynb b/demos/Wiener Filter.ipynb new file mode 100644 index 00000000..eddcb03d --- /dev/null +++ b/demos/Wiener Filter.ipynb @@ -0,0 +1,879 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# A NIFTy demonstration" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## IFT: Big Picture\n", + "IFT starting point:\n", + "\n", + "$$d = Rs+n$$\n", + "\n", + "Typically, $s$ continuous field, $d$ discrete data vector. Particularily, $R$ is not invertible.\n", + "\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, en. raffiniert) is a Python framework in which IFT problems can be tackeled easily.\n", + "\n", + "Main Interfaces:\n", + "\n", + "- **Spaces**: Cartesian, 2-Spheres (Healpix, Gauss-Legendre) and their respective harmonic spaces.\n", + "- **Fields**: Defined on spaces.\n", + "- **Operators**: Acting on spaces." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## Wiener Filter: Formulae\n", + "\n", + "### Assumptions\n", + "\n", + "- $d=Rs+n$, $R$ linear operator.\n", + "- $\\mathcal P (s) = \\mathcal G (s,S)$, $\\mathcal P (n) = \\mathcal G (n,N)$ where $S, N$ are positive definite matrices.\n", + "\n", + "### Posterior\n", + "The Posterior is given by:\n", + "\n", + "$$\\mathcal P (s|d) \\propto P(s,d) = \\mathcal G(d-Rs,N) \\,\\mathcal G(s,S) \\propto \\mathcal G (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", + "\n", + "Let us implement this in NIFTy!" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## Wiener Filter: Example\n", + "\n", + "- One-dimensional signal with powerspectrum: $$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$. Recall: $P(k)$ defines an isotropic and homogeneous $S$.\n", + "- $N = 0.5 \\cdot \\text{id}$.\n", + "- Number data points $N_{pix} = 512$.\n", + "- Response operator:\n", + "$$R_x=\\begin{pmatrix} \\delta(x-0)\\\\\\delta(x-1)\\\\\\ldots\\\\ \\delta(x-511) \\end{pmatrix}.$$\n", + "However, the signal space is also discrete on the computer and $R = \\text{id}$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "N_pixels = 512 # Number of pixels\n", + "sigma2 = .5 # Noise variance\n", + "\n", + "def pow_spec(k):\n", + " P0, k0, gamma = [.2, 5, 6]\n", + " return P0 * (1. + (k/k0)**2)**(- gamma / 2)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Wiener Filter: Implementation" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "-" + } + }, + "source": [ + "### Import Modules" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from nifty import (DiagonalOperator, EndomorphicOperator, FFTOperator, Field,\n", + " InvertibleOperatorMixin, PowerSpace, RGSpace,\n", + " create_power_operator, SmoothingOperator, DiagonalProberMixin, Prober)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Implement Propagator" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "class PropagatorOperator(InvertibleOperatorMixin, EndomorphicOperator):\n", + " def __init__(self, R, N, Sh, default_spaces=None):\n", + " super(PropagatorOperator, self).__init__(default_spaces=default_spaces,\n", + " preconditioner=lambda x : fft.adjoint_times(Sh.times(fft.times(x))))\n", + "\n", + " self.R = R\n", + " self.N = N\n", + " self.Sh = Sh\n", + " self._domain = R.domain\n", + " self.fft = FFTOperator(domain=R.domain, target=Sh.domain)\n", + "\n", + " def _inverse_times(self, x, spaces, x0=None):\n", + " return self.R.adjoint_times(self.N.inverse_times(self.R(x))) \\\n", + " + self.fft.adjoint_times(self.Sh.inverse_times(self.fft(x)))\n", + "\n", + " @property\n", + " def domain(self):\n", + " return self._domain\n", + "\n", + " @property\n", + " def unitary(self):\n", + " return False\n", + "\n", + " @property\n", + " def symmetric(self):\n", + " return False\n", + "\n", + " @property\n", + " def self_adjoint(self):\n", + " return True" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "### Conjugate Gradient Preconditioning\n", + "\n", + "- $D$ is defined via:\n", + "$$D^{-1} = \\mathcal F^\\dagger S_h^{-1}\\mathcal F + R^\\dagger N^{-1} R.$$\n", + "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*). \n", + "\n", + "- One can define the *condition number* of a non-singular and normal matrix $A$:\n", + "$$\\kappa (A) := \\frac{|\\lambda_{\\text{max}}|}{|\\lambda_{\\text{min}}|},$$\n", + "where $\\lambda_{\\text{max}}$ and $\\lambda_{\\text{min}}$ are the largest and smallest eigenvalue of $A$, respectively.\n", + "\n", + "- The larger $\\kappa$ the slower Conjugate Gradient.\n", + "\n", + "- By default, conjugate gradient solves: $D^{-1} m = j$ for $m$, where $D^{-1}$ can be bad conditioned. If one knows a non-singular matrix $T$ for which $TD^{-1}$ is better conditioned, one can solve the equivalent problem:\n", + "$$\\tilde A m = \\tilde j,$$\n", + "where $\\tilde A = T D^{-1}$ and $\\tilde j = Tj$.\n", + "\n", + "- In our case $S^{-1}$ is responsible for the bad conditioning of $D$ depending on the chosen power spectrum. Thus, we choose\n", + "\n", + "$$T = \\mathcal F^\\dagger S_h^{-1} \\mathcal F.$$" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Generate Mock data\n", + "\n", + "- Generate a field $s$ and $n$ with given covariances.\n", + "- Calculate $d$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "s_space = RGSpace(N_pixels)\n", + "fft = FFTOperator(s_space)\n", + "h_space = fft.target[0]\n", + "p_space = PowerSpace(h_space)\n", + "\n", + "\n", + "# Operators\n", + "Sh = create_power_operator(h_space, power_spectrum=pow_spec)\n", + "N = DiagonalOperator(s_space, diagonal=sigma2, bare=True)\n", + "R = DiagonalOperator(s_space, diagonal=1.)\n", + "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", + "\n", + "# Fields and data\n", + "sh = Field(p_space, val=pow_spec).power_synthesize(real_signal=True)\n", + "s = fft.adjoint_times(sh)\n", + "n = Field.from_random(domain=s_space, random_type='normal',\n", + " std=np.sqrt(sigma2), mean=0)\n", + "d = R(s) + n\n", + "j = R.adjoint_times(N.inverse_times(d))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Run Wiener Filter" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "m = D(j)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Create Power Spectra of Signal and Reconstruction" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "s_power = sh.power_analyze()\n", + "m_power = fft(m).power_analyze()\n", + "s_power_data = s_power.val.get_full_data().real\n", + "m_power_data = m_power.val.get_full_data().real\n", + "\n", + "# Get signal data and reconstruction data\n", + "s_data = s.val.get_full_data().real\n", + "m_data = m.val.get_full_data().real\n", + "\n", + "d_data = d.val.get_full_data().real" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Signal Reconstruction" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "plt.plot(s_data, 'k', label=\"Signal\", alpha=.5, linewidth=.5)\n", + "plt.plot(d_data, 'k+', label=\"Data\")\n", + "plt.plot(m_data, 'r', label=\"Reconstruction\")\n", + "plt.title(\"Reconstruction\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.plot(s_data - s_data, 'k', label=\"Signal\", alpha=.5, linewidth=.5)\n", + "plt.plot(d_data - s_data, 'k+', label=\"Data\")\n", + "plt.plot(m_data - s_data, 'r', label=\"Reconstruction\")\n", + "plt.axhspan(-np.sqrt(sigma2),np.sqrt(sigma2), facecolor='0.9', alpha=.5)\n", + "plt.title(\"Residuals\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Power Spectrum" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "plt.loglog()\n", + "plt.xlim(1, int(N_pixels/2))\n", + "ymin = min(m_power_data)\n", + "plt.ylim(ymin, 1)\n", + "xs = np.arange(1,int(N_pixels/2),.1)\n", + "plt.plot(xs, pow_spec(xs), label=\"True Power Spectrum\", linewidth=.7, color='k')\n", + "plt.plot(s_power_data, 'k', label=\"Signal\", alpha=.5, linewidth=.5)\n", + "plt.plot(m_power_data, 'r', label=\"Reconstruction\")\n", + "plt.axhline(sigma2 / N_pixels, color=\"k\", linestyle='--', label=\"Noise level\", alpha=.5)\n", + "plt.axhspan(sigma2 / N_pixels, ymin, facecolor='0.9', alpha=.5)\n", + "plt.title(\"Power Spectrum\")\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Wiener Filter on Incomplete Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "# Operators\n", + "Sh = create_power_operator(h_space, power_spectrum=pow_spec)\n", + "N = DiagonalOperator(s_space, diagonal=sigma2, bare=True)\n", + "# R is defined below\n", + "\n", + "# Fields\n", + "sh = Field(p_space, val=pow_spec).power_synthesize(real_signal=True)\n", + "s = fft.adjoint_times(sh)\n", + "n = Field.from_random(domain=s_space, random_type='normal',\n", + " std=np.sqrt(sigma2), mean=0)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "### Partially Lose Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "l = int(N_pixels * 0.2)\n", + "h = int(N_pixels * 0.2 * 4)\n", + "\n", + "mask = Field(s_space, val=1)\n", + "mask.val[ l : h] = 0\n", + "\n", + "R = DiagonalOperator(s_space, diagonal = mask)\n", + "n.val[l:h] = 0\n", + "\n", + "d = R(s) + n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", + "j = R.adjoint_times(N.inverse_times(d))\n", + "m = D(j)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Compute Uncertainty\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "class DiagonalProber(DiagonalProberMixin, Prober):\n", + " def __init__(self, *args, **kwargs):\n", + " super(DiagonalProber,self).__init__(*args, **kwargs)\n", + "\n", + "diagProber = DiagonalProber(domain=s_space, probe_dtype=np.complex, probe_count=200)\n", + "diagProber(D)\n", + "m_var = Field(s_space,val=diagProber.diagonal.val).weight(-1)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "### Get data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "s_power = sh.power_analyze()\n", + "m_power = fft(m).power_analyze()\n", + "s_power_data = s_power.val.get_full_data().real\n", + "m_power_data = m_power.val.get_full_data().real\n", + "\n", + "# Get signal data and reconstruction data\n", + "s_data = s.val.get_full_data().real\n", + "m_data = m.val.get_full_data().real\n", + "m_var_data = m_var.val.get_full_data().real\n", + "uncertainty = np.sqrt(np.abs(m_var_data))\n", + "\n", + "d_data = d.val.get_full_data().real\n", + "\n", + "# Set lost data to NaN for proper plotting\n", + "d_data[d_data == 0] = np.nan" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(15,10))\n", + "plt.plot(s_data, 'k', label=\"Signal\", alpha=.5, linewidth=1)\n", + "plt.plot(d_data, 'k+', label=\"Data\", alpha=1)\n", + "plt.axvspan(l, h, facecolor='0.8', alpha=.5)\n", + "plt.title(\"Incomplete Data\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "fig" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(15,10))\n", + "plt.plot(s_data, 'k', label=\"Signal\", alpha=1, linewidth=1)\n", + "plt.plot(d_data, 'k+', label=\"Data\", alpha=.5)\n", + "plt.plot(m_data, 'r', label=\"Reconstruction\")\n", + "plt.axvspan(l, h, facecolor='0.8', alpha=.5)\n", + "plt.fill_between(range(N_pixels), m_data - uncertainty, m_data + uncertainty, facecolor='0')\n", + "plt.title(\"Reconstruction of incomplete data\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# 2d Example" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "N_pixels = 256 # Number of pixels\n", + "sigma2 = 1000 # Noise variance\n", + "\n", + "\n", + "def pow_spec(k):\n", + " P0, k0, gamma = [.2, 20, 4]\n", + " return P0 * (1. + (k/k0)**2)**(- gamma / 2)\n", + "\n", + "\n", + "s_space = RGSpace([N_pixels, N_pixels])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "fft = FFTOperator(s_space)\n", + "h_space = fft.target[0]\n", + "p_space = PowerSpace(h_space)\n", + "\n", + "# Operators\n", + "Sh = create_power_operator(h_space, power_spectrum=pow_spec)\n", + "N = DiagonalOperator(s_space, diagonal=sigma2, bare=True)\n", + "R = SmoothingOperator(s_space, sigma=.01)\n", + "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", + "\n", + "# Fields and data\n", + "sh = Field(p_space, val=pow_spec).power_synthesize(real_signal=True)\n", + "s = fft.adjoint_times(sh)\n", + "n = Field.from_random(domain=s_space, random_type='normal',\n", + " std=np.sqrt(sigma2), mean=0)\n", + "\n", + "# Lose some data\n", + "\n", + "l = int(N_pixels * 0.2)\n", + "h = int(N_pixels * 0.2 * 2)\n", + "\n", + "mask = Field(s_space, val=1)\n", + "mask.val[l:h,l:h] = 0\n", + "\n", + "R = DiagonalOperator(s_space, diagonal = mask)\n", + "n.val[l:h, l:h] = 0\n", + "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", + "\n", + "d = R(s) + n\n", + "j = R.adjoint_times(N.inverse_times(d))\n", + "\n", + "# Run Wiener filter\n", + "m = D(j)\n", + "\n", + "# Uncertainty\n", + "diagProber = DiagonalProber(domain=s_space, probe_dtype=np.complex, probe_count=10)\n", + "diagProber(D)\n", + "m_var = Field(s_space, val=diagProber.diagonal.val).weight(-1)\n", + "\n", + "# Get data\n", + "s_power = sh.power_analyze()\n", + "m_power = fft(m).power_analyze()\n", + "s_power_data = s_power.val.get_full_data().real\n", + "m_power_data = m_power.val.get_full_data().real\n", + "s_data = s.val.get_full_data().real\n", + "m_data = m.val.get_full_data().real\n", + "m_var_data = m_var.val.get_full_data().real\n", + "d_data = d.val.get_full_data().real\n", + "\n", + "uncertainty = np.sqrt(np.abs(m_var_data))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "cm = ['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", + "\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", + " vmax=ma)\n", + " ax.set_title(caption.pop(0))\n", + "\n", + "fig.subplots_adjust(right=0.8)\n", + "cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n", + "fig.colorbar(im, cax=cbar_ax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "mi = np.min(s_data)\n", + "ma = np.max(s_data)\n", + "\n", + "fig, axes = plt.subplots(2, 2, figsize=(15, 15))\n", + "\n", + "data = [s_data, m_data, s_data - m_data, uncertainty]\n", + "caption = [\"Signal\", \"Reconstruction\", \"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", + " ax.set_title(caption.pop(0))\n", + "\n", + "fig.subplots_adjust(right=0.8)\n", + "cbar_ax = fig.add_axes([.85, 0.15, 0.05, 0.7])\n", + "fig.colorbar(im, cax=cbar_ax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Is the uncertainty map reliable?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "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", + "fig = plt.figure()\n", + "plt.imshow(precise.astype(float), cmap=\"brg\")\n", + "plt.colorbar()\n", + "fig" + ] + }, + { + "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", + "\n", + "commit 1d10be4674a42945f8548f3b68688bf0f0d753fe\n", + "\n", + "NIFTy v3 **not (yet) stable!**" + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} -- GitLab From 35e632b032eec6cdb1ad04f2effd052cbb29b441 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Sun, 4 Feb 2018 16:10:59 +0100 Subject: [PATCH 2/8] make the 1D part (sort of) work with NIFTy4 --- demos/Wiener Filter.ipynb | 187 ++++++++++++++++---------------------- 1 file changed, 80 insertions(+), 107 deletions(-) diff --git a/demos/Wiener Filter.ipynb b/demos/Wiener Filter.ipynb index eddcb03d..0ad954ca 100644 --- a/demos/Wiener Filter.ipynb +++ b/demos/Wiener Filter.ipynb @@ -93,7 +93,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "slideshow": { "slide_type": "-" } @@ -101,7 +100,6 @@ "outputs": [], "source": [ "N_pixels = 512 # Number of pixels\n", - "sigma2 = .5 # Noise variance\n", "\n", "def pow_spec(k):\n", " P0, k0, gamma = [.2, 5, 6]\n", @@ -140,11 +138,11 @@ }, "outputs": [], "source": [ - "import matplotlib.pyplot as plt\n", "import numpy as np\n", - "from nifty import (DiagonalOperator, EndomorphicOperator, FFTOperator, Field,\n", - " InvertibleOperatorMixin, PowerSpace, RGSpace,\n", - " create_power_operator, SmoothingOperator, DiagonalProberMixin, Prober)" + "np.random.seed(42)\n", + "import nifty4 as ift\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline" ] }, { @@ -162,43 +160,20 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ - "class PropagatorOperator(InvertibleOperatorMixin, EndomorphicOperator):\n", - " def __init__(self, R, N, Sh, default_spaces=None):\n", - " super(PropagatorOperator, self).__init__(default_spaces=default_spaces,\n", - " preconditioner=lambda x : fft.adjoint_times(Sh.times(fft.times(x))))\n", - "\n", - " self.R = R\n", - " self.N = N\n", - " self.Sh = Sh\n", - " self._domain = R.domain\n", - " self.fft = FFTOperator(domain=R.domain, target=Sh.domain)\n", - "\n", - " def _inverse_times(self, x, spaces, x0=None):\n", - " return self.R.adjoint_times(self.N.inverse_times(self.R(x))) \\\n", - " + self.fft.adjoint_times(self.Sh.inverse_times(self.fft(x)))\n", - "\n", - " @property\n", - " def domain(self):\n", - " return self._domain\n", - "\n", - " @property\n", - " def unitary(self):\n", - " return False\n", - "\n", - " @property\n", - " def symmetric(self):\n", - " return False\n", - "\n", - " @property\n", - " def self_adjoint(self):\n", - " return True" + "def PropagatorOperator(R, N, Sh):\n", + " IC = ift.GradientNormController(name=\"inverter\", iteration_limit=50000,\n", + " tol_abs_gradnorm=0.1)\n", + " inverter = ift.ConjugateGradient(controller=IC)\n", + " D = (R.adjoint*N.inverse*R + Sh.inverse).inverse\n", + " # MR FIXME: we can/should provide a preconditioner here as well!\n", + " return ift.InversionEnabler(D, inverter)\n", + " #return ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter).inverse\n" ] }, { @@ -247,30 +222,33 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ - "s_space = RGSpace(N_pixels)\n", - "fft = FFTOperator(s_space)\n", - "h_space = fft.target[0]\n", - "p_space = PowerSpace(h_space)\n", - "\n", + "s_space = ift.RGSpace(N_pixels)\n", + "h_space = s_space.get_default_codomain()\n", + "HT = ift.HarmonicTransformOperator(h_space, target=s_space)\n", + "p_space = ift.PowerSpace(h_space)\n", "\n", "# Operators\n", - "Sh = create_power_operator(h_space, power_spectrum=pow_spec)\n", - "N = DiagonalOperator(s_space, diagonal=sigma2, bare=True)\n", - "R = DiagonalOperator(s_space, diagonal=1.)\n", - "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", + "Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)\n", + "R = HT #*ift.create_harmonic_smoothing_operator((h_space,), 0, 0.02)\n", "\n", "# Fields and data\n", - "sh = Field(p_space, val=pow_spec).power_synthesize(real_signal=True)\n", - "s = fft.adjoint_times(sh)\n", - "n = Field.from_random(domain=s_space, random_type='normal',\n", - " std=np.sqrt(sigma2), mean=0)\n", - "d = R(s) + n\n", - "j = R.adjoint_times(N.inverse_times(d))" + "sh = ift.power_synthesize(ift.PS_field(p_space, pow_spec),real_signal=True)\n", + "noiseless_data=R(sh)\n", + "signal_to_noise = 5\n", + "noise_amplitude = noiseless_data.std()/signal_to_noise\n", + "N = ift.ScalingOperator(noise_amplitude**2, s_space)\n", + "\n", + "n = ift.Field.from_random(domain=s_space, random_type='normal',\n", + " std=noise_amplitude, mean=0)\n", + "ift.plot(n)\n", + "d = noiseless_data + n\n", + "ift.plot(d)\n", + "j = R.adjoint_times(N.inverse_times(d))\n", + "ift.plot(HT(j))\n", + "D = PropagatorOperator(R=R, N=N, Sh=Sh)" ] }, { @@ -288,7 +266,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "slideshow": { "slide_type": "-" } @@ -313,23 +290,22 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ - "s_power = sh.power_analyze()\n", - "m_power = fft(m).power_analyze()\n", - "s_power_data = s_power.val.get_full_data().real\n", - "m_power_data = m_power.val.get_full_data().real\n", + "s_power = ift.power_analyze(sh)\n", + "m_power = ift.power_analyze(m)\n", + "s_power_data = s_power.val.real\n", + "m_power_data = m_power.val.real\n", "\n", "# Get signal data and reconstruction data\n", - "s_data = s.val.get_full_data().real\n", - "m_data = m.val.get_full_data().real\n", + "s_data = HT(sh).val.real\n", + "m_data = HT(m).val.real\n", "\n", - "d_data = d.val.get_full_data().real" + "d_data = d.val.real" ] }, { @@ -375,7 +351,7 @@ "plt.plot(s_data - s_data, 'k', label=\"Signal\", alpha=.5, linewidth=.5)\n", "plt.plot(d_data - s_data, 'k+', label=\"Data\")\n", "plt.plot(m_data - s_data, 'r', label=\"Reconstruction\")\n", - "plt.axhspan(-np.sqrt(sigma2),np.sqrt(sigma2), facecolor='0.9', alpha=.5)\n", + "plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5)\n", "plt.title(\"Residuals\")\n", "plt.legend()\n", "plt.show()" @@ -410,8 +386,8 @@ "plt.plot(xs, pow_spec(xs), label=\"True Power Spectrum\", linewidth=.7, color='k')\n", "plt.plot(s_power_data, 'k', label=\"Signal\", alpha=.5, linewidth=.5)\n", "plt.plot(m_power_data, 'r', label=\"Reconstruction\")\n", - "plt.axhline(sigma2 / N_pixels, color=\"k\", linestyle='--', label=\"Noise level\", alpha=.5)\n", - "plt.axhspan(sigma2 / N_pixels, ymin, facecolor='0.9', alpha=.5)\n", + "plt.axhline(noise_amplitude**2 / N_pixels, color=\"k\", linestyle='--', label=\"Noise level\", alpha=.5)\n", + "plt.axhspan(noise_amplitude**2 / N_pixels, ymin, facecolor='0.9', alpha=.5)\n", "plt.title(\"Power Spectrum\")\n", "plt.legend()\n", "plt.show()" @@ -432,7 +408,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "slideshow": { "slide_type": "skip" } @@ -440,15 +415,15 @@ "outputs": [], "source": [ "# Operators\n", - "Sh = create_power_operator(h_space, power_spectrum=pow_spec)\n", - "N = DiagonalOperator(s_space, diagonal=sigma2, bare=True)\n", + "Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)\n", + "N = ift.ScalingOperator(noise_amplitude**2,s_space)\n", "# R is defined below\n", "\n", "# Fields\n", - "sh = Field(p_space, val=pow_spec).power_synthesize(real_signal=True)\n", - "s = fft.adjoint_times(sh)\n", - "n = Field.from_random(domain=s_space, random_type='normal',\n", - " std=np.sqrt(sigma2), mean=0)" + "sh = ift.power_synthesize(ift.PS_field(p_space,pow_spec),real_signal=True)\n", + "s = HT(sh)\n", + "n = ift.Field.from_random(domain=s_space, random_type='normal',\n", + " std=noise_amplitude, mean=0)" ] }, { @@ -466,7 +441,6 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "slideshow": { "slide_type": "-" } @@ -474,22 +448,21 @@ "outputs": [], "source": [ "l = int(N_pixels * 0.2)\n", - "h = int(N_pixels * 0.2 * 4)\n", + "h = int(N_pixels * 0.2 * 2)\n", "\n", - "mask = Field(s_space, val=1)\n", + "mask = ift.Field(s_space, val=1)\n", "mask.val[ l : h] = 0\n", "\n", - "R = DiagonalOperator(s_space, diagonal = mask)\n", + "R = ift.DiagonalOperator(mask)*HT\n", "n.val[l:h] = 0\n", "\n", - "d = R(s) + n" + "d = R(sh) + n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "slideshow": { "slide_type": "skip" } @@ -516,17 +489,21 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true + "scrolled": true }, "outputs": [], "source": [ - "class DiagonalProber(DiagonalProberMixin, Prober):\n", - " def __init__(self, *args, **kwargs):\n", - " super(DiagonalProber,self).__init__(*args, **kwargs)\n", + "sc = ift.probing.utils.StatCalculator()\n", "\n", - "diagProber = DiagonalProber(domain=s_space, probe_dtype=np.complex, probe_count=200)\n", - "diagProber(D)\n", - "m_var = Field(s_space,val=diagProber.diagonal.val).weight(-1)" + "IC = ift.GradientNormController(name=\"inverter\", iteration_limit=50000,\n", + " tol_abs_gradnorm=0.1)\n", + "inverter = ift.ConjugateGradient(controller=IC)\n", + "curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter)\n", + "\n", + "for i in range(200):\n", + " sc.add(HT(curv.generate_posterior_sample()))\n", + "\n", + "m_var = sc.var" ] }, { @@ -544,25 +521,24 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ - "s_power = sh.power_analyze()\n", - "m_power = fft(m).power_analyze()\n", - "s_power_data = s_power.val.get_full_data().real\n", - "m_power_data = m_power.val.get_full_data().real\n", + "s_power = ift.power_analyze(sh)\n", + "m_power = ift.power_analyze(m)\n", + "s_power_data = s_power.val.real\n", + "m_power_data = m_power.val.real\n", "\n", "# Get signal data and reconstruction data\n", - "s_data = s.val.get_full_data().real\n", - "m_data = m.val.get_full_data().real\n", - "m_var_data = m_var.val.get_full_data().real\n", + "s_data = s.val.real\n", + "m_data = HT(m).val.real\n", + "m_var_data = m_var.val.real\n", "uncertainty = np.sqrt(np.abs(m_var_data))\n", - "\n", - "d_data = d.val.get_full_data().real\n", + "ift.plot(ift.sqrt(m_var))\n", + "d_data = d.val.real\n", "\n", "# Set lost data to NaN for proper plotting\n", "d_data[d_data == 0] = np.nan" @@ -646,9 +622,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "N_pixels = 256 # Number of pixels\n", @@ -660,14 +634,13 @@ " return P0 * (1. + (k/k0)**2)**(- gamma / 2)\n", "\n", "\n", - "s_space = RGSpace([N_pixels, N_pixels])" + "s_space = ift.RGSpace([N_pixels, N_pixels])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": true, "slideshow": { "slide_type": "skip" } @@ -857,21 +830,21 @@ "metadata": { "celltoolbar": "Slideshow", "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 2", "language": "python", - "name": "python3" + "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 3 + "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.6.4" + "pygments_lexer": "ipython2", + "version": "2.7.12" } }, "nbformat": 4, -- GitLab From c92db8accc9bb87bee671cf3d4fc986e0d541660 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Sun, 4 Feb 2018 17:57:05 +0100 Subject: [PATCH 3/8] tweaks --- demos/Wiener Filter.ipynb | 93 +++++++++++++++------------------------ 1 file changed, 36 insertions(+), 57 deletions(-) diff --git a/demos/Wiener Filter.ipynb b/demos/Wiener Filter.ipynb index 0ad954ca..f8cf5d3e 100644 --- a/demos/Wiener Filter.ipynb +++ b/demos/Wiener Filter.ipynb @@ -24,20 +24,20 @@ "\n", "$$d = Rs+n$$\n", "\n", - "Typically, $s$ continuous field, $d$ discrete data vector. Particularily, $R$ is not invertible.\n", + "Typically, $s$ is a continuous field, $d$ a discrete data vector. Particularly, $R$ is not invertible.\n", "\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, en. raffiniert) is a Python framework in which IFT problems can be tackeled easily.\n", + "NIFTy (Numerical Information Field Theory) is a Python framework in which IFT problems can be tackled easily.\n", "\n", "Main Interfaces:\n", "\n", "- **Spaces**: Cartesian, 2-Spheres (Healpix, Gauss-Legendre) and their respective harmonic spaces.\n", "- **Fields**: Defined on spaces.\n", - "- **Operators**: Acting on spaces." + "- **Operators**: Acting on fields." ] }, { @@ -80,7 +80,7 @@ "source": [ "## Wiener Filter: Example\n", "\n", - "- One-dimensional signal with powerspectrum: $$P(k) = P_0\\,\\left(1+\\left(\\frac{k}{k_0}\\right)^2\\right)^{-\\gamma /2},$$\n", + "- One-dimensional signal with power spectrum: $$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$. Recall: $P(k)$ defines an isotropic and homogeneous $S$.\n", "- $N = 0.5 \\cdot \\text{id}$.\n", "- Number data points $N_{pix} = 512$.\n", @@ -102,8 +102,8 @@ "N_pixels = 512 # Number of pixels\n", "\n", "def pow_spec(k):\n", - " P0, k0, gamma = [.2, 5, 6]\n", - " return P0 * (1. + (k/k0)**2)**(- gamma / 2)" + " P0, k0, gamma = [.2, 5, 4]\n", + " return P0 / ((1. + (k/k0)**2)**(gamma / 2))" ] }, { @@ -139,7 +139,7 @@ "outputs": [], "source": [ "import numpy as np\n", - "np.random.seed(42)\n", + "np.random.seed(40)\n", "import nifty4 as ift\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" @@ -172,8 +172,7 @@ " inverter = ift.ConjugateGradient(controller=IC)\n", " D = (R.adjoint*N.inverse*R + Sh.inverse).inverse\n", " # MR FIXME: we can/should provide a preconditioner here as well!\n", - " return ift.InversionEnabler(D, inverter)\n", - " #return ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter).inverse\n" + " return ift.InversionEnabler(D, inverter)\n" ] }, { @@ -196,7 +195,7 @@ "\n", "- The larger $\\kappa$ the slower Conjugate Gradient.\n", "\n", - "- By default, conjugate gradient solves: $D^{-1} m = j$ for $m$, where $D^{-1}$ can be bad conditioned. If one knows a non-singular matrix $T$ for which $TD^{-1}$ is better conditioned, one can solve the equivalent problem:\n", + "- By default, conjugate gradient solves: $D^{-1} m = j$ for $m$, where $D^{-1}$ can be badly conditioned. If one knows a non-singular matrix $T$ for which $TD^{-1}$ is better conditioned, one can solve the equivalent problem:\n", "$$\\tilde A m = \\tilde j,$$\n", "where $\\tilde A = T D^{-1}$ and $\\tilde j = Tj$.\n", "\n", @@ -237,17 +236,13 @@ "# Fields and data\n", "sh = ift.power_synthesize(ift.PS_field(p_space, pow_spec),real_signal=True)\n", "noiseless_data=R(sh)\n", - "signal_to_noise = 5\n", - "noise_amplitude = noiseless_data.std()/signal_to_noise\n", + "noise_amplitude = np.sqrt(0.05)\n", "N = ift.ScalingOperator(noise_amplitude**2, s_space)\n", "\n", "n = ift.Field.from_random(domain=s_space, random_type='normal',\n", - " std=noise_amplitude, mean=0)\n", - "ift.plot(n)\n", + " std=noise_amplitude, mean=0)\n", "d = noiseless_data + n\n", - "ift.plot(d)\n", "j = R.adjoint_times(N.inverse_times(d))\n", - "ift.plot(HT(j))\n", "D = PropagatorOperator(R=R, N=N, Sh=Sh)" ] }, @@ -329,7 +324,7 @@ }, "outputs": [], "source": [ - "plt.plot(s_data, 'k', label=\"Signal\", alpha=.5, linewidth=.5)\n", + "plt.plot(s_data, 'g', label=\"Signal\")\n", "plt.plot(d_data, 'k+', label=\"Data\")\n", "plt.plot(m_data, 'r', label=\"Reconstruction\")\n", "plt.title(\"Reconstruction\")\n", @@ -348,7 +343,7 @@ "outputs": [], "source": [ "plt.figure()\n", - "plt.plot(s_data - s_data, 'k', label=\"Signal\", alpha=.5, linewidth=.5)\n", + "plt.plot(s_data - s_data, 'g', label=\"Signal\")\n", "plt.plot(d_data - s_data, 'k+', label=\"Data\")\n", "plt.plot(m_data - s_data, 'r', label=\"Reconstruction\")\n", "plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5)\n", @@ -537,7 +532,6 @@ "m_data = HT(m).val.real\n", "m_var_data = m_var.val.real\n", "uncertainty = np.sqrt(np.abs(m_var_data))\n", - "ift.plot(ift.sqrt(m_var))\n", "d_data = d.val.real\n", "\n", "# Set lost data to NaN for proper plotting\n", @@ -562,19 +556,6 @@ "plt.legend()" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "outputs": [], - "source": [ - "fig" - ] - }, { "cell_type": "code", "execution_count": null, @@ -595,19 +576,6 @@ "plt.legend()" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "outputs": [], - "source": [ - "fig" - ] - }, { "cell_type": "markdown", "metadata": { @@ -647,20 +615,20 @@ }, "outputs": [], "source": [ - "fft = FFTOperator(s_space)\n", - "h_space = fft.target[0]\n", - "p_space = PowerSpace(h_space)\n", + "h_space = s_space.get_default_codomain()\n", + "fft = ift.FFTOperator(s_space,h_space)\n", + "p_space = ift.PowerSpace(h_space)\n", "\n", "# Operators\n", - "Sh = create_power_operator(h_space, power_spectrum=pow_spec)\n", - "N = DiagonalOperator(s_space, diagonal=sigma2, bare=True)\n", - "R = SmoothingOperator(s_space, sigma=.01)\n", - "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", + "Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)\n", + "N = ift.ScalingOperator(sigma2,s_space)\n", + "R = ift.FFTSmoothingOperator(s_space, sigma=.01)\n", + "#D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", "\n", "# Fields and data\n", - "sh = Field(p_space, val=pow_spec).power_synthesize(real_signal=True)\n", + "sh = ift.power_synthesize(ift.PS_field(p_space,pow_spec),real_signal=True)\n", "s = fft.adjoint_times(sh)\n", - "n = Field.from_random(domain=s_space, random_type='normal',\n", + "n = ift.Field.from_random(domain=s_space, random_type='normal',\n", " std=np.sqrt(sigma2), mean=0)\n", "\n", "# Lose some data\n", @@ -668,12 +636,12 @@ "l = int(N_pixels * 0.2)\n", "h = int(N_pixels * 0.2 * 2)\n", "\n", - "mask = Field(s_space, val=1)\n", + "mask = ift.Field(s_space, val=1)\n", "mask.val[l:h,l:h] = 0\n", "\n", - "R = DiagonalOperator(s_space, diagonal = mask)\n", + "R = ift.DiagonalOperator(mask)\n", "n.val[l:h, l:h] = 0\n", - "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", + "D = PropagatorOperator(R=R, N=N, Sh=fft.inverse*Sh*fft)\n", "\n", "d = R(s) + n\n", "j = R.adjoint_times(N.inverse_times(d))\n", @@ -682,6 +650,17 @@ "m = D(j)\n", "\n", "# Uncertainty\n", + "sc = ift.probing.utils.StatCalculator()\n", + "\n", + "IC = ift.GradientNormController(name=\"inverter\", iteration_limit=50000,\n", + " tol_abs_gradnorm=0.1)\n", + "inverter = ift.ConjugateGradient(controller=IC)\n", + "curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,fft.inverse*Sh*fft,inverter)\n", + "\n", + "for i in range(20):\n", + " sc.add(HT(curv.generate_posterior_sample()))\n", + "\n", + "m_var = sc.var\n", "diagProber = DiagonalProber(domain=s_space, probe_dtype=np.complex, probe_count=10)\n", "diagProber(D)\n", "m_var = Field(s_space, val=diagProber.diagonal.val).weight(-1)\n", -- GitLab From a4c2996b4bdc85e01a36c076d55f84cdfbc33dd3 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Sun, 4 Feb 2018 21:59:48 +0100 Subject: [PATCH 4/8] fix 2D part of notebook --- demos/Wiener Filter.ipynb | 67 ++++++++++++--------------------------- 1 file changed, 21 insertions(+), 46 deletions(-) diff --git a/demos/Wiener Filter.ipynb b/demos/Wiener Filter.ipynb index f8cf5d3e..7f2d61c6 100644 --- a/demos/Wiener Filter.ipynb +++ b/demos/Wiener Filter.ipynb @@ -490,12 +490,13 @@ "source": [ "sc = ift.probing.utils.StatCalculator()\n", "\n", - "IC = ift.GradientNormController(name=\"inverter\", iteration_limit=50000,\n", + "IC = ift.GradientNormController(iteration_limit=50000,\n", " tol_abs_gradnorm=0.1)\n", "inverter = ift.ConjugateGradient(controller=IC)\n", "curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter)\n", "\n", "for i in range(200):\n", + " print i\n", " sc.add(HT(curv.generate_posterior_sample()))\n", "\n", "m_var = sc.var" @@ -594,11 +595,11 @@ "outputs": [], "source": [ "N_pixels = 256 # Number of pixels\n", - "sigma2 = 1000 # Noise variance\n", + "sigma2 = 10. # Noise variance\n", "\n", "\n", "def pow_spec(k):\n", - " P0, k0, gamma = [.2, 20, 4]\n", + " P0, k0, gamma = [.2, 5, 4]\n", " return P0 * (1. + (k/k0)**2)**(- gamma / 2)\n", "\n", "\n", @@ -616,7 +617,7 @@ "outputs": [], "source": [ "h_space = s_space.get_default_codomain()\n", - "fft = ift.FFTOperator(s_space,h_space)\n", + "HT = ift.HarmonicTransformOperator(h_space,s_space)\n", "p_space = ift.PowerSpace(h_space)\n", "\n", "# Operators\n", @@ -627,7 +628,6 @@ "\n", "# Fields and data\n", "sh = ift.power_synthesize(ift.PS_field(p_space,pow_spec),real_signal=True)\n", - "s = fft.adjoint_times(sh)\n", "n = ift.Field.from_random(domain=s_space, random_type='normal',\n", " std=np.sqrt(sigma2), mean=0)\n", "\n", @@ -639,11 +639,11 @@ "mask = ift.Field(s_space, val=1)\n", "mask.val[l:h,l:h] = 0\n", "\n", - "R = ift.DiagonalOperator(mask)\n", + "R = ift.DiagonalOperator(mask)*HT\n", "n.val[l:h, l:h] = 0\n", - "D = PropagatorOperator(R=R, N=N, Sh=fft.inverse*Sh*fft)\n", + "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", "\n", - "d = R(s) + n\n", + "d = R(sh) + n\n", "j = R.adjoint_times(N.inverse_times(d))\n", "\n", "# Run Wiener filter\n", @@ -652,28 +652,26 @@ "# Uncertainty\n", "sc = ift.probing.utils.StatCalculator()\n", "\n", - "IC = ift.GradientNormController(name=\"inverter\", iteration_limit=50000,\n", - " tol_abs_gradnorm=0.1)\n", + "IC = ift.GradientNormController(iteration_limit=50000,\n", + " tol_abs_gradnorm=0.1)\n", "inverter = ift.ConjugateGradient(controller=IC)\n", - "curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,fft.inverse*Sh*fft,inverter)\n", + "curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter)\n", "\n", "for i in range(20):\n", + " print i\n", " sc.add(HT(curv.generate_posterior_sample()))\n", "\n", "m_var = sc.var\n", - "diagProber = DiagonalProber(domain=s_space, probe_dtype=np.complex, probe_count=10)\n", - "diagProber(D)\n", - "m_var = Field(s_space, val=diagProber.diagonal.val).weight(-1)\n", "\n", "# Get data\n", - "s_power = sh.power_analyze()\n", - "m_power = fft(m).power_analyze()\n", - "s_power_data = s_power.val.get_full_data().real\n", - "m_power_data = m_power.val.get_full_data().real\n", - "s_data = s.val.get_full_data().real\n", - "m_data = m.val.get_full_data().real\n", - "m_var_data = m_var.val.get_full_data().real\n", - "d_data = d.val.get_full_data().real\n", + "s_power = ift.power_analyze(sh)\n", + "m_power = ift.power_analyze(m)\n", + "s_power_data = s_power.val.real\n", + "m_power_data = m_power.val.real\n", + "s_data = HT(sh).val.real\n", + "m_data = HT(m).val.real\n", + "m_var_data = m_var.val.real\n", + "d_data = d.val.real\n", "\n", "uncertainty = np.sqrt(np.abs(m_var_data))" ] @@ -708,15 +706,6 @@ "fig.colorbar(im, cax=cbar_ax)" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig" - ] - }, { "cell_type": "code", "execution_count": null, @@ -744,19 +733,6 @@ "fig.colorbar(im, cax=cbar_ax)" ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "outputs": [], - "source": [ - "fig" - ] - }, { "cell_type": "markdown", "metadata": { @@ -783,8 +759,7 @@ "\n", "fig = plt.figure()\n", "plt.imshow(precise.astype(float), cmap=\"brg\")\n", - "plt.colorbar()\n", - "fig" + "plt.colorbar()" ] }, { -- GitLab From 759a13ce46c7b21bcf267b77b39b3951b74790ad Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Sun, 4 Feb 2018 22:09:54 +0100 Subject: [PATCH 5/8] tweak docs --- demos/Wiener Filter.ipynb | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/demos/Wiener Filter.ipynb b/demos/Wiener Filter.ipynb index 7f2d61c6..90cdf2f4 100644 --- a/demos/Wiener Filter.ipynb +++ b/demos/Wiener Filter.ipynb @@ -82,11 +82,11 @@ "\n", "- One-dimensional signal with power spectrum: $$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$. Recall: $P(k)$ defines an isotropic and homogeneous $S$.\n", - "- $N = 0.5 \\cdot \\text{id}$.\n", - "- Number data points $N_{pix} = 512$.\n", + "- $N = 0.05 \\cdot \\text{id}$.\n", + "- Number of data points $N_{pix} = 512$.\n", + "- reconstruction in harmonic space.\n", "- Response operator:\n", - "$$R_x=\\begin{pmatrix} \\delta(x-0)\\\\\\delta(x-1)\\\\\\ldots\\\\ \\delta(x-511) \\end{pmatrix}.$$\n", - "However, the signal space is also discrete on the computer and $R = \\text{id}$." + "$$R = FFT(\\text{harmonic} \\rightarrow \\text{position})$$\n" ] }, { @@ -167,7 +167,7 @@ "outputs": [], "source": [ "def PropagatorOperator(R, N, Sh):\n", - " IC = ift.GradientNormController(name=\"inverter\", iteration_limit=50000,\n", + " IC = ift.GradientNormController(iteration_limit=50000,\n", " tol_abs_gradnorm=0.1)\n", " inverter = ift.ConjugateGradient(controller=IC)\n", " D = (R.adjoint*N.inverse*R + Sh.inverse).inverse\n", @@ -186,9 +186,10 @@ "### Conjugate Gradient Preconditioning\n", "\n", "- $D$ is defined via:\n", - "$$D^{-1} = \\mathcal F^\\dagger S_h^{-1}\\mathcal F + R^\\dagger N^{-1} R.$$\n", + "$$D^{-1} = \\mathcal S_h^{-1} + R^\\dagger N^{-1} R.$$\n", "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*). \n", "\n", + "" ] }, { -- GitLab From 6437bc49d1e50be35b80288e1d8111c3ec28ba76 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Sun, 4 Feb 2018 22:44:15 +0100 Subject: [PATCH 6/8] tweaks --- demos/Wiener Filter.ipynb | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/demos/Wiener Filter.ipynb b/demos/Wiener Filter.ipynb index 90cdf2f4..cf92a9ea 100644 --- a/demos/Wiener Filter.ipynb +++ b/demos/Wiener Filter.ipynb @@ -82,11 +82,11 @@ "\n", "- One-dimensional signal with power spectrum: $$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$. Recall: $P(k)$ defines an isotropic and homogeneous $S$.\n", - "- $N = 0.05 \\cdot \\text{id}$.\n", + "- $N = 0.2 \\cdot \\mathbb{1}$.\n", "- Number of data points $N_{pix} = 512$.\n", "- reconstruction in harmonic space.\n", "- Response operator:\n", - "$$R = FFT(\\text{harmonic} \\rightarrow \\text{position})$$\n" + "$$R = FFT_{\\text{harmonic} \\rightarrow \\text{position}}$$\n" ] }, { @@ -238,7 +238,7 @@ "# Fields and data\n", "sh = ift.power_synthesize(ift.PS_field(p_space, pow_spec),real_signal=True)\n", "noiseless_data=R(sh)\n", - "noise_amplitude = np.sqrt(0.05)\n", + "noise_amplitude = np.sqrt(0.2)\n", "N = ift.ScalingOperator(noise_amplitude**2, s_space)\n", "\n", "n = ift.Field.from_random(domain=s_space, random_type='normal',\n", @@ -326,6 +326,7 @@ }, "outputs": [], "source": [ + "plt.figure(figsize=(15,10))\n", "plt.plot(s_data, 'g', label=\"Signal\")\n", "plt.plot(d_data, 'k+', label=\"Data\")\n", "plt.plot(m_data, 'r', label=\"Reconstruction\")\n", @@ -344,7 +345,7 @@ }, "outputs": [], "source": [ - "plt.figure()\n", + "plt.figure(figsize=(15,10))\n", "plt.plot(s_data - s_data, 'g', label=\"Signal\")\n", "plt.plot(d_data - s_data, 'k+', label=\"Data\")\n", "plt.plot(m_data - s_data, 'r', label=\"Reconstruction\")\n", @@ -375,13 +376,14 @@ }, "outputs": [], "source": [ + "plt.figure(figsize=(15,10))\n", "plt.loglog()\n", "plt.xlim(1, int(N_pixels/2))\n", "ymin = min(m_power_data)\n", "plt.ylim(ymin, 1)\n", "xs = np.arange(1,int(N_pixels/2),.1)\n", "plt.plot(xs, pow_spec(xs), label=\"True Power Spectrum\", linewidth=.7, color='k')\n", - "plt.plot(s_power_data, 'k', label=\"Signal\", alpha=.5, linewidth=.5)\n", + "plt.plot(s_power_data, 'g', label=\"Signal\")\n", "plt.plot(m_power_data, 'r', label=\"Reconstruction\")\n", "plt.axhline(noise_amplitude**2 / N_pixels, color=\"k\", linestyle='--', label=\"Noise level\", alpha=.5)\n", "plt.axhspan(noise_amplitude**2 / N_pixels, ymin, facecolor='0.9', alpha=.5)\n", @@ -551,8 +553,8 @@ }, "outputs": [], "source": [ - "fig = plt.figure(figsize=(15,10))\n", - "plt.plot(s_data, 'k', label=\"Signal\", alpha=.5, linewidth=1)\n", + "plt.figure(figsize=(15,10))\n", + "plt.plot(s_data, 'g', label=\"Signal\", linewidth=1)\n", "plt.plot(d_data, 'k+', label=\"Data\", alpha=1)\n", "plt.axvspan(l, h, facecolor='0.8', alpha=.5)\n", "plt.title(\"Incomplete Data\")\n", @@ -570,7 +572,7 @@ "outputs": [], "source": [ "fig = plt.figure(figsize=(15,10))\n", - "plt.plot(s_data, 'k', label=\"Signal\", alpha=1, linewidth=1)\n", + "plt.plot(s_data, 'g', label=\"Signal\", alpha=1, linewidth=1)\n", "plt.plot(d_data, 'k+', label=\"Data\", alpha=.5)\n", "plt.plot(m_data, 'r', label=\"Reconstruction\")\n", "plt.axvspan(l, h, facecolor='0.8', alpha=.5)\n", @@ -597,11 +599,11 @@ "outputs": [], "source": [ "N_pixels = 256 # Number of pixels\n", - "sigma2 = 10. # Noise variance\n", + "sigma2 = 2. # Noise variance\n", "\n", "\n", "def pow_spec(k):\n", - " P0, k0, gamma = [.2, 5, 4]\n", + " P0, k0, gamma = [.2, 2, 4]\n", " return P0 * (1. + (k/k0)**2)**(- gamma / 2)\n", "\n", "\n", @@ -635,8 +637,8 @@ "\n", "# Lose some data\n", "\n", - "l = int(N_pixels * 0.2)\n", - "h = int(N_pixels * 0.2 * 2)\n", + "l = int(N_pixels * 0.33)\n", + "h = int(N_pixels * 0.33 * 2)\n", "\n", "mask = ift.Field(s_space, val=1)\n", "mask.val[l:h,l:h] = 0\n", @@ -759,7 +761,7 @@ "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", - "fig = plt.figure()\n", + "plt.figure(figsize=(15,10))\n", "plt.imshow(precise.astype(float), cmap=\"brg\")\n", "plt.colorbar()" ] @@ -777,9 +779,7 @@ "\n", "https://gitlab.mpcdf.mpg.de/ift/NIFTy\n", "\n", - "commit 1d10be4674a42945f8548f3b68688bf0f0d753fe\n", - "\n", - "NIFTy v3 **not (yet) stable!**" + "NIFTy v4 **more or less stable!**" ] } ], -- GitLab From 9b072104b3e96a648d6025794a6b4ede218abfde Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Mon, 5 Feb 2018 10:56:59 +0100 Subject: [PATCH 7/8] rename --- demos/{Wiener Filter.ipynb => Wiener_Filter.ipynb} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename demos/{Wiener Filter.ipynb => Wiener_Filter.ipynb} (100%) diff --git a/demos/Wiener Filter.ipynb b/demos/Wiener_Filter.ipynb similarity index 100% rename from demos/Wiener Filter.ipynb rename to demos/Wiener_Filter.ipynb -- GitLab From d3f0bf6c0bfeaf2259c7accf3020b799e8e6c7cc Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Mon, 5 Feb 2018 11:06:52 +0100 Subject: [PATCH 8/8] cleanup --- demos/Wiener_Filter.ipynb | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/demos/Wiener_Filter.ipynb b/demos/Wiener_Filter.ipynb index cf92a9ea..a8da9671 100644 --- a/demos/Wiener_Filter.ipynb +++ b/demos/Wiener_Filter.ipynb @@ -166,13 +166,13 @@ }, "outputs": [], "source": [ - "def PropagatorOperator(R, N, Sh):\n", + "def Curvature(R, N, Sh):\n", " IC = ift.GradientNormController(iteration_limit=50000,\n", " tol_abs_gradnorm=0.1)\n", " inverter = ift.ConjugateGradient(controller=IC)\n", - " D = (R.adjoint*N.inverse*R + Sh.inverse).inverse\n", - " # MR FIXME: we can/should provide a preconditioner here as well!\n", - " return ift.InversionEnabler(D, inverter)\n" + " # WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy\n", + " # helper methods.\n", + " return ift.library.WienerFilterCurvature(R,N,Sh,inverter)\n" ] }, { @@ -245,7 +245,8 @@ " std=noise_amplitude, mean=0)\n", "d = noiseless_data + n\n", "j = R.adjoint_times(N.inverse_times(d))\n", - "D = PropagatorOperator(R=R, N=N, Sh=Sh)" + "curv = Curvature(R=R, N=N, Sh=Sh)\n", + "D = curv.inverse" ] }, { @@ -468,7 +469,8 @@ }, "outputs": [], "source": [ - "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", + "curv = Curvature(R=R, N=N, Sh=Sh)\n", + "D = curv.inverse\n", "j = R.adjoint_times(N.inverse_times(d))\n", "m = D(j)" ] @@ -493,12 +495,6 @@ "outputs": [], "source": [ "sc = ift.probing.utils.StatCalculator()\n", - "\n", - "IC = ift.GradientNormController(iteration_limit=50000,\n", - " tol_abs_gradnorm=0.1)\n", - "inverter = ift.ConjugateGradient(controller=IC)\n", - "curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter)\n", - "\n", "for i in range(200):\n", " print i\n", " sc.add(HT(curv.generate_posterior_sample()))\n", @@ -627,8 +623,6 @@ "# Operators\n", "Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)\n", "N = ift.ScalingOperator(sigma2,s_space)\n", - "R = ift.FFTSmoothingOperator(s_space, sigma=.01)\n", - "#D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", "\n", "# Fields and data\n", "sh = ift.power_synthesize(ift.PS_field(p_space,pow_spec),real_signal=True)\n", @@ -645,7 +639,8 @@ "\n", "R = ift.DiagonalOperator(mask)*HT\n", "n.val[l:h, l:h] = 0\n", - "D = PropagatorOperator(R=R, N=N, Sh=Sh)\n", + "curv = Curvature(R=R, N=N, Sh=Sh)\n", + "D = curv.inverse\n", "\n", "d = R(sh) + n\n", "j = R.adjoint_times(N.inverse_times(d))\n", -- GitLab