Finding the prime integral for a spring mass system

 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}