range.hpp

Go to the documentation of this file.

#ifndef TATAMI_STATS_SKIP_NAN_RANGE_HPP
#define TATAMI_STATS_SKIP_NAN_RANGE_HPP
 
#include "../utils.hpp"
 
#include <vector>
#include <algorithm>
#include <type_traits>
#include <limits>
 
#include "tatami/tatami.hpp"
 
namespace tatami_stats {
 
namespace skip_nan { 
 
template<typename Output_>
constexpr Output_ default_minimum_placeholder() {
    if constexpr(std::numeric_limits<Output_>::has_infinity) {
        return std::numeric_limits<Output_>::infinity();
    } else {
        return std::numeric_limits<Output_>::max();
    }
}
 
template<typename Output_>
constexpr Output_ default_maximum_placeholder() {
    if constexpr(std::numeric_limits<Output_>::has_infinity) {
        return -std::numeric_limits<Output_>::infinity();
    } else {
        return std::numeric_limits<Output_>::lowest();
    }
}
 
template<typename Output_ = double>
struct RangeOptions {
    int num_threads = 1;
 
    Output_ minimum_placeholder = default_minimum_placeholder<Output_>();
 
    Output_ maximum_placeholder = default_maximum_placeholder<Output_>();
};
 
template<typename Output_, typename Index_>
struct RangeDirectResult {
    Output_ minimum;
    Output_ maximum;
    Index_ count;
};
 
template<typename Value_, typename Index_, typename Output_>
RangeDirectResult<Output_, Index_> range_direct(const Value_* const ptr, const Index_ num, const RangeOptions<Output_>& opt) {
    RangeDirectResult<Output_, Index_> output;
    output.minimum = opt.minimum_placeholder;
    output.maximum = opt.maximum_placeholder;
    output.count = 0;
 
    // First loop to get to the first non-NA value.
    Index_ i = 0;
    for (; i < num; ++i) {
        auto val = ptr[i];
        if (!std::isnan(val)) {
            output.minimum = val;
            output.maximum = val;
            output.count = 1;
            ++i;
            break;
        }
    }
 
    // Second loop to actually get the minimum.
    for (; i < num; ++i) {
        auto val = ptr[i];
        if (!std::isnan(val)) {
            ++output.count;
            output.minimum = std::min(output.minimum, val);
            output.maximum = std::max(output.maximum, val);
        }
    }
 
    return output;
}
 
template<typename Value_, typename Index_, typename Output_>
RangeDirectResult<Output_, Index_> range_direct(const Value_* value, const Index_ num_nonzero, const Index_ num_all, const RangeOptions<Output_>& opt) {
    if (num_nonzero) {
        auto candidate = range_direct(value, num_nonzero, opt);
        if (num_nonzero < num_all) {
            candidate.minimum = std::min(candidate.minimum, static_cast<Output_>(0));
            candidate.maximum = std::max(candidate.maximum, static_cast<Output_>(0));
            candidate.count += num_all - num_nonzero;
        }
        return candidate;
    } else if (num_all) {
        RangeDirectResult<Output_, Index_> output;
        output.minimum = 0;
        output.maximum = 0;
        output.count = num_all;
        return output;
    } else {
        RangeDirectResult<Output_, Index_> output;
        output.minimum = opt.minimum_placeholder;
        output.maximum = opt.maximum_placeholder;
        output.count = 0;
        return output;
    }
}
template<typename Output_, typename Count_>
struct RangeBuffers {
    Output_* minimum;
 
    Output_* maximum;
 
    Count_* count;
};
 
template<typename Value_, typename Index_, typename Output_, typename Count_>
void range_direct(bool row, const tatami::Matrix<Value_, Index_>& mat, RangeBuffers<Output_, Count_>& output, const RangeOptions<Output_>& opt) {
    const auto dim = (row ? mat.nrow() : mat.ncol());
    const auto otherdim = (row ? mat.ncol() : mat.nrow());
 
    if (mat.is_sparse()) {
        tatami::Options topt;
        topt.sparse_extract_index = false;
        tatami::parallelize([&](int, Index_ s, Index_ l) -> void {
            auto ext = tatami::consecutive_extractor<true>(mat, row, s, l, topt);
            auto vbuffer = tatami::create_container_of_Index_size<std::vector<Value_> >(otherdim);
            for (Index_ x = 0; x < l; ++x) {
                auto out = ext->fetch(vbuffer.data(), NULL);
                auto res = range_direct(out.value, out.number, otherdim, opt);
                output.minimum[x + s] = res.minimum;
                output.maximum[x + s] = res.maximum;
                output.count[x + s] = res.count;
            }
        }, dim, opt.num_threads);
 
    } else {
        tatami::parallelize([&](int, Index_ s, Index_ l) -> void {
            auto ext = tatami::consecutive_extractor<false>(mat, row, s, l);
            auto buffer = tatami::create_container_of_Index_size<std::vector<Value_> >(otherdim);
            for (Index_ x = 0; x < l; ++x) {
                auto ptr = ext->fetch(buffer.data());
                auto res = range_direct(ptr, otherdim, opt);
                output.minimum[x + s] = res.minimum;
                output.maximum[x + s] = res.maximum;
                output.count[x + s] = res.count;
            }
        }, dim, opt.num_threads);
    }
}
 
template<typename Value_, typename Index_, typename Output_, typename Count_>
void range_running(bool row, const tatami::Matrix<Value_, Index_>& mat, RangeBuffers<Output_, Count_>& output, const RangeOptions<Output_>& opt) {
    const auto dim = (row ? mat.nrow() : mat.ncol());
    const auto otherdim = (row ? mat.ncol() : mat.nrow());
    const bool is_sparse = mat.is_sparse();
 
    const bool do_parallel = opt.num_threads > 1; 
    std::optional<std::vector<std::optional<std::vector<Output_> > > > all_partial_min, all_partial_max;
    std::optional<std::vector<std::optional<std::vector<Count_> > > > all_partial_count;
    if (do_parallel) {
        all_partial_min.emplace(sanisizer::cast<I<decltype(all_partial_min->size())> >(opt.num_threads - 1));
        all_partial_max.emplace(sanisizer::cast<I<decltype(all_partial_max->size())> >(opt.num_threads - 1));
        all_partial_count.emplace(sanisizer::cast<I<decltype(all_partial_count->size())> >(opt.num_threads - 1));
    }
 
    std::fill_n(output.count, dim, 0);
 
    // If we're not skipping NaNs and we have at least one dimension element,
    // the output arrays will be fully populated when thread 0 processes the first dimension element.
    if (otherdim == 0) {
        std::fill_n(output.minimum, dim, opt.minimum_placeholder);
        std::fill_n(output.maximum, dim, opt.maximum_placeholder);
        return;
    }
 
    // No need to wipe dirty output buffers in the dense case, as we already set each entry of the output buffers.
    if (is_sparse) {
        std::fill_n(output.minimum, dim, 0);
        std::fill_n(output.maximum, dim, 0);
    }
 
    const auto nused = tatami::parallelize([&](int thread, Index_ s, Index_ l) -> void {
        Output_* min_ptr;
        Output_* max_ptr;
        Count_* count_ptr;
        std::optional<std::vector<Output_> > cur_min, cur_max;
        std::optional<std::vector<Count_> > cur_count;
        if (!do_parallel) {
            min_ptr = output.minimum;
            max_ptr = output.maximum;
            count_ptr = output.count;
        } else {
            if (thread == 0) {
                min_ptr = output.minimum;
                max_ptr = output.maximum;
                count_ptr = output.count;
            } else {
                cur_min.emplace(tatami::cast_Index_to_container_size<std::vector<Output_> >(dim));
                cur_max.emplace(tatami::cast_Index_to_container_size<std::vector<Output_> >(dim));
                cur_count.emplace(tatami::cast_Index_to_container_size<std::vector<Count_> >(dim));
                min_ptr = cur_min->data();
                max_ptr = cur_max->data();
                count_ptr = cur_count->data();
            }
        }
 
        if (is_sparse) {
            tatami::Options topt;
            topt.sparse_ordered_index = false;
            auto ext = tatami::consecutive_extractor<true>(mat, !row, s, l, topt);
            auto vbuffer = tatami::create_container_of_Index_size<std::vector<Value_> >(dim);
            auto ibuffer = tatami::create_container_of_Index_size<std::vector<Index_> >(dim);
            auto nonzeros = tatami::create_container_of_Index_size<std::vector<Index_> >(dim);
 
            for (Index_ x = 0; x < l; ++x) {
                auto out = ext->fetch(vbuffer.data(), ibuffer.data());
 
                // For the first observed vector in each thread, we can optimize it a little as we don't need to read existing min/max.
                if (x == 0) {
                    for (Index_ i = 0; i < out.number; ++i) {
                        const auto val = out.value[i];
                        const auto idx = out.index[i];
                        if (!std::isnan(val)) {
                            min_ptr[idx] = val;
                            max_ptr[idx] = val;
                            ++count_ptr[idx];
                        } else {
                            min_ptr[idx] = opt.minimum_placeholder;
                            max_ptr[idx] = opt.maximum_placeholder;
                        }
                        ++nonzeros[idx];
                    }
                } else {
                    for (Index_ i = 0; i < out.number; ++i) {
                        const auto val = out.value[i];
                        const auto idx = out.index[i];
                        if (!std::isnan(val)) {
                            auto& min_current = min_ptr[idx];
                            auto& max_current = max_ptr[idx];
                            if (count_ptr[idx] == 0) {
                                min_current = val;
                                max_current = val;
                            } else {
                                min_current = std::min(min_current, val);
                                max_current = std::max(max_current, val);
                            }
                            ++count_ptr[idx];
                        }
                        ++nonzeros[idx];
                    }
                }
            }
 
            for (Index_ d = 0; d < dim; ++d) {
                if (l > nonzeros[d]) {
                    count_ptr[d] += l - nonzeros[d];
                    auto& min_current = min_ptr[d];
                    min_current = std::min(min_current, static_cast<Output_>(0));
                    auto& max_current = max_ptr[d];
                    max_current = std::max(max_current, static_cast<Output_>(0));
                }
            }
 
        } else {
            auto ext = tatami::consecutive_extractor<false>(mat, !row, s, l);
            auto buffer = tatami::create_container_of_Index_size<std::vector<Value_> >(dim);
 
            for (Index_ x = 0; x < l; ++x) {
                auto ptr = ext->fetch(buffer.data());
 
                if (x == 0) {
                    // For the first observed vector in each thread,
                    // we can optimize it a little as we don't need to read existing min/max.
                    for (Index_ i = 0; i < dim; ++i) {
                        const auto val = ptr[i];
                        if (!std::isnan(val)) {
                            min_ptr[i] = val;
                            max_ptr[i] = val;
                            ++count_ptr[i];
                        } else {
                            min_ptr[i] = opt.minimum_placeholder;
                            max_ptr[i] = opt.maximum_placeholder;
                        }
                    }
                } else {
                    for (Index_ i = 0; i < dim; ++i) {
                        const auto val = ptr[i];
                        if (!std::isnan(val)) {
                            auto& min_current = min_ptr[i];
                            auto& max_current = max_ptr[i];
                            if (count_ptr[i] == 0) {
                                min_current = val;
                                max_current = val;
                            } else {
                                min_current = std::min(min_current, val);
                                max_current = std::max(max_current, val);
                            }
                            ++count_ptr[i];
                        }
                    }
                }
            }
        }
 
        if (do_parallel) {
            if (thread > 0) {
                (*all_partial_min)[thread - 1] = std::move(cur_min);
                (*all_partial_max)[thread - 1] = std::move(cur_max);
                (*all_partial_count)[thread - 1] = std::move(cur_count);
            }
        }
    }, otherdim, opt.num_threads);
 
    if (do_parallel) {
        for (int u = 1; u < nused; ++u) {
            const auto& cur_min = *((*all_partial_min)[u - 1]);
            const auto& cur_max = *((*all_partial_max)[u - 1]);
            const auto& cur_count = *((*all_partial_count)[u - 1]);
            for (Index_ d = 0; d < dim; ++d) {
                if (!cur_count[d]) {
                    continue;
                }
                if (output.count[d]) {
                    output.minimum[d] = std::min(cur_min[d], output.minimum[d]);
                    output.maximum[d] = std::max(cur_max[d], output.maximum[d]);
                } else {
                    output.minimum[d] = cur_min[d];
                    output.maximum[d] = cur_max[d];
                }
                output.count[d] += cur_count[d];
            }
        }
    }
}
template<typename Value_, typename Index_, typename Output_, typename Count_>
void range(bool row, const tatami::Matrix<Value_, Index_>& mat, RangeBuffers<Output_, Count_>& output, const RangeOptions<Output_>& opt) {
    if (mat.prefer_rows() == row) {
        range_direct(row, mat, output, opt);
    } else {
        range_running(row, mat, output, opt);
    }
}
 
template<typename Output_, typename Count_>
struct RangeResult {
    std::vector<Output_> minimum;
 
    std::vector<Output_> maximum;
 
    std::vector<Count_> count;
};
 
template<typename Value_, typename Index_, typename Output_ = Value_, typename Count_ = Index_>
RangeResult<Output_, Count_> range(bool row, const tatami::Matrix<Value_, Index_>& mat, const RangeOptions<Output_>& opt) {
    RangeResult<Output_, Count_> output;
    const auto dim = (row ? mat.nrow() : mat.ncol());
    tatami::resize_container_to_Index_size(output.minimum, dim
#ifdef TATAMI_STATS_TEST_DIRTY
        , -1
#endif
    );
    tatami::resize_container_to_Index_size(output.maximum, dim
#ifdef TATAMI_STATS_TEST_DIRTY
        , -1
#endif
    );
    tatami::resize_container_to_Index_size(output.count, dim
#ifdef TATAMI_STATS_TEST_DIRTY
        , -1
#endif
    );
 
    RangeBuffers<Output_, Count_> buffers;
    buffers.minimum = output.minimum.data();
    buffers.maximum = output.maximum.data();
    buffers.count = output.count.data();
    range(row, mat, buffers, opt);
 
    return output;
}
 
}
 
}
 
#endif