Commit a0b26666 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'master' into issue77_2

parents 93db476d 42951308
......@@ -10,10 +10,10 @@ rank = comm.rank
if __name__ == "__main__":
distribution_strategy = 'fftw'
distribution_strategy = 'not'
# Setting up the geometry
s_space = RGSpace([512, 512], dtype=np.float64)
s_space = RGSpace([512, 512])
fft = FFTOperator(s_space)
h_space = fft.target[0]
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
......@@ -54,4 +54,3 @@ if __name__ == "__main__":
# pl.plot([go.Heatmap(z=d_data)], filename='data.html')
# pl.plot([go.Heatmap(z=m_data)], filename='map.html')
# pl.plot([go.Heatmap(z=ss_data)], filename='map_orig.html')
#
......@@ -53,10 +53,10 @@ class WienerFilterEnergy(Energy):
if __name__ == "__main__":
distribution_strategy = 'fftw'
distribution_strategy = 'not'
# Set up spaces and fft transformation
s_space = RGSpace([512, 512], dtype=np.float)
s_space = RGSpace([512, 512])
fft = FFTOperator(s_space)
h_space = fft.target[0]
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
......
......@@ -18,8 +18,6 @@
import abc
import numpy as np
from keepers import Loggable,\
Versionable
......@@ -27,8 +25,7 @@ from keepers import Loggable,\
class DomainObject(Versionable, Loggable, object):
__metaclass__ = abc.ABCMeta
def __init__(self, dtype):
self._dtype = np.dtype(dtype)
def __init__(self):
self._ignore_for_hash = []
def __hash__(self):
......@@ -57,10 +54,6 @@ class DomainObject(Versionable, Loggable, object):
def __ne__(self, x):
return not self.__eq__(x)
@property
def dtype(self):
return self._dtype
@abc.abstractproperty
def shape(self):
raise NotImplementedError(
......@@ -85,10 +78,9 @@ class DomainObject(Versionable, Loggable, object):
# ---Serialization---
def _to_hdf5(self, hdf5_group):
hdf5_group.attrs['dtype'] = self.dtype.name
return None
@classmethod
def _from_hdf5(cls, hdf5_group, repository):
result = cls(dtype=np.dtype(hdf5_group.attrs['dtype']))
result = cls()
return result
......@@ -45,8 +45,7 @@ class Field(Loggable, Versionable, object):
self.domain_axes = self._get_axes_tuple(self.domain)
self.dtype = self._infer_dtype(dtype=dtype,
val=val,
domain=self.domain)
val=val)
self.distribution_strategy = self._parse_distribution_strategy(
distribution_strategy=distribution_strategy,
......@@ -86,18 +85,19 @@ class Field(Loggable, Versionable, object):
axes_list += [tuple(l)]
return tuple(axes_list)
def _infer_dtype(self, dtype, val, domain):
def _infer_dtype(self, dtype, val):
if dtype is None:
if isinstance(val, Field) or \
isinstance(val, distributed_data_object):
try:
dtype = val.dtype
dtype_tuple = (np.dtype(gc['default_field_dtype']),)
except AttributeError:
try:
if val is None:
raise TypeError
dtype = np.result_type(val)
except(TypeError):
dtype = np.dtype(gc['default_field_dtype'])
else:
dtype_tuple = (np.dtype(dtype),)
if domain is not None:
dtype_tuple += tuple(np.dtype(sp.dtype) for sp in domain)
dtype = reduce(lambda x, y: np.result_type(x, y), dtype_tuple)
dtype = np.dtype(dtype)
return dtype
......@@ -303,59 +303,41 @@ class Field(Loggable, Versionable, object):
return result_obj
def power_synthesize(self, spaces=None, real_signal=True,
def power_synthesize(self, spaces=None, real_power=True, real_signal=True,
mean=None, std=None):
# check if all spaces in `self.domain` are either of signal-type or
# power_space instances
for sp in self.domain:
if not sp.harmonic and not isinstance(sp, PowerSpace):
self.logger.info(
"Field has a space in `domain` which is neither "
"harmonic nor a PowerSpace.")
# check if the `spaces` input is valid
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
if spaces is None:
if len(self.domain) == 1:
spaces = (0,)
else:
raise ValueError(
"Field has multiple spaces as domain "
"but `spaces` is None.")
if len(spaces) == 0:
raise ValueError(
"No space for synthesis specified.")
elif len(spaces) > 1:
raise ValueError(
"Conversion of only one space at a time is allowed.")
if spaces is None:
spaces = range(len(self.domain))
power_space_index = spaces[0]
power_domain = self.domain[power_space_index]
if not isinstance(power_domain, PowerSpace):
raise ValueError(
"A PowerSpace is needed for field synthetization.")
for power_space_index in spaces:
power_space = self.domain[power_space_index]
if not isinstance(power_space, PowerSpace):
raise ValueError("A PowerSpace is needed for field "
"synthetization.")
# create the result domain
result_domain = list(self.domain)
harmonic_domain = power_domain.harmonic_domain
result_domain[power_space_index] = harmonic_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
# create random samples: one or two, depending on whether the
# power spectrum is real or complex
if issubclass(power_domain.dtype.type, np.complexfloating):
result_list = [None, None]
else:
if real_power:
result_list = [None]
else:
result_list = [None, None]
result_list = [self.__class__.from_random(
'normal',
mean=mean,
std=std,
domain=result_domain,
dtype=harmonic_domain.dtype,
dtype=np.complex,
distribution_strategy=self.distribution_strategy)
for x in result_list]
......@@ -363,18 +345,51 @@ class Field(Loggable, Versionable, object):
# processing without killing the fields.
# if the signal-space field should be real, hermitianize the field
# components
spec = self.val.get_full_data()
for power_space_index in spaces:
spec = self._spec_to_rescaler(spec, result_list, power_space_index)
local_rescaler = spec
result_val_list = [x.val for x in result_list]
# apply the rescaler to the random fields
result_val_list[0].apply_scalar_function(
lambda x: x * local_rescaler.real,
inplace=True)
if not real_power:
result_val_list[1].apply_scalar_function(
lambda x: x * local_rescaler.imag,
inplace=True)
if real_signal:
result_val_list = [harmonic_domain.hermitian_decomposition(
x.val,
axes=x.domain_axes[power_space_index],
for power_space_index in spaces:
harmonic_domain = result_domain[power_space_index]
result_val_list = [harmonic_domain.hermitian_decomposition(
result_val,
axes=result.domain_axes[power_space_index],
preserve_gaussian_variance=True)[0]
for x in result_list]
for (result, result_val)
in zip(result_list, result_val_list)]
# store the result into the fields
[x.set_val(new_val=y, copy=False) for x, y in
zip(result_list, result_val_list)]
if real_power:
result = result_list[0]
else:
result_val_list = [x.val for x in result_list]
result = result_list[0] + 1j*result_list[1]
return result
def _spec_to_rescaler(self, spec, result_list, power_space_index):
power_space = self.domain[power_space_index]
# weight the random fields with the power spectrum
# therefore get the pindex from the power space
pindex = power_domain.pindex
pindex = power_space.pindex
# take the local data from pindex. This data must be compatible to the
# local data of the field given the slice of the PowerSpace
local_distribution_strategy = \
......@@ -390,34 +405,12 @@ class Field(Loggable, Versionable, object):
# power spectrum into the appropriate places of the pindex array.
# Do this for every 'pindex-slice' in parallel using the 'slice(None)'s
local_pindex = pindex.get_local_data(copy=False)
full_spec = self.val.get_full_data()
local_blow_up = [slice(None)]*len(self.shape)
local_blow_up[self.domain_axes[power_space_index][0]] = local_pindex
# here, the power_spectrum is distributed into the new shape
local_rescaler = full_spec[local_blow_up]
# apply the rescaler to the random fields
result_val_list[0].apply_scalar_function(
lambda x: x * local_rescaler.real,
inplace=True)
if issubclass(power_domain.dtype.type, np.complexfloating):
result_val_list[1].apply_scalar_function(
lambda x: x * local_rescaler.imag,
inplace=True)
# store the result into the fields
[x.set_val(new_val=y, copy=False) for x, y in
zip(result_list, result_val_list)]
if issubclass(power_domain.dtype.type, np.complexfloating):
result = result_list[0] + 1j*result_list[1]
else:
result = result_list[0]
return result
local_rescaler = spec[local_blow_up]
return local_rescaler
# ---Properties---
......
......@@ -105,10 +105,10 @@ class DiagonalOperator(EndomorphicOperator):
return self._implemented
@property
def symmetric(self):
if self._symmetric is None:
self._symmetric = (self._diagonal.val.imag == 0).all()
return self._symmetric
def self_adjoint(self):
if self._self_adjoint is None:
self._self_adjoint = (self._diagonal.val.imag == 0).all()
return self._self_adjoint
@property
def unitary(self):
......@@ -155,8 +155,8 @@ class DiagonalOperator(EndomorphicOperator):
# Otherwise, inplace weightening would change the external field
f.weight(inplace=copy, power=-1)
# Reset the symmetric property:
self._symmetric = None
# Reset the self_adjoint property:
self._self_adjoint = None
# Reset the unitarity property
self._unitary = None
......
......@@ -26,7 +26,7 @@ class EndomorphicOperator(LinearOperator):
# ---Overwritten properties and methods---
def inverse_times(self, x, spaces=None):
if self.symmetric and self.unitary:
if self.self_adjoint and self.unitary:
return self.times(x, spaces)
else:
return super(EndomorphicOperator, self).inverse_times(
......@@ -34,7 +34,7 @@ class EndomorphicOperator(LinearOperator):
spaces=spaces)
def adjoint_times(self, x, spaces=None):
if self.symmetric:
if self.self_adjoint:
return self.times(x, spaces)
else:
return super(EndomorphicOperator, self).adjoint_times(
......@@ -42,7 +42,7 @@ class EndomorphicOperator(LinearOperator):
spaces=spaces)
def adjoint_inverse_times(self, x, spaces=None):
if self.symmetric:
if self.self_adjoint:
return self.inverse_times(x, spaces)
else:
return super(EndomorphicOperator, self).adjoint_inverse_times(
......@@ -50,7 +50,7 @@ class EndomorphicOperator(LinearOperator):
spaces=spaces)
def inverse_adjoint_times(self, x, spaces=None):
if self.symmetric:
if self.self_adjoint:
return self.inverse_times(x, spaces)
else:
return super(EndomorphicOperator, self).inverse_adjoint_times(
......@@ -66,5 +66,5 @@ class EndomorphicOperator(LinearOperator):
# ---Added properties and methods---
@abc.abstractproperty
def symmetric(self):
def self_adjoint(self):
raise NotImplementedError
......@@ -83,12 +83,14 @@ class FFTOperator(LinearOperator):
# Store the dtype information
if domain_dtype is None:
self.domain_dtype = None
self.logger.info("Setting domain_dtype to np.float.")
self.domain_dtype = np.float
else:
self.domain_dtype = np.dtype(domain_dtype)
if target_dtype is None:
self.target_dtype = None
self.logger.info("Setting target_dtype to np.complex.")
self.target_dtype = np.complex
else:
self.target_dtype = np.dtype(target_dtype)
......
......@@ -21,7 +21,7 @@ import numpy as np
from nifty.config import dependency_injector as gdi
from nifty import GLSpace, LMSpace
from slicing_transformation import SlicingTransformation
import lm_transformation_factory
import lm_transformation_helper
pyHealpix = gdi.get('pyHealpix')
......@@ -31,12 +31,17 @@ class GLLMTransformation(SlicingTransformation):
# ---Overwritten properties and methods---
def __init__(self, domain, codomain=None, module=None):
if module is None:
module = 'pyHealpix'
if module != 'pyHealpix':
raise ValueError("Unsupported SHT module.")
if 'pyHealpix' not in gdi:
raise ImportError(
"The module pyHealpix is needed but not available.")
super(GLLMTransformation, self).__init__(domain, codomain,
module=module)
super(GLLMTransformation, self).__init__(domain, codomain, module)
# ---Mandatory properties and methods---
......@@ -63,7 +68,7 @@ class GLLMTransformation(SlicingTransformation):
nlat = domain.nlat
lmax = nlat - 1
result = LMSpace(lmax=lmax, dtype=domain.dtype)
result = LMSpace(lmax=lmax)
return result
@classmethod
......@@ -91,6 +96,11 @@ class GLLMTransformation(SlicingTransformation):
super(GLLMTransformation, cls).check_codomain(domain, codomain)
def _transformation_of_slice(self, inp, **kwargs):
if inp.dtype not in (np.float, np.complex):
self.logger.warn("The input array has dtype: %s. The FFT will "
"be performed at double precision." %
str(inp.dtype))
nlat = self.domain.nlat
nlon = self.domain.nlon
lmax = self.codomain.lmax
......@@ -104,12 +114,12 @@ class GLLMTransformation(SlicingTransformation):
for x in (inp.real, inp.imag)]
[resultReal,
resultImag] = [lm_transformation_factory.buildIdx(x, lmax=lmax)
resultImag] = [lm_transformation_helper.buildIdx(x, lmax=lmax)
for x in [resultReal, resultImag]]
result = self._combine_complex_result(resultReal, resultImag)
else:
result = sjob.map2alm(inp)
result = lm_transformation_factory.buildIdx(result, lmax=lmax)
result = lm_transformation_helper.buildIdx(result, lmax=lmax)
return result
......@@ -22,7 +22,7 @@ from nifty.config import dependency_injector as gdi
from nifty import HPSpace, LMSpace
from slicing_transformation import SlicingTransformation
import lm_transformation_factory
import lm_transformation_helper
pyHealpix = gdi.get('pyHealpix')
......@@ -32,12 +32,17 @@ class HPLMTransformation(SlicingTransformation):
# ---Overwritten properties and methods---
def __init__(self, domain, codomain=None, module=None):
if module is None:
module = 'pyHealpix'
if module != 'pyHealpix':
raise ValueError("Unsupported SHT module.")
if 'pyHealpix' not in gdi:
raise ImportError(
"The module pyHealpix is needed but not available")
super(HPLMTransformation, self).__init__(domain, codomain,
module=module)
super(HPLMTransformation, self).__init__(domain, codomain, module)
# ---Mandatory properties and methods---
......@@ -63,7 +68,7 @@ class HPLMTransformation(SlicingTransformation):
lmax = 2*domain.nside
result = LMSpace(lmax=lmax, dtype=domain.dtype)
result = LMSpace(lmax=lmax)
return result
@classmethod
......@@ -83,6 +88,11 @@ class HPLMTransformation(SlicingTransformation):
super(HPLMTransformation, cls).check_codomain(domain, codomain)
def _transformation_of_slice(self, inp, **kwargs):
if inp.dtype not in (np.float, np.complex):
self.logger.warn("The input array has dtype: %s. The FFT will "
"be performed at double precision." %
str(inp.dtype))
lmax = self.codomain.lmax
mmax = lmax
......@@ -92,13 +102,13 @@ class HPLMTransformation(SlicingTransformation):
for x in (inp.real, inp.imag)]
[resultReal,
resultImag] = [lm_transformation_factory.buildIdx(x, lmax=lmax)
resultImag] = [lm_transformation_helper.buildIdx(x, lmax=lmax)
for x in [resultReal, resultImag]]
result = self._combine_complex_result(resultReal, resultImag)
else:
result = pyHealpix.map2alm_iter(inp, lmax, mmax, 3)
result = lm_transformation_factory.buildIdx(result, lmax=lmax)
result = lm_transformation_helper.buildIdx(result, lmax=lmax)
return result
......@@ -21,7 +21,7 @@ from nifty.config import dependency_injector as gdi
from nifty import GLSpace, LMSpace
from slicing_transformation import SlicingTransformation
import lm_transformation_factory
import lm_transformation_helper
pyHealpix = gdi.get('pyHealpix')
......@@ -31,12 +31,17 @@ class LMGLTransformation(SlicingTransformation):
# ---Overwritten properties and methods---
def __init__(self, domain, codomain=None, module=None):
if module is None:
module = 'pyHealpix'
if module != 'pyHealpix':
raise ValueError("Unsupported SHT module.")
if 'pyHealpix' not in gdi:
raise ImportError(
"The module pyHealpix is needed but not available.")
super(LMGLTransformation, self).__init__(domain, codomain,
module=module)
super(LMGLTransformation, self).__init__(domain, codomain, module)
# ---Mandatory properties and methods---
......@@ -69,7 +74,7 @@ class LMGLTransformation(SlicingTransformation):
nlat = domain.lmax + 1
nlon = domain.lmax*2 + 1
result = GLSpace(nlat=nlat, nlon=nlon, dtype=domain.dtype)
result = GLSpace(nlat=nlat, nlon=nlon)
return result
@classmethod
......@@ -97,6 +102,11 @@ class LMGLTransformation(SlicingTransformation):
super(LMGLTransformation, cls).check_codomain(domain, codomain)
def _transformation_of_slice(self, inp, **kwargs):
if inp.dtype not in (np.float, np.complex):
self.logger.warn("The input array has dtype: %s. The FFT will "
"be performed at double precision." %
str(inp.dtype))
nlat = self.codomain.nlat
nlon = self.codomain.nlon
lmax = self.domain.lmax
......@@ -107,7 +117,7 @@ class LMGLTransformation(SlicingTransformation):
sjob.set_triangular_alm_info(lmax, mmax)
if issubclass(inp.dtype.type, np.complexfloating):
[resultReal,
resultImag] = [lm_transformation_factory.buildLm(x, lmax=lmax)
resultImag] = [lm_transformation_helper.buildLm(x, lmax=lmax)
for x in (inp.real, inp.imag)]
[resultReal, resultImag] = [sjob.alm2map(x)
......@@ -116,7 +126,7 @@ class LMGLTransformation(SlicingTransformation):
result = self._combine_complex_result(resultReal, resultImag)
else:
result = lm_transformation_factory.buildLm(inp, lmax=lmax)
result = lm_transformation_helper.buildLm(inp, lmax=lmax)
result = sjob.alm2map(result)
return result
......@@ -20,7 +20,7 @@ import numpy as np
from nifty.config import dependency_injector as gdi
from nifty import HPSpace, LMSpace
from slicing_transformation import SlicingTransformation
import lm_transformation_factory
import lm_transformation_helper
pyHealpix = gdi.get('pyHealpix')
......@@ -30,12 +30,17 @@ class LMHPTransformation(SlicingTransformation):
# ---Overwritten properties and methods---
def __init__(self, domain, codomain=None, module=None):
if module is None:
module = 'pyHealpix'
if module != 'pyHealpix':
raise ValueError("Unsupported SHT module.")
if gdi.get('pyHealpix') is None:
raise ImportError(
"The module pyHealpix is needed but not available.")
super(LMHPTransformation, self).__init__(domain, codomain,
module=module)
super(LMHPTransformation, self).__init__(domain, codomain, module)
# ---Mandatory properties and methods---
......@@ -65,7 +70,7 @@ class LMHPTransformation(SlicingTransformation):
raise TypeError("domain needs to be a LMSpace.")
nside = max((domain.lmax + 1)//2, 1)
result = HPSpace(nside=nside, dtype=domain.dtype)
result = HPSpace(nside=nside)
return result
@classmethod
......@@ -85,13 +90,18 @@ class LMHPTransformation(SlicingTransformation):
super(LMHPTransformation, cls).check_codomain(domain, codomain)
def _transformation_of_slice(self, inp, **kwargs):