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!
|
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.
|
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?
|
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.
|
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.
|
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:
- The issue reported originally, when passed an array with columns of complex pairs, the results for all but the first column are incorrect
- When passsing a tolerance value of 1 or greater, matlab returns an error
- Values that cplxpair determines to be real (i.e. have an imaginary part less than tol * abs(z)) are forced to be real
- 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.
- 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)
|
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)
|
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)
|
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.
|
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.
|
Fri 27 Jan 2017 04:30:40 AM UTC, comment #15:
|
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?
|
Fri 27 Jan 2017 03:45:33 AM UTC, comment #13:
Re: the FIXME question in Comment #7
matlab 2016b:
|
Fri 27 Jan 2017 03:42:14 AM UTC, comment #12:
Re: Comment #7
Matlab 2016b:
|
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).
|
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.
|
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
|
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).
|
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?
|
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.
|
Mon 23 Jan 2017 03:55:58 AM UTC, comment #5:
Matlab 2016b output:
|
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:
|
Mon 23 Jan 2017 03:41:47 AM UTC, comment #3:
matlab 2016b output for the test script:
|
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)
|
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:
|
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.
|