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 * http://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 package org.apache.commons.math4.legacy.random; 18 19 import java.util.function.Supplier; 20 21 import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 22 import org.apache.commons.math4.legacy.exception.NotPositiveException; 23 import org.apache.commons.math4.legacy.exception.NullArgumentException; 24 import org.apache.commons.math4.legacy.exception.OutOfRangeException; 25 26 /** 27 * Implementation of a Halton sequence. 28 * <p> 29 * A Halton sequence is a low-discrepancy sequence generating points in the interval [0, 1] according to 30 * <pre> 31 * H(n) = d_0 / b + d_1 / b^2 .... d_j / b^j+1 32 * 33 * with 34 * 35 * n = d_j * b^j-1 + ... d_1 * b + d_0 * b^0 36 * </pre> 37 * For higher dimensions, subsequent prime numbers are used as base, e.g. { 2, 3, 5 } for a Halton sequence in R^3. 38 * <p> 39 * Halton sequences are known to suffer from linear correlation for larger prime numbers, thus the individual digits 40 * are usually scrambled. This implementation already comes with support for up to 40 dimensions with optimal weight 41 * numbers from <a href="http://etd.lib.fsu.edu/theses/available/etd-07062004-140409/unrestricted/dissertation1.pdf"> 42 * H. Chi: Scrambled quasirandom sequences and their applications</a>. 43 * <p> 44 * The generator supports two modes: 45 * <ul> 46 * <li>sequential generation of points: {@link #get()}</li> 47 * <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li> 48 * </ul> 49 * 50 * @see <a href="http://en.wikipedia.org/wiki/Halton_sequence">Halton sequence (Wikipedia)</a> 51 * @see <a href="https://lirias.kuleuven.be/bitstream/123456789/131168/1/mcm2005_bartv.pdf"> 52 * On the Halton sequence and its scramblings</a> 53 * @since 3.3 54 */ 55 public class HaltonSequenceGenerator implements Supplier<double[]> { 56 57 /** The first 40 primes. */ 58 private static final int[] PRIMES = new int[] { 59 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 60 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 61 149, 151, 157, 163, 167, 173 62 }; 63 64 /** The optimal weights used for scrambling of the first 40 dimension. */ 65 private static final int[] WEIGHTS = new int[] { 66 1, 2, 3, 3, 8, 11, 12, 14, 7, 18, 12, 13, 17, 18, 29, 14, 18, 43, 41, 67 44, 40, 30, 47, 65, 71, 28, 40, 60, 79, 89, 56, 50, 52, 61, 108, 56, 68 66, 63, 60, 66 69 }; 70 71 /** Space dimension. */ 72 private final int dimension; 73 74 /** The current index in the sequence. */ 75 private int count; 76 77 /** The base numbers for each component. */ 78 private final int[] base; 79 80 /** The scrambling weights for each component. */ 81 private final int[] weight; 82 83 /** 84 * Construct a new Halton sequence generator for the given space dimension. 85 * 86 * @param dimension the space dimension 87 * @throws OutOfRangeException if the space dimension is outside the allowed range of [1, 40] 88 */ 89 public HaltonSequenceGenerator(final int dimension) { 90 this(dimension, PRIMES, WEIGHTS); 91 } 92 93 /** 94 * Construct a new Halton sequence generator with the given base numbers and weights for each dimension. 95 * The length of the bases array defines the space dimension and is required to be > 0. 96 * 97 * @param dimension the space dimension 98 * @param bases the base number for each dimension, entries should be (pairwise) prime, may not be null 99 * @param weights the weights used during scrambling, may be null in which case no scrambling will be performed 100 * @throws NullArgumentException if base is null 101 * @throws OutOfRangeException if the space dimension is outside the range [1, len], where 102 * len refers to the length of the bases array 103 * @throws DimensionMismatchException if weights is non-null and the length of the input arrays differ 104 */ 105 public HaltonSequenceGenerator(final int dimension, final int[] bases, final int[] weights) { 106 NullArgumentException.check(bases); 107 108 if (dimension < 1 || dimension > bases.length) { 109 throw new OutOfRangeException(dimension, 1, PRIMES.length); 110 } 111 112 if (weights != null && weights.length != bases.length) { 113 throw new DimensionMismatchException(weights.length, bases.length); 114 } 115 116 this.dimension = dimension; 117 this.base = bases.clone(); 118 this.weight = weights == null ? null : weights.clone(); 119 count = 0; 120 } 121 122 /** {@inheritDoc} */ 123 @Override 124 public double[] get() { 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 {@code index < 0}. 167 */ 168 public double[] skipTo(final int index) { 169 if (index < 0) { 170 throw new NotPositiveException(index); 171 } 172 173 count = index; 174 return get(); 175 } 176 177 /** 178 * Returns the index i of the next point in the Halton sequence that will be returned 179 * by calling {@link #get()}. 180 * 181 * @return the index of the next point 182 */ 183 public int getNextIndex() { 184 return count; 185 } 186 }