gsl_fit.c

Go to the documentation of this file.

/*============================================================================*/
/*==========================================================================*/
#include "gsl_fit.h"
/*==========================================================================*/
gsl_fit *gsl_fit_alloc(unsigned int method, \
                       unsigned int algotrithm, \
                       unsigned int max_iter, \
                       unsigned int type)
{
  /* 
   * This instance will be created in memory and the 
   * pointer to it passed back to the user 
   */
  gsl_fit *gft;
  
  /* Sanity checks */
  if ((method<MULTIDIM_MINIMIZAION)||(method>SIMULATED_ANNEAL)) {
    GSL_ERROR_VAL("Not a valid method", \
                   GSL_EINVAL, NULL); }
  if ((algotrithm<CONJUGATE_FR)||(algotrithm>LMDR)) {
    GSL_ERROR_VAL("Not a valid algortithm", \
                   GSL_EINVAL, NULL); }
  if ((method == MULTIDIM_MINIMIZAION) && (algotrithm > STEEPEST_DESCENT)) {
    GSL_ERROR_VAL("Wrong method/algotrithm pair.", \
                   GSL_EINVAL, NULL); }
  else if (method == LEAST_SQUARES)  {
    if ((algotrithm <= STEEPEST_DESCENT) || (algotrithm > LMDR)) {
      GSL_ERROR_VAL("Wrong method/algotrithm pair.", \
                     GSL_EINVAL, NULL); } }
  if ((method!=SIMULATED_ANNEAL)&&(max_iter<=0)) {
    GSL_ERROR_VAL("Not enough iterations", \
               GSL_EINVAL, NULL); }
  if ((type!=SINGLE)&&(type!=MULTIPLE)) {
    GSL_ERROR_VAL("Not valid type", \
               GSL_EINVAL, NULL); }
  /* Allocate memory for the gsl_fit */
  gft = (gsl_fit *)malloc(sizeof(gsl_fit));
  if (gft==NULL) {
    GSL_ERROR_VAL("Failed to allocate a gsl_fit", \
                   GSL_ENOMEM, NULL); }
  
  /* Initialize the data struct */
  gft->numd = 0;
  gft->data = NULL;
  
  /* Initialize stuff */
  gft->chis = 0.0;  
  gft->maxi = max_iter;
  gft->indx = 0;
  gft->type = type;
  
  /* Initialize the paramter struct */
  gft->para.num = 0;
  gft->para.idx = 0;
  gft->para.par = NULL;
  gft->para.fix = NULL;
  gft->para.err = NULL;
  gft->para.tmp = NULL;

  /* Initialize the function struct */
  gft->numf = 0;
  gft->func = NULL;

  /*------------------------------------------------------*/
  /* Set the fitting method */
  switch (method) {
    case LEAST_SQUARES:
      gft->meth.solver_function_call = &lstsq_main;
      break;
    case MULTIDIM_MINIMIZAION: 
      gft->meth.solver_function_call = &mdmin_main;
      break;
    case SIMULATED_ANNEAL: 
      gft->meth.solver_function_call = &siman_main;
      break;
    default: 
      gft->meth.solver_function_call = &lstsq_main;
  }

  switch (algotrithm) {
    case CONJUGATE_FR: 
      gft->meth.mdmin_type = gsl_multimin_fdfminimizer_conjugate_fr;
      break;
    case CONJUGATE_PR: 
      gft->meth.mdmin_type = gsl_multimin_fdfminimizer_conjugate_pr;
      break;
    case VECTOR_BFGS: 
      gft->meth.mdmin_type = gsl_multimin_fdfminimizer_vector_bfgs;
      break;
    case STEEPEST_DESCENT: 
      gft->meth.mdmin_type = \
                     gsl_multimin_fdfminimizer_steepest_descent;
      break;
    case LMSDR: 
      gft->meth.lstsq_type = gsl_multifit_fdfsolver_lmsder;
      break;
    case LMDR: 
      gft->meth.lstsq_type = gsl_multifit_fdfsolver_lmder;
      break;
    default:
      gft->meth.lstsq_type = gsl_multifit_fdfsolver_lmder;
      break;
  }
  
  return gft;
}
/*==========================================================================*/
int gsl_fit_add_dataset(gsl_fit *gft, \
                        double *xpt, double *ypt, double *sig, \
                        unsigned int num)
{
  int i;
  
  /* Sanity checks */
  if (gft==NULL) {
    GSL_ERROR("Attempt to add data to null gsl_fit", \
               GSL_EFAULT); }
  if ((gft->type==SINGLE)&&(gft->numd!=0)) {
    GSL_ERROR("Already have data to fit", \
               GSL_EINVAL); }
  if (num<=2) {
    GSL_ERROR("Not enough data points", \
               GSL_EINVAL); }
               
  /*==============*/
  /* allocate the memory for the */
  gft->data = realloc(gft->data, (gft->numd+1)*sizeof(struct _data));
  if (gft->data == NULL) {
      GSL_ERROR("No memory for data set",\
                 GSL_ENOMEM); }
    
  /*==============*/
  gft->data[gft->numd].xpt = (double *)malloc(num*sizeof(double));
  gft->data[gft->numd].ypt = (double *)malloc(num*sizeof(double));
  gft->data[gft->numd].sig = (double *)malloc(num*sizeof(double));
  if ((gft->data[gft->numd].xpt == NULL) || \
      (gft->data[gft->numd].ypt == NULL) || \
      (gft->data[gft->numd].sig == NULL)) {
      GSL_ERROR("No memory for data set",\
                 GSL_ENOMEM); }
  
  /*==============*/
  gft->data[gft->numd].num = num;
  gft->data[gft->numd].idx = 0;
  
  /*==============*/
  for (i=0; i<num; i++)
  {
    gft->data[gft->numd].xpt[i] = xpt[i];
    gft->data[gft->numd].ypt[i] = ypt[i];
    gft->data[gft->numd].sig[i] = sig[i];
  }
  
  /*==============*/
  gft->numd += 1;

  /*==============*/
  return GSL_SUCCESS;
}
/*==========================================================================*/
int gsl_fit_add_model_multiple(gsl_fit *gft, \
                               double (*fun)(double, double *, \
                                             unsigned int, unsigned int), \
                               double *par, unsigned int *fix, unsigned int num)
{
  unsigned int i;
  
  /* Sanity checks */
  if (gft==NULL) {
    GSL_ERROR("Attempt to add data to null gsl_fit", \
               GSL_EFAULT); }
  if (num<=0) {
    GSL_ERROR("Not enough data points", \
               GSL_EINVAL); }
  if (gft->type==SINGLE) {
    GSL_ERROR("Wrong model type", \
               GSL_EINVAL); } 
                             
  gft->func     = realloc(gft->func,     (gft->numf+1)*sizeof(struct _func));
  gft->para.par = realloc(gft->para.par, (gft->para.num+num)*sizeof(double));
  gft->para.fix = realloc(gft->para.fix, \
                                   (gft->para.num+num)*sizeof(unsigned int));
  gft->para.err = realloc(gft->para.err, (gft->para.num+num)*sizeof(double));
  gft->para.tmp = realloc(gft->para.tmp, (gft->para.num+num)*sizeof(double));
  
  gft->func[gft->numf].fun_m = fun;
  gft->func[gft->numf].fun_s = NULL;
  gft->func[gft->numf].num = num;
  gft->func[gft->numf].idx = gft->para.num;
  
  for (i=0; i<num; i++)
  {
    gft->para.par[i+gft->para.num] = par[i];
    gft->para.fix[i+gft->para.num] = fix[i];
    gft->para.err[i+gft->para.num] = 0.0;
    gft->para.tmp[i+gft->para.num] = 0.0;
  } 

  gft->numf += 1;
  gft->para.num += num;
  
  return GSL_SUCCESS;
}
/*==========================================================================*/
int gsl_fit_add_model_single(gsl_fit *gft, \
                             double (*fun)(double, double *), \
                             double *par, unsigned int *fix, unsigned int num)
{
  unsigned int i;
  
  /* Sanity checks */
  if (gft==NULL) {
    GSL_ERROR("Attempt to add data to null gsl_fit", \
               GSL_EFAULT); }
  if (num<=0) {
    GSL_ERROR("Not enough data points", \
               GSL_EINVAL); }
  if (gft->type==MULTIPLE) {
    GSL_ERROR("Wrong model type", \
               GSL_EINVAL); }       
       
  gft->func     = realloc(gft->func,     (gft->numf+1)*sizeof(struct _func));
  gft->para.par = realloc(gft->para.par, (gft->para.num+num)*sizeof(double));
  gft->para.fix = realloc(gft->para.fix, \
                                   (gft->para.num+num)*sizeof(unsigned int));
  gft->para.err = realloc(gft->para.err, (gft->para.num+num)*sizeof(double));
  gft->para.tmp = realloc(gft->para.tmp, (gft->para.num+num)*sizeof(double));
  
  gft->func[gft->numf].fun_m = NULL;
  gft->func[gft->numf].fun_s = fun;
  gft->func[gft->numf].num = num;
  gft->func[gft->numf].idx = gft->para.num;
  
  for (i=0; i<num; i++)
  {
    gft->para.par[i+gft->para.num] = par[i];
    gft->para.fix[i+gft->para.num] = fix[i];
    gft->para.err[i+gft->para.num] = 0.0;
    gft->para.tmp[i+gft->para.num] = 0.0;
  } 

  gft->numf += 1;
  gft->para.num += num;
  
  return GSL_SUCCESS;
}
/*==========================================================================*/
int gsl_fit_print(gsl_fit *gft, char *filename)
{
  int    i, j, k;
  double tmp;
  FILE  *fp;
        
  /* Sanity checks */
  if (gft==NULL) {
    GSL_ERROR("Attempt to print data from a null gsl_fit", \
               GSL_EFAULT); }

  fp = fopen(filename, "w");
  if (fp==NULL) {
    GSL_ERROR("Cannot open file for writing", \
               GSL_ENOMEM); }

  /* loop over the data sets */
  for (k=0; k<gft->numd; k++)
  {
    /* loop over the data points in data set k */
    for (i=0; i<gft->data[k].num; i++)
    {
      /* loop over the model functions */
      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);
        }
      }
      fprintf(fp,"%f\t%f\n", gft->data[k].xpt[i], tmp);
    }
    fprintf(fp,"\n");
  }
    
  fclose(fp);

  return GSL_SUCCESS;
}
/*==========================================================================*/
int gsl_fit_get_result(gsl_fit *gft, \
                       unsigned int function, \
                       double *par, double *err, unsigned int num)
{
  int i;
  
  /* Sanity checks */
  if (gft==NULL) {
    GSL_ERROR(" attempt to get results from a null gsl_fit", \
               GSL_EINVAL); }
  if (function>=gft->numf) {
    GSL_ERROR(" attempt to get results from non-existant function", \
               GSL_EINVAL); }
  if (num!=gft->func[function].num) {
    GSL_ERROR(" wrong number of parameters to get", \
               GSL_EINVAL); }

  for (i=0; i<num; i++)
  {
    par[i] = (gft->para.par+gft->func[function].idx)[i];
    err[i] = (gft->para.err+gft->func[function].idx)[i];
  }
  return GSL_SUCCESS;
}
/*==========================================================================*/
int gsl_fit_do(gsl_fit *gft)
{
  int status;
  
  status = gft->meth.solver_function_call(gft);

  return status;
}
/*==========================================================================*/
void gsl_fit_free(gsl_fit *gft)
{
  int k;
  
  if (gft == NULL) return;

  for(k=0; k<gft->numd; k++)
  {
    free(gft->data[k].xpt);
    free(gft->data[k].ypt);
    free(gft->data[k].sig);
  }
  if (gft->numd) free(gft->data);
  
  if (gft->para.num)
  {
    free(gft->para.par);
    free(gft->para.fix);
    free(gft->para.err);
    free(gft->para.tmp);
  }

  if (gft->numf)
  {
    free(gft->func);
  }
  
  free(gft);
  /*==============*/
  return;
}

---