1 /* Copyright 2017-2021 PaGMO development team
2 
3 This file is part of the PaGMO library.
4 
5 The PaGMO library is free software; you can redistribute it and/or modify
6 it under the terms of either:
7 
8   * the GNU Lesser General Public License as published by the Free
9     Software Foundation; either version 3 of the License, or (at your
10     option) any later version.
11 
12 or
13 
14   * the GNU General Public License as published by the Free Software
15     Foundation; either version 3 of the License, or (at your option) any
16     later version.
17 
18 or both in parallel, as here.
19 
20 The PaGMO library is distributed in the hope that it will be useful, but
21 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
22 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
23 for more details.
24 
25 You should have received copies of the GNU General Public License and the
26 GNU Lesser General Public License along with the PaGMO library.  If not,
27 see https://www.gnu.org/licenses/. */
28 
29 #include <cmath>
30 #include <iomanip>
31 #include <iostream>
32 #include <numeric>
33 #include <random>
34 #include <sstream>
35 #include <stdexcept>
36 #include <string>
37 #include <utility>
38 #include <vector>
39 
40 #include <pagmo/algorithm.hpp>
41 #include <pagmo/algorithms/de1220.hpp>
42 #include <pagmo/exceptions.hpp>
43 #include <pagmo/io.hpp>
44 #include <pagmo/population.hpp>
45 #include <pagmo/s11n.hpp>
46 #include <pagmo/types.hpp>
47 #include <pagmo/utils/generic.hpp>
48 
49 namespace pagmo
50 {
51 
de1220(unsigned gen,std::vector<unsigned> allowed_variants,unsigned variant_adptv,double ftol,double xtol,bool memory,unsigned seed)52 de1220::de1220(unsigned gen, std::vector<unsigned> allowed_variants, unsigned variant_adptv, double ftol, double xtol,
53                bool memory, unsigned seed)
54     : m_gen(gen), m_F(), m_CR(), m_variant(), m_allowed_variants(allowed_variants), m_variant_adptv(variant_adptv),
55       m_ftol(ftol), m_xtol(xtol), m_memory(memory), m_e(seed), m_seed(seed), m_verbosity(0u)
56 {
57     for (auto variant : allowed_variants) {
58         if (variant < 1u || variant > 18u) {
59             pagmo_throw(std::invalid_argument,
60                         "All mutation variants considered must be in [1, .., 18], while a value of "
61                             + std::to_string(variant) + " was detected.");
62         }
63     }
64     if (variant_adptv < 1u || variant_adptv > 2u) {
65         pagmo_throw(std::invalid_argument, "The variant for self-adaptation must be in [1,2], while a value of "
66                                                + std::to_string(variant_adptv) + " was detected.");
67     }
68 }
69 
70 /// Algorithm evolve method
71 /**
72  * Evolves the population for a maximum number of generations, until one of
73  * tolerances set on the population flatness (x_tol, f_tol) are met.
74  *
75  * @param pop population to be evolved
76  * @return evolved population
77  * @throws std::invalid_argument if the problem is multi-objective or constrained or stochastic
78  * @throws std::invalid_argument if the population size is not at least 7
79  */
evolve(population pop) const80 population de1220::evolve(population pop) const
81 {
82     // We store some useful variables
83     const auto &prob = pop.get_problem(); // This is a const reference, so using set_seed for example will not be
84                                           // allowed
85     auto dim = prob.get_nx();             // This getter does not return a const reference but a copy
86     const auto bounds = prob.get_bounds();
87     const auto &lb = bounds.first;
88     const auto &ub = bounds.second;
89     auto NP = pop.size();
90     auto prob_f_dimension = prob.get_nf();
91     auto fevals0 = prob.get_fevals(); // disount for the already made fevals
92     unsigned count = 1u;              // regulates the screen output
93 
94     // PREAMBLE-------------------------------------------------------------------------------------------------
95     // We start by checking that the problem is suitable for this
96     // particular algorithm.
97     if (prob.get_nc() != 0u) {
98         pagmo_throw(std::invalid_argument, "Non linear constraints detected in " + prob.get_name() + " instance. "
99                                                + get_name() + " cannot deal with them");
100     }
101     if (prob_f_dimension != 1u) {
102         pagmo_throw(std::invalid_argument, "Multiple objectives detected in " + prob.get_name() + " instance. "
103                                                + get_name() + " cannot deal with them");
104     }
105     if (prob.is_stochastic()) {
106         pagmo_throw(std::invalid_argument,
107                     "The problem appears to be stochastic " + get_name() + " cannot deal with it");
108     }
109     // Get out if there is nothing to do.
110     if (m_gen == 0u) {
111         return pop;
112     }
113     if (pop.size() < 7u) {
114         pagmo_throw(std::invalid_argument, get_name() + " needs at least 7 individuals in the population, "
115                                                + std::to_string(pop.size()) + " detected");
116     }
117     // ---------------------------------------------------------------------------------------------------------
118 
119     // No throws, all valid: we clear the logs
120     m_log.clear();
121 
122     // Some vectors used during evolution are declared.
123     vector_double tmp(dim);                              // contains the mutated candidate
124     std::uniform_real_distribution<double> drng(0., 1.); // to generate a number in [0, 1)
125     std::normal_distribution<double> n_dist(0., 1.);     // to generate a normally distributed number
126     std::uniform_int_distribution<vector_double::size_type> c_idx(
127         0u, dim - 1u); // to generate a random index in the chromosome
128     std::uniform_int_distribution<vector_double::size_type> p_idx(0u, NP - 1u); // to generate a random index in pop
129     std::uniform_int_distribution<vector_double::size_type> v_idx(0u, m_allowed_variants.size()
130                                                                           - 1u); // to generate a random variant
131 
132     // We extract from pop the chromosomes and fitness associated
133     auto popold = pop.get_x();
134     auto fit = pop.get_f();
135     auto popnew = popold;
136 
137     // Initialise the global bests
138     auto best_idx = pop.best_idx();
139     vector_double::size_type worst_idx = 0u;
140     auto gbX = popnew[best_idx];
141     auto gbfit = fit[best_idx];
142     // the best decision vector of a generation
143     auto gbIter = gbX;
144     std::vector<vector_double::size_type> r(7); // indexes of 7 selected population members
145 
146     // Initialize the F and CR vectors
147     if ((m_CR.size() != NP) || (m_F.size() != NP) || (m_variant.size() != NP) || (!m_memory)) {
148         m_CR.resize(NP);
149         m_F.resize(NP);
150         m_variant.resize(NP);
151         if (m_variant_adptv == 1u) {
152             for (decltype(NP) i = 0u; i < NP; ++i) {
153                 m_CR[i] = drng(m_e);
154                 m_F[i] = drng(m_e) * 0.9 + 0.1;
155             }
156         } else if (m_variant_adptv == 2u) {
157             for (decltype(NP) i = 0u; i < NP; ++i) {
158                 m_CR[i] = n_dist(m_e) * 0.15 + 0.5;
159                 m_F[i] = n_dist(m_e) * 0.15 + 0.5;
160             }
161         }
162         for (auto &variant : m_variant) {
163             variant = m_allowed_variants[v_idx(m_e)];
164         }
165     }
166     // Initialize the global and iteration bests for F and CR
167     double gbF = m_F[0];   // initialization to the 0 ind, will soon be forgotten
168     double gbCR = m_CR[0]; // initialization to the 0 ind, will soon be forgotten
169     unsigned gbVariant = m_variant[0];
170     double gbIterF = gbF;
171     double gbIterCR = gbCR;
172     unsigned gbIterVariant;
173 
174     // We initialize the global best for F and CR as the first individual (this will soon be forgotten)
175 
176     // Main DE iterations
177     for (decltype(m_gen) gen = 1u; gen <= m_gen; ++gen) {
178         // Start of the loop through the population
179         for (decltype(NP) i = 0u; i < NP; ++i) {
180             /*-----We select at random 5 indexes from the population---------------------------------*/
181             std::vector<vector_double::size_type> idxs(NP);
182             std::iota(idxs.begin(), idxs.end(), vector_double::size_type(0u));
183             for (auto j = 0u; j < 7u; ++j) { // Durstenfeld's algorithm to select 7 indexes at random
184                 auto idx = std::uniform_int_distribution<vector_double::size_type>(0u, NP - 1u - j)(m_e);
185                 r[j] = idxs[idx];
186                 std::swap(idxs[idx], idxs[NP - 1u - j]);
187             }
188 
189             // Adapt amplification factor, crossover probability and mutation variant for DE 1220
190             double F = 0., CR = 0.;
191             unsigned VARIANT = 0u;
192             VARIANT = (drng(m_e) < 0.9) ? m_variant[i] : m_allowed_variants[v_idx(m_e)];
193             if (m_variant_adptv == 1u) {
194                 F = (drng(m_e) < 0.9) ? m_F[i] : drng(m_e) * 0.9 + 0.1;
195                 CR = (drng(m_e) < 0.9) ? m_CR[i] : drng(m_e);
196             }
197 
198             /*-------DE/best/1/exp--------------------------------------------------------------------*/
199             if (VARIANT == 1u) {
200                 if (m_variant_adptv == 2u) {
201                     F = gbIterF + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[r[2]]);
202                     CR = gbIterCR + n_dist(m_e) * 0.5 * (m_CR[r[1]] - m_CR[r[2]]);
203                 }
204                 tmp = popold[i];
205                 auto n = c_idx(m_e);
206                 auto L = 0u;
207                 do {
208                     tmp[n] = gbIter[n] + F * (popold[r[1]][n] - popold[r[2]][n]);
209                     n = (n + 1u) % dim;
210                     ++L;
211                 } while ((drng(m_e) < CR) && (L < dim));
212             }
213 
214             /*-------DE/rand/1/exp-------------------------------------------------------------------*/
215             else if (VARIANT == 2u) {
216                 if (m_variant_adptv == 2u) {
217                     F = m_F[r[0]] + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[r[2]]);
218                     CR = m_CR[r[0]] + n_dist(m_e) * 0.5 * (m_CR[r[1]] - m_CR[r[2]]);
219                 }
220                 tmp = popold[i];
221                 auto n = c_idx(m_e);
222                 decltype(dim) L = 0u;
223                 do {
224                     tmp[n] = popold[r[0]][n] + F * (popold[r[1]][n] - popold[r[2]][n]);
225                     n = (n + 1u) % dim;
226                     ++L;
227                 } while ((drng(m_e) < CR) && (L < dim));
228             }
229             /*-------DE/rand-to-best/1/exp-----------------------------------------------------------*/
230             else if (VARIANT == 3u) {
231                 if (m_variant_adptv == 2u) {
232                     F = m_F[i] + n_dist(m_e) * 0.5 * (gbIterF - m_F[i]) + n_dist(m_e) * 0.5 * (m_F[r[0]] - m_F[r[1]]);
233                     CR = m_CR[i] + n_dist(m_e) * 0.5 * (gbIterCR - m_CR[i])
234                          + n_dist(m_e) * 0.5 * (m_CR[r[0]] - m_CR[r[1]]);
235                 }
236                 tmp = popold[i];
237                 auto n = c_idx(m_e);
238                 auto L = 0u;
239                 do {
240                     tmp[n] = tmp[n] + F * (gbIter[n] - tmp[n]) + F * (popold[r[0]][n] - popold[r[1]][n]);
241                     n = (n + 1u) % dim;
242                     ++L;
243                 } while ((drng(m_e) < CR) && (L < dim));
244             }
245             /*-------DE/best/2/exp is another powerful variant worth trying--------------------------*/
246             else if (VARIANT == 4u) {
247                 if (m_variant_adptv == 2u) {
248                     F = gbIterF + n_dist(m_e) * 0.5 * (m_F[r[0]] - m_F[r[1]])
249                         + n_dist(m_e) * 0.5 * (m_F[r[2]] - m_F[r[3]]);
250                     CR = gbIterCR + n_dist(m_e) * 0.5 * (m_CR[r[0]] - m_CR[r[1]])
251                          + n_dist(m_e) * 0.5 * (m_CR[r[2]] - m_CR[r[3]]);
252                 }
253                 tmp = popold[i];
254                 auto n = c_idx(m_e);
255                 auto L = 0u;
256                 do {
257                     tmp[n]
258                         = gbIter[n] + (popold[r[0]][n] - popold[r[1]][n]) * F + (popold[r[2]][n] - popold[r[3]][n]) * F;
259                     n = (n + 1u) % dim;
260                     ++L;
261                 } while ((drng(m_e) < CR) && (L < dim));
262             }
263             /*-------DE/rand/2/exp seems to be a robust optimizer for many functions-------------------*/
264             else if (VARIANT == 5u) {
265                 if (m_variant_adptv == 2u) {
266                     F = m_F[r[4]] + n_dist(m_e) * 0.5 * (m_F[r[0]] - m_F[r[1]])
267                         + n_dist(m_e) * 0.5 * (m_F[r[2]] - m_F[r[3]]);
268                     CR = m_CR[r[4]] + n_dist(m_e) * 0.5 * (m_CR[r[0]] - m_CR[r[1]])
269                          + n_dist(m_e) * 0.5 * (m_CR[r[2]] - m_CR[r[3]]);
270                 }
271                 tmp = popold[i];
272                 auto n = c_idx(m_e);
273                 auto L = 0u;
274                 do {
275                     tmp[n] = popold[r[4]][n] + (popold[r[0]][n] - popold[r[1]][n]) * F
276                              + (popold[r[2]][n] - popold[r[3]][n]) * F;
277                     n = (n + 1u) % dim;
278                     ++L;
279                 } while ((drng(m_e) < CR) && (L < dim));
280             }
281 
282             /*=======Essentially same strategies but BINOMIAL CROSSOVER===============================*/
283             /*-------DE/best/1/bin--------------------------------------------------------------------*/
284             else if (VARIANT == 6u) {
285                 if (m_variant_adptv == 2u) {
286                     F = gbIterF + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[r[2]]);
287                     CR = gbIterCR + n_dist(m_e) * 0.5 * (m_CR[r[1]] - m_CR[r[2]]);
288                 }
289                 tmp = popold[i];
290                 auto n = c_idx(m_e);
291                 for (decltype(dim) L = 0u; L < dim; ++L) {   /* perform Dc binomial trials */
292                     if ((drng(m_e) < CR) || L + 1u == dim) { /* change at least one parameter */
293                         tmp[n] = gbIter[n] + F * (popold[r[1]][n] - popold[r[2]][n]);
294                     }
295                     n = (n + 1u) % dim;
296                 }
297             }
298             /*-------DE/rand/1/bin-------------------------------------------------------------------*/
299             else if (VARIANT == 7u) {
300                 if (m_variant_adptv == 2u) {
301                     F = m_F[r[0]] + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[r[2]]);
302                     CR = m_CR[r[0]] + n_dist(m_e) * 0.5 * (m_CR[r[1]] - m_CR[r[2]]);
303                 }
304                 tmp = popold[i];
305                 auto n = c_idx(m_e);
306                 for (decltype(dim) L = 0u; L < dim; ++L) {   /* perform Dc binomial trials */
307                     if ((drng(m_e) < CR) || L + 1u == dim) { /* change at least one parameter */
308                         tmp[n] = popold[r[0]][n] + F * (popold[r[1]][n] - popold[r[2]][n]);
309                     }
310                     n = (n + 1u) % dim;
311                 }
312             }
313             /*-------DE/rand-to-best/1/bin-----------------------------------------------------------*/
314             else if (VARIANT == 8u) {
315                 if (m_variant_adptv == 2u) {
316                     F = m_F[i] + n_dist(m_e) * 0.5 * (gbIterF - m_F[i]) + n_dist(m_e) * 0.5 * (m_F[r[0]] - m_F[r[1]]);
317                     CR = m_CR[i] + n_dist(m_e) * 0.5 * (gbIterCR - m_CR[i])
318                          + n_dist(m_e) * 0.5 * (m_CR[r[0]] - m_CR[r[1]]);
319                 }
320                 tmp = popold[i];
321                 auto n = c_idx(m_e);
322                 for (decltype(dim) L = 0u; L < dim; ++L) {   /* perform Dc binomial trials */
323                     if ((drng(m_e) < CR) || L + 1u == dim) { /* change at least one parameter */
324                         tmp[n] = tmp[n] + F * (gbIter[n] - tmp[n]) + F * (popold[r[0]][n] - popold[r[1]][n]);
325                     }
326                     n = (n + 1u) % dim;
327                 }
328             }
329             /*-------DE/best/2/bin--------------------------------------------------------------------*/
330             else if (VARIANT == 9u) {
331                 if (m_variant_adptv == 2u) {
332                     F = gbIterF + n_dist(m_e) * 0.5 * (m_F[r[0]] - m_F[r[1]])
333                         + n_dist(m_e) * 0.5 * (m_F[r[2]] - m_F[r[3]]);
334                     CR = gbIterCR + n_dist(m_e) * 0.5 * (m_CR[r[0]] - m_CR[r[1]])
335                          + n_dist(m_e) * 0.5 * (m_CR[r[2]] - m_CR[r[3]]);
336                 }
337                 tmp = popold[i];
338                 auto n = c_idx(m_e);
339                 for (decltype(dim) L = 0u; L < dim; ++L) {   /* perform Dc binomial trials */
340                     if ((drng(m_e) < CR) || L + 1u == dim) { /* change at least one parameter */
341                         tmp[n] = gbIter[n] + (popold[r[0]][n] - popold[r[1]][n]) * F
342                                  + (popold[r[2]][n] - popold[r[3]][n]) * F;
343                     }
344                     n = (n + 1u) % dim;
345                 }
346             }
347             /*-------DE/rand/2/bin--------------------------------------------------------------------*/
348             else if (VARIANT == 10u) {
349                 if (m_variant_adptv == 2u) {
350                     F = m_F[r[4]] + n_dist(m_e) * 0.5 * (m_F[r[0]] - m_F[r[1]])
351                         + n_dist(m_e) * 0.5 * (m_F[r[2]] - m_F[r[3]]);
352                     CR = m_CR[r[4]] + n_dist(m_e) * 0.5 * (m_CR[r[0]] - m_CR[r[1]])
353                          + n_dist(m_e) * 0.5 * (m_CR[r[2]] - m_CR[r[3]]);
354                 }
355                 tmp = popold[i];
356                 auto n = c_idx(m_e);
357                 for (decltype(dim) L = 0u; L < dim; ++L) {   /* perform Dc binomial trials */
358                     if ((drng(m_e) < CR) || L + 1u == dim) { /* change at least one parameter */
359                         tmp[n] = popold[r[4]][n] + (popold[r[0]][n] - popold[r[1]][n]) * F
360                                  + (popold[r[2]][n] - popold[r[3]][n]) * F;
361                     }
362                     n = (n + 1u) % dim;
363                 }
364             }
365             /*-------DE/rand/3/exp --------------------------------------------------------------------*/
366             else if (VARIANT == 11u) {
367                 if (m_variant_adptv == 2u) {
368                     F = m_F[r[0]] + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[r[2]])
369                         + n_dist(m_e) * 0.5 * (m_F[r[3]] - m_F[r[4]]) + n_dist(m_e) * 0.5 * (m_F[r[5]] - m_F[r[6]]);
370                     CR = m_CR[r[4]] + n_dist(m_e) * 0.5 * (m_CR[r[0]] + m_CR[r[1]] - m_CR[r[2]] - m_CR[r[3]]);
371                 }
372                 tmp = popold[i];
373                 auto n = c_idx(m_e);
374                 auto L = 0u;
375                 do {
376                     tmp[n] = popold[r[0]][n] + (popold[r[1]][n] - popold[r[2]][n]) * F
377                              + (popold[r[3]][n] - popold[r[4]][n]) * F + (popold[r[5]][n] - popold[r[6]][n]) * F;
378                     n = (n + 1u) % dim;
379                     ++L;
380                 } while ((drng(m_e) < CR) && (L < dim));
381             }
382             /*-------DE/rand/3/bin --------------------------------------------------------------------*/
383             else if (VARIANT == 12u) {
384                 if (m_variant_adptv == 2u) {
385                     F = m_F[r[0]] + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[r[2]])
386                         + n_dist(m_e) * 0.5 * (m_F[r[3]] - m_F[r[4]]) + n_dist(m_e) * 0.5 * (m_F[r[5]] - m_F[r[6]]);
387                     CR = m_CR[r[4]] + n_dist(m_e) * 0.5 * (m_CR[r[0]] + m_CR[r[1]] - m_CR[r[2]] - m_CR[r[3]]);
388                 }
389                 tmp = popold[i];
390                 auto n = c_idx(m_e);
391                 for (decltype(dim) L = 0u; L < dim; ++L) {   /* perform Dc binomial trials */
392                     if ((drng(m_e) < CR) || L + 1u == dim) { /* change at least one parameter */
393                         tmp[n] = popold[r[0]][n] + (popold[r[1]][n] - popold[r[2]][n]) * F
394                                  + (popold[r[3]][n] - popold[r[4]][n]) * F + (popold[r[5]][n] - popold[r[6]][n]) * F;
395                     }
396                     n = (n + 1u) % dim;
397                 }
398             }
399             /*-------DE/best/3/exp --------------------------------------------------------------------*/
400             else if (VARIANT == 13u) {
401                 if (m_variant_adptv == 2u) {
402                     F = gbIterF + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[r[2]])
403                         + n_dist(m_e) * 0.5 * (m_F[r[3]] - m_F[r[4]]) + n_dist(m_e) * 0.5 * (m_F[r[5]] - m_F[r[6]]);
404                     CR = gbIterCR + n_dist(m_e) * 0.5 * (m_CR[r[0]] + m_CR[r[1]] - m_CR[r[2]] - m_CR[r[3]]);
405                 }
406                 tmp = popold[i];
407                 auto n = c_idx(m_e);
408                 auto L = 0u;
409                 do {
410                     tmp[n] = gbIter[n] + (popold[r[1]][n] - popold[r[2]][n]) * F
411                              + (popold[r[3]][n] - popold[r[4]][n]) * F + (popold[r[5]][n] - popold[r[6]][n]) * F;
412                     n = (n + 1u) % dim;
413                     ++L;
414                 } while ((drng(m_e) < CR) && (L < dim));
415             }
416             /*-------DE/best/3/bin --------------------------------------------------------------------*/
417             else if (VARIANT == 14u) {
418                 if (m_variant_adptv == 2u) {
419                     F = gbIterF + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[r[2]])
420                         + n_dist(m_e) * 0.5 * (m_F[r[3]] - m_F[r[4]]) + n_dist(m_e) * 0.5 * (m_F[r[5]] - m_F[r[6]]);
421                     CR = gbIterCR + n_dist(m_e) * 0.5 * (m_CR[r[0]] + m_CR[r[1]] - m_CR[r[2]] - m_CR[r[3]]);
422                 }
423                 tmp = popold[i];
424                 auto n = c_idx(m_e);
425                 for (decltype(dim) L = 0u; L < dim; ++L) {   /* perform Dc binomial trials */
426                     if ((drng(m_e) < CR) || L + 1u == dim) { /* change at least one parameter */
427                         tmp[n] = gbIter[n] + (popold[r[1]][n] - popold[r[2]][n]) * F
428                                  + (popold[r[3]][n] - popold[r[4]][n]) * F + (popold[r[5]][n] - popold[r[6]][n]) * F;
429                     }
430                     n = (n + 1u) % dim;
431                 }
432             }
433             /*-------DE/rand-to-current/2/exp --------------------------------------------------------------------*/
434             else if (VARIANT == 15u) {
435                 if (m_variant_adptv == 2u) {
436                     F = m_F[r[0]] + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[i])
437                         + n_dist(m_e) * 0.5 * (m_F[r[3]] - m_F[r[4]]);
438                     CR = m_CR[r[0]] + n_dist(m_e) * 0.5 * (m_CR[r[1]] - m_CR[i])
439                          + n_dist(m_e) * 0.5 * (m_CR[r[3]] - m_CR[r[4]]);
440                 }
441                 tmp = popold[i];
442                 auto n = c_idx(m_e);
443                 auto L = 0u;
444                 do {
445                     tmp[n] = popold[r[0]][n] + (popold[r[1]][n] - popold[i][n]) * F
446                              + (popold[r[2]][n] - popold[r[3]][n]) * F;
447                     n = (n + 1u) % dim;
448                     ++L;
449                 } while ((drng(m_e) < CR) && (L < dim));
450             }
451             /*-------DE/rand-to-current/2/bin --------------------------------------------------------------------*/
452             else if (VARIANT == 16u) {
453                 if (m_variant_adptv == 2u) {
454                     F = m_F[r[0]] + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[i])
455                         + n_dist(m_e) * 0.5 * (m_F[r[3]] - m_F[r[4]]);
456                     CR = m_CR[r[0]] + n_dist(m_e) * 0.5 * (m_CR[r[1]] - m_CR[i])
457                          + n_dist(m_e) * 0.5 * (m_CR[r[3]] - m_CR[r[4]]);
458                 }
459                 tmp = popold[i];
460                 auto n = c_idx(m_e);
461                 for (decltype(dim) L = 0u; L < dim; ++L) {   /* perform Dc binomial trials */
462                     if ((drng(m_e) < CR) || L + 1u == dim) { /* change at least one parameter */
463                         tmp[n] = popold[r[0]][n] + (popold[r[1]][n] - popold[i][n]) * F
464                                  + (popold[r[2]][n] - popold[r[3]][n]) * F;
465                     }
466                     n = (n + 1u) % dim;
467                 }
468             }
469             /*-------DE/rand-to-best-and-current/2/exp
470                --------------------------------------------------------------------*/
471             else if (VARIANT == 17u) {
472                 if (m_variant_adptv == 2u) {
473                     F = m_F[r[0]] + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[i])
474                         - n_dist(m_e) * 0.5 * (m_F[r[2]] - gbIterF);
475                     CR = m_CR[r[0]] + n_dist(m_e) * 0.5 * (m_CR[r[1]] - m_CR[i])
476                          - n_dist(m_e) * 0.5 * (m_CR[r[3]] - gbIterCR);
477                 }
478                 tmp = popold[i];
479                 auto n = c_idx(m_e);
480                 auto L = 0u;
481                 do {
482                     tmp[n] = popold[r[0]][n] + (popold[r[1]][n] - popold[i][n]) * F - (popold[r[2]][n] - gbIter[n]) * F;
483                     n = (n + 1u) % dim;
484                     ++L;
485                 } while ((drng(m_e) < CR) && (L < dim));
486             }
487             /*-------DE/rand-to-best-and-current/2/bin
488                --------------------------------------------------------------------*/
489             else if (VARIANT == 18u) {
490                 if (m_variant_adptv == 2u) {
491                     F = m_F[r[0]] + n_dist(m_e) * 0.5 * (m_F[r[1]] - m_F[i])
492                         - n_dist(m_e) * 0.5 * (m_F[r[2]] - gbIterF);
493                     CR = m_CR[r[0]] + n_dist(m_e) * 0.5 * (m_CR[r[1]] - m_CR[i])
494                          - n_dist(m_e) * 0.5 * (m_CR[r[3]] - gbIterCR);
495                 }
496                 tmp = popold[i];
497                 auto n = c_idx(m_e);
498                 for (decltype(dim) L = 0u; L < dim; ++L) {   /* perform Dc binomial trials */
499                     if ((drng(m_e) < CR) || L + 1u == dim) { /* change at least one parameter */
500                         tmp[n] = popold[r[0]][n] + (popold[r[1]][n] - popold[i][n]) * F
501                                  - (popold[r[2]][n] - gbIter[n]) * F;
502                     }
503                     n = (n + 1u) % dim;
504                 }
505             }
506 
507             /*==Trial mutation now in tmp. force feasibility and see how good this choice really was.==*/
508             // a) feasibility
509             for (decltype(dim) j = 0u; j < dim; ++j) {
510                 if ((tmp[j] < lb[j]) || (tmp[j] > ub[j])) {
511                     tmp[j] = uniform_real_from_range(lb[j], ub[j], m_e);
512                 }
513             }
514             // b) how good?
515             auto newfitness = prob.fitness(tmp); /* Evaluates tmp[] */
516             if (newfitness[0] <= fit[i][0]) {    /* improved objective function value ? */
517                 fit[i] = newfitness;
518                 popnew[i] = tmp;
519                 // updates the individual in pop (avoiding to recompute the objective function)
520                 pop.set_xf(i, popnew[i], newfitness);
521                 // Update the adapted parameters
522                 m_CR[i] = CR;
523                 m_F[i] = F;
524                 m_variant[i] = VARIANT;
525 
526                 if (newfitness[0] <= gbfit[0]) {
527                     /* if so...*/
528                     gbfit = newfitness; /* reset gbfit to new low...*/
529                     gbX = popnew[i];
530                     gbF = F;   /* these were forgotten in PaGMOlegacy */
531                     gbCR = CR; /* these were forgotten in PaGMOlegacy */
532                     gbVariant = VARIANT;
533                 }
534             } else {
535                 popnew[i] = popold[i];
536             }
537         } // End of one generation
538         /* Save best population member of current iteration */
539         gbIter = gbX;
540         gbIterF = gbF;
541         gbIterCR = gbCR;
542         gbIterVariant = gbVariant; // the gbIterVariant is not really needed and is only kept for consistency
543         /* swap population arrays. New generation becomes old one */
544         std::swap(popold, popnew);
545 
546         // Check the exit conditions
547         double dx = 0., df = 0.;
548 
549         best_idx = pop.best_idx();
550         worst_idx = pop.worst_idx();
551         for (decltype(dim) i = 0u; i < dim; ++i) {
552             dx += std::abs(pop.get_x()[worst_idx][i] - pop.get_x()[best_idx][i]);
553         }
554         if (dx < m_xtol) {
555             if (m_verbosity > 0u) {
556                 std::cout << "Exit condition -- xtol < " << m_xtol << std::endl;
557             }
558             return pop;
559         }
560 
561         df = std::abs(pop.get_f()[worst_idx][0] - pop.get_f()[best_idx][0]);
562         if (df < m_ftol) {
563             if (m_verbosity > 0u) {
564                 std::cout << "Exit condition -- ftol < " << m_ftol << std::endl;
565             }
566             return pop;
567         }
568 
569         // Logs and prints (verbosity modes > 1: a line is added every m_verbosity generations)
570         if (m_verbosity > 0u) {
571             // Every m_verbosity generations print a log line
572             if (gen % m_verbosity == 1u || m_verbosity == 1u) {
573                 best_idx = pop.best_idx();
574                 worst_idx = pop.worst_idx();
575                 dx = 0.;
576                 // The population flattness in chromosome
577                 for (decltype(dim) i = 0u; i < dim; ++i) {
578                     dx += std::abs(pop.get_x()[worst_idx][i] - pop.get_x()[best_idx][i]);
579                 }
580                 // The population flattness in fitness
581                 df = std::abs(pop.get_f()[worst_idx][0] - pop.get_f()[best_idx][0]);
582                 // Every 50 lines print the column names
583                 if (count % 50u == 1u) {
584                     print("\n", std::setw(7), "Gen:", std::setw(15), "Fevals:", std::setw(15), "Best:", std::setw(15),
585                           "F:", std::setw(15), "CR:", std::setw(15), "Variant:", std::setw(15), "dx:", std::setw(15),
586                           std::setw(15), "df:", '\n');
587                 }
588                 print(std::setw(7), gen, std::setw(15), prob.get_fevals() - fevals0, std::setw(15),
589                       pop.get_f()[best_idx][0], std::setw(15), gbIterF, std::setw(15), gbIterCR, std::setw(15),
590                       gbIterVariant, std::setw(15), dx, std::setw(15), df, '\n');
591                 ++count;
592                 // Logs
593                 m_log.emplace_back(gen, prob.get_fevals() - fevals0, pop.get_f()[best_idx][0], gbIterF, gbIterCR,
594                                    gbIterVariant, dx, df);
595             }
596         }
597     } // end main DE iterations
598     if (m_verbosity) {
599         std::cout << "Exit condition -- generations = " << m_gen << std::endl;
600     }
601     return pop;
602 }
603 
604 /// Sets the seed
605 /**
606  * @param seed the seed controlling the algorithm stochastic behaviour
607  */
set_seed(unsigned seed)608 void de1220::set_seed(unsigned seed)
609 {
610     m_e.seed(seed);
611     m_seed = seed;
612 }
613 
614 /// Extra info
615 /**
616  * One of the optional methods of any user-defined algorithm (UDA).
617  *
618  * @return a string containing extra info on the algorithm
619  */
get_extra_info() const620 std::string de1220::get_extra_info() const
621 {
622     std::ostringstream ss;
623     stream(ss, "\tGenerations: ", m_gen);
624     stream(ss, "\n\tAllowed variants: ", m_allowed_variants);
625     stream(ss, "\n\tSelf adaptation variant: ", m_variant_adptv);
626     stream(ss, "\n\tStopping xtol: ", m_xtol);
627     stream(ss, "\n\tStopping ftol: ", m_ftol);
628     stream(ss, "\n\tMemory: ", m_memory);
629     stream(ss, "\n\tVerbosity: ", m_verbosity);
630     stream(ss, "\n\tSeed: ", m_seed);
631     return ss.str();
632 }
633 
634 // Object serialization
635 template <typename Archive>
serialize(Archive & ar,unsigned)636 void de1220::serialize(Archive &ar, unsigned)
637 {
638     detail::archive(ar, m_gen, m_F, m_CR, m_allowed_variants, m_variant_adptv, m_ftol, m_xtol, m_memory, m_e, m_seed,
639                     m_verbosity, m_log);
640 }
641 
642 } // namespace pagmo
643 
644 PAGMO_S11N_ALGORITHM_IMPLEMENT(pagmo::de1220)
645