FastLoadedDiceRollerDiscreteSampler.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.rng.sampling.distribution;
import java.math.BigInteger;
import java.util.Arrays;
import org.apache.commons.rng.UniformRandomProvider;
/**
* Distribution sampler that uses the Fast Loaded Dice Roller (FLDR). It can be used to
* sample from {@code n} values each with an associated relative weight. If all unique items
* are assigned the same weight it is more efficient to use the {@link DiscreteUniformSampler}.
*
* <p>Given a list {@code L} of {@code n} positive numbers,
* where {@code L[i]} represents the relative weight of the {@code i}th side, FLDR returns
* integer {@code i} with relative probability {@code L[i]}.
*
* <p>FLDR produces <em>exact</em> samples from the specified probability distribution.
* <ul>
* <li>For integer weights, the probability of returning {@code i} is precisely equal to the
* rational number {@code L[i] / m}, where {@code m} is the sum of {@code L}.
* <li>For floating-points weights, each weight {@code L[i]} is converted to the
* corresponding rational number {@code p[i] / q[i]} where {@code p[i]} is a positive integer and
* {@code q[i]} is a power of 2. The rational weights are then normalized (exactly) to sum to unity.
* </ul>
*
* <p>Note that if <em>exact</em> samples are not required then an alternative sampler that
* ignores very small relative weights may have improved sampling performance.
*
* <p>This implementation is based on the algorithm in:
*
* <blockquote>
* Feras A. Saad, Cameron E. Freer, Martin C. Rinard, and Vikash K. Mansinghka.
* The Fast Loaded Dice Roller: A Near-Optimal Exact Sampler for Discrete Probability
* Distributions. In AISTATS 2020: Proceedings of the 23rd International Conference on
* Artificial Intelligence and Statistics, Proceedings of Machine Learning Research 108,
* Palermo, Sicily, Italy, 2020.
* </blockquote>
*
* <p>Sampling uses {@link UniformRandomProvider#nextInt()} as the source of random bits.
*
* @see <a href="https://arxiv.org/abs/2003.03830">Saad et al (2020)
* Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics,
* PMLR 108:1036-1046.</a>
* @since 1.5
*/
public abstract class FastLoadedDiceRollerDiscreteSampler
implements SharedStateDiscreteSampler {
/**
* The maximum size of an array.
*
* <p>This value is taken from the limit in Open JDK 8 {@code java.util.ArrayList}.
* It allows VMs to reserve some header words in an array.
*/
private static final int MAX_ARRAY_SIZE = Integer.MAX_VALUE - 8;
/** The maximum biased exponent for a finite double.
* This is offset by 1023 from {@code Math.getExponent(Double.MAX_VALUE)}. */
private static final int MAX_BIASED_EXPONENT = 2046;
/** Size of the mantissa of a double. Equal to 52 bits. */
private static final int MANTISSA_SIZE = 52;
/** Mask to extract the 52-bit mantissa from a long representation of a double. */
private static final long MANTISSA_MASK = 0x000f_ffff_ffff_ffffL;
/** BigInteger representation of {@link Long#MAX_VALUE}. */
private static final BigInteger MAX_LONG = BigInteger.valueOf(Long.MAX_VALUE);
/** The maximum offset that will avoid loss of bits for a left shift of a 53-bit value.
* The value will remain positive for any shift {@code <=} this value. */
private static final int MAX_OFFSET = 10;
/** Initial value for no leaf node label. */
private static final int NO_LABEL = Integer.MAX_VALUE;
/** Name of the sampler. */
private static final String SAMPLER_NAME = "Fast Loaded Dice Roller";
/**
* Class to handle the edge case of observations in only one category.
*/
private static class FixedValueDiscreteSampler extends FastLoadedDiceRollerDiscreteSampler {
/** The sample value. */
private final int sampleValue;
/**
* @param sampleValue Sample value.
*/
FixedValueDiscreteSampler(int sampleValue) {
this.sampleValue = sampleValue;
}
@Override
public int sample() {
return sampleValue;
}
@Override
public FastLoadedDiceRollerDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
return this;
}
@Override
public String toString() {
return SAMPLER_NAME;
}
}
/**
* Class to implement the FLDR sample algorithm.
*/
private static class FLDRSampler extends FastLoadedDiceRollerDiscreteSampler {
/** Empty boolean source. This is the location of the sign-bit after 31 right shifts on
* the boolean source. */
private static final int EMPTY_BOOL_SOURCE = 1;
/** Underlying source of randomness. */
private final UniformRandomProvider rng;
/** Number of categories. */
private final int n;
/** Number of levels in the discrete distribution generating (DDG) tree.
* Equal to {@code ceil(log2(m))} where {@code m} is the sum of observations. */
private final int k;
/** Number of leaf nodes at each level. */
private final int[] h;
/** Stores the leaf node labels in increasing order. Named {@code H} in the FLDR paper. */
private final int[] lH;
/**
* Provides a bit source for booleans.
*
* <p>A cached value from a call to {@link UniformRandomProvider#nextInt()}.
*
* <p>Only stores 31-bits when full as 1 bit has already been consumed.
* The sign bit is a flag that shifts down so the source eventually equals 1
* when all bits are consumed and will trigger a refill.
*/
private int booleanSource = EMPTY_BOOL_SOURCE;
/**
* Creates a sampler.
*
* <p>The input parameters are not validated and must be correctly computed tables.
*
* @param rng Generator of uniformly distributed random numbers.
* @param n Number of categories
* @param k Number of levels in the discrete distribution generating (DDG) tree.
* Equal to {@code ceil(log2(m))} where {@code m} is the sum of observations.
* @param h Number of leaf nodes at each level.
* @param lH Stores the leaf node labels in increasing order.
*/
FLDRSampler(UniformRandomProvider rng,
int n,
int k,
int[] h,
int[] lH) {
this.rng = rng;
this.n = n;
this.k = k;
// Deliberate direct storage of input arrays
this.h = h;
this.lH = lH;
}
/**
* Creates a copy with a new source of randomness.
*
* @param rng Generator of uniformly distributed random numbers.
* @param source Source to copy.
*/
private FLDRSampler(UniformRandomProvider rng,
FLDRSampler source) {
this.rng = rng;
this.n = source.n;
this.k = source.k;
this.h = source.h;
this.lH = source.lH;
}
/** {@inheritDoc} */
@Override
public int sample() {
// ALGORITHM 5: SAMPLE
int c = 0;
int d = 0;
for (;;) {
// b = flip()
// d = 2 * d + (1 - b)
d = (d << 1) + flip();
if (d < h[c]) {
// z = H[d][c]
final int z = lH[d * k + c];
// assert z != NO_LABEL
if (z < n) {
return z;
}
d = 0;
c = 0;
} else {
d = d - h[c];
c++;
}
}
}
/**
* Provides a source of boolean bits.
*
* <p>Note: This replicates the boolean cache functionality of
* {@code o.a.c.rng.core.source32.IntProvider}. The method has been simplified to return
* an {@code int} value rather than a {@code boolean}.
*
* @return the bit (0 or 1)
*/
private int flip() {
int bits = booleanSource;
if (bits == 1) {
// Refill
bits = rng.nextInt();
// Store a refill flag in the sign bit and the unused 31 bits, return lowest bit
booleanSource = Integer.MIN_VALUE | (bits >>> 1);
return bits & 0x1;
}
// Shift down eventually triggering refill, return current lowest bit
booleanSource = bits >>> 1;
return bits & 0x1;
}
/** {@inheritDoc} */
@Override
public String toString() {
return SAMPLER_NAME + " [" + rng.toString() + "]";
}
/** {@inheritDoc} */
@Override
public FastLoadedDiceRollerDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
return new FLDRSampler(rng, this);
}
}
/** Package-private constructor. */
FastLoadedDiceRollerDiscreteSampler() {
// Intentionally empty
}
/** {@inheritDoc} */
// Redeclare the signature to return a FastLoadedDiceRollerSampler not a SharedStateLongSampler
@Override
public abstract FastLoadedDiceRollerDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng);
/**
* Creates a sampler.
*
* <p>Note: The discrete distribution generating (DDG) tree requires {@code (n + 1) * k} entries
* where {@code n} is the number of categories, {@code k == ceil(log2(m))} and {@code m}
* is the sum of the observed frequencies. An exception is raised if this cannot be allocated
* as a single array.
*
* <p>For reference the sum is limited to {@link Long#MAX_VALUE} and the value {@code k} to 63.
* The number of categories is limited to approximately {@code ((2^31 - 1) / k) = 34,087,042}
* when the sum of frequencies is large enough to create k=63.
*
* @param rng Generator of uniformly distributed random numbers.
* @param frequencies Observed frequencies of the discrete distribution.
* @return the sampler
* @throws IllegalArgumentException if {@code frequencies} is null or empty, a
* frequency is negative, the sum of all frequencies is either zero or
* above {@link Long#MAX_VALUE}, or the size of the discrete distribution generating tree
* is too large.
*/
public static FastLoadedDiceRollerDiscreteSampler of(UniformRandomProvider rng,
long[] frequencies) {
final long m = sum(frequencies);
// Obtain indices of non-zero frequencies
final int[] indices = indicesOfNonZero(frequencies);
// Edge case for 1 non-zero weight. This also handles edge case for 1 observation
// (as log2(m) == 0 will break the computation of the DDG tree).
if (indices.length == 1) {
return new FixedValueDiscreteSampler(indexOfNonZero(frequencies));
}
return createSampler(rng, frequencies, indices, m);
}
/**
* Creates a sampler.
*
* <p>Weights are converted to rational numbers {@code p / q} where {@code q} is a power of 2.
* The numerators {@code p} are scaled to use a common denominator before summing.
*
* <p>All weights are used to create the sampler. Weights with a small magnitude relative
* to the largest weight can be excluded using the constructor method with the
* relative magnitude parameter {@code alpha} (see {@link #of(UniformRandomProvider, double[], int)}).
*
* @param rng Generator of uniformly distributed random numbers.
* @param weights Weights of the discrete distribution.
* @return the sampler
* @throws IllegalArgumentException if {@code weights} is null or empty, a
* weight is negative, infinite or {@code NaN}, the sum of all weights is zero, or the size
* of the discrete distribution generating tree is too large.
* @see #of(UniformRandomProvider, double[], int)
*/
public static FastLoadedDiceRollerDiscreteSampler of(UniformRandomProvider rng,
double[] weights) {
return of(rng, weights, 0);
}
/**
* Creates a sampler.
*
* <p>Weights are converted to rational numbers {@code p / q} where {@code q} is
* a power of 2. The numerators {@code p} are scaled to use a common
* denominator before summing.
*
* <p>Note: The discrete distribution generating (DDG) tree requires
* {@code (n + 1) * k} entries where {@code n} is the number of categories,
* {@code k == ceil(log2(m))} and {@code m} is the sum of the weight numerators
* {@code q}. An exception is raised if this cannot be allocated as a single
* array.
*
* <p>For reference the value {@code k} is equal to or greater than the ratio of
* the largest to the smallest weight expressed as a power of 2. For
* {@code Double.MAX_VALUE / Double.MIN_VALUE} this is ~2098. The value
* {@code k} increases with the sum of the weight numerators. A number of
* weights in excess of 1,000,000 with values equal to {@link Double#MAX_VALUE}
* would be required to raise an exception when the minimum weight is
* {@link Double#MIN_VALUE}.
*
* <p>Weights with a small magnitude relative to the largest weight can be
* excluded using the relative magnitude parameter {@code alpha}. This will set
* any weight to zero if the magnitude is approximately 2<sup>alpha</sup>
* <em>smaller</em> than the largest weight. This comparison is made using only
* the exponent of the input weights. The {@code alpha} parameter is ignored if
* not above zero. Note that a small {@code alpha} parameter will exclude more
* weights than a large {@code alpha} parameter.
*
* <p>The alpha parameter can be used to exclude categories that
* have a very low probability of occurrence and will improve the construction
* performance of the sampler. The effect on sampling performance depends on
* the relative weights of the excluded categories; typically a high {@code alpha}
* is used to exclude categories that would be visited with a very low probability
* and the sampling performance is unchanged.
*
* <p><b>Implementation Note</b>
*
* <p>This method creates a sampler with <em>exact</em> samples from the
* specified probability distribution. It is recommended to use this method:
* <ul>
* <li>if the weights are computed, for example from a probability mass function; or
* <li>if the weights sum to an infinite value.
* </ul>
*
* <p>If the weights are computed from empirical observations then it is
* recommended to use the factory method
* {@link #of(UniformRandomProvider, long[]) accepting frequencies}. This
* requires the total number of observations to be representable as a long
* integer.
*
* <p>Note that if all weights are scaled by a power of 2 to be integers, and
* each integer can be represented as a positive 64-bit long value, then the
* sampler created using this method will match the output from a sampler
* created with the scaled weights converted to long values for the factory
* method {@link #of(UniformRandomProvider, long[]) accepting frequencies}. This
* assumes the sum of the integer values does not overflow.
*
* <p>It should be noted that the conversion of weights to rational numbers has
* a performance overhead during construction (sampling performance is not
* affected). This may be avoided by first converting them to integer values
* that can be summed without overflow. For example by scaling values by
* {@code 2^62 / sum} and converting to long by casting or rounding.
*
* <p>This approach may increase the efficiency of construction. The resulting
* sampler may no longer produce <em>exact</em> samples from the distribution.
* In particular any weights with a converted frequency of zero cannot be
* sampled.
*
* @param rng Generator of uniformly distributed random numbers.
* @param weights Weights of the discrete distribution.
* @param alpha Alpha parameter.
* @return the sampler
* @throws IllegalArgumentException if {@code weights} is null or empty, a
* weight is negative, infinite or {@code NaN}, the sum of all weights is zero,
* or the size of the discrete distribution generating tree is too large.
* @see #of(UniformRandomProvider, long[])
*/
public static FastLoadedDiceRollerDiscreteSampler of(UniformRandomProvider rng,
double[] weights,
int alpha) {
final int n = checkWeightsNonZeroLength(weights);
// Convert floating-point double to a relative weight
// using a shifted integer representation
final long[] frequencies = new long[n];
final int[] offsets = new int[n];
convertToIntegers(weights, frequencies, offsets, alpha);
// Obtain indices of non-zero weights
final int[] indices = indicesOfNonZero(frequencies);
// Edge case for 1 non-zero weight.
if (indices.length == 1) {
return new FixedValueDiscreteSampler(indexOfNonZero(frequencies));
}
final BigInteger m = sum(frequencies, offsets, indices);
// Use long arithmetic if possible. This occurs when the weights are similar in magnitude.
if (m.compareTo(MAX_LONG) <= 0) {
// Apply the offset
for (int i = 0; i < n; i++) {
frequencies[i] <<= offsets[i];
}
return createSampler(rng, frequencies, indices, m.longValue());
}
return createSampler(rng, frequencies, offsets, indices, m);
}
/**
* Sum the frequencies.
*
* @param frequencies Frequencies.
* @return the sum
* @throws IllegalArgumentException if {@code frequencies} is null or empty, a
* frequency is negative, or the sum of all frequencies is either zero or above
* {@link Long#MAX_VALUE}
*/
private static long sum(long[] frequencies) {
// Validate
if (frequencies == null || frequencies.length == 0) {
throw new IllegalArgumentException("frequencies must contain at least 1 value");
}
// Sum the values.
// Combine all the sign bits in the observations and the intermediate sum in a flag.
long m = 0;
long signFlag = 0;
for (final long o : frequencies) {
m += o;
signFlag |= o | m;
}
// Check for a sign-bit.
if (signFlag < 0) {
// One or more observations were negative, or the sum overflowed.
for (final long o : frequencies) {
if (o < 0) {
throw new IllegalArgumentException("frequencies must contain positive values: " + o);
}
}
throw new IllegalArgumentException("Overflow when summing frequencies");
}
if (m == 0) {
throw new IllegalArgumentException("Sum of frequencies is zero");
}
return m;
}
/**
* Convert the floating-point weights to relative weights represented as
* integers {@code value * 2^exponent}. The relative weight as an integer is:
*
* <pre>
* BigInteger.valueOf(value).shiftLeft(exponent)
* </pre>
*
* <p>Note that the weights are created using a common power-of-2 scaling
* operation so the minimum exponent is zero.
*
* <p>A positive {@code alpha} parameter is used to set any weight to zero if
* the magnitude is approximately 2<sup>alpha</sup> <em>smaller</em> than the
* largest weight. This comparison is made using only the exponent of the input
* weights.
*
* @param weights Weights of the discrete distribution.
* @param values Output floating-point mantissas converted to 53-bit integers.
* @param exponents Output power of 2 exponent.
* @param alpha Alpha parameter.
* @throws IllegalArgumentException if a weight is negative, infinite or
* {@code NaN}, or the sum of all weights is zero.
*/
private static void convertToIntegers(double[] weights, long[] values, int[] exponents, int alpha) {
int maxExponent = Integer.MIN_VALUE;
for (int i = 0; i < weights.length; i++) {
final double weight = weights[i];
// Ignore zero.
// When creating the integer value later using bit shifts the result will remain zero.
if (weight == 0) {
continue;
}
final long bits = Double.doubleToRawLongBits(weight);
// For the IEEE 754 format see Double.longBitsToDouble(long).
// Extract the exponent (with the sign bit)
int exp = (int) (bits >>> MANTISSA_SIZE);
// Detect negative, infinite or NaN.
// Note: Negative values sign bit will cause the exponent to be too high.
if (exp > MAX_BIASED_EXPONENT) {
throw new IllegalArgumentException("Invalid weight: " + weight);
}
long mantissa;
if (exp == 0) {
// Sub-normal number:
mantissa = (bits & MANTISSA_MASK) << 1;
// Here we convert to a normalised number by counting the leading zeros
// to obtain the number of shifts of the most significant bit in
// the mantissa that is required to get a 1 at position 53 (i.e. as
// if it were a normal number with assumed leading bit).
final int shift = Long.numberOfLeadingZeros(mantissa << 11);
mantissa <<= shift;
exp -= shift;
} else {
// Normal number. Add the implicit leading 1-bit.
mantissa = (bits & MANTISSA_MASK) | (1L << MANTISSA_SIZE);
}
// Here the floating-point number is equal to:
// mantissa * 2^(exp-1075)
values[i] = mantissa;
exponents[i] = exp;
maxExponent = Math.max(maxExponent, exp);
}
// No exponent indicates that all weights are zero
if (maxExponent == Integer.MIN_VALUE) {
throw new IllegalArgumentException("Sum of weights is zero");
}
filterWeights(values, exponents, alpha, maxExponent);
scaleWeights(values, exponents);
}
/**
* Filters small weights using the {@code alpha} parameter.
* A positive {@code alpha} parameter is used to set any weight to zero if
* the magnitude is approximately 2<sup>alpha</sup> <em>smaller</em> than the
* largest weight. This comparison is made using only the exponent of the input
* weights.
*
* @param values 53-bit values.
* @param exponents Power of 2 exponent.
* @param alpha Alpha parameter.
* @param maxExponent Maximum exponent.
*/
private static void filterWeights(long[] values, int[] exponents, int alpha, int maxExponent) {
if (alpha > 0) {
// Filter weights. This must be done before the values are shifted so
// the exponent represents the approximate magnitude of the value.
for (int i = 0; i < exponents.length; i++) {
if (maxExponent - exponents[i] > alpha) {
values[i] = 0;
}
}
}
}
/**
* Scale the weights represented as integers {@code value * 2^exponent} to use a
* minimum exponent of zero. The values are scaled to remove any common trailing zeros
* in their representation. This ultimately reduces the size of the discrete distribution
* generating (DGG) tree.
*
* @param values 53-bit values.
* @param exponents Power of 2 exponent.
*/
private static void scaleWeights(long[] values, int[] exponents) {
// Find the minimum exponent and common trailing zeros.
int minExponent = Integer.MAX_VALUE;
for (int i = 0; i < exponents.length; i++) {
if (values[i] != 0) {
minExponent = Math.min(minExponent, exponents[i]);
}
}
// Trailing zeros occur when the original weights have a representation with
// less than 52 binary digits, e.g. {1.5, 0.5, 0.25}.
int trailingZeros = Long.SIZE;
for (int i = 0; i < values.length && trailingZeros != 0; i++) {
trailingZeros = Math.min(trailingZeros, Long.numberOfTrailingZeros(values[i]));
}
// Scale by a power of 2 so the minimum exponent is zero.
for (int i = 0; i < exponents.length; i++) {
exponents[i] -= minExponent;
}
// Remove common trailing zeros.
if (trailingZeros != 0) {
for (int i = 0; i < values.length; i++) {
values[i] >>>= trailingZeros;
}
}
}
/**
* Sum the integers at the specified indices.
* Integers are represented as {@code value * 2^exponent}.
*
* @param values 53-bit values.
* @param exponents Power of 2 exponent.
* @param indices Indices to sum.
* @return the sum
*/
private static BigInteger sum(long[] values, int[] exponents, int[] indices) {
BigInteger m = BigInteger.ZERO;
for (final int i : indices) {
m = m.add(toBigInteger(values[i], exponents[i]));
}
return m;
}
/**
* Convert the value and left shift offset to a BigInteger.
* It is assumed the value is at most 53-bits. This allows optimising the left
* shift if it is below 11 bits.
*
* @param value 53-bit value.
* @param offset Left shift offset (must be positive).
* @return the BigInteger
*/
private static BigInteger toBigInteger(long value, int offset) {
// Ignore zeros. The sum method uses indices of non-zero values.
if (offset <= MAX_OFFSET) {
// Assume (value << offset) <= Long.MAX_VALUE
return BigInteger.valueOf(value << offset);
}
return BigInteger.valueOf(value).shiftLeft(offset);
}
/**
* Creates the sampler.
*
* <p>It is assumed the frequencies are all positive and the sum does not
* overflow.
*
* @param rng Generator of uniformly distributed random numbers.
* @param frequencies Observed frequencies of the discrete distribution.
* @param indices Indices of non-zero frequencies.
* @param m Sum of the frequencies.
* @return the sampler
*/
private static FastLoadedDiceRollerDiscreteSampler createSampler(UniformRandomProvider rng,
long[] frequencies,
int[] indices,
long m) {
// ALGORITHM 5: PREPROCESS
// a == frequencies
// m = sum(a)
// h = leaf node count
// H = leaf node label (lH)
final int n = frequencies.length;
// k = ceil(log2(m))
final int k = 64 - Long.numberOfLeadingZeros(m - 1);
// r = a(n+1) = 2^k - m
final long r = (1L << k) - m;
// Note:
// A sparse matrix can often be used for H, as most of its entries are empty.
// This implementation uses a 1D array for efficiency at the cost of memory.
// This is limited to approximately ((2^31 - 1) / k), e.g. 34087042 when the sum of
// observations is large enough to create k=63.
// This could be handled using a 2D array. In practice a number of categories this
// large is not expected and is currently not supported.
final int[] h = new int[k];
final int[] lH = new int[checkArraySize((n + 1L) * k)];
Arrays.fill(lH, NO_LABEL);
int d;
for (int j = 0; j < k; j++) {
final int shift = (k - 1) - j;
final long bitMask = 1L << shift;
d = 0;
for (final int i : indices) {
// bool w ← (a[i] >> (k − 1) − j)) & 1
// h[j] = h[j] + w
// if w then:
if ((frequencies[i] & bitMask) != 0) {
h[j]++;
// H[d][j] = i
lH[d * k + j] = i;
d++;
}
}
// process a(n+1) without extending the input frequencies array by 1
if ((r & bitMask) != 0) {
h[j]++;
lH[d * k + j] = n;
}
}
return new FLDRSampler(rng, n, k, h, lH);
}
/**
* Creates the sampler. Frequencies are represented as a 53-bit value with a
* left-shift offset.
* <pre>
* BigInteger.valueOf(value).shiftLeft(offset)
* </pre>
*
* <p>It is assumed the frequencies are all positive.
*
* @param rng Generator of uniformly distributed random numbers.
* @param frequencies Observed frequencies of the discrete distribution.
* @param offsets Left shift offsets (must be positive).
* @param indices Indices of non-zero frequencies.
* @param m Sum of the frequencies.
* @return the sampler
*/
private static FastLoadedDiceRollerDiscreteSampler createSampler(UniformRandomProvider rng,
long[] frequencies,
int[] offsets,
int[] indices,
BigInteger m) {
// Repeat the logic from createSampler(...) using extended arithmetic to test the bits
// ALGORITHM 5: PREPROCESS
// a == frequencies
// m = sum(a)
// h = leaf node count
// H = leaf node label (lH)
final int n = frequencies.length;
// k = ceil(log2(m))
final int k = m.subtract(BigInteger.ONE).bitLength();
// r = a(n+1) = 2^k - m
final BigInteger r = BigInteger.ONE.shiftLeft(k).subtract(m);
final int[] h = new int[k];
final int[] lH = new int[checkArraySize((n + 1L) * k)];
Arrays.fill(lH, NO_LABEL);
int d;
for (int j = 0; j < k; j++) {
final int shift = (k - 1) - j;
d = 0;
for (final int i : indices) {
// bool w ← (a[i] >> (k − 1) − j)) & 1
// h[j] = h[j] + w
// if w then:
if (testBit(frequencies[i], offsets[i], shift)) {
h[j]++;
// H[d][j] = i
lH[d * k + j] = i;
d++;
}
}
// process a(n+1) without extending the input frequencies array by 1
if (r.testBit(shift)) {
h[j]++;
lH[d * k + j] = n;
}
}
return new FLDRSampler(rng, n, k, h, lH);
}
/**
* Test the logical bit of the shifted integer representation.
* The value is assumed to have at most 53-bits of information. The offset
* is assumed to be positive. This is functionally equivalent to:
* <pre>
* BigInteger.valueOf(value).shiftLeft(offset).testBit(n)
* </pre>
*
* @param value 53-bit value.
* @param offset Left shift offset.
* @param n Index of bit to test.
* @return true if the bit is 1
*/
private static boolean testBit(long value, int offset, int n) {
if (n < offset) {
// All logical trailing bits are zero
return false;
}
// Test if outside the 53-bit value (note that the implicit 1 bit
// has been added to the 52-bit mantissas for 'normal' floating-point numbers).
final int bit = n - offset;
return bit <= MANTISSA_SIZE && (value & (1L << bit)) != 0;
}
/**
* Check the weights have a non-zero length.
*
* @param weights Weights.
* @return the length
*/
private static int checkWeightsNonZeroLength(double[] weights) {
if (weights == null || weights.length == 0) {
throw new IllegalArgumentException("weights must contain at least 1 value");
}
return weights.length;
}
/**
* Create the indices of non-zero values.
*
* @param values Values.
* @return the indices
*/
private static int[] indicesOfNonZero(long[] values) {
int n = 0;
final int[] indices = new int[values.length];
for (int i = 0; i < values.length; i++) {
if (values[i] != 0) {
indices[n++] = i;
}
}
return Arrays.copyOf(indices, n);
}
/**
* Find the index of the first non-zero frequency.
*
* @param frequencies Frequencies.
* @return the index
* @throws IllegalStateException if all frequencies are zero.
*/
static int indexOfNonZero(long[] frequencies) {
for (int i = 0; i < frequencies.length; i++) {
if (frequencies[i] != 0) {
return i;
}
}
throw new IllegalStateException("All frequencies are zero");
}
/**
* Check the size is valid for a 1D array.
*
* @param size Size
* @return the size as an {@code int}
* @throws IllegalArgumentException if the size is too large for a 1D array.
*/
static int checkArraySize(long size) {
if (size > MAX_ARRAY_SIZE) {
throw new IllegalArgumentException("Unable to allocate array of size: " + size);
}
return (int) size;
}
}