Opened 9 years ago

Closed 9 years ago

#30 closed defect (fixed)

Results under OS X & Linux differ from "known good"

Reported by: flip Owned by: flip
Priority: major Milestone:
Component: hlsvd Version:
Keywords: Cc: bsoher

Description (last modified by flip)

Some history is warranted. HLSVD's origins in this project go back to the IDL version of Vespa. Brian and/or I cajoled two sets of output from IDL Vespa's HLSVD function and matched them up with their respective inputs. They became cp0.xml and mystery.xml in the vespa/hlsvd/demo directory. We assumed that the HLSVD executed via IDL Vespa gave correct results.

Fast forward to today where we have re-run the original Fortran code through f2c, recompiled it with modern compilers, and added a Python wrapper. Running the demo code under OS X and Windows gives similar results (to one another), with one exception. Using mystery.xml as input, OS X gets back one fewer singular values (25 instead of the expected 26) and a final iteration count of 500 versus the expected 550.

Although the results agree (mostly) under OS X & Windows, they diverge from the "known good" results produced by IDL Vespa. In some cases the difference is many orders of magnitude. (A random example -- in cp0, amplitude #6 is expected to be 670.377, the actual value is 0.272528.)

Furthermore, there was a bug in the IDL Vespa code that wrapped the HLSVD DLL. HLSVD's C function signature looked like this:

int hlsvd_( real        *yr,
            real        *yi,
            integer     *ndata,
            integer     *ndp,
            real        *tstep,
            integer     *m,
            integer     *nssv,
            doublereal  *rootr,
            doublereal  *rooti,
            real        *fre,
            real        *dam,
            real        *amp,
            real        *phase,
            doublereal  *sv,
            doublereal  *snois,
            doublereal  *rms,
            integer     *ndiv,
            integer     *nit,
            integer     *nitf)

But in the IDL code, sv, snoise and rms were all declared as fltarr or float when they should be doubles. There was not enough space allocated for the C DLL to write to. In practice, I saw snoise and rms get trashed by this bug, and I'm not sure if our "known good" values were affected by it.

Even taking into account that bug, it doesn't explain why the results from our modern version of HLSVD diverge so wildly from the old version. One (or both) must be wrong.

Change History (4)

comment:1 Changed 9 years ago by flip

  • Description modified (diff)

Some history is warranted. HLSVD's origins in this project go back to the IDL version of Vespa. Brian and/or I cajoled two sets of output from IDL Vespa's HLSVD function and matched them up with their respective inputs. They became cp0.xml and mystery.xml in the vespa/hlsvd/demo directory. We assumed that the HLSVD executed via IDL Vespa gave correct results.

Fast forward to today. We have re-run the original Fortran code through f2c, recompiled it with modern compilers, and added a Python wrapper. Running the demo code under OS X and Windows gives results similar to one another, with one exception. Using mystery.xml as input, OS X gets back one fewer singular values (25 instead of the expected 26) and a final iteration count of 500 versus the expected 550.

Although the results agree (mostly) under OS X & Windows, they diverge from the "known good" results produced by IDL Vespa. In some cases the difference is many orders of magnitude. (A random example -- in cp0, amplitude #6 is expected to be 670.377, the actual value is 0.272528.)

Furthermore, there was a bug in the IDL Vespa code that wrapped the HLSVD DLL. HLSVD's C function signature looks like this:

int hlsvd_( real        *yr,
            real        *yi,
            integer     *ndata,
            integer     *ndp,
            real        *tstep,
            integer     *m,
            integer     *nssv,
            doublereal  *rootr,
            doublereal  *rooti,
            real        *fre,
            real        *dam,
            real        *amp,
            real        *phase,
            doublereal  *sv,
            doublereal  *snois,
            doublereal  *rms,
            integer     *ndiv,
            integer     *nit,
            integer     *nitf)

But in the IDL code, sv, snoise and rms were all declared as fltarr or float when they should have been doubles. There was not enough space allocated for the C DLL to write to. In practice, I saw snoise and rms get trashed by this bug, and I'm not sure if our "known good" values were affected by it.

Even taking into account that bug, it doesn't explain why the results from our modern version of HLSVD diverge so wildly from the old version. One (or both) must be wrong.

comment:2 Changed 9 years ago by flip

Brian and I have discussed this, and what we (mostly Brian) concluded is
that since HLSVD present as simplified picture of data, that simplification
can happen differently on different platforms depending on minor math library
differences, floating point noise, etc.

In short, we can expect answers in a certain range, but there's no single
correct answer.

The "known good" results in cp0.xml are almost certainly wrong, since no
version of HLSVD that we have (f2c or pure Fortran, 32- or 64-bit, on three
different platforms) can reproduce the results in that file. What to replace
them with I'm not sure since the whole idea of the demo as it is currently
constructed requires known good values.

I would discard the entire demo except that it also makes a handy unit test.
Probably the best thing to do is to convert the demo into a proper test
and include a note to the effect that the results should be close but
might not necessarily match perfectly.

comment:3 Changed 9 years ago by flip

  • Owner set to flip
  • Status changed from new to assigned

comment:4 Changed 9 years ago by flip

  • Resolution set to fixed
  • Status changed from assigned to closed

Fixed in r2669 by changing the demo into a test and better advertising its limitations.

Note: See TracTickets for help on using tickets.