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

Question

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

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

кода, который будет производить репрезентативные данные: (. Моя проблема очень похож на Function to compute 3D gradient with unevenly spaced sample locations но, казалось, не было никакого решения, поэтому я думал, что я снова спрашиваю)

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? 

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

Answer 1

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

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(g-y)/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 

Update
Я не знаю, 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(g-y)/np.linalg.norm(g))) 

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

 
`Given a set of M vectors in n-dimensions (call them b_k), find a set of 
`coeffs (call them a_k) which yields the best estimate of the identity 
`matrix and the zero vector 
` 
`         M 
` (1) min ||E - I||, where E = sum a_k b_k b_k 
`  a_k      k=1 
` 
`         M 
` (2) min ||z - 0||, where z = sum a_k b_k 
`  a_k      k=1 
` 
` 
`Note that the basis vectors {b_k} are not required 
`to be normalized, orthogonal, or even linearly independent. 
` 
`First, define the following quantities: 
` 
` B    ==> matrix whose columns are the b_k 
` G = B'.B  ==> transpose of B times B 
` F = G @ G  ==> @ represents the hadamard product 
` v = diag(G) ==> vector composed of diag elements of G 
` 
`The above minimizations are equivalent to this linearly constrained problem 
` 
` Solve F.a = v 
` s.t. G.a = 0 
` 
`Let {X} denote the Moore-Penrose inverse of X. 
`Then the solution of the linear problem can be written: 
` 
` N = I - {G}.G  ==> projector into nullspace of G 
` a = N . {F.N} . v 
` 
`The utility of these coeffs is that they allow you to write 
`very simple expressions for the derivatives of a tensor field. 
` 
` 
`Let D be the del (or nabla) operator 
`and d be the difference operator wrt the central (aka 0th) node, 
`so that, for any scalar/vector/tensor quantity Y, we have: 
` dY = Y - Y_0 
` 
`Let x_k be the position vector of the kth node. 
`And for our basis vectors, take 
` b_k = dx_k = x_k - x_0. 
` 
`Assume that each node has a field value associated with it 
` (e.g. temperature), and assume a quadratic model [about x = x_0] 
` for the field [g=gradient, H=hessian, ":" is the double-dot product] 
` 
`  Y = Y_0 + (x-x_0).g + (x-x_0)(x-x_0):H 
` dY = dx.g + dxdx:H 
` D2Y = 2*I:H   ==> Laplacian of Y 
` 
` 
`Evaluate the model at the kth node 
` 
` dY_k = dx_k.g + dx_k dx_k:H 
` 
`Multiply by a_k and sum 
` 
`  M    M     M 
` 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 
` 
`Thus, we have a second order estimate of the Laplacian 
` 
`    M 
` Lap(Y_0) = sum 2 a_k dY_k 
`    k=1 
` 
` 
`Now play the same game with a linear model 
` dY_k = dx_k.g 
` 
`But this time multiply by (a_k dx_k) and sum 
` 
`  M     M 
` sum a_k dx_k dY_k = sum a_k dx_k dx_k.g 
` k=1     k=1 
` 
`      = I.g 
`      = g 
` 
` 
`In general, the derivatives at the central node can be estimated as 
` 
`   M 
` D#Y = sum a_k dx_k#dY_k 
`   k=1 
` 
`   M 
` D2Y = sum 2 a_k dY_k 
`   k=1 
` 
` where 
` # stands for the {dot, cross, or tensor} product 
`  yielding the {div, curl, or grad} of Y 
` and 
` D2Y stands for the Laplacian of Y 
` D2Y = D.DY = Lap(Y) 

Answer 2

Наглядно для деривата WRT один DataPoint, я хотел бы сделать следующий

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

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

EDIT:

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

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

Answer 3

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

1. я предлагаю использовать query(x,distance_upper_bound) вместо query(x,k), так как это, вероятно, предотвратит некоторые неустойчивости из-за кластеризации

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

Answer 4

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

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

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

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