On Mon, Feb 2, 2015 at 3:40 PM, Paulo van Breugel <[email protected]> wrote:
> > > On Mon, Feb 2, 2015 at 1:38 PM, Glynn Clements <[email protected]> > wrote: > >> >> Paulo van Breugel wrote: >> >> > > The main speed-up will come from "inlining" distance.mahalanobis(), >> > > which is essentially: >> > > >> > > delta = u - v >> > > m = np.dot(np.dot(delta, VI), delta) >> > > return np.sqrt(m) >> > > >> > > np.dot(v, m) is equivalent to np.sum(v * m,axis=0), and np.dot(u,v) is >> > > equivalent to np.sum(u * v), so the second line above is equivalent to >> > > >> > > m = np.sum(np.sum(delta[:,None] * VI * delta[None,:], axis=-1), >> > > axis=-1) >> > >> > >> > Smart! I can follow the logic, but I am not sure how to solve the >> problem >> > below: >> > >> > Traceback (most recent call last): >> > File "<stdin>", line 1, in <module> >> > ValueError: operands could not be broadcast together with shapes >> > (3,1,77,78) (3,3) >> > >> > Which refers to the different dimensions of delta and VI? >> >> The first version of the code (which is quoted above) will only work >> with 1-D arrays. It's just a fairly direct translation of the existing >> distance.mahalanobis() function, given to explain how it can then be >> extened to 3-D arrays (i.e. a 2-D array of 1-D vectors). >> >> The second version: >> >> > > VI = np.linalg.inv(covar) >> > > delta = dat_ref - stat_mean[:,None,None] >> > > m = np.sum(np.sum(delta[:,:,:,None] * delta[:,:,None,:] * >> VI[None,None,:,:],axis=-1),axis=-1) >> > > stat_mah = np.sqrt(m) >> >> should work with delta being a 3-D array. >> > > > Yes, but when running the lines above, I am getting the same error message: > > Traceback (most recent call last): > File "<stdin>", line 1, in <module> > ValueError: operands could not be broadcast together with shapes > (3,77,78,78) (1,1,3,3) > > 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? > > delta[:,:,:,None].shape > (3, 77, 78, 1) > > delta[:,:,None,:].shape > (3, 77, 1, 78) > > VI[None,None,:,:].shape > (1, 1, 3, 3) > > 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. > > m = np.sum(np.sum(delta[None,:,:,:] * VI[:,:,None,None],axis=1) * delta, > axis=0) > > > This does, in my limited testing, gives the same result as using the loop > with > > pow(distance.mahalanobis(cell_ref, stat_mean, VI),2). > > > to be sure, does the above makes sense to you? > 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 :-) > > Paulo > > > > > > > > > > >> -- >> Glynn Clements <[email protected]> >> > >
_______________________________________________ grass-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-dev
