use iterative map analysis for Healpix

nifty/operators/fft_operator/transformations/hplmtransformation.py

@@ -88,11 +88,8 @@ class HPLMTransformation(SlicingTransformation):
         mmax = lmax
         nside= self.domain.nside
 
-        sjob = pyHealpix.sharpjob_d()
-        sjob.set_Healpix_geometry(nside)
-        sjob.set_triangular_alm_info(lmax, mmax)
         if issubclass(inp.dtype.type, np.complexfloating):
-            [resultReal, resultImag] = [sjob.map2alm(x)
+            [resultReal, resultImag] = [pyHealpix.map2alm_iter(x,lmax,mmax,3)
                                         for x in (inp.real, inp.imag)]
 
             [resultReal,
@@ -102,7 +99,7 @@ class HPLMTransformation(SlicingTransformation):
             result = self._combine_complex_result(resultReal, resultImag)
 
         else:
-            result = sjob.map2alm(inp)
+            result = pyHealpix.map2alm_iter(inp,lmax,mmax,3)
             result = lm_transformation_factory.buildIdx(result, lmax=lmax)
 
         return result

nifty/operators/fft_operator/transformations/lmhptransformation.py

@@ -64,7 +64,7 @@ class LMHPTransformation(SlicingTransformation):
         if not isinstance(domain, LMSpace):
             raise TypeError("domain needs to be a LMSpace.")
 
-        nside = np.max(domain.lmax//2, 1)
+        nside = max((domain.lmax+1) // 2, 1)
         result = HPSpace(nside=nside, dtype=domain.dtype)
         return result
 
@@ -89,21 +89,18 @@ class LMHPTransformation(SlicingTransformation):
         lmax = self.domain.lmax
         mmax = lmax
 
-        sjob = pyHealpix.sharpjob_d()
-        sjob.set_Healpix_geometry(nside)
-        sjob.set_triangular_alm_info(lmax, mmax)
         if issubclass(inp.dtype.type, np.complexfloating):
             [resultReal,
              resultImag] = [lm_transformation_factory.buildLm(x, lmax=lmax)
                             for x in (inp.real, inp.imag)]
 
-            [resultReal, resultImag] = [sjob.alm2map(x)
+            [resultReal, resultImag] = [pyHealpix.alm2map(x,lmax,mmax,nside)
                                         for x in [resultReal, resultImag]]
 
             result = self._combine_complex_result(resultReal, resultImag)
 
         else:
             result = lm_transformation_factory.buildLm(inp, lmax=lmax)
-            result = sjob.alm2map(result)
+            result = pyHealpix.alm2map(result,lmax,mmax,nside)
 
         return result

test/test_misc.py

@@ -31,6 +31,7 @@ from nifty import Field,\
                   FieldArray, \
                   RGRGTransformation, \
                   LMGLTransformation, \
+                  LMHPTransformation, \
                   FFTOperator
 
 from nose.plugins.skip import SkipTest
@@ -103,3 +104,15 @@ class Misc_Tests(unittest.TestCase):
                 inp = Field.from_random(domain=a,random_type='normal',std=7,mean=3,dtype=tp)
                 out = fft.inverse_times(fft.times(inp))
                 assert_allclose(inp.val, out.val)
+
+    def test_sht2(self):
+        if 'pyHealpix' not in di:
+            raise SkipTest
+        for lm in [128,256]:
+            for tp in [np.float64,np.complex128,np.float32,np.complex64]:
+                a = LMSpace(lmax=lm)
+                b = LMHPTransformation.get_codomain(a)
+                fft = FFTOperator(domain=a, target=b, domain_dtype=tp, target_dtype=tp)
+                inp = Field.from_random(domain=a,random_type='normal',std=1,mean=0,dtype=tp)
+                out = fft.inverse_times(fft.times(inp))
+                assert_allclose(inp.val, out.val,rtol=1e-3,atol=1e-3)