numpy - Python粒子模拟器:核心外处理

Question

问题描述

用python / numpy编写蒙特卡洛粒子模拟器（布朗运动和光子发射）。我需要将模拟输出（>> 10GB）保存到文件中，然后在第二步中处理数据。与Windows和Linux的兼容性很重要。

颗粒数（n_particles）是10-100。时间步数（time_size）为〜10 ^ 9。

仿真包括3个步骤（以下代码适用于全内存版本）：


模拟（并存储）一个emission速率数组（包含许多几乎为0的元素）：


形状（n_particles x time_size），float32，大小80GB

计算counts数组，（来自泊松过程的随机值，具有先前计算的速率）：


形状（n_particles x time_size），uint8，大小20GB

counts = np.random.poisson(lam=emission).astype(np.uint8)

查找计数的时间戳（或索引）。计数几乎总是为0，因此时间戳数组将适合RAM。

# Loop across the particles
timestamps = [np.nonzero(c) for c in counts]

我只执行一次步骤1，然后重复执行步骤2-3多次（约100次）。将来我可能需要在计算 emission之前对 cumsum进行预处理（应用 counts或其他函数）。

题

我有一个有效的内存中实现，并且我试图了解什么是实现可以扩展到（很多）更长的仿真的核外版本的最佳方法。

我希望它存在

我需要将数组保存到文件中，并且我想使用单个文件进行仿真。我还需要一种“简单”的方式来存储和调用仿真参数字典（标量）。

理想情况下，我想要一个文件支持的numpy数组，它可以预分配并填充块。然后，我希望numpy数组方法（ max， cumsum，...）透明地工作，只需要 chunksize关键字来指定每次迭代要加载多少数组。

更好的是，我希望Numexpr不在缓存和RAM之间运行，而在RAM和硬盘驱动器之间运行。

实际的选择是什么

作为首选
我开始尝试pyTables，但对其复杂性和抽象（与numpy不同）感到不满意。而且，我当前的解决方案（阅读下文）非常丑陋，效率不高。

所以我寻求答案的选择是


实现具有所需功能的numpy数组（如何？）
以更聪明的方式使用pytable（不同的数据结构/方法）
使用另一个库：h5py，blaze，pandas ...（到目前为止，我还没有尝试过）。


临时解决方案（pyTables）

我将模拟参数保存在 '/parameters'组中：每个参数都转换为numpy数组标量。详细的解决方案，但它可以工作。

我将 emission保存为可扩展数组（ EArray），因为我以块的形式生成数据，并且需要附加每个新块（尽管我知道最终的大小）。保存 counts更有问题。如果将其保存为pytable数组，则很难执行“ counts> = 2”之类的查询。因此，我将计数保存为多个表（每个粒子一个）[UGLY]，并使用 .get_where_list('counts >= 2')查询。我不确定这是否节省空间，并且
生成所有这些表而不是使用单个数组，会使HDF5文件变得更加混乱。此外，奇怪的是，创建这些表需要创建一个自定义dtype（即使对于标准numpy dtype）：

    dt = np.dtype([('counts', 'u1')])        
    for ip in xrange(n_particles):
        name = "particle_%d" % ip
        data_file.create_table(
                    group, name, description=dt, chunkshape=chunksize,
                    expectedrows=time_size,
                    title='Binned timetrace of emitted ph (bin = t_step)'
                        ' - particle_%d' % particle)

每个粒子计数“表”都有一个不同的名称（ name = "particle_%d" % ip），为了方便迭代，我需要将它们放在python列表中。

编辑：这个问题的结果是一个称为 PyBroMo的布朗运动模拟器。

Answer 1

Dask.array可以在磁盘阵列（如PyTables或h5py）上执行诸如max，cumsum等的分块操作。

import h5py
d = h5py.File('myfile.hdf5')['/data']
import dask.array as da
x = da.from_array(d, chunks=(1000, 1000))

X的外观就像一个numpy数组，并复制了许多API。 x上的操作将创建一个内存中操作的DAG，将在必要时使用磁盘中的多核流有效执行

da.exp(x).mean(axis=0).compute()

http://dask.pydata.org/en/latest/

conda install dask
or 
pip install dask