asa_usr.c

simulated annealing code ASA

   * %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;
  } 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 */
  /* 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)

#if ASA_LIB
/***********************************************************************
* LONG_INT asa_seed - returns initial random seed
***********************************************************************/

#if HAVE_ANSI
LONG_INT
asa_seed (LONG_INT seed)
#else
LONG_INT
asa_seed (seed)
     LONG_INT seed;
#endif
{
  static LONG_INT rand_seed;

  if (fabs (seed) > 0) {
    asa_rand_seed = &rand_seed;
    rand_seed = seed;
  }

  return (rand_seed);
}
#endif /* ASA_LIB */

/***********************************************************************
* double myrand - returns random number between 0 and 1
*	This routine returns the random number generator between 0 and 1
***********************************************************************/

#if HAVE_ANSI
double
myrand (LONG_INT * rand_seed)
#else
double
myrand (rand_seed)
     LONG_INT *rand_seed;
#endif
  /* 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
***********************************************************************/

#if HAVE_ANSI
double
randflt (LONG_INT * rand_seed)
#else
double
randflt (rand_seed)
     LONG_INT *rand_seed;
#endif
{
  return (resettable_randflt (rand_seed, 0));
}

/***********************************************************************
* double resettable_randflt
***********************************************************************/

#if HAVE_ANSI
double
resettable_randflt (LONG_INT * rand_seed, int reset)
#else
double
resettable_randflt (rand_seed, reset)
     LONG_INT *rand_seed;
     int reset;
#endif
  /* 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;
#if ASA_SAVE
  /* random_array[] local to all of asa_usr.c set at top of file */
#else
  static double random_array[SHUFFLE];  /* random variables */
#endif

  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 */