Cal*_*leb 82 python numpy hdf5 pytables h5py
我正在处理大型3D阵列,我经常需要以各种方式进行切片以进行各种数据分析.典型的"立方体"可以是~100GB(将来可能会变大)
似乎python中大型数据集的典型推荐文件格式是使用HDF5(h5py或pytables).我的问题是:使用HDF5存储和分析这些立方体而不是将它们存储在简单的平面二进制文件中是否有任何速度或内存使用效益？HDF5是否更适合表格数据,而不像我正在使用的大型数组？我看到HDF5可以提供很好的压缩,但我对处理速度和处理内存溢出更感兴趣.
我经常只想分析立方体的一个大的子集.pytables和h5py的一个缺点是,当我拿一个数组时,我总是得到一个numpy数组,用尽内存.但是,如果我切片平面二进制文件的numpy memmap,我可以得到一个视图,它将数据保存在磁盘上.因此,似乎我可以更轻松地分析我的数据的特定部分而不会超出我的记忆.
我已经探讨了pytables和h5py,到目前为止还没有看到我的目的的好处.
Joe*_*ton 142
HDF5的一些主要优点是其层次结构(类似于文件夹/文件),与每个项目一起存储的可选任意元数据,以及其灵活性(例如压缩).这种组织结构和元数据存储可能听起来微不足道,但在实践中它非常有用.
HDF的另一个优点是数据集可以是固定大小或灵活大小.因此,无需创建整个新副本,就可以轻松地将数据附加到大型数据集.
此外,HDF5是一种标准化格式,其库可用于几乎所有语言,因此使用HDF可以轻松地在Matlab,Fortran,R,C和Python之间共享磁盘数据.(公平地说,只要你知道C和F的顺序并知道存储数组的形状,类型等,它对于大二进制数组也不会太难.)
就像TL/DR一样:对于~8GB的3D阵列,沿着任意轴读取"完整"切片需要大约20秒使用分块的HDF5数据集,并且0.3秒(最佳情况)到超过3小时(最差情况)一个相同数据的memmapped数组.
除了上面列出的内容之外,还有一个很大的优势就是"分块"*磁盘数据格式,如HDF5:读取任意切片(强调任意)通常要快得多,因为磁盘上的数据更加连续平均.
*(HDF5不一定是分块数据格式.它支持分块,但不需要它.事实上h5py,如果我没记错的话,创建数据集的默认值不是块.)
基本上,对于给定切片的数据集,最佳情况下的磁盘读取速度和最坏情况下的磁盘读取速度将与分块的HDF数据集相当接近(假设您选择了合理的块大小或让库为您选择一个).用一个简单的二进制数组,最好的情况也比较快,但最坏的情况是多糟糕.
有一点需要注意,如果您有SSD,您可能不会注意到读/写速度的巨大差异.但是,对于常规硬盘驱动器,顺序读取比随机读取快得多.(即常规硬盘驱动器有很seek长时间.)HDF在SSD上仍然具有优势,但更多的是由于其他功能(例如元数据,组织等)而不是原始速度.
首先,为了消除混淆,访问h5py数据集会返回一个与numpy数组非常相似的对象,但在切片之前不会将数据加载到内存中.(与memmap类似,但不完全相同.)查看h5py介绍以获取更多信息.
切片数据集会将数据的一个子集加载到内存中,但可能你想用它做一些事情,此时你无论如何都需要它在内存中.
如果您确实想进行核外计算,可以使用pandas或者相当容易地获取表格数据pytables.有可能h5py(对于大型ND阵列更好),但您需要下拉到较低级别并自己处理迭代.
然而,类似numpy的核外计算的未来是Blaze.如果你真的想采取这条路线,请看看它.
首先,考虑写入磁盘的3D C排序数组(我将通过调用arr.ravel()和打印结果来模拟它,以使事物更加明显):
In [1]: import numpy as np
In [2]: arr = np.arange(4*6*6).reshape(4,6,6)
In [3]: arr
Out[3]:
array([[[ 0, 1, 2, 3, 4, 5],
[ 6, 7, 8, 9, 10, 11],
[ 12, 13, 14, 15, 16, 17],
[ 18, 19, 20, 21, 22, 23],
[ 24, 25, 26, 27, 28, 29],
[ 30, 31, 32, 33, 34, 35]],
[[ 36, 37, 38, 39, 40, 41],
[ 42, 43, 44, 45, 46, 47],
[ 48, 49, 50, 51, 52, 53],
[ 54, 55, 56, 57, 58, 59],
[ 60, 61, 62, 63, 64, 65],
[ 66, 67, 68, 69, 70, 71]],
[[ 72, 73, 74, 75, 76, 77],
[ 78, 79, 80, 81, 82, 83],
[ 84, 85, 86, 87, 88, 89],
[ 90, 91, 92, 93, 94, 95],
[ 96, 97, 98, 99, 100, 101],
[102, 103, 104, 105, 106, 107]],
[[108, 109, 110, 111, 112, 113],
[114, 115, 116, 117, 118, 119],
[120, 121, 122, 123, 124, 125],
[126, 127, 128, 129, 130, 131],
[132, 133, 134, 135, 136, 137],
[138, 139, 140, 141, 142, 143]]])
这些值将按顺序存储在磁盘上,如下面第4行所示.(暂时忽略文件系统细节和碎片.)
In [4]: arr.ravel(order='C')
Out[4]:
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,
78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143])
在最好的情况下,让我们沿第一轴切片.请注意,这些只是数组的前36个值.这将是一个非常快速的阅读!(一个寻求,一个阅读)
In [5]: arr[0,:,:]
Out[5]:
array([[ 0, 1, 2, 3, 4, 5],
[ 6, 7, 8, 9, 10, 11],
[12, 13, 14, 15, 16, 17],
[18, 19, 20, 21, 22, 23],
[24, 25, 26, 27, 28, 29],
[30, 31, 32, 33, 34, 35]])
类似地,沿第一轴的下一个切片将仅是接下来的36个值.要沿此轴读取整个切片,我们只需要一个seek操作.如果我们要读的是沿着这个轴的各种切片,那么这就是完美的文件结构.
但是,让我们考虑最坏情况:沿最后一个轴的切片.
In [6]: arr[:,:,0]
Out[6]:
array([[ 0, 6, 12, 18, 24, 30],
[ 36, 42, 48, 54, 60, 66],
[ 72, 78, 84, 90, 96, 102],
[108, 114, 120, 126, 132, 138]])
要读取此片段,我们需要36次搜索和36次读取,因为所有值都在磁盘上分开.它们都不相邻!
这可能看起来很小,但随着我们越来越大的阵列,操作的数量和规模seek迅速增长.对于以这种方式存储并读入通道的大型(~10Gb)3D阵列memmap,即使使用现代硬件,沿"最差"轴读取整个切片也可能需要数十分钟.同时,沿最佳轴的切片可能需要不到一秒钟.为简单起见,我只是沿着单个轴显示"完整"切片,但是对于任何数据子集的任意切片都会发生完全相同的事情.
顺便提一下,有几种文件格式可以利用这一点,并且基本上在磁盘上存储三个巨大的 3D阵列副本:一个是C顺序,一个是F顺序,一个是中间的两个.(这方面的一个例子是Geoprobe的D3D格式,虽然我不确定它是否记录在任何地方.)谁在乎最终文件大小是4TB,存储是否便宜!关于这一点的疯狂之处在于,因为主要用例是在每个方向上提取单个子切片,所以要进行的读取非常非常快.它运作得很好!
假设我们将3D阵列的2x2x2"块"存储为磁盘上的连续块.换句话说,像:
nx, ny, nz = arr.shape
slices = []
for i in range(0, nx, 2):
for j in range(0, ny, 2):
for k in range(0, nz, 2):
slices.append((slice(i, i+2), slice(j, j+2), slice(k, k+2)))
chunked = np.hstack([arr[chunk].ravel() for chunk in slices])
所以磁盘上的数据看起来像chunked:
array([ 0, 1, 6, 7, 36, 37, 42, 43, 2, 3, 8, 9, 38,
39, 44, 45, 4, 5, 10, 11, 40, 41, 46, 47, 12, 13,
18, 19, 48, 49, 54, 55, 14, 15, 20, 21, 50, 51, 56,
57, 16, 17, 22, 23, 52, 53, 58, 59, 24, 25, 30, 31,
60, 61, 66, 67, 26, 27, 32, 33, 62, 63, 68, 69, 28,
29, 34, 35, 64, 65, 70, 71, 72, 73, 78, 79, 108, 109,
114, 115, 74, 75, 80, 81, 110, 111, 116, 117, 76, 77, 82,
83, 112, 113, 118, 119, 84, 85, 90, 91, 120, 121, 126, 127,
86, 87, 92, 93, 122, 123, 128, 129, 88, 89, 94, 95, 124,
125, 130, 131, 96, 97, 102, 103, 132, 133, 138, 139, 98, 99,
104, 105, 134, 135, 140, 141, 100, 101, 106, 107, 136, 137, 142, 143])
只是为了表明它们是2x2x2块arr,注意这些是前8个值chunked:
In [9]: arr[:2, :2, :2]
Out[9]:
array([[[ 0, 1],
[ 6, 7]],
[[36, 37],
[42, 43]]])
要读入沿轴的任何切片,我们将读取6或9个连续的块(数据量是我们需要的两倍),然后只保留我们想要的部分.这是最糟糕的情况下最多9次寻求,而最多只寻求36次寻找非分块版本.(但对于memmapped数组,最好的情况仍然是6次搜索与1次.)因为顺序读取与搜索相比非常快,这大大减少了将任意子集读入内存所需的时间.再次,对于更大的阵列,这种效果会变得更大.
HDF5更进一步.块不必连续存储,并且它们由B树索引.此外,它们不必在磁盘上具有相同的大小,因此可以对每个块应用压缩.
h5py默认情况下,h5py不会在磁盘上创建分块的HDF文件(pytables相反,我认为这样做).chunks=True但是,如果在创建数据集时指定,则会在磁盘上获得一个分块数组.
作为一个快速,最小的例子:
import numpy as np
import h5py
data = np.random.random((100, 100, 100))
with h5py.File('test.hdf', 'w') as outfile:
dset = outfile.create_dataset('a_descriptive_name', data=data, chunks=True)
dset.attrs['some key'] = 'Did you want some metadata?'
请注意,chunks=True告诉h5py我们为我们自动选择一个块大小.如果您对最常见的用例有更多了解,可以通过指定形状元组来优化块大小/形状(例如(2,2,2),在上面的简单示例中).这使您可以更有效地沿特定轴读取或优化特定大小的读/写.
为了强调这一点,让我们比较分块中的读数,从分块的HDF5数据集和包含相同精确数据的大型(~8GB)Fortran排序的3D阵列.
我已经清除了每次运行之间的所有操作系统缓存,因此我们看到了"冷"性能.
对于每种文件类型,我们将测试沿第一轴的"完整"x切片和沿最后一轴的"完整"z-slize的读数.对于Fortran排序的memmapped数组,"x"切片是最坏的情况,"z"切片是最好的情况.
使用的代码是一个要点(包括创建hdf文件).我无法轻松共享此处使用的数据,但您可以通过相同形状(621, 4991, 2600)和类型)的零数组来模拟它np.uint8.
该chunked_hdf.py如下所示:
import sys
import h5py
def main():
data = read()
if sys.argv[1] == 'x':
x_slice(data)
elif sys.argv[1] == 'z':
z_slice(data)
def read():
f = h5py.File('/tmp/test.hdf5', 'r')
return f['seismic_volume']
def z_slice(data):
return data[:,:,0]
def x_slice(data):
return data[0,:,:]
main()
memmapped_array.py是类似的,但触摸更复杂,以确保切片实际上加载到内存中(默认情况下,memmapped将返回另一个数组,这不是一个苹果对苹果的比较).
import numpy as np
import sys
def main():
data = read()
if sys.argv[1] == 'x':
x_slice(data)
elif sys.argv[1] == 'z':
z_slice(data)
def read():
big_binary_filename = '/data/nankai/data/Volumes/kumdep01_flipY.3dv.vol'
shape = 621, 4991, 2600
header_len = 3072
data = np.memmap(filename=big_binary_filename, mode='r', offset=header_len,
order='F', shape=shape, dtype=np.uint8)
return data
def z_slice(data):
dat = np.empty(data.shape[:2], dtype=data.dtype)
dat[:] = data[:,:,0]
return dat
def x_slice(data):
dat = np.empty(data.shape[1:], dtype=data.dtype)
dat[:] = data[0,:,:]
return dat
main()
我们先来看看HDF性能:
jofer at cornbread in ~
$ sudo ./clear_cache.sh
jofer at cornbread in ~
$ time python chunked_hdf.py z
python chunked_hdf.py z 0.64s user 0.28s system 3% cpu 23.800 total
jofer at cornbread in ~
$ sudo ./clear_cache.sh
jofer at cornbread in ~
$ time python chunked_hdf.py x
python chunked_hdf.py x 0.12s user 0.30s system 1% cpu 21.856 total
"完整"x切片和"完整"z切片花费大约相同的时间量(~20秒).考虑到这是一个8GB阵列,这还不错.大多数时候
如果我们将它与memmapped数组时间进行比较(它是Fortran有序的:"z-slice"是最好的情况,"x-slice"是最坏的情况.):
jofer at cornbread in ~
$ sudo ./clear_cache.sh
jofer at cornbread in ~
$ time python memmapped_array.py z
python memmapped_array.py z 0.07s user 0.04s system 28% cpu 0.385 total
jofer at cornbread in ~
$ sudo ./clear_cache.sh
jofer at cornbread in ~
$ time python memmapped_array.py x
python memmapped_array.py x 2.46s user 37.24s system 0% cpu 3:35:26.85 total
是的,你没有看错.一个切片方向为0.3秒,另一个为3.5 小时.
在"x"方向切片的时间远远长于将整个8GB阵列加载到内存并选择我们想要的切片所需的时间!(同样,这是一个Fortran排序的数组.相反的x/z切片时序就是C有序数组的情况.)
但是,如果我们总是希望沿着最好的方向进行切片,那么磁盘上的大二进制数组非常好.(~0.3秒!)
使用memmapped数组,你会遇到这种I/O差异(或者各向异性可能是一个更好的术语).但是,对于分块的HDF数据集,您可以选择chunksize,以使访问相等或针对特定用例进行优化.它为您提供了更多的灵活性.
希望无论如何,这有助于澄清问题的一部分.与"原始"memmaps相比,HDF5还有许多其他优点,但我没有足够的空间在这里展开所有这些.压缩可以加快一些东西(我使用的数据并没有太大的压缩受益,所以我很少使用它),和OS级缓存经常与HDF5文件播放更漂亮比"原始" memmaps.除此之外,HDF5是一种非常出色的容器格式.它为您提供了很多管理数据的灵活性,并且可以或多或少地使用任何编程语言.
总的来说,尝试一下,看看它是否适用于您的用例.我想你可能会感到惊讶.
| 归档时间： |
|
| 查看次数： |
20598 次 |
| 最近记录： |