diff options
Diffstat (limited to 'gsl-1.9/siman/siman.c')
| -rw-r--r-- | gsl-1.9/siman/siman.c | 266 |
1 files changed, 266 insertions, 0 deletions
diff --git a/gsl-1.9/siman/siman.c b/gsl-1.9/siman/siman.c new file mode 100644 index 0000000..051c235 --- /dev/null +++ b/gsl-1.9/siman/siman.c @@ -0,0 +1,266 @@ +/* siman/siman.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <stdio.h> +#include <math.h> +#include <stdlib.h> +#include <string.h> +#include <assert.h> + +#include <gsl/gsl_machine.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_siman.h> + +static inline +double safe_exp (double x) /* avoid underflow errors for large uphill steps */ +{ + return (x < GSL_LOG_DBL_MIN) ? 0.0 : exp(x); +} + +/* implementation of a basic simulated annealing algorithm */ + +void +gsl_siman_solve (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef, + gsl_siman_step_t take_step, + gsl_siman_metric_t distance, + gsl_siman_print_t print_position, + gsl_siman_copy_t copyfunc, + gsl_siman_copy_construct_t copy_constructor, + gsl_siman_destroy_t destructor, + size_t element_size, + gsl_siman_params_t params) +{ + void *x, *new_x, *best_x; + double E, new_E, best_E; + int i, done; + double T; + int n_evals = 1, n_iter = 0, n_accepts, n_rejects, n_eless; + + /* this function requires that either the dynamic functions (copy, + copy_constructor and destrcutor) are passed, or that an element + size is given */ + assert((copyfunc != NULL && copy_constructor != NULL && destructor != NULL) + || (element_size != 0)); + + distance = 0 ; /* This parameter is not currently used */ + E = Ef(x0_p); + + if (copyfunc) { + x = copy_constructor(x0_p); + new_x = copy_constructor(x0_p); + best_x = copy_constructor(x0_p); + } else { + x = (void *) malloc (element_size); + memcpy (x, x0_p, element_size); + new_x = (void *) malloc (element_size); + best_x = (void *) malloc (element_size); + memcpy (best_x, x0_p, element_size); + } + + best_E = E; + + T = params.t_initial; + done = 0; + + if (print_position) { + printf ("#-iter #-evals temperature position energy\n"); + } + + while (!done) { + + 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); + } + + 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 */ + if (new_E < E) { + /* yay! take a step */ + if (copyfunc) { + copyfunc(new_x, x); + } else { + memcpy (x, new_x, element_size); + } + E = new_E; + ++n_eless; + } else if (gsl_rng_uniform(r) < safe_exp (-(new_E - E)/(params.k * T)) ) { + /* yay! take a step */ + if (copyfunc) { + copyfunc(new_x, x); + } else { + memcpy(x, new_x, element_size); + } + E = new_E; + ++n_accepts; + } else { + ++n_rejects; + } + } + + if (print_position) { + /* see if we need to print stuff as we go */ + /* printf("%5d %12g %5d %3d %3d %3d", n_iter, T, n_evals, */ + /* 100*n_eless/n_steps, 100*n_accepts/n_steps, */ + /* 100*n_rejects/n_steps); */ + printf ("%5d %7d %12g", n_iter, n_evals, T); + print_position (x); + printf (" %12g\n", E); + } + + /* apply the cooling schedule to the temperature */ + /* FIXME: I should also introduce a cooling schedule for the iters */ + T /= params.mu_t; + ++n_iter; + if (T < params.t_min) { + done = 1; + } + } + + /* 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); + } + + if (copyfunc) { + destructor(x); + destructor(new_x); + destructor(best_x); + } else { + free (x); + free (new_x); + free (best_x); + } +} + +/* implementation of a simulated annealing algorithm with many tries */ + +void +gsl_siman_solve_many (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef, + gsl_siman_step_t take_step, + gsl_siman_metric_t distance, + gsl_siman_print_t print_position, + size_t element_size, + gsl_siman_params_t params) +{ + /* the new set of trial points, and their energies and probabilities */ + void *x, *new_x; + double *energies, *probs, *sum_probs; + double Ex; /* energy of the chosen point */ + double T; /* the temperature */ + int i, done; + double u; /* throw the die to choose a new "x" */ + int n_iter; + + if (print_position) { + printf ("#-iter temperature position"); + printf (" delta_pos energy\n"); + } + + x = (void *) malloc (params.n_tries * element_size); + new_x = (void *) 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); */ + memcpy (x, x0_p, element_size); + done = 0; + + n_iter = 0; + while (!done) + { + Ex = Ef (x); + for (i = 0; i < params.n_tries - 1; ++i) + { /* only go to N_TRIES-2 */ + /* center the new_x[] around x, then pass it to take_step() */ + sum_probs[i] = 0; + 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)); + } + /* 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)); + + /* now throw biased die to see which new_x[i] we choose */ + sum_probs[0] = probs[0]; + for (i = 1; i < params.n_tries; ++i) + { + sum_probs[i] = sum_probs[i - 1] + probs[i]; + } + u = gsl_rng_uniform (r) * sum_probs[params.n_tries - 1]; + for (i = 0; i < params.n_tries; ++i) + { + if (u < sum_probs[i]) + { + memcpy (x, (char *)new_x + i * element_size, element_size); + break; + } + } + if (print_position) + { + printf ("%5d\t%12g\t", n_iter, T); + print_position (x); + printf ("\t%12g\t%12g\n", distance (x, x0_p), Ex); + } + T /= params.mu_t; + ++n_iter; + if (T < params.t_min) + { + done = 1; + } + } + + /* now return the value via x0_p */ + memcpy (x0_p, x, element_size); + + /* printf("the result is: %g (E=%g)\n", x, Ex); */ + + free (x); + free (new_x); + free (energies); + free (probs); + free (sum_probs); +} |