rtbulksort.cpp

此文件包含了在linux下实现tpr-tree索引的源代码

//--------------------------------------------------------------------
//      RTbulksort.cpp
//      --------------
//      Implements the array that can be sorted as required by 
//      STR bulk loading algorithms 
//
//      TPR-tree - Index for continuously moving objects
//      July 2001 release, Aalborg University
//

//#define NDEBUG     // if NDEBUG is defined asserts do not work

#include "rt.h"
#include <stdlib.h>

RTperiod RTparam;  
int      RToffset = 0;

int      RTsortedDim;
int      RTsortedSpeed;

// For HilbertSort
//
static       int  RTnumSortedDims;
static       int  RTmaxCoordExponent;
static const int  BITS_IN_MANTISSA = 22;    

// Slice exponent off from IEEE 754 single precision floating-point number
//
inline int floatExp(float f) { return (*(unsigned*)&f << 1) >> 24; }

//------------------------------------------------
extern "C" 
int CompareCoord  (const void *i, const void *j)
{
   if (((RTsortEntry*)i)->coords[RToffset + RTsortedDim] < 
       ((RTsortEntry*)j)->coords[RToffset + RTsortedDim]) return -1;
   if (((RTsortEntry*)i)->coords[RToffset + RTsortedDim] > 
       ((RTsortEntry*)j)->coords[RToffset + RTsortedDim]) return  1;
 
   return 0;
}

//------------------------------------------------
extern "C" 
int CompareCoordParam (const void *i, const void *j)
{
   RTcoord ic = ((RTsortEntry*)i)->coords[RToffset + RTsortedDim] + 
                ((RTsortEntry*)i)->coords[RToffset + RTsortedSpeed] * RTparam;
   RTcoord jc = ((RTsortEntry*)j)->coords[RToffset + RTsortedDim] + 
                ((RTsortEntry*)j)->coords[RToffset + RTsortedSpeed] * RTparam;

   if (ic < jc) return -1;
   if (ic > jc) return  1;
 
   return 0;
}

//------------------------------------------------
extern "C" 
int ComparePid    (const void *i, const void *j)
{
   if (((RTsortEntry*)i)->pid < ((RTsortEntry*)j)->pid) return -1;
   if (((RTsortEntry*)i)->pid > ((RTsortEntry*)j)->pid) return  1;
 
   return 0;
}

//---------------------------------------------------------------
//   CompareHilbert with auxilary routines

//------------------------------------------------
//   bitFromFloat returns the pos bit of x. 
//   Bit 0 is rightmost (least-significant) bit 
//

int bitFromFloat(int pos, float x)
{
   int mantPos;
   union {int i; float f;} u; 

   if (x == 0.0) return 0;

   mantPos = BITS_IN_MANTISSA + 1 + floatExp(x);

   if (pos >  mantPos) return 0;        // "Phantom" Leading Zeroes
   if (pos == mantPos) return 1;        // "Phantom" Leading One from normalized float 
   if (pos > mantPos - BITS_IN_MANTISSA - 1)
   {
      // Return the appropriate bit from the mantissa
      pos = mantPos - pos - 1;
      pos = BITS_IN_MANTISSA - pos;
      u.f = x;
      return u.i / (1 << pos ) % 2;
   }
   return 0;                            // Must be referencing a Trailing 0
}

//------------------------------------------------
inline unsigned char X_OR( unsigned char a, unsigned char b ) 
{ return ((a+b) & 1) ? 1: 0; }

//------------------------------------------------
inline void ROTATE( unsigned char a[], unsigned char len ) 
{
   int i; 
   static unsigned char b[TPR_MAX_DIMS*2];

   for (i = 0; i < RTnumSortedDims; i++) 
      b[(RTnumSortedDims + i + len) % RTnumSortedDims] = a[i]; 
   for (i = 0; i < RTnumSortedDims; i++) a[i] = b[i]; 
}

//------------------------------------------------
void HilbertDvalue( unsigned char in[], unsigned char h[],
                    unsigned char w[], unsigned char *rotation )
{
   unsigned char x, i, o[TPR_MAX_DIMS*2];
   char principal;

   /* Collect the j-digits over all coordinates */
   for (i = 0; i < RTnumSortedDims; i++) o[i] = X_OR(w[i], in[i]);

   /* rotate to the left by accumulated principals */
   if (*rotation) ROTATE(o, -(*rotation)); 

   /* decode the actual hilbert digits */
   for (i = 1, h[0] = o[0]; i < RTnumSortedDims; i++) h[i] = X_OR(h[i-1], o[i]);

   /* get principal pos so as to prepare the mask for next */
   for (principal = RTnumSortedDims - 2;
        principal >= 0 && h[principal] == h[RTnumSortedDims-1];
        principal--);
   if (principal == -1) principal += RTnumSortedDims;

   /* adjust the final mask */
   o[RTnumSortedDims - 1] = 1 - o[RTnumSortedDims - 1];
   for (i = 0, x = 0; i < RTnumSortedDims; i++) 
      if (o[i]) x = 1 - x;
   if (x == 1) o[principal] = 1 - o[principal];

   /* takes the rotation back (to the right) */
   if (*rotation) ROTATE(o, *rotation);
   *rotation = (*rotation + principal) % RTnumSortedDims;

   /* finally, prepare w for next run of j-digits */
   for (i = 0; i < RTnumSortedDims; i++) w[i] = X_OR(w[i], o[i]);
}

//------------------------------------------------
extern "C" 
int CompareHilbert (const void *i, const void *j)
{
   int d, k, h, Hdiff = 0;

   int limit;
   unsigned char Ibit[2][TPR_MAX_DIMS*2];
   unsigned char Hbit[2][TPR_MAX_DIMS*2];
   unsigned char mask[2][TPR_MAX_DIMS*2];
   unsigned char rotation[2];
   RTsortEntry *in[2] = { (RTsortEntry*)i, (RTsortEntry*)j };

   for (k = 0; k < 2; k++)
   {
      rotation[k] = 0;
      for (d = 0; d < RTnumSortedDims; d++)
      mask[k][d] = 0;
   } 
   limit = BITS_IN_MANTISSA + 1 + RTmaxCoordExponent;

   for (k = limit; !Hdiff && k >= 0; k--)
   {
      for (h = 0; h < 2; h++)
      {
         for (d = 0; d < RTnumSortedDims; d++)
            Ibit[h][d] = bitFromFloat(k, in[h]->coords[RToffset + d]);
         HilbertDvalue(Ibit[h], Hbit[h], mask[h], &rotation[h]);
      }

      for(h = 0; Hdiff == 0 && h < RTnumSortedDims; h++)
         Hdiff = Hbit[0][h] - Hbit[1][h];
   }

   return Hdiff;   
}

//----------------------------------------------------------------
//    RTsortArray
//

RTsortArray::RTsortArray(int ln, int rc_size) : 
             rec_size(rc_size), len(ln), num(0)
{
   // Do some tricks to achieve a proper address alignment for data
   data = (char*) new RTsortEntry[len * rec_size / sizeof(RTsortEntry) + 1];
   assert (data); 
}

//------------------------------------------------
RTsortArray::~RTsortArray()
{
   if (data) delete[] data; 
}

//------------------------------------------------
void* RTsortArray::operator[] (int n) const
{
   assert(data && n < num);
   return data + n * rec_size;
}

//------------------------------------------------
void* RTsortArray::NewRec ()
{
   assert(data && num < len);
   return data + (num++) * rec_size;
}

//------------------------------------------------
void  RTsortArray::Sort (int (*cmpf)(const void *, const void *))
{
   assert(data);
   qsort(data, num, rec_size, cmpf);
}

//------------------------------------------------
void  RTsortArray::SortNTimes (int from, int numRecs, int maxRecs, double alpha,
                               int (*cmpf)(const void *, const void *))
{
   int   numSorts;
   int   recsPerSlice;
   int   nodesPerSlice;
   int   numNodes;
   int   nRoot;    // nth root of the number of pages where n = dimension
   int   i, j;        
   int   numDims = Settings.GetLoadAlg() == RTsettings::aNormal || 
                   Settings.GetLoadAlg() == RTsettings::aVelocity ?
                   Settings.GetDims() : Settings.GetDims()*2;
   RTcoord  uexts[TPR_MAX_DIMS*2];
   RTcoord  s;    // length of an edge of the spatial part of a BR 
   int      dimOrder[TPR_MAX_DIMS*2]; 
   char*    dataptr = data + from * rec_size; 
 
   assert(data);

   // Prepare dimOrder - order of dimensions
   if (Settings.GetLoadAlg() == RTsettings::aDualSpeed || 
       Settings.GetLoadAlg() == RTsettings::aVelocity) 
   {
      for (i = 0, j = Settings.GetDims(); i < Settings.GetDims(); j++, i++) dimOrder[i] = j;
      for (j = 0; j < Settings.GetDims(); j++, i++) dimOrder[i] = j;
   }
   else if (Settings.GetLoadAlg() == RTsettings::aCSD)
   {  
      for (i = 0; i < Settings.GetDims(); i++) dimOrder[i] = i;
      dimOrder[i++] = Settings.GetDims() * 3 - 1;
      for (j = Settings.GetDims()*2; j < Settings.GetDims()*3 - 1; i++, j++) 
         dimOrder[i] = j; 
   }
   else if (Settings.GetLoadAlg() == RTsettings::aCDS)
   {  
      for (i = 0; i < Settings.GetDims(); i++) dimOrder[i] = i;
      for (j = Settings.GetDims()*2; j < Settings.GetDims()*3; i++, j++) dimOrder[i] = j; 
   }
   else if (Settings.GetLoadAlg() == RTsettings::aDSC)
   {  
      for (i = 0, j = Settings.GetDims()*2; j < Settings.GetDims()*3; i++, j++) 
         dimOrder[i] = j; 
      for (j = 0; j < Settings.GetDims(); j++, i++) dimOrder[i] = j;
   }
   else if (Settings.GetLoadAlg() == RTsettings::aSDC)
   {  
      dimOrder[0] = Settings.GetDims() * 3 - 1;
      for (i = 1, j = Settings.GetDims()*2; j < Settings.GetDims()*3 - 1; i++, j++) 
         dimOrder[i] = j; 
      for (j = 0; j < Settings.GetDims(); j++, i++) dimOrder[i] = j;
   }
   else   // RTsettings::aDual, RTsettings::aNormal, RTsettings::aRatio
      for (i = 0; i < Settings.GetDims()*2; i++) dimOrder[i] = i;  

   numNodes = (int)ceil((double)numRecs/maxRecs);
   nRoot    = (int)ceil(pow(numNodes, 1.0/numDims));

   // For ratio based algorithm, we have to find the extents of the universe  
   if (Settings.GetLoadAlg() == RTsettings::aRatio)
   {
      RTsortEntry *mincp, *maxcp;
      char        *cp;

      for (RTsortedDim = 0, s = 1; RTsortedDim < numDims; RTsortedDim++) 
      {
         mincp = maxcp = (RTsortEntry*) dataptr;
         for (i = 0, cp = dataptr; i < numRecs; i++, cp += rec_size)
	      {
            if (cmpf(cp, mincp) < 0) mincp = (RTsortEntry*) cp;
            if (cmpf(cp, maxcp) > 0) maxcp = (RTsortEntry*) cp; 
     	   }
         uexts[RTsortedDim] = maxcp->coords[RToffset + RTsortedDim] - 
                              mincp->coords[RToffset + RTsortedDim];
         s *= uexts[RTsortedDim];
      }
      s = pow (s / numNodes, 1.0 / numDims) / sqrt (alpha);

      // prepare uexts, to store the # of bands for each dimension (# of cuts + 1)
      for (j = 0; j < numDims; j++)
      {
         uexts[j] /= j < Settings.GetDims() ? s : s * alpha;
         if (uexts[j] < 1) uexts[j] = 1;
      }   
   }

   // Sort on the First Coordinate of the Center Point
   RTsortedDim = dimOrder[0];
   RTsortedSpeed = Settings.GetDims() + RTsortedDim;
   qsort(dataptr, numRecs, rec_size, cmpf);
  
   nodesPerSlice = numNodes;
   // Perform the sorts on each coordinate
   for(j = 1; j < numDims; j++)
   {
      RTsortedDim   = dimOrder[j];
      RTsortedSpeed = Settings.GetDims() + RTsortedDim;

      if (Settings.GetLoadAlg() == RTsettings::aRatio) 
      {
         nodesPerSlice = (int) floor(nodesPerSlice / uexts[dimOrder[j - 1]] + 0.5);
         if (!nodesPerSlice) nodesPerSlice = 1;
      }
      else nodesPerSlice = (int) pow(nRoot, numDims - j);
      recsPerSlice = nodesPerSlice * maxRecs ;
      numSorts     = (int) ceil((double)numRecs / recsPerSlice);

      for (i = 0; i < numSorts - 1; i++)
         qsort(dataptr + i*recsPerSlice*rec_size, recsPerSlice, rec_size, cmpf);
 
      // Sort the last column
      qsort(dataptr + i*recsPerSlice*rec_size, numRecs - i*recsPerSlice,
           rec_size, cmpf);
   }
}

//------------------------------------------------
//  Sorting along Hilbert space filling curve
//

void RTsortArray::HilbertSort()
{
    RTcoord  shifts[TPR_MAX_DIMS*2];
    RTcoord  c; 
    bool     needMove = false;     // need for a shift of a coordinate center
    int      i, j, ex;   
   
    RTnumSortedDims = Settings.GetLoadAlg() == RTsettings::aNormalHilbert ? 
                      Settings.GetDims() : Settings.GetDims()*2;
    
    for (i = 0; i < RTnumSortedDims; i++) shifts[i] = 0;
    RTmaxCoordExponent = 0;   
  
    // Move the coordinate center to a position were all coordinates are positive
    // and find maximum coordinate exponent
    for (i = 0; i < num; i++)
       for (j = 0; j < RTnumSortedDims; j++)
       {
          c = ((RTsortEntry*)(*this)[i])->coords[RToffset + j];
          if (c < shifts[j]) shifts[j] = c;
          ex = floatExp(c); 
          assert(ex || c == 0.0); 
          if (ex > RTmaxCoordExponent) RTmaxCoordExponent = ex;
       }
    for (i = 0; i < RTnumSortedDims && !needMove; i++) 
       if (shifts[i] != 0.0) needMove = true;
    if (needMove)
       for (i = 0; i < num; i++)
          for (j = 0; j < RTnumSortedDims; j++)
             ((RTsortEntry*)(*this)[i])->coords[RToffset + j] -= shifts[j];
    
    qsort(data, num, rec_size, CompareHilbert);
   
    // Move the coordinate center to the original position
    if (needMove)
       for (i = 0; i < num; i++)
          for (j = 0; j < RTnumSortedDims; j++)
             ((RTsortEntry*)(*this)[i])->coords[RToffset + j] += shifts[j];    
}

//------------------------------------------------
void RTsortArray::ReadPoints (ifstream& is)
{
   if (data) delete[] data;
   
   rec_size = SIZEOF_LEAFREC; 
   
   is.seekg(0, ios::end);
   len = num = is.tellg() / rec_size;
   is.seekg(0, ios::beg);

   // Do some tricks to achieve a proper address alignment for data
   data = (char*) new RTsortEntry[len * rec_size / sizeof(RTsortEntry) + 1];
   is.read(data, len * rec_size);   
}

//------------------------------------------------
void RTsortArray::AdjustPoints (RTperiod deltaT)
{
   int i, j;
  
   for (i = 0; i < num; i++)
      for (j = 0; j < Settings.GetDims(); j++) 
        ((RTsortEntry*)(*this)[i])->coords[j] += 
        ((RTsortEntry*)(*this)[i])->Velocity(j) * deltaT;
}