WIP: Further transformation adaptions

nifty/operators/fft_operator/transformations/gllmtransformation.py

@@ -6,7 +6,6 @@ import nifty.nifty_utilities as utilities
 from nifty import GLSpace, LMSpace
 import lm_transformation_factory as ltf
 
-hp = gdi.get('healpy')
 gl = gdi.get('libsharp_wrapper_gl')
 
 
@@ -119,10 +118,8 @@ class GLLMTransformation(Transformation):
                 inpImg = gl.map2alm(
                     np.imag(inp).astype(np.float64, copy=False), nlat=nlat,
                     nlon=nlon, lmax=lmax, mmax=mmax)
-                #TODO gl shouldn't depend on hp
-                lmaxArray, mmaxArray = hp.Alm.getlm(lmax=lmax)
-                inpReal = ltf.buildIdx(inpReal, lmaxArray, mmaxArray)
-                inpImg = ltf.buildIdx(inpImg, lmaxArray, mmaxArray)
+                inpReal = ltf.buildIdx(inpReal, lmax=lmax)
+                inpImg = ltf.buildIdx(inpImg, lmax=lmax)
                 inp = inpReal + inpImg * 1j
             else:
                 if self.domain.dtype == np.dtype('float32'):
@@ -133,6 +130,7 @@ class GLLMTransformation(Transformation):
                     inp = gl.map2alm(inp,
                                      nlat=nlat, nlon=nlon,
                                      lmax=lmax, mmax=mmax)
+                inp = ltf.buildIdx(inp, lmax=lmax)
 
             if slice_list == [slice(None, None)]:
                 return_val = inp

nifty/operators/fft_operator/transformations/hplmtransformation.py

@@ -117,15 +117,14 @@ class HPLMTransformation(Transformation):
                     np.imag(inp).astype(np.float64, copy=False),
                     lmax=lmax, mmax=mmax, iter=niter, pol=True,
                     use_weights=False, datapath=None)
-                lmaxArray, mmaxArray = hp.Alm.getlm(lmax=lmax)
-                inpReal = ltf.buildIdx(inpReal,lmaxArray, mmaxArray)
-                inpImg = ltf.buildIdx(inpImg,lmaxArray, mmaxArray)
+                inpReal = ltf.buildIdx(inpReal, lmax=lmax)
+                inpImg = ltf.buildIdx(inpImg, lmax=lmax)
                 inp = inpReal + inpImg * 1j
             else:
                 inp = hp.map2alm(inp.astype(np.float64, copy=False),
                                  lmax=lmax, mmax=mmax, iter=niter, pol=True,
                                  use_weights=False, datapath=None)
-
+                inp = ltf.buildIdx(inp, lmax=lmax)
             if slice_list == [slice(None, None)]:
                 return_val = inp
             else:

nifty/operators/fft_operator/transformations/lm_transformation_factory.pyx

 import numpy as np
 cimport numpy as np
 
-cpdef buildIdx(np.ndarray[np.complex128_t, ndim=1] nr, np.ndarray[np.int_t] l, np.ndarray[np.int_t] m, np.int_t lmax):
+cpdef buildIdx(np.ndarray[np.complex128_t, ndim=1] nr, np.int_t lmax):
     cdef np.int size = (lmax+1)*(lmax+1)
 
     cdef np.ndarray res=np.zeros([size], dtype=np.complex128)
     cdef np.ndarray final=np.zeros([size], dtype=np.float64)
     res[0:lmax+1] = nr[0:lmax+1]
 
-    cdef np.ndarray resL=np.zeros([size], dtype=np.int)
-    resL[0:lmax+1] = np.arange(lmax+1)
-
-    cdef np.ndarray resM=np.zeros([size], dtype=np.int)
-
     for i in xrange(len(nr)-lmax-1):
         res[i*2+lmax+1] = nr[i+lmax+1]
         res[i*2+lmax+2] = np.conjugate(nr[i+lmax+1])
-        resL[i*2+lmax+1] = l[i+lmax+1]
-        resL[i*2+lmax+2] = l[i+lmax+1]
-        resM[i*2+lmax+1] = m[i+lmax+1]
-        resM[i*2+lmax+2] = -m[i+lmax+1]
-    final = realify(res,resL, resM)
-    return final, resL, resM
+    final = realify(res, lmax)
+    return final
+
+cpdef np.ndarray[np.float64_t, ndim=1] realify(np.ndarray[np.complex128_t, ndim=1] nr, np.int_t lmax):
+    cdef np.ndarray resi=np.zeros([len(nr)], dtype=np.float64)
+
+    resi[0:lmax+1] = np.real(nr[0:lmax+1])
+
+    for i in xrange(lmax+1,len(nr),2):
+        # m calculation print i,(i-lmax)/2+1,(2*lmax+1)*(2*lmax+1)-4*(i-lmax)+1
+        mi =  int(np.ceil(((2*lmax+1)-np.sqrt((2*lmax+1)*(2*lmax+1)-4*(i-lmax)+1))/2))
+        resi[i]=np.sqrt(2)*np.real(nr[i])*(-1)**(mi*mi)
+        resi[i+1]=np.sqrt(2)*np.imag(nr[i])*(-1)**(mi*mi)
+    return resi
 
-cpdef buildLm(np.ndarray[np.float64_t, ndim=1] nr, np.ndarray[np.int_t] l, np.ndarray[np.int_t] m, np.int_t lmax):
+cpdef buildLm(np.ndarray[np.float64_t, ndim=1] nr, np.int_t lmax):
     cdef np.int size = (len(nr)-lmax-1)/2+lmax+1
 
     cdef np.ndarray res=np.zeros([size], dtype=np.complex128)
     cdef np.ndarray temp=np.zeros([len(nr)], dtype=np.complex128)
     res[0:lmax+1] = nr[0:lmax+1]
 
-    cdef np.ndarray resL=np.zeros([size], dtype=np.int)
-    resL[0:lmax+1] = np.arange(lmax+1)
-
-    cdef np.ndarray resM=np.zeros([size], dtype=np.int)
-
-    temp = inverseRealify(nr, l, m)
+    temp = inverseRealify(nr, lmax)
 
     for i in xrange(0,len(temp)-lmax-1,2):
         res[i/2+lmax+1] = temp[i+lmax+1]
-        resL[i/2+lmax+1] = l[i+lmax+1]
-        resM[i/2+lmax+1] = m[i+lmax+1]
-    return res,resL,resM
-
-cpdef np.ndarray[np.float64_t, ndim=1] realify(np.ndarray[np.complex128_t, ndim=1] nr, np.ndarray[np.int_t] l, np.ndarray[np.int_t] m):
-    cdef np.ndarray resi=np.zeros([len(nr)], dtype=np.float64)
-
-    for i in xrange(len(nr)):
-        if m[i]<0:
-            resi[i]=np.sqrt(2)*np.imag(nr[i-1])*(-1)**(m[i]*m[i])
-        elif m[i]>0:
-            resi[i]=np.sqrt(2)*np.real(nr[i])*(-1)**(m[i]*m[i])
-        else:
-            resi[i]=np.real(nr[i])
-    return resi
-
+    return res
 
-cpdef np.ndarray[np.complex128_t, ndim=1] inverseRealify(np.ndarray[np.float64_t, ndim=1] nr, np.ndarray[np.int_t] l, np.ndarray[np.int_t] m):
+cpdef np.ndarray[np.complex128_t, ndim=1] inverseRealify(np.ndarray[np.float64_t, ndim=1] nr, np.int_t lmax):
     cdef np.ndarray resi=np.zeros([len(nr)], dtype=np.complex128)
+    resi[0:lmax+1] = np.real(nr[0:lmax+1])
 
-    for i in xrange(len(nr)):
-        if m[i]<0:
-            resi[i]=1/np.sqrt(2)*(nr[i-1]-1j*nr[i])
-        elif m[i]>0:
-            resi[i]=(-1)**m[i]/np.sqrt(2)*(nr[i]+1j*nr[i+1])
-        else:
-            resi[i]=np.real(nr[i])
+    for i in xrange(lmax+1,len(nr),2):
+        # m calculation print i,(i-lmax)/2+1,(2*lmax+1)*(2*lmax+1)-4*(i-lmax)+1
+        mi =  int(np.ceil(((2*lmax+1)-np.sqrt((2*lmax+1)*(2*lmax+1)-4*(i-lmax)+1))/2))
+        resi[i]=(-1)**mi/np.sqrt(2)*(nr[i]+1j*nr[i+1])
+        resi[i+1]=1/np.sqrt(2)*(nr[i]-1j*nr[i+1])
     return resi
 

nifty/operators/fft_operator/transformations/lmgltransformation.py

@@ -5,6 +5,7 @@ from nifty.config import dependency_injector as gdi
 import nifty.nifty_utilities as utilities
 from nifty import GLSpace, LMSpace
 
+import lm_transformation_factory as ltf
 gl = gdi.get('libsharp_wrapper_gl')
 
 
@@ -117,12 +118,25 @@ class LMGLTransformation(Transformation):
             lmax = self.domain.lmax
             mmax = self.mmax
 
-            if self.domain.dtype == np.dtype('complex64'):
-                inp = gl.alm2map_f(inp, nlat=nlat, nlon=nlon,
-                                   lmax=lmax, mmax=mmax, cl=False)
-            else:
-                inp = gl.alm2map(inp, nlat=nlat, nlon=nlon,
+            if self.domain.dtype == np.dtype('complex128'):
+                inpReal = np.real(inp)
+                inpImag = np.imag(inp)
+                inpReal = ltf.buildLm(inpReal,lmax=lmax)
+                inpImag = ltf.buildLm(inpImag,lmax=lmax)
+                inpReal = gl.alm2map(inpReal, nlat=nlat, nlon=nlon,
+                                 lmax=lmax, mmax=mmax, cl=False)
+                inpImag = gl.alm2map(inpImag, nlat=nlat, nlon=nlon,
                                  lmax=lmax, mmax=mmax, cl=False)
+                inp = inpReal+inpImag*(1j)
+            else:
+                inp = ltf.buildLm(inp, lmax=lmax)
+
+                if self.domain.dtype == np.dtype('complex64'):
+                    inp = gl.alm2map_f(inp, nlat=nlat, nlon=nlon,
+                                       lmax=lmax, mmax=mmax, cl=False)
+                else:
+                    inp = gl.alm2map(inp, nlat=nlat, nlon=nlon,
+                                     lmax=lmax, mmax=mmax, cl=False)
 
             if slice_list == [slice(None, None)]:
                 return_val = inp

nifty/operators/fft_operator/transformations/lmhptransformation.py

@@ -5,6 +5,8 @@ from nifty.config import dependency_injector as gdi
 import nifty.nifty_utilities as utilities
 from nifty import HPSpace, LMSpace
 
+import lm_transformation_factory as ltf
+
 hp = gdi.get('healpy')
 
 
@@ -102,10 +104,26 @@ class LMHPTransformation(Transformation):
             lmax = self.domain.lmax
             mmax = self.domain.mmax
 
-            inp = inp.astype(np.complex128, copy=False)
-            inp = hp.alm2map(inp, nside, lmax=lmax, mmax=mmax,
-                             pixwin=False, fwhm=0.0, sigma=None,
-                             pol=True, inplace=False)
+            if self.domain.dtype == np.dtype('complex128'):
+                inpReal = np.real(inp)
+                inpImag = np.imag(inp)
+                inpReal = ltf.buildLm(inpReal,lmax=lmax)
+                inpImag = ltf.buildLm(inpImag,lmax=lmax)
+                inpReal = hp.alm2map(inpReal.astype(np.complex128, copy=False), nside,
+                                 lmax=lmax, mmax=mmax,
+                                 pixwin=False, fwhm=0.0, sigma=None,
+                                 pol=True, inplace=False)
+                inpImag = hp.alm2map(inpImag.astype(np.complex128, copy=False), nside,
+                                 lmax=lmax, mmax=mmax,
+                                 pixwin=False, fwhm=0.0, sigma=None,
+                                 pol=True, inplace=False)
+                inp = inpReal+inpImag*(1j)
+            else:
+                inp = ltf.buildLm(inp, lmax=lmax)
+                inp = hp.alm2map(inp.astype(np.complex128, copy=False), nside,
+                                 lmax=lmax, mmax=mmax,
+                                 pixwin=False, fwhm=0.0, sigma=None,
+                                 pol=True, inplace=False)
 
             if slice_list == [slice(None, None)]:
                 return_val = inp