1#include <iostream>
2
3#include <dcgp/expression.hpp>
4#include <dcgp/kernel_set.hpp>
5
6// Here we search for first integrals of a mass spring sistem (one dimension) using Lipson method
7// The hamiltonian is H = p^2 + q^2 and is consistently found by the evolution
8
9double fitness(const dcgp::expression<audi::gdual_d> &ex, const std::vector<std::vector<audi::gdual_d>> &in)
10{
11 double retval = 0;
12 for (auto i = 0u; i < in.size(); ++i) {
13 auto T = ex(in[i]); // We compute all the derivatives up to order one
14 double dFp = T[0].get_derivative({{"dp", 1}});
15 double dFq = T[0].get_derivative({{"dq", 1}});
16 double p = in[i][0].constant_cf();
17 double q = in[i][1].constant_cf();
18 double err = dFp / dFq - p / q; // Here we set (dp/dt) / (dq/dt) = dp/dq
19 retval += std::log(1 + std::abs(err));
20 }
21
22 return retval / static_cast<double>(in.size());
23}
24
25using namespace dcgp;
26
27int main()
28{
29 // Random seed
30 std::random_device rd;
31
32 // Function set
33 dcgp::kernel_set<audi::gdual_d> basic_set({"sum", "diff", "mul", "div"});
34
35 // d-CGP expression
36 dcgp::expression<audi::gdual_d> ex(2, 1, 1, 15, 16, 2, basic_set(), rd());
37
38 // Symbols
39 std::vector<std::string> in_sym({"p", "q"});
40
41 // We create the grid over x
42 std::vector<std::vector<audi::gdual_d>> in(10u);
43 for (auto i = 0u; i < in.size(); ++i) {
44 audi::gdual_d p_var(0.12 + 0.9 / static_cast<double>((in.size() - 1)) * i, "p", 1u);
45 audi::gdual_d q_var(1. - 0.143 / static_cast<double>((in.size() - 1)) * i, "q", 1u);
46 in[i] = std::vector<audi::gdual_d>{p_var, q_var};
47 }
48 // We run the (1-4)-ES
49 double best_fit = 1e32;
50 std::vector<double> newfits(4, 0.);
51 std::vector<std::vector<unsigned int>> newchromosomes(4);
52 std::vector<unsigned int> best_chromosome(ex.get());
53 unsigned int gen = 0;
54 do {
55 gen++;
56 for (auto i = 0u; i < newfits.size(); ++i) {
57 ex.set(best_chromosome);
58 ex.mutate_active(6);
59 newfits[i] = fitness(ex, in); // Total fitness
60 newchromosomes[i] = ex.get();
61 }
62 for (auto i = 0u; i < newfits.size(); ++i) {
63 if (newfits[i] <= best_fit) {
64 if (newfits[i] != best_fit) {
65 std::cout << "New best found: gen: " << std::setw(7) << gen << "\t value: " << newfits[i]
66 << std::endl;
67 // std::cout << "Expression: " << ex(in_sym) << std::endl;
68 }
69 best_fit = newfits[i];
70 best_chromosome = newchromosomes[i];
71 ex.set(best_chromosome);
72 }
73 }
74 } while (best_fit > 1e-12 && gen < 10000);
75
76 audi::print("Number of generations: ", gen, "\n");
77 audi::print("Expression: ", ex, "\n");
78 audi::print("Expression: ", ex(in_sym), "\n");
79 audi::print("Point: ", in[2], "\n");
80 audi::print("Taylor: ", ex(in[2]), "\n");
81}