Commit 95de2e53 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

step 1

parent 85b5d01c
...@@ -17,4 +17,3 @@ ...@@ -17,4 +17,3 @@
# and financially supported by the Studienstiftung des deutschen Volkes. # and financially supported by the Studienstiftung des deutschen Volkes.
from power_space import PowerSpace from power_space import PowerSpace
from power_index_factory import PowerIndexFactory
\ No newline at end of file
# 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()
...@@ -17,19 +17,14 @@ ...@@ -17,19 +17,14 @@
# and financially supported by the Studienstiftung des deutschen Volkes. # and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np import numpy as np
from d2o import distributed_data_object,\ 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): class PowerIndices(object):
"""Computes helpful quantities to deal with power spectra. """Computes helpful quantities to deal with power spectra.
Given the shape and the density of a underlying rectangular grid this Given the shape and the density of a underlying rectangular grid this
class provides the user class provides the user
with the pindex, kindex, rho and pundex. The indices are binned with the pindex, kindex and rho. The indices are binned
according to the supplied parameter scheme. If wanted, computed according to the supplied parameter scheme. If wanted, computed
results are stored for future reuse. results are stored for future reuse.
...@@ -150,7 +145,7 @@ class PowerIndices(object): ...@@ -150,7 +145,7 @@ class PowerIndices(object):
def get_index_dict(self, **kwargs): def get_index_dict(self, **kwargs):
""" """
Returns a dictionary containing the pindex, kindex, rho and pundex Returns a dictionary containing the pindex, kindex and rho
binned according to the supplied parameter scheme and a binned according to the supplied parameter scheme and a
configuration dict containing this scheme. configuration dict containing this scheme.
...@@ -170,8 +165,7 @@ class PowerIndices(object): ...@@ -170,8 +165,7 @@ class PowerIndices(object):
Returns Returns
------- -------
index_dict : dict index_dict : dict
Contains the keys: 'config', 'pindex', 'kindex', 'rho' and Contains the keys: 'config', 'pindex', 'kindex' and 'rho'
'pundex'
""" """
# Cast the input arguments # Cast the input arguments
temp_config_dict = self._cast_config(**kwargs) temp_config_dict = self._cast_config(**kwargs)
...@@ -207,18 +201,18 @@ class PowerIndices(object): ...@@ -207,18 +201,18 @@ class PowerIndices(object):
def _compute_index_dict(self, config_dict): def _compute_index_dict(self, config_dict):
""" """
Internal helper function which takes a config_dict, asks for the Internal helper function which takes a config_dict, asks for the
pindex/kindex/rho/pundex set, and bins them according to the config pindex/kindex/rho set, and bins them according to the config
""" """
# if no binning is requested, compute the indices, build the dict, # if no binning is requested, compute the indices, build the dict,
# and return it straight. # and return it straight.
if not config_dict["logarithmic"] and config_dict["nbin"] is None and \ if not config_dict["logarithmic"] and config_dict["nbin"] is None and \
config_dict["binbounds"] is None: config_dict["binbounds"] is None:
(temp_pindex, temp_kindex, temp_rho, temp_pundex) =\ (temp_pindex, temp_kindex, temp_rho) =\
self._compute_indices(self.k_array) self._compute_indices(self.k_array)
temp_k_array = self.k_array temp_k_array = self.k_array
# if binning is required, make a recursive call to get the unbinned # if binning is required, make a recursive call to get the unbinned
# indices, bin them, compute the pundex and then return everything. # indices, bin them, and then return everything.
else: else:
# Get the unbinned indices # Get the unbinned indices
temp_unbinned_indices = self.get_index_dict(nbin=None, temp_unbinned_indices = self.get_index_dict(nbin=None,
...@@ -226,7 +220,7 @@ class PowerIndices(object): ...@@ -226,7 +220,7 @@ class PowerIndices(object):
logarithmic=False, logarithmic=False,
store=False) store=False)
# Bin them # Bin them
(temp_pindex, temp_kindex, temp_rho, temp_pundex) = \ (temp_pindex, temp_kindex, temp_rho) = \
self._bin_power_indices( self._bin_power_indices(
temp_unbinned_indices, **config_dict) temp_unbinned_indices, **config_dict)
# Make a binned version of k_array # Make a binned version of k_array
...@@ -237,7 +231,6 @@ class PowerIndices(object): ...@@ -237,7 +231,6 @@ class PowerIndices(object):
"pindex": temp_pindex, "pindex": temp_pindex,
"kindex": temp_kindex, "kindex": temp_kindex,
"rho": temp_rho, "rho": temp_rho,
"pundex": temp_pundex,
"k_array": temp_k_array} "k_array": temp_k_array}
return temp_index_dict return temp_index_dict
...@@ -249,7 +242,7 @@ class PowerIndices(object): ...@@ -249,7 +242,7 @@ class PowerIndices(object):
def _compute_indices(self, k_array): def _compute_indices(self, k_array):
""" """
Internal helper function which computes pindex, kindex, rho and pundex Internal helper function which computes pindex, kindex and rho
from a given k_array from a given k_array
""" """
########## ##########
...@@ -275,63 +268,7 @@ class PowerIndices(object): ...@@ -275,63 +268,7 @@ class PowerIndices(object):
####### #######
global_rho = global_pindex.bincount().get_full_data() global_rho = global_pindex.bincount().get_full_data()
########## return global_pindex, global_kindex, global_rho
# 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): def _bin_power_indices(self, index_dict, **kwargs):
""" """
...@@ -356,7 +293,7 @@ class PowerIndices(object): ...@@ -356,7 +293,7 @@ class PowerIndices(object):
Returns Returns
------- -------
pindex : distributed_data_object pindex : distributed_data_object
kindex, rho, pundex : ndarrays kindex, rho : ndarrays
The (re)binned power indices. The (re)binned power indices.
""" """
...@@ -411,6 +348,5 @@ class PowerIndices(object): ...@@ -411,6 +348,5 @@ class PowerIndices(object):
pindex_ = pindex.copy_empty() pindex_ = pindex.copy_empty()
pindex_.set_local_data(reorder[pindex.get_local_data()]) pindex_.set_local_data(reorder[pindex.get_local_data()])
pundex_ = self._compute_pundex(pindex_, kindex_)
return pindex_, kindex_, rho_, pundex_ return pindex_, kindex_, rho_
...@@ -20,9 +20,8 @@ import numpy as np ...@@ -20,9 +20,8 @@ import numpy as np
import d2o import d2o
from power_index_factory import PowerIndexFactory
from nifty.spaces.space import Space from nifty.spaces.space import Space
from power_indices import PowerIndices
class PowerSpace(Space): class PowerSpace(Space):
...@@ -57,9 +56,6 @@ class PowerSpace(Space): ...@@ -57,9 +56,6 @@ class PowerSpace(Space):
mapped to which power bin mapped to which power bin
kindex : numpy.ndarray kindex : numpy.ndarray
Sorted array of all k-modes. 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 rho : numpy.ndarray
The amount of k-modes that get mapped to one power bin is given by The amount of k-modes that get mapped to one power bin is given by
rho. rho.
...@@ -88,8 +84,7 @@ class PowerSpace(Space): ...@@ -88,8 +84,7 @@ class PowerSpace(Space):
distribution_strategy='not', distribution_strategy='not',
logarithmic=False, nbin=None, binbounds=None): logarithmic=False, nbin=None, binbounds=None):
super(PowerSpace, self).__init__() super(PowerSpace, self).__init__()
self._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex', self._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_k_array']
'_k_array']
if not isinstance(harmonic_partner, Space): if not isinstance(harmonic_partner, Space):
raise ValueError( raise ValueError(
...@@ -99,19 +94,19 @@ class PowerSpace(Space): ...@@ -99,19 +94,19 @@ class PowerSpace(Space):
"harmonic_partner must be a harmonic space.") "harmonic_partner must be a harmonic space.")
self._harmonic_partner = harmonic_partner self._harmonic_partner = harmonic_partner
power_index = PowerIndexFactory.get_power_index( tmp = PowerIndices(self.harmonic_partner, distribution_strategy,
domain=self.harmonic_partner, logarithmic=logarithmic,
distribution_strategy=distribution_strategy, nbin=nbin,
logarithmic=logarithmic, binbounds=binbounds)
nbin=nbin, power_index = tmp.get_index_dict(logarithmic=logarithmic,
binbounds=binbounds) nbin=nbin,
binbounds=binbounds)
self._config = power_index['config'] self._config = power_index['config']
self._pindex = power_index['pindex'] self._pindex = power_index['pindex']
self._kindex = power_index['kindex'] self._kindex = power_index['kindex']
self._rho = power_index['rho'] self._rho = power_index['rho']
self._pundex = power_index['pundex']
self._k_array = power_index['k_array'] self._k_array = power_index['k_array']
if self.config['nbin'] is not None: if self.config['nbin'] is not None:
...@@ -242,14 +237,6 @@ class PowerSpace(Space): ...@@ -242,14 +237,6 @@ class PowerSpace(Space):
""" """
return self._rho 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 @property
def k_array(self): def k_array(self):
""" An array containing distances to the grid center (i.e. zero-mode) """ An array containing distances to the grid center (i.e. zero-mode)
...@@ -262,7 +249,6 @@ class PowerSpace(Space): ...@@ -262,7 +249,6 @@ class PowerSpace(Space):
def _to_hdf5(self, hdf5_group): def _to_hdf5(self, hdf5_group):
hdf5_group['kindex'] = self.kindex hdf5_group['kindex'] = self.kindex
hdf5_group['rho'] = self.rho hdf5_group['rho'] = self.rho
hdf5_group['pundex'] = self.pundex
hdf5_group['logarithmic'] = self.config["logarithmic"] hdf5_group['logarithmic'] = self.config["logarithmic"]
# Store nbin as string, since it can be None # Store nbin as string, since it can be None
hdf5_group.attrs['nbin'] = str(self.config["nbin"]) hdf5_group.attrs['nbin'] = str(self.config["nbin"])
...@@ -295,10 +281,8 @@ class PowerSpace(Space): ...@@ -295,10 +281,8 @@ class PowerSpace(Space):
new_ps._pindex = repository.get('pindex', hdf5_group) new_ps._pindex = repository.get('pindex', hdf5_group)
new_ps._kindex = hdf5_group['kindex'][:] new_ps._kindex = hdf5_group['kindex'][:]
new_ps._rho = hdf5_group['rho'][:] new_ps._rho = hdf5_group['rho'][:]
new_ps._pundex = hdf5_group['pundex'][:]
new_ps._k_array = repository.get('k_array', hdf5_group) new_ps._k_array = repository.get('k_array', hdf5_group)
new_ps._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex', new_ps._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_k_array']
'_k_array']
return new_ps return new_ps
......
...@@ -32,20 +32,20 @@ from itertools import product, chain ...@@ -32,20 +32,20 @@ from itertools import product, chain
from d2o.config import dependency_injector as gdi from d2o.config import dependency_injector as gdi
HARMONIC_SPACES = [RGSpace((8,), harmonic=True), HARMONIC_SPACES = [RGSpace((8,), harmonic=True),
RGSpace((7,), harmonic=True,zerocenter=True), RGSpace((7,), harmonic=True,zerocenter=True),
RGSpace((8,), harmonic=True,zerocenter=True), RGSpace((8,), harmonic=True,zerocenter=True),
RGSpace((7,8), harmonic=True), RGSpace((7,8), harmonic=True),
RGSpace((7,8), harmonic=True, zerocenter=True), RGSpace((7,8), harmonic=True, zerocenter=True),
RGSpace((6,6), harmonic=True, zerocenter=True), RGSpace((6,6), harmonic=True, zerocenter=True),
RGSpace((7,5), harmonic=True, zerocenter=True), RGSpace((7,5), harmonic=True, zerocenter=True),
RGSpace((5,5), harmonic=True), RGSpace((5,5), harmonic=True),
RGSpace((4,5,7), harmonic=True), RGSpace((4,5,7), harmonic=True),
RGSpace((4,5,7), harmonic=True, zerocenter=True), RGSpace((4,5,7), harmonic=True, zerocenter=True),
LMSpace(6), LMSpace(6),
LMSpace(9)] LMSpace(9)]
#Try all sensible kinds of combinations of spaces, distributuion strategy and #Try all sensible kinds of combinations of spaces, distributuion strategy and
#binning parameters #binning parameters
_maybe_fftw = ["fftw"] if ('pyfftw' in gdi) else [] _maybe_fftw = ["fftw"] if ('pyfftw' in gdi) else []
...@@ -68,7 +68,6 @@ CONSTRUCTOR_CONFIGS = [ ...@@ -68,7 +68,6 @@ CONSTRUCTOR_CONFIGS = [
'pindex': distributed_data_object([0, 1, 2, 3, 4, 3, 2, 1]), 'pindex': distributed_data_object([0, 1, 2, 3, 4, 3, 2, 1]),
'kindex': np.array([0., 1., 2., 3., 4.]), 'kindex': np.array([0., 1., 2., 3., 4.]),
'rho': np.array([1, 2, 2, 2, 1]), '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.]), 'k_array': np.array([0., 1., 2., 3., 4., 3., 2., 1.]),
}], }],
[RGSpace((8,), harmonic=True), 'not', True, None, None, { [RGSpace((8,), harmonic=True), 'not', True, None, None, {
...@@ -81,7 +80,6 @@ CONSTRUCTOR_CONFIGS = [ ...@@ -81,7 +80,6 @@ CONSTRUCTOR_CONFIGS = [
'pindex': distributed_data_object([0, 1, 1, 1, 1, 1, 1, 1]), 'pindex': distributed_data_object([0, 1, 1, 1, 1, 1, 1, 1]),
'kindex': np.array([0., 2.28571429]), 'kindex': np.array([0., 2.28571429]),
'rho': np.array([1, 7]), 'rho': np.array([1, 7]),
'pundex': np.array([0, 1]),
'k_array': np.array([0., 2.28571429, 2.28571429, 2.28571429, 'k_array': np.array([0., 2.28571429, 2.28571429, 2.28571429,
2.28571429, 2.28571429, 2.28571429, 2.28571429]), 2.28571429, 2.28571429, 2.28571429, 2.28571429]),
}], }],
...@@ -120,7 +118,6 @@ class PowerSpaceInterfaceTest(unittest.TestCase): ...@@ -120,7 +118,6 @@ class PowerSpaceInterfaceTest(unittest.TestCase):
['pindex', distributed_data_object], ['pindex', distributed_data_object],
['kindex', np.ndarray], ['kindex', np.ndarray],
['rho', np.ndarray], ['rho', np.ndarray],
['pundex', np.ndarray],
['k_array', distributed_data_object], ['k_array', distributed_data_object],
]) ])
def test_property_ret_type(self, attribute, expected_type): def test_property_ret_type(self, attribute, expected_type):
...@@ -129,22 +126,16 @@ class PowerSpaceInterfaceTest(unittest.TestCase): ...@@ -129,22 +126,16 @@ class PowerSpaceInterfaceTest(unittest.TestCase):
assert_(isinstance(getattr(p, attribute), expected_type)) assert_(isinstance(getattr(p, attribute), expected_type))
class PowerSpaceConsistencyCheck(unittest.TestCase): class PowerSpaceConsistencyCheck(unittest.TestCase):
@expand(CONSISTENCY_CONFIGS)
def test_pipundexInversion(self, harmonic_partner, distribution_strategy,
binbounds, nbin,logarithmic):
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) @expand(CONSISTENCY_CONFIGS)
def test_rhopindexConsistency(self, harmonic_partner, distribution_strategy, def test_rhopindexConsistency(self, harmonic_partner, distribution_strategy,
binbounds, nbin,logarithmic): binbounds, nbin,logarithmic):
p = PowerSpace(harmonic_partner=harmonic_partner,
distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
assert_equal(p.pindex.flatten().bincount(), p.rho, assert_equal(p.pindex.flatten().bincount(), p.rho,
err_msg='rho is not equal to pindex degeneracy') err_msg='rho is not equal to pindex degeneracy')
class PowerSpaceFunctionalityTest(unittest.TestCase): class PowerSpaceFunctionalityTest(unittest.TestCase):
@expand(CONSISTENCY_CONFIGS) @expand(CONSISTENCY_CONFIGS)
def test_constructor(self, harmonic_partner, distribution_strategy, def test_constructor(self, harmonic_partner, distribution_strategy,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment