yanera_mdmin.c

Go to the documentation of this file.

/*============================================================================
Copyright 2006 Thad Harroun, Kevin Yager

This file is part of "yanera 2.0".

"yanera 2.0" is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

"yanera 2.0" is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with "yanera 2.0"; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA 
=============================================================================*/
#include "yanera.h"
#include "yanera_mdmin.h"
#include "yanera_quadrature.h"
#include "yanera_reflectivity.h"
#include "yanera_util.h"
/*============================================================================*/
int yanera_mdmin_main(yanera_container *yanera)
{
  int                         i, count;
  int                         status;
  gsl_vector                 *guess;
  gsl_multimin_function_fdf   mdmin_func;
  gsl_multimin_fdfminimizer  *solvr;
  
  /* Inital guess is the paramters as they were passed into here */
  guess = gsl_vector_alloc(yanera->parameters.n);
  for (i=0; i<yanera->parameters.n; i++) 
  {
    gsl_vector_set(guess, i, yanera->parameters.p[i]);
  }

  mdmin_func.f      = &mdmin_f;
  mdmin_func.df     = &mdmin_df;
  mdmin_func.fdf    = &mdmin_fdf;
  mdmin_func.n      = yanera->parameters.n;
  mdmin_func.params = (void *)yanera;
  
  solvr = gsl_multimin_fdfminimizer_alloc( \
             yanera->misc.fit_method_mdmin, yanera->parameters.n);
  
  /*==================================================*/
  count = 0;
  if ((yanera->misc.fit_method == STEEP_DESC) ||
      (yanera->misc.fit_method == VECTOR_BFGS))
    gsl_multimin_fdfminimizer_set(solvr, &mdmin_func, guess, 1.1, 1E-2);
  else
    gsl_multimin_fdfminimizer_set(solvr, &mdmin_func, guess, 1.1, 1e-4);
  yanera_report(yanera, count, solvr);

  do
  {
    count++;
    
    status = gsl_multimin_fdfminimizer_iterate(solvr);
    if (status != GSL_SUCCESS) 
    {
      fprintf(stderr,"STATUS: %s\n", gsl_strerror(status));
      break;
    }
    
    yanera_report(yanera, 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 < yanera->misc.fit_iterations));

  yanera_report(yanera, -1, solvr);

  gsl_multimin_fdfminimizer_free(solvr);
  gsl_vector_free(guess);
  
  return GSL_SUCCESS;
}
/*============================================================================*/
void mdmin_fdf(const gsl_vector *parameters, \
               void *yanera_pointer, \
               double *f, gsl_vector *df)
{
  *f = mdmin_f(parameters, yanera_pointer);
  mdmin_df(parameters, yanera_pointer, df);
  return;
}
/*============================================================================*/
double mdmin_f(const gsl_vector *parameters, \
               void *yanera_pointer)
{
  int     i;
  double  tmp;
  yanera_container  *yanera;
  /*
  ** r =>> residual
  */
  yanera = (yanera_container *)yanera_pointer;
  /*
  **  Load the params back into the parameter array.
  */
  for(i=0; i<yanera->parameters.n; i++)
  {
    yanera->parameters.p[i] = gsl_vector_get(parameters, i);
  }

  tmp = mdmin_chisq(yanera);
  
  return tmp;
}
/*============================================================================*/
void mdmin_df(const gsl_vector *parameters, \
              void *yanera_pointer, \
              gsl_vector *df)
{
  int                i;
  double             res, err;
  gsl_function       F;
  yanera_container  *yanera;
  /*
  **
  */
  yanera = (yanera_container *)yanera_pointer;
  /*
  **  Load the params into the temporary parameter array.
  */
  for(i=0; i<yanera->parameters.n; i++)
  {
    yanera->parameters.tmp[i] = gsl_vector_get(parameters, i);
  }
  /*
   *
   */
  F.function = &mdmin_derivatives;
  F.params   = yanera_pointer;
  /*
  **  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<yanera->parameters.n; i++) 
  {
    yanera->parameters.idx = i;
    if (yanera->parameters.f[i] == YES) res=0.0;
    else gsl_deriv_central(&F, yanera->parameters.tmp[i], 1e-6, &res, &err);
    gsl_vector_set(df, i, res);
  }

  return;
}
/*============================================================================*/
double mdmin_derivatives(double new_p, void *yanera_pointer)
{
  int                i;
  double             tmp;
  yanera_container  *yanera;
  /*
  ** p =>> varying parameter for derivative
  */
  yanera = (yanera_container  *)yanera_pointer; 
  /*
   *  We need move the parameters from 'tmp', back into 'p'
   *  with the derivative parameter 'new_p' 
   *  replacing the parameter at index 'y->parameters.idx'.
   *
   *  model and data index  = yanera->idx;
   *  parameter index       = yanera->parameters.idx;
   *  data point index      = yanera->data[yanera->idx].idx;
   */
  for(i=0; i<yanera->parameters.n; i++) 
  {
    yanera->parameters.p[i] = yanera->parameters.tmp[i];
  }
  yanera->parameters.p[yanera->parameters.idx] = new_p;

  tmp  = mdmin_chisq(yanera);
  
  return tmp;
}
/*============================================================================*/
double mdmin_chisq(yanera_container *yanera)
{
  int    i, j;
  double tmp;
  double residual;
  int    numdata;
  
  /*
  ** Loop over the data sets with index k.
  */
  residual = 0.0;
  numdata  = 0;
  for (j=0; j<yanera->number_of_models; j++)
  {
    /*
    ** Redo quadrature for this moel
    */
    if ((yanera->type == MODEL_COMPONENT) || \
        (yanera->type == MODEL_SLAB)      || \
        (yanera->type == MODEL_FUNCTION))
    { doQuadrature(&yanera->models_complete[yanera->idx], yanera); }
    /*
    ** Loop over the data points
    */ 
    for (i=0; i<yanera->data[j].n; i++)
    {
      tmp = calcReflectivity(yanera->data[j].q[i], \
                             &yanera->data[j].resolution, \
                             &yanera->models_complete[j], \
                             yanera);
      tmp = (tmp - yanera->data[j].R[i]) / yanera->data[j].e[i];
      residual += tmp * tmp;
      numdata++;
    }
    /*
    ** 
    */
  }
  
  tmp =  residual / (double)(numdata - yanera->parameters.n);
  return tmp;
}

---