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.math4.legacy.random;
018
019import java.util.function.Supplier;
020
021import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
022import org.apache.commons.math4.legacy.exception.NotPositiveException;
023import org.apache.commons.math4.legacy.exception.NullArgumentException;
024import org.apache.commons.math4.legacy.exception.OutOfRangeException;
025
026/**
027 * Implementation of a Halton sequence.
028 * <p>
029 * A Halton sequence is a low-discrepancy sequence generating points in the interval [0, 1] according to
030 * <pre>
031 *   H(n) = d_0 / b + d_1 / b^2 .... d_j / b^j+1
032 *
033 *   with
034 *
035 *   n = d_j * b^j-1 + ... d_1 * b + d_0 * b^0
036 * </pre>
037 * For higher dimensions, subsequent prime numbers are used as base, e.g. { 2, 3, 5 } for a Halton sequence in R^3.
038 * <p>
039 * Halton sequences are known to suffer from linear correlation for larger prime numbers, thus the individual digits
040 * are usually scrambled. This implementation already comes with support for up to 40 dimensions with optimal weight
041 * numbers from <a href="http://etd.lib.fsu.edu/theses/available/etd-07062004-140409/unrestricted/dissertation1.pdf">
042 * H. Chi: Scrambled quasirandom sequences and their applications</a>.
043 * <p>
044 * The generator supports two modes:
045 * <ul>
046 *   <li>sequential generation of points: {@link #get()}</li>
047 *   <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li>
048 * </ul>
049 *
050 * @see <a href="http://en.wikipedia.org/wiki/Halton_sequence">Halton sequence (Wikipedia)</a>
051 * @see <a href="https://lirias.kuleuven.be/bitstream/123456789/131168/1/mcm2005_bartv.pdf">
052 * On the Halton sequence and its scramblings</a>
053 * @since 3.3
054 */
055public class HaltonSequenceGenerator implements Supplier<double[]> {
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;
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) {
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        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}