src/util/hilbert.C

   1 #include "hilbert.h"
   2
   3 #ifdef __cplusplus
   4 #include <queue>
   5 #include <vector>
   6 #include <iostream>
   7 #include <sstream>
   8 #include <algorithm>
   9 #include <math.h>
  10
  11 using namespace std;
  12
  13 #ifndef PARTITION_TOPOLOGY_VERBOSE
  14 #define PARTITION_TOPOLOGY_VERBOSE 0
  15 #endif
  16
  17 /**
  18  *  Author: Harshitha Menon, Nikhil Jain, Yanhua Sun
  19  *  Contact: gplkrsh2@illinois.edu, nikhil@illinois.edu, sun51@illinois.edu
  20  *
  21  *  More details about this implementation of the Hilbert curve can be found
  22  *  from https://github.com/straup/gae-spacetimeid/blob/master/hilbert.py
  23  *  and this is a C++ implementation of what is given there.
  24  */
  25 static int gray_encode(int i) {
  26   //cout << "gray_encode " << i << " " << (i^(i/2)) << endl;
  27   return i ^ (i/2);
  28 }
  29
  30 static int gray_decode(int n) {
  31   int sh = 1;
  32   int div;
  33   while (true) {
  34     div = n >> sh;
  35     n ^= div;
  36     if (div <= 1) {
  37       return n;
  38     }
  39     sh <<=1;
  40   }
  41 }
  42
  43 static void initial_start_end(int nChunks, int dim, int& start, int& end) {
  44   start = 0;
  45   ////cout << "-nchunks - 1 mod dim " << ((-nChunks-1)%dim) << endl;
  46   int mod_val = (-nChunks - 1) % dim;
  47   if (mod_val < 0) {
  48     mod_val += dim;
  49   }
  50   end = pow((double)2, (double)mod_val);
  51 }
  52
  53 static int pack_index(vector<int> chunks, int dim) {
  54   int p = pow((double)2, (double)dim);
  55   int chunk_size = chunks.size();
  56   int val = 0;
  57   for (int i = 0;i < chunk_size; i++) {
  58       val = val*p + chunks[i] ;
  59   }
  60   return val;
  61 }
  62
  63 static vector<int> transpose_bits(vector<int> srcs, int nDests) {
  64   int nSrcs = srcs.size();
  65   vector<int> dests;
  66   dests.resize(nDests);
  67   int dest = 0;
  68   for (int j = nDests - 1; j > -1; j--) {
  69     dest = 0;
  70     ////cout << "nSrcs " << nSrcs << endl;
  71     for (int k = 0; k < nSrcs; k++) {
  72       dest = dest * 2 + srcs[k] % 2;
  73       srcs[k] /= 2;
  74       ////cout << "dest " << dest << endl;
  75     }
  76     dests[j] = dest;
  77   }
  78   return dests;
  79 }
  80
  81 static vector<int> pack_coords(vector<int> coord_chunks, int dim) {
  82   return transpose_bits(coord_chunks, dim);
  83 }
  84
  85 static vector<int> unpack_coords(const vector<int> &coords, int dim) {
  86     int biggest = 0;
  87     for(int i=0; i<coords.size(); i++)
  88     {
  89         if(coords[i] > biggest)
  90             biggest = coords[i];
  91     }
  92     int nChunks = max(1, int( ceil( log( (double)biggest + 1)/log(2.0f) ) ) );
  93     return transpose_bits( coords, nChunks );
  94 }
  95
  96
  97 static void unpack_index(int i, int dim, vector<int>& chunks) {
  98   int p = pow((double)2, (double)dim);
  99   int nChunks = max(1, int(ceil(double(log((double)i+1))/log((double)p))));
 100   //cout << "num chunks " << nChunks << endl;
 101   chunks.resize(nChunks);
 102   for (int j = nChunks-1; j > -1; j--) {
 103     chunks[j] = i % p;
 104     i /= p;
 105     ////cout << "chunks[" << j << "] " << chunks[j] << endl;
 106   }
 107   //cout << "chunk size " << chunks.size() << endl;
 108 }
 109
 110 static int gray_encode_travel(int start, int end, int mask, int i) {
 111   int travel_bit = start ^ end;
 112   int modulus = mask + 1;
 113   int g = gray_encode(i) * (travel_bit * 2);
 114   //cout << "start " << start << " end " << end << "travel_bits " << travel_bit << " modulus " << modulus << " g " << g << endl;
 115   return ((g | (g/modulus) ) & mask) ^ start;
 116 }
 117
 118 static int gray_decode_travel(int start, int end, int mask, int i) {
 119   int travel_bit = start ^ end;
 120   int modulus = mask + 1;
 121   int rg = (i ^ start) * ( modulus/(travel_bit * 2));
 122   return gray_decode((rg | (rg / modulus)) & mask);
 123 }
 124
 125 static void child_start_end(int parent_start, int parent_end, int mask, int i, int&
 126     child_start, int& child_end) {
 127   int start_i = max(0, (i-1)&~1);
 128   //cout << "start_i " << start_i << endl;
 129   int end_i = min(mask, (i+1)|1);
 130   //cout << "end_i " << end_i << endl;
 131   child_start = gray_encode_travel(parent_start, parent_end, mask, start_i);
 132   child_end = gray_encode_travel(parent_start, parent_end, mask, end_i);
 133   //cout << "child_start " << child_start << " child end " << child_end << endl;
 134 }
 135
 136 vector<int> int_to_Hilbert(int i, int dim) {
 137   int nChunks, mask, start, end;
 138   vector<int> index_chunks;
 139   unpack_index(i, dim, index_chunks);
 140   nChunks = index_chunks.size();
 141   //cout << "int to hilbert of " << i << " in dim " << dim << " size " << nChunks << endl;
 142   mask = pow((double)2, (double)dim) - 1;
 143   //cout << "mask " << mask << endl;
 144   initial_start_end(nChunks, dim, start, end);
 145   //cout << "start " << start << " end " << end << endl;
 146   vector<int> coord_chunks;
 147   coord_chunks.resize(nChunks);
 148   for (int j = 0; j < nChunks; j++) {
 149     i = index_chunks[j];
 150     coord_chunks[j] = gray_encode_travel(start, end, mask, i);
 151     //cout << "coord_chuunk " << j << " : " << coord_chunks[j] << endl;
 152     ////cout << "going for child start end" << endl;
 153     child_start_end(start, end, mask, i, start, end);
 154     //cout << "child start end " << start << "  " << end << endl;
 155   }
 156   return pack_coords(coord_chunks, dim);
 157 }
 158
 159 int Hilbert_to_int(const vector<int>& coords, int dim)
 160 {
 161     vector<int> coord_chunks =  unpack_coords( coords, dim );
 162     int i;
 163     int nChunks = coord_chunks.size();
 164     int mask = pow((double)2, (double)dim) - 1;
 165     int start, end;
 166     initial_start_end(nChunks, dim, start, end);
 167     vector<int> index_chunks;
 168     index_chunks.resize(nChunks);
 169     for (int j = 0; j < nChunks; j++) {
 170         i = gray_decode_travel( start, end, mask, coord_chunks[ j ] );
 171         index_chunks[ j ] = i;
 172         child_start_end(start, end, mask, i, start, end);
 173     }
 174
 175     return pack_index( index_chunks, dim );
 176 }
 177
 178 #endif