yanera_ool.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_OOL

#include "yanera.h"
#include "yanera_ool.h"
#include "yanera_mdmin.h"
#include "yanera_util.h"
/*============================================================================*/
ool_diff_Hv_accel *Hv_accel;
/*============================================================================*/
int yanera_ool_main(yanera_container *yanera)
{ 
  int                   i, count, status;
  ool_conmin_function     ool_func;
  ool_conmin_constraint   ool_constraint;
  ool_conmin_minimizer   *solvr;
  gsl_vector             *guess;
  ool_conmin_spg_parameters    spgP;
  ool_conmin_gencan_parameters gencanP;
 
  ool_func.n   = yanera->parameters.n;
  ool_func.f   = &ool_f;
  ool_func.df  = &ool_df;
  ool_func.fdf = &ool_fdf;
  ool_func.Hv  = &ool_Hv;
  ool_func.params = (void *)yanera;
  Hv_accel = ool_diff_Hv_accel_alloc(yanera->parameters.n);
  
  /* Set the constraints */
  ool_constraint.n = yanera->parameters.n;
  ool_constraint.L = gsl_vector_alloc( ool_constraint.n );
  ool_constraint.U = gsl_vector_alloc( ool_constraint.n );
  for (i=0; i<yanera->parameters.n; i++) 
  {
    gsl_vector_set(ool_constraint.L, i, yanera->parameters.min[i]);
    gsl_vector_set(ool_constraint.U, i, yanera->parameters.max[i]);
    if (yanera->parameters.f[i] == YES)
    {
      gsl_vector_set(ool_constraint.L, i, yanera->parameters.p[i]);
      gsl_vector_set(ool_constraint.U, i, yanera->parameters.p[i]);
    }
  }
  
  /* 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]);
  }

  solvr = ool_conmin_minimizer_alloc(yanera->misc.fit_method_ool, \
                                     yanera->parameters.n);
  if (yanera->misc.fit_method == OOL_SPG)
  {
    ool_conmin_parameters_default(yanera->misc.fit_method_ool, \
                                 (void*)(&spgP));
    /*spgP.alphamax = 1.0e+3;
    spgP.alphamin = 1.0e-3;
    spgP.tol      = 1.0e-5;*/
  }  
  else if (yanera->misc.fit_method == OOL_GENCAN)
  {
    ool_conmin_parameters_default(yanera->misc.fit_method_ool, \
                                 (void*)(&gencanP));
  }
  
  /*==================================================*/
  count = 0;
  if (yanera->misc.fit_method == OOL_SPG)
  {
    ool_conmin_minimizer_set(solvr, &ool_func, &ool_constraint, guess, \
                             (void*)(&spgP) );
  }
  else if (yanera->misc.fit_method == OOL_GENCAN)
  {
    ool_conmin_minimizer_set(solvr, &ool_func, &ool_constraint, guess, \
                             (void*)(&gencanP) );
  }
  yanera_report(yanera,count,solvr);

  do 
  {
    count++;
     
    ool_conmin_minimizer_iterate(solvr);
    status = ool_conmin_is_optimal(solvr);
    if (status == OOL_SUCCESS) 
    {
      fprintf(stderr,"STATUS: SUCCESS\n");
    }
    
    yanera_report(yanera, count, solvr);

  } while((status == GSL_CONTINUE) && (count < yanera->misc.fit_iterations));
  
  yanera_report(yanera, -1, solvr);

  ool_conmin_minimizer_free(solvr);
  gsl_vector_free(guess);
  gsl_vector_free(ool_constraint.L);
  gsl_vector_free(ool_constraint.U);
  ool_diff_Hv_accel_free (Hv_accel);
  return OOL_SUCCESS;
}
/*============================================================================*/
void ool_fdf(const gsl_vector *parameters, void *yanera_pointer,
             double *f, gsl_vector *df )
{
  *f = ool_f(parameters, yanera_pointer);
  ool_diff_g(ool_f, parameters, yanera_pointer, df, 1.E-6);
}
/*============================================================================*/
double ool_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 ool_df(const gsl_vector *parameters, void *yanera_pointer, gsl_vector *df)
{
  ool_diff_g(ool_f, parameters, yanera_pointer, df, 1.E-6);
}
/*============================================================================*/
void ool_Hv(const gsl_vector *parameters, void *yanera_pointer, \
            const gsl_vector *V, gsl_vector *hv)
{
  ool_diff_Hv(Hv_accel, ool_df, parameters, yanera_pointer, V, hv, 1.E-6);
}

#endif

---