yanera_siman.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 Parratt; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA 
============================================================================*/ 
/*============================================================================*/
#include "yanera.h"
#include "yanera_siman.h"
#include "yanera_mdmin.h"
#include "yanera_util.h"
#include "yanera_xml_util.h"
/*============================================================================*/
int yanera_siman_main(yanera_container *yanera)
{
  double              start_temp, final_temp;
  const gsl_rng      *rng;
  gsl_siman_params_t  params;

  rng = gsl_rng_alloc(gsl_rng_rand);
  gsl_rng_set(rng, time(NULL));

  
  start_temp = 4.0*mdmin_chisq(yanera);
  final_temp = start_temp / 1e6;

  params.n_tries       = 200;
  params.iters_fixed_T =   yanera->misc.fit_write_iterations;
  params.step_size     =   0.1;
  params.k             =   1.0;
  params.t_initial     =   start_temp;
  params.mu_t          =   exp(log(1.0E6)/yanera->misc.fit_iterations);
  params.t_min         =   final_temp;
  
  gsl_siman_solve(rng, (void *)yanera, \
                  siman_energy, \
                  siman_step,   \
                  siman_metric, \
                  siman_print,  \
                  siman_copy, siman_copy_construct, siman_destroy,
                  0, params);
                  
  yanera_report(yanera,-1,NULL);
  return GSL_SUCCESS;
}
/*============================================================================*/
double siman_energy(void *yanera_pointer)
{
  double tmp;
  tmp = mdmin_chisq((yanera_container *)yanera_pointer);
  /*fprintf(stderr,"%f\n", tmp);*/
  return tmp;
}
/*============================================================================*/
void siman_step(const gsl_rng *rng, void *yanera_pointer, double step_size)
{
  int i;
  double s;
  yanera_container  *yanera;
  
  yanera = (yanera_container *)yanera_pointer;
  
  for (i=0; i<yanera->parameters.n; i++)
  {
    if (yanera->parameters.f[i] == NO) 
    {
      s = gsl_rng_uniform(rng);
      yanera->parameters.p[i] += step_size * (2.0*s-1.0);
      
      /* ======= */
      if ((yanera->parameters.con[i] == CONSTRAIN_MAX) || \
          (yanera->parameters.con[i] == CONSTRAIN_BOTH))
      {
        if (yanera->parameters.p[i] > yanera->parameters.max[i])
        { yanera->parameters.p[i] = yanera->parameters.max[i]; }
      }
      if ((yanera->parameters.con[i] == CONSTRAIN_MIN) || \
          (yanera->parameters.con[i] == CONSTRAIN_BOTH))
      {
        if (yanera->parameters.p[i] < yanera->parameters.min[i])
        { yanera->parameters.p[i] = yanera->parameters.min[i]; }
      }
      /* ======= */
    }
  }
  return;
}
/*============================================================================*/
double siman_metric(void *yanera_pointer1, void *yanera_pointer2)
{
  yanera_container  *yanera1;
  yanera_container  *yanera2;
  double sum;
  int i;

  yanera1 = (yanera_container *)yanera_pointer1;
  yanera2 = (yanera_container *)yanera_pointer2;
  sum = 0.0;
  for (i=0; i<yanera1->parameters.n; i++)
  {
    sum += fabs(yanera1->parameters.p[i] - yanera2->parameters.p[i]);
  }
  return sum;
}
/*============================================================================*/
void siman_print(void *yanera_pointer)
{
  int               i;
  yanera_container *yanera;
  
  yanera = (yanera_container *)yanera_pointer;
  
  for (i=0; i<yanera->parameters.n; i++)
  {
    if (yanera->parameters.f[i] == NO) 
    {
      printf("%8.4g", yanera->parameters.p[i]);
    }
  }
  return;
}
/*============================================================================*/
void siman_copy(void *source_yanera_pointer, void *dest_yanera_pointer)
{
  int               i;
  yanera_container *new_yanera, *yanera;
  new_yanera = (yanera_container *)dest_yanera_pointer;
      yanera = (yanera_container *)source_yanera_pointer;
  /*
   *
   */
  for (i=0; i<yanera->parameters.n; i++)
  {
    new_yanera->parameters.p[i] = yanera->parameters.p[i];
  }
  
  
  return;
}
/*============================================================================*/
void *siman_copy_construct(void *yanera_pointer)
{
  int               i, j;
  yanera_container *new_yanera, *yanera;
  yanera_layer     *new_layer,  *layer;
  /*
   *
   */
  new_yanera = yanera_initialize(NULL);
  yanera = (yanera_container *)yanera_pointer;
  /* 
   * 
   */
  new_yanera->type                = yanera->type;
  new_yanera->parrattFunction     = yanera->parrattFunction;
  new_yanera->profileFunctionReal = yanera->profileFunctionReal;
  new_yanera->profileFunctionImag = yanera->profileFunctionImag;
  /* 
   * 
   */
  new_yanera->number_of_models = yanera->number_of_models;
  for (i=0; i<yanera->number_of_models; i++)
  {
    new_yanera->data[i].q = malloc(yanera->data[i].n * sizeof(double));
    new_yanera->data[i].R = malloc(yanera->data[i].n * sizeof(double));
    new_yanera->data[i].e = malloc(yanera->data[i].n * sizeof(double));
    new_yanera->data[i].resolution.value = \
                         malloc(yanera->data[i].resolution.n * sizeof(double));
    new_yanera->data[i].resolution.q = \
                         malloc(yanera->data[i].resolution.n * sizeof(double));
    /* ----- */
    layer = yanera->models_complete[i].first_layer;
    while (layer != NULL)
    {
      new_layer = copyLayer(layer);
      linkLayerIntoModel(new_layer, &new_yanera->models_complete[i]);
      layer = layer->next;
    }
    /* ----- */
    new_yanera->models_complete[i].background = \
            yanera->models_complete[i].background;
    new_yanera->models_complete[i].polarized = \
            yanera->models_complete[i].polarized;
    /* --- */      
    new_yanera->data[i].n = yanera->data[i].n;
    for (j=0; j<yanera->data[i].n; j++)
    {
      new_yanera->data[i].q[j] = yanera->data[i].q[j];
      new_yanera->data[i].R[j] = yanera->data[i].R[j];
      new_yanera->data[i].e[j] = yanera->data[i].e[j];
    }
    new_yanera->data[i].resolution.type = yanera->data[i].resolution.type;
    new_yanera->data[i].resolution.n    = yanera->data[i].resolution.n;
    for (j=0; j<yanera->data[i].resolution.n; j++)
    {
      new_yanera->data[i].resolution.value[j] = \
                  yanera->data[i].resolution.value[j];
      new_yanera->data[i].resolution.q[j]     = \
                  yanera->data[i].resolution.q[j];
    }
  }
  /*
   *
   */
  new_yanera->parameters.p   = malloc(yanera->parameters.n * sizeof(double));
  new_yanera->parameters.f   = malloc(yanera->parameters.n * \
                                      sizeof(unsigned short));
  new_yanera->parameters.con = malloc(yanera->parameters.n * \
                               sizeof(unsigned short));
  new_yanera->parameters.min = malloc(yanera->parameters.n * sizeof(double));
  new_yanera->parameters.max = malloc(yanera->parameters.n * sizeof(double));  
  new_yanera->parameters.n   = yanera->parameters.n;
  for (i=0; i<yanera->parameters.n; i++)
  {
    new_yanera->parameters.p[i] = yanera->parameters.p[i];
    new_yanera->parameters.f[i] = yanera->parameters.f[i];
    new_yanera->parameters.min[i] = yanera->parameters.min[i];
    new_yanera->parameters.max[i] = yanera->parameters.max[i];
    new_yanera->parameters.con[i] = yanera->parameters.con[i];
  }
  /*
   *
   */
  new_yanera->misc.quadrature_error = yanera->misc.quadrature_error;
  new_yanera->misc.quadrature_thik = yanera->misc.quadrature_thik;
  new_yanera->misc.fit_weighting = yanera->misc.fit_weighting;
  new_yanera->misc.fit_iterations = yanera->misc.fit_iterations;
  new_yanera->misc.fit_iterations = yanera->misc.fit_write_iterations;  /*
   *
   */
  siman_copy(yanera_pointer, (void *)new_yanera);
  return (void *)new_yanera;
}
/*============================================================================*/
void siman_destroy(void *yanera_pointer)
{
  yanera_free((yanera_container *)yanera_pointer);
  return;
}

---