1

标题可能含糊不清，不知道该怎么说。

我在 python 中使用 numpy 和 matplotlib 的粒子模拟器已经有点远了，我已经设法实现了 coloumb、重力和风，现在我只想添加温度和压力，但我有一个优化前的问题（万恶之源）。我想看看粒子何时崩溃：

问：在 numpy 中是否可以根据 bool 条件来获取数组与每个元素的差异？我想避免循环。

例如：(x - any element in x) < a 应该返回类似

[True, True, False, True]

如果 x 中的元素 0,1 和 3 满足条件。

编辑：

等效循环是：

for i in len(x):
  for j in in len(x):
    #!= not so important
    ##earlier question I asked lets me figure that one out
    if i!=j: 
      if x[j] - x[i] < a:
       True

我注意到 numpy 操作比 if 测试快得多，这帮助我加快了速度。

如果有人想玩它，这是一个示例代码。

#Simple circular box simulator, part of part_sim
#Restructure to import into gravity() or coloumb () or wind() or pressure()
#Or to use all forces: sim_full()
#Note: Implement crashing as backbone to all forces

import numpy as np
import matplotlib.pyplot as plt

N = 1000        #Number of particles
R = 8000        #Radius of box          
r = np.random.randint(0,R/2,2*N).reshape(N,2)
v = np.random.randint(-200,200,r.shape)
v_limit = 10000 #Speedlimit


plt.ion()
line, = plt.plot([],'o')
plt.axis([-10000,10000,-10000,10000])

while True:
    r_hit = np.sqrt(np.sum(r**2,axis=1))>R   #Who let the dogs out, who, who?
    r_nhit = ~r_hit                     
    N_rhit = r_hit[r_hit].shape[0]
    r[r_hit] = r[r_hit] - 0.1*v[r_hit]       #Get the dogs back inside
    r[r_nhit] = r[r_nhit] +0.1*v[r_nhit]
    #Dogs should turn tail before they crash!
    #---
    #---crash code here....
    #---crash end
    #---
    vmin, vmax = np.min(v), np.max(v)        
    #Give the particles a random kick when they hit the wall
    v[r_hit]  = -v[r_hit] + np.random.randint(vmin, vmax, (N_rhit,2))
    #Slow down honey
    v_abs = np.abs(v) > v_limit
    #Hit the wall at too high v honey? You are getting a speed reduction
    v[v_abs] *=0.5
    line.set_ydata(r[:,1])
    line.set_xdata(r[:,0])
    plt.draw()

一旦我弄清楚如何......以便可以在更大的盒子中轻松区分高速粒子，我计划为上面的数据点添加颜色。

4

1 回答 1

5

例如：x - x < a 中的任何元素应该返回类似

【对，对，对，对，对】

如果 x 中的元素 0,1 和 3 满足条件。我注意到 numpy 操作比 if 测试快得多，这帮助我加快了速度。

是的，它只是m < a。例如：

>>> m = np.array((1, 3, 10, 5))
>>> a = 6
>>> m2 = m < a
>>> m2
array([ True,  True, False,  True], dtype=bool)

现在，对于这个问题：

问：在 numpy 中是否可以根据 bool 条件来获取数组与每个元素的差异？我想避免循环。

我不确定您在这里要求什么，但它似乎与下面的示例不匹配。例如，您是否尝试从满足谓词的每个元素中减去 1？在这种情况下，您可以依赖以下事实，False==0并且True==1只需减去布尔数组：

>>> m3 = m - m2
>>> m3
>>> array([ 0,  2, 10,  4])

根据您的说明，您需要此伪代码循环的等效项：

for i in len(x):
  for j in in len(x):
    #!= not so important
    ##earlier question I asked lets me figure that one out
    if i!=j: 
      if x[j] - x[i] < a:
        True

我认为这里的困惑在于这与您所说的完全相反：您不希望“基于布尔条件的数组与其每个元素的差异”，而是“基于差异的布尔条件具有每个元素的数组”。即使这样也只会让你真正得到一个 len(m)*len(m) bools 的方阵，但我认为剩下的部分是“任何”。

无论如何，您要求的是隐式笛卡尔积，将 m 的每个元素与 m 的每个元素进行比较。

您可以轻松地将其从两个循环减少到一个（或者，更确切地说，隐式矢量化其中一个，获得通常的 numpy 性能优势）。对于每个值，通过从每个元素中减去该值并将结果与​​ 进行比较来创建一个新数组a，然后将它们连接起来：

>>> a = -2
>>> comparisons = np.array([m - x < a for x in m])
>>> flattened = np.any(comparisons, 0)
>>> flattened
array([ True,  True, False,  True], dtype=bool)

但是你也可以很容易地把它变成一个简单的矩阵运算。m从 的每个其他元素中减去 的每个元素m就是m - m.T。（您可以使乘积更明确，但numpy处理添加行向量和列向量的方式不是必需的。）然后您只需将其中的每个元素与 scalar 进行比较a，然后 reduce with any，就完成了：

>>> a = -2
>>> m = np.matrix((1, 3, 10, 5))
>>> subtractions = m - m.T
>>> subtractions
matrix([[ 0,  2,  9,  4],
        [-2,  0,  7,  2],
        [-9, -7,  0, -5],
        [-4, -2,  5,  0]])
>>> comparisons = subtractions < a
>>> comparisons
matrix([[False, False, False, False],
        [False, False, False, False],
        [ True,  True, False,  True],
        [ True, False, False, False]], dtype=bool)
>>> np.any(comparisons, 0)
matrix([[ True,  True, False,  True]], dtype=bool)

或者，将它们放在一行中：

>>> np.any((m - m.T) < a, 0)
matrix([[ True,  True,  True,  True]], dtype=bool)

如果您需要m成为数组而不是矩阵，则可以将减法线替换为m - np.matrix(m).T.

对于更高维度，您实际上确实需要在数组中工作，因为您正在尝试将 2D 数组与自身进行笛卡尔积以获得 4D 数组，并且numpy不做 4D 矩阵。因此，您不能使用简单的“行向量 - 列向量 = 矩阵”技巧。但是您可以手动执行此操作：

>>> m = np.array([[1,2], [3,4]]) # 2x2
>>> m4d = m.reshape(1, 1, 2, 2)  # 1x1x2x2
>>> m4d
array([[[[1, 2],
         [3, 4]]]])
>>> mt4d = m4d.T # 2x2x1x1
>>> mt4d
array([[[[1]],
        [[3]]],
       [[[2]],
        [[4]]]])
>>> subtractions = m - mt4d # 2x2x2x2
>>> subtractions
array([[[[ 0,  1],
         [ 2,  3]],
        [[-2, -1],
         [ 0,  1]]],
       [[[-1,  0],
         [ 1,  2]],
        [[-3, -2],
         [-1,  0]]]])

从那里开始，其余部分与以前相同。把它放在一行中：

>>> np.any((m - m.reshape(1, 1, 2, 2).T) < a, 0)

（如果您还记得我最初的答案，我会不知何故reshape通过乘以m1 的列向量来做同样的事情，这显然是一种更愚蠢的方法。）

最后一个快速的想法：如果您的算法真的是“（对于 , 的任何元素）对于每个元素的布尔结果，那么y您实际上不需要“对于任何元素”，您可以只使用“对于最大元素” . 所以你可以从 O(N^2) 简化为 O(N)：mx - y < axmyy

>>> (m - m.max()) < a

或者，如果a是肯定的，那总是错误的，所以你可以简化为 O(1)：

>>> np.zeros(m.shape, dtype=bool)

但我猜你真正的算法实际上是在使用abs(x - y)，或者更复杂的东西，不能以这种方式简化。

于 2012-11-29T23:28:47.450 回答