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<bool minimum_, typename Value_>
constexpr auto choose_placeholder() {
    if constexpr(minimum_) {
        // 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();
        }
    } else {
        // 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<bool minimum_, typename Output_, typename Value_>
bool is_better(Output_ best, Value_ alt) {
    if constexpr(minimum_) {
        return best > static_cast<Output_>(alt);
    } else {
        return best < static_cast<Output_>(alt);
    }
}
 
}
template<bool minimum_, typename Value_, typename Index_>
Value_ direct(const Value_* ptr, Index_ num, bool skip_nan) {
    return ::tatami_stats::internal::nanable_ifelse_with_value<Value_>(
        skip_nan,
        [&]() -> Value_ {
            auto current = internal::choose_placeholder<minimum_, Value_>(); 
            for (Index_ i = 0; i < num; ++i) {
                auto val = ptr[i];
                if (internal::is_better<minimum_>(current, val)) { // no need to explicitly handle NaNs, as any comparison with NaNs is always false.
                    current = val;
                }
            }
            return current;
        },
        [&]() -> Value_ {
            if (num) {
                if constexpr(minimum_) {
                    return *std::min_element(ptr, ptr + num);
                } else {
                    return *std::max_element(ptr, ptr + num);
                }
            } else {
                return internal::choose_placeholder<minimum_, Value_>(); 
            }
        }
    );
}
 
template<bool minimum_, typename Value_, typename Index_>
Value_ direct(const Value_* value, Index_ num_nonzero, Index_ num_all, bool skip_nan) {
    if (num_nonzero) {
        auto candidate = direct<minimum_>(value, num_nonzero, skip_nan);
        if (num_nonzero < num_all && internal::is_better<minimum_>(candidate, 0)) {
            candidate = 0;
        }
        return candidate;
    } else if (num_all) {
        return 0;
    } else {
        return internal::choose_placeholder<minimum_, Value_>();
    }
}
 
template<bool minimum_, typename Output_, typename Value_, typename Index_>
class RunningDense {
public:
    RunningDense(Index_ num, Output_* store, bool skip_nan) : my_num(num), my_store(store), 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,
                [&]() {
                    for (Index_ i = 0; i < my_num; ++i, ++ptr) {
                        auto val = *ptr;
                        if (std::isnan(val)) {
                            my_store[i] = internal::choose_placeholder<minimum_, Value_>();
                        } else {
                            my_store[i] = val;
                        }
                    }
                },
                [&]() {
                    std::copy_n(ptr, my_num, my_store);
                }
            );
 
        } else {
            for (Index_ i = 0; i < my_num; ++i, ++ptr) {
                auto val = *ptr;
                if (internal::is_better<minimum_>(my_store[i], val)) { // this should implicitly skip NaNs, any NaN comparison will be false.
                    my_store[i] = val;
                }
            }
        }
    }
 
    void finish() {
        if (my_init) {
            std::fill_n(my_store, my_num, internal::choose_placeholder<minimum_, Value_>());
        }
    }
 
private:
    bool my_init = true;
    Index_ my_num;
    Output_* my_store;
    bool my_skip_nan;
};
 
template<bool minimum_, typename Output_, typename Value_, typename Index_>
class RunningSparse {
public:
    RunningSparse(Index_ num, Output_* store, bool skip_nan, Index_ subtract = 0) : 
        my_num(num), my_store(store), 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);
            std::fill_n(my_store, my_num, internal::choose_placeholder<minimum_, Value_>());
 
            if (!my_skip_nan) {
                for (Index_ i = 0; i < number; ++i, ++value, ++index) {
                    auto val = *value;
                    auto idx = *index - my_subtract;
                    my_store[idx] = val;
                    ++my_nonzero[idx];
                }
                my_count = 1;
                return;
            }
        }
 
        for (Index_ i = 0; i < number; ++i, ++value, ++index) {
            auto val = *value;
            auto idx = *index - my_subtract;
            auto& current = my_store[idx];
            if (internal::is_better<minimum_>(current, val)) { // this should implicitly skip NaNs, any NaN comparison will be false.
                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]) {
                    auto& current = my_store[i];
                    if (internal::is_better<minimum_>(current, 0)) {
                        current = 0;
                    }
                }
            }
        } else {
            std::fill_n(my_store, my_num, internal::choose_placeholder<minimum_, Value_>());
        }
    }
 
private:
    Index_ my_num;
    Output_* my_store;
    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_>* p, Output_* min_out, Output_* max_out, const Options& ropt) {
    auto dim = (row ? p->nrow() : p->ncol());
    auto otherdim = (row ? p->ncol() : p->nrow());
    const bool direct = p->prefer_rows() == row;
 
    bool store_min = min_out != NULL;
    bool store_max = max_out != NULL;
 
    if (p->sparse()) {
        tatami::Options opt;
        opt.sparse_ordered_index = false;
 
        if (direct) {
            opt.sparse_extract_index = false;
            tatami::parallelize([&](int, Index_ s, Index_ l) {
                auto ext = tatami::consecutive_extractor<true>(p, 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<true>(out.value, out.number, otherdim, ropt.skip_nan);
                    }
                    if (store_max) {
                        max_out[x + s] = ranges::direct<false>(out.value, out.number, otherdim, ropt.skip_nan);
                    }
                }
            }, dim, ropt.num_threads);
 
        } else {
            tatami::parallelize([&](int thread, Index_ s, Index_ l) {
                auto ext = tatami::consecutive_extractor<true>(p, !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<true, Output_, Value_, Index_> runmin(l, local_min.data(), ropt.skip_nan, s);
                ranges::RunningSparse<false, Output_, Value_, Index_> runmax(l, local_max.data(), ropt.skip_nan, s);
 
                for (Index_ x = 0; x < otherdim; ++x) {
                    auto out = ext->fetch(vbuffer.data(), ibuffer.data());
                    if (store_min) {
                        runmin.add(out.value, out.index, out.number);
                    }
                    if (store_max) {
                        runmax.add(out.value, out.index, out.number);
                    }
                }
 
                if (store_min) {
                    runmin.finish();
                    local_min.transfer();
                }
                if (store_max) {
                    runmax.finish();
                    local_max.transfer();
                }
            }, dim, ropt.num_threads);
        }
 
    } else {
        if (direct) {
            tatami::parallelize([&](int, Index_ s, Index_ l) {
                auto ext = tatami::consecutive_extractor<false>(p, 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<true>(ptr, otherdim, ropt.skip_nan);
                    }
                    if (store_max) {
                        max_out[x + s] = ranges::direct<false>(ptr, otherdim, ropt.skip_nan);
                    }
                }
            }, dim, ropt.num_threads);
 
        } else {
            tatami::parallelize([&](int thread, Index_ s, Index_ l) {
                auto ext = tatami::consecutive_extractor<false>(p, !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<true, Output_, Value_, Index_> runmin(l, local_min.data(), ropt.skip_nan);
                ranges::RunningDense<false, Output_, Value_, Index_> runmax(l, local_max.data(), ropt.skip_nan);
 
                for (Index_ x = 0; x < otherdim; ++x) {
                    auto ptr = ext->fetch(buffer.data());
                    if (store_min) {
                        runmin.add(ptr);
                    }
                    if (store_max) {
                        runmax.add(ptr);
                    }
                }
 
                if (store_min) {
                    runmin.finish();
                    local_min.transfer();
                }
                if (store_max) {
                    runmax.finish();
                    local_max.transfer();
                }
            }, dim, ropt.num_threads);
        }
    }
 
    return;
}
 
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) {
    std::vector<Output_> mins(p->ncol()), maxs(p->ncol());
    apply(false, p, mins.data(), maxs.data(), ropt);
    return std::make_pair(std::move(mins), std::move(maxs));
}
 
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, Options());
}
 
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) {
    std::vector<Output_> mins(p->nrow()), maxs(p->nrow());
    apply(true, p, mins.data(), maxs.data(), ropt);
    return std::make_pair(std::move(mins), std::move(maxs));
}
 
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, Options());
}
 
}
 
}
 
#endif