TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
84 / 92
91.3%
Branches:
75 / 89
84.3%
Functions:
11 / 11
100.0%
distribution/size_frequency.cpp
Line Branch Exec Source
1 #include "terraingen/distribution/size_frequency.h"
2
3 #include "terraingen/util/checks.h"
4
5 #include <algorithm>
6 #include <cmath>
7 #include <stdexcept>
8 #include <string>
9
10 namespace terraingen {
11
12 namespace {
13
14 102604 void CheckRange(const DiameterRange& range) {
15
8/10
✓ Branch 1 taken 102604 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 102603 times.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 102601 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 102601 times.
✓ Branch 9 taken 3 times.
✓ Branch 10 taken 102601 times.
 205205 if (!(std::isfinite(range.min) && range.min >= 0.0 && range.max > range.min) ||
16 102601 std::isnan(range.max)) {
17
1/1
✓ Branch 2 taken 3 times.
 3 throw std::invalid_argument("Diameter range must satisfy 0 <= min < max");
18 }
19 102601 }
20
21 5003968 void CheckInterval(double dMin, double dMax) {
22
7/8
✓ Branch 1 taken 5003967 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 5003964 times.
✓ Branch 4 taken 3 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5003964 times.
✓ Branch 8 taken 4 times.
✓ Branch 9 taken 5003964 times.
 5003968 if (!(IsPositive(dMin) && dMax > dMin) || std::isnan(dMax)) {
23
1/1
✓ Branch 2 taken 4 times.
 4 throw std::invalid_argument("Diameter interval must satisfy 0 < dMin < dMax");
24 }
25 5003964 }
26
27 } // namespace
28
29 102604 SizeFrequencyDistribution::SizeFrequencyDistribution(SfdKind kind, DiameterRange range)
30 102604 : kind_(kind), range_(range) {
31 102604 CheckRange(range_);
32 102601 }
33
34 21749909 void SizeFrequencyDistribution::CheckDiameter(double diameter) const {
35
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 21749905 times.
 21749909 if (!IsPositive(diameter)) {
36
1/1
✓ Branch 2 taken 4 times.
 4 throw std::domain_error("Diameter must be positive and finite");
37 }
38
6/6
✓ Branch 0 taken 21618133 times.
✓ Branch 1 taken 131772 times.
✓ Branch 3 taken 19 times.
✓ Branch 4 taken 21618114 times.
✓ Branch 5 taken 19 times.
✓ Branch 6 taken 21749886 times.
 21749905 if (!extrapolate_ && !range_.Contains(diameter)) {
39
2/2
✓ Branch 1 taken 19 times.
✓ Branch 4 taken 19 times.
 38 throw std::domain_error("Diameter " + std::to_string(diameter) +
40
2/2
✓ Branch 2 taken 19 times.
✓ Branch 5 taken 19 times.
 57 " m is outside the distribution's valid range");
41 }
42 21749886 }
43
44 14184236 double SizeFrequencyDistribution::Cumulative(double diameter) const {
45 14184236 CheckDiameter(diameter);
46 14184217 return CumulativeAt(diameter);
47 }
48
49 7565673 double SizeFrequencyDistribution::Incremental(double diameter) const {
50 7565673 CheckDiameter(diameter);
51 7565669 return IncrementalAt(diameter);
52 }
53
54 2 double SizeFrequencyDistribution::IncrementalAt(double diameter) const {
55 2 constexpr double kStep = 1e-5;
56 2 const double low = diameter * std::exp(-kStep);
57 2 const double high = diameter * std::exp(kStep);
58 2 return (CumulativeAt(low) - CumulativeAt(high)) / (high - low);
59 }
60
61 1893683 double SizeFrequencyDistribution::InvertCumulative(double target, double dMin, double dMax) const {
62 1893683 double low = std::log(dMin);
63
2/2
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 1893663 times.
 1893683 double high = std::isinf(dMax) ? low + 1.0 : std::log(dMax);
64
6/6
✓ Branch 1 taken 39 times.
✓ Branch 2 taken 1893663 times.
✓ Branch 4 taken 19 times.
✓ Branch 5 taken 20 times.
✓ Branch 6 taken 19 times.
✓ Branch 7 taken 1893683 times.
 1893702 while (std::isinf(dMax) && CumulativeAt(std::exp(high)) > target) {
65 19 high += 1.0;
66 }
67
1/2
✓ Branch 0 taken 1893683 times.
✗ Branch 1 not taken.
 1893683 if (target > 0.0) {
68 1893683 const double logTarget = std::log(target);
69 10026747 const auto excess = [&](double u) {
70 10026747 return std::log(CumulativeAt(std::exp(u))) - logTarget;
71 1893683 };
72
1/1
✓ Branch 1 taken 1893683 times.
 1893683 double gLow = excess(low);
73
1/1
✓ Branch 1 taken 1893683 times.
 1893683 double gHigh = excess(high);
74
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1893683 times.
 1893683 if (!(gLow > 0.0)) {
75 1893683 return std::exp(low);
76 }
77
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1893683 times.
 1893683 if (!(gHigh < 0.0)) {
78 ✗ return std::exp(high);
79 }
80 1893683 int side = 0;
81 1893683 double u = low;
82
1/2
✓ Branch 0 taken 6239381 times.
✗ Branch 1 not taken.
 6239381 for (int i = 0; i < 100; ++i) {
83 6239381 const double previous = u;
84 6239381 u = (low * gHigh - high * gLow) / (gHigh - gLow);
85
4/4
✓ Branch 0 taken 6030982 times.
✓ Branch 1 taken 208399 times.
✓ Branch 2 taken 229451 times.
✓ Branch 3 taken 5801531 times.
 6239381 if (!(u > low && u < high)) {
86 437850 u = 0.5 * (low + high);
87 }
88
1/1
✓ Branch 1 taken 6239381 times.
 6239381 const double g = excess(u);
89
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 6239381 times.
 6239381 if (std::isnan(g)) {
90 ✗ break;
91 }
92
4/4
✓ Branch 0 taken 4924593 times.
✓ Branch 1 taken 1314788 times.
✓ Branch 2 taken 578895 times.
✓ Branch 3 taken 4345698 times.
 6239381 if (g == 0.0 || std::fabs(u - previous) <= 1e-14 * std::fmax(1.0, std::fabs(u))) {
93 1893683 return std::exp(u);
94 }
95
2/2
✓ Branch 0 taken 2226150 times.
✓ Branch 1 taken 2119548 times.
 4345698 if (g > 0.0) {
96 2226150 low = u;
97 2226150 gLow = g;
98
2/2
✓ Branch 0 taken 763612 times.
✓ Branch 1 taken 1462538 times.
 2226150 if (side == 1) {
99 763612 gHigh *= 0.5;
100 }
101 2226150 side = 1;
102 } else {
103 2119548 high = u;
104 2119548 gHigh = g;
105
2/2
✓ Branch 0 taken 712536 times.
✓ Branch 1 taken 1407012 times.
 2119548 if (side == -1) {
106 712536 gLow *= 0.5;
107 }
108 2119548 side = -1;
109 }
110 }
111 }
112
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 ✗ for (int i = 0; i < 200 && high - low > 1e-15 * std::fmax(1.0, std::fabs(high)); ++i) {
113 ✗ const double middle = 0.5 * (low + high);
114
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 ✗ if (CumulativeAt(std::exp(middle)) > target) {
115 ✗ low = middle;
116 } else {
117 ✗ high = middle;
118 }
119 }
120 ✗ return std::exp(0.5 * (low + high));
121 }
122
123 5003968 double SizeFrequencyDistribution::ExpectedCount(double area, double dMin, double dMax) const {
124 5003968 CheckInterval(dMin, dMax);
125
1/1
✓ Branch 2 taken 5003962 times.
 5003964 RequireNonNegative(area, "Area");
126 5003962 double upper = 0.0;
127
2/2
✓ Branch 1 taken 40183 times.
✓ Branch 2 taken 4963779 times.
 5003962 if (std::isinf(dMax)) {
128
6/6
✓ Branch 0 taken 40181 times.
✓ Branch 1 taken 2 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 40180 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 40182 times.
 40183 if (!extrapolate_ && !std::isinf(range_.max)) {
129
1/1
✓ Branch 2 taken 1 time.
 1 throw std::domain_error("An infinite dMax needs an unbounded valid range");
130 }
131 } else {
132 4963779 upper = Cumulative(dMax);
133 }
134 5003960 return area * (Cumulative(dMin) - upper);
135 }
136
137 3561571 double SizeFrequencyDistribution::SampleDiameter(RandomStream& random, double dMin,
138 double dMax) const {
139 3561571 const double count = ExpectedCount(1.0, dMin, dMax);
140
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 3561568 times.
 3561570 if (!(count > 0.0)) {
141
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("The distribution has no mass in the diameter interval");
142 }
143 3561568 const double top = Cumulative(dMin);
144 3561568 const double target = top - random.NextDouble() * count;
145
2/2
✓ Branch 1 taken 3561568 times.
✓ Branch 4 taken 3561568 times.
 3561568 return std::clamp(InvertCumulative(target, dMin, dMax), dMin, dMax);
146 }
147
148 } // namespace terraingen
149