1#include <audi/audi.hpp>
2#include <iostream>
3
4#include <dcgp/expression.hpp>
5#include <dcgp/kernel_set.hpp>
6
7// Here we solve the differential equation d^2y dy = - 4 / x^3 (NLODE3) 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 dy = T[0].get_derivative({1});
16 double ddy = T[0].get_derivative({2});
17 double x = in[i][0].constant_cf();
18 double ode1 = -4 / x / x / x;
19 retval += (ode1 - ddy * dy) * (ode1 - ddy * 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", "log"});
31
32 // d-CGP expression
33 dcgp::expression<audi::gdual_d> ex(1, 1, 1, 15, 16, 2, basic_set(), rd());
34
35 // Symbols for streaming out the hr expression
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 in[i].push_back(audi::gdual_d(1. + 1. / static_cast<double>((in.size() - 1)) * i, "x", 2)); // 1, .., 2
42 }
43
44 // We run the (1-4)-ES
45 double best_fit = 1e32;
46 std::vector<double> newfits(4, 0.);
47 std::vector<std::vector<unsigned int>> newchromosomes(4);
48 std::vector<unsigned int> best_chromosome(ex.get());
49 unsigned int gen = 0;
50
51 do {
52 gen++;
53 for (auto i = 0u; i < newfits.size(); ++i) {
54 ex.set(best_chromosome);
55 ex.mutate_active(2);
56 auto fitness_ic
57 = ex(std::vector<audi::gdual_d>{audi::gdual_d(1.)})[0]; // 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 < 10000);
75 audi::stream(std::cout, "Number of generations: ", gen, "\n");
76 audi::stream(std::cout, "Expression: ", ex(in_sym), "\n");
77}