Commit 935bc62b authored by Theo Steininger's avatar Theo Steininger

Simplified _hermitian_decomposition according to Reimars suggestion.

parent aab0727b
...@@ -596,22 +596,15 @@ class Field(Loggable, Versionable, object): ...@@ -596,22 +596,15 @@ class Field(Loggable, Versionable, object):
@staticmethod @staticmethod
def _hermitian_decomposition(domain, val, spaces, domain_axes, def _hermitian_decomposition(domain, val, spaces, domain_axes,
preserve_gaussian_variance=False): preserve_gaussian_variance=False):
# hermitianize for the first space
(h, a) = domain[spaces[0]].hermitian_decomposition( flipped_val = val
val, for space in spaces:
domain_axes[spaces[0]]) flipped_val = domain[space].hermitianize_inverter(
# hermitianize all remaining spaces using the iterative formula x=flipped_val,
for space in spaces[1:]: axes=domain_axes[space])
(hh, ha) = domain[space].hermitian_decomposition( flipped_val = flipped_val.conjugate()
h, h = (val + flipped_val)/2.
domain_axes[space]) a = val - h
(ah, aa) = domain[space].hermitian_decomposition(
a,
domain_axes[space])
c = (hh - ha - ah + aa).conjugate()
full = (hh + ha + ah + aa)
h = (full + c)/2.
a = (full - c)/2.
# correct variance # correct variance
if preserve_gaussian_variance: if preserve_gaussian_variance:
......
...@@ -117,7 +117,6 @@ class FFTOperator(LinearOperator): ...@@ -117,7 +117,6 @@ class FFTOperator(LinearOperator):
super(FFTOperator, self).__init__(default_spaces) super(FFTOperator, self).__init__(default_spaces)
# Initialize domain and target # Initialize domain and target
self._domain = self._parse_domain(domain) self._domain = self._parse_domain(domain)
if len(self.domain) != 1: if len(self.domain) != 1:
raise ValueError("TransformationOperator accepts only exactly one " raise ValueError("TransformationOperator accepts only exactly one "
......
...@@ -89,22 +89,6 @@ class LMSpace(Space): ...@@ -89,22 +89,6 @@ class LMSpace(Space):
super(LMSpace, self).__init__() super(LMSpace, self).__init__()
self._lmax = self._parse_lmax(lmax) self._lmax = self._parse_lmax(lmax)
def hermitian_decomposition(self, x, axes=None):
if issubclass(x.dtype.type, np.complexfloating):
hermitian_part = x.copy_empty()
anti_hermitian_part = x.copy_empty()
hermitian_part[:] = x.real
anti_hermitian_part[:] = x.imag * 1j
else:
hermitian_part = x.copy()
anti_hermitian_part = x.copy_empty()
anti_hermitian_part[:] = 0
return (hermitian_part, anti_hermitian_part)
def hermitian_fixed_points(self):
return None
# ---Mandatory properties and methods--- # ---Mandatory properties and methods---
def __repr__(self): def __repr__(self):
......
...@@ -100,26 +100,6 @@ class RGSpace(Space): ...@@ -100,26 +100,6 @@ class RGSpace(Space):
self._distances = self._parse_distances(distances) self._distances = self._parse_distances(distances)
self._zerocenter = self._parse_zerocenter(zerocenter) self._zerocenter = self._parse_zerocenter(zerocenter)
def hermitian_decomposition(self, x, axes=None,
preserve_gaussian_variance=False):
# check axes
if axes is None:
axes = range(len(self.shape))
assert len(x.shape) >= len(self.shape), "shapes mismatch"
assert len(axes) == len(self.shape), "axes mismatch"
# compute the hermitian part
flipped_x = self._hermitianize_inverter(x, axes=axes)
flipped_x = flipped_x.conjugate()
# average x and flipped_x.
hermitian_part = x + flipped_x
hermitian_part /= 2.
# use subtraction since it is faster than flipping another time
anti_hermitian_part = (x-hermitian_part)
return (hermitian_part, anti_hermitian_part)
def hermitian_fixed_points(self): def hermitian_fixed_points(self):
dimensions = len(self.shape) dimensions = len(self.shape)
mid_index = np.array(self.shape)//2 mid_index = np.array(self.shape)//2
...@@ -138,7 +118,7 @@ class RGSpace(Space): ...@@ -138,7 +118,7 @@ class RGSpace(Space):
fixed_points += [tuple(index * mid_index)] fixed_points += [tuple(index * mid_index)]
return fixed_points return fixed_points
def _hermitianize_inverter(self, x, axes): def hermitianize_inverter(self, x, axes):
# calculate the number of dimensions the input array has # calculate the number of dimensions the input array has
dimensions = len(x.shape) dimensions = len(x.shape)
# prepare the slicing object which will be used for mirroring # prepare the slicing object which will be used for mirroring
......
...@@ -147,35 +147,34 @@ class Space(DomainObject): ...@@ -147,35 +147,34 @@ class Space(DomainObject):
raise NotImplementedError( raise NotImplementedError(
"There is no generic co-smoothing kernel for Space base class.") "There is no generic co-smoothing kernel for Space base class.")
def hermitian_decomposition(self, x, axes, def hermitian_fixed_points(self):
preserve_gaussian_variance=False): """ Returns the array points which remain invariant under the action
""" Decomposes x into its hermitian and anti-hermitian constituents. of `hermitianize_inverter`
This method decomposes a field's array x into its hermitian and
antihermitian part, which corresponds to real and imaginary part Returns
in a corresponding harmonic partner space. This is an internal function -------
that is mainly used for power-synthesizing and -analyzing Fields. list of index-tuples
The list contains the index-coordinates of the invariant points.
"""
return None
def hermitianize_inverter(self, x, axes):
""" Inverts/flips x in the context of Hermitian decomposition.
This method is mainly used for power-synthesizing and -analyzing
Fields.
Parameters Parameters
---------- ----------
x : distributed_data_object
The field's val object.
axes : tuple of ints axes : tuple of ints
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*
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
they lose half of their variance since the real(/imaginary) part
gets lost.
Returns Returns
------- -------
(distributed_data_object, distributed_data_object) distributed_data_object
A tuple of two distributed_data_objects, the first being the The Hermitian-flipped of x.
hermitian and the second the anti-hermitian part of x.
""" """
raise NotImplementedError return x
...@@ -155,56 +155,25 @@ class RGSpaceFunctionalityTests(unittest.TestCase): ...@@ -155,56 +155,25 @@ class RGSpaceFunctionalityTests(unittest.TestCase):
@expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16), @expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16),
(4, 6, 8), (17, 5, 3)], (4, 6, 8), (17, 5, 3)],
[True, False])) [True, False]))
def test_hermitian_decomposition(self, shape, zerocenter): def test_hermitianize_inverter(self, shape, zerocenter):
r = RGSpace(shape, harmonic=True, zerocenter=zerocenter) r = RGSpace(shape, harmonic=True, zerocenter=zerocenter)
v = np.empty(shape, dtype=np.complex128) v = np.empty(shape, dtype=np.complex128)
v.real = np.random.random(shape) v.real = np.random.random(shape)
v.imag = np.random.random(shape) v.imag = np.random.random(shape)
h, a = r.hermitian_decomposition(v) inverted = r.hermitianize_inverter(v, axes=range(len(shape)))
# make sure that data == h + a
# NOTE: this is only correct for preserve_gaussian_variance==False,
# but I consider this an intrinsic property of a hermitian
# decomposition.
assert_almost_equal(v, h+a)
print (h, a)
# test hermitianity of h
it = np.nditer(h, flags=['multi_index'])
while not it.finished:
i1 = it.multi_index
i2 = []
for i in range(len(i1)):
if r.zerocenter[i] and r.shape[i] % 2 != 0:
i2.append(h.shape[i]-i1[i]-1)
else:
i2.append(h.shape[i]-i1[i] if i1[i] > 0 else 0)
i2 = tuple(i2)
assert_almost_equal(h[i1], np.conj(h[i2]))
assert_almost_equal(a[i1], -np.conj(a[i2]))
it.iternext()
@expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16), # test hermitian flipping of `inverted`
(4, 6, 8), (17, 5, 3)], it = np.nditer(v, flags=['multi_index'])
[True, False]))
def test_hermitian_decomposition2(self, shape, zerocenter):
r = RGSpace(shape, harmonic=True, zerocenter=zerocenter)
v = np.random.random(shape)
h, a = r.hermitian_decomposition(v)
# make sure that data == h + a
assert_almost_equal(v, h+a)
# test hermitianity of h
it = np.nditer(h, flags=['multi_index'])
while not it.finished: while not it.finished:
i1 = it.multi_index i1 = it.multi_index
i2 = [] i2 = []
for i in range(len(i1)): for i in range(len(i1)):
if r.zerocenter[i] and r.shape[i] % 2 != 0: if r.zerocenter[i] and r.shape[i] % 2 != 0:
i2.append(h.shape[i]-i1[i]-1) i2.append(v.shape[i]-i1[i]-1)
else: else:
i2.append(h.shape[i]-i1[i] if i1[i] > 0 else 0) i2.append(v.shape[i]-i1[i] if i1[i] > 0 else 0)
i2 = tuple(i2) i2 = tuple(i2)
assert_almost_equal(h[i1], np.conj(h[i2])) assert_almost_equal(inverted[i1], v[i2])
assert_almost_equal(a[i1], -np.conj(a[i2]))
it.iternext() it.iternext()
@expand(get_distance_array_configs()) @expand(get_distance_array_configs())
......
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