1 /* Copyright (c) 2008-2013 Carlos Henrique da Silva Santos
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:
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
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.
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;
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);
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) );
43 cauchy_mit = -cauchy_mit;
45 cauchy_mit = cauchy_mit + (band*0.5);
46 valor = cauchy_mit/band;
47 valor = min+(max-min)*valor;
51 // ---------------------------------------------------------------------------
53 // main Individual representation type
54 typedef struct IndividualStructure {
59 static int CompareIndividuals(void *unused, const void *a_, const void *b_) {
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);
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 */
74 nlopt_stopping* stop, /* nlopt stop condition */
75 unsigned np, /* Number of Parents */
76 unsigned no) { /* Number of Offsprings */
78 // variables from nlopt
79 nlopt_result ret = NLOPT_SUCCESS;
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. */
93 // -------------------------------
94 // controling the population size
95 // -------------------------------
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;
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;
115 // main vector of parameters to randcauchy
116 vetor[0] = 4; // ignored
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;
132 for (item=0; item<nparameters; item++) {
136 esparents[id].parameters[item] = randcauchy(vetor);
139 memcpy(esparents[0].parameters, x, nparameters * sizeof(double));
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;
151 for (item=0; item<nparameters; item++) {
155 esoffsprings[id].parameters[item] = randcauchy(vetor);
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;
166 if (*minf > esparents[id].fitness) {
167 *minf = esparents[id].fitness;
168 memcpy(x, esparents[id].parameters,
169 nparameters * sizeof(double));
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;
177 // ------------------------------------
178 // Main Loop - Generations
179 // ------------------------------------
181 // ------------------------------------
183 // ------------------------------------
184 for (id=0; id < no; id++)
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];
196 // ------------------------------------
198 // ------------------------------------
199 totalmutation = (int) ((no * nparameters) / 10);
200 if (totalmutation < 1) totalmutation = 1;
201 for (contmutation=0; contmutation < totalmutation;
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]
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;
219 if (*minf > esoffsprings[id].fitness) {
220 *minf = esoffsprings[id].fitness;
221 memcpy(x, esoffsprings[id].parameters,
222 nparameters * sizeof(double));
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;
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];
240 nlopt_qsort_r(estotal, no+np, sizeof(Individual), NULL,
242 // copy after sorting:
243 for (i=0; i < no+np; i++) {
245 esparents[i] = estotal[i];
247 esoffsprings[i-np] = estotal[i];
249 } // generations loop
252 for (id=0; id < np; id++) free(esparents[id].parameters);
253 for (id=0; id < no; id++) free(esoffsprings[id].parameters);
255 if (esparents) free(esparents);
256 if (esoffsprings) free(esoffsprings);
257 if (estotal) free(estotal);