From da193c1450cd7151745646c8062060e2a904287f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 19 Apr 2017 17:43:24 +0300 Subject: [PATCH] something already works Former-commit-id: 9b1cfe39ae1453d26db40ad36d649b820d091ee4 --- qpms/qpms_c.pyx | 161 ++++++++++++++++++++++++++++++++++++++++++-- qpms/translations.c | 11 ++- qpms/translations.h | 6 ++ setup.py | 18 ++++- 4 files changed, 189 insertions(+), 7 deletions(-) diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index e9c5e66..1a91542 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -1,6 +1,8 @@ -# cythonized functions here -cimport numpy as np +# Cythonized parts of QPMS here +# ----------------------------- + 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 = (data)#[0] + #cdef char *func_name= (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 = (func)( + (ip0)[0], + (ip1)[0], + (ip2)[0], + (ip3)[0], + (ip4)[0], + (ip5)[0], + (ip6)[0], + (ip7)[0], + (ip8)[0], + ) + (op0)[0] = 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] = qpms_trans_single_A_Taylor_ext +elementwise_funcs_B[0] = 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 + ) + diff --git a/qpms/translations.c b/qpms/translations.c index 8b35505..ffddc28 100644 --- a/qpms/translations.c +++ b/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); +} diff --git a/qpms/translations.h b/qpms/translations.h index 78cb325..2bbd7d2 100644 --- a/qpms/translations.h +++ b/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 diff --git a/setup.py b/setup.py index 9779442..7bb4b17 100644 --- a/setup.py +++ b/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'],