ranges.hpp

Go to the documentation of this file.

#ifndef TATAMI_STATS_RANGES_HPP
#define TATAMI_STATS_RANGES_HPP
 
#include "tatami/tatami.hpp"
#include "utils.hpp"
 
#include <vector>
#include <algorithm>
#include <type_traits>
 
namespace tatami_stats {
 
namespace ranges {
 
struct Options {
    bool skip_nan = false;
 
    int num_threads = 1;
};
 
namespace internal {
 
template<typename Value_>
constexpr Value_ choose_minimum_placeholder() {
    // Placeholder value 'x' is such that 'x >= y' is always true for any non-NaN 'y'.
    if constexpr(std::numeric_limits<Value_>::has_infinity) {
        return std::numeric_limits<Value_>::infinity();
    } else {
        return std::numeric_limits<Value_>::max();
    }
}
 
template<typename Value_>
constexpr Value_ choose_maximum_placeholder() {
    // Placeholder value 'x' is such that 'x <= y' is always true for any non-NaN 'y'.
    if constexpr(std::numeric_limits<Value_>::has_infinity) {
        return -std::numeric_limits<Value_>::infinity();
    } else {
        return std::numeric_limits<Value_>::lowest();
    }
}
 
}
template<typename Value_, typename Index_>
Value_ direct(const Value_* ptr, Index_ num, bool minimum, bool skip_nan) {
    return ::tatami_stats::internal::nanable_ifelse_with_value<Value_>(
        skip_nan,
 
        [&]() -> Value_ {
            if (minimum) {
                auto current = internal::choose_minimum_placeholder<Value_>(); 
                for (Index_ i = 0; i < num; ++i) {
                    auto val = ptr[i];
                    if (val < current) { // no need to explicitly handle NaNs, as any comparison with NaNs is always false.
                        current = val;
                    }
                }
                return current;
            } else {
                auto current = internal::choose_maximum_placeholder<Value_>(); 
                for (Index_ i = 0; i < num; ++i) {
                    auto val = ptr[i];
                    if (val > current) { // again, no need to explicitly handle NaNs, as any comparison with NaNs is always false.
                        current = val;
                    }
                }
                return current;
            }
        },
 
        [&]() -> Value_ {
            if (num) {
                if (minimum) {
                    return *std::min_element(ptr, ptr + num);
                } else {
                    return *std::max_element(ptr, ptr + num);
                }
            } else {
                if (minimum) {
                    return internal::choose_minimum_placeholder<Value_>();
                } else {
                    return internal::choose_maximum_placeholder<Value_>();
                }
            }
        }
    );
}
 
template<typename Value_, typename Index_>
Value_ direct(const Value_* value, Index_ num_nonzero, Index_ num_all, bool minimum, bool skip_nan) {
    if (num_nonzero) {
        auto candidate = direct(value, num_nonzero, minimum, skip_nan);
        if (num_nonzero < num_all) {
            if (minimum) {
                if (candidate > 0) {
                    candidate = 0;
                }
            } else {
                if (candidate < 0) {
                    candidate = 0;
                }
            }
        }
        return candidate;
 
    } else if (num_all) {
        return 0;
 
    } else {
        if (minimum) {
            return internal::choose_minimum_placeholder<Value_>();
        } else {
            return internal::choose_maximum_placeholder<Value_>();
        }
    }
}
 
template<typename Output_, typename Value_, typename Index_>
class RunningDense {
public:
    RunningDense(Index_ num, Output_* store_min, Output_* store_max, bool skip_nan) :
        my_num(num), my_store_min(store_min), my_store_max(store_max), my_skip_nan(skip_nan) {}
 
    void add(const Value_* ptr) {
        if (my_init) {
            my_init = false;
            ::tatami_stats::internal::nanable_ifelse<Value_>(
                my_skip_nan,
 
                [&]() -> void {
                    if (my_store_min) {
                        for (Index_ i = 0; i < my_num; ++i) {
                            auto val = ptr[i];
                            if (std::isnan(val)) {
                                my_store_min[i] = internal::choose_minimum_placeholder<Value_>();
                            } else {
                                my_store_min[i] = val;
                            }
                        }
                    }
                    if (my_store_max) {
                        for (Index_ i = 0; i < my_num; ++i) {
                            auto val = ptr[i];
                            if (std::isnan(val)) {
                                my_store_max[i] = internal::choose_maximum_placeholder<Value_>();
                            } else {
                                my_store_max[i] = val;
                            }
                        }
                    }
                },
 
                [&]() -> void {
                    if (my_store_min) {
                        std::copy_n(ptr, my_num, my_store_min);
                    }
                    if (my_store_max) {
                        std::copy_n(ptr, my_num, my_store_max);
                    }
                }
            );
 
        } else {
            if (my_store_min) {
                for (Index_ i = 0; i < my_num; ++i) {
                    auto val = ptr[i];
                    auto& current = my_store_min[i];
                    if (val < current) { // this should implicitly skip val=NaNs, as any NaN comparison will be false.
                        current = val;
                    }
                }
            }
            if (my_store_max) {
                for (Index_ i = 0; i < my_num; ++i) {
                    auto val = ptr[i];
                    auto& current = my_store_max[i];
                    if (val > current) { // this should implicitly skip val=NaNs, as any NaN comparison will be false.
                        current = val;
                    }
                }
            }
        }
    }
 
    void finish() {
        if (my_init) {
            if (my_store_min) {
                std::fill_n(my_store_min, my_num, internal::choose_minimum_placeholder<Value_>());
            }
            if (my_store_max) {
                std::fill_n(my_store_max, my_num, internal::choose_maximum_placeholder<Value_>());
            }
        }
    }
 
private:
    bool my_init = true;
    Index_ my_num;
    Output_* my_store_min;
    Output_* my_store_max;
    bool my_skip_nan;
};
 
template<typename Output_, typename Value_, typename Index_>
class RunningSparse {
public:
    RunningSparse(Index_ num, Output_* store_min, Output_* store_max, bool skip_nan, Index_ subtract = 0) : 
        my_num(num), my_store_min(store_min), my_store_max(store_max), my_skip_nan(skip_nan), my_subtract(subtract) {}
 
    void add(const Value_* value, const Index_* index, Index_ number) {
        if (my_count == 0) {
            my_nonzero.resize(my_num);
 
            if (my_store_min) {
                std::fill_n(my_store_min, my_num, internal::choose_minimum_placeholder<Value_>());
            }
            if (my_store_max) {
                std::fill_n(my_store_max, my_num, internal::choose_maximum_placeholder<Value_>());
            }
 
            if (!my_skip_nan) {
                for (Index_ i = 0; i < number; ++i) {
                    auto val = value[i];
                    auto idx = index[i] - my_subtract;
                    if (my_store_min) {
                        my_store_min[idx] = val;
                    }
                    if (my_store_max) {
                        my_store_max[idx] = val;
                    }
                    ++my_nonzero[idx];
                }
                my_count = 1;
                return;
            }
        }
 
        for (Index_ i = 0; i < number; ++i) {
            auto val = value[i];
            auto idx = index[i] - my_subtract;
            if (my_store_min) { // this should implicitly skip NaNs, any NaN comparison will be false.
                auto& current = my_store_min[idx];
                if (current > val) {
                    current = val;
                }
            }
            if (my_store_max) {
                auto& current = my_store_max[idx];
                if (current < val) {
                    current = val;
                }
            }
            ++my_nonzero[idx];
        }
 
        ++my_count;
    }
 
    void finish() {
        if (my_count) {
            for (Index_ i = 0; i < my_num; ++i) {
                if (my_count > my_nonzero[i]) {
                    if (my_store_min) {
                        auto& current = my_store_min[i];
                        if (current > 0) {
                            current = 0;
                        }
                    }
                    if (my_store_max) {
                        auto& current = my_store_max[i];
                        if (current < 0) {
                            current = 0;
                        }
                    }
                }
            }
        } else {
            if (my_store_min) {
                std::fill_n(my_store_min, my_num, internal::choose_minimum_placeholder<Value_>());
            }
            if (my_store_max) {
                std::fill_n(my_store_max, my_num, internal::choose_maximum_placeholder<Value_>());
            }
        }
    }
 
private:
    Index_ my_num;
    Output_* my_store_min;
    Output_* my_store_max;
    bool my_skip_nan;
    Index_ my_subtract;
    Index_ my_count = 0;
    std::vector<Index_> my_nonzero;
};
 
template<typename Value_, typename Index_, typename Output_>
void apply(bool row, const tatami::Matrix<Value_, Index_>& mat, Output_* min_out, Output_* max_out, const Options& ropt) {
    auto dim = (row ? mat.nrow() : mat.ncol());
    auto otherdim = (row ? mat.ncol() : mat.nrow());
    const bool direct = mat.prefer_rows() == row;
 
    bool store_min = min_out != NULL;
    bool store_max = max_out != NULL;
 
    if (mat.sparse()) {
        tatami::Options opt;
        opt.sparse_ordered_index = false;
 
        if (direct) {
            opt.sparse_extract_index = false;
            tatami::parallelize([&](int, Index_ s, Index_ l) -> void {
                auto ext = tatami::consecutive_extractor<true>(mat, row, s, l, opt);
                std::vector<Value_> vbuffer(otherdim);
                for (Index_ x = 0; x < l; ++x) {
                    auto out = ext->fetch(vbuffer.data(), NULL);
                    if (store_min) {
                        min_out[x + s] = ranges::direct(out.value, out.number, otherdim, true, ropt.skip_nan);
                    }
                    if (store_max) {
                        max_out[x + s] = ranges::direct(out.value, out.number, otherdim, false, ropt.skip_nan);
                    }
                }
            }, dim, ropt.num_threads);
 
        } else {
            tatami::parallelize([&](int thread, Index_ s, Index_ l) -> void {
                auto ext = tatami::consecutive_extractor<true>(mat, !row, static_cast<Index_>(0), otherdim, s, l, opt);
                std::vector<Value_> vbuffer(l);
                std::vector<Index_> ibuffer(l);
 
                auto local_min = (store_min ? LocalOutputBuffer<Output_>(thread, s, l, min_out) : LocalOutputBuffer<Output_>());
                auto local_max = (store_max ? LocalOutputBuffer<Output_>(thread, s, l, max_out) : LocalOutputBuffer<Output_>());
                ranges::RunningSparse<Output_, Value_, Index_> runner(l, local_min.data(), local_max.data(), ropt.skip_nan, s);
 
                for (Index_ x = 0; x < otherdim; ++x) {
                    auto out = ext->fetch(vbuffer.data(), ibuffer.data());
                    runner.add(out.value, out.index, out.number);
                }
 
                runner.finish();
                if (store_min) {
                    local_min.transfer();
                }
                if (store_max) {
                    local_max.transfer();
                }
            }, dim, ropt.num_threads);
        }
 
    } else {
        if (direct) {
            tatami::parallelize([&](int, Index_ s, Index_ l) -> void {
                auto ext = tatami::consecutive_extractor<false>(mat, row, s, l);
                std::vector<Value_> buffer(otherdim);
                for (Index_ x = 0; x < l; ++x) {
                    auto ptr = ext->fetch(buffer.data());
                    if (store_min) {
                        min_out[x + s] = ranges::direct(ptr, otherdim, true, ropt.skip_nan);
                    }
                    if (store_max) {
                        max_out[x + s] = ranges::direct(ptr, otherdim, false, ropt.skip_nan);
                    }
                }
            }, dim, ropt.num_threads);
 
        } else {
            tatami::parallelize([&](int thread, Index_ s, Index_ l) -> void {
                auto ext = tatami::consecutive_extractor<false>(mat, !row, static_cast<Index_>(0), otherdim, s, l);
                std::vector<Value_> buffer(l);
 
                auto local_min = (store_min ? LocalOutputBuffer<Output_>(thread, s, l, min_out) : LocalOutputBuffer<Output_>());
                auto local_max = (store_max ? LocalOutputBuffer<Output_>(thread, s, l, max_out) : LocalOutputBuffer<Output_>());
                ranges::RunningDense<Output_, Value_, Index_> runner(l, local_min.data(), local_max.data(), ropt.skip_nan);
 
                for (Index_ x = 0; x < otherdim; ++x) {
                    auto ptr = ext->fetch(buffer.data());
                    runner.add(ptr);
                }
 
                runner.finish();
                if (store_min) {
                    local_min.transfer();
                }
                if (store_max) {
                    local_max.transfer();
                }
            }, dim, ropt.num_threads);
        }
    }
 
    return;
}
 
// Back-compatibility.
template<typename Value_, typename Index_, typename Output_>
void apply(bool row, const tatami::Matrix<Value_, Index_>* p, Output_* min_out, Output_* max_out, const Options& ropt) {
    apply(row, *p, min_out, max_out, ropt); 
}
template<typename Output_ = double, typename Value_, typename Index_>
std::pair<std::vector<Output_>, std::vector<Output_> > by_column(const tatami::Matrix<Value_, Index_>& mat, const Options& ropt) {
    std::vector<Output_> mins(mat.ncol()), maxs(mat.ncol());
    apply(false, mat, mins.data(), maxs.data(), ropt);
    return std::make_pair(std::move(mins), std::move(maxs));
}
 
// Back-compatibility.
template<typename Output_ = double, typename Value_, typename Index_>
std::pair<std::vector<Output_>, std::vector<Output_> > by_column(const tatami::Matrix<Value_, Index_>* p, const Options& ropt) {
    return by_column<Output_>(*p, ropt);
}
 
template<typename Output_ = double, typename Value_, typename Index_>
std::pair<std::vector<Output_>, std::vector<Output_> > by_column(const tatami::Matrix<Value_, Index_>& mat) {
    return by_column<Output_>(mat, {});
}
 
template<typename Output_ = double, typename Value_, typename Index_>
std::pair<std::vector<Output_>, std::vector<Output_> > by_column(const tatami::Matrix<Value_, Index_>* p) {
    return by_column<Output_>(*p);
}
template<typename Output_ = double, typename Value_, typename Index_>
std::pair<std::vector<Output_>, std::vector<Output_> > by_row(const tatami::Matrix<Value_, Index_>& mat, const Options& ropt) {
    std::vector<Output_> mins(mat.nrow()), maxs(mat.nrow());
    apply(true, mat, mins.data(), maxs.data(), ropt);
    return std::make_pair(std::move(mins), std::move(maxs));
}
 
// Back-compatibility.
template<typename Output_ = double, typename Value_, typename Index_>
std::pair<std::vector<Output_>, std::vector<Output_> > by_row(const tatami::Matrix<Value_, Index_>* p, const Options& ropt) {
    return by_row<Output_>(*p, ropt);
}
 
template<typename Output_ = double, typename Value_, typename Index_>
std::pair<std::vector<Output_>, std::vector<Output_> > by_row(const tatami::Matrix<Value_, Index_>& mat) {
    return by_row<Output_>(mat, {}); 
}
 
template<typename Output_ = double, typename Value_, typename Index_>
std::pair<std::vector<Output_>, std::vector<Output_> > by_row(const tatami::Matrix<Value_, Index_>* p) {
    return by_row<Output_>(*p);
}
}
 
}
 
#endif