0

My aim is to create an nxn array of the angle between one particle and every other particle in a simulation, for all particles. Then a masked array can select all particles in a particular field of view.

My question is how to take two nxn arrays (dx/distance and dy/distance) row by row and take the dot product with a row of a nx2 vector to result in a nx1 array of angles between one particle and every other particle. The repeat this for all rows and result in a nxn array of all the angles.

Context - There are n particles of position (x, y) and velocity (x, y). The offset vector between each particle can be calculated by creating an n x n array dx and n x n array dy. The offset vector (from particle i to particle j) is (xi - xj, yi - yj), which we can get from (dx, dy). Create unit vector component dx/distance and dy/distance.

n = 4
k = 10
width = 50
boid_radius = 8
dim = 2

position = np.random.rand(n, dim) * width  # random dataset
velocity = 0.5 * np.random.random_sample((n, dim)) + 1

from sklearn import preprocessing as pp
velocity_normalized = pp.normalize(velocity)

dx = np.subtract.outer(position[:, 0], position[:, 0])
dy = np.subtract.outer(position[:, 1], position[:, 1])
distance = np.hypot(dx, dy)
# mask zeros
ox = dx/distance
ox = dy/distance

Example Data:

position
Out[233]: 
array([[  6.68625116,  34.35642605],
       [ 18.96766714,  45.61291941],
       [ 49.49921981,  37.95450382],
       [ 28.22272906,  42.90652135]])

dx
Out[234]: 
array([[  0.        , -12.28141597, -42.81296865, -21.5364779 ],
       [ 12.28141597,   0.        , -30.53155268,  -9.25506192],
       [ 42.81296865,  30.53155268,   0.        ,  21.27649075],
       [ 21.5364779 ,   9.25506192, -21.27649075,   0.        ]])

dy
Out[235]: 
array([[  0.        , -11.25649336,  -3.59807777,  -8.5500953 ],
       [ 11.25649336,   0.        ,   7.65841559,   2.70639806],
       [  3.59807777,  -7.65841559,   0.        ,  -4.95201753],
       [  8.5500953 ,  -2.70639806,   4.95201753,   0.        ]])

Form the offset unit vector components:

distance = np.hypot(dx, dy).round()
array([[  0.,  17.,  43.,  23.],
       [ 17.,   0.,  31.,  10.],
       [ 43.,  31.,   0.,  22.],
       [ 23.,  10.,  22.,   0.]])

zeros = ma.masked_where(distance==0, distance)
masked_array(data =
 [[-- 17.0 43.0 23.0]
 [17.0 -- 31.0 10.0]
 [43.0 31.0 -- 22.0]
 [23.0 10.0 22.0 --]],
             mask =
 [[ True False False False]
 [False  True False False]
 [False False  True False]
 [False False False  True]],
       fill_value = 1e+20)

ox = dx / zeros
masked_array(data =
 [[-- -0.7224362336046727 -0.9956504337130146 -0.9363686041759272]
 [0.7224362336046727 -- -0.9848887960767806 -0.9255061924766889]
 [0.9956504337130146 0.9848887960767806 -- 0.9671132160733322]
 [0.9363686041759272 0.9255061924766889 -0.9671132160733322 --]],
             mask =
 [[ True False False False]
 [False  True False False]
 [False False  True False]
 [False False False  True]],
       fill_value = 1e+20)

The next step is to take the dot product between the unit velocity vector and the unit offset vector. That's where I'm stuck.

velocity
Out[239]: 
array([[ 1.08980931,  1.10142992],
       [ 1.42378512,  1.31445528],
       [ 1.4599431 ,  1.34567875],
       [ 1.45934809,  1.03997269]])

pp.normalize(velocity)
Out[242]: 
array([[ 0.70334695,  0.71084672],
       [ 0.73475404,  0.67833363],
       [ 0.73529551,  0.67774665],
       [ 0.81437139,  0.58034407]])

Next step? I can't work out how to index the vectors to get it all done in the way that numpy makes it look like magic.

Blue Shrapnel
  • 95
  • 1
  • 13
  • What is the question? – MB-F Nov 23 '17 at 07:46
  • I would not use sklearn's normalize function for this. If I understand correctly all you need is `dx/distance` and `dy/distance` to get unit directions. – MB-F Nov 23 '17 at 08:02
  • Good point, thank you. My question though is about converting this into the numpythonic approach where vectors are treated as a whole instead of with for loops to iterate through them. – Blue Shrapnel Nov 23 '17 at 08:37
  • What for loops? `dx/distance` for example computes all normalized x directions without any loop in Python. – MB-F Nov 23 '17 at 08:52
  • Thank you yes, I've got that step. It's computing the dot product of the velocity vector and the offset vector between each particle and every other particle, resulting in a nxn array of angles between all particles and every other particle. – Blue Shrapnel Nov 23 '17 at 09:25
  • I see. Can you show what the expected result would be in the example? You will most likely need something like `np.matmul(np.stack([dx, dy]).T, vel[..., None])[..., 0]`, or `np.einsum('kij,jk->ji', np.stack([dx, dy]), vel)`, or simply `dx*vel[:,0] + dy*vel[:,1]`. – MB-F Nov 23 '17 at 11:53

1 Answers1

0

Thanks @kazemakase, answer staring me in the face... trying to look for the over complicated. 

import numpy as np
import sklearn.preprocessing as pp
import numpy.ma as ma

n = 4 # number of particles
dim = 2 # x and y coords

pos = np.random.randint(0,9,(n,2))
vel = np.random.randint(0,9,(n,2))
vel_unit = pp.normalize(vel)

dx = np.subtract.outer(pos[:,0], pos[:,0])
dy = np.subtract.outer(pos[:,1], pos[:,1])
distance = np.hypot(dx, dy).round()
m = ma.masked_where(distance==0, distance)
o_unit_x = dx/m
o_unit_y = dy/m
angles = o_unit_y * vel_unit[:,1] + o_unit_x * vel_unit[:,0]

Just need to test that's right and it's doing what I think it is.

Sample Output:

pos
Out[154]: 
array([[6, 1],
       [5, 7],
       [0, 7],
       [5, 2]])

vel
Out[155]: 
array([[6, 3],
       [7, 5],
       [5, 4],
       [2, 2]])

distance
Out[156]: 
array([[ 0.,  6.,  8.,  1.],
       [ 6.,  0.,  5.,  5.],
       [ 8.,  5.,  0.,  7.],
       [ 1.,  5.,  7.,  0.]])

angles.round(3)
Out[158]: 
masked_array(data =
 [[-- -0.446 0.117 0.0]
 [0.298 -- 0.781 0.707]
 [-0.335 -0.814 -- 0.0]
 [-0.447 -0.581 0.112 --]],
             mask =
 [[ True False False False]
 [False  True False False]
 [False False  True False]
 [False False False  True]],
       fill_value = 1e+20)
Blue Shrapnel
  • 95
  • 1
  • 13