[Help-gsl] Non-linear least squares fits.

[ja]

Hello,
I have a problem where I'm trying to fit five parameters using a
non-linear least squares fit to two spectral lines simultaneously. One
spectral line is held in one array and one in another. The code is written
in C++ using static class members. I'm assuming at the moment that there's
a problem in the differentials as the code compiles but the solver exits
immediately, with no change in the parameters. I'm trying to check the
values of the differentials using mathematica.

I would be really grateful if anyone could check if my code is roughly
correct given what I'm trying to achieve.

Thanks in advance!

 Following the example given in the documentation, my structure looks like
this;

struct data{
  size_t n;
  double * xl;
  double * yl;
  double *xu;
  double *yu;
  double sigma1;
  double sigma2;
  double * param;
};

My actual fitting code is as follows;
int nh3fitting::fit()
{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *s;

  int status;
  size_t i, iter=0;
 std::cout << n << std::endl;
  const size_t ni =n;
  const size_t p =np;

  gsl_matrix *covar = gsl_matrix_alloc(p,p);
  struct data d ={ni,xlower,ylower,xupper,yupper,noise1,noise2,param};
//be careful these need to be in the order specified in the struct
  gsl_multifit_function_fdf f;
  double x_init[np]={param[0],param[1],param[2],param[3],param[4]};
//initial guesses
  gsl_vector_view x = gsl_vector_view_array(x_init,p);

  f.f = &expb_f; //loads f with functions from expfit
  f.df = &expb_df;
  f.fdf = &expb_fdf;
  f.n =2*ni; //number of data points
  f.p=p; // number of parameters to fit.
  f.params=&d;  //loads data

  T=gsl_multifit_fdfsolver_lmsder;
  s=gsl_multifit_fdfsolver_alloc(T,2*ni,p);
  std::cout << "here " << std::endl;
  gsl_multifit_fdfsolver_set(s,&f,&x.vector);

  do
    {
      iter++;
      status = gsl_multifit_fdfsolver_iterate(s);

      if(status)
        break;

      status = gsl_multifit_test_delta(s->dx,s->x,1e-10,1e-10);
    }
  while (status == GSL_CONTINUE && iter <500);

  gsl_multifit_covar(s->J,0.0, covar);

#define FIT(i) gsl_vector_get(s->x,i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))

 {
   for(i=0;i<np;i++)
     {
       param[i]=FIT(i);
       error[i]=ERR(i);
     }
 }
 gsl_multifit_fdfsolver_free (s);
 std::cout << "Parameters are (fit) " << param[0] << " " << param[1] << "
"  << param[2] << " " << param[3] << "  " << param[4] << std::endl;
 return 0;
}

The helper functions providing the function and the differential are
below.calculatep* gives the differential and calcTaFull gives the values
of the functions at that value.

int nh3fitting::expb_f(const gsl_vector * x, void *data,  gsl_vector *f )
{
        std::cout << "starting expb_f " << std::endl;

        size_t n=((struct data *)data)->n;
        double *xl =((struct data *)data)->xl;
        double *xu =((struct data *)data)->xu;
        double *yl =((struct data *)data)->yl;
        double *yu =((struct data *)data)->yu;
        double sigma1 =((struct data *)data)->sigma1;
        double sigma2 =((struct data *)data)->sigma2;
        double *param=((struct data *)data)->param;

        double tr,tau2,vlsrfreqlower,vlsrfrequpper,lwfreqlower,lwfrequpper;
        nh3fitting::initialise(
param,tr,tau2,vlsrfreqlower,vlsrfrequpper,lwfreqlower,lwfrequpper);

        size_t i;

        for(i=0;i<n;i++)
    {
        double Yi =calcTaFull(xl[i], param, tau2,vlsrfreqlower,
vlsrfrequpper, lwfreqlower, lwfrequpper);
                gsl_vector_set(f,i,(Yi-yl[i])/sigma1);
    }

    for(i=0;i<n;i++)
    {
        double Yi =calcTaFull(xu[i], param, tau2,vlsrfreqlower,
vlsrfrequpper, lwfreqlower, lwfrequpper);
                gsl_vector_set(f,i+n,(Yi-yu[i])/sigma2);
    }
        return GSL_SUCCESS;
}


int nh3fitting::expb_df(const gsl_vector * x, void *data, gsl_matrix * J)
{
        size_t n=((struct data *)data)->n;
        double *xl =((struct data *)data)->xl;
        double *xu =((struct data *)data)->xu;
        double s1 =((struct data *)data)->sigma1;
        double s2 =((struct data *)data)->sigma2;
        double *param=((struct data *)data)->param;

        size_t i;

        double
tr,tau2,vlsrfreqlower,vlsrfrequpper,lwfreqlower,lwfrequpper,tau,dtdtau;
        
nh3fitting::initialise(param,tr,tau2,vlsrfreqlower,vlsrfrequpper,lwfreqlower,lwfrequpper);

  for(i=0;i<n;i++)
    {
        tau=calcTau(xl[i],param[4],  tau2, vlsrfreqlower,  vlsrfrequpper,
lwfreqlower, lwfrequpper);
        dtdtau=calculatedtadtau(xl[i], param,tau);
        gsl_matrix_set(J,i,0,calculatedp0(xl[i],param,tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper, lwfreqlower)/s1);
        gsl_matrix_set(J,i,1,calculatedp1(xl[i], param, tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper, lwfreqlower)/s1);
        gsl_matrix_set(J,i,2,calculatedp2(xl[i], param, tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper, lwfreqlower)/s1);
        gsl_matrix_set(J,i,3,calculatedp3(xl[i], param, tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper, lwfreqlower)/s1);
        gsl_matrix_set(J,i,4,calculatedp4(xl[i], param, tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper, lwfreqlower)/s1);
    }

     for(i=0;i<n;i++)
    {
        tau=calcTau(xu[i],param[4],  tau2, vlsrfreqlower,  vlsrfrequpper,
lwfreqlower, lwfrequpper);
        dtdtau=calculatedtadtau(xu[i], param,tau);
        gsl_matrix_set(J,i+n,0,calculatedp0(xu[i], param, tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper, lwfreqlower)/s2);
        gsl_matrix_set(J,i+n,1,calculatedp1(xu[i], param, tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper, lwfreqlower)/s2);//this all
needs to  be i+1;
        gsl_matrix_set(J,i+n,2,calculatedp2(xu[i], param, tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper,  lwfreqlower)/s2);
        gsl_matrix_set(J,i+n,3,calculatedp3(xu[i], param, tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper,  lwfreqlower)/s2);
        gsl_matrix_set(J,i+n,4,calculatedp4(xu[i], param, tau,dtdtau,tau2,
vlsrfreqlower, vlsrfrequpper, lwfrequpper,  lwfreqlower)/s2);
    }

        return GSL_SUCCESS;
}