Commit 54885033 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'mmm5' into 'master'

Martin's monster merge part 5/N: final cleanups, fixes and new testcases

See merge request !68
parents fe6bf4e5 99b6031e
Pipeline #11742 failed with stages
in 10 minutes and 7 seconds
......@@ -99,7 +99,7 @@ Starting with a fresh Ubuntu installation move to a folder like
- Install pyHealpix:
git clone https://gitlab.mpcdf.mpg.de/mtr/pyHealpix.git
git clone https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
cd pyHealpix
autoreconf -i && ./configure && sudo make install
cd ..
......@@ -133,7 +133,7 @@ your local user directory.
For pyHealpix, use:
git clone https://gitlab.mpcdf.mpg.de/mtr/pyHealpix.git
git clone https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
cd pyHealpix
autoreconf -i && ./configure --prefix=$HOME/.local && make install
cd ..
......@@ -154,7 +154,7 @@ may cause trouble.
- Install pyHealpix:
git clone https://gitlab.mpcdf.mpg.de/mtr/pyHealpix.git
git clone https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git
cd pyHealpix
autoreconf -i && ./configure && sudo make install
cd ..
......
#!/bin/bash
git clone -b mpi https://github.com/fredRos/pyFFTW.git
git clone -b mpi https://github.com/ultimanet/pyFFTW.git
cd pyFFTW/
CC=mpicc python setup.py build_ext install
cd ..
......
......@@ -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
......@@ -83,16 +83,13 @@ class HPLMTransformation(SlicingTransformation):
super(HPLMTransformation, cls).check_codomain(domain, codomain)
def _transformation_of_slice(self, inp, **kwargs):
nside = self.domain.nside
lmax = self.codomain.lmax
mmax = lmax
sjob = pyHealpix.sharpjob_d()
sjob.set_Healpix_geometry(nside)
sjob.set_triangular_alm_info(lmax, mmax)
if issubclass(inp.dtype.type, np.complexfloating):
[resultReal, resultImag] = [sjob.map2alm(x)
for x in (inp.real, inp.imag)]
[resultReal,
resultImag] = [pyHealpix.map2alm_iter(x, lmax, mmax, 3)
for x in (inp.real, inp.imag)]
[resultReal,
resultImag] = [lm_transformation_factory.buildIdx(x, lmax=lmax)
......@@ -101,7 +98,7 @@ class HPLMTransformation(SlicingTransformation):
result = self._combine_complex_result(resultReal, resultImag)
else:
result = sjob.map2alm(inp)
result = pyHealpix.map2alm_iter(inp, lmax, mmax, 3)
result = lm_transformation_factory.buildIdx(result, lmax=lmax)
return result
......@@ -64,7 +64,7 @@ class LMHPTransformation(SlicingTransformation):
if not isinstance(domain, LMSpace):
raise TypeError("domain needs to be a LMSpace.")
nside = np.max(domain.lmax//2, 1)
nside = max((domain.lmax + 1)//2, 1)
result = HPSpace(nside=nside, dtype=domain.dtype)
return result
......@@ -89,21 +89,18 @@ class LMHPTransformation(SlicingTransformation):
lmax = self.domain.lmax
mmax = lmax
sjob = pyHealpix.sharpjob_d()
sjob.set_Healpix_geometry(nside)
sjob.set_triangular_alm_info(lmax, mmax)
if issubclass(inp.dtype.type, np.complexfloating):
[resultReal,
resultImag] = [lm_transformation_factory.buildLm(x, lmax=lmax)
for x in (inp.real, inp.imag)]
[resultReal, resultImag] = [sjob.alm2map(x)
[resultReal, resultImag] = [pyHealpix.alm2map(x, lmax, mmax, nside)
for x in [resultReal, resultImag]]
result = self._combine_complex_result(resultReal, resultImag)
else:
result = lm_transformation_factory.buildLm(inp, lmax=lmax)
result = sjob.alm2map(result)
result = pyHealpix.alm2map(result, lmax, mmax, nside)
return result
......@@ -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
......
......@@ -50,56 +50,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---
......@@ -111,23 +73,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).
......@@ -175,7 +130,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
......@@ -271,7 +226,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']:
......@@ -298,7 +253,7 @@ class RGSpace(Space):
cords = np.ogrid[inds]
dists = ((np.float128(0) + cords[0] - shape[0] // 2) * dk[0])**2
dists = ((cords[0] - shape[0]//2)*dk[0])**2
# apply zerocenterQ shift
if not self.zerocenter[0]:
dists = np.fft.ifftshift(dists)
......
# NIFTy
# Copyright (C) 2017 Theo Steininger
#
# Author: Theo Steininger
#
# 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/>.
import unittest
import numpy as np
from numpy.testing import assert_,\
assert_equal, \
assert_allclose
from nifty.config import dependency_injector as di
from nifty import Field,\
RGSpace,\
LMSpace,\
GLSpace,\
FieldArray, \
RGRGTransformation, \
LMGLTransformation, \
LMHPTransformation, \
FFTOperator
from nose.plugins.skip import SkipTest
def _harmonic_type(itp):
otp=itp
if otp==np.float64:
otp=np.complex128
elif otp==np.float32:
otp=np.complex64
return otp
class Misc_Tests(unittest.TestCase):
def test_RG_distance_1D(self):
for dim1 in [10,11]:
for zc1 in [False,True]:
for d in [0.1,1,3.7]:
foo = RGSpace([dim1],zerocenter=zc1)
res = foo.get_distance_array('not')
assert_equal(res[zc1*(dim1//2)],0.)
def test_RG_distance_2D(self):
for dim1 in [10,11]:
for dim2 in [9,28]:
for zc1 in [False,True]:
for zc2 in [False,True]:
for d in [0.1,1,3.7]:
foo = RGSpace([dim1,dim2],zerocenter=[zc1,zc2])
res = foo.get_distance_array('not')
assert_equal(res[zc1*(dim1//2),zc2*(dim2//2)],0.)
def test_fft1D(self):
for dim1 in [10,11]:
for zc1 in [False,True]:
for zc2 in [False,True]:
for d in [0.1,1,3.7]:
for itp in [np.float64,np.complex128,np.float32,np.complex64]:
a = RGSpace(dim1, zerocenter=zc1, distances=d)
b = RGRGTransformation.get_codomain(a, zerocenter=zc2)
fft = FFTOperator(domain=a, target=b, domain_dtype=itp, target_dtype=_harmonic_type(itp))
inp = Field.from_random(domain=a,random_type='normal',std=7,mean=3,dtype=itp)
out = fft.inverse_times(fft.times(inp))
assert_allclose(inp.val, out.val)
def test_fft2D(self):
for dim1 in [10,11]:
for dim2 in [9,12]:
for zc1 in [False,True]:
for zc2 in [False,True]:
for zc3 in [False,True]:
for zc4 in [False,True]:
for d in [0.1,1,3.7]:
for itp in [np.float64,np.complex128,np.float32,np.complex64]:
a = RGSpace([dim1,dim2], zerocenter=[zc1,zc2], distances=d)
b = RGRGTransformation.get_codomain(a, zerocenter=[zc3,zc4])
fft = FFTOperator(domain=a, target=b, domain_dtype=itp, target_dtype=_harmonic_type(itp))
inp = Field.from_random(domain=a,random_type='normal',std=7,mean=3,dtype=itp)
out = fft.inverse_times(fft.times(inp))
assert_allclose(inp.val, out.val)
def test_sht(self):
if 'pyHealpix' not in di:
raise SkipTest
for lm in [0,3,6,11,30]:
for tp in [np.float64,np.complex128,np.float32,np.complex64]:
a = LMSpace(lmax=lm)
b = LMGLTransformation.get_codomain(a)
fft = FFTOperator(domain=a, target=b, domain_dtype=tp, target_dtype=tp)
inp = Field.from_random(domain=a,random_type='normal',std=7,mean=3,dtype=tp)
out = fft.inverse_times(fft.times(inp))
assert_allclose(inp.val, out.val)
def test_sht2(self):
if 'pyHealpix' not in di:
raise SkipTest
for lm in [128,256]:
for tp in [np.float64,np.complex128,np.float32,np.complex64]:
a = LMSpace(lmax=lm)
b = LMHPTransformation.get_codomain(a)
fft = FFTOperator(domain=a, target=b, domain_dtype=tp, target_dtype=tp)
inp = Field.from_random(domain=a,random_type='normal',std=1,mean=0,dtype=tp)
out = fft.inverse_times(fft.times(inp))
assert_allclose(inp.val, out.val,rtol=1e-3,atol=1e-3)
......@@ -21,6 +21,7 @@ import unittest
from numpy.testing import assert_equal
from keepers import Repository
from test.common import expand, generate_spaces
import os
try:
import h5py
......@@ -44,3 +45,7 @@ if h5py_available:
self._repo.commit()
assert_equal(space, self._repo.get('space'))
@classmethod
def tearDownClass(cls):
os.remove('test.h5')
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