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