1#include <audi/io.hpp>
2#include <iostream>
3
4#include <dcgp/expression.hpp>
5#include <dcgp/kernel_set.hpp>
6
7// Here we solve the differential equation dy = (2x - y) / x from Tsoulos paper
8// Tsoulos and Lagaris: "Solving Differential equations with genetic programming"
9
10double fitness(const dcgp::expression<audi::gdual_d> &ex, const std::vector<std::vector<audi::gdual_d>> &in)
11{
12 double retval = 0;
13 for (auto i = 0u; i < in.size(); ++i) {
14 auto T = ex(in[i]); // We compute the expression and thus the derivatives
15 double y = T[0].get_derivative({0});
16 double dy = T[0].get_derivative({1});
17 double x = in[i][0].constant_cf();
18 double ode1 = (2. * x - y) / x;
19 retval += (ode1 - dy) * (ode1 - dy); // We compute the quadratic error
20 }
21 return retval;
22}
23
24int main()
25{
26 // Random seed
27 std::random_device rd;
28
29 // Function set
30 dcgp::kernel_set<audi::gdual_d> basic_set({"sum", "diff", "mul", "div", "exp", "log", "sin", "cos"});
31
32 // d-CGP expression
33 dcgp::expression<audi::gdual_d> ex(1, 1, 1, 15, 16, 2, basic_set(), rd());
34
35 // Symbols
36 std::vector<std::string> in_sym({"x"});
37
38 // We create the grid over x
39 std::vector<std::vector<audi::gdual_d>> in(10u);
40 for (auto i = 0u; i < in.size(); ++i) {
41 audi::gdual_d point(0.1 + 0.9 / static_cast<double>((in.size() - 1)) * i, "x", 1);
42 in[i].push_back(point); // 1, .., 2
43 }
44
45 // We run the (1-4)-ES
46 auto best_fit = 1e32;
47 std::vector<double> newfits(4, 0.);
48 std::vector<std::vector<unsigned int>> newchromosomes(4);
49 auto best_chromosome(ex.get());
50 auto gen = 0u;
51
52 do {
53 gen++;
54 for (auto i = 0u; i < newfits.size(); ++i) {
55 ex.set(best_chromosome);
56 ex.mutate_active(2);
57 auto fitness_ic = ex({audi::gdual_d(1.)})[0] - 3.; // Penalty term to enforce the initial conditions
58 newfits[i] = fitness(ex, in) + fitness_ic.constant_cf() * fitness_ic.constant_cf(); // Total fitness
59 newchromosomes[i] = ex.get();
60 }
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-3 && gen < 3000);
75 audi::stream(std::cout, "Number of generations: ", gen, "\n");
76 audi::stream(std::cout, "Expression: ", ex(in_sym), "\n");
77}