[Octave-bug-tracker] [bug #52919] comparison of complex values

[Rik]

Update of bug #52919 (project octave):

                 Release:                   4.0.3 => dev                    

    _______________________________________________________

Follow-up Comment #14:

I think we all agree that sort and sortrows have performance issues that
should be addressed.  That is really its own topic and should be segregated
into a new bug report.  Could someone create that and port the comments over
from this report?

As for Matlab's behavior, it really was quite bad.  As far as I know, the
initial implementation was *not* a stable sort.  They added an extra option to
commands like unique() so you could specify whether you wanted a "stable" or
"sorted" output (http://www.mathworks.com/help/matlab/ref/unique.html).

Real numbers are not a special case of complex numbers in Matlab.  They are a
different class.  The functions isreal, iscomplex, and isinteger all test the
class that is used to represent the number, rather than testing the actual
value of the number.  For example


x = complex (1,0)
x =  1 + 0i
octave:9> isreal (x)
ans = 0
octave:10> iscomplex (x)
ans = 1
octave:11> isinteger (x)
ans = 0


The number 1.0 is clearly an integer, but Octave says it isn't because it is
being stored using 8 bytes and an IEEE double format.

The issue that you ran in to with indexing of a complex array producing a real
result is known as narrowing.  Because complex numbers take twice as much
memory to store, Octave and Matlab attempt to reduce complex arrays to real
arrays whenever they can.  For example, you cannot creat a complex scalar with
a zero magnitude imaginary component on the command line.  Try


octave:12> x = 1+0i
x =  1


In this case Octave's parser notices that you are trying to create a complex
number, but the interpreter realizes that it can narrow this complex number to
a real number and it does so.  The way around that is to use the complex()
function which avoids automatic narrowing.


octave:13> x = complex (1,0)
x =  1 + 0i


Even so, Octave will narrow this to a real result whenever it has a chance. 
An operation like indexing, or addition, will do this.  For example,


octave:13> x = complex (1,0)
x =  1 + 0i
octave:14> x(1)
ans =  1
octave:15> x
x =  1 + 0i
octave:16> x + 0
ans =  1


Matlab has some documentation to help out on this:
https://www.mathworks.com/help/coder/ug/code-generation-for-complex-data.html.

Note:


For code generation, complex data that has all zero-valued imaginary parts
remains complex. This data does not become real. This behavior has the
following implications:

    In some cases, results from functions that sort complex data by absolute
value can differ from the MATLAB® results. See Functions That Sort Complex
Values by Absolute Value.

    For functions that require that complex inputs are sorted by absolute
value, complex inputs with zero-valued imaginary parts must be sorted by
absolute value. These functions include ismember, union, intersect, setdiff,
and setxor.


Because Octave was second, we have to follow Matlab's rules which are that one
sort algorithm is used for real numbers (value) and another algorithm for
complex numbers (absolute value).

The extra complication around complex numbers on the negative real axis is not
a buglet.  Matlab is (now) breaking ties in absolute value by using the phase
angle which they get by calling the underlying atan2 library function.  This
function returns an angle in the range [-pi, pi] (both endpoints are
included).  But this causes problems because it can assign two values to the
same underlying real number.  For example,


octave:21> x = [3-4i, complex(-5,0), 3+4i]
x =

   3 - 4i  -5 + 0i   3 + 4i

octave:22> abs (x)
ans =

   5   5   5

octave:23> angle (x)
ans =

  -0.92730   3.14159   0.92730


All numbers have the same magnitude, so the tie is broken by phase angle
resulting in


[3-4i, 3+4i, -5]


But the IEEE standard for floating points includes both +0 and -0.  What
happens when the imaginary part is -0


x(2) = complex (-5,-0)
x =

   3 - 4i  -5 - 0i   3 + 4i

octave:25> abs (x)
ans =

   5   5   5

octave:26> angle (x)
ans =

  -0.92730  -3.14159   0.92730


Now the sort order is


[-5, 3-4i, 3+4i]


That bothers most people, so Octave takes the additional step of remapping the
atan2 function on to the range (-pi, pi] (first endpoint is not included).

The question to me is whether, given that we have to use abs() as the first
test, it would be useful to prefer lexigraphic order during ties.  The current
sort order puts negative real values after positive real values which is
consistent, but doesn't generally accord with people's natural sorting
algorithm.  For example,


sort ([complex(-5,0), complex(5,0), 3+4i, 3-4i])
ans =

   3 - 4i   5 + 0i   3 + 4i  -5 + 0i



We could have -5+0i appear first which might make more sense.  There are
potentially other benefits as well.  It may be faster to simply compare
magnitudes of real and imaginary parts rather than calculating the atan2
function.  Finally, there is an issue on Mac OSX platforms that makes the
remapping of the range [-pi,pi] to (-pi, pi] cumbersome because they use a
different value of pi than anyone else.  If the code were just comparing
magnitudes it should be platform independent.

It appears that there could be some benefits to re-defining the sort
algorithm, but first the Octave Maintainer's should decide whether it is worth
deviating from Matlab compatibility in this regard.



    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?52919>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/