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 org.apache.commons.math3.exception.DimensionMismatchException; 020import org.apache.commons.math3.exception.NotPositiveException; 021import org.apache.commons.math3.exception.NullArgumentException; 022import org.apache.commons.math3.exception.OutOfRangeException; 023import org.apache.commons.math3.util.MathUtils; 024 025/** 026 * Implementation of a Halton sequence. 027 * <p> 028 * A Halton sequence is a low-discrepancy sequence generating points in the interval [0, 1] according to 029 * <pre> 030 * H(n) = d_0 / b + d_1 / b^2 .... d_j / b^j+1 031 * 032 * with 033 * 034 * n = d_j * b^j-1 + ... d_1 * b + d_0 * b^0 035 * </pre> 036 * For higher dimensions, subsequent prime numbers are used as base, e.g. { 2, 3, 5 } for a Halton sequence in R^3. 037 * <p> 038 * Halton sequences are known to suffer from linear correlation for larger prime numbers, thus the individual digits 039 * are usually scrambled. This implementation already comes with support for up to 40 dimensions with optimal weight 040 * numbers from <a href="http://etd.lib.fsu.edu/theses/available/etd-07062004-140409/unrestricted/dissertation1.pdf"> 041 * H. Chi: Scrambled quasirandom sequences and their applications</a>. 042 * <p> 043 * The generator supports two modes: 044 * <ul> 045 * <li>sequential generation of points: {@link #nextVector()}</li> 046 * <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li> 047 * </ul> 048 * 049 * @see <a href="http://en.wikipedia.org/wiki/Halton_sequence">Halton sequence (Wikipedia)</a> 050 * @see <a href="https://lirias.kuleuven.be/bitstream/123456789/131168/1/mcm2005_bartv.pdf"> 051 * On the Halton sequence and its scramblings</a> 052 * @since 3.3 053 */ 054public class HaltonSequenceGenerator implements RandomVectorGenerator { 055 056 /** The first 40 primes. */ 057 private static final int[] PRIMES = new int[] { 058 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 059 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 060 149, 151, 157, 163, 167, 173 061 }; 062 063 /** The optimal weights used for scrambling of the first 40 dimension. */ 064 private static final int[] WEIGHTS = new int[] { 065 1, 2, 3, 3, 8, 11, 12, 14, 7, 18, 12, 13, 17, 18, 29, 14, 18, 43, 41, 066 44, 40, 30, 47, 65, 71, 28, 40, 60, 79, 89, 56, 50, 52, 61, 108, 56, 067 66, 63, 60, 66 068 }; 069 070 /** Space dimension. */ 071 private final int dimension; 072 073 /** The current index in the sequence. */ 074 private int count = 0; 075 076 /** The base numbers for each component. */ 077 private final int[] base; 078 079 /** The scrambling weights for each component. */ 080 private final int[] weight; 081 082 /** 083 * Construct a new Halton sequence generator for the given space dimension. 084 * 085 * @param dimension the space dimension 086 * @throws OutOfRangeException if the space dimension is outside the allowed range of [1, 40] 087 */ 088 public HaltonSequenceGenerator(final int dimension) throws OutOfRangeException { 089 this(dimension, PRIMES, WEIGHTS); 090 } 091 092 /** 093 * Construct a new Halton sequence generator with the given base numbers and weights for each dimension. 094 * The length of the bases array defines the space dimension and is required to be > 0. 095 * 096 * @param dimension the space dimension 097 * @param bases the base number for each dimension, entries should be (pairwise) prime, may not be null 098 * @param weights the weights used during scrambling, may be null in which case no scrambling will be performed 099 * @throws NullArgumentException if base is null 100 * @throws OutOfRangeException if the space dimension is outside the range [1, len], where 101 * len refers to the length of the bases array 102 * @throws DimensionMismatchException if weights is non-null and the length of the input arrays differ 103 */ 104 public HaltonSequenceGenerator(final int dimension, final int[] bases, final int[] weights) 105 throws NullArgumentException, OutOfRangeException, DimensionMismatchException { 106 107 MathUtils.checkNotNull(bases); 108 109 if (dimension < 1 || dimension > bases.length) { 110 throw new OutOfRangeException(dimension, 1, PRIMES.length); 111 } 112 113 if (weights != null && weights.length != bases.length) { 114 throw new DimensionMismatchException(weights.length, bases.length); 115 } 116 117 this.dimension = dimension; 118 this.base = bases.clone(); 119 this.weight = weights == null ? null : weights.clone(); 120 count = 0; 121 } 122 123 /** {@inheritDoc} */ 124 public double[] nextVector() { 125 final double[] v = new double[dimension]; 126 for (int i = 0; i < dimension; i++) { 127 int index = count; 128 double f = 1.0 / base[i]; 129 130 int j = 0; 131 while (index > 0) { 132 final int digit = scramble(i, j, base[i], index % base[i]); 133 v[i] += f * digit; 134 index /= base[i]; // floor( index / base ) 135 f /= base[i]; 136 } 137 } 138 count++; 139 return v; 140 } 141 142 /** 143 * Performs scrambling of digit {@code d_j} according to the formula: 144 * <pre> 145 * ( weight_i * d_j ) mod base 146 * </pre> 147 * Implementations can override this method to do a different scrambling. 148 * 149 * @param i the dimension index 150 * @param j the digit index 151 * @param b the base for this dimension 152 * @param digit the j-th digit 153 * @return the scrambled digit 154 */ 155 protected int scramble(final int i, final int j, final int b, final int digit) { 156 return weight != null ? (weight[i] * digit) % b : digit; 157 } 158 159 /** 160 * Skip to the i-th point in the Halton sequence. 161 * <p> 162 * This operation can be performed in O(1). 163 * 164 * @param index the index in the sequence to skip to 165 * @return the i-th point in the Halton sequence 166 * @throws NotPositiveException if index < 0 167 */ 168 public double[] skipTo(final int index) throws NotPositiveException { 169 count = index; 170 return nextVector(); 171 } 172 173 /** 174 * Returns the index i of the next point in the Halton sequence that will be returned 175 * by calling {@link #nextVector()}. 176 * 177 * @return the index of the next point 178 */ 179 public int getNextIndex() { 180 return count; 181 } 182 183}