something already works

Former-commit-id: 9b1cfe39ae1453d26db40ad36d649b820d091ee4

qpms/qpms_c.pyx

@@ -1,6 +1,8 @@
-# cythonized functions here
+# Cythonized parts of QPMS here
-cimport numpy as np
+# -----------------------------
+
 import numpy as np
+cimport numpy as np
 cimport cython
 
 ## Auxillary function for retrieving the "meshgrid-like" indices; inc. nmax
@@ -68,7 +70,158 @@ def get_y_mn_unsigned(int nmax):
     return(ymn_plus, ymn_minus)
 
 cdef int q_max(int m, int n, int mu, int nu):
-    return min(n,nu,(n+nu-abs(m+mu)//2)
+    return min(n,nu,(n+nu-abs(m+mu)//2))
+
+"""
+Now we generate our own universal functions to be used with numpy.
+
+Good way to see how this is done is to look at scipy/scipy/special/generate_ufuncs.py
+and scipy/scipy/special/generate_ufuncs.py
+
+In simple words, it works like this:
+- Let's have a single element function. This can be function which returns or a "subroutine".
+- Then we need a loop function; this is a wrapper that gets bunch of pointers from numpy and
+  has to properly call the single element function.
+- From those two, we build a python object using PyUFunc_FromFuncAndData.
+  * If the ufunc is supposed to work on different kinds of input/output types,
+    then a pair of single-element and loop functions is o be provided for
+    each combination of types. However, the single-element function can be reused if
+    the corresponding loop functions do the proper casting.
+"""
+
+## as in scipy/special/_ufuncs_cxx.pyx
+##-------------------------------------
+#cdef extern from "numpy/ufuncobject.h":
+#    int PyUFunc_getfperr() nogil
+#
+#cdef public int wrap_PyUFunc_getfperr() nogil:
+#    """
+#    Call PyUFunc_getfperr in a context where PyUFunc_API array is initialized;
+#
+#    """
+#    return PyUFunc_getfperr()
+#
+#cimport sf_error
+#-------------------------------------
 
 
-#cdef translation_coefficient_A(
+ctypedef double complex cdouble
+
+cdef void loop_D_iiiidddii_As_D_lllldddbl(char **args, np.npy_intp *dims, np.npy_intp *steps, void *data) nogil:
+    cdef np.npy_intp i, n = dims[0]
+    cdef void *func = (<void**>data)#[0]
+    #cdef char *func_name= <char*>(<void**>data)[1] # i am not using this, nor have I saved func_name to data
+    cdef char *ip0 = args[0]
+    cdef char *ip1 = args[1]
+    cdef char *ip2 = args[2]
+    cdef char *ip3 = args[3]
+    cdef char *ip4 = args[4]
+    cdef char *ip5 = args[5]
+    cdef char *ip6 = args[6]
+    cdef char *ip7 = args[7]
+    cdef char *ip8 = args[8]
+    cdef char *op0 = args[9]
+    cdef cdouble ov0
+    for i in range(n): # iterating over dimensions
+        ov0 = (<double complex(*)(int, int, int, int, double, double, double, int, int) nogil>func)(
+            <int>(<np.npy_long*>ip0)[0],
+            <int>(<np.npy_long*>ip1)[0],
+            <int>(<np.npy_long*>ip2)[0],
+            <int>(<np.npy_long*>ip3)[0],
+            <double>(<np.npy_double*>ip4)[0],
+            <double>(<np.npy_double*>ip5)[0],
+            <double>(<np.npy_double*>ip6)[0],
+            <int>(<np.npy_bool*>ip7)[0],
+            <int>(<np.npy_long*>ip8)[0],
+        )
+        (<double *>op0)[0] = <double>ov0
+        ip0 += steps[0]
+        ip1 += steps[1]
+        ip2 += steps[2]
+        ip3 += steps[3]
+        ip4 += steps[4]
+        ip5 += steps[5]
+        ip6 += steps[6]
+        ip7 += steps[7]
+        ip8 += steps[8]
+        op0 += steps[9]
+# FIXME ERROR HANDLING!!! requires correct import and different data passed (see scipy's generated ufuncs)
+#    sf_error.check_fpe(func_name)
+
+
+#cdef extern from "numpy/arrayobject.h":
+#    cdef enum NPY_TYPES:
+#        NPY_DOUBLE
+#        NPY_CDOUBLE # complex double
+#        NPY_LONG # int
+#    ctypedef int npy_intp
+
+
+cdef extern from "translations.h":
+    cdouble qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu,
+        double r, double th, double ph, int r_ge_d, int J) nogil
+    cdouble qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu,
+        double r, double th, double ph, int r_ge_d, int J) nogil
+
+# Module initialisation
+# ---------------------
+
+np.import_array()  # not sure whether this is really needed
+np.import_ufunc()
+
+# Arrays passed to PyUFunc_FromFuncAndData()
+# ------------------------------------------
+
+# BTW, aren't there anonymous arrays in cython?
+
+cdef np.PyUFuncGenericFunction loop_func[1]
+cdef char input_output_types[10]
+cdef void *elementwise_funcs_A[1]
+cdef void *elementwise_funcs_B[1]
+
+loop_func[0] = loop_D_iiiidddii_As_D_lllldddbl
+
+input_output_types[0] = np.NPY_LONG
+input_output_types[1] = np.NPY_LONG
+input_output_types[2] = np.NPY_LONG
+input_output_types[3] = np.NPY_LONG
+input_output_types[4] = np.NPY_DOUBLE
+input_output_types[5] = np.NPY_DOUBLE
+input_output_types[6] = np.NPY_DOUBLE
+input_output_types[7] = np.NPY_BOOL
+input_output_types[8] = np.NPY_LONG
+input_output_types[9] = np.NPY_CDOUBLE
+
+elementwise_funcs_A[0] = <void*> qpms_trans_single_A_Taylor_ext
+elementwise_funcs_B[0] = <void*> qpms_trans_single_B_Taylor_ext
+
+trans_A_Taylor = np.PyUFunc_FromFuncAndData(
+        loop_func, # func
+        elementwise_funcs_A, #data 
+        input_output_types, # types
+        1, # ntypes: number of supported input types
+        9, # nin: number of input args
+        1, # nout: number of output args
+        0, # identity element, unused
+        "trans_A_Taylor", # name
+        """
+        TODO computes the E-E or M-M translation coefficient in Taylor's normalisation
+        """, # doc
+        0 # unused, for backward compatibility of numpy c api
+        )
+
+trans_B_Taylor = np.PyUFunc_FromFuncAndData(
+        loop_func,
+        elementwise_funcs_B,
+        input_output_types,
+        1, # number of supported input types
+        9, # number of input args
+        1, # number of output args
+        0, # identity element, unused
+        "trans_B_Taylor",
+        """
+        TODO computes the E-E or M-M translation coefficient in Taylor's normalisation
+        """,
+        0 # unused
+        )
+

qpms/translations.c

@@ -160,5 +160,14 @@ complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kd
     return (presum / prenormratio) * sum;
 }
 
+complex double qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu, 
+		double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) {
+	sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi};
+	return qpms_trans_single_A_Taylor(m,n,mu,nu,kdlj,r_ge_d,J);
+}
 
-
+complex double qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu, 
+		double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J) {
+	sph_t kdlj = {kdlj_r, kdlj_theta, kdlj_phi};
+	return qpms_trans_single_B_Taylor(m,n,mu,nu,kdlj,r_ge_d,J);
+}

qpms/translations.h

@@ -23,4 +23,10 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd
 complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
                 bool r_ge_d, qpms_bessel_t J);
 
+complex double qpms_trans_single_A_Taylor_ext(int m, int n, int mu, int nu, double kdlj_r,
+		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
+
+complex double qpms_trans_single_B_Taylor_ext(int m, int n, int mu, int nu, double kdlj_r,
+		double kdlj_th, double kdlj_phi, int r_ge_d, int J);
+
 #endif // QPMS_TRANSLATIONS_H

setup.py

@@ -8,11 +8,25 @@ from distutils.extension import Extension
 #        m = sys.modules['setuptools.extension']
 #        m.Extension.__dict__ = m._Extension.__dict__
 
+# TODO CHECK THIS OUT http://stackoverflow.com/questions/4056657/what-is-the-easiest-way-to-make-an-optional-c-extension-for-a-python-package
+# also this: https://docs.python.org/2/extending/building.html
+import os
+
+print("You might want to add additional library path to LD_LIBRARY_PATH (especially if you are not using"
+        " GNU GSL in your system library path) and if import fails. ")
+print(os.environ['LD_LIBRARY_PATH'].split(':'))
+
 qpms_c = Extension('qpms_c',
-        sources = ['qpms/qpms_c.pyx'])
+        sources = ['qpms/qpms_c.pyx','qpms/gaunt.c',#'qpms/gaunt.h','qpms/vectors.h','qpms/translations.h',
+            # FIXME http://stackoverflow.com/questions/4259170/python-setup-script-extensions-how-do-you-include-a-h-file
+            'qpms/translations.c'],
+        extra_compile_args=['-std=c99','-ggdb','-O3'],
+        libraries=['gsl', 'blas'],
+        runtime_library_dirs=os.environ['LD_LIBRARY_PATH'].split(':')
+        )
 
 setup(name='qpms',
-        version = "0.1.9",
+        version = "0.2.0",
         packages=['qpms'],
 #        setup_requires=['setuptools_cython'],
         install_requires=['cython>=0.21','quaternion','spherical_functions','py_gmm'],