我有3个Numpy阵列,需要在它们之间形成笛卡尔产品。阵列的尺寸不是固定的,因此它们可以采用不同的值,一个示例可以是a =(10000,50(,b =(40,50(,c =(10000,50(。
然后,我在下面执行一些处理(例如A B-C(是我用于产品的功能。
def cartesian_2d(arrays, out=None):
arrays = [np.asarray(x) for x in arrays]
dtype = arrays[0].dtype
n = np.prod([x.shape[0] for x in arrays])
if out is None:
out = np.empty([n, len(arrays), arrays[0].shape[1]], dtype=dtype)
m = n // arrays[0].shape[0]
out[:, 0] = np.repeat(arrays[0], m, axis=0)
if arrays[1:]:
cartesian_2d(arrays[1:], out=out[0:m, 1:, :])
for j in range(1, arrays[0].shape[0]):
out[j * m:(j + 1) * m, 1:] = out[0:m, 1:]
return out
a = [[ 0, -0.02], [1, -0.15]]
b = [[0, 0.03]]
result = cartesian_2d([a,b,a])
// array([[[ 0. , -0.02],
[ 0. , 0.03],
[ 0. , -0.02]],
[[ 0. , -0.02],
[ 0. , 0.03],
[ 1. , -0.15]],
[[ 1. , -0.15],
[ 0. , 0.03],
[ 0. , -0.02]],
[[ 1. , -0.15],
[ 0. , 0.03],
[ 1. , -0.15]]])
输出与itertools.product
相同。但是,我正在使用我的自定义功能来利用Numpy Vectorized操作,与Itertool相比,它可以正常工作。
之后,我做
result[:, 0, :] + result[:, 1, :] - result[:, 2, :]
//array([[ 0. , 0.03],
[-1. , 0.16],
[ 1. , -0.1 ],
[ 0. , 0.03]])
这是最终的预期结果。
只要我的数组适合内存,该功能就可以按预期工作。但是我的用户酶要求我使用大量数据,并且由于无法分配所需的内存,因此我在np.empty()
行中获得了MemoryError。目前,我正在使用大约20GB数据,这可能会在将来增加。
这些数组代表向量,必须存储在float
中,因此我不能使用int
。另外,它们是密集的数组,因此使用sparse
不是一个选项。
我将使用这些数组进行进一步处理,理想情况下,我不想将它们存储在此阶段。因此memmap
/h5py
格式可能无济于事,尽管我不确定这一点。
如果还有其他形成该产品的方法,那也可以。
我敢肯定,有比这更大的数据集的应用程序,我希望有人以前遇到过此类问题,并想知道如何解决这个问题。请帮助。
如果至少您的结果适合内存
以下会产生您的预期结果,而无需依赖中间体的三倍。它使用广播。
请注意,几乎所有的numpy操作都是可以播放的,因此实际上,可能不需要明确的笛卡尔产品:
#shared dimensions:
sh = a.shape[1:]
aba = (a[:, None, None] + b[None, :, None] - a[None, None, :]).reshape(-1, *sh)
aba
#array([[ 0. , 0.03],
# [-1. , 0.16],
# [ 1. , -0.1 ],
# [ 0. , 0.03]])
通过'id'
解决结果行您可以考虑遗漏reshape
。这将使您可以通过组合索引在结果中解决这些行。如果您的组件ID仅为0,1,2,那么...就像您的示例中一样,这将与合并ID相同。例如,ABA [1,0,0]将对应于A。
一些解释
广播:例如,当添加两个阵列时,它们的形状并不一定是相同的,仅由于广播而兼容。从某种意义上说,广播是将标量添加到数组的概括:
[[2], [[7], [[2],
7 + [3], equiv to [7], + [3],
[4]] [7]] [4]]
广播:
[[4], [[1, 2, 3], [[4, 4, 4],
[[1, 2, 3]] + [5], equiv to [1, 2, 3], + [5, 5, 5],
[6]] [1, 2, 3]] [6, 6, 6]]
为此,每个操作数的每个维度都必须是1或等于彼此操作数的相应维度,除非是1。如果操作数的尺寸少于其他尺寸,则其形状在<em上的形状填充。>左。请注意,图中所示的Equiv数组未明确创建。
如果结果也不适合
在那种情况下,我看不出如何避免使用存储,所以h5py
或类似的东西。
从每个操作数中删除第一列
这只是切片的问题:
a_no_id = a[:, 1:]
等。请注意,与Python列表不同,切成薄片时,Numpy阵列不会返回副本,而是视图。因此,效率(内存或运行时(不是这里的问题。
替代解决方案是创建指数的笛卡尔产品(这更容易,因为解决方案更容易对于存在1D数组的笛卡尔产品(:
idx = cartesian_product(
np.arange(len(a)),
np.arange(len(b)) + len(a),
np.arange(len(a))
)
然后使用花哨的索引创建输出数组:
x = np.concatenate((a, b))
result = x[idx.ravel(), :].reshape(*idx.shape, -1)
在磁盘上有效地写作结果
起初关于结果数据大小的几个思维。
结果数据的大小
size_in_GB = A.shape[0]**2*A.shape[1]*B.shape[0]*(size_of_datatype)/1e9
在您提到的问题中,A.Shape =(10000,50(,B =(40,50(。使用float64,您的结果将为1600 GB。如果您有足够的磁盘空间,可以在没有问题的情况下完成此操作,但是您必须考虑接下来与数据无关。也许这只是一个中间结果,可以在块中处理数据。
如果不是这种情况,这是一个示例,如何有效地处理1600GB的数据(RAM使用率约为200 MB(。在现实数据上,水槽应约为200 mb/s。
计算结果的代码来自@paulpanzer。
import numpy as np
import tables #register blosc
import h5py as h5
import h5py_cache as h5c
a=np.arange(500*50).reshape(500, 50)
b=np.arange(40*50).reshape(40, 50)
# isn't well documented, have a look at https://github.com/Blosc/hdf5-blosc
compression_opts=(0, 0, 0, 0, 5, 1, 1)
compression_opts[4]=9 #compression level 0...9
compression_opts[5]=1 #shuffle
compression_opts[6]=1 #compressor (I guess that's lz4)
File_Name_HDF5='Test.h5'
f = h5.File(File_Name_HDF5, 'w',chunk_cache_mem_size=1024**2*300)
dset = f.create_dataset('Data', shape=(a.shape[0]**2*b.shape[0],a.shape[1]),dtype='d',chunks=(a.shape[0]*b.shape[0],1),compression=32001,compression_opts=(0, 0, 0, 0, 9, 1, 1), shuffle=False)
#Write the data
for i in range(a.shape[0]):
sh = a.shape[1:]
aba = (a[i] + b[:, None] - a).reshape(-1, *sh)
dset[i*a.shape[0]*b.shape[0]:(i+1)*a.shape[0]*b.shape[0]]=aba
f.close()
阅读数据
File_Name_HDF5='Test.h5'
f = h5c.File(File_Name_HDF5, 'r',chunk_cache_mem_size=1024**2*300)
dset=f['Data']
chunks_size=500
for i in range(0,dset.shape[0],chunks_size):
#Iterate over the first column
data=dset[i:i+chunks_size,:] #avoid excessive calls to the hdf5 library
#Do something with the data
f.close()
f = h5c.File(File_Name_HDF5, 'r',chunk_cache_mem_size=1024**2*300)
dset=f['Data']
for i in range(dset.shape[1]):
# Iterate over the second dimension
# fancy indexing e.g.[:,i] will be much slower
# use np.expand_dims or in this case np.squeeze after the read operation from the dset
# if you wan't to have the same result than [:,i] (1 dim array)
data=dset[:,i:i+1]
#Do something with the data
f.close()
在此测试示例中,我得到了约550 m/s的写入吞吐量,大约(500 m/s的第一个昏暗,1000m/s秒的读取(和50的压缩比。可接受的速度如果您沿最快的变化方向读取或写入数据(在C的最后一个维度中(,则使用HDF5使用的块数据格式,这根本不是问题。Numpy Memmap也无法进行压缩,导致更高的文件尺寸和较慢的速度。
请注意,必须设置压缩过滤器和块状形状。这取决于您之后如何阅读数据和实际数据。
如果您做了完全错误的事情,则与正确的方法相比,perfornance可以慢10-100倍(例如,可以针对第一个或第二个读取示例优化块状(。