[Numpy-discussion] Trick for fast

[josef.pktd@gmai...]

On Fri, Feb 3, 2012 at 8:44 AM, Alan G Isaac <alan.isaac@gmail.com> wrote:
> On 2/3/2012 5:16 AM, santhu kumar wrote:
>> x = nX3 vector.
>> mass = nX1 vector
>> inert = zeros((3,3))
>> for i in range(n):
>>        ri = x[i,:].reshape(1,3)
>>        inert = inert + mass[i,]*(sum(ri*ri)*eye(3) - dot(ri.T,ri))
>>
>
>
> This should buy you a bit.
>
> xdot = (x*x).sum(axis=1)
> for (massi,xi,xdoti) in zip(mass.flat,x,xdot):
>       temp = -np.outer(xi,xi)
>       temp.flat[slice(0,None,4)] += xdoti
>       inert += massi*temp
>
> Alan Isaac

maybe something like this, (self contained example and name spaces to
make running it easier)

import numpy as np
n = 15
x = np.arange(n*3.).reshape(-1,3) #nX3 vector.
mass = np.linspace(1,2,n)[:,None] #nX1 vector
inert = np.zeros((3,3))
for i in range(n):
      ri = x[i,:].reshape(1,3)
      inert = inert + mass[i,]*(sum(ri*ri)*np.eye(3) - np.dot(ri.T,ri))
print inert

print  np.diag((mass * x**2).sum(0))  - np.dot(x.T, mass*x)

[[     0.  -16755.  -17287.5]
 [-16755.       0.  -17865. ]
 [-17287.5 -17865.       0. ]]
[[     0.  -16755.  -17287.5]
 [-16755.       0.  -17865. ]
 [-17287.5 -17865.       0. ]]

Josef


>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion