为每个网格单元创建粒子位置列表

Question

我正在创建电荷涂抹功能。我有一个矩阵，每行是一个带电荷和位置的粒子。然后我查看每个粒子在网格中的位置，以计算每个网格单元中有多少粒子，但我需要知道每个粒子在哪个单元中，以便我可以找到特定网格单元中每个粒子位置的平均值。我的解决方案是创建一个列表，其中行数是矩阵中网格单元的数量，并让列在x、y和z方向上的位置，但显然我不能在每个索引中附加一个以上的数字，但可能会有一些变化？对不起，我的问题是开放式的。先谢谢你

import matplotlib.pyplot as plt
import random
import numpy as np


###Initalize particle lists 
particle_arrayX=[]
particle_arrayY=[]


###The resolution
N = 10

###Number of particles
M = 1000
col=None
row=None
###Size of box
Box_size=100
###gridsize
Grid_size=Box_size/N

###Initalize particles
for i in range(M):
    particle_arrayX.append(random.random()*Box_size)
    particle_arrayY.append(random.random()*Box_size)
    
###Intialize matrix
ParticleMatrix_Inital=[[0 for i in range(N)]]*N


###Measure density in each cell
for i in range(M):
    
    col=None
    row=None
    #The x and y components are diveded by the gridsize
    #Then they are converted to integers and then asigned either to a row or column
    #If value is float with decimal 0 EX 2.0, then 1 is substracted before converted to int
    coln=particle_arrayX[i]/Grid_size
    rown=particle_arrayY[i]/Grid_size
    
    if coln.is_integer()==True:
        col=int(coln)-1
    else:
        col=int(coln)
    if rown.is_integer()==True:
        row=int(rown)-1
    else:
        row=int(rown)
    ParticleMatrix_Inital=np.array(ParticleMatrix_Inital)
    
    ParticleMatrix_Inital[row,col]=ParticleMatrix_Inital[row,col]+1
    ParticleMatrix_Inital=ParticleMatrix_Inital.tolist()

#Plot matrix
plt.imshow(ParticleMatrix_Inital)
plt.colorbar()
plt.show()

Answer 1

欢迎来到SO

有很多方法可以解决“装箱”经验数据的问题。我在下面提出了一个面向对象（OO）解决方案，因为（在我的主观观点中）它提供了干净、整洁和高度可读的代码。另一方面，如果要模拟巨大的多粒子系统，OO解决方案可能不是最有效的。如果下面的代码不能完全解决您的问题，我仍然希望它的部分内容能对您有所帮助

也就是说，我建议将您的网格实现为class。为了让我们自己的生活更轻松，我们可以应用所有粒子都有正坐标的惯例。也就是说x、y甚至z（如果引入）从0延伸到您定义的任何box_size。然而，class Grid并不真正关心实际的box_size，只关心网格的分辨率

class Grid:
    def __init__(self, _delta_x, _delta_y):
        self.delta_x = _delta_x
        self.delta_y = _delta_y
    def bounding_cell(self, x, y):
        # Note that y-coordinates corresponds to matrix rows,
        # and that x-coordinates corresponds to matrix columns
        return (int(y/self.delta_y), int(x/self.delta_x))

是，这可能是一个简单的函数。然而，作为一个类，它很容易扩展。此外，函数必须依赖于全局变量（糟糕！），或者在每个维度方向上显式地给出网格跨度（delta），以便每次确定给定坐标(x,y)所属的矩阵单元（或bin）

现在，它是如何工作的？想象一下最简单的情况，网格分辨率为1。然后，位置(x,y) = (1.2, 4,9)处的粒子应放置在(row,col) = (4,1)处的矩阵中。也就是row = int(y/delta_y)，同样地x。分辨率越高（越小delta），矩阵的行数和列数就越大

现在我们有了一个Grid，让我们也对粒子进行对象定向！直截了当地说：

class Particle:
    def __init__(self, _id, _charge, _pos_x, _pos_y):
        self.id = _id
        self.charge = _charge
        self.pos_x = _pos_x
        self.pos_y = _pos_y
    def __str__(self):
        # implementing the __str__ method let's us 'print(a_particle)'
        return "{%d, %d, %3.1f, %3.1f}" % (self.id, self.charge, self.pos_x, self.pos_y)
    def get_position(self):
        return self.pos_x, self.pos_y
    def get_charge(self):
        return self.charge

这个类或多或少只是一个数据的集合，很容易被dict所取代。然而，这个类清楚地表达了它的意图，它干净整洁，也很容易扩展

现在，让我们创建一些粒子实例！下面是一个函数，它通过列表理解创建一个具有id、charge和位置(x,y)的粒子列表：

import random
def create_n_particles(n_particles, max_pos):
    return [Particle(id, # unique ID
                     random.randint(-1,1), # charge [-1, 0, 1]
                     random.random()*max_pos, # x coord
                     random.random()*max_pos) # y coord
            for id in range(n_particles)]

最后，我们来看看有趣的部分：把它们放在一起：

import numpy as np

if __name__ == "__main__":

    n_particles = 1000
    box_size = 100
    grid_resolution = 10
    grid_size = int(box_size / grid_resolution)

    grid = Grid(grid_resolution, grid_resolution)

    particles = create_n_particles(n_particles, box_size)

    charge_matrix = np.zeros((grid_size, grid_size))

    id_matrix = [[ [] for i in range(grid_size)] for j in range(grid_size)]

    for particle in particles:
        x, y = particle.get_position()
        row, col = grid.bounding_cell(x, y)
        charge_matrix[row][col] += particle.get_charge()
        # The ID-matrix is similar to the charge-matrix,
        # but every cell contains a new list of particle IDs
        id_matrix[row][col].append(particle.id)

请注意ID矩阵的初始化：这是您要求的每个网格单元的粒子位置列表。它是一个矩阵，表示粒子容器，每个单元都包含一个要用粒子ID填充的列表。您还可以使用整个粒子实例（不仅仅是它们的ID）填充这些列表：id_matrix[row][col].append(particle)

最后一个for循环完成了真正的工作，这里的面向对象策略向我们展示了它的魅力：循环很短，很容易阅读和理解正在发生的事情：charge_matrix中的一个单元包含这个网格单元/箱中的总电荷。同时，id_matrix被包含在此网格单元/箱中的粒子ID填充

从我们构建粒子列表的方式来看，particles，我们看到一个粒子的ID与列表中该粒子的索引是等价的。因此，它们可以像这样被取回

for i,row in enumerate(id_matrix):
    for j,col in enumerate(row):
        print("[%d][%d] : " % (i, j), end="")
        for particle_id in id_matrix[i][j]:
            p = particles[particle_id]
            # do some work with 'p', or just print it:
            print(p, end=", ")
        print("") # print new line


# Output:
# [0][0] : {32, -1, 0.2, 0.4}, ... <  all data of that particle
# ....

我把这个检索的优化留给你们，因为我真的不知道你们需要什么样的数据以及你们将如何处理它。也许最好将所有粒子都包含在一个dict中，而不是列表中；我不知道（？）。你选择

最后，我建议您使用matshow来显示矩阵，而imshow更适合图像

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)
cax = ax.matshow(charge_matrix)
fig.colorbar(cax)
ax.invert_yaxis() # flip the matrix such that the y-axis points upwards
fig.savefig("charge_matrix.png")

我们还可以散点绘制粒子，并在上面的matshow中添加与我们的网格对应的网格线。我们给散点图上色，使负电荷为蓝色，中性电荷为灰色，正电荷为红色

def charge_color(charge):
    if charge > 0: return 'red'
    elif charge < 0: return 'blue'
    else: return 'gray'

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_aspect('equal')
ax.set_xticks(np.arange(0, 101, grid_resolution))
ax.set_yticks(np.arange(0, 101, grid_resolution))
ax.grid()
ax.scatter([p.pos_x for p in particles],
            [p.pos_y for p in particles],
            c=[charge_color(p.get_charge()) for p in particles])
fig.savefig("particle_dist.png")