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
