bugGNU Octave - Bugs: bug #37297, Meaningless test in fftfilt.m,...

 
 

bug #37297: Meaningless test in fftfilt.m, patch to address this issue as well enhance final touchup of routine and fix rounding bug

Submitter:  Dan Sebald <sebald>
Submitted:  Fri 07 Sep 2012 09:08:33 PM UTC
   
 
Category:  None Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  jordigh
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 04 Apr 2013 01:31:39 PM UTC, comment #7: 

Ok, Dan, I've applied and merged your patch:

http://hg.savannah.gnu.org/hgweb/octave/graph/1766d8655006

Ed, can you use this as a starting point for fixing the test tolerances?

Thanks to you both for the work on this.

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Wed 03 Apr 2013 11:29:51 PM UTC, comment #6: 

Here is an update to the patch.  Same code, comments altered a tad, but mostly a much reduced changeset comment reformatted with line breaks.

I suggest applying this changeset, then closing this bug report.  After that, Ed can create a separate patch or bug report with tolerance changes based on norm(), starting with the diffs in my.patch.

(file #27766)

Dan Sebald <sebald>
Wed 03 Apr 2013 08:29:31 PM UTC, comment #5: 

Sorry that these patches fell through the cracks. Ed, I tried looking at yours, but it doesn't apply against the current sources anymore. Dan, your patch could work, but the commit message needs to be cleaned up. Overall, I do like the approach about fixing checking for real and imaginary parts.

Can you two agree on an updated patch, and I'll apply it?

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Fri 19 Oct 2012 04:03:26 AM UTC, comment #4: 

This patch is an update of octave-fftfilt-2012sep11.patch, having only a minor change of fix for computing the index for columns to be rounded.  The index should be an integer array, not logical array.  Note that this patch was generated using "hg export" and has the prescribed documentation header.

I didn't include eem2314's patch, so that will have to be re-added.  Let's try this a bit differently than using patch/modify/diff, and see if we can create a series of commits in a file as though it were done using push/pull.  That is,

1) use "hg import" to apply the latest patch
2) make the modifications for the new tolerance settings
3) commit those changes with the proper header documentation
4) "hg export" one previous to the tip (look at the "hg log")
5) attach the exported diff file here

This way, we'll have a series of diffs in just one file, easily applied by maintainers and giving credit of mods where appropriate.  Rather than redoing the patch for any tweaks, we'll just add more change sets.  We'll see how that works.

As for the problems that someone on the help-octave list is seeing for a compilation that might have a different math library, I just don't see how that can affect the rounding behavior and cause a test failure.  The result should be exact integer values for the first set of tests, and there isn't much that affects this:

  if (! any (b - fix (b)))
    idx = find (! any (x - fix (x)));
    y (:, idx) = round (y (:, idx));
  endif

That's it.  Could it be that the fix() routine fails under some other library so that one of the conditional tests fails?  I don't know.  I think the person who is seeing the problem will have to simply cut and past the test script and walk through the test looking at intermediate results.

(file #26788)

Dan Sebald <sebald>
Wed 10 Oct 2012 12:14:47 AM UTC, comment #3: 

This problem is identical to #35959
I believe this is another case where we need to account for the
size of the numbers being tested (see #34461). I sent a note to
the maintainers list with some tests I ran - this problem apparently arises when octave is built without fftw3 and fftpack is used instead.
I'm posting a modified version of Dan's patch from 11 Sep; either will fix the problem but this one uses norms instead of large multiples of eps.

(file #26729)

Edward Meyer <eem2314>
Wed 12 Sep 2012 02:08:04 AM UTC, comment #2: 

Here is a revamped version of the patch that is likely more robust using the all(imag(x(:,i))==0) construct rather than isreal().

In the case of strictly imaginary numbers (which are stored as complex) I've chosen to maintain the signs of the real component when set to zero, i.e.:

    y (:,xisimag) = complex (real (y (:,xisimag)) * 0, imag (y (:,xisimag)));

because that is more consistent with current behavior generally in Octave.  If it is desired that -0 should not be used in the real field, then one could do:

    y (:,xisimag) = complex (0, imag (y (:,xisimag)));

which probably requires a bit less CPU.  However, there may be some final tweak for this patch yet to address that.  So hold off until a couple weeks after this date in case the handling of -0 is changed in Octave.

(file #26537)

Dan Sebald <sebald>
Fri 07 Sep 2012 10:50:12 PM UTC, comment #1: 

It makes sense that we could allow the zeroing of real and imaginary parts on a column by column basis.  Here is a patch that does that and includes an additional test.


(file #26503)

Dan Sebald <sebald>
Fri 07 Sep 2012 09:08:33 PM UTC, original submission:  

One of the tests in fftfilt.m failed on my system:

  *** test
 r = sqrt (1/2) * (1+i);
 b = b*r;
 assert (fftfilt (b, x  ), r*[1 1 0 0 0 0 0 0 0 0]  , eps);
 assert (fftfilt (b, r*x), r*r*[1 1 0 0 0 0 0 0 0 0], eps);
 assert (fftfilt (b, x.'), r*[1 1 0 0 0 0 0 0 0 0].', eps);
!!!!! test failed
assert (fftfilt (b, r x),r r * [1, 1, 0, 0, 0, 0, 0, 0, 0, 0],eps) expected
 Columns 1 through 3:
   0.00000 + 1.00000i   0.00000 + 1.00000i   0.00000 + 0.00000i
 Columns 4 through 6:
   0.00000 + 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i
 Columns 7 through 9:
   0.00000 + 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i
 Column 10:
   0.00000 + 0.00000i
but got
 Columns 1 through 3:
   0.00000 + 1.00000i   0.00000 + 1.00000i   0.00000 - 0.00000i
 Columns 4 through 6:
  -0.00000 + 0.00000i  -0.00000 + 0.00000i   0.00000 + 0.00000i
 Columns 7 through 9:
  -0.00000 + 0.00000i   0.00000 - 0.00000i  -0.00000 + 0.00000i
 Column 10:
  -0.00000 + 0.00000i
maximum absolute error 2.22478e-16 exceeds tolerance 2.22045e-16
shared variables   scalar structure containing the fields:
    b =
       1   1
    x =
       1   0   0   0   0   0   0   0   0   0
    r = [](0x0)

This brought to light the fact that another test is somewhat meaningless:

%!test
%! b  = rand (10, 1);
%! x  = rand (10, 1);
%! y0 = filter (b, 1, x);
%! y  = filter (b, 1, x);
%! assert (y, y0);

And then came forth a bug in the code itself.


I've attached a patch that

1) Relaxes the tolerance of the test that failed from eps to 2*eps.  (Explained later.)
2) Changes the meaningless test so that convolution and fftfilt are compared for both real and imaginary inputs, with appropriate tolerance.
3) In the final touchup of the results, handles the cases of b real and x imaginary, b imaginary and x real, and b and x imaginary.  (Note that the existing code for integerizing the output works appropriately for both x and b integer requirement, which is nice.)
4) Adds another test in which both x and b are complex.
5) Moves the transpose to the very end of the routine because rounding was failing.

Here is greater details on the tests and choice of tolerance:

%!test
%! b = [1 1];
%! x = [1, zeros(1,9)];
%! assert (fftfilt (b,  x  ), [1 1 0 0 0 0 0 0 0 0]  );
%! assert (fftfilt (b,  x.'), [1 1 0 0 0 0 0 0 0 0].');
%! assert (fftfilt (b.',x  ), [1 1 0 0 0 0 0 0 0 0]  );
%! assert (fftfilt (b.',x.'), [1 1 0 0 0 0 0 0 0 0].');

The tolerance of epsilon was removed.  Because the inputs are integers, the output is cast to integer and these results should be an exact comparison.  (I do wonder if there should be an option to remove the final touchups.)  NOTE: THIS CHANGE BROUGHT OUT A BUG IN rounding.  y was transposed before final touchup so y could be a ROW vector in which case:

    idx = !any (x - fix (x));

is a value of 1.  Thus for a row vector

    y(:, idx) = round (y(:, idx));

only was rounding the first element of the vector.

%! r = sqrt (1/2) * (1+i);
%! b = b*r;
%! assert (fftfilt (b, x  ), r*[1 1 0 0 0 0 0 0 0 0]  , eps  );

This test passes because although b is complex, x is not complex.  So given that x is a unit pulse (all imaginary components are 0), the mixing is such that there is really only one multiply summed throughout the whole transform.

%! assert (fftfilt (b, r*x), r*r*[1 1 0 0 0 0 0 0 0 0], 2*eps);

This is the test that failed.  Notice it is similar to the previous, but now both b and r*x are complex.  No longer is the imaginary component all zeros, with the implication that butterflies contribute "on both wings".  Hence, with the unit pulse, but complex operations, I'm thinking that works out to two multiplications within summations, hence 2*eps can be lost.

%! assert (fftfilt (b, x.'), r*[1 1 0 0 0 0 0 0 0 0].', eps  );

Similar to a couple tests previous, just different orientation.

%! b  = [1 1];
%! x  = zeros (10,3); x(1,1)=-1; x(1,2)=1;
%! y0 = zeros (10,3); y0(1:2,1)=-1; y0(1:2,2)=1;
%! y  = fftfilt (b, x);
%! assert (y0, y);
%! y  = fftfilt (b*i, x);
%! assert (y0*i, y);
%! y  = fftfilt (b, x*i);
%! assert (y0*i, y);
%! y  = fftfilt (b*i, x*i);
%! assert (-y0, y);

I simply added all the combinations of real and imaginary inputs to test the enhancements to the final touchups.

%!test
%! b  = rand (10, 1);
%! x  = rand (10, 1);
%! y0 = filter (b, 1, x);
%! y  = fftfilt (b, x);
%! assert (y0, y, 15*eps);
%! y0 = filter (b*i, 1, x*i);
%! y  = fftfilt (b*i, x*i);
%! assert (y0, y, 20*eps);

This is a test that really checks the integrity of the algorithm.  (The meaningless test prior to the patch.)  After some initial guesses and trial and error, I tested 1e6 trials, relaxing tolerance so that they all passed.  A tolerance of 15*eps worked for real inputs, but 15*eps failed long into tests for the imaginary inputs.  I tried that test with 15*eps tolerance again and all 1e6 computations were within tolerance.  So, I bumped that test to 16*eps for both cases.

%!test
%! b  = rand (10, 1) + i*rand (10, 1);
%! x  = rand (10, 1) + i*rand (10, 1);
%! y0 = filter (b, 1, x);
%! y  = fftfilt (b, x);
%! assert (y0, y, 55*eps);

OK, this test pushes the limit: both b and x are complex, non-integer so there is no zeroing out either the real or imaginary component.  The 50*eps tolerance did not pass 1e6 trials, so I changed that to 55*eps and that passed.

Dan Sebald <sebald>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #27766:  octave-fftfilt-2013apr03.patch added by sebald (4KiB - application/octet-stream)
file #26729:  my.patch added by eem2314 (4KiB - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by eem2314 (Updated the item)
  • -email is unavailable- added by sebald (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 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2013-04-04 jordigh StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2013-04-03 sebald Attached File- Added octave-fftfilt-2013apr03.patch, #27766
    2013-04-03 jordigh StatusNone In Progress
        Assigned toNone jordigh
    2012-10-19 sebald Attached File- Added octave-fftfilt-2012oct18.patch, #26788
    2012-10-10 eem2314 Attached File- Added my.patch, #26729
    2012-09-12 sebald Attached File- Added octave-fftfilt-2012sep11.patch, #26537
    2012-09-07 sebald Attached File- Added octave-fftfilt-2012sep07_6pm.patch, #26503
    2012-09-07 sebald Attached File- Added octave-fftfilt-2012sep07.patch, #26501

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code