Python粒子模拟器:核外处理

Question

问题描述

在python/numpy中编写蒙特卡罗粒子模拟器(布朗运动和光子发射).我需要将模拟输出(>> 10GB)保存到文件中,然后在第二步中处理数据.与Windows和Linux的兼容性非常重要.

粒子数(n_particles)为10-100.时间步数(time_size)的数量是~10 ^ 9.

模拟有3个步骤(下面的代码是针对全内存版本):

1. 模拟(并存储)一个emission速率数组(包含许多几乎为0的元素):

   • shape(n_particlesx time_size),float32,大小80GB

2. 计算counts数组,(来自泊松过程的随机值,具有先前计算的速率):

   • shape(n_particlesx time_size),uint8,大小20GB

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

3. 查找计数的时间戳(或索引).计数几乎总是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不同)不满意.此外,我目前的解决方案(如下所示)是UGLY,效率不高.

所以我寻求答案的选择是

1. 实现一个具有所需功能的numpy数组(如何？)

2. 以更智能的方式使用pytable(不同的数据结构/方法)

3. 使用另一个库:h5py,blaze,pandas ......(到目前为止我还没有尝试过任何一个).

暂定解决方案(pyTables)

我将模拟参数保存在'/parameters'组中:每个参数都转换为numpy数组标量.详细的解决方案,但它的工作原理.

我保存emission为Extensible array(EArray),因为我以块的形式生成数据,我需要追加每个新块(我知道最终的大小).储蓄counts更成问题.如果将它保存为pytable数组,则很难执行"计数> = 2"之类的查询.因此,我将计数保存为多个表(每个粒子一个)[UGLY],我查询.get_where_list('counts >= 2').我不确定这是否节省空间,并且生成所有这些表而不是使用单个数组,显着地破坏了HDF5文件.而且,奇怪的是,创建这些表需要创建自定义dtype(即使对于标准的numpy dtypes):

    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

PyTable 解决方案

由于不需要 Pandas 提供的功能，并且处理速度要慢得多（请参见下面的笔记本），因此最好的方法是直接使用 PyTables 或 h5py。到目前为止我只尝试过 pytables 方法。

所有测试均在此笔记本中进行：

• Python粒子模拟器：numpy核外处理

pytables 数据结构简介

参考：PyTables 官方文档

Pytables 允许以两种格式将数据存储在 HDF5 文件中：数组和表。

数组

有 3 种类型的数组Array、CArray和EArray。它们都允许使用类似于 numpy 切片的表示法来存储和检索（多维）切片。

# Write data to store (broadcasting works)
array1[:]  = 3

# Read data from store
in_ram_array = array1[:]

为了在某​​些用例中进行优化，保存在“块”中，可以在创建时CArray选择其大小。chunk_shape

Array并且CArray大小在创建时是固定的。不过，您可以在创建后逐块填充/写入数组。反之EArray可以用该.append()方法进行扩展。

表格

这table是一个完全不同的野兽。它基本上是一张“桌子”。您只有一维索引，每个元素都是一行。每行内部都有“列”数据类型，每列可以有不同的类型。如果您熟悉numpy record-arrays，那么表基本上是一个一维记录数组，每个元素都有许多字段作为列。

一维或二维 numpy 数组可以存储在表中，但这有点棘手：我们需要创建一个行数据类型。例如，要存储 1D uint8 numpy 数组，我们需要执行以下操作：

table_uint8 = np.dtype([('field1', 'u1')])
table_1d = data_file.create_table('/', 'array_1d', description=table_uint8)

那么为什么要使用表格呢？因为与数组不同，表可以高效地查询。例如，如果我们想在一个巨大的基于磁盘的表中搜索> 3的元素，我们可以这样做：

index = table_1d.get_where_list('field1 > 3')

它不仅简单（与我们需要分块扫描整个文件并index在循环中构建的数组相比），而且速度也非常快。

如何存储仿真参数

存储模拟参数的最佳方法是使用组（即/parameters），将每个标量转换为 numpy 数组并将其存储为CArray.

数组“ emission”

emission是按顺序生成和读取的最大数组。对于这种使用模式，一个好的数据结构是EArray。在具有约 50% 零元素的“模拟”数据上，blosc压缩 ( level=5) 实现了 2.2 倍的压缩比。我发现 2^18 (256k) 的块大小具有最短的处理时间。

正在存储“ counts”

同时存储“ counts”将使文件大小增加 10%，并且计算时间戳的时间将增加 40%。存储counts本身并不是一个优势，因为最终只需要时间戳。

优点是重建索引（时间戳）更简单，因为我们在单个命令（）中查询完整时间轴.get_where_list('counts >= 1')。相反，通过分块处理，我们需要执行一些索引算术，这有点棘手，并且可能是维护的负担。

然而，与这两种情况下所需的所有其他操作（排序和合并）相比，代码复杂性可能很小。

正在存储“ timestamps”

时间戳可以在 RAM 中累积。然而，我们在开始之前不知道数组的大小，并且hstack()需要最终调用来“合并”存储在列表中的不同块。这会使内存需求加倍，因此 RAM 可能会不足。

我们可以使用将即时时间戳存储到表中.append()。最后我们可以使用 将该表加载到内存中.read()。这仅比全内存计算慢 10%，但避免了“双 RAM”要求。此外，我们可以避免最终的满载并具有最小的 RAM 使用量。

吡咯

H5py是一个比pytables简单得多的库。对于（主要）顺序处理的用例来说，似乎比 pytables 更适合。唯一缺少的功能是缺乏“blosc”压缩。这是否会导致很大的性能损失还有待测试。