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?
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
#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);
}