gsl_mdmin.c

Go to the documentation of this file.

/*============================================================================*/
#include "gsl_fit.h"
/*============================================================================*/
int mdmin_main(void *void_gsl_fit_struct)
{
  int                         count;
  int                         status;
  gsl_vector                 *guess;
  gsl_multimin_function_fdf   mdmin_func;
  gsl_multimin_fdfminimizer  *solvr;
  gsl_fit                    *gft;
  
  gft = (gsl_fit *)void_gsl_fit_struct;
  
  /* Inital guess */
  guess = gsl_vector_alloc(gft->para.num);
  for (count=0; count<gft->para.num; count++) 
  {
    gsl_vector_set(guess, count, gft->para.par[count]);
  }

  mdmin_func.f      = &mdmin_f;
  mdmin_func.df     = &mdmin_df;
  mdmin_func.fdf    = &mdmin_fdf;
  mdmin_func.n      = gft->para.num;
  mdmin_func.params = (void *)gft;
  
  solvr = gsl_multimin_fdfminimizer_alloc(gft->meth.mdmin_type, \
                                          gft->para.num);
  
  /*==================================================*/
  count = 0;
  gsl_multimin_fdfminimizer_set(solvr, &mdmin_func, guess, 0.01, 1e-6);
  mdmin_status(gft, count, solvr);

  do
  {
    count++;
    
    status = gsl_multimin_fdfminimizer_iterate(solvr);
    if (status != GSL_SUCCESS) 
    {
      fprintf(stderr,"STATUS: %s\n", gsl_strerror(status));
      break;
    }
    
    mdmin_status(gft, count, solvr);
    
    status = gsl_multimin_test_gradient(solvr->gradient, 1e-6);
    
    if (status == GSL_SUCCESS) 
    {
      fprintf(stderr,"STATUS: %s\n", gsl_strerror(status));
    }
    
  } while ((status == GSL_CONTINUE) && (count < gft->maxi));

  mdmin_status(gft, -1, solvr);

  gsl_multimin_fdfminimizer_free(solvr);
  gsl_vector_free(guess);
  
  return GSL_SUCCESS;
}
/*============================================================================*/
void mdmin_fdf(const gsl_vector *parameters, \
               void *void_gsl_fit_struct, \
               double *f, gsl_vector *df)
{
  *f = mdmin_f(parameters, void_gsl_fit_struct);
  mdmin_df(parameters, void_gsl_fit_struct, df);
  return;
}
/*============================================================================*/
double mdmin_f(const gsl_vector *parameters, \
               void *void_gsl_fit_struct)
{
  int     i;
  double  tmp;
  gsl_fit   *gft;
  
  gft = (gsl_fit *)void_gsl_fit_struct;
  
  /*
  **  Load the params back into the tmp parameter array
  */
  for(i=0; i<gft->para.num; i++)
  {
    gft->para.tmp[i] = gsl_vector_get(parameters, i);
  }

  tmp = mdmin_chisq(gft);
  
  return tmp;
}
/*============================================================================*/
void mdmin_df(const gsl_vector *parameters, \
              void *void_gsl_fit_struct, \
              gsl_vector *df)
{
  int                i;
  double             res, err;
  gsl_function       F;
  gsl_fit *gft;
  
  gft = (gsl_fit *)void_gsl_fit_struct;
  
  /*
  **  Load the params back into the parameter array
  */
  for(i=0; i<gft->para.num; i++)
  {
    gft->para.par[i] = gsl_vector_get(parameters, i);
  }

  /*
  ** Pass the gsl_fit as a void pointer to the
  ** function calculating the derivative.
  */
  F.params   = void_gsl_fit_struct;
  F.function = &mdmin_derivatives;
  
  /*
  **  d(chisq)/d(p_i)
  **  => i is the parameter index
  ** If the parameter is to be held fixed, then
  ** set the derivative to zero.
  */
  for (i=0; i<gft->para.num; i++) 
  {
    gft->para.idx = i;
    if (gft->para.fix[i]) res=0.0;
    else gsl_deriv_central(&F, gft->para.par[i], 1e-4, &res, &err);
    gsl_vector_set(df, i, res);
  }

  return;
}
/*============================================================================*/
double mdmin_derivatives(double p, void *void_gsl_fit_struct)
{
  int                i;
  double             tmp;
  gsl_fit *gft;
  
  gft = (gsl_fit *)void_gsl_fit_struct;
  
  /*
  **  We need a complete set of paramters, 
  **  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  = mdmin_chisq(gft);
  
  return tmp;
}
/*============================================================================*/
double mdmin_chisq(gsl_fit *gft)
{
  int    i, j, k;
  double tmp;
  double residual = 0.0;
  int    numdata;
  
  /*
  ** Loop over the data sets with index k.
  */
  numdata = 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.tmp + gft->func[j].idx, \
                                    k, j);    
        }
        else
        {
          tmp += gft->func[j].fun_s(gft->data[k].xpt[i], \
                                    gft->para.tmp + gft->func[j].idx);    
        }
      }
      tmp = (tmp - gft->data[k].ypt[i]) / gft->data[k].sig[i];    
      residual += tmp * tmp;
      numdata++;
    }
  } 
  
  tmp =  residual / (numdata - gft->para.num);
  return tmp;
}
/*============================================================================*/
void mdmin_status(gsl_fit *gft, \
                  int iteration, \
                  gsl_multimin_fdfminimizer *s)
{
  int      j;
  double   chisq;

  if (iteration >= 0)
  {
    fprintf(stderr,"ITERATION:%3u ", iteration);
    fprintf(stderr,"===================\n");
  }
  else
  {
    fprintf(stderr,"COMPLETE ");
    fprintf(stderr,"===================\n");
  }
  
  /*===========*/
  for (j=0; j<gft->para.num; j++)
  {
    gft->para.par[j] = gsl_vector_get(s->x, j);
    gft->para.tmp[j] = gsl_vector_get(s->x, j);
  }
  chisq = mdmin_chisq(gft);
  fprintf(stderr,"ChiSq: %12.6e\n", chisq);
  for (j=0; j<gft->para.num; j++)
  {
    fprintf(stderr,"%3d: %12.6f\n", j, gsl_vector_get(s->x, 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] = 0.0;
    }
    gft->chis = chisq;
  } 
  
  return;
}

---