asa.c

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

FALSE *valid_state_generated_flag\n");
#endif
  } else {
    /* let asa generate valid initial parameters */
    repeated_invalid_states = 0;
    OPTIONS->Sequential_Parameters = *start_sequence - 1;
    do {
      ++(*number_invalid_generated_states);
      generate_new_state (user_random_generator,
                          seed,
                          parameter_minimum,
                          parameter_maximum, current_user_parameter_temp,
#if USER_GENERATING_FUNCTION
                          initial_user_parameter_temp,
                          temperature_scale_parameters,
#endif
                          number_parameters,
                          parameter_type,
                          current_generated_state, last_saved_state, OPTIONS);
      *valid_state_generated_flag = TRUE;
#if USER_ACCEPTANCE_TEST
      OPTIONS->User_Acceptance_Flag = TRUE;
      OPTIONS->Cost_Acceptance_Flag = FALSE;
#endif
      current_generated_state->cost =
        user_cost_function (current_generated_state->parameter,
                            parameter_minimum,
                            parameter_maximum,
                            tangents,
                            curvature,
                            number_parameters,
                            parameter_type,
                            valid_state_generated_flag, exit_status, OPTIONS);
      if (cost_function_test
          (current_generated_state->cost,
           current_generated_state->parameter, parameter_minimum,
           parameter_maximum, number_parameters, xnumber_parameters) == 0) {
        *exit_status = INVALID_COST_FUNCTION;
        goto EXIT_ASA;
      }
      ++repeated_invalid_states;
      if (repeated_invalid_states > OPTIONS->Limit_Invalid_Generated_States) {
        *exit_status = TOO_MANY_INVALID_STATES;
        goto EXIT_ASA;
      }
    }
    while (*valid_state_generated_flag == FALSE);
    --(*number_invalid_generated_states);
  }                             /* OPTIONS->User_Initial_Parameters */

  /* set all states to the last one generated */
  VFOR (index_v) {
#if DROPPED_PARAMETERS
    /* ignore parameters that have too small a range */
    if (PARAMETER_RANGE_TOO_SMALL (index_v))
      continue;
#endif
    best_generated_state->parameter[index_v] =
      last_saved_state->parameter[index_v] =
      current_generated_state->parameter[index_v];
  }

  /* set all costs to the last one generated */
  best_generated_state->cost = last_saved_state->cost =
    current_generated_state->cost;

  *accepted_to_generated_ratio = ONE;

  /* do not calculate curvatures initially */
  *curvature_flag = FALSE;

#if ASA_PRINT
  fprintf (ptr_asa_out,
           "temperature_scale = %*.*g\n",
           G_FIELD, G_PRECISION, temperature_scale);
#if RATIO_TEMPERATURE_SCALES
#if ASA_PRINT_INTERMED
  VFOR (index_v) {
    fprintf (ptr_asa_out,
#if INT_ALLOC
             "temperature_scale_parameters[%d] = %*.*g\n",
#else
#if INT_LONG
             "temperature_scale_parameters[%ld] = %*.*g\n",
#else
             "temperature_scale_parameters[%d] = %*.*g\n",
#endif
#endif
             index_v,
             G_FIELD, G_PRECISION, temperature_scale_parameters[index_v]);
  }
#endif
#else
  fprintf (ptr_asa_out,
           "temperature_scale_parameters[0] = %*.*g\n",
           G_FIELD, G_PRECISION, temperature_scale_parameters[0]);
#endif /* RATIO_TEMPERATURE_SCALES */
  fprintf (ptr_asa_out,
           "*temperature_scale_cost = %*.*g\n",
           G_FIELD, G_PRECISION, *temperature_scale_cost);
  fprintf (ptr_asa_out, "\n\n");

#if ASA_PRINT_INTERMED
  print_state (parameter_minimum,
               parameter_maximum,
               tangents,
               curvature,
               current_cost_temperature,
               current_user_parameter_temp,
               accepted_to_generated_ratio,
               number_parameters,
               curvature_flag,
               number_accepted,
               index_cost_acceptances,
               number_generated,
               number_invalid_generated_states,
               last_saved_state, best_generated_state, ptr_asa_out, OPTIONS);
#endif
  fprintf (ptr_asa_out, "\n");

  fflush (ptr_asa_out);
#endif

#if ASA_SAMPLE
#if ASA_PRINT
  fprintf (ptr_asa_out,
           ":SAMPLE:   n_accept   cost        cost_temp    bias_accept    \
 aver_weight\n");
  fprintf (ptr_asa_out,
           ":SAMPLE:   index      param[]     temp[]       bias_gener[]   \
 range[]\n");
#endif
#endif

  /* reset the current cost and the number of generations performed */
  *number_invalid_generated_states = 0;
  *best_number_generated_saved =
    *number_generated = *recent_number_generated = 0;
  OPTIONS->N_Generated = *number_generated;
  VFOR (index_v) {
    /* ignore parameters that have too small a range */
    if (PARAMETER_RANGE_TOO_SMALL (index_v))
      continue;
    index_parameter_generations[index_v] = 1;
  }
#if USER_ACCEPTANCE_TEST
  OPTIONS->User_Acceptance_Flag = FALSE;
  OPTIONS->Cost_Acceptance_Flag = FALSE;
#endif

#if ASA_QUEUE
#if ASA_PRINT
#if INT_ALLOC
  fprintf (ptr_asa_out, "OPTIONS->Queue_Size = %d\n", OPTIONS->Queue_Size);
#else
#if INT_LONG
  fprintf (ptr_asa_out, "OPTIONS->Queue_Size = %ld\n", OPTIONS->Queue_Size);
#else
  fprintf (ptr_asa_out, "OPTIONS->Queue_Size = %d\n", OPTIONS->Queue_Size);
#endif
#endif
  VFOR (index_v) {
    fprintf (ptr_asa_out,
#if INT_ALLOC
             "Queue_Resolution[%d] = %*.*g\n",
#else
#if INT_LONG
             "Queue_Resolution[%ld] = %*.*g\n",
#else
             "Queue_Resolution[%d] = %*.*g\n",
#endif
#endif
             index_v,
             G_FIELD, G_PRECISION, OPTIONS->Queue_Resolution[index_v]);
  }
#endif /* ASA_PRINT */

  /* fill arrays to check allocated memory */
  for (queue = 0; queue < queue_size_tmp; ++queue) {
    VFOR (index_v) {
      if (PARAMETER_RANGE_TOO_SMALL (index_v)) {
        continue;
      }
      queue_v = index_v + queue * (LONG_INT) (*number_parameters);
      save_queue_param[queue_v] = current_generated_state->parameter[index_v];
    }
    save_queue_cost[queue] = current_generated_state->cost;
    save_queue_flag[queue] = *valid_state_generated_flag;
  }
  save_queue = save_queue_indx = 0;
#endif /* ASA_QUEUE */

#if ASA_RESOLUTION
#if ASA_PRINT
  VFOR (index_v) {
    fprintf (ptr_asa_out,
#if INT_ALLOC
             "Coarse_Resolution[%d] = %*.*g\n",
#else
#if INT_LONG
             "Coarse_Resolution[%ld] = %*.*g\n",
#else
             "Coarse_Resolution[%d] = %*.*g\n",
#endif
#endif
             index_v,
             G_FIELD, G_PRECISION, OPTIONS->Coarse_Resolution[index_v]);
  }
#endif /* ASA_PRINT */
#endif /* ASA_RESOLUTION */

#if MULTI_MIN
  multi_sort[OPTIONS->Multi_Number] = OPTIONS->Multi_Number;
  multi_cost[OPTIONS->Multi_Number] = current_generated_state->cost;
  VFOR (index_v) {
    multi_params[OPTIONS->Multi_Number][index_v] =
      current_generated_state->parameter[index_v];
  }
  for (multi_index = 0; multi_index < OPTIONS->Multi_Number; ++multi_index) {
    multi_sort[multi_index] = multi_index;
    multi_cost[multi_index] = OPTIONS->Multi_Cost[multi_index] =
      current_generated_state->cost;
    VFOR (index_v) {
      multi_params[multi_index][index_v] =
        OPTIONS->Multi_Params[multi_index][index_v] =
        current_generated_state->parameter[index_v];
    }
  }
#endif /* MULTI_MIN */

  OPTIONS->Sequential_Parameters = *start_sequence - 1;

  /* MAIN ANNEALING LOOP */
  while (((*number_accepted <= OPTIONS->Limit_Acceptances)
          || (OPTIONS->Limit_Acceptances == 0))
         && ((*number_generated <= OPTIONS->Limit_Generated)
             || (OPTIONS->Limit_Generated == 0))) {

    tmp_var_db1 = -F_LOG ((OPTIONS->Temperature_Ratio_Scale));

    /* compute temperature scales */
    tmp_var_db2 = F_LOG (OPTIONS->Temperature_Anneal_Scale);
    temperature_scale = tmp_var_db1 *
      F_EXP (-tmp_var_db2 / *xnumber_parameters);

#if QUENCH_PARAMETERS
#if RATIO_TEMPERATURE_SCALES
    VFOR (index_v)
      temperature_scale_parameters[index_v] = tmp_var_db1 * F_EXP
#if QUENCH_PARAMETERS_SCALE
      (-(tmp_var_db2 * OPTIONS->User_Quench_Param_Scale[index_v])
#else
      (-(tmp_var_db2)
#endif
       / *xnumber_parameters)
      * OPTIONS->User_Temperature_Ratio[index_v];
#else
    VFOR (index_v)
      temperature_scale_parameters[index_v] = tmp_var_db1 * F_EXP
#if QUENCH_PARAMETERS_SCALE
      (-(tmp_var_db2 * OPTIONS->User_Quench_Param_Scale[index_v])
#else
      (-(tmp_var_db2)
#endif
       / *xnumber_parameters);
#endif /* RATIO_TEMPERATURE_SCALES */
#else /* QUENCH_PARAMETERS */
#if RATIO_TEMPERATURE_SCALES
    VFOR (index_v)
      temperature_scale_parameters[index_v] =
      tmp_var_db1 * F_EXP (-(tmp_var_db2) / *xnumber_parameters)
      * OPTIONS->User_Temperature_Ratio[index_v];
#else
    VFOR (index_v)
      temperature_scale_parameters[index_v] =
      tmp_var_db1 * F_EXP (-(tmp_var_db2) / *xnumber_parameters);
#endif /* RATIO_TEMPERATURE_SCALES */
#endif /* QUENCH_PARAMETERS */

#if USER_ACCEPTANCE_TEST
    OPTIONS->Cost_Temp_Scale =
#endif
      *temperature_scale_cost =
#if QUENCH_COST
#if QUENCH_COST_SCALE
      tmp_var_db1 * F_EXP (-(tmp_var_db2 * OPTIONS->User_Quench_Cost_Scale[0])
#else
      tmp_var_db1 * F_EXP (-(tmp_var_db2)
#endif
                           / *xnumber_parameters) *
      OPTIONS->Cost_Parameter_Scale_Ratio;
#else /* QUENCH_COST */
      tmp_var_db1 * F_EXP (-(tmp_var_db2)
                           / *xnumber_parameters) *
      OPTIONS->Cost_Parameter_Scale_Ratio;
#endif /* QUENCH_COST */

    /* CALCULATE NEW TEMPERATURES */

    /* calculate new parameter temperatures */
    VFOR (index_v) {
      /* skip parameters with too small a range */
      if (PARAMETER_RANGE_TOO_SMALL (index_v))
        continue;

      log_new_temperature_ratio =
        -temperature_scale_parameters[index_v] *
        F_POW ((double) index_parameter_generations[index_v],
#if QUENCH_PARAMETERS
               OPTIONS->User_Quench_Param_Scale[index_v]
#else /* QUENCH_PARAMETERS */
               ONE
#endif /* QUENCH_PARAMETERS */
               / *xnumber_parameters);
      /* check (and correct) for too large an exponent */
      log_new_temperature_ratio = EXPONENT_CHECK (log_new_temperature_ratio);
      current_user_parameter_temp[index_v] =
        initial_user_parameter_temp[index_v]
        * F_EXP (log_new_temperature_ratio);

#if NO_PARAM_TEMP_TEST
      if (current_user_parameter_temp[index_v] < (double) EPS_DOUBLE)
        current_user_parameter_temp[index_v] = (double) EPS_DOUBLE;
#else
      /* check for too small a parameter temperature */
      if (current_user_parameter_temp[index_v] < (double) EPS_DOUBLE) {
        *exit_status = P_TEMP_TOO_SMALL;
        *index_exit_v = index_v;
        goto EXIT_ASA;
      }
#endif
    }

    /* calculate new cost temperature */
    log_new_temperature_ratio =
      -*temperature_scale_cost * F_POW ((double) *index_cost_acceptances,
#if QUENCH_COST
                                        OPTIONS->User_Quench_Cost_Scale[0]
#else
                                        ONE
#endif
                                        / *xnumber_parameters);
    log_new_temperature_ratio = EXPONENT_CHECK (log_new_temperature_ratio);
#if USER_ACCEPTANCE_TEST
    OPTIONS->Cost_Temp_Curr = OPTIONS->Cost_Temp_Init =
#endif
      *current_cost_temperature = *initial_cost_temperature
      * F_EXP (log_new_temperature_ratio);

#if NO_COST_TEMP_TEST
    if (*current_cost_temperature < (double) EPS_DOUBLE)
#if USER_ACCEPTANCE_TEST
      OPTIONS->Cost_Temp_Curr =
#endif
        *current_cost_temperature = (double) EPS_DOUBLE;
#else
    /* check for too small a cost temperature */
    if (*current_cost_temperature < (double) EPS_DOUBLE) {
      *exit_status = C_TEMP_TOO_SMALL;
      goto EXIT_ASA;
    }
#endif

#if ASA_SAVE
    if (asa_read == TRUE && OPTIONS->Asa_Recursive_Level == asa_recursive_max) {
      if (OPTIONS->Asa_Recursive_Level > 0)
        sprintf (asa_save_comm, "asa_save_%d", OPTIONS->Asa_Recursive_Level);
      else