Commit ca11700c authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'NIFTy_5' into more_tests

parents 96105fd5 bab4d8c1
FROM debian:testing-slim
RUN apt-get update && apt-get install -y \
# Needed for gitlab tests
git \
# Needed for setup
git python3-pip \
# Packages needed for NIFTy
libfftw3-dev \
python3 python3-pip python3-dev python3-scipy \
python3-scipy \
# Documentation build dependencies
python3-sphinx python3-sphinx-rtd-theme \
python3-sphinx-rtd-theme \
# Testing dependencies
python3-coverage python3-pytest python3-pytest-cov \
python3-pytest-cov jupyter \
# Optional NIFTy dependencies
openmpi-bin libopenmpi-dev python3-mpi4py \
# Packages needed for NIFTy
libfftw3-dev python3-mpi4py python3-matplotlib \
# more optional NIFTy dependencies
&& pip3 install pyfftw \
# Optional NIFTy dependencies
&& pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
# Testing dependencies
&& rm -rf /var/lib/apt/lists/*
# Needed for demos to be running
RUN apt-get update && apt-get install -y python3-matplotlib \
&& python3 -m pip install --upgrade pip && python3 -m pip install jupyter \
&& pip3 install jupyter \
&& rm -rf /var/lib/apt/lists/*
# Set matplotlib backend
......
......@@ -29,4 +29,4 @@ add_module_names = False
html_theme = "sphinx_rtd_theme"
html_logo = 'nifty_logo_black.png'
exclude_patterns = ['mod/modules.rst']
exclude_patterns = ['mod/modules.rst', 'mod/*version.rst']
......@@ -219,14 +219,14 @@ Here, the uncertainty of the field and the power spectrum of its generating proc
| **Output of tomography demo getting_started_3.py** |
+----------------------------------------------------+
| .. image:: images/getting_started_3_setup.png |
| :width: 50 % |
| |
+----------------------------------------------------+
| Non-Gaussian signal field, |
| data backprojected into the image domain, power |
| spectrum of underlying Gausssian process. |
+----------------------------------------------------+
| .. image:: images/getting_started_3_results.png |
| :width: 50 % |
| |
+----------------------------------------------------+
| Posterior mean field signal |
| reconstruction, its uncertainty, and the power |
......@@ -235,10 +235,10 @@ Here, the uncertainty of the field and the power spectrum of its generating proc
| orange line). |
+----------------------------------------------------+
Maximim a Posteriori
Maximum a Posteriori
--------------------
One popular field estimation method is Maximim a Posteriori (MAP).
One popular field estimation method is Maximum a Posteriori (MAP).
It only requires to minimize the information Hamiltonian, e.g by a gradient descent method that stops when
......@@ -246,7 +246,7 @@ It only requires to minimize the information Hamiltonian, e.g by a gradient desc
\frac{\partial \mathcal{H}(d,\xi)}{\partial \xi} = 0.
NIFTy5 automatically calculates the necessary gradient from a generative model of the signal and the data and to minimize the Hamiltonian.
NIFTy5 automatically calculates the necessary gradient from a generative model of the signal and the data and uses this to minimize the Hamiltonian.
However, MAP often provides unsatisfactory results in cases of deep hirachical Bayesian networks.
The reason for this is that MAP ignores the volume factors in parameter space, which are not to be neglected in deciding whether a solution is reasonable or not.
......@@ -270,7 +270,7 @@ As a compromise between being optimal and being computationally affordable, the
\int \mathcal{D}\xi \,\mathcal{Q}(\xi) \log \left( \frac{\mathcal{Q}(\xi)}{\mathcal{P}(\xi)} \right)
Minimizing this with respect to all entries of the covariance :math:`D` is unfeasible for fields.
Therefore, Metric Gaussian Variational Inference (MGVI) approximates the precision matrix at the location of the current mean :math:`M=D^{-1}` by the Bayesian Fisher information metric,
Therefore, Metric Gaussian Variational Inference (MGVI) approximates the posterior precision matrix :math:`D^{-1}` at the location of the current mean :math:`m` by the Bayesian Fisher information metric,
.. math::
......
......@@ -64,7 +64,7 @@ from .minimization.nonlinear_cg import NonlinearCG
from .minimization.descent_minimizers import (
DescentMinimizer, SteepestDescent, VL_BFGS, L_BFGS, RelaxedNewton,
NewtonCG)
from .minimization.scipy_minimizer import (ScipyMinimizer, L_BFGS_B, ScipyCG)
from .minimization.scipy_minimizer import L_BFGS_B
from .minimization.energy import Energy
from .minimization.quadratic_energy import QuadraticEnergy
from .minimization.energy_adapter import EnergyAdapter
......@@ -73,7 +73,7 @@ from .minimization.metric_gaussian_kl import MetricGaussianKL
from .sugar import *
from .plot import Plot
from .library.smooth_linear_amplitude import (SLAmplitude, CepstrumOperator)
from .library.smooth_linear_amplitude import SLAmplitude, CepstrumOperator
from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse
from .library.dynamic_operator import (dynamic_operator,
......
......@@ -19,6 +19,7 @@ import numpy as np
from .field import Field
from .linearization import Linearization
from .operators.linear_operator import LinearOperator
from .sugar import from_random
__all__ = ["consistency_check", "check_value_gradient_consistency",
......@@ -62,8 +63,40 @@ def _full_implementation(op, domain_dtype, target_dtype, atol, rtol):
_inverse_implementation(op, domain_dtype, target_dtype, atol, rtol)
def _check_linearity(op, domain_dtype, atol, rtol):
fld1 = from_random("normal", op.domain, dtype=domain_dtype)
fld2 = from_random("normal", op.domain, dtype=domain_dtype)
alpha = np.random.random()
val1 = op(alpha*fld1+fld2)
val2 = alpha*op(fld1)+op(fld2)
_assert_allclose(val1, val2, atol=atol, rtol=rtol)
def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
atol=0, rtol=1e-7):
"""Checks whether times(), adjoint_times(), inverse_times() and
adjoint_inverse_times() (if in capability list) is implemented
consistently. Additionally, it checks whether the operator is linear
actually.
Parameters
----------
op : LinearOperator
Operator which shall be checked.
domain_dtype : FIXME
The data type of the random vectors in the operator's domain. Default
is `np.float64`.
target_dtype : FIXME
The data type of the random vectors in the operator's target. Default
is `np.float64`.
atol : float
FIXME. Default is 0.
rtol : float
FIXME. Default is 0.
"""
if not isinstance(op, LinearOperator):
raise TypeError('This test tests only linear operators.')
_check_linearity(op, domain_dtype, atol, rtol)
_full_implementation(op, domain_dtype, target_dtype, atol, rtol)
_full_implementation(op.adjoint, target_dtype, domain_dtype, atol, rtol)
_full_implementation(op.inverse, target_dtype, domain_dtype, atol, rtol)
......@@ -124,8 +157,10 @@ def _check_consistency(op, loc, tol, ntries, do_metric):
def check_value_gradient_consistency(op, loc, tol=1e-8, ntries=100):
"""FIXME"""
_check_consistency(op, loc, tol, ntries, False)
def check_value_gradient_metric_consistency(op, loc, tol=1e-8, ntries=100):
"""FIXME"""
_check_consistency(op, loc, tol, ntries, True)
......@@ -31,13 +31,13 @@ class Field(object):
Parameters
----------
domain : DomainTuple
the domain of the new Field
The domain of the new Field.
val : data_object
This object's global shape must match the domain shape
After construction, the object will no longer be writeable!
Notes
-----
Note
----
If possible, do not invoke the constructor directly, but use one of the
many convenience functions for instantiation!
"""
......@@ -76,14 +76,14 @@ class Field(object):
Parameters
----------
domain : Domain, tuple of Domain, or DomainTuple
domain of the new Field
Domain of the new Field.
val : float/complex/int scalar
fill value. Data type of the field is inferred from val.
Fill value. Data type of the field is inferred from val.
Returns
-------
Field
the newly created field
The newly created Field.
"""
if not np.isscalar(val):
raise TypeError("val must be a scalar")
......@@ -99,7 +99,7 @@ class Field(object):
Parameters
----------
domain : DomainTuple, tuple of Domain, or Domain
the domain of the new Field
The domain of the new Field.
arr : numpy.ndarray
The data content to be used for the new Field.
Its shape must match the shape of `domain`.
......@@ -132,8 +132,9 @@ class Field(object):
Returns
-------
numpy.ndarray : array containing all field entries, which can be
modified. Its shape is identical to `self.shape`.
numpy.ndarray
Array containing all field entries, which can be modified. Its
shape is identical to `self.shape`.
"""
return dobj.to_global_data_rw(self._val)
......@@ -171,9 +172,9 @@ class Field(object):
random_type : 'pm1', 'normal', or 'uniform'
The random distribution to use.
domain : DomainTuple
The domain of the output random field
The domain of the output random Field.
dtype : type
The datatype of the output random field
The datatype of the output random Field.
Returns
-------
......@@ -187,10 +188,10 @@ class Field(object):
@property
def val(self):
"""dobj.data_object : the data object storing the field's entries
"""dobj.data_object : the data object storing the field's entries.
Notes
-----
Note
----
This property is intended for low-level, internal use only. Do not use
from outside of NIFTy's core; there should be better alternatives.
"""
......@@ -236,13 +237,13 @@ class Field(object):
Parameters
----------
spaces : int, tuple of int or None
indices of the sub-domains of the field's domain to be considered.
Indices of the sub-domains of the field's domain to be considered.
If `None`, the entire domain is used.
Returns
-------
float or None
if the requested sub-domain has a uniform volume element, it is
If the requested sub-domain has a uniform volume element, it is
returned. Otherwise, `None` is returned.
"""
if np.isscalar(spaces):
......@@ -264,7 +265,7 @@ class Field(object):
Parameters
----------
spaces : int, tuple of int or None
indices of the sub-domains of the field's domain to be considered.
Indices of the sub-domains of the field's domain to be considered.
If `None`, the entire domain is used.
Returns
......@@ -331,8 +332,9 @@ class Field(object):
x : Field
Returns
----------
Field, defined on the product space of self.domain and x.domain
-------
Field
Defined on the product space of self.domain and x.domain.
"""
if not isinstance(x, Field):
raise TypeError("The multiplier must be an instance of " +
......
......@@ -18,7 +18,8 @@
from ..minimization.energy_adapter import EnergyAdapter
from ..multi_field import MultiField
from ..operators.distributors import PowerDistributor
from ..operators.energy_operators import StandardHamiltonian, InverseGammaLikelihood
from ..operators.energy_operators import (StandardHamiltonian,
InverseGammaLikelihood)
from ..operators.scaling_operator import ScalingOperator
from ..operators.simple_linear_operators import ducktape
......@@ -72,7 +73,8 @@ def make_adjust_variances(a,
if scaling is not None:
x = ScalingOperator(scaling, x.target)(x)
return StandardHamiltonian(InverseGammaLikelihood(d_eval)(x), ic_samp=ic_samp)
return StandardHamiltonian(InverseGammaLikelihood(d_eval/2.)(x),
ic_samp=ic_samp)
def do_adjust_variances(position,
......
......@@ -16,6 +16,7 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
from operator import mul
from ..domain_tuple import DomainTuple
from ..operators.contraction_operator import ContractionOperator
......@@ -61,6 +62,10 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
power_distributor = PowerDistributor(h_space, p_space)
A = power_distributor(amplitude_operator)
vol = h_space.scalar_dvol**-0.5
# When doubling the resolution of `tgt` the value of the highest k-mode
# will scale with a square root. `vol` cancels this effect such that the
# same power spectrum can be used for the spaces with the same volume,
# different resolutions and the same object in them.
return ht(vol*A*ducktape(h_space, None, name))
......@@ -107,7 +112,8 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
d = [dd0, dd1]
a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
a = reduce(lambda x, y: x*y, a)
a = reduce(mul, a)
A = pd @ a
vol = reduce(lambda x, y: x*y, [sp.scalar_dvol**-0.5 for sp in hsp])
# For `vol` see comment in `CorrelatedField`
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
return ht(vol*A*ducktape(hsp, None, name))
......@@ -26,12 +26,13 @@ from ..sugar import makeOp
class InverseGammaOperator(Operator):
"""Operator which transforms a Gaussian into an inverse gamma distribution.
"""Transforms a Gaussian into an inverse gamma distribution.
The pdf of the inverse gamma distribution is defined as follows:
.. math ::
\\frac{q^\\alpha}{\\Gamma(\\alpha)}x^{-\\alpha -1}\\exp \\left(-\\frac{q}{x}\\right)
.. math::
\\frac{q^\\alpha}{\\Gamma(\\alpha)}x^{-\\alpha -1}
\\exp \\left(-\\frac{q}{x}\\right)
That means that for large x the pdf falls off like :math:`x^(-\\alpha -1)`.
The mean of the pdf is at :math:`q / (\\alpha - 1)` if :math:`\\alpha > 1`.
......@@ -54,7 +55,8 @@ class InverseGammaOperator(Operator):
"""
def __init__(self, domain, alpha, q, delta=0.001):
self._domain = self._target = DomainTuple.make(domain)
self._alpha, self._q, self._delta = float(alpha), float(q), float(delta)
self._alpha, self._q, self._delta = \
float(alpha), float(q), float(delta)
self._xmin, self._xmax = -8.2, 8.2
# Precompute
xs = np.arange(self._xmin, self._xmax+2*delta, delta)
......
......@@ -138,7 +138,8 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
of the linear function.
The prior on the smooth component is parametrized by two real numbers: the
strength and the cutoff of the smoothness prior (see :class:`CepstrumOperator`).
strength and the cutoff of the smoothness prior
(see :class:`CepstrumOperator`).
Parameters
----------
......
......@@ -130,4 +130,4 @@ class MetricGaussianKL(Energy):
def __repr__(self):
return 'KL ({} samples):\n'.format(len(
self._samples)) + utilities.indent(self._ham.__repr__())
self._samples)) + utilities.indent(self._hamiltonian.__repr__())
......@@ -93,7 +93,7 @@ class _MinHelper(object):
return _toArray_rw(res)
class ScipyMinimizer(Minimizer):
class _ScipyMinimizer(Minimizer):
"""Scipy-based minimizer
Parameters
......@@ -136,19 +136,19 @@ class ScipyMinimizer(Minimizer):
def L_BFGS_B(ftol, gtol, maxiter, maxcor=10, disp=False, bounds=None):
"""Returns a ScipyMinimizer object carrying out the L-BFGS-B algorithm.
"""Returns a _ScipyMinimizer object carrying out the L-BFGS-B algorithm.
See Also
--------
ScipyMinimizer
_ScipyMinimizer
"""
options = {"ftol": ftol, "gtol": gtol, "maxiter": maxiter,
"maxcor": maxcor, "disp": disp}
return ScipyMinimizer("L-BFGS-B", options, False, bounds)
return _ScipyMinimizer("L-BFGS-B", options, False, bounds)
class ScipyCG(Minimizer):
"""Returns a ScipyMinimizer object carrying out the conjugate gradient
class _ScipyCG(Minimizer):
"""Returns a _ScipyMinimizer object carrying out the conjugate gradient
algorithm as implemented by SciPy.
This class is only intended for double-checking NIFTy's own conjugate
......
......@@ -24,7 +24,8 @@ from .linear_operator import LinearOperator
class ContractionOperator(LinearOperator):
"""A :class:`LinearOperator` which sums up fields into the direction of subspaces.
"""A :class:`LinearOperator` which sums up fields into the direction of
subspaces.
This Operator sums up a field with is defined on a :class:`DomainTuple`
to a :class:`DomainTuple` which contains the former as a subset.
......
......@@ -23,31 +23,38 @@ from .linear_operator import LinearOperator
class DomainTupleFieldInserter(LinearOperator):
"""Writes the content of a :class:`Field` into one slice of a :class:`DomainTuple`.
"""Writes the content of a :class:`Field` into one slice of a
:class:`DomainTuple`.
Parameters
----------
domain : Domain, tuple of Domain or DomainTuple
new_space : Domain, tuple of Domain or DomainTuple
index : Integer
Index at which new_space shall be added to domain.
position : tuple
Slice in new_space in which the input field shall be written into.
target : Domain, tuple of Domain or DomainTuple
space : int
The index of the sub-domain which is inserted.
index : tuple
Slice in new sub-domain in which the input field shall be written into.
"""
def __init__(self, domain, new_space, index, position):
self._domain = DomainTuple.make(domain)
tgt = list(self.domain)
tgt.insert(index, new_space)
self._target = DomainTuple.make(tgt)
def __init__(self, target, space, pos):
if not space <= len(target) or space < 0:
raise ValueError
self._target = DomainTuple.make(target)
dom = list(self.target)
dom.pop(space)
self._domain = DomainTuple.make(dom)
self._capability = self.TIMES | self.ADJOINT_TIMES
fst_dims = sum(len(dd.shape) for dd in self.domain[:index])
new_space = target[space]
nshp = new_space.shape
if len(position) != len(nshp):
fst_dims = sum(len(dd.shape) for dd in self.target[:space])
if len(pos) != len(nshp):
raise ValueError("shape mismatch between new_space and position")
for s, p in zip(nshp, position):
for s, p in zip(nshp, pos):
if p < 0 or p >= s:
raise ValueError("bad position value")
self._slc = (slice(None),)*fst_dims + position
self._slc = (slice(None),)*fst_dims + pos
def apply(self, x, mode):
self._check_input(x, mode)
......
......@@ -21,8 +21,8 @@ from .linear_operator import LinearOperator
class EndomorphicOperator(LinearOperator):
"""Represents a :class:`LinearOperator` which is endomorphic, i.e. one which
has identical domain and target.
"""Represents a :class:`LinearOperator` which is endomorphic, i.e. one
which has identical domain and target.
"""
@property
def target(self):
......
......@@ -95,7 +95,7 @@ class QuadraticFormOperator(EnergyOperator):
class GaussianEnergy(EnergyOperator):
"""Class for energies of fields with Gaussian probability distribution.
"""Computes a negative-log Gaussian.
Represents up to constants in :math:`m`:
......@@ -162,8 +162,8 @@ class GaussianEnergy(EnergyOperator):
class PoissonianEnergy(EnergyOperator):
"""Class for likelihood Hamiltonians of expected count field constrained
by Poissonian count data.
"""Computes likelihood Hamiltonians of expected count field constrained by
Poissonian count data.
Represents up to an f-independent term :math:`log(d!)`:
......@@ -200,24 +200,46 @@ class PoissonianEnergy(EnergyOperator):
class InverseGammaLikelihood(EnergyOperator):
"""
FIXME
"""Computes the negative log-likelihood of the inverse gamma distribution.
It negative log-pdf(x) is given by
.. math ::
\\sum_i (\\alpha_i+1)*\\ln(x_i) + \\beta_i/x_i
This is the likelihood for the variance :math:`x=S_k` given data
:math:`\\beta = 0.5 |s_k|^2` where the Field :math:`s` is known to have
the covariance :math:`S_k`.
Parameters
----------
beta : Field
beta parameter of the inverse gamma distribution
alpha : Scalar, Field, optional
alpha parameter of the inverse gamma distribution
"""
def __init__(self, d):
if not isinstance(d, Field):
def __init__(self, beta, alpha=-0.5):
if not isinstance(beta, Field):
raise TypeError
self._d = d
self._domain = DomainTuple.make(d.domain)
self._beta = beta
if np.isscalar(alpha):
alpha = Field.from_local_data(
beta.domain, np.full(beta.local_data.shape, alpha))
elif not isinstance(alpha, Field):
raise TypeError
self._alphap1 = alpha+1
self._domain = DomainTuple.make(beta.domain)
def apply(self, x):
self._check_input(x)
res = 0.5*(x.log().sum() + (1./x).vdot(self._d))
res = x.log().vdot(self._alphap1) + (1./x).vdot(self._beta)
if not isinstance(x, Linearization):
return Field.scalar(res)
if not x.want_metric:
return res
metric = SandwichOperator.make(x.jac, makeOp(0.5/(x.val**2)))
metric = SandwichOperator.make(x.jac, makeOp(self._alphap1/(x.val**2)))
return res.add_metric(metric)
......@@ -318,7 +340,7 @@ class StandardHamiltonian(EnergyOperator):
class AveragedEnergy(EnergyOperator):
"""Averages an energy over samples
"""Averages an energy over samples.
Parameters
----------
......@@ -328,15 +350,15 @@ class AveragedEnergy(EnergyOperator):
Set of residual sample points to be added to mean field for
approximate estimation of the KL.
Note
----
Having symmetrized residual samples, with both v_i and -v_i being
present, ensures that the distribution mean is exactly represented.
Notes
-----
- Having symmetrized residual samples, with both :math:`v_i` and
:math:`-v_i` being present, ensures that the distribution mean is
exactly represented.
:class:`AveragedEnergy(h)` approximates
:math:`\\left< H(f) \\right>_{G(f-m,D)}` if the residuals
:math:`f-m` are drawn from a Gaussian distribution with covariance
:math:`D`.
- :class:`AveragedEnergy(h)` approximates
:math:`\\left< H(f) \\right>_{G(f-m,D)}` if the residuals :math:`f-m`
are drawn from a Gaussian distribution with covariance :math:`D`.
"""
def __init__(self, h, res_samples):
......
......@@ -53,9 +53,8 @@ class LinearInterpolator(LinearOperator):
# FIXME This needs to be removed as soon as the bug below is fixed.
if dims.count(dims[0]) != len(dims):
raise TypeError(
'This is a bug. Please extend LinearInterpolators functionality!'
)
raise TypeError("This is a bug. Please extend"
"LinearInterpolator's functionality!")
shp = sampling_points.shape
if not (isinstance(sampling_points, np.ndarray) and len(shp) == 2):
......
......@@ -38,10 +38,10 @@ class SamplingEnabler(EndomorphicOperator):
The iteration controller to use for the iterative numerical inversion
done by a :class:`ConjugateGradient` object.
approximation : :class:`LinearOperator`, optional
if not None, this linear operator should be an approximation to the operator, which
supports the operation modes that the operator doesn't have. It is used as a
preconditioner during the iterative inversion, to accelerate
convergence.
if not None, this linear operator should be an approximation to the
operator, which supports the operation modes that the operator doesn't
have. It is used as a preconditioner during the iterative inversion,
to accelerate convergence.
"""
def __init__(self, likelihood, prior, iteration_controller,
......
......@@ -187,7 +187,7 @@ class GeometryRemover(LinearOperator):
domain: Domain, tuple of Domain, or DomainTuple:
the full input domain of the operator.
space: int, optional
The index of the subdomain on which the operator should act. Default is None.
The index of the subdomain on which the operator should act.
If None, it acts on all spaces.
Notes
......
......@@ -22,10 +22,9 @@ import numpy as np
from . import dobj
from .domains.gl_space import GLSpace
from .domains.hp_space import HPSpace
from .domains.log_rg_space import LogRGSpace
from .domains.power_space import PowerSpace
from .domains.rg_space import RGSpace
from .domains.log_rg_space import LogRGSpace
from .domain_tuple import DomainTuple
from .field import Field
# relevant properties:
......@@ -237,6 +236,8 @@ def _plot2D(f, ax, **kwargs):
foo = kwargs.pop("norm", None)
norm = {} if foo is None else {'norm': foo}