From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 28725 invoked by alias); 25 May 2007 19:46:54 -0000 Received: (qmail 28717 invoked by uid 22791); 25 May 2007 19:46:53 -0000 X-Spam-Check-By: sourceware.org Received: from mail.cs.tut.fi (HELO mail.cs.tut.fi) (130.230.4.42) by sourceware.org (qpsmtpd/0.31) with ESMTP; Fri, 25 May 2007 19:46:50 +0000 Received: from spam2.cs.tut.fi (spam2.cs.tut.fi [130.230.4.7]) by mail.cs.tut.fi (Postfix) with ESMTP id E8E3EAF6F for ; Fri, 25 May 2007 22:46:43 +0300 (EEST) Received: from mail.cs.tut.fi ([130.230.4.42]) by spam2.cs.tut.fi (spam2.cs.tut.fi [130.230.4.7]) (amavisd-maia, port 10024) with ESMTP id 01867-13-3 for ; Fri, 25 May 2007 22:46:42 +0300 (EEST) Received: from modeemi.modeemi.cs.tut.fi (modeemi.modeemi.cs.tut.fi [130.230.72.134]) by mail.cs.tut.fi (Postfix) with ESMTP id 81727AC7C for ; Fri, 25 May 2007 22:43:34 +0300 (EEST) Received: from jolt.modeemi.cs.tut.fi (jolt.modeemi.cs.tut.fi [130.230.72.144]) by modeemi.modeemi.cs.tut.fi (Postfix) with ESMTP id 72922248BC for ; Fri, 25 May 2007 22:43:34 +0300 (EEST) Received: by jolt.modeemi.cs.tut.fi (Postfix, from userid 16311) id 57A7F4F966; Fri, 25 May 2007 22:43:34 +0300 (EEST) Date: Fri, 25 May 2007 19:46:00 -0000 To: gsl-discuss@sources.redhat.com Subject: Simulated annealing cleanup Message-ID: <20070525194334.GF24297@jolt.modeemi.cs.tut.fi> MIME-Version: 1.0 Content-Type: multipart/mixed; boundary="UugvWAfsgieZRqgk" Content-Disposition: inline User-Agent: Mutt/1.5.11 From: shd@jolt.modeemi.cs.tut.fi (Heikki Orsila) Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2007-q2/txt/msg00022.txt.bz2 --UugvWAfsgieZRqgk Content-Type: text/plain; charset=iso-8859-1 Content-Disposition: inline Content-length: 1578 Hello. I propose the attached patch to cleanup and optimize the simulated annealing algorithm. Diffstat: siman.c | 112 ++++++++++++++++++++++++++++++------------------------------ siman_tsp.c | 11 ++++- 2 files changed, 66 insertions(+), 57 deletions(-) Here's the description of changes: siman/siman.c: - remove useless "done" variable - optimization: new_E <= best_E is only necessary when new_E < E - new_E < best_E is the proper comparison, not new_E <= best_E. This reduces useless copying of state. - fix whitespace use - replace repeated occurence of: if (copyfunc) { copyfunc(src, dst); } else { memcpy(dst, src, size); } with copy_state(). - spelling: boltzman -> boltzmann - make boltzmann function an inline function -> easier to modularize the transition probability function later - Use T *= T_factor instead of T /= params.mu_t. Avoids division -> faster. - malloc() returns void * always, remove useless cast - Print the best energy so far on progress. In the TSP example it reveals that optimization continues long after the optimum point has been found. On one test the optimum was reached with just 28% of total iterations. siman/siman_tsp.c: - initialize parameters using C99 initialization style. This way params structure may change inside and siman_tsp.c need not be updated. It is also more clear. You should attribute the patch to "Heikki Orsila " if you apply it. -- Heikki Orsila Barbie's law: heikki.orsila@iki.fi "Math is hard, let's go shopping!" http://www.iki.fi/shd --UugvWAfsgieZRqgk Content-Type: text/plain; charset=iso-8859-1 Content-Disposition: attachment; filename="simulated_annealing_cleanup.diff" Content-length: 7778 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; --UugvWAfsgieZRqgk--