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 | * http://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 | package agents.anac.y2019.harddealer.math3.random;
|
---|
18 |
|
---|
19 | import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
|
---|
20 | import agents.anac.y2019.harddealer.math3.exception.NotPositiveException;
|
---|
21 | import agents.anac.y2019.harddealer.math3.exception.NullArgumentException;
|
---|
22 | import agents.anac.y2019.harddealer.math3.exception.OutOfRangeException;
|
---|
23 | import agents.anac.y2019.harddealer.math3.util.MathUtils;
|
---|
24 |
|
---|
25 | /**
|
---|
26 | * Implementation of a Halton sequence.
|
---|
27 | * <p>
|
---|
28 | * A Halton sequence is a low-discrepancy sequence generating points in the interval [0, 1] according to
|
---|
29 | * <pre>
|
---|
30 | * H(n) = d_0 / b + d_1 / b^2 .... d_j / b^j+1
|
---|
31 | *
|
---|
32 | * with
|
---|
33 | *
|
---|
34 | * n = d_j * b^j-1 + ... d_1 * b + d_0 * b^0
|
---|
35 | * </pre>
|
---|
36 | * For higher dimensions, subsequent prime numbers are used as base, e.g. { 2, 3, 5 } for a Halton sequence in R^3.
|
---|
37 | * <p>
|
---|
38 | * Halton sequences are known to suffer from linear correlation for larger prime numbers, thus the individual digits
|
---|
39 | * are usually scrambled. This implementation already comes with support for up to 40 dimensions with optimal weight
|
---|
40 | * numbers from <a href="http://etd.lib.fsu.edu/theses/available/etd-07062004-140409/unrestricted/dissertation1.pdf">
|
---|
41 | * H. Chi: Scrambled quasirandom sequences and their applications</a>.
|
---|
42 | * <p>
|
---|
43 | * The generator supports two modes:
|
---|
44 | * <ul>
|
---|
45 | * <li>sequential generation of points: {@link #nextVector()}</li>
|
---|
46 | * <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li>
|
---|
47 | * </ul>
|
---|
48 | *
|
---|
49 | * @see <a href="http://en.wikipedia.org/wiki/Halton_sequence">Halton sequence (Wikipedia)</a>
|
---|
50 | * @see <a href="https://lirias.kuleuven.be/bitstream/123456789/131168/1/mcm2005_bartv.pdf">
|
---|
51 | * On the Halton sequence and its scramblings</a>
|
---|
52 | * @since 3.3
|
---|
53 | */
|
---|
54 | public class HaltonSequenceGenerator implements RandomVectorGenerator {
|
---|
55 |
|
---|
56 | /** The first 40 primes. */
|
---|
57 | private static final int[] PRIMES = new int[] {
|
---|
58 | 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
|
---|
59 | 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
|
---|
60 | 149, 151, 157, 163, 167, 173
|
---|
61 | };
|
---|
62 |
|
---|
63 | /** The optimal weights used for scrambling of the first 40 dimension. */
|
---|
64 | private static final int[] WEIGHTS = new int[] {
|
---|
65 | 1, 2, 3, 3, 8, 11, 12, 14, 7, 18, 12, 13, 17, 18, 29, 14, 18, 43, 41,
|
---|
66 | 44, 40, 30, 47, 65, 71, 28, 40, 60, 79, 89, 56, 50, 52, 61, 108, 56,
|
---|
67 | 66, 63, 60, 66
|
---|
68 | };
|
---|
69 |
|
---|
70 | /** Space dimension. */
|
---|
71 | private final int dimension;
|
---|
72 |
|
---|
73 | /** The current index in the sequence. */
|
---|
74 | private int count = 0;
|
---|
75 |
|
---|
76 | /** The base numbers for each component. */
|
---|
77 | private final int[] base;
|
---|
78 |
|
---|
79 | /** The scrambling weights for each component. */
|
---|
80 | private final int[] weight;
|
---|
81 |
|
---|
82 | /**
|
---|
83 | * Construct a new Halton sequence generator for the given space dimension.
|
---|
84 | *
|
---|
85 | * @param dimension the space dimension
|
---|
86 | * @throws OutOfRangeException if the space dimension is outside the allowed range of [1, 40]
|
---|
87 | */
|
---|
88 | public HaltonSequenceGenerator(final int dimension) throws OutOfRangeException {
|
---|
89 | this(dimension, PRIMES, WEIGHTS);
|
---|
90 | }
|
---|
91 |
|
---|
92 | /**
|
---|
93 | * Construct a new Halton sequence generator with the given base numbers and weights for each dimension.
|
---|
94 | * The length of the bases array defines the space dimension and is required to be > 0.
|
---|
95 | *
|
---|
96 | * @param dimension the space dimension
|
---|
97 | * @param bases the base number for each dimension, entries should be (pairwise) prime, may not be null
|
---|
98 | * @param weights the weights used during scrambling, may be null in which case no scrambling will be performed
|
---|
99 | * @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 < 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 | }
|
---|