still one failing test

nifty/spaces/power_space/power_indices.py

-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-# Copyright(C) 2013-2017 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
-# and financially supported by the Studienstiftung des deutschen Volkes.
-
-import numpy as np
-from d2o import distributed_data_object
-
-class PowerIndices(object):
-    """Computes helpful quantities to deal with power spectra.
-
-    Given the shape and the density of a underlying grid this class provides
-    the user with the pindex, kindex and rho. The indices are binned
-    according to the supplied parameter scheme.
-    """
-    @staticmethod
-    def get_arrays(domain, distribution_strategy, logarithmic, nbin, binbounds):
-        """
-            Returns a dictionary containing the pindex, kindex and rho
-            binned according to the supplied parameter scheme and a
-            configuration dict containing this scheme.
-
-            Parameters
-            ----------
-            domain : NIFTy harmonic space
-                The space for which the power indices get computed
-            distribution_strategy : str
-                The distribution_strategy that will be used for the k_array and
-                pindex distributed_data_object.
-            logarithmic : bool
-                Flag specifying if the binning is performed on logarithmic
-                scale.
-            nbin : integer
-                Number of used bins.
-            binbounds : {list, array}
-                Array-like inner boundaries of the used bins.
-
-        """
-        # if no binning is requested, compute the indices, build the dict,
-        # and return it straight.
-        if not logarithmic and nbin is None and binbounds is None:
-            k_array = domain.get_distance_array(distribution_strategy)
-            temp_kindex = k_array.unique()
-            # compute the local pindex slice on basis of the local k_array data
-            local_pindex = np.searchsorted(temp_kindex, k_array.get_local_data())
-            # prepare the distributed_data_object
-            temp_pindex = distributed_data_object(
-                            global_shape=k_array.shape,
-                            dtype=local_pindex.dtype,
-                            distribution_strategy=distribution_strategy)
-            # store the local pindex data in the global_pindex d2o
-            temp_pindex.set_local_data(local_pindex)
-            temp_rho = temp_pindex.bincount().get_full_data()
-            temp_k_array = k_array
-
-        # if binning is required, make a recursive call to get the unbinned
-        # indices, bin them, and then return everything.
-        else:
-            # Get the unbinned indices
-            pindex, kindex, rho, dummy = PowerIndices.get_arrays(domain,distribution_strategy,nbin=None,
-                                                        binbounds=None,
-                                                        logarithmic=False)
-            # Bin them
-            (temp_pindex, temp_kindex, temp_rho) = \
-                PowerIndices._bin_power_indices(
-                    pindex, kindex, rho, logarithmic, nbin, binbounds)
-            # Make a binned version of k_array
-            tempindex = temp_pindex.copy(dtype=temp_kindex.dtype)
-            temp_k_array = tempindex.apply_scalar_function(
-                            lambda x: temp_kindex[x.astype(np.dtype('int'))])
-
-        return temp_pindex, temp_kindex, temp_rho, temp_k_array
-
-    @staticmethod
-    def _bin_power_indices(pindex, kindex, rho, logarithmic, nbin, binbounds):
-        """
-            Returns the binned power indices associated with the Fourier grid.
-
-            Parameters
-            ----------
-            pindex : distributed_data_object
-                Index of the Fourier grid points in a distributed_data_object.
-            kindex : ndarray
-                Array of all k-vector lengths.
-            rho : ndarray
-                Degeneracy factor of the individual k-vectors.
-            logarithmic : bool
-                Flag specifying if the binning is performed on logarithmic
-                scale.
-            nbin : integer
-                Number of used bins.
-            binbounds : {list, array}
-                Array-like inner boundaries of the used bins.
-
-            Returns
-            -------
-            pindex : distributed_data_object
-            kindex, rho : ndarrays
-                The (re)binned power indices.
-
-        """
-
-        # boundaries
-        if(binbounds is not None):
-            binbounds = np.sort(binbounds)
-        # equal binning
-        else:
-            if(logarithmic is None):
-                logarithmic = False
-            if(logarithmic):
-                k = np.r_[0, np.log(kindex[1:])]
-            else:
-                k = kindex
-            dk = np.max(k[2:] - k[1:-1])  # maximal dk
-            if(nbin is None):
-                nbin = int((k[-1] - 0.5 * (k[2] + k[1])) /
-                           dk - 0.5)  # maximal nbin
-            else:
-                nbin = min(int(nbin), int(
-                    (k[-1] - 0.5 * (k[2] + k[1])) / dk + 2.5))
-                dk = (k[-1] - 0.5 * (k[2] + k[1])) / (nbin - 2.5)
-            binbounds = np.r_[0.5 * (3 * k[1] - k[2]),
-                              0.5 * (k[1] + k[2]) + dk * np.arange(nbin - 2)]
-            if(logarithmic):
-                binbounds = np.exp(binbounds)
-        # reordering
-        reorder = np.searchsorted(binbounds, kindex)
-        rho_ = np.bincount(reorder,weights=rho).astype(rho.dtype)
-        kindex_ = np.bincount(reorder,weights=kindex*rho)/rho_
-        pindex_ = pindex.copy_empty()
-        pindex_.set_local_data(reorder[pindex.get_local_data()])
-
-        return pindex_, kindex_, rho_

nifty/spaces/power_space/power_space.py

@@ -18,10 +18,9 @@
 
 import numpy as np
 
-import d2o
+from d2o import distributed_data_object
 
 from nifty.spaces.space import Space
-from power_indices import PowerIndices
 
 
 class PowerSpace(Space):
@@ -42,8 +41,8 @@ class PowerSpace(Space):
         (default : None).
         if nbin == None, then nbin is set to the length of kindex.
     binbounds :  {list, array-like} *optional*
-        Array-like inner boundaries of the used bins of the default
-        indices.
+        Boundaries between the power spectrum bins.
+        (If binbounds has n entries, there will be n+1 bins.)
         (default : None)
         if binbounds == None :
             Calculates the bounds from the kindex while applying the
@@ -67,19 +66,8 @@ class PowerSpace(Space):
         The total volume of the space.
     shape : tuple of np.ints
         The shape of the space's data array.
-    logarithmic : bool
-        True if logarithmic binning should be used.
-    nbin : {int, None}
-        The number of bins that should be used for power spectrum binning
-        (default : None).
-        if nbin == None, then nbin is set to the length of kindex.
     binbounds :  {list, array-like}
-        Array-like inner boundaries of the used bins of the default
-        indices.
-        (default : None)
-        if binbounds == None :
-            Calculates the bounds from the kindex while applying the
-            logarithmic and nbin keywords.
+        Boundaries between the power spectrum bins
 
     Notes
     -----
@@ -96,34 +84,15 @@ class PowerSpace(Space):
         super(PowerSpace, self).__init__()
         self._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_k_array']
 
-        if not isinstance(harmonic_partner, Space):
-            raise ValueError(
-                "harmonic_partner must be a Space.")
-        if not harmonic_partner.harmonic:
-            raise ValueError(
-                "harmonic_partner must be a harmonic space.")
+        if not (isinstance(harmonic_partner, Space) and \
+                harmonic_partner.harmonic):
+            raise ValueError("harmonic_partner must be a harmonic space.")
         self._harmonic_partner = harmonic_partner
-        try:
-            self._logarithmic = bool(logarithmic)
-        except(TypeError):
-            self._logarithmic = False
-        try:
-            self._nbin = int(nbin)
-        except(TypeError):
-            self._nbin = None
-        try:
-            self._binbounds = tuple(np.array(binbounds))
-        except(TypeError):
-            self._binbounds = None
 
-#        tmp = PowerIndices(self.harmonic_partner, distribution_strategy)
-        self._pindex, self._kindex, self._rho, self._k_array = PowerIndices.get_arrays(self.harmonic_partner, distribution_strategy, self._logarithmic,
-                                                   self._nbin, self._binbounds)
+        self._k_array = self._harmonic_partner.get_distance_array(distribution_strategy)
+        self._compute_binbounds (binbounds, logarithmic, nbin)
 
-        if self._nbin is not None:
-            if self._nbin > len(self.kindex):
-                self.logger.warn("nbin was set to a value being larger than "
-                                 "the length of kindex!")
+        self._do_binning()
 
     def pre_cast(self, x, axes):
         """ Casts power spectrum functions to discretized power spectra.
@@ -154,9 +123,9 @@ class PowerSpace(Space):
 
     def __repr__(self):
         return ("PowerSpace(harmonic_partner=%r, distribution_strategy=%r, "
-                "logarithmic=%r, nbin=%r, binbounds=%r)"
+                "binbounds=%r)"
                 % (self.harmonic_partner, self.pindex.distribution_strategy,
-                   self._logarithmic, self._nbin, self._binbounds))
+                   self._binbounds))
 
     @property
     def harmonic(self):
@@ -179,8 +148,6 @@ class PowerSpace(Space):
         distribution_strategy = self.pindex.distribution_strategy
         return self.__class__(harmonic_partner=self.harmonic_partner,
                               distribution_strategy=distribution_strategy,
-                              logarithmic=self._logarithmic,
-                              nbin=self._nbin,
                               binbounds=self._binbounds)
 
     def weight(self, x, power=1, axes=None, inplace=False):
@@ -201,10 +168,9 @@ class PowerSpace(Space):
         return result_x
 
     def get_distance_array(self, distribution_strategy):
-        result = d2o.distributed_data_object(
+        return distributed_data_object(
                                 self.kindex, dtype=np.float64,
                                 distribution_strategy=distribution_strategy)
-        return result
 
     def get_fft_smoothing_kernel_function(self, sigma):
         raise NotImplementedError(
@@ -218,14 +184,6 @@ class PowerSpace(Space):
         """
         return self._harmonic_partner
 
-    @property
-    def logarithmic(self):
-        return self._logarithmic
-
-    @property
-    def nbin(self):
-        return self._nbin
-
     @property
     def binbounds(self):
         return self._binbounds
@@ -261,9 +219,6 @@ class PowerSpace(Space):
     def _to_hdf5(self, hdf5_group):
         hdf5_group['kindex'] = self.kindex
         hdf5_group['rho'] = self.rho
-        hdf5_group['logarithmic'] = self._logarithmic
-        # Store nbin as string, since it can be None
-        hdf5_group.attrs['nbin'] = str(self._nbin)
         hdf5_group.attrs['binbounds'] = str(self._binbounds)
 
         #MR FIXME: why not "return None" as happens everywhere else?
@@ -285,8 +240,6 @@ class PowerSpace(Space):
         new_ps._harmonic_partner = repository.get('harmonic_partner',
                                                   hdf5_group)
 
-        new_ps._logarithmic = hdf5_group['logarithmic'][()]
-        exec("new_ps._nbin = " + hdf5_group.attrs['nbin'])
         exec("new_ps._binbounds = " + hdf5_group.attrs['binbounds'])
 
         new_ps._pindex = repository.get('pindex', hdf5_group)
@@ -297,6 +250,74 @@ class PowerSpace(Space):
 
         return new_ps
 
+    def _do_binning(self):
+        """ Computes pindex, kindex and rho according to self._binbounds.
+        """
+        # if no binning is requested, compute the "natural" binning, i.e. there
+        # is one bin for every distinct k-vector length
+        if self._binbounds is None:
+            self._kindex = self._k_array.unique()
+            # prepare the distributed_data_object
+            self._pindex = distributed_data_object(
+                global_shape=self._k_array.shape,
+                dtype=np.int,
+                distribution_strategy=self._k_array.distribution_strategy)
+            # store the local pindex data in the global_pindex d2o
+            self._pindex.set_local_data(
+                np.searchsorted(self._kindex, self._k_array.get_local_data()))
+            self._rho = self._pindex.bincount().get_full_data()
+
+        # use the provided binbounds
+        else:
+            self._pindex = distributed_data_object(
+                global_shape=self._k_array.shape,
+                dtype=np.int,
+                distribution_strategy=self._k_array.distribution_strategy)
+            self._pindex.set_local_data(np.searchsorted(
+                self._binbounds, self._k_array.get_local_data()))
+            self._rho = self._pindex.bincount().get_full_data()
+            self._kindex = self._pindex.bincount(
+                weights=self._k_array).get_full_data()/self._rho
+
+    def _compute_binbounds (self, binbounds, logarithmic, nbin):
+        try:
+            self._binbounds = np.sort(tuple(np.array(binbounds)))
+        except(TypeError):
+            self._binbounds = None
+
+        if(self._binbounds is None):
+            try:
+                logarithmic = bool(logarithmic)
+            except(TypeError):
+                logarithmic = False
+            try:
+                nbin = int(nbin)
+            except(TypeError):
+                nbin = None
+
+            if not logarithmic and nbin is None:
+                return # no binbounds, use "natural" binning
+
+            # equal binning
+            # MR FIXME: this needs to improve
+            kindex = self._k_array.unique()
+            if (logarithmic):
+                k = np.r_[0, np.log(kindex[1:])]
+            else:
+                k = kindex
+            dk = np.max(k[2:] - k[1:-1])  # minimum dk to avoid empty bins
+            if(nbin is None):
+                nbin = int((k[-1] - 0.5 * (k[2] + k[1])) /
+                           dk - 0.5)  # maximal nbin
+            else:
+                nbin = min(int(nbin), int(
+                    (k[-1] - 0.5 * (k[2] + k[1])) / dk + 2.5))
+                dk = (k[-1] - 0.5 * (k[2] + k[1])) / (nbin - 2.5)
+            self._binbounds = np.r_[0.5 * (3 * k[1] - k[2]),
+                              0.5 * (k[1] + k[2]) + dk * np.arange(nbin - 2)]
+            if(logarithmic):
+                self._binbounds = np.exp(self._binbounds)
+
 
 class EmptyPowerSpace(PowerSpace):
     def __init__(self):

test/test_spaces/test_power_space.py

@@ -64,8 +64,6 @@ CONSTRUCTOR_CONFIGS = [
         'dim': 5,
         'total_volume': 8.0,
         'harmonic_partner': RGSpace((8,), harmonic=True),
-        'logarithmic': False,
-        'nbin': None,
         'binbounds': None,
         'pindex': distributed_data_object([0, 1, 2, 3, 4, 3, 2, 1]),
         'kindex': np.array([0., 1., 2., 3., 4.]),
@@ -78,8 +76,6 @@ CONSTRUCTOR_CONFIGS = [
         'dim': 2,
         'total_volume': 8.0,
         'harmonic_partner': RGSpace((8,), harmonic=True),
-        'logarithmic': True,
-        'nbin': None,
         'binbounds': None,
         'pindex': distributed_data_object([0, 1, 1, 1, 1, 1, 1, 1]),
         'kindex': np.array([0., 2.28571429]),
@@ -118,8 +114,6 @@ def get_weight_configs():
 class PowerSpaceInterfaceTest(unittest.TestCase):
     @expand([
         ['harmonic_partner', Space],
-        ['logarithmic', bool],
-        ['nbin', NoneType],
         ['binbounds', NoneType],
         ['pindex', distributed_data_object],
         ['kindex', np.ndarray],
@@ -143,7 +137,7 @@ class PowerSpaceConsistencyCheck(unittest.TestCase):
             err_msg='rho is not equal to pindex degeneracy')
 
 class PowerSpaceFunctionalityTest(unittest.TestCase):
-    @expand(CONSISTENCY_CONFIGS)
+    @expand(CONSTRUCTOR_CONFIGS)
     def test_constructor(self, harmonic_partner, distribution_strategy,
                          logarithmic, nbin, binbounds, expected):
         if 'error' in expected: