1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.numerics.math;
39
40 import org.apache.commons.math3.util.FastMath;
41
42 import java.math.BigInteger;
43 import java.util.Arrays;
44
45
46
47
48
49
50
51
52
53 public class HilbertCurveTransforms {
54
55 private HilbertCurveTransforms() {
56
57 }
58
59 private static int adjust_rotation(int rotation, int nDims, int bits) {
60 long nd1Ones = (ones(nDims) >> 1);
61 bits &= (int) (-bits & nd1Ones);
62 while (bits != 0) {
63 bits >>= 1;
64 ++rotation;
65 }
66 if (rotation >= nDims) {
67 rotation -= nDims;
68 }
69 return rotation;
70 }
71
72 private static long ones(int k) {
73 return ((2L << (k - 1)) - 1);
74 }
75
76 private static long rotateLeft(long arg, int nRots, int nDims) {
77 return (((arg << nRots) | (arg >> (nDims - nRots))) & ones(nDims));
78 }
79
80 private static long rotateRight(long arg, int nRots, int nDims) {
81 return ((arg >> nRots) | (arg << (nDims - nRots))) & ones(nDims);
82 }
83
84 private static long bitTranspose(int nDims, int nBits, long inCoords) {
85 int nDims1 = nDims - 1;
86 int inB = nBits;
87 int utB;
88 long inFieldEnds = 1;
89 long inMask = ones(inB);
90 long coords = 0;
91
92 while ((utB = inB / 2) != 0) {
93 int shiftAmt = nDims1 * utB;
94 long utFieldEnds = inFieldEnds | (inFieldEnds << (shiftAmt + utB));
95 long utMask = (utFieldEnds << utB) - utFieldEnds;
96 long utCoords = 0;
97 int d;
98 if (inB % 2 == 1) {
99 long inFieldStarts = inFieldEnds << (inB - 1);
100 int oddShift = 2 * shiftAmt;
101 for (d = 0; d < nDims; ++d) {
102 long in = inCoords & inMask;
103 inCoords >>= inB;
104 coords |= (in & inFieldStarts) << oddShift++;
105 in &= ~inFieldStarts;
106 in = (in | (in << shiftAmt)) & utMask;
107 utCoords |= in << (d * utB);
108 }
109 } else {
110 for (d = 0; d < nDims; ++d) {
111 long in = inCoords & inMask;
112 inCoords >>= inB;
113 in = (in | (in << shiftAmt)) & utMask;
114 utCoords |= in << (d * utB);
115 }
116 }
117 inCoords = utCoords;
118 inB = utB;
119 inFieldEnds = utFieldEnds;
120 inMask = utMask;
121 }
122 coords |= inCoords;
123 return coords;
124 }
125
126
127
128
129
130
131
132
133
134 public static long[] hilbertIndexToCoordinates(int nBonds, int nBitsPerDim, long index) {
135 long[] coord = new long[nBonds];
136
137 if (nBonds > 1) {
138 long coords;
139 int nbOnes = (int) ones(nBitsPerDim);
140 int d;
141
142 if (nBitsPerDim > 1) {
143 int nDimsBits = nBonds * nBitsPerDim;
144 int ndOnes = (int) ones(nBonds);
145 int nd1Ones = ndOnes >> 1;
146 int b = nDimsBits;
147 int rotation = 0;
148 int flipBit = 0;
149 long nthbits = ones(nDimsBits) / ndOnes;
150 index ^= (index ^ nthbits) >> 1;
151 coords = 0;
152
153 do {
154 int bits = (int) ((index >> (b -= nBonds)) & ndOnes);
155 coords <<= nBonds;
156 coords |= rotateLeft(bits, rotation, nBonds) ^ flipBit;
157 flipBit = 1 << rotation;
158 rotation = adjust_rotation(rotation, nBonds, bits);
159 } while (b > 0);
160
161 for (b = nBonds; b < nDimsBits; b *= 2) {
162 coords ^= coords >> b;
163 }
164 coords = bitTranspose(nBitsPerDim, nBonds, coords);
165 } else {
166 coords = index ^ (index >> 1);
167 }
168
169 for (d = 0; d < nBonds; ++d) {
170 coord[d] = coords & nbOnes;
171 coords >>= nBitsPerDim;
172 }
173 } else {
174 coord[0] = index;
175 }
176
177 return coord;
178 }
179
180 private static long coordinatesToHilbertIndex(int nDims, int nBits, long[] coord) {
181 long index;
182 if (nDims > 1) {
183 int nDimsBits = nDims * nBits;
184 long coords = 0;
185 for (int d = nDims; d-- > 0; ) {
186 coords <<= nBits;
187 coords |= coord[d];
188 }
189
190 if (nBits > 1) {
191 long ndOnes = ones(nDims);
192 long nd1Ones = ndOnes >> 1;
193 int b = nDimsBits;
194 int rotation = 0;
195 long flipBit = 0;
196 long nthbits = ones(nDimsBits) / ndOnes;
197 coords = bitTranspose(nDims, nBits, coords);
198 coords ^= coords >> nDims;
199 index = 0;
200
201 do {
202 long bits = (coords >> (b -= nDims)) & ndOnes;
203 bits = rotateRight(flipBit ^ bits, rotation, nDims);
204 index <<= nDims;
205 index |= bits;
206 flipBit = 1L << rotation;
207 rotation = adjust_rotation(rotation, nDims, (int) bits);
208 } while (b > 0);
209
210 index ^= nthbits >> 1;
211 } else {
212 index = coords;
213 }
214
215 for (int d = 1; d < nDimsBits; d *= 2) {
216 index ^= index >> d;
217 }
218
219 } else {
220 index = coord[0];
221 }
222
223 return index;
224 }
225
226
227
228
229
230
231 public static void main(String[] args) {
232 int nBonds = 2;
233 int nTorsions = 4;
234 int nBits = (int) Math.ceil(FastMath.log(2, nTorsions));
235
236
237
238
239
240 BigInteger maxIndex = BigInteger.valueOf(2).pow(nBonds * nBits).subtract(BigInteger.ONE);
241 BigInteger numConfigs = BigInteger.valueOf(nTorsions).pow(nBonds);
242
243 System.out.println("Maximum index: " + maxIndex);
244 System.out.println("Number of configurations: " + numConfigs);
245
246 BigInteger index = BigInteger.ZERO;
247 int counter = 0;
248 while (index.longValue() <= maxIndex.longValue()) {
249 long[] coordinates = hilbertIndexToCoordinates(nBonds, nBits, index.longValue());
250 long convertedIndex = coordinatesToHilbertIndex(nBonds, nBits, coordinates);
251 boolean valid = true;
252
253 for (long coord : coordinates) {
254 if (coord > nTorsions - 1) {
255 valid = false;
256 break;
257 }
258 }
259
260
261 System.out.println("Hilbert index: " + counter + "; Coordinates: " + Arrays.toString(coordinates));
262 System.out.println("Converted index: " + convertedIndex);
263 counter++;
264 index = index.add(BigInteger.ONE);
265 }
266 System.out.println("Number of valid configurations: " + counter);
267 }
268 }