[Numpy-discussion] Trick for fast

[josef.pktd at gmail.com]

On Fri, Feb 3, 2012 at 2:33 PM,  <josef.pktd at gmail.com> wrote:
> On Fri, Feb 3, 2012 at 1:58 PM, santhu kumar <mesanthu at gmail.com> wrote:
>> Hi Josef,
>>
>> I am unclear on what you want to say, but all I am doing in the code is
>> getting inertia tensor for a bunch of particle masses.
>> (http://en.wikipedia.org/wiki/Moment_of_inertia#Moment_of_inertia_tensor)
>>
>> So the diagonals are not actually zeros but would have z^2 + y^2 ..
>> The reason which I said 3secs could be misunderstood .. This code is called
>> many times over a loop and the bulk time is taken in computing this inertial
>> tensor.
>> After code change, the entire loop finishes off in 3 ses.
>
> ok, I still had a python sum instead of the numpy sum in there  (I
> really really like namespaces :)
>
>>>> np.sum(ri*ri)
> 5549.0
>>>> sum(ri*ri)
> array([ 1764.,  1849.,  1936.])
>
> this might match better.
>
> print  (mass * x**2).sum() * np.eye(3)  - np.dot(x.T, mass*x)
>
> or
> res = - np.dot(x.T, mass*x)
> res[np.arange(3), np.arange(3)] += (mass * x**2).sum()
> print res

if my interpretation is now correct, and because I have seen many traces lately

>>> res = - np.dot(x.T, mass*x)
>>> res[np.arange(3), np.arange(3)] -= np.trace(res)
>>> res
array([[ 35752.5, -16755. , -17287.5],
       [-16755. ,  34665. , -17865. ],
       [-17287.5, -17865. ,  33532.5]])


>>> res - inert
array([[  7.27595761e-12,   0.00000000e+00,   3.63797881e-12],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00]])


Josef


>
> Josef
>>
>> Thanks for alertness,
>> Santhosh
>>
>> On Fri, Feb 3, 2012 at 12:47 PM, <numpy-discussion-request at scipy.org> wrote:
>>>
>>> Send NumPy-Discussion mailing list submissions to
>>>        numpy-discussion at scipy.org
>>>
>>> To subscribe or unsubscribe via the World Wide Web, visit
>>>        http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>> or, via email, send a message with subject or body 'help' to
>>>        numpy-discussion-request at scipy.org
>>>
>>> You can reach the person managing the list at
>>>        numpy-discussion-owner at scipy.org
>>>
>>> When replying, please edit your Subject line so it is more specific
>>> than "Re: Contents of NumPy-Discussion digest..."
>>>
>>>
>>> Today's Topics:
>>>
>>>   1. Re: Trick for fast (santhu kumar)
>>>   2. Re: Trick for fast (josef.pktd at gmail.com)
>>>
>>>
>>> ----------------------------------------------------------------------
>>>
>>> Message: 1
>>> Date: Fri, 3 Feb 2012 12:29:26 -0600
>>> From: santhu kumar <mesanthu at gmail.com>
>>>
>>> Subject: Re: [Numpy-discussion] Trick for fast
>>> To: numpy-discussion at scipy.org
>>> Message-ID:
>>>
>>>  <CA+7TRsszzxPLAOz_9pgU5e9mYcHTDQZMzRoMPzpfhrRkcStZEQ at mail.gmail.com>
>>> Content-Type: text/plain; charset="iso-8859-1"
>>>
>>>
>>> On Fri, Feb 3, 2012 at 1:29 PM, santhu kumar <mesanthu at gmail.com> wrote:
>>> > Hello all,
>>> >
>>> > Thanks for lovely solutions. I have sat on it for some time and wrote it
>>> > myself :
>>> >
>>> > n =x.shape[0]
>>> > ea = np.array([1,0,0,0,1,0,0,0,1])
>>> > inert = ((np.tile(ea,(n,1))*((x*x).sum(axis=1)[:,np.newaxis]) -
>>> >
>>> > np.hstack([x*x[:,0][:,np.newaxis],x*x[:,1][:,np.newaxis],x*x[:,2][:,np.newaxis]]))*mass[:,np.newaxis]).sum(axis=0)
>>> > inert.shape = 3,3
>>> >
>>> > Does the trick and reduces the time from over 45 secs to 3 secs.
>>> > I do want to try einsum but my numpy is little old and it does not have
>>> > it.
>>> >
>>> > Thanks Sebastian (it was tricky to understand your code for me) and
>>> > Josef
>>> > (clean).
>>>
>>> Isn't the entire substraction of the first term just to set the
>>> diagonal of the result to zero.
>>>
>>> It looks to me now just like the weighted dot product and setting the
>>> diagonal to zero. That shouldn't take 3 secs unless you actual
>>> dimensions are huge.
>>>
>>> Josef
>>>
>>>
>>>
>>>
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>