Commit 37ec065b authored by Philipp Arras's avatar Philipp Arras
Browse files

Restructure ipynb

parent 11683784
......@@ -2,33 +2,30 @@
# 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.
......@@ -37,21 +34,19 @@
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$.
......@@ -70,15 +65,15 @@
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
......@@ -87,11 +82,11 @@
%matplotlib inline
```
%% Cell type:markdown id: tags:
### Implement Propagator
#### Implement Propagator
%% Cell type:code id: tags:
``` python
def Curvature(R, N, Sh):
......@@ -102,11 +97,11 @@
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*).
......@@ -126,11 +121,11 @@
$$T = \mathcal F^\dagger S_h^{-1} \mathcal F.$$
-->
%% 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:
......@@ -158,21 +153,21 @@
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
......@@ -202,11 +197,11 @@
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
......@@ -319,11 +314,11 @@
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
......
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