asa_usr.c

simulated annealing code ASA

{
  double uniform_test, curr_cost_temp;
#if USER_ACCEPT_ASYMP_EXP
  double x, q, delta_cost;
#endif

#if ASA_TEMPLATE                /* ASA cost index */
  /* Calculate the current ASA cost index.  This could be useful
     to define a new schedule for the cost temperature, beyond
     simple changes that can be made using USER_COST_SCHEDULE. */

  int index;
  double k_temperature, quench, y;
  double xparameter_dimension;

#if QUENCH_COST
  quench = USER_OPTIONS->User_Quench_Cost_Scale[0];
#else
  quench = 1.0;
#endif /* QUENCH_COST */
  xparameter_dimension = (double) *parameter_dimension;
  for (index = 0; index < *parameter_dimension; ++index)
    if (fabs (parameter_upper_bound[index] - parameter_lower_bound[index]) <
        (double) EPS_DOUBLE)
      *xparameter_dimension -= 1.0;

  y = -F_LOG (USER_OPTIONS->Cost_Temp_Curr
              / USER_OPTIONS->Cost_Temp_Init) / USER_OPTIONS->Cost_Temp_Scale;

  k_temperature = F_POW (y, xparameter_dimension / quench);
#endif /* ASA cost index */

  uniform_test = randflt (USER_OPTIONS->Random_Seed);
  curr_cost_temp = USER_OPTIONS->Cost_Temp_Curr;

#if ASA_TEMPLATE
#if USER_COST_SCHEDULE
  curr_cost_temp =
    (USER_OPTIONS->Cost_Schedule (USER_OPTIONS->Cost_Temp_Curr,
                                  USER_OPTIONS) + (double) EPS_DOUBLE);
#else
  curr_cost_temp = USER_OPTIONS->Cost_Temp_Curr;
#endif
#endif /* ASA_TEMPLATE */

  /* You must add in your own test here.  If USER_ACCEPT_ASYMP_EXP
     also is TRUE here, then you can use the default
     Asymp_Exp_Param=1 to replicate the code in asa.c. */

#if USER_ACCEPT_ASYMP_EXP
#if USER_COST_SCHEDULE
  curr_cost_temp =
    (USER_OPTIONS->Cost_Schedule (USER_OPTIONS->Cost_Temp_Curr,
                                  USER_OPTIONS) + (double) EPS_DOUBLE);
#endif

  delta_cost = (current_cost - *(USER_OPTIONS->Last_Cost))
    / (curr_cost_temp + (double) EPS_DOUBLE);

  /* The following asymptotic approximation to the exponential
   * function, "Tsallis statistics," was proposed in
   * %A T.J.P. Penna
   * %T Traveling salesman problem and Tsallis statistics
   * %J Phys. Rev. E
   * %V 50
   * %N 6
   * %P R1-R3
   * %D 1994
   * While the use of the TSP for a test case is of dubious value (since
   * there are many special algorithms for this problem), the use of this
   * function is another example of how to control the rate of annealing
   * of the acceptance criteria.  E.g., if you require a more moderate
   * acceptance test, then negative q may be helpful. */

  q = USER_OPTIONS->Asymp_Exp_Param;
  if (fabs (1.0 - q) < (double) EPS_DOUBLE)
    x = MIN (1.0, (F_EXP (-delta_cost)));       /* Boltzmann test */
  else if ((1.0 - (1.0 - q) * delta_cost) < (double) EPS_DOUBLE)
    x = MIN (1.0, (F_EXP (-delta_cost)));       /* Boltzmann test */
  else
    x = MIN (1.0, F_POW ((1.0 - (1.0 - q) * delta_cost), (1.0 / (1.0 - q))));

  USER_OPTIONS->Prob_Bias = x;
  if (x >= uniform_test)
    USER_OPTIONS->User_Acceptance_Flag = TRUE;
  else
    USER_OPTIONS->User_Acceptance_Flag = FALSE;

#endif /* USER_ACCEPT_ASYMP_EXP */
}
#endif /* USER_ACCEPTANCE_TEST */

#if USER_GENERATING_FUNCTION
#if HAVE_ANSI
double
user_generating_distrib (LONG_INT * seed,
                         ALLOC_INT * parameter_dimension,
                         ALLOC_INT index_v,
                         double temperature_v,
                         double init_param_temp_v,
                         double temp_scale_params_v,
                         double parameter_v,
                         double parameter_range_v,
                         double *last_saved_parameter,
                         USER_DEFINES * USER_OPTIONS)
#else
double
user_generating_distrib (seed,
                         parameter_dimension,
                         index_v,
                         temperature_v,
                         init_param_temp_v,
                         temp_scale_params_v,
                         parameter_v,
                         parameter_range_v,
                         last_saved_parameter, USER_OPTIONS)
     LONG_INT *seed;
     ALLOC_INT *parameter_dimension;
     ALLOC_INT index_v;
     double temperature_v;
     double init_param_temp_v;
     double temp_scale_params_v;
     double parameter_v;
     double parameter_range_v;
     double *last_saved_parameter;
     USER_DEFINES *USER_OPTIONS;
#endif
{
#if ASA_TEMPLATE
  double x, y, z;

  /* This is the ASA distribution.  A slower temperature schedule can be
     obtained here, e.g., temperature_v = pow(temperature_v, 0.5); */

  x = randflt (seed);
  y = x < 0.5 ? -1.0 : 1.0;
  z = y * temperature_v * (F_POW ((1.0 + 1.0 / temperature_v),
                                  fabs (2.0 * x - 1.0)) - 1.0);

  x = parameter_v + z * parameter_range_v;

  return (x);
#endif /* ASA_TEMPLATE */
}
#endif /* USER_GENERATING_FUNCTION */

#if USER_REANNEAL_COST
#if HAVE_ANSI
int
user_reanneal_cost (double *cost_best,
                    double *cost_last,
                    double *initial_cost_temperature,
                    double *current_cost_temperature,
                    USER_DEFINES * USER_OPTIONS)
#else
int
user_reanneal_cost (cost_best,
                    cost_last,
                    initial_cost_temperature,
                    current_cost_temperature, USER_OPTIONS)
     double *cost_best;
     double *cost_last;
     double *initial_cost_temperature;
     double *current_cost_temperature;
     USER_DEFINES *USER_OPTIONS;
#endif /* HAVE_ANSI */
{
  int cost_test;
  double tmp_dbl;

#if ASA_TEMPLATE
  static int first_time = 1;
  static double save_last[3];
  double average_cost_last;

  if (first_time == 1) {
    first_time = 0;
    save_last[0] = save_last[1] = save_last[2] = *cost_last;
  }

  save_last[2] = save_last[1];
  save_last[1] = save_last[0];
  save_last[0] = *cost_last;
  average_cost_last =
    fabs ((save_last[0] + save_last[1] + save_last[2]) / 3.0);

  tmp_dbl = MAX (fabs (*cost_best), average_cost_last);
  tmp_dbl = MAX ((double) EPS_DOUBLE, tmp_dbl);
  *initial_cost_temperature = MIN (*initial_cost_temperature, tmp_dbl);

  /* This test can be useful if your cost function goes from a positive
     to a negative value, and you do not want to get get stuck in a local
     minima around zero due to the default in reanneal().  Pick any
     number instead of 0.0001 */
  tmp_dbl = MIN (fabs (*cost_last), fabs (*cost_best));
  if (tmp_dbl < 0.0001)
    cost_test = FALSE;
  else
    cost_test = TRUE;
#endif /* ASA_TEMPLATE */

  tmp_dbl = MAX (fabs (cost_last), fabs (cost_best));
  tmp_dbl = MAX ((double) EPS_DOUBLE, tmp_dbl);
  *initial_cost_temperature = MIN (*initial_cost_temperature, tmp_dbl);

  *current_cost_temperature =
    MAX (fabs (cost_last - cost_best), *current_cost_temperature);
  *current_cost_temperature =
    MAX ((double) EPS_DOUBLE, *current_cost_temperature);
  *current_cost_temperature =
    MIN (*current_cost_temperature, *initial_cost_temperature);

  cost_test = TRUE;

  return (cost_test);
}
#endif /* USER_REANNEAL_COST */

#if USER_REANNEAL_PARAMETERS
#if HAVE_ANSI
double
user_reanneal_params (double current_temp,
                      double tangent,
                      double max_tangent, USER_DEFINES * USER_OPTIONS)
#else
double
user_reanneal_params (current_temp, tangent, max_tangent, USER_OPTIONS)
     double current_temp;
     double tangent;
     double max_tangent;
     USER_DEFINES *USER_OPTIONS;
#endif /* HAVE_ANSI */
{
#if ASA_TEMPLATE
  double x;

  x = current_temp * (max_tangent / tangent);

  return (x);
#endif
}
#endif /* USER_REANNEAL_PARAMETERS */

#if SELF_OPTIMIZE

/***********************************************************************
* main
*	This is a sample calling program to self-optimize ASA
***********************************************************************/
#if HAVE_ANSI

#if ASA_LIB
int
asa_main (
#if ASA_TEMPLATE_LIB
           double *main_recur_cost_value,
           double *main_recur_cost_parameters, int *main_recur_exit_code
#endif
  )
#else /* ASA_LIB */
int
main (int argc, char **argv)
#endif                          /* ASA_LIB */
#else /* HAVE_ANSI */

#if ASA_LIB
int
asa_main (
#if ASA_TEMPLATE_LIB
           main_recur_cost_value,
           main_recur_cost_parameters, main_recur_exit_code
#endif
  )
#if ASA_TEMPLATE_LIB
     double *main_recur_cost_value;
     double *main_recur_cost_parameters;
     int *main_recur_exit_code;
#endif

#else /* ASA_LIB */
int
main (argc, argv)
     int argc;
     char **argv;
#endif /* ASA_LIB */

#endif /* HAVE_ANSI */
{

  /* seed for random number generator */
  LONG_INT *recur_rand_seed;
  char user_exit_msg[160];
  FILE *ptr_out;

#if RECUR_OPTIONS_FILE
  FILE *recur_ptr_options;
  char read_option[80];
  char read_if[4], read_FALSE[6], read_comm1[3], read_ASA_SAVE[9],
    read_comm2[3];
  int read_int;
#if INT_LONG
  LONG_INT read_long;
#endif
  double read_double;
#endif

  int *recur_exit_code;
#if ASA_LIB
#else
  int compile_cnt;
#endif
#if MULTI_MIN
  int multi_index;
  ALLOC_INT n_param;
#endif

  double *recur_parameter_lower_bound, *recur_parameter_upper_bound;
  double *recur_cost_parameters, *recur_cost_tangents, *recur_cost_curvature;
  double recur_cost_value;

  ALLOC_INT *recur_parameter_dimension;
  int *recur_parameter_int_real;
  int *recur_cost_flag;
  int recur_initialize_params_value;
  ALLOC_INT recur_v;
#if MY_TEMPLATE                 /* MY_TEMPLATE_recur_main_decl */
  /* add some declarations if required */
#endif

  USER_DEFINES *RECUR_USER_OPTIONS;

  if ((recur_parameter_dimension =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (user_exit_msg, "main()/asa_main(): recur_parameter_dimension");
    Exit_USER (user_exit_msg);
    return (-2);
  }
  if ((recur_exit_code = (int *) calloc (1, sizeof (int))) == NULL) {
    strcpy (user_exit_msg, "main()/asa_main(): recur_exit_code");
    Exit_USER (user_exit_msg);
    return (-2);
  }
  if ((recur_cost_flag = (int *) calloc (1, sizeof (int))) == NULL) {
    strcpy (user_exit_msg, "main()/asa_main(): recur_cost_flag");
    Exit_USER (user_exit_msg);
    return (-2);
  }

  if ((RECUR_USER_OPTIONS =
       (USER_DEFINES *) calloc (1, sizeof (USER_DEFINES))) == NULL) {
    strcpy (user_exit_msg, "main()/asa_main(): RECUR_USER_OPTIONS");
    Exit_USER (user_exit_msg);
    return (-2);
  }
#if RECUR_OPTIONS_FILE
  recur_ptr_options = fopen ("asa_opt_recur", "r");

  fscanf (recur_ptr_options, "%s%s%s%s%s",
          read_if, read_FALSE, read_comm1, read_ASA_SAVE, read_comm2);
  if (strcmp (read_if, "#if") || strcmp (read_FALSE, "FALSE") ||
      strcmp (read_comm1, "/*") || strcmp (read_ASA_SAVE, "ASA_SAVE") ||
      strcmp (read_comm2, "*/")) {
#if INCL_STDOUT
    printf ("\n\n*** EXIT not asa_opt_recur for this version *** \n\n");
#endif /* INCL_STDOUT */
    return (-6);
  }
#if INT_LONG
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%ld", &read_long);
  RECUR_USER_OPTIONS->Limit_Acceptances = read_long;
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%ld", &read_long);
  RECUR_USER_OPTIONS->Limit_Generated = read_long;
#else
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%d", &read_int);
  RECUR_USER_OPTIONS->Limit_Acceptances = read_int;
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%d", &read_int);
  RECUR_USER_OPTIONS->Limit_Generated = read_int;
#endif
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%d", &read_int);
  RECUR_USER_OPTIONS->Limit_Invalid_Generated_States = read_int;
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%lf", &read_double);
  RECUR_USER_OPTIONS->Accepted_To_Generated_Ratio = read_double;

  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%lf", &read_double);
  RECUR_USER_OPTIONS->Cost_Precision = read_double;
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%d", &read_int);
  RECUR_USER_OPTIONS->Maximum_Cost_Repeat = read_int;
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%d", &read_int);
  RECUR_USER_OPTIONS->Number_Cost_Samples = read_int;
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%lf", &read_double);
  RECUR_USER_OPTIONS->Temperature_Ratio_Scale = read_double;
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%lf", &read_double);
  RECUR_USER_OPTIONS->Cost_Parameter_Scale_Ratio = read_double;
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%lf", &read_double);
  RECUR_USER_OPTIONS->Temperature_Anneal_Scale = read_double;

  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%d", &read_int);
  RECUR_USER_OPTIONS->Include_Integer_Parameters = read_int;
  fscanf (recur_ptr_options, "%s", read_option);
  fscanf (recur_ptr_options, "%d", &read_int);
  RECUR_USER_OPTIONS->User_Initial_Parameters = read_int;
#if INT_AL