python执行效率有多低,基于numpy的计算效率低下的多重处理

I'm trying to parallelize some calculations that use numpy with the help of Python's multiprocessing module. Consider this simplified example:

import time

import numpy

from multiprocessing import Pool

def test_func(i):

a = numpy.random.normal(size=1000000)

b = numpy.random.normal(size=1000000)

for i in range(2000):

a = a + b

b = a - b

a = a - b

return 1

t1 = time.time()

test_func(0)

single_time = time.time() - t1

print("Single time:", single_time)

n_par = 4

pool = Pool()

t1 = time.time()

results_async = [

pool.apply_async(test_func, [i])

for i in range(n_par)]

results = [r.get() for r in results_async]

multicore_time = time.time() - t1

print("Multicore time:", multicore_time)

print("Efficiency:", single_time / multicore_time)

When I execute it, the multicore_time is roughly equal to single_time * n_par, while I would expect it to be close to single_time. Indeed, if I replace numpy calculations with just time.sleep(10), this is what I get — perfect efficiency. But for some reason it does not work with numpy. Can this be solved, or is it some internal limitation of numpy?

Some additional info which may be useful:

I'm using OSX 10.9.5, Python 3.4.2 and the CPU is Core i7 with (as reported by the system info) 4 cores (although the above program only takes 50% of CPU time in total, so the system info may not be taking into account hyperthreading).

when I run this I see n_par processes in top working at 100% CPU

if I replace numpy array operations with a loop and per-index operations, the efficiency rises significantly (to about 75% for n_par = 4).

解决方案

It looks like the test function you're using is memory bound. That means that the run time you're seeing is limited by how fast the computer can pull the arrays from memory into cache. For example, the line a = a + b is actually using 3 arrays, a, b and a new array that will replace a. These three arrays are about 8MB each (1e6 floats * 8 bytes per floats). I believe the different i7s have something like 3MB - 8MB of shared L3 cache so you cannot fit all 3 arrays in cache at once. Your cpu adds the floats faster than the array can be loaded into cache so most of the time is spent waiting on the array to be read from memory. Because the cache is shared between the cores, you don't see any speedup by spreading the work onto multiple cores.

Memory bound operations are an issue for numpy in general and the only way I know to deal with them is to use something like cython or numba.

import rasterio import numpy as np import json import glob from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import accuracy_score, cohen_kappa_score # 国家一级分类编码(文件名需包含类别编号class-6.tif) CLASS_MAP = { 1: "耕地", 2: "林地", 3: "草地", 4: "工矿仓储用地", 5: "住宅用地", 6: "交通用地", 7: "水域", 8: "未利用地" } def load_labels(img_shape, class_dir='labels/'): """ 合并类别标签文件 :param img_shape: 原始影像的(height, width) :param class_dir: 存放类别TIFF的目录 :return: 合并后的标签矩阵 (height, width) """ labels = np.zeros(img_shape, dtype=np.uint8) conflict_count = 0 # 遍历所有类别文件 for tif_path in glob.glob(f"{class_dir}class_*.tif"): class_id = int(tif_path.split('_')[-1].split('.')[0]) # 从文件名提取编号 with rasterio.open(tif_path) as src: class_mask = src.read(1) # 读取单波段 # 验证空间对齐 if class_mask.shape != img_shape: raise ValueError(f"标签文件{tif_path}尺寸{class_mask.shape}与影像不匹配") # 合并逻辑(假设标签文件为二值掩膜) overlap = (labels > 0) & (class_mask > 0) conflict_count += np.sum(overlap) labels[class_mask > 0] = class_id # 覆盖写入 print(f"警告:发现{conflict_count}个像素的重标注") return labels # 1. 数据预处理 with rasterio.open('JL1_data.tif') as src: img = src.read([3,2,1,4]) # B,G,R,NIR波段 meta = src.profile height, width = src.shape # 计算NDVI特征 nir, red = img[3].astype(float), img[0].astype(float) ndvi = (nir - red) / (nir + red + 1e-6) features = np.stack([img[0], img[1], img[2], img[3], ndvi], axis=0) labels = load_labels((height, width), class_dir='train_labels/') # 提取有效样本 mask = labels > 0 X = features[:, mask].T # 特征矩阵 (n_samples, 5) y = labels[mask] # 标签向量 # 3. 监督分类模型 X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3) model = RandomForestClassifier(n_estimators=100, oob_score=True) model.fit(X_train, y_train) # 4. 全图预测与保存 height, width = labels.shape full_pred = np.zeros((height, width), dtype=np.uint8) full_pred[mask] = model.predict(X) # 仅预测已标注区域 # 保存分类结果TIFF meta.update(dtype='uint8', count=1) with rasterio.open('classification.tif', 'w', **meta) as dst: dst.write(full_pred, 1) # 5. 精度评估 y_pred = model.predict(X_test) report = { "overall_accuracy": accuracy_score(y_test, y_pred), "kappa": cohen_kappa_score(y_test, y_pred), "class_mapping": CLASS_MAP } with open('accuracy_report.json', 'w') as f: json.dump(report, f, ensure_ascii=False)要用上述代码实现土地利用分类,看一下上面代码有没有错误,保证他能运行,并告诉我需要自己修改的输入路径在哪里
最新发布
03-30
<think>好的,我现在需要帮助用户检查他们的土地利用分类代码是否存在错误,并指出需要修改的输入路径。首先,我需要通读用户提供的代码,逐步分析可能的问题。 首先,用户使用了rasterio来读取遥感影像,这可能涉及到波段数据的处理。在数据预处理部分,他们用`src.read([3,2,1,4])`来读取波段,这里可能需要确认波段的顺序是否正确。例如,常见的顺序可能是红、绿、蓝、近红外,但不同数据源可能不同,用户需要确认自己的数据波段排列是否符合预期。 接下来是加载标签的部分,用户定义了一个`load_labels`函数。这里需要注意的是,标签文件的命名方式。代码中通过`glob.glob(f"{class_dir}class_*.tif")`查找文件,并从中提取类别ID。用户提到文件名应为“class-6.tif”,但代码中使用的是下划线分割,例如“class_6.tif”。这里可能存在文件名解析错误,需要检查文件名是否符合代码中的解析逻辑,否则会导致类别ID提取错误,从而影响标签的合并。 在合并标签时,用户使用覆盖写入的方式,这可能导致标签冲突。例如,如果一个像素在个标签文件中被标记,后面的类别会覆盖前面的,这可能不符合实际情况。用户需要确认标签是否互斥,或者是否需要处理这种冲突情况。当前代码统计了冲突数量,但处理方式是直接覆盖,可能需要用户注意这一点是否符合他们的数据情况。 在特征提取部分,用户计算了NDVI,这里需要注意分母是否为0的情况。代码中添加了1e-6来避免除以零,这应该是正确的。特征堆叠时,用户使用了`np.stack`,可能需要确认各波段的顺序是否正确,以及NDVI是否正确添加到了特征矩阵中。 数据拆分部分,用户使用`train_test_split`,但随机森林通常不需要分训练集和测试集来训练,因为可以使用袋外分数(OOB)。这里用户同时使用了测试集和OOB分数,可能重复评估,但不算错误,只是效率问题。可以建议用户直接使用OOB分数,减少计算时间。 全图预测时,用户使用`model.predict(X)`,这里的X是所有的有效样本,包括训练集和测试集。这可能存在问题,因为模型在预测时应该使用未见过的数据。但根据代码,用户可能在预测时仅使用已标注区域,这需要确认是否符合实际需求。此外,全图预测可能需要处理未标注区域,但用户代码中`mask`仅包含已标注区域,可能实际应用中需要预测整个影像,而不仅仅是已标注部分,这里可能存在逻辑错误。 保存结果时,用户更新了meta信息,并写入了分类结果。需要确认meta的dtype是否正确设置为uint8,以及count是否为1,这些与输出数据的一致性相关。 精度评估部分,用户计算了总体精度和Kappa系数,但缺少了每个类别的生产者精度和用户精度,这可能会影响评估的全面性。可以建议用户添加更详细的评估指标,比如分类报告。 关于输入路径,用户需要修改的地方包括输入的影像文件路径'JL1_data.tif'和标签目录'class_dir='train_labels/''。需要确认这些路径是否存在,以及标签文件是否存放在指定目录下,并且文件名符合要求。 此外,代码中可能存在一些潜在问题,比如标签合并时的覆盖逻辑、测试集的使用方式、全图预测的范围等,需要逐一指出并给出修改建议。同时,确保用户了解需要根据实际情况调整这些路径,并检查标签文件的命名和内容是否正确。</think>以下是对代码的逐项检查与修正建议(需特别注意⚠️标记部分): --- ### **一、路径设置检查(用户必须修改部分)** ```python # ⚠️需要用户自行确认的路径(红色为关键修改点): 1. 输入影像路径: 'JL1_data.tif' 2. 标签目录路径: class_dir='train_labels/' (需确保该目录包含class_*.tif文件) 3. 输出文件路径: 'classification.tif' 和 'accuracy_report.json' ``` --- ### **二、代码逻辑修正** #### 1. **标签加载问题** ```python # 原代码文件名解析逻辑 tif_path.split('_')[-1].split('.')[0] # 要求文件名必须为 class_6.tif # ⚠️若实际文件名为class-6.tif需改为: tif_path.split('-')[-1].split('.')[0] # 匹配连字符分隔的命名 ``` #### 2. **波段索引错误** ```python # 原代码波段读取顺序 img = src.read([3,2,1,4]) # 假设原始波段顺序为 [B1,B2,B3,B4] # ⚠️验证实际波段对应关系: # 常见卫星数据波段对照表 # | 索引 | Landsat-8波段 | 典型用途 | # |------|--------------|------------| # | 0 | Coastal | 非标准配置 | # | 1 | Blue | 水体识别 | # | 2 | Green | 植被 | # | 3 | Red | 地物边界 | # | 4 | NIR | 植被指数 | ``` **建议修改**:根据实际数据文档确认波段顺序 #### 3. **NDVI计算优化** ```python # 原分母处理存在精度风险 nir + red + 1e-6 → 建议改为 np.where((nir + red) == 0, 1e-6, nir + red) ``` #### 4. **全图预测逻辑错误** ```python # 原代码仅预测已标注区域(失去分类意义) full_pred[mask] = model.predict(X) # ❌ 错误逻辑 # ✅ 正确做法:预测全图所有像素 X_full = features.reshape(5, -1).T # 将(height, width, 5)转为(n_pixels, 5) full_pred = model.predict(X_full).reshape(height, width) ``` --- ### **三、潜在问题预警** | 问题类型 | 具体表现 | 解决方案 | |-----------------|-------------------------|---------------------------| | 内存不足 | 大型影像导致程序崩溃 | 分块处理数据 | | 类别不平衡 | 小样本类别精度低下 | 添加class_weight='balanced' | | 坐标系统不匹配 | 标签与影像无法对齐 | 检查rasterio的transform属性 | --- ### **四、改进建议** ```python # 在模型初始化处添加类别平衡 model = RandomForestClassifier( n_estimators=100, class_weight='balanced', # ← 新增 oob_score=True ) # 增加更详细的评估指标 from sklearn.metrics import classification_report print(classification_report(y_test, y_pred, target_names=CLASS_MAP.values())) ``` --- ### **五、代码验证检查表** 1. [ ] 确认`train_labels/`目录存在且包含class_1.tif到class_8.tif 2. [ ] 验证影像文件JL1_data.tif的波段顺序 3. [ ] 检查输出目录的写入权限 4. [ ] 确保numpy版本≥1.18(支持stack函数) **执行顺序建议**:先在小范围测试(可通过设置`height, width = 1000, 1000`截取部分数据),确认无误后再进行全图处理。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值