Commit b760b9dd authored by Martin Reinecke's avatar Martin Reinecke

merge NIFTy_5

parents c9cf46c9 d85c3e6f
NIFTy - Numerical Information Field Theory
==========================================
[![build status](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/NIFTy_5/build.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/NIFTy_5)
[![coverage report](https://gitlab.mpcdf.mpg.de/ift/NIFTy/badges/NIFTy_5/coverage.svg)](https://gitlab.mpcdf.mpg.de/ift/NIFTy/commits/NIFTy_5)
[![build status](https://gitlab.mpcdf.mpg.de/ift/nifty-dev/badges/NIFTy_5/build.svg)](https://gitlab.mpcdf.mpg.de/ift/nifty-dev/commits/NIFTy_5)
[![coverage report](https://gitlab.mpcdf.mpg.de/ift/nifty-dev/badges/NIFTy_5/coverage.svg)](https://gitlab.mpcdf.mpg.de/ift/nifty-dev/commits/NIFTy_5)
**NIFTy** project homepage:
[http://ift.pages.mpcdf.de/NIFTy](http://ift.pages.mpcdf.de/NIFTy)
......@@ -108,7 +108,7 @@ For a quick start, you can browse through the [informal
introduction](http://ift.pages.mpcdf.de/NIFTy/code.html) or
dive into NIFTy by running one of the demonstrations, e.g.:
python demos/wiener_filter_via_curvature.py
python demos/getting_started_1.py
### Acknowledgement
......
......@@ -74,7 +74,7 @@ if __name__ == '__main__':
N_samples = 20
for i in range(5):
H = H.at(position)
samples = [H.curvature.draw_sample(from_inverse=True)
samples = [H.metric.draw_sample(from_inverse=True)
for _ in range(N_samples)]
KL = ift.SampledKullbachLeiblerDivergence(H, samples)
......
......@@ -306,16 +306,16 @@ Energy functionals
In NIFTy5 such functions are represented by objects of type :class:`Energy`.
These hold the prescription how to calculate the function's
:attr:`~Energy.value`, :attr:`~Energy.gradient` and
(optionally) :attr:`~Energy.curvature` at any given :attr:`~Energy.position`
(optionally) :attr:`~Energy.metric` at any given :attr:`~Energy.position`
in parameter space.
Function values are floating-point scalars, gradients have the form of fields
living on the energy's position domain, and curvatures are represented by
living on the energy's position domain, and metrics are represented by
linear operator objects.
Energies are classes that typically have to be provided by the user when
tackling new IFT problems.
Some examples of concrete energy classes delivered with NIFTy5 are
:class:`QuadraticEnergy` (with position-independent curvature, mainly used with
:class:`QuadraticEnergy` (with position-independent metric, mainly used with
conjugate gradient minimization) and :class:`~nifty5.library.WienerFilterEnergy`.
......@@ -367,7 +367,7 @@ This family of algorithms is encapsulated in NIFTy's :class:`DescentMinimizer`
class, which currently has three concrete implementations:
:class:`SteepestDescent`, :class:`VL_BFGS`, and :class:`RelaxedNewton`.
Of these algorithms, only :class:`RelaxedNewton` requires the energy object to
provide a :attr:`~Energy.curvature` property, the others only need energy
provide a :attr:`~Energy.metric` property, the others only need energy
values and gradients.
The flexibility of NIFTy's design allows using externally provided
......
......@@ -21,7 +21,7 @@ from .nonlinearities import Exponential, Linear, PositiveTanh, Tanh
from .models.constant import Constant
from .models.linear_model import LinearModel
from .models.local_nonlinearity import (LocalModel, PointwiseExponential,
PointwisePositiveTanh, PointwiseTanh)
PointwisePositiveTanh, PointwiseTanh)
from .models.model import Model
from .models.multi_model import MultiModel
from .models.variable import Variable
......@@ -65,7 +65,8 @@ from .minimization.steepest_descent import SteepestDescent
from .minimization.vl_bfgs import VL_BFGS
from .minimization.l_bfgs import L_BFGS
from .minimization.relaxed_newton import RelaxedNewton
from .minimization.scipy_minimizer import ScipyMinimizer, NewtonCG, L_BFGS_B, ScipyCG
from .minimization.scipy_minimizer import (ScipyMinimizer, NewtonCG, L_BFGS_B,
ScipyCG)
from .minimization.energy import Energy
from .minimization.quadratic_energy import QuadraticEnergy
from .minimization.line_energy import LineEnergy
......@@ -82,7 +83,8 @@ from .library.point_sources import PointSources
from .library.poissonian_energy import PoissonianEnergy
from .library.wiener_filter_curvature import WienerFilterCurvature
from .library.wiener_filter_energy import WienerFilterEnergy
from .library.correlated_fields import make_correlated_field, make_mf_correlated_field
from .library.correlated_fields import (make_correlated_field,
make_mf_correlated_field)
from .library.bernoulli_energy import BernoulliEnergy
from . import extra
......
......@@ -49,13 +49,13 @@ class Hamiltonian(Energy):
@property
@memo
def curvature(self):
prior_curv = self._prior.curvature
def metric(self):
prior_mtr = self._prior.metric
if self._ic_samp is None:
return self._lh.curvature + prior_curv
return self._lh.metric + prior_mtr
else:
return SamplingEnabler(self._lh.curvature, prior_curv.inverse,
self._ic_samp, prior_curv.inverse)
return SamplingEnabler(self._lh.metric, prior_mtr.inverse,
self._ic_samp, prior_mtr.inverse)
def __str__(self):
res = 'Likelihood:\t{:.2E}\n'.format(self._lh.value)
......
......@@ -33,6 +33,6 @@ class SampledKullbachLeiblerDivergence(Energy):
@property
@memo
def curvature(self):
return (my_sum(map(lambda v: v.curvature, self._energy_list)) *
(1./len(self._energy_list)))
def metric(self):
return (my_sum(map(lambda v: v.metric, self._energy_list)) *
(1./len(self._energy_list)))
......@@ -22,7 +22,7 @@ from ..minimization.energy import Energy
from ..models.model import Model
__all__ = ["check_value_gradient_consistency",
"check_value_gradient_curvature_consistency"]
"check_value_gradient_metric_consistency"]
def _get_acceptable_model(M):
......@@ -30,7 +30,7 @@ def _get_acceptable_model(M):
if not np.isfinite(val.sum()):
raise ValueError('Initial Model value must be finite')
dir = from_random("normal", M.position.domain)
dirder = M.gradient(dir)
dirder = M.jacobian(dir)
dir *= val/(dirder).norm()*1e-5
# Find a step length that leads to a "reasonable" Model
for i in range(50):
......@@ -82,7 +82,7 @@ def check_value_gradient_consistency(E, tol=1e-8, ntries=100):
if isinstance(E, Energy):
dirder = Emid.gradient.vdot(dir)/dirnorm
else:
dirder = Emid.gradient(dir)/dirnorm
dirder = Emid.jacobian(dir)/dirnorm
numgrad = (E2.value-val)/dirnorm
if isinstance(E, Model):
xtol = tol * dirder.norm() / np.sqrt(dirder.size)
......@@ -100,9 +100,9 @@ def check_value_gradient_consistency(E, tol=1e-8, ntries=100):
E = Enext
def check_value_gradient_curvature_consistency(E, tol=1e-8, ntries=100):
def check_value_gradient_metric_consistency(E, tol=1e-8, ntries=100):
if isinstance(E, Model):
raise ValueError('Models have no curvature, thus it cannot be tested.')
raise ValueError('Models have no metric, thus it cannot be tested.')
for _ in range(ntries):
E2 = _get_acceptable_energy(E)
val = E.value
......@@ -112,7 +112,7 @@ def check_value_gradient_curvature_consistency(E, tol=1e-8, ntries=100):
for i in range(50):
Emid = E.at(E.position + 0.5*dir)
dirder = Emid.gradient.vdot(dir)/dirnorm
dgrad = Emid.curvature(dir)/dirnorm
dgrad = Emid.metric(dir)/dirnorm
xtol = tol*Emid.gradient_norm
if abs((E2.value-val)/dirnorm - dirder) < xtol and \
(abs((E2.gradient-E.gradient)/dirnorm-dgrad) < xtol).all():
......@@ -121,5 +121,5 @@ def check_value_gradient_curvature_consistency(E, tol=1e-8, ntries=100):
dirnorm *= 0.5
E2 = Emid
else:
raise ValueError("gradient, value and curvature seem inconsistent")
raise ValueError("gradient, value and metric seem inconsistent")
E = Enext
......@@ -39,9 +39,9 @@ class BernoulliEnergy(Energy):
if isnan(self._value):
self._value = inf
metric = makeOp(1./((p_val) * (1.-p_val)))
self._gradient = self._p.gradient.adjoint_times(metric(p_val-d))
self._gradient = self._p.jacobian.adjoint_times(metric(p_val-d))
self._curvature = SandwichOperator.make(self._p.gradient, metric)
self._metric = SandwichOperator.make(self._p.jacobian, metric)
def at(self, position):
return self.__class__(self._p.at(position), self._d)
......@@ -55,5 +55,5 @@ class BernoulliEnergy(Energy):
return self._gradient
@property
def curvature(self):
return self._curvature
def metric(self):
return self._metric
......@@ -47,7 +47,8 @@ def make_mf_correlated_field(s_space_spatial, s_space_energy,
from ..models.variable import Variable
from ..domain_tuple import DomainTuple
from ..operators.domain_distributor import DomainDistributor
from ..operators.harmonic_transform_operator import HarmonicTransformOperator
from ..operators.harmonic_transform_operator \
import HarmonicTransformOperator
h_space_spatial = s_space_spatial.get_default_codomain()
h_space_energy = s_space_energy.get_default_codomain()
h_space = DomainTuple.make((h_space_spatial, h_space_energy))
......
......@@ -49,18 +49,20 @@ class GaussianEnergy(Energy):
def value(self):
if self._cov is None:
return .5 * self.residual.vdot(self.residual).real
return .5 * self.residual.vdot(self._cov.inverse(self.residual)).real
return .5 * self.residual.vdot(
self._cov.inverse_times(self.residual)).real
@property
@memo
def gradient(self):
if self._cov is None:
return self._inp.gradient.adjoint(self.residual)
return self._inp.gradient.adjoint(self._cov.inverse(self.residual))
return self._inp.jacobian.adjoint_times(self.residual)
return self._inp.jacobian.adjoint_times(
self._cov.inverse_times(self.residual))
@property
@memo
def curvature(self):
def metric(self):
if self._cov is None:
return SandwichOperator.make(self._inp.gradient, None)
return SandwichOperator.make(self._inp.gradient, self._cov.inverse)
return SandwichOperator.make(self._inp.jacobian, None)
return SandwichOperator.make(self._inp.jacobian, self._cov.inverse)
......@@ -28,7 +28,7 @@ class PointSources(Model):
@property
@memo
def gradient(self):
def jacobian(self):
u = self.position['points'].local_data
inner = norm.pdf(u)
outer_inv = invgamma.pdf(invgamma.ppf(norm.cdf(u),
......@@ -51,7 +51,8 @@ class PointSources(Model):
@staticmethod
def IG_prime(field, alpha, q):
inner = norm.pdf(field.local_data)
outer = invgamma.pdf(invgamma.ppf(norm.cdf(field.local_data), alpha, scale=q), alpha, scale=q)
outer = invgamma.pdf(invgamma.ppf(norm.cdf(field.local_data), alpha,
scale=q), alpha, scale=q)
# # FIXME
# outer = np.clip(outer, 1e-20, None)
outer = 1/outer
......
......@@ -40,11 +40,11 @@ class PoissonianEnergy(Energy):
self._value = lamb_val.sum() - d.vdot(log(lamb_val))
if isnan(self._value):
self._value = inf
self._gradient = self._lamb.gradient.adjoint_times(1 - d/lamb_val)
self._gradient = self._lamb.jacobian.adjoint_times(1 - d/lamb_val)
# metric = makeOp(d/lamb_val/lamb_val)
metric = makeOp(1./lamb_val)
self._curvature = SandwichOperator.make(self._lamb.gradient, metric)
self._metric = SandwichOperator.make(self._lamb.jacobian, metric)
def at(self, position):
return self.__class__(self._lamb.at(position), self._d)
......@@ -58,5 +58,5 @@ class PoissonianEnergy(Energy):
return self._gradient
@property
def curvature(self):
return self._curvature
def metric(self):
return self._metric
......@@ -47,7 +47,7 @@ class ConjugateGradient(Minimizer):
Parameters
----------
energy : Energy object at the starting point of the iteration.
Its curvature operator must be independent of position, otherwise
Its metric operator must be independent of position, otherwise
linear conjugate gradient minimization will fail.
preconditioner : Operator *optional*
This operator can be provided which transforms the variables of the
......@@ -73,7 +73,7 @@ class ConjugateGradient(Minimizer):
return energy, controller.CONVERGED
while True:
q = energy.curvature(d)
q = energy.metric(d)
ddotq = d.vdot(q).real
if ddotq == 0.:
logger.error("Error: ConjugateGradient: ddotq==0.")
......
......@@ -51,7 +51,7 @@ class DescentMinimizer(Minimizer):
Parameters
----------
energy : Energy
Energy object which provides value, gradient and curvature at a
Energy object which provides value, gradient and metric at a
specific position in parameter space.
Returns
......
......@@ -25,7 +25,7 @@ class Energy(NiftyMetaBase()):
""" Provides the functional used by minimization schemes.
The Energy object is an implementation of a scalar function including its
gradient and curvature at some position.
gradient and metric at some position.
Parameters
----------
......@@ -35,11 +35,11 @@ class Energy(NiftyMetaBase()):
Notes
-----
An instance of the Energy class is defined at a certain location. If one
is interested in the value, gradient or curvature of the abstract energy
is interested in the value, gradient or metric of the abstract energy
functional one has to 'jump' to the new position using the `at` method.
This method returns a new energy instance residing at the new position. By
this approach, intermediate results from computing e.g. the gradient can
safely be reused for e.g. the value or the curvature.
safely be reused for e.g. the value or the metric.
Memorizing the evaluations of some quantities (using the memo decorator)
minimizes the computational effort for multiple calls.
......@@ -75,7 +75,7 @@ class Energy(NiftyMetaBase()):
Field : selected location in parameter space.
The Field location in parameter space where value, gradient and
curvature are evaluated.
metric are evaluated.
"""
return self._position
......@@ -104,11 +104,11 @@ class Energy(NiftyMetaBase()):
return self.gradient.norm()
@property
def curvature(self):
def metric(self):
"""
LinearOperator : implicitly defined curvature.
LinearOperator : implicitly defined metric.
A positive semi-definite operator or function describing the
curvature of the potential at the given `position`.
metric of the potential at the given `position`.
"""
raise NotImplementedError
......@@ -132,7 +132,7 @@ class Energy(NiftyMetaBase()):
from .iteration_controller import IterationController
if not isinstance(controller, IterationController):
raise TypeError
return CurvatureInversionEnabler(self, controller, preconditioner)
return MetricInversionEnabler(self, controller, preconditioner)
def __mul__(self, factor):
from .energy_sum import EnergySum
......@@ -160,9 +160,9 @@ class Energy(NiftyMetaBase()):
return EnergySum.make([self], [-1.])
class CurvatureInversionEnabler(Energy):
class MetricInversionEnabler(Energy):
def __init__(self, ene, controller, preconditioner):
super(CurvatureInversionEnabler, self).__init__(ene.position)
super(MetricInversionEnabler, self).__init__(ene.position)
self._energy = ene
self._controller = controller
self._preconditioner = preconditioner
......@@ -170,7 +170,7 @@ class CurvatureInversionEnabler(Energy):
def at(self, position):
if self._position.isSubsetOf(position):
return self
return CurvatureInversionEnabler(
return MetricInversionEnabler(
self._energy.at(position), self._controller, self._preconditioner)
@property
......@@ -186,16 +186,16 @@ class CurvatureInversionEnabler(Energy):
return self._energy.gradient
@property
def curvature(self):
def metric(self):
from ..operators.linear_operator import LinearOperator
from ..operators.inversion_enabler import InversionEnabler
curv = self._energy.curvature
curv = self._energy.metric
if self._preconditioner is None:
precond = None
elif isinstance(self._preconditioner, LinearOperator):
precond = self._preconditioner
elif isinstance(self._preconditioner, Energy):
precond = self._preconditioner.at(self.position).curvature
precond = self._preconditioner.at(self.position).metric
return InversionEnabler(curv, self._controller, precond)
def longest_step(self, dir):
......
......@@ -67,6 +67,6 @@ class EnergySum(Energy):
@property
@memo
def curvature(self):
return my_lincomb(map(lambda v: v.curvature, self._energies),
def metric(self):
return my_lincomb(map(lambda v: v.metric, self._energies),
self._factors)
......@@ -21,7 +21,7 @@ from .energy import Energy
class QuadraticEnergy(Energy):
"""The Energy for a quadratic form.
The most important aspect of this energy is that its curvature must be
The most important aspect of this energy is that its metric must be
position-independent.
"""
......@@ -74,5 +74,5 @@ class QuadraticEnergy(Energy):
return self._grad
@property
def curvature(self):
def metric(self):
return self._A
......@@ -24,7 +24,7 @@ class RelaxedNewton(DescentMinimizer):
""" Calculates the descent direction according to a Newton scheme.
The descent direction is determined by weighting the gradient at the
current parameter position with the inverse local curvature.
current parameter position with the inverse local metric.
"""
def __init__(self, controller, line_searcher=None):
if line_searcher is None:
......@@ -34,4 +34,4 @@ class RelaxedNewton(DescentMinimizer):
line_searcher=line_searcher)
def get_descent_direction(self, energy):
return -energy.curvature.inverse_times(energy.gradient)
return -energy.metric.inverse_times(energy.gradient)
......@@ -56,7 +56,7 @@ class _MinHelper(object):
def hessp(self, x, p):
self._update(x)
res = self._energy.curvature(_toField(p, self._domain))
res = self._energy.metric(_toField(p, self._domain))
return _toFlatNdarray(res)
......
......@@ -42,7 +42,7 @@ class ScalarMul(Model):
self._factor = factor
self._value = self._factor * self._model.value
self._gradient = self._factor * self._model.gradient
self._jacobian = self._factor * self._model.jacobian
def at(self, position):
return self.__class__(self._factor, self._model.at(position))
......@@ -57,7 +57,7 @@ class Add(Model):
self._model2 = model2.at(position)
self._value = self._model1.value + self._model2.value
self._gradient = self._model1.gradient + self._model2.gradient
self._jacobian = self._model1.jacobian + self._model2.jacobian
@staticmethod
def make(model1, model2):
......@@ -87,8 +87,8 @@ class Mul(Model):
self._model2 = model2.at(position)
self._value = self._model1.value * self._model2.value
self._gradient = (makeOp(self._model1.value) * self._model2.gradient +
makeOp(self._model2.value) * self._model1.gradient)
self._jacobian = (makeOp(self._model1.value) * self._model2.jacobian +
makeOp(self._model2.value) * self._model1.jacobian)
@staticmethod
def make(model1, model2):
......
......@@ -32,7 +32,7 @@ class Constant(Model):
-----
Since there is no model-function associated:
- Position has no influence on value.
- There is no gradient.
- There is no Jacobian.
"""
# TODO Remove position
def __init__(self, position, constant):
......@@ -40,7 +40,7 @@ class Constant(Model):
self._constant = constant
self._value = self._constant
self._gradient = 0.
self._jacobian = 0.
def at(self, position):
return self.__class__(position, self._constant)
......@@ -34,7 +34,7 @@ class LinearModel(Model):
-------
Model with linear Operator applied:
- Model.value = LinOp (inp.value) [key-wise]
- Gradient = LinOp * inp.gradient
- Jacobian = LinOp * inp.jacobian
"""
from ..operators.linear_operator import LinearOperator
super(LinearModel, self).__init__(inp.position)
......@@ -49,7 +49,7 @@ class LinearModel(Model):
self._lin_op._key)
self._value = self._lin_op(self._inp.value)
self._gradient = self._lin_op*self._inp.gradient
self._jacobian = self._lin_op*self._inp.jacobian
def at(self, position):
return self.__class__(self._inp.at(position), self._lin_op)
......@@ -27,7 +27,7 @@ class LocalModel(Model):
"""
Computes nonlinearity(inp)
- LocalModel.value = nonlinearity(value) (pointwise)
- LocalModel.gradient = Outer Product of gradients
- LocalModel.jacobian = Outer Product of Jacobians
Parameters
----------
......@@ -40,9 +40,9 @@ class LocalModel(Model):
self._inp = inp
self._nonlinearity = nonlinearity
self._value = nonlinearity(self._inp.value)
d_inner = self._inp.gradient
d_inner = self._inp.jacobian
d_outer = makeOp(self._nonlinearity.derivative(self._inp.value))
self._gradient = d_outer * d_inner
self._jacobian = d_outer * d_inner
def at(self, position):
return self.__class__(self._inp.at(position), self._nonlinearity)
......
......@@ -28,7 +28,7 @@ class Model(NiftyMetaBase()):
The Model object is an implementation of a * which knows:
- position in parameterspace. (Field, MulitField)
- value according to its modelfunction A. A(position)
- gradient of the modelfunction at the current position.
- Jacobian of the model function at the current position.
Parameters
----------
......@@ -37,9 +37,9 @@ class Model(NiftyMetaBase()):
Notes
-----
An instance of the model class knows its position, value and gradient.
An instance of the model class knows its position, value and Jacobian.
One can 'jump' to a new position, with the help of the 'at' method, whereby
one automatically gets the value and gradient of the model. The 'at' method
one automatically gets the value and Jacobian of the model. The 'at' method
creates a new instance of the class.
"""
def __init__(self, position):
......@@ -65,7 +65,7 @@ class Model(NiftyMetaBase()):
"""
Field or MultiField: selected location in parameter space.
The location in parameter space where value and gradient are
The location in parameter space where value and Jacobian are
evaluated.
"""
return self._position
......@@ -80,11 +80,11 @@ class Model(NiftyMetaBase()):
return self._value
@property
def gradient(self):
def jacobian(self):
"""
LinearOperator : The derivative of the model at given `position`.
"""
return self._gradient
return self._jacobian
def __getitem__(self, key):
sel = SelectionOperator(self.value.domain, key)
......
......@@ -34,7 +34,8 @@ class MultiModel(Model):
if not isinstance(val.domain, DomainTuple):
raise TypeError
self._value = MultiField({key: val})
self._gradient = MultiAdaptor(self.value.domain) * self._model.gradient
self._jacobian = (MultiAdaptor(self.value.domain) *
self._model.jacobian)
def at(self, position):
return self.__class__(self._model.at(position), self._key)
......@@ -32,7 +32,7 @@ class Variable(Model):
super(Variable, self).__init__(position)
self._value = position
self._gradient = ScalingOperator(1., position.domain)
self._jacobian = ScalingOperator(1., position.domain)
def at(self, position):
return self.__class__(position)
......@@ -94,9 +94,8 @@ class ChainOperator(LinearOperator):
@staticmethod
def make(ops):
"""Build a ChainOperator (or something simpler if possible),
"""Build a ChainOperator (or something simpler if possible),
a sequence of concatenated LinearOperators.
Parameters
----------
......
......@@ -277,7 +277,8 @@ class LinearOperator(NiftyMetaBase()):
if not self._validMode[mode]:
raise NotImplementedError("invalid operator mode specified")
if mode & self.capability == 0:
raise NotImplementedError("requested operator mode is not supported")
raise NotImplementedError(
"requested operator mode is not supported")
def _check_input(self, x, mode):
self._check_mode(mode)
......
......@@ -32,9 +32,9 @@ class SamplingEnabler(EndomorphicOperator):
Parameters
----------
likelihood : :class:`EndomorphicOperator`
Curvature of the likelihood
Metric of the likelihood
prior : :class:`EndomorphicOperator`
Inverse curvature of the prior
Inverse metric of the prior
iteration_controller : :class:`IterationController`
The iteration controller to use for the iterative numerical inversion
done by a :class:`ConjugateGradient` object.
......
......@@ -59,7 +59,7 @@ class Energy_Tests(unittest.TestCase):
S = ift.create_power_operator(hspace, power_spectrum=_flat_PS)
energy = ift.WienerFilterEnergy(
position=s0, d=d, R=R, N=N, S=S, iteration_controller=IC)