8

我有两个 numpy 数组“元素”和“节点”。我的目标是收集这些数组的一些数据。我需要用“节点”数组中包含的两个坐标替换最后两列的“元素”数据。这两个数组非常庞大,我必须自动化它。

这篇文章引用了一篇旧文章:Replace data of an array by 2 values of a second array

不同的是,数组非常大(元素:(3342558,5)和节点:(581589,4)),以前的出路不起作用。

一个例子 :

    import numpy as np
    
    Elements = np.array([[1.,11.,14.],[2.,12.,13.]])
    
    nodes = np.array([[11.,0.,0.],[12.,1.,1.],[13.,2.,2.],[14.,3.,3.]])
    
    results = np.array([[1., 0., 0., 3., 3.],
    [2., 1., 1., 2., 2.]])

hpaulj 之前提出的出路

    e = Elements[:,1:].ravel().astype(int)
    n=nodes[:,0].astype(int)
    
    I, J = np.where(e==n[:,None])
    
    results = np.zeros((e.shape[0],2),nodes.dtype)
    results[J] = nodes[I,:1]
    results = results.reshape(2,4)

但是对于巨大的数组,这个脚本不起作用:
DepreciationWarning: elementwise comparison failed; this will raise an error in the future......

4

4 回答 4

2

Elements大部分游戏都是从in 中找出相应的匹配索引nodes

方法#1

由于您似乎可以转换为整数,因此假设我们可以将它们视为整数。有了这个,我们可以使用基于array-assignment+mapping的方法,如下所示:

ar = Elements.astype(int)
a = ar[:,1:].ravel()
nd = nodes[:,0].astype(int)

n = a.max()+1
# for generalized case of neagtive ints in a or nodes having non-matching values:
# n = max(a.max()-min(0,a.min()), nd.max()-min(0,nd.min()))+1

lookup = np.empty(n, dtype=int)
lookup[nd] = np.arange(len(nd))
indices = lookup[a]

nc = (Elements.shape[1]-1)*(nodes.shape[1]-1) # 4 for given setup
out = np.concatenate((ar[:,0,None], nodes[indices,1:].reshape(-1,nc)),axis=1)

方法#2

我们也可以np.searchsorted用来获取那些indices.

对于根据第一个 col 和匹配大小写对行进行排序的节点,我们可以简单地使用:

indices = np.searchsorted(nd, a)

对于不必要的排序案例和匹配案例:

sidx = nd.argsort()
idx = np.searchsorted(nd, a, sorter=sidx)
indices = sidx[idx]

对于不匹配的情况,使用无效的 bool 数组:

invalid = idx==len(nd)
idx[invalid] = 0
indices = sidx[idx]

方法#3

另一个带concatenation+ sorting-

b = np.concatenate((nd,a))
sidx = b.argsort(kind='stable')

n = len(nd)
v = sidx<n
counts = np.diff(np.flatnonzero(np.r_[v,True]))
r = np.repeat(sidx[v], counts)

indices = np.empty(len(a), dtype=int)
indices[sidx[~v]-n] = r[sidx>=n]

要检测不匹配的,请使用:

nd[indices] != a

将这里的想法移植到numba

from numba import njit

def numba1(Elements, nodes):
    a = Elements[:,1:].ravel()
    nd = nodes[:,0]
    b = np.concatenate((nd,a))
    sidx = b.argsort(kind='stable')
    
    n = len(nodes)        
    ncols = Elements.shape[1]-1
    size = nodes.shape[1]-1        
    dt = np.result_type(Elements.dtype, nodes.dtype)
    nc = ncols*size
    
    out = np.empty((len(Elements),1+nc), dtype=dt)
    out[:,0] = Elements[:,0]
    return numba1_func(out, sidx, nodes, n, ncols, size)

@njit
def numba1_func(out, sidx, nodes, n, ncols, size):
    N = len(sidx)    
    for i in range(N):
        if sidx[i]<n:
            cur_id = sidx[i]
            continue
        else:
            idx = sidx[i]-n        
            row = idx//ncols
            col = idx-row*ncols        
            cc = col*size+1
            for ii in range(size):
                out[row, cc+ii] = nodes[cur_id,ii+1]
    return out
于 2020-10-23T10:06:15.947 回答
1

你会考虑使用pandas吗?

import pandas as pd
Elements = np.array([[1.,11.,14.],[2.,12.,13.]])
nodes = np.array([[11.,0.,0.],[12.,1.,1.],[13.,2.,2.],[14.,3.,3.]])

df_elements = pd.DataFrame(Elements,columns = ['idx','node1','node2'])
df_nodes = pd.DataFrame(nodes, columns = ['node_id','x','y'])

#Double merge to get the coordinates from df_nodes
results = df_elements.merge(df_nodes, left_on = 'node1', right_on="node_id", how='left').merge(df_nodes, left_on="node2",right_on = "node_id", how='left')[['idx',"x_x",'y_x','x_y','y_y']].values

输出

array([[1., 0., 0., 3., 3.],
       [2., 1., 1., 2., 2.]])
于 2020-10-23T10:19:52.627 回答
0

首先,让我们估计一下数组的大小,看看是否会遇到内存错误

from sys import getsizeof

Element_size = getsizeof(np.random.randint(0,100,(3342558,5))) / (1024**3)
nodes_size = getsizeof(np.random.randint(0,100,(581589,4))) / (1024**3)
result_size = getsizeof(np.random.randint(0,100,(3342558,13))) / (1024**3)

total_size = Element_size + nodes_size + result_size

运行13=(5-1)*(4-1)+1这个脚本total_size0.46 GB

我们首先创建要使用的数组

elements = np.random.randint(0,100,(100,5))
elements[:,0] = np.arange(100)
nodes = np.random.randint(0,100,(300,4))

# create an empty result array 
results = np.empty((100,13)).astype(elements.dtype)
results[:,:5] = elements

如您所见,我们results首先创建了数组,一开始创建此数组有两个好处

  1. 大多数操作可以是在results.
  2. 如果内存空间不够,创建的时候就知道了results

使用这些数组,您可以解决您的问题

aux_inds = np.arange(4)
def argmax_with_exception(row):
    
    mask = row[1:5][:,None] == nodes[:,0]
    indices = np.argmax(mask,axis=1)
    node_slices = nodes[indices][:,1:]

    # if a node in Element is not found in the array nodes
    not_found = aux_inds[~np.any(mask,axis=1)]
    node_slices[not_found] = np.ones(3) * -999
    row[1:] = node_slices.flatten()
    
np.apply_along_axis(argmax_with_exception,1,results)

其中,如果在Element中找不到节点,则将nodes其值分配给(-999,-999,-999)

在这种方法中,np.apply_along_axis(argmax_with_exception,1, results)将对数组执行就地操作results,因此,只要可以首先创建数组,就不太可能遇到内存错误。但是,如果您正在使用的机器的 RAM 非常小,您可以Elements首先将阵列保存到磁盘,然后将其加载resultsresults[:,:5] = np.load('Elements.npy')

于 2020-10-24T06:37:54.250 回答
0

为了理解pythonic的解决方案先看sgnfis在老帖子上提供的解决方案: Old solution

import numpy as np
# I used numpy 1.10.1 here

Elements = np.array([[1.,11.,14.],[2.,12.,13.]])
nodes = np.array([[11.,0.,0.],[12.,1.,1.],[13.,2.,2.],[14.,3.,3.]])

# Create an array with enough rows and five columns
res = np.zeros((np.shape(Elements)[0],5))

for i in range(np.shape(Elements)[0]):
    res[i,0] = Elements[i,0] # The first column stays the same

    # Find the Value of the 2nd column of Elements in the first column of nodes.
    nodesindex = np.where(nodes[:,0]==Elements[i,1])
    # Replace second and third row of the results with the ventries from nodes.
    res[i,1:3]=nodes[nodesindex,1:3]

    #Do the same for the 3rd column of Elements
    nodesindex = np.where(nodes[:,0]==Elements[i,2])
    res[i,3:5]=nodes[nodesindex,1:3]

print(res)

上面的解决方案现在变成了 pythonic 解决方案,如下所示: 新解决方案:

import numpy as np

Elements = np.array([[1.,11.,14.],[2.,12.,13.]])
nodes = np.array([[11.,0.,0.],[12.,1.,1.],[13.,2.,2.],[14.,3.,3.]])

# Create an array with enough rows and five columns
res = np.zeros((np.shape(Elements)[0],5))
res[:,0] = Elements[:,0]  # The first column stays the same
res[:,1:3]=[nodes[np.where(nodes[:,0]==Elements[i,1]),1:3] for i in range(np.shape(Elements)[0])]
res[:,3:5]=[nodes[np.where(nodes[:,0]==Elements[i,2]),1:3] for i in range(np.shape(Elements)[0])]
print(res)
于 2020-10-26T02:33:50.607 回答