numpy-discussion

[Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[<pe...@ce...>]

Hi,

I have recieved the following note from a user:

"""
In SciPy 0.3.x the ufuncs were overloaded by more "intelligent" versions.
A very attractive feature was that sqrt(-1) would yield 1j as in Matlab.
Then you can program formulas directly (e.g., roots of a 2nd order
polynomial) and the right answer is always achieved. In the Matlab-Python
battle in mathematics education, this feature is important.

Now in SciPy 0.5.x sqrt(-1) yields nan. A lot of code we have, especially
for introductory numerics and physics courses, is now broken.
This has already made my colleagues at the University skeptical to
Python as "this lack of backward compatibility would never happen in Matlab".

Another problem related to Numeric and numpy is that in these courses we
use ScientificPython several places, which applies Numeric and will
continue to do so. You then easily get a mix of numpy and Numeric
in scripts, which may cause problems and at least extra overhead.
Just converting to numpy in your own scripts isn't enough if you call
up libraries using and returning Numeric.
"""

I wonder, what are the reasons that numpy.sqrt(-1) returns nan?
Could sqrt(-1) made to return 1j again? If not, shouldn't
numpy.sqrt(-1) raise a ValueError instead of returning silently nan?

Thanks,
Pearu

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

pe...@ce... wrote:

>Hi,
>
>I have recieved the following note from a user:
>
>"""
>In SciPy 0.3.x the ufuncs were overloaded by more "intelligent" versions.
>A very attractive feature was that sqrt(-1) would yield 1j as in Matlab.
>Then you can program formulas directly (e.g., roots of a 2nd order
>polynomial) and the right answer is always achieved. In the Matlab-Python
>battle in mathematics education, this feature is important.
>
>Now in SciPy 0.5.x sqrt(-1) yields nan. A lot of code we have, especially
>for introductory numerics and physics courses, is now broken.
>This has already made my colleagues at the University skeptical to
>Python as "this lack of backward compatibility would never happen in Matlab".
>  
>
This was a consequence of moving scipy_base into NumPy but not exposing 
the scimath library in NumPy.   It would be a very easy thing to put 
from numpy.lib.scimath import *
into the scipy name-space.

I'm supportive of that as a backward-compatibility measure.

>Another problem related to Numeric and numpy is that in these courses we
>use ScientificPython several places, which applies Numeric and will
>continue to do so. You then easily get a mix of numpy and Numeric
>in scripts, which may cause problems and at least extra overhead.
>Just converting to numpy in your own scripts isn't enough if you call
>up libraries using and returning Numeric.
>  
>
>"""
>
>I wonder, what are the reasons that numpy.sqrt(-1) returns nan?
>  
>
Because that is backwards compatible.  You have to construct a 
function-wrapper in order to handle the negative case correctly.  The 
function wrapper is going to be slower.  Thus, it is placed in a 
separate library.

>Could sqrt(-1) made to return 1j again? 
>
Not in NumPy.  But, in scipy it could.

>If not, shouldn't
>  
>
>numpy.sqrt(-1) raise a ValueError instead of returning silently nan?
>  
>
This is user adjustable.  You change the error mode to raise on 
'invalid' instead of pass silently which is now the default.

-Travis

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[Sven S. <sve...@gm...>]

Travis Oliphant schrieb:

> 
>> If not, shouldn't
>>  
>>
>> numpy.sqrt(-1) raise a ValueError instead of returning silently nan?
>>  
>>
> This is user adjustable.  You change the error mode to raise on 
> 'invalid' instead of pass silently which is now the default.
> 
> -Travis
> 

Could you please explain how this adjustment is done, or point to the
relevant documentation.
Thank you,
Sven

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

Sven Schreiber wrote:

>>This is user adjustable.  You change the error mode to raise on 
>>'invalid' instead of pass silently which is now the default.
>>
>>-Travis
>>
>>    
>>
>
>Could you please explain how this adjustment is done, or point to the
>relevant documentation.
>  
>

numpy.sqrt(-1)

old = seterr(invalid='raise')
numpy.sqrt(-1)  # should raise an error

seterr(**old)  # restores error-modes for current thread
numpy.sqrt(-1)

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

Travis Oliphant wrote:
> Sven Schreiber wrote:
>
>   
>>> This is user adjustable.  You change the error mode to raise on 
>>> 'invalid' instead of pass silently which is now the default.
>>>
>>> -Travis
>>>
>>>    
>>>
>>>       
>> Could you please explain how this adjustment is done, or point to the
>> relevant documentation.
>>  
>>
>>     
>
> numpy.sqrt(-1)
>
> old = seterr(invalid='raise')
> numpy.sqrt(-1)  # should raise an error
>
> seterr(**old)  # restores error-modes for current thread
> numpy.sqrt(-1)
>
>   
With python 2.5 out now, perhaps it's time to come up with a with 
statement context manager. Something like:

    from __future__ import with_statement
    import numpy

    class errstate(object):
        def __init__(self, **kwargs):
            self.kwargs = kwargs
        def __enter__(self):
            self.oldstate = numpy.seterr(**self.kwargs)
        def __exit__(self, *exc_info):
            numpy.seterr(**self.oldstate)
           
    a = numpy.arange(10)
    a/a # ignores divide by zero
    with errstate(divide='raise'):
        a/a # raise exception on divide by zer
    # Would ignore divide by zero again if we got here.

-tim

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

Tim Hochberg wrote:

>With python 2.5 out now, perhaps it's time to come up with a with 
>statement context manager. Something like:
>
>    from __future__ import with_statement
>    import numpy
>
>    class errstate(object):
>        def __init__(self, **kwargs):
>            self.kwargs = kwargs
>        def __enter__(self):
>            self.oldstate = numpy.seterr(**self.kwargs)
>        def __exit__(self, *exc_info):
>            numpy.seterr(**self.oldstate)
>           
>    a = numpy.arange(10)
>    a/a # ignores divide by zero
>    with errstate(divide='raise'):
>        a/a # raise exception on divide by zer
>    # Would ignore divide by zero again if we got here.
>
>-tim
>
>  
>

This looks great.  I think most people aren't aware of the with 
statement and what it can do (I'm only aware because of your posts, for 
example). 

So, what needs to be added to your example in order to just add it to 
numpy?

-Travis

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

Travis Oliphant wrote:
> Tim Hochberg wrote:
>
>   
>> With python 2.5 out now, perhaps it's time to come up with a with 
>> statement context manager. Something like:
>>
>>    from __future__ import with_statement
>>    import numpy
>>
>>    class errstate(object):
>>        def __init__(self, **kwargs):
>>            self.kwargs = kwargs
>>        def __enter__(self):
>>            self.oldstate = numpy.seterr(**self.kwargs)
>>        def __exit__(self, *exc_info):
>>            numpy.seterr(**self.oldstate)
>>           
>>    a = numpy.arange(10)
>>    a/a # ignores divide by zero
>>    with errstate(divide='raise'):
>>        a/a # raise exception on divide by zer
>>    # Would ignore divide by zero again if we got here.
>>
>> -tim
>>
>>  
>>
>>     
>
> This looks great.  I think most people aren't aware of the with 
> statement and what it can do (I'm only aware because of your posts, for 
> example). 
>
> So, what needs to be added to your example in order to just add it to 
> numpy?
>   
As far as I know, just testing  and documentation -- however testing was 
so minimal that I may find some other stuff. I'll try to clean it up 
tomorrow so that I'm a little more confident that it works correctly and 
I'll send another note out then.


-tim

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[Andrew J. <a.h...@gm...>]

Tim Hochberg wrote:
> Travis Oliphant wrote:
>> Tim Hochberg wrote:
>>
>>   
>>> With python 2.5 out now, perhaps it's time to come up with a with 
>>> statement context manager. Something like:
>>>
>>>    a = numpy.arange(10)
>>>    a/a # ignores divide by zero
>>>    with errstate(divide='raise'):
>>>        a/a # raise exception on divide by zer
>>>    # Would ignore divide by zero again if we got here.
>>>
>>> -tim
>>>
>> This looks great.  I think most people aren't aware of the with 
>> statement and what it can do (I'm only aware because of your posts, for 
>> example). 
>>
>> So, what needs to be added to your example in order to just add it to 
>> numpy?
>>   
> As far as I know, just testing  and documentation -- however testing was 
> so minimal that I may find some other stuff. I'll try to clean it up 
> tomorrow so that I'm a little more confident that it works correctly and 
> I'll send another note out then.

For this particular application, why not a context manager which just 
substitutes in the appropriately-optimized version of sqrt? That is, 
don't change the error state, but actually change the value of the 
object pointed at by the name sqrt?

Andrew

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[Fernando P. <fpe...@gm...>]

On 10/11/06, Travis Oliphant <oli...@ee...> wrote:

> pe...@ce... wrote:
> >Could sqrt(-1) made to return 1j again?
> >
> Not in NumPy.  But, in scipy it could.

Without taking sides on which way to go, I'd like to -1 the idea of a
difference in behavior between numpy and scipy.

IMHO, scipy should be within reason a strict superset of numpy.
Gratuitious differences in behavior like this one are going to drive
us all mad.

There are people who import scipy for everything, others distinguish
between numpy and scipy, others use numpy alone and at some point in
their life's code they do

import numpy as N -> import scipy as N

because they start needing stuff not in plain numpy.  Having different
APIs and behaviors appear there is, I think, a Seriously Bad Idea
(TM).

Just my 1e-2j,

Cheers,

f

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[Stefan v. d. W. <st...@su...>]

On Wed, Oct 11, 2006 at 03:37:34PM -0600, Fernando Perez wrote:
> On 10/11/06, Travis Oliphant <oli...@ee...> wrote:
>=20
> > pe...@ce... wrote:
> > >Could sqrt(-1) made to return 1j again?
> > >
> > Not in NumPy.  But, in scipy it could.
>=20
> Without taking sides on which way to go, I'd like to -1 the idea of a
> difference in behavior between numpy and scipy.
>=20
> IMHO, scipy should be within reason a strict superset of numpy.
> Gratuitious differences in behavior like this one are going to drive
> us all mad.
>=20
> There are people who import scipy for everything, others distinguish
> between numpy and scipy, others use numpy alone and at some point in
> their life's code they do
>=20
> import numpy as N -> import scipy as N
>=20
> because they start needing stuff not in plain numpy.  Having different
> APIs and behaviors appear there is, I think, a Seriously Bad Idea
> (TM).

I agree with Fernando on this one.

Further, if I understand correctly, changing sqrt and power to give
the right answer by default will slow things down somewhat.  But is it
worth sacrificing intuitive usage for speed?

N.power(2,-2) =3D=3D 0

and

N.sqrt(-1) =3D=3D nan

just doesn't feel right.  Why not then have

N.power(2,-2) =3D=3D 0.24
N.sqrt(-1) =3D=3D 1j

and write a special function that does fast calculation of
square-roots for positive values?

Cheers
St=E9fan

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

Stefan van der Walt wrote:

>I agree with Fernando on this one.
>
>Further, if I understand correctly, changing sqrt and power to give
>the right answer by default will slow things down somewhat.  But is it
>worth sacrificing intuitive usage for speed?
>  
>
For NumPy, yes. 

This is one reason that NumPy by itself is not a MATLAB replacement. 

>N.power(2,-2) == 0
>
>and
>
>N.sqrt(-1) == nan
>
>just doesn't feel right.  
>

Only because your expectations are that NumPy *be* a MATLAB 
replacement.  The problem is that it sacrifices too much for that to be 
the case.   And we all realize that NumPy needs more stuff added to it 
to be like IDL/MATLAB such as SciPy, Matplotlib, IPython, etc.

>Why not then have
>
>N.power(2,-2) == 0.24
>N.sqrt(-1) == 1j
>
>and write a special function that does fast calculation of
>square-roots for positive values?
>  
>

We've already done this.  The special functions are called

numpy.power
numpy.sqrt

(notice that if you do numpy.sqrt(-1+0j) you get the "expected" answer 
emphasizing that numpy does no "argument" checking to determine the output).

The "intuitive" functions (which must do argument checking) are (in 
numpy.lib.scimath) but exported as

scipy.power (actually I need to check that one...)
scipy.sqrt

What could be simpler?  ;-)

-Travis

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[Stefan v. d. W. <st...@su...>]

On Wed, Oct 11, 2006 at 05:21:44PM -0600, Travis Oliphant wrote:
> Stefan van der Walt wrote:
>=20
> >I agree with Fernando on this one.
> >
> >Further, if I understand correctly, changing sqrt and power to give
> >the right answer by default will slow things down somewhat.  But is it
> >worth sacrificing intuitive usage for speed?
> > =20
> >
> For NumPy, yes.=20
>=20
> This is one reason that NumPy by itself is not a MATLAB
> replacement.=20

Intuitive usage is hopefully not a MATLAB-only feature.

> >N.power(2,-2) =3D=3D 0
> >
> >and
> >
> >N.sqrt(-1) =3D=3D nan
> >
> >just doesn't feel right. =20
> >
>=20
> Only because your expectations are that NumPy *be* a MATLAB=20
> replacement.  The problem is that it sacrifices too much for that to be=
=20
> the case.   And we all realize that NumPy needs more stuff added to it=20
> to be like IDL/MATLAB such as SciPy, Matplotlib, IPython, etc.

I have none such expectations -- I havn't used MATLAB in over 5 years.
All I'm saying is that, since the value of the square-root of -1 is
not nan and 2^(-2) is not 0, it doesn't surprise me that this
behaviour confuses people.

> The "intuitive" functions (which must do argument checking) are (in=20
> numpy.lib.scimath) but exported as
>=20
> scipy.power (actually I need to check that one...)
> scipy.sqrt
>=20
> What could be simpler?  ;-)

I'm sure this is going to come back and haunt us.  We have two
libraries, one depends on and exposes the API of the other, yet it
also overrides some functions with its own behaviour, while keeping
the same names.

I'll shut up now :)

Cheers
St=E9fan

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[Greg W. <gre...@gm...>]

On 10/11/06, Travis Oliphant <oli...@ee...> wrote:
>
> Stefan van der Walt wrote:
> >Further, if I understand correctly, changing sqrt and power to give
> >the right answer by default will slow things down somewhat.  But is it
> >worth sacrificing intuitive usage for speed?
> >
> For NumPy, yes.
>
> This is one reason that NumPy by itself is not a MATLAB replacement.


This is not about being a Matlab replacement.
This is about correctness.
Numpy purports to handle complex numbers.
Mathematically, sqrt(-1) is a complex number.
Therefore Numpy *must* return a complex number.
Speed should not take precedence over correctness.

If Numpy doesn't return a complex number then it shouldn't pretend to
support complex numbers.

Greg

-- 
Linux.  Because rebooting is for adding hardware.

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

On 10/12/06, Greg Willden <gre...@gm...> wrote:
> On 10/11/06, Travis Oliphant <oli...@ee...> wrote:
> > Stefan van der Walt wrote:
> > >Further, if I understand correctly, changing sqrt and power to give
> > >the right answer by default will slow things down somewhat.  But is it
> > >worth sacrificing intuitive usage for speed?
> > >
> > For NumPy, yes.
> >
> > This is one reason that NumPy by itself is not a MATLAB replacement.
>
> This is not about being a Matlab replacement.
> This is about correctness.
> Numpy purports to handle complex numbers.
> Mathematically, sqrt(-1) is a complex number.
> Therefore Numpy *must* return a complex number.
> Speed should not take precedence over correctness.

Unless your goal is speed.  Then speed should take precedence over correctness.

And unless you're a fan of quaternions, in which case *which* square
root of -1 should it return?

It's interesting to note that although python has had complex numbers
pretty much from the beginning, math.sqrt(-1) returns an error.  If
you want to work with complex square roots you need to use
cmath.sqrt().  Basically you have to tell python that compex numbers
are something you care about by using the module designed for complex
math.  This scimath module is a similar deal.

But perhaps the name could be a little more obvious / short?  Right
now it seems it only deals with complex numbers, so maybe having
"complex" or "cmath" in the name would make it clearer.  Hmm there is
a numpy.math, why not a numpy.cmath.

> If Numpy doesn't return a complex number then it shouldn't pretend to
> support complex numbers.

That's certainly being overdramatic.  Lot's of folks are doing nifty
stuff with complex FFTs every day using numpy/scipy.  And I'm sure
they will continue to no matter what numpy.sqrt(-1) returns.

--bb

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[Greg W. <gre...@gm...>]

On 10/11/06, Bill Baxter <wb...@gm...> wrote:
>
> On 10/12/06, Greg Willden <gre...@gm...> wrote:
> > Speed should not take precedence over correctness.
>
> Unless your goal is speed.  Then speed should take precedence over
> correctness.


Huh?
Person 1:  "Hey you should use function X."
Person 2:  "No, it doesn't give me the correct answer."
Person 1: "Who cares?  It's fast!"

What kind of logic is that?

I'm having serious doubts about my conversion to Numpy.
Greg

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[James G. <jg...@ca...>]

Greg Willden wrote:
> On 10/11/06, *Bill Baxter* <wb...@gm... 
> <mailto:wb...@gm...>> wrote:
> 
>     On 10/12/06, Greg Willden <gre...@gm...
>     <mailto:gre...@gm...>> wrote:
>      > Speed should not take precedence over correctness.
> 
>     Unless your goal is speed.  Then speed should take precedence over
>     correctness. 
> 
> 
> Huh?
> Person 1:  "Hey you should use function X."
> Person 2:  "No, it doesn't give me the correct answer."
> Person 1: "Who cares?  It's fast!"
> 
> What kind of logic is that?
> 
> I'm having serious doubts about my conversion to Numpy.

It's more like
Person 1: Function x is designed to be fast, have predictable memory usage and 
so be suitable for a wide variety of applications. Alternatively there is 
function y which is slower but handles a larger numerical domain. y Might be 
more sutiable for your needs; you should use it.
Person 2: Numpy is broken! Function x should just do what I want!

Really, given that both options are available, I don't see a big issue here. 
Claims that returning j for sqrt(-1) is "more correct" seem highly dubious as a) 
NumPy is a numerical library, not an algebraic one, hence the fact that it 
behaves in a type-dependent way is not so surprising, and b) no one complains 
that functions like sqrt return only positive floats, despite the fact that, 
mathematically, this function is multivalued.

-- 
"Eternity's a terrible thought. I mean, where's it all going to end?"
  -- Tom Stoppard, Rosencrantz and Guildenstern are Dead

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[Stefan v. d. W. <st...@su...>]

On Thu, Oct 12, 2006 at 08:58:21AM -0500, Greg Willden wrote:
> On 10/11/06, Bill Baxter <wb...@gm...> wrote:
>=20
>     On 10/12/06, Greg Willden <gre...@gm...> wrote:
>     > Speed should not take precedence over correctness.
>=20
>     Unless your goal is speed.  Then speed should take precedence over
>     correctness.
>=20
>=20
> Huh?
> Person 1:  "Hey you should use function X."
> Person 2:  "No, it doesn't give me the correct answer."
> Person 1: "Who cares?  It's fast!"
>=20
> What kind of logic is that?

I tried to explain the argument at

http://www.scipy.org/NegativeSquareRoot

Cheers
St=E9fan

Re: [Numpy-discussion] =?iso-8859-1?q?Should_numpy=2Esqrt=28-1=29_ret?= =?iso-8859-1?q?urn_1j_rather_than=09nan=3F?=

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

On Thu, 12 Oct 2006, Stefan van der Walt apparently wrote:=20
> I tried to explain the argument at=20
> http://www.scipy.org/NegativeSquareRoot=20

Helpful.  But you start off by saying:
        In mathematics, the above assumption is true -- that=20
        the square root of -1 is 1j.

Since square root is (here) a function, and part of the=20
function definition is the domain and codomain, this=20
statement is not correct.  If the codomain is real numbers,=20
the domain must correspondingly be (a subset of) the=20
nonnegative reals.  The nan output is the result of an=20
**invalid** input.

So the question is really: "Why is a negative real number an=20
invalid input in this implementation", or "Why in this=20
implementation is the type of the output restricted by the=20
type of the input?" You get a good start on answering that.

Cheers,
Alan Isaac

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than?nan?

[Stefan v. d. W. <st...@su...>]

On Thu, Oct 12, 2006 at 10:53:12AM -0400, Alan G Isaac wrote:
> On Thu, 12 Oct 2006, Stefan van der Walt apparently wrote:=20
> > I tried to explain the argument at=20
> > http://www.scipy.org/NegativeSquareRoot=20
>=20
> Helpful.  But you start off by saying:
>         In mathematics, the above assumption is true -- that=20
>         the square root of -1 is 1j.
>=20
> Since square root is (here) a function, and part of the=20
> function definition is the domain and codomain, this=20
> statement is not correct.  If the codomain is real numbers,=20
> the domain must correspondingly be (a subset of) the=20
> nonnegative reals.  The nan output is the result of an=20
> **invalid** input.
>=20
> So the question is really: "Why is a negative real number an=20
> invalid input in this implementation", or "Why in this=20
> implementation is the type of the output restricted by the=20
> type of the input?" You get a good start on answering that.

You are more than welcome to change the wiki if you can think of a
simpler way to explain the problem... In fact, I would encourage that
-- I just put it there as a starting point.

Cheers
St=E9fan

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

Greg Willden wrote:
> On 10/11/06, *Travis Oliphant* <oli...@ee... 
> <mailto:oli...@ee...>> wrote:
>
>     Stefan van der Walt wrote:
>     >Further, if I understand correctly, changing sqrt and power to give
>     >the right answer by default will slow things down somewhat.  But
>     is it
>     >worth sacrificing intuitive usage for speed?
>     >
>     For NumPy, yes.
>
>     This is one reason that NumPy by itself is not a MATLAB replacement.
>
>
> This is not about being a Matlab replacement.
> This is about correctness.
I disagree.   NumPy does the "correct" thing when you realize that sqrt 
is a function that returns the same type as it's input.   The field 
over-which the operation takes place is defined by the input data-type 
and not the input "values".  Either way can be considered correct 
mathematically.   As Paul said it was a design decision not to go 
searching through the array to determine whether or not there are 
negative numbers in the array.

Of course you can do that if you want and that's what scipy.sqrt does. 
> Numpy purports to handle complex numbers.
> Mathematically, sqrt(-1) is a complex number.
Or, maybe it's undefined if you are in the field of real numbers.  It 
all depends.
> Therefore Numpy *must* return a complex number.
Only if the input is complex.  That is a reasonable alternative to your 
specification.
>
> If Numpy doesn't return a complex number then it shouldn't pretend to 
> support complex numbers.
Of course it supports complex numbers, it just doesn't support automatic 
conversion to complex numbers.    It supports complex numbers the same 
way Python supports them (i.e. you have to use cmath to get sqrt(-1) == 1j)

People can look at this many ways without calling the other way of 
looking at it unreasonable. 

I don't see a pressing need to change this in NumPy, and in fact see 
many reasons to leave it the way it is.   This discussion should move to 
the scipy list because that is the only place where a change could occur.

-Travis.

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

Greg Willden wrote:
> On 10/11/06, *Travis Oliphant* <oli...@ee... 
> <mailto:oli...@ee...>> wrote:
>
>     Stefan van der Walt wrote:
>     >Further, if I understand correctly, changing sqrt and power to give
>     >the right answer by default will slow things down somewhat.  But
>     is it
>     >worth sacrificing intuitive usage for speed?
>     >
>     For NumPy, yes.
>
>     This is one reason that NumPy by itself is not a MATLAB replacement.
>
>
> This is not about being a Matlab replacement.
> This is about correctness.
> Numpy purports to handle complex numbers.
> Mathematically, sqrt(-1) is a complex number.
That's vastly oversimplified. If you are working with the reals, then 
sqrt(-1) is undefined (AKA nan). If you are working in the complex 
plane, then sqrt(-1) is indeed *a* complex number; of course you don't 
know *which* complex number it is unless you also specify the branch.
> Therefore Numpy *must* return a complex number.
No I don't think that it must. I've found it a very useful tool for the 
past decade plus without it returning complex numbers from sqrt.
> Speed should not take precedence over correctness.
The current behavior is not incorrect.
>
> If Numpy doesn't return a complex number then it shouldn't pretend to 
> support complex numbers.
Please relax.

Personally I think that the default error mode should be tightened up. 
Then people would only see these sort of things if they really care 
about them. Using Python 2.5 and the errstate class I posted earlier:

    # This is what I like for the default error state
    numpy.seterr(invalid='raise', divide='raise', over='raise',
    under='ignore')

    a = -numpy.arange(10)
    with errstate(invalid='ignore'):
        print numpy.sqrt(a) # This happily returns a bunch of NANs; and
    one zero.
    print numpy.sqrt(a.astype(complex)) # This returns a bunch of
    complex values
    print numpy.sqrt(a) # This raise a floating point error. No silent
    NANs returned

This same error state make the vagaries of dividing by zero less 
surprising as well.

-tim

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

On 10/11/06, Tim Hochberg <tim...@ie...> wrote:
>
> Greg Willden wrote:
> > On 10/11/06, *Travis Oliphant* <oli...@ee...
> > <mailto:oli...@ee...>> wrote:
> >
> >     Stefan van der Walt wrote:
> >     >Further, if I understand correctly, changing sqrt and power to give
> >     >the right answer by default will slow things down somewhat.  But
> >     is it
> >     >worth sacrificing intuitive usage for speed?
> >     >
> >     For NumPy, yes.
> >
> >     This is one reason that NumPy by itself is not a MATLAB replacement.
> >
> >
> > This is not about being a Matlab replacement.
> > This is about correctness.
> > Numpy purports to handle complex numbers.
> > Mathematically, sqrt(-1) is a complex number.
> That's vastly oversimplified. If you are working with the reals, then
> sqrt(-1) is undefined (AKA nan). If you are working in the complex
> plane, then sqrt(-1) is indeed *a* complex number;


And if you are working over the rationals, sqrt(-1) and sqrt(-2) lead to
different field extensions ;) Of course, numpy doesn't *have* rationals, so
I'm just being cute.


<snip>

Personally I think that the default error mode should be tightened up.
> Then people would only see these sort of things if they really care
> about them. Using Python 2.5 and the errstate class I posted earlier:
>
>     # This is what I like for the default error state
>     numpy.seterr(invalid='raise', divide='raise', over='raise',
>     under='ignore')


I like these choices too. Overflow, division by zero, and sqrt(-x) are
usually errors, indicating bad data or programming bugs. Underflowed floats,
OTOH, are just really, really small numbers and can be treated as zero. An
exception might be if the result is used in division and no error is raised,
resulting in a loss of accuracy.

If complex results are *expected*, then this should be made explicit by
using complex numbers. Numpy allows finegrained control of data types and
array ordering, it is a bit closer to the metal than Matlab. This extra
control allows greater efficiency, both in speed and in storage, at the cost
of a bit more care on the programmers side.

Chuck

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

>
>     Personally I think that the default error mode should be tightened
>     up.
>     Then people would only see these sort of things if they really care
>     about them. Using Python 2.5 and the errstate class I posted earlier:
>
>         # This is what I like for the default error state
>         numpy.seterr (invalid='raise', divide='raise', over='raise',
>         under='ignore')
>
>
> I like these choices too. Overflow, division by zero, and sqrt(-x) are 
> usually errors, indicating bad data or programming bugs. Underflowed 
> floats, OTOH, are just really, really small numbers and can be treated 
> as zero. An exception might be if the result is used in division and 
> no error is raised, resulting in a loss of accuracy.
>

I'm fine with this.  I've hesitated because error checking is by default 
slower.  But, I can agree that it is "less surprising" to new-comers.  
People that don't mind no-checking can simply set their defaults back to 
'ignore'


-Travis

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

[<pe...@ce...>]

PS: I am still sending this message to numpy list only because the
    proposal below affects numpy code, not scipy one.

I think Fernando points make sense, numpy.foo(x) != scipy.foo(x) can
cause confusion and frustration both among new numpy/scipy users and
developers (who need to find explanations for the choises made).

So, let me propose the following solution so that all parties will get
the same results without sacrifying numpy.sqrt speed on non-negative input
and scipy.sqrt backward compatibility:

Define numpy.sqrt as follows:

def sqrt(x): 
    r = nx.sqrt(x)
    if nx.nan in r:
        i = nx.where(nx.isnan(r))
        r = _tocomplex(r)
        r[i] = nx.sqrt(_tocomplex(x[i]))
    return r

and define

  numpy.sqrtp

that takes only non-negative input, this is for those users who expect
sqrt to fail on negative input (as Numeric.sqrt and math.sqrt do).

Pearu

Re: [Numpy-discussion] Should numpy.sqrt(-1) return 1j rather than nan?

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

Travis Oliphant wrote:
>>     Personally I think that the default error mode should be tightened
>>     up.
>>     Then people would only see these sort of things if they really care
>>     about them. Using Python 2.5 and the errstate class I posted earlier:
>>
>>         # This is what I like for the default error state
>>         numpy.seterr (invalid='raise', divide='raise', over='raise',
>>         under='ignore')
>>
>>
>> I like these choices too. Overflow, division by zero, and sqrt(-x) are 
>> usually errors, indicating bad data or programming bugs. Underflowed 
>> floats, OTOH, are just really, really small numbers and can be treated 
>> as zero. An exception might be if the result is used in division and 
>> no error is raised, resulting in a loss of accuracy.
>>
>>     
>
> I'm fine with this.  I've hesitated because error checking is by default 
> slower.  But, I can agree that it is "less surprising" to new-comers.  
> People that don't mind no-checking can simply set their defaults back to 
> 'ignore'
>
>   
Great.

One thing we may want to do (numarray had this), was add a pseudo 
argument 'all', that allows you to set all of the values at once. Then 
if you want the full-bore, ignore-all-errors computation (and your using 
2.5 and "from __future__ import with_statement") you can just do:

with errstate(all='ignore'):
    # computation here
# back to being picky

-tim