| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.commons.math4.legacy.random; |
| |
| import java.util.function.Supplier; |
| |
| import org.apache.commons.math4.legacy.exception.DimensionMismatchException; |
| import org.apache.commons.math4.legacy.exception.NullArgumentException; |
| import org.apache.commons.math4.legacy.exception.OutOfRangeException; |
| |
| /** |
| * Implementation of a Halton sequence. |
| * <p> |
| * A Halton sequence is a low-discrepancy sequence generating points in the interval [0, 1] according to |
| * <pre> |
| * H(n) = d_0 / b + d_1 / b^2 .... d_j / b^j+1 |
| * |
| * with |
| * |
| * n = d_j * b^j-1 + ... d_1 * b + d_0 * b^0 |
| * </pre> |
| * For higher dimensions, subsequent prime numbers are used as base, e.g. { 2, 3, 5 } for a Halton sequence in R^3. |
| * <p> |
| * Halton sequences are known to suffer from linear correlation for larger prime numbers, thus the individual digits |
| * are usually scrambled. This implementation already comes with support for up to 40 dimensions with optimal weight |
| * numbers from <a href="http://etd.lib.fsu.edu/theses/available/etd-07062004-140409/unrestricted/dissertation1.pdf"> |
| * H. Chi: Scrambled quasirandom sequences and their applications</a>. |
| * <p> |
| * The generator supports two modes: |
| * <ul> |
| * <li>sequential generation of points: {@link #get()}</li> |
| * <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li> |
| * </ul> |
| * |
| * @see <a href="http://en.wikipedia.org/wiki/Halton_sequence">Halton sequence (Wikipedia)</a> |
| * @see <a href="https://lirias.kuleuven.be/bitstream/123456789/131168/1/mcm2005_bartv.pdf"> |
| * On the Halton sequence and its scramblings</a> |
| * @since 3.3 |
| */ |
| public class HaltonSequenceGenerator implements Supplier<double[]> { |
| |
| /** The first 40 primes. */ |
| private static final int[] PRIMES = new int[] { |
| 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, |
| 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, |
| 149, 151, 157, 163, 167, 173 |
| }; |
| |
| /** The optimal weights used for scrambling of the first 40 dimension. */ |
| private static final int[] WEIGHTS = new int[] { |
| 1, 2, 3, 3, 8, 11, 12, 14, 7, 18, 12, 13, 17, 18, 29, 14, 18, 43, 41, |
| 44, 40, 30, 47, 65, 71, 28, 40, 60, 79, 89, 56, 50, 52, 61, 108, 56, |
| 66, 63, 60, 66 |
| }; |
| |
| /** Space dimension. */ |
| private final int dimension; |
| |
| /** The current index in the sequence. */ |
| private int count; |
| |
| /** The base numbers for each component. */ |
| private final int[] base; |
| |
| /** The scrambling weights for each component. */ |
| private final int[] weight; |
| |
| /** |
| * Construct a new Halton sequence generator for the given space dimension. |
| * |
| * @param dimension the space dimension |
| * @throws OutOfRangeException if the space dimension is outside the allowed range of [1, 40] |
| */ |
| public HaltonSequenceGenerator(final int dimension) { |
| this(dimension, PRIMES, WEIGHTS); |
| } |
| |
| /** |
| * Construct a new Halton sequence generator with the given base numbers and weights for each dimension. |
| * The length of the bases array defines the space dimension and is required to be > 0. |
| * |
| * @param dimension the space dimension |
| * @param bases the base number for each dimension, entries should be (pairwise) prime, may not be null |
| * @param weights the weights used during scrambling, may be null in which case no scrambling will be performed |
| * @throws NullArgumentException if base is null |
| * @throws OutOfRangeException if the space dimension is outside the range [1, len], where |
| * len refers to the length of the bases array |
| * @throws DimensionMismatchException if weights is non-null and the length of the input arrays differ |
| */ |
| public HaltonSequenceGenerator(final int dimension, final int[] bases, final int[] weights) { |
| NullArgumentException.check(bases); |
| |
| if (dimension < 1 || dimension > bases.length) { |
| throw new OutOfRangeException(dimension, 1, PRIMES.length); |
| } |
| |
| if (weights != null && weights.length != bases.length) { |
| throw new DimensionMismatchException(weights.length, bases.length); |
| } |
| |
| this.dimension = dimension; |
| this.base = bases.clone(); |
| this.weight = weights == null ? null : weights.clone(); |
| count = 0; |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public double[] get() { |
| final double[] v = new double[dimension]; |
| for (int i = 0; i < dimension; i++) { |
| int index = count; |
| double f = 1.0 / base[i]; |
| |
| int j = 0; |
| while (index > 0) { |
| final int digit = scramble(i, j, base[i], index % base[i]); |
| v[i] += f * digit; |
| index /= base[i]; // floor( index / base ) |
| f /= base[i]; |
| } |
| } |
| count++; |
| return v; |
| } |
| |
| /** |
| * Performs scrambling of digit {@code d_j} according to the formula: |
| * <pre> |
| * ( weight_i * d_j ) mod base |
| * </pre> |
| * Implementations can override this method to do a different scrambling. |
| * |
| * @param i the dimension index |
| * @param j the digit index |
| * @param b the base for this dimension |
| * @param digit the j-th digit |
| * @return the scrambled digit |
| */ |
| protected int scramble(final int i, final int j, final int b, final int digit) { |
| return weight != null ? (weight[i] * digit) % b : digit; |
| } |
| |
| /** |
| * Skip to the i-th point in the Halton sequence. |
| * <p> |
| * This operation can be performed in O(1). |
| * |
| * @param index the index in the sequence to skip to |
| * @return the i-th point in the Halton sequence |
| * @throws org.apache.commons.math4.legacy.exception.NotPositiveException NotPositiveException if index < 0 |
| */ |
| public double[] skipTo(final int index) { |
| count = index; |
| return get(); |
| } |
| |
| /** |
| * Returns the index i of the next point in the Halton sequence that will be returned |
| * by calling {@link #get()}. |
| * |
| * @return the index of the next point |
| */ |
| public int getNextIndex() { |
| return count; |
| } |
| } |