Merge branch 'mmm2' into 'master'

Martin's monster merge part 2/N: partial de-cythonization

See merge request !65

nifty/operators/fft_operator/transformations/lm_transformation_factory.py

+import numpy as np
+
+
+def buildLm(nr, lmax):
+    new_dtype = np.result_type(nr.dtype, np.complex64)
+
+    size = (len(nr)-lmax-1)/2+lmax+1
+    res = np.zeros([size], dtype=new_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):
+    if nr.dtype == np.dtype('complex64'):
+        new_dtype = np.float32
+    elif nr.dtype == np.dtype('complex128'):
+        new_dtype = np.float64
+    elif nr.dtype == np.dtype('complex256'):
+        new_dtype = np.float128
+    else:
+        raise TypeError("dtype of nr not supported.")
+
+    size = (lmax+1)*(lmax+1)
+    final = np.zeros([size], dtype=new_dtype)
+    final[0:lmax+1] = nr[0:lmax+1].real
+    final[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
+    final[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
+    return final

nifty/operators/fft_operator/transformations/lm_transformation_factory.pyx

-
-import numpy as np
-cimport numpy as np
-
-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 final=np.zeros([size], dtype=np.float32)
-    final[0:lmax+1] = nr[0:lmax+1].real
-    final[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
-    final[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
-    return final
-
-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 final=np.zeros([size], dtype=np.float64)
-    final[0:lmax+1] = nr[0:lmax+1].real
-    final[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
-    final[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
-    return final
-
-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)
-    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
-
-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)
-    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

nifty/spaces/lm_space/lm_helper.pyx

-import numpy as np
-cimport numpy as np
-
-cpdef _distance_array_helper(np.ndarray[np.int_t] index_array, np.int_t lmax):
-    cdef np.int size = (lmax + 1) * (lmax  + 1)
-
-    cdef np.ndarray res = np.zeros([len(index_array)], dtype=np.float128)
-
-    for idx, index_full in enumerate(index_array):
-        if index_full <= lmax:
-            index_half = index_full
-        else:
-            if (index_full - lmax) % 2 == 0:
-                index_half = (index_full + lmax) / 2;
-            else:
-                index_half = (index_full + lmax + 1) / 2;
-
-        m = (np.ceil(((2 * lmax + 1) -
-                      np.sqrt((2 * lmax + 1)**2 -
-                              8 * (index_half - lmax))) / 2)).astype(int)
-        res[idx] = index_half - m * (2 * lmax + 1 - m) // 2
-
-    return res

nifty/spaces/lm_space/lm_space.py

@@ -25,8 +25,6 @@ from nifty.spaces.space import Space
 from nifty.config import nifty_configuration as gc,\
                          dependency_injector as gdi
 
-from lm_helper import _distance_array_helper
-
 from d2o import arange
 
 gl = gdi.get('libsharp_wrapper_gl')
@@ -177,11 +175,19 @@ class LMSpace(Space):
                        distribution_strategy=distribution_strategy)
 
         dists = dists.apply_scalar_function(
-            lambda x: _distance_array_helper(x, self.lmax),
+            lambda x: self._distance_array_helper(x, self.lmax),
             dtype=np.float)
 
         return dists
 
+    @staticmethod
+    def _distance_array_helper(index_array, lmax):
+        u = 2*lmax + 1
+        index_half = (index_array+np.minimum(lmax, index_array)+1)//2
+        m = (np.ceil((u - np.sqrt(u*u - 8*(index_half - lmax)))/2)).astype(int)
+        res = (index_half - m*(u - m)//2).astype(np.float64)
+        return res
+
     def get_fft_smoothing_kernel_function(self, sigma):
         if sigma is None:
             sigma = np.sqrt(2) * np.pi / (self.lmax + 1)

setup.py

@@ -17,7 +17,6 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 from setuptools import setup, find_packages
-import os
 
 from Cython.Build import cythonize
 import numpy
@@ -35,9 +34,6 @@ setup(name="ift_nifty",
       package_dir={"nifty": "nifty"},
       zip_safe=False,
       ext_modules=cythonize([
-          "nifty/operators/fft_operator/transformations/"
-              "lm_transformation_factory.pyx",
-          "nifty/spaces/lm_space/lm_helper.pyx",
           "nifty/operators/smoothing_operator/smooth_util.pyx"
       ]),
       include_dirs=[numpy.get_include()],