HilbertCurveBatchSort.java

/*
 * Class:        HilbertCurveBatchSort
 * Description:  performs a batch sort on the arrays
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2015  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since
 
 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.
 
 * SSJ 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.
 
 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
package umontreal.ssj.util.multidimsort;

public class HilbertCurveBatchSort<T extends MultiDimComparable<? super T>> implements MultiDimSortComparable<T> {
   int dimension; // Dimension d of the points used for the sort.
   int m; // Number of bit retained for each coordinate. Must be <= 31.
   HilbertCurveMap hcMap;
 
   BatchSortPow2<T> bsort; // Batch sort used.
   long[][] indexH; // Table whose each line gives the number of an object after the
                    // batch sort in first coord. and its Hilbert index in second coord.
                    // Here, the second coordinates are all distinct!
                    // For a given n, this index is computed once for all and saved.
   int nSavedIndex = 0; // Value of n for which indexH is now computed.

   public HilbertCurveBatchSort(double[] batchExponents, int m) {
      hcMap = new HilbertCurveMap(batchExponents.length, m);
      dimension = batchExponents.length;
      this.m = m;
      bsort = new BatchSortPow2<T>(batchExponents);
   }

   public HilbertCurveBatchSort(double[] batchExponents, int m, HilbertCurveMap map) {
      dimension = batchExponents.length;
      this.m = m;
      if (map.dimension() != dimension)
         throw new IllegalArgumentException("HilbertCurveMap has a different dimension! Expecting: " + dimension);
      this.hcMap = map;
      bsort = new BatchSortPow2<T>(batchExponents);
   }

   public void computeIndexH(int n) {
      if (n == nSavedIndex)
         return;
      nSavedIndex = n;
      indexH = new long[n][2]; // To put Hilbert indexes of points.
      // For the following, we need e_0, ..., e_{d-1}.
      int twom = (1 << m); // 2^m.
      int mdim = dimension - 1;
      int[] ej = bsort.getBitNumbers();
      int stepj[] = new int[dimension]; // The integers 2^{m-e_j}.
      int[] icoord = new int[dimension]; // The integer coordinates of current point.
      for (int j = 0; j < dimension; ++j) {
         if (ej[j] > m)
            throw new RuntimeException("ej[j] is larger than m");
         icoord[j] = 0;
         stepj[j] = 1 << (m - ej[j]);
      }
      indexH[0][0] = 0;
      indexH[0][1] = 0;
      for (int i = 1; i < n; ++i) {
         icoord[mdim] += stepj[mdim]; // Increment last coordinate.
         if (icoord[mdim] == twom) { // Must propagate the carry.
            for (int j = mdim; j >= 0; --j) {
               if (icoord[j] < twom)
                  break; // Exit for loop if no carry.
               icoord[j] = 0;
               ++icoord[j - 1];
            }
         }
         indexH[i][0] = i;
         indexH[i][1] = hcMap.coordinatesToIndex(icoord);
      }
 
      // use static method of HilbertCurveSort to sort
      HilbertCurveSort.sortIndexOfLong2(indexH, 0, n); // First coordinate of indexH gives the ordering
      // that must be applied.
   }

   public void sort(T[] a, int iMin, int iMax) {
      if (iMin + 1 == iMax)
         return;
      if (iMax - iMin != nSavedIndex)
         computeIndexH(iMax - iMin); // This changes nSavedIndex.
      bsort.sort(a, iMin, iMax); // This changes nSaved if needed.
      // Now use the sorted indexH to sort a.
      // First make a copy, and pick the elements from this unmodified copy.
      T[] aclone = a.clone(); // We should not clone all the objects in a,
                              // but only the array of pointers.
      // T[] aclone = new Object[iMax];
 
      for (int i = iMin; i < iMax; ++i) {
         a[i] = aclone[(int) indexH[i][0]];
      }
   }

   public void sort(T[] a) {
      sort(a, 0, a.length);
   }
 
   public void sort(double[][] a, int iMin, int iMax) {
      if (iMin + 1 == iMax)
         return;
      if (iMax - iMin != nSavedIndex)
         computeIndexH(iMax - iMin); // This changes nSavedIndex.
      bsort.sort(a, iMin, iMax); // This changes nSaved if needed.
      // Now use the sorted indexH to sort a.
      // First make a copy, and pick the elements from this unmodified copy.
      double[][] aclone = a.clone();
      for (int i = iMin; i < iMax; ++i) {
         a[i] = aclone[(int) indexH[i][0]];
      }
   }
 
   public void sort(double[][] a) {
      sort(a, 0, a.length);
   }

   public int dimension() {
      return dimension;
   }

   public HilbertCurveMap getHilbertCurveMap() {
      return hcMap;
   }
 
}