View Javadoc
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    *      https://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  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.random;
23  
24  import org.hipparchus.exception.LocalizedCoreFormats;
25  import org.hipparchus.exception.MathIllegalArgumentException;
26  import org.hipparchus.exception.NullArgumentException;
27  import org.hipparchus.util.MathUtils;
28  
29  /**
30   * Implementation of a Halton sequence.
31   * <p>
32   * A Halton sequence is a low-discrepancy sequence generating points in the interval [0, 1] according to
33   * <pre>
34   *   H(n) = d_0 / b + d_1 / b^2 .... d_j / b^j+1
35   *
36   *   with
37   *
38   *   n = d_j * b^j-1 + ... d_1 * b + d_0 * b^0
39   * </pre>
40   * For higher dimensions, subsequent prime numbers are used as base, e.g. { 2, 3, 5 } for a Halton sequence in R^3.
41   * <p>
42   * Halton sequences are known to suffer from linear correlation for larger prime numbers, thus the individual digits
43   * are usually scrambled. This implementation already comes with support for up to 40 dimensions with optimal weight
44   * numbers from <a href="http://etd.lib.fsu.edu/theses/available/etd-07062004-140409/unrestricted/dissertation1.pdf">
45   * H. Chi: Scrambled quasirandom sequences and their applications</a>.
46   * <p>
47   * The generator supports two modes:
48   * <ul>
49   *   <li>sequential generation of points: {@link #nextVector()}</li>
50   *   <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li>
51   * </ul>
52   *
53   * @see <a href="http://en.wikipedia.org/wiki/Halton_sequence">Halton sequence (Wikipedia)</a>
54   * @see <a href="https://lirias.kuleuven.be/bitstream/123456789/131168/1/mcm2005_bartv.pdf">
55   * On the Halton sequence and its scramblings</a>
56   */
57  public class HaltonSequenceGenerator implements RandomVectorGenerator {
58  
59      /** The first 40 primes. */
60      private static final int[] PRIMES = {
61          2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
62          71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
63          149, 151, 157, 163, 167, 173
64      };
65  
66      /** The optimal weights used for scrambling of the first 40 dimension. */
67      private static final int[] WEIGHTS = {
68          1, 2, 3, 3, 8, 11, 12, 14, 7, 18, 12, 13, 17, 18, 29, 14, 18, 43, 41,
69          44, 40, 30, 47, 65, 71, 28, 40, 60, 79, 89, 56, 50, 52, 61, 108, 56,
70          66, 63, 60, 66
71      };
72  
73      /** Space dimension. */
74      private final int dimension;
75  
76      /** The current index in the sequence. */
77      private int count;
78  
79      /** The base numbers for each component. */
80      private final int[] base;
81  
82      /** The scrambling weights for each component. */
83      private final int[] weight;
84  
85      /**
86       * Construct a new Halton sequence generator for the given space dimension.
87       *
88       * @param dimension the space dimension
89       * @throws MathIllegalArgumentException if the space dimension is outside the allowed range of [1, 40]
90       */
91      public HaltonSequenceGenerator(final int dimension) throws MathIllegalArgumentException {
92          this(dimension, PRIMES, WEIGHTS);
93      }
94  
95      /**
96       * Construct a new Halton sequence generator with the given base numbers and weights for each dimension.
97       * The length of the bases array defines the space dimension and is required to be &gt; 0.
98       *
99       * @param dimension the space dimension
100      * @param bases the base number for each dimension, entries should be (pairwise) prime, may not be null
101      * @param weights the weights used during scrambling, may be null in which case no scrambling will be performed
102      * @throws NullArgumentException if base is null
103      * @throws MathIllegalArgumentException if the space dimension is outside the range [1, len], where
104      *   len refers to the length of the bases array
105      * @throws MathIllegalArgumentException if weights is non-null and the length of the input arrays differ
106      */
107     public HaltonSequenceGenerator(final int dimension, final int[] bases, final int[] weights)
108             throws MathIllegalArgumentException, NullArgumentException {
109 
110         MathUtils.checkNotNull(bases);
111         MathUtils.checkRangeInclusive(dimension, 1, bases.length);
112 
113         if (weights != null && weights.length != bases.length) {
114             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
115                                                    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     @Override
126     public double[] nextVector() {
127         final double[] v = new double[dimension];
128         for (int i = 0; i < dimension; i++) {
129             int index = count;
130             double f = 1.0 / base[i];
131 
132             int j = 0;
133             while (index > 0) {
134                 final int digit = scramble(i, j, base[i], index % base[i]);
135                 v[i] += f * digit;
136                 index /= base[i]; // floor( index / base )
137                 f /= base[i];
138             }
139         }
140         count++;
141         return v;
142     }
143 
144     /**
145      * Performs scrambling of digit {@code d_j} according to the formula:
146      * <pre>
147      *   ( weight_i * d_j ) mod base
148      * </pre>
149      * Implementations can override this method to do a different scrambling.
150      *
151      * @param i the dimension index
152      * @param j the digit index
153      * @param b the base for this dimension
154      * @param digit the j-th digit
155      * @return the scrambled digit
156      */
157     protected int scramble(final int i, final int j, final int b, final int digit) {
158         return weight != null ? (weight[i] * digit) % b : digit;
159     }
160 
161     /**
162      * Skip to the i-th point in the Halton sequence.
163      * <p>
164      * This operation can be performed in O(1).
165      *
166      * @param index the index in the sequence to skip to
167      * @return the i-th point in the Halton sequence
168      * @throws MathIllegalArgumentException if index &lt; 0
169      */
170     public double[] skipTo(final int index) throws MathIllegalArgumentException {
171         count = index;
172         return nextVector();
173     }
174 
175     /**
176      * Returns the index i of the next point in the Halton sequence that will be returned
177      * by calling {@link #nextVector()}.
178      *
179      * @return the index of the next point
180      */
181     public int getNextIndex() {
182         return count;
183     }
184 
185 }