Skip to content
Snippets Groups Projects
Commit 231eeda7 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

remaining cleanups

parent 09bba27a
No related branches found
No related tags found
1 merge request!68Martin's monster merge part 5/N: final cleanups, fixes and new testcases
Pipeline #
......@@ -19,7 +19,6 @@
import numpy as np
from itertools import product
def get_slice_list(shape, axes):
"""
Helper function which generates slice list(s) to traverse over all
......@@ -66,174 +65,6 @@ def get_slice_list(shape, axes):
return
def hermitianize_gaussian(x, axes=None):
# make the point inversions
flipped_x = _hermitianize_inverter(x, axes=axes)
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.
if axes is None:
axes = xrange(dimensions)
ndlist = [2 if i in axes else 1 for i in xrange(dimensions)]
ndlist = tuple(ndlist)
for i in np.ndindex(ndlist):
temp_index = tuple(i * mid_index)
x[temp_index] *= np.sqrt(0.5)
try:
x.hermitian = True
except(AttributeError):
pass
return x
def hermitianize(x, axes=None):
# make the point inversions
flipped_x = _hermitianize_inverter(x, axes=axes)
flipped_x = flipped_x.conjugate()
# average x and flipped_x.
# x = (x + flipped_x) / 2.
result_x = x + flipped_x
result_x /= 2.
# try:
# x.hermitian = True
# except(AttributeError):
# pass
return result_x
def _hermitianize_inverter(x, axes):
# 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()
if axes is None:
axes = xrange(dimensions)
# flip in the desired directions
for i in axes:
slice_picker = slice_primitive[:]
slice_picker[i] = slice(1, None, None)
slice_picker = tuple(slice_picker)
slice_inverter = slice_primitive[:]
slice_inverter[i] = slice(None, 0, -1)
slice_inverter = tuple(slice_inverter)
try:
y.set_data(to_key=slice_picker, data=y,
from_key=slice_inverter)
except(AttributeError):
y[slice_picker] = y[slice_inverter]
return y
def direct_vdot(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 not isinstance(current_extract, list) and \
not isinstance(current_extract, tuple):
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(reduce(lambda x, y: x * y, 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(reduce(lambda x, y: x * y, 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(reduce(lambda x, y: x * y, 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_axis_to_tuple(axis, length):
if axis is None:
......@@ -261,41 +92,6 @@ def cast_axis_to_tuple(axis, length):
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)
def get_default_codomain(domain):
from nifty.spaces import RGSpace, HPSpace, GLSpace, LMSpace
from nifty.operators.fft_operator.transformations import RGRGTransformation, \
HPLMTransformation, GLLMTransformation, LMGLTransformation
if isinstance(domain, RGSpace):
return RGRGTransformation.get_codomain(domain)
elif isinstance(domain, HPSpace):
return HPLMTransformation.get_codomain(domain)
elif isinstance(domain, GLSpace):
return GLLMTransformation.get_codomain(domain)
elif isinstance(domain, LMSpace):
# TODO: get the preferred transformation path from config
return LMGLTransformation.get_codomain(domain)
else:
raise TypeError('ERROR: unknown domain')
def parse_domain(domain):
from nifty.domain_object import DomainObject
if domain is None:
......@@ -308,6 +104,6 @@ def parse_domain(domain):
for d in domain:
if not isinstance(d, DomainObject):
raise TypeError(
"Given object contains something that is not a "
"Given object contains something that is not an "
"instance of DomainObject-class.")
return domain
......@@ -119,7 +119,7 @@ class PowerSpace(Space):
def get_distance_array(self, distribution_strategy):
result = d2o.distributed_data_object(
self.kindex,
self.kindex, dtype=np.float64,
distribution_strategy=distribution_strategy)
return result
......
......@@ -48,56 +48,18 @@ class RGSpace(Space):
NIFTY subclass for spaces of regular Cartesian grids.
Parameters
----------
num : {int, numpy.ndarray}
Number of gridpoints or numbers of gridpoints along each axis.
naxes : int, *optional*
Number of axes (default: None).
zerocenter : {bool, numpy.ndarray}, *optional*
Whether the Fourier zero-mode is located in the center of the grid
(or the center of each axis speparately) or not (default: True).
hermitian : bool, *optional*
Whether the fields living in the space follow hermitian symmetry or
not (default: True).
purelyreal : bool, *optional*
Whether the field values are purely real (default: True).
dist : {float, numpy.ndarray}, *optional*
Distance between two grid points along each axis (default: None).
fourier : bool, *optional*
Whether the space represents a Fourier or a position grid
(default: False).
Notes
-----
Only even numbers of grid points per axis are supported.
The basis transformations between position `x` and Fourier mode `k`
rely on (inverse) fast Fourier transformations using the
:math:`exp(2 \pi i k^\dagger x)`-formulation.
Attributes
----------
para : numpy.ndarray
One-dimensional array containing information on the axes of the
space in the following form: The first entries give the grid-points
along each axis in reverse order; the next entry is 0 if the
fields defined on the space are purely real-valued, 1 if they are
hermitian and complex, and 2 if they are not hermitian, but
complex-valued; the last entries hold the information on whether
the axes are centered on zero or not, containing a one for each
zero-centered axis and a zero for each other one, in reverse order.
dtype : numpy.dtype
Data type of the field values for a field defined on this space,
either ``numpy.float64`` or ``numpy.complex128``.
discrete : bool
Whether or not the underlying space is discrete, always ``False``
for regular grids.
vol : numpy.ndarray
One-dimensional array containing the distances between two grid
points along each axis, in reverse order. By default, the total
length of each axis is assumed to be one.
fourier : bool
harmonic : bool
Whether or not the grid represents a Fourier basis.
zerocenter : {bool, numpy.ndarray}, *optional*
Whether the Fourier zero-mode is located in the center of the grid
(or the center of each axis speparately) or not (default: True).
distances : {float, numpy.ndarray}, *optional*
Distance between two grid points along each axis (default: None).
"""
# ---Overwritten properties and methods---
......@@ -109,23 +71,16 @@ class RGSpace(Space):
Parameters
----------
num : {int, numpy.ndarray}
shape : {int, numpy.ndarray}
Number of gridpoints or numbers of gridpoints along each axis.
naxes : int, *optional*
Number of axes (default: None).
zerocenter : {bool, numpy.ndarray}, *optional*
Whether the Fourier zero-mode is located in the center of the
grid (or the center of each axis speparately) or not
(default: False).
hermitian : bool, *optional*
Whether the fields living in the space follow hermitian
symmetry or not (default: True).
purelyreal : bool, *optional*
Whether the field values are purely real (default: True).
dist : {float, numpy.ndarray}, *optional*
distances : {float, numpy.ndarray}, *optional*
Distance between two grid points along each axis
(default: None).
fourier : bool, *optional*
harmonic : bool, *optional*
Whether the space represents a Fourier or a position grid
(default: False).
......@@ -173,7 +128,7 @@ class RGSpace(Space):
hermitian_part = hermitian_part * np.sqrt(2)
anti_hermitian_part = anti_hermitian_part * np.sqrt(2)
# The fixed points of the point inversion must not be avaraged.
# The fixed points of the point inversion must not be averaged.
# Hence one must divide out the sqrt(2) again
# -> Get the middle index of the array
mid_index = np.array(hermitian_part.shape, dtype=np.int) // 2
......@@ -269,7 +224,7 @@ class RGSpace(Space):
shape = self.shape
# prepare the distributed_data_object
nkdict = distributed_data_object(
global_shape=shape,
global_shape=shape, dtype=np.float64,
distribution_strategy=distribution_strategy)
if distribution_strategy in DISTRIBUTION_STRATEGIES['slicing']:
......@@ -296,7 +251,7 @@ class RGSpace(Space):
cords = np.ogrid[inds]
dists = ((np.float128(0) + cords[0] - shape[0] // 2) * dk[0])**2
dists = ((np.float64(0) + cords[0] - shape[0] // 2) * dk[0])**2
# apply zerocenterQ shift
if not self.zerocenter[0]:
dists = np.fft.ifftshift(dists)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment