Commit 41a5faf9 by Theo Steininger

### Merge branch 'fix_hermitianizer' into 'master'

```Fix hermitianizer

Closes #157

See merge request !157```
parents 50d6a969 b3a06798
Pipeline #14518 passed with stages
in 12 minutes and 10 seconds
 ... @@ -599,58 +599,55 @@ class Field(Loggable, Versionable, object): ... @@ -599,58 +599,55 @@ class Field(Loggable, Versionable, object): # 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=preserve_gaussian_variance) # 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 spaces[1:]: (hh, ha) = domain[space].hermitian_decomposition( (hh, ha) = domain[space].hermitian_decomposition( h, h, domain_axes[space], domain_axes[space]) preserve_gaussian_variance=False) (ah, aa) = domain[space].hermitian_decomposition( (ah, aa) = domain[space].hermitian_decomposition( a, a, domain_axes[space], domain_axes[space]) preserve_gaussian_variance=False) c = (hh - ha - ah + aa).conjugate() c = (hh - ha - ah + aa).conjugate() full = (hh + ha + ah + aa) full = (hh + ha + ah + aa) h = (full + c)/2. h = (full + c)/2. a = (full - c)/2. a = (full - c)/2. # correct variance # correct variance if preserve_gaussian_variance: # in principle one must not correct the variance for the fixed h *= np.sqrt(2) # points of the hermitianization. However, for a complex field a *= np.sqrt(2) # the input field loses half of its power at its fixed points # in the `hermitian` part. Hence, here a factor of sqrt(2) is if not issubclass(val.dtype.type, np.complexfloating): # also necessary! # in principle one must not correct the variance for the fixed # => The hermitianization can be done on a space level since either # points of the hermitianization. However, for a complex field # nothing must be done (LMSpace) or ALL points need a factor of sqrt(2) # the input field loses half of its power at its fixed points # => use the preserve_gaussian_variance flag in the # in the `hermitian` part. Hence, here a factor of sqrt(2) is # hermitian_decomposition method above. # also necessary! # => The hermitianization can be done on a space level since # This code is for educational purposes: # either nothing must be done (LMSpace) or ALL points need a # fixed_points = [domain[i].hermitian_fixed_points() for i in spaces] # factor of sqrt(2) # # check if there was at least one flipping during hermitianization # => use the preserve_gaussian_variance flag in the # flipped_Q = np.any([fp is not None for fp in fixed_points]) # hermitian_decomposition method above. # # if the array got flipped, correct the variance # if flipped_Q: # This code is for educational purposes: # h *= np.sqrt(2) fixed_points = [domain[i].hermitian_fixed_points() # a *= np.sqrt(2) for i in spaces] # fixed_points = [[fp] if fp is None else fp # fixed_points = [[fp] if fp is None else fp for fp in fixed_points] for fp in fixed_points] # for product_point in itertools.product(*fixed_points): # slice_object = np.array((slice(None), )*len(val.shape), for product_point in itertools.product(*fixed_points): # dtype=np.object) slice_object = np.array((slice(None), )*len(val.shape), # for i, sp in enumerate(spaces): dtype=np.object) # point_component = product_point[i] for i, sp in enumerate(spaces): # if point_component is None: point_component = product_point[i] # point_component = slice(None) if point_component is None: # slice_object[list(domain_axes[sp])] = point_component point_component = slice(None) # slice_object[list(domain_axes[sp])] = point_component # slice_object = tuple(slice_object) # h[slice_object] /= np.sqrt(2) slice_object = tuple(slice_object) # a[slice_object] /= np.sqrt(2) h[slice_object] /= np.sqrt(2) a[slice_object] /= np.sqrt(2) return (h, a) return (h, a) def _spec_to_rescaler(self, spec, result_list, power_space_index): def _spec_to_rescaler(self, spec, result_list, power_space_index): ... @@ -667,7 +664,7 @@ class Field(Loggable, Versionable, object): ... @@ -667,7 +664,7 @@ class Field(Loggable, Versionable, object): if pindex.distribution_strategy is not local_distribution_strategy: if pindex.distribution_strategy is not local_distribution_strategy: self.logger.warn( self.logger.warn( "The distribution_stragey of pindex does not fit the " "The distribution_strategy of pindex does not fit the " "slice_local distribution strategy of the synthesized field.") "slice_local distribution strategy of the synthesized field.") # Now use numpy advanced indexing in order to put the entries of the # Now use numpy advanced indexing in order to put the entries of the ... ...
 ... @@ -89,25 +89,21 @@ class LMSpace(Space): ... @@ -89,25 +89,21 @@ 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, def hermitian_decomposition(self, x, axes=None): preserve_gaussian_variance=False): if issubclass(x.dtype.type, np.complexfloating): if issubclass(x.dtype.type, np.complexfloating): hermitian_part = x.copy_empty() hermitian_part = x.copy_empty() anti_hermitian_part = x.copy_empty() anti_hermitian_part = x.copy_empty() hermitian_part[:] = x.real hermitian_part[:] = x.real anti_hermitian_part[:] = x.imag * 1j anti_hermitian_part[:] = x.imag * 1j if preserve_gaussian_variance: hermitian_part *= np.sqrt(2) anti_hermitian_part *= np.sqrt(2) else: else: hermitian_part = x.copy() hermitian_part = x.copy() anti_hermitian_part = x.copy_empty() anti_hermitian_part = x.copy_empty() anti_hermitian_part.val[:] = 0 anti_hermitian_part[:] = 0 return (hermitian_part, anti_hermitian_part) return (hermitian_part, anti_hermitian_part) # def hermitian_fixed_points(self): def hermitian_fixed_points(self): # return None return None # ---Mandatory properties and methods--- # ---Mandatory properties and methods--- ... ...
 ... @@ -102,6 +102,12 @@ class RGSpace(Space): ... @@ -102,6 +102,12 @@ class RGSpace(Space): def hermitian_decomposition(self, x, axes=None, def hermitian_decomposition(self, x, axes=None, preserve_gaussian_variance=False): 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 # compute the hermitian part flipped_x = self._hermitianize_inverter(x, axes=axes) flipped_x = self._hermitianize_inverter(x, axes=axes) flipped_x = flipped_x.conjugate() flipped_x = flipped_x.conjugate() ... @@ -112,68 +118,46 @@ class RGSpace(Space): ... @@ -112,68 +118,46 @@ class RGSpace(Space): # use subtraction since it is faster than flipping another time # use subtraction since it is faster than flipping another time anti_hermitian_part = (x-hermitian_part) anti_hermitian_part = (x-hermitian_part) if preserve_gaussian_variance: hermitian_part, anti_hermitian_part = \ self._hermitianize_correct_variance(hermitian_part, anti_hermitian_part, axes=axes) return (hermitian_part, anti_hermitian_part) return (hermitian_part, anti_hermitian_part) def _hermitianize_correct_variance(self, hermitian_part, def hermitian_fixed_points(self): anti_hermitian_part, axes): dimensions = len(self.shape) # Correct the variance by multiplying sqrt(2) mid_index = np.array(self.shape)//2 hermitian_part = hermitian_part * np.sqrt(2) ndlist = [1]*dimensions anti_hermitian_part = anti_hermitian_part * np.sqrt(2) for k in range(dimensions): if self.shape[k] % 2 == 0: # If the dtype of the input is complex, the fixed points lose the power ndlist[k] = 2 # of their imaginary-part (or real-part, respectively). Therefore ndlist = tuple(ndlist) # the factor of sqrt(2) also applies there fixed_points = [] if not issubclass(hermitian_part.dtype.type, np.complexfloating): for index in np.ndindex(ndlist): # The fixed points of the point inversion must not be averaged. for k in range(dimensions): # Hence one must divide out the sqrt(2) again if self.shape[k] % 2 != 0 and self.zerocenter[k]: # -> Get the middle index of the array index = list(index) mid_index = np.array(hermitian_part.shape, dtype=np.int) // 2 index[k] = 1 dimensions = mid_index.size index = tuple(index) # Use ndindex to iterate over all combinations of zeros and the fixed_points += [tuple(index * mid_index)] # mid_index in order to correct all fixed points. return 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) hermitian_part[temp_index] /= np.sqrt(2) anti_hermitian_part[temp_index] /= np.sqrt(2) return hermitian_part, anti_hermitian_part def _hermitianize_inverter(self, x, axes): def _hermitianize_inverter(self, x, axes): shape = x.shape # calculate the number of dimensions the input array has # calculate the number of dimensions the input array has dimensions = len(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 slice_primitive = [slice(None), ] * dimensions slice_primitive = [slice(None), ] * dimensions # copy the input data # copy the input data y = x.copy() y = x.copy() if axes is None: axes = xrange(dimensions) # flip in the desired directions # flip in the desired directions for i in axes: for k in range(len(axes)): i = axes[k] slice_picker = slice_primitive[:] slice_picker = slice_primitive[:] if shape[i] % 2 == 0: slice_picker[i] = slice(1, None, None) else: slice_picker[i] = slice(None) slice_picker = tuple(slice_picker) slice_inverter = slice_primitive[:] slice_inverter = slice_primitive[:] if shape[i] % 2 == 0: if (not self.zerocenter[k]) or self.shape[k] % 2 == 0: slice_picker[i] = slice(1, None, None) slice_inverter[i] = slice(None, 0, -1) slice_inverter[i] = slice(None, 0, -1) else: else: slice_picker[i] = slice(None) slice_inverter[i] = slice(None, None, -1) slice_inverter[i] = slice(None, None, -1) slice_picker = tuple(slice_picker) slice_inverter = tuple(slice_inverter) slice_inverter = tuple(slice_inverter) try: try: ... ...
 ... @@ -167,7 +167,7 @@ class Space(DomainObject): ... @@ -167,7 +167,7 @@ class Space(DomainObject): If the hermitian decomposition is done via computing the half If the hermitian decomposition is done via computing the half sums and differences of `x` and mirrored `x`, all points except the sums and differences of `x` and mirrored `x`, all points except the fixed points lose half of their variance. If `x` is complex also fixed points lose half of their variance. If `x` is complex also the lose half of their variance since the real(/imaginary) part they lose half of their variance since the real(/imaginary) part gets lost. gets lost. Returns Returns ... ...
 ... @@ -20,20 +20,17 @@ import unittest ... @@ -20,20 +20,17 @@ import unittest import numpy as np import numpy as np from numpy.testing import assert_,\ from numpy.testing import assert_,\ assert_equal assert_almost_equal from itertools import product from itertools import product from nifty import Field,\ from nifty import Field,\ RGSpace,\ RGSpace FieldArray from d2o import distributed_data_object,\ from d2o import distributed_data_object STRATEGIES from test.common import expand from test.common import expand np.random.seed(123) SPACES = [RGSpace((4,)), RGSpace((5))] SPACES = [RGSpace((4,)), RGSpace((5))] SPACE_COMBINATIONS = [(), SPACES[0], SPACES[1], SPACES] SPACE_COMBINATIONS = [(), SPACES[0], SPACES[1], SPACES] ... @@ -55,10 +52,31 @@ class Test_Interface(unittest.TestCase): ... @@ -55,10 +52,31 @@ class Test_Interface(unittest.TestCase): f = Field(domain=domain) f = Field(domain=domain) assert_(isinstance(getattr(f, attribute), desired_type)) assert_(isinstance(getattr(f, attribute), desired_type)) #class Test_Initialization(unittest.TestCase): # class Test_Functionality(unittest.TestCase): # @parameterized.expand( @expand(product([True, False], [True, False], # itertools.product(SPACE_COMBINATIONS, [True, False], [True, False], # [] [(1,), (4,), (5,)], [(1,), (6,), (7,)])) # ) def test_hermitian_decomposition(self, z1, z2, preserve, complexdata, # def test_ s1, s2): np.random.seed(123) r1 = RGSpace(s1, harmonic=True, zerocenter=(z1,)) r2 = RGSpace(s2, harmonic=True, zerocenter=(z2,)) ra = RGSpace(s1+s2, harmonic=True, zerocenter=(z1, z2)) v = np.random.random(s1+s2) if complexdata: v = v + 1j*np.random.random(s1+s2) f1 = Field(ra, val=v, copy=True) f2 = Field((r1, r2), val=v, copy=True) h1, a1 = Field._hermitian_decomposition((ra,), f1.val, (0,), ((0, 1,),), preserve) h2, a2 = Field._hermitian_decomposition((r1, r2), f2.val, (0, 1), ((0,), (1,)), preserve) h3, a3 = Field._hermitian_decomposition((r1, r2), f2.val, (1, 0), ((0,), (1,)), preserve) assert_almost_equal(h1.get_full_data(), h2.get_full_data()) assert_almost_equal(a1.get_full_data(), a2.get_full_data()) assert_almost_equal(h1.get_full_data(), h3.get_full_data()) assert_almost_equal(a1.get_full_data(), a3.get_full_data())
 ... @@ -130,3 +130,7 @@ class LMSpaceFunctionalityTests(unittest.TestCase): ... @@ -130,3 +130,7 @@ class LMSpaceFunctionalityTests(unittest.TestCase): def test_distance_array(self, lmax, expected): def test_distance_array(self, lmax, expected): l = LMSpace(lmax) l = LMSpace(lmax) assert_almost_equal(l.get_distance_array('not').data, expected) assert_almost_equal(l.get_distance_array('not').data, expected) def test_hermitian_fixed_points(self): x = LMSpace(5) assert_equal(x.hermitian_fixed_points(), None)
 ... @@ -24,6 +24,7 @@ import numpy as np ... @@ -24,6 +24,7 @@ import numpy as np from numpy.testing import assert_, assert_equal, assert_almost_equal from numpy.testing import assert_, assert_equal, assert_almost_equal from nifty import RGSpace from nifty import RGSpace from test.common import expand from test.common import expand from itertools import product # [shape, zerocenter, distances, harmonic, expected] # [shape, zerocenter, distances, harmonic, expected] CONSTRUCTOR_CONFIGS = [ CONSTRUCTOR_CONFIGS = [ ... @@ -135,83 +136,6 @@ def get_weight_configs(): ... @@ -135,83 +136,6 @@ def get_weight_configs(): ] ] def get_hermitian_configs(): h_0_x = np.array([ [0.88250339+0.12102381j, 0.54293435+0.7345584j, 0.87057998+0.20515315j, 0.16602950+0.09396132j], [0.83853902+0.17974696j, 0.79735933+0.37104425j, 0.22057732+0.9498977j, 0.14329183+0.47899678j], [0.96934284+0.3792878j, 0.13118669+0.45643055j, 0.16372149+0.48235714j, 0.66141537+0.20383357j], [0.49168197+0.77572178j, 0.09570420+0.14219071j, 0.69735595+0.33017333j, 0.83692452+0.18544449j]]) h_0_res_real = np.array([ [0.88250339+0.j, 0.35448193+0.32029854j, 0.87057998+0.j, 0.35448193-0.32029854j], [0.66511049-0.29798741j, 0.81714193+0.09279988j, 0.45896664+0.30986218j, 0.11949801+0.16840303j], [0.96934284+0.j, 0.39630103+0.12629849j, 0.16372149+0.j, 0.39630103-0.12629849j], [0.66511049+0.29798741j, 0.11949801-0.16840303j, 0.45896664-0.30986218j, 0.81714193-0.09279988j]]) h_0_res_imag = np.array([ [0.12102381+0.j, 0.41425986-0.18845242j, 0.20515315+0.j, 0.41425986+0.18845242j], [0.47773437-0.17342852j, 0.27824437+0.0197826j, 0.64003551+0.23838932j, 0.31059374-0.02379381j], [0.37928780+0.j, 0.33013206+0.26511434j, 0.48235714+0.j, 0.33013206-0.26511434j], [0.47773437+0.17342852j, 0.31059374+0.02379381j, 0.64003551-0.23838932j, 0.27824437-0.0197826j]])*1j h_1_x = np.array([ [[0.23987021+0.41617749j, 0.34605012+0.55462234j, 0.07947035+0.73360723j, 0.22853748+0.39275304j], [0.90254910+0.02107809j, 0.28195470+0.56031588j, 0.23004043+0.33873536j, 0.56398377+0.68913034j], [0.81897406+0.2050369j, 0.88724852+0.8137488j, 0.84645004+0.0059284j, 0.14950377+0.50013099j]], [[0.93491597+0.73251066j, 0.74764790+0.11539037j, 0.48090736+0.04352568j, 0.49363732+0.97233093j], [0.72761881+0.74636216j, 0.46390134+0.4343401j, 0.88436859+0.79415269j, 0.67027606+0.85498234j], [0.86318727+0.19076379j, 0.36859448+0.89842333j, 0.73407193+0.85091112j, 0.44187657+0.08936409j]] ]) h_1_res_real = np.array([ [[0.23987021+0.j, 0.28729380+0.08093465j, 0.07947035+0.j, 0.28729380-0.08093465j], [0.90254910+0.j, 0.42296924-0.06440723j, 0.23004043+0.j, 0.42296924+0.06440723j], [0.81897406+0.j, 0.51837614+0.1568089j, 0.84645004+0.j, 0.51837614-0.1568089j]], [[0.93491597+0.j, 0.62064261-0.42847028j, 0.48090736+0.j, 0.62064261+0.42847028j], [0.72761881+0.j, 0.56708870-0.21032112j, 0.88436859+0.j, 0.56708870+0.21032112j], [0.86318727+0.j, 0.40523552+0.40452962j, 0.73407193+0.j, 0.40523552-0.40452962j]] ]) h_1_res_imag = np.array([ [[0.41617749+0.j, 0.47368769-0.05875632j, 0.73360723+0.j, 0.47368769+0.05875632j], [0.02107809+0.j, 0.62472311+0.14101454j, 0.33873536+0.j, 0.62472311-0.14101454j], [0.20503690+0.j, 0.65693990-0.36887238j, 0.00592840+0.j, 0.65693990+0.36887238j]], [[0.73251066+0.j, 0.54386065-0.12700529j, 0.04352568+0.j, 0.54386065+0.12700529j], [0.74636216+0.j, 0.64466122+0.10318736j, 0.79415269+0.j, 0.64466122-0.10318736j], [0.19076379+0.j, 0.49389371+0.03664104j, 0.85091112+0.j, 0.49389371-0.03664104j]] ])*1j return [ [h_0_x, None, h_0_res_real, h_0_res_imag], [h_1_x, (2,), h_1_res_real, h_1_res_imag] ] class RGSpaceInterfaceTests(unittest.TestCase): class RGSpaceInterfaceTests(unittest.TestCase): @expand([['distances', tuple], @expand([['distances', tuple], ['zerocenter', tuple]]) ['zerocenter', tuple]]) ... @@ -228,11 +152,60 @@ class RGSpaceFunctionalityTests(unittest.TestCase): ... @@ -228,11 +152,60 @@ class RGSpaceFunctionalityTests(unittest.TestCase): for key, value in expected.iteritems(): for key, value in expected.iteritems(): assert_equal(getattr(x, key), value) assert_equal(getattr(x, key), value) @expand(get_hermitian_configs()) @expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16), def test_hermitian_decomposition(self, x, axes, real, imag): (4, 6, 8), (17, 5, 3)], r = RGSpace(5) [True, False])) assert_almost_equal(r.hermitian_decomposition(x, axes=axes)[0], real) def test_hermitian_decomposition(self, shape, zerocenter): assert_almost_equal(r.hermitian_decomposition(x, axes=axes)[1], imag) 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() @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']) 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(get_distance_array_configs()) @expand(get_distance_array_configs()) def test_distance_array(self, shape, distances, zerocenter, expected): def test_distance_array(self, shape, distances, zerocenter, expected): ... @@ -247,3 +220,8 @@ class RGSpaceFunctionalityTests(unittest.TestCase): ... @@ -247,3 +220,8 @@ class RGSpaceFunctionalityTests(unittest.TestCase): assert_almost_equal(res, expected) assert_almost_equal(res, expected) if inplace: if inplace: assert_(x is res) assert_(x is res) def test_hermitian_fixed_points(self): x = RGSpace((5, 6, 5, 6), zerocenter=[False, False, True, True]) assert_equal(x.hermitian_fixed_points(), [(0, 0, 2, 0), (0, 0, 2, 3), (0, 3, 2, 0), (0, 3, 2, 3)])
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!