Cartesian.h

#ifndef DASH__CARTESIAN_H_
#define DASH__CARTESIAN_H_

#include <dash/Types.h>
#include <dash/Dimensional.h>
#include <dash/Exception.h>
#include <dash/internal/Logging.h>

#include <array>
#include <algorithm>
#include <sstream>
#include <iostream>
#include <cstring>
#include <type_traits>

namespace dash {

template<
  dim_t    NumDimensions,
  typename SizeType = dash::default_size_t >
class CartesianSpace
{
private:
  typedef typename std::make_signed<SizeType>::type
    IndexType;
  typedef CartesianSpace<NumDimensions, SizeType>
    self_t;

public:
  typedef IndexType                           index_type;
  typedef SizeType                            size_type;
  typedef std::array<SizeType, NumDimensions> extents_type;

public:
  template<dim_t NDim_, typename SizeType_>
  friend std::ostream & operator<<(
    std::ostream                           & os,
    const CartesianSpace<NDim_, SizeType_> & cartesian_space);

protected:
  SizeType     _size    = 0;
  SizeType     _rank    = NumDimensions;
  extents_type _extents = {{  }};

public:
  typedef std::integral_constant<dim_t, NumDimensions> ndim;

  constexpr CartesianSpace()
  : _size(0),
    _rank(NumDimensions)
  { }

  template <typename... Args>
  CartesianSpace(SizeType arg, Args... args)
  : _size(0),
    _rank(NumDimensions) {
    resize(arg, args...);
  }

  CartesianSpace(
    const extents_type & extents)
  : _size(0),
    _rank(NumDimensions) {
    resize(extents);
  }

  bool operator==(const self_t & other) const {
    if (this == &other) {
      return true;
    }
    for(auto i = 0; i < NumDimensions; i++) {
      if (_extents[i] != other._extents[i]) {
        return false;
      }
    }
    // No need to compare _size as it is derived from _extents.
    return true;
  }

  constexpr bool operator!=(const self_t & other) const {
    return !(*this == other);
  }

  template<typename... Args>
  void resize(SizeType arg, Args... args) {
    static_assert(
      sizeof...(Args) == NumDimensions-1,
      "Invalid number of arguments");
    std::array<SizeType, NumDimensions> extents =
      {{ arg, (SizeType)(args)... }};
    resize(extents);
  }

  template<typename SizeType_>
  void resize(const std::array<SizeType_, NumDimensions> & extents) {
    // Update size:
    _size = 1;
    for(auto i = 0; i < NumDimensions; i++ ) {
      _extents[i] = static_cast<SizeType>(extents[i]);
      _size      *= _extents[i];
    }
  }

  void resize(dim_t dim, SizeType extent) {
    _extents[dim] = extent;
    resize(_extents);
  }

  constexpr dim_t rank() const noexcept {
    return _rank;
  }

  constexpr SizeType size() const noexcept {
    return _size;
  }

  constexpr const extents_type & extents() const noexcept {
    return _extents;
  }

  SizeType extent(dim_t dim) const {
    DASH_ASSERT_RANGE(
      0, dim, NumDimensions-1,
      "Dimension for CartesianSpace::extent(dim) is out of bounds");
    return _extents[dim];
  }
}; // class CartesianSpace

template<
  dim_t    NumDimensions,
  typename SizeType = dash::default_size_t >
class SizeSpec : public CartesianSpace<NumDimensions, SizeType>
{
private:
  typedef CartesianSpace<NumDimensions, SizeType>
    parent_t;

public:
  constexpr SizeSpec() : parent_t() {
  }

  template <typename... Args>
  constexpr SizeSpec(SizeType arg, Args... args)
  : parent_t(arg, args...) {
  }

  constexpr SizeSpec(
    const ::std::array<SizeType, NumDimensions> & extents)
  : parent_t(extents) {
  }
};

template<
  dim_t      NumDimensions,
  MemArrange Arrangement    = ROW_MAJOR,
  typename   IndexType      = dash::default_index_t >
class CartesianIndexSpace
{
private:
  typedef typename std::make_unsigned<IndexType>::type
    SizeType;
  typedef CartesianIndexSpace<NumDimensions, Arrangement, IndexType>
    self_t;
  typedef ViewSpec<NumDimensions, IndexType>
    ViewSpec_t;

public:
  typedef IndexType                           index_type;
  typedef SizeType                            size_type;
  typedef std::array<SizeType, NumDimensions> extents_type;

  template<dim_t NDim_>
  friend std::ostream & operator<<(
    std::ostream                     & os,
    const CartesianIndexSpace<NDim_> & cartesian_space);

protected:
  SizeType     _size             = 0;
  SizeType     _rank             = NumDimensions;
  extents_type _extents          = {  };
  extents_type _offset_row_major = { };
  extents_type _offset_col_major = { };

public:
  typedef std::integral_constant<dim_t, NumDimensions> ndim;

  constexpr CartesianIndexSpace() = default;

  CartesianIndexSpace(
    const extents_type & extents)
  : _extents(extents)
  {
    resize(extents);
  }

  template<typename... Args>
  CartesianIndexSpace(SizeType arg, Args... args)
  : _extents({{ }})
  {
    resize(arg, args...);
  }

  constexpr dim_t rank() const noexcept {
    return _rank;
  }

  bool operator==(const self_t & other) const {
    if (this == &other) {
      return true;
    }
    for(auto i = 0; i < NumDimensions; i++) {
      if (_extents[i] != other._extents[i]) {
        return false;
      }
      // Comparing either row- or column major offsets suffices:
      if (_offset_row_major[i] != other._offset_row_major[i]) {
        return false;
      }
    }
    // No need to compare _size as it is derived from _extents.
    return true;
  }

  constexpr bool operator!=(const self_t & other) const {
    return !(*this == other);
  }

  template<typename... Args>
  void resize(SizeType arg, Args... args) {
    static_assert(
      sizeof...(Args) == NumDimensions-1,
      "Invalid number of arguments");
    std::array<SizeType, NumDimensions> extents =
      {{ arg, (SizeType)(args)... }};
    resize(extents);
  }

  template<typename SizeType_>
  void resize(const std::array<SizeType_, NumDimensions> & extents) {
    // Update size:
    _size = 1;
    for(auto i = 0; i < NumDimensions; i++ ) {
      _extents[i] = static_cast<SizeType>(extents[i]);
      _size      *= _extents[i];
    }
    // Update offsets:
    _offset_row_major[NumDimensions-1] = 1;
    for(auto i = NumDimensions-2; i >= 0; --i) {
      _offset_row_major[i] = _offset_row_major[i+1] * _extents[i+1];
    }
    _offset_col_major[0] = 1;
    for(auto i = 1; i < NumDimensions; ++i) {
      _offset_col_major[i] = _offset_col_major[i-1] * _extents[i-1];
    }
  }

  void resize(dim_t dim, SizeType extent) {
    _extents[dim] = extent;
    resize(_extents);
  }

  constexpr SizeType size() const noexcept {
    return _size;
  }

  constexpr const extents_type & extents() const noexcept {
    return _extents;
  }

  SizeType extent(dim_t dim) const {
    DASH_ASSERT_RANGE(
      0, dim, NumDimensions-1,
      "Given dimension " << dim <<
      " for CartesianIndexSpace::extent(dim) is out of bounds");
    return _extents[dim];
  }

  template<
    typename... Args,
    MemArrange AtArrangement = Arrangement>
  constexpr IndexType at(
      IndexType arg, Args... args) const {
    static_assert(
      sizeof...(Args) == NumDimensions-1,
      "Invalid number of arguments");
    return at<AtArrangement>(
             std::array<IndexType, NumDimensions> {{
               arg, (IndexType)(args) ... }}
           );
  }

  template<
    MemArrange AtArrangement = Arrangement,
    typename OffsetType>
  IndexType at(
    const std::array<OffsetType, NumDimensions> & point) const {
    SizeType offs = 0;
    DASH_ASSERT_GT(_size, 0, "CartesianIndexSpace has size 0");
    for (auto i = 0; i < NumDimensions; i++) {
      DASH_ASSERT_RANGE(
        0,
        static_cast<IndexType>(point[i]),
        static_cast<IndexType>(_extents[i]-1),
        "Given coordinate for CartesianIndexSpace::at() exceeds extent");
      SizeType offset_dim = 0;
      if (AtArrangement == ROW_MAJOR) {
        offset_dim = _offset_row_major[i];
      } else if (AtArrangement == COL_MAJOR) {
        offset_dim = _offset_col_major[i];
      }
      offs += offset_dim * point[i];
    }
    return offs;
  }

  template<
    MemArrange AtArrangement = Arrangement,
    typename OffsetType>
  IndexType at(
    const std::array<OffsetType, NumDimensions> & point,
    const ViewSpec_t & viewspec) const {
    std::array<OffsetType, NumDimensions> coords{};
    for (auto d = 0; d < NumDimensions; ++d) {
      coords[d] = point[d] + viewspec.offset(d);
    }
    return at(coords);
  }

  template<MemArrange CoordArrangement = Arrangement>
  std::array<IndexType, NumDimensions> coords(
    IndexType index) const
  {
    DASH_ASSERT_GT(_size, 0, "CartesianIndexSpace has size 0");
    DASH_ASSERT_RANGE(
      0, index, static_cast<IndexType>(_size-1),
      "Given index for CartesianIndexSpace::coords() is out of bounds");

    ::std::array<IndexType, NumDimensions> pos{};
    if (CoordArrangement == ROW_MAJOR) {
      for(auto i = 0; i < NumDimensions; ++i) {
        pos[i] = index / _offset_row_major[i];
        index  = index % _offset_row_major[i];
      }
    } else if (CoordArrangement == COL_MAJOR) {
      for(auto i = NumDimensions-1; i >= 0; --i) {
        pos[i] = index / _offset_col_major[i];
        index  = index % _offset_col_major[i];
      }
    }
    return pos;
  }

  template<MemArrange CoordArrangement = Arrangement>
  std::array<IndexType, NumDimensions> coords(
    IndexType          index,
    const ViewSpec_t & viewspec) const
  {
    std::array<IndexType, NumDimensions> pos{};
    extents_type offset;
    if (CoordArrangement == ROW_MAJOR)
    {
      offset[NumDimensions-1] = 1;
      for(auto i = NumDimensions-2; i >= 0; --i)
        offset[i] = offset[i+1] * viewspec.extent(i+1);

      for(auto i = 0; i < NumDimensions; ++i)
      {
        pos[i] = index / offset[i] + viewspec.offset(i);
        index  = index % offset[i];
      }
    }
    else if (CoordArrangement == COL_MAJOR)
    {
      offset[0] = 1;
      for(auto i = 1; i < NumDimensions; ++i)
        offset[i] = offset[i-1] * viewspec.extent(i-1);

      for(auto i = NumDimensions-1; i >= 0; --i)
      {
        pos[i] = index / offset[i] + viewspec.offset(i);
        index  = index % offset[i];
      }
    }
    return pos;
  }

#if 0

  template<MemArrange CoordArrangement = Arrangement>
  bool includes_index(
    IndexType index,
    dim_t dimension,
    IndexType dim_offset) const {
    if (_rank == 1) {
      // Shortcut for trivial case
      return (index >= 0 && index < this->size());
    }
    auto base_offset = 0;
    if (CoordArrangement == COL_MAJOR) {
      base_offset = _offset_col_major[dimension];
    } else if (CoordArrangement == ROW_MAJOR) {
      base_offset = _offset_row_major[dimension];
    }
    for (auto d = 0; d < NumDimensions; ++d) {
      // TODO
    }
    return true;
  }
#endif

  template<dim_t U = NumDimensions>
  constexpr typename std::enable_if< (U > 0), SizeType >::type
  x(SizeType offs) const {
    return coords(offs)[0];
  }

  template<dim_t U = NumDimensions>
  constexpr typename std::enable_if< (U > 1), SizeType >::type
  y(SizeType offs) const {
    return coords(offs)[1];
  }

  template<dim_t U = NumDimensions>
  constexpr typename std::enable_if< (U > 2), SizeType >::type
  z(SizeType offs) const {
    return coords(offs)[2];
  }

}; // class CartesianIndexSpace

template<
  size_t     NumDimensions,
  MemArrange Arrangement   = ROW_MAJOR,
  typename   IndexType     = dash::default_index_t >
class LocalMemoryLayout :
  public CartesianIndexSpace<NumDimensions, Arrangement, IndexType>
{
private:
  typedef typename std::make_unsigned<IndexType>::type
    SizeType;
  typedef LocalMemoryLayout<NumDimensions, Arrangement, IndexType>
    self_t;
  typedef CartesianIndexSpace<NumDimensions, Arrangement, IndexType>
    parent_t;
  typedef ViewSpec<NumDimensions, IndexType>
    ViewSpec_t;
public:
  constexpr LocalMemoryLayout(
    const SizeSpec<NumDimensions> & sizespec,
    const DistributionSpec<NumDimensions> & distspec)
  : parent_t(sizespec),
    _distspec(distspec) {
  }

  constexpr LocalMemoryLayout(
    const DistributionSpec<NumDimensions> & distspec)
  : parent_t(SizeSpec<NumDimensions>()),
    _distspec(distspec) {
  }

  bool operator==(const self_t & other) const {
    if (!(parent_t::operator==(other))) {
      return false;
    }
    return _distspec == other._distspec;
  }

  constexpr bool operator!=(const self_t & other) const {
    return !(*this == other);
  }

  template<typename... Args>
  void resize(SizeType arg, Args... args) {
    static_assert(
      sizeof...(Args) == NumDimensions-1,
      "Invalid number of arguments");
    std::array<SizeType, NumDimensions> extents =
      {{ arg, (SizeType)(args)... }};
    resize(extents);
  }

  template<typename SizeType_>
  void resize(std::array<SizeType_, NumDimensions> extents) {
    if (!_distspec.is_tiled()) {
      parent_t::resize(extents);
    }
    // Tiles in at least one dimension
    // TODO
    DASH_THROW(dash::exception::NotImplemented,
               "CartesianIndexSpace::resize(IndexType) not implemented for tiles!");
  }

  template<
    typename... Args,
    MemArrange AtArrangement = Arrangement>
  constexpr IndexType at(
      IndexType arg, Args... args) const {
    static_assert(
      sizeof...(Args) == NumDimensions-1,
      "Invalid number of arguments");
    return at<AtArrangement>(
             std::array<IndexType, NumDimensions> {
               arg, (IndexType)(args) ... }
           );
  }

  template<
    MemArrange AtArrangement = Arrangement,
    typename OffsetType>
  IndexType at(
    const std::array<OffsetType, NumDimensions> & point) const {
    if (!_distspec.is_tiled()) {
      // Default case, no tiles
      return parent_t::at(point);
    }
    // Tiles in at least one dimension
    // TODO
    DASH_THROW(dash::exception::NotImplemented,
               "CartesianIndexSpace::at(IndexType) not implemented for tiles!");
  }

  template<
    MemArrange AtArrangement = Arrangement,
    typename OffsetType>
  IndexType at(
    const std::array<OffsetType, NumDimensions> & point,
    const ViewSpec_t & viewspec) const {
    std::array<OffsetType, NumDimensions> coords;
    for (auto d = 0; d < NumDimensions; ++d) {
      coords[d] = point[d] + viewspec[d].offset;
    }
    return at(coords);
  }

  template<MemArrange CoordArrangement = Arrangement>
  std::array<IndexType, NumDimensions> coords(IndexType index) const {
    if (!_distspec.is_tiled()) {
      // Default case, no tiles
      return parent_t::coords(index);
    }
    // Tiles in at least one dimension
    // TODO
    DASH_THROW(dash::exception::NotImplemented,
               "CartesianIndexSpace::coords(IndexType) not implemented for tiles!");
  }

private:
  DistributionSpec<NumDimensions> _distspec;
}; // class LocalMemoryLayout

template <
  dash::dim_t NumDimensions,
  typename    SizeType >
std::ostream & operator<<(
  std::ostream & os,
  const dash::CartesianSpace<NumDimensions, SizeType> & cartesian_space)
{
  std::ostringstream ss;
  ss << dash::typestr(cartesian_space)
     << ": "
     << "extents(";
  for (auto dim = 0; dim < NumDimensions; ++dim) {
    if (dim > 0) {
      ss << ",";
    }
    ss << cartesian_space.extents()[dim];
  }
  ss << ")";
  return operator<<(os, ss.str());
}

template <
  dash::dim_t NumDimensions >
std::ostream & operator<<(
  std::ostream & os,
  const dash::CartesianIndexSpace<NumDimensions> & cartesian_space)
{
  std::ostringstream ss;
  ss << dash::typestr(cartesian_space)
     << ": "
     << "extents(";
  for (auto dim = 0; dim < NumDimensions; ++dim) {
    if (dim > 0) {
      ss << ",";
    }
    ss << cartesian_space.extents()[dim];
  }
  ss << ")";
  return operator<<(os, ss.str());
}

}  // namespace dash

#endif // DASH__CARTESIAN_H_