summary_statistics.h

Go to the documentation of this file.

#pragma once

#include "interop/model/model_exceptions.h"
#include "interop/util/statistics.h"
#include "interop/model/summary/run_summary.h"


namespace illumina { namespace interop { namespace logic { namespace summary
{
    template<typename I>
    void reserve(I beg, I end, const ptrdiff_t n)
    {
        for (; beg != end; ++beg) beg->reserve(n);
    }

    template<typename T>
    class summary_by_lane_read
    {
    public:
        typedef std::vector<T> vector_t; // TODO: Speed up with 2d vector?
        typedef std::vector<vector_t> vector2d_t;
        typedef std::vector<vector2d_t> vector3d_t;
        typedef typename vector_t::const_iterator const_iterator;
        typedef typename vector3d_t::iterator read_iterator;
        typedef typename vector_t::const_reference const_reference;


    public:
        summary_by_lane_read(const model::summary::run_summary &run, const ptrdiff_t n, const size_t surface_count=1) :
                m_summary_by_lane_read(run.size(), vector2d_t(run.lane_count()*std::max(static_cast<size_t>(1),surface_count))),
                m_lane_count(run.lane_count()),
                m_surface_count(std::max(static_cast<size_t>(1),surface_count))
        {
            for (read_iterator beg = m_summary_by_lane_read.begin(), end = m_summary_by_lane_read.end();
                 beg != end; ++beg)
                reserve(beg->begin(), beg->end(), n);
        }

    public:
        vector_t &operator()(const size_t read, const size_t lane, const size_t surface=0)
        {
            INTEROP_ASSERTMSG(surface < m_surface_count, surface << " < " << m_surface_count);
            INTEROP_ASSERTMSG(read < m_summary_by_lane_read.size(), read << " < " << m_summary_by_lane_read.size());
            const size_t offset = lane*m_surface_count+surface;
            INTEROP_ASSERTMSG(offset < m_summary_by_lane_read[read].size(),
                              offset << " < " << m_summary_by_lane_read[read].size());
            return m_summary_by_lane_read[read][offset];
        }

        void clear()
        {
            for (read_iterator beg = m_summary_by_lane_read.begin(), end = m_summary_by_lane_read.end();
                 beg != end; ++beg)
                clear(beg->begin(), beg->end());
        }

        size_t size() const
        {
            return m_summary_by_lane_read.size();
        }

        size_t read_count() const
        {
            return m_summary_by_lane_read.size();
        }

        size_t lane_count() const
        {
            return m_lane_count;
        }
        size_t surface_count() const
        {
            return m_surface_count;
        }

    private:
        template<typename I>
        static void clear(I beg, I end)
        {
            for (; beg != end; ++beg) beg->clear();
        }

    private:
        vector3d_t m_summary_by_lane_read;
        size_t m_lane_count;
        size_t m_surface_count;
    };

    template<typename I, typename S>
    void summarize(I beg, I end, S &stat, const bool skip_median)
    {
        stat.clear();
        if (beg == end) return;
        stat.mean(util::mean<float>(beg, end));
        stat.stddev(std::sqrt(util::variance_with_mean<float>(beg, end, stat.mean())));
        if(!skip_median) stat.median(util::median_interpolated<float>(beg, end));
    }

    template<typename I, typename S, typename BinaryOp, typename Compare>
    void summarize(I beg, I end, S &stat, BinaryOp op, Compare comp, const bool skip_median)
    {
        if (beg == end) return;
        stat.mean(util::mean<float>(beg, end, op));
        stat.stddev(std::sqrt(util::variance_with_mean<float>(beg, end, stat.mean(), op)));
        if(!skip_median) stat.median(util::median_interpolated<float>(beg, end, comp, op));
    }

    template<typename I, typename S, typename BinaryOp, typename Compare>
    size_t nan_summarize(I beg, I end, S &stat, BinaryOp op, Compare comp, const bool skip_median)
    {
        stat.clear();
        if (beg == end) return 0;
        end = util::remove_nan(beg, end, op);
        if (beg == end) return 0;
        stat.mean(util::mean<float>(beg, end, op));
        INTEROP_ASSERT(!std::isnan(stat.mean()));
        stat.stddev(std::sqrt(util::variance_with_mean<float>(beg, end, stat.mean(), op)));
        if(!skip_median) stat.median(util::median_interpolated<float>(beg, end, comp, op));
        return size_t(std::distance(beg, end));
    }

    template<typename I, typename S, typename Op>
    S nan_accumulate(I beg, I end, const S init, Op op)
    {
        S return_value = init;
        if (beg == end) return init;
        const float temp_value = 0;
        for (; beg != end; ++beg)
        {
            const float current_value = op(temp_value, *beg);
            if(!std::isnan(current_value))
            {
                return_value += static_cast<S>(current_value);
            }
        }
        return return_value;
    }

    inline float divide(const float num, const float div)
    {
        static const float eps = 1e-9f;
        return (div < eps) ? 0 : num / div;
    }
    inline double divide(const double num, const double div)
    {
        static const double eps = 1e-9;
        return (div < eps) ? 0 : num / div;
    }

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

            template<class T>
            size_t operator()(const size_t last, const T &obj) const
            {
                return last + m_op(obj).total_cycles();
            }

        private:
            UnaryOp m_op;
        };

        struct default_get_read
        {
            model::run::read_info &operator()(model::summary::read_summary &summary) const
            {
                return summary.read();
            }

            model::run::read_info &operator()(model::run::read_info &read) const
            {
                return read;
            }

            const model::run::read_info &operator()(const model::summary::read_summary &summary) const
            {
                return summary.read();
            }

            const model::run::read_info &operator()(const model::run::read_info &read) const
            {
                return read;
            }
        };
    }

}}}}