Рассчитайте астрономическое расстояние двух наборов положения неба с помощью theano

Question

Я хочу вычислить угловое расстояние между всеми точками в двух разных наборах, что-то вроде cdist из scipy, но с другим алгоритмом расстояния и с использованием theano.Угловое расстояние между двумя источниками с прямым восхождением (ra) в (0,2pi) и с склонением (dec) в (-pi / 2, pi / 2) составляет:

theta = arccos(sin(dec1)*sin(dec2)+cos(dec1)*cos(dec2)*cos(ra1-ra2))

предположим, что X - это матрица, состоящая из N источников с их положением (ra, dec):

#RA     DEC 
54.29   -35.19  
54.62   -35.45
...

и W - другой набор источников M различных источников.Как я могу определить угловое разделение всех X источников со всеми W источниками?

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

edist = T.sqrt((X ** 2).sum(1).reshape((X.shape[0], 1)) + (W ** 2).sum(1).reshape((1, W.shape[0])) - 2 * X.dot(W.T))

Я пробовал с:

d = T.arccos(\\
T.sin(X.reshape((X.shape[0], 1, -1))[...,1])*T.sin(W.reshape((1, W.shape[0], -1))[..., 1])+\\
T.cos(X.reshape((X.shape[0], 1, -1))[...,1])*T.cos(W.reshape((1, W.shape[0], -1))[..., 1])*\\
T.cos(X.reshape((X.shape[0], 1, -1))[...,0] -W.reshape((1, W.shape[0], -1))[...,0]))

эта результирующая матрица d имеет форму (N, M) вместо (N, M, 2), поскольку я ожидал суммирования по третьей оси;далее численный результат неверен (я сравнил его с TOPCAT , который ориентирован на астрономию программного обеспечения. Любое предложение?

Answer 1

Я решил проблему: просто мне нужно преобразовать правильное восхождение и склонение из градуса в радиан.Теперь метод работает.

Answer 2

Вам нужно отладить выражение по частям - сначала рассчитайте sin(dec1) и убедитесь, что вы получите правильную форму и правильный числовой результат. Затем умножьте на sin(dec2) и так далее, пока не получите полное выражение arccos.

Одной из идей о том, что, возможно, неправильно в вашем коде, является использование * для умножения - если вы хотите умножить матрицы, вы должны использовать T.multiply() вместо *.