ComputeMedianMad.hpp

Go to the documentation of this file.

#ifndef SCRAN_COMPUTE_MEDIAN_MAD_H
#define SCRAN_COMPUTE_MEDIAN_MAD_H
 
#include "../utils/macros.hpp"
 
#include <vector>
#include <limits>
#include <cmath>
#include <algorithm>
#include <cstdint>
 
#include "utils.hpp"
 
namespace scran {
 
class ComputeMedianMad {
public:
    struct Defaults {
        static constexpr bool log = false;
 
        static constexpr bool median_only = false;
    };
 
private:
    bool log = Defaults::log;
 
    bool median_only = Defaults::median_only;
 
public:
    ComputeMedianMad& set_log(bool l = Defaults::log) {
        log = l;
        return *this;
    }
 
    ComputeMedianMad& set_median_only(bool m = Defaults::median_only) {
        median_only = m;
        return *this;
    }
 
public:
    struct Results {
        Results(int nblocks=0) : medians(nblocks), mads(nblocks) {}
        std::vector<double> medians;
 
        std::vector<double> mads;
 
        bool log = false;
    };
 
public:
    template<typename Block>
    static std::vector<int> compute_block_starts(size_t n, const Block* block) {
        std::vector<int> starts;
 
        for (size_t i = 0; i < n; ++i) {
            size_t candidate = block[i] + 1;
            if (candidate > starts.size()) {
                starts.resize(candidate);
            }
            ++starts[block[i]];
        }
        
        int sofar = 0;
        for (auto& s : starts) {
            int last = sofar;
            sofar += s;
            s = last;
        }
 
        return starts;
    }
 
public:
    template<typename Input, typename Buffer> 
    Results run(size_t n, const Input* metrics, Buffer* buffer) const {
        Results output(1);
        output.log = log;
 
        auto copy = buffer;
        for (size_t i = 0; i < n; ++i) {
            if (!std::isnan(metrics[i])) {
                *copy = metrics[i];
                ++copy;
            }
        }
 
        compute(copy - buffer, buffer, output.medians[0], output.mads[0]);
        return output;
    }
 
    template<typename Input>
    Results run(size_t n, const Input* metrics) const {
        std::vector<double> buffer(n);
        return run(n, metrics, buffer.data());
    } 
 
public:
    template<typename Input, typename Block, typename Buffer>
    Results run_blocked(size_t n, const Block* block, std::vector<int> starts, const Input* metrics, Buffer* buffer) const {
        if (block) {
            // Unscrambling into the buffer.
            auto ends = starts;
            for (size_t i = 0; i < n; ++i) {
                if (!std::isnan(metrics[i])) {
                    auto& pos = ends[block[i]];
                    buffer[pos] = metrics[i];
                    ++pos;
                }
            }
 
            // Using the ranges on the buffer.
            size_t nblocks = starts.size();
            Results output(nblocks);
            output.log = log;
            for (size_t g = 0; g < nblocks; ++g) {
                compute(ends[g] - starts[g], buffer + starts[g], output.medians[g], output.mads[g]);
            }
 
            return output;
        } else {
            return run(n, metrics, buffer);
        }
    }
 
    template<typename Block, typename Input>
    Results run_blocked(size_t n, const Block* block, const Input* metrics) const {
        std::vector<double> buffer(n);
        return run_blocked(n, block, compute_block_starts(n, block), metrics, buffer.data());
    }
 
private:
    template<typename Buffer>
    void compute(size_t n, Buffer* ptr, double& median, double& mad) const {
        if (!n) {
            median = std::numeric_limits<double>::quiet_NaN();
            mad = std::numeric_limits<double>::quiet_NaN();
            return;
        }
 
        if (log) {
            auto copy = ptr;
            for (size_t i = 0; i < n; ++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");
                }
            }
        }
 
        // First, getting the median.
        size_t halfway = n / 2;
        median = compute_median(ptr, halfway, n);
 
        if (median_only || std::isnan(median)) {
            // Giving up.
            mad = std::numeric_limits<double>::quiet_NaN();
            return;
        } 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.
            mad = 0;
            return;
        }
 
        // Now getting the MAD. 
        auto copy = ptr;
        for (size_t i = 0; i < n; ++i, ++copy) {
            *copy = std::abs(*copy - median);
        }
        mad = compute_median(ptr, halfway, n);
        mad *= 1.4826; // for equivalence with the standard deviation under normality.
        return;
    }
 
    static double compute_median (double* ptr, size_t halfway, size_t n) {
        std::nth_element(ptr, ptr + halfway, ptr + n);
        double medtmp = *(ptr + halfway);
 
        if (n % 2 == 0) {
            std::nth_element(ptr, ptr + halfway - 1, ptr + n);
            double left = *(ptr + halfway - 1);
 
            if (std::isinf(left) && std::isinf(medtmp)) {
                if ((left > 0) != (medtmp > 0)) {
                    return std::numeric_limits<double>::quiet_NaN();
                } else {
                    return medtmp;
                }
            } 
 
            medtmp = (medtmp + left)/2;
        }
 
        return medtmp;
    }
};
 
}
 
#endif