Solving the ODE1 problem from Tsoulos

 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}