diff --git a/ci/requirements.txt b/ci/requirements.txt
index 638b653efb5f31e6cf852aac378c95015804ca5e..424f78f5a6757773b236896f08586eead7ed6c8a 100644
--- a/ci/requirements.txt
+++ b/ci/requirements.txt
@@ -1,4 +1,5 @@
numpy
+cython
mpi4py
matplotlib
plotly
diff --git a/ci/requirements_base.txt b/ci/requirements_base.txt
index 597d5963d0196f59628b63e25a4cc8681a490f59..5bf7be86d73f06787e4c438cd590bd480e8805e2 100644
--- a/ci/requirements_base.txt
+++ b/ci/requirements_base.txt
@@ -1,4 +1,5 @@
numpy
+cython
nose
parameterized
coverage
diff --git a/nifty/field.py b/nifty/field.py
index ab647ab8fdee61da005869afe749e8bf5783c6f9..fe01bca857a21b12d0dd0a7ac389aa36f4e88e46 100644
--- a/nifty/field.py
+++ b/nifty/field.py
@@ -17,6 +17,8 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
from __future__ import division
+
+import itertools
import numpy as np
from keepers import Versionable,\
@@ -104,6 +106,7 @@ class Field(Loggable, Versionable, object):
distributed_data_object
"""
+
# ---Initialization methods---
def __init__(self, domain=None, val=None, dtype=None,
@@ -512,7 +515,7 @@ class Field(Loggable, Versionable, object):
result_domain = list(self.domain)
for power_space_index in spaces:
power_space = self.domain[power_space_index]
- harmonic_domain = power_space.harmonic_domain
+ harmonic_domain = power_space.harmonic_partner
result_domain[power_space_index] = harmonic_domain
# create random samples: one or two, depending on whether the
@@ -554,14 +557,12 @@ class Field(Loggable, Versionable, object):
inplace=True)
if real_signal:
- for power_space_index in spaces:
- harmonic_domain = result_domain[power_space_index]
- result_val_list = [harmonic_domain.hermitian_decomposition(
- result_val,
- axes=result.domain_axes[power_space_index],
- preserve_gaussian_variance=True)[0]
- for (result, result_val)
- in zip(result_list, result_val_list)]
+ result_val_list = [self._hermitian_decomposition(
+ result_domain,
+ result_val,
+ spaces,
+ result_list[0].domain_axes)[0]
+ for result_val in result_val_list]
# store the result into the fields
[x.set_val(new_val=y, copy=False) for x, y in
@@ -574,6 +575,46 @@ class Field(Loggable, Versionable, object):
return result
+ @staticmethod
+ def _hermitian_decomposition(domain, val, spaces, domain_axes):
+ # 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 xrange(1, len(spaces)):
+ (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()
+ h = (val + c)/2.
+ a = (val - c)/2.
+
+ # correct variance
+ 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)
+
+ return (h, a)
+
def _spec_to_rescaler(self, spec, result_list, power_space_index):
power_space = self.domain[power_space_index]
diff --git a/nifty/spaces/lm_space/lm_space.py b/nifty/spaces/lm_space/lm_space.py
index 9ba38e7902198c6b0a530d263cf9e46ff5eb2212..5cadf785bdb0bd89947d704242fb0fea21db56ce 100644
--- a/nifty/spaces/lm_space/lm_space.py
+++ b/nifty/spaces/lm_space/lm_space.py
@@ -94,9 +94,12 @@ class LMSpace(Space):
hermitian_part = x.copy_empty()
anti_hermitian_part = x.copy_empty()
hermitian_part[:] = x.real
- anti_hermitian_part[:] = x.imag
+ anti_hermitian_part[:] = x.imag * 1j
return (hermitian_part, anti_hermitian_part)
+ def hermitian_fixed_points(self):
+ return None
+
# ---Mandatory properties and methods---
@property
diff --git a/nifty/spaces/rg_space/rg_space.py b/nifty/spaces/rg_space/rg_space.py
index cb70c66452be8bc762b53ea91801ac077ee91304..887e52de8431296ed9d8de35e8f7dd5259f5afde 100644
--- a/nifty/spaces/rg_space/rg_space.py
+++ b/nifty/spaces/rg_space/rg_space.py
@@ -110,7 +110,7 @@ class RGSpace(Space):
hermitian_part /= 2.
# use subtraction since it is faster than flipping another time
- anti_hermitian_part = (x-hermitian_part)/1j
+ anti_hermitian_part = (x-hermitian_part)
if preserve_gaussian_variance:
hermitian_part, anti_hermitian_part = \
@@ -120,6 +120,18 @@ class RGSpace(Space):
return (hermitian_part, anti_hermitian_part)
+ def hermitian_fixed_points(self):
+ shape = self.shape
+ mid_index = np.array(shape)//2
+ ndlist = [2 if (shape[i] % 2 == 0) else 1 for i in xrange(len(shape))]
+ ndlist = tuple(ndlist)
+ odd_axes_list = np.array([1 if (shape[i] % 2 == 1) else 0
+ for i in xrange(len(shape))])
+ fixed_points = []
+ for i in np.ndindex(ndlist):
+ fixed_points += [tuple((i+odd_axes_list) * mid_index)]
+ return fixed_points
+
def _hermitianize_correct_variance(self, hermitian_part,
anti_hermitian_part, axes):
# Correct the variance by multiplying sqrt(2)
@@ -145,8 +157,9 @@ class RGSpace(Space):
return hermitian_part, anti_hermitian_part
def _hermitianize_inverter(self, x, axes):
+ shape = x.shape
# calculate the number of dimensions the input array has
- dimensions = len(x.shape)
+ dimensions = len(shape)
# prepare the slicing object which will be used for mirroring
slice_primitive = [slice(None), ] * dimensions
# copy the input data
@@ -158,11 +171,17 @@ class RGSpace(Space):
# flip in the desired directions
for i in axes:
slice_picker = slice_primitive[:]
- slice_picker[i] = slice(1, None, None)
+ 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[i] = slice(None, 0, -1)
+ if shape[i] % 2 == 0:
+ slice_inverter[i] = slice(None, 0, -1)
+ else:
+ slice_inverter[i] = slice(None, None, -1)
slice_inverter = tuple(slice_inverter)
try: