Brightway中高效存储采样矩阵值的技术方案问询

看起来你已经把LCA蒙特卡洛分析的基础铺垫工作都搞定了，接下来咱们一步步实现你要的这个可复用采样数组，完美适配敏感性分析这类场景：

核心思路 #

我们的目标是把每个Monte Carlo迭代的CSR稀疏矩阵的非零元素值提取出来，按「行=固定非零元素位置，列=迭代次数」的结构堆叠成二维数组——这样后续做敏感性分析时，直接取某一行就能拿到对应矩阵元素的所有采样值，效率拉满。

步骤1：验证稀疏矩阵结构一致性 #

首先要确认所有Monte Carlo迭代的CSR矩阵非零元素的位置是固定的（LCA的蒙特卡洛通常是参数波动，稀疏结构不会变，但还是要检查一下）：

import numpy as np
from scipy.sparse import csr_matrix

# 假设你的所有迭代CSR矩阵存在这个列表里
mc_csr_matrices = [your_iter_1_csr, your_iter_2_csr, ...]

# 检查所有矩阵的indptr和indices是否和第一个矩阵一致
first_mat = mc_csr_matrices[0]
structure_consistent = all(
    (mat.indptr == first_mat.indptr).all() and (mat.indices == first_mat.indices).all()
    for mat in mc_csr_matrices
)

if not structure_consistent:
    raise ValueError("部分迭代的CSR矩阵稀疏结构不一致，需要先对齐非零元素位置！")

步骤2：提取并堆叠采样值 #

直接提取每个CSR矩阵的data属性（就是非零元素值的数组），然后用column_stack把它们拼成二维数组：

# 提取所有迭代的非零元素值数组
mc_nonzero_values = [mat.data for mat in mc_csr_matrices]

# 水平堆叠成目标数组：行=非零元素位置，列=迭代次数
tech_matrix_samples = np.column_stack(mc_nonzero_values)

现在tech_matrix_samples完全符合你的需求：每一行对应技术圈矩阵中一个固定的非零元素，每一列是该元素在对应Monte Carlo迭代中的采样值。

步骤3：关联元素元数据（可选但实用） #

结合你已有的activity_dict、product_dict和项目映射字典，给每个非零元素加上业务层面的标识，后续分析时能快速定位：

# 从第一个矩阵获取所有非零元素的行列索引
rows, cols = first_mat.nonzero()

# 生成元数据列表
element_metadata = []
for row_idx, col_idx in zip(rows, cols):
    # 映射到实际的活动和产品
    activity = activity_dict.get(row_idx, f"未知活动_{row_idx}")
    product = product_dict.get(col_idx, f"未知产品_{col_idx}")
    # 用项目映射字典做自定义转换
    mapped_activity = your_project_mapping.get(activity, activity)
    mapped_product = your_project_mapping.get(product, product)
    
    element_metadata.append({
        "row_idx": row_idx,
        "col_idx": col_idx,
        "activity": mapped_activity,
        "product": mapped_product
    })

步骤4：高效保存与复用 #

用numpy的npz压缩格式保存数组和元数据，既省空间又方便加载：

# 保存数据
np.savez(
    "tech_matrix_mc_samples.npz",
    samples=tech_matrix_samples,
    element_metadata=np.array(element_metadata, dtype=object)
)

# 后续复用加载
loaded_data = np.load("tech_matrix_mc_samples.npz", allow_pickle=True)
tech_matrix_samples = loaded_data["samples"]
element_metadata = loaded_data["element_metadata"].tolist()

# 举个敏感性分析的例子：提取指定活动-产品对的所有采样值
target_item = next(
    item for item in element_metadata 
    if item["activity"] == "你的目标活动ID/名称"
)
target_samples = tech_matrix_samples[target_item["row_idx"], :]

注意事项 #

• 如果遇到稀疏结构不一致的迭代矩阵，需要先对齐所有矩阵到同一个稀疏模板（比如用scipy.sparse.csr_matrix的tocoo()格式统一后再转换回来）；
• 若迭代次数极多或非零元素数量巨大，可以考虑用dask.array做分块处理，避免内存溢出。

内容的提问来源于stack exchange，提问作者MPa