<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Feb 2, 2015 at 3:40 PM, Paulo van Breugel <span dir="ltr"><<a href="mailto:p.vanbreugel@gmail.com" target="_blank">p.vanbreugel@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote"><div><div>On Mon, Feb 2, 2015 at 1:38 PM, Glynn Clements <span dir="ltr"><<a href="mailto:glynn@gclements.plus.com" target="_blank">glynn@gclements.plus.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><span><br>
Paulo van Breugel wrote:<br>
<br>
> > The main speed-up will come from "inlining" distance.mahalanobis(),<br>
> > which is essentially:<br>
> ><br>
> > delta = u - v<br>
> > m = np.dot(np.dot(delta, VI), delta)<br>
> > return np.sqrt(m)<br>
> ><br>
> > np.dot(v, m) is equivalent to np.sum(v * m,axis=0), and np.dot(u,v) is<br>
> > equivalent to np.sum(u * v), so the second line above is equivalent to<br>
> ><br>
> > m = np.sum(np.sum(delta[:,None] * VI * delta[None,:], axis=-1),<br>
> > axis=-1)<br>
><br>
><br>
> Smart! I can follow the logic, but I am not sure how to solve the problem<br>
> below:<br>
><br>
> Traceback (most recent call last):<br>
> File "<stdin>", line 1, in <module><br>
> ValueError: operands could not be broadcast together with shapes<br>
> (3,1,77,78) (3,3)<br>
><br>
> Which refers to the different dimensions of delta and VI?<br>
<br>
</span>The first version of the code (which is quoted above) will only work<br>
with 1-D arrays. It's just a fairly direct translation of the existing<br>
distance.mahalanobis() function, given to explain how it can then be<br>
extened to 3-D arrays (i.e. a 2-D array of 1-D vectors).<br>
<br>
The second version:<br>
<span><br>
> > VI = np.linalg.inv(covar)<br>
> > delta = dat_ref - stat_mean[:,None,None]<br>
> > m = np.sum(np.sum(delta[:,:,:,None] * delta[:,:,None,:] * VI[None,None,:,:],axis=-1),axis=-1)<br>
> > stat_mah = np.sqrt(m)<br>
<br>
</span>should work with delta being a 3-D array.<br></blockquote><div><br><br></div></div></div><div>Yes, but when running the lines above, I am getting the same error message:<br><br><span style="font-family:monospace,monospace"><span>Traceback (most recent call last):<br> File "<stdin>", line 1, in <module><br></span>ValueError: operands could not be broadcast together with shapes (3,77,78,78) (1,1,3,3)</span><br><br></div><div>From what I can see from broadcasting rules, there is mismatch in the last and fore-last dimensions of VI[None,None,:,:] compared to the other two?<br><br><span style="font-family:monospace,monospace">delta[:,:,:,None].shape<br>(3, 77, 78, 1)<br><br>delta[:,:,None,:].shape<br>(3, 77, 1, 78)<br><br>VI[None,None,:,:].shape<br>(1, 1, 3, 3)</span><br><br></div><div>I am really have to get a better grasp of working with these arrays, but in any case, after a bit trial and error, I came up with the following to compute m.<br><br>
<p style="margin:0px;text-indent:0px">
</p><p style="margin:0px;text-indent:0px"><span style="font-family:monospace,monospace">m = np.sum(np.sum(delta[None,:,:,:] * VI[:,:,None,None],axis=1) * delta, axis=0)</span></p><p style="margin:0px;text-indent:0px"><br></p><p style="margin:0px;text-indent:0px">This does, in my limited testing, gives the same result as using the loop with
</p><p style="margin:0px;text-indent:0px">pow(distance.mahalanobis(cell_ref, stat_mean, VI),2).</p><p style="margin:0px;text-indent:0px"><br></p><p style="margin:0px;text-indent:0px">to be sure, does the above makes sense to you?</p></div></div></div></div></blockquote><div><br><br></div><div>About the speed up.. that is indeed rather significant. With three layers as input (cells: 1122582) my original function took 30 seconds, the one above 0.23 seconds. Not bad :-)<br><br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><span><font color="#888888"><p style="margin:0px;text-indent:0px"><br></p><p style="margin:0px;text-indent:0px">Paulo<br></p><p style="margin:0px;text-indent:0px"><br></p><p style="margin:0px;text-indent:0px"><br></p> <br><br><br><br></font></span></div><span><font color="#888888"><div><br><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<span><font color="#888888"><br>
--<br>
Glynn Clements <<a href="mailto:glynn@gclements.plus.com" target="_blank">glynn@gclements.plus.com</a>><br>
</font></span></blockquote></font></span></div><br></div></div>
</blockquote></div><br></div></div>