Jacky's blog
Jacky
2023-09-30
目录

numpy

# 基础知识

NumPy是Python科学计算的基础包，提供了高性能的多维数组对象和用于处理这些数组的工具。

# 数组创建与操作

import numpy as np

# 创建数组
arr1 = np.array([1, 2, 3, 4, 5])  # 从列表创建
arr2 = np.zeros((3, 3))           # 创建3x3的零矩阵
arr3 = np.ones((2, 4))            # 创建2x4的全1矩阵
arr4 = np.empty((2, 3))           # 创建未初始化的数组
arr5 = np.arange(0, 10, 2)        # 创建等差数列 [0,2,4,6,8]
arr6 = np.linspace(0, 1, 5)       # 创建5个等间距的数 [0, 0.25, 0.5, 0.75, 1]

# 数组形状操作
arr = np.array([[1, 2, 3], [4, 5, 6]])
reshaped = arr.reshape(3, 2)      # 改变形状为3x2
flattened = arr.flatten()         # 展平为一维数组

# 数组切片与索引
arr = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(arr[0, 1])      # 获取第一行第二列的元素
print(arr[1:3, :])    # 获取第2-3行的所有列
print(arr[:, 1])      # 获取所有行的第二列

# 数组拼接与分割
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
c = np.concatenate([a, b])        # 拼接数组。 连接操作，将数组在指定轴上连接起来
                                  # + 是加法操作，对对应位置的元素进行相加
arr = np.array([[1, 2], [3, 4], [5, 6]])
split_arr = np.split(arr, 3)      # 将数组分割成3份
split_arr_col = np.split(arr, 2, axis=1)  # 沿着列分割， 第一列: array([[1], [3], [5]]); 第二列: array([[2], [4], [6]])

# 网格生成
x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
X, Y = np.meshgrid(x, y)         # 生成二维网格
print("X网格:")
print(X)  # [[1 2 3]
          #  [1 2 3]
          #  [1 2 3]]
print("Y网格:")
print(Y)  # [[4 4 4]
          #  [5 5 5]
          #  [6 6 6]]

# meshgrid 的高级用法
x = np.linspace(-2, 2, 5)
y = np.linspace(-2, 2, 5)
X, Y = np.meshgrid(x, y, indexing='ij')  # 使用 'ij' 索引方式
Z = X**2 + Y**2  # 计算每个网格点的值

# 三维网格
x = np.array([1, 2])
y = np.array([3, 4])
z = np.array([5, 6])
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
print("三维网格形状:", X.shape, Y.shape, Z.shape)  # (2, 2, 2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

# 数学运算

NumPy提供了丰富的数学运算功能，支持向量化操作，大大提高了计算效率。

# 数学运算模块概览

  • np.linalg: 线性代数相关函数。 示例 (opens new window)
  • np.random: 随机数生成
  • np.fft: 快速傅里叶变换
  • np.polynomial: 多项式运算
  • np.ma: 掩码数组操作

# 基本运算

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

# 基本运算
print(a + b)          # 加法 [5, 7, 9]
print(a - b)          # 减法 [-3, -3, -3]
print(a * b)          # 乘法 [4, 10, 18]
print(a / b)          # 除法 [0.25, 0.4, 0.5]
print(a ** 2)         # 幂运算 [1, 4, 9]

# 矩阵运算
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
print(np.dot(A, B))   # 矩阵乘法
print(np.matmul(A, B))# 矩阵乘法（另一种方式）

# 统计函数
arr = np.array([1, 2, 3, 4, 5])
print(np.mean(arr))   # 平均值
print(np.std(arr))    # 标准差
print(np.var(arr))    # 方差

# 三角函数
angles = np.array([0, 30, 45, 60, 90])
radians = np.deg2rad(angles)
print(np.sin(radians))# 正弦值
print(np.cos(radians))# 余弦值
print(np.tan(radians))# 正切值
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

# 高级数学函数

# 1. 指数和对数函数

arr = np.array([1, 2, 3, 4, 5])
# arr.tolist()

# 指数函数
print(np.exp(arr))    # e^x
print(np.exp2(arr))   # 2^x
print(np.power(arr, 2)) # x^2

# 对数函数
print(np.log(arr))    # ln(x)
print(np.log2(arr))   # log2(x)
print(np.log10(arr))  # log10(x)

# 双曲函数
print(np.sinh(arr))   # 双曲正弦
print(np.cosh(arr))   # 双曲余弦
print(np.tanh(arr))   # 双曲正切
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# 2. 舍入和取整函数

arr = np.array([1.2, 2.7, -1.8, 3.4])

print(np.floor(arr))  # 向下取整
print(np.ceil(arr))   # 向上取整
print(np.round(arr))  # 四舍五入
print(np.trunc(arr))  # 截断小数部分
print(np.fix(arr))    # 向零取整
1
2
3
4
5
6
7

# 3. 复数运算

# 创建复数数组
z = np.array([1+2j, 3+4j, 5+6j])

print(np.real(z))     # 实部
print(np.imag(z))     # 虚部
print(np.conj(z))     # 共轭复数
print(np.angle(z))    # 相位角
print(np.abs(z))      # 模长
1
2
3
4
5
6
7
8

# 4. 累积和差分

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

# 累积函数
print(np.cumsum(arr))     # 累积和 [1, 3, 6, 10, 15]
print(np.cumprod(arr))    # 累积积 [1, 2, 6, 24, 120]

# 差分
print(np.diff(arr))       # 一阶差分 [1, 1, 1, 1]
print(np.diff(arr, n=2))  # 二阶差分 [0, 0, 0]
1
2
3
4
5
6
7
8
9

# 5. 梯度计算

# 计算梯度
x = np.array([1, 2, 4, 7, 11, 16])
gradient = np.gradient(x)
print(gradient)  # 计算相邻元素的差分

# 二维梯度
X, Y = np.meshgrid(np.linspace(0, 1, 5), np.linspace(0, 1, 5))
Z = X**2 + Y**2
grad_x, grad_y = np.gradient(Z)
1
2
3
4
5
6
7
8
9

# 统计函数详解

# 1. 描述性统计

arr = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

# 基本统计量
print(np.mean(arr))       # 均值
print(np.median(arr))     # 中位数
print(np.std(arr))        # 标准差； 衡量数据分散程度的统计量，表示数据点与平均值的平均距离。
print(np.var(arr))        # 方差
print(np.min(arr))        # 最小值
print(np.max(arr))        # 最大值
print(np.ptp(arr))        # 极差 (max - min)

# 分位数
print(np.percentile(arr, [25, 50, 75]))  # 四分位数
print(np.quantile(arr, [0.1, 0.5, 0.9])) # 分位数
1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 2. 相关性分析

# 协方差和相关系数
x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 4, 5, 4, 5])

print(np.cov(x, y))       # 协方差矩阵
print(np.corrcoef(x, y))  # 相关系数矩阵

# 自相关
def autocorr(x, lag=1):
    return np.corrcoef(x[:-lag], x[lag:])[0, 1]

print(autocorr(x, lag=1))
1
2
3
4
5
6
7
8
9
10
11
12

# 3. 分布函数

# 正态分布
from scipy import stats

# 生成正态分布数据
normal_data = np.random.normal(0, 1, 1000)

# 计算分布参数
mean, std = stats.norm.fit(normal_data)
print(f"拟合的均值: {mean:.3f}, 标准差: {std:.3f}")

# 其他分布
uniform_data = np.random.uniform(0, 1, 1000)
exponential_data = np.random.exponential(1, 1000)
1
2
3
4
5
6
7
8
9
10
11
12
13

# 信号处理函数

# 1. 卷积和滤波

# 一维卷积
signal = np.array([1, 2, 3, 4, 5])
kernel = np.array([0.5, 1, 0.5])
convolved = np.convolve(signal, kernel, mode='same')

# 二维卷积（图像滤波）
image = np.random.rand(10, 10)
filter_kernel = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]]) / 9
filtered_image = np.convolve2d(image, filter_kernel, mode='same')
1
2
3
4
5
6
7
8
9

# 2. 傅里叶变换

# 一维FFT
signal = np.sin(2 * np.pi * 10 * np.linspace(0, 1, 1000))
fft_result = np.fft.fft(signal)
freq = np.fft.fftfreq(len(signal))

# 功率谱密度
psd = np.abs(fft_result)**2

# 逆FFT
reconstructed = np.fft.ifft(fft_result)
1
2
3
4
5
6
7
8
9
10

# 3. 窗函数

# 常用窗函数
N = 100
hann_window = np.hanning(N)
hamming_window = np.hamming(N)
blackman_window = np.blackman(N)

# 应用窗函数
signal = np.sin(2 * np.pi * 10 * np.linspace(0, 1, N))
windowed_signal = signal * hann_window
1
2
3
4
5
6
7
8
9

# 优化和求解函数

# 1. 线性规划

from scipy.optimize import linprog

# 线性规划示例
# 最小化: -x1 - 2x2
# 约束: x1 + x2 <= 3, x1 >= 0, x2 >= 0
c = [-1, -2]  # 目标函数系数
A = [[1, 1]]  # 约束矩阵
b = [3]       # 约束右端
bounds = [(0, None), (0, None)]  # 变量边界

result = linprog(c, A_ub=A, b_ub=b, bounds=bounds)
print(f"最优解: {result.x}")
print(f"最优值: {result.fun}")
1
2
3
4
5
6
7
8
9
10
11
12
13

# 2. 非线性优化

from scipy.optimize import minimize

# 最小化函数 f(x) = x^2 + 2x + 1
def objective(x):
    return x[0]**2 + 2*x[0] + 1

# 初始猜测
x0 = [0]

# 优化
result = minimize(objective, x0, method='BFGS')
print(f"最优解: {result.x[0]}")
print(f"最优值: {result.fun}")
1
2
3
4
5
6
7
8
9
10
11
12
13

# 3. 方程求解

from scipy.optimize import fsolve

# 求解非线性方程组
def equations(vars):
    x, y = vars
    eq1 = x**2 + y**2 - 1
    eq2 = x - y
    return [eq1, eq2]

# 初始猜测
initial_guess = [0.5, 0.5]

# 求解
solution = fsolve(equations, initial_guess)
print(f"解: x={solution[0]:.3f}, y={solution[1]:.3f}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 数值积分和微分

# 1. 数值积分

from scipy import integrate

# 定积分 ∫(0 to 1) x^2 dx
def integrand(x):
    return x**2

result, error = integrate.quad(integrand, 0, 1)
print(f"积分结果: {result:.3f}, 误差: {error:.2e}")

# 二重积分
def integrand_2d(x, y):
    return x**2 + y**2

result_2d, error_2d = integrate.dblquad(integrand_2d, 0, 1, lambda x: 0, lambda x: 1)
print(f"二重积分结果: {result_2d:.3f}")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 2. 数值微分

# 数值微分
def f(x):
    return x**2

# 中心差分
def numerical_derivative(f, x, h=1e-6):
    return (f(x + h) - f(x - h)) / (2 * h)

x = 2.0
exact_derivative = 2 * x  # f'(x) = 2x
numerical_result = numerical_derivative(f, x)
print(f"精确导数: {exact_derivative}")
print(f"数值导数: {numerical_result:.6f}")
1
2
3
4
5
6
7
8
9
10
11
12
13

# 插值和拟合

# 1. 多项式插值

from scipy.interpolate import interp1d, lagrange

# 数据点
x = np.array([0, 1, 2, 3, 4])
y = np.array([0, 1, 4, 9, 16])

# 线性插值
linear_interp = interp1d(x, y, kind='linear')
x_new = np.linspace(0, 4, 100)
y_linear = linear_interp(x_new)

# 拉格朗日插值
poly = lagrange(x, y)
y_lagrange = poly(x_new)
1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 2. 样条插值

from scipy.interpolate import CubicSpline

# 三次样条插值
cs = CubicSpline(x, y)
y_spline = cs(x_new)

# 样条导数
derivative = cs.derivative()
y_derivative = derivative(x_new)
1
2
3
4
5
6
7
8
9

# 3. 曲线拟合

from scipy.optimize import curve_fit

# 定义拟合函数
def func(x, a, b, c):
    return a * np.exp(-b * x) + c

# 生成带噪声的数据
x_data = np.linspace(0, 4, 50)
y_data = func(x_data, 2.5, 1.3, 0.5) + 0.2 * np.random.normal(size=len(x_data))

# 拟合
popt, pcov = curve_fit(func, x_data, y_data)
print(f"拟合参数: a={popt[0]:.3f}, b={popt[1]:.3f}, c={popt[2]:.3f}")
1
2
3
4
5
6
7
8
9
10
11
12
13

# 特殊函数

# 1. 贝塞尔函数

from scipy.special import j0, j1, y0, y1

x = np.linspace(0, 10, 100)
bessel_j0 = j0(x)  # 第一类贝塞尔函数
bessel_j1 = j1(x)  # 第一类贝塞尔函数
bessel_y0 = y0(x)  # 第二类贝塞尔函数
bessel_y1 = y1(x)  # 第二类贝塞尔函数
1
2
3
4
5
6
7

# 2. 伽马函数

from scipy.special import gamma, loggamma

x = np.linspace(1, 5, 100)
gamma_values = gamma(x)
log_gamma_values = loggamma(x)
1
2
3
4
5

# 3. 误差函数

from scipy.special import erf, erfc

x = np.linspace(-3, 3, 100)
error_function = erf(x)
complementary_error = erfc(x)
1
2
3
4
5

# 数组操作

# 广播机制(Broadcasting)

广播是NumPy中一种强大的机制，它允许不同形状的数组进行运算。广播遵循以下规则：

  1. 广播规则：

    • 从右向左比较数组的维度
    • 如果维度相等或其中一个为1，则可以广播
    • 如果维度不相等且都不为1，则报错
    • 缺失的维度会被视为1

  2. 广播步骤：

    • 确定最终输出数组的形状
    • 将较小的数组扩展到与较大数组相同的形状
    • 进行元素级运算

  3. 常见广播示例：

# 示例1：向量与矩阵相加
a = np.array([[1, 2, 3], [4, 5, 6]])  # 形状: (2, 3)
b = np.array([1, 2, 3])               # 形状: (3,)
print(a + b)  # b被广播为 [[1,2,3], [1,2,3]]
# 输出: [[2, 4, 6], [5, 7, 9]]

# 示例2：不同维度的数组运算
a = np.array([1, 2, 3])      # 形状: (3,)
b = np.array([[1], [2]])     # 形状: (2, 1)
print(a + b)  # a被广播为 [[1,2,3], [1,2,3]], b被广播为 [[1,1,1], [2,2,2]]
# 输出: [[2, 3, 4], [3, 4, 5]]

# 示例3：三维数组的广播
a = np.array([[[1], [2]], [[3], [4]]])  # 形状: (2, 2, 1)
b = np.array([1, 2, 3])                 # 形状: (3,)
print(a + b)  # b被广播为 [[[1,2,3], [1,2,3]], [[1,2,3], [1,2,3]]]
# 输出: [[[2,3,4], [3,4,5]], [[4,5,6], [5,6,7]]]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
  1. 广播的应用场景：
# 标准化数据
data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
mean = np.mean(data, axis=0)  # 计算每列的均值
std = np.std(data, axis=0)    # 计算每列的标准差
normalized = (data - mean) / std  # 广播用于标准化

# 图像处理
image = np.array([[1, 2], [3, 4]])  # 2x2图像
filter = np.array([1, -1])          # 1x2滤波器
# 滤波器被广播到每一行
filtered = image * filter

# 机器学习中的批量运算
weights = np.array([0.1, 0.2, 0.3])  # 权重向量
batch = np.array([[1, 2, 3], [4, 5, 6]])  # 批量数据
# 权重被广播到每个样本
predictions = np.dot(batch, weights)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

  1. 广播的性能考虑：

    • 广播不会实际复制数据，而是创建视图
    • 对于大型数组，广播操作可能比显式复制更高效
    • 但要注意内存使用，避免创建过大的临时数组

  2. 常见错误和解决方案：

# 错误示例：维度不兼容
a = np.array([1, 2, 3])
b = np.array([1, 2])
try:
    print(a + b)  # 会报错：operands could not be broadcast together
except ValueError as e:
    print("错误：", e)

# 解决方案1：调整数组形状
b = np.array([1, 2, 3])  # 确保维度匹配

# 解决方案2：使用reshape
b = np.array([1, 2]).reshape(2, 1)  # 改变形状以匹配

# 解决方案3：使用expand_dims
b = np.array([1, 2])
b = np.expand_dims(b, axis=0)  # 添加新维度
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
  1. 最佳实践：
    • 尽量使用广播而不是显式循环
    • 注意数组的形状，确保广播规则满足
    • 对于复杂的广播操作，使用 np.broadcast_to() 显式控制广播
    • 在性能关键代码中，考虑使用 np.einsum() 进行复杂的广播运算

# 轴(axis)概念

arr = np.array([[1, 2, 3], [4, 5, 6]])
print(np.sum(arr, axis=0))  # 沿着0轴（行）求和。 沿着行方向（垂直方向）求和。 结果 [5 7 9]
print(np.sum(arr, axis=1))  # 沿着1轴（列）求和。沿着列方向（水平方向）求和 。 结果 [6, 16]
print(np.sum(arr))  # 不指定axis（所有元素求和）。结果: 21
1
2
3
4

# 条件筛选

arr = np.array([1, 2, 3, 4, 5])
mask = arr > 3
print(arr[mask])  # 输出大于3的元素
1
2
3

# 排序和去重

arr = np.array([3, 1, 4, 1, 5, 9, 2, 6])
print(np.sort(arr))           # 排序
print(np.unique(arr))         # 去重
1
2
3

# 线性代数

# 矩阵转置
A = np.array([[1, 2], [3, 4]])
print(A.T)  # 转置

# 矩阵求逆
A = np.array([[1, 2], [3, 4]])
print(np.linalg.inv(A))  # 求逆矩阵

# 特征值与特征向量
A = np.array([[1, 2], [2, 1]])
eigenvalues, eigenvectors = np.linalg.eig(A)

# 解线性方程组
A = np.array([[1, 2], [3, 4]])
b = np.array([5, 6])
x = np.linalg.solve(A, b)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

# 文件操作

# 保存数据
arr = np.array([1, 2, 3, 4, 5])
np.save('array.npy', arr)           # 保存单个数组
np.savez('arrays.npz', a=arr, b=arr) # 保存多个数组

# 加载数据
loaded_arr = np.load('array.npy')
data = np.load('arrays.npz')
a = data['a']
b = data['b']
1
2
3
4
5
6
7
8
9
10

# 随机数生成

# 均匀分布
uniform = np.random.rand(3, 3)  # 生成3x3的均匀分布随机数

# 正态分布
normal = np.random.randn(3, 3)  # 生成3x3的正态分布随机数

# 随机整数
integers = np.random.randint(0, 10, size=(3, 3))  # 生成3x3的随机整数矩阵
1
2
3
4
5
6
7
8

# 高级特性

# 通用函数(ufunc)

# 自定义ufunc
def my_func(x):
    return x * 2

my_ufunc = np.frompyfunc(my_func, 1, 1)
arr = np.array([1, 2, 3])
result = my_ufunc(arr)
1
2
3
4
5
6
7

# 结构化数组

# 创建结构化数组
dtype = [('name', 'U10'), ('age', 'i4'), ('height', 'f4')]
values = [('Alice', 25, 1.75), ('Bob', 30, 1.85)]
arr = np.array(values, dtype=dtype)
1
2
3
4

# 内存视图与拷贝

# 视图（共享内存）
arr = np.array([1, 2, 3])
view = arr.view()
view[0] = 10  # 会改变原始数组

# 拷贝（独立内存）
arr = np.array([1, 2, 3])
copy = arr.copy()
copy[0] = 10  # 不会改变原始数组
1
2
3
4
5
6
7
8
9

# 掩码数组

# 创建掩码数组
arr = np.array([1, 2, 3, 4, 5])
mask = np.array([True, False, True, False, True])
masked_arr = np.ma.masked_array(arr, mask=~mask)
1
2
3
4

# 网格生成 (meshgrid)

meshgrid 是 NumPy 中用于生成坐标网格的重要函数，广泛应用于科学计算、数据可视化和数值分析。

# 基本概念

import numpy as np

# 基本用法
x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
X, Y = np.meshgrid(x, y)

print("X网格:")
print(X)
# 输出:
# [[1 2 3]
#  [1 2 3]
#  [1 2 3]]

print("Y网格:")
print(Y)
# 输出:
# [[4 4 4]
#  [5 5 5]
#  [6 6 6]]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# 索引方式

# 默认索引方式 (xy)
x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
X_xy, Y_xy = np.meshgrid(x, y, indexing='xy')
print("xy索引方式:")
print("X形状:", X_xy.shape)  # (3, 3)
print("Y形状:", Y_xy.shape)  # (3, 3)

# ij索引方式
X_ij, Y_ij = np.meshgrid(x, y, indexing='ij')
print("ij索引方式:")
print("X形状:", X_ij.shape)  # (3, 3)
print("Y形状:", Y_ij.shape)  # (3, 3)

# 稀疏网格
X_sparse, Y_sparse = np.meshgrid(x, y, sparse=True)
print("稀疏网格形状:", X_sparse.shape, Y_sparse.shape)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# 实际应用示例

# 1. 函数可视化
def visualize_function():
    # 创建网格
    x = np.linspace(-3, 3, 50)
    y = np.linspace(-3, 3, 50)
    X, Y = np.meshgrid(x, y)
    
    # 计算函数值
    Z = np.sin(X) * np.cos(Y)
    
    # 可视化
    import matplotlib.pyplot as plt
    
    plt.figure(figsize=(15, 5))
    
    # 等高线图
    plt.subplot(1, 3, 1)
    plt.contour(X, Y, Z, levels=20)
    plt.colorbar()
    plt.title('等高线图')
    plt.xlabel('X')
    plt.ylabel('Y')
    
    # 3D表面图
    from mpl_toolkits.mplot3d import Axes3D
    ax = plt.subplot(1, 3, 2, projection='3d')
    surf = ax.plot_surface(X, Y, Z, cmap='viridis')
    plt.colorbar(surf)
    plt.title('3D表面图')
    
    # 热力图
    plt.subplot(1, 3, 3)
    plt.imshow(Z, extent=[-3, 3, -3, 3], origin='lower', cmap='hot')
    plt.colorbar()
    plt.title('热力图')
    
    plt.tight_layout()
    plt.show()

# 2. 向量场计算
def vector_field_example():
    # 创建网格
    x = np.linspace(-2, 2, 20)
    y = np.linspace(-2, 2, 20)
    X, Y = np.meshgrid(x, y)
    
    # 计算向量场
    U = -Y  # x方向分量
    V = X   # y方向分量
    
    # 计算场强
    magnitude = np.sqrt(U**2 + V**2)
    
    # 可视化
    import matplotlib.pyplot as plt
    
    plt.figure(figsize=(12, 5))
    
    # 向量场图
    plt.subplot(1, 2, 1)
    plt.quiver(X, Y, U, V, magnitude, cmap='viridis', scale=50)
    plt.colorbar()
    plt.title('向量场')
    plt.xlabel('X')
    plt.ylabel('Y')
    
    # 流线图
    plt.subplot(1, 2, 2)
    plt.streamplot(X, Y, U, V, density=1.5, color=magnitude, cmap='viridis')
    plt.colorbar()
    plt.title('流线图')
    plt.xlabel('X')
    plt.ylabel('Y')
    
    plt.tight_layout()
    plt.show()

# 3. 数值积分
def numerical_integration():
    # 定义积分区域
    x = np.linspace(0, 1, 100)
    y = np.linspace(0, 1, 100)
    X, Y = np.meshgrid(x, y)
    
    # 定义被积函数
    def f(x, y):
        return x**2 + y**2
    
    # 计算函数值
    Z = f(X, Y)
    
    # 数值积分（矩形法则）
    dx = x[1] - x[0]
    dy = y[1] - y[0]
    integral = np.sum(Z) * dx * dy
    
    print(f"数值积分结果: {integral:.6f}")
    print(f"理论值: {2/3:.6f}")
    
    return integral

# 4. 三维网格
def three_dimensional_grid():
    # 创建三维网格
    x = np.linspace(-1, 1, 10)
    y = np.linspace(-1, 1, 10)
    z = np.linspace(-1, 1, 10)
    X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
    
    # 计算三维函数
    R = np.sqrt(X**2 + Y**2 + Z**2)
    V = np.exp(-R)
    
    print("三维网格形状:", X.shape, Y.shape, Z.shape)
    print("函数值范围:", np.min(V), "到", np.max(V))
    
    return X, Y, Z, V

# 5. 有限差分法
def finite_difference_method():
    # 网格参数
    nx, ny = 30, 30
    dx = dy = 0.1
    x = np.linspace(0, (nx-1)*dx, nx)
    y = np.linspace(0, (ny-1)*dy, ny)
    X, Y = np.meshgrid(x, y)
    
    # 初始化解
    u = np.zeros((nx, ny))
    
    # 边界条件
    u[0, :] = 0      # 左边界
    u[-1, :] = 0     # 右边界
    u[:, 0] = 0      # 下边界
    u[:, -1] = 100   # 上边界
    
    # 迭代求解拉普拉斯方程
    for iteration in range(500):
        u_old = u.copy()
        for i in range(1, nx-1):
            for j in range(1, ny-1):
                u[i, j] = 0.25 * (u_old[i+1, j] + u_old[i-1, j] + 
                                 u_old[i, j+1] + u_old[i, j-1])
        
        # 检查收敛
        if np.max(np.abs(u - u_old)) < 1e-6:
            print(f"收敛于第 {iteration} 次迭代")
            break
    
    return X, Y, u
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

# 性能优化技巧

# 1. 使用稀疏网格减少内存使用
def sparse_grid_example():
    x = np.linspace(-5, 5, 1000)
    y = np.linspace(-5, 5, 1000)
    
    # 密集网格（内存密集）
    X_dense, Y_dense = np.meshgrid(x, y)
    print("密集网格内存使用:", X_dense.nbytes / 1024 / 1024, "MB")
    
    # 稀疏网格（内存节省）
    X_sparse, Y_sparse = np.meshgrid(x, y, sparse=True)
    print("稀疏网格内存使用:", X_sparse.nbytes / 1024 / 1024, "MB")
    
    return X_sparse, Y_sparse

# 2. 分块处理大网格
def chunked_grid_processing():
    # 大网格分块处理
    nx, ny = 1000, 1000
    chunk_size = 100
    
    result = np.zeros((nx, ny))
    
    for i in range(0, nx, chunk_size):
        for j in range(0, ny, chunk_size):
            # 处理当前块
            x_chunk = np.linspace(i*0.01, (i+chunk_size)*0.01, chunk_size)
            y_chunk = np.linspace(j*0.01, (j+chunk_size)*0.01, chunk_size)
            X_chunk, Y_chunk = np.meshgrid(x_chunk, y_chunk)
            
            # 计算块内函数值
            Z_chunk = np.sin(X_chunk) * np.cos(Y_chunk)
            
            # 存储结果
            result[i:i+chunk_size, j:j+chunk_size] = Z_chunk
    
    return result
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

# 性能优化建议

  1. 使用向量化操作代替循环
  2. 避免频繁的数组复制
  3. 使用适当的数据类型
  4. 利用广播机制
  5. 使用NumPy的内置函数

# 常见问题解决

  1. 内存溢出：使用适当的数据类型和分块处理
  2. 性能问题：使用向量化操作和NumPy内置函数
  3. 精度问题：注意浮点数计算的精度限制

# 最佳实践

  1. 始终检查数组的形状和数据类型
  2. 使用适当的错误处理机制
  3. 注意内存使用
  4. 保持代码的可读性和可维护性

# 实际应用案例

# 1. 图像处理应用

# 图像滤波和增强

import numpy as np
from PIL import Image

# 加载图像
image = np.array(Image.open('image.jpg'))

# 灰度化
def rgb_to_gray(rgb_image):
    return np.dot(rgb_image[..., :3], [0.299, 0.587, 0.114])

# 高斯滤波
def gaussian_filter(image, sigma=1.0):
    from scipy.ndimage import gaussian_filter
    return gaussian_filter(image, sigma=sigma)

# 边缘检测
def sobel_edge_detection(image):
    sobel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    sobel_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
    
    grad_x = np.convolve2d(image, sobel_x, mode='same')
    grad_y = np.convolve2d(image, sobel_y, mode='same')
    
    return np.sqrt(grad_x**2 + grad_y**2)

# 应用示例
gray_image = rgb_to_gray(image)
filtered_image = gaussian_filter(gray_image, sigma=2.0)
edge_image = sobel_edge_detection(gray_image)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

# 图像变换

# 图像旋转
def rotate_image(image, angle):
    from scipy.ndimage import rotate
    return rotate(image, angle, reshape=False)

# 图像缩放
def resize_image(image, scale_factor):
    from scipy.ndimage import zoom
    return zoom(image, scale_factor)

# 图像平移
def translate_image(image, dx, dy):
    from scipy.ndimage import shift
    return shift(image, (dy, dx))
1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 2. 数据分析应用

# 时间序列分析

# 移动平均
def moving_average(data, window_size):
    return np.convolve(data, np.ones(window_size)/window_size, mode='valid')

# 指数平滑
def exponential_smoothing(data, alpha=0.3):
    smoothed = np.zeros_like(data)
    smoothed[0] = data[0]
    for i in range(1, len(data)):
        smoothed[i] = alpha * data[i] + (1 - alpha) * smoothed[i-1]
    return smoothed

# 季节性分解
def seasonal_decomposition(data, period):
    # 简单季节性分解
    trend = moving_average(data, period)
    seasonal = data - trend
    residual = data - trend - seasonal
    return trend, seasonal, residual

# 应用示例
time_series = np.random.randn(100) + np.sin(np.linspace(0, 4*np.pi, 100))
ma_result = moving_average(time_series, 10)
es_result = exponential_smoothing(time_series, alpha=0.3)
trend, seasonal, residual = seasonal_decomposition(time_series, 20)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

# 数据清洗和预处理

# 异常值检测
def detect_outliers(data, threshold=3):
    z_scores = np.abs((data - np.mean(data)) / np.std(data))
    return z_scores > threshold

# 数据标准化
def standardize_data(data):
    return (data - np.mean(data)) / np.std(data)

# 数据归一化
def normalize_data(data):
    return (data - np.min(data)) / (np.max(data) - np.min(data))

# 缺失值处理
def handle_missing_values(data, method='mean'):
    if method == 'mean':
        return np.where(np.isnan(data), np.nanmean(data), data)
    elif method == 'median':
        return np.where(np.isnan(data), np.nanmedian(data), data)
    elif method == 'interpolate':
        from scipy.interpolate import interp1d
        valid_indices = ~np.isnan(data)
        if np.sum(valid_indices) > 1:
            f = interp1d(np.where(valid_indices)[0], data[valid_indices])
            return f(np.arange(len(data)))
    return data
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

# 3. 机器学习应用

# 特征工程

# 特征缩放
def feature_scaling(X, method='standard'):
    if method == 'standard':
        return (X - np.mean(X, axis=0)) / np.std(X, axis=0)
    elif method == 'minmax':
        return (X - np.min(X, axis=0)) / (np.max(X, axis=0) - np.min(X, axis=0))
    elif method == 'robust':
        from scipy.stats import iqr
        q75, q25 = np.percentile(X, [75, 25], axis=0)
        return (X - np.median(X, axis=0)) / iqr(X, axis=0)

# 特征选择
def feature_selection(X, y, method='correlation', threshold=0.5):
    if method == 'correlation':
        correlations = np.abs(np.corrcoef(X.T, y)[:-1, -1])
        return correlations > threshold
    elif method == 'variance':
        variances = np.var(X, axis=0)
        return variances > threshold

# 多项式特征
def polynomial_features(X, degree=2):
    from sklearn.preprocessing import PolynomialFeatures
    poly = PolynomialFeatures(degree=degree, include_bias=False)
    return poly.fit_transform(X)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

# 模型评估

# 交叉验证
def cross_validation(X, y, model, k=5):
    from sklearn.model_selection import KFold
    kf = KFold(n_splits=k, shuffle=True, random_state=42)
    scores = []
    
    for train_idx, test_idx in kf.split(X):
        X_train, X_test = X[train_idx], X[test_idx]
        y_train, y_test = y[train_idx], y[test_idx]
        
        model.fit(X_train, y_train)
        score = model.score(X_test, y_test)
        scores.append(score)
    
    return np.mean(scores), np.std(scores)

# 混淆矩阵
def confusion_matrix(y_true, y_pred):
    from sklearn.metrics import confusion_matrix
    return confusion_matrix(y_true, y_pred)

# ROC曲线
def roc_curve(y_true, y_pred_proba):
    from sklearn.metrics import roc_curve, auc
    fpr, tpr, _ = roc_curve(y_true, y_pred_proba)
    roc_auc = auc(fpr, tpr)
    return fpr, tpr, roc_auc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

# 4. 科学计算应用

# 数值积分和微分方程

# 常微分方程求解
from scipy.integrate import odeint

def model(y, t, k):
    dydt = -k * y
    return dydt

# 求解ODE
t = np.linspace(0, 20, 100)
y0 = 1.0
k = 0.1
solution = odeint(model, y0, t, args=(k,))

# 偏微分方程（热传导方程）
def heat_equation_solver(nx=50, nt=100, dx=0.1, dt=0.001, alpha=1.0):
    # 一维热传导方程
    x = np.linspace(0, nx*dx, nx)
    t = np.linspace(0, nt*dt, nt)
    
    # 初始条件
    u = np.zeros((nt, nx))
    u[0, :] = np.sin(np.pi * x / (nx * dx))
    
    # 边界条件
    u[:, 0] = 0
    u[:, -1] = 0
    
    # 时间步进
    for n in range(0, nt-1):
        for i in range(1, nx-1):
            u[n+1, i] = u[n, i] + alpha * dt / dx**2 * (u[n, i+1] - 2*u[n, i] + u[n, i-1])
    
    return x, t, u
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

# 网格生成和可视化应用

# 二维函数可视化
def plot_2d_function():
    # 创建网格
    x = np.linspace(-3, 3, 100)
    y = np.linspace(-3, 3, 100)
    X, Y = np.meshgrid(x, y)
    
    # 计算函数值
    Z = np.sin(np.sqrt(X**2 + Y**2))
    
    # 绘制等高线图
    import matplotlib.pyplot as plt
    plt.figure(figsize=(10, 8))
    
    # 等高线图
    plt.subplot(2, 2, 1)
    plt.contour(X, Y, Z, levels=20)
    plt.colorbar()
    plt.title('等高线图')
    
    # 填充等高线图
    plt.subplot(2, 2, 2)
    plt.contourf(X, Y, Z, levels=20)
    plt.colorbar()
    plt.title('填充等高线图')
    
    # 3D表面图
    from mpl_toolkits.mplot3d import Axes3D
    ax = plt.subplot(2, 2, 3, projection='3d')
    surf = ax.plot_surface(X, Y, Z, cmap='viridis')
    plt.colorbar(surf)
    plt.title('3D表面图')
    
    # 热力图
    plt.subplot(2, 2, 4)
    plt.imshow(Z, extent=[-3, 3, -3, 3], origin='lower', cmap='hot')
    plt.colorbar()
    plt.title('热力图')
    
    plt.tight_layout()
    plt.show()

# 三维网格应用
def plot_3d_function():
    # 创建三维网格
    x = np.linspace(-2, 2, 20)
    y = np.linspace(-2, 2, 20)
    z = np.linspace(-2, 2, 20)
    X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
    
    # 计算三维函数值
    R = np.sqrt(X**2 + Y**2 + Z**2)
    V = np.sin(R) / (R + 1e-10)  # 避免除零
    
    # 绘制等值面
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # 选择特定等值面
    level = 0.5
    ax.contour3D(X, Y, Z, V, levels=[level], cmap='viridis')
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title(f'等值面 V = {level}')
    
    plt.show()

# 向量场可视化
def plot_vector_field():
    # 创建网格
    x = np.linspace(-2, 2, 20)
    y = np.linspace(-2, 2, 20)
    X, Y = np.meshgrid(x, y)
    
    # 计算向量场
    U = -Y  # x方向分量
    V = X   # y方向分量
    
    # 绘制向量场
    import matplotlib.pyplot as plt
    
    plt.figure(figsize=(10, 8))
    plt.quiver(X, Y, U, V, scale=50)
    plt.streamplot(X, Y, U, V, density=1.5, color='red', linewidth=1)
    plt.title('向量场可视化')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.axis('equal')
    plt.show()

# 数值积分中的网格应用
def numerical_integration_2d():
    # 使用网格进行二重积分
    def f(x, y):
        return x**2 + y**2
    
    # 创建网格
    x = np.linspace(0, 1, 100)
    y = np.linspace(0, 1, 100)
    X, Y = np.meshgrid(x, y)
    
    # 计算函数值
    Z = f(X, Y)
    
    # 数值积分（矩形法则）
    dx = x[1] - x[0]
    dy = y[1] - y[0]
    integral = np.sum(Z) * dx * dy
    
    print(f"数值积分结果: {integral:.6f}")
    print(f"理论值: {2/3:.6f}")  # ∫∫(x²+y²)dxdy = 2/3
    
    return integral

# 有限差分法中的网格应用
def finite_difference_2d():
    # 二维拉普拉斯方程求解
    # ∇²u = 0
    
    # 网格参数
    nx, ny = 50, 50
    dx = dy = 0.1
    x = np.linspace(0, (nx-1)*dx, nx)
    y = np.linspace(0, (ny-1)*dy, ny)
    X, Y = np.meshgrid(x, y)
    
    # 初始化解
    u = np.zeros((nx, ny))
    
    # 边界条件
    u[0, :] = 0      # 左边界
    u[-1, :] = 0     # 右边界
    u[:, 0] = 0      # 下边界
    u[:, -1] = 100   # 上边界
    
    # 迭代求解
    for iteration in range(1000):
        u_old = u.copy()
        for i in range(1, nx-1):
            for j in range(1, ny-1):
                u[i, j] = 0.25 * (u_old[i+1, j] + u_old[i-1, j] + 
                                 u_old[i, j+1] + u_old[i, j-1])
        
        # 检查收敛
        if np.max(np.abs(u - u_old)) < 1e-6:
            print(f"收敛于第 {iteration} 次迭代")
            break
    
    # 可视化结果
    import matplotlib.pyplot as plt
    
    plt.figure(figsize=(12, 5))
    
    plt.subplot(1, 2, 1)
    plt.contourf(X, Y, u, levels=20, cmap='viridis')
    plt.colorbar()
    plt.title('拉普拉斯方程解')
    plt.xlabel('X')
    plt.ylabel('Y')
    
    plt.subplot(1, 2, 2)
    plt.plot(x, u[:, ny//2], 'b-', linewidth=2)
    plt.title('中心线剖面')
    plt.xlabel('X')
    plt.ylabel('u')
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()
    
    return u
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176

# 性能优化技巧

# 1. 内存优化

# 内存使用监控

import psutil
import os

def get_memory_usage():
    process = psutil.Process(os.getpid())
    return process.memory_info().rss / 1024 / 1024  # MB

def monitor_memory(func):
    def wrapper(*args, **kwargs):
        before = get_memory_usage()
        result = func(*args, **kwargs)
        after = get_memory_usage()
        print(f"内存使用: {after - before:.2f} MB")
        return result
    return wrapper

# 使用示例
@monitor_memory
def large_array_operation():
    arr = np.random.rand(10000, 10000)
    return np.sum(arr)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

# 内存映射

# 使用内存映射处理大文件
def process_large_file(filename, chunk_size=1000):
    data = np.memmap(filename, dtype='float64', mode='r', shape=(10000, 1000))
    
    results = []
    for i in range(0, len(data), chunk_size):
        chunk = data[i:i+chunk_size]
        result = np.mean(chunk, axis=0)
        results.append(result)
    
    return np.array(results)
1
2
3
4
5
6
7
8
9
10
11

# 2. 计算优化

# 向量化优化

# 避免循环，使用向量化操作
def vectorized_operation_example():
    # 慢的方法（循环）
    def slow_method(arr):
        result = np.zeros_like(arr)
        for i in range(len(arr)):
            result[i] = arr[i] * 2 + 1
        return result
    
    # 快的方法（向量化）
    def fast_method(arr):
        return arr * 2 + 1
    
    # 测试性能
    arr = np.random.rand(1000000)
    
    import time
    start = time.time()
    result1 = slow_method(arr)
    slow_time = time.time() - start
    
    start = time.time()
    result2 = fast_method(arr)
    fast_time = time.time() - start
    
    print(f"循环方法: {slow_time:.4f}秒")
    print(f"向量化方法: {fast_time:.4f}秒")
    print(f"加速比: {slow_time/fast_time:.1f}x")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

# 并行计算

# 使用多进程
from multiprocessing import Pool
import numpy as np

def parallel_processing_example():
    def process_chunk(chunk):
        return np.sum(chunk**2)
    
    # 创建大数据
    data = np.random.rand(1000000)
    chunk_size = 100000
    chunks = [data[i:i+chunk_size] for i in range(0, len(data), chunk_size)]
    
    # 并行处理
    with Pool() as pool:
        results = pool.map(process_chunk, chunks)
    
    return np.sum(results)

# 使用Numba加速
from numba import jit

@jit(nopython=True)
def numba_optimized_function(arr):
    result = np.zeros_like(arr)
    for i in range(len(arr)):
        result[i] = arr[i] * 2 + 1
    return result
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

# 3. 算法优化

# 分块计算

def block_matrix_operations(A, B, block_size=1000):
    """分块矩阵运算，减少内存使用"""
    m, n = A.shape
    p = B.shape[1]
    result = np.zeros((m, p))
    
    for i in range(0, m, block_size):
        for j in range(0, p, block_size):
            for k in range(0, n, block_size):
                result[i:i+block_size, j:j+block_size] += np.matmul(
                    A[i:i+block_size, k:k+block_size],
                    B[k:k+block_size, j:j+block_size]
                )
    
    return result
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

# 缓存优化

from functools import lru_cache

@lru_cache(maxsize=128)
def cached_computation(x):
    """缓存计算结果"""
    return np.sum(x**2)

# 使用示例
for i in range(100):
    result = cached_computation(tuple(range(i)))  # 需要转换为不可变类型
1
2
3
4
5
6
7
8
9
10

# 调试和测试技巧

# 1. 调试工具

# 数组检查工具

def array_debug_info(arr, name="array"):
    """打印数组的详细信息"""
    print(f"=== {name} 调试信息 ===")
    print(f"形状: {arr.shape}")
    print(f"数据类型: {arr.dtype}")
    print(f"内存大小: {arr.nbytes / 1024 / 1024:.2f} MB")
    print(f"最小值: {np.min(arr)}")
    print(f"最大值: {np.max(arr)}")
    print(f"均值: {np.mean(arr):.4f}")
    print(f"标准差: {np.std(arr):.4f}")
    print(f"是否包含NaN: {np.any(np.isnan(arr))}")
    print(f"是否包含Inf: {np.any(np.isinf(arr))}")
    print("=" * 30)

# 使用示例
arr = np.random.randn(1000, 1000)
array_debug_info(arr, "测试数组")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# 性能分析

import time
import cProfile
import pstats

def profile_function(func, *args, **kwargs):
    """分析函数性能"""
    profiler = cProfile.Profile()
    profiler.enable()
    
    start_time = time.time()
    result = func(*args, **kwargs)
    end_time = time.time()
    
    profiler.disable()
    stats = pstats.Stats(profiler)
    stats.sort_stats('cumulative')
    stats.print_stats(10)
    
    print(f"总执行时间: {end_time - start_time:.4f}秒")
    return result
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# 2. 单元测试

# NumPy测试框架

import numpy.testing as npt

def test_array_operations():
    """测试数组操作"""
    # 测试基本运算
    a = np.array([1, 2, 3])
    b = np.array([4, 5, 6])
    
    # 测试加法
    expected = np.array([5, 7, 9])
    npt.assert_array_equal(a + b, expected)
    
    # 测试乘法
    expected = np.array([4, 10, 18])
    npt.assert_array_equal(a * b, expected)
    
    # 测试统计函数
    npt.assert_almost_equal(np.mean(a), 2.0)
    npt.assert_almost_equal(np.std(a), 0.816496580927726)

def test_matrix_operations():
    """测试矩阵操作"""
    A = np.array([[1, 2], [3, 4]])
    B = np.array([[5, 6], [7, 8]])
    
    # 测试矩阵乘法
    expected = np.array([[19, 22], [43, 50]])
    npt.assert_array_equal(np.dot(A, B), expected)
    
    # 测试矩阵求逆
    A_inv = np.linalg.inv(A)
    identity = np.dot(A, A_inv)
    npt.assert_array_almost_equal(identity, np.eye(2))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

# 3. 错误处理

# 自定义异常

class NumPyError(Exception):
    """NumPy相关错误的基类"""
    pass

class ShapeMismatchError(NumPyError):
    """形状不匹配错误"""
    pass

class InvalidDataTypeError(NumPyError):
    """无效数据类型错误"""
    pass

def safe_array_operation(arr1, arr2, operation='add'):
    """安全的数组操作"""
    try:
        # 检查形状
        if arr1.shape != arr2.shape:
            raise ShapeMismatchError(f"形状不匹配: {arr1.shape} vs {arr2.shape}")
        
        # 检查数据类型
        if arr1.dtype != arr2.dtype:
            print(f"警告: 数据类型不同 {arr1.dtype} vs {arr2.dtype}")
        
        # 执行操作
        if operation == 'add':
            return arr1 + arr2
        elif operation == 'multiply':
            return arr1 * arr2
        else:
            raise ValueError(f"不支持的操作: {operation}")
            
    except Exception as e:
        print(f"操作失败: {e}")
        return None
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

# 常见陷阱和解决方案

# 1. 浮点数精度问题

# 问题：浮点数比较
a = 0.1 + 0.2
b = 0.3
print(a == b)  # False!

# 解决方案：使用numpy.testing
npt.assert_almost_equal(a, b, decimal=15)

# 或者使用相对误差
def is_close(a, b, rtol=1e-5, atol=1e-8):
    return np.abs(a - b) <= (atol + rtol * np.abs(b))

print(is_close(a, b))  # True
1
2
3
4
5
6
7
8
9
10
11
12
13

# 2. 内存泄漏

# 问题：大数组没有及时释放
def memory_leak_example():
    large_arrays = []
    for i in range(100):
        arr = np.random.rand(10000, 10000)
        large_arrays.append(arr)  # 内存泄漏！
    
    # 解决方案：及时释放
    del large_arrays
    import gc
    gc.collect()
1
2
3
4
5
6
7
8
9
10
11

# 3. 广播陷阱

# 问题：意外的广播
a = np.array([1, 2, 3])
b = np.array([[1], [2]])
try:
    result = a + b  # 可能产生意外结果
except ValueError as e:
    print(f"广播错误: {e}")

# 解决方案：显式控制广播
a_expanded = np.expand_dims(a, axis=0)  # 形状变为 (1, 3)
result = a_expanded + b  # 明确的广播
1
2
3
4
5
6
7
8
9
10
11

# 扩展资源

# 1. 相关库

  • SciPy: 科学计算库，基于NumPy
  • Pandas: 数据分析库，使用NumPy作为底层
  • Matplotlib: 绘图库，与NumPy紧密集成
  • Scikit-learn: 机器学习库，大量使用NumPy
  • Numba: JIT编译器，加速NumPy代码
  • Cython: 编译Python代码，提高性能

# 2. 学习资源

  • 官方文档: https://numpy.org/doc/
  • NumPy教程: https://numpy.org/doc/stable/user/quickstart.html
  • NumPy参考: https://numpy.org/doc/stable/reference/
  • NumPy中文文档: https://www.numpy.org.cn/

# 3. 社区资源

  • Stack Overflow: NumPy相关问题
  • GitHub: NumPy源码和示例
  • Reddit: r/Python和r/DataScience
  • Discord: Python和科学计算社区
#numpy#Python#data-science
上次更新: 2025/10/08, 16:24:59
pandas
matplotlib

← matplotlib→

最近更新
01
npx 使用指南
10-12
02
cursor
09-28
03
inspect
07-20
更多文章>
Theme by Vdoing | Copyright © 2019-2025 Jacky | MIT License