range.h

Go to the documentation of this file.

/*
 *  This file is a part of TiledArray.
 *  Copyright (C) 2013  Virginia Tech
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef TILEDARRAY_RANGE_H__INCLUDED
#define TILEDARRAY_RANGE_H__INCLUDED

#include <TiledArray/range_iterator.h>
#include <TiledArray/permutation.h>
#include <TiledArray/size_array.h>

namespace TiledArray {


  class Range {
  public:
    typedef Range Range_; 
    typedef std::size_t size_type; 
    typedef std::vector<size_type> index; 
    typedef index index_type; 
    typedef detail::SizeArray<const size_type> size_array; 
    typedef size_array extent_type;    
    typedef std::size_t ordinal_type; 
    typedef detail::RangeIterator<size_type, Range_> const_iterator; 
    friend class detail::RangeIterator<size_type, Range_>;

  protected:

    size_type* data_ = nullptr;
    size_type offset_ = 0ul; 
    size_type volume_ = 0ul; 
    unsigned int rank_ = 0u; 

  private:


    template <typename Index>
    void init_range_data(const Index& lower_bound, const Index& upper_bound) {
      // Construct temp pointers
      size_type* MADNESS_RESTRICT const lower  = data_;
      size_type* MADNESS_RESTRICT const upper  = lower + rank_;
      size_type* MADNESS_RESTRICT const extent = upper + rank_;
      size_type* MADNESS_RESTRICT const stride = extent + rank_;
      const auto* MADNESS_RESTRICT const lower_data = detail::data(lower_bound);
      const auto* MADNESS_RESTRICT const upper_data = detail::data(upper_bound);

      // Set the volume seed
      volume_ = 1ul;
      offset_ = 0ul;

      // Compute range data
      for(int i = int(rank_) - 1; i >= 0; --i) {
        // Check input dimensions
        TA_ASSERT(lower_data[i] >= 0);
        TA_ASSERT(lower_data[i] < upper_data[i]);

        // Compute data for element i of lower, upper, and extent
        const size_type lower_bound_i = lower_data[i];
        const size_type upper_bound_i = upper_data[i];
        const size_type extent_i      = upper_bound_i - lower_bound_i;

        lower[i]  = lower_bound_i;
        upper[i]  = upper_bound_i;
        extent[i] = extent_i;
        stride[i] = volume_;
        offset_  += lower_bound_i * volume_;
        volume_  *= extent_i;
      }
    }


    template <typename Index,
              typename std::enable_if<detail::is_pair<
                  typename Index::value_type>::value>::type* = nullptr>
    void init_range_data(const Index& bound) {
      // Construct temp pointers
      size_type* MADNESS_RESTRICT const lower  = data_;
      size_type* MADNESS_RESTRICT const upper  = lower + rank_;
      size_type* MADNESS_RESTRICT const extent = upper + rank_;
      size_type* MADNESS_RESTRICT const stride = extent + rank_;
      const auto* MADNESS_RESTRICT const bound_data = detail::data(bound);

      // Set the volume seed
      volume_ = 1ul;
      offset_ = 0ul;

      // Compute range data
      for(int i = int(rank_) - 1; i >= 0; --i) {
        // Compute data for element i of lower, upper, and extent
        const size_type lower_bound_i = bound_data[i].first;
        const size_type upper_bound_i = bound_data[i].second;
        const size_type extent_i      = upper_bound_i - lower_bound_i;

        // Check input dimensions
        TA_ASSERT(lower_bound_i >= 0ul);
        TA_ASSERT(lower_bound_i < upper_bound_i);

        lower[i]  = lower_bound_i;
        upper[i]  = upper_bound_i;
        extent[i] = extent_i;
        stride[i] = volume_;
        offset_  += lower_bound_i * volume_;
        volume_  *= extent_i;
      }
    }


    template <typename Index,
              typename std::enable_if<detail::is_integral_list<
                  typename Index::value_type>::value>::type* = nullptr>
    void init_range_data(const Index& extents) {
      // Construct temp pointers
      size_type* MADNESS_RESTRICT const lower  = data_;
      size_type* MADNESS_RESTRICT const upper  = lower + rank_;
      size_type* MADNESS_RESTRICT const extent = upper + rank_;
      size_type* MADNESS_RESTRICT const stride = extent + rank_;
      const auto* MADNESS_RESTRICT const extent_data = detail::data(extents);

      // Set the offset and volume initial values
      volume_ = 1ul;
      offset_ = 0ul;

      // Compute range data
      for(int i = int(rank_) - 1; i >= 0; --i) {
        // Check bounds of the input extent
        TA_ASSERT(extent_data[i] > 0);

        // Get extent i
        const size_type extent_i = extent_data[i];

        lower[i]  = 0ul;
        upper[i]  = extent_i;
        extent[i] = extent_i;
        stride[i] = volume_;
        volume_  *= extent_i;
      }
    }


    template <typename ... Indices,
              typename std::enable_if<detail::is_integral_list<Indices...>::value>::type* = nullptr>
    void init_range_data(const std::tuple<Indices...>& extents) {
      const constexpr std::size_t rank = std::tuple_size<std::tuple<Indices...>>::value;
      TA_ASSERT(rank_ == rank);

      // Set the offset and volume initial values
      volume_ = 1ul;
      offset_ = 0ul;

      // initialize by recursion
      init_range_data_helper(extents, std::make_index_sequence<rank>{});
    }

    template <typename ... Indices, std::size_t ... Is>
    void init_range_data_helper(const std::tuple<Indices...>& extents, std::index_sequence<Is...>) {
      int workers[] = {0, (init_range_data_helper_iter<Is>(extents),0)...};
      ++workers[0];
    }

    template <std::size_t I, typename ... Indices>
    void init_range_data_helper_iter(const std::tuple<Indices...>& extents) {
      // Check bounds of the input extent
      TA_ASSERT(std::get<I>(extents) > 0ul);

      // Get extent i
      const size_type extent_i = std::get<I>(extents);

      size_type* MADNESS_RESTRICT const lower  = data_;
      size_type* MADNESS_RESTRICT const upper  = lower + rank_;
      size_type* MADNESS_RESTRICT const extent = upper + rank_;
      size_type* MADNESS_RESTRICT const stride = extent + rank_;

      lower[I]  = 0ul;
      upper[I]  = extent_i;
      extent[I] = extent_i;
      stride[I] = volume_;
      volume_  *= extent_i;
    }


    void init_range_data(const Permutation& perm,
        const size_type* MADNESS_RESTRICT const other_lower_bound,
        const size_type* MADNESS_RESTRICT const other_upper_bound)
    {
      // Create temporary pointers to this range data
      auto* MADNESS_RESTRICT const lower  = data_;
      auto* MADNESS_RESTRICT const upper  = lower + rank_;
      auto* MADNESS_RESTRICT const extent = upper + rank_;
      auto* MADNESS_RESTRICT const stride = extent + rank_;

      // Copy the permuted lower, upper, and extent into this range.
      for(unsigned int i = 0u; i < rank_; ++i) {
        const auto perm_i = perm[i];

        // Get the lower bound, upper bound, and extent from other for rank i.
        const auto other_lower_bound_i = other_lower_bound[i];
        const auto other_upper_bound_i = other_upper_bound[i];
        const auto other_extent_i = other_upper_bound_i - other_lower_bound_i;

        // Store the permuted lower bound, upper bound, and extent
        lower[perm_i]  = other_lower_bound_i;
        upper[perm_i]  = other_upper_bound_i;
        extent[perm_i] = other_extent_i;
      }

      // Recompute stride, offset, and volume
      volume_ = 1ul;
      offset_ = 0ul;
      for(int i = int(rank_) - 1; i >= 0; --i) {
        const auto lower_i = lower[i];
        const auto extent_i = extent[i];
        stride[i] = volume_;
        offset_ += lower_i * volume_;
        volume_ *= extent_i;
      }
    }

  public:


    Range() { }


    template <typename Index,
        typename std::enable_if<! std::is_integral<Index>::value>::type* = nullptr>
    Range(const Index& lower_bound, const Index& upper_bound) {
      const size_type n = detail::size(lower_bound);
      TA_ASSERT(n == detail::size(upper_bound));
      if(n) {
        // Initialize array memory
        data_ = new size_type[n << 2];
        rank_ = n;
        init_range_data(lower_bound, upper_bound);
      }
    }


    template <typename Index1,
              typename std::enable_if<std::is_integral<Index1>::value>::type* = nullptr>
    Range(const std::initializer_list<Index1>& lower_bound,
        const std::initializer_list<Index1>& upper_bound)
    {
      const size_type n = detail::size(lower_bound);
      TA_ASSERT(n == detail::size(upper_bound));
      if(n) {
        // Initialize array memory
        data_ = new size_type[n << 2];
        rank_ = n;
        init_range_data(lower_bound, upper_bound);
      }
    }


    template <typename Index,
        typename std::enable_if<! std::is_integral<Index>::value &&
        std::is_integral<typename Index::value_type>::value>::type* = nullptr>
    explicit Range(const Index& extent) {
      const size_type n = detail::size(extent);
      if(n) {
        // Initialize array memory
        data_ = new size_type[n << 2];
        rank_ = n;
        init_range_data(extent);
      }
    }


    template <typename Index1,
              typename std::enable_if<std::is_integral<Index1>::value>::type* = nullptr>
    explicit Range(const std::initializer_list<Index1>& extent) {
      const size_type n = detail::size(extent);
      if(n) {
        // Initialize array memory
        data_ = new size_type[n << 2];
        rank_ = n;
        init_range_data(extent);
      }
    }


    template <typename Index,
        typename std::enable_if<! std::is_integral<Index>::value &&
        detail::is_pair<typename Index::value_type>::value>::type* = nullptr>
    Range(const Index& bounds) {
      const size_type n = detail::size(bounds);
      if(n) {
        // Initialize array memory
        data_ = new size_type[n << 2];
        rank_ = n;
        init_range_data(bounds);
      }
    }


    template <typename Index1, typename Index2>
    Range(const std::initializer_list<std::pair<Index1,Index2>>& bounds) {
      const size_type n = detail::size(bounds);
      if(n) {
        // Initialize array memory
        data_ = new size_type[n << 2];
        rank_ = n;
        init_range_data(bounds);
      }
    }


    template<typename... Index,
        typename std::enable_if<detail::is_integral_list<Index...>::value>::type* = nullptr>
    explicit Range(const Index... extents) :
      Range(std::array<size_t, sizeof...(Index)>{{static_cast<std::size_t>(extents)...}})
    { }


    template <typename ... IndexPairs,
              typename std::enable_if<detail::is_integral_pair_list<IndexPairs...>::value>::type* = nullptr
             >
    explicit Range(const IndexPairs... bounds) :
    Range(std::array<std::pair<std::size_t,std::size_t>, sizeof...(IndexPairs)>{{static_cast<std::pair<std::size_t,std::size_t>>(bounds)...}})
    { }



    Range(const Range_& other) {
      if(other.rank_ > 0ul) {
        data_ = new size_type[other.rank_ << 2];
        offset_ = other.offset_;
        volume_ = other.volume_;
        rank_ = other.rank_;
        memcpy(data_, other.data_, (sizeof(size_type) << 2) * other.rank_);
      }
    }


    Range(Range_&& other) :
      data_(other.data_), offset_(other.offset_), volume_(other.volume_),
      rank_(other.rank_)
    {
      other.data_ = nullptr;
      other.offset_ = 0ul;
      other.volume_ = 0ul;
      other.rank_ = 0u;
    }


    Range(const Permutation& perm, const Range_& other) {
      TA_ASSERT(perm.dim() == other.rank_);

      if(other.rank_ > 0ul) {
        data_ = new size_type[other.rank_ << 2];
        rank_ = other.rank_;

        if(perm) {
          init_range_data(perm, other.data_, other.data_ + rank_);
        } else {
          // Simple copy will due.
          memcpy(data_, other.data_, (sizeof(size_type) << 2) * rank_);
          offset_ = other.offset_;
          volume_ = other.volume_;
        }
      }
    }

    ~Range() { delete [] data_; }


    Range_& operator=(const Range_& other) {
      if(rank_ != other.rank_) {
        delete [] data_;
        data_ = (other.rank_ > 0ul ? new size_type[other.rank_ << 2] : nullptr);
        rank_ = other.rank_;
      }
      memcpy(data_, other.data_, (sizeof(size_type) << 2) * rank_);
      offset_ = other.offset_;
      volume_ = other.volume_;

      return *this;
    }


    Range_& operator=(Range_&& other) {
      data_ = other.data_;
      offset_ = other.offset_;
      volume_ = other.volume_;
      rank_ = other.rank_;

      other.data_ = nullptr;
      other.offset_ = 0ul;
      other.volume_ = 0ul;
      other.rank_ = 0u;

      return *this;
    }


    unsigned int rank() const { return rank_; }


    const size_type* lobound_data() const { return data_; }


    size_array lobound() const { return size_array(lobound_data(), rank_); }


    size_type lobound(size_t dim) const {
      TA_ASSERT(dim < rank_);
      return *(lobound_data() + dim);
    }


    const size_type* upbound_data() const { return data_ + rank_; }


    size_array upbound() const {
      return size_array(upbound_data(), rank_);
    }


    size_type upbound(size_t dim) const {
      TA_ASSERT(dim < rank_);
      return *(upbound_data() + dim);
    }


    const size_type* extent_data() const { return data_ + (rank_ + rank_); }


    extent_type extent() const {
      return size_array(extent_data(), rank_);
    }


    size_type extent(size_t dim) const {
      TA_ASSERT(dim < rank_);
      return *(extent_data() + dim);
    }


    const size_type* stride_data() const { return data_ + (rank_ + rank_ + rank_); }


    size_array stride() const {
      return size_array(stride_data(), rank_);
    }


    size_type stride(size_t dim) const {
      TA_ASSERT(dim < rank_);
      return *(stride_data() + dim);
    }


    ordinal_type volume() const { return volume_; }

    ordinal_type area() const { return volume_; }


    ordinal_type offset() const { return offset_; }


    const_iterator begin() const { return const_iterator(data_, this); }


    const_iterator end() const { return const_iterator(data_ + rank_, this); }


    template <typename Index,
        typename std::enable_if<! std::is_integral<Index>::value, bool>::type* = nullptr>
    bool includes(const Index& index) const {
      TA_ASSERT(detail::size(index) == rank_);
      const size_type* MADNESS_RESTRICT const lower  = data_;
      const size_type* MADNESS_RESTRICT const upper = lower + rank_;

      bool result = (rank_ > 0u);
      auto it = std::begin(index); // TODO C++14 switch to std::cbegin
      for(unsigned int i = 0u; result && (i < rank_); ++i, ++it) {
        const size_type index_i = *it;
        const size_type lower_i = lower[i];
        const size_type upper_i = upper[i];
        result = result && (index_i >= lower_i) && (index_i < upper_i);
      }

      return result;
    }


    template <typename Integer>
    bool includes(const std::initializer_list<Integer>& index) const {
      return includes<std::initializer_list<Integer>>(index);
    }



    template <typename Ordinal>
    typename std::enable_if<std::is_integral<Ordinal>::value, bool>::type
    includes(Ordinal i) const {
      return include_ordinal_(i);
    }

    template <typename... Index>
    typename std::enable_if<(sizeof...(Index) > 1ul), size_type>::type
    includes(const Index&... index) const {
      const size_type i[sizeof...(Index)] = {static_cast<size_type>(index)...};
      return includes(i);
    }




    Range_& operator *=(const Permutation& perm);


    template <typename Index>
    Range_& resize(const Index& lower_bound, const Index& upper_bound) {
      const size_type n = detail::size(lower_bound);
      TA_ASSERT(n == detail::size(upper_bound));

      // Reallocate memory for range arrays
      if(rank_ != n) {
        delete [] data_;
        data_ = (n > 0ul ? new size_type[n << 2] : nullptr);
        rank_ = n;
      }
      if(n > 0ul)
        init_range_data(lower_bound, upper_bound);
      else
        volume_ = 0ul;

      return *this;
    }


    template <typename Index>
    Range_& inplace_shift(const Index& bound_shift) {
      TA_ASSERT(detail::size(bound_shift) == rank_);

      const auto* MADNESS_RESTRICT const bound_shift_data = detail::data(bound_shift);
      size_type* MADNESS_RESTRICT const lower = data_;
      size_type* MADNESS_RESTRICT const upper = data_ + rank_;
      const size_type* MADNESS_RESTRICT const stride = upper + rank_ + rank_;

      offset_ = 0ul;
      for(unsigned i = 0u; i < rank_; ++i) {
        // Load range data
        const auto bound_shift_i = bound_shift_data[i];
        auto lower_i = lower[i];
        auto upper_i = upper[i];
        const auto stride_i = stride[i];

        // Compute new range bounds
        lower_i += bound_shift_i;
        upper_i += bound_shift_i;

        // Update range data
        offset_ += lower_i * stride_i;
        lower[i] = lower_i;
        upper[i] = upper_i;
      }

      return *this;
    }


    template <typename Index>
    Range_ shift(const Index& bound_shift) {
      Range_ result(*this);
      result.inplace_shift(bound_shift);
      return result;
    }


    ordinal_type ordinal(const ordinal_type index) const {
      TA_ASSERT(includes(index));
      return index;
    }


    template <typename Index,
        typename std::enable_if<! std::is_integral<Index>::value>::type* = nullptr>
    ordinal_type ordinal(const Index& index) const {
      TA_ASSERT(detail::size(index) == rank_);
      TA_ASSERT(includes(index));

      size_type* MADNESS_RESTRICT const stride = data_ + rank_ + rank_ + rank_;

      size_type result = 0ul;
      auto index_it = std::begin(index);
      for(unsigned int i = 0u; i < rank_; ++i, ++index_it) {
        const size_type stride_i = stride[i];
        result += *(index_it) * stride_i;
      }

      return result - offset_;
    }


    template <typename... Index,
        typename std::enable_if<(sizeof...(Index) > 1ul)>::type* = nullptr>
    size_type ordinal(const Index&... index) const {
      const size_type temp_index[sizeof...(Index)] = { static_cast<size_type>(index)... };
      return ordinal(temp_index);
    }


    index idx(size_type index) const {
      // Check that index is contained by range.
      TA_ASSERT(includes(index));

      // Construct result coordinate index object and allocate its memory.
      Range_::index result(rank_, 0);

      // Get pointers to the data
      size_type * MADNESS_RESTRICT const result_data = result.data();
      size_type const * MADNESS_RESTRICT const lower = data_;
      size_type const * MADNESS_RESTRICT const size = data_ + rank_ + rank_;

      // Compute the coordinate index of index in range.
      for(int i = int(rank_) - 1; i >= 0; --i) {
        const size_type lower_i = lower[i];
        const size_type size_i = size[i];

        // Compute result index element i
        const size_type result_i = (index % size_i) + lower_i;
        index /= size_i;

        // Store result
        result_data[i] = result_i;
      }

      return result;
    }


    template <typename Index,
        typename std::enable_if<! std::is_integral<Index>::value>::type* = nullptr>
    const Index& idx(const Index& i) const {
      TA_ASSERT(includes(i));
      return i;
    }

    template <typename Archive,
        typename std::enable_if<madness::archive::is_input_archive<Archive>::value>::type* = nullptr>
    void serialize(const Archive& ar) {
      // Get rank
      unsigned int rank = 0ul;
      ar & rank;

      // Reallocate the array
      const unsigned int four_x_rank = rank << 2;
      if(rank_ != rank) {
        delete [] data_;
        data_ = (rank > 0u ? new size_type[four_x_rank] : nullptr);
        rank_ = rank;
      }

      // Get range data
      ar & madness::archive::wrap(data_, four_x_rank) & offset_ & volume_;
    }

    template <typename Archive,
        typename std::enable_if<madness::archive::is_output_archive<Archive>::value>::type* = nullptr>
    void serialize(const Archive& ar) const {
      ar & rank_ & madness::archive::wrap(data_, rank_ << 2) & offset_ & volume_;
    }

    void swap(Range_& other) {
      // Get temp data
      std::swap(data_, other.data_);
      std::swap(offset_, other.offset_);
      std::swap(volume_, other.volume_);
      std::swap(rank_, other.rank_);
    }

  private:


    template <typename Index>
    typename std::enable_if<std::is_signed<Index>::value, bool>::type
    include_ordinal_(Index i) const { return (i >= Index(0)) && (i < Index(volume_)); }


    template <typename Index>
    typename std::enable_if<! std::is_signed<Index>::value, bool>::type
    include_ordinal_(Index i) const { return i < volume_; }


    void increment(index& i) const {
      TA_ASSERT(includes(i));

      size_type const * MADNESS_RESTRICT const lower = data_;
      size_type const * MADNESS_RESTRICT const upper = data_ + rank_;

      for(int d = int(rank_) - 1; d >= 0; --d) {
        // increment coordinate
        ++i[d];

        // break if done
        if(i[d] < upper[d])
          return;

        // Reset current index to lower bound.
        i[d] = lower[d];
      }

      // if the current location was set to lower then it was at the end and
      // needs to be reset to equal upper.
      std::copy(upper, upper + rank_, i.begin());
    }


    void advance(index& i, std::ptrdiff_t n) const {
      TA_ASSERT(includes(i));
      const size_type o = ordinal(i) + n;
      TA_ASSERT(includes(o));
      i = idx(o);
    }


    std::ptrdiff_t distance_to(const index& first, const index& last) const {
      TA_ASSERT(includes(first));
      TA_ASSERT(includes(last));
      return ordinal(last) - ordinal(first);
    }

  }; // class Range

  inline Range& Range::operator *=(const Permutation& perm) {
    TA_ASSERT(perm.dim() == rank_);
    if(rank_ > 1ul) {
      // Copy the lower and upper bound data into a temporary array
      size_type* MADNESS_RESTRICT const temp_lower = new size_type[rank_ << 1];
      const size_type* MADNESS_RESTRICT const temp_upper = temp_lower + rank_;
      std::memcpy(temp_lower, data_, (sizeof(size_type) << 1) * rank_);

      init_range_data(perm, temp_lower, temp_upper);

      // Cleanup old memory.
      delete[] temp_lower;
    }
    return *this;
  }

  inline void swap(Range& r0, Range& r1) { // no throw
    r0.swap(r1);
  }



  inline Range operator*(const Permutation& perm, const Range& r) {
    return Range(perm, r);
  }


  inline bool operator ==(const Range& r1, const Range& r2) {
    return (r1.rank() == r2.rank()) && !std::memcmp(r1.lobound_data(), r2.lobound_data(),
        r1.rank() * (2u * sizeof(Range::size_type)));
  }

  inline bool operator !=(const Range& r1, const Range& r2) {
    return (r1.rank() != r2.rank()) || std::memcmp(r1.lobound_data(), r2.lobound_data(),
        r1.rank() * (2u * sizeof(Range::size_type)));
  }


  inline std::ostream& operator<<(std::ostream& os, const Range& r) {
    os << "[ ";
    detail::print_array(os, r.lobound_data(), r.rank());
    os << ", ";
    detail::print_array(os, r.upbound_data(), r.rank());
    os << " )";
    return os;
  }

} // namespace TiledArray
#endif // TILEDARRAY_RANGE_H__INCLUDED