diff --git a/nvbench/detail/measure_cold.cu b/nvbench/detail/measure_cold.cu index ed22854..eb63d40 100644 --- a/nvbench/detail/measure_cold.cu +++ b/nvbench/detail/measure_cold.cu @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -33,7 +34,6 @@ #include #include #include -#include #include #include @@ -201,6 +201,25 @@ bool measure_cold_base::is_finished() void measure_cold_base::run_trials_epilogue() { m_walltime_timer.stop(); } +void measure_cold_base::log_timeout_warnings(nvbench::printer_base &printer, + std::optional cuda_stdev_noise) +{ + if (!m_max_time_exceeded) + { + return; + } + + const auto timeout = m_walltime_timer.get_duration(); + + log_measurement_timeout_warnings(printer, + m_criterion_params, + timeout, + m_total_samples, + m_min_samples, + m_total_cuda_time, + cuda_stdev_noise); +} + void measure_cold_base::generate_summaries() { { @@ -211,6 +230,30 @@ void measure_cold_base::generate_summaries() summ.set_int64("value", m_total_samples); } + { + auto &summ = m_state.add_summary("nv/cold/walltime"); + summ.set_string("name", "Walltime"); + summ.set_string("hint", "duration"); + summ.set_string("description", "Walltime used for isolated measurements"); + summ.set_float64("value", m_walltime_timer.get_duration()); + summ.set_string("hide", "Hidden by default."); + } + + if (m_total_samples == 0) + { + // Throttling can discard every trial before a timeout stops collection. + // Keep timeout diagnostics, but skip sample-derived summaries. + if (auto printer_ptr = m_state.get_benchmark().get_printer()) + { + auto &printer = *printer_ptr; + this->log_timeout_warnings(printer); + printer.log(nvbench::log_level::warn, + fmt::format("Cold: no accepted samples recorded in {:0.2f}s walltime.", + m_walltime_timer.get_duration())); + } + return; + } + // cpu time statistics { auto &summ = m_state.add_summary("nv/cold/time/cpu/min"); @@ -257,15 +300,15 @@ void measure_cold_base::generate_summaries() summ.set_string("hide", "Hidden by default."); } - const auto cpu_stdev_noise = statistics::compute_relative_dispersion(cpu_stdev, cpu_mean); - if (cpu_stdev_noise) + const auto cpu_stdev_noise = + statistics::compute_standard_deviation_noise(m_total_samples, cpu_stdev, cpu_mean); { auto &summ = m_state.add_summary("nv/cold/time/cpu/stdev/relative"); summ.set_string("name", "Noise"); summ.set_string("hint", "percentage"); summ.set_string("description", "Relative standard deviation of isolated kernel execution CPU times"); - summ.set_float64("value", *cpu_stdev_noise); + summ.set_float64("value", statistics::stdev_noise_or_sentinel(cpu_stdev_noise)); } const auto [cpu_time_first_quartile, cpu_time_median, cpu_time_third_quartile] = @@ -363,15 +406,15 @@ void measure_cold_base::generate_summaries() summ.set_string("hide", "Hidden by default."); } - const auto cuda_stdev_noise = statistics::compute_relative_dispersion(cuda_stdev, cuda_mean); - if (cuda_stdev_noise) + const auto cuda_stdev_noise = + statistics::compute_standard_deviation_noise(m_total_samples, cuda_stdev, cuda_mean); { auto &summ = m_state.add_summary("nv/cold/time/gpu/stdev/relative"); summ.set_string("name", "Noise"); summ.set_string("hint", "percentage"); summ.set_string("description", "Relative standard deviation of isolated kernel execution GPU times"); - summ.set_float64("value", *cuda_stdev_noise); + summ.set_float64("value", statistics::stdev_noise_or_sentinel(cuda_stdev_noise)); } const auto [cuda_time_first_quartile, cuda_time_median, cuda_time_third_quartile] = @@ -460,15 +503,6 @@ void measure_cold_base::generate_summaries() } } // bandwidth - { - auto &summ = m_state.add_summary("nv/cold/walltime"); - summ.set_string("name", "Walltime"); - summ.set_string("hint", "duration"); - summ.set_string("description", "Walltime used for isolated measurements"); - summ.set_float64("value", m_walltime_timer.get_duration()); - summ.set_string("hide", "Hidden by default."); - } - if (m_sm_clock_rate_accumulator != 0.) { const auto clock_mean = m_sm_clock_rate_accumulator / d_samples; @@ -500,53 +534,7 @@ void measure_cold_base::generate_summaries() { auto &printer = *printer_ptr; - if (m_max_time_exceeded) - { - const auto timeout = m_walltime_timer.get_duration(); - - auto get_param = [this](std::optional ¶m, const std::string &name) { - if (m_criterion_params.has_value(name)) - { - param = m_criterion_params.get_float64(name); - } - }; - - std::optional max_noise; - get_param(max_noise, "max-noise"); - - std::optional min_time; - get_param(min_time, "min-time"); - - if (max_noise && cuda_stdev_noise && *cuda_stdev_noise > *max_noise) - { - printer.log(nvbench::log_level::warn, - fmt::format("Current measurement timed out ({:0.2f}s) " - "while over noise threshold ({:0.2f}% > " - "{:0.2f}%)", - timeout, - *cuda_stdev_noise * 100, - *max_noise * 100)); - } - if (m_total_samples < m_min_samples) - { - printer.log(nvbench::log_level::warn, - fmt::format("Current measurement timed out ({:0.2f}s) " - "before accumulating min_samples ({} < {})", - timeout, - m_total_samples, - m_min_samples)); - } - if (min_time && m_total_cuda_time < *min_time) - { - printer.log(nvbench::log_level::warn, - fmt::format("Current measurement timed out ({:0.2f}s) " - "before accumulating min_time ({:0.2f}s < " - "{:0.2f}s)", - timeout, - m_total_cuda_time, - *min_time)); - } - } + this->log_timeout_warnings(printer, cuda_stdev_noise); // Log to stdout: printer.log(nvbench::log_level::pass, diff --git a/nvbench/detail/measure_cold.cuh b/nvbench/detail/measure_cold.cuh index 5f5ce0b..df24519 100644 --- a/nvbench/detail/measure_cold.cuh +++ b/nvbench/detail/measure_cold.cuh @@ -46,12 +46,14 @@ #include #include +#include #include #include namespace nvbench { +struct printer_base; struct state; namespace detail @@ -96,6 +98,11 @@ protected: __forceinline__ void unblock_stream() { m_blocker.unblock(); } __forceinline__ void unblock_stream_noexcept() noexcept { m_blocker.unblock_noexcept(); } +private: + void log_timeout_warnings(nvbench::printer_base &printer, + std::optional cuda_stdev_noise = std::nullopt); + +protected: nvbench::state &m_state; nvbench::launch m_launch; diff --git a/nvbench/detail/measure_cpu_only.cxx b/nvbench/detail/measure_cpu_only.cxx index c291ba6..d8efdf3 100644 --- a/nvbench/detail/measure_cpu_only.cxx +++ b/nvbench/detail/measure_cpu_only.cxx @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -29,9 +30,7 @@ #include #include #include -#include #include -#include #include namespace nvbench::detail @@ -168,14 +167,14 @@ void measure_cpu_only_base::generate_summaries() summ.set_string("hide", "Hidden by default."); } - const auto cpu_stdev_noise = statistics::compute_relative_dispersion(cpu_stdev, cpu_mean); - if (cpu_stdev_noise) + const auto cpu_stdev_noise = + statistics::compute_standard_deviation_noise(m_total_samples, cpu_stdev, cpu_mean); { auto &summ = m_state.add_summary("nv/cpu_only/time/cpu/stdev/relative"); summ.set_string("name", "Noise"); summ.set_string("hint", "percentage"); summ.set_string("description", "Relative standard deviation of isolated CPU times"); - summ.set_float64("value", *cpu_stdev_noise); + summ.set_float64("value", statistics::stdev_noise_or_sentinel(cpu_stdev_noise)); } const auto [cpu_first_quartile, cpu_median, cpu_third_quartile] = @@ -268,48 +267,13 @@ void measure_cpu_only_base::generate_summaries() { const auto timeout = m_walltime_timer.get_duration(); - auto get_param = [this](std::optional ¶m, const std::string &name) { - if (m_criterion_params.has_value(name)) - { - param = m_criterion_params.get_float64(name); - } - }; - - std::optional max_noise; - get_param(max_noise, "max-noise"); - - std::optional min_time; - get_param(min_time, "min-time"); - - if (max_noise && cpu_stdev_noise && *cpu_stdev_noise > *max_noise) - { - printer.log(nvbench::log_level::warn, - fmt::format("Current measurement timed out ({:0.2f}s) " - "while over noise threshold ({:0.2f}% > " - "{:0.2f}%)", - timeout, - *cpu_stdev_noise * 100, - *max_noise * 100)); - } - if (m_total_samples < m_min_samples) - { - printer.log(nvbench::log_level::warn, - fmt::format("Current measurement timed out ({:0.2f}s) " - "before accumulating min_samples ({} < {})", - timeout, - m_total_samples, - m_min_samples)); - } - if (min_time && m_total_cpu_time < *min_time) - { - printer.log(nvbench::log_level::warn, - fmt::format("Current measurement timed out ({:0.2f}s) " - "before accumulating min_time ({:0.2f}s < " - "{:0.2f}s)", - timeout, - m_total_cpu_time, - *min_time)); - } + log_measurement_timeout_warnings(printer, + m_criterion_params, + timeout, + m_total_samples, + m_min_samples, + m_total_cpu_time, + cpu_stdev_noise); } // Log to stdout: diff --git a/nvbench/detail/measure_timeout_warnings.cuh b/nvbench/detail/measure_timeout_warnings.cuh new file mode 100644 index 0000000..7eb1451 --- /dev/null +++ b/nvbench/detail/measure_timeout_warnings.cuh @@ -0,0 +1,117 @@ +/* + * Copyright 2026 NVIDIA Corporation + * + * Licensed under the Apache License, Version 2.0 with the LLVM exception + * (the "License"); you may not use this file except in compliance with + * the License. + * + * You may obtain a copy of the License at + * + * http://llvm.org/foundation/relicensing/LICENSE.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once + +#include + +#if defined(NVBENCH_IMPLICIT_SYSTEM_HEADER_GCC) +#pragma GCC system_header +#elif defined(NVBENCH_IMPLICIT_SYSTEM_HEADER_CLANG) +#pragma clang system_header +#elif defined(NVBENCH_IMPLICIT_SYSTEM_HEADER_MSVC) +#pragma system_header +#endif + +#include +#include +#include +#include + +#include + +#include +#include +#include + +namespace nvbench::detail +{ + +inline std::optional +get_float64_criterion_param(const nvbench::criterion_params ¶ms, const std::string &name) +{ + if (!params.has_value(name)) + { + return std::nullopt; + } + return params.get_float64(name); +} + +inline void log_measurement_timeout_warnings(nvbench::printer_base &printer, + const nvbench::criterion_params &criterion_params, + nvbench::float64_t timeout, + nvbench::int64_t total_samples, + nvbench::int64_t min_samples, + nvbench::float64_t accumulated_time, + std::optional stdev_noise) +{ + const auto max_noise = get_float64_criterion_param(criterion_params, "max-noise"); + const auto min_time = get_float64_criterion_param(criterion_params, "min-time"); + + const auto enough_samples_for_noise = + statistics::has_enough_samples_for_noise_estimate(total_samples); + const auto stdev_noise_unavailable = !stdev_noise || std::isnan(*stdev_noise) || + *stdev_noise < nvbench::float64_t{}; + if (max_noise && !enough_samples_for_noise) + { + printer.log(nvbench::log_level::warn, + fmt::format("Current measurement timed out ({:0.2f}s) " + "before accumulating enough samples to estimate noise ({} < {})", + timeout, + total_samples, + statistics::min_samples_for_noise_estimate)); + } + else if (max_noise && stdev_noise_unavailable) + { + printer.log(nvbench::log_level::warn, + fmt::format("Current measurement timed out ({:0.2f}s) " + "while unable to estimate noise for max-noise", + timeout)); + } + else if (max_noise && *stdev_noise > *max_noise) + { + printer.log(nvbench::log_level::warn, + fmt::format("Current measurement timed out ({:0.2f}s) " + "while over noise threshold ({:0.2f}% > " + "{:0.2f}%)", + timeout, + *stdev_noise * 100, + *max_noise * 100)); + } + if (total_samples < min_samples) + { + printer.log(nvbench::log_level::warn, + fmt::format("Current measurement timed out ({:0.2f}s) " + "before accumulating min_samples ({} < {})", + timeout, + total_samples, + min_samples)); + } + if (min_time && accumulated_time < *min_time) + { + printer.log(nvbench::log_level::warn, + fmt::format("Current measurement timed out ({:0.2f}s) " + "before accumulating min_time ({:0.2f}s < " + "{:0.2f}s)", + timeout, + accumulated_time, + *min_time)); + } +} + +} // namespace nvbench::detail diff --git a/nvbench/detail/statistics.cuh b/nvbench/detail/statistics.cuh index ab12fa0..33bc376 100644 --- a/nvbench/detail/statistics.cuh +++ b/nvbench/detail/statistics.cuh @@ -54,6 +54,21 @@ namespace nvbench::detail::statistics inline constexpr nvbench::int64_t min_samples_for_noise_estimate = 5; +// Heuristic crossover near L1-sized float64 data; tune if sorting remains faster. +inline constexpr std::size_t quartile_selection_threshold = 4096; + +inline constexpr bool has_enough_samples_for_noise_estimate(nvbench::int64_t num_samples) +{ + return num_samples >= min_samples_for_noise_estimate; +} + +template +constexpr ValueType standard_deviation_unavailable_sentinel() +{ + static_assert(std::is_floating_point_v); + return std::numeric_limits::infinity(); +} + /** * Computes and returns the unbiased sample standard deviation. * @@ -66,9 +81,9 @@ ValueType standard_deviation(Iter first, Iter last, ValueType mean) const auto num = std::distance(first, last); - if (num < min_samples_for_noise_estimate) // don't bother with low sample sizes. + if (!has_enough_samples_for_noise_estimate(num)) // don't bother with low sample sizes. { - return std::numeric_limits::infinity(); + return standard_deviation_unavailable_sentinel(); } const auto variance = nvbench::detail::transform_reduce(first, @@ -205,10 +220,13 @@ public: } }; -// Compute percentile rank using nearest rank method +// Use a rounded zero-based percentile rank inline std::size_t percentile_rank(int percentile, std::size_t size) { + // Precondition: sample_size > 0. Public percentile helpers handle empty + // inputs before calling this internal rank helper. assert(size > 0 && "percentile_rank requires non-empty sample set"); + const auto p = std::clamp(percentile, 0, 100); const auto q = static_cast(p) / 100.0; @@ -216,12 +234,23 @@ inline std::size_t percentile_rank(int percentile, std::size_t size) return static_cast(std::round(q * max_rank)); } +template +bool contains_nan(const std::vector &samples) +{ + return std::any_of(samples.cbegin(), samples.cend(), [](ValueType value) { + return std::isnan(value); + }); +} + template std::array compute_percentiles_by_sorting(std::vector &&samples, const std::array &percentiles) { + static_assert(std::is_floating_point_v, + "compute_percentiles_by_sorting requires a floating-point value type."); + std::array result{}; - if (samples.empty()) + if (samples.empty() || contains_nan(samples)) { result.fill(std::numeric_limits::quiet_NaN()); return result; @@ -242,7 +271,8 @@ std::array compute_percentiles_by_sorting(std::vector & * Computes exact percentile values using rank round(p / 100 * (S - 1)). * * The input range is copied before sorting, so const iterators are supported. - * If the input has fewer than 1 sample, all percentiles are returned as quiet NaNs. + * If the input has fewer than 1 sample or contains a NaN, all percentiles are returned as quiet + * NaNs. */ template quartiles_t compute_quartiles_by_sorting(std::vector &&samples) { + static_assert(std::is_floating_point_v, + "compute_quartiles_by_sorting requires a floating-point value type."); + constexpr std::array qs{25, 50, 75}; - const auto r = ::nvbench::detail::statistics::compute_percentiles_by_sorting( - std::forward>(samples), - qs); + const auto r = ::nvbench::detail::statistics::compute_percentiles_by_sorting(std::move(samples), + qs); return {r[0], r[1], r[2]}; } template quartiles_t compute_quartiles_by_selection(std::vector &&samples) { - if (samples.empty()) + static_assert(std::is_floating_point_v, + "compute_quartiles_by_selection requires a floating-point value type."); + + if (samples.empty() || contains_nan(samples)) { constexpr auto nan = std::numeric_limits::quiet_NaN(); return {nan, nan, nan}; @@ -319,9 +354,8 @@ quartiles_t compute_quartiles(Iter first, Iter last) static_assert(std::is_floating_point_v); std::vector samples(first, last); - constexpr std::size_t selection_threshold = 4096; - if (samples.size() >= selection_threshold) + if (samples.size() >= quartile_selection_threshold) { return ::nvbench::detail::statistics::compute_quartiles_by_selection(std::move(samples)); } @@ -357,13 +391,36 @@ compute_relative_interquartile_range(nvbench::float64_t first_quartile, return ::nvbench::detail::statistics::compute_relative_dispersion(interquartile_range, median); } +// Returns nullopt until there are enough samples for a meaningful standard deviation estimate. +inline std::optional +compute_standard_deviation_noise(nvbench::int64_t num_samples, + nvbench::float64_t standard_deviation, + nvbench::float64_t center) +{ + if (!has_enough_samples_for_noise_estimate(num_samples)) + { + return std::nullopt; + } + if (!std::isfinite(standard_deviation)) + { + return std::nullopt; + } + + return ::nvbench::detail::statistics::compute_relative_dispersion(standard_deviation, center); +} + +inline nvbench::float64_t stdev_noise_or_sentinel(std::optional noise) +{ + return noise.value_or(standard_deviation_unavailable_sentinel()); +} + // Returns nullopt until there are enough samples for a meaningful robust noise estimate. inline std::optional compute_robust_noise(nvbench::int64_t num_samples, nvbench::float64_t first_quartile, nvbench::float64_t median, nvbench::float64_t third_quartile) { - if (num_samples < min_samples_for_noise_estimate) + if (!has_enough_samples_for_noise_estimate(num_samples)) { return std::nullopt; } diff --git a/testing/CMakeLists.txt b/testing/CMakeLists.txt index 7dc6cef..b390edf 100644 --- a/testing/CMakeLists.txt +++ b/testing/CMakeLists.txt @@ -16,6 +16,7 @@ set(test_srcs exception_safety.cu float64_axis.cu int64_axis.cu + measure_timeout_warnings.cu named_values.cu option_parser.cu range.cu diff --git a/testing/measure_timeout_warnings.cu b/testing/measure_timeout_warnings.cu new file mode 100644 index 0000000..58372f3 --- /dev/null +++ b/testing/measure_timeout_warnings.cu @@ -0,0 +1,118 @@ +/* + * Copyright 2026 NVIDIA Corporation + * + * Licensed under the Apache License, Version 2.0 with the LLVM exception + * (the "License"); you may not use this file except in compliance with + * the License. + * + * You may obtain a copy of the License at + * + * http://llvm.org/foundation/relicensing/LICENSE.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "test_asserts.cuh" + +struct recording_printer : nvbench::printer_base +{ + explicit recording_printer(std::ostream &stream) + : nvbench::printer_base{stream} + {} + + std::vector> logs; + +protected: + void do_log(nvbench::log_level level, const std::string &message) override + { + logs.emplace_back(level, message); + } +}; + +void check_noise_warning( + std::optional stdev_noise, + const std::string &expected_message, + nvbench::int64_t total_samples = nvbench::detail::statistics::min_samples_for_noise_estimate) +{ + std::ostringstream stream; + recording_printer printer{stream}; + nvbench::criterion_params params; + params.set_float64("max-noise", 0.01); + + nvbench::detail::log_measurement_timeout_warnings(printer, + params, + 1.0, + total_samples, + 1, + 1.0, + stdev_noise); + + ASSERT(printer.logs.size() == 1); + ASSERT(printer.logs[0].first == nvbench::log_level::warn); + ASSERT(printer.logs[0].second.find(expected_message) != std::string::npos); +} + +void test_non_finite_or_invalid_stdev_noise_timeout_warning() +{ + check_noise_warning(std::nullopt, + "before accumulating enough samples to estimate noise", + nvbench::detail::statistics::min_samples_for_noise_estimate - 1); + check_noise_warning(std::nullopt, "unable to estimate noise"); + check_noise_warning(std::numeric_limits::quiet_NaN(), + "unable to estimate noise"); + check_noise_warning(-1.0, "unable to estimate noise"); + check_noise_warning(std::numeric_limits::infinity(), "over noise threshold"); +} + +void test_min_samples_timeout_warning() +{ + std::ostringstream stream; + recording_printer printer{stream}; + nvbench::criterion_params params; + + nvbench::detail::log_measurement_timeout_warnings(printer, params, 1.0, 4, 5, 1.0, std::nullopt); + + ASSERT(printer.logs.size() == 1); + ASSERT(printer.logs[0].first == nvbench::log_level::warn); + ASSERT(printer.logs[0].second.find("before accumulating min_samples") != std::string::npos); +} + +void test_min_time_timeout_warning() +{ + std::ostringstream stream; + recording_printer printer{stream}; + nvbench::criterion_params params; + params.set_float64("min-time", 2.0); + + nvbench::detail::log_measurement_timeout_warnings(printer, params, 1.0, 5, 1, 1.5, std::nullopt); + + ASSERT(printer.logs.size() == 1); + ASSERT(printer.logs[0].first == nvbench::log_level::warn); + ASSERT(printer.logs[0].second.find("before accumulating min_time") != std::string::npos); +} + +int main() +{ + test_non_finite_or_invalid_stdev_noise_timeout_warning(); + test_min_samples_timeout_warning(); + test_min_time_timeout_warning(); + return 0; +} diff --git a/testing/statistics.cu b/testing/statistics.cu index 7aa995a..9c8d2b9 100644 --- a/testing/statistics.cu +++ b/testing/statistics.cu @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -63,6 +64,27 @@ void assert_quartiles_nan(statistics::quartiles_t actual) ASSERT(std::isnan(actual.third_quartile)); } +statistics::quartiles_t expected_rank_quartiles(std::size_t num_samples) +{ + const auto expected_value = [num_samples](int percentile) { + const auto q = static_cast(percentile) / 100.0; + return std::round(q * static_cast(num_samples - 1)); + }; + return {expected_value(25), expected_value(50), expected_value(75)}; +} + +statistics::quartiles_t +expected_duplicate_heavy_quartiles(std::size_t num_samples) +{ + const auto value_at_percentile = [num_samples](int percentile) { + const auto q = static_cast(percentile) / 100.0; + const auto rank = + static_cast(std::round(q * static_cast(num_samples - 1))); + return static_cast((4 * rank) / num_samples); + }; + return {value_at_percentile(25), value_at_percentile(50), value_at_percentile(75)}; +} + void test_mean() { { @@ -259,6 +281,16 @@ void test_percentiles() ASSERT(std::isnan(actual[1])); ASSERT(std::isnan(actual[2])); } + + { + constexpr auto nan = std::numeric_limits::quiet_NaN(); + const std::vector data{10.0, nan, 30.0, 20.0}; + const auto actual = + statistics::compute_percentiles(data.cbegin(), data.cend(), std::array{25, 50, 75}); + ASSERT(std::isnan(actual[0])); + ASSERT(std::isnan(actual[1])); + ASSERT(std::isnan(actual[2])); + } } void test_quartiles_methods_agree() @@ -282,14 +314,32 @@ void test_quartiles_methods_agree() assert_quartiles_equal(selection, sorting); } + { + constexpr auto nan = std::numeric_limits::quiet_NaN(); + const std::vector data{40.0, 10.0, nan, 20.0}; + assert_quartiles_nan( + statistics::compute_quartiles_by_sorting(std::vector(data))); + assert_quartiles_nan( + statistics::compute_quartiles_by_selection(std::vector(data))); + assert_quartiles_nan(statistics::compute_quartiles(data.cbegin(), data.cend())); + } + // test around threshold when public API switches between implementations - for (const auto n : std::array{4095, 4096, 4097}) + constexpr auto threshold = statistics::quartile_selection_threshold; + if constexpr (threshold < 2) + { + return; + } + + for (const auto n : std::array{threshold - 1, threshold, threshold + 1}) { std::vector data(n); for (std::size_t i = 0; i < data.size(); ++i) { - data[i] = static_cast((i * 37) % data.size()); + data[i] = static_cast(i); } + std::mt19937 rng{37u}; + std::shuffle(data.begin(), data.end(), rng); const auto public_api = statistics::compute_quartiles(data.cbegin(), data.cend()); const auto sorting = @@ -298,6 +348,41 @@ void test_quartiles_methods_agree() statistics::compute_quartiles_by_selection(std::vector(data)); assert_quartiles_equal(selection, sorting); assert_quartiles_equal(public_api, sorting); + assert_quartiles_equal(public_api, expected_rank_quartiles(n)); + } +} + +void test_quartiles_methods_agree_with_duplicate_heavy_inputs() +{ + // Test around threshold when public API switches between implementations. + constexpr auto threshold = statistics::quartile_selection_threshold; + if constexpr (threshold < 2) + { + return; + } + + for (const auto n : std::array{threshold - 1, threshold, threshold + 1}) + { + for (const auto seed : std::array{17u, 12345u, 987654321u}) + { + std::vector data(n); + for (std::size_t i = 0; i < data.size(); ++i) + { + data[i] = static_cast((4 * i) / data.size()); + } + + std::mt19937 rng{seed}; + std::shuffle(data.begin(), data.end(), rng); + + const auto public_api = statistics::compute_quartiles(data.cbegin(), data.cend()); + const auto sorting = + statistics::compute_quartiles_by_sorting(std::vector(data)); + const auto selection = + statistics::compute_quartiles_by_selection(std::vector(data)); + assert_quartiles_equal(selection, sorting); + assert_quartiles_equal(public_api, sorting); + assert_quartiles_equal(public_api, expected_duplicate_heavy_quartiles(n)); + } } } @@ -385,6 +470,119 @@ void test_compute_robust_noise() ASSERT(actual); ASSERT(is_close(*actual, 1.0)); } + + { + const auto actual = + statistics::compute_robust_noise(statistics::min_samples_for_noise_estimate, 0.0, 0.0, 1.0); + ASSERT(!actual); + } + + { + const auto actual = + statistics::compute_robust_noise(statistics::min_samples_for_noise_estimate, -2.0, -1.0, 0.0); + ASSERT(!actual); + } + + { + const auto actual = + statistics::compute_robust_noise(statistics::min_samples_for_noise_estimate, + std::numeric_limits::quiet_NaN(), + 4.0, + 6.0); + ASSERT(!actual); + } + + { + const auto actual = + statistics::compute_robust_noise(statistics::min_samples_for_noise_estimate, + 2.0, + 4.0, + std::numeric_limits::infinity()); + ASSERT(!actual); + } +} + +void test_compute_standard_deviation_noise() +{ + ASSERT(!statistics::has_enough_samples_for_noise_estimate( + statistics::min_samples_for_noise_estimate - 1)); + ASSERT( + statistics::has_enough_samples_for_noise_estimate(statistics::min_samples_for_noise_estimate)); + + { + const auto actual = + statistics::compute_standard_deviation_noise(statistics::min_samples_for_noise_estimate - 1, + 2.0, + 1.0); + ASSERT(!actual); + } + + { + const auto actual = statistics::compute_standard_deviation_noise( + statistics::min_samples_for_noise_estimate, + std::numeric_limits::quiet_NaN(), + 1.0); + ASSERT(!actual); + } + + { + const auto actual = statistics::compute_standard_deviation_noise( + statistics::min_samples_for_noise_estimate, + std::numeric_limits::infinity(), + 1.0); + ASSERT(!actual); + } + + { + const auto actual = + statistics::compute_standard_deviation_noise(statistics::min_samples_for_noise_estimate, + 1.0, + 0.0); + ASSERT(!actual); + } + + { + const auto actual = + statistics::compute_standard_deviation_noise(statistics::min_samples_for_noise_estimate, + 1.0, + -1.0); + ASSERT(!actual); + } + + { + const auto actual = + statistics::compute_standard_deviation_noise(statistics::min_samples_for_noise_estimate, + -1.0, + 1.0); + ASSERT(!actual); + } + + { + const auto actual = + statistics::compute_standard_deviation_noise(statistics::min_samples_for_noise_estimate, + 2.0, + 4.0); + ASSERT(actual); + ASSERT(is_close(*actual, 0.5)); + } +} + +void test_stdev_noise_or_sentinel() +{ + { + const auto actual = statistics::standard_deviation_unavailable_sentinel(); + ASSERT(std::isinf(actual)); + } + + { + const auto actual = statistics::stdev_noise_or_sentinel(nvbench::float64_t{0.25}); + ASSERT(is_close(actual, 0.25)); + } + + { + const auto actual = statistics::stdev_noise_or_sentinel(std::nullopt); + ASSERT(actual == statistics::standard_deviation_unavailable_sentinel()); + } } void test_relative_interquartile_range() @@ -496,9 +694,12 @@ int main() test_percentiles(); test_quartiles(); test_quartiles_methods_agree(); + test_quartiles_methods_agree_with_duplicate_heavy_inputs(); test_compute_relative_dispersion_nominal_input(); test_compute_relative_dispersion_invalid_inputs(); test_relative_interquartile_range(); + test_compute_standard_deviation_noise(); + test_stdev_noise_or_sentinel(); test_compute_robust_noise(); test_lin_regression(); test_r2();