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