TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
79 / 79
100.0%
Branches:
87 / 88
98.9%
Functions:
9 / 9
100.0%
random/samplers.cpp
Line Branch Exec Source
1 #include "terraingen/random/samplers.h"
2
3 #include "terraingen/util/checks.h"
4
5 #include <algorithm>
6 #include <cmath>
7 #include <iterator>
8 #include <stdexcept>
9 #include <utility>
10
11 namespace terraingen {
12
13 namespace {
14
15 /// Poisson count by multiplying uniforms until the product falls below exp(-`mean`); for small
16 /// means.
17 5581567 std::uint64_t PoissonMultiplication(RandomStream& random, double mean) {
18 5581567 const double limit = std::exp(-mean);
19 5581567 std::uint64_t count = 0;
20 5581567 double product = random.NextDouble();
21
2/2
✓ Branch 0 taken 13886660 times.
✓ Branch 1 taken 5581567 times.
 19468227 while (product > limit) {
22 13886660 ++count;
23 13886660 product *= random.NextDouble();
24 }
25 5581567 return count;
26 }
27
28 /// Poisson count by Hormann's transformed rejection with squeeze (PTRS); for means of 10 or more.
29 314801 std::uint64_t PoissonPtrs(RandomStream& random, double mean) {
30 314801 const double sqrtMean = std::sqrt(mean);
31 314801 const double logMean = std::log(mean);
32 314801 const double b = 0.931 + 2.53 * sqrtMean;
33 314801 const double a = -0.059 + 0.02483 * b;
34 314801 const double inverseAlpha = 1.1239 + 1.1328 / (b - 3.4);
35 314801 const double vr = 0.9277 - 3.6224 / (b - 2.0);
36 while (true) {
37 385365 const double u = random.NextDouble() - 0.5;
38 385365 const double v = random.NextDouble();
39 385365 const double us = 0.5 - std::fabs(u);
40 385365 const double k = std::floor((2.0 * a / us + b) * u + mean + 0.43);
41
4/4
✓ Branch 0 taken 333620 times.
✓ Branch 1 taken 51745 times.
✓ Branch 2 taken 215373 times.
✓ Branch 3 taken 118247 times.
 385365 if (us >= 0.07 && v <= vr) {
42 215373 return static_cast<std::uint64_t>(k);
43 }
44
6/6
✓ Branch 0 taken 165403 times.
✓ Branch 1 taken 4589 times.
✓ Branch 2 taken 5696 times.
✓ Branch 3 taken 159707 times.
✓ Branch 4 taken 5668 times.
✓ Branch 5 taken 28 times.
 169992 if (k < 0.0 || (us < 0.013 && v > us)) {
45 10257 continue;
46 }
47 159735 if (std::log(v) + std::log(inverseAlpha) - std::log(a / (us * us) + b) <=
48
2/2
✓ Branch 0 taken 99428 times.
✓ Branch 1 taken 60307 times.
 159735 -mean + k * logMean - std::lgamma(k + 1.0)) {
49 99428 return static_cast<std::uint64_t>(k);
50 }
51 70564 }
52 }
53
54 } // namespace
55
56 5896394 std::uint64_t SamplePoisson(RandomStream& random, double mean) {
57
8/8
✓ Branch 1 taken 5896392 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 5896390 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 5896389 times.
✓ Branch 7 taken 5 times.
✓ Branch 8 taken 5896389 times.
 5896394 if (!(std::isfinite(mean) && mean >= 0.0 && mean <= 1e9)) {
58
1/1
✓ Branch 2 taken 5 times.
 5 throw std::invalid_argument("Poisson mean must be finite and in [0, 1e9]");
59 }
60
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 5896368 times.
 5896389 if (mean == 0.0) {
61 21 return 0;
62 }
63
2/2
✓ Branch 0 taken 5581567 times.
✓ Branch 1 taken 314801 times.
 5896368 return mean < 10.0 ? PoissonMultiplication(random, mean) : PoissonPtrs(random, mean);
64 }
65
66 72007 double SampleTruncatedPowerLaw(RandomStream& random, double min, double max, double exponent) {
67
7/8
✓ Branch 1 taken 72004 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 72002 times.
✓ Branch 4 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 72002 times.
✓ Branch 8 taken 5 times.
✓ Branch 9 taken 72002 times.
 72007 if (!(IsPositive(min) && max > min && !std::isnan(max))) {
68
1/1
✓ Branch 2 taken 5 times.
 5 throw std::invalid_argument("Power law bounds must satisfy 0 < min < max");
69 }
70
1/1
✓ Branch 2 taken 72000 times.
 72002 RequirePositive(exponent, "Power law exponent");
71 72000 const double low = std::pow(min, -exponent);
72
2/2
✓ Branch 1 taken 1000 times.
✓ Branch 2 taken 71000 times.
 72000 const double high = std::isinf(max) ? 0.0 : std::pow(max, -exponent);
73
1/1
✓ Branch 1 taken 72000 times.
 72000 const double u = random.NextDouble();
74 72000 const double x = std::pow(low - u * (low - high), -1.0 / exponent);
75
1/1
✓ Branch 1 taken 72000 times.
 144000 return std::clamp(x, min, max);
76 }
77
78 50007 double SampleExponential(RandomStream& random, double rate) {
79
1/1
✓ Branch 2 taken 50004 times.
 50007 RequirePositive(rate, "Exponential rate");
80 50004 return -std::log(random.NextDoubleOpenBelow()) / rate;
81 }
82
83 50006 double SampleNormal(RandomStream& random, double mean, double sigma) {
84
6/6
✓ Branch 1 taken 50005 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 50002 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 50002 times.
 50006 if (!(std::isfinite(mean) && IsNonNegative(sigma))) {
85
1/1
✓ Branch 2 taken 4 times.
 4 throw std::invalid_argument("Normal mean must be finite and sigma finite and not negative");
86 }
87 50002 constexpr double kTwoPi = 6.283185307179586477;
88 50002 const double radius = std::sqrt(-2.0 * std::log(random.NextDoubleOpenBelow()));
89 50002 return mean + sigma * radius * std::cos(kTwoPi * random.NextDouble());
90 }
91
92 14 TabulatedDistribution::TabulatedDistribution(std::vector<double> values,
93 14 std::vector<double> cumulative)
94 14 : values_(std::move(values)), cumulative_(std::move(cumulative)) {
95
6/6
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 12 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 12 times.
 14 if (values_.size() < 2 || values_.size() != cumulative_.size()) {
96 throw std::invalid_argument(
97
1/1
✓ Branch 2 taken 2 times.
 2 "A tabulated distribution needs at least two values and one CDF value per value");
98 }
99
2/2
✓ Branch 1 taken 32 times.
✓ Branch 2 taken 8 times.
 40 for (std::size_t i = 0; i < values_.size(); ++i) {
100
6/6
✓ Branch 2 taken 31 times.
✓ Branch 3 taken 1 time.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 30 times.
✓ Branch 8 taken 2 times.
✓ Branch 9 taken 30 times.
 32 if (!std::isfinite(values_[i]) || !std::isfinite(cumulative_[i])) {
101
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Tabulated values and CDF must be finite");
102 }
103
8/8
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 12 times.
✓ Branch 4 taken 17 times.
✓ Branch 5 taken 1 time.
✓ Branch 8 taken 1 time.
✓ Branch 9 taken 16 times.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 28 times.
 30 if (i > 0 && (values_[i] <= values_[i - 1] || cumulative_[i] < cumulative_[i - 1])) {
104 throw std::invalid_argument(
105
1/1
✓ Branch 2 taken 2 times.
 2 "Tabulated values must increase strictly and the CDF must not decrease");
106 }
107 }
108
6/6
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 5 times.
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 5 times.
 8 if (cumulative_.front() != 0.0 || cumulative_.back() != 1.0) {
109
1/1
✓ Branch 2 taken 3 times.
 3 throw std::invalid_argument("Tabulated CDF must start at 0 and end at 1");
110 }
111 23 }
112
113 70009 double TabulatedDistribution::Cdf(double x) const {
114
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 70006 times.
 70009 if (x <= values_.front()) {
115 3 return 0.0;
116 }
117
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 70003 times.
 70006 if (x >= values_.back()) {
118 3 return 1.0;
119 }
120
1/1
✓ Branch 3 taken 70003 times.
 70003 const auto upper = std::upper_bound(values_.begin(), values_.end(), x);
121 70003 const auto i = static_cast<std::size_t>(std::distance(values_.begin(), upper)) - 1;
122 70003 const double t = (x - values_[i]) / (values_[i + 1] - values_[i]);
123 70003 return cumulative_[i] + t * (cumulative_[i + 1] - cumulative_[i]);
124 }
125
126 75001 double TabulatedDistribution::Sample(RandomStream& random) const {
127
1/1
✓ Branch 1 taken 75001 times.
 75001 const double u = random.NextDouble();
128
1/1
✓ Branch 3 taken 75001 times.
 75001 const auto upper = std::upper_bound(cumulative_.begin(), cumulative_.end(), u);
129 75001 const auto i = static_cast<std::size_t>(std::distance(cumulative_.begin(), upper)) - 1;
130 75001 const double t = (u - cumulative_[i]) / (cumulative_[i + 1] - cumulative_[i]);
131 150002 return values_[i] + t * (values_[i + 1] - values_[i]);
132 }
133
134 } // namespace terraingen
135