Commit 95df1f22 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'master' into now_DescentMinimizer_and_bugfix

parents 2e48be5f 54f203dd
Pipeline #12262 failed with stage
in 5 minutes and 42 seconds
......@@ -107,7 +107,7 @@ if __name__ == "__main__":
# callback=distance_measure,
# max_history_length=3)
m0 = Field(s_space, val=1)
m0 = Field(s_space, val=1.)
energy = WienerFilterEnergy(position=m0, D=D, j=j)
......
......@@ -41,6 +41,20 @@ class DomainObject(Versionable, Loggable, object):
return result_hash
def __eq__(self, x):
""" Checks if two domain_objects are equal.
Parameters
----------
x: domain_object
The domain_object `self` is compared to.
Returns
-------
bool
True if `self` and x describe the same manifold.
"""
if isinstance(x, type(self)):
for key in vars(self).keys():
item1 = vars(self)[key]
......@@ -58,23 +72,137 @@ class DomainObject(Versionable, Loggable, object):
@abc.abstractproperty
def shape(self):
""" Returns the shape of the underlying array-like object.
Returns
-------
tuple of ints
The shape of the underlying array-like object.
Raises
------
NotImplementedError
If called for this abstract class.
"""
raise NotImplementedError(
"There is no generic shape for DomainObject.")
@abc.abstractproperty
def dim(self):
""" Returns the number of pixel-dimensions the object has.
Returns
-------
int
An Integer representing the number of pixels the discretized
manifold has.
Raises
------
NotImplementedError
If called for this abstract class.
"""
raise NotImplementedError(
"There is no generic dim for DomainObject.")
@abc.abstractmethod
def weight(self, x, power=1, axes=None, inplace=False):
""" Weights the field on this domain with the space's volume-weights.
Weights hereby refer to integration weights, as they appear in
discretized integrals. Per default, this function mutliplies each bin
of the field x by its volume, which lets it behave like a density
(top form). However, different powers of the volume can be applied
with the power parameter. The axes parameter specifies which of the
field array's indices correspond to this domain.
Parameters
----------
x : distributed_data_object
The fields data array.
power : int, *optional*
The power to which the volume-weight is raised (default: 1).
axes : {int, tuple}, *optional*
Specifies the axes of x which represent this domain
(default: None).
If axes==None:
weighting is applied with respect to all axes
inplace : bool, *optional*
If this is True, the weighting is done on the values of x,
if it is False, x is not modified and this method returns a
weighted copy of x (default: False).
Returns
-------
distributed_data_object
A weighted version of x, with volume-weights raised to the
given power.
Raises
------
NotImplementedError
If called for this abstract class.
"""
raise NotImplementedError(
"There is no generic weight-method for DomainObject.")
def pre_cast(self, x, axes=None):
def pre_cast(self, x, axes):
""" Casts input for Field.val before Field performs the cast.
Parameters
----------
x : {array-like, castable}
an array-like object or anything that can be cast to arrays.
axes : tuple of ints
Specifies the axes of x which correspond to this domain.
Returns
-------
{array-like, castable}
Processed input where casting that needs Space-specific knowledge
(for example location of pixels on the manifold) was performed.
See Also
--------
post_cast
Notes
-----
Usually returns x, except if a power spectrum is given to a
PowerSpace, where this spectrum is evaluated at the power indices.
"""
return x
def post_cast(self, x, axes=None):
def post_cast(self, x, axes):
""" Performs casting operations that are done after Field's cast.
Parameters
----------
x : {array-like, castable}
an array-like object or anything that can be cast to arrays.
axes : tuple of ints
Specifies the axes of x which correspond to this domain.
See Also
--------
pre_cast
Returns
-------
distributed_data_object
Processed input where casting that needs Space-specific knowledge
(for example location of pixels on the manifold) was performed.
"""
return x
# ---Serialization---
......
......@@ -167,8 +167,8 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods---
def power_analyze(self, spaces=None, log=False, nbin=None, binbounds=None,
real_signal=True):
def power_analyze(self, spaces=None, logarithmic=False, nbin=None, binbounds=None,
decompose_power=False):
# check if all spaces in `self.domain` are either harmonic or
# power_space instances
for sp in self.domain:
......@@ -210,18 +210,18 @@ class Field(Loggable, Versionable, object):
self.val.get_axes_local_distribution_strategy(
self.domain_axes[space_index])
harmonic_domain = self.domain[space_index]
power_domain = PowerSpace(harmonic_domain=harmonic_domain,
harmonic_partner = self.domain[space_index]
power_domain = PowerSpace(harmonic_partner=harmonic_partner,
distribution_strategy=distribution_strategy,
log=log, nbin=nbin, binbounds=binbounds)
logarithmic=logarithmic, nbin=nbin, binbounds=binbounds)
# extract pindex and rho from power_domain
pindex = power_domain.pindex
rho = power_domain.rho
if real_signal:
if decompose_power:
hermitian_part, anti_hermitian_part = \
harmonic_domain.hermitian_decomposition(
harmonic_partner.hermitian_decomposition(
self.val,
axes=self.domain_axes[space_index])
......@@ -245,7 +245,7 @@ class Field(Loggable, Versionable, object):
result_domain = list(self.domain)
result_domain[space_index] = power_domain
if real_signal:
if decompose_power:
result_dtype = np.complex
else:
result_dtype = np.float
......@@ -303,8 +303,8 @@ class Field(Loggable, Versionable, object):
return result_obj
def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
mean=None, std=None):
def power_synthesize(self, spaces=None, real_power=True,
real_signal=False, mean=None, std=None):
# check if the `spaces` input is valid
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
......@@ -322,8 +322,8 @@ class Field(Loggable, Versionable, object):
result_domain = list(self.domain)
for power_space_index in spaces:
power_space = self.domain[power_space_index]
harmonic_domain = power_space.harmonic_domain
result_domain[power_space_index] = harmonic_domain
harmonic_partner = power_space.harmonic_partner
result_domain[power_space_index] = harmonic_partner
# create random samples: one or two, depending on whether the
# power spectrum is real or complex
......@@ -365,8 +365,8 @@ class Field(Loggable, Versionable, object):
if real_signal:
for power_space_index in spaces:
harmonic_domain = result_domain[power_space_index]
result_val_list = [harmonic_domain.hermitian_decomposition(
harmonic_partner = result_domain[power_space_index]
result_val_list = [harmonic_partner.hermitian_decomposition(
result_val,
axes=result.domain_axes[power_space_index],
preserve_gaussian_variance=True)[0]
......
......@@ -19,6 +19,7 @@
import numpy as np
from itertools import product
def get_slice_list(shape, axes):
"""
Helper function which generates slice list(s) to traverse over all
......@@ -65,8 +66,7 @@ def get_slice_list(shape, axes):
return
def cast_axis_to_tuple(axis, length):
def cast_axis_to_tuple(axis, length=None):
if axis is None:
return None
try:
......@@ -78,16 +78,17 @@ def cast_axis_to_tuple(axis, length):
raise TypeError(
"Could not convert axis-input to tuple of ints")
# shift negative indices to positive ones
axis = tuple(item if (item >= 0) else (item + length) for item in axis)
if length is not None:
# shift negative indices to positive ones
axis = tuple(item if (item >= 0) else (item + length) for item in axis)
# Deactivated this, in order to allow for the ComposedOperator
# remove duplicate entries
# axis = tuple(set(axis))
# Deactivated this, in order to allow for the ComposedOperator
# remove duplicate entries
# axis = tuple(set(axis))
# assert that all entries are elements in [0, length]
for elem in axis:
assert (0 <= elem < length)
# assert that all entries are elements in [0, length]
for elem in axis:
assert (0 <= elem < length)
return axis
......
......@@ -21,7 +21,9 @@ from nifty.operators.linear_operator import LinearOperator
class ComposedOperator(LinearOperator):
# ---Overwritten properties and methods---
def __init__(self, operators):
def __init__(self, operators, default_spaces=None):
super(ComposedOperator, self).__init__(default_spaces)
self._operator_store = ()
for op in operators:
if not isinstance(op, LinearOperator):
......
......@@ -30,9 +30,10 @@ class DiagonalOperator(EndomorphicOperator):
# ---Overwritten properties and methods---
def __init__(self, domain=(),
diagonal=None, bare=False, copy=True,
distribution_strategy=None):
def __init__(self, domain=(), diagonal=None, bare=False, copy=True,
distribution_strategy=None, default_spaces=None):
super(DiagonalOperator, self).__init__(default_spaces)
self._domain = self._parse_domain(domain)
if distribution_strategy is None:
......
......@@ -34,13 +34,72 @@ from transformations import RGRGTransformation,\
class FFTOperator(LinearOperator):
"""Transforms between a pair of position and harmonic domains.
Built-in domain pairs are
- a harmonic and a non-harmonic RGSpace (with matching distances)
- a HPSpace and a LMSpace
- a GLSpace and a LMSpace
Within a domain pair, both orderings are possible.
The operator provides a "times" and an "adjoint_times" operation.
For a pair of RGSpaces, the "adjoint_times" operation is equivalent to
"inverse_times"; for the sphere-related domains this is not the case, since
the operator matrix is not square.
Parameters
----------
domain: Space or single-element tuple of Spaces
The domain of the data that is input by "times" and output by
"adjoint_times".
target: Space or single-element tuple of Spaces (optional)
The domain of the data that is output by "times" and input by
"adjoint_times".
If omitted, a co-domain will be chosen automatically.
Whenever "domain" is an RGSpace, the codomain (and its parameters) are
uniquely determined (except for "zerocenter").
For GLSpace, HPSpace, and LMSpace, a sensible (but not unique)
co-domain is chosen that should work satisfactorily in most situations,
but for full control, the user should explicitly specify a codomain.
module: String (optional)
Software module employed for carrying out the transform operations.
For RGSpace pairs this can be "numpy" or "fftw", where "numpy" is
always available, but "fftw" offers higher performance and
parallelization. For sphere-related domains, only "pyHealpix" is
available. If omitted, "fftw" is selected for RGSpaces if available,
else "numpy"; on the sphere the default is "pyHealpix".
domain_dtype: data type (optional)
Data type of the fields that go into "times" and come out of
"adjoint_times". Default is "numpy.complex".
target_dtype: data type (optional)
Data type of the fields that go into "adjoint_times" and come out of
"times". Default is "numpy.complex".
Attributes
----------
domain: Tuple of Spaces (with one entry)
The domain of the data that is input by "times" and output by
"adjoint_times".
target: Tuple of Spaces (with one entry)
The domain of the data that is output by "times" and input by
"adjoint_times".
unitary: bool
Returns True if the operator is unitary (currently only the case if
the domain and codomain are RGSpaces), else False.
Raises
------
ValueError:
if "domain" or "target" are not of the proper type.
"""
# ---Class attributes---
default_codomain_dictionary = {RGSpace: RGSpace,
HPSpace: LMSpace,
GLSpace: LMSpace,
LMSpace: HPSpace,
LMSpace: GLSpace,
}
transformation_dictionary = {(RGSpace, RGSpace): RGRGTransformation,
......@@ -52,8 +111,9 @@ class FFTOperator(LinearOperator):
# ---Overwritten properties and methods---
def __init__(self, domain=(), target=None, module=None,
domain_dtype=None, target_dtype=None):
def __init__(self, domain, target=None, module=None,
domain_dtype=None, target_dtype=None, default_spaces=None):
super(FFTOperator, self).__init__(default_spaces)
# Initialize domain and target
......@@ -83,8 +143,8 @@ class FFTOperator(LinearOperator):
# Store the dtype information
if domain_dtype is None:
self.logger.info("Setting domain_dtype to np.float.")
self.domain_dtype = np.float
self.logger.info("Setting domain_dtype to np.complex.")
self.domain_dtype = np.complex
else:
self.domain_dtype = np.dtype(domain_dtype)
......@@ -114,11 +174,11 @@ class FFTOperator(LinearOperator):
result_field = x.copy_empty(domain=result_domain,
dtype=self.target_dtype)
result_field.set_val(new_val=new_val, copy=False)
result_field.set_val(new_val=new_val, copy=True)
return result_field
def _inverse_times(self, x, spaces):
def _adjoint_times(self, x, spaces):
spaces = utilities.cast_axis_to_tuple(spaces, len(x.domain))
if spaces is None:
# this case means that x lives on only one space, which is
......@@ -138,7 +198,7 @@ class FFTOperator(LinearOperator):
result_field = x.copy_empty(domain=result_domain,
dtype=self.domain_dtype)
result_field.set_val(new_val=new_val, copy=False)
result_field.set_val(new_val=new_val, copy=True)
return result_field
......@@ -154,12 +214,38 @@ class FFTOperator(LinearOperator):
@property
def unitary(self):
return True
return (self._forward_transformation.unitary and
self._backward_transformation.unitary)
# ---Added properties and methods---
@classmethod
def get_default_codomain(cls, domain):
"""Returns a codomain to the given domain.
Parameters
----------
domain: Space
An instance of RGSpace, HPSpace, GLSpace or LMSpace.
Returns
-------
target: Space
A (more or less perfect) counterpart to "domain" with respect
to a FFT operation.
Whenever "domain" is an RGSpace, the codomain (and its parameters)
are uniquely determined (except for "zerocenter").
For GLSpace, HPSpace, and LMSpace, a sensible (but not unique)
co-domain is chosen that should work satisfactorily in most
situations. For full control however, the user should not rely on
this method.
Raises
------
ValueError:
if no default codomain is defined for "domain".
"""
domain_class = domain.__class__
try:
codomain_class = cls.default_codomain_dictionary[domain_class]
......
......@@ -45,6 +45,10 @@ class GLLMTransformation(SlicingTransformation):
# ---Mandatory properties and methods---
@property
def unitary(self):
return False
@classmethod
def get_codomain(cls, domain):
"""
......@@ -85,13 +89,13 @@ class GLLMTransformation(SlicingTransformation):
mmax = codomain.mmax
if lmax != mmax:
cls.Logger.warn("Unrecommended: codomain has lmax != mmax.")
cls.logger.warn("Unrecommended: codomain has lmax != mmax.")
if lmax != nlat - 1:
cls.Logger.warn("Unrecommended: codomain has lmax != nlat - 1.")
cls.logger.warn("Unrecommended: codomain has lmax != nlat - 1.")
if nlon != 2*nlat - 1:
cls.Logger.warn("Unrecommended: domain has nlon != 2*nlat - 1.")
cls.logger.warn("Unrecommended: domain has nlon != 2*nlat - 1.")
super(GLLMTransformation, cls).check_codomain(domain, codomain)
......
......@@ -46,6 +46,10 @@ class HPLMTransformation(SlicingTransformation):
# ---Mandatory properties and methods---
@property
def unitary(self):
return False
@classmethod
def get_codomain(cls, domain):
"""
......@@ -83,7 +87,7 @@ class HPLMTransformation(SlicingTransformation):
nside = domain.nside
if lmax != 2*nside:
cls.Logger.warn("Unrecommended: lmax != 2*nside.")
cls.logger.warn("Unrecommended: lmax != 2*nside.")
super(HPLMTransformation, cls).check_codomain(domain, codomain)
......@@ -98,7 +102,7 @@ class HPLMTransformation(SlicingTransformation):
if issubclass(inp.dtype.type, np.complexfloating):
[resultReal,
resultImag] = [pyHealpix.map2alm_iter(x, lmax, mmax, 3)
resultImag] = [pyHealpix.map2alm(x, lmax, mmax)
for x in (inp.real, inp.imag)]
[resultReal,
......@@ -108,7 +112,7 @@ class HPLMTransformation(SlicingTransformation):
result = self._combine_complex_result(resultReal, resultImag)
else:
result = pyHealpix.map2alm_iter(inp, lmax, mmax, 3)
result = pyHealpix.map2alm(inp, lmax, mmax)
result = lm_transformation_helper.buildIdx(result, lmax=lmax)
return result
......@@ -45,6 +45,10 @@ class LMGLTransformation(SlicingTransformation):
# ---Mandatory properties and methods---
@property
def unitary(self):
return False
@classmethod
def get_codomain(cls, domain):
"""
......@@ -91,13 +95,13 @@ class LMGLTransformation(SlicingTransformation):
mmax = domain.mmax
if lmax != mmax:
cls.Logger.warn("Unrecommended: codomain has lmax != mmax.")
cls.logger.warn("Unrecommended: codomain has lmax != mmax.")
if nlat != lmax + 1:
cls.Logger.warn("Unrecommended: codomain has nlat != lmax + 1.")
cls.logger.warn("Unrecommended: codomain has nlat != lmax + 1.")
if nlon != 2*lmax + 1:
cls.Logger.warn("Unrecommended: domain has nlon != 2*lmax + 1.")
cls.logger.warn("Unrecommended: domain has nlon != 2*lmax + 1.")
super(LMGLTransformation, cls).check_codomain(domain, codomain)
......
......@@ -44,6 +44,10 @@ class LMHPTransformation(SlicingTransformation):
# ---Mandatory properties and methods---
@property
def unitary(self):
return False
@classmethod
def get_codomain(cls, domain):
"""
......@@ -85,7 +89,7 @@ class LMHPTransformation(SlicingTransformation):
lmax = domain.lmax
if lmax != 2*nside:
cls.Logger.warn("Unrecommended: lmax != 2*nside.")
cls.logger.warn("Unrecommended: lmax != 2*nside.")
super(LMHPTransformation, cls).check_codomain(domain, codomain)
......
......@@ -23,6 +23,9 @@ from nifty import RGSpace, nifty_configuration
class RGRGTransformation(Transformation):
# ---Overwritten properties and methods---
def __init__(self, domain, codomain=None, module=None):
super(RGRGTransformation, self).__init__(domain, codomain, module)
......@@ -42,6 +45,12 @@ class RGRGTransformation(Transformation):
else:
raise ValueError('Unsupported FFT module:' + module)
# ---Mandatory properties and methods---
@property
def unitary(self):
return True
@classmethod
def get_codomain(cls, domain, zerocenter=None):
"""
......
......@@ -37,6 +37,10 @@ class Transformation(Loggable, object):
self.domain = domain
self.codomain = codomain
@abc.abstractproperty