以不同的方式切割不同的numpy阵列行

Question

我正在研究蒙特卡洛辐射传输代码，它通过介质模拟射击光子并统计建模随机游走。它一次缓慢地发射一个光子，所以我想对它进行矢量化并一次运行1000个光子。

我已将光子通过的平板分成nlayers切片，其光学深度为0和depth。实际上，这意味着我有nlayers + 2区域（nlayers加上板块上方的区域和板块下方的区域）。在每一步中，我都必须跟踪每个光子通过哪些层。

我们假设我已经知道两个光子从0层开始。一个步骤结束于第2层，另一个步骤结束于第6层。这由一个数组pastpresent表示看起来像这样：

[[ 0 2] 
[ 0 6]] 

我想与(nlayers + 2)列和2行，以产生一个阵列traveled_through，描述光子i是否通过层j传递（端点包）。它看起来像这样（与nlayers = 10）：

[[ 1 1 1 0 0 0 0 0 0 0 0 0] 
[ 1 1 1 1 1 1 1 0 0 0 0 0]] 

我可以遍历光子产生的traveled_through每个单独的行做到这一点，但这是相当缓慢的，和那种失败运行多个光子点一次，所以我宁愿不这样做。

我试图定义数组如下：

traveled_through = np.zeros((2, nlayers)).astype(int) 
traveled_through[ : , np.min(pastpresent, axis = 1) : np.max(pastpresent, axis = 1) + ] = 1 

的想法是，在给定的光子的行，从开始层直到并且包括结束层索引将被设置为1，其中所有其他剩余0。但是，我得到以下错误：

traveled_through[ : , np.min(pastpresent, axis = 1) : np.max(pastpresent, axis = 1) + 1 ] = 1 
IndexError: invalid slice 

我最好的猜测是，numpy的不允许阵列的不同行以不同的方式使用此方法进行索引。有没有人有如何为任意数量的光子和任意数量的图层生成traveled_through的建议？

Answer 1

如果两个光子始终从0开始，那么可以按如下构建阵列。

首先设置变量...

>>> pastpresent = np.array([[0, 2], [0, 6]]) 
>>> nlayers = 10 

...然后构建数组：

>>> (pastpresent[:,1][:,np.newaxis] + 1 > np.arange(nlayers+2)).astype(int) 
array([[1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0]]) 

或者如果光子具有的任意起始点层：

>>> pastpresent2 = np.array([[1, 7], [3, 9]]) 
>>> (pastpresent2[:,0][:,np.newaxis] < np.arange(nlayers+2)) & 
    (pastpresent2[:,1][:,np.newaxis] + 1 > np.arange(nlayers+2)).astype(int)  
array([[0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0], 
     [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0]]) 

Answer 2

一个小技巧我的那种像这种事情涉及logical_xor ufunc的accumulate方法：

>>> a = np.zeros(10, dtype=int) 
>>> b = [3, 7] 
>>> a[b] = 1 
>>> a 
array([0, 0, 0, 1, 0, 0, 0, 1, 0, 0]) 
>>> np.logical_xor.accumulate(a, out=a) 
array([0, 0, 0, 1, 1, 1, 1, 0, 0, 0]) 

请注意，此设置为1在b位置之间的条目，第一个索引包括最后一个索引，所以你必须处理1个错误，这取决于你究竟在做什么。

随着几行，你可以把它作为工作：

>>> a = np.zeros((3, 10), dtype=int) 
>>> b = np.array([[1, 7], [0, 4], [3, 8]]) 
>>> b[:, 1] += 1 # handle the off by 1 error 
>>> a[np.arange(len(b))[:, None], b] = 1 
>>> a 
array([[0, 1, 0, 0, 0, 0, 0, 0, 1, 0], 
     [1, 0, 0, 0, 0, 1, 0, 0, 0, 0], 
     [0, 0, 0, 1, 0, 0, 0, 0, 0, 1]]) 
>>> np.logical_xor.accumulate(a, axis=1, out=a) 
array([[0, 1, 1, 1, 1, 1, 1, 1, 0, 0], 
     [1, 1, 1, 1, 1, 0, 0, 0, 0, 0], 
     [0, 0, 0, 1, 1, 1, 1, 1, 1, 0]])