numpy-discussion — Discussion list for all users of Numerical Python

RE: [Numpy-discussion] numarray.where confusion

[Perry G. <pe...@st...>]

Francesc Alted va escriure:

> A Dimecres 26 Maig 2004 21:01, Perry Greenfield va escriure:
> > correct. You'd have to break apart the m1 tuple and
> > index all the components, e.g.,
> > 
> > m11, m12 = m1
> > x[m11[m2],m12[m2]] = ...
> > 
> > This gets clumsier with the more dimensions that must
> > be handled, but you still can do it. It would be most
> > useful if the indexed array is very large, the number
> > of items selected is relatively small and one
> > doesn't want to incur the memory overhead of all the
> > mask arrays of the admittedly much nicer notational
> > approach that Francesc illustrated.
> 
> Well, boolean arrays have the property that they use very little memory
> (only 1 byte / element), and normally perform quite well doing indexing.
> Some timings:
> 
> >>> import timeit
> >>> t1 = 
> timeit.Timer("m1=where(x>4);m2=where(x[m1]<7);m11,m12=m1;x[m11[m2]
> ,m12[m2]]","from numarray import 
> arange,where;dim=3;x=arange(dim*dim);x.shape=(dim,dim)")
> >>> t2 = timeit.Timer("x[(x>4) & (x<7)]","from numarray import 
> arange,where;dim=3;x=arange(dim*dim);x.shape=(dim,dim)")
> >>> t1.repeat(3,1000)
> [3.1320240497589111, 3.1235389709472656, 3.1198310852050781]
> >>> t2.repeat(3,1000)
> [1.1218469142913818, 1.117638111114502, 1.1156759262084961]
> 
> i.e. using boolean arrays for indexing is roughly 3 times faster.
> 
> For larger arrays this difference is even more noticeable:
> 
> >>> t3 = 
> timeit.Timer("m1=where(x>4);m2=where(x[m1]<7);m11,m12=m1;x[m11[m2]
> ,m12[m2]]","from numarray import 
> arange,where;dim=1000;x=arange(dim*dim);x.shape=(dim,dim)")
> >>> t4 = timeit.Timer("x[(x>4) & (x<7)]","from numarray import 
> arange,where;dim=1000;x=arange(dim*dim);x.shape=(dim,dim)")
> >>> t3.repeat(3,10)
> [3.1818649768829346, 3.20477294921875, 3.190640926361084]
> >>> t4.repeat(3,10)
> [0.42328095436096191, 0.42140507698059082, 0.41979002952575684]
> 
> as you see, now the difference is almost an order of magnitude (!).
> 
> So, perhaps assuming the small memory overhead, in most of cases it is
> better to use boolean selections. However, it would be nice to know the
> ultimate reason of why this happens, because the Perry approach seems
> intuitively faster.
>
Yes I agree. It was good of you to post these timings. I don't 
think we had actually compared the two approaches though the
results don't surprise me (though I suspect the results may change
if the first mask has a very small percentage of elements; the
large timing test has nearly all elements selected for the first
mask).

Perry

Re: [Numpy-discussion] numarray.where confusion

[Francesc A. <fa...@py...>]

A Dimecres 26 Maig 2004 21:01, Perry Greenfield va escriure:
> correct. You'd have to break apart the m1 tuple and
> index all the components, e.g.,
> 
> m11, m12 = m1
> x[m11[m2],m12[m2]] = ...
> 
> This gets clumsier with the more dimensions that must
> be handled, but you still can do it. It would be most
> useful if the indexed array is very large, the number
> of items selected is relatively small and one
> doesn't want to incur the memory overhead of all the
> mask arrays of the admittedly much nicer notational
> approach that Francesc illustrated.

Well, boolean arrays have the property that they use very little memory
(only 1 byte / element), and normally perform quite well doing indexing.
Some timings:

>>> import timeit
>>> t1 = timeit.Timer("m1=where(x>4);m2=where(x[m1]<7);m11,m12=m1;x[m11[m2],m12[m2]]","from numarray import arange,where;dim=3;x=arange(dim*dim);x.shape=(dim,dim)")
>>> t2 = timeit.Timer("x[(x>4) & (x<7)]","from numarray import arange,where;dim=3;x=arange(dim*dim);x.shape=(dim,dim)")
>>> t1.repeat(3,1000)
[3.1320240497589111, 3.1235389709472656, 3.1198310852050781]
>>> t2.repeat(3,1000)
[1.1218469142913818, 1.117638111114502, 1.1156759262084961]

i.e. using boolean arrays for indexing is roughly 3 times faster.

For larger arrays this difference is even more noticeable:

>>> t3 = timeit.Timer("m1=where(x>4);m2=where(x[m1]<7);m11,m12=m1;x[m11[m2],m12[m2]]","from numarray import arange,where;dim=1000;x=arange(dim*dim);x.shape=(dim,dim)")
>>> t4 = timeit.Timer("x[(x>4) & (x<7)]","from numarray import arange,where;dim=1000;x=arange(dim*dim);x.shape=(dim,dim)")
>>> t3.repeat(3,10)
[3.1818649768829346, 3.20477294921875, 3.190640926361084]
>>> t4.repeat(3,10)
[0.42328095436096191, 0.42140507698059082, 0.41979002952575684]

as you see, now the difference is almost an order of magnitude (!).

So, perhaps assuming the small memory overhead, in most of cases it is
better to use boolean selections. However, it would be nice to know the
ultimate reason of why this happens, because the Perry approach seems
intuitively faster.

-- 
Francesc Alted

RE: [Numpy-discussion] numarray.where confusion

[Perry G. <pe...@st...>]

> > try 
> > 
> > x[m1[0][m2]] = array([10,20])
> > 
> > instead. The intent here is to provide x with the net index array
> > by indexing m1 first rather than indexing x first.
> > (note the odd use of m1[0]; this is necessary since where() will
> > return a tuple of index arrays (to allow use in multidimensional 
> > cases as indices, so the m1[0] extracts the array from the tuple;
> > Since m1 is a tuple, indexing it with another index array (well, 
> > tuple containing an index array) doesn't work).
> 
> This works, but for the fact that in my real code I *am* dealing with
> multidimensional arrays.  But this is a nice trick to remember.
> 
> (So, the following "does not work":
> 
>     x = arange(9)
>     x.shape=(3,3)
>     m1 = where(x > 4)
>     m2 = where(x[m1] < 7)
>     x[m1[0][m2]]
> )
 
correct. You'd have to break apart the m1 tuple and
index all the components, e.g.,

m11, m12 = m1
x[m11[m2],m12[m2]] = ...

This gets clumsier with the more dimensions that must
be handled, but you still can do it. It would be most
useful if the indexed array is very large, the number
of items selected is relatively small and one
doesn't want to incur the memory overhead of all the
mask arrays of the admittedly much nicer notational
approach that Francesc illustrated.

Perry
 

Re: [Numpy-discussion] numarray.where confusion

[Alok S. <as...@vi...>]

On 26/05/04: 10:43, Andrew Straw wrote:
> Todd Miller wrote:
> >On Wed, 2004-05-26 at 12:06, Francesc Alted wrote:
> >
> >>>>>a = arange(10)
> >>>>>a[(a>5) & (a<8)] = array([10, 20])
> >>>>>
> Is there an equivalently slick way to accomplish to what I'm trying 
> below?  (the the values in c[:,1] get changed based on the same-row 
> values in c[:,0]?)
> 
> from numarray import *
> a=arange(10)
> b=arange(10)+20
> c=concatenate((a[:,NewAxis],b[:,NewAxis]),axis=1)
> c[c[:,0]>7][:,1] = 0 # doesn't work because it makes a copy and 
> therefore doesn't modify c

Well, for your case, the following works:

>>> print c
[[ 0 20]
 [ 1 21]
 [ 2 22]
 [ 3 23]
 [ 4 24]
 [ 5 25]
 [ 6 26]
 [ 7 27]
 [ 8 28]
 [ 9 29]]
>>> t0 = c[:, 0]
>>> t1 = c[:, 1]
>>> t1[t0 > 7] = 0
>>> print c
[[ 0 20]
 [ 1 21]
 [ 2 22]
 [ 3 23]
 [ 4 24]
 [ 5 25]
 [ 6 26]
 [ 7 27]
 [ 8  0]
 [ 9  0]]

Not sure this helps in your real code though.

Alok

-- 
Alok Singhal (as...@vi...)          *   *          
Graduate Student, dept. of Astronomy   *           *     *
University of Virginia                                    
http://www.astro.virginia.edu/~as8ca/              *    * 

Re: [Numpy-discussion] numarray.where confusion

[Andrew S. <str...@as...>]

Todd Miller wrote:

>On Wed, 2004-05-26 at 12:06, Francesc Alted wrote:
>  
>
>>A Dimecres 26 Maig 2004 17:41, Todd Miller va escriure:
>>    
>>
>>>Here's how I did it (there was an easier way I overlooked):
>>>
>>>a = arange(10)
>>>m1 = where(a > 5, 1, 0).astype('Bool')
>>>m2 = where(a < 8, 1, 0).astype('Bool')
>>>a[m1 & m2] = array([10, 20])
>>>      
>>>
>>Perhaps the easier way looks like this?
>>
>>    
>>
>>>>>a = arange(10)
>>>>>a[(a>5) & (a<8)] = array([10, 20])
>>>>>          
>>>>>
Is there an equivalently slick way to accomplish to what I'm trying 
below?  (the the values in c[:,1] get changed based on the same-row 
values in c[:,0]?)

from numarray import *
a=arange(10)
b=arange(10)+20
c=concatenate((a[:,NewAxis],b[:,NewAxis]),axis=1)
c[c[:,0]>7][:,1] = 0 # doesn't work because it makes a copy and 
therefore doesn't modify c

Cheers!
Andrew

Re: [Numpy-discussion] numarray.where confusion

[Alok S. <as...@vi...>]

On 26/05/04: 11:24, Perry Greenfield wrote:
> (due to confusions with "a" in text I'll use x in place of "a")
> I believe the problem you are seeing (I'm not 100% certain yet)
> is that although it is possible to assign to an array-indexed
> array, that doing that twice over doesn't work since Python is,
> in effect, treating x[m1] as an expression even though it is
> on the left side. That expression results in a new array that the
> second indexing updates, but then is thrown away since it is not
> assigned to anything else.
> 
> Your second try creates a temporary t which is also not a view into
> a so when you update t, a is not updated.

Thanks or this info.  It makes sense now.  I suspected earlier that t
was not a view but a copy, but didn't realise that the same thing was
happening with x[m1][m2].

> try 
> 
> x[m1[0][m2]] = array([10,20])
> 
> instead. The intent here is to provide x with the net index array
> by indexing m1 first rather than indexing x first.
> (note the odd use of m1[0]; this is necessary since where() will
> return a tuple of index arrays (to allow use in multidimensional 
> cases as indices, so the m1[0] extracts the array from the tuple;
> Since m1 is a tuple, indexing it with another index array (well, 
> tuple containing an index array) doesn't work).

This works, but for the fact that in my real code I *am* dealing with
multidimensional arrays.  But this is a nice trick to remember.

(So, the following "does not work":

    x = arange(9)
    x.shape=(3,3)
    m1 = where(x > 4)
    m2 = where(x[m1] < 7)
    x[m1[0][m2]]
)

On 26/05/04: 11:41, Todd Miller wrote:
> Here's how I did it (there was an easier way I overlooked):
> 
> a = arange(10)
> m1 = where(a > 5, 1, 0).astype('Bool')
> m2 = where(a < 8, 1, 0).astype('Bool')
> a[m1 & m2] = array([10, 20])

Ah.  This works!  Even for multidimensional arrays.

On 26/05/04: 18:06, Francesc Alted wrote:
> Perhaps the easier way looks like this?
> 
> >>> a = arange(10)
> >>> a[(a>5) & (a<8)] = array([10, 20])
> >>> a
> array([ 0,  1,  2,  3,  4,  5, 10, 20,  8,  9])
> 
> Indexing is a very powerful (and fun) thing, indeed :)

I like this too.  Thank you all for the help!

Alok

-- 
Alok Singhal (as...@vi...)       __
Graduate Student, dept. of Astronomy   /  _
University of Virginia                 \_O \
http://www.astro.virginia.edu/~as8ca/    __/

Re: [Numpy-discussion] numarray.where confusion

[Todd M. <jm...@st...>]

On Wed, 2004-05-26 at 12:06, Francesc Alted wrote:
> A Dimecres 26 Maig 2004 17:41, Todd Miller va escriure:
> > Here's how I did it (there was an easier way I overlooked):
> > 
> > a = arange(10)
> > m1 = where(a > 5, 1, 0).astype('Bool')
> > m2 = where(a < 8, 1, 0).astype('Bool')
> > a[m1 & m2] = array([10, 20])
> 
> Perhaps the easier way looks like this?
> 
> >>> a = arange(10)
> >>> a[(a>5) & (a<8)] = array([10, 20])

Much, much better.  Thanks!

Todd

> >>> a
> array([ 0,  1,  2,  3,  4,  5, 10, 20,  8,  9])
> 
> Indexing is a very powerful (and fun) thing, indeed :)
-- 
Todd Miller <jm...@st...>

Re: [Numpy-discussion] numarray.where confusion

[Francesc A. <fa...@py...>]

A Dimecres 26 Maig 2004 17:41, Todd Miller va escriure:
> Here's how I did it (there was an easier way I overlooked):
> 
> a = arange(10)
> m1 = where(a > 5, 1, 0).astype('Bool')
> m2 = where(a < 8, 1, 0).astype('Bool')
> a[m1 & m2] = array([10, 20])

Perhaps the easier way looks like this?

>>> a = arange(10)
>>> a[(a>5) & (a<8)] = array([10, 20])
>>> a
array([ 0,  1,  2,  3,  4,  5, 10, 20,  8,  9])

Indexing is a very powerful (and fun) thing, indeed :)

-- 
Francesc Alted

Re: [Numpy-discussion] numarray.where confusion

[Todd M. <jm...@st...>]

On Wed, 2004-05-26 at 10:48, Alok Singhal wrote:
> Hi,
> 
> I am having trouble understanding how exactly "where" works in
> numarray.
> 
> What I am trying to do:
> 
> I am preparing a two-level mask in an array and then assign values to
> the array where both masks are true:
> 
> >>> from numarray import *
> >>> a = arange(10)
> >>> # First mask
> >>> m1 = where(a > 5)
> >>> a[m1]
> array([6, 7, 8, 9])
> >>> # Second mask
> >>> m2 = where(a[m1] < 8)
> >>> a[m1][m2]

a[m1] is a new array here.

> array([6, 7])
> >>> # So far so good
> >>> # Now change some values
> >>> a[m1][m2] = array([10, 20])

And here too.  This does a write into what is effectively a temporary
variable returned by the expression a[m1].  Although the write occurs,
it is lost.

> >>> a[m1][m2]
> array([6, 7])
> >>> a
> array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

Here's how I did it (there was an easier way I overlooked):

a = arange(10)
m1 = where(a > 5, 1, 0).astype('Bool')
m2 = where(a < 8, 1, 0).astype('Bool')
a[m1 & m2] = array([10, 20])

The principle here is to keep the masks as "full sized" boolean arrays
rather than index arrays so they can be combined using the bitwise and
operator.   The resulting mask can be used to index just once
eliminating the temporary.

Regards,
Todd

RE: [Numpy-discussion] numarray.where confusion

[Perry G. <pe...@st...>]

Alok Singhal wrote:
> Hi,
> 
> I am having trouble understanding how exactly "where" works in
> numarray.
> 
> What I am trying to do:
> 
> I am preparing a two-level mask in an array and then assign values to
> the array where both masks are true:
> 
> >>> from numarray import *
> >>> a = arange(10)
> >>> # First mask
> >>> m1 = where(a > 5)
> >>> a[m1]
> array([6, 7, 8, 9])
> >>> # Second mask
> >>> m2 = where(a[m1] < 8)
> >>> a[m1][m2]
> array([6, 7])
> >>> # So far so good
> >>> # Now change some values
> >>> a[m1][m2] = array([10, 20])
> >>> a[m1][m2]
> array([6, 7])
> >>> a
> array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
> >>> # Didn't  work
> >>> # Let's try a temporary variable
> >>> t = a[m1]
> >>> t[m2]
> array([6, 7])
> >>> t[m2] = array([10, 20])
> >>> t[m2], t
> (array([10, 20]), array([10, 20,  8,  9]))
> >>> a
> array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
> 
> So, my assignment to a[m1][m2] seems to work (no messages), but it
> doesn't produce the effect I want it to.
> 
> I have read the documentation but I couldn't find something that would
> explain this behavior.
> 
> So my questions:
> 
> - did I miss something important in the documentation,
> - I am expecting something I shouldn't, or
> - there is a bug in numarray?
> 
(due to confusions with "a" in text I'll use x in place of "a")
I believe the problem you are seeing (I'm not 100% certain yet)
is that although it is possible to assign to an array-indexed
array, that doing that twice over doesn't work since Python is,
in effect, treating x[m1] as an expression even though it is
on the left side. That expression results in a new array that the
second indexing updates, but then is thrown away since it is not
assigned to anything else.

Your second try creates a temporary t which is also not a view into
a so when you update t, a is not updated.

try 

x[m1[0][m2]] = array([10,20])

instead. The intent here is to provide x with the net index array
by indexing m1 first rather than indexing x first.
(note the odd use of m1[0]; this is necessary since where() will
return a tuple of index arrays (to allow use in multidimensional 
cases as indices, so the m1[0] extracts the array from the tuple;
Since m1 is a tuple, indexing it with another index array (well, 
tuple containing an index array) doesn't work).

Perry Greenfield

[Numpy-discussion] numarray.where confusion

[Alok S. <as...@vi...>]

Hi,

I am having trouble understanding how exactly "where" works in
numarray.

What I am trying to do:

I am preparing a two-level mask in an array and then assign values to
the array where both masks are true:

>>> from numarray import *
>>> a = arange(10)
>>> # First mask
>>> m1 = where(a > 5)
>>> a[m1]
array([6, 7, 8, 9])
>>> # Second mask
>>> m2 = where(a[m1] < 8)
>>> a[m1][m2]
array([6, 7])
>>> # So far so good
>>> # Now change some values
>>> a[m1][m2] = array([10, 20])
>>> a[m1][m2]
array([6, 7])
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> # Didn't  work
>>> # Let's try a temporary variable
>>> t = a[m1]
>>> t[m2]
array([6, 7])
>>> t[m2] = array([10, 20])
>>> t[m2], t
(array([10, 20]), array([10, 20,  8,  9]))
>>> a
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

So, my assignment to a[m1][m2] seems to work (no messages), but it
doesn't produce the effect I want it to.

I have read the documentation but I couldn't find something that would
explain this behavior.

So my questions:

- did I miss something important in the documentation,
- I am expecting something I shouldn't, or
- there is a bug in numarray?

Thanks,
Alok

-- 
Alok Singhal (as...@vi...)       __
Graduate Student, dept. of Astronomy   /  _
University of Virginia                 \_O \
http://www.astro.virginia.edu/~as8ca/    __/

RE: [Numpy-discussion] image correlations

[Perry G. <pe...@st...>]

John Hunter writes:
> >>>>> "Todd" == Todd Miller <jm...@st...> writes:
> 
>     >> I have a series of luminance images that I want to do some
>     >> correlation analyses on.  Each image is an MxN frame of a
>     >> movie.  I want to compute the correlation between a given pixel
>     >> i,j, and every other pixel in the image over each frame.  That
>     >> is, if I assume xij is a numFrames length time series, and xkl
>     >> is another numFrames length time series (the pixel intensities
>     >> at two points in the image over time), I want to compute the
>     
>     >> corrcoeff(xij, xkl) for every kl with ij fixed.
>     
>     >> I know I could do this by looping over the pixels of the image,
>     >> but I'm hoping for something a bit faster.
> 
>     Todd> For numarray try numarray.convolve.correlate2d and set
>     Todd> fft=1.
> 
Something that Todd and I have talked about is making general
functions (like fft, convolve, correlate) broadcastable to
other dimensions. Not much has been done to do that but I think 
a general mechanism could be developed to do just that. In that
way one could do a 1-d correlation over a 2 or 3-d array without
looping over the extra dimensions (all that would be done implicitly
in C) with the option of specifying which dimension the function
should be applied to while looping over all the other dimensions.

In the meantime, it would seemt that one possible way of dealing
with this is to simply recast your array as a 1-d array to use
the 1-d correlation. This means some memory copying and explicit
padding (how do you want the correlation to handle points beyond
the ends?). Supposing you want it to wraparound (in effect a
circular correlation) you could do something like this (untested
and probably has some mistakes :-)

if imagecube is a T x M x N array where there are T time samples.

data = concatenate((imagecube, imagecube, imagecube))
    # pad the time series on both ends with identical copies
reference = imagecube[:, i, j]
flatdata = ravel(transpose(data))
flatresult = correlate(flatdata, reference)
result = transpose(reshape(flatresult, (N,M,T)))[T:2*T,:,:]

This is admittedly a bit clumsy, but should be much, much faster
than iterating over all pixels.

It would be much nicer just to have

result = correlate(imagecube, imagecube[:,i,j], axis=0)

which the broadcasting facility would permit.



Perry

Re: [Numpy-discussion] image correlations

[John H. <jdh...@ac...>]

>>>>> "Todd" == Todd Miller <jm...@st...> writes:

    >> I have a series of luminance images that I want to do some
    >> correlation analyses on.  Each image is an MxN frame of a
    >> movie.  I want to compute the correlation between a given pixel
    >> i,j, and every other pixel in the image over each frame.  That
    >> is, if I assume xij is a numFrames length time series, and xkl
    >> is another numFrames length time series (the pixel intensities
    >> at two points in the image over time), I want to compute the
    
    >> corrcoeff(xij, xkl) for every kl with ij fixed.
    
    >> I know I could do this by looping over the pixels of the image,
    >> but I'm hoping for something a bit faster.

    Todd> For numarray try numarray.convolve.correlate2d and set
    Todd> fft=1.

I've looked at this and don't know if I'm missing something or if this
isn't the right function for me.  I want to correlate a given pixel
intensity with the intensity at all other pixels over a series of
images.   Is this possible with correlate2d?

JDH

[Numpy-discussion] numarray equivalent to sign() function

[Mike Z. <zi...@uc...>]

Hi, I am in the process of converting some code from Numeric to numarray,
and it seems that numarray no longer has the sign() function -- is that
so?

ex:

sign(-30.0) = -1
sign(0) = 0
sign(1000) = 1

Is there a replacement?

Thanks,

Mike

Re: [Numpy-discussion] Read only array?

[Todd M. <jm...@st...>]

On Fri, 2004-05-21 at 23:05, Shin, Daehyok wrote:
> I am wondering if there is any way to set an array as read-only,
> so that any trial to modify values of elements in the array raises some
> warning.
> Is it a cool feature to prevent unexpected side effect?
> 
> Daehyok Shin (Peter)
> 

In numarray you can do this:

>>> from numarray import *
>>> def make_readonly(a):
..     v = a.view()
..     s = a.tostring()
..     v._data = s
..     return v
..
>>> a = arange(100)
>>> b = make_readonly(a)
>>> b[0] = 1
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
ValueError: NA_setFromPythonScalar: assigment to readonly array buffer

This works because of the buffer protocol and the fact that a string is
a read-only buffer.  Using the buffer protocol is a numarray feature so
I'm pretty sure Numeric can't do this.  Also note that in C, this
read-only array is actually writable.

Regards,
Todd

[Numpy-discussion] Read only array?

[Shin, D. <sd...@em...>]

I am wondering if there is any way to set an array as read-only,
so that any trial to modify values of elements in the array raises some
warning.
Is it a cool feature to prevent unexpected side effect?

Daehyok Shin (Peter)

Re: [Numpy-discussion] image correlations

[Todd M. <jm...@st...>]

On Fri, 2004-05-21 at 15:01, John Hunter wrote:
> I have a series of luminance images that I want to do some correlation
> analyses on.  Each image is an MxN frame of a movie.  I want to
> compute the correlation between a given pixel i,j, and every other
> pixel in the image over each frame.  That is, if I assume
> xij is a numFrames length time series, and xkl is another numFrames
> length time series (the pixel intensities at two points in the image
> over time), I want to compute the 
> 
> corrcoeff(xij, xkl) for every kl with ij fixed.
> 
> I know I could do this by looping over the pixels of the image, but
> I'm hoping for something a bit faster.
> 
> Any suggestions?

For numarray try numarray.convolve.correlate2d and set fft=1.

Todd

[Numpy-discussion] image correlations

[John H. <jdh...@ac...>]

I have a series of luminance images that I want to do some correlation
analyses on.  Each image is an MxN frame of a movie.  I want to
compute the correlation between a given pixel i,j, and every other
pixel in the image over each frame.  That is, if I assume
xij is a numFrames length time series, and xkl is another numFrames
length time series (the pixel intensities at two points in the image
over time), I want to compute the 

corrcoeff(xij, xkl) for every kl with ij fixed.

I know I could do this by looping over the pixels of the image, but
I'm hoping for something a bit faster.

Any suggestions?

John Hunter

Re: [Numpy-discussion] Performance issues with numarray.searchsorted()

[Francesc A. <fa...@py...>]

A Dijous 20 Maig 2004 19:13, vareu escriure:
> The first case agrees with your results of ~10x difference in favor of
> Numeric at 10**6 elements.  The last case shows a ~3x numarray advantage
> given large numbers of values.
> 
> My analysis is this:  since searchsorted runs in O(log2(N)) time,  even
> with 10**6 elements there are only 20 iterations or so.  This is just
> not enough to overcome numarray's Python ufunc overhead.  I'll see what
> I can do,  but I think we're up against the standard numarray
> performance wall for small arrays.	

I see. What about creating a special case for when the item to be searched
is an scalar (and not an array)? In such a case, it would be theoretically
possible to get at least the bisect.bisect() performance.

Besides, if this special case would be implemented, users would be able to
call repeatedly this scalar version of searchsorted() in order to deal with
very small "values" arrays (len()<5) and still having a gain.

Just a thought,

-- 
Francesc Alted

Re: [Numpy-discussion] Performance issues with numarray.searchsorted()

[Todd M. <jm...@st...>]

On Thu, 2004-05-20 at 05:26, Robert Kern wrote:
> Francesc Alted wrote:
> > Hi,
> > 
> > I'm willing to use a lot the searchsorted function in numarray, but I'm a
> > bit surprised about the poor performance of it. Look at that:
> > 
> > 
> >>>>from time import time
> >>>>import numarray
> >>>>import Numeric
> >>>>na=numarray.arange(1000*1000)
> >>>>nu=Numeric.arange(1000*1000)
> >>>>t1=time();numarray.searchsorted(na,200*1000);time()-t1
> > 
> > 200000
> > 0.00055098533630371094
> > 
> >>>>t1=time();Numeric.searchsorted(nu,200*1000);time()-t1
> > 
> > 200000
> > 7.7962875366210938e-05
> > 
> > It may seem that Numeric is better optimised, but the standard python module
> > bisect is even faster than numarray.searchsorted:
> > 
> > 
> >>>>import bisect
> >>>>t1=time();bisect.bisect_left(na,200*1000);time()-t1
> > 
> > 200000
> > 8.8930130004882812e-05
> > 
> >>>>t1=time();bisect.bisect_left(nu,200*1000);time()-t1
> > 
> > 200000
> > 8.6069107055664062e-05
> 
> A better timing (IMHO), but with similar conclusions:
> 
> Python 2.3 (#1, Sep 13 2003, 00:49:11)
> [GCC 3.3 20030304 (Apple Computer, Inc. build 1495)] on darwin
> Type "help", "copyright", "credits" or "license" for more information.
>  >>> import timeit
>  >>> t1 = timeit.Timer("Numeric.searchsorted(a,200000)",
> "import Numeric;a=Numeric.arange(1000000)")
>  >>> t2 = timeit.Timer("numarray.searchsorted(a,200000)",
> "import numarray;a=numarray.arange(1000000)")
>  >>> t3 = timeit.Timer("bisect.bisect_left(a,200000)",
> "import Numeric;import bisect;a=Numeric.arange(1000000)")
>  >>> t4 = timeit.Timer("bisect.bisect_left(a,200000)",
> "import numarray;import bisect;a=numarray.arange(1000000)")
>  >>> t1.repeat(3,10000)
> [0.15758609771728516, 0.17469501495361328, 0.15456986427307129]
>  >>> t2.repeat(3,10000)
> [6.7581729888916016, 6.9644770622253418, 6.6776731014251709]
>  >>> t3.repeat(3,10000)
> [0.41335701942443848, 0.45698308944702148, 0.39665889739990234]
>  >>> t4.repeat(3,10000)
> [0.49930000305175781, 0.48063492774963379, 0.52067780494689941]
> 
> [Apologies for the linewraps.]
> 
> I also get similar results with double arrays. Weird.
> 
> Python 2.3 on Mac OS X 10.3.sumthin', latest CVS checkout of numarray, 
> Numeric 23.1.
Here's what I got with the numarray benchmarks after adding an extra
case:

benchmark                    i    numarray (usec)    Numeric (usec)  numarray:Numeric

searchsorted(a,p/5)        0.0               446                12              36.7
                           1.0               438                11              36.8
                           2.0               450                12              35.0
                           3.0               459                14              32.6
                           4.0               511                25              19.7
                           5.0               636                56              11.2
                           6.0               653                64              10.1

searchsorted(a,a)          0.0               285                 5              48.0
                           1.0               283                 7              39.6
                           2.0               291                29               9.8
                           3.0               368               308               1.2
                           4.0              1771              4120               0.4
                           5.0             17335             52127               0.3
                           6.0            201277            605787               0.3

p = 10**i
a = arange(p)

The first case agrees with your results of ~10x difference in favor of
Numeric at 10**6 elements.  The last case shows a ~3x numarray advantage
given large numbers of values.

My analysis is this:  since searchsorted runs in O(log2(N)) time,  even
with 10**6 elements there are only 20 iterations or so.  This is just
not enough to overcome numarray's Python ufunc overhead.  I'll see what
I can do,  but I think we're up against the standard numarray
performance wall for small arrays.	

Regards,
Todd

Re: [Numpy-discussion] Performance issues with numarray.searchsorted()

[Robert K. <rk...@uc...>]

Francesc Alted wrote:
> Hi,
> 
> I'm willing to use a lot the searchsorted function in numarray, but I'm a
> bit surprised about the poor performance of it. Look at that:
> 
> 
>>>>from time import time
>>>>import numarray
>>>>import Numeric
>>>>na=numarray.arange(1000*1000)
>>>>nu=Numeric.arange(1000*1000)
>>>>t1=time();numarray.searchsorted(na,200*1000);time()-t1
> 
> 200000
> 0.00055098533630371094
> 
>>>>t1=time();Numeric.searchsorted(nu,200*1000);time()-t1
> 
> 200000
> 7.7962875366210938e-05
> 
> It may seem that Numeric is better optimised, but the standard python module
> bisect is even faster than numarray.searchsorted:
> 
> 
>>>>import bisect
>>>>t1=time();bisect.bisect_left(na,200*1000);time()-t1
> 
> 200000
> 8.8930130004882812e-05
> 
>>>>t1=time();bisect.bisect_left(nu,200*1000);time()-t1
> 
> 200000
> 8.6069107055664062e-05

A better timing (IMHO), but with similar conclusions:

Python 2.3 (#1, Sep 13 2003, 00:49:11)
[GCC 3.3 20030304 (Apple Computer, Inc. build 1495)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
 >>> import timeit
 >>> t1 = timeit.Timer("Numeric.searchsorted(a,200000)",
"import Numeric;a=Numeric.arange(1000000)")
 >>> t2 = timeit.Timer("numarray.searchsorted(a,200000)",
"import numarray;a=numarray.arange(1000000)")
 >>> t3 = timeit.Timer("bisect.bisect_left(a,200000)",
"import Numeric;import bisect;a=Numeric.arange(1000000)")
 >>> t4 = timeit.Timer("bisect.bisect_left(a,200000)",
"import numarray;import bisect;a=numarray.arange(1000000)")
 >>> t1.repeat(3,10000)
[0.15758609771728516, 0.17469501495361328, 0.15456986427307129]
 >>> t2.repeat(3,10000)
[6.7581729888916016, 6.9644770622253418, 6.6776731014251709]
 >>> t3.repeat(3,10000)
[0.41335701942443848, 0.45698308944702148, 0.39665889739990234]
 >>> t4.repeat(3,10000)
[0.49930000305175781, 0.48063492774963379, 0.52067780494689941]

[Apologies for the linewraps.]

I also get similar results with double arrays. Weird.

Python 2.3 on Mac OS X 10.3.sumthin', latest CVS checkout of numarray, 
Numeric 23.1.

-- 
Robert Kern
rk...@uc...

"In the fields of hell where the grass grows high
  Are the graves of dreams allowed to die."
   -- Richard Harter

[Numpy-discussion] Performance issues with numarray.searchsorted()

[Francesc A. <fa...@py...>]

Hi,

I'm willing to use a lot the searchsorted function in numarray, but I'm a
bit surprised about the poor performance of it. Look at that:

>>> from time import time
>>> import numarray
>>> import Numeric
>>> na=3Dnumarray.arange(1000*1000)
>>> nu=3DNumeric.arange(1000*1000)
>>> t1=3Dtime();numarray.searchsorted(na,200*1000);time()-t1
200000
0.00055098533630371094
>>> t1=3Dtime();Numeric.searchsorted(nu,200*1000);time()-t1
200000
7.7962875366210938e-05

It may seem that Numeric is better optimised, but the standard python module
bisect is even faster than numarray.searchsorted:

>>> import bisect
>>> t1=3Dtime();bisect.bisect_left(na,200*1000);time()-t1
200000
8.8930130004882812e-05
>>> t1=3Dtime();bisect.bisect_left(nu,200*1000);time()-t1
200000
8.6069107055664062e-05

So, bisect performance is similar to that of Numeric searchsorted and both
are almost an order of magnitude better than numarray searchsorted. This a
little bit surprising as the bisect_left module is written in plain python.
=46rom the python 2.3.3 sources:

def bisect_left(a, x, lo=3D0, hi=3DNone):
    if hi is None:
        hi =3D len(a)
    while lo < hi:
        mid =3D (lo+hi)//2
        if a[mid] < x: lo =3D mid+1
        else: hi =3D mid
    return lo

I'm using python2.3.3, numarray 0.9.1 (latest CVS) and Debian Linux (sid).

Cheers,

=2D-=20
=46rancesc Alted

Re: [Numpy-discussion] random_array.poisson(0) gives -1 and -2

[Sebastian H. <ha...@ms...>]

On Wednesday 19 May 2004 12:03 pm, you wrote:
> On Wed, 2004-05-19 at 14:25, Sebastian Haase wrote:
> > Hi,
> >
> > the random_array poisson functions returns negative values if mean=0:
> > >>>from numarray import random_array as ra
> > >>> ra.seed(x=1, y=1)
> > >>> ra.poisson(0)
> >
> > 5
> >
> > >>> ra.poisson(0)
> >
> > -2
> >
> > My "math book" tells me that it should be always zero.
> > This seems to be a constructed case, but I'm using this to put
> > "quantum statistic" into a simulated image:
> > obj = na.array( something )
> > imageFromDetector = ra.poisson( obj ) + gaussianNoiseArray
> > The object array might have lots of zeros surrounding the "actual
> > object". Thinking of a fluorescent object sending out photons it makes
> > sense to not get any photons at all from 'empty' regions.
> > I'm using numarray 0.8;
>
> I tried this on Fedora-1 i386 with Python-2.3.3 and it returned zero
> consistently.  What platform are you on?
>
> Regards,
> Todd

I running debian (Woody)
$ uname -a
Linux baboon 2.4.18 #1 Tue Dec 16 14:11:01 PST 2003 i686 unknown
$python -v
<snip>
Python 2.2.1 (#1, Feb 28 2004, 00:52:10)
[GCC 2.95.4 20011002 (Debian prerelease)] on linux2

and I get this:
>>> ra.poisson(0, 100000).min()
-4
>>> ra.poisson(0, 100000).min()
-4
>>> ra.poisson(0, 100000).mean()
1.9383
>>> ra.poisson(0, 100000).mean()
1.9607
>>> ra.poisson(0, 100000).max()
29
>>> ra.poisson(0, 100000).max()
28

Thanks for checking,
Sebastian

Re: [Numpy-discussion] random_array.poisson(0) gives -1 and -2

[Todd M. <jm...@st...>]

On Wed, 2004-05-19 at 14:25, Sebastian Haase wrote:
> Hi,
> the random_array poisson functions returns negative values if mean=0:
> >>>from numarray import random_array as ra
> >>> ra.seed(x=1, y=1)
> >>> ra.poisson(0)
> 5
> >>> ra.poisson(0)
> -2
> 
> My "math book" tells me that it should be always zero.
> This seems to be a constructed case, but I'm using this to put
> "quantum statistic" into a simulated image:
> obj = na.array( something )
> imageFromDetector = ra.poisson( obj ) + gaussianNoiseArray
> The object array might have lots of zeros surrounding the "actual object". 
> Thinking of a fluorescent object sending out photons it makes sense to not get 
> any photons at all from 'empty' regions.
> I'm using numarray 0.8;

I tried this on Fedora-1 i386 with Python-2.3.3 and it returned zero
consistently.  What platform are you on?

Regards,
Todd

[Numpy-discussion] random_array.poisson(0) gives -1 and -2

[Sebastian H. <ha...@ms...>]

Hi,
the random_array poisson functions returns negative values if mean=0:
>>>from numarray import random_array as ra
>>> ra.seed(x=1, y=1)
>>> ra.poisson(0)
5
>>> ra.poisson(0)
-2

My "math book" tells me that it should be always zero.
This seems to be a constructed case, but I'm using this to put
"quantum statistic" into a simulated image:
obj = na.array( something )
imageFromDetector = ra.poisson( obj ) + gaussianNoiseArray
The object array might have lots of zeros surrounding the "actual object". 
Thinking of a fluorescent object sending out photons it makes sense to not get 
any photons at all from 'empty' regions.
I'm using numarray 0.8;

Thanks for numarray,
Sebastian Haase