/*********************************************************************** * Adaptive Simulated Annealing (ASA) * Lester Ingber * Copyright (c) 1993-2004 Lester Ingber. All Rights Reserved. * The LICENSE file must be included with ASA code. * Modified by John Meacham for Haskell interface ***********************************************************************/ #include "asa_usr.h" char user_exit_msg[160]; /* temp storage for exit messages */ FILE *ptr_out; static double resettable_randflt (LONG_INT * rand_seed, int reset); static double randflt (LONG_INT * rand_seed); /*********************************************************************** * main * This is a sample calling program to optimize using ASA ***********************************************************************/ int asa_main ( hs_cost_func *func, int number_parameters, double *upper_bounds, double *lower_bounds, int *type, double *main_cost_value, double *main_cost_parameters, int *main_exit_code, long int initial_rand_seed ) { int i; int *exit_code; ALLOC_INT n_param; #if ASA_TEMPLATE_SAMPLE FILE *ptr_asa; #endif #if MULTI_MIN int multi_index; #endif /* pointer to array storage for asa arguments */ double *parameter_lower_bound, *parameter_upper_bound, *cost_parameters, *cost_tangents, *cost_curvature; double cost_value; int initialize_parameters_value; /* the number of parameters to optimize */ ALLOC_INT *parameter_dimension; /* pointer to array storage for parameter type flags */ int *parameter_int_real; /* valid flag for cost function */ int *cost_flag; /* seed for random number generator */ LONG_INT *rand_seed; USER_DEFINES *USER_OPTIONS; #if MY_TEMPLATE /* MY_TEMPLATE_main_decl */ /* add some declarations if required */ #endif #if ASA_TEMPLATE_MULTIPLE int n_asa, n_trajectory; ALLOC_INT index; #if HAVE_ANSI char asa_file[8] = "asa_x_y"; #else char asa_file[8]; #endif /* HAVE_ANSI */ #endif /* ASA_TEMPLATE_MULTIPLE */ if ((USER_OPTIONS = (USER_DEFINES *) calloc (1, sizeof (USER_DEFINES))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): USER_DEFINES"); Exit_USER (user_exit_msg); return (-2); } #if OPTIONAL_DATA_PTR #if ASA_TEMPLATE USER_OPTIONS->Asa_Data_Dim_Ptr = 256; if ((USER_OPTIONS->Asa_Data_Ptr = (OPTIONAL_PTR_TYPE *) calloc (USER_OPTIONS->Asa_Data_Dim_Ptr, sizeof (OPTIONAL_PTR_TYPE))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): USER_OPTIONS->Asa_Data_Ptr"); Exit_USER (user_exit_msg); return (-2); } #endif /* ASA_TEMPLATE */ #endif /* OPTIONAL_DATA_PTR */ if (!strcmp (USER_OUT, "STDOUT")) { #if INCL_STDOUT ptr_out = stdout; #endif /* INCL_STDOUT */ } else { ptr_out = fopen (USER_OUT, "w"); } fflush (ptr_out); if ((rand_seed = (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): rand_seed"); Exit_USER (user_exit_msg); return (-2); } *rand_seed = initial_rand_seed; /* initialize random number generator with first call */ resettable_randflt (rand_seed, 1); /* Initialize the users parameters, allocating space, etc. Note that the default is to have asa generate the initial cost_parameters that satisfy the user's constraints. */ if ((parameter_dimension = (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): parameter_dimension"); Exit_USER (user_exit_msg); return (-2); } if ((exit_code = (int *) calloc (1, sizeof (int))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): exit_code"); Exit_USER (user_exit_msg); return (-2); } if ((cost_flag = (int *) calloc (1, sizeof (int))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): cost_flag"); Exit_USER (user_exit_msg); return (-2); } //USER_OPTIONS->Limit_Acceptances = 10000; USER_OPTIONS->Limit_Acceptances = 1000; USER_OPTIONS->Limit_Generated = 99999; USER_OPTIONS->Limit_Invalid_Generated_States = 1000; /* USER_OPTIONS->Accepted_To_Generated_Ratio = 1.0E-6; */ USER_OPTIONS->Accepted_To_Generated_Ratio = 1.0E-4; USER_OPTIONS->Cost_Precision = 1.0E-18; USER_OPTIONS->Maximum_Cost_Repeat = 5; USER_OPTIONS->Number_Cost_Samples = 5; USER_OPTIONS->Temperature_Ratio_Scale = 1.0E-5; USER_OPTIONS->Cost_Parameter_Scale_Ratio = 1.0; USER_OPTIONS->Temperature_Anneal_Scale = 100.0; USER_OPTIONS->Include_Integer_Parameters = FALSE; USER_OPTIONS->User_Initial_Parameters = FALSE; USER_OPTIONS->Sequential_Parameters = -1; USER_OPTIONS->Initial_Parameter_Temperature = 1.0; USER_OPTIONS->Acceptance_Frequency_Modulus = 100; USER_OPTIONS->Generated_Frequency_Modulus = 10000; USER_OPTIONS->Reanneal_Cost = 1; USER_OPTIONS->Reanneal_Parameters = TRUE; USER_OPTIONS->Delta_X = 0.001; USER_OPTIONS->User_Tangents = FALSE; USER_OPTIONS->Curvature_0 = FALSE; /* ALLOCATE STORAGE */ #if USER_ASA_OUT if ((USER_OPTIONS->Asa_Out_File = (char *) calloc (80, sizeof (char))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): USER_OPTIONS->Asa_Out_File"); Exit_USER (user_exit_msg); return (-2); } #endif /* the number of parameters for the cost function */ #if OPTIONS_FILE_DATA fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%s", read_option); #if INT_ALLOC fscanf (ptr_options, "%d", &read_int); *parameter_dimension = read_int; #else #if INT_LONG fscanf (ptr_options, "%ld", &read_long); *parameter_dimension = read_long; #else fscanf (ptr_options, "%d", &read_int); *parameter_dimension = read_int; #endif #endif #else /* OPTIONS_FILE_DATA */ #endif /* OPTIONS_FILE_DATA */ #if MY_TEMPLATE /* MY_TEMPLATE_dim */ *parameter_dimension = number_parameters; /* If not using OPTIONS_FILE_DATA or data read from asa_opt, insert the number of parameters for the cost_function */ #endif /* MY_TEMPLATE dim */ #if ASA_TEMPLATE_SAMPLE *parameter_dimension = 2; USER_OPTIONS->Limit_Acceptances = 2000; USER_OPTIONS->User_Tangents = TRUE; USER_OPTIONS->Limit_Weights = 1.0E-7; #endif #if ASA_TEMPLATE_PARALLEL USER_OPTIONS->Gener_Block = 100; USER_OPTIONS->Gener_Block_Max = 512; USER_OPTIONS->Gener_Mov_Avr = 3; #endif /* allocate parameter minimum space */ if ((parameter_lower_bound = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): parameter_lower_bound"); Exit_USER (user_exit_msg); return (-2); } /* allocate parameter maximum space */ if ((parameter_upper_bound = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): parameter_upper_bound"); Exit_USER (user_exit_msg); return (-2); } /* allocate parameter initial values; the parameter final values will be stored here later */ if ((cost_parameters = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): cost_parameters"); Exit_USER (user_exit_msg); return (-2); } /* allocate the parameter types, real or integer */ if ((parameter_int_real = (int *) calloc (*parameter_dimension, sizeof (int))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): parameter_int_real"); Exit_USER (user_exit_msg); return (-2); } /* allocate space for parameter cost_tangents - used for reannealing */ if ((cost_tangents = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): cost_tangents"); Exit_USER (user_exit_msg); return (-2); } if (USER_OPTIONS->Curvature_0 == FALSE || USER_OPTIONS->Curvature_0 == -1) { /* allocate space for parameter cost_curvatures/covariance */ if ((cost_curvature = (double *) calloc ((*parameter_dimension) * (*parameter_dimension), sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): cost_curvature"); Exit_USER (user_exit_msg); return (-2); } } else { cost_curvature = (double *) NULL; } #if USER_COST_SCHEDULE USER_OPTIONS->Cost_Schedule = user_cost_schedule; #endif #if USER_ACCEPTANCE_TEST USER_OPTIONS->Acceptance_Test = user_acceptance_test; #endif #if USER_ACCEPT_ASYMP_EXP USER_OPTIONS->Asymp_Exp_Param = 1.0; #endif #if USER_GENERATING_FUNCTION USER_OPTIONS->Generating_Distrib = user_generating_distrib; #endif #if USER_REANNEAL_COST USER_OPTIONS->Reanneal_Cost_Function = user_reanneal_cost; #endif #if USER_REANNEAL_PARAMETERS USER_OPTIONS->Reanneal_Params_Function = user_reanneal_params; #endif #if MY_TEMPLATE /* MY_TEMPLATE_pre_initialize */ /* last changes before entering initialize_parameters() */ USER_OPTIONS->Asa_Data_Ptr = func; USER_OPTIONS->Asa_Data_Dim_Ptr = 1; memcpy(parameter_lower_bound,lower_bounds,sizeof(double)*number_parameters); memcpy(cost_parameters,lower_bounds,sizeof(double)*number_parameters); memcpy(parameter_upper_bound,upper_bounds,sizeof(double)*number_parameters); memcpy(parameter_int_real, type, sizeof(int)*number_parameters); #endif initialize_parameters_value = initialize_parameters (cost_parameters, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, #if OPTIONS_FILE_DATA ptr_options, #endif USER_OPTIONS); if (initialize_parameters_value == -2) return (initialize_parameters_value); for(i = 0; i < number_parameters; i++) { USER_OPTIONS->User_Quench_Param_Scale[i] = 1.0; } USER_OPTIONS->User_Quench_Cost_Scale[0] = 1.0; /* optimize the cost_function, returning the results in cost_value and cost_parameters */ #if ASA_TEMPLATE_MULTIPLE /* multiple asa() quenched calls + multiple asa_out files (To get longer quenched runs, decrease SMALL_FLOAT.) */ for (n_asa = 1; n_asa <= *parameter_dimension; n_asa++) { asa_file[4] = 'A' + n_asa - 1; USER_OPTIONS->User_Quench_Cost_Scale[0] = (double) n_asa; for (index = 0; index < *parameter_dimension; ++index) USER_OPTIONS->User_Quench_Param_Scale[index] = (double) n_asa; for (n_trajectory = 0; n_trajectory < 3; ++n_trajectory) { asa_file[6] = 'a' + n_trajectory; strcpy (USER_OPTIONS->Asa_Out_File, asa_file); #endif #if ASA_TEMPLATE_ASA_OUT_PID pid_file[0] = 'a'; pid_file[1] = 's'; pid_file[2] = 'a'; pid_file[3] = '_'; pid_file[4] = 'o'; pid_file[5] = 'u'; pid_file[6] = 't'; pid_file[7] = '_'; pid_int = getpid (); if (pid_int < 0) { pid_file[7] = '0'; pid_int = -pid_int; } strcpy (USER_OPTIONS->Asa_Out_File, pid_file); #endif cost_value = asa (USER_COST_FUNCTION, randflt, rand_seed, cost_parameters, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, USER_OPTIONS); if (*exit_code == -1) { #if INCL_STDOUT printf ("\n\n*** error in calloc in ASA ***\n\n"); #endif /* INCL_STDOUT */ fprintf (ptr_out, "\n\n*** error in calloc in ASA ***\n\n"); fflush (ptr_out); return (-1); } #if MULTI_MIN fprintf (ptr_out, "Multi_Specify = %d\n", USER_OPTIONS->Multi_Specify); #if INT_LONG fprintf (ptr_out, "N_Accepted = %ld\n", USER_OPTIONS->N_Accepted); #else fprintf (ptr_out, "N_Accepted = %d\n", USER_OPTIONS->N_Accepted); #endif #if ASA_RESOLUTION for (n_param = 0; n_param < *parameter_dimension; ++n_param) { fprintf (ptr_out, #if INT_ALLOC "Coarse_Resolution[%d] = %12.7g\n", #else #if INT_LONG "Coarse_Resolution[%ld] = %12.7g\n", #else "Coarse_Resolution[%d] = %12.7g\n", #endif #endif n_param, USER_OPTIONS->Coarse_Resolution[n_param]); } #else /* ASA_RESOLUTION */ for (n_param = 0; n_param < *parameter_dimension; ++n_param) { fprintf (ptr_out, #if INT_ALLOC "Multi_Grid[%d] = %12.7g\n", #else #if INT_LONG "Multi_Grid[%ld] = %12.7g\n", #else "Multi_Grid[%d] = %12.7g\n", #endif #endif n_param, USER_OPTIONS->Multi_Grid[n_param]); } #endif /* ASA_RESOLUTION */ fprintf (ptr_out, "\n"); for (multi_index = 0; multi_index < USER_OPTIONS->Multi_Number; ++multi_index) { fprintf (ptr_out, "\n"); fprintf (ptr_out, "Multi_Cost[%d] = %12.7g\n", multi_index, USER_OPTIONS->Multi_Cost[multi_index]); for (n_param = 0; n_param < *parameter_dimension; ++n_param) { fprintf (ptr_out, #if INT_ALLOC "Multi_Params[%d][%d] = %12.7g\n", #else #if INT_LONG "Multi_Params[%d][%ld] = %12.7g\n", #else "Multi_Params[%d][%d] = %12.7g\n", #endif #endif multi_index, n_param, USER_OPTIONS->Multi_Params[multi_index][n_param]); } } fprintf (ptr_out, "\n"); fflush (ptr_out); #endif /* MULTI_MIN */ #if FITLOC /* Fit_Local, Iter_Max and Penalty may be set adaptively */ USER_OPTIONS->Penalty = 1000; USER_OPTIONS->Fit_Local = 0; USER_OPTIONS->Iter_Max = 500; if (USER_OPTIONS->Fit_Local >= 1) { cost_value = fitloc (USER_COST_FUNCTION, cost_parameters, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, USER_OPTIONS, ptr_out); } #endif /* FITLOC */ #if MY_TEMPLATE /* MY_TEMPLATE_post_asa */ #endif *main_cost_value = cost_value; for (n_param = 0; n_param < *parameter_dimension; ++n_param) { main_cost_parameters[n_param] = cost_parameters[n_param]; } *main_exit_code = *exit_code; fprintf (ptr_out, "exit code = %d\n", *exit_code); fprintf (ptr_out, "final cost value = %-12.7g\n", cost_value); fprintf (ptr_out, "%12s %12s\n","parameter","value"); for (n_param = 0; n_param < *parameter_dimension; ++n_param) { fprintf (ptr_out, #if INT_ALLOC "%12d %12.7g\n", #else #if INT_LONG "%12ld %12.7g\n", #else "%12d %12.7g\n", #endif #endif n_param, cost_parameters[n_param]); } #if ASA_TEMPLATE_MULTIPLE } } #endif #if ASA_TEMPLATE_SAMPLE ptr_asa = fopen ("asa_out", "r"); sample (ptr_out, ptr_asa); #endif /* close all files */ ptr_out != stdout && fclose (ptr_out); #if OPTIONAL_DATA_DBL free (USER_OPTIONS->Asa_Data_Dbl); #endif #if OPTIONAL_DATA_INT free (USER_OPTIONS->Asa_Data_Int); #endif #if OPTIONAL_DATA_PTR #if MY_TEMPLATE /* Instead of freeing Asa_Data_Ptr, if memory has been allocated * outside ASA, e.g., by the use of ASA_LIB, use the following: */ USER_OPTIONS->Asa_Data_Ptr = NULL; #endif /* MY_TEMPLATE */ free (USER_OPTIONS->Asa_Data_Ptr); #endif #if USER_ASA_OUT free (USER_OPTIONS->Asa_Out_File); #endif #if ASA_SAMPLE free (USER_OPTIONS->Bias_Generated); #endif #if ASA_QUEUE #if ASA_RESOLUTION #else free (USER_OPTIONS->Queue_Resolution); #endif #endif #if ASA_RESOLUTION free (USER_OPTIONS->Coarse_Resolution); #endif if (USER_OPTIONS->Curvature_0 == FALSE || USER_OPTIONS->Curvature_0 == -1) free (cost_curvature); #if USER_INITIAL_PARAMETERS_TEMPS free (USER_OPTIONS->User_Parameter_Temperature); #endif #if USER_INITIAL_COST_TEMP free (USER_OPTIONS->User_Cost_Temperature); #endif #if DELTA_PARAMETERS free (USER_OPTIONS->User_Delta_Parameter); #endif #if QUENCH_PARAMETERS free (USER_OPTIONS->User_Quench_Param_Scale); #endif #if QUENCH_COST free (USER_OPTIONS->User_Quench_Cost_Scale); #endif #if RATIO_TEMPERATURE_SCALES free (USER_OPTIONS->User_Temperature_Ratio); #endif #if MULTI_MIN free (USER_OPTIONS->Multi_Cost); free (USER_OPTIONS->Multi_Grid); for (multi_index = 0; multi_index < USER_OPTIONS->Multi_Number; ++multi_index) { free (USER_OPTIONS->Multi_Params[multi_index]); } free (USER_OPTIONS->Multi_Params); #endif /* MULTI_MIN */ free (USER_OPTIONS); free (parameter_dimension); free (exit_code); free (cost_flag); free (parameter_lower_bound); free (parameter_upper_bound); free (cost_parameters); free (parameter_int_real); free (cost_tangents); free (rand_seed); return (0); /* NOTREACHED */ } /*********************************************************************** * initialize_parameters - sample parameter initialization function * This depends on the users cost function to optimize (minimum). * The routine allocates storage needed for asa. The user should * define the number of parameters and their ranges, * and make sure the initial parameters are within * the minimum and maximum ranges. The array * parameter_int_real should be REAL_TYPE (-1) for real parameters, * and INTEGER_TYPE (1) for integer values ***********************************************************************/ #if HAVE_ANSI int initialize_parameters (double *cost_parameters, double *parameter_lower_bound, double *parameter_upper_bound, double *cost_tangents, double *cost_curvature, ALLOC_INT * parameter_dimension, int *parameter_int_real, #if OPTIONS_FILE_DATA FILE * ptr_options, #endif USER_DEFINES * USER_OPTIONS) #else int initialize_parameters (cost_parameters, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, #if OPTIONS_FILE_DATA ptr_options, #endif USER_OPTIONS) double *cost_parameters; double *parameter_lower_bound; double *parameter_upper_bound; double *cost_tangents; double *cost_curvature; ALLOC_INT *parameter_dimension; int *parameter_int_real; #if OPTIONS_FILE_DATA FILE *ptr_options; #endif USER_DEFINES *USER_OPTIONS; #endif { ALLOC_INT index; #if OPTIONS_FILE_DATA char read_option[80]; ALLOC_INT read_index; #endif #if MULTI_MIN int multi_index; #endif #if MY_TEMPLATE /* MY_TEMPLATE_init_decl */ /* add some declarations if required */ #endif index = 0; #if OPTIONS_FILE_DATA fscanf (ptr_options, "%s", read_option); for (index = 0; index < *parameter_dimension; ++index) { #if MY_TEMPLATE /* MY_TEMPLATE_read_opt */ /* put in some code as required to alter lines read from asa_opt */ #endif #if INT_ALLOC fscanf (ptr_options, "%d", &read_index); #else #if INT_LONG fscanf (ptr_options, "%ld", &read_index); #else fscanf (ptr_options, "%d", &read_index); #endif #endif fscanf (ptr_options, "%lf%lf%lf%d", &(parameter_lower_bound[read_index]), &(parameter_upper_bound[read_index]), &(cost_parameters[read_index]), &(parameter_int_real[read_index])); } #else /* OPTIONS_FILE_DATA */ #if ASA_TEST /* store the parameter ranges */ for (index = 0; index < *parameter_dimension; ++index) parameter_lower_bound[index] = -10000.0; for (index = 0; index < *parameter_dimension; ++index) parameter_upper_bound[index] = 10000.0; /* store the initial parameter types */ for (index = 0; index < *parameter_dimension; ++index) parameter_int_real[index] = REAL_TYPE; /* store the initial parameter values */ for (index = 0; index < *parameter_dimension / 4.0; ++index) { cost_parameters[4 * (index + 1) - 4] = 999.0; cost_parameters[4 * (index + 1) - 3] = -1007.0; cost_parameters[4 * (index + 1) - 2] = 1001.0; cost_parameters[4 * (index + 1) - 1] = -903.0; } #endif /* ASA_TEST */ #endif /* OPTIONS_FILE_DATA */ #if ASA_TEMPLATE_SAMPLE for (index = 0; index < *parameter_dimension; ++index) parameter_lower_bound[index] = 0; for (index = 0; index < *parameter_dimension; ++index) parameter_upper_bound[index] = 2.0; for (index = 0; index < *parameter_dimension; ++index) parameter_int_real[index] = REAL_TYPE; for (index = 0; index < *parameter_dimension; ++index) cost_parameters[index] = 0.5; #endif #if USER_INITIAL_PARAMETERS_TEMPS if ((USER_OPTIONS->User_Parameter_Temperature = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->User_Parameter_Temperature"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE for (index = 0; index < *parameter_dimension; ++index) USER_OPTIONS->User_Parameter_Temperature[index] = 1.0; #endif #endif /* USER_INITIAL_PARAMETERS_TEMPS */ #if USER_INITIAL_COST_TEMP if ((USER_OPTIONS->User_Cost_Temperature = (double *) calloc (1, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->User_Cost_Temperature"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE USER_OPTIONS->User_Cost_Temperature[0] = 5.936648E+09; #endif #endif /* USER_INITIAL_COST_TEMP */ #if DELTA_PARAMETERS if ((USER_OPTIONS->User_Delta_Parameter = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->User_Delta_Parameter"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE for (index = 0; index < *parameter_dimension; ++index) USER_OPTIONS->User_Delta_Parameter[index] = 0.001; #endif #endif /* DELTA_PARAMETERS */ #if QUENCH_PARAMETERS if ((USER_OPTIONS->User_Quench_Param_Scale = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->User_Quench_Param_Scale"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE for (index = 0; index < *parameter_dimension; ++index) USER_OPTIONS->User_Quench_Param_Scale[index] = 1.0; #endif #if ASA_TEMPLATE_MULTIPLE for (index = 0; index < *parameter_dimension; ++index) USER_OPTIONS->User_Quench_Param_Scale[index] = 1.0; #endif #if ASA_TEMPLATE_SAVE for (index = 0; index < *parameter_dimension; ++index) USER_OPTIONS->User_Quench_Param_Scale[index] = 1.0; #endif #endif /* QUENCH_PARAMETERS */ #if QUENCH_COST if ((USER_OPTIONS->User_Quench_Cost_Scale = (double *) calloc (1, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->User_Quench_Cost_Scale"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE USER_OPTIONS->User_Quench_Cost_Scale[0] = 1.0; #endif #if ASA_TEMPLATE_MULTIPLE USER_OPTIONS->User_Quench_Cost_Scale[0] = 1.0; #endif #if ASA_TEMPLATE_SAVE USER_OPTIONS->User_Quench_Cost_Scale[0] = 1.0; #endif #endif /* QUENCH_COST */ /* use asa_opt to read in QUENCH USER_OPTIONS */ #if OPTIONS_FILE_DATA #if QUENCH_COST fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%lf", &(USER_OPTIONS->User_Quench_Cost_Scale[0])); #if QUENCH_PARAMETERS fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%s", read_option); for (index = 0; index < *parameter_dimension; ++index) { #if INT_ALLOC fscanf (ptr_options, "%d", &read_index); #else #if INT_LONG fscanf (ptr_options, "%ld", &read_index); #else fscanf (ptr_options, "%d", &read_index); #endif #endif fscanf (ptr_options, "%lf", &(USER_OPTIONS->User_Quench_Param_Scale[read_index])); } #endif /* QUENCH_PARAMETERS */ #endif /* QUENCH_COST */ #endif /* OPTIONS_FILE_DATA */ #if RATIO_TEMPERATURE_SCALES if ((USER_OPTIONS->User_Temperature_Ratio = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->User_Temperature_Ratio"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE for (index = 0; index < *parameter_dimension; ++index) USER_OPTIONS->User_Temperature_Ratio[index] = 1.0; #endif #endif /* RATIO_TEMPERATURE_SCALES */ /* Defines the limit of collection of sampled data by asa */ #if ASA_SAMPLE /* create memory for Bias_Generated[] */ if ((USER_OPTIONS->Bias_Generated = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->Bias_Generated"); Exit_USER (user_exit_msg); return (-2); } #endif #if ASA_RESOLUTION if ((USER_OPTIONS->Coarse_Resolution = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->Coarse_Resolution"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE for (index = 0; index < *parameter_dimension; ++index) USER_OPTIONS->Coarse_Resolution[index] = 1.0; #endif #endif /* ASA_RESOLUTION */ #if ASA_QUEUE #if ASA_RESOLUTION USER_OPTIONS->Queue_Resolution = USER_OPTIONS->Coarse_Resolution; #else /* ASA_RESOLUTION */ if ((USER_OPTIONS->Queue_Resolution = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->Queue_Resolution"); Exit_USER (user_exit_msg); return (-2); } #endif /* ASA_RESOLUTION */ #if ASA_TEMPLATE_QUEUE USER_OPTIONS->Queue_Size = 100; for (index = 0; index < *parameter_dimension; ++index) USER_OPTIONS->Queue_Resolution[index] = 0.001; #endif #endif /* ASA_QUEUE */ #if MULTI_MIN #if ASA_TEMPLATE USER_OPTIONS->Multi_Number = 2; #endif if ((USER_OPTIONS->Multi_Cost = (double *) calloc (USER_OPTIONS->Multi_Number, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->Multi_Cost"); Exit_USER (user_exit_msg); return (-2); } if ((USER_OPTIONS->Multi_Grid = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->Multi_Grid"); Exit_USER (user_exit_msg); return (-2); } if ((USER_OPTIONS->Multi_Params = (double **) calloc (USER_OPTIONS->Multi_Number, sizeof (double *))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->Multi_Params"); Exit_USER (user_exit_msg); return (-2); } for (multi_index = 0; multi_index < USER_OPTIONS->Multi_Number; ++multi_index) { if ((USER_OPTIONS->Multi_Params[multi_index] = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): USER_OPTIONS->Multi_Params[multi_index]"); Exit_USER (user_exit_msg); return (-2); } } #if ASA_TEST for (index = 0; index < *parameter_dimension; ++index) { USER_OPTIONS->Multi_Grid[index] = 0.05; } USER_OPTIONS->Multi_Specify = 0; #endif #if ASA_TEMPLATE for (index = 0; index < *parameter_dimension; ++index) { USER_OPTIONS->Multi_Grid[index] = (parameter_upper_bound[index] - parameter_lower_bound[index]) / 100.0; } USER_OPTIONS->Multi_Specify = 0; #endif /* ASA_TEMPLATE */ #endif /* MULTI_MIN */ USER_OPTIONS->Asa_Recursive_Level = 0; #if MY_TEMPLATE /* MY_TEMPLATE_params */ /* If not using RECUR_OPTIONS_FILE_DATA or data read from asa_opt, store the parameter ranges store the parameter types store the initial parameter values other changes needed for initialization */ #endif /* MY_TEMPLATE params */ return (0); } #if COST_FILE #else /*********************************************************************** * double cost_function * This is the users cost function to optimize * (find the minimum). * cost_flag is set to TRUE if the parameter set * does not violates any constraints * parameter_lower_bound and parameter_upper_bound may be * adaptively changed during the search. ***********************************************************************/ #if HAVE_ANSI double cost_function (double *x, double *parameter_lower_bound, double *parameter_upper_bound, double *cost_tangents, double *cost_curvature, ALLOC_INT * parameter_dimension, int *parameter_int_real, int *cost_flag, int *exit_code, USER_DEFINES * USER_OPTIONS) #else double cost_function (x, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, USER_OPTIONS) double *x; double *parameter_lower_bound; double *parameter_upper_bound; double *cost_tangents; double *cost_curvature; ALLOC_INT *parameter_dimension; int *parameter_int_real; int *cost_flag; int *exit_code; USER_DEFINES *USER_OPTIONS; #endif { #if ASA_TEST /* ASA test problem */ /* Objective function from * %A A. Corana * %A M. Marchesi * %A C. Martini * %A S. Ridella * %T Minimizing multimodal functions of continuous variables * with the "simulated annealing" algorithm * %J ACM Trans. Mathl. Software * %V 13 * %N 3 * %P 262-279 * %D 1987 * * This function, when used with ASA_TEST_POINT set to TRUE, contains * 1.0E20 local minima. When *parameter_dimension is equal to 4, visiting * each minimum for a millisecond would take about the present age of the * universe to visit all these minima. */ /* defines for the test problem, which assume *parameter_dimension is a multiple of 4. If this is set to a large number, you likely should set Curvature_0 to TRUE. */ double q_n, d_i, s_i, t_i, z_i, c_r; int k_i; #if ASA_TEST_POINT ALLOC_INT k_flag; #endif ALLOC_INT i, j; #if SELF_OPTIMIZE #else static LONG_INT funevals = 0; #endif #if ASA_TEMPLATE_SAVE static int read_test = 0; FILE *ptr_read_test; #endif #if MY_TEMPLATE /* MY_TEMPLATE_diminishing_ranges */ /* insert code to automate changing ranges of parameters */ #endif #if ASA_TEMPLATE /* example of diminishing ranges */ if (USER_OPTIONS->Locate_Cost == 12 && *(USER_OPTIONS->Best_Cost) < 1.0) { fprintf (ptr_out, "best_cost = %g\n", *(USER_OPTIONS->Best_Cost)); for (i = 0; i < *parameter_dimension; ++i) { parameter_lower_bound[i] = USER_OPTIONS->Best_Parameters[i] - 0.5 * fabs (parameter_lower_bound[i] - USER_OPTIONS->Best_Parameters[i]); parameter_upper_bound[i] = USER_OPTIONS->Best_Parameters[i] + 0.5 * fabs (parameter_upper_bound[i] - USER_OPTIONS->Best_Parameters[i]); parameter_lower_bound[i] = MIN (parameter_lower_bound[i], USER_OPTIONS->Best_Parameters[i] - 0.01); parameter_upper_bound[i] = MAX (parameter_upper_bound[i], USER_OPTIONS->Best_Parameters[i] + 0.01); } } #endif /* ASA_TEMPLATE */ /* a_i = parameter_upper_bound[i] */ s_i = 0.2; t_i = 0.05; c_r = 0.15; #if ASA_TEST_POINT k_flag = 0; for (i = 0; i < *parameter_dimension; ++i) { if (fabs (parameter_upper_bound[i] - parameter_lower_bound[i]) < (double) EPS_DOUBLE) continue; if (x[i] > 0.0) { k_i = (int) (x[i] / s_i + 0.5); } else if (x[i] < 0.0) { k_i = (int) (x[i] / s_i - 0.5); } else { k_i = 0; } if (k_i == 0) ++k_flag; } #endif /* ASA_TEST_POINT */ q_n = 0.0; for (i = 0; i < *parameter_dimension; ++i) { if (fabs (parameter_upper_bound[i] - parameter_lower_bound[i]) < (double) EPS_DOUBLE) continue; j = i % 4; switch (j) { case 0: d_i = 1.0; break; case 1: d_i = 1000.0; break; case 2: d_i = 10.0; break; default: d_i = 100.0; } if (x[i] > 0.0) { k_i = (int) (x[i] / s_i + 0.5); } else if (x[i] < 0.0) { k_i = (int) (x[i] / s_i - 0.5); } else { k_i = 0; } #if ASA_TEST_POINT if (fabs (k_i * s_i - x[i]) < t_i && k_flag != *parameter_dimension) #else if (fabs (k_i * s_i - x[i]) < t_i) #endif { if (k_i < 0) { z_i = k_i * s_i + t_i; } else if (k_i > 0) { z_i = k_i * s_i - t_i; } else { z_i = 0.0; } q_n += c_r * d_i * z_i * z_i; } else { q_n += d_i * x[i] * x[i]; } } funevals = funevals + 1; #if ASA_TEMPLATE_SAVE /* cause a crash */ if ((ptr_read_test = fopen ("asa_save", "r")) == NULL) { read_test = 1; fclose (ptr_read_test); } else { fclose (ptr_read_test); } /* will need a few hundred if testing ASA_PARALLEL to get an asa_save */ if (funevals == 50 && read_test == 1) { fprintf (ptr_out, "\n\n*** intended crash to test ASA_SAVE *** \n\n"); fflush (ptr_out); #if INCL_STDOUT printf ("\n\n*** intended crash to test ASA_SAVE *** \n\n"); #endif /* INCL_STDOUT */ exit (2); } #endif *cost_flag = TRUE; #if SELF_OPTIMIZE #else #if TIME_CALC /* print the time every PRINT_FREQUENCY evaluations */ if ((PRINT_FREQUENCY > 0) && ((funevals % PRINT_FREQUENCY) == 0)) { fprintf (ptr_out, "funevals = %ld ", funevals); #if INCL_STDOUT print_time ("", ptr_out); #endif /* INCL_STDOUT */ } #endif #endif #if ASA_TEMPLATE_SAMPLE USER_OPTIONS->Cost_Acceptance_Flag = TRUE; if (USER_OPTIONS->User_Acceptance_Flag == FALSE && *cost_flag == TRUE) USER_OPTIONS->Acceptance_Test (q_n, parameter_lower_bound, parameter_upper_bound, *parameter_dimension, USER_OPTIONS); #endif /* ASA_TEMPLATE_SAMPLE */ return (q_n); #endif /* ASA_TEST */ #if ASA_TEMPLATE_SAMPLE int n; double cost; if (*cost_flag == FALSE) { for (n = 0; n < *parameter_dimension; ++n) if (fabs (parameter_upper_bound[n] - parameter_lower_bound[n]) < (double) EPS_DOUBLE) continue; cost_tangents[n] = 2.0 * x[n]; } cost = 0.0; for (n = 0; n < *parameter_dimension; ++n) { if (fabs (parameter_upper_bound[n] - parameter_lower_bound[n]) < (double) EPS_DOUBLE) continue; cost += (x[n] * x[n]); } *cost_flag = TRUE; USER_OPTIONS->Cost_Acceptance_Flag = TRUE; if (USER_OPTIONS->User_Acceptance_Flag == FALSE && *cost_flag == TRUE) USER_OPTIONS->Acceptance_Test (cost, parameter_lower_bound, parameter_upper_bound, *parameter_dimension, USER_OPTIONS); return (cost); #endif /* ASA_TEMPLATE_SAMPLE */ #if MY_TEMPLATE /* MY_TEMPLATE_cost */ return USER_OPTIONS->Asa_Data_Ptr(x,cost_flag); /* Use the parameter values x[] and define your cost_function. The {} brackets around this function are already in place. */ #endif /* MY_TEMPLATE cost */ } #endif /* COST_FILE */ /* Here is a good random number generator */ #define MULT ((LONG_INT) 25173) #define MOD ((LONG_INT) 65536) #define INCR ((LONG_INT) 13849) #define FMOD ((double) 65536.0) /*********************************************************************** * double myrand - returns random number between 0 and 1 * This routine returns the random number generator between 0 and 1 ***********************************************************************/ static double myrand (LONG_INT * rand_seed) /* returns random number in {0,1} */ { #if TRUE /* (change to FALSE for alternative RNG) */ *rand_seed = (LONG_INT) ((MULT * (*rand_seed) + INCR) % MOD); return ((double) (*rand_seed) / FMOD); #else /* See "Random Number Generators: Good Ones Are Hard To Find," Park & Miller, CACM 31 (10) (October 1988) pp. 1192-1201. *********************************************************** THIS IMPLEMENTATION REQUIRES AT LEAST 32 BIT INTEGERS *********************************************************** */ #define _A_MULTIPLIER 16807L #define _M_MODULUS 2147483647L /* (2**31)-1 */ #define _Q_QUOTIENT 127773L /* 2147483647 / 16807 */ #define _R_REMAINDER 2836L /* 2147483647 % 16807 */ long lo; long hi; long test; hi = *rand_seed / _Q_QUOTIENT; lo = *rand_seed % _Q_QUOTIENT; test = _A_MULTIPLIER * lo - _R_REMAINDER * hi; if (test > 0) { *rand_seed = test; } else { *rand_seed = test + _M_MODULUS; } return ((double) *rand_seed / _M_MODULUS); #endif /* alternative RNG */ } /*********************************************************************** * double randflt ***********************************************************************/ static double randflt (LONG_INT * rand_seed) { return (resettable_randflt (rand_seed, 0)); } /*********************************************************************** * double resettable_randflt ***********************************************************************/ static double resettable_randflt (LONG_INT * rand_seed, int reset) /* shuffles random numbers in random_array[SHUFFLE] array */ { /* This RNG is a modified algorithm of that presented in * %A K. Binder * %A D. Stauffer * %T A simple introduction to Monte Carlo simulations and some * specialized topics * %B Applications of the Monte Carlo Method in statistical physics * %E K. Binder * %I Springer-Verlag * %C Berlin * %D 1985 * %P 1-36 * where it is stated that such algorithms have been found to be * quite satisfactory in many statistical physics applications. */ double rranf; unsigned kranf; int n; static int initial_flag = 0; LONG_INT initial_seed; static double random_array[SHUFFLE]; /* random variables */ if (*rand_seed < 0) *rand_seed = -*rand_seed; if ((initial_flag == 0) || reset) { initial_seed = *rand_seed; for (n = 0; n < SHUFFLE; ++n) random_array[n] = myrand (&initial_seed); initial_flag = 1; for (n = 0; n < 1000; ++n) /* warm up random generator */ rranf = randflt (&initial_seed); rranf = randflt (rand_seed); return (rranf); } kranf = (unsigned) (myrand (rand_seed) * SHUFFLE) % SHUFFLE; rranf = *(random_array + kranf); *(random_array + kranf) = myrand (rand_seed); return (rranf); } #if USER_COST_SCHEDULE #if HAVE_ANSI double user_cost_schedule (double test_temperature, USER_DEFINES * USER_OPTIONS) #else double user_cost_schedule (test_temperature, USER_OPTIONS) double test_temperature; USER_DEFINES *USER_OPTIONS; #endif /* HAVE_ANSI */ { double x; #if ASA_TEMPLATE_SAMPLE x = F_POW (test_temperature, 0.15); #endif #if ASA_TEMPLATE x = test_temperature; #endif return (x); } #endif /* USER_COST_SCHEDULE */ #if USER_ACCEPTANCE_TEST #if HAVE_ANSI void user_acceptance_test (double current_cost, double *parameter_lower_bound, double *parameter_upper_bound, ALLOC_INT * parameter_dimension, USER_DEFINES * USER_OPTIONS) #else void user_acceptance_test (current_cost, parameter_lower_bound, parameter_upper_bound, parameter_dimension, USER_OPTIONS) double current_cost; double *parameter_lower_bound; double *parameter_upper_bound; ALLOC_INT *parameter_dimension; USER_DEFINES *USER_OPTIONS; #endif /* HAVE_ANSI */ { double uniform_test, curr_cost_temp; #if USER_ACCEPT_ASYMP_EXP double x, q, delta_cost; #endif #if ASA_TEMPLATE /* ASA cost index */ /* Calculate the current ASA cost index. This could be useful to define a new schedule for the cost temperature, beyond simple changes that can be made using USER_COST_SCHEDULE. */ int index; double k_temperature, quench, y; double xparameter_dimension; #if QUENCH_COST quench = USER_OPTIONS->User_Quench_Cost_Scale[0]; #else quench = 1.0; #endif /* QUENCH_COST */ xparameter_dimension = (double) *parameter_dimension; for (index = 0; index < *parameter_dimension; ++index) if (fabs (parameter_upper_bound[index] - parameter_lower_bound[index]) < (double) EPS_DOUBLE) *xparameter_dimension -= 1.0; y = -F_LOG (USER_OPTIONS->Cost_Temp_Curr / USER_OPTIONS->Cost_Temp_Init) / USER_OPTIONS->Cost_Temp_Scale; k_temperature = F_POW (y, xparameter_dimension / quench); #endif /* ASA cost index */ uniform_test = randflt (USER_OPTIONS->Random_Seed); curr_cost_temp = USER_OPTIONS->Cost_Temp_Curr; #if ASA_TEMPLATE #if USER_COST_SCHEDULE curr_cost_temp = (USER_OPTIONS->Cost_Schedule (USER_OPTIONS->Cost_Temp_Curr, USER_OPTIONS) + (double) EPS_DOUBLE); #else curr_cost_temp = USER_OPTIONS->Cost_Temp_Curr; #endif #endif /* ASA_TEMPLATE */ /* You must add in your own test here. If USER_ACCEPT_ASYMP_EXP also is TRUE here, then you can use the default Asymp_Exp_Param=1 to replicate the code in asa.c. */ #if USER_ACCEPT_ASYMP_EXP #if USER_COST_SCHEDULE curr_cost_temp = (USER_OPTIONS->Cost_Schedule (USER_OPTIONS->Cost_Temp_Curr, USER_OPTIONS) + (double) EPS_DOUBLE); #endif delta_cost = (current_cost - *(USER_OPTIONS->Last_Cost)) / (curr_cost_temp + (double) EPS_DOUBLE); /* The following asymptotic approximation to the exponential * function, "Tsallis statistics," was proposed in * %A T.J.P. Penna * %T Traveling salesman problem and Tsallis statistics * %J Phys. Rev. E * %V 50 * %N 6 * %P R1-R3 * %D 1994 * While the use of the TSP for a test case is of dubious value (since * there are many special algorithms for this problem), the use of this * function is another example of how to control the rate of annealing * of the acceptance criteria. E.g., if you require a more moderate * acceptance test, then negative q may be helpful. */ q = USER_OPTIONS->Asymp_Exp_Param; if (fabs (1.0 - q) < (double) EPS_DOUBLE) x = MIN (1.0, (F_EXP (-delta_cost))); /* Boltzmann test */ else if ((1.0 - (1.0 - q) * delta_cost) < (double) EPS_DOUBLE) x = MIN (1.0, (F_EXP (-delta_cost))); /* Boltzmann test */ else x = MIN (1.0, F_POW ((1.0 - (1.0 - q) * delta_cost), (1.0 / (1.0 - q)))); USER_OPTIONS->Prob_Bias = x; if (x >= uniform_test) USER_OPTIONS->User_Acceptance_Flag = TRUE; else USER_OPTIONS->User_Acceptance_Flag = FALSE; #endif /* USER_ACCEPT_ASYMP_EXP */ } #endif /* USER_ACCEPTANCE_TEST */ #if USER_GENERATING_FUNCTION #if HAVE_ANSI double user_generating_distrib (LONG_INT * seed, ALLOC_INT * parameter_dimension, ALLOC_INT index_v, double temperature_v, double init_param_temp_v, double temp_scale_params_v, double parameter_v, double parameter_range_v, double *last_saved_parameter, USER_DEFINES * USER_OPTIONS) #else double user_generating_distrib (seed, parameter_dimension, index_v, temperature_v, init_param_temp_v, temp_scale_params_v, parameter_v, parameter_range_v, last_saved_parameter, USER_OPTIONS) LONG_INT *seed; ALLOC_INT *parameter_dimension; ALLOC_INT index_v; double temperature_v; double init_param_temp_v; double temp_scale_params_v; double parameter_v; double parameter_range_v; double *last_saved_parameter; USER_DEFINES *USER_OPTIONS; #endif { #if ASA_TEMPLATE double x, y, z; /* This is the ASA distribution. A slower temperature schedule can be obtained here, e.g., temperature_v = pow(temperature_v, 0.5); */ x = randflt (seed); y = x < 0.5 ? -1.0 : 1.0; z = y * temperature_v * (F_POW ((1.0 + 1.0 / temperature_v), fabs (2.0 * x - 1.0)) - 1.0); x = parameter_v + z * parameter_range_v; return (x); #endif /* ASA_TEMPLATE */ } #endif /* USER_GENERATING_FUNCTION */ #if USER_REANNEAL_COST #if HAVE_ANSI int user_reanneal_cost (double *cost_best, double *cost_last, double *initial_cost_temperature, double *current_cost_temperature, USER_DEFINES * USER_OPTIONS) #else int user_reanneal_cost (cost_best, cost_last, initial_cost_temperature, current_cost_temperature, USER_OPTIONS) double *cost_best; double *cost_last; double *initial_cost_temperature; double *current_cost_temperature; USER_DEFINES *USER_OPTIONS; #endif /* HAVE_ANSI */ { int cost_test; double tmp_dbl; #if ASA_TEMPLATE static int first_time = 1; static double save_last[3]; double average_cost_last; if (first_time == 1) { first_time = 0; save_last[0] = save_last[1] = save_last[2] = *cost_last; } save_last[2] = save_last[1]; save_last[1] = save_last[0]; save_last[0] = *cost_last; average_cost_last = fabs ((save_last[0] + save_last[1] + save_last[2]) / 3.0); tmp_dbl = MAX (fabs (*cost_best), average_cost_last); tmp_dbl = MAX ((double) EPS_DOUBLE, tmp_dbl); *initial_cost_temperature = MIN (*initial_cost_temperature, tmp_dbl); /* This test can be useful if your cost function goes from a positive to a negative value, and you do not want to get get stuck in a local minima around zero due to the default in reanneal(). Pick any number instead of 0.0001 */ tmp_dbl = MIN (fabs (*cost_last), fabs (*cost_best)); if (tmp_dbl < 0.0001) cost_test = FALSE; else cost_test = TRUE; #endif /* ASA_TEMPLATE */ tmp_dbl = MAX (fabs (cost_last), fabs (cost_best)); tmp_dbl = MAX ((double) EPS_DOUBLE, tmp_dbl); *initial_cost_temperature = MIN (*initial_cost_temperature, tmp_dbl); *current_cost_temperature = MAX (fabs (cost_last - cost_best), *current_cost_temperature); *current_cost_temperature = MAX ((double) EPS_DOUBLE, *current_cost_temperature); *current_cost_temperature = MIN (*current_cost_temperature, *initial_cost_temperature); cost_test = TRUE; return (cost_test); } #endif /* USER_REANNEAL_COST */ #if USER_REANNEAL_PARAMETERS #if HAVE_ANSI double user_reanneal_params (double current_temp, double tangent, double max_tangent, USER_DEFINES * USER_OPTIONS) #else double user_reanneal_params (current_temp, tangent, max_tangent, USER_OPTIONS) double current_temp; double tangent; double max_tangent; USER_DEFINES *USER_OPTIONS; #endif /* HAVE_ANSI */ { #if ASA_TEMPLATE double x; x = current_temp * (max_tangent / tangent); return (x); #endif } #endif /* USER_REANNEAL_PARAMETERS */ #if SELF_OPTIMIZE /*********************************************************************** * main * This is a sample calling program to self-optimize ASA ***********************************************************************/ #if HAVE_ANSI #if ASA_LIB int asa_main ( #if ASA_TEMPLATE_LIB double *main_recur_cost_value, double *main_recur_cost_parameters, int *main_recur_exit_code #endif ) #else /* ASA_LIB */ int main (int argc, char **argv) #endif /* ASA_LIB */ #else /* HAVE_ANSI */ #if ASA_LIB int asa_main ( #if ASA_TEMPLATE_LIB main_recur_cost_value, main_recur_cost_parameters, main_recur_exit_code #endif ) #if ASA_TEMPLATE_LIB double *main_recur_cost_value; double *main_recur_cost_parameters; int *main_recur_exit_code; #endif #else /* ASA_LIB */ int main (argc, argv) int argc; char **argv; #endif /* ASA_LIB */ #endif /* HAVE_ANSI */ { /* seed for random number generator */ LONG_INT *recur_rand_seed; #if RECUR_OPTIONS_FILE FILE *recur_ptr_options; char read_option[80]; char read_if[4], read_FALSE[6], read_comm1[3], read_ASA_SAVE[9], read_comm2[3]; int read_int; #if INT_LONG LONG_INT read_long; #endif double read_double; #endif int *recur_exit_code; #if MULTI_MIN int multi_index; ALLOC_INT n_param; #endif double *recur_parameter_lower_bound, *recur_parameter_upper_bound; double *recur_cost_parameters, *recur_cost_tangents, *recur_cost_curvature; double recur_cost_value; ALLOC_INT *recur_parameter_dimension; int *recur_parameter_int_real; int *recur_cost_flag; int recur_initialize_parameters_value; ALLOC_INT recur_v; #if MY_TEMPLATE /* MY_TEMPLATE_recur_main_decl */ /* add some declarations if required */ #endif USER_DEFINES *RECUR_USER_OPTIONS; if ((recur_parameter_dimension = (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_parameter_dimension"); Exit_USER (user_exit_msg); return (-2); } if ((recur_exit_code = (int *) calloc (1, sizeof (int))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_exit_code"); Exit_USER (user_exit_msg); return (-2); } if ((recur_cost_flag = (int *) calloc (1, sizeof (int))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_cost_flag"); Exit_USER (user_exit_msg); return (-2); } if ((RECUR_USER_OPTIONS = (USER_DEFINES *) calloc (1, sizeof (USER_DEFINES))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): RECUR_USER_OPTIONS"); Exit_USER (user_exit_msg); return (-2); } #if RECUR_OPTIONS_FILE recur_ptr_options = fopen ("asa_opt_recur", "r"); fscanf (recur_ptr_options, "%s%s%s%s%s", read_if, read_FALSE, read_comm1, read_ASA_SAVE, read_comm2); if (strcmp (read_if, "#if") || strcmp (read_FALSE, "FALSE") || strcmp (read_comm1, "/*") || strcmp (read_ASA_SAVE, "ASA_SAVE") || strcmp (read_comm2, "*/")) { fprintf (ptr_out, "\n\n*** not asa_opt_recur for this version *** \n\n"); fflush (ptr_out); #if INCL_STDOUT printf ("\n\n*** EXIT not asa_opt_recur for this version *** \n\n"); #endif /* INCL_STDOUT */ return (-6); } #if INT_LONG fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%ld", &read_long); RECUR_USER_OPTIONS->Limit_Acceptances = read_long; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%ld", &read_long); RECUR_USER_OPTIONS->Limit_Generated = read_long; #else fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Limit_Acceptances = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Limit_Generated = read_int; #endif fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Limit_Invalid_Generated_States = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%lf", &read_double); RECUR_USER_OPTIONS->Accepted_To_Generated_Ratio = read_double; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%lf", &read_double); RECUR_USER_OPTIONS->Cost_Precision = read_double; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Maximum_Cost_Repeat = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Number_Cost_Samples = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%lf", &read_double); RECUR_USER_OPTIONS->Temperature_Ratio_Scale = read_double; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%lf", &read_double); RECUR_USER_OPTIONS->Cost_Parameter_Scale_Ratio = read_double; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%lf", &read_double); RECUR_USER_OPTIONS->Temperature_Anneal_Scale = read_double; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Include_Integer_Parameters = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->User_Initial_Parameters = read_int; #if INT_ALLOC fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Sequential_Parameters = read_int; #else #if INT_LONG fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%ld", &read_long); RECUR_USER_OPTIONS->Sequential_Parameters = read_long; #else fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Sequential_Parameters = read_int; #endif #endif fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%lf", &read_double); RECUR_USER_OPTIONS->Initial_Parameter_Temperature = read_double; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Acceptance_Frequency_Modulus = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Generated_Frequency_Modulus = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Reanneal_Cost = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Reanneal_Parameters = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%lf", &read_double); RECUR_USER_OPTIONS->Delta_X = read_double; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->User_Tangents = read_int; fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%d", &read_int); RECUR_USER_OPTIONS->Curvature_0 = read_int; #else /* RECUR_OPTIONS_FILE */ RECUR_USER_OPTIONS->Limit_Acceptances = 100; RECUR_USER_OPTIONS->Limit_Generated = 1000; RECUR_USER_OPTIONS->Limit_Invalid_Generated_States = 1000; RECUR_USER_OPTIONS->Accepted_To_Generated_Ratio = 1.0E-4; RECUR_USER_OPTIONS->Cost_Precision = 1.0E-18; RECUR_USER_OPTIONS->Maximum_Cost_Repeat = 2; RECUR_USER_OPTIONS->Number_Cost_Samples = 2; RECUR_USER_OPTIONS->Temperature_Ratio_Scale = 1.0E-5; RECUR_USER_OPTIONS->Cost_Parameter_Scale_Ratio = 1.0; RECUR_USER_OPTIONS->Temperature_Anneal_Scale = 100.0; RECUR_USER_OPTIONS->Include_Integer_Parameters = FALSE; RECUR_USER_OPTIONS->User_Initial_Parameters = FALSE; RECUR_USER_OPTIONS->Sequential_Parameters = -1; RECUR_USER_OPTIONS->Initial_Parameter_Temperature = 1.0; RECUR_USER_OPTIONS->Acceptance_Frequency_Modulus = 15; RECUR_USER_OPTIONS->Generated_Frequency_Modulus = 10000; RECUR_USER_OPTIONS->Reanneal_Cost = FALSE; RECUR_USER_OPTIONS->Reanneal_Parameters = FALSE; RECUR_USER_OPTIONS->Delta_X = 1.0E-6; RECUR_USER_OPTIONS->User_Tangents = FALSE; RECUR_USER_OPTIONS->Curvature_0 = TRUE; #endif /* RECUR_OPTIONS_FILE */ /* the number of parameters for the recur_cost_function */ #if RECUR_OPTIONS_FILE_DATA fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%s", read_option); #if INT_ALLOC fscanf (recur_ptr_options, "%d", &read_int); *recur_parameter_dimension = read_int; #else #if INT_LONG fscanf (recur_ptr_options, "%ld", &read_long); *recur_parameter_dimension = read_long; #else fscanf (recur_ptr_options, "%d", &read_int); *recur_parameter_dimension = read_int; #endif #endif #else /* RECUR_OPTIONS_FILE_DATA */ #if ASA_TEMPLATE_SELFOPT *recur_parameter_dimension = 2; #endif #if MY_TEMPLATE /* MY_TEMPLATE_recur_dim */ /* If not using RECUR_OPTIONS_FILE_DATA or data read from recur_asa_opt, insert the number of parameters for the recur_cost_function */ #endif /* MY_TEMPLATE recur_dim */ #endif /* RECUR_OPTIONS_FILE_DATA */ if ((recur_parameter_lower_bound = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_parameter_lower_bound"); Exit_USER (user_exit_msg); return (-2); } if ((recur_parameter_upper_bound = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_parameter_upper_bound"); Exit_USER (user_exit_msg); return (-2); } if ((recur_cost_parameters = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_cost_parameters"); Exit_USER (user_exit_msg); return (-2); } if ((recur_parameter_int_real = (int *) calloc (*recur_parameter_dimension, sizeof (int))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_parameter_int_real"); Exit_USER (user_exit_msg); return (-2); } if ((recur_cost_tangents = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_cost_tangents"); Exit_USER (user_exit_msg); return (-2); } if (RECUR_USER_OPTIONS->Curvature_0 == FALSE || RECUR_USER_OPTIONS->Curvature_0 == -1) { if ((recur_cost_curvature = (double *) calloc ((*recur_parameter_dimension) * (*recur_parameter_dimension), sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_cost_curvature"); Exit_USER (user_exit_msg); return (-2); } } else { recur_cost_curvature = (double *) NULL; } #if ASA_TEMPLATE_SELFOPT /* Set memory to that required for use. */ RECUR_USER_OPTIONS->Asa_Data_Dim_Dbl = 1; if ((RECUR_USER_OPTIONS->Asa_Data_Dbl = (double *) calloc (RECUR_USER_OPTIONS->Asa_Data_Dim_Dbl, sizeof (double))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): RECUR_USER_OPTIONS->Asa_Data_Dbl"); Exit_USER (user_exit_msg); return (-2); } /* Use Asa_Data[0] as flag, e.g., if used with SELF_OPTIMIZE. */ RECUR_USER_OPTIONS->Asa_Data_Dbl[0] = 0; #endif /* ASA_TEMPLATE_SELFOPT */ #if OPTIONAL_DATA_PTR #if ASA_TEMPLATE RECUR_USER_OPTIONS->Asa_Data_Dim_Ptr = 1; if ((RECUR_USER_OPTIONS->Asa_Data_Ptr = (OPTIONAL_PTR_TYPE *) calloc (RECUR_USER_OPTIONS->Asa_Data_Dim_Ptr, sizeof (OPTIONAL_PTR_TYPE))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): RECUR_USER_OPTIONS->Asa_Data_Ptr"); Exit_USER (user_exit_msg); return (-2); } #endif /* ASA_TEMPLATE */ #endif /* OPTIONAL_DATA_PTR */ #if ASA_SAVE /* Such data could be saved in a user_save file, but for convenience here everything is saved in asa_save. */ RECUR_USER_OPTIONS->Random_Array_Dim = SHUFFLE; RECUR_USER_OPTIONS->Random_Array = random_array; #endif /* ASA_SAVE */ /* open the output file */ #if ASA_SAVE if (!strcmp (USER_OUT, "STDOUT")) { #if INCL_STDOUT ptr_out = stdout; #endif /* INCL_STDOUT */ } else { ptr_out = fopen (USER_OUT, "a"); } #else if (!strcmp (USER_OUT, "STDOUT")) { #if INCL_STDOUT ptr_out = stdout; #endif /* INCL_STDOUT */ } else { ptr_out = fopen (USER_OUT, "w"); } #endif // fprintf (ptr_out, "%s\n\n", USER_ID); fflush (ptr_out); if ((recur_rand_seed = (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): recur_rand_seed"); Exit_USER (user_exit_msg); return (-2); } /* first value of *recur_rand_seed */ #if ASA_LIB *recur_rand_seed = (asa_rand_seed ? *asa_rand_seed : (LONG_INT) 696969); #else *recur_rand_seed = 696969; #endif randflt (recur_rand_seed); #if USER_COST_SCHEDULE RECUR_USER_OPTIONS->Cost_Schedule = recur_user_cost_schedule; #endif #if USER_ACCEPTANCE_TEST RECUR_USER_OPTIONS->Acceptance_Test = recur_user_acceptance_test; #endif #if USER_ACCEPT_ASYMP_EXP RECUR_USER_OPTIONS->Asymp_Exp_Param = 1.0; #endif #if USER_GENERATING_FUNCTION RECUR_USER_OPTIONS->Generating_Distrib = recur_user_generating_distrib; #endif #if USER_REANNEAL_COST RECUR_USER_OPTIONS->Reanneal_Cost_Function = recur_user_reanneal_cost; #endif #if USER_REANNEAL_PARAMETERS RECUR_USER_OPTIONS->Reanneal_Params_Function = recur_user_reanneal_params; #endif #if MY_TEMPLATE /* MY_TEMPLATE_recur_pre_initialize */ /* last changes before entering recur_initialize_parameters() */ #endif /* initialize the users parameters, allocating space, etc. Note that the default is to have asa generate the initial recur_cost_parameters that satisfy the user's constraints. */ recur_initialize_parameters_value = recur_initialize_parameters (recur_cost_parameters, recur_parameter_lower_bound, recur_parameter_upper_bound, recur_cost_tangents, recur_cost_curvature, recur_parameter_dimension, recur_parameter_int_real, #if RECUR_OPTIONS_FILE_DATA recur_ptr_options, #endif RECUR_USER_OPTIONS); #if RECUR_OPTIONS_FILE fclose (recur_ptr_options); #endif if (recur_initialize_parameters_value == -2) return (recur_initialize_parameters_value); #if USER_ASA_OUT if ((RECUR_USER_OPTIONS->Asa_Out_File = (char *) calloc (80, sizeof (char))) == NULL) { strcpy (user_exit_msg, "main()/asa_main(): RECUR_USER_OPTIONS->Asa_Out_File"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE_SELFOPT strcpy (RECUR_USER_OPTIONS->Asa_Out_File, "asa_sfop"); #endif #endif recur_cost_value = asa (RECUR_USER_COST_FUNCTION, randflt, recur_rand_seed, recur_cost_parameters, recur_parameter_lower_bound, recur_parameter_upper_bound, recur_cost_tangents, recur_cost_curvature, recur_parameter_dimension, recur_parameter_int_real, recur_cost_flag, recur_exit_code, RECUR_USER_OPTIONS); if (*recur_exit_code == -1) { #if INCL_STDOUT printf ("\n\n*** error in calloc in ASA ***\n\n"); #endif /* INCL_STDOUT */ fprintf (ptr_out, "\n\n*** error in calloc in ASA ***\n\n"); fflush (ptr_out); return (-1); } #if MULTI_MIN fprintf (ptr_out, "Multi_Specify = %d\n", RECUR_USER_OPTIONS->Multi_Specify); for (n_param = 0; n_param < *recur_parameter_dimension; ++n_param) { fprintf (ptr_out, #if INT_ALLOC "Multi_Grid[%d] = %12.7g\n", #else #if INT_LONG "Multi_Grid[%ld] = %12.7g\n", #else "Multi_Grid[%d] = %12.7g\n", #endif #endif n_param, RECUR_USER_OPTIONS->Multi_Grid[n_param]); } fprintf (ptr_out, "\n"); for (multi_index = 0; multi_index < RECUR_USER_OPTIONS->Multi_Number; ++multi_index) { fprintf (ptr_out, "\n"); fprintf (ptr_out, "Multi_Cost[%d] = %12.7g\n", multi_index, RECUR_USER_OPTIONS->Multi_Cost[multi_index]); for (n_param = 0; n_param < *recur_parameter_dimension; ++n_param) { fprintf (ptr_out, #if INT_ALLOC "Multi_Params[%d][%d] = %12.7g\n", #else #if INT_LONG "Multi_Params[%d][%ld] = %12.7g\n", #else "Multi_Params[%d][%d] = %12.7g\n", #endif #endif multi_index, n_param, RECUR_USER_OPTIONS->Multi_Params[multi_index][n_param]); } } fprintf (ptr_out, "\n"); fflush (ptr_out); #endif /* MULTI_MIN */ #if FITLOC /* Fit_Local and Penalty may be set adaptively */ RECUR_USER_OPTIONS->Penalty = 1000; RECUR_USER_OPTIONS->Fit_Local = 0; RECUR_USER_OPTIONS->Iter_Max = 500; if (RECUR_USER_OPTIONS->Fit_Local >= 1) { recur_cost_value = fitloc (RECUR_USER_COST_FUNCTION, recur_cost_parameters, recur_parameter_lower_bound, recur_parameter_upper_bound, recur_cost_tangents, recur_cost_curvature, recur_parameter_dimension, recur_parameter_int_real, recur_cost_flag, recur_exit_code, RECUR_USER_OPTIONS, ptr_out); } #endif /* FITLOC */ fprintf (ptr_out, "\n\n recur_cost_value = %12.7g\n", recur_cost_value); #if MY_TEMPLATE /* MY_TEMPLATE_recur_post_recur_asa */ #endif #if ASA_TEMPLATE_LIB *main_recur_cost_value = recur_cost_value; for (recur_v = 0; recur_v < *recur_parameter_dimension; ++recur_v) { main_recur_cost_parameters[recur_v] = recur_cost_parameters[recur_v]; } *main_recur_exit_code = *recur_exit_code; #endif for (recur_v = 0; recur_v < *recur_parameter_dimension; ++recur_v) #if INT_ALLOC fprintf (ptr_out, "recur_cost_parameters[%d] = %12.7g\n", #else #if INT_LONG fprintf (ptr_out, "recur_cost_parameters[%ld] = %12.7g\n", #else fprintf (ptr_out, "recur_cost_parameters[%d] = %12.7g\n", #endif #endif recur_v, recur_cost_parameters[recur_v]); fprintf (ptr_out, "\n\n"); /* close all files */ ptr_out != stdout && fclose (ptr_out); #if OPTIONAL_DATA_DBL free (RECUR_USER_OPTIONS->Asa_Data_Dbl); #endif #if OPTIONAL_DATA_INT free (RECUR_USER_OPTIONS->Asa_Data_Int); #endif #if OPTIONAL_DATA_PTR free (RECUR_USER_OPTIONS->Asa_Data_Ptr); #endif #if USER_ASA_OUT free (RECUR_USER_OPTIONS->Asa_Out_File); #endif #if ASA_SAMPLE free (RECUR_USER_OPTIONS->Bias_Generated); #endif #if ASA_QUEUE #if ASA_RESOLUTION #else free (RECUR_USER_OPTIONS->Queue_Resolution); #endif #endif #if ASA_RESOLUTION free (RECUR_USER_OPTIONS->Coarse_Resolution); #endif if (RECUR_USER_OPTIONS->Curvature_0 == FALSE || RECUR_USER_OPTIONS->Curvature_0 == -1) free (recur_cost_curvature); #if USER_INITIAL_PARAMETERS_TEMPS free (RECUR_USER_OPTIONS->User_Parameter_Temperature); #endif #if USER_INITIAL_COST_TEMP free (RECUR_USER_OPTIONS->User_Cost_Temperature); #endif #if DELTA_PARAMETERS free (RECUR_USER_OPTIONS->User_Delta_Parameter); #endif #if QUENCH_PARAMETERS free (RECUR_USER_OPTIONS->User_Quench_Param_Scale); #endif #if QUENCH_COST free (RECUR_USER_OPTIONS->User_Quench_Cost_Scale); #endif #if RATIO_TEMPERATURE_SCALES free (RECUR_USER_OPTIONS->User_Temperature_Ratio); #endif #if MULTI_MIN free (RECUR_USER_OPTIONS->Multi_Cost); free (RECUR_USER_OPTIONS->Multi_Grid); for (multi_index = 0; multi_index < RECUR_USER_OPTIONS->Multi_Number; ++multi_index) { free (RECUR_USER_OPTIONS->Multi_Params[multi_index]); } free (RECUR_USER_OPTIONS->Multi_Params); #endif /* MULTI_MIN */ free (RECUR_USER_OPTIONS); free (recur_parameter_dimension); free (recur_exit_code); free (recur_cost_flag); free (recur_parameter_lower_bound); free (recur_parameter_upper_bound); free (recur_cost_parameters); free (recur_parameter_int_real); free (recur_cost_tangents); free (recur_rand_seed); return (0); /* NOTREACHED */ } /*********************************************************************** * recur_initialize_parameters * This depends on the users cost function to optimize (minimum). * The routine allocates storage needed for asa. The user should * define the number of parameters and their ranges, * and make sure the initial parameters are within * the minimum and maximum ranges. The array * recur_parameter_int_real should be REAL_TYPE (-1) * for real parameters, ***********************************************************************/ #if HAVE_ANSI int recur_initialize_parameters (double *recur_cost_parameters, double *recur_parameter_lower_bound, double *recur_parameter_upper_bound, double *recur_cost_tangents, double *recur_cost_curvature, ALLOC_INT * recur_parameter_dimension, int *recur_parameter_int_real, #if RECUR_OPTIONS_FILE_DATA FILE * recur_ptr_options, #endif USER_DEFINES * RECUR_USER_OPTIONS) #else int recur_initialize_parameters (recur_cost_parameters, recur_parameter_lower_bound, recur_parameter_upper_bound, recur_cost_tangents, recur_cost_curvature, recur_parameter_dimension, recur_parameter_int_real, #if RECUR_OPTIONS_FILE_DATA recur_ptr_options, #endif RECUR_USER_OPTIONS) double *recur_parameter_lower_bound; double *recur_parameter_upper_bound; double *recur_cost_parameters; double *recur_cost_tangents; double *recur_cost_curvature; ALLOC_INT *recur_parameter_dimension; int *recur_parameter_int_real; #if RECUR_OPTIONS_FILE_DATA FILE *recur_ptr_options; #endif USER_DEFINES *RECUR_USER_OPTIONS; #endif { ALLOC_INT index; #if RECUR_OPTIONS_FILE_DATA char read_option[80]; ALLOC_INT read_index; #endif #if MY_TEMPLATE /* MY_TEMPLATE_recur_init_decl */ /* add some declarations if required */ #endif #if MULTI_MIN int multi_index; #endif #if RECUR_OPTIONS_FILE_DATA fscanf (recur_ptr_options, "%s", read_option); for (index = 0; index < *recur_parameter_dimension; ++index) { #if MY_TEMPLATE /* MY_TEMPLATE_recur_read_opt */ /* put in some code as required to alter lines read from recur_asa_opt */ #endif #if INT_ALLOC fscanf (recur_ptr_options, "%d", &read_index); #else #if INT_LONG fscanf (recur_ptr_options, "%ld", &read_index); #else fscanf (recur_ptr_options, "%d", &read_index); #endif #endif fscanf (recur_ptr_options, "%lf%lf%lf%d", &(recur_parameter_lower_bound[read_index]), &(recur_parameter_upper_bound[read_index]), &(recur_cost_parameters[read_index]), &(recur_parameter_int_real[read_index])); } #else /* RECUR_OPTIONS_FILE_DATA */ #if ASA_TEMPLATE_SELFOPT /* NOTE: USER_OPTIONS->Temperature_Ratio_Scale = x[0]; USER_OPTIONS->Cost_Parameter_Scale_Ratio = x[1]; */ /* store the initial parameter values */ recur_cost_parameters[0] = 1.0E-5; recur_cost_parameters[1] = 1.0; recur_parameter_lower_bound[0] = 1.0E-6; recur_parameter_upper_bound[0] = 1.0E-4; recur_parameter_lower_bound[1] = 0.5; recur_parameter_upper_bound[1] = 3.0; /* store the initial parameter types */ for (index = 0; index < *recur_parameter_dimension; ++index) recur_parameter_int_real[index] = REAL_TYPE; #endif #endif /* RECUR_OPTIONS_FILE_DATA */ #if USER_INITIAL_PARAMETERS_TEMPS if ((RECUR_USER_OPTIONS->User_Parameter_Temperature = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->User_Parameter_Temperature"); Exit_USER (user_exit_msg); return (-2); } for (index = 0; index < *recur_parameter_dimension; ++index) RECUR_USER_OPTIONS->User_Parameter_Temperature[index] = 1.0; #endif /* USER_INITIAL_PARAMETERS_TEMPS */ #if USER_INITIAL_COST_TEMP if ((RECUR_USER_OPTIONS->User_Cost_Temperature = (double *) calloc (1, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->User_Cost_Temperature"); Exit_USER (user_exit_msg); return (-2); } RECUR_USER_OPTIONS->User_Cost_Temperature[0] = 5.936648E+09; #endif /* USER_INITIAL_COST_TEMP */ #if DELTA_PARAMETERS if ((RECUR_USER_OPTIONS->User_Delta_Parameter = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->User_Delta_Parameter"); Exit_USER (user_exit_msg); return (-2); } for (index = 0; index < *recur_parameter_dimension; ++index) RECUR_USER_OPTIONS->User_Delta_Parameter[index] = 0.001; #endif /* DELTA_PARAMETERS */ #if QUENCH_PARAMETERS if ((RECUR_USER_OPTIONS->User_Quench_Param_Scale = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->User_Quench_Param_Scale"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE for (index = 0; index < *recur_parameter_dimension; ++index) RECUR_USER_OPTIONS->User_Quench_Param_Scale[index] = 1.0; #endif #endif /* QUENCH_PARAMETERS */ #if QUENCH_COST if ((RECUR_USER_OPTIONS->User_Quench_Cost_Scale = (double *) calloc (1, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->User_Quench_Cost_Scale"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE RECUR_USER_OPTIONS->User_Quench_Cost_Scale[0] = 1.0; #endif #endif /* QUENCH_COST */ /* use asa_opt_recur to read in QUENCH RECUR_USER_OPTIONS */ #if RECUR_OPTIONS_FILE_DATA #if QUENCH_COST fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%lf", &(RECUR_USER_OPTIONS->User_Quench_Cost_Scale[0])); #if QUENCH_PARAMETERS fscanf (recur_ptr_options, "%s", read_option); fscanf (recur_ptr_options, "%s", read_option); for (index = 0; index < *recur_parameter_dimension; ++index) { #if INT_ALLOC fscanf (recur_ptr_options, "%d", &read_index); #else #if INT_LONG fscanf (recur_ptr_options, "%ld", &read_index); #else fscanf (recur_ptr_options, "%d", &read_index); #endif #endif fscanf (recur_ptr_options, "%lf", &(RECUR_USER_OPTIONS->User_Quench_Param_Scale[read_index])); } #endif /* QUENCH_PARAMETERS */ #endif /* QUENCH_COST */ #endif /* RECUR_OPTIONS_FILE_DATA */ #if RATIO_TEMPERATURE_SCALES if ((RECUR_USER_OPTIONS->User_Temperature_Ratio = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->User_Temperature_Ratio"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE for (index = 0; index < *recur_parameter_dimension; ++index) RECUR_USER_OPTIONS->User_Temperature_Ratio[index] = 1.0; #endif #endif /* RATIO_TEMPERATURE_SCALES */ /* Defines the limit of collection of sampled data by asa */ #if ASA_SAMPLE /* create memory for Bias_Generated[] */ if ((RECUR_USER_OPTIONS->Bias_Generated = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->Bias_Generated"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE RECUR_USER_OPTIONS->Limit_Weights = 1.0E-7; #if QUENCH_COST RECUR_USER_OPTIONS->User_Quench_Cost_Scale[0] = 1.0; #endif #if QUENCH_PARAMETERS for (index = 0; index < *recur_parameter_dimension; ++index) RECUR_USER_OPTIONS->User_Quench_Param_Scale[index] = 1.0; #endif #endif /* ASA_TEMPLATE */ #endif /* ASA_SAMPLE */ #if ASA_TEMPLATE #if ASA_PARALLEL RECUR_USER_OPTIONS->Gener_Block = 1; RECUR_USER_OPTIONS->Gener_Block_Max = 1; RECUR_USER_OPTIONS->Gener_Mov_Avr = 1; #endif #endif #if ASA_RESOLUTION if ((RECUR_USER_OPTIONS->Coarse_Resolution = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->Coarse_Resolution"); Exit_USER (user_exit_msg); return (-2); } #endif #if MULTI_MIN #if ASA_TEMPLATE RECUR_USER_OPTIONS->Multi_Number = 2; #endif if ((RECUR_USER_OPTIONS->Multi_Cost = (double *) calloc (RECUR_USER_OPTIONS->Multi_Number, sizeof (double))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): RECUR_USER_OPTIONS->Multi_Cost"); Exit_USER (user_exit_msg); return (-2); } if ((RECUR_USER_OPTIONS->Multi_Grid = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->Multi_Grid"); Exit_USER (user_exit_msg); return (-2); } if ((RECUR_USER_OPTIONS->Multi_Params = (double **) calloc (RECUR_USER_OPTIONS->Multi_Number, sizeof (double *))) == NULL) { strcpy (user_exit_msg, "initialize_parameters(): RECUR_USER_OPTIONS->Multi_Params"); Exit_USER (user_exit_msg); return (-2); } for (multi_index = 0; multi_index < RECUR_USER_OPTIONS->Multi_Number; ++multi_index) { if ((RECUR_USER_OPTIONS->Multi_Params[multi_index] = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->Multi_Params[multi_index]"); Exit_USER (user_exit_msg); return (-2); } } #if ASA_TEST for (index = 0; index < *recur_parameter_dimension; ++index) { RECUR_USER_OPTIONS->Multi_Grid[index] = 0.05; } RECUR_USER_OPTIONS->Multi_Specify = 0; #endif #if ASA_TEMPLATE for (index = 0; index < *recur_parameter_dimension; ++index) { RECUR_USER_OPTIONS->Multi_Grid[index] = (recur_parameter_upper_bound[index] - recur_parameter_lower_bound[index]) / 100.0; } RECUR_USER_OPTIONS->Multi_Specify = 0; #endif /* ASA_TEMPLATE */ #endif /* MULTI_MIN */ #if ASA_TEMPLATE_QUEUE RECUR_USER_OPTIONS->Queue_Size = 0; #endif #if ASA_QUEUE #if ASA_RESOLUTION RECUR_USER_OPTIONS->Queue_Resolution = RECUR_USER_OPTIONS->Coarse_Resolution; #else /* ASA_RESOLUTION */ if ((RECUR_USER_OPTIONS->Queue_Resolution = (double *) calloc (*recur_parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_initialize_parameters(): RECUR_USER_OPTIONS->Queue_Resolution"); Exit_USER (user_exit_msg); return (-2); } #endif /* ASA_RESOLUTION */ #if ASA_TEMPLATE_QUEUE RECUR_USER_OPTIONS->Queue_Size = 0; #endif #endif /* ASA_QUEUE */ #if MY_TEMPLATE /* MY_TEMPLATE_recur_params */ /* If not using RECUR_OPTIONS_FILE_DATA or data read from recur_asa_opt, store the recur_parameter ranges store the recur_parameter types store the initial recur_parameter values other changes needed for initialization */ #endif /* MY_TEMPLATE recur_params */ RECUR_USER_OPTIONS->Asa_Recursive_Level = 1; return (0); } /*********************************************************************** * double recur_cost_function * This is the users cost function to optimize * (find the minimum). * cost_flag is set to TRUE if the parameter set * does not violates any constraints * recur_parameter_lower_bound and recur_parameter_upper_bound * may be adaptively changed during the search. ***********************************************************************/ #if HAVE_ANSI double recur_cost_function (double *x, double *recur_parameter_lower_bound, double *recur_parameter_upper_bound, double *recur_cost_tangents, double *recur_cost_curvature, ALLOC_INT * recur_parameter_dimension, int *recur_parameter_int_real, int *recur_cost_flag, int *recur_exit_code, USER_DEFINES * RECUR_USER_OPTIONS) #else double recur_cost_function (x, recur_parameter_lower_bound, recur_parameter_upper_bound, recur_cost_tangents, recur_cost_curvature, recur_parameter_dimension, recur_parameter_int_real, recur_cost_flag, recur_exit_code, RECUR_USER_OPTIONS) double *x; double *recur_parameter_lower_bound; double *recur_parameter_upper_bound; double *recur_cost_tangents; double *recur_cost_curvature; ALLOC_INT *recur_parameter_dimension; int *recur_parameter_int_real; int *recur_cost_flag; int *recur_exit_code; USER_DEFINES *RECUR_USER_OPTIONS; #endif { double cost_value; static LONG_INT recur_funevals = 0; int *exit_code; #if OPTIONAL_DATA_PTR int data_ptr_flg; #endif #if OPTIONS_FILE FILE *ptr_options; char read_option[80]; char read_if[4], read_FALSE[6], read_comm1[3], read_ASA_SAVE[9], read_comm2[3]; int read_int; #if INT_LONG LONG_INT read_long; #endif double read_double; #endif #if MY_TEMPLATE /* MY_TEMPLATE_recur_cost_decl */ /* add some declarations if required */ #endif double *parameter_lower_bound, *parameter_upper_bound; double *cost_parameters; double *cost_tangents, *cost_curvature; ALLOC_INT *parameter_dimension; int *parameter_int_real; int *cost_flag; static LONG_INT *rand_seed; static int initial_flag = 0; #if MULTI_MIN int multi_index; #endif USER_DEFINES *USER_OPTIONS; recur_funevals = recur_funevals + 1; if ((rand_seed = (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): rand_seed"); Exit_USER (user_exit_msg); return (-2); } if ((USER_OPTIONS = (USER_DEFINES *) calloc (1, sizeof (USER_DEFINES))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): USER_OPTIONS"); Exit_USER (user_exit_msg); return (-2); } #if OPTIONS_FILE /* Test to see if asa_opt is in correct directory. This is useful for some PC and Mac compilers. */ if ((ptr_options = fopen ("asa_opt", "r")) == NULL) { fprintf (ptr_out, "\n\n*** fopen asa_opt failed *** \n\n"); fflush (ptr_out); #if INCL_STDOUT printf ("\n\n*** EXIT fopen asa_opt failed *** \n\n"); #endif /* INCL_STDOUT */ return (6); } fscanf (ptr_options, "%s%s%s%s%s", read_if, read_FALSE, read_comm1, read_ASA_SAVE, read_comm2); if (strcmp (read_if, "#if") || strcmp (read_FALSE, "FALSE") || strcmp (read_comm1, "/*") || strcmp (read_ASA_SAVE, "ASA_SAVE") || strcmp (read_comm2, "*/")) { fprintf (ptr_out, "\n\n*** not asa_opt for this version *** \n\n"); fflush (ptr_out); #if INCL_STDOUT printf ("\n\n*** EXIT not asa_opt for this version *** \n\n"); #endif /* INCL_STDOUT */ return (-6); } #if INT_LONG fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%ld", &read_long); USER_OPTIONS->Limit_Acceptances = read_long; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%ld", &read_long); USER_OPTIONS->Limit_Generated = read_long; #else fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Limit_Acceptances = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Limit_Generated = read_int; #endif fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Limit_Invalid_Generated_States = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%lf", &read_double); USER_OPTIONS->Accepted_To_Generated_Ratio = read_double; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%lf", &read_double); USER_OPTIONS->Cost_Precision = read_double; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Maximum_Cost_Repeat = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Number_Cost_Samples = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%lf", &read_double); USER_OPTIONS->Temperature_Ratio_Scale = read_double; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%lf", &read_double); USER_OPTIONS->Cost_Parameter_Scale_Ratio = read_double; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%lf", &read_double); USER_OPTIONS->Temperature_Anneal_Scale = read_double; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Include_Integer_Parameters = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->User_Initial_Parameters = read_int; #if INT_ALLOC fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Sequential_Parameters = read_int; #else #if INT_LONG fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%ld", &read_long); USER_OPTIONS->Sequential_Parameters = read_long; #else fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Sequential_Parameters = read_int; #endif #endif fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%lf", &read_double); USER_OPTIONS->Initial_Parameter_Temperature = read_double; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Acceptance_Frequency_Modulus = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Generated_Frequency_Modulus = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Reanneal_Cost = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Reanneal_Parameters = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%lf", &read_double); USER_OPTIONS->Delta_X = read_double; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->User_Tangents = read_int; fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%d", &read_int); USER_OPTIONS->Curvature_0 = read_int; #else /* OPTIONS_FILE */ /* USER_OPTIONS->Limit_Acceptances = 10000; */ USER_OPTIONS->Limit_Acceptances = 1000; USER_OPTIONS->Limit_Generated = 99999; USER_OPTIONS->Limit_Invalid_Generated_States = 1000; USER_OPTIONS->Accepted_To_Generated_Ratio = 1.0E-6; USER_OPTIONS->Cost_Precision = 1.0E-18; USER_OPTIONS->Maximum_Cost_Repeat = 2; USER_OPTIONS->Number_Cost_Samples = 2; /* These variables are set below in x[.] */ /* USER_OPTIONS->Temperature_Ratio_Scale = 1.0E-5; */ /* USER_OPTIONS->Cost_Parameter_Scale_Ratio = 1.0; */ USER_OPTIONS->Temperature_Anneal_Scale = 100.; USER_OPTIONS->Include_Integer_Parameters = FALSE; USER_OPTIONS->User_Initial_Parameters = FALSE; USER_OPTIONS->Sequential_Parameters = -1; USER_OPTIONS->Initial_Parameter_Temperature = 1.0; USER_OPTIONS->Acceptance_Frequency_Modulus = 100; USER_OPTIONS->Generated_Frequency_Modulus = 10000; USER_OPTIONS->Reanneal_Cost = 1; USER_OPTIONS->Reanneal_Parameters = TRUE; USER_OPTIONS->Delta_X = 0.001; USER_OPTIONS->User_Tangents = FALSE; USER_OPTIONS->Curvature_0 = TRUE; #endif /* OPTIONS_FILE */ USER_OPTIONS->Temperature_Ratio_Scale = x[0]; USER_OPTIONS->Cost_Parameter_Scale_Ratio = x[1]; if (initial_flag == 0) { /* first value of *rand_seed */ #if ASA_LIB *rand_seed = (asa_rand_seed ? *asa_rand_seed : (LONG_INT) 696969); #else *rand_seed = 696969; #endif } if ((parameter_dimension = (ALLOC_INT *) calloc (1, sizeof (ALLOC_INT))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): parameter_dimension"); Exit_USER (user_exit_msg); return (-2); } if ((exit_code = (int *) calloc (1, sizeof (int))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): exit_code"); Exit_USER (user_exit_msg); return (-2); } if ((cost_flag = (int *) calloc (1, sizeof (int))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): cost_flag"); Exit_USER (user_exit_msg); return (-2); } /* the number of parameters for the cost function */ #if OPTIONS_FILE_DATA fscanf (ptr_options, "%s", read_option); fscanf (ptr_options, "%s", read_option); #if INT_ALLOC fscanf (ptr_options, "%d", &read_int); *parameter_dimension = read_int; #else #if INT_LONG fscanf (ptr_options, "%ld", &read_long); *parameter_dimension = read_long; #else fscanf (ptr_options, "%d", &read_int); *parameter_dimension = read_int; #endif #endif #else /* OPTIONS_FILE_DATA */ #if ASA_TEST /* set parameter dimension if SELF_OPTIMIZE=TRUE */ *parameter_dimension = 4; #endif /* ASA_TEST */ #endif /* OPTIONS_FILE_DATA */ #if MY_TEMPLATE /* MY_TEMPLATE_recur_dim */ /* If not using OPTIONS_FILE_DATA or data read from asa_opt, set parameter dimension if SELF_OPTIMIZE=TRUE */ #endif /* MY_TEMPLATE recur_dim */ if ((parameter_lower_bound = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): parameter_lower_bound"); Exit_USER (user_exit_msg); return (-2); } if ((parameter_upper_bound = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): parameter_upper_bound"); Exit_USER (user_exit_msg); return (-2); } if ((cost_parameters = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): cost_parameters"); Exit_USER (user_exit_msg); return (-2); } if ((parameter_int_real = (int *) calloc (*parameter_dimension, sizeof (int))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): parameter_int_real"); Exit_USER (user_exit_msg); return (-2); } if ((cost_tangents = (double *) calloc (*parameter_dimension, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): cost_tangents"); Exit_USER (user_exit_msg); return (-2); } if (USER_OPTIONS->Curvature_0 == FALSE || USER_OPTIONS->Curvature_0 == -1) { if ((cost_curvature = (double *) calloc ((*parameter_dimension) * (*parameter_dimension), sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): cost_curvature"); Exit_USER (user_exit_msg); return (-2); } } else { cost_curvature = (double *) NULL; } #if ASA_TEMPLATE_SELFOPT /* Set memory to that required for use. */ USER_OPTIONS->Asa_Data_Dim_Dbl = 2; if ((USER_OPTIONS->Asa_Data_Dbl = (double *) calloc (USER_OPTIONS->Asa_Data_Dim_Dbl, sizeof (double))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): USER_OPTIONS->Asa_Data_Dbl"); Exit_USER (user_exit_msg); return (-2); } /* Use Asa_Data_Dbl[0] as flag, e.g., if used with SELF_OPTIMIZE. */ USER_OPTIONS->Asa_Data_Dbl[0] = 1.0; #endif /* ASA_TEMPLATE_SELFOPT */ #if USER_COST_SCHEDULE USER_OPTIONS->Cost_Schedule = user_cost_schedule; #endif #if USER_ACCEPTANCE_TEST USER_OPTIONS->Acceptance_Test = user_acceptance_test; #endif #if USER_ACCEPT_ASYMP_EXP USER_OPTIONS->Asymp_Exp_Param = 1.0; #endif #if USER_GENERATING_FUNCTION USER_OPTIONS->Generating_Distrib = user_generating_distrib; #endif #if USER_REANNEAL_COST USER_OPTIONS->Reanneal_Cost_Function = user_reanneal_cost; #endif #if USER_REANNEAL_PARAMETERS USER_OPTIONS->Reanneal_Params_Function = user_reanneal_params; #endif initialize_parameters (cost_parameters, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, #if OPTIONS_FILE_DATA ptr_options, #endif USER_OPTIONS); #if OPTIONS_FILE fclose (ptr_options); #endif #if ASA_SAVE USER_OPTIONS->Random_Array_Dim = SHUFFLE; USER_OPTIONS->Random_Array = random_array; #endif /* ASA_SAVE */ /* It might be a good idea to place a loop around this call, and to average over several values of funevals returned by trajectories of cost_value. */ funevals = 0; #if USER_ASA_OUT if ((USER_OPTIONS->Asa_Out_File = (char *) calloc (80, sizeof (char))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): USER_OPTIONS->Asa_Out_File"); Exit_USER (user_exit_msg); return (-2); } #if ASA_TEMPLATE_SELFOPT strcpy (USER_OPTIONS->Asa_Out_File, "asa_rcur"); #endif #endif #if OPTIONAL_DATA_PTR data_ptr_flg = 1; #if ASA_TEMPLATE /* N.b.: If OPTIONAL_DATA_PTR is being used for RECUR_USER_OPTIONS * as well as for USER_OPTIONS, do not create (or free) additional memory * in recur_cost_function() for Asa_Data_Dim_Ptr and Asa_Data_Ptr to * be passed to the inner cost_function(), but rather link pointers to * those in RECUR_USER_OPTIONS. Typically, define separate structures * within the structure defined by Asa_Data_Ptr to access info depending * on whether the run in a particular level of cost function in this * recursive operation. In this case, set * #if TRUE to #if FALSE just * below. See the ASA-README for more discussion. */ #if TRUE USER_OPTIONS->Asa_Data_Dim_Ptr = 1; if ((USER_OPTIONS->Asa_Data_Ptr = (OPTIONAL_PTR_TYPE *) calloc (USER_OPTIONS->Asa_Data_Dim_Ptr, sizeof (OPTIONAL_PTR_TYPE))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): USER_OPTIONS->Asa_Data_Ptr"); Exit_USER (user_exit_msg); return (-2); } #else USER_OPTIONS->Asa_Data_Dim_Ptr = RECUR_USER_OPTIONS->Asa_Data_Dim_Ptr; USER_OPTIONS->Asa_Data_Ptr = RECUR_USER_OPTIONS->Asa_Data_Ptr; data_ptr_flg = 0; #endif #endif /* ASA_TEMPLATE */ USER_OPTIONS->Asa_Data_Dim_Ptr = 1; if ((USER_OPTIONS->Asa_Data_Ptr = (OPTIONAL_PTR_TYPE *) calloc (USER_OPTIONS->Asa_Data_Dim_Ptr, sizeof (OPTIONAL_PTR_TYPE))) == NULL) { strcpy (user_exit_msg, "recur_cost_function(): USER_OPTIONS->Asa_Data_Ptr"); Exit_USER (user_exit_msg); return (-2); } #endif /* OPTIONAL_DATA_PTR */ cost_value = asa (USER_COST_FUNCTION, randflt, rand_seed, cost_parameters, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, USER_OPTIONS); if (*exit_code == -1) { #if INCL_STDOUT printf ("\n\n*** error in calloc in ASA ***\n\n"); #endif /* INCL_STDOUT */ fprintf (ptr_out, "\n\n*** error in calloc in ASA ***\n\n"); fflush (ptr_out); return (-1); } #if MY_TEMPLATE /* MY_TEMPLATE_recur_post_asa */ #endif if (cost_value > .001) { *recur_cost_flag = FALSE; } else { *recur_cost_flag = TRUE; } #if FALSE /* set to 1 to activate FAST EXIT */ /* Make a quick exit */ if (recur_funevals >= 10) { *recur_cost_flag = FALSE; RECUR_USER_OPTIONS->Limit_Invalid_Generated_States = 0; fprintf (ptr_out, "FAST EXIT set at recur_funevals = 10\n\n"); } #endif #if TIME_CALC /* print every RECUR_PRINT_FREQUENCY evaluations */ if ((RECUR_PRINT_FREQUENCY > 0) && ((recur_funevals % RECUR_PRINT_FREQUENCY) == 0)) { USER_OPTIONS->Temperature_Ratio_Scale = x[0]; fprintf (ptr_out, "USER_OPTIONS->Temperature_Ratio_Scale = %12.7g\n", USER_OPTIONS->Temperature_Ratio_Scale); USER_OPTIONS->Cost_Parameter_Scale_Ratio = x[1]; fprintf (ptr_out, "USER_OPTIONS->Cost_Parameter_Scale_Ratio = %12.7g\n", USER_OPTIONS->Cost_Parameter_Scale_Ratio); } print_time ("", ptr_out); #endif fprintf (ptr_out, "recur_funevals = %ld, *recur_cost_flag = %d\n", recur_funevals, *recur_cost_flag); /* cost function = number generated at best cost */ #if ASA_TEMPLATE_SELFOPT funevals = (LONG_INT) (USER_OPTIONS->Asa_Data_Dbl[1]); fprintf (ptr_out, "\tbest_funevals = %ld, cost_value = %12.7g\n\n", funevals, cost_value); /* cost function = total number generated during run */ #endif /* ASA_TEMPLATE_SELFOPT */ #if ASA_SAMPLE fprintf (ptr_out, "\tfunevals = %ld, cost_value = %12.7g\n\n", funevals, cost_value); #endif fflush (ptr_out); #if ASA_TEMPLATE_SAMPLE ptr_asa = fopen ("asa_out", "r"); sample (ptr_out, ptr_asa); #endif #if OPTIONAL_DATA_DBL free (USER_OPTIONS->Asa_Data_Dbl); #endif #if OPTIONAL_DATA_INT free (USER_OPTIONS->Asa_Data_Int); #endif #if OPTIONAL_DATA_PTR if (data_ptr_flg == 1) { free (USER_OPTIONS->Asa_Data_Ptr); } #endif #if USER_ASA_OUT free (USER_OPTIONS->Asa_Out_File); #endif #if ASA_SAMPLE free (USER_OPTIONS->Bias_Generated); #endif #if ASA_QUEUE #if ASA_RESOLUTION #else free (USER_OPTIONS->Queue_Resolution); #endif #endif #if ASA_RESOLUTION free (USER_OPTIONS->Coarse_Resolution); #endif if (USER_OPTIONS->Curvature_0 == FALSE || USER_OPTIONS->Curvature_0 == -1) free (cost_curvature); #if USER_INITIAL_PARAMETERS_TEMPS free (USER_OPTIONS->User_Parameter_Temperature); #endif #if USER_INITIAL_COST_TEMP free (USER_OPTIONS->User_Cost_Temperature); #endif #if DELTA_PARAMETERS free (USER_OPTIONS->User_Delta_Parameter); #endif #if QUENCH_PARAMETERS free (USER_OPTIONS->User_Quench_Param_Scale); #endif #if QUENCH_COST free (USER_OPTIONS->User_Quench_Cost_Scale); #endif #if RATIO_TEMPERATURE_SCALES free (USER_OPTIONS->User_Temperature_Ratio); #endif #if MULTI_MIN free (USER_OPTIONS->Multi_Grid); for (multi_index = 0; multi_index < USER_OPTIONS->Multi_Number; ++multi_index) { free (USER_OPTIONS->Multi_Params[multi_index]); } #endif /* MULTI_MIN */ #if OPTIONAL_DATA_PTR if (data_ptr_flg == 0) { USER_OPTIONS = NULL; } #endif free (USER_OPTIONS); free (parameter_dimension); free (exit_code); free (cost_flag); free (parameter_lower_bound); free (parameter_upper_bound); free (cost_parameters); free (parameter_int_real); free (cost_tangents); free (rand_seed); return ((double) funevals); } #if USER_COST_SCHEDULE #if HAVE_ANSI double recur_user_cost_schedule (double test_temperature, USER_DEFINES * RECUR_USER_OPTIONS) #else double recur_user_cost_schedule (test_temperature, RECUR_USER_OPTIONS) double test_temperature; USER_DEFINES *RECUR_USER_OPTIONS; #endif /* HAVE_ANSI */ { #if ASA_TEMPLATE double x; x = test_temperature; return (x); #endif } #endif /* USER_COST_SCHEDULE */ #if USER_ACCEPTANCE_TEST #if HAVE_ANSI void recur_user_acceptance_test (double current_cost, double *recur_parameter_lower_bound, double *recur_parameter_upper_bound, ALLOC_INT * recur_parameter_dimension, USER_DEFINES * RECUR_USER_OPTIONS) #else void recur_user_acceptance_test (current_cost, recur_parameter_lower_bound, recur_parameter_upper_bound, recur_parameter_dimension, RECUR_USER_OPTIONS) double current_cost; double *recur_parameter_lower_bound; double *recur_parameter_upper_bound; ALLOC_INT *recur_parameter_dimension; USER_DEFINES *RECUR_USER_OPTIONS; #endif /* HAVE_ANSI */ { double uniform_test, curr_cost_temp; #if USER_ACCEPT_ASYMP_EXP double x, q, delta_cost; #endif #if ASA_TEMPLATE /* ASA cost index */ /* Calculate the current ASA cost index. This could be useful to define a new schedule for the cost temperature, beyond simple changes that can be made using USER_COST_SCHEDULE. */ int index; double k_temperature, quench, y; double xrecur_parameter_dimension; #if QUENCH_COST quench = RECUR_USER_OPTIONS->User_Quench_Cost_Scale[0]; #else quench = 1.0; #endif /* QUENCH_COST */ xrecur_parameter_dimension = (double) *recur_parameter_dimension; for (index = 0; index < *recur_parameter_dimension; ++index) if (fabs (recur_parameter_upper_bound[index] - recur_parameter_lower_bound[index]) < (double) EPS_DOUBLE) *xrecur_parameter_dimension -= 1.0; y = -F_LOG (RECUR_USER_OPTIONS->Cost_Temp_Curr / RECUR_USER_OPTIONS->Cost_Temp_Init) / RECUR_USER_OPTIONS->Cost_Temp_Scale; k_temperature = F_POW (y, xrecur_parameter_dimension / quench); #endif /* ASA cost index */ uniform_test = randflt (RECUR_USER_OPTIONS->Random_Seed); curr_cost_temp = RECUR_USER_OPTIONS->Cost_Temp_Curr; #if ASA_TEMPLATE #if USER_COST_SCHEDULE curr_cost_temp = (RECUR_USER_OPTIONS->Cost_Schedule (RECUR_USER_OPTIONS->Cost_Temp_Curr, RECUR_USER_OPTIONS) + (double) EPS_DOUBLE); #else curr_cost_temp = RECUR_USER_OPTIONS->Cost_Temp_Curr; #endif #endif /* ASA_TEMPLATE */ #if USER_ACCEPT_ASYMP_EXP #if USER_COST_SCHEDULE curr_cost_temp = (RECUR_USER_OPTIONS->Cost_Schedule (RECUR_USER_OPTIONS->Cost_Temp_Curr, RECUR_USER_OPTIONS) + (double) EPS_DOUBLE); #endif delta_cost = (current_cost - *(RECUR_USER_OPTIONS->Last_Cost)) / (curr_cost_temp + (double) EPS_DOUBLE); q = RECUR_USER_OPTIONS->Asymp_Exp_Param; if (fabs (1.0 - q) < (double) EPS_DOUBLE) x = MIN (1.0, (F_EXP (-delta_cost))); /* Boltzmann test */ else if ((1.0 - (1.0 - q) * delta_cost) < (double) EPS_DOUBLE) x = MIN (1.0, (F_EXP (-delta_cost))); /* Boltzmann test */ else x = MIN (1.0, F_POW ((1.0 - (1.0 - q) * delta_cost), (1.0 / (1.0 - q)))); RECUR_USER_OPTIONS->Prob_Bias = x; if (x >= uniform_test) RECUR_USER_OPTIONS->User_Acceptance_Flag = TRUE; else RECUR_USER_OPTIONS->User_Acceptance_Flag = FALSE; #endif /* USER_ACCEPT_ASYMP_EXP */ } #endif /* USER_ACCEPTANCE_TEST */ #if USER_GENERATING_FUNCTION #if HAVE_ANSI double recur_user_generating_distrib (LONG_INT * seed, ALLOC_INT * recur_parameter_dimension, ALLOC_INT index_v, double temperature_v, double init_param_temp_v, double temp_scale_params_v, double parameter_v, double parameter_range_v, double *last_saved_parameter, USER_DEFINES * RECUR_USER_OPTIONS) #else double recur_user_generating_distrib (seed, recur_parameter_dimension, index_v, temperature_v, init_param_temp_v, temp_scale_params_v, parameter_v, parameter_range_v, last_saved_parameter, RECUR_USER_OPTIONS) LONG_INT *seed; ALLOC_INT *recur_parameter_dimension; ALLOC_INT index_v; double temperature_v; double init_param_temp_v; double temp_scale_params_v; double parameter_v; double parameter_range_v; double *last_saved_parameter; USER_DEFINES *RECUR_USER_OPTIONS; #endif { #if ASA_TEMPLATE double x, y, z; /* This is the ASA distribution. A slower temperature schedule can be obtained here, e.g., temperature_v = pow(temperature_v, 0.5); */ x = randflt (seed); y = x < 0.5 ? -1.0 : 1.0; z = y * temperature_v * (F_POW ((1.0 + 1.0 / temperature_v), fabs (2.0 * x - 1.0)) - 1.0); x = parameter_v + z * parameter_range_v; return (x); #endif /* ASA_TEMPLATE */ } #endif /* USER_GENERATING_FUNCTION */ #if USER_REANNEAL_COST #if HAVE_ANSI int recur_user_reanneal_cost (double *cost_best, double *cost_last, double *initial_cost_temperature, double *current_cost_temperature, USER_DEFINES * RECUR_USER_OPTIONS) #else int recur_user_reanneal_cost (cost_best, cost_last, initial_cost_temperature, current_cost_temperature, RECUR_USER_OPTIONS) double *cost_best; double *cost_last; double *initial_cost_temperature; double *current_cost_temperature; USER_DEFINES *RECUR_USER_OPTIONS; #endif /* HAVE_ANSI */ { #if ASA_TEMPLATE double tmp_dbl; tmp_dbl = MAX (fabs (*cost_last), fabs (*cost_best)); tmp_dbl = MAX ((double) EPS_DOUBLE, tmp_dbl); *initial_cost_temperature = MIN (*initial_cost_temperature, tmp_dbl); return (TRUE); #endif } #endif /* USER_REANNEAL_COST */ #if USER_REANNEAL_PARAMETERS #if HAVE_ANSI double recur_user_reanneal_params (double current_temp, double tangent, double max_tangent, USER_DEFINES * RECUR_USER_OPTIONS) #else double recur_user_reanneal_params (current_temp, tangent, max_tangent, RECUR_USER_OPTIONS) double current_temp; double tangent; double max_tangent; USER_DEFINES *RECUR_USER_OPTIONS; #endif /* HAVE_ANSI */ { #if ASA_TEMPLATE double x; x = current_temp * (max_tangent / tangent); return (x); #endif } #endif /* USER_REANNEAL_PARAMETERS */ #endif /* SELF_OPTIMIZE */ #if FITLOC #if HAVE_ANSI double calcf (double (*user_cost_function) (double *, double *, double *, double *, double *, ALLOC_INT *, int *, int *, int *, USER_DEFINES *), double *xloc, double *parameter_lower_bound, double *parameter_upper_bound, double *cost_tangents, double *cost_curvature, ALLOC_INT * parameter_dimension, int *parameter_int_real, int *cost_flag, int *exit_code, USER_DEFINES * OPTIONS, FILE * ptr_out) #else double calcf (user_cost_function, xloc, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out) double (*user_cost_function) (); double *x; double *parameter_lower_bound; double *parameter_upper_bound; double *cost_tangents; double *cost_curvature; ALLOC_INT *parameter_dimension; int *parameter_int_real; int *cost_flag; int *exit_code; USER_DEFINES *OPTIONS; FILE *ptr_out; #endif { ALLOC_INT index_v; #if FITLOC_ROUND double x, min_parameter_v, max_parameter_v, parameter_range_v; #endif double floc; #if ASA_RESOLUTION double xres, xint, xplus, xminus, dx, dxminus, dxplus; #endif #if FITLOC_ROUND /* The following section for adjustments of parameters is taken from generate_new_state() in asa.c */ for (index_v = 0; index_v < *parameter_dimension; ++index_v) { if (fabs (parameter_lower_bound[index_v] - parameter_upper_bound[index_v]) < EPS_DOUBLE) continue; x = xloc[index_v]; min_parameter_v = parameter_lower_bound[index_v]; max_parameter_v = parameter_upper_bound[index_v]; parameter_range_v = max_parameter_v - min_parameter_v; /* Handle discrete parameters. */ #if ASA_RESOLUTION xres = OPTIONS->Coarse_Resolution[index_v]; if (xres > EPS_DOUBLE) { min_parameter_v -= (xres / 2.0); max_parameter_v += (xres / 2.0); parameter_range_v = max_parameter_v - min_parameter_v; } #endif /* ASA_RESOLUTION */ if (parameter_int_real[index_v] > 0) { #if ASA_RESOLUTION if (xres > EPS_DOUBLE) { ; } else { #endif /* ASA_RESOLUTION */ min_parameter_v -= 0.5; max_parameter_v += 0.5; parameter_range_v = max_parameter_v - min_parameter_v; } #if ASA_RESOLUTION } #endif #if ASA_RESOLUTION if (xres > EPS_DOUBLE) { xint = xres * (double) ((LONG_INT) (x / xres)); xplus = xint + xres; xminus = xint - xres; dx = fabs (xint - x); dxminus = fabs (xminus - x); dxplus = fabs (xplus - x); if (dx < dxminus && dx < dxplus) x = xint; else if (dxminus < dxplus) x = xminus; else x = xplus; } #endif /* ASA_RESOLUTION */ /* Handle discrete parameters. You might have to check rounding on your machine. */ if (parameter_int_real[index_v] > 0) { #if ASA_RESOLUTION if (xres > EPS_DOUBLE) { ; } else { #endif /* ASA_RESOLUTION */ if (x < min_parameter_v + 0.5) x = min_parameter_v + 0.5 + (double) EPS_DOUBLE; if (x > max_parameter_v - 0.5) x = max_parameter_v - 0.5 + (double) EPS_DOUBLE; if (x + 0.5 > 0.0) { x = (double) ((LONG_INT) (x + 0.5)); } else { x = (double) ((LONG_INT) (x - 0.5)); } if (x > parameter_upper_bound[index_v]) x = parameter_upper_bound[index_v]; if (x < parameter_lower_bound[index_v]) x = parameter_lower_bound[index_v]; } #if ASA_RESOLUTION } if (xres > EPS_DOUBLE) { if (x < min_parameter_v + xres / 2.0) x = min_parameter_v + xres / 2.0 + (double) EPS_DOUBLE; if (x > max_parameter_v - xres / 2.0) x = max_parameter_v - xres / 2.0 + (double) EPS_DOUBLE; if (x > parameter_upper_bound[index_v]) x = parameter_upper_bound[index_v]; if (x < parameter_lower_bound[index_v]) x = parameter_lower_bound[index_v]; } #endif /* ASA_RESOLUTION */ if ((x < parameter_lower_bound[index_v]) || (x > parameter_upper_bound[index_v])) { ; } else { xloc[index_v] = x; } } #endif /* FITLOC_ROUND */ floc = user_cost_function (xloc, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS); if (*cost_flag == FALSE) { floc += OPTIONS->Penalty; } for (index_v = 0; index_v < *parameter_dimension; ++index_v) { if (parameter_upper_bound[index_v] - xloc[index_v] < EPS_DOUBLE) floc += OPTIONS->Penalty; else if (xloc[index_v] - parameter_lower_bound[index_v] < EPS_DOUBLE) floc += OPTIONS->Penalty; } return (floc); } #if HAVE_ANSI double fitloc (double (*user_cost_function) (double *, double *, double *, double *, double *, ALLOC_INT *, int *, int *, int *, USER_DEFINES *), double *xloc, double *parameter_lower_bound, double *parameter_upper_bound, double *cost_tangents, double *cost_curvature, ALLOC_INT * parameter_dimension, int *parameter_int_real, int *cost_flag, int *exit_code, USER_DEFINES * OPTIONS, FILE * ptr_out) #else double fitloc (user_cost_function, xloc, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out) double (*user_cost_function) (); double *xloc; double *parameter_lower_bound; double *parameter_upper_bound; double *cost_tangents; double *cost_curvature; ALLOC_INT *parameter_dimension; int *parameter_int_real; int *cost_flag; int *exit_code; USER_DEFINES *OPTIONS; FILE *ptr_out; #endif { double x; ALLOC_INT index_v; #if FITLOC_ROUND double min_parameter_v, max_parameter_v, parameter_range_v; #endif double *xsave; double tol1, tol2, alpha, beta1, beta2, gamma, delta, floc, fsave, ffinal; int no_progress, tot_iters, locflg, bndflg; #if ASA_RESOLUTION double xres, xint, xminus, xplus, dx, dxminus, dxplus; #endif #if FITLOC_PRINT if (OPTIONS->Fit_Local >= 1) { fprintf (ptr_out, "\n\nSTART LOCAL FIT\n"); } else { fprintf (ptr_out, "\n\nSTART LOCAL FIT Independent of ASA\n"); } fflush (ptr_out); #endif /* FITLOC_PRINT */ xsave = (double *) calloc (*parameter_dimension, sizeof (double)); bndflg = 0; /* The following simplex parameters may need adjustments for your system. */ tol1 = EPS_DOUBLE; tol2 = EPS_DOUBLE * 100.; no_progress = 4; alpha = 1.0; beta1 = 0.75; beta2 = 0.75; gamma = 1.25; delta = 2.50; for (index_v = 0; index_v < *parameter_dimension; ++index_v) { xsave[index_v] = xloc[index_v]; } fsave = user_cost_function (xloc, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS); tot_iters = simplex (user_cost_function, xloc, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out, tol1, tol2, no_progress, alpha, beta1, beta2, gamma, delta); fflush (ptr_out); for (index_v = 0; index_v < *parameter_dimension; ++index_v) { x = xloc[index_v]; if ((x < parameter_lower_bound[index_v]) || (x > parameter_upper_bound[index_v])) { bndflg = 1; } } /* The following section for adjustments of parameters is taken from generate_new_state() in asa.c */ #if FITLOC_ROUND for (index_v = 0; index_v < *parameter_dimension; ++index_v) { if (fabs (parameter_lower_bound[index_v] - parameter_upper_bound[index_v]) < EPS_DOUBLE) continue; x = xloc[index_v]; min_parameter_v = parameter_lower_bound[index_v]; max_parameter_v = parameter_upper_bound[index_v]; parameter_range_v = max_parameter_v - min_parameter_v; /* Handle discrete parameters. */ #if ASA_RESOLUTION xres = OPTIONS->Coarse_Resolution[index_v]; if (xres > EPS_DOUBLE) { min_parameter_v -= (xres / 2.0); max_parameter_v += (xres / 2.0); parameter_range_v = max_parameter_v - min_parameter_v; } #endif /* ASA_RESOLUTION */ if (parameter_int_real[index_v] > 0) { #if ASA_RESOLUTION if (xres > EPS_DOUBLE) { ; } else { #endif /* ASA_RESOLUTION */ min_parameter_v -= 0.5; max_parameter_v += 0.5; parameter_range_v = max_parameter_v - min_parameter_v; } #if ASA_RESOLUTION } #endif #if ASA_RESOLUTION if (xres > EPS_DOUBLE) { xint = xres * (double) ((LONG_INT) (x / xres)); xplus = xint + xres; xminus = xint - xres; dx = fabs (xint - x); dxminus = fabs (xminus - x); dxplus = fabs (xplus - x); if (dx < dxminus && dx < dxplus) x = xint; else if (dxminus < dxplus) x = xminus; else x = xplus; } #endif /* ASA_RESOLUTION */ /* Handle discrete parameters. You might have to check rounding on your machine. */ if (parameter_int_real[index_v] > 0) { #if ASA_RESOLUTION if (xres > EPS_DOUBLE) { ; } else { #endif /* ASA_RESOLUTION */ if (x < min_parameter_v + 0.5) x = min_parameter_v + 0.5 + (double) EPS_DOUBLE; if (x > max_parameter_v - 0.5) x = max_parameter_v - 0.5 + (double) EPS_DOUBLE; if (x + 0.5 > 0.0) { x = (double) ((LONG_INT) (x + 0.5)); } else { x = (double) ((LONG_INT) (x - 0.5)); } if (x > parameter_upper_bound[index_v]) x = parameter_upper_bound[index_v]; if (x < parameter_lower_bound[index_v]) x = parameter_lower_bound[index_v]; } #if ASA_RESOLUTION } if (xres > EPS_DOUBLE) { if (x < min_parameter_v + xres / 2.0) x = min_parameter_v + xres / 2.0 + (double) EPS_DOUBLE; if (x > max_parameter_v - xres / 2.0) x = max_parameter_v - xres / 2.0 + (double) EPS_DOUBLE; if (x > parameter_upper_bound[index_v]) x = parameter_upper_bound[index_v]; if (x < parameter_lower_bound[index_v]) x = parameter_lower_bound[index_v]; } #endif /* ASA_RESOLUTION */ if ((x < parameter_lower_bound[index_v]) || (x > parameter_upper_bound[index_v])) { bndflg = 1; #if FITLOC_PRINT if (OPTIONS->Fit_Local == 2) fprintf (ptr_out, "IGNORE FITLOC: OUT OF BOUNDS xloc[%ld] = %g\n", index_v, xloc[index_v]); else fprintf (ptr_out, "OUT OF BOUNDS xloc[%ld] = %g\n", index_v, xloc[index_v]); #else ; #endif /* FITLOC_PRINT */ } else { xloc[index_v] = x; } } #endif /* FITLOC_ROUND */ floc = user_cost_function (xloc, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS); if (fabs (floc - fsave) < (double) EPS_DOUBLE) { locflg = 1; ffinal = fsave; #if FITLOC_PRINT fprintf (ptr_out, "\nsame global cost = %g\tlocal cost = %g\n\n", fsave, floc); #endif /* FITLOC_PRINT */ } else { if (floc < fsave) { if (OPTIONS->Fit_Local == 2 && bndflg == 1) { locflg = 1; ffinal = fsave; } else { locflg = 0; ffinal = floc; } } else { locflg = 1; ffinal = fsave; } #if FITLOC_PRINT fprintf (ptr_out, "\nDIFF global cost = %g\tlocal cost = %g\n\n", fsave, floc); #endif /* FITLOC_PRINT */ } for (index_v = 0; index_v < *parameter_dimension; ++index_v) { if (fabs (xloc[index_v] - xsave[index_v]) < (double) EPS_DOUBLE) { #if FITLOC_PRINT fprintf (ptr_out, "same global param[%ld] = %g\tlocal param = %g\n", index_v, xsave[index_v], xloc[index_v]); #else ; #endif /* FITLOC_PRINT */ } else { #if FITLOC_PRINT fprintf (ptr_out, "DIFF global param[%ld] = %g\tlocal param = %g\n", index_v, xsave[index_v], xloc[index_v]); #else ; #endif /* FITLOC_PRINT */ if (locflg == 1) { xloc[index_v] = xsave[index_v]; } } } #if FITLOC_PRINT fprintf (ptr_out, "\n"); fflush (ptr_out); #endif /* FITLOC_PRINT */ free (xsave); return (ffinal); } /* Written by Mark Johnson , based on %A J.A. Nelder %A R. Mead %T A simplex method for function minimization %J Computer J. (UK) %V 7 %D 1964 %P 308-313 with improvements from %A G.P. Barabino %A G.S. Barabino %A B. Bianco %A M. Marchesi %T A study on the performances of simplex methods for function minimization %B Proc. IEEE Int. Conf. Circuits and Computers %D 1980 %P 1150-1153 adapted for use in ASA by Lester Ingber */ #if HAVE_ANSI int simplex (double (*user_cost_function) (double *, double *, double *, double *, double *, ALLOC_INT *, int *, int *, int *, USER_DEFINES *), double *x, double *parameter_lower_bound, double *parameter_upper_bound, double *cost_tangents, double *cost_curvature, ALLOC_INT * parameter_dimension, int *parameter_int_real, int *cost_flag, int *exit_code, USER_DEFINES * OPTIONS, FILE * ptr_out, double tol1, double tol2, int no_progress, double alpha, double beta1, double beta2, double gamma, double delta) #else int simplex (user_cost_function, x, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out, tol1, tol2, no_progress, alpha, beta1, beta2, gamma, delta) double (*user_cost_function) (); double *x; double *parameter_lower_bound; double *parameter_upper_bound; double *cost_tangents; double *cost_curvature; ALLOC_INT *parameter_dimension; int *parameter_int_real; int *cost_flag; int *exit_code; USER_DEFINES *OPTIONS; FILE *ptr_out; double tol1; double tol2; int no_progress; double alpha; double beta1; double beta2; double gamma; double delta; #endif { double fs, fl, fh, fr, fe, fc1, fc2, ftmp, flast; double err1; double *fvals; double **splx; /* the simplex of points */ double *x0; /* centroid of simplex */ double *xr; /* point for a reflection */ double *xe; /* point for an expansion */ double *xc1; /* point for a minor contraction */ double *xc2; /* point for a major contraction */ int s, l, h; int i, j, iters, futility; int lastprint; fvals = (double *) calloc (*parameter_dimension + 1, sizeof (double)); splx = (double **) calloc (*parameter_dimension + 1, sizeof (double *)); for (i = 0; i <= *parameter_dimension; i++) splx[i] = (double *) calloc (*parameter_dimension, sizeof (double)); x0 = (double *) calloc (*parameter_dimension, sizeof (double)); xr = (double *) calloc (*parameter_dimension, sizeof (double)); xe = (double *) calloc (*parameter_dimension, sizeof (double)); xc1 = (double *) calloc (*parameter_dimension, sizeof (double)); xc2 = (double *) calloc (*parameter_dimension, sizeof (double)); /* build the initial simplex */ for (i = 0; i < *parameter_dimension; i++) { splx[0][i] = x[i]; } for (i = 1; i <= *parameter_dimension; i++) { for (j = 0; j < *parameter_dimension; j++) { if ((j + 1) == i) splx[i][j] = (x[j] * 2.25) + tol2; else splx[i][j] = x[j]; xr[j] = splx[i][j]; } fvals[i] = calcf (user_cost_function, xr, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out); } /* and of course compute function at starting point */ fvals[0] = calcf (user_cost_function, x, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out); /* now find the largest, 2nd largest, smallest f values */ if (fvals[0] > fvals[1]) { h = 0; s = 1; l = 1; } else { h = 1; s = 0; l = 0; } fh = fvals[h]; fs = fvals[s]; fl = fvals[l]; for (i = 2; i <= *parameter_dimension; i++) { if (fvals[i] <= fvals[l]) { l = i; fl = fvals[i]; } else { if (fvals[i] >= fvals[h]) { s = h; fs = fh; h = i; fh = fvals[i]; } else if (fvals[i] >= fvals[s]) { s = i; fs = fvals[i]; } } } #if FITLOC_PRINT if ((s == h) || (s == l) || (h == l)) fprintf (ptr_out, "\nPANIC: s,l,h not unique %d %d %d\n", s, h, l); fprintf (ptr_out, "INITIAL SIMPLEX:\n"); for (i = 0; i <= *parameter_dimension; i++) { for (j = 0; j < *parameter_dimension; j++) { fprintf (ptr_out, " %11.4g", splx[i][j]); } fprintf (ptr_out, " f = %12.5g", fvals[i]); if (i == h) fprintf (ptr_out, " HIGHEST"); if (i == s) fprintf (ptr_out, " SECOND HIGHEST"); if (i == l) fprintf (ptr_out, " LOWEST"); fprintf (ptr_out, "\n"); } #endif /* FITLOC_PRINT */ /* MAJOR LOOP */ flast = fl; futility = 0; lastprint = 0; iters = 0; err1 = 1.1 + (1.1 * tol1); while ((err1 > tol1) && (iters < OPTIONS->Iter_Max) && (futility < (*parameter_dimension * no_progress))) { iters++; /* now find the largest, 2nd largest, smallest f values */ if (fvals[0] > fvals[1]) { h = 0; s = 1; l = 1; } else { h = 1; s = 0; l = 0; } fh = fvals[h]; fs = fvals[s]; fl = fvals[l]; for (i = 2; i <= *parameter_dimension; i++) { if (fvals[i] <= fvals[l]) { l = i; fl = fvals[i]; } else { if (fvals[i] >= fvals[h]) { s = h; fs = fh; h = i; fh = fvals[i]; } else if (fvals[i] >= fvals[s]) { s = i; fs = fvals[i]; } } } #if FITLOC_PRINT if ((s == h) || (s == l) || (h == l)) fprintf (ptr_out, "\nPANIC: s,l,h not unique %d %d %d\n", s, h, l); #endif /* compute the centroid */ for (j = 0; j < *parameter_dimension; j++) { x0[j] = 0.0; for (i = 0; i <= *parameter_dimension; i++) { if (i != h) x0[j] += splx[i][j]; } x0[j] /= ((double) *parameter_dimension); } if (fl < flast) { flast = fl; futility = 0; } else futility += 1; #if FITLOC_PRINT fprintf (ptr_out, "Iteration %3d f(best) = %12.6g halt? = %11.5g\n", iters, fl, err1); if ((iters - lastprint) >= 100) { fprintf (ptr_out, "\n Best point seen so far:\n"); for (i = 0; i < *parameter_dimension; i++) { fprintf (ptr_out, " x[%3d] = %15.7g\n", i, splx[l][i]); } lastprint = iters; fprintf (ptr_out, "\n"); } fflush (ptr_out); #endif /* FITLOC_PRINT */ /* STEP 1: compute a reflected point xr */ for (i = 0; i < *parameter_dimension; i++) { xr[i] = ((1.0 + alpha) * x0[i]) - (alpha * splx[h][i]); } fr = calcf (user_cost_function, xr, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out); /* typical outcome: <2nd-biggest , >lowest . Go again */ if ((fr < fs) && (fr > fl)) { for (i = 0; i < *parameter_dimension; i++) { splx[h][i] = xr[i]; } fvals[h] = fr; goto more_iterations; } /* STEP 2: if reflected point is favorable, expand the simplex */ if (fr < fl) { for (i = 0; i < *parameter_dimension; i++) { xe[i] = (gamma * xr[i]) + ((1.0 - gamma) * x0[i]); } fe = calcf (user_cost_function, xe, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out); if (fe < fr) { /* win big; expansion point tiny */ for (i = 0; i < *parameter_dimension; i++) { splx[h][i] = xe[i]; } fvals[h] = fh = fe; } else /* still ok; reflection point a winner */ { for (i = 0; i < *parameter_dimension; i++) { splx[h][i] = xr[i]; } fvals[h] = fh = fr; } goto more_iterations; } /* STEP 3: if reflected point is unfavorable, contract simplex */ if (fr > fs) { if (fr < fh) { /* may as well replace highest pt */ for (i = 0; i < *parameter_dimension; i++) { splx[h][i] = xr[i]; } fvals[h] = fh = fr; } for (i = 0; i < *parameter_dimension; i++) { xc1[i] = (beta1 * xr[i]) + ((1.0 - beta1) * x0[i]); } fc1 = calcf (user_cost_function, xc1, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out); if (fc1 < fh) { /* slight contraction worked */ for (i = 0; i < *parameter_dimension; i++) { splx[h][i] = xc1[i]; } fvals[h] = fh = fc1; goto more_iterations; } /* now have to try strong contraction */ for (i = 0; i < *parameter_dimension; i++) { xc2[i] = (beta2 * splx[h][i]) + ((1.0 - beta2) * x0[i]); } fc2 = calcf (user_cost_function, xc2, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out); if (fc2 < fh) { /* strong contraction worked */ for (i = 0; i < *parameter_dimension; i++) { splx[h][i] = xc2[i]; } fvals[h] = fh = fc2; goto more_iterations; } } /* STEP 4: nothing worked. collapse the simplex around xl */ for (i = 0; i <= *parameter_dimension; i++) { if (i != l) { for (j = 0; j < *parameter_dimension; j++) { splx[i][j] = (splx[i][j] + splx[l][j]) / delta; xr[j] = splx[i][j]; } fvals[i] = calcf (user_cost_function, xr, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, OPTIONS, ptr_out); } } more_iterations: ftmp = 0.00; for (i = 0; i <= *parameter_dimension; i++) { ftmp += fvals[i]; } ftmp /= ((double) (*parameter_dimension + 1)); err1 = 0.00; for (i = 0; i <= *parameter_dimension; i++) { err1 += ((fvals[i] - ftmp) * (fvals[i] - ftmp)); } err1 /= ((double) (*parameter_dimension + 1)); err1 = sqrt (err1); } /* end of major while loop */ /* find the smallest f value */ l = 0; fl = fvals[0]; for (i = 1; i <= *parameter_dimension; i++) { if (fvals[i] < fvals[l]) l = i; } /* give it back to the user */ for (i = 0; i < *parameter_dimension; i++) { x[i] = splx[l][i]; } free (fvals); for (i = 0; i <= *parameter_dimension; i++) free (splx[i]); free (splx); free (x0); free (xr); free (xe); free (xc1); free (xc2); return (iters); } #else #endif /* FITLOC */ #if ASA_TEMPLATE_SAMPLE #if HAVE_ANSI void sample (FILE * ptr_out, FILE * ptr_asa) #else void sample (ptr_out, ptr_asa) FILE *ptr_out; FILE *ptr_asa; #endif { int ind, n_samples, n_accept, index, dim; double cost, cost_temp, bias_accept; double param, temp, bias_gener, aver_weight, range; double sum, norm, answer, prod, binsize; char ch[80], sample[8]; /* This is a demonstration of using ASA_SAMPLE to perform the double integral of exp(-x^2 - y^2) for x and y between 0 and 2. The mesh is quite crude. The temperature-dependent acceptance and generated biases factor are divided out, and the actual cost function weights each point. */ dim = 2; norm = sum = 0.; n_samples = 0; fprintf (ptr_out, ":SAMPLE: n_accept cost cost_temp bias_accept \ aver_weight\n"); fprintf (ptr_out, ":SAMPLE: index param[] temp[] bias_gener[] \ range[]\n"); for (;;) { fscanf (ptr_asa, "%s", ch); if (!strcmp (ch, "exit_status")) { break; } if (strcmp (ch, ":SAMPLE#")) { continue; } ++n_samples; fprintf (ptr_out, "%s\n", ch); fflush (ptr_out); fscanf (ptr_asa, "%s%d%lf%lf%lf%lf", sample, &n_accept, &cost, &cost_temp, &bias_accept, &aver_weight); if (strcmp (sample, ":SAMPLE+")) { fprintf (ptr_out, "%s %11d %12.7g %12.7g %12.7g %12.7g\n", sample, n_accept, cost, cost_temp, bias_accept, aver_weight); } else { fprintf (ptr_out, "%s %10d %12.7g %12.7g %12.7g %12.7g\n", sample, n_accept, cost, cost_temp, bias_accept, aver_weight); } prod = bias_accept; binsize = 1.0; for (ind = 0; ind < dim; ++ind) { fscanf (ptr_asa, "%s%d%lf%lf%lf%lf", sample, &index, ¶m, &temp, &bias_gener, &range); fprintf (ptr_out, "%s %11d %12.7g %12.7g %12.7g %12.7g\n", sample, index, param, temp, bias_gener, range); prod *= bias_gener; binsize *= range; } /* In this example, retrieve integrand from sampling function */ sum += ((F_EXP (-cost) * binsize) / prod); norm += (binsize / prod); } sum /= norm; answer = 1.0; for (ind = 0; ind < dim; ++ind) { answer *= (0.5 * sqrt (3.14159265) * erf (2.0)); } fprintf (ptr_out, "\n"); fprintf (ptr_out, "sum = %12.7g, answer = %12.7g\n", sum, answer); fprintf (ptr_out, "n_samples = %d, norm = %12.7g\n", n_samples, norm); fflush (ptr_out); } #endif /* ASA_TEMPLATE_SAMPLE */ #if ASA_TEMPLATE_LIB int main () { double main_cost_value; double *main_cost_parameters; int main_exit_code; LONG_INT number_params; ALLOC_INT n_param; FILE *ptr_main; #if INCL_STDOUT ptr_main = stdout; #endif /* INCL_STDOUT */ /* Note this assumes the *parameter_dimension = 4 */ number_params = 4; if ((main_cost_parameters = (double *) calloc (number_params, sizeof (double))) == NULL) { strcpy (user_exit_msg, "ASA_TEMPLATE_LIB main(): main_cost_parameters"); Exit_USER (user_exit_msg); return (-2); } asa_seed (696969); /* This is the default random seed. */ asa_main (&main_cost_value, main_cost_parameters, &main_exit_code); fprintf (ptr_main, "main_exit_code = %d\n", main_exit_code); fprintf (ptr_main, "main_cost_value = %12.7g\n", main_cost_value); fprintf (ptr_main, "parameter\tvalue\n"); for (n_param = 0; n_param < number_params; ++n_param) { fprintf (ptr_main, #if INT_ALLOC "%d\t\t%12.7g\n", #else #if INT_LONG "%ld\t\t%12.7g\n", #else "%d\t\t%12.7g\n", #endif #endif n_param, main_cost_parameters[n_param]); } free (main_cost_parameters); return (0); /* NOTREACHED */ } #endif /* ASA_TEMPLATE_LIB */ void Exit_USER (char *statement) { #if INCL_STDOUT printf ("\n\n*** EXIT calloc failed *** %s\n\n", statement); #else ; #endif /* INCL_STDOUT */ }