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}