Skip to content

GitLab

N
NIFTy
Commit 470d9aab authored by Martin Reinecke's avatar Martin Reinecke
Browse files
Options

Merge branch 'master' into vl_bfgs

parents 5d972816 41a5faf9
Pipeline #14521 passed with stage
    in 7 minutes and 3 seconds
    Hide whitespace changes
    Inline Side-by-side
    Showing with 169 additions and 192 deletions
    nifty/field.py
    View file @ 470d9aab
    ......@@ -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):
    ......@@ -667,7 +664,7 @@ class Field(Loggable, Versionable, object):
    if pindex.distribution_strategy is not local_distribution_strategy:
    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.")
    # Now use numpy advanced indexing in order to put the entries of the
    ......
    nifty/spaces/lm_space/lm_space.py
    View file @ 470d9aab
    ......@@ -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---
    ......
    nifty/spaces/rg_space/rg_space.py
    View file @ 470d9aab
    ......@@ -102,6 +102,12 @@ class RGSpace(Space):
    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()
    ......@@ -112,68 +118,46 @@ 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.
    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 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):
    shape = x.shape
    # 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
    slice_primitive = [slice(None), ] * dimensions
    # copy the input data
    y = x.copy()
    if axes is None:
    axes = xrange(dimensions)
    # flip in the desired directions
    for i in axes:
    for k in range(len(axes)):
    i = axes[k]
    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[:]
    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)
    else:
    slice_picker[i] = slice(None)
    slice_inverter[i] = slice(None, None, -1)
    slice_picker = tuple(slice_picker)
    slice_inverter = tuple(slice_inverter)
    try:
    ......
    nifty/spaces/space/space.py
    View file @ 470d9aab
    ......@@ -167,7 +167,7 @@ class Space(DomainObject):
    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
    they lose half of their variance since the real(/imaginary) part
    gets lost.
    Returns
    ......
    test/test_field.py
    View file @ 470d9aab
    ......@@ -20,20 +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]
    ......@@ -55,10 +52,31 @@ class Test_Interface(unittest.TestCase):
    f = Field(domain=domain)
    assert_(isinstance(getattr(f, attribute), desired_type))
    #class Test_Initialization(unittest.TestCase):
    #
    # @parameterized.expand(
    # itertools.product(SPACE_COMBINATIONS,
    # []
    # )
    # def test_
    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))
    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())
    test/test_spaces/test_lm_space.py
    View file @ 470d9aab
    ......@@ -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)
    test/test_spaces/test_rg_space.py
    View file @ 470d9aab
    ......@@ -24,6 +24,7 @@ import numpy as np
    from numpy.testing import assert_, assert_equal, assert_almost_equal
    from nifty import RGSpace
    from test.common import expand
    from itertools import product
    # [shape, zerocenter, distances, harmonic, expected]
    CONSTRUCTOR_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):
    @expand([['distances', tuple],
    ['zerocenter', tuple]])
    ......@@ -228,11 +152,60 @@ class RGSpaceFunctionalityTests(unittest.TestCase):
    for key, value in expected.iteritems():
    assert_equal(getattr(x, key), value)
    @expand(get_hermitian_configs())
    def test_hermitian_decomposition(self, x, axes, real, imag):
    r = RGSpace(5)
    assert_almost_equal(r.hermitian_decomposition(x, axes=axes)[0], real)
    assert_almost_equal(r.hermitian_decomposition(x, axes=axes)[1], imag)
    @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)
    # 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())
    def test_distance_array(self, shape, distances, zerocenter, expected):
    ......@@ -247,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