У меня есть несколько точек (координаты x, y, z) в 3D-окне с соответствующими массами. Я хочу нарисовать гистограмму плотности массы, которая находится в сферах данного радиуса R.
Я написал код, который при условии, что я не допустил ошибок, которые, как мне кажется, могут возникнуть, работает следующим образом:
Мои "настоящие" данные - это что-то огромное, поэтому я написал небольшой код, чтобы случайным образом генерировать непересекающиеся точки с произвольной массой в коробке.
Я вычисляю трехмерную гистограмму (взвешенную по массе) с биннингом, примерно в 10 раз меньшим, чем радиус моих сфер.
Я беру БПФ своей гистограммы, вычисляю волновые моды (kx, ky и kz) и использую их для умножения моей гистограммы в пространстве Фурье на аналитическое выражение трехмерной вершины -это функция окна (сферическая фильтрация) в пространстве Фурье.
Я инвертирую БПФ мою недавно вычисленную сетку.
Таким образом, построение 1D-гистограммы значений на каждом бине даст мне то, что я хочу.
Моя проблема заключается в следующем: учитывая то, что я делаю, в моей инвертированной сетке БПФ не должно быть никаких отрицательных значений (шаг 4), но я получаю некоторые и со значениями, намного превышающими числовую ошибку.
Если я запускаю свой код на маленькой коробке (300x300x300 см 3 и точки, разделенные как минимум на 1 см), я не получаю проблему. Хотя я получаю его за 600x600x600 см3.
Если я установлю все массы на 0, таким образом работая на пустой сетке, я получу обратно 0 без каких-либо замеченных проблем.
Я даю свой код в полном блоке, чтобы его можно было легко скопировать.
import numpy as np
import matplotlib.pyplot as plt
import random
from numba import njit
# 1. Generate a bunch of points with masses from 1 to 3 separated by a radius of 1 cm
radius = 1
rangeX = (0, 100)
rangeY = (0, 100)
rangeZ = (0, 100)
rangem = (1,3)
qty = 20000 # or however many points you want
# Generate a set of all points within 1 of the origin, to be used as offsets later
deltas = set()
for x in range(-radius, radius+1):
for y in range(-radius, radius+1):
for z in range(-radius, radius+1):
if x*x + y*y + z*z<= radius*radius:
deltas.add((x,y,z))
X = []
Y = []
Z = []
M = []
excluded = set()
for i in range(qty):
x = random.randrange(*rangeX)
y = random.randrange(*rangeY)
z = random.randrange(*rangeZ)
m = random.uniform(*rangem)
if (x,y,z) in excluded: continue
X.append(x)
Y.append(y)
Z.append(z)
M.append(m)
excluded.update((x+dx, y+dy, z+dz) for (dx,dy,dz) in deltas)
print("There is ",len(X)," points in the box")
# Compute the 3D histogram
a = np.vstack((X, Y, Z)).T
b = 200
H, edges = np.histogramdd(a, weights=M, bins = b)
# Compute the FFT of the grid
Fh = np.fft.fftn(H, axes=(-3,-2, -1))
# Compute the different wave-modes
kx = 2*np.pi*np.fft.fftfreq(len(edges[0][:-1]))*len(edges[0][:-1])/(np.amax(X)-np.amin(X))
ky = 2*np.pi*np.fft.fftfreq(len(edges[1][:-1]))*len(edges[1][:-1])/(np.amax(Y)-np.amin(Y))
kz = 2*np.pi*np.fft.fftfreq(len(edges[2][:-1]))*len(edges[2][:-1])/(np.amax(Z)-np.amin(Z))
# I create a matrix containing the values of the filter in each point of the grid in Fourier space
R = 5
Kh = np.empty((len(kx),len(ky),len(kz)))
@njit(parallel=True)
def func_njit(kx, ky, kz, Kh):
for i in range(len(kx)):
for j in range(len(ky)):
for k in range(len(kz)):
if np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2) != 0:
Kh[i][j][k] = (np.sin((np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2))*R)-(np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2))*R*np.cos((np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2))*R))*3/((np.sqrt(kx[i]**2+ky[j]**2+kz[k]**2))*R)**3
else:
Kh[i][j][k] = 1
return Kh
Kh = func_njit(kx, ky, kz, Kh)
# I multiply each point of my grid by the associated value of the filter (multiplication in Fourier space = convolution in real space)
Gh = np.multiply(Fh, Kh)
# I take the inverse FFT of my filtered grid. I take the real part to get back floats but there should only be zeros for the imaginary part.
Density = np.real(np.fft.ifftn(Gh,axes=(-3,-2, -1)))
# Here it shows if there are negative values the magnitude of the error
print(np.min(Density))
D = Density.flatten()
N = np.mean(D)
# I then compute the histogram I want
hist, bins = np.histogram(D/N, bins='auto', density=True)
bin_centers = (bins[1:]+bins[:-1])*0.5
plt.plot(bin_centers, hist)
plt.xlabel('rho/rhom')
plt.ylabel('P(rho)')
plt.show()
Знаете ли вы, почему я получаю эти отрицательные значения? Как вы думаете, есть более простой способ продолжить?
Извините, если это очень длинный пост, я постарался сделать его очень понятным и отредактирую его с вашими комментариями, большое спасибо!
-EDIT-
Дополнительный вопрос по этому вопросу можно найти [здесь]. 1