[SciPy-User] Get array of separation vectors from an array a vectors

[josef.pktd at gmail.com]

On Wed, Jan 13, 2010 at 7:47 PM, davide_fiocco at yahoo.it
<davide_fiocco at yahoo.it> wrote:
> Hi folks,
> I'm new to Python and I'm trying to implement a basic molecular
> dynamics code.
>
> The problem I have is the following:
> Suppose you have an array of N vectors in R^3 like:
> A = [ [x1,y1,z1], [x2,y2,z2], ..., [xN,yN,zN] ]
>
> what I need is to get N!/(2! (N-2)!) separation vectors between the
> vectors in A, i.e.
> D = [ [x1-x2,y1-y2,z1-z2], [x1-x3,y1-y3,z1-z3], ..., [x2-x3,y2-y3,z2-
> z3], ...,  [x_i-x_j,y_i-y_j,z_i-z_j],...]
>
> and I need the code to be FAST! Else I think I'll switch to a Fortran/
> F2Py implementation.
>
> I'd say this task is not too different to what
> scipy.spatial.distance.pdist() does, with the difference that i don't
> need (the euclidean, say) distance but the differences between all the
> pairs of vectors in A.
>
> All suggestions will be very welcome, and I apologize if this is too
> trivial! Thank you.


Something along the following, is the only thing I can come up with.
Still requires intermediate arrays, and I thought I saw somewhere in
numpy a function that creates the indices for a triu (but don't
remember where)


import numpy as np

n = 5 #4
a = np.arange(n*3).reshape(n,3)
print a

#full
ind0, ind1 = np.mgrid[0:n,0:n]
ind0, ind1 = ind0.ravel(), ind1.ravel()
d = a[ind1,:]-a[ind0,:]
print d

#reduced
triuind0, triuind1 = np.nonzero(np.triu(np.ones((n,n)),k=1))
dr = a[triuind0,:]-a[triuind1,:]
print dr

'''
>>> import scipy
>>> scipy.comb(4,2,exact=1)
6L
>>> scipy.comb(5,2,exact=1)
10L
'''

Warning quickly written and untested.

>>> a
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11],
       [12, 13, 14]])
>>> dr
array([[ -3,  -3,  -3],
       [ -6,  -6,  -6],
       [ -9,  -9,  -9],
       [-12, -12, -12],
       [ -3,  -3,  -3],
       [ -6,  -6,  -6],
       [ -9,  -9,  -9],
       [ -3,  -3,  -3],
       [ -6,  -6,  -6],
       [ -3,  -3,  -3]])

Josef

>
> Davide
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>