001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE 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 org.apache.commons.math3.random;
018
019import java.io.BufferedReader;
020import java.io.IOException;
021import java.io.InputStream;
022import java.io.InputStreamReader;
023import java.nio.charset.Charset;
024import java.util.Arrays;
025import java.util.NoSuchElementException;
026import java.util.StringTokenizer;
027
028import org.apache.commons.math3.exception.MathInternalError;
029import org.apache.commons.math3.exception.MathParseException;
030import org.apache.commons.math3.exception.NotPositiveException;
031import org.apache.commons.math3.exception.NotStrictlyPositiveException;
032import org.apache.commons.math3.exception.OutOfRangeException;
033import org.apache.commons.math3.util.FastMath;
034
035/**
036 * Implementation of a Sobol sequence.
037 * <p>
038 * A Sobol sequence is a low-discrepancy sequence with the property that for all values of N,
039 * its subsequence (x1, ... xN) has a low discrepancy. It can be used to generate pseudo-random
040 * points in a space S, which are equi-distributed.
041 * <p>
042 * The implementation already comes with support for up to 1000 dimensions with direction numbers
043 * calculated from <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Stephen Joe and Frances Kuo</a>.
044 * <p>
045 * The generator supports two modes:
046 * <ul>
047 *   <li>sequential generation of points: {@link #nextVector()}</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 * @version $Id: SobolSequenceGenerator.html 908881 2014-05-15 07:10:28Z luc $
055 * @since 3.3
056 */
057public class SobolSequenceGenerator implements RandomVectorGenerator {
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 = FastMath.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 = "/assets/org/apache/commons/math3/random/new-joe-kuo-6.1000";
070
071    /** Character set for file input. */
072    private static final String FILE_CHARSET = "US-ASCII";
073
074    /** Space dimension. */
075    private final int dimension;
076
077    /** The current index in the sequence. */
078    private int count = 0;
079
080    /** The direction vector for each component. */
081    private final long[][] direction;
082
083    /** The current state. */
084    private final long[] x;
085
086    /**
087     * Construct a new Sobol sequence generator for the given space dimension.
088     *
089     * @param dimension the space dimension
090     * @throws OutOfRangeException if the space dimension is outside the allowed range of [1, 1000]
091     */
092    public SobolSequenceGenerator(final int dimension) throws OutOfRangeException {
093        if (dimension < 1 || dimension > MAX_DIMENSION) {
094            throw new OutOfRangeException(dimension, 1, MAX_DIMENSION);
095        }
096
097        // initialize the other dimensions with direction numbers from a resource
098        final InputStream is = getClass().getResourceAsStream(RESOURCE_NAME);
099        if (is == null) {
100            throw new MathInternalError();
101        }
102
103        this.dimension = dimension;
104
105        // init data structures
106        direction = new long[dimension][BITS + 1];
107        x = new long[dimension];
108
109        try {
110            initFromStream(is);
111        } catch (IOException e) {
112            // the internal resource file could not be read -> should not happen
113            throw new MathInternalError();
114        } catch (MathParseException e) {
115            // the internal resource file could not be parsed -> should not happen
116            throw new MathInternalError();
117        } finally {
118            try {
119                is.close();
120            } catch (IOException e) { // NOPMD
121                // ignore
122            }
123        }
124    }
125
126    /**
127     * Construct a new Sobol sequence generator for the given space dimension with
128     * direction vectors loaded from the given stream.
129     * <p>
130     * The expected format is identical to the files available from
131     * <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Stephen Joe and Frances Kuo</a>.
132     * The first line will be ignored as it is assumed to contain only the column headers.
133     * The columns are:
134     * <ul>
135     *  <li>d: the dimension</li>
136     *  <li>s: the degree of the primitive polynomial</li>
137     *  <li>a: the number representing the coefficients</li>
138     *  <li>m: the list of initial direction numbers</li>
139     * </ul>
140     * Example:
141     * <pre>
142     * d       s       a       m_i
143     * 2       1       0       1
144     * 3       2       1       1 3
145     * </pre>
146     * <p>
147     * The input stream <i>must</i> be an ASCII text containing one valid direction vector per line.
148     *
149     * @param dimension the space dimension
150     * @param is the stream to read the direction vectors from
151     * @throws NotStrictlyPositiveException if the space dimension is &lt; 1
152     * @throws OutOfRangeException if the space dimension is outside the range [1, max], where
153     *   max refers to the maximum dimension found in the input stream
154     * @throws MathParseException if the content in the stream could not be parsed successfully
155     * @throws IOException if an error occurs while reading from the input stream
156     */
157    public SobolSequenceGenerator(final int dimension, final InputStream is)
158            throws NotStrictlyPositiveException, MathParseException, IOException {
159
160        if (dimension < 1) {
161            throw new NotStrictlyPositiveException(dimension);
162        }
163
164        this.dimension = dimension;
165
166        // init data structures
167        direction = new long[dimension][BITS + 1];
168        x = new long[dimension];
169
170        // initialize the other dimensions with direction numbers from the stream
171        int lastDimension = initFromStream(is);
172        if (lastDimension < dimension) {
173            throw new OutOfRangeException(dimension, 1, lastDimension);
174        }
175    }
176
177    /**
178     * Load the direction vector for each dimension from the given stream.
179     * <p>
180     * The input stream <i>must</i> be an ASCII text containing one
181     * valid direction vector per line.
182     *
183     * @param is the input stream to read the direction vector from
184     * @return the last dimension that has been read from the input stream
185     * @throws IOException if the stream could not be read
186     * @throws MathParseException if the content could not be parsed successfully
187     */
188    private int initFromStream(final InputStream is) throws MathParseException, IOException {
189
190        // special case: dimension 1 -> use unit initialization
191        for (int i = 1; i <= BITS; i++) {
192            direction[0][i] = 1l << (BITS - i);
193        }
194
195        final Charset charset = Charset.forName(FILE_CHARSET);
196        final BufferedReader reader = new BufferedReader(new InputStreamReader(is, charset));
197        int dim = -1;
198
199        try {
200            // ignore first line
201            reader.readLine();
202
203            int lineNumber = 2;
204            int index = 1;
205            String line = null;
206            while ( (line = reader.readLine()) != null) {
207                StringTokenizer st = new StringTokenizer(line, " ");
208                try {
209                    dim = Integer.parseInt(st.nextToken());
210                    if (dim >= 2 && dim <= dimension) { // we have found the right dimension
211                        final int s = Integer.parseInt(st.nextToken());
212                        final int a = Integer.parseInt(st.nextToken());
213                        final int[] m = new int[s + 1];
214                        for (int i = 1; i <= s; i++) {
215                            m[i] = Integer.parseInt(st.nextToken());
216                        }
217                        initDirectionVector(index++, a, m);
218                    }
219
220                    if (dim > dimension) {
221                        return dim;
222                    }
223                } catch (NoSuchElementException e) {
224                    throw new MathParseException(line, lineNumber);
225                } catch (NumberFormatException e) {
226                    throw new MathParseException(line, lineNumber);
227                }
228                lineNumber++;
229            }
230        } finally {
231            reader.close();
232        }
233
234        return dim;
235    }
236
237    /**
238     * Calculate the direction numbers from the given polynomial.
239     *
240     * @param d the dimension, zero-based
241     * @param a the coefficients of the primitive polynomial
242     * @param m the initial direction numbers
243     */
244    private void initDirectionVector(final int d, final int a, final int[] m) {
245        final int s = m.length - 1;
246        for (int i = 1; i <= s; i++) {
247            direction[d][i] = ((long) m[i]) << (BITS - i);
248        }
249        for (int i = s + 1; i <= BITS; i++) {
250            direction[d][i] = direction[d][i - s] ^ (direction[d][i - s] >> s);
251            for (int k = 1; k <= s - 1; k++) {
252                direction[d][i] ^= ((a >> (s - 1 - k)) & 1) * direction[d][i - k];
253            }
254        }
255    }
256
257    /** {@inheritDoc} */
258    public double[] nextVector() {
259        final double[] v = new double[dimension];
260        if (count == 0) {
261            count++;
262            return v;
263        }
264
265        // find the index c of the rightmost 0
266        int c = 1;
267        int value = count - 1;
268        while ((value & 1) == 1) {
269            value >>= 1;
270            c++;
271        }
272
273        for (int i = 0; i < dimension; i++) {
274            x[i] ^= direction[i][c];
275            v[i] = (double) x[i] / SCALE;
276        }
277        count++;
278        return v;
279    }
280
281    /**
282     * Skip to the i-th point in the Sobol sequence.
283     * <p>
284     * This operation can be performed in O(1).
285     *
286     * @param index the index in the sequence to skip to
287     * @return the i-th point in the Sobol sequence
288     * @throws NotPositiveException if index &lt; 0
289     */
290    public double[] skipTo(final int index) throws NotPositiveException {
291        if (index == 0) {
292            // reset x vector
293            Arrays.fill(x, 0);
294        } else {
295            final int i = index - 1;
296            final long grayCode = i ^ (i >> 1); // compute the gray code of i = i XOR floor(i / 2)
297            for (int j = 0; j < dimension; j++) {
298                long result = 0;
299                for (int k = 1; k <= BITS; k++) {
300                    final long shift = grayCode >> (k - 1);
301                    if (shift == 0) {
302                        // stop, as all remaining bits will be zero
303                        break;
304                    }
305                    // the k-th bit of i
306                    final long ik = shift & 1;
307                    result ^= ik * direction[j][k];
308                }
309                x[j] = result;
310            }
311        }
312        count = index;
313        return nextVector();
314    }
315
316    /**
317     * Returns the index i of the next point in the Sobol sequence that will be returned
318     * by calling {@link #nextVector()}.
319     *
320     * @return the index of the next point
321     */
322    public int getNextIndex() {
323        return count;
324    }
325
326}