diff --git a/include/boost/multi/algorithms/ordering.hpp b/include/boost/multi/algorithms/ordering.hpp new file mode 100644 index 000000000..377554cbe --- /dev/null +++ b/include/boost/multi/algorithms/ordering.hpp @@ -0,0 +1,115 @@ +// Copyright 2026 Alfredo A. Correa +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +// `ordering` computes the permutation of indices that would sort a (multi)dimensional +// range *without moving the data* (a.k.a. argsort). The result is written into a +// caller-provided output range (no internal allocation). Because it orders by index, +// it handles non-zero-based arrays naturally and never copies/moves the (possibly +// proxy-referenced) elements. +// +// The sort backend is not hard-coded to `std::sort`: it dispatches to an ADL-found +// `sort` (e.g. `thrust::sort` when the output iterator is a thrust iterator) and falls +// back to `std::sort` otherwise. Overloads taking an execution policy as the first +// argument forward it (e.g. `std::execution::par`, or a thrust execution policy). +// NOTE: because dispatch is by ADL, `boost::multi` must not declare an iterator-triple +// `sort(first, last, comp)`, which would otherwise be selected for Multi iterators. + +// #pragma once +#ifndef BOOST_MULTI_ALGORITHMS_ORDERING_HPP +#define BOOST_MULTI_ALGORITHMS_ORDERING_HPP + +#include // for std::sort, std::copy +#include // for std::less +#include // for std::enable_if_t, std::void_t +#include // for std::forward, std::declval + +namespace boost::multi { + +namespace detail { +// detects a Multi array/subarray (has `.extension()`), used to tell an array argument +// apart from an execution-policy argument in the policy overloads below. +template struct has_extension : std::false_type {}; +template +struct has_extension().extension())>> : std::true_type {}; + +// overload-priority tag: priority_tag<1> is preferred over priority_tag<0>. +template struct priority_tag : priority_tag {}; +template<> struct priority_tag<0> {}; + +// Prefer an ADL-found `sort` (e.g. `thrust::sort` for thrust iterators) and fall back to +// `std::sort`. Crucially, `std::sort` is NOT brought into scope in the ADL overload, so a +// same-signature `thrust::sort` wins unambiguously instead of tying with `std::sort`. +template +auto sort_dispatch(priority_tag<1> /*prefer ADL*/, It first, It last, Compare comp) + -> decltype(sort(first, last, comp)) { + return sort(first, last, comp); // ADL only (no `using std::sort`) +} +template +auto sort_dispatch(priority_tag<0> /*fallback*/, It first, It last, Compare comp) -> void { + std::sort(first, last, comp); +} + +template +auto sort_dispatch(priority_tag<1> /*prefer ADL*/, Policy&& policy, It first, It last, Compare comp) + -> decltype(sort(std::forward(policy), first, last, comp)) { + return sort(std::forward(policy), first, last, comp); // ADL only (e.g. thrust policies) +} +template +auto sort_dispatch(priority_tag<0> /*fallback*/, Policy&& policy, It first, It last, Compare comp) -> void { + std::sort(std::forward(policy), first, last, comp); // std execution policies +} +} // namespace detail + +// Writes into [first, ...) the permutation of `arr`'s indices such that +// `arr[result[0]], arr[result[1]], ...` is non-decreasing according to `comp`. +// `first` must point to a mutable random-access range of at least `arr.size()` elements. +// `arr` is not modified. Returns the end of the written range. +template::value, int> = 0> // NOLINT(modernize-use-constraints) for C++17 +auto ordering(Array1D const& arr, RandomAccessIt first, Compare comp) -> RandomAccessIt { + auto const last = std::copy(arr.extension().begin(), arr.extension().end(), first); // seed output with the index values of `arr` + + detail::sort_dispatch( + detail::priority_tag<1>{}, first, last, + [&arr, comp](auto idx1, auto idx2) { return comp(arr[idx1], arr[idx2]); } + ); + + return last; +} + +template +auto ordering(Array1D const& arr, RandomAccessIt first) -> RandomAccessIt { + return ordering(arr, first, std::less<>{}); +} + +// Policy-aware overloads: `policy` (e.g. `std::execution::par`, or a thrust policy) is +// forwarded to the sort. Disambiguated from the overloads above by requiring the second +// argument to be a Multi array (an execution policy is not). +template::value, int> = 0> // NOLINT(modernize-use-constraints) for C++17 +auto ordering(Policy&& policy, Array1D const& arr, RandomAccessIt first, Compare comp) -> RandomAccessIt { + auto const last = std::copy(arr.extension().begin(), arr.extension().end(), first); + + detail::sort_dispatch( + detail::priority_tag<1>{}, std::forward(policy), first, last, + [&arr, comp](auto idx1, auto idx2) { return comp(arr[idx1], arr[idx2]); } + ); + + return last; +} + +template::value, int> = 0> // NOLINT(modernize-use-constraints) for C++17 +auto ordering(Policy&& policy, Array1D const& arr, RandomAccessIt first) -> RandomAccessIt { + return ordering(std::forward(policy), arr, first, std::less<>{}); +} + +// To *apply* the resulting `order` (or any index permutation) there is no need for a +// dedicated algorithm here: +// - out of place (gather): dst[k] = src[order[k]], i.e. +// std::transform(order.begin(), order.end(), dst.begin(), [&src](auto idx) { return src[idx]; }); +// (or, lazily, a view: order | std::views::transform([&src](auto idx) { return src[idx]; })) +// - in place: boost::algorithm::apply_permutation(src.begin(), src.end(), order.begin(), order.end()); +// Both work for N-dimensional `src`, where `src[idx]` is a whole slice. + +} // end namespace boost::multi + +#endif // BOOST_MULTI_ALGORITHMS_ORDERING_HPP diff --git a/test/ordering.cpp b/test/ordering.cpp new file mode 100644 index 000000000..3de18cafc --- /dev/null +++ b/test/ordering.cpp @@ -0,0 +1,133 @@ +// Copyright 2026 Alfredo A. Correa +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include // for ordering +#include // for array + +#include + +#include // for sort +#include // for greater +#include // for next + +namespace multi = boost::multi; + +namespace { +// a minimal user-defined "execution policy" with an ADL-findable sort, used to exercise +// the policy overload and its ADL dispatch portably (no /TBB, nvcc/nvhpc-safe). +struct seq_policy {}; +template +void sort(seq_policy /*policy*/, It first, It last, Compare comp) { std::sort(first, last, comp); } +} // namespace + +auto main() -> int { // NOLINT(readability-function-cognitive-complexity,bugprone-exception-escape) + // basic ascending order into a caller-provided output buffer + { + multi::array const arr = {3, 1, 2}; + + multi::array order(arr.extents()); + + BOOST_TEST( multi::ordering(arr, order.begin()) == order.end() ); + + BOOST_TEST( order[0] == 1 ); + BOOST_TEST( order[1] == 2 ); + BOOST_TEST( order[2] == 0 ); + + // the data must be untouched + BOOST_TEST( arr[0] == 3 ); + BOOST_TEST( arr[1] == 1 ); + BOOST_TEST( arr[2] == 2 ); + + // reading through the order gives a sorted view + BOOST_TEST( arr[order[0]] == 1 ); + BOOST_TEST( arr[order[1]] == 2 ); + BOOST_TEST( arr[order[2]] == 3 ); + } + + // descending order via custom comparator + { + multi::array const arr = {3, 1, 2}; + + multi::array order(arr.extents()); + + multi::ordering(arr, order.begin(), std::greater<>{}); + + BOOST_TEST( arr[order[0]] == 3 ); + BOOST_TEST( arr[order[1]] == 2 ); + BOOST_TEST( arr[order[2]] == 1 ); + } + + // already sorted -> identity permutation + { + multi::array const arr = {1, 2, 3, 4}; + + multi::array order(arr.extents()); + + multi::ordering(arr, order.begin()); + + BOOST_TEST( order[0] == 0 ); + BOOST_TEST( order[1] == 1 ); + BOOST_TEST( order[2] == 2 ); + BOOST_TEST( order[3] == 3 ); + } + + // floating-point elements; verify the order yields a sorted sequence + { + multi::array const arr = {2.5, -1.0, 0.0, 9.9, 3.3}; + + multi::array order(arr.extents()); + + BOOST_TEST( multi::ordering(arr, order.begin()) == order.end() ); + + for(multi::index k = 0; k + 1 != order.size(); ++k) { // NOLINT(altera-unroll-loops,altera-id-dependent-backward-branch) + BOOST_TEST( arr[order[k]] <= arr[order[k + 1]] ); + } + } + + // ties: std::sort is not stable, so only require the result to be sorted and a valid permutation + { + multi::array const arr = {2, 2, 1, 3, 1}; + + multi::array order(arr.extents()); + + multi::ordering(arr, order.begin()); + + for(multi::index k = 0; k + 1 != order.size(); ++k) { // NOLINT(altera-unroll-loops,altera-id-dependent-backward-branch) + BOOST_TEST( arr[order[k]] <= arr[order[k + 1]] ); + } + + multi::index sum = 0; + for(auto idx : order) { // NOLINT(altera-unroll-loops) + sum += idx; + } + BOOST_TEST( sum == (0 + 1 + 2 + 3 + 4) ); // it is a permutation of {0..4} + } + + // caller-provided raw pointer (the most general output) via data_elements() + { + multi::array const arr = {5, 4, 6}; + + multi::array buf(arr.extents()); + + BOOST_TEST( multi::ordering(arr, buf.data_elements()) == std::next(buf.data_elements(), 3) ); + + BOOST_TEST( arr[buf[0]] == 4 ); + BOOST_TEST( arr[buf[1]] == 5 ); + BOOST_TEST( arr[buf[2]] == 6 ); + } + + // policy overload + ADL sort dispatch (also exercises disambiguation from the policy-free overloads) + { + multi::array const arr = {3, 1, 2, 5, 4}; + + multi::array order(arr.extents()); + multi::ordering(seq_policy{}, arr, order.begin()); + + for(multi::index k = 0; k + 1 != order.size(); ++k) { // NOLINT(altera-unroll-loops,altera-id-dependent-backward-branch) + BOOST_TEST( arr[order[k]] <= arr[order[k + 1]] ); + } + } + + return boost::report_errors(); +}