Commit e01ad3a5 authored by Philipp Arras's avatar Philipp Arras
Browse files

Merge branch 'NIFTy_5' into normalized_amplitudes_pp

parents 61a550d1 ec4c6f1a
Pipeline #63204 passed with stages
in 5 minutes and 51 seconds
...@@ -211,3 +211,18 @@ Consequently, the inverse covariance operator will automatically have lower indi ...@@ -211,3 +211,18 @@ Consequently, the inverse covariance operator will automatically have lower indi
ensuring that the whole log-likelihood expression does not scale with resolution. ensuring that the whole log-likelihood expression does not scale with resolution.
**This upper-lower index convention is not coded into NIFTy**, in order to not reduce user freedom. **This upper-lower index convention is not coded into NIFTy**, in order to not reduce user freedom.
One should however have this in mind when constructing log-likelihoods in order to ensure resolution independence. One should however have this in mind when constructing log-likelihoods in order to ensure resolution independence.
Harmonic Transform Convention
.............................
In NIFTy the convention for the harmonic transformations is set by requiring the zero mode of the transformed field to be the integral over the original field.
This choice is convenient for the Fourier transformation and used throughout the literature.
Note that for the spherical harmonics this convention is only rarely used and might result in unexpected factors in the transformed field.
To be specific, for the spherical harmonics transformation this means that a monopole of unit amplitude in position-space which is transformed to the spherical harmonics domain and back to the original position space via the adjoint transformation will have a non-unit amplitude afterwards.
The factor between the input field and the twice transformed field is :math:`1/4\pi`.
In comparison to the convention used in HEALPix, this corresponds to dividing the output of the HEALPix transformed field by :math:`\sqrt{4\pi}` in each transformation.
Depending on the use-case, additional volume factors must be accounted for.
This is especially true if one wants to define the inverse transformation.
Note that the convention for the harmonic transformations used in NIFTy results in non-unitary operators.
...@@ -19,6 +19,7 @@ from .multi_field import MultiField ...@@ -19,6 +19,7 @@ from .multi_field import MultiField
from .operators.operator import Operator from .operators.operator import Operator
from .operators.adder import Adder from .operators.adder import Adder
from .operators.log1p import Log1p
from .operators.diagonal_operator import DiagonalOperator from .operators.diagonal_operator import DiagonalOperator
from .operators.distributors import DOFDistributor, PowerDistributor from .operators.distributors import DOFDistributor, PowerDistributor
from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter
...@@ -46,7 +47,7 @@ from .operators.value_inserter import ValueInserter ...@@ -46,7 +47,7 @@ from .operators.value_inserter import ValueInserter
from .operators.energy_operators import ( from .operators.energy_operators import (
EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood, EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
BernoulliEnergy, StandardHamiltonian, AveragedEnergy, QuadraticFormOperator, BernoulliEnergy, StandardHamiltonian, AveragedEnergy, QuadraticFormOperator,
Squared2NormOperator) Squared2NormOperator, StudentTEnergy)
from .operators.convolution_operators import FuncConvolutionOperator from .operators.convolution_operators import FuncConvolutionOperator
from .probing import probe_with_posterior_samples, probe_diagonal, \ from .probing import probe_with_posterior_samples, probe_diagonal, \
......
...@@ -32,7 +32,7 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full", ...@@ -32,7 +32,7 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full",
"redistribute", "default_distaxis", "is_numpy", "absmax", "norm", "redistribute", "default_distaxis", "is_numpy", "absmax", "norm",
"lock", "locked", "uniform_full", "transpose", "to_global_data_rw", "lock", "locked", "uniform_full", "transpose", "to_global_data_rw",
"ensure_not_distributed", "ensure_default_distributed", "ensure_not_distributed", "ensure_default_distributed",
"tanh", "conjugate", "sin", "cos", "tan", "tanh", "conjugate", "sin", "cos", "tan", "log10",
"sinh", "cosh", "sinc", "absolute", "sign", "clip"] "sinh", "cosh", "sinc", "absolute", "sign", "clip"]
_comm = MPI.COMM_WORLD _comm = MPI.COMM_WORLD
...@@ -297,7 +297,7 @@ def _math_helper(x, function, out): ...@@ -297,7 +297,7 @@ def _math_helper(x, function, out):
_current_module = sys.modules[__name__] _current_module = sys.modules[__name__]
for f in ["sqrt", "exp", "log", "tanh", "conjugate", "sin", "cos", "tan", for f in ["sqrt", "exp", "log", "tanh", "conjugate", "sin", "cos", "tan",
"sinh", "cosh", "sinc", "absolute", "sign"]: "sinh", "cosh", "sinc", "absolute", "sign", "log10"]:
def func(f): def func(f):
def func2(x, out=None): def func2(x, out=None):
return _math_helper(x, f, out) return _math_helper(x, f, out)
......
...@@ -18,9 +18,11 @@ ...@@ -18,9 +18,11 @@
# Data object module that uses simple numpy ndarrays. # Data object module that uses simple numpy ndarrays.
import numpy as np import numpy as np
from numpy import absolute, clip, cos, cosh, empty, empty_like, exp, full, log
from numpy import ndarray as data_object from numpy import ndarray as data_object
from numpy import ones, sign, sin, sinc, sinh, sqrt, tan, tanh, vdot, zeros from numpy import empty, empty_like, ones, zeros, full
from numpy import absolute, sign, clip, vdot
from numpy import sin, cos, sinh, cosh, tan, tanh
from numpy import exp, log, log10, sqrt, sinc
from .random import Random from .random import Random
...@@ -34,7 +36,7 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full", ...@@ -34,7 +36,7 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full",
"lock", "locked", "uniform_full", "to_global_data_rw", "lock", "locked", "uniform_full", "to_global_data_rw",
"ensure_not_distributed", "ensure_default_distributed", "ensure_not_distributed", "ensure_default_distributed",
"clip", "sin", "cos", "tan", "sinh", "clip", "sin", "cos", "tan", "sinh",
"cosh", "absolute", "sign", "sinc"] "cosh", "absolute", "sign", "sinc", "log10"]
ntask = 1 ntask = 1
rank = 0 rank = 0
......
...@@ -19,6 +19,8 @@ from functools import reduce ...@@ -19,6 +19,8 @@ from functools import reduce
from . import utilities from . import utilities
from .domains.domain import Domain from .domains.domain import Domain
import numpy as np
class DomainTuple(object): class DomainTuple(object):
"""Ordered sequence of Domain objects. """Ordered sequence of Domain objects.
...@@ -125,6 +127,58 @@ class DomainTuple(object): ...@@ -125,6 +127,58 @@ class DomainTuple(object):
""" """
return self._size return self._size
def scalar_weight(self, spaces=None):
"""Returns the uniform volume element for a sub-domain of `self`.
Parameters
----------
spaces : int, tuple of int or None
Indices of the sub-domains 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
returned. Otherwise, `None` is returned.
"""
if np.isscalar(spaces):
return self._dom[spaces].scalar_dvol
if spaces is None:
spaces = range(len(self._dom))
res = 1.
for i in spaces:
tmp = self._dom[i].scalar_dvol
if tmp is None:
return None
res *= tmp
return res
def total_volume(self, spaces=None):
"""Returns the total volume of `self` or of a subspace of it.
Parameters
----------
spaces : int, tuple of int or None
Indices of the sub-domains of the domain to be considered.
If `None`, the total volume of the whole domain is returned.
Returns
-------
float
the total volume of the requested (sub-)domain.
"""
if np.isscalar(spaces):
return self._dom[spaces].total_volume
if spaces is None:
spaces = range(len(self._dom))
res = 1.
for i in spaces:
res *= self._dom[i].total_volume
return res
@property @property
def axes(self): def axes(self):
"""tuple of tuple of int : axis indices of the underlying domains""" """tuple of tuple of int : axis indices of the underlying domains"""
......
...@@ -16,11 +16,13 @@ ...@@ -16,11 +16,13 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np import numpy as np
from numpy.testing import assert_
from .domain_tuple import DomainTuple from .domain_tuple import DomainTuple
from .field import Field from .field import Field
from .linearization import Linearization from .linearization import Linearization
from .multi_domain import MultiDomain from .multi_domain import MultiDomain
from .multi_field import MultiField
from .operators.linear_operator import LinearOperator from .operators.linear_operator import LinearOperator
from .sugar import from_random from .sugar import from_random
...@@ -81,6 +83,38 @@ def _check_linearity(op, domain_dtype, atol, rtol): ...@@ -81,6 +83,38 @@ def _check_linearity(op, domain_dtype, atol, rtol):
_assert_allclose(val1, val2, atol=atol, rtol=rtol) _assert_allclose(val1, val2, atol=atol, rtol=rtol)
def _actual_domain_check(op, domain_dtype=None, inp=None):
needed_cap = op.TIMES
if (op.capability & needed_cap) != needed_cap:
return
if domain_dtype is not None:
inp = from_random("normal", op.domain, dtype=domain_dtype)
elif inp is None:
raise ValueError('Need to specify either dtype or inp')
assert_(inp.domain is op.domain)
assert_(op(inp).domain is op.target)
def _actual_domain_check_nonlinear(op, loc):
assert isinstance(loc, (Field, MultiField))
assert_(loc.domain is op.domain)
lin = Linearization.make_var(loc, False)
reslin = op(lin)
assert_(lin.domain is op.domain)
assert_(lin.target is op.domain)
assert_(lin.val.domain is lin.domain)
assert_(reslin.domain is op.domain)
assert_(reslin.target is op.target)
assert_(reslin.val.domain is reslin.target)
assert_(reslin.target is op.target)
assert_(reslin.jac.domain is reslin.domain)
assert_(reslin.jac.target is reslin.target)
_actual_domain_check(reslin.jac, inp=loc)
_actual_domain_check(reslin.jac.adjoint, inp=reslin.jac(loc))
def _domain_check(op): def _domain_check(op):
for dd in [op.domain, op.target]: for dd in [op.domain, op.target]:
if not isinstance(dd, (DomainTuple, MultiDomain)): if not isinstance(dd, (DomainTuple, MultiDomain)):
...@@ -123,6 +157,10 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64, ...@@ -123,6 +157,10 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
if not isinstance(op, LinearOperator): if not isinstance(op, LinearOperator):
raise TypeError('This test tests only linear operators.') raise TypeError('This test tests only linear operators.')
_domain_check(op) _domain_check(op)
_actual_domain_check(op, domain_dtype)
_actual_domain_check(op.adjoint, target_dtype)
_actual_domain_check(op.inverse, target_dtype)
_actual_domain_check(op.adjoint.inverse, domain_dtype)
_check_linearity(op, domain_dtype, atol, rtol) _check_linearity(op, domain_dtype, atol, rtol)
_check_linearity(op.adjoint, target_dtype, atol, rtol) _check_linearity(op.adjoint, target_dtype, atol, rtol)
_check_linearity(op.inverse, target_dtype, atol, rtol) _check_linearity(op.inverse, target_dtype, atol, rtol)
...@@ -180,6 +218,7 @@ def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100): ...@@ -180,6 +218,7 @@ def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100):
Tolerance for the check. Tolerance for the check.
""" """
_domain_check(op) _domain_check(op)
_actual_domain_check_nonlinear(op, loc)
for _ in range(ntries): for _ in range(ntries):
lin = op(Linearization.make_var(loc)) lin = op(Linearization.make_var(loc))
loc2, lin2 = _get_acceptable_location(op, loc, lin) loc2, lin2 = _get_acceptable_location(op, loc, lin)
......
...@@ -246,42 +246,23 @@ class Field(object): ...@@ -246,42 +246,23 @@ class Field(object):
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. returned. Otherwise, `None` is returned.
""" """
if np.isscalar(spaces): return self._domain.scalar_weight(spaces)
return self._domain[spaces].scalar_dvol
if spaces is None:
spaces = range(len(self._domain))
res = 1.
for i in spaces:
tmp = self._domain[i].scalar_dvol
if tmp is None:
return None
res *= tmp
return res
def total_volume(self, spaces=None): def total_volume(self, spaces=None):
"""Returns the total volume of a sub-domain of `self`. """Returns the total volume of the field's domain or of a subspace of it.
Parameters Parameters
---------- ----------
spaces : int, tuple of int or None 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. If `None`, the total volume of the whole domain is returned.
Returns Returns
------- -------
float float
the total volume of the requested sub-domain. the total volume of the requested (sub-)domain.
""" """
if np.isscalar(spaces): return self._domain.total_volume(spaces)
return self._domain[spaces].total_volume
if spaces is None:
spaces = range(len(self._domain))
res = 1.
for i in spaces:
res *= self._domain[i].total_volume
return res
def weight(self, power=1, spaces=None): def weight(self, power=1, spaces=None):
"""Weights the pixels of `self` with their invidual pixel volumes. """Weights the pixels of `self` with their invidual pixel volumes.
...@@ -682,7 +663,7 @@ for op in ["__iadd__", "__isub__", "__imul__", "__idiv__", ...@@ -682,7 +663,7 @@ for op in ["__iadd__", "__isub__", "__imul__", "__idiv__",
return func2 return func2
setattr(Field, op, func(op)) setattr(Field, op, func(op))
for f in ["sqrt", "exp", "log", "tanh", for f in ["sqrt", "exp", "log", "log10", "tanh",
"sin", "cos", "tan", "cosh", "sinh", "sin", "cos", "tan", "cosh", "sinh",
"absolute", "sinc", "sign"]: "absolute", "sinc", "sign"]:
def func(f): def func(f):
......
...@@ -63,7 +63,7 @@ class _LightConeDerivative(LinearOperator): ...@@ -63,7 +63,7 @@ class _LightConeDerivative(LinearOperator):
if mode == self.TIMES: if mode == self.TIMES:
res += self._derivatives[i]*x[i] res += self._derivatives[i]*x[i]
else: else:
res[i] = np.sum(self._derivatives[i]*x) res[i] = np.sum(self._derivatives[i]*x.real)
return Field.from_global_data(self._tgt(mode), res) return Field.from_global_data(self._tgt(mode), res)
......
...@@ -330,6 +330,11 @@ class Linearization(object): ...@@ -330,6 +330,11 @@ class Linearization(object):
tmp = self._val.log() tmp = self._val.log()
return self.new(tmp, makeOp(1./self._val)(self._jac)) return self.new(tmp, makeOp(1./self._val)(self._jac))
def log10(self):
tmp = self._val.log10()
tmp2 = 1. / (self._val * np.log(10))
return self.new(tmp, makeOp(tmp2)(self._jac))
def sinh(self): def sinh(self):
tmp = self._val.sinh() tmp = self._val.sinh()
tmp2 = self._val.cosh() tmp2 = self._val.cosh()
......
...@@ -73,7 +73,6 @@ class KLMetric(EndomorphicOperator): ...@@ -73,7 +73,6 @@ class KLMetric(EndomorphicOperator):
return self._KL.metric_sample(from_inverse, dtype) return self._KL.metric_sample(from_inverse, dtype)
class MetricGaussianKL_MPI(Energy): class MetricGaussianKL_MPI(Energy):
"""Provides the sampled Kullback-Leibler divergence between a distribution """Provides the sampled Kullback-Leibler divergence between a distribution
and a Metric Gaussian. and a Metric Gaussian.
...@@ -116,7 +115,7 @@ class MetricGaussianKL_MPI(Energy): ...@@ -116,7 +115,7 @@ class MetricGaussianKL_MPI(Energy):
_samples : None _samples : None
Only a parameter for internal uses. Typically not to be set by users. Only a parameter for internal uses. Typically not to be set by users.
seed_offset : int seed_offset : int
A parameter with which one can controll from which seed the samples A parameter with which one can controll from which seed the samples
are drawn. Per default, the seed is different for MPI tasks, but the are drawn. Per default, the seed is different for MPI tasks, but the
same every time this class is initialized. same every time this class is initialized.
...@@ -132,7 +131,6 @@ class MetricGaussianKL_MPI(Energy): ...@@ -132,7 +131,6 @@ class MetricGaussianKL_MPI(Energy):
Torsten A. Enßlin, `<https://arxiv.org/abs/1901.11033>`_ Torsten A. Enßlin, `<https://arxiv.org/abs/1901.11033>`_
""" """
def __init__(self, mean, hamiltonian, n_samples, constants=[], def __init__(self, mean, hamiltonian, n_samples, constants=[],
point_estimates=[], mirror_samples=False, point_estimates=[], mirror_samples=False,
napprox=0, _samples=None, seed_offset=0): napprox=0, _samples=None, seed_offset=0):
...@@ -161,6 +159,7 @@ class MetricGaussianKL_MPI(Energy): ...@@ -161,6 +159,7 @@ class MetricGaussianKL_MPI(Energy):
if napprox > 1: if napprox > 1:
met._approximation = makeOp(approximation2endo(met, napprox)) met._approximation = makeOp(approximation2endo(met, napprox))
_samples = [] _samples = []
rand_state = np.random.get_state()
for i in range(lo, hi): for i in range(lo, hi):
if mirror_samples: if mirror_samples:
np.random.seed(i//2+seed_offset) np.random.seed(i//2+seed_offset)
...@@ -171,8 +170,9 @@ class MetricGaussianKL_MPI(Energy): ...@@ -171,8 +170,9 @@ class MetricGaussianKL_MPI(Energy):
_samples.append(((i % 2)*2-1) * _samples.append(((i % 2)*2-1) *
met.draw_sample(from_inverse=True)) met.draw_sample(from_inverse=True))
else: else:
np.random.seed(i) np.random.seed(i+seed_offset)
_samples.append(met.draw_sample(from_inverse=True)) _samples.append(met.draw_sample(from_inverse=True))
np.random.set_state(rand_state)
_samples = tuple(_samples) _samples = tuple(_samples)
if mirror_samples: if mirror_samples:
n_samples *= 2 n_samples *= 2
...@@ -242,8 +242,11 @@ class MetricGaussianKL_MPI(Energy): ...@@ -242,8 +242,11 @@ class MetricGaussianKL_MPI(Energy):
raise NotImplementedError() raise NotImplementedError()
lin = self._lin.with_want_metric() lin = self._lin.with_want_metric()
samp = full(self._hamiltonian.domain, 0.) samp = full(self._hamiltonian.domain, 0.)
rand_state = np.random.get_state()
np.random.seed(rank+np.random.randint(99999))
for v in self._samples: for v in self._samples:
samp = samp + self._hamiltonian(lin+v).metric.draw_sample(from_inverse=False, dtype=dtype) samp = samp + self._hamiltonian(lin+v).metric.draw_sample(from_inverse=False, dtype=dtype)
np.random.set_state(rand_state)
return allreduce_sum_field(samp) return allreduce_sum_field(samp)
def metric_sample(self, from_inverse=False, dtype=np.float64): def metric_sample(self, from_inverse=False, dtype=np.float64):
......
...@@ -27,6 +27,7 @@ from .linear_operator import LinearOperator ...@@ -27,6 +27,7 @@ from .linear_operator import LinearOperator
from .operator import Operator from .operator import Operator
from .sampling_enabler import SamplingEnabler from .sampling_enabler import SamplingEnabler
from .sandwich_operator import SandwichOperator from .sandwich_operator import SandwichOperator
from .scaling_operator import ScalingOperator
from .simple_linear_operators import VdotOperator from .simple_linear_operators import VdotOperator
...@@ -248,6 +249,43 @@ class InverseGammaLikelihood(EnergyOperator): ...@@ -248,6 +249,43 @@ class InverseGammaLikelihood(EnergyOperator):
return res.add_metric(metric) return res.add_metric(metric)
class StudentTEnergy(EnergyOperator):
"""Computes likelihood energy of expected event frequency constrained by
event data.
.. math ::
E(f) = -\\log \\text{Bernoulli}(d|f)
= -d^\\dagger \\log f - (1-d)^\\dagger \\log(1-f),
where f is a field defined on `d.domain` with the expected
frequencies of events.
Parameters
----------
d : Field
Data field with events (1) or non-events (0).
theta : Scalar
Degree of freedom parameter for the student t distribution
"""
def __init__(self, domain, theta):
self._domain = DomainTuple.make(domain)
self._theta = theta
from .log1p import Log1p
self._l1p = Log1p(domain)
def apply(self, x):
self._check_input(x)
v = ((self._theta+1)/2)*self._l1p(x**2/self._theta).sum()
if not isinstance(x, Linearization):
return Field.scalar(v)
if not x.want_metric:
return v
met = ScalingOperator((self._theta+1)/(self._theta+3), self.domain)
met = SandwichOperator.make(x.jac, met)
return v.add_metric(met)
class BernoulliEnergy(EnergyOperator): class BernoulliEnergy(EnergyOperator):
"""Computes likelihood energy of expected event frequency constrained by """Computes likelihood energy of expected event frequency constrained by
event data. event data.
......
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..field import Field
from ..multi_field import MultiField
from .operator import Operator
from .diagonal_operator import DiagonalOperator
from ..linearization import Linearization
from ..sugar import from_local_data
from numpy import log1p
class Log1p(Operator):
"""computes x -> log(1+x)
"""
def __init__(self, dom):
self._domain = dom
self._target = dom
def apply(self, x):
lin = isinstance(x, Linearization)
xval = x.val if lin else x
xlval = xval.local_data
res = from_local_data(xval.domain, log1p(xlval))
if not lin:
return res
jac = DiagonalOperator(1/(1+xval))
return x.new(res, jac@x.jac)
...@@ -159,7 +159,7 @@ class Operator(metaclass=NiftyMeta): ...@@ -159,7 +159,7 @@ class Operator(metaclass=NiftyMeta):
for f in ["sqrt", "exp", "log", "tanh", "sigmoid", 'sin', 'cos', 'tan', for f in ["sqrt", "exp", "log", "tanh", "sigmoid", 'sin', 'cos', 'tan',
'sinh', 'cosh', 'absolute', 'sinc', 'one_over']: 'sinh', 'cosh', 'absolute', 'sinc', 'one_over', 'log10']:
def func(f): def func(f):
def func2(self): def func2(self):
fa = _FunctionApplier(self.target, f) fa = _FunctionApplier(self.target, f)
......
...@@ -187,9 +187,7 @@ class _SlowFieldAdapter(LinearOperator): ...@@ -187,9 +187,7 @@ class _SlowFieldAdapter(LinearOperator):
self._check_input(x, mode) self._check_input(x, mode)
if isinstance(x, MultiField): if isinstance(x, MultiField):