yanera_levmar.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 
=============================================================================*/
#ifdef HAVE_LEVMAR

#include "yanera.h"
#include "lm.h"
#include "yanera_levmar.h"
#include "yanera_quadrature.h"
#include "yanera_reflectivity.h"
#include "yanera_util.h"
/* ========================================================================== */
void levmar_f( double *p, double *y, int np, int ny, void *yanera_pointer)
{
  int                i, j;
  int                y_index;
  double             tmp;
  yanera_container  *yanera;
  /*
  ** 
  */
  yanera = (yanera_container *)yanera_pointer;
  /*
  **  Load the params back into the parameter array.
  */
  if (yanera->parameters.n != np) 
  {
    yanera_error("Error in levmar.", __FILE__, __LINE__, yanera);
  }
  
  for(i=0; i<yanera->parameters.n; i++)
  {
    yanera->parameters.p[i] = p[i];
  }
  /*
  ** Loop over the data sets and models
  */
  y_index = 0;
  for (j=0; j<yanera->number_of_models; j++)
  {
    /*
    ** Redo quadrature for these model types
    */
    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].e[i];
      y[y_index] = tmp;
      y_index++; 
      if (y_index > ny) 
      {
        yanera_error("Error in levmar.", __FILE__, __LINE__, yanera);
      }
    }
    /*
    ** 
    */
  }
  
  return;
}
/* ========================================================================== */
double levmar_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].e[i];
  
  return tmp;
}
/* ========================================================================== */
void levmar_df(double *p, double *j, int np, int ny, void *yanera_pointer)
{
  static short       iter = 0;
  int                i, l, k;
  int                j_index;
  double             res, err;
  gsl_function       F;
  yanera_container  *yanera;
  /*
   * Assuming 1 iteration is equal to 1 jacobian evaluation, call
   * yanera_report from here
   */
  yanera = (yanera_container *)yanera_pointer;
  yanera_report(yanera, iter++, NULL);
  /*
   *  Load the params into the temporary parameter array.
   */
  if (yanera->parameters.n != np)
  {
    yanera_error("Error in levmar.", __FILE__, __LINE__, yanera);
  }
  for(i=0; i<yanera->parameters.n; i++)
  {
    yanera->parameters.tmp[i] = p[i];
  }
  /*
   * Functions for gsl_deriv_central
   */
  F.function = &levmar_derivatives;
  F.params   = yanera_pointer;
  /*
   * Loop over the data sets
   */
  j_index = 0;
  for (k=0; k<yanera->number_of_models; k++)
  {
    /*
    **  The jacobian mtirx J_il = d(f_i)/d(a_l)
    **  => i is the data index of f
    **  => l 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 (l=0; l<yanera->parameters.n; l++) 
      {
        yanera->parameters.idx = l;
        /*
        ** If fixed, J_ij = 0
        */
        if (yanera->parameters.f[l] == YES) res=0.0;
        else gsl_deriv_central(&F, yanera->parameters.tmp[l], 1e-5, &res, &err);
        j[j_index] = res;
        j_index++;
      }
    }
  }
    
  return;
}
/* ========================================================================== */
int yanera_levmar_main(yanera_container *yanera)
{
  int           i, j, count;
  double        opts[4], info[LM_INFO_SZ];
  double        tmp, *y, *p, *work;


  /*
   * Levmar expects a single vector of all the data.
   * I have no idea what option \mu is, but LEVMAR works well
   * by calculating the smallest data point times error = tmp. 
   */
  count = 0;
  for (i=0; i<yanera->number_of_models; i++)
  { count += yanera->data[i].n; }
  y = (double *)malloc(count*sizeof(double));
  if (y == NULL)
  {
    yanera_error("Error in levmar.", __FILE__, __LINE__, yanera);
  }
  count = 0;
  tmp = 10000.5;
  for (i=0; i<yanera->number_of_models; i++)
  {
    for (j=0; j<yanera->data[i].n; j++)
    {
      y[count] = yanera->data[i].R[j] / yanera->data[i].e[j];
      if (tmp > yanera->data[i].R[j] * yanera->data[i].e[j])
      {
        tmp = yanera->data[i].R[j] * yanera->data[i].e[j];
      }
      count++; 
    }
  }
  /*
   * Vector p is the starting guess, and is also used as workspace,
   * so it was a mistake before to just pass yanera->parameters.p
   * as the inital guess.
   */
  p = (double *)malloc(yanera->parameters.n*sizeof(double));
  for (i=0; i<yanera->parameters.n; i++)
  {
    p[i] = yanera->parameters.p[i];
  }
  /*
   * Some user settable options. See the Levmar web page for what
   * they are. I don't know if setting opts[0] really helps or not.
   */
  opts[0] = tmp*10.0;
  opts[1] = 1E-8;
  opts[2] = 1E-8;
  opts[3] = 1E-8;
  /*
   *  Workspace provided to Levmar.
   */
  work = (double *)malloc(LM_DIF_WORKSZ(yanera->parameters.n,count)* \
                          sizeof(double));
  if (work == NULL)
  {
    yanera_error("Error in levmar.", __FILE__, __LINE__, yanera);
  }
  
  /*
   * Levmar routine. Returns number of iterations = jacobian evaluations
   * as far as I can tell.
   * Not sure how to get modulous for reporting.
   */
  tmp = dlevmar_bc_der(levmar_f, levmar_df,                            \
                       p, y,                        \
                       yanera->parameters.n, count,                    \
                       yanera->parameters.min, yanera->parameters.max, \
                       yanera->misc.fit_iterations, opts, info, \
                       work, NULL, (void *)yanera);
  /*
   * Get the answers to the fit and report completed.
   */                  
  for (i=0; i<yanera->parameters.n; i++)
  {
    yanera->parameters.p[i] = p[i];
  }
  yanera_report(yanera, -1, NULL);
  /*
   * Optional information
  fprintf(stderr, "RET:%g\n", tmp);
  fprintf(stderr, "I 0:%g\n", info[0]);
  fprintf(stderr, "I 1:%g\n", info[1]);
  fprintf(stderr, "I 2:%g\n", info[2]);
  fprintf(stderr, "I 3:%g\n", info[3]);
  fprintf(stderr, "I 4:%g\n", info[4]);
  fprintf(stderr, "I 5:%g\n", info[5]);
  fprintf(stderr, "I 6:%g\n", info[6]);
  fprintf(stderr, "I 7:%g\n", info[7]);
  fprintf(stderr, "I 8:%g\n", info[8]);
  fprintf(stderr, "c  :%d\n", count);
  fprintf(stderr, "n  :%d\n", yanera->parameters.n);
  */
  
  free(p);
  free(y);
  free(work);
  
  return (GSL_SUCCESS);
}

#endif

---