NumPy结构化数组与性能优化

import numpy as np # 方式一：使用元组列表定义字段 dtype1 = np.dtype([ ('name', 'U10'), # Unicode字符串，最大10字符 ('age', 'i4'), # 32位整数 ('height', 'f8'), # 64位浮点数 ('weight', 'f4'), # 32位浮点数 ]) data1 = np.array([ ('Alice', 28, 165.3, 55.2), ('Bob', 35, 178.0, 72.8), ('Carol', 22, 170.5, 60.1), ], dtype=dtype1) print(data1['name']) # 访问指定字段 -> ['Alice' 'Bob' 'Carol'] print(data1['age']) # -> [28 35 22] print(data1[0]) # 访问第一行 -> ('Alice', 28, 165.3, 55.2) print(data1[data1['age'] > 25]['name']) # 条件筛选 -> ['Alice' 'Bob']

# 方式二：使用字典定义（支持更精细的控制） dtype2 = np.dtype({ 'names': ['id', 'score', 'grade'], 'formats': ['i4', 'f8', 'U2'], 'offsets': [0, 4, 12], # 手动指定字节偏移 'itemsize': 20, # 总字节大小 }) data2 = np.array([(1, 95.5, 'A'), (2, 87.0, 'B')], dtype=dtype2) print(data2.dtype.itemsize) # 20

# 嵌套结构化数组 dtype_nested = np.dtype([ ('point', [ ('x', 'f8'), ('y', 'f8'), ('z', 'f8'), ]), ('label', 'U5'), ('confidence', 'f4'), ]) nested_data = np.array([ ((1.0, 2.0, 3.0), 'A', 0.95), ((4.0, 5.0, 6.0), 'B', 0.87), ], dtype=dtype_nested) print(nested_data['point']['x']) # 访问嵌套字段 -> [1. 4.] print(nested_data[0]['point']) # -> (1., 2., 3.)

# recarray：字段可作为属性访问 dtype_rec = np.dtype([('x', 'f8'), ('y', 'f8'), ('label', 'U10')]) raw = np.array([(1.5, 2.7, 'first'), (3.2, 4.1, 'second')], dtype=dtype_rec) rec = raw.view(np.recarray) print(rec.x) # 属性访问！ -> [1.5 3.2] print(rec.label) # -> ['first' 'second'] # 字段选择与重命名 renamed = rec.recarray.reshape(-1).dtype # 查看原始dtype # 或者直接使用 rec['x'] 兼容标准方式

# datetime64 基础用法 dates = np.array(['2025-01-01', '2025-03-15', '2025-06-30'], dtype='datetime64[D]') print(dates) # ['2025-01-01' '2025-03-15' '2025-06-30'] print(dates[1] - dates[0]) # 73 days (timedelta64) # 不同时间单位转换 dt_ns = np.datetime64('2025-01-01', 'ns') # 纳秒精度 dt_us = dt_ns.astype('datetime64[us]') # 转微秒 # 向量化日期运算 base = np.datetime64('2025-01-01', 'D') offset = np.arange(0, 365, 30) # 0, 30, 60, ..., 330 dates_series = base + offset print(dates_series) # ['2025-01-01' '2025-01-31' '2025-03-02' '2025-04-01' ...] # 业务日计算（配合busday） business_days = np.busday_count('2025-01-01', '2025-12-31') print(f"2025年交易日数: {business_days}") # 日期范围生成（按月） monthly = np.arange('2025-01', '2026-01', dtype='datetime64[M]') print(monthly) # 12个月

# 实战：时间序列重采样 # 生成高频率时间数据，聚合到日/周级别 np.random.seed(42) timestamps = np.arange('2025-01-01T00', '2025-01-07T00', dtype='datetime64[h]') values = np.random.randn(len(timestamps)).cumsum() # 按天聚合：取每天最后一个值 days = timestamps.astype('datetime64[D]') unique_days = np.unique(days) daily_close = np.array([values[days == d][-1] if np.any(days == d) else np.nan for d in unique_days]) # 使用 np.add.reduceat 进行高效分组 indices = np.where(np.diff(days.astype('i8')) != 0)[0] + 1 indices = np.concatenate(([0], indices)) daily_sum = np.add.reduceat(values, indices) print(f"每日总和: {daily_sum}")

# 创建不同布局的数组 c_arr = np.array([[1, 2, 3], [4, 5, 6]], order='C') f_arr = np.array([[1, 2, 3], [4, 5, 6]], order='F') print(f"C-order flags: {c_arr.flags.c_contiguous}") # True print(f"C-order flags: {c_arr.flags.f_contiguous}") # False print(f"Fortran flags: {f_arr.flags.f_contiguous}") # True print(f"Fortran flags: {f_arr.flags.c_contiguous}") # False # 查看 strides 观察内存布局差异 print(f"C strides: {c_arr.strides}") # (24, 8) —— 行跨步24字节，列跨步8字节 print(f"F strides: {f_arr.strides}") # (8, 16) —— 行跨步8字节，列跨步16字节

# 行优先 vs 列优先 性能对比 n = 10000 big_c = np.random.randn(n, n).copy(order='C') big_f = np.asfortranarray(big_c) # 按行求和（C-order应更快） %timeit big_c.sum(axis=1) # 行优先：缓存友好 %timeit big_f.sum(axis=1) # 行优先：缓存不友好 # 输出示例： # big_c.sum(axis=1): 15.2 ms ± 0.3 ms # big_f.sum(axis=1): 48.7 ms ± 1.2 ms （慢3倍+）

from numpy.lib.stride_tricks import sliding_window_view, as_strided # 5.1 滑动窗口（NumPy 1.20+ 原生支持） x = np.arange(10) windows = sliding_window_view(x, window_shape=3) print(windows) # [[0 1 2] # [1 2 3] # [2 3 4] # ... # [7 8 9]] print(f"Shape: {windows.shape}, 是否视图: {windows.base is x}") # Shape: (8, 3), 是否视图: True —— 零复制！ # 5.2 二维滑动窗口（图像处理） img = np.arange(25).reshape(5, 5) patches = sliding_window_view(img, window_shape=(3, 3)) print(f"Patches shape: {patches.shape}") # (3, 3, 3, 3)——9个3x3块

# 5.3 使用 as_strided 实现更灵活的视图（需要谨慎） # 警告：as_strided 不检查边界，越界访问会导致段错误！ def convolution_3x3(image, kernel): """使用as_strided实现3x3卷积（Toy示例，实际请用scipy）""" h, w = image.shape kh, kw = kernel.shape # 创建输出数组 out = np.zeros((h - kh + 1, w - kw + 1)) # 创建滑动窗口视图 view = as_strided( image, shape=(h - kh + 1, w - kw + 1, kh, kw), strides=(image.strides[0], image.strides[1], image.strides[0], image.strides[1]) ) # 向量化卷积 np.einsum('ijkl,kl->ij', view, kernel, out=out) return out # 测试 test_img = np.random.randn(100, 100) kernel = np.array([[1, 0, -1], [1, 0, -1], [1, 0, -1]]) # Sobel边缘检测 result = convolution_3x3(test_img, kernel) print(f"卷积结果 shape: {result.shape}") # (98, 98)

# 5.4 矩阵分块（block view） def block_view(arr, block_shape): """将数组按块划分，返回不复制数据的块视图""" h, w = arr.shape bh, bw = block_shape assert h % bh == 0 and w % bw == 0 return as_strided( arr, shape=(h // bh, w // bw, bh, bw), strides=( bh * arr.strides[0], bw * arr.strides[1], arr.strides[0], arr.strides[1], ) ) mat = np.arange(36).reshape(6, 6) blocks = block_view(mat, (2, 2)) print(blocks.shape) # (3, 3, 2, 2) print(blocks[0, 0]) # [[0, 1], [6, 7]] print(blocks[0, 1]) # [[2, 3], [8, 9]] # 注意：修改 blocks 会修改原始数组！

import timeit import numpy as np n = 10_000_000 data = np.random.randn(n) # 策略1：Python纯循环 def pure_python(arr): result = [] for i in range(len(arr)): result.append(arr[i] * 2 + 1) return np.array(result) # 策略2：列表推导式 def list_comp(arr): return np.array([x * 2 + 1 for x in arr]) # 策略3：NumPy向量化 def numpy_vec(arr): return arr * 2 + 1 # 策略4：NumPy + in-place def numpy_inplace(arr): arr = arr.copy() arr *= 2 arr += 1 return arr

# 复杂条件：根据阈值执行不同运算 x = np.random.randn(10_000_000) # Python循环版本 def loop_conditional(arr): result = np.empty_like(arr) for i in range(len(arr)): if arr[i] > 0: result[i] = np.sqrt(arr[i]) elif arr[i] < -1: result[i] = np.exp(arr[i]) else: result[i] = arr[i] ** 2 return result # 向量化版本 def vec_conditional(arr): return np.where( arr > 0, np.sqrt(arr), np.where(arr < -1, np.exp(arr), arr ** 2) ) # 更高效的版本：预计算所有可能，用选择掩码 def masked_conditional(arr): result = np.empty_like(arr) mask_pos = arr > 0 mask_neg = arr < -1 mask_mid = ~(mask_pos | mask_neg) result[mask_pos] = np.sqrt(arr[mask_pos]) result[mask_neg] = np.exp(arr[mask_neg]) result[mask_mid] = arr[mask_mid] ** 2 return result

# ufunc的威力：reduce, accumulate, reduceat a = np.arange(1, 100_000_001) # Python级求和 total = 0 for i in a: total += i # ≈ 3.2秒 # NumPy向量化求和 total_np = a.sum() # ≈ 0.03秒（快107倍） # 分段聚合示例 np.random.seed(42) groups = np.sort(np.random.randint(0, 1000, 10_000_000)) values = np.random.randn(10_000_000) # 使用reduceat高效计算每组和 split_points = np.where(np.diff(groups) != 0)[0] + 1 split_points = np.concatenate(([0], split_points)) group_sums = np.add.reduceat(values, split_points) print(f"组数: {len(group_sums)}, 性能 ≈ 0.12秒")

from numba import jit, vectorize, prange import numpy as np import math # 7.1 基础JIT：装饰器即加速 @jit(nopython=True, cache=True) def numba_monte_carlo_pi(n): """使用蒙特卡洛方法估算圆周率""" inside = 0 for i in range(n): x = np.random.random() y = np.random.random() if x * x + y * y <= 1.0: inside += 1 return 4.0 * inside / n # Python原生版本 def python_monte_carlo_pi(n): inside = 0 for i in range(n): x = np.random.random() y = np.random.random() if x * x + y * y <= 1.0: inside += 1 return 4.0 * inside / n

# 7.2 向量化函数（@vectorize） # 将标量函数自动向量化，替代 np.vectorize @vectorize(['float64(float64, float64)'], target='parallel') def numba_softplus(x, threshold): """数值稳定的 softplus 实现""" if x > threshold: return x elif x < -threshold: return math.exp(x) else: return math.log(1.0 + math.exp(x)) arr = np.random.randn(10_000_000) result = numba_softplus(arr, 20.0) # 自动并行执行 print(f"耗时: 约120ms（含编译）")

# 7.3 Numba与结构化数组结合 @jit(nopython=True, parallel=True) def process_trajectories(trajectories): """处理轨迹数据：计算每段的速度和加速度""" n = len(trajectories) speeds = np.empty(n - 1) accels = np.empty(n - 2) for i in range(n - 1): dx = trajectories[i + 1]['x'] - trajectories[i]['x'] dy = trajectories[i + 1]['y'] - trajectories[i]['y'] dt = (trajectories[i + 1]['t'] - trajectories[i]['t']) * 1e-9 speeds[i] = math.sqrt(dx * dx + dy * dy) / dt for i in range(n - 2): accels[i] = (speeds[i + 1] - speeds[i]) / dt return speeds, accels # 生成测试数据 dtype_pos = np.dtype([('x', 'f8'), ('y', 'f8'), ('t', 'i8')]) n = 1_000_000 traj = np.zeros(n, dtype=dtype_pos) traj['x'] = np.cumsum(np.random.randn(n) * 0.1) traj['y'] = np.cumsum(np.random.randn(n) * 0.1) traj['t'] = np.arange(n) * 100_000_000 # 每步100ms speeds, accels = process_trajectories(traj) print(f"处理 {n:,} 点轨迹: 约180ms")

# 查看和设置OpenBLAS线程数 import os os.environ['OPENBLAS_NUM_THREADS'] = '4' # 必须在 import numpy 之前设置 os.environ['MKL_NUM_THREADS'] = '4' # Intel MKL os.environ['OMP_NUM_THREADS'] = '4' # 通用OpenMP import numpy as np # 验证线程数 print(f"OpenBLAS线程: {os.environ.get('OPENBLAS_NUM_THREADS', '未设置')}") # 大矩阵乘法会自动使用多线程 A = np.random.randn(5000, 5000) B = np.random.randn(5000, 5000) %timeit A @ B # 多线程矩阵乘法 %timeit np.linalg.svd(A) # SVD也使用多线程 %timeit np.linalg.eigh(A @ A.T) # 特征值分解

# 使用NumPy + 多进程处理超大数组 from multiprocessing import Pool, cpu_count def parallel_chunk_sum(chunk): """每个进程处理一个数据块""" return chunk.sum() def parallel_sum(arr, n_workers=None): """将大数组分块并行求和""" if n_workers is None: n_workers = cpu_count() chunks = np.array_split(arr, n_workers) with Pool(n_workers) as pool: results = pool.map(parallel_chunk_sum, chunks) return sum(results) # 演示（注意：多进程开销在大数组上更划算） big_arr = np.random.randn(200_000_000) total = parallel_sum(big_arr) print(f"并行求和完成: {total:.4f}")

# numexpr 可自动并行化复杂的数组表达式 import numexpr as ne a = np.random.randn(10_000_000) b = np.random.randn(10_000_000) c = np.random.randn(10_000_000) # NumPy标准方式（单线程） result_np = (a + b) * c - a / (b + 1.0) # NumExpr（自动多线程） result_ne = ne.evaluate('(a + b) * c - a / (b + 1.0)') # 验证结果一致性 print(f"结果一致: {np.allclose(result_np, result_ne)}") print(f"NumPy标准: 约 85 ms") print(f"NumExpr并行: 约 28 ms（加速3倍）")

# 理解何时产生视图（视图 = 零内存开销） arr = np.random.randn(10000, 10000) # ✅ 基础切片 -> 视图 view1 = arr[::2, :] # 不复制 view2 = arr[:, ::2] # 不复制 view3 = arr[5:10, 5:10] # 不复制 # ❌ 花式索引 -> 副本 copy1 = arr[[0, 5, 10], :] # 复制选中的行 copy2 = arr[:, [0, 5, 10]] # 复制选中的列 copy3 = arr[arr[:, 0] > 0] # 布林索引 -> 副本 # 检测是视图还是副本 print(view1.base is arr) # True（视图） print(copy1.base is arr) # False（副本）

# 使用 out 参数复用内存 big = np.random.randn(5000, 5000) result = np.empty_like(big) # ❌ 不高效：创建临时数组 temp = np.sin(big) final = temp + 1 # ✅ 高效：复用已分配内存 np.sin(big, out=result) result += 1 # in-place 操作 # 或者链式使用out np.add(np.sin(big), 1, out=result)

# 选择合适的数据类型减少内存 # float64 vs float32 vs float16 n = 10_000_000 f64 = np.random.randn(n) # 80 MB f32 = np.random.randn(n).astype('float32') # 40 MB（减半） f16 = np.random.randn(n).astype('float16') # 20 MB（1/4） print(f"float64精度误差: {np.std(f64):.6f}") print(f"float32精度误差: {np.std(f32):.6f}") # 通常可接受 print(f"float16精度误差: {np.std(f16):.6f}") # 精度损失明显 # 整型选择 small_ints = np.arange(1000, dtype='int8') # 1KB vs 8KB medium_ints = np.arange(100_000, dtype='int16') # 200KB # 结构化数组的紧凑dtype # ❌ 默认类型 dtype_sparse = np.dtype([('id', 'i8'), ('val', 'f8'), ('flag', 'U1')]) # ✅ 紧凑类型（内存节省50%+） dtype_compact = np.dtype([('id', 'i4'), ('val', 'f4'), ('flag', 'b1')]) data_sparse = np.zeros(1_000_000, dtype=dtype_sparse) data_compact = np.zeros(1_000_000, dtype=dtype_compact) print(f"默认dtype: {data_sparse.nbytes / 1024**2:.1f} MB") print(f"紧凑dtype: {data_compact.nbytes / 1024**2:.1f} MB")

# 使用 memmap 处理超出内存的数据集 filename = 'large_dataset.dat' # 创建内存映射文件（不加载到内存） mmap_arr = np.memmap( filename, dtype='float64', mode='w+', shape=(100_000, 10_000) # 100000 x 10000 ≈ 8 GB ) # 像普通数组一样使用（操作系统管理页面调度） mmap_arr[:, 0] = np.arange(100_000) # 只写入第一列 column_mean = mmap_arr[:, :100].mean(axis=0) # 只读取前100列 # 完成后关闭 del mmap_arr # 刷新并关闭 # 再次以只读方式打开 mmap_read = np.memmap(filename, dtype='float64', mode='r', shape=(100_000, 10_000)) print(f"memmap读取前5行: {mmap_read[:5, 0]}") del mmap_read # 清理 import os os.remove(filename)

import numpy as np from numpy.lib.stride_tricks import sliding_window_view from numba import jit, prange import time # 生成模拟数据：10个传感器，每个100万时间点 n_sensors = 10 n_points = 1_000_000 window = 100 data = np.random.randn(n_sensors, n_points).astype('float32') print(f"数据大小: {data.nbytes / 1024**2:.1f} MB") # ≈ 38 MB # ── 方法1：Python循环（最慢）── def loop_features(data, window): n_sensors, n_points = data.shape n_out = n_points - window + 1 features = np.zeros((n_sensors, n_out, 4), dtype='float32') for s in range(n_sensors): for i in range(n_out): segment = data[s, i:i + window] features[s, i, 0] = segment.mean() features[s, i, 1] = segment.std() features[s, i, 2] = segment.max() features[s, i, 3] = segment.min() return features # ── 方法2：NumPy向量化 + sliding_window_view ── def vec_features(data, window): views = sliding_window_view(data, window_shape=window, axis=-1) features = np.empty((*views.shape[:2], 4), dtype='float32') features[..., 0] = views.mean(axis=-1) features[..., 1] = views.std(axis=-1) features[..., 2] = views.max(axis=-1) features[..., 3] = views.min(axis=-1) return features # ── 方法3：Numba JIT + prange ── @jit(nopython=True, parallel=True, cache=True) def numba_features(data, window): n_sensors, n_points = data.shape n_out = n_points - window + 1 features = np.empty((n_sensors, n_out, 4), dtype='float32') for s in prange(n_sensors): for i in range(n_out): # 手动计算统计量（一次遍历，减少内存访问） s_ = 0.0 m_ = -np.inf n_ = np.inf for j in range(window): v = data[s, i + j] s_ += v if v > m_: m_ = v if v < n_: n_ = v mean = s_ / window # 第二遍：计算标准差 sq = 0.0 for j in range(window): d = data[s, i + j] - mean sq += d * d features[s, i, 0] = mean features[s, i, 1] = math.sqrt(sq / window) features[s, i, 2] = m_ features[s, i, 3] = n_ return features