merge master

demos/wiener_filter.py

@@ -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')
-#

demos/wiener_filter_hamiltonian.py

@@ -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)

nifty/domain_object.py

@@ -18,8 +18,6 @@
 
 import abc
 
-import numpy as np
-
 from keepers import Loggable,\
                     Versionable
 
@@ -27,9 +25,9 @@ from keepers import Loggable,\
 class DomainObject(Versionable, Loggable, object):
     __metaclass__ = abc.ABCMeta
 
-    def __init__(self, dtype):
-        self._dtype = np.dtype(dtype)
-        self._ignore_for_hash = []
+    def __init__(self):
+        # _global_id is used in the Versioning module from keepers
+        self._ignore_for_hash = ['_global_id']
 
     def __hash__(self):
         # Extract the identifying parts from the vars(self) dict.
@@ -43,17 +41,20 @@ class DomainObject(Versionable, Loggable, object):
 
     def __eq__(self, x):
         if isinstance(x, type(self)):
-            return hash(self) == hash(x)
+            for key in vars(self).keys():
+                item1 = vars(self)[key]
+                if key in self._ignore_for_hash or key == '_ignore_for_hash':
+                    continue
+                item2 = vars(x)[key]
+                if item1 != item2:
+                    return False
+            return True
         else:
             return False
 
     def __ne__(self, x):
         return not self.__eq__(x)
 
-    @property
-    def dtype(self):
-        return self._dtype
-
     @abc.abstractproperty
     def shape(self):
         raise NotImplementedError(
@@ -78,10 +79,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

nifty/field.py

@@ -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
 
@@ -210,16 +210,10 @@ class Field(Loggable, Versionable, object):
             self.val.get_axes_local_distribution_strategy(
                 self.domain_axes[space_index])
 
-        if real_signal:
-            power_dtype = np.dtype('complex')
-        else:
-            power_dtype = np.dtype('float')
-
         harmonic_domain = self.domain[space_index]
         power_domain = PowerSpace(harmonic_domain=harmonic_domain,
                                   distribution_strategy=distribution_strategy,
-                                  log=log, nbin=nbin, binbounds=binbounds,
-                                  dtype=power_dtype)
+                                  log=log, nbin=nbin, binbounds=binbounds)
 
         # extract pindex and rho from power_domain
         pindex = power_domain.pindex
@@ -251,8 +245,14 @@ class Field(Loggable, Versionable, object):
         result_domain = list(self.domain)
         result_domain[space_index] = power_domain
 
+        if real_signal:
+            result_dtype = np.complex
+        else:
+            result_dtype = np.float
+
         result_field = self.copy_empty(
                    domain=result_domain,
+                   dtype=result_dtype,
                    distribution_strategy=power_spectrum.distribution_strategy)
         result_field.set_val(new_val=power_spectrum, copy=False)
 
@@ -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---
 

nifty/field_types/field_array.py

@@ -16,20 +16,18 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-import numpy as np
-
 from field_type import FieldType
 
 
 class FieldArray(FieldType):
 
-    def __init__(self, shape, dtype=np.float):
+    def __init__(self, shape):
         try:
             new_shape = tuple([int(i) for i in shape])
         except TypeError:
             new_shape = (int(shape), )
         self._shape = new_shape
-        super(FieldArray, self).__init__(dtype=dtype)
+        super(FieldArray, self).__init__()
 
     @property
     def shape(self):
@@ -43,14 +41,9 @@ class FieldArray(FieldType):
 
     def _to_hdf5(self, hdf5_group):
         hdf5_group['shape'] = self.shape
-        hdf5_group.attrs['dtype'] = self.dtype.name
-
         return None
 
     @classmethod
     def _from_hdf5(cls, hdf5_group, loopback_get):
-        result = cls(
-            shape=hdf5_group['shape'][:],
-            dtype=np.dtype(hdf5_group.attrs['dtype'])
-            )
+        result = cls(shape=hdf5_group['shape'][:])
         return result

nifty/operators/diagonal_operator/diagonal_operator.py

@@ -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

nifty/operators/endomorphic_operator/endomorphic_operator.py

@@ -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

nifty/operators/fft_operator/fft_operator.py

@@ -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)
 

nifty/operators/fft_operator/transformations/gllmtransformation.py

@@ -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

nifty/operators/fft_operator/transformations/hplmtransformation.py

@@ -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

nifty/operators/fft_operator/transformations/lmgltransformation.py

@@ -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.")
+