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