From e47faf2db7a2e3c9b7d75b2efa0f97f7a6399004 Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 7 Jul 2017 16:04:29 +0200
Subject: [PATCH] Revert "Merge branch 'index_games' into working_on_demos"

This reverts commit df5326d210ad33d34a46b106c194241e4dcdd509, reversing
changes made to 91f06b354ba2748c67f0c1566833c3288b374aa9.
---
 nifty/field.py                                |  18 +-
 .../response_operator/response_operator.py    |   3 +-
 nifty/spaces/lm_space/lm_space.py             |  25 +-
 nifty/spaces/power_space/__init__.py          |   1 +
 .../spaces/power_space/power_index_factory.py |  43 ++
 nifty/spaces/power_space/power_indices.py     | 416 ++++++++++++++++++
 nifty/spaces/power_space/power_space.py       | 208 +++++----
 nifty/spaces/rg_space/rg_space.py             |  30 --
 nifty/spaces/space/space.py                   |  14 -
 test/test_spaces/test_power_space.py          |  41 +-
 10 files changed, 619 insertions(+), 180 deletions(-)
 create mode 100644 nifty/spaces/power_space/power_index_factory.py
 create mode 100644 nifty/spaces/power_space/power_indices.py

diff --git a/nifty/field.py b/nifty/field.py
index 533eefc61..fc697139f 100644
--- a/nifty/field.py
+++ b/nifty/field.py
@@ -294,8 +294,7 @@ class Field(Loggable, Versionable, object):
             if nbin==None : maximum number of bins is used
         binbounds : array-like *optional*
             Inner bounds of the bins (default : None).
-            Overrides nbin and logarithmic.
-            if binbounds==None : bins are inferred.
+            if binbounds==None : bins are inferred. Overwrites nbins and log
         keep_phase_information : boolean, *optional*
             If False, return a real-valued result containing the power spectrum
             of the input Field.
@@ -398,9 +397,14 @@ class Field(Loggable, Versionable, object):
                                   logarithmic=logarithmic, nbin=nbin,
                                   binbounds=binbounds)
 
+        # extract pindex and rho from power_domain
+        pindex = power_domain.pindex
+        rho = power_domain.rho
+
         power_spectrum = cls._calculate_power_spectrum(
                                 field_val=work_field.val,
-                                pdomain=power_domain,
+                                pindex=pindex,
+                                rho=rho,
                                 axes=work_field.domain_axes[space_index])
 
         # create the result field and put power_spectrum into it
@@ -417,11 +421,8 @@ class Field(Loggable, Versionable, object):
         return result_field
 
     @classmethod
-    def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
+    def _calculate_power_spectrum(cls, field_val, pindex, rho, axes=None):
 
-        pindex = pdomain.pindex
-        # MR FIXME: how about iterating over slices, instead of replicating
-        # pindex? Would save memory and probably isn't slower.
         if axes is not None:
             pindex = cls._shape_up_pindex(
                             pindex=pindex,
@@ -430,7 +431,6 @@ class Field(Loggable, Versionable, object):
                             axes=axes)
         power_spectrum = pindex.bincount(weights=field_val,
                                          axis=axes)
-        rho = pdomain.rho
         if axes is not None:
             new_rho_shape = [1, ] * len(power_spectrum.shape)
             new_rho_shape[axes[0]] = len(rho)
@@ -762,7 +762,7 @@ class Field(Loggable, Versionable, object):
         Returns
         -------
         out : tuple
-            The output object. The tuple contains the dimensions of the spaces
+            The output object. The tuple contains the dimansions of the spaces
             in domain.
 
         See Also
diff --git a/nifty/operators/response_operator/response_operator.py b/nifty/operators/response_operator/response_operator.py
index cb3172077..75361469d 100644
--- a/nifty/operators/response_operator/response_operator.py
+++ b/nifty/operators/response_operator/response_operator.py
@@ -33,7 +33,8 @@ class ResponseOperator(LinearOperator):
     domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
         The domain on which the Operator's input Field lives.
     target : tuple of DomainObjects, i.e. Spaces and FieldTypes
-        The domain in which the outcome of the operator lives.
+        The domain in which the outcome of the operator lives. As the Operator
+        is endomorphic this is the same as its domain.
     unitary : boolean
         Indicates whether the Operator is unitary or not.
 
diff --git a/nifty/spaces/lm_space/lm_space.py b/nifty/spaces/lm_space/lm_space.py
index d131dea69..c542339fc 100644
--- a/nifty/spaces/lm_space/lm_space.py
+++ b/nifty/spaces/lm_space/lm_space.py
@@ -106,6 +106,9 @@ class LMSpace(Space):
 
         return (hermitian_part, anti_hermitian_part)
 
+#    def hermitian_fixed_points(self):
+#        return None
+
     # ---Mandatory properties and methods---
 
     def __repr__(self):
@@ -144,22 +147,6 @@ class LMSpace(Space):
             return x.copy()
 
     def get_distance_array(self, distribution_strategy):
-        if distribution_strategy == 'not':  # short cut
-            lmax = self.lmax
-            ldist = np.empty((self.dim,), dtype=np.float64)
-            ldist[0:lmax+1] = np.arange(lmax+1, dtype=np.float64)
-            tmp = np.empty((2*lmax+2), dtype=np.float64)
-            tmp[0::2] = np.arange(lmax+1)
-            tmp[1::2] = np.arange(lmax+1)
-            idx = lmax+1
-            for l in range(1, lmax+1):
-                ldist[idx:idx+2*(lmax+1-l)] = tmp[2*l:]
-                idx += 2*(lmax+1-l)
-            dists = arange(start=0, stop=self.shape[0],
-                           distribution_strategy=distribution_strategy)
-            dists.set_local_data(ldist)
-            return dists
-
         dists = arange(start=0, stop=self.shape[0],
                        distribution_strategy=distribution_strategy)
 
@@ -169,12 +156,6 @@ class LMSpace(Space):
 
         return dists
 
-    def get_unique_distances(self):
-        return np.arange(self.lmax+1, dtype=np.float64)
-
-    def get_natural_binbounds(self):
-        return np.arange(self.lmax, dtype=np.float64) + 0.5
-
     @staticmethod
     def _distance_array_helper(index_array, lmax):
         u = 2*lmax + 1
diff --git a/nifty/spaces/power_space/__init__.py b/nifty/spaces/power_space/__init__.py
index 4801bdbe1..cad1d25cf 100644
--- a/nifty/spaces/power_space/__init__.py
+++ b/nifty/spaces/power_space/__init__.py
@@ -17,3 +17,4 @@
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
 from power_space import PowerSpace
+from power_index_factory import PowerIndexFactory
\ No newline at end of file
diff --git a/nifty/spaces/power_space/power_index_factory.py b/nifty/spaces/power_space/power_index_factory.py
new file mode 100644
index 000000000..8caeee21d
--- /dev/null
+++ b/nifty/spaces/power_space/power_index_factory.py
@@ -0,0 +1,43 @@
+# 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.
+
+from power_indices import PowerIndices
+
+
+class _PowerIndexFactory(object):
+    def __init__(self):
+        self.power_indices_storage = {}
+
+    def get_power_index(self, domain, distribution_strategy,
+                        logarithmic=False, nbin=None, binbounds=None):
+        key = (domain, distribution_strategy)
+
+        if key not in self.power_indices_storage:
+            self.power_indices_storage[key] = \
+                PowerIndices(domain, distribution_strategy,
+                             logarithmic=logarithmic,
+                             nbin=nbin,
+                             binbounds=binbounds)
+        power_indices = self.power_indices_storage[key]
+        power_index = power_indices.get_index_dict(logarithmic=logarithmic,
+                                                   nbin=nbin,
+                                                   binbounds=binbounds)
+        return power_index
+
+
+PowerIndexFactory = _PowerIndexFactory()
diff --git a/nifty/spaces/power_space/power_indices.py b/nifty/spaces/power_space/power_indices.py
new file mode 100644
index 000000000..9841a005c
--- /dev/null
+++ b/nifty/spaces/power_space/power_indices.py
@@ -0,0 +1,416 @@
+# 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,\
+                STRATEGIES as DISTRIBUTION_STRATEGIES
+
+from d2o.config import dependency_injector as d2o_di
+from d2o.config import configuration as d2o_config
+
+
+class PowerIndices(object):
+    """Computes helpful quantities to deal with power spectra.
+
+    Given the shape and the density of a underlying rectangular grid this
+    class provides the user
+    with the pindex, kindex, rho and pundex. The indices are binned
+    according to the supplied parameter scheme. If wanted, computed
+    results are stored for future reuse.
+
+    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 *optional*
+        Flag specifying if the binning of the default indices is
+        performed on logarithmic scale.
+    nbin : integer *optional*
+        Number of used bins for the binning of the default indices.
+    binbounds : {list, array}
+        Array-like inner boundaries of the used bins of the default
+        indices.
+    """
+    def __init__(self, domain, distribution_strategy,
+                 logarithmic=False, nbin=None, binbounds=None):
+        self.domain = domain
+        self.distribution_strategy = distribution_strategy
+
+        # Compute the global k_array
+        self.k_array = self.domain.get_distance_array(distribution_strategy)
+        # Initialize the dictionary which stores all individual index-dicts
+        self.global_dict = {}
+        # Set self.default_parameters
+        self.set_default(config_dict={'logarithmic': logarithmic,
+                                      'nbin': nbin,
+                                      'binbounds': binbounds})
+
+    # Redirect the direct calls approaching a power_index instance to the
+    # default_indices dict
+    @property
+    def default_indices(self):
+        return self.get_index_dict(**self.default_parameters)
+
+    def __getitem__(self, x):
+        return self.default_indices.get(x)
+
+    def __contains__(self, x):
+        return self.default_indices.__contains__(x)
+
+    def __iter__(self):
+        return self.default_indices.__iter__()
+
+    def __getattr__(self, x):
+        return self.default_indices.__getattribute__(x)
+
+    def set_default(self, **kwargs):
+        """
+            Sets the index-set which is specified by the parameters as the
+            default for the power_index instance.
+
+            Parameters
+            ----------
+            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
+            -------
+            None
+        """
+        parsed_kwargs = self._cast_config(**kwargs)
+        self.default_parameters = parsed_kwargs
+
+    def _cast_config(self, **kwargs):
+        """
+            internal helper function which casts the various combinations of
+            possible parameters into a properly defaulted dictionary
+        """
+        temp_config_dict = kwargs.get('config_dict', None)
+        if temp_config_dict is not None:
+            return self._cast_config_helper(**temp_config_dict)
+        else:
+            defaults = self.default_parameters
+            temp_logarithmic = kwargs.get("logarithmic",
+                                          defaults['logarithmic'])
+            temp_nbin = kwargs.get("nbin", defaults['nbin'])
+            temp_binbounds = kwargs.get("binbounds", defaults['binbounds'])
+
+            return self._cast_config_helper(logarithmic=temp_logarithmic,
+                                            nbin=temp_nbin,
+                                            binbounds=temp_binbounds)
+
+    def _cast_config_helper(self, logarithmic, nbin, binbounds):
+        """
+            internal helper function which sets the defaults for the
+            _cast_config function
+        """
+
+        try:
+            temp_logarithmic = bool(logarithmic)
+        except(TypeError):
+            temp_logarithmic = False
+
+        try:
+            temp_nbin = int(nbin)
+        except(TypeError):
+            temp_nbin = None
+
+        try:
+            temp_binbounds = tuple(np.array(binbounds))
+        except(TypeError):
+            temp_binbounds = None
+
+        temp_dict = {"logarithmic": temp_logarithmic,
+                     "nbin": temp_nbin,
+                     "binbounds": temp_binbounds}
+        return temp_dict
+
+    def get_index_dict(self, **kwargs):
+        """
+            Returns a dictionary containing the pindex, kindex, rho and pundex
+            binned according to the supplied parameter scheme and a
+            configuration dict containing this scheme.
+
+            Parameters
+            ----------
+            store : bool
+                Flag specifying if  the calculated index dictionary should be
+                stored in the global_dict for future use.
+            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
+            -------
+            index_dict : dict
+                Contains the keys: 'config', 'pindex', 'kindex', 'rho' and
+                'pundex'
+        """
+        # Cast the input arguments
+        temp_config_dict = self._cast_config(**kwargs)
+        # Compute a hashable identifier from the config which will be used
+        # as dict key
+        temp_key = self._freeze_config(temp_config_dict)
+        # Check if the result should be stored for future use.
+        storeQ = kwargs.get("store", True)
+        # Try to find the requested index dict in the global_dict
+        try:
+            return self.global_dict[temp_key]
+        except(KeyError):
+            # If it is not found, calculate it.
+            temp_index_dict = self._compute_index_dict(temp_config_dict)
+            # Store it, if required
+            if storeQ:
+                self.global_dict[temp_key] = temp_index_dict
+                # Important: If the result is stored, return a reference to
+                # the dictionary entry, not anly a plain copy. Otherwise,
+                # set_default breaks!
+                return self.global_dict[temp_key]
+            else:
+                # Return the plain result.
+                return temp_index_dict
+
+    def _freeze_config(self, config_dict):
+        """
+            a helper function which forms a hashable identifying object from
+            a config dictionary which can be used as key of a dict
+        """
+        return frozenset(config_dict.items())
+
+    def _compute_index_dict(self, config_dict):
+        """
+            Internal helper function which takes a config_dict, asks for the
+            pindex/kindex/rho/pundex set, and bins them according to the config
+        """
+        # if no binning is requested, compute the indices, build the dict,
+        # and return it straight.
+        if not config_dict["logarithmic"] and config_dict["nbin"] is None and \
+                config_dict["binbounds"] is None:
+            (temp_pindex, temp_kindex, temp_rho, temp_pundex) =\
+                self._compute_indices(self.k_array)
+            temp_k_array = self.k_array
+
+        # if binning is required, make a recursive call to get the unbinned
+        # indices, bin them, compute the pundex and then return everything.
+        else:
+            # Get the unbinned indices
+            temp_unbinned_indices = self.get_index_dict(nbin=None,
+                                                        binbounds=None,
+                                                        logarithmic=False,
+                                                        store=False)
+            # Bin them
+            (temp_pindex, temp_kindex, temp_rho, temp_pundex) = \
+                self._bin_power_indices(
+                    temp_unbinned_indices, **config_dict)
+            # Make a binned version of k_array
+            temp_k_array = self._compute_k_array_from_pindex_kindex(
+                               temp_pindex, temp_kindex)
+
+        temp_index_dict = {"config": config_dict,
+                           "pindex": temp_pindex,
+                           "kindex": temp_kindex,
+                           "rho": temp_rho,
+                           "pundex": temp_pundex,
+                           "k_array": temp_k_array}
+        return temp_index_dict
+
+    def _compute_k_array_from_pindex_kindex(self, pindex, kindex):
+        tempindex = pindex.copy(dtype=kindex.dtype)
+        result = tempindex.apply_scalar_function(
+                            lambda x: kindex[x.astype(np.dtype('int'))])
+        return result
+
+    def _compute_indices(self, k_array):
+        """
+        Internal helper function which computes pindex, kindex, rho and pundex
+        from a given k_array
+        """
+        ##########
+        # kindex #
+        ##########
+        global_kindex = k_array.unique()
+
+        ##########
+        # pindex #
+        ##########
+        # compute the local pindex slice on basis of the local k_array data
+        local_pindex = np.searchsorted(global_kindex, k_array.get_local_data())
+        # prepare the distributed_data_object
+        global_pindex = distributed_data_object(
+                            global_shape=k_array.shape,
+                            dtype=local_pindex.dtype,
+                            distribution_strategy=self.distribution_strategy)
+        # store the local pindex data in the global_pindex d2o
+        global_pindex.set_local_data(local_pindex)
+
+        #######
+        # rho #
+        #######
+        global_rho = global_pindex.bincount().get_full_data()
+
+        ##########
+        # pundex #
+        ##########
+        global_pundex = self._compute_pundex(global_pindex,
+                                             global_kindex)
+
+        return global_pindex, global_kindex, global_rho, global_pundex
+
+    def _compute_pundex(self, global_pindex, global_kindex):
+        """
+        Internal helper function which computes the pundex array from a
+        pindex and a kindex array. This function is separated from the
+        _compute_indices function as it is needed in _bin_power_indices,
+        too.
+        """
+        if self.distribution_strategy in DISTRIBUTION_STRATEGIES['slicing']:
+            ##########
+            # pundex #
+            ##########
+            # Prepare the local data
+            local_pindex = global_pindex.get_local_data()
+            # Compute the local pundices for the local pindices
+            (temp_uniqued_pindex, local_temp_pundex) = np.unique(
+                                                            local_pindex,
+                                                            return_index=True)
+            # Shift the local pundices by the nodes' local_dim_offset
+            local_temp_pundex += global_pindex.distributor.local_dim_offset
+
+            # Prepare the pundex arrays used for the Allreduce operation
+            # pundex has the same length as the kindex array
+            local_pundex = np.zeros(shape=global_kindex.shape, dtype=np.int)
+            # Set the default value higher than the maximal possible pundex
+            # value so that MPI.MIN can sort out the default
+            local_pundex += np.prod(global_pindex.shape) + 1
+            # Set the default value higher than the length
+            global_pundex = np.empty_like(local_pundex)
+            # Store the individual pundices in the local_pundex array
+            local_pundex[temp_uniqued_pindex] = local_temp_pundex
+            # Extract the MPI module from the global_pindex d2o
+            MPI = d2o_di[d2o_config['mpi_module']]
+            # Use Allreduce to find the first occurences/smallest pundices
+            global_pindex.comm.Allreduce(local_pundex,
+                                         global_pundex,
+                                         op=MPI.MIN)
+            return global_pundex
+
+        elif self.distribution_strategy in DISTRIBUTION_STRATEGIES['not']:
+            ##########
+            # pundex #
+            ##########
+            pundex = np.unique(global_pindex.get_local_data(),
+                               return_index=True)[1]
+            return pundex
+        else:
+            raise NotImplementedError(
+                "_compute_pundex_d2o not implemented for given "
+                "distribution_strategy")
+
+    def _bin_power_indices(self, index_dict, **kwargs):
+        """
+            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, pundex : ndarrays
+                The (re)binned power indices.
+
+        """
+        # Cast the given config
+        temp_config_dict = self._cast_config(**kwargs)
+        logarithmic = temp_config_dict['logarithmic']
+        nbin = temp_config_dict['nbin']
+        binbounds = temp_config_dict['binbounds']
+
+        # Extract the necessary indices from the supplied index dict
+        pindex = index_dict["pindex"]
+        kindex = index_dict["kindex"]
+        rho = index_dict["rho"]
+
+        # 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])  # minimal 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.zeros(len(binbounds) + 1, dtype=rho.dtype)
+        kindex_ = np.empty(len(binbounds) + 1, dtype=kindex.dtype)
+        for ii in range(len(reorder)):
+            if(rho_[reorder[ii]] == 0):
+                kindex_[reorder[ii]] = kindex[ii]
+                rho_[reorder[ii]] += rho[ii]
+            else:
+                kindex_[reorder[ii]] = ((kindex_[reorder[ii]] *
+                                         rho_[reorder[ii]] +
+                                         kindex[ii] * rho[ii]) /
+                                        (rho_[reorder[ii]] + rho[ii]))
+                rho_[reorder[ii]] += rho[ii]
+
+        pindex_ = pindex.copy_empty()
+        pindex_.set_local_data(reorder[pindex.get_local_data()])
+        pundex_ = self._compute_pundex(pindex_, kindex_)
+
+        return pindex_, kindex_, rho_, pundex_
diff --git a/nifty/spaces/power_space/power_space.py b/nifty/spaces/power_space/power_space.py
index 4841e958e..9be622643 100644
--- a/nifty/spaces/power_space/power_space.py
+++ b/nifty/spaces/power_space/power_space.py
@@ -18,11 +18,11 @@
 
 import numpy as np
 
-from d2o import distributed_data_object
+import d2o
 
-from nifty.spaces.space import Space
+from power_index_factory import PowerIndexFactory
 
-_PSCache = {}
+from nifty.spaces.space import Space
 
 
 class PowerSpace(Space):
@@ -37,24 +37,18 @@ class PowerSpace(Space):
         derived from this PowerSpace, e.g. the pindex.
         (default : 'not')
     logarithmic : bool *optional*
-        True if logarithmic binning should be used (default : None).
+        True if logarithmic binning should be used (default : False).
     nbin : {int, None} *optional*
         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} *optional*
-        Boundaries between the power spectrum bins.
-        (If binbounds has n entries, there will be n+1 bins, the first bin
-        starting at -inf, the last bin ending at +inf.)
+        Array-like inner boundaries of the used bins of the default
+        indices.
         (default : None)
-        if binbounds == None:
+        if binbounds == None :
             Calculates the bounds from the kindex while applying the
             logarithmic and nbin keywords.
-    Note: if "bindounds" is not None, both "logarithmic" and "nbin" must be
-        None!
-    Note: if "binbounds", "logarithmic", and "nbin" are all None, then
-        "natural" binning is performed, i.e. there will be one bin for every
-        distinct k-vector length.
 
     Attributes
     ----------
@@ -63,6 +57,9 @@ class PowerSpace(Space):
         mapped to which power bin
     kindex : numpy.ndarray
         Sorted array of all k-modes.
+    pundex : numpy.ndarray
+        Flat index of the first occurence of a k-vector with length==kindex[n]
+        in the k_array.
     rho : numpy.ndarray
         The amount of k-modes that get mapped to one power bin is given by
         rho.
@@ -74,8 +71,9 @@ class PowerSpace(Space):
         The total volume of the space.
     shape : tuple of np.ints
         The shape of the space's data array.
-    binbounds :  tuple or None
-        Boundaries between the power spectrum bins
+    config : {logarithmic, nbin, binbounds}
+        Dictionary storing the values for `logarithmic`, `nbin`, and
+        `binbounds` that were used during initialization.
 
     Notes
     -----
@@ -88,79 +86,38 @@ class PowerSpace(Space):
 
     def __init__(self, harmonic_partner,
                  distribution_strategy='not',
-                 logarithmic=None, nbin=None, binbounds=None):
+                 logarithmic=False, nbin=None, binbounds=None):
         super(PowerSpace, self).__init__()
-        self._ignore_for_hash += ['_pindex', '_kindex', '_rho']
+        self._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',
+                                  '_k_array']
 
-        if not (isinstance(harmonic_partner, Space) and
-                harmonic_partner.harmonic):
-            raise ValueError("harmonic_partner must be a harmonic space.")
+        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.")
         self._harmonic_partner = harmonic_partner
 
-        # sanity check
-        if binbounds is not None and not(nbin is None and logarithmic is None):
-            raise ValueError(
-                "if binbounds is defined, nbin and logarithmic must be None")
+        power_index = PowerIndexFactory.get_power_index(
+                        domain=self.harmonic_partner,
+                        distribution_strategy=distribution_strategy,
+                        logarithmic=logarithmic,
+                        nbin=nbin,
+                        binbounds=binbounds)
 
-        if binbounds is not None:
-            binbounds = tuple(binbounds)
+        self._config = power_index['config']
 
-        key = (harmonic_partner, distribution_strategy, logarithmic, nbin,
-               binbounds)
-        if _PSCache.get(key) is not None:
-            (self._binbounds, self._pindex, self._kindex, self._rho) \
-                = _PSCache[key]
-            return
+        self._pindex = power_index['pindex']
+        self._kindex = power_index['kindex']
+        self._rho = power_index['rho']
+        self._pundex = power_index['pundex']
+        self._k_array = power_index['k_array']
 
-        self._binbounds = None
-        if logarithmic is None and nbin is None and binbounds is None:
-            bb = self._harmonic_partner.get_natural_binbounds()
-        else:
-            if binbounds is not None:
-                bb = np.sort(np.array(binbounds))
-            else:
-                if logarithmic is not None:
-                    logarithmic = bool(logarithmic)
-                if nbin is not None:
-                    nbin = int(nbin)
-
-                # equidistant binning (linear or log)
-                # MR FIXME: this needs to improve
-                kindex = self._harmonic_partner.get_unique_distances()
-                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)
-                bb = np.r_[0.5 * (3 * k[1] - k[2]),
-                           0.5 * (k[1] + k[2]) + dk * np.arange(nbin-2)]
-                if(logarithmic):
-                    bb = np.exp(bb)
-            self._binbounds = tuple(bb)
-
-        dists = self._harmonic_partner.get_distance_array(
-                distribution_strategy)
-        # Compute pindex, kindex and rho according to bb
-        self._pindex = distributed_data_object(
-            global_shape=dists.shape,
-            dtype=np.int,
-            distribution_strategy=distribution_strategy)
-
-        self._pindex.set_local_data(np.searchsorted(
-            bb, dists.get_local_data()))  # also expensive!
-        self._rho = self._pindex.bincount().get_full_data()
-        self._kindex = self._pindex.bincount(
-            weights=dists).get_full_data()/self._rho
-
-        _PSCache[key] = \
-            (self._binbounds, self._pindex, self._kindex, self._rho)
+        if self.config['nbin'] is not None:
+            if self.config['nbin'] > len(self.kindex):
+                self.logger.warn("nbin was set to a value being larger than "
+                                 "the length of kindex!")
 
     def pre_cast(self, x, axes):
         """ Casts power spectrum functions to discretized power spectra.
@@ -185,15 +142,19 @@ class PowerSpace(Space):
 
         """
 
-        return x(self.kindex) if callable(x) else x
+        if callable(x):
+            return x(self.kindex)
+        else:
+            return x
 
     # ---Mandatory properties and methods---
 
     def __repr__(self):
         return ("PowerSpace(harmonic_partner=%r, distribution_strategy=%r, "
-                "binbounds=%r)"
+                "logarithmic=%r, nbin=%r, binbounds=%r)"
                 % (self.harmonic_partner, self.pindex.distribution_strategy,
-                   self._binbounds))
+                   self.config['logarithmic'], self.config['nbin'],
+                   self.config['binbounds']))
 
     @property
     def harmonic(self):
@@ -216,9 +177,11 @@ class PowerSpace(Space):
         distribution_strategy = self.pindex.distribution_strategy
         return self.__class__(harmonic_partner=self.harmonic_partner,
                               distribution_strategy=distribution_strategy,
-                              binbounds=self._binbounds)
+                              logarithmic=self.config["logarithmic"],
+                              nbin=self.config["nbin"],
+                              binbounds=self.config["binbounds"])
 
-    def weight(self, x, power, axes, inplace=False):
+    def weight(self, x, power=1, axes=None, inplace=False):
         reshaper = [1, ] * len(x.shape)
         # we know len(axes) is always 1
         reshaper[axes[0]] = self.shape[0]
@@ -236,9 +199,10 @@ class PowerSpace(Space):
         return result_x
 
     def get_distance_array(self, distribution_strategy):
-        return distributed_data_object(
+        result = d2o.distributed_data_object(
                                 self.kindex, dtype=np.float64,
                                 distribution_strategy=distribution_strategy)
+        return result
 
     def get_fft_smoothing_kernel_function(self, sigma):
         raise NotImplementedError(
@@ -253,8 +217,11 @@ class PowerSpace(Space):
         return self._harmonic_partner
 
     @property
-    def binbounds(self):
-        return self._binbounds
+    def config(self):
+        """ Returns the configuration which was used for `logarithmic`, `nbin`
+        and `binbounds` during initialization.
+        """
+        return self._config
 
     @property
     def pindex(self):
@@ -275,20 +242,67 @@ class PowerSpace(Space):
         """
         return self._rho
 
+    @property
+    def pundex(self):
+        """ An array for which the n-th entry gives the flat index of the
+        first occurence of a k-vector with length==kindex[n] in the
+        k_array.
+        """
+        return self._pundex
+
+    @property
+    def k_array(self):
+        """ An array containing distances to the grid center (i.e. zero-mode)
+        for every k-mode in the grid of the harmonic partner space.
+        """
+        return self._k_array
+
     # ---Serialization---
 
     def _to_hdf5(self, hdf5_group):
-        hdf5_group.attrs['binbounds'] = str(self._binbounds)
-        hdf5_group.attrs['distribution_strategy'] = \
-            self._pindex.distribution_strategy
-
+        hdf5_group['kindex'] = self.kindex
+        hdf5_group['rho'] = self.rho
+        hdf5_group['pundex'] = self.pundex
+        hdf5_group['logarithmic'] = self.config["logarithmic"]
+        # Store nbin as string, since it can be None
+        hdf5_group.attrs['nbin'] = str(self.config["nbin"])
+        hdf5_group.attrs['binbounds'] = str(self.config["binbounds"])
+
+        #MR FIXME: why not "return None" as happens everywhere else?
         return {
             'harmonic_partner': self.harmonic_partner,
+            'pindex': self.pindex,
+            'k_array': self.k_array
         }
 
     @classmethod
     def _from_hdf5(cls, hdf5_group, repository):
-        hp = repository.get('harmonic_partner', hdf5_group)
-        exec("bb = " + hdf5_group.attrs['binbounds'])
-        ds = hdf5_group.attrs['distribution_strategy']
-        return PowerSpace(hp, ds, binbounds=bb)
+        # make an empty PowerSpace object
+        new_ps = EmptyPowerSpace()
+        # reset class
+        new_ps.__class__ = cls
+        # call instructor so that classes are properly setup
+        super(PowerSpace, new_ps).__init__()
+        # set all values
+        new_ps._harmonic_partner = repository.get('harmonic_partner',
+                                                  hdf5_group)
+
+        new_ps._config = {}
+        new_ps._config['logarithmic'] = hdf5_group['logarithmic'][()]
+        exec("new_ps._config['nbin'] = " + hdf5_group.attrs['nbin'])
+        exec("new_ps._config['binbounds'] = " + hdf5_group.attrs['binbounds'])
+
+        new_ps._pindex = repository.get('pindex', hdf5_group)
+        new_ps._kindex = hdf5_group['kindex'][:]
+        new_ps._rho = hdf5_group['rho'][:]
+        new_ps._pundex = hdf5_group['pundex'][:]
+        new_ps._k_array = repository.get('k_array', hdf5_group)
+        new_ps._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',
+                                    '_k_array']
+
+        return new_ps
+
+
+class EmptyPowerSpace(PowerSpace):
+    def __init__(self):
+        pass
diff --git a/nifty/spaces/rg_space/rg_space.py b/nifty/spaces/rg_space/rg_space.py
index 8ca0f45ed..26e9f8bd6 100644
--- a/nifty/spaces/rg_space/rg_space.py
+++ b/nifty/spaces/rg_space/rg_space.py
@@ -283,36 +283,6 @@ class RGSpace(Space):
         dists = np.sqrt(dists)
         return dists
 
-    def get_unique_distances(self):
-        dimensions = len(self.shape)
-        if dimensions == 1:  # extra easy
-            maxdist = self.shape[0]//2
-            return np.arange(maxdist+1, dtype=np.float64) * self.distances[0]
-        if np.all(self.distances == self.distances[0]):  # shortcut
-            maxdist = np.asarray(self.shape)//2
-            tmp = np.sum(maxdist*maxdist)
-            tmp = np.zeros(tmp+1, dtype=np.bool)
-            t2 = np.arange(maxdist[0]+1, dtype=np.int64)
-            t2 *= t2
-            for i in range(1, dimensions):
-                t3 = np.arange(maxdist[i]+1, dtype=np.int64)
-                t3 *= t3
-                t2 = np.add.outer(t2, t3)
-            tmp[t2] = True
-            return np.sqrt(np.nonzero(tmp)[0])*self.distances[0]
-        else:  # do it the hard way
-            tmp = self.get_distance_array('not').unique()  # expensive!
-            tol = 1e-12*tmp[-1]
-            # remove all points that are closer than tol to their right
-            # neighbors.
-            # I'm appending the last value*2 to the array to treat the
-            # rightmost point correctly.
-            return tmp[np.diff(np.r_[tmp, 2*tmp[-1]]) > tol]
-
-    def get_natural_binbounds(self):
-        tmp = self.get_unique_distances()
-        return 0.5*(tmp[:-1]+tmp[1:])
-
     def get_fft_smoothing_kernel_function(self, sigma):
         return lambda x: np.exp(-2. * np.pi*np.pi * x*x * sigma*sigma)
 
diff --git a/nifty/spaces/space/space.py b/nifty/spaces/space/space.py
index 12b2b4a18..dc49c0c4c 100644
--- a/nifty/spaces/space/space.py
+++ b/nifty/spaces/space/space.py
@@ -117,20 +117,6 @@ class Space(DomainObject):
         raise NotImplementedError(
             "There is no generic distance structure for Space base class.")
 
-    def get_unique_distances(self):
-        raise NotImplementedError
-
-    def get_natural_binbounds(self):
-        """ The boundaries for natural power spectrum binning.
-
-        Returns
-        -------
-        distributed_data_object
-            A numpy array containing the binbounds
-
-        """
-        raise NotImplementedError
-
     def get_fft_smoothing_kernel_function(self, sigma):
         """ This method returns a smoothing kernel function.
 
diff --git a/test/test_spaces/test_power_space.py b/test/test_spaces/test_power_space.py
index 04f9a7452..256b38647 100644
--- a/test/test_spaces/test_power_space.py
+++ b/test/test_spaces/test_power_space.py
@@ -48,12 +48,14 @@ HARMONIC_SPACES = [RGSpace((8,), harmonic=True),
 
 #Try all sensible kinds of combinations of spaces, distributuion strategy and
 #binning parameters
+_maybe_fftw = ["fftw"] if ('pyfftw' in gdi) else []
+
 CONSISTENCY_CONFIGS_IMPLICIT = product(HARMONIC_SPACES,
                                        ["not", "equal", "fftw"],
                                        [None], [None, 3, 4], [True, False])
 CONSISTENCY_CONFIGS_EXPLICIT = product(HARMONIC_SPACES,
                                        ["not", "equal", "fftw"],
-                                       [[0., 1.3]], [None], [None])
+                                       [[0., 1.3]], [None], [False])
 CONSISTENCY_CONFIGS = chain(CONSISTENCY_CONFIGS_IMPLICIT,
                             CONSISTENCY_CONFIGS_EXPLICIT)
 
@@ -62,16 +64,18 @@ CONSISTENCY_CONFIGS = chain(CONSISTENCY_CONFIGS_IMPLICIT,
 CONSTRUCTOR_CONFIGS = [
     [1, 'not', False, None, None, {'error': ValueError}],
     [RGSpace((8,)), 'not', False, None, None, {'error': ValueError}],
-    [RGSpace((8,), harmonic=True), 'not', None, None, None, {
+    [RGSpace((8,), harmonic=True), 'not', False, None, None, {
         'harmonic': True,
         'shape': (5,),
         'dim': 5,
         'total_volume': 8.0,
         'harmonic_partner': RGSpace((8,), harmonic=True),
-        'binbounds': None,
+        'config': {'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.]),
         'rho': np.array([1, 2, 2, 2, 1]),
+        'pundex': np.array([0, 1, 2, 3, 4]),
+        'k_array': np.array([0., 1., 2., 3., 4., 3., 2., 1.]),
         }],
     [RGSpace((8,), harmonic=True), 'not', True, None, None, {
         'harmonic': True,
@@ -79,10 +83,13 @@ CONSTRUCTOR_CONFIGS = [
         'dim': 2,
         'total_volume': 8.0,
         'harmonic_partner': RGSpace((8,), harmonic=True),
-        'binbounds': (0.70710678118654757,),
+        'config': {'logarithmic': True, 'nbin': None, 'binbounds': None},
         'pindex': distributed_data_object([0, 1, 1, 1, 1, 1, 1, 1]),
         'kindex': np.array([0., 2.28571429]),
         'rho': np.array([1, 7]),
+        'pundex': np.array([0, 1]),
+        'k_array': np.array([0., 2.28571429, 2.28571429, 2.28571429,
+                             2.28571429, 2.28571429, 2.28571429, 2.28571429]),
         }],
     ]
 
@@ -115,10 +122,12 @@ def get_weight_configs():
 class PowerSpaceInterfaceTest(unittest.TestCase):
     @expand([
         ['harmonic_partner', Space],
-        ['binbounds', NoneType],
+        ['config', dict],
         ['pindex', distributed_data_object],
         ['kindex', np.ndarray],
         ['rho', np.ndarray],
+        ['pundex', np.ndarray],
+        ['k_array', distributed_data_object],
         ])
     def test_property_ret_type(self, attribute, expected_type):
         r = RGSpace((4, 4), harmonic=True)
@@ -128,15 +137,32 @@ class PowerSpaceInterfaceTest(unittest.TestCase):
 
 class PowerSpaceConsistencyCheck(unittest.TestCase):
     @expand(CONSISTENCY_CONFIGS)
-    def test_rhopindexConsistency(self, harmonic_partner, distribution_strategy,
-                         binbounds, nbin,logarithmic):
+    def test_pipundexInversion(self, harmonic_partner, distribution_strategy,
+                               binbounds, nbin, logarithmic):
+        if distribution_strategy == "fftw":
+            if not hasattr(gdi.get('fftw'), 'FFTW_MPI'):
+                raise SkipTest
+        p = PowerSpace(harmonic_partner=harmonic_partner,
+                       distribution_strategy=distribution_strategy,
+                       logarithmic=logarithmic, nbin=nbin,
+                       binbounds=binbounds)
+        assert_equal(p.pindex.flatten()[p.pundex], np.arange(p.dim),
+                     err_msg='pundex is not right-inverse of pindex!')
+
+    @expand(CONSISTENCY_CONFIGS)
+    def test_rhopindexConsistency(self, harmonic_partner,
+                                  distribution_strategy, binbounds, nbin,
+                                  logarithmic):
         if distribution_strategy == "fftw":
             if not hasattr(gdi.get('fftw'), 'FFTW_MPI'):
+                print (gdi.get('fftw'), "blub \n\n\n")
                 raise SkipTest
         p = PowerSpace(harmonic_partner=harmonic_partner,
                            distribution_strategy=distribution_strategy,
                            logarithmic=logarithmic, nbin=nbin,
                            binbounds=binbounds)
+        assert_equal(p.pindex.flatten()[p.pundex],np.arange(p.dim),
+            err_msg='pundex is not right-inverse of pindex!')
 
         assert_equal(p.pindex.flatten().bincount(), p.rho,
             err_msg='rho is not equal to pindex degeneracy')
@@ -148,6 +174,7 @@ class PowerSpaceFunctionalityTest(unittest.TestCase):
         if distribution_strategy == "fftw":
             if not hasattr(gdi.get('fftw'), 'FFTW_MPI'):
                 raise SkipTest
+            raise SkipTest
         if 'error' in expected:
             with assert_raises(expected['error']):
                 PowerSpace(harmonic_partner=harmonic_partner,
-- 
GitLab