IndexUtil.H

Go to the documentation of this file.

#ifndef __INDEX_UTIL_H__
#define __INDEX_UTIL_H__
#include <vector>
#include <istream>
#include <ostream>
#include <cstdlib>

namespace ix { namespace util {

    template<typename DataType>
    bool IntervalOverlaps(const DataType &i1, const DataType &i2,const DataType &j1, const DataType &j2)
    {
      return((i1 <= j2) && (i2 >= j1));
    }

    class sizeextent : public std::vector< std::pair<size_t,size_t> > 
    {      
      friend std::ostream &operator<<(std::ostream &outStream,const sizeextent &sizeExtent);
      friend std::istream &operator>>(std::istream &inStream,sizeextent &sizeExtent);
    private:
      size_t _nd;
      std::vector<size_t> _Np;
      std::vector<size_t> _N;
    public:
      sizeextent()
      {
        this->resize(0);
        _nd = 0;
        _Np.resize(0);
        _N.resize(0);
      }
      
      sizeextent(const sizeextent &inExtent)
      {
        if(!inExtent.empty())
          Copy(inExtent);
      }
      
      template<typename IntType>
      explicit sizeextent(const std::vector<IntType> &inStart,
                          const std::vector<IntType> &inSize)
      {
        this->resize(0);
        _nd = 0;
        _Np.resize(0);
        _N.resize(0);
        int nStart = inStart.size();
        int nSize = inSize.size();
        if(nStart != nSize)
          return;
        _nd = inStart.size();
        this->resize(_nd);
        for(size_t i = 0;i < _nd;i++){
          (*this)[i].first  = inStart[i];
          (*this)[i].second = inStart[i]+inSize[i]-1; 
        }
        Sync();
      };
      
      template<typename IntType>
      explicit sizeextent(const std::vector<std::vector<IntType> > &inextent)
      {
        _nd = inextent.size();
        this->resize(_nd); 
        for(size_t i = 0; i < _nd; i++){
          typename std::vector<IntType>::const_iterator iei = inextent[i].begin();
          while(iei != inextent[i].end()){
            (*this)[i].first  = *iei++;
            (*this)[i].second = *iei++;
          } 
        }     
        Sync();
      };
      
      template<typename ContainerType>
      explicit sizeextent(const ContainerType &inflatextent)
      {
        _nd = inflatextent.size()/2;
        this->resize(_nd);
        typename ContainerType::const_iterator ifIt = inflatextent.begin();
        for(int i = 0;i < _nd;i++){
          (*this)[i].first  = *ifIt++;
          (*this)[i].second = *ifIt++;
        }
        Sync();
      };
      
      // Constructor for flat extent [1:2*nd]
      template<typename T>
      explicit sizeextent(const T *src,int nd=3){
        _nd = nd;
        this->resize(_nd);
        unsigned int nindex = 0;
        typename sizeextent::iterator bsi = this->begin();
        while(bsi != this->end()){
          bsi->first  = src[nindex++];
          bsi->second = src[nindex++];
          bsi++;
        }
        Sync();
      };
      
      // Constructor for simple extent with sizes src [1:nd]
      template<typename T>
      explicit sizeextent(int nd,const T *src){
        _nd = nd;
        this->resize(_nd);
        unsigned int nindex = 0;
        typename sizeextent::iterator bsi = this->begin();
        while(bsi != this->end()){
          bsi->first  = 0;
          bsi->second = src[nindex++]-1;
          bsi++;
        }
        Sync();
      };

      void destroy()
        
      {
        this->resize(0);
        _nd = 0;
        _Np.resize(0);
        _N.resize(0);
      };
      
      template<typename ContainerType>
      void Init(const ContainerType &inflatextent){
        destroy();
        _nd = inflatextent.size()/2;
        this->resize(_nd);
        size_t nd = 0;
        typename sizeextent::iterator bsi = this->begin();
        typename ContainerType::const_iterator ifi = inflatextent.begin();
        while(bsi != this->end()){
          bsi->first  = *ifi++;
          bsi->second = *ifi++;
          bsi++;
        }
        Sync();
      };

      template<typename IntType>
      void Init(const std::vector<IntType> &inStart,
                const std::vector<IntType> &inSize)
      {
        destroy();
        _Np.resize(0);
        _N.resize(0);
        int nStart = inStart.size();
        int nSize = inSize.size();
        if(nStart != nSize)
          return;
        _nd = inStart.size();
        this->resize(_nd);
        for(size_t i = 0;i < _nd;i++){
          (*this)[i].first  = inStart[i];
          (*this)[i].second = inStart[i]+inSize[i]-1; 
        }
        Sync();
      };
      
      template<typename ContainerType>
      void InitSimple(const ContainerType &inSize){
        destroy();
        _nd = inSize.size();
        this->resize(_nd);
        typename sizeextent::iterator bsi = this->begin();
        typename ContainerType::const_iterator ifi = inSize.begin();
        while(bsi != this->end()){
          bsi->first  = 0;
          bsi->second = *ifi++-1;
          bsi++;
        }
        Sync();
      };
      

      template<typename T>
      void InitSimple(int numDim,const T *inSizes){
        destroy();
        _nd = numDim;
        this->resize(_nd);
        size_t nd = 0;
        typename sizeextent::iterator bsi = this->begin();
        while(bsi != this->end()){
          bsi->first  = 0;
          bsi->second = inSizes[nd++];
          bsi++;
        }
        Sync();
      };
      
      std::vector<size_t> Starts() const
      {
        std::vector<size_t> startArray;
        sizeextent::const_iterator myIt = this->begin();
        while(myIt != this->end())
          startArray.push_back(myIt++->first);
        return(startArray);
      };

      std::vector<size_t> Ends() const
      {
        std::vector<size_t> startArray;
        sizeextent::const_iterator myIt = this->begin();
        while(myIt != this->end())
          startArray.push_back(myIt++->second);
        return(startArray);        
      };

      std::vector<size_t> Sizes() const
      {
        std::vector<size_t> startArray;
        sizeextent::const_iterator myIt = this->begin();
        while(myIt != this->end()){
          startArray.push_back(myIt->second - myIt->first + 1);
          myIt++;
        }
        return(startArray);                
      };

      void Copy(const sizeextent &inExtent)
      {
        this->resize(inExtent.size());
        _nd = inExtent._nd;
        _Np = inExtent._Np;
        _N  = inExtent._N;
        sizeextent::iterator myIt = this->begin();
        sizeextent::const_iterator inIt = inExtent.begin();
        while(myIt != this->end())
          *myIt++ = *inIt++;
      };

      void Sync()
      {
        _nd = this->size();
        _N.resize(_nd);
        _Np.resize(_nd);
        for(size_t i = 0; i < _nd; i++){
          _N[i] = (*this)[i].second - (*this)[i].first + 1;
        }     
        _Np[0] = 1;
        for(size_t i = 1; i < _nd;i++){
          _Np[i] = _N[i-1] * _Np[i-1];
        }
      };
      
      size_t NNodes() const {
        if(this->empty())
          return(0);
        size_t nnodes = 1;
        for(size_t i = 0;i < _nd;i++)
          nnodes *= ((*this)[i].second - (*this)[i].first+1);
        return(nnodes);
      };

      size_t NNodes(int inDir) const {
        if(this->empty())
          return(0);
        if(inDir > _nd)
          return(0);
        return((*this)[inDir].second - (*this)[inDir].first+1);
      };
      
      size_t ND() const {return(this->size());};
      
      template<typename ContainerType>
      void Flatten(ContainerType &output) const
      {
        output.resize(0);
        for(size_t i = 0;i < _nd;i++){
          output.push_back((*this)[i].first);
          output.push_back((*this)[i].second);
        }
      };
      
      template<typename ContainerType>
      void dir_loop(size_t nd,size_t indoff,std::vector<size_t> &NP,
                    const sizeextent &inExtent, ContainerType &indices) const
      {
        size_t dindex = nd - 1;
        size_t start  = inExtent[dindex].first;
        size_t end    = inExtent[dindex].second;
        if(dindex == 0){
          for(size_t dd = start;dd <= end;dd++)
            indices.push_back(dd - (*this)[dindex].first+indoff);
        } else{
          for(size_t dd = start;dd <= end;dd++){          
            size_t ndoff = indoff + (dd - (*this)[dindex].first)*NP[dindex];
            dir_loop(dindex,ndoff,NP,inExtent,indices);
          }
        }
      };

      template<typename ContainerType>
      void GetFlatIndices(const sizeextent &extent,ContainerType &indices) const
      {
        int nd = ND();
        indices.resize(0);
        std::vector<size_t> N(nd);
        std::vector<size_t> NP(nd);
        for(int nc = 0;nc < nd;nc++){
          N[nc] = (*this)[nc].second - (*this)[nc].first + 1;
          NP[nc] = (nc == 0 ? 1 : N[nc-1] * NP[nc-1]);
        }
        NP[0] = 0;
        dir_loop(nd,0,NP,extent,indices);
      };
 
      template<typename ContainerType>
      ContainerType FlatIndices(const sizeextent &extent) const
      {
        ContainerType indices;
        GetFlatIndices<ContainerType>(extent,indices);
        return(indices);
      };

      void Overlap(const sizeextent &inextent,sizeextent &outextent) const
      {
        size_t nd = ND();
        outextent.resize(0);
        if(inextent.empty() || this->empty())
          return;
        bool match = true;
        for(size_t j = 0;j < nd && match;j++){
          if(!((((*this)[j].first >= inextent[j].first &&
                (*this)[j].first <= inextent[j].second) ||
               ((*this)[j].second >= inextent[j].first &&
                (*this)[j].second <= inextent[j].second)) || 
               ((inextent[j].first > (*this)[j].first) &&
                (inextent[j].second < (*this)[j].second))))
            match = false;
        }
        if(match){ // then this extent has some overlap with the searched
          outextent.resize(nd);
          for(size_t i = 0;i < nd;i++){
            outextent[i].first = std::max((*this)[i].first,inextent[i].first);
            outextent[i].second = std::min((*this)[i].second,inextent[i].second);
          }
          outextent.Sync();
        }
      };

     
      sizeextent Neighboring(const sizeextent &inextent,int nbrDim) const
      {
        sizeextent emptyExtent;
        const sizeextent &self(*this);
        if(inextent.empty() || self.empty())
          return(emptyExtent);
        size_t nd = ND();
        sizeextent nbrExtent(inextent);
        bool possible = true;
        int numNbrDims = 0;
        for(size_t j = 0;j < nd && possible;j++){
          if(IntervalOverlaps(self[j].first,self[j].second,nbrExtent[j].first,nbrExtent[j].second)){
            if(nbrExtent[j].first < self[j].first) nbrExtent[j].first = self[j].first;
            if(nbrExtent[j].second > self[j].second) nbrExtent[j].second = self[j].second;
          } else {
            if((nbrExtent[j].first != (self[j].second + 1)) &&
               (self[j].first != (nbrExtent[j].second + 1))){
              possible = false;
            } else {
              numNbrDims++;
            }
          }
        }
        if(possible && (numNbrDims == nbrDim))
          return(nbrExtent);
        return(emptyExtent);
      };

      sizeextent Overlap(const sizeextent &inextent) const
      {
        sizeextent outExtent;
        if(inextent.empty() || this->empty())
          return(outExtent);
        size_t nd = ND();
        bool match = true;
        for(size_t j = 0;j < nd && match;j++){
          if(!((((*this)[j].first >= inextent[j].first &&
                (*this)[j].first <= inextent[j].second) ||
               ((*this)[j].second >= inextent[j].first &&
                (*this)[j].second <= inextent[j].second)) || 
               ((inextent[j].first > (*this)[j].first) &&
                (inextent[j].second < (*this)[j].second))))
            match = false;
        }
        if(match){ // then this extent has some overlap with the searched
          outExtent.resize(nd);
          for(size_t i = 0;i < nd;i++){
            outExtent[i].first = std::max((*this)[i].first,inextent[i].first);
            outExtent[i].second = std::min((*this)[i].second,inextent[i].second);
          }
          outExtent.Sync();
        }
        return(outExtent);
      };

      int Protrude(const sizeextent &inExtent)
      {
        if(inExtent.size() != this->size())
          return(1);
        sizeextent::const_iterator inExtIt = inExtent.begin();
        sizeextent::iterator thisIt = this->begin();
        while(inExtIt != inExtent.end()){
          if(thisIt++->first < inExtIt++->first)
            return(1);
        }
        inExtIt = inExtent.begin();
        thisIt  = this->begin();
        while(inExtIt != inExtent.end()){
          thisIt->first -= inExtIt->first;
          thisIt->second += inExtIt->second;
          inExtIt++;
          thisIt++;
        }
        Sync();
        return(0);
      };

      template<typename ContainerType>
      int ProtrudeLimit(const ContainerType &inOffset,const ContainerType &inLimits)
      {
        typename ContainerType::const_iterator inLimIt = inLimits.begin();
        typename ContainerType::const_iterator inExtIt = inOffset.begin();
        sizeextent::iterator thisIt = this->begin();
        while(inExtIt != inOffset.end()){
          if(*inExtIt > thisIt->first)
            thisIt->first = 0;
          else
            thisIt->first -= *inExtIt;
          inExtIt++;
          thisIt->second = thisIt->second + *inExtIt++;
          if(thisIt->second >= *inLimIt)
            thisIt->second = *inLimIt-1;
          inLimIt++;
          thisIt++;
        }
        Sync();
        return(0);
      };

      sizeextent Protruded(const sizeextent &inExtent)
      {
        sizeextent returnExtent;
        if(inExtent.size() != this->size())
          return(returnExtent);
        sizeextent::const_iterator inExtIt = inExtent.begin();
        sizeextent::iterator thisIt = this->begin();
        while(inExtIt != inExtent.end()){
          if(thisIt++->first < inExtIt++->first)
            return(returnExtent);
        }
        returnExtent.Copy(*this);
        inExtIt = inExtent.begin();
        thisIt  = returnExtent.begin();
        while(inExtIt != inExtent.end()){
          thisIt->first -= inExtIt->first;
          thisIt->second += inExtIt->second;
          inExtIt++;
          thisIt++;
        }
        returnExtent.Sync();
        return(returnExtent);
      };

      template<typename ContainerType>
      int Protrude(const ContainerType &protrusionOffsets)
      {
        if(protrusionOffsets.size() != _nd*2)
          return(1);
        typename ContainerType::const_iterator offsetIt = protrusionOffsets.begin();
        sizeextent::iterator extentIt = this->begin();
        while(extentIt != this->end()){
          std::pair<size_t,size_t> &curPair(*extentIt++);
          if((*offsetIt < 0) && (-*offsetIt > curPair.first))
            return(1);
          offsetIt++;
          if((*offsetIt < 0) && (-*offsetIt > curPair.second))
            return(1);
          offsetIt++;
        }
        offsetIt = protrusionOffsets.begin();
        extentIt = this->begin();
        while(extentIt != this->end()){
          std::pair<size_t,size_t> &curPair(*extentIt++);
          curPair.first  += *offsetIt++;
          curPair.second += *offsetIt++;
        }
        Sync();
        return(0);
      };

      sizeextent operator-(const sizeextent &inExtent) const 
      {
        sizeextent returnExtent;
        if(inExtent.size() != this->size())
          return(returnExtent);
        returnExtent.Copy(*this);
        sizeextent::const_iterator inExtIt = inExtent.begin();
        sizeextent::iterator myIt          = returnExtent.begin();
        while(myIt != returnExtent.end()){
          std::pair<size_t,size_t> &curPair(*myIt++);
          const std::pair<size_t,size_t> &inPair(*inExtIt++);
          if(curPair.first < inPair.first)
            return(*this);
          curPair.first  -= inPair.first;
          curPair.second -= inPair.first;
        }
        return(returnExtent);
      };

      sizeextent &operator+=(const sizeextent &inExtent) 
      {
        if(inExtent.size() != this->size())
          return(*this);
        sizeextent::const_iterator inExtIt = inExtent.begin();
        sizeextent::iterator myIt          = this->begin();
        while(myIt != this->end()){
          std::pair<size_t,size_t> &curPair(*myIt++);
          const std::pair<size_t,size_t> &inPair(*inExtIt++);
          curPair.first  += inPair.first;
          curPair.second += inPair.first;
        }
        return(*this);
      };

      sizeextent RelativeTranslation(const sizeextent &inOrigin,const sizeextent &inDestination) const
      {
        sizeextent returnExtent;

        if(inOrigin.size() != inDestination.size())
          return(returnExtent);
        if(inOrigin.size() != this->size())
          return(returnExtent);

        returnExtent.Copy(*this);

        sizeextent::const_iterator originIt = inOrigin.begin();
        sizeextent::const_iterator destIt   = inDestination.begin();
        sizeextent::iterator returnIt       = returnExtent.begin();

        while(originIt != inOrigin.end()){
  
          const std::pair<size_t,size_t> &originPair(*originIt++);
          const std::pair<size_t,size_t> &destPair(*destIt++);
          std::pair<size_t,size_t> &returnPair(*returnIt++);
          
          returnPair.first  -= originPair.first;
          returnPair.second -= originPair.first;
          returnPair.first  += destPair.first;
          returnPair.second += destPair.first;

        }

        return(returnExtent);

      };

      // on success, returnExtent = (*this*[] - inOrigin[].first + inDestination[].first)\
      // the intent is to resolve an interval within a larger interval - useful in things like
      // domain decomposition where the partitioned interval is located inside a larger buffer 
      int RelativeTranslation(const sizeextent &inOrigin,const sizeextent &inDestination,sizeextent &returnExtent) const
      {

        if(inOrigin.size() != inDestination.size())
          return(1);
        if(inOrigin.size() != this->size())
          return(1);

        returnExtent.Copy(*this);

        sizeextent::const_iterator originIt = inOrigin.begin();
        sizeextent::const_iterator destIt   = inDestination.begin();
        sizeextent::iterator returnIt       = returnExtent.begin();

        while(originIt != inOrigin.end()){
  
          const std::pair<size_t,size_t> &originPair(*originIt++);
          const std::pair<size_t,size_t> &destPair(*destIt++);
          std::pair<size_t,size_t> &returnPair(*returnIt++);
          
          returnPair.first  -= originPair.first;
          returnPair.second -= originPair.first;
          returnPair.first  += destPair.first;
          returnPair.second += destPair.first;

        }

        return(0);

      };

      template<typename ContainerType>
      sizeextent Protruded(const ContainerType &protrusionOffsets)
      {
        sizeextent returnExtent;
        if(protrusionOffsets.size() != _nd*2)
          return(returnExtent);
        typename ContainerType::const_iterator offsetIt = protrusionOffsets.begin();
        sizeextent::iterator extentIt = this->begin();
        while(extentIt != this->end()){
          std::pair<size_t,size_t> &curPair(*extentIt++);
          if((*offsetIt < 0) && (-*offsetIt > curPair.first))
            return(returnExtent);
          offsetIt++;
          if((*offsetIt < 0) && (-*offsetIt > curPair.second))
            return(returnExtent);
          offsetIt++;
        }
        returnExtent.Copy(*this);
        offsetIt = protrusionOffsets.begin();
        extentIt = returnExtent.begin();
        while(extentIt != returnExtent.end()){
          std::pair<size_t,size_t> &curPair(*extentIt++);
          curPair.first  += *offsetIt++;
          curPair.second += *offsetIt++;
        }
        returnExtent.Sync();
        return(returnExtent);
      };

      
      template<typename ContainerType>
      int Translate(const ContainerType &translationOffsets,bool reverseOffset=false) 
      {
        if(translationOffsets.size() != _nd)
          return(1);
        sizeextent::iterator thisIt = this->begin();
        typename ContainerType::const_iterator transIt = translationOffsets.begin();
        while(thisIt != this->end()){
          std::pair<size_t,size_t> &curPair(*thisIt++);
          if(*transIt < 0){
            int transVal = std::abs(*transIt);
            if(transVal > curPair.first || transVal > curPair.second)
              return(1);
          }
          transIt++;
        }
        thisIt = this->begin();
        transIt = translationOffsets.begin();
        if(!reverseOffset){
          while(thisIt != this->end()){
            std::pair<size_t,size_t> &curPair(*thisIt++);
            curPair.first += *transIt;
            curPair.second += *transIt++;
          }
        } else {
          while(thisIt != this->end()){
            std::pair<size_t,size_t> &curPair(*thisIt++);
            curPair.first -= *transIt;
            curPair.second -= *transIt++;
          }
        }
        Sync();
        return(0);
      };

      template<typename ContainerType>
      void FindSharedNodes(std::vector<sizeextent> &extent_pool,std::vector<sizeextent> &shared_extents,
                           ContainerType &neighbors) const
      {
        shared_extents.resize(0);
        size_t nindex = 0;
        typename std::vector<sizeextent>::iterator epi = extent_pool.begin();
        while(epi != extent_pool.end()){
          if(*epi != *this){
            sizeextent overlap;
            this->Overlap(*epi,overlap);
            if(!overlap.empty()){
              neighbors.push_back(nindex);
              shared_extents.push_back(overlap);
            }
          }
          nindex++;
          epi++;
        }
      };

      // Count number of dimensions with no overlap
      int Skewness(const sizeextent &inExtent)
      {
        if(inExtent._nd != _nd)
          return(_nd);
        sizeextent::const_iterator extIt = inExtent.begin();
        sizeextent::iterator myIt = begin();
        int skewNum = 0;
        while(extIt != inExtent.end()){
          if((extIt->first > myIt->second) ||
             (extIt->second < myIt->first)) 
            skewNum++;
          extIt++;
          myIt++;
        }
        return(skewNum);
      };

      template<typename ContainerType>
      inline std::vector<size_t> Coordinates(size_t nodeIndex)
      {
        std::vector<size_t> retVal(_nd,0);
        std::vector<size_t>::reverse_iterator rvIt = retVal.rbegin();
        size_t newIndex = nodeIndex;
        for(int i = (_nd-1);i >= 0;i--){
          *rvIt = newIndex/_Np[i];
          newIndex = nodeIndex - (*rvIt * _Np[i]);
          rvIt++;
        }
      };
 
      template<typename ContainerType>
      inline size_t NodeIndex(ContainerType &index) const
      {
        size_t l0 = 0;
        typename ContainerType::reverse_iterator indIt = index.rbegin();
        for(int i = (_nd-1);i >= 0;i--)
          l0 += (*indIt++ - (*this)[i].first)*_Np[i]; 
        return(l0);
      };

      std::ostream &PrettyPrint(std::ostream &outStream) const;
    };
    
    std::ostream &operator<<(std::ostream &outStream,const ix::util::sizeextent &sizeExtent);
    std::istream &operator>>(std::istream &inStream,ix::util::sizeextent &sizeExtent);

  }
}
#endif