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