Tim, Thank you very much, the code does what's it expected to do. Unfortunately the thing is still pretty slow on large data sets. I will probably now look for ways to calculate the variogram from some random samples of my data. Thanks for the observation regarding the square array, that would have bitten me, later.
Again, thanks, Hanno Timothy Hochberg <[EMAIL PROTECTED]> said: > ------=_Part_160363_13970610.1182535377481 > Content-Type: text/plain; charset=ISO-8859-1; format=flowed > Content-Transfer-Encoding: 7bit > Content-Disposition: inline > > OK, generally in code like this I leave the outer loops alone and try to > vectorize just the inner loop.I have some ideas in this direction, but > first, there seems to be some problems with the code at well. The code looks > like it is written to take non-square 'data' arrays. However, > > for i in range(data.shape[0]): > datasquared = (data - data[:,i])**2 > > This is looping over shape[0], but indexing on axis-1, which doesn't work > for non-square arrays. One suggestion is make a function to compute the > variogram along a given axis and then calling it twice instead of computing > them both independently. Can you try the following code and see if this > correctly implements a variogram? I don't have time to check that it really > implements a variogram, but I'm hoping it's close: > > def variogram(data, binsize, axis=-1): > data = data.swapaxes(-1, axis) > n = data.shape[-1] > resultsize = int(N.ceil(n / float(binsize))) > result = N.zeros([resultsize], data.dtype) > for i in range(resultsize): > j0 = max(i*binsize, 1) > j1 = min(j0+binsize, n) > denominator = 0 > for j in range(j0, j1): > d2 = (data[...,j:] - data[...,:-j])**2 > result[i] += d2.sum() > denominator += N.prod(d2.shape) > result[i] /= denominator > return result > > > > > On 6/22/07, Hanno Klemm <[EMAIL PROTECTED] > wrote: > > > > Tim, > > > > this is the best I could come up with until now: > > > > > > import numpy as N > > > > def naive_variogram(data, binsize=100., stepsize=5.): > > """calculates variograms along the rows and columns of the given > > array which is supposed to contain equally spaced data with > > stepsize stepsize""" > > > > # how many elements do fit in one bin? > > > > binlength = int(binsize/stepsize) > > > > #bins in x- and y- direction (+1 for the possible > > #elements larger than int(binsize/stepsize): > > x_bins = (data.shape[1])/binlength+1 > > y_bins = (data.shape[0])/binlength+1 > > > > #arrays to store the reuslts in > > x_result = N.zeros(x_bins, dtype = float) > > y_result = N.zeros(y_bins, dtype = float) > > > > #arrays to get teh denominators right > > x_denominators = N.zeros(x_bins, dtype=float) > > y_denominators = N.zeros(x_bins, dtype=float) > > > > #what is the last index? > > xlast = data.shape[1] > > ylast = data.shape[0] > > for i in range(data.shape[0]): > > datasquared = (data - data[:,i])**2 > > #number of bins to fill until the end of the array: > > numbins = 1 + (xlast - i)/binlength > > for j in range(numbins): > > x_result[j]+=\ > > datasquared[:,i+1+j*binlength:i+1+(j+1)*binlength].sum() > > x_denominators[j] +=\ > > datasquared[:,i+1+j*binlength:i+1+(j+1)*binlength].size > > try: > > #Is there a rest? > > x_result[numbins] += > > datasquared[:,i+1+numbins*binlength:].sum() > > x_denominators[numbins] += > > datasquared[:,i+1+numbins*binlength:].size > > except IndexError: > > pass > > > > x_result /= x_denominators > > > > > > for i in range(data.shape[1]): > > datasquared = (data - data[i])**2 > > #number of bins to fill until the end of the array: > > numbins = 1 + (ylast - i)/binlength > > #Fill the bins > > for j in range(numbins): > > > > y_result[j]+=datasquared[i+1+j*binlength:i+1+(j+1)*binlength].sum() > > y_denominators[j] += > > datasquared[i+1+j*binlength:i+1+(j+1)*binlength].size > > try: > > #Is there a rest? > > y_result[numbins] += > > datasquared[:,i+1+numbins*binlength:].sum() > > y_denominators[numbins] += > > datasquared[:,i+1+numbins*binlength:].size > > except IndexError: > > pass > > > > y_result /= y_denominators > > > > return x_result, y_result > > > > Thanks, > > Hanno > > > > > > Timothy Hochberg < [EMAIL PROTECTED]> said: > > > > > ------=_Part_157389_1558912.1182523880067 > > > Content-Type: text/plain; charset=ISO-8859-1; format=flowed > > > Content-Transfer-Encoding: 7bit > > > Content-Disposition: inline > > > > > > On 6/22/07, Hanno Klemm <[EMAIL PROTECTED]> wrote: > > > > > > > > > > > > Hi, > > > > > > > > I have an array which represents regularly spaced spatial data. I now > > > > would like to compute the (semi-)variogram, i.e. > > > > > > > > gamma(h) = 1/N(h) \sum_{i,j\in N(h)} (z_i - z_j)**2, > > > > > > > > where h is the (approximate) spatial difference between the > > > > measurements z_i, and z_j, and N(h) is the number of measurements with > > > > distance h. > > > > > > > > However, I only want to calculate the thing along the rows and > > > > columns. The naive approach involves two for loops and a lot of > > > > searching, which becomes painfully slow on large data sets. Are there > > > > better implementations around in numpy/scipy or does anyone have a > > > > good idea of how to do that more efficient? I looked around a bit but > > > > couldn't find anything. > > > > > > > > > Can you send the naive code as well. Its often easier to see what's > > going on > > > with code in addition to the equations. > > > > > > Regards. > > > > > > -tim > > > > > > > > > > > > -- > > > . __ > > > . |-\ > > > . > > > . [EMAIL PROTECTED] > > > > > > ------=_Part_157389_1558912.1182523880067 > > > Content-Type: text/html; charset=ISO-8859-1 > > > Content-Transfer-Encoding: 7bit > > > Content-Disposition: inline > > > > > > <br><br><div><span class="gmail_quote">On 6/22/07, <b > > class="gmail_sendername">Hanno Klemm</b> <<a > > href="mailto:[EMAIL PROTECTED] ">[EMAIL PROTECTED] </a>> > > wrote:</span><blockquote class="gmail_quote" style="border-left: 1px > > solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"> > > > <br>Hi,<br><br>I have an array which represents regularly spaced > > spatial data. I now<br>would like to compute the (semi-)variogram, > > i.e.<br><br>gamma(h) = 1/N(h) \sum_{i,j\in N(h)} (z_i - > > z_j)**2,<br><br>where h is the (approximate) spatial difference > > between the > > > <br>measurements z_i, and z_j, and N(h) is the number of > > measurements with<br>distance h.<br><br>However, I only want to > > calculate the thing along the rows and<br>columns. The naive approach > > involves two for loops and a lot of > > > <br>searching, which becomes painfully slow on large data sets. Are > > there<br>better implementations around in numpy/scipy or does anyone > > have a<br>good idea of how to do that more efficient? I looked around > > a bit but<br>couldn't find anything. > > > </blockquote><div><br>Can you send the naive code as well. Its often > > easier to see what's going on with code in addition to the > > equations.<br><br>Regards.<br><br>-tim<br><br></div></div><br > > clear="all"><br>-- <br>. __ > > > <br>. |-\<br>.<br>. <a > > href="mailto:[EMAIL PROTECTED]"> [EMAIL PROTECTED]</a> > > > > > > ------=_Part_157389_1558912.1182523880067-- > > > > > > > > > > > -- > > Hanno Klemm > > [EMAIL PROTECTED] > > > > > > _______________________________________________ > > Numpy-discussion mailing list > > Numpy-discussion@scipy.org > > http://projects.scipy.org/mailman/listinfo/numpy-discussion > > > > > > -- > . __ > . |-\ > . > . [EMAIL PROTECTED] > > ------=_Part_160363_13970610.1182535377481 > Content-Type: text/html; charset=ISO-8859-1 > Content-Transfer-Encoding: 7bit > Content-Disposition: inline > > <br>OK, generally in code like this I leave the outer loops alone and try to vectorize just the inner loop.I have some ideas in this direction, but first, there seems to be some problems with the code at well. The code looks like it is written to take non-square 'data' arrays. However, > <br><br> for i in range(data.shape[0]):<br> datasquared = (data - data[:,i])**2<br><div style="margin-left: 40px;"><br></div>This is looping over shape[0], but indexing on axis-1, which doesn't work for non-square arrays. One suggestion is make a function to compute the variogram along a given axis and then calling it twice instead of computing them both independently. Can you try the following code and see if this correctly implements a variogram? I don't have time to check that it really implements a variogram, but I'm hoping it's close: > <br><br> def variogram(data, binsize, axis=-1):<br> data = data.swapaxes(-1, axis)<br> n = data.shape[-1]<br> resultsize = int(N.ceil(n / float(binsize)))<br> result = N.zeros([resultsize], > data.dtype)<br> for i in range(resultsize):<br> j0 = max(i*binsize, 1)<br> j1 = min(j0+binsize, n)<br> denominator = 0<br> for j in range(j0, j1):<br> d2 = (data[...,j:] - data[...,:-j])**2 > <br> result[i] += d2.sum()<br> denominator += N.prod(d2.shape)<br> result[i] /= denominator<br> return result<br><br><br><br><br><div><span class="gmail_quote">On 6/22/07, <b class="gmail_sendername"> > Hanno Klemm</b> <<a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">[EMAIL PROTECTED] > </a>> wrote:</span><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"> > Tim,<br><br>this is the best I could come up with until now:<br><br><br>import numpy as N<br><br>def naive_variogram(data, binsize=100., stepsize=5.):<br> """calculates variograms along the rows and columns of the given > <br> array which is supposed to contain equally spaced data with<br> stepsize stepsize"""<br><br> # how many elements do fit in one bin?<br><br> binlength = int(binsize/stepsize)<br><br> #bins in x- and y- direction (+1 for the possible > <br> #elements larger than int(binsize/stepsize):<br> x_bins = (data.shape[1])/binlength+1<br> y_bins = (data.shape[0])/binlength+1<br><br> #arrays to store the reuslts in<br> x_result = N.zeros(x_bins, dtype = float) > <br> y_result = N.zeros(y_bins, dtype = float)<br><br> #arrays to get teh denominators right<br> x_denominators = N.zeros(x_bins, dtype=float)<br> y_denominators = N.zeros(x_bins, dtype=float)<br><br> #what is the last index? > <br> xlast = data.shape[1]<br> ylast = data.shape[0]<br> for i in range(data.shape[0]):<br> datasquared = (data - data[:,i])**2<br> #number of bins to fill until the end of the array:<br> numbins = 1 + (xlast - i)/binlength > <br> for j in range(numbins):<br> x_result[j]+=\<br>datasquared[:,i+1+j*binlength:i+1+(j+1)*binlength].sum()<br> x_denominators[j] +=\<br>datasquared[:,i+1+j*binlength:i+1+(j+1)*binlength].size > <br> try:<br> #Is there a rest?<br> x_result[numbins] +=<br>datasquared[:,i+1+numbins*binlength:].sum()<br> x_denominators[numbins] +=<br>datasquared[:,i+1+numbins*binlength:].size > <br> > except IndexError:<br> pass<br><br> x_result /= x_denominators<br><br><br> for i in range(data.shape[1]):<br> datasquared = (data - data[i])**2<br> #number of bins to fill until the end of the array: > <br> numbins = 1 + (ylast - i)/binlength<br> #Fill the bins<br> for j in range(numbins):<br><br>y_result[j]+=datasquared[i+1+j*binlength:i+1+(j+1)*binlength].sum()<br> y_denominators[j] += > <br> > datasquared[i+1+j*binlength:i+1+(j+1)*binlength].size<br> try:<br> #Is there a rest?<br> y_result[numbins] +=<br>datasquared[:,i+1+numbins*binlength:].sum()<br> y_denominators[numbins] += > <br>datasquared[:,i+1+numbins*binlength:].size<br> except IndexError:<br> pass<br><br> y_result /= y_denominators<br><br> return x_result, y_result<br><br>Thanks,<br>Hanno<br><br><br>Timothy Hochberg < > <a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">[EMAIL PROTECTED]</a>> said:<br><br>> ------=_Part_157389_1558912.1182523880067<br>> Content-Type: text/plain; charset=ISO-8859-1; format=flowed > <br>> Content-Transfer-Encoding: 7bit > <br>> Content-Disposition: inline<br>><br>> On 6/22/07, Hanno Klemm <<a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">[EMAIL PROTECTED]</a>> wrote: > <br>> ><br>> ><br>> > Hi,<br>> ><br>> > I have an array which represents regularly spaced spatial data. I now > <br>> > would like to compute the (semi-)variogram, i.e.<br>> ><br>> > gamma(h) = 1/N(h) \sum_{i,j\in N(h)} (z_i - z_j)**2,<br>> ><br>> > where h is the (approximate) spatial difference between the > <br>> > measurements z_i, and z_j, and N(h) is the number of measurements with<br>> > distance h.<br>> ><br>> > However, I only want to calculate the thing along the rows and<br>> > columns. The naive approach involves two for loops and a lot of > <br>> > searching, which becomes painfully slow on large data sets. Are there<br>> > better implementations around in numpy/scipy or does anyone have a<br>> > good idea of how to do that more efficient? I looked around a bit but > <br>> > couldn't find anything.<br>><br>><br>> Can you send the naive code as well. Its often easier to see what's<br>going on<br>> with code in addition to the equations.<br>><br>> Regards. > <br>><br>> -tim<br>><br>><br>><br>> --<br>> __<br>> . |-\<br>> .<br>> <a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">[EMAIL PROTECTED] > </a><br>><br>> ------=_Part_157389_1558912.1182523880067 > <br>> Content-Type: text/html; charset=ISO-8859-1<br>> Content-Transfer-Encoding: 7bit<br>> Content-Disposition: inline<br>><br>> <br><br><div><span class="gmail_quote">On 6/22/07, <b > <br>class="gmail_sendername">Hanno Klemm</b> &lt;<a<br>href="mailto:<a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">[EMAIL PROTECTED] > </a>"><a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">[EMAIL PROTECTED] > </a></a>&gt;<br>wrote:</span><blockquote class="gmail_quote" style="border-left: 1px<br>solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"><br>> <br>Hi,<br><br>I have an array which represents regularly spaced > <br>spatial data. I now<br>would like to compute the (semi-)variogram,<br>i.e.<br><br>gamma(h) = 1/N(h) \sum_{i,j\in N(h)} (z_i -<br>z_j)**2,<br><br>where h is the (approximate) spatial difference > <br>between the<br>> <br>measurements z_i, and z_j, and N(h) is the number of<br>measurements with<br>distance h.<br><br>However, I only want to<br>calculate the thing along the rows and<br>columns. The naive approach > <br>involves two for loops and a lot of<br>> <br>searching, which becomes painfully slow on large data sets. Are<br>there<br>better implementations around in numpy/scipy or does anyone<br>have a<br>good idea of how to do that more efficient? I looked around > <br>a bit but<br>couldn&#39;t find anything.<br>> </blockquote><div><br>Can you send the naive code as well. Its often<br>easier to see what&#39;s going on with code in addition to the<br> > > > equations.<br><br>Regards.<br><br>-tim<br><br></div></div><br<br>clear="all"><br>-- <br>.&nbsp;&nbsp;__<br>> <br>.&nbsp;&nbsp; |-\<br>.<br>.&nbsp;&nbsp;<a > <br>href="mailto:<a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">[EMAIL PROTECTED]</a>"><a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)"> > > [EMAIL PROTECTED]</a></a><br>><br>> ------=_Part_157389_1558912.1182523880067-- > <br>><br><br><br><br>--<br>Hanno Klemm<br><a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">[EMAIL PROTECTED]</a><br><br><br>_______________________________________________ > <br>Numpy-discussion mailing list<br><a href="mailto:Numpy-discussion@scipy.org" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)"> > Numpy-discussion@scipy.org</a><br><a href="http://projects.scipy.org/mailman/listinfo/numpy-discussion" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">http://projects.scipy.org/mailman/listinfo/numpy-discussion > </a><br></blockquote></div><br><br clear="all"><br>-- <br> > . __<br>. |-\<br>.<br>. <a href="mailto:[EMAIL PROTECTED]" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">[EMAIL PROTECTED]</a> > > ------=_Part_160363_13970610.1182535377481-- > -- Hanno Klemm [EMAIL PROTECTED] _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion