如何加速 scipy.weave 中的多维数组访问？

Question

我正在 python 中编写我的 c 代码来加速循环：

from scipy import weave
from numpy import *

#1) create the array
a=zeros((200,300,400),int)
for i in range(200):
    for j in range(300):
        for k in range(400):    
            a[i,j,k]=i*300*400+j*400+k
#2) test on c code to access the array
code="""
for(int i=0;i<200;++i){
for(int j=0;j<300;++j){
for(int k=0;k<400;++k){
printf("%ld,",a[i*300*400+j*400+k]);    
}
printf("\\n");
}
printf("\\n\\n");
}
"""
test =weave.inline(code, ['a'])

它工作得很好，但是当数组很大时它仍然很昂贵。 有人建议我使用 a.strides 而不是令人讨厌的“a[i*300*400+j*400+k]” 我无法理解有关 .strides 的文档。

任何想法

提前致谢

Answer 1

您可以将 3 个 for 循环替换为

grid=np.ogrid[0:200,0:300,0:400]
a=grid[0]*300*400+grid[1]*400+grid[2]

以下建议这可能会导致约 68 倍（或更好？见下文）加速：

% python -mtimeit -s"import test" "test.m1()"
100 loops, best of 3: 17.5 msec per loop
% python -mtimeit -s"import test" "test.m2()"
1000 loops, best of 3: 247 usec per loop

test.py:

import numpy as np

n1,n2,n3=20,30,40
def m1():
    a=np.zeros((n1,n2,n3),int)
    for i in range(n1):
        for j in range(n2):
            for k in range(n3):    
                a[i,j,k]=i*300*400+j*400+k
    return a

def m2():    
    grid=np.ogrid[0:n1,0:n2,0:n3]
    b=grid[0]*300*400+grid[1]*400+grid[2]
    return b 

if __name__=='__main__':
    assert(np.all(m1()==m2()))

With n1,n2,n3 = 200,300,400，

python -mtimeit -s"import test" "test.m2()"

在我的机器上花费了 182 毫秒，并且

python -mtimeit -s"import test" "test.m1()"

有尚未完成。

Answer 2

问题是您在 C 代码中将 240 万个数字打印到屏幕上。这当然需要一段时间，因为必须将数字转换为字符串，然后打印到屏幕上。您真的需要将它们全部打印到屏幕上吗？您的最终目标是什么？

为了进行比较，我尝试将另一个数组设置为 a 中的每个元素。这个过程在编织中花费了大约 0.05 秒。我放弃了在 30 秒左右将所有元素打印到屏幕上的计时。

Answer 3

在 C 中没有办法加速访问多维数组。你必须计算数组索引，并且必须取消引用它，这就是最简单的。

Answer 4

我真的希望，您没有使用所有打印语句运行循环，正如贾斯汀已经指出的那样。除此之外：

from scipy import weave
n1, n2, n3 = 200, 300, 400

def m1():
    a = np.zeros((n1,n2,n3), int)
    for i in xrange(n1):
        for j in xrange(n2):
            for k in xrange(n3):
                a[i,j,k] = i*300*400 + j*400 + k
    return a

def m2():    
    grid = np.ogrid[0:n1,0:n2,0:n3]
    b = grid[0]*300*400 + grid[1]*400 + grid[2]
    return b 

def m3():
    a = np.zeros((n1,n2,n3), int)
    code = """
    int rows = Na[0];
    int cols = Na[1];
    int depth = Na[2];
    int val = 0;      
    for (int i=0; i<rows; i++) {
        for (int j=0; j<cols; j++) {
            for (int k=0; k<depth; k++) {
                val = (i*cols + j)*depth + k;
                a[val] = val;
            }
        }
    }"""
    weave.inline(code, ['a'])
    return a

%timeit m1()
%timeit m2()
%timeit m3()
np.all(m1() == m2())
np.all(m2() == m3())

给我：

1 loops, best of 3: 19.6 s per loop
1 loops, best of 3: 248 ms per loop
10 loops, best of 3: 144 ms per loop

这似乎很合理。如果您想进一步加快速度，您可能需要开始使用 GPU，它非常适合此类数字运算。

在这种特殊情况下，您甚至可以这样做：

def m4():
    a = np.zeros((n1,n2,n3), int)
    code = """
    int rows = Na[0];
    int cols = Na[1];
    int depth = Na[2];
    for (int i=0; i<rows*cols*depth; i++) {
        a[i] = i;
    }"""
    weave.inline(code, ['a'])
    return a

但这并没有变得更好，因为 np.zeros() 已经占用了大部分时间：

%timeit np.zeros((n1,n2,n3), int)
10 loops, best of 3: 113 ms per loop