yanera_util.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 
============================================================================*/ 
/*============================================================================*/
#define _CRT_SECURE_NO_DEPRECATE

#include "yanera.h"
#include "yanera_util.h"
#include "yanera_xml_input.h"
#include "yanera_mdmin.h"
#include "yanera_lstsq.h"
#include "yanera_profile.h"
#include "yanera_quadrature.h"
#include "yanera_reflectivity.h"
#include "yanera_xml_output.h"
/*============================================================================*/
yanera_container *yanera_initialize(const char *filename)
{
  unsigned short    i, j;
  yanera_container *yanera;
  
  fprintf(stdout,"YANERA         : 2.0 beta\n");
  fprintf(stdout,"   Date        : %s\n", __DATE__ );
  fprintf(stdout,"               :\n");

  /*---------------------------------
   *  Allocate a new yanera struct.
   */
  yanera = (yanera_container *)malloc(sizeof(yanera_container));
  if (yanera == NULL) 
  {
    yanera_error("No memory for new yanera_container", \
                 __FILE__, __LINE__, NULL);
  }
  /*---------------------------------
   *
   */
  yanera->filename = (char *)calloc((1+strlen(filename)),sizeof(char));
  if (yanera->filename == NULL) 
  {
    yanera_error("No memory for new yanera XML filename", \
                 __FILE__, __LINE__, yanera);
  }
  strcpy(yanera->filename, filename);

  /*---------------------------------
   *  
   */
  yanera->type = MODEL_LAYER;
  yanera->parrattFunction = &calcParrattFromModel;
  yanera->profileFunctionReal = &layerProfileReal;
  yanera->profileFunctionImag = &layerProfileImag;
  /*---------------------------------
   *  
   */
  for (i=0; i<10; i++)
  {
    yanera->models_xml[i].data_idref  = NULL;
    yanera->models_xml[i].name        = NULL;
    yanera->models_xml[i].background  = -1;
    yanera->models_xml[i].polarized   = POLARIZED_NONE;
    yanera->models_xml[i].first_layer = NULL;
    yanera->models_xml[i].last_layer  = NULL;
    yanera->models_xml[i].slabs.number_of_slab_edges = 0;
    yanera->models_xml[i].slabs.number_of_slabs      = 0;
    for (j=0; j<USHRT_MAX; j++)
    {
      yanera->models_xml[i].slabs.zpos[j] = 0.0;
      yanera->models_xml[i].slabs.thik[j] = 0.0;
      yanera->models_xml[i].slabs.rsld[j] = 0.0;
      yanera->models_xml[i].slabs.isld[j] = 0.0;
    }
    
    yanera->models_complete[i].data_idref  = NULL;
    yanera->models_complete[i].name        = NULL;
    yanera->models_complete[i].background  = -1;
    yanera->models_complete[i].polarized   = POLARIZED_NONE;
    yanera->models_complete[i].first_layer = NULL;
    yanera->models_complete[i].last_layer  = NULL;
    yanera->models_complete[i].slabs.number_of_slab_edges = 0;
    yanera->models_complete[i].slabs.number_of_slabs      = 0;
    for (j=0; j<USHRT_MAX; j++)
    {
      yanera->models_complete[i].slabs.zpos[j] = 0.0;
      yanera->models_complete[i].slabs.thik[j] = 0.0;
      yanera->models_complete[i].slabs.rsld[j] = 0.0;
      yanera->models_complete[i].slabs.isld[j] = 0.0;
    }
    
    yanera->data[i].data_id   = NULL;
    yanera->data[i].file_name = NULL;
    yanera->data[i].n         = 0;
    yanera->data[i].q         = NULL;
    yanera->data[i].R         = NULL;
    yanera->data[i].e         = NULL;
    yanera->data[i].idx = 0;
    yanera->data[i].resolution.type      = RESOLUTION_NONE;
    yanera->data[i].resolution.value     = NULL;
    yanera->data[i].resolution.q         = NULL;
    yanera->data[i].resolution.n         = 0;
    yanera->data[i].resolution.file_name = NULL;
  }
  yanera->number_of_models = 0;
  yanera->number_of_data   = 0;
  yanera->idx              = 0;

  /*---------------------------------
   *  zero the paramter struct
   */
  yanera->parameters.p   = NULL;
  yanera->parameters.f   = NULL;
  yanera->parameters.min = NULL;
  yanera->parameters.max = NULL;
  yanera->parameters.con = NULL;
  yanera->parameters.n   = 0;
  yanera->parameters.tmp = NULL;
  yanera->parameters.idx = 0;
  
  /*---------------------------------
   *  zero the components struct
   */
  yanera->components  = NULL;
  
  /*---------------------------------
   *  Default misc options/values
   */
  yanera->misc.write_results           = NO;
  yanera->misc.write_profile           = YES;
  yanera->misc.write_reflectivity      = YES;
  yanera->misc.write_slabs             = NO;
  yanera->misc.q_min                   = 1.0E-05;
  yanera->misc.q_max                   = 0.2;
  yanera->misc.q_num                   = 201;
  yanera->misc.z_num                   = 201;
  yanera->misc.quadrature_error        = 0.5E-02;
  yanera->misc.quadrature_thik         = 100.0;
  yanera->misc.fit_do                  = NO;
  yanera->misc.fit_weighting           = WEIGHT_NONE;
  yanera->misc.fit_iterations          = 25;
  yanera->misc.fit_write_iterations    = 10;
  yanera->misc.fit_method              = LMSDR;
  yanera->misc.fit_method_mdmin        = \
                                  gsl_multimin_fdfminimizer_conjugate_fr; 
  yanera->misc.fit_method_mdmin_simplex = \
                                  gsl_multimin_fminimizer_nmsimplex; 
  yanera->misc.fit_method_lstsq        = gsl_multifit_fdfsolver_lmsder;
#ifdef HAVE_OOL
  yanera->misc.fit_method_ool          = ool_conmin_minimizer_spg;
#endif
  return yanera;
}
/*=========================================================================*/
void yanera_free(yanera_container *yanera)
{
  struct _postfix_queue_element *postfix1, *postfix2;
  short         i;
  yanera_layer *layer_pointer_1, *layer_pointer_2;

  
  free(yanera->filename);
  
  free(yanera->parameters.p);
  free(yanera->parameters.f);
  free(yanera->parameters.tmp);
  free(yanera->parameters.min);
  free(yanera->parameters.max);
  free(yanera->parameters.con);
  
  for (i=0; i<10; i++)
  {
    /* free layers */
    layer_pointer_1 = yanera->models_xml[i].first_layer;
    while (layer_pointer_1 != NULL)
    {
      free(layer_pointer_1->layer_id);
      free(layer_pointer_1->layer_idref);
      free(layer_pointer_1->component_id);
      free(layer_pointer_1->component_idref);
      free(layer_pointer_1->name);
      free(layer_pointer_1->func);
      free(layer_pointer_1->parm_names);
      free(layer_pointer_1->parm_values);
        
      layer_pointer_2 = layer_pointer_1;
      layer_pointer_1 = layer_pointer_1->next;
      free(layer_pointer_2);
    };
    
    free(yanera->models_xml[i].data_idref);
    free(yanera->models_xml[i].name);
  };

  for (i=0; i<10; i++)
  {
    /* free layers */
    layer_pointer_1 = yanera->models_complete[i].first_layer;
    while (layer_pointer_1 != NULL)
    {
      free(layer_pointer_1->layer_id);
      free(layer_pointer_1->layer_idref);
      free(layer_pointer_1->component_id);
      free(layer_pointer_1->component_idref);
      free(layer_pointer_1->name);
      free(layer_pointer_1->func);
      if (layer_pointer_1->type == LAYER_FUNCTION)
      {
        postfix1 = layer_pointer_1->func_postfix->front;
        while (postfix1 != NULL)
        { 
          postfix2 = postfix1;
          postfix1 = postfix1->next;
          free(postfix2);
        }
      }
      free(layer_pointer_1->func_postfix);
      free(layer_pointer_1->parm_names);
      free(layer_pointer_1->parm_values);
      layer_pointer_2 = layer_pointer_1;
      layer_pointer_1 = layer_pointer_1->next;
      free(layer_pointer_2);
    };
    
    free(yanera->models_complete[i].data_idref);
    free(yanera->models_complete[i].name);
  };


  /* free components */
  layer_pointer_1 = yanera->components;
  while (layer_pointer_1 != NULL)
  {
      free(layer_pointer_1->layer_id);
      free(layer_pointer_1->layer_idref);
      free(layer_pointer_1->component_id);
      free(layer_pointer_1->component_idref);
      free(layer_pointer_1->name);
      free(layer_pointer_1->func);
      free(layer_pointer_1->parm_names);
      free(layer_pointer_1->parm_values);
    
    layer_pointer_2 = layer_pointer_1;
    layer_pointer_1 = layer_pointer_1->next;
    free(layer_pointer_2);
  };

  /* free datafiles */
  for (i=0; i<10; i++)
  {
    free(yanera->data[i].file_name);
    free(yanera->data[i].data_id);
    free(yanera->data[i].resolution.value);
    free(yanera->data[i].resolution.q);
    free(yanera->data[i].resolution.file_name);
    free(yanera->data[i].q);
    free(yanera->data[i].R);
    free(yanera->data[i].e);
  };
 
  free(yanera);
  yanera = NULL;
}
/*============================================================================*/
void yanera_error(const char       *reason, \
                  const char       *file,   \
                  int               line,   \
                  yanera_container *yanera)
{
  fprintf(stdout, "               :\n");
  fprintf(stdout, "YANERA ERROR   : %s\n", reason);
  fprintf(stdout, "   In file     : %s\n", file);
  fprintf(stdout, "   At line     : %d\n", line);
  fprintf(stdout, "\n");
  
  if (yanera != NULL) yanera_free(yanera);
  
  exit(EXIT_FAILURE);
}
/*============================================================================*/
void yanera_report(yanera_container *yanera, \
                  int iteration, \
                  void *solver)
{
  int j, i;
  gsl_multifit_fdfsolver    *mfs;
  gsl_multimin_fdfminimizer *mms;
  gsl_multimin_fminimizer   *mms2; 
#ifdef HAVE_OOL
  ool_conmin_minimizer      *oos;
#endif

  /*===========*/
  if (iteration >= 0)
  {
    fprintf(stderr,"   Iteration   :%3u \n", iteration);
    /*===========*/
    if ((yanera->misc.fit_method == LMSDR) || \
        (yanera->misc.fit_method == LMDR)  || \
        (yanera->misc.fit_method == CLMDR) || \
        (yanera->misc.fit_method == CLMSDR))
    {
      if (solver != NULL)
      { 
        mfs = (gsl_multifit_fdfsolver *)solver;
        for (j=0; j<yanera->parameters.n; j++)
        { yanera->parameters.p[j] = gsl_vector_get(mfs->x, j); }
        fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera));
      }
    }
    if ((yanera->misc.fit_method == CONJUGATE_FR) || \
        (yanera->misc.fit_method == CONJUGATE_PR) || \
        (yanera->misc.fit_method == VECTOR_BFGS)  || \
        (yanera->misc.fit_method == STEEP_DESC)   || \
        (yanera->misc.fit_method == NMSIMPLEX))
    { fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera)); }
#ifdef HAVE_OOL
    if ((yanera->misc.fit_method == OOL_SPG)      || \
        (yanera->misc.fit_method == OOL_GENCAN))
    { fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera)); }
#endif
#ifdef HAVE_LEVMAR       
    if (yanera->misc.fit_method == LEVMAR)
    { fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera)); }
#endif        
  }

  /*===========*/
  if (iteration == -1)
  {
    fprintf(stderr,"YANERA FINISHED: \n");
    /*===========*/
    if (solver != NULL)
    {
      if ((yanera->misc.fit_method == LMSDR) || \
          (yanera->misc.fit_method == LMDR)  || \
          (yanera->misc.fit_method == CLMDR) || \
          (yanera->misc.fit_method == CLMSDR))
      { 
        mfs = (gsl_multifit_fdfsolver *)solver;
        /*
        fprintf(stderr,"CHISQ     : %12.6e\n", gsl_blas_dnrm2(mfs->f));
        */
        for (j=0; j<yanera->parameters.n; j++)
        { yanera->parameters.p[j] = gsl_vector_get(mfs->x, j); }
        fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera));
      }
      if ((yanera->misc.fit_method == CONJUGATE_FR) || \
          (yanera->misc.fit_method == CONJUGATE_PR) || \
          (yanera->misc.fit_method == VECTOR_BFGS)  || \
          (yanera->misc.fit_method == STEEP_DESC))
      { 
        mms = (gsl_multimin_fdfminimizer *)solver;
        for (j=0; j<yanera->parameters.n; j++)
        { yanera->parameters.p[j] = gsl_vector_get(mms->x, j); }
        fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera));
      }  
      if (yanera->misc.fit_method == NMSIMPLEX)
      { 
        mms2 = (gsl_multimin_fminimizer *)solver;
        for (j=0, i=0; j<yanera->parameters.n; j++)
        { 
          if (yanera->parameters.f[j] == NO)
          {
            yanera->parameters.p[j] = gsl_vector_get(mms2->x, i);
            i++;
          }
        }
        fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera));
      }
#ifdef HAVE_OOL
      if ((yanera->misc.fit_method == OOL_SPG) || \
          (yanera->misc.fit_method == OOL_GENCAN))
      { 
        oos = (ool_conmin_minimizer *)solver;
        for (j=0; j<yanera->parameters.n; j++)
        { yanera->parameters.p[j] = gsl_vector_get(oos->x, j); }
        fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera));
      } 
#endif
    }
    if (yanera->misc.fit_method == SIMAN)
    {
      fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera));
    }
#ifdef HAVE_LEVMAR
    if (yanera->misc.fit_method == LEVMAR)
    {
      fprintf(stderr,"   Chi-squared : %12.6e\n", mdmin_chisq(yanera));
    }
#endif
  }
  if (((fmod(iteration,yanera->misc.fit_write_iterations) == 0) && \
      (iteration > 0)) || \
      (iteration == -1))
  {
    fprintf(stderr,"   Writing out :");
    /*===========*/
    /*---------------------------------
     * Write the profiles to file.
     */
    if (yanera->misc.write_profile == YES)
    {
      fprintf(stderr,"...profile");
      yanera_write_profile(yanera); 
    }
    /*---------------------------------
     * Do slab profiles to file.
     */
    if (yanera->misc.write_slabs == YES)
    {
      if (yanera->type != MODEL_LAYER)
      {
        yanera_write_slabs(yanera); 
        fprintf(stderr,"...slabs");
      }
    }
   /*---------------------------------
    * Write reflectivity to file
    */
    if (yanera->misc.write_reflectivity == YES)
    {
      yanera_write_reflectivity(yanera);
      fprintf(stderr,"...reflectivity");
    }
    /*---------------------------------
    * Write an XML file of the results
    */
    if (yanera->misc.write_results == YES)  
    {
      if (yanera->misc.fit_do == YES) yanera_find_errors(yanera);
      yanera_write_xml(yanera);
      fprintf(stderr,"...results\n");
    }
  }
  
  return;
}
/*============================================================================*/
int yanera_find_errors(yanera_container *yanera)
{
  int          i, count;
  gsl_matrix  *covar;
  gsl_matrix  *jacob;
  gsl_vector  *guess;
  
  for (i=0, count=0; i<yanera->number_of_models; i++) 
  {
    count += yanera->data[i].n;
  }
  jacob = gsl_matrix_alloc(count,yanera->parameters.n);
  
  covar = gsl_matrix_alloc(yanera->parameters.n, yanera->parameters.n);
  
  guess = gsl_vector_alloc(yanera->parameters.n);
  for (i=0; i<yanera->parameters.n; i++) 
  {
    gsl_vector_set(guess, i, yanera->parameters.p[i]);
  }
    
  lstsq_df(guess, (void *)yanera, jacob);
  gsl_multifit_covar(jacob, 0.0, covar);
  
  for (i=0; i<yanera->parameters.n; i++)
  {
    yanera->parameters.tmp[i] = gsl_matrix_get(covar,i,i);
  }
  
  gsl_matrix_free(covar);
  gsl_matrix_free(jacob);
  gsl_vector_free(guess);
  
  return GSL_SUCCESS;
}

---