Является ли двумерный numpy.take быстрым?

numpy.take можно применять в 2 измерениях с помощью

np.take(np.take(T,ix,axis=0), iy,axis=1 )

Я протестировал шаблон дискретного двумерного лапласиана

ΔT = T[ix-1,iy] + T[ix+1, iy] + T[ix,iy-1] + T[ix,iy+1] - 4 * T[ix,iy]

с 2 схемами взятия и обычной схемой numpy.array. Функции p и q введены для более компактного написания кода и обращаются к оси 0 и 1 в другом порядке. Это код:

nx = 300; ny= 300
T  = np.arange(nx*ny).reshape(nx, ny)
ix = np.linspace(1,nx-2,nx-2,dtype=int) 
iy = np.linspace(1,ny-2,ny-2,dtype=int)
#------------------------------------------------------------
def p(Φ,kx,ky):
    return np.take(np.take(Φ,ky,axis=1), kx,axis=0 )
#------------------------------------------------------------
def q(Φ,kx,ky):
    return np.take(np.take(Φ,kx,axis=0), ky,axis=1 )
#------------------------------------------------------------
%timeit ΔT_n = T[0:nx-2,1:ny-1] + T[2:nx,1:ny-1] + T[1:nx-1,0:ny-2]  + T[1:nx-1,2:ny] - 4.0 * T[1:nx-1,1:ny-1] 
%timeit ΔT_t = p(T,ix-1,iy)  + p(T,ix+1,iy)  + p(T,ix,iy-1)  + p(T,ix,iy+1)  - 4.0 * p(T,ix,iy)
%timeit ΔT_t = q(T,ix-1,iy)  + q(T,ix+1,iy)  + q(T,ix,iy-1)  + q(T,ix,iy+1)  - 4.0 * q(T,ix,iy)
.
1000 loops, best of 3: 944 µs per loop
100 loops, best of 3: 3.11 ms per loop
100 loops, best of 3: 2.02 ms per loop

Результаты кажутся очевидными:

  1. обычная арифметика индекса numpy быстрее всего
  2. схема взятия q занимает на 100 % больше времени (= C-заказ?)
  3. схема взятия p занимает на 200% больше времени (= Fortran-упорядочивание?)

Даже одномерный пример руководства по scipy указывает, что numpy.take работает быстро:

a = np.array([4, 3, 5, 7, 6, 8])
indices = [0, 1, 4]
%timeit np.take(a, indices)
%timeit a[indices]
.
The slowest run took 6.58 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 4.32 µs per loop
The slowest run took 7.34 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 3.87 µs per loop

У кого-нибудь есть опыт, как сделать numpy.take быстро? Это был бы гибкий и привлекательный способ написания бережливого кода, быстрый в написании кода и
также считается быстрым в исполнении. Спасибо за некоторые подсказки, чтобы улучшить мой подход!


python numpy performance take
person pyano    schedule 24.07.2017    source источник
comment
Как насчет np.ix_ : T[np.ix_(ix,iy)]?   -  person Divakar    schedule 25.07.2017
comment
Насколько я помню из прошлых тестов, np.take немного быстрее, чем нотация индексации. Но преимущество достаточно маленькое, чтобы обернуть его в вызов функции, как вы это делаете, может разрушить его. stackoverflow.com/ вопросы/44487889/   -  person hpaulj    schedule 25.07.2017
comment
@Divakar: да, я тоже пробовал np.ix_ (но пропустил его из-за краткости в моем вопросе): в моих тестах np.ix_ был медленнее, чем лучший np.take   -  person pyano    schedule 25.07.2017
comment
Будут ли ix и iy всегда следовать такому шаблону постоянного размера шага в своих индексах?   -  person Divakar    schedule 25.07.2017
comment
@hpailj: вы правы: я тоже должен попробовать без переноса функций. Но я хотел бы написать довольно сложный код CFD (вычисление гидродинамики). Таким образом, бережливое письмо имеет важное значение, соответственно. небережливый код очень подвержен ошибкам.   -  person pyano    schedule 25.07.2017
comment
T[0:nx-2,1:ny-1] — это представление (базовая индексация с 2 срезами). p(T,ix-1,iy) — это копия, индексация с двумя массивами. Суммы будут, конечно, копией. take может быть быстрее, чем прямое "причудливое" индексирование, T[ix-1,:][:,iy]   -  person hpaulj    schedule 25.07.2017
comment
@Divakar: я могу попробовать np.ix_, добавив в мой код выше следующее: from numpy import ix_ as q; ΔT_ix_ = T[q(ix-1,iy)] + T[q(ix+1,iy)] + T[q(ix,iy-1)] + T[q(ix,iy+1)] - 4.0 * T[q(ix,iy)]   -  person pyano    schedule 25.07.2017
comment
@hpaulj: прямая причудливая индексация не работает в двух измерениях - я пробовал. Так что T[ix-1,iy] не дает желаемого результата.   -  person pyano    schedule 25.07.2017
comment
@Divakar: извините, я неправильно истолковал ваш вопрос: да, в CFD (вычислительная гидродинамика) со структурированными сетками модели индексов ix и iy обычно не слишком дикие. Но есть большие массивы.   -  person pyano    schedule 25.07.2017
comment
Обратите внимание, что я использовал двойную индексацию. В качестве альтернативы T[(ix-1)[:,None], iy] или пусть ix_ сделает это за вас.   -  person hpaulj    schedule 25.07.2017
comment
@pyano Думаю, в этом случае вы можете использовать slicing : T[ix_start:ix_end:ix_stepsize,..]`.   -  person Divakar    schedule 25.07.2017

Ответы (2)


arrow_upward
2
arrow_downward

Индексированную версию можно очистить с помощью таких объектов среза:

T[0:nx-2,1:ny-1] + T[2:nx,1:ny-1] + T[1:nx-1,0:ny-2]  + T[1:nx-1,2:ny] - 4.0 * T[1:nx-1,1:ny-1]

sy1 = slice(1,ny-1)
sx1 = slice(1,nx-1)
sy2 = slice(2,ny)
sy_2 = slice(0,ny-2)
T[0:nx-2,sy1] + T[2:nx,sy1] + T[sx1,xy_2]  + T[sx1,sy2] - 4.0 * T[sx1,sy1]
person hpaulj    schedule 24.07.2017

arrow_upward
1
arrow_downward

Спасибо @Divakar и @hpaulj! Да, работать с slice тоже можно. Сравнение всех 4 подходов дает:

  1. самый быстрый ex aequo: t(usual np) и t(slice)
  2. t(take) = 2 * t(slice)
  3. t(ix_) = 3 * t(slice)

Вот код и результаты:

import numpy as np
from numpy import ix_ as r
nx = 500;    ny = 500
T = np.arange(nx*ny).reshape(nx, ny)

ix = np.arange(1,nx-1); 
iy = np.arange(1,ny-1);

jx = slice(1,nx-1); jxm = slice(0,nx-2); jxp = slice(2,nx)
jy = slice(1,ny-1); jym = slice(0,ny-2); jyp = slice(2,ny)

#------------------------------------------------------------
def p(U,kx,ky):
    return np.take(np.take(U,kx, axis=0), ky,axis=1)
#------------------------------------------------------------

%timeit ΔT_slice= -T[jxm,jy]     + T[jxp,jy]     - T[jx,jym]     + T[jx,jyp]     - 0.0 * T[jx,jy]
%timeit ΔT_npy  = -T[0:nx-2,1:ny-1] + T[2:nx,1:ny-1] - T[1:nx-1,0:ny-2]  + T[1:nx-1,2:ny] - 0.0 * T[1:nx-1,1:ny-1]
%timeit ΔT_take = -p(T,ix-1,iy)  + p(T,ix+1,iy)  - p(T,ix,iy-1)  + p(T,ix,iy+1)  - 0.0 * p(T,ix,iy)
%timeit ΔT_ix_  = -T[r(ix-1,iy)] + T[r(ix+1,iy)] - T[r(ix,iy-1)] + T[r(ix,iy+1)] - 0.0 * T[r(ix,iy)]
.
100 loops, best of 3: 3.14 ms per loop
100 loops, best of 3: 3.13 ms per loop
100 loops, best of 3: 7.03 ms per loop
100 loops, best of 3: 9.58 ms per loop

Относительно обсуждения просмотра и копирования может быть поучительным следующее:

print("if False --> a view ;   if True --> a copy"  )
print("_slice_ :", T[jx,jy].base is None)
print("_npy_   :", T[1:nx-1,1:ny-1].base is None)
print("_take_  :", p(T,ix,iy).base is None)
print("_ix_    :", T[r(ix,iy)].base is None)
.
if False --> a view ;   if True --> a copy
_slice_ : False
_npy_   : False
_take_  : True
_ix_    : True
person pyano    schedule 25.07.2017