Commit 550f3aa1 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'master' into fix_parameterized

# Conflicts:
#	test/common.py
parents dce9dcb5 dd68d8d0
Pipeline #13328 passed with stage
in 5 minutes and 21 seconds
.. currentmodule:: nifty
The ``ProjectionOperator`` class -- ...
.......................................
.. autoclass:: ProjectionOperator
:show-inheritance:
:members:
.. currentmodule:: nifty
The ``PropagatorOperator`` class -- ...
.......................................
.. autoclass:: PropagatorOperator
:show-inheritance:
:members:
.. currentmodule:: nifty
The ``ResponseOperator`` class -- A possible response implementation
....................................................................
.. autoclass:: ResponseOperator
:show-inheritance:
:members:
.. currentmodule:: nifty
The ``RGSpace`` class -- Regular Cartesian grids
................................................
.. autoclass:: RGSpace
:show-inheritance:
:members:
.. currentmodule:: nifty
The ``RGRGTransformation`` class -- A transformation routine
............................................................
.. autoclass:: RGRGTransformation
:show-inheritance:
:members:
.. currentmodule:: nifty
The ``SmoothingOperator`` class -- Smoothing fields
...................................................
.. autoclass:: SmoothingOperator
:show-inheritance:
:members:
Spaces
======
The :py:class:`Space` classes of NIFTY represent geometrical spaces approximated by grids in the computer environment. Each subclass of the base class corresponds to a specific grid type and replaces some of the inherited methods with its own methods that are unique to the respective grid. This framework ensures an abstract handling of spaces independent of the underlying geometrical grid and the grid's resolution.
Each instance of a :py:class:`Space` needs to capture all structural and dimensional specifics of the grid and all computationally relevant quantities such as the data type of associated field values. These parameters are stored as properties of an instance of the class at its initialization, and they do not need to be accessed explicitly by the user thereafter. This prevents the writing of grid or resolution dependent code.
Spatial symmetries of a system can be exploited by corresponding coordinate transformations. Often, transformations from one basis to its harmonic counterpart can greatly reduce the computational complexity of algorithms. The harmonic basis is defined by the eigenbasis of the Laplace operator; e.g., for a flat position space it is the Fourier basis. This conjugation of bases is implemented in NIFTY by distinguishing conjugate space classes, which can be obtained by the instance method *get_codomain* (and checked for by *check_codomain*). Moreover, transformations between conjugate spaces are performed automatically if required.
Space classes
-------------
Next to the generic :py:class:`Space` class, NIFTY has implementations of five subclasses, representing specific geometrical spaces and their discretizations.
.. toctree::
:maxdepth: 1
rg_space
hp_space
gl_space
lm_space
power_space
.. currentmodule:: nifty
The ``Space`` class -- The base Space object
--------------------------------------------
.. autoclass:: Space
:show-inheritance:
:members:
Transformations
---------------
NIFTY provides transformations
.. toctree::
:maxdepth: 1
rgrg_transformation
gllm_transformation
hplm_transformation
lmgl_transformation
lmhp_transformation
...@@ -269,7 +269,7 @@ class Field(Loggable, Versionable, object): ...@@ -269,7 +269,7 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods--- # ---Powerspectral methods---
def power_analyze(self, spaces=None, logarithmic=False, nbin=None, def power_analyze(self, spaces=None, logarithmic=False, nbin=None,
binbounds=None, decompose_power=True): binbounds=None, keep_phase_information=False):
""" Computes the square root power spectrum for a subspace of `self`. """ Computes the square root power spectrum for a subspace of `self`.
Creates a PowerSpace for the space addressed by `spaces` with the given Creates a PowerSpace for the space addressed by `spaces` with the given
...@@ -283,7 +283,6 @@ class Field(Loggable, Versionable, object): ...@@ -283,7 +283,6 @@ class Field(Loggable, Versionable, object):
spaces : int *optional* spaces : int *optional*
The subspace for which the powerspectrum shall be computed The subspace for which the powerspectrum shall be computed
(default : None). (default : None).
if spaces==None : Tries to synthesize for the whole domain
logarithmic : boolean *optional* logarithmic : boolean *optional*
True if the output PowerSpace should use logarithmic binning. True if the output PowerSpace should use logarithmic binning.
{default : False} {default : False}
...@@ -294,11 +293,16 @@ class Field(Loggable, Versionable, object): ...@@ -294,11 +293,16 @@ class Field(Loggable, Versionable, object):
binbounds : array-like *optional* binbounds : array-like *optional*
Inner bounds of the bins (default : None). Inner bounds of the bins (default : None).
if binbounds==None : bins are inferred. Overwrites nbins and log if binbounds==None : bins are inferred. Overwrites nbins and log
decompose_power : boolean, *optional* keep_phase_information : boolean, *optional*
Whether the analysed signal-space Field is intrinsically real or If False, return a real-valued result containing the power spectrum
complex and if the power spectrum shall therefore be computed of the input Field.
for the real and the imaginary part of the Field separately If True, return a complex-valued result whose real component
(default : True). contains the power spectrum computed from the real part of the
input Field, and whose imaginary component contains the power
spectrum computed from the imaginary part of the input Field.
The absolute value of this result should be identical to the output
of power_analyze with keep_phase_information=False.
(default : False).
Raise Raise
----- -----
...@@ -331,23 +335,47 @@ class Field(Loggable, Versionable, object): ...@@ -331,23 +335,47 @@ class Field(Loggable, Versionable, object):
# check if the `spaces` input is valid # check if the `spaces` input is valid
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain)) spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None: if spaces is None:
if len(self.domain) == 1: spaces = range(len(self.domain))
spaces = (0,)
else:
raise ValueError(
"Field has multiple spaces as domain "
"but `spaces` is None.")
if len(spaces) == 0: if len(spaces) == 0:
raise ValueError( raise ValueError(
"No space for analysis specified.") "No space for analysis specified.")
elif len(spaces) > 1:
raise ValueError(
"Conversion of only one space at a time is allowed.")
space_index = spaces[0] if keep_phase_information:
parts_val = self._hermitian_decomposition(
domain=self.domain,
val=self.val,
spaces=spaces,
domain_axes=self.domain_axes,
preserve_gaussian_variance=False)
parts = [self.copy_empty().set_val(part_val, copy=False)
for part_val in parts_val]
else:
parts = [self]
parts = [abs(part)**2 for part in parts]
for space_index in spaces:
parts = [self._single_power_analyze(
work_field=part,
space_index=space_index,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds)
for part in parts]
if keep_phase_information:
result_field = parts[0] + 1j*parts[1]
else:
result_field = parts[0]
return result_field
if not self.domain[space_index].harmonic: @classmethod
def _single_power_analyze(cls, work_field, space_index, logarithmic, nbin,
binbounds):
if not work_field.domain[space_index].harmonic:
raise ValueError( raise ValueError(
"The analyzed space must be harmonic.") "The analyzed space must be harmonic.")
...@@ -358,10 +386,10 @@ class Field(Loggable, Versionable, object): ...@@ -358,10 +386,10 @@ class Field(Loggable, Versionable, object):
# If it was complex, all the power is put into a real power spectrum. # If it was complex, all the power is put into a real power spectrum.
distribution_strategy = \ distribution_strategy = \
self.val.get_axes_local_distribution_strategy( work_field.val.get_axes_local_distribution_strategy(
self.domain_axes[space_index]) work_field.domain_axes[space_index])
harmonic_domain = self.domain[space_index] harmonic_domain = work_field.domain[space_index]
power_domain = PowerSpace(harmonic_partner=harmonic_domain, power_domain = PowerSpace(harmonic_partner=harmonic_domain,
distribution_strategy=distribution_strategy, distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin, logarithmic=logarithmic, nbin=nbin,
...@@ -371,35 +399,18 @@ class Field(Loggable, Versionable, object): ...@@ -371,35 +399,18 @@ class Field(Loggable, Versionable, object):
pindex = power_domain.pindex pindex = power_domain.pindex
rho = power_domain.rho rho = power_domain.rho
if decompose_power: power_spectrum = cls._calculate_power_spectrum(
hermitian_part, anti_hermitian_part = \ field_val=work_field.val,
harmonic_domain.hermitian_decomposition( pindex=pindex,
self.val, rho=rho,
axes=self.domain_axes[space_index]) axes=work_field.domain_axes[space_index])
[hermitian_power, anti_hermitian_power] = \
[self._calculate_power_spectrum(
x=part,
pindex=pindex,
rho=rho,
axes=self.domain_axes[space_index])
for part in [hermitian_part, anti_hermitian_part]]
power_spectrum = hermitian_power + 1j * anti_hermitian_power
else:
power_spectrum = self._calculate_power_spectrum(
x=self.val,
pindex=pindex,
rho=rho,
axes=self.domain_axes[space_index])
# create the result field and put power_spectrum into it # create the result field and put power_spectrum into it
result_domain = list(self.domain) result_domain = list(work_field.domain)
result_domain[space_index] = power_domain result_domain[space_index] = power_domain
result_dtype = power_spectrum.dtype result_dtype = power_spectrum.dtype
result_field = self.copy_empty( result_field = work_field.copy_empty(
domain=result_domain, domain=result_domain,
dtype=result_dtype, dtype=result_dtype,
distribution_strategy=power_spectrum.distribution_strategy) distribution_strategy=power_spectrum.distribution_strategy)
...@@ -407,17 +418,16 @@ class Field(Loggable, Versionable, object): ...@@ -407,17 +418,16 @@ class Field(Loggable, Versionable, object):
return result_field return result_field
def _calculate_power_spectrum(self, x, pindex, rho, axes=None): @classmethod
fieldabs = abs(x) def _calculate_power_spectrum(cls, field_val, pindex, rho, axes=None):
fieldabs **= 2
if axes is not None: if axes is not None:
pindex = self._shape_up_pindex( pindex = cls._shape_up_pindex(
pindex=pindex, pindex=pindex,
target_shape=x.shape, target_shape=field_val.shape,
target_strategy=x.distribution_strategy, target_strategy=field_val.distribution_strategy,
axes=axes) axes=axes)
power_spectrum = pindex.bincount(weights=fieldabs, power_spectrum = pindex.bincount(weights=field_val,
axis=axes) axis=axes)
if axes is not None: if axes is not None:
new_rho_shape = [1, ] * len(power_spectrum.shape) new_rho_shape = [1, ] * len(power_spectrum.shape)
...@@ -425,10 +435,10 @@ class Field(Loggable, Versionable, object): ...@@ -425,10 +435,10 @@ class Field(Loggable, Versionable, object):
rho = rho.reshape(new_rho_shape) rho = rho.reshape(new_rho_shape)
power_spectrum /= rho power_spectrum /= rho
power_spectrum **= 0.5
return power_spectrum return power_spectrum
def _shape_up_pindex(self, pindex, target_shape, target_strategy, axes): @staticmethod
def _shape_up_pindex(pindex, target_shape, target_strategy, axes):
if pindex.distribution_strategy not in \ if pindex.distribution_strategy not in \
DISTRIBUTION_STRATEGIES['global']: DISTRIBUTION_STRATEGIES['global']:
raise ValueError("pindex's distribution strategy must be " raise ValueError("pindex's distribution strategy must be "
...@@ -543,6 +553,8 @@ class Field(Loggable, Versionable, object): ...@@ -543,6 +553,8 @@ class Field(Loggable, Versionable, object):
# components # components
spec = self.val.get_full_data() spec = self.val.get_full_data()
spec = np.sqrt(spec)
for power_space_index in spaces: for power_space_index in spaces:
spec = self._spec_to_rescaler(spec, result_list, power_space_index) spec = self._spec_to_rescaler(spec, result_list, power_space_index)
local_rescaler = spec local_rescaler = spec
...@@ -561,10 +573,11 @@ class Field(Loggable, Versionable, object): ...@@ -561,10 +573,11 @@ class Field(Loggable, Versionable, object):
if real_signal: if real_signal:
result_val_list = [self._hermitian_decomposition( result_val_list = [self._hermitian_decomposition(
result_domain, result_domain,
result_val, result_val,
spaces, spaces,
result_list[0].domain_axes)[0] result_list[0].domain_axes,
preserve_gaussian_variance=True)[0]
for result_val in result_val_list] for result_val in result_val_list]
# store the result into the fields # store the result into the fields
...@@ -579,22 +592,23 @@ class Field(Loggable, Versionable, object): ...@@ -579,22 +592,23 @@ class Field(Loggable, Versionable, object):
return result return result
@staticmethod @staticmethod
def _hermitian_decomposition(domain, val, spaces, domain_axes): def _hermitian_decomposition(domain, val, spaces, domain_axes,
preserve_gaussian_variance=False):
# hermitianize for the first space # hermitianize for the first space
(h, a) = domain[spaces[0]].hermitian_decomposition( (h, a) = domain[spaces[0]].hermitian_decomposition(
val, val,
domain_axes[spaces[0]], domain_axes[spaces[0]],
preserve_gaussian_variance=True) preserve_gaussian_variance=preserve_gaussian_variance)
# hermitianize all remaining spaces using the iterative formula # hermitianize all remaining spaces using the iterative formula
for space in xrange(1, len(spaces)): for space in xrange(1, len(spaces)):
(hh, ha) = domain[space].hermitian_decomposition( (hh, ha) = domain[space].hermitian_decomposition(
h, h,
domain_axes[space], domain_axes[space],
preserve_gaussian_variance=True) preserve_gaussian_variance=False)
(ah, aa) = domain[space].hermitian_decomposition( (ah, aa) = domain[space].hermitian_decomposition(
a, a,
domain_axes[space], domain_axes[space],
preserve_gaussian_variance=True) preserve_gaussian_variance=False)
c = (hh - ha - ah + aa).conjugate() c = (hh - ha - ah + aa).conjugate()
full = (hh + ha + ah + aa) full = (hh + ha + ah + aa)
h = (full + c)/2. h = (full + c)/2.
...@@ -1142,8 +1156,8 @@ class Field(Loggable, Versionable, object): ...@@ -1142,8 +1156,8 @@ class Field(Loggable, Versionable, object):
""" """
return_field = self.copy_empty()
new_val = abs(self.get_val(copy=False)) new_val = abs(self.get_val(copy=False))
return_field = self.copy_empty(dtype=new_val.dtype)
return_field.set_val(new_val, copy=False) return_field.set_val(new_val, copy=False)
return return_field return return_field
......
...@@ -37,7 +37,7 @@ class GLLMTransformation(SlicingTransformation): ...@@ -37,7 +37,7 @@ class GLLMTransformation(SlicingTransformation):
if module != 'pyHealpix': if module != 'pyHealpix':
raise ValueError("Unsupported SHT module.") raise ValueError("Unsupported SHT module.")
if 'pyHealpix' not in gdi: if pyHealpix is None:
raise ImportError( raise ImportError(
"The module pyHealpix is needed but not available.") "The module pyHealpix is needed but not available.")
......
...@@ -38,7 +38,7 @@ class HPLMTransformation(SlicingTransformation): ...@@ -38,7 +38,7 @@ class HPLMTransformation(SlicingTransformation):
if module != 'pyHealpix': if module != 'pyHealpix':
raise ValueError("Unsupported SHT module.") raise ValueError("Unsupported SHT module.")
if 'pyHealpix' not in gdi: if pyHealpix is None:
raise ImportError( raise ImportError(
"The module pyHealpix is needed but not available") "The module pyHealpix is needed but not available")
......
...@@ -37,7 +37,7 @@ class LMGLTransformation(SlicingTransformation): ...@@ -37,7 +37,7 @@ class LMGLTransformation(SlicingTransformation):
if module != 'pyHealpix': if module != 'pyHealpix':
raise ValueError("Unsupported SHT module.") raise ValueError("Unsupported SHT module.")
if 'pyHealpix' not in gdi: if pyHealpix is None:
raise ImportError( raise ImportError(
"The module pyHealpix is needed but not available.") "The module pyHealpix is needed but not available.")
......
...@@ -36,7 +36,7 @@ class LMHPTransformation(SlicingTransformation): ...@@ -36,7 +36,7 @@ class LMHPTransformation(SlicingTransformation):
if module != 'pyHealpix': if module != 'pyHealpix':
raise ValueError("Unsupported SHT module.") raise ValueError("Unsupported SHT module.")
if gdi.get('pyHealpix') is None: if pyHealpix is None:
raise ImportError( raise ImportError(
"The module pyHealpix is needed but not available.") "The module pyHealpix is needed but not available.")
......
...@@ -13,7 +13,7 @@ class GLMollweide(Heatmap): ...@@ -13,7 +13,7 @@ class GLMollweide(Heatmap):
def __init__(self, data, xsize=800, color_map=None, def __init__(self, data, xsize=800, color_map=None,
webgl=False, smoothing=False): webgl=False, smoothing=False):
# smoothing 'best', 'fast', False # smoothing 'best', 'fast', False
if 'pyHealpix' not in gdi: if pyHealpix is None:
raise ImportError( raise ImportError(
"The module pyHealpix is needed but not available.") "The module pyHealpix is needed but not available.")
self.xsize = xsize self.xsize = xsize
......
...@@ -12,7 +12,7 @@ pyHealpix = gdi.get('pyHealpix') ...@@ -12,7 +12,7 @@ pyHealpix = gdi.get('pyHealpix')
class HPMollweide(Heatmap): class HPMollweide(Heatmap):
def __init__(self, data, xsize=800, color_map=None, webgl=False, def __init__(self, data, xsize=800, color_map=None, webgl=False,
smoothing=False): # smoothing 'best', 'fast', False smoothing=False): # smoothing 'best', 'fast', False
if 'pyHealpix' not in gdi: if pyHealpix is None:
raise ImportError( raise ImportError(
"The module pyHealpix is needed but not available.") "The module pyHealpix is needed but not available.")
self.xsize = xsize self.xsize = xsize
......
...@@ -31,7 +31,7 @@ class PlotterBase(Loggable, object): ...@@ -31,7 +31,7 @@ class PlotterBase(Loggable, object):
__metaclass__ = abc.ABCMeta __metaclass__ = abc.ABCMeta
def __init__(self, interactive=False, path='.', title=""): def __init__(self, interactive=False, path='.', title=""):
if 'plotly' not in gdi: if plotly is None:
raise ImportError("The module plotly is needed but not available.") raise ImportError("The module plotly is needed but not available.")
self.interactive = interactive self.interactive = interactive
self.path = path self.path = path
......
...@@ -45,7 +45,7 @@ class GLSpace(Space): ...@@ -45,7 +45,7 @@ class GLSpace(Space):
Number of latitudinal bins (or rings) that are used for this Number of latitudinal bins (or rings) that are used for this
pixelization. pixelization.
nlon : int, *optional* nlon : int, *optional*
Number of longditudinal bins that are used for this pixelization. Number of longitudinal bins that are used for this pixelization.
Attributes Attributes
---------- ----------
...@@ -57,7 +57,7 @@ class GLSpace(Space): ...@@ -57,7 +57,7 @@ class GLSpace(Space):
Number of latitudinal bins (or rings) that are used for this Number of latitudinal bins (or rings) that are used for this
pixelization. pixelization.
nlon : int nlon : int
Number of longditudinal bins that are used for this pixelization. Number of longitudinal bins that are used for this pixelization.
total_volume : np.float total_volume : np.float
The total volume of the space. The total volume of the space.
shape : tuple of np.ints shape : tuple of np.ints
...@@ -89,7 +89,7 @@ class GLSpace(Space): ...@@ -89,7 +89,7 @@ class GLSpace(Space):
# ---Overwritten properties and methods--- # ---Overwritten properties and methods---
def __init__(self, nlat, nlon=None): def __init__(self, nlat, nlon=None):
if 'pyHealpix' not in gdi: if pyHealpix is None:
raise ImportError( raise ImportError(
"The module pyHealpix is needed but not available.") "The module pyHealpix is needed but not available.")
...@@ -163,7 +163,7 @@ class GLSpace(Space): ...@@ -163,7 +163,7 @@ class GLSpace(Space):
@property @property
def nlon(self): def nlon(self):
""" Number of longditudinal bins that are used for this pixelization. """ Number of longitudinal bins that are used for this pixelization.
""" """
return self._nlon return self._nlon
......
...@@ -32,7 +32,7 @@ def create_power_operator(domain, power_spectrum, dtype=None, ...@@ -32,7 +32,7 @@ def create_power_operator(domain, power_spectrum, dtype=None,
Parameters Parameters
---------- ----------
domain : DomainObject domain : DomainObject
Domain over which the power operator shall live. Domain over which the power operator shall live.
power_spectrum : (array-like, method) power_spectrum : (array-like, method)
An array-like object, or a method that implements the square root An array-like object, or a method that implements the square root
of a power spectrum as a function of k. of a power spectrum as a function of k.
...@@ -54,7 +54,6 @@ def create_power_operator(domain, power_spectrum, dtype=None, ...@@ -54,7 +54,6 @@ def create_power_operator(domain, power_spectrum, dtype=None,
distribution_strategy=distribution_strategy) distribution_strategy=distribution_strategy)
fp = Field(power_domain, val=power_spectrum, dtype=dtype, fp = Field(power_domain, val=power_spectrum, dtype=dtype,
distribution_strategy=distribution_strategy) distribution_strategy=distribution_strategy)
fp **= 2
f = fp.power_synthesize(mean=1, std=0, real_signal=False) f = fp.power_synthesize(mean=1, std=0, real_signal=False)
return DiagonalOperator(domain, diagonal=f, bare=True) return DiagonalOperator(domain, diagonal=f, bare=True)
...@@ -16,47 +16,15 @@ ...@@ -1