variance.hpp

Go to the documentation of this file.

#ifndef TATAMI_STATS_VARIANCE_HPP
#define TATAMI_STATS_VARIANCE_HPP
 
#include <vector>
#include <cmath>
#include <numeric>
#include <limits>
#include <algorithm>
#include <cstddef>
#include <optional>
 
#include "tatami/tatami.hpp"
#include "sanisizer/sanisizer.hpp"
#include "quickstats/quickstats.hpp"
 
#include "rss.hpp"
#include "skip_nan/rss.hpp"
#include "utils.hpp"
 
namespace tatami_stats {
 
template<typename Output_ = double>
struct VarianceOptions {
    bool skip_nan = false;
 
    int num_threads = 1;
 
    Output_ mean_placeholder = quickstats::nan_if_available_else_zero<Output_>();
 
    Output_ variance_placeholder = quickstats::nan_if_available_else_zero<Output_>();
};
 
template<typename Output_>
struct VarianceBuffers {
    Output_* mean;
 
    Output_* variance;
};
 
template<typename Value_, typename Index_, typename Output_>
void variance(bool row, const tatami::Matrix<Value_, Index_>& mat, VarianceBuffers<Output_>& output, const VarianceOptions<Output_>& opt) {
    const auto dim = (row ? mat.nrow() : mat.ncol());
    const auto otherdim = (row ? mat.ncol() : mat.nrow());
 
    nanable_ifelse<Value_>(
        opt.skip_nan,
 
        [&]() -> void {
            skip_nan::RssBuffers<Output_, Index_> tmp;
            tmp.mean = output.mean;
            tmp.rss = output.variance;
            auto count = tatami::create_container_of_Index_size<std::vector<Index_> >(row ? mat.nrow() : mat.ncol());
            tmp.count = count.data();
 
            skip_nan::RssOptions ropt;
            ropt.num_threads = opt.num_threads;
            ropt.mean_placeholder = opt.mean_placeholder;
            skip_nan::rss(row, mat, tmp, ropt);
 
            for (Index_ i = 0; i < dim; ++i) {
                if (count[i] <= 1) {
                    output.variance[i] = opt.variance_placeholder;
                } else {
                    output.variance[i] /= count[i] - 1;
                }
            }
        },
 
        [&]() -> void {
            RssBuffers<Output_> tmp;
            tmp.mean = output.mean;
            tmp.rss = output.variance;
 
            RssOptions ropt;
            ropt.num_threads = opt.num_threads;
            ropt.mean_placeholder = opt.mean_placeholder;
            rss(row, mat, tmp, ropt);
 
            if (otherdim <= 1) {
                std::fill_n(output.variance, dim, opt.variance_placeholder);
            } else {
                for (Index_ i = 0; i < dim; ++i) {
                    output.variance[i] /= otherdim - 1;
                }
            }
        }
    );
}
 
template<typename Output_>
struct VarianceResult {
    std::vector<Output_> mean;
 
    std::vector<Output_> variance;
};
 
template<typename Output_ = double, typename Value_, typename Index_>
VarianceResult<Output_> variance(bool row, const tatami::Matrix<Value_, Index_>& mat, const VarianceOptions<Output_>& opt) {
    VarianceResult<Output_> output;
    const auto dim = (row ? mat.nrow() : mat.ncol());
    tatami::resize_container_to_Index_size(output.mean, dim
#ifdef TATAMI_STATS_TEST_DIRTY
        , -1
#endif
    );
    tatami::resize_container_to_Index_size(output.variance, dim
#ifdef TATAMI_STATS_TEST_DIRTY
        , -1
#endif
    );
 
    VarianceBuffers<Output_> buffers;
    buffers.mean = output.mean.data();
    buffers.variance = output.variance.data();
 
    variance(row, mat, buffers, opt);
    return output;
}
 
}
 
#endif