merge step 1

nifty/field.py

@@ -270,7 +270,7 @@ class Field(Loggable, Versionable, object):
 
     # ---Powerspectral methods---
 
-    def power_analyze(self, spaces=None, logarithmic=False, nbin=None,
+    def power_analyze(self, spaces=None, logarithmic=None, nbin=None,
                       binbounds=None, keep_phase_information=False):
         """ Computes the square root power spectrum for a subspace of `self`.
 
@@ -287,14 +287,15 @@ class Field(Loggable, Versionable, object):
             (default : None).
         logarithmic : boolean *optional*
             True if the output PowerSpace should use logarithmic binning.
-            {default : False}
+            {default : None}
         nbin : int *optional*
             The number of bins the resulting PowerSpace shall have
             (default : None).
             if nbin==None : maximum number of bins is used
         binbounds : array-like *optional*
             Inner bounds of the bins (default : None).
-            if binbounds==None : bins are inferred. Overwrites nbins and log
+            Overrides nbin and logarithmic.
+            if binbounds==None : bins are inferred.
         keep_phase_information : boolean, *optional*
             If False, return a real-valued result containing the power spectrum
             of the input Field.
@@ -397,14 +398,9 @@ 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,
-                                pindex=pindex,
-                                rho=rho,
+                                pdomain=power_domain,
                                 axes=work_field.domain_axes[space_index])
 
         # create the result field and put power_spectrum into it
@@ -421,8 +417,11 @@ class Field(Loggable, Versionable, object):
         return result_field
 
     @classmethod
-    def _calculate_power_spectrum(cls, field_val, pindex, rho, axes=None):
+    def _calculate_power_spectrum(cls, field_val, pdomain, 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,
@@ -431,6 +430,7 @@ 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 dimansions of the spaces
+            The output object. The tuple contains the dimensions of the spaces
             in domain.
 
         See Also

nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py

@@ -38,7 +38,7 @@ class InvertibleOperatorMixin(object):
         (default: ConjugateGradient)
 
     preconditioner : LinearOperator
-        Preconditioner that is used by ConjugateGraduent if no minimizer was
+        Preconditioner that is used by ConjugateGradient if no minimizer was
         given.
 
     Attributes

nifty/operators/response_operator/response_operator.py

@@ -33,8 +33,7 @@ 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. As the Operator
-        is endomorphic this is the same as its domain.
+        The domain in which the outcome of the operator lives.
     unitary : boolean
         Indicates whether the Operator is unitary or not.
 

nifty/spaces/lm_space/lm_space.py

@@ -143,6 +143,22 @@ 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)
 
@@ -152,6 +168,12 @@ 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

nifty/spaces/power_space/__init__.py

@@ -17,4 +17,3 @@
 # 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

nifty/spaces/power_space/power_index_factory.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.
-
-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()

nifty/spaces/rg_space/rg_space.py

@@ -267,6 +267,36 @@ 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)
 

nifty/spaces/space/space.py

@@ -117,6 +117,20 @@ 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.
 

test/test_spaces/test_power_space.py

@@ -48,14 +48,12 @@ 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], [False])
+                                       [[0., 1.3]], [None], [None])
 CONSISTENCY_CONFIGS = chain(CONSISTENCY_CONFIGS_IMPLICIT,
                             CONSISTENCY_CONFIGS_EXPLICIT)
 
@@ -64,18 +62,16 @@ 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', False, None, None, {
+    [RGSpace((8,), harmonic=True), 'not', None, None, None, {
         'harmonic': True,
         'shape': (5,),
         'dim': 5,
         'total_volume': 8.0,
         'harmonic_partner': RGSpace((8,), harmonic=True),
-        'config': {'logarithmic': False, 'nbin': None, 'binbounds': 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,
@@ -83,13 +79,10 @@ CONSTRUCTOR_CONFIGS = [
         'dim': 2,
         'total_volume': 8.0,
         'harmonic_partner': RGSpace((8,), harmonic=True),
-        'config': {'logarithmic': True, 'nbin': None, 'binbounds': None},
+        'binbounds': (0.70710678118654757,),
         '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]),
         }],
     ]
 
@@ -122,12 +115,10 @@ def get_weight_configs():
 class PowerSpaceInterfaceTest(unittest.TestCase):
     @expand([
         ['harmonic_partner', Space],
-        ['config', dict],
+        ['binbounds', NoneType],
         ['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)
@@ -137,32 +128,15 @@ class PowerSpaceInterfaceTest(unittest.TestCase):
 
 class PowerSpaceConsistencyCheck(unittest.TestCase):
     @expand(CONSISTENCY_CONFIGS)
-    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):
+    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')
@@ -174,7 +148,6 @@ 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,