TensorIndices.h

Go to the documentation of this file.

#ifndef LARDATA_UTILITIES_TENSORINDICES_H
#define LARDATA_UTILITIES_TENSORINDICES_H

// C/C++ standard libraries
#include <cstddef>     // std::size_t
#include <stdexcept>   // std::out_of_range
#include <string>      // std::to_string()
#include <type_traits> // std::enable_if

namespace util {

  struct TensorIndicesBasicTypes {

    using Index_t = std::ptrdiff_t;

    using DimSize_t = std::size_t;

    using LinIndex_t = std::size_t;

  }; // TensorIndicesBasicTypes

  // This is the class declaration; it has a specialisation for rank 1 and
  // a generic implementation (valid for rank larger than 1).
  template <unsigned int RANK>
  class TensorIndices;

  namespace details {
    template <unsigned int RANK, unsigned int DIM>
    struct ExtractTensorDimension;
  } // namespace details

  template <>
  class TensorIndices<1U> {
  public:
    using Index_t = TensorIndicesBasicTypes::Index_t;

    using DimSize_t = TensorIndicesBasicTypes::DimSize_t;

    using LinIndex_t = TensorIndicesBasicTypes::LinIndex_t;

    static constexpr unsigned int rank() { return 1U; }

    TensorIndices(DimSize_t dim) : dimSize(dim) {}

    template <
      typename ITER,
      typename = std::enable_if_t<std::is_convertible<decltype(*(ITER())), DimSize_t>::value, void>>
    TensorIndices(ITER dimIter) : TensorIndices(*dimIter)
    {}

    LinIndex_t operator()(Index_t index) const { return index; }

    template <typename ITER>
    std::enable_if_t<std::is_convertible<decltype(*(ITER())), DimSize_t>::value, LinIndex_t>
    operator()(ITER indexIter) const
    {
      return *indexIter;
    }

    LinIndex_t at(Index_t index) const { return checkOuterIndex(index); }

    template <typename ITER>
    std::enable_if_t<std::is_convertible<decltype(*(ITER())), DimSize_t>::value, LinIndex_t> at(
      ITER indexIter) const
    {
      return checkOuterIndex(*indexIter);
    }

    bool has(Index_t index) const { return (index >= 0) && ((DimSize_t)index < size()); }

    template <typename ITER>
    std::enable_if_t<std::is_convertible<decltype(*(ITER())), DimSize_t>::value, bool> has(
      ITER indexIter) const
    {
      return (*indexIter >= 0) && ((DimSize_t)*indexIter < size());
    }

    template <unsigned int DIM = 0>
    DimSize_t size() const
    {
      return (DIM == 0) ? totalSize() : dim<DIM>();
    }

    template <unsigned int DIM>
    DimSize_t dim() const
    {
      static_assert(DIM == 0, "Invalid dimension requested");
      return dimSize;
    } // dim()

    template <unsigned int DIM>
    bool hasIndex(Index_t index) const
    {
      static_assert(DIM == 0, "Invalid dimension requested");
      return has(index);
    }

    bool hasLinIndex(LinIndex_t linIndex) const { return has((Index_t)linIndex); }

    bool operator==(TensorIndices<1U> const& t) const { return t.size() == size(); }
    bool operator!=(TensorIndices<1U> const& t) const { return t.size() != size(); }

  protected:
    Index_t checkOuterIndex(Index_t index) const
    {
      if (has(index)) return index; // good
      throw std::out_of_range("Requested index " + std::to_string(index) +
                              " for a dimension of size " + std::to_string(size()));
    }

  protected:
    DimSize_t dim0() const { return dimSize; }

  private:
    template <unsigned int R, unsigned int D>
    friend struct details::ExtractTensorDimension;

    DimSize_t dimSize; 

    DimSize_t totalSize() const { return dim0(); }

  }; // class TensorIndices<1>

  template <unsigned int RANK>
  class TensorIndices : private TensorIndices<1U> {
    static_assert(RANK > 1, "TensorIndices must have rank 1 or higher");

    using Base_t = TensorIndices<1U>; 

  public:
    using Index_t = Base_t::Index_t;

    using DimSize_t = Base_t::DimSize_t;

    using LinIndex_t = Base_t::LinIndex_t;

    static constexpr unsigned int rank() { return RANK; }

    using MinorTensor_t = TensorIndices<rank() - 1>;

    template <typename... OTHERDIMS>
    TensorIndices(DimSize_t first, OTHERDIMS... others)
      : Base_t(first), m(others...), totSize(dim0() * m.size())
    {}

    // only if ITER can be dereferenced into something convertible to DimSize_t
    template <
      typename ITER,
      typename = std::enable_if_t<std::is_convertible<decltype(*(ITER())), DimSize_t>::value, void>>
    TensorIndices(ITER dimIter) : Base_t(dimIter), m(++dimIter), totSize(dim0() * m.size())
    {}

    template <typename... OTHERINDICES>
    LinIndex_t operator()(Index_t first, OTHERINDICES... others) const
    {
      return first * minorTensor().size() + minorTensor()(others...);
    }

    template <typename ITER>
    std::enable_if_t<std::is_convertible<decltype(*(ITER())), DimSize_t>::value, LinIndex_t>
    operator()(ITER indexIter) const
    {
      auto const baseSize = (*indexIter) * minorTensor().size();
      return baseSize + minorTensor()(++indexIter);
    }

    template <typename... OTHERINDICES>
    LinIndex_t at(Index_t first, OTHERINDICES... others) const
    {
      return Base_t::checkOuterIndex(first) * minorTensor().size() + minorTensor().at(others...);
    }

    template <typename ITER>
    std::enable_if_t<std::is_convertible<decltype(*(ITER())), DimSize_t>::value, LinIndex_t> at(
      ITER indexIter) const
    {
      auto const baseSize = Base_t::checkOuterIndex(*indexIter) * minorTensor().size();
      return baseSize + minorTensor()(++indexIter);
    }

    template <typename... OTHERINDICES>
    bool has(Index_t first, OTHERINDICES... others) const
    {
      return Base_t::has(first) && minorTensor().has(others...);
    }

    template <typename ITER>
    std::enable_if_t<std::is_convertible<decltype(*(ITER())), DimSize_t>::value, bool> has(
      ITER indexIter) const
    {
      return Base_t::has(*indexIter) ? minorTensor().has(++indexIter) : false;
    }

    template <unsigned int DIM>
    DimSize_t dim() const
    {
      return details::ExtractTensorDimension<rank(), DIM>::dim(*this);
    }

    template <unsigned int DIM>
    bool hasIndex(Index_t index) const
    {
      return (index >= 0U) && ((DimSize_t)index < dim<DIM>());
    }

    template <unsigned int DIM = 0>
    DimSize_t size() const
    {
      return (DIM == 0) ? totalSize() : details::ExtractTensorDimension<rank(), DIM>::size(*this);
    }

    bool hasLinIndex(LinIndex_t linIndex) const
    {
      return (linIndex >= 0U) && ((DimSize_t)linIndex < size());
    }

    MinorTensor_t const& minorTensor() const { return m; }

    bool operator==(TensorIndices<RANK> const& t) const
    {
      return Base_t::operator==(t) && (minorTensor() == t.minorTensor());
    }

    bool operator!=(TensorIndices<RANK> const& t) const
    {
      return Base_t::operator!=(t) || (minorTensor() != t.minorTensor());
    }

  protected:
    // need to be friend of the dimension extractor
    template <unsigned int R, unsigned int D>
    friend struct details::ExtractTensorDimension;

    MinorTensor_t m;   
    DimSize_t totSize; 

    DimSize_t totalSize() const { return totSize; }

  }; // class TensorIndices<>

  template <unsigned int RANK1,
            unsigned int RANK2,
            typename = std::enable_if_t<(RANK1 != RANK2), bool>>
  bool operator==(TensorIndices<RANK1> const& a, TensorIndices<RANK2> const& b)
  {
    return false;
  }

  template <unsigned int RANK1,
            unsigned int RANK2,
            typename = std::enable_if_t<(RANK1 != RANK2), bool>>
  bool operator!=(TensorIndices<RANK1> const& a, TensorIndices<RANK2> const& b)
  {
    return true;
  }

  template <typename... DIMS>
  auto makeTensorIndices(DIMS... dims)
  {
    return TensorIndices<sizeof...(DIMS)>{TensorIndicesBasicTypes::DimSize_t(dims)...};
  }

  using MatrixIndices = TensorIndices<2U>;

} // namespace util

//------------------------------------------------------------------------------
//--- details implementation
//---
namespace util {
  namespace details {

    template <unsigned int RANK, unsigned int DIM>
    struct ExtractTensorDimension {
      static_assert(RANK > DIM, "Invalid dimension requested");
      static TensorIndicesBasicTypes::DimSize_t dim(TensorIndices<RANK> const& t)
      {
        return ExtractTensorDimension<RANK - 1U, DIM - 1U>::dim(t.minorTensor());
      }

      static TensorIndicesBasicTypes::DimSize_t size(TensorIndices<RANK> const& t)
      {
        return ExtractTensorDimension<RANK - 1U, DIM - 1U>::size(t.minorTensor());
      }

    }; // ExtractTensorDimension<>()

    template <unsigned int RANK>
    struct ExtractTensorDimension<RANK, 0U> {
      static_assert(RANK > 0, "Invalid rank 0 for TensorIndices");
      static TensorIndicesBasicTypes::DimSize_t dim(TensorIndices<RANK> const& t)
      {
        return t.dim0();
      }

      static TensorIndicesBasicTypes::DimSize_t size(TensorIndices<RANK> const& t)
      {
        return t.size();
      }

    }; // ExtractTensorDimension<RANK, 0>()

  } // namespace details
} // namespace util

#endif // LARDATA_UTILITIES_TENSORINDICES_H