Source code

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