# [Numpy-discussion] effectively computing variograms with numpy

Timothy Hochberg tim.hochberg at ieee.org
Fri Jun 22 14:02:57 EDT 2007

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 <klemm at phys.ethz.ch > 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 < tim.hochberg at ieee.org> 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 <klemm at phys.ethz.ch> 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
> >
> >
> >
> > --
> > .  __
> > .   |-\
> > .
> > .  tim.hochberg at ieee.org
> >
> > ------=_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:klemm at phys.ethz.ch ">klemm at 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;">
> > <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:tim.hochberg at ieee.org"> tim.hochberg at ieee.org</a>
> >
> > ------=_Part_157389_1558912.1182523880067--
> >
>
>
>
> --
> Hanno Klemm
> klemm at phys.ethz.ch
>
>
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>

--
.  __
.   |-\
.
.  tim.hochberg at ieee.org
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20070622/59bdea77/attachment.html>