1
votes

In my program I use boost::multi_array and sometimes it is important to convert a pointer to an element inside the container into its indices, for example, to check that the element is not on the boundary of the array by any dimension. Here is oversimplified code just to explain the intension:

    boost::multi_array<int, 2> arr2d(boost::extents[10][10]);
    const auto data = arr2d.data();
    const auto end = data + arr2d.num_elements();

    for ( auto it = data; it != end; ++it )
    {
        [[maybe_unused]] int v = *it;
        // how to convert it into x and y indices in arr2d?
    } 

Is there any build-in valid means in boost::multi_array for pointer-to-indices conversion? Or the only possibility is manual division of the offset relative to data start on the dimensions?

1
The normal way to get the number of elements between two pointers (both pointing into the same array), one would subtract them. To get the index, one of the pointers should be the pointer to the first element. I don't know how this translates to a boost::multi_array.Some programmer dude

1 Answers

1
votes

This facility does not seem to be in the library. Which is a bit of a shame, but makes sense from the interface/use cases for the library.

I came up with this solution which takes into account all storage orderings/directions, and base indexes:

template <typename MultiArray,
          typename Element = typename MultiArray::element_type>
auto get_coords(Element const* el, MultiArray const& arr)
{
    using index            = typename MultiArray::index;
    using size_type        = typename MultiArray::size_type;
    auto constexpr NumDims = MultiArray::dimensionality;

    auto raw_index = std::distance(arr.data(), el);
    assert(raw_index >= 0 && size_t(raw_index) < arr.num_elements());

    std::array<index, NumDims> coord {0};
    auto shape    = arr.shape();
    auto strides  = arr.strides();
    auto bases    = arr.index_bases();
    auto& storage = arr.storage_order();

    for (size_type n = 0; n < NumDims; ++n) {
        auto dim = storage.ordering(NumDims - 1 - n); // reverse order
        //fmt::print("dim: {} stride: {}\n", dim, strides[dim]);
        coord[dim] = raw_index / strides[dim];
        raw_index -= coord[dim] * strides[dim];
    }

    for (size_type dim = 0; dim < NumDims; ++dim) {
        coord[dim] += bases[dim];
        if (!storage.ascending(dim))
            coord[dim] += shape[dim] - 1;
    }

    return coord;
}

Because this is way too much logic not to test it I came up with a rigorous test:

void test(auto arr) {
    iota_n(arr.data(), 1, arr.num_elements());
    fmt::print("\n{} -> {}\n", arr, map(arr));
    check_all(arr);
}

This causes the mapped coordinates to be printed for every element in memory order:

template <typename MultiArray> auto map(MultiArray const& arr)
{
    using Coord =
        std::array<typename MultiArray::index, MultiArray::dimensionality>;

    std::vector<Coord> v;
    for (size_t n = 0; n<arr.num_elements(); ++n) {
        v.push_back(get_coords(arr.data() + n, arr));
    }

    return v;
}

And also checks all elements:

void check_element(auto const& el, auto& arr) {
    auto coord = get_coords(&el, arr);
    //fmt::print("{} at {}\n", el, coord);

    // assert same value at same address
    auto& check = deref(arr, coord);
    assert(el == check);
    assert(&el == &check);
    s_elements_checked += 1;
}

void check_all(int    const& el, auto& arr) { check_element(el, arr); }
void check_all(double const& el, auto& arr) { check_element(el, arr); }

void check_all(auto const& rank, auto& arr) {
    for (auto&& rank_or_val : rank) check_all(rank_or_val, arr);
}
void check_all(auto const& arr) {
    for (auto&& rank : arr) check_all(rank, arr);
}

Then I ran it on a variety of test matrices:

int main() {
    // default
    test(boost::multi_array<int, 2>{boost::extents[3][2], boost::c_storage_order()});
    // FORTRAN
    test(boost::multi_array<int, 2>{boost::extents[3][2], boost::fortran_storage_order()});

    {
        boost::multi_array<int, 3> based{boost::extents[3][2][3], boost::fortran_storage_order()};
        // using standard bases (0)
        test(based);

        // using non-standard bases
        based.reindex(std::array<int, 3> {50,-20,100});
        test(based);
    }

    {
        // with custom storage ordering of dimensions, including descending
        std::array<int, 3> order{2, 0, 1};
        std::array<bool, 3> ascending {false,true,true};
        boost::multi_array<int, 3> custom{
            boost::extents[3][2][4],
            boost::general_storage_order<3>(order.data(), ascending.data())};
        custom.reindex(std::array<int, 3> {100,100,100});
        test(custom);
    }

    fmt::print("Checked a total of {} elements, verified ok\n", s_elements_checked);
}

This prints

{{1, 2}, {3, 4}, {5, 6}} -> {{0, 0}, {0, 1}, {1, 0}, {1, 1}, {2, 0}, {2, 1}}

{{1, 4}, {2, 5}, {3, 6}} -> {{0, 0}, {1, 0}, {2, 0}, {0, 1}, {1, 1}, {2, 1}}

{{{1, 7, 13}, {4, 10, 16}}, {{2, 8, 14}, {5, 11, 17}}, {{3, 9, 15}, {6, 12, 18}}} -> {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {0, 1, 0}, {1, 1, 0}, {2, 1, 0}, {0, 0, 1}, {1, 0, 1}, {2, 0, 1}, {0, 1, 1}, {1, 1, 1}, {2, 1, 1}, {0, 0, 2}, {1, 0, 2}, {2, 0, 2}, {0, 1, 2}, {1, 1, 2}, {2, 1, 2}}

{{{1, 7, 13}, {4, 10, 16}}, {{2, 8, 14}, {5, 11, 17}}, {{3, 9, 15}, {6, 12, 18}}} -> {{50, -20, 100}, {51, -20, 100}, {52, -20, 100}, {50, -19, 100}, {51, -19, 100}, {52, -19, 100}, {50, -20, 101}, {51, -20, 101}, {52, -20, 101}, {50, -19, 101}, {51, -19, 101}, {52, -19, 101}, {50, -20, 102}, {51, -20, 102}, {52, -20, 102}, {50, -19, 102}, {51, -19, 102}, {52, -19, 102}}

{{{9, 10, 11, 12}, {21, 22, 23, 24}}, {{5, 6, 7, 8}, {17, 18, 19, 20}}, {{1, 2, 3, 4}, {13, 14, 15, 16}}} -> {{102, 100, 100}, {102, 100, 101}, {102, 100, 102}, {102, 100, 103}, {101, 100, 100}, {101, 100, 101}, {101, 100, 102}, {101, 100, 103}, {100, 100, 100}, {100, 100, 101}, {100, 100, 102}, {100, 100, 103}, {102, 101, 100}, {102, 101, 101}, {102, 101, 102}, {102, 101, 103}, {101, 101, 100}, {101, 101, 101}, {101, 101, 102}, {101, 101, 103}, {100, 101, 100}, {100, 101, 101}, {100, 101, 102}, {100, 101, 103}} Checked a total of 72 elements, verified ok

Full Listing

Live On Compiler Explorer

#include <boost/multi_array.hpp>
#include <fmt/ranges.h>

template <typename MultiArray,
        typename Element = typename MultiArray::element_type>
auto get_coords(Element const* el, MultiArray const& arr)
{
    using index            = typename MultiArray::index;
    using size_type        = typename MultiArray::size_type;
    auto constexpr NumDims = MultiArray::dimensionality;

    auto raw_index = std::distance(arr.data(), el);
    assert(raw_index >= 0 && size_t(raw_index) < arr.num_elements());

    std::array<index, NumDims> coord {0};
    auto shape    = arr.shape();
    auto strides  = arr.strides();
    auto bases    = arr.index_bases();
    auto& storage = arr.storage_order();

    for (size_type n = 0; n < NumDims; ++n) {
        auto dim = storage.ordering(NumDims - 1 - n); // reverse order
        //fmt::print("dim: {} stride: {}\n", dim, strides[dim]);
        coord[dim] = raw_index / strides[dim];
        raw_index -= coord[dim] * strides[dim];
    }

    for (size_type dim = 0; dim < NumDims; ++dim) {
        coord[dim] += bases[dim];
        if (!storage.ascending(dim))
            coord[dim] += shape[dim] - 1;
    }

    return coord;
}

template <typename MultiArray> auto map(MultiArray const& arr)
{
    using Coord =
        std::array<typename MultiArray::index, MultiArray::dimensionality>;

    std::vector<Coord> v;
    for (size_t n = 0; n<arr.num_elements(); ++n) {
        v.push_back(get_coords(arr.data() + n, arr));
    }

    return v;
}

#include <boost/algorithm/cxx11/iota.hpp>
using boost::algorithm::iota_n;

auto& deref(auto const& arr, std::array<long, 1> coord) { return arr[coord[0]]; }
auto& deref(auto const& arr, std::array<long, 2> coord) { return arr[coord[0]][coord[1]]; }
auto& deref(auto const& arr, std::array<long, 3> coord) { return arr[coord[0]][coord[1]][coord[2]]; }
auto& deref(auto const& arr, std::array<long, 4> coord) { return arr[coord[0]][coord[1]][coord[2]][coord[3]]; }

static int s_elements_checked = 0;

void check_element(auto const& el, auto& arr) {
    auto coord = get_coords(&el, arr);
    //fmt::print("{} at {}\n", el, coord);

    // assert same value at same address
    auto& check = deref(arr, coord);
    assert(el == check);
    assert(&el == &check);
    s_elements_checked += 1;
}

void check_all(int    const& el, auto& arr) { check_element(el, arr); }
void check_all(double const& el, auto& arr) { check_element(el, arr); }

void check_all(auto const& rank, auto& arr) {
    for (auto&& rank_or_val : rank) check_all(rank_or_val, arr);
}
void check_all(auto const& arr) {
    for (auto&& rank : arr) check_all(rank, arr);
}

void test(auto arr) {
    iota_n(arr.data(), 1, arr.num_elements());
    fmt::print("\n{} -> {}\n", arr, map(arr));
    check_all(arr);
}

int main() {
    // default
    test(boost::multi_array<int, 2>{boost::extents[3][2], boost::c_storage_order()});
    // FORTRAN
    test(boost::multi_array<int, 2>{boost::extents[3][2], boost::fortran_storage_order()});

    {
        boost::multi_array<int, 3> based{boost::extents[3][2][3], boost::fortran_storage_order()};
        // using standard bases (0)
        test(based);

        // using non-standard bases
        based.reindex(std::array<int, 3> {50,-20,100});
        test(based);
    }

    {
        // with custom storage ordering of dimensions, including descending
        std::array<int, 3> order{2, 0, 1};
        std::array<bool, 3> ascending {false,true,true};
        boost::multi_array<int, 3> custom{
            boost::extents[3][2][4],
            boost::general_storage_order<3>(order.data(), ascending.data())};
        custom.reindex(std::array<int, 3> {100,100,100});
        test(custom);
    }

    fmt::print("Checked a total of {} elements, verified ok\n", s_elements_checked);
}