View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.random;
23  
24  import java.io.BufferedReader;
25  import java.io.IOException;
26  import java.io.InputStream;
27  import java.io.InputStreamReader;
28  import java.nio.charset.Charset;
29  import java.util.Arrays;
30  import java.util.NoSuchElementException;
31  import java.util.StringTokenizer;
32  
33  import org.hipparchus.exception.LocalizedCoreFormats;
34  import org.hipparchus.exception.MathIllegalArgumentException;
35  import org.hipparchus.exception.MathIllegalStateException;
36  import org.hipparchus.exception.MathRuntimeException;
37  import org.hipparchus.util.FastMath;
38  import org.hipparchus.util.MathUtils;
39  
40  /**
41   * Implementation of a Sobol sequence.
42   * <p>
43   * A Sobol sequence is a low-discrepancy sequence with the property that for all values of N,
44   * its subsequence (x1, ... xN) has a low discrepancy. It can be used to generate pseudo-random
45   * points in a space S, which are equi-distributed.
46   * <p>
47   * The implementation already comes with support for up to 1000 dimensions with direction numbers
48   * calculated from <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Stephen Joe and Frances Kuo</a>.
49   * <p>
50   * The generator supports two modes:
51   * <ul>
52   *   <li>sequential generation of points: {@link #nextVector()}</li>
53   *   <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li>
54   * </ul>
55   *
56   * @see <a href="http://en.wikipedia.org/wiki/Sobol_sequence">Sobol sequence (Wikipedia)</a>
57   * @see <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Sobol sequence direction numbers</a>
58   *
59   */
60  public class SobolSequenceGenerator implements RandomVectorGenerator {
61  
62      /** The number of bits to use. */
63      private static final int BITS = 52;
64  
65      /** The scaling factor. */
66      private static final double SCALE = FastMath.pow(2, BITS);
67  
68      /** The maximum supported space dimension. */
69      private static final int MAX_DIMENSION = 1000;
70  
71      /** The resource containing the direction numbers. */
72      private static final String RESOURCE_NAME = "/assets/org/hipparchus/random/new-joe-kuo-6.1000";
73  
74      /** Character set for file input. */
75      private static final String FILE_CHARSET = "US-ASCII";
76  
77      /** Space dimension. */
78      private final int dimension;
79  
80      /** The current index in the sequence. */
81      private int count;
82  
83      /** The direction vector for each component. */
84      private final long[][] direction;
85  
86      /** The current state. */
87      private final long[] x;
88  
89      /**
90       * Construct a new Sobol sequence generator for the given space dimension.
91       *
92       * @param dimension the space dimension
93       * @throws MathIllegalArgumentException if the space dimension is outside the allowed range of [1, 1000]
94       */
95      public SobolSequenceGenerator(final int dimension) throws MathIllegalArgumentException {
96          MathUtils.checkRangeInclusive(dimension, 1, MAX_DIMENSION);
97  
98          // initialize the other dimensions with direction numbers from a resource
99          try (InputStream is = getClass().getResourceAsStream(RESOURCE_NAME)) {
100             if (is == null) {
101                 throw MathRuntimeException.createInternalError();
102             }
103 
104             this.dimension = dimension;
105 
106             // init data structures
107             direction = new long[dimension][BITS + 1];
108             x = new long[dimension];
109 
110             initFromStream(is);
111         } catch (IOException | MathIllegalStateException e) {
112             // the internal resource file could not be parsed -> should not happen
113             throw MathRuntimeException.createInternalError(e);
114         }
115     }
116 
117     /**
118      * Construct a new Sobol sequence generator for the given space dimension with
119      * direction vectors loaded from the given stream.
120      * <p>
121      * The expected format is identical to the files available from
122      * <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Stephen Joe and Frances Kuo</a>.
123      * The first line will be ignored as it is assumed to contain only the column headers.
124      * The columns are:
125      * <ul>
126      *  <li>d: the dimension</li>
127      *  <li>s: the degree of the primitive polynomial</li>
128      *  <li>a: the number representing the coefficients</li>
129      *  <li>m: the list of initial direction numbers</li>
130      * </ul>
131      * Example:
132      * <pre>
133      * d       s       a       m_i
134      * 2       1       0       1
135      * 3       2       1       1 3
136      * </pre>
137      * <p>
138      * The input stream <i>must</i> be an ASCII text containing one valid direction vector per line.
139      *
140      * @param dimension the space dimension
141      * @param is the stream to read the direction vectors from
142      * @throws MathIllegalArgumentException if the space dimension is &lt; 1
143      * @throws MathIllegalArgumentException if the space dimension is outside the range [1, max], where
144      *   max refers to the maximum dimension found in the input stream
145      * @throws MathIllegalStateException if the content in the stream could not be parsed successfully
146      * @throws IOException if an error occurs while reading from the input stream
147      */
148     public SobolSequenceGenerator(final int dimension, final InputStream is)
149             throws MathIllegalArgumentException, MathIllegalStateException, IOException {
150 
151         if (dimension < 1) {
152             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
153                                                    dimension, 1);
154         }
155 
156         this.dimension = dimension;
157 
158         // init data structures
159         direction = new long[dimension][BITS + 1];
160         x = new long[dimension];
161 
162         // initialize the other dimensions with direction numbers from the stream
163         int lastDimension = initFromStream(is);
164         MathUtils.checkRangeInclusive(dimension, 1, lastDimension);
165     }
166 
167     /**
168      * Load the direction vector for each dimension from the given stream.
169      * <p>
170      * The input stream <i>must</i> be an ASCII text containing one
171      * valid direction vector per line.
172      *
173      * @param is the input stream to read the direction vector from
174      * @return the last dimension that has been read from the input stream
175      * @throws IOException if the stream could not be read
176      * @throws MathIllegalStateException if the content could not be parsed successfully
177      */
178     private int initFromStream(final InputStream is) throws MathIllegalStateException, IOException {
179 
180         // special case: dimension 1 -> use unit initialization
181         for (int i = 1; i <= BITS; i++) {
182             direction[0][i] = 1l << (BITS - i);
183         }
184 
185         final Charset charset = Charset.forName(FILE_CHARSET);
186         int dim = -1;
187 
188         try (BufferedReader reader = new BufferedReader(new InputStreamReader(is, charset))) {
189             // ignore first line
190             reader.readLine();
191 
192             int lineNumber = 2;
193             int index = 1;
194             for (String line = reader.readLine(); line != null; line = reader.readLine()) {
195                 StringTokenizer st = new StringTokenizer(line, " ");
196                 try {
197                     dim = Integer.parseInt(st.nextToken());
198                     if (dim >= 2 && dim <= dimension) { // we have found the right dimension
199                         final int s = Integer.parseInt(st.nextToken());
200                         final int a = Integer.parseInt(st.nextToken());
201                         final int[] m = new int[s + 1];
202                         for (int i = 1; i <= s; i++) {
203                             m[i] = Integer.parseInt(st.nextToken());
204                         }
205                         initDirectionVector(index++, a, m);
206                     }
207 
208                     if (dim > dimension) {
209                         return dim;
210                     }
211                 } catch (NoSuchElementException|NumberFormatException e) {
212                     throw new MathIllegalStateException(e, LocalizedCoreFormats.CANNOT_PARSE,
213                                                         line, lineNumber);
214                 }
215                 lineNumber++;
216             }
217         }
218 
219         return dim;
220     }
221 
222     /**
223      * Calculate the direction numbers from the given polynomial.
224      *
225      * @param d the dimension, zero-based
226      * @param a the coefficients of the primitive polynomial
227      * @param m the initial direction numbers
228      */
229     private void initDirectionVector(final int d, final int a, final int[] m) {
230         final int s = m.length - 1;
231         for (int i = 1; i <= s; i++) {
232             direction[d][i] = ((long) m[i]) << (BITS - i);
233         }
234         for (int i = s + 1; i <= BITS; i++) {
235             direction[d][i] = direction[d][i - s] ^ (direction[d][i - s] >> s);
236             for (int k = 1; k <= s - 1; k++) {
237                 direction[d][i] ^= ((a >> (s - 1 - k)) & 1) * direction[d][i - k];
238             }
239         }
240     }
241 
242     /** {@inheritDoc} */
243     @Override
244     public double[] nextVector() {
245         final double[] v = new double[dimension];
246         if (count == 0) {
247             count++;
248             return v;
249         }
250 
251         // find the index c of the rightmost 0
252         int c = 1;
253         int value = count - 1;
254         while ((value & 1) == 1) {
255             value >>= 1;
256             c++;
257         }
258 
259         for (int i = 0; i < dimension; i++) {
260             x[i] ^= direction[i][c];
261             v[i] = x[i] / SCALE;
262         }
263         count++;
264         return v;
265     }
266 
267     /**
268      * Skip to the i-th point in the Sobol sequence.
269      * <p>
270      * This operation can be performed in O(1).
271      *
272      * @param index the index in the sequence to skip to
273      * @return the i-th point in the Sobol sequence
274      * @throws MathIllegalArgumentException if index &lt; 0
275      */
276     public double[] skipTo(final int index) throws MathIllegalArgumentException {
277         if (index == 0) {
278             // reset x vector
279             Arrays.fill(x, 0);
280         } else {
281             final int i = index - 1;
282             final long grayCode = i ^ (i >> 1); // compute the gray code of i = i XOR floor(i / 2)
283             for (int j = 0; j < dimension; j++) {
284                 long result = 0;
285                 for (int k = 1; k <= BITS; k++) {
286                     final long shift = grayCode >> (k - 1);
287                     if (shift == 0) {
288                         // stop, as all remaining bits will be zero
289                         break;
290                     }
291                     // the k-th bit of i
292                     final long ik = shift & 1;
293                     result ^= ik * direction[j][k];
294                 }
295                 x[j] = result;
296             }
297         }
298         count = index;
299         return nextVector();
300     }
301 
302     /**
303      * Returns the index i of the next point in the Sobol sequence that will be returned
304      * by calling {@link #nextVector()}.
305      *
306      * @return the index of the next point
307      */
308     public int getNextIndex() {
309         return count;
310     }
311 
312 }