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 * @version $Id: HaltonSequenceGenerator.html 885258 2013-11-03 02:46:49Z tn $
053 * @since 3.3
054 */
055public class HaltonSequenceGenerator implements RandomVectorGenerator {
056
057    /** The first 40 primes. */
058    private static final int[] PRIMES = new int[] {
059        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
060        71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
061        149, 151, 157, 163, 167, 173
062    };
063
064    /** The optimal weights used for scrambling of the first 40 dimension. */
065    private static final int[] WEIGHTS = new int[] {
066        1, 2, 3, 3, 8, 11, 12, 14, 7, 18, 12, 13, 17, 18, 29, 14, 18, 43, 41,
067        44, 40, 30, 47, 65, 71, 28, 40, 60, 79, 89, 56, 50, 52, 61, 108, 56,
068        66, 63, 60, 66
069    };
070
071    /** Space dimension. */
072    private final int dimension;
073
074    /** The current index in the sequence. */
075    private int count = 0;
076
077    /** The base numbers for each component. */
078    private final int[] base;
079
080    /** The scrambling weights for each component. */
081    private final int[] weight;
082
083    /**
084     * Construct a new Halton sequence generator for the given space dimension.
085     *
086     * @param dimension the space dimension
087     * @throws OutOfRangeException if the space dimension is outside the allowed range of [1, 40]
088     */
089    public HaltonSequenceGenerator(final int dimension) throws OutOfRangeException {
090        this(dimension, PRIMES, WEIGHTS);
091    }
092
093    /**
094     * Construct a new Halton sequence generator with the given base numbers and weights for each dimension.
095     * The length of the bases array defines the space dimension and is required to be &gt; 0.
096     *
097     * @param dimension the space dimension
098     * @param bases the base number for each dimension, entries should be (pairwise) prime, may not be null
099     * @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            throws NullArgumentException, OutOfRangeException, DimensionMismatchException {
107
108        MathUtils.checkNotNull(bases);
109
110        if (dimension < 1 || dimension > bases.length) {
111            throw new OutOfRangeException(dimension, 1, PRIMES.length);
112        }
113
114        if (weights != null && weights.length != bases.length) {
115            throw new DimensionMismatchException(weights.length, bases.length);
116        }
117
118        this.dimension = dimension;
119        this.base = bases.clone();
120        this.weight = weights == null ? null : weights.clone();
121        count = 0;
122    }
123
124    /** {@inheritDoc} */
125    public double[] nextVector() {
126        final double[] v = new double[dimension];
127        for (int i = 0; i < dimension; i++) {
128            int index = count;
129            double f = 1.0 / base[i];
130
131            int j = 0;
132            while (index > 0) {
133                final int digit = scramble(i, j, base[i], index % base[i]);
134                v[i] += f * digit;
135                index /= base[i]; // floor( index / base )
136                f /= base[i];
137            }
138        }
139        count++;
140        return v;
141    }
142
143    /**
144     * Performs scrambling of digit {@code d_j} according to the formula:
145     * <pre>
146     *   ( weight_i * d_j ) mod base
147     * </pre>
148     * Implementations can override this method to do a different scrambling.
149     *
150     * @param i the dimension index
151     * @param j the digit index
152     * @param b the base for this dimension
153     * @param digit the j-th digit
154     * @return the scrambled digit
155     */
156    protected int scramble(final int i, final int j, final int b, final int digit) {
157        return weight != null ? (weight[i] * digit) % b : digit;
158    }
159
160    /**
161     * Skip to the i-th point in the Halton sequence.
162     * <p>
163     * This operation can be performed in O(1).
164     *
165     * @param index the index in the sequence to skip to
166     * @return the i-th point in the Halton sequence
167     * @throws NotPositiveException if index &lt; 0
168     */
169    public double[] skipTo(final int index) throws NotPositiveException {
170        count = index;
171        return nextVector();
172    }
173
174    /**
175     * Returns the index i of the next point in the Halton sequence that will be returned
176     * by calling {@link #nextVector()}.
177     *
178     * @return the index of the next point
179     */
180    public int getNextIndex() {
181        return count;
182    }
183
184}