From a929944b2661f2f75db8585eb690c20e4f24f26f Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Fri, 2 Jun 2017 21:55:14 +0200
Subject: [PATCH] still one failing test
---
nifty/spaces/power_space/power_indices.py | 146 ----------------------
nifty/spaces/power_space/power_space.py | 141 ++++++++++++---------
test/test_spaces/test_power_space.py | 8 +-
3 files changed, 82 insertions(+), 213 deletions(-)
delete mode 100644 nifty/spaces/power_space/power_indices.py
diff --git a/nifty/spaces/power_space/power_indices.py b/nifty/spaces/power_space/power_indices.py
deleted file mode 100644
index df67b9bae..000000000
--- a/nifty/spaces/power_space/power_indices.py
+++ /dev/null
@@ -1,146 +0,0 @@
-# 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_
diff --git a/nifty/spaces/power_space/power_space.py b/nifty/spaces/power_space/power_space.py
index b6244f4a6..3e0d5f851 100644
--- a/nifty/spaces/power_space/power_space.py
+++ b/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):
diff --git a/test/test_spaces/test_power_space.py b/test/test_spaces/test_power_space.py
index 94e21bc3f..ea704dd24 100644
--- a/test/test_spaces/test_power_space.py
+++ b/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:
--
GitLab