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)) return NLOPT_INVALID_ARGS;
  99      esparents    = (Individual*) malloc(sizeof(Individual) * np);
 100      esoffsprings = (Individual*) malloc(sizeof(Individual) * no);
 101      estotal     = (Individual*) malloc(sizeof(Individual) * (np+no));
 102      if ((!esparents)||(!esoffsprings)||(!estotal)) {
 103           free(esparents); free(esoffsprings); free(estotal);
 104           return NLOPT_OUT_OF_MEMORY;
 105      }
 106      for (id=0; id < np; id++) esparents[id].parameters = NULL;
 107      for (id=0; id < no; id++) esoffsprings[id].parameters = NULL;
 108      // From here the population is initialized
 109      /* we don't handle unbounded search regions;
 110         this check is unnecessary since it is performed in nlopt_optimize.
 111         for (j = 0; j < nparameters; ++j)
 112           if (nlopt_isinf(low[j]) || nlopt_isinf(up[j]))
 113             return NLOPT_INVALID_ARGS;
 114      */
 115      // main vector of parameters to randcauchy
 116      vetor[0] = 4; // ignored
 117      vetor[3] = 0;
 118      vetor[4] = 1;
 119      vetor[5] = 10;
 120      vetor[6] = 1;
 121      vetor[7] = 0; // ignored
 122      // ------------------------------------
 123      // Initializing parents population
 124      // ------------------------------------
 125      for (id=0; id < np; id++) {
 126           esparents[id].parameters =
 127                (double*) malloc(sizeof(double) * nparameters);
 128           if (!esparents[id].parameters) {
 129                ret = NLOPT_OUT_OF_MEMORY;
 130                goto done;
 131           }
 132           for (item=0; item<nparameters; item++) {
 133                vetor[1] = lb[item];
 134                vetor[2] = ub[item];
 135                vetor[7]=vetor[7]+1;
 136                esparents[id].parameters[item] = randcauchy(vetor);
 137           }
 138      }
 139      memcpy(esparents[0].parameters, x, nparameters * sizeof(double));
 140
 141      // ------------------------------------
 142      // Initializing offsprings population
 143      // ------------------------------------
 144      for (id=0; id < no; id++) {
 145           esoffsprings[id].parameters =
 146                (double*) malloc(sizeof(double) * nparameters);
 147           if (!esoffsprings[id].parameters) {
 148                ret = NLOPT_OUT_OF_MEMORY;
 149                goto done;
 150           }
 151           for (item=0; item<nparameters; item++) {
 152                vetor[1] = lb[item];
 153                vetor[2] = ub[item];
 154                vetor[7]=vetor[7]+1;
 155                esoffsprings[id].parameters[item] = randcauchy(vetor);
 156           }
 157      }
 158      // ------------------------------------
 159      // Parents fitness evaluation
 160      // ------------------------------------
 161      for (id=0; id < np; id++) {
 162           esparents[id].fitness =
 163                f(nparameters, esparents[id].parameters, NULL, data_f);
 164           estotal[id].fitness = esparents[id].fitness;
 165           stop->nevals++;
 166           if (*minf > esparents[id].fitness) {
 167                *minf = esparents[id].fitness;
 168                memcpy(x, esparents[id].parameters,
 169                       nparameters * sizeof(double));
 170           }
 171           if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 172           else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
 173           else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 174           else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 175           if (ret != NLOPT_SUCCESS) goto done;
 176      }
 177      // ------------------------------------
 178      // Main Loop - Generations
 179      // ------------------------------------
 180      while (1) {
 181           // ------------------------------------
 182           // Crossover
 183           // ------------------------------------
 184           for (id=0; id < no; id++)
 185           {
 186                parent1  = nlopt_iurand((int) np);
 187                parent2  = nlopt_iurand((int) np);
 188                crosspoint = (unsigned) nlopt_iurand((int) nparameters);
 189                for (item=0; item < crosspoint; item++)
 190                     esoffsprings[id].parameters[item]
 191                          = esparents[parent1].parameters[item];
 192                for (item=crosspoint; item < nparameters; item++)
 193                     esoffsprings[id].parameters[item]
 194                          = esparents[parent2].parameters[item];
 195           }
 196           // ------------------------------------
 197           // Gaussian Mutation
 198           // ------------------------------------
 199           totalmutation = (int) ((no * nparameters) / 10);
 200           if (totalmutation < 1) totalmutation = 1;
 201           for (contmutation=0; contmutation < totalmutation;
 202                contmutation++) {
 203                idoffmutation = nlopt_iurand((int) no);
 204                paramoffmutation = nlopt_iurand((int) nparameters);
 205                vetor[1] = lb[paramoffmutation];
 206                vetor[2] = ub[paramoffmutation];
 207                vetor[7] = vetor[7]+contmutation;
 208                esoffsprings[idoffmutation].parameters[paramoffmutation]
 209                     = randcauchy(vetor);
 210           }
 211           // ------------------------------------
 212           // Offsprings fitness evaluation
 213           // ------------------------------------
 214           for (id=0; id < no; id++){
 215                //esoffsprings[id].fitness = (double)fitness(esoffsprings[id].parameters, nparameters,fittype);
 216                esoffsprings[id].fitness = f(nparameters, esoffsprings[id].parameters, NULL, data_f);
 217                estotal[id+np].fitness = esoffsprings[id].fitness;
 218                stop->nevals++;
 219                if (*minf > esoffsprings[id].fitness) {
 220                     *minf = esoffsprings[id].fitness;
 221                     memcpy(x, esoffsprings[id].parameters,
 222                            nparameters * sizeof(double));
 223                }
 224                if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 225                else if (*minf < stop->minf_max)
 226                     ret = NLOPT_MINF_MAX_REACHED;
 227                else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 228                else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
 229                if (ret != NLOPT_SUCCESS) goto done;
 230           }
 231           // ------------------------------------
 232           // Individual selection
 233           // ------------------------------------
 234           // all the individuals are copied to one vector to easily identify best solutions
 235           for (i=0; i < np; i++)
 236                estotal[i] = esparents[i];
 237           for (i=0; i < no; i++)
 238                estotal[np+i] = esoffsprings[i];
 239           // Sorting
 240           nlopt_qsort_r(estotal, no+np, sizeof(Individual), NULL,
 241                         CompareIndividuals);
 242           // copy after sorting:
 243           for (i=0; i < no+np; i++) {
 244                if (i<np)
 245                     esparents[i] = estotal[i];
 246                else
 247                     esoffsprings[i-np] = estotal[i];
 248           }
 249      } // generations loop
 250
 251 done:
 252      for (id=0; id < np; id++) free(esparents[id].parameters);
 253      for (id=0; id < no; id++) free(esoffsprings[id].parameters);
 254
 255      if (esparents)     free(esparents);
 256      if (esoffsprings)  free(esoffsprings);
 257      if (estotal)               free(estotal);
 258      return ret;
 259 }