Skip to content

Commit

Permalink
Merge pull request #14 from darioizzo/main
Browse files Browse the repository at this point in the history 
round of clang format and throws removed
  • Loading branch information
@darioizzo
darioizzo authored Oct 4, 2023
2 parents 91707cc + e171a4a commit 2c02d90
Show file tree
Hide file tree
Showing 48 changed files with 18,537 additions and 16,619 deletions.
142 changes: 71 additions & 71 deletions benchmark/convert_anomalies_benchmark.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,86 +24,86 @@ using std::chrono::microseconds;
// In this benchmark we test the speed and accuracy of the Kepler's equation
// solvers

void perform_test_speed(double min_ecc, double max_ecc, unsigned N) {
//
// Engines
//
// NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
std::mt19937 rng_engine(122012203u);
//
// Distributions
//
std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
std::uniform_real_distribution<double> M_d(-1e8, 1e8);
void perform_test_speed(double min_ecc, double max_ecc, unsigned N)
{
//
// Engines
//
// NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
std::mt19937 rng_engine(122012203u);
//
// Distributions
//
std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
std::uniform_real_distribution<double> M_d(-1e8, 1e8);

// We generate the random dataset
std::vector<double> eccenricities(N);
std::vector<double> mean_anomalies(N);
// We generate the random dataset
std::vector<double> eccenricities(N);
std::vector<double> mean_anomalies(N);

for (auto i = 0u; i < N; ++i) {
mean_anomalies[i] = M_d(rng_engine);
eccenricities[i] = ecc_d(rng_engine);
}
for (auto i = 0u; i < N; ++i) {
mean_anomalies[i] = M_d(rng_engine);
eccenricities[i] = ecc_d(rng_engine);
}

// We log progress
fmt::print("{:.2f} min_ecc, {:.2f} max_ecc, on {} data points: ", min_ecc,
max_ecc, N);
// We log progress
fmt::print("{:.2f} min_ecc, {:.2f} max_ecc, on {} data points: ", min_ecc, max_ecc, N);

auto start = high_resolution_clock::now();
for (auto i = 0u; i < N; ++i) {
m2e(mean_anomalies[i], eccenricities[i]);
}
auto stop = high_resolution_clock::now();
auto duration = duration_cast<microseconds>(stop - start);
fmt::print("{:.3f}s\n", (static_cast<double>(duration.count()) / 1e6));
auto start = high_resolution_clock::now();
for (auto i = 0u; i < N; ++i) {
m2e(mean_anomalies[i], eccenricities[i]);
}
auto stop = high_resolution_clock::now();
auto duration = duration_cast<microseconds>(stop - start);
fmt::print("{:.3f}s\n", (static_cast<double>(duration.count()) / 1e6));
}

void perform_test_accuracy(double min_ecc, double max_ecc, unsigned N) {
//
// Engines
//
// NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
std::mt19937 rng_engine(122012203u);
//
// Distributions
//
std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
std::uniform_real_distribution<double> M_d(-1e8, 1e8);
void perform_test_accuracy(double min_ecc, double max_ecc, unsigned N)
{
//
// Engines
//
// NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
std::mt19937 rng_engine(122012203u);
//
// Distributions
//
std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
std::uniform_real_distribution<double> M_d(-1e8, 1e8);

// We generate the random dataset
std::vector<double> eccenricities(N);
std::vector<double> mean_anomalies(N);
// We generate the random dataset
std::vector<double> eccenricities(N);
std::vector<double> mean_anomalies(N);

for (auto i = 0u; i < N; ++i) {
mean_anomalies[i] = M_d(rng_engine);
eccenricities[i] = ecc_d(rng_engine);
}
for (auto i = 0u; i < N; ++i) {
mean_anomalies[i] = M_d(rng_engine);
eccenricities[i] = ecc_d(rng_engine);
}

// We log progress
fmt::print("{:.2f} min_ecc, {:.2f} max_ecc, on {} data points: ", min_ecc,
max_ecc, N);
std::vector<double> err(N);
for (auto i = 0u; i < N; ++i) {
auto res = e2m(m2e(mean_anomalies[i], eccenricities[i]), eccenricities[i]);
// error is arbitrarily: (|sinM-sinMtrue| +|cosM-cosMtrue|)/2
err[i] = (std::abs(std::sin(res) - std::sin(mean_anomalies[i])) +
std::abs(std::cos(res) - std::cos(mean_anomalies[i]))) /
2.;
}
auto max_it = max_element(std::begin(err), std::end(err));
auto min_it = min_element(std::begin(err), std::end(err));
auto avg = std::accumulate(err.begin(), err.end(), 0.0) /
static_cast<double>(err.size());
fmt::print("{:.3e} avg, {:.3e} min, {:.3e} max\n", avg, *min_it, *max_it);
// We log progress
fmt::print("{:.2f} min_ecc, {:.2f} max_ecc, on {} data points: ", min_ecc, max_ecc, N);
std::vector<double> err(N);
for (auto i = 0u; i < N; ++i) {
auto res = e2m(m2e(mean_anomalies[i], eccenricities[i]), eccenricities[i]);
// error is arbitrarily: (|sinM-sinMtrue| +|cosM-cosMtrue|)/2
err[i] = (std::abs(std::sin(res) - std::sin(mean_anomalies[i]))
+ std::abs(std::cos(res) - std::cos(mean_anomalies[i])))
/ 2.;
}
auto max_it = max_element(std::begin(err), std::end(err));
auto min_it = min_element(std::begin(err), std::end(err));
auto avg = std::accumulate(err.begin(), err.end(), 0.0) / static_cast<double>(err.size());
fmt::print("{:.3e} avg, {:.3e} min, {:.3e} max\n", avg, *min_it, *max_it);
}

int main() {
fmt::print("\nComputes speed at different eccentricity ranges:\n");
perform_test_speed(0, 0.5, 1000000);
perform_test_speed(0.5, 0.9, 1000000);
perform_test_speed(0.9, 0.99, 1000000);
fmt::print("\nComputes error at different eccentricity ranges:\n");
perform_test_accuracy(0, 0.5, 100000);
perform_test_accuracy(0.5, 0.9, 100000);
perform_test_accuracy(0.9, 0.99, 100000);
int main()
{
fmt::print("\nComputes speed at different eccentricity ranges:\n");
perform_test_speed(0, 0.5, 1000000);
perform_test_speed(0.5, 0.9, 1000000);
perform_test_speed(0.9, 0.99, 1000000);
fmt::print("\nComputes error at different eccentricity ranges:\n");
perform_test_accuracy(0, 0.5, 100000);
perform_test_accuracy(0.5, 0.9, 100000);
perform_test_accuracy(0.9, 0.99, 100000);
}
104 changes: 52 additions & 52 deletions benchmark/lambert_problem_benchmark.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,61 +41,61 @@ using std::chrono::duration_cast;
using std::chrono::high_resolution_clock;
using std::chrono::microseconds;

int main() {
// Number of trials
const unsigned trials = 50000u;
int main()
{
// Number of trials
const unsigned trials = 50000u;

std::array<std::array<double, 3>, trials> r1s{}, r2s{};
std::array<double, trials> tof{};
std::array<bool, trials> cw{};
std::array<double, trials> mu{};
std::array<std::array<double, 3>, trials> r1s{}, r2s{};
std::array<double, trials> tof{};
std::array<bool, trials> cw{};
std::array<double, trials> mu{};

// NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
std::mt19937 rng_engine(122012203u);
std::uniform_int_distribution<unsigned> cw_d(0, 1);
std::uniform_real_distribution<double> r_d(-2, 2);
std::uniform_real_distribution<double> tof_d(2., 40.);
std::uniform_real_distribution<double> mu_d(0.9, 1.1);
unsigned revs_max = 20u;
unsigned long count = 0lu;
// NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
std::mt19937 rng_engine(122012203u);
std::uniform_int_distribution<unsigned> cw_d(0, 1);
std::uniform_real_distribution<double> r_d(-2, 2);
std::uniform_real_distribution<double> tof_d(2., 40.);
std::uniform_real_distribution<double> mu_d(0.9, 1.1);
unsigned revs_max = 20u;
unsigned long count = 0lu;

for (auto i = 0u; i < trials; ++i) {
// 1 - generate a random problem geometry
r1s[i][0] = r_d(rng_engine);
r1s[i][1] = r_d(rng_engine);
r1s[i][2] = r_d(rng_engine);
r2s[i][0] = r_d(rng_engine);
r2s[i][1] = r_d(rng_engine);
r2s[i][2] = r_d(rng_engine);
tof[i] = tof_d(rng_engine);
cw[i] = static_cast<bool>(cw_d(rng_engine));
mu[i] = mu_d(rng_engine);
}
for (auto i = 0u; i < trials; ++i) {
// 1 - generate a random problem geometry
r1s[i][0] = r_d(rng_engine);
r1s[i][1] = r_d(rng_engine);
r1s[i][2] = r_d(rng_engine);
r2s[i][0] = r_d(rng_engine);
r2s[i][1] = r_d(rng_engine);
r2s[i][2] = r_d(rng_engine);
tof[i] = tof_d(rng_engine);
cw[i] = static_cast<bool>(cw_d(rng_engine));
mu[i] = mu_d(rng_engine);
}

auto start = high_resolution_clock::now();
for (auto i = 0u; i < trials; ++i) {
// 2 - Solve the lambert problem
kep3::lambert_problem lp(r1s[i], r2s[i], tof[i], mu[i], cw[i], revs_max);
count = count + lp.get_v1().size();
}
auto stop = high_resolution_clock::now();
auto duration = duration_cast<microseconds>(stop - start);
fmt::print("Lambert:\n{} solutions computed in {:.3f}s\n", count,
(static_cast<double>(duration.count()) / 1e6));
fmt::print("Projected number of solutions per second: {}\n",
static_cast<double>(count) / ((static_cast<double>(duration.count()) / 1e6)));
auto start = high_resolution_clock::now();
for (auto i = 0u; i < trials; ++i) {
// 2 - Solve the lambert problem
kep3::lambert_problem lp(r1s[i], r2s[i], tof[i], mu[i], cw[i], revs_max);
count = count + lp.get_v1().size();
}
auto stop = high_resolution_clock::now();
auto duration = duration_cast<microseconds>(stop - start);
fmt::print("Lambert:\n{} solutions computed in {:.3f}s\n", count, (static_cast<double>(duration.count()) / 1e6));
fmt::print("Projected number of solutions per second: {}\n",
static_cast<double>(count) / ((static_cast<double>(duration.count()) / 1e6)));

count = 0; //reset counter
start = high_resolution_clock::now();
for (auto i = 0u; i < trials; ++i) {
// 3 - Solve the lambert problem for the singe rev case
kep3::lambert_problem lp(r1s[i], r2s[i], tof[i], mu[i], cw[i], 0u);
count += 1u;
}
stop = high_resolution_clock::now();
duration = duration_cast<microseconds>(stop - start);
fmt::print("\nLambert (0 revs only):\n{} solutions computed in {:.3f}s\n", count,
(static_cast<double>(duration.count()) / 1e6));
fmt::print("Projected number of solutions per second: {}\n",
static_cast<double>(count) / ((static_cast<double>(duration.count()) / 1e6)));
count = 0; // reset counter
start = high_resolution_clock::now();
for (auto i = 0u; i < trials; ++i) {
// 3 - Solve the lambert problem for the singe rev case
kep3::lambert_problem lp(r1s[i], r2s[i], tof[i], mu[i], cw[i], 0u);
count += 1u;
}
stop = high_resolution_clock::now();
duration = duration_cast<microseconds>(stop - start);
fmt::print("\nLambert (0 revs only):\n{} solutions computed in {:.3f}s\n", count,
(static_cast<double>(duration.count()) / 1e6));
fmt::print("Projected number of solutions per second: {}\n",
static_cast<double>(count) / ((static_cast<double>(duration.count()) / 1e6)));
}
Loading

0 comments on commit 2c02d90

Please sign in to comment.