[Numpy-discussion] Trick for fast

[Sebastian Berg]

I guess Einsum is much cleaner, but I already had started with this and
maybe someone likes it, this is fully vectorized and uses a bit of funny
stuff too:

# The dot product(s), written using broadcasting rules:
a = -(x.reshape(-1,1,3) * x[...,None])

# Magic, to avoid the eye thing, takes all diagonal elements as view,
maybe there is a cooler way for it:
diagonals = np.lib.stride_tricks.as_strided(a, (a.shape[0], 3), (a.dtype.itemsize*9, a.dtype.itemsize*4))

# Add the x**2 (s is a view on the diagonals), the sum is broadcasted.
diagonals += (sum(x**2, 1))[:,None]

# And multiply by mass using broadcasting:
a *= mass[...,None]

# And sum up all the intermediat results:
inert = a.sum(0)

print inert

Regards,

Sebastian

On Fri, 2012-02-03 at 15:43 +0100, Søren Gammelmark wrote:
> What about this?
> 
> 
> A = einsum("i,ij->", mass, x ** 2)
> B = einsum("i,ij,ik->jk", mass, x, x)
> I = A * eye(3) - B
> 
> 
> /Søren
> 
> On 3 February 2012 15:10, <josef.pktd at gmail.com> wrote:
>         On Fri, Feb 3, 2012 at 8:44 AM, Alan G Isaac
>         <alan.isaac at 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 at scipy.org
>         > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>         _______________________________________________
>         NumPy-Discussion mailing list
>         NumPy-Discussion at scipy.org
>         http://mail.scipy.org/mailman/listinfo/numpy-discussion
>         
> 
> 
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion