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
18 package org.apache.commons.rng.sampling.distribution;
19
20 import org.apache.commons.rng.UniformRandomProvider;
21 import org.apache.commons.rng.sampling.SharedStateSampler;
22
23 /**
24 * Functions used by some of the samplers.
25 * This class is not part of the public API, as it would be
26 * better to group these utilities in a dedicated component.
27 */
28 final class InternalUtils { // Class is package-private on purpose; do not make it public.
29 /** All long-representable factorials. */
30 private static final long[] FACTORIALS = {
31 1L, 1L, 2L,
32 6L, 24L, 120L,
33 720L, 5040L, 40320L,
34 362880L, 3628800L, 39916800L,
35 479001600L, 6227020800L, 87178291200L,
36 1307674368000L, 20922789888000L, 355687428096000L,
37 6402373705728000L, 121645100408832000L, 2432902008176640000L };
38
39 /** The first array index with a non-zero log factorial. */
40 private static final int BEGIN_LOG_FACTORIALS = 2;
41
42 /**
43 * The multiplier to convert the least significant 53-bits of a {@code long} to a {@code double}.
44 * Taken from org.apache.commons.rng.core.util.NumberFactory.
45 */
46 private static final double DOUBLE_MULTIPLIER = 0x1.0p-53d;
47
48 /** Utility class. */
49 private InternalUtils() {}
50
51 /**
52 * @param n Argument.
53 * @return {@code n!}
54 * @throws IndexOutOfBoundsException if the result is too large to be represented
55 * by a {@code long} (i.e. if {@code n > 20}), or {@code n} is negative.
56 */
57 static long factorial(int n) {
58 return FACTORIALS[n];
59 }
60
61 /**
62 * Validate the probabilities sum to a finite positive number.
63 *
64 * @param probabilities the probabilities
65 * @return the sum
66 * @throws IllegalArgumentException if {@code probabilities} is null or empty, a
67 * probability is negative, infinite or {@code NaN}, or the sum of all
68 * probabilities is not strictly positive.
69 */
70 static double validateProbabilities(double[] probabilities) {
71 if (probabilities == null || probabilities.length == 0) {
72 throw new IllegalArgumentException("Probabilities must not be empty.");
73 }
74
75 double sumProb = 0;
76 for (final double prob : probabilities) {
77 validateProbability(prob);
78 sumProb += prob;
79 }
80
81 if (Double.isInfinite(sumProb) || sumProb <= 0) {
82 throw new IllegalArgumentException("Invalid sum of probabilities: " + sumProb);
83 }
84 return sumProb;
85 }
86
87 /**
88 * Validate the probability is a finite positive number.
89 *
90 * @param probability Probability.
91 * @throws IllegalArgumentException if {@code probability} is negative, infinite or {@code NaN}.
92 */
93 static void validateProbability(double probability) {
94 if (probability < 0 ||
95 Double.isInfinite(probability) ||
96 Double.isNaN(probability)) {
97 throw new IllegalArgumentException("Invalid probability: " +
98 probability);
99 }
100 }
101
102 /**
103 * Create a new instance of the given sampler using
104 * {@link SharedStateSampler#withUniformRandomProvider(UniformRandomProvider)}.
105 *
106 * @param sampler Source sampler.
107 * @param rng Generator of uniformly distributed random numbers.
108 * @return the new sampler
109 * @throws UnsupportedOperationException if the underlying sampler is not a
110 * {@link SharedStateSampler} or does not return a {@link NormalizedGaussianSampler} when
111 * sharing state.
112 */
113 static NormalizedGaussianSampler newNormalizedGaussianSampler(
114 NormalizedGaussianSampler sampler,
115 UniformRandomProvider rng) {
116 if (!(sampler instanceof SharedStateSampler<?>)) {
117 throw new UnsupportedOperationException("The underlying sampler cannot share state");
118 }
119 final Object newSampler = ((SharedStateSampler<?>) sampler).withUniformRandomProvider(rng);
120 if (!(newSampler instanceof NormalizedGaussianSampler)) {
121 throw new UnsupportedOperationException(
122 "The underlying sampler did not create a normalized Gaussian sampler");
123 }
124 return (NormalizedGaussianSampler) newSampler;
125 }
126
127 /**
128 * Creates a {@code double} in the interval {@code (0, 1]} from a {@code long} value.
129 *
130 * @param v Number.
131 * @return a {@code double} value in the interval {@code (0, 1]}.
132 */
133 static double makeNonZeroDouble(long v) {
134 // This matches the method in o.a.c.rng.core.util.NumberFactory.makeDouble(long)
135 // but shifts the range from [0, 1) to (0, 1].
136 return ((v >>> 11) + 1L) * DOUBLE_MULTIPLIER;
137 }
138
139 /**
140 * Class for computing the natural logarithm of the factorial of {@code n}.
141 * It allows to allocate a cache of precomputed values.
142 * In case of cache miss, computation is performed by a call to
143 * {@link InternalGamma#logGamma(double)}.
144 */
145 public static final class FactorialLog {
146 /**
147 * Precomputed values of the function:
148 * {@code LOG_FACTORIALS[i] = log(i!)}.
149 */
150 private final double[] logFactorials;
151
152 /**
153 * Creates an instance, reusing the already computed values if available.
154 *
155 * @param numValues Number of values of the function to compute.
156 * @param cache Existing cache.
157 * @throws NegativeArraySizeException if {@code numValues < 0}.
158 */
159 private FactorialLog(int numValues,
160 double[] cache) {
161 logFactorials = new double[numValues];
162
163 int endCopy;
164 if (cache != null && cache.length > BEGIN_LOG_FACTORIALS) {
165 // Copy available values.
166 endCopy = Math.min(cache.length, numValues);
167 System.arraycopy(cache, BEGIN_LOG_FACTORIALS, logFactorials, BEGIN_LOG_FACTORIALS,
168 endCopy - BEGIN_LOG_FACTORIALS);
169 } else {
170 // All values to be computed
171 endCopy = BEGIN_LOG_FACTORIALS;
172 }
173
174 // Compute remaining values.
175 for (int i = endCopy; i < numValues; i++) {
176 if (i < FACTORIALS.length) {
177 logFactorials[i] = Math.log(FACTORIALS[i]);
178 } else {
179 logFactorials[i] = logFactorials[i - 1] + Math.log(i);
180 }
181 }
182 }
183
184 /**
185 * Creates an instance with no precomputed values.
186 *
187 * @return an instance with no precomputed values.
188 */
189 public static FactorialLog create() {
190 return new FactorialLog(0, null);
191 }
192
193 /**
194 * Creates an instance with the specified cache size.
195 *
196 * @param cacheSize Number of precomputed values of the function.
197 * @return a new instance where {@code cacheSize} values have been
198 * precomputed.
199 * @throws IllegalArgumentException if {@code n < 0}.
200 */
201 public FactorialLog withCache(final int cacheSize) {
202 return new FactorialLog(cacheSize, logFactorials);
203 }
204
205 /**
206 * Computes {@code log(n!)}.
207 *
208 * @param n Argument.
209 * @return {@code log(n!)}.
210 * @throws IndexOutOfBoundsException if {@code numValues < 0}.
211 */
212 public double value(final int n) {
213 // Use cache of precomputed values.
214 if (n < logFactorials.length) {
215 return logFactorials[n];
216 }
217
218 // Use cache of precomputed factorial values.
219 if (n < FACTORIALS.length) {
220 return Math.log(FACTORIALS[n]);
221 }
222
223 // Delegate.
224 return InternalGamma.logGamma(n + 1.0);
225 }
226 }
227 }