[Bug-gsl] Bug in lm(s)der_set - cannot reuse fitter

[Christian Holm Christensen]

Hi there,

There's a bug in the lm(s)der_set functions (called by
gsl_multifit_fdfsolver_set), that makes it impossible to reuse a
fitter. 

Suppose you fit one data set by doing 

        struct data { ... }
        int func(const gsl_vector* x, void* data, gsl_vector* f, gsl_matrix* j)
        {
          ...
        }
        
        int main()
        {
          const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
          gsl_multifit_fdfsolver            *s = gsl_matrix_alloc (p, p);
          gsl_multifit_function_fdf          f;
          f.fdf = &func;
          ....
          make_data();
          ...
          gsl_multifit_fdfsolver_set (s, &f, &x.vector);
          do
          {
            status = gsl_multifit_fdfsolver_iterate (s);
            if (status) break;
            status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
          } while (status == GSL_CONTINUE && iter < 500);
          ...
          make_data();
          gsl_multifit_fdfsolver_set (s, &f, &x.vector);
          do
          {
            status = gsl_multifit_fdfsolver_iterate (s);
            if (status) break;
            status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
          } while (status == GSL_CONTINUE && iter < 500); 
          ...
        }

Then, the second fit may fail.  The reason is, that the internal
structures of the lm(s)der algorithm are not reset properly.  

To see this with real data, please download "rcugsl" (C++ wrappers of
GSL histograms & fitters - not big) and some data.  Both can be found
at 

        http://fmd.nbi.dk/fmd/fee/rcugsl/
        Direct links:
                http://fmd.nbi.dk/fmd/fee/rcugsl/rcugsl-0.2.tar.gz
                http://fmd.nbi.dk/fmd/fee/rcugsl/data.tar.gz 
        
and then run the program "testBug" with option '-e' (enable bug fix) and
option '-d' (disable bug fix) on the extracted data 

        ~/rcugsl-0.2/rcugsl/testBug [-e|-d] *.hist
        
The bug fix is in `NonLinearFitter::Set(Function& f, DataSet& d, bool
bugfix)'. 

The data is "real" pedestal data from a silicon sensor attached to a
10bit ADC.   As you can see, the pedestals are rather narrow - however,
that is not the problem (I filter out spectra with an sigma less than
0.01 ADC counts). 

I'm using GSL 1.9-3 from Debian GNU/Linux unstable (sid) on an i386.
I'm using GCC 4.1.3 (prerelease). 

Thank you for your attention and a great library.

Yours,

-- 
 ___  |  Christian Holm Christensen 
  |_| |  -------------------------------------------------------------
    | |  Address: Sankt Hansgade 23, 1. th.  Phone:  (+45) 35 35 96 91
     _|           DK-2200 Copenhagen N       Cell:   (+45) 24 61 85 91
    _|            Denmark                    Office: (+45) 353  25 404
 ____|   Email:   address@hidden               Web:    www.nbi.dk/~cholm
 | |