001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the "math NOTICE.txt" file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package squidpony.squidmath; 018 019import squidpony.annotation.GwtIncompatible; 020 021import java.io.BufferedReader; 022import java.io.IOException; 023import java.io.InputStream; 024import java.io.InputStreamReader; 025import java.nio.charset.Charset; 026import java.util.Arrays; 027import java.util.NoSuchElementException; 028import java.util.StringTokenizer; 029 030/** 031 * Implementation of a Sobol sequence as a Quasi-Random Number Generator. 032 * <p> 033 * A Sobol sequence is a low-discrepancy sequence with the property that for all values of N, 034 * its subsequence (x1, ... xN) has a low discrepancy. It can be used to generate quasi-random 035 * points in a space S, which are equi-distributed. This is not a true random number generator, 036 * and is not even a pseudo-random number generator; the sequence generated from identical 037 * starting points with identical dimensions will be exactly the same. Calling this class' 038 * nextVector, nextIntVector, and nextLongVector methods all increment the position in the 039 * sequence, and do not use separate sequences for separate types. 040 * <p> 041 * The implementation already comes with support for up to 1000 dimensions with direction numbers 042 * calculated from <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Stephen Joe and Frances Kuo</a>. 043 * <p> 044 * The generator supports two modes: 045 * <ul> 046 * <li>sequential generation of points: {@link #nextVector()}, {@link #nextIntVector()}, 047 * {@link #nextLongVector()}, and the bounded variants on each of those</li> 048 * <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li> 049 * </ul> 050 * 051 * @see <a href="http://en.wikipedia.org/wiki/Sobol_sequence">Sobol sequence (Wikipedia)</a> 052 * @see <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Sobol sequence direction numbers</a> 053 * 054 * Created by Tommy Ettinger on 5/2/2015 based off Apache Commons Math 4. 055 */ 056@GwtIncompatible /* Because of getResourceAsStream */ 057public class SobolQRNG implements RandomnessSource { 058 059 /** The number of bits to use. */ 060 private static final int BITS = 52; 061 062 /** The scaling factor. */ 063 private static final double SCALE = Math.pow(2, BITS); 064 065 /** The maximum supported space dimension. */ 066 private static final int MAX_DIMENSION = 1000; 067 068 /** The resource containing the direction numbers. */ 069 private static final String RESOURCE_NAME = "/qrng/new-joe-kuo-6.1000.txt"; 070 071 /** Character set for file input. */ 072 private static final String FILE_CHARSET = "US-ASCII"; 073 074 private static final long serialVersionUID = -6759002780425873173L; 075 076 /** Space dimension. */ 077 private final int dimension; 078 079 /** The current index in the sequence. Starts at 1, not 0, because 0 acts differently and shouldn't be typical.*/ 080 private int count = 1; 081 082 /** The direction vector for each component. */ 083 private final long[][] direction; 084 085 /** The current state. */ 086 private final long[] x; 087 088 /** 089 * Construct a new Sobol sequence generator for the given space dimension. 090 * You should call {@link #skipTo(int)} with a fairly large number (over 1000) to ensure the results aren't 091 * too obviously non-random. If you skipTo(1), all doubles in that result will be 0.5, and if you skipTo(0), 092 * all will be 0 (this class starts at index 1 instead of 0 for that reason). This is true for all dimensions. 093 * 094 * @param dimension the space dimension 095 * @throws ArithmeticException if the space dimension is outside the allowed range of [1, 1000] 096 */ 097 public SobolQRNG(final int dimension) throws ArithmeticException { 098 if (dimension < 1 || dimension > MAX_DIMENSION) { 099 throw new ArithmeticException("Dimension " + dimension + "is outside the allowed range of [1, 1000]"); 100 } 101 102 // initialize the other dimensions with direction numbers from a resource 103 final InputStream is = getClass().getResourceAsStream(RESOURCE_NAME); 104 if (is == null) { 105 throw new ArithmeticException("Could not load embedded resource 'new-joe-kuo-6.1000'"); 106 } 107 108 this.dimension = dimension; 109 110 // init data structures 111 direction = new long[dimension][BITS + 1]; 112 x = new long[dimension]; 113 114 try { 115 initFromStream(is); 116 } catch (IOException e) { 117 // the internal resource file could not be read -> should not happen 118 throw new InternalError("Resource is unreadable, abort!\n" + e.getMessage()); 119 } catch (ArithmeticException e) { 120 // the internal resource file could not be parsed -> should not happen 121 throw new InternalError("Resource is unreadable, abort!\n" + e.getMessage()); 122 } finally { 123 try { 124 is.close(); 125 } catch (IOException e) { // NOPMD 126 // ignore 127 } 128 } 129 } 130/* 131 / ** 132 * Construct a new Sobol sequence generator for the given space dimension with 133 * direction vectors loaded from the given stream. 134 * <p> 135 * The expected format is identical to the files available from 136 * <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Stephen Joe and Frances Kuo</a>. 137 * The first line will be ignored as it is assumed to contain only the column headers. 138 * The columns are: 139 * <ul> 140 * <li>d: the dimension</li> 141 * <li>s: the degree of the primitive polynomial</li> 142 * <li>a: the number representing the coefficients</li> 143 * <li>m: the list of initial direction numbers</li> 144 * </ul> 145 * Example: 146 * <pre> 147 * d s a m_i 148 * 2 1 0 1 149 * 3 2 1 1 3 150 * </pre> 151 * <p> 152 * The input stream <i>must</i> be an ASCII text containing one valid direction vector per line. 153 * 154 * @param dimension the space dimension 155 * @param is the stream to read the direction vectors from 156 * @throws ArithmeticException if the space dimension is outside the range [1, max], where 157 * max refers to the maximum dimension found in the input stream 158 * @throws IOException if an error occurs while reading from the input stream 159 * / 160 public SobolQRNG(final int dimension, final InputStream is) 161 throws ArithmeticException, IOException { 162 163 if (dimension < 1) { 164 throw new ArithmeticException("The given dimension, " + dimension + ", is too low (must be greater than 0)"); 165 } 166 167 this.dimension = dimension; 168 169 // init data structures 170 direction = new long[dimension][BITS + 1]; 171 x = new long[dimension]; 172 173 // initialize the other dimensions with direction numbers from the stream 174 int lastDimension = initFromStream(is); 175 if (lastDimension < dimension) { 176 throw new ArithmeticException("Dimension " + dimension + "is outside the allowed range of [1, 1000]"); 177 } 178 } 179*/ 180 /** 181 * Load the direction vector for each dimension from the given stream. 182 * <p> 183 * The input stream <i>must</i> be an ASCII text containing one 184 * valid direction vector per line. 185 * 186 * @param is the input stream to read the direction vector from 187 * @return the last dimension that has been read from the input stream 188 * @throws IOException if the stream could not be read 189 * @throws ArithmeticException if the content could not be parsed successfully 190 */ 191 private int initFromStream(final InputStream is) throws ArithmeticException, IOException { 192 193 // special case: dimension 1 -> use unit initialization 194 for (int i = 1; i <= BITS; i++) { 195 direction[0][i] = 1l << (BITS - i); 196 } 197 198 final Charset charset = Charset.forName(FILE_CHARSET); 199 final BufferedReader reader = new BufferedReader(new InputStreamReader(is, charset)); 200 int dim = -1; 201 202 try { 203 // ignore first line 204 reader.readLine(); 205 206 int lineNumber = 2; 207 int index = 1; 208 String line = null; 209 while ( (line = reader.readLine()) != null) { 210 StringTokenizer st = new StringTokenizer(line, " "); 211 try { 212 dim = Integer.parseInt(st.nextToken()); 213 if (dim >= 2 && dim <= dimension) { // we have found the right dimension 214 final int s = Integer.parseInt(st.nextToken()); 215 final int a = Integer.parseInt(st.nextToken()); 216 final int[] m = new int[s + 1]; 217 for (int i = 1; i <= s; i++) { 218 m[i] = Integer.parseInt(st.nextToken()); 219 } 220 initDirectionVector(index++, a, m); 221 } 222 223 if (dim > dimension) { 224 return dim; 225 } 226 } catch (NoSuchElementException e) { 227 throw new ArithmeticException("Resource error at line " + lineNumber + ":\n " + line); 228 } catch (NumberFormatException e) { 229 throw new ArithmeticException("Resource error at line " + lineNumber + ":\n " + line); 230 } 231 lineNumber++; 232 } 233 } finally { 234 reader.close(); 235 } 236 237 return dim; 238 } 239 240 /** 241 * Calculate the direction numbers from the given polynomial. 242 * 243 * @param d the dimension, zero-based 244 * @param a the coefficients of the primitive polynomial 245 * @param m the initial direction numbers 246 */ 247 private void initDirectionVector(final int d, final int a, final int[] m) { 248 final int s = m.length - 1; 249 for (int i = 1; i <= s; i++) { 250 direction[d][i] = ((long) m[i]) << (BITS - i); 251 } 252 for (int i = s + 1; i <= BITS; i++) { 253 direction[d][i] = direction[d][i - s] ^ (direction[d][i - s] >> s); 254 for (int k = 1; k <= s - 1; k++) { 255 direction[d][i] ^= ((a >> (s - 1 - k)) & 1) * direction[d][i - k]; 256 } 257 } 258 } 259 260 /** Generate a random vector. 261 * @return a random vector as an array of double in the range [0.0, 1.0). 262 */ 263 public double[] nextVector() { 264 final double[] v = new double[dimension]; 265 if (count == 0) { 266 count++; 267 return v; 268 } 269 270 // find the index c of the rightmost 0 271 int c = 1; 272 int value = count - 1; 273 while ((value & 1) == 1) { 274 value >>= 1; 275 c++; 276 } 277 278 for (int i = 0; i < dimension; i++) { 279 x[i] ^= direction[i][c]; 280 v[i] = (double) x[i] / SCALE; 281 } 282 count++; 283 return v; 284 } 285 286 /** Generate a random vector. 287 * @param max the maximum exclusive value for the elements of the desired vector; minimum is 0 inclusive. 288 * @return a random vector as an array of double. 289 */ 290 public double[] nextVector(final double max) { 291 final double[] v = new double[dimension]; 292 if (count == 0) { 293 count++; 294 return v; 295 } 296 297 // find the index c of the rightmost 0 298 int c = 1; 299 int value = count - 1; 300 while ((value & 1) == 1) { 301 value >>= 1; 302 c++; 303 } 304 305 for (int i = 0; i < dimension; i++) { 306 x[i] ^= direction[i][c]; 307 v[i] = (double) x[i] / SCALE * max; 308 } 309 count++; 310 return v; 311 } 312 313 /** Generate a random vector. 314 * @return a random vector as an array of long (only 52 bits are actually used for the result, plus sign bit). 315 */ 316 public long[] nextLongVector() { 317 final long[] v = new long[dimension]; 318 if (count == 0) { 319 count++; 320 return v; 321 } 322 323 // find the index c of the rightmost 0 324 int c = 1 + Integer.numberOfTrailingZeros(count); 325 326 for (int i = 0; i < dimension; i++) { 327 x[i] ^= direction[i][c]; 328 v[i] = x[i]; 329 } 330 count++; 331 return v; 332 } 333 /** Generate a random vector. 334 * @param max the maximum exclusive value for the elements of the desired vector; minimum is 0 inclusive. 335 * @return a random vector as an array of long (only 52 bits are actually used for the result, plus sign bit). 336 */ 337 public long[] nextLongVector(final long max) { 338 final long[] v = new long[dimension]; 339 340 if (count == 0) { 341 count++; 342 return v; 343 } 344 345 // find the index c of the rightmost 0 346 int c = 1 + Integer.numberOfTrailingZeros(count); 347 348 for (int i = 0; i < dimension; i++) { 349 x[i] ^= direction[i][c]; 350 //suboptimal, but this isn't meant for quality of randomness, actually the opposite. 351 v[i] = (long)(x[i] / SCALE * max) % max; 352 } 353 count++; 354 return v; 355 } 356 357 /** Generate a random vector. 358 * @return a random vector as an array of int. 359 */ 360 public int[] nextIntVector() { 361 final int[] v = new int[dimension]; 362 if (count == 0) { 363 count++; 364 return v; 365 } 366 367 // find the index c of the rightmost 0 368 int c = 1 + Integer.numberOfTrailingZeros(count); 369 370 for (int i = 0; i < dimension; i++) { 371 x[i] ^= direction[i][c]; 372 v[i] = (int) (x[i] & 0x7fffffff); 373 } 374 count++; 375 return v; 376 } 377 378 /** Generate a random vector. 379 * @param max the maximum exclusive value for the elements of the desired vector; minimum is 0 inclusive. 380 * @return a random vector as an array of int. 381 */ 382 public int[] nextIntVector(final int max) { 383 final int[] v = new int[dimension]; 384 if (count == 0) { 385 count++; 386 return v; 387 } 388 389 // find the index c of the rightmost 0 390 int c = 1 + Integer.numberOfTrailingZeros(count); 391 392 for (int i = 0; i < dimension; i++) { 393 x[i] ^= direction[i][c]; 394 //suboptimal, but this isn't meant for quality of randomness, actually the opposite. 395 v[i] = (int)(x[i] / SCALE * max) % max; 396 } 397 count++; 398 return v; 399 } 400 401 /** 402 * Skip to the i-th point in the Sobol sequence. 403 * <p> 404 * This operation can be performed in O(1). 405 * If index is somehow negative, this uses its absolute value instead of throwing an exception. 406 * If index is 0, the result will always be entirely 0. 407 * You should skipTo a number greater than 1000 if you want random-seeming individual numbers in each vector. 408 * 409 * @param index the index in the sequence to skip to 410 * @return the i-th point in the Sobol sequence 411 */ 412 public double[] skipTo(final int index) { 413 if (index == 0) { 414 // reset x vector 415 Arrays.fill(x, 0); 416 } else { 417 final int i = (index > 0) ? (index - 1) : (-index - 1); 418 final long grayCode = i ^ (i >> 1); // compute the gray code of i = i XOR floor(i / 2) 419 for (int j = 0; j < dimension; j++) { 420 long result = 0; 421 for (int k = 1; k <= BITS; k++) { 422 final long shift = grayCode >> (k - 1); 423 if (shift == 0) { 424 // stop, as all remaining bits will be zero 425 break; 426 } 427 // the k-th bit of i 428 final long ik = shift & 1; 429 result ^= ik * direction[j][k]; 430 } 431 x[j] = result; 432 } 433 } 434 count = (index >= 0) ? index : -index; 435 return nextVector(); 436 } 437 438 /** 439 * Returns the index i of the next point in the Sobol sequence that will be returned 440 * by calling {@link #nextVector()}. 441 * 442 * @return the index of the next point 443 */ 444 public int getNextIndex() { 445 return count; 446 } 447 448 /** 449 * 450 * @param bits the number of bits to be returned 451 * @return 452 */ 453 @Override 454 public int next(int bits) { 455 return (int) (nextIntVector()[0] & (1L << bits) - 1); 456 } 457 458 /** 459 * Using this method, any algorithm that needs to efficiently generate more 460 * than 32 bits of random data can interface with this randomness source. 461 * <p/> 462 * Get a random long between Long.MIN_VALUE and Long.MAX_VALUE (both inclusive). 463 * 464 * @return a random long between Long.MIN_VALUE and Long.MAX_VALUE (both inclusive) 465 */ 466 @Override 467 public long nextLong() { 468 if (dimension > 1) { 469 long[] l = nextLongVector(); 470 return (l[0] << 32) ^ (l[1]); 471 } 472 return (nextLongVector()[0] << 32) ^ (nextLongVector()[0]); 473 } 474 475 /** 476 * Produces a copy of this RandomnessSource that, if next() and/or nextLong() are called on this object and the 477 * copy, both will generate the same sequence of random numbers from the point copy() was called. This just need to 478 * copy the state so it isn't shared, usually, and produce a new value with the same exact state. 479 * 480 * @return a copy of this RandomnessSource 481 */ 482 @Override 483 public RandomnessSource copy() { 484 SobolQRNG next = new SobolQRNG(dimension); 485 next.count = count; 486 return next; 487 } 488}