Commit 9f31d6a2 authored by Philipp Arras's avatar Philipp Arras

Merge remote-tracking branch 'dev/NIFTy_5' into sebastian_review

parents 279a9d51 4f8fa350
......@@ -39,9 +39,9 @@ Installation
- [Python 3](https://www.python.org/) (3.5.x or later)
- [SciPy](https://www.scipy.org/)
- [pyFFTW](https://pypi.python.org/pypi/pyFFTW)
Optional dependencies:
- [pyFFTW](https://pypi.python.org/pypi/pyFFTW) for faster Fourier transforms
- [pyHealpix](https://gitlab.mpcdf.mpg.de/ift/pyHealpix) (for harmonic
transforms involving domains on the sphere)
- [mpi4py](https://mpi4py.scipy.org) (for MPI-parallel execution)
......@@ -61,18 +61,29 @@ distributions, the "apt" lines will need slight changes.
NIFTy5 and its mandatory dependencies can be installed via:
sudo apt-get install git libfftw3-dev 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_5
Plotting support is added via:
pip3 install --user matplotlib
FFTW support is added via:
sudo apt-get install libfftw3-dev
pip3 install --user pyfftw
To actually use FFTW in your Nifty calculations, you need to call
nifty5.fft.enable_fftw()
at the beginning of your code.
(Note: If you encounter problems related to `pyFFTW`, make sure that you are
using a pip-installed `pyFFTW` package. Unfortunately, some distributions are
shipping an incorrectly configured `pyFFTW` package, which does not cooperate
with the installed `FFTW3` libraries.)
Plotting support is added via:
pip3 install --user matplotlib
Support for spherical harmonic transforms is added via:
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
......@@ -86,7 +97,7 @@ MPI support is added via:
To run the tests, additional packages are required:
sudo apt-get install python3-coverage python3-parameterized python3-pytest python3-pytest-cov
sudo apt-get install python3-coverage python3-pytest python3-pytest-cov
Afterwards the tests (including a coverage report) can be run using the
following command in the repository root:
......
......@@ -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']
......@@ -9,19 +9,20 @@ Theoretical Background
IFT is fully Bayesian. How else could infinitely many field degrees of freedom be constrained by finite 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).
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, renormalization, and more. IFT reproduces many known well working algorithms. This should be reassuring. Also, 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). (*FIXME* does it?)
.. tip:: *In-a-nutshell introductions to information field theory* can be found in [2]_, [3]_, [4]_, and [5]_, with the latter probably being the most didactically.
.. tip:: *In-a-nutshell introductions to information field theory* can be found in [2]_, [3]_, [4]_, and [5]_, with the latter probably being the most didactical.
.. [1] T.A. Enßlin et al. (2009), "Information field theory for cosmological perturbation reconstruction and nonlinear signal analysis", PhysRevD.80.105005, 09/2009; `arXiv:0806.3474 <http://www.arxiv.org/abs/0806.3474>`_
.. [1] T.A. Enßlin et al. (2009), "Information field theory for cosmological perturbation reconstruction and nonlinear signal analysis", PhysRevD.80.105005, 09/2009; `[arXiv:0806.3474] <http://www.arxiv.org/abs/0806.3474>`_
.. [2] T.A. Enßlin (2013), "Information field theory", proceedings of MaxEnt 2012 -- the 32nd International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering; AIP Conference Proceedings, Volume 1553, Issue 1, p.184; `arXiv:1301.2556 <http://arxiv.org/abs/1301.2556>`_
.. [2] T.A. Enßlin (2013), "Information field theory", proceedings of MaxEnt 2012 -- the 32nd International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering; AIP Conference Proceedings, Volume 1553, Issue 1, p.184; `[arXiv:1301.2556] <http://arxiv.org/abs/1301.2556>`_
.. [3] T.A. Enßlin (2014), "Astrophysical data analysis with information field theory", AIP Conference Proceedings, Volume 1636, Issue 1, p.49; `arXiv:1405.7701 <http://arxiv.org/abs/1405.7701>`_
.. [3] T.A. Enßlin (2014), "Astrophysical data analysis with information field theory", AIP Conference Proceedings, Volume 1636, Issue 1, p.49; `[arXiv:1405.7701] <http://arxiv.org/abs/1405.7701>`_
.. [4] Wikipedia contributors (2018), `"Information field theory" <https://en.wikipedia.org/w/index.php?title=Information_field_theory&oldid=876731720>`_, Wikipedia, The Free Encyclopedia.
.. [5] T.A. Enßlin (2019), "Information theory for fields", accepted by Annalen der Physik; `arXiv:1804.03350 <http://arxiv.org/abs/1804.03350>`_
.. [5] T.A. Enßlin (2019), "Information theory for fields", accepted by Annalen der Physik; `[DOI] <https://doi.org/10.1002/andp.201800127>`_, `[arXiv:1804.03350] <http://arxiv.org/abs/1804.03350>`_
Discretized continuum
---------------------
......@@ -83,7 +84,7 @@ The above line of argumentation analogously applies to the discretization of ope
A(x \in \Omega_p, y \in \Omega_q) \quad\mapsto\quad A_{pq} \quad=\quad \frac{\iint_{\Omega_p \Omega_q} \mathrm{d}x \, \mathrm{d}y \; A(x,y)}{\iint_{\Omega_p \Omega_q} \mathrm{d}x \, \mathrm{d}y} \quad=\quad \big< \big< A(x,y) \big>_{\Omega_p} \big>_{\Omega_q}
.
The proper discretization of spaces, fields, and operators, as well as the normalization of position integrals, is essential for the conservation of the continuum limit. Their consistent implementation in NIFTY allows a pixelization independent coding of algorithms.
The proper discretization of spaces, fields, and operators, as well as the normalization of position integrals, is essential for the conservation of the continuum limit. Their consistent implementation in NIFTy allows a pixelization independent coding of algorithms.
Free Theory & Implicit Operators
--------------------------------
......@@ -100,10 +101,9 @@ and the measurement equation is linear in both signal and noise,
d= R\, s + n,
with :math:`{R}` the measurement response, which maps the continous signal field into the discrete data space.
with :math:`{R}` being the measurement response, which maps the continous signal field into the discrete data space.
This is called a free theory, as the information Hamiltonian
associate professor
.. math::
......@@ -135,9 +135,9 @@ the posterior covariance operator, and
j = R^\dagger N^{-1} d
the information source. The operation in :math:`{d= D\,R^\dagger N^{-1} d}` is also called the generalized Wiener filter.
the information source. The operation in :math:`{m = D\,R^\dagger N^{-1} d}` is also called the generalized Wiener filter.
NIFTy permits to define the involved operators :math:`{R}`, :math:`{R^\dagger}`, :math:`{S}`, and :math:`{N}` implicitely, as routines that can be applied to vectors, but which do not require the explicit storage of the matrix elements of the operators.
NIFTy permits to define the involved operators :math:`{R}`, :math:`{R^\dagger}`, :math:`{S}`, and :math:`{N}` implicitly, as routines that can be applied to vectors, but which do not require the explicit storage of the matrix elements of the operators.
Some of these operators are diagonal in harmonic (Fourier) basis, and therefore only require the specification of a (power) spectrum and :math:`{S= F\,\widehat{P_s} F^\dagger}`. Here :math:`{F = \mathrm{HarmonicTransformOperator}}`, :math:`{\widehat{P_s} = \mathrm{DiagonalOperator}(P_s)}`, and :math:`{P_s(k)}` is the power spectrum of the process that generated :math:`{s}` as a function of the (absolute value of the) harmonic (Fourier) space koordinate :math:`{k}`. For those, NIFTy can easily also provide inverse operators, as :math:`{S^{-1}= F\,\widehat{\frac{1}{P_s}} F^\dagger}` in case :math:`{F}` is unitary, :math:`{F^\dagger=F^{-1}}`.
......@@ -175,10 +175,10 @@ NIFTy takes advantage of this formulation in several ways:
1) All prior degrees of freedom have unit covariance which improves the condition number of operators which need to be inverted.
2) The amplitude operator can be regarded as part of the response, :math:`{R'=R\,A}`. In general, more sophisticated responses can be constructed out of the composition of simpler operators.
3) The response can be non-linear, e.g. :math:`{R'(s)=R \exp(A\,\xi)}`, see demos/getting_started_2.py.
4) The amplitude operator can be made dependent on unknowns as well, e.g. :math:`A=A(\tau)= F\, \widehat{e^\tau}` represents an amplitude operator with a positive definite, unknown spectrum defined in Fourier domain. The amplitude field :math:`{\tau}` would get its own amplitude operator, with a cepstrum (spectrum of a log spectrum) defined in quefrency space (harmonic space of a logarithmically binned harmonic space) to regularize its degrees of freedom by imposing some (user-defined level of) spectral smoothness.
5) NIFTy can calculate the gradient of the information Hamiltonian and the Fischer information metric with respect to all unknown parameters, here :math:`{\xi}` and :math:`{\tau}`, by automatic differentiation. The gradients are used for MAP and HMCF estimates, and the Fischer matrix is required in addition to the gradient by Metric Gaussian Variational Inference (MGVI), which is available in NIFTy as well. MGVI is an implicit operator extension of Automatic Differentiation Variational Inference (ADVI).
4) The amplitude operator can be made dependent on unknowns as well, e.g. :math:`A=A(\tau)= F\, \widehat{e^\tau}` represents an amplitude operator with a positive definite, unknown spectrum defined in the Fourier domain. The amplitude field :math:`{\tau}` would get its own amplitude operator, with a cepstrum (spectrum of a log spectrum) defined in quefrency space (harmonic space of a logarithmically binned harmonic space) to regularize its degrees of freedom by imposing some (user-defined degree of) spectral smoothness.
5) NIFTy can calculate the gradient of the information Hamiltonian and the Fisher information metric with respect to all unknown parameters, here :math:`{\xi}` and :math:`{\tau}`, by automatic differentiation. The gradients are used for MAP and HMCF estimates, and the Fisher matrix is required in addition to the gradient by Metric Gaussian Variational Inference (MGVI), which is available in NIFTy as well. MGVI is an implicit operator extension of Automatic Differentiation Variational Inference (ADVI).
The reconstruction of a non-Gaussian signal with unknown covarinance from a non-trivial (tomographic) response is demonstrated in demos/getting_started_3.py. Here, the uncertainty of the field and the power spectrum of its generating process are probed via posterior samples provided by the MGVI algorithm.
The reconstruction of a non-Gaussian signal with unknown covariance from a non-trivial (tomographic) response is demonstrated in demos/getting_started_3.py. Here, the uncertainty of the field and the power spectrum of its generating process are probed via posterior samples provided by the MGVI algorithm.
+----------------------------------------------------+
| **Output of tomography demo getting_started_3.py** |
......
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>
......@@ -7,18 +7,29 @@ distributions, the "apt" lines will need slight changes.
NIFTy5 and its mandatory dependencies can be installed via::
sudo apt-get install git libfftw3-dev 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_5
Plotting support is added via::
pip3 install --user matplotlib
FFTW support is added via:
sudo apt-get install libfftw3-dev
pip3 install --user pyfftw
To actually use FFTW in your Nifty calculations, you need to call
nifty5.fft.enable_fftw()
at the beginning of your code.
(Note: If you encounter problems related to `pyFFTW`, make sure that you are
using a pip-installed `pyFFTW` package. Unfortunately, some distributions are
shipping an incorrectly configured `pyFFTW` package, which does not cooperate
with the installed `FFTW3` libraries.)
Plotting support is added via::
pip3 install --user matplotlib
Support for spherical harmonic transforms is added via::
pip3 install --user git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
......
......@@ -16,10 +16,10 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
from ..utilities import NiftyMetaBase
from ..utilities import NiftyMeta
class Domain(NiftyMetaBase()):
class Domain(metaclass=NiftyMeta):
"""The abstract class repesenting a (structured or unstructured) domain.
"""
def __repr__(self):
......
......@@ -19,23 +19,63 @@ from .utilities import iscomplextype
import numpy as np
_use_fftw = True
_use_fftw = False
_fftw_prepped = False
_fft_extra_args = {}
if _use_fftw:
import pyfftw
from pyfftw.interfaces.numpy_fft import fftn, rfftn, ifftn
pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(1000.)
# Optional extra arguments for the FFT calls
# if exact reproducibility is needed,
# set "planner_effort" to "FFTW_ESTIMATE"
import os
nthreads = int(os.getenv("OMP_NUM_THREADS", "1"))
_fft_extra_args = dict(planner_effort='FFTW_ESTIMATE', threads=nthreads)
else:
from numpy.fft import fftn, rfftn, ifftn
_fft_extra_args = {}
def enable_fftw():
global _use_fftw
_use_fftw = True
def disable_fftw():
global _use_fftw
_use_fftw = False
def _init_pyfftw():
global _fft_extra_args, _fftw_prepped
if not _fftw_prepped:
import pyfftw
from pyfftw.interfaces.numpy_fft import fftn, rfftn, ifftn
pyfftw.interfaces.cache.enable()
pyfftw.interfaces.cache.set_keepalive_time(1000.)
# Optional extra arguments for the FFT calls
# if exact reproducibility is needed,
# set "planner_effort" to "FFTW_ESTIMATE"
import os
nthreads = int(os.getenv("OMP_NUM_THREADS", "1"))
_fft_extra_args = dict(planner_effort='FFTW_ESTIMATE',
threads=nthreads)
_fftw_prepped = True
def fftn(a, axes=None):
if _use_fftw:
from pyfftw.interfaces.numpy_fft import fftn
_init_pyfftw()
return fftn(a, axes=axes, **_fft_extra_args)
else:
return np.fft.fftn(a, axes=axes)
def rfftn(a, axes=None):
if _use_fftw:
from pyfftw.interfaces.numpy_fft import rfftn
_init_pyfftw()
return rfftn(a, axes=axes, **_fft_extra_args)
else:
return np.fft.rfftn(a, axes=axes)
def ifftn(a, axes=None):
if _use_fftw:
from pyfftw.interfaces.numpy_fft import ifftn
_init_pyfftw()
return ifftn(a, axes=axes, **_fft_extra_args)
else:
return np.fft.ifftn(a, axes=axes)
def hartley(a, axes=None):
......@@ -46,7 +86,7 @@ def hartley(a, axes=None):
if iscomplextype(a.dtype):
raise TypeError("Hartley transform requires real-valued arrays.")
tmp = rfftn(a, axes=axes, **_fft_extra_args)
tmp = rfftn(a, axes=axes)
def _fill_array(tmp, res, axes):
if axes is None:
......@@ -89,7 +129,7 @@ def my_fftn_r2c(a, axes=None):
if iscomplextype(a.dtype):
raise TypeError("Transform requires real-valued input arrays.")
tmp = rfftn(a, axes=axes, **_fft_extra_args)
tmp = rfftn(a, axes=axes)
def _fill_complex_array(tmp, res, axes):
if axes is None:
......@@ -123,4 +163,4 @@ def my_fftn_r2c(a, axes=None):
def my_fftn(a, axes=None):
return fftn(a, axes=axes, **_fft_extra_args)
return fftn(a, axes=axes)
......@@ -43,7 +43,7 @@ def CepstrumOperator(target, a, k0):
and ceps is the so-called cepstrum:
.. math::
\\mathrm{sqrt\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
\\mathrm{sqrt\\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
These operators are combined in this fashion in order to generate:
......
......@@ -461,7 +461,7 @@ class Linearization(object):
if len(constants) == 0:
return Linearization.make_var(field, want_metric)
else:
ops = [ScalingOperator(0. if key in constants else 1., dom)
for key, dom in field.domain.items()]
bdop = BlockDiagonalOperator(field.domain, tuple(ops))
ops = {key: ScalingOperator(0. if key in constants else 1., dom)
for key, dom in field.domain.items()}
bdop = BlockDiagonalOperator(field.domain, ops)
return Linearization(field, bdop, want_metric=want_metric)
......@@ -15,10 +15,10 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..utilities import NiftyMetaBase
from ..utilities import NiftyMeta
class Energy(NiftyMetaBase()):
class Energy(metaclass=NiftyMeta):
"""Provides the functional used by minimization schemes.
The Energy object is an implementation of a scalar function including its
......
......@@ -16,11 +16,11 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..logger import logger
from ..utilities import NiftyMetaBase
from ..utilities import NiftyMeta
import numpy as np
class IterationController(NiftyMetaBase()):
class IterationController(metaclass=NiftyMeta):
"""The abstract base class for all iteration controllers.
An iteration controller is an object that monitors the progress of a
minimization iteration. At the begin of the minimization, its start()
......
......@@ -15,7 +15,7 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..utilities import NiftyMetaBase
from ..utilities import NiftyMeta
from ..logger import logger
import numpy as np
......@@ -103,7 +103,7 @@ class LineEnergy(object):
return res.real
class LineSearch(NiftyMetaBase()):
class LineSearch(metaclass=NiftyMeta):
"""Class for finding a step size that satisfies the strong Wolfe
conditions.
......
......@@ -15,13 +15,12 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..utilities import NiftyMetaBase
from ..utilities import NiftyMeta
class Minimizer(NiftyMetaBase()):
class Minimizer(metaclass=NiftyMeta):
"""A base class used by all minimizers."""
# MR FIXME: the docstring is partially ignored by Sphinx. Why?
def __call__(self, energy, preconditioner=None):
"""Performs the minimization of the provided Energy functional.
......
......@@ -24,18 +24,17 @@ class BlockDiagonalOperator(EndomorphicOperator):
"""
Parameters
----------
domain : :class:`MultiDomain`
domain : MultiDomain
Domain and target of the operator.
operators : dict
Dictionary with subdomain names as keys and :class:`LinearOperators` as items.
Dictionary with subdomain names as keys and :class:`LinearOperator`s
as items.
"""
def __init__(self, domain, operators):
if not isinstance(domain, MultiDomain):
raise TypeError("MultiDomain expected")
if not isinstance(operators, tuple):
raise TypeError("tuple expected")
self._domain = domain
self._ops = operators
self._ops = tuple(operators[key] for key in domain.keys())
self._capability = self._all_ops
for op in self._ops:
if op is not None:
......@@ -56,13 +55,14 @@ class BlockDiagonalOperator(EndomorphicOperator):
def _combine_chain(self, op):
if self._domain != op._domain:
raise ValueError("domain mismatch")
res = tuple(v1(v2) for v1, v2 in zip(self._ops, op._ops))
res = {key: v1(v2)
for key, v1, v2 in zip(self._domain.keys(), self._ops, op._ops)}
return BlockDiagonalOperator(self._domain, res)
def _combine_sum(self, op, selfneg, opneg):
from ..operators.sum_operator import SumOperator
if self._domain != op._domain:
raise ValueError("domain mismatch")
res = tuple(SumOperator.make([v1, v2], [selfneg, opneg])
for v1, v2 in zip(self._ops, op._ops))
res = {key: SumOperator.make([v1, v2], [selfneg, opneg])
for key, v1, v2 in zip(self._domain.keys(), self._ops, op._ops)}
return BlockDiagonalOperator(self._domain, res)
......@@ -318,14 +318,14 @@ class Hamiltonian(EnergyOperator):
class AveragedEnergy(EnergyOperator):
"""Computes Kullbach-Leibler (KL) divergence or Gibbs free energies.
"""Computes Kullback-Leibler (KL) divergence or Gibbs free energies.
A sample-averaged energy, e.g. an Hamiltonian, approximates the relevant
part of a KL to be used in Variational Bayes inference if the samples are
drawn from the approximating Gaussian:
.. math ::
\\text{KL}(m) = \\frac1{\\#\{v_i\}} \\sum_{v_i} H(m+v_i),
\\text{KL}(m) = \\frac1{\\#\\{v_i\\}} \\sum_{v_i} H(m+v_i),
where :math:`v_i` are the residual samples and :math:`m` is the mean field
around which the samples are drawn.
......
......@@ -16,10 +16,10 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from ..utilities import NiftyMetaBase, indent
from ..utilities import NiftyMeta, indent
class Operator(NiftyMetaBase()):
class Operator(metaclass=NiftyMeta):
"""Transforms values defined on one domain into values defined on another
domain, and can also provide the Jacobian.
"""
......
......@@ -15,6 +15,9 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
from operator import mul
import numpy as np
from ..domain_tuple import DomainTuple
......@@ -25,7 +28,7 @@ from .linear_operator import LinearOperator
class ValueInserter(LinearOperator):
"""Inserts one value into a field which is zero otherwise.
"""Inserts one value into a field which is constant otherwise.
Parameters
----------
......@@ -33,25 +36,37 @@ class ValueInserter(LinearOperator):
index : iterable of int
The index of the target into which the value of the domain shall be
inserted.
default_value : float
Constant value which is inserted everywhere where the input operator
is not inserted. Default is 0.
"""
def __init__(self, target, index):
def __init__(self, target, index, default_value=0.):
self._domain = makeDomain(UnstructuredDomain(1))
self._target = DomainTuple.make(target)
# Type and value checks
index = tuple(index)
if not all([isinstance(n, int) and n>=0 and n<self.target.shape[i] for i, n in enumerate(index)]):
if not all([
isinstance(n, int) and n >= 0 and n < self.target.shape[i]
for i, n in enumerate(index)
]):
raise TypeError
if not len(index) == len(self.target.shape):
raise ValueError
np.empty(self.target.shape)[index]
self._index = index
self._dv = float(default_value)
self._dvsum = self._dv*(reduce(mul, self.target.shape) - 1)
self._capability = self.TIMES | self.ADJOINT_TIMES
# Check whether index is in bounds
np.empty(self.target.shape)[self._index]
def apply(self, x, mode):
self._check_input(x, mode)
x = x.to_global_data()
if mode == self.TIMES:
res = np.zeros(self.target.shape, dtype=x.dtype)
res = np.full(self.target.shape, self._dv, dtype=x.dtype)
res[self._index] = x
else:
res = np.full((1,), x[self._index], dtype=x.dtype)
res = np.full((1,), x[self._index] + self._dvsum, dtype=x.dtype)
return Field.from_global_data(self._tgt(mode), res)
......@@ -16,6 +16,8 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from .field import Field
from .operators.endomorphic_operator import EndomorphicOperator
from .operators.operator import Operator
class StatCalculator(object):
......@@ -69,6 +71,29 @@ class StatCalculator(object):
def probe_with_posterior_samples(op, post_op, nprobes):
'''FIXME
Parameters
----------
op : EndomorphicOperator
FIXME
post_op : Operator
FIXME
nprobes : int
Number of samples which shall be drawn.