parallel_means.hpp

Go to the documentation of this file.

#ifndef SCRAN_BLOCKS_PARALLEL_MEANS_HPP
#define SCRAN_BLOCKS_PARALLEL_MEANS_HPP
 
#include <vector>
#include <limits>
#include <algorithm>
#include <cmath>
#include <numeric>
#include <cstddef>
 
#include "sanisizer/sanisizer.hpp"
 
#include "utils.hpp"
 
namespace scran_blocks {
 
template<bool weighted_, typename Stat_, typename Weight_, typename Output_>
void parallel_means_internal(const std::size_t n, std::vector<Stat_*> in, const Weight_* const w, Output_* const out, const bool skip_nan) {
    const auto nblocks = in.size();
    if (nblocks == 0) {
        std::fill_n(out, n, std::numeric_limits<Output_>::quiet_NaN());
        return;
    } else if (nblocks == 1) {
        if constexpr(weighted_) {
            if (w[0] == 0) {
                std::fill_n(out, n, std::numeric_limits<Output_>::quiet_NaN());
                return;
            }
        } 
        std::copy(in[0], in[0] + n, out);
        return;
    }
 
    if (skip_nan) {
        for (I<decltype(n)> i = 0; i < n; ++i) {
            Output_ prod = 0;
            Output_ denom = 0;
 
            for (I<decltype(nblocks)> b = 0; b < nblocks; ++b) {
                const auto val = in[b][i];
                if (!std::isnan(val)) {
                    if constexpr(weighted_) {
                        const auto curw = w[b];
                        prod += val * curw;
                        denom += curw;
                    } else {
                        prod += val;
                        denom += 1;
                    }
                }
            }
 
            out[i] = prod/denom;
        }
 
    } else {
        std::fill_n(out, n, 0);
        auto wcopy = w;
        for (const auto current : in) {
            if constexpr(weighted_) {
                const Weight_ weight = *(wcopy++);
                if (weight != 1) { 
                    for (I<decltype(n)> i = 0; i < n; ++i) {
                        out[i] += current[i] * weight;
                    }
                    continue;
                }
            }
            for (I<decltype(n)> i = 0; i < n; ++i) {
                out[i] += current[i];
            }
        }
 
        const Output_ denom = [&]() {
            if constexpr(weighted_) {
                return std::accumulate(w, w + in.size(), static_cast<Output_>(0));
            } else {
                return in.size();
            }
        }();
        for (I<decltype(n)> i = 0; i < n; ++i) {
            out[i] /= denom;
        }
    }
}
template<typename Stat_, typename Output_>
void parallel_means(const std::size_t n, std::vector<Stat_*> in, Output_* const out, const bool skip_nan) {
    parallel_means_internal<false>(n, std::move(in), static_cast<int*>(NULL), out, skip_nan);
    return;
}
 
template<typename Output_ = double, typename Stat_>
std::vector<Output_> parallel_means(const std::size_t n, std::vector<Stat_*> in, const bool skip_nan) {
    auto out = sanisizer::create<std::vector<Output_> >(n);
    parallel_means(n, std::move(in), out.data(), skip_nan);
    return out;
}
 
template<typename Stat_, typename Weight_, typename Output_>
void parallel_weighted_means(const std::size_t n, std::vector<Stat_*> in, const Weight_* const w, Output_* const out, const bool skip_nan) {
    if (!in.empty()) {
        bool same = true;
        const auto numin = in.size();
        for (I<decltype(numin)> i = 1; i < numin; ++i) {
            if (w[i] != w[0]) {
                same = false;
                break;
            }
        }
 
        if (same) {
            if (w[0] == 0) {
                std::fill_n(out, n, std::numeric_limits<Output_>::quiet_NaN());
            } else {
                parallel_means(n, std::move(in), out, skip_nan);
            }
            return;
        }
    }
 
    parallel_means_internal<true>(n, std::move(in), w, out, skip_nan);
    return;
}
 
template<typename Output_ = double, typename Stat_, typename Weight_>
std::vector<Output_> parallel_weighted_means(const std::size_t n, std::vector<Stat_*> in, const Weight_* const w, const bool skip_nan) {
    auto out = sanisizer::create<std::vector<Output_> >(n);
    parallel_weighted_means(n, std::move(in), w, out.data(), skip_nan);
    return out;
}
 
// Methods for back-compatibility.
template<typename Stat_, typename Output_>
void average_vectors(const std::size_t n, std::vector<Stat_*> in, Output_* const out, const bool skip_nan) {
    parallel_means(n, in, out, skip_nan);
}
 
template<typename Output_ = double, typename Stat_>
std::vector<Output_> average_vectors(const std::size_t n, std::vector<Stat_*> in, const bool skip_nan) {
    return parallel_means(n, in, skip_nan);
}
 
template<typename Stat_, typename Weight_, typename Output_>
void average_vectors_weighted(const std::size_t n, std::vector<Stat_*> in, const Weight_* const w, Output_* const out, const bool skip_nan) {
    parallel_weighted_means(n, in, w, out, skip_nan);
}
 
template<typename Output_ = double, typename Stat_, typename Weight_>
std::vector<Output_> average_vectors_weighted(const std::size_t n, std::vector<Stat_*> in, const Weight_* const w, const bool skip_nan) {
    return parallel_weighted_means(n, in, w, skip_nan);
}
}
 
#endif