N dimensional array index utility

Question

I created the below class to help working with ND Arrays, mapping based on this question. This will help in implementing the code for handling convolutions.

How can I improve upon this? Is there something essential I'm missing from the interface?

header:

class NDArrayIndex{
public:
  NDArrayIndex(
    std::initializer_list<std::uint32_t> dimensions, std::int32_t padding = 0, 
    std::initializer_list<std::uint32_t> position = {}
  );

  NDArrayIndex& set(const std::vector<std::uint32_t>& position);
  NDArrayIndex& step();
  NDArrayIndex& step(std::uint32_t dimension, std::int32_t delta);
  const std::vector<std::uint32_t>& position() const{
    return m_position;
  }
  std::optional<std::uint32_t> calculate_mapped_position(const std::vector<std::uint32_t>& position) const;
  std::optional<std::uint32_t> mapped_position() const{
    return m_mappedIndex;
  }
  bool inside_bounds(const std::vector<std::uint32_t>& position, std::uint32_t dimension = 0u, std::int32_t delta = 0) const;
  bool inside_bounds(std::uint32_t dimension = 0u, std::int32_t delta = 0) const{
    return inside_bounds(m_position, dimension, delta);
  }
  bool inside_bounds(const NDArrayIndex& index, std::uint32_t dimension = 0u, std::int32_t delta = 0) const{
    return inside_bounds(index.position(), dimension, delta);
  }
  bool inside_content(const std::vector<std::uint32_t>& position, std::uint32_t dimension = 0u, std::int32_t delta = 0) const;
  bool inside_content(std::uint32_t dimension = 0u, std::int32_t delta = 0) const{
    return inside_content(m_position, dimension, delta);
  }
  bool inside_content(const NDArrayIndex& index, std::uint32_t dimension = 0u, std::int32_t delta = 0) const{
    return inside_content(index.position(), dimension, delta);
  }

  using IntervalPart = std::pair<std::uint32_t, std::uint32_t>;
  std::vector<IntervalPart> mappable_parts_of(std::uint32_t dimension, std::int32_t delta) const{
    return mappable_parts_of(m_position, dimension, delta);
  }
  std::vector<IntervalPart> mappable_parts_of(
    const std::vector<std::uint32_t>& position, std::uint32_t dimension, std::int32_t delta
  ) const;

  std::uint32_t buffer_size(){
    return m_bufferSize;
  }

private:
  const std::vector<std::uint32_t> m_dimensions;
  const std::int32_t m_padding;
  const std::vector<std::uint32_t> m_strides;
  const std::uint32_t m_bufferSize;
  std::vector<std::uint32_t> m_position;
  std::optional<std::uint32_t> m_mappedIndex;
};

source:

NDArrayIndex::NDArrayIndex(
  std::initializer_list<std::uint32_t> dimensions, std::int32_t padding, 
  std::initializer_list<std::uint32_t> position
)
: m_dimensions(dimensions)
, m_padding(padding)
, m_strides(init_strides(dimensions, m_padding))
, m_bufferSize(std::accumulate(m_dimensions.begin(), m_dimensions.end(), 1.0, 
  [](const std::uint32_t& partial, const std::uint32_t& element){ return partial * element; }
))
, m_position(init_position(m_dimensions, position))
, m_mappedIndex(calculate_mapped_position(m_position))
{
  assert(0 == std::count(m_dimensions.begin(), m_dimensions.end(), 0));
  assert(inside_bounds(m_position));
}

NDArrayIndex& NDArrayIndex::set(const std::vector<std::uint32_t>& position){
  assert(position.size() == m_position.size());
  assert(inside_bounds(position));
  m_position = position;
  m_mappedIndex = calculate_mapped_position(m_position);
  assert( (!m_mappedIndex.has_value())||(m_mappedIndex.value() < m_bufferSize) );
  return *this;
}

NDArrayIndex& NDArrayIndex::step(){
  std::uint32_t dim = 0;
  bool changed = false;
  while(dim < m_dimensions.size()){
    if(inside_bounds(dim, 1)){
      step(dim, 1);
      break;
    }else{
      changed = true;
      m_position[dim] = 0;
    }
    ++dim;
  }
  if(dim >= m_dimensions.size()){
    m_mappedIndex = 0; /* Overflow happened, start from the beginning */
  }else{
    if(changed)m_mappedIndex = calculate_mapped_position(m_position);
    assert(m_mappedIndex < m_bufferSize);
  }
  return *this;
}

NDArrayIndex& NDArrayIndex::step(std::uint32_t dimension, std::int32_t delta){
  const std::int32_t new_position = static_cast<std::int32_t>(m_position[dimension]) + delta;
  assert(0 <= new_position);
  assert((m_dimensions[dimension] + (2 * std::max(0, m_padding))) > static_cast<std::uint32_t>(new_position));

  m_position[dimension] = new_position;

  bool new_position_is_inside_content = inside_content(m_position);
  if(m_mappedIndex.has_value() && new_position_is_inside_content){ /* m_mappedIndex has a value if the previous position was valid */
    m_mappedIndex.value() += m_strides[dimension] * delta;
    assert(m_mappedIndex < m_bufferSize);
  }else if(new_position_is_inside_content){ /* if the new position is inside bounds, then the mapped index can be caluclated */
    m_mappedIndex = calculate_mapped_position(m_position);
  }else m_mappedIndex = {}; /* No mapped index for positions inside the padding */
  return *this;
}

std::optional<std::uint32_t> NDArrayIndex::calculate_mapped_position(const std::vector<std::uint32_t>& position) const{
  assert(position.size() == m_strides.size());
  if(!inside_content(position))
    return {};

  std::uint32_t result_index = 0u;
  for(std::uint32_t dim = 0; dim < position.size(); ++dim){
    result_index += (position[dim] - std::max(m_padding, -m_padding)) * m_strides[dim];
  }
  return result_index;
}

bool NDArrayIndex::inside_bounds(const std::vector<std::uint32_t>& position, std::uint32_t dimension, std::int32_t delta) const{
  std::uint32_t dimension_index = 0;
  return std::all_of(position.begin(), position.end(), 
    [this, &dimension_index, dimension, delta](const std::uint32_t& pos){
      std::int32_t position = static_cast<std::int32_t>(pos);
      if(dimension_index == dimension) position += delta;
      return( (0 <= position)&&(position < static_cast<int32_t>(2 * std::max(0, m_padding) + m_dimensions[dimension_index++])) );
    }
  );
}

bool NDArrayIndex::inside_content(const std::vector<std::uint32_t>& position, std::uint32_t dimension, std::int32_t delta) const{
  std::uint32_t dimension_index = 0;
  return std::all_of(position.begin(), position.end(), 
    [this, &dimension_index, dimension, delta](const std::uint32_t& pos){
      std::int32_t actual_position = static_cast<std::int32_t>(pos);
      if(dimension_index == dimension) actual_position += delta;
      return( 
        (std::max(m_padding, -m_padding) <= actual_position)
        &&(actual_position < static_cast<std::int32_t>(m_dimensions[dimension_index++] + m_padding)) 
      );
    }
  );
}

std::vector<NDArrayIndex::IntervalPart> NDArrayIndex::mappable_parts_of(
  const std::vector<std::uint32_t>& position, std::uint32_t dimension, std::int32_t delta
) const{
  std::vector<NDArrayIndex::IntervalPart> result;
  bool part_in_progress = false;
  for(std::int32_t delta_index = 0; delta_index < delta; delta_index += std::copysign(1, delta)){
    const bool current_position_in_inside_content = inside_content(position, dimension, delta_index);
    if(current_position_in_inside_content && part_in_progress){
      assert(0 < result.size());
      ++std::get<1>(result.back()); /* Increase the size of the current part of the interval */
    }else if(current_position_in_inside_content){ /* If the interval iteration became inside bounds */
      result.push_back({(position[dimension] + delta_index), 1}); /* Add the new part as a result */
      part_in_progress = true;
    }else part_in_progress = false;
  }
  return result;
}

and with the following tests:

TEST_CASE("Testing NDArray Indexing with a 2D array without padding", "[NDArray]"){
  std::uint32_t width = rand()%100;
  std::uint32_t height = rand()%100;
  NDArrayIndex idx({width, height});

  for(std::uint32_t variant = 0; variant < 5; ++variant){
    std::uint32_t x = rand()%width;
    std::uint32_t y = rand()%height;
    idx.set({x,y});
    REQUIRE(idx.inside_bounds());
    REQUIRE(idx.mapped_position().has_value());
    REQUIRE(idx.mapped_position().value() == (x + (y * width)));
    std::uint32_t elements_after_x_row = width - x;
    REQUIRE(1 == idx.mappable_parts_of(0,width).size());
    REQUIRE(x == std::get<0>(idx.mappable_parts_of(0,width)[0]));
    REQUIRE(elements_after_x_row == std::get<1>(idx.mappable_parts_of(0,width)[0]));
    /*!Note: using width in the above interfaces because it is guaranteed
     * that an interval of that size spans over the relevant dimension
     * */
  }

  REQUIRE(idx.buffer_size() == (width * height));
  idx.set({0,0});
  for(std::uint32_t i = 0; i < idx.buffer_size(); ++i){
    REQUIRE(idx.inside_bounds());
    REQUIRE(idx.inside_content());
    REQUIRE(idx.mapped_position().has_value() == true);
    REQUIRE(idx.mapped_position().value() == i);
    idx.step();
  }
}

TEST_CASE("Testing NDArray Indexing with a 2D array with positive padding", "[NDArray][padding]"){
  std::uint32_t width = 1 + rand()%20;
  std::uint32_t height = 1 + rand()%20;
  std::int32_t padding = 5;
  NDArrayIndex idx({width, height}, padding);

  for(std::uint32_t variant = 0; variant < 5; ++variant){
    std::uint32_t x = padding + rand()%(width);
    std::uint32_t y = padding + rand()%(height);
    idx.set({x,y});
    REQUIRE(idx.inside_bounds());
    REQUIRE(idx.mapped_position().has_value());
    REQUIRE( idx.mapped_position().value() == (x - padding + ((y - padding) * width)) );
    std::uint32_t elements_after_x_row = padding + width - x;
    REQUIRE(1 == idx.mappable_parts_of(0,width).size());
    REQUIRE(x == std::get<0>(idx.mappable_parts_of(0,width)[0]));
    REQUIRE(elements_after_x_row == std::get<1>(idx.mappable_parts_of(0,width)[0]));
  }

  REQUIRE(idx.buffer_size() == (width * height));
  std::uint32_t x = 0u;
  std::uint32_t y = 0u;
  std::uint32_t reference_mapped_position = 0u;
  idx.set({0,0});
  for(std::uint32_t i = 0; i < idx.buffer_size(); ++i){
    if(
      (padding <= static_cast<std::int32_t>(x) && x < (padding + width))
      &&(padding <= static_cast<std::int32_t>(y) && y < (padding + height))      
    ){
      REQUIRE(idx.inside_bounds());
      REQUIRE(idx.inside_content());
      REQUIRE(idx.mapped_position().has_value() == true);
      REQUIRE(idx.mapped_position().value() == reference_mapped_position);
      ++reference_mapped_position;
    }else{
      REQUIRE(idx.inside_bounds());
      REQUIRE(idx.mapped_position().has_value() == false);
    } 
    idx.step();
    if(x < padding + width + padding - 1){
      ++x;
    }else{
      x = 0;
      ++y;
    }
  }
}

TEST_CASE("Testing NDArray Indexing with a 2D array with negative padding", "[NDArray][padding]"){
  std::uint32_t width = 11 + rand()%20;
  std::uint32_t height = 11 + rand()%20;
  std::int32_t padding = -5;
  NDArrayIndex idx({width, height}, padding);

  for(std::uint32_t variant = 0; variant < 5; ++variant){
    std::uint32_t x = -padding + rand()%(width + 2 * padding);
    std::uint32_t y = -padding + rand()%(height + 2 * padding);
    idx.set({x,y});

    REQUIRE(idx.inside_bounds());
    REQUIRE(idx.mapped_position().has_value());
    REQUIRE( idx.mapped_position().value() == (x + padding + ((y + padding) * (width + 2 * padding))) );
    std::uint32_t elements_after_x_row = padding + width - x;
    REQUIRE(1 == idx.mappable_parts_of(0,width).size());
    REQUIRE(x == std::get<0>(idx.mappable_parts_of(0,width)[0]));
    REQUIRE(elements_after_x_row == std::get<1>(idx.mappable_parts_of(0,width)[0]));
  }

  REQUIRE(idx.buffer_size() == (width * height));
  std::uint32_t x = 0u;
  std::uint32_t y = 0u;
  std::uint32_t reference_mapped_position = 0u;
  idx.set({0,0});
  for(std::uint32_t i = 0; i < idx.buffer_size(); ++i){
    if(
      (-padding <= static_cast<std::int32_t>(x) && x < (padding + width))
      &&(-padding <= static_cast<std::int32_t>(y) && y < (padding + height))      
    ){
      REQUIRE(idx.inside_bounds());
      REQUIRE(idx.inside_content());
      REQUIRE(idx.mapped_position().has_value() == true);
      REQUIRE(idx.mapped_position().value() == reference_mapped_position);
      ++reference_mapped_position;
    }else{
      REQUIRE(idx.inside_bounds());
      REQUIRE(idx.mapped_position().has_value() == false);
    } 
    idx.step();
    if(x < (width - 1)){
      ++x;
    }else{
      x = 0;
      ++y;
    }
  }
}

Answer 1

Add Doxygen documentation

I'm having a hard time understanding what the purpose is of all those member functions. Adding Doxygen documentation for the class and all its members would be of great help.

Do you need a variable number of dimensions?

Usually when dealing with data, you already know its dimensionality. It would make more sense then to make NDArrayIndex be a template, with the template parameter being the number of dimensions, and then use std::array instead of std::vector to store the things that are now stored in std::vectors. This will allow the compiler to optimize the code much better, and avoids all the memory allocations.

Even if you don't know the number of dimensions up front, maybe you can make the template parameter be the maximum number of dimensions, so you can still use std::arrays. Consider that, at least on Linux on AMD64, the size of an empty std::vector is 24 bytes, which is the same size as a std::array<uint32_t, 6>.

Why is padding a scalar?

It's weird that m_dimensions and m_strides are vectors, but m_padding is a scalar. There is no reason why you could not have a different padding size for each dimension.

Prefer declaring a struct instead of using std::pair

While std::pair and std::tuple are sometimes helpful in generic code, if you can just declare a struct instead, prefer the latter. This allows you to give names to the two elements of the pair, and avoids the ugly calls to std::get<>():

struct IntervalPart {
    std::uint32_t start;
    std::uint32_t size;
};

Unsafe conversions between signed and unsigned integers

The API makes it seem like any unsigned 32-bit value is safe to be used, but the cast to std::int32_t inside inside_bounds() makes large values unsafe. You should handle this somehow.

Consider adding iterators

It looks to me like in the end, you want to iterate over the mappable part of a multi-dimensional array. Now you have to call a mix of mappable_parts_of(), step() and mapped_position(). But all of this is slow; if you have the mappable parts, you shouldn't need the fancyness of step(), and you don't need the bounds checking done by mapped_position(). I think that an interface that provides an iterator to efficiently iterate over a mappable part would be best. For example:

NDArray array = ...;
NDArrayIndex index = ...;

for (auto position: index.mappable_parts_range(...)) {
    do_something_with(array[position]);
}

Of course, for convolutions you need two positions, one for each array you want to convolve, that move together, so you either want to be able to derive one from the other quickly, or have an even more fancy API that allows provides you with both positions simultaneously, or perhaps even pass in the arrays directly and have it return references to array elements:

NDArray array1, array2;
float sum = {};
...
for (auto& [el1, el2]: overlapping_range(array1, array2, ...)) {
    sum += el1 * el2;
}

Answer 2

Do you really need exactly 32 bits where you use std::uint32_t? This prevents your code compiling where a 32-bit type is not available. Perhaps std::uint_fast32_t would be a better choice?