PCPPIntervalUtils.C

Go to the documentation of this file.

#include "PCPPIntervalUtils.H"

namespace pcpp {
  namespace util {


    int SubIntervalFromStream(std::istream &inStream,const pcpp::IndexIntervalType &inInterval,
                              pcpp::IndexIntervalType &subInterval)
    {
      subInterval.Copy(inInterval);
      int numDim = inInterval.size();
      int ibDir = 0;
      std::vector<int> iStart(numDim,0);
      std::vector<int> iEnd(numDim,0);
      inStream >> ibDir;
      if(ibDir == 0)
        return(1);
      for(int iDim = 0;iDim < numDim;iDim++){
        inStream >> iStart[iDim] >> iEnd[iDim];
        if(iStart[iDim] == 0 || iEnd[iDim] == 0)
          return(1);
      }
      int direction = ibDir;
      int leftBoundary = 1;
      if(ibDir < 0){
        direction = -1*direction;
        leftBoundary = 0;
      }
      direction = direction - 1;
      if(leftBoundary){
        subInterval[direction].first = iStart[direction]-1;
        subInterval[direction].second = iEnd[direction]-1;
      } else {
        int startIndex = iStart[direction];
        int endIndex   = iEnd[direction];
        if(iStart[direction] < 0)
          startIndex *= -1;
        if(iEnd[direction] < 0)
          endIndex *= -1;
        startIndex -= 1;
        endIndex -= 1;
        subInterval[direction].first  -= startIndex;
        subInterval[direction].second -= endIndex;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        if(iDim != direction){
          int startIndex = iStart[iDim];
          int endIndex   = iEnd[iDim];
          if(startIndex < 0){
            startIndex += 1;
            startIndex *= -1;
            subInterval[iDim].first = inInterval[iDim].second - startIndex;
          } else {
            startIndex -= 1;
            subInterval[iDim].first = inInterval[iDim].first + startIndex;
          }
          if(endIndex < 0){
            endIndex += 1;
            endIndex *= 1;
            subInterval[iDim].second = inInterval[iDim].second - endIndex;
          } else {
            endIndex -= 1;
            subInterval[iDim].second = inInterval[iDim].first + endIndex;
          }
        }
      }
      // Check to make sure it actually overlaps
      pcpp::IndexIntervalType overlapInterval(inInterval.Overlap(subInterval));
      if(overlapInterval != subInterval)
        return(1);
      return(0);
    };

    std::vector<pcpp::IndexIntervalType> 
    UniqueUnion(const std::vector<pcpp::IndexIntervalType> &intervalSet)
    {
      std::vector<pcpp::IndexIntervalType> uniqueUnion(intervalSet);
      int numIntervals = intervalSet.size();
      if(numIntervals <= 1)
        return(uniqueUnion);
      for(int iSet = 0;iSet < numIntervals-1;iSet++){
        pcpp::IndexIntervalType interval1(uniqueUnion[iSet]);
        for(int iCompare = iSet+1;iCompare < numIntervals;iSet++){
          pcpp::IndexIntervalType compareInterval(interval1.Overlap(uniqueUnion[iCompare]));
          if(compareInterval.NNodes() > 0){
            
          }
        }
      }
      
    }

    void CollapseInterval(size_t &xSize,size_t &ySize,size_t &zSize,
                          const pcpp::IndexIntervalType &inInterval,
                          pcpp::IndexIntervalType &outInterval)
    {
      outInterval = inInterval;
      if(xSize == 1){
        if(ySize == 1){
          outInterval[0] = outInterval[2];
          xSize = zSize;
          outInterval[1].first = 0;
          outInterval[1].second = 0;
          outInterval[2] = outInterval[1];
        } else {
          outInterval[0] = outInterval[1];
          xSize = ySize;
          outInterval[1] = outInterval[2];
          ySize = zSize;
          outInterval[2].first = 0;
          outInterval[2].second = 0;
          zSize = 1;
        } 
      } else if (ySize == 1){
        if(zSize == 1){
          outInterval[1].first = 0;
          outInterval[1].second = 0;
          outInterval[2] = outInterval[1];
        } else {
          outInterval[1] = outInterval[2];
          outInterval[2].first = 0;
          outInterval[2].second = 0;
          ySize = zSize;
          zSize = 1;
        }
      } else if (zSize == 1){
        outInterval[2].first = 0;
        outInterval[2].second = 0;
      }
    };

    std::deque<size_t> PrimeFactors(size_t inNumber)
    {
      std::deque<size_t> primeFactors;
      size_t baseFac = 2;
      while(baseFac*baseFac <= inNumber){
        if(inNumber%baseFac == 0){
          primeFactors.push_back(baseFac);
          inNumber /= baseFac;
        } else {
          baseFac++;
        }
      }
      if(inNumber >= 1)
        primeFactors.push_back(inNumber);
      return(primeFactors);
    };

    int SimplePartitionInterval(const pcpp::IndexIntervalType &inInterval,
                                std::vector<bool> partDirection,
                                int partID,int numPart,
                                pcpp::IndexIntervalType &outInterval,
                                std::vector<int> &numPartitions)
    {

      if(numPart == 1){
        outInterval = inInterval;
        return(0);
      }

      outInterval.resize(0);

      if(partID < 0 || numPart < 1){
        std::cout << "numpart prob" << std::endl;
        return(-1);
      }

      int numDim = inInterval.size();

      if(numDim <= 0){
        std::cout << "numdim prob" << std::endl;
        return(-1);
      }

      if(partDirection.empty()){
        partDirection.resize(numDim,true);
      } else if (partDirection.size() != numDim){
        std::cout << "partdir prob" << std::endl;
        return(-1);
      }
      
      size_t maxTotalPossibleParts = 1;
      std::vector<int> candidateDims;
      std::vector<size_t> candidateSize;
      numPartitions.resize(numDim,0);

      // Initialize some needed data structures 
      for(int iDim = 0;iDim < numDim;iDim++){
        if(partDirection[iDim]){
          size_t mySize = inInterval[iDim].second - inInterval[iDim].first + 1;
          candidateDims.push_back(iDim);
          candidateSize.push_back(mySize);
          maxTotalPossibleParts *= mySize;
        } else {
          numPartitions[iDim] = 1;
        }
      }

      if(numPart > maxTotalPossibleParts){
        std::cout << "maxpart" << std::endl; 
        return(-1);
      }

      int partDims = candidateDims.size();
      outInterval = inInterval;

      // prime factorization of desired number of partitions
      std::deque<size_t> primeFactors(PrimeFactors(numPart));
      int numFac = primeFactors.size();

      if(numFac == 0)
        return(0);
      if(numFac == 1 && primeFactors.front() == 1)
        return(0);

      // If the number of prime factors is larger than the 
      // number of possible dimenions, reduce by recursively
      // combining (i.e. multiplying) the two smallest factors.
      // Resulting deque is sorted smallest to largest
      while(numFac > partDims){
        primeFactors[1] *= primeFactors.front();
        primeFactors.pop_front();
        std::sort(primeFactors.begin(),primeFactors.end());
        numFac--;
      }

      // Next, match the largest partitions to the dimensions with
      // largest size until all the dimensions are partitioned.
      while(!primeFactors.empty()){
        int dimPart = primeFactors.back();
        primeFactors.pop_back();
        int dimPlane = 1;
        if(!primeFactors.empty()){
          std::deque<size_t>::iterator facIt = primeFactors.begin();
          while(facIt != primeFactors.end())
            dimPlane *= *facIt++;
        }

        int maxSize = 0;
        int maxDim = -1;
        for(int iCand = 0;iCand < partDims;iCand++){
          if((candidateSize[iCand] >= maxSize) && (numPartitions[candidateDims[iCand]] == 0)){
            maxDim = candidateDims[iCand];
            maxSize = candidateSize[iCand];
          }
        }
        if(maxSize == 0 || maxDim < 0){
          std::cout << "maxsize = " << maxSize << "maxDim = " << maxDim << std::endl;
          return(-1);
        }

        numPartitions[maxDim] = dimPart;
        size_t iStart = outInterval[maxDim].first;
        size_t iEnd   = outInterval[maxDim].second;
        size_t partRank = (partID/dimPlane)%dimPart;
        Part1D(iStart,iEnd,dimPart,partRank,outInterval[maxDim].first,
               outInterval[maxDim].second);
      }

      return(0);

    };
    

    int Part1D(const size_t iStart,const size_t iEnd,const size_t numPart,const size_t partIndex,
               size_t &partStart,size_t &partEnd)
    {
      size_t intervalSize = iEnd - iStart + 1;
      size_t partSize = intervalSize/numPart;
      if(partSize == 0)
        return(1);
      size_t leftOver = intervalSize%numPart;
      if(partIndex < leftOver)
        partSize++;
      size_t intervalStart = partIndex*partSize;
      if(partIndex >= leftOver)
        intervalStart += leftOver;
      partStart = intervalStart + iStart;
      partEnd = partStart + partSize - 1;
      return(0);
    };

    int PartitionCartesianExtent(const pcpp::IndexIntervalType &globalExtent,
                                 const std::vector<int> &cartDims,
                                 const std::vector<int> &cartCoords,
                                 pcpp::IndexIntervalType &partExtent,
                                 std::ostream &messageStream)
    {

      if(globalExtent.empty()){
        messageStream << "PartitionCartesianExtent:Error: Global extent empty." 
                      << std::endl;
        return(1);
      }
      int numDim = cartDims.size();
      if(globalExtent.ND() < numDim){
        messageStream << "PartitionCartesianExtent:Error: Global extent dimension deficient, "
                      << globalExtent.ND() << " < " << numDim << "." << std::endl;
        return(1);
      }
      std::vector<size_t> localStart(numDim,0);
      std::vector<size_t> localSize(numDim,0);
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t numPart = cartDims[iDim];
        size_t intervalStart = globalExtent[iDim].first;
        size_t intervalEnd   = globalExtent[iDim].second;
        size_t partStart = 0;
        size_t partEnd   = 0;
        size_t partIndex = cartCoords[iDim];
        if(Part1D(intervalStart,intervalEnd,numPart,partIndex,partStart,partEnd)){
          messageStream << "PartitionCartesianExtent:Error: Global extent of insufficient size for number of partitions."
                        << std::endl;
          return(1);
        }
        localSize[iDim]  = partEnd - partStart + 1;
        localStart[iDim] = partStart;
      }
      partExtent.Init(localStart,localSize);
      return(0);
    };
    int PartitionCartesianInterval(const pcpp::IndexIntervalType &globalInterval,
                                   const std::vector<int> &cartDims,
                                   const std::vector<int> &cartCoords,
                                   pcpp::IndexIntervalType &partInterval,
                                   std::ostream &messageStream)
    {

      if(globalInterval.empty()){
        messageStream << "PartitionCartesianInterval:Error: Global interval empty." 
                      << std::endl;
        return(1);
      }
      int numDim = cartDims.size();
      if(globalInterval.ND() < numDim){
        messageStream << "PartitionCartesianInterval:Error: Global interval dimension deficient, "
                      << globalInterval.ND() << " < " << numDim << "." << std::endl;
        return(1);
      }
      std::vector<size_t> localStart(numDim,0);
      std::vector<size_t> localSize(numDim,0);
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t numPart = cartDims[iDim];
        size_t intervalStart = globalInterval[iDim].first;
        size_t intervalEnd   = globalInterval[iDim].second;
        size_t partStart = 0;
        size_t partEnd   = 0;
        size_t partIndex = cartCoords[iDim];
        if(Part1D(intervalStart,intervalEnd,numPart,partIndex,partStart,partEnd)){
          messageStream << "PartitionCartesianInterval:Error: Global interval of insufficient size for number of partitions."
                        << std::endl;
          return(1);
        }
        localSize[iDim]  = partEnd - partStart + 1;
        localStart[iDim] = partStart;
      }
      partInterval.Init(localStart,localSize);
      return(0);
    };
  };
};