Вычисление 3D-градиента с неравномерно разнесенными точками

В настоящее время у меня есть объем, наложенный несколькими миллионами неравномерно распределенных частиц, и каждая частица имеет атрибут (потенциальный, для любопытных), который я хочу вычислить для локальной силы (ускорения).

np.gradient работает только с равномерно распределенными данными, и я смотрел здесь: Градиент второго порядка в numpy, где требуется интерполяция, но я не смог найти реализацию 3D сплайна в Numpy.

Некоторый код, который будет выдавать репрезентативные данные:

import numpy as np from scipy.spatial import cKDTree x = np.random.uniform(-10, 10, 10000) y = np.random.uniform(-10, 10, 10000) z = np.random.uniform(-10, 10, 10000) phi = np.random.uniform(-10**9, 0, 10000) kdtree = cKDTree(np.c_[x,y,z]) _, index = kdtree.query([0,0,0], 32) #find 32 nearest particles to the origin #find the gradient at (0,0,0) by considering the 32 nearest particles? 

(Моя проблема очень похожа на функцию для вычисления 3D-градиента с неравномерно расположенными точками выборки, но там, казалось, не было никакого решения, поэтому я подумал, что снова попрошу.)

Любая помощь будет оценена по достоинству.

Вот реализация Джулии, которая делает то, что вы просите

 using NearestNeighbors n = 3; k = 32; # for stability use k > n*(n+3)/2 # Take a point near the center of cube point = 0.5 + rand(n)*1e-3; data = rand(n, 10^4); kdtree = KDTree(data); idxs, dists = knn(kdtree, point, k, true); # Coords of the k-Nearest Neighbors X = data[:,idxs]; # Least-squares recipe for coefficients C = point * ones(1,k); # central node dX = X - C; # diffs from central node G = dX' * dX; F = G .* G; v = diag(G); N = pinv(G) * G; N = eye(N) - N; a = N * pinv(F*N) * v; # ...these are the coeffs # Use a temperature distribution of T = 25.4 * r^2 # whose analytical gradient is gradT = 25.4 * 2*x X2 = X .* X; C2 = C .* C; T = 25.4 * n * mean(X2, 1)'; Tc = 25.4 * n * mean(C2, 1)'; # central node dT = T - Tc; # diffs from central node y = dX * (a .* dT); # Estimated gradient g = 2 * 25.4 * point; # Analytical # print results @printf "Estimated Grad = %s\n" string(y') @printf "Analytical Grad = %s\n" string(g') @printf "Relative Error = %.8f\n" vecnorm(gy)/vecnorm(g) 

Этот метод имеет относительную погрешность 1%. Вот результаты нескольких запусков …

 Estimated Grad = [25.51670916224472 25.421038632006926 25.6711949674633] Analytical Grad = [25.41499027802736 25.44913042322385 25.448202594123806] Relative Error = 0.00559934 Estimated Grad = [25.310574056859014 25.549736360607493 25.368056350800604] Analytical Grad = [25.43200914200516 25.43243178887198 25.45061497749628] Relative Error = 0.00426558 

Обновить
Я не очень хорошо знаю Python, но вот перевод, который, кажется, работает

 import numpy as np from scipy.spatial import KDTree n = 3; k = 32; # fill the cube with random points data = np.random.rand(10000,n) kdtree = KDTree(data) # pick a point (at the center of the cube) point = 0.5 * np.ones((1,n)) # Coords of k-Nearest Neighbors dists, idxs = kdtree.query(point, k) idxs = idxs[0] X = data[idxs,:] # Calculate coefficients C = (np.dot(point.T, np.ones((1,k)))).T # central node dX= X - C # diffs from central node G = np.dot(dX, dX.T) F = np.multiply(G, G) v = np.diag(G); N = np.dot(np.linalg.pinv(G), G) N = np.eye(k) - N; a = np.dot(np.dot(N, np.linalg.pinv(np.dot(F,N))), v) # these are the coeffs # Temperature distribution is T = 25.4 * r^2 X2 = np.multiply(X, X) C2 = np.multiply(C, C) T = 25.4 * n * np.mean(X2, 1).T Tc = 25.4 * n * np.mean(C2, 1).T # central node dT = T - Tc; # diffs from central node # Analytical gradient ==> gradT = 2*25.4* x g = 2 * 25.4 * point; print( "g[]: %s" % (g) ) # Estimated gradient y = np.dot(dX.T, np.multiply(a, dT)) print( "y[]: %s, Relative Error = %.8f" % (y, np.linalg.norm(gy)/np.linalg.norm(g)) ) 

Обновление # 2
Я думаю, что я могу написать что-то понятное, используя форматированный ASCII вместо LaTeX …

 `Учитывая набор M векторов в n-измерениях (назовем их b_k), найдите набор
 `coeffs (назовем их a_k), что дает наилучшую оценку тождества
 `матрица и нулевой вектор
 `
 `M
 `(1) min || E - I ||, где E = sum a_k b_k b_k
 `a_k k = 1
 `
 `M
 (2) min || z - 0 ||, где z = sum a_k b_k
 `a_k k = 1
 `
 `
 `Обратите внимание, что базисные векторы {b_k} не требуются
 `быть нормированным, ортогональным или даже линейно независимым.
 `
 `Сначала определите следующие величины:
 `
 `B ==> матрица, столбцами которой являются b_k
 `G = B'.B ==> транспонирование B раз B
 `F = G @ G ==> @ представляет продукт hadamard
 `v = diag (G) ==> вектор, состоящий из элементов diag из G
 `
 `Вышеуказанные минимизации эквивалентны этой линейно ограниченной задаче
 `
 `Решите Fa = v
 `st Ga = 0
 `
 «Пусть {X} обозначает Мура-Пенроуз, обратный X.
 «Тогда решение линейной задачи можно записать:
 `
 `N = I - {G} .G ==> проектор в пустое пространство G
 `a = N.  {FN}.  v
 `
 «Полезность этих коэффициентов заключается в том, что они позволяют писать
 `очень простые выражения для производных тензорного поля.
 `
 `
 `Пусть D - оператор del (или nabla)
 `и d - оператор разности по центральному (aka 0) узлу,
 так что для любой скалярной / векторной / тензорной величины Y имеем:
 `dY = Y - Y_0
 `
 `Пусть x_k - вектор положения k-го узла.
 `И для наших базисных векторов возьмем
 `b_k = dx_k = x_k - x_0.
 `
 `Предположим, что каждый узел имеет связанное с ним значение поля
 (например, температура) и принять квадратичную модель [о х = х_0]
 `для поля [g = градиент, H = hessian,": "- произведение с двойной точкой)
 `
 `Y = Y_0 + (x-x_0) .g + (x-x_0) (x-x_0): H
 `dY = dx.g + dxdx: H
 `D2Y = 2 * I: H ==> Лапласиан Y
 `
 `
 `Оценить модель на k-м узле 
 `
 `dY_k = dx_k.g + dx_k dx_k: H
 `
 `Умножение на a_k и sum
 `
 `MMM
 `sum a_k dY_k = sum a_k dx_k.g + sum a_k dx_k dx_k: H
 `k = 1 k = 1 k = 1
 `
 `= 0.g + I: H
 `= D2Y / 2
 `
 Таким образом, мы имеем оценку второго порядка лапласиана
 `
 `M
 `Lap (Y_0) = сумма 2 a_k dY_k
 `k = 1
 `
 `
 `Теперь играйте в ту же игру с линейной моделью
 `dY_k = dx_k.g
 `
 `Но на этот раз умножим на (a_k dx_k) и сумму
 `
 `ММ
 `sum a_k dx_k dY_k = sum a_k dx_k dx_k.g
 `k = 1 k = 1
 `
 `= Ig
 `= g
 `
 `
 `В общем случае производные в центральном узле могут быть оценены как
 `
 `M
 `D # Y = sum a_k dx_k # dY_k
 `k = 1
 `
 `M
 `D2Y = сумма 2 a_k dY_k
 `k = 1
 `
 `где
 `# обозначает произведение {точка, крест или тензор}
 `уступая {div, curl, или grad} Y
 и
 `D2Y означает лапласиан Y
 `D2Y = D.DY = Lap (Y)

Интуитивно, для деривации по одному datapoint, я бы сделал следующее

  • Возьмите фрагмент окружающих данных: data=phi[x_id-1:x_id+1, y_id-1:y_id+1, z_id-1:z_id+1] . Подход с kdTre выглядит очень красиво, конечно, вы также можете использовать его для подмножества данных.
  • Установите 3D-полином, вы можете посмотреть на polyvander3D . Определите точку в середине среза как центр. Вычислите смещения в других точках. Передайте их в виде координат для полифита.
  • Выведите полином в свое положение.

Это будет простым решением вашей проблемы. Однако это, вероятно, будет очень медленным.

РЕДАКТИРОВАТЬ:

На самом деле это, по-видимому, обычный метод: https://scicomp.stackexchange.com/questions/480/how-can-i-numerically-differentiate-an-unevenly-sampled-function

В принятом ответе говорится о выводе интерполяционного многочлена. Хотя очевидно, что полином должен охватывать все данные (матрица Вандермонда). Для вас это невозможно, слишком много данных. Взятие локального подмножества кажется очень разумным.

Многое зависит от соотношения сигнал / шум ваших потенциальных данных. Ваш пример – это все шумы, поэтому «подгонка» чего-либо к нему всегда будет «чрезмерной». Степень шума будет определять степень, в которой вы хотите быть поли-фитингом (как с ответом lhk), и насколько вы хотите быть Kriging (используя pyKriging или иначе)

  1. Я бы предложил использовать query(x,distance_upper_bound) вместо query(x,k) , так как это, вероятно, предотвратит некоторые нестабильности из-за кластеризации

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

Я дам свои два цента поздно. В случае, когда пространство равномерно растянуто и велико, обычно извлекается только локальная информация для каждой частицы.

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

  1. N ближайшего соседа, используя, например, дерево KD. Это определяет местность динамически, что может быть или не быть хорошей идеей.
  2. Случайно разбиваем пространство с плоскостями на группы частиц. В основном тестирование на N-неравенство сокращает пространство N раз.

После определения локальности вы можете интерполировать полином, который аналитически дифференцируется. Я призываю больше думать в разных определениях мест. (может дать интересное различие)