average_vectors.hpp

Go to the documentation of this file.

#ifndef SCRAN_BLOCKS_AVERAGE_VECTORS_HPP
#define SCRAN_BLOCKS_AVERAGE_VECTORS_HPP
 
#include <vector>
#include <limits>
#include <algorithm>
#include <cmath>
#include <numeric>
 
namespace scran_blocks {
 
namespace internal {
 
template<bool weighted_, typename Stat_, typename Weight_, typename Output_>
void average_vectors(size_t n, std::vector<Stat_*> in, const Weight_* w, Output_* out, bool skip_nan) {
    if (in.empty()) {
        std::fill_n(out, n, std::numeric_limits<Output_>::quiet_NaN());
        return;
    } else if (in.size() == 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;
    }
 
    std::fill_n(out, n, 0);
    std::vector<Weight_> accumulated(skip_nan ? n : 0);
 
    auto wcopy = w;
    for (auto current : in) {
        auto copy = out;
 
        if constexpr(weighted_) {
            // Don't skip if weight = 0, as we need to still compute the product, e.g.,
            // if the value is Inf, we'd end up with 0 * Inf => NaN that can't be skipped.
            Weight_ weight = *(wcopy++);
 
            // Use the other loop and skip an unnecessary multiplication when the weight is 1.
            if (weight != 1) { 
                if (skip_nan) {
                    for (size_t i = 0; i < n; ++i, ++current, ++copy) {
                        auto x = *current * weight;
                         if (!std::isnan(x)) {
                            *copy += x; 
                            accumulated[i] += weight;
                         }
                    }
                } else {
                    for (size_t i = 0; i < n; ++i, ++current, ++copy) {
                        *copy += *current * weight;
                    }
                }
                continue;
            }
        }
 
        if (skip_nan) {
            for (size_t i = 0; i < n; ++i, ++current, ++copy) {
                auto x = *current;
                if (!std::isnan(x)) {
                    *copy += x; 
                    ++accumulated[i];
                }
            }
        } else {
            for (size_t i = 0; i < n; ++i, ++current, ++copy) {
                *copy += *current;
            }
        }
    }
 
    if (skip_nan) {
        for (size_t i = 0; i < n; ++i, ++out) {
            *out /= accumulated[i];
        }
    } else {
        double denom = 1;
        if constexpr(weighted_) {
            denom /= std::accumulate(w, w + in.size(), 0.0);
        } else {
            denom /= in.size();
        }
        for (size_t i = 0; i < n; ++i, ++out) {
            *out *= denom;
        }
    }
}
 
}
template<typename Stat_, typename Output_>
void average_vectors(size_t n, std::vector<Stat_*> in, Output_* out, bool skip_nan) {
    internal::average_vectors<false>(n, std::move(in), static_cast<int*>(NULL), out, skip_nan);
    return;
}
 
template<typename Output_ = double, typename Stat_>
std::vector<Output_> average_vectors(size_t n, std::vector<Stat_*> in, bool skip_nan) {
    std::vector<Output_> out(n);
    average_vectors(n, std::move(in), out.data(), skip_nan);
    return out;
}
 
template<typename Stat_, typename Weight_, typename Output_>
void average_vectors_weighted(size_t n, std::vector<Stat_*> in, const Weight_* w, Output_* out, bool skip_nan) {
    if (!in.empty()) {
        bool same = true;
        for (size_t i = 1, end = in.size(); i < end; ++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 {
                average_vectors(n, std::move(in), out, skip_nan);
            }
            return;
        }
    }
 
    internal::average_vectors<true>(n, std::move(in), w, out, skip_nan);
    return;
}
 
template<typename Output_ = double, typename Stat_, typename Weight_>
std::vector<Output_> average_vectors_weighted(size_t n, std::vector<Stat_*> in, const Weight_* w, bool skip_nan) {
    std::vector<Output_> out(n);
    average_vectors_weighted(n, std::move(in), w, out.data(), skip_nan);
    return out;
}
 
}
 
#endif