gsl_lstsq.c

Go to the documentation of this file.

/*============================================================================*/
#include "gsl_fit.h"

int lstsq_main(void *void_gsl_fit_struct)
{
  unsigned int                        count, i;
  int                                 stop1, stop2;
  gsl_vector                         *guess;
  gsl_vector                         *gradt;
  gsl_matrix                         *covar;
  gsl_multifit_function_fdf           lstsq_func;
  gsl_multifit_fdfsolver             *solvr;
  gsl_fit                            *gft;
  
  gft = (gsl_fit *)void_gsl_fit_struct;
  
  /* Inital guess is the paramters as they were passed into here */
  guess = gsl_vector_alloc(gft->para.num);
  for (i=0; i<gft->para.num; i++) 
  {
    gsl_vector_set(guess, i, gft->para.par[i]);
  }

  /* count the number of data points */
  count = 0;
  for (i=0; i<gft->numd; i++) count += gft->data[i].num;

  /* need this vector for stop test #2 */
  gradt = gsl_vector_alloc(gft->para.num);
  
  /* covariant matrix  */
  covar = gsl_matrix_alloc(gft->para.num,gft->para.num);

  lstsq_func.f      = &lstsq_f;
  lstsq_func.df     = &lstsq_df;
  lstsq_func.fdf    = &lstsq_fdf;
  lstsq_func.n      = count;     /* number of data points */
  lstsq_func.p      = gft->para.num; /* number of paramters   */
  lstsq_func.params = (void *)gft; /* pass the full gsl_fit struct */
  
  solvr = gsl_multifit_fdfsolver_alloc(gft->meth.lstsq_type, \
                                       lstsq_func.n, lstsq_func.p);

  /*==================================================*/
  count = 0;
  gsl_multifit_fdfsolver_set(solvr, &lstsq_func, guess);
  gsl_multifit_covar(solvr->J, 0.0, covar);
  lstsq_status(gft, count, solvr, covar);

  do
  { 
    count++;
    stop1 = gsl_multifit_fdfsolver_iterate(solvr);    
    gsl_multifit_covar(solvr->J, 0.0, covar);
    lstsq_status(gft, count, solvr, covar);
    
    if (stop1) 
    {
      fprintf(stderr,"STATUS: %s\n", gsl_strerror(stop1));
      break;
    }
    
    /* Stop test #1 */
    stop1 = gsl_multifit_test_delta(solvr->dx, solvr->x, 1e-6, 1e-6);

    /* Stop test #2 */    
    gsl_multifit_gradient(solvr->J, solvr->f, gradt);
    stop2 = gsl_multifit_test_gradient(gradt, 1e-6);
    
    if ((stop1==GSL_SUCCESS) && (stop2==GSL_SUCCESS))
    {
      fprintf(stderr,"STATUS: %s\n", gsl_strerror(stop1));
    }

  } while (((stop1 == GSL_CONTINUE) || (stop2 == GSL_CONTINUE)) \
            && (count < gft->maxi));

  /*==================================================*/
  gsl_multifit_covar(solvr->J, 0.0, covar);
  lstsq_status(gft, -1, solvr, covar);
  
  /*==================================================*/

  gsl_multifit_fdfsolver_free(solvr);
  gsl_matrix_free(covar);
  gsl_vector_free(guess);
  gsl_vector_free(gradt);

  return GSL_SUCCESS;
}
/*===================================================================
 *    lstsq_fdf 
 */
int lstsq_fdf(const gsl_vector *parameters, \
              void *void_gsl_fit_struct, \
              gsl_vector *r, gsl_matrix *J)
{
  lstsq_f (parameters, void_gsl_fit_struct, r);
  lstsq_df(parameters, void_gsl_fit_struct, J);
  return GSL_SUCCESS;
}
/*===================================================================
 *    lstsq_f 
 */
int lstsq_f(const gsl_vector *parameters, \
            void *void_gsl_fit_struct, \
            gsl_vector *r)
{
  int                i, j, k;
  int                r_index;
  double             tmp;
  gsl_fit *gft;
  
  gft = (gsl_fit *)void_gsl_fit_struct;
  /*
  **  Load the params back into the parameter array.
  **  Why should we do this step?
  **  We could say below that p = params.data + d->func[j].idx
  */
  for(i=0; i<gft->para.num; i++)
  {
    gft->para.par[i] = gsl_vector_get(parameters, i);
  }
  /*
  ** Loop over the data sets with index k.
  */
  r_index = 0;
  for (k=0; k<gft->numd; k++)
  {
    /*
    ** Loop over the data points with index i
    */ 
    for (i=0; i<gft->data[k].num; i++)
    {
      /*
      **   Loop over the functions with index j
      */
      tmp = 0.0;
      for (j=0; j<gft->numf; j++) 
      {
        if (gft->type==MULTIPLE)
        {
          tmp += gft->func[j].fun_m(gft->data[k].xpt[i], \
                                    gft->para.par + gft->func[j].idx, \
                                    k, j);
        }
        else
        {
          tmp += gft->func[j].fun_s(gft->data[k].xpt[i], \
                                    gft->para.par + gft->func[j].idx);
        }
      }
      tmp = (tmp - gft->data[k].ypt[i]) / gft->data[k].sig[i];    
      gsl_vector_set(r, r_index, tmp);
      r_index++; 
    }
  } 
  return GSL_SUCCESS;
}
/*===================================================================
 *    lstsq_df 
 */
int lstsq_df(const gsl_vector *parameters, \
             void *void_gsl_fit_struct, \
             gsl_matrix *J)
{
  int                i, j, k;
  int                r_index;
  double             res, err;
  gsl_function       F;
  gsl_fit *gft;
  
  gft = (gsl_fit *)void_gsl_fit_struct;
  /*
  **  Load the params back into the parameter array.
  **  Why should we do this step?
  **  We could say below that p = params.data + d->func[j].idx
  */
  for(i=0; i<gft->para.num; i++)
  {
    gft->para.par[i] = gsl_vector_get(parameters,i);
  }

  F.params   = void_gsl_fit_struct;
  F.function = &lstsq_derivatives;

  /*
  ** Loop over the data sets with index k.
  */
  r_index = 0;
  for (k=0; k<gft->numd; k++)
  {
    gft->indx = k;
    /*
    **  The jacobian mtirx J_ij = d(r_i)/d(a_j)
    **  => i is the data index of residual r
    **  => j is the parameter index
    */
    /*
    ** Loop over the data points with index i
    */ 
    for (i=0; i<gft->data[k].num; i++)
    {
      gft->data[k].idx = i;
      /*
      **   Loop over the parameters with index j
      */

      for (j=0; j<gft->para.num; j++) 
      {
        gft->para.idx = j;
        /*
        ** If fixed, J_ij = 0
        */
        if (gft->para.fix[j]) res=0.0;
        else gsl_deriv_central(&F, gft->para.par[j], 1e-4, &res, &err);
        gsl_matrix_set(J, r_index, j, res);
      }
      r_index += 1;
    }
  } 
    
  return GSL_SUCCESS;
}
/*===================================================================
 *    lstsq_derivatives 
 */
double lstsq_derivatives(double p, void *void_gsl_fit_struct)
{
  int                i, j;
  double             tmp;
  gsl_fit *gft;
  
  gft = (gsl_fit *)void_gsl_fit_struct; 
  /*
  **  We need a complete set of paramters 'tmp',
  **  with the derivative parameter 'p' 
  **  replacing the parameter at index 'd->para.idx'.
  */
  for(i=0; i<gft->para.num; i++) gft->para.tmp[i] = gft->para.par[i];
  gft->para.tmp[gft->para.idx] = p;
  
  tmp = 0.0;
  i  = gft->data[gft->indx].idx;
  for (j=0; j<gft->numf; j++) 
  {
    if (gft->type == MULTIPLE)
    {
      tmp += gft->func[j].fun_m(gft->data[gft->indx].xpt[i], \
                                gft->para.tmp+gft->func[j].idx, \
                                gft->indx, j);
    }
    else
    {
      tmp += gft->func[j].fun_s(gft->data[gft->indx].xpt[i], \
                                gft->para.tmp+gft->func[j].idx);
    }
  }
  tmp = (tmp - gft->data[gft->indx].ypt[i]) / gft->data[gft->indx].sig[i];    
  return tmp;
}
/*===================================================================
 *    lstsq_status 
 */
int lstsq_status(gsl_fit *gft, \
                 int iteration, \
                 gsl_multifit_fdfsolver *s, \
                 gsl_matrix *c)
{
  int      j;
  double   chisq;
  
  if (iteration >= 0)
  {
    fprintf(stderr,"ITERATION:%3u ", iteration);
    fprintf(stderr,"===================\n");
  }
  else
  {
    fprintf(stderr,"COMPLETE ");
    fprintf(stderr,"===================\n");
  }
  chisq = gsl_blas_dnrm2(s->f)/(double)((gft->numd)-(gft->para.num));
  fprintf(stderr,"ChiSq: %12.6e\n", chisq);

  for (j=0; j<gft->para.num; j++)
  {
    fprintf(stderr,"%3d: %12.6f +/- %8.6f\n",
            j,
            gsl_vector_get(s->x, j),
            sqrt(gsl_matrix_get(c, j, j)));
  }
  fprintf(stderr,"\n");
  
  if (iteration < 0)
  {
    for (j=0; j<gft->para.num; j++)
    {
      gft->para.par[j] = gsl_vector_get(s->x, j);
      gft->para.err[j] = sqrt(gsl_matrix_get(c, j, j));
    }
    gft->chis = chisq;
  }
  
  return GSL_SUCCESS;
}

---