Commit 39a0aa47 authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'cosm' into 'NIFTy_5'

Various stuff by Philipp

See merge request ift/nifty-dev!194
parents a4b34dce e6f7a95d
FROM debian:testing-slim FROM debian:testing-slim
RUN apt-get update && apt-get install -y \ RUN apt-get update && apt-get install -y \
# Needed for gitlab tests # Needed for setup
git \ git python3-pip \
# Packages needed for NIFTy # Packages needed for NIFTy
libfftw3-dev \ python3-scipy \
python3 python3-pip python3-dev python3-scipy \
# Documentation build dependencies # Documentation build dependencies
python3-sphinx python3-sphinx-rtd-theme \ python3-sphinx-rtd-theme \
# Testing dependencies # Testing dependencies
python3-coverage python3-pytest python3-pytest-cov \ python3-pytest-cov jupyter \
# Optional NIFTy dependencies # Optional NIFTy dependencies
openmpi-bin libopenmpi-dev python3-mpi4py \ libfftw3-dev python3-mpi4py python3-matplotlib \
# Packages needed for NIFTy # more optional NIFTy dependencies
&& pip3 install pyfftw \ && pip3 install pyfftw \
# Optional NIFTy dependencies
&& pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \ && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \
# Testing dependencies && pip3 install jupyter \
&& 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 \
&& rm -rf /var/lib/apt/lists/* && rm -rf /var/lib/apt/lists/*
# Set matplotlib backend # Set matplotlib backend
......
...@@ -29,4 +29,4 @@ add_module_names = False ...@@ -29,4 +29,4 @@ add_module_names = False
html_theme = "sphinx_rtd_theme" html_theme = "sphinx_rtd_theme"
html_logo = 'nifty_logo_black.png' html_logo = 'nifty_logo_black.png'
exclude_patterns = ['mod/modules.rst'] exclude_patterns = ['mod/modules.rst', 'mod/*version.rst']
...@@ -219,14 +219,14 @@ Here, the uncertainty of the field and the power spectrum of its generating proc ...@@ -219,14 +219,14 @@ Here, the uncertainty of the field and the power spectrum of its generating proc
| **Output of tomography demo getting_started_3.py** | | **Output of tomography demo getting_started_3.py** |
+----------------------------------------------------+ +----------------------------------------------------+
| .. image:: images/getting_started_3_setup.png | | .. image:: images/getting_started_3_setup.png |
| :width: 50 % | | |
+----------------------------------------------------+ +----------------------------------------------------+
| Non-Gaussian signal field, | | Non-Gaussian signal field, |
| data backprojected into the image domain, power | | data backprojected into the image domain, power |
| spectrum of underlying Gausssian process. | | spectrum of underlying Gausssian process. |
+----------------------------------------------------+ +----------------------------------------------------+
| .. image:: images/getting_started_3_results.png | | .. image:: images/getting_started_3_results.png |
| :width: 50 % | | |
+----------------------------------------------------+ +----------------------------------------------------+
| Posterior mean field signal | | Posterior mean field signal |
| reconstruction, its uncertainty, and the power | | reconstruction, its uncertainty, and the power |
...@@ -238,7 +238,7 @@ Here, the uncertainty of the field and the power spectrum of its generating proc ...@@ -238,7 +238,7 @@ Here, the uncertainty of the field and the power spectrum of its generating proc
Maximum a Posteriori Maximum a Posteriori
-------------------- --------------------
One popular field estimation method is Maximim a Posteriori (MAP). One popular field estimation method is Maximum a Posteriori (MAP).
It only requires to minimize the information Hamiltonian, e.g by a gradient descent method that stops when It only requires to minimize the information Hamiltonian, e.g by a gradient descent method that stops when
......
...@@ -19,6 +19,7 @@ import numpy as np ...@@ -19,6 +19,7 @@ import numpy as np
from .field import Field from .field import Field
from .linearization import Linearization from .linearization import Linearization
from .operators.linear_operator import LinearOperator
from .sugar import from_random from .sugar import from_random
__all__ = ["consistency_check", "check_value_gradient_consistency", __all__ = ["consistency_check", "check_value_gradient_consistency",
...@@ -62,8 +63,40 @@ def _full_implementation(op, domain_dtype, target_dtype, atol, rtol): ...@@ -62,8 +63,40 @@ def _full_implementation(op, domain_dtype, target_dtype, atol, rtol):
_inverse_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, def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
atol=0, rtol=1e-7): 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, domain_dtype, target_dtype, atol, rtol)
_full_implementation(op.adjoint, target_dtype, domain_dtype, atol, rtol) _full_implementation(op.adjoint, target_dtype, domain_dtype, atol, rtol)
_full_implementation(op.inverse, 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): ...@@ -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): def check_value_gradient_consistency(op, loc, tol=1e-8, ntries=100):
"""FIXME"""
_check_consistency(op, loc, tol, ntries, False) _check_consistency(op, loc, tol, ntries, False)
def check_value_gradient_metric_consistency(op, loc, tol=1e-8, ntries=100): def check_value_gradient_metric_consistency(op, loc, tol=1e-8, ntries=100):
"""FIXME"""
_check_consistency(op, loc, tol, ntries, True) _check_consistency(op, loc, tol, ntries, True)
...@@ -31,13 +31,13 @@ class Field(object): ...@@ -31,13 +31,13 @@ class Field(object):
Parameters Parameters
---------- ----------
domain : DomainTuple domain : DomainTuple
the domain of the new Field The domain of the new Field.
val : data_object val : data_object
This object's global shape must match the domain shape This object's global shape must match the domain shape
After construction, the object will no longer be writeable! After construction, the object will no longer be writeable!
Notes Note
----- ----
If possible, do not invoke the constructor directly, but use one of the If possible, do not invoke the constructor directly, but use one of the
many convenience functions for instantiation! many convenience functions for instantiation!
""" """
...@@ -76,14 +76,14 @@ class Field(object): ...@@ -76,14 +76,14 @@ class Field(object):
Parameters Parameters
---------- ----------
domain : Domain, tuple of Domain, or DomainTuple domain : Domain, tuple of Domain, or DomainTuple
domain of the new Field Domain of the new Field.
val : float/complex/int scalar 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 Returns
------- -------
Field Field
the newly created field The newly created Field.
""" """
if not np.isscalar(val): if not np.isscalar(val):
raise TypeError("val must be a scalar") raise TypeError("val must be a scalar")
...@@ -99,7 +99,7 @@ class Field(object): ...@@ -99,7 +99,7 @@ class Field(object):
Parameters Parameters
---------- ----------
domain : DomainTuple, tuple of Domain, or Domain domain : DomainTuple, tuple of Domain, or Domain
the domain of the new Field The domain of the new Field.
arr : numpy.ndarray arr : numpy.ndarray
The data content to be used for the new Field. The data content to be used for the new Field.
Its shape must match the shape of `domain`. Its shape must match the shape of `domain`.
...@@ -132,8 +132,9 @@ class Field(object): ...@@ -132,8 +132,9 @@ class Field(object):
Returns Returns
------- -------
numpy.ndarray : array containing all field entries, which can be numpy.ndarray
modified. Its shape is identical to `self.shape`. Array containing all field entries, which can be modified. Its
shape is identical to `self.shape`.
""" """
return dobj.to_global_data_rw(self._val) return dobj.to_global_data_rw(self._val)
...@@ -171,9 +172,9 @@ class Field(object): ...@@ -171,9 +172,9 @@ class Field(object):
random_type : 'pm1', 'normal', or 'uniform' random_type : 'pm1', 'normal', or 'uniform'
The random distribution to use. The random distribution to use.
domain : DomainTuple domain : DomainTuple
The domain of the output random field The domain of the output random Field.
dtype : type dtype : type
The datatype of the output random field The datatype of the output random Field.
Returns Returns
------- -------
...@@ -187,10 +188,10 @@ class Field(object): ...@@ -187,10 +188,10 @@ class Field(object):
@property @property
def val(self): 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 This property is intended for low-level, internal use only. Do not use
from outside of NIFTy's core; there should be better alternatives. from outside of NIFTy's core; there should be better alternatives.
""" """
...@@ -236,13 +237,13 @@ class Field(object): ...@@ -236,13 +237,13 @@ class Field(object):
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 entire domain is used.
Returns Returns
------- -------
float or None 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. returned. Otherwise, `None` is returned.
""" """
if np.isscalar(spaces): if np.isscalar(spaces):
...@@ -264,7 +265,7 @@ class Field(object): ...@@ -264,7 +265,7 @@ class Field(object):
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 entire domain is used.
Returns Returns
...@@ -331,8 +332,9 @@ class Field(object): ...@@ -331,8 +332,9 @@ class Field(object):
x : Field x : Field
Returns 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): if not isinstance(x, Field):
raise TypeError("The multiplier must be an instance of " + raise TypeError("The multiplier must be an instance of " +
......
...@@ -16,6 +16,7 @@ ...@@ -16,6 +16,7 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce from functools import reduce
from operator import mul
from ..domain_tuple import DomainTuple from ..domain_tuple import DomainTuple
from ..operators.contraction_operator import ContractionOperator from ..operators.contraction_operator import ContractionOperator
...@@ -61,6 +62,10 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None): ...@@ -61,6 +62,10 @@ def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
power_distributor = PowerDistributor(h_space, p_space) power_distributor = PowerDistributor(h_space, p_space)
A = power_distributor(amplitude_operator) A = power_distributor(amplitude_operator)
vol = h_space.scalar_dvol**-0.5 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)) return ht(vol*A*ducktape(h_space, None, name))
...@@ -107,7 +112,8 @@ def MfCorrelatedField(target, amplitudes, name='xi'): ...@@ -107,7 +112,8 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
d = [dd0, dd1] d = [dd0, dd1]
a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)] a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
a = reduce(lambda x, y: x*y, a) a = reduce(mul, a)
A = pd @ a A = pd @ a
vol = reduce(lambda x, y: x*y, [sp.scalar_dvol**-0.5 for sp in hsp]) # For `vol` see comment in `CorrelatedField`
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
return ht(vol*A*ducktape(hsp, None, name)) return ht(vol*A*ducktape(hsp, None, name))
...@@ -130,4 +130,4 @@ class MetricGaussianKL(Energy): ...@@ -130,4 +130,4 @@ class MetricGaussianKL(Energy):
def __repr__(self): def __repr__(self):
return 'KL ({} samples):\n'.format(len( return 'KL ({} samples):\n'.format(len(
self._samples)) + utilities.indent(self._ham.__repr__()) self._samples)) + utilities.indent(self._hamiltonian.__repr__())
...@@ -23,31 +23,38 @@ from .linear_operator import LinearOperator ...@@ -23,31 +23,38 @@ from .linear_operator import LinearOperator
class DomainTupleFieldInserter(LinearOperator): class DomainTupleFieldInserter(LinearOperator):
"""Writes the content of a :class:`Field` into one slice of a :class:`DomainTuple`. """Writes the content of a :class:`Field` into one slice of a
:class:`DomainTuple`.
Parameters Parameters
---------- ----------
domain : Domain, tuple of Domain or DomainTuple target : Domain, tuple of Domain or DomainTuple
new_space : Domain, tuple of Domain or DomainTuple space : int
index : Integer The index of the sub-domain which is inserted.
Index at which new_space shall be added to domain. index : tuple
position : tuple Slice in new sub-domain in which the input field shall be written into.
Slice in new_space in which the input field shall be written into.
""" """
def __init__(self, domain, new_space, index, position):
self._domain = DomainTuple.make(domain) def __init__(self, target, space, pos):
tgt = list(self.domain) if not space <= len(target) or space < 0:
tgt.insert(index, new_space) raise ValueError
self._target = DomainTuple.make(tgt) self._target = DomainTuple.make(target)
dom = list(self.target)
dom.pop(space)
self._domain = DomainTuple.make(dom)
self._capability = self.TIMES | self.ADJOINT_TIMES self._capability = self.TIMES | self.ADJOINT_TIMES
fst_dims = sum(len(dd.shape) for dd in self.domain[:index])
new_space = target[space]
nshp = new_space.shape nshp = new_space.shape
if len(position) != len(nshp): fst_dims = sum(len(dd.shape) for dd in self.target[:space])
if len(pos) != len(nshp):
raise ValueError("shape mismatch between new_space and position") raise ValueError("shape mismatch between new_space and position")
for s, p in zip(nshp, position): for s, p in zip(nshp, pos):
if p < 0 or p >= s: if p < 0 or p >= s:
raise ValueError("bad position value") raise ValueError("bad position value")
self._slc = (slice(None),)*fst_dims + position
self._slc = (slice(None),)*fst_dims + pos
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
......
...@@ -22,10 +22,9 @@ import numpy as np ...@@ -22,10 +22,9 @@ import numpy as np
from . import dobj from . import dobj
from .domains.gl_space import GLSpace from .domains.gl_space import GLSpace
from .domains.hp_space import HPSpace from .domains.hp_space import HPSpace
from .domains.log_rg_space import LogRGSpace
from .domains.power_space import PowerSpace from .domains.power_space import PowerSpace
from .domains.rg_space import RGSpace from .domains.rg_space import RGSpace
from .domains.log_rg_space import LogRGSpace
from .domain_tuple import DomainTuple
from .field import Field from .field import Field
# relevant properties: # relevant properties:
...@@ -237,6 +236,8 @@ def _plot2D(f, ax, **kwargs): ...@@ -237,6 +236,8 @@ def _plot2D(f, ax, **kwargs):
foo = kwargs.pop("norm", None) foo = kwargs.pop("norm", None)
norm = {} if foo is None else {'norm': foo} norm = {} if foo is None else {'norm': foo}
aspect = kwargs.pop("aspect", None)
aspect = {} if foo is None else {'aspect': foo}
ax.set_title(kwargs.pop("title", "")) ax.set_title(kwargs.pop("title", ""))
ax.set_xlabel(kwargs.pop("xlabel", "")) ax.set_xlabel(kwargs.pop("xlabel", ""))
...@@ -250,7 +251,7 @@ def _plot2D(f, ax, **kwargs): ...@@ -250,7 +251,7 @@ def _plot2D(f, ax, **kwargs):
im = ax.imshow( im = ax.imshow(
f.to_global_data().T, extent=[0, nx*dx, 0, ny*dy], f.to_global_data().T, extent=[0, nx*dx, 0, ny*dy],
vmin=kwargs.get("zmin"), vmax=kwargs.get("zmax"), vmin=kwargs.get("zmin"), vmax=kwargs.get("zmax"),
cmap=cmap, origin="lower", **norm) cmap=cmap, origin="lower", **norm, **aspect)
plt.colorbar(im) plt.colorbar(im)
_limit_xy(**kwargs) _limit_xy(**kwargs)
return return
...@@ -313,8 +314,8 @@ class Plot(object): ...@@ -313,8 +314,8 @@ class Plot(object):
Notes Notes
----- -----
After doing one or more calls `plot()`, one also needs to call After doing one or more calls `add()`, one needs to call `output()` to
`plot_finish()` to output the result. show or save the plot.
Parameters Parameters
---------- ----------
...@@ -324,7 +325,7 @@ class Plot(object): ...@@ -324,7 +325,7 @@ class Plot(object):
If it is a list, all list members must be Fields defined over the If it is a list, all list members must be Fields defined over the
same one-dimensional `RGSpace` or `PowerSpace`. same one-dimensional `RGSpace` or `PowerSpace`.
title: string title: string
title of the plot. Title of the plot.
xlabel: string xlabel: string
Label for the x axis. Label for the x axis.
ylabel: string ylabel: string
...@@ -338,7 +339,7 @@ class Plot(object): ...@@ -338,7 +339,7 @@ class Plot(object):
label: string of list of strings label: string of list of strings
Annotation string. Annotation string.
alpha: float or list of floats alpha: float or list of floats
transparency value Transparency value.
""" """
self._plots.append(f) self._plots.append(f)
self._kwargs.append(kwargs) self._kwargs.append(kwargs)
......
...@@ -25,9 +25,9 @@ class StatCalculator(object): ...@@ -25,9 +25,9 @@ class StatCalculator(object):
Notes Notes
----- -----
- the memory usage of this object is constant, i.e. it does not increase - The memory usage of this object is constant, i.e. it does not increase
with the number of samples added with the number of samples added.
- the code computes the unbiased variance (which contains a `1./(n-1)` - The code computes the unbiased variance (which contains a `1./(n-1)`
term for `n` samples). term for `n` samples).
""" """
def __init__(self): def __init__(self):
......
...@@ -189,11 +189,10 @@ def testContractionOperator(spaces, wgt, dtype): ...@@ -189,11 +189,10 @@ def testContractionOperator(spaces, wgt, dtype):
def testDomainTupleFieldInserter(): def testDomainTupleFieldInserter():
domain = ift.DomainTuple.make((ift.UnstructuredDomain(12), target = ift.DomainTuple.make((ift.UnstructuredDomain([3, 2]),
ift.UnstructuredDomain(7),
ift.RGSpace([4, 22]))) ift.RGSpace([4, 22])))
new_space = ift.UnstructuredDomain(7) op = ift.DomainTupleFieldInserter(target, 1, (5,))
pos = (5,)
op = ift.DomainTupleFieldInserter(domain, new_space, 0, pos)
ift.extra.consistency_check(op) ift.extra.consistency_check(op)
......
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