scran_blocks
Blocking utilities for libscran
Loading...
Searching...
No Matches
average_vectors.hpp
Go to the documentation of this file.
1#ifndef SCRAN_BLOCKS_AVERAGE_VECTORS_HPP
2#define SCRAN_BLOCKS_AVERAGE_VECTORS_HPP
3
4#include <vector>
5#include <limits>
6#include <algorithm>
7#include <cmath>
8#include <numeric>
9#include <cstddef>
10
11#include "sanisizer/sanisizer.hpp"
12
13#include "utils.hpp"
14
20namespace scran_blocks {
21
25namespace internal {
26
27template<bool weighted_, typename Stat_, typename Weight_, typename Output_>
28void average_vectors(const std::size_t n, std::vector<Stat_*> in, const Weight_* const w, Output_* const out, const bool skip_nan) {
29 if (in.empty()) {
30 std::fill_n(out, n, std::numeric_limits<Output_>::quiet_NaN());
31 return;
32 } else if (in.size() == 1) {
33 if constexpr(weighted_) {
34 if (w[0] == 0) {
35 std::fill_n(out, n, std::numeric_limits<Output_>::quiet_NaN());
36 return;
37 }
38 }
39 std::copy(in[0], in[0] + n, out);
40 return;
41 }
42
43 std::fill_n(out, n, 0);
44 std::vector<Weight_> accumulated(skip_nan ? n : 0);
45
46 auto wcopy = w;
47 for (const auto current : in) {
48 if constexpr(weighted_) {
49 // Don't skip if weight = 0, as we need to still compute the product, e.g.,
50 // if the value is Inf, we'd end up with 0 * Inf => NaN that can't be skipped.
51 const Weight_ weight = *(wcopy++);
52
53 // Use the other loop and skip an unnecessary multiplication when the weight is 1.
54 if (weight != 1) {
55 if (skip_nan) {
56 for (decltype(I(n)) i = 0; i < n; ++i) {
57 const auto x = current[i] * weight;
58 if (!std::isnan(x)) {
59 out[i] += x;
60 accumulated[i] += weight;
61 }
62 }
63 } else {
64 for (decltype(I(n)) i = 0; i < n; ++i) {
65 out[i] += current[i] * weight;
66 }
67 }
68 continue;
69 }
70 }
71
72 if (skip_nan) {
73 for (decltype(I(n)) i = 0; i < n; ++i) {
74 const auto x = current[i];
75 if (!std::isnan(x)) {
76 out[i] += x;
77 ++accumulated[i];
78 }
79 }
80 } else {
81 for (decltype(I(n)) i = 0; i < n; ++i) {
82 out[i] += current[i];
83 }
84 }
85 }
86
87 if (skip_nan) {
88 for (decltype(I(n)) i = 0; i < n; ++i) {
89 out[i] /= accumulated[i];
90 }
91 } else {
92 Output_ denom = 1;
93 if constexpr(weighted_) {
94 denom /= std::accumulate(w, w + in.size(), static_cast<Output_>(0));
95 } else {
96 denom /= in.size();
97 }
98 for (decltype(I(n)) i = 0; i < n; ++i) {
99 out[i] *= denom;
100 }
101 }
102}
103
104}
124template<typename Stat_, typename Output_>
125void average_vectors(const std::size_t n, std::vector<Stat_*> in, Output_* const out, const bool skip_nan) {
126 internal::average_vectors<false>(n, std::move(in), static_cast<int*>(NULL), out, skip_nan);
127 return;
128}
129
143template<typename Output_ = double, typename Stat_>
144std::vector<Output_> average_vectors(const std::size_t n, std::vector<Stat_*> in, const bool skip_nan) {
145 auto out = sanisizer::create<std::vector<Output_> >(n);
146 average_vectors(n, std::move(in), out.data(), skip_nan);
147 return out;
148}
149
168template<typename Stat_, typename Weight_, typename Output_>
169void average_vectors_weighted(const std::size_t n, std::vector<Stat_*> in, const Weight_* const w, Output_* const out, const bool skip_nan) {
170 if (!in.empty()) {
171 bool same = true;
172 const auto numin = in.size();
173 for (decltype(I(numin)) i = 1; i < numin; ++i) {
174 if (w[i] != w[0]) {
175 same = false;
176 break;
177 }
178 }
179
180 if (same) {
181 if (w[0] == 0) {
182 std::fill_n(out, n, std::numeric_limits<Output_>::quiet_NaN());
183 } else {
184 average_vectors(n, std::move(in), out, skip_nan);
185 }
186 return;
187 }
188 }
189
190 internal::average_vectors<true>(n, std::move(in), w, out, skip_nan);
191 return;
192}
193
209template<typename Output_ = double, typename Stat_, typename Weight_>
210std::vector<Output_> average_vectors_weighted(const std::size_t n, std::vector<Stat_*> in, const Weight_* const w, const bool skip_nan) {
211 auto out = sanisizer::create<std::vector<Output_> >(n);
212 average_vectors_weighted(n, std::move(in), w, out.data(), skip_nan);
213 return out;
214}
215
216}
217
218#endif
Blocking utilities for libscran.
Definition average_vectors.hpp:20
void average_vectors(const std::size_t n, std::vector< Stat_ * > in, Output_ *const out, const bool skip_nan)
Definition average_vectors.hpp:125
void average_vectors_weighted(const std::size_t n, std::vector< Stat_ * > in, const Weight_ *const w, Output_ *const out, const bool skip_nan)
Definition average_vectors.hpp:169