Блог

Доброе утро, я внедряю фильтр Крессмана для вычисления средневзвешенных значений расстояния в Numpy.. Я использую имплиментацию шарового дерева (благодаря Джейку Вандерпласу), чтобы возвращать список местоположений для каждой точки в массиве запросов. массив запросов (q) имеет форму [n, 3] и в каждой точке имеет x, y, z в точке, в которой я хочу получить средневзвешенное значение точек, хранящихся в дереве.. код, обернутый вокруг дерева, возвращает точки на определенном расстоянии, поэтому я получаю массивы массивов переменной длины.. Я использую where для поиска непустых записей (т. Е. позиций, где было по крайней мере несколько точек в радиусе влияния), создавая массив isgood…

Затем я перебираю все точки запроса, чтобы вернуть средневзвешенное значение значений self.z (обратите внимание, что это может быть либо dims = 1, либо dims = 2, чтобы разрешить множественную совместную привязку)

итак, что усложняет использование map или других более быстрых методов, так это неравномерность длин массивов в пределах self.distances и self.locations… Я все еще довольно новичок в numpy / python, но я не могу придумать, как сделать этот массив разумным (т. Е. Не возвращаться к циклам)

 self.locations, self.distances = self.tree.query_radius( q, r, return_distance=True)
t2=time()
if debug: print "Removing voids"
isgood=np.where( np.array([len(x) for x in self.locations])!=0)[0]
interpol = np.zeros( (len(self.locations),)   np.shape(self.z[0]) )
interpol.fill(np.nan)
for dist, ix, posn, roi in zip(self.distances[isgood], self.locations[isgood], isgood, r[isgood]):
    interpol[isgood[jinterpol]] = np.average(self.z[ix], weights=(roi**2-dist**2) / (roi**2   dist**2), axis=0)
    jinterpol  = 1
  

итак… Есть какие-либо подсказки о том, как ускорить цикл?..

Для типичного отображения, применяемого для отображения данных метеорологического радара из сетки дальности, азимута, высоты в декартову сетку, где у меня 240x240x34 точек и 4 переменных, требуется 99 секунд для запроса дерева (написанного Джейком на C и cython .. это сложный шаг, поскольку вам нужно выполнить поиск данных!) и 100 секунд для выполнения вычисления… что, на мой взгляд, является медленным?? где мои накладные расходы? эффективен ли np.mean или, как его называют в миллионы раз, можно ли здесь добиться ускорения? выиграю ли я, используя float32, а не default64 … или даже масштабирование до целых чисел (чего было бы очень сложно избежать при взвешивании… любые подсказки с благодарностью получены!

Вы можете найти обсуждение относительных преимуществ схемы Крессмана по сравнению с использованием весовой функции Гаусса на:

http://www.flame.org /~cdoswell/публикации/radar_oa_00.pdf

Главное — сопоставить параметр сглаживания с данными (я рекомендую использовать значение, близкое к среднему расстоянию между точками данных). Как только вы узнаете параметр сглаживания, вы можете установить «радиус влияния», равный радиусу, при котором весовая функция падает до 0,01 (или любого другого значения).

Насколько важна скорость? Если вы хотите, вместо вызова экспоненциальной функции для определения веса, вы можете составить дискретную таблицу весов для некоторого фиксированного числа приращений радиуса, что значительно ускоряет вычисление. В идеале у вас должны быть данные за пределами границ сетки, которые можно использовать при отображении значений, окружающих точки сетки (даже в граничных точках сетки). Обратите внимание, что это не настоящая схема интерполяции — она не вернет точно наблюдаемые значения в точках данных. Как и схема Крессмана, это файлер нижних частот.