Inconsistent numpy.angle behavior

[rubenvb]

Following up on a Stackoverflow question of mine:
http://stackoverflow.com/questions/15137003/strange-complex-phase-displayed

I am bringing this here. The issue is that angle called on a zero-matrix causes inconsistent values of the phase of a complex number. I understand the underlying float implementation is the cause, but I think it is important to understand e.g. Matlab does not have this problem, and instead does the more sensible thing and returns 0.

For completeness, the numpy script:

from __future__ import division
from numpy import *
from matplotlib.pyplot import *

Radius=10
N=1024
dx=2*Radius/N
dy=dx
x=r_[-Radius:Radius:dx]
y=r_[-Radius:Radius:dy]
X, Y = meshgrid(x,y)
R = sqrt(X**2+Y**2)
PHI = arctan2(Y,X)

ringthing = R < Radius
ring = zeros((2,N,N),dtype=complex)
ring[0] = ringthing
ring[1] = ringthing*exp(1j*PHI)

f=fig()
p1=f.add_subplot(121)
p1.imshow(angle(ring[0]))
p2=f.add_subplot(122)
p2.imshow(angle(ring[1]))

f.show()

And the corresponding Matlab code:

N=64;
Radius=10;
dx=2*Radius/(N-1);
x=-Radius:dx:Radius
y=-Radius:dx:Radius;
[X,Y] = meshgrid(x,y);
[X,Y] = meshgrid(x,y);
[PHI,R] = cart2pol(X,Y);
ringthing = (R<Radius).*(R>.75*Radius);
imshow(angle(exp(1j*PHI).*ringthing))

I agree the phase of a zero isn't strictly defined (and GNU Octave... yuck, handles it even worse than numpy), but when showing these types of images, this is a hindrance and having to code around this seems artificial (and adds an extra line of code ;-)). If it has to be undefined, make it at least numerically handy undefined.

[ewmoore]

Internally angle calls numpy.arctan2, which really just is the atan2 function from libm. This wrapping to pi for -0+0j is a special value of atan2 as defined in the C standard, (c.f. www.open-std.org/jtc1/sc22/wg14/www/docs/n1124.pdf, section F.9.1.4.)

I'll note that python's cmath.phase, which is the python equivalent to numpy.angle, also returns, pi. I don't think this is really an error.

[njsmith]

[I added from __future__ import division to the test code, since otherwise it doesn't work on python 2]

A simpler example showing what's going on here is:

def signed_complex_zero(neg_real, neg_imag):
    c = np.empty((), dtype=complex)
    c.real = np.copysign(0, 1 - 2 * neg_real)
    c.imag = np.copysign(0, 1 - 2 * neg_imag)
    return c[()]
signed_zeros = [signed_complex_zero(False, False), signed_complex_zero(True, False), signed_complex_zero(True, True), signed_complex_zero(False, True)]
for zero in signed_zeros:
    print("angle(%r) = %r" % (zero, np.angle(zero)))

prints:

angle((-0-0j)) = -3.1415926535897931
angle((-0+0j)) = 3.1415926535897931
angle(-0j) = -0
angle(0j) = 0.0

I don't think np.angle(0) can really return NaN, because that breaks the equivalence x = abs(x) * exp(angle(x) * 1j); it means that you can't actually reliably round-trip from real/imag to magnitude/angle representations. (Well, this round-tripping won't preserve the signs on zeros even now, but oh well.) And as @ewmoore notes, the proper return values for atan2 given any combination of positive or negative zeros is defined by Annex F.

[ewmoore]

Let's close this. It's not an error.

[pv]

Closing, the current behavior is correct.

[Yadunandanaacharya]

print(np.angle((-0-0j)))

This is returning output : 0.0.
Even I'm trying get angle from FFT signal, but returning zero.

[eric-wieser]

@Yadunandanaacharya, your issue is that -0 is 0, you have to use -0.0:

In [5]: np.angle((-0-0j))
Out[5]: 0.0

In [6]: np.angle((-0.0-0j))
Out[6]: 3.141592653589793

This is because only python floats distinguish -0 and 0, but 0 is a python int - this is a python issue, not a numpy one,

[eric-wieser]

I opened https://bugs.python.org/issue41485 about this, the builtin complex.__repr__ is being very misleading.