NumPy与数组优化

概述

本章讲解Cython与NumPy的深度集成。NumPy数组是Python科学计算的基础，Cython通过typed memoryview提供接近C的性能访问NumPy数据。

学习路径：NumPy基础 → typed memoryview → 数组操作 → 高级索引 → C互操作

核心概念：

• typed memoryview：零拷贝访问NumPy数组
• cnp.ndarray：C级NumPy数组声明
• 广播机制：标量与数组的自动扩展

---

8.1 NumPy基础回顾

ndarray结构

功能说明：NumPy数组的核心属性：shape、dtype、strides。

# NumPy数组创建
import numpy as np

arr = np.array([1, 2, 3, 4, 5], dtype=np.int64)
print(arr.shape)   # (5,)
print(arr.dtype)   # int64
print(arr.strides) # (8,) - 字节步长

输出示例：

>>> arr.shape
(5,)
>>> arr.dtype
int64
>>> arr.strides
(8,)

数据类型

功能说明：声明NumPy数组类型时指定dtype和维度。

cimport numpy as cnp

# 声明NumPy数组类型
cdef cnp.ndarray[cnp.int64_t, ndim=1] int_arr
cdef cnp.ndarray[cnp.float64_t, ndim=2] float_arr

# dtype代码对照表
# b1 - boolean
# i1, i2, i4, i8 - signed integers (1/2/4/8字节)
# u1, u2, u4, u8 - unsigned integers
# f4, f8 - floats (32/64位)
# c8, c16 - complex

内存布局

功能说明：C连续（行优先）和Fortran连续（列优先）的内存布局差异。

# C连续（行优先，Python默认）
c_arr = np.array([[1, 2], [3, 4]], order='C')
# 内存布局: [1, 2, 3, 4] - 行内连续

# Fortran连续（列优先）
f_arr = np.array([[1, 2], [3, 4]], order='F')
# 内存布局: [1, 3, 2, 4] - 列内连续

输出示例：

>>> c_arr.strides  # 行优先
(16, 8)
>>> f_arr.strides  # 列优先
(8, 16)

---

8.2 typed memoryview

memoryview语法

功能说明：声明memoryview类型的数组参数。

# 声明memoryview
cdef double[:] vec1d    # 一维double数组
cdef int[:, :] mat2d    # 二维int矩阵
cdef double[:, :, :] vol3d  # 三维double体积数据

# 函数参数 - nogil安全的数组操作
cpdef double mean(double[:] arr) nogil:
    cdef int n = arr.shape[0]
    cdef int i
    cdef double total = 0.0
    for i in range(n):
        total += arr[i]
    return total / n

输出示例：

>>> arr = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> mean(arr)
3.0

切片操作

功能说明：memoryview切片创建新视图，修改影响原数组。

import numpy as np

cdef double[:] arr = np.array([1.0, 2.0, 3.0, 4.0, 5.0])

# 获取切片（创建新view）
cdef double[:] sub = arr[1:4]  # [2.0, 3.0, 4.0]

# 修改会影响原数组
sub[0] = 100.0
# arr[1] 变为 100.0

# 多维切片
cdef double[:, :] matrix = np.array([[1, 2], [3, 4], [5, 6]])
cdef double[:, :] submat = matrix[0:2, 0:2]  # 前两行两列

输出示例：

>>> arr[1]
100.0
>>> submat
array([[1, 2],
       [3, 4]])

连续性与副本

功能说明：检查数组连续性，创建连续副本。

cdef double[:] arr = np.array([1, 2, 3, 4, 5])

# 检查连续性
print(arr.c_contiguous)  # True
print(arr.f_contiguous)  # 可能False

# 强制C连续副本
cdef double[::1] c_arr = arr.copy()

# 强制F连续
cdef double[::1, ::1] f_arr = np.asfortranarray(arr.reshape((5, 1)))

---

8.3 数组操作优化

批量元素操作

功能说明：Cython循环替代NumPy操作，获得更高性能。

# Python循环（慢）
def add_py(a, b):
    result = np.empty_like(a)
    for i in range(len(a)):
        result[i] = a[i] + b[i]
    return result

# Cython版本（快）
import numpy as np

cpdef np.ndarray[np.float64_t, ndim=1] add_cy(
    np.ndarray[np.float64_t, ndim=1] a,
    np.ndarray[np.float64_t, ndim=1] b
):
    cdef np.ndarray[np.float64_t, ndim=1] result = np.empty_like(a)
    cdef np.float64_t[:] a_view = a
    cdef np.float64_t[:] b_view = b
    cdef np.float64_t[:] result_view = result
    cdef int i
    cdef int n = a.shape[0]
    for i in range(n):
        result_view[i] = a_view[i] + b_view[i]
    return result

性能对比（1000万元素）：

方式	耗时
Python循环	~2000ms
NumPy原生	~50ms
Cython typed view	~30ms

广播机制

功能说明：标量操作自动广播到所有元素。

# 标量操作 - Cython自动广播
cdef void add_scalar(double[:] arr, double scalar) nogil:
    cdef int i
    for i in range(arr.shape[0]):
        arr[i] += scalar

# 行广播
cdef void add_to_rows(double[:, :] matrix, double[:] row_vec) nogil:
    cdef int i, j
    for i in range(matrix.shape[0]):
        for j in range(matrix.shape[1]):
            matrix[i, j] += row_vec[j]

输出示例：

>>> data = np.array([1.0, 2.0, 3.0])
>>> add_scalar(data, 10)
>>> data
array([11., 12., 13.])

原地操作

功能说明：原地修改数组节省内存。

# 原地操作节省内存
cdef void square_inplace(double[:] arr) nogil:
    cdef int i
    for i in range(arr.shape[0]):
        arr[i] = arr[i] * arr[i]

# 使用
cdef double[:] data = np.array([1.0, 2.0, 3.0, 4.0])
square_inplace(data)
# data 变为 [1.0, 4.0, 9.0, 16.0]

输出示例：

>>> data
array([ 1.,  4.,  9., 16.])

---

8.4 高级索引

Fancy indexing

功能说明：使用整数数组索引，返回副本。

import numpy as np

cdef double[:] arr = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
cdef int[:] indices = np.array([0, 2, 4])

# Fancy索引（返回副本）
cdef double[:] selected = arr[indices]
print(selected[0])  # 1.0

输出示例：

>>> selected
array([1., 3., 5.])

布尔掩码

功能说明：布尔数组过滤元素。

cdef double[:] arr = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
cdef bint[:] mask = np.array([True, False, True, False, True])

# 掩码过滤
cdef int count = 0
cdef int i
for i in range(arr.shape[0]):
    if mask[i]:
        count += 1

# 生成掩码
cdef bint[:] create_mask(double[:] arr, double threshold) nogil:
    cdef int n = arr.shape[0]
    cdef bint[:] mask = np.zeros(n, dtype=np.bool_)
    cdef int i
    for i in range(n):
        mask[i] = arr[i] > threshold
    return mask

输出示例：

>>> create_mask(arr, 2.5)
array([False, False,  True,  True,  True])

多维索引优化

功能说明：矩阵的行列索引和切片。

cdef double[:, :] matrix = np.array([
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
])

# 行列索引
cdef int row = 1
cdef int col = 2
cdef double val = matrix[row, col]  # 6.0

# 行切片
cdef double[:] first_row = matrix[0, :]  # [1, 2, 3]
cdef double[:] last_col = matrix[:, 2]   # [3, 6, 9]

# 遍历矩阵
cdef int i, j
cdef double total = 0.0
for i in range(matrix.shape[0]):
    for j in range(matrix.shape[1]):
        total += matrix[i, j]

输出示例：

>>> val
6.0
>>> first_row
array([1., 2., 3.])
>>> total
45.0

---

8.5 NumPy与C互操作

直接缓冲区接口

功能说明：任何支持buffer协议的对象都可用memoryview处理。

from cython cimport view

# buffer protocol支持
cdef void process_buffer(object arr):
    cdef double[:] view = arr
    # 可以直接处理任何支持buffer协议的对象
    cdef int i
    for i in range(view.shape[0]):
        view[i] *= 2.0

PyArray_DATA使用

功能说明：获取NumPy数组的原始数据指针。

cimport numpy as cnp

# 获取数据指针
cdef np.ndarray[cnp.float64_t, ndim=1] arr = np.array([1.0, 2.0, 3.0])
cdef double* data_ptr = <double*>cnp.PyArray_DATA(arr)

# 通过指针访问
data_ptr[0] = 100.0
# arr 变为 [100.0, 2.0, 3.0]

输出示例：

>>> arr
array([100.,   2.,   3.])

零拷贝技术

功能说明：memoryview不复制数据，只创建新视图。

# memoryview实现零拷贝
cdef double[:] fast_view(np.ndarray[cnp.float64_t, ndim=1] arr):
    # 不复制数据，创建新视图
    return arr

# 使用
cdef np.ndarray[cnp.float64_t, ndim=1] np_arr = np.array([1.0, 2.0, 3.0])
cdef double[:] mv = fast_view(np_arr)
mv[0] = 999.0
# np_arr[0] 也变为 999.0（零拷贝）

输出示例：

>>> np_arr
array([999.,   2.,   3.])
>>> mv[0] = 0
>>> np_arr[0]
0.0

---

本章小结

数组技术对比

技术	性能	适用场景
ndarray + typed view	极高	大数组运算
memoryview	零拷贝	数组切片
buffer protocol	高	通用数组接口
C指针	最快	危险，需小心

最佳实践清单

1. 函数参数用typed memoryview（double[:]）
2. 批量操作用nogil循环
3. 原地修改节省内存
4. 检查c_contiguous确保最优访问模式
5. 使用PyArray_DATA直接操作指针（需谨慎）

---

练习

1. 实现矩阵乘法（使用typed memoryview）
2. 比较NumPy操作和Cython循环的性能差异
3. 实现带掩码的数组过滤函数
4. 创建零拷贝的memoryview切片
5. 实现批量广播操作（标量+数组）
6. 测试C连续和Fortran连续的访问性能差异