2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
38 /**H*O*C**************************************************************
40 ** No.!Version! Date ! Request ! Modification ! Author **
41 ** ---+-------+------+---------+---------------------------+------- **
42 ** 1 + 3.1 +5/18/95+ - + strategy DE/rand-to-best/1+ Storn **
43 ** + + + + included + **
44 ** 1 + 3.2 +6/06/95+C.Fleiner+ change loops into memcpy + Storn **
45 ** 2 + 3.2 +6/06/95+ - + update comments + Storn **
46 ** 1 + 3.3 +6/15/95+ K.Price + strategy DE/best/2 incl. + Storn **
47 ** 2 + 3.3 +6/16/95+ - + comments and beautifying + Storn **
48 ** 3 + 3.3 +7/13/95+ - + upper and lower bound for + Storn **
49 ** + + + + initialization + **
50 ** 1 + 3.4 +2/12/96+ - + increased printout prec. + Storn **
51 ** 1 + 3.5 +5/28/96+ - + strategies revisited + Storn **
52 ** 2 + 3.5 +5/28/96+ - + strategy DE/rand/2 incl. + Storn **
53 ** 1 + 3.6 +8/06/96+ K.Price + Binomial Crossover added + Storn **
54 ** 2 + 3.6 +9/30/96+ K.Price + cost variance output + Storn **
55 ** 3 + 3.6 +9/30/96+ - + alternative to ASSIGND + Storn **
56 ** 4 + 3.6 +10/1/96+ - + variable checking inserted + Storn **
57 ** 5 + 3.6 +10/1/96+ - + strategy indic. improved + Storn **
59 ***H*O*C*E***********************************************************/
61 /* Adopted for use in GROMACS by David van der Spoel, Oct 2001 */
74 #include "gmx_fatal.h"
81 static const char *strat[] = {
82 "DE/best/1/exp", "DE/rand/1/exp",
83 "DE/rand-to-best/1/exp", "DE/best/2/exp",
84 "DE/rand/2/exp", "DE/best/1/bin",
85 "DE/rand/1/bin", "DE/rand-to-best/1/bin",
86 "DE/best/2/bin", "DE/rand/2/bin"
89 /*------------------------Macros----------------------------------------*/
90 static void assignd(int D,real a[],real b[])
98 static real **make2d(int n,int m)
109 static void copy2range(int D,real c[],t_range r[])
113 for(i=0; (i<D); i++) {
115 while (c[i] < r[i].rmin)
117 while (c[i] > r[i].rmax)
119 /* if (c[i] < r[i].rmin)
121 if (c[i] > r[i].rmax)
128 t_genalg *init_ga(FILE *fplog,const char *infile,int D,t_range range[])
135 /*------Initializations----------------------------*/
137 /*-----Read input data------------------------------------------------*/
138 fpin_ptr = gmx_fio_fopen(infile,"r");
139 if ( fscanf(fpin_ptr,"%d",&ga->NP) != 1 || /*---choice of strategy---*/
140 fscanf(fpin_ptr,"%d",&ga->strategy) != 1 || /*---choice of strategy---*/
141 fscanf(fpin_ptr,"%lf",&ff) != 1 || /*---weight factor------------*/
142 fscanf(fpin_ptr,"%lf",&cr) != 1 || /*---crossing over factor-----*/
143 fscanf(fpin_ptr,"%d",&ga->seed) != 1 || /*---random seed----------*/
144 gmx_fio_fclose(fpin_ptr) != 1)
146 gmx_fatal(FARGS,"Error reading from file %s",infile);
155 /* Allocate memory */
156 ga->pold = make2d(ga->NP,ga->D);
157 ga->pnew = make2d(ga->NP,ga->D);
159 snew(ga->best,ga->D);
160 snew(ga->bestit,ga->D);
161 snew(ga->cost,ga->NP);
162 snew(ga->msf,ga->NP);
163 snew(ga->pres,ga->NP);
164 snew(ga->scale,ga->NP);
165 snew(ga->energy,ga->NP);
167 /*-----Checking input variables for proper range--------------*/
168 if ((ga->CR < 0) || (ga->CR > 1.0))
169 gmx_fatal(FARGS,"CR=%f, should be ex [0,1]",ga->CR);
171 gmx_fatal(FARGS,"seed=%d, should be > 0",ga->seed);
172 if ((ga->strategy < 0) || (ga->strategy > 10))
173 gmx_fatal(FARGS,"strategy=%d, should be ex {1-10}",ga->strategy);
175 /* spread initial population members */
176 for (i=0; (i<ga->NP); i++) {
177 for (j=0; (j<ga->D); j++) {
178 ga->pold[i][j] = value_rand(&(range[j]),&ga->seed);
182 fprintf(fplog,"-----------------------------------------------\n");
183 fprintf(fplog,"Genetic algorithm parameters\n");
184 fprintf(fplog,"-----------------------------------------------\n");
185 fprintf(fplog,"Number of variables: %d\n",ga->D);
186 fprintf(fplog,"Population size: %d\n",ga->NP);
187 fprintf(fplog,"Strategy: %s\n",strat[ga->strategy]);
188 fprintf(fplog,"Weight factor: %g\n",ga->FF);
189 fprintf(fplog,"Crossing over factor: %g\n",ga->CR);
190 fprintf(fplog,"Random seed: %d\n",ga->seed);
191 fprintf(fplog,"-----------------------------------------------\n");
196 void update_ga(FILE *fpout_ptr,t_range range[],t_genalg *ga)
198 static int i_init=0; /* Initialisation related stuff */
199 int i, j, L, n; /* counting variables */
200 int r1,r2,r3,r4,r5; /* placeholders for random indexes */
202 if (i_init < ga->NP) {
203 /* Copy data for first force evaluation to range array */
204 copy2range(ga->D,ga->pold[i_init],range);
210 /* Now starts real genetic stuff, a new trial set is made */
211 if (ga->ipop == ga->NP) {
218 do { /* Pick a random population member */
219 /* Endless loop for ga->NP < 2 !!! */
220 r1 = (int)(rando(&ga->seed)*ga->NP);
223 do { /* Pick a random population member */
224 /* Endless loop for ga->NP < 3 !!! */
225 r2 = (int)(rando(&ga->seed)*ga->NP);
226 } while((r2==i) || (r2==r1));
229 /* Pick a random population member */
230 /* Endless loop for ga->NP < 4 !!! */
231 r3 = (int)(rando(&ga->seed)*ga->NP);
232 } while((r3==i) || (r3==r1) || (r3==r2));
235 /* Pick a random population member */
236 /* Endless loop for ga->NP < 5 !!! */
237 r4 = (int)(rando(&ga->seed)*ga->NP);
238 } while((r4==i) || (r4==r1) || (r4==r2) || (r4==r3));
241 /* Pick a random population member */
242 /* Endless loop for ga->NP < 6 !!! */
243 r5 = (int)(rando(&ga->seed)*ga->NP);
244 } while((r5==i) || (r5==r1) || (r5==r2) || (r5==r3) || (r5==r4));
247 /* Choice of strategy
248 * We have tried to come up with a sensible naming-convention: DE/x/y/z
249 * DE : stands for Differential Evolution
250 * x : a string which denotes the vector to be perturbed
251 * y : number of difference vectors taken for perturbation of x
252 * z : crossover method (exp = exponential, bin = binomial)
254 * There are some simple rules which are worth following:
255 * 1) ga->FF is usually between 0.5 and 1 (in rare cases > 1)
256 * 2) ga->CR is between 0 and 1 with 0., 0.3, 0.7 and 1. being worth to
258 * 3) To start off ga->NP = 10*ga->D is a reasonable choice. Increase ga->NP if
259 * misconvergence happens.
260 * 4) If you increase ga->NP, ga->FF usually has to be decreased
261 * 5) When the DE/ga->best... schemes fail DE/rand... usually works and
263 * EXPONENTIAL ga->CROSSOVER
264 *-------DE/ga->best/1/exp-------
265 *-------Our oldest strategy but still not bad. However, we have found several
266 *-------optimization problems where misconvergence occurs.
268 assignd(ga->D,ga->tmp,ga->pold[i]);
269 n = (int)(rando(&ga->seed)*ga->D);
272 switch (ga->strategy) {
274 /* strategy DE0 (not in our paper) */
276 ga->tmp[n] = ga->bestit[n] + ga->FF*(ga->pold[r2][n]-ga->pold[r3][n]);
279 } while((rando(&ga->seed) < ga->CR) && (L < ga->D));
283 * This is one of my favourite strategies. It works especially
284 * well when the ga->bestit[]"-schemes experience misconvergence.
285 * Try e.g. ga->FF=0.7 and ga->CR=0.5 * as a first guess.
288 /* strategy DE1 in the techreport */
290 ga->tmp[n] = ga->pold[r1][n] + ga->FF*(ga->pold[r2][n]-ga->pold[r3][n]);
293 } while((rando(&ga->seed) < ga->CR) && (L < ga->D));
296 /* DE/rand-to-ga->best/1/exp
297 * This strategy seems to be one of the ga->best strategies.
298 * Try ga->FF=0.85 and ga->CR=1.
299 * If you get misconvergence try to increase ga->NP.
300 * If this doesn't help you should play around with all three
304 /* similar to DE2 but generally better */
306 ga->tmp[n] = ga->tmp[n] + ga->FF*(ga->bestit[n] - ga->tmp[n]) +
307 ga->FF*(ga->pold[r1][n]-ga->pold[r2][n]);
310 } while((rando(&ga->seed) < ga->CR) && (L < ga->D));
313 /* DE/ga->best/2/exp is another powerful strategy worth trying */
316 ga->tmp[n] = ga->bestit[n] +
317 (ga->pold[r1][n]+ga->pold[r2][n]-ga->pold[r3][n]-ga->pold[r4][n])*ga->FF;
320 } while((rando(&ga->seed) < ga->CR) && (L < ga->D));
323 /*----DE/rand/2/exp seems to be a robust optimizer for many functions-----*/
326 ga->tmp[n] = ga->pold[r5][n] +
327 (ga->pold[r1][n]+ga->pold[r2][n]-ga->pold[r3][n]-ga->pold[r4][n])*ga->FF;
330 } while((rando(&ga->seed) < ga->CR) && (L < ga->D));
333 /*===Essentially same strategies but BINOMIAL ga->CROSSOVER===*/
335 /*-------DE/ga->best/1/bin------*/
337 for (L=0; L<ga->D; L++) {
338 /* perform D binomial trials */
339 if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1))) {
340 /* change at least one parameter */
341 ga->tmp[n] = ga->bestit[n] + ga->FF*(ga->pold[r2][n]-ga->pold[r3][n]);
347 /*-------DE/rand/1/bin------*/
349 for (L=0; L<ga->D; L++) {
350 /* perform D binomial trials */
351 if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1))) {
352 /* change at least one parameter */
353 ga->tmp[n] = ga->pold[r1][n] + ga->FF*(ga->pold[r2][n]-ga->pold[r3][n]);
359 /*-------DE/rand-to-ga->best/1/bin------*/
361 for (L=0; L<ga->D; L++) {
362 /* perform ga->D binomial trials */
363 if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1))) {
364 /* change at least one parameter */
365 ga->tmp[n] = ga->tmp[n] + ga->FF*(ga->bestit[n] - ga->tmp[n]) +
366 ga->FF*(ga->pold[r1][n]-ga->pold[r2][n]);
372 /*-------DE/ga->best/2/bin------*/
374 for (L=0; L<ga->D; L++) {
375 /* perform ga->D binomial trials */
376 if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1))) {
377 /* change at least one parameter */
378 ga->tmp[n] = ga->bestit[n] +
379 (ga->pold[r1][n]+ga->pold[r2][n]-ga->pold[r3][n]-ga->pold[r4][n])*ga->FF;
385 /*-------DE/rand/2/bin-------*/
387 for (L=0; L<ga->D; L++) {
388 /* perform ga->D binomial trials */
389 if ((rando(&ga->seed) < ga->CR) || (L == (ga->D-1))) {
390 /* change at least one parameter */
391 ga->tmp[n] = ga->pold[r5][n] +
392 (ga->pold[r1][n]+ga->pold[r2][n]-ga->pold[r3][n]-ga->pold[r4][n])*ga->FF;
399 /*===Trial mutation now in ga->tmp[]. Test how good this choice really was.==*/
400 copy2range(ga->D,ga->tmp,range);
404 gmx_bool print_ga(FILE *fp,t_genalg *ga,real msf,tensor pres,rvec scale,
405 real energy,t_range range[],real tol)
407 static int nfeval=0; /* number of function evaluations */
408 static gmx_bool bImproved;
410 real cvar; /* computes the cost variance */
411 real cmean; /* mean cost */
415 trial_cost = cost(pres,msf,energy);
416 if (nfeval < ga->NP) {
417 ga->cost[nfeval] = trial_cost;
418 ga->msf[nfeval] = msf;
419 ga->energy[nfeval] = energy;
420 copy_mat(pres,ga->pres[nfeval]);
421 copy_rvec(scale,ga->scale[nfeval]);
423 pr_rvec(debug,0,"scale",scale,DIM,TRUE);
424 pr_rvec(debug,0,"pold ",ga->pold[nfeval]+4,DIM,TRUE);
429 /* When we get here we have done an initial evaluation for all
430 * animals in the population
435 /* First iteration after first round of trials */
436 if (nfeval == ga->NP) {
437 /* Evaluate who is ga->best */
439 for (j=1; (j<ga->NP); j++) {
440 if (ga->cost[j] < ga->cost[ga->imin])
443 assignd(ga->D,ga->best,ga->pold[ga->imin]);
444 /* save best member ever */
445 assignd(ga->D,ga->bestit,ga->pold[ga->imin]);
446 /* save best member of generation */
449 if (trial_cost < ga->cost[ga->ipop]) {
450 if (trial_cost < ga->cost[ga->imin]) {
451 /* Was this a new minimum? */
452 /* if so, reset cmin to new low...*/
454 assignd(ga->D,ga->best,ga->tmp);
457 /* improved objective function value ? */
458 ga->cost[ga->ipop] = trial_cost;
460 ga->msf[ga->ipop] = msf;
461 ga->energy[ga->ipop] = energy;
462 copy_mat(pres,ga->pres[ga->ipop]);
463 copy_rvec(scale,ga->scale[ga->ipop]);
465 assignd(ga->D,ga->pnew[ga->ipop],ga->tmp);
469 /* replace target with old value */
470 assignd(ga->D,ga->pnew[ga->ipop],ga->pold[ga->ipop]);
472 /* #define SCALE_DEBUG*/
476 rvec_sub(ga->scale[ga->imin],&(ga->best[ga->D-3]),dscale);
477 if (norm(dscale) > 0) {
478 pr_rvec(fp,0,"scale",scale,DIM,TRUE);
479 pr_rvec(fp,0,"best ",&(ga->best[ga->D-3]),DIM,TRUE);
480 fprintf(fp,"imin = %d, ipop = %d, nfeval = %d\n",ga->imin,
482 gmx_fatal(FARGS,"Scale inconsistency");
487 /* Increase population member count */
490 /* End mutation loop through population? */
491 if (ga->ipop == ga->NP) {
492 /* Save ga->best population member of current iteration */
493 assignd(ga->D,ga->bestit,ga->best);
495 /* swap population arrays. New generation becomes old one */
500 /*----Compute the energy variance (just for monitoring purposes)-----------*/
501 /* compute the mean value first */
503 for (j=0; (j<ga->NP); j++)
504 cmean += ga->cost[j];
505 cmean = cmean/ga->NP;
507 /* now the variance */
509 for (j=0; (j<ga->NP); j++)
510 cvar += sqr(ga->cost[j] - cmean);
511 cvar = cvar/(ga->NP-1);
513 /*----Output part----------------------------------------------------------*/
514 if (1 || bImproved || (nfeval == ga->NP)) {
515 fprintf(fp,"\nGen: %5d\n Cost: %12.3e <Cost>: %12.3e\n"
516 " Ener: %8.4f RMSF: %8.3f Pres: %8.1f %8.1f %8.1f (%8.1f)\n"
517 " Box-Scale: %15.10f %15.10f %15.10f\n",
518 ga->gen,ga->cost[ga->imin],cmean,ga->energy[ga->imin],
519 sqrt(ga->msf[ga->imin]),ga->pres[ga->imin][XX][XX],
520 ga->pres[ga->imin][YY][YY],ga->pres[ga->imin][ZZ][ZZ],
521 trace(ga->pres[ga->imin])/3.0,
522 ga->scale[ga->imin][XX],ga->scale[ga->imin][YY],
523 ga->scale[ga->imin][ZZ]);
525 for (j=0; (j<ga->D); j++)
526 fprintf(fp,"\tbest[%d]=%-15.10f\n",j,ga->best[j]);
527 if (ga->cost[ga->imin] < tol) {
528 for (i=0; (i<ga->NP); i++) {
529 fprintf(fp,"Animal: %3d Cost:%8.3f Ener: %8.3f RMSF: %8.3f%s\n",
530 i,ga->cost[i],ga->energy[i],sqrt(ga->msf[i]),
531 (i == ga->imin) ? " ***" : "");
532 for (j=0; (j<ga->D); j++)
533 fprintf(fp,"\tParam[%d][%d]=%15.10g\n",i,j,ga->pold[i][j]);