/*********************************************************************** * Lester Ingber (copyright) (c) * See COPYING License in this directory * Date 1 Jan 93 ***********************************************************************/ #define ASA_ID \ "/* $Id: asa.c,v 2.3 1993/12/22 04:47:38 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 (*user_random_generator) (), 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, DEFINES * OPTIONS) #else double asa(user_cost_function, user_random_generator, 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) (); 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; DEFINES *OPTIONS; #endif { int index_cost_constraint, /* index initial cost functions averaged to set initial temperature */ temp_var_int, /* integer temporary value */ index_cost_repeat; /* test OPTIONS->COST_PRECISION when = OPTIONS->MAXIMUM_COST_REPEAT */ ALLOC_INT index_v; /* iteration index */ double sampled_cost_sum, /* temporary initial costs */ final_cost, /* best cost to return to user */ tmp_var_db1, tmp_var_db2; /* temporary doubles */ int *curvature_flag; int *quench_param_flag; int *quench_cost_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; /* Holds values of parameter increments */ double *delta_parameter; /* Holds values of quench_param scales */ double *quench_param; /* Holds value of quench_cost scales */ double *quench_cost; /* 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; /* used by EXPONENT_CHECK */ /* 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) exit(9); if ((quench_param_flag = (int *) calloc(1, sizeof(int))) == NULL) exit(9); if ((quench_cost_flag = (int *) calloc(1, sizeof(int))) == NULL) exit(9); if ((maximum_tangent = (double *) calloc(1, sizeof(double))) == NULL) exit(9); if ((accepted_to_generated_ratio = (double *) calloc(1, sizeof(double))) == NULL) exit(9); if ((temperature_scale_cost = (double *) calloc(1, sizeof(double))) == NULL) exit(9); if ((current_cost_temperature = (double *) calloc(1, sizeof(double))) == NULL) exit(9); if ((initial_cost_temperature = (double *) calloc(1, sizeof(double))) == NULL) exit(9); if ((index_exit_v = (ALLOC_INT *) calloc(1, sizeof(ALLOC_INT))) == NULL) exit(9); if ((number_generated = (LONG_INT *) calloc(1, sizeof(LONG_INT))) == NULL) exit(9); if ((best_number_generated_saved = (LONG_INT *) calloc(1, sizeof(LONG_INT))) == NULL) exit(9); if ((recent_number_generated = (LONG_INT *) calloc(1, sizeof(LONG_INT))) == NULL) exit(9); if ((number_accepted = (LONG_INT *) calloc(1, sizeof(LONG_INT))) == NULL) exit(9); if ((recent_number_acceptances = (LONG_INT *) calloc(1, sizeof(LONG_INT))) == NULL) exit(9); if ((index_cost_acceptances = (LONG_INT *) calloc(1, sizeof(LONG_INT))) == NULL) exit(9); if ((number_acceptances_saved = (LONG_INT *) calloc(1, sizeof(LONG_INT))) == NULL) exit(9); if ((best_number_accepted_saved = (LONG_INT *) calloc(1, sizeof(LONG_INT))) == NULL) exit(9); if ((number_invalid_generated_states = (LONG_INT *) calloc(1, sizeof(LONG_INT))) == NULL) exit(9); if ((current_generated_state = (STATE *) calloc(1, sizeof(STATE))) == NULL) exit(9); if ((last_saved_state = (STATE *) calloc(1, sizeof(STATE))) == NULL) exit(9); if ((best_generated_state = (STATE *) calloc(1, sizeof(STATE))) == NULL) exit(9); 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 ptr_asa_out = fopen(OPTIONS->asa_out_file, "w"); #else ptr_asa_out = fopen(ASA_OUT, "w"); #endif } else { #if USER_ASA_OUT ptr_asa_out = fopen(OPTIONS->asa_out_file, "a"); #else ptr_asa_out = fopen(ASA_OUT, "a"); #endif fprintf(ptr_asa_out, "\n\n\t\t number_asa_open = %d\n", number_asa_open); } #endif } else { ++recursive_asa_open; #if ASA_PRINT if (recursive_asa_open == 1) { /* open the output file */ #if USER_ASA_OUT ptr_asa_out = fopen(OPTIONS->asa_out_file, "w"); #else ptr_asa_out = fopen(ASA_OUT, "w"); #endif } else { #if USER_ASA_OUT ptr_asa_out = fopen(OPTIONS->asa_out_file, "a"); #else ptr_asa_out = fopen(ASA_OUT, "a"); #endif fprintf(ptr_asa_out, "\n\n\t\t recursive_asa_open = %d\n", recursive_asa_open); } #endif } #if ASA_PRINT /* print header information as defined by user */ fprintf(ptr_asa_out, "\t\tADAPTIVE SIMULATED ANNEALING\n\n"); fprintf(ptr_asa_out, "%s\n\n", ASA_ID); fprintf(ptr_asa_out, "OPTIONS_FILE = %d\n", (int) OPTIONS_FILE); fprintf(ptr_asa_out, "HAVE_ANSI = %d\n", (int) HAVE_ANSI); fprintf(ptr_asa_out, "IO_PROTOTYPES = %d\n", (int) IO_PROTOTYPES); fprintf(ptr_asa_out, "TIME_CALC = %d\n", (int) TIME_CALC); fprintf(ptr_asa_out, "TIME_STD = %d\n", (int) TIME_STD); fprintf(ptr_asa_out, "INT_LONG = %d\n", (int) INT_LONG); fprintf(ptr_asa_out, "INT_ALLOC = %d\n", (int) INT_ALLOC); fprintf(ptr_asa_out, "SMALL_FLOAT = %g\n", (double) SMALL_FLOAT); fprintf(ptr_asa_out, "MIN_DOUBLE = %g\n", (double) MIN_DOUBLE); fprintf(ptr_asa_out, "MAX_DOUBLE = %g\n", (double) MAX_DOUBLE); fprintf(ptr_asa_out, "EPS_DOUBLE = %g\n", (double) EPS_DOUBLE); fprintf(ptr_asa_out, "NO_PARAM_TEMP_TEST = %d\n", (int) NO_PARAM_TEMP_TEST); fprintf(ptr_asa_out, "NO_COST_TEMP_TEST = %d\n", (int) NO_COST_TEMP_TEST); fprintf(ptr_asa_out, "SELF_OPTIMIZE = %d\n", (int) SELF_OPTIMIZE); fprintf(ptr_asa_out, "ASA_TEST = %d\n", (int) ASA_TEST); fprintf(ptr_asa_out, "ASA_TEMPLATE = %d\n", (int) ASA_TEMPLATE); fprintf(ptr_asa_out, "OPTIONAL_DATA = %d\n\n", (int) OPTIONAL_DATA); fprintf(ptr_asa_out, "ASA_PRINT = %d\n", (int) ASA_PRINT); fprintf(ptr_asa_out, "ASA_OUT = %s\n", #if USER_ASA_OUT OPTIONS->asa_out_file); #else ASA_OUT); #endif fprintf(ptr_asa_out, "USER_ASA_OUT = %d\n", (int) USER_ASA_OUT); fprintf(ptr_asa_out, "ASA_PRINT_INTERMED = %d\n", (int) ASA_PRINT_INTERMED); fprintf(ptr_asa_out, "ASA_PRINT_MORE = %d\n\n", (int) ASA_PRINT_MORE); fprintf(ptr_asa_out, "OPTIONS->LIMIT_ACCEPTANCES = %d\n", OPTIONS->LIMIT_ACCEPTANCES); fprintf(ptr_asa_out, "OPTIONS->LIMIT_INVALID_GENERATED_STATES = %d\n", OPTIONS->LIMIT_INVALID_GENERATED_STATES); fprintf(ptr_asa_out, "OPTIONS->ACCEPTED_TO_GENERATED_RATIO = %g\n\n", OPTIONS->ACCEPTED_TO_GENERATED_RATIO); fprintf(ptr_asa_out, "OPTIONS->COST_PRECISION = %g\n", OPTIONS->COST_PRECISION); fprintf(ptr_asa_out, "OPTIONS->MAXIMUM_COST_REPEAT = %d\n", OPTIONS->MAXIMUM_COST_REPEAT); fprintf(ptr_asa_out, "OPTIONS->NUMBER_COST_SAMPLES = %d\n", OPTIONS->NUMBER_COST_SAMPLES); fprintf(ptr_asa_out, "OPTIONS->TEMPERATURE_RATIO_SCALE = %g\n", OPTIONS->TEMPERATURE_RATIO_SCALE); fprintf(ptr_asa_out, "OPTIONS->COST_PARAMETER_SCALE = %g\n", OPTIONS->COST_PARAMETER_SCALE); fprintf(ptr_asa_out, "OPTIONS->TEMPERATURE_ANNEAL_SCALE = %g\n", OPTIONS->TEMPERATURE_ANNEAL_SCALE); fprintf(ptr_asa_out, "OPTIONS->USER_INITIAL_COST_TEMP = %d\n\n", OPTIONS->USER_INITIAL_COST_TEMP); fprintf(ptr_asa_out, "OPTIONS->INCLUDE_INTEGER_PARAMETERS = %d\n", OPTIONS->INCLUDE_INTEGER_PARAMETERS); fprintf(ptr_asa_out, "OPTIONS->USER_INITIAL_PARAMETERS = %d\n", OPTIONS->USER_INITIAL_PARAMETERS); fprintf(ptr_asa_out, "OPTIONS->INITIAL_PARAMETER_TEMPERATURE = %g\n", OPTIONS->INITIAL_PARAMETER_TEMPERATURE); fprintf(ptr_asa_out, "OPTIONS->RATIO_TEMPERATURE_SCALES = %d\n", OPTIONS->RATIO_TEMPERATURE_SCALES); fprintf(ptr_asa_out, "OPTIONS->USER_INITIAL_PARAMETERS_TEMPS = %d\n\n", OPTIONS->USER_INITIAL_PARAMETERS_TEMPS); fprintf(ptr_asa_out, "OPTIONS->TESTING_FREQUENCY_MODULUS = %d\n", OPTIONS->TESTING_FREQUENCY_MODULUS); fprintf(ptr_asa_out, "OPTIONS->ACTIVATE_REANNEAL = %d\n", OPTIONS->ACTIVATE_REANNEAL); fprintf(ptr_asa_out, "OPTIONS->REANNEAL_RESCALE = %g\n", OPTIONS->REANNEAL_RESCALE); #if INT_LONG fprintf(ptr_asa_out, "OPTIONS->MAXIMUM_REANNEAL_INDEX = %ld\n", (LONG_INT) OPTIONS->MAXIMUM_REANNEAL_INDEX); #else fprintf(ptr_asa_out, "OPTIONS->MAXIMUM_REANNEAL_INDEX = %d\n\n", OPTIONS->MAXIMUM_REANNEAL_INDEX); #endif fprintf(ptr_asa_out, "OPTIONS->DELTA_X = %g\n", OPTIONS->DELTA_X); fprintf(ptr_asa_out, "OPTIONS->DELTA_PARAMETERS = %d\n", OPTIONS->DELTA_PARAMETERS); fprintf(ptr_asa_out, "OPTIONS->CURVATURE_0 = %d\n\n", OPTIONS->CURVATURE_0); fprintf(ptr_asa_out, "OPTIONS->QUENCH_PARAMETERS = %d\n", OPTIONS->QUENCH_PARAMETERS); fprintf(ptr_asa_out, "OPTIONS->QUENCH_COST = %d\n\n", OPTIONS->QUENCH_COST); fprintf(ptr_asa_out, "\n"); #if TIME_CALC /* print starting time */ print_asa_time("start_asa", ptr_asa_out); #endif fflush(ptr_asa_out); #endif /* set indices and counts to 0 */ *best_number_generated_saved = *number_generated = *recent_number_generated = *recent_number_acceptances = 0; *index_cost_acceptances = *best_number_accepted_saved = *number_accepted = *number_acceptances_saved = 0; index_cost_repeat = 0; /* do not calculate curvatures initially */ *curvature_flag = FALSE; /* allocate storage for all parameters */ if ((current_generated_state->parameter = (double *) calloc(*number_parameters, sizeof(double)) ) == NULL) exit(9); if ((last_saved_state->parameter = (double *) calloc(*number_parameters, sizeof(double)) ) == NULL) exit(9); if ((best_generated_state->parameter = (double *) calloc(*number_parameters, sizeof(double)) ) == NULL) exit(9); if ((initial_user_parameter_temp = (double *) calloc(*number_parameters, sizeof(double)) ) == NULL) exit(9); if ((index_parameter_generations = (LONG_INT *) calloc(*number_parameters, sizeof(LONG_INT)) ) == NULL) exit(9); /* set all delta_parameter[] */ if (OPTIONS->DELTA_PARAMETERS == TRUE) { delta_parameter = OPTIONS->user_delta_parameter; VFOR(index_v) delta_parameter[index_v] = OPTIONS->user_delta_parameter[index_v]; } else { if ((delta_parameter = (double *) calloc(*number_parameters, sizeof(double)) ) == NULL) exit(9); VFOR(index_v) delta_parameter[index_v] = OPTIONS->DELTA_X; } /* set all quench_param[] */ if (OPTIONS->QUENCH_PARAMETERS == TRUE) { quench_param = OPTIONS->user_quench_param_scale; VFOR(index_v) quench_param[index_v] = OPTIONS->user_quench_param_scale[index_v]; } else { if ((quench_param = (double *) calloc(*number_parameters, sizeof(double)) ) == NULL) exit(9); VFOR(index_v) quench_param[index_v] = ONE; } /* set quench_cost[] */ if (OPTIONS->QUENCH_COST == TRUE) { quench_cost = OPTIONS->user_quench_cost_scale; *quench_cost = OPTIONS->user_quench_cost_scale[0]; } else { if ((quench_cost = (double *) calloc(1, sizeof(double)) ) == NULL) exit(9); *quench_cost = ONE; } /* set all temperatures */ if ((current_user_parameter_temp = (double *) calloc(*number_parameters, sizeof(double)) ) == NULL) exit(9); if (OPTIONS->USER_INITIAL_PARAMETERS_TEMPS == TRUE) { VFOR(index_v) current_user_parameter_temp[index_v] = initial_user_parameter_temp[index_v] = OPTIONS->user_parameter_temperature[index_v]; } else { VFOR(index_v) current_user_parameter_temp[index_v] = initial_user_parameter_temp[index_v] = OPTIONS->INITIAL_PARAMETER_TEMPERATURE; } if (OPTIONS->RATIO_TEMPERATURE_SCALES == TRUE) { temperature_scale_parameters = OPTIONS->user_temperature_ratio; VFOR(index_v) temperature_scale_parameters[index_v] = OPTIONS->user_temperature_ratio[index_v]; } else { if ((temperature_scale_parameters = (double *) calloc(*number_parameters, sizeof(double)) ) == NULL) exit(9); } if (OPTIONS->USER_INITIAL_COST_TEMP == TRUE) *initial_cost_temperature = *current_cost_temperature = OPTIONS->user_cost_temperature[0]; /* set parameters to the initial parameter values */ VFOR(index_v) last_saved_state->parameter[index_v] = current_generated_state->parameter[index_v] = parameter_initial_final[index_v]; #if ASA_PRINT #if INT_ALLOC fprintf(ptr_asa_out, "*number_parameters = %d\n\n", *number_parameters); #else fprintf(ptr_asa_out, #if INT_LONG "*number_parameters = %ld\n\n", *number_parameters); #else "*number_parameters = %d\n\n", *number_parameters); #endif #endif /* print the min, max, current values, and types of parameters */ fprintf(ptr_asa_out, "index_v parameter_minimum parameter_maximum\ parameter_value parameter_type \n"); #if ASA_PRINT_INTERMED VFOR(index_v) #if INT_ALLOC fprintf(ptr_asa_out, " %-8d %-18.7g %-17.7g %-15.7g %-7d\n", #else #if INT_LONG fprintf(ptr_asa_out, " %-8ld %-18.7g %-17.7g %-15.7g %-7d\n", #else fprintf(ptr_asa_out, " %-8d %-18.7g %-17.7g %-15.7g %-7d\n", #endif #endif index_v, parameter_minimum[index_v], parameter_maximum[index_v], current_generated_state->parameter[index_v], parameter_type[index_v]); fprintf(ptr_asa_out, "\n\n"); #endif /* Print out user-defined OPTIONS */ if (OPTIONS->DELTA_PARAMETERS) VFOR(index_v) fprintf(ptr_asa_out, #if INT_ALLOC "OPTIONS->user_delta_parameter[%d] = %12.7g\n", #else #if INT_LONG "OPTIONS->user_delta_parameter[%ld] = %12.7g\n", #else "OPTIONS->user_delta_parameter[%d] = %12.7g\n", #endif #endif index_v, delta_parameter[index_v]); if (OPTIONS->QUENCH_PARAMETERS) VFOR(index_v) fprintf(ptr_asa_out, #if INT_ALLOC "OPTIONS->user_quench_param_scale[%d] = %12.7g\n", #else #if INT_LONG "OPTIONS->user_quench_param_scale[%ld] = %12.7g\n", #else "OPTIONS->user_quench_param_scale[%d] = %12.7g\n", #endif #endif index_v, quench_param[index_v]); if (OPTIONS->QUENCH_COST) fprintf(ptr_asa_out, "OPTIONS->user_quench_cost_scale = %12.7g\n", *quench_cost); if (OPTIONS->USER_INITIAL_PARAMETERS_TEMPS) VFOR(index_v) fprintf(ptr_asa_out, #if INT_ALLOC "OPTIONS->user_parameter_temperature[%d] = %12.7g\n", #else #if INT_LONG "OPTIONS->user_parameter_temperature[%ld] = %12.7g\n", #else "OPTIONS->user_parameter_temperature[%d] = %12.7g\n", #endif #endif index_v, initial_user_parameter_temp[index_v]); if (OPTIONS->RATIO_TEMPERATURE_SCALES) VFOR(index_v) fprintf(ptr_asa_out, #if INT_ALLOC "OPTIONS->user_temperature_ratio[%d] = %12.7g\n", #else #if INT_LONG "OPTIONS->user_temperature_ratio[%ld] = %12.7g\n", #else "OPTIONS->user_temperature_ratio[%d] = %12.7g\n", #endif #endif index_v, temperature_scale_parameters[index_v]); if (OPTIONS->USER_INITIAL_COST_TEMP) fprintf(ptr_asa_out, "OPTIONS->user_cost_temperature[0] = %12.7g\n", *initial_cost_temperature); fflush(ptr_asa_out); #endif /* set the initial index of parameter generations to 1 */ VFOR(index_v) index_parameter_generations[index_v] = 1; if (OPTIONS->USER_INITIAL_COST_TEMP == FALSE) { /* calculate the average cost over samplings of the cost function */ sampled_cost_sum = ZERO; for (index_cost_constraint = 0; index_cost_constraint < (OPTIONS->NUMBER_COST_SAMPLES); ++index_cost_constraint) { *number_invalid_generated_states = 0; repeated_invalid_states = 0; do { ++(*number_invalid_generated_states); ++repeated_invalid_states; *valid_state_generated_flag = TRUE; generate_new_state(user_random_generator, parameter_minimum, parameter_maximum, current_user_parameter_temp, number_parameters, parameter_type, current_generated_state, last_saved_state); 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 (repeated_invalid_states > OPTIONS->LIMIT_INVALID_GENERATED_STATES) { *exit_status = (int) TOO_MANY_INVALID_STATES; asa_exit(user_cost_function, &final_cost, parameter_initial_final, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, current_cost_temperature, initial_user_parameter_temp, current_user_parameter_temp, temperature_scale_parameters, accepted_to_generated_ratio, number_parameters, parameter_type, valid_state_generated_flag, exit_status, index_exit_v, number_accepted, best_number_accepted_saved, index_cost_acceptances, number_generated, number_invalid_generated_states, index_parameter_generations, best_number_generated_saved, current_generated_state, last_saved_state, best_generated_state, ptr_asa_out, OPTIONS); free(curvature_flag); free(quench_param_flag); free(quench_cost_flag); free(maximum_tangent); free(accepted_to_generated_ratio); free(temperature_scale_cost); free(current_cost_temperature); free(initial_cost_temperature); free(number_generated); free(best_number_generated_saved); free(recent_number_generated); free(number_accepted); free(recent_number_acceptances); free(index_cost_acceptances); free(number_acceptances_saved); free(best_number_accepted_saved); free(number_invalid_generated_states); free(current_generated_state->parameter); free(last_saved_state->parameter); free(best_generated_state->parameter); free(current_generated_state); free(last_saved_state); free(best_generated_state); free(initial_user_parameter_temp); free(index_exit_v); free(index_parameter_generations); if (OPTIONS->DELTA_PARAMETERS == FALSE) free(delta_parameter); if (OPTIONS->QUENCH_PARAMETERS == FALSE) free(quench_param); if (OPTIONS->QUENCH_COST == FALSE) free(quench_cost); free(current_user_parameter_temp); if (OPTIONS->RATIO_TEMPERATURE_SCALES == FALSE) free(temperature_scale_parameters); if (recursive_asa_open == 0) asa_open = FALSE; return (final_cost); } } while (*valid_state_generated_flag == FALSE); --(*number_invalid_generated_states); sampled_cost_sum += fabs(current_generated_state->cost); } /* set all cost temperatures to the average cost */ *initial_cost_temperature = *current_cost_temperature = sampled_cost_sum / (OPTIONS->NUMBER_COST_SAMPLES); } /* set all parameters to the initial parameter values */ VFOR(index_v) best_generated_state->parameter[index_v] = last_saved_state->parameter[index_v] = current_generated_state->parameter[index_v] = parameter_initial_final[index_v]; /* if using user's initial parameters */ if (OPTIONS->USER_INITIAL_PARAMETERS == TRUE) { *valid_state_generated_flag = TRUE; 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 ASA_PRINT if (*valid_state_generated_flag == FALSE) fprintf(ptr_asa_out, "user's initial parameters generated \ FALSE *valid_state_generated_flag\n"); #endif } else { /* let asa generate valid initial parameters */ repeated_invalid_states = 0; do { ++(*number_invalid_generated_states); ++repeated_invalid_states; *valid_state_generated_flag = TRUE; generate_new_state(user_random_generator, parameter_minimum, parameter_maximum, current_user_parameter_temp, number_parameters, parameter_type, current_generated_state, last_saved_state); 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 (repeated_invalid_states > OPTIONS->LIMIT_INVALID_GENERATED_STATES) { *exit_status = (int) TOO_MANY_INVALID_STATES; asa_exit(user_cost_function, &final_cost, parameter_initial_final, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, current_cost_temperature, initial_user_parameter_temp, current_user_parameter_temp, temperature_scale_parameters, accepted_to_generated_ratio, number_parameters, parameter_type, valid_state_generated_flag, exit_status, index_exit_v, number_accepted, best_number_accepted_saved, index_cost_acceptances, number_generated, number_invalid_generated_states, index_parameter_generations, best_number_generated_saved, current_generated_state, last_saved_state, best_generated_state, ptr_asa_out, OPTIONS); free(curvature_flag); free(quench_param_flag); free(quench_cost_flag); free(maximum_tangent); free(accepted_to_generated_ratio); free(temperature_scale_cost); free(current_cost_temperature); free(initial_cost_temperature); free(number_generated); free(best_number_generated_saved); free(recent_number_generated); free(number_accepted); free(recent_number_acceptances); free(index_cost_acceptances); free(number_acceptances_saved); free(best_number_accepted_saved); free(number_invalid_generated_states); free(current_generated_state->parameter); free(last_saved_state->parameter); free(best_generated_state->parameter); free(current_generated_state); free(last_saved_state); free(best_generated_state); free(initial_user_parameter_temp); free(index_exit_v); free(index_parameter_generations); if (OPTIONS->DELTA_PARAMETERS == FALSE) free(delta_parameter); if (OPTIONS->QUENCH_PARAMETERS == FALSE) free(quench_param); if (OPTIONS->QUENCH_COST == FALSE) free(quench_cost); free(current_user_parameter_temp); if (OPTIONS->RATIO_TEMPERATURE_SCALES == FALSE) free(temperature_scale_parameters); if (recursive_asa_open == 0) asa_open = FALSE; return (final_cost); } } 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) { /* ignore parameters that have too small a range */ if (PARAMETER_RANGE_TOO_SMALL(index_v)) continue; 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; tmp_var_db1 = -log((OPTIONS->TEMPERATURE_RATIO_SCALE)); /* compute temperature scales */ tmp_var_db2 = log(OPTIONS->TEMPERATURE_ANNEAL_SCALE); temperature_scale = tmp_var_db1 * exp(-tmp_var_db2 / (double) *number_parameters); if (OPTIONS->RATIO_TEMPERATURE_SCALES) { VFOR(index_v) temperature_scale_parameters[index_v] = tmp_var_db1 * exp(-(tmp_var_db2 * quench_param[index_v]) / (double) *number_parameters) * OPTIONS->user_temperature_ratio[index_v]; } else { VFOR(index_v) temperature_scale_parameters[index_v] = tmp_var_db1 * exp(-(tmp_var_db2 * quench_param[index_v]) / (double) *number_parameters); } *temperature_scale_cost = tmp_var_db1 * exp(-(tmp_var_db2 * quench_cost[0]) / (double) *number_parameters) * OPTIONS->COST_PARAMETER_SCALE; /* do not calculate curvatures initially */ *curvature_flag = FALSE; #if ASA_PRINT fprintf(ptr_asa_out, "temperature_scale = %12.7g\n", temperature_scale); if (OPTIONS->RATIO_TEMPERATURE_SCALES) { #if ASA_PRINT_INTERMED VFOR(index_v) { fprintf(ptr_asa_out, #if INT_ALLOC "temperature_scale_parameters[%d] = %12.7g\n", #else #if INT_LONG "temperature_scale_parameters[%ld] = %12.7g\n", #else "temperature_scale_parameters[%d] = %12.7g\n", #endif #endif index_v, temperature_scale_parameters[index_v]); } #endif } else { fprintf(ptr_asa_out, "temperature_scale_parameters[0] = %12.7g\n", temperature_scale_parameters[0]); } fprintf(ptr_asa_out, "*temperature_scale_cost = %12.7g\n", *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 /* 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; 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; } /* MAIN ANNEALING LOOP */ while (*number_accepted <= ((LONG_INT) OPTIONS->LIMIT_ACCEPTANCES)) { /* 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] * pow((double) index_parameter_generations[index_v], quench_param[index_v] / (double) *number_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] * 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 = (int) P_TEMP_TOO_SMALL; *index_exit_v = index_v; /* Note that this abridged exit can hold and print old values of some variables, such as *current_cost_temperature. */ asa_exit(user_cost_function, &final_cost, parameter_initial_final, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, current_cost_temperature, initial_user_parameter_temp, current_user_parameter_temp, temperature_scale_parameters, accepted_to_generated_ratio, number_parameters, parameter_type, valid_state_generated_flag, exit_status, index_exit_v, number_accepted, best_number_accepted_saved, index_cost_acceptances, number_generated, number_invalid_generated_states, index_parameter_generations, best_number_generated_saved, current_generated_state, last_saved_state, best_generated_state, ptr_asa_out, OPTIONS); free(curvature_flag); free(quench_param_flag); free(quench_cost_flag); free(maximum_tangent); free(accepted_to_generated_ratio); free(temperature_scale_cost); free(current_cost_temperature); free(initial_cost_temperature); free(number_generated); free(best_number_generated_saved); free(recent_number_generated); free(number_accepted); free(recent_number_acceptances); free(index_cost_acceptances); free(number_acceptances_saved); free(best_number_accepted_saved); free(number_invalid_generated_states); free(current_generated_state->parameter); free(last_saved_state->parameter); free(best_generated_state->parameter); free(current_generated_state); free(last_saved_state); free(best_generated_state); free(initial_user_parameter_temp); free(index_exit_v); free(index_parameter_generations); if (OPTIONS->DELTA_PARAMETERS == FALSE) free(delta_parameter); if (OPTIONS->QUENCH_PARAMETERS == FALSE) free(quench_param); if (OPTIONS->QUENCH_COST == FALSE) free(quench_cost); free(current_user_parameter_temp); if (OPTIONS->RATIO_TEMPERATURE_SCALES == FALSE) free(temperature_scale_parameters); if (recursive_asa_open == 0) asa_open = FALSE; return (final_cost); } #endif } /* calculate new cost temperature */ log_new_temperature_ratio = -*temperature_scale_cost * pow((double) *index_cost_acceptances, (*quench_cost) / (double) *number_parameters); log_new_temperature_ratio = EXPONENT_CHECK(log_new_temperature_ratio); *current_cost_temperature = *initial_cost_temperature * exp(log_new_temperature_ratio); #if NO_COST_TEMP_TEST if (*current_cost_temperature < (double) EPS_DOUBLE) *current_cost_temperature = (double) EPS_DOUBLE; #else /* check for too small a cost temperature */ if (*current_cost_temperature < (double) EPS_DOUBLE) { *exit_status = (int) C_TEMP_TOO_SMALL; /* Note that this abridged exit can hold and print old values of some variables. */ asa_exit(user_cost_function, &final_cost, parameter_initial_final, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, current_cost_temperature, initial_user_parameter_temp, current_user_parameter_temp, temperature_scale_parameters, accepted_to_generated_ratio, number_parameters, parameter_type, valid_state_generated_flag, exit_status, index_exit_v, number_accepted, best_number_accepted_saved, index_cost_acceptances, number_generated, number_invalid_generated_states, index_parameter_generations, best_number_generated_saved, current_generated_state, last_saved_state, best_generated_state, ptr_asa_out, OPTIONS); free(curvature_flag); free(quench_param_flag); free(quench_cost_flag); free(maximum_tangent); free(accepted_to_generated_ratio); free(temperature_scale_cost); free(current_cost_temperature); free(initial_cost_temperature); free(number_generated); free(best_number_generated_saved); free(recent_number_generated); free(number_accepted); free(recent_number_acceptances); free(index_cost_acceptances); free(number_acceptances_saved); free(best_number_accepted_saved); free(number_invalid_generated_states); free(current_generated_state->parameter); free(last_saved_state->parameter); free(best_generated_state->parameter); free(current_generated_state); free(last_saved_state); free(best_generated_state); free(initial_user_parameter_temp); free(index_exit_v); free(index_parameter_generations); if (OPTIONS->DELTA_PARAMETERS == FALSE) free(delta_parameter); if (OPTIONS->QUENCH_PARAMETERS == FALSE) free(quench_param); if (OPTIONS->QUENCH_COST == FALSE) free(quench_cost); free(current_user_parameter_temp); if (OPTIONS->RATIO_TEMPERATURE_SCALES == FALSE) free(temperature_scale_parameters); if (recursive_asa_open == 0) asa_open = FALSE; return (final_cost); } #endif /* GENERATE NEW PARAMETERS */ /* generate a new valid set of parameters */ repeated_invalid_states = 0; do { ++(*number_invalid_generated_states); ++repeated_invalid_states; *valid_state_generated_flag = TRUE; generate_new_state(user_random_generator, parameter_minimum, parameter_maximum, current_user_parameter_temp, number_parameters, parameter_type, current_generated_state, last_saved_state); 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 (repeated_invalid_states > OPTIONS->LIMIT_INVALID_GENERATED_STATES) { *exit_status = (int) TOO_MANY_INVALID_STATES; asa_exit(user_cost_function, &final_cost, parameter_initial_final, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, current_cost_temperature, initial_user_parameter_temp, current_user_parameter_temp, temperature_scale_parameters, accepted_to_generated_ratio, number_parameters, parameter_type, valid_state_generated_flag, exit_status, index_exit_v, number_accepted, best_number_accepted_saved, index_cost_acceptances, number_generated, number_invalid_generated_states, index_parameter_generations, best_number_generated_saved, current_generated_state, last_saved_state, best_generated_state, ptr_asa_out, OPTIONS); free(curvature_flag); free(quench_param_flag); free(quench_cost_flag); free(maximum_tangent); free(accepted_to_generated_ratio); free(temperature_scale_cost); free(current_cost_temperature); free(initial_cost_temperature); free(number_generated); free(best_number_generated_saved); free(recent_number_generated); free(number_accepted); free(recent_number_acceptances); free(index_cost_acceptances); free(number_acceptances_saved); free(best_number_accepted_saved); free(number_invalid_generated_states); free(current_generated_state->parameter); free(last_saved_state->parameter); free(best_generated_state->parameter); free(current_generated_state); free(last_saved_state); free(best_generated_state); free(initial_user_parameter_temp); free(index_exit_v); free(index_parameter_generations); if (OPTIONS->DELTA_PARAMETERS == FALSE) free(delta_parameter); if (OPTIONS->QUENCH_PARAMETERS == FALSE) free(quench_param); if (OPTIONS->QUENCH_COST == FALSE) free(quench_cost); free(current_user_parameter_temp); if (OPTIONS->RATIO_TEMPERATURE_SCALES == FALSE) free(temperature_scale_parameters); if (recursive_asa_open == 0) asa_open = FALSE; return (final_cost); } } while (*valid_state_generated_flag == FALSE); --(*number_invalid_generated_states); /* ACCEPT/REJECT NEW PARAMETERS */ /* decide to accept/reject the new state */ accept_new_state(user_random_generator, parameter_minimum, parameter_maximum, current_cost_temperature, number_parameters, recent_number_acceptances, number_accepted, index_cost_acceptances, number_acceptances_saved, recent_number_generated, number_generated, index_parameter_generations, current_generated_state, last_saved_state); /* calculate the ratio of acceptances to generated states */ *accepted_to_generated_ratio = (double) (*recent_number_acceptances + 1) / (double) (*recent_number_generated + 1); /* CHECK FOR NEW MINIMUM */ if (current_generated_state->cost < best_generated_state->cost) { /* NEW MINIMUM FOUND */ /* reset the recent acceptances and generated counts */ *recent_number_acceptances = *recent_number_generated = 0; *best_number_generated_saved = *number_generated; *best_number_accepted_saved = *number_accepted; index_cost_repeat = 0; /* copy the current state into the best_generated state */ best_generated_state->cost = current_generated_state->cost; VFOR(index_v) { /* ignore parameters that have too small a range */ if (PARAMETER_RANGE_TOO_SMALL(index_v)) continue; best_generated_state->parameter[index_v] = current_generated_state->parameter[index_v]; } /* printout the new minimum state and value */ #if ASA_PRINT #if INT_LONG fprintf(ptr_asa_out, "best...->cost=%-12.7g \ *number_accepted=%ld *number_generated=%ld\n", best_generated_state->cost, *number_accepted, *number_generated); #else fprintf(ptr_asa_out, "best...->cost=%-12.7g \ *number_accepted=%d *number_generated=%d\n", best_generated_state->cost, *number_accepted, *number_generated); #endif #if ASA_PRINT_MORE 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 fflush(ptr_asa_out); #endif } /* PERIODIC TESTING/REANNEALING/PRINTING SECTION */ temp_var_int = (int) (*number_accepted % ((LONG_INT) OPTIONS->TESTING_FREQUENCY_MODULUS)); if ((temp_var_int == 0 && *number_acceptances_saved == *number_accepted) || *accepted_to_generated_ratio < (OPTIONS->ACCEPTED_TO_GENERATED_RATIO)) { if (*accepted_to_generated_ratio < (OPTIONS->ACCEPTED_TO_GENERATED_RATIO)) *recent_number_acceptances = *recent_number_generated = 0; /* if best.cost repeats OPTIONS->MAXIMUM_COST_REPEAT then exit */ if (fabs(last_saved_state->cost - best_generated_state->cost) < OPTIONS->COST_PRECISION) { ++index_cost_repeat; if (index_cost_repeat == (OPTIONS->MAXIMUM_COST_REPEAT)) { *exit_status = (int) COST_REPEATING; /* Note that this abridged exit can hold and print old values of some variables. */ asa_exit(user_cost_function, &final_cost, parameter_initial_final, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, current_cost_temperature, initial_user_parameter_temp, current_user_parameter_temp, temperature_scale_parameters, accepted_to_generated_ratio, number_parameters, parameter_type, valid_state_generated_flag, exit_status, index_exit_v, number_accepted, best_number_accepted_saved, index_cost_acceptances, number_generated, number_invalid_generated_states, index_parameter_generations, best_number_generated_saved, current_generated_state, last_saved_state, best_generated_state, ptr_asa_out, OPTIONS); free(curvature_flag); free(quench_param_flag); free(quench_cost_flag); free(maximum_tangent); free(accepted_to_generated_ratio); free(temperature_scale_cost); free(current_cost_temperature); free(initial_cost_temperature); free(number_generated); free(best_number_generated_saved); free(recent_number_generated); free(number_accepted); free(recent_number_acceptances); free(index_cost_acceptances); free(number_acceptances_saved); free(best_number_accepted_saved); free(number_invalid_generated_states); free(current_generated_state->parameter); free(last_saved_state->parameter); free(best_generated_state->parameter); free(current_generated_state); free(last_saved_state); free(best_generated_state); free(initial_user_parameter_temp); free(index_exit_v); free(index_parameter_generations); if (OPTIONS->DELTA_PARAMETERS == FALSE) free(delta_parameter); if (OPTIONS->QUENCH_PARAMETERS == FALSE) free(quench_param); if (OPTIONS->QUENCH_COST == FALSE) free(quench_cost); free(current_user_parameter_temp); if (OPTIONS->RATIO_TEMPERATURE_SCALES == FALSE) free(temperature_scale_parameters); if (recursive_asa_open == 0) asa_open = FALSE; return (final_cost); } } else { index_cost_repeat = 0; } /* set to FALSE to skip reannealing */ if (OPTIONS->ACTIVATE_REANNEAL == TRUE) { /* calculate tangents, not curvatures, to reanneal */ *curvature_flag = FALSE; cost_derivatives(user_cost_function, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, number_parameters, parameter_type, exit_status, curvature_flag, valid_state_generated_flag, number_invalid_generated_states, current_generated_state, best_generated_state, ptr_asa_out, OPTIONS); reanneal(parameter_minimum, parameter_maximum, tangents, maximum_tangent, current_cost_temperature, initial_cost_temperature, temperature_scale_cost, current_user_parameter_temp, initial_user_parameter_temp, temperature_scale_parameters, quench_param, quench_cost, number_parameters, parameter_type, index_cost_acceptances, index_parameter_generations, last_saved_state, best_generated_state, OPTIONS); } #if ASA_PRINT_INTERMED #if ASA_PRINT 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); fprintf(ptr_asa_out, "\n"); fflush(ptr_asa_out); #endif #endif } } /* FINISHED ANNEALING and MINIMIZATION */ *exit_status = (int) NORMAL_EXIT; asa_exit(user_cost_function, &final_cost, parameter_initial_final, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, current_cost_temperature, initial_user_parameter_temp, current_user_parameter_temp, temperature_scale_parameters, accepted_to_generated_ratio, number_parameters, parameter_type, valid_state_generated_flag, exit_status, index_exit_v, number_accepted, best_number_accepted_saved, index_cost_acceptances, number_generated, number_invalid_generated_states, index_parameter_generations, best_number_generated_saved, current_generated_state, last_saved_state, best_generated_state, ptr_asa_out, OPTIONS); free(curvature_flag); free(quench_param_flag); free(quench_cost_flag); free(maximum_tangent); free(accepted_to_generated_ratio); free(temperature_scale_cost); free(current_cost_temperature); free(initial_cost_temperature); free(number_generated); free(best_number_generated_saved); free(recent_number_generated); free(number_accepted); free(recent_number_acceptances); free(index_cost_acceptances); free(number_acceptances_saved); free(best_number_accepted_saved); free(number_invalid_generated_states); free(current_generated_state->parameter); free(last_saved_state->parameter); free(best_generated_state->parameter); free(current_generated_state); free(last_saved_state); free(best_generated_state); free(initial_user_parameter_temp); free(index_exit_v); free(index_parameter_generations); if (OPTIONS->DELTA_PARAMETERS == FALSE) free(delta_parameter); if (OPTIONS->QUENCH_PARAMETERS == FALSE) free(quench_param); if (OPTIONS->QUENCH_COST == FALSE) free(quench_cost); free(current_user_parameter_temp); if (OPTIONS->RATIO_TEMPERATURE_SCALES == FALSE) free(temperature_scale_parameters); if (recursive_asa_open == 0) asa_open = FALSE; return (final_cost); } /*********************************************************************** * asa_exit * This procedures copies the best parameters and cost into * final_cost and parameter_initial_final ***********************************************************************/ #if HAVE_ANSI void asa_exit(double (*user_cost_function) (), double *final_cost, double *parameter_initial_final, double *parameter_minimum, double *parameter_maximum, double *tangents, double *curvature, double *maximum_tangent, double *delta_parameter, double *current_cost_temperature, double *initial_user_parameter_temp, double *current_user_parameter_temp, double *temperature_scale_parameters, double *accepted_to_generated_ratio, ALLOC_INT * number_parameters, int *parameter_type, int *valid_state_generated_flag, int *exit_status, ALLOC_INT * index_exit_v, LONG_INT * number_accepted, LONG_INT * best_number_accepted_saved, LONG_INT * index_cost_acceptances, LONG_INT * number_generated, LONG_INT * number_invalid_generated_states, LONG_INT * index_parameter_generations, LONG_INT * best_number_generated_saved, STATE * current_generated_state, STATE * last_saved_state, STATE * best_generated_state, FILE * ptr_asa_out, DEFINES * OPTIONS) #else void asa_exit(user_cost_function, final_cost, parameter_initial_final, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, current_cost_temperature, initial_user_parameter_temp, current_user_parameter_temp, temperature_scale_parameters, accepted_to_generated_ratio, number_parameters, parameter_type, valid_state_generated_flag, exit_status, index_exit_v, number_accepted, best_number_accepted_saved, index_cost_acceptances, number_generated, number_invalid_generated_states, index_parameter_generations, best_number_generated_saved, current_generated_state, last_saved_state, best_generated_state, ptr_asa_out, OPTIONS) double (*user_cost_function) (); double *final_cost; double *parameter_initial_final; double *parameter_minimum; double *parameter_maximum; double *tangents; double *curvature; double *maximum_tangent; double *delta_parameter; double *current_cost_temperature; double *initial_user_parameter_temp; double *current_user_parameter_temp; double *temperature_scale_parameters; double *accepted_to_generated_ratio; ALLOC_INT *number_parameters; int *parameter_type; int *valid_state_generated_flag; int *exit_status; ALLOC_INT *index_exit_v; LONG_INT *number_accepted; LONG_INT *best_number_accepted_saved; LONG_INT *index_cost_acceptances; LONG_INT *number_generated; LONG_INT *number_invalid_generated_states; LONG_INT *index_parameter_generations; LONG_INT *best_number_generated_saved; STATE *current_generated_state; STATE *last_saved_state; STATE *best_generated_state; FILE *ptr_asa_out; DEFINES *OPTIONS; #endif { ALLOC_INT index_v; /* iteration index */ int *curvature_flag; if ((curvature_flag = (int *) calloc(1, sizeof(int))) == NULL) exit(9); if (*exit_status != TOO_MANY_INVALID_STATES) { /* calculate curvatures and tangents at best point */ *curvature_flag = TRUE; cost_derivatives(user_cost_function, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, number_parameters, parameter_type, exit_status, curvature_flag, valid_state_generated_flag, number_invalid_generated_states, current_generated_state, best_generated_state, ptr_asa_out, OPTIONS); } /* return final function minimum and associated parameters */ *final_cost = best_generated_state->cost; VFOR(index_v) { parameter_initial_final[index_v] = best_generated_state->parameter[index_v]; } #if ASA_PRINT 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); switch (*exit_status) { case NORMAL_EXIT: fprintf(ptr_asa_out, "\n\n NORMAL_EXIT exit_status = %d\n", *exit_status); break; case P_TEMP_TOO_SMALL: fprintf(ptr_asa_out, "\n\n P_TEMP_TOO_SMALL exit_status = %d\n", *exit_status); fprintf(ptr_asa_out, #if INT_ALLOC "current_user_parameter_temp[%d] too small = %12.7g\n", #else #if INT_LONG "current_user_parameter_temp[%ld] too small = %12.7g\n", #else "current_user_parameter_temp[%d] too small = %12.7g\n", #endif #endif *index_exit_v, current_user_parameter_temp[*index_exit_v]); break; case C_TEMP_TOO_SMALL: fprintf(ptr_asa_out, "\n\n C_TEMP_TOO_SMALL exit_status = %d\n", *exit_status); fprintf(ptr_asa_out, "*current_cost_temperature too small = %12.7g\n", *current_cost_temperature); break; case COST_REPEATING: fprintf(ptr_asa_out, "\n\n COST_REPEATING exit_status = %d\n", *exit_status); break; case TOO_MANY_INVALID_STATES: fprintf(ptr_asa_out, "\n\n TOO_MANY_INVALID_STATES exit_status = %d\n", *exit_status); break; default: fprintf(ptr_asa_out, "\n\n ERR: no exit code available = %d\n", *exit_status); } fprintf(ptr_asa_out, "final_cost = best_generated_state->cost = %-12.7g\n", *final_cost); #if INT_LONG fprintf(ptr_asa_out, "*number_accepted at best_generated_state->cost = %ld\n", *best_number_accepted_saved); fprintf(ptr_asa_out, "*number_generated at best_generated_state->cost = %ld\n", *best_number_generated_saved); #else fprintf(ptr_asa_out, "*number_accepted at best_generated_state->cost = %d\n", *best_number_accepted_saved); fprintf(ptr_asa_out, "*number_generated at best_generated_state->cost = %d\n", *best_number_generated_saved); #endif #endif #if ASA_TEMPLATE #if OPTIONAL_DATA if (OPTIONS->asa_data[0] > (double) MIN_DOUBLE) OPTIONS->asa_data[1] = (double) (*best_number_generated_saved); #endif #endif /* return unused space */ free(curvature_flag); #if ASA_PRINT #if TIME_CALC /* print ending time */ print_asa_time("asa_end", ptr_asa_out); #endif fprintf(ptr_asa_out, "\n\n\n"); fflush(ptr_asa_out); fclose(ptr_asa_out); #endif } /*********************************************************************** * generate_new_state * Generates a valid new state from the old state ***********************************************************************/ /* generate a new state from the old state that lies between the min and max parameter values */ #if HAVE_ANSI void generate_new_state(double (*user_random_generator) (), double *parameter_minimum, double *parameter_maximum, double *current_user_parameter_temp, ALLOC_INT * number_parameters, int *parameter_type, STATE * current_generated_state, STATE * last_saved_state) #else void generate_new_state(user_random_generator, parameter_minimum, parameter_maximum, current_user_parameter_temp, number_parameters, parameter_type, current_generated_state, last_saved_state) double (*user_random_generator) (); double *parameter_minimum; double *parameter_maximum; double *current_user_parameter_temp; ALLOC_INT *number_parameters; int *parameter_type; STATE *current_generated_state; STATE *last_saved_state; #endif { ALLOC_INT index_v; double x; double parameter_v, min_parameter_v, max_parameter_v, temperature_v, parameter_range_v; /* generate a new value for each parameter */ VFOR(index_v) { min_parameter_v = parameter_minimum[index_v]; max_parameter_v = parameter_maximum[index_v]; parameter_range_v = max_parameter_v - min_parameter_v; /* ignore parameters that have too small a range */ if (fabs(parameter_range_v) < (double) EPS_DOUBLE) continue; temperature_v = current_user_parameter_temp[index_v]; parameter_v = last_saved_state->parameter[index_v]; /* Handle discrete integer parameters. */ if (INTEGER_PARAMETER(index_v)) { min_parameter_v -= HALF; max_parameter_v += HALF; parameter_range_v = max_parameter_v - min_parameter_v; } /* generate a new state x until within the parameter bounds */ for (;;) { x = parameter_v + generate_asa_state(user_random_generator, &temperature_v) * parameter_range_v; /* exit the loop if within its valid parameter range */ if (x <= max_parameter_v - (double) EPS_DOUBLE && x >= min_parameter_v + (double) EPS_DOUBLE) break; } /* Handle discrete integer parameters. You might have to check rounding on your machine. */ if (INTEGER_PARAMETER(index_v)) { if (x < min_parameter_v + HALF) x = min_parameter_v + HALF + (double) EPS_DOUBLE; if (x > max_parameter_v - HALF) x = max_parameter_v - HALF + (double) EPS_DOUBLE; if (x + HALF > ZERO) { x = (int) (x + HALF); } else { x = (int) (x - HALF); } if (x > parameter_maximum[index_v]) x = parameter_maximum[index_v]; if (x < parameter_minimum[index_v]) x = parameter_minimum[index_v]; } /* save the newly generated value */ current_generated_state->parameter[index_v] = x; } } /*********************************************************************** * generate_asa_state * This function generates a single value according to the * ASA generating function and the passed temperature ***********************************************************************/ #if HAVE_ANSI double generate_asa_state(double (*user_random_generator) (), double *temp) #else double generate_asa_state(user_random_generator, temp) double (*user_random_generator) (); double *temp; #endif { double x, y, z; x = (*user_random_generator) (); y = x < HALF ? -ONE : ONE; z = y * *temp * (pow((ONE + ONE / *temp), fabs(TWO * x - ONE)) - ONE); return (z); } /*********************************************************************** * accept_new_state * This procedure accepts or rejects a newly generated state, * depending on whether the difference between new and old * cost functions passes a statistical test. If accepted, * the current state is updated. ***********************************************************************/ /* determine whether to accept or reject the new state */ #if HAVE_ANSI void accept_new_state(double (*user_random_generator) (), double *parameter_minimum, double *parameter_maximum, double *current_cost_temperature, ALLOC_INT * number_parameters, LONG_INT * recent_number_acceptances, LONG_INT * number_accepted, LONG_INT * index_cost_acceptances, LONG_INT * number_acceptances_saved, LONG_INT * recent_number_generated, LONG_INT * number_generated, LONG_INT * index_parameter_generations, STATE * current_generated_state, STATE * last_saved_state) #else void accept_new_state(user_random_generator, parameter_minimum, parameter_maximum, current_cost_temperature, number_parameters, recent_number_acceptances, number_accepted, index_cost_acceptances, number_acceptances_saved, recent_number_generated, number_generated, index_parameter_generations, current_generated_state, last_saved_state) double (*user_random_generator) (); double *parameter_minimum; double *parameter_maximum; double *current_cost_temperature; ALLOC_INT *number_parameters; LONG_INT *recent_number_acceptances; LONG_INT *number_accepted; LONG_INT *index_cost_acceptances; LONG_INT *number_acceptances_saved; LONG_INT *recent_number_generated; LONG_INT *number_generated; LONG_INT *index_parameter_generations; STATE *current_generated_state; STATE *last_saved_state; #endif { double delta_cost; ALLOC_INT index_v; /* update accepted and generated count */ ++*number_acceptances_saved; ++*recent_number_generated; ++*number_generated; /* increment the parameter index generation for each parameter */ VFOR(index_v) { /* ignore parameters with too small a range */ if (PARAMETER_RANGE_TOO_SMALL(index_v)) continue; ++(index_parameter_generations[index_v]); } /* effective cost function for testing acceptance criteria, calculate the cost difference and divide by the temperature */ delta_cost = (current_generated_state->cost - last_saved_state->cost) / (*current_cost_temperature + EPS_DOUBLE); /* decide to accept/reject the new state */ if (exp(EXPONENT_CHECK(-delta_cost)) > (*user_random_generator) ()) { /* copy the current generated parameters to the last saved state */ last_saved_state->cost = current_generated_state->cost; VFOR(index_v) { /* ignore parameters with too small a range */ if (PARAMETER_RANGE_TOO_SMALL(index_v)) continue; last_saved_state->parameter[index_v] = current_generated_state->parameter[index_v]; } /* update acceptance counts */ ++*recent_number_acceptances; ++*number_accepted; ++*index_cost_acceptances; *number_acceptances_saved = *number_accepted; } } /*********************************************************************** * reanneal * Readjust temperatures of generating and acceptance functions ***********************************************************************/ /* REANNEALING - readjust the temperatures */ #if HAVE_ANSI void reanneal(double *parameter_minimum, double *parameter_maximum, double *tangents, double *maximum_tangent, double *current_cost_temperature, double *initial_cost_temperature, double *temperature_scale_cost, double *current_user_parameter_temp, double *initial_user_parameter_temp, double *temperature_scale_parameters, double *quench_param, double *quench_cost, ALLOC_INT * number_parameters, int *parameter_type, LONG_INT * index_cost_acceptances, LONG_INT * index_parameter_generations, STATE * last_saved_state, STATE * best_generated_state, DEFINES * OPTIONS) #else void reanneal(parameter_minimum, parameter_maximum, tangents, maximum_tangent, current_cost_temperature, initial_cost_temperature, temperature_scale_cost, current_user_parameter_temp, initial_user_parameter_temp, temperature_scale_parameters, quench_param, quench_cost, number_parameters, parameter_type, index_cost_acceptances, index_parameter_generations, last_saved_state, best_generated_state, OPTIONS) double *parameter_minimum; double *parameter_maximum; double *tangents; double *maximum_tangent; double *current_cost_temperature; double *initial_cost_temperature; double *temperature_scale_cost; double *current_user_parameter_temp; double *initial_user_parameter_temp; double *temperature_scale_parameters; double *quench_param; double *quench_cost; ALLOC_INT *number_parameters; int *parameter_type; LONG_INT *index_cost_acceptances; LONG_INT *index_parameter_generations; STATE *last_saved_state; STATE *best_generated_state; DEFINES *OPTIONS; #endif { ALLOC_INT index_v; double tmp_var_db3; double new_temperature; /* the temperature */ double log_new_temperature_ratio; double log_init_cur_temp_ratio; double temperature_rescale_power; double temp_dbl; VFOR(index_v) { if (NO_REANNEAL(index_v)) continue; /* use the temp double to prevent overflow */ temp_dbl = (double) index_parameter_generations[index_v]; /* skip parameters with too small range or integer parameters */ if (OPTIONS->INCLUDE_INTEGER_PARAMETERS == TRUE) { if (PARAMETER_RANGE_TOO_SMALL(index_v)) continue; } else { if (PARAMETER_RANGE_TOO_SMALL(index_v) || INTEGER_PARAMETER(index_v)) continue; } /* ignore parameters with too small tangents */ if (fabs(tangents[index_v]) < (double) EPS_DOUBLE) continue; /* reset the index of parameter generations appropriately */ new_temperature = fabs(REANNEAL_FUNCTION(current_user_parameter_temp[index_v], tangents[index_v], *maximum_tangent)); if (new_temperature < initial_user_parameter_temp[index_v]) { log_init_cur_temp_ratio = fabs(log(((double) EPS_DOUBLE + initial_user_parameter_temp[index_v]) / ((double) EPS_DOUBLE + new_temperature))); temp_dbl = (double) EPS_DOUBLE + pow(log_init_cur_temp_ratio / temperature_scale_parameters[index_v], (double) *number_parameters / quench_param[index_v]); } else { temp_dbl = 1.0; } /* Reset index_parameter_generations if index reset too large, and also reset the initial_user_parameter_temp, to achieve the same new temperature. */ while (temp_dbl > ((double) OPTIONS->MAXIMUM_REANNEAL_INDEX)) { log_new_temperature_ratio = -temperature_scale_parameters[index_v] * pow(temp_dbl, quench_param[index_v] / (double) *number_parameters); log_new_temperature_ratio = EXPONENT_CHECK(log_new_temperature_ratio); new_temperature = initial_user_parameter_temp[index_v] * exp(log_new_temperature_ratio); temp_dbl /= OPTIONS->REANNEAL_RESCALE; temperature_rescale_power = ONE / pow(OPTIONS->REANNEAL_RESCALE, quench_param[index_v] / (double) *number_parameters); initial_user_parameter_temp[index_v] = new_temperature * pow(initial_user_parameter_temp [index_v] / new_temperature, temperature_rescale_power); } /* restore from temporary double */ index_parameter_generations[index_v] = (LONG_INT) temp_dbl; } /* reanneal : Reset the index of cost acceptances to take into account finer detail in cost terrain. */ /* set the starting cost_temperature to last cost found so far */ if (*initial_cost_temperature > fabs(last_saved_state->cost)) *initial_cost_temperature = fabs(last_saved_state->cost); temp_dbl = (double) *index_cost_acceptances; if (*current_cost_temperature > fabs(best_generated_state->cost)) { tmp_var_db3 = fabs(log(((double) EPS_DOUBLE + *initial_cost_temperature) / ((double) EPS_DOUBLE + fabs(best_generated_state->cost)))); temp_dbl = (double) EPS_DOUBLE + pow(tmp_var_db3 / *temperature_scale_cost, (double) *number_parameters / (*quench_cost)); } else { log_init_cur_temp_ratio = fabs(log(((double) EPS_DOUBLE + *initial_cost_temperature) / ((double) EPS_DOUBLE + *current_cost_temperature))); temp_dbl = (double) EPS_DOUBLE + pow(log_init_cur_temp_ratio / *temperature_scale_cost, (double) *number_parameters / (*quench_cost)); } /* reset index_cost_temperature if index reset too large */ while (temp_dbl > ((double) OPTIONS->MAXIMUM_REANNEAL_INDEX)) { log_new_temperature_ratio = -*temperature_scale_cost * pow((double) temp_dbl, (*quench_cost) / (double) *number_parameters); log_new_temperature_ratio = EXPONENT_CHECK(log_new_temperature_ratio); new_temperature = *initial_cost_temperature * exp(log_new_temperature_ratio); temp_dbl /= OPTIONS->REANNEAL_RESCALE; temperature_rescale_power = ONE / pow(OPTIONS->REANNEAL_RESCALE, (*quench_cost) / (double) *number_parameters); *initial_cost_temperature = new_temperature * pow(*initial_cost_temperature / new_temperature, temperature_rescale_power); } *index_cost_acceptances = (LONG_INT) temp_dbl; } /*********************************************************************** * cost_derivatives * This procedure calculates the derivatives of the cost function * with respect to its parameters. The first derivatives are * used as a sensitivity measure for reannealing. The second * derivatives are calculated only if *curvature_flag=TRUE; * these are a measure of the covariance of the fit when a * minimum is found. ***********************************************************************/ /* Calculate the numerical derivatives of the best generated state found so far */ /* In this implementation of ASA, no checks are made for *valid_state_generated_flag=FALSE for differential neighbors to the current best state. */ /* Assuming no information is given about the metric of the parameter space, use simple Cartesian space to calculate curvatures. */ #if HAVE_ANSI void cost_derivatives(double (*user_cost_function) (), double *parameter_minimum, double *parameter_maximum, double *tangents, double *curvature, double *maximum_tangent, double *delta_parameter, ALLOC_INT * number_parameters, int *parameter_type, int *exit_status, int *curvature_flag, int *valid_state_generated_flag, LONG_INT * number_invalid_generated_states, STATE * current_generated_state, STATE * best_generated_state, FILE * ptr_asa_out, DEFINES * OPTIONS) #else void cost_derivatives(user_cost_function, parameter_minimum, parameter_maximum, tangents, curvature, maximum_tangent, delta_parameter, number_parameters, parameter_type, exit_status, curvature_flag, valid_state_generated_flag, number_invalid_generated_states, current_generated_state, best_generated_state, ptr_asa_out, OPTIONS) double (*user_cost_function) (); double *parameter_minimum; double *parameter_maximum; double *tangents; double *curvature; double *maximum_tangent; double *delta_parameter; ALLOC_INT *number_parameters; int *parameter_type; int *exit_status; int *curvature_flag; int *valid_state_generated_flag; LONG_INT *number_invalid_generated_states; STATE *current_generated_state; STATE *best_generated_state; FILE *ptr_asa_out; DEFINES *OPTIONS; #endif { ALLOC_INT index_v, index_vv, index_v_vv, index_vv_v; LONG_INT saved_num_invalid_gen_states; double parameter_v, parameter_vv, parameter_v_offset; double recent_best_cost; double new_cost_state_1, new_cost_state_2, new_cost_state_3; double delta_parameter_v, delta_parameter_vv; if (OPTIONS->CURVATURE_0 == TRUE) *curvature_flag = FALSE; /* save the best cost */ recent_best_cost = best_generated_state->cost; /* copy the best state into the current state */ VFOR(index_v) { /* ignore parameters with too small ranges */ if (PARAMETER_RANGE_TOO_SMALL(index_v)) continue; current_generated_state->parameter[index_v] = best_generated_state->parameter[index_v]; } /* set parameters (& possibly constraints) to best state */ saved_num_invalid_gen_states = (*number_invalid_generated_states); *valid_state_generated_flag = TRUE; 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 ASA_PRINT if (*valid_state_generated_flag == FALSE) /* extra check */ fprintf(ptr_asa_out, "Generated an invalid best state when calculating \ the derivatives\n"); #endif (*number_invalid_generated_states) = saved_num_invalid_gen_states; *valid_state_generated_flag = TRUE; VFOR(index_v) { if (NO_REANNEAL(index_v)) { tangents[index_v] = ZERO; index_v_vv = ROW_COL_INDEX(index_v, index_v); curvature[index_v_vv] = ZERO; continue; } /* skip parameters with too small range or integer parameters */ if (OPTIONS->INCLUDE_INTEGER_PARAMETERS == TRUE) { if (PARAMETER_RANGE_TOO_SMALL(index_v)) { tangents[index_v] = ZERO; if (*curvature_flag == TRUE) { index_v_vv = ROW_COL_INDEX(index_v, index_v); curvature[index_v_vv] = ZERO; continue; } } } else { if (PARAMETER_RANGE_TOO_SMALL(index_v) || INTEGER_PARAMETER(index_v)) { tangents[index_v] = ZERO; if (*curvature_flag == TRUE) { index_v_vv = ROW_COL_INDEX(index_v, index_v); curvature[index_v_vv] = ZERO; continue; } } } /* save the v_th parameter and delta_parameter */ parameter_v = best_generated_state->parameter[index_v]; delta_parameter_v = delta_parameter[index_v]; /* This does _not_ fix going out of bounds by delta_parameter_v or delta_parameter_vv during a curvature calculation, but does address such problems during reannealing. */ parameter_v_offset = (ONE + delta_parameter_v) * parameter_v; if (parameter_v_offset > parameter_maximum[index_v] || parameter_v_offset < parameter_minimum[index_v]) delta_parameter_v = -delta_parameter_v; /* generate the first sample point and calculate cost at current point + OPTIONS->DELTA_X */ current_generated_state->parameter[index_v] = parameter_v_offset; /* generate the first sample point */ 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); new_cost_state_1 = current_generated_state->cost; *valid_state_generated_flag = TRUE; /* restore the parameter state */ current_generated_state->parameter[index_v] = parameter_v; /* calculate the numerical derivative */ tangents[index_v] = (new_cost_state_1 - recent_best_cost) / (delta_parameter_v * parameter_v + (double) EPS_DOUBLE); /* calculate the diagonal curvatures */ if (*curvature_flag == TRUE) { /* generate the second sample point and calculate cost at current point - OPTIONS->DELTA_X */ current_generated_state->parameter[index_v] = (ONE - delta_parameter_v) * parameter_v; 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); new_cost_state_2 = current_generated_state->cost; *valid_state_generated_flag = TRUE; /* restore the parameter state */ current_generated_state->parameter[index_v] = parameter_v; /* index_v_vv: row index_v, column index_v */ index_v_vv = ROW_COL_INDEX(index_v, index_v); /* calculate and store the curvature */ curvature[index_v_vv] = (new_cost_state_2 - TWO * recent_best_cost + new_cost_state_1) / (delta_parameter_v * delta_parameter_v * parameter_v * parameter_v + (double) EPS_DOUBLE); } } /* calculate off-diagonal curvatures */ if (*curvature_flag == TRUE) { VFOR(index_v) { if (NO_REANNEAL(index_v)) continue; /* save the v_th parameter and delta_x */ parameter_v = current_generated_state->parameter[index_v]; delta_parameter_v = delta_parameter[index_v]; VFOR(index_vv) { if (NO_REANNEAL(index_vv)) continue; /* calculate only the upper diagonal */ if (index_v <= index_vv) continue; /* index_v_vv: row index_v, column index_vv */ index_v_vv = ROW_COL_INDEX(index_v, index_vv); index_vv_v = ROW_COL_INDEX(index_vv, index_v); /* skip parms with too small range or integer parameters */ if (OPTIONS->INCLUDE_INTEGER_PARAMETERS == TRUE) { if (PARAMETER_RANGE_TOO_SMALL(index_v) || PARAMETER_RANGE_TOO_SMALL(index_vv)) { curvature[index_vv_v] = curvature[index_v_vv] = ZERO; continue; } } else { if (INTEGER_PARAMETER(index_v) || INTEGER_PARAMETER(index_vv) || PARAMETER_RANGE_TOO_SMALL(index_v) || PARAMETER_RANGE_TOO_SMALL(index_vv)) { curvature[index_vv_v] = curvature[index_v_vv] = ZERO; continue; } } /* save the vv_th parameter and delta_parameter */ parameter_vv = current_generated_state->parameter[index_vv]; delta_parameter_vv = delta_parameter[index_vv]; /* generate first sample point */ current_generated_state->parameter[index_v] = (ONE + delta_parameter_v) * parameter_v; current_generated_state->parameter[index_vv] = (ONE + delta_parameter_vv) * parameter_vv; 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); new_cost_state_1 = current_generated_state->cost; *valid_state_generated_flag = TRUE; /* restore the v_th parameter */ current_generated_state->parameter[index_v] = parameter_v; /* generate second sample point */ 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); new_cost_state_2 = current_generated_state->cost; *valid_state_generated_flag = TRUE; /* restore the vv_th parameter */ current_generated_state->parameter[index_vv] = parameter_vv; /* generate third sample point */ current_generated_state->parameter[index_v] = (ONE + delta_parameter_v) * best_generated_state->parameter[index_v]; 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); new_cost_state_3 = current_generated_state->cost; *valid_state_generated_flag = TRUE; /* restore the v_th parameter */ current_generated_state->parameter[index_v] = parameter_v; /* calculate and store the curvature */ curvature[index_vv_v] = curvature[index_v_vv] = (new_cost_state_1 - new_cost_state_2 - new_cost_state_3 + recent_best_cost) / (delta_parameter_v * delta_parameter_vv * parameter_v * parameter_vv + (double) EPS_DOUBLE); } } } /* restore the best cost function value */ current_generated_state->cost = recent_best_cost; /* find the maximum |tangent| from all tangents */ *maximum_tangent = 0; VFOR(index_v) { if (NO_REANNEAL(index_v)) continue; /* ignore too small ranges and integer parameters types */ if (OPTIONS->INCLUDE_INTEGER_PARAMETERS == TRUE) { if (PARAMETER_RANGE_TOO_SMALL(index_v)) continue; } else { if (PARAMETER_RANGE_TOO_SMALL(index_v) || INTEGER_PARAMETER(index_v)) continue; } /* find the maximum |tangent| (from all tangents) */ if (fabs(tangents[index_v]) > *maximum_tangent) { *maximum_tangent = fabs(tangents[index_v]); } } } #if ASA_PRINT /*********************************************************************** * print_state * Prints a description of the current state of the system ***********************************************************************/ /* print the state of the system to the user */ #if HAVE_ANSI void print_state(double *parameter_minimum, double *parameter_maximum, double *tangents, double *curvature, double *current_cost_temperature, double *current_user_parameter_temp, double *accepted_to_generated_ratio, ALLOC_INT * number_parameters, int *curvature_flag, LONG_INT * number_accepted, LONG_INT * index_cost_acceptances, LONG_INT * number_generated, LONG_INT * number_invalid_generated_states, STATE * last_saved_state, STATE * best_generated_state, FILE * ptr_asa_out, DEFINES * OPTIONS) #else void 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) double *parameter_minimum; double *parameter_maximum; double *tangents; double *curvature; double *current_cost_temperature; double *current_user_parameter_temp; double *accepted_to_generated_ratio; ALLOC_INT *number_parameters; int *curvature_flag; LONG_INT *number_accepted; LONG_INT *index_cost_acceptances; LONG_INT *number_generated; LONG_INT *number_invalid_generated_states; STATE *last_saved_state; STATE *best_generated_state; FILE *ptr_asa_out; DEFINES *OPTIONS; #endif { ALLOC_INT index_v; ALLOC_INT index_vv, index_v_vv; fprintf(ptr_asa_out, "\n"); #if TIME_CALC print_asa_time("", ptr_asa_out); #endif if (OPTIONS->CURVATURE_0 == TRUE) *curvature_flag = FALSE; #if INT_LONG fprintf(ptr_asa_out, "*index_cost_acceptances = %ld, *current_cost_temperature = %12.7g\n", *index_cost_acceptances, *current_cost_temperature); fprintf(ptr_asa_out, "*accepted_to_generated_ratio = %12.7g,\ *number_invalid... = %ld\n", *accepted_to_generated_ratio, (*number_invalid_generated_states)); fprintf(ptr_asa_out, "*number_generated = %ld, *number_accepted = %ld\n", *number_generated, *number_accepted); #else fprintf(ptr_asa_out, "*index_cost_acceptances = %d, *current_cost_temperature = %12.7g\n", *index_cost_acceptances, *current_cost_temperature); fprintf(ptr_asa_out, "*accepted_to_generated_ratio = %12.7g,\ *number_invalid... = %d\n", *accepted_to_generated_ratio, *number_invalid_generated_states); fprintf(ptr_asa_out, "*number_generated = %d, *number_accepted = %d\n", *number_generated, *number_accepted); #endif fprintf(ptr_asa_out, "best...->cost = %12.7g,\ last...->cost = %12.7g\n", best_generated_state->cost, last_saved_state->cost); /* Note that tangents will not be calculated until reanneal is called, and therefore their listing in the printout only is relevant then */ VFOR(index_v) { /* ignore too small ranges */ if (PARAMETER_RANGE_TOO_SMALL(index_v)) continue; fprintf(ptr_asa_out, #if INT_ALLOC "best_generated_state->parameter[%d] = %12.7g\n\ current_user_parameter_temp[%d] = %12.7g\n\ tangents[%d]: %12.7g\n", #else #if INT_LONG "best_generated_state->parameter[%ld] = %12.7g\n\ current_user_parameter_temp[%ld] = %12.7g\n\ tangents[%ld]: %12.7g\n", #else "best_generated_state->parameter[%d] = %12.7g\n\ current_user_parameter_temp[%d] = %12.7g\n\ tangents[%d]: %12.7g\n", #endif #endif index_v, best_generated_state->parameter[index_v], index_v, current_user_parameter_temp[index_v], index_v, tangents[index_v]); } if (*curvature_flag == TRUE) { /* print curvatures */ VFOR(index_v) { /* ignore too small ranges */ if (PARAMETER_RANGE_TOO_SMALL(index_v)) continue; VFOR(index_vv) { /* only print upper diagonal of matrix */ if (index_v < index_vv) continue; /* ignore too small ranges (index_vv) */ if (PARAMETER_RANGE_TOO_SMALL(index_vv)) continue; /* index_v_vv: row index_v, column index_vv */ index_v_vv = ROW_COL_INDEX(index_v, index_vv); #if INT_ALLOC fprintf(ptr_asa_out, "curvature[%d][%d] = %12.7g\n", #else #if INT_LONG fprintf(ptr_asa_out, "curvature[%ld][%ld] = %12.7g\n", #else fprintf(ptr_asa_out, "curvature[%d][%d] = %12.7g\n", #endif #endif index_v, index_vv, curvature[index_v_vv]); } } } fprintf(ptr_asa_out, "\n"); fflush(ptr_asa_out); } #endif #if ASA_PRINT #if TIME_CALC /*********************************************************************** * print_asa_time * This calculates the time and runtime and prints it. ***********************************************************************/ #if HAVE_ANSI void print_asa_time(char *message, FILE * ptr_asa_out) #else void print_asa_time(message, ptr_asa_out) char *message; FILE *ptr_asa_out; #endif { int who = RUSAGE_SELF; /* Check our own time */ struct rusage usage; /* get the resource usage information */ #if TIME_STD syscall(SYS_GETRUSAGE, who, &usage); #else getrusage(who, &usage); #endif /* print the usage time in reasonable form */ aux_print_asa_time(&usage.ru_utime, message, ptr_asa_out); } /*********************************************************************** * aux_print_asa_time * auxiliary print the time routine ***********************************************************************/ #if HAVE_ANSI void aux_print_asa_time(struct timeval *time, char *message, FILE * ptr_asa_out) #else void aux_print_asa_time(time, message, ptr_asa_out) struct timeval *time; char *message; FILE *ptr_asa_out; #endif { static double sx; double us, s, m, h; double ds, dm, dh; /* calculate the new microseconds, seconds, minutes, hours and the differences since the last call */ us = (double) ((int) ((double) EPS_DOUBLE + time->tv_usec)) / 1.E6; s = (double) ((int) ((double) EPS_DOUBLE + time->tv_sec)) + us; ds = s - sx; sx = s; h = (int) ((double) EPS_DOUBLE + s / 3600.); m = (int) ((double) EPS_DOUBLE + s / 60.) - 60. * h; s -= (3600. * h + 60. * m); dh = (int) ((double) EPS_DOUBLE + ds / 3600.); dm = (int) ((double) EPS_DOUBLE + ds / 60.) - 60. * dh; ds -= (3600. * dh + 60. * dm); /* print the statistics */ fprintf(ptr_asa_out, "%s:time: %gh %gm %gs; incr: %gh %gm %gs\n", message, h, m, s, dh, dm, ds); } #endif /* TIME_CALC */ #endif /* ASA_PRINT */