python numpy log_Cython教程与代码之——Cython_tutorial之cython和numpy以及pandas

前面我们探讨的都是一些比较简单的应用情况，下面来介绍一下我们更常用的优化场景。

一般在数据挖掘的工程中我们最常用的两个python库，numpy和pandas，下面我们来学习一下怎么使用cython来优化numpy和pandas的代码。

纯python版

%load_ext Cython
import numpy as np


def shannon_entropy_py(p_x):
    return - np.sum(p_x * np.log(p_x))

（补充一下这里的参数-a，很nice的使用参数，具体的可以显示代码中可优化的部分：

图中的黄色部分表示这部分代码都是全部或者部分用到了python的解释器，所以这部分代码是拖慢程序运行速度的罪魁祸首，我们需要在这部分上进行cython的改造。）

首先，最简单也最有效的改造方式，把python中涉及到变量的部分 用cython定义一遍，这里就要用到 cimport numpy as cnp了，numpy给cython留了调用的c-level的接口，使用cimport可以在cython中导入c模块从而通过调用c模块来加快程序运行的速度。

numpy真是个好东西！

我们这里把函数的参数进行类型声明，“cnp.ndarray p_x” （和c完全一样的变量声明方式）

%%cython 

import numpy as np
cimport numpy as cnp

def shannon_entropy_cy(cnp.ndarray p_x):
    return - np.sum(p_x * np.log(p_x))

比较一下二者的运行速度区别：

辣鸡，才提升了这么一点。

接下来在上一部的cython定义里面，我们把np.log改为从c库中导入（没错，使用cimport允许你在cython中导入标准的c函数）

%%cython 

cimport numpy as cnp
from libc.math cimport log as clog

def shannon_entropy_v1(cnp.ndarray p_x):
    cdef double res = 0.0
    cdef int n = p_x.shape[0]
    cdef int i
    for i in range(n):
        res += p_x[i] * clog(p_x[i])
    return -res

wtf？？？为什么速度慢了这么多？？？？

这里实际上涉及到矢量化的概念，矢量化的大概的含义就是将一些比较复杂的运算转化为张量（比如二维条件下的矩阵运算）的运算，最典型的例子就是通过矩阵乘的方法代替了for循环，比如，针对两个5*5的数组进行求和操作，如果使用for循环我们需要循环一大圈，但是从矩阵的角度出发就方便多了，对应位置的元素相加即可。虽然这二者听起来运算的目的是一样的，但是在cpu上运行的效率是不可比拟的。

为了比较，我们把循环的写法写一遍然后计算一下时间：

def shannon_entropy_v1(p_x):
    res=0.0
    n = p_x.shape[0]
    for i in range(n):
        res += p_x[i] * np.log(p_x[i])
    return -res
%%timeit
shannon_entropy_v1(pmf)

下面我们来使用神器的numpy数据类型的指定。相对于上面再上面的那个cython实现，下面的代码仅仅改变了cnp.ndarray[double]，从原来的cnp.ndarray改成了cnp.ndarray[double]，

%%cython 

cimport numpy as cnp
from libc.math cimport log as clog

def shannon_entropy_v2(cnp.ndarray[double] p_x):
    cdef double res = 0.0
    cdef int n = p_x.shape[0]
    cdef int i
    for i in range(n):
        res += p_x[i] * clog(p_x[i])
    return -res

可以看到，速度比矢量化的numpy版还要快5倍左右！相对于python版本的循环函数快了161.73倍。

在上面的基础上，我们继续努力！

%%cython -a

cimport cython
cimport numpy as cnp
from libc.math cimport log

@cython.boundscheck(False)
@cython.wraparound(False)
def shannon_entropy_v3(cnp.ndarray[double] p_x):
    cdef double res = 0.0
    cdef int n = p_x.shape[0]
    cdef int i
    for i in range(n):
        res += p_x[i] * log(p_x[i])
    return -res

%%timeit
shannon_entropy_v3(pmf)

可以看到，在取消了边界检查和环绕检查之后速度又提升了一丢丢。

我们再来看看有没有办法针对于循环做一个简单的优化，重新生成一个更大的pmf来比较：

from scipy.stats import poisson
poi = poisson(10.0)
n = 10000000
pmf = poi.pmf(np.arange(n))

原来的最快的方式耗时：

我们针对循环体做了一些改变，具体参照：

https://stackoverflow.com/questions/21382180/cython-pure-c-loop-optimization​stackoverflow.com

%%cython -a

cimport cython
cimport numpy as cnp
from libc.math cimport log

@cython.boundscheck(False)
@cython.wraparound(False)
def shannon_entropy_v3(cnp.ndarray[double] p_x):
    cdef double res = 0.0
    cdef int n = p_x.shape[0]
    cdef int i
    for i from 0<=i<n:
        res += p_x[i] * log(p_x[i])
    return -res

可以看到,改变了循环体的结构不使用python的range之后运行的效率又提升了不少。具体原因可以参考上面的stack overflow的讨论。

对于pandas，tutorial给出的建议是，用numpy代替pandas，因为pandas是在numpy的基础上进行的高级的封装，所以，不可避免的要带来一定的性能的损失。因此建议在使用pandas的场合都替换为numpy，一方面使用numpy在大部分场景下都可以直接加速，一方面也便于扩展，julia、cython、numba都是可以直接和numpy进行交互的。

---

介绍完了函数，该介绍一些cython中的类如何来定义了。其实也挺简单的。

下面是纯python的类定义。

class PyLCG(object):
    
    def __init__(self, a=1664525, c=1013904223, m=2**3, seed=0):
        self.a = a
        self.c = c
        if m <= 0:
            raise ValueError("m must be > 0, given {}".format(m))
        self.m = m
        # The RNG state.
        self.x = seed
        
    def _advance(self):
        r = self.x
        self.x = (self.a * self.x + self.c) % self.m
        return r
        
    def randint(self, size=None):
        if size is None:
            return self._advance()
        return np.asarray([self._advance() for _ in range(size)])

时间花销为6.21ms

然后是cython下的类定义：可以看到比较特别的地方在于：

（1）cython下的类在初始化__init__的时候用的是__cinit__

（2）__cinit__中的参数要提前在_cinit_之间进行声明否则类无法直接识别

（3）针对与如何在cython中定义list的问题，google上给出很多解决方案，不过对于存放纯数字的list来说，numpy的c api不失为一个很方便的选择,或者如果有可能，在cython中使用numcpp？（最快的速度应该是使用c的malloc来开辟动态数组，或者c++的vector等方式，当然如果通过list来append dict，tuple之类的会比较麻烦这里暂不赘述，因为我还没有搞清楚，c和c++的语法知识我需要足够的时间来慢慢研究），我觉得list、dict、tuple这些python的数据类型要在cython中通过纯c或c++的方式表达出来也是比较难，这一块后续单独写一篇总结的文章放到当前的专栏下，因为涉及到的东西实在太多了。

%%cython -a

import numpy as np
cimport numpy as cnp
cimport cython

# Creates a new extension type: https://docs.python.org/3/extending/newtypes.html
cdef class CyLCG:
    
    # We declare the compile-time types of our *instance* attributes here.
    # This is similar to C++ class declaration syntax.
    cdef long a, c, m, x
    
    
    # Special Cython-defined initializer.
    # Called before __init__ to initialize all C-level attributes.
    def __cinit__(self, long a=1664525, long c=1013904223, long m=2**3, long seed=0):
        self.a = a
        self.c = c
        if m <= 0:
            raise ValueError("m must be > 0, given {}".format(m))
        self.m = m
        self.x = seed
    
    # cdef / cpdef methods are supported
    @cython.cdivision(True)
    cdef long _advance(self):
        cdef long r = self.x
        self.x = (self.a * self.x + self.c) % self.m
        return r
    
    # Regular def method
    @cython.boundscheck(False)
    @cython.wraparound(False)
    def randint(self, size=None):
        cdef long r
        if size is None:
            # Call to self._advance() here is efficient and at the C level.
            r = self._advance()
            return r
        cdef long[:] a = np.empty((size,), dtype='i4')
        cdef int i
        cdef int n = int(size)
        for i in range(n):
            a[i] = self._advance()
        return np.asarray(a)

速度提高了42.53倍。

除此之后，使用cython定义的类的内存占用也要远小于python下的类：

### Pure-python内存占用

*`PyLCG（）`对象本身

*实例`__dict__`

*和实例`__dict__`中的每个键/值

python下的类确实多了不少属性，后面会在《python高性能编程》的专栏里添加一篇关于python中的数据的属性的删减的问题。

import sys
(sys.getsizeof(pyrng) # the object itself
 + sys.getsizeof(pyrng.__dict__)  # the instance __dict__
 + sum(sys.getsizeof(k) + sys.getsizeof(v) for k, v in pyrng.__dict__.items())) # k/v memory use

而cython下定义的类型仅包含了这个类本身所占用的内存。

(sys.getsizeof(CyLCG()) # the object itself
 + 4 * 8) # The 4 8-byte longs (a, c, m, x)

可以看到两个类属性上的差距。（CyLCG中cdef的_advance在python中是不可见的）

同时根据上图我们可以知道，凡是在cython定义的类中用cdef定义的方法和属性，在python中都是不可见的，

同时也不可以给这个类创造新的属性（而这一切在python中都是可以实现的）

这个时候我们可以使用cython中的 cdef public和cdef readonly来避免这个问题(很遗憾只能使得属性变得可见，cdef的方法无法变得可见只能通过别的手段比如cpdef或者是def中的嵌套cdef的方式)

%%cython -a

import numpy as np
cimport cython

cdef class CyLCGOpen:
    
    cdef public long x
    cdef readonly long a, c, m
    
    def __cinit__(self, long a=1664525, long c=1013904223, long m=2**32, long seed=0):
        self.a = a
        self.c = c
        if m <= 0:
            raise ValueError("m must be > 0, given {}".format(m))
        self.m = m
        self.x = seed
        
    # cdef / cpdef methods are supported
    @cython.cdivision(True)        
    cdef  long _advance(self):
        cdef long r = self.x
        self.x = (self.a * self.x + self.c) % self.m
        return r
    
    # Regular def method
    @cython.boundscheck(False)
    @cython.wraparound(False)
    def randint(self, size=None):
        cdef long r
        if size is None:
            # Call to self._advance() here is efficient and at the C level.
            r = self._advance()
            return r
        cdef long[::1] a = np.empty((size,), dtype='i8')
        cdef int i
        cdef int n = int(size)
        for i in range(n):
            a[i] = self._advance()
        return np.asarray(a)