• Python 调用 C/C++实现卷积


    scipy与numpy中很多函数用c实现。

    本篇以最简单的convolve为例做一个相应的动态库。

    1.[转载] Python 调用 C/C++(基础篇)

    下面转载自Jerry Jho在知乎“如何实现 C/C++ 与 Python 的通信?”问题下的回复。
    示例简单,说明清晰,在我电脑上测试无误。

    这种做法称为Python扩展。比如说,我们有一个功能强大的C函数:
    int great_function(int a) {
        return a + 1;
    }
    期望在Python里这样使用:
    >>> from great_module import great_function 
    >>> great_function(2)
    3
    考虑最简单的情况。我们把功能强大的函数放入C文件 great_module.c 中。
    #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外,这个文件中还有以下部分:
    • 包裹函数_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)

    下面转载自Jerry Jho在知乎“如何实现 C/C++ 与 Python 的通信?”问题下的回复。
    示例简单,说明清晰,在我电脑上测试无误。
    用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

    swig -c++ -python great_class.i
    Windows下编译:
    cl /LD great_class_wrap.cxx /o _great_class.pyd -IC:Python27include C:Python27libspython27.lib

    Linux,使用C++的编译器

    g++ -fPIC -shared great_class_wrap.cxx -o _great_class.so  -I/usr/include/python2.7/ -lpython2.7
    在Python交互模式下测试:
    >>> import great_class
    >>> c = great_class.Great()
    >>> c.setWall(5)
    >>> c.getWall()
    5
    也就是说C++的class会直接映射到Python classSWIG非常强大,对于Python接口而言,简单类型,甚至指针,都无需人工干涉即可自动转换,而复杂类型,尤其是自定义类型,SWIG提供了typemap供转换。而一旦使用了typemap,配置文件将不再在各个语言当中通用。

    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
    """
  • 相关阅读:
    oracle将blob转为varchar2
    根据给定年份和周数获取指定周的开始结束日期
    获取指定月份的第一个周五
    找回误删的表和数据
    oracle大数据库
    Java数据库编程
    流的使用
    (转)支持 PS/2 与 USB 的键盘过滤驱动(可卸载)
    Unity编辑器扩展-Custom List, displaying data your way
    值得推荐的C/C++框架和库 (真的很强大)〔转〕
  • 原文地址:https://www.cnblogs.com/qw12/p/6227691.html
Copyright © 2020-2023  润新知