Taking care of types

nifty/operators/fft_operator/transformations/gllmtransformation.py

@@ -111,23 +111,18 @@ class GLLMTransformation(Transformation):
                     return_val = np.empty_like(temp_val)
                 inp = temp_val[slice_list]
 
-            if self.domain.dtype == np.dtype('complex128'):
-                inpReal = gl.map2alm(
+            if inp.dtype >= np.dtype('complex64'):
+                inpReal = self.GlMap2Alm(
                     np.real(inp).astype(np.float64, copy=False), nlat=nlat,
                     nlon=nlon, lmax=lmax, mmax=mmax)
-                inpImg = gl.map2alm(
+                inpImg = self.GlMap2Alm(
                     np.imag(inp).astype(np.float64, copy=False), nlat=nlat,
                     nlon=nlon, lmax=lmax, mmax=mmax)
                 inpReal = ltf.buildIdx(inpReal, lmax=lmax)
                 inpImg = ltf.buildIdx(inpImg, lmax=lmax)
                 inp = inpReal + inpImg * 1j
             else:
-                if self.domain.dtype == np.dtype('float32'):
-                    inp = gl.map2alm_f(inp,
-                                       nlat=nlat, nlon=nlon,
-                                       lmax=lmax, mmax=mmax)
-                else:
-                    inp = gl.map2alm(inp,
+                inp = self.GlMap2Alm(inp,
                                      nlat=nlat, nlon=nlon,
                                      lmax=lmax, mmax=mmax)
                 inp = ltf.buildIdx(inp, lmax=lmax)
@@ -144,3 +139,9 @@ class GLLMTransformation(Transformation):
             return_val = return_val.astype(self.codomain.dtype, copy=False)
 
         return return_val
+
+    def GlMap2Alm(self, inp, **kwargs):
+        if inp.dtype == np.dtype('float32'):
+            return gl.map2alm_f(inp, kwargs)
+        else:
+            return gl.map.alm(inp, kwargs)

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.int_t lmax):
+def buildLm(inp, **kwargs):
+    if inp.dtype == np.dtype('float32'):
+        return _buildLm_f(inp, **kwargs)
+    else:
+        return _buildLm(inp, **kwargs)
+
+def buildIdx(inp, **kwargs):
+    if inp.dtype == np.dtype('complex64'):
+        return _buildIdx_f(inp, **kwargs)
+    else:
+        return _buildIdx(inp, **kwargs)
+
+cpdef np.ndarray[np.float32_t, ndim=1]  _buildIdx_f(np.ndarray[np
+.complex64_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.complex64)
+    cdef np.ndarray final=np.zeros([size], dtype=np.float32)
+    res[0:lmax+1] = nr[0:lmax+1]
+
+    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])
+    final = _realify_f(res, lmax)
+    return final
+
+cpdef np.ndarray[np.float32_t, ndim=1] _realify_f(np.ndarray[np.complex64_t, ndim=1] nr, np.int_t lmax):
+    cdef np.ndarray resi=np.zeros([len(nr)], dtype=np.float32)
+
+    resi[0:lmax+1] = np.real(nr[0:lmax+1])
+
+    for i in xrange(lmax+1,len(nr),2):
+        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 np.ndarray[np.float64_t, ndim=1]  _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)
@@ -11,40 +49,64 @@ cpdef buildIdx(np.ndarray[np.complex128_t, ndim=1] nr, np.int_t lmax):
     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])
-    final = realify(res, lmax)
+    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):
+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.int_t lmax):
+cpdef np.ndarray[np.complex64_t, ndim=1] _buildLm_f(np.ndarray[np.float32_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.complex64)
+    cdef np.ndarray temp=np.zeros([len(nr)], dtype=np.complex64)
+    res[0:lmax+1] = nr[0:lmax+1]
+
+    temp = _inverseRealify_f(nr, lmax)
+
+    for i in xrange(0,len(temp)-lmax-1,2):
+        res[i/2+lmax+1] = temp[i+lmax+1]
+    return res
+
+cpdef np.ndarray[np.complex64_t, ndim=1] _inverseRealify_f(np.ndarray[np.float32_t, ndim=1] nr, np.int_t lmax):
+    cdef np.ndarray resi=np.zeros([len(nr)], dtype=np.complex64)
+    resi[0:lmax+1] = np.real(nr[0:lmax+1])
+
+    for i in xrange(lmax+1,len(nr),2):
+        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
+
+
+cpdef np.ndarray[np.complex128_t, ndim=1] _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]
 
-    temp = inverseRealify(nr, lmax)
+    temp = _inverseRealify(nr, lmax)
 
     for i in xrange(0,len(temp)-lmax-1,2):
         res[i/2+lmax+1] = temp[i+lmax+1]
     return res
 
-cpdef np.ndarray[np.complex128_t, ndim=1] inverseRealify(np.ndarray[np.float64_t, ndim=1] nr, np.int_t lmax):
+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(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])

nifty/operators/fft_operator/transformations/lmgltransformation.py

@@ -118,25 +118,20 @@ class LMGLTransformation(Transformation):
             lmax = self.domain.lmax
             mmax = self.mmax
 
-            if self.domain.dtype == np.dtype('complex128'):
+            if inp.dtype >= np.dtype('complex64'):
                 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,
+                inpReal = self.GlAlm2Map(inpReal, nlat=nlat, nlon=nlon,
                                  lmax=lmax, mmax=mmax, cl=False)
-                inpImag = gl.alm2map(inpImag, nlat=nlat, nlon=nlon,
+                inpImag = self.GlAlm2Map(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)
+                inp = self.GlAlm2Map(inp, nlat=nlat, nlon=nlon,
+                                   lmax=lmax, mmax=mmax, cl=False)
 
             if slice_list == [slice(None, None)]:
                 return_val = inp
@@ -154,3 +149,9 @@ class LMGLTransformation(Transformation):
             return_val = return_val.astype(self.codomain.dtype, copy=False)
 
         return return_val
+
+    def GlAlm2Map(self, inp, **kwargs):
+        if inp.dtype == np.dtype('complex64'):
+            return gl.alm2map_f(inp, kwargs)
+        else:
+            return gl.alm2map(inp, kwargs)