Skip to content

GitLab

Commit c92db8ac authored by Martin Reinecke's avatar Martin Reinecke
Browse files

tweaks

parent 35e632b0
Pipeline #24368 passed with stage
in 24 minutes and 15 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 36 additions and 57 deletions
demos/Wiener Filter.ipynb
View file @ c92db8ac
......@@ -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",
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment