Проецирование скорости на линии прямой видимости - РЕДАКТИРОВАННЫЙ

Question

У меня есть огромный текстовый файл, который содержит положение (x, y, z) и компонент скорости (vx, vy, vz) миллиона звезд. После некоторых вращений и проекций я получаю новые компоненты положения и скорости (x', y', z', vx', vy', vz') звезд.

My последний шаг - вычислить скорость вдоль луча зрения, это как будто я должен «усреднить» компонент vz, и для этого я пытаюсь создать файл FITS, в котором каждый пиксель содержит среднее значение vz component.

Вот часть моего кода:

mod = np.genfromtxt('data_bar_region.txt')

x = list(mod[:,0])
y = list(mod[:,1])
vz = mod[:,5]

x_rang_1 = np.arange(-40, 41, 1)
y_rang_1 = np.arange(-40, 41, 1)

fake_data_1 = np.array((len(x_rang_1),len(x_rang_1)))

for i in range(len(x_rang_1)-1):
    for j in range(len(y_rang_1)-1):
        vel_tmp = []
        for index in range(len(x)):
            if x_rang_1[i] <= x[index] <= x_rang_1[i+1]:
                if y_rang_1[j] <= y[index] <= y_rang_1[j+1]:
                    vel_tmp.append(vz[index])
        fake_data_1[j,i] = np.mean(vel_tmp)
            

hdu1 = fits.PrimaryHDU(fake_data_1)
hdu1.writeto('TEST.fits')

Этот код слишком медленный (на моем ноутбуке на это ушло около 8 часов), и я не знаю, как я может ускориться.

Были ли у вас предложения или другие способы вычисления v_LOS лучше и быстрее?

РЕДАКТИРОВАТЬ: Перед выполнением «усреднения» я должен разделить изображение на части нескольких форм и размеров (такие части называются «ячейками»).

Здесь есть [изображение ячеек ( на правой панели есть то же изображение ящиков, но оно увеличено, чтобы лучше понять, что такое ящики)] 1 .

Итак, у меня есть еще один файл FITS (называемый bins.fits) с тот же размер fake_data_1, и я просто хочу найти соответствие между этими двумя файлами, потому что я хочу вычислить среднее и стандартное распределение звезд в нескольких ячейках.

В качестве альтернативы я иметь текстовый файл, содержащий информацию о том, какой пиксель принадлежит указанному c bin, например:

xy bin 1 1 34 1 2 34 1 3 34 . . . 34 56 37 34 57 37 34 58 37

и так далее. Файл bins.file имеет размер (564,585), как и fake_data_1, изменяя возможность начала и конца диапазона x и y. Прикрепил весь скрипт:

mod = np.genfromtxt('data_new_bar_scaled.txt')

# to match the correct position and size of the observation,
# I have to multiply by a factor equal to the semi-size
x = mod[:, 0]*(585-1)/200
y = mod[:, 1]*(564-1)/200
vz = mod[:,5]

A = fits.open('NGC4277_TESIkinematic.fits')
bins = A[7].data.T


start_x = -(585-1)/2
stop_x = (585-1)/2
step_x = step  # step in x_rang_1
x_rang = np.arange(start_x, stop_x + step_x, step_x)    

start_y = -(564-1)/2
stop_y = (564-1)/2
step_y = step  # step in y_rang_1
y_rang = np.arange(start_y, stop_y + step_y, step_y)

fake_data_1 = np.empty((len(x_rang), len(y_rang)))
fake_data_1[:] = np.NaN  # initialize with NaN

 
print(fake_data_1.shape)
print(bins.shape)


d = {}
for i in range(len(x)):
    index_for_x = math.floor((x[i] - start_x) / step_x)
    index_for_y = math.floor((y[i] - start_y) / step_y)
    if 0 <= index_for_x < len(x_rang) and 0 <= index_for_y < len(y_rang):
        key = (x_rang[index_for_x], y_rang[index_for_y]) 
        
        if key in d:
            d[key].append(vz[i])
        else:
            d[key] = [vz[i]]

bb = np.unique(bins)
print(len(bb))

for i, x in enumerate(x_rang):
    for j, y in enumerate(y_rang):
        key = (x, y) 
    
        for z in range(len(bb)):
            j,k = np.where(bb[z]==bins)
            print('index :', z)
            if key in d:
                fake_data_1[j,k] = np.mean(d[key])

Answer 1

Ваш код настолько медленный, поскольку вложенные циклы в вашем коде повторяются более миллиона звезд 1600 (80 * 80) раз. Вы можете улучшить производительность, используя словарь и перебирая более миллиона звезд только один раз.

Вы можете попробовать следующий код, который примерно в 1600 раз быстрее:

import numpy as np
import math

mod = np.genfromtxt('data_bar_region.txt')

x = list(mod[:, 0])
y = list(mod[:, 1])
vz = mod[:, 5]

x_rang_1 = np.arange(-40, 41, 1)
y_rang_1 = np.arange(-40, 41, 1)

fake_data_1 = np.empty((len(x_rang_1), len(y_rang_1)))
fake_data_1[:] = np.NaN  # initialize with NaN

d = {}
for i in range(len(x)):
    key = (math.floor(x[i]), math.floor(y[i]))
    if key in d:
        d[key].append(vz[i])
    else:
        d[key] = [vz[i]]

for i, x in enumerate(x_rang_1):
    for j, y in enumerate(y_rang_1):
        key = (x, y)
        if key in d:
            fake_data_1[i, j] = np.mean(d[key])

hdu1 = fits.PrimaryHDU(fake_data_1)
hdu1.writeto('TEST.fits')

ОБНОВЛЕНИЕ

Для обобщенной версии для step в x_rang_1 (или y_rang_1) вы можете попробовать следующий код:

import numpy as np
import math

mod = np.genfromtxt('data_bar_region.txt')

x = list(mod[:, 0])
y = list(mod[:, 1])
vz = mod[:, 5]

start_x_rang_1 = -40
stop_x_rang_1 = 40
step_x_rang_1 = 0.5  # step in x_rang_1
x_rang_1 = np.arange(start_x_rang_1, stop_x_rang_1 + step_x_rang_1, step_x_rang_1)

start_y_rang_1 = -40
stop_y_rang_1 = 40
step_y_rang_1 = 1  # step in y_rang_1
y_rang_1 = np.arange(start_y_rang_1, stop_y_rang_1 + step_y_rang_1, step_y_rang_1)

fake_data_1 = np.empty((len(x_rang_1), len(y_rang_1)))
fake_data_1[:] = np.NaN  # initialize with NaN

d = {}
for i in range(len(x)):
    index_for_x_rang_1 = math.floor((x[i] - start_x_rang_1) / step_x_rang_1)
    index_for_y_rang_1 = math.floor((y[i] - start_y_rang_1) / step_y_rang_1)
    if 0 <= index_for_x_rang_1 < len(x_rang_1) and 0 <= index_for_y_rang_1 < len(y_rang_1):
        key = (x_rang_1[index_for_x_rang_1], y_rang_1[index_for_y_rang_1])
        if key in d:
            d[key].append(vz[i])
        else:
            d[key] = [vz[i]]

for i, x in enumerate(x_rang_1):
    for j, y in enumerate(y_rang_1):
        key = (x, y)
        if key in d:
            fake_data_1[i, j] = np.mean(d[key])

hdu1 = fits.PrimaryHDU(fake_data_1)
hdu1.writeto('TEST.fits')

ОБНОВЛЕНИЕ 2

Может быть, как показано ниже?

Когда я предполагал, что входы

x    y    vz
0    0.1  10
1.8  0    4
1.2  1.9  5.2

bins = np.array(
    [[34, 35, 34, 34, 36],
     [37, 36, 34, 35, 36],
     [34, 35, 37, 36, 34]])  # shape: (5, 3)

Вам нужны следующие код?

import numpy as np
import math

x = np.array([0, 1.8, 1.2, ])
y = np.array([0.1, 0, 1.9, ])
vz = np.array([10, 4, 5.2])

start_x_rang_1 = 0
stop_x_rang_1 = 2
step_x_rang_1 = 1  # step in x_rang_1
x_rang_1 = np.arange(start_x_rang_1, stop_x_rang_1 + step_x_rang_1, step_x_rang_1)

start_y_rang_1 = 0
stop_y_rang_1 = 0.5
step_y_rang_1 = 2  # step in y_rang_1
y_rang_1 = np.arange(start_y_rang_1, stop_y_rang_1 + step_y_rang_1, step_y_rang_1)

fake_data_1 = np.empty((len(x_rang_1), len(y_rang_1)))  # shape: (3, 5)
fake_data_1[:] = np.NaN  # initialize with NaN
bins = np.array(
    [[34, 35, 34, 34, 36],
     [37, 36, 34, 35, 36],
     [34, 35, 37, 36, 34]])  # shape: (3, 5)

d_bins = {}
for i in range(len(x)):
    index_for_x_rang_1 = math.floor((x[i] - start_x_rang_1) / step_x_rang_1)
    index_for_y_rang_1 = math.floor((y[i] - start_y_rang_1) / step_y_rang_1)
    if 0 <= index_for_x_rang_1 < len(x_rang_1) and 0 <= index_for_y_rang_1 < len(y_rang_1):
        key = bins[index_for_x_rang_1, index_for_y_rang_1]
        if key in d_bins:
            d_bins[key].append(vz[i])
        else:
            d_bins[key] = [vz[i]]

d_bins_mean = {}
for bin in d_bins:
    d_bins_mean[bin] = np.mean(d_bins[bin])

get_corresponding_mean = np.vectorize(lambda x: d_bins_mean.get(x, np.NaN))
result = get_corresponding_mean(bins)
print(result)

который печатает

[[10.   nan 10.  10.   nan]
 [ 4.6  nan 10.   nan  nan]
 [10.   nan  4.6  nan 10. ]]