test/bench: measure standard deviation for benchmarks (#6390)

* Add Welfords square of sums algorithm for variance and std deviation calculations.
* Add standard deviation to benchmarks
* Simplify Welford
* Add geometric mean and standard deviation to Welford Sums
This commit is contained in:
Rot127 2026-05-24 21:08:06 +00:00 committed by GitHub
parent 6be10c2428
commit 702250eb4f
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
9 changed files with 263 additions and 4 deletions

30
librz/include/rz_math.h Normal file
View file

@ -0,0 +1,30 @@
// SPDX-FileCopyrightText: 2026 Rot127 <rot127@posteo.com>
// SPDX-License-Identifier: LGPL-3.0-only
#ifndef RZ_MATH_H
#define RZ_MATH_H
#include <rz_types.h>
/**
* \brief A Welford Variance implementation taken from
* https://www.johndcook.com/blog/standard_deviation/
*/
typedef struct {
ut64 n; ///< Number of variables
double old_mean; ///< The mean before a variable is added.
double new_mean; ///< The mean after a variable was added.
double oldS; ///< Sum of squares before a variable is added.
double newS; ///< Sum of squares after a variable was added.
} RzMathWelfordSums;
RZ_API void rz_math_welford_init(RZ_BORROW RzMathWelfordSums *wf);
RZ_API void rz_math_welford_clear(RZ_BORROW RzMathWelfordSums *wf);
RZ_API void rz_math_welford_push(RZ_BORROW RzMathWelfordSums *wf, double var);
RZ_API ut64 rz_math_welford_n(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_mean(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_variance(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_std_deviation(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_sum_of_squares(const RzMathWelfordSums *wf);
#endif // RZ_MATH_H

View file

@ -48,6 +48,7 @@
#include <rz_util/rz_lang_byte_array.h>
#include <rz_util/rz_log.h>
#include <rz_util/rz_luhn.h>
#include <rz_util/rz_math.h>
#include <rz_util/rz_mem.h>
#include <rz_util/rz_name.h>
#include <rz_util/rz_num.h>

View file

@ -0,0 +1,44 @@
// SPDX-FileCopyrightText: 2026 Rot127 <rot127@posteo.com>
// SPDX-License-Identifier: LGPL-3.0-only
#ifndef RZ_MATH_H
#define RZ_MATH_H
#include <rz_types.h>
#ifdef __cplusplus
extern "C" {
#endif
/**
* \brief A Welford Sums of Squares implementation take from
* https://www.johndcook.com/blog/standard_deviation/
* doi: http://dx.doi.org/10.1080/00401706.1962.10490022
*/
typedef struct {
ut64 n; ///< Number of variables
double amean; ///< The arithmetic mean of the variables.
double asums; ///< Sum of squares.
double gmean; ///< Geometric mean
double ln_v_sums; ///< Geometric sum of ln(x_i)
double gsums; ///< Geometric sums of squares
} RzMathWelfordSums;
RZ_API void rz_math_welford_init(RZ_BORROW RzMathWelfordSums *wf);
RZ_API void rz_math_welford_clear(RZ_BORROW RzMathWelfordSums *wf);
RZ_API void rz_math_welford_push(RZ_BORROW RzMathWelfordSums *wf, double var);
RZ_API ut64 rz_math_welford_n(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_amean(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_avar(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_astddev(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_asums(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_gmean(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_gvar(const RzMathWelfordSums *wf);
RZ_API double rz_math_welford_gstddev(const RzMathWelfordSums *wf);
#ifdef __cplusplus
}
#endif
#endif // RZ_MATH_H

70
librz/util/math.c Normal file
View file

@ -0,0 +1,70 @@
// SPDX-FileCopyrightText: 2026 Rot127 <rot127@posteo.com>
// SPDX-License-Identifier: LGPL-3.0-only
#include <rz_util.h>
RZ_API void rz_math_welford_init(RZ_BORROW RzMathWelfordSums *wf) {
rz_return_if_fail(wf);
memset(wf, 0, sizeof(RzMathWelfordSums));
}
RZ_API void rz_math_welford_clear(RZ_BORROW RzMathWelfordSums *wf) {
rz_return_if_fail(wf);
memset(wf, 0, sizeof(RzMathWelfordSums));
}
RZ_API void rz_math_welford_push(RZ_BORROW RzMathWelfordSums *wf, double var) {
rz_return_if_fail(wf);
wf->n++;
// See Knuth TAOCP vol 2, 3rd edition, page 232
if (wf->n == 1) {
wf->amean = var;
wf->asums = 0.0;
wf->gmean = var;
wf->ln_v_sums = log(var);
wf->gsums = 0.0;
} else {
// Arithmetic
double old_amean = wf->amean;
wf->amean = old_amean + (var - old_amean) / wf->n;
wf->asums += (var - old_amean) * (var - wf->amean);
// Geometric
double old_gmean = wf->gmean;
wf->ln_v_sums += log(var);
wf->gmean = exp(wf->ln_v_sums / wf->n);
wf->gsums += log(var / old_gmean) * log(var / wf->gmean);
}
}
RZ_API ut64 rz_math_welford_n(const RzMathWelfordSums *wf) {
rz_return_val_if_fail(wf, 0);
return wf->n;
}
RZ_API double rz_math_welford_amean(const RzMathWelfordSums *wf) {
rz_return_val_if_fail(wf, 0.0);
return (wf->n > 0) ? wf->amean : 0.0;
}
RZ_API double rz_math_welford_gmean(const RzMathWelfordSums *wf) {
rz_return_val_if_fail(wf, 0.0);
return (wf->n > 0) ? wf->gmean : 0.0;
}
RZ_API double rz_math_welford_avar(const RzMathWelfordSums *wf) {
rz_return_val_if_fail(wf, 0.0);
return ((wf->n > 1) ? wf->asums / wf->n : 0.0);
}
RZ_API double rz_math_welford_astddev(const RzMathWelfordSums *wf) {
rz_return_val_if_fail(wf, 0.0);
return sqrt(rz_math_welford_avar(wf));
}
RZ_API double rz_math_welford_gstddev(const RzMathWelfordSums *wf) {
rz_return_val_if_fail(wf, 0.0);
return (wf->n > 1) ? exp(sqrt(wf->gsums / (wf->n - 1))) : 0.0;
}

View file

@ -35,6 +35,7 @@ rz_util_common_sources = [
'list.c',
'log.c',
'luhn.c',
'math.c',
'mem.c',
'name.c',
'path.c',

View file

@ -2,7 +2,6 @@
// SPDX-License-Identifier: LGPL-3.0-only
#include "bench_utils.h"
#include <stdio.h>
/**
* \brief Initialize benchmark context
@ -38,8 +37,7 @@ RZ_API void rz_bench_report(RZ_NONNULL RzBenchCtx *ctx, RZ_NONNULL RzTable *t) {
rz_return_if_fail(ctx && t);
double total_ms = ctx->total_time / 1000.0;
double avg_us = ctx->iterations != 0 ? (double)ctx->total_time / ctx->iterations : 0;
double ops_per_sec = ctx->total_time != 0 ? (ctx->iterations * 1000000.0) / ctx->total_time : 0;
rz_table_add_rowf(t, "sdfff", ctx->name, (int)ctx->iterations, total_ms, avg_us, ops_per_sec);
rz_table_add_rowf(t, "sdffff", ctx->name, (int)ctx->iterations, total_ms, ctx->arith_mean_us, ops_per_sec, ctx->arith_std_dev);
}

View file

@ -15,6 +15,10 @@ typedef struct rz_bench_ctx_t {
ut64 iterations; ///< number of iterations
ut64 start_time; ///< start time of the benchmark in microseconds
ut64 total_time; ///< total elapsed time of the benchmark in microseconds
double arith_mean_us; ///< The average time needed per iteration in microseconds (arithmetic mean).
double arith_std_dev; ///< Arithmetic standard deviation of benchmark samples in microseconds.
double geo_mean_us; ///< The average time needed per iteration in microseconds (geometric mean).
double geo_std_dev; ///< Geometric standard deviation of benchmark samples.
} RzBenchCtx;
RZ_API void rz_bench_init(RZ_NONNULL RzBenchCtx *ctx, RZ_NONNULL const char *name, ut64 iterations);
@ -34,25 +38,39 @@ RZ_API void rz_bench_report(RZ_NONNULL RzBenchCtx *ctx, RZ_NONNULL RzTable *t);
*/
#define RZ_BENCH_RUN(name, table, iterations, code) \
do { \
RzMathWelfordSums wf = { 0 }; \
RzBenchCtx ctx; \
rz_bench_init(&ctx, name, iterations); \
rz_bench_start(&ctx); \
for (ut64 i = 0; i < iterations; i++) { \
ut64 spl = rz_time_now_mono(); \
code; \
rz_math_welford_push(&wf, (double)(rz_time_now_mono() - spl)); \
} \
ctx.arith_std_dev = rz_math_welford_astddev(&wf); \
ctx.arith_mean_us = rz_math_welford_amean(&wf); \
ctx.geo_std_dev = rz_math_welford_gstddev(&wf); \
ctx.geo_mean_us = rz_math_welford_gmean(&wf); \
rz_bench_end(&ctx); \
rz_bench_report(&ctx, table); \
} while (0)
#define RZ_BENCH_RUN_I(name, i, table, iterations, code) \
do { \
RzMathWelfordSums wf = { 0 }; \
RzBenchCtx ctx; \
rz_bench_init(&ctx, name, iterations); \
rz_bench_start(&ctx); \
for (ut64(i) = 0; (i) < iterations; (i)++) { \
ut64 spl = rz_time_now_mono(); \
code; \
rz_math_welford_push(&wf, (double)(rz_time_now_mono() - spl)); \
} \
rz_bench_end(&ctx); \
ctx.arith_std_dev = rz_math_welford_astddev(&wf); \
ctx.arith_mean_us = rz_math_welford_amean(&wf); \
ctx.geo_std_dev = rz_math_welford_gstddev(&wf); \
ctx.geo_mean_us = rz_math_welford_gmean(&wf); \
rz_bench_report(&ctx, table); \
} while (0)
@ -61,7 +79,7 @@ RZ_API void rz_bench_report(RZ_NONNULL RzBenchCtx *ctx, RZ_NONNULL RzTable *t);
* \param T table to initialize.
*/
#define RZ_BENCH_TABLE_INIT(T) \
rz_table_set_columnsf(T, "snnnn", "Benchmark", "Iterations", "Total time [ms]", "Avg iter time [us/iteration]", "Throughput [iterations/sec]");
rz_table_set_columnsf(T, "sdnnnn", "Benchmark", "Iterations", "Total time [ms]", "Avg iter time [us/iteration]", "Throughput [iterations/sec]", "Std Deviation");
/**
* \brief Prints microbenchmark results and frees the RzTable \p T. Should be called at end of a benchmark suite.

View file

@ -94,6 +94,7 @@ if get_option('enable_tests')
'list',
'log',
'lzma',
'math',
'mem',
'ovf',
'pj',

96
test/unit/test_math.c Normal file
View file

@ -0,0 +1,96 @@
// SPDX-FileCopyrightText: 2026 Rot127 <rot127@posteo.com>
// SPDX-License-Identifier: LGPL-3.0-only
#include <rz_util.h>
#include <rz_vector.h>
#include "minunit.h"
static bool test_welford_arithmetic(void) {
RzMathWelfordSums wf = { 0 };
rz_math_welford_init(&wf);
rz_math_welford_push(&wf, 0.49671415);
rz_math_welford_push(&wf, -0.1382643);
rz_math_welford_push(&wf, 0.64768854);
rz_math_welford_push(&wf, 1.52302986);
rz_math_welford_push(&wf, -0.23415337);
rz_math_welford_push(&wf, -0.23413696);
rz_math_welford_push(&wf, 1.57921282);
rz_math_welford_push(&wf, 0.76743473);
rz_math_welford_push(&wf, -0.46947439);
rz_math_welford_push(&wf, 0.54256004);
rz_math_welford_push(&wf, -0.46341769);
rz_math_welford_push(&wf, -0.46572975);
rz_math_welford_push(&wf, 0.24196227);
rz_math_welford_push(&wf, -1.91328024);
rz_math_welford_push(&wf, -1.72491783);
rz_math_welford_push(&wf, -0.56228753);
rz_math_welford_push(&wf, -1.01283112);
rz_math_welford_push(&wf, 0.31424733);
rz_math_welford_push(&wf, -0.90802408);
rz_math_welford_push(&wf, -1.4123037);
mu_assert_eq(rz_math_welford_n(&wf), 20, "n");
char val[16] = { 0 };
rz_strf(val, "%.6f", rz_math_welford_avar(&wf));
mu_assert_streq(val, "0.875572", "variance");
rz_strf(val, "%.6f", rz_math_welford_astddev(&wf));
mu_assert_streq(val, "0.935720", "std deviation");
rz_strf(val, "%.6f", rz_math_welford_amean(&wf));
mu_assert_streq(val, "-0.171299", "mean");
rz_math_welford_clear(&wf);
mu_end;
}
static bool test_welford_geometric(void) {
RzMathWelfordSums wf = { 0 };
rz_math_welford_init(&wf);
rz_math_welford_push(&wf, 131.5502965);
rz_math_welford_push(&wf, 47.68105836);
rz_math_welford_push(&wf, 56.85572523);
rz_math_welford_push(&wf, 23.22318396);
rz_math_welford_push(&wf, 39.38442231);
rz_math_welford_push(&wf, 58.35549654);
rz_math_welford_push(&wf, 27.36880479);
rz_math_welford_push(&wf, 68.40314554);
rz_math_welford_push(&wf, 38.0772422);
rz_math_welford_push(&wf, 45.8320556);
rz_math_welford_push(&wf, 38.05285189);
rz_math_welford_push(&wf, 165.8969663);
rz_math_welford_push(&wf, 54.15778147);
rz_math_welford_push(&wf, 28.94430432);
rz_math_welford_push(&wf, 89.43632748);
rz_math_welford_push(&wf, 26.24548069);
rz_math_welford_push(&wf, 61.88749606);
rz_math_welford_push(&wf, 16.84742668);
rz_math_welford_push(&wf, 24.60841286);
rz_math_welford_push(&wf, 61.4434194);
mu_assert_eq(rz_math_welford_n(&wf), 20, "n");
char val[16] = { 0 };
rz_strf(val, "%.6f", rz_math_welford_gmean(&wf));
mu_assert_streq(val, "46.544783", "mean");
rz_strf(val, "%.6f", rz_math_welford_gstddev(&wf));
mu_assert_streq(val, "1.787509", "std deviation");
rz_math_welford_clear(&wf);
mu_end;
}
static int all_tests(void) {
mu_run_test(test_welford_arithmetic);
mu_run_test(test_welford_geometric);
return tests_passed != tests_run;
}
mu_main(all_tests)