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 <algorithm>
30 #include <cmath>
31 #include <iomanip>
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/not_population_based.hpp>
42 #include <pagmo/algorithms/simulated_annealing.hpp>
43 #include <pagmo/exceptions.hpp>
44 #include <pagmo/io.hpp>
45 #include <pagmo/population.hpp>
46 #include <pagmo/s11n.hpp>
47 #include <pagmo/types.hpp>
48 #include <pagmo/utils/generic.hpp>
49
50 namespace pagmo
51 {
52
simulated_annealing(double Ts,double Tf,unsigned n_T_adj,unsigned n_range_adj,unsigned bin_size,double start_range,unsigned seed)53 simulated_annealing::simulated_annealing(double Ts, double Tf, unsigned n_T_adj, unsigned n_range_adj,
54 unsigned bin_size, double start_range, unsigned seed)
55 : m_Ts(Ts), m_Tf(Tf), m_n_T_adj(n_T_adj), m_n_range_adj(n_range_adj), m_bin_size(bin_size),
56 m_start_range(start_range), m_e(seed), m_seed(seed), m_verbosity(0u), m_log()
57 {
58 if (Ts <= 0. || !std::isfinite(Ts)) {
59 pagmo_throw(std::invalid_argument, "The starting temperature must be finite and positive, while a value of "
60 + std::to_string(Ts) + " was detected.");
61 }
62 if (Tf <= 0. || !std::isfinite(Tf)) {
63 pagmo_throw(std::invalid_argument, "The final temperature must be finite and positive, while a value of "
64 + std::to_string(Tf) + " was detected.");
65 }
66 if (Tf > Ts) {
67 pagmo_throw(std::invalid_argument,
68 "The final temperature must be smaller than the initial temperature, while a value of "
69 + std::to_string(Tf) + " >= " + std::to_string(Ts) + " was detected.");
70 }
71 if (start_range <= 0. || start_range > 1. || !std::isfinite(start_range)) {
72 pagmo_throw(std::invalid_argument, "The start range must be in the (0, 1] range, while a value of "
73 + std::to_string(start_range) + " was detected.");
74 }
75 if (n_T_adj == 0u) {
76 pagmo_throw(std::invalid_argument,
77 "The number of temperature adjustments must be strictly positive, while a value of "
78 + std::to_string(n_T_adj) + " was detected.");
79 }
80 if (n_range_adj == 0u) {
81 pagmo_throw(std::invalid_argument,
82 "The number of range adjustments must be strictly positive, while a value of "
83 + std::to_string(n_range_adj) + " was detected.");
84 }
85 }
86
87 /// Algorithm evolve method
88 /**
89 * @param pop population to be evolved
90 * @return evolved population
91 * @throws std::invalid_argument if the problem is multi-objective, constrained or stochastic
92 * @throws std::invalid_argument if the population size is < 1u
93 */
evolve(population pop) const94 population simulated_annealing::evolve(population pop) const
95 {
96 // We store some useful properties
97 const auto &prob = pop.get_problem(); // This is a const reference, so using set_seed for example will not be
98 // allowed
99 auto dim = prob.get_nx(); // not const as used type for counters
100 const auto bounds = prob.get_bounds();
101 const auto &lb = bounds.first;
102 const auto &ub = bounds.second;
103 auto fevals0 = prob.get_fevals(); // disount for the already made fevals
104 unsigned count = 1u; // regulates the screen output
105
106 // PREAMBLE-------------------------------------------------------------------------------------------------
107 // We start by checking that the problem is suitable for this particular algorithm.
108 if (prob.get_nc() != 0u) {
109 pagmo_throw(std::invalid_argument, "Non linear constraints detected in " + prob.get_name() + " instance. "
110 + get_name() + " cannot deal with them");
111 }
112 if (prob.get_nf() != 1u) {
113 pagmo_throw(std::invalid_argument, "Multiple objectives detected in " + prob.get_name() + " instance. "
114 + get_name() + " cannot deal with them");
115 }
116 if (prob.is_stochastic()) {
117 pagmo_throw(std::invalid_argument,
118 "The problem appears to be stochastic " + get_name() + " cannot deal with it");
119 }
120 if (!pop.size()) {
121 pagmo_throw(std::invalid_argument, get_name() + " does not work on an empty population");
122 }
123 // ---------------------------------------------------------------------------------------------------------
124
125 // No throws, all valid: we clear the logs
126 m_log.clear();
127
128 std::uniform_real_distribution<double> drng(0., 1.); // to generate a number in [0, 1)
129
130 // We init the starting point
131 auto sel_xf = select_individual(pop);
132 vector_double x0(std::move(sel_xf.first)), fit0(std::move(sel_xf.second));
133 // Determines the coefficient to decrease the temperature
134 const double Tcoeff = std::pow(m_Tf / m_Ts, 1.0 / static_cast<double>(m_n_T_adj));
135 // Stores the current and new points
136 auto xNEW = x0;
137 auto xOLD = x0;
138 auto best_x = x0;
139 auto fNEW = fit0;
140 auto fOLD = fit0;
141 auto best_f = fit0;
142
143 // Stores the adaptive ranges for each component
144 vector_double step(dim, m_start_range);
145
146 // Stores the number of accepted points for each component
147 std::vector<int> acp(dim, 0u);
148 double ratio = 0., currentT = m_Ts, probab = 0.;
149
150 // Main SA loops
151 for (decltype(m_n_T_adj) jter = 0u; jter < m_n_T_adj; ++jter) {
152 for (decltype(m_n_range_adj) mter = 0u; mter < m_n_range_adj; ++mter) {
153 // 1 - Annealing
154 for (decltype(m_bin_size) kter = 0u; kter < m_bin_size; ++kter) {
155 auto nter = std::uniform_int_distribution<vector_double::size_type>(0u, dim - 1u)(m_e);
156 for (decltype(dim) numb = 0u; numb < dim; ++numb) {
157 nter = (nter + 1u) % dim;
158 // We modify the current point by mutating its nter component within the adaptive step
159 auto width = step[nter] * (ub[nter] - lb[nter]);
160 xNEW[nter] = uniform_real_from_range(std::max(xOLD[nter] - width, lb[nter]),
161 std::min(xOLD[nter] + width, ub[nter]), m_e);
162 // And we valuate the objective function for the new point
163 fNEW = prob.fitness(xNEW);
164 // We decide wether to accept or discard the point
165 if (fNEW[0] <= fOLD[0]) {
166 // accept
167 xOLD[nter] = xNEW[nter];
168 fOLD = fNEW;
169 acp[nter]++; // Increase the number of accepted values
170 // We update the best
171 if (fNEW[0] <= best_f[0]) {
172 best_f = fNEW;
173 best_x = xNEW;
174 }
175 } else {
176 // test it with Boltzmann to decide the acceptance
177 probab = std::exp(-std::abs(fOLD[0] - fNEW[0]) / currentT);
178 // we compare prob with a random probability.
179 if (probab > drng(m_e)) {
180 xOLD[nter] = xNEW[nter];
181 fOLD = fNEW;
182 acp[nter]++; // Increase the number of accepted values
183 } else {
184 xNEW[nter] = xOLD[nter];
185 }
186 }
187 // 2 - We log to screen
188 if (m_verbosity > 0u) {
189 // Prints a log line every m_verbosity fitness evaluations
190 auto fevals_count = prob.get_fevals() - fevals0;
191 if (fevals_count >= (count - 1u) * m_verbosity) {
192 // 1 - Every 50 lines print the column names
193 if (count % 50u == 1u) {
194 print("\n", std::setw(7), "Fevals:", std::setw(15), "Best:", std::setw(15),
195 "Current:", std::setw(15), "Mean range:", std::setw(15), "Temperature:", '\n');
196 }
197 auto avg_range
198 = std::accumulate(step.begin(), step.end(), 0.) / static_cast<double>(step.size());
199 // 2 - Print
200 print(std::setw(7), fevals_count, std::setw(15), best_f[0], std::setw(15), fOLD[0],
201 std::setw(15), avg_range, std::setw(15), currentT);
202 ++count;
203 std::cout << std::endl; // we flush here as we want the user to read in real time ...
204 // Logs
205 m_log.emplace_back(fevals_count, best_f[0], fOLD[0], avg_range, currentT);
206 }
207 }
208 } // end for(nter = 0; ...
209 } // end for(kter = 0; ...
210 // adjust the step (adaptively)
211 for (decltype(dim) iter = 0u; iter < dim; ++iter) {
212 ratio = static_cast<double>(acp[iter]) / static_cast<double>(m_bin_size);
213 acp[iter] = 0u; // reset the counter
214 if (ratio > .6) {
215 // too many acceptances, increase the step by a factor 3 maximum
216 step[iter] = step[iter] * (1. + 2. * (ratio - .6) / .4);
217 } else {
218 if (ratio < .4) {
219 // too few acceptance, decrease the step by a factor 3 maximum
220 step[iter] = step[iter] / (1. + 2. * ((.4 - ratio) / .4));
221 };
222 };
223 // And if it becomes too large, reset it to its initial value
224 if (step[iter] > m_start_range) step[iter] = m_start_range;
225 }
226 }
227 // Cooling schedule
228 currentT *= Tcoeff;
229 }
230 // We update the decision vector in pop, but only if things have improved
231 if (best_f[0] <= fit0[0]) {
232 replace_individual(pop, best_x, best_f);
233 }
234 return pop;
235 }
236
237 /// Sets the seed
238 /**
239 * @param seed the seed controlling the algorithm stochastic behaviour
240 */
set_seed(unsigned seed)241 void simulated_annealing::set_seed(unsigned seed)
242 {
243 m_e.seed(seed);
244 m_seed = seed;
245 }
246
247 /// Extra info
248 /**
249 * One of the optional methods of any user-defined algorithm (UDA).
250 *
251 * @return a string containing extra info on the algorithm
252 */
get_extra_info() const253 std::string simulated_annealing::get_extra_info() const
254 {
255 std::ostringstream ss;
256 stream(ss, "\tStarting temperature: ", m_Ts);
257 stream(ss, "\n\tFinal temperature: ", m_Tf);
258 stream(ss, "\n\tNumber of temperature adjustments: ", m_n_T_adj);
259 stream(ss, "\n\tNumber of range adjustments: ", m_n_range_adj);
260 stream(ss, "\n\tBin size: ", m_bin_size);
261 stream(ss, "\n\tStarting range: ", m_start_range);
262 stream(ss, "\n\tSeed: ", m_seed);
263 stream(ss, "\n\tVerbosity: ", m_verbosity);
264 return ss.str();
265 }
266
267 // Object serialization
268 template <typename Archive>
serialize(Archive & ar,unsigned)269 void simulated_annealing::serialize(Archive &ar, unsigned)
270 {
271 detail::archive(ar, boost::serialization::base_object<not_population_based>(*this), m_Ts, m_Tf, m_n_T_adj,
272 m_n_range_adj, m_bin_size, m_start_range, m_e, m_seed, m_verbosity, m_log);
273 }
274
275 } // namespace pagmo
276
277 PAGMO_S11N_ALGORITHM_IMPLEMENT(pagmo::simulated_annealing)
278