gsl_siman.c

Go to the documentation of this file.

/*============================================================================*/
#include "gsl_fit.h"
/*===================================================================
 *    siman_energy 
 */
int siman_main(void *void_gsl_fit_struct)
{
  double              st, ft; 
  gsl_rng            *rng;
  gsl_siman_params_t  params;
  gsl_fit            *gft;
  
  gft = (gsl_fit *)void_gsl_fit_struct;

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

  /* starting and final temperatures */
  st = 2.0 * siman_energy(void_gsl_fit_struct);
  ft = st / 1.0e6;

  /*
  **  The number of steps in temperature, n_T, is given by
  **  n_T = ln(st/ft)/ln(mu_t)
  **
  **  The total number of trials, N, is given by
  **  N=n_T*iters_fixed_T
  */

  params.iters_fixed_T = 500;
  params.k             =   1.0;
  params.t_initial     =   st;
  params.t_min         =   ft;
  params.mu_t          =   1.01;
  params.step_size     =   0.01;
    
  fprintf(stderr, "Starting %d simulated annealing steps.\n\n",
          (int)floor(log(st/ft)/log(params.mu_t)));
          
  gsl_siman_solve(rng, void_gsl_fit_struct, \
                  siman_energy,    \
                  siman_step,      \
                  siman_metric,    \
                  siman_print,     \
                  siman_copy,      \
                  siman_construct, \
                  siman_destroy,   \
                  0, params);


  gft->chis = siman_energy(void_gsl_fit_struct);
  
  gsl_rng_free(rng);
  
  return GSL_SUCCESS;
}
/*===================================================================
 *    siman_energy 
 */
double siman_energy(void *void_gsl_fit_struct)
{
  int i, j, k;
  double tmp;
  double residual = 0.0;
  int    numdata;
  gsl_fit  *gft;

  gft = (gsl_fit *)void_gsl_fit_struct;
  
  /*
  ** Loop over the data sets with index k.
  */
  numdata = 0;
  for (k=0; k<gft->numd; k++)
  {
    /*
    ** Loop over the data points with index i
    */ 
    for (i=0; i<gft->data[k].num; i++)
    {
      /*
      **   Loop over the functions with index j
      */
      tmp = 0.0;
      for (j=0; j<gft->numf; j++) 
      {
        if (gft->type == MULTIPLE)
        {
          tmp += gft->func[j].fun_m(gft->data[k].xpt[i], \
                                    gft->para.par + gft->func[j].idx, \
                                    k, j);    
        }
        else
        {
          tmp += gft->func[j].fun_s(gft->data[k].xpt[i], \
                                    gft->para.par + gft->func[j].idx);    
        }
      }
      tmp = (tmp - gft->data[k].ypt[i]) / gft->data[k].sig[i];    
      residual += tmp * tmp;
      numdata++;
    }
  } 
  tmp =  residual / (numdata - gft->para.num);
  return tmp;
}
/*===================================================================
 *    siman_step 
 */
void siman_step(const gsl_rng *rng, void *void_gsl_fit_struct, double step_size)
{
  gsl_fit *gft;
  int i;
  double s;
  
  gft = (gsl_fit *)void_gsl_fit_struct;
  
  for (i=0; i<gft->para.num; i++)
  {
    if (gft->para.fix[i] == 0) 
    {
      s = gsl_rng_uniform(rng);
      gft->para.par[i] = step_size * (2.0*s-1.0) + gft->para.par[i];
    }
  }
  return;
}
/*===================================================================
 *    siman_copy 
 */
void siman_copy(void *source, void *dest)
{
  int i, j;
  
  gsl_fit  *gft_old, *gft_new;
  
  gft_old = (gsl_fit *)source;
  gft_new = (gsl_fit *)dest;

  gft_new->numd = gft_old->numd;
  for (i=0; i<gft_old->numd; i++)
  {
    gft_new->data[i].idx = gft_old->data[i].idx;
    gft_new->data[i].num = gft_old->data[i].num;
    for (j=0; j<gft_old->data[i].num; j++)
    {
      gft_new->data[i].xpt[j] = gft_old->data[i].xpt[j];
      gft_new->data[i].ypt[j] = gft_old->data[i].ypt[j];
      gft_new->data[i].sig[j] = gft_old->data[i].sig[j];
    }
  }
  
  gft_new->para.num = gft_old->para.num;
  for (i=0; i<gft_old->para.num; i++)
  {
    gft_new->para.par[i] = gft_old->para.par[i];
    gft_new->para.fix[i] = gft_old->para.fix[i];
    gft_new->para.err[i] = gft_old->para.err[i];
    gft_new->para.tmp[i] = gft_old->para.tmp[i];
  }
  gft_new->para.idx = gft_old->para.idx;

  gft_new->numf = gft_old->numf;
  for (i=0; i<gft_old->numf; i++)
  {
    gft_new->func[i].idx = gft_old->func[i].idx;
    gft_new->func[i].num = gft_old->func[i].num;
    gft_new->func[i].fun_s = gft_old->func[i].fun_s;
    gft_new->func[i].fun_m = gft_old->func[i].fun_m;
  }
  
  gft_new->meth.solver_function_call = gft_old->meth.solver_function_call;
  gft_new->meth.mdmin_type           = gft_old->meth.mdmin_type;
  gft_new->meth.lstsq_type           = gft_old->meth.lstsq_type;

  gft_new->chis = gft_old->chis;  
  gft_new->maxi = gft_old->maxi;  
  gft_new->indx = gft_old->indx;  
  gft_new->type = gft_old->type;
  return; 

}
/*===================================================================
 *    siman_construct 
 */
void *siman_construct(void *void_gsl_fit_struct)
{
  int i,j;
  
  gsl_fit  *gft_old, *gft_new;
  
  gft_old = (gsl_fit *)void_gsl_fit_struct;
  gft_new = (gsl_fit *)malloc(sizeof(gsl_fit));

  gft_new->numd = gft_old->numd;
  gft_new->data = (struct _data *)malloc(gft_old->numd * sizeof(struct _data));

  for (i=0; i<gft_old->numd; i++)
  {
    gft_new->data[i].xpt = (double *)malloc(gft_old->data[i].num * \
                                            sizeof(double));
    gft_new->data[i].ypt = (double *)malloc(gft_old->data[i].num * \
                                            sizeof(double));
    gft_new->data[i].sig = (double *)malloc(gft_old->data[i].num * \
                                            sizeof(double));
    gft_new->data[i].idx = gft_old->data[i].idx;
    gft_new->data[i].num = gft_old->data[i].num;
    for (j=0; j<gft_old->data[i].num; j++)
    {
      gft_new->data[i].xpt[j] = gft_old->data[i].xpt[j];
      gft_new->data[i].ypt[j] = gft_old->data[i].ypt[j];
      gft_new->data[i].sig[j] = gft_old->data[i].sig[j];
    }
  }
  
  gft_new->para.num = gft_old->para.num;
  gft_new->para.par = (double *)malloc(gft_old->para.num*sizeof(double));
  gft_new->para.fix = (unsigned int *)malloc(gft_old->para.num*\
                                             sizeof(unsigned int));
  gft_new->para.err = (double *)malloc(gft_old->para.num*sizeof(double));
  gft_new->para.tmp = (double *)malloc(gft_old->para.num*sizeof(double));
  for (i=0; i<gft_old->para.num; i++)
  {
    gft_new->para.par[i] = gft_old->para.par[i];
    gft_new->para.fix[i] = gft_old->para.fix[i];
    gft_new->para.err[i] = gft_old->para.err[i];
    gft_new->para.tmp[i] = gft_old->para.tmp[i];
  }
  gft_new->para.idx = gft_old->para.idx;

  gft_new->numf = gft_old->numf;
  gft_new->func     = malloc(gft_old->numf*sizeof(struct _func));
  for (i=0; i<gft_old->numf; i++)
  {
    gft_new->func[i].idx = gft_old->func[i].idx;
    gft_new->func[i].num = gft_old->func[i].num;
    gft_new->func[i].fun_s = gft_old->func[i].fun_s;
    gft_new->func[i].fun_m = gft_old->func[i].fun_m;
  }
  
  gft_new->meth.solver_function_call = gft_old->meth.solver_function_call;
  gft_new->meth.mdmin_type           = gft_old->meth.mdmin_type;
  gft_new->meth.lstsq_type           = gft_old->meth.lstsq_type;

  gft_new->chis = gft_old->chis;  
  gft_new->maxi = gft_old->maxi;  
  gft_new->indx = gft_old->indx;  
  gft_new->type = gft_old->type;
  
  return (void *)gft_new; 
}
/*===================================================================
 *    siman_destroy 
 */
void siman_destroy(void *void_gsl_fit_struct)
{
  int i;
  gsl_fit  *gft;

  if (void_gsl_fit_struct == NULL) return;

  gft = (gsl_fit *)void_gsl_fit_struct;

  for(i=0; i<gft->numd; i++)
  {
    free(gft->data[i].xpt);
    free(gft->data[i].ypt);
    free(gft->data[i].sig);
  }
  if (gft->numd) free(gft->data);
  
  if (gft->para.num)
  {
    free(gft->para.par);
    free(gft->para.fix);
    free(gft->para.err);
    free(gft->para.tmp);
  }

  if (gft->numf)
  {
    free(gft->func);
  }
  
  free(void_gsl_fit_struct);
  
  return;
}
/*===================================================================
 *    siman_metric 
 */
double siman_metric(void *void_gsl_fit_struct1, void *void_gsl_fit_struct2)
{
  gsl_fit *gft1;
  gsl_fit *gft2;
  double sum, tmp;
  int i;

  gft1 = (gsl_fit *)void_gsl_fit_struct1;
  gft2 = (gsl_fit *)void_gsl_fit_struct2;
  sum = 0.0;
  for (i=0; i<gft1->para.num; i++)
  {
    /*sum += fabs(gft1->para.par[i] - gft2->para.par[i]);*/
    tmp = (gft1->para.par[i] - gft2->para.par[i]);
    sum += tmp * tmp;
  }
  return sqrt(sum);
}
/*===================================================================
 *    siman_print 
 */
void siman_print(void *void_gsl_fit_struct)
{
  gsl_fit *gft;
  int i;
  
  gft = (gsl_fit *)void_gsl_fit_struct;
  
  for (i=0; i<gft->para.num; i++)
  {
    printf(" %10.5g ", gft->para.par[i]);
  }
  return;
}

---