重點介紹廣播機制以及針對高維數組的軸操作，最后對 NumPy 的 C 語言擴展作了介紹。

1廣播

NumPy 運算通常是在兩個數組的元素級別上進行的。最簡單情況就是，兩個具有完全相同 shape 的數組運算，如下面例子所示，

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

numpy 的廣播機制是指在執行算術運算時處理不同 shape 的數組的方式。在一定規則下，較小的數組在較大的數組上廣播，從而使得數組具有兼容的 shape。

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

發現這兩個計算的結果是一樣的，但第二個是有廣播機制在發揮作用。

廣播規則

在兩個數組上執行運算時，NumPy 比較它們的形狀。它從 shape 的最右邊開始往左一一比較。如果所有位子比較下來都是下面兩種情況之一，

• 相同位子上的兩個數字相等
• 或者其中之一是 1

那么這兩個數組可以運算。如果不滿足這些條件，則將引發 ValueError，表明數組的 shape 不兼容。

可見，數組的 shape 好比人的八字，兩個人如果八字不合，那是不能在一起滴。

在下面這些示例中，A 和 B 數組中長度為 1 的那些軸（缺失的軸自動補 1），在廣播期間會擴展為另一個數組相同位子上更大的長度，

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

下面例子中第一個數組的 shape 為 (3,3)，第二個數組的 shape 為 (3,)，此時相當于 (1,3)，因此先將第二個數組的 shape 改為 (3,3)，相當于原來數組沿著 0 軸再復制 2 份。

 

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

為了更好地理解這個機制，下面再給出幾個例子。下圖共三行，分別對應三種廣播方式，請對照后面代碼。

 

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))

# 對應第一種情況
A + B

# 對應第二種情況
A + b

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

# 對應第三種情況
a1 + b

而下面例子不滿足廣播規則，因而不能執行運算。

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

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

 

廣播機制小結

• 廣播機制為數組運算提供了一種便捷方式。
• 話雖如此，它并非在所有情況下都有效，并且實際上強加了執行廣播必須滿足的嚴格規則。
• 僅當數組中每個維的形狀相等或維的大小為 1 時，才能執行算術運算。

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

維度壓縮

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

numpy.squeeze( )

• 從數組中刪除單維度的軸，即把 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數組轉置（換軸）

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])

請對照下圖理解這個三維數組在內存中的樣子以及對它的不同視圖（view）。關于這點，文末附上的進階篇有詳細解讀。

 

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])

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

 

這是怎么實現對內存中同一份數據使用不同的軸序呢？實際上，數據還是那些數據，更改的是各個軸上的步長 stride。

x.strides, y1.strides, y2.strides

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

再看一個例子，三維數組有三個軸，注意換軸后每個軸的步長。

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

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

 

兩個數組三個軸對應的步長不同了。

 

軸更換后，下標也跟著換了，所以換軸前后相同下標指向的數據是不同的。

 

其實，軸的意義主要體現在步長上，所以換軸一定意義上就是更換了步長。

實際例子

RGB 圖像數據

• 每張圖像由紅綠藍三個通道組成，每個通道對應一個     32×32 的二維數組

 

看下圖，從左到右，分別對應圖像數據在內存中的存放，將一維數組轉化為三維數組，更換軸。

 

那么，為什么要換軸呢？因為不同程序包對數據的要求不同，我們為了使用它們，需要按照它們對參數的要求來對數據作相應調整。

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

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

4通用函數

ufunc 函數

• ufunc 是 universal function 的縮寫，是指能對數組的每個元素進行操作的函數，而不是針對 narray 對象操作。
• NumPy 提供了大量的 ufunc 函數。這些函數在對 narray 進行運算的速度比使用循環或者列表推導式要快得多。
• NumPy 內置的許多 ufunc 函數是 C 語言實現的，因此計算效率很高。

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

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

# 將結果直接傳給輸入 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 內置的通用函數，如果需要這樣的計算，建議直接使用，不要自己實現。
• 沿著軸對數組進行操作，相當于將運算符       <op>      插入到沿軸的所有子數組或者元素當中。
• 格式為：     <op>.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 內置的通用函數，如果需要這樣的計算，建議直接使用，不要自己實現。

• 與 reduce 類似，只是它返回的數組和輸入的數組的 shape 相同，保存所有的中間計算結果。

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

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

自定義 ufunc 函數

# 定義一個 python 函數
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 函數

• 將一個計算單個元素的函數轉換成 ufunc 函數

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

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

NumPy 之 C 擴展

本文主要介紹兩種擴展方式，

• ctypes
• Cython

ctypes

• ctypes 是 Python 的一個外部庫，提供和 C 語言兼容的數據類型，可以很方便地調用 dll/so 中輸出的 C 接口函數。

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

    double xx,r,c=0.618,c0=0.518,hc=1.0;
    for(int i=0;i<size;i++) {
        xx = x[i]-(int)(x[i]);
        if (xx>=c) r=0.0;
        else if (xx<c0) r=xx/c0*hc;
        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 的一個超集，可以編譯為 C，Cython 結合了 Python 的易用性和原生 C 代碼的高效率。

# ufunc_diy.h
void ufunc_diy(double <em> in_array, double </em> out_array, int size);

# ufunc_diy.c
void ufunc_diy(double <em>x, double </em>y, int size) {

    double xx,r,c=0.618,c0=0.518,hc=1.0;
    for(int i=0;i<size;i++) {
        xx = x[i]-(int)(x[i]);
        if (xx>=c) r=0.0;
        else if (xx<c0) r=xx/c0*hc;
        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，使用函數。

#_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 <em> in_array, double </em> 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(<double*> np.PyArray_DATA(in_array),
              <double*> 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

可以看到多了兩個文件，一個是 _ufunc_cython.c，一個是 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)