getting_started_0.ipynb 17.6 KB
Newer Older
Philipp Arras's avatar
Philipp Arras committed
1
2
3
4
5
6
7
8
9
10
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
11
    "# Code example: Wiener filter"
Philipp Arras's avatar
Philipp Arras committed
12
13
14
15
16
17
18
19
20
21
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
22
    "## Introduction\n",
Philipp Arras's avatar
Philipp Arras committed
23
24
25
26
    "IFT starting point:\n",
    "\n",
    "$$d = Rs+n$$\n",
    "\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
27
    "Typically, $s$ is a continuous field, $d$ a discrete data vector. Particularly, $R$ is not invertible.\n",
Philipp Arras's avatar
Philipp Arras committed
28
29
30
    "\n",
    "IFT aims at **inverting** the above uninvertible problem in the **best possible way** using Bayesian statistics.\n",
    "\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
31
    "NIFTy (Numerical Information Field Theory) is a Python framework in which IFT problems can be tackled easily.\n",
Philipp Arras's avatar
Philipp Arras committed
32
33
34
35
36
    "\n",
    "Main Interfaces:\n",
    "\n",
    "- **Spaces**: Cartesian, 2-Spheres (Healpix, Gauss-Legendre) and their respective harmonic spaces.\n",
    "- **Fields**: Defined on spaces.\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
37
    "- **Operators**: Acting on fields."
Philipp Arras's avatar
Philipp Arras committed
38
39
40
41
42
43
44
45
46
47
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
48
    "## Wiener filter on one-dimensional fields\n",
Philipp Arras's avatar
Philipp Arras committed
49
50
51
52
53
54
55
56
57
    "\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",
Philipp Arras's avatar
Philipp Arras committed
58
    "$$\\mathcal P (s|d) \\propto P(s,d) = \\mathcal G(d-Rs,N) \\,\\mathcal G(s,S) \\propto \\mathcal G (s-m,D) $$\n",
Philipp Arras's avatar
Philipp Arras committed
59
60
    "\n",
    "where\n",
Philipp Arras's avatar
Philipp Arras committed
61
62
63
    "$$m = Dj$$\n",
    "with\n",
    "$$D = (S^{-1} +R^\\dagger N^{-1} R)^{-1} , \\quad j = R^\\dagger N^{-1} d.$$\n",
Philipp Arras's avatar
Philipp Arras committed
64
65
66
67
68
69
70
71
72
73
74
75
    "\n",
    "Let us implement this in NIFTy!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
76
    "### In NIFTy\n",
Philipp Arras's avatar
Philipp Arras committed
77
    "\n",
Martin Reinecke's avatar
Martin Reinecke committed
78
79
    "- 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",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
80
    "- $N = 0.2 \\cdot \\mathbb{1}$.\n",
Martin Reinecke's avatar
Martin Reinecke committed
81
82
    "- Number of data points $N_{pix} = 512$.\n",
    "- reconstruction in harmonic space.\n",
Philipp Arras's avatar
Philipp Arras committed
83
    "- Response operator:\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
84
    "$$R = FFT_{\\text{harmonic} \\rightarrow \\text{position}}$$\n"
Philipp Arras's avatar
Philipp Arras committed
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N_pixels = 512     # Number of pixels\n",
    "\n",
    "def pow_spec(k):\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
100
101
    "    P0, k0, gamma = [.2, 5, 4]\n",
    "    return P0 / ((1. + (k/k0)**2)**(gamma / 2))"
Philipp Arras's avatar
Philipp Arras committed
102
103
104
105
106
107
108
109
110
111
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
112
    "### Implementation"
Philipp Arras's avatar
Philipp Arras committed
113
114
115
116
117
118
119
120
121
122
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
123
    "#### Import Modules"
Philipp Arras's avatar
Philipp Arras committed
124
125
126
127
128
129
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
Martin Reinecke's avatar
Martin Reinecke committed
130
    "scrolled": true,
Philipp Arras's avatar
Philipp Arras committed
131
132
133
134
135
136
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
Philipp Arras's avatar
Philipp Arras committed
137
    "%matplotlib inline\n",
Philipp Arras's avatar
Philipp Arras committed
138
    "import numpy as np\n",
Philipp Arras's avatar
Philipp Arras committed
139
    "import nifty8 as ift\n",
140
    "import matplotlib.pyplot as plt\n",
Martin Reinecke's avatar
Martin Reinecke committed
141
    "plt.rcParams['figure.dpi'] = 100\n",
Philipp Arras's avatar
Philipp Arras committed
142
    "plt.style.use(\"seaborn-notebook\")"
Philipp Arras's avatar
Philipp Arras committed
143
144
145
146
147
148
149
150
151
152
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
153
    "#### Implement Propagator"
Philipp Arras's avatar
Philipp Arras committed
154
155
156
157
158
159
160
161
162
163
164
165
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
166
    "def Curvature(R, N, Sh):\n",
Martin Reinecke's avatar
Martin Reinecke committed
167
    "    IC = ift.GradientNormController(iteration_limit=50000,\n",
168
    "                                    tol_abs_gradnorm=0.1)\n",
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
169
170
    "    # WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy\n",
    "    # helper methods.\n",
171
    "    return ift.WienerFilterCurvature(R,N,Sh,iteration_controller=IC,iteration_controller_sampling=IC)"
Philipp Arras's avatar
Philipp Arras committed
172
173
174
175
176
177
178
179
180
181
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
182
    "#### Conjugate Gradient Preconditioning\n",
Philipp Arras's avatar
Philipp Arras committed
183
184
    "\n",
    "- $D$ is defined via:\n",
Martin Reinecke's avatar
Martin Reinecke committed
185
    "$$D^{-1} = \\mathcal S_h^{-1} + R^\\dagger N^{-1} R.$$\n",
Philipp Arras's avatar
Philipp Arras committed
186
187
    "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",
Martin Reinecke's avatar
Martin Reinecke committed
188
    "<!--\n",
Philipp Arras's avatar
Philipp Arras committed
189
190
191
192
193
194
    "- 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",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
195
    "- 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",
Philipp Arras's avatar
Philipp Arras committed
196
197
198
199
200
    "$$\\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",
Martin Reinecke's avatar
Martin Reinecke committed
201
202
    "$$T = \\mathcal F^\\dagger S_h^{-1} \\mathcal F.$$\n",
    "-->"
Philipp Arras's avatar
Philipp Arras committed
203
204
205
206
207
208
209
210
211
212
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
213
    "#### Generate Mock data\n",
Philipp Arras's avatar
Philipp Arras committed
214
215
216
217
218
219
220
221
    "\n",
    "- Generate a field $s$ and $n$ with given covariances.\n",
    "- Calculate $d$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
Martin Reinecke's avatar
Martin Reinecke committed
222
223
224
   "metadata": {
    "scrolled": true
   },
Philipp Arras's avatar
Philipp Arras committed
225
226
   "outputs": [],
   "source": [
227
228
229
    "s_space = ift.RGSpace(N_pixels)\n",
    "h_space = s_space.get_default_codomain()\n",
    "HT = ift.HarmonicTransformOperator(h_space, target=s_space)\n",
Philipp Arras's avatar
Philipp Arras committed
230
231
    "\n",
    "# Operators\n",
232
    "Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)\n",
Philipp Arras's avatar
Philipp Arras committed
233
    "R = HT # @ ift.create_harmonic_smoothing_operator((h_space,), 0, 0.02)\n",
Philipp Arras's avatar
Philipp Arras committed
234
235
    "\n",
    "# Fields and data\n",
Philipp Arras's avatar
Philipp Arras committed
236
    "sh = Sh.draw_sample_with_dtype(dtype=np.float64)\n",
237
    "noiseless_data=R(sh)\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
238
    "noise_amplitude = np.sqrt(0.2)\n",
239
    "N = ift.ScalingOperator(s_space, noise_amplitude**2)\n",
240
241
    "\n",
    "n = ift.Field.from_random(domain=s_space, random_type='normal',\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
242
    "                          std=noise_amplitude, mean=0)\n",
243
244
    "d = noiseless_data + n\n",
    "j = R.adjoint_times(N.inverse_times(d))\n",
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
245
246
    "curv = Curvature(R=R, N=N, Sh=Sh)\n",
    "D = curv.inverse"
Philipp Arras's avatar
Philipp Arras committed
247
248
249
250
251
252
253
254
255
256
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
257
    "#### Run Wiener Filter"
Philipp Arras's avatar
Philipp Arras committed
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "m = D(j)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
281
    "#### Results"
Philipp Arras's avatar
Philipp Arras committed
282
283
284
285
286
287
288
289
290
291
292
293
294
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Get signal data and reconstruction data\n",
Martin Reinecke's avatar
Martin Reinecke committed
295
296
297
    "s_data = HT(sh).val\n",
    "m_data = HT(m).val\n",
    "d_data = d.val\n",
Philipp Arras's avatar
Philipp Arras committed
298
    "\n",
299
    "plt.plot(s_data, 'r', label=\"Signal\", linewidth=2)\n",
Martin Reinecke's avatar
Martin Reinecke committed
300
    "plt.plot(d_data, 'k.', label=\"Data\")\n",
301
    "plt.plot(m_data, 'k', label=\"Reconstruction\",linewidth=2)\n",
Philipp Arras's avatar
Philipp Arras committed
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
    "plt.title(\"Reconstruction\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
317
    "plt.plot(s_data - s_data, 'r', label=\"Signal\", linewidth=2)\n",
Martin Reinecke's avatar
Martin Reinecke committed
318
    "plt.plot(d_data - s_data, 'k.', label=\"Data\")\n",
319
    "plt.plot(m_data - s_data, 'k', label=\"Reconstruction\",linewidth=2)\n",
320
    "plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5)\n",
Philipp Arras's avatar
Philipp Arras committed
321
322
323
324
325
326
327
328
329
330
331
332
333
    "plt.title(\"Residuals\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
334
    "#### Power Spectrum"
Philipp Arras's avatar
Philipp Arras committed
335
336
337
338
339
340
341
342
343
344
345
346
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
Martin Reinecke's avatar
Martin Reinecke committed
347
348
    "s_power_data = ift.power_analyze(sh).val\n",
    "m_power_data = ift.power_analyze(m).val\n",
Philipp Arras's avatar
Philipp Arras committed
349
350
351
352
353
    "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",
Martin Reinecke's avatar
Martin Reinecke committed
354
355
356
    "plt.plot(xs, pow_spec(xs), label=\"True Power Spectrum\", color='k',alpha=0.5)\n",
    "plt.plot(s_power_data, 'r', label=\"Signal\")\n",
    "plt.plot(m_power_data, 'k', label=\"Reconstruction\")\n",
357
358
    "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",
Philipp Arras's avatar
Philipp Arras committed
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
    "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": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Operators\n",
386
    "Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)\n",
387
    "N = ift.ScalingOperator(s_space, noise_amplitude**2)\n",
Philipp Arras's avatar
Philipp Arras committed
388
389
390
    "# R is defined below\n",
    "\n",
    "# Fields\n",
Philipp Arras's avatar
Philipp Arras committed
391
    "sh = Sh.draw_sample_with_dtype(dtype=np.float64)\n",
392
393
394
    "s = HT(sh)\n",
    "n = ift.Field.from_random(domain=s_space, random_type='normal',\n",
    "                      std=noise_amplitude, mean=0)"
Philipp Arras's avatar
Philipp Arras committed
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Partially Lose Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "l = int(N_pixels * 0.2)\n",
419
    "h = int(N_pixels * 0.2 * 2)\n",
Philipp Arras's avatar
Philipp Arras committed
420
    "\n",
Martin Reinecke's avatar
Martin Reinecke committed
421
422
    "mask = np.full(s_space.shape, 1.)\n",
    "mask[l:h] = 0\n",
Martin Reinecke's avatar
Martin Reinecke committed
423
    "mask = ift.Field.from_raw(s_space, mask)\n",
Philipp Arras's avatar
Philipp Arras committed
424
    "\n",
Martin Reinecke's avatar
Martin Reinecke committed
425
    "R = ift.DiagonalOperator(mask)(HT)\n",
Martin Reinecke's avatar
Martin Reinecke committed
426
    "n = n.val_rw()\n",
Martin Reinecke's avatar
Martin Reinecke committed
427
    "n[l:h] = 0\n",
Martin Reinecke's avatar
Martin Reinecke committed
428
    "n = ift.Field.from_raw(s_space, n)\n",
Philipp Arras's avatar
Philipp Arras committed
429
    "\n",
430
    "d = R(sh) + n"
Philipp Arras's avatar
Philipp Arras committed
431
432
433
434
435
436
437
438
439
440
441
442
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
443
444
    "curv = Curvature(R=R, N=N, Sh=Sh)\n",
    "D = curv.inverse\n",
Philipp Arras's avatar
Philipp Arras committed
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
    "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": {
464
    "scrolled": true
Philipp Arras's avatar
Philipp Arras committed
465
466
467
   },
   "outputs": [],
   "source": [
Martin Reinecke's avatar
Martin Reinecke committed
468
    "m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 200, np.float64)"
Philipp Arras's avatar
Philipp Arras committed
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Get data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Get signal data and reconstruction data\n",
Martin Reinecke's avatar
Martin Reinecke committed
493
494
495
    "s_data = s.val\n",
    "m_data = HT(m).val\n",
    "m_var_data = m_var.val\n",
Martin Reinecke's avatar
Martin Reinecke committed
496
    "uncertainty = np.sqrt(m_var_data)\n",
Martin Reinecke's avatar
Martin Reinecke committed
497
    "d_data = d.val_rw()\n",
Philipp Arras's avatar
Philipp Arras committed
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
    "\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": [
Martin Reinecke's avatar
Martin Reinecke committed
513
514
    "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",
515
    "plt.plot(s_data, 'r', label=\"Signal\", alpha=1, linewidth=2)\n",
Martin Reinecke's avatar
Martin Reinecke committed
516
    "plt.plot(d_data, 'k.', label=\"Data\")\n",
517
    "plt.plot(m_data, 'k', label=\"Reconstruction\", linewidth=2)\n",
Philipp Arras's avatar
Philipp Arras committed
518
519
520
521
522
523
524
525
526
527
528
529
    "plt.title(\"Reconstruction of incomplete data\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
Philipp Arras's avatar
Philipp Arras committed
530
    "## Wiener filter on two-dimensional field"
Philipp Arras's avatar
Philipp Arras committed
531
532
533
534
535
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
Martin Reinecke's avatar
Martin Reinecke committed
536
   "metadata": {},
Philipp Arras's avatar
Philipp Arras committed
537
538
539
   "outputs": [],
   "source": [
    "N_pixels = 256      # Number of pixels\n",
Martin Reinecke's avatar
Martin Reinecke committed
540
    "sigma2 = 2.         # Noise variance\n",
Philipp Arras's avatar
Philipp Arras committed
541
542
    "\n",
    "def pow_spec(k):\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
543
    "    P0, k0, gamma = [.2, 2, 4]\n",
Martin Reinecke's avatar
Martin Reinecke committed
544
    "    return P0 * (1. + (k/k0)**2)**(-gamma/2)\n",
Philipp Arras's avatar
Philipp Arras committed
545
    "\n",
546
    "s_space = ift.RGSpace([N_pixels, N_pixels])"
Philipp Arras's avatar
Philipp Arras committed
547
548
549
550
551
552
553
554
555
556
557
558
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
559
    "h_space = s_space.get_default_codomain()\n",
Martin Reinecke's avatar
Martin Reinecke committed
560
    "HT = ift.HarmonicTransformOperator(h_space,s_space)\n",
Philipp Arras's avatar
Philipp Arras committed
561
562
    "\n",
    "# Operators\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
563
    "Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)\n",
564
    "N = ift.ScalingOperator(s_space, sigma2)\n",
Philipp Arras's avatar
Philipp Arras committed
565
566
    "\n",
    "# Fields and data\n",
Philipp Arras's avatar
Philipp Arras committed
567
    "sh = Sh.draw_sample_with_dtype(dtype=np.float64)\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
568
    "n = ift.Field.from_random(domain=s_space, random_type='normal',\n",
Philipp Arras's avatar
Philipp Arras committed
569
570
571
572
    "                      std=np.sqrt(sigma2), mean=0)\n",
    "\n",
    "# Lose some data\n",
    "\n",
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
573
574
    "l = int(N_pixels * 0.33)\n",
    "h = int(N_pixels * 0.33 * 2)\n",
Philipp Arras's avatar
Philipp Arras committed
575
    "\n",
Martin Reinecke's avatar
Martin Reinecke committed
576
577
    "mask = np.full(s_space.shape, 1.)\n",
    "mask[l:h,l:h] = 0.\n",
Martin Reinecke's avatar
Martin Reinecke committed
578
    "mask = ift.Field.from_raw(s_space, mask)\n",
Philipp Arras's avatar
Philipp Arras committed
579
    "\n",
Martin Reinecke's avatar
Martin Reinecke committed
580
    "R = ift.DiagonalOperator(mask)(HT)\n",
Martin Reinecke's avatar
Martin Reinecke committed
581
    "n = n.val_rw()\n",
Martin Reinecke's avatar
Martin Reinecke committed
582
    "n[l:h, l:h] = 0\n",
Martin Reinecke's avatar
Martin Reinecke committed
583
    "n = ift.Field.from_raw(s_space, n)\n",
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
584
585
    "curv = Curvature(R=R, N=N, Sh=Sh)\n",
    "D = curv.inverse\n",
Philipp Arras's avatar
Philipp Arras committed
586
    "\n",
Martin Reinecke's avatar
Martin Reinecke committed
587
    "d = R(sh) + n\n",
Philipp Arras's avatar
Philipp Arras committed
588
589
590
591
592
593
    "j = R.adjoint_times(N.inverse_times(d))\n",
    "\n",
    "# Run Wiener filter\n",
    "m = D(j)\n",
    "\n",
    "# Uncertainty\n",
Martin Reinecke's avatar
Martin Reinecke committed
594
    "m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 20, np.float64)\n",
Philipp Arras's avatar
Philipp Arras committed
595
596
    "\n",
    "# Get data\n",
Martin Reinecke's avatar
Martin Reinecke committed
597
598
599
600
    "s_data = HT(sh).val\n",
    "m_data = HT(m).val\n",
    "m_var_data = m_var.val\n",
    "d_data = d.val\n",
Philipp Arras's avatar
Philipp Arras committed
601
602
603
604
605
606
607
608
609
610
611
612
613
    "uncertainty = np.sqrt(np.abs(m_var_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
614
    "cmap = ['magma', 'inferno', 'plasma', 'viridis'][1]\n",
Philipp Arras's avatar
Philipp Arras committed
615
616
617
618
    "\n",
    "mi = np.min(s_data)\n",
    "ma = np.max(s_data)\n",
    "\n",
Philipp Arras's avatar
Philipp Arras committed
619
    "fig, axes = plt.subplots(1, 2)\n",
Philipp Arras's avatar
Philipp Arras committed
620
621
622
623
624
    "\n",
    "data = [s_data, d_data]\n",
    "caption = [\"Signal\", \"Data\"]\n",
    "\n",
    "for ax in axes.flat:\n",
625
    "    im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cmap, vmin=mi,\n",
Philipp Arras's avatar
Philipp Arras committed
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
    "                   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": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "mi = np.min(s_data)\n",
    "ma = np.max(s_data)\n",
    "\n",
Philipp Arras's avatar
Philipp Arras committed
647
    "fig, axes = plt.subplots(3, 2, figsize=(10, 15))\n",
Philipp Arras's avatar
Philipp Arras committed
648
    "sample = HT(curv.draw_sample(from_inverse=True)+m).val\n",
Martin Reinecke's avatar
Martin Reinecke committed
649
    "post_mean = (m_mean + HT(m)).val\n",
Philipp Arras's avatar
Philipp Arras committed
650
    "\n",
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
651
652
    "data = [s_data, m_data, post_mean, sample, s_data - m_data, uncertainty]\n",
    "caption = [\"Signal\", \"Reconstruction\", \"Posterior mean\", \"Sample\", \"Residuals\", \"Uncertainty Map\"]\n",
Philipp Arras's avatar
Philipp Arras committed
653
654
    "\n",
    "for ax in axes.flat:\n",
655
    "    im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cmap, vmin=mi, vmax=ma)\n",
Philipp Arras's avatar
Philipp Arras committed
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
    "    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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Is the uncertainty map reliable?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
Martin Reinecke's avatar
Martin Reinecke committed
684
    "precise = (np.abs(s_data-m_data) < uncertainty)\n",
Philipp Arras's avatar
Philipp Arras committed
685
686
687
    "print(\"Error within uncertainty map bounds: \" + str(np.sum(precise) * 100 / N_pixels**2) + \"%\")\n",
    "\n",
    "plt.imshow(precise.astype(float), cmap=\"brg\")\n",
Martin Reinecke's avatar
Martin Reinecke committed
688
    "plt.colorbar()"
Philipp Arras's avatar
Philipp Arras committed
689
690
691
692
693
694
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
695
   "display_name": "Python 3",
Philipp Arras's avatar
Philipp Arras committed
696
   "language": "python",
697
   "name": "python3"
Philipp Arras's avatar
Philipp Arras committed
698
699
700
701
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
702
    "version": 3
Philipp Arras's avatar
Philipp Arras committed
703
704
705
706
707
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
708
   "pygments_lexer": "ipython3",
709
   "version": "3.9.2"
Philipp Arras's avatar
Philipp Arras committed
710
711
712
713
714
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}