Merge branch 'Real-Field_Synthesizer' into 'master'

Real field synthesizer

See merge request !125

ci/requirements.txt

 numpy
+cython
 mpi4py
 matplotlib
 plotly

ci/requirements_base.txt

 numpy
+cython
 nose
 parameterized
 coverage

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,
@@ -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]
 

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

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:

test/test_spaces/test_lm_space.py

@@ -116,7 +116,7 @@ class LMSpaceFunctionalityTests(unittest.TestCase):
             real)
         assert_almost_equal(
             l.hermitian_decomposition(distributed_data_object(x))[1],
-            imag)
+            imag*1j)
 
     @expand(get_weight_configs())
     def test_weight(self, x, power, axes, inplace, expected):

test/test_spaces/test_rg_space.py

@@ -162,7 +162,7 @@ def get_hermitian_configs():
         [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]])
+            0.27824437-0.0197826j]])*1j
 
     h_1_x = np.array([
         [[0.23987021+0.41617749j, 0.34605012+0.55462234j, 0.07947035+0.73360723j,
@@ -205,7 +205,7 @@ def get_hermitian_configs():
              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]