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 &gt; 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 &lt; 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}