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

Re: [Numpy-discussion] Re: selecting random array element

[Alan G I. <ai...@am...>]

On Mon, 20 Feb 2006, Robert Kern apparently wrote: 
> You may want to consider numpy.random.permutation() 

Yes indeed.
Thanks,
Alan

Re: [Numpy-discussion] Re: Some missing linalg things (wanted: LU decomposition)

[Bill B. <wb...@gm...>]

On 2/20/06, Bill Baxter <wb...@gm...> wrote:

> > Should have mentioned -- I was using numpy 0.9.4 / scipy 0.4.4.
> > Looks like it works in numpy 0.9.5 / scipy 0.4.6
> >
> > But matplotlib, which I also need, hasn't been updated for numpy 0.9.5y=
et.
> > :-(
> >
>
>
On 2/21/06, Jonathan Taylor wrote:
>
> For matplotlib, I just use tolist() like
>
> a =3D array([1,3,2,3])
>
> ...
>
> pylab.plot(a.tolist())
>
> Maybe that will work for you until you can fix your problem.
> J.



Excellent idea!  That does the trick for now (if I take the numerix: numpy
line out of my .matplotlibrc to stop it from crashing on import).

--bb

[Numpy-discussion] Re: selecting random array element

[Robert K. <rob...@gm...>]

Alan G Isaac wrote:
> On Mon, 20 Feb 2006, Robert Kern apparently wrote: 
> 
>>        length = numpy.multiply.reduce(x.shape) 
> 
> Can this be different from x.size?

No, it's just that old habits die hard. I knew there was a clean way to get that
information, I just didn't remember it or bother myself to look for it.

-- 
Robert Kern
rob...@gm...

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

Re: [Numpy-discussion] Re: selecting random array element

[Alan G I. <ai...@am...>]

On Mon, 20 Feb 2006, Robert Kern apparently wrote:=20
>         length =3D numpy.multiply.reduce(x.shape)=20


Can this be different from x.size?

Thanks,
Alan Isaac

[Numpy-discussion] Re: selecting random array element

[Robert K. <rob...@gm...>]

Chris Fonnesbeck wrote:
> What is the best way to select a random element from a numpy array? I
> know I could index by a random integer, but was wondering if there was
> a built-in method or function.

Generating a random index is what I do.

I think there's certainly room for a RandomState.choice() method. I think
something like this covers most of the use cases:

import numpy
from numpy import random
def choice(x, axis=None):
    """Select an element or subarray uniformly randomly.
    If axis is None, then a single element is chosen from the entire array.
    Otherwise, a subarray is chosen from the given axis.
    """
    x = numpy.asarray(x)
    if axis is None:
        length = numpy.multiply.reduce(x.shape)
        n = random.randint(length)
        return x.flat[n]
    else:
        n = random.randint(x.shape[axis])
        # I'm sure there's a better way of doing this
        idx = map(slice, x.shape)
        idx[axis] = n
        return x[tuple(idx)]

-- 
Robert Kern
rob...@gm...

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

[Numpy-discussion] Re: selecting random array element

[Robert K. <rob...@gm...>]

Alan G Isaac wrote:
> At http://www.american.edu/econ/pytrix/pytrix.py find
> def permute(x):
> 	'''Return a permutation of a sequence or array.
> 
> 	:note: Also consider numpy.random.shuffle
> 	       (to permute *inplace* 1-d arrays)
> 	'''
>         x = numpy.asarray(x)
> 	xshape = x.shape
> 	pidx = numpy.random.random(x.size).argsort()
> 	return x.flat[pidx].reshape(xshape)

You may want to consider numpy.random.permutation()

In [22]: numpy.random.permutation?
Type:           builtin_function_or_method
Base Class:     <type 'builtin_function_or_method'>
String Form:    <built-in method permutation of mtrand.RandomState object at
0x3f0f0>
Namespace:      Interactive
Docstring:
    Given an integer, return a shuffled sequence of integers >= 0 and
    < x; given a sequence, return a shuffled array copy.

    permutation(x)

-- 
Robert Kern
rob...@gm...

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

Re: [Numpy-discussion] [nd_image] histograms of RGB images

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

See this example. You don't really need pylab/matplotlib -- you could
just use numpy.histogram.

import pylab
import matplotlib.numerix as nx
import Image

im = Image.open('data/lena.jpg')
imbuf = im.tostring('raw','RGB',0,-1)
imnx = nx.fromstring(imbuf,nx.UInt8)
imnx.shape = im.size[1], im.size[0], 3

bins = nx.arange(0,256)
pylab.hist( nx.ravel(imnx[:,:,0]), bins=bins, facecolor='r', edgecolor='r' )
pylab.hist( nx.ravel(imnx[:,:,1]), bins=bins, facecolor='g', edgecolor='g' )
pylab.hist( nx.ravel(imnx[:,:,2]), bins=bins, facecolor='b', edgecolor='b' )
pylab.show()

Curzio Basso wrote:

>Hello everybody!
>
>I was wondering if someone already had the problem of computing
>histograms of RGB images for all channels simultaneously (that is
>getting a rank-3 array) rather than on the three channels separately.
>Just looking for a way to avoid writing the C function :-)
>
>cheers
>curzio
>
>
>-------------------------------------------------------
>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=k&kid3432&bid#0486&dat1642
>_______________________________________________
>Numpy-discussion mailing list
>Num...@li...
>https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>  
>

Re: [Numpy-discussion] selecting random array element

[Alan G I. <ai...@am...>]

At http://www.american.edu/econ/pytrix/pytrix.py find
def permute(x):
=09'''Return a permutation of a sequence or array.

=09:note: Also consider numpy.random.shuffle
=09       (to permute *inplace* 1-d arrays)
=09'''
        x =3D numpy.asarray(x)
=09xshape =3D x.shape
=09pidx =3D numpy.random.random(x.size).argsort()
=09return x.flat[pidx].reshape(xshape)

Note the note.  ;-)

Cheers,
Alan Isaac

[Numpy-discussion] selecting random array element

[Chris F. <fon...@gm...>]

What is the best way to select a random element from a numpy array? I
know I could index by a random integer, but was wondering if there was
a built-in method or function.

Thanks,
C.

--
Chris Fonnesbeck + Atlanta, GA + http://trichech.us

[Numpy-discussion] [nd_image] histograms of RGB images

[Curzio B. <cur...@gm...>]

Hello everybody!

I was wondering if someone already had the problem of computing
histograms of RGB images for all channels simultaneously (that is
getting a rank-3 array) rather than on the three channels separately.
Just looking for a way to avoid writing the C function :-)

cheers
curzio

Re: [Numpy-discussion] Problems compiling on Cygwin

[Jose L. G. D. <jos...@gm...>]

Hi Travis,
> I looked in to how people at cygwin ports got the IEEE math stuff done.  
> They borrowed it from BSD basically.  So, I've taken their patch and 
> placed it in the main tree.
> 
> Jose, could you check out the latest SVN version of numpy and try to 
> build and install it on cygwin to see if I made the right changes?


I can´t access remote SVNs/CVSs servers from work, but I will download it at
home and try it tonight.

Many thanks!
Jose

-- 
Lust, ein paar Euro nebenbei zu verdienen? Ohne Kosten, ohne Risiko!
Satte Provisionen für GMX Partner: http://www.gmx.net/de/go/partner

Re: [Numpy-discussion] Re: Some missing linalg things (wanted: LU decomposition)

[Bill B. <wb...@gm...>]

Should have mentioned -- I was using numpy 0.9.4 / scipy 0.4.4.
Looks like it works in numpy 0.9.5 / scipy 0.4.6

But matplotlib, which I also need, hasn't been updated for numpy 0.9.5 yet.
:-(

It's also still pretty weird to me that you have to do "from
scipy.linalgimport lu" specifically.  And then after doing that one
import, then all the
other scipy.linalg.* functions magically spring into existence too.  Is tha=
t
sort of hing expected behavior from Python imports?

>>> import numpy as N
>>> import scipy as S
>>> S.linalg.lu
Traceback (most recent call last):
  File "<input>", line 1, in ?
AttributeError: 'module' object has no attribute 'lu'
>>> from scipy.linalg import lu
>>> S.linalg.lu(N.rand(2,2))
(array([[ 0.,  1.],
       [ 1.,  0.]]), array([[ 1.        ,  0.        ],
       [ 0.18553085,  1.        ]]), array([[ 0.71732168,  0.48540043],
       [ 0.        ,  0.61379118]]))
>>> (N.__version__, S.__version__)
('0.9.5', '0.4.6')

--bb

On 2/20/06, Nils Wagner <nw...@me...> wrote:
>
> Bill Baxter wrote:
> > Ack.  I may be able to get references to lu, lu_factor, et al, but
> > they don't actually work with numpy arrays:
> >
> > from scipy.linalg import lu,lu_factor,lu_solve
> > import scipy as S
> > A =3D S.rand(2,2)
> > lu(A)
> > Traceback (most recent call last):
> >   File "<input>", line 1, in ?
> >   File "C:\Python24\Lib\site-packages\scipy\linalg\decomp.py", line
> > 249, in lu
> >     flu, =3D get_flinalg_funcs(('lu',),(a1,))
> >   File "C:\Python24\Lib\site-packages\scipy\linalg\flinalg.py", line
> > 30, in get_flinalg_funcs
> >     t =3D arrays[i].dtypechar
> > AttributeError: 'numpy.ndarray' object has no attribute 'dtypechar'
> >
> >
> > Ok, so, once again, does anyone have an lu_factor / lu_solve
> > implementation in python that I could borrow?
> >
> > Apologies for the monologue.
> >
> > --bb
> >
> >
> > On 2/20/06, *Bill Baxter* <wb...@gm...
> > <mailto:wb...@gm...>> wrote:
> >
> >     Upon further inspection I find that if I call 'from scipy import
> >     *' then linalg.lu <http://linalg.lu> etc are defined.
> >     But if I do anything else to import scipy like 'import scipy' or
> >     'import scipy as S' or 'from scipy import linalg', then lu, cg etc
> >     are not defined.
> >
> >     Why is that?
> >
> >     I can get at them without importing * by doing 'from scipy.linalg
> >     import lu', but that's kind of odd to have to do that.
> >
> >     --bb
> >
> >
> >     On 2/20/06, * Bill Baxter* <wb...@gm...
> >     <mailto:wb...@gm...>> wrote:
> >
> >         This url http://www.rexx.com/~dkuhlman/scipy_course_01.html
> >         <http://www.rexx.com/%7Edkuhlman/scipy_course_01.html> seems
> >         to keep turning up in my searches for numpy and scipy things,
> >         but many of the linalg operations it lists don't seem to exist
> >         in recent versions of numpy (or scipy).
> >
> >         Some of them are:
> >
> >         * norm
> >         * factorizations: lu, lu_factor, lu_solve, qr
> >         * iterative solvers: cg, cgs, gmres etc.
> >
> >         Did these things used to exist in Numeric but they haven't
> >         been ported over?  Will they be re-introduced sometime?
> >
> >         In the short term, the one I'm after right now is LU decompose
> >         and solve functionality.  Anyone have a numpy implementation?
> >
> >         --Bill Baxter
> >
> No problem here.
>
> >>> from scipy.linalg import lu,lu_factor,lu_solve
> >>> import scipy as S
> >>> A =3D S.rand(2,2)
> >>> lu(A)
> (array([[ 0.,  1.],
>        [ 1.,  0.]]), array([[ 1.        ,  0.        ],
>        [ 0.81367315,  1.        ]]), array([[ 0.49886054,  0.57065709],
>        [ 0.        , -0.30862809]]))
> >>> S.__version__
> '0.4.7.1614'
>
>
> Nils
>
>

[Numpy-discussion] Re: Some missing linalg things (wanted: LU decomposition)

[Bill B. <wb...@gm...>]

Ack.  I may be able to get references to lu, lu_factor, et al, but they
don't actually work with numpy arrays:

from scipy.linalg import lu,lu_factor,lu_solve
import scipy as S
A =3D S.rand(2,2)
lu(A)
Traceback (most recent call last):
  File "<input>", line 1, in ?
  File "C:\Python24\Lib\site-packages\scipy\linalg\decomp.py", line 249, in
lu
    flu, =3D get_flinalg_funcs(('lu',),(a1,))
  File "C:\Python24\Lib\site-packages\scipy\linalg\flinalg.py", line 30, in
get_flinalg_funcs
    t =3D arrays[i].dtypechar
AttributeError: 'numpy.ndarray' object has no attribute 'dtypechar'


Ok, so, once again, does anyone have an lu_factor / lu_solve implementation
in python that I could borrow?

Apologies for the monologue.

--bb


On 2/20/06, Bill Baxter <wb...@gm...> wrote:
>
> Upon further inspection I find that if I call 'from scipy import *' then
> linalg.lu etc are defined.
> But if I do anything else to import scipy like 'import scipy' or 'import
> scipy as S' or 'from scipy import linalg', then lu, cg etc are not define=
d.
>
>
> Why is that?
>
> I can get at them without importing * by doing 'from scipy.linalg import
> lu', but that's kind of odd to have to do that.
>
> --bb
>
> On 2/20/06, Bill Baxter <wb...@gm...> wrote:
> >
> > This url http://www.rexx.com/~dkuhlman/scipy_course_01.html<http://www.=
rexx.com/%7Edkuhlman/scipy_course_01.html>seems to keep turning up in my se=
arches for numpy and scipy things,
> > but many of the linalg operations it lists don't seem to exist in recen=
t
> > versions of numpy (or scipy).
> >
> > Some of them are:
> >
> > * norm
> > * factorizations: lu, lu_factor, lu_solve, qr
> > * iterative solvers: cg, cgs, gmres etc.
> >
> > Did these things used to exist in Numeric but they haven't been ported
> > over?  Will they be re-introduced sometime?
> >
> > In the short term, the one I'm after right now is LU decompose and solv=
e
> > functionality.  Anyone have a numpy implementation?
> >
> > --Bill Baxter
> >
>

[Numpy-discussion] Re: Some missing linalg things (wanted: LU decomposition)

[Bill B. <wb...@gm...>]

Upon further inspection I find that if I call 'from scipy import *' then
linalg.lu etc are defined.
But if I do anything else to import scipy like 'import scipy' or 'import
scipy as S' or 'from scipy import linalg', then lu, cg etc are not defined.


Why is that?

I can get at them without importing * by doing 'from scipy.linalg import
lu', but that's kind of odd to have to do that.

--bb

On 2/20/06, Bill Baxter <wb...@gm...> wrote:
>
> This url http://www.rexx.com/~dkuhlman/scipy_course_01.html<http://www.re=
xx.com/%7Edkuhlman/scipy_course_01.html>seems to keep turning up in my sear=
ches for numpy and scipy things,
> but many of the linalg operations it lists don't seem to exist in recent
> versions of numpy (or scipy).
>
> Some of them are:
>
> * norm
> * factorizations: lu, lu_factor, lu_solve, qr
> * iterative solvers: cg, cgs, gmres etc.
>
> Did these things used to exist in Numeric but they haven't been ported
> over?  Will they be re-introduced sometime?
>
> In the short term, the one I'm after right now is LU decompose and solve
> functionality.  Anyone have a numpy implementation?
>
> --Bill Baxter
>



--
William V. Baxter III
OLM Digital
Kono Dens Building Rm 302
1-8-8 Wakabayashi Setagaya-ku
Tokyo, Japan  154-0023
+81 (3) 3422-3380

[Numpy-discussion] Some missing linalg things (wanted: LU decomposition)

[Bill B. <wb...@gm...>]

This url http://www.rexx.com/~dkuhlman/scipy_course_01.html seems to keep
turning up in my searches for numpy and scipy things,
but many of the linalg operations it lists don't seem to exist in recent
versions of numpy (or scipy).

Some of them are:

* norm
* factorizations: lu, lu_factor, lu_solve, qr
* iterative solvers: cg, cgs, gmres etc.

Did these things used to exist in Numeric but they haven't been ported
over?  Will they be re-introduced sometime?

In the short term, the one I'm after right now is LU decompose and solve
functionality.  Anyone have a numpy implementation?

--Bill Baxter

[Numpy-discussion] What is the status of the multidimensional arrays PEP?

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

What is the status of the multidimensional arrays PEP?
<http://numeric.scipy.org/PEP.txt>

It seems to me that there is one part of the PEP that can be easily
separated into a rather uncontroversial PEP.  This is the part that
defines array protocol.  Python already has a (1-dimensional) array
object in the standard library. Python array already supports buffer
protocol and it looks like implementing full array protocol is
straightforward.

I believe that having an object that supports array protocol even
without multiple dimensions will be immediately useful.

Re: [Numpy-discussion] optimizing power() for complex and real cases

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

David M. Cooke wrote:

>On Sat, Feb 18, 2006 at 06:17:47PM -0700, Tim Hochberg wrote:
>  
>
>>OK, I now have a faily clean implementation in C 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)
>>
>>
>>First a couple of technical questions, then on to the philosophical portion 
>>of this note.
>>
>>1. Is there a nice fast way to get a matrix filled with ones from C. I've 
>>been tempted to write a ufunc 'ones_like', but I'm afraid that might be 
>>considered inappropriate.
>>
>>2. Are people aware that array_power is sometimes passed non arrays as its 
>>first argument? Despite having the signature:
>>
>>array_power(PyArrayObject *a1, PyObject *o2) 
>>
>>This caused me almost no end of headaches, not to mention crashes during 
>>numpy.test().
>>    
>>
>
>Yes; because it's the implementation of __pow__, the second argument can
>be anything.
>  
>
No, you misunderstand.. What I was talking about was that the *first* 
argument can also be something that's not a PyArrayObject, despite the 
functions signature.

>  
>
>>I'll check this into the power_optimization branch RSN, hopefully with a 
>>fix for the zero power case. Possibly also after extending it to inplace 
>>power as well.
>>
>>
>>OK, now on to more important stuff. As I've been playing with this my 
>>opinion has gone in circles a couple of times. I now think the issue of 
>>optimizing integer powers of complex numbers and integer powers of floats 
>>are almost completely different. Because complex powers are quite slow and 
>>relatively inaccurate, it is appropriate to optimize them for integer 
>>powers at the level of nc_pow. This should be just a matter of liberal 
>>borrowing from complexobject.c, but I haven't tried it yet.
>>    
>>
>
>Ok.
>
>  
>
>>On the other hand, real powers are fast enough that doing anything at the 
>>single element level is unlikely to help. So in that case we're left with 
>>either optimizing the cases where the dimension is zero as David has done, 
>>or optimizing at the __pow__ (AKA array_power) level as I've done now based 
>>on David's original suggestion. This second approach is faster because it 
>>avoids the mysterious python scalar -> zero-D array conversion overhead. 
>>However, it suffers if we want to optimize lots of different powers since 
>>one needs a ufunc for each one. So the question becomes, which powers 
>>should we optimize?
>>    
>>
>
>Hmm, ufuncs are passed a void* argument for passing info to them. Now,
>what that argument is defined when the ufunc is created, but maybe
>there's a way to piggy-back on it.
>  
>
Yeah, I really felt like I was fighting the ufuncs when I was playing 
with this. On the one hand, you really want to use the ufunc machinery. 
On the other hand that forces you into using the same types for both 
arguments. That really wouldn't be a problem, since we could just define 
an integer_power that took doubles, but did integer powers, except for 
the conversion overhead of Python_Integers into arrays. It looks like 
you started down this road and I played with this as well.   I can think 
a of at least one (horrible) way around the matrix overhead, but the 
real fix would be to dig into PyArray_EnsureArray and see why it's slow 
for Python_Ints. It is much faster for numarray scalars.

Another approach is to actually compute (x*x)*(x*x) for pow(x,4) at the 
level of array_power. I think I could make this work. It would probably 
work well for medium size arrays, but might well make things worse for 
large arrays that are limited by memory bandwidth since it would need to 
move the array from memory into the cache multiple times.

>>My latest thinking on this is that we should optimize only those cases 
>>where the optimized result is no less accurate than that produced by pow. 
>>I'm going to assume that all C operations are equivalently accurate, so 
>>pow(x,2) has roughly the same amount of error as x*x. (Something on the 
>>order of 0.5 ULP I'd guess). In that case:
>>  pow(x, -1) -> 1 / x
>>  pow(x, 0) -> 1
>>  pow(x, 0.5) -> sqrt(x)
>>  pow(x, 1) -> x
>>  pow(x, 2) -> x*x
>>can all be implemented in terms of multiply or divide with the same 
>>accuracy as the original power methods. Once we get beyond these, the error 
>>will go up progressively.
>>
>>The minimal set described above seems like it should be relatively 
>>uncontroversial and it's what I favor. Once we get beyond this basic set, 
>>we would need to reach some sort of consensus on how much additional error 
>>we are willing to tolerate for optimizing these extra cases. You'll notice 
>>that I've changed my mind, yet again, over whether to optimize A**0.5. 
>>Since the set of additional ufuncs needed in this case is relatively small, 
>>just square and inverse (==1/x), this minimal set works well if optimizing 
>>in pow as I've done.
>>    
>>

Just to add a little more confusion to the mix. I did a little testing 
to see how close pow(x,n) and x*x*... actually are. They are slightly 
less close for small values of N and slightly closer for large values of 
N than I would have expected. The upshot of this is that integer powers 
between  -2 and +4 all seem to vary by the same amount when computed 
using pow(x,n) versus multiplies. I'm including the test code at the 
end. Assuming that this result is not a fluke that expands the 
noncontroversial set by at least 3 more values. That's starting to 
strain the ufunc aproach, so perhaps optimizing in @TYP@_power is the 
way to go after all. Or,  more likely, adding @TYP@_int_power or maybe 
@TYP@_fast_power (so as to be able to include some half integer powers) 
and dispatching appropriately from array_power.

The problem here, of course, is the overhead that PyArray_EnsureArray 
runs into. I'm not sure if the ufuncs actually call that, but I was 
using that to convert things to arrays at one point and I saw the 
slowdown, so I suspect that the slowdown is in something 
PyArray_EnsureArray calls if not in that routine itself. I'm afraid to 
dig into that stuff though.. On the other hand, it would probably speed 
up all kinds of stuff if that was sped up.

>
>Ok. I'm still not happy with the speed of pow(), though. I'll have to
>sit and look at. We may be able to optimize integer powers better.
>And there's another area: integer powers of integers. Right now that
>uses pow(), whereas we might be able to do better. 
>
So that's why that's so slow. I assumed it was doing some sort of 
sucessive multiplication. For this, the code that complexobject uses for 
integer powers might be helpful.

>I'm looking into
>that. A framework for that could be helpful for the complex powers too.
>
>Too bad we couldn't make a function generator :-) [Well, we could using
>weave...]\
>  
>
Yaigh!

-tim




def check(n=100000):
    import math
    sqrt = math.sqrt
    failures = {}
    for x in [math.pi, math.e, 1.1]+[1.0+1.0/y for y in range(1,1+n)]:
        for e, expr in [
                        (-5, "1/((x*x)*(x*x)*x)"), 
(-4,"1/((x*x)*(x*x))"), (3, "1/((x*x)*x)"), (1, "x"),
                        (-2, "1/(x*x)"), (-1,"1/x"), (0, "1"), (1, "x"),
                        (2, "x*x"), (3, "x*x*x"), (4, "(x*x)*(x*x)"),
                        (4, "x*x*x*x"),
                        (5, '(x*x)*(x*x)*x'),
                        (6, '(x*x)*(x*x)*(x*x)'),
                        (7, '(x*x)*(x*x)*(x*x)*x'),
                        (8, '((x*x)*(x*x))*((x*x)*(x*x))'),
                        (-1.5, "1/(sqrt(x)*x)"), (-0.5, "1/sqrt(x)"),
                        (0.5, "sqrt(x)"), (1.5, "x*sqrt(x)")]:
            delta =  abs(pow(x,e) - eval(expr, locals())) / pow(x,e)
            if delta:
                key = (e, expr)
                if key not in failures:
                    failures[key] = [(delta, x)]
                failures[key].append((delta, x))
    for key in sorted(failures.keys()):
        e, expr = key
        fstring = ', '.join(str(x) for x in 
list(reversed(sorted(failures[key])))[:1])
        if len(failures[key]) > 1:
            fstring += ', ...'
        print "Failures for x**%s (%s): %s" % (e, expr, fstring)

[Numpy-discussion] Multiple inheritance from ndarray

[Robert L. <rh...@as...>]

I have a swig extension that defines a class that inherits from
both a personal C-coded image struct (actImage), and also from
Numeric's UserArray.  This works very nicely, but I thought that
it was about time to upgrade to numpy.

The code looks like:

from UserArray import *

class Image(UserArray, actImage):
     def __init__(self, *args):
         actImage.__init__(self, *args)
         UserArray.__init__(self, self.getArray(), 'd', copy=False,  
savespace=False)

I can't figure out how to convert this to use ndarray, as ndarray  
doesn't
seem to have an __init__ method, merely a __new__.

So what's the approved numpy way to handle multiple inheritance?   
I've a nasty
idea that this is a python question that I should know the answer to,  
but I'm
afraid that I don't...

				R

[Numpy-discussion] numpy 0.9.4 on cygwin

[Tariq M. <mah...@gm...>]

Hi,
Has anyone been successful at installing numpy 0.9.4 on cygwin?

Details:
os.name: posix
uname -a: CYGWIN_NT-5.1
sys.platform: cygwin
sys.version: 2.4.1
numpy.version: 0.9.4

Steps taken:
1. unpacked numpy-0.9.4.tar.gz
2. changed to numpy-0.9.4 directory
3. python setup.py install

Major error messages while compiling C source:

1. undefined reference to `_feclearexcept' in umathmodule
2. Command "gcc -shared -Wl,--enable-auto-image-base
build/temp.cygwin-1.5.18-i686-2.4/build/src/numpy/core/src/umathmodule.o
-L/usr/lib/python2.4/config -lpython2.4 -o
build/lib.cygwin-1.5.18-i686-2.4/numpy/core/umath.dll" failed with
exit status 1

Other information that might be useful:
1. successfully installed both Numeric 24.2 and numarray 1.5.0

Any help would be appreciated.

Tariq

[Numpy-discussion] number ranges

[Colin J. W. <cj...@sy...>]

<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
  <meta content="text/html;charset=ISO-8859-1" http-equiv="Content-Type">
</head>
<body bgcolor="#ffffff" text="#000000">
Alan G Isaac wrote:
<blockquote
 cite="mid...@am..."
 type="cite">
  <pre wrap="">On Sun, 19 Feb 2006, "Colin J. Williams" apparently wrote: 
  </pre>
  <blockquote type="cite">
    <pre wrap="">The conditional expression (PEP 308) was rejected but is currently being 
implemented in Python version 2.5.  It will have the syntax A if C else B. 
    </pre>
  </blockquote>
  <pre wrap=""><!---->
It's coming
<a class="moz-txt-link-freetext" href="http://www.python.org/peps/pep-0308.html">http://www.python.org/peps/pep-0308.html</a>
But in 2.5?
<a class="moz-txt-link-freetext" href="http://www.python.org/dev/doc/devel/whatsnew/whatsnew25.html">http://www.python.org/dev/doc/devel/whatsnew/whatsnew25.html</a>

Thank you,
Alan Isaac
  </pre>
</blockquote>
It's true that Andrew Kuchling doesn't list it in his What's New.<br>
<br>
I was quoting Guido van Rossum.&nbsp; Please see below.<br>
<br>
My purpose in mentioning PEP 308 was to illustrate that there can be
second thoughts and to support the idea that there should be second
thoughts with respect to PEP 204.<br>
<pre wrap="">see <a class="moz-txt-link-freetext"
 href="http://www.python.org/peps/pep-0204.html">http://www.python.org/peps/pep-0204.html</a></pre>
Colin W.<br>
<br>
<h1>[Python-Dev] Conditional Expression Resolution</h1>
<b>Guido van Rossum</b> <a
 href="mailto:python-dev%40python.org?Subject=%5BPython-Dev%5D%20Conditional%20Expression%20Resolution&amp;In-Reply-To="
 title="[Python-Dev] Conditional Expression Resolution">guido at
python.org </a><br>
<i>Fri Sep 30 03:21:53 CEST 2005</i>
<ul>
  <li>Previous message: <a
 href="http://mail.python.org/pipermail/python-dev/2005-September/056853.html">[Python-Dev]
[PATCH][BUG] Segmentation Fault in xml.dom.minidom.parse
    </a></li>
  <li>Next message: <a
 href="http://mail.python.org/pipermail/python-dev/2005-September/056847.html">[Python-Dev]
Conditional Expression Resolution
    </a></li>
  <li> <b>Messages sorted by:</b> <a
 href="http://mail.python.org/pipermail/python-dev/2005-September/date.html#56846">[
date ]</a> <a
 href="http://mail.python.org/pipermail/python-dev/2005-September/thread.html#56846">[
thread ]</a> <a
 href="http://mail.python.org/pipermail/python-dev/2005-September/subject.html#56846">[
subject ]</a> <a
 href="http://mail.python.org/pipermail/python-dev/2005-September/author.html#56846">[
author ]</a> </li>
</ul>
<hr><!--beginarticle-->
<pre>After a long discussion I've decided to add a shortcut conditional
expression to Python 2.5.

The syntax will be ...

    A if C else B</pre>
<br>
</body>
</html>

Re: [Numpy-discussion] optimizing power() for complex and real cases

[David M. C. <co...@ph...>]

On Sat, Feb 18, 2006 at 06:17:47PM -0700, Tim Hochberg wrote:
> 
> OK, I now have a faily clean implementation in C 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)
> 
> 
> First a couple of technical questions, then on to the philosophical portion 
> of this note.
> 
> 1. Is there a nice fast way to get a matrix filled with ones from C. I've 
> been tempted to write a ufunc 'ones_like', but I'm afraid that might be 
> considered inappropriate.
> 
> 2. Are people aware that array_power is sometimes passed non arrays as its 
> first argument? Despite having the signature:
> 
> array_power(PyArrayObject *a1, PyObject *o2) 
> 
> This caused me almost no end of headaches, not to mention crashes during 
> numpy.test().

Yes; because it's the implementation of __pow__, the second argument can
be anything.

> I'll check this into the power_optimization branch RSN, hopefully with a 
> fix for the zero power case. Possibly also after extending it to inplace 
> power as well.
> 
> 
> OK, now on to more important stuff. As I've been playing with this my 
> opinion has gone in circles a couple of times. I now think the issue of 
> optimizing integer powers of complex numbers and integer powers of floats 
> are almost completely different. Because complex powers are quite slow and 
> relatively inaccurate, it is appropriate to optimize them for integer 
> powers at the level of nc_pow. This should be just a matter of liberal 
> borrowing from complexobject.c, but I haven't tried it yet.

Ok.

> On the other hand, real powers are fast enough that doing anything at the 
> single element level is unlikely to help. So in that case we're left with 
> either optimizing the cases where the dimension is zero as David has done, 
> or optimizing at the __pow__ (AKA array_power) level as I've done now based 
> on David's original suggestion. This second approach is faster because it 
> avoids the mysterious python scalar -> zero-D array conversion overhead. 
> However, it suffers if we want to optimize lots of different powers since 
> one needs a ufunc for each one. So the question becomes, which powers 
> should we optimize?

Hmm, ufuncs are passed a void* argument for passing info to them. Now,
what that argument is defined when the ufunc is created, but maybe
there's a way to piggy-back on it.

> My latest thinking on this is that we should optimize only those cases 
> where the optimized result is no less accurate than that produced by pow. 
> I'm going to assume that all C operations are equivalently accurate, so 
> pow(x,2) has roughly the same amount of error as x*x. (Something on the 
> order of 0.5 ULP I'd guess). In that case:
>   pow(x, -1) -> 1 / x
>   pow(x, 0) -> 1
>   pow(x, 0.5) -> sqrt(x)
>   pow(x, 1) -> x
>   pow(x, 2) -> x*x
> can all be implemented in terms of multiply or divide with the same 
> accuracy as the original power methods. Once we get beyond these, the error 
> will go up progressively.
> 
> The minimal set described above seems like it should be relatively 
> uncontroversial and it's what I favor. Once we get beyond this basic set, 
> we would need to reach some sort of consensus on how much additional error 
> we are willing to tolerate for optimizing these extra cases. You'll notice 
> that I've changed my mind, yet again, over whether to optimize A**0.5. 
> Since the set of additional ufuncs needed in this case is relatively small, 
> just square and inverse (==1/x), this minimal set works well if optimizing 
> in pow as I've done.

Ok. I'm still not happy with the speed of pow(), though. I'll have to
sit and look at. We may be able to optimize integer powers better.
And there's another area: integer powers of integers. Right now that
uses pow(), whereas we might be able to do better. I'm looking into
that. A framework for that could be helpful for the complex powers too.

Too bad we couldn't make a function generator :-) [Well, we could using
weave...]

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

Re: [Numpy-discussion] complex division

[David M. C. <co...@ph...>]

On Sun, Feb 19, 2006 at 03:32:34PM -0700, Tim Hochberg wrote:
> 
> While rummaging around Python's complexobject.c looking for code to 
> steal for complex power, I came across the following comment relating to 
> complex division:
> 
>        /******************************************************************
>        This was the original algorithm.  It's grossly prone to spurious
>        overflow and underflow errors.  It also merrily divides by 0 despite
>        checking for that(!).  The code still serves a doc purpose here, as
>        the algorithm following is a simple by-cases transformation of this
>        one:
> 
>        Py_complex r;
>        double d = b.real*b.real + b.imag*b.imag;
>        if (d == 0.)
>            errno = EDOM;
>        r.real = (a.real*b.real + a.imag*b.imag)/d;
>        r.imag = (a.imag*b.real - a.real*b.imag)/d;
>        return r;
>        ******************************************************************/
> 
>        /* This algorithm is better, and is pretty obvious:  first
>    divide the
>         * numerators and denominator by whichever of {b.real, b.imag} has
>         * larger magnitude.  The earliest reference I found was to CACM
>         * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
>         * University).  As usual, though, we're still ignoring all IEEE
>         * endcases.
>         */
> 
> The algorithm shown, and maligned, in this comment is pretty much 
> exactly what is done in numpy at present. The function goes on to use 
> the improved algorithm, which I will include at the bottom of the post. 
> It seems nearly certain that using this algorithm will result in some 
> speed hit, although I'm not certain how much. I will probably try this 
> out at some point and see what the speed hit, but in case I drop the 
> ball I thought I'd throw this out there as something we should at least 
> look at. In most cases, I'll take accuracy over raw speed (within reason).
> 
> -tim

The condition for accuracy on this is
|Z - z| < epsilon |z|

where I'm using Z for the computed value of z=a/b, and epsilon is on the
order of machine accuracy.

As pointed out by Steward (ACM TOMS, v. 11, pg 238 (1985)), this doesn't
mean that the real and imaginary components are accurate. The example he
gives is a = 1e70 + 1e-70i and b=1e56+1e-56i, where z=a/b=1e14 + 1e-99i,
which is susceptible to underflow for a machine with 10 decimal digits
and a exponent range of +-99.

Priest (ACM TOMS v30, pg 389 (2004)) gives an alternative, which I won't
show here, b/c it does bit-twiddling with the double representation. But
it does a better job of handling overflow and underflow in intermediate
calculations, is competitive in terms of accuracy, and is faster (at
least on a 750 MHz UltraSPARC-III ;) than the other algorithms except
for the textbook version. One problem is the sample code is for double
precision; for single or longdouble, we'd have to figure out some magic
constants.

Maybe I'll look into it later, but for now Smith's algorithm is better
than the textbook one we were using :-)

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

Re: [Numpy-discussion] complex division

[Charles R H. <cha...@gm...>]

Hmm...

The new algorithm does look better with respect to overflow and
underflow, but I wonder if it is not a bit of overkill. It seems to me
that the same underflow/overflow problems attend complex
multiplication, which is pretty much all that goes on in the original
algorithm.

One thing I do know is that division is expensive. I wonder if one
division and two multiplications might be cheaper than two divisions.
I'll have to check that out.

Chuck

On 2/19/06, Tim Hochberg <tim...@co...> wrote:
>
> While rummaging around Python's complexobject.c looking for code to
> steal for complex power, I came across the following comment relating to
> complex division:
>
>         /****************************************************************=
**
>         This was the original algorithm.  It's grossly prone to spurious
>         overflow and underflow errors.  It also merrily divides by 0 desp=
ite
>         checking for that(!).  The code still serves a doc purpose here, =
as
>         the algorithm following is a simple by-cases transformation of th=
is
>         one:
>
>         Py_complex r;
>         double d =3D b.real*b.real + b.imag*b.imag;
>         if (d =3D=3D 0.)
>             errno =3D EDOM;
>         r.real =3D (a.real*b.real + a.imag*b.imag)/d;
>         r.imag =3D (a.imag*b.real - a.real*b.imag)/d;
>         return r;
>         *****************************************************************=
*/
>
>         /* This algorithm is better, and is pretty obvious:  first
>     divide the
>          * numerators and denominator by whichever of {b.real, b.imag} ha=
s
>          * larger magnitude.  The earliest reference I found was to CACM
>          * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
>          * University).  As usual, though, we're still ignoring all IEEE
>          * endcases.
>          */
>
> The algorithm shown, and maligned, in this comment is pretty much
> exactly what is done in numpy at present. The function goes on to use
> the improved algorithm, which I will include at the bottom of the post.
> It seems nearly certain that using this algorithm will result in some
> speed hit, although I'm not certain how much. I will probably try this
> out at some point and see what the speed hit, but in case I drop the
> ball I thought I'd throw this out there as something we should at least
> look at. In most cases, I'll take accuracy over raw speed (within reason)=
.
>
> -tim
>
>
>
>
>
>
>
>          Py_complex r;    /* the result */
>           const double abs_breal =3D b.real < 0 ? -b.real : b.real;
>          const double abs_bimag =3D b.imag < 0 ? -b.imag : b.imag;
>
>          if (abs_breal >=3D abs_bimag) {
>              /* divide tops and bottom by b.real */
>              if (abs_breal =3D=3D 0.0) {
>                  errno =3D EDOM;
>                  r.real =3D r.imag =3D 0.0;
>              }
>              else {
>                  const double ratio =3D b.imag / b.real;
>                  const double denom =3D b.real + b.imag * ratio;
>                  r.real =3D (a.real + a.imag * ratio) / denom;
>                  r.imag =3D (a.imag - a.real * ratio) / denom;
>              }
>         }
>         else {
>             /* divide tops and bottom by b.imag */
>             const double ratio =3D b.real / b.imag;
>             const double denom =3D b.real * ratio + b.imag;
>             assert(b.imag !=3D 0.0);
>             r.real =3D (a.real * ratio + a.imag) / denom;
>             r.imag =3D (a.imag * ratio - a.real) / denom;
>         }
>         return r;
>
>
>
>
> -------------------------------------------------------
> 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] complex division

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

While rummaging around Python's complexobject.c looking for code to 
steal for complex power, I came across the following comment relating to 
complex division:

        /******************************************************************
        This was the original algorithm.  It's grossly prone to spurious
        overflow and underflow errors.  It also merrily divides by 0 despite
        checking for that(!).  The code still serves a doc purpose here, as
        the algorithm following is a simple by-cases transformation of this
        one:

        Py_complex r;
        double d = b.real*b.real + b.imag*b.imag;
        if (d == 0.)
            errno = EDOM;
        r.real = (a.real*b.real + a.imag*b.imag)/d;
        r.imag = (a.imag*b.real - a.real*b.imag)/d;
        return r;
        ******************************************************************/

        /* This algorithm is better, and is pretty obvious:  first
    divide the
         * numerators and denominator by whichever of {b.real, b.imag} has
         * larger magnitude.  The earliest reference I found was to CACM
         * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
         * University).  As usual, though, we're still ignoring all IEEE
         * endcases.
         */

The algorithm shown, and maligned, in this comment is pretty much 
exactly what is done in numpy at present. The function goes on to use 
the improved algorithm, which I will include at the bottom of the post. 
It seems nearly certain that using this algorithm will result in some 
speed hit, although I'm not certain how much. I will probably try this 
out at some point and see what the speed hit, but in case I drop the 
ball I thought I'd throw this out there as something we should at least 
look at. In most cases, I'll take accuracy over raw speed (within reason).

-tim







         Py_complex r;    /* the result */
          const double abs_breal = b.real < 0 ? -b.real : b.real;
         const double abs_bimag = b.imag < 0 ? -b.imag : b.imag;

         if (abs_breal >= abs_bimag) {
             /* divide tops and bottom by b.real */
             if (abs_breal == 0.0) {
                 errno = EDOM;
                 r.real = r.imag = 0.0;
             }
             else {
                 const double ratio = b.imag / b.real;
                 const double denom = b.real + b.imag * ratio;
                 r.real = (a.real + a.imag * ratio) / denom;
                 r.imag = (a.imag - a.real * ratio) / denom;
             }
        }
        else {
            /* divide tops and bottom by b.imag */
            const double ratio = b.real / b.imag;
            const double denom = b.real * ratio + b.imag;
            assert(b.imag != 0.0);
            r.real = (a.real * ratio + a.imag) / denom;
            r.imag = (a.imag * ratio - a.real) / denom;
        }
        return r;

Re: [Numpy-discussion] Re: number ranges (was Re: Matlab page on scipy wiki)

[Alan G I. <ai...@am...>]

On Sun, 19 Feb 2006, "Colin J. Williams" apparently wrote: 
> The conditional expression (PEP 308) was rejected but is currently being 
> implemented in Python version 2.5.  It will have the syntax A if C else B. 

It's coming
http://www.python.org/peps/pep-0308.html
But in 2.5?
http://www.python.org/dev/doc/devel/whatsnew/whatsnew25.html

Thank you,
Alan Isaac