libfgen  0.1.15
Library for optimization using a genetic algorithm or particle swarm optimization
pso.c
00001 
00002 #include <stdlib.h>
00003 #include <stdio.h>
00004 #include <string.h>
00005 #include <math.h>
00006 #include <float.h>
00007 #include "fgen.h"
00008 #include "win32_compat.h"
00009 
00010 static double *LocalBestKnownPosition(FpsoPopulation *pop, int index);
00011 
00023 FpsoPopulation *fpso_create(int population_size, int nu_parameters, FpsoGenerationCallbackFunc
00024 fpso_generation_callback_func, FpsoCalculateErrorFunc fpso_calculate_error_func) {
00025         int i;
00026         FpsoPopulation *pop = (FpsoPopulation *)malloc(sizeof(FpsoPopulation));
00027         pop->size = population_size;
00028         pop->nu_params = nu_parameters;
00029         pop->fpso_generation_callback_func = fpso_generation_callback_func;
00030         pop->fpso_calculate_error_func = fpso_calculate_error_func;
00031         pop->ind = (FpsoIndividual *)malloc(sizeof(FpsoIndividual) * population_size);
00032         for (i = 0; i < population_size; i++) {
00033                 pop->ind[i].position = (double *)malloc(sizeof(double) * nu_parameters);
00034                 pop->ind[i].velocity = (double *)malloc(sizeof(double) * nu_parameters);
00035                 pop->ind[i].best_known_position = (double *)malloc(sizeof(double) * nu_parameters);
00036         }
00037         pop->lower_bound = (double *)malloc(sizeof(double) * nu_parameters);
00038         pop->upper_bound = (double *)malloc(sizeof(double) * nu_parameters);
00039         pop->best_known_position = (double *)malloc(sizeof(double) * nu_parameters);
00040         pop->stop_signalled = 0;
00041         pop->rng = fgen_random_create_rng();
00042         return pop;
00043 }
00044 
00060 void fpso_set_parameters(FpsoPopulation *pop, int topology, int bounding_strategy, double omega, double phi1,
00061 double phi2) {
00062         pop->topology = topology;
00063         pop->bounding_strategy = bounding_strategy;
00064         pop->omega = omega;
00065         pop->phi1 = phi1;
00066         pop->phi2 = phi2;
00067 }
00068 
00074 void fpso_destroy(FpsoPopulation *pop) {
00075         int i;
00076         fgen_random_destroy_rng(pop->rng);
00077         for (i = 0; i < pop->size; i++) {
00078                 free(pop->ind[i].position);
00079                 free(pop->ind[i].velocity);
00080                 free(pop->ind[i].best_known_position);
00081         }
00082         free(pop->ind);
00083         free(pop->lower_bound);
00084         free(pop->upper_bound);
00085         free(pop->best_known_position);
00086         free(pop);
00087 }
00088 
00089 static double bound(double x, double l, double u) {
00090         if (x < l)
00091                 return l;
00092         if (x > u)
00093                 return u;
00094         return x;
00095 }
00096 
00104 void fpso_bound_position(const FpsoPopulation *pop, const double *src, double *dest) {
00105         int i;
00106         for (i = 0; i < pop->nu_params; i++)
00107                 dest[i] =
00108                     bound(src[i], pop->lower_bound[i], pop->upper_bound[i]);
00109 }
00110 
00111 static void CopyPosition(FpsoPopulation *pop, const double *src, double *dest) {
00112         memcpy(dest, src, sizeof(double) * pop->nu_params);
00113 }
00114 
00115 static void InitializePopulation(FpsoPopulation *pop) {
00116         int i;
00117         pop->best_known_error = POSITIVE_INFINITY_DOUBLE;
00118         for (i = 0; i < pop->size; i++) {
00119                 int j;
00120                 // Initialize position.
00121                 for (j = 0; j < pop->nu_params; j++)
00122                         pop->ind[i].position[j] = fgen_random_from_range_d(pop->rng, pop->lower_bound[j],
00123                                 pop->upper_bound[j]);
00124                 CopyPosition(pop, pop->ind[i].position, pop->ind[i].best_known_position);
00125                 pop->ind[i].best_known_error = pop->fpso_calculate_error_func(pop, pop->ind[i].best_known_position);
00126                 if (pop->topology == FPSO_TOPOLOGY_GBEST) {
00127                         // Update the best known position of the population.
00128                         double f = pop->ind[i].best_known_error;
00129                         if (f < pop->best_known_error) {
00130                                 pop->best_known_error = f;
00131                                 CopyPosition(pop, pop->ind[i].position, pop->best_known_position);
00132                         }
00133                 }
00134                 /* Initialize velocity. */
00135                 for (j = 0; j < pop->nu_params; j++)
00136                         pop->ind[i].velocity[j] = fgen_random_from_range_d(pop->rng,
00137                                 - fabs(pop->upper_bound[j] - pop->lower_bound[j]),
00138                                 fabs(pop->upper_bound[j] - pop->lower_bound[j]));
00139         }
00140 }
00141 
00142 static void UpdateVelocity(FpsoPopulation *pop, int i, const double *local_best_known_position) {
00143         int j;
00144         for (j = 0; j < pop->nu_params; j++) {
00145                 double r_1 = fgen_random_from_range_d(pop->rng, 0, pop->phi1);
00146                 double r_2 = fgen_random_from_range_d(pop->rng, 0, pop->phi2);
00147                 if (pop->topology == FPSO_TOPOLOGY_GBEST)
00148                         pop->ind[i].velocity[j] = pop->omega * pop->ind[i].velocity[j] +
00149                                 r_1 * (pop->ind[i].best_known_position[j] - pop->ind[i].position[j]) +
00150                                 r_2 * (pop->best_known_position[j] - pop->ind[i].position[j]);
00151                 else    // topology == FPSO_TOPOLOGY_LBEST
00152                         pop->ind[i].velocity[j] = pop->omega * pop->ind[i].velocity[j] +
00153                                 r_1 * (pop->ind[i].best_known_position[j] - pop->ind[i].position[j]) +
00154                                 r_2 * (local_best_known_position[j] - pop->ind[i].position[j]);
00155                 if (pop->bounding_strategy & FPSO_BOUND_VELOCITY_ELEMENT)
00156                         pop->ind[i].velocity[j] = bound(pop->ind[i].velocity[j],
00157                                 -fabs(pop->upper_bound[j] - pop->lower_bound[j]),
00158                                 fabs(pop->upper_bound[j] - pop->lower_bound[j]));
00159         }
00160 }
00161 
00162 static void UpdatePosition(FpsoPopulation *pop, int i) {
00163         int j;
00164         for (j = 0; j < pop->nu_params; j++) {
00165                 pop->ind[i].position[j] += pop->ind[i].velocity[j];
00166                 if (pop->bounding_strategy & FPSO_BOUND_POSITION_ELEMENT)
00167                         pop->ind[i].position[j] = bound(pop->ind[i].position[j],
00168                                   pop->lower_bound[j], pop->upper_bound[j]);
00169         }
00170 }
00171 
00172 static void UpdateLBESTBestKnownPosition(FpsoPopulation *pop, int i, double *scratch_position) {
00173         double error_position;
00174         // If the position is not bound, do bound it for the error calculation only.
00175         if (!(pop->bounding_strategy & FPSO_BOUND_POSITION_ELEMENT)) {
00176                 fpso_bound_position(pop,pop->ind[i].position, scratch_position);
00177                 error_position = pop->fpso_calculate_error_func(pop, scratch_position);
00178                 if (error_position < pop->ind[i].best_known_error) {
00179                         CopyPosition(pop, scratch_position, pop->ind[i].best_known_position);
00180                         pop->ind[i].best_known_error = error_position;
00181                 }
00182         } else {
00183                 error_position = pop->fpso_calculate_error_func(pop, pop->ind[i].position);
00184                 if (error_position < pop->ind[i].best_known_error) {
00185                         CopyPosition(pop, pop->ind[i].position, pop->ind[i].best_known_position);
00186                         pop->ind[i].best_known_error = error_position;
00187                 }
00188         }
00189 }
00190 
00191 static void UpdateGBESTBestKnownPositions(FpsoPopulation *pop, double *scratch_position) {
00192         int i;
00193         for (i = 0; i < pop->size; i++) {
00194                 double error_position;
00195                 // If the position is not bound, do bound it for the error calculation only.
00196                 if (!(pop->bounding_strategy & FPSO_BOUND_POSITION_ELEMENT)) {
00197                         fpso_bound_position(pop, pop->ind[i].position, scratch_position);
00198                         error_position = pop->fpso_calculate_error_func(pop, scratch_position);
00199                         if (error_position < pop->ind[i].best_known_error) {
00200                                 CopyPosition(pop, scratch_position, pop->ind[i].best_known_position);
00201                                 pop->ind[i].best_known_error = error_position;
00202                         }
00203                         if (error_position < pop->best_known_error) {
00204                                 pop->best_known_error = error_position;
00205                                 CopyPosition(pop, scratch_position, pop->best_known_position);
00206                         }
00207                 } else {
00208                         error_position = pop->fpso_calculate_error_func(pop, pop->ind[i].position);
00209                         if (error_position < pop->ind[i].best_known_error) {
00210                                 CopyPosition(pop, pop->ind[i].position, pop->ind[i].best_known_position);
00211                                 pop->ind[i].best_known_error = error_position;
00212                         }
00213                         if (error_position < pop->best_known_error) {
00214                                 pop->best_known_error = error_position;
00215                                 CopyPosition(pop, pop->ind[i].position, pop->best_known_position);
00216                         }
00217                 }
00218         }
00219 }
00220 
00221 static void UpdateLBESTGlobalBestKnownPosition(FpsoPopulation *pop, double *scratch_position) {
00222         int i;
00223         pop->best_known_error = POSITIVE_INFINITY_DOUBLE;
00224         for (i = 0; i < pop->size; i++)
00225                 if (pop->ind[i].best_known_error < pop->best_known_error) {
00226                         if (!(pop->bounding_strategy & FPSO_BOUND_POSITION_ELEMENT)) {
00227 //                              fpso_bound_position(pop, pop->ind[i].position, scratch_position); Wrong?
00228                                 fpso_bound_position(pop, pop->ind[i].best_known_position, scratch_position);
00229                                 CopyPosition(pop, scratch_position, pop->best_known_position);
00230                         } else
00231                                 CopyPosition(pop, pop->ind[i].best_known_position, pop->best_known_position);
00232                         pop->best_known_error = pop->ind[i].best_known_error;
00233                 }
00234 }
00235 
00236 /* Main algorithm. */
00237 
00245 void fpso_run(FpsoPopulation *pop, int max_generation) {
00246         double *scratch_position;
00247         /* Initialize the population. */
00248         InitializePopulation(pop);
00249 
00250         int generation = 0;
00251         pop->fpso_generation_callback_func(pop, 0);
00252         if (pop->stop_signalled)
00253                 return;
00254 
00255         scratch_position = (double *)malloc(sizeof(double) * pop->nu_params);
00256 
00257         /* Run the algorithm. */
00258         for (;;) {
00259                 int i;
00260                 for (i = 0; i < pop->size; i++) {
00261                         double *local_best_known_position;
00262                         if (pop->topology == FPSO_TOPOLOGY_LBEST)
00263                                 local_best_known_position = LocalBestKnownPosition(pop, i);
00264 
00265                         UpdateVelocity(pop, i, local_best_known_position);
00266 
00267                         UpdatePosition(pop, i);
00268 
00269                         /* For LBEST topology, immediately update the best known position. */
00270                         if (pop->topology == FPSO_TOPOLOGY_LBEST)
00271                                 UpdateLBESTBestKnownPosition(pop, i, scratch_position);
00272 
00273                 }
00274                 /* For GBEST topology, do sychronous updates (update best postions after all positions and */
00275                 /* velocities have been updated. */
00276                 if (pop->topology == FPSO_TOPOLOGY_GBEST)
00277                         UpdateGBESTBestKnownPositions(pop, scratch_position);
00278 
00279                 /* For LBEST topology, the best known position of the population is undefined at this stage. */
00280                 /* Before calling the generation_callback_func, update the global best known position. */
00281                 if (pop->topology == FPSO_TOPOLOGY_LBEST)
00282                         UpdateLBESTGlobalBestKnownPosition(pop, scratch_position);
00283 
00284                 generation++;
00285                 pop->fpso_generation_callback_func(pop, generation);
00286                 if ((max_generation != - 1 && generation == max_generation) || pop->stop_signalled)
00287                         break;
00288         }
00289         free(scratch_position);
00290 }
00291 
00292 double *LocalBestKnownPosition(FpsoPopulation *pop, int index) {
00293 /* Ring topology. */
00294         double *best_position = pop->ind[index].best_known_position;
00295         double best_error = pop->ind[index].best_known_error;
00296         /* Check the left neighbour. */
00297         int index2;
00298         if (index == 0)
00299                 index2 = pop->size - 1;
00300         else
00301                 index2 = index - 1;
00302         if (pop->ind[index2].best_known_error < best_error) {
00303                 best_position = pop->ind[index2].best_known_position;
00304                 best_error = pop->ind[index2].best_known_error;
00305         }
00306         /* Check the right neighbour. */
00307         int index3 = (index + 1) % pop->size;
00308         if (pop->ind[index3].best_known_error < best_error) {
00309                 best_position = pop->ind[index3].best_known_position;
00310 /*        best_error = ind[index3].best_known_error; */
00311         }
00312         return best_position;
00313 }
00314 
00315 /* Parameter setting. */
00316 
00325 void fpso_set_parameter_bounds(FpsoPopulation *pop, int index, double min_bound, double max_bound) {
00326         pop->lower_bound[index] = min_bound;
00327         pop->upper_bound[index] = max_bound;
00328 }
00329 
00336 void fpso_set_topology(FpsoPopulation *pop, int type) {
00337         pop->topology = type;
00338 }
00339 
00353 void fpso_set_bounding_strategy(FpsoPopulation *pop, int type) {
00354         pop->bounding_strategy = type;
00355 }
00356 
00357 
00362 void fpso_set_user_data(FpsoPopulation *pop, void *user_data) {
00363         pop->user_data = user_data;
00364 }
00365 
00372 double *fpso_get_best_known_position(const FpsoPopulation *pop) {
00373         return pop->best_known_position;
00374 }
00375 
00382 double fpso_get_best_known_error(const FpsoPopulation *pop) {
00383         return pop->best_known_error;
00384 }
00385 
00392 void fpso_signal_stop(FpsoPopulation *pop) {
00393         pop->stop_signalled = 1;
00394 }
00395 
00400 FgenRNG *fpso_get_rng(const FpsoPopulation *pop) {
00401         return pop->rng;
00402 }
00403 
 All Data Structures Functions Variables