1#ifndef SCRAN_MARKERS_SUMMARIZE_COMPARISONS_HPP
2#define SCRAN_MARKERS_SUMMARIZE_COMPARISONS_HPP
10#include "tatami_stats/tatami_stats.hpp"
11#include "sanisizer/sanisizer.hpp"
28template<
typename Stat_ =
double,
typename Rank_ =
int>
72template<
typename Stat_ =
double,
typename Rank_ =
int>
78 std::vector<Stat_>
min;
96 std::vector<Stat_>
max;
110template<
typename Stat_,
typename Gene_,
typename Rank_>
111void summarize_comparisons(
112 const std::size_t ngroups,
113 const Stat_*
const effects,
114 const std::size_t group,
117 std::vector<Stat_>& buffer)
120 std::size_t ncomps = 0;
121 for (
decltype(I(ngroups)) r = 0; r < ngroups; ++r) {
122 if (r == group || std::isnan(effects[r])) {
125 buffer[ncomps] = effects[r];
130 Stat_ val = (ncomps == 0 ? std::numeric_limits<Stat_>::quiet_NaN() : buffer[0]);
132 output.
min[gene] = val;
135 output.
mean[gene] = val;
138 output.
max[gene] = val;
141 output.
median[gene] = val;
145 const auto ebegin = buffer.data(), elast = ebegin + ncomps;
147 output.
min[gene] = *std::min_element(ebegin, elast);
150 output.
mean[gene] = std::accumulate(ebegin, elast,
static_cast<Stat_
>(0)) / ncomps;
153 output.
max[gene] = *std::max_element(ebegin, elast);
156 output.
median[gene] = tatami_stats::medians::direct(ebegin, ncomps,
false);
161template<
typename Gene_,
typename Stat_,
typename Rank_>
162void summarize_comparisons(
164 const std::size_t ngroups,
165 const Stat_*
const effects,
166 const std::vector<SummaryBuffers<Stat_, Rank_> >& output,
170 auto buffer = sanisizer::create<std::vector<Stat_> >(ngroups);
171 for (Gene_ gene = start, end = start + length; gene < end; ++gene) {
172 for (
decltype(I(ngroups)) l = 0; l < ngroups; ++l) {
173 const auto current_effects = effects + sanisizer::nd_offset<std::size_t>(0, ngroups, l, ngroups, gene);
174 summarize_comparisons(ngroups, current_effects, l, gene, output[l], buffer);
180template<
typename Stat_,
typename Gene_>
181Gene_ fill_and_sort_rank_buffer(
const Stat_*
const effects,
const std::size_t stride, std::vector<std::pair<Stat_, Gene_> >& buffer) {
183 for (Gene_ i = 0, end = buffer.size(); i < end; ++i) {
184 const auto cureffect = effects[sanisizer::product_unsafe<std::size_t>(i, stride)];
185 if (!std::isnan(cureffect)) {
186 auto& current = buffer[counter];
187 current.first = cureffect;
195 buffer.begin() + counter,
196 [&](
const std::pair<Stat_, Gene_>& left,
const std::pair<Stat_, Gene_>& right) ->
bool {
198 if (left.first == right.first) {
199 return left.second < right.second;
201 return left.first > right.first;
209template<
typename Stat_,
typename Gene_,
typename Rank_>
210void compute_min_rank_internal(
const Gene_ use,
const std::vector<std::pair<Stat_, Gene_> >& buffer, Rank_*
const output) {
212 for (Gene_ i = 0; i < use; ++i) {
213 auto& current = output[buffer[i].second];
214 if (counter < current) {
221template<
typename Stat_,
typename Gene_,
typename Rank_>
222void compute_min_rank_for_group(
const Gene_ ngenes,
const std::size_t ngroups,
const std::size_t group,
const Stat_*
const effects, Rank_*
const output,
const int threads) {
223 std::vector<std::vector<Rank_> > stores(threads - 1);
224 const auto maxrank_placeholder = ngenes;
225 std::fill_n(output, ngenes, maxrank_placeholder);
227 tatami::parallelize([&](
const int t,
const std::size_t start,
const std::size_t length) ->
void {
232 auto& curstore = stores[t - 1];
233 if (curstore.empty()) {
234 sanisizer::resize(curstore, ngenes, maxrank_placeholder);
236 curoutput = curstore.data();
239 auto buffer = sanisizer::create<std::vector<std::pair<Stat_, Gene_> > >(ngenes);
240 for (
decltype(I(start)) g = start, end = start + length; g < end; ++g) {
244 const auto used = fill_and_sort_rank_buffer(effects + g, ngroups, buffer);
245 compute_min_rank_internal(used, buffer, curoutput);
247 }, ngroups, threads);
249 for (
const auto& curstore : stores) {
251 for (
auto x : curstore) {
260template<
typename Stat_,
typename Gene_,
typename Rank_>
261void compute_min_rank_pairwise(
263 const std::size_t ngroups,
264 const Stat_*
const effects,
265 const std::vector<SummaryBuffers<Stat_, Rank_> >& output,
268 const auto ngroups2 = sanisizer::product_unsafe<std::size_t>(ngroups, ngroups);
270 tatami::parallelize([&](
const int,
const std::size_t start,
const std::size_t length) ->
void {
271 auto buffer = sanisizer::create<std::vector<std::pair<Stat_, Gene_> > >(ngenes);
272 for (
decltype(I(start)) g = start, end = start + length; g < end; ++g) {
273 const auto target = output[g].min_rank;
274 if (target == NULL) {
278 std::fill_n(target, ngenes, ngenes);
280 for (
decltype(I(ngroups)) g2 = 0; g2 < ngroups; ++g2) {
284 const auto offset = sanisizer::nd_offset<std::size_t>(g2, ngroups, g);
285 const auto used = fill_and_sort_rank_buffer(effects + offset, ngroups2, buffer);
286 compute_min_rank_internal(used, buffer, target);
289 }, ngroups, threads);
292template<
typename Gene_,
typename Stat_,
typename Rank_>
293SummaryBuffers<Stat_, Rank_> fill_summary_results(
295 SummaryResults<Stat_, Rank_>& out,
296 const bool compute_min,
297 const bool compute_mean,
298 const bool compute_median,
299 const bool compute_max,
300 const bool compute_min_rank)
302 SummaryBuffers<Stat_, Rank_> ptr;
303 const auto out_len = sanisizer::cast<typename std::vector<Stat_>::size_type>(ngenes);
306 out.min.resize(out_len
307#ifdef SCRAN_MARKERS_TEST_INIT
308 , SCRAN_MARKERS_TEST_INIT
311 ptr.min = out.min.data();
314 out.mean.resize(out_len
315#ifdef SCRAN_MARKERS_TEST_INIT
316 , SCRAN_MARKERS_TEST_INIT
319 ptr.mean = out.mean.data();
321 if (compute_median) {
322 out.median.resize(out_len
323#ifdef SCRAN_MARKERS_TEST_INIT
324 , SCRAN_MARKERS_TEST_INIT
327 ptr.median = out.median.data();
330 out.max.resize(out_len
331#ifdef SCRAN_MARKERS_TEST_INIT
332 , SCRAN_MARKERS_TEST_INIT
335 ptr.max = out.max.data();
337 if (compute_min_rank) {
338 out.min_rank.resize(out_len
339#ifdef SCRAN_MARKERS_TEST_INIT
340 , SCRAN_MARKERS_TEST_INIT
343 ptr.min_rank = out.min_rank.data();
349template<
typename Gene_,
typename Stat_,
typename Rank_>
350std::vector<SummaryBuffers<Stat_, Rank_> > fill_summary_results(
352 const std::size_t ngroups,
353 std::vector<SummaryResults<Stat_, Rank_> >& outputs,
354 const bool compute_min,
355 const bool compute_mean,
356 const bool compute_median,
357 const bool compute_max,
358 const bool compute_min_rank)
360 sanisizer::resize(outputs, ngroups);
361 std::vector<SummaryBuffers<Stat_, Rank_> > ptrs;
362 ptrs.reserve(ngroups);
363 for (
decltype(I(ngroups)) g = 0; g < ngroups; ++g) {
364 ptrs.emplace_back(fill_summary_results(ngenes, outputs[g], compute_min, compute_mean, compute_median, compute_max, compute_min_rank));
Marker detection for single-cell data.
Definition score_markers_pairwise.hpp:25
void parallelize(Function_ fun, Index_ tasks, int threads)
Pointers to arrays to hold the summary statistics.
Definition summarize_comparisons.hpp:29
Stat_ * mean
Definition summarize_comparisons.hpp:42
Stat_ * median
Definition summarize_comparisons.hpp:49
Stat_ * min
Definition summarize_comparisons.hpp:35
Stat_ * max
Definition summarize_comparisons.hpp:56
Rank_ * min_rank
Definition summarize_comparisons.hpp:63
Container for the summary statistics.
Definition summarize_comparisons.hpp:73
std::vector< Stat_ > mean
Definition summarize_comparisons.hpp:84
std::vector< Stat_ > median
Definition summarize_comparisons.hpp:90
std::vector< Stat_ > min
Definition summarize_comparisons.hpp:78
std::vector< Rank_ > min_rank
Definition summarize_comparisons.hpp:102
std::vector< Stat_ > max
Definition summarize_comparisons.hpp:96