1 /**********************************************************************
2 ga_simplex.c
3 **********************************************************************
4
5 ga_simplex - A simplex search algorithm for comparison and local search.
6 Copyright ©2002-2005, Stewart Adcock <stewart@linux-domain.com>
7 All rights reserved.
8
9 The latest version of this program should be available at:
10 http://gaul.sourceforge.net/
11
12 This program is free software; you can redistribute it and/or modify
13 it under the terms of the GNU General Public License as published by
14 the Free Software Foundation; either version 2 of the License, or
15 (at your option) any later version. Alternatively, if your project
16 is incompatible with the GPL, I will probably agree to requests
17 for permission to use the terms of any other license.
18
19 This program is distributed in the hope that it will be useful, but
20 WITHOUT ANY WARRANTY WHATSOEVER.
21
22 A full copy of the GNU General Public License should be in the file
23 "COPYING" provided with this distribution; if not, see:
24 http://www.gnu.org/
25
26 **********************************************************************
27
28 Synopsis: A simplex search algorithm for comparison and local search.
29
30 Note that, this algorithm requires that chromosomes
31 may be reversibly mapped to arrays of double-precision
32 floating-point array chromsomes. If this is not
33 possible then, hmmm, tough luck.
34
35 You might want to think carefully about your convergence
36 criteria.
37
38 References: Press, Flannery, Teukolsky, and Vetterling,
39 "Numerical Recipes in C: The Art of Scientific Computing"
40 Cambridge University Press, 2nd edition (1992) pp. 408-412.
41
42 Nelder, J.A., and Mead, R. Computer Journal, 7:308-313 (1965)
43
44 Yarbro, L.A., and Deming, S.N. Analytica Chim. Acta,
45 73:391-398 (1974)
46
47 To do: Make alpha, beta and gamma parameters.
48
49 **********************************************************************/
50
51 #include "gaul/ga_simplex.h"
52
53 /**********************************************************************
54 ga_population_set_simplex_parameters()
55 synopsis: Sets the simplex-search parameters for a population.
56 parameters: population *pop Population to set parameters of.
57 const int Number of dimensions for double array (Needn't match dimensions of chromosome.)
58 const double Initial step size.
59 const GAto_double Map chromosomal data to array of doubles.
60 const GAfrom_double Map array of doubles to chromosomal data.
61 return: none
62 last updated: 29 Mar 2004
63 **********************************************************************/
64
ga_population_set_simplex_parameters(population * pop,const int dimensions,const double step,const GAto_double to_double,const GAfrom_double from_double)65 void ga_population_set_simplex_parameters( population *pop,
66 const int dimensions,
67 const double step,
68 const GAto_double to_double,
69 const GAfrom_double from_double)
70 {
71
72 if ( !pop ) die("Null pointer to population structure passed.");
73 /*
74 if ( !to_double ) die("Null pointer to GAto_double callback passed.");
75 if ( !from_double ) die("Null pointer to GAfrom_double callback passed.");
76 */
77
78 plog( LOG_VERBOSE, "Population's simplex-search parameters set" );
79
80 if (pop->simplex_params == NULL)
81 pop->simplex_params = s_malloc(sizeof(ga_simplex_t));
82
83 pop->simplex_params->to_double = to_double;
84 pop->simplex_params->from_double = from_double;
85 pop->simplex_params->dimensions = dimensions;
86
87 pop->simplex_params->step = step; /* range: >0, 1=unit step randomisation, higher OK. */
88
89 pop->simplex_params->alpha = 1.50; /* range: 0=no extrap, 1=unit step extrap, higher OK. */
90 pop->simplex_params->beta = 0.75; /* range: 0=no contraction, 1=full contraction. */
91 pop->simplex_params->gamma = 0.25; /* range: 0=no contraction, 1=full contraction. */
92
93 return;
94 }
95
96
97 /**********************************************************************
98 ga_simplex()
99 synopsis: Performs optimisation on the passed entity by using a
100 simplistic simplex-search. The fitness evaluations
101 are performed using the standard and evaluation
102 callback mechanism.
103 The passed entity will have its data overwritten. The
104 remainder of the population will be let untouched.
105 Note that it is safe to pass a NULL initial structure,
106 in which case a random starting structure will be
107 generated, however the final solution will not be
108 available to the caller in any obvious way.
109 parameters:
110 return:
111 last updated: 18 Feb 2005
112 **********************************************************************/
113
ga_simplex(population * pop,entity * initial,const int max_iterations)114 int ga_simplex( population *pop,
115 entity *initial,
116 const int max_iterations )
117 {
118 int iteration=0; /* Current iteration number. */
119 int i, j; /* Index into putative solution array. */
120 entity **putative; /* Current working solutions. */
121 entity *new1, *new2; /* New putative solutions. */
122 entity *tmpentity; /* Used to swap working solutions. */
123 double *tmpdoubleptr; /* Used to swap working solutions. */
124 int num_points; /* Number of search points. */
125 double **putative_d, *putative_d_buffer; /* Storage for double arrays. */
126 double *average; /* Vector average of solutions. */
127 double *new1_d, *new2_d; /* New putative solutions. */
128 int first=0, last; /* Indices into solution arrays. */
129 boolean done=FALSE; /* Whether the shuffle sort is complete. */
130 boolean did_replace; /* Whether worst solution was replaced. */
131 boolean restart_needed; /* Whether the search needs restarting. */
132
133 /*
134 * Checks.
135 */
136 if (!pop) die("NULL pointer to population structure passed.");
137 if (!pop->evaluate) die("Population's evaluation callback is undefined.");
138 if (!pop->simplex_params) die("ga_population_set_simplex_params(), or similar, must be used prior to ga_simplex().");
139 if (!pop->simplex_params->to_double) die("Population's genome to double callback is undefined.");
140 if (!pop->simplex_params->from_double) die("Population's genome from double callback is undefined.");
141
142 /*
143 * Prepare working entities and double arrays.
144 * The space for the average and new arrays are allocated simultaneously.
145 */
146 num_points = pop->simplex_params->dimensions+1;
147 putative = s_malloc(sizeof(entity *)*num_points);
148 putative_d = s_malloc(sizeof(double *)*num_points);
149 putative_d_buffer = s_malloc(sizeof(double)*pop->simplex_params->dimensions*num_points*3);
150
151 putative_d[0] = putative_d_buffer;
152 average = &(putative_d_buffer[num_points*pop->simplex_params->dimensions]);
153 new1_d = &(putative_d_buffer[(num_points+1)*pop->simplex_params->dimensions]);
154 new2_d = &(putative_d_buffer[(num_points+2)*pop->simplex_params->dimensions]);
155
156 for (i=1; i<num_points; i++)
157 {
158 putative[i] = ga_get_free_entity(pop); /* The 'working' solutions. */
159 putative_d[i] = &(putative_d_buffer[i*pop->simplex_params->dimensions]);
160 }
161
162 new1 = ga_get_free_entity(pop);
163 new2 = ga_get_free_entity(pop);
164
165 /* Do we need to generate a random starting solution? */
166 if (!initial)
167 {
168 plog(LOG_VERBOSE, "Will perform simplex search with random starting solution.");
169
170 putative[0] = ga_get_free_entity(pop);
171 ga_entity_seed(pop, putative[0]);
172 initial = ga_get_free_entity(pop);
173 }
174 else
175 {
176 plog(LOG_VERBOSE, "Will perform simplex search with specified starting solution.");
177
178 putative[0] = ga_get_free_entity(pop);
179 ga_entity_copy(pop, putative[0], initial);
180 }
181
182 /*
183 * Generate sample points.
184 * Ensure that these initial solutions are scored.
185 *
186 * NOTE: Only perturb each solution by one dimension, by an
187 * amount specified by the step parameter; it might be better to perturb
188 * all dimensions and/or by a randomized amount.
189 */
190 #pragma omp parallel \
191 if (GAUL_DETERMINISTIC_OPENMP==0) \
192 shared(pop,num_points,putative_d,putative) private(i,j)
193 {
194 #pragma omp single \
195 nowait
196 pop->simplex_params->to_double(pop, putative[0], putative_d[0]);
197 pop->evaluate(pop, putative[0]);
198
199 #pragma omp for \
200 schedule(static) nowait
201 for (i=1; i<num_points; i++)
202 {
203 for (j=0; j<pop->simplex_params->dimensions; j++)
204 putative_d[i][j] = putative_d[0][j] +
205 random_double_range(-pop->simplex_params->step,pop->simplex_params->step);
206
207 pop->simplex_params->from_double(pop, putative[i], putative_d[i]);
208 pop->evaluate(pop, putative[i]);
209 }
210 } /* End of parallel block. */
211
212 /*
213 * Sort the initial solutions by fitness.
214 * We use a bi-directional bubble sort algorithm (which is
215 * called shuffle sort, apparently).
216 */
217 last = pop->simplex_params->dimensions-1;
218 while (done == FALSE && first < last)
219 {
220 for (j = last ; j > first ; j--)
221 {
222 if ( putative[j]->fitness > putative[j-1]->fitness )
223 { /* Swap! */
224 tmpentity = putative[j];
225 putative[j] = putative[j-1];
226 putative[j-1] = tmpentity;
227 tmpdoubleptr = putative_d[j];
228 putative_d[j] = putative_d[j-1];
229 putative_d[j-1] = tmpdoubleptr;
230 }
231 }
232 first++; /* The first one is definitely correct now. */
233
234 done = TRUE;
235
236 for (j = first ; j < last ; j++)
237 {
238 if ( putative[j]->fitness < putative[j+1]->fitness )
239 { /* Swap! */
240 tmpentity = putative[j];
241 putative[j] = putative[j+1];
242 putative[j+1] = tmpentity;
243 tmpdoubleptr = putative_d[j];
244 putative_d[j] = putative_d[j+1];
245 putative_d[j+1] = tmpdoubleptr;
246 done = FALSE;
247 }
248 }
249 last--; /* The last one is definitely correct now. */
250 }
251
252 plog( LOG_VERBOSE,
253 "Prior to the first iteration, the current solution has fitness score of %f",
254 putative[0]->fitness );
255
256 /*
257 * Do all the iterations:
258 *
259 * Stop when (a) max_iterations reached, or
260 * (b) "pop->iteration_hook" returns FALSE.
261 */
262 while ( (pop->iteration_hook?pop->iteration_hook(iteration, putative[0]):TRUE) &&
263 iteration<max_iterations )
264 {
265 iteration++;
266
267 /*
268 * Compute the vector average of all solutions except the least fit.
269 * Exploration will proceed along the vector from the least fit point
270 * to that vector average.
271 */
272 for (j = 0; j < pop->simplex_params->dimensions; j++)
273 {
274 average[j] = 0.0;
275
276 for (i = 0; i < num_points-1; i++)
277 {
278 average[j] += putative_d[i][j];
279 }
280
281 average[j] /= num_points-1;
282 }
283
284 /*
285 * Check for convergence and restart if needed.
286 * Reduce step, alpha, beta and gamma each time this happens.
287 */
288 restart_needed = TRUE;
289 /*printf("DEBUG: average = ");*/
290 for (j = 0; j < pop->simplex_params->dimensions; j++)
291 {
292 if ( average[j]-TINY > putative_d[pop->simplex_params->dimensions][j] ||
293 average[j]+TINY < putative_d[pop->simplex_params->dimensions][j] )
294 restart_needed = FALSE;
295
296 /*printf("%f ", average[j]/pop->simplex_params->dimensions);*/
297 }
298 /*printf("\n");*/
299
300 if (restart_needed != FALSE)
301 {
302 /*printf("DEBUG: restarting search.\n");*/
303 pop->simplex_params->step *= 0.50;
304 pop->simplex_params->alpha *= 0.75;
305 pop->simplex_params->beta *= 0.75;
306 pop->simplex_params->gamma *= 0.75;
307
308 for (i=1; i<num_points; i++)
309 {
310 for (j=0; j<pop->simplex_params->dimensions; j++)
311 putative_d[i][j] = putative_d[0][j] +
312 random_double_range(-pop->simplex_params->step,pop->simplex_params->step);
313
314 pop->simplex_params->from_double(pop, putative[i], putative_d[i]);
315 pop->evaluate(pop, putative[i]);
316 }
317 }
318
319 /*
320 * Simplex reflection - Extrapolate by a factor alpha away from worst point.
321 */
322 for (j = 0; j < pop->simplex_params->dimensions; j++)
323 {
324 new1_d[j] = (1.0 + pop->simplex_params->alpha) * average[j] -
325 pop->simplex_params->alpha * putative_d[num_points-1][j];
326 }
327
328 /*
329 * Evaluate the function at this reflected point.
330 */
331 pop->simplex_params->from_double(pop, new1, new1_d);
332 pop->evaluate(pop, new1);
333
334 if (new1->fitness > putative[0]->fitness)
335 {
336 /*
337 * The new solution is fitter than the previously fittest solution, so attempt an
338 * additional extrapolation by a factor alpha.
339 */
340 /*printf("DEBUG: new1 (%f) is fitter than p0 ( %f )\n", new1->fitness, putative[0]->fitness);*/
341
342 for (j = 0; j < pop->simplex_params->dimensions; j++)
343 new2_d[j] = (1.0 + pop->simplex_params->alpha) * new1_d[j] -
344 pop->simplex_params->alpha * putative_d[num_points-1][j];
345
346 pop->simplex_params->from_double(pop, new2, new2_d);
347 pop->evaluate(pop, new2);
348
349 if (new2->fitness > putative[0]->fitness)
350 {
351 /*
352 * This additional extrapolation succeeded, so replace the least fit solution
353 * by inserting new solution in correct position.
354 */
355 /*printf("DEBUG: new2 (%f) is fitter than p0 ( %f )\n", new2->fitness, putative[0]->fitness);*/
356
357 tmpentity = putative[pop->simplex_params->dimensions];
358 tmpdoubleptr = putative_d[pop->simplex_params->dimensions];
359
360 for (j = pop->simplex_params->dimensions; j > 0; j--)
361 {
362 putative[j]=putative[j-1];
363 putative_d[j]=putative_d[j-1];
364 }
365
366 putative[0] = new2;
367 putative_d[0] = new2_d;
368
369 new2 = tmpentity;
370 new2_d = tmpdoubleptr;
371 }
372 else
373 {
374 /*
375 * This additional extrapolation failed, so use the original
376 * reflected solution.
377 */
378 tmpentity = putative[pop->simplex_params->dimensions];
379 tmpdoubleptr = putative_d[pop->simplex_params->dimensions];
380
381 for (j = pop->simplex_params->dimensions; j > 0; j--)
382 {
383 putative[j]=putative[j-1];
384 putative_d[j]=putative_d[j-1];
385 }
386
387 putative[0] = new1;
388 putative_d[0] = new1_d;
389
390 new1 = tmpentity;
391 new1_d = tmpdoubleptr;
392 }
393 }
394 else if (new1->fitness < putative[pop->simplex_params->dimensions-1]->fitness)
395 {
396 /*
397 * The reflected point is worse than the second-least fit.
398 */
399 /*printf("DEBUG: new1 (%f) is less fit than p(n-1) ( %f )\n", new1->fitness, putative[pop->simplex_params->dimensions-1]->fitness);*/
400
401 did_replace = FALSE;
402
403 if (new1->fitness > putative[pop->simplex_params->dimensions]->fitness)
404 {
405 /*
406 * It is better than the least fit, so use it to replace the
407 * least fit.
408 */
409 /*printf("DEBUG: but fitter than p(n) ( %f )\n", putative[pop->simplex_params->dimensions]->fitness);*/
410 did_replace = TRUE;
411
412 tmpentity = putative[pop->simplex_params->dimensions];
413 tmpdoubleptr = putative_d[pop->simplex_params->dimensions];
414
415 putative[pop->simplex_params->dimensions] = new1;
416 putative_d[pop->simplex_params->dimensions] = new1_d;
417
418 new1 = tmpentity;
419 new1_d = tmpdoubleptr;
420 }
421 /*
422 * Perform a contraction of the simplex along one dimension, away from worst point.
423 */
424 for (j = 0; j < pop->simplex_params->dimensions; j++)
425 new1_d[j] = (1.0 - pop->simplex_params->beta) * average[j] +
426 pop->simplex_params->beta * putative_d[num_points-1][j];
427
428 pop->simplex_params->from_double(pop, new1, new1_d);
429 pop->evaluate(pop, new1);
430
431 if (new1->fitness > putative[pop->simplex_params->dimensions]->fitness)
432 {
433 /*
434 * The contraction gave an improvement, so accept it by
435 * inserting the new solution at the correct position.
436 */
437 did_replace = TRUE;
438
439 /*printf("DEBUG: contracted new1 (%f) is fitter than p(n) ( %f )\n", new1->fitness, putative[pop->simplex_params->dimensions]->fitness);*/
440
441 i = 0;
442 while (putative[i]->fitness > new1->fitness) i++;
443
444 tmpentity = putative[pop->simplex_params->dimensions];
445 tmpdoubleptr = putative_d[pop->simplex_params->dimensions];
446
447 for (j = pop->simplex_params->dimensions; j > i; j--)
448 {
449 putative[j]=putative[j-1];
450 putative_d[j]=putative_d[j-1];
451 }
452
453 putative[i] = new1;
454 putative_d[i] = new1_d;
455
456 new1 = tmpentity;
457 new1_d = tmpdoubleptr;
458 }
459
460 if (did_replace == FALSE)
461 {
462 /*
463 * The new solution is worse than the previous worse. So, contract
464 * toward the average point.
465 */
466 /*printf("DEBUG: new1 (%f) is worse than all.\n", new1->fitness);*/
467
468 for (i = 1; i < num_points; i++)
469 {
470 for (j = 0; j < pop->simplex_params->dimensions; j++)
471 putative_d[i][j] = average[j] +
472 pop->simplex_params->gamma * (putative_d[i][j] - average[j]);
473
474 pop->simplex_params->from_double(pop, putative[i], putative_d[i]);
475 pop->evaluate(pop, putative[i]);
476 }
477
478 /*
479 * Alternative is to contact toward the most fit point.
480 for (i = 1; i < num_points; i++)
481 {
482 for (j = 0; j < pop->simplex_params->dimensions; j++)
483 putative_d[i][j] = putative_d[0][j] +
484 pop->simplex_params->gamma * (putative_d[i][j] - putative_d[0][j]);
485
486 pop->simplex_params->from_double(pop, putative[i], putative_d[i]);
487 pop->evaluate(pop, putative[i]);
488 }
489 */
490
491 }
492
493 }
494 else
495 {
496 /*
497 * The reflection gave a solution which was better than the worst two
498 * solutions, but worse than the best solution.
499 * Replace the old worst solution by inserting the new solution at the
500 * correct position.
501 */
502 /*printf("DEBUG: new1 (%f) is fitter than worst 2\n", new1->fitness);
503 for (j=0; j < pop->simplex_params->dimensions; j++)
504 printf("%d fitness = %f\n", j, putative[j]->fitness);
505 */
506
507 i = 0;
508 while (putative[i]->fitness > new1->fitness) i++;
509
510 /*printf("DEBUG: new1 inserted at position %d\n", i);*/
511
512 tmpentity = putative[pop->simplex_params->dimensions];
513 tmpdoubleptr = putative_d[pop->simplex_params->dimensions];
514
515 for (j = pop->simplex_params->dimensions; j > i; j--)
516 {
517 putative[j]=putative[j-1];
518 putative_d[j]=putative_d[j-1];
519 }
520
521 putative[i] = new1;
522 putative_d[i] = new1_d;
523
524 new1 = tmpentity;
525 new1_d = tmpdoubleptr;
526 }
527
528 /*
529 * Use the iteration callback.
530 */
531 plog( LOG_VERBOSE,
532 "After iteration %d, the current solution has fitness score of %f",
533 iteration,
534 putative[0]->fitness );
535
536 } /* Iteration loop. */
537
538 /*
539 * Store best solution.
540 */
541 ga_entity_copy(pop, initial, putative[0]);
542
543 /*
544 * Cleanup.
545 */
546 ga_entity_dereference(pop, new1);
547 ga_entity_dereference(pop, new2);
548
549 for (i=0; i<num_points; i++)
550 {
551 ga_entity_dereference(pop, putative[i]);
552 }
553
554 s_free(putative);
555 s_free(putative_d);
556 s_free(putative_d_buffer);
557
558 return iteration;
559 }
560
561
562 /**********************************************************************
563 ga_simplex_double()
564 synopsis: Performs optimisation on the passed entity by using a
565 simplistic simplex-search. The fitness evaluations
566 are performed using the standard and evaluation
567 callback mechanism.
568 The passed entity will have its data overwritten. The
569 remainder of the population will be let untouched.
570 Note that it is safe to pass a NULL initial structure,
571 in which case a random starting structure will be
572 generated, however the final solution will not be
573 available to the caller in any obvious way.
574
575 Only double chromosomes may be used in this optimised
576 version of the algorithm.
577 parameters:
578 return:
579 last updated: 13 Apr 2005
580 **********************************************************************/
581
ga_simplex_double(population * pop,entity * initial,const int max_iterations)582 int ga_simplex_double( population *pop,
583 entity *initial,
584 const int max_iterations )
585 {
586 int iteration=0; /* Current iteration number. */
587 int i, j; /* Index into putative solution array. */
588 entity **putative; /* Current working solutions. */
589 entity *new1, *new2; /* New putative solutions. */
590 entity *tmpentity; /* Used to swap working solutions. */
591 int num_points; /* Number of search points. */
592 double *average; /* Vector average of solutions. */
593 int first=0, last; /* Indices into solution arrays. */
594 boolean done=FALSE; /* Whether the shuffle sort is complete. */
595 boolean did_replace; /* Whether worst solution was replaced. */
596 boolean restart_needed; /* Whether the search needs restarting. */
597
598 /*
599 * Checks.
600 */
601 if (!pop) die("NULL pointer to population structure passed.");
602 if (!pop->evaluate) die("Population's evaluation callback is undefined.");
603 if (!pop->simplex_params) die("ga_population_set_simplex_params(), or similar, must be used prior to ga_simplex().");
604
605 /*
606 * Prepare working entities and double arrays.
607 * The space for the average and new arrays are allocated simultaneously.
608 */
609 num_points = pop->len_chromosomes+1;
610 putative = s_malloc(sizeof(entity *)*num_points);
611 average = s_malloc(sizeof(double)*pop->len_chromosomes);
612
613 for (i=1; i<num_points; i++)
614 {
615 putative[i] = ga_get_free_entity(pop); /* The 'working' solutions. */
616 }
617
618 new1 = ga_get_free_entity(pop);
619 new2 = ga_get_free_entity(pop);
620
621 /* Do we need to generate a random starting solution? */
622 if (!initial)
623 {
624 plog(LOG_VERBOSE, "Will perform simplex search with random starting solution.");
625
626 putative[0] = ga_get_free_entity(pop);
627 ga_entity_seed(pop, putative[0]);
628 initial = ga_get_free_entity(pop);
629 }
630 else
631 {
632 plog(LOG_VERBOSE, "Will perform simplex search with specified starting solution.");
633
634 putative[0] = ga_get_free_entity(pop);
635 ga_entity_copy(pop, putative[0], initial);
636 }
637
638 /*
639 * Generate sample points.
640 * Ensure that these initial solutions are scored.
641 *
642 * NOTE: Only perturb each solution by one dimension, by an
643 * amount specified by the step parameter; it might be better to perturb
644 * all dimensions and/or by a randomized amount.
645 */
646 #pragma omp parallel \
647 if (GAUL_DETERMINISTIC_OPENMP==0) \
648 shared(pop,num_points,putative) private(i,j)
649 {
650 #pragma omp single \
651 nowait
652 pop->evaluate(pop, putative[0]);
653
654 #pragma omp for \
655 schedule(static) nowait
656 for (i=1; i<num_points; i++)
657 {
658 for (j=0; j<pop->len_chromosomes; j++)
659 ((double *)putative[i]->chromosome[0])[j]
660 = ((double *)putative[0]->chromosome[0])[j] +
661 random_double_range(-pop->simplex_params->step,pop->simplex_params->step);
662
663 pop->evaluate(pop, putative[i]);
664 }
665 } /* End of parallel block. */
666
667 /*
668 * Sort the initial solutions by fitness.
669 * We use a bi-directional bubble sort algorithm (which is
670 * called shuffle sort, apparently).
671 */
672 last = pop->len_chromosomes-1;
673 while (done == FALSE && first < last)
674 {
675 for (j = last ; j > first ; j--)
676 {
677 if ( putative[j]->fitness > putative[j-1]->fitness )
678 { /* Swap! */
679 tmpentity = putative[j];
680 putative[j] = putative[j-1];
681 putative[j-1] = tmpentity;
682 }
683 }
684 first++; /* The first one is definitely correct now. */
685
686 done = TRUE;
687
688 for (j = first ; j < last ; j++)
689 {
690 if ( putative[j]->fitness < putative[j+1]->fitness )
691 { /* Swap! */
692 tmpentity = putative[j];
693 putative[j] = putative[j+1];
694 putative[j+1] = tmpentity;
695 done = FALSE;
696 }
697 }
698 last--; /* The last one is definitely correct now. */
699 }
700
701 plog( LOG_VERBOSE,
702 "Prior to the first iteration, the current solution has fitness score of %f",
703 putative[0]->fitness );
704
705 /*
706 * Do all the iterations:
707 *
708 * Stop when (a) max_iterations reached, or
709 * (b) "pop->iteration_hook" returns FALSE.
710 */
711 while ( (pop->iteration_hook?pop->iteration_hook(iteration, putative[0]):TRUE) &&
712 iteration<max_iterations )
713 {
714 iteration++;
715
716 /*
717 * Compute the vector average of all solutions except the least fit.
718 * Exploration will proceed along the vector from the least fit point
719 * to that vector average.
720 */
721 for (j = 0; j < pop->len_chromosomes; j++)
722 {
723 average[j] = 0.0;
724
725 for (i = 0; i < num_points-1; i++)
726 average[j] += ((double *)putative[i]->chromosome[0])[j];
727
728 average[j] /= num_points-1;
729 }
730
731 /*
732 * Check for convergence and restart if needed.
733 * Reduce step, alpha, beta and gamma each time this happens.
734 */
735 restart_needed = TRUE;
736 /*printf("DEBUG: average = ");*/
737
738 for (j = 0; j < pop->len_chromosomes; j++)
739 {
740 if ( average[j]-TINY > ((double *)putative[pop->len_chromosomes]->chromosome[0])[j] ||
741 average[j]+TINY < ((double *)putative[pop->len_chromosomes]->chromosome[0])[j] )
742 restart_needed = FALSE;
743
744 /*printf("%f ", average[j]/pop->len_chromosomes);*/
745 }
746 /*printf("\n");*/
747
748 if (restart_needed != FALSE)
749 {
750 /*printf("DEBUG: restarting search.\n");*/
751 pop->simplex_params->step *= 0.50;
752 pop->simplex_params->alpha *= 0.75;
753 pop->simplex_params->beta *= 0.75;
754 pop->simplex_params->gamma *= 0.75;
755
756 for (i=1; i<num_points; i++)
757 {
758 for (j=0; j<pop->len_chromosomes; j++)
759 ((double *)putative[i]->chromosome[0])[j]
760 = ((double *)putative[0]->chromosome[0])[j] +
761 random_double_range(-pop->simplex_params->step,pop->simplex_params->step);
762
763 pop->evaluate(pop, putative[i]);
764 }
765 }
766
767 /*
768 * Simplex reflection - Extrapolate by a factor alpha away from worst point.
769 */
770 for (j = 0; j < pop->len_chromosomes; j++)
771 ((double *)new1->chromosome[0])[j]
772 = (1.0 + pop->simplex_params->alpha) * average[j] -
773 pop->simplex_params->alpha * ((double *)putative[num_points-1]->chromosome[0])[j];
774
775 /*
776 * Evaluate the function at this reflected point.
777 */
778 pop->evaluate(pop, new1);
779
780 if (new1->fitness > putative[0]->fitness)
781 {
782 /*
783 * The new solution is fitter than the previously fittest solution, so attempt an
784 * additional extrapolation by a factor alpha.
785 */
786 /*printf("DEBUG: new1 (%f) is fitter than p0 ( %f )\n", new1->fitness, putative[0]->fitness);*/
787
788 for (j = 0; j < pop->len_chromosomes; j++)
789 ((double *)new2->chromosome[0])[j]
790 = (1.0 + pop->simplex_params->alpha) * ((double *)new1->chromosome[0])[j] -
791 pop->simplex_params->alpha * ((double *)putative[num_points-1]->chromosome[0])[j];
792
793 pop->evaluate(pop, new2);
794
795 if (new2->fitness > putative[0]->fitness)
796 {
797 /*
798 * This additional extrapolation succeeded, so replace the least fit solution
799 * by inserting new solution in correct position.
800 */
801 /*printf("DEBUG: new2 (%f) is fitter than p0 ( %f )\n", new2->fitness, putative[0]->fitness);*/
802
803 tmpentity = putative[pop->len_chromosomes];
804
805 for (j = pop->len_chromosomes; j > 0; j--)
806 putative[j]=putative[j-1];
807
808 putative[0] = new2;
809
810 new2 = tmpentity;
811 }
812 else
813 {
814 /*
815 * This additional extrapolation failed, so use the original
816 * reflected solution.
817 */
818 tmpentity = putative[pop->len_chromosomes];
819
820 for (j = pop->len_chromosomes; j > 0; j--)
821 {
822 putative[j]=putative[j-1];
823 }
824
825 putative[0] = new1;
826
827 new1 = tmpentity;
828 }
829 }
830 else if (new1->fitness < putative[pop->len_chromosomes-1]->fitness)
831 {
832 /*
833 * The reflected point is worse than the second-least fit.
834 */
835 /*printf("DEBUG: new1 (%f) is less fit than p(n-1) ( %f )\n", new1->fitness, putative[pop->len_chromosomes-1]->fitness);*/
836
837 did_replace = FALSE;
838
839 if (new1->fitness > putative[pop->len_chromosomes]->fitness)
840 {
841 /*
842 * It is better than the least fit, so use it to replace the
843 * least fit.
844 */
845 /*printf("DEBUG: but fitter than p(n) ( %f )\n", putative[pop->len_chromosomes]->fitness);*/
846
847 did_replace = TRUE;
848
849 tmpentity = putative[pop->len_chromosomes];
850
851 putative[pop->len_chromosomes] = new1;
852
853 new1 = tmpentity;
854 }
855 /*
856 * Perform a contraction of the simplex along one dimension, away from worst point.
857 */
858 for (j = 0; j < pop->len_chromosomes; j++)
859 ((double *)new1->chromosome[0])[j]
860 = (1.0 - pop->simplex_params->beta) * average[j] +
861 pop->simplex_params->beta * ((double *)putative[num_points-1]->chromosome[0])[j];
862
863 pop->evaluate(pop, new1);
864
865 if (new1->fitness > putative[pop->len_chromosomes]->fitness)
866 {
867 /*
868 * The contraction gave an improvement, so accept it by
869 * inserting the new solution at the correct position.
870 */
871 did_replace = TRUE;
872
873 /*printf("DEBUG: contracted new1 (%f) is fitter than p(n) ( %f )\n", new1->fitness, putative[pop->len_chromosomes]->fitness);*/
874
875 i = 0;
876 while (putative[i]->fitness > new1->fitness) i++;
877
878 tmpentity = putative[pop->len_chromosomes];
879
880 for (j = pop->len_chromosomes; j > i; j--)
881 putative[j]=putative[j-1];
882
883 putative[i] = new1;
884
885 new1 = tmpentity;
886 }
887
888 if (did_replace == FALSE)
889 {
890 /*
891 * The new solution is worse than the previous worse. So, contract
892 * toward the average point.
893 */
894 /*printf("DEBUG: new1 (%f) is worse than all.\n", new1->fitness);*/
895
896 for (i = 1; i < num_points; i++)
897 {
898 for (j = 0; j < pop->len_chromosomes; j++)
899 ((double *)putative[i]->chromosome[0])[j]
900 = average[j] +
901 pop->simplex_params->gamma
902 * (((double *)putative[i]->chromosome[0])[j] - average[j]);
903
904 pop->evaluate(pop, putative[i]);
905 }
906
907 /*
908 * Alternative is to contact toward the most fit point.
909 for (i = 1; i < num_points; i++)
910 {
911 for (j = 0; j < pop->len_chromosomes; j++)
912 ((double *)putative[i]->chromosome[0])[j]
913 = ((double *)putative[0]->chromosome[0])[j] +
914 pop->simplex_params->gamma
915 * (((double *)putative[i]->chromosome[0])[j] - ((double *)putative[0]->chromosome[0])[j]);
916
917 pop->evaluate(pop, putative[i]);
918 }
919 */
920
921 }
922
923 }
924 else
925 {
926 /*
927 * The reflection gave a solution which was better than the worst two
928 * solutions, but worse than the best solution.
929 * Replace the old worst solution by inserting the new solution at the
930 * correct position.
931 */
932 /*printf("DEBUG: new1 (%f) is fitter than worst 2\n", new1->fitness);
933 for (j=0; j < pop->len_chromosomes; j++)
934 printf("%d fitness = %f\n", j, putative[j]->fitness);
935 */
936
937 i = 0;
938 while (putative[i]->fitness > new1->fitness) i++;
939
940 /*printf("DEBUG: new1 inserted at position %d\n", i);*/
941
942 tmpentity = putative[pop->len_chromosomes];
943
944 for (j = pop->len_chromosomes; j > i; j--)
945 putative[j]=putative[j-1];
946
947 putative[i] = new1;
948
949 new1 = tmpentity;
950 }
951
952 /*
953 * Use the iteration callback.
954 */
955 plog( LOG_VERBOSE,
956 "After iteration %d, the current solution has fitness score of %f",
957 iteration,
958 putative[0]->fitness );
959
960 } /* Iteration loop. */
961
962 /*
963 * Store best solution.
964 */
965 ga_entity_copy(pop, initial, putative[0]);
966
967 /*
968 * Cleanup.
969 */
970 ga_entity_dereference(pop, new1);
971 ga_entity_dereference(pop, new2);
972
973 for (i=0; i<num_points; i++)
974 {
975 ga_entity_dereference(pop, putative[i]);
976 }
977
978 s_free(putative);
979
980 return iteration;
981 }
982
983