Commit f4cc3157 authored by theos's avatar theos
Browse files

Intermediate commit for merging.

parent 36082174
......@@ -43,16 +43,17 @@ class space_paradict(object):
temp_hash = hash(item)
except TypeError:
temp_hash = hash(tuple(item))
result_hash ^= temp_hash ^ 131*hash(key)
result_hash ^= temp_hash ^ int(hash(key)/131)
return result_hash
class rg_space_paradict(space_paradict):
def __init__(self, shape, zerocenter, distances):
def __init__(self, shape, zerocenter, distances, harmonic):
if not hasattr(self, 'parameters'):
self.parameters = {}
self.parameters.__setitem__('shape', shape)
self.parameters.__setitem__('harmonic', harmonic)
space_paradict.__init__(
self, zerocenter=zerocenter, distances=distances)
......@@ -64,17 +65,24 @@ class rg_space_paradict(space_paradict):
if key == 'shape':
temp = np.empty(len(self['shape']), dtype=np.int)
temp[:] = arg
temp = tuple(temp)
elif key == 'zerocenter':
temp = np.empty(len(self['shape']), dtype=bool)
temp[:] = arg
temp = tuple(temp)
elif key == 'distances':
if arg is None:
temp = 1 / np.array(self['shape'], dtype=np.float)
if self['harmonic']:
temp = np.ones_like(self['shape'], dtype=np.float)
else:
temp = 1 / np.array(self['shape'], dtype=np.float)
else:
temp = np.empty(len(self['shape']), dtype=np.float)
temp[:] = arg
temp = tuple(temp)
elif key == 'harmonic':
temp = bool(arg)
temp = tuple(temp)
self.parameters.__setitem__(key, temp)
......@@ -185,20 +193,32 @@ class hp_space_paradict(space_paradict):
class power_space_paradict(space_paradict):
def __init__(self, power_indices):
def __init__(self, pindex, kindex, rho, pundex, k_array, config,
harmonic_domain):
space_paradict.__init__(self,
power_indices=power_indices)
pindex=pindex,
kindex=kindex,
rho=rho,
pundex=pundex,
k_array=k_array,
config=config,
harmonic_domain=harmonic_domain)
def __setitem__(self, key, arg):
if key not in ['power_indices']:
if key not in ['pindex', 'kindex', 'rho', 'config', 'harmonic_domain']:
raise ValueError(about._errors.cstring(
"ERROR: Unsupported PowerSpace parameter: " + key))
if key == 'power_indices':
if key == 'harmonic_domain':
if not arg.harmonic:
raise ValueError(about._errors.cstring(
"ERROR: harmonic_domain must be harmonic."))
temp = arg
else:
temp = arg
self.parameters.__setitem__(key, temp)
def __hash__(self):
return (self.parameters['power_indices']['hash'] ^
131*hash('power_indices'))
return (hash(frozenset(self.parameters['config'].items())) ^
(hash(self.parameters['domain'])/131))
# -*- coding: utf-8 -*-
from nifty.power import PowerSpace
class LMPowerSpace(PowerSpace):
pass
......@@ -6,8 +6,8 @@ class _PowerIndexFactory(object):
def __init__(self):
self.power_indices_storage = {}
def get_power_indices(self, log, nbin, binbounds,
domain, distribution_strategy):
def get_power_indices(self, domain, distribution_strategy,
log=False, nbin=None, binbounds=None):
current_hash = domain.__hash__() ^ (111*hash(distribution_strategy))
if current_hash not in self.power_indices_storage:
......
......@@ -44,8 +44,8 @@ class PowerIndices(object):
self.domain = domain
self.distribution_strategy = distribution_strategy
# Compute the global kdict
self.kdict = self.domain.compute_k_array(distribution_strategy)
# Compute the global k_array
self.k_array = self.domain.compute_k_array(distribution_strategy)
# Initialize the dictonary which stores all individual index-dicts
self.global_dict = {}
# Set self.default_parameters
......@@ -137,10 +137,10 @@ class PowerIndices(object):
"binbounds": temp_binbounds}
return temp_dict
def compute_kdict(self):
def compute_k_array(self):
raise NotImplementedError(
about._errors.cstring(
"ERROR: No generic compute_kdict method implemented."))
"ERROR: No generic compute_k_array method implemented."))
def get_index_dict(self, **kwargs):
"""
......@@ -208,8 +208,8 @@ class PowerIndices(object):
if not config_dict["log"] and config_dict["nbin"] is None and \
config_dict["binbounds"] is None:
(temp_pindex, temp_kindex, temp_rho, temp_pundex) =\
self._compute_indices(self.kdict)
temp_kdict = self.kdict
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.
......@@ -223,46 +223,42 @@ class PowerIndices(object):
(temp_pindex, temp_kindex, temp_rho, temp_pundex) = \
self._bin_power_indices(
temp_unbinned_indices, **config_dict)
# Make a binned version of kdict
temp_kdict = self._compute_kdict_from_pindex_kindex(temp_pindex,
temp_kindex)
# Make a binned version of k_array
temp_k_array = self._compute_k_array_from_pindex_kindex(
temp_pindex, temp_kindex)
index_dict_hash = (hash(self.domain) ^
127*hash(self._freeze_config(config_dict)))
temp_index_dict = {"config": config_dict,
"pindex": temp_pindex,
"kindex": temp_kindex,
"rho": temp_rho,
"pundex": temp_pundex,
"kdict": temp_kdict,
"domain": self.domain,
"hash": index_dict_hash}
"k_array": temp_k_array}
return temp_index_dict
def _compute_kdict_from_pindex_kindex(self, pindex, kindex):
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, nkdict):
def _compute_indices(self, k_array):
"""
Internal helper function which computes pindex, kindex, rho and pundex
from a given nkdict
from a given k_array
"""
##########
# kindex #
##########
global_kindex = nkdict.unique()
global_kindex = k_array.unique()
##########
# pindex #
##########
# compute the local pindex slice on basis of the local nkdict data
local_pindex = np.searchsorted(global_kindex, nkdict.get_local_data())
# 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=nkdict.shape,
global_shape=k_array.shape,
dtype=local_pindex.dtype,
distribution_strategy=self.distribution_strategy)
# store the local pindex data in the global_pindex d2o
......@@ -411,133 +407,133 @@ class PowerIndices(object):
return pindex_, kindex_, rho_, pundex_
class LMPowerIndices(PowerIndices):
def __init__(self, lmax, dim, distribution_strategy,
log=False, nbin=None, binbounds=None):
"""
Returns an instance of the PowerIndices class. Given the shape and
the density of a underlying rectangular grid it provides the user
with the pindex, kindex, rho and pundex. The indices are bined
according to the supplied parameter scheme. If wanted, computed
results are stored for future reuse.
Parameters
----------
shape : tuple, list, ndarray
Array-like object which specifies the shape of the underlying
rectangular grid
dgrid : tuple, list, ndarray
Array-like object which specifies the step-width of the
underlying grid
zerocenter : boolean, tuple/list/ndarray of boolean *optional*
Specifies which dimensions are zerocentered. (default:False)
log : 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.
"""
# Basic inits and consistency checks
self.lmax = np.uint(lmax)
self.dim = np.uint(dim)
super(LMPowerIndices, self).__init__(
distribution_strategy=distribution_strategy,
log=log,
nbin=nbin,
binbounds=binbounds)
def compute_kdict(self):
"""
Calculates an n-dimensional array with its entries being the
lengths of the k-vectors from the zero point of the grid.
Parameters
----------
None : All information is taken from the parent object.
Returns
-------
nkdict : distributed_data_object
"""
if self.distribution_strategy != 'not':
about.warnings.cprint(
"WARNING: full kdict is temporarily stored on every node " +
"altough disribution strategy != 'not'!")
if 'healpy' in gdi:
nkdict = hp.Alm.getlm(self.lmax, i=None)[0]
else:
nkdict = self._getlm()[0]
nkdict = distributed_data_object(
nkdict,
distribution_strategy=self.distribution_strategy,
comm=self._comm)
return nkdict
def _getlm(self): # > compute all (l,m)
index = np.arange(self.dim)
n = 2 * self.lmax + 1
m = np.ceil((n - np.sqrt(n**2 - 8 * (index - self.lmax))) / 2
).astype(np.int)
l = index - self.lmax * m + m * (m - 1) // 2
return l, m
def _compute_indices_d2o(self, nkdict):
"""
Internal helper function which computes pindex, kindex, rho and pundex
from a given nkdict
"""
##########
# kindex #
##########
kindex = np.arange(self.lmax + 1, dtype=np.float)
##########
# pindex #
##########
pindex = nkdict.copy(dtype=np.int)
#######
# rho #
#######
rho = (2 * kindex + 1).astype(np.int)
##########
# pundex #
##########
pundex = self._compute_pundex_d2o(pindex, kindex)
return pindex, kindex, rho, pundex
def _compute_indices_np(self, nkdict):
"""
Internal helper function which computes pindex, kindex, rho and pundex
from a given nkdict
"""
##########
# kindex #
##########
kindex = np.arange(self.lmax + 1, dtype=np.float)
##########
# pindex #
##########
pindex = nkdict.astype(dtype=np.int, copy=True)
#######
# rho #
#######
rho = (2 * kindex + 1).astype(np.int)
##########
# pundex #
##########
pundex = self._compute_pundex_np(pindex, kindex)
return pindex, kindex, rho, pundex
#
#class LMPowerIndices(PowerIndices):
#
# def __init__(self, lmax, dim, distribution_strategy,
# log=False, nbin=None, binbounds=None):
# """
# Returns an instance of the PowerIndices class. Given the shape and
# the density of a underlying rectangular grid it provides the user
# with the pindex, kindex, rho and pundex. The indices are bined
# according to the supplied parameter scheme. If wanted, computed
# results are stored for future reuse.
#
# Parameters
# ----------
# shape : tuple, list, ndarray
# Array-like object which specifies the shape of the underlying
# rectangular grid
# dgrid : tuple, list, ndarray
# Array-like object which specifies the step-width of the
# underlying grid
# zerocenter : boolean, tuple/list/ndarray of boolean *optional*
# Specifies which dimensions are zerocentered. (default:False)
# log : 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.
# """
# # Basic inits and consistency checks
# self.lmax = np.uint(lmax)
# self.dim = np.uint(dim)
# super(LMPowerIndices, self).__init__(
# distribution_strategy=distribution_strategy,
# log=log,
# nbin=nbin,
# binbounds=binbounds)
#
# def compute_kdict(self):
# """
# Calculates an n-dimensional array with its entries being the
# lengths of the k-vectors from the zero point of the grid.
#
# Parameters
# ----------
# None : All information is taken from the parent object.
#
# Returns
# -------
# nkdict : distributed_data_object
# """
#
# if self.distribution_strategy != 'not':
# about.warnings.cprint(
# "WARNING: full kdict is temporarily stored on every node " +
# "altough disribution strategy != 'not'!")
#
# if 'healpy' in gdi:
# nkdict = hp.Alm.getlm(self.lmax, i=None)[0]
# else:
# nkdict = self._getlm()[0]
# nkdict = distributed_data_object(
# nkdict,
# distribution_strategy=self.distribution_strategy,
# comm=self._comm)
# return nkdict
#
# def _getlm(self): # > compute all (l,m)
# index = np.arange(self.dim)
# n = 2 * self.lmax + 1
# m = np.ceil((n - np.sqrt(n**2 - 8 * (index - self.lmax))) / 2
# ).astype(np.int)
# l = index - self.lmax * m + m * (m - 1) // 2
# return l, m
#
# def _compute_indices_d2o(self, nkdict):
# """
# Internal helper function which computes pindex, kindex, rho and pundex
# from a given nkdict
# """
# ##########
# # kindex #
# ##########
# kindex = np.arange(self.lmax + 1, dtype=np.float)
#
# ##########
# # pindex #
# ##########
# pindex = nkdict.copy(dtype=np.int)
#
# #######
# # rho #
# #######
# rho = (2 * kindex + 1).astype(np.int)
#
# ##########
# # pundex #
# ##########
# pundex = self._compute_pundex_d2o(pindex, kindex)
#
# return pindex, kindex, rho, pundex
#
# def _compute_indices_np(self, nkdict):
# """
# Internal helper function which computes pindex, kindex, rho and pundex
# from a given nkdict
# """
# ##########
# # kindex #
# ##########
# kindex = np.arange(self.lmax + 1, dtype=np.float)
#
# ##########
# # pindex #
# ##########
# pindex = nkdict.astype(dtype=np.int, copy=True)
#
# #######
# # rho #
# #######
# rho = (2 * kindex + 1).astype(np.int)
#
# ##########
# # pundex #
# ##########
# pundex = self._compute_pundex_np(pindex, kindex)
#
# return pindex, kindex, rho, pundex
......@@ -8,19 +8,25 @@ from nifty.nifty_paradict import power_space_paradict
class PowerSpace(Space):
def __init__(self, power_indices, dtype=np.dtype('float')):
def __init__(self, pindex, kindex, rho, config,
harmonic_domain, dtype=np.dtype('float'), **kwargs):
# the **kwargs is in the __init__ in order to enable a
# PowerSpace(**power_index) initialization
self.dtype = np.dtype(dtype)
self.paradict = power_space_paradict(power_indices=power_indices)
self.harmonic = True
self.paradict = power_space_paradict(pindex=pindex,
kindex=kindex,
rho=rho,
config=config,
harmonic_domain=harmonic_domain)
self._harmonic = True
@property
def shape(self):
return tuple(self.paradict['shape'])
return self.paradict['kindex'].shape
def calculate_power_spectrum(self, x, axes=None):
fieldabs = abs(x)**2
pindex = self.power_indices['pindex']
pindex = self.paradict['pindex']
if axes is not None:
pindex = self._shape_up_pindex(
pindex=pindex,
......@@ -30,7 +36,7 @@ class PowerSpace(Space):
power_spectrum = pindex.bincount(weights=fieldabs,
axis=axes)
rho = self.power_indices['rho']
rho = self.paradict['rho']
if axes is not None:
new_rho_shape = [1, ] * len(power_spectrum.shape)
new_rho_shape[axes[0]] = len(rho)
......
# -*- coding: utf-8 -*-
import numpy as np
from d2o import STRATEGIES
from nifty.power import PowerSpace
from nifty.nifty_paradict import rg_power_space_paradict
from power_index_factory import RGPowerIndexFactory
class RGPowerSpace(PowerSpace):
def __init__(self, distribution_strategy, dtype=np.dtype('float'),
log=False, nbin=None, binbounds=None):
self.dtype = np.dtype(dtype)
self.paradict = rg_power_space_paradict(
distribution_strategy=distribution_strategy,
log=log,
nbin=nbin,
binbounds=binbounds)
temp_dict = self.paradict.parameters.copy()
del temp_dict['complexity']
self.power_indices = RGPowerIndexFactory.get_power_indices(**temp_dict)
self.distances = (tuple(self.power_indices['rho']),)
self.harmonic = True
@property
def shape(self):
return tuple(self.paradict['shape'])
def calculate_power_spectrum(self, x, axes=None):
fieldabs = abs(x)**2
pindex = self.power_indices['pindex']
if axes is not None:
pindex = self._shape_up_pindex(
pindex=pindex,
target_shape=x.shape,
target_strategy=x.distribution_strategy,
axes=axes)
power_spectrum = pindex.bincount(weights=fieldabs,
axis=axes)
rho = self.power_indices['rho']
if axes is not None:
new_rho_shape = [1, ] * len(power_spectrum.shape)
new_rho_shape[axes[0]] = len(rho)
rho = rho.reshape(new_rho_shape)
power_spectrum /= rho
return power_spectrum
def _shape_up_pindex(self, pindex, target_shape, target_strategy, axes):
if pindex.distribution_strategy not in STRATEGIES['global']:
raise ValueError("ERROR: pindex's distribution strategy must be "
"global-type")
if pindex.distribution_strategy in STRATEGIES['slicing']:
if ((0 not in axes) or
(target_strategy is not pindex.distribution_strategy)):
raise ValueError(
"ERROR: A slicing distributor shall not be reshaped to "
"something non-sliced.")
semiscaled_shape = [1, ] * len(target_shape)
for i in axes:
semiscaled_shape[i] = target_shape[i]
local_data = pindex.get_local_data(copy=False)
semiscaled_local_data = local_data.reshape(semiscaled_shape)
result_obj = pindex.copy_empty(global_shape=target_shape,
distribution_strategy=target_strategy)
result_obj.set_full_data(semiscaled_local_data, copy=False)
return result_obj
......@@ -152,9 +152,13 @@ class RGSpace(Space):
self.paradict = rg_space_paradict(shape=shape,
zerocenter=zerocenter,
distances=distances)
distances=distances,
harmonic=harmonic)
self.dtype = np.dtype(dtype)
self.harmonic = bool(harmonic)
@property
def harmonic(self):
return self.paradict['harmonic']
def copy(self):
return RGSpace(dtype=self.dtype, harmonic=self.harmonic,
......
......@@ -203,8 +203,11 @@ class Space(object):
# parse dtype
dtype = np.dtype(dtype)
self.dtype = dtype
self._harmonic = None
self.harmonic = None
@property
def harmonic(self):
return self._harmonic
def __hash__(self):
# Extract the identifying parts from the vars(self) dict.
......@@ -212,21 +215,17 @@ class Space(object):
for (key, item) in vars(self).items():
if key in []:
continue
result_hash ^= item.__hash__() ^ 113*hash