Commit 0fd92060 authored by Martin Reinecke's avatar Martin Reinecke

Merge remote-tracking branch 'origin/NIFTy_5' into simplify_for_const

parents 09f9e0bf 1a9b8429
......@@ -39,9 +39,9 @@ test_serial:
script:
- pytest-3 -q --cov=nifty5 test
- >
python3 -m coverage report --omit "*plotting*,*distributed_do*"
python3 -m coverage report --omit "*plot*,*distributed_do*"
- >
python3 -m coverage report --omit "*plotting*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }'
python3 -m coverage report --omit "*plot*,*distributed_do*" | grep TOTAL | awk '{ print "TOTAL: "$4; }'
test_mpi:
stage: test
......
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
......
......@@ -74,7 +74,7 @@ if __name__ == '__main__':
ic_sampling = ift.GradientNormController(iteration_limit=100)
# Minimize the Hamiltonian
H = ift.Hamiltonian(likelihood, ic_sampling)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
H = ift.EnergyAdapter(position, H, want_metric=True)
# minimizer = ift.L_BFGS(ic_newton)
H, convergence = minimizer(H)
......
......@@ -99,7 +99,7 @@ if __name__ == '__main__':
minimizer = ift.NewtonCG(ic_newton)
# Compute MAP solution by minimizing the information Hamiltonian
H = ift.Hamiltonian(likelihood)
H = ift.StandardHamiltonian(likelihood)
initial_position = ift.from_random('normal', domain)
H = ift.EnergyAdapter(initial_position, H, want_metric=True)
H, convergence = minimizer(H)
......
......@@ -100,10 +100,10 @@ if __name__ == '__main__':
# Set up likelihood and information Hamiltonian
likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
H = ift.Hamiltonian(likelihood, ic_sampling)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
initial_position = ift.MultiField.full(H.domain, 0.)
position = initial_position
initial_mean = ift.MultiField.full(H.domain, 0.)
mean = initial_mean
plot = ift.Plot()
plot.add(signal(mock_position), title='Ground Truth')
......@@ -117,9 +117,9 @@ if __name__ == '__main__':
# 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 = ift.MetricGaussianKL(mean, H, N_samples)
KL, convergence = minimizer(KL)
position = KL.position
mean = KL.position
# Plot current reconstruction
plot = ift.Plot()
......@@ -128,7 +128,7 @@ if __name__ == '__main__':
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)
KL = ift.MetricGaussianKL(mean, H, N_samples)
sc = ift.StatCalculator()
for sample in KL.samples:
sc.add(signal(sample + KL.position))
......
......@@ -103,7 +103,7 @@ N = ift.DiagonalOperator(ift.from_global_data(d_space, var))
IC = ift.DeltaEnergyController(tol_rel_deltaE=1e-12, iteration_limit=200)
likelihood = ift.GaussianEnergy(d, N)(R)
Ham = ift.Hamiltonian(likelihood, IC)
Ham = ift.StandardHamiltonian(likelihood, IC)
H = ift.EnergyAdapter(params, Ham, want_metric=True)
# Minimize
......
# rm -rf docs/build docs/source/mod
sphinx-apidoc -e -o docs/source/mod nifty5
sphinx-build -b html docs/source/ docs/build/
......@@ -264,13 +264,13 @@ This functionality is provided by NIFTy's
:class:`~inversion_enabler.InversionEnabler` class, which is itself a linear
operator.
.. currentmodule:: nifty5.operators.linear_operator
.. currentmodule:: nifty5.operators.operator
Direct multiplication and adjoint inverse multiplication transform a field
defined on the operator's :attr:`~LinearOperator.domain` to one defined on the
operator's :attr:`~LinearOperator.target`, whereas adjoint multiplication and
inverse multiplication transform from :attr:`~LinearOperator.target` to
:attr:`~LinearOperator.domain`.
defined on the operator's :attr:`~Operator.domain` to one defined on the
operator's :attr:`~Operator.target`, whereas adjoint multiplication and inverse
multiplication transform from :attr:`~Operator.target` to
:attr:`~Operator.domain`.
.. currentmodule:: nifty5.operators
......@@ -379,7 +379,7 @@ Minimization algorithms
All minimization algorithms in NIFTy inherit from the abstract
:class:`~minimizer.Minimizer` class, which presents a minimalistic interface
consisting only of a :meth:`~minimizer.Minimizer.__call__()` method taking an
consisting only of a :meth:`~minimizer.Minimizer.__call__` method taking an
:class:`~energy.Energy` object and optionally a preconditioning operator, and
returning the energy at the discovered minimum and a status code.
......@@ -399,17 +399,16 @@ Many minimizers for nonlinear problems can be characterized as
This family of algorithms is encapsulated in NIFTy's
:class:`~descent_minimizers.DescentMinimizer` class, which currently has three
concrete implementations: :class:`~descent_minimizers.SteepestDescent`,
:class:`~descent_minimizers.RelaxedNewton`,
:class:`~descent_minimizers.NewtonCG`, :class:`~descent_minimizers.L_BFGS` and
:class:`~descent_minimizers.VL_BFGS`. Of these algorithms, only
:class:`~descent_minimizers.RelaxedNewton` requires the energy object to provide
:class:`~descent_minimizers.NewtonCG` requires the energy object to provide
a :attr:`~energy.Energy.metric` property, the others only need energy values and
gradients.
The flexibility of NIFTy's design allows using externally provided minimizers.
With only small effort, adapters for two SciPy minimizers were written; they are
available under the names :class:`~scipy_minimizer.ScipyCG` and
:class:`L_BFGS_B`.
:class:`~scipy_minimizer.L_BFGS_B`.
Application to operator inversion
......@@ -432,4 +431,4 @@ performs a minimization of a
:class:`~minimization.quadratic_energy.QuadraticEnergy` by means of the
:class:`~minimization.conjugate_gradient.ConjugateGradient` algorithm. An
example is provided in
:func:`~ļibrary.wiener_filter_curvature.WienerFilterCurvature`.
:func:`~library.wiener_filter_curvature.WienerFilterCurvature`.
......@@ -13,6 +13,7 @@ napoleon_use_ivar = True
napoleon_use_admonition_for_notes = True
napoleon_use_admonition_for_examples = True
napoleon_use_admonition_for_references = True
napoleon_include_special_with_doc = True
project = u'NIFTy5'
copyright = u'2013-2019, Max-Planck-Society'
......@@ -27,3 +28,5 @@ add_module_names = False
html_theme = "sphinx_rtd_theme"
html_logo = 'nifty_logo_black.png'
exclude_patterns = ['mod/modules.rst', 'mod/*version.rst']
This diff is collapsed.
NIFTy -- Numerical Information Field Theory
===========================================
**NIFTy** [1]_, "\ **N**\umerical **I**\nformation **F**\ield **T**\heor\ **y**\ ", is a versatile library designed to enable the development of signal inference algorithms that are independent of the underlying spatial grid and its resolution.
**NIFTy** [1]_, [2]_, "\ **N**\umerical **I**\nformation **F**\ield **T**\heor\ **y**\ ", is a versatile library designed to enable the development of signal inference algorithms that are independent of the underlying spatial grid and its resolution.
Its object-oriented framework is written in Python, although it accesses libraries written in C++ and C for efficiency.
NIFTy offers a toolkit that abstracts discretized representations of continuous spaces, fields in these spaces, and operators acting on fields into classes.
......@@ -13,22 +13,18 @@ The set of spaces on which NIFTy operates comprises point sets, *n*-dimensional
References
----------
.. [1] Steininger et al., "NIFTy 3 - Numerical Information Field Theory - A Python framework for multicomponent signal inference on HPC clusters", 2017, submitted to PLOS One; `[arXiv:1708.01073] <https://arxiv.org/abs/1708.01073>`_
.. [1] Selig et al., "NIFTY - Numerical Information Field Theory. A versatile PYTHON library for signal inference ", 2013, Astronmy and Astrophysics 554, 26; `[DOI] <https://ui.adsabs.harvard.edu/link_gateway/2013A&A...554A..26S/doi:10.1051/0004-6361/201321236>`_, `[arXiv:1301.4499] <https://arxiv.org/abs/1301.4499>`_
.. [2] Steininger et al., "NIFTy 3 - Numerical Information Field Theory - A Python framework for multicomponent signal inference on HPC clusters", 2017, accepted by Annalen der Physik; `[arXiv:1708.01073] <https://arxiv.org/abs/1708.01073>`_
Contents
........
.. toctree::
:maxdepth: 2
ift
Gallery <http://wwwmpa.mpa-garching.mpg.de/ift/nifty/gallery/>
installation
code
citations
Indices and tables
..................
* :any:`Module Index <mod/modules>`
* :ref:`search`
Package Documentation <mod/nifty5>
......@@ -19,6 +19,7 @@ from .field import Field
from .multi_field import MultiField
from .operators.operator import Operator
from .operators.adder import Adder
from .operators.diagonal_operator import DiagonalOperator
from .operators.distributors import DOFDistributor, PowerDistributor
from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter
......@@ -33,7 +34,6 @@ from .operators.field_zero_padder import FieldZeroPadder
from .operators.inversion_enabler import InversionEnabler
from .operators.linear_operator import LinearOperator
from .operators.mask_operator import MaskOperator
from .operators.offset_operator import OffsetOperator
from .operators.qht_operator import QHTOperator
from .operators.regridding_operator import RegriddingOperator
from .operators.sampling_enabler import SamplingEnabler
......@@ -49,7 +49,7 @@ from .operators.simple_linear_operators import (
from .operators.value_inserter import ValueInserter
from .operators.energy_operators import (
EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
BernoulliEnergy, Hamiltonian, AveragedEnergy)
BernoulliEnergy, StandardHamiltonian, AveragedEnergy)
from .probing import probe_with_posterior_samples, probe_diagonal, \
StatCalculator
......@@ -64,16 +64,16 @@ 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
from .minimization.kl_energy import KL_Energy
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,
......
......@@ -307,7 +307,7 @@ for f in ["sqrt", "exp", "log", "tanh", "conjugate", "sin", "cos", "tan",
def clip(x, a_min=None, a_max=None):
return data_object(x.shape, np.clip(x.data, a_min, a_max), x.distaxis)
return data_object(x.shape, np.clip(x._data, a_min, a_max), x._distaxis)
def from_object(object, dtype, copy, set_locked):
......@@ -387,10 +387,14 @@ def distaxis(arr):
def from_local_data(shape, arr, distaxis=0):
if arr.dtype.kind not in "fciub":
raise TypeError
return data_object(shape, arr, distaxis)
def from_global_data(arr, sum_up=False, distaxis=0):
if arr.dtype.kind not in "fciub":
raise TypeError
if sum_up:
arr = np_allreduce_sum(arr)
if distaxis == -1:
......
......@@ -97,10 +97,14 @@ def distaxis(arr):
def from_local_data(shape, arr, distaxis=-1):
if tuple(shape) != arr.shape:
raise ValueError
if arr.dtype.kind not in "fciub":
raise TypeError
return arr
def from_global_data(arr, sum_up=False, distaxis=-1):
if arr.dtype.kind not in "fciub":
raise TypeError
return arr
......
......@@ -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)
......@@ -25,19 +25,19 @@ from .domain_tuple import DomainTuple
class Field(object):
"""The discrete representation of a continuous field over multiple spaces.
Stores data arrays and carries all the needed metainformation (i.e. the
Stores data arrays and carries all the needed meta-information (i.e. the
domain) for operators to be able to operate on them.
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 " +
......
......@@ -16,10 +16,10 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..minimization.energy_adapter import EnergyAdapter
from ..multi_domain import MultiDomain
from ..multi_field import MultiField
from ..operators.distributors import PowerDistributor
from ..operators.energy_operators import Hamiltonian, InverseGammaLikelihood
from ..operators.energy_operators import (StandardHamiltonian,
InverseGammaLikelihood)
from ..operators.scaling_operator import ScalingOperator
from ..operators.simple_linear_operators import ducktape
......@@ -35,25 +35,27 @@ def make_adjust_variances(a,
Constructs a Hamiltonian to solve constant likelihood optimizations of the
form phi = a * xi under the constraint that phi remains constant.
FIXME xi is white.
Parameters
----------
a : Operator
Operator which gives the amplitude when evaluated at a position
Gives the amplitude when evaluated at a position.
xi : Operator
Operator which gives the excitation when evaluated at a position
postion : Field, MultiField
Position of the whole problem
Gives the excitation when evaluated at a position.
position : Field, MultiField
Position of the entire problem.
samples : Field, MultiField
Residual samples of the whole problem
Residual samples of the whole problem.
scaling : Float
Optional rescaling of the Likelihood
Optional rescaling of the Likelihood.
ic_samp : Controller
Iteration Controller for Hamiltonian
Iteration Controller for Hamiltonian.
Returns
-------
Hamiltonian
A Hamiltonian that can be used for further minimization
StandardHamiltonian
A Hamiltonian that can be used for further minimization.
"""
d = a*xi
......@@ -71,7 +73,8 @@ def make_adjust_variances(a,
if scaling is not None:
x = ScalingOperator(scaling, x.target)(x)
return Hamiltonian(InverseGammaLikelihood(d_eval)(x), ic_samp=ic_samp)
return StandardHamiltonian(InverseGammaLikelihood(d_eval/2.)(x),
ic_samp=ic_samp)
def do_adjust_variances(position,
......@@ -79,6 +82,9 @@ def do_adjust_variances(position,
minimizer,
xi_key='xi',
samples=[]):
'''
FIXME
'''
h_space = position[xi_key].domain[0]
pd = PowerDistributor(h_space, amplitude_operator.target[0])
......
......@@ -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
......@@ -24,16 +25,16 @@ from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import ducktape
def CorrelatedField(target, amplitude_operator, name='xi'):
'''Constructs an operator which turns a white Gaussian excitation field
def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
"""Constructs an operator which turns a white Gaussian excitation field
into a correlated field.
This function returns an operator which implements:
ht @ (vol * A * xi),
where `ht` is a harmonic transform operator, `A` is the sqare root of the
prior covariance an `xi` is the excitation field.
where `ht` is a harmonic transform operator, `A` is the square root of the
prior covariance and `xi` is the excitation field.
Parameters
----------
......@@ -41,26 +42,35 @@ def CorrelatedField(target, amplitude_operator, name='xi'):
Target of the operator. Must contain exactly one space.
amplitude_operator: Operator
name : string
:class:`MultiField` key for xi-field.
:class:`MultiField` key for the xi-field.
codomain : Domain
The codomain for target[0]. If not supplied, it is inferred.
Returns
-------
Correlated field : Operator
'''
Operator
Correlated field
"""
tgt = DomainTuple.make(target)
if len(tgt) > 1:
raise ValueError
h_space = tgt[0].get_default_codomain()
ht = HarmonicTransformOperator(h_space, tgt[0])
if codomain is None:
codomain = tgt[0].get_default_codomain()
h_space = codomain
ht = HarmonicTransformOperator(h_space, target=tgt[0])
p_space = amplitude_operator.target[0]
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))
def MfCorrelatedField(target, amplitudes, name='xi'):
'''Constructs an operator which turns white Gaussian excitation fields
"""Constructs an operator which turns white Gaussian excitation fields
into a correlated field defined on a DomainTuple with two entries and two
separate correlation structures.
......@@ -70,7 +80,7 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
Parameters
----------
target : Domain, DomainTuple or tuple of Domain
Target of the operator. Must contain exactly one space.
Target of the operator. Must contain exactly two spaces.
amplitudes: iterable of Operator
List of two amplitude operators.
name : string
......@@ -78,8 +88,9 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
Returns
-------
Correlated field : Operator
'''
Operator
Correlated field
"""
tgt = DomainTuple.make(target)