merge from nifty2go

nifty/operators/fft_operator_support.py

@@ -74,52 +74,15 @@ class RGRGTransformation(Transformation):
         return Tval
 
 
-class SlicingTransformation(Transformation):
-    def transform(self, x):
-        axes = x.domain.axes[self.space]
-        if dobj.dist_axis(x.val) in axes:
-            raise NotImplementedError
-        p2h = x.domain == self.pdom
-        idat = dobj.local_data(x.val)
-        if p2h:
-            res = Field(self.hdom, dtype=x.dtype)
-            odat = dobj.local_data(res.val)
-
-            for slice in utilities.get_slice_list(idat.shape, axes):
-                odat[slice] = self._slice_p2h(idat[slice])
-        else:
-            res = Field(self.pdom, dtype=x.dtype)
-            odat = dobj.local_data(res.val)
-
-            for slice in utilities.get_slice_list(idat.shape, axes):
-                odat[slice] = self._slice_h2p(idat[slice])
-        return res
-
-
-def buildLm(nr, lmax):
-    res = np.empty((len(nr)+lmax+1)//2, dtype=(nr[0]*1j).dtype)
-    res[0:lmax+1] = nr[0:lmax+1]
-    res[lmax+1:] = np.sqrt(0.5)*(nr[lmax+1::2] + 1j*nr[lmax+2::2])
-    return res
-
-
-def buildIdx(nr, lmax):
-    res = np.empty((lmax+1)*(lmax+1), dtype=nr[0].real.dtype)
-    res[0:lmax+1] = nr[0:lmax+1].real
-    res[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
-    res[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
-    return res
-
-
-class SphericalTransformation(SlicingTransformation):
+class SphericalTransformation(Transformation):
     def __init__(self, pdom, hdom, space):
         super(SphericalTransformation, self).__init__(pdom, hdom, space)
         from pyHealpix import sharpjob_d
 
         self.lmax = self.hdom[self.space].lmax
+        self.mmax = self.hdom[self.space].mmax
         self.sjob = sharpjob_d()
-        self.sjob.set_triangular_alm_info(self.hdom[self.space].lmax,
-                                          self.hdom[self.space].mmax)
+        self.sjob.set_triangular_alm_info(self.lmax, self.mmax)
         if isinstance(self.pdom[self.space], GLSpace):
             self.sjob.set_Gauss_geometry(self.pdom[self.space].nlat,
                                          self.pdom[self.space].nlon)
@@ -132,8 +95,39 @@ class SphericalTransformation(SlicingTransformation):
 
     def _slice_p2h(self, inp):
         rr = self.sjob.map2alm(inp)
-        return buildIdx(rr, lmax=self.lmax)
+        assert len(rr) == ((self.mmax+1)*(self.mmax+2))//2 + \
+                          (self.mmax+1)*(self.lmax-self.mmax)
+        res = np.empty(2*len(rr)-self.lmax-1, dtype=rr[0].real.dtype)
+        res[0:self.lmax+1] = rr[0:self.lmax+1].real
+        res[self.lmax+1::2] = np.sqrt(2)*rr[self.lmax+1:].real
+        res[self.lmax+2::2] = np.sqrt(2)*rr[self.lmax+1:].imag
+        return res
 
     def _slice_h2p(self, inp):
-        result = buildLm(inp, lmax=self.lmax)
-        return self.sjob.alm2map(result)
+        res = np.empty((len(inp)+self.lmax+1)//2, dtype=(inp[0]*1j).dtype)
+        assert len(res) == ((self.mmax+1)*(self.mmax+2))//2 + \
+                           (self.mmax+1)*(self.lmax-self.mmax)
+        res[0:self.lmax+1] = inp[0:self.lmax+1]
+        res[self.lmax+1:] = np.sqrt(0.5)*(inp[self.lmax+1::2] +
+                                          1j*inp[self.lmax+2::2])
+        return self.sjob.alm2map(res)
+
+    def transform(self, x):
+        axes = x.domain.axes[self.space]
+        if dobj.dist_axis(x.val) in axes:
+            raise NotImplementedError
+        p2h = x.domain == self.pdom
+        idat = dobj.local_data(x.val)
+        if p2h:
+            res = Field(self.hdom, dtype=x.dtype)
+            odat = dobj.local_data(res.val)
+
+            for slice in utilities.get_slice_list(idat.shape, axes):
+                odat[slice] = self._slice_p2h(idat[slice])
+        else:
+            res = Field(self.pdom, dtype=x.dtype)
+            odat = dobj.local_data(res.val)
+
+            for slice in utilities.get_slice_list(idat.shape, axes):
+                odat[slice] = self._slice_h2p(idat[slice])
+        return res

nifty/spaces/gl_space.py

@@ -112,7 +112,7 @@ class GLSpace(Space):
 
     def get_default_codomain(self):
         from .. import LMSpace
-        return LMSpace(lmax=self.nlat-1)
+        return LMSpace(lmax=self._nlat-1, mmax=self._nlon//2)
 
     def check_codomain(self, codomain):
         from .. import LMSpace

nifty/spaces/lm_space.py

@@ -33,6 +33,13 @@ class LMSpace(Space):
     lmax : int
         The maximum :math:`l` value of any spherical harmonics
         :math:`Y_{lm}` that is represented in this Space.
+        Must be >=0.
+
+    mmax : int *optional*
+        The maximum :math:`m` value of any spherical harmonics
+        :math:`Y_{lm}` that is represented in this Space.
+        If not supplied, it is set to lmax.
+        Must be >=0 and <=lmax.
 
     See Also
     --------
@@ -40,15 +47,6 @@ class LMSpace(Space):
     GLSpace : A class for the Gauss-Legendre discretization of the
         sphere [#]_.
 
-    Raises
-    ------
-    ValueError
-        If given lmax is negative.
-
-    Notes
-    -----
-        This implementation implicitly sets the mmax parameter to lmax.
-
     References
     ----------
     .. [#] K.M. Gorski et al., 2005, "HEALPix: A Framework for
@@ -59,13 +57,20 @@ class LMSpace(Space):
            `arXiv:1303.4945 <http://www.arxiv.org/abs/1303.4945>`_
     """
 
-    def __init__(self, lmax):
+    def __init__(self, lmax, mmax=None):
         super(LMSpace, self).__init__()
-        self._needed_for_hash += ["_lmax"]
-        self._lmax = self._parse_lmax(lmax)
+        self._needed_for_hash += ["_lmax", "_mmax"]
+        self._lmax = np.int(lmax)
+        if self._lmax < 0:
+            raise ValueError("lmax must be >=0.")
+        if mmax is None:
+            mmax = self._lmax
+        self._mmax = np.int(mmax)
+        if self._mmax < 0 or self._mmax > self._lmax:
+            raise ValueError("mmax must be >=0 and <=lmax.")
 
     def __repr__(self):
-        return ("LMSpace(lmax=%r)" % self.lmax)
+        return ("LMSpace(lmax=%r, mmax=%r)" % (self.lmax, self.mmax))
 
     @property
     def harmonic(self):
@@ -77,28 +82,27 @@ class LMSpace(Space):
 
     @property
     def dim(self):
-        l = self.lmax
+        l = self._lmax
+        m = self._mmax
         # the LMSpace consists of the full triangle (including -m's!),
         # minus two little triangles if mmax < lmax
-        # dim = (((2*(l+1)-1)+1)**2/4 - 2 * (l-m)(l-m+1)/2
-        # dim = np.int((l+1)**2 - (l-m)*(l-m+1.))
-        # We fix l == m
-        return np.int((l+1)*(l+1))
+        return (l+1)**2 - (l-m)*(l-m+1)
 
     def scalar_dvol(self):
         return 1.
 
     def get_k_length_array(self):
-        lmax = self.lmax
+        lmax = self._lmax
+        mmax = self._mmax
         ldist = np.empty((self.dim,), dtype=np.float64)
         ldist[0:lmax+1] = np.arange(lmax+1, dtype=np.float64)
         tmp = np.empty((2*lmax+2), dtype=np.float64)
         tmp[0::2] = np.arange(lmax+1)
         tmp[1::2] = np.arange(lmax+1)
         idx = lmax+1
-        for l in range(1, lmax+1):
-            ldist[idx:idx+2*(lmax+1-l)] = tmp[2*l:]
-            idx += 2*(lmax+1-l)
+        for m in range(1, mmax+1):
+            ldist[idx:idx+2*(lmax+1-m)] = tmp[2*m:]
+            idx += 2*(lmax+1-m)
         return Field((self,), from_np(ldist))
 
     def get_unique_k_lengths(self):
@@ -106,6 +110,9 @@ class LMSpace(Space):
 
     @staticmethod
     def _kernel(x, sigma):
+        # cf. "All-sky convolution for polarimetry experiments"
+        # by Challinor et al.
+        # http://arxiv.org/abs/astro-ph/0008228
         res = x+1.
         res *= x
         res *= -0.5*sigma*sigma
@@ -113,9 +120,6 @@ class LMSpace(Space):
         return res
 
     def get_fft_smoothing_kernel_function(self, sigma):
-        # cf. "All-sky convolution for polarimetry experiments"
-        # by Challinor et al.
-        # http://arxiv.org/abs/astro-ph/0008228
         return lambda x: self._kernel(x, sigma)
 
     @property
@@ -129,15 +133,8 @@ class LMSpace(Space):
     def mmax(self):
         """ Returns the maximum :math:`m` value of any spherical harmonic
         :math:`Y_{lm}` that is represented in this Space.
-        Currently this is identical to lmax.
         """
-        return self._lmax
-
-    def _parse_lmax(self, lmax):
-        lmax = np.int(lmax)
-        if lmax < 0:
-            raise ValueError("lmax must be >=0.")
-        return lmax
+        return self._mmax
 
     def get_default_codomain(self):
         from .. import GLSpace
@@ -146,4 +143,4 @@ class LMSpace(Space):
     def check_codomain(self, codomain):
         from .. import GLSpace, HPSpace
         if not isinstance(codomain, (GLSpace, HPSpace)):
-            raise TypeError("codomain must be a GLSpace.")
+            raise TypeError("codomain must be a GLSpace or HPSpace.")

test/test_spaces/test_lm_space.py

@@ -27,21 +27,21 @@ from nifty2go.dobj import to_ndarray as to_np
 
 # [lmax, expected]
 CONSTRUCTOR_CONFIGS = [
-        [5, {
+        [5, None, {
             'lmax': 5,
             'mmax': 5,
             'shape': (36,),
             'harmonic': True,
             'dim': 36,
             }],
-        [7, {
+        [7, 4, {
             'lmax': 7,
-            'mmax': 7,
-            'shape': (64,),
+            'mmax': 4,
+            'shape': (52,),
             'harmonic': True,
-            'dim': 64,
+            'dim': 52,
             }],
-        [-1, {
+        [-1, 28, {
             'error': ValueError
             }]
     ]
@@ -72,18 +72,18 @@ class LMSpaceInterfaceTests(unittest.TestCase):
             ['mmax', int],
             ['dim', int]])
     def test_property_ret_type(self, attribute, expected_type):
-        l = ift.LMSpace(7)
+        l = ift.LMSpace(7, 5)
         assert_(isinstance(getattr(l, attribute), expected_type))
 
 
 class LMSpaceFunctionalityTests(unittest.TestCase):
     @expand(CONSTRUCTOR_CONFIGS)
-    def test_constructor(self, lmax, expected):
+    def test_constructor(self, lmax, mmax, expected):
         if 'error' in expected:
             with assert_raises(expected['error']):
-                ift.LMSpace(lmax)
+                ift.LMSpace(lmax, mmax)
         else:
-            l = ift.LMSpace(lmax)
+            l = ift.LMSpace(lmax, mmax)
             for key, value in expected.items():
                 assert_equal(getattr(l, key), value)