重點(diǎn)介紹廣播機(jī)制以及針對(duì)高維數(shù)組的軸操作，最后對(duì) NumPy 的 C 語(yǔ)言擴(kuò)展作了介紹。

1廣播

NumPy 運(yùn)算通常是在兩個(gè)數(shù)組的元素級(jí)別上進(jìn)行的。最簡(jiǎn)單情況就是，兩個(gè)具有完全相同 shape 的數(shù)組運(yùn)算，如下面例子所示，

a = np.array([1.0, 2.0, 3.0])
b = np.array([2.0, 2.0, 2.0])
a * b

numpy 的廣播機(jī)制是指在執(zhí)行算術(shù)運(yùn)算時(shí)處理不同 shape 的數(shù)組的方式。在一定規(guī)則下，較小的數(shù)組在較大的數(shù)組上廣播，從而使得數(shù)組具有兼容的 shape。

a = np.array([1.0, 2.0, 3.0])
b = 2.0
a * b

發(fā)現(xiàn)這兩個(gè)計(jì)算的結(jié)果是一樣的，但第二個(gè)是有廣播機(jī)制在發(fā)揮作用。

廣播規(guī)則

在兩個(gè)數(shù)組上執(zhí)行運(yùn)算時(shí)，NumPy 比較它們的形狀。它從 shape 的最右邊開(kāi)始往左一一比較。如果所有位子比較下來(lái)都是下面兩種情況之一，

• 相同位子上的兩個(gè)數(shù)字相等
• 或者其中之一是 1

那么這兩個(gè)數(shù)組可以運(yùn)算。如果不滿足這些條件，則將引發(fā) ValueError，表明數(shù)組的 shape 不兼容。

可見(jiàn)，數(shù)組的 shape 好比人的八字，兩個(gè)人如果八字不合，那是不能在一起滴。

在下面這些示例中，A 和 B 數(shù)組中長(zhǎng)度為 1 的那些軸（缺失的軸自動(dòng)補(bǔ) 1），在廣播期間會(huì)擴(kuò)展為另一個(gè)數(shù)組相同位子上更大的長(zhǎng)度，

A      (3d array):  15 x 3 x 5
B      (3d array):  15 x 1 x 5
Result (3d array):  15 x 3 x 5

A      (3d array):  15 x 3 x 5
B      (2d array):       3 x 5
Result (3d array):  15 x 3 x 5

A      (3d array):  15 x 3 x 5
B      (2d array):       3 x 1
Result (3d array):  15 x 3 x 5
    
A      (4d array):  8 x 1 x 6 x 1
B      (3d array):      7 x 1 x 5
Result (4d array):  8 x 7 x 6 x 5

下面例子中第一個(gè)數(shù)組的 shape 為 (3,3)，第二個(gè)數(shù)組的 shape 為 (3,)，此時(shí)相當(dāng)于 (1,3)，因此先將第二個(gè)數(shù)組的 shape 改為 (3,3)，相當(dāng)于原來(lái)數(shù)組沿著 0 軸再?gòu)?fù)制 2 份。

 

MatA = np.array([[1, 2, 3],[4,5,6],[7,8,9]])
MatB = np.array([1, 2, 3])
MatA + MatB

為了更好地理解這個(gè)機(jī)制，下面再給出幾個(gè)例子。下圖共三行，分別對(duì)應(yīng)三種廣播方式，請(qǐng)對(duì)照后面代碼。

 

a = np.array([0,10,20,30])
b = np.array([0,1,2])

A = np.stack((a,a,a), axis=1)

B = np.stack((b,b,b,b))

# 對(duì)應(yīng)第一種情況
A + B

# 對(duì)應(yīng)第二種情況
A + b

a1 = np.array([[0,10,20,30]]).T

# 對(duì)應(yīng)第三種情況
a1 + b

而下面例子不滿足廣播規(guī)則，因而不能執(zhí)行運(yùn)算。

A      (1d array):  3
B      (1d array):  4 # 倒數(shù)最后的軸長(zhǎng)度不兼容

A      (2d array):  4 x 3
B      (1d array):      4 # 倒數(shù)最后的軸長(zhǎng)度不兼容
    
A      (2d array):      2 x 1
B      (3d array):  8 x 4 x 3 # 倒數(shù)第二個(gè)軸長(zhǎng)度不兼容

 

廣播機(jī)制小結(jié)

• 廣播機(jī)制為數(shù)組運(yùn)算提供了一種便捷方式。
• 話雖如此，它并非在所有情況下都有效，并且實(shí)際上強(qiáng)加了執(zhí)行廣播必須滿足的嚴(yán)格規(guī)則。
• 僅當(dāng)數(shù)組中每個(gè)維的形狀相等或維的大小為 1 時(shí)，才能執(zhí)行算術(shù)運(yùn)算。

2維度增減

維度增加

在需要增加軸的位子使用 np.newaxis 或者 None。

x = np.arange(6).reshape(2,3)
x, x.shape

x1 = x[:,np.newaxis,:]
x1, x1.shape

# 或者
x2 = x[:,None,:]
x2, x2.shape

維度壓縮

有時(shí)候需要將數(shù)組中多余的軸去掉，以降低數(shù)組維度的目的。

numpy.squeeze( )

• 從數(shù)組中刪除單維度的軸，即把 shape 中為 1 的維度去掉。

x = np.arange(6).reshape(2,1,3)

y = x.squeeze()

xd = x.__array_interface__['data'][0] 
yd = y.__array_interface__['data'][0] 

 

3數(shù)組轉(zhuǎn)置（換軸）

x = np.arange(9).reshape(3, 3)

y = np.transpose(x) # 或者 y = x.transpose() 或者 x.T

y = np.transpose(x, [1, 0])

x = np.array([3,2,1,0,4,5,6,7,8,9,10,11,12,13,14,15]).reshape(2, 2, 4)

y1 = np.transpose(x, [1, 0, 2])

請(qǐng)對(duì)照下圖理解這個(gè)三維數(shù)組在內(nèi)存中的樣子以及對(duì)它的不同視圖（view）。關(guān)于這點(diǎn)，文末附上的進(jìn)階篇有詳細(xì)解讀。

 

y2 = np.transpose(x, [2, 0, 1])

# 代碼放一起
x = np.array([3,2,1,0,4,5,6,7,8,9,10,11,12,13,14,15]).reshape(2, 2, 4)
y0 = np.transpose(x, [1, 2, 0])
y1 = np.transpose(x, [1, 0, 2])
y2 = np.transpose(x, [2, 0, 1])

看看變軸后各個(gè)數(shù)組的元素具體是怎樣的，注意，它們都指向同一份數(shù)據(jù)。

 

這是怎么實(shí)現(xiàn)對(duì)內(nèi)存中同一份數(shù)據(jù)使用不同的軸序呢？實(shí)際上，數(shù)據(jù)還是那些數(shù)據(jù)，更改的是各個(gè)軸上的步長(zhǎng) stride。

x.strides, y1.strides, y2.strides

# 數(shù)據(jù)還是同一份
id(x.data), id(y1.data), id(y2.data)

再看一個(gè)例子，三維數(shù)組有三個(gè)軸，注意換軸后每個(gè)軸的步長(zhǎng)。

x = np.arange(16).reshape(2, 2, 4)

y = x.transpose((1, 0, 2))

 

兩個(gè)數(shù)組三個(gè)軸對(duì)應(yīng)的步長(zhǎng)不同了。

 

軸更換后，下標(biāo)也跟著換了，所以換軸前后相同下標(biāo)指向的數(shù)據(jù)是不同的。

 

其實(shí)，軸的意義主要體現(xiàn)在步長(zhǎng)上，所以換軸一定意義上就是更換了步長(zhǎng)。

實(shí)際例子

RGB 圖像數(shù)據(jù)

• 每張圖像由紅綠藍(lán)三個(gè)通道組成，每個(gè)通道對(duì)應(yīng)一個(gè)     32×32 的二維數(shù)組

 

看下圖，從左到右，分別對(duì)應(yīng)圖像數(shù)據(jù)在內(nèi)存中的存放，將一維數(shù)組轉(zhuǎn)化為三維數(shù)組，更換軸。

 

那么，為什么要換軸呢？因?yàn)椴煌绦虬鼘?duì)數(shù)據(jù)的要求不同，我們?yōu)榱耸褂盟鼈?，需要按照它們?duì)參數(shù)的要求來(lái)對(duì)數(shù)據(jù)作相應(yīng)調(diào)整。

而有時(shí)候，并不需要換軸，只需要更換某個(gè)軸上元素的次序即可，例如，

# 變換某個(gè)軸上元素的次序
z = x[..., (3, 2, 1, 0)]

4通用函數(shù)

ufunc 函數(shù)

• ufunc 是 universal function 的縮寫(xiě)，是指能對(duì)數(shù)組的每個(gè)元素進(jìn)行操作的函數(shù)，而不是針對(duì) narray 對(duì)象操作。
• NumPy 提供了大量的 ufunc 函數(shù)。這些函數(shù)在對(duì) narray 進(jìn)行運(yùn)算的速度比使用循環(huán)或者列表推導(dǎo)式要快得多。
• NumPy 內(nèi)置的許多 ufunc 函數(shù)是 C 語(yǔ)言實(shí)現(xiàn)的，因此計(jì)算效率很高。

x = np.linspace(0, 2*np.pi, 5)

y, z = np.sin(x), np.cos(x)

# 將結(jié)果直接傳給輸入 x
np.sin(x, x)

性能比較

import time
import math
import numpy as np

x = [i for i in range(1000000)]

# math.sin
start = time.process_time()
for i, t in enumerate(x):
    x[i] = math.sin(t)
math_time = time.process_time()-start    

# numpy.sin
x = np.array(x, dtype=np.float64)
start = time.process_time()
np.sin(x, x)
numpy_time = time.process_time()-start

# comparison
math_time, numpy_time, math_time/numpy_time

reduce 操作

• 這是 NumPy 內(nèi)置的通用函數(shù)，如果需要這樣的計(jì)算，建議直接使用，不要自己實(shí)現(xiàn)。
• 沿著軸對(duì)數(shù)組進(jìn)行操作，相當(dāng)于將運(yùn)算符            插入到沿軸的所有子數(shù)組或者元素當(dāng)中。
• 格式為：     .reduce (array=, axis=0, dtype=None)

np.add.reduce([1,2,3])

np.add.reduce([[1,2,3],[4,5,6]], axis=1)

np.multiply.reduce([[1,2,3],[4,5,6]], axis=1)

accumulate 操作

• 這也是 NumPy 內(nèi)置的通用函數(shù)，如果需要這樣的計(jì)算，建議直接使用，不要自己實(shí)現(xiàn)。

• 與 reduce 類似，只是它返回的數(shù)組和輸入的數(shù)組的 shape 相同，保存所有的中間計(jì)算結(jié)果。

np.add.accumulate([1,2,3])

np.add.accumulate([[1,2,3],[4,5,6]], axis=1)

自定義 ufunc 函數(shù)

# 定義一個(gè) python 函數(shù)
def ufunc_diy(x):
    c, c0, hc = 0.618, 0.518, 1.0
    x = x - int(x) 
    if x >= c: 
        r = 0.0
    elif x < c0: 
        r = x / c0 * hc
    else: 
        r = (c-x) / (c-c0) * hc
    return r

x = np.linspace(0, 2, 1000000)

ufunc_diy(x)

start = time.process_time()
y1 = np.array([ufunc_diy(t) for t in x])
time_1 = time.process_time()-start
time_1

np.frompyfunc 函數(shù)

• 將一個(gè)計(jì)算單個(gè)元素的函數(shù)轉(zhuǎn)換成 ufunc 函數(shù)

ufunc = np.frompyfunc(ufunc_diy, 1, 1)

start = time.process_time()
y2 = ufunc(x)
time_2 = time.process_time()-start
time_2

NumPy 之 C 擴(kuò)展

本文主要介紹兩種擴(kuò)展方式，

• ctypes
• Cython

ctypes

• ctypes 是 Python 的一個(gè)外部庫(kù)，提供和 C 語(yǔ)言兼容的數(shù)據(jù)類型，可以很方便地調(diào)用 dll/so 中輸出的 C 接口函數(shù)。

#ufunc.c
'''
void ufunc_diy(double x, double y, int size) {

    double xx,r,c=0.618,c0=0.518,hc=1.0;
    for(int i=0;i        xx = x[i]-(int)(x[i]);
        if (xx>=c) r=0.0;
        else if (xx        else r=(c-xx)/(c-c0)*hc;
        y[i]=r;
    }
}
'''

#ufunc.py
""" Example of wrapping a C library function that accepts a C double array as
    input using the numpy.ctypeslib. """

import numpy as np
import numpy.ctypeslib as npct
from ctypes import c_int

array_1d_double = npct.ndpointer(dtype=np.double, ndim=1, flags='CONTIGUOUS')

# load the library, using numpy mechanisms
lib = npct.load_library("lib_ufunc", ".")

# setup the return types and argument types
lib.ufunc_diy.restype = None
lib.ufunc_diy.argtypes = [array_1d_double, array_1d_double, c_int]

def ufunc_diy_func(in_array, out_array):
    return lib.ufunc_diy(in_array, out_array, len(in_array))

# 編譯
# gcc -shared -fPIC -O2 ufunc.c -ldl -o lib_ufunc.so

import time
import numpy as np
import ufunc

start = time.process_time() 
ufunc.ufunc_diy_func(x, x)
end = time.process_time()
print("ufunc_diy time: ", end-start)

# python test_ufunc.py 
# ufunc_diy time:  0.003 - 0.008

Cython

• Cython 是 Python 的一個(gè)超集，可以編譯為 C，Cython 結(jié)合了 Python 的易用性和原生 C 代碼的高效率。

# ufunc_diy.h
void ufunc_diy(double  in_array, double  out_array, int size);

# ufunc_diy.c
void ufunc_diy(double x, double y, int size) {

    double xx,r,c=0.618,c0=0.518,hc=1.0;
    for(int i=0;i        xx = x[i]-(int)(x[i]);
        if (xx>=c) r=0.0;
        else if (xx        else r=(c-xx)/(c-c0)*hc;
        y[i]=r;
    }
}

# Cython支持 NumPy
# 在代碼中聲明 a = np.array([0,10,20,30])
b = np.array([0,1,2])cimport numpy，使用函數(shù)。

#_ufunc_cython.pyx_
""" Example of wrapping a C function that takes C double arrays as input using
    the Numpy declarations from Cython """

# cimport the Cython declarations for numpy
cimport numpy as np

# if you want to use the Numpy-C-API from Cython
# (not strictly necessary for this example, but good practice)
np.import_array()

# cdefine the signature of our c function
cdef extern from "ufunc_diy.h":
    void ufunc_diy (double  in_array, double  out_array, int size)

# create the wrapper code, with numpy type annotations
def ufunc_diy_func(np.ndarray[double, ndim=1, mode="c"] in_array not Noa = np.array([0,10,20,30])
b = np.array([0,1,2])ne,
                   np.ndarray[double, ndim=1, mode="c"] out_array not None):
    ufunc_diy( np.PyArray_DATA(in_array),
               np.PyArray_DATA(out_array),
              in_array.shape[0])

# setup.py
from distutils.core import setup, Extension
import numpy
from Cython.Distutils import build_ext

setup(
    cmdclass={'build_ext': build_ext},
    ext_modules=[Extension("ufunc_cython",
                 sources=["_ufunc_cython.pyx", "ufunc_diy.c"],
                 include_dirs=[numpy.get_include()])],
)
# 或者
from distutils.core import setup
import numpy
from Cython.Build import cythonize
 
setup(
    ext_modules=cythonize("_ufunc_cython.pyx", annotate=True),
    include_dirs=[numpy.get_include()]
)

# 編譯
python setup.py build_ext --inplace

可以看到多了兩個(gè)文件，一個(gè)是 _ufunc_cython.c，一個(gè)是 ufunc_cython.so（如果是 windows，則是 .pyd）。

import time
import numpy as np
import ufunc_cython

start = time.process_time() 
ufunc_cython.ufunc_diy_func(x, x)
end = time.process_time()
print("ufunc_diy time: ", end-start)