<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):<br>           datasquared = (data - data[:,i])**2<br><div style="margin-left: 40px;"><br></div>This is looping over shape, 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:klemm@phys.ethz.ch" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">klemm@phys.ethz.ch
</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)/binlength+1<br>    y_bins = (data.shape)/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<br>    ylast = data.shape<br>    for i in range(data.shape):<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):<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:tim.hochberg@ieee.org" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">tim.hochberg@ieee.org</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:klemm@phys.ethz.ch" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">klemm@phys.ethz.ch</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:tim.hochberg@ieee.org" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">tim.hochberg@ieee.org
</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:klemm@phys.ethz.ch" target="_blank" onclick="return top.js.OpenExtLink(window,event,this)">klemm@phys.ethz.ch
</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