yanera_lstsq.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_lstsq.h"
#include "yanera_quadrature.h"
#include "yanera_reflectivity.h"
#include "yanera_util.h"
/*===========================================================================*/
int yanera_lstsq_main(yanera_container *yanera)
{
  unsigned int                        i;
  int                                 count, stop1, stop2;
  gsl_vector                         *guess;
  gsl_vector                         *gradt;
  gsl_matrix                         *covar;
  gsl_multifit_function_fdf           lstsq_func;
  gsl_multifit_fdfsolver             *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]);
  }

  /* count the number of data points */
  count = 0;
  for (i=0; i<yanera->number_of_models; i++)
  {
    count += yanera->data[i].n;
  }

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

  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      = yanera->parameters.n; /* number of paramters   */
  lstsq_func.params = (void *)yanera; /* pass the full yanera_container */
  
  solvr = gsl_multifit_fdfsolver_alloc(\
          yanera->misc.fit_method_lstsq, lstsq_func.n, lstsq_func.p);

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

  do
  { 
    count++;
    
    stop1 = gsl_multifit_fdfsolver_iterate(solvr);    
    /*gsl_multifit_covar(solvr->J, 0.0, covar);*/
    yanera_report(yanera, count, solvr);
    
    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 < yanera->misc.fit_iterations));

  /*==================================================*/
  /*gsl_multifit_covar(solvr->J, 0.0, covar);*/
  yanera_report(yanera, -1, solvr);
  
  /*==================================================*/

  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 *yanera_pointer, \
              gsl_vector *r, gsl_matrix *J)
{
  lstsq_f (parameters, yanera_pointer, r);
  lstsq_df(parameters, yanera_pointer, J);
  return GSL_SUCCESS;
}
/*===================================================================
 *    lstsq_f 
 */
int lstsq_f(const gsl_vector *parameters, \
            void *yanera_pointer, \
            gsl_vector *r)
{
  int                i, j;
  int                r_index;
  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);
  }
  /*
  ** Loop over the data sets and models
  */
  r_index = 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[j], 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];
      gsl_vector_set(r, r_index, tmp);
      r_index++; 
    }
    /*
    ** 
    */
  }
  
  return GSL_SUCCESS;
}
/*===================================================================
 *    lstsq_df 
 */
int lstsq_df(const gsl_vector *parameters, \
             void *yanera_pointer, \
             gsl_matrix *J)
{
  int                i, j, k;
  int                r_index;
  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 = &lstsq_derivatives;
  F.params   = yanera_pointer;
  /*
  ** Loop over the data sets
  */
  r_index = 0;
  for (k=0; k<yanera->number_of_models; 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
    */
    yanera->idx = k;
    /*
    ** Loop over the data points
    */
    for (i=0; i<yanera->data[k].n; i++)
    {
      yanera->data[k].idx = i;
      /*
      ** Loop over the parameters 
      */
      for (j=0; j<yanera->parameters.n; j++) 
      {
        yanera->parameters.idx = j;
        /*
        ** If fixed, J_ij = 0
        */
        if (yanera->parameters.f[j] == YES) res=0.0;
        else gsl_deriv_central(&F, yanera->parameters.tmp[j], 1e-5, &res, &err);
        gsl_matrix_set(J, r_index, j, res);
      }
      r_index++; 
    }
    /*
    **
    */
  }
    
  return GSL_SUCCESS;
}
/*===================================================================
 *    lstsq_derivatives 
 */
double lstsq_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;
  if ((yanera->type == MODEL_COMPONENT) ||
      (yanera->type == MODEL_SLAB) ||
      (yanera->type == MODEL_FUNCTION))
  { doQuadrature(&yanera->models_complete[yanera->idx], yanera); }
  
  tmp = 0.0;
  i = yanera->data[yanera->idx].idx;
  tmp = calcReflectivity(yanera->data[yanera->idx].q[i], \
                         &yanera->data[yanera->idx].resolution, \
                         &yanera->models_complete[yanera->idx], \
                          yanera);
  tmp = (tmp - yanera->data[yanera->idx].R[i]) /  \
               yanera->data[yanera->idx].e[i];
  return tmp;
}

---