distribution/tabulated.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/distribution/tabulated.h" | ||
| 2 | |||
| 3 | #include "terraingen/util/checks.h" | ||
| 4 | |||
| 5 | #include <algorithm> | ||
| 6 | #include <cmath> | ||
| 7 | #include <stdexcept> | ||
| 8 | #include <utility> | ||
| 9 | |||
| 10 | namespace terraingen { | ||
| 11 | |||
| 12 | namespace { | ||
| 13 | |||
| 14 | /// The table's diameter range in metres, after checking the table and units. | ||
| 15 | 70 | DiameterRange TableRange(const std::vector<double>& diameters, | |
| 16 | const std::vector<double>& cumulative, double diameterUnit, | ||
| 17 | double areaUnit) { | ||
| 18 | 70 | if (diameters.size() < 2 || diameters.size() != cumulative.size()) { | |
| 19 | 5 | throw std::invalid_argument("A table needs at least two diameters with one value each"); | |
| 20 | } | ||
| 21 | 65 | if (!IsPositive(diameterUnit) || !IsPositive(areaUnit)) { | |
| 22 | 2 | throw std::invalid_argument("Diameter and area units must be positive and finite"); | |
| 23 | } | ||
| 24 | 640 | for (std::size_t i = 0; i < diameters.size(); ++i) { | |
| 25 | 584 | if (!IsPositive(diameters[i]) || !IsPositive(cumulative[i])) { | |
| 26 | 4 | throw std::invalid_argument("Table diameters and values must be positive and finite"); | |
| 27 | } | ||
| 28 | 580 | if (i > 0 && !(diameters[i] > diameters[i - 1])) { | |
| 29 | 2 | throw std::invalid_argument("Table diameters must strictly increase"); | |
| 30 | } | ||
| 31 | 578 | if (i > 0 && cumulative[i] > cumulative[i - 1]) { | |
| 32 | 1 | throw std::invalid_argument("A cumulative distribution cannot increase with diameter"); | |
| 33 | } | ||
| 34 | } | ||
| 35 | 56 | return DiameterRange{diameters.front() * diameterUnit, diameters.back() * diameterUnit}; | |
| 36 | } | ||
| 37 | |||
| 38 | } // namespace | ||
| 39 | |||
| 40 | 70 | TabulatedSfd::TabulatedSfd(std::vector<double> diameters, std::vector<double> cumulative, | |
| 41 | 70 | SfdKind kind, double diameterUnit, double areaUnit) | |
| 42 | : SizeFrequencyDistribution(kind, TableRange(diameters, cumulative, diameterUnit, areaUnit)), | ||
| 43 | 70 | diameters_(std::move(diameters)), cumulative_(std::move(cumulative)) { | |
| 44 | 628 | for (std::size_t i = 0; i < diameters_.size(); ++i) { | |
| 45 | 572 | diameters_[i] *= diameterUnit; | |
| 46 | 572 | cumulative_[i] /= areaUnit; | |
| 47 | 572 | if (i > 0) { | |
| 48 | 1032 | slopes_.push_back(std::log(cumulative_[i] / cumulative_[i - 1]) / | |
| 49 | 516 | std::log(diameters_[i] / diameters_[i - 1])); | |
| 50 | } | ||
| 51 | } | ||
| 52 | 56 | } | |
| 53 | |||
| 54 | 1025074 | std::size_t TabulatedSfd::IntervalAt(double diameter) const { | |
| 55 | 1025074 | const auto upper = std::upper_bound(diameters_.begin(), diameters_.end(), diameter); | |
| 56 | 1025074 | const auto index = static_cast<std::size_t>(upper - diameters_.begin()); | |
| 57 | 1025074 | return std::clamp<std::size_t>(index, 1, diameters_.size() - 1) - 1; | |
| 58 | } | ||
| 59 | |||
| 60 | 917903 | double TabulatedSfd::CumulativeAt(double diameter) const { | |
| 61 | 917903 | const std::size_t i = IntervalAt(diameter); | |
| 62 | 917903 | return cumulative_[i] * std::pow(diameter / diameters_[i], slopes_[i]); | |
| 63 | } | ||
| 64 | |||
| 65 | 107171 | double TabulatedSfd::IncrementalAt(double diameter) const { | |
| 66 | 107171 | return -slopes_[IntervalAt(diameter)] * CumulativeAt(diameter) / diameter; | |
| 67 | } | ||
| 68 | |||
| 69 | 42988 | double TabulatedSfd::InvertCumulative(double target, double /*dMin*/, double /*dMax*/) const { | |
| 70 | 42988 | const auto above = std::partition_point(cumulative_.begin(), cumulative_.end(), | |
| 71 | 121017 | [target](double value) { return value > target; }); | |
| 72 | 42988 | const auto index = static_cast<std::size_t>(above - cumulative_.begin()); | |
| 73 | 42988 | const std::size_t i = std::clamp<std::size_t>(index, 1, cumulative_.size() - 1) - 1; | |
| 74 | 42988 | if (slopes_[i] == 0.0) { | |
| 75 | 2 | return diameters_[i]; | |
| 76 | } | ||
| 77 | 42986 | return diameters_[i] * std::pow(target / cumulative_[i], 1.0 / slopes_[i]); | |
| 78 | } | ||
| 79 | |||
| 80 | } // namespace terraingen | ||
| 81 |