Python粒子模拟器:内核外处理
问题描述
在python/numpy中编写蒙特卡罗粒子模拟器(布朗运动和光子发射)。我需要将模拟输出(>>10GB)保存到一个文件中，并在第二步中处理数据。与Windows和Linux的兼容性都很重要。
粒子数(n_particles)为10-100。时间步数(time_size)约为10^9。
该模拟有3个步骤(以下代码用于all-in-RAM版本)：
- 模拟(并存储)
emission速率数组(包含许多几乎为-0的元素)：
- shape (`n_particles` x `time_size`), float32, size **80GB**
泊松计算counts数组，(来自泊松过程的随机值和之前计算的速率)：
- shape (`n_particles` x `time_size`), uint8, size **20GB**计数= np.random.poisson(lam=emission).astype(np.uint8)
- 查找计数的时间戳(或索引)。Counts几乎总是0，因此时间戳数组将适合粒子的时间戳= RAM.Loop (C) for c in counts
我做了一次步骤1，然后重复步骤2-3多次(~100)次。将来，我可能需要在计算counts.之前对emission进行预处理(应用cumsum或其他函数
问题
我有一个在内存中工作的实现，我正在尝试理解实现一个可以扩展到(更)更长模拟的内核外版本的最佳方法。
我希望它存在的东西
我需要将数组保存到一个文件中，并且我想使用单个文件进行模拟。我还需要一种“简单”的方式来存储和调用模拟参数(标量)的字典。
理想情况下，我想要一个文件支持的numpy数组，我可以预先分配和填充块。然后，我想要numpy数组方法(max，cumsum，...)为了透明地工作，只需要一个chunksize关键字来指定在每次迭代中加载多少数组。
更好的是，我想要一个Numexpr，它不是在缓存和内存之间运行，而是在内存和硬盘驱动器之间运行。
实用的选项是什么？
作为第一个选择，我开始尝试pyTables，但我对它的复杂性和抽象性不满意(与numpy如此不同)。此外，我目前的解决方案(阅读下面)是丑陋的，不是很有效。
因此，我寻求答案的选择是
- 实现一个具有所需功能的numpy数组(如何实现？)
- 以更智能的方式使用pytable (不同于另一个库: h5py，blaze，pandas...(到目前为止，我还没有尝试过任何一个)。
暂定解决方案(pyTables)
我将模拟参数保存在'/parameters'组中:每个参数都被转换为一个数值数组标量。详细的解决方案，但它是有效的。
我将emission保存为可扩展数组(EArray)，因为我是以块为单位生成数据的，并且需要追加每个新块(虽然我知道最终的大小)。保存counts更成问题。如果像pytable数组一样保存它，就很难执行像"counts >= 2“这样的查询。因此，我将计数保存为多个表(每个粒子一个)丑陋，并使用.get_where_list('counts >= 2')进行查询。我不确定这是否节省空间，而且生成所有这些表而不是使用单个数组，会显著地占用HDF5文件。此外，奇怪的是，创建这些表需要创建一个自定义数据类型(即使是标准的numpy数据类型)：
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列表中，以便于迭代。
EDIT：这个问题的结果是一个名为PyBroMo的布朗运动模拟器。
回答 5
Stack Overflow用户
发布于 2014-01-09 08:25:48
PyTable解决方案
由于不需要Pandas提供的功能，并且处理速度要慢得多(请参阅下面的notebook )，因此最好的方法是直接使用PyTables或h5py。到目前为止，我只尝试了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大小在创建时是固定的。不过，您可以在创建后逐个块填充/写入数组。相反，可以使用.append()方法扩展EArray。
表格
table是一个完全不同的野兽。它基本上就是一个“表”。您只有一维索引，并且每个元素都是一行。在每一行中有“列”数据类型，每列可以有不同的类型。如果您熟悉numpy record-arrays，表基本上是一个一维记录数组，每个元素都有许多字段作为列。
一维或二维numpy数组可以存储在表中，但这有点棘手:我们需要创建一个row数据类型。例如，要存储一个1DNumpy数组，我们需要这样做: uint8：
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
H5py是一个比pytables简单得多的库。对于这种(主要)顺序处理的用例，似乎比pytable更适合。唯一缺少的功能是缺少“blosc”压缩。这是否会导致较大的性能损失还有待测试。
Stack Overflow用户
发布于 2015-05-02 03:16:01
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，该DAG将根据需要使用来自磁盘的多核数据流高效地执行
da.exp(x).mean(axis=0).compute()http://dask.pydata.org/en/latest/
conda install dask
or
pip install daskStack Overflow用户
发布于 2014-01-06 20:53:40
有关如何在HDF5文件中存储参数的信息，请参阅here (它会对参数进行存储，因此您可以按自己的方式存储参数；它们的大小限制为64kb )。
import pandas as pd
import numpy as np
n_particles = 10
chunk_size = 1000
# 1) create a new emission file, compressing as we go
emission = pd.HDFStore('emission.hdf',mode='w',complib='blosc')
# generate simulated data
for i in range(10):
df = pd.DataFrame(np.abs(np.random.randn(chunk_size,n_particles)),dtype='float32')
# create a globally unique index (time)
# http://stackoverflow.com/questions/16997048/how-does-one-append-large-amounts-of-
data-to-a-pandas-hdfstore-and-get-a-natural/16999397#16999397
try:
nrows = emission.get_storer('df').nrows
except:
nrows = 0
df.index = pd.Series(df.index) + nrows
emission.append('df',df)
emission.close()
# 2) create counts
cs = pd.HDFStore('counts.hdf',mode='w',complib='blosc')
# this is an iterator, can be any size
for df in pd.read_hdf('emission.hdf','df',chunksize=200):
counts = pd.DataFrame(np.random.poisson(lam=df).astype(np.uint8))
# set the index as the same
counts.index = df.index
# store the sum across all particles (as most are zero this will be a
# nice sub-selector
# better maybe to have multiple of these sums that divide the particle space
# you don't have to do this but prob more efficient
# you can do this in another file if you want/need
counts['particles_0_4'] = counts.iloc[:,0:4].sum(1)
counts['particles_5_9'] = counts.iloc[:,5:9].sum(1)
# make the non_zero column indexable
cs.append('df',counts,data_columns=['particles_0_4','particles_5_9'])
cs.close()
# 3) find interesting counts
print pd.read_hdf('counts.hdf','df',where='particles_0_4>0')
print pd.read_hdf('counts.hdf','df',where='particles_5_9>0') 或者，也可以使每个粒子成为data_column，然后分别对其进行选择。
和一些输出(在本例中是相当活跃的发射:)
0 1 2 3 4 5 6 7 8 9 particles_0_4 particles_5_9
0 2 2 2 3 2 1 0 2 1 0 9 4
1 1 0 0 0 1 0 1 0 3 0 1 4
2 0 2 0 0 2 0 0 1 2 0 2 3
3 0 0 0 1 1 0 0 2 0 3 1 2
4 3 1 0 2 1 0 0 1 0 0 6 1
5 1 0 0 1 0 0 0 3 0 0 2 3
6 0 0 0 1 1 0 1 0 0 0 1 1
7 0 2 0 2 0 0 0 0 2 0 4 2
8 0 0 0 1 3 0 0 0 0 1 1 0
10 1 0 0 0 0 0 0 0 0 1 1 0
11 0 0 1 1 0 2 0 1 2 1 2 5
12 0 2 2 4 0 0 1 1 0 1 8 2
13 0 2 1 0 0 0 0 1 1 0 3 2
14 1 0 0 0 0 3 0 0 0 0 1 3
15 0 0 0 1 1 0 0 0 0 0 1 0
16 0 0 0 4 3 0 3 0 1 0 4 4
17 0 2 2 3 0 0 2 2 0 2 7 4
18 0 1 2 1 0 0 3 2 1 2 4 6
19 1 1 0 0 0 0 1 2 1 1 2 4
20 0 0 2 1 2 2 1 0 0 1 3 3
22 0 1 2 2 0 0 0 0 1 0 5 1
23 0 2 4 1 0 1 2 0 0 2 7 3
24 1 1 1 0 1 0 0 1 2 0 3 3
26 1 3 0 4 1 0 0 0 2 1 8 2
27 0 1 1 4 0 1 2 0 0 0 6 3
28 0 1 0 0 0 0 0 0 0 0 1 0
29 0 2 0 0 1 0 1 0 0 0 2 1
30 0 1 0 2 1 2 0 2 1 1 3 5
31 0 0 1 1 1 1 1 0 1 1 2 3
32 3 0 2 1 0 0 1 0 1 0 6 2
33 1 3 1 0 4 1 1 0 1 4 5 3
34 1 1 0 0 0 0 0 3 0 1 2 3
35 0 1 0 0 1 1 2 0 1 0 1 4
36 1 0 1 0 1 2 1 2 0 1 2 5
37 0 0 0 1 0 0 0 0 3 0 1 3
38 2 5 0 0 0 3 0 1 0 0 7 4
39 1 0 0 2 1 1 3 0 0 1 3 4
40 0 1 0 0 1 0 0 4 2 2 1 6
41 0 3 3 1 1 2 0 0 2 0 7 4
42 0 1 0 2 0 0 0 0 0 1 3 0
43 0 0 2 0 5 0 3 2 1 1 2 6
44 0 2 0 1 0 0 1 0 0 0 3 1
45 1 0 0 2 0 0 0 1 4 0 3 5
46 0 2 0 0 0 0 0 1 1 0 2 2
48 3 0 0 0 0 1 1 0 0 0 3 2
50 0 1 0 1 0 1 0 0 2 1 2 3
51 0 0 2 0 0 0 2 3 1 1 2 6
52 0 0 2 3 2 3 1 0 1 5 5 5
53 0 0 0 2 1 1 0 0 1 1 2 2
54 0 1 2 2 2 0 1 0 2 0 5 3
55 0 2 1 0 0 0 0 0 3 2 3 3
56 0 1 0 0 0 2 2 0 1 1 1 5
57 0 0 0 1 1 0 0 1 0 0 1 1
58 6 1 2 0 2 2 0 0 0 0 9 2
59 0 1 1 0 0 0 0 0 2 0 2 2
60 2 0 0 0 1 0 0 1 0 1 2 1
61 0 0 3 1 1 2 0 0 1 0 4 3
62 2 0 1 0 0 0 0 1 2 1 3 3
63 2 0 1 0 1 0 1 0 0 0 3 1
65 0 0 1 0 0 0 1 5 0 1 1 6
.. .. .. .. .. .. .. .. .. .. ... ...
[9269 rows x 12 columns]https://stackoverflow.com/questions/20940805
复制相似问题
腾讯云开发者
