audio_utils/include/audio_utils/Statistics.h

/*
 * Copyright 2018 The Android Open Source Project
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * 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.
 */

#ifndef ANDROID_AUDIO_UTILS_STATISTICS_H
#define ANDROID_AUDIO_UTILS_STATISTICS_H

#ifdef __cplusplus

#include "variadic_utils.h"

// variadic_utils already contains stl headers; in addition:
#include <deque> // for ReferenceStatistics implementation
#include <sstream>

namespace android {
namespace audio_utils {

/**
 * Compensated summation is used to accumulate a sequence of floating point
 * values, with "compensation" information to help preserve precision lost
 * due to catastrophic cancellation, e.g. (BIG) + (SMALL) - (BIG) = 0.
 *
 * We provide two forms of compensated summation:
 * the Kahan variant (which has better properties if the sum is generally
 * larger than the data added; and the Neumaier variant which is better if
 * the sum or delta may alternatively be larger.
 *
 * https://en.wikipedia.org/wiki/Kahan_summation_algorithm
 *
 * Alternative approaches include divide-and-conquer summation
 * which provides increased accuracy with log n stack depth (recursion).
 *
 * https://en.wikipedia.org/wiki/Pairwise_summation
 */

template <typename T>
struct KahanSum {
    T mSum{};
    T mCorrection{}; // negative low order bits of mSum.

    constexpr KahanSum<T>() = default;

    explicit constexpr KahanSum<T>(const T& value)
        : mSum{value}
    { }

    // takes T not KahanSum<T>
    friend constexpr KahanSum<T> operator+(KahanSum<T> lhs, const T& rhs) {
        const T y = rhs - lhs.mCorrection;
        const T t = lhs.mSum + y;

#ifdef __FAST_MATH__
#warning "fast math enabled, could optimize out KahanSum correction"
#endif

        lhs.mCorrection = (t - lhs.mSum) - y; // compiler, please do not optimize with /fp:fast
        lhs.mSum = t;
        return lhs;
    }

    constexpr KahanSum<T>& operator+=(const T& rhs) { // takes T not KahanSum<T>
        *this = *this + rhs;
        return *this;
    }

    constexpr operator T() const {
        return mSum;
    }

    constexpr void reset() {
        mSum = {};
        mCorrection = {};
    }
};

// A more robust version of Kahan summation for input greater than sum.
// TODO: investigate variants that reincorporate mCorrection bits into mSum if possible.
template <typename T>
struct NeumaierSum {
    T mSum{};
    T mCorrection{}; // low order bits of mSum.

    constexpr NeumaierSum<T>() = default;

    explicit constexpr NeumaierSum<T>(const T& value)
        : mSum{value}
    { }

    friend constexpr NeumaierSum<T> operator+(NeumaierSum<T> lhs, const T& rhs) {
        const T t = lhs.mSum + rhs;

        if (const_abs(lhs.mSum) >= const_abs(rhs)) {
            lhs.mCorrection += (lhs.mSum - t) + rhs;
        } else {
            lhs.mCorrection += (rhs - t) + lhs.mSum;
        }
        lhs.mSum = t;
        return lhs;
    }

    constexpr NeumaierSum<T>& operator+=(const T& rhs) { // takes T not NeumaierSum<T>
        *this = *this + rhs;
        return *this;
    }

    static constexpr T const_abs(T x) {
        return x < T{} ? -x : x;
    }

    constexpr operator T() const {
        return mSum + mCorrection;
    }

    constexpr void reset() {
        mSum = {};
        mCorrection = {};
    }
};

//-------------------------------------------------------------------
// Constants and limits

template <typename T, typename T2=void>  struct StatisticsConstants;

template <typename T>
struct StatisticsConstants<T, std::enable_if_t<std::is_arithmetic<T>::value>> {
    // value closest to negative infinity for type T
    static constexpr T negativeInfinity() {
        return std::numeric_limits<T>::has_infinity ?
                -std::numeric_limits<T>::infinity() : std::numeric_limits<T>::min();
    };

    static constexpr T mNegativeInfinity = negativeInfinity();

    // value closest to positive infinity for type T
    static constexpr T positiveInfinity() {
        return std::numeric_limits<T>::has_infinity ?
                std::numeric_limits<T>::infinity() : std::numeric_limits<T>::max();
    }

    static constexpr T mPositiveInfinity = positiveInfinity();
};

// specialize for tuple and pair
template <typename T>
struct StatisticsConstants<T, std::enable_if_t<!std::is_arithmetic<T>::value>> {
private:
    template <std::size_t... I >
    static constexpr auto negativeInfinity(std::index_sequence<I...>) {
       return T{StatisticsConstants<
               typename std::tuple_element<I, T>::type>::mNegativeInfinity...};
    }
    template <std::size_t... I >
    static constexpr auto positiveInfinity(std::index_sequence<I...>) {
       return T{StatisticsConstants<
               typename std::tuple_element<I, T>::type>::mPositiveInfinity...};
    }
public:
    static constexpr auto negativeInfinity() {
       return negativeInfinity(std::make_index_sequence<std::tuple_size<T>::value>());
    }
    static constexpr auto mNegativeInfinity =
        negativeInfinity(std::make_index_sequence<std::tuple_size<T>::value>());
    static constexpr auto positiveInfinity() {
       return positiveInfinity(std::make_index_sequence<std::tuple_size<T>::value>());
    }
    static constexpr auto mPositiveInfinity =
        positiveInfinity(std::make_index_sequence<std::tuple_size<T>::value>());
};

/**
 * Statistics provides a running weighted average, variance, and standard deviation of
 * a sample stream. It is more numerically stable for floating point computation than a
 * naive sum of values, sum of values squared.
 *
 * The weighting is like an IIR filter, with the most recent sample weighted as 1, and decaying
 * by alpha (between 0 and 1).  With alpha == 1. this is rectangular weighting, reducing to
 * Welford's algorithm.
 *
 * The IIR filter weighting emphasizes more recent samples, has low overhead updates,
 * constant storage, and constant computation (per update or variance read).
 *
 * This is a variant of the weighted mean and variance algorithms described here:
 * https://en.wikipedia.org/wiki/Moving_average
 * https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
 * https://en.wikipedia.org/wiki/Weighted_arithmetic_mean
 *
 * weight = sum_{i=1}^n \alpha^{n-i}
 * mean = 1/weight * sum_{i=1}^n \alpha^{n-i}x_i
 * var = 1/weight * sum_{i=1}^n alpha^{n-i}(x_i- mean)^2
 *
 * The Statistics class is safe to call from a SCHED_FIFO thread with the exception of
 * the toString() method, which uses std::stringstream to format data for printing.
 *
 * Long term data accumulation and constant alpha:
 * If the alpha weight is 1 (or not specified) then statistics objects with float
 * summation types (D, S) should NOT add more than the mantissa-bits elements
 * without reset to prevent variance increases due to weight precision underflow.
 * This is 1 << 23 elements for float and 1 << 52 elements for double.
 *
 * Setting alpha less than 1 avoids this underflow problem.
 * Alpha < 1 - (epsilon * 32), where epsilon is std::numeric_limits<D>::epsilon()
 * is recommended for continuously running statistics (alpha <= 0.999996
 * for float summation precision).
 *
 * Alpha may also change on-the-fly, based on the reliability of
 * new information.  In that case, alpha may be set temporarily greater
 * than 1.
 *
 * https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Reliability_weights_2
 *
 * Statistics may also be collected on variadic "vector" object instead of
 * scalars, where the variance may be computed as an inner product radial squared
 * distance from the mean; or as an outer product where the variance returned
 * is a covariance matrix.
 *
 * TODO:
 * 1) Alternative versions of Kahan/Neumaier sum that better preserve precision.
 * 2) Add binary math ops to corrected sum classes for better precision in lieu of long double.
 * 3) Add Cholesky decomposition to ensure positive definite covariance matrices if
 *    the input is a variadic object.
 */

/**
 * Mean may have catastrophic cancellation of positive and negative sample values,
 * so we use Kahan summation in the algorithms below (or substitute "D" if not needed).
 *
 * https://en.wikipedia.org/wiki/Kahan_summation_algorithm
 */

template <
    typename T,               // input data type
    typename D = double,      // output mean data type
    typename S = KahanSum<D>, // compensated mean summation type, if any
    typename A = double,      // weight type
    typename D2 = double,     // output variance "D^2" type
    typename PRODUCT = std::multiplies<D> // how the output variance is computed
    >
class Statistics {
public:
    /** alpha is the weight (if alpha == 1. we use a rectangular window) */
    explicit constexpr Statistics(A alpha = A(1.))
        : mAlpha(alpha)
    { }

    template <size_t N>
    explicit constexpr Statistics(const T (&a)[N], A alpha = A(1.))
        : mAlpha(alpha)
    {
        for (const auto &data : a) {
            add(data);
        }
    }

    constexpr void setAlpha(A alpha) {
        mAlpha = alpha;
    }

    constexpr void add(const T &value) {
        // Note: fastest implementation uses fmin fminf but would not be constexpr

        mMax = audio_utils::max(mMax, value); // order important: reject NaN
        mMin = audio_utils::min(mMin, value); // order important: reject NaN
        ++mN;
        const D delta = value - mMean;
        /* if (mAlpha == 1.) we have Welford's algorithm
            ++mN;
            mMean += delta / mN;
            mM2 += delta * (value - mMean);

            Note delta * (value - mMean) should be non-negative.
        */
        mWeight = A(1.) + mAlpha * mWeight;
        mWeight2 = A(1.) + mAlpha * mAlpha * mWeight2;
        D meanDelta = delta / mWeight;
        mMean += meanDelta;
        mM2 = mAlpha * mM2 + PRODUCT()(delta, (value - mMean));

        /*
           Alternate variant related to:
           http://mathworld.wolfram.com/SampleVarianceComputation.html

           const double sweight = mAlpha * mWeight;
           mWeight = 1. + sweight;
           const double dmean = delta / mWeight;
           mMean += dmean;
           mM2 = mAlpha * mM2 + mWeight * sweight * dmean * dmean;

           The update is slightly different than Welford's algorithm
           showing a by-construction non-negative update to M2.
        */
    }

    constexpr int64_t getN() const {
        return mN;
    }

    constexpr void reset() {
        mMin = StatisticsConstants<T>::positiveInfinity();
        mMax = StatisticsConstants<T>::negativeInfinity();
        mN = 0;
        mWeight = {};
        mWeight2 = {};
        mMean = {};
        mM2 = {};
    }

    constexpr A getWeight() const {
        return mWeight;
    }

    constexpr D getMean() const {
        return mMean;
    }

    constexpr D2 getVariance() const {
        if (mN < 2) {
            // must have 2 samples for sample variance.
            return {};
        } else {
            return mM2 / getSampleWeight();
        }
    }

    constexpr D2 getPopVariance() const {
        if (mN < 1) {
            return {};
        } else {
            return mM2 / mWeight;
        }
    }

    // explicitly use sqrt_constexpr if you need a constexpr version
    D2 getStdDev() const {
        return android::audio_utils::sqrt(getVariance());
    }

    D2 getPopStdDev() const {
        return android::audio_utils::sqrt(getPopVariance());
    }

    constexpr T getMin() const {
        return mMin;
    }

    constexpr T getMax() const {
        return mMax;
    }

    std::string toString() const {
        const int64_t N = getN();
        if (N == 0) return "unavail";

        std::stringstream ss;
        ss << "ave=" << getMean();
        if (N > 1) {
            // we use the sample standard deviation (not entirely unbiased,
            // though the sample variance is unbiased).
            ss << " std=" << getStdDev();
        }
        ss << " min=" << getMin();
        ss << " max=" << getMax();
        return ss.str();
    }

private:
    A mAlpha;
    T mMin{StatisticsConstants<T>::positiveInfinity()};
    T mMax{StatisticsConstants<T>::negativeInfinity()};

    int64_t mN = 0;  // running count of samples.
    A mWeight{};     // sum of weights.
    A mWeight2{};    // sum of weights squared.
    S mMean{};       // running mean.
    D2 mM2{};         // running unnormalized variance.

    // Reliability correction for unbiasing variance, since mean is estimated
    // from same sample stream as variance.
    // if mAlpha == 1 this is mWeight - 1;
    //
    // TODO: consider exposing the correction factor.
    constexpr A getSampleWeight() const {
        // if mAlpha is constant then the mWeight2 member variable is not
        // needed, one can use instead:
        // return (mWeight - D(1.)) * D(2.) / (D(1.) + mAlpha);

        return mWeight - mWeight2 / mWeight;
    }
};

/**
 * ReferenceStatistics is a naive implementation of the weighted running variance,
 * which consumes more space and is slower than Statistics.  It is provided for
 * comparison and testing purposes.  Do not call from a SCHED_FIFO thread!
 *
 * Note: Common code not combined for implementation clarity.
 *       We don't invoke Kahan summation or other tricks.
 */
template <
    typename T, // input data type
    typename D = double // output mean/variance data type
    >
class ReferenceStatistics {
public:
    /** alpha is the weight (alpha == 1. is rectangular window) */
    explicit ReferenceStatistics(D alpha = D(1.))
        : mAlpha(alpha)
    { }

    constexpr void setAlpha(D alpha) {
        mAlpha = alpha;
    }

    // For independent testing, have intentionally slightly different behavior
    // of min and max than Statistics with respect to Nan.
    constexpr void add(const T &value) {
        if (getN() == 0) {
            mMax = value;
            mMin = value;
        } else if (value > mMax) {
            mMax = value;
        } else if (value < mMin) {
            mMin = value;
        }

        mData.push_front(value);
        mAlphaList.push_front(mAlpha);
    }

    int64_t getN() const {
        return mData.size();
    }

    void reset() {
        mMin = {};
        mMax = {};
        mData.clear();
        mAlphaList.clear();
    }

    D getWeight() const {
        D weight{};
        D alpha_i(1.);
        for (size_t i = 0; i < mData.size(); ++i) {
            weight += alpha_i;
            alpha_i *= mAlphaList[i];
        }
        return weight;
    }

    D getWeight2() const {
        D weight2{};
        D alpha2_i(1.);
        for (size_t i = 0; i < mData.size(); ++i) {
            weight2 += alpha2_i;
            alpha2_i *= mAlphaList[i] * mAlphaList[i];
        }
        return weight2;
    }

    D getMean() const {
        D wsum{};
        D alpha_i(1.);
        for (size_t i = 0; i < mData.size(); ++i) {
            wsum += alpha_i * mData[i];
            alpha_i *= mAlphaList[i];
        }
        return wsum / getWeight();
    }

    // Should always return a positive value.
    D getVariance() const {
        return getUnweightedVariance() / (getWeight() - getWeight2() / getWeight());
    }

    // Should always return a positive value.
    D getPopVariance() const {
        return getUnweightedVariance() / getWeight();
    }

    D getStdDev() const {
        return sqrt(getVariance());
    }

    D getPopStdDev() const {
        return sqrt(getPopVariance());
    }

    T getMin() const {
        return mMin;
    }

    T getMax() const {
        return mMax;
    }

    std::string toString() const {
        const auto N = getN();
        if (N == 0) return "unavail";

        std::stringstream ss;
        ss << "ave=" << getMean();
        if (N > 1) {
            // we use the sample standard deviation (not entirely unbiased,
            // though the sample variance is unbiased).
            ss << " std=" << getStdDev();
        }
        ss << " min=" << getMin();
        ss << " max=" << getMax();
        return ss.str();
    }

private:
    T mMin{};
    T mMax{};

    D mAlpha;                 // current alpha value
    std::deque<T> mData;      // store all the data for exact summation, mData[0] most recent.
    std::deque<D> mAlphaList; // alpha value for the data added.

    D getUnweightedVariance() const {
        const D mean = getMean();
        D wsum{};
        D alpha_i(1.);
        for (size_t i = 0; i < mData.size(); ++i) {
            const D diff = mData[i] - mean;
            wsum += alpha_i * diff * diff;
            alpha_i *= mAlphaList[i];
        }
        return wsum;
    }
};

/**
 * Least squares fitting of a 2D straight line based on the covariance matrix.
 *
 * See formula from:
 * http://mathworld.wolfram.com/LeastSquaresFitting.html
 *
 * y = a + b*x
 *
 * returns a: y intercept
 *         b: slope
 *         r2: correlation coefficient (1.0 means great fit, 0.0 means no fit.)
 *
 * For better numerical stability, it is suggested to use the slope b only:
 * as the least squares fit line intersects the mean.
 *
 * (y - mean_y) = b * (x - mean_x).
 *
 */
template <typename T>
constexpr void computeYLineFromStatistics(
        T &a, T& b, T &r2,
        const T& mean_x,
        const T& mean_y,
        const T& var_x,
        const T& cov_xy,
        const T& var_y) {

    // Dimensionally r2 is unitless.  If there is no correlation
    // then r2 is clearly 0 as cov_xy == 0.  If x and y are identical up to a scale
    // and shift, then r2 is 1.
    r2 = cov_xy * cov_xy / (var_x * var_y);

    // The least squares solution to the overconstrained matrix equation requires
    // the pseudo-inverse. In 2D, the best-fit slope is the mean removed
    // (via covariance and variance) dy/dx derived from the joint expectation
    // (this is also dimensionally correct).
    b = cov_xy / var_x;

    // The best fit line goes through the mean, and can be used to find the y intercept.
    a = mean_y - b * mean_x;
}

/**
 * LinearLeastSquaresFit<> class is derived from the Statistics<> class, with a 2 element array.
 * Arrays are preferred over tuples or pairs because copy assignment is constexpr and
 * arrays are trivially copyable.
 */
template <typename T>
class LinearLeastSquaresFit : public
    Statistics<std::array<T, 2>, // input
               std::array<T, 2>, // mean data output
               std::array<T, 2>, // compensated mean sum
               T,                // weight type
               std::array<T, 3>, // covariance_ut
               audio_utils::outerProduct_UT_array<std::array<T, 2>>>
{
public:
    constexpr explicit LinearLeastSquaresFit(const T &alpha = T(1.))
        : Statistics<std::array<T, 2>,
             std::array<T, 2>,
             std::array<T, 2>,
             T,
             std::array<T, 3>, // covariance_ut
             audio_utils::outerProduct_UT_array<std::array<T, 2>>>(alpha) { }

    /* Note: base class method: add(value)

    constexpr void add(const std::array<T, 2>& value);

       use:
          add({1., 2.});
       or
          add(to_array(myTuple));
    */

    /**
     * y = a + b*x
     *
     * returns a: y intercept
     *         b: y slope (dy / dx)
     *         r2: correlation coefficient (1.0 means great fit, 0.0 means no fit.)
     */
    constexpr void computeYLine(T &a, T &b, T &r2) const {
        computeYLineFromStatistics(a, b, r2,
                std::get<0>(this->getMean()), /* mean_x */
                std::get<1>(this->getMean()), /* mean_y */
                std::get<0>(this->getPopVariance()), /* var_x */
                std::get<1>(this->getPopVariance()), /* cov_xy */
                std::get<2>(this->getPopVariance())); /* var_y */
    }

    /**
     * x = a + b*y
     *
     * returns a: x intercept
     *         b: x slope (dx / dy)
     *         r2: correlation coefficient (1.0 means great fit, 0.0 means no fit.)
     */
    constexpr void computeXLine(T &a, T &b, T &r2) const {
        // reverse x and y for X line computation
        computeYLineFromStatistics(a, b, r2,
                std::get<1>(this->getMean()), /* mean_x */
                std::get<0>(this->getMean()), /* mean_y */
                std::get<2>(this->getPopVariance()), /* var_x */
                std::get<1>(this->getPopVariance()), /* cov_xy */
                std::get<0>(this->getPopVariance())); /* var_y */
    }

    /**
     * this returns the estimate of y from a given x
     */
    constexpr T getYFromX(const T &x) const {
        const T var_x = std::get<0>(this->getPopVariance());
        const T cov_xy = std::get<1>(this->getPopVariance());
        const T b = cov_xy / var_x;  // dy / dx

        const T mean_x = std::get<0>(this->getMean());
        const T mean_y = std::get<1>(this->getMean());
        return /* y = */ b * (x - mean_x) + mean_y;
    }

    /**
     * this returns the estimate of x from a given y
     */
    constexpr T getXFromY(const T &y) const {
        const T cov_xy = std::get<1>(this->getPopVariance());
        const T var_y = std::get<2>(this->getPopVariance());
        const T b = cov_xy / var_y;  // dx / dy

        const T mean_x = std::get<0>(this->getMean());
        const T mean_y = std::get<1>(this->getMean());
        return /* x = */ b * (y - mean_y) + mean_x;
    }

    constexpr T getR2() const {
        const T var_x = std::get<0>(this->getPopVariance());
        const T cov_xy = std::get<1>(this->getPopVariance());
        const T var_y = std::get<2>(this->getPopVariance());
        return cov_xy * cov_xy / (var_x * var_y);
    }
};

/**
 * constexpr statistics functions of form:
 * algorithm(forward_iterator begin, forward_iterator end)
 *
 * These check that the input looks like an iterator, but doesn't
 * check if __is_forward_iterator<>.
 *
 * divide-and-conquer pairwise summation forms will require
 * __is_random_access_iterator<>.
 */

// returns max of elements, or if no elements negative infinity.
template <typename T,
          std::enable_if_t<is_iterator<T>::value, int> = 0>
constexpr auto max(T begin, T end) {
    using S = std::remove_cv_t<std::remove_reference_t<
            decltype(*begin)>>;
    S maxValue = StatisticsConstants<S>::mNegativeInfinity;
    for (auto it = begin; it != end; ++it) {
        maxValue = std::max(maxValue, *it);
    }
    return maxValue;
}

// returns min of elements, or if no elements positive infinity.
template <typename T,
          std::enable_if_t<is_iterator<T>::value, int> = 0>
constexpr auto min(T begin, T end) {
    using S = std::remove_cv_t<std::remove_reference_t<
            decltype(*begin)>>;
    S minValue = StatisticsConstants<S>::mPositiveInfinity;
    for (auto it = begin; it != end; ++it) {
        minValue = std::min(minValue, *it);
    }
    return minValue;
}

template <typename D = double, typename S = KahanSum<D>, typename T,
          std::enable_if_t<is_iterator<T>::value, int> = 0>
constexpr auto sum(T begin, T end) {
    S sum{};
    for (auto it = begin; it != end; ++it) {
        sum += D(*it);
    }
    return sum;
}

template <typename D = double, typename S = KahanSum<D>, typename T,
          std::enable_if_t<is_iterator<T>::value, int> = 0>
constexpr auto sumSqDiff(T begin, T end, D x = {}) {
    S sum{};
    for (auto it = begin; it != end; ++it) {
        const D diff = *it - x;
        sum += diff * diff;
    }
    return sum;
}

// Form: algorithm(array[]), where array size is known to the compiler.
template <typename T, size_t N>
constexpr T max(const T (&a)[N]) {
    return max(&a[0], &a[N]);
}

template <typename T, size_t N>
constexpr T min(const T (&a)[N]) {
    return min(&a[0], &a[N]);
}

template <typename D = double, typename S = KahanSum<D>, typename T, size_t N>
constexpr D sum(const T (&a)[N]) {
    return sum<D, S>(&a[0], &a[N]);
}

template <typename D = double, typename S = KahanSum<D>, typename T, size_t N>
constexpr D sumSqDiff(const T (&a)[N], D x = {}) {
    return sumSqDiff<D, S>(&a[0], &a[N], x);
}

// TODO: remove when std::isnan is constexpr
template <typename T>
constexpr T isnan(T x) {
    return __builtin_isnan(x);
}

// constexpr sqrt computed by the Babylonian (Newton's) method.
// Please use math libraries for non-constexpr cases.
// TODO: remove when there is some std::sqrt which is constexpr.
//
// https://en.wikipedia.org/wiki/Methods_of_computing_square_roots

// watch out using the unchecked version, use the checked version below.
template <typename T>
constexpr T sqrt_constexpr_unchecked(T x, T prev) {
    static_assert(std::is_floating_point<T>::value, "must be floating point type");
    const T next = T(0.5) * (prev + x / prev);
    return next == prev ? next : sqrt_constexpr_unchecked(x, next);
}

// checked sqrt
template <typename T>
constexpr T sqrt_constexpr(T x) {
    static_assert(std::is_floating_point<T>::value, "must be floating point type");
    if (x < T{}) { // negative values return nan
        return std::numeric_limits<T>::quiet_NaN();
    } else if (isnan(x)
            || x == std::numeric_limits<T>::infinity()
            || x == T{}) {
        return x;
    } else { // good to go.
        return sqrt_constexpr_unchecked(x, T(1.));
    }
}

} // namespace audio_utils
} // namespace android

#endif // __cplusplus

/** \cond */
 __BEGIN_DECLS
/** \endcond */

/** Simple stats structure for low overhead statistics gathering.
 * Designed to be accessed by C (with no functional getters).
 * Zero initialize {} to clear or reset.
 */
typedef struct {
   int64_t n;
   double min;
   double max;
   double last;
   double mean;
} simple_stats_t;

/** logs new value to the simple_stats_t */
static inline void simple_stats_log(simple_stats_t *stats, double value) {
    if (++stats->n == 1) {
        stats->min = stats->max = stats->last = stats->mean = value;
    } else {
        stats->last = value;
        if (value < stats->min) {
            stats->min = value;
        } else if (value > stats->max) {
            stats->max = value;
        }
        // Welford's algorithm for mean
        const double delta = value - stats->mean;
        stats->mean += delta / stats->n;
    }
}

/** dumps statistics to a string, returns the length of string excluding null termination. */
static inline size_t simple_stats_to_string(simple_stats_t *stats, char *buffer, size_t size) {
    if (size == 0) {
        return 0;
    } else if (stats->n == 0) {
        return snprintf(buffer, size, "none");
    } else {
        return snprintf(buffer, size, "(mean: %lf  min: %lf  max: %lf  last: %lf  n: %lld)",
                stats->mean, stats->min, stats->max, stats->last, (long long)stats->n);
    }
}

/** \cond */
__END_DECLS
/** \endcond */

#endif // !ANDROID_AUDIO_UTILS_STATISTICS_H