Эффективный способ создания множества случайных координат палочки с numpy

Question

В моделировании методом Монте-Карло я создаю множество списков случайных координат палочки (фактически два списка координат для повторения, представляющих два разных типа палок) в форме [[x0,y0,x1,y1]*N]. Используя векторизованные методы numpy, я попытался свести к минимуму время создания. Однако при определенных условиях длина массивов превышает 10 миллионов, а генерация становится узким местом.

Следующий код дает минимальный пример с некоторым тестом значением

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, вместо того, чтобы делать какое-то случайное сравнение чисел для каждой палочки. Тем не менее, сохраняется более длинная часть создания случайных начальных координат и направления плюс расчет конечных координат.

Кто-нибудь видит дальнейшие оптимизации, чтобы ускорить создание массивов?

Answer 1

Как показывают тайминги x1 & y1 Расчеты являются самыми медленными частями кода. В нем мы имеем cosine и sine вычисления, масштабирование с помощью length, а затем округление и преобразование в int32. Теперь один из способов, которым люди используют для повышения производительности NumPy, - это модуль numexpr.

В нашей самой медленной части операциями, которые могут быть вычислены с помощью numexpr, являются, cosine и масштабирование. Таким образом, 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 

Далее, у нас есть два brances для целей тестирования во время выполнения - один с исходным кодом, а другой с предложил 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% сокращение во время выполнения там, неплохо я угадываю!