Commit 1762da4a authored by Martin Reinecke's avatar Martin Reinecke
Browse files

changes for NIFTy_7, round 1

parent fe178ade
...@@ -37,7 +37,7 @@ build_docker_from_cache: ...@@ -37,7 +37,7 @@ build_docker_from_cache:
test_serial: test_serial:
stage: test stage: test
script: script:
- pytest-3 -q --cov=nifty6 test - pytest-3 -q --cov=nifty7 test
- > - >
python3 -m coverage report --omit "*plot*" | tee coverage.txt python3 -m coverage report --omit "*plot*" | tee coverage.txt
- > - >
...@@ -59,7 +59,7 @@ pages: ...@@ -59,7 +59,7 @@ pages:
paths: paths:
- public - public
only: only:
- NIFTy_5 - NIFTy_6
before_script: before_script:
......
NIFTy - Numerical Information Field Theory NIFTy - Numerical Information Field Theory
========================================== ==========================================
[![build status](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/NIFTy_6/build.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/NIFTy_6) [![build status](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/NIFTy_7/build.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/NIFTy_7)
[![coverage report](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/NIFTy_6/coverage.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/NIFTy_6) [![coverage report](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/NIFTy_7/coverage.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/NIFTy_7)
**NIFTy** project homepage: **NIFTy** project homepage:
[http://ift.pages.mpcdf.de/nifty](http://ift.pages.mpcdf.de/nifty) [http://ift.pages.mpcdf.de/nifty](http://ift.pages.mpcdf.de/nifty)
...@@ -59,8 +59,8 @@ Optional dependencies: ...@@ -59,8 +59,8 @@ Optional dependencies:
### Sources ### Sources
The current version of NIFTy6 can be obtained by cloning the repository and The current version of NIFTy7 can be obtained by cloning the repository and
switching to the NIFTy_6 branch: switching to the NIFTy_7 branch:
git clone https://gitlab.mpcdf.mpg.de/ift/nifty.git git clone https://gitlab.mpcdf.mpg.de/ift/nifty.git
...@@ -69,10 +69,10 @@ switching to the NIFTy_6 branch: ...@@ -69,10 +69,10 @@ switching to the NIFTy_6 branch:
In the following, we assume a Debian-based distribution. For other In the following, we assume a Debian-based distribution. For other
distributions, the "apt" lines will need slight changes. distributions, the "apt" lines will need slight changes.
NIFTy6 and its mandatory dependencies can be installed via: NIFTy7 and its mandatory dependencies can be installed via:
sudo apt-get install git python3 python3-pip python3-dev sudo apt-get install git python3 python3-pip python3-dev
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_6 pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/nifty.git@NIFTy_7
Plotting support is added via: Plotting support is added via:
...@@ -107,7 +107,7 @@ To run the tests, additional packages are required: ...@@ -107,7 +107,7 @@ To run the tests, additional packages are required:
Afterwards the tests (including a coverage report) can be run using the Afterwards the tests (including a coverage report) can be run using the
following command in the repository root: following command in the repository root:
pytest-3 --cov=nifty6 test pytest-3 --cov=nifty7 test
### First Steps ### First Steps
......
...@@ -3,7 +3,7 @@ from time import time ...@@ -3,7 +3,7 @@ from time import time
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
import nifty6 as ift import nifty7 as ift
N0s, a0s, b0s, c0s = [], [], [], [] N0s, a0s, b0s, c0s = [], [], [], []
......
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
import numpy as np import numpy as np
import nifty6 as ift import nifty7 as ift
if __name__ == '__main__': if __name__ == '__main__':
# Set up the position space of the signal # Set up the position space of the signal
......
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# A NIFTy demonstration # A NIFTy demonstration
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## IFT: Big Picture ## IFT: Big Picture
IFT starting point: IFT starting point:
$$d = Rs+n$$ $$d = Rs+n$$
Typically, $s$ is a continuous field, $d$ a discrete data vector. Particularly, $R$ is not invertible. 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. IFT aims at **inverting** the above uninvertible problem in the **best possible way** using Bayesian statistics.
## NIFTy ## NIFTy
NIFTy (Numerical Information Field Theory) is a Python framework in which IFT problems can be tackled easily. NIFTy (Numerical Information Field Theory) is a Python framework in which IFT problems can be tackled easily.
Main Interfaces: Main Interfaces:
- **Spaces**: Cartesian, 2-Spheres (Healpix, Gauss-Legendre) and their respective harmonic spaces. - **Spaces**: Cartesian, 2-Spheres (Healpix, Gauss-Legendre) and their respective harmonic spaces.
- **Fields**: Defined on spaces. - **Fields**: Defined on spaces.
- **Operators**: Acting on fields. - **Operators**: Acting on fields.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Wiener Filter: Formulae ## Wiener Filter: Formulae
### Assumptions ### Assumptions
- $d=Rs+n$, $R$ linear operator. - $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. - $\mathcal P (s) = \mathcal G (s,S)$, $\mathcal P (n) = \mathcal G (n,N)$ where $S, N$ are positive definite matrices.
### Posterior ### Posterior
The Posterior is given by: 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) $$ $$\mathcal P (s|d) \propto P(s,d) = \mathcal G(d-Rs,N) \,\mathcal G(s,S) \propto \mathcal G (s-m,D) $$
where where
$$\begin{align} $$\begin{align}
m &= Dj \\ m &= Dj \\
D^{-1}&= (S^{-1} +R^\dagger N^{-1} R )\\ D^{-1}&= (S^{-1} +R^\dagger N^{-1} R )\\
j &= R^\dagger N^{-1} d j &= R^\dagger N^{-1} d
\end{align}$$ \end{align}$$
Let us implement this in NIFTy! Let us implement this in NIFTy!
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Wiener Filter: Example ## 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},$$ - 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$. with $P_0 = 0.2, k_0 = 5, \gamma = 4$.
- $N = 0.2 \cdot \mathbb{1}$. - $N = 0.2 \cdot \mathbb{1}$.
- Number of data points $N_{pix} = 512$. - Number of data points $N_{pix} = 512$.
- reconstruction in harmonic space. - reconstruction in harmonic space.
- Response operator: - Response operator:
$$R = FFT_{\text{harmonic} \rightarrow \text{position}}$$ $$R = FFT_{\text{harmonic} \rightarrow \text{position}}$$
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
N_pixels = 512 # Number of pixels N_pixels = 512 # Number of pixels
def pow_spec(k): def pow_spec(k):
P0, k0, gamma = [.2, 5, 4] P0, k0, gamma = [.2, 5, 4]
return P0 / ((1. + (k/k0)**2)**(gamma / 2)) return P0 / ((1. + (k/k0)**2)**(gamma / 2))
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Wiener Filter: Implementation ## Wiener Filter: Implementation
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Import Modules ### Import Modules
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import numpy as np import numpy as np
import nifty6 as ift import nifty7 as ift
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
%matplotlib inline %matplotlib inline
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Implement Propagator ### Implement Propagator
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def Curvature(R, N, Sh): def Curvature(R, N, Sh):
IC = ift.GradientNormController(iteration_limit=50000, IC = ift.GradientNormController(iteration_limit=50000,
tol_abs_gradnorm=0.1) tol_abs_gradnorm=0.1)
# WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy # WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy
# helper methods. # helper methods.
return ift.WienerFilterCurvature(R,N,Sh,iteration_controller=IC,iteration_controller_sampling=IC) return ift.WienerFilterCurvature(R,N,Sh,iteration_controller=IC,iteration_controller_sampling=IC)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Conjugate Gradient Preconditioning ### Conjugate Gradient Preconditioning
- $D$ is defined via: - $D$ is defined via:
$$D^{-1} = \mathcal S_h^{-1} + R^\dagger N^{-1} R.$$ $$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*). 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$: - One can define the *condition number* of a non-singular and normal matrix $A$:
$$\kappa (A) := \frac{|\lambda_{\text{max}}|}{|\lambda_{\text{min}}|},$$ $$\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. where $\lambda_{\text{max}}$ and $\lambda_{\text{min}}$ are the largest and smallest eigenvalue of $A$, respectively.
- The larger $\kappa$ the slower Conjugate Gradient. - 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: - 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,$$ $$\tilde A m = \tilde j,$$
where $\tilde A = T D^{-1}$ and $\tilde j = Tj$. 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 - 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.$$ $$T = \mathcal F^\dagger S_h^{-1} \mathcal F.$$
--> -->
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Generate Mock data ### Generate Mock data
- Generate a field $s$ and $n$ with given covariances. - Generate a field $s$ and $n$ with given covariances.
- Calculate $d$. - Calculate $d$.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
s_space = ift.RGSpace(N_pixels) s_space = ift.RGSpace(N_pixels)
h_space = s_space.get_default_codomain() h_space = s_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(h_space, target=s_space) HT = ift.HarmonicTransformOperator(h_space, target=s_space)
# Operators # Operators
Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec) Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
R = HT #*ift.create_harmonic_smoothing_operator((h_space,), 0, 0.02) R = HT #*ift.create_harmonic_smoothing_operator((h_space,), 0, 0.02)
# Fields and data # Fields and data
sh = Sh.draw_sample_with_dtype(dtype=np.float64) sh = Sh.draw_sample_with_dtype(dtype=np.float64)
noiseless_data=R(sh) noiseless_data=R(sh)
noise_amplitude = np.sqrt(0.2) noise_amplitude = np.sqrt(0.2)
N = ift.ScalingOperator(s_space, noise_amplitude**2) N = ift.ScalingOperator(s_space, noise_amplitude**2)
n = ift.Field.from_random(domain=s_space, random_type='normal', n = ift.Field.from_random(domain=s_space, random_type='normal',
std=noise_amplitude, mean=0) std=noise_amplitude, mean=0)
d = noiseless_data + n d = noiseless_data + n
j = R.adjoint_times(N.inverse_times(d)) j = R.adjoint_times(N.inverse_times(d))
curv = Curvature(R=R, N=N, Sh=Sh) curv = Curvature(R=R, N=N, Sh=Sh)
D = curv.inverse D = curv.inverse
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Run Wiener Filter ### Run Wiener Filter
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
m = D(j) m = D(j)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Signal Reconstruction ### Signal Reconstruction
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Get signal data and reconstruction data # Get signal data and reconstruction data
s_data = HT(sh).val s_data = HT(sh).val
m_data = HT(m).val m_data = HT(m).val
d_data = d.val d_data = d.val
plt.figure(figsize=(15,10)) plt.figure(figsize=(15,10))
plt.plot(s_data, 'r', label="Signal", linewidth=3) plt.plot(s_data, 'r', label="Signal", linewidth=3)
plt.plot(d_data, 'k.', label="Data") plt.plot(d_data, 'k.', label="Data")
plt.plot(m_data, 'k', label="Reconstruction",linewidth=3) plt.plot(m_data, 'k', label="Reconstruction",linewidth=3)
plt.title("Reconstruction") plt.title("Reconstruction")
plt.legend() plt.legend()
plt.show() plt.show()
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.figure(figsize=(15,10)) plt.figure(figsize=(15,10))
plt.plot(s_data - s_data, 'r', label="Signal", linewidth=3) plt.plot(s_data - s_data, 'r', label="Signal", linewidth=3)
plt.plot(d_data - s_data, 'k.', label="Data") plt.plot(d_data - s_data, 'k.', label="Data")
plt.plot(m_data - s_data, 'k', label="Reconstruction",linewidth=3) plt.plot(m_data - s_data, 'k', label="Reconstruction",linewidth=3)
plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5) plt.axhspan(-noise_amplitude,noise_amplitude, facecolor='0.9', alpha=.5)
plt.title("Residuals") plt.title("Residuals")
plt.legend() plt.legend()
plt.show() plt.show()
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Power Spectrum ### Power Spectrum
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
s_power_data = ift.power_analyze(sh).val s_power_data = ift.power_analyze(sh).val
m_power_data = ift.power_analyze(m).val m_power_data = ift.power_analyze(m).val
plt.figure(figsize=(15,10)) plt.figure(figsize=(15,10))
plt.loglog() plt.loglog()
plt.xlim(1, int(N_pixels/2)) plt.xlim(1, int(N_pixels/2))
ymin = min(m_power_data) ymin = min(m_power_data)
plt.ylim(ymin, 1) plt.ylim(ymin, 1)
xs = np.arange(1,int(N_pixels/2),.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(xs, pow_spec(xs), label="True Power Spectrum", color='k',alpha=0.5)
plt.plot(s_power_data, 'r', label="Signal") plt.plot(s_power_data, 'r', label="Signal")
plt.plot(m_power_data, 'k', label="Reconstruction") plt.plot(m_power_data, 'k', label="Reconstruction")
plt.axhline(noise_amplitude**2 / N_pixels, color="k", linestyle='--', label="Noise level", alpha=.5) 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.axhspan(noise_amplitude**2 / N_pixels, ymin, facecolor='0.9', alpha=.5)
plt.title("Power Spectrum") plt.title("Power Spectrum")
plt.legend() plt.legend()
plt.show() plt.show()
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Wiener Filter on Incomplete Data ## Wiener Filter on Incomplete Data
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Operators # Operators
Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec) Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
N = ift.ScalingOperator(s_space, noise_amplitude**2) N = ift.ScalingOperator(s_space, noise_amplitude**2)
# R is defined below # R is defined below
# Fields # Fields
sh = Sh.draw_sample_with_dtype(dtype=np.float64) sh = Sh.draw_sample_with_dtype(dtype=np.float64)
s = HT(sh) s = HT(sh)
n = ift.Field.from_random(domain=s_space, random_type='normal', n = ift.Field.from_random(domain=s_space, random_type='normal',
std=noise_amplitude, mean=0) std=noise_amplitude, mean=0)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Partially Lose Data ### Partially Lose Data
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
l = int(N_pixels * 0.2) l = int(N_pixels * 0.2)
h = int(N_pixels * 0.2 * 2) h = int(N_pixels * 0.2 * 2)
mask = np.full(s_space.shape, 1.) mask = np.full(s_space.shape, 1.)
mask[l:h] = 0 mask[l:h] = 0
mask = ift.Field.from_raw(s_space, mask) mask = ift.Field.from_raw(s_space, mask)
R = ift.DiagonalOperator(mask)(HT) R = ift.DiagonalOperator(mask)(HT)
n = n.val_rw() n = n.val_rw()
n[l:h] = 0 n[l:h] = 0
n = ift.Field.from_raw(s_space, n) n = ift.Field.from_raw(s_space, n)
d = R(sh) + n d = R(sh) + n
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
curv = Curvature(R=R, N=N, Sh=Sh) curv = Curvature(R=R, N=N, Sh=Sh)
D = curv.inverse D = curv.inverse
j = R.adjoint_times(N.inverse_times(d)) j = R.adjoint_times(N.inverse_times(d))
m = D(j) m = D(j)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Compute Uncertainty ### Compute Uncertainty
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 200, np.float64) m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 200, np.float64)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Get data ### Get data
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
# Get signal data and reconstruction data # Get signal data and reconstruction data
s_data = s.val s_data = s.val
m_data = HT(m).val m_data = HT(m).val
m_var_data = m_var.val m_var_data = m_var.val
uncertainty = np.sqrt(m_var_data) uncertainty = np.sqrt(m_var_data)
d_data = d.val_rw() d_data = d.val_rw()
# Set lost data to NaN for proper plotting # Set lost data to NaN for proper plotting
d_data[d_data == 0] = np.nan d_data[d_data == 0] = np.nan
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
fig = plt.figure(figsize=(15,10)) fig = plt.figure(figsize=(15,10))
plt.axvspan(l, h, facecolor='0.8',alpha=0.5) 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.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(s_data, 'r', label="Signal", alpha=1, linewidth=3)
plt.plot(d_data, 'k.', label="Data") plt.plot(d_data, 'k.', label="Data")
plt.plot(m_data, 'k', label="Reconstruction", linewidth=3) plt.plot(m_data, 'k', label="Reconstruction", linewidth=3)
plt.title("Reconstruction of incomplete data") plt.title("Reconstruction of incomplete data")
plt.legend() plt.legend()
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# 2d Example # 2d Example
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
N_pixels = 256 # Number of pixels N_pixels = 256 # Number of pixels
sigma2 = 2. # Noise variance sigma2 = 2. # Noise variance
def pow_spec(k): def pow_spec(k):
P0, k0, gamma = [.2, 2, 4] P0, k0, gamma = [.2, 2, 4]
return P0 * (1. + (k/k0)**2)**(-gamma/2) return P0 * (1. + (k/k0)**2)**(-gamma/2)
s_space = ift.RGSpace([N_pixels, N_pixels]) s_space = ift.RGSpace([N_pixels, N_pixels])
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
h_space = s_space.get_default_codomain() h_space = s_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(h_space,s_space) HT = ift.HarmonicTransformOperator(h_space,s_space)
# Operators # Operators
Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec) Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
N = ift.ScalingOperator(s_space, sigma2) N = ift.ScalingOperator(s_space, sigma2)
# Fields and data # Fields and data
sh = Sh.draw_sample_with_dtype(dtype=np.float64) sh = Sh.draw_sample_with_dtype(dtype=np.float64)
n = ift.Field.from_random(domain=s_space, random_type='normal', n = ift.Field.from_random(domain=s_space, random_type='normal',
std=np.sqrt(sigma2), mean=0) std=np.sqrt(sigma2), mean=0)
# Lose some data # Lose some data
l = int(N_pixels * 0.33) l = int(N_pixels * 0.33)
h = int(N_pixels * 0.33 * 2) h = int(N_pixels * 0.33 * 2)
mask = np.full(s_space.shape, 1.) mask = np.full(s_space.shape, 1.)
mask[l:h,l:h] = 0. mask[l:h,l:h] = 0.
mask = ift.Field.from_raw(s_space, mask) mask = ift.Field.from_raw(s_space, mask)
R = ift.DiagonalOperator(mask)(HT) R = ift.DiagonalOperator(mask)(HT)
n = n.val_rw() n = n.val_rw()
n[l:h, l:h] = 0 n[l:h, l:h] = 0
n = ift.Field.from_raw(s_space, n) n = ift.Field.from_raw(s_space, n)
curv = Curvature(R=R, N=N, Sh=Sh) curv = Curvature(R=R, N=N, Sh=Sh)
D = curv.inverse D = curv.inverse
d = R(sh) + n d = R(sh) + n
j = R.adjoint_times(N.inverse_times(d)) j = R.adjoint_times(N.inverse_times(d))
# Run Wiener filter # Run Wiener filter
m = D(j) m = D(j)
# Uncertainty # Uncertainty
m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 20, np.float64) m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 20, np.float64)
# Get data # Get data
s_data = HT(sh).val s_data = HT(sh).val
m_data = HT(m).val m_data = HT(m).val
m_var_data = m_var.val m_var_data = m_var.val
d_data = d.val d_data = d.val
uncertainty = np.sqrt(np.abs(m_var_data)) uncertainty = np.sqrt(np.abs(m_var_data))
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
cm = ['magma', 'inferno', 'plasma', 'viridis'][1] cm = ['magma', 'inferno', 'plasma', 'viridis'][1]
mi = np.min(s_data) mi = np.min(s_data)
ma = np.max(s_data) ma = np.max(s_data)
fig, axes = plt.subplots(1, 2, figsize=(15, 7)) fig, axes = plt.subplots(1, 2, figsize=(15, 7))
data = [s_data, d_data] data = [s_data, d_data]
caption = ["Signal", "Data"] caption = ["Signal", "Data"]
for ax in axes.flat: for ax in axes.flat:
im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi, im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi,
vmax=ma) vmax=ma)
ax.set_title(caption.pop(0)) ax.set_title(caption.pop(0))
fig.subplots_adjust(right=0.8) fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7]) cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
fig.colorbar(im, cax=cbar_ax) fig.colorbar(im, cax=cbar_ax)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
mi = np.min(s_data) mi = np.min(s_data)
ma = np.max(s_data) ma = np.max(s_data)
fig, axes = plt.subplots(3, 2, figsize=(15, 22.5)) fig, axes = plt.subplots(3, 2, figsize=(15, 22.5))
sample = HT(curv.draw_sample(from_inverse=True)+m).val sample = HT(curv.draw_sample(from_inverse=True)+m).val
post_mean = (m_mean + HT(m)).val post_mean = (m_mean + HT(m)).val
data = [s_data, m_data, post_mean, sample, s_data - m_data, uncertainty] data = [s_data, m_data, post_mean, sample, s_data - m_data, uncertainty]
caption = ["Signal", "Reconstruction", "Posterior mean", "Sample", "Residuals", "Uncertainty Map"] caption = ["Signal", "Reconstruction", "Posterior mean", "Sample", "Residuals", "Uncertainty Map"]
for ax in axes.flat: for ax in axes.flat:
im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi, vmax=ma) im = ax.imshow(data.pop(0), interpolation='nearest', cmap=cm, vmin=mi, vmax=ma)
ax.set_title(caption.pop(0)) ax.set_title(caption.pop(0))
fig.subplots_adjust(right=0.8) fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([.85, 0.15, 0.05, 0.7]) cbar_ax = fig.add_axes([.85, 0.15, 0.05, 0.7])
fig.colorbar(im, cax=cbar_ax) fig.colorbar(im, cax=cbar_ax)
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Is the uncertainty map reliable? ### Is the uncertainty map reliable?
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
precise = (np.abs(s_data-m_data) < uncertainty) precise = (np.abs(s_data-m_data) < uncertainty)
print("Error within uncertainty map bounds: " + str(np.sum(precise) * 100 / N_pixels**2) + "%") print("Error within uncertainty map bounds: " + str(np.sum(precise) * 100 / N_pixels**2) + "%")
plt.figure(figsize=(15,10)) plt.figure(figsize=(15,10))
plt.imshow(precise.astype(float), cmap="brg") plt.imshow(precise.astype(float), cmap="brg")
plt.colorbar() plt.colorbar()
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Start Coding # Start Coding
## NIFTy Repository + Installation guide ## NIFTy Repository + Installation guide
https://gitlab.mpcdf.mpg.de/ift/NIFTy https://gitlab.mpcdf.mpg.de/ift/NIFTy
NIFTy v6 **more or less stable!** NIFTy v6 **more or less stable!**
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
......
...@@ -25,7 +25,7 @@ import sys ...@@ -25,7 +25,7 @@ import sys
import numpy as np import numpy as np
import nifty6 as ift import nifty7 as ift
def make_checkerboard_mask(position_space): def make_checkerboard_mask(position_space):
......
...@@ -25,7 +25,7 @@ import sys ...@@ -25,7 +25,7 @@ import sys
import numpy as np import numpy as np
import nifty6 as ift import nifty7 as ift
def exposure_2d(): def exposure_2d():
......
...@@ -29,7 +29,7 @@ import sys ...@@ -29,7 +29,7 @@ import sys
import numpy as np import numpy as np
import nifty6 as ift import nifty7 as ift
def random_los(n_los): def random_los(n_los):
......
...@@ -29,7 +29,7 @@ import sys ...@@ -29,7 +29,7 @@ import sys
import numpy as np import numpy as np
import nifty6 as ift import nifty7 as ift
class SingleDomain(ift.LinearOperator): class SingleDomain(ift.LinearOperator):
......
import numpy as np import numpy as np
import nifty6 as ift import nifty7 as ift
def convtest(test_signal, delta, func): def convtest(test_signal, delta, func):
......
...@@ -18,7 +18,7 @@ ...@@ -18,7 +18,7 @@
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
import nifty6 as ift import nifty7 as ift
def polynomial(coefficients, sampling_points): def polynomial(coefficients, sampling_points):
......
rm -rf docs/build docs/source/mod rm -rf docs/build docs/source/mod
sphinx-apidoc -e -o docs/source/mod nifty6 sphinx-apidoc -e -o docs/source/mod nifty7
sphinx-build -b html docs/source/ docs/build/ sphinx-build -b html docs/source/ docs/build/
...@@ -13,7 +13,7 @@ NIFTy-related publications ...@@ -13,7 +13,7 @@ NIFTy-related publications
author = {{Martin Reinecke, Theo Steininger, Marco Selig}}, author = {{Martin Reinecke, Theo Steininger, Marco Selig}},
title = {NIFTy -- Numerical Information Field TheorY}, title = {NIFTy -- Numerical Information Field TheorY},
url = {https://gitlab.mpcdf.mpg.de/ift/NIFTy}, url = {https://gitlab.mpcdf.mpg.de/ift/NIFTy},
version = {nifty6}, version = {nifty7},
date = {2018-04-05}, date = {2018-04-05},
} }
......
...@@ -36,9 +36,9 @@ Domains ...@@ -36,9 +36,9 @@ Domains
Abstract base class Abstract base class
------------------- -------------------
.. currentmodule:: nifty6.domains.domain .. currentmodule:: nifty7.domains.domain
One of the fundamental building blocks of the NIFTy6 framework is the *domain*. One of the fundamental building blocks of the NIFTy7 framework is the *domain*.
Its required capabilities are expressed by the abstract :py:class:`Domain` class. Its required capabilities are expressed by the abstract :py:class:`Domain` class.
A domain must be able to answer the following queries: A domain must be able to answer the following queries:
...@@ -52,7 +52,7 @@ A domain must be able to answer the following queries: ...@@ -52,7 +52,7 @@ A domain must be able to answer the following queries:
Unstructured domains Unstructured domains
-------------------- --------------------
.. currentmodule:: nifty6.domains.unstructured_domain .. currentmodule:: nifty7.domains.unstructured_domain
Domains can be either *structured* (i.e. there is geometrical information Domains can be either *structured* (i.e. there is geometrical information
associated with them, like position in space and volume factors), associated with them, like position in space and volume factors),
...@@ -65,7 +65,7 @@ Unstructured domains can be described by instances of NIFTy's ...@@ -65,7 +65,7 @@ Unstructured domains can be described by instances of NIFTy's
Structured domains Structured domains
------------------ ------------------
.. currentmodule:: nifty6.domains.structured_domain .. currentmodule:: nifty7.domains.structured_domain
In contrast to unstructured domains, these domains have an assigned geometry. In contrast to unstructured domains, these domains have an assigned geometry.
NIFTy requires them to provide the volume elements of their grid cells. NIFTy requires them to provide the volume elements of their grid cells.
...@@ -86,7 +86,7 @@ The additional methods are specified in the abstract class ...@@ -86,7 +86,7 @@ The additional methods are specified in the abstract class
NIFTy comes with several concrete subclasses of :class:`StructuredDomain`: NIFTy comes with several concrete subclasses of :class:`StructuredDomain`:
.. currentmodule:: nifty6.domains .. currentmodule:: nifty7.domains
- :class:`~rg_space.RGSpace` represents a regular Cartesian grid with an arbitrary - :class:`~rg_space.RGSpace` represents a regular Cartesian grid with an arbitrary
number of dimensions, which is supposed to be periodic in each dimension. number of dimensions, which is supposed to be periodic in each dimension.
...@@ -123,7 +123,7 @@ Some examples are: ...@@ -123,7 +123,7 @@ Some examples are:
(on a harmonic domain) and a few inferred model parameters (e.g. on an (on a harmonic domain) and a few inferred model parameters (e.g. on an
unstructured grid). unstructured grid).
.. currentmodule:: nifty6 .. currentmodule:: nifty7
Consequently, NIFTy defines a class called :class:`~domain_tuple.DomainTuple` Consequently, NIFTy defines a class called :class:`~domain_tuple.DomainTuple`
holding a sequence of :class:`~domains.domain.Domain` objects. The full domain is holding a sequence of :class:`~domains.domain.Domain` objects. The full domain is
...@@ -228,7 +228,7 @@ Advanced operators ...@@ -228,7 +228,7 @@ Advanced operators
NIFTy provides a library of commonly employed operators which can be used for NIFTy provides a library of commonly employed operators which can be used for
specific inference problems. Currently these are: specific inference problems. Currently these are:
.. currentmodule:: nifty6.library .. currentmodule:: nifty7.library
- :class:`~smooth_linear_amplitude.SLAmplitude`, which returns a smooth power spectrum. - :class:`~smooth_linear_amplitude.SLAmplitude`, which returns a smooth power spectrum.
- :class:`~inverse_gamma_operator.InverseGammaOperator`, which models point sources which are - :class:`~inverse_gamma_operator.InverseGammaOperator`, which models point sources which are
...@@ -240,9 +240,9 @@ specific inference problems. Currently these are: ...@@ -240,9 +240,9 @@ specific inference problems. Currently these are:
Linear Operators Linear Operators
================ ================
.. currentmodule:: nifty6.operators .. currentmodule:: nifty7.operators
A linear operator (represented by NIFTy6's abstract A linear operator (represented by NIFTy7's abstract
:class:`~linear_operator.LinearOperator` class) is derived from :class:`~linear_operator.LinearOperator` class) is derived from
:class:`~operator.Operator` and can be interpreted as an (implicitly defined) :class:`~operator.Operator` and can be interpreted as an (implicitly defined)
matrix. Since its operation is linear, it can provide some additional matrix. Since its operation is linear, it can provide some additional
...@@ -282,7 +282,7 @@ This functionality is provided by NIFTy's ...@@ -282,7 +282,7 @@ This functionality is provided by NIFTy's
:class:`~inversion_enabler.InversionEnabler` class, which is itself a linear :class:`~inversion_enabler.InversionEnabler` class, which is itself a linear
operator. operator.
.. currentmodule:: nifty6.operators.operator .. currentmodule:: nifty7.operators.operator
Direct multiplication and adjoint inverse multiplication transform a field Direct multiplication and adjoint inverse multiplication transform a field
defined on the operator's :attr:`~Operator.domain` to one defined on the defined on the operator's :attr:`~Operator.domain` to one defined on the
...@@ -290,7 +290,7 @@ operator's :attr:`~Operator.target`, whereas adjoint multiplication and inverse ...@@ -290,7 +290,7 @@ operator's :attr:`~Operator.target`, whereas adjoint multiplication and inverse
multiplication transform from :attr:`~Operator.target` to multiplication transform from :attr:`~Operator.target` to
:attr:`~Operator.domain`. :attr:`~Operator.domain`.
.. currentmodule:: nifty6.operators .. currentmodule:: nifty7.operators
Operators with identical domain and target can be derived from Operators with identical domain and target can be derived from
:class:`~endomorphic_operator.EndomorphicOperator`. Typical examples for this :class:`~endomorphic_operator.EndomorphicOperator`. Typical examples for this
...@@ -299,7 +299,7 @@ multiplies its input by a scalar value, and ...@@ -299,7 +299,7 @@ multiplies its input by a scalar value, and
:class:`~diagonal_operator.DiagonalOperator`, which multiplies every value of :class:`~diagonal_operator.DiagonalOperator`, which multiplies every value of
its input field with potentially different values. its input field with potentially different values.
.. currentmodule:: nifty6 .. currentmodule:: nifty7
Further operator classes provided by NIFTy are Further operator classes provided by NIFTy are
...@@ -325,7 +325,7 @@ and ``f1`` and ``f2`` are of type :class:`~field.Field`, writing:: ...@@ -325,7 +325,7 @@ and ``f1`` and ``f2`` are of type :class:`~field.Field`, writing::
X = A(B.inverse(A.adjoint)) + C X = A(B.inverse(A.adjoint)) + C
f2 = X(f1) f2 = X(f1)
.. currentmodule:: nifty6.operators.linear_operator .. currentmodule:: nifty7.operators.linear_operator
will perform the operation suggested intuitively by the notation, checking will perform the operation suggested intuitively by the notation, checking
domain compatibility while building the composed operator. domain compatibility while building the composed operator.
...@@ -349,12 +349,12 @@ Minimization ...@@ -349,12 +349,12 @@ Minimization
Most problems in IFT are solved by (possibly nested) minimizations of Most problems in IFT are solved by (possibly nested) minimizations of
high-dimensional functions, which are often nonlinear. high-dimensional functions, which are often nonlinear.
.. currentmodule:: nifty6.minimization .. currentmodule:: nifty7.minimization
Energy functionals Energy functionals
------------------ ------------------
In NIFTy6 such functions are represented by objects of type In NIFTy7 such functions are represented by objects of type
:class:`~energy.Energy`. These hold the prescription how to calculate the :class:`~energy.Energy`. These hold the prescription how to calculate the
function's :attr:`~energy.Energy.value`, :attr:`~energy.Energy.gradient` and function's :attr:`~energy.Energy.value`, :attr:`~energy.Energy.gradient` and
(optionally) :attr:`~energy.Energy.metric` at any given (optionally) :attr:`~energy.Energy.metric` at any given
...@@ -362,11 +362,11 @@ function's :attr:`~energy.Energy.value`, :attr:`~energy.Energy.gradient` and ...@@ -362,11 +362,11 @@ function's :attr:`~energy.Energy.value`, :attr:`~energy.Energy.gradient` and
floating-point scalars, gradients have the form of fields defined on the energy's floating-point scalars, gradients have the form of fields defined on the energy's
position domain, and metrics are represented by linear operator objects. position domain, and metrics are represented by linear operator objects.
.. currentmodule:: nifty6 .. currentmodule:: nifty7
Energies are classes that typically have to be provided by the user when Energies are classes that typically have to be provided by the user when
tackling new IFT problems. An example of concrete energy classes delivered with tackling new IFT problems. An example of concrete energy classes delivered with
NIFTy6 is :class:`~minimization.quadratic_energy.QuadraticEnergy` (with NIFTy7 is :class:`~minimization.quadratic_energy.QuadraticEnergy` (with
position-independent metric, mainly used with conjugate gradient minimization). position-independent metric, mainly used with conjugate gradient minimization).
For MGVI, NIFTy provides the :class:`~energy.Energy` subclass For MGVI, NIFTy provides the :class:`~energy.Energy` subclass
...@@ -394,14 +394,14 @@ functional for MGVI and minimize it. ...@@ -394,14 +394,14 @@ functional for MGVI and minimize it.
Iteration control Iteration control
----------------- -----------------
.. currentmodule:: nifty6.minimization.iteration_controllers .. currentmodule:: nifty7.minimization.iteration_controllers
Iterative minimization of an energy requires some means of Iterative minimization of an energy requires some means of
checking the quality of the current solution estimate and stopping once checking the quality of the current solution estimate and stopping once
it is sufficiently accurate. In case of numerical problems, the iteration needs it is sufficiently accurate. In case of numerical problems, the iteration needs
to be terminated as well, returning a suitable error description. to be terminated as well, returning a suitable error description.
In NIFTy6, this functionality is encapsulated in the abstract In NIFTy7, this functionality is encapsulated in the abstract
:class:`IterationController` class, which is provided with the initial energy :class:`IterationController` class, which is provided with the initial energy
object before starting the minimization, and is updated with the improved