diff -urp gsl-1.9-org/siman/siman.c gsl-1.9/siman/siman.c --- gsl-1.9-org/siman/siman.c 2005-11-15 18:22:47.000000000 +0200 +++ gsl-1.9/siman/siman.c 2007-05-25 22:25:02.000000000 +0300 @@ -33,6 +33,22 @@ double safe_exp (double x) /* avoid unde { return (x < GSL_LOG_DBL_MIN) ? 0.0 : exp(x); } + +static inline double +boltzmann(double E, double new_E, double T, gsl_siman_params_t *params) +{ + return safe_exp (-(new_E - E) / (params->k * T)); +} + +static inline void +copy_state(void *src, void *dst, size_t size, gsl_siman_copy_t copyfunc) +{ + if (copyfunc) { + copyfunc(src, dst); + } else { + memcpy(dst, src, size); + } +} /* implementation of a basic simulated annealing algorithm */ @@ -49,8 +65,8 @@ gsl_siman_solve (const gsl_rng * r, void { void *x, *new_x, *best_x; double E, new_E, best_E; - int i, done; - double T; + int i; + double T, T_factor; int n_evals = 1, n_iter = 0, n_accepts, n_rejects, n_eless; /* this function requires that either the dynamic functions (copy, @@ -67,67 +83,57 @@ gsl_siman_solve (const gsl_rng * r, void new_x = copy_constructor(x0_p); best_x = copy_constructor(x0_p); } else { - x = (void *) malloc (element_size); + x = malloc (element_size); memcpy (x, x0_p, element_size); - new_x = (void *) malloc (element_size); - best_x = (void *) malloc (element_size); + new_x = malloc (element_size); + best_x = malloc (element_size); memcpy (best_x, x0_p, element_size); } best_E = E; T = params.t_initial; - done = 0; + T_factor = 1.0 / params.mu_t; if (print_position) { - printf ("#-iter #-evals temperature position energy\n"); + printf ("#-iter #-evals temperature position energy bestenergy\n"); } - while (!done) { + while (1) { n_accepts = 0; n_rejects = 0; n_eless = 0; + for (i = 0; i < params.iters_fixed_T; ++i) { - if (copyfunc) { - copyfunc(x, new_x); - } else { - memcpy (new_x, x, element_size); - } + + copy_state(x, new_x, element_size, copyfunc); take_step (r, new_x, params.step_size); new_E = Ef (new_x); - if(new_E <= best_E){ - if (copyfunc) { - copyfunc(new_x,best_x); - } else { - memcpy (best_x, new_x, element_size); - } - best_E=new_E; - } - ++n_evals; /* keep track of Ef() evaluations */ + /* now take the crucial step: see if the new point is accepted - or not, as determined by the boltzman probability */ + or not, as determined by the boltzmann probability */ if (new_E < E) { + + if (new_E < best_E) { + copy_state(new_x, best_x, element_size, copyfunc); + best_E = new_E; + } + /* yay! take a step */ - if (copyfunc) { - copyfunc(new_x, x); - } else { - memcpy (x, new_x, element_size); - } + copy_state(new_x, x, element_size, copyfunc); E = new_E; ++n_eless; - } else if (gsl_rng_uniform(r) < safe_exp (-(new_E - E)/(params.k * T)) ) { + + } else if (gsl_rng_uniform(r) < boltzmann(E, new_E, T, ¶ms)) { /* yay! take a step */ - if (copyfunc) { - copyfunc(new_x, x); - } else { - memcpy(x, new_x, element_size); - } + copy_state(new_x, x, element_size, copyfunc); E = new_E; ++n_accepts; + } else { ++n_rejects; } @@ -140,25 +146,21 @@ gsl_siman_solve (const gsl_rng * r, void /* 100*n_rejects/n_steps); */ printf ("%5d %7d %12g", n_iter, n_evals, T); print_position (x); - printf (" %12g\n", E); + printf (" %12g %12g\n", E, best_E); } /* apply the cooling schedule to the temperature */ /* FIXME: I should also introduce a cooling schedule for the iters */ - T /= params.mu_t; + T *= T_factor; ++n_iter; if (T < params.t_min) { - done = 1; + break; } } /* at the end, copy the result onto the initial point, so we pass it back to the caller */ - if (copyfunc) { - copyfunc(best_x, x0_p); - } else { - memcpy (x0_p, best_x, element_size); - } + copy_state(best_x, x0_p, element_size, copyfunc); if (copyfunc) { destructor(x); @@ -185,8 +187,8 @@ gsl_siman_solve_many (const gsl_rng * r, void *x, *new_x; double *energies, *probs, *sum_probs; double Ex; /* energy of the chosen point */ - double T; /* the temperature */ - int i, done; + double T, T_factor; /* the temperature and a step multiplier */ + int i; double u; /* throw the die to choose a new "x" */ int n_iter; @@ -195,19 +197,19 @@ gsl_siman_solve_many (const gsl_rng * r, printf (" delta_pos energy\n"); } - x = (void *) malloc (params.n_tries * element_size); - new_x = (void *) malloc (params.n_tries * element_size); + x = malloc (params.n_tries * element_size); + new_x = malloc (params.n_tries * element_size); energies = (double *) malloc (params.n_tries * sizeof (double)); probs = (double *) malloc (params.n_tries * sizeof (double)); sum_probs = (double *) malloc (params.n_tries * sizeof (double)); T = params.t_initial; -/* memcpy (x, x0_p, element_size); */ + T_factor = 1.0 / params.mu_t; + memcpy (x, x0_p, element_size); - done = 0; n_iter = 0; - while (!done) + while (1) { Ex = Ef (x); for (i = 0; i < params.n_tries - 1; ++i) @@ -217,12 +219,12 @@ gsl_siman_solve_many (const gsl_rng * r, memcpy ((char *)new_x + i * element_size, x, element_size); take_step (r, (char *)new_x + i * element_size, params.step_size); energies[i] = Ef ((char *)new_x + i * element_size); - probs[i] = safe_exp (-(energies[i] - Ex) / (params.k * T)); + probs[i] = boltzmann(Ex, energies[i], T, ¶ms); } /* now add in the old value of "x", so it is a contendor */ memcpy ((char *)new_x + (params.n_tries - 1) * element_size, x, element_size); energies[params.n_tries - 1] = Ex; - probs[params.n_tries - 1] = safe_exp (-(energies[i] - Ex) / (params.k * T)); + probs[params.n_tries - 1] = boltzmann(Ex, energies[i], T, ¶ms); /* now throw biased die to see which new_x[i] we choose */ sum_probs[0] = probs[0]; @@ -235,7 +237,7 @@ gsl_siman_solve_many (const gsl_rng * r, { if (u < sum_probs[i]) { - memcpy (x, (char *)new_x + i * element_size, element_size); + memcpy (x, (char *) new_x + i * element_size, element_size); break; } } @@ -245,11 +247,11 @@ gsl_siman_solve_many (const gsl_rng * r, print_position (x); printf ("\t%12g\t%12g\n", distance (x, x0_p), Ex); } - T /= params.mu_t; + T *= T_factor; ++n_iter; if (T < params.t_min) - { - done = 1; + { + break; } } diff -urp gsl-1.9-org/siman/siman_tsp.c gsl-1.9/siman/siman_tsp.c --- gsl-1.9-org/siman/siman_tsp.c 2005-06-26 16:25:35.000000000 +0300 +++ gsl-1.9/siman/siman_tsp.c 2007-05-25 21:56:30.000000000 +0300 @@ -36,8 +36,15 @@ #define MU_T 1.002 /* damping factor for temperature */ #define T_MIN 5.0e-1 -gsl_siman_params_t params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE, - K, T_INITIAL, MU_T, T_MIN}; +gsl_siman_params_t params = { + .n_tries = N_TRIES, + .iters_fixed_T = ITERS_FIXED_T, + .step_size = STEP_SIZE, + .k = K, + .t_initial = T_INITIAL, + .mu_t = MU_T, + .t_min = T_MIN +}; struct s_tsp_city { const char * name;