首页
学习
活动
专区
工具
TVP
发布
社区首页 >问答首页 >Python粒子模拟器:内核外处理

Python粒子模拟器:内核外处理
EN

Stack Overflow用户
提问于 2014-01-06 07:55:16
回答 5查看 1.8K关注 0票数 17

问题描述

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

粒子数(n_particles)为10-100。时间步数(time_size)约为10^9。

该模拟有3个步骤(以下代码用于all-in-RAM版本):

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

代码语言:javascript
复制
- shape (`n_particles` x `time_size`), float32, size **80GB**

泊松计算counts数组,(来自泊松过程的随机值和之前计算的速率):

代码语言:javascript
复制
- shape (`n_particles` x `time_size`), uint8, size **20GB**

计数= np.random.poisson(lam=emission).astype(np.uint8)

  1. 查找计数的时间戳(或索引)。Counts几乎总是0,因此时间戳数组将适合粒子的时间戳= RAM.Loop (C) for c in counts

我做了一次步骤1,然后重复步骤2-3多次(~100)次。将来,我可能需要在计算counts.之前对emission进行预处理(应用cumsum或其他函数

问题

我有一个在内存中工作的实现,我正在尝试理解实现一个可以扩展到(更)更长模拟的内核外版本的最佳方法。

我希望它存在的东西

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

理想情况下,我想要一个文件支持的numpy数组,我可以预先分配和填充块。然后,我想要numpy数组方法(maxcumsum,...)为了透明地工作,只需要一个chunksize关键字来指定在每次迭代中加载多少数组。

更好的是,我想要一个Numexpr,它不是在缓存和内存之间运行,而是在内存和硬盘驱动器之间运行。

实用的选项是什么?

作为第一个选择,我开始尝试pyTables,但我对它的复杂性和抽象性不满意(与numpy如此不同)。此外,我目前的解决方案(阅读下面)是丑陋的,不是很有效。

因此,我寻求答案的选择是

  1. 实现一个具有所需功能的numpy数组(如何实现?)
  2. 以更智能的方式使用pytable (不同于另一个库: h5py,blaze,pandas...(到目前为止,我还没有尝试过任何一个)。

暂定解决方案(pyTables)

我将模拟参数保存在'/parameters'组中:每个参数都被转换为一个数值数组标量。详细的解决方案,但它是有效的。

我将emission保存为可扩展数组(EArray),因为我是以块为单位生成数据的,并且需要追加每个新块(虽然我知道最终的大小)。保存counts更成问题。如果像pytable数组一样保存它,就很难执行像"counts >= 2“这样的查询。因此,我将计数保存为多个表(每个粒子一个)丑陋,并使用.get_where_list('counts >= 2')进行查询。我不确定这是否节省空间,而且生成所有这些表而不是使用单个数组,会显著地占用HDF5文件。此外,奇怪的是,创建这些表需要创建一个自定义数据类型(即使是标准的numpy数据类型):

代码语言:javascript
复制
    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的布朗运动模拟器。

EN

回答 5

Stack Overflow用户

回答已采纳

发布于 2014-01-09 08:25:48

PyTable解决方案

由于不需要Pandas提供的功能,并且处理速度要慢得多(请参阅下面的notebook ),因此最好的方法是直接使用PyTables或h5py。到目前为止,我只尝试了pytables方法。

所有测试都在此笔记本中执行:

pytables数据结构简介

参考:Official PyTables Docs

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

数组

有3种类型的阵列ArrayCArrayEArray。它们都允许使用类似于numpy切片的符号来存储和检索(多维)切片。

代码语言:javascript
复制
# Write data to store (broadcasting works)
array1[:]  = 3

# Read data from store
in_ram_array = array1[:]

对于某些用例中的优化,块被保存在“CArray”中,其大小可以在创建时使用chunk_shape选择。

ArrayCArray大小在创建时是固定的。不过,您可以在创建后逐个块填充/写入数组。相反,可以使用.append()方法扩展EArray

表格

table是一个完全不同的野兽。它基本上就是一个“表”。您只有一维索引,并且每个元素都是一行。在每一行中有“列”数据类型,每列可以有不同的类型。如果您熟悉numpy record-arrays,表基本上是一个一维记录数组,每个元素都有许多字段作为列。

一维或二维numpy数组可以存储在表中,但这有点棘手:我们需要创建一个row数据类型。例如,要存储一个1DNumpy数组,我们需要这样做: uint8:

代码语言:javascript
复制
table_uint8 = np.dtype([('field1', 'u1')])
table_1d = data_file.create_table('/', 'array_1d', description=table_uint8)

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

代码语言:javascript
复制
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”压缩。这是否会导致较大的性能损失还有待测试。

票数 2
EN

Stack Overflow用户

发布于 2015-05-02 03:16:01

Dask.array可以在PyTables或h5py等磁盘阵列上执行maxcumsum等分块操作。

代码语言:javascript
复制
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将根据需要使用来自磁盘的多核数据流高效地执行

代码语言:javascript
复制
da.exp(x).mean(axis=0).compute()

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

代码语言:javascript
复制
conda install dask
or 
pip install dask
票数 3
EN

Stack Overflow用户

发布于 2014-01-06 20:53:40

有关如何在HDF5文件中存储参数的信息,请参阅here (它会对参数进行存储,因此您可以按自己的方式存储参数;它们的大小限制为64kb )。

代码语言:javascript
复制
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,然后分别对其进行选择。

和一些输出(在本例中是相当活跃的发射:)

代码语言:javascript
复制
    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]
票数 2
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/20940805

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档