scipy与numpy中很多函数用c实现。
本篇以最简单的convolve为例做一个相应的动态库。
1.[转载] Python 调用 C/C++(基础篇)
这种做法称为Python扩展。比如说,我们有一个功能强大的C函数:期望在Python里这样使用:int great_function(int a) { return a + 1; }考虑最简单的情况。我们把功能强大的函数放入C文件 great_module.c 中。>>> from great_module import great_function >>> great_function(2) 3除了功能强大的函数great_function外,这个文件中还有以下部分:#include <Python.h> int great_function(int a) { return a + 1; } static PyObject * _great_function(PyObject *self, PyObject *args) { int _a; int res; if (!PyArg_ParseTuple(args, "i", &_a)) return NULL; res = great_function(_a); return PyLong_FromLong(res); } static PyMethodDef GreateModuleMethods[] = { { "great_function", _great_function, METH_VARARGS, "" }, {NULL, NULL, 0, NULL} }; PyMODINIT_FUNC initgreat_module(void) { (void) Py_InitModule("great_module", GreateModuleMethods); }
- 包裹函数_great_function。它负责将Python的参数转化为C的参数(PyArg_ParseTuple),调用实际的great_function,并处理great_function的返回值,最终返回给Python环境。
- 导出表GreateModuleMethods。它负责告诉Python这个模块里有哪些函数可以被Python调用。导出表的名字可以随便起,每一项有4个参数:第一个参数是提供给Python环境的函数名称,第二个参数是_great_function,即包裹函数。第三个参数的含义是参数变长,第四个参数是一个说明性的字符串。导出表总是以{NULL, NULL, 0, NULL}结束。
- 导出函数initgreat_module。这个的名字不是任取的,是你的module名称添加前缀init。导出函数中将模块名称与导出表进行连接。
在Windows下面,在Visual Studio命令提示符下编译这个文件的命令是
cl /LD great_module.c /o great_module.pyd -IC:Python27include C:Python27libspython27.lib/LD 即生成动态链接库。编译成功后在当前目录可以得到 great_module.pyd(实际上是dll)。这个pyd可以在Python环境下直接当作module使用。
在Linux下面,则用gcc编译:
gcc -fPIC -shared great_module.c -o great_module.so -I/usr/include/python2.7/ -lpython2.7在当前目录下得到great_module.so,同理可以在Python中直接使用。
测试结果如下
2.[转载] Python调用C/C++(使用SWIG)
用C/C++对脚本语言的功能扩展是非常常见的事情,Python也不例外。除了SWIG,市面上还有若干用于Python扩展的工具包,比较知名的还有Boost.Python、SIP等,此外,Cython由于可以直接集成C/C++代码,并方便的生成Python模块,故也可以完成扩展Python的任务。答主在这里选用SWIG的一个重要原因是,它不仅可以用于Python,也可以用于其他语言。如今SWIG已经支持C/C++的好基友Java,主流脚本语言Python、Perl、Ruby、PHP、JavaScript、tcl、Lua,还有Go、C#,以及R。SWIG是基于配置的,也就是说,原则上一套配置改变不同的编译方法就能适用各种语言(当然,这是理想情况了……)SWIG的安装方便,有Windows的预编译包,解压即用,绿色健康。主流Linux通常集成swig的包,也可以下载源代码自己编译,SWIG非常小巧,通常安装不会出什么问题。还记得本文第2节的那个great_function吗?有了SWIG,事情就会变得如此简单:
/* great_module.i */ %module great_module %{ int great_function(int a) { return a + 1; } %} int great_function(int a);换句话说,SWIG自动完成了诸如Python类型转换、module初始化、导出代码表生成的诸多工作。
对于C++,SWIG也可以应对。例如以下代码有C++类的定义:
//great_class.h #ifndef GREAT_CLASS #define GREAT_CLASS class Great { private: int s; public: void setWall (int _s) {s = _s;}; int getWall () {return s;}; }; #endif // GREAT_CLASS对应的SWIG配置文件
/* great_class.i */ %module great_class %{ #include "great_class.h" %} %include "great_class.h"这里不再重新敲一遍class的定义了,直接使用SWIG的%include指令
SWIG编译时要加-c++这个选项,生成的扩展名为cxx
Windows下编译:swig -c++ -python great_class.icl /LD great_class_wrap.cxx /o _great_class.pyd -IC:Python27include C:Python27libspython27.libLinux,使用C++的编译器
在Python交互模式下测试:g++ -fPIC -shared great_class_wrap.cxx -o _great_class.so -I/usr/include/python2.7/ -lpython2.7也就是说C++的class会直接映射到Python classSWIG非常强大,对于Python接口而言,简单类型,甚至指针,都无需人工干涉即可自动转换,而复杂类型,尤其是自定义类型,SWIG提供了typemap供转换。而一旦使用了typemap,配置文件将不再在各个语言当中通用。>>> import great_class >>> c = great_class.Great() >>> c.setWall(5) >>> c.getWall() 5
3.Python调用C/C++(convolve函数,使用SWIG)
由上可见,使用swig相比第一种方法要简单很多。
需要实现的函数在Python中的代码是
def conv1d(a,b): m=len(a) n=len(b) t = m + n - 1 r = [0] * t for i in range(t): for j in range(min(m,i+1)): if (i - j < n): r[i] += a[j] * b[i-j] return r
使用swig时遇到几个问题:
1. 输入的Python数组无法被c/c++识别
这里采用
%typemap(in) int[ANY](int temp[$1_dim0])
将参数中的PySequence转化为c中的数组。
$1_dim0指的是数组传入大小,如下例中为20,PyObject_Length($input)为输入数据中PyList的长度。
通过一个循环,将每个数据用PySequence_GetItem($input,i)获取并传递给 $1作为转化后的参数。
2. 输出的c数组指针无法被Python识别
这里采用
%typemap(out) int* conv1d
将返回值c数组转化为PyList。
通过一个循环,将每个数据用PyList_SetItem($input,i)赋值并传递给 $result作为转化后的返回值。
3.将c数组转化为PyList时,不知道长度。
这里采用
%typemap(in,numinputs=0,noblock=1) size_t *len
将长度传递给上一步中的typemap。
numinputs=0 指不需要数据输入。
Typemaps支持一个名为noblock的特殊属性,其中可以使用{...}分隔符,但是分隔符实际上不会生成到代码中。 效果类似于使用“”或%{%}分隔符,但代码通过预处理器运行。
%module t %include <stdint.i> %typemap(in,numinputs=0,noblock=1) size_t *len { size_t templen; $1 = &templen; } %typemap(in) int[ANY](int temp[$1_dim0]) { int i; if (!PySequence_Check($input)) { PyErr_SetString(PyExc_TypeError,"Expecting a sequence"); return NULL; } if (PyObject_Length($input) > $1_dim0) { PyErr_SetString(PyExc_ValueError,"Expecting a sequence with less elements"); return NULL; } int j=PyObject_Length($input); for (i =0; i < j; i++) { PyObject *o = PySequence_GetItem($input,i); if (!PyInt_Check(o)) { Py_XDECREF(o); PyErr_SetString(PyExc_ValueError,"Expecting a sequence of ints"); return NULL; } temp[i] = PyInt_AsLong(o); Py_DECREF(o); } $1 = &temp[0]; } %typemap(out) int* conv1d { int i; $result = PyList_New(templen); for (i = 0; i < templen; i++) { PyObject *o = PyInt_FromLong((long)$1[i]); PyList_SetItem($result,i,o); } } %inline %{ int *conv1d(int a[20],int b[20],int n, int m, size_t *len) { int j,i; const int t = n+m-1; int *r = (int *) malloc((t)*sizeof(int)); for (i = 0; i < t; i++) for (r[i]=0,j = 0; j < m && j <= i; j++) if (i - j < n) r[i] += a[j] * b[i - j]; *len = t; return r; } %}
4.Python调用C/C++(convolve函数,使用Cython)
以下介绍引用自用Cython来提高Python代码速度 [一]
什么是Cython?
Cython并不是C,也不是Python,他是一个将介于C和Python之间的语言结构编译成C,然后在Python中运行的接口。它保持了大多数Python的语法结构,支持所有Python的命令,并通过一定方法加入C的结构,从而绕过Python的重重转译和类型检查,来实现快速的程序。
大多数Python整合包都会包括Cython。你也可以通过pip来安装:
pip install Cython与Python不一样,你并不能直接使用Cython的源代码。需要通过编译来完成从Python调用C程序的目的。此外,没有错,如果你有用C直接编写的程序,可以用Cython来调用纯C程序!
Cython同样支持C++
Cython自带一些C和C++的标准库
ipython支持直接在命令行编译和使用cython程序!然而我的实践发现可能出现奇怪bug。
Cython的参考书非常少,最常用的一本是Cython A GUIDE FOR PYTHON PROGRAMMERS。然而……此书非常不详尽,绝对不是noob friendly…… 建议有什么问题还是stackoverflow求助大牛。
什么时候用Cython?
切记:Cython并不是所有python速度问题的良药。它牺牲了flexibility来换取速度,你必须要考虑所有细节问题,包括很多python新手甚至老手忽略了的、Python帮你暗地里处理好的问题。
-
如果你需要大量使用python的类,请慎重考虑是否使用。
-
如果并没有大量循环,或者最内层的循环有一些复杂python的类,请慎重使用。
-
如果不需要重复使用脚本,请慎重使用。
-
如果提升只有不到20%,请慎重使用。
-
如果这个函数或者类会被很多人debug和改进,请慎重使用。除非很多人熟悉Cython,不然maintain程序会非常麻烦。
-
如果不能确定性能瓶颈是你要改写为cython的函数,请慎重使用。
-
如果你不会从200%的速度提升中获得快感,请慎重使用。
根据上面的需求,我们制作了如下pyx文件。
import numpy as np cimport cython cimport numpy as np @cython.boundscheck(False) @cython.wraparound(False) def _cython_convolve(np.ndarray[np.int_t, ndim=1] a, np.ndarray[np.int_t, ndim=1] b): cdef int m = a.shape[0] cdef int n = b.shape[0] cdef int i,j,t = m + n - 1 cdef np.ndarray[np.int_t, ndim=1] r = np.zeros(t, dtype=np.int) for i in range(t): for j in range(min(m,i+1)): if (i - j < n): r[i] += a[j] * b[i-j] return r def cython_convolve(a, b): return _cython_convolve(np.array(a), np.array(b))
以及如下的setup.py文件
# setup.py from distutils.core import setup, Extension from Cython.Build import cythonize import numpy setup(ext_modules = cythonize(Extension( 'cython_convolve', sources=['cython_convolve.pyx'], language='c', include_dirs=[numpy.get_include()], library_dirs=[], libraries=[], extra_compile_args=[], extra_link_args=[] )))
在终端运行
python setup.py build_ext --inplace
生成的cython_convolve.so即可用于调用。
5.在Python中的使用效果
自定义的conv1d最慢,使用swig实现的最快。使用cython实现的比numpy自带的略慢。
import scipy.signal.signaltools as sg import numpy as np import timeit import t import cython_convolve a = [[1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4]] b = [[1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4]] def conv1d(a,b): m=len(a) n=len(b) t = m + n - 1 r = [0] * t for i in range(t): for j in range(min(m,i+1)): if (i - j < n): r[i] += a[j] * b[i-j] return r def npconv(): np.convolve(a[0],b[0]) def sgconv(): sg.convolve2d(a, b) def newconv(): conv1d(a[0], b[0]) def tconv(): t.conv1d(a[0],b[0],16,16) def cconv(): cython_convolve.cython_convolve(a[0], b[0]) if __name__=='__main__': print sg.convolve2d(a, b) print np.convolve(a[0], b[0]) print conv1d(a[0], b[0]) print t.conv1d(a[0],b[0],16,16) print cython_convolve.cython_convolve(a[0], b[0]) print timeit.timeit('sgconv()',setup='from __main__ import sgconv',number=100000) print timeit.timeit('npconv()',setup='from __main__ import npconv',number=100000) print timeit.timeit('newconv()',setup='from __main__ import newconv',number=100000) print timeit.timeit('tconv()', setup='from __main__ import tconv', number=100000) print timeit.timeit('cconv()', setup='from __main__ import cconv', number=100000) """ [[ 1 4 10 20 27 32 36 40 53 60 62 60 79 88 88 80 103 108 94 60 77 80 68 40 51 52 42 20 25 24 16]] [ 1 4 10 20 27 32 36 40 53 60 62 60 79 88 88 80 103 108 94 60 77 80 68 40 51 52 42 20 25 24 16] [1, 4, 10, 20, 27, 32, 36, 40, 53, 60, 62, 60, 79, 88, 88, 80, 103, 108, 94, 60, 77, 80, 68, 40, 51, 52, 42, 20, 25, 24, 16] [1, 4, 10, 20, 27, 32, 36, 40, 53, 60, 62, 60, 79, 88, 88, 80, 103, 108, 94, 60, 77, 80, 68, 40, 51, 52, 42, 20, 25, 24, 16] [ 1 4 10 20 27 32 36 40 53 60 62 60 79 88 88 80 103 108 94 60 77 80 68 40 51 52 42 20 25 24 16] 1.34557986259 0.578319072723 5.11206579208 0.260015964508 0.591851968765 """