Есть ли способ, которым я могу перебирать массив быстрее, чем цикл for?

Question

Я пишу код, который сравнивает потоки пикселей на астрономической карте с соответствующей областью на другой.Обе карты представляют собой массивы данных.

Для этого мне нужно преобразовать индексы пикселей на первой карте (Av) в их эквивалент по координатам неба, а затем преобразовать эти координаты неба в эквивалент их индексов пикселей.на второй карте (CO).Затем я масштабирую потоки второй карты, чтобы соответствовать значениям первой карты.После этого я должен продолжать обрабатывать данные.

Проблема в том, что при наличии тысяч пикселей на первой карте коду требуется очень много времени, чтобы выполнить то, что он должен делать, что являетсяхлопоты для устранения неполадок.Я понял, что самой медленной вещью в этой части кода является цикл for.

Есть ли способ, которым я могу перебирать массивный массив, работая с индексами и вычисляяданные из каждого пикселя, быстрее, чем цикл for? Есть ли лучший способ сделать это вообще?

В псевдокоде мой код выглядит примерно так:

for pixel i,j in 1st map:
    sky_x1,sky_y1 = pixel_2_skycoord(i,j)
    i2,j2 = skycoord_2_pixel(sky_x1,sky_y1)

    Avmap.append(Avflux[i,j])
    COmap.append(COflux[i2,j2]*scale)

Фактический код:

for i in xrange(0,sAv_y-1):
    for j in xrange(0,sAv_x-1):
        if not np.isnan(Avdata[i,j]):           

            y,x=wcs.utils.skycoord_to_pixel(wcs.utils.pixel_to_skycoord(i,j,wAv,0),wcs=wCO)

            x=x.astype(int)+0 #the zero is because i don't understand the problem with numpy but it fixes it anyway
            y=y.astype(int)+0 #i couldn't get the number from an array with 1 value but adding zero resolves it somehow
            COflux=COdata[x,y]

            ylist.append(Avdata[i,j])
            xlist.append(COflux*(AvArea/COArea))

Answer 1

Виновником здесь являются две петли.У Numpy есть много функций, которые не позволяют использовать циклы for для быстрого компилирования кода.Хитрость заключается в векторизации вашего кода.

Вы можете заглянуть в функцию numpy meshgrid, чтобы преобразовать эти данные в векторизованную форму, чтобы затем вы могли использовать что-то вроде этот вопрос SO , чтобы применитьпроизвольная функция для этого вектора.

Что-то вроде:

x_width = 15
y_width = 10

x, y = np.meshgrid(range(x_width), range(y_width))

def translate(x, y, x_o, y_o):
    x_new = x + x_o
    y_new = y + y_o
    return x_new, y_new

x_new, y_new = translate(x, y, 3, 3)

x_new[4,5], y[4,5]

(8, 4)

Answer 2

Вы должны избегать циклов и выполнять тяжелые вычисления в базовом коде C, в Numpy или в Astropy для преобразования небо / пиксель. Есть несколько вариантов сделать это с astropy.wcs.

Первый с SkyCoord. Давайте сначала создадим сетку значений для ваших индексов пикселей:

In [30]: xx, yy = np.mgrid[:5, :5] 
    ...: xx, yy 
Out[30]: 
(array([[0, 0, 0, 0, 0],
        [1, 1, 1, 1, 1],
        [2, 2, 2, 2, 2],
        [3, 3, 3, 3, 3],
        [4, 4, 4, 4, 4]]), array([[0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4]]))

Теперь мы можем создать объект SkyCoord (который является подклассом массива Numpy) из индексов пикселей с использованием wcs:

In [33]: from astropy.coordinates import SkyCoord 
    ...: sky = SkyCoord.from_pixel(xx, yy, wcs) 
    ...: sky                                                                                   
Out[33]: 
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
    [[(53.17127889, -27.78771333), (53.17127889, -27.78765778),
      (53.17127889, -27.78760222), (53.17127889, -27.78754667),
      (53.17127889, -27.78749111)],
      ....

Обратите внимание, что используется wcs.utils.skycoord_to_pixel. Этот объект также имеет метод проецирования в пиксели с помощью wcs. Я сделаю то же самое здесь для практических целей:

In [34]: sky.to_pixel(wcs)                                                                     
Out[34]: 
(array([[ 0.00000000e+00, -1.11022302e-16, -2.22044605e-16,
         -3.33066907e-16,  1.13149046e-10],
        ...
        [ 4.00000000e+00,  4.00000000e+00,  4.00000000e+00,
          4.00000000e+00,  4.00000000e+00]]),
 array([[-6.31503738e-11,  1.00000000e+00,  2.00000000e+00,
          3.00000000e+00,  4.00000000e+00],
        ...
        [-1.11457732e-10,  1.00000000e+00,  2.00000000e+00,
          3.00000000e+00,  4.00000000e+00]]))

Мы получаем набор значений с плавающей запятой для новых индексов x и y. Поэтому вам нужно будет округлить эти значения и преобразовать их в int, чтобы использовать их в качестве индексов массива.

Второй вариант - использовать функции более низкого уровня, например, wcs.pixel_to_world_values и wcs.world_to_pixel_values, который принимает массивы Nx2 и также возвращает это:

In [37]: wcs.pixel_to_world_values(np.array([xx.ravel(), yy.ravel()]).T)                       
Out[37]: 
array([[ 53.17127889, -27.78771333],
       [ 53.17127889, -27.78765778],
       [ 53.17127889, -27.78760222],
       [ 53.17127889, -27.78754667],
       ...

Answer 3

Мое предложение состоит в том, чтобы использовать многопроцессорность или многопоточность вместо цикла for, который делает параллельное выполнение и делает код быстрее ..