Commit 58d4b593 authored by Philipp Arras's avatar Philipp Arras
Browse files

Rename nifty4 -> nifty5

parent 02f153a5
......@@ -34,7 +34,7 @@ test_python2_with_coverage:
script:
- python setup.py install --user -f
- mpiexec -n 2 --bind-to none nosetests -q 2> /dev/null
- nosetests -q --with-coverage --cover-package=nifty4 --cover-erase
- nosetests -q --with-coverage --cover-package=nifty5 --cover-erase
- >
coverage report --omit "*plotting*,*distributed_do*"
- >
......
......@@ -49,8 +49,8 @@ Optional dependencies:
### Sources
The current version of Nifty4 can be obtained by cloning the repository and
switching to the NIFTy_4 branch:
The current version of Nifty5 can be obtained by cloning the repository and
switching to the NIFTy_5 branch:
git clone https://gitlab.mpcdf.mpg.de/ift/NIFTy.git
......@@ -59,7 +59,7 @@ switching to the NIFTy_4 branch:
In the following, we assume a Debian-based distribution. For other
distributions, the "apt" lines will need slight changes.
NIFTy4 and its mandatory dependencies can be installed via:
NIFTy5 and its mandatory dependencies can be installed via:
sudo apt-get install git libfftw3-dev python python-pip python-dev
pip install --user git+https://gitlab.mpcdf.mpg.de/ift/NIFTy.git@NIFTy_4
......@@ -99,7 +99,7 @@ In oder to run the tests one needs two additional packages:
Afterwards the tests (including a coverage report) can be run using the
following command in the repository root:
nosetests -x --with-coverage --cover-html --cover-package=nifty4
nosetests -x --with-coverage --cover-html --cover-package=nifty5
### First Steps
......
......@@ -140,7 +140,7 @@
"source": [
"import numpy as np\n",
"np.random.seed(40)\n",
"import nifty4 as ift\n",
"import nifty5 as ift\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
......
%% Cell type:markdown id: tags:
# A NIFTy demonstration
%% Cell type:markdown id: tags:
## IFT: Big Picture
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
### Assumptions
- $d=Rs+n$, $R$ linear operator.
- $\mathcal P (s) = \mathcal G (s,S)$, $\mathcal P (n) = \mathcal G (n,N)$ where $S, N$ are positive definite matrices.
### Posterior
The Posterior is given by:
$$\mathcal P (s|d) \propto P(s,d) = \mathcal G(d-Rs,N) \,\mathcal G(s,S) \propto \mathcal G (m,D) $$
where
$$\begin{align}
m &= Dj \\
D^{-1}&= (S^{-1} +R^\dagger N^{-1} R )\\
j &= R^\dagger N^{-1} d
\end{align}$$
Let us implement this in NIFTy!
%% Cell type:markdown id: tags:
## Wiener Filter: Example
- We assume statistical homogeneity and isotropy. Therefore the signal covariance $S$ is diagonal in harmonic space, and is described by a one-dimensional power spectrum, assumed here as $$P(k) = P_0\,\left(1+\left(\frac{k}{k_0}\right)^2\right)^{-\gamma /2},$$
with $P_0 = 0.2, k_0 = 5, \gamma = 4$.
- $N = 0.2 \cdot \mathbb{1}$.
- Number of data points $N_{pix} = 512$.
- reconstruction in harmonic space.
- Response operator:
$$R = FFT_{\text{harmonic} \rightarrow \text{position}}$$
%% Cell type:code id: tags:
``` python
N_pixels = 512 # Number of pixels
def pow_spec(k):
P0, k0, gamma = [.2, 5, 4]
return P0 / ((1. + (k/k0)**2)**(gamma / 2))
```
%% Cell type:markdown id: tags:
## Wiener Filter: Implementation
%% Cell type:markdown id: tags:
### Import Modules
%% Cell type:code id: tags:
``` python
import numpy as np
np.random.seed(40)
import nifty4 as ift
import nifty5 as ift
import matplotlib.pyplot as plt
%matplotlib inline
```
%% Cell type:markdown id: tags:
### Implement Propagator
%% Cell type:code id: tags:
``` python
def Curvature(R, N, Sh):
IC = ift.GradientNormController(iteration_limit=50000,
tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=IC)
# WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy
# helper methods.
return ift.library.WienerFilterCurvature(R,N,Sh,inverter)
```
%% Cell type:markdown id: tags:
### 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*).
<!--
- One can define the *condition number* of a non-singular and normal matrix $A$:
$$\kappa (A) := \frac{|\lambda_{\text{max}}|}{|\lambda_{\text{min}}|},$$
where $\lambda_{\text{max}}$ and $\lambda_{\text{min}}$ are the largest and smallest eigenvalue of $A$, respectively.
- The larger $\kappa$ the slower Conjugate Gradient.
- 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:
$$\tilde A m = \tilde j,$$
where $\tilde A = T D^{-1}$ and $\tilde j = Tj$.
- In our case $S^{-1}$ is responsible for the bad conditioning of $D$ depending on the chosen power spectrum. Thus, we choose
$$T = \mathcal F^\dagger S_h^{-1} \mathcal F.$$
-->
%% Cell type:markdown id: tags:
### Generate Mock data
- Generate a field $s$ and $n$ with given covariances.
- Calculate $d$.
%% Cell type:code id: tags:
``` python
s_space = ift.RGSpace(N_pixels)
h_space = s_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(h_space, target=s_space)
# Operators
Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
R = HT #*ift.create_harmonic_smoothing_operator((h_space,), 0, 0.02)
# Fields and data
sh = Sh.draw_sample()
noiseless_data=R(sh)
noise_amplitude = np.sqrt(0.2)
N = ift.ScalingOperator(noise_amplitude**2, s_space)
n = ift.Field.from_random(domain=s_space, random_type='normal',
std=noise_amplitude, mean=0)
d = noiseless_data + n
j = R.adjoint_times(N.inverse_times(d))
curv = Curvature(R=R, N=N, Sh=Sh)
D = curv.inverse
```
%% Cell type:markdown id: tags:
### Run Wiener Filter
%% Cell type:code id: tags:
``` python
m = D(j)
```
%% Cell type:markdown id: tags:
### Signal Reconstruction
%% Cell type:code id: tags:
``` python
# Get signal data and reconstruction data
s_data = HT(sh).to_global_data()
m_data = HT(m).to_global_data()
d_data = d.to_global_data()
plt.figure(figsize=(15,10))
plt.plot(s_data, 'r', label="Signal", linewidth=3)
plt.plot(d_data, 'k.', label="Data")
plt.plot(m_data, 'k', label="Reconstruction",linewidth=3)
plt.title("Reconstruction")
plt.legend()
plt.show()
```
%% Cell type:code id: tags:
``` python
plt.figure(figsize=(15,10))
plt.plot(s_data - s_data, 'r', label="Signal", linewidth=3)
plt.plot(d_data - s_data, 'k.', label="Data")
plt.plot(m_data - s_data, 'k', label="Reconstruction",linewidth=3)
plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5)
plt.title("Residuals")
plt.legend()
plt.show()
```
%% Cell type:markdown id: tags:
### Power Spectrum
%% Cell type:code id: tags:
``` python
s_power_data = ift.power_analyze(sh).to_global_data()
m_power_data = ift.power_analyze(m).to_global_data()
plt.figure(figsize=(15,10))
plt.loglog()
plt.xlim(1, int(N_pixels/2))
ymin = min(m_power_data)
plt.ylim(ymin, 1)
xs = np.arange(1,int(N_pixels/2),.1)
plt.plot(xs, pow_spec(xs), label="True Power Spectrum", color='k',alpha=0.5)
plt.plot(s_power_data, 'r', label="Signal")
plt.plot(m_power_data, 'k', label="Reconstruction")
plt.axhline(noise_amplitude**2 / N_pixels, color="k", linestyle='--', label="Noise level", alpha=.5)
plt.axhspan(noise_amplitude**2 / N_pixels, ymin, facecolor='0.9', alpha=.5)
plt.title("Power Spectrum")
plt.legend()
plt.show()
```
%% Cell type:markdown id: tags:
## Wiener Filter on Incomplete Data
%% Cell type:code id: tags:
``` python
# Operators
Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
N = ift.ScalingOperator(noise_amplitude**2,s_space)
# R is defined below
# Fields
sh = Sh.draw_sample()
s = HT(sh)
n = ift.Field.from_random(domain=s_space, random_type='normal',
std=noise_amplitude, mean=0)
```
%% Cell type:markdown id: tags:
### Partially Lose Data
%% Cell type:code id: tags:
``` python
l = int(N_pixels * 0.2)
h = int(N_pixels * 0.2 * 2)
mask = np.full(s_space.shape, 1.)
mask[l:h] = 0
mask = ift.Field.from_global_data(s_space, mask)
R = ift.DiagonalOperator(mask)*HT
n = n.to_global_data()
n[l:h] = 0
n = ift.Field.from_global_data(s_space, n)
d = R(sh) + n
```
%% Cell type:code id: tags:
``` python
curv = Curvature(R=R, N=N, Sh=Sh)
D = curv.inverse
j = R.adjoint_times(N.inverse_times(d))
m = D(j)
```
%% Cell type:markdown id: tags:
### Compute Uncertainty
%% Cell type:code id: tags:
``` python
m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 200)
```
%% Cell type:markdown id: tags:
### Get data
%% Cell type:code id: tags:
``` python
# Get signal data and reconstruction data
s_data = s.to_global_data()
m_data = HT(m).to_global_data()
m_var_data = m_var.to_global_data()
uncertainty = np.sqrt(m_var_data)
d_data = d.to_global_data()
# Set lost data to NaN for proper plotting
d_data[d_data == 0] = np.nan
```
%% Cell type:code id: tags:
``` python
fig = plt.figure(figsize=(15,10))
plt.axvspan(l, h, facecolor='0.8',alpha=0.5)
plt.fill_between(range(N_pixels), m_data - uncertainty, m_data + uncertainty, facecolor='0.5', alpha=0.5)
plt.plot(s_data, 'r', label="Signal", alpha=1, linewidth=3)
plt.plot(d_data, 'k.', label="Data")
plt.plot(m_data, 'k', label="Reconstruction", linewidth=3)
plt.title("Reconstruction of incomplete data")
plt.legend()
```
%% Cell type:markdown id: tags:
# 2d Example
%% Cell type:code id: tags:
``` python
N_pixels = 256 # Number of pixels
sigma2 = 2. # Noise variance
def pow_spec(k):
P0, k0, gamma = [.2, 2, 4]
return P0 * (1. + (k/k0)**2)**(-gamma/2)
s_space = ift.RGSpace([N_pixels, N_pixels])
```
%% Cell type:code id: tags:
``` python
h_space = s_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(h_space,s_space)
# Operators
Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
N = ift.ScalingOperator(sigma2,s_space)
# Fields and data
sh = Sh.draw_sample()
n = ift.Field.from_random(domain=s_space, random_type='normal',
std=np.sqrt(sigma2), mean=0)
# Lose some data
l = int(N_pixels * 0.33)
h = int(N_pixels * 0.33 * 2)
mask = np.full(s_space.shape, 1.)
mask[l:h,l:h] = 0.
mask = ift.Field.from_global_data(s_space, mask)
R = ift.DiagonalOperator(mask)*HT
n = n.to_global_data()
n[l:h, l:h] = 0
n = ift.Field.from_global_data(s_space, n)
curv = Curvature(R=R, N=N, Sh=Sh)
D = curv.inverse
d = R(sh) + n
j = R.adjoint_times(N.inverse_times(d))
# Run Wiener filter
m = D(j)
# Uncertainty
m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 20)
# Get data
s_data = HT(sh).to_global_data()
m_data = HT(m).to_global_data()
m_var_data = m_var.to_global_data()
d_data = d.to_global_data()
uncertainty = np.sqrt(np.abs(m_var_data))
```
%% Cell type:code id: tags:
``` python
cm = ['magma', 'inferno', 'plasma', 'viridis'][1]
mi = np.min(s_data)
ma = np.max(s_data)
fig, axes = plt.subplots(1, 2, figsize=(15, 7))
data = [s_data, d_data]
caption = ["Signal", "Data"]
for ax in axes.flat:
im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi,
vmax=ma)
ax.set_title(caption.pop(0))
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
fig.colorbar(im, cax=cbar_ax)
```
%% Cell type:code id: tags:
``` python
mi = np.min(s_data)
ma = np.max(s_data)
fig, axes = plt.subplots(3, 2, figsize=(15, 22.5))
sample = HT(curv.draw_sample(from_inverse=True)+m).to_global_data()
post_mean = (m_mean + HT(m)).to_global_data()
data = [s_data, m_data, post_mean, sample, s_data - m_data, uncertainty]
caption = ["Signal", "Reconstruction", "Posterior mean", "Sample", "Residuals", "Uncertainty Map"]
for ax in axes.flat:
im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi, vmax=ma)
ax.set_title(caption.pop(0))
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([.85, 0.15, 0.05, 0.7])
fig.colorbar(im, cax=cbar_ax)
```
%% Cell type:markdown id: tags:
### Is the uncertainty map reliable?
%% Cell type:code id: tags:
``` python
precise = (np.abs(s_data-m_data) < uncertainty)
print("Error within uncertainty map bounds: " + str(np.sum(precise) * 100 / N_pixels**2) + "%")
plt.figure(figsize=(15,10))
plt.imshow(precise.astype(float), cmap="brg")
plt.colorbar()
```
%% Cell type:markdown id: tags:
# Start Coding
## NIFTy Repository + Installation guide
https://gitlab.mpcdf.mpg.de/ift/NIFTy
NIFTy v4 **more or less stable!**
......
import nifty4 as ift
from nifty4.library.nonlinearities import Linear, Tanh, Exponential
import nifty5 as ift
from nifty5.library.nonlinearities import Linear, Tanh, Exponential
import numpy as np
np.random.seed(42)
......
import nifty4 as ift
from nifty4.library.nonlinearities import Exponential
import nifty5 as ift
from nifty5.library.nonlinearities import Exponential
import numpy as np
np.random.seed(42)
......
import nifty4 as ift
from nifty4.library.nonlinearities import Linear, Exponential, Tanh
import nifty5 as ift
from nifty5.library.nonlinearities import Linear, Exponential, Tanh
import numpy as np
np.random.seed(20)
......
import numpy as np
import nifty4 as ift
import nifty5 as ift
if __name__ == "__main__":
signal_to_noise = 0.5 # The signal to noise ratio
......
import nifty4 as ift
import nifty5 as ift
import numpy as np
......
......@@ -4,7 +4,7 @@
# to get exact figure of paper comment out lines marked by "COMMENT OUT"
import numpy as np
import nifty4 as ift
import nifty5 as ift
import matplotlib.pyplot as plt
......
import nifty4 as ift
import nifty5 as ift
import numpy as np
......
import nifty4 as ift
import nifty5 as ift
import numpy as np
import matplotlib.pyplot as plt
from nifty4.sugar import create_power_operator
from nifty5.sugar import create_power_operator
np.random.seed(42)
......
import nifty4 as ift
import nifty5 as ift
import numpy as np
......
import numpy as np
import nifty4 as ift
import nifty5 as ift
# TODO: MAKE RESPONSE MPI COMPATIBLE OR USE LOS RESPONSE INSTEAD
......
import numpy as np
import nifty4 as ift
import nifty5 as ift
if __name__ == "__main__":
......
import numpy as np
import nifty4 as ift
import nifty5 as ift
import numericalunits as nu
......
import nifty4 as ift
import nifty5 as ift
import numpy as np
np.random.seed(42)
......
rm -rf docs/build docs/source/mod
sphinx-apidoc -l -e -d 2 -o docs/source/mod nifty4
sphinx-apidoc -l -e -d 2 -o docs/source/mod nifty5
sphinx-build -b html docs/source/ docs/build/
......@@ -7,7 +7,7 @@ NIFTy-related publications
author = {{Martin Reinecke, Theo Steininger, Marco Selig}},
title = {NIFTy -- Numerical Information Field TheorY},
url = {https://gitlab.mpcdf.mpg.de/ift/NIFTy},
version = {nifty4},
version = {nifty5},
date = {2018-04-05},
}
......
.. currentmodule:: nifty4
.. currentmodule:: nifty5
=============
Code Overview
......@@ -35,7 +35,7 @@ Domains
Abstract base class
-------------------
One of the fundamental building blocks of the NIFTy4 framework is the *domain*.
One of the fundamental building blocks of the NIFTy5 framework is the *domain*.
Its required capabilities are expressed by the abstract :class:`Domain` class.
A domain must be able to answer the following queries:
......@@ -151,7 +151,7 @@ performance benefits.
Linear Operators
================
A linear operator (represented by NIFTy4's abstract :class:`LinearOperator`
A linear operator (represented by NIFTy5's abstract :class:`LinearOperator`
class) can be interpreted as an (implicitly defined) matrix.
It can be applied to :class:`Field` instances, resulting in other :class:`Field`
instances that potentially live on other domains.
......@@ -245,7 +245,7 @@ high-dimensional functions, which are often nonlinear.
Energy functionals
------------------
In NIFTy4 such functions are represented by objects of type :class:`Energy`.
In NIFTy5 such functions are represented by objects of type :class:`Energy`.
These hold the prescription how to calculate the function's
:attr:`~Energy.value`, :attr:`~Energy.gradient` and
(optionally) :attr:`~Energy.curvature` at any given :attr:`~Energy.position`
......@@ -256,9 +256,9 @@ linear operator objects.
Energies are classes that typically have to be provided by the user when
tackling new IFT problems.
Some examples of concrete energy classes delivered with NIFTy4 are
Some examples of concrete energy classes delivered with NIFTy5 are
:class:`QuadraticEnergy` (with position-independent curvature, mainly used with
conjugate gradient minimization) and :class:`~nifty4.library.WienerFilterEnergy`.
conjugate gradient minimization) and :class:`~nifty5.library.WienerFilterEnergy`.
Iteration control
......@@ -269,7 +269,7 @@ checking the quality of the current solution estimate and stopping once
it is sufficiently accurate. In case of numerical problems, the iteration needs
to be terminated as well, returning a suitable error description.
In NIFTy4, this functionality is encapsulated in the abstract
In NIFTy5, this functionality is encapsulated in the abstract
:class:`IterationController` class, which is provided with the initial energy
object before starting the minimization, and is updated with the improved
energy after every iteration. Based on this information, it can either continue
......
......@@ -73,7 +73,7 @@ source_suffix = '.rst'
master_doc = 'index'
# General information about the project.
project = u'NIFTy4'
project = u'NIFTy5'
copyright = u'2013-2018, Max-Planck-Society'
author = u'Theo Steininger / Martin Reinecke'
......
Supports Markdown
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