Re: [Numpy-discussion] all elements equal

[Zachary Pincus]

How about the following?
exact: numpy.all(a == a[0])
inexact: numpy.allclose(a, a[0])

On Mar 5, 2012, at 2:19 PM, Keith Goodman wrote:

> On Mon, Mar 5, 2012 at 11:14 AM, Neal Becker  wrote:
>> What is a simple, efficient way to determine if all elements in an array (in 
>> my
>> case, 1D) are equal?  How about close?
> 
> For the exactly equal case, how about:
> 
> I[1] a = np.array([1,1,1,1])
> I[2] np.unique(a).size
> O[2] 1# All equal
> 
> I[3] a = np.array([1,1,1,2])
> I[4] np.unique(a).size
> O[4] 2   # All not equal
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Viewing a float64 array with a float32 array

[Zachary Pincus]

> Hi,
> 
> Is it possible to have a view of a float64 array that is itself float32?
> So that:
> 
 A = np.arange(5, dtype='d')
 A.view(dtype='f')
> 
> would return a size 5 float32 array looking at A's data?

I think not. The memory layout of a 32-bit IEEE float is not a subset of that 
of a 64-bit float -- see e.g. the first table in:
http://steve.hollasch.net/cgindex/coding/ieeefloat.html

This would work for int8/int16 or other int types (so long as the number 
doesn't exceed the range of the smaller type), but AFAIK not floats.

Note how the subset relationship works for the int8/int16 case, but not 
float32/float64:

str(numpy.array(100,dtype=numpy.int8).data)
'd'

str(numpy.array(100,dtype=numpy.int16).data)
'd\x00'

str(numpy.array(-100,dtype=numpy.int16).data)
'\x9c\xff'

str(numpy.array(-100,dtype=numpy.int8).data)
'\x9c'

str(numpy.array(100,dtype=numpy.float32).data)
'\x00\x00\xc8B'

str(numpy.array(100,dtype=numpy.float64).data)
'\x00\x00\x00\x00\x00\x00Y@'


Zach



___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Viewing a float64 array with a float32 array

[Zachary Pincus]

> I'm not sure if you are referring to rounding errors but that's OK with
> me.
> 
> I was thinking something along the lines of changing how numpy looks at
> the data of A's view by modifying say the stride attribute, etc.

Yes, so was I. As you can see in my example with ints below, you could skip 
every other byte of the int16 array to "look" at an int8 array. This is because 
the memory layout of an int8 is a proper subset of int16. (Modulo 
endian-concerns of course...)

But looking at the link I provided, you can see that taking the first 32 bits 
of an float64 (or the last 32 or any 32) does not yield something that can be 
interpreted as a float32. So there's no subset relationship, and you can't do 
the strides-trick.

To be extra clear, look at the memory layout of a float that's expressible 
without rounding error:
str(numpy.array(128,dtype=numpy.float64).data)
'\x00\x00\x00\x00\x00\x00`@'

str(numpy.array(128,dtype=numpy.float32).data)
'\x00\x00\x00C'

There's obviously no stride trick whereby one will "look" like the other. 

Zach


> 
> On Wed, Mar 21, 2012, at 11:19, Zachary Pincus wrote:
>>> Hi,
>>> 
>>> Is it possible to have a view of a float64 array that is itself float32?
>>> So that:
>>> 
>>>>>> A = np.arange(5, dtype='d')
>>>>>> A.view(dtype='f')
>>> 
>>> would return a size 5 float32 array looking at A's data?
>> 
>> I think not. The memory layout of a 32-bit IEEE float is not a subset of
>> that of a 64-bit float -- see e.g. the first table in:
>> http://steve.hollasch.net/cgindex/coding/ieeefloat.html
>> 
>> This would work for int8/int16 or other int types (so long as the number
>> doesn't exceed the range of the smaller type), but AFAIK not floats.
>> 
>> Note how the subset relationship works for the int8/int16 case, but not
>> float32/float64:
>> 
>> str(numpy.array(100,dtype=numpy.int8).data)
>> 'd'
>> 
>> str(numpy.array(100,dtype=numpy.int16).data)
>> 'd\x00'
>> 
>> str(numpy.array(-100,dtype=numpy.int16).data)
>> '\x9c\xff'
>> 
>> str(numpy.array(-100,dtype=numpy.int8).data)
>> '\x9c'
>> 
>> str(numpy.array(100,dtype=numpy.float32).data)
>> '\x00\x00\xc8B'
>> 
>> str(numpy.array(100,dtype=numpy.float64).data)
>> '\x00\x00\x00\x00\x00\x00Y@'
>> 
>> 
>> Zach
>> 
>> 
>> 
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] creating/working NumPy-ndarrays in C++

[Zachary Pincus]

> That all sounds like no option -- sad.
> Cython is no solution cause, all I want is to leave Python Syntax in
> favor for strong OOP design patterns.

What about ctypes?

For straight numerical work where sometimes all one needs to hand across the 
python-to-C/C++/Fortran boundary is a pointer to some memory and the size of 
the memory region. So I will often just write a very thin C interface to 
whatever I'm using and then call into it with ctypes.

So you could just design your algorithm in C++ with all the "strong OOP design 
patterns" you want, and then just write a minimal C interface on top with one 
or two functions that receive pointers to memory. Then allocate numpy arrays in 
python with whatever memory layout you need, and use the a "ctypes" attribute 
to find the pointer data etc. that you need to pass over. 

Zach

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Slices from an index list

[Zachary Pincus]

> Here's one way you could do it:
> 
> In [43]: indices = [0,1,2,3,5,7,8,9,10,12,13,14]
> 
> In [44]: jumps = where(diff(indices) != 1)[0] + 1
> 
> In [45]: starts = hstack((0, jumps))
> 
> In [46]: ends = hstack((jumps, len(indices)))
> 
> In [47]: slices = [slice(start, end) for start, end in zip(starts, ends)]
> 
> In [48]: slices
> Out[48]: [slice(0, 4, None), slice(4, 5, None), slice(5, 9, None), slice(9, 
> 12, None)]

If you're only going to use the slices to divide up the list, you could use 
numpy.split and skip creating the slice objects:

indices = [0,1,2,3,5,7,8,9,10,12,13,14]
jumps = numpy.where(numpy.diff(indices) != 1)[0] + 1
numpy.split(indices, jumps)

giving:
[array([0, 1, 2, 3]), array([5]), array([ 7,  8,  9, 10]), array([12, 13, 14])]

Zach

(btw, Warren, the method to calculate the jumps is cute. I'll have to remember 
that.)

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] Fancy-indexing reorders output in corner cases?

[Zachary Pincus]

Hello all,

The below seems to be a bug, but perhaps it's unavoidably part of the indexing 
mechanism?

It's easiest to show via example... note that using "[0,1]" to pull two columns 
out of the array gives the same shape as using ":2" in the simple case, but 
when there's additional slicing happening, the shapes get transposed or 
something.

In [2]: numpy.version.version # latest git version
Out[2]: '1.7.0.dev-3d4'

In [3]: d = numpy.empty((10, 9, 8, 7))

In [4]: d[:,:,:,[0,1]].shape
Out[4]: (10, 9, 8, 2)

In [5]: d[:,:,:,:2].shape
Out[5]: (10, 9, 8, 2)

In [6]: d[:,0,:,[0,1]].shape
Out[6]: (2, 10, 8)

In [7]: d[:,0,:,:2].shape
Out[7]: (10, 8, 2)

In [8]: d[0,:,:,[0,1]].shape
Out[8]: (2, 9, 8)

In [9]: d[0,:,:,:2].shape
Out[9]: (9, 8, 2)

Oddly, this error can appear/disappear depending on the position of the other 
axis sliced:
In [14]: d = numpy.empty((10, 9, 8))

In [15]: d[:,:,[0,1]].shape
Out[15]: (10, 9, 2)

In [16]: d[:,0,[0,1]].shape
Out[16]: (10, 2)

In [17]: d[0,:,[0,1]].shape
Out[17]: (2, 9)

This cannot be the expected behavior, right?
Zach

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Fancy-indexing reorders output in corner cases?

[Zachary Pincus]

> On Mon, May 14, 2012 at 4:33 PM, Zachary Pincus  
> wrote:
>> The below seems to be a bug, but perhaps it's unavoidably part of the 
>> indexing mechanism?
>> 
>> It's easiest to show via example... note that using "[0,1]" to pull two 
>> columns out of the array gives the same shape as using ":2" in the simple 
>> case, but when there's additional slicing happening, the shapes get 
>> transposed or something.
> 
> When fancy indexing and slicing is mixed, the resulting shape is
> essentially unpredictable.

Aah, right -- this does come up on the list not infrequently, doesn't it. I'd 
always thought it was more exotic usages that raised these issues. Good to know.

>  The "correct" way to do it is to only use
> fancy indexing, i.e. generate the indices of the sliced dimension as
> well.
> 

Excellent -- thanks!


> Stéfan
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Changes in PyArray_FromAny between 1.5.x and 1.6.x

[Zachary Pincus]

> It isn't just the array() calls which end up getting problems.  For
> example, in 1.5.x
> 
> sage: f = 10; type(f)
> 
> sage: numpy.arange(f)
> array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) #int64
> 
> while in 1.6.x
> 
> sage: numpy.arange(f)
> array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=object)
> 
> We also see problems with calls like
> 
> sage: scipy.stats.uniform(0,15).ppf([0.5,0.7])
> array([  7.5,  10.5])
> 
> which work in 1.5.x, but fail with a traceback "TypeError: array
> cannot be safely cast to required type" in 1.6.x.

I'm getting problems like this after a 1.6 upgrade as well. Lots of object 
arrays being created when previously there would either be an error, or an 
array of floats. 

Also, lots of the "TypeError: array cannot be safely cast to required type" are 
cropping up.

Honestly, most of these are in places where my code was lax and so I just 
cleaned things up to use the right dtypes etc. But still a bit unexpected in 
terms of having more code to fix than I was used to for 0.X numpy revisions.

Just another data-point, though. Not really a complaint.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Changes in PyArray_FromAny between 1.5.x and 1.6.x

[Zachary Pincus]

> There is a fine line here. We do need to make people clean up lax code in 
> order to improve numpy, but hopefully we can keep the cleanups reasonable.

Oh agreed. Somehow, though, I was surprised by this, even though I keep tabs on 
the numpy lists -- at no point did it become clear that "big changes in how 
arrays get constructed and typecast are ahead that may require code fixes". 
That was my main point, but probably a PEBCAK issue more than anything.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Changes in PyArray_FromAny between 1.5.x and 1.6.x

[Zachary Pincus]

>> On Tue, Jun 5, 2012 at 8:41 PM, Zachary Pincus 
>> wrote:
>>> 
>>>> There is a fine line here. We do need to make people clean up lax code
>>>> in order to improve numpy, but hopefully we can keep the cleanups
>>>> reasonable.
>>> 
>>> Oh agreed. Somehow, though, I was surprised by this, even though I keep
>>> tabs on the numpy lists -- at no point did it become clear that "big changes
>>> in how arrays get constructed and typecast are ahead that may require code
>>> fixes". That was my main point, but probably a PEBCAK issue more than
>>> anything.
>> 
>> 
>> It was fairly extensively discussed when introduced,
>> http://thread.gmane.org/gmane.comp.python.numeric.general/44206, and again
>> at some later point.
> 
> Those are the not-yet-finalized changes in 1.7; Zachary (I think) is
> talking about problems upgrading from ~1.5 to 1.6.

Yes, unless I'm wrong I experienced these problems from 1.5.something to 1.6.1. 
I didn't take notes as it was in the middle of a deadline-crunch so I just 
fixed the code and moved on (long, stupid story about why the upgrade before a 
deadline...). It's just that the issues mentioned above seem to have hit me too 
and I wanted to mention that. But unhelpfully, I think, without code, and now 
I've hijacked this thread! Sorry.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] [numpy] ENH: Initial implementation of a 'neighbor' calculation (#303)

[Zachary Pincus]

> On 10.10.2012 15:42, Nathaniel Smith wrote:
>> This PR submitted a few months ago adds a substantial new API to numpy,
>> so it'd be great to get more review. No-one's replied yet, though...
>> 
>> Any thoughts, anyone? Is it useful, could it be better...?
> 
> Fast neighbor search is what scipy.spatial.cKDTree is designed for. 
> There is an brand new version in SciPy master.

It's not neighbor *search*, it's applying a function over an (arbitrary chosen 
and weighted) moving neighborhood in an nd array.

It would be useful for the author of the PR to post a detailed comparison of 
this functionality with scipy.ndimage.generic_filter, which appears to have 
very similar functionality.

Zach

 

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] [numpy] ENH: Initial implementation of a 'neighbor' calculation (#303)

[Zachary Pincus]

Hi Tim,

It's tricky to find functionality sometimes because as you've seen numpy and 
especially scipy are spread over very diverse domains with very diverse 
terminology! Usually someone on one or the other of the lists can help folks 
find some functionality, if not by name than by description...

In any case, though, I hope you'll keep your code around and accessible to 
people in standalone form. Despite being a bit slower than the ndimage stuff, 
it seems like you've got a better set of boundary conditions, and some other 
niceties that may appeal to others too.

Zach



> On Sun, Oct 14, 2012 at 8:24 PM, Zachary Pincus  
> wrote:
>> It would be useful for the author of the PR to post a detailed comparison of 
>> this functionality with scipy.ndimage.generic_filter, which appears to have 
>> very similar functionality.
>> 
> I'll be durned.   I created neighbor because I didn't find what I wanted, and 
> to find now that I just didn't look in the right place is well ...  Let's 
> just say that I went for a long run last night.
> 
> Searching for ndimage, I found that is has been around a long, long time.  
> First in numarray, then moved to scipy.
> 
> Really I could only nitpick about minor differences - kinda like a primary 
> political campaign.  On the face of it though, generic_filter looks better.  
> First off it is written in C so likely will be faster and more efficient 
> memory use.  I didn't look at optimizing neighbor at all and at least my part 
> of it is pure Python.  Of course for all of the small differences, I like my 
> choices better.  :-)
> 
> I would like to make a mild suggestion.  Emphasis on mild.  Maybe ndimage, in 
> all or in part, should be brought into (back into?) Numpy and renamed.  
> 
> About the PR.  Given that the neighbor functionality exists already, I will 
> close the PR later today.  Move along, nothing to see here...
> 
> Side note:  I wrote arraypad with the future idea that it would become 
> easyish to include that functionality in other places, for example in 
> neighbor.  A Don't Repeat Yourself idea.  Previously I had only seen Fortran 
> pad capabilities in some of the Fast Fourier Transform functions. The 
> generic_filter function includes several padding functions - written in C.  
> This means that if arraypad needs be more efficient we have C code to base a 
> better arraypad.
> 
> Another side node:  The measurements functions in ndimage are called zonal 
> functions in the GIS field.
> 
> Kindest regards,
> Tim
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Manipulate neighboring points in 2D array

[Zachary Pincus]

> I have 2D array, let's say: `np.random.random((100,100))` and I want to do 
> simple manipulation on each point neighbors, like divide their values by 3.
> 
> So for each array value, x, and it neighbors n:
> 
> n n nn/3 n/3 n/3
> n x n -> n/3  x  n/3
> n n nn/3 n/3 n/3
> 
> I searched a bit, and found about scipy ndimage filters, but if I'm not 
> wrong, there is no such function. Of course me being wrong is quite possible, 
> as I did not comprehend whole ndimage module, but I tried generic filter for 
> example and browser other functions.
> 
> Is there better way to make above manipulation, instead using for loop over 
> every array element?

I am not sure I understand the above manipulation... typically neighborhood 
operators take an array element and the its neighborhood and then give a single 
output that becomes the value of the new array at that point. That is, a 3x3 
neighborhood filter would act as a function F(R^{3x3}) -> R. It appears that 
what you're talking about above is a function F(R^{3x3}) -> R^{3x3}. But how is 
this output to map onto the original array positions? Is the function to be 
applied to non-overlapping neighborhoods? Is it to be applied to all 
neighborhoods and then summed at each position to give the output array?

If you can describe the problem in a bit more detail, with perhaps some sample 
input and output for what you desire (and/or with some pseudocode describing 
how it would work in a looping-over-each-element approach), I'm sure folks can 
figure out how best to do this in numpy.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Manipulate neighboring points in 2D array

[Zachary Pincus]

> You are right. I needed generic filter - to update current point, and not the 
> neighbors as I wrote.
> Initial code is slow loop over 2D python lists, which I'm trying to convert 
> to numpy and make it useful. In that loop there is inner loop for calculating 
> neighbors properties, which confused me yesterday, and mislead to search for 
> something that probably does not make sense.
> 
> It's clear now :)

It's possible that some generic filter operations can be cast in terms of 
pure-numpy operations, or composed out of existing filters available in 
scipy.ndimage. If you can describe the filter operation you wish to perform, 
perhaps someone can make some suggestions.

Alternately, scipy.ndimage.generic_filter can take an arbitrary python 
function. Though it's not really fast...

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] create a numpy array of images

[Zachary Pincus]

>>> I am using python for a while now and I have a requirement of  
>>> creating a
>>> numpy array of microscopic tiff images ( this data is 3d, meaning  
>>> there are
>>> 100 z slices of 512 X 512 pixels.) How can I create an array of  
>>> images?
>>
>> It's quite straightforward to create a 3-d array to hold this kind  
>> of data:
>>
>> image_block = np.empty((100, 512, 512), dtype=??)
>>
>> now you can load it up by using some lib (PIL, or ???) to load the  
>> tif
>> images, and then:
>>
>> for i in images:
>> image_block[i,:,:] = i
>
> Notice that since PIL 1.1.6, PIL Image objects support the numpy
> interface: http://effbot.org/zone/pil-changes-116.htm

For even longer than this, PIL has been somewhat broken with regard to  
16-bit images (very common in microscopy); you may run into strange  
byte-ordering issues that scramble the data on reading or writing.  
Also, PIL's numpy interface is somewhat broken in similar ways.  
(Numerous people have provided patches to PIL, but these are never  
incorporated into any releases, as far as I can tell.)

So try PIL, but if the images come out all wrong, you might want to  
check out the scikits.image package, which has hooks for various other  
image read/write tools.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] how to do this efficiently?

[Zachary Pincus]

> In a 1-d array, find the first point where all subsequent points  
> have values
> less than a threshold, T.

Maybe something like:

last_greater = numpy.arange(arr.shape)[arr >= T][-1]
first_lower = last_greater + 1

There's probably a better way to do it, without the arange, though...



___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] how to do this efficiently?

[Zachary Pincus]

>>> In a 1-d array, find the first point where all subsequent points
>>> have values
>>> less than a threshold, T.
>>
>> Maybe something like:
>>
>> last_greater = numpy.arange(arr.shape)[arr >= T][-1]
>> first_lower = last_greater + 1
>>
>> There's probably a better way to do it, without the arange, though...
>
> I'm trying to find the first point in a power spectrum such that all  
> subsequent
> points are below some level.  I've started with:
>
> db is my power spectrum in dB,   It is already reversed.
>
> mag = np.maximum.accumulate (db) - db[-1]
>
> Now all I need is to find the first point such that mag < -50.  How  
> to do this
> efficiently?

Right -- that's what I showed above. Find the last point in mag that  
is >= -50, and by definition the next point is the first point such  
that the remainder of mag is < -50.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] how to do this efficiently?

[Zachary Pincus]

On Feb 9, 2011, at 10:58 AM, Neal Becker wrote:

> Zachary Pincus wrote:
>
>>>>> In a 1-d array, find the first point where all subsequent points
>>>>> have values
>>>>> less than a threshold, T.
>>>>
>>>> Maybe something like:
>>>>
>>>> last_greater = numpy.arange(arr.shape)[arr >= T][-1]
>>>> first_lower = last_greater + 1
>>>>
>>>> There's probably a better way to do it, without the arange,  
>>>> though...
>>>
>>> I'm trying to find the first point in a power spectrum such that all
>>> subsequent
>>> points are below some level.  I've started with:
>>>
>>> db is my power spectrum in dB,   It is already reversed.
>>>
>>> mag = np.maximum.accumulate (db) - db[-1]
>>>
>>> Now all I need is to find the first point such that mag < -50.  How
>>> to do this
>>> efficiently?
>>
>> Right -- that's what I showed above. Find the last point in mag that
>> is >= -50, and by definition the next point is the first point such
>> that the remainder of mag is < -50.
>
> But where is numpy's 'find_first' function?  I can't seem to find it  
> so I had to
> make my own in C++.


As before, the line below does what you said you need, though not  
maximally efficiently. (Try it in an interpreter...) There may be  
another way in numpy that doesn't rely on constructing the index  
array, but this is the first thing that came to mind.

last_greater = numpy.arange(arr.shape)[arr >= T][-1]

Let's unpack that dense line a bit:

mask = arr >= T
indices = numpy.arange(arr.shape)
above_threshold_indices = indices[mask]
last_above_threshold_index = above_threshold_indices[-1]

Does this make sense?











___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] how to do this efficiently?

[Zachary Pincus]

>> As before, the line below does what you said you need, though not
>> maximally efficiently. (Try it in an interpreter...) There may be
>> another way in numpy that doesn't rely on constructing the index
>> array, but this is the first thing that came to mind.
>>
>> last_greater = numpy.arange(arr.shape)[arr >= T][-1]
>>
>> Let's unpack that dense line a bit:
>>
>> mask = arr >= T
>> indices = numpy.arange(arr.shape)
>> above_threshold_indices = indices[mask]
>> last_above_threshold_index = above_threshold_indices[-1]
>>
>> Does this make sense?
>
> This assumes monotonicity. Is that allowed?

The twice-stated problem was:

> In a 1-d array, find the first point where all subsequent points  
> have values less than a threshold, T.

So that should do the trick... Though Alan's argmax solution is  
definitely a better one than indexing an indices array. Same logic and  
result though, just more compact.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] how to do this efficiently?

[Zachary Pincus]

>> This assumes monotonicity. Is that allowed?
>
> The twice-stated problem was:

[Note to avert email-miscommunications] BTW, I wasn't trying to snipe  
at you with that comment, Josef...

I just meant to say that this solution solves the problem as Neal  
posed it, though that might not be the exact problem he has.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] how to do this efficiently?

[Zachary Pincus]

>
>>> In a 1-d array, find the first point where all subsequent points  
>>> have values
>>> less than a threshold.
>
> This doesn't imply monotonicity.
> Suppose with have a sin curve, and I want to find the last trough. Or
> a business cycle and I want to find the last recession.
>
> Unless my english deteriorated recently.
>

Not sure that I follow? I read the statement as:
find the smallest non-negative index i such that array[j] < t for all  
i <= j < len(array)

for which the various solutions proposed work, except for the corner  
case where array.min() >= t.

I'm not sure where monotonicity comes into play, unless one interprets  
the "all subsequent points" clause as "some number of subsequent  
points" or something...

For those sort of problems, I usually wind up smoothing the array  
[optional], and then using scipy.ndimage.maximum to calculate peak  
values within a given window, and then find the points that equal the  
peak values -- this works pretty robustly for peak/trough detection in  
arbitrary dimension.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] [OT] any image io module that works with python3?

[Zachary Pincus]

Here's a ctypes interface to FreeImage that I wrote a while back and  
was since cleaned up (and maintained) by the scikits.image folk:

https://github.com/stefanv/scikits.image/blob/master/scikits/image/io/_plugins/freeimage_plugin.py

If it doesn't work out of the box on python 3, then it should be  
pretty simple to fix.

Zach



On Mar 12, 2011, at 4:40 AM, Christoph Gohlke wrote:

>
>
> On 3/12/2011 1:08 AM, Nadav Horesh wrote:
>> Having numpy, scipy, and matplotlib working reasonably with  
>> python3, a
>> major piece of code I miss for a major python3 migration is an  
>> image IO.
>> I found that pylab's imread works fine for png image, but I need to  
>> read
>> all the other image format as well as png and jpeg output.
>> Any hints (including advices how easyly construct my own module) are
>> appreciated.
>> Nadav.
>>
>
> On Windows, PIL (private port at
> ), PythonMagick
> , and pygame 1.9.2pre
>  are working reasonably well for image IO. Also
> the FreeImage library  is easy to  
> use
> with ctypes .
>
> Christoph
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] strange behavior of np.minimum and np.maximum

[Zachary Pincus]

>
 a, b, c = np.array([10]), np.array([2]), np.array([7])
 min_val = np.minimum(a, b, c)
 min_val
> array([2])
 max_val = np.maximum(a, b, c)
 max_val
> array([10])
 min_val
> array([10])
>
> (I'm using numpy 1.4, and I observed the same behavior with numpy
> 2.0.0.dev8600 on another machine). I'm quite surprised by this  
> behavior
> (It took me quite a long time to figure out what was happening in a
> script of mine that wasn't giving what I expected, because of
> np.maximum changing the output of np.minimum). Is it a bug, or am I
> missing something?

Read the documentation for numpy.minimum and numpy.maximum: they give  
you element-wise minimum values from two arrays passed as arguments.  
E.g.:

 >>> numpy.minimum([1,2,3],[3,2,1])
array([1, 2, 1])

The optional third parameter to numpy.minimum is an "out" array - an  
array to place the results into instead of making a new array for that  
purpose. (This can save time / memory in various cases.)

This should therefore be enough to explain the above behavior. (That  
is, min_val and max_val wind up being just other names for the array  
'c', which gets modified in-place by the numpy.minimum and  
numpy.maximum.)

If you want the minimum value of a sequence of arbitrary length, use  
the python min() function. If you have a numpy array already and you  
want the minimum (global, or along a particular axis), use the min()  
method of the array, or numpy.min(arr).

Zach

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] reading tiff images

[Zachary Pincus]

On Apr 26, 2011, at 2:31 PM, Daniel Lepage wrote:

> You need PIL no matter what; scipy.misc.imread, scipy.ndimage.imread,
> and scikits.image.io.imread all call PIL.

scikits.image.io also has a ctypes wrapper for the freeimage library.  
I prefer these (well, I wrote them), though apparently there are some  
64-bit issues (crashes?). I haven't been working on a 64-bit system so  
I haven't been able to address them, but I will be soon. It's a very  
thin wrapper around a simple image IO library, so there's lots of room  
to add and extend as need be...

All of the PIL wrappers are kluges around serious flaws in how PIL  
reads images, particularly non-8-bit images and in particular non- 
native-endian 16-bit images.

Zach


> Theoretically there's no difference between any of them, although in
> actuality some use "import Image" and others use "from PIL import
> Image"; one of these may fail depending on how you installed PIL. (I'm
> not sure which is supposed to be standard - the PIL docs use both
> interchangeably, and I think the latest version of PIL on pypi sets it
> up so that both will work).
>
> I'd use whichever tool you're already importing - if you're using
> ndimage anyway, just use ndimage.imread rather than adding more
> imports.
>
> Note that using PIL directly is easy, but does require adding an extra
> step; OTOH, if you're familiar with PIL, you can use some of its
> transformations from the start, e.g.
>
> def imread(fname, mode='RGBA'):
>return np.asarray(Image.open(fname).convert(mode))
>
> to ensure that you always get 4-channel images, even for images that
> were initially RGB or grayscale.
>
> HTH,
> Dan
>
> On Tue, Apr 26, 2011 at 2:00 PM, Mathew Yeates  
>  wrote:
>> Hi
>> What is current method of using ndiimage on a Tiff file? I've seen
>> different methods using ndimage itself, scipy.misc and Pil.
>>
>> Mathew
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] fast grayscale conversion

[Zachary Pincus]

You could try:
src_mono = src_rgb.astype(float).sum(axis=-1) / 3.

But that speed does seem slow. Here are the relevant timings on my machine (a 
recent MacBook Pro) for a 3.1-megapixel-size array:
In [16]: a = numpy.empty((2048, 1536, 3), dtype=numpy.uint8)

In [17]: timeit numpy.dot(a.astype(float), numpy.ones(3)/3.)
10 loops, best of 3: 116 ms per loop

In [18]: timeit a.astype(float).sum(axis=-1)/3.
10 loops, best of 3: 85.3 ms per loop

In [19]: timeit a.astype(float)
10 loops, best of 3: 23.3 ms per loop




On Jun 20, 2011, at 4:15 PM, Alex Flint wrote:

> At the moment I'm using numpy.dot to convert a WxHx3 RGB image to a grayscale 
> image:
> 
> src_mono = np.dot(src_rgb.astype(np.float), np.ones(3)/3.);
> 
> This seems quite slow though (several seconds for a 3 megapixel image) - is 
> there a more specialized routine better suited to this?
> 
> Cheers,
> Alex
> 
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] poor performance of sum with sub-machine-word integer types

[Zachary Pincus]

Hello all,

As a result of the "fast greyscale conversion" thread, I noticed an anomaly 
with numpy.ndararray.sum(): summing along certain axes is much slower with 
sum() than versus doing it explicitly, but only with integer dtypes and when 
the size of the dtype is less than the machine word. I checked in 32-bit and 
64-bit modes and in both cases only once the dtype got as large as that did the 
speed difference go away. See below...

Is this something to do with numpy or something inexorable about machine / 
memory architecture?

Zach

Timings -- 64-bit mode:
--
In [2]: i = numpy.ones((1024,1024,4), numpy.int8)
In [3]: timeit i.sum(axis=-1)
10 loops, best of 3: 131 ms per loop
In [4]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 2.57 ms per loop

In [5]: i = numpy.ones((1024,1024,4), numpy.int16)
In [6]: timeit i.sum(axis=-1)
10 loops, best of 3: 131 ms per loop
In [7]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 4.75 ms per loop

In [8]: i = numpy.ones((1024,1024,4), numpy.int32)
In [9]: timeit i.sum(axis=-1)
10 loops, best of 3: 131 ms per loop
In [10]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 6.37 ms per loop

In [11]: i = numpy.ones((1024,1024,4), numpy.int64)
In [12]: timeit i.sum(axis=-1)
100 loops, best of 3: 16.6 ms per loop
In [13]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 15.1 ms per loop



Timings -- 32-bit mode:
--
In [2]: i = numpy.ones((1024,1024,4), numpy.int8)
In [3]: timeit i.sum(axis=-1)
10 loops, best of 3: 138 ms per loop
In [4]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 3.68 ms per loop

In [5]: i = numpy.ones((1024,1024,4), numpy.int16)
In [6]: timeit i.sum(axis=-1)
10 loops, best of 3: 140 ms per loop
In [7]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 4.17 ms per loop

In [8]: i = numpy.ones((1024,1024,4), numpy.int32)
In [9]: timeit i.sum(axis=-1)
10 loops, best of 3: 22.4 ms per loop
In [10]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
100 loops, best of 3: 12.2 ms per loop

In [11]: i = numpy.ones((1024,1024,4), numpy.int64)
In [12]: timeit i.sum(axis=-1)
10 loops, best of 3: 29.2 ms per loop
In [13]: timeit i[...,0]+i[...,1]+i[...,2]+i[...,3]
10 loops, best of 3: 23.8 ms per loop

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] poor performance of sum with sub-machine-word integer types

[Zachary Pincus]

On Jun 21, 2011, at 1:16 PM, Charles R Harris wrote:

> It's because of the type conversion sum uses by default for greater precision.

Aah, makes sense. Thanks for the detailed explanations and timings!
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Comparing NumPy/IDL Performance

[Zachary Pincus]

Hello Keith,

While I also echo Johann's points about the arbitrariness and non-utility of 
benchmarking I'll briefly comment just on just a few tests to help out with 
getting things into idiomatic python/numpy:

Tests 1 and 2 are fairly pointless (empty for loop and empty procedure) that 
won't actually influence the running time of well-written non-pathological code.

Test 3: 
#Test 3 - Add 20 scalar ints
nrep = 200 * scale_factor
for i in range(nrep):
a = i + 1

well, python looping is slow... one doesn't do such loops in idiomatic code if 
the underlying intent can be re-cast into array operations in numpy. But here 
the test is on such a simple operation that it's not clear how to recast in a 
way that would remain reasonable. Ideally you'd test something like:
i = numpy.arange(20)
for j in range(scale_factor):
  a = i + 1

but that sort of changes what the test is testing.


Finally, test 21:
#Test 21 - Smooth 512 by 512 byte array, 5x5 boxcar
for i in range(nrep):
b = scipy.ndimage.filters.median_filter(a, size=(5, 5))
timer.log('Smooth 512 by 512 byte array, 5x5 boxcar, %d times' % nrep)

A median filter is definitely NOT a boxcar filter! You want "uniform_filter":

In [4]: a = numpy.empty((1000,1000))

In [5]: timeit scipy.ndimage.filters.median_filter(a, size=(5, 5))
10 loops, best of 3: 93.2 ms per loop

In [6]: timeit scipy.ndimage.filters.uniform_filter(a, size=(5, 5))
10 loops, best of 3: 27.7 ms per loop

Zach


On Sep 26, 2011, at 10:19 AM, Keith Hughitt wrote:

> Hi all,
> 
> Myself and several colleagues have recently started work on a Python library 
> for solar physics, in order to provide an alternative to the current mainstay 
> for solar physics, which is written in IDL.
> 
> One of the first steps we have taken is to create a Python port of a popular 
> benchmark for IDL (time_test3) which measures performance for a variety of 
> (primarily matrix) operations. In our initial attempt, however, Python 
> performs significantly poorer than IDL for several of the tests. I have 
> attached a graph which shows the results for one machine: the x-axis is the 
> test # being compared, and the y-axis is the time it took to complete the 
> test, in milliseconds. While it is possible that this is simply due to 
> limitations in Python/Numpy, I suspect that this is due at least in part to 
> our lack in familiarity with NumPy and SciPy.
> 
> So my question is, does anyone see any places where we are doing things very 
> inefficiently in Python?
> 
> In order to try and ensure a fair comparison between IDL and Python there are 
> some things (e.g. the style of timing and output) which we have deliberately 
> chosen to do a certain way. In other cases, however, it is likely that we 
> just didn't know a better method.
> 
> Any feedback or suggestions people have would be greatly appreciated. 
> Unfortunately, due to the proprietary nature of IDL, we cannot share the 
> original version of time_test3, but hopefully the comments in time_test3.py 
> will be clear enough.
> 
> Thanks!
> Keith
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Comparing NumPy/IDL Performance

[Zachary Pincus]

I think the remaining delta between the integer and float "boxcar" smoothing is 
that the integer version (test 21) still uses median_filter(), while the float 
one (test 22) is using uniform_filter(), which is a boxcar.

Other than that and the slow roll() implementation in numpy, things look pretty 
solid, yes?

Zach


On Sep 29, 2011, at 12:11 PM, Keith Hughitt wrote:

> Thank you all for the comments and suggestions.
> 
> First off, I would like to say that I entirely agree with people's 
> suggestions about lack of objectiveness in the test design, and the caveat 
> about optimizing early. The main reason we put together the Python version of 
> the benchmark was as a quick "sanity check" to make sure that there are no 
> major show-stoppers before we began work on the library. We also wanted to 
> put together something to show other people who are firmly in the IDL camp 
> that this is a viable option.
> 
> We did in fact put together another short test-suite (test_testr.py & 
> time_testr.pro) which consists of operations that would are frequently used 
> by us, but it also is testing a very small portion of the kinds of things our 
> library will eventually do.
> 
> That said, I made a few small changes to the original benchmark, based on 
> people's feedback, and put together a new plot.
> 
> The changes made include:
> 
> 1. Using xrange instead of range
> 2. Using uniform filter instead of median filter
> 3. Fixed a typo for tests 2 & 3 which resulted in slower Python results
> 
> Again, note that some of the tests are testing non-numpy functionality. 
> Several of the results still stand out,  but overall the results are much 
> more reasonable than before.
> 
> Cheers,
> Keith
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] wanted: decent matplotlib alternative

[Zachary Pincus]

I keep meaning to use matplotlib as well, but every time I try I also get 
really turned off by the matlabish interface in the examples. I get that it's a 
selling point for matlab refugees, but I find it counterintuitive in the same 
way Christoph seems to.

I'm glad to hear the OO interface isn't as clunky as it looks on some of the 
doc pages, though. This is good news. Can anyone point out any good 
tutorials/docs on using matplotlib idiomatically via its OO interface?

Zach



On Oct 13, 2011, at 3:21 PM, Joe Kington wrote:

> Have a look at Chaco: http://code.enthought.com/chaco/  If you're wanting a 
> more pythonic api, it's a good choice.
> 
> Personally, I still prefer matplotlib.
> 
> You don't every need to touch the state machine interface. 
> 
> The OO interface is slighly un-pythonic, but it's hardly clunky. I think 
> you're referring to one of the webpage examples of it which avoids _any_ 
> convenience functions.  You can still use the convenience functions without 
> having to rely on the state machine in any way. E.g.:
> 
> import matplotlib.pyplot as plt
> 
> fig, axes = plt.subplots(ncols=4)
> 
> for ax in axes:
> ax.plot(range(10))
> 
> plt.show()
> 
> All in all, matplotlib deliberately tries to mimic matlab for a lot of the 
> conventions.  This is mostly to make it easier to switch if you're already 
> familiar with matlab.
> 
> To each his own, but for better or worse, matplotlib is the most widely used 
> plotting library for python.  It's worth getting a bit more familiar with, if 
> nothing else just to see past some of the rough edges.
> 
> -Joe
> 
> 
> 
> On Thu, Oct 13, 2011 at 1:58 PM, Christoph Groth  wrote:
> Hello,
> 
> Is it just me who thinks that matplotlib is ugly and a pain to use?  So
> far I haven't found a decent alternative usable from within python.  (I
> haven't tried all the packages out there.)  I'm mostly interested in 2d
> plots.  Who is happy enough with a numpy-compatible plotting package to
> recommend it?
> 
> Thanks,
> 
> Christoph
> 
> 
> 
> 
> A few things I can't stand about matplotlib:
> 
>  * It works as a state machine.  (There is an OO-API, too, but it's ugly
>   and cumbersome to use, and most examples use the state machine mode
>   or, even worse, a mixture of OO and global state.)
> 
>  * It uses inches by default. (I propose to switch to nails: 1 nail = 3
>   digits = 2 1/4 inches = 1/16 yard.)
> 
>  * subplot(211)  (ugh!)
> 
>  * Concepts are named in a confusing way. ("ax = subplot(112)" anyone?)
> 
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> 
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] NumPy nogil API

[Zachary Pincus]

> As an example, it'd be nice to have scipy.ndimage available without the GIL:
> http://docs.scipy.org/doc/scipy/reference/ndimage.html
> 
> Now, this *can* easily be done as the core is written in C++. I'm just 
> pointing out that some people may wish more for calling scipy.ndimage 
> inside their prange than for some parts of NumPy.

Not exactly to your larger point wrt the GIL, but I *think* some skimage (née 
scikits.image) folks are trying to rewrite most of ndimage's functionality in 
cython. I don't know what the status of this effort is though...

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] dedicated function for resize with averaging or rebin 2d arrays?

[Zachary Pincus]

Hi Andrea,

scipy.ndimage.zoom will do this nicely for magnification. (Just set the spline 
order to 0 to get nearest-neighbor interpolation; otherwise you can use higher 
orders for better smoothing.)

For decimation (zooming out) scipy.ndimage.zoom also works, but it's not as 
nice as a dedicated decimation filter that would average properly over the area 
that's being squeezed into a single output pixel. (You'd have to choose the 
spline order manually to approximate that.) I'm afraid I don't have enough 
signal-processing background to know how to write a proper general-purpose 
decimation filter -- basically, you convolve with whatever bandlimiting filter 
(e.g. a gaussian, or do it in the Fourier domain), then just do 
nearest-neighbor downsampling, but I'm never sure how to properly choose the 
filter parameters!

Between this and ndimage.zoom for magnifying, one could get together a much 
better "rebin" function that in the edge cases of integer 
magnification/minification should work the same as the IDL one. But the 
participants in the old discussion you highlighted seemed unhappy with the 
time/space used for proper decimation, so I'm not sure what really would be 
best. 

Zach


On Nov 11, 2011, at 1:41 AM, Andrea Zonca wrote:

> hi,
> I work in astrophysics where the most common programming language is
> currently IDL.
> A common request of people switching from IDL to python is the
> implementation of the REBIN function, which either downsizes a 2d
> array by averaging or increases its dimension by repeating its
> elements. In both cases the new shape must be an integer factor of the
> old shape.
> 
> I believe it is a very handy function for quick smoothing of 2 dimensional 
> data.
> 
> I found a discussion about this topic in the archives:
> http://thread.gmane.org/gmane.comp.python.numeric.general/885/focus=894
> 
> Do you think it would be useful to add such function to numpy?
> 
> I created a simple implementation to help in the discussion:
> https://gist.github.com/1348792
> 
> thanks,
> Andrea Zonca
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] bug in numpy.mean() ?

[Zachary Pincus]

On Jan 24, 2012, at 1:33 PM, K.-Michael Aye wrote:

> I know I know, that's pretty outrageous to even suggest, but please 
> bear with me, I am stumped as you may be:
> 
> 2-D data file here:
> http://dl.dropbox.com/u/139035/data.npy
> 
> Then:
> In [3]: data.mean()
> Out[3]: 3067.024383998
> 
> In [4]: data.max()
> Out[4]: 3052.4343
> 
> In [5]: data.shape
> Out[5]: (1000, 1000)
> 
> In [6]: data.min()
> Out[6]: 3040.498
> 
> In [7]: data.dtype
> Out[7]: dtype('float32')
> 
> 
> A mean value calculated per loop over the data gives me 3045.747251076416
> I first thought I still misunderstand how data.mean() works, per axis 
> and so on, but did the same with a flattenend version with the same 
> results.
> 
> Am I really soo tired that I can't see what I am doing wrong here?
> For completion, the data was read by a osgeo.gdal dataset method called 
> ReadAsArray()
> My numpy.__version__ gives me 1.6.1 and my whole setup is based on 
> Enthought's EPD.


I get the same result:

In [1]: import numpy

In [2]: data = numpy.load('data.npy')

In [3]: data.mean()
Out[3]: 3067.024383998

In [4]: data.max()
Out[4]: 3052.4343

In [5]: data.min()
Out[5]: 3040.498

In [6]: numpy.version.version
Out[6]: '2.0.0.dev-433b02a'

This on OS X 10.7.2 with Python 2.7.1, on an intel Core i7. Running python as a 
32 vs. 64-bit process doesn't make a difference.

The data matrix doesn't look too strange when I view it as an image -- all 
pretty smooth variation around the (min, max) range. But maybe it's still 
somehow floating-point pathological?

This is fun too:
In [12]: data.mean()
Out[12]: 3067.024383998

In [13]: (data/3000).mean()*3000
Out[13]: 3020.807437501

In [15]: (data/2).mean()*2
Out[15]: 3067.024383998

In [16]: (data/200).mean()*200
Out[16]: 3013.67541


Zach


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] bug in numpy.mean() ?

[Zachary Pincus]

> You have a million 32-bit floating point numbers that are in the 
> thousands. Thus you are exceeding the 32-bitfloat precision and, if you 
> can, you need to increase precision of the accumulator in np.mean() or 
> change the input dtype:
 a.mean(dtype=np.float32) # default and lacks precision
> 3067.024383998
 a.mean(dtype=np.float64)
> 3045.747251076416
 a.mean(dtype=np.float128)
> 3045.7472510764160156
 b=a.astype(np.float128)
 b.mean()
> 3045.7472510764160156
> 
> Otherwise you are left to using some alternative approach to calculate 
> the mean.
> 
> Bruce

Interesting -- I knew that float64 accumulators were used with integer arrays, 
and I had just assumed that 64-bit or higher accumulators would be used with 
floating-point arrays too, instead of the array's dtype. This is actually quite 
a bit of a gotcha for floating-point imaging-type tasks -- good to know!

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Addressing arrays

[Zachary Pincus]

a[x,y,:]

Read the slicing part of the tutorial:
http://www.scipy.org/Tentative_NumPy_Tutorial 
(section 1.6)

And the documentation:
http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html



On Jan 30, 2012, at 10:25 AM, Ted To wrote:

> Hi,
> 
> Is there some straightforward way to access an array by values across a
> subset of its dimensions?  For example, if I have a three dimensional
> array a=(x,y,z), can I look at the values of z given particular values
> for x and y?
> 
> Thanks,
> Ted
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Addressing arrays

[Zachary Pincus]

Ted, can you clarify what you're asking for? Maybe give a trivial example of an 
array and the desired output?

I'm pretty sure this is a slicing question though:
> If I have a three dimensional array a=(x,y,z), can I look at the values of z 
> given particular values for x and y?
Given that element values are scalars in this case, and indices are (x,y,z) 
triples, it seems likely that looking for "values of z" given an (x,y) pair is 
an slicing-by-index question, no?

For indexing-by-value, "fancy indexing" with boolean masks is usually the way 
to go... again, Ted (or Chao), if you can describe your indexing needs in a bit 
more detail, it's often easy to find a compact slicing and/or fancy-indexing 
strategy that works well and reasonably efficiently.

Zach



On Jan 30, 2012, at 10:33 AM, Chao YUE wrote:

> he is not asking for slicing. he is asking for how to index array by element 
> value but not element index.
> 
> 2012/1/30 Zachary Pincus 
> a[x,y,:]
> 
> Read the slicing part of the tutorial:
> http://www.scipy.org/Tentative_NumPy_Tutorial
> (section 1.6)
> 
> And the documentation:
> http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html
> 
> 
> 
> On Jan 30, 2012, at 10:25 AM, Ted To wrote:
> 
> > Hi,
> >
> > Is there some straightforward way to access an array by values across a
> > subset of its dimensions?  For example, if I have a three dimensional
> > array a=(x,y,z), can I look at the values of z given particular values
> > for x and y?
> >
> > Thanks,
> > Ted
> > ___
> > NumPy-Discussion mailing list
> > NumPy-Discussion@scipy.org
> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> 
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> 
> 
> 
> -- 
> ***
> Chao YUE
> Laboratoire des Sciences du Climat et de l'Environnement (LSCE-IPSL)
> UMR 1572 CEA-CNRS-UVSQ
> Batiment 712 - Pe 119
> 91191 GIF Sur YVETTE Cedex
> Tel: (33) 01 69 08 29 02; Fax:01.69.08.77.16
> 
> 
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Addressing arrays

[Zachary Pincus]

> Thanks!  That works great if I only want to search over one index but I
> can't quite figure out what to do with more than a single index.  So
> suppose I have a labeled, multidimensional array with labels 'month',
> 'year' and 'quantity'.  a[['month','year']] gives me an array of indices
> but "a[['month','year']]==(1,1960)" produces "False".  I'm sure I simply
> don't know the proper syntax and I apologize for that -- I'm kind of new
> to numpy.

I think that your best bet is to form the boolean masks independently and then 
logical-and them together:

mask = (a['month'] == 1) & (a['year'] == 1960)
jan_60 = a[mask]

Someone might have more insight here. Though I should note that if you have 
large data and are doing lots of "queries" like this, a more database-ish 
approach might be better. Something like sqlite's python bindings, or PyTables. 
Alternately, if your data are all time-series based things, PANDAS might be 
worth looking at. 

But the above approach should be just fine for non-huge datasets...

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] sampling arrays

[Zachary Pincus]

> Thank you, 
>  No if the location ( space time or depth) of choice is not 
> available then the function I was looking for should give an interpolated 
> value at the choice.
> with best regards,
> Sudheer

scipy.ndimage.map_coordinates may be exactly what you want.


> - Original Message -
> 
>> From: Henry Gomersall 
>> To: Discussion of Numerical Python 
>> Cc: 
>> Sent: Sunday, 16 June 2013 2:49 PM
>> Subject: Re: [Numpy-discussion] sampling arrays
>> 
>> On Sun, 2013-06-16 at 14:48 +0800, Sudheer Joseph wrote:
>>> Is it possible to sample a 4D array of numpy at given dimensions with
>>> out writing loops? ie a smart python way?
>> 
>> It's not clear how what you want to do is different from simply indexing
>> the array...?
>> 
>> Henry
>> 
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>> 
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Relative speed

[Zachary Pincus]

> And as you pointed out, 
> most of the time for non-trivial datasets the numpy operations will be 
> faster. (I'm daunted by the notion of trying to do linear algebra on 
> lists of tuples, assuming that's the relevant set of operations given 
> the comparison to the matrix class.)

Note the important and pretty common exception of building up a list one 
element (or row of elements) at a time. Here, python lists usually rule, unless 
the final size is known in advance.

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] u in [u+1]

[Zachary Pincus]

> I'm having some trouble here.  I have a list of numpy arrays.  I  
> want to
> know if an array 'u' is in the list.

Try:

any(numpy.all(u == l) for l in array_list)

standard caveats about float comparisons apply; perhaps
any(numpy.allclose(u, l) for l in array_list)
is more appropriate in certain circumstances.

Can of course replace the first 'any' with 'all' or 'sum' to get  
different kinds of information, but using 'any' is equivalent to the  
'in' query that you wanted.

Why the 'in' operator below fails is that behind the scenes, 'u not in  
[u+1]' causes Python to iterate through the list testing each element  
for equality with u. Except that as the error states, arrays don't  
support testing for equality because such tests are ambiguous. (cf.  
many threads about this.)

Zach


On Feb 5, 2010, at 6:47 AM, Neal Becker wrote:

> I'm having some trouble here.  I have a list of numpy arrays.  I  
> want to
> know if an array 'u' is in the list.
>
> As an example,
> u = np.arange(10)
>
> : u not in [u+1]
> ---
> ValueErrorTraceback (most recent  
> call last)
>
> /home/nbecker/raysat/test/ in ()
>
> ValueError: The truth value of an array with more than one element is
> ambiguous. Use a.any() or a.all()
>
> What would be the way to do this?
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Loading bit strings

[Zachary Pincus]

> Is there a good way in NumPy to convert from a bit string to a boolean
> array?
>
> For example, if I have a 2-byte string s='\xfd\x32', I want to get a
> 16-length boolean array out of it.

numpy.unpackbits(numpy.fromstring('\xfd\x32', dtype=numpy.uint8))
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] numpy.array(arr.flat) mutates arr if arr.flags.fortran: bug?

[Zachary Pincus]

Hello,

I assume it is a bug that calling numpy.array() on a flatiter of a  
fortran-strided array that owns its own data causes that array to be  
rearranged somehow?

Not sure what happens with a fancier-strided array that also owns its  
own data (because I'm not sure how to create one of those in python).

This is from the latest svn version (2.0.0.dev8302) but was also  
present in a previous version too.

Zach


In [9]: a = numpy.array([[1,2],[3,4]]).copy('F')

In [10]: a
Out[10]:
array([[1, 2],
[3, 4]])

In [11]: list(a.flat)
Out[11]: [1, 2, 3, 4]

In [12]: a # no problem
Out[12]:
array([[1, 2],
[3, 4]])

In [13]: numpy.array(a.flat)
Out[13]: array([1, 2, 3, 4])

In [14]: a # this ain't right!
Out[14]:
array([[1, 3],
[2, 4]])

In [15]: a = numpy.array([[1,2],[3,4]]).copy('C')

In [16]: numpy.array(a.flat)
Out[16]: array([1, 2, 3, 4])

In [17]: a
Out[17]:
array([[1, 2],
[3, 4]])

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] numpy.array(arr.flat) mutates arr if arr.flags.fortran: bug?

[Zachary Pincus]

> You should open a ticket for this.

http://projects.scipy.org/numpy/ticket/1439



On Mar 26, 2010, at 11:26 AM, Charles R Harris wrote:

>
>
> On Wed, Mar 24, 2010 at 1:13 PM, Zachary Pincus  > wrote:
> Hello,
>
> I assume it is a bug that calling numpy.array() on a flatiter of a
> fortran-strided array that owns its own data causes that array to be
> rearranged somehow?
>
> Not sure what happens with a fancier-strided array that also owns its
> own data (because I'm not sure how to create one of those in python).
>
> This is from the latest svn version (2.0.0.dev8302) but was also
> present in a previous version too.
>
> Zach
>
>
> In [9]: a = numpy.array([[1,2],[3,4]]).copy('F')
>
> In [10]: a
> Out[10]:
> array([[1, 2],
>[3, 4]])
>
> In [11]: list(a.flat)
> Out[11]: [1, 2, 3, 4]
>
> In [12]: a # no problem
> Out[12]:
> array([[1, 2],
>[3, 4]])
>
> In [13]: numpy.array(a.flat)
> Out[13]: array([1, 2, 3, 4])
>
> In [14]: a # this ain't right!
> Out[14]:
> array([[1, 3],
>[2, 4]])
>
> In [15]: a = numpy.array([[1,2],[3,4]]).copy('C')
>
> In [16]: numpy.array(a.flat)
> Out[16]: array([1, 2, 3, 4])
>
> In [17]: a
> Out[17]:
> array([[1, 2],
>[3, 4]])
>
>
> You should open a ticket for this.
>
> Chuck
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Replacing NANs

[Zachary Pincus]

> In an array I want to replace all NANs with some number say 100, I  
> found a method nan_to_num but it only replaces with zero.
> Any solution for this?

Indexing with a mask is one approach here:

a[numpy.isnan(a)] = 100

also cf. numpy.isfinite as well in case you want the same with infs.

Zach


On Mar 30, 2010, at 2:59 PM, Vishal Rana wrote:

> Hi,
>
> In an array I want to replace all NANs with some number say 100, I  
> found a method nan_to_num but it only replaces with zero.
> Any solution for this?
>
> Thanks
> Vishal
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] How to combine a pair of 1D arrays?

[Zachary Pincus]

> On Wed, Apr 14, 2010 at 10:25, Peter Shinners   
> wrote:
>> Is there a way to combine two 1D arrays with the same size into a 2D
>> array? It seems like the internal pointers and strides could be
>> combined. My primary goal is to not make any copies of the data.
>
> There is absolutely no way to get around that, I am afraid.

You could start with the 2d array... instead of:
>  >>> a = np.array((1,2,3,4))
>  >>> b = np.array((11,12,13,14))
>  >>> c = np.magical_fuse(a, b)   # what goes here?

perhaps something like:

  >>> c = np.empty((2,4), int)
  >>> a = c[0]
  >>> b = c[1]

Now a and b are views on c. So if you (say) pass them to, say, some  
(strides-aware) C routine that fills them in element-by-element, c  
will get filled in at the same time.

Zach


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] 2D binning

[Zachary Pincus]

> Hi
> Can anyone think of a clever (non-lopping) solution to the following?
>
> A have a list of latitudes, a list of longitudes, and list of data  
> values. All lists are the same length.
>
> I want to compute an average  of data values for each lat/lon pair.  
> e.g. if lat[1001] lon[1001] = lat[2001] [lon [2001] then
> data[1001] = (data[1001] + data[2001])/2
>
> Looping is going to take wa to long.

As a start, are the "equal" lat/lon pairs exactly equal (i.e. either  
not floating-point, or floats that will always compare equal, that is,  
the floating-point bit-patterns will be guaranteed to be identical) or  
approximately equal to float tolerance?

If you're in the approx-equal case, then look at the KD-tree in scipy  
for doing near-neighbors queries.

If you're in the exact-equal case, you could consider hashing the lat/ 
lon pairs or something. At least then the looping is O(N) and not  
O(N^2):

import collections
grouped = collections.defaultdict(list)
for lt, ln, da in zip(lat, lon, data):
   grouped[(lt, ln)].append(da)

averaged = dict((ltln, numpy.mean(da)) for ltln, da in grouped.items())

Is that fast enough?

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] 2D binning

[Zachary Pincus]

> I guess it's as fast as I'm going to get. I don't really see any  
> other way. BTW, the lat/lons are integers)

You could (in c or cython) try a brain-dead "hashtable" with no  
collision detection:

for lat, long, data in dataset:
   bin = (lat ^ long) % num_bins
   hashtable[bin] = update_incremental_mean(hashtable[bin], data)

you'll of course want to do some experiments to see if your data are  
sufficiently sparse and/or you can afford a large enough hashtable  
array that you won't get spurious hash collisions. Adding error- 
checking to ensure that there are no collisions would be pretty  
trivial (just keep a table of the lat/long for each hash value, which  
you'll need anyway, and check that different lat/long pairs don't get  
assigned the same bin).

Zach



> -Mathew
>
> On Tue, Jun 1, 2010 at 1:49 PM, Zachary Pincus  > wrote:
> > Hi
> > Can anyone think of a clever (non-lopping) solution to the  
> following?
> >
> > A have a list of latitudes, a list of longitudes, and list of data
> > values. All lists are the same length.
> >
> > I want to compute an average  of data values for each lat/lon pair.
> > e.g. if lat[1001] lon[1001] = lat[2001] [lon [2001] then
> > data[1001] = (data[1001] + data[2001])/2
> >
> > Looping is going to take wa to long.
>
> As a start, are the "equal" lat/lon pairs exactly equal (i.e. either
> not floating-point, or floats that will always compare equal, that is,
> the floating-point bit-patterns will be guaranteed to be identical) or
> approximately equal to float tolerance?
>
> If you're in the approx-equal case, then look at the KD-tree in scipy
> for doing near-neighbors queries.
>
> If you're in the exact-equal case, you could consider hashing the lat/
> lon pairs or something. At least then the looping is O(N) and not
> O(N^2):
>
> import collections
> grouped = collections.defaultdict(list)
> for lt, ln, da in zip(lat, lon, data):
>   grouped[(lt, ln)].append(da)
>
> averaged = dict((ltln, numpy.mean(da)) for ltln, da in  
> grouped.items())
>
> Is that fast enough?
>
> Zach
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Installing numpy on py 3.1.2 , osx

[Zachary Pincus]

Hi Vincent,

(1) Fortran compiler isn't necessary for numpy, but is for scipy,  
which isn't ported to python 3 yet.
(2) Could you put up on pastebin or somewhere online the full error  
you got?

The problem isn't one of not finding the Python.h header file, which  
will be present in the Python 3 framework that you installed, but some  
sort of obscure compilation error that got misinterpreted as a missing  
header by the python build process. Unfortunately, the [CLIP] in your  
previous email may have obscured the source of the error.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Installing numpy on py 3.1.2 , osx

[Zachary Pincus]

On Jun 7, 2010, at 5:19 PM, Vincent Davis wrote:

> Here is a link to the full output after typing python setup.py build.
> https://docs.google.com/Doc?docid=0AVQgwG2qUDgdZGYyaGo0NjNfMjI5Z3BraHd6ZDg&hl=en

that's just bringing up an empty document page for me...
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Installing numpy on py 3.1.2 , osx

[Zachary Pincus]

This is unexpected, from the error log:
> /Library/Frameworks/Python.framework/Versions/3.1/include/python3.1/ 
> Python.h:11:20: error: limits.h: No such file or directory

No good... it can't find basic system headers. Perhaps it's due to the  
MACOSX_DEPLOYMENT_TARGET environment variable that was warned about  
earlier? (I think I'd seen that sort of thing before... regardless of  
whether you have the 10.4 SDK, if MACOSX_DEPLOYMENT_TARGET is 10.3 I  
think it doesn't get used? Not sure.) Try setting  
MACOSX_DEPLOYMENT_TARGET to 10.4:
MACOSX_DEPLOYMENT_TARGET=10.4
export MACOSX_DEPLOYMENT_TARGET
(if you're using bash, the default OS X shell; for tcsh use "setenv  
MACOSX_DEPLOYMENT_TARGET 10.4")

If this fixes things, I suspect that the build process should probably  
give an error about MACOSX_DEPLOYMENT_TARGET instead of a warning...

Zach



On Jun 7, 2010, at 11:49 PM, Vincent Davis wrote:

> On Mon, Jun 7, 2010 at 7:59 PM, Vincent Davis  > wrote:
>> On Mon, Jun 7, 2010 at 7:26 PM, Ralf Gommers
>>  wrote:
>>>
>>>
>>> On Tue, Jun 8, 2010 at 4:46 AM, Vincent Davis >> >
>>> wrote:

 On Mon, Jun 7, 2010 at 2:11 PM, Skipper Seabold >>> >
 wrote:
> On Mon, Jun 7, 2010 at 4:02 PM, Vincent Davis  >
> wrote:
>>   "Cannot compile 'Python.h'. Perhaps you need to "\
>> SystemError: Cannot compile 'Python.h'. Perhaps you need to  
>> install
>> python-dev|python-devel.
>
> Not sure about this.  Might be the compiler issue, but I know on  
> Linux
> you need to install the python-dev (or python-devel) package,  
> but I
> don't see that as a pre-req on the Mac page.

 So it seems to be seeing the newly installed fortran, but I do  
 get an
 error at the end.
  File "numpy/core/setup.py", line 260, in check_types
"Cannot compile 'Python.h'. Perhaps you need to "\
 SystemError: Cannot compile 'Python.h'. Perhaps you need to install
 python-dev|python-devel.

 Any ideas, I looked into how python-dev|python-devel applied to  
 python
 but could make much of what I found. I tried to build python 3.1.2
 from source but that didn't work out. If I need to go that  
 direction I
 guess I need to.

>>> You are probably missing the 10.4 SDK. It's an optional install of  
>>> XCode on
>>> your Snow Leopard DVD.
>>
>> True I have 10.5 and 10.6 SDK I will install the 10.4 and try again.
>> Thanks
>> Vincent
>
> Still no go after installing 10.4 sdk
> Attached terminal output.
>
>>>
>>> Cheers,
>>> Ralf
>>>
>>>
>>>
>>> ___
>>> NumPy-Discussion mailing list
>>> NumPy-Discussion@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>>
>>>
>>
>  Output.txt.zip>___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Installing numpy on py 3.1.2 , osx

[Zachary Pincus]

> On Tue, Jun 8, 2010 at 7:58 AM, Zachary Pincus  > wrote:
>> This is unexpected, from the error log:
>>> /Library/Frameworks/Python.framework/Versions/3.1/include/python3.1/
>>> Python.h:11:20: error: limits.h: No such file or directory
>>
>> No good... it can't find basic system headers. Perhaps it's due to  
>> the
>> MACOSX_DEPLOYMENT_TARGET environment variable that was warned about
>> earlier? (I think I'd seen that sort of thing before... regardless of
>> whether you have the 10.4 SDK, if MACOSX_DEPLOYMENT_TARGET is 10.3 I
>> think it doesn't get used? Not sure.) Try setting
>> MACOSX_DEPLOYMENT_TARGET to 10.4:
>> MACOSX_DEPLOYMENT_TARGET=10.4
>> export MACOSX_DEPLOYMENT_TARGET
>> (if you're using bash, the default OS X shell; for tcsh use "setenv
>> MACOSX_DEPLOYMENT_TARGET 10.4")
>>
> Before I make any changes I wanted to make sure I kinda know what I am
> doing and how :)
> I have several versions of python installed. The python 3.1.2 is clean
> except for nose being installed. to access python 3 on my system I
> have a terminal alias to py 3.1.2 as python3

Basic summary of environment variables:
http://en.wikipedia.org/wiki/Environment_variable

In a unix shell, setting an environment variable will apply to the  
processes that are started from within that shell session. To  
permanently set environment variables, you'd need to add the relevant  
commands to the bash .profile or other startup files (a topic in and  
of itself). So the good news is that nothing you do will w/r/t  
environment vars will persist past exiting that particular shell.

Anyhow, from the below, you don't have MACOSX_DEPLOYMENT_TARGET set at  
all, which is what would be expected. I wonder if perhaps the python/ 
numpy build process is setting MACOSX_DEPLOYMENT_TARGET to 10.3 when  
there is none already set? If so this is not good.

Anyhow, given that you're using bash, here's what to do to try to build:

cd /path/to/numpy/src
MACOSX_DEPLOYMENT_TARGET=10.4
export MACOSX_DEPLOYMENT_TARGET
python3 setup.py build

and see if that fixes the issue? If so there's some bug in the numpy  
build process you've stumbled across.

Zach



>
> So when in the terminal typing env returns
> MacBookPro-new-2:numpy vmd$ env
> TERM_PROGRAM=Apple_Terminal
> TERM=xterm-color
> SHELL=/bin/bash
> TMPDIR=/var/folders/2f/2fiXYQSSE+CgAzDQPp9+kTI/-Tmp-/
> Apple_PubSub_Socket_Render=/tmp/launch-4g8LJN/Render
> TERM_PROGRAM_VERSION=273
> OLDPWD=/Users/vmd/DropBox/numpy
> USER=vmd
> COMMAND_MODE=unix2003
> SSH_AUTH_SOCK=/tmp/launch-fxsj0I/Listeners
> __CF_USER_TEXT_ENCODING=0x1F5:0:0
> PATH=/Library/Frameworks/EPD64.framework/Versions/Current/bin:/ 
> Library/Frameworks/EPD64.framework/Versions/Current/bin:/usr/bin:/ 
> bin:/usr/sbin:/sbin:/usr/local/bin:/usr/local/git/bin:/usr/X11/bin:/ 
> opt/local/bin
> MKL_NUM_THREADS=1
> PWD=/Users/vmd/DropBox/numpy/numpy
> LANG=en_US.UTF-8
> SHLVL=1
> HOME=/Users/vmd
> LOGNAME=vmd
> DISPLAY=/tmp/launch-HlM7T8/org.x:0
> _=/usr/bin/env
> MacBookPro-new-2:numpy vmd$
>
> So looking into this a bit more.
> MacBookPro-new-2:numpy vmd$ python3
> Python 3.1.2 (r312:79360M, Mar 24 2010, 01:33:18)
> [GCC 4.0.1 (Apple Inc. build 5493)] on darwin
> Type "help", "copyright", "credits" or "license" for more information.
>
>>>> import os
>>>> for x in os.environ.items(): print(x)
> ...
> ('MKL_NUM_THREADS', '1')
> ('LANG', 'en_US.UTF-8')
> ('TERM', 'xterm-color')
> ('Apple_PubSub_Socket_Render', '/tmp/launch-4g8LJN/Render')
> ('TMPDIR', '/var/folders/2f/2fiXYQSSE+CgAzDQPp9+kTI/-Tmp-/')
> ('SHLVL', '1')
> ('OLDPWD', '/Users/vmd/DropBox/numpy')
> ('SSH_AUTH_SOCK', '/tmp/launch-fxsj0I/Listeners')
> ('TERM_PROGRAM_VERSION', '273')
> ('__CF_USER_TEXT_ENCODING', '0x1F5:0:0')
> ('PWD', '/Users/vmd/DropBox/numpy/numpy')
> ('SHELL', '/bin/bash')
> ('LOGNAME', 'vmd')
> ('USER', 'vmd')
> ('PATH', '/Library/Frameworks/EPD64.framework/Versions/Current/bin:/ 
> Library/Frameworks/EPD64.framework/Versions/Current/bin:/usr/bin:/ 
> bin:/usr/sbin:/sbin:/usr/local/bin:/usr/local/git/bin:/usr/X11/bin:/ 
> opt/local/bin')
> ('HOME', '/Users/vmd')
> ('TERM_PROGRAM', 'Apple_Terminal')
> ('DISPLAY', '/tmp/launch-HlM7T8/org.x:0')
> ('_', '/usr/local/b

Re: [Numpy-discussion] Installing numpy on py 3.1.2 , osx

[Zachary Pincus]

> Failed again, I have attached the output including the execution of
> the above commands.
> Thanks for link to the environment variables, I need to read that.

In the attached file (and the one from the next email too) I didn't  
see the

MACOSX_DEPLOYMENT_TARGET=10.4
export MACOSX_DEPLOYMENT_TARGET

lines. Those need to be re-run with each new shell (terminal window)  
you start up. (Unless you put them in ~/.profile, etc.)

You can see from the error messages (after the barf from the 2to3 tool  
finishes) that there's still the warning about  
MACOSX_DEPLOYMENT_TARGET being set to 10.3.

If you did in fact run those commands immediately prior to "python3  
setup.py build" and still see the MACOSX_DEPLOYMENT_TARGET warning,  
then there is definitely some bug in the build process that's  
obliterating that environment variable.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Installing numpy on py 3.1.2 , osx

[Zachary Pincus]

Hi Vincent,

I'm not really sure -- now the build is using the 10.4 SDK but still  
not finding the right headers.

Could you check to see that the following file exists:
/Developer/SDKs/MacOSX10.4u.sdk/usr/include/limits.h

If it does, then I'm really not sure what's going on.

Zach



On Jun 8, 2010, at 3:35 PM, Vincent Davis wrote:

> On Tue, Jun 8, 2010 at 12:26 PM, Zachary Pincus  > wrote:
>>
>>> Failed again, I have attached the output including the execution of
>>> the above commands.
>>> Thanks for link to the environment variables, I need to read that.
>>
>> In the attached file (and the one from the next email too) I didn't
>> see the
>>
>> MACOSX_DEPLOYMENT_TARGET=10.4
>> export MACOSX_DEPLOYMENT_TARGET
>>
>> lines. Those need to be re-run with each new shell (terminal window)
>> you start up. (Unless you put them in ~/.profile, etc.)
>>
>> You can see from the error messages (after the barf from the 2to3  
>> tool
>> finishes) that there's still the warning about
>> MACOSX_DEPLOYMENT_TARGET being set to 10.3.
>>
>> If you did in fact run those commands immediately prior to "python3
>> setup.py build" and still see the MACOSX_DEPLOYMENT_TARGET warning,
>> then there is definitely some bug in the build process that's
>> obliterating that environment variable.
>
> I did but somehow what was being saved from the terminal was from  
> yesterday ???
> So to be sure I started over.
> I have run it again and for sure I have entered the commands, (see
> attached evidence)
> MACOSX_DEPLOYMENT_TARGET=10.4
> export MACOSX_DEPLOYMENT_TARGET
>
> I then search the terminal output for "MACOSX_DEPLOYMENT_TARGET" and
> it only occurred at the beginning where I entered the commands.
>
> Thanks
> Vincent
>
>>
>> Zach
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>  V3.txt.zip>___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] match or vectorized in-type function.

[Zachary Pincus]

> match(v1, v2) => returns a boolean array of length len(v1) indicating
> whether element i in v1 is in v2.

You want numpy.in1d (and friends, probably, like numpy.unique and the  
others that are all collected in numpy.lib.arraysetops...)


Definition:   numpy.in1d(ar1, ar2, assume_unique=False)
Docstring:
 Test whether each element of a 1D array is also present in a  
second array.

 Returns a boolean array the same length as `ar1` that is True
 where an element of `ar1` is in `ar2` and False otherwise.

 Parameters
 --
 ar1 : array_like, shape (M,)
 Input array.
 ar2 : array_like
 The values against which to test each value of `ar1`.
 assume_unique : bool, optional
 If True, the input arrays are both assumed to be unique, which
 can speed up the calculation.  Default is False.

 Returns
 ---
 mask : ndarray of bools, shape(M,)
 The values `ar1[mask]` are in `ar2`.

 See Also
 
 numpy.lib.arraysetops : Module with a number of other functions for
 performing set operations on arrays.

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Crosstabulation

[Zachary Pincus]

Hi Ionut,

Check out the "tabular" package:
http://parsemydata.com/tabular/index.html

It seems to be basically what you want... it does "pivot tables" (aka  
crosstabulation), it's built on top of numpy, and has simple data IO  
tools.

Also check out this discussion on "pivot tables" from the numpy list a  
while ago:
http://mail.scipy.org/pipermail/numpy-discussion/2007-August/028739.html

One question -- what do you mean by "raster data"? In my arena, that  
usually means "images"... and I'm not sure what a crosstabulation on  
image data would mean!

Zach


On Jul 14, 2010, at 2:28 PM, Ionut Sandric wrote:

>
> Sorry, the first email was sent before to finish it...
>
>
> Hi:
>
> I have two raster data and I would like to do a crosstabulation  
> between them and export the results to a table in a text file. Is it  
> possible to do it with NumPy? Does someone have an example?
>
> Thank you,
>
> Ionut
>
>
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] ndarray.resize that preserve view content ?

[Zachary Pincus]

> A[:5,:5] shows the data I want, but it's not contiguous in memory.
> A.resize(5,5) is contiguous, but do not contains the data I want.
>
> How to get both efficiently?

A[:5,:5].copy()
will give you a new, contiguous array that has the same contents as  
A[5:,5:], but in a new chunk of memory. Is this what you need?




On Aug 4, 2010, at 11:17 AM, Antoine Dechaume wrote:

> I forgot to refer to resize, sorry about that.
>
> A[:5,:5] shows the data I want, but it's not contiguous in memory.
> A.resize(5,5) is contiguous, but do not contains the data I want.
>
> How to get both efficiently?
>
>
> On Wed, Aug 4, 2010 at 5:01 PM, Robert Kern   
> wrote:
> On Wed, Aug 4, 2010 at 09:29, Antoine Dechaume   
> wrote:
> > Hi,
> >
> > given A=empty([10,10]), I would like to keep A[:5,:5] as a  
> contiguous memory
> > segment.
> >
> > How to do it efficiently?
>
> I'm not sure I understand what you want. Your Subject line and the
> body of your email conflict with each other. Can you try to explain
> what you want in other words?
>
> --
> Robert Kern
>
> "I have come to believe that the whole world is an enigma, a harmless
> enigma that is made terrible by our own mad attempt to interpret it as
> though it had an underlying truth."
>   -- Umberto Eco
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] ndarray.resize that preserve view content ?

[Zachary Pincus]

> Yes it is, but is there a way to do it in-place?

So you want the first 25 elements of the array (in a flat "contiguous"  
view) to contain the 25 elements of A[:5,:5]? This will do that, but  
having to do stuff like this (rather than just copying the memory  
region) might be indicative that maybe your code design might not  
really be right. (Why does it absolutely have to be in-place? Memory  
pressure?)

In [34]: a = numpy.array(numpy.arange(100, dtype=int).reshape(10,10))

In [35]: a
Out[35]:
array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
[10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
[20, 21, 22, 23, 24, 25, 26, 27, 28, 29],
[30, 31, 32, 33, 34, 35, 36, 37, 38, 39],
[40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
[50, 51, 52, 53, 54, 55, 56, 57, 58, 59],
[60, 61, 62, 63, 64, 65, 66, 67, 68, 69],
[70, 71, 72, 73, 74, 75, 76, 77, 78, 79],
[80, 81, 82, 83, 84, 85, 86, 87, 88, 89],
[90, 91, 92, 93, 94, 95, 96, 97, 98, 99]])

In [36]: a[:5,:5]
Out[36]:
array([[ 0,  1,  2,  3,  4],
[10, 11, 12, 13, 14],
[20, 21, 22, 23, 24],
[30, 31, 32, 33, 34],
[40, 41, 42, 43, 44]])

In [37]: a.flat[:25] = a[:5,:5]

In [38]: a.resize((5,5))
---
ValueErrorTraceback (most recent call  
last)

/Users/zpincus/ in ()

ValueError: cannot resize an array that has been referenced or is  
referencing
another array in this way.  Use the resize function

In [39]: b = numpy.ndarray(buffer=a, shape=(5,5), dtype=int)
In [40]: b

Out[40]:
array([[ 0,  1,  2,  3,  4],
[10, 11, 12, 13, 14],
[20, 21, 22, 23, 24],
[30, 31, 32, 33, 34],
[40, 41, 42, 43, 44]])


In [41]: b.flags.c_contiguous
Out[41]: True

In [42]: b.flags.owndata
Out[42]: False


Zach


> On Wed, Aug 4, 2010 at 5:20 PM, Zachary Pincus  > wrote:
> > A[:5,:5] shows the data I want, but it's not contiguous in memory.
> > A.resize(5,5) is contiguous, but do not contains the data I want.
> >
> > How to get both efficiently?
>
> A[:5,:5].copy()
> will give you a new, contiguous array that has the same contents as
> A[5:,5:], but in a new chunk of memory. Is this what you need?
>
>
>
>
> On Aug 4, 2010, at 11:17 AM, Antoine Dechaume wrote:
>
> > I forgot to refer to resize, sorry about that.
> >
> > A[:5,:5] shows the data I want, but it's not contiguous in memory.
> > A.resize(5,5) is contiguous, but do not contains the data I want.
> >
> > How to get both efficiently?
> >
> >
> > On Wed, Aug 4, 2010 at 5:01 PM, Robert Kern 
> > wrote:
> > On Wed, Aug 4, 2010 at 09:29, Antoine Dechaume 
> > wrote:
> > > Hi,
> > >
> > > given A=empty([10,10]), I would like to keep A[:5,:5] as a
> > contiguous memory
> > > segment.
> > >
> > > How to do it efficiently?
> >
> > I'm not sure I understand what you want. Your Subject line and the
> > body of your email conflict with each other. Can you try to explain
> > what you want in other words?
> >
> > --
> > Robert Kern
> >
> > "I have come to believe that the whole world is an enigma, a  
> harmless
> > enigma that is made terrible by our own mad attempt to interpret  
> it as
> > though it had an underlying truth."
> >   -- Umberto Eco
> > ___
> > NumPy-Discussion mailing list
> > NumPy-Discussion@scipy.org
> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >
> > ___
> > NumPy-Discussion mailing list
> > NumPy-Discussion@scipy.org
> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] ndarray.resize that preserve view content ?

[Zachary Pincus]

Oh and PS. Robert's right that there's no general way to do this! What  
I have only works because the data existing in the first 25 elements  
of A that get clobbered by the copy operation aren't the same data  
that are being copied (or where they are, the new copy is identical to  
the old one). Other slices won't have this property... A[:] = A[::-1]  
e.g. will fail totally.



On Aug 4, 2010, at 11:52 AM, Zachary Pincus wrote:

>> Yes it is, but is there a way to do it in-place?
>
> So you want the first 25 elements of the array (in a flat "contiguous"
> view) to contain the 25 elements of A[:5,:5]? This will do that, but
> having to do stuff like this (rather than just copying the memory
> region) might be indicative that maybe your code design might not
> really be right. (Why does it absolutely have to be in-place? Memory
> pressure?)
>
> In [34]: a = numpy.array(numpy.arange(100, dtype=int).reshape(10,10))
>
> In [35]: a
> Out[35]:
> array([[ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9],
>[10, 11, 12, 13, 14, 15, 16, 17, 18, 19],
>[20, 21, 22, 23, 24, 25, 26, 27, 28, 29],
>[30, 31, 32, 33, 34, 35, 36, 37, 38, 39],
>[40, 41, 42, 43, 44, 45, 46, 47, 48, 49],
>[50, 51, 52, 53, 54, 55, 56, 57, 58, 59],
>[60, 61, 62, 63, 64, 65, 66, 67, 68, 69],
>[70, 71, 72, 73, 74, 75, 76, 77, 78, 79],
>[80, 81, 82, 83, 84, 85, 86, 87, 88, 89],
>[90, 91, 92, 93, 94, 95, 96, 97, 98, 99]])
>
> In [36]: a[:5,:5]
> Out[36]:
> array([[ 0,  1,  2,  3,  4],
>[10, 11, 12, 13, 14],
>[20, 21, 22, 23, 24],
>[30, 31, 32, 33, 34],
>[40, 41, 42, 43, 44]])
>
> In [37]: a.flat[:25] = a[:5,:5]
>
> In [38]: a.resize((5,5))
> ---
> ValueErrorTraceback (most recent call
> last)
>
> /Users/zpincus/ in ()
>
> ValueError: cannot resize an array that has been referenced or is
> referencing
> another array in this way.  Use the resize function
>
> In [39]: b = numpy.ndarray(buffer=a, shape=(5,5), dtype=int)
> In [40]: b
>
> Out[40]:
> array([[ 0,  1,  2,  3,  4],
>[10, 11, 12, 13, 14],
>[20, 21, 22, 23, 24],
>[30, 31, 32, 33, 34],
>[40, 41, 42, 43, 44]])
>
>
> In [41]: b.flags.c_contiguous
> Out[41]: True
>
> In [42]: b.flags.owndata
> Out[42]: False
>
>
> Zach
>
>
>> On Wed, Aug 4, 2010 at 5:20 PM, Zachary Pincus >> wrote:
>>> A[:5,:5] shows the data I want, but it's not contiguous in memory.
>>> A.resize(5,5) is contiguous, but do not contains the data I want.
>>>
>>> How to get both efficiently?
>>
>> A[:5,:5].copy()
>> will give you a new, contiguous array that has the same contents as
>> A[5:,5:], but in a new chunk of memory. Is this what you need?
>>
>>
>>
>>
>> On Aug 4, 2010, at 11:17 AM, Antoine Dechaume wrote:
>>
>>> I forgot to refer to resize, sorry about that.
>>>
>>> A[:5,:5] shows the data I want, but it's not contiguous in memory.
>>> A.resize(5,5) is contiguous, but do not contains the data I want.
>>>
>>> How to get both efficiently?
>>>
>>>
>>> On Wed, Aug 4, 2010 at 5:01 PM, Robert Kern 
>>> wrote:
>>> On Wed, Aug 4, 2010 at 09:29, Antoine Dechaume 
>>> wrote:
>>>> Hi,
>>>>
>>>> given A=empty([10,10]), I would like to keep A[:5,:5] as a
>>> contiguous memory
>>>> segment.
>>>>
>>>> How to do it efficiently?
>>>
>>> I'm not sure I understand what you want. Your Subject line and the
>>> body of your email conflict with each other. Can you try to explain
>>> what you want in other words?
>>>
>>> --
>>> Robert Kern
>>>
>>> "I have come to believe that the whole world is an enigma, a
>> harmless
>>> enigma that is made terrible by our own mad attempt to interpret
>> it as
>>> though it had an underlying truth."
>>>  -- Umberto Eco
>>> ___
>>> NumPy-Discussion mailing list
>>> NumPy-Discussion@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>>
>>> ___
>>> NumPy-Discussion mailing list
>>> NumPy-Discussion@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Inverting argsort(a, axis=0) to obtain column-wise ranks

[Zachary Pincus]

> indices = argsort(a1)
> ranks = zeros_like(indices)
> ranks[indices] = arange(len(indices))

Doesn't answer your original question directly, but I only recently  
learned from this list that the following does the same as the above:
ranks = a1.argsort().argsort()
Will wonders never cease...

So does ranks=a2.argsort(axis=0).argsort(axis=0) then do the trick?

Zach

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] loadtxt stop

[Zachary Pincus]

>> In the end, the question was; is worth adding start= and stop=  
>> markers
>> into loadtxt to allow grabbing sections of a file between two known
>> headers?  I imagine it's something that people come up against  
>> regularly.


Simple enough to wrap your file in a new file-like object that stops  
coughing up lines when the delimiter is found, no?

class TruncatingFile(object):
   def __init__(self, fh, delimiter='END'):
 self.fh = fh
 self.delimiter=delimiter
 self.done = False
   def readline(self):
 if self.done: return ''
 line = self.fh.readline()
 if line.strip() == self.delimiter:
   self.done = True
   return ''
 return line
   def __iter__(self):
 return self
   def next(self):
 line = self.fh.next()
 if line.strip() == self.delimiter:
   self.done = True
   raise StopIteration()
 return line

from StringIO import StringIO
c = StringIO("0 1\n2 3\nEND")
numpy.loadtxt(TruncatingFile(c))

Though, really, it's annoying that numpy.loadtxt needs both the  
readline function *and* the iterator protocol. If it just used  
iterators, you could do:

def truncator(fh, delimiter='END'):
   for line in fh:
 if line.strip() == delimiter:
   break
 yield line

numpy.loadtxt(truncator(c))

Maybe I'll try to work up a patch for this.

Zach



On Sep 17, 2010, at 2:51 PM, Christopher Barker wrote:

> Neil Hodgson wrote:
>> In the end, the question was; is worth adding start= and stop=  
>> markers
>> into loadtxt to allow grabbing sections of a file between two known
>> headers?  I imagine it's something that people come up against  
>> regularly.
>
> maybe not so regular. However, a common use would be to be able load
> only n rows, which also does not appear to be supported. That would  
> be nice.
>
> -Chris
>
>
>
>> Thanks,
>> Neil
>>
>> 
>> *From:* Neil Hodgson 
>> *To:* numpy-discussion@scipy.org
>> *Sent:* Fri, 17 September, 2010 14:17:12
>> *Subject:* loadtxt stop
>>
>> Hi,
>>
>> I been looking around and could spot anything on this.  Quite often I
>> want to read a homogeneous block of data from within a file.  The
>> skiprows option is great for missing out the section before the data
>> starts, but if there is anything below then loadtxt will choke.  I
>> wondered if there was a possibility to put an endmarker= ?
>>
>> For example, if I want to load text from a large! file that looks  
>> like this
>>
>> header line
>> header line
>> 1 2.0 3.0
>> 2 4.5 5.7
>> ...
>> 500 4.3 5.4
>> END
>> more headers
>> more headers
>> 1 2.0 3.0 3.14 1.1414
>> 2 4.5 5.7 1.14 3.1459
>> ...
>> 500 4.3 5.4 0.000 0.001
>> END
>>
>> Then I can use skiprows=2, but loadtxt will choke when it gets to
>> 'END'.  To read t
>>
>>
>>
>> 
>>
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
> -- 
> Christopher Barker, Ph.D.
> Oceanographer
>
> Emergency Response Division
> NOAA/NOS/OR&R(206) 526-6959   voice
> 7600 Sand Point Way NE   (206) 526-6329   fax
> Seattle, WA  98115   (206) 526-6317   main reception
>
> chris.bar...@noaa.gov
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] loadtxt stop

[Zachary Pincus]

> Though, really, it's annoying that numpy.loadtxt needs both the
> readline function *and* the iterator protocol. If it just used
> iterators, you could do:
>
> def truncator(fh, delimiter='END'):
>   for line in fh:
> if line.strip() == delimiter:
>   break
> yield line
>
> numpy.loadtxt(truncator(c))
>
> Maybe I'll try to work up a patch for this.


That seemed easy... worth applying? Won't break compatibility, because  
the previous loadtxt required both fname.readline and fname.__iter__,  
while this requires only the latter.


Index: numpy/lib/npyio.py
===
--- numpy/lib/npyio.py  (revision 8716)
+++ numpy/lib/npyio.py  (working copy)
@@ -597,10 +597,11 @@
  fh = bz2.BZ2File(fname)
  else:
  fh = open(fname, 'U')
-elif hasattr(fname, 'readline'):
-fh = fname
  else:
-raise ValueError('fname must be a string or file handle')
+  try:
+  fh = iter(fname)
+  except:
+  raise ValueError('fname must be a string or file handle')
  X = []

  def flatten_dtype(dt):
@@ -633,14 +634,18 @@

  # Skip the first `skiprows` lines
  for i in xrange(skiprows):
-fh.readline()
+try:
+fh.next()
+except StopIteration:
+raise IOError('End-of-file reached before  
encountering data.')

  # Read until we find a line with some values, and use
  # it to estimate the number of columns, N.
  first_vals = None
  while not first_vals:
-first_line = fh.readline()
-if not first_line: # EOF reached
+try:
+first_line = fh.next()
+except StopIteration:
  raise IOError('End-of-file reached before  
encountering data.')
  first_vals = split_line(first_line)
  N = len(usecols or first_vals)

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] loadtxt stop

[Zachary Pincus]

On Sep 17, 2010, at 3:59 PM, Benjamin Root wrote:

> So, this code will still raise an error for an empty file.   
> Personally, I consider that a bug because I would expect to receive  
> an empty array.  I could understand raising an error for a non-empty  
> file that does not contain anything useful.  For comparison, Matlab  
> returns an empty matrix for loading an emtpy text file.
>
> This has been a long-standing annoyance for me, along with the  
> behavior with a single-line data file.

Agreed... I just wanted to make the patch as identical in behavior to  
the old version as possible. Though again, simple shims around loadtxt  
(as in my previous examples) can yield the desired behavior easily  
enough.

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] loadtxt stop

[Zachary Pincus]

>> Though, really, it's annoying that numpy.loadtxt needs both the
>> readline function *and* the iterator protocol. If it just used
>> iterators, you could do:
>>
>> def truncator(fh, delimiter='END'):
>>  for line in fh:
>>if line.strip() == delimiter:
>>  break
>>yield line
>>
>> numpy.loadtxt(truncator(c))
>>
>> Maybe I'll try to work up a patch for this.


http://projects.scipy.org/numpy/ticket/1616

Zach

>
>
> That seemed easy... worth applying? Won't break compatibility, because
> the previous loadtxt required both fname.readline and fname.__iter__,
> while this requires only the latter.
>
>
> Index: numpy/lib/npyio.py
> ===
> --- numpy/lib/npyio.py(revision 8716)
> +++ numpy/lib/npyio.py(working copy)
> @@ -597,10 +597,11 @@
>  fh = bz2.BZ2File(fname)
>  else:
>  fh = open(fname, 'U')
> -elif hasattr(fname, 'readline'):
> -fh = fname
>  else:
> -raise ValueError('fname must be a string or file handle')
> +  try:
> +  fh = iter(fname)
> +  except:
> +  raise ValueError('fname must be a string or file handle')
>  X = []
>
>  def flatten_dtype(dt):
> @@ -633,14 +634,18 @@
>
>  # Skip the first `skiprows` lines
>  for i in xrange(skiprows):
> -fh.readline()
> +try:
> +fh.next()
> +except StopIteration:
> +raise IOError('End-of-file reached before
> encountering data.')
>
>  # Read until we find a line with some values, and use
>  # it to estimate the number of columns, N.
>  first_vals = None
>  while not first_vals:
> -first_line = fh.readline()
> -if not first_line: # EOF reached
> +try:
> +first_line = fh.next()
> +except StopIteration:
>  raise IOError('End-of-file reached before
> encountering data.')
>  first_vals = split_line(first_line)
>  N = len(usecols or first_vals)
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] constructing an array from memory

[Zachary Pincus]

> I'm trying to do something ... unusual.
>
> gdb support scripting with Python. From within my python script, I can
> get the address of a contiguous area of memory that stores a  fortran
> array. I want to creat a NumPy array using "frombuffer". I see that
> the CPython API supports the creation of a buffer, but, is there an
> easier, more direct, way?

Here's how I do a similar task:
numpy.ndarray(shape, dtype=dtype,  
buffer=(ctypes.c_char*size_in_bytes).from_address(address))

You may need the strides or order parameters, as well.

Perhaps there's an easier way to create a buffer from an integer  
memory address, but this seems pretty straightforward.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] str/bytes object from arr.data?

[Zachary Pincus]

As str objects are supposed to be immutable, I think anything  
"official" that makes a string from a numpy array is supposed to copy  
the data. But I think you can use ctypes to wrap a pointer and a  
length as a python string.

Zach


On Sep 27, 2010, at 8:28 AM, Francesc Alted wrote:

> Hi,
>
> Anybody knows a way to get a str object (bytes for Python >= 2.6)  
> out of
> a buffer object (i.e. the .data attribute of ndarrays) without copying
> data?
>
> I need this for avoid creating a new function that deals with buffer
> objects (instead of reusing the one for str/byte objects that I  
> already
> have).  So, a matter of laziness :-)
>
> Thanks,
>
> -- 
> Francesc Alted
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] point line distance

[Zachary Pincus]

Here's a good list of basic geometry algorithms:
http://www.softsurfer.com/algorithms.htm

Zach


On Oct 6, 2010, at 5:08 PM, Renato Fabbri wrote:

> supose you have a line defined by two points and a point. you want  
> the distance
>
> what are easiest possibilities? i am doing it, but its nasty'n ugly
>
> -- 
> GNU/Linux User #479299
> skype: fabbri.renato
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] function name as parameter

[Zachary Pincus]

> I'm trying to write an implementation of the amoeba function from
> numerical recipes and need to be able to pass a function name and
> parameter list to be called from within the amoeba function.  Simply
> passing the name as a string doesn't work since python doesn't know it
> is a function and throws a typeerror.  Is there something similar to
> IDL's 'call_function' routine in python/numpy or a pythonic/numpy  
> means
> of passing function names?

Just pass the function itself! For example:

def foo():
   print 6

def call_function_repeatedly(func, count):
   for i in range(count):
 func()

call_function_repeatedly(foo, 2) # calls foo twice

bar = foo
bar() # still calls foo... we've just assigned the function to a  
different name


In python, functions (and classes, and everything else) are first- 
class objects and can be assigned to variables, passed around, etc,  
etc, just as anything else.

However, note that scipy.optimize.fmin implements the Nelder-Mead  
simplex algorithm, which is (I think) the same as the "amoeba"  
optimizer. Also you might be interested in the openopt package, which  
implements more optimizers a bit more consistently than scipy.optimize.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] function name as parameter

[Zachary Pincus]

> Hi Robert,
> so in a big data analysis framework, that is essentially written in C 
> ++,
> exposed to python with SWIG, plus dedicated python modules, the user
> performs an analysis choosing some given modules by name,as in :
> myOpt="foo"
> my_analyse.perform(use_optimizer=myOpt)
>
> The attribute use_optimizer is then checked to perform the right
> calls/instantiations of python but also C++ objects. and the  
> actual
> crunching number is in the C++ part.
> But then I realize that I need to tweak this optimizer's state, and  
> the
> optimizer object is accessible from a library pyOpt that has been
> swigified from a C++ library.
> Then I access the object by calling optObject =
> eval("pyOpt.%s(some_args)"%myOpt)
> where myOpt would be "foo" in this particular analysis. This is  
> because
> what the attribute use_optimizer expects is also the object name in  
> the
> library, of course.
> It goes without saying that I could do :
> if myOpt=="foo":
> optObject=pyOpt.foo(some_args)
> else:
> 
> and if you guys tell me it is way safer, I will do that instead of the
> use of eval that I liked because of the compactness.

Well, optObject=getattr(pyOpt, myOpt) is a bit nicer and marginally  
more likely to return sane error messages than optObject=eval("pyOpt. 
%s"%myOpt). Then you could do result=optObject(some_args)...

But why not have the user just pass in the relevant optObject from the  
pyOpt namespace (or some restricted namespace that just has the  
relevant functions exposed)? E.g.
my_analysis.perform(optimizer=pyOpt.Amoeba)
rather than
my_analysis.perform(optimizer='Amoeba')

This lets users do introspection on the pyOpt namespace to see what  
functions they can choose from, which is rather friendlier in an  
interactive environment.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] short circuit != ?

[Zachary Pincus]

This is silly: the structure of the python language prevents  
meaningful short-circuiting in the case of
np.any(a!=b)

While it's true that np.any itself may short-circuit, the 'a!=b'  
statement itself will be evaluated in its entirety before the result  
(a boolean array) is passed to np.any. This is just how python (like  
most other languages) works. This means that the big-O complexity  
isn't reduced any: all elements of a and b have to be examined. If any  
short-circuits, then at least all the elements don't need to be  
examined twice!

If performance actually matters here (unlikely?), a quick bit of  
cython code could do this easily enough in the simple case where a and  
b are the same shape (so no broadcasting required). Note that the  
below hard-codes the dtype also.

import cython
cimport numpy
import numpy

TYPE_T = numpy.int64_t
TYPE = numpy.int64

@cython.boundscheck(False)
def equal(a, b):
   cdef:
 numpy.ndarray[TYPE, ndim=1, negative_indices=False] a_flat =  
a.astype(TYPE).flatten()
 numpy.ndarray[TYPE, ndim=1, negative_indices=False] b_flat =  
b.astype(TYPE).flatten()
 unsigned int i, l
   assert a_flat.shape[0] == b_flat.shape[0]
   l = a_flat.shape[0]
   for i in range(l):
 if a_flat[i] != b_flat[i]:
   return False
   return True

Zach



On Oct 27, 2010, at 9:44 AM, Skipper Seabold wrote:

> On Wed, Oct 27, 2010 at 9:37 AM, Johann Cohen-Tanugi
>  wrote:
>>
>>
>> On 10/27/2010 03:31 PM, Neal Becker wrote:
>>> Johann Cohen-Tanugi wrote:
>>>
>>>
 how about np.any(a!=b)  ??

 On 10/27/2010 12:25 PM, Neal Becker wrote:

> Is there a way to get a short circuit != ?
>
> That is, compare 2 arrays but stop as soon as the first element
> comparison fails?
>
> I'm assuming that np.all (a != b) will _not_ do this, but will  
> first
> compare all elements.
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
>
>>> I don't think that will do short ciruit, will it?  I think that  
>>> will compare
>>> each element, returning a bool array, then short-circuit eval that  
>>> bool
>>> array.
>>>
>>>
>> In [3]: a=np.array([2,3,
>>
>> In [4]: b=np.array([2,5,
>>
>> In [5]: np.any(a!=b)
>> Out[5]: True
>>
>> it does not return a bool array, it seems. I do not see how you would
>> "broadcast" the notion of "any" along axes maybe?
>>
>
> Not definitive by any means, but
>
> In [28]: a = np.arange(1,1000)
>
> In [29]: b = np.arange(1,1000)
>
> In [30]: timeit np.any(a!=b)
> 10 loops, best of 3: 37.8 ms per loop
>
> In [31]: a[0] = 10.
>
> In [32]: timeit np.any(a!=b)
> 10 loops, best of 3: 24.7 ms per loop
>
> In [33]: a[0] = 1
>
> In [34]: a[-1] = 1
>
> In [35]: timeit np.any(a!=b)
> 10 loops, best of 3: 37.7 ms per loop
>
> It seems to at least take less time when the difference is at the
> "beginning," though I'm sure there could be exceptions.
>
> Skipper
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

[Zachary Pincus]

> Help! I'm having a problem in searching through the *elements* if a  
> 2d array. I have a loop over a numpy array:
>
> n,m = G.shape
> print n,m
> for i in xrange(n):
> for j in xrange(m):
> print type(G), type(G[i,j]), type(float(G[i,j]))
> g = float(abs(G[i,j]))
> if g < cut:
> print i,j,G[i,j]
>
> However, I'm getting the error message:
>
> The truth value of an array with more than one element is  
> ambiguous. Use a.any() or a.all()
>
> despite the fact that I'm not computing a truth value of an array.  
> I'd like to just compare to G[i,j] or abs(G[i,j]), which should be a  
> *single* value, and *not* an array. I even call 'float' here to make  
> sure that I cast this to a normal python float. But I still get this  
> error message.

Which line is raising the error?
My guess is it's the only truth-value testing line: "if g < cut". (I'm  
not sure what else could error here in that way...) It looks like g is  
probably a scalar, but your code isn't showing where cut comes from,  
nor have you printed out it's type... is it an array?

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] ndarray __getattr__ to perform __getitem__

[Zachary Pincus]

>>> But wouldn't the performance hit only come when I use it in this  
>>> way?
>>> __getattr__ is only called if the named attribute is *not* found (I
>>> guess it falls off the end of the case statement, or is the result  
>>> of
>>> the attribute hash table "miss").
>> That's why I said that __getattr__ would perhaps work better.
>
>
> So do you want me to try out an implementation and supply a patch?  If
> so, where should I send the patch?
>
> Ian

Note that there are various extant projects that I think attempt to  
provide similar functionality to what you're wanting (unless I badly  
misread your original email, in which case apologies):
http://projects.scipy.org/numpy/wiki/NdarrayWithNamedAxes

You might want to check these out (or pitch in with those efforts)  
before starting up your own variant. (Though if your idea has more  
modest goals, then it might be a good complement to these more  
extensive/intrusive solutions.)

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] N dimensional dichotomy optimization

[Zachary Pincus]

On Nov 23, 2010, at 10:57 AM, Gael Varoquaux wrote:

> On Tue, Nov 23, 2010 at 04:33:00PM +0100, Sebastian Walter wrote:
>> At first glance it looks as if a relaxation is simply not possible:
>> either there are additional rows or not.
>> But with some technical transformations it is possible to reformulate
>> the problem into a form that allows the relaxation of the integer
>> constraint in a natural way.
>
>> Maybe this is also possible in your case?
>
> Well, given that it is a cross-validation score that I am optimizing,
> there is not simple algorithm giving this score, so it's not obvious  
> at
> all that there is a possible relaxation. A road to follow would be to
> find an oracle giving empirical risk after estimation of the penalized
> problem, and try to relax this oracle. That's two steps further than  
> I am
> (I apologize if the above paragraph is incomprehensible, I am  
> getting too
> much in the technivalities of my problem.
>
>> Otherwise, well, let me know if you find a working solution ;)
>
> Nelder-Mead seems to be working fine, so far. It will take a few weeks
> (or more) to have a real insight on what works and what doesn't.

Jumping in a little late, but it seems that simulated annealing might  
be a decent method here: take random steps (drawing from a  
distribution of integer step sizes), reject steps that fall outside  
the fitting range, and accept steps according to the standard  
annealing formula.

Something with a global optimum but spikes along the way is pretty  
well-suited to SA in general, and it's also an easy algorithm to make  
work on a lattice. If you're in high dimensions, there are also bolt- 
on methods for biasing the steps toward "good directions" as opposed  
to just taking isotropic random steps. Again, pretty easy to think of  
discrete implementations of this...

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Threshold

[Zachary Pincus]

> mask = numpy.zeros(medical_image.shape, dtype="uint16")
> mask[ numpy.logical_and( medical_image >= lower, medical_image <=  
> upper)] = 255
>
> Where lower and upper are the threshold bounds. Here I' m marking the
> array positions where medical_image is between the threshold bounds
> with 255, where isn' t with 0. The question is: Is there a better  
> way to do that?

This will give you a True/False boolean mask:
mask = numpy.logical_and( medical_image >= lower, medical_image <=  
upper)

And this a 0/255 mask:
mask = 255*numpy.logical_and( medical_image >= lower, medical_image <=  
upper)

You can make the code a bit more terse/idiomatic by using the bitwise  
operators, which do logical operations on boolean arrays:
mask = 255*((medical_image >= lower) & (medical_image <= upper))

Though this is a bit annoying as the bitwise ops (& | ^ ~) have higher  
precedence than the comparison ops (< <= > >=), so you need to  
parenthesize carefully, as above.

Zach


On Dec 2, 2010, at 7:35 AM, totonixs...@gmail.com wrote:

> Hi all,
>
> I' m developing a medical software named InVesalius [1], it is a free
> software. It uses numpy arrays to store the medical images (CT and
> MRI) and the mask, the mask is used to mark the region of interest and
> to create 3D surfaces. Those array generally have 512x512 elements.
> The mask is created based in threshold, with lower and upper bound,
> this way:
>
> mask = numpy.zeros(medical_image.shape, dtype="uint16")
> mask[ numpy.logical_and( medical_image >= lower, medical_image <=  
> upper)] = 255
>
> Where lower and upper are the threshold bounds. Here I' m marking the
> array positions where medical_image is between the threshold bounds
> with 255, where isn' t with 0. The question is: Is there a better way
> to do that?
>
> Thank!
>
> [1] - svn.softwarepublico.gov.br/trac/invesalius
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Arrays with aliased elements?

[Zachary Pincus]

def repeat(arr, num):
   arr = numpy.asarray(arr)
   return numpy.ndarray(arr.shape+(num,), dtype=arr.dtype,
 buffer=arr, strides=arr.strides+(0,))

There are limits to what these sort of stride tricks can accomplish,  
but repeating as above, or similar, is feasible.


On Jan 1, 2011, at 8:42 PM, Enzo Michelangeli wrote:

> Is there any way, not involving compilation of C code, to define  
> ndarrays
> where some rows or columns share the same data buffers? For example,
> something built by a hypothetical variant of the np.repeat()  
> function, such
> that, if a = array([2,3]), calling:
>
>   b = np.aliasedrepeat(x, [1, 2], axis=0)
>
> would return in b:
>
>   array([[2, 3],
>  [2, 3],
>  [2, 3]])
>
> ...with the understanding that the three rows would actually share  
> the same
> data, so setting e.g.:
>
>   b[0,1] = 5
>
> ...would change b into:
>
>   array([[2, 5],
>  [2, 5],
>  [2, 5]])
>
> In other words, something with a behaviour similar to a list of lists:
>
 a = [2,3]
 b = [a,a,a]
 b
> [[2, 3], [2, 3], [2, 3]]
 b[0][1] = 5
 b
> [[2, 5], [2, 5], [2, 5]]
>
> This would save memory (and time spent in unnecessary copying) in some
> applications with large arrays, and would allow to cope with the  
> current
> inability of weave.blitz to understand broadcasting rules, e.g. for
> calculating outer products (I mentioned this in a previous thread).
>
> Enzo
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Variable in an array name?

[Zachary Pincus]

> Is it possible to use a variable in an array name?  I am looping  
> through a
> bunch of calculations, and need to have each array as a separate  
> entity.
> I'm pretty new to python and numpy, so forgive my ignorance.  I'm  
> sure there
> is a simple answer, but I can't seem to find it.
>
> let's say i have a variable 'i':
>
> i = 5
>
> I would like my array to have the name array5
>
> I know how I could do this manually, but not in a loop where i is  
> redefined
> several times.

There are ways to do this, but what you likely actually want is just  
to put several arrays in a python list and then index into the list,  
instead of constructing numbered names.

e.g.:

array_list = []

for whatever:
   array_list.append(numpy.array(whatever))

for array in array_list:
   do_something(array)

given_array = array_list[i]
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Variable in an array name?

[Zachary Pincus]

> Thank you very much for the prompt response.  I have already done  
> what you
> have suggested, but there are a few cases where I do need to have an  
> array
> named with a variable (looping through large numbers of unrelated  
> files and
> calculations that need to be dumped into different analyses).  It  
> would be
> extraordinarily helpful if someone could post a solution to this  
> problem,
> regardless of inefficiency of the method.  Thanks a ton for any  
> additional
> help.

You could store arrays associated with string names, or other  
identifiers, (as opposed to integer indices) in a python dict.

Global and local namespaces are also just dicts that you can grab with  
globals() and locals(), if you really want to look up variable names  
algorithmically, but I promise you that this is really not what you  
want to be doing.

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] NOOB Alert: Looping Through Text Files...

[Zachary Pincus]

> textlist = ["test1.txt", "test2.txt", "test3.txt"]
>
> for i in textlist:
>   text_file = open(textlist, "a")
>   text_file.write("\nI suck at Python and need help")
>   text_file.close()
>
> But, this doesn't work.  It gives me the error:
>
> coercing to Unicode: need string or buffer, list found

Yeah, it's probably the wrong list; still, easy enough to answer...

You want this:
   text_file = open(i, "a")

to open the individual file, as opposed to:
   text_file = open(textlist, "a")

which tries to open the whole list of filenames. Python cannot figure  
out how to turn the list into a single (unicode) filename, hence the  
error.

As for your wanting to write files with new names:

for txtfile in txtlist:
   f = open(txtfile, 'r')
   data = f.read()
   f.close()
   fnew = open(get_new_name(txtfile), 'w')
   fnew.write(data)
   fnew.write('\nHelp is on the way.')
   fnew.close()

where get_new_name() or whatever is defined appropriately to transform  
the old name ('test1.txt', say) into 'test1_appended.txt'...

Zach


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] intersect1d and setmember1d

[Zachary Pincus]

Hi,

> intersect1d and setmember1d doesn't give expected results in case  
> there are duplicate values in either array becuase it works by  
> sorting data and substracting previous value. Is there an  
> alternative in numpy to get indices of intersected values.

 From the docstring for setmember1d (and other set operations), you  
are only supposed to pass it arrays with unique values (i.e. arrays  
that represent sets in the mathematical sense):

 >>> print numpy.setmember1d.__doc__
 Return a boolean array set True where first element is in second  
array.

 Boolean array is the shape of `ar1` containing True where the  
elements
 of `ar1` are in `ar2` and False otherwise.

 Use unique1d() to generate arrays with only unique elements to  
use as
 inputs to this function.
 [...]

As stated, use unique1d to generate set-arrays from your input.

On the other hand, intersect1d is supposed to work with repeated  
elements:
 >>> print numpy.intersect1d.__doc__
 Intersection returning repeated or unique elements common to both  
arrays.

 Parameters
 --
 ar1,ar2 : array_like
 Input arrays.

 Returns
 ---
 out : ndarray, shape(N,)
 Sorted 1D array of common elements with repeating elements.

 See Also
 
 intersect1d_nu : Returns only unique common elements.
 [...]

Do you have an example of intersect1d not working right? If so, what  
version of numpy are you using (print numpy.version.version)?

Zach




On Feb 26, 2009, at 12:48 PM, mudit sharma wrote:

>
> intersect1d and setmember1d doesn't give expected results in case  
> there are duplicate values in either array becuase it works by  
> sorting data and substracting previous value. Is there an  
> alternative in numpy to get indices of intersected values.
>
> In [31]: p nonzero(setmember1d(v1.Id, v2.Id))[0]
> [ 0  1  2  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22  
> 23 24 25<-- index 2 shouldn't be  
> here look at the data below.
> 26 27 28 29]
>
> In [32]: p v1.Id[:10]
> [ 232.  232.  233.  233.  234.  234.  235.  235.  237.  237.]
>
> In [33]: p v2.Id[:10]
> [ 232.  232.  234.  234.  235.  235.  236.  236.  237.  237.]
>
> In [34]: p setmember1d(v1.Id, v2.Id)
> [ True  True  True False  True  True  True  True  True  True  True   
> True   <-- index 2 shouldn't be True
>  True  True  True  True  True  True  True  True  True  True  True   
> True
>  True  True  True  True  True  True]
>
> In [35]: p setmember1d(v1.Id[:10], v2.Id[:10])
> [ True  True  True False  True  True  True  True  True  True]
>
>
>
>
> ___
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] image processing using numpy-scipy?

[Zachary Pincus]

Hello,

>> This a little wiered problem. I am having a black and white image.  
>> (black
>> background)
>> Entire image is filled with noisy white patterns of different size  
>> and
>> shape. I need to fill the
>> white patches if there area is more then given one. Logically this  
>> could
>> possible to use a quickfill algorithm
>> for every pixel and if the filled area generated by that pixel is  
>> more then
>> given area(in pixel) then paint that
>> patch (filled area) with black.
>>
>> I have read the docs of PIL but there is no function for this. Can  
>> I use
>> numpy-scipy for the matter?
>
> Sure, there are a couple of options.  First, look at scipy.ndimage if
> there is anything there you can use (i.e. maybe binary dilation is
> sufficient).
>
> Otherwise, I've got some connected component code at:
>
> http://mentat.za.net/source/connected_components.tar.bz2
>
> (The repository is at http://mentat.za.net/hg/ccomp if you prefer)

I think that scipy.ndimage.label also provides connected-component  
labeling, with arbitrary connectivity. But this problem could also be  
solvable via dilation / erosion, if the area constraint is loose --  
e.g. if you want to quash all blobs smaller than, say, 5x5 pixels, you  
could just erode for 5 iterations, and then dilate the eroded image  
back for '-1' iterations (which causes the ndimage algorithm to  
iterate until no pixels change), using a mask of the original image  
(so that no pixels outside of the original blobs are turned on). This  
basically means that any object that survives the original erosion  
will be dilated back to its initial size. (Similar tricks can also be  
used to find / kill objects touching the edges -- use a non-zero  
constant out-of-bounds value and dilate an all-zeros array for -1  
iterations, using the original array as the mask. Then only objects  
touching the edge will get filled in...)

Alternately, if you need a very stringent area threshold -- e.g.  
remove all objects with five or fewer total pixels, regardless of  
their configuration -- then the connected-components approach is  
required. Below is a stab at it, though note that there's a slow step  
that I can't right now figure out how to avoid, short of coding just  
that in C or with cython or something...

Zach



import numpy
import scipy.ndimage as ndimage

_4_connected = numpy.array([[0, 1, 0],
 [1, 1, 1],
 [0, 1, 0]], dtype=bool)

def kill_small(binary_image, min_size, structure=_4_connected):
   label_image, num_objects = ndimage.label(binary_image, structure)
   # Label 0 is the background...
   object_labels = numpy.arange(1, num_objects+1)
   # Get the area of each labeled object by summing up the pixel  
values in the
   # binary image. (Divide the binary image by its max val to ensure  
that the image
   # consists of 1s and 0s, so that the sum equals the area in pixels.)
   areas = numpy.array(nd.sum(binary_image / binary_image.max(),  
label_image, object_labels))
   big_objects = object_labels[areas >= min_size]
   # This part will be pretty slow! But I can't think right now how to  
speed it up.
   # (If there are more big objects than small objects, it would be  
faster to
   # reverse this process and xor out the small objects from the  
binary image,
   # rather than the below, which or's up a new image of just the  
large objects.)
   big_object_image = numpy.zeros(binary_image.shape, dtype=bool)
   for bo in big_objects:
 big_object_image |= label_image == bo
   return big_object_image
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Bilateral filter

[Zachary Pincus]

Hi all,

I just grabbed the latest bilateral filter from Stéfan's repository,  
but I can't get it to work! I'm using a recent numpy SVN and the  
latest release of cython...

In [10]: bl = bilateral.bilateral(image, 2, 150)
---
NameError Traceback (most recent call  
last)

/Users/zpincus/Documents/Research/Slack Lab/Experimental Data/Big Data/ 
2009-2-17/ in ()

/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site- 
packages/bilateral/bilateral.pyc in bilateral(mat, xy_sig, z_sig,  
filter_size, mode)
  36 '''
  37
---> 38 filter_fcn = _BB.Bilat_fcn(xy_sig, z_sig, filter_size)
  39 size = filter_fcn.xy_size
  40 return generic_filter(mat, filter_fcn.cfilter, size=size,  
mode=mode)

/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site- 
packages/bilateral/bilateral_base.so in  
bilateral.bilateral_base.Bilat_fcn.__init__ (bilateral/ 
bilateral_base.c:590)()

NameError: np

I know very little about the internals of cython, so I can't figure  
this error out... it seems to be a problem with the 'cimport numpy as  
np' line in the pyx file, but beyond that I'm lost.


Let me know if I can provide any additional information,

Zach




On Feb 14, 2009, at 1:29 PM, Stéfan van der Walt wrote:

> Hi Nadav
>
> 2009/2/13 Nadav Horesh :
>> I tried the installer and it did not copy bilateral.py. I tried to  
>> improve it and the result is in the attached file.  I hope it would  
>> pass the mail filter, if not contact me directly to the email  
>> address below.
>
> Thanks!  I applied your changes and modified the setup.py to support
> in-place building.  Again available here:
>
> http://github.com/stefanv/bilateral/tree/master
>
> Cheers
> Stéfan
> ___
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Bilateral filter

[Zachary Pincus]

Hi all,

So, I re-grabbed the latest bilateral code from the repository, and  
did the following:

[axlotl:~/Developer/Python/stefanv-bilateral] zpincus% python setup.py  
build_ext -i
[axlotl:~/Developer/Python/stefanv-bilateral] zpincus% ipython
In [1]: import numpy
In [2]: import bilateral
In [3]: bilateral.bilateral(numpy.ones((10,10)), 4, 4)
---
NameError Traceback (most recent call  
last)

/Users/zpincus/Developer/Python/stefanv-bilateral/ in  
()

/Users/zpincus/Developer/Python/stefanv-bilateral/bilateral/ 
bilateral.pyc in bilateral(mat, xy_sig, z_sig, filter_size, mode)
  36 '''
  37
---> 38 filter_fcn = _BB.Bilat_fcn(xy_sig, z_sig, filter_size)
  39 size = filter_fcn.xy_size
  40 return generic_filter(mat, filter_fcn.cfilter, size=size,  
mode=mode)

/Users/zpincus/Developer/Python/stefanv-bilateral/bilateral/ 
bilateral_base.so in bilateral.bilateral_base.Bilat_fcn.__init__  
(bilateral/bilateral_base.c:590)()

NameError: np

In [4]: print numpy.version.version
1.3.0.dev6359
In [5]: import Cython.Compiler.Version
In [6]: print Cython.Compiler.Version.version
0.10.3

I'm using Python 2.5.2 on OS X 10.5.6. Oddly, when I run the cython  
tests, numpy_test fails to even compile because it can't find 'numpy/ 
arrayobject.h' -- I assume this is some problem with their test  
harness. However, I noticed that the numpy_test.pyx file does both:
'cimport numpy as np'
and
'import numpy as np'

so I added the second line to the bilateral_base.pyx file... re- 
running the above then gets a bit farther, but things choke again,  
with an "AttributeError: 'numpy.ndarray' object has no attribute  
'dimensions'", which comes from bilateral_base.c:1060, i.e.  
bilateral_base.pyx line 73.

This is perplexing. I wonder if the problem is using a dev numpy  
version, but a few-months-old cython. (Unfortunately, I don't have  
mercurial installed, so I can't easily grab the dev cython -- there's  
no simple way to download a mercurial repository over HTTP with wget  
or whatnot, which is annoying.)

Any thoughts? In the mean time, I guess I'll install mercurial so I  
can try the latest cython.

Zach


On Feb 28, 2009, at 1:53 PM, Nadav Horesh wrote:

> I could not reproduce your error. I am using latest official  
> releases of numpy and cython on a linux box (do you use Mac?).  I am  
> attaching the package I have on my PC, for the small chance it would  
> help.
>
>  Nadav.
>
>
> -הודעה מקורית-
> מאת: numpy-discussion-boun...@scipy.org בשם Zachary Pincus
> נשלח: ו 27-פברואר-09 22:26
> אל: Discussion of Numerical Python
> נושא: Re: [Numpy-discussion] Bilateral filter
>
> Hi all,
>
> I just grabbed the latest bilateral filter from St?fan's repository,
> but I can't get it to work! I'm using a recent numpy SVN and the
> latest release of cython...
>
> In [10]: bl = bilateral.bilateral(image, 2, 150)
> ---
> NameError Traceback (most recent call
> last)
>
> /Users/zpincus/Documents/Research/Slack Lab/Experimental Data/Big  
> Data/
> 2009-2-17/ in ()
>
> /Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
> packages/bilateral/bilateral.pyc in bilateral(mat, xy_sig, z_sig,
> filter_size, mode)
>  36 '''
>  37
> ---> 38 filter_fcn = _BB.Bilat_fcn(xy_sig, z_sig, filter_size)
>  39 size = filter_fcn.xy_size
>  40 return generic_filter(mat, filter_fcn.cfilter, size=size,
> mode=mode)
>
> /Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
> packages/bilateral/bilateral_base.so in
> bilateral.bilateral_base.Bilat_fcn.__init__ (bilateral/
> bilateral_base.c:590)()
>
> NameError: np
>
> I know very little about the internals of cython, so I can't figure
> this error out... it seems to be a problem with the 'cimport numpy as
> np' line in the pyx file, but beyond that I'm lost.
>
>
> Let me know if I can provide any additional information,
>
> Zach
>
>
>
>
> On Feb 14, 2009, at 1:29 PM, St?fan van der Walt wrote:
>
>> Hi Nadav
>>
>> 2009/2/13 Nadav Horesh :
>>> I tried the installer and it did not copy bilateral.py. I tried to
>>> improve it and the result is in the attached file.  I hope it would
>>> pass the mail filter, if not contact me directly to the email
>>> address below.
>>
>> Thanks!  I applied your changes and modified the setup.py to support
>> in-place building.  Again available here:
>

Re: [Numpy-discussion] Bilateral filter

[Zachary Pincus]

Well, the latest cython doesn't help -- both errors still appear as  
below. (Also, the latest cython can't run the numpy tests either.) I'm  
befuddled.

Zach

On Feb 28, 2009, at 5:18 PM, Zachary Pincus wrote:

> Hi all,
>
> So, I re-grabbed the latest bilateral code from the repository, and
> did the following:
>
> [axlotl:~/Developer/Python/stefanv-bilateral] zpincus% python setup.py
> build_ext -i
> [axlotl:~/Developer/Python/stefanv-bilateral] zpincus% ipython
> In [1]: import numpy
> In [2]: import bilateral
> In [3]: bilateral.bilateral(numpy.ones((10,10)), 4, 4)
> ---
> NameError Traceback (most recent call
> last)
>
> /Users/zpincus/Developer/Python/stefanv-bilateral/ in
> ()
>
> /Users/zpincus/Developer/Python/stefanv-bilateral/bilateral/
> bilateral.pyc in bilateral(mat, xy_sig, z_sig, filter_size, mode)
>  36 '''
>  37
> ---> 38 filter_fcn = _BB.Bilat_fcn(xy_sig, z_sig, filter_size)
>  39 size = filter_fcn.xy_size
>  40 return generic_filter(mat, filter_fcn.cfilter, size=size,
> mode=mode)
>
> /Users/zpincus/Developer/Python/stefanv-bilateral/bilateral/
> bilateral_base.so in bilateral.bilateral_base.Bilat_fcn.__init__
> (bilateral/bilateral_base.c:590)()
>
> NameError: np
>
> In [4]: print numpy.version.version
> 1.3.0.dev6359
> In [5]: import Cython.Compiler.Version
> In [6]: print Cython.Compiler.Version.version
> 0.10.3
>
> I'm using Python 2.5.2 on OS X 10.5.6. Oddly, when I run the cython
> tests, numpy_test fails to even compile because it can't find 'numpy/
> arrayobject.h' -- I assume this is some problem with their test
> harness. However, I noticed that the numpy_test.pyx file does both:
> 'cimport numpy as np'
> and
> 'import numpy as np'
>
> so I added the second line to the bilateral_base.pyx file... re-
> running the above then gets a bit farther, but things choke again,
> with an "AttributeError: 'numpy.ndarray' object has no attribute
> 'dimensions'", which comes from bilateral_base.c:1060, i.e.
> bilateral_base.pyx line 73.
>
> This is perplexing. I wonder if the problem is using a dev numpy
> version, but a few-months-old cython. (Unfortunately, I don't have
> mercurial installed, so I can't easily grab the dev cython -- there's
> no simple way to download a mercurial repository over HTTP with wget
> or whatnot, which is annoying.)
>
> Any thoughts? In the mean time, I guess I'll install mercurial so I
> can try the latest cython.
>
> Zach
>
>
> On Feb 28, 2009, at 1:53 PM, Nadav Horesh wrote:
>
>> I could not reproduce your error. I am using latest official
>> releases of numpy and cython on a linux box (do you use Mac?).  I am
>> attaching the package I have on my PC, for the small chance it would
>> help.
>>
>> Nadav.
>>
>>
>> -הודעה מקורית-
>> מאת: numpy-discussion-boun...@scipy.org בשם Zachary Pincus
>> נשלח: ו 27-פברואר-09 22:26
>> אל: Discussion of Numerical Python
>> נושא: Re: [Numpy-discussion] Bilateral filter
>>
>> Hi all,
>>
>> I just grabbed the latest bilateral filter from St?fan's repository,
>> but I can't get it to work! I'm using a recent numpy SVN and the
>> latest release of cython...
>>
>> In [10]: bl = bilateral.bilateral(image, 2, 150)
>> ---
>> NameError Traceback (most recent call
>> last)
>>
>> /Users/zpincus/Documents/Research/Slack Lab/Experimental Data/Big
>> Data/
>> 2009-2-17/ in ()
>>
>> /Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
>> packages/bilateral/bilateral.pyc in bilateral(mat, xy_sig, z_sig,
>> filter_size, mode)
>> 36 '''
>> 37
>> ---> 38 filter_fcn = _BB.Bilat_fcn(xy_sig, z_sig, filter_size)
>> 39 size = filter_fcn.xy_size
>> 40 return generic_filter(mat, filter_fcn.cfilter, size=size,
>> mode=mode)
>>
>> /Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
>> packages/bilateral/bilateral_base.so in
>> bilateral.bilateral_base.Bilat_fcn.__init__ (bilateral/
>> bilateral_base.c:590)()
>>
>> NameError: np
>>
>> I know very little about the internals of cython, so I can't figure
>> this error out... it seems to be a problem with the 'cimport numpy as
>> np' line

Re: [Numpy-discussion] Bilateral filter

[Zachary Pincus]

>> Well, the latest cython doesn't help -- both errors still appear as
>> below. (Also, the latest cython can't run the numpy tests either.)  
>> I'm
>> befuddled.
>
> That's pretty weird.  Did you remove the .so that was build as well as
> any source files, before doing build_ext with the latest Cython?  Also
> good to make sure that the latest Cython is, in fact, the one being
> used.

Yeah... and I just tried that again, with the same results. I have no  
idea what could be going wrong. E.g. why would 'cimport numpy as np'  
not add np to the namespace on my machine whereas it does so on yours...

Also, I assume that constructs like:
 cdef int i, dim = data.dimensions[0]
are some special numpy-support syntax that's supposed to be added by  
the cimport numpy line? (Because numpy arrays don't expose a  
'dimensions' attribute to python code...) It's like for some reason on  
my machine, cython isn't building its numpy support correctly. Which  
is understandable in the light of the fact that cython can't pass the  
numpy tests on my machine either. Odd indeed. Maybe I'll try on the  
cython list, since you guys seem to have demonstrated that the problem  
isn't in the bilateral code!

Zach
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Bilateral filter

[Zachary Pincus]

Hi guys,

Dag, the cython person who seems to deal with the numpy stuff, had  
this to say:
> - cimport and import are different things; you need both.
> - The "dimensions" field is in Cython renamed "shape" to be closer  
> to the Python interface. This is done in Cython/Includes/numpy.pxd

After including both the 'cimport' and 'import' lines, and changing  
'dimensions' to 'shape', things work perfectly. And now the shoe is on  
the other foot -- I wonder why you guys were getting proper results  
with what Dag claims to be buggy code! Odd indeed.

Zach






On Mar 1, 2009, at 2:36 PM, Nadav Horesh wrote:

> 1. "dimensions" is a field in the C struct, that describes the array  
> object.
> 2. Is there a chance that the header file numpy/arrayobject.h  
> belongs to a
>   different numpy version that you run?
>
> I am not very experienced with cython (I suppose that Stefan has  
> some experience).
> As you said, probably the cython list is a  better place to look for  
> an answer. I would be happy to see how this issue resolved.
>
>  Nadav
>
>
>
> -הודעה מקורית-
> מאת: numpy-discussion-boun...@scipy.org בשם Zachary Pincus
> נשלח: א 01-מרץ-09 20:59
> אל: Discussion of Numerical Python
> נושא: Re: [Numpy-discussion] Bilateral filter
>
>>> Well, the latest cython doesn't help -- both errors still appear as
>>> below. (Also, the latest cython can't run the numpy tests either.)
>>> I'm
>>> befuddled.
>>
>> That's pretty weird.  Did you remove the .so that was build as well  
>> as
>> any source files, before doing build_ext with the latest Cython?   
>> Also
>> good to make sure that the latest Cython is, in fact, the one being
>> used.
>
> Yeah... and I just tried that again, with the same results. I have no
> idea what could be going wrong. E.g. why would 'cimport numpy as np'
> not add np to the namespace on my machine whereas it does so on  
> yours...
>
> Also, I assume that constructs like:
> cdef int i, dim = data.dimensions[0]
> are some special numpy-support syntax that's supposed to be added by
> the cimport numpy line? (Because numpy arrays don't expose a
> 'dimensions' attribute to python code...) It's like for some reason on
> my machine, cython isn't building its numpy support correctly. Which
> is understandable in the light of the fact that cython can't pass the
> numpy tests on my machine either. Odd indeed. Maybe I'll try on the
> cython list, since you guys seem to have demonstrated that the problem
> isn't in the bilateral code!
>
> Zach
> ___
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
> ___
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Bilateral filter

[Zachary Pincus]

> 2009/3/1 Zachary Pincus :
>> Dag, the cython person who seems to deal with the numpy stuff, had
>> this to say:
>>> - cimport and import are different things; you need both.
>>> - The "dimensions" field is in Cython renamed "shape" to be closer
>>> to the Python interface. This is done in Cython/Includes/numpy.pxd
>
> Thanks for following up.  I made the fixes in:
>
> http://github.com/stefanv/bilateral.git

Cool! Does this, out of curiosity, break things for you? (Or Nadav?)

> I think we should combine all these image processing algorithms (and
> the ones you sent to the list) into an image processing scikit.  We've
> certainly got enough algorithms lying around!
>
> If you think that's a good idea, I'll set it up this coming week.

I'm all for it. I've got a few other bits lying around that might be  
good there too:
  - 2D iso-contour finding (sub-pixel precision)
  - 2D image warping via thin-plate splines

I also have some code for various geometric algorithms lying around:
  - calculating optimal rigid alignments of point-sets ("Procrustes  
Analysis")
  - line intersections, closest points to lines, distance to lines, etc.

if that would be of any use to anyone.

Zach
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Bilateral filter

[Zachary Pincus]

Hi Stéfan,

>>> http://github.com/stefanv/bilateral.git
>>
>> Cool! Does this, out of curiosity, break things for you? (Or Nadav?)
>
> I wish I had some way to test.  Do you maybe have a short example that
> I can convert to a test?

Here's my test case for basic working-ness (e.g. non exception- 
throwing) of that bilateral code:

In [7]: bilateral.bilateral(numpy.arange(25).reshape((5,5)), 4, 10)

Out[7]:
array([[ 7,  7,  7,  8,  8],
[ 9,  9,  9, 10, 10],
[11, 11, 12, 12, 12],
[13, 13, 14, 14, 14],
[15, 15, 16, 16, 16]])

That's all I'd been using to provoke the errors before, so presumably  
if you get that far with the fixed code, then things should be good as  
far as cython's concerned?

>> I'm all for it. I've got a few other bits lying around that might be
>> good there too:
>>  - 2D iso-contour finding (sub-pixel precision)
>>  - 2D image warping via thin-plate splines
>
>> I also have some code for various geometric algorithms lying around:
>>  - calculating optimal rigid alignments of point-sets ("Procrustes
>> Analysis")
>>  - line intersections, closest points to lines, distance to lines,  
>> etc.
>>
>> if that would be of any use to anyone.
>
> Definitely.  In addition I have code for polygon clipping, hough
> transforms, grey-level co-occurrence matrices, connected components,
> shortest paths and linear position-invariant filtering.

Aah, fantastic. The co-occurrence matrix stuff will be very useful to  
me!

Zach

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] image processing using numpy-scipy?

[Zachary Pincus]

> did you have a look at OpenCV?
>
> http://sourceforge.net/projects/opencvlibrary
>
> Since a couple of weeks, we have implemented the numpy array  
> interface so data exchange is easy [check out from SVN].

Oh fantastic! That is great news indeed.

Zach
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] 1.3.0 release notes erratum

[Zachary Pincus]

Hi David,

Thanks again for bundling in the architecture-specification flag into  
the numpy superpack installers: being able to choose sse vs. nosse  
installs is really helpful to me, and from what I hear, many others as  
well!

Anyhow, I just noticed (sorry I didn't see this before the release) a  
minor erratum in the release note wrt this feature, which reads:
> Automatic arch detection can now be bypassed from the command line  
> for the superpack installed:
>
> numpy-1.3.0-superpack-win32.exe /arch=nosse
>
> will install a numpy which works on any x86, even if the running  
> computer
> supports SSE set.
>

However, the actual command-line should be:
numpy-1.3.0-superpack-win32.exe /arch nosse

Probably anyone advanced enough to use this feature would figure that  
out, but no harm being clear...

Thanks,
Zach
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] replace voids in 2d dat with nearest neighbour value

[Zachary Pincus]

Hi Christian,

Check out this discussion from a little while ago on a very similar  
issue (but in 3d):
http://www.nabble.com/Help-with-interpolating-missing-values-from-a-3D-scanner-td21489413.html

Most of the suggestions should be directly applicable.

Zach


On Apr 6, 2009, at 9:01 AM, Christian K. wrote:

> Hi,
>
> I am looking for an elegant and fast way to fill the voids of a 2d  
> array with
> neighbouring values. The array's size can be up to (1000, 1000) and  
> its values
> are slowly varying around a mean value. What I call voids are values  
> which are
> far from the mean value (+- 80%). A void usually extends over some  
> adjacent
> coordinates.
>
> Currently I am using
>
> tmp = N.ma.array(tmp, tmp data[tmp.mask] = tmp.mean()
>
> which moves the voids closer to the mean value but which is still  
> far from
> beeing a smooth interpolation.
>
> Regards, Christian
>
>
> ___
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] shift for optimal superimposition of two 3D matrices according to correlation computed using FFT

[Zachary Pincus]

> I have two 3D density maps (meaning volumetric data, each one
> essentially a IxJxK matrix containing real values that sum to one) and
> want to find translation of between the two that maximises  
> correlation.
> This can be done by computing the correlation between the two
> (correlation theorem -> I use FFT for this) and then shifting one of  
> the
> two maps by the vector that corresponds to the index of the maximum in
> the 3D correlation matrix.
>
> My problem with this method is, however, that the shift wraps around  
> the
> edges of the matrix. e.g: for a 30x30x30 matrix that needs to be
> slightly shifted to fit the seconds 3D matrix, sometimes shift vectors
> like (0, 29, 0) appear, that actually mean a shift by (0, -1, 0). This
> means that the size of the matrix needs to be subtracted from some
> components of the shift vector in special cases.
>
> So far, I have been unable to determine the cutoff value, outside of
> which I need to subtract the matrix size from a component of the shift
> vector, as this cutoff seems to vary between different data sets.

Does it work to use a cutoff of half the size of the input arrays in  
each dimension? This is equivalent to calculating both shifts (the  
positive and negative) and using whichever has a smaller absolute value.

Alternately, you could use numpy.roll to shift the data (one axis at a  
time). Since roll wraps around, you wouldn't need to bother figuring  
out which shift is "correct".

Finally, you could not use FFTs but instead directly optimize a  
transformation between the two, using scipy.ndimage's affine  
transforms and scipy.optimize's numerical optimizers.

Zach

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] shift for optimal superimposition of two 3D matrices according to correlation computed using FFT

[Zachary Pincus]

>> Does it work to use a cutoff of half the size of the input arrays in
>> each dimension? This is equivalent to calculating both shifts (the
>> positive and negative) and using whichever has a smaller absolute  
>> value.
> no, unfortunately the cutoff is not half of the dimensions.

Explain more about what you need from this "cutoff" then. Given that  
the peak is just an index into the array, it seems like the "cutoff"  
for using the peak index as the shift vs. subtracting the size of the  
array from the index should be a function of the size of the arrays.  
How is the cutoff value data-dependent? Can you provide some simple 1D  
examples?

>> Alternately, you could use numpy.roll to shift the data (one axis  
>> at a
>> time). Since roll wraps around, you wouldn't need to bother figuring
>> out which shift is "correct".
> Maybe I do not understand correctly what you mean by this, but as  
> far as
> I get it this will not help me, given that I do not know the number of
> positions I need to roll. I think the parameter  to the roll
> function would be the same value as the cutoff I need.
> If you wanted to suggest that I can check some property for every  
> shift
> and hence could try every possible shift in each axis: this is
> computationally not feasable for me. I need to run time consuming
> operations to ensure the correctness and with a 200x200x200 matrix
> (which is very roughly as large as it gets for me), so this would ruin
> the speed benfits of using the FFT.

1D: say that you have 30-element arrays and a correlation peak of 25.  
I assume this means you need to shift one array 5 positions backward  
or the other 5 forward. However (and this ought to be justifiable  
since the FFT assumes wrap-around as well) you could also just shift  
one array 25 positions forward, wrapping around as needed. No?  
numpy.roll shifts arrays with wrap-around.

If this isn't what you need, you'll have to explain a bit more about  
your problem, and about what it is you need from the cutoff value.

Zach
___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] help getting started

[Zachary Pincus]

Hi John,

First, did you build your own Python 2.6 or install from a binary?  
When you type "python" at the command prompt, which python runs? (You  
can find this out by running "which python" from the command line.)

Second, it appears that numpy is *already installed* for a non-apple  
python 2.5 on your computer. (Either you manually built that python,  
or installed it from an installer, and then you built or installed  
numpy for that python?) I can tell this because you have a directory:

/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site- 
packages/numpy

which normally gets created during the process of installing numpy for  
Python 2.5. Unless you manually created this directory yourself  
somehow? Anyhow, that directory is *not* the source-code directory  
that you should be running numpy's setup.py from; it is the directory  
that the installation process for numpy copies all of the necessary  
files into!

Anyhow, how did you come to be executing commands from that directory?  
The numpy install instructions definitely don't have a "cd /Library/ 
Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/ 
numpy" step...

So, if you want to use numpy with python 2.5, you should be all set  
already. If you need to use it with 2.6, then you will want to  
download the source code into a directory somewhere -- say in your  
home folder. Then execute "python setup.py install" from *this*  
directory. (Note: if typing "python" does not run python 2.6, but  
instead runs python 2.5, you will want to do "python2.6 setup.py  
install"...

Feel free to ask any more questions -- I'm heading out the door, so  
this might be a bit terse.

Zach


On Apr 10, 2009, at 1:30 PM, John Seales wrote:

> Intel mac.
> Python 2.6.
>
> I'm trying to start using numpy and scipy, and am having trouble.
>
> I'm trying to install numpy following the instructions on 
> http://scipy.org/Installing_SciPy/Mac_OS_X 
>  .
>
> when I give the python setup.py build command, it doesn't work.  
> Here's my interaction at the terminal:
>
> john-seales-computer:/Library/Frameworks/Python.framework/Versions/ 
> 2.5/lib/python2.5/site-packages/numpy johnseales$ python setup.py  
> build
>
> This is the wrong setup.py file to run
> john-seales-computer:/Library/Frameworks/Python.framework/Versions/ 
> 2.5/lib/python2.5/site-packages/numpy
>
> Anybody have any ideas how I can get up and runnng?
>
>
>
> Rediscover Hotmail®: Now available on your iPhone or BlackBerry  
> Check it out.___
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Creating a RGB-image from BW

[Zachary Pincus]

Hi Johannes,

According to http://www.pygtk.org/pygtk2reference/class- 
gdkpixbuf.html , the pixels_array is a numeric python array (a  
predecessor to numpy). The upshot is that perhaps the nice  
broadcasting machinery will work fine:

pb_pixels[...] = fits_pixels[..., numpy.newaxis]

This might not work, though. But perhaps it would be possible to make  
a numpy array that's really just a view onto the memory of pb_pixels,  
or perhaps one could convert fits_pixels into a numeric array...  
Hopefully someone on the list can make a suggestion about dealing with  
numeric arrays.

Alternately, there are pixbuf methods for reading image data from  
strings. You'd just need to get fits_pixels set up properly, then call  
tostring() on it, and pass that to the pixbuf. The trick is in setting  
up fits_pixels so that its memory layout corresponds to what gtk  
wants. Usually, images are stored in memory as (r,g,b) tuples packed  
by rows and then columns; this is I assume what GTK wants. So you'd do  
something like:

fits_color = numpy.empty((height, width, 3), dtype=numpy.uint8)
fits_color[...] = fits_pixels[..., numpy.newaxis]
fits_string = fits_color.tostring()
pb = gtk.gdk.pixbuf_new_from_data(fits_string, gtk.gdk.COLORSPACE_RGB,  
False, 8, 640, 480, 3*640)

Zach

On Apr 28, 2009, at 2:36 AM, Johannes Bauer wrote:

> Hello group,
>
> I've been redicted from usenet ("Convert numpy.ndarray into "normal"
> array", <75dgm1f16hqn...@mid.dfncis.de>) here and hope this is the  
> right
> place.
>
> Basically, what I have is a numpy-Array which I got from a FITS-file
> (it's black/white). I want to display that using GTK. Therefore every
> color needs to appear three times (to make it look gray R = G = B).
>
> The basic framework looks like
>
> [...]
> pb = gtk.gdk.Pixbuf(gtk.gdk.COLORSPACE_RGB, False, 8, width, height)
> pb_pixels = pb.get_pixels_array()
>
>
> print(type(pb_pixels), pb_pixels.shape, pb_pixels.typecode())
> print(type(fits_pixels), fits_pixels.shape, fits_pixels.dtype)
>
> which gives
>
> (, (480, 640, 3), 'b')
> (, (480, 640), dtype('uint8'))
>
> so now I need to assign values. Naively I started out with
>
> for x in range(width):
>   for y in range(height):
>   pb_pixels[y, x] = fits_pixels[y, x]
>
> which was horribly slow (around 3 seconds). Thanks to the usenet  
> help, I
> now got somewhat better:
>
> fits_colors = numpy.zeros((height, width, 3), dtype="uint8")
> for y in range(height):
>   for x in range(width):
>   fits_colors[height -  y - 1, x] = fits_pixels[y, x]
> pb_pixels[:, :] = fits_colors
>
> This also works, and is a lot faster (around 0.7 seconds). However,
> there seems to be a better way to do it. I played around with
>
> fits_colors = numpy.fromfunction(lambda y, x, z: fits_pixels[y, x],
> (height, width, 3), dtype="uint8")
> pb_pixels[:, :] = fits_colors
>
> Which worked somewhat - but gives weird results: The picture is
> rotatated 90° to the right and the lower left part is displayed
> repeatedly after 256 pixels... (I can make a screenshot if that's
> easier). The fromfunction Function is quite fast in my context (around
> 0.2 second).
>
> How should I solve this problem the right way?
>
> Kind regards,
> Johannes
> ___
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Efficient scaling of array

[Zachary Pincus]

scipy.ndimage.zoom (and related interpolation functions) would be a  
good bet -- different orders of interpolation are available, too,  
which can be useful.

Zach



On May 4, 2009, at 11:40 AM, Johannes Bauer wrote:

> Hello list,
>
> is there a possibility to scale an array by interpolation,
> automatically? For illustration a 1D-example would be an array of size
> 5, which is scaled to size 3:
>
> before: [ 1,   2,   3,   4,   5 ]
>  1/1  2/3
>   1/3  1   1/3
>2/3   1
> after : [ 2.33, 5, 7.66 ]
>
>
> The same thing should be possible in nD, with the obvious analogy. Is
> there such a function in numpy?
>
> Kind regards,
> Johannes
> ___
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] interpolation in numpy

[Zachary Pincus]

You might want also to look into scipy.ndimage.zoom.

Zach

On Jul 9, 2009, at 9:42 AM, Thomas Hrabe wrote:

> Hi all,
>
> I am not a newbie to python and numpy, but however, I kinda do not
> find a proper solution for my interpolation problem without coding it
> explicitly myself.
>
> All I want to do is to increase the resolution of an tree  
> dimensional array.
>
> I have a volume 'a'
>
> a = numpy.random.rand(3,3,3)
>
> now, I would like to expand a to a volume b of size say 10x10x10
>
> I want b to have interpolated values from a. One could think of this
> procedure as "zooming" into the volume, similar for images or so.
>
> numpy.interp does such a thing for 1d, but is there such function  
> for 3d, too?
>
> Thank you in advance for your help,
>
> Thomas
>
> FYI: I do not have scipy installed
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] mean of two or more arrays while ignoring a specific value

[Zachary Pincus]

Might want to look into masked arrays: numpy.ma.array.
a = numpy.array([1,5,4,99])
b = numpy.array([3,7,2,8])
arr = numpy.array([a, b])
masked = numpy.ma.array(arr, mask = arr==99)
masked.mean(axis=0)

masked_array(data = [2.0 6.0 3.0 8.0],
  mask = [False False False False],
fill_value = 1e+20)

Zach


On Jul 15, 2009, at 9:40 AM, Greg Fiske wrote:

> Dear list,
>
> I’m learning to work with numpy arrays.  Can somebody explain how to  
> get the average of two separate arrays while ignoring a user defined  
> value in one array?
>
> For example:
> >>>a = numpy.array([1,5,4,99])
> >>>b = numpy.array([3,7,2,8])
>
> Ignoring the value 99, the result should be an array like c=  
> ([2,6,3,8])
>
> Thanks for any advice,
>
> Greg
>
>
> Greg Fiske
> Research Associate/GIS Manager
> The Woods Hole Research Center
> 149 Woods Hole Road
> Falmouth MA, 02540
> 508-540-9900x139
> gfi...@whrc.org
> http://whrc.org
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Functions to pack/unpack bytes?

[Zachary Pincus]

> Does numpy have functions to convert between e.g. an array of uint32  
> and
> uint8, where the uint32 array is a packed version of the uint8 array
> (selecting little/big endian)?


You could use the ndarray constructor to look at the memory differently:

In : a = numpy.arange(240, 260, dtype=numpy.uint32)
In : a
Out:
array([240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252,
253, 254, 255, 256, 257, 258, 259], dtype=uint32)

In : b = numpy.ndarray(shape=(len(a)*4,), dtype=numpy.uint8, buffer=a)

In : b
Out:
array([240,   0,   0,   0, 241,   0,   0,   0, 242,   0,   0,   0, 243,
  0,   0,   0, 244,   0,   0,   0, 245,   0,   0,   0, 246,   0,
  0,   0, 247,   0,   0,   0, 248,   0,   0,   0, 249,   0,   0,
  0, 250,   0,   0,   0, 251,   0,   0,   0, 252,   0,   0,   0,
253,   0,   0,   0, 254,   0,   0,   0, 255,   0,   0,   0,   0,
  1,   0,   0,   1,   1,   0,   0,   2,   1,   0,   0,   3,   1,
  0,   0], dtype=uint8)


I assume for selecting little/big-endian going the other way, you  
could use the other methods of specifying dtypes that allow for byte- 
order descriptors. (Like dtype objects or the format strings.)

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Optimized half-sizing of images?

[Zachary Pincus]

> We have a need to to generate half-size version of RGB images as  
> quickly
> as possible.

How good do these need to look? You could just throw away every other  
pixel... image[::2, ::2].

Failing that, you could also try using ndimage's convolve routines to  
run a 2x2 box filter over the image, and then throw away half of the  
pixels. But this would be slower than optimal, because the kernel  
would be convolved over every pixel, not just the ones you intend to  
keep.

Really though, I'd just bite the bullet and write a C extension (or  
cython, whatever, an extension to work for a defined-dimensionality,  
defined-dtype array is pretty simple), or as suggested before, do it  
on the GPU. (Though I find that readback from the GPU can be slow  
enough that C code can beat it in some cases.)

Zach
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] [SciPy-User] Simple pattern recognition

[Zachary Pincus]

I believe that pretty generic connected-component finding is already  
available with scipy.ndimage.label, as David suggested at the  
beginning of the thread...

This function takes a binary array (e.g. zeros where the background  
is, non-zero where foreground is) and outputs an array where each  
connected component of non-background pixels has a unique non-zero  
"label" value.

ndimage.find_objects will then give slices (e.g. bounding boxes) for  
each labeled object (or a subset of them as specified). There are also  
a ton of statistics you can calculate based on the labeled objects --  
look at the entire ndimage.measurements namespace.

Zach

On Sep 21, 2009, at 1:45 PM, Gökhan Sever wrote:

> I asked this question at 
> http://stackoverflow.com/questions/1449139/simple-object-recognition 
>  and get lots of nice feedback, and finally I have managed to  
> implement what I wanted.
>
> What I was looking for is named "connected component labelling or  
> analysis" for my "connected component extraction"
>
> I have put the code (lab2.py) and the image (particles.png) under:
> http://code.google.com/p/ccnworks/source/browse/#svn/trunk/AtSc450/ 
> labs
>
> What do you think of improving that code and adding into scipy's  
> ndimage library (like connected_components())  ?
>
> Comments and suggestions are welcome :)
>
>
> On Wed, Sep 16, 2009 at 7:22 PM, Gökhan Sever  
>  wrote:
> Hello all,
>
> I want to be able to count predefined simple rectangle shapes on an  
> image as shown like in this one: 
> http://img7.imageshack.us/img7/2327/particles.png
>
> Which is in my case to count all the blue pixels (they are ice-snow  
> flake shadows in reality) in one of the column.
>
> What is the way to automate this task, which library or technique  
> should I study to tackle it.
>
> Thanks.
>
> -- 
> Gökhan
>
>
>
> -- 
> Gökhan
> ___
> SciPy-User mailing list
> scipy-u...@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] difficulty with numpy.where

[Zachary Pincus]

Hello,

a < b < c (or any equivalent expression) is python syntactic sugar for  
(a < b) and (b < c).

Now, for numpy arrays, a < b gives an array with boolean True or False  
where the elements of a are less than those of b. So this gives us two  
arrays that python now wants to "and" together. To do this, python  
tries to convert the array "a < b" to a single True or False value,  
and the array "b < c" to a single True or False value, which it then  
knows how to "and" together. Except that "a < b" could contain many  
True or False elements, so how to convert them to a single one?  
There's no obvious way to guess -- typically, one uses "any" or "all"  
to convert a boolean array to a single true or false value, depending,  
obviously, on what one needs.

So this explains the error you see, but has nothing to do with the  
results you desire... you need to and-together two boolean arrays  
*element-wise* -- which is something Python doesn't know how to do  
with the builtin "and" operator (which cannot be overridden). To do  
this, you need to use the bitwise logic operators:
(a < b) & (b < c).

So:

def sin_half_period(x): return where((0.0 <= x) & (x <= pi), sin(x),  
0.0)

Zach



On Oct 1, 2009, at 12:55 PM, Dr. Phillip M. Feldman wrote:

>
> I've defined the following one-line function that uses numpy.where:
>
> def sin_half_period(x): return where(0.0 <= x <= pi, sin(x), 0.0)
>
> When I try to use this function, I get an error message:
>
> In [4]: z=linspace(0,2*pi,9)
> In [5]: sin_half_period(z)
> ---
> ValueErrorTraceback (most recent  
> call last)
>
> The truth value of an array with more than one element is ambiguous.  
> Use
> a.any
> () or a.all()
>
> Any suggestions will be appreciated.
> -- 
> View this message in context: 
> http://www.nabble.com/difficulty-with-numpy.where-tp25702676p25702676.html
> Sent from the Numpy-discussion mailing list archive at Nabble.com.
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion