Commit b3a06798 authored by Theo Steininger's avatar Theo Steininger

Fixed hermitianization. Moved preserve_gaussian_covariance completely to the Field.

Made tests more complete.
parent 4a36f41c
Pipeline #14517 passed with stage
in 6 minutes and 14 seconds
......@@ -599,58 +599,55 @@ class Field(Loggable, Versionable, object):
# hermitianize for the first space
(h, a) = domain[spaces[0]].hermitian_decomposition(
val,
domain_axes[spaces[0]],
preserve_gaussian_variance=preserve_gaussian_variance)
domain_axes[spaces[0]])
# 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(
h,
domain_axes[space],
preserve_gaussian_variance=False)
domain_axes[space])
(ah, aa) = domain[space].hermitian_decomposition(
a,
domain_axes[space],
preserve_gaussian_variance=False)
domain_axes[space])
c = (hh - ha - ah + aa).conjugate()
full = (hh + ha + ah + aa)
h = (full + c)/2.
a = (full - c)/2.
# correct variance
# 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]
# for product_point in itertools.product(*fixed_points):
# slice_object = np.array((slice(None), )*len(val.shape),
# dtype=np.object)
# for i, sp in enumerate(spaces):
# point_component = product_point[i]
# if point_component is None:
# point_component = slice(None)
# slice_object[list(domain_axes[sp])] = point_component
#
# slice_object = tuple(slice_object)
# h[slice_object] /= np.sqrt(2)
# a[slice_object] /= np.sqrt(2)
if preserve_gaussian_variance:
h *= np.sqrt(2)
a *= np.sqrt(2)
if not issubclass(val.dtype.type, np.complexfloating):
# 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]
fixed_points = [[fp] if fp is None else fp
for fp in fixed_points]
for product_point in itertools.product(*fixed_points):
slice_object = np.array((slice(None), )*len(val.shape),
dtype=np.object)
for i, sp in enumerate(spaces):
point_component = product_point[i]
if point_component is None:
point_component = slice(None)
slice_object[list(domain_axes[sp])] = point_component
slice_object = tuple(slice_object)
h[slice_object] /= np.sqrt(2)
a[slice_object] /= np.sqrt(2)
return (h, a)
def _spec_to_rescaler(self, spec, result_list, power_space_index):
......
......@@ -89,25 +89,21 @@ class LMSpace(Space):
super(LMSpace, self).__init__()
self._lmax = self._parse_lmax(lmax)
def hermitian_decomposition(self, x, axes=None,
preserve_gaussian_variance=False):
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
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
anti_hermitian_part[:] = 0
return (hermitian_part, anti_hermitian_part)
# def hermitian_fixed_points(self):
# return None
def hermitian_fixed_points(self):
return None
# ---Mandatory properties and methods---
......
......@@ -118,43 +118,25 @@ class RGSpace(Space):
# use subtraction since it is faster than flipping another time
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)
def _hermitianize_correct_variance(self, hermitian_part,
anti_hermitian_part, axes):
# Correct the variance by multiplying sqrt(2)
hermitian_part = hermitian_part * np.sqrt(2)
anti_hermitian_part = anti_hermitian_part * np.sqrt(2)
# If the dtype of the input is complex, the fixed points lose the power
# of their imaginary-part (or real-part, respectively). Therefore
# the factor of sqrt(2) also applies there
if not issubclass(hermitian_part.dtype.type, np.complexfloating):
# 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
dimensions = mid_index.size
# Use ndindex to iterate over all combinations of zeros and the
# mid_index in order to correct all fixed points.
ndlist=[1]*dimensions
for k in range(len(axes)):
i = axes[k]
if self.shape[k]%2 == 0:
ndlist[i] = 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
def hermitian_fixed_points(self):
dimensions = len(self.shape)
mid_index = np.array(self.shape)//2
ndlist = [1]*dimensions
for k in range(dimensions):
if self.shape[k] % 2 == 0:
ndlist[k] = 2
ndlist = tuple(ndlist)
fixed_points = []
for index in np.ndindex(ndlist):
for k in range(dimensions):
if self.shape[k] % 2 != 0 and self.zerocenter[k]:
index = list(index)
index[k] = 1
index = tuple(index)
fixed_points += [tuple(index * mid_index)]
return fixed_points
def _hermitianize_inverter(self, x, axes):
# calculate the number of dimensions the input array has
......
......@@ -20,21 +20,17 @@ import unittest
import numpy as np
from numpy.testing import assert_,\
assert_equal,\
assert_almost_equal
from itertools import product
from nifty import Field,\
RGSpace,\
FieldArray
RGSpace
from d2o import distributed_data_object,\
STRATEGIES
from d2o import distributed_data_object
from test.common import expand
np.random.seed(123)
SPACES = [RGSpace((4,)), RGSpace((5))]
SPACE_COMBINATIONS = [(), SPACES[0], SPACES[1], SPACES]
......@@ -56,46 +52,31 @@ class Test_Interface(unittest.TestCase):
f = Field(domain=domain)
assert_(isinstance(getattr(f, attribute), desired_type))
def test_hermitian_decomposition0(self):
s1=(25,)
s2=(16,)
r1 = RGSpace(s1, harmonic=True)
r2 = RGSpace(s2, harmonic=True)
ra = RGSpace(s1+s2, harmonic=True)
v = np.random.random(s1+s2) + 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,),),False)
h2,a2 = Field._hermitian_decomposition((r1,r2),f2.val,(0,1),((0,),(1,)),False)
h3,a3 = Field._hermitian_decomposition((r1,r2),f2.val,(1,0),((0,),(1,)),False)
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())
@expand(product([False,True],[False,True]))
def test_hermitian_decomposition1(self, complexdata, preserve):
s0=(1,)
s1=(56,25)
r0 = RGSpace(s0, harmonic=True)
r1 = RGSpace(s1, harmonic=True)
ra = RGSpace(s0+s1, harmonic=True)
v = np.random.random(s1)
if (complexdata):
v = v + 1j*np.random.random(s1)
f1=Field(r1,val=v,copy=True)
f2=Field((r0,r1),val=v,copy=True)
h1,a1 = Field._hermitian_decomposition((r1,),f1.val,(0,),((0,1,),),preserve)
h2,a2 = Field._hermitian_decomposition((r0,r1),f2.val,(0,1),((0,),(1,2)),preserve)
h2=h2[0,:,:]
a2=a2[0,:,:]
assert_almost_equal(h1.get_full_data(),h2.get_full_data())
assert_almost_equal(a1.get_full_data(),a2.get_full_data())
class Test_Functionality(unittest.TestCase):
@expand(product([True, False], [True, False],
[True, False], [True, False],
[(1,), (4,), (5,)], [(1,), (6,), (7,)]))
def test_hermitian_decomposition(self, z1, z2, preserve, complexdata,
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))
#class Test_Initialization(unittest.TestCase):
#
# @parameterized.expand(
# itertools.product(SPACE_COMBINATIONS,
# []
# )
# def test_
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):
def test_distance_array(self, lmax, expected):
l = LMSpace(lmax)
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)
......@@ -23,7 +23,6 @@ import numpy as np
from numpy.testing import assert_, assert_equal, assert_almost_equal
from nifty import RGSpace
from nifty import Field
from test.common import expand
from itertools import product
......@@ -153,47 +152,59 @@ class RGSpaceFunctionalityTests(unittest.TestCase):
for key, value in expected.iteritems():
assert_equal(getattr(x, key), value)
@expand(product([(10,),(11,),(1,1),(4,4),(5,7),(8,12),(7,16),(4,6,8),
(17,5,3)],))
def test_hermitian_decomposition(self, shape):
r = RGSpace(shape, harmonic=True)
v = np.empty(shape,dtype=np.complex128)
@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):
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)
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)
# 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'])
it = np.nditer(h, flags=['multi_index'])
while not it.finished:
i1 = it.multi_index
i2 = []
for i in range(len(i1)):
i2.append(h.shape[i]-i1[i] if i1[i]>0 else 0)
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]))
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)],))
def test_hermitian_decomposition2(self, shape):
r = RGSpace(shape, harmonic=True)
@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)
h, a = r.hermitian_decomposition(v)
# make sure that data == h + a
assert_almost_equal(v,h+a)
assert_almost_equal(v, h+a)
# test hermitianity of h
it = np.nditer (h, flags=['multi_index'])
it = np.nditer(h, flags=['multi_index'])
while not it.finished:
i1 = it.multi_index
i2 = []
for i in range(len(i1)):
i2.append(h.shape[i]-i1[i] if i1[i]>0 else 0)
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]))
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())
......@@ -209,3 +220,8 @@ class RGSpaceFunctionalityTests(unittest.TestCase):
assert_almost_equal(res, expected)
if inplace:
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)])
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