Commit 5b4de458 authored by Jait Dixit's avatar Jait Dixit
Browse files

Merge branch 'feature/field_multiple_space' of gitlab.mpcdf.mpg.de:ift/NIFTy...

Merge branch 'feature/field_multiple_space' of gitlab.mpcdf.mpg.de:ift/NIFTy into feature/field_multiple_space
parents 447d785c 36082174
......@@ -55,8 +55,7 @@ from field_types import FieldType,\
from operators import *
from power import PowerSpace,\
RGPowerSpace,\
LMPowerSpace
PowerIndexFactory
## optional submodule `rg`
try:
......
......@@ -43,40 +43,38 @@ class space_paradict(object):
temp_hash = hash(item)
except TypeError:
temp_hash = hash(tuple(item))
result_hash ^= temp_hash * hash(key)
result_hash ^= temp_hash ^ 131*hash(key)
return result_hash
class rg_space_paradict(space_paradict):
def __init__(self, shape, complexity, zerocenter):
self.ndim = len(np.array(shape).flatten())
def __init__(self, shape, zerocenter, distances):
if not hasattr(self, 'parameters'):
self.parameters = {}
self.parameters.__setitem__('shape', shape)
space_paradict.__init__(
self, shape=shape, complexity=complexity, zerocenter=zerocenter)
self, zerocenter=zerocenter, distances=distances)
def __setitem__(self, key, arg):
if key not in ['shape', 'complexity', 'zerocenter']:
if key not in ['shape', 'zerocenter', 'distances']:
raise ValueError(about._errors.cstring(
"ERROR: Unsupported RGSpace parameter:" + key))
if key == 'shape':
temp = np.array(arg, dtype=np.int).flatten()
if np.any(temp < 0):
raise ValueError("ERROR: negative number in shape.")
temp = tuple(temp)
if len(temp) != self.ndim:
raise ValueError(about._errors.cstring(
"ERROR: Number of dimensions does not match the init " +
"value."))
elif key == 'complexity':
temp = int(arg)
if temp not in [0, 1, 2]:
raise ValueError(about._errors.cstring(
"ERROR: Unsupported complexity parameter: " + str(temp)))
temp = np.empty(len(self['shape']), dtype=np.int)
temp[:] = arg
elif key == 'zerocenter':
temp = np.empty(self.ndim, dtype=bool)
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)
else:
temp = np.empty(len(self['shape']), dtype=np.float)
temp[:] = arg
temp = tuple(temp)
self.parameters.__setitem__(key, temp)
......@@ -187,82 +185,20 @@ class hp_space_paradict(space_paradict):
class power_space_paradict(space_paradict):
def __init__(self, distribution_strategy, log, nbin, binbounds):
def __init__(self, power_indices):
space_paradict.__init__(self,
distribution_strategy=distribution_strategy,
log=log,
nbin=nbin,
binbounds=binbounds)
power_indices=power_indices)
def __setitem__(self, key, arg):
if key not in ['distribution_strategy', 'log', 'nbin', 'binbounds']:
if key not in ['power_indices']:
raise ValueError(about._errors.cstring(
"ERROR: Unsupported PowerSpace parameter: " + key))
if key == 'log':
try:
temp = bool(arg)
except(TypeError):
temp = False
elif key == 'nbin':
try:
temp = int(arg)
except(TypeError):
temp = None
elif key == 'binbounds':
try:
temp = tuple(np.array(arg))
except(TypeError):
temp = None
elif key == 'distribution_strategy':
temp = str(arg)
if key == 'power_indices':
temp = arg
self.parameters.__setitem__(key, temp)
class rg_power_space_paradict(power_space_paradict, rg_space_paradict):
def __init__(self, shape, dgrid, zerocenter, distribution_strategy,
log, nbin, binbounds):
rg_space_paradict.__init__(self,
shape=shape,
complexity=0,
zerocenter=zerocenter)
power_space_paradict.__init__(
self,
distribution_strategy=distribution_strategy,
log=log,
nbin=nbin,
binbounds=binbounds)
self['dgrid'] = dgrid
def __setitem__(self, key, arg):
if key not in ['shape', 'complexity', 'zerocenter',
'distribution_strategy', 'log', 'nbin', 'binbounds',
'dgrid']:
raise ValueError(about._errors.cstring(
"ERROR: Unsupported RGPowerSpace parameter: " + key))
if key in ['distribution_strategy', 'log', 'nbin', 'binbounds']:
power_space_paradict.__setitem__(self, key, arg)
elif key == 'dgrid':
temp = np.array(arg, dtype=np.float).flatten()
if np.size(temp) != self.ndim:
temp = temp * np.ones(self.ndim, dtype=np.float)
temp = tuple(temp)
self.parameters.__setitem__(key, temp)
else:
rg_space_paradict.__setitem__(self, key, arg)
def __hash__(self):
return (self.parameters['power_indices']['hash'] ^
131*hash('power_indices'))
......@@ -236,4 +236,21 @@ def cast_axis_to_tuple(axis, length):
# shift negative indices to positive ones
axis = tuple(item if (item >= 0) else (item + length) for item in axis)
return axis
\ No newline at end of file
return axis
def complex_bincount(x, weights=None, minlength=None):
try:
complex_weights_Q = issubclass(weights.dtype.type,
np.complexfloating)
except AttributeError:
complex_weights_Q = False
if complex_weights_Q:
real_bincount = x.bincount(weights=weights.real,
minlength=minlength)
imag_bincount = x.bincount(weights=weights.imag,
minlength=minlength)
return real_bincount + imag_bincount
else:
return x.bincount(weights=weights, minlength=minlength)
# -*- coding: utf-8 -*-
from power_space import PowerSpace
from rg_power_space import RGPowerSpace
from lm_power_space import LMPowerSpace
from power_index_factory import PowerIndexFactory,\
RGPowerIndexFactory,\
LMPowerIndexFactory
\ No newline at end of file
from power_index_factory import PowerIndexFactory
\ No newline at end of file
# -*- coding: utf-8 -*-
from power_indices import PowerIndices,\
RGPowerIndices,\
LMPowerIndices
from power_indices import PowerIndices
class _PowerIndexFactory(object):
def __init__(self):
self.power_indices_storage = {}
def _get_power_index_class(self):
return PowerIndices
def get_power_indices(self, log, nbin, binbounds,
domain, distribution_strategy):
current_hash = domain.__hash__() ^ (111*hash(distribution_strategy))
def hash_arguments(self, **kwargs):
return frozenset(kwargs.items())
def get_power_indices(self, log, nbin, binbounds, **kwargs):
current_hash = self.hash_arguments(**kwargs)
if current_hash not in self.power_indices_storage:
power_class = self._get_power_index_class()
self.power_indices_storage[current_hash] = power_class(
log=log,
nbin=nbin,
binbounds=binbounds,
**kwargs)
self.power_indices_storage[current_hash] = \
PowerIndices(domain, distribution_strategy,
log=log, nbin=nbin, binbounds=binbounds)
power_indices = self.power_indices_storage[current_hash]
power_index = power_indices.get_index_dict(log=log,
nbin=nbin,
......@@ -30,15 +21,4 @@ class _PowerIndexFactory(object):
return power_index
class _RGPowerIndexFactory(_PowerIndexFactory):
def _get_power_index_class(self):
return RGPowerIndices
class _LMPowerIndexFactory(_PowerIndexFactory):
def _get_power_index_class(self):
return LMPowerIndices
PowerIndexFactory = _PowerIndexFactory()
RGPowerIndexFactory = _RGPowerIndexFactory()
LMPowerIndexFactory = _LMPowerIndexFactory()
......@@ -13,8 +13,8 @@ hp = gdi.get('healpy')
class PowerIndices(object):
def __init__(self, distribution_strategy, log=False, nbin=None,
binbounds=None):
def __init__(self, domain, 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
......@@ -41,12 +41,11 @@ class PowerIndices(object):
Array-like inner boundaries of the used bins of the default
indices.
"""
self._comm = getattr(gdi[gc['mpi_module']], gc['default_comm'])
self.domain = domain
self.distribution_strategy = distribution_strategy
# Compute the global kdict
self.kdict = self.compute_kdict()
self.kdict = self.domain.compute_k_array(distribution_strategy)
# Initialize the dictonary which stores all individual index-dicts
self.global_dict = {}
# Set self.default_parameters
......@@ -228,39 +227,61 @@ class PowerIndices(object):
temp_kdict = self._compute_kdict_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}
"kdict": temp_kdict,
"domain": self.domain,
"hash": index_dict_hash}
return temp_index_dict
def _compute_kdict_from_pindex_kindex(self, pindex, kindex):
if isinstance(pindex, distributed_data_object):
tempindex = pindex.copy(dtype=kindex.dtype)
result = tempindex.apply_scalar_function(
lambda x: kindex[x.astype(np.dtype('int'))])
else:
result = kindex[pindex].astype(dtype=kindex.dtype)
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):
if self.distribution_strategy in ['np']:
return self._compute_indices_np(nkdict)
else:
return self._compute_indices_d2o(nkdict)
def _compute_indices_d2o(self, nkdict):
"""
Internal helper function which computes pindex, kindex, rho and pundex
from a given nkdict
"""
raise NotImplementedError(
about._errors.cstring(
"ERROR: No generic _compute_indices_d2o method implemented."))
##########
# kindex #
##########
global_kindex = nkdict.unique()
def _compute_pundex_d2o(self, global_pindex, global_kindex):
##########
# pindex #
##########
# compute the local pindex slice on basis of the local nkdict data
local_pindex = np.searchsorted(global_kindex, nkdict.get_local_data())
# prepare the distributed_data_object
global_pindex = distributed_data_object(
global_shape=nkdict.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
......@@ -291,7 +312,9 @@ class PowerIndices(object):
# Store the individual pundices in the local_pundex array
local_pundex[temp_uniqued_pindex] = local_temp_pundex
# Use Allreduce to find the first occurences/smallest pundices
self._comm.Allreduce(local_pundex, global_pundex, op=MPI.MIN)
global_pindex.comm.Allreduce(local_pundex,
global_pundex,
op=MPI.MIN)
return global_pundex
elif self.distribution_strategy in DISTRIBUTION_STRATEGIES['not']:
......@@ -306,28 +329,6 @@ class PowerIndices(object):
"ERROR: _compute_pundex_d2o not implemented for given "
"distribution_strategy"))
def _compute_indices_np(self, nkdict):
"""
Internal helper function which computes pindex, kindex, rho and pundex
from a given nkdict
"""
raise NotImplementedError(
about._errors.cstring(
"ERROR: No generic _compute_indices_np method implemented."))
def _compute_pundex_np(self, pindex, 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.
"""
##########
# pundex #
##########
pundex = np.unique(pindex, return_index=True)[1]
return pundex
def _bin_power_indices(self, index_dict, **kwargs):
"""
Returns the binned power indices associated with the Fourier grid.
......@@ -404,203 +405,13 @@ class PowerIndices(object):
(rho_[reorder[ii]] + rho[ii]))
rho_[reorder[ii]] += rho[ii]
if self.distribution_strategy == 'np':
pindex_ = reorder[pindex]
pundex_ = self._compute_pundex_np(pindex_, kindex_)
else:
pindex_ = pindex.copy_empty()
pindex_.set_local_data(reorder[pindex.get_local_data()])
pundex_ = self._compute_pundex_d2o(pindex_, kindex_)
pindex_ = pindex.copy_empty()
pindex_.set_local_data(reorder[pindex.get_local_data()])
pundex_ = self._compute_pundex(pindex_, kindex_)
return pindex_, kindex_, rho_, pundex_
class RGPowerIndices(PowerIndices):
def __init__(self, shape, dgrid, distribution_strategy, zerocenter=False,
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.shape = np.array(shape, dtype=int)
self.dgrid = np.abs(np.array(dgrid))
if self.shape.shape != self.dgrid.shape:
raise ValueError(about._errors.cstring("ERROR: The supplied shape\
and dgrid have not the same dimensionality"))
self.zerocenter = self._cast_zerocenter(zerocenter)
super(RGPowerIndices, self).__init__(
distribution_strategy=distribution_strategy,
log=log,
nbin=nbin,
binbounds=binbounds)
def _cast_zerocenter(self, zerocenter=False):
"""
internal helper function which brings the zerocenter input in
the form of a boolean-tuple
"""
zc = np.array(zerocenter).astype(bool)
if zc.shape == self.shape.shape:
return tuple(zc)
else:
temp = np.empty(shape=self.shape.shape, dtype=bool)
temp[:] = zc
return tuple(temp)
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
"""
shape = self.shape
if self.distribution_strategy == 'np':
slice_of_first_dimension = slice(0, shape[0])
nkdict = self._compute_kdict_helper(slice_of_first_dimension)
else:
# prepare the distributed_data_object
nkdict = distributed_data_object(
global_shape=shape,
dtype=np.float128,
distribution_strategy=self.distribution_strategy,
comm=self._comm)
if self.distribution_strategy in \
DISTRIBUTION_STRATEGIES['slicing']:
# get the node's individual slice of the first dimension
slice_of_first_dimension = slice(
*nkdict.distributor.local_slice[0:2])
elif self.distribution_strategy in DISTRIBUTION_STRATEGIES['not']:
slice_of_first_dimension = slice(0, shape[0])
else:
raise ValueError(about._errors.cstring(
"ERROR: Unsupported distribution strategy"))
dists = self._compute_kdict_helper(slice_of_first_dimension)
nkdict.set_local_data(dists)
return nkdict
def _compute_kdict_helper(self, slice_of_first_dimension):
dk = self.dgrid
shape = self.shape
inds = []
for a in shape:
inds += [slice(0, a)]
cords = np.ogrid[inds]
dists = ((np.float128(0) + cords[0] - shape[0] // 2) * dk[0])**2
# apply zerocenterQ shift
if self.zerocenter[0] == False:
dists = np.fft.fftshift(dists)
# only save the individual slice
dists = dists[slice_of_first_dimension]
for ii in range(1, len(shape)):
temp = ((cords[ii] - shape[ii] // 2) * dk[ii])**2
if self.zerocenter[ii] == False:
temp = np.fft.fftshift(temp)
dists = dists + temp
dists = np.sqrt(dists)
return dists
def _compute_indices_d2o(self, nkdict):
"""
Internal helper function which computes pindex, kindex, rho and pundex
from a given nkdict
"""
##########
# kindex #
##########
global_kindex = nkdict.unique()
##########
# pindex #
##########
# compute the local pindex slice on basis of the local nkdict data
local_pindex = np.searchsorted(global_kindex, nkdict.get_local_data())
# prepare the distributed_data_object
global_pindex = distributed_data_object(
global_shape=nkdict.shape,
dtype=local_pindex.dtype,
distribution_strategy=self.distribution_strategy,
comm=self._comm)
# 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_d2o(global_pindex,
global_kindex)
return global_pindex, global_kindex, global_rho, global_pundex
def _compute_indices_np(self, nkdict):
"""
Internal helper function which computes pindex, kindex, rho and pundex
from a given nkdict
"""
##########
# kindex #
##########
kindex = np.unique(nkdict)
##########
# pindex #
##########
pindex = np.searchsorted(kindex, nkdict)
#######
# rho #
#######
rho = np.bincount(pindex.flatten())
##########
# pundex #
##########
pundex = self._compute_pundex_np(pindex, kindex)
return pindex, kindex, rho, pundex
class LMPowerIndices(PowerIndices):
def __init__(self, lmax, dim, distribution_strategy,
......
# -*- coding: utf-8 -*-
import numpy as np
from d2o import STRATEGIES
from nifty.space import Space
from nifty.nifty_paradict import power_space_paradict
class PowerSpace(Space):
def __init__(self, distribution_strategy, dtype=np.dtype('float'),
log=False, nbin=None, binbounds=None):
def __init__(self, power_indices, dtype=np.dtype('float')):
self.dtype = np.dtype(dtype)
self.paradict = power_space_paradict(
distribution_strategy=distribution_strategy,
log=log,