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 回答