Commit 0676944f authored by Ultima's avatar Ultima
Browse files

Added nifty_utilities.py.

Started rework of the d2o slicing: introduced 0-dimensional advanced indexing capabilities.
parent 23e0643e
...@@ -35,6 +35,7 @@ from nifty_mpi_data import distributed_data_object ...@@ -35,6 +35,7 @@ from nifty_mpi_data import distributed_data_object
from nifty_power import * from nifty_power import *
from nifty_random import random from nifty_random import random
from nifty_simple_math import * from nifty_simple_math import *
from nifty_utilities import *
from nifty_paradict import space_paradict,\ from nifty_paradict import space_paradict,\
point_space_paradict,\ point_space_paradict,\
......
...@@ -170,7 +170,9 @@ class notification(switch): ...@@ -170,7 +170,9 @@ class notification(switch):
""" """
_code = "\033[0m" ## "\033[39;49m" _code = "\033[0m" ## "\033[39;49m"
_ccode_red = "\033[31;1m"
_ccode_yellow = "\033[33;1m"
_ccode_green = "\033[32;1m"
def __init__(self,default=True,ccode="\033[0m"): def __init__(self,default=True,ccode="\033[0m"):
""" """
Initializes the notification and sets `status` and `ccode` Initializes the notification and sets `status` and `ccode`
...@@ -387,9 +389,12 @@ class _about(object): ## nifty support class for global settings ...@@ -387,9 +389,12 @@ class _about(object): ## nifty support class for global settings
self._version = str(__version__) self._version = str(__version__)
## switches and notifications ## switches and notifications
self._errors = notification(default=True,ccode=notification._code) self._errors = notification(default=True,
self.warnings = notification(default=True,ccode=notification._code) ccode=notification._code)
self.infos = notification(default=False,ccode=notification._code) self.warnings = notification(default=True,
ccode=notification._code)
self.infos = notification(default=False,
ccode=notification._code)
self.multiprocessing = switch(default=True) self.multiprocessing = switch(default=True)
self.hermitianize = switch(default=True) self.hermitianize = switch(default=True)
self.lm2gl = switch(default=True) self.lm2gl = switch(default=True)
......
## nifty configuration ## nifty configuration
## ##
## errors colour code ## errors colour code
0 31
## warnings ## warnings
1 1
## warnings colour code ## warnings colour code
0 33;1
## infos ## infos
0 0
## infos colour code ## infos colour code
0 32
## multiprocessing ## multiprocessing
1 1
## hermitianize ## hermitianize
......
...@@ -152,6 +152,7 @@ from nifty_about import about ...@@ -152,6 +152,7 @@ from nifty_about import about
from nifty_random import random from nifty_random import random
from nifty.nifty_mpi_data import distributed_data_object from nifty.nifty_mpi_data import distributed_data_object
import nifty.nifty_utilities as utilities
pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
...@@ -3046,32 +3047,7 @@ class field(object): ...@@ -3046,32 +3047,7 @@ class field(object):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _map(self, function, *args): def _map(self, function, *args):
if self.ishape == (): return utilities.field_map(self.ishape, function, *args)
return function(*args)
else:
if args == ():
result = np.empty(self.ishape, dtype=np.object)
for i in xrange(np.prod(self.ishape)):
ii = np.unravel_index(i, self.ishape)
result[ii] = function()
return result
else:
## define a helper function in order to clip the get-indices
## to be suitable for the foreign arrays in args.
## This allows you to do operations, like adding to fields
## with ishape (3,4,3) and (3,4,1)
def get_clipped(w, ind):
w_shape = np.array(np.shape(w))
get_tuple = tuple(np.clip(ind, 0, w_shape-1))
return w[get_tuple]
result = np.empty_like(args[0])
for i in xrange(np.prod(result.shape)):
ii = np.unravel_index(i, result.shape)
result[ii] = function(*map(
lambda z: get_clipped(z, ii), args)
)
#result[ii] = function(*map(lambda z: z[ii], args))
return result
def cast(self, x = None, ishape = None): def cast(self, x = None, ishape = None):
if ishape is None: if ishape is None:
......
This diff is collapsed.
...@@ -483,72 +483,7 @@ def conjugate(x): ...@@ -483,72 +483,7 @@ def conjugate(x):
""" """
return _math_helper(x, np.conjugate) return _math_helper(x, np.conjugate)
def direct_dot(x, y):
## the input could be fields. Try to extract the data
try:
x = x.get_val()
except(AttributeError):
pass
try:
y = y.get_val()
except(AttributeError):
pass
## try to make a direct vdot
try:
return x.vdot(y)
except(AttributeError):
pass
try:
return y.vdot(x)
except(AttributeError):
pass
## fallback to numpy
return np.vdot(x, y)
def convert_nested_list_to_object_array(x):
## if x is a nested_list full of ndarrays all having the same size,
## np.shape returns the shape of the ndarrays, too, i.e. too many
## dimensions
possible_shape = np.shape(x)
## Check if possible_shape goes too deep.
dimension_counter = 0
current_extract = x
for i in xrange(len(possible_shape)):
if isinstance(current_extract, list) == False and\
isinstance(current_extract, tuple) == False:
break
current_extract = current_extract[0]
dimension_counter += 1
real_shape = possible_shape[:dimension_counter]
## if the numpy array was not encapsulated at all, return x directly
if real_shape == ():
return x
## Prepare the carrier-object
carrier = np.empty(real_shape, dtype = np.object)
for i in xrange(np.prod(real_shape)):
ii = np.unravel_index(i, real_shape)
try:
carrier[ii] = x[ii]
except(TypeError):
extracted = x
for j in xrange(len(ii)):
extracted = extracted[ii[j]]
carrier[ii] = extracted
return carrier
......
# -*- coding: utf-8 -*-
import numpy as np
def hermitianize(x):
## make the point inversions
flipped_x = _hermitianize_inverter(x)
flipped_x = flipped_x.conjugate()
## check if x was already hermitian
if (x == flipped_x).all():
return x
## average x and flipped_x.
## Correct the variance by multiplying sqrt(0.5)
x = (x + flipped_x) * np.sqrt(0.5)
## The fixed points of the point inversion must not be avaraged.
## Hence one must multiply them again with sqrt(0.5)
## -> Get the middle index of the array
mid_index = np.array(x.shape, dtype=np.int)//2
dimensions = mid_index.size
## Use ndindex to iterate over all combinations of zeros and the
## mid_index in order to correct all fixed points.
for i in np.ndindex((2,)*dimensions):
temp_index = tuple(i*mid_index)
x[temp_index] *= np.sqrt(0.5)
try:
x.hermitian = True
except(AttributeError):
pass
return x
def _hermitianize_inverter(x):
## calculate the number of dimensions the input array has
dimensions = len(x.shape)
## prepare the slicing object which will be used for mirroring
slice_primitive = [slice(None),]*dimensions
## copy the input data
y = x.copy()
## flip in every direction
for i in xrange(dimensions):
slice_picker = slice_primitive[:]
slice_picker[i] = slice(1, None,None)
slice_inverter = slice_primitive[:]
slice_inverter[i] = slice(None, 0, -1)
try:
y.inject(to_slices=slice_picker, data=y,
from_slices=slice_inverter)
except(AttributeError):
y[slice_picker] = y[slice_inverter]
return y
def direct_dot(x, y):
## the input could be fields. Try to extract the data
try:
x = x.get_val()
except(AttributeError):
pass
try:
y = y.get_val()
except(AttributeError):
pass
## try to make a direct vdot
try:
return x.vdot(y)
except(AttributeError):
pass
try:
return y.vdot(x)
except(AttributeError):
pass
## fallback to numpy
return np.vdot(x, y)
def convert_nested_list_to_object_array(x):
## if x is a nested_list full of ndarrays all having the same size,
## np.shape returns the shape of the ndarrays, too, i.e. too many
## dimensions
possible_shape = np.shape(x)
## Check if possible_shape goes too deep.
dimension_counter = 0
current_extract = x
for i in xrange(len(possible_shape)):
if isinstance(current_extract, list) == False and\
isinstance(current_extract, tuple) == False:
break
current_extract = current_extract[0]
dimension_counter += 1
real_shape = possible_shape[:dimension_counter]
## if the numpy array was not encapsulated at all, return x directly
if real_shape == ():
return x
## Prepare the carrier-object
carrier = np.empty(real_shape, dtype = np.object)
for i in xrange(np.prod(real_shape)):
ii = np.unravel_index(i, real_shape)
try:
carrier[ii] = x[ii]
except(TypeError):
extracted = x
for j in xrange(len(ii)):
extracted = extracted[ii[j]]
carrier[ii] = extracted
return carrier
def field_map(ishape, function, *args):
if ishape == ():
return function(*args)
else:
if args == ():
result = np.empty(ishape, dtype=np.object)
for i in xrange(np.prod(ishape)):
ii = np.unravel_index(i, ishape)
result[ii] = function()
return result
else:
## define a helper function in order to clip the get-indices
## to be suitable for the foreign arrays in args.
## This allows you to do operations, like adding to fields
## with ishape (3,4,3) and (3,4,1)
def get_clipped(w, ind):
w_shape = np.array(np.shape(w))
get_tuple = tuple(np.clip(ind, 0, w_shape-1))
return w[get_tuple]
result = np.empty_like(args[0])
for i in xrange(np.prod(result.shape)):
ii = np.unravel_index(i, result.shape)
result[ii] = function(*map(
lambda z: get_clipped(z, ii), args)
)
#result[ii] = function(*map(lambda z: z[ii], args))
return result
\ No newline at end of file
...@@ -31,7 +31,7 @@ from nifty_probing import trace_prober,\ ...@@ -31,7 +31,7 @@ from nifty_probing import trace_prober,\
inverse_trace_prober,\ inverse_trace_prober,\
diagonal_prober,\ diagonal_prober,\
inverse_diagonal_prober inverse_diagonal_prober
import nifty.nifty_utilities as utilities
import nifty_simple_math import nifty_simple_math
...@@ -2159,17 +2159,26 @@ class power_operator(diagonal_operator): ...@@ -2159,17 +2159,26 @@ class power_operator(diagonal_operator):
self.val = diag self.val = diag
""" """
@property
def val(self):
return self._val
## The domain is used for calculations of the power-spectrum, not for ## The domain is used for calculations of the power-spectrum, not for
## actual field values. Therefore the casting of self.val must be switched ## actual field values. Therefore the casting of self.val must be switched
## off. ## off.
@val.setter def set_val(self, new_val):
def val(self, x): """
self._val = x Resets the field values.
Parameters
----------
new_val : {scalar, ndarray}
New field values either as a constant or an arbitrary array.
"""
self.val = new_val
return self.val
def get_val(self):
return self.val
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def set_power(self, new_spec, bare=True, pindex=None, **kwargs): def set_power(self, new_spec, bare=True, pindex=None, **kwargs):
...@@ -2374,7 +2383,9 @@ class power_operator(diagonal_operator): ...@@ -2374,7 +2383,9 @@ class power_operator(diagonal_operator):
""" """
pindex = self._cast_pindex(pindex, **kwargs) pindex = self._cast_pindex(pindex, **kwargs)
return projection_operator(self.domain, assign=pindex) return projection_operator(self.domain,
codomain = self.codomain,
assign=pindex)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
...@@ -2668,7 +2679,7 @@ class projection_operator(operator): ...@@ -2668,7 +2679,7 @@ class projection_operator(operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#@profile
def pseudo_tr(self, x, axis=(), **kwargs): def pseudo_tr(self, x, axis=(), **kwargs):
""" """
Computes the pseudo trace of a given object for all projection bands Computes the pseudo trace of a given object for all projection bands
...@@ -2739,13 +2750,16 @@ class projection_operator(operator): ...@@ -2739,13 +2750,16 @@ class projection_operator(operator):
self.domain.get_meta_volume(total=False), self.domain.get_meta_volume(total=False),
power = -1) power = -1)
## prepare the result object ## prepare the result object
big_shape = working_field.ishape + (len(self.ind),) #big_shape = working_field.ishape + (len(self.ind),)
projection_result = np.empty(big_shape, dtype = np.object_) # working_data = working_field.get_val()
for i in xrange(len(self.ind)): #projection_result = np.empty(big_shape,
temp_projected = self(working_field, bands = self.ind[i]) # dtype = working_field.domain.datatype)
temp_dotted = temp_projected.dot(1, axis=(), bare=True)
projection_result[(slice(None),)*len(working_field.ishape) + (i,)]\ projection_result = utilities.field_map(
= temp_dotted working_field.ishape,
lambda z: self.domain.calc_bincount(self.assign, weights = z),
working_field.get_val())
projection_result = np.sum(projection_result, axis=axis) projection_result = np.sum(projection_result, axis=axis)
return projection_result return projection_result
......
...@@ -24,7 +24,7 @@ from __future__ import division ...@@ -24,7 +24,7 @@ from __future__ import division
from nifty.nifty_about import about from nifty.nifty_about import about
from nifty.nifty_core import space, \ from nifty.nifty_core import space, \
field field
from nifty.nifty_simple_math import direct_dot from nifty.nifty_utilities import direct_dot
##============================================================================= ##=============================================================================
......
...@@ -48,7 +48,7 @@ from nifty.nifty_core import point_space, \ ...@@ -48,7 +48,7 @@ from nifty.nifty_core import point_space, \
from nifty.nifty_random import random from nifty.nifty_random import random
from nifty.nifty_mpi_data import distributed_data_object from nifty.nifty_mpi_data import distributed_data_object
from nifty.nifty_paradict import rg_space_paradict from nifty.nifty_paradict import rg_space_paradict
import nifty.nifty_utilities as utilities
import fft_rg import fft_rg
''' '''
...@@ -2383,59 +2383,3 @@ class power_indices(object): ...@@ -2383,59 +2383,3 @@ class power_indices(object):
pundex_ = self.__compute_pundex__(pindex_, kindex_) pundex_ = self.__compute_pundex__(pindex_, kindex_)
return pindex_, kindex_, rho_, pundex_ return pindex_, kindex_, rho_, pundex_
class utilities(object):
def __init__(self):
pass
@staticmethod
def hermitianize(x):
## make the point inversions
flipped_x = utilities._hermitianize_inverter(x)
flipped_x = flipped_x.conjugate()
## check if x was already hermitian
if (x == flipped_x).all():
return x
## average x and flipped_x.
## Correct the variance by multiplying sqrt(0.5)
x = (x + flipped_x) * np.sqrt(0.5)
## The fixed points of the point inversion must not be avaraged.
## Hence one must multiply them again with sqrt(0.5)
## -> Get the middle index of the array
mid_index = np.array(x.shape, dtype=np.int)//2
dimensions = mid_index.size
## Use ndindex to iterate over all combinations of zeros and the
## mid_index in order to correct all fixed points.
for i in np.ndindex((2,)*dimensions):
temp_index = tuple(i*mid_index)
x[temp_index] *= np.sqrt(0.5)
try:
x.hermitian = True
except(AttributeError):
pass
return x
@staticmethod
def _hermitianize_inverter(x):
## calculate the number of dimensions the input array has
dimensions = len(x.shape)
## prepare the slicing object which will be used for mirroring
slice_primitive = [slice(None),]*dimensions
## copy the input data
y = x.copy()
## flip in every direction
for i in xrange(dimensions):
slice_picker = slice_primitive[:]
slice_picker[i] = slice(1, None,None)
slice_inverter = slice_primitive[:]
slice_inverter[i] = slice(None, 0, -1)
try:
y.inject(to_slices=slice_picker, data=y,
from_slices=slice_inverter)
except(AttributeError):
y[slice_picker] = y[slice_inverter]
return y
...@@ -19,9 +19,9 @@ except(ImportError): ...@@ -19,9 +19,9 @@ except(ImportError):
found['h5py'] = False found['h5py'] = False
all_datatypes = [np.bool_, np.int16, np.uint16, np.uint32, np.int32, np.int, all_datatypes = [np.bool_, np.int16, np.uint16, np.uint32, np.int32, np.int_,
np.int64, np.uint64, np.int64, np.uint64, np.float32, np.float, np.int64, np.uint64, np.int64, np.uint64, np.float32, np.float_,
np.float64, np.float128, np.complex64, np.complex, np.complex128] np.float64, np.float128, np.complex64, np.complex_, np.complex128]
all_distribution_strategies = ['not', 'equal', 'fftw'] all_distribution_strategies = ['not', 'equal', 'fftw']
binary_non_inplace_operators = ['__add__', '__radd__', '__sub__', '__rsub__', binary_non_inplace_operators = ['__add__', '__radd__', '__sub__', '__rsub__',
...@@ -110,7 +110,7 @@ class Test_Initialization(unittest.TestCase): ...@@ -110,7 +110,7 @@ class Test_Initialization(unittest.TestCase):
@parameterized.expand(itertools.product([ @parameterized.expand(itertools.product([
[1, (13,7), np.float64, (13,7), np.float64], [1, (13,7), np.float64, (13,7), np.float64],
[np.array([1]), (13,7), np.float64, (1,), np.float64], [np.array([1]), (13,7), np.float64, (1,), np.float64],
[np.array([[1.,2.],[3.,4.]]), (13,7), np.int, (2,2), np.int], [np.array([[1.,2.],[3.,4.]]), (13,7), np.int_, (2,2), np.int_],
], ['not', 'equal', 'fftw']), ], ['not', 'equal', 'fftw']),
testcase_func_name=custom_name_func) testcase_func_name=custom_name_func)
def test_special_init_cases(self, def test_special_init_cases(self,
...@@ -298,7 +298,7 @@ def scalar_only_square(x): ...@@ -298,7 +298,7 @@ def scalar_only_square(x):
class Test_copy_and_copy_empty(unittest.TestCase): class Test_copy_and_copy_empty(unittest.TestCase):
@parameterized.expand([ @parameterized.expand([
((8,7), None, None), ((8,7), None, None),
(None, np.float, None), (None, np.float_, None),
(None, None, 'not') (None, None, 'not')
], testcase_func_name=custom_name_func) ], testcase_func_name=custom_name_func)
def test_copy_empty(self, def test_copy_empty(self,
...@@ -322,7 +322,7 @@ class Test_copy_and_copy_empty(unittest.TestCase): ...@@ -322,7 +322,7 @@ class Test_copy_and_copy_empty(unittest.TestCase):
assert_equal(p.distribution_strategy, new_distribution_strategy) assert_equal(p.distribution_strategy, new_distribution_strategy)
@parameterized.expand([ @parameterized.expand([
(np.float, None), (np.float_, None),
(None, 'not') (None, 'not')
], testcase_func_name=custom_name_func) ], testcase_func_name=custom_name_func)
def test_copy(self, def test_copy(self,
......
Markdown is supported
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