Fix standard deviation

Author: leandrolcampos

libc/benchmarks/gpu/LibcGpuBenchmark.cpp

@@ -166,15 +166,13 @@ benchmark(const BenchmarkOptions &options,
           const cpp::function<uint64_t(uint32_t)> &wrapper_func) {
   BenchmarkResult result;
   RuntimeEstimationProgression rep;
-  uint32_t total_iterations = 0;
   uint32_t iterations = options.initial_iterations;
+
   if (iterations < 1u)
     iterations = 1;
 
   uint32_t samples = 0;
   uint64_t total_time = 0;
-  uint64_t best_guess = 0;
-  uint64_t cycles_squared = 0;
   uint64_t min = UINT64_MAX;
   uint64_t max = 0;
 
@@ -186,46 +184,55 @@ benchmark(const BenchmarkOptions &options,
   uint32_t call_index = 0;
 
   for (int64_t time_budget = options.max_duration; time_budget >= 0;) {
-    uint64_t sample_cycles = 0;
-    const clock_t start = static_cast<double>(clock());
-    for (uint32_t i = 0; i < iterations; i++) {
+    RefinableRuntimeEstimator sample_estimator;
+
+    const clock_t start = clock();
+    while (sample_estimator.get_iterations() < iterations) {
       auto wrapper_intermediate = wrapper_func(call_index++);
-      uint64_t current_result = wrapper_intermediate - overhead;
+      uint64_t current_result =
+          wrapper_intermediate < overhead ? 0 : wrapper_intermediate - overhead;
       max = cpp::max(max, current_result);
       min = cpp::min(min, current_result);
-      sample_cycles += current_result;
+      sample_estimator.update(current_result);
     }
     const clock_t end = clock();
+
     const clock_t duration_ns =
         ((end - start) * 1000 * 1000 * 1000) / CLOCKS_PER_SEC;
     total_time += duration_ns;
     time_budget -= duration_ns;
     samples++;
-    cycles_squared += sample_cycles * sample_cycles;
 
-    total_iterations += iterations;
-    const double change_ratio =
-        rep.compute_improvement({iterations, sample_cycles});
-    best_guess = rep.current_estimation;
+    const double change_ratio = rep.compute_improvement(sample_estimator);
 
     if (samples >= options.max_samples || iterations >= options.max_iterations)
       break;
+
+    const auto total_iterations = rep.get_estimator().get_iterations();
+
     if (total_time >= options.min_duration && samples >= options.min_samples &&
         total_iterations >= options.min_iterations &&
         change_ratio < options.epsilon)
       break;
 
-    iterations *= options.scaling_factor;
+    iterations = static_cast<uint32_t>(iterations * options.scaling_factor);
   }
-  result.cycles = best_guess;
-  result.standard_deviation = fputil::sqrt<double>(
-      static_cast<double>(cycles_squared) / total_iterations -
-      static_cast<double>(best_guess * best_guess));
+
+  const auto &estimator = rep.get_estimator();
+  result.cycles = static_cast<uint64_t>(estimator.get_mean());
+  result.standard_deviation = estimator.get_stddev();
+
   result.min = min;
   result.max = max;
   result.samples = samples;
-  result.total_iterations = total_iterations;
-  result.total_time = total_time / total_iterations;
+
+  result.total_iterations = estimator.get_iterations();
+  if (result.total_iterations > 0) {
+    result.total_time = total_time / result.total_iterations;
+  } else {
+    result.total_time = 0;
+  }
+
   return result;
 }
 

libc/benchmarks/gpu/LibcGpuBenchmark.h

@@ -11,6 +11,7 @@
 #include "src/__support/CPP/string_view.h"
 #include "src/__support/CPP/type_traits.h"
 #include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/sqrt.h"
 #include "src/__support/macros/config.h"
 #include "src/time/clock.h"
 
@@ -30,40 +31,82 @@ struct BenchmarkOptions {
   double scaling_factor = 1.4;
 };
 
-struct Measurement {
+class RefinableRuntimeEstimator {
   uint32_t iterations = 0;
-  uint64_t elapsed_cycles = 0;
-};
-
-class RefinableRuntimeEstimation {
-  uint64_t total_cycles = 0;
-  uint32_t total_iterations = 0;
+  uint64_t sum_of_cycles = 0;
+  uint64_t sum_of_squared_cycles = 0;
 
 public:
-  uint64_t update(const Measurement &M) {
-    total_cycles += M.elapsed_cycles;
-    total_iterations += M.iterations;
-    return total_cycles / total_iterations;
+  void update(uint64_t cycles) noexcept {
+    iterations += 1;
+    sum_of_cycles += cycles;
+    sum_of_squared_cycles += cycles * cycles;
+  }
+
+  double get_mean() const noexcept {
+    if (iterations == 0)
+      return 0.0;
+
+    return static_cast<double>(sum_of_cycles) / iterations;
+  }
+
+  void update(const RefinableRuntimeEstimator &other) noexcept {
+    iterations += other.iterations;
+    sum_of_cycles += other.sum_of_cycles;
+    sum_of_squared_cycles += other.sum_of_squared_cycles;
+  }
+
+  double get_variance() const noexcept {
+    if (iterations == 0)
+      return 0.0;
+
+    const double num = static_cast<double>(iterations);
+    const double sum_x = static_cast<double>(sum_of_cycles);
+    const double sum_x2 = static_cast<double>(sum_of_squared_cycles);
+
+    const double mean_of_squares = sum_x2 / num;
+    const double mean = sum_x / num;
+    const double mean_squared = mean * mean;
+    const double variance = mean_of_squares - mean_squared;
+
+    return variance < 0.0 ? 0.0 : variance;
+  }
+
+  double get_stddev() const noexcept {
+    return fputil::sqrt<double>(get_variance());
   }
+
+  uint32_t get_iterations() const noexcept { return iterations; }
 };
 
 // Tracks the progression of the runtime estimation
 class RuntimeEstimationProgression {
-  RefinableRuntimeEstimation rre;
+  RefinableRuntimeEstimator estimator;
+  double current_mean = 0.0;
 
 public:
-  uint64_t current_estimation = 0;
+  const RefinableRuntimeEstimator &get_estimator() const noexcept {
+    return estimator;
+  }
+
+  double
+  compute_improvement(const RefinableRuntimeEstimator &sample_estimator) {
+    if (sample_estimator.get_iterations() == 0)
+      return 1.0;
 
-  double compute_improvement(const Measurement &M) {
-    const uint64_t new_estimation = rre.update(M);
-    double ratio =
-        (static_cast<double>(current_estimation) / new_estimation) - 1.0;
+    estimator.update(sample_estimator);
+
+    const double new_mean = estimator.get_mean();
+    if (current_mean == 0.0 || new_mean == 0.0) {
+      current_mean = new_mean;
+      return 1.0;
+    }
 
-    // Get absolute value
+    double ratio = (current_mean / new_mean) - 1.0;
     if (ratio < 0)
-      ratio *= -1;
+      ratio = -ratio;
 
-    current_estimation = new_estimation;
+    current_mean = new_mean;
     return ratio;
   }
 };