statistics.h

Go to the documentation of this file.

#pragma once

#include <cstddef>
#include <limits>
#include <numeric>
#include <algorithm>
#include "interop/util/assert.h"
#include "interop/util/math.h"

namespace illumina { namespace interop { namespace util
{

    namespace op
    {
        struct parameter_none_type
        {
        };

        template<class T, typename R, typename P1=parameter_none_type>
        struct const_member_function_w
        {
            const_member_function_w(P1 param1, R (T::*func )(P1) const) : m_param1(param1), m_function(func)
            { }

            template<class F>
            F operator()(const F val, const T &obj) const
            {
                return val + operator()(obj);
            }
            float operator()(const float val, const T &obj) const
            {
                const float ret = static_cast<float>(operator()(obj));
                if(std::isnan(ret)) return val;
                return val + ret;
            }
            double operator()(const double val, const T &obj) const
            {
                const double ret = static_cast<double>(operator()(obj));
                if(std::isnan(ret)) return val;
                return val + ret;
            }

            R operator()(const T &obj) const
            {
                return (obj.*m_function)(m_param1);
            }

        private:
            P1 m_param1;

            R (T::*m_function )(P1) const;
        };

        template<class T, typename R>
        struct const_member_function_w<T, R, parameter_none_type>
        {
            const_member_function_w(R (T::*func )() const) : m_function(func)
            { }

            template<class F>
            F operator()(const F val, const T &obj) const
            {
                const F ret = static_cast<F>(operator()(obj));
                return val + ret;
            }
            float operator()(const float val, const T &obj) const
            {
                const float ret = static_cast<float>(operator()(obj));
                if(std::isnan(ret)) return val;
                return val + ret;
            }
            double operator()(const double val, const T &obj) const
            {
                const double ret = static_cast<double>(operator()(obj));
                if(std::isnan(ret)) return val;
                return val + ret;
            }

            R operator()(const T &obj) const
            {
                return (obj.*m_function)();
            }

        private:
            R (T::*m_function )() const;
        };

        template<class T, typename R, typename P1=parameter_none_type>
        struct const_member_function_less_w
        {
            const_member_function_less_w(const const_member_function_w<T, R, P1> &func) : m_func(func)
            { }

            bool operator()(const T &lhs, const T &rhs) const
            {
                return m_func(lhs) < m_func(rhs);
            }

        private:
            const_member_function_w<T, R, P1> m_func;
        };
        struct dummy_arg { dummy_arg(){}};

        template<class T, typename R, typename P2, typename P1>
        const_member_function_w<T, R, P1> const_member_function(const P2& param1, R (T::*func )(P1) const)
        {
            return const_member_function_w<T, R, P1>(param1, func);
        }

        template<class T, typename R>
        const_member_function_w<T, R> const_member_function(R (T::*func )() const)
        {
            return const_member_function_w<T, R>(func);
        }


        template<class T, typename R>
        const_member_function_w<T, R> const_member_function(dummy_arg, R (T::*func )() const)
        {
            return const_member_function_w<T, R>(func);
        }

        template<class T, typename R, typename P2, typename P1>
        const_member_function_less_w<T, R, P1> const_member_function_less(const P2& param1, R (T::*func )(P1) const)
        {
            return const_member_function_less_w<T, R, P1>(const_member_function(param1, func));
        }

        template<class T, typename R>
        const_member_function_less_w<T, R> const_member_function_less(R (T::*func )() const)
        {
            return const_member_function_less_w<T, R>(const_member_function(func));
        }

        struct operator_none
        {
            template<typename F, typename T>
            F operator()(const F val1, const T &val2)
            {
                return static_cast<F>(val1 + val2);
            }

            template<class T>
            const T &operator()(const T &val)
            {
                return val;
            }
        };

        template<typename UnaryOp>
        struct nan_check
        {
            nan_check(const UnaryOp &op) : m_op(op)
            { }

            template<class T>
            size_t operator()(const T &obj) const
            {
                return !std::isnan(m_op(obj));
            }

        private:
            UnaryOp m_op;
        };
    }

    template<typename I, typename O, typename F>
    void outliers_lower(I beg, I end, const F bound, O out)
    {
        for (; beg != end; ++beg)
        {
            if (*beg < bound)
            {
                *out = *beg;
                ++out;
            }
            else break;
        }
    }

    template<typename I, typename O, typename F>
    void outliers_upper(I beg, I end, const F bound, O out)
    {
        if (end == beg) return;
        for (--end; beg != end; --end)
        {
            if (*end > bound)
            {
                *out = *end;
                ++out;
            }
            else return;
        }
        if (beg == end) return;
        *out = *end;
        ++out;

    }

    template<typename I>
    I percentile(I beg, I end, const size_t percentile)
    {
        INTEROP_ASSERT(percentile > 0 && percentile <= 100);
        const size_t n = static_cast<size_t>(std::distance(beg, end));
        const size_t nth_index = static_cast<size_t>(std::ceil(percentile * n / 100.0) - 1);
        I nth_element_iterator = beg + nth_index;
        std::nth_element(beg, nth_element_iterator, end);
        return nth_element_iterator;
    }

    template<typename I, typename Compare>
    I percentile(I beg, I end, const size_t percentile, Compare comp)
    {
        INTEROP_ASSERT(percentile > 0 && percentile <= 100);
        const size_t n = static_cast<size_t>(std::distance(beg, end));
        const size_t nth_index = static_cast<size_t>(std::ceil(percentile * n / 100.0) - 1);
        I nth_element_iterator = beg + nth_index;
        std::nth_element(beg, nth_element_iterator, end, comp);
        return nth_element_iterator;
    }

    template<typename F>
    F interpolate_linear(const F y1, const F y2, const F x1, const F x2, const F xt)
    {
        return y1 + (y2 - y1) / (x2 - x1) * (xt - x1);
    }

    template<typename F, typename I>
    F percentile_sorted(I beg, I end, const size_t percentile)
    {
        INTEROP_ASSERT(beg != end);
        INTEROP_ASSERT(percentile > 0 && percentile <= 100);
        const size_t n = static_cast<size_t>(std::distance(beg, end));
        if (n == 0) return std::numeric_limits<F>::quiet_NaN();
        size_t nth_index = percentile * n / 100;
        if ((n * percentile / 100.0f - nth_index) < 0.5f)
        {
            if (nth_index == 0) return *beg;
            nth_index--;
        }
        if (nth_index >= (n - 1)) return *(end - 1);
        const F y1 = *(beg + nth_index);
        const F y2 = *(beg + nth_index + 1);
        const F x1 = 100.0f * (nth_index + 0.5f) / n;
        const F x2 = 100.0f * (nth_index + 0.5f + 1) / n;
        return interpolate_linear(y1, y2, x1, x2, static_cast<float>(percentile));
    }
    template<typename F, typename I, typename Op>
    F percentile_sorted(I beg, I end, const size_t percentile, Op op)
    {
        INTEROP_ASSERT(beg != end);
        INTEROP_ASSERT(percentile > 0 && percentile <= 100);
        const size_t n = static_cast<size_t>(std::distance(beg, end));
        if (n == 0) return std::numeric_limits<F>::quiet_NaN();
        size_t nth_index = percentile * n / 100;
        if ((n * percentile / 100.0f - nth_index) < 0.5f)
        {
            if (nth_index == 0) return op(*beg);
            nth_index--;
        }
        if (nth_index >= (n - 1)) return op(*(end - 1));
        const F y1 = op(*(beg + nth_index));
        const F y2 = op(*(beg + nth_index + 1));
        const F x1 = 100.0f * (nth_index + 0.5f) / n;
        const F x2 = 100.0f * (nth_index + 0.5f + 1) / n;
        return interpolate_linear(y1, y2, x1, x2, static_cast<float>(percentile));
    }

// TODO: median using nth_element
    template<typename I, typename UnaryOp>
    I remove_nan(I beg, I end, UnaryOp op)
    {
        return std::stable_partition(beg, end, op::nan_check<UnaryOp>(op));
    }

    template<typename I>
    I median(I beg, I end)
    {
        return percentile(beg, end, 50);
    }

    template<typename I, typename Compare>
    I median(I beg, I end, Compare comp)
    {
        return percentile(beg, end, 50, comp);
    }
    template<typename F, typename I>
    F median_interpolated(I beg, I end)
    {
        std::stable_sort(beg, end);
        return percentile_sorted<F>(beg, end, 50);
    }
    template<typename F, typename I, typename Compare>
    F median_interpolated(I beg, I end, Compare comp)
    {
        std::stable_sort(beg, end, comp);
        return percentile_sorted<F>(beg, end, 50);
    }
    template<typename F, typename I, typename Compare, typename Op>
    F median_interpolated(I beg, I end, Compare comp, Op op)
    {
        std::stable_sort(beg, end, comp);
        return percentile_sorted<F>(beg, end, 50, op);
    }

    //TODO: remove_nan
    // TODO:  nan_median using nth_element

    template<typename R, typename I, typename BinaryOp>
    R nan_mean(I beg, I end, BinaryOp op)
    {
        ptrdiff_t n = 0;
        R sum = 0;
        for (; beg != end; ++beg)
        {
            R val = op(*beg);
            if (std::isnan(val)) continue;
            sum += val;
            ++n;
        }
        if (n == 0)return std::numeric_limits<R>::quiet_NaN();
        return sum / n;
    }

    template<typename R, typename I, typename BinaryOp>
    R nan_variance_with_mean(I beg, I end, const R mean_val, BinaryOp op)
    {
        ptrdiff_t n = 0;
        R sum2 = 0;
        R sum3 = 0;
        for (; beg != end; ++beg)
        {
            const R val = op(*beg) - mean_val;
            if (std::isnan(val)) continue;
            sum2 += val * val;
            sum3 += val;
            ++n;
        }
        if (n <= 1) return std::numeric_limits<R>::quiet_NaN();
        return (sum2 - sum3 * sum3 / n) / (n - 1);
    }

    template<typename R, typename I, typename BinaryOp>
    R nan_variance(I beg, I end, BinaryOp op)
    {
        const R mean_val = nan_mean<R>(beg, end, op);
        return nan_variance_with_mean<R>(beg, end, mean_val, op);
    }

    template<typename R, typename I, typename BinaryOp>
    R mean(I beg, I end, BinaryOp op)
    {
        ptrdiff_t n = std::distance(beg, end);
        if (n == 0) return 0;
        return std::accumulate(beg, end, R(0), op) / R(n);
    }

    template<typename R, typename I, typename BinaryOp>
    R variance_with_mean(I beg, I end, const R mean_val, BinaryOp op)
    {
        ptrdiff_t n = std::distance(beg, end);
        R sum2 = 0;
        R sum3 = 0;
        for (; beg != end; ++beg)
        {
            const R val = op(*beg) - mean_val;
            sum2 += val * val;
            sum3 += val;
        }
        if (n <= 1) return 0;
        return (sum2 - sum3 * sum3 / n) / (n - 1);
    }

    template<typename R, typename I, typename BinaryOp>
    R variance(I beg, I end, BinaryOp op)
    {
        const R mean_val = mean<R>(beg, end, op);
        return variance_with_mean<R>(beg, end, mean_val, op);
    }


    template<typename R, typename I>
    R mean(I beg, I end)
    {
        return mean<R>(beg, end, op::operator_none());
    }

    template<typename R, typename I>
    R variance(I beg, I end)
    {
        return variance<R>(beg, end, op::operator_none());
    }

    template<typename R, typename I>
    R variance_with_mean(I beg, I end, const R mean)
    {
        return variance_with_mean<R>(beg, end, mean, op::operator_none());
    }

}}}