From 935bc62baa8031d5750757c5bd2ca378b3018893 Mon Sep 17 00:00:00 2001 From: Theo Steininger Date: Thu, 13 Jul 2017 02:03:14 +0200 Subject: [PATCH] Simplified _hermitian_decomposition according to Reimars suggestion. --- nifty/field.py | 25 ++++------- nifty/operators/fft_operator/fft_operator.py | 1 - nifty/spaces/lm_space/lm_space.py | 16 ------- nifty/spaces/rg_space/rg_space.py | 22 +--------- nifty/spaces/space/space.py | 41 +++++++++--------- test/test_spaces/test_rg_space.py | 45 +++----------------- 6 files changed, 37 insertions(+), 113 deletions(-) diff --git a/nifty/field.py b/nifty/field.py index 67e092d0..9dce821f 100644 --- a/nifty/field.py +++ b/nifty/field.py @@ -596,22 +596,15 @@ class Field(Loggable, Versionable, object): @staticmethod def _hermitian_decomposition(domain, val, spaces, domain_axes, preserve_gaussian_variance=False): - # hermitianize for the first space - (h, a) = domain[spaces[0]].hermitian_decomposition( - val, - domain_axes[spaces[0]]) - # hermitianize all remaining spaces using the iterative formula - for space in spaces[1:]: - (hh, ha) = domain[space].hermitian_decomposition( - h, - domain_axes[space]) - (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. + + flipped_val = val + for space in spaces: + flipped_val = domain[space].hermitianize_inverter( + x=flipped_val, + axes=domain_axes[space]) + flipped_val = flipped_val.conjugate() + h = (val + flipped_val)/2. + a = val - h # correct variance if preserve_gaussian_variance: diff --git a/nifty/operators/fft_operator/fft_operator.py b/nifty/operators/fft_operator/fft_operator.py index c04bfbbb..c5dd2c3d 100644 --- a/nifty/operators/fft_operator/fft_operator.py +++ b/nifty/operators/fft_operator/fft_operator.py @@ -117,7 +117,6 @@ class FFTOperator(LinearOperator): super(FFTOperator, self).__init__(default_spaces) # Initialize domain and target - self._domain = self._parse_domain(domain) if len(self.domain) != 1: raise ValueError("TransformationOperator accepts only exactly one " diff --git a/nifty/spaces/lm_space/lm_space.py b/nifty/spaces/lm_space/lm_space.py index 32bf0690..42022f05 100644 --- a/nifty/spaces/lm_space/lm_space.py +++ b/nifty/spaces/lm_space/lm_space.py @@ -89,22 +89,6 @@ class LMSpace(Space): super(LMSpace, self).__init__() 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--- def __repr__(self): diff --git a/nifty/spaces/rg_space/rg_space.py b/nifty/spaces/rg_space/rg_space.py index fdb17107..a78136db 100644 --- a/nifty/spaces/rg_space/rg_space.py +++ b/nifty/spaces/rg_space/rg_space.py @@ -100,26 +100,6 @@ class RGSpace(Space): self._distances = self._parse_distances(distances) 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): dimensions = len(self.shape) mid_index = np.array(self.shape)//2 @@ -138,7 +118,7 @@ class RGSpace(Space): fixed_points += [tuple(index * mid_index)] return fixed_points - def _hermitianize_inverter(self, x, axes): + def hermitianize_inverter(self, 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 diff --git a/nifty/spaces/space/space.py b/nifty/spaces/space/space.py index 4b239141..d9b37ddc 100644 --- a/nifty/spaces/space/space.py +++ b/nifty/spaces/space/space.py @@ -147,35 +147,34 @@ class Space(DomainObject): raise NotImplementedError( "There is no generic co-smoothing kernel for Space base class.") - def hermitian_decomposition(self, x, axes, - preserve_gaussian_variance=False): - """ Decomposes x into its hermitian and anti-hermitian constituents. + def hermitian_fixed_points(self): + """ Returns the array points which remain invariant under the action + of `hermitianize_inverter` - This method decomposes a field's array x into its hermitian and - antihermitian part, which corresponds to real and imaginary part - in a corresponding harmonic partner space. This is an internal function - that is mainly used for power-synthesizing and -analyzing Fields. + + Returns + ------- + 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 ---------- - x : distributed_data_object - The field's val object. axes : tuple of ints 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 ------- - (distributed_data_object, distributed_data_object) - A tuple of two distributed_data_objects, the first being the - hermitian and the second the anti-hermitian part of x. - + distributed_data_object + The Hermitian-flipped of x. """ - raise NotImplementedError + return x diff --git a/test/test_spaces/test_rg_space.py b/test/test_spaces/test_rg_space.py index 8e434d06..2b3e348e 100644 --- a/test/test_spaces/test_rg_space.py +++ b/test/test_spaces/test_rg_space.py @@ -155,56 +155,25 @@ class RGSpaceFunctionalityTests(unittest.TestCase): @expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16), (4, 6, 8), (17, 5, 3)], [True, False])) - def test_hermitian_decomposition(self, shape, zerocenter): + def test_hermitianize_inverter(self, shape, zerocenter): r = RGSpace(shape, harmonic=True, zerocenter=zerocenter) v = np.empty(shape, dtype=np.complex128) v.real = np.random.random(shape) v.imag = np.random.random(shape) - h, a = r.hermitian_decomposition(v) - # 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() + inverted = r.hermitianize_inverter(v, axes=range(len(shape))) - @expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16), - (4, 6, 8), (17, 5, 3)], - [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']) + # test hermitian flipping of `inverted` + it = np.nditer(v, 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) + i2.append(v.shape[i]-i1[i]-1) 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) - assert_almost_equal(h[i1], np.conj(h[i2])) - assert_almost_equal(a[i1], -np.conj(a[i2])) + assert_almost_equal(inverted[i1], v[i2]) it.iternext() @expand(get_distance_array_configs()) -- GitLab