subArray1D.h

Go to the documentation of this file.

#ifndef _DLR_SUBARRAY1D_H_
#define _DLR_SUBARRAY1D_H_

#include <dlrNumeric/array1D.h>
#include <dlrNumeric/slice.h>

namespace dlr {

  namespace numeric {
    
    template <class Type>
    class SubArray1D {
    public:

      SubArray1D(const Array1D<Type>& source);

    
      SubArray1D(const Array1D<Type>& source, const Slice& slice0);

    
      SubArray1D(const SubArray1D<Type> &other);

    
      virtual
      ~SubArray1D() {}

    
      operator Array1D<Type>() const;

    
      inline size_t
      size() const {return this->m_size;}

    
      inline size_t
      start() const {return this->m_start;}


      inline size_t
      stride() const {return this->m_stride;}

    
      inline Array1D<Type>
      getArray() const {return this->m_source;}
  

      SubArray1D<Type>&
      operator=(const SubArray1D<Type>& other);

    
      SubArray1D<Type>&
      operator=(Type value);

      SubArray1D<Type>&
      operator+=(const SubArray1D<Type>& other);

      SubArray1D<Type>&
      operator-=(const SubArray1D<Type>& other);

    private:
      inline int
      abs(int argument) {return (argument >= 0) ? argument : -argument;}

      inline void
      checkArray1DSize(const Array1D<Type>& other) const;

      inline void
      checkSubArray1DSize(const SubArray1D<Type>& other) const;

      SubArray1D<Type>&
      copy(const SubArray1D<Type>& other);

      Array1D<Type> m_source;
      int m_start;
      int m_stop;
      int m_stride;
      size_t m_size;
    };


    /* Non-member functions */

    template <class Type>
    inline SubArray1D<Type>
    subArray(const Array1D<Type>& source) {return SubArray1D<Type>(source);}

  
    template <class Type>
    inline SubArray1D<Type>
    subArray(const Array1D<Type>& source,
             const Slice& rowSlice) {
      return SubArray1D<Type>(source, rowSlice);
    }

  
    template <class Type>
    std::ostream&
    operator<<(std::ostream& stream, const SubArray1D<Type>& subArray0);

  } // namespace numeric

} // namespace dlr


/* ======= Declarations to maintain compatibility with legacy code. ======= */

namespace dlr {

  using numeric::SubArray1D;
  using numeric::subArray;

} // namespace dlr


/*******************************************************************
 * Function definitions follow.  This would be a .C file
 * if SubArray1D weren't templated.
 *******************************************************************/

#include <cmath>
#include <dlrCommon/stridedPointer.h>

namespace dlr {

  namespace numeric {
    
    template <class Type>
    SubArray1D<Type>::
    SubArray1D(const Array1D<Type>& source)
      : m_source(source),
        m_start(0),
        m_stop(source.size()),
        m_stride(1),
        m_size(source.size())
    {
      // Empty
    }

    template <class Type>
    SubArray1D<Type>::
    SubArray1D(const Array1D<Type>& source, const Slice& slice0)
      : m_source(source),
        m_start(slice0.start()),
        m_stop(slice0.stop()),
        m_stride(slice0.stride()),
        m_size(0)
    {
      // It's convenient to be able to specify "last element" somehow.
      // Setting the stop value to zero will wrap to source.size()
      // if appropriate.
      if((this->m_stop == 0) && (this->m_stride > 0)) {
        this->m_stop = source.size();
      }
    
      // Negative indexing is also super-convenient
      while(this->m_start < 0) {
        this->m_start += source.size();
      }
      while(this->m_stop < 0) {
        this->m_stop += source.size();
      }

      // Now compute sizes (yuck).
      this->m_size = ((this->m_stop - this->m_start)
                      / this->m_stride); // integer division
      if(this->m_size < 0) {
        this->m_size = 0; 
      } else {
        // Can't think of a better way to do this.
        if(this->abs(this->m_size * this->m_stride)
           < this->abs(this->m_stop - this->m_start)) {
          ++this->m_size;
        }
      }

      // Make sure we won't read/write outside of the array.
      this->checkArray1DSize(source);
    }

  
    template <class Type>
    SubArray1D<Type>::
    SubArray1D(const SubArray1D<Type> &other)
      : m_source(other.m_source),
        m_start(other.m_start),
        m_stop(other.m_stop),
        m_stride(other.m_stride),
        m_size(other.m_size)
    {
      // Empty
    }

  
    template <class Type>
    SubArray1D<Type>::operator Array1D<Type>() const
    {
      Array1D<Type> returnVal(this->size());
      subArray(returnVal) = *this;
      return returnVal;
    }

  
    template <class Type>
    SubArray1D<Type>& SubArray1D<Type>::
    operator=(const SubArray1D<Type>& other)
    {
      this->checkSubArray1DSize(other);
      return this->copy(other);
    }

  
    template <class Type>
    SubArray1D<Type>& SubArray1D<Type>::
    operator=(Type value)
    {
      StridedPointer<Type> outPtr0(this->m_source.data(m_start), this->m_stride);
      StridedPointer<Type> outPtr1 = outPtr0 + this->size();
      std::fill(outPtr0, outPtr1, value);
      return *this;
    }

  
    template <class Type>
    SubArray1D<Type>& SubArray1D<Type>::
    operator+=(const SubArray1D<Type>& other)
    {
      this->checkSubArray1DSize(other);
      StridedPointer<Type> thisDataPtr0(this->m_source.data(this->m_start),
                                        this->m_stride);
      StridedPointer<Type> thisDataPtr1 = thisDataPtr0 + this->size();
      StridedPointer<const Type>
        otherDataPtr0(other.m_source.data(other.m_start), other.m_stride);
      std::transform(thisDataPtr0, thisDataPtr1, otherDataPtr0, thisDataPtr0,
                     std::plus<Type>());
      return *this;
    }

  
    template <class Type>
    SubArray1D<Type>& SubArray1D<Type>::
    operator-=(const SubArray1D<Type>& other)
    {
      this->checkSubArray1DSize(other);
      StridedPointer<Type> thisDataPtr0(this->m_source.data(this->m_start),
                                        this->m_stride);
      StridedPointer<Type> thisDataPtr1 = thisDataPtr0 + this->size();
      StridedPointer<const Type>
        otherDataPtr0(other.m_source.data(other.m_start), other.m_stride);
      std::transform(thisDataPtr0, thisDataPtr1, otherDataPtr0, thisDataPtr0,
                     std::minus<Type>());
      return *this;
    }
  
    template <class Type>
    inline void SubArray1D<Type>::
    checkArray1DSize(const Array1D<Type>& other) const
    {
#ifdef _DLRNUMERIC_CHECKBOUNDS_
      if((m_start < 0) || (m_start >= static_cast<int>(other.size()))) {
        std::ostringstream message;
        message << "Invalid start index: " << m_start << std::endl;
        DLR_THROW(IndexException, "SubArray1D::checkArray1DSize()",
                  message.str().c_str());
      }
      if(m_stride > 0) {
        if(m_stop > static_cast<int>(other.size())) {
          std::ostringstream message;
          message << "Invalid stop index: " << m_stop << std::endl;
          DLR_THROW(IndexException, "SubArray1D::checkArray1DSize()",
                    message.str().c_str());
        }
      } else if(m_stride < 0) {
        if(m_stop < -1) {
          std::ostringstream message;
          message << "Invalid stop index: " << m_stop << std::endl;
          DLR_THROW(IndexException, "SubArray1D::checkArray1DSize()",
                    message.str().c_str());
        }
      } else {
        // m_stride == 0
        DLR_THROW(IndexException, "SubArray1D::checkArray1DSize()",
                  "Invalid stride: 0");
      }
#endif
    }

  
    template <class Type>
    inline void SubArray1D<Type>::
    checkSubArray1DSize(const SubArray1D<Type>& other) const
    {
#ifdef _DLRNUMERIC_CHECKBOUNDS_
      if(other.size() != this->size()) {
        std::ostringstream message;
        message << "Size mismatch: " << other.size() << " vs. "
                << this->size() << std::endl;
        DLR_THROW(IndexException, "SubArray1D::checkSubArray1DSize()",
                  message.str().c_str());
      }
#endif
    }

  
    template <class Type>
    SubArray1D<Type>& SubArray1D<Type>::
    copy(const SubArray1D<Type>& other)
    {
      StridedPointer<const Type>
        inPtr0(other.m_source.data(other.m_start), other.stride());
      StridedPointer<const Type> inPtr1 = inPtr0 + other.size();
      StridedPointer<Type>
        outPtr0(this->m_source.data(this->m_start), this->m_stride);
      std::copy(inPtr0, inPtr1, outPtr0);
      return *this;
    }

    /* Non-member functions */

    template <class Type>
    std::ostream& operator<<(std::ostream& stream,
                             const SubArray1D<Type>& subArray)
    {
      Array1D<Type> array(subArray.size());
      subArray1D(array) = subArray;
      stream << "SubArray1D([";
      for(int index0 = 0; index < array.size() - 1; ++index0) {
        stream << array(index0) << ", ";
      }
      stream << array(array.size() - 1) << "])";
      stream.flush();
      return stream;
    }

  } // namespace numeric

} // namespace dlr

#endif // #ifdef _DLR_SUBARRAY1D_H_