Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2019 Jakub Rydzewski (jr@fizyka.umk.pl). All rights reserved.
3 :
4 : See http://www.maze-code.github.io for more information.
5 :
6 : This file is part of maze.
7 :
8 : maze is free software: you can redistribute it and/or modify it under the
9 : terms of the GNU Lesser General Public License as published by the Free
10 : Software Foundation, either version 3 of the License, or (at your option)
11 : any later version.
12 :
13 : maze is distributed in the hope that it will be useful, but WITHOUT ANY
14 : WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15 : FOR A PARTICULAR PURPOSE.
16 :
17 : See the GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with maze. If not, see <https://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #ifndef __PLUMED_maze_Memetic_h
23 : #define __PLUMED_maze_Memetic_h
24 :
25 : /**
26 : * @file Memetic.h
27 : *
28 : * @author J. Rydzewski (jr@fizyka.umk.pl)
29 : */
30 :
31 : #include "core/ActionRegister.h"
32 :
33 : #include "Core.h"
34 : #include "Member.h"
35 : #include "Optimizer.h"
36 :
37 : namespace PLMD {
38 : namespace maze {
39 :
40 : /**
41 : * @class Memetic Memetic.h "maze/Memetic.h"
42 : *
43 : * @brief Memetic algorithms for the optimization of the loss function.
44 : */
45 : class Memetic: public Optimizer {
46 : public:
47 : /**
48 : * PLMD constructor.
49 : *
50 : * @param[in] ao PLMD::ActionOptions&.
51 : */
52 : explicit Memetic(const ActionOptions& ao);
53 :
54 : /**
55 : * Registers PLMD keywords.
56 : *
57 : * @param[in] keys Keywords.
58 : */
59 : static void registerKeywords(Keywords& keys);
60 :
61 : /**
62 : * Each class deriving from Optimizer needs to override this function.
63 : */
64 : void optimize() override;
65 :
66 : private:
67 : /**
68 : * Create a set of translations relative to the ligand, each translation
69 : * encodes a probable biasing direction.
70 : */
71 : void initialize_members();
72 :
73 : /**
74 : * Score each translation by the loss function.
75 : */
76 : void score_members();
77 :
78 : /**
79 : * Calculate the mean score.
80 : */
81 : double score_mean();
82 :
83 : /**
84 : * Sort the population using heaps, required for finding a minimum of the loss
85 : * function.
86 : */
87 : void sort_members();
88 :
89 : /**
90 : * Encode a ligand conformation.
91 : */
92 : Vector create_coding();
93 :
94 : /**
95 : * Check if the vector length is out of bounds.
96 : *
97 : * @param[in] v Vector length.
98 : */
99 : bool out_of_bounds(double v);
100 :
101 : /**
102 : * Score a single member.
103 : *
104 : * @param[in] v Member's translation.
105 : */
106 : double score_member(const Vector& v);
107 :
108 : /**
109 : * Print a status.
110 : */
111 : void print_status() const;
112 :
113 : /**
114 : * Select a new population using the roulette selection.
115 : */
116 : void selection_roulette();
117 :
118 : /**
119 : * Perform mating in the population.
120 : *
121 : * @param[in,out] members Population.
122 : */
123 : void mating(std::vector<Member>& members);
124 :
125 : /**
126 : * Mate two members.
127 : *
128 : * @param[in,out] m1 1st member.
129 : * @param[in,out] m2 2nd member.
130 : */
131 : void crossover(Member& m1, Member& m2);
132 :
133 : /**
134 : * Mutate a member.
135 : *
136 : * @param[in,out] m Member.
137 : */
138 : void mutation(Member& m);
139 :
140 : /**
141 : * Mutate the population.
142 : *
143 : * @param[in,out] members Population.
144 : */
145 : void mutation(std::vector<Member>& members);
146 :
147 : protected:
148 : /**
149 : * Local searches to improve global solutions.
150 : */
151 :
152 : /**
153 : * Stochastic hill climbing.
154 : *
155 : * Neighbors with better or equal cost should be accepted, allowing the
156 : * technique to navigate across plateaus in the response surface.
157 : *
158 : * @param[in,out] m Member.
159 : * @param[in] params None.
160 : */
161 : void stochastic_hill_climbing(
162 : Member& m,
163 : /* none */ const std::vector<double>& params = {}
164 : );
165 :
166 : /**
167 : * Random-restart hill climbing.
168 : *
169 : * The algorithm can be restarted and repeated a number of times after it
170 : * converges to provide an improved result.
171 : *
172 : * @param[in,out] m Member.
173 : * @param[in] params Number of restarts.
174 : */
175 : void random_restart_hill_climbing(
176 : Member& m,
177 : /* n_restarts */ const std::vector<double>& params = {10}
178 : );
179 :
180 : /**
181 : * Solis-Wets random walk.
182 : *
183 : * Adaptive random search algorithm was designed to address the limitations of
184 : * the fixed step size. The strategy for adaptive random search is to
185 : * continually approximate the optimal step size required to reach the global
186 : * optimum in the search space. This is achieved by trialling and adopting
187 : * smaller or larger step sizes only if they result in an improvement in the
188 : * search performance.
189 : *
190 : * Very large step sizes are trialled with a much lower frequency. This
191 : * strategy of preferring large moves is intended to allow the technique to
192 : * escape local optima. Smaller step sizes are adopted if no improvement is
193 : * made for an extended period.
194 : *
195 : * @param[in,out] m Member.
196 : * @param[in] params
197 : */
198 : void adaptive_random_search(
199 : Member& m,
200 : /* */ const std::vector<double>& params = {1.0, 1.e-5, 2.0, 2.0, 3.0, 3.0}
201 : );
202 :
203 : /**
204 : * Luus–Jaakola heuristics.
205 : *
206 : * @param[in,out] m Member.
207 : * @param[in] params Bounds.
208 : */
209 : void luus_jaakola(
210 : Member& m,
211 : /* bounds */ const std::vector<double>& params
212 : );
213 :
214 : /**
215 : * Local annealing.
216 : *
217 : * @param[in,out] m Member.
218 : * @param[in] params T, alpha.
219 : */
220 : void annealing(
221 : Member& m,
222 : /* T, alpha */const std::vector<double>& params = {300.0, 0.95}
223 : );
224 :
225 : /**
226 : * Apply local search to members.
227 : *
228 : * @param[in,out] members Population.
229 : */
230 : void local_search(std::vector<Member>& members);
231 :
232 : protected:
233 : /**
234 : * Return an optimal biasing direction.
235 : */
236 : Vector solve();
237 :
238 : public:
239 : /**
240 : * Setters and getters.
241 : */
242 :
243 : unsigned int get_capacity() const;
244 : void set_capacity(unsigned int);
245 :
246 : unsigned int get_coding_len() const;
247 : void set_coding_len(unsigned int);
248 :
249 : unsigned int get_n_local_iterations() const;
250 : void set_n_local_iterations(unsigned int);
251 :
252 : unsigned int get_n_global_iterations() const;
253 : void set_n_global_iterations(unsigned int);
254 :
255 : double get_mutation_rate() const;
256 : void set_mutation_rate(double);
257 :
258 : double get_mating_rate() const;
259 : void set_mating_rate(double);
260 :
261 : double get_cauchy_mean() const;
262 : void set_cauchy_mean(double);
263 :
264 : double get_cauchy_spread() const;
265 : void set_cauchy_spread(double);
266 :
267 : bool is_local_search_on() const;
268 : void local_search_on();
269 : void local_search_off();
270 :
271 : double get_local_search_rate() const;
272 : void set_local_search_rate(double);
273 :
274 : std::string get_local_search_type() const;
275 : void set_local_search_type(const std::string&);
276 :
277 : protected:
278 : //! Population
279 : std::vector<Member> members_;
280 :
281 : //! Bound
282 : double bound_;
283 :
284 : //! Scores
285 : double score_worst_;
286 : double score_best_;
287 : Member member_best_;
288 :
289 : //! If a local search is performed
290 : bool adaptation_;
291 :
292 : protected:
293 : //! Size of population
294 : unsigned int capacity_;
295 : //! Length of coding
296 : unsigned int coding_len_;
297 :
298 : //! Number of local search iterations
299 : unsigned int n_local_iterations_;
300 : //! Number of global search iterations, doomsday
301 : unsigned int n_global_iterations_;
302 :
303 : //! Probability of mutation
304 : double mutation_rate_;
305 : //! Probability of mating
306 : double mating_rate_;
307 :
308 : //! Mean and spread of cauchy sampler
309 : double cauchy_mean_alpha_;
310 : double cauchy_mean_beta_;
311 :
312 : //! If local search is employed in sampling
313 : bool local_search_on_;
314 : //! Rate of local mutation
315 : double local_search_rate_;
316 : //! Type of local search, stochastic_hill_climbing or adaptive_random_search
317 : std::string local_search_type_;
318 : };
319 :
320 6 : void Memetic::initialize_members() {
321 6 : members_.clear();
322 6 : members_.resize(capacity_);
323 :
324 42 : for (size_t i = 0; i < capacity_; ++i) {
325 36 : Member m{};
326 :
327 36 : m.score=0;
328 36 : m.translation=create_coding();
329 :
330 36 : members_.at(i) = m;
331 : }
332 6 : }
333 :
334 30 : void Memetic::score_members() {
335 210 : for (size_t i = 0; i < members_.size(); ++i) {
336 180 : double s = score_member(members_[i].translation);
337 180 : members_[i].score=s;
338 : }
339 30 : }
340 :
341 30 : void Memetic::sort_members() {
342 30 : std::make_heap(
343 : members_.begin(),
344 : members_.end(),
345 : compare
346 : );
347 :
348 30 : std::sort_heap(
349 : members_.begin(),
350 : members_.end(),
351 : compare
352 : );
353 :
354 30 : member_best_ = members_[capacity_ - 1];
355 30 : score_best_ = members_[capacity_ - 1].score;
356 30 : score_worst_ = members_[0].score;
357 30 : }
358 :
359 0 : double Memetic::score_mean() {
360 : auto acc = [](double s, const Member& m) {
361 0 : return s + m.score;
362 : };
363 :
364 : return std::accumulate(
365 : members_.begin(),
366 : members_.end(),
367 : 0.0,
368 0 : acc) / capacity_;
369 : }
370 :
371 30 : void Memetic::selection_roulette() {
372 30 : std::vector<Member> sel(members_);
373 30 : std::vector<double> rel_scores(capacity_, 0.0);
374 :
375 210 : for (std::size_t i = 0; i < capacity_; ++i) {
376 180 : double r = 1.0 / (members_[i].score + 0.01);
377 180 : rel_scores.at(i) = r;
378 : }
379 :
380 30 : std::vector<double> cum_sum(capacity_, 0.0);
381 30 : std::partial_sum(
382 : rel_scores.begin(),
383 : rel_scores.end(),
384 : cum_sum.begin(),
385 : std::plus<double>()
386 : );
387 :
388 30 : double sum = cum_sum.back();
389 : members_.clear();
390 30 : members_.resize(capacity_);
391 : int chosen = -1;
392 :
393 210 : for (size_t j = 0; j < capacity_; ++j) {
394 : double probability=rnd::next_double(sum);
395 639 : for (size_t i = 0; i < cum_sum.size(); ++i) {
396 639 : if (cum_sum[i] > probability) {
397 180 : chosen = i;
398 :
399 180 : members_.at(j).score = sel.at(chosen).score;
400 180 : members_.at(j).translation = sel.at(chosen).translation;
401 :
402 180 : break;
403 : }
404 : }
405 : }
406 30 : }
407 :
408 0 : void Memetic::crossover(Member& s1, Member& s2) {
409 0 : size_t i = rnd::next_int(1, coding_len_ - 1);
410 :
411 0 : Member z1(s1);
412 0 : Member z2(s2);
413 :
414 0 : for (size_t j = i; j < coding_len_; ++j) {
415 0 : z1.translation[j] = s2.translation[j];
416 0 : z2.translation[j] = s1.translation[j];
417 : }
418 :
419 0 : if (!out_of_bounds(z1.translation.modulo()) && !out_of_bounds(z2.translation.modulo())) {
420 0 : s1 = z1;
421 0 : s2 = z2;
422 : }
423 0 : }
424 :
425 113 : void Memetic::mutation(Member& m) {
426 113 : int which = rnd::next_int(coding_len_);
427 113 : double v = rnd::next_cauchy(cauchy_mean_alpha_, cauchy_mean_beta_);
428 113 : m.translation[which] += v;
429 113 : if (out_of_bounds(m.translation.modulo())) {
430 35 : m.translation[which] -= v;
431 : }
432 113 : }
433 :
434 30 : void Memetic::mutation(std::vector<Member>& m) {
435 210 : for (std::vector<Member>::iterator it = m.begin(), end = m.end(); it != end; ++it) {
436 180 : double r = rnd::next_double();
437 180 : if (r < mutation_rate_) {
438 3 : mutation(*it);
439 : }
440 : }
441 30 : }
442 :
443 11 : void Memetic::stochastic_hill_climbing(
444 : Member& m,
445 : const std::vector<double>& params) {
446 121 : for (std::size_t i = 0; i < n_local_iterations_; ++i) {
447 110 : Member n;
448 110 : n.translation = m.translation;
449 110 : mutation(n);
450 110 : double score_n = score_member(n.translation);
451 :
452 110 : if (m.score > score_n) {
453 52 : m.translation = n.translation;
454 : }
455 : }
456 11 : }
457 :
458 0 : void Memetic::random_restart_hill_climbing(
459 : Member& m,
460 : const std::vector<double>& params) {
461 0 : unsigned int n_restarts = static_cast<int>(params[0]);
462 0 : std::vector<Member> s(n_restarts);
463 : unsigned int ndx = 0;
464 0 : double min = m.score;
465 :
466 0 : for (std::size_t r = 0; r < n_restarts; ++r) {
467 0 : Member n = m;
468 0 : stochastic_hill_climbing(n);
469 0 : s[r] = n;
470 :
471 0 : if (min > n.score) {
472 : min = n.score;
473 0 : ndx = r;
474 : }
475 : }
476 :
477 0 : m = s[ndx];
478 0 : }
479 :
480 0 : void Memetic::annealing(
481 : Member& m,
482 : const std::vector<double>& params) {
483 0 : double T = params[0];
484 0 : double alpha = params[1];
485 :
486 0 : for (std::size_t i = 0; i < n_local_iterations_; ++i) {
487 0 : Member n = m;
488 0 : mutation(n);
489 0 : double score_n = score_member(n.translation);
490 :
491 : double probability = std::min(
492 0 : 1.0,
493 0 : std::exp(-(score_n - m.score) / T)
494 0 : );
495 :
496 0 : double r = rnd::next_double();
497 :
498 0 : if (r < probability) {
499 0 : m = n;
500 : }
501 :
502 0 : T *= alpha;
503 : }
504 0 : }
505 :
506 0 : void Memetic::luus_jaakola(
507 : Member& m,
508 : const std::vector<double>& params) {
509 : /* TODO */
510 0 : }
511 :
512 0 : void Memetic::adaptive_random_search(
513 : Member& m,
514 : const std::vector<double>& params) {
515 0 : double rho_start = params[0];
516 0 : double rho_lower_bound = params[1];
517 0 : double ex = params[2];
518 0 : double ct = params[3];
519 0 : int s_ex = static_cast<int>(params[4]);
520 0 : int f_ct = static_cast<int>(params[5]);
521 :
522 : unsigned int k = 0;
523 0 : Vector xk = m.translation;
524 : int ns = 0;
525 : int nf = 0;
526 0 : Vector bk;
527 0 : bk.zero();
528 : double rho_k = rho_start;
529 :
530 0 : while (rho_k > rho_lower_bound && k < n_local_iterations_) {
531 0 : if (ns >= s_ex) {
532 0 : rho_k *= ex;
533 0 : } else if (nf >= f_ct) {
534 0 : rho_k *= ct;
535 : }
536 :
537 0 : std::vector<double> chiv = rnd::next_double(-rho_k, rho_k, 3);
538 0 : Vector chi = tls::vector2Vector(chiv);
539 0 : Vector tmp;
540 :
541 0 : for (unsigned int i = 0; i < 3; ++i) {
542 0 : tmp[i] = 2.0 * (xk[i] - chi[i]);
543 : }
544 :
545 0 : double score_chi = score_member(chi);
546 0 : double score_xk = score_member(xk);
547 0 : double score_tmp = score_member(tmp);
548 :
549 0 : if (score_chi < score_xk) {
550 0 : ns++;
551 : nf = 0;
552 :
553 0 : for (unsigned int i = 0; i < 3; i++) {
554 0 : bk[i] = 0.2 * bk[i] + 0.4 * (chi[i] - xk[i]);
555 0 : xk[i] = chi[i];
556 : }
557 0 : } else if (score_tmp < score_xk && score_xk <= score_chi) {
558 0 : ns++;
559 : nf = 0;
560 :
561 0 : for (unsigned int i = 0; i < 3; i++) {
562 0 : bk[i] = bk[i] - 0.4 * (chi[i] - xk[i]);
563 0 : xk[i] = 2.0 * xk[i] - chi[i];
564 : }
565 : } else {
566 : ns = 0;
567 0 : nf++;
568 :
569 0 : for (unsigned int i = 0; i < 3; i++) {
570 0 : bk[i] = 0.5 * bk[i];
571 : }
572 : }
573 :
574 0 : k++;
575 : }
576 :
577 0 : m.translation = xk;
578 0 : }
579 :
580 30 : void Memetic::local_search(std::vector<Member>& m) {
581 30 : adaptation_ = true;
582 :
583 30 : if (local_search_on_) {
584 105 : for (std::size_t i = 0; i < capacity_; ++i) {
585 90 : double probability = rnd::next_double();
586 :
587 90 : if (probability < local_search_rate_) {
588 11 : if (local_search_type_ == "stochastic_hill_climbing") {
589 22 : stochastic_hill_climbing(m[i]);
590 0 : } else if (local_search_type_ == "adaptive_random_search") {
591 0 : adaptive_random_search(m[i]);
592 0 : } else if (local_search_type_ == "random_restart_hill_climbing") {
593 0 : random_restart_hill_climbing(m[i]);
594 : }
595 : }
596 : }
597 : }
598 :
599 30 : adaptation_ = false;
600 30 : }
601 :
602 30 : void Memetic::mating(std::vector<Member>& m) {
603 : std::vector<Member> offspring;
604 :
605 120 : while (m.size() != 0) {
606 90 : int j = rnd::next_int(m.size());
607 90 : int i = rnd::next_int(m.size());
608 :
609 139 : while (i == j) {
610 49 : j=rnd::next_int(m.size());
611 : }
612 :
613 90 : Member m1 = m[i];
614 90 : Member m2 = m[j];
615 :
616 90 : if (i > j) {
617 : m.erase(m.begin() + i);
618 : m.erase(m.begin() + j);
619 58 : } else if (j > i) {
620 : m.erase(m.begin() + j);
621 : m.erase(m.begin() + i);
622 : }
623 :
624 90 : double probability = rnd::next_double();
625 90 : if (probability < mating_rate_) {
626 0 : crossover(m1, m2);
627 : }
628 :
629 90 : offspring.push_back(m1);
630 90 : offspring.push_back(m2);
631 : }
632 :
633 30 : m = offspring;
634 : offspring.clear();
635 30 : }
636 :
637 36 : Vector Memetic::create_coding() {
638 36 : double s = sampling_radius();
639 : double r = rnd::next_double(s);
640 :
641 36 : return rnd::next_plmd_vector(r);
642 : }
643 :
644 113 : bool Memetic::out_of_bounds(double v) {
645 113 : double s = sampling_radius();
646 :
647 113 : return v > s;
648 : }
649 :
650 290 : double Memetic::score_member(const Vector& coding) {
651 : double action = 0;
652 290 : Vector distance;
653 290 : const unsigned nl_size = neighbor_list_->size();
654 290 : Vector dev = coding;
655 :
656 290 : #pragma omp parallel num_threads(get_n_threads_openmp())
657 : {
658 : #pragma omp for reduction(+:action)
659 : for (unsigned int i = 0; i < nl_size; i++) {
660 : unsigned i0 = neighbor_list_->getClosePair(i).first;
661 : unsigned i1 = neighbor_list_->getClosePair(i).second;
662 :
663 : if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
664 : continue;
665 : }
666 :
667 : if (pbc_) {
668 : distance = pbcDistance(getPosition(i0) + dev, getPosition(i1));
669 : } else {
670 : distance = delta(getPosition(i0) + dev, getPosition(i1));
671 : }
672 :
673 : action += pairing(distance.modulo());
674 : }
675 : }
676 :
677 290 : return action;
678 : }
679 :
680 0 : void Memetic::print_status() const {
681 0 : log.printf("Lowest score: %f \n", score_best_);
682 0 : }
683 :
684 6 : Vector Memetic::solve() {
685 6 : initialize_members();
686 :
687 : size_t epoch = 0;
688 36 : while (epoch < n_global_iterations_) {
689 30 : score_members();
690 :
691 30 : selection_roulette();
692 30 : mating(members_);
693 30 : mutation(members_);
694 30 : local_search(members_);
695 :
696 30 : sort_members();
697 :
698 30 : epoch++;
699 : }
700 :
701 6 : return member_best_.translation / member_best_.translation.modulo();
702 : }
703 :
704 : inline unsigned int Memetic::get_capacity() const {
705 : return capacity_;
706 : }
707 :
708 : inline void Memetic::set_capacity(unsigned int capacity) {
709 : capacity_ = capacity;
710 : }
711 :
712 : inline unsigned int Memetic::get_coding_len() const {
713 : return coding_len_;
714 : }
715 :
716 : inline void Memetic::set_coding_len(unsigned int coding_len) {
717 : coding_len_ = coding_len;
718 : }
719 :
720 : inline unsigned int Memetic::get_n_global_iterations() const {
721 : return n_global_iterations_;
722 : }
723 :
724 : inline void Memetic::set_n_global_iterations(unsigned int n_global_iterations) {
725 2 : n_global_iterations_ = n_global_iterations;
726 : }
727 :
728 : inline unsigned int Memetic::get_n_local_iterations() const {
729 : return n_local_iterations_;
730 : }
731 :
732 : inline void Memetic::set_n_local_iterations(unsigned int n_local_iterations) {
733 : n_local_iterations_ = n_local_iterations;
734 : }
735 :
736 : inline double Memetic::get_mating_rate() const {
737 : return mating_rate_;
738 : }
739 :
740 : inline void Memetic::set_mating_rate(double mating_rate) {
741 : mating_rate_ = mating_rate;
742 : }
743 :
744 : inline double Memetic::get_mutation_rate() const {
745 : return mutation_rate_;
746 : }
747 :
748 : inline void Memetic::set_mutation_rate(double mutation_rate) {
749 : mutation_rate_ = mutation_rate;
750 : }
751 :
752 : inline double Memetic::get_cauchy_mean() const {
753 : return cauchy_mean_alpha_;
754 : }
755 :
756 : inline void Memetic::set_cauchy_mean(double cauchy_mean_alpha) {
757 : cauchy_mean_alpha_ = cauchy_mean_alpha;
758 : }
759 :
760 : inline double Memetic::get_cauchy_spread() const {
761 : return cauchy_mean_beta_;
762 : }
763 :
764 : inline void Memetic::set_cauchy_spread(double cauchy_mean_beta) {
765 : cauchy_mean_beta_ = cauchy_mean_beta;
766 : }
767 :
768 : inline bool Memetic::is_local_search_on() const {
769 : return local_search_on_;
770 : }
771 :
772 : inline void Memetic::local_search_on() {
773 : local_search_on_ = true;
774 : }
775 :
776 : inline void Memetic::local_search_off() {
777 : local_search_on_ = false;
778 : }
779 :
780 : inline double Memetic::get_local_search_rate() const {
781 : return local_search_rate_;
782 : }
783 :
784 : inline void Memetic::set_local_search_rate(double local_search_rate) {
785 : local_search_rate_ = local_search_rate;
786 : }
787 :
788 : inline std::string Memetic::get_local_search_type() const {
789 : return local_search_type_;
790 : }
791 :
792 : inline void Memetic::set_local_search_type(const std::string& local_search_type) {
793 : local_search_type_ = local_search_type;
794 : }
795 :
796 : } // namespace maze
797 : } // namespace PLMD
798 :
799 : #endif // __PLUMED_maze_Memetic_h
|
