diff --git a/yoshimi-testrunner/src/util/statistic.cpp b/yoshimi-testrunner/src/util/statistic.cpp deleted file mode 100644 index 22c3d21f0..000000000 --- a/yoshimi-testrunner/src/util/statistic.cpp +++ /dev/null @@ -1,57 +0,0 @@ -/* - * statistic - helpers for generic statistics calculations - * - * Copyright 2021, Hermann Vosseler - * - * This file is part of the Yoshimi-Testsuite, which is free software: - * you can redistribute and/or modify it under the terms of the GNU - * General Public License as published by the Free Software Foundation, - * either version 3 of the License, or (at your option) any later version. - * - * Yoshimi-Testsuite is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with yoshimi. If not, see . - ***************************************************************/ - - -/** @file statistic.cpp - ** Implementation details for generic statistics calculations. - ** - ** @todo WIP as of 9/21 - ** - */ - - -#include "util/statistic.hpp" -//#include "util/format.hpp" -#include "util/error.hpp" - -//#include -//#include -//#include - - -namespace util { - -//using util::formatVal; -//using std::log10; -//using std::fabs; -//using std::max; -//using std::min; - - -namespace { // Implementation details -}//(End)Implementation namespace - - - - -/** - */ - - -}//(End)namespace util diff --git a/yoshimi-testrunner/src/util/statistic.hpp b/yoshimi-testrunner/src/util/statistic.hpp index c71a3dfa2..9fe37e6e2 100644 --- a/yoshimi-testrunner/src/util/statistic.hpp +++ b/yoshimi-testrunner/src/util/statistic.hpp @@ -20,9 +20,10 @@ /** @file statistic.cpp ** Support for generic statistics calculations. - ** - ** @todo WIP as of 9/21 - ** + ** - average over the N last elements in a data sequence + ** - simple linear regression with weights (single predictor variable) + ** - also over a time series with zero-based indices + ** */ @@ -32,9 +33,10 @@ #include "util/error.hpp" +#include "util/nocopy.hpp" +#include "util/format.hpp" +#include "util/utils.hpp" -//#include -//#include #include #include #include @@ -42,26 +44,110 @@ namespace util { -//using util::formatVal; -//using std::log10; -//using std::fabs; -//using std::max; -//using std::min; +using std::fabs; +using std::array; +using std::tuple; +using std::make_tuple; + + +/** + * Read-only view into a segment within a sequence of data + * @tparam D value type of the data series + * @remark simplistic workaround since we don't support C++20 yet + * @todo replace by const std::span + */ +template +class DataSpan + : util::Cloneable +{ + const D* const b_{nullptr}; + const D* const e_{nullptr}; + +public: + DataSpan() = default; + DataSpan(D const& begin, D const& end) + : b_{&begin} + , e_{&end} + { if (e_ < b_) throw error::Invalid("End point before begin."); } + + template + DataSpan(CON const& container) + : DataSpan{*std::begin(container), *std::end(container)} + { } + + + using iterator = const D*; + + size_t size() const { return e_ - b_; } + bool empty() const { return b_ == e_; } + + iterator begin() const { return b_; } + iterator end() const { return e_; } + + D const& operator[](size_t i) const { return b_ + i; } + D const& at(size_t i) const + { + if (i >= size()) throw error::Invalid("Index "+str(i)+" beyond size="+str(size())); + return this->operator[](i); + } +}; + + +template +using Span = DataSpan const&; using VecD = std::vector; -namespace { // Implementation details -}//(End)Implementation namespace +/** helper to unpack a std::tuple into a homogeneous std::array */ +template +constexpr auto array_from_tuple(TUP&& tuple) +{ + constexpr auto makeArray = [](auto&& ... x){ return std::array{std::forward(x) ... }; }; + return std::apply(makeArray, std::forward(tuple)); +} -inline double averageLastN(VecD const& data, size_t n) + +template +inline double average(Span data) +{ + if (isnil(data)) return 0.0; + double sum = 0.0; + for (auto val : data) + sum += val; + return sum / data.size(); +} + + +inline DataSpan lastN(VecD const& data, size_t n) { n = std::min(n, data.size()); size_t oldest = data.size() - n; - double sum = 0.0; - for (size_t i=data.size(); oldest < i; --i) - sum += data[i-1]; - return 0{data[oldest], *data.end()}; +} + +inline double averageLastN(VecD const& data, size_t n) +{ + return average(lastN(data,n)); +} + + +/** "building blocks" for mean, variance and covariance of time series data */ +template +inline auto computeStatSums(Span series) +{ + double ysum = 0.0; + double yysum = 0.0; + double xysum = 0.0; + size_t x = 0; + for (auto& y : series) + { + ysum += y; + yysum += y*y; + xysum += x*y; + ++x; + } + return make_tuple(ysum,yysum, xysum); } @@ -81,9 +167,10 @@ using RegressionData = std::vector; /** "building blocks" for weighted mean, weighted variance and covariance */ -inline auto computeWeightedStatSums(RegressionData const& points) +inline auto computeWeightedStatSums(Span points) { - std::array sums = {0.0,0.0,0.0,0.0,0.0,0.0}; + std::array sums; + sums.fill(0.0); auto& [wsum, wxsum, wysum, wxxsum, wyysum, wxysum] = sums; for (auto& p : points) { @@ -99,7 +186,7 @@ inline auto computeWeightedStatSums(RegressionData const& points) /** * Compute simple linear regression with a single predictor variable (x). - * @param points data to fit the linear model with + * @param points 2D data to fit the linear model with, including weights. * @return the computed linear model `b + a·x`, and the resulting fit * - socket (constant offset `b`) * - gradient (linear factor `a`) @@ -109,7 +196,7 @@ inline auto computeWeightedStatSums(RegressionData const& points) * - maximum absolute delta * - delta standard deviation */ -inline std::tuple computeLinearRegression(RegressionData const& points) +inline auto computeLinearRegression(Span points) { auto [wsum, wxsum, wysum, wxxsum, wyysum, wxysum] = computeWeightedStatSums(points); @@ -140,15 +227,54 @@ inline std::tuple computeLinearR maxDelta = std::max(maxDelta, fabs(delta)); variance += p.w * delta*delta; } - return std::make_tuple(socket,gradient - ,move(predicted) - ,move(deltas) - ,correlation - ,maxDelta - ,sqrt(variance/wsum) - ); + return make_tuple(socket,gradient + ,move(predicted) + ,move(deltas) + ,correlation + ,maxDelta + ,sqrt(variance/wsum) + ); } +inline auto computeLinearRegression(RegressionData const& points) +{ return computeLinearRegression(DataSpan{points}); } + + + +/** + * Compute linear regression over a time series with zero-based indices. + * @remark using the indices as x-values, the calculations for a regression line + * can be simplified, using the known closed formula for a sum of integers, + * shifting the indices to 0…n-1 (leaving out the 0 and 0² term) + * - `1+…+n = n·(n+1)/2` + * - `1+…+n² = n·(n+1)·(2n+1)/6` + * @return `(socket,gradient)` to describe the regression line y = socket + gradient · i + */ +template +inline auto computeTimeSeriesLinearRegression(Span series) +{ + auto [ysum,yysum, xysum] = computeStatSums(series); + + size_t n = series.size(); + double im = (n-1)/2.0; // mean of zero-based indices i ∈ {0 … n-1} + double ym = ysum / n; // mean y + double varx = (n-1)*(n+1)/12.0; // variance of zero-based indices Σ(i-im)² / n + double vary = yysum/n - ym*ym; // variance of data values Σ(y-ym)² / n + double cova = xysum - ysum *(n-1)/2; // Time series Covariance = Σ(i-im)(y-ym) = Σiy + im·ym·n - ym·Σi - im·Σy; use n*ym = Σy + + // Linear Regression minimising σ² + double gradient = cova / (n*varx); // Gradient = Correlation · σy / σx ; σ = √Variance; Correlation = Covariance /(√Σx √Σy) + double socket = ym - gradient * im; // Regression line: Y-ym = Gradient · (i-im) ; set i≔0 yields socket + + // Correlation (Pearson's r) + double correlation = gradient * sqrt(varx/vary); + return make_tuple(socket,gradient,correlation); +} + +inline auto computeTimeSeriesLinearRegression(VecD const& series) +{ return computeTimeSeriesLinearRegression(Span{series}); } + + }//(End)namespace util #endif /*TESTRUNNER_UTIL_STATISTIC_HPP_*/