1

The title might be ambiguous, didn't know how else to word it.

I have gotten a bit far with my particle simulator in python using numpy and matplotlib, I have managed to implement coloumb, gravity and wind, now I just want to add temperature and pressure but I have a pre-optimization question (root of all evil). I want to see when particles crash:

Q: Is it in numpy possible to take the difference of an array with each of its own element based on a bool condition? I want to avoid looping.

Eg: (x - any element in x) < a Should return something like

[True, True, False, True]

If element 0,1 and 3 in x meets the condition.

Edit:

The loop quivalent would be:

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

I notice numpy operations are far faster than if tests and this has helped me speed up things ALOT.

Here is a sample code if anyone wants to play with it.

#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()

I plan to add colors to the datapoints above once I figure out how...such that high velocity particles can easily be distinguished in larger boxes.

3
  • you might want to use scipy.spatial.cKDTree instead... though I guess for a 1-D case its maybe not worth it. Commented Nov 30, 2012 at 1:12
  • @seberg: Actually, it sounds like his actual work is 2D arrays; he just gave the 1D case in the question to keep things simple. But I don't think this is relevant to him—he wants to compare each value against every other value, not against neighbors. Commented Nov 30, 2012 at 18:33
  • @abarnert I don't know, but if he uses those distances to find all those within (or not within) a certain range, then KDTree is exactly made for that (ok, you can find N-nearest neighbors too). Commented Nov 30, 2012 at 20:56

1 Answer 1

6

Eg: x - any element in x < a Should return something like

[True, True, False, True]

If element 0,1 and 3 in x meets the condition. I notice numpy operations are far faster than if tests and this has helped me speed up things ALOT.

Yes, it's just m < a. For example:

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

Now, to the question:

Q: Is it in numpy possible to take the difference of an array with each of its own element based on a bool condition? I want to avoid looping.

I'm not sure what you're asking for here, but it doesn't seem to match the example directly below it. Are you trying to, e.g., subtract 1 from each element that satisfies the predicate? In that case, you can rely on the fact that False==0 and True==1 and just subtract the boolean array:

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

From your clarification, you want the equivalent of this pseudocode loop:

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

I think the confusion here is that this is the exact opposite of what you said: you don't want "the difference of an array with each of its own element based on a bool condition", but "a bool condition based on the difference of an array with each of its own elements". And even that only really gets you to a square matrix of len(m)*len(m) bools, but I think the part left over is that the "any".

At any rate, you're asking for an implicit cartesian product, comparing each element of m to each element of m.

You can easily reduce this from two loops to one (or, rather, implicitly vectorize one of them, gaining the usual numpy performance benefits). For each value, create a new array by subtracting that value from each element and comparing the result with a, and then join those up:

>>> 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)

But you can also turn this into a simple matrix operation pretty easily. Subtracting every element of m from every other element of m is just m - m.T. (You can make the product more explicit, but the way numpy handles adding row and column vectors, it isn't necessary.) And then you just compare every element of that to the scalar a, and reduce with any, and you're done:

>>> 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)

Or, putting it all together in one line:

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

If you need m to be an array rather than a matrix, you can replace the subtraction line with m - np.matrix(m).T.

For higher dimensions, you actually do need to work in arrays, because you're trying to cartesian-product a 2D array with itself to get a 4D array, and numpy doesn't do 4D matrices. So, you can't use the simple "row vector - column vector = matrix" trick. But you can do it manually:

>>> 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]]]])

And from there, the remainder is the same as before. Putting it together into one line:

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

(If you remember my original answer, I'd somehow blanked on reshape and was doing the same thing by multiplying m by a column vector of 1s, which obviously is a much stupider way to proceed.)

One last quick thought: If your algorithm really is "the bool result of (for any element y of m, x - y < a) for each element x of m", you don't actually need "for any element y", you can just use "for the maximal element y". So you can simplify from O(N^2) to O(N):

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

Or, if a is positive, that's always false, so you can simplify to O(1):

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

But I'm guessing your real algorithm is actually using abs(x - y), or something more complicated, which can't be simplified in this way.

Sign up to request clarification or add additional context in comments.

6 Comments

Essentially I want to take the distance of an array with each of its elements and check if the result fits <a. Please see edit :)
@user948652: OK, I think I know what you want now. See the new version.
Clever way! Works well for 1D arrays, the problem is if I have [x,y,z] in each row then the shapes won't agree but at least this cuts a normal python loop. Anyway I learned a lot and will try and apply this, hopefully I can find a workaround to disable both loops.
The same idea extends to higher dimensions… except that the shortcut trick of "row vector - column vector = 2D matrix" doesn't work; you have to manually raise the dimension of the array before transposing so you can subtract across to get the higher-D array.
@user948652: Reshape it from N dimensions to 2N dimensions, with the extras having a size of 1, then transpose, then subtract that. See my updated answer for details.
|

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.