bugGNU Octave - Bugs: bug #50124, cplxpair produces incorrect result...

 
 

bug #50124: cplxpair produces incorrect result for second column of complex pairs

Submitted by:  Bill Lash <welash>
Submitted on:  Mon 23 Jan 2017 12:07:27 AM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: FixedAssigned to: None
Originator Name: William LashOpen/Closed: Closed
Release: 4.0.3Operating System: GNU/Linux

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission Jump to the original submission)

Wed 30 Aug 2017 12:37:24 AM UTC, comment #25:

I tested the script with the changes that Rik made, and all of my test cases pass. I think it can be closed. Thanks!

Bill Lash <welash>
Tue 29 Aug 2017 10:46:41 PM UTC, comment #24:

I made several changes to update the documentation and follow Octave coding conventions here (http://hg.savannah.gnu.org/hgweb/octave/rev/a315ac23dc6c).

The only thing of note was this change

This was done to be analogous to

If we don't need the call to ones(), because we are relying on implicit broadcasting, then it can be removed from this line as well.

Rik <rik5>
Project Administrator
Tue 29 Aug 2017 09:08:37 PM UTC, comment #23:

I pushed the change here:

http://hg.savannah.gnu.org/hgweb/octave/rev/b4cf41d173dd

Can this report be closed?

John W. Eaton <jwe>
Project Administrator
Mon 17 Jul 2017 05:50:55 PM UTC, comment #22:

Nothing to be done. It just needs someone to do the final review and push. But it's the middle of summer and even Octave volunteers take time off. Also, there has been some really deep re-structuring of the Octave code base around using C++11 features and that has sucked up a lot of developer time.

Rik <rik5>
Project Administrator
Sun 16 Jul 2017 01:31:12 AM UTC, comment #21:

I happened to think about this the other day, and notice that the patch has not yet been applied. Is there anything I need to do? I think the latest version obeys coding standards and has tests.

Bill Lash <welash>
Sun 29 Jan 2017 05:38:29 AM UTC, comment #20:

Thanks for matlab-izing the test script, I always forget that octave has a superset of the matlab syntax.

So as far as I can tell, these are the issues where matlab and octave disagree:

  1. The issue reported originally, when passed an array with columns of complex pairs, the results for all but the first column are incorrect
  1. When passsing a tolerance value of 1 or greater, matlab returns an error
  1. Values that cplxpair determines to be real (i.e. have an imaginary part less than tol * abs(z)) are forced to be real
  1. Complex pairs are forced to be true complex conjugates. It appears that Matlab chooses the member of the pair that has a positive imaginary part, and conjugates it to get the member of the pair with the negative part. The pairs are sorted with the member that has negative real part listed first.
  1. The way that matlab seems to be comparing values to determine if they are complex pairs is by comparing the conjugate of each value to the other values, and finding the one with the minimum magnitude of the difference. If this minimum difference is less than tol * abs(z) then it is the pair.

I tried my hand at adressing these issues, and I think I have version that works. I made the changes against the latest version in mercurial. I am attaching a patch (new_cplxpair.patch) as well as the results of running your version of the script (latest.out). There are some formatting differences (matlab gives one more digit than octave in format long, and in some cases octave prints 1 or 0 instead of 1.000000000000000), but it seems to match up other than that.

I did have to change some of the tests, due to forcing of complex pairs to be conjugates causing test failures.

(file #39599, file #39600)

Bill Lash <welash>
Fri 27 Jan 2017 02:25:39 PM UTC, comment #19:

for future reference, had to cleanup the script to run in Matlab:

no printf in matlab.

printf --> fprintf

" --> '

endfor --> end

end_try_catch --> end

"Undefined function 'printf' for input arguments of type 'char'"

can't do abs(z)(1)

absz = abs(z)
absz(1)

anyway, after that it ran and I don't think I changed any intended function:

i attached the matlab compatible version as tcplxpair_matlab.m

(file #39570)

Nicholas Jankowski <nrjank>
Fri 27 Jan 2017 06:33:42 AM UTC, comment #18:

Should there also be a dependency on bug #47865 https://savannah.gnu.org/bugs/?47865 that was mentioned earlier, which seems to be trying to address the tolerance issues.

If matlab supports try/catch I came up with a test script the ranges over a set of "error" values for real numbers and complex pairs to try to determine how the tolerance argument is used in each case.

It also does a few tests to determine how the results are sorted and which value of the complex pair is used when forcing the complex pairs to be equal. I am attaching the script, called tcplxpair, and the results of running it as:

octave --no-gui tcplxpair > output

I ran it with 4.0.3, and I also downloaded the current version of the cplxpair.m file from mercurial (which seems to also be in octave 4.2) and got the results from that.

If someone could try it with matlab and share the results, I'll try to analyze the differences.

(file #39566, file #39567, file #39568)

Bill Lash <welash>
Fri 27 Jan 2017 05:18:48 AM UTC, comment #17:

I added a dependency on bug #45810. That was the bug that motivated the last serious change to the tolerance calculation.

Rik <rik5>
Project Administrator
Fri 27 Jan 2017 05:01:53 AM UTC, comment #16:

So it now seems that there are multiple incompatibilities. First, the fact that the second column of values is not sorted correctly. Second, for Matlab compatibility, results should be exact conjugates. This is documented by The Mathworks:

Third, the tolerance needs to be in the range 0 <= tol <= 1 and probably the algorithm in cplxpair.m updated as well.

Rik <rik5>
Project Administrator
Fri 27 Jan 2017 04:30:40 AM UTC, comment #15:
Nicholas Jankowski <nrjank>
Fri 27 Jan 2017 04:22:04 AM UTC, comment #14:

Can you also try:

From your results it looks like tolerance needs to be between 0 and 1, and is not multiplied by eps internally (since (1-eps)epsabs(z) would not be greater 0.01)

Also, for the FIXME, can you try:

Also, some of this is getting away from the original issue that I reported. Should I create another bug report?

Bill Lash <welash>
Fri 27 Jan 2017 03:45:33 AM UTC, comment #13:

Re: the FIXME question in Comment #7

matlab 2016b:

Nicholas Jankowski <nrjank>
Fri 27 Jan 2017 03:42:14 AM UTC, comment #12:

Re: Comment #7

Matlab 2016b:

Nicholas Jankowski <nrjank>
Fri 27 Jan 2017 01:53:42 AM UTC, comment #11:

is this a recent change in behavior for octave? With 4.0.3 I see:

So for 4.0.3 it looks like cplxpair is comparing to tol*abs(z).

Bill Lash <welash>
Thu 26 Jan 2017 11:10:45 PM UTC, comment #10:

To address Mike's comment: if we don't know Matlab's exact tolerance calculation, then start with a low tolerance. With a low tolerance, the two complex numbers won't pair and Matlab will emit an error. Increase the tolerance until Matlab is able to pair the two complex numbers and then see whether the result is an exactly conjugated pair or the original two complex numbers paired.

To address Bill's comment: The calculated tolerance used by Octave is

For Z == 1+0.5i,

The test metric is

or

Since .01 < 0.022204, they will pair.

If you introduce an absolutely crazy tolerance, like 1e25, then the numbers will be treated as all real. However, this is easily detected because the true real number (1+0i) will not be at the end of the list separated from the complex numbers at the start. Instead, the list ordering will be the same as the input in Octave.

Rik <rik5>
Project Administrator
Thu 26 Jan 2017 10:25:28 PM UTC, comment #9:

Won't setting the tolerance that high mean that all of the data you have given will be treated as real? I would think using a tolerance of 0.1 would be more appropriate. So

Bill Lash <welash>
Thu 26 Jan 2017 10:21:52 PM UTC, comment #8:

We also have the still open question about how that tolerance argument is actually interpreted, since Matlab's docs are very confusing about it (bug #47865).

Mike Miller <mtmiller>
Project Administrator
Thu 26 Jan 2017 10:14:41 PM UTC, comment #7:

The function cplxpair is defined in an m-file. If you look at the start of the code there is a series of FIXME notes.

This seems relevant. Also, with a high enough tolerance you can force Octave to match pairs which are fairly far apart. Octave does pair them, but returns each half of the pair exactly as it appeared in the input. If Matlab is choosing to keep only one half of a pair and conjugate it to get the second half, then this test will immediately show it.

The results for Octave are

@Nicholas: Can you run the test code in Matlab to see what happens?

Rik <rik5>
Project Administrator
Mon 23 Jan 2017 04:07:58 AM UTC, comment #6:

Thanks, so it looks like Matlab is choosing one of the two values of the set of complex pairs, and conjugating it to get the other value.

I wonder if that is a big enough difference to file a matlab compatability bug against.

Bill Lash <welash>
Mon 23 Jan 2017 03:55:58 AM UTC, comment #5:

Matlab 2016b output:

Nicholas Jankowski <nrjank>
Mon 23 Jan 2017 03:52:06 AM UTC, comment #4:

Thanks, that is more what I expected. It is strange that matlab shows a slight difference between the before and after for the first set of complex pairs. Would it be possible to do that again with "format long"?

It seems like it is possible that matlab is forcing the complex pairs to be conjugates of each other. With my proposed patch, and adding a format long at the front of the script, I get:

Bill Lash <welash>
Mon 23 Jan 2017 03:41:47 AM UTC, comment #3:

matlab 2016b output for the test script:

Nicholas Jankowski <nrjank>
Mon 23 Jan 2017 12:43:07 AM UTC, comment #2:

I think the patch below fixes the issue. I am curious if there is a better way to rearrange the z matrix than the way I did it. Basically idx contains columns of numbers between 1 and 7 in my test case, and the intent is to rearrange each column of z according to the corresponding column of idx. I created midx which creates an array of pairs of coordinates.

Comments are welcome.

(file #39534)

Bill Lash <welash>
Mon 23 Jan 2017 12:32:11 AM UTC, comment #1:

oops, the cut and paste lost something in the markup. The test script is actually:

Bill Lash <welash>
Mon 23 Jan 2017 12:07:27 AM UTC, original submission:

If I understand the documentation for cplxpair correctly, if you pass in a matrix with 2 columns, the function should sort the pairs in each column (unless you pass the DIM argument to tell it to sort along the specified dimension).

I am seeing that the second column is filled in with values from the first column. Here is my test script:

z = exp (2ipi[4; 3; 5; 2; 6; 1; 0]/7);
y = [ -1-1i; -1+1i;-3; -2; 1; 2; 3];
[z,y]
cplxpair([z,y])
cplxpair([z,y]) - [z,y]

I get:

>> z = exp (2ipi[4; 3; 5; 2; 6; 1; 0]/7);
>> y = [ -1-1i; -1+1i;-3; -2; 1; 2; 3];
>> [z,y]

ans =

-0.90097 - 0.43388i -1.00000 - 1.00000i
-0.90097 + 0.43388i -1.00000 + 1.00000i
-0.22252 - 0.97493i -3.00000 + 0.00000i
-0.22252 + 0.97493i -2.00000 + 0.00000i
0.62349 - 0.78183i 1.00000 + 0.00000i
0.62349 + 0.78183i 2.00000 + 0.00000i
1.00000 + 0.00000i 3.00000 + 0.00000i

>> cplxpair([z,y])

ans =

-0.90097 - 0.43388i -0.22252 - 0.97493i
-0.90097 + 0.43388i -0.22252 + 0.97493i
-0.22252 - 0.97493i -0.90097 - 0.43388i
-0.22252 + 0.97493i -0.90097 + 0.43388i
0.62349 - 0.78183i 0.62349 - 0.78183i
0.62349 + 0.78183i 0.62349 + 0.78183i
1.00000 + 0.00000i 1.00000 + 0.00000i

>> cplxpair([z,y]) - [z,y]

ans =

0.00000 + 0.00000i 0.77748 + 0.02507i
0.00000 + 0.00000i 0.77748 - 0.02507i
0.00000 + 0.00000i 2.09903 - 0.43388i
0.00000 + 0.00000i 1.09903 + 0.43388i
0.00000 + 0.00000i -0.37651 - 0.78183i
0.00000 + 0.00000i -1.37651 + 0.78183i
0.00000 + 0.00000i -2.00000 + 0.00000i

>>


looking at the code in cplxpair.m, I am suspicious of a couple of things. In:

[idxi, idxj] = find (abs (imag (z)) ./ (abs (z) + realmin (cls)) < tol);
q = sparse (idxi, idxj, 1, n, m);
nr = sum (q, 1);
[q, idx] = sort (q, 1);
z = z(idx);
y = z;

if I take out the semicolons, I see

idx =

1 3
2 4
3 1
4 2
5 5
6 6
7 7

and then

z(idx) is

z =

-0.90097 - 0.43388i -0.22252 - 0.97493i
-0.90097 + 0.43388i -0.22252 + 0.97493i
-0.22252 - 0.97493i -0.90097 - 0.43388i
-0.22252 + 0.97493i -0.90097 + 0.43388i
0.62349 - 0.78183i 0.62349 - 0.78183i
0.62349 + 0.78183i 0.62349 + 0.78183i
1.00000 + 0.00000i 1.00000 + 0.00000i

I think the indexing for the second column needs to be relative to the second column.

Additionally in the loop that does the pairing, the variable j is not used in indexing the y or z array.

Bill Lash <welash>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):
   
   
Comment:
   

Attached Files
file #39599:  new_cplxpair.patch added by welash (3KiB - text/x-patch - new version of patch to address all issues, and test results using this patch)
file #39600:  latest.out added by welash (3KiB - application/octet-stream - new version of patch to address all issues, and test results using this patch)
file #39570:  tcplxpair_matlab.m added by nrjank (4KiB - application/octet-stream - matlab compatible version of testscript)
file #39566:  tcplxpair added by welash (4KiB - application/octet-stream)
file #39567:  results4.2 added by welash (3KiB - application/octet-stream)
file #39568:  results_4.0.3 added by welash (3KiB - application/octet-stream)
file #39534:  cplxpair.patch added by welash (2KiB - text/x-patch)

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by jwe (Posted a comment)
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by nrjank (Posted a comment)
  • -unavailable- added by welash (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 13 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Wed 30 Aug 2017 05:41:57 PM UTCrik5StatusReady For Test=>Fixed
      Open/ClosedOpen=>Closed
    Tue 29 Aug 2017 09:08:37 PM UTCjweStatusPatch Submitted=>Ready For Test
    Sun 29 Jan 2017 05:38:29 AM UTCwelashAttached File-=>Added new_cplxpair.patch, #39599
      Attached File-=>Added latest.out, #39600
    Fri 27 Jan 2017 02:25:39 PM UTCnrjankAttached File-=>Added tcplxpair_matlab.m, #39570
    Fri 27 Jan 2017 06:59:57 AM UTCsiko1056Carbon-CopyRemoved 91577=>-
    Fri 27 Jan 2017 06:33:42 AM UTCwelashAttached File-=>Added tcplxpair, #39566
      Attached File-=>Added results4.2, #39567
      Attached File-=>Added results_4.0.3, #39568
    Fri 27 Jan 2017 05:18:48 AM UTCrik5Dependencies-=>Depends on bugs #45810
    Mon 23 Jan 2017 08:24:36 AM UTCsiko1056StatusNone=>Patch Submitted
    Mon 23 Jan 2017 12:43:07 AM UTCwelashAttached File-=>Added cplxpair.patch, #39534

    Back to the top


    Powered by Savane 3.1-cleanup1