1//! @file optimal.cpp
2//! @author ryftchen
3//! @brief The definitions (optimal) in the algorithm module.
4//! @version 0.1.0
5//! @copyright Copyright (c) 2022-2026 ryftchen. All rights reserved.
6
7#include "optimal.hpp"
8
9#include <algorithm>
10
11namespace algorithm::optimal
12{
13//! @brief Function version number.
14//! @return version number (major.minor.patch)
15const char* version() noexcept
16{
17 static const char* const ver = "0.1.0";
18 return ver;
19}
20
21Result Gradient::operator()(const double left, const double right, const double eps)
22{
23 if (!func)
24 {
25 return std::nullopt;
26 }
27
28 const auto climbing = createClimbers(left, right);
29 double xBest = *climbing.cbegin();
30 double yBest = func(xBest);
31
32 for (std::uint32_t iteration = 0; auto x : climbing)
33 {
34 iteration = 0;
35 double learningRate = initialLR;
36 double gradient = calculateFirstDerivative(x, eps);
37 double dx = learningRate * gradient;
38 while ((std::fabs(x: dx) > eps) && ((x - dx) >= left) && ((x - dx) <= right))
39 {
40 x -= dx;
41 ++iteration;
42 learningRate = initialLR / (1.0 + decay * iteration);
43 gradient = calculateFirstDerivative(x, eps);
44 dx = learningRate * gradient;
45 }
46
47 if (const double y = func(x); y < yBest)
48 {
49 xBest = x;
50 yBest = y;
51 }
52 }
53 return std::make_optional(t: std::make_tuple(args&: yBest, args&: xBest));
54}
55
56std::unordered_multiset<double> Gradient::createClimbers(const double left, const double right) const
57{
58 std::unordered_multiset<double> climbing{};
59 climbing.reserve(n: loopTime);
60 std::generate_n(
61 first: std::inserter(x&: climbing, i: climbing.end()),
62 n: loopTime,
63 gen: [engine = std::mt19937_64(std::random_device{}()),
64 candidate = std::uniform_real_distribution<double>(left, right)]() mutable { return candidate(engine); });
65 return climbing;
66}
67
68double Gradient::calculateFirstDerivative(const double x, const double eps) const
69{
70 const double differential = eps / 2.0;
71 return (func(x + differential) - func(x - differential)) / eps;
72}
73
74Result Tabu::operator()(const double left, const double right, const double eps)
75{
76 if (!func)
77 {
78 return std::nullopt;
79 }
80
81 std::mt19937_64 engine(std::random_device{}());
82 double solution = std::uniform_real_distribution<double>{left, right}(engine);
83 double xBest = solution;
84 double yBest = func(xBest);
85 double stepLen = initialStep;
86
87 std::vector<double> neighborhood{};
88 std::vector<double> tabuList{};
89 neighborhood.reserve(n: neighborSize);
90 tabuList.reserve(n: tabuTenure);
91 tabuList.emplace_back(args&: solution);
92 for (std::uint32_t i = 0; i < maxIterations; ++i)
93 {
94 updateNeighborhood(neighborhood, solution, stepLen, left, right);
95 if (const auto [yBestNbr, xBestNbr] = neighborhoodSearch(neighborhood, solution, tabuList); yBestNbr < yBest)
96 {
97 xBest = xBestNbr;
98 yBest = yBestNbr;
99 }
100 neighborhood.clear();
101
102 solution = xBest;
103 if (tabuList.size() == tabuTenure)
104 {
105 tabuList.erase(position: tabuList.cbegin());
106 }
107 tabuList.emplace_back(args&: solution);
108 stepLen *= expDecay;
109 if (stepLen < eps)
110 {
111 break;
112 }
113 }
114 return std::make_optional(t: std::make_tuple(args&: yBest, args&: xBest));
115}
116
117void Tabu::updateNeighborhood(
118 std::vector<double>& neighborhood,
119 const double solution,
120 const double stepLen,
121 const double left,
122 const double right) const
123{
124 for (std::uint32_t i = 0; i < (neighborSize / 2); ++i)
125 {
126 const double neighbor1 = std::min(a: solution + (i * stepLen), b: right);
127 neighborhood.emplace_back(args: neighbor1);
128
129 const double neighbor2 = std::max(a: solution - (i * stepLen), b: left);
130 neighborhood.emplace_back(args: neighbor2);
131 }
132}
133
134std::tuple<double, double> Tabu::neighborhoodSearch(
135 const std::vector<double>& neighborhood, const double solution, const std::vector<double>& tabuList)
136{
137 double xBestNbr = solution;
138 double yBestNbr = func(xBestNbr);
139 for (const auto neighbor : neighborhood)
140 {
141 if (std::ranges::find(tabuList, neighbor) == tabuList.cend())
142 {
143 if (const double fitness = func(neighbor); fitness < yBestNbr)
144 {
145 xBestNbr = neighbor;
146 yBestNbr = fitness;
147 }
148 }
149 }
150 return std::make_tuple(args&: yBestNbr, args&: xBestNbr);
151}
152
153Result Annealing::operator()(const double left, const double right, const double eps)
154{
155 if (!func)
156 {
157 return std::nullopt;
158 }
159
160 std::mt19937_64 engine(std::random_device{}());
161 std::uniform_real_distribution<double> perturbation(left, right);
162 std::uniform_real_distribution<double> pr(0.0, 1.0);
163 double temperature = initialT;
164 double x = std::round(x: perturbation(engine) / eps) * eps;
165 double y = func(x);
166 double xBest = x;
167 double yBest = y;
168
169 while (temperature > minimalT)
170 {
171 x = xBest;
172 y = yBest;
173 for (std::uint32_t i = 0; i < markovChainLen; ++i)
174 {
175 double xNbr = cauchyLikeDistribution(prev: x, min: left, max: right, temp: temperature, xi: pr(engine));
176 if ((xNbr < left) || (xNbr > right))
177 {
178 xNbr = std::round(x: perturbation(engine) / eps) * eps;
179 }
180 if (const double yNbr = func(xNbr); metropolisAcceptanceCriterion(deltaE: yNbr - y, temp: temperature, xi: pr(engine)))
181 {
182 x = xNbr;
183 y = yNbr;
184 if (y < yBest)
185 {
186 xBest = x;
187 yBest = y;
188 }
189 }
190 }
191 temperature *= coolingRate;
192 }
193 return std::make_optional(t: std::make_tuple(args&: yBest, args&: xBest));
194}
195
196double Annealing::cauchyLikeDistribution(
197 const double prev, const double min, const double max, const double temp, const double xi)
198{
199 const double sgn = static_cast<int>((xi - 0.5) > 0.0) - static_cast<int>((xi - 0.5) < 0.0);
200 return prev + ((temp * sgn * (1.0 + 1.0 / (std::pow(x: temp, y: (2.0 * xi) - 1.0) - 1.0))) * (max - min));
201}
202
203bool Annealing::metropolisAcceptanceCriterion(const double deltaE, const double temp, const double xi)
204{
205 return (deltaE < 0.0) || (std::exp(x: -deltaE / temp) > xi);
206}
207
208Result Particle::operator()(const double left, const double right, const double eps)
209{
210 if (!func)
211 {
212 return std::nullopt;
213 }
214
215 auto swarm = swarmInit(left, right);
216 const auto initialBest =
217 std::ranges::min_element(swarm, [](const auto& min1, const auto& min2) { return min1.fitness < min2.fitness; });
218 double gloBest = initialBest->x;
219 double gloBestFitness = initialBest->fitness;
220
221 for (std::uint32_t i = 0; i < maxIterations; ++i)
222 {
223 updateParticles(swarm, iteration: i, gloBest, left, right, eps);
224 updateBests(swarm, gloBest, gloBestFitness);
225 }
226 return std::make_optional(t: std::make_tuple(args&: gloBestFitness, args&: gloBest));
227}
228
229Particle::Swarm Particle::swarmInit(const double left, const double right)
230{
231 Swarm swarm{};
232 swarm.reserve(n: swarmSize);
233 std::generate_n(
234 first: std::back_inserter(x&: swarm),
235 n: swarmSize,
236 gen: [this,
237 candidate = std::uniform_real_distribution<double>(left, right),
238 velocity = std::uniform_real_distribution<double>(vMin, vMax)]() mutable
239 {
240 const double x = candidate(engine);
241 return Individual{.x: x, .v: velocity(engine), .persBest: x, .fitness: func(x), .persBestFitness: func(x)};
242 });
243 return swarm;
244}
245
246void Particle::updateParticles(
247 Swarm& swarm,
248 const std::uint32_t iteration,
249 const double gloBest,
250 const double left,
251 const double right,
252 const double eps)
253{
254 for (const double w = nonlinearDecreasingWeight(iteration); auto& ind : swarm)
255 {
256 const double rand1 = std::round(x: perturbation(engine) / eps) * eps;
257 const double rand2 = std::round(x: perturbation(engine) / eps) * eps;
258 ind.v = w * ind.v + c1 * rand1 * (ind.persBest - ind.x) + c2 * rand2 * (gloBest - ind.x);
259 ind.v = std::clamp(val: ind.v, lo: vMin, hi: vMax);
260
261 ind.x += ind.v;
262 ind.x = std::clamp(val: ind.x, lo: left, hi: right);
263
264 ind.fitness = func(ind.x);
265 }
266}
267
268double Particle::nonlinearDecreasingWeight(const std::uint32_t iteration) const
269{
270 return wBegin - ((wBegin - wEnd) * std::pow(x: (1.0 + iteration) / maxIterations, y: 2.0));
271}
272
273void Particle::updateBests(Swarm& swarm, double& gloBest, double& gloBestFitness)
274{
275 for (auto& ind : swarm)
276 {
277 if (ind.fitness < ind.persBestFitness)
278 {
279 ind.persBest = ind.x;
280 ind.persBestFitness = ind.fitness;
281 }
282 if (ind.fitness < gloBestFitness)
283 {
284 gloBest = ind.x;
285 gloBestFitness = ind.fitness;
286 }
287 }
288}
289
290Result Ant::operator()(const double left, const double right, const double eps)
291{
292 if (!func)
293 {
294 return std::nullopt;
295 }
296
297 auto colony = colonyInit(left, right);
298 double xBest = colony.front().position;
299 double yBest = func(xBest);
300 double stepLen = initialStep;
301
302 for (std::uint32_t i = 0; i < maxIterations; ++i)
303 {
304 stateTransition(colony, eps);
305 pathConstruction(colony, stepLen, left, right);
306 updatePheromones(colony);
307
308 const auto currBest = std::ranges::min_element(
309 colony, [](const auto& min1, const auto& min2) { return min1.pheromone < min2.pheromone; });
310 if (const double x = currBest->position, y = func(x); y < yBest)
311 {
312 xBest = x;
313 yBest = y;
314 }
315 stepLen = initialStep / (1.0 + i);
316 if (stepLen < eps)
317 {
318 break;
319 }
320 }
321 return std::make_optional(t: std::make_tuple(args&: yBest, args&: xBest));
322}
323
324Ant::Colony Ant::colonyInit(const double left, const double right)
325{
326 Colony colony{};
327 colony.reserve(n: numOfAnts);
328 std::generate_n(
329 first: std::back_inserter(x&: colony),
330 n: numOfAnts,
331 gen: [this, candidate = std::uniform_real_distribution<double>(left, right)]() mutable
332 {
333 const double x = candidate(engine);
334 return State{.position: x, .pheromone: func(x), .transPr: 0.0};
335 });
336 return colony;
337}
338
339void Ant::stateTransition(Colony& colony, const double eps)
340{
341 const double minTau =
342 std::ranges::min_element(
343 colony, [](const auto& min1, const auto& min2) { return min1.pheromone < min2.pheromone; })
344 ->pheromone;
345 const double maxTau =
346 std::ranges::max_element(
347 colony, [](const auto& max1, const auto& max2) { return max1.pheromone < max2.pheromone; })
348 ->pheromone;
349 std::ranges::for_each(
350 colony,
351 [minTau, maxTau, eps](auto& state) { state.transPr = (maxTau - state.pheromone) / (maxTau - minTau + eps); });
352}
353
354void Ant::pathConstruction(Colony& colony, const double stepLen, const double left, const double right)
355{
356 for (auto& state : colony)
357 {
358 double pos = state.position
359 + ((state.transPr < p0) ? stepLen * localCoeff(engine) : (right - left) * globalCoeff(engine));
360 pos = std::clamp(val: pos, lo: left, hi: right);
361 if (func(pos) < func(state.position))
362 {
363 state.position = pos;
364 }
365 }
366}
367
368void Ant::updatePheromones(Colony& colony)
369{
370 std::ranges::for_each(
371 colony, [this](auto& state) { state.pheromone = (1.0 - rho) * state.pheromone + func(state.position); });
372}
373
374auto Genetic::getBestIndividual(const Population& pop)
375{
376 std::vector<double> fitness{};
377 fitness.reserve(n: pop.size());
378 std::ranges::transform(pop, std::back_inserter(x&: fitness), [this](const auto& ind) { return calculateFitness(chr: ind); });
379 return std::next(x: pop.cbegin(), n: std::ranges::distance(fitness.cbegin(), std::ranges::max_element(fitness)));
380}
381
382Result Genetic::operator()(const double left, const double right, const double eps)
383{
384 if (!func || !updateSpecies(left, right, eps))
385 {
386 return std::nullopt;
387 }
388
389 auto pop = populationInit();
390 for (std::uint32_t i = 0; i < numOfGens; ++i)
391 {
392 auto elite = extractElite(pop);
393 select(pop);
394 crossover(pop);
395 mutate(pop);
396 pop.emplace_back(args: std::move(elite));
397 }
398 const double xBest = geneticDecode(chr: *getBestIndividual(pop));
399 return std::make_optional(t: std::make_tuple(args: func(xBest), args: xBest));
400}
401
402bool Genetic::updateSpecies(const double left, const double right, const double eps)
403{
404 const double stateNum = ((right - left) / eps) + 1.0;
405 const auto estimatedLen = static_cast<std::uint32_t>(std::ceil(x: std::log2(x: stateNum)));
406 if (estimatedLen < minChrLen)
407 {
408 return false;
409 }
410
411 chromosomeLen = estimatedLen;
412 property = Property{.lower: left, .upper: right, .prec: eps};
413 return true;
414}
415
416double Genetic::geneticDecode(const Chromosome& chr) const
417{
418 std::uint32_t currDecoded = 0;
419 std::ranges::for_each(chr, [&currDecoded](const auto bit) { currDecoded = (currDecoded << 1) | bit; });
420 const auto maxDecoded = static_cast<double>((1ULL << chromosomeLen) - 1ULL);
421 return property.lower + ((property.upper - property.lower) * currDecoded / maxDecoded);
422}
423
424Genetic::Population Genetic::populationInit()
425{
426 Population pop{};
427 pop.reserve(n: popSize);
428 std::generate_n(
429 first: std::back_inserter(x&: pop),
430 n: popSize,
431 gen: [this]()
432 {
433 Chromosome chr{};
434 chr.reserve(n: chromosomeLen);
435 std::generate_n(
436 first: std::back_inserter(x&: chr),
437 n: chromosomeLen,
438 gen: [this, bit = std::uniform_int_distribution<std::uint8_t>(0, 1)]() mutable { return bit(engine); });
439 return chr;
440 });
441 return pop;
442}
443
444void Genetic::geneticCross(Chromosome& chr1, Chromosome& chr2)
445{
446 std::uniform_int_distribution<std::uint32_t> randomPos(1, chromosomeLen - 2);
447 std::uint32_t point1 = randomPos(engine);
448 std::uint32_t point2 = randomPos(engine);
449 if (point1 > point2)
450 {
451 std::swap(a&: point1, b&: point2);
452 }
453 std::swap_ranges(first1: chr1.begin() + point1, last1: chr1.begin() + point2 + 1, first2: chr2.begin() + point1);
454}
455
456void Genetic::crossover(Population& pop)
457{
458 std::vector<std::reference_wrapper<Chromosome>> selector(pop.begin(), pop.end());
459 std::ranges::shuffle(selector, engine);
460 for (auto chrIter = selector.cbegin(); (chrIter != selector.cend()) && (std::next(x: chrIter, n: 1) != selector.cend());
461 std::advance(i&: chrIter, n: 2))
462 {
463 if (probability(engine) < crossPr)
464 {
465 auto& parent1 = chrIter->get();
466 auto& parent2 = std::next(x: chrIter, n: 1)->get();
467 geneticCross(chr1&: parent1, chr2&: parent2);
468 }
469 }
470}
471
472void Genetic::geneticMutation(Chromosome& chr)
473{
474 std::ranges::for_each(
475 chr,
476 [this](auto& bit)
477 {
478 if (probability(engine) < mutatePr)
479 {
480 bit = !bit;
481 }
482 });
483}
484
485void Genetic::mutate(Population& pop)
486{
487 std::ranges::for_each(pop, [this](auto& ind) { geneticMutation(chr&: ind); });
488}
489
490double Genetic::calculateFitness(const Chromosome& chr)
491{
492 return -func(geneticDecode(chr));
493}
494
495std::optional<std::pair<double, double>> Genetic::goldbergLinearScaling(
496 const std::vector<double>& fitness, const double eps)
497{
498 const double fitMax = *std::ranges::max_element(fitness);
499 const double fitAvg = std::reduce(first: fitness.cbegin(), last: fitness.cend(), init: 0.0) / static_cast<double>(fitness.size());
500 if (std::fabs(x: fitMax - fitAvg) < eps)
501 {
502 return std::nullopt;
503 }
504 const double alpha = (cMult - 1.0) * fitAvg / (fitMax - fitAvg);
505 const double beta = fitAvg * (1.0 - alpha);
506 return std::make_optional(t: std::pair<double, double>(alpha, beta));
507}
508
509auto Genetic::rouletteWheelSelection(const Population& pop, const std::vector<double>& cumFitness)
510{
511 const auto cumIter = std::ranges::lower_bound(cumFitness, probability(engine));
512 return (cumIter != cumFitness.cend()) ? std::next(x: pop.cbegin(), n: std::distance(first: cumFitness.cbegin(), last: cumIter))
513 : std::prev(x: pop.cend());
514}
515
516void Genetic::select(Population& pop)
517{
518 std::vector<double> fitness{};
519 fitness.reserve(n: pop.size());
520 std::ranges::transform(pop, std::back_inserter(x&: fitness), [this](const auto& ind) { return calculateFitness(chr: ind); });
521 if (const double min = *std::ranges::min_element(fitness); min <= 0.0)
522 {
523 std::ranges::for_each(fitness, [min, eps = property.prec](auto& fit) { fit = fit - min + eps; });
524 }
525
526 const auto coeff = goldbergLinearScaling(fitness, eps: property.prec)
527 .value_or(u: std::pair<double, double>(0.0, 1.0 / static_cast<double>(pop.size())));
528 const double sum = std::accumulate(
529 first: fitness.begin(),
530 last: fitness.end(),
531 init: 0.0,
532 binary_op: [alpha = coeff.first, beta = coeff.second](const auto acc, auto& fit)
533 {
534 fit = alpha * fit + beta;
535 return acc + fit;
536 });
537 std::ranges::for_each(fitness, [sum](auto& fit) { fit /= sum; });
538 std::partial_sum(first: fitness.begin(), last: fitness.end(), result: fitness.begin());
539
540 Population selected{};
541 selected.reserve(n: pop.size());
542 std::generate_n(
543 first: std::back_inserter(x&: selected),
544 n: pop.size(),
545 gen: [this, &pop, &fitness]() { return *rouletteWheelSelection(pop, cumFitness: fitness); });
546 pop = std::move(selected);
547}
548
549Genetic::Chromosome Genetic::extractElite(Population& pop)
550{
551 const auto bestIter = getBestIndividual(pop);
552 auto elite = *bestIter;
553 pop.erase(position: bestIter);
554 return elite;
555}
556} // namespace algorithm::optimal
557