esch/esch.c

   1 /* Copyright (c) 2008-2013 Carlos Henrique da Silva Santos
   2  *
   3  * Permission is hereby granted, free of charge, to any person obtaining
   4  * a copy of this software and associated documentation files (the
   5  * "Software"), to deal in the Software without restriction, including
   6  * without limitation the rights to use, copy, modify, merge, publish,
   7  * distribute, sublicense, and/or sell copies of the Software, and to
   8  * permit persons to whom the Software is furnished to do so, subject to
   9  * the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be
  12  * included in all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  16  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  17  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  18  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  19  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  20  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  21  */
  22
  23 #include <stdlib.h>
  24 #include <string.h>
  25 #include "esch.h"
  26
  27 // ---------------------------------------------------------------------------
  28 // Cauchy random number distribution
  29 static double randcauchy(const double params[7]) {
  30      /* double min, double max, double mi, double t, double band */
  31      double na_unif, cauchy_mit, limit_inf, limit_sup;
  32      double valor;
  33      double min = params[1], max = params[2], mi = params[3],
  34           t = params[4], band = params[5];
  35      limit_inf = mi - (band*0.5);
  36      limit_sup = mi + (band*0.5);
  37      do {
  38           na_unif = nlopt_urand(0,1); // ran2(0,1);
  39           cauchy_mit = t*tan((na_unif-0.5)*3.14159265358979323846) + mi;
  40      } while ( (cauchy_mit<limit_inf) || (cauchy_mit>limit_sup) );
  41
  42      if (cauchy_mit < 0)
  43           cauchy_mit = -cauchy_mit;
  44      else
  45           cauchy_mit = cauchy_mit + (band*0.5);
  46      valor  = cauchy_mit/band;
  47      valor = min+(max-min)*valor;
  48      return valor;
  49 }
  50
  51 // ---------------------------------------------------------------------------
  52
  53 // main Individual representation type
  54 typedef struct IndividualStructure {
  55      double * parameters;
  56      double fitness;
  57 } Individual;
  58
  59 static int CompareIndividuals(void *unused, const void *a_, const void *b_) {
  60      // (void) unused;
  61      const Individual *a = (const Individual *) a_;
  62      const Individual *b = (const Individual *) b_;
  63      return a->fitness < b->fitness ? -1 : (a->fitness > b->fitness ? +1 : 0);
  64 }
  65
  66 nlopt_result chevolutionarystrategy(
  67      unsigned nparameters, /* Number of input parameters */
  68      nlopt_func f,      /* Recursive Objective Funtion Call */
  69      void * data_f,     /* Data to Objective Function */
  70      const double* lb,                  /* Lower bound values */
  71      const double* ub,                  /* Upper bound values */
  72      double* x,                         /*in: initial guess, out: minimizer */
  73      double* minf,
  74      nlopt_stopping* stop,              /* nlopt stop condition */
  75      unsigned np,                       /* Number of Parents */
  76      unsigned no) {                     /* Number of Offsprings */
  77
  78      // variables from nlopt
  79      nlopt_result ret = NLOPT_SUCCESS;
  80      double vetor[8];
  81      unsigned  i, id, item;
  82      int  parent1, parent2;
  83      unsigned crosspoint;  // crossover parameteres
  84      int  contmutation, totalmutation;  // mutation parameters
  85      int  idoffmutation, paramoffmutation;      // mutation parameters
  86      Individual * esparents;                    // Parents population
  87      Individual * esoffsprings;         // Offsprings population
  88      Individual * estotal;// copy containing Parents and Offsprings pops
  89      /* It is interesting to maintain the parents and offsprings
  90       * populations stablished and sorted; when the final iterations
  91       * is achieved, they are ranked and updated. */
  92
  93      // -------------------------------
  94      // controling the population size
  95      // -------------------------------
  96      if (!np) np = 40;
  97      if (!no) no = 60;
  98      if ((np < 1)||(no<1)) {
  99          nlopt_stop_msg(stop, "populations %d, %d are too small", np, no);
 100          return NLOPT_INVALID_ARGS;
 101      }
 102      esparents    = (Individual*) malloc(sizeof(Individual) * np);
 103      esoffsprings = (Individual*) malloc(sizeof(Individual) * no);
 104      estotal     = (Individual*) malloc(sizeof(Individual) * (np+no));
 105      if ((!esparents)||(!esoffsprings)||(!estotal)) {
 106           free(esparents); free(esoffsprings); free(estotal);
 107           return NLOPT_OUT_OF_MEMORY;
 108      }
 109      for (id=0; id < np; id++) esparents[id].parameters = NULL;
 110      for (id=0; id < no; id++) esoffsprings[id].parameters = NULL;
 111      // From here the population is initialized
 112      /* we don't handle unbounded search regions;
 113         this check is unnecessary since it is performed in nlopt_optimize.
 114         for (j = 0; j < nparameters; ++j)
 115           if (nlopt_isinf(low[j]) || nlopt_isinf(up[j]))
 116             return NLOPT_INVALID_ARGS;
 117      */
 118      // main vector of parameters to randcauchy
 119      vetor[0] = 4; // ignored
 120      vetor[3] = 0;
 121      vetor[4] = 1;
 122      vetor[5] = 10;
 123      vetor[6] = 1;
 124      vetor[7] = 0; // ignored
 125      // ------------------------------------
 126      // Initializing parents population
 127      // ------------------------------------
 128      for (id=0; id < np; id++) {
 129           esparents[id].parameters =
 130                (double*) malloc(sizeof(double) * nparameters);
 131           if (!esparents[id].parameters) {
 132                ret = NLOPT_OUT_OF_MEMORY;
 133                goto done;
 134           }
 135           for (item=0; item<nparameters; item++) {
 136                vetor[1] = lb[item];
 137                vetor[2] = ub[item];
 138                vetor[7]=vetor[7]+1;
 139                esparents[id].parameters[item] = randcauchy(vetor);
 140           }
 141      }
 142      memcpy(esparents[0].parameters, x, nparameters * sizeof(double));
 143
 144      // ------------------------------------
 145      // Initializing offsprings population
 146      // ------------------------------------
 147      for (id=0; id < no; id++) {
 148           esoffsprings[id].parameters =
 149                (double*) malloc(sizeof(double) * nparameters);
 150           if (!esoffsprings[id].parameters) {
 151                ret = NLOPT_OUT_OF_MEMORY;
 152                goto done;
 153           }
 154           for (item=0; item<nparameters; item++) {
 155                vetor[1] = lb[item];
 156                vetor[2] = ub[item];
 157                vetor[7]=vetor[7]+1;
 158                esoffsprings[id].parameters[item] = randcauchy(vetor);
 159           }
 160      }
 161      // ------------------------------------
 162      // Parents fitness evaluation
 163      // ------------------------------------
 164      for (id=0; id < np; id++) {
 165           esparents[id].fitness =
 166                f(nparameters, esparents[id].parameters, NULL, data_f);
 167           estotal[id].fitness = esparents[id].fitness;
 168           stop->nevals++;
 169           if (*minf > esparents[id].fitness) {
 170                *minf = esparents[id].fitness;
 171                memcpy(x, esparents[id].parameters,
 172                       nparameters * sizeof(double));
 173           }
 174           if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 175           else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 176           else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 177           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 178           if (ret != NLOPT_SUCCESS) goto done;
 179      }
 180      // ------------------------------------
 181      // Main Loop - Generations
 182      // ------------------------------------
 183      while (1) {
 184           // ------------------------------------
 185           // Crossover
 186           // ------------------------------------
 187           for (id=0; id < no; id++)
 188           {
 189                parent1  = nlopt_iurand((int) np);
 190                parent2  = nlopt_iurand((int) np);
 191                crosspoint = (unsigned) nlopt_iurand((int) nparameters);
 192                for (item=0; item < crosspoint; item++)
 193                     esoffsprings[id].parameters[item]
 194                          = esparents[parent1].parameters[item];
 195                for (item=crosspoint; item < nparameters; item++)
 196                     esoffsprings[id].parameters[item]
 197                          = esparents[parent2].parameters[item];
 198           }
 199           // ------------------------------------
 200           // Gaussian Mutation
 201           // ------------------------------------
 202           totalmutation = (int) ((no * nparameters) / 10);
 203           if (totalmutation < 1) totalmutation = 1;
 204           for (contmutation=0; contmutation < totalmutation;
 205                contmutation++) {
 206                idoffmutation = nlopt_iurand((int) no);
 207                paramoffmutation = nlopt_iurand((int) nparameters);
 208                vetor[1] = lb[paramoffmutation];
 209                vetor[2] = ub[paramoffmutation];
 210                vetor[7] = vetor[7]+contmutation;
 211                esoffsprings[idoffmutation].parameters[paramoffmutation]
 212                     = randcauchy(vetor);
 213           }
 214           // ------------------------------------
 215           // Offsprings fitness evaluation
 216           // ------------------------------------
 217           for (id=0; id < no; id++){
 218                //esoffsprings[id].fitness = (double)fitness(esoffsprings[id].parameters, nparameters,fittype);
 219                esoffsprings[id].fitness = f(nparameters, esoffsprings[id].parameters, NULL, data_f);
 220                estotal[id+np].fitness = esoffsprings[id].fitness;
 221                stop->nevals++;
 222                if (*minf > esoffsprings[id].fitness) {
 223                     *minf = esoffsprings[id].fitness;
 224                     memcpy(x, esoffsprings[id].parameters,
 225                            nparameters * sizeof(double));
 226                }
 227                if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 228                else if (*minf < stop->minf_max)
 229                     ret = NLOPT_MINF_MAX_REACHED;
 230                else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 231                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 232                if (ret != NLOPT_SUCCESS) goto done;
 233           }
 234           // ------------------------------------
 235           // Individual selection
 236           // ------------------------------------
 237           // all the individuals are copied to one vector to easily identify best solutions
 238           for (i=0; i < np; i++)
 239                estotal[i] = esparents[i];
 240           for (i=0; i < no; i++)
 241                estotal[np+i] = esoffsprings[i];
 242           // Sorting
 243           nlopt_qsort_r(estotal, no+np, sizeof(Individual), NULL,
 244                         CompareIndividuals);
 245           // copy after sorting:
 246           for (i=0; i < no+np; i++) {
 247                if (i<np)
 248                     esparents[i] = estotal[i];
 249                else
 250                     esoffsprings[i-np] = estotal[i];
 251           }
 252      } // generations loop
 253
 254 done:
 255      for (id=0; id < np; id++) free(esparents[id].parameters);
 256      for (id=0; id < no; id++) free(esoffsprings[id].parameters);
 257
 258      if (esparents)     free(esparents);
 259      if (esoffsprings)  free(esoffsprings);
 260      if (estotal)               free(estotal);
 261      return ret;
 262 }