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 <cstdint>
 
#include "tatami_stats/tatami_stats.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 Index_, typename Float_> 
FindMedianMadResults<Float_> find_median_mad(Index_ 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.
    Index_ lost = 0;
    for (Index_ i = 0; i < num; ++i) {
        if (std::isnan(metrics[i])) {
            std::swap(metrics[i], metrics[lost]);
            ++lost;
        }
    }
    metrics += lost;
    num -= lost;
 
    if (options.log) {
        auto copy = metrics;
        for (Index_ i = 0; i < num; ++i, ++copy) {
            if (*copy > 0) {
                *copy = std::log(*copy);
            } else if (*copy == 0) {
                *copy = -std::numeric_limits<double>::infinity();
            } else {
                throw std::runtime_error("cannot log-transform negative values");
            }
        }
    }
 
    // No need to skip the NaNs again.
    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.
 
    auto copy = metrics;
    for (Index_ i = 0; i < num; ++i, ++copy) {
        *copy = std::abs(*copy - 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 Index_, typename Value_> 
FindMedianMadResults<Float_> find_median_mad(Index_ num, const Value_* metrics, Float_* buffer, const FindMedianMadOptions& options) {
    std::unique_ptr<std::vector<Float_> > xbuffer;
    if (buffer == NULL) {
        xbuffer = std::make_unique<std::vector<Float_> >(num);
        buffer = xbuffer->data();
    }
    std::copy_n(metrics, num, buffer);
    return find_median_mad(num, buffer, options);
}
 
template<typename Float_, typename Index_>
class FindMedianMadWorkspace {
public:
    template<typename Block_>
    FindMedianMadWorkspace(Index_ num, const Block_* block) : my_buffer(num) {
        set(num, block);
    }
 
    FindMedianMadWorkspace() = default;
 
    template<typename Block_>
    void set(Index_ num, const Block_* block) {
        my_block_starts.clear();
 
        if (block) { 
            for (Index_ i = 0; i < num; ++i) {
                size_t candidate = block[i];
                if (candidate >= my_block_starts.size()) {
                    my_block_starts.resize(candidate + 1);
                }
                ++my_block_starts[candidate];
            }
 
            Index_ sofar = 0;
            for (auto& s : my_block_starts) {
                Index_ last = sofar;
                sofar += s;
                s = last;
            }
        }
 
        my_buffer.resize(num);
        my_block_ends.resize(my_block_starts.size());
    }
 
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<Index_> my_block_starts;
    std::vector<Index_> my_block_ends;
};
 
template<typename Output_ = double, typename Index_, typename Value_, typename Block_>
std::vector<FindMedianMadResults<Output_> > find_median_mad_blocked(
    Index_ num,
    const Value_* metrics, 
    const Block_* block,
    FindMedianMadWorkspace<Output_, Index_>* workspace,
    const FindMedianMadOptions& options)
{
    std::unique_ptr<FindMedianMadWorkspace<Output_, Index_> > xworkspace;
    if (workspace == NULL) {
        xworkspace = std::make_unique<FindMedianMadWorkspace<Output_, Index_> >(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 (Index_ i = 0; i < num; ++i) {
        auto& pos = ends[block[i]];
        buffer[pos] = metrics[i];
        ++pos;
    }
 
    // Using the ranges on the buffer.
    size_t nblocks = starts.size();
    output.reserve(nblocks);
    for (size_t g = 0; g < nblocks; ++g) {
        output.push_back(find_median_mad(ends[g] - starts[g], buffer.data() + starts[g], options));
    }
 
    return output;
}
 
}
 
#endif