Commit d985e8ee authored by Theo Steininger's avatar Theo Steininger

Fixed power for power_synthesize on LMSpace.

Added first methods for __repr__
parent dd80d59d
Pipeline #13043 failed with stages
in 5 minutes and 16 seconds
...@@ -45,6 +45,10 @@ class DomainObject(Versionable, Loggable, object): ...@@ -45,6 +45,10 @@ class DomainObject(Versionable, Loggable, object):
# _global_id is used in the Versioning module from keepers # _global_id is used in the Versioning module from keepers
self._ignore_for_hash = ['_global_id'] self._ignore_for_hash = ['_global_id']
@abc.abstractmethod
def __repr__(self):
raise NotImplementedError
def __hash__(self): def __hash__(self):
# Extract the identifying parts from the vars(self) dict. # Extract the identifying parts from the vars(self) dict.
result_hash = 0 result_hash = 0
......
...@@ -386,6 +386,7 @@ class Field(Loggable, Versionable, object): ...@@ -386,6 +386,7 @@ class Field(Loggable, Versionable, object):
for part in [hermitian_part, anti_hermitian_part]] for part in [hermitian_part, anti_hermitian_part]]
power_spectrum = hermitian_power + 1j * anti_hermitian_power power_spectrum = hermitian_power + 1j * anti_hermitian_power
else: else:
power_spectrum = self._calculate_power_spectrum( power_spectrum = self._calculate_power_spectrum(
x=self.val, x=self.val,
...@@ -396,11 +397,7 @@ class Field(Loggable, Versionable, object): ...@@ -396,11 +397,7 @@ class Field(Loggable, Versionable, object):
# create the result field and put power_spectrum into it # create the result field and put power_spectrum into it
result_domain = list(self.domain) result_domain = list(self.domain)
result_domain[space_index] = power_domain result_domain[space_index] = power_domain
result_dtype = power_spectrum.dtype
if decompose_power:
result_dtype = np.complex
else:
result_dtype = np.float
result_field = self.copy_empty( result_field = self.copy_empty(
domain=result_domain, domain=result_domain,
...@@ -585,33 +582,44 @@ class Field(Loggable, Versionable, object): ...@@ -585,33 +582,44 @@ class Field(Loggable, Versionable, object):
def _hermitian_decomposition(domain, val, spaces, domain_axes): def _hermitian_decomposition(domain, val, spaces, domain_axes):
# hermitianize for the first space # hermitianize for the first space
(h, a) = domain[spaces[0]].hermitian_decomposition( (h, a) = domain[spaces[0]].hermitian_decomposition(
val, val,
domain_axes[spaces[0]]) domain_axes[spaces[0]],
preserve_gaussian_variance=True)
# hermitianize all remaining spaces using the iterative formula # hermitianize all remaining spaces using the iterative formula
for space in xrange(1, len(spaces)): for space in xrange(1, len(spaces)):
(hh, ha) = \ (hh, ha) = domain[space].hermitian_decomposition(
domain[space].hermitian_decomposition(h, domain_axes[space]) h,
(ah, aa) = \ domain_axes[space],
domain[space].hermitian_decomposition(a, domain_axes[space]) preserve_gaussian_variance=True)
(ah, aa) = domain[space].hermitian_decomposition(
a,
domain_axes[space],
preserve_gaussian_variance=True)
c = (hh - ha - ah + aa).conjugate() c = (hh - ha - ah + aa).conjugate()
h = (val + c)/2. h = (val + c)/2.
a = (val - c)/2. a = (val - c)/2.
# correct variance # correct variance
fixed_points = [domain[i].hermitian_fixed_points() for i in spaces]
# check if there was at least one flipping during hermitianization
flipped_Q = np.any([fp is not None for fp in fixed_points])
# if the array got flipped, correct the variance
if flipped_Q:
h *= np.sqrt(2)
a *= np.sqrt(2)
# in principle one must not correct the variance for the fixed
# points of the hermitianization. However, for a complex field
# the input field looses half of its power at its fixed points
# in the `hermitian` part. Hence, here a factor of sqrt(2) is
# also necessary!
# in principle one must not correct the variance for the fixed
# points of the hermitianization. However, for a complex field
# the input field loses half of its power at its fixed points
# in the `hermitian` part. Hence, here a factor of sqrt(2) is
# also necessary!
# => The hermitianization can be done on a space level since either
# nothing must be done (LMSpace) or ALL points need a factor of sqrt(2)
# => use the preserve_gaussian_variance flag in the
# hermitian_decomposition method above.
# This code is for educational purposes:
# fixed_points = [domain[i].hermitian_fixed_points() for i in spaces]
# # check if there was at least one flipping during hermitianization
# flipped_Q = np.any([fp is not None for fp in fixed_points])
# # if the array got flipped, correct the variance
# if flipped_Q:
# h *= np.sqrt(2)
# a *= np.sqrt(2)
#
# fixed_points = [[fp] if fp is None else fp for fp in fixed_points] # fixed_points = [[fp] if fp is None else fp for fp in fixed_points]
# for product_point in itertools.product(*fixed_points): # for product_point in itertools.product(*fixed_points):
# slice_object = np.array((slice(None), )*len(val.shape), # slice_object = np.array((slice(None), )*len(val.shape),
......
...@@ -100,6 +100,9 @@ class GLSpace(Space): ...@@ -100,6 +100,9 @@ class GLSpace(Space):
# ---Mandatory properties and methods--- # ---Mandatory properties and methods---
def __repr__(self):
return ("GLSpace(nlat=%r, nlon=%r)" % (self.nlat, self.nlon))
@property @property
def harmonic(self): def harmonic(self):
return False return False
......
...@@ -84,6 +84,9 @@ class HPSpace(Space): ...@@ -84,6 +84,9 @@ class HPSpace(Space):
# ---Mandatory properties and methods--- # ---Mandatory properties and methods---
def __repr__(self):
return ("HPSpace(nside=%r)" % self.nside)
@property @property
def harmonic(self): def harmonic(self):
return False return False
......
...@@ -91,10 +91,19 @@ class LMSpace(Space): ...@@ -91,10 +91,19 @@ class LMSpace(Space):
def hermitian_decomposition(self, x, axes=None, def hermitian_decomposition(self, x, axes=None,
preserve_gaussian_variance=False): preserve_gaussian_variance=False):
hermitian_part = x.copy_empty() if issubclass(x.dtype.type, np.complexfloating):
anti_hermitian_part = x.copy_empty() hermitian_part = x.copy_empty()
hermitian_part[:] = x.real anti_hermitian_part = x.copy_empty()
anti_hermitian_part[:] = x.imag * 1j hermitian_part[:] = x.real
anti_hermitian_part[:] = x.imag * 1j
if preserve_gaussian_variance:
hermitian_part *= np.sqrt(2)
anti_hermitian_part *= np.sqrt(2)
else:
hermitian_part = x.copy()
anti_hermitian_part = x.copy_empty()
anti_hermitian_part.val[:] = 0
return (hermitian_part, anti_hermitian_part) return (hermitian_part, anti_hermitian_part)
# def hermitian_fixed_points(self): # def hermitian_fixed_points(self):
...@@ -102,6 +111,9 @@ class LMSpace(Space): ...@@ -102,6 +111,9 @@ class LMSpace(Space):
# ---Mandatory properties and methods--- # ---Mandatory properties and methods---
def __repr__(self):
return ("LMSpace(lmax=%r)" % self.lmax)
@property @property
def harmonic(self): def harmonic(self):
return True return True
......
...@@ -149,6 +149,13 @@ class PowerSpace(Space): ...@@ -149,6 +149,13 @@ class PowerSpace(Space):
# ---Mandatory properties and methods--- # ---Mandatory properties and methods---
def __repr__(self):
return ("PowerSpace(harmonic_partner=%r, distribution_strategy=%r, "
"logarithmic=%r, nbin=%r, binbounds=%r)"
% (self.harmonic_partner, self.pindex.distribution_strategy,
self.config['logarithmic'], self.config['nbin'],
self.config['binbounds']))
@property @property
def harmonic(self): def harmonic(self):
return True return True
......
...@@ -52,10 +52,10 @@ class RGSpace(Space): ...@@ -52,10 +52,10 @@ class RGSpace(Space):
---------- ----------
shape : {int, numpy.ndarray} shape : {int, numpy.ndarray}
Number of grid points or numbers of gridpoints along each axis. Number of grid points or numbers of gridpoints along each axis.
zerocenter : {bool, numpy.ndarray}, *optional* zerocenter : {bool, numpy.ndarray} *optional*
Whether x==0 (or k==0, respectively) is located in the center of Whether x==0 (or k==0, respectively) is located in the center of
the grid (or the center of each axis speparately) or not. the grid (or the center of each axis speparately) or not.
(default: False). (default: False).
distances : {float, numpy.ndarray}, *optional* distances : {float, numpy.ndarray}, *optional*
Distance between two grid points along each axis Distance between two grid points along each axis
(default: None). (default: None).
...@@ -120,40 +120,32 @@ class RGSpace(Space): ...@@ -120,40 +120,32 @@ class RGSpace(Space):
return (hermitian_part, anti_hermitian_part) return (hermitian_part, anti_hermitian_part)
# def hermitian_fixed_points(self):
# shape = self.shape
# mid_index = np.array(shape)//2
# ndlist = [2 if (shape[i] % 2 == 0) else 1 for i in xrange(len(shape))]
# ndlist = tuple(ndlist)
# odd_axes_list = np.array([1 if (shape[i] % 2 == 1) else 0
# for i in xrange(len(shape))])
# fixed_points = []
# for i in np.ndindex(ndlist):
# fixed_points += [tuple((i+odd_axes_list) * mid_index)]
# return fixed_points
def _hermitianize_correct_variance(self, hermitian_part, def _hermitianize_correct_variance(self, hermitian_part,
anti_hermitian_part, axes): anti_hermitian_part, axes):
# Correct the variance by multiplying sqrt(2) # Correct the variance by multiplying sqrt(2)
hermitian_part = hermitian_part * np.sqrt(2) hermitian_part = hermitian_part * np.sqrt(2)
anti_hermitian_part = anti_hermitian_part * np.sqrt(2) anti_hermitian_part = anti_hermitian_part * np.sqrt(2)
# The fixed points of the point inversion must not be averaged. # If the dtype of the input is complex, the fixed points lose the power
# Hence one must divide out the sqrt(2) again # of their imaginary-part (or real-part, respectively). Therefore
# -> Get the middle index of the array # the factor of sqrt(2) also applies there
mid_index = np.array(hermitian_part.shape, dtype=np.int) // 2 if not issubclass(hermitian_part.dtype.type, np.complexfloating):
dimensions = mid_index.size # The fixed points of the point inversion must not be averaged.
# Use ndindex to iterate over all combinations of zeros and the # Hence one must divide out the sqrt(2) again
# mid_index in order to correct all fixed points. # -> Get the middle index of the array
if axes is None: mid_index = np.array(hermitian_part.shape, dtype=np.int) // 2
axes = xrange(dimensions) dimensions = mid_index.size
# Use ndindex to iterate over all combinations of zeros and the
ndlist = [2 if i in axes else 1 for i in xrange(dimensions)] # mid_index in order to correct all fixed points.
ndlist = tuple(ndlist) if axes is None:
for i in np.ndindex(ndlist): axes = xrange(dimensions)
temp_index = tuple(i * mid_index)
hermitian_part[temp_index] /= np.sqrt(2) ndlist = [2 if i in axes else 1 for i in xrange(dimensions)]
anti_hermitian_part[temp_index] /= np.sqrt(2) ndlist = tuple(ndlist)
for i in np.ndindex(ndlist):
temp_index = tuple(i * mid_index)
hermitian_part[temp_index] /= np.sqrt(2)
anti_hermitian_part[temp_index] /= np.sqrt(2)
return hermitian_part, anti_hermitian_part return hermitian_part, anti_hermitian_part
def _hermitianize_inverter(self, x, axes): def _hermitianize_inverter(self, x, axes):
...@@ -193,6 +185,10 @@ class RGSpace(Space): ...@@ -193,6 +185,10 @@ class RGSpace(Space):
# ---Mandatory properties and methods--- # ---Mandatory properties and methods---
def __repr__(self):
return ("RGSpace(shape=%r, zerocenter=%r, distances=%r, harmonic=%r)"
% (self.shape, self.zerocenter, self.distances, self.harmonic))
@property @property
def harmonic(self): def harmonic(self):
return self._harmonic return self._harmonic
......
...@@ -64,12 +64,6 @@ class Space(DomainObject): ...@@ -64,12 +64,6 @@ class Space(DomainObject):
@abc.abstractproperty @abc.abstractproperty
def harmonic(self): def harmonic(self):
""" Returns True if this space is a harmonic space. """ Returns True if this space is a harmonic space.
Raises
------
NotImplementedError
If called for this abstract class.
""" """
raise NotImplementedError raise NotImplementedError
...@@ -83,12 +77,8 @@ class Space(DomainObject): ...@@ -83,12 +77,8 @@ class Space(DomainObject):
float float
A real number representing the sum of all pixel volumes. A real number representing the sum of all pixel volumes.
Raises
------
NotImplementedError
If called for this abstract class.
""" """
raise NotImplementedError( raise NotImplementedError(
"There is no generic volume for the Space base class.") "There is no generic volume for the Space base class.")
...@@ -122,11 +112,6 @@ class Space(DomainObject): ...@@ -122,11 +112,6 @@ class Space(DomainObject):
distributed_data_object distributed_data_object
A d2o containing the distances A d2o containing the distances
Raises
------
NotImplementedError
If called for this abstract class.
""" """
raise NotImplementedError( raise NotImplementedError(
...@@ -157,11 +142,6 @@ class Space(DomainObject): ...@@ -157,11 +142,6 @@ class Space(DomainObject):
A smoothing operation that multiplies values with a Gaussian A smoothing operation that multiplies values with a Gaussian
kernel. kernel.
Raises
------
NotImplementedError :
If called for this abstract class.
""" """
raise NotImplementedError( raise NotImplementedError(
...@@ -184,7 +164,11 @@ class Space(DomainObject): ...@@ -184,7 +164,11 @@ class Space(DomainObject):
Specifies the axes of x which correspond to this space. Specifies the axes of x which correspond to this space.
preserve_gaussian_variance : bool *optional* preserve_gaussian_variance : bool *optional*
FIXME: figure out what this does If the hermitian decomposition is done via computing the half
sums and differences of `x` and mirrored `x`, all points except the
fixed points lose half of their variance. If `x` is complex also
the lose half of their variance since the real(/imaginary) part
gets lost.
Returns Returns
------- -------
...@@ -192,14 +176,6 @@ class Space(DomainObject): ...@@ -192,14 +176,6 @@ class Space(DomainObject):
A tuple of two distributed_data_objects, the first being the A tuple of two distributed_data_objects, the first being the
hermitian and the second the anti-hermitian part of x. hermitian and the second the anti-hermitian part of x.
Raises
------
NotImplementedError
If called for this abstract class.
""" """
raise NotImplementedError raise NotImplementedError
def __repr__(self):
return str(type(self)) + "\n"
...@@ -18,14 +18,15 @@ ...@@ -18,14 +18,15 @@
import unittest import unittest
import numpy as np import numpy as np
from numpy.testing import assert_, assert_equal
from itertools import product from itertools import product
from types import LambdaType from types import LambdaType
from numpy.testing import assert_, assert_raises, assert_equal
from nifty import LMSpace, GLSpace, HPSpace
from nifty.config import dependency_injector as di
from test.common import expand, generate_spaces, generate_harmonic_spaces from test.common import expand, generate_spaces, generate_harmonic_spaces
from d2o import distributed_data_object
from nifty.spaces import *
class SpaceInterfaceTests(unittest.TestCase): class SpaceInterfaceTests(unittest.TestCase):
@expand(product(generate_spaces(), [ @expand(product(generate_spaces(), [
...@@ -34,28 +35,17 @@ class SpaceInterfaceTests(unittest.TestCase): ...@@ -34,28 +35,17 @@ class SpaceInterfaceTests(unittest.TestCase):
['dim', int], ['dim', int],
['total_volume', np.float]])) ['total_volume', np.float]]))
def test_property_ret_type(self, space, attr_expected_type): def test_property_ret_type(self, space, attr_expected_type):
assert_( assert_(isinstance(getattr(space, attr_expected_type[0]),
isinstance(getattr( attr_expected_type[1]))
space,
attr_expected_type[0]
), attr_expected_type[1])
)
@expand(product(generate_harmonic_spaces(), [ @expand(product(generate_harmonic_spaces(), [
['get_fft_smoothing_kernel_function', None, LambdaType], ['get_distance_array', 'not', distributed_data_object],
['get_fft_smoothing_kernel_function', 2.0, LambdaType], ['get_fft_smoothing_kernel_function', 2.0, LambdaType],
])) ]))
def test_method_ret_type(self, space, method_expected_type): def test_method_ret_type(self, space, method_expected_type):
getattr( assert_(type(getattr(space, method_expected_type[0])(
space, method_expected_type[0])(*method_expected_type[1:-1]) *method_expected_type[1:-1])) is
method_expected_type[-1])
assert_equal(
type(getattr(
space,
method_expected_type[0])(*method_expected_type[1:-1])
),
method_expected_type[-1]
)
@expand([[space] for space in generate_spaces()]) @expand([[space] for space in generate_spaces()])
def test_copy(self, space): def test_copy(self, space):
...@@ -63,3 +53,7 @@ class SpaceInterfaceTests(unittest.TestCase): ...@@ -63,3 +53,7 @@ class SpaceInterfaceTests(unittest.TestCase):
assert_(space is not space.copy()) assert_(space is not space.copy())
# make sure contents are the same # make sure contents are the same
assert_equal(space, space.copy()) assert_equal(space, space.copy())
@expand([[space] for space in generate_spaces()])
def test_repr(self, space):
assert_(space == eval(space.__repr__()))
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