Search for first integrals of Kepler’s problem

  1#include <iostream>
  2
  3#include <dcgp/expression.hpp>
  4#include <dcgp/kernel_set.hpp>
  5
  6// Here we search for first integrals of the Kepler's problem) using our "mutation suppression" method
  7// The hamiltonian is H = 1/(2m) (pr^2+pt^2 / r^2) + mu / r
  8
  9double fitness(const dcgp::expression<audi::gdual_d> &ex, const std::vector<std::vector<audi::gdual_d>> &in, double &check)
 10{
 11    double retval = 0;
 12    check = 0;
 13    for (auto i = 0u; i < in.size(); ++i) {
 14        auto T = ex(in[i]); // We compute all the derivatives up to order one
 15        std::vector<std::string> symbol_list{"pr", "pt", "r", "th", "m", "mu"};
 16        for (auto sym : symbol_list) {
 17            T[0] += audi::gdual_d(0, sym, 0); // We make sure that the symbols are all in the final expression
 18        }
 19        double dFpr = T[0].get_derivative({{"dpr", 1}});
 20        // double dFpt= T[0].get_derivative({{"dpt", 1}});
 21        double dFqr = T[0].get_derivative({{"dr", 1}});
 22        double dFqt = T[0].get_derivative({{"dth", 1}});
 23        double pr = in[i][0].constant_cf();
 24        double pt = in[i][1].constant_cf();
 25        double qr = in[i][2].constant_cf();
 26        // double qt = in[i][3].constant_cf();
 27        double m = in[i][4].constant_cf();
 28        double mu = in[i][5].constant_cf();
 29        double err = dFpr * (pt * pt / m / qr / qr / qr - mu / qr / qr) + dFqr * (pr / m) + dFqt * (pt / qr / qr / m);
 30        retval += (err) * (err); // We compute the quadratic error
 31        check += dFpr * dFpr + dFqr * dFqr
 32                 + dFqt * dFqt; // If check is 0 then ex represent a constant or an ignorable variable
 33    }
 34
 35    return retval;
 36}
 37
 38using namespace dcgp;
 39
 40int main()
 41{
 42    // Random seed
 43    std::random_device rd;
 44
 45    // Function set
 46    dcgp::kernel_set<audi::gdual_d> basic_set({"sum", "diff", "mul", "div"});
 47
 48    // d-CGP expression
 49    dcgp::expression<audi::gdual_d> ex(6, 1, 1, 100, 50, 2, basic_set(), rd());
 50
 51    // Symbols
 52    std::vector<std::string> in_sym({"pr", "pt", "r", "th", "m", "mu"});
 53
 54    // We create the grid over x
 55    std::vector<std::vector<audi::gdual_d>> in(50u);
 56    for (auto i = 0u; i < in.size(); ++i) {
 57        audi::gdual_d pr_var(0.12 + 20. / static_cast<double>((in.size() - 1)) * i, "pr", 1u);
 58        audi::gdual_d pt_var(1. + 20. / static_cast<double>((in.size() - 1)) * i, "pt", 1u);
 59        audi::gdual_d r_var(0.12 + 20. / static_cast<double>((in.size() - 1)) * i, "r", 1u);
 60        audi::gdual_d th_var(1. + 20. / static_cast<double>((in.size() - 1)) * i, "th", 1u);
 61        audi::gdual_d m_var(0.12 + 20. / static_cast<double>((in.size() - 1)) * i, "m", 1u);
 62        audi::gdual_d mu_var(1. + 20. / static_cast<double>((in.size() - 1)) * i, "mu", 1u);
 63        in[i] = std::vector<audi::gdual_d>{pr_var, pt_var, r_var, th_var, m_var, mu_var};
 64    }
 65
 66    // We run the (1-4)-ES
 67    double best_fit = 1e32;
 68    std::vector<double> newfits(4, 0.);
 69    std::vector<std::vector<unsigned int>> newchromosomes(4);
 70    std::vector<unsigned int> best_chromosome(ex.get());
 71    unsigned int gen = 0;
 72    do {
 73        gen++;
 74        for (auto i = 0u; i < newfits.size(); ++i) {
 75            ex.set(best_chromosome);
 76            double check = 0;
 77            while (check < 1e-3) {
 78                ex.mutate_active(6);
 79                newfits[i] = fitness(ex, in, check); // Total fitness
 80            }
 81            newchromosomes[i] = ex.get();
 82        }
 83        for (auto i = 0u; i < newfits.size(); ++i) {
 84            if (newfits[i] <= best_fit) {
 85                if (newfits[i] != best_fit) {
 86                    audi::print("New best found: gen: ", std::setw(7), gen, "\t value: ", newfits[i], "\n");
 87                    // std::cout << "Expression: " << ex(in_sym), "\n");
 88                }
 89                best_fit = newfits[i];
 90                best_chromosome = newchromosomes[i];
 91                ex.set(best_chromosome);
 92            }
 93        }
 94    } while (best_fit > 1e-12 && gen < 10000);
 95
 96    audi::print("Number of generations: ", gen, "\n");
 97    audi::print("Expression: ", ex, "\n");
 98    audi::print("Expression: ", ex(in_sym), "\n");
 99    audi::print("Point: ", in[2], "\n");
100    audi::print("Taylor: ", ex(in[2]), "\n");
101}