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

[Numpy-discussion] distutils env variables

[Charlie M. <cw...@gm...>]

So numpy's distutils highjacking doesn't seem to respect environment
variables.  When I try something like...

CPPFLAGS=3D"-I/extra/include/path" python setup.py build

... it seems ignored.  Is this intentional, or just missing functionality.

Thanks,
    Charlie

[Numpy-discussion] field method in recarray

[Christopher H. <ch...@st...>]

Hi Travis,

I have added a field method to recarray.  This allows field access via 
either field name or index number.

Example below:

In [1]: from numpy import rec

In [2]: r = 
rec.fromrecords([[456,'dbe',1.2],[2,'de',1.3]],names='col1,col2,col3')

In [3]: r.field('col1')
Out[3]: array([456,   2])

In [4]: r.field(0)
Out[4]: array([456,   2])

In [5]: r.field(0)[1]=1000

In [6]: r.field(0)
Out[6]: array([ 456, 1000])


Chris

[Numpy-discussion] compressed method in doc

[manouchk <man...@gm...>]

Hi,

First of all, I'm quite new to python and numpy (maybe to new!)... I was 
loooking for a convenient way to convert a (1D) masked array to a numeric 
array which only contain the not masked values.
I spent several hours to figure out that the method "compressed" was the exact 
thing I was looking for! 

So I'm wondering what is the common way to find the right method (if there is) 
for a beginner ? Look at help of all methods (using completion to find all 
methods with tab key)? I hope the question is not too much basic for numpy 
mailing-list!


Emmanuel

[Numpy-discussion] numpy/scipy transition

[Arnd B. <arn...@we...>]

Hi,

concerning the transition from Numeric/scipy to the new numpy/scipy
I have a couple of points for which I would be very interested
in any advice/suggestions:

As some of you might know we are running a computational physics
course using python+Numeric+scipy+ipython.
In six weeks the course will be run again and
we are facing the question whether to switch to
numpy/scipy right now or to delay this for one more year.

Reasons to switch
+ numpy is better in several respects compared to Numeric
+ numpy/scipy installs more easily on newer machines
+ students will learn the most recent tools
+ extensive and alive documentation on www.scipy.org

Reasons not to switch
- there is no enthought edition yet (right?)
- there are only packages for a few platforms/distribution
- we need scipy.xplt
  (matplotlib is still no option at this point)

Discussion/Background:

To us the two main show-stoppers are scipy.xplt and the question
about an Enthought Edition  for Windows.
For the Pool of PCs where the tutorial groups are to be held,
it won't be a problem to install numpy/scipy in such
a way that scipy.sandbox.xplt is visible as scipy.xplt
(at least I hope).
However, the students will either have windows (around 80%)
or Linux at home. For windows users we have used
the Enthought Edition (http://code.enthought.com/enthon/)
and linux users were pointed to
available packages for their machines or to install Numeric/scipy
themselves.

Concerning xplt another option might be to
install scipy.sandbox.xplt in such a way
that a `import xplt` would work. If that is possible we could
try to supply `xplt` separately for some of the distributions,
and maybe also for windows (which I don't use, so I have
no idea how difficult that would be).
If something like this was possible, the main question is
whether a new enthon distribution with new numpy/scipy/ipython
and all the other niceties of mayavi/VTK/wxPython/....
will come out in the near future?

I would really love to use the new numpy/scipy - so any ideas
are very welcome!

Best, Arnd

Re: [Numpy-discussion] New release of NumPy coming

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

A Dimecres 15 Febrer 2006 11:18, Pearu Peterson va escriure:
> On Wed, 15 Feb 2006, Francesc Altet wrote:
> > I've been trying to see how to correctly load the unicode tests, but
> > failed miserably. Perhaps Pearu can tell us about the correct way to
> > do that.
>
> I have fixed it in svn. When importing modules from tests/ directory, one
> must surround the corresponding import statements with set_local_path()
> and restore_path() calls.

Ah, ok. Is there any place where this is explained or we have to use
the source to figure out these sort of things?

Thanks,

=2D-=20
>0,0<   Francesc Altet =A0 =A0 http://www.carabos.com/
V   V   C=E1rabos Coop. V. =A0=A0Enjoy Data
 "-"

Re: [Numpy-discussion] New release of NumPy coming

[Pearu P. <pe...@sc...>]

On Wed, 15 Feb 2006, Francesc Altet wrote:

> I've been trying to see how to correctly load the unicode tests, but
> failed miserably. Perhaps Pearu can tell us about the correct way to
> do that.

I have fixed it in svn. When importing modules from tests/ directory, one 
must surround the corresponding import statements with set_local_path() 
and restore_path() calls.

Regards,
Pearu

Re: [Numpy-discussion] avoiding the matrix copy performance hit

[Ed S. <sch...@ft...>]

Travis Oliphant wrote:

> I think I originally tried to make mat *not* return a copy, but this
> actually broke code in SciPy.  So, I left the default as it was as a
> copy on input.  There is an *asmatrix* command that does not return a
> copy...

All SciPy's unit tests actually pass with a default of copy=False in the
matrix constructor.  So SciPy needn't be the blocker here.  I'd like to
cast a vote for not copying by default, in the interests of efficiency
and, as Bill Baxter argued, usefulness.

-- Ed

Re: [Numpy-discussion] New release of NumPy coming

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

A Dimecres 15 Febrer 2006 00:01, Travis Oliphant va escriure:
> I'd like to make a new release of NumPy in the next day or two.   If
> there are any outstanding issues, please let me know.

Well, I've run all the recently added tests for unicode in both UCS4
and UCS2 platforms and all passes flawlessly. Very good :-)

However, there are some problem when trying to load the unicode tests
when running the complete suite:

In [1]: import numpy

In [2]: numpy.test(1)
  Found 11 tests for numpy.core.umath
  Found 8 tests for numpy.lib.arraysetops
  Found 26 tests for numpy.core.ma
  Found 6 tests for numpy.core.records
  Found 14 tests for numpy.core.numeric
  Found 4 tests for numpy.distutils.misc_util
  Found 3 tests for numpy.lib.getlimits
  Found 30 tests for numpy.core.numerictypes
  Found 9 tests for numpy.lib.twodim_base
  Found 1 tests for numpy.core.oldnumeric
  Found 44 tests for numpy.lib.shape_base
  Found 4 tests for numpy.lib.index_tricks
  Found 42 tests for numpy.lib.type_check
  Found 3 tests for numpy.dft.helper
Warning:    !! FAILURE importing tests for <module 'numpy.core.multiarray'=
=20
from '...kages/numpy/core/multiarray.so'>
/usr/lib/python2.3/site-packages/numpy/core/tests/test_multiarray.py:195:=20
ImportError: No module named test_unicode (in ?)
  Found 7 tests for numpy.core.defmatrix
  Found 33 tests for numpy.lib.function_base
  Found 0 tests for __main__
=2E........................................................................=
=2E........................................................................=
=2E........................................................................=
=2E...........................
=2D---------------------------------------------------------------------
Ran 247 tests in 0.951s

OK

I've been trying to see how to correctly load the unicode tests, but
failed miserably. Perhaps Pearu can tell us about the correct way to
do that.

Thanks,

>0,0<   Francesc Altet =A0 =A0 http://www.carabos.com/
V   V   C=E1rabos Coop. V. =A0=A0Enjoy Data
 "-"

Re: [Numpy-discussion] New release of NumPy coming

[Arnd B. <arn...@we...>]

On Tue, 14 Feb 2006, Travis Oliphant wrote:

> I'd like to make a new release of NumPy in the next day or two.   If
> there are any outstanding issues, please let me know.

It seems that the icc stuff has fallen through the cracks,
  http://article.gmane.org/gmane.comp.python.numeric.general/3517/

(it is not relevant to me at this point - it was only a test
after your request for compilations with compilers other than gcc ;-).

Best, Arnd

Re: [Numpy-discussion] New release of NumPy coming

[Pearu P. <pe...@sc...>]

On Tue, 14 Feb 2006, Andrew Straw wrote:

> Travis Oliphant wrote:
>
>> 
>> I'd like to make a new release of NumPy in the next day or two.   If there 
>> are any outstanding issues, please let me know.
>> 
> Here's one:
> http://projects.scipy.org/scipy/numpy/ticket/4

Fixed in svn.

Pearu

Re: [Numpy-discussion] Re: NumPy Glossary Was:Matlab page on scipy wiki

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

Bill Baxter wrote:

> On the point of professionalism, I'd like to change the matlab page's
> title from "NumPy for Matlab Addicts" to simply "NumPy for Matlab
> Users".  It's been bugging me since I put it up there initially... but
> I'm not really sure how to chage the name of a page in the wiki.

I went ahead and renamed the page and created a redirect from the old page.

Re: [Numpy-discussion] New release of NumPy coming

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

Travis Oliphant wrote:

>
> I'd like to make a new release of NumPy in the next day or two.   If 
> there are any outstanding issues, please let me know.
>
Here's one:
http://projects.scipy.org/scipy/numpy/ticket/4

Re: [Numpy-discussion] New release of NumPy coming

[Sasha <nd...@ma...>]

I was deliquent checking in ma code that takes advantage of context in
__array__ .  Will try to do it tomorrow. (Need to add more tests.)

On 2/14/06, Travis Oliphant <oli...@ie...> wrote:
>
> I'd like to make a new release of NumPy in the next day or two.   If
> there are any outstanding issues, please let me know.
>
> -Travis
>
>
>
> -------------------------------------------------------
> This SF.net email is sponsored by: Splunk Inc. Do you grep through log fi=
les
> for problems?  Stop!  Download the new AJAX search engine that makes
> searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
> http://sel.as-us.falkag.net/sel?cmd=3Dlnk&kid=3D103432&bid=3D230486&dat=
=3D121642
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>

[Numpy-discussion] optimizing power() for complex and real cases (was: indexing problem)

[<co...@ph...>]

[changed subject to reflect this better]

Tim Hochberg <tim...@co...> writes:
> David M. Cooke wrote:
>>Tim Hochberg <tim...@co...> writes:
>>>David M. Cooke wrote:
>>>  2, Don't distinguish between scalars and arrays -- that just makes
>>>things harder to explain.
>>Makes the optimizations better, though.
>>
>>
> Ah, Because you can hoist all the checks for what type of optimization
> to do, if any, out of the core loop, right? That's a good point. Still
> I'm not keen on a**b having different performance *and* different
> results depending on whether b is a scalar or matrix. The first thing
> to do is to measure how much overhead doing the optimization element
> by element is going to add. Assuming that it's signifigant that leaves
> us with the familiar dilema: fast, simple or general purpose; pick any
> two.
>
> 1. Do what I've proposed: optimize things at the c_pow level. This is
> general purpose and relatively simple to implement (since we can steal
> most of the code from complexobject.c). It may have a signifigant
> speed penalty versus 2 though:
>
> 2. Do what you've proposed: optimize things at the ufunc level. This
> fast and relatively simple to implement. It's more limited in scope
> and a bit harder to explain than 2.
>
> 3. Do both. This is straightforward, but adds a bunch of extra code
> paths with all the attendant required testing and possibility for
> bugs. So, fast, general purpose, but not simple.

Start with #1, then try #2. The problem with #2 is that you still have
to include #1: if you're doing x**y when y is an array, then you have
do if (y==2) etc. checks in your inner loop anyways. In that case, you
might as well do it in nc_pow. At that point, it may be better to move
the #1 optimization to the level of x.__pow__ (see below).

>>The speed comes from the inlining of how to do the power _outside_ the
>>inner loop. The reason x**2, etc. are slower currently is there is a
>>function call in the inner loop. Your's and mine C library's pow()
>>function mostly likely does something like I have above, for a single
>>case: pow(x, 2.0) is calculated as x*x. However, each time through it
>>has do decide _how_ to do it.
>>
>>
> Part of our difference in perspective comes from the fact that I've
> just been staring at the guts of complex power. In this case you
> always have function calls at present, even for s*s. (At least I'm
> fairly certain that doesn't get inlined although I haven't checked).
> Since much of the work I do is with complex matrices,  it's
> appropriate that I focus on this.

Ah, ok, now things are clicking. Complex power is going to be harder,
because making sure that going from x**2.001 to x**2 doesn't do some
funny complex branch cut stuff (I work in reals all the time :-)

For the real numbers, these type of optimizations *are* a big win, and
don't have the same type of continuity problems. I'll put them into numpy soon.

> Have you measured the effect of a function call on the speed here, or
> is that just an educated guess. If it's an educated guess, it's
> probably worth determining how much of speed hit the function call
> actually causes.  I was going to try to get a handle on this by
> comparing multiplication of Complex numbers (which requires a function
> call plus more math), with multiplication of Floats which does not.
> Perversly, the Complex multiplication came out marginally faster,
> which is hard to explain whichever way you look at it.
>
>>>> timeit.Timer("a*b", "from numpy import arange; a =
> arange(10000)+0j; b = arange(10000)+0j").time
> it(10000)
> 3.2974959107959876
>>>> timeit.Timer("a*b", "from numpy import arange; a = arange(10000);
>>>> b
> = arange(10000)").timeit(100
> 00)
> 3.4541194481425919

You're not multiplying floats in the last one: you're multiplying
integers. You either need to use a = arange(10000.0), or a =
arange(10000.0, dtype=float) (to be more specific).

Your integer numbers are about 3x better than mine, though (difference
in architecture, maybe? I'm on an Athlon64).

>>That's why I limited the optimization to scalar exponents: array
>>exponents would mean it's about as slow as the pow() call, even if the
>>checks were inlined into the loop. It would probably be even slower
>>for the non-optimized case, as you'd check for the special exponents,
>>then call pow() if it fails (which would likely recheck the exponents).
>>
>>
> Again, here I'm thinking of the complex case. In that case at least, I
> don't think that the non-optimized case would take a noticeable speed
> hit. I would put it into pow itself, which already special cases a==0
> and b==0. For float pow it might, but that's already slow, so I doubt
> that it would make much difference.

It does make a bit of difference with float pow: the general case
slows down a bit.

>>Maybe a simple way to add this is to rewrite x.__pow__() as something
>>like the C equivalent of
>>
>>def __pow__(self, p):
>>    if p is not a scalar:
>>       return power(self, p)
>>    elif p == 1:
>>       return p
>>    elif p == 2:
>>       return square(self)
>>    elif p == 3:
>>       return cube(self)
>>    elif p == 4:
>>       return power_4(self)
>>    elif p == 0:
>>       return ones(self.shape, dtype=self.dtype)
>>    elif p == -1:
>>       return 1.0/self
>>    elif p == 0.5:
>>       return sqrt(self)
>>
>>and add ufuncs square, cube, power_4 (etc.).
>
> It sounds like we need to benchmark some stuff and see what we come up
> with. One approach would be for each of us to implement this for one
> time (say float) and see how the approaches compare speed wise. That's
> not entirely fair as my approach will do much better at complex than
> float I believe, but it's certainly easier.

The way the ufuncs are templated, we can split out the complex
routines easily enough.

Here's what I propose:

- add a square() ufunc, where square(x) == x*x (but faster of course)
- I'll fiddle with the floats
- you fiddle with the complex numbers :-)

I've created a new branch in svn, at
http://svn.scipy.org/svn/numpy/branches/power_optimization
to do this fiddling. The changes below I mention are all checked in as
revision 2104 (http://projects.scipy.org/scipy/numpy/changeset/2104).

I've added a square() ufunc to the power_optimization branch because
I'd argue that it's probably *the* most common use of **. I've
implemented it, and it's as fast as a*a for reals, and runs in 2/3 the
time as a*a for complex (which makes sense: squaring a complex numbers
has 3 real multiplications, while multiplying has 4 in the (simple)
scheme [1]).

At least with square(), there's no argument about continuity, as it
only squares :-).

The next step I'd suggest is special-casing x.__pow__, like I suggest
above. We could just test for integer scalar exponents (0, 1, 2), and
just special-case those (returning ones(), x.copy(), square(x)), and
leave all the rest to power().

I've also checked in code to the power_optimization branch that
special cases power(x, <scalar exponent>), or anytime the basic ufunc
gets called with a stride of 0 for the exponent. It doesn't do complex
x, so no problems on your side, but it's a good chunk faster for this
case than what we've got now. One reason I'm also looking at adding
square() is because my optimization of power() makes x**2 run (only)
1.5 slower than x*x (and I can't for the life of me see where that 0.5
is coming from! It should be 1.0 like square()!).

[1] which brings up another point. Would using the 3-multiplication
version for complex multiplication be good? There might be some
effects with cancellation errors due to the extra subtractions...

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|co...@ph...

[Numpy-discussion] New release of NumPy coming

[Travis O. <oli...@ie...>]

I'd like to make a new release of NumPy in the next day or two.   If 
there are any outstanding issues, please let me know.

-Travis

Re: [Numpy-discussion] NumPy Glossary: a request for review

[Christopher B. <Chr...@no...>]

Sasha wrote:
> However, if others think a link to more detailed explanation belongs
> to glossary entries, the natural destination of the link would be a
> page in Travis' book.

Not unless that glossary is part of the book. Links in the Wiki should 
point to the wiki, or maybe to other openly available sources on the web.

A link isn't critical, but a page about broadcasting would still be 
nice, it's a great feature of numpy.

-Chris


-- 
Christopher Barker, Ph.D.
Oceanographer
                                     		
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chr...@no...

Re: [Numpy-discussion] Re: indexing problem

[Tim H. <tim...@co...>]

David M. Cooke wrote:

>Tim Hochberg <tim...@co...> writes:
>
>  
>
>>David M. Cooke wrote:
>>
>>    
>>
>>>Gary Ruben <gr...@bi...> writes:
>>>
>>>
>>>
>>>      
>>>
>>>>Tim Hochberg wrote:
>>>><snip>
>>>>
>>>>
>>>>        
>>>>
>>>>>However, I'm not convinced this is a good idea for numpy. This would
>>>>>introduce a discontinuity in a**b that could cause problems in some
>>>>>cases. If, for instance, one were running an iterative solver of
>>>>>some sort (something I've been known to do), and b was a free
>>>>>variable, it could get stuck at b = 2 since things would go
>>>>>nonmonotonic there.
>>>>>
>>>>>
>>>>>          
>>>>>
>>>>I don't quite understand the problem here. Tim says Python special
>>>>cases integer powers but then talks about the problem when b is a
>>>>floating type. I think special casing x**2 and maybe even x**3 when
>>>>the power is an integer is still a good idea.
>>>>
>>>>
>>>>        
>>>>
>>>Well, what I had done with Numeric did special case x**0, x**1,
>>>x**(-1), x**0.5, x**2, x**3, x**4, and x**5, and only when the
>>>exponent was a scalar (so x**y where y was an array wouldn't be). I
>>>think this is very useful, as I don't want to microoptimize my code to
>>>x*x instead of x**2. The reason for just scalar exponents was so
>>>choosing how to do the power was lifted out of the inner loop. With
>>>that, x**2 was as fast as x*x.
>>>
>>>
>>>      
>>>
>>This is getting harder to object to since, try as I might I can't get
>>a**b to go nonmontonic in the vicinity of b==2. I run out of floating
>>point resolution before the slight shift due to special casing at 2
>>results in nonmonoticity. I suspect that I could manage it with enough
>>work, but it would require some unlikely function of a**b. I'm not
>>sure if I'm really on board with this, but let me float a slightly
>>modified proposal anyway:
>>
>>   1. numpy.power stays as it is now. That way in the rare case that
>>someone runs into trouble they can drop back to power. Alternatively
>>there could be rawpower and power where rawpower has the current
>>behaviour. While the name rawpower sounds cool/cheesy, power is used
>>infrequently enough that I doubt it matters whether it has these
>>special case optimazations.
>>    
>>
>
>+1
>
>  
>
>>  2, Don't distinguish between scalars and arrays -- that just makes
>>things harder to explain.
>>    
>>
>
>Makes the optimizations better, though.
>  
>
Ah, Because you can hoist all the checks for what type of optimization 
to do, if any, out of the core loop, right? That's a good point. Still 
I'm not keen on a**b having different performance *and* different 
results depending on whether b is a scalar or matrix. The first thing to 
do is to measure how much overhead doing the optimization element  by 
element is going to add. Assuming that it's signifigant that leaves us 
with the familiar dilema: fast, simple or general purpose; pick any two.

1. Do what I've proposed: optimize things at the c_pow level. This is 
general purpose and relatively simple to implement (since we can steal 
most of the code from complexobject.c). It may have a signifigant speed 
penalty versus 2 though:

2. Do what you've proposed: optimize things at the ufunc level. This 
fast and relatively simple to implement. It's more limited in scope and 
a bit harder to explain than 2.

3. Do both. This is straightforward, but adds a bunch of extra code 
paths with all the attendant required testing and possibility for bugs. 
So, fast, general purpose, but not simple.


>  
>
>>  3. Python itself special cases all integral powers between -100 and
>>100. Beg/borrow/steal their code. This makes it easier to explain
>>since all smallish integer powers are just automagically faster.
>>
>>  4. Is the performance advantage of special casing a**0.5
>>signifigant? If so use the above trick to special case all half
>>integral  and integral powers between -N and N. Since sqrt probably
>>chews up some time the cutoff. The cutoff probably shifts somewhat if
>>we're optimizing half integral as well as integral powers. Perhaps N
>>would be 32 or 64.
>>
>>The net result of this is that a**b would be computed using a
>>combination of repeated multiplication and sqrt for  real integral and
>>half integral values of b between -N and N. That seems simpler to
>>explain and somewhat more useful as well.
>>
>>It sounds like a fun project although I'm not certain yet that it's a
>>good idea.
>>    
>>
>
>Basically, my Numeric code looked like this:
>
>#define POWER_UFUNC3(prefix, basetype, exptype, outtype)                \
>static void prefix##_power(char **args, int *dimensions,                \
>                           int *steps, void *func) {                    \
>    int i, cis1=steps[0], cis2=steps[1], cos=steps[2], n=dimensions[0]; \
>    int is1=cis1/sizeof(basetype);                                      \
>    int is2=cis2/sizeof(exptype);                                       \
>    int os=cos/sizeof(outtype);                                         \
>    basetype *i1 = (basetype *)(args[0]);                               \
>    exptype *i2=(exptype *)(args[1]);                                   \
>    outtype *op=(outtype *)(args[2]);                                   \
>    if (is2 == 0) {                                                     \
>        exptype exponent = i2[0];                                       \
>        if (POWER_equal(exponent, 0.0)) {                               \
>            for (i = 0; i < n; i++, op += os) {                         \
>                POWER_one((*op))                                        \
>            }                                                           \
>        } else if (POWER_equal(exponent, 1.0)) {                        \
>            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
>                *op = *i1;                                              \
>            }                                                           \
>        } else if (POWER_equal(exponent, 2.0)) {                        \
>            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
>                POWER_square((*op),(*i1))                               \
>            }                                                           \
>        } else if (POWER_equal(exponent, -1.0)) {                       \
>            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
>                POWER_inverse((*op),(*i1))                              \
>            }                                                           \
>        } else if (POWER_equal(exponent, 3.0)) {                        \
>            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
>                POWER_cube((*op),(*i1))                                 \
>            }                                                           \
>        } else if (POWER_equal(exponent, 4.0)) {                        \
>            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
>                POWER_fourth((*op),(*i1))                               \
>            }                                                           \
>        } else if (POWER_equal(exponent, 0.5)) {                        \
>            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
>                POWER_sqrt((*op),(*i1))                                 \
>            }                                                           \
>        } else {                                                        \
>            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
>                POWER_pow((*op), (*i1), (exponent))                     \
>            }                                                           \
>        }                                                               \
>    } else {                                                            \
>        for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {                 \
>            POWER_pow((*op), (*i1), (*i2))                              \
>        }                                                               \
>    }                                                                   \
>}
>#define POWER_UFUNC(prefix, type) POWER_UFUNC3(prefix, type, type, type)
>
>#define FTYPE float
>#define POWER_equal(x,y)   x == y
>#define POWER_one(o)       o = 1.0;
>#define POWER_square(o,x)  o = x*x;
>#define POWER_inverse(o,x) o = 1.0 / x;
>#define POWER_cube(o,x)    FTYPE y=x; o = y*y*y;
>#define POWER_fourth(o,x)  FTYPE y=x, s = y*y; o = s * s;
>#define POWER_sqrt(o,x)    o = sqrt(x);
>#define POWER_pow(o,x,n)   o = pow(x, n);
>POWER_UFUNC(FLOAT, float)
>POWER_UFUNC3(FLOATD, float, double, float)
>
>plus similiar definitions for float, double, complex float, and
>complex double. Using the POWER_square, etc. macros means the complex
>case was easy to add.
>
>The speed comes from the inlining of how to do the power _outside_ the
>inner loop. The reason x**2, etc. are slower currently is there is a
>function call in the inner loop. Your's and mine C library's pow()
>function mostly likely does something like I have above, for a single
>case: pow(x, 2.0) is calculated as x*x. However, each time through it
>has do decide _how_ to do it.
>  
>
Part of our difference in perspective comes from the fact that I've just 
been staring at the guts of complex power. In this case you always have 
function calls at present, even for s*s. (At least I'm fairly certain 
that doesn't get inlined although I haven't checked). Since much of the 
work I do is with complex matrices,  it's appropriate that I focus on this.

Have you measured the effect of a function call on the speed here, or is 
that just an educated guess. If it's an educated guess, it's probably 
worth determining how much of speed hit the function call actually 
causes.  I was going to try to get a handle on this by comparing 
multiplication of Complex numbers (which requires a function call plus 
more math), with multiplication of Floats which does not. Perversly, the 
Complex multiplication came out marginally faster, which is hard to 
explain whichever way you look at it.

 >>> timeit.Timer("a*b", "from numpy import arange; a = 
arange(10000)+0j; b = arange(10000)+0j").time
it(10000)
3.2974959107959876
 >>> timeit.Timer("a*b", "from numpy import arange; a = arange(10000); b 
= arange(10000)").timeit(100
00)
3.4541194481425919

>That's why I limited the optimization to scalar exponents: array
>exponents would mean it's about as slow as the pow() call, even if the
>checks were inlined into the loop. It would probably be even slower
>for the non-optimized case, as you'd check for the special exponents,
>then call pow() if it fails (which would likely recheck the exponents).
>  
>
Again, here I'm thinking of the complex case. In that case at least, I 
don't think that the non-optimized case would take a noticeable speed 
hit. I would put it into pow itself, which already special cases a==0 
and b==0. For float pow it might, but that's already slow, so I doubt 
that it would make much difference.

>Maybe a simple way to add this is to rewrite x.__pow__() as something
>like the C equivalent of
>
>def __pow__(self, p):
>    if p is not a scalar:
>       return power(self, p)
>    elif p == 1:
>       return p
>    elif p == 2:
>       return square(self)
>    elif p == 3:
>       return cube(self)
>    elif p == 4:
>       return power_4(self)
>    elif p == 0:
>       return ones(self.shape, dtype=self.dtype)
>    elif p == -1:
>       return 1.0/self
>    elif p == 0.5:
>       return sqrt(self)
>
>and add ufuncs square, cube, power_4 (etc.).
>  
>

It sounds like we need to benchmark some stuff and see what we come up 
with. One approach would be for each of us to implement this for one 
time (say float) and see how the approaches compare speed wise. That's 
not entirely fair as my approach will do much better at complex than 
float I believe, but it's certainly easier.


regards,

-tim

[Numpy-discussion] transpose function

[Robert C. <cim...@nt...>]

Just now I have stumbled on a non-intuitive thing with the transpose 
function.

The help says:
"""
Help on function transpose in module numpy.core.oldnumeric:
transpose(a, axes=None)
     transpose(a, axes=None) returns array with dimensions permuted
     according to axes.  If axes is None (default) returns array with
     dimensions reversed.
"""

There are many functions in scipy that accept the 'axis' argument which 
is a single integer number. I have overlooked that here it is 'axes' and 
see what happens (I expected 'flipping' the array around the given 
single axis, well...):

import scipy as nm

In [26]:b = nm.zeros( (3,4) )
In [27]:b
Out[27]:
array([[0, 0, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 0]])

In [29]:nm.transpose( b, 0 )
Out[29]:array([0, 0, 0])

In [30]:nm.transpose( b, 1 )
Out[30]:array([0, 0, 0, 0])

So I propose either to replace 'axes' with 'order' or give an example in 
the docstring. It would be also good to raise an exception when the 
length of the 'axes' argument does not match the array rank and/or does 
not contain a permutation (no repetitions) of relevant indices.

What do the gurus think?

r.

Re: [Numpy-discussion] Re: indexing problem

[Arnd B. <arn...@we...>]

On Mon, 13 Feb 2006, David M. Cooke wrote:

> Gary Ruben <gr...@bi...> writes:
>
> > Tim Hochberg wrote:
> > <snip>
> >> However, I'm not convinced this is a good idea for numpy. This would
> >> introduce a discontinuity in a**b that could cause problems in some
> >> cases. If, for instance, one were running an iterative solver of
> >> some sort (something I've been known to do), and b was a free
> >> variable, it could get stuck at b = 2 since things would go
> >> nonmonotonic there.
> >
> > I don't quite understand the problem here. Tim says Python special
> > cases integer powers but then talks about the problem when b is a
> > floating type. I think special casing x**2 and maybe even x**3 when
> > the power is an integer is still a good idea.
>
> Well, what I had done with Numeric did special case x**0, x**1,
> x**(-1), x**0.5, x**2, x**3, x**4, and x**5, and only when the
> exponent was a scalar (so x**y where y was an array wouldn't be). I
> think this is very useful, as I don't want to microoptimize my code to
> x*x instead of x**2. The reason for just scalar exponents was so
> choosing how to do the power was lifted out of the inner loop. With
> that, x**2 was as fast as x*x.

+1 from me for the special casing.
A speed improvement of more than a factor 16 (from David's numbers)
is something relevant!

Best, Arnd

[Numpy-discussion] Python Equivalent of Matlab lpc

[Curtis S. <th...@gm...>]

Hi,

I am trying to pull cepstral coefficients from wav files for a speech 
recognizer and I am wondering if there is a python equivalent of the lpc 
function in matlab? If not, anyone know of any other good ways to 
featurize speech vectors with python included functions?

Thanks,
Curtis

Re: [Numpy-discussion] Re: indexing problem

[<co...@ph...>]

Tim Hochberg <tim...@co...> writes:

> David M. Cooke wrote:
>
>>Gary Ruben <gr...@bi...> writes:
>>
>>
>>
>>>Tim Hochberg wrote:
>>><snip>
>>>
>>>
>>>>However, I'm not convinced this is a good idea for numpy. This would
>>>>introduce a discontinuity in a**b that could cause problems in some
>>>>cases. If, for instance, one were running an iterative solver of
>>>>some sort (something I've been known to do), and b was a free
>>>>variable, it could get stuck at b = 2 since things would go
>>>>nonmonotonic there.
>>>>
>>>>
>>>I don't quite understand the problem here. Tim says Python special
>>>cases integer powers but then talks about the problem when b is a
>>>floating type. I think special casing x**2 and maybe even x**3 when
>>>the power is an integer is still a good idea.
>>>
>>>
>>
>>Well, what I had done with Numeric did special case x**0, x**1,
>>x**(-1), x**0.5, x**2, x**3, x**4, and x**5, and only when the
>>exponent was a scalar (so x**y where y was an array wouldn't be). I
>>think this is very useful, as I don't want to microoptimize my code to
>>x*x instead of x**2. The reason for just scalar exponents was so
>>choosing how to do the power was lifted out of the inner loop. With
>>that, x**2 was as fast as x*x.
>>
>>
> This is getting harder to object to since, try as I might I can't get
> a**b to go nonmontonic in the vicinity of b==2. I run out of floating
> point resolution before the slight shift due to special casing at 2
> results in nonmonoticity. I suspect that I could manage it with enough
> work, but it would require some unlikely function of a**b. I'm not
> sure if I'm really on board with this, but let me float a slightly
> modified proposal anyway:
>
>    1. numpy.power stays as it is now. That way in the rare case that
> someone runs into trouble they can drop back to power. Alternatively
> there could be rawpower and power where rawpower has the current
> behaviour. While the name rawpower sounds cool/cheesy, power is used
> infrequently enough that I doubt it matters whether it has these
> special case optimazations.

+1

>
>   2, Don't distinguish between scalars and arrays -- that just makes
> things harder to explain.

Makes the optimizations better, though.

>   3. Python itself special cases all integral powers between -100 and
> 100. Beg/borrow/steal their code. This makes it easier to explain
> since all smallish integer powers are just automagically faster.
>
>   4. Is the performance advantage of special casing a**0.5
> signifigant? If so use the above trick to special case all half
> integral  and integral powers between -N and N. Since sqrt probably
> chews up some time the cutoff. The cutoff probably shifts somewhat if
> we're optimizing half integral as well as integral powers. Perhaps N
> would be 32 or 64.
>
> The net result of this is that a**b would be computed using a
> combination of repeated multiplication and sqrt for  real integral and
> half integral values of b between -N and N. That seems simpler to
> explain and somewhat more useful as well.
>
> It sounds like a fun project although I'm not certain yet that it's a
> good idea.

Basically, my Numeric code looked like this:

#define POWER_UFUNC3(prefix, basetype, exptype, outtype)                \
static void prefix##_power(char **args, int *dimensions,                \
                           int *steps, void *func) {                    \
    int i, cis1=steps[0], cis2=steps[1], cos=steps[2], n=dimensions[0]; \
    int is1=cis1/sizeof(basetype);                                      \
    int is2=cis2/sizeof(exptype);                                       \
    int os=cos/sizeof(outtype);                                         \
    basetype *i1 = (basetype *)(args[0]);                               \
    exptype *i2=(exptype *)(args[1]);                                   \
    outtype *op=(outtype *)(args[2]);                                   \
    if (is2 == 0) {                                                     \
        exptype exponent = i2[0];                                       \
        if (POWER_equal(exponent, 0.0)) {                               \
            for (i = 0; i < n; i++, op += os) {                         \
                POWER_one((*op))                                        \
            }                                                           \
        } else if (POWER_equal(exponent, 1.0)) {                        \
            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
                *op = *i1;                                              \
            }                                                           \
        } else if (POWER_equal(exponent, 2.0)) {                        \
            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
                POWER_square((*op),(*i1))                               \
            }                                                           \
        } else if (POWER_equal(exponent, -1.0)) {                       \
            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
                POWER_inverse((*op),(*i1))                              \
            }                                                           \
        } else if (POWER_equal(exponent, 3.0)) {                        \
            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
                POWER_cube((*op),(*i1))                                 \
            }                                                           \
        } else if (POWER_equal(exponent, 4.0)) {                        \
            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
                POWER_fourth((*op),(*i1))                               \
            }                                                           \
        } else if (POWER_equal(exponent, 0.5)) {                        \
            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
                POWER_sqrt((*op),(*i1))                                 \
            }                                                           \
        } else {                                                        \
            for (i = 0; i < n; i++, i1 += is1, op += os) {              \
                POWER_pow((*op), (*i1), (exponent))                     \
            }                                                           \
        }                                                               \
    } else {                                                            \
        for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {                 \
            POWER_pow((*op), (*i1), (*i2))                              \
        }                                                               \
    }                                                                   \
}
#define POWER_UFUNC(prefix, type) POWER_UFUNC3(prefix, type, type, type)

#define FTYPE float
#define POWER_equal(x,y)   x == y
#define POWER_one(o)       o = 1.0;
#define POWER_square(o,x)  o = x*x;
#define POWER_inverse(o,x) o = 1.0 / x;
#define POWER_cube(o,x)    FTYPE y=x; o = y*y*y;
#define POWER_fourth(o,x)  FTYPE y=x, s = y*y; o = s * s;
#define POWER_sqrt(o,x)    o = sqrt(x);
#define POWER_pow(o,x,n)   o = pow(x, n);
POWER_UFUNC(FLOAT, float)
POWER_UFUNC3(FLOATD, float, double, float)

plus similiar definitions for float, double, complex float, and
complex double. Using the POWER_square, etc. macros means the complex
case was easy to add.

The speed comes from the inlining of how to do the power _outside_ the
inner loop. The reason x**2, etc. are slower currently is there is a
function call in the inner loop. Your's and mine C library's pow()
function mostly likely does something like I have above, for a single
case: pow(x, 2.0) is calculated as x*x. However, each time through it
has do decide _how_ to do it.

That's why I limited the optimization to scalar exponents: array
exponents would mean it's about as slow as the pow() call, even if the
checks were inlined into the loop. It would probably be even slower
for the non-optimized case, as you'd check for the special exponents,
then call pow() if it fails (which would likely recheck the exponents).

Maybe a simple way to add this is to rewrite x.__pow__() as something
like the C equivalent of

def __pow__(self, p):
    if p is not a scalar:
       return power(self, p)
    elif p == 1:
       return p
    elif p == 2:
       return square(self)
    elif p == 3:
       return cube(self)
    elif p == 4:
       return power_4(self)
    elif p == 0:
       return ones(self.shape, dtype=self.dtype)
    elif p == -1:
       return 1.0/self
    elif p == 0.5:
       return sqrt(self)

and add ufuncs square, cube, power_4 (etc.).

-- 
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|co...@ph...

Re: [Numpy-discussion] Re: indexing problem

[Ryan K. <rya...@gm...>]

I think that would be great, but is there any chance there would be a
problem with the scenario Tim posted earlier:  If a script was running
some sort of optimization on x**y, could the y value ever actually be
returned as an integer and could that throw off the optimization if
round off error caused the float version returned a significantly
different value than the integer version?

Ryan

On 2/13/06, Gary Ruben <gr...@bi...> wrote:
> Hi David,
> So, I think what you had done would be OK provided you removed the
> x**0.5 case to avoid the problem Tim raised and checked that the
> exponent is an integer, not just a scalar.
> Does anyone see a problem with this approach.
> Gary R.
>
> David M. Cooke wrote:
> > Gary Ruben <gr...@bi...> writes:
> >
> >> Tim Hochberg wrote:
> >> <snip>
> >>> However, I'm not convinced this is a good idea for numpy. This would
> >>> introduce a discontinuity in a**b that could cause problems in some
> >>> cases. If, for instance, one were running an iterative solver of
> >>> some sort (something I've been known to do), and b was a free
> >>> variable, it could get stuck at b =3D 2 since things would go
> >>> nonmonotonic there.
> >> I don't quite understand the problem here. Tim says Python special
> >> cases integer powers but then talks about the problem when b is a
> >> floating type. I think special casing x**2 and maybe even x**3 when
> >> the power is an integer is still a good idea.
> >
> > Well, what I had done with Numeric did special case x**0, x**1,
> > x**(-1), x**0.5, x**2, x**3, x**4, and x**5, and only when the
> > exponent was a scalar (so x**y where y was an array wouldn't be). I
> > think this is very useful, as I don't want to microoptimize my code to
> > x*x instead of x**2. The reason for just scalar exponents was so
> > choosing how to do the power was lifted out of the inner loop. With
> > that, x**2 was as fast as x*x.
> >
>
>
> -------------------------------------------------------
> This SF.net email is sponsored by: Splunk Inc. Do you grep through log fi=
les
> for problems?  Stop!  Download the new AJAX search engine that makes
> searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
> http://sel.as-us.falkag.net/sel?cmd=3Dlnk&kid=3D103432&bid=3D230486&dat=
=3D121642
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>

Re: [Numpy-discussion] Re: indexing problem

[Tim H. <tim...@co...>]

David M. Cooke wrote:

>Gary Ruben <gr...@bi...> writes:
>
>  
>
>>Tim Hochberg wrote:
>><snip>
>>    
>>
>>>However, I'm not convinced this is a good idea for numpy. This would
>>>introduce a discontinuity in a**b that could cause problems in some
>>>cases. If, for instance, one were running an iterative solver of
>>>some sort (something I've been known to do), and b was a free
>>>variable, it could get stuck at b = 2 since things would go
>>>nonmonotonic there.
>>>      
>>>
>>I don't quite understand the problem here. Tim says Python special
>>cases integer powers but then talks about the problem when b is a
>>floating type. I think special casing x**2 and maybe even x**3 when
>>the power is an integer is still a good idea.
>>    
>>
>
>Well, what I had done with Numeric did special case x**0, x**1,
>x**(-1), x**0.5, x**2, x**3, x**4, and x**5, and only when the
>exponent was a scalar (so x**y where y was an array wouldn't be). I
>think this is very useful, as I don't want to microoptimize my code to
>x*x instead of x**2. The reason for just scalar exponents was so
>choosing how to do the power was lifted out of the inner loop. With
>that, x**2 was as fast as x*x.
>  
>
This is getting harder to object to since, try as I might I can't get 
a**b to go nonmontonic in the vicinity of b==2. I run out of floating 
point resolution before the slight shift due to special casing at 2 
results in nonmonoticity. I suspect that I could manage it with enough 
work, but it would require some unlikely function of a**b. I'm not sure 
if I'm really on board with this, but let me float a slightly modified 
proposal anyway:

    1. numpy.power stays as it is now. That way in the rare case that 
someone runs into trouble they can drop back to power. Alternatively 
there could be rawpower and power where rawpower has the current 
behaviour. While the name rawpower sounds cool/cheesy, power is used 
infrequently enough that I doubt it matters whether it has these special 
case optimazations.

   2, Don't distinguish between scalars and arrays -- that just makes 
things harder to explain.

   3. Python itself special cases all integral powers between -100 and 
100. Beg/borrow/steal their code. This makes it easier to explain since 
all smallish integer powers are just automagically faster.

   4. Is the performance advantage of special casing a**0.5 signifigant? 
If so use the above trick to special case all half integral  and 
integral powers between -N and N. Since sqrt probably chews up some time 
the cutoff. The cutoff probably shifts somewhat if we're optimizing half 
integral as well as integral powers. Perhaps N would be 32 or 64.

The net result of this is that a**b would be computed using a 
combination of repeated multiplication and sqrt for  real integral and 
half integral values of b between -N and N. That seems simpler to 
explain and somewhat more useful as well.

It sounds like a fun project although I'm not certain yet that it's a 
good idea.

-tim

Re: [Numpy-discussion] Re: indexing problem

[Gary R. <gr...@bi...>]

Hi David,
So, I think what you had done would be OK provided you removed the 
x**0.5 case to avoid the problem Tim raised and checked that the 
exponent is an integer, not just a scalar.
Does anyone see a problem with this approach.
Gary R.

David M. Cooke wrote:
> Gary Ruben <gr...@bi...> writes:
> 
>> Tim Hochberg wrote:
>> <snip>
>>> However, I'm not convinced this is a good idea for numpy. This would
>>> introduce a discontinuity in a**b that could cause problems in some
>>> cases. If, for instance, one were running an iterative solver of
>>> some sort (something I've been known to do), and b was a free
>>> variable, it could get stuck at b = 2 since things would go
>>> nonmonotonic there.
>> I don't quite understand the problem here. Tim says Python special
>> cases integer powers but then talks about the problem when b is a
>> floating type. I think special casing x**2 and maybe even x**3 when
>> the power is an integer is still a good idea.
> 
> Well, what I had done with Numeric did special case x**0, x**1,
> x**(-1), x**0.5, x**2, x**3, x**4, and x**5, and only when the
> exponent was a scalar (so x**y where y was an array wouldn't be). I
> think this is very useful, as I don't want to microoptimize my code to
> x*x instead of x**2. The reason for just scalar exponents was so
> choosing how to do the power was lifted out of the inner loop. With
> that, x**2 was as fast as x*x.
> 

Re: [Numpy-discussion] String array equality test does not broadcast

[Travis O. <oli...@ie...>]

Tim Hochberg wrote:

> Travis Oliphant wrote:
>
>> Russel Howe wrote:
>>
>>> I am converting some numarray code to numpy and I noticed this 
>>> behavior:
>>>
>>> >>> from numpy import *
>>> >>> sta=array(['abc', 'def', 'ghi'])
>>> >>> stb=array(['abc', 'jkl', 'ghi'])
>>> >>> sta==stb
>>> False
>>>
>>> I expected the same as this:
>>> >>> a1=array([1,2,3])
>>> >>> a2=array([1,4,3])
>>> >>> a1==a2
>>> array([True, False, True], dtype=bool)
>>>
>>> I am trying to figure out how to fix this now...
>>
>>
>>
>>
>> Equality testing on string arrays does not work (equality testing 
>> uses ufuncs internally which are not supported generally for flexible 
>> arrays).   You must use chararray's.
>
>
> Should string arrays then perhaps raise an exception here to keep 
> people out of trouble?


Probably.  The equal (not_equal) rich comparison code has some left-over 
stuff from Numeric which is implemented so that if the ufunc equal 
(not_equal) failed False (True) was returned.

I did not special-case the string arrays in this code.

-Travis

>
> -tim
>
>
>>
>> Thus,
>>
>> sta.view(chararray) == stb.view(chararray)
>>
>> Or create chararray's from the beginning:
>>
>> sta = char.array(['abc','def','ghi'])
>> stb = char.array(['abc','jkl','ghi'])
>>
>> Char arrays are a special subclass of the ndarray that give arrays 
>> all the methods of strings (and unicode) elements and allow (rich) 
>> comparison operations.
>>
>> -Travis
>>
>>
>>
>>
>>
>>
>>
>>
>> -------------------------------------------------------
>> This SF.net email is sponsored by: Splunk Inc. Do you grep through 
>> log files
>> for problems?  Stop!  Download the new AJAX search engine that makes
>> searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
>> http://sel.as-us.falkag.net/sel?cmd=lnk&kid=103432&bid=230486&dat=121642
>> _______________________________________________
>> Numpy-discussion mailing list
>> Num...@li...
>> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>>
>>
>
>
>
>
> -------------------------------------------------------
> This SF.net email is sponsored by: Splunk Inc. Do you grep through log 
> files
> for problems?  Stop!  Download the new AJAX search engine that makes
> searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
> http://sel.as-us.falkag.net/sel?cmd=lnk&kid=103432&bid=230486&dat=121642
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion