bugGNU Octave - Bugs: bug #51148, residue function potentially...

 
 

bug #51148: residue function potentially inaccurate

Submitter:  None
Submitted:  Wed 31 May 2017 12:08:58 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Frank Wiedmann Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 4.2.1
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 09 Jun 2017 06:44:24 PM UTC, comment #22: 

Agreed.  If someone can open a new bug report with a small example that shows clearly whether the problem is with deconv or filter that would be good.

Rik <rik5>
Group administrator
Fri 09 Jun 2017 11:38:18 AM UTC, comment #21: 

@Anonymous: yes, I can confirm the behaviour in Matlab. This is why I said that there is maybe a problem either in deconv or in filter.
Anyway, the problem in residue should be fixed and now a new bug report should be opened.

Marco Caliari <caliari>
Group Member
Fri 09 Jun 2017 11:21:49 AM UTC, comment #20: 

By the way, with MATLAB, I get two leading zeros for both of your experiments.

Anonymous
Fri 09 Jun 2017 10:47:37 AM UTC, comment #19: 

I agree that comparing to 0 does invite problems (although the code has worked fine for me so far in MATLAB). However, even if I replace the line by the one that you suggested, I still get an error:


>> foster2cauer_test
error: 'foster2cauer_test' undefined near line 1 column 1
>> foster2cauer_test
error: div(2): out of bound 1
error: called from
    foster2cauer at line 31 column 11
    foster2cauer_test at line 7 column 22
>>


Anonymous
Fri 09 Jun 2017 10:33:04 AM UTC, comment #18: 

Although in general it is not safe something like


rem = rem(find(rem~=0,1):end);


and should be replaced by (e.g.)


rem = rem(find(abs(rem)>eps/2,1):end);


now the problem is maybe in deconv (or even in filter). For instance


[div,rem]=deconv([1.92306958582241e+15,3.20449986572221e+24,1.34271290136344e+32,2.32739765751038e+38],[7.33727670161595e+27,1.05919311870816e+36,4.56169848520627e+42])


gives rem with two leading zeros, while


[div,rem]=deconv([1.92306958582241e+15,3.204499865722215e+24,1.34271290136344e+32,2.32739765751038e+38],[7.33727670161595e+27,1.05919311870816e+36,4.56169848520627e+42])


(I added a 5 as last digit of second element of first input) gives rem with one leading zero.

Marco Caliari <caliari>
Group Member
Fri 09 Jun 2017 09:23:15 AM UTC, comment #17: 

Thank you all for working on this. I just edited my residue.m file (by hand) to apply the changes from residue.diff. Unfortunately, when I run my original testcase, I still get an error:


>> foster2cauer_test
error: foster2cauer: operator +: nonconformant arguments (op1 is 1x4, op2 is 1x5)
error: called from
    foster2cauer at line 33 column 7
    foster2cauer_test at line 7 column 22
>>


Could someone please try to run my testcase with the patched version of residue.m in order to see if there is still a problem remaining or if I made a mistake when I edited my version of residue.m?

Anonymous
Thu 08 Jun 2017 07:00:12 PM UTC, comment #16: 

@Marco: I pushed your patch to the stable branch here (http://hg.savannah.gnu.org/hgweb/octave/rev/5faf1767f9a9).  I changed the BIST syntax to


%!test <51148>


so that Octave will automatically recognize that this test is associated with bug #51148.

I also did some benchmarking and the new method is slower by about 20%.  However, the average running time increased from 3.4 milliseconds to 4.0 milliseconds.  Almost no human will notice the 0.6 milliseconds, while they will notice the incorrectly calculated values.

Rik <rik5>
Group administrator
Thu 08 Jun 2017 07:00:04 AM UTC, comment #15: 

@Rik: I tried different normalizations and to use the residual from deconv, but no success. I attach my patch wrt stable.

(file #40879)

Marco Caliari <caliari>
Group Member
Tue 06 Jun 2017 09:41:45 PM UTC, comment #14: 

There are two outputs from deconv, the polynomial and the residual.  Would keeping track of the residual and adding it back in at the next iteration of the for loop help?

Rik <rik5>
Group administrator
Tue 06 Jun 2017 09:35:07 PM UTC, comment #13: 

No better idea at the moment, this is not really my mathematical area.

However, as a general rule Octave prefers accuracy over performance.  Once the calculation is correct one can go back and see if there is a way to make it faster.

Three for loops is definitely going to be slow.  The current conv() and deconv() functions make use of filter() which is written in C++ and is quite fast.  Is there a way to scale/normalize the pden and pm inputs to deconv so that


pn = deconv (pden, pm);


computes correctly?


Rik <rik5>
Group administrator
Tue 06 Jun 2017 12:51:12 PM UTC, comment #12: 

The problem is in the subfunction rresidue in residue.m. It is dangerous to compute pn as deconvolution of pden and current pm (around line 340). In fact, pm is the product of all the other denominators. If you replace


pn = deconv (pden, pm);


with


pn = 1;
for j = indx(1):n - 1
  pn = conv (pn, [1, -p(j)]);
end
for j = n + 1:indx(end)
  pn = conv (pn, [1, -p(j)]);
end
for j = 1:e(n) - 1
  pn = deconv (pn, [1, -p(n)]);
end


you get the right solution. I'm still not satisfied with the last deconvolution step. Any better idea?

Marco Caliari <caliari>
Group Member
Thu 01 Jun 2017 09:31:19 PM UTC, comment #11: 

I just examined the effect of different permutations of the input data on the result of the residue function (by removing the semicolon after the residue function and commenting out the rest of the foster2cauer function). I found that in two cases, I do not get


p =

  1.0000e+018  9.0036e+027  1.1289e+036  1.7010e+047

q =

  1.0000e+000  1.9231e+015  3.2045e+024  1.3427e+032  2.3274e+038


but


p =

  1.0000e+018  9.0036e+027  1.1294e+036  -7.5446e+047

q =

  1.0000e+000  1.9231e+015  3.2045e+024  1.3427e+032  2.3274e+038


This reinforces the impression that there is something wrong here.

Anonymous
Thu 01 Jun 2017 04:46:52 PM UTC, comment #10: 

Re-titling bug report to reflect what appears to be an underlying problem in residue.m.  I verified things by changing foster2cauer.m to use the p and q values calculated by Matlab:


% calculate polynomial coefficients of impedance (numerator and denominator)
%[p,q] = residue(1./cf,-1./(rf.*cf),0);

p = [1.0000e+18   9.0036e+27   1.1290e+36   4.6827e+42];
q = [ 1.0000e+00   1.9231e+15   3.2045e+24   1.3427e+32   2.3274e+38];


With this change, the test case foster2cauer_test runs through to completion.

I stopped the script during execution and have saved the input values to the residue function in a file residue.var which is attached.  To test things run


load residue.var
[p,q] = residue (r, pin, k)


This has the bad result for p(4) which is in error by 10^5.


(file #40836)

Rik <rik5>
Group administrator
Thu 01 Jun 2017 01:02:58 PM UTC, comment #9: 

Thanks for looking into this. Shouldn't we set the status to "Confirmed" instead of "Fixed"?

Anonymous
Thu 01 Jun 2017 12:53:26 PM UTC, comment #8: 

I did spend 2hrs on this and my conclusion is that residue is broken! I have not looked any deeper.

Doug Stewart <dastew>
Thu 01 Jun 2017 07:43:13 AM UTC, comment #7: 

Matlab returns the following for the residue calculation:

p =

   1.0000e+18   9.0036e+27   1.1290e+36   4.6827e+42


q =

   1.0000e+00   1.9231e+15   3.2045e+24   1.3427e+32   2.3274e+38

p(4) is clearly different, which is probably the cause for the problem.


Anonymous
Wed 31 May 2017 08:34:52 PM UTC, comment #6: 

For a quick test, I removed the semicolon from the call to residue and I get


p =

   1.0000e+18   9.0036e+27   1.1289e+36   1.7010e+47

q =

   1.0000e+00   1.9231e+15   3.2045e+24   1.3427e+32   2.3274e+38


It would be interesting to see what Matlab returns for the residue calculation.

The actual error in Octave is caused by the unequal sizes of div and rem after the call to deconv.


q = div(2)*p+[zeros(1,length(p)-length(rem)) rem];


The variable p is 1x4 while the variable rem is 1x5.

Rik <rik5>
Group administrator
Wed 31 May 2017 08:26:58 PM UTC, comment #5: 

I'll let Doug look at this more when he has time.

Rik <rik5>
Group administrator
Wed 31 May 2017 08:15:58 PM UTC, comment #4: 

Did you run my testcase with the fixed residue.m? Like I wrote in comment #2, when I commented out the trimming in residue.m, I still got an error for my testcase in Octave, albeit a different one.

Anonymous
Wed 31 May 2017 06:16:07 PM UTC, comment #3: 

I grafted the changeset for bug #49291 on to the stable branch (http://hg.savannah.gnu.org/hgweb/octave/rev/7eba3cd5b6e0).  It will be a part of a bug fix release 4.2.2 if there is one, or it will be a part of the next major release 4.4.0.

In the meantime, this is just a change to an m-file so you can grab the residue.m code from the Mercurial repository and overwrite the old m-file in your local directory to see the change.

Rik <rik5>
Group administrator
Wed 31 May 2017 01:47:00 PM UTC, comment #2: 

The problem you mentioned was fixed as bug #49291 (also reported by me). However, the fix is apparently not yet available in the official release. When I comment out the trimming in residue.m, I still get an error:

>> foster2cauer_test

error: foster2cauer: operator +: nonconformant arguments (op1 is 1x4, op2 is 1x5)
error: called from
    foster2cauer at line 33 column 7
    foster2cauer_test at line 7 column 22

>>

Anonymous
Wed 31 May 2017 01:33:18 PM UTC, comment #1: 

I tried you test and also got the same error.
The problem is in residue.m
Residue thinks that the numbers are too big or small and discards some.
when I re-scaled your numbers I got a better answer.

I don't have time now but will look at it more latter.

Doug Stewart <dastew>
Wed 31 May 2017 12:08:58 PM UTC, original submission:  

I have written a short function that calculates the equivalent Cauer RC network for a given Foster RC network. The test script shows that, as expected, the result is independent (within machine precision) of the order of the RC pairs in the Foster network. When I try to run the test script in Octave, I get an error.

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #40879:  residue.diff added by caliari (2KiB - application/x-tex)
file #40836:  residue.var added by rik5 (175B - application/octet-stream)
file #40824:  foster2cauer.m added by None (2KiB - text/plain)
file #40825:  foster2cauer_test.m added by None (289B - text/plain)
file #40826:  foster2cauer_test_matlab.txt added by None (6KiB - text/plain)
file #40827:  foster2cauer_test_octave.txt added by None (829B - text/plain)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 15 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-06-08 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2017-06-08 caliari Attached File- Added residue.diff, #40879
    2017-06-01 rik5 Attached File- Added residue.var, #40836
        Item GroupMatlab Compatibility Inaccurate Result
        StatusFixed Confirmed
        Open/ClosedClosed Open
        Operating SystemMicrosoft Windows Any
        SummaryFunction that runs in Matlab produces error in Octave residue function potentially inaccurate
    2017-05-31 rik5 StatusNone Fixed
        Open/ClosedOpen Closed
    2017-05-31 None Attached File- Added foster2cauer.m, #40824
        Attached File- Added foster2cauer_test.m, #40825
        Attached File- Added foster2cauer_test_matlab.txt, #40826
        Attached File- Added foster2cauer_test_octave.txt, #40827

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code