find_median_mad.hpp

Go to the documentation of this file.

#ifndef SCRAN_QC_FIND_MEDIAN_MAD_H
#define SCRAN_QC_FIND_MEDIAN_MAD_H
 
#include <vector>
#include <limits>
#include <cmath>
#include <algorithm>
#include <cstddef>
 
#include "tatami_stats/tatami_stats.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "utils.hpp"
 
namespace scran_qc {
 
struct FindMedianMadOptions {
    bool log = false;
 
    bool median_only = false;
};
 
template<typename Float_>
struct FindMedianMadResults {
    FindMedianMadResults(Float_ m1, Float_ m2) : median(m1), mad(m2) {}
    FindMedianMadResults() = default;
    Float_ median = 0;
 
    Float_ mad = 0;
};
 
template<typename Float_> 
FindMedianMadResults<Float_> find_median_mad(std::size_t num, Float_* metrics, const FindMedianMadOptions& options) {
    static_assert(std::is_floating_point<Float_>::value);
 
    // Rotate all the NaNs to the front of the buffer and ignore them.
    decltype(I(num)) lost = 0;
    for (decltype(I(num)) i = 0; i < num; ++i) {
        if (std::isnan(metrics[i])) {
            std::swap(metrics[i], metrics[lost]);
            ++lost;
        }
    }
    metrics += lost;
    num -= lost;
 
    if (options.log) {
        for (decltype(I(num)) i = 0; i < num; ++i) {
            auto& val = metrics[i];
            if (val > 0) {
                val = std::log(val);
            } else if (val == 0) {
                val = -std::numeric_limits<double>::infinity();
            } else {
                throw std::runtime_error("cannot log-transform negative values");
            }
        }
    }
 
    // No need to skip the NaNs again.
    const auto median = tatami_stats::medians::direct<Float_>(metrics, num, /* skip_nan = */ false);
 
    if (options.median_only || std::isnan(median)) {
        // Giving up.
        return FindMedianMadResults<Float_>(median, std::numeric_limits<Float_>::quiet_NaN());
    } else if (std::isinf(median)) {
        // MADs should be no-ops when added/subtracted from infinity. Any
        // finite value will do here, so might as well keep it simple.
        return FindMedianMadResults<Float_>(median, static_cast<Float_>(0));
    }
 
    // As an aside, there's no way to avoid passing in 'metrics' as a Float_,
    // even if the original values were integers, because we need to do this
    // subtraction here that could cast integers to floats. So at some point we
    // will need a floating-point buffer, and so we might as well just pass the
    // metrics in as floats in the first place. Technically the first sort
    // could be done with an integer buffer but then we'd need an extra argument.
 
    for (decltype(I(num)) i = 0; i < num; ++i) {
        metrics[i] = std::abs(metrics[i] - median);
    }
    auto mad = tatami_stats::medians::direct<Float_>(metrics, num, /* skip_nan = */ false);
    mad *= 1.4826; // for equivalence with the standard deviation under normality.
 
    return FindMedianMadResults<Float_>(median, mad);
}
 
template<typename Float_ = double, typename Value_> 
FindMedianMadResults<Float_> find_median_mad(const std::size_t num, const Value_* const metrics, Float_* buffer, const FindMedianMadOptions& options) {
    std::vector<Float_> xbuffer;
    if (buffer == NULL) {
        sanisizer::resize(xbuffer, num
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        buffer = xbuffer.data();
    }
    std::copy_n(metrics, num, buffer);
    return find_median_mad(num, buffer, options);
}
 
template<typename Float_>
class FindMedianMadWorkspace {
public:
    template<typename Block_>
    FindMedianMadWorkspace(const std::size_t num, const Block_* const block) :
        my_buffer(sanisizer::cast<decltype(I(my_buffer.size()))>(num))
    {
        set(num, block);
    }
 
    FindMedianMadWorkspace() = default;
 
    template<typename Block_>
    void set(const std::size_t num, const Block_* const block) {
        my_block_starts.clear();
 
        if (block) { 
            for (decltype(I(num)) i = 0; i < num; ++i) {
                const auto candidate = block[i];
                if (sanisizer::is_greater_than_or_equal(candidate, my_block_starts.size())) {
                    my_block_starts.resize(sanisizer::sum<decltype(I(my_block_starts.size()))>(candidate, 1));
                }
                ++my_block_starts[candidate];
            }
 
            std::size_t sofar = 0;
            for (auto& s : my_block_starts) {
                const auto last = sofar;
                sofar += s;
                s = last;
            }
        }
 
        sanisizer::resize(my_buffer, num
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
        sanisizer::resize(my_block_ends, num
#ifdef SCRAN_QC_TEST_INIT
            , SCRAN_QC_TEST_INIT
#endif
        );
    }
 
public:
    // Can't figure out how to make compute_blocked() a friend,
    // so these puppies are public for simplicity.
    std::vector<Float_> my_buffer;
    std::vector<std::size_t> my_block_starts;
    std::vector<std::size_t> my_block_ends;
};
 
template<typename Output_ = double, typename Value_, typename Block_>
std::vector<FindMedianMadResults<Output_> > find_median_mad_blocked(
    const std::size_t num,
    const Value_* const metrics, 
    const Block_* const block,
    FindMedianMadWorkspace<Output_>* workspace,
    const FindMedianMadOptions& options)
{
    std::unique_ptr<FindMedianMadWorkspace<Output_> > xworkspace;
    if (workspace == NULL) {
        xworkspace = std::make_unique<FindMedianMadWorkspace<Output_> >(num, block);
        workspace = xworkspace.get();
    }
 
    std::vector<FindMedianMadResults<Output_> > output;
 
    auto& buffer = workspace->my_buffer;
    if (!block) {
        std::copy_n(metrics, num, buffer.begin());
        output.push_back(find_median_mad(num, buffer.data(), options));
        return output;
    }
 
    const auto& starts = workspace->my_block_starts;
    auto& ends = workspace->my_block_ends;
    std::copy(starts.begin(), starts.end(), ends.begin());
    for (decltype(I(num)) i = 0; i < num; ++i) {
        auto& pos = ends[block[i]];
        buffer[pos] = metrics[i];
        ++pos;
    }
 
    // Using the ranges on the buffer.
    const auto nblocks = starts.size();
    output.reserve(nblocks);
    for (decltype(I(nblocks)) g = 0; g < nblocks; ++g) {
        output.push_back(find_median_mad(ends[g] - starts[g], buffer.data() + starts[g], options));
    }
 
    return output;
}
 
}
 
#endif