asa.c

模拟退火的源代码。关于matlab的大家好好学习以下吧

/***********************************************************************
* Adaptive Simulated Annealing (ASA)
* Lester Ingber <[e-mail]>
* Copyright (c) 1993-2004 Lester Ingber.  All Rights Reserved.
* The LICENSE file must be included with ASA code.
***********************************************************************/

#define ASA_ID \
"/* $Id: asa.c,v 25.15 2004/09/23 18:10:46 ingber Exp ingber $ */"

#include "asa.h"

/***********************************************************************
* asa
*       This procedure implements the full ASA function optimization.
***********************************************************************/
#if HAVE_ANSI
double
asa (double (*user_cost_function)

      
     (double *, double *, double *, double *, double *, ALLOC_INT *, int *,
      int *, int *, USER_DEFINES *),
     double (*user_random_generator) (LONG_INT *), LONG_INT * seed,
     double *parameter_initial_final, double *parameter_minimum,
     double *parameter_maximum, double *tangents, double *curvature,
     ALLOC_INT * number_parameters, int *parameter_type,
     int *valid_state_generated_flag, int *exit_status,
     USER_DEFINES * OPTIONS)
#else
double
asa (user_cost_function,
     user_random_generator,
     seed,
     parameter_initial_final,
     parameter_minimum,
     parameter_maximum,
     tangents,
     curvature,
     number_parameters,
     parameter_type, valid_state_generated_flag, exit_status, OPTIONS)
     double (*user_cost_function) ();
     double (*user_random_generator) ();
     LONG_INT *seed;
     double *parameter_initial_final;
     double *parameter_minimum;
     double *parameter_maximum;
     double *tangents;
     double *curvature;
     ALLOC_INT *number_parameters;
     int *parameter_type;
     int *valid_state_generated_flag;
     int *exit_status;
     USER_DEFINES *OPTIONS;
#endif /* HAVE_ANSI */
{
#if USER_INITIAL_COST_TEMP
#if USER_REANNEAL_COST
#else
  int index_cost_constraint;    /* index cost functions averaged */
#endif /* USER_REANNEAL_COST */
#else /* USER_INITIAL_COST_TEMP */
  int index_cost_constraint;    /* index cost functions averaged */
#endif /* USER_INITIAL_COST_TEMP */

  int index_cost_repeat,        /* test OPTIONS->Cost_Precision when =
                                   OPTIONS->Maximum_Cost_Repeat */
    tmp_var_int, tmp_var_int1, tmp_var_int2;    /* temporary integers */

  ALLOC_INT index_v,            /* iteration index */
   *start_sequence;             /* initial OPTIONS->Sequential_Parameters
                                   used if >= 0 */
  double final_cost,            /* best cost to return to user */
    tmp_var_db, tmp_var_db1, tmp_var_db2;       /* temporary doubles */
  int *curvature_flag;
  FILE *ptr_asa_out;            /* file ptr to output file */

  /* The 3 states that are kept track of during the annealing process */
  STATE *current_generated_state, *last_saved_state, *best_generated_state;

#if ASA_SAVE
  FILE *ptr_save, *ptr_comm;
  int asa_read;
  char asa_save_comm[100];
#if ASA_SAVE_OPT
  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;
  FILE *ptr_save_opt;
#endif
#endif /* ASA_SAVE */

#if ASA_PIPE_FILE
  FILE *ptr_asa_pipe;
#endif

  int immediate_flag;           /* save Immediate_Exit */
  int asa_exit_value;

  double xnumber_parameters[1];

  /* The array of tangents (absolute value of the numerical derivatives),
     and the maximum |tangent| of the array */
  double *maximum_tangent;

  /* ratio of acceptances to generated points - determines when to
     test/reanneal */
  double *accepted_to_generated_ratio;

  /* temperature parameters */
  double temperature_scale, *temperature_scale_parameters;
  /* relative scalings of cost and parameters to temperature_scale */
  double *temperature_scale_cost;
  double *current_user_parameter_temp;
  double *initial_user_parameter_temp;
  double *current_cost_temperature;
  double *initial_cost_temperature;
  double log_new_temperature_ratio;     /* current *temp = initial *temp *
                                           exp(log_new_temperature_ratio) */
  ALLOC_INT *index_exit_v;      /* information for asa_exit */

  /* counts of generated states and acceptances */
  LONG_INT *index_parameter_generations;
  LONG_INT *number_generated, *best_number_generated_saved;
  LONG_INT *recent_number_generated, *number_accepted;
  LONG_INT *recent_number_acceptances, *index_cost_acceptances;
  LONG_INT *number_acceptances_saved, *best_number_accepted_saved;

  /* Flag indicates that the parameters generated were
     invalid according to the cost function validity criteria. */
  LONG_INT *number_invalid_generated_states;
  LONG_INT repeated_invalid_states;

#if ASA_QUEUE
  int queue_new;                /* flag to add new entry */
  int *save_queue_flag;         /* save valid_state_generated_flag */
  LONG_INT queue;               /* index of queue */
  LONG_INT queue_v;             /* index of parameters in queue */
  LONG_INT save_queue_test;     /* test if all parameters are present */
  LONG_INT save_queue;          /* last filled position in queue */
  LONG_INT save_queue_indx;     /* current position in queue */
  double *save_queue_cost, *save_queue_param;   /* saved states */
  ALLOC_INT queue_size_tmp;
#endif

#if MULTI_MIN
  int multi_index;
  int multi_test, multi_test_cmp, multi_test_dim;
  int *multi_sort;
  double *multi_cost;
  double **multi_params;
#endif /* MULTI_MIN */

#if ASA_PARALLEL
  LONG_INT *parallel_sort;
  LONG_INT index_parallel, sort_index;  /* count of parallel generated states */
  LONG_INT parallel_generated;  /* saved *recent_number_generated */
  LONG_INT parallel_block_max;  /* saved OPTIONS->Gener_Block_Max */
  STATE *gener_block_state;
#endif

  /* used to index repeated and recursive calls to asa */
  /* This assumes that multiple calls (>= 1) _or_ recursive
     calls are being made to asa */
  static int asa_open = FALSE;
  static int number_asa_open = 0;
  static int recursive_asa_open = 0;

  /* initializations */

  if ((curvature_flag = (int *) calloc (1, sizeof (int))) == NULL) {
    strcpy (exit_msg, "asa(): curvature_flag");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((maximum_tangent = (double *) calloc (1, sizeof (double))) == NULL) {
    strcpy (exit_msg, "asa(): maximum_tangent");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((accepted_to_generated_ratio =
       (double *) calloc (1, sizeof (double))) == NULL) {
    strcpy (exit_msg, "asa(): accepted_to_generated_ratio");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((temperature_scale_cost =
       (double *) calloc (1, sizeof (double))) == NULL) {
    strcpy (exit_msg, "asa(): temperature_scale_cost");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((current_cost_temperature =
       (double *) calloc (1, sizeof (double))) == NULL) {
    strcpy (exit_msg, "asa(): current_cost_temperature");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((initial_cost_temperature =
       (double *) calloc (1, sizeof (double))) == NULL) {
    strcpy (exit_msg, "asa(): initial_cost_temperature");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((index_exit_v = (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): index_exit_v");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((start_sequence = (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): start_sequence");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((number_generated =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): number_generated");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((best_number_generated_saved =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): best_number_generated_saved");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((recent_number_generated =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): recent_number_generated");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((number_accepted =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): number_accepted");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((recent_number_acceptances =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): recent_number_acceptances");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((index_cost_acceptances =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): index_cost_acceptances");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((number_acceptances_saved =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): number_acceptances_saved");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((best_number_accepted_saved =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): best_number_accepted_saved");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((number_invalid_generated_states =
       (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) {
    strcpy (exit_msg, "asa(): number_invalid_generated_states");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }

  if ((current_generated_state =
       (STATE *) calloc (1, sizeof (STATE))) == NULL) {
    strcpy (exit_msg, "asa(): current_generated_state");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((last_saved_state = (STATE *) calloc (1, sizeof (STATE))) == NULL) {
    strcpy (exit_msg, "asa(): last_saved_state");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  if ((best_generated_state = (STATE *) calloc (1, sizeof (STATE))) == NULL) {
    strcpy (exit_msg, "asa(): best_generated_state");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
#if ASA_PARALLEL
  if ((gener_block_state =
       (STATE *) calloc (OPTIONS->Gener_Block_Max, sizeof (STATE))) == NULL) {
    strcpy (exit_msg, "asa(): gener_block_state");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
  gener_block_state_qsort = gener_block_state;
  if ((parallel_sort =
       (LONG_INT *) calloc (OPTIONS->Gener_Block_Max,
                            sizeof (LONG_INT))) == NULL) {
    strcpy (exit_msg, "asa(): parallel_sort");
    Exit_ASA (exit_msg);
    *exit_status = CALLOC_FAILED;
    return (-1);
  }
#endif

  /* set default */
  ptr_asa_out = (FILE *) NULL;

  OPTIONS->Immediate_Exit = FALSE;

  if (asa_open == FALSE) {
    asa_open = TRUE;
    ++number_asa_open;
#if ASA_PRINT
    if (number_asa_open == 1) {
      /* open the output file */
#if USER_ASA_OUT
      if (!strcmp (OPTIONS->Asa_Out_File, "STDOUT")) {
#if INCL_STDOUT
        ptr_asa_out = stdout;
#endif /* INCL_STDOUT */
      } else {
#if ASA_SAVE
        ptr_asa_out = fopen (OPTIONS->Asa_Out_File, "a");
#else
        ptr_asa_out = fopen (OPTIONS->Asa_Out_File, "w");
#endif
      }
#else /* USER_ASA_OUT */
      if (!strcmp (ASA_OUT, "STDOUT")) {
#if INCL_STDOUT
        ptr_asa_out = stdout;
#endif /* INCL_STDOUT */
      } else {
#if ASA_SAVE
        ptr_asa_out = fopen (ASA_OUT, "a");
#else
        ptr_asa_out = fopen (ASA_OUT, "w");
#endif
      }
#endif /* USER_ASA_OUT */
    } else {
#if USER_ASA_OUT
      if (!strcmp (OPTIONS->Asa_Out_File, "STDOUT")) {