range.hpp

Go to the documentation of this file.

#ifndef TATAMI_STATS_RANGE_HPP
#define TATAMI_STATS_RANGE_HPP
 
#include "utils.hpp"
 
#include <vector>
#include <algorithm>
#include <type_traits>
#include <limits>
 
#include "tatami/tatami.hpp"
 
namespace tatami_stats {
 
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 Value_, typename Index_, typename Output_>
Output_ min_direct(const Value_* const ptr, const Index_ num, const RangeOptions<Output_>& opt) {
    if (num) {
        return *std::min_element(ptr, ptr + num);
    } else {
        return opt.minimum_placeholder;
    }
}
 
template<typename Value_, typename Index_, typename Output_>
Output_ max_direct(const Value_* ptr, const Index_ num, const RangeOptions<Output_>& opt) {
    if (num) {
        return *std::max_element(ptr, ptr + num);
    } else {
        return opt.maximum_placeholder;
    }
}
 
template<typename Value_, typename Index_, typename Output_>
Output_ min_direct(const Value_* value, const Index_ num_nonzero, const Index_ num_all, const RangeOptions<Output_>& opt) {
    if (num_nonzero) {
        auto candidate = min_direct(value, num_nonzero, opt);
        if (num_nonzero < num_all) {
            if (candidate > 0) {
                candidate = 0;
            }
        }
        return candidate;
    } else if (num_all) {
        return 0;
    } else {
        return opt.minimum_placeholder;
    }
}
 
template<typename Value_, typename Index_, typename Output_>
Output_ max_direct(const Value_* value, const Index_ num_nonzero, const Index_ num_all, const RangeOptions<Output_>& opt) {
    if (num_nonzero) {
        auto candidate = max_direct(value, num_nonzero, opt);
        if (num_nonzero < num_all) {
            if (candidate < 0) {
                candidate = 0;
            }
        }
        return candidate;
    } else if (num_all) {
        return 0;
    } else {
        return opt.maximum_placeholder;
    }
}
template<typename Output_>
struct RangeBuffers {
    Output_* minimum;
 
    Output_* maximum;
};
 
template<typename Value_, typename Index_, typename Output_>
void range_direct(bool row, const tatami::Matrix<Value_, Index_>& mat, RangeBuffers<Output_>& 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);
                output.minimum[x + s] = min_direct(out.value, out.number, otherdim, opt);
                output.maximum[x + s] = max_direct(out.value, out.number, otherdim, opt);
            }
        }, 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());
                output.minimum[x + s] = min_direct(ptr, otherdim, opt);
                output.maximum[x + s] = max_direct(ptr, otherdim, opt);
            }
        }, dim, opt.num_threads);
    }
}
 
template<typename Value_, typename Index_, typename Output_>
void range_running(bool row, const tatami::Matrix<Value_, Index_>& mat, RangeBuffers<Output_>& 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;
    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));
    }
 
    if (otherdim == 0) {
        std::fill_n(output.minimum, dim, opt.minimum_placeholder);
        std::fill_n(output.maximum, dim, opt.maximum_placeholder);
        return;
    }
 
    const auto nused = tatami::parallelize([&](int thread, Index_ s, Index_ l) -> void {
        Output_* min_ptr;
        Output_* max_ptr;
        std::optional<std::vector<Output_> > cur_min, cur_max;
        if (!do_parallel) {
            min_ptr = output.minimum;
            max_ptr = output.maximum;
        } else {
            if (thread == 0) {
                min_ptr = output.minimum;
                max_ptr = output.maximum;
            } 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));
                min_ptr = cur_min->data();
                max_ptr = cur_max->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);
 
            // We pretend to start at one structural non-zero for every dimension element. 
            // This is because the first iteration effectively populates 'min_ptr/max_ptr' as a dense vector.
            // So even a structural zero is already considered here, being treated as a structural non-zero with a value of zero.
            auto nonzeros = tatami::create_container_of_Index_size<std::vector<Index_> >(dim, 1);
 
            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) {
                    // We treat the first extracted row/column as a dense vector, expanding it with all of the zeros.
                    // For non-main threads, our thread-local buffer is already zeroed so no need to explicitly do this. 
                    if (!do_parallel || thread == 0) {
                        std::fill_n(min_ptr, dim, 0);
                        std::fill_n(max_ptr, dim, 0);
                    }
                    for (Index_ i = 0; i < out.number; ++i) {
                        const auto val = out.value[i];
                        const auto idx = out.index[i];
                        min_ptr[idx] = val;
                        max_ptr[idx] = val;
                    }
                } else {
                    for (Index_ i = 0; i < out.number; ++i) {
                        const auto val = out.value[i];
                        const auto idx = out.index[i];
                        auto& min_current = min_ptr[idx];
                        min_current = std::min(min_current, val);
                        auto& max_current = max_ptr[idx];
                        max_current = std::max(max_current, val);
                        ++nonzeros[idx];
                    }
                }
            }
 
            for (Index_ d = 0; d < dim; ++d) {
                if (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());
 
                // 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) {
                    std::copy_n(ptr, dim, min_ptr);
                    std::copy_n(ptr, dim, max_ptr);
                } else {
                    for (Index_ i = 0; i < dim; ++i) {
                        const auto val = ptr[i];
                        auto& min_current = min_ptr[i];
                        min_current = std::min(min_current, val);
                        auto& max_current = max_ptr[i];
                        max_current = std::max(max_current, val);
                    }
                }
            }
        }
 
        if (do_parallel) {
            if (thread > 0) {
                (*all_partial_min)[thread - 1] = std::move(cur_min);
                (*all_partial_max)[thread - 1] = std::move(cur_max);
            }
        }
    }, 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]);
            for (Index_ d = 0; d < dim; ++d) {
                // All threads would have processed at least one element,
                // so we don't have to worry about dirty input buffers.
                output.minimum[d] = std::min(output.minimum[d], cur_min[d]);
                output.maximum[d] = std::max(output.maximum[d], cur_max[d]);
            }
        }
    }
}
template<typename Value_, typename Index_, typename Output_>
void range(bool row, const tatami::Matrix<Value_, Index_>& mat, RangeBuffers<Output_>& 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_>
struct RangeResult {
    std::vector<Output_> minimum;
 
    std::vector<Output_> maximum;
};
 
template<typename Value_, typename Index_, typename Output_ = Value_>
RangeResult<Output_> range(bool row, const tatami::Matrix<Value_, Index_>& mat, const RangeOptions<Output_>& opt) {
    RangeResult<Output_> 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
    );
 
    RangeBuffers<Output_> buffers;
    buffers.minimum = output.minimum.data();
    buffers.maximum = output.maximum.data();
    range(row, mat, buffers, opt);
 
    return output;
}
 
}
 
#endif