使用 numpy 创建许多随机棒坐标的省时方法

Time-efficient way for creation of many random stick coordinates with numpy

在蒙特卡洛模拟中,我以 [[x0,y0,x1,y1]*N] 的形式创建了许多随机棒坐标列表(实际上每次重复两个坐标列表,代表两种不同的棒类型)。通过使用向量化的 numpy 方法,我尝试最小化创建时间。然而,在某些情况下,数组的长度超过 10mio,生成成为瓶颈。

下面的代码给出了一个带有一些测试值的最小示例

import numpy as np

def create_coordinates_vect(dimensions=[1500,2500], length=50, count=12000000, type1_content=0.001):
    # two arrays with random start coordinates in area of dimensions
    x0 = np.random.randint(dimensions[0], size=count)
    y0 = np.random.randint(dimensions[1], size=count)
    # random direction of each stick
    dirrad = 2 * np.pi * np.random.rand(count)
    # to destinguish between type1 and type2 sticks based on random values
    stick_type = np.random.rand(count)   
    is_type1 = np.zeros_like(stick_type)
    is_type1[stick_type < type1_content] = True
    # calculate end coordinates
    x1 = x0 + np.rint(np.cos(dirrad) * length).astype(np.int32)
    y1 = y0 + np.rint(np.sin(dirrad) * length).astype(np.int32)
    # stack together start and end coordinates
    coordinates = np.vstack((x0, y0, x1, y1)).T.astype(np.int32)
    # split array according to type
    coords_type1 = coordinates[is_type1 == True]
    coords_type2 = coordinates[is_type1 == False]
    return ([coords_type1, coords_type2])

list1, list2 = create_coordinates_vect()

时序分析给出了不同部分的以下结果

=> x0, y0:                       477.3640632629945 ms
=> dirrad, stick_type:           317.4648284911094 ms
=> is_type1:                      27.3699760437172 ms
=> x1, y1:                      1184.7038269042969 ms
=> vstack:                       189.0783309965234 ms
=> coords_type1, coords_type2:   309.9758625035176 ms

我仍然可以通过预先定义 type1 和 type2 棒的数量来节省一些时间,而不是为每个棒做一些随机数比较。但是,创建随机起始坐标和方向以及计算结束坐标的较长部分将保留。

有人看到进一步的优化来加速数组的创建吗?

如时间所示,x1y1 计算是代码中最慢的部分。在其中,我们有 cosinesine 计算,用 length 缩放,然后四舍五入并转换为 int32。现在,人们用来提升 NumPy 性能的方法之一是使用 numexpr 模块。

在我们最慢的部分,可以用 numexpr 计算的操作是 sinecosine 和缩放。因此,代码的 numexpr 修改版本看起来像这样 -

import numexpr as ne

x1 = x0 + np.rint(ne.evaluate("cos(dirrad) * length")).astype(np.int32)
y1 = y0 + np.rint(ne.evaluate("sin(dirrad) * length")).astype(np.int32)

运行时测试 -

让我们考虑第 (1/100) 个形状到原始数组形状。因此,我们有 -

dimensions=[15,25]
length=50
count=120000
type1_content=0.001

代码的初始部分保持不变 -

# two arrays with random start coordinates in area of dimensions
x0 = np.random.randint(dimensions[0], size=count)
y0 = np.random.randint(dimensions[1], size=count)

# random direction of each stick
dirrad = 2 * np.pi * np.random.rand(count)
# to destinguish between type1 and type2 sticks based on random values
stick_type = np.random.rand(count)   
is_type1 = np.zeros_like(stick_type)
is_type1[stick_type < type1_content] = True

接下来,我们有两个分支用于运行时测试 - 一个是原始代码,另一个是建议的基于 numexpr 的方法 -

def org_app(x0,y0,dirrad,length):
    x1 = x0 + np.rint(np.cos(dirrad) * length).astype(np.int32)
    y1 = y0 + np.rint(np.sin(dirrad) * length).astype(np.int32)

def new_app(x0,y0,dirrad,length):
    x1 = x0 + np.rint(ne.evaluate("cos(dirrad) * length")).astype(np.int32)
    y1 = y0 + np.rint(ne.evaluate("sin(dirrad) * length")).astype(np.int32)

最后,运行时测试自己-

In [149]: %timeit org_app(x0,y0,dirrad,length)
10 loops, best of 3: 23.5 ms per loop

In [150]: %timeit new_app(x0,y0,dirrad,length)
100 loops, best of 3: 14.6 ms per loop    

所以,我们正在考虑 40% 运行时间的减少,我想这还不错!