Commit 28fe4cf9 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'docu_fixes' into 'NIFTy_5'

High-level docu

See merge request ift/nifty-dev!153
parents addd0938 17041b3f
......@@ -15,28 +15,33 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import nifty5 as ift
#####################################################################
# Bernoulli reconstruction
# Reconstruct an event probability field with values between 0 and 1
# from the observed events
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
#####################################################################
import numpy as np
import nifty5 as ift
if __name__ == '__main__':
# FIXME ABOUT THIS CODE
np.random.seed(41)
# Set up the position space of the signal
#
# # One dimensional regular grid with uniform exposure
# position_space = ift.RGSpace(1024)
# exposure = np.ones(position_space.shape)
# Two-dimensional regular grid with inhomogeneous exposure
mode = 2
if mode == 0:
# One-dimensional regular grid
position_space = ift.RGSpace(1024)
elif mode == 1:
# Two-dimensional regular grid
position_space = ift.RGSpace([512, 512])
else:
# Sphere
position_space = ift.HPSpace(128)
# Sphere with uniform exposure
# position_space = ift.HPSpace(128)
# exposure = ift.Field.full(position_space, 1.)
# Defining harmonic space and transform
# Define harmonic space and transform
harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
......@@ -44,15 +49,13 @@ if __name__ == '__main__':
# Define power spectrum and amplitudes
def sqrtpspec(k):
return 1. / (20. + k**2)
return 1./(20. + k**2)
A = ift.create_power_operator(harmonic_space, sqrtpspec)
# Set up a sky model
# Set up a sky model and instrumental response
sky = ift.positive_tanh(HT(A))
GR = ift.GeometryRemover(position_space)
# Set up instrumental response
R = GR
# Generate mock data
......@@ -65,8 +68,8 @@ if __name__ == '__main__':
# Compute likelihood and Hamiltonian
position = ift.from_random('normal', harmonic_space)
likelihood = ift.BernoulliEnergy(data)(p)
ic_newton = ift.DeltaEnergyController(name='Newton', iteration_limit=100,
tol_rel_deltaE=1e-8)
ic_newton = ift.DeltaEnergyController(
name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
minimizer = ift.NewtonCG(ic_newton)
ic_sampling = ift.GradientNormController(iteration_limit=100)
......@@ -82,5 +85,4 @@ if __name__ == '__main__':
plot.add(reconstruction, title='reconstruction')
plot.add(GR.adjoint_times(data), title='data')
plot.add(sky(mock_position), title='truth')
plot.output(nx=3, xsize=16, ysize=5, title="results",
name="bernoulli.png")
plot.output(nx=3, xsize=16, ysize=9, title="results", name="bernoulli.png")
......@@ -15,26 +15,36 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import nifty5 as ift
###############################################################################
# Compute a Wiener filter solution with NIFTy
# Shows how measurement gaps are filled in
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
###############################################################################
import numpy as np
import nifty5 as ift
def make_chess_mask(position_space):
def make_checkerboard_mask(position_space):
# Checkerboard mask for 2D mode
mask = np.ones(position_space.shape)
for i in range(4):
for j in range(4):
if (i+j) % 2 == 0:
mask[i*128//4:(i+1)*128//4, j*128//4:(j+1)*128//4] = 0
if (i + j) % 2 == 0:
mask[i*128//4:(i + 1)*128//4, j*128//4:(j + 1)*128//4] = 0
return mask
def make_random_mask():
# Random mask for spherical mode
mask = ift.from_random('pm1', position_space)
mask = (mask+1)/2
mask = (mask + 1)/2
return mask.to_global_data()
def mask_to_nan(mask, field):
# Set masked pixels to nan for plotting
masked_data = field.local_data.copy()
masked_data[mask.local_data == 0] = np.nan
return ift.from_local_data(field.domain, masked_data)
......@@ -42,49 +52,68 @@ def mask_to_nan(mask, field):
if __name__ == '__main__':
np.random.seed(42)
# FIXME description of the tutorial
# Choose problem geometry and masking
# Choose space on which the signal field is defined
mode = 1
if mode == 0:
# One dimensional regular grid
# One-dimensional regular grid
position_space = ift.RGSpace([1024])
mask = np.ones(position_space.shape)
elif mode == 1:
# Two dimensional regular grid with chess mask
# Two-dimensional regular grid with checkerboard mask
position_space = ift.RGSpace([128, 128])
mask = make_chess_mask(position_space)
mask = make_checkerboard_mask(position_space)
else:
# Sphere with half of its locations randomly masked
# Sphere with half of its pixels randomly masked
position_space = ift.HPSpace(128)
mask = make_random_mask()
# Specify harmonic space corresponding to signal
harmonic_space = position_space.get_default_codomain()
# Harmonic transform from harmonic space to position space
HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
# Set correlation structure with a power spectrum and build
# prior correlation covariance
# Set prior correlation covariance with a power spectrum leading to
# homogeneous and isotropic statistics
def power_spectrum(k):
return 100. / (20.+k**3)
return 100./(20. + k**3)
# 1D spectral space on which the power spectrum is defined
power_space = ift.PowerSpace(harmonic_space)
# Mapping to (higher dimensional) harmonic space
PD = ift.PowerDistributor(harmonic_space, power_space)
# Apply the mapping
prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
# Insert the result into the diagonal of an harmonic space operator
S = ift.DiagonalOperator(prior_correlation_structure)
# S is the prior field covariance
# Build instrument response consisting of a discretization, mask
# and harmonic transformaion
# Data is defined on a geometry-free space, thus the geometry is removed
GR = ift.GeometryRemover(position_space)
# Masking operator to model that parts of the field have not been observed
mask = ift.Field.from_global_data(position_space, mask)
Mask = ift.DiagonalOperator(mask)
# Operators can be composed either with paranthesis
# The response operator consists of
# - an harmonic transform (to get to image space)
# - the application of the mask
# - the removal of geometric information
# Operators can be composed either with parenthesis
R = GR(Mask(HT))
# or with @
R = GR @ Mask @ HT
data_space = GR.target
# Set the noise covariance
# Set the noise covariance N
noise = 5.
N = ift.ScalingOperator(noise, data_space)
......@@ -93,17 +122,17 @@ if __name__ == '__main__':
MOCK_NOISE = N.draw_sample()
data = R(MOCK_SIGNAL) + MOCK_NOISE
# Build propagator D and information source j
j = R.adjoint_times(N.inverse_times(data))
# Build inverse propagator D and information source j
D_inv = R.adjoint @ N.inverse @ R + S.inverse
# Make it invertible
j = R.adjoint_times(N.inverse_times(data))
# Make D_inv invertible (via Conjugate Gradient)
IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse
# WIENER FILTER
# Calculate WIENER FILTER solution
m = D(j)
# PLOTTING
# Plotting
rg = isinstance(position_space, ift.RGSpace)
plot = ift.Plot()
if rg and len(position_space.shape) == 1:
......@@ -118,5 +147,5 @@ if __name__ == '__main__':
plot.add(HT(MOCK_SIGNAL), title='Mock Signal')
plot.add(mask_to_nan(mask, (GR(Mask)).adjoint(data)), title='Data')
plot.add(HT(m), title='Reconstruction')
plot.add(mask_to_nan(mask, HT(m-MOCK_SIGNAL)), title='Residuals')
plot.add(mask_to_nan(mask, HT(m - MOCK_SIGNAL)), title='Residuals')
plot.output(nx=2, ny=2, xsize=10, ysize=10, title="getting_started_1")
......@@ -15,11 +15,19 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import nifty5 as ift
###############################################################################
# Log-normal field reconstruction from Poissonian data with inhomogenous
# exposure (in case for 2D mode)
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
###############################################################################
import numpy as np
import nifty5 as ift
def get_2D_exposure():
def exposure_2d():
# Structured exposure for 2D mode
x_shape, y_shape = position_space.shape
exposure = np.ones(position_space.shape)
......@@ -30,72 +38,73 @@ def get_2D_exposure():
exposure[:, x_shape*4//5:x_shape] *= .1
exposure[:, x_shape//2:x_shape*3//2] *= 3.
exposure = ift.Field.from_global_data(position_space, exposure)
return exposure
return ift.Field.from_global_data(position_space, exposure)
if __name__ == '__main__':
# FIXME description of the tutorial
# FIXME All random seeds to 42
np.random.seed(41)
# Set up the position space of the signal
#
# # One dimensional regular grid with uniform exposure
# position_space = ift.RGSpace(1024)
# exposure = ift.Field.full(position_space, 1.)
# Choose space on which the signal field is defined
mode = 2
if mode == 0:
# One-dimensional regular grid with uniform exposure
position_space = ift.RGSpace(1024)
exposure = ift.Field.full(position_space, 1.)
elif mode == 1:
# Two-dimensional regular grid with inhomogeneous exposure
position_space = ift.RGSpace([512, 512])
exposure = get_2D_exposure()
exposure = exposure_2d()
else:
# Sphere with uniform exposure
# position_space = ift.HPSpace(128)
# exposure = ift.Field.full(position_space, 1.)
position_space = ift.HPSpace(128)
exposure = ift.Field.full(position_space, 1.)
# Defining harmonic space and transform
# Define harmonic space and harmonic transform
harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
# Domain on which the field's degrees of freedom are defined
domain = ift.DomainTuple.make(harmonic_space)
position = ift.from_random('normal', domain)
# Define power spectrum and amplitudes
# Define amplitude (square root of power spectrum)
def sqrtpspec(k):
return 1. / (20. + k**2)
return 1./(20. + k**2)
p_space = ift.PowerSpace(harmonic_space)
pd = ift.PowerDistributor(harmonic_space, p_space)
a = ift.PS_field(p_space, sqrtpspec)
A = pd(a)
# Set up a sky model
# Define sky model
sky = ift.exp(HT(ift.makeOp(A)))
M = ift.DiagonalOperator(exposure)
GR = ift.GeometryRemover(position_space)
# Set up instrumental response
# Define instrumental response
R = GR(M)
# Generate mock data
# Generate mock data and define likelihood operator
d_space = R.target[0]
lamb = R(sky)
mock_position = ift.from_random('normal', domain)
data = lamb(mock_position)
data = np.random.poisson(data.to_global_data().astype(np.float64))
data = ift.Field.from_global_data(d_space, data)
# Compute likelihood and Hamiltonian
ic_newton = ift.DeltaEnergyController(name='Newton', iteration_limit=100,
tol_rel_deltaE=1e-8)
likelihood = ift.PoissonianEnergy(data)(lamb)
# Settings for minimization
ic_newton = ift.DeltaEnergyController(
name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
minimizer = ift.NewtonCG(ic_newton)
# Minimize the Hamiltonian
# Compute MAP solution by minimizing the information Hamiltonian
H = ift.Hamiltonian(likelihood)
H = ift.EnergyAdapter(position, H, want_metric=True)
initial_position = ift.from_random('normal', domain)
H = ift.EnergyAdapter(initial_position, H, want_metric=True)
H, convergence = minimizer(H)
# Plot results
# Plotting
signal = sky(mock_position)
reconst = sky(H.position)
plot = ift.Plot()
......@@ -103,4 +112,4 @@ if __name__ == '__main__':
plot.add(GR.adjoint(data), title='Data')
plot.add(reconst, title='Reconstruction')
plot.add(reconst - signal, title='Residuals')
plot.output(name='getting_started_2.png', xsize=16, ysize=16)
plot.output(name='getting_started_2.pdf', xsize=16, ysize=16)
......@@ -15,99 +15,133 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import nifty5 as ift
############################################################
# Non-linear tomography
# The data is integrated lines of sight
# Random lines (set mode=0), radial lines (mode=1)
#############################################################
import numpy as np
import nifty5 as ift
def random_los(n_los):
starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
ends = list(0.5 + 0*np.random.uniform(0, 1, (n_los, 2)).T)
return starts, ends
def get_random_LOS(n_los):
def radial_los(n_los):
starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
return starts, ends
if __name__ == '__main__':
# FIXME description of the tutorial
np.random.seed(42)
np.seterr(all='raise')
position_space = ift.RGSpace([128, 128])
np.random.seed(420)
# Setting up an amplitude model
A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -4., 1, 1., 1.)
# Choose between random line-of-sight response (mode=1) and radial lines
# of sight (mode=2)
mode = 1
# Building the model for a correlated signal
position_space = ift.RGSpace([128, 128])
# Set up an amplitude model for the field
# The parameters mean:
# 64 spectral bins
#
# Spectral smoothness (affects Gaussian process part)
# 3 = relatively high variance of spectral curbvature
# 0.4 = quefrency mode below which cepstrum flattens
#
# Power-law part of spectrum:
# -5 = preferred power-law slope
# 0.5 = low variance of power-law slope
# 0.4 = y-intercept mean
# 0.3 = relatively high y-intercept variance
A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3)
# Build the model for a correlated signal
harmonic_space = position_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
power_space = A.target[0]
power_distributor = ift.PowerDistributor(harmonic_space, power_space)
vol = harmonic_space.scalar_dvol
vol = ift.ScalingOperator(vol**(-0.5), harmonic_space)
vol = ift.ScalingOperator(harmonic_space.scalar_dvol**(-0.5),
harmonic_space)
correlated_field = ht(
vol(power_distributor(A))*ift.ducktape(harmonic_space, None, 'xi'))
# alternatively to the block above one can do:
#correlated_field = ift.CorrelatedField(position_space, A)
# Alternatively, one can use:
# correlated_field = ift.CorrelatedField(position_space, A)
# apply some nonlinearity
# Apply a nonlinearity
signal = ift.positive_tanh(correlated_field)
# Building the Line of Sight response
LOS_starts, LOS_ends = get_random_LOS(100)
# Build the line-of-sight response and define signal response
LOS_starts, LOS_ends = random_los(100) if mode == 1 else radial_los(100)
R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
# build signal response model and model likelihood
signal_response = R(signal)
# specify noise
# Specify noise
data_space = R.target
noise = .001
N = ift.ScalingOperator(noise, data_space)
# generate mock data
MOCK_POSITION = ift.from_random('normal', signal_response.domain)
data = signal_response(MOCK_POSITION) + N.draw_sample()
# Generate mock signal and data
mock_position = ift.from_random('normal', signal_response.domain)
data = signal_response(mock_position) + N.draw_sample()
# set up model likelihood
likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
# set up minimization and inversion schemes
# Minimization parameters
ic_sampling = ift.GradientNormController(iteration_limit=100)
ic_newton = ift.GradInfNormController(
name='Newton', tol=1e-7, iteration_limit=35)
minimizer = ift.NewtonCG(ic_newton)
# build model Hamiltonian
# Set up model likelihood and information Hamiltonian
likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
H = ift.Hamiltonian(likelihood, ic_sampling)
INITIAL_POSITION = ift.MultiField.full(H.domain, 0.)
position = INITIAL_POSITION
initial_position = ift.MultiField.full(H.domain, 0.)
position = initial_position
plot = ift.Plot()
plot.add(signal(MOCK_POSITION), title='Ground Truth')
plot.add(signal(mock_position), title='Ground Truth')
plot.add(R.adjoint_times(data), title='Data')
plot.add([A.force(MOCK_POSITION)], title='Power Spectrum')
plot.add([A.force(mock_position)], title='Power Spectrum')
plot.output(ny=1, nx=3, xsize=24, ysize=6, name="setup.png")
# number of samples used to estimate the KL
N_samples = 20
# Draw new samples to approximate the KL five times
for i in range(5):
# Draw new samples and minimize KL
KL = ift.KL_Energy(position, H, N_samples)
KL, convergence = minimizer(KL)
position = KL.position
# Plot current reconstruction
plot = ift.Plot()
plot.add(signal(KL.position), title="reconstruction")
plot.add([A.force(KL.position), A.force(MOCK_POSITION)], title="power")
plot.add([A.force(KL.position), A.force(mock_position)], title="power")
plot.output(ny=1, ysize=6, xsize=16, name="loop-{:02}.png".format(i))
# Draw posterior samples
KL = ift.KL_Energy(position, H, N_samples)
plot = ift.Plot()
sc = ift.StatCalculator()
for sample in KL.samples:
sc.add(signal(sample + KL.position))
# Plotting
plot = ift.Plot()
plot.add(sc.mean, title="Posterior Mean")
plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
powers = [A.force(s + KL.position) for s in KL.samples]
plot.add(
powers + [A.force(KL.position), A.force(MOCK_POSITION)],
powers + [A.force(KL.position),
A.force(mock_position)],
title="Sampled Posterior Power Spectrum",
linewidth=[1.]*len(powers) + [3., 3.])
plot.output(ny=1, nx=3, xsize=24, ysize=6, name="results.png")
......@@ -97,8 +97,8 @@ Combinations of domains
=======================
The fundamental classes described above are often sufficient to specify the
domain of a field. In some cases, however, it will be necessary to have the
field live on a product of elementary domains instead of a single one.
domain of a field. In some cases, however, it will be necessary to define the
field on a product of elementary domains instead of a single one.
More sophisticated models also require a set of several such fields.
Some examples are:
......@@ -121,7 +121,7 @@ A :class:`DomainTuple` supports iteration and indexing, and also provides the
properties :attr:`~DomainTuple.shape`, :attr:`~DomainTuple.size` in analogy to
the elementary :class:`Domain`.
An aggregation of several :class:`DomainTuple`s, each member identified by a
An aggregation of several :class:`DomainTuple` s, each member identified by a
name, is described by the :class:`MultiDomain` class.
Fields
......@@ -157,11 +157,11 @@ that are not covered by the provided standard operations, its data content must
be extracted first, then changed, and a new field has to be created from the
result.
Fields living on a MultiDomain
Fields defined on a MultiDomain
------------------------------
The :class:`MultiField` class can be seen as a dictionary of individual
:class:`Field`s, each identified by a name, which lives on an associated
:class:`Field` s, each identified by a name, which is defined on a
:class:`MultiDomain`.
......@@ -170,21 +170,21 @@ Operators
All transformations between different NIFTy fields are expressed (explicitly
or implicitly) in the form of :class:`Operator` objects. The interface of this
class is very minimalistic: it has a property called `domain` which returns
a `Domaintuple` or `MultiDomain` object specifying the structure of the
`Field`s or `MultiField`s it expects as input, another property `target`
class is very minimalistic: it has a property called :class:`domain` which returns
a :class:`DomainTuple` or :class:`MultiDomain` object specifying the structure of the
:class:`Field` s or :class:`MultiField` s it expects as input, another property :class:`target`
describing its output, and finally an overloaded `apply` method, which can
take
- a `Field`/`MultiField`object, in which case it returns the transformed
`Field`/`MultiField`
- a `Linearization` object, in which case it returns the transformed
`Linearization`
- a :class:`Field`/:class:`MultiField`object, in which case it returns the transformed
:class:`Field`/:class:`MultiField`
- a :class:`Linearization` object, in which case it returns the transformed
:class:`Linearization`
This is the interface that all objects derived from `Operator` must implement.
In addition, `Operator` objects can be added/subtracted, multiplied, chained
(via the `__call__` method) and support pointwise application of functions like
`exp()`, `log()`, `sqrt()`, `conjugate()` etc.
This is the interface that all objects derived from :class:`Operator` must implement.
In addition, :class:`Operator` objects can be added/subtracted, multiplied, chained
(via the :class:`__call__` method) and support pointwise application of functions like
:class:`exp()`, :class:`log()`, :class:`sqrt()`, :class:`conjugate()` etc.
Linear Operators
......@@ -193,7 +193,7 @@ Linear Operators
A linear operator (represented by NIFTy5's abstract :class:`LinearOperator`
class) is derived from `Operator` and can be interpreted as an
(implicitly defined) matrix. Since its operation is linear, it can provide some
additional functionality which is not available for the more generic `Operator`
additional functionality which is not available for the more generic :class:`Operator`
class.
......
......@@ -9,19 +9,19 @@ Theoretical Background
IFT is fully Bayesian. How else could infinitely many field degrees of freedom be constrained by finite data?
It can be used without the knowledge of Feynman diagrams. There is a full toolbox of methods. It reproduces many known well working algorithms. This should be reassuring. And, there were certainly previous works in a similar spirit. Anyhow, in many cases IFT provides novel rigorous ways to extract information from data.
There is a full toolbox of methods that can be used, like the classical approximation (= Maximum a posteriori = MAP), effective action (= Variational Bayes = VI), Feynman diagrams, renormalitation, and more. IFT reproduces many known well working algorithms. This should be reassuring. And, there were certainly previous works in a similar spirit. Anyhow, in many cases IFT provides novel rigorous ways to extract information from data. NIFTy comes with reimplemented MAP and VI estimators. It also provides a Hamiltonian Monte Carlo sampler for Fields (HMCF).