<commit_before>﻿
#include <random>

#include "astronomy/frames.hpp"
#include "base/bundle.hpp"
#include "base/file.hpp"
#include "base/not_null.hpp"
#include "base/status.hpp"
#include "geometry/named_quantities.hpp"
#include "geometry/sign.hpp"
#include "gtest/gtest.h"
#include "integrators/methods.hpp"
#include "integrators/symmetric_linear_multistep_integrator.hpp"
#include "mathematica/mathematica.hpp"
#include "numerics/root_finders.hpp"
#include "physics/degrees_of_freedom.hpp"
#include "physics/ephemeris.hpp"
#include "physics/kepler_orbit.hpp"
#include "physics/massive_body.hpp"
#include "physics/solar_system.hpp"
#include "quantities/astronomy.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/si.hpp"

namespace principia {

using base::Bundle;
using base::not_null;
using base::OFStream;
using base::Status;
using geometry::Instant;
using geometry::Position;
using geometry::Sign;
using integrators::SymmetricLinearMultistepIntegrator;
using integrators::methods::Quinlan1999Order8A;
using numerics::Bisect;
using physics::Ephemeris;
using physics::KeplerianElements;
using physics::KeplerOrbit;
using physics::MassiveBody;
using physics::RelativeDegreesOfFreedom;
using physics::SolarSystem;
using quantities::Square;
using quantities::Sqrt;
using quantities::Time;
using quantities::astronomy::JulianYear;
using quantities::si::Day;
using quantities::si::Degree;
using quantities::si::Hour;
using quantities::si::Metre;
using quantities::si::Milli;
using quantities::si::Second;

namespace astronomy {

using Transits = std::vector<Instant>;
using TransitsByPlanet = std::map<std::string, Transits>;

// The random number generator used by the optimisation.

// The description of the characteristics of an individual, i.e., a
// configuration of the Trappist system.
class Genome {
 public:
  explicit Genome(std::vector<KeplerianElements<Trappist>> const& elements);

  std::vector<KeplerianElements<Trappist>> const& elements() const;

  void Mutate(std::mt19937_64& engine, int generation);

  static Genome OnePointCrossover(Genome const& g1,
                                  Genome const& g2,
                                  std::mt19937_64& engine);
  static Genome TwoPointCrossover(Genome const& g1,
                                  Genome const& g2,
                                  std::mt19937_64& engine);
  static Genome Blend(Genome const& g1,
                      Genome const& g2,
                      std::mt19937_64& engine);

 private:
  std::vector<KeplerianElements<Trappist>> elements_;
};

// A set of genomes which can reproduce based on their fitness.
class Population {
 public:
  Population(Genome const& luca,
             int const size,
             std::function<double(Genome const&)> compute_fitness,
             std::function<std::string(Genome const&)> residual_trace);

  void ComputeAllFitnesses();

  void BegetChildren();

  Genome best_genome() const;

 private:
  Genome const* Pick() const;

  std::function<double(Genome const&)> const compute_fitness_;
  std::function<std::string(Genome const&)> const residual_trace_;
  mutable std::mt19937_64 engine_;
  std::vector<Genome> current_;
  std::vector<Genome> next_;
  std::vector<double> fitnesses_;
  std::vector<double> cumulative_fitnesses_;

  int generation_ = 0;
  double best_fitness_ = 0.0;
  std::optional<Genome> best_genome_;
};

Genome::Genome(std::vector<KeplerianElements<Trappist>> const& elements)
    : elements_(elements) {}

std::vector<KeplerianElements<Trappist>> const& Genome::elements() const {
  return elements_;
}

void Genome::Mutate(std::mt19937_64& engine, int generation)  {
  for (auto& element : elements_) {
    element.asymptotic_true_anomaly = std::nullopt;
    element.turning_angle = std::nullopt;
    element.semimajor_axis = std::nullopt;
    element.specific_energy = std::nullopt;
    element.characteristic_energy = std::nullopt;
    element.mean_motion = std::nullopt;
    element.hyperbolic_mean_motion = std::nullopt;
    element.hyperbolic_excess_velocity = std::nullopt;
    element.semiminor_axis = std::nullopt;
    element.impact_parameter = std::nullopt;
    element.semilatus_rectum = std::nullopt;
    element.specific_angular_momentum = std::nullopt;
    element.periapsis_distance = std::nullopt;
    element.apoapsis_distance = std::nullopt;
    element.longitude_of_periapsis = std::nullopt;
    element.true_anomaly = std::nullopt;
    element.hyperbolic_mean_anomaly = std::nullopt;
    // The standard deviation of the distribution below has a strong effect on
    // the convergence of the algorithm: if it's too small we do not explore the
    // genomic space efficiently and it takes forever to find decent solutions;
    // if it's too large we explore the genomic space haphazardly and suffer
    // from deleterious mutations.
    double multiplicator =
        std::exp2(-2 - (generation % 200) / 20);
    if (generation == -1) multiplicator = 1;
    std::student_t_distribution<> distribution(1);
    *element.argument_of_periapsis +=
        distribution(engine) * 10 * Degree * multiplicator;
    *element.argument_of_periapsis =
        std::fmod(*element.argument_of_periapsis / quantities::si::Radian,
                  2 * π) *
        quantities::si::Radian;
    if (*element.argument_of_periapsis < 0 * quantities::si::Radian) {
      *element.argument_of_periapsis += 2 * π * quantities::si::Radian;
    }
    *element.mean_anomaly += distribution(engine) * 10 * Degree * multiplicator;
    *element.mean_anomaly =
        std::fmod(*element.mean_anomaly / quantities::si::Radian, 2 * π) *
        quantities::si::Radian;
    if (*element.mean_anomaly < 0 * quantities::si::Radian) {
      *element.mean_anomaly += 2 * π * quantities::si::Radian;
    }
    *element.period += distribution(engine) * 4 * Second * multiplicator;
    element.eccentricity =
        std::max(0.0,
                 std::min(*element.eccentricity +
                              distribution(engine) * 1e-3 * multiplicator,
                          0.2));
  }
}

Genome Genome::OnePointCrossover(Genome const& g1,
                                 Genome const& g2,
                                 std::mt19937_64& engine) {
  CHECK_EQ(g1.elements_.size(), g2.elements_.size());
  std::vector<KeplerianElements<Trappist>> new_elements;
  std::uniform_int_distribution<> order_distribution(0, 1);
  std::uniform_int_distribution<> split_distribution(0, g1.elements_.size());
  bool const reverse = order_distribution(engine) == 1;
  int const split = split_distribution(engine);<commit_msg>faster cooling, but slower for periods<commit_after>﻿
#include <random>

#include "astronomy/frames.hpp"
#include "base/bundle.hpp"
#include "base/file.hpp"
#include "base/not_null.hpp"
#include "base/status.hpp"
#include "geometry/named_quantities.hpp"
#include "geometry/sign.hpp"
#include "gtest/gtest.h"
#include "integrators/methods.hpp"
#include "integrators/symmetric_linear_multistep_integrator.hpp"
#include "mathematica/mathematica.hpp"
#include "numerics/root_finders.hpp"
#include "physics/degrees_of_freedom.hpp"
#include "physics/ephemeris.hpp"
#include "physics/kepler_orbit.hpp"
#include "physics/massive_body.hpp"
#include "physics/solar_system.hpp"
#include "quantities/astronomy.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/si.hpp"

namespace principia {

using base::Bundle;
using base::not_null;
using base::OFStream;
using base::Status;
using geometry::Instant;
using geometry::Position;
using geometry::Sign;
using integrators::SymmetricLinearMultistepIntegrator;
using integrators::methods::Quinlan1999Order8A;
using numerics::Bisect;
using physics::Ephemeris;
using physics::KeplerianElements;
using physics::KeplerOrbit;
using physics::MassiveBody;
using physics::RelativeDegreesOfFreedom;
using physics::SolarSystem;
using quantities::Square;
using quantities::Sqrt;
using quantities::Time;
using quantities::astronomy::JulianYear;
using quantities::si::Day;
using quantities::si::Degree;
using quantities::si::Hour;
using quantities::si::Metre;
using quantities::si::Milli;
using quantities::si::Second;

namespace astronomy {

using Transits = std::vector<Instant>;
using TransitsByPlanet = std::map<std::string, Transits>;

// The random number generator used by the optimisation.

// The description of the characteristics of an individual, i.e., a
// configuration of the Trappist system.
class Genome {
 public:
  explicit Genome(std::vector<KeplerianElements<Trappist>> const& elements);

  std::vector<KeplerianElements<Trappist>> const& elements() const;

  void Mutate(std::mt19937_64& engine, int generation);

  static Genome OnePointCrossover(Genome const& g1,
                                  Genome const& g2,
                                  std::mt19937_64& engine);
  static Genome TwoPointCrossover(Genome const& g1,
                                  Genome const& g2,
                                  std::mt19937_64& engine);
  static Genome Blend(Genome const& g1,
                      Genome const& g2,
                      std::mt19937_64& engine);

 private:
  std::vector<KeplerianElements<Trappist>> elements_;
};

// A set of genomes which can reproduce based on their fitness.
class Population {
 public:
  Population(Genome const& luca,
             int const size,
             std::function<double(Genome const&)> compute_fitness,
             std::function<std::string(Genome const&)> residual_trace);

  void ComputeAllFitnesses();

  void BegetChildren();

  Genome best_genome() const;

 private:
  Genome const* Pick() const;

  std::function<double(Genome const&)> const compute_fitness_;
  std::function<std::string(Genome const&)> const residual_trace_;
  mutable std::mt19937_64 engine_;
  std::vector<Genome> current_;
  std::vector<Genome> next_;
  std::vector<double> fitnesses_;
  std::vector<double> cumulative_fitnesses_;

  int generation_ = 0;
  double best_fitness_ = 0.0;
  std::optional<Genome> best_genome_;
};

Genome::Genome(std::vector<KeplerianElements<Trappist>> const& elements)
    : elements_(elements) {}

std::vector<KeplerianElements<Trappist>> const& Genome::elements() const {
  return elements_;
}

void Genome::Mutate(std::mt19937_64& engine, int generation)  {
  for (auto& element : elements_) {
    element.asymptotic_true_anomaly = std::nullopt;
    element.turning_angle = std::nullopt;
    element.semimajor_axis = std::nullopt;
    element.specific_energy = std::nullopt;
    element.characteristic_energy = std::nullopt;
    element.mean_motion = std::nullopt;
    element.hyperbolic_mean_motion = std::nullopt;
    element.hyperbolic_excess_velocity = std::nullopt;
    element.semiminor_axis = std::nullopt;
    element.impact_parameter = std::nullopt;
    element.semilatus_rectum = std::nullopt;
    element.specific_angular_momentum = std::nullopt;
    element.periapsis_distance = std::nullopt;
    element.apoapsis_distance = std::nullopt;
    element.longitude_of_periapsis = std::nullopt;
    element.true_anomaly = std::nullopt;
    element.hyperbolic_mean_anomaly = std::nullopt;
    // The standard deviation of the distribution below has a strong effect on
    // the convergence of the algorithm: if it's too small we do not explore the
    // genomic space efficiently and it takes forever to find decent solutions;
    // if it's too large we explore the genomic space haphazardly and suffer
    // from deleterious mutations.
    double multiplicator =
        std::exp2(-2 - (generation % 100) / 15);
    if (generation == -1) multiplicator = 1;
    std::student_t_distribution<> distribution(1);
    *element.argument_of_periapsis +=
        distribution(engine) * 10 * Degree * multiplicator;
    *element.argument_of_periapsis =
        std::fmod(*element.argument_of_periapsis / quantities::si::Radian,
                  2 * π) *
        quantities::si::Radian;
    if (*element.argument_of_periapsis < 0 * quantities::si::Radian) {
      *element.argument_of_periapsis += 2 * π * quantities::si::Radian;
    }
    *element.mean_anomaly += distribution(engine) * 10 * Degree * multiplicator;
    *element.mean_anomaly =
        std::fmod(*element.mean_anomaly / quantities::si::Radian, 2 * π) *
        quantities::si::Radian;
    if (*element.mean_anomaly < 0 * quantities::si::Radian) {
      *element.mean_anomaly += 2 * π * quantities::si::Radian;
    }
    *element.period += distribution(engine) * 4 * Second * std::sqrt(multiplicator);
    element.eccentricity =
        std::max(0.0,
                 std::min(*element.eccentricity +
                              distribution(engine) * 1e-3 * multiplicator,
                          0.2));
  }
}

Genome Genome::OnePointCrossover(Genome const& g1,
                                 Genome const& g2,
                                 std::mt19937_64& engine) {
  CHECK_EQ(g1.elements_.size(), g2.elements_.size());
  std::vector<KeplerianElements<Trappist>> new_elements;
  std::uniform_int_distribution<> order_distribution(0, 1);
  std::uniform_int_distribution<> split_distribution(0, g1.elements_.size());
  bool const reverse = order_distribution(engine) == 1;
  int const split = split_distribution(engine);
{'fraction_non_alphanumeric': 0.11684700931770364, 'fraction_numerical': 0.014878268710550046, 'mean_word_length': 4.022264150943396, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5537148', 'n_tokens_mistral': 4322, 'n_tokens_neox': 4045, 'n_words': 886}