Re: modified residues() for matlab compatibility

[Ben Abbott]

On Oct 8, 2007, at 3:41 PM, John W. Eaton wrote:

> On  6-Oct-2007, Ben Abbott wrote:
>
> | I've attached a modified residue.m that includes modifications to  
> the
> | doc-string and the  addition of test scripts. In addition, the
> | multiplicity vector, "e", which the prior change dropped has been
> | reinstated.
>
> OK, I applied this patch, but had to do some of it by hand since there
> had been intervening formatting changes.  Would you please check the
> current CVS version and make sure that it is correct?  In particular,
> I am seeing the failure:
>
>   octave:1> test residue
>     ***** test
>    b = [1, 0, 1];
>    a = [1, 0, 18, 0, 81];
>    [r, p, k, e] = residue(b, a);
>    assert ((abs (54*r - [-5i; 12; 5i; 12]) < 1e-6
>    && abs (p - [3i; 3i; -3i; -3i]) < 1e-7
>    && isempty (k)
>    && e == [1; 2; 1; 2]));
>    [br, ar] = residue (r, p, k);
>    assert ((abs (br - b) < 1e-12
>    && abs (ar - a) < 1e-12));
>   !!!!! test failed
>   error: assert ((abs (54 * r - [-5i; 12; 5i; 12]) < 1e-6 && abs (p  
> - [3i; 3i; -3i; -3i]) < 1e-7 && isempty (k) && e == [1; 2; 1; 2]))  
> failed
>   octave:2> b = [1, 0, 1];
>   octave:3> a = [1, 0, 18, 0, 81];
>   octave:4>  [r, p, k, e] = residue(b, a);
>   octave:5> abs (54*r - [-5i; 12; 5i; 12])
>   ans =
>
>      1.0000e+01
>      4.4460e-08
>      1.0000e+01
>      4.4460e-08
>
> and I'm not sure whether it is the test that is bad, or residue that
> is computing the wrong result.
>
> | The test scripts include the examples give in the doc-string as well
> | as the original error which recently prompted much discussion and
> | lead to my involvement. The accuracy of the residues may have been
> | impacted. I had to lower the increase the permitted error in the  
> test
> | scripts in order for the scripts to pass ... at least relative to  
> the
> | test script in your prior email. Specifically
> |
> |       -%! assert ((abs (r - [-2; 7; 3]) < 1e-7
> |       +%! assert ((abs (r - [-2; 7; 3]) < 1e-5
> |
> | It occurs to me that I need to get set up with CVS so I can see the
> | present state of things. That way I can send patches instead of
> | sources ... unfortunately the web-site, http://www.octave.org/
> | download.html, is down. I'll get to that task when it's back up.
>
> I should be up again.
>
> | In any event, I've attached a new residue.m as well as the  
> equivalent
> | patch (which may have a unix/dos EOL problem). I've also attached a
> | ChangeLog.
> |
> | Some questions,
> |
> | (1) Is there a naming convention for patches?
>
> No the name doesn't matter, but it is helpful if you send them as
> plain text attachments so I can apply them directly from Emacs instead
> of having to save them to a file.
>
> | (2) Which variety of EOL is desired, or does it matter?
>
> The files should have Unix line endings (LF only), but I don't think
> that matters if you send text/plain attachments.
>
> thanks,
>
> jwe

John,

I did not get the same error. However, I did get an error from the  
final test script placed on the bottom of the file. Since I was  
unable to connect to the CVS over the weekend and I did now know if  
there is a protocol for test-scripts, I'd placed my testing between  
the doc-strip and remainder of the function.

Since having access to the CVS, I've moved the scripts to the bottom  
and reconciled the duplicate.

However, since this is my 1st time using the "test" functionality, it  
might be wise to double check my work.

Meanwhile, there still remains the error you've brought to my  
attention. The problem appears to be with the order in which the  
poles/roots are sorted. The script residue.m calls mpoles.m, which  
uses the internal sort() routine, to sort (reorder) the roots.

It appears that your octave reorders the roots differently than mine.  
I've verified that the mpoles.m checked into the CVS is the same as  
the one I'm using (however, there was a prior version I had posted  
that would produce this problem).

I haven't yet had the time to mirror the entire CVS to my system. So  
I am unable to isolate the problem, and be certain of my guess.

However, I think we can verify the hypothesis quickly. Using my FInk  
install of  Octave, I get the following

        b = [1, 0, 1];
        a = [1, 0, 18, 0, 81];
        [r, p, k, e] = residue(b, a)

r =

   0.00000 - 0.09259i
   0.22222 + 0.00000i
   0.00000 + 0.09259i
   0.22222 - 0.00000i

p =

   0.0000 + 3.0000i
   0.0000 + 3.0000i
  -0.0000 - 3.0000i
  -0.0000 - 3.0000i

k = [](0x0)
e =

   1
   2
   1
   2

On your system, I suspect the poles are ordered as

p =

   0.0000 - 3.0000i
   0.0000 - 3.0000i
  -0.0000 + 3.0000i
  -0.0000 + 3.0000i

i.e. the -3i pair comes ahead of the +3i pair.

Can you verify?

If the problem is that simple, then I'll venture a guess that the  
sort function has been modified since the 2.9.14 release ... or ...  
your CPU is giving slightly different answers than mine. In either  
event, it makes since to include both results in the test-script  
(since they are each actually correct).

Assuming my assumptions are correct, I've modified the script. The  
patch is below. Please let me know if the syntax for creating the  
patch is not what you prefer.

diff -Naur residue_cvs.m residue.m

--- residue_cvs.m       2007-10-08 21:10:06.000000000 -0400
+++ residue.m   2007-10-08 22:06:07.000000000 -0400
@@ -122,31 +122,6 @@
 ## Created: June 1994
 ## Adapted-By: jwe

-%!test
-%! b = [1, 1, 1];
-%! a = [1, -5, 8, -4];
-%! [r, p, k, e] = residue (b, a);
-%! assert ((abs (r - [-2; 7; 3]) < 1e-5
-%! && abs (p - [2; 2; 1]) < 1e-7
-%! && isempty (k)
-%! && e == [1; 2; 1]));
-%! k = [1 0];
-%! [b, a] = residue (r, p, k);
-%! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12
-%! && abs (a - [1, -5, 8, -4]) < 1e-12));
-
-%!test
-%! b = [1, 0, 1];
-%! a = [1, 0, 18, 0, 81];
-%! [r, p, k, e] = residue(b, a);
-%! assert ((abs (54*r - [-5i; 12; 5i; 12]) < 1e-6
-%! && abs (p - [3i; 3i; -3i; -3i]) < 1e-7
-%! && isempty (k)
-%! && e == [1; 2; 1; 2]));
-%! [br, ar] = residue (r, p, k);
-%! assert ((abs (br - b) < 1e-12
-%! && abs (ar - a) < 1e-12));
-
 function [r, p, k, e] = residue (b, a, varargin)

   if (nargin < 2 || nargin > 3)
@@ -323,12 +298,27 @@

 endfunction

-%% test/octave.test/poly/residue-1.m
 %!test
 %! b = [1, 1, 1];
 %! a = [1, -5, 8, -4];
 %! [r, p, k, e] = residue (b, a);
-%! assert((abs (r - [-2; 7; 3]) < 1e-6
+%! assert ((abs (r - [-2; 7; 3]) < 1e-5
 %! && abs (p - [2; 2; 1]) < 1e-7
 %! && isempty (k)
 %! && e == [1; 2; 1]));
+%! k = [1 0];
+%! [b, a] = residue (r, p, k);
+%! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12
+%! && abs (a - [1, -5, 8, -4]) < 1e-12));
+
+%!test
+%! b = [1, 0, 1];
+%! a = [1, 0, 18, 0, 81];
+%! [r, p, k, e] = residue(b, a);
+%! assert ((abs (54*r - [-5i; 12; 5i; 12]) < 1e-6
+%! && abs (p - [3i; 3i; -3i; -3i]) < 1e-7
+%! && isempty (k)
+%! && e == [1; 2; 1; 2]));
+%! [br, ar] = residue (r, p, k);
+%! assert ((abs (br - b) < 1e-12
+%! && abs (ar - a) < 1e-12));

Ben