1#include <iostream>
2
3#include <dcgp/expression.hpp>
4#include <dcgp/kernel_set.hpp>
5
6// Here we search for first integrals of the Kepler's problem) using our "mutation suppression" method
7// The hamiltonian is H = 1/(2m) (pr^2+pt^2 / r^2) + mu / r
8
9double fitness(const dcgp::expression<audi::gdual_d> &ex, const std::vector<std::vector<audi::gdual_d>> &in, double &check)
10{
11 double retval = 0;
12 check = 0;
13 for (auto i = 0u; i < in.size(); ++i) {
14 auto T = ex(in[i]); // We compute all the derivatives up to order one
15 std::vector<std::string> symbol_list{"pr", "pt", "r", "th", "m", "mu"};
16 for (auto sym : symbol_list) {
17 T[0] += audi::gdual_d(0, sym, 0); // We make sure that the symbols are all in the final expression
18 }
19 double dFpr = T[0].get_derivative({{"dpr", 1}});
20 // double dFpt= T[0].get_derivative({{"dpt", 1}});
21 double dFqr = T[0].get_derivative({{"dr", 1}});
22 double dFqt = T[0].get_derivative({{"dth", 1}});
23 double pr = in[i][0].constant_cf();
24 double pt = in[i][1].constant_cf();
25 double qr = in[i][2].constant_cf();
26 // double qt = in[i][3].constant_cf();
27 double m = in[i][4].constant_cf();
28 double mu = in[i][5].constant_cf();
29 double err = dFpr * (pt * pt / m / qr / qr / qr - mu / qr / qr) + dFqr * (pr / m) + dFqt * (pt / qr / qr / m);
30 retval += (err) * (err); // We compute the quadratic error
31 check += dFpr * dFpr + dFqr * dFqr
32 + dFqt * dFqt; // If check is 0 then ex represent a constant or an ignorable variable
33 }
34
35 return retval;
36}
37
38using namespace dcgp;
39
40int main()
41{
42 // Random seed
43 std::random_device rd;
44
45 // Function set
46 dcgp::kernel_set<audi::gdual_d> basic_set({"sum", "diff", "mul", "div"});
47
48 // d-CGP expression
49 dcgp::expression<audi::gdual_d> ex(6, 1, 1, 100, 50, 2, basic_set(), rd());
50
51 // Symbols
52 std::vector<std::string> in_sym({"pr", "pt", "r", "th", "m", "mu"});
53
54 // We create the grid over x
55 std::vector<std::vector<audi::gdual_d>> in(50u);
56 for (auto i = 0u; i < in.size(); ++i) {
57 audi::gdual_d pr_var(0.12 + 20. / static_cast<double>((in.size() - 1)) * i, "pr", 1u);
58 audi::gdual_d pt_var(1. + 20. / static_cast<double>((in.size() - 1)) * i, "pt", 1u);
59 audi::gdual_d r_var(0.12 + 20. / static_cast<double>((in.size() - 1)) * i, "r", 1u);
60 audi::gdual_d th_var(1. + 20. / static_cast<double>((in.size() - 1)) * i, "th", 1u);
61 audi::gdual_d m_var(0.12 + 20. / static_cast<double>((in.size() - 1)) * i, "m", 1u);
62 audi::gdual_d mu_var(1. + 20. / static_cast<double>((in.size() - 1)) * i, "mu", 1u);
63 in[i] = std::vector<audi::gdual_d>{pr_var, pt_var, r_var, th_var, m_var, mu_var};
64 }
65
66 // We run the (1-4)-ES
67 double best_fit = 1e32;
68 std::vector<double> newfits(4, 0.);
69 std::vector<std::vector<unsigned int>> newchromosomes(4);
70 std::vector<unsigned int> best_chromosome(ex.get());
71 unsigned int gen = 0;
72 do {
73 gen++;
74 for (auto i = 0u; i < newfits.size(); ++i) {
75 ex.set(best_chromosome);
76 double check = 0;
77 while (check < 1e-3) {
78 ex.mutate_active(6);
79 newfits[i] = fitness(ex, in, check); // Total fitness
80 }
81 newchromosomes[i] = ex.get();
82 }
83 for (auto i = 0u; i < newfits.size(); ++i) {
84 if (newfits[i] <= best_fit) {
85 if (newfits[i] != best_fit) {
86 audi::print("New best found: gen: ", std::setw(7), gen, "\t value: ", newfits[i], "\n");
87 // std::cout << "Expression: " << ex(in_sym), "\n");
88 }
89 best_fit = newfits[i];
90 best_chromosome = newchromosomes[i];
91 ex.set(best_chromosome);
92 }
93 }
94 } while (best_fit > 1e-12 && gen < 10000);
95
96 audi::print("Number of generations: ", gen, "\n");
97 audi::print("Expression: ", ex, "\n");
98 audi::print("Expression: ", ex(in_sym), "\n");
99 audi::print("Point: ", in[2], "\n");
100 audi::print("Taylor: ", ex(in[2]), "\n");
101}