Finding the prime integral for a spring mass system using Lipson method

 1#include <iostream>
 2
 3#include <dcgp/expression.hpp>
 4#include <dcgp/kernel_set.hpp>
 5
 6// Here we search for first integrals of a mass spring sistem (one dimension) using Lipson method
 7// The hamiltonian is H = p^2 + q^2 and is consistently found by the evolution
 8
 9double fitness(const dcgp::expression<audi::gdual_d> &ex, const std::vector<std::vector<audi::gdual_d>> &in)
10{
11    double retval = 0;
12    for (auto i = 0u; i < in.size(); ++i) {
13        auto T = ex(in[i]); // We compute all the derivatives up to order one
14        double dFp = T[0].get_derivative({{"dp", 1}});
15        double dFq = T[0].get_derivative({{"dq", 1}});
16        double p = in[i][0].constant_cf();
17        double q = in[i][1].constant_cf();
18        double err = dFp / dFq - p / q; // Here we set (dp/dt) / (dq/dt) = dp/dq
19        retval += std::log(1 + std::abs(err));
20    }
21
22    return retval / static_cast<double>(in.size());
23}
24
25using namespace dcgp;
26
27int main()
28{
29    // Random seed
30    std::random_device rd;
31
32    // Function set
33    dcgp::kernel_set<audi::gdual_d> basic_set({"sum", "diff", "mul", "div"});
34
35    // d-CGP expression
36    dcgp::expression<audi::gdual_d> ex(2, 1, 1, 15, 16, 2, basic_set(), rd());
37
38    // Symbols
39    std::vector<std::string> in_sym({"p", "q"});
40
41    // We create the grid over x
42    std::vector<std::vector<audi::gdual_d>> in(10u);
43    for (auto i = 0u; i < in.size(); ++i) {
44        audi::gdual_d p_var(0.12 + 0.9 / static_cast<double>((in.size() - 1)) * i, "p", 1u);
45        audi::gdual_d q_var(1. - 0.143 / static_cast<double>((in.size() - 1)) * i, "q", 1u);
46        in[i] = std::vector<audi::gdual_d>{p_var, q_var};
47    }
48    // We run the (1-4)-ES
49    double best_fit = 1e32;
50    std::vector<double> newfits(4, 0.);
51    std::vector<std::vector<unsigned int>> newchromosomes(4);
52    std::vector<unsigned int> best_chromosome(ex.get());
53    unsigned int gen = 0;
54    do {
55        gen++;
56        for (auto i = 0u; i < newfits.size(); ++i) {
57            ex.set(best_chromosome);
58            ex.mutate_active(6);
59            newfits[i] = fitness(ex, in); // Total fitness
60            newchromosomes[i] = ex.get();
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-12 && gen < 10000);
75
76    audi::print("Number of generations: ", gen, "\n");
77    audi::print("Expression: ", ex, "\n");
78    audi::print("Expression: ", ex(in_sym), "\n");
79    audi::print("Point: ", in[2], "\n");
80    audi::print("Taylor: ", ex(in[2]), "\n");
81}