Commit acdf2265 authored by Reimar H Leike's avatar Reimar H Leike
Browse files

Merge branch 'NIFTy_5' into fix_mpi_kl

parents b43228de 046d074c
Pipeline #63008 passed with stages
in 7 minutes and 7 seconds
......@@ -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.
**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.
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.
......@@ -32,7 +32,7 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full",
"redistribute", "default_distaxis", "is_numpy", "absmax", "norm",
"lock", "locked", "uniform_full", "transpose", "to_global_data_rw",
"ensure_not_distributed", "ensure_default_distributed",
"tanh", "conjugate", "sin", "cos", "tan",
"tanh", "conjugate", "sin", "cos", "tan", "log10",
"sinh", "cosh", "sinc", "absolute", "sign", "clip"]
_comm = MPI.COMM_WORLD
......@@ -297,7 +297,7 @@ def _math_helper(x, function, out):
_current_module = sys.modules[__name__]
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 func2(x, out=None):
return _math_helper(x, f, out)
......
......@@ -18,9 +18,11 @@
# Data object module that uses simple numpy ndarrays.
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 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
......@@ -34,7 +36,7 @@ __all__ = ["ntask", "rank", "master", "local_shape", "data_object", "full",
"lock", "locked", "uniform_full", "to_global_data_rw",
"ensure_not_distributed", "ensure_default_distributed",
"clip", "sin", "cos", "tan", "sinh",
"cosh", "absolute", "sign", "sinc"]
"cosh", "absolute", "sign", "sinc", "log10"]
ntask = 1
rank = 0
......
......@@ -19,6 +19,8 @@ from functools import reduce
from . import utilities
from .domains.domain import Domain
import numpy as np
class DomainTuple(object):
"""Ordered sequence of Domain objects.
......@@ -125,6 +127,58 @@ class DomainTuple(object):
"""
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
def axes(self):
"""tuple of tuple of int : axis indices of the underlying domains"""
......
......@@ -246,42 +246,23 @@ class Field(object):
If the requested sub-domain has a uniform volume element, it is
returned. Otherwise, `None` is returned.
"""
if np.isscalar(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
return self._domain.scalar_weight(spaces)
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
----------
spaces : int, tuple of int or None
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
-------
float
the total volume of the requested sub-domain.
the total volume of the requested (sub-)domain.
"""
if np.isscalar(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
return self._domain.total_volume(spaces)
def weight(self, power=1, spaces=None):
"""Weights the pixels of `self` with their invidual pixel volumes.
......@@ -682,7 +663,7 @@ for op in ["__iadd__", "__isub__", "__imul__", "__idiv__",
return func2
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",
"absolute", "sinc", "sign"]:
def func(f):
......
......@@ -330,6 +330,11 @@ class Linearization(object):
tmp = self._val.log()
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):
tmp = self._val.sinh()
tmp2 = self._val.cosh()
......
......@@ -73,7 +73,6 @@ class KLMetric(EndomorphicOperator):
return self._KL.metric_sample(from_inverse, dtype)
class MetricGaussianKL_MPI(Energy):
"""Provides the sampled Kullback-Leibler divergence between a distribution
and a Metric Gaussian.
......@@ -116,7 +115,7 @@ class MetricGaussianKL_MPI(Energy):
_samples : None
Only a parameter for internal uses. Typically not to be set by users.
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
same every time this class is initialized.
......@@ -132,7 +131,6 @@ class MetricGaussianKL_MPI(Energy):
Torsten A. Enßlin, `<https://arxiv.org/abs/1901.11033>`_
"""
def __init__(self, mean, hamiltonian, n_samples, constants=[],
point_estimates=[], mirror_samples=False,
napprox=0, _samples=None, seed_offset=0):
......
......@@ -159,7 +159,7 @@ class Operator(metaclass=NiftyMeta):
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 func2(self):
fa = _FunctionApplier(self.target, f)
......
......@@ -171,14 +171,15 @@ def _find_closest(A, target):
return idx
def _makeplot(name):
def _makeplot(name, block=True):
import matplotlib.pyplot as plt
if dobj.rank != 0:
plt.close()
return
if name is None:
plt.show()
plt.close()
plt.show(block=block)
if block:
plt.close()
return
extension = os.path.splitext(name)[1]
if extension in (".pdf", ".png", ".svg"):
......@@ -511,6 +512,9 @@ class Plot(object):
If left empty, the plot will be shown on the screen,
otherwise it will be written to a file with the given name.
Supported extensions: .png and .pdf. Default: None.
block: bool
Override the blocking behavior of the non-interactive plotting
mode. The plot will not be closed in this case but is left open!
"""
import matplotlib.pyplot as plt
nplot = len(self._plots)
......@@ -537,4 +541,4 @@ class Plot(object):
ax = fig.add_subplot(ny, nx, i+1)
_plot(self._plots[i], ax, **self._kwargs[i])
fig.tight_layout()
_makeplot(kwargs.pop("name", None))
_makeplot(kwargs.pop("name", None), block=kwargs.pop("block", True))
......@@ -35,7 +35,7 @@ __all__ = ['PS_field', 'power_analyze', 'create_power_operator',
'create_harmonic_smoothing_operator', 'from_random',
'full', 'from_global_data', 'from_local_data',
'makeDomain', 'sqrt', 'exp', 'log', 'tanh', 'sigmoid',
'sin', 'cos', 'tan', 'sinh', 'cosh',
'sin', 'cos', 'tan', 'sinh', 'cosh', 'log10',
'absolute', 'one_over', 'clip', 'sinc',
'conjugate', 'get_signal_variance', 'makeOp', 'domain_union',
'get_default_codomain', 'single_plot']
......@@ -389,7 +389,7 @@ def domain_union(domains):
_current_module = sys.modules[__name__]
for f in ["sqrt", "exp", "log", "tanh", "sigmoid",
for f in ["sqrt", "exp", "log", "log10", "tanh", "sigmoid",
"conjugate", 'sin', 'cos', 'tan', 'sinh', 'cosh',
'absolute', 'one_over', 'sinc']:
def func(f):
......
......@@ -219,17 +219,23 @@ def test_weight():
f = ift.Field.full(s1, 10.)
f2 = f.weight(1)
assert_equal(f.weight(1).local_data, f2.local_data)
assert_equal(f.domain.total_volume(), 1)
assert_equal(f.domain.total_volume(0), 1)
assert_equal(f.domain.total_volume((0,)), 1)
assert_equal(f.total_volume(), 1)
assert_equal(f.total_volume(0), 1)
assert_equal(f.total_volume((0,)), 1)
assert_equal(f.domain.scalar_weight(), 0.1)
assert_equal(f.domain.scalar_weight(0), 0.1)
assert_equal(f.domain.scalar_weight((0,)), 0.1)
assert_equal(f.scalar_weight(), 0.1)
assert_equal(f.scalar_weight(0), 0.1)
assert_equal(f.scalar_weight((0,)), 0.1)
s1 = ift.GLSpace(10)
f = ift.Field.full(s1, 10.)
assert_equal(f.scalar_weight(), None)
assert_equal(f.scalar_weight(0), None)
assert_equal(f.scalar_weight((0,)), None)
assert_equal(f.domain.scalar_weight(), None)
assert_equal(f.domain.scalar_weight(0), None)
assert_equal(f.domain.scalar_weight((0,)), None)
@pmp('dom', [ift.RGSpace(10), ift.GLSpace(10)])
......@@ -329,7 +335,7 @@ def test_emptydomain():
@pmp('dom', [ift.RGSpace((8,), harmonic=True), ()])
@pmp('func', [
"exp", "log", "sin", "cos", "tan", "sinh", "cosh", "sinc", "absolute",
"sign"
"sign", "log10"
])
def test_funcs(num, dom, func):
num = 5
......
......@@ -54,7 +54,7 @@ def test_special_gradients():
@pmp('f', [
'log', 'exp', 'sqrt', 'sin', 'cos', 'tan', 'sinc', 'sinh', 'cosh', 'tanh',
'absolute', 'one_over', 'sigmoid'
'absolute', 'one_over', 'sigmoid', 'log10'
])
def test_actual_gradients(f):
dom = ift.UnstructuredDomain((1,))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment