diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/random')
32 files changed, 6739 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/random/AbstractRandomGenerator.java b/src/main/java/org/apache/commons/math3/random/AbstractRandomGenerator.java new file mode 100644 index 0000000..ce8ad85 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/AbstractRandomGenerator.java @@ -0,0 +1,250 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.util.FastMath; + +/** + * Abstract class implementing the {@link RandomGenerator} interface. Default implementations for + * all methods other than {@link #nextDouble()} and {@link #setSeed(long)} are provided. + * + * <p>All data generation methods are based on {@code code nextDouble()}. Concrete implementations + * <strong>must</strong> override this method and <strong>should</strong> provide better / more + * performant implementations of the other methods if the underlying PRNG supplies them. + * + * @since 1.1 + */ +public abstract class AbstractRandomGenerator implements RandomGenerator { + + /** + * Cached random normal value. The default implementation for {@link #nextGaussian} generates + * pairs of values and this field caches the second value so that the full algorithm is not + * executed for every activation. The value {@code Double.NaN} signals that there is no cached + * value. Use {@link #clear} to clear the cached value. + */ + private double cachedNormalDeviate = Double.NaN; + + /** Construct a RandomGenerator. */ + public AbstractRandomGenerator() { + super(); + } + + /** + * Clears the cache used by the default implementation of {@link #nextGaussian}. Implementations + * that do not override the default implementation of {@code nextGaussian} should call this + * method in the implementation of {@link #setSeed(long)} + */ + public void clear() { + cachedNormalDeviate = Double.NaN; + } + + /** {@inheritDoc} */ + public void setSeed(int seed) { + setSeed((long) seed); + } + + /** {@inheritDoc} */ + public void setSeed(int[] seed) { + // the following number is the largest prime that fits in 32 bits (it is 2^32 - 5) + final long prime = 4294967291l; + + long combined = 0l; + for (int s : seed) { + combined = combined * prime + s; + } + setSeed(combined); + } + + /** + * Sets the seed of the underlying random number generator using a {@code long} seed. Sequences + * of values generated starting with the same seeds should be identical. + * + * <p>Implementations that do not override the default implementation of {@code nextGaussian} + * should include a call to {@link #clear} in the implementation of this method. + * + * @param seed the seed value + */ + public abstract void setSeed(long seed); + + /** + * Generates random bytes and places them into a user-supplied byte array. The number of random + * bytes produced is equal to the length of the byte array. + * + * <p>The default implementation fills the array with bytes extracted from random integers + * generated using {@link #nextInt}. + * + * @param bytes the non-null byte array in which to put the random bytes + */ + public void nextBytes(byte[] bytes) { + int bytesOut = 0; + while (bytesOut < bytes.length) { + int randInt = nextInt(); + for (int i = 0; i < 3; i++) { + if (i > 0) { + randInt >>= 8; + } + bytes[bytesOut++] = (byte) randInt; + if (bytesOut == bytes.length) { + return; + } + } + } + } + + /** + * Returns the next pseudorandom, uniformly distributed {@code int} value from this random + * number generator's sequence. All 2<font size="-1"><sup>32</sup></font> possible {@code int} + * values should be produced with (approximately) equal probability. + * + * <p>The default implementation provided here returns + * + * <pre> + * <code>(int) (nextDouble() * Integer.MAX_VALUE)</code> + * </pre> + * + * @return the next pseudorandom, uniformly distributed {@code int} value from this random + * number generator's sequence + */ + public int nextInt() { + return (int) ((2d * nextDouble() - 1d) * Integer.MAX_VALUE); + } + + /** + * Returns a pseudorandom, uniformly distributed {@code int} value between 0 (inclusive) and the + * specified value (exclusive), drawn from this random number generator's sequence. + * + * <p>The default implementation returns + * + * <pre> + * <code>(int) (nextDouble() * n</code> + * </pre> + * + * @param n the bound on the random number to be returned. Must be positive. + * @return a pseudorandom, uniformly distributed {@code int} value between 0 (inclusive) and n + * (exclusive). + * @throws NotStrictlyPositiveException if {@code n <= 0}. + */ + public int nextInt(int n) { + if (n <= 0) { + throw new NotStrictlyPositiveException(n); + } + int result = (int) (nextDouble() * n); + return result < n ? result : n - 1; + } + + /** + * Returns the next pseudorandom, uniformly distributed {@code long} value from this random + * number generator's sequence. All 2<font size="-1"><sup>64</sup></font> possible {@code long} + * values should be produced with (approximately) equal probability. + * + * <p>The default implementation returns + * + * <pre> + * <code>(long) (nextDouble() * Long.MAX_VALUE)</code> + * </pre> + * + * @return the next pseudorandom, uniformly distributed {@code long} value from this random + * number generator's sequence + */ + public long nextLong() { + return (long) ((2d * nextDouble() - 1d) * Long.MAX_VALUE); + } + + /** + * Returns the next pseudorandom, uniformly distributed {@code boolean} value from this random + * number generator's sequence. + * + * <p>The default implementation returns + * + * <pre> + * <code>nextDouble() <= 0.5</code> + * </pre> + * + * @return the next pseudorandom, uniformly distributed {@code boolean} value from this random + * number generator's sequence + */ + public boolean nextBoolean() { + return nextDouble() <= 0.5; + } + + /** + * Returns the next pseudorandom, uniformly distributed {@code float} value between {@code 0.0} + * and {@code 1.0} from this random number generator's sequence. + * + * <p>The default implementation returns + * + * <pre> + * <code>(float) nextDouble() </code> + * </pre> + * + * @return the next pseudorandom, uniformly distributed {@code float} value between {@code 0.0} + * and {@code 1.0} from this random number generator's sequence + */ + public float nextFloat() { + return (float) nextDouble(); + } + + /** + * Returns the next pseudorandom, uniformly distributed {@code double} value between {@code 0.0} + * and {@code 1.0} from this random number generator's sequence. + * + * <p>This method provides the underlying source of random data used by the other methods. + * + * @return the next pseudorandom, uniformly distributed {@code double} value between {@code 0.0} + * and {@code 1.0} from this random number generator's sequence + */ + public abstract double nextDouble(); + + /** + * Returns the next pseudorandom, Gaussian ("normally") distributed {@code double} value with + * mean {@code 0.0} and standard deviation {@code 1.0} from this random number generator's + * sequence. + * + * <p>The default implementation uses the <em>Polar Method</em> due to G.E.P. Box, M.E. Muller + * and G. Marsaglia, as described in D. Knuth, <u>The Art of Computer Programming</u>, 3.4.1C. + * + * <p>The algorithm generates a pair of independent random values. One of these is cached for + * reuse, so the full algorithm is not executed on each activation. Implementations that do not + * override this method should make sure to call {@link #clear} to clear the cached value in the + * implementation of {@link #setSeed(long)}. + * + * @return the next pseudorandom, Gaussian ("normally") distributed {@code double} value with + * mean {@code 0.0} and standard deviation {@code 1.0} from this random number generator's + * sequence + */ + public double nextGaussian() { + if (!Double.isNaN(cachedNormalDeviate)) { + double dev = cachedNormalDeviate; + cachedNormalDeviate = Double.NaN; + return dev; + } + double v1 = 0; + double v2 = 0; + double s = 1; + while (s >= 1) { + v1 = 2 * nextDouble() - 1; + v2 = 2 * nextDouble() - 1; + s = v1 * v1 + v2 * v2; + } + if (s != 0) { + s = FastMath.sqrt(-2 * FastMath.log(s) / s); + } + cachedNormalDeviate = v2 * s; + return v1 * s; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/AbstractWell.java b/src/main/java/org/apache/commons/math3/random/AbstractWell.java new file mode 100644 index 0000000..87f2ee1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/AbstractWell.java @@ -0,0 +1,216 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.util.FastMath; + +import java.io.Serializable; + +/** + * This abstract class implements the WELL class of pseudo-random number generator from + * François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + * + * <p>This generator is described in a paper by François Panneton, Pierre L'Ecuyer and Makoto + * Matsumoto <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM Transactions on Mathematical + * Software, 32, 1 (2006). The errata for the paper are in <a + * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>. + * + * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number + * generator</a> + * @since 2.2 + */ +public abstract class AbstractWell extends BitsStreamGenerator implements Serializable { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -817701723016583596L; + + /** Current index in the bytes pool. */ + protected int index; + + /** Bytes pool. */ + protected final int[] v; + + /** + * Index indirection table giving for each index its predecessor taking table size into account. + */ + protected final int[] iRm1; + + /** + * Index indirection table giving for each index its second predecessor taking table size into + * account. + */ + protected final int[] iRm2; + + /** + * Index indirection table giving for each index the value index + m1 taking table size into + * account. + */ + protected final int[] i1; + + /** + * Index indirection table giving for each index the value index + m2 taking table size into + * account. + */ + protected final int[] i2; + + /** + * Index indirection table giving for each index the value index + m3 taking table size into + * account. + */ + protected final int[] i3; + + /** + * Creates a new random number generator. + * + * <p>The instance is initialized using the current time plus the system identity hash code of + * this instance as the seed. + * + * @param k number of bits in the pool (not necessarily a multiple of 32) + * @param m1 first parameter of the algorithm + * @param m2 second parameter of the algorithm + * @param m3 third parameter of the algorithm + */ + protected AbstractWell(final int k, final int m1, final int m2, final int m3) { + this(k, m1, m2, m3, null); + } + + /** + * Creates a new random number generator using a single int seed. + * + * @param k number of bits in the pool (not necessarily a multiple of 32) + * @param m1 first parameter of the algorithm + * @param m2 second parameter of the algorithm + * @param m3 third parameter of the algorithm + * @param seed the initial seed (32 bits integer) + */ + protected AbstractWell(final int k, final int m1, final int m2, final int m3, final int seed) { + this(k, m1, m2, m3, new int[] {seed}); + } + + /** + * Creates a new random number generator using an int array seed. + * + * @param k number of bits in the pool (not necessarily a multiple of 32) + * @param m1 first parameter of the algorithm + * @param m2 second parameter of the algorithm + * @param m3 third parameter of the algorithm + * @param seed the initial seed (32 bits integers array), if null the seed of the generator will + * be related to the current time + */ + protected AbstractWell( + final int k, final int m1, final int m2, final int m3, final int[] seed) { + + // the bits pool contains k bits, k = r w - p where r is the number + // of w bits blocks, w is the block size (always 32 in the original paper) + // and p is the number of unused bits in the last block + final int w = 32; + final int r = (k + w - 1) / w; + this.v = new int[r]; + this.index = 0; + + // precompute indirection index tables. These tables are used for optimizing access + // they allow saving computations like "(j + r - 2) % r" with costly modulo operations + iRm1 = new int[r]; + iRm2 = new int[r]; + i1 = new int[r]; + i2 = new int[r]; + i3 = new int[r]; + for (int j = 0; j < r; ++j) { + iRm1[j] = (j + r - 1) % r; + iRm2[j] = (j + r - 2) % r; + i1[j] = (j + m1) % r; + i2[j] = (j + m2) % r; + i3[j] = (j + m3) % r; + } + + // initialize the pool content + setSeed(seed); + } + + /** + * Creates a new random number generator using a single long seed. + * + * @param k number of bits in the pool (not necessarily a multiple of 32) + * @param m1 first parameter of the algorithm + * @param m2 second parameter of the algorithm + * @param m3 third parameter of the algorithm + * @param seed the initial seed (64 bits integer) + */ + protected AbstractWell(final int k, final int m1, final int m2, final int m3, final long seed) { + this(k, m1, m2, m3, new int[] {(int) (seed >>> 32), (int) (seed & 0xffffffffl)}); + } + + /** + * Reinitialize the generator as if just built with the given int seed. + * + * <p>The state of the generator is exactly the same as a new generator built with the same + * seed. + * + * @param seed the initial seed (32 bits integer) + */ + @Override + public void setSeed(final int seed) { + setSeed(new int[] {seed}); + } + + /** + * Reinitialize the generator as if just built with the given int array seed. + * + * <p>The state of the generator is exactly the same as a new generator built with the same + * seed. + * + * @param seed the initial seed (32 bits integers array). If null the seed of the generator will + * be the system time plus the system identity hash code of the instance. + */ + @Override + public void setSeed(final int[] seed) { + if (seed == null) { + setSeed(System.currentTimeMillis() + System.identityHashCode(this)); + return; + } + + System.arraycopy(seed, 0, v, 0, FastMath.min(seed.length, v.length)); + + if (seed.length < v.length) { + for (int i = seed.length; i < v.length; ++i) { + final long l = v[i - seed.length]; + v[i] = (int) ((1812433253l * (l ^ (l >> 30)) + i) & 0xffffffffL); + } + } + + index = 0; + clear(); // Clear normal deviate cache + } + + /** + * Reinitialize the generator as if just built with the given long seed. + * + * <p>The state of the generator is exactly the same as a new generator built with the same + * seed. + * + * @param seed the initial seed (64 bits integer) + */ + @Override + public void setSeed(final long seed) { + setSeed(new int[] {(int) (seed >>> 32), (int) (seed & 0xffffffffl)}); + } + + /** {@inheritDoc} */ + @Override + protected abstract int next(final int bits); +} diff --git a/src/main/java/org/apache/commons/math3/random/BitsStreamGenerator.java b/src/main/java/org/apache/commons/math3/random/BitsStreamGenerator.java new file mode 100644 index 0000000..07ab156 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/BitsStreamGenerator.java @@ -0,0 +1,249 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.util.FastMath; + +import java.io.Serializable; + +/** + * Base class for random number generators that generates bits streams. + * + * @since 2.0 + */ +public abstract class BitsStreamGenerator implements RandomGenerator, Serializable { + /** Serializable version identifier */ + private static final long serialVersionUID = 20130104L; + + /** Next gaussian. */ + private double nextGaussian; + + /** Creates a new random number generator. */ + public BitsStreamGenerator() { + nextGaussian = Double.NaN; + } + + /** {@inheritDoc} */ + public abstract void setSeed(int seed); + + /** {@inheritDoc} */ + public abstract void setSeed(int[] seed); + + /** {@inheritDoc} */ + public abstract void setSeed(long seed); + + /** + * Generate next pseudorandom number. + * + * <p>This method is the core generation algorithm. It is used by all the public generation + * methods for the various primitive types {@link #nextBoolean()}, {@link #nextBytes(byte[])}, + * {@link #nextDouble()}, {@link #nextFloat()}, {@link #nextGaussian()}, {@link #nextInt()}, + * {@link #next(int)} and {@link #nextLong()}. + * + * @param bits number of random bits to produce + * @return random bits generated + */ + protected abstract int next(int bits); + + /** {@inheritDoc} */ + public boolean nextBoolean() { + return next(1) != 0; + } + + /** {@inheritDoc} */ + public double nextDouble() { + final long high = ((long) next(26)) << 26; + final int low = next(26); + return (high | low) * 0x1.0p-52d; + } + + /** {@inheritDoc} */ + public float nextFloat() { + return next(23) * 0x1.0p-23f; + } + + /** {@inheritDoc} */ + public double nextGaussian() { + + final double random; + if (Double.isNaN(nextGaussian)) { + // generate a new pair of gaussian numbers + final double x = nextDouble(); + final double y = nextDouble(); + final double alpha = 2 * FastMath.PI * x; + final double r = FastMath.sqrt(-2 * FastMath.log(y)); + random = r * FastMath.cos(alpha); + nextGaussian = r * FastMath.sin(alpha); + } else { + // use the second element of the pair already generated + random = nextGaussian; + nextGaussian = Double.NaN; + } + + return random; + } + + /** {@inheritDoc} */ + public int nextInt() { + return next(32); + } + + /** + * {@inheritDoc} + * + * <p>This default implementation is copied from Apache Harmony java.util.Random (r929253). + * + * <p>Implementation notes: + * + * <ul> + * <li>If n is a power of 2, this method returns {@code (int) ((n * (long) next(31)) >> 31)}. + * <li>If n is not a power of 2, what is returned is {@code next(31) % n} with {@code + * next(31)} values rejected (i.e. regenerated) until a value that is larger than the + * remainder of {@code Integer.MAX_VALUE / n} is generated. Rejection of this initial + * segment is necessary to ensure a uniform distribution. + * </ul> + */ + public int nextInt(int n) throws IllegalArgumentException { + if (n > 0) { + if ((n & -n) == n) { + return (int) ((n * (long) next(31)) >> 31); + } + int bits; + int val; + do { + bits = next(31); + val = bits % n; + } while (bits - val + (n - 1) < 0); + return val; + } + throw new NotStrictlyPositiveException(n); + } + + /** {@inheritDoc} */ + public long nextLong() { + final long high = ((long) next(32)) << 32; + final long low = ((long) next(32)) & 0xffffffffL; + return high | low; + } + + /** + * Returns a pseudorandom, uniformly distributed {@code long} value between 0 (inclusive) and + * the specified value (exclusive), drawn from this random number generator's sequence. + * + * @param n the bound on the random number to be returned. Must be positive. + * @return a pseudorandom, uniformly distributed {@code long} value between 0 (inclusive) and n + * (exclusive). + * @throws IllegalArgumentException if n is not positive. + */ + public long nextLong(long n) throws IllegalArgumentException { + if (n > 0) { + long bits; + long val; + do { + bits = ((long) next(31)) << 32; + bits |= ((long) next(32)) & 0xffffffffL; + val = bits % n; + } while (bits - val + (n - 1) < 0); + return val; + } + throw new NotStrictlyPositiveException(n); + } + + /** Clears the cache used by the default implementation of {@link #nextGaussian}. */ + public void clear() { + nextGaussian = Double.NaN; + } + + /** + * Generates random bytes and places them into a user-supplied array. + * + * <p>The array is filled with bytes extracted from random integers. This implies that the + * number of random bytes generated may be larger than the length of the byte array. + * + * @param bytes Array in which to put the generated bytes. Cannot be {@code null}. + */ + public void nextBytes(byte[] bytes) { + nextBytesFill(bytes, 0, bytes.length); + } + + /** + * Generates random bytes and places them into a user-supplied array. + * + * <p>The array is filled with bytes extracted from random integers. This implies that the + * number of random bytes generated may be larger than the length of the byte array. + * + * @param bytes Array in which to put the generated bytes. Cannot be {@code null}. + * @param start Index at which to start inserting the generated bytes. + * @param len Number of bytes to insert. + * @throws OutOfRangeException if {@code start < 0} or {@code start >= bytes.length}. + * @throws OutOfRangeException if {@code len < 0} or {@code len > bytes.length - start}. + */ + public void nextBytes(byte[] bytes, int start, int len) { + if (start < 0 || start >= bytes.length) { + throw new OutOfRangeException(start, 0, bytes.length); + } + if (len < 0 || len > bytes.length - start) { + throw new OutOfRangeException(len, 0, bytes.length - start); + } + + nextBytesFill(bytes, start, len); + } + + /** + * Generates random bytes and places them into a user-supplied array. + * + * <p>The array is filled with bytes extracted from random integers. This implies that the + * number of random bytes generated may be larger than the length of the byte array. + * + * @param bytes Array in which to put the generated bytes. Cannot be {@code null}. + * @param start Index at which to start inserting the generated bytes. + * @param len Number of bytes to insert. + */ + private void nextBytesFill(byte[] bytes, int start, int len) { + int index = start; // Index of first insertion. + + // Index of first insertion plus multiple 4 part of length (i.e. length + // with two least significant bits unset). + final int indexLoopLimit = index + (len & 0x7ffffffc); + + // Start filling in the byte array, 4 bytes at a time. + while (index < indexLoopLimit) { + final int random = next(32); + bytes[index++] = (byte) random; + bytes[index++] = (byte) (random >>> 8); + bytes[index++] = (byte) (random >>> 16); + bytes[index++] = (byte) (random >>> 24); + } + + final int indexLimit = start + len; // Index of last insertion + 1. + + // Fill in the remaining bytes. + if (index < indexLimit) { + int random = next(32); + while (true) { + bytes[index++] = (byte) random; + if (index < indexLimit) { + random >>>= 8; + } else { + break; + } + } + } + } +} diff --git a/src/main/java/org/apache/commons/math3/random/CorrelatedRandomVectorGenerator.java b/src/main/java/org/apache/commons/math3/random/CorrelatedRandomVectorGenerator.java new file mode 100644 index 0000000..6668356 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/CorrelatedRandomVectorGenerator.java @@ -0,0 +1,178 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.RectangularCholeskyDecomposition; + +/** + * A {@link RandomVectorGenerator} that generates vectors with with correlated components. + * + * <p>Random vectors with correlated components are built by combining the uncorrelated components + * of another random vector in such a way that the resulting correlations are the ones specified by + * a positive definite covariance matrix. + * + * <p>The main use for correlated random vector generation is for Monte-Carlo simulation of physical + * problems with several variables, for example to generate error vectors to be added to a nominal + * vector. A particularly interesting case is when the generated vector should be drawn from a <a + * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">Multivariate Normal + * Distribution</a>. The approach using a Cholesky decomposition is quite usual in this case. + * However, it can be extended to other cases as long as the underlying random generator provides + * {@link NormalizedRandomGenerator normalized values} like {@link GaussianRandomGenerator} or + * {@link UniformRandomGenerator}. + * + * <p>Sometimes, the covariance matrix for a given simulation is not strictly positive definite. + * This means that the correlations are not all independent from each other. In this case, however, + * the non strictly positive elements found during the Cholesky decomposition of the covariance + * matrix should not be negative either, they should be null. Another non-conventional extension + * handling this case is used here. Rather than computing <code>C = U<sup>T</sup>.U</code> where + * <code>C</code> is the covariance matrix and <code>U</code> is an upper-triangular matrix, we + * compute <code>C = B.B<sup>T</sup></code> where <code>B</code> is a rectangular matrix having more + * rows than columns. The number of columns of <code>B</code> is the rank of the covariance matrix, + * and it is the dimension of the uncorrelated random vector that is needed to compute the component + * of the correlated vector. This class handles this situation automatically. + * + * @since 1.2 + */ +public class CorrelatedRandomVectorGenerator implements RandomVectorGenerator { + /** Mean vector. */ + private final double[] mean; + + /** Underlying generator. */ + private final NormalizedRandomGenerator generator; + + /** Storage for the normalized vector. */ + private final double[] normalized; + + /** Root of the covariance matrix. */ + private final RealMatrix root; + + /** + * Builds a correlated random vector generator from its mean vector and covariance matrix. + * + * @param mean Expected mean values for all components. + * @param covariance Covariance matrix. + * @param small Diagonal elements threshold under which column are considered to be dependent on + * previous ones and are discarded + * @param generator underlying generator for uncorrelated normalized components. + * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException if the covariance + * matrix is not strictly positive definite. + * @throws DimensionMismatchException if the mean and covariance arrays dimensions do not match. + */ + public CorrelatedRandomVectorGenerator( + double[] mean, + RealMatrix covariance, + double small, + NormalizedRandomGenerator generator) { + int order = covariance.getRowDimension(); + if (mean.length != order) { + throw new DimensionMismatchException(mean.length, order); + } + this.mean = mean.clone(); + + final RectangularCholeskyDecomposition decomposition = + new RectangularCholeskyDecomposition(covariance, small); + root = decomposition.getRootMatrix(); + + this.generator = generator; + normalized = new double[decomposition.getRank()]; + } + + /** + * Builds a null mean random correlated vector generator from its covariance matrix. + * + * @param covariance Covariance matrix. + * @param small Diagonal elements threshold under which column are considered to be dependent on + * previous ones and are discarded. + * @param generator Underlying generator for uncorrelated normalized components. + * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException if the covariance + * matrix is not strictly positive definite. + */ + public CorrelatedRandomVectorGenerator( + RealMatrix covariance, double small, NormalizedRandomGenerator generator) { + int order = covariance.getRowDimension(); + mean = new double[order]; + for (int i = 0; i < order; ++i) { + mean[i] = 0; + } + + final RectangularCholeskyDecomposition decomposition = + new RectangularCholeskyDecomposition(covariance, small); + root = decomposition.getRootMatrix(); + + this.generator = generator; + normalized = new double[decomposition.getRank()]; + } + + /** + * Get the underlying normalized components generator. + * + * @return underlying uncorrelated components generator + */ + public NormalizedRandomGenerator getGenerator() { + return generator; + } + + /** + * Get the rank of the covariance matrix. The rank is the number of independent rows in the + * covariance matrix, it is also the number of columns of the root matrix. + * + * @return rank of the square matrix. + * @see #getRootMatrix() + */ + public int getRank() { + return normalized.length; + } + + /** + * Get the root of the covariance matrix. The root is the rectangular matrix <code>B</code> such + * that the covariance matrix is equal to <code>B.B<sup>T</sup></code> + * + * @return root of the square matrix + * @see #getRank() + */ + public RealMatrix getRootMatrix() { + return root; + } + + /** + * Generate a correlated random vector. + * + * @return a random vector as an array of double. The returned array is created at each call, + * the caller can do what it wants with it. + */ + public double[] nextVector() { + + // generate uncorrelated vector + for (int i = 0; i < normalized.length; ++i) { + normalized[i] = generator.nextNormalizedDouble(); + } + + // compute correlated vector + double[] correlated = new double[mean.length]; + for (int i = 0; i < correlated.length; ++i) { + correlated[i] = mean[i]; + for (int j = 0; j < root.getColumnDimension(); ++j) { + correlated[i] += root.getEntry(i, j) * normalized[j]; + } + } + + return correlated; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/EmpiricalDistribution.java b/src/main/java/org/apache/commons/math3/random/EmpiricalDistribution.java new file mode 100644 index 0000000..9ed3f4a --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/EmpiricalDistribution.java @@ -0,0 +1,866 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.distribution.AbstractRealDistribution; +import org.apache.commons.math3.distribution.ConstantRealDistribution; +import org.apache.commons.math3.distribution.NormalDistribution; +import org.apache.commons.math3.distribution.RealDistribution; +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.exception.ZeroException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.stat.descriptive.StatisticalSummary; +import org.apache.commons.math3.stat.descriptive.SummaryStatistics; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.MathUtils; + +import java.io.BufferedReader; +import java.io.File; +import java.io.FileInputStream; +import java.io.IOException; +import java.io.InputStream; +import java.io.InputStreamReader; +import java.net.URL; +import java.nio.charset.Charset; +import java.util.ArrayList; +import java.util.List; + +/** + * Represents an <a href="http://http://en.wikipedia.org/wiki/Empirical_distribution_function"> + * empirical probability distribution</a> -- a probability distribution derived from observed data + * without making any assumptions about the functional form of the population distribution that the + * data come from. + * + * <p>An <code>EmpiricalDistribution</code> maintains data structures, called <i>distribution + * digests</i>, that describe empirical distributions and support the following operations: + * + * <ul> + * <li>loading the distribution from a file of observed data values + * <li>dividing the input data into "bin ranges" and reporting bin frequency counts (data for + * histogram) + * <li>reporting univariate statistics describing the full set of data values as well as the + * observations within each bin + * <li>generating random values from the distribution + * </ul> + * + * Applications can use <code>EmpiricalDistribution</code> to build grouped frequency histograms + * representing the input data or to generate random values "like" those in the input file -- i.e., + * the values generated will follow the distribution of the values in the file. + * + * <p>The implementation uses what amounts to the <a + * href="http://nedwww.ipac.caltech.edu/level5/March02/Silverman/Silver2_6.html">Variable Kernel + * Method</a> with Gaussian smoothing: + * + * <p><strong>Digesting the input file</strong> + * + * <ol> + * <li>Pass the file once to compute min and max. + * <li>Divide the range from min-max into <code>binCount</code> "bins." + * <li>Pass the data file again, computing bin counts and univariate statistics (mean, std dev.) + * for each of the bins + * <li>Divide the interval (0,1) into subintervals associated with the bins, with the length of a + * bin's subinterval proportional to its count. + * </ol> + * + * <strong>Generating random values from the distribution</strong> + * + * <ol> + * <li>Generate a uniformly distributed value in (0,1) + * <li>Select the subinterval to which the value belongs. + * <li>Generate a random Gaussian value with mean = mean of the associated bin and std dev = std + * dev of associated bin. + * </ol> + * + * <p>EmpiricalDistribution implements the {@link RealDistribution} interface as follows. Given x + * within the range of values in the dataset, let B be the bin containing x and let K be the + * within-bin kernel for B. Let P(B-) be the sum of the probabilities of the bins below B and let + * K(B) be the mass of B under K (i.e., the integral of the kernel density over B). Then set P(X < + * x) = P(B-) + P(B) * K(x) / K(B) where K(x) is the kernel distribution evaluated at x. This + * results in a cdf that matches the grouped frequency distribution at the bin endpoints and + * interpolates within bins using within-bin kernels. <strong>USAGE NOTES:</strong> + * + * <ul> + * <li>The <code>binCount</code> is set by default to 1000. A good rule of thumb is to set the bin + * count to approximately the length of the input file divided by 10. + * <li>The input file <i>must</i> be a plain text file containing one valid numeric entry per + * line. + * </ul> + */ +public class EmpiricalDistribution extends AbstractRealDistribution { + + /** Default bin count */ + public static final int DEFAULT_BIN_COUNT = 1000; + + /** Character set for file input */ + private static final String FILE_CHARSET = "US-ASCII"; + + /** Serializable version identifier */ + private static final long serialVersionUID = 5729073523949762654L; + + /** RandomDataGenerator instance to use in repeated calls to getNext() */ + protected final RandomDataGenerator randomData; + + /** List of SummaryStatistics objects characterizing the bins */ + private final List<SummaryStatistics> binStats; + + /** Sample statistics */ + private SummaryStatistics sampleStats = null; + + /** Max loaded value */ + private double max = Double.NEGATIVE_INFINITY; + + /** Min loaded value */ + private double min = Double.POSITIVE_INFINITY; + + /** Grid size */ + private double delta = 0d; + + /** number of bins */ + private final int binCount; + + /** is the distribution loaded? */ + private boolean loaded = false; + + /** upper bounds of subintervals in (0,1) "belonging" to the bins */ + private double[] upperBounds = null; + + /** Creates a new EmpiricalDistribution with the default bin count. */ + public EmpiricalDistribution() { + this(DEFAULT_BIN_COUNT); + } + + /** + * Creates a new EmpiricalDistribution with the specified bin count. + * + * @param binCount number of bins. Must be strictly positive. + * @throws NotStrictlyPositiveException if {@code binCount <= 0}. + */ + public EmpiricalDistribution(int binCount) { + this(binCount, new RandomDataGenerator()); + } + + /** + * Creates a new EmpiricalDistribution with the specified bin count using the provided {@link + * RandomGenerator} as the source of random data. + * + * @param binCount number of bins. Must be strictly positive. + * @param generator random data generator (may be null, resulting in default JDK generator) + * @throws NotStrictlyPositiveException if {@code binCount <= 0}. + * @since 3.0 + */ + public EmpiricalDistribution(int binCount, RandomGenerator generator) { + this(binCount, new RandomDataGenerator(generator)); + } + + /** + * Creates a new EmpiricalDistribution with default bin count using the provided {@link + * RandomGenerator} as the source of random data. + * + * @param generator random data generator (may be null, resulting in default JDK generator) + * @since 3.0 + */ + public EmpiricalDistribution(RandomGenerator generator) { + this(DEFAULT_BIN_COUNT, generator); + } + + /** + * Creates a new EmpiricalDistribution with the specified bin count using the provided {@link + * RandomDataImpl} instance as the source of random data. + * + * @param binCount number of bins + * @param randomData random data generator (may be null, resulting in default JDK generator) + * @since 3.0 + * @deprecated As of 3.1. Please use {@link #EmpiricalDistribution(int,RandomGenerator)} + * instead. + */ + @Deprecated + public EmpiricalDistribution(int binCount, RandomDataImpl randomData) { + this(binCount, randomData.getDelegate()); + } + + /** + * Creates a new EmpiricalDistribution with default bin count using the provided {@link + * RandomDataImpl} as the source of random data. + * + * @param randomData random data generator (may be null, resulting in default JDK generator) + * @since 3.0 + * @deprecated As of 3.1. Please use {@link #EmpiricalDistribution(RandomGenerator)} instead. + */ + @Deprecated + public EmpiricalDistribution(RandomDataImpl randomData) { + this(DEFAULT_BIN_COUNT, randomData); + } + + /** + * Private constructor to allow lazy initialisation of the RNG contained in the {@link + * #randomData} instance variable. + * + * @param binCount number of bins. Must be strictly positive. + * @param randomData Random data generator. + * @throws NotStrictlyPositiveException if {@code binCount <= 0}. + */ + private EmpiricalDistribution(int binCount, RandomDataGenerator randomData) { + super(randomData.getRandomGenerator()); + if (binCount <= 0) { + throw new NotStrictlyPositiveException(binCount); + } + this.binCount = binCount; + this.randomData = randomData; + binStats = new ArrayList<SummaryStatistics>(); + } + + /** + * Computes the empirical distribution from the provided array of numbers. + * + * @param in the input data array + * @exception NullArgumentException if in is null + */ + public void load(double[] in) throws NullArgumentException { + DataAdapter da = new ArrayDataAdapter(in); + try { + da.computeStats(); + // new adapter for the second pass + fillBinStats(new ArrayDataAdapter(in)); + } catch (IOException ex) { + // Can't happen + throw new MathInternalError(); + } + loaded = true; + } + + /** + * Computes the empirical distribution using data read from a URL. + * + * <p>The input file <i>must</i> be an ASCII text file containing one valid numeric entry per + * line. + * + * @param url url of the input file + * @throws IOException if an IO error occurs + * @throws NullArgumentException if url is null + * @throws ZeroException if URL contains no data + */ + public void load(URL url) throws IOException, NullArgumentException, ZeroException { + MathUtils.checkNotNull(url); + Charset charset = Charset.forName(FILE_CHARSET); + BufferedReader in = new BufferedReader(new InputStreamReader(url.openStream(), charset)); + try { + DataAdapter da = new StreamDataAdapter(in); + da.computeStats(); + if (sampleStats.getN() == 0) { + throw new ZeroException(LocalizedFormats.URL_CONTAINS_NO_DATA, url); + } + // new adapter for the second pass + in = new BufferedReader(new InputStreamReader(url.openStream(), charset)); + fillBinStats(new StreamDataAdapter(in)); + loaded = true; + } finally { + try { + in.close(); + } catch (IOException ex) { // NOPMD + // ignore + } + } + } + + /** + * Computes the empirical distribution from the input file. + * + * <p>The input file <i>must</i> be an ASCII text file containing one valid numeric entry per + * line. + * + * @param file the input file + * @throws IOException if an IO error occurs + * @throws NullArgumentException if file is null + */ + public void load(File file) throws IOException, NullArgumentException { + MathUtils.checkNotNull(file); + Charset charset = Charset.forName(FILE_CHARSET); + InputStream is = new FileInputStream(file); + BufferedReader in = new BufferedReader(new InputStreamReader(is, charset)); + try { + DataAdapter da = new StreamDataAdapter(in); + da.computeStats(); + // new adapter for second pass + is = new FileInputStream(file); + in = new BufferedReader(new InputStreamReader(is, charset)); + fillBinStats(new StreamDataAdapter(in)); + loaded = true; + } finally { + try { + in.close(); + } catch (IOException ex) { // NOPMD + // ignore + } + } + } + + /** + * Provides methods for computing <code>sampleStats</code> and <code>beanStats</code> + * abstracting the source of data. + */ + private abstract class DataAdapter { + + /** + * Compute bin stats. + * + * @throws IOException if an error occurs computing bin stats + */ + public abstract void computeBinStats() throws IOException; + + /** + * Compute sample statistics. + * + * @throws IOException if an error occurs computing sample stats + */ + public abstract void computeStats() throws IOException; + } + + /** <code>DataAdapter</code> for data provided through some input stream */ + private class StreamDataAdapter extends DataAdapter { + + /** Input stream providing access to the data */ + private BufferedReader inputStream; + + /** + * Create a StreamDataAdapter from a BufferedReader + * + * @param in BufferedReader input stream + */ + StreamDataAdapter(BufferedReader in) { + super(); + inputStream = in; + } + + /** {@inheritDoc} */ + @Override + public void computeBinStats() throws IOException { + String str = null; + double val = 0.0d; + while ((str = inputStream.readLine()) != null) { + val = Double.parseDouble(str); + SummaryStatistics stats = binStats.get(findBin(val)); + stats.addValue(val); + } + + inputStream.close(); + inputStream = null; + } + + /** {@inheritDoc} */ + @Override + public void computeStats() throws IOException { + String str = null; + double val = 0.0; + sampleStats = new SummaryStatistics(); + while ((str = inputStream.readLine()) != null) { + val = Double.parseDouble(str); + sampleStats.addValue(val); + } + inputStream.close(); + inputStream = null; + } + } + + /** <code>DataAdapter</code> for data provided as array of doubles. */ + private class ArrayDataAdapter extends DataAdapter { + + /** Array of input data values */ + private double[] inputArray; + + /** + * Construct an ArrayDataAdapter from a double[] array + * + * @param in double[] array holding the data + * @throws NullArgumentException if in is null + */ + ArrayDataAdapter(double[] in) throws NullArgumentException { + super(); + MathUtils.checkNotNull(in); + inputArray = in; + } + + /** {@inheritDoc} */ + @Override + public void computeStats() throws IOException { + sampleStats = new SummaryStatistics(); + for (int i = 0; i < inputArray.length; i++) { + sampleStats.addValue(inputArray[i]); + } + } + + /** {@inheritDoc} */ + @Override + public void computeBinStats() throws IOException { + for (int i = 0; i < inputArray.length; i++) { + SummaryStatistics stats = binStats.get(findBin(inputArray[i])); + stats.addValue(inputArray[i]); + } + } + } + + /** + * Fills binStats array (second pass through data file). + * + * @param da object providing access to the data + * @throws IOException if an IO error occurs + */ + private void fillBinStats(final DataAdapter da) throws IOException { + // Set up grid + min = sampleStats.getMin(); + max = sampleStats.getMax(); + delta = (max - min) / ((double) binCount); + + // Initialize binStats ArrayList + if (!binStats.isEmpty()) { + binStats.clear(); + } + for (int i = 0; i < binCount; i++) { + SummaryStatistics stats = new SummaryStatistics(); + binStats.add(i, stats); + } + + // Filling data in binStats Array + da.computeBinStats(); + + // Assign upperBounds based on bin counts + upperBounds = new double[binCount]; + upperBounds[0] = ((double) binStats.get(0).getN()) / (double) sampleStats.getN(); + for (int i = 1; i < binCount - 1; i++) { + upperBounds[i] = + upperBounds[i - 1] + + ((double) binStats.get(i).getN()) / (double) sampleStats.getN(); + } + upperBounds[binCount - 1] = 1.0d; + } + + /** + * Returns the index of the bin to which the given value belongs + * + * @param value the value whose bin we are trying to find + * @return the index of the bin containing the value + */ + private int findBin(double value) { + return FastMath.min( + FastMath.max((int) FastMath.ceil((value - min) / delta) - 1, 0), binCount - 1); + } + + /** + * Generates a random value from this distribution. <strong>Preconditions:</strong> + * + * <ul> + * <li>the distribution must be loaded before invoking this method + * </ul> + * + * @return the random value. + * @throws MathIllegalStateException if the distribution has not been loaded + */ + public double getNextValue() throws MathIllegalStateException { + + if (!loaded) { + throw new MathIllegalStateException(LocalizedFormats.DISTRIBUTION_NOT_LOADED); + } + + return sample(); + } + + /** + * Returns a {@link StatisticalSummary} describing this distribution. + * <strong>Preconditions:</strong> + * + * <ul> + * <li>the distribution must be loaded before invoking this method + * </ul> + * + * @return the sample statistics + * @throws IllegalStateException if the distribution has not been loaded + */ + public StatisticalSummary getSampleStats() { + return sampleStats; + } + + /** + * Returns the number of bins. + * + * @return the number of bins. + */ + public int getBinCount() { + return binCount; + } + + /** + * Returns a List of {@link SummaryStatistics} instances containing statistics describing the + * values in each of the bins. The list is indexed on the bin number. + * + * @return List of bin statistics. + */ + public List<SummaryStatistics> getBinStats() { + return binStats; + } + + /** + * Returns a fresh copy of the array of upper bounds for the bins. Bins are: <br> + * [min,upperBounds[0]],(upperBounds[0],upperBounds[1]],..., (upperBounds[binCount-2], + * upperBounds[binCount-1] = max]. + * + * <p>Note: In versions 1.0-2.0 of commons-math, this method incorrectly returned the array of + * probability generator upper bounds now returned by {@link #getGeneratorUpperBounds()}. + * + * @return array of bin upper bounds + * @since 2.1 + */ + public double[] getUpperBounds() { + double[] binUpperBounds = new double[binCount]; + for (int i = 0; i < binCount - 1; i++) { + binUpperBounds[i] = min + delta * (i + 1); + } + binUpperBounds[binCount - 1] = max; + return binUpperBounds; + } + + /** + * Returns a fresh copy of the array of upper bounds of the subintervals of [0,1] used in + * generating data from the empirical distribution. Subintervals correspond to bins with lengths + * proportional to bin counts. <strong>Preconditions:</strong> + * + * <ul> + * <li>the distribution must be loaded before invoking this method + * </ul> + * + * <p>In versions 1.0-2.0 of commons-math, this array was (incorrectly) returned by {@link + * #getUpperBounds()}. + * + * @since 2.1 + * @return array of upper bounds of subintervals used in data generation + * @throws NullPointerException unless a {@code load} method has been called beforehand. + */ + public double[] getGeneratorUpperBounds() { + int len = upperBounds.length; + double[] out = new double[len]; + System.arraycopy(upperBounds, 0, out, 0, len); + return out; + } + + /** + * Property indicating whether or not the distribution has been loaded. + * + * @return true if the distribution has been loaded + */ + public boolean isLoaded() { + return loaded; + } + + /** + * Reseeds the random number generator used by {@link #getNextValue()}. + * + * @param seed random generator seed + * @since 3.0 + */ + public void reSeed(long seed) { + randomData.reSeed(seed); + } + + // Distribution methods --------------------------- + + /** + * {@inheritDoc} + * + * @since 3.1 + */ + @Override + public double probability(double x) { + return 0; + } + + /** + * {@inheritDoc} + * + * <p>Returns the kernel density normalized so that its integral over each bin equals the bin + * mass. + * + * <p>Algorithm description: + * + * <ol> + * <li>Find the bin B that x belongs to. + * <li>Compute K(B) = the mass of B with respect to the within-bin kernel (i.e., the integral + * of the kernel density over B). + * <li>Return k(x) * P(B) / K(B), where k is the within-bin kernel density and P(B) is the + * mass of B. + * </ol> + * + * @since 3.1 + */ + public double density(double x) { + if (x < min || x > max) { + return 0d; + } + final int binIndex = findBin(x); + final RealDistribution kernel = getKernel(binStats.get(binIndex)); + return kernel.density(x) * pB(binIndex) / kB(binIndex); + } + + /** + * {@inheritDoc} + * + * <p>Algorithm description: + * + * <ol> + * <li>Find the bin B that x belongs to. + * <li>Compute P(B) = the mass of B and P(B-) = the combined mass of the bins below B. + * <li>Compute K(B) = the probability mass of B with respect to the within-bin kernel and + * K(B-) = the kernel distribution evaluated at the lower endpoint of B + * <li>Return P(B-) + P(B) * [K(x) - K(B-)] / K(B) where K(x) is the within-bin kernel + * distribution function evaluated at x. + * </ol> + * + * If K is a constant distribution, we return P(B-) + P(B) (counting the full mass of B). + * + * @since 3.1 + */ + public double cumulativeProbability(double x) { + if (x < min) { + return 0d; + } else if (x >= max) { + return 1d; + } + final int binIndex = findBin(x); + final double pBminus = pBminus(binIndex); + final double pB = pB(binIndex); + final RealDistribution kernel = k(x); + if (kernel instanceof ConstantRealDistribution) { + if (x < kernel.getNumericalMean()) { + return pBminus; + } else { + return pBminus + pB; + } + } + final double[] binBounds = getUpperBounds(); + final double kB = kB(binIndex); + final double lower = binIndex == 0 ? min : binBounds[binIndex - 1]; + final double withinBinCum = + (kernel.cumulativeProbability(x) - kernel.cumulativeProbability(lower)) / kB; + return pBminus + pB * withinBinCum; + } + + /** + * {@inheritDoc} + * + * <p>Algorithm description: + * + * <ol> + * <li>Find the smallest i such that the sum of the masses of the bins through i is at least + * p. + * <li>Let K be the within-bin kernel distribution for bin i.</br> Let K(B) be the mass of B + * under K. <br> + * Let K(B-) be K evaluated at the lower endpoint of B (the combined mass of the bins + * below B under K).<br> + * Let P(B) be the probability of bin i.<br> + * Let P(B-) be the sum of the bin masses below bin i. <br> + * Let pCrit = p - P(B-)<br> + * <li>Return the inverse of K evaluated at <br> + * K(B-) + pCrit * K(B) / P(B) + * </ol> + * + * @since 3.1 + */ + @Override + public double inverseCumulativeProbability(final double p) throws OutOfRangeException { + if (p < 0.0 || p > 1.0) { + throw new OutOfRangeException(p, 0, 1); + } + + if (p == 0.0) { + return getSupportLowerBound(); + } + + if (p == 1.0) { + return getSupportUpperBound(); + } + + int i = 0; + while (cumBinP(i) < p) { + i++; + } + + final RealDistribution kernel = getKernel(binStats.get(i)); + final double kB = kB(i); + final double[] binBounds = getUpperBounds(); + final double lower = i == 0 ? min : binBounds[i - 1]; + final double kBminus = kernel.cumulativeProbability(lower); + final double pB = pB(i); + final double pBminus = pBminus(i); + final double pCrit = p - pBminus; + if (pCrit <= 0) { + return lower; + } + return kernel.inverseCumulativeProbability(kBminus + pCrit * kB / pB); + } + + /** + * {@inheritDoc} + * + * @since 3.1 + */ + public double getNumericalMean() { + return sampleStats.getMean(); + } + + /** + * {@inheritDoc} + * + * @since 3.1 + */ + public double getNumericalVariance() { + return sampleStats.getVariance(); + } + + /** + * {@inheritDoc} + * + * @since 3.1 + */ + public double getSupportLowerBound() { + return min; + } + + /** + * {@inheritDoc} + * + * @since 3.1 + */ + public double getSupportUpperBound() { + return max; + } + + /** + * {@inheritDoc} + * + * @since 3.1 + */ + public boolean isSupportLowerBoundInclusive() { + return true; + } + + /** + * {@inheritDoc} + * + * @since 3.1 + */ + public boolean isSupportUpperBoundInclusive() { + return true; + } + + /** + * {@inheritDoc} + * + * @since 3.1 + */ + public boolean isSupportConnected() { + return true; + } + + /** + * {@inheritDoc} + * + * @since 3.1 + */ + @Override + public void reseedRandomGenerator(long seed) { + randomData.reSeed(seed); + } + + /** + * The probability of bin i. + * + * @param i the index of the bin + * @return the probability that selection begins in bin i + */ + private double pB(int i) { + return i == 0 ? upperBounds[0] : upperBounds[i] - upperBounds[i - 1]; + } + + /** + * The combined probability of the bins up to but not including bin i. + * + * @param i the index of the bin + * @return the probability that selection begins in a bin below bin i. + */ + private double pBminus(int i) { + return i == 0 ? 0 : upperBounds[i - 1]; + } + + /** + * Mass of bin i under the within-bin kernel of the bin. + * + * @param i index of the bin + * @return the difference in the within-bin kernel cdf between the upper and lower endpoints of + * bin i + */ + @SuppressWarnings("deprecation") + private double kB(int i) { + final double[] binBounds = getUpperBounds(); + final RealDistribution kernel = getKernel(binStats.get(i)); + return i == 0 + ? kernel.cumulativeProbability(min, binBounds[0]) + : kernel.cumulativeProbability(binBounds[i - 1], binBounds[i]); + } + + /** + * The within-bin kernel of the bin that x belongs to. + * + * @param x the value to locate within a bin + * @return the within-bin kernel of the bin containing x + */ + private RealDistribution k(double x) { + final int binIndex = findBin(x); + return getKernel(binStats.get(binIndex)); + } + + /** + * The combined probability of the bins up to and including binIndex. + * + * @param binIndex maximum bin index + * @return sum of the probabilities of bins through binIndex + */ + private double cumBinP(int binIndex) { + return upperBounds[binIndex]; + } + + /** + * The within-bin smoothing kernel. Returns a Gaussian distribution parameterized by {@code + * bStats}, unless the bin contains only one observation, in which case a constant distribution + * is returned. + * + * @param bStats summary statistics for the bin + * @return within-bin kernel parameterized by bStats + */ + protected RealDistribution getKernel(SummaryStatistics bStats) { + if (bStats.getN() == 1 || bStats.getVariance() == 0) { + return new ConstantRealDistribution(bStats.getMean()); + } else { + return new NormalDistribution( + randomData.getRandomGenerator(), + bStats.getMean(), + bStats.getStandardDeviation(), + NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + } + } +} diff --git a/src/main/java/org/apache/commons/math3/random/GaussianRandomGenerator.java b/src/main/java/org/apache/commons/math3/random/GaussianRandomGenerator.java new file mode 100644 index 0000000..33eec56 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/GaussianRandomGenerator.java @@ -0,0 +1,49 @@ +/* + * 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.math3.random; + +/** + * This class is a gaussian normalized random generator for scalars. + * + * <p>This class is a simple wrapper around the {@link RandomGenerator#nextGaussian} method. + * + * @since 1.2 + */ +public class GaussianRandomGenerator implements NormalizedRandomGenerator { + + /** Underlying generator. */ + private final RandomGenerator generator; + + /** + * Create a new generator. + * + * @param generator underlying random generator to use + */ + public GaussianRandomGenerator(final RandomGenerator generator) { + this.generator = generator; + } + + /** + * Generate a random scalar with null mean and unit standard deviation. + * + * @return a random scalar with null mean and unit standard deviation + */ + public double nextNormalizedDouble() { + return generator.nextGaussian(); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/HaltonSequenceGenerator.java b/src/main/java/org/apache/commons/math3/random/HaltonSequenceGenerator.java new file mode 100644 index 0000000..2b1c623 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/HaltonSequenceGenerator.java @@ -0,0 +1,195 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.util.MathUtils; + +/** + * Implementation of a Halton sequence. + * + * <p>A Halton sequence is a low-discrepancy sequence generating points in the interval [0, 1] + * according to + * + * <pre> + * H(n) = d_0 / b + d_1 / b^2 .... d_j / b^j+1 + * + * with + * + * n = d_j * b^j-1 + ... d_1 * b + d_0 * b^0 + * </pre> + * + * For higher dimensions, subsequent prime numbers are used as base, e.g. { 2, 3, 5 } for a Halton + * sequence in R^3. + * + * <p>Halton sequences are known to suffer from linear correlation for larger prime numbers, thus + * the individual digits are usually scrambled. This implementation already comes with support for + * up to 40 dimensions with optimal weight numbers from <a + * href="http://etd.lib.fsu.edu/theses/available/etd-07062004-140409/unrestricted/dissertation1.pdf"> + * H. Chi: Scrambled quasirandom sequences and their applications</a>. + * + * <p>The generator supports two modes: + * + * <ul> + * <li>sequential generation of points: {@link #nextVector()} + * <li>random access to the i-th point in the sequence: {@link #skipTo(int)} + * </ul> + * + * @see <a href="http://en.wikipedia.org/wiki/Halton_sequence">Halton sequence (Wikipedia)</a> + * @see <a href="https://lirias.kuleuven.be/bitstream/123456789/131168/1/mcm2005_bartv.pdf">On the + * Halton sequence and its scramblings</a> + * @since 3.3 + */ +public class HaltonSequenceGenerator implements RandomVectorGenerator { + + /** The first 40 primes. */ + private static final int[] PRIMES = + new int[] { + 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, + 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, + 173 + }; + + /** The optimal weights used for scrambling of the first 40 dimension. */ + private static final int[] WEIGHTS = + new int[] { + 1, 2, 3, 3, 8, 11, 12, 14, 7, 18, 12, 13, 17, 18, 29, 14, 18, 43, 41, 44, 40, 30, + 47, 65, 71, 28, 40, 60, 79, 89, 56, 50, 52, 61, 108, 56, 66, 63, 60, 66 + }; + + /** Space dimension. */ + private final int dimension; + + /** The current index in the sequence. */ + private int count = 0; + + /** The base numbers for each component. */ + private final int[] base; + + /** The scrambling weights for each component. */ + private final int[] weight; + + /** + * Construct a new Halton sequence generator for the given space dimension. + * + * @param dimension the space dimension + * @throws OutOfRangeException if the space dimension is outside the allowed range of [1, 40] + */ + public HaltonSequenceGenerator(final int dimension) throws OutOfRangeException { + this(dimension, PRIMES, WEIGHTS); + } + + /** + * Construct a new Halton sequence generator with the given base numbers and weights for each + * dimension. The length of the bases array defines the space dimension and is required to be + * > 0. + * + * @param dimension the space dimension + * @param bases the base number for each dimension, entries should be (pairwise) prime, may not + * be null + * @param weights the weights used during scrambling, may be null in which case no scrambling + * will be performed + * @throws NullArgumentException if base is null + * @throws OutOfRangeException if the space dimension is outside the range [1, len], where len + * refers to the length of the bases array + * @throws DimensionMismatchException if weights is non-null and the length of the input arrays + * differ + */ + public HaltonSequenceGenerator(final int dimension, final int[] bases, final int[] weights) + throws NullArgumentException, OutOfRangeException, DimensionMismatchException { + + MathUtils.checkNotNull(bases); + + if (dimension < 1 || dimension > bases.length) { + throw new OutOfRangeException(dimension, 1, PRIMES.length); + } + + if (weights != null && weights.length != bases.length) { + throw new DimensionMismatchException(weights.length, bases.length); + } + + this.dimension = dimension; + this.base = bases.clone(); + this.weight = weights == null ? null : weights.clone(); + count = 0; + } + + /** {@inheritDoc} */ + public double[] nextVector() { + final double[] v = new double[dimension]; + for (int i = 0; i < dimension; i++) { + int index = count; + double f = 1.0 / base[i]; + + int j = 0; + while (index > 0) { + final int digit = scramble(i, j, base[i], index % base[i]); + v[i] += f * digit; + index /= base[i]; // floor( index / base ) + f /= base[i]; + } + } + count++; + return v; + } + + /** + * Performs scrambling of digit {@code d_j} according to the formula: + * + * <pre> + * ( weight_i * d_j ) mod base + * </pre> + * + * Implementations can override this method to do a different scrambling. + * + * @param i the dimension index + * @param j the digit index + * @param b the base for this dimension + * @param digit the j-th digit + * @return the scrambled digit + */ + protected int scramble(final int i, final int j, final int b, final int digit) { + return weight != null ? (weight[i] * digit) % b : digit; + } + + /** + * Skip to the i-th point in the Halton sequence. + * + * <p>This operation can be performed in O(1). + * + * @param index the index in the sequence to skip to + * @return the i-th point in the Halton sequence + * @throws NotPositiveException if index < 0 + */ + public double[] skipTo(final int index) throws NotPositiveException { + count = index; + return nextVector(); + } + + /** + * Returns the index i of the next point in the Halton sequence that will be returned by calling + * {@link #nextVector()}. + * + * @return the index of the next point + */ + public int getNextIndex() { + return count; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/ISAACRandom.java b/src/main/java/org/apache/commons/math3/random/ISAACRandom.java new file mode 100644 index 0000000..8d10af5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/ISAACRandom.java @@ -0,0 +1,279 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.util.FastMath; + +import java.io.Serializable; + +/** + * <a href="http://burtleburtle.net/bob/rand/isaacafa.html">ISAAC: a fast cryptographic + * pseudo-random number generator</a> <br> + * ISAAC (Indirection, Shift, Accumulate, Add, and Count) generates 32-bit random numbers. ISAAC has + * been designed to be cryptographically secure and is inspired by RC4. Cycles are guaranteed to be + * at least 2<sup>40</sup> values long, and they are 2<sup>8295</sup> values long on average. The + * results are uniformly distributed, unbiased, and unpredictable unless you know the seed. <br> + * This code is based (with minor changes and improvements) on the original implementation of the + * algorithm by Bob Jenkins. <br> + * + * @since 3.0 + */ +public class ISAACRandom extends BitsStreamGenerator implements Serializable { + /** Serializable version identifier */ + private static final long serialVersionUID = 7288197941165002400L; + + /** Log of size of rsl[] and mem[] */ + private static final int SIZE_L = 8; + + /** Size of rsl[] and mem[] */ + private static final int SIZE = 1 << SIZE_L; + + /** Half-size of rsl[] and mem[] */ + private static final int H_SIZE = SIZE >> 1; + + /** For pseudo-random lookup */ + private static final int MASK = SIZE - 1 << 2; + + /** The golden ratio */ + private static final int GLD_RATIO = 0x9e3779b9; + + /** The results given to the user */ + private final int[] rsl = new int[SIZE]; + + /** The internal state */ + private final int[] mem = new int[SIZE]; + + /** Count through the results in rsl[] */ + private int count; + + /** Accumulator */ + private int isaacA; + + /** The last result */ + private int isaacB; + + /** Counter, guarantees cycle is at least 2^40 */ + private int isaacC; + + /** Service variable. */ + private final int[] arr = new int[8]; + + /** Service variable. */ + private int isaacX; + + /** Service variable. */ + private int isaacI; + + /** Service variable. */ + private int isaacJ; + + /** + * Creates a new ISAAC random number generator. <br> + * The instance is initialized using a combination of the current time and system hash code of + * the instance as the seed. + */ + public ISAACRandom() { + setSeed(System.currentTimeMillis() + System.identityHashCode(this)); + } + + /** + * Creates a new ISAAC random number generator using a single long seed. + * + * @param seed Initial seed. + */ + public ISAACRandom(long seed) { + setSeed(seed); + } + + /** + * Creates a new ISAAC random number generator using an int array seed. + * + * @param seed Initial seed. If {@code null}, the seed will be related to the current time. + */ + public ISAACRandom(int[] seed) { + setSeed(seed); + } + + /** {@inheritDoc} */ + @Override + public void setSeed(int seed) { + setSeed(new int[] {seed}); + } + + /** {@inheritDoc} */ + @Override + public void setSeed(long seed) { + setSeed(new int[] {(int) (seed >>> 32), (int) (seed & 0xffffffffL)}); + } + + /** {@inheritDoc} */ + @Override + public void setSeed(int[] seed) { + if (seed == null) { + setSeed(System.currentTimeMillis() + System.identityHashCode(this)); + return; + } + final int seedLen = seed.length; + final int rslLen = rsl.length; + System.arraycopy(seed, 0, rsl, 0, FastMath.min(seedLen, rslLen)); + if (seedLen < rslLen) { + for (int j = seedLen; j < rslLen; j++) { + long k = rsl[j - seedLen]; + rsl[j] = (int) (0x6c078965L * (k ^ k >> 30) + j & 0xffffffffL); + } + } + initState(); + } + + /** {@inheritDoc} */ + @Override + protected int next(int bits) { + if (count < 0) { + isaac(); + count = SIZE - 1; + } + return rsl[count--] >>> 32 - bits; + } + + /** Generate 256 results */ + private void isaac() { + isaacI = 0; + isaacJ = H_SIZE; + isaacB += ++isaacC; + while (isaacI < H_SIZE) { + isaac2(); + } + isaacJ = 0; + while (isaacJ < H_SIZE) { + isaac2(); + } + } + + /** Intermediate internal loop. */ + private void isaac2() { + isaacX = mem[isaacI]; + isaacA ^= isaacA << 13; + isaacA += mem[isaacJ++]; + isaac3(); + isaacX = mem[isaacI]; + isaacA ^= isaacA >>> 6; + isaacA += mem[isaacJ++]; + isaac3(); + isaacX = mem[isaacI]; + isaacA ^= isaacA << 2; + isaacA += mem[isaacJ++]; + isaac3(); + isaacX = mem[isaacI]; + isaacA ^= isaacA >>> 16; + isaacA += mem[isaacJ++]; + isaac3(); + } + + /** Lowest level internal loop. */ + private void isaac3() { + mem[isaacI] = mem[(isaacX & MASK) >> 2] + isaacA + isaacB; + isaacB = mem[(mem[isaacI] >> SIZE_L & MASK) >> 2] + isaacX; + rsl[isaacI++] = isaacB; + } + + /** Initialize, or reinitialize, this instance of rand. */ + private void initState() { + isaacA = 0; + isaacB = 0; + isaacC = 0; + for (int j = 0; j < arr.length; j++) { + arr[j] = GLD_RATIO; + } + for (int j = 0; j < 4; j++) { + shuffle(); + } + // fill in mem[] with messy stuff + for (int j = 0; j < SIZE; j += 8) { + arr[0] += rsl[j]; + arr[1] += rsl[j + 1]; + arr[2] += rsl[j + 2]; + arr[3] += rsl[j + 3]; + arr[4] += rsl[j + 4]; + arr[5] += rsl[j + 5]; + arr[6] += rsl[j + 6]; + arr[7] += rsl[j + 7]; + shuffle(); + setState(j); + } + // second pass makes all of seed affect all of mem + for (int j = 0; j < SIZE; j += 8) { + arr[0] += mem[j]; + arr[1] += mem[j + 1]; + arr[2] += mem[j + 2]; + arr[3] += mem[j + 3]; + arr[4] += mem[j + 4]; + arr[5] += mem[j + 5]; + arr[6] += mem[j + 6]; + arr[7] += mem[j + 7]; + shuffle(); + setState(j); + } + isaac(); + count = SIZE - 1; + clear(); + } + + /** Shuffle array. */ + private void shuffle() { + arr[0] ^= arr[1] << 11; + arr[3] += arr[0]; + arr[1] += arr[2]; + arr[1] ^= arr[2] >>> 2; + arr[4] += arr[1]; + arr[2] += arr[3]; + arr[2] ^= arr[3] << 8; + arr[5] += arr[2]; + arr[3] += arr[4]; + arr[3] ^= arr[4] >>> 16; + arr[6] += arr[3]; + arr[4] += arr[5]; + arr[4] ^= arr[5] << 10; + arr[7] += arr[4]; + arr[5] += arr[6]; + arr[5] ^= arr[6] >>> 4; + arr[0] += arr[5]; + arr[6] += arr[7]; + arr[6] ^= arr[7] << 8; + arr[1] += arr[6]; + arr[7] += arr[0]; + arr[7] ^= arr[0] >>> 9; + arr[2] += arr[7]; + arr[0] += arr[1]; + } + + /** + * Set the state by copying the internal arrays. + * + * @param start First index into {@link #mem} array. + */ + private void setState(int start) { + mem[start] = arr[0]; + mem[start + 1] = arr[1]; + mem[start + 2] = arr[2]; + mem[start + 3] = arr[3]; + mem[start + 4] = arr[4]; + mem[start + 5] = arr[5]; + mem[start + 6] = arr[6]; + mem[start + 7] = arr[7]; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/JDKRandomGenerator.java b/src/main/java/org/apache/commons/math3/random/JDKRandomGenerator.java new file mode 100644 index 0000000..e9fe954 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/JDKRandomGenerator.java @@ -0,0 +1,55 @@ +/* + * 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.math3.random; + +import java.util.Random; + +/** + * Extension of <code>java.util.Random</code> to implement {@link RandomGenerator}. + * + * @since 1.1 + */ +public class JDKRandomGenerator extends Random implements RandomGenerator { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -7745277476784028798L; + + /** Create a new JDKRandomGenerator with a default seed. */ + public JDKRandomGenerator() { + super(); + } + + /** + * Create a new JDKRandomGenerator with the given seed. + * + * @param seed initial seed + * @since 3.6 + */ + public JDKRandomGenerator(int seed) { + setSeed(seed); + } + + /** {@inheritDoc} */ + public void setSeed(int seed) { + setSeed((long) seed); + } + + /** {@inheritDoc} */ + public void setSeed(int[] seed) { + setSeed(RandomGeneratorFactory.convertToLong(seed)); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/MersenneTwister.java b/src/main/java/org/apache/commons/math3/random/MersenneTwister.java new file mode 100644 index 0000000..44ce85c --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/MersenneTwister.java @@ -0,0 +1,275 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.util.FastMath; + +import java.io.Serializable; + +/** + * This class implements a powerful pseudo-random number generator developed by Makoto Matsumoto and + * Takuji Nishimura during 1996-1997. + * + * <p>This generator features an extremely long period (2<sup>19937</sup>-1) and 623-dimensional + * equidistribution up to 32 bits accuracy. The home page for this generator is located at <a + * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html"> + * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html</a>. + * + * <p>This generator is described in a paper by Makoto Matsumoto and Takuji Nishimura in 1998: <a + * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">Mersenne Twister: A + * 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator</a>, ACM Transactions on + * Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3--30 + * + * <p>This class is mainly a Java port of the 2002-01-26 version of the generator written in C by + * Makoto Matsumoto and Takuji Nishimura. Here is their original copyright: + * + * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0"> + * <tr><td>Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + * All rights reserved.</td></tr> + * + * <tr><td>Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * <ol> + * <li>Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer.</li> + * <li>Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution.</li> + * <li>The names of its contributors may not be used to endorse or promote + * products derived from this software without specific prior written + * permission.</li> + * </ol></td></tr> + * + * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND + * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, + * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS + * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, + * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY + * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE + * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH + * DAMAGE.</strong></td></tr> + * </table> + * + * @since 2.0 + */ +public class MersenneTwister extends BitsStreamGenerator implements Serializable { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 8661194735290153518L; + + /** Size of the bytes pool. */ + private static final int N = 624; + + /** Period second parameter. */ + private static final int M = 397; + + /** X * MATRIX_A for X = {0, 1}. */ + private static final int[] MAG01 = {0x0, 0x9908b0df}; + + /** Bytes pool. */ + private int[] mt; + + /** Current index in the bytes pool. */ + private int mti; + + /** + * Creates a new random number generator. + * + * <p>The instance is initialized using the current time plus the system identity hash code of + * this instance as the seed. + */ + public MersenneTwister() { + mt = new int[N]; + setSeed(System.currentTimeMillis() + System.identityHashCode(this)); + } + + /** + * Creates a new random number generator using a single int seed. + * + * @param seed the initial seed (32 bits integer) + */ + public MersenneTwister(int seed) { + mt = new int[N]; + setSeed(seed); + } + + /** + * Creates a new random number generator using an int array seed. + * + * @param seed the initial seed (32 bits integers array), if null the seed of the generator will + * be related to the current time + */ + public MersenneTwister(int[] seed) { + mt = new int[N]; + setSeed(seed); + } + + /** + * Creates a new random number generator using a single long seed. + * + * @param seed the initial seed (64 bits integer) + */ + public MersenneTwister(long seed) { + mt = new int[N]; + setSeed(seed); + } + + /** + * Reinitialize the generator as if just built with the given int seed. + * + * <p>The state of the generator is exactly the same as a new generator built with the same + * seed. + * + * @param seed the initial seed (32 bits integer) + */ + @Override + public void setSeed(int seed) { + // we use a long masked by 0xffffffffL as a poor man unsigned int + long longMT = seed; + // NB: unlike original C code, we are working with java longs, the cast below makes masking + // unnecessary + mt[0] = (int) longMT; + for (mti = 1; mti < N; ++mti) { + // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. + // initializer from the 2002-01-09 C version by Makoto Matsumoto + longMT = (1812433253l * (longMT ^ (longMT >> 30)) + mti) & 0xffffffffL; + mt[mti] = (int) longMT; + } + + clear(); // Clear normal deviate cache + } + + /** + * Reinitialize the generator as if just built with the given int array seed. + * + * <p>The state of the generator is exactly the same as a new generator built with the same + * seed. + * + * @param seed the initial seed (32 bits integers array), if null the seed of the generator will + * be the current system time plus the system identity hash code of this instance + */ + @Override + public void setSeed(int[] seed) { + + if (seed == null) { + setSeed(System.currentTimeMillis() + System.identityHashCode(this)); + return; + } + + setSeed(19650218); + int i = 1; + int j = 0; + + for (int k = FastMath.max(N, seed.length); k != 0; k--) { + long l0 = (mt[i] & 0x7fffffffl) | ((mt[i] < 0) ? 0x80000000l : 0x0l); + long l1 = (mt[i - 1] & 0x7fffffffl) | ((mt[i - 1] < 0) ? 0x80000000l : 0x0l); + long l = (l0 ^ ((l1 ^ (l1 >> 30)) * 1664525l)) + seed[j] + j; // non linear + mt[i] = (int) (l & 0xffffffffl); + i++; + j++; + if (i >= N) { + mt[0] = mt[N - 1]; + i = 1; + } + if (j >= seed.length) { + j = 0; + } + } + + for (int k = N - 1; k != 0; k--) { + long l0 = (mt[i] & 0x7fffffffl) | ((mt[i] < 0) ? 0x80000000l : 0x0l); + long l1 = (mt[i - 1] & 0x7fffffffl) | ((mt[i - 1] < 0) ? 0x80000000l : 0x0l); + long l = (l0 ^ ((l1 ^ (l1 >> 30)) * 1566083941l)) - i; // non linear + mt[i] = (int) (l & 0xffffffffL); + i++; + if (i >= N) { + mt[0] = mt[N - 1]; + i = 1; + } + } + + mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array + + clear(); // Clear normal deviate cache + } + + /** + * Reinitialize the generator as if just built with the given long seed. + * + * <p>The state of the generator is exactly the same as a new generator built with the same + * seed. + * + * @param seed the initial seed (64 bits integer) + */ + @Override + public void setSeed(long seed) { + setSeed(new int[] {(int) (seed >>> 32), (int) (seed & 0xffffffffl)}); + } + + /** + * Generate next pseudorandom number. + * + * <p>This method is the core generation algorithm. It is used by all the public generation + * methods for the various primitive types {@link #nextBoolean()}, {@link #nextBytes(byte[])}, + * {@link #nextDouble()}, {@link #nextFloat()}, {@link #nextGaussian()}, {@link #nextInt()}, + * {@link #next(int)} and {@link #nextLong()}. + * + * @param bits number of random bits to produce + * @return random bits generated + */ + @Override + protected int next(int bits) { + + int y; + + if (mti >= N) { // generate N words at one time + int mtNext = mt[0]; + for (int k = 0; k < N - M; ++k) { + int mtCurr = mtNext; + mtNext = mt[k + 1]; + y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff); + mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 0x1]; + } + for (int k = N - M; k < N - 1; ++k) { + int mtCurr = mtNext; + mtNext = mt[k + 1]; + y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff); + mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 0x1]; + } + y = (mtNext & 0x80000000) | (mt[0] & 0x7fffffff); + mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 0x1]; + + mti = 0; + } + + y = mt[mti++]; + + // tempering + y ^= y >>> 11; + y ^= (y << 7) & 0x9d2c5680; + y ^= (y << 15) & 0xefc60000; + y ^= y >>> 18; + + return y >>> (32 - bits); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/NormalizedRandomGenerator.java b/src/main/java/org/apache/commons/math3/random/NormalizedRandomGenerator.java new file mode 100644 index 0000000..2abcc72 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/NormalizedRandomGenerator.java @@ -0,0 +1,38 @@ +/* + * 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.math3.random; + +/** + * This interface represent a normalized random generator for scalars. Normalized generator provide + * null mean and unit standard deviation scalars. + * + * @since 1.2 + */ +public interface NormalizedRandomGenerator { + + /** + * Generate a random scalar with null mean and unit standard deviation. + * + * <p>This method does <strong>not</strong> specify the shape of the distribution, it is the + * implementing class that provides it. The only contract here is to generate numbers with null + * mean and unit standard deviation. + * + * @return a random scalar with null mean and unit standard deviation + */ + double nextNormalizedDouble(); +} diff --git a/src/main/java/org/apache/commons/math3/random/RandomAdaptor.java b/src/main/java/org/apache/commons/math3/random/RandomAdaptor.java new file mode 100644 index 0000000..e7030d0 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/RandomAdaptor.java @@ -0,0 +1,182 @@ +/* + * 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.math3.random; + +import java.util.Random; + +/** + * Extension of <code>java.util.Random</code> wrapping a {@link RandomGenerator}. + * + * @since 1.1 + */ +public class RandomAdaptor extends Random implements RandomGenerator { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 2306581345647615033L; + + /** Wrapped randomGenerator instance */ + private final RandomGenerator randomGenerator; + + /** Prevent instantiation without a generator argument */ + @SuppressWarnings("unused") + private RandomAdaptor() { + randomGenerator = null; + } + + /** + * Construct a RandomAdaptor wrapping the supplied RandomGenerator. + * + * @param randomGenerator the wrapped generator + */ + public RandomAdaptor(RandomGenerator randomGenerator) { + this.randomGenerator = randomGenerator; + } + + /** + * Factory method to create a <code>Random</code> using the supplied <code>RandomGenerator + * </code>. + * + * @param randomGenerator wrapped RandomGenerator instance + * @return a Random instance wrapping the RandomGenerator + */ + public static Random createAdaptor(RandomGenerator randomGenerator) { + return new RandomAdaptor(randomGenerator); + } + + /** + * Returns the next pseudorandom, uniformly distributed <code>boolean</code> value from this + * random number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed <code>boolean</code> value from this + * random number generator's sequence + */ + @Override + public boolean nextBoolean() { + return randomGenerator.nextBoolean(); + } + + /** + * Generates random bytes and places them into a user-supplied byte array. The number of random + * bytes produced is equal to the length of the byte array. + * + * @param bytes the non-null byte array in which to put the random bytes + */ + @Override + public void nextBytes(byte[] bytes) { + randomGenerator.nextBytes(bytes); + } + + /** + * Returns the next pseudorandom, uniformly distributed <code>double</code> value between <code> + * 0.0</code> and <code>1.0</code> from this random number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed <code>double</code> value between <code> + * 0.0</code> and <code>1.0</code> from this random number generator's sequence + */ + @Override + public double nextDouble() { + return randomGenerator.nextDouble(); + } + + /** + * Returns the next pseudorandom, uniformly distributed <code>float</code> value between <code> + * 0.0</code> and <code>1.0</code> from this random number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed <code>float</code> value between <code> + * 0.0</code> and <code>1.0</code> from this random number generator's sequence + */ + @Override + public float nextFloat() { + return randomGenerator.nextFloat(); + } + + /** + * Returns the next pseudorandom, Gaussian ("normally") distributed <code>double</code> value + * with mean <code>0.0</code> and standard deviation <code>1.0</code> from this random number + * generator's sequence. + * + * @return the next pseudorandom, Gaussian ("normally") distributed <code>double</code> value + * with mean <code>0.0</code> and standard deviation <code>1.0</code> from this random + * number generator's sequence + */ + @Override + public double nextGaussian() { + return randomGenerator.nextGaussian(); + } + + /** + * Returns the next pseudorandom, uniformly distributed <code>int</code> value from this random + * number generator's sequence. All 2<font size="-1"><sup>32</sup></font> possible {@code int} + * values should be produced with (approximately) equal probability. + * + * @return the next pseudorandom, uniformly distributed <code>int</code> value from this random + * number generator's sequence + */ + @Override + public int nextInt() { + return randomGenerator.nextInt(); + } + + /** + * Returns a pseudorandom, uniformly distributed {@code int} value between 0 (inclusive) and the + * specified value (exclusive), drawn from this random number generator's sequence. + * + * @param n the bound on the random number to be returned. Must be positive. + * @return a pseudorandom, uniformly distributed {@code int} value between 0 (inclusive) and n + * (exclusive). + * @throws IllegalArgumentException if n is not positive. + */ + @Override + public int nextInt(int n) { + return randomGenerator.nextInt(n); + } + + /** + * Returns the next pseudorandom, uniformly distributed <code>long</code> value from this random + * number generator's sequence. All 2<font size="-1"><sup>64</sup></font> possible {@code long} + * values should be produced with (approximately) equal probability. + * + * @return the next pseudorandom, uniformly distributed <code>long</code> value from this random + * number generator's sequence + */ + @Override + public long nextLong() { + return randomGenerator.nextLong(); + } + + /** {@inheritDoc} */ + public void setSeed(int seed) { + if (randomGenerator != null) { // required to avoid NPE in constructor + randomGenerator.setSeed(seed); + } + } + + /** {@inheritDoc} */ + public void setSeed(int[] seed) { + if (randomGenerator != null) { // required to avoid NPE in constructor + randomGenerator.setSeed(seed); + } + } + + /** {@inheritDoc} */ + @Override + public void setSeed(long seed) { + if (randomGenerator != null) { // required to avoid NPE in constructor + randomGenerator.setSeed(seed); + } + } +} diff --git a/src/main/java/org/apache/commons/math3/random/RandomData.java b/src/main/java/org/apache/commons/math3/random/RandomData.java new file mode 100644 index 0000000..98662e3 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/RandomData.java @@ -0,0 +1,245 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.NotANumberException; +import org.apache.commons.math3.exception.NotFiniteNumberException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; + +import java.util.Collection; + +/** + * Random data generation utilities. + * + * @deprecated to be removed in 4.0. Use {@link RandomDataGenerator} directly + */ +@Deprecated +public interface RandomData { + /** + * Generates a random string of hex characters of length {@code len}. + * + * <p>The generated string will be random, but not cryptographically secure. To generate + * cryptographically secure strings, use {@link #nextSecureHexString(int)}. + * + * @param len the length of the string to be generated + * @return a random string of hex characters of length {@code len} + * @throws NotStrictlyPositiveException if {@code len <= 0} + */ + String nextHexString(int len) throws NotStrictlyPositiveException; + + /** + * Generates a uniformly distributed random integer between {@code lower} and {@code upper} + * (endpoints included). + * + * <p>The generated integer will be random, but not cryptographically secure. To generate + * cryptographically secure integer sequences, use {@link #nextSecureInt(int, int)}. + * + * @param lower lower bound for generated integer + * @param upper upper bound for generated integer + * @return a random integer greater than or equal to {@code lower} and less than or equal to + * {@code upper} + * @throws NumberIsTooLargeException if {@code lower >= upper} + */ + int nextInt(int lower, int upper) throws NumberIsTooLargeException; + + /** + * Generates a uniformly distributed random long integer between {@code lower} and {@code upper} + * (endpoints included). + * + * <p>The generated long integer values will be random, but not cryptographically secure. To + * generate cryptographically secure sequences of longs, use {@link #nextSecureLong(long, + * long)}. + * + * @param lower lower bound for generated long integer + * @param upper upper bound for generated long integer + * @return a random long integer greater than or equal to {@code lower} and less than or equal + * to {@code upper} + * @throws NumberIsTooLargeException if {@code lower >= upper} + */ + long nextLong(long lower, long upper) throws NumberIsTooLargeException; + + /** + * Generates a random string of hex characters from a secure random sequence. + * + * <p>If cryptographic security is not required, use {@link #nextHexString(int)}. + * + * @param len the length of the string to be generated + * @return a random string of hex characters of length {@code len} + * @throws NotStrictlyPositiveException if {@code len <= 0} + */ + String nextSecureHexString(int len) throws NotStrictlyPositiveException; + + /** + * Generates a uniformly distributed random integer between {@code lower} and {@code upper} + * (endpoints included) from a secure random sequence. + * + * <p>Sequences of integers generated using this method will be cryptographically secure. If + * cryptographic security is not required, {@link #nextInt(int, int)} should be used instead of + * this method. + * + * <p><strong>Definition</strong>: <a + * href="http://en.wikipedia.org/wiki/Cryptographically_secure_pseudo-random_number_generator"> + * Secure Random Sequence</a> + * + * @param lower lower bound for generated integer + * @param upper upper bound for generated integer + * @return a random integer greater than or equal to {@code lower} and less than or equal to + * {@code upper}. + * @throws NumberIsTooLargeException if {@code lower >= upper}. + */ + int nextSecureInt(int lower, int upper) throws NumberIsTooLargeException; + + /** + * Generates a uniformly distributed random long integer between {@code lower} and {@code upper} + * (endpoints included) from a secure random sequence. + * + * <p>Sequences of long values generated using this method will be cryptographically secure. If + * cryptographic security is not required, {@link #nextLong(long, long)} should be used instead + * of this method. + * + * <p><strong>Definition</strong>: <a + * href="http://en.wikipedia.org/wiki/Cryptographically_secure_pseudo-random_number_generator"> + * Secure Random Sequence</a> + * + * @param lower lower bound for generated integer + * @param upper upper bound for generated integer + * @return a random long integer greater than or equal to {@code lower} and less than or equal + * to {@code upper}. + * @throws NumberIsTooLargeException if {@code lower >= upper}. + */ + long nextSecureLong(long lower, long upper) throws NumberIsTooLargeException; + + /** + * Generates a random value from the Poisson distribution with the given mean. + * + * <p><strong>Definition</strong>: <a + * href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda366j.htm">Poisson + * Distribution</a> + * + * @param mean the mean of the Poisson distribution + * @return a random value following the specified Poisson distribution + * @throws NotStrictlyPositiveException if {@code mean <= 0}. + */ + long nextPoisson(double mean) throws NotStrictlyPositiveException; + + /** + * Generates a random value from the Normal (or Gaussian) distribution with specified mean and + * standard deviation. + * + * <p><strong>Definition</strong>: <a + * href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm">Normal + * Distribution</a> + * + * @param mu the mean of the distribution + * @param sigma the standard deviation of the distribution + * @return a random value following the specified Gaussian distribution + * @throws NotStrictlyPositiveException if {@code sigma <= 0}. + */ + double nextGaussian(double mu, double sigma) throws NotStrictlyPositiveException; + + /** + * Generates a random value from the exponential distribution with specified mean. + * + * <p><strong>Definition</strong>: <a + * href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda3667.htm">Exponential + * Distribution</a> + * + * @param mean the mean of the distribution + * @return a random value following the specified exponential distribution + * @throws NotStrictlyPositiveException if {@code mean <= 0}. + */ + double nextExponential(double mean) throws NotStrictlyPositiveException; + + /** + * Generates a uniformly distributed random value from the open interval {@code (lower, upper)} + * (i.e., endpoints excluded). + * + * <p><strong>Definition</strong>: <a + * href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda3662.htm">Uniform + * Distribution</a> {@code lower} and {@code upper - lower} are the <a href = + * "http://www.itl.nist.gov/div898/handbook/eda/section3/eda364.htm"> location and scale + * parameters</a>, respectively. + * + * @param lower the exclusive lower bound of the support + * @param upper the exclusive upper bound of the support + * @return a uniformly distributed random value between lower and upper (exclusive) + * @throws NumberIsTooLargeException if {@code lower >= upper} + * @throws NotFiniteNumberException if one of the bounds is infinite + * @throws NotANumberException if one of the bounds is NaN + */ + double nextUniform(double lower, double upper) + throws NumberIsTooLargeException, NotFiniteNumberException, NotANumberException; + + /** + * Generates a uniformly distributed random value from the interval {@code (lower, upper)} or + * the interval {@code [lower, upper)}. The lower bound is thus optionally included, while the + * upper bound is always excluded. + * + * <p><strong>Definition</strong>: <a + * href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda3662.htm">Uniform + * Distribution</a> {@code lower} and {@code upper - lower} are the <a href = + * "http://www.itl.nist.gov/div898/handbook/eda/section3/eda364.htm"> location and scale + * parameters</a>, respectively. + * + * @param lower the lower bound of the support + * @param upper the exclusive upper bound of the support + * @param lowerInclusive {@code true} if the lower bound is inclusive + * @return uniformly distributed random value in the {@code (lower, upper)} interval, if {@code + * lowerInclusive} is {@code false}, or in the {@code [lower, upper)} interval, if {@code + * lowerInclusive} is {@code true} + * @throws NumberIsTooLargeException if {@code lower >= upper} + * @throws NotFiniteNumberException if one of the bounds is infinite + * @throws NotANumberException if one of the bounds is NaN + */ + double nextUniform(double lower, double upper, boolean lowerInclusive) + throws NumberIsTooLargeException, NotFiniteNumberException, NotANumberException; + + /** + * Generates an integer array of length {@code k} whose entries are selected randomly, without + * repetition, from the integers {@code 0, ..., n - 1} (inclusive). + * + * <p>Generated arrays represent permutations of {@code n} taken {@code k} at a time. + * + * @param n the domain of the permutation + * @param k the size of the permutation + * @return a random {@code k}-permutation of {@code n}, as an array of integers + * @throws NumberIsTooLargeException if {@code k > n}. + * @throws NotStrictlyPositiveException if {@code k <= 0}. + */ + int[] nextPermutation(int n, int k) + throws NumberIsTooLargeException, NotStrictlyPositiveException; + + /** + * Returns an array of {@code k} objects selected randomly from the Collection {@code c}. + * + * <p>Sampling from {@code c} is without replacement; but if {@code c} contains identical + * objects, the sample may include repeats. If all elements of {@code c} are distinct, the + * resulting object array represents a <a + * href="http://rkb.home.cern.ch/rkb/AN16pp/node250.html#SECTION0002500000000000000000">Simple + * Random Sample</a> of size {@code k} from the elements of {@code c}. + * + * @param c the collection to be sampled + * @param k the size of the sample + * @return a random sample of {@code k} elements from {@code c} + * @throws NumberIsTooLargeException if {@code k > c.size()}. + * @throws NotStrictlyPositiveException if {@code k <= 0}. + */ + Object[] nextSample(Collection<?> c, int k) + throws NumberIsTooLargeException, NotStrictlyPositiveException; +} diff --git a/src/main/java/org/apache/commons/math3/random/RandomDataGenerator.java b/src/main/java/org/apache/commons/math3/random/RandomDataGenerator.java new file mode 100644 index 0000000..5d48580 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/RandomDataGenerator.java @@ -0,0 +1,780 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.distribution.BetaDistribution; +import org.apache.commons.math3.distribution.BinomialDistribution; +import org.apache.commons.math3.distribution.CauchyDistribution; +import org.apache.commons.math3.distribution.ChiSquaredDistribution; +import org.apache.commons.math3.distribution.ExponentialDistribution; +import org.apache.commons.math3.distribution.FDistribution; +import org.apache.commons.math3.distribution.GammaDistribution; +import org.apache.commons.math3.distribution.HypergeometricDistribution; +import org.apache.commons.math3.distribution.PascalDistribution; +import org.apache.commons.math3.distribution.PoissonDistribution; +import org.apache.commons.math3.distribution.TDistribution; +import org.apache.commons.math3.distribution.UniformIntegerDistribution; +import org.apache.commons.math3.distribution.WeibullDistribution; +import org.apache.commons.math3.distribution.ZipfDistribution; +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.exception.NotANumberException; +import org.apache.commons.math3.exception.NotFiniteNumberException; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.MathArrays; + +import java.io.Serializable; +import java.security.MessageDigest; +import java.security.NoSuchAlgorithmException; +import java.security.NoSuchProviderException; +import java.security.SecureRandom; +import java.util.Collection; + +/** + * Implements the {@link RandomData} interface using a {@link RandomGenerator} instance to generate + * non-secure data and a {@link java.security.SecureRandom} instance to provide data for the <code> + * nextSecureXxx</code> methods. If no <code>RandomGenerator</code> is provided in the constructor, + * the default is to use a {@link Well19937c} generator. To plug in a different implementation, + * either implement <code>RandomGenerator</code> directly or extend {@link AbstractRandomGenerator}. + * + * <p>Supports reseeding the underlying pseudo-random number generator (PRNG). The <code> + * SecurityProvider</code> and <code>Algorithm</code> used by the <code>SecureRandom</code> instance + * can also be reset. + * + * <p>For details on the default PRNGs, see {@link java.util.Random} and {@link + * java.security.SecureRandom}. + * + * <p><strong>Usage Notes</strong>: + * + * <ul> + * <li>Instance variables are used to maintain <code>RandomGenerator</code> and <code>SecureRandom + * </code> instances used in data generation. Therefore, to generate a random sequence of + * values or strings, you should use just <strong>one</strong> <code>RandomDataImpl</code> + * instance repeatedly. + * <li>The "secure" methods are *much* slower. These should be used only when a cryptographically + * secure random sequence is required. A secure random sequence is a sequence of pseudo-random + * values which, in addition to being well-dispersed (so no subsequence of values is an any + * more likely than other subsequence of the the same length), also has the additional + * property that knowledge of values generated up to any point in the sequence does not make + * it any easier to predict subsequent values. + * <li>When a new <code>RandomDataImpl</code> is created, the underlying random number generators + * are <strong>not</strong> initialized. If you do not explicitly seed the default non-secure + * generator, it is seeded with the current time in milliseconds plus the system identity hash + * code on first use. The same holds for the secure generator. If you provide a <code> + * RandomGenerator</code> to the constructor, however, this generator is not reseeded by the + * constructor nor is it reseeded on first use. + * <li>The <code>reSeed</code> and <code>reSeedSecure</code> methods delegate to the corresponding + * methods on the underlying <code>RandomGenerator</code> and <code>SecureRandom</code> + * instances. Therefore, <code>reSeed(long)</code> fully resets the initial state of the + * non-secure random number generator (so that reseeding with a specific value always results + * in the same subsequent random sequence); whereas reSeedSecure(long) does + * <strong>not</strong> reinitialize the secure random number generator (so secure sequences + * started with calls to reseedSecure(long) won't be identical). + * <li>This implementation is not synchronized. The underlying <code>RandomGenerator</code> or + * <code>SecureRandom</code> instances are not protected by synchronization and are not + * guaranteed to be thread-safe. Therefore, if an instance of this class is concurrently + * utilized by multiple threads, it is the responsibility of client code to synchronize access + * to seeding and data generation methods. + * </ul> + * + * @since 3.1 + */ +public class RandomDataGenerator implements RandomData, Serializable { + + /** Serializable version identifier */ + private static final long serialVersionUID = -626730818244969716L; + + /** underlying random number generator */ + private RandomGenerator rand = null; + + /** underlying secure random number generator */ + private RandomGenerator secRand = null; + + /** + * Construct a RandomDataGenerator, using a default random generator as the source of + * randomness. + * + * <p>The default generator is a {@link Well19937c} seeded with {@code + * System.currentTimeMillis() + System.identityHashCode(this))}. The generator is initialized + * and seeded on first use. + */ + public RandomDataGenerator() {} + + /** + * Construct a RandomDataGenerator using the supplied {@link RandomGenerator} as the source of + * (non-secure) random data. + * + * @param rand the source of (non-secure) random data (may be null, resulting in the default + * generator) + */ + public RandomDataGenerator(RandomGenerator rand) { + this.rand = rand; + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description:</strong> hex strings are generated using a 2-step process. + * + * <ol> + * <li>{@code len / 2 + 1} binary bytes are generated using the underlying Random + * <li>Each binary byte is translated into 2 hex digits + * </ol> + * + * @param len the desired string length. + * @return the random string. + * @throws NotStrictlyPositiveException if {@code len <= 0}. + */ + public String nextHexString(int len) throws NotStrictlyPositiveException { + if (len <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.LENGTH, len); + } + + // Get a random number generator + RandomGenerator ran = getRandomGenerator(); + + // Initialize output buffer + StringBuilder outBuffer = new StringBuilder(); + + // Get int(len/2)+1 random bytes + byte[] randomBytes = new byte[(len / 2) + 1]; + ran.nextBytes(randomBytes); + + // Convert each byte to 2 hex digits + for (int i = 0; i < randomBytes.length; i++) { + Integer c = Integer.valueOf(randomBytes[i]); + + /* + * Add 128 to byte value to make interval 0-255 before doing hex + * conversion. This guarantees <= 2 hex digits from toHexString() + * toHexString would otherwise add 2^32 to negative arguments. + */ + String hex = Integer.toHexString(c.intValue() + 128); + + // Make sure we add 2 hex digits for each byte + if (hex.length() == 1) { + hex = "0" + hex; + } + outBuffer.append(hex); + } + return outBuffer.toString().substring(0, len); + } + + /** {@inheritDoc} */ + public int nextInt(final int lower, final int upper) throws NumberIsTooLargeException { + return new UniformIntegerDistribution(getRandomGenerator(), lower, upper).sample(); + } + + /** {@inheritDoc} */ + public long nextLong(final long lower, final long upper) throws NumberIsTooLargeException { + if (lower >= upper) { + throw new NumberIsTooLargeException( + LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND, lower, upper, false); + } + final long max = (upper - lower) + 1; + if (max <= 0) { + // the range is too wide to fit in a positive long (larger than 2^63); as it covers + // more than half the long range, we use directly a simple rejection method + final RandomGenerator rng = getRandomGenerator(); + while (true) { + final long r = rng.nextLong(); + if (r >= lower && r <= upper) { + return r; + } + } + } else if (max < Integer.MAX_VALUE) { + // we can shift the range and generate directly a positive int + return lower + getRandomGenerator().nextInt((int) max); + } else { + // we can shift the range and generate directly a positive long + return lower + nextLong(getRandomGenerator(), max); + } + } + + /** + * Returns a pseudorandom, uniformly distributed {@code long} value between 0 (inclusive) and + * the specified value (exclusive), drawn from this random number generator's sequence. + * + * @param rng random generator to use + * @param n the bound on the random number to be returned. Must be positive. + * @return a pseudorandom, uniformly distributed {@code long} value between 0 (inclusive) and n + * (exclusive). + * @throws IllegalArgumentException if n is not positive. + */ + private static long nextLong(final RandomGenerator rng, final long n) + throws IllegalArgumentException { + if (n > 0) { + final byte[] byteArray = new byte[8]; + long bits; + long val; + do { + rng.nextBytes(byteArray); + bits = 0; + for (final byte b : byteArray) { + bits = (bits << 8) | (((long) b) & 0xffL); + } + bits &= 0x7fffffffffffffffL; + val = bits % n; + } while (bits - val + (n - 1) < 0); + return val; + } + throw new NotStrictlyPositiveException(n); + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description:</strong> hex strings are generated in 40-byte segments + * using a 3-step process. + * + * <ol> + * <li>20 random bytes are generated using the underlying <code>SecureRandom</code>. + * <li>SHA-1 hash is applied to yield a 20-byte binary digest. + * <li>Each byte of the binary digest is converted to 2 hex digits. + * </ol> + * + * @throws NotStrictlyPositiveException if {@code len <= 0} + */ + public String nextSecureHexString(int len) throws NotStrictlyPositiveException { + if (len <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.LENGTH, len); + } + + // Get SecureRandom and setup Digest provider + final RandomGenerator secRan = getSecRan(); + MessageDigest alg = null; + try { + alg = MessageDigest.getInstance("SHA-1"); + } catch (NoSuchAlgorithmException ex) { + // this should never happen + throw new MathInternalError(ex); + } + alg.reset(); + + // Compute number of iterations required (40 bytes each) + int numIter = (len / 40) + 1; + + StringBuilder outBuffer = new StringBuilder(); + for (int iter = 1; iter < numIter + 1; iter++) { + byte[] randomBytes = new byte[40]; + secRan.nextBytes(randomBytes); + alg.update(randomBytes); + + // Compute hash -- will create 20-byte binary hash + byte[] hash = alg.digest(); + + // Loop over the hash, converting each byte to 2 hex digits + for (int i = 0; i < hash.length; i++) { + Integer c = Integer.valueOf(hash[i]); + + /* + * Add 128 to byte value to make interval 0-255 This guarantees + * <= 2 hex digits from toHexString() toHexString would + * otherwise add 2^32 to negative arguments + */ + String hex = Integer.toHexString(c.intValue() + 128); + + // Keep strings uniform length -- guarantees 40 bytes + if (hex.length() == 1) { + hex = "0" + hex; + } + outBuffer.append(hex); + } + } + return outBuffer.toString().substring(0, len); + } + + /** {@inheritDoc} */ + public int nextSecureInt(final int lower, final int upper) throws NumberIsTooLargeException { + return new UniformIntegerDistribution(getSecRan(), lower, upper).sample(); + } + + /** {@inheritDoc} */ + public long nextSecureLong(final long lower, final long upper) + throws NumberIsTooLargeException { + if (lower >= upper) { + throw new NumberIsTooLargeException( + LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND, lower, upper, false); + } + final RandomGenerator rng = getSecRan(); + final long max = (upper - lower) + 1; + if (max <= 0) { + // the range is too wide to fit in a positive long (larger than 2^63); as it covers + // more than half the long range, we use directly a simple rejection method + while (true) { + final long r = rng.nextLong(); + if (r >= lower && r <= upper) { + return r; + } + } + } else if (max < Integer.MAX_VALUE) { + // we can shift the range and generate directly a positive int + return lower + rng.nextInt((int) max); + } else { + // we can shift the range and generate directly a positive long + return lower + nextLong(rng, max); + } + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description</strong>: + * + * <ul> + * <li>For small means, uses simulation of a Poisson process using Uniform deviates, as + * described <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm">here.</a> The + * Poisson process (and hence value returned) is bounded by 1000 * mean. + * <li>For large means, uses the rejection algorithm described in <br> + * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i> + * <strong>Computing</strong> vol. 26 pp. 197-207. + * </ul> + * + * @throws NotStrictlyPositiveException if {@code len <= 0} + */ + public long nextPoisson(double mean) throws NotStrictlyPositiveException { + return new PoissonDistribution( + getRandomGenerator(), + mean, + PoissonDistribution.DEFAULT_EPSILON, + PoissonDistribution.DEFAULT_MAX_ITERATIONS) + .sample(); + } + + /** {@inheritDoc} */ + public double nextGaussian(double mu, double sigma) throws NotStrictlyPositiveException { + if (sigma <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.STANDARD_DEVIATION, sigma); + } + return sigma * getRandomGenerator().nextGaussian() + mu; + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description</strong>: Uses the Algorithm SA (Ahrens) from p. 876 in: + * [1]: Ahrens, J. H. and Dieter, U. (1972). Computer methods for sampling from the exponential + * and normal distributions. Communications of the ACM, 15, 873-882. + */ + public double nextExponential(double mean) throws NotStrictlyPositiveException { + return new ExponentialDistribution( + getRandomGenerator(), + mean, + ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY) + .sample(); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.GammaDistribution Gamma Distribution}. + * + * <p>This implementation uses the following algorithms: + * + * <p>For 0 < shape < 1: <br> + * Ahrens, J. H. and Dieter, U., <i>Computer methods for sampling from gamma, beta, Poisson and + * binomial distributions.</i> Computing, 12, 223-246, 1974. + * + * <p>For shape >= 1: <br> + * Marsaglia and Tsang, <i>A Simple Method for Generating Gamma Variables.</i> ACM Transactions + * on Mathematical Software, Volume 26 Issue 3, September, 2000. + * + * @param shape the median of the Gamma distribution + * @param scale the scale parameter of the Gamma distribution + * @return random value sampled from the Gamma(shape, scale) distribution + * @throws NotStrictlyPositiveException if {@code shape <= 0} or {@code scale <= 0}. + */ + public double nextGamma(double shape, double scale) throws NotStrictlyPositiveException { + return new GammaDistribution( + getRandomGenerator(), + shape, + scale, + GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY) + .sample(); + } + + /** + * Generates a random value from the {@link HypergeometricDistribution Hypergeometric + * Distribution}. + * + * @param populationSize the population size of the Hypergeometric distribution + * @param numberOfSuccesses number of successes in the population of the Hypergeometric + * distribution + * @param sampleSize the sample size of the Hypergeometric distribution + * @return random value sampled from the Hypergeometric(numberOfSuccesses, sampleSize) + * distribution + * @throws NumberIsTooLargeException if {@code numberOfSuccesses > populationSize}, or {@code + * sampleSize > populationSize}. + * @throws NotStrictlyPositiveException if {@code populationSize <= 0}. + * @throws NotPositiveException if {@code numberOfSuccesses < 0}. + */ + public int nextHypergeometric(int populationSize, int numberOfSuccesses, int sampleSize) + throws NotPositiveException, NotStrictlyPositiveException, NumberIsTooLargeException { + return new HypergeometricDistribution( + getRandomGenerator(), populationSize, numberOfSuccesses, sampleSize) + .sample(); + } + + /** + * Generates a random value from the {@link PascalDistribution Pascal Distribution}. + * + * @param r the number of successes of the Pascal distribution + * @param p the probability of success of the Pascal distribution + * @return random value sampled from the Pascal(r, p) distribution + * @throws NotStrictlyPositiveException if the number of successes is not positive + * @throws OutOfRangeException if the probability of success is not in the range {@code [0, 1]}. + */ + public int nextPascal(int r, double p) + throws NotStrictlyPositiveException, OutOfRangeException { + return new PascalDistribution(getRandomGenerator(), r, p).sample(); + } + + /** + * Generates a random value from the {@link TDistribution T Distribution}. + * + * @param df the degrees of freedom of the T distribution + * @return random value from the T(df) distribution + * @throws NotStrictlyPositiveException if {@code df <= 0} + */ + public double nextT(double df) throws NotStrictlyPositiveException { + return new TDistribution( + getRandomGenerator(), df, TDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY) + .sample(); + } + + /** + * Generates a random value from the {@link WeibullDistribution Weibull Distribution}. + * + * @param shape the shape parameter of the Weibull distribution + * @param scale the scale parameter of the Weibull distribution + * @return random value sampled from the Weibull(shape, size) distribution + * @throws NotStrictlyPositiveException if {@code shape <= 0} or {@code scale <= 0}. + */ + public double nextWeibull(double shape, double scale) throws NotStrictlyPositiveException { + return new WeibullDistribution( + getRandomGenerator(), + shape, + scale, + WeibullDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY) + .sample(); + } + + /** + * Generates a random value from the {@link ZipfDistribution Zipf Distribution}. + * + * @param numberOfElements the number of elements of the ZipfDistribution + * @param exponent the exponent of the ZipfDistribution + * @return random value sampled from the Zipf(numberOfElements, exponent) distribution + * @exception NotStrictlyPositiveException if {@code numberOfElements <= 0} or {@code exponent + * <= 0}. + */ + public int nextZipf(int numberOfElements, double exponent) throws NotStrictlyPositiveException { + return new ZipfDistribution(getRandomGenerator(), numberOfElements, exponent).sample(); + } + + /** + * Generates a random value from the {@link BetaDistribution Beta Distribution}. + * + * @param alpha first distribution shape parameter + * @param beta second distribution shape parameter + * @return random value sampled from the beta(alpha, beta) distribution + */ + public double nextBeta(double alpha, double beta) { + return new BetaDistribution( + getRandomGenerator(), + alpha, + beta, + BetaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY) + .sample(); + } + + /** + * Generates a random value from the {@link BinomialDistribution Binomial Distribution}. + * + * @param numberOfTrials number of trials of the Binomial distribution + * @param probabilityOfSuccess probability of success of the Binomial distribution + * @return random value sampled from the Binomial(numberOfTrials, probabilityOfSuccess) + * distribution + */ + public int nextBinomial(int numberOfTrials, double probabilityOfSuccess) { + return new BinomialDistribution(getRandomGenerator(), numberOfTrials, probabilityOfSuccess) + .sample(); + } + + /** + * Generates a random value from the {@link CauchyDistribution Cauchy Distribution}. + * + * @param median the median of the Cauchy distribution + * @param scale the scale parameter of the Cauchy distribution + * @return random value sampled from the Cauchy(median, scale) distribution + */ + public double nextCauchy(double median, double scale) { + return new CauchyDistribution( + getRandomGenerator(), + median, + scale, + CauchyDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY) + .sample(); + } + + /** + * Generates a random value from the {@link ChiSquaredDistribution ChiSquare Distribution}. + * + * @param df the degrees of freedom of the ChiSquare distribution + * @return random value sampled from the ChiSquare(df) distribution + */ + public double nextChiSquare(double df) { + return new ChiSquaredDistribution( + getRandomGenerator(), + df, + ChiSquaredDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY) + .sample(); + } + + /** + * Generates a random value from the {@link FDistribution F Distribution}. + * + * @param numeratorDf the numerator degrees of freedom of the F distribution + * @param denominatorDf the denominator degrees of freedom of the F distribution + * @return random value sampled from the F(numeratorDf, denominatorDf) distribution + * @throws NotStrictlyPositiveException if {@code numeratorDf <= 0} or {@code denominatorDf <= + * 0}. + */ + public double nextF(double numeratorDf, double denominatorDf) + throws NotStrictlyPositiveException { + return new FDistribution( + getRandomGenerator(), + numeratorDf, + denominatorDf, + FDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY) + .sample(); + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description</strong>: scales the output of Random.nextDouble(), but + * rejects 0 values (i.e., will generate another random double if Random.nextDouble() returns + * 0). This is necessary to provide a symmetric output interval (both endpoints excluded). + * + * @throws NumberIsTooLargeException if {@code lower >= upper} + * @throws NotFiniteNumberException if one of the bounds is infinite + * @throws NotANumberException if one of the bounds is NaN + */ + public double nextUniform(double lower, double upper) + throws NumberIsTooLargeException, NotFiniteNumberException, NotANumberException { + return nextUniform(lower, upper, false); + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description</strong>: if the lower bound is excluded, scales the output + * of Random.nextDouble(), but rejects 0 values (i.e., will generate another random double if + * Random.nextDouble() returns 0). This is necessary to provide a symmetric output interval + * (both endpoints excluded). + * + * @throws NumberIsTooLargeException if {@code lower >= upper} + * @throws NotFiniteNumberException if one of the bounds is infinite + * @throws NotANumberException if one of the bounds is NaN + */ + public double nextUniform(double lower, double upper, boolean lowerInclusive) + throws NumberIsTooLargeException, NotFiniteNumberException, NotANumberException { + + if (lower >= upper) { + throw new NumberIsTooLargeException( + LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND, lower, upper, false); + } + + if (Double.isInfinite(lower)) { + throw new NotFiniteNumberException(LocalizedFormats.INFINITE_BOUND, lower); + } + if (Double.isInfinite(upper)) { + throw new NotFiniteNumberException(LocalizedFormats.INFINITE_BOUND, upper); + } + + if (Double.isNaN(lower) || Double.isNaN(upper)) { + throw new NotANumberException(); + } + + final RandomGenerator generator = getRandomGenerator(); + + // ensure nextDouble() isn't 0.0 + double u = generator.nextDouble(); + while (!lowerInclusive && u <= 0.0) { + u = generator.nextDouble(); + } + + return u * upper + (1.0 - u) * lower; + } + + /** + * {@inheritDoc} + * + * <p>This method calls {@link MathArrays#shuffle(int[],RandomGenerator) MathArrays.shuffle} in + * order to create a random shuffle of the set of natural numbers {@code { 0, 1, ..., n - 1 }}. + * + * @throws NumberIsTooLargeException if {@code k > n}. + * @throws NotStrictlyPositiveException if {@code k <= 0}. + */ + public int[] nextPermutation(int n, int k) + throws NumberIsTooLargeException, NotStrictlyPositiveException { + if (k > n) { + throw new NumberIsTooLargeException(LocalizedFormats.PERMUTATION_EXCEEDS_N, k, n, true); + } + if (k <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.PERMUTATION_SIZE, k); + } + + int[] index = MathArrays.natural(n); + MathArrays.shuffle(index, getRandomGenerator()); + + // Return a new array containing the first "k" entries of "index". + return MathArrays.copyOf(index, k); + } + + /** + * {@inheritDoc} + * + * <p>This method calls {@link #nextPermutation(int,int) nextPermutation(c.size(), k)} in order + * to sample the collection. + */ + public Object[] nextSample(Collection<?> c, int k) + throws NumberIsTooLargeException, NotStrictlyPositiveException { + + int len = c.size(); + if (k > len) { + throw new NumberIsTooLargeException( + LocalizedFormats.SAMPLE_SIZE_EXCEEDS_COLLECTION_SIZE, k, len, true); + } + if (k <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES, k); + } + + Object[] objects = c.toArray(); + int[] index = nextPermutation(len, k); + Object[] result = new Object[k]; + for (int i = 0; i < k; i++) { + result[i] = objects[index[i]]; + } + return result; + } + + /** + * Reseeds the random number generator with the supplied seed. + * + * <p>Will create and initialize if null. + * + * @param seed the seed value to use + */ + public void reSeed(long seed) { + getRandomGenerator().setSeed(seed); + } + + /** + * Reseeds the secure random number generator with the current time in milliseconds. + * + * <p>Will create and initialize if null. + */ + public void reSeedSecure() { + getSecRan().setSeed(System.currentTimeMillis()); + } + + /** + * Reseeds the secure random number generator with the supplied seed. + * + * <p>Will create and initialize if null. + * + * @param seed the seed value to use + */ + public void reSeedSecure(long seed) { + getSecRan().setSeed(seed); + } + + /** + * Reseeds the random number generator with {@code System.currentTimeMillis() + + * System.identityHashCode(this))}. + */ + public void reSeed() { + getRandomGenerator().setSeed(System.currentTimeMillis() + System.identityHashCode(this)); + } + + /** + * Sets the PRNG algorithm for the underlying SecureRandom instance using the Security Provider + * API. The Security Provider API is defined in <a href = + * "http://java.sun.com/j2se/1.3/docs/guide/security/CryptoSpec.html#AppA"> Java Cryptography + * Architecture API Specification & Reference.</a> + * + * <p><strong>USAGE NOTE:</strong> This method carries <i>significant</i> overhead and may take + * several seconds to execute. + * + * @param algorithm the name of the PRNG algorithm + * @param provider the name of the provider + * @throws NoSuchAlgorithmException if the specified algorithm is not available + * @throws NoSuchProviderException if the specified provider is not installed + */ + public void setSecureAlgorithm(String algorithm, String provider) + throws NoSuchAlgorithmException, NoSuchProviderException { + secRand = + RandomGeneratorFactory.createRandomGenerator( + SecureRandom.getInstance(algorithm, provider)); + } + + /** + * Returns the RandomGenerator used to generate non-secure random data. + * + * <p>Creates and initializes a default generator if null. Uses a {@link Well19937c} generator + * with {@code System.currentTimeMillis() + System.identityHashCode(this))} as the default seed. + * + * @return the Random used to generate random data + * @since 3.2 + */ + public RandomGenerator getRandomGenerator() { + if (rand == null) { + initRan(); + } + return rand; + } + + /** + * Sets the default generator to a {@link Well19937c} generator seeded with {@code + * System.currentTimeMillis() + System.identityHashCode(this))}. + */ + private void initRan() { + rand = new Well19937c(System.currentTimeMillis() + System.identityHashCode(this)); + } + + /** + * Returns the SecureRandom used to generate secure random data. + * + * <p>Creates and initializes if null. Uses {@code System.currentTimeMillis() + + * System.identityHashCode(this)} as the default seed. + * + * @return the SecureRandom used to generate secure random data, wrapped in a {@link + * RandomGenerator}. + */ + private RandomGenerator getSecRan() { + if (secRand == null) { + secRand = RandomGeneratorFactory.createRandomGenerator(new SecureRandom()); + secRand.setSeed(System.currentTimeMillis() + System.identityHashCode(this)); + } + return secRand; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/RandomDataImpl.java b/src/main/java/org/apache/commons/math3/random/RandomDataImpl.java new file mode 100644 index 0000000..d5749e9 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/RandomDataImpl.java @@ -0,0 +1,544 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.distribution.IntegerDistribution; +import org.apache.commons.math3.distribution.RealDistribution; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.NotANumberException; +import org.apache.commons.math3.exception.NotFiniteNumberException; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.OutOfRangeException; + +import java.io.Serializable; +import java.security.NoSuchAlgorithmException; +import java.security.NoSuchProviderException; +import java.util.Collection; + +/** + * Generates random deviates and other random data using a {@link RandomGenerator} instance to + * generate non-secure data and a {@link java.security.SecureRandom} instance to provide data for + * the <code>nextSecureXxx</code> methods. If no <code>RandomGenerator</code> is provided in the + * constructor, the default is to use a {@link Well19937c} generator. To plug in a different + * implementation, either implement <code>RandomGenerator</code> directly or extend {@link + * AbstractRandomGenerator}. + * + * <p>Supports reseeding the underlying pseudo-random number generator (PRNG). The <code> + * SecurityProvider</code> and <code>Algorithm</code> used by the <code>SecureRandom</code> instance + * can also be reset. + * + * <p>For details on the default PRNGs, see {@link java.util.Random} and {@link + * java.security.SecureRandom}. + * + * <p><strong>Usage Notes</strong>: + * + * <ul> + * <li>Instance variables are used to maintain <code>RandomGenerator</code> and <code>SecureRandom + * </code> instances used in data generation. Therefore, to generate a random sequence of + * values or strings, you should use just <strong>one</strong> <code>RandomDataGenerator + * </code> instance repeatedly. + * <li>The "secure" methods are *much* slower. These should be used only when a cryptographically + * secure random sequence is required. A secure random sequence is a sequence of pseudo-random + * values which, in addition to being well-dispersed (so no subsequence of values is an any + * more likely than other subsequence of the the same length), also has the additional + * property that knowledge of values generated up to any point in the sequence does not make + * it any easier to predict subsequent values. + * <li>When a new <code>RandomDataGenerator</code> is created, the underlying random number + * generators are <strong>not</strong> initialized. If you do not explicitly seed the default + * non-secure generator, it is seeded with the current time in milliseconds plus the system + * identity hash code on first use. The same holds for the secure generator. If you provide a + * <code>RandomGenerator</code> to the constructor, however, this generator is not reseeded by + * the constructor nor is it reseeded on first use. + * <li>The <code>reSeed</code> and <code>reSeedSecure</code> methods delegate to the corresponding + * methods on the underlying <code>RandomGenerator</code> and <code>SecureRandom</code> + * instances. Therefore, <code>reSeed(long)</code> fully resets the initial state of the + * non-secure random number generator (so that reseeding with a specific value always results + * in the same subsequent random sequence); whereas reSeedSecure(long) does + * <strong>not</strong> reinitialize the secure random number generator (so secure sequences + * started with calls to reseedSecure(long) won't be identical). + * <li>This implementation is not synchronized. The underlying <code>RandomGenerator</code> or + * <code>SecureRandom</code> instances are not protected by synchronization and are not + * guaranteed to be thread-safe. Therefore, if an instance of this class is concurrently + * utilized by multiple threads, it is the responsibility of client code to synchronize access + * to seeding and data generation methods. + * </ul> + * + * @deprecated to be removed in 4.0. Use {@link RandomDataGenerator} instead + */ +@Deprecated +public class RandomDataImpl implements RandomData, Serializable { + + /** Serializable version identifier */ + private static final long serialVersionUID = -626730818244969716L; + + /** RandomDataGenerator delegate */ + private final RandomDataGenerator delegate; + + /** + * Construct a RandomDataImpl, using a default random generator as the source of randomness. + * + * <p>The default generator is a {@link Well19937c} seeded with {@code + * System.currentTimeMillis() + System.identityHashCode(this))}. The generator is initialized + * and seeded on first use. + */ + public RandomDataImpl() { + delegate = new RandomDataGenerator(); + } + + /** + * Construct a RandomDataImpl using the supplied {@link RandomGenerator} as the source of + * (non-secure) random data. + * + * @param rand the source of (non-secure) random data (may be null, resulting in the default + * generator) + * @since 1.1 + */ + public RandomDataImpl(RandomGenerator rand) { + delegate = new RandomDataGenerator(rand); + } + + /** + * @return the delegate object. + * @deprecated To be removed in 4.0. + */ + @Deprecated + RandomDataGenerator getDelegate() { + return delegate; + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description:</strong> hex strings are generated using a 2-step process. + * + * <ol> + * <li>{@code len / 2 + 1} binary bytes are generated using the underlying Random + * <li>Each binary byte is translated into 2 hex digits + * </ol> + * + * @param len the desired string length. + * @return the random string. + * @throws NotStrictlyPositiveException if {@code len <= 0}. + */ + public String nextHexString(int len) throws NotStrictlyPositiveException { + return delegate.nextHexString(len); + } + + /** {@inheritDoc} */ + public int nextInt(int lower, int upper) throws NumberIsTooLargeException { + return delegate.nextInt(lower, upper); + } + + /** {@inheritDoc} */ + public long nextLong(long lower, long upper) throws NumberIsTooLargeException { + return delegate.nextLong(lower, upper); + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description:</strong> hex strings are generated in 40-byte segments + * using a 3-step process. + * + * <ol> + * <li>20 random bytes are generated using the underlying <code>SecureRandom</code>. + * <li>SHA-1 hash is applied to yield a 20-byte binary digest. + * <li>Each byte of the binary digest is converted to 2 hex digits. + * </ol> + */ + public String nextSecureHexString(int len) throws NotStrictlyPositiveException { + return delegate.nextSecureHexString(len); + } + + /** {@inheritDoc} */ + public int nextSecureInt(int lower, int upper) throws NumberIsTooLargeException { + return delegate.nextSecureInt(lower, upper); + } + + /** {@inheritDoc} */ + public long nextSecureLong(long lower, long upper) throws NumberIsTooLargeException { + return delegate.nextSecureLong(lower, upper); + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description</strong>: + * + * <ul> + * <li>For small means, uses simulation of a Poisson process using Uniform deviates, as + * described <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm">here.</a> The + * Poisson process (and hence value returned) is bounded by 1000 * mean. + * <li>For large means, uses the rejection algorithm described in <br> + * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i> + * <strong>Computing</strong> vol. 26 pp. 197-207. + * </ul> + */ + public long nextPoisson(double mean) throws NotStrictlyPositiveException { + return delegate.nextPoisson(mean); + } + + /** {@inheritDoc} */ + public double nextGaussian(double mu, double sigma) throws NotStrictlyPositiveException { + return delegate.nextGaussian(mu, sigma); + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description</strong>: Uses the Algorithm SA (Ahrens) from p. 876 in: + * [1]: Ahrens, J. H. and Dieter, U. (1972). Computer methods for sampling from the exponential + * and normal distributions. Communications of the ACM, 15, 873-882. + */ + public double nextExponential(double mean) throws NotStrictlyPositiveException { + return delegate.nextExponential(mean); + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description</strong>: scales the output of Random.nextDouble(), but + * rejects 0 values (i.e., will generate another random double if Random.nextDouble() returns + * 0). This is necessary to provide a symmetric output interval (both endpoints excluded). + */ + public double nextUniform(double lower, double upper) + throws NumberIsTooLargeException, NotFiniteNumberException, NotANumberException { + return delegate.nextUniform(lower, upper); + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description</strong>: if the lower bound is excluded, scales the output + * of Random.nextDouble(), but rejects 0 values (i.e., will generate another random double if + * Random.nextDouble() returns 0). This is necessary to provide a symmetric output interval + * (both endpoints excluded). + * + * @since 3.0 + */ + public double nextUniform(double lower, double upper, boolean lowerInclusive) + throws NumberIsTooLargeException, NotFiniteNumberException, NotANumberException { + return delegate.nextUniform(lower, upper, lowerInclusive); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.BetaDistribution Beta Distribution}. This + * implementation uses {@link #nextInversionDeviate(RealDistribution) inversion} to generate + * random values. + * + * @param alpha first distribution shape parameter + * @param beta second distribution shape parameter + * @return random value sampled from the beta(alpha, beta) distribution + * @since 2.2 + */ + public double nextBeta(double alpha, double beta) { + return delegate.nextBeta(alpha, beta); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.BinomialDistribution Binomial Distribution}. This + * implementation uses {@link #nextInversionDeviate(RealDistribution) inversion} to generate + * random values. + * + * @param numberOfTrials number of trials of the Binomial distribution + * @param probabilityOfSuccess probability of success of the Binomial distribution + * @return random value sampled from the Binomial(numberOfTrials, probabilityOfSuccess) + * distribution + * @since 2.2 + */ + public int nextBinomial(int numberOfTrials, double probabilityOfSuccess) { + return delegate.nextBinomial(numberOfTrials, probabilityOfSuccess); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.CauchyDistribution Cauchy Distribution}. This + * implementation uses {@link #nextInversionDeviate(RealDistribution) inversion} to generate + * random values. + * + * @param median the median of the Cauchy distribution + * @param scale the scale parameter of the Cauchy distribution + * @return random value sampled from the Cauchy(median, scale) distribution + * @since 2.2 + */ + public double nextCauchy(double median, double scale) { + return delegate.nextCauchy(median, scale); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.ChiSquaredDistribution ChiSquare Distribution}. This + * implementation uses {@link #nextInversionDeviate(RealDistribution) inversion} to generate + * random values. + * + * @param df the degrees of freedom of the ChiSquare distribution + * @return random value sampled from the ChiSquare(df) distribution + * @since 2.2 + */ + public double nextChiSquare(double df) { + return delegate.nextChiSquare(df); + } + + /** + * Generates a random value from the {@link org.apache.commons.math3.distribution.FDistribution + * F Distribution}. This implementation uses {@link #nextInversionDeviate(RealDistribution) + * inversion} to generate random values. + * + * @param numeratorDf the numerator degrees of freedom of the F distribution + * @param denominatorDf the denominator degrees of freedom of the F distribution + * @return random value sampled from the F(numeratorDf, denominatorDf) distribution + * @throws NotStrictlyPositiveException if {@code numeratorDf <= 0} or {@code denominatorDf <= + * 0}. + * @since 2.2 + */ + public double nextF(double numeratorDf, double denominatorDf) + throws NotStrictlyPositiveException { + return delegate.nextF(numeratorDf, denominatorDf); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.GammaDistribution Gamma Distribution}. + * + * <p>This implementation uses the following algorithms: + * + * <p>For 0 < shape < 1: <br> + * Ahrens, J. H. and Dieter, U., <i>Computer methods for sampling from gamma, beta, Poisson and + * binomial distributions.</i> Computing, 12, 223-246, 1974. + * + * <p>For shape >= 1: <br> + * Marsaglia and Tsang, <i>A Simple Method for Generating Gamma Variables.</i> ACM Transactions + * on Mathematical Software, Volume 26 Issue 3, September, 2000. + * + * @param shape the median of the Gamma distribution + * @param scale the scale parameter of the Gamma distribution + * @return random value sampled from the Gamma(shape, scale) distribution + * @throws NotStrictlyPositiveException if {@code shape <= 0} or {@code scale <= 0}. + * @since 2.2 + */ + public double nextGamma(double shape, double scale) throws NotStrictlyPositiveException { + return delegate.nextGamma(shape, scale); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.HypergeometricDistribution Hypergeometric + * Distribution}. This implementation uses {@link #nextInversionDeviate(IntegerDistribution) + * inversion} to generate random values. + * + * @param populationSize the population size of the Hypergeometric distribution + * @param numberOfSuccesses number of successes in the population of the Hypergeometric + * distribution + * @param sampleSize the sample size of the Hypergeometric distribution + * @return random value sampled from the Hypergeometric(numberOfSuccesses, sampleSize) + * distribution + * @throws NumberIsTooLargeException if {@code numberOfSuccesses > populationSize}, or {@code + * sampleSize > populationSize}. + * @throws NotStrictlyPositiveException if {@code populationSize <= 0}. + * @throws NotPositiveException if {@code numberOfSuccesses < 0}. + * @since 2.2 + */ + public int nextHypergeometric(int populationSize, int numberOfSuccesses, int sampleSize) + throws NotPositiveException, NotStrictlyPositiveException, NumberIsTooLargeException { + return delegate.nextHypergeometric(populationSize, numberOfSuccesses, sampleSize); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.PascalDistribution Pascal Distribution}. This + * implementation uses {@link #nextInversionDeviate(IntegerDistribution) inversion} to generate + * random values. + * + * @param r the number of successes of the Pascal distribution + * @param p the probability of success of the Pascal distribution + * @return random value sampled from the Pascal(r, p) distribution + * @since 2.2 + * @throws NotStrictlyPositiveException if the number of successes is not positive + * @throws OutOfRangeException if the probability of success is not in the range {@code [0, 1]}. + */ + public int nextPascal(int r, double p) + throws NotStrictlyPositiveException, OutOfRangeException { + return delegate.nextPascal(r, p); + } + + /** + * Generates a random value from the {@link org.apache.commons.math3.distribution.TDistribution + * T Distribution}. This implementation uses {@link #nextInversionDeviate(RealDistribution) + * inversion} to generate random values. + * + * @param df the degrees of freedom of the T distribution + * @return random value from the T(df) distribution + * @since 2.2 + * @throws NotStrictlyPositiveException if {@code df <= 0} + */ + public double nextT(double df) throws NotStrictlyPositiveException { + return delegate.nextT(df); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.WeibullDistribution Weibull Distribution}. This + * implementation uses {@link #nextInversionDeviate(RealDistribution) inversion} to generate + * random values. + * + * @param shape the shape parameter of the Weibull distribution + * @param scale the scale parameter of the Weibull distribution + * @return random value sampled from the Weibull(shape, size) distribution + * @since 2.2 + * @throws NotStrictlyPositiveException if {@code shape <= 0} or {@code scale <= 0}. + */ + public double nextWeibull(double shape, double scale) throws NotStrictlyPositiveException { + return delegate.nextWeibull(shape, scale); + } + + /** + * Generates a random value from the {@link + * org.apache.commons.math3.distribution.ZipfDistribution Zipf Distribution}. This + * implementation uses {@link #nextInversionDeviate(IntegerDistribution) inversion} to generate + * random values. + * + * @param numberOfElements the number of elements of the ZipfDistribution + * @param exponent the exponent of the ZipfDistribution + * @return random value sampled from the Zipf(numberOfElements, exponent) distribution + * @since 2.2 + * @exception NotStrictlyPositiveException if {@code numberOfElements <= 0} or {@code exponent + * <= 0}. + */ + public int nextZipf(int numberOfElements, double exponent) throws NotStrictlyPositiveException { + return delegate.nextZipf(numberOfElements, exponent); + } + + /** + * Reseeds the random number generator with the supplied seed. + * + * <p>Will create and initialize if null. + * + * @param seed the seed value to use + */ + public void reSeed(long seed) { + delegate.reSeed(seed); + } + + /** + * Reseeds the secure random number generator with the current time in milliseconds. + * + * <p>Will create and initialize if null. + */ + public void reSeedSecure() { + delegate.reSeedSecure(); + } + + /** + * Reseeds the secure random number generator with the supplied seed. + * + * <p>Will create and initialize if null. + * + * @param seed the seed value to use + */ + public void reSeedSecure(long seed) { + delegate.reSeedSecure(seed); + } + + /** + * Reseeds the random number generator with {@code System.currentTimeMillis() + + * System.identityHashCode(this))}. + */ + public void reSeed() { + delegate.reSeed(); + } + + /** + * Sets the PRNG algorithm for the underlying SecureRandom instance using the Security Provider + * API. The Security Provider API is defined in <a href = + * "http://java.sun.com/j2se/1.3/docs/guide/security/CryptoSpec.html#AppA"> Java Cryptography + * Architecture API Specification & Reference.</a> + * + * <p><strong>USAGE NOTE:</strong> This method carries <i>significant</i> overhead and may take + * several seconds to execute. + * + * @param algorithm the name of the PRNG algorithm + * @param provider the name of the provider + * @throws NoSuchAlgorithmException if the specified algorithm is not available + * @throws NoSuchProviderException if the specified provider is not installed + */ + public void setSecureAlgorithm(String algorithm, String provider) + throws NoSuchAlgorithmException, NoSuchProviderException { + delegate.setSecureAlgorithm(algorithm, provider); + } + + /** + * {@inheritDoc} + * + * <p>Uses a 2-cycle permutation shuffle. The shuffling process is described <a + * href="http://www.maths.abdn.ac.uk/~igc/tch/mx4002/notes/node83.html">here</a>. + */ + public int[] nextPermutation(int n, int k) + throws NotStrictlyPositiveException, NumberIsTooLargeException { + return delegate.nextPermutation(n, k); + } + + /** + * {@inheritDoc} + * + * <p><strong>Algorithm Description</strong>: Uses a 2-cycle permutation shuffle to generate a + * random permutation of <code>c.size()</code> and then returns the elements whose indexes + * correspond to the elements of the generated permutation. This technique is described, and + * proven to generate random samples <a + * href="http://www.maths.abdn.ac.uk/~igc/tch/mx4002/notes/node83.html">here</a> + */ + public Object[] nextSample(Collection<?> c, int k) + throws NotStrictlyPositiveException, NumberIsTooLargeException { + return delegate.nextSample(c, k); + } + + /** + * Generate a random deviate from the given distribution using the <a + * href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">inversion method.</a> + * + * @param distribution Continuous distribution to generate a random value from + * @return a random value sampled from the given distribution + * @throws MathIllegalArgumentException if the underlynig distribution throws one + * @since 2.2 + * @deprecated use the distribution's sample() method + */ + @Deprecated + public double nextInversionDeviate(RealDistribution distribution) + throws MathIllegalArgumentException { + return distribution.inverseCumulativeProbability(nextUniform(0, 1)); + } + + /** + * Generate a random deviate from the given distribution using the <a + * href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">inversion method.</a> + * + * @param distribution Integer distribution to generate a random value from + * @return a random value sampled from the given distribution + * @throws MathIllegalArgumentException if the underlynig distribution throws one + * @since 2.2 + * @deprecated use the distribution's sample() method + */ + @Deprecated + public int nextInversionDeviate(IntegerDistribution distribution) + throws MathIllegalArgumentException { + return distribution.inverseCumulativeProbability(nextUniform(0, 1)); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/RandomGenerator.java b/src/main/java/org/apache/commons/math3/random/RandomGenerator.java new file mode 100644 index 0000000..aeba5e3 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/RandomGenerator.java @@ -0,0 +1,130 @@ +/* + * 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.math3.random; + +/** + * Interface extracted from <code>java.util.Random</code>. This interface is implemented by {@link + * AbstractRandomGenerator}. + * + * @since 1.1 + */ +public interface RandomGenerator { + + /** + * Sets the seed of the underlying random number generator using an <code>int</code> seed. + * + * <p>Sequences of values generated starting with the same seeds should be identical. + * + * @param seed the seed value + */ + void setSeed(int seed); + + /** + * Sets the seed of the underlying random number generator using an <code>int</code> array seed. + * + * <p>Sequences of values generated starting with the same seeds should be identical. + * + * @param seed the seed value + */ + void setSeed(int[] seed); + + /** + * Sets the seed of the underlying random number generator using a <code>long</code> seed. + * + * <p>Sequences of values generated starting with the same seeds should be identical. + * + * @param seed the seed value + */ + void setSeed(long seed); + + /** + * Generates random bytes and places them into a user-supplied byte array. The number of random + * bytes produced is equal to the length of the byte array. + * + * @param bytes the non-null byte array in which to put the random bytes + */ + void nextBytes(byte[] bytes); + + /** + * Returns the next pseudorandom, uniformly distributed <code>int</code> value from this random + * number generator's sequence. All 2<font size="-1"><sup>32</sup></font> possible {@code int} + * values should be produced with (approximately) equal probability. + * + * @return the next pseudorandom, uniformly distributed <code>int</code> value from this random + * number generator's sequence + */ + int nextInt(); + + /** + * Returns a pseudorandom, uniformly distributed {@code int} value between 0 (inclusive) and the + * specified value (exclusive), drawn from this random number generator's sequence. + * + * @param n the bound on the random number to be returned. Must be positive. + * @return a pseudorandom, uniformly distributed {@code int} value between 0 (inclusive) and n + * (exclusive). + * @throws IllegalArgumentException if n is not positive. + */ + int nextInt(int n); + + /** + * Returns the next pseudorandom, uniformly distributed <code>long</code> value from this random + * number generator's sequence. All 2<font size="-1"><sup>64</sup></font> possible {@code long} + * values should be produced with (approximately) equal probability. + * + * @return the next pseudorandom, uniformly distributed <code>long</code> value from this random + * number generator's sequence + */ + long nextLong(); + + /** + * Returns the next pseudorandom, uniformly distributed <code>boolean</code> value from this + * random number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed <code>boolean</code> value from this + * random number generator's sequence + */ + boolean nextBoolean(); + + /** + * Returns the next pseudorandom, uniformly distributed <code>float</code> value between <code> + * 0.0</code> and <code>1.0</code> from this random number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed <code>float</code> value between <code> + * 0.0</code> and <code>1.0</code> from this random number generator's sequence + */ + float nextFloat(); + + /** + * Returns the next pseudorandom, uniformly distributed <code>double</code> value between <code> + * 0.0</code> and <code>1.0</code> from this random number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed <code>double</code> value between <code> + * 0.0</code> and <code>1.0</code> from this random number generator's sequence + */ + double nextDouble(); + + /** + * Returns the next pseudorandom, Gaussian ("normally") distributed <code>double</code> value + * with mean <code>0.0</code> and standard deviation <code>1.0</code> from this random number + * generator's sequence. + * + * @return the next pseudorandom, Gaussian ("normally") distributed <code>double</code> value + * with mean <code>0.0</code> and standard deviation <code>1.0</code> from this random + * number generator's sequence + */ + double nextGaussian(); +} diff --git a/src/main/java/org/apache/commons/math3/random/RandomGeneratorFactory.java b/src/main/java/org/apache/commons/math3/random/RandomGeneratorFactory.java new file mode 100644 index 0000000..b610b44 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/RandomGeneratorFactory.java @@ -0,0 +1,118 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.NotStrictlyPositiveException; + +import java.util.Random; + +/** + * Utilities for creating {@link RandomGenerator} instances. + * + * @since 3.3 + */ +public class RandomGeneratorFactory { + /** Class contains only static methods. */ + private RandomGeneratorFactory() {} + + /** + * Creates a {@link RandomDataGenerator} instance that wraps a {@link Random} instance. + * + * @param rng JDK {@link Random} instance that will generate the the random data. + * @return the given RNG, wrapped in a {@link RandomGenerator}. + */ + public static RandomGenerator createRandomGenerator(final Random rng) { + return new RandomGenerator() { + /** {@inheritDoc} */ + public void setSeed(int seed) { + rng.setSeed((long) seed); + } + + /** {@inheritDoc} */ + public void setSeed(int[] seed) { + rng.setSeed(convertToLong(seed)); + } + + /** {@inheritDoc} */ + public void setSeed(long seed) { + rng.setSeed(seed); + } + + /** {@inheritDoc} */ + public void nextBytes(byte[] bytes) { + rng.nextBytes(bytes); + } + + /** {@inheritDoc} */ + public int nextInt() { + return rng.nextInt(); + } + + /** {@inheritDoc} */ + public int nextInt(int n) { + if (n <= 0) { + throw new NotStrictlyPositiveException(n); + } + return rng.nextInt(n); + } + + /** {@inheritDoc} */ + public long nextLong() { + return rng.nextLong(); + } + + /** {@inheritDoc} */ + public boolean nextBoolean() { + return rng.nextBoolean(); + } + + /** {@inheritDoc} */ + public float nextFloat() { + return rng.nextFloat(); + } + + /** {@inheritDoc} */ + public double nextDouble() { + return rng.nextDouble(); + } + + /** {@inheritDoc} */ + public double nextGaussian() { + return rng.nextGaussian(); + } + }; + } + + /** + * Converts seed from one representation to another. + * + * @param seed Original seed. + * @return the converted seed. + */ + public static long convertToLong(int[] seed) { + // The following number is the largest prime that fits + // in 32 bits (i.e. 2^32 - 5). + final long prime = 4294967291l; + + long combined = 0l; + for (int s : seed) { + combined = combined * prime + s; + } + + return combined; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/RandomVectorGenerator.java b/src/main/java/org/apache/commons/math3/random/RandomVectorGenerator.java new file mode 100644 index 0000000..55571ef --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/RandomVectorGenerator.java @@ -0,0 +1,33 @@ +/* + * 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.math3.random; + +/** + * This interface represents a random generator for whole vectors. + * + * @since 1.2 + */ +public interface RandomVectorGenerator { + + /** + * Generate a random vector. + * + * @return a random vector as an array of double. + */ + double[] nextVector(); +} diff --git a/src/main/java/org/apache/commons/math3/random/SobolSequenceGenerator.java b/src/main/java/org/apache/commons/math3/random/SobolSequenceGenerator.java new file mode 100644 index 0000000..7b91c7b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/SobolSequenceGenerator.java @@ -0,0 +1,329 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.exception.MathParseException; +import org.apache.commons.math3.exception.NotPositiveException; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.util.FastMath; + +import java.io.BufferedReader; +import java.io.IOException; +import java.io.InputStream; +import java.io.InputStreamReader; +import java.nio.charset.Charset; +import java.util.Arrays; +import java.util.NoSuchElementException; +import java.util.StringTokenizer; + +/** + * Implementation of a Sobol sequence. + * + * <p>A Sobol sequence is a low-discrepancy sequence with the property that for all values of N, its + * subsequence (x1, ... xN) has a low discrepancy. It can be used to generate pseudo-random points + * in a space S, which are equi-distributed. + * + * <p>The implementation already comes with support for up to 1000 dimensions with direction numbers + * calculated from <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Stephen Joe and Frances + * Kuo</a>. + * + * <p>The generator supports two modes: + * + * <ul> + * <li>sequential generation of points: {@link #nextVector()} + * <li>random access to the i-th point in the sequence: {@link #skipTo(int)} + * </ul> + * + * @see <a href="http://en.wikipedia.org/wiki/Sobol_sequence">Sobol sequence (Wikipedia)</a> + * @see <a href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Sobol sequence direction numbers</a> + * @since 3.3 + */ +public class SobolSequenceGenerator implements RandomVectorGenerator { + + /** The number of bits to use. */ + private static final int BITS = 52; + + /** The scaling factor. */ + private static final double SCALE = FastMath.pow(2, BITS); + + /** The maximum supported space dimension. */ + private static final int MAX_DIMENSION = 1000; + + /** The resource containing the direction numbers. */ + private static final String RESOURCE_NAME = + "/assets/org/apache/commons/math3/random/new-joe-kuo-6.1000"; + + /** Character set for file input. */ + private static final String FILE_CHARSET = "US-ASCII"; + + /** Space dimension. */ + private final int dimension; + + /** The current index in the sequence. */ + private int count = 0; + + /** The direction vector for each component. */ + private final long[][] direction; + + /** The current state. */ + private final long[] x; + + /** + * Construct a new Sobol sequence generator for the given space dimension. + * + * @param dimension the space dimension + * @throws OutOfRangeException if the space dimension is outside the allowed range of [1, 1000] + */ + public SobolSequenceGenerator(final int dimension) throws OutOfRangeException { + if (dimension < 1 || dimension > MAX_DIMENSION) { + throw new OutOfRangeException(dimension, 1, MAX_DIMENSION); + } + + // initialize the other dimensions with direction numbers from a resource + final InputStream is = getClass().getResourceAsStream(RESOURCE_NAME); + if (is == null) { + throw new MathInternalError(); + } + + this.dimension = dimension; + + // init data structures + direction = new long[dimension][BITS + 1]; + x = new long[dimension]; + + try { + initFromStream(is); + } catch (IOException e) { + // the internal resource file could not be read -> should not happen + throw new MathInternalError(); + } catch (MathParseException e) { + // the internal resource file could not be parsed -> should not happen + throw new MathInternalError(); + } finally { + try { + is.close(); + } catch (IOException e) { // NOPMD + // ignore + } + } + } + + /** + * Construct a new Sobol sequence generator for the given space dimension with direction vectors + * loaded from the given stream. + * + * <p>The expected format is identical to the files available from <a + * href="http://web.maths.unsw.edu.au/~fkuo/sobol/">Stephen Joe and Frances Kuo</a>. The first + * line will be ignored as it is assumed to contain only the column headers. The columns are: + * + * <ul> + * <li>d: the dimension + * <li>s: the degree of the primitive polynomial + * <li>a: the number representing the coefficients + * <li>m: the list of initial direction numbers + * </ul> + * + * Example: + * + * <pre> + * d s a m_i + * 2 1 0 1 + * 3 2 1 1 3 + * </pre> + * + * <p>The input stream <i>must</i> be an ASCII text containing one valid direction vector per + * line. + * + * @param dimension the space dimension + * @param is the stream to read the direction vectors from + * @throws NotStrictlyPositiveException if the space dimension is < 1 + * @throws OutOfRangeException if the space dimension is outside the range [1, max], where max + * refers to the maximum dimension found in the input stream + * @throws MathParseException if the content in the stream could not be parsed successfully + * @throws IOException if an error occurs while reading from the input stream + */ + public SobolSequenceGenerator(final int dimension, final InputStream is) + throws NotStrictlyPositiveException, MathParseException, IOException { + + if (dimension < 1) { + throw new NotStrictlyPositiveException(dimension); + } + + this.dimension = dimension; + + // init data structures + direction = new long[dimension][BITS + 1]; + x = new long[dimension]; + + // initialize the other dimensions with direction numbers from the stream + int lastDimension = initFromStream(is); + if (lastDimension < dimension) { + throw new OutOfRangeException(dimension, 1, lastDimension); + } + } + + /** + * Load the direction vector for each dimension from the given stream. + * + * <p>The input stream <i>must</i> be an ASCII text containing one valid direction vector per + * line. + * + * @param is the input stream to read the direction vector from + * @return the last dimension that has been read from the input stream + * @throws IOException if the stream could not be read + * @throws MathParseException if the content could not be parsed successfully + */ + private int initFromStream(final InputStream is) throws MathParseException, IOException { + + // special case: dimension 1 -> use unit initialization + for (int i = 1; i <= BITS; i++) { + direction[0][i] = 1l << (BITS - i); + } + + final Charset charset = Charset.forName(FILE_CHARSET); + final BufferedReader reader = new BufferedReader(new InputStreamReader(is, charset)); + int dim = -1; + + try { + // ignore first line + reader.readLine(); + + int lineNumber = 2; + int index = 1; + String line = null; + while ((line = reader.readLine()) != null) { + StringTokenizer st = new StringTokenizer(line, " "); + try { + dim = Integer.parseInt(st.nextToken()); + if (dim >= 2 && dim <= dimension) { // we have found the right dimension + final int s = Integer.parseInt(st.nextToken()); + final int a = Integer.parseInt(st.nextToken()); + final int[] m = new int[s + 1]; + for (int i = 1; i <= s; i++) { + m[i] = Integer.parseInt(st.nextToken()); + } + initDirectionVector(index++, a, m); + } + + if (dim > dimension) { + return dim; + } + } catch (NoSuchElementException e) { + throw new MathParseException(line, lineNumber); + } catch (NumberFormatException e) { + throw new MathParseException(line, lineNumber); + } + lineNumber++; + } + } finally { + reader.close(); + } + + return dim; + } + + /** + * Calculate the direction numbers from the given polynomial. + * + * @param d the dimension, zero-based + * @param a the coefficients of the primitive polynomial + * @param m the initial direction numbers + */ + private void initDirectionVector(final int d, final int a, final int[] m) { + final int s = m.length - 1; + for (int i = 1; i <= s; i++) { + direction[d][i] = ((long) m[i]) << (BITS - i); + } + for (int i = s + 1; i <= BITS; i++) { + direction[d][i] = direction[d][i - s] ^ (direction[d][i - s] >> s); + for (int k = 1; k <= s - 1; k++) { + direction[d][i] ^= ((a >> (s - 1 - k)) & 1) * direction[d][i - k]; + } + } + } + + /** {@inheritDoc} */ + public double[] nextVector() { + final double[] v = new double[dimension]; + if (count == 0) { + count++; + return v; + } + + // find the index c of the rightmost 0 + int c = 1; + int value = count - 1; + while ((value & 1) == 1) { + value >>= 1; + c++; + } + + for (int i = 0; i < dimension; i++) { + x[i] ^= direction[i][c]; + v[i] = (double) x[i] / SCALE; + } + count++; + return v; + } + + /** + * Skip to the i-th point in the Sobol sequence. + * + * <p>This operation can be performed in O(1). + * + * @param index the index in the sequence to skip to + * @return the i-th point in the Sobol sequence + * @throws NotPositiveException if index < 0 + */ + public double[] skipTo(final int index) throws NotPositiveException { + if (index == 0) { + // reset x vector + Arrays.fill(x, 0); + } else { + final int i = index - 1; + final long grayCode = i ^ (i >> 1); // compute the gray code of i = i XOR floor(i / 2) + for (int j = 0; j < dimension; j++) { + long result = 0; + for (int k = 1; k <= BITS; k++) { + final long shift = grayCode >> (k - 1); + if (shift == 0) { + // stop, as all remaining bits will be zero + break; + } + // the k-th bit of i + final long ik = shift & 1; + result ^= ik * direction[j][k]; + } + x[j] = result; + } + } + count = index; + return nextVector(); + } + + /** + * Returns the index i of the next point in the Sobol sequence that will be returned by calling + * {@link #nextVector()}. + * + * @return the index of the next point + */ + public int getNextIndex() { + return count; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/StableRandomGenerator.java b/src/main/java/org/apache/commons/math3/random/StableRandomGenerator.java new file mode 100644 index 0000000..b43a770 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/StableRandomGenerator.java @@ -0,0 +1,143 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.FastMath; + +/** + * This class provides a stable normalized random generator. It samples from a stable distribution + * with location parameter 0 and scale 1. + * + * <p>The implementation uses the Chambers-Mallows-Stuck method as described in <i>Handbook of + * computational statistics: concepts and methods</i> by James E. Gentle, Wolfgang Härdle, + * Yuichi Mori. + * + * @since 3.0 + */ +public class StableRandomGenerator implements NormalizedRandomGenerator { + /** Underlying generator. */ + private final RandomGenerator generator; + + /** stability parameter */ + private final double alpha; + + /** skewness parameter */ + private final double beta; + + /** cache of expression value used in generation */ + private final double zeta; + + /** + * Create a new generator. + * + * @param generator underlying random generator to use + * @param alpha Stability parameter. Must be in range (0, 2] + * @param beta Skewness parameter. Must be in range [-1, 1] + * @throws NullArgumentException if generator is null + * @throws OutOfRangeException if {@code alpha <= 0} or {@code alpha > 2} or {@code beta < -1} + * or {@code beta > 1} + */ + public StableRandomGenerator( + final RandomGenerator generator, final double alpha, final double beta) + throws NullArgumentException, OutOfRangeException { + if (generator == null) { + throw new NullArgumentException(); + } + + if (!(alpha > 0d && alpha <= 2d)) { + throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_LEFT, alpha, 0, 2); + } + + if (!(beta >= -1d && beta <= 1d)) { + throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_SIMPLE, beta, -1, 1); + } + + this.generator = generator; + this.alpha = alpha; + this.beta = beta; + if (alpha < 2d && beta != 0d) { + zeta = beta * FastMath.tan(FastMath.PI * alpha / 2); + } else { + zeta = 0d; + } + } + + /** + * Generate a random scalar with zero location and unit scale. + * + * @return a random scalar with zero location and unit scale + */ + public double nextNormalizedDouble() { + // we need 2 uniform random numbers to calculate omega and phi + double omega = -FastMath.log(generator.nextDouble()); + double phi = FastMath.PI * (generator.nextDouble() - 0.5); + + // Normal distribution case (Box-Muller algorithm) + if (alpha == 2d) { + return FastMath.sqrt(2d * omega) * FastMath.sin(phi); + } + + double x; + // when beta = 0, zeta is zero as well + // Thus we can exclude it from the formula + if (beta == 0d) { + // Cauchy distribution case + if (alpha == 1d) { + x = FastMath.tan(phi); + } else { + x = + FastMath.pow(omega * FastMath.cos((1 - alpha) * phi), 1d / alpha - 1d) + * FastMath.sin(alpha * phi) + / FastMath.pow(FastMath.cos(phi), 1d / alpha); + } + } else { + // Generic stable distribution + double cosPhi = FastMath.cos(phi); + // to avoid rounding errors around alpha = 1 + if (FastMath.abs(alpha - 1d) > 1e-8) { + double alphaPhi = alpha * phi; + double invAlphaPhi = phi - alphaPhi; + x = + (FastMath.sin(alphaPhi) + zeta * FastMath.cos(alphaPhi)) + / cosPhi + * (FastMath.cos(invAlphaPhi) + zeta * FastMath.sin(invAlphaPhi)) + / FastMath.pow(omega * cosPhi, (1 - alpha) / alpha); + } else { + double betaPhi = FastMath.PI / 2 + beta * phi; + x = + 2d + / FastMath.PI + * (betaPhi * FastMath.tan(phi) + - beta + * FastMath.log( + FastMath.PI + / 2d + * omega + * cosPhi + / betaPhi)); + + if (alpha != 1d) { + x += beta * FastMath.tan(FastMath.PI * alpha / 2); + } + } + } + return x; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/SynchronizedRandomGenerator.java b/src/main/java/org/apache/commons/math3/random/SynchronizedRandomGenerator.java new file mode 100644 index 0000000..410ec3c --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/SynchronizedRandomGenerator.java @@ -0,0 +1,95 @@ +/* + * 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.math3.random; + +/** + * Any {@link RandomGenerator} implementation can be thread-safe if it is used through an instance + * of this class. This is achieved by enclosing calls to the methods of the actual generator inside + * the overridden {@code synchronized} methods of this class. + * + * @since 3.1 + */ +public class SynchronizedRandomGenerator implements RandomGenerator { + /** Object to which all calls will be delegated. */ + private final RandomGenerator wrapped; + + /** + * Creates a synchronized wrapper for the given {@code RandomGenerator} instance. + * + * @param rng Generator whose methods will be called through their corresponding overridden + * synchronized version. To ensure thread-safety, the wrapped generator <em>must</em> not be + * used directly. + */ + public SynchronizedRandomGenerator(RandomGenerator rng) { + wrapped = rng; + } + + /** {@inheritDoc} */ + public synchronized void setSeed(int seed) { + wrapped.setSeed(seed); + } + + /** {@inheritDoc} */ + public synchronized void setSeed(int[] seed) { + wrapped.setSeed(seed); + } + + /** {@inheritDoc} */ + public synchronized void setSeed(long seed) { + wrapped.setSeed(seed); + } + + /** {@inheritDoc} */ + public synchronized void nextBytes(byte[] bytes) { + wrapped.nextBytes(bytes); + } + + /** {@inheritDoc} */ + public synchronized int nextInt() { + return wrapped.nextInt(); + } + + /** {@inheritDoc} */ + public synchronized int nextInt(int n) { + return wrapped.nextInt(n); + } + + /** {@inheritDoc} */ + public synchronized long nextLong() { + return wrapped.nextLong(); + } + + /** {@inheritDoc} */ + public synchronized boolean nextBoolean() { + return wrapped.nextBoolean(); + } + + /** {@inheritDoc} */ + public synchronized float nextFloat() { + return wrapped.nextFloat(); + } + + /** {@inheritDoc} */ + public synchronized double nextDouble() { + return wrapped.nextDouble(); + } + + /** {@inheritDoc} */ + public synchronized double nextGaussian() { + return wrapped.nextGaussian(); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/UncorrelatedRandomVectorGenerator.java b/src/main/java/org/apache/commons/math3/random/UncorrelatedRandomVectorGenerator.java new file mode 100644 index 0000000..ec38c4d --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/UncorrelatedRandomVectorGenerator.java @@ -0,0 +1,91 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.DimensionMismatchException; + +import java.util.Arrays; + +/** + * A {@link RandomVectorGenerator} that generates vectors with uncorrelated components. Components + * of generated vectors follow (independent) Gaussian distributions, with parameters supplied in the + * constructor. + * + * @since 1.2 + */ +public class UncorrelatedRandomVectorGenerator implements RandomVectorGenerator { + + /** Underlying scalar generator. */ + private final NormalizedRandomGenerator generator; + + /** Mean vector. */ + private final double[] mean; + + /** Standard deviation vector. */ + private final double[] standardDeviation; + + /** + * Simple constructor. + * + * <p>Build an uncorrelated random vector generator from its mean and standard deviation + * vectors. + * + * @param mean expected mean values for each component + * @param standardDeviation standard deviation for each component + * @param generator underlying generator for uncorrelated normalized components + */ + public UncorrelatedRandomVectorGenerator( + double[] mean, double[] standardDeviation, NormalizedRandomGenerator generator) { + if (mean.length != standardDeviation.length) { + throw new DimensionMismatchException(mean.length, standardDeviation.length); + } + this.mean = mean.clone(); + this.standardDeviation = standardDeviation.clone(); + this.generator = generator; + } + + /** + * Simple constructor. + * + * <p>Build a null mean random and unit standard deviation uncorrelated vector generator + * + * @param dimension dimension of the vectors to generate + * @param generator underlying generator for uncorrelated normalized components + */ + public UncorrelatedRandomVectorGenerator(int dimension, NormalizedRandomGenerator generator) { + mean = new double[dimension]; + standardDeviation = new double[dimension]; + Arrays.fill(standardDeviation, 1.0); + this.generator = generator; + } + + /** + * Generate an uncorrelated random vector. + * + * @return a random vector as a newly built array of double + */ + public double[] nextVector() { + + double[] random = new double[mean.length]; + for (int i = 0; i < random.length; ++i) { + random[i] = mean[i] + standardDeviation[i] * generator.nextNormalizedDouble(); + } + + return random; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/UniformRandomGenerator.java b/src/main/java/org/apache/commons/math3/random/UniformRandomGenerator.java new file mode 100644 index 0000000..4db7f8e --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/UniformRandomGenerator.java @@ -0,0 +1,58 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.util.FastMath; + +/** + * This class implements a normalized uniform random generator. + * + * <p>Since it is a normalized random generator, it generates values from a uniform distribution + * with mean equal to 0 and standard deviation equal to 1. Generated values fall in the range + * [-√3, +√3]. + * + * @since 1.2 + */ +public class UniformRandomGenerator implements NormalizedRandomGenerator { + + /** Square root of three. */ + private static final double SQRT3 = FastMath.sqrt(3.0); + + /** Underlying generator. */ + private final RandomGenerator generator; + + /** + * Create a new generator. + * + * @param generator underlying random generator to use + */ + public UniformRandomGenerator(RandomGenerator generator) { + this.generator = generator; + } + + /** + * Generate a random scalar with null mean and unit standard deviation. + * + * <p>The number generated is uniformly distributed between -&sqrt;(3) and +&sqrt;(3). + * + * @return a random scalar with null mean and unit standard deviation + */ + public double nextNormalizedDouble() { + return SQRT3 * (2 * generator.nextDouble() - 1.0); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/UnitSphereRandomVectorGenerator.java b/src/main/java/org/apache/commons/math3/random/UnitSphereRandomVectorGenerator.java new file mode 100644 index 0000000..ddcaa97 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/UnitSphereRandomVectorGenerator.java @@ -0,0 +1,74 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.util.FastMath; + +/** + * Generate random vectors isotropically located on the surface of a sphere. + * + * @since 2.1 + */ +public class UnitSphereRandomVectorGenerator implements RandomVectorGenerator { + /** RNG used for generating the individual components of the vectors. */ + private final RandomGenerator rand; + + /** Space dimension. */ + private final int dimension; + + /** + * @param dimension Space dimension. + * @param rand RNG for the individual components of the vectors. + */ + public UnitSphereRandomVectorGenerator(final int dimension, final RandomGenerator rand) { + this.dimension = dimension; + this.rand = rand; + } + + /** + * Create an object that will use a default RNG ({@link MersenneTwister}), in order to generate + * the individual components. + * + * @param dimension Space dimension. + */ + public UnitSphereRandomVectorGenerator(final int dimension) { + this(dimension, new MersenneTwister()); + } + + /** {@inheritDoc} */ + public double[] nextVector() { + final double[] v = new double[dimension]; + + // See http://mathworld.wolfram.com/SpherePointPicking.html for example. + // Pick a point by choosing a standard Gaussian for each element, and then + // normalizing to unit length. + double normSq = 0; + for (int i = 0; i < dimension; i++) { + final double comp = rand.nextGaussian(); + v[i] = comp; + normSq += comp * comp; + } + + final double f = 1 / FastMath.sqrt(normSq); + for (int i = 0; i < dimension; i++) { + v[i] *= f; + } + + return v; + } +} diff --git a/src/main/java/org/apache/commons/math3/random/ValueServer.java b/src/main/java/org/apache/commons/math3/random/ValueServer.java new file mode 100644 index 0000000..9c15292 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/ValueServer.java @@ -0,0 +1,465 @@ +/* + * 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.math3.random; + +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.ZeroException; +import org.apache.commons.math3.exception.util.LocalizedFormats; + +import java.io.BufferedReader; +import java.io.IOException; +import java.io.InputStreamReader; +import java.net.MalformedURLException; +import java.net.URL; + +/** + * Generates values for use in simulation applications. + * + * <p>How values are generated is determined by the <code>mode</code> property. + * + * <p>Supported <code>mode</code> values are: + * + * <ul> + * <li>DIGEST_MODE -- uses an empirical distribution + * <li>REPLAY_MODE -- replays data from <code>valuesFileURL</code> + * <li>UNIFORM_MODE -- generates uniformly distributed random values with mean = <code>mu</code> + * <li>EXPONENTIAL_MODE -- generates exponentially distributed random values with mean = <code>mu + * </code> + * <li>GAUSSIAN_MODE -- generates Gaussian distributed random values with mean = <code>mu</code> + * and standard deviation = <code>sigma</code> + * <li>CONSTANT_MODE -- returns <code>mu</code> every time. + * </ul> + */ +public class ValueServer { + + /** Use empirical distribution. */ + public static final int DIGEST_MODE = 0; + + /** Replay data from valuesFilePath. */ + public static final int REPLAY_MODE = 1; + + /** Uniform random deviates with mean = μ. */ + public static final int UNIFORM_MODE = 2; + + /** Exponential random deviates with mean = μ. */ + public static final int EXPONENTIAL_MODE = 3; + + /** Gaussian random deviates with mean = μ, std dev = σ. */ + public static final int GAUSSIAN_MODE = 4; + + /** Always return mu */ + public static final int CONSTANT_MODE = 5; + + /** mode determines how values are generated. */ + private int mode = 5; + + /** URI to raw data values. */ + private URL valuesFileURL = null; + + /** Mean for use with non-data-driven modes. */ + private double mu = 0.0; + + /** Standard deviation for use with GAUSSIAN_MODE. */ + private double sigma = 0.0; + + /** Empirical probability distribution for use with DIGEST_MODE. */ + private EmpiricalDistribution empiricalDistribution = null; + + /** File pointer for REPLAY_MODE. */ + private BufferedReader filePointer = null; + + /** RandomDataImpl to use for random data generation. */ + private final RandomDataGenerator randomData; + + // Data generation modes ====================================== + + /** Creates new ValueServer */ + public ValueServer() { + randomData = new RandomDataGenerator(); + } + + /** + * Construct a ValueServer instance using a RandomDataImpl as its source of random data. + * + * @param randomData the RandomDataImpl instance used to source random data + * @since 3.0 + * @deprecated use {@link #ValueServer(RandomGenerator)} + */ + @Deprecated + public ValueServer(RandomDataImpl randomData) { + this.randomData = randomData.getDelegate(); + } + + /** + * Construct a ValueServer instance using a RandomGenerator as its source of random data. + * + * @since 3.1 + * @param generator source of random data + */ + public ValueServer(RandomGenerator generator) { + this.randomData = new RandomDataGenerator(generator); + } + + /** + * Returns the next generated value, generated according to the mode value (see MODE constants). + * + * @return generated value + * @throws IOException in REPLAY_MODE if a file I/O error occurs + * @throws MathIllegalStateException if mode is not recognized + * @throws MathIllegalArgumentException if the underlying random generator thwrows one + */ + public double getNext() + throws IOException, MathIllegalStateException, MathIllegalArgumentException { + switch (mode) { + case DIGEST_MODE: + return getNextDigest(); + case REPLAY_MODE: + return getNextReplay(); + case UNIFORM_MODE: + return getNextUniform(); + case EXPONENTIAL_MODE: + return getNextExponential(); + case GAUSSIAN_MODE: + return getNextGaussian(); + case CONSTANT_MODE: + return mu; + default: + throw new MathIllegalStateException( + LocalizedFormats.UNKNOWN_MODE, + mode, + "DIGEST_MODE", + DIGEST_MODE, + "REPLAY_MODE", + REPLAY_MODE, + "UNIFORM_MODE", + UNIFORM_MODE, + "EXPONENTIAL_MODE", + EXPONENTIAL_MODE, + "GAUSSIAN_MODE", + GAUSSIAN_MODE, + "CONSTANT_MODE", + CONSTANT_MODE); + } + } + + /** + * Fills the input array with values generated using getNext() repeatedly. + * + * @param values array to be filled + * @throws IOException in REPLAY_MODE if a file I/O error occurs + * @throws MathIllegalStateException if mode is not recognized + * @throws MathIllegalArgumentException if the underlying random generator thwrows one + */ + public void fill(double[] values) + throws IOException, MathIllegalStateException, MathIllegalArgumentException { + for (int i = 0; i < values.length; i++) { + values[i] = getNext(); + } + } + + /** + * Returns an array of length <code>length</code> with values generated using getNext() + * repeatedly. + * + * @param length length of output array + * @return array of generated values + * @throws IOException in REPLAY_MODE if a file I/O error occurs + * @throws MathIllegalStateException if mode is not recognized + * @throws MathIllegalArgumentException if the underlying random generator thwrows one + */ + public double[] fill(int length) + throws IOException, MathIllegalStateException, MathIllegalArgumentException { + double[] out = new double[length]; + for (int i = 0; i < length; i++) { + out[i] = getNext(); + } + return out; + } + + /** + * Computes the empirical distribution using values from the file in <code>valuesFileURL</code>, + * using the default number of bins. + * + * <p><code>valuesFileURL</code> must exist and be readable by *this at runtime. + * + * <p>This method must be called before using <code>getNext()</code> with <code> + * mode = DIGEST_MODE</code> + * + * @throws IOException if an I/O error occurs reading the input file + * @throws NullArgumentException if the {@code valuesFileURL} has not been set + * @throws ZeroException if URL contains no data + */ + public void computeDistribution() throws IOException, ZeroException, NullArgumentException { + computeDistribution(EmpiricalDistribution.DEFAULT_BIN_COUNT); + } + + /** + * Computes the empirical distribution using values from the file in <code>valuesFileURL</code> + * and <code>binCount</code> bins. + * + * <p><code>valuesFileURL</code> must exist and be readable by this process at runtime. + * + * <p>This method must be called before using <code>getNext()</code> with <code> + * mode = DIGEST_MODE</code> + * + * @param binCount the number of bins used in computing the empirical distribution + * @throws NullArgumentException if the {@code valuesFileURL} has not been set + * @throws IOException if an error occurs reading the input file + * @throws ZeroException if URL contains no data + */ + public void computeDistribution(int binCount) + throws NullArgumentException, IOException, ZeroException { + empiricalDistribution = + new EmpiricalDistribution(binCount, randomData.getRandomGenerator()); + empiricalDistribution.load(valuesFileURL); + mu = empiricalDistribution.getSampleStats().getMean(); + sigma = empiricalDistribution.getSampleStats().getStandardDeviation(); + } + + /** + * Returns the data generation mode. See {@link ValueServer the class javadoc} for description + * of the valid values of this property. + * + * @return Value of property mode. + */ + public int getMode() { + return mode; + } + + /** + * Sets the data generation mode. + * + * @param mode New value of the data generation mode. + */ + public void setMode(int mode) { + this.mode = mode; + } + + /** + * Returns the URL for the file used to build the empirical distribution when using {@link + * #DIGEST_MODE}. + * + * @return Values file URL. + */ + public URL getValuesFileURL() { + return valuesFileURL; + } + + /** + * Sets the {@link #getValuesFileURL() values file URL} using a string URL representation. + * + * @param url String representation for new valuesFileURL. + * @throws MalformedURLException if url is not well formed + */ + public void setValuesFileURL(String url) throws MalformedURLException { + this.valuesFileURL = new URL(url); + } + + /** + * Sets the the {@link #getValuesFileURL() values file URL}. + * + * <p>The values file <i>must</i> be an ASCII text file containing one valid numeric entry per + * line. + * + * @param url URL of the values file. + */ + public void setValuesFileURL(URL url) { + this.valuesFileURL = url; + } + + /** + * Returns the {@link EmpiricalDistribution} used when operating in {@value #DIGEST_MODE}. + * + * @return EmpircalDistribution built by {@link #computeDistribution()} + */ + public EmpiricalDistribution getEmpiricalDistribution() { + return empiricalDistribution; + } + + /** + * Resets REPLAY_MODE file pointer to the beginning of the <code>valuesFileURL</code>. + * + * @throws IOException if an error occurs opening the file + * @throws NullPointerException if the {@code valuesFileURL} has not been set. + */ + public void resetReplayFile() throws IOException { + if (filePointer != null) { + try { + filePointer.close(); + filePointer = null; + } catch (IOException ex) { // NOPMD + // ignore + } + } + filePointer = + new BufferedReader(new InputStreamReader(valuesFileURL.openStream(), "UTF-8")); + } + + /** + * Closes {@code valuesFileURL} after use in REPLAY_MODE. + * + * @throws IOException if an error occurs closing the file + */ + public void closeReplayFile() throws IOException { + if (filePointer != null) { + filePointer.close(); + filePointer = null; + } + } + + /** + * Returns the mean used when operating in {@link #GAUSSIAN_MODE}, {@link #EXPONENTIAL_MODE} or + * {@link #UNIFORM_MODE}. When operating in {@link #CONSTANT_MODE}, this is the constant value + * always returned. Calling {@link #computeDistribution()} sets this value to the overall mean + * of the values in the {@link #getValuesFileURL() values file}. + * + * @return Mean used in data generation. + */ + public double getMu() { + return mu; + } + + /** + * Sets the {@link #getMu() mean} used in data generation. Note that calling this method after + * {@link #computeDistribution()} has been called will have no effect on data generated in + * {@link #DIGEST_MODE}. + * + * @param mu new Mean value. + */ + public void setMu(double mu) { + this.mu = mu; + } + + /** + * Returns the standard deviation used when operating in {@link #GAUSSIAN_MODE}. Calling {@link + * #computeDistribution()} sets this value to the overall standard deviation of the values in + * the {@link #getValuesFileURL() values file}. This property has no effect when the data + * generation mode is not {@link #GAUSSIAN_MODE}. + * + * @return Standard deviation used when operating in {@link #GAUSSIAN_MODE}. + */ + public double getSigma() { + return sigma; + } + + /** + * Sets the {@link #getSigma() standard deviation} used in {@link #GAUSSIAN_MODE}. + * + * @param sigma New standard deviation. + */ + public void setSigma(double sigma) { + this.sigma = sigma; + } + + /** + * Reseeds the random data generator. + * + * @param seed Value with which to reseed the {@link RandomDataImpl} used to generate random + * data. + */ + public void reSeed(long seed) { + randomData.reSeed(seed); + } + + // ------------- private methods --------------------------------- + + /** + * Gets a random value in DIGEST_MODE. + * + * <p><strong>Preconditions</strong>: + * + * <ul> + * <li>Before this method is called, <code>computeDistribution()</code> must have completed + * successfully; otherwise an <code>IllegalStateException</code> will be thrown + * </ul> + * + * @return next random value from the empirical distribution digest + * @throws MathIllegalStateException if digest has not been initialized + */ + private double getNextDigest() throws MathIllegalStateException { + if ((empiricalDistribution == null) || (empiricalDistribution.getBinStats().size() == 0)) { + throw new MathIllegalStateException(LocalizedFormats.DIGEST_NOT_INITIALIZED); + } + return empiricalDistribution.getNextValue(); + } + + /** + * Gets next sequential value from the <code>valuesFileURL</code>. + * + * <p>Throws an IOException if the read fails. + * + * <p>This method will open the <code>valuesFileURL</code> if there is no replay file open. + * + * <p>The <code>valuesFileURL</code> will be closed and reopened to wrap around from EOF to BOF + * if EOF is encountered. EOFException (which is a kind of IOException) may still be thrown if + * the <code>valuesFileURL</code> is empty. + * + * @return next value from the replay file + * @throws IOException if there is a problem reading from the file + * @throws MathIllegalStateException if URL contains no data + * @throws NumberFormatException if an invalid numeric string is encountered in the file + */ + private double getNextReplay() throws IOException, MathIllegalStateException { + String str = null; + if (filePointer == null) { + resetReplayFile(); + } + if ((str = filePointer.readLine()) == null) { + // we have probably reached end of file, wrap around from EOF to BOF + closeReplayFile(); + resetReplayFile(); + if ((str = filePointer.readLine()) == null) { + throw new MathIllegalStateException( + LocalizedFormats.URL_CONTAINS_NO_DATA, valuesFileURL); + } + } + return Double.parseDouble(str); + } + + /** + * Gets a uniformly distributed random value with mean = mu. + * + * @return random uniform value + * @throws MathIllegalArgumentException if the underlying random generator thwrows one + */ + private double getNextUniform() throws MathIllegalArgumentException { + return randomData.nextUniform(0, 2 * mu); + } + + /** + * Gets an exponentially distributed random value with mean = mu. + * + * @return random exponential value + * @throws MathIllegalArgumentException if the underlying random generator thwrows one + */ + private double getNextExponential() throws MathIllegalArgumentException { + return randomData.nextExponential(mu); + } + + /** + * Gets a Gaussian distributed random value with mean = mu and standard deviation = sigma. + * + * @return random Gaussian value + * @throws MathIllegalArgumentException if the underlying random generator thwrows one + */ + private double getNextGaussian() throws MathIllegalArgumentException { + return randomData.nextGaussian(mu, sigma); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/Well1024a.java b/src/main/java/org/apache/commons/math3/random/Well1024a.java new file mode 100644 index 0000000..72c7f16 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/Well1024a.java @@ -0,0 +1,110 @@ +/* + * 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.math3.random; + +/** + * This class implements the WELL1024a pseudo-random number generator from François Panneton, + * Pierre L'Ecuyer and Makoto Matsumoto. + * + * <p>This generator is described in a paper by François Panneton, Pierre L'Ecuyer and Makoto + * Matsumoto <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM Transactions on Mathematical + * Software, 32, 1 (2006). The errata for the paper are in <a + * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>. + * + * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number + * generator</a> + * @since 2.2 + */ +public class Well1024a extends AbstractWell { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 5680173464174485492L; + + /** Number of bits in the pool. */ + private static final int K = 1024; + + /** First parameter of the algorithm. */ + private static final int M1 = 3; + + /** Second parameter of the algorithm. */ + private static final int M2 = 24; + + /** Third parameter of the algorithm. */ + private static final int M3 = 10; + + /** + * Creates a new random number generator. + * + * <p>The instance is initialized using the current time as the seed. + */ + public Well1024a() { + super(K, M1, M2, M3); + } + + /** + * Creates a new random number generator using a single int seed. + * + * @param seed the initial seed (32 bits integer) + */ + public Well1024a(int seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using an int array seed. + * + * @param seed the initial seed (32 bits integers array), if null the seed of the generator will + * be related to the current time + */ + public Well1024a(int[] seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using a single long seed. + * + * @param seed the initial seed (64 bits integer) + */ + public Well1024a(long seed) { + super(K, M1, M2, M3, seed); + } + + /** {@inheritDoc} */ + @Override + protected int next(final int bits) { + + final int indexRm1 = iRm1[index]; + + final int v0 = v[index]; + final int vM1 = v[i1[index]]; + final int vM2 = v[i2[index]]; + final int vM3 = v[i3[index]]; + + final int z0 = v[indexRm1]; + final int z1 = v0 ^ (vM1 ^ (vM1 >>> 8)); + final int z2 = (vM2 ^ (vM2 << 19)) ^ (vM3 ^ (vM3 << 14)); + final int z3 = z1 ^ z2; + final int z4 = (z0 ^ (z0 << 11)) ^ (z1 ^ (z1 << 7)) ^ (z2 ^ (z2 << 13)); + + v[index] = z3; + v[indexRm1] = z4; + index = indexRm1; + + return z4 >>> (32 - bits); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/Well19937a.java b/src/main/java/org/apache/commons/math3/random/Well19937a.java new file mode 100644 index 0000000..a93358b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/Well19937a.java @@ -0,0 +1,112 @@ +/* + * 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.math3.random; + +/** + * This class implements the WELL19937a pseudo-random number generator from François + * Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + * + * <p>This generator is described in a paper by François Panneton, Pierre L'Ecuyer and Makoto + * Matsumoto <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM Transactions on Mathematical + * Software, 32, 1 (2006). The errata for the paper are in <a + * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>. + * + * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number + * generator</a> + * @since 2.2 + */ +public class Well19937a extends AbstractWell { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -7462102162223815419L; + + /** Number of bits in the pool. */ + private static final int K = 19937; + + /** First parameter of the algorithm. */ + private static final int M1 = 70; + + /** Second parameter of the algorithm. */ + private static final int M2 = 179; + + /** Third parameter of the algorithm. */ + private static final int M3 = 449; + + /** + * Creates a new random number generator. + * + * <p>The instance is initialized using the current time as the seed. + */ + public Well19937a() { + super(K, M1, M2, M3); + } + + /** + * Creates a new random number generator using a single int seed. + * + * @param seed the initial seed (32 bits integer) + */ + public Well19937a(int seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using an int array seed. + * + * @param seed the initial seed (32 bits integers array), if null the seed of the generator will + * be related to the current time + */ + public Well19937a(int[] seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using a single long seed. + * + * @param seed the initial seed (64 bits integer) + */ + public Well19937a(long seed) { + super(K, M1, M2, M3, seed); + } + + /** {@inheritDoc} */ + @Override + protected int next(final int bits) { + + final int indexRm1 = iRm1[index]; + final int indexRm2 = iRm2[index]; + + final int v0 = v[index]; + final int vM1 = v[i1[index]]; + final int vM2 = v[i2[index]]; + final int vM3 = v[i3[index]]; + + final int z0 = (0x80000000 & v[indexRm1]) ^ (0x7FFFFFFF & v[indexRm2]); + final int z1 = (v0 ^ (v0 << 25)) ^ (vM1 ^ (vM1 >>> 27)); + final int z2 = (vM2 >>> 9) ^ (vM3 ^ (vM3 >>> 1)); + final int z3 = z1 ^ z2; + final int z4 = z0 ^ (z1 ^ (z1 << 9)) ^ (z2 ^ (z2 << 21)) ^ (z3 ^ (z3 >>> 21)); + + v[index] = z3; + v[indexRm1] = z4; + v[indexRm2] &= 0x80000000; + index = indexRm1; + + return z4 >>> (32 - bits); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/Well19937c.java b/src/main/java/org/apache/commons/math3/random/Well19937c.java new file mode 100644 index 0000000..a699195 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/Well19937c.java @@ -0,0 +1,117 @@ +/* + * 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.math3.random; + +/** + * This class implements the WELL19937c pseudo-random number generator from François + * Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + * + * <p>This generator is described in a paper by François Panneton, Pierre L'Ecuyer and Makoto + * Matsumoto <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM Transactions on Mathematical + * Software, 32, 1 (2006). The errata for the paper are in <a + * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>. + * + * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number + * generator</a> + * @since 2.2 + */ +public class Well19937c extends AbstractWell { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -7203498180754925124L; + + /** Number of bits in the pool. */ + private static final int K = 19937; + + /** First parameter of the algorithm. */ + private static final int M1 = 70; + + /** Second parameter of the algorithm. */ + private static final int M2 = 179; + + /** Third parameter of the algorithm. */ + private static final int M3 = 449; + + /** + * Creates a new random number generator. + * + * <p>The instance is initialized using the current time as the seed. + */ + public Well19937c() { + super(K, M1, M2, M3); + } + + /** + * Creates a new random number generator using a single int seed. + * + * @param seed the initial seed (32 bits integer) + */ + public Well19937c(int seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using an int array seed. + * + * @param seed the initial seed (32 bits integers array), if null the seed of the generator will + * be related to the current time + */ + public Well19937c(int[] seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using a single long seed. + * + * @param seed the initial seed (64 bits integer) + */ + public Well19937c(long seed) { + super(K, M1, M2, M3, seed); + } + + /** {@inheritDoc} */ + @Override + protected int next(final int bits) { + + final int indexRm1 = iRm1[index]; + final int indexRm2 = iRm2[index]; + + final int v0 = v[index]; + final int vM1 = v[i1[index]]; + final int vM2 = v[i2[index]]; + final int vM3 = v[i3[index]]; + + final int z0 = (0x80000000 & v[indexRm1]) ^ (0x7FFFFFFF & v[indexRm2]); + final int z1 = (v0 ^ (v0 << 25)) ^ (vM1 ^ (vM1 >>> 27)); + final int z2 = (vM2 >>> 9) ^ (vM3 ^ (vM3 >>> 1)); + final int z3 = z1 ^ z2; + int z4 = z0 ^ (z1 ^ (z1 << 9)) ^ (z2 ^ (z2 << 21)) ^ (z3 ^ (z3 >>> 21)); + + v[index] = z3; + v[indexRm1] = z4; + v[indexRm2] &= 0x80000000; + index = indexRm1; + + // add Matsumoto-Kurita tempering + // to get a maximally-equidistributed generator + z4 ^= (z4 << 7) & 0xe46e1700; + z4 ^= (z4 << 15) & 0x9b868000; + + return z4 >>> (32 - bits); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/Well44497a.java b/src/main/java/org/apache/commons/math3/random/Well44497a.java new file mode 100644 index 0000000..7a4f34c --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/Well44497a.java @@ -0,0 +1,115 @@ +/* + * 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.math3.random; + +/** + * This class implements the WELL44497a pseudo-random number generator from François + * Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + * + * <p>This generator is described in a paper by François Panneton, Pierre L'Ecuyer and Makoto + * Matsumoto <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM Transactions on Mathematical + * Software, 32, 1 (2006). The errata for the paper are in <a + * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>. + * + * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number + * generator</a> + * @since 2.2 + */ +public class Well44497a extends AbstractWell { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -3859207588353972099L; + + /** Number of bits in the pool. */ + private static final int K = 44497; + + /** First parameter of the algorithm. */ + private static final int M1 = 23; + + /** Second parameter of the algorithm. */ + private static final int M2 = 481; + + /** Third parameter of the algorithm. */ + private static final int M3 = 229; + + /** + * Creates a new random number generator. + * + * <p>The instance is initialized using the current time as the seed. + */ + public Well44497a() { + super(K, M1, M2, M3); + } + + /** + * Creates a new random number generator using a single int seed. + * + * @param seed the initial seed (32 bits integer) + */ + public Well44497a(int seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using an int array seed. + * + * @param seed the initial seed (32 bits integers array), if null the seed of the generator will + * be related to the current time + */ + public Well44497a(int[] seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using a single long seed. + * + * @param seed the initial seed (64 bits integer) + */ + public Well44497a(long seed) { + super(K, M1, M2, M3, seed); + } + + /** {@inheritDoc} */ + @Override + protected int next(final int bits) { + + final int indexRm1 = iRm1[index]; + final int indexRm2 = iRm2[index]; + + final int v0 = v[index]; + final int vM1 = v[i1[index]]; + final int vM2 = v[i2[index]]; + final int vM3 = v[i3[index]]; + + // the values below include the errata of the original article + final int z0 = (0xFFFF8000 & v[indexRm1]) ^ (0x00007FFF & v[indexRm2]); + final int z1 = (v0 ^ (v0 << 24)) ^ (vM1 ^ (vM1 >>> 30)); + final int z2 = (vM2 ^ (vM2 << 10)) ^ (vM3 << 26); + final int z3 = z1 ^ z2; + final int z2Prime = ((z2 << 9) ^ (z2 >>> 23)) & 0xfbffffff; + final int z2Second = ((z2 & 0x00020000) != 0) ? (z2Prime ^ 0xb729fcec) : z2Prime; + final int z4 = z0 ^ (z1 ^ (z1 >>> 20)) ^ z2Second ^ z3; + + v[index] = z3; + v[indexRm1] = z4; + v[indexRm2] &= 0xFFFF8000; + index = indexRm1; + + return z4 >>> (32 - bits); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/Well44497b.java b/src/main/java/org/apache/commons/math3/random/Well44497b.java new file mode 100644 index 0000000..9da51a5 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/Well44497b.java @@ -0,0 +1,122 @@ +/* + * 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.math3.random; + +/** + * This class implements the WELL44497b pseudo-random number generator from François + * Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + * + * <p>This generator is described in a paper by François Panneton, Pierre L'Ecuyer and Makoto + * Matsumoto <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM Transactions on Mathematical + * Software, 32, 1 (2006). The errata for the paper are in <a + * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>. + * + * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number + * generator</a> + * @since 2.2 + */ +public class Well44497b extends AbstractWell { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 4032007538246675492L; + + /** Number of bits in the pool. */ + private static final int K = 44497; + + /** First parameter of the algorithm. */ + private static final int M1 = 23; + + /** Second parameter of the algorithm. */ + private static final int M2 = 481; + + /** Third parameter of the algorithm. */ + private static final int M3 = 229; + + /** + * Creates a new random number generator. + * + * <p>The instance is initialized using the current time as the seed. + */ + public Well44497b() { + super(K, M1, M2, M3); + } + + /** + * Creates a new random number generator using a single int seed. + * + * @param seed the initial seed (32 bits integer) + */ + public Well44497b(int seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using an int array seed. + * + * @param seed the initial seed (32 bits integers array), if null the seed of the generator will + * be related to the current time + */ + public Well44497b(int[] seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using a single long seed. + * + * @param seed the initial seed (64 bits integer) + */ + public Well44497b(long seed) { + super(K, M1, M2, M3, seed); + } + + /** {@inheritDoc} */ + @Override + protected int next(final int bits) { + + // compute raw value given by WELL44497a generator + // which is NOT maximally-equidistributed + final int indexRm1 = iRm1[index]; + final int indexRm2 = iRm2[index]; + + final int v0 = v[index]; + final int vM1 = v[i1[index]]; + final int vM2 = v[i2[index]]; + final int vM3 = v[i3[index]]; + + // the values below include the errata of the original article + final int z0 = (0xFFFF8000 & v[indexRm1]) ^ (0x00007FFF & v[indexRm2]); + final int z1 = (v0 ^ (v0 << 24)) ^ (vM1 ^ (vM1 >>> 30)); + final int z2 = (vM2 ^ (vM2 << 10)) ^ (vM3 << 26); + final int z3 = z1 ^ z2; + final int z2Prime = ((z2 << 9) ^ (z2 >>> 23)) & 0xfbffffff; + final int z2Second = ((z2 & 0x00020000) != 0) ? (z2Prime ^ 0xb729fcec) : z2Prime; + int z4 = z0 ^ (z1 ^ (z1 >>> 20)) ^ z2Second ^ z3; + + v[index] = z3; + v[indexRm1] = z4; + v[indexRm2] &= 0xFFFF8000; + index = indexRm1; + + // add Matsumoto-Kurita tempering + // to get a maximally-equidistributed generator + z4 ^= (z4 << 7) & 0x93dd1400; + z4 ^= (z4 << 15) & 0xfa118000; + + return z4 >>> (32 - bits); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/Well512a.java b/src/main/java/org/apache/commons/math3/random/Well512a.java new file mode 100644 index 0000000..d09cbb1 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/Well512a.java @@ -0,0 +1,111 @@ +/* + * 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.math3.random; + +/** + * This class implements the WELL512a pseudo-random number generator from François Panneton, + * Pierre L'Ecuyer and Makoto Matsumoto. + * + * <p>This generator is described in a paper by François Panneton, Pierre L'Ecuyer and Makoto + * Matsumoto <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM Transactions on Mathematical + * Software, 32, 1 (2006). The errata for the paper are in <a + * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>. + * + * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number + * generator</a> + * @since 2.2 + */ +public class Well512a extends AbstractWell { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -6104179812103820574L; + + /** Number of bits in the pool. */ + private static final int K = 512; + + /** First parameter of the algorithm. */ + private static final int M1 = 13; + + /** Second parameter of the algorithm. */ + private static final int M2 = 9; + + /** Third parameter of the algorithm. */ + private static final int M3 = 5; + + /** + * Creates a new random number generator. + * + * <p>The instance is initialized using the current time as the seed. + */ + public Well512a() { + super(K, M1, M2, M3); + } + + /** + * Creates a new random number generator using a single int seed. + * + * @param seed the initial seed (32 bits integer) + */ + public Well512a(int seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using an int array seed. + * + * @param seed the initial seed (32 bits integers array), if null the seed of the generator will + * be related to the current time + */ + public Well512a(int[] seed) { + super(K, M1, M2, M3, seed); + } + + /** + * Creates a new random number generator using a single long seed. + * + * @param seed the initial seed (64 bits integer) + */ + public Well512a(long seed) { + super(K, M1, M2, M3, seed); + } + + /** {@inheritDoc} */ + @Override + protected int next(final int bits) { + + final int indexRm1 = iRm1[index]; + + final int vi = v[index]; + final int vi1 = v[i1[index]]; + final int vi2 = v[i2[index]]; + final int z0 = v[indexRm1]; + + // the values below include the errata of the original article + final int z1 = (vi ^ (vi << 16)) ^ (vi1 ^ (vi1 << 15)); + final int z2 = vi2 ^ (vi2 >>> 11); + final int z3 = z1 ^ z2; + final int z4 = + (z0 ^ (z0 << 2)) ^ (z1 ^ (z1 << 18)) ^ (z2 << 28) ^ (z3 ^ ((z3 << 5) & 0xda442d24)); + + v[index] = z3; + v[indexRm1] = z4; + index = indexRm1; + + return z4 >>> (32 - bits); + } +} diff --git a/src/main/java/org/apache/commons/math3/random/package-info.java b/src/main/java/org/apache/commons/math3/random/package-info.java new file mode 100644 index 0000000..212b018 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/random/package-info.java @@ -0,0 +1,115 @@ +/* + * 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. + */ +/** + * Random number and random data generators. + * + * <p>Commons-math provides a few pseudo random number generators. The top level interface is + * RandomGenerator. It is implemented by three classes: + * + * <ul> + * <li>{@link org.apache.commons.math3.random.JDKRandomGenerator JDKRandomGenerator} that extends + * the JDK provided generator + * <li>AbstractRandomGenerator as a helper for users generators + * <li>BitStreamGenerator which is an abstract class for several generators and which in turn is + * extended by: + * <ul> + * <li>{@link org.apache.commons.math3.random.MersenneTwister MersenneTwister} + * <li>{@link org.apache.commons.math3.random.Well512a Well512a} + * <li>{@link org.apache.commons.math3.random.Well1024a Well1024a} + * <li>{@link org.apache.commons.math3.random.Well19937a Well19937a} + * <li>{@link org.apache.commons.math3.random.Well19937c Well19937c} + * <li>{@link org.apache.commons.math3.random.Well44497a Well44497a} + * <li>{@link org.apache.commons.math3.random.Well44497b Well44497b} + * </ul> + * </ul> + * + * <p>The JDK provided generator is a simple one that can be used only for very simple needs. The + * Mersenne Twister is a fast generator with very good properties well suited for Monte-Carlo + * simulation. It is equidistributed for generating vectors up to dimension 623 and has a huge + * period: 2<sup>19937</sup> - 1 (which is a Mersenne prime). This generator is described in a paper + * by Makoto Matsumoto and Takuji Nishimura in 1998: <a + * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">Mersenne Twister: A + * 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator</a>, ACM Transactions on + * Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3--30. The WELL generators are + * a family of generators with period ranging from 2<sup>512</sup> - 1 to 2<sup>44497</sup> - 1 + * (this last one is also a Mersenne prime) with even better properties than Mersenne Twister. These + * generators are described in a paper by François Panneton, Pierre L'Ecuyer and Makoto + * Matsumoto <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM Transactions on Mathematical + * Software, 32, 1 (2006). The errata for the paper are in <a + * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>. + * + * <p>For simple sampling, any of these generators is sufficient. For Monte-Carlo simulations the + * JDK generator does not have any of the good mathematical properties of the other generators, so + * it should be avoided. The Mersenne twister and WELL generators have equidistribution properties + * proven according to their bits pool size which is directly linked to their period (all of them + * have maximal period, i.e. a generator with size n pool has a period 2<sup>n</sup>-1). They also + * have equidistribution properties for 32 bits blocks up to s/32 dimension where s is their pool + * size. So WELL19937c for exemple is equidistributed up to dimension 623 (19937/32). This means a + * Monte-Carlo simulation generating a vector of n variables at each iteration has some guarantees + * on the properties of the vector as long as its dimension does not exceed the limit. However, + * since we use bits from two successive 32 bits generated integers to create one double, this limit + * is smaller when the variables are of type double. so for Monte-Carlo simulation where less the 16 + * doubles are generated at each round, WELL1024 may be sufficient. If a larger number of doubles + * are needed a generator with a larger pool would be useful. + * + * <p>The WELL generators are more modern then MersenneTwister (the paper describing than has been + * published in 2006 instead of 1998) and fix some of its (few) drawbacks. If initialization array + * contains many zero bits, MersenneTwister may take a very long time (several hundreds of thousands + * of iterations to reach a steady state with a balanced number of zero and one in its bits pool). + * So the WELL generators are better to <i>escape zeroland</i> as explained by the WELL generators + * creators. The Well19937a and Well44497a generator are not maximally equidistributed (i.e. there + * are some dimensions or bits blocks size for which they are not equidistributed). The Well512a, + * Well1024a, Well19937c and Well44497b are maximally equidistributed for blocks size up to 32 bits + * (they should behave correctly also for double based on more than 32 bits blocks, but + * equidistribution is not proven at these blocks sizes). + * + * <p>The MersenneTwister generator uses a 624 elements integer array, so it consumes less than 2.5 + * kilobytes. The WELL generators use 6 integer arrays with a size equal to the pool size, so for + * example the WELL44497b generator uses about 33 kilobytes. This may be important if a very large + * number of generator instances were used at the same time. + * + * <p>All generators are quite fast. As an example, here are some comparisons, obtained on a 64 bits + * JVM on a linux computer with a 2008 processor (AMD phenom Quad 9550 at 2.2 GHz). The generation + * rate for MersenneTwister was about 27 millions doubles per second (remember we generate two 32 + * bits integers for each double). Generation rates for other PRNG, relative to MersenneTwister: + * + * <p> + * + * <table border="1" align="center"> + * <tr BGCOLOR="#CCCCFF"><td colspan="2"><font size="+2">Example of performances</font></td></tr> + * <tr BGCOLOR="#EEEEFF"><font size="+1"><td>Name</td><td>generation rate (relative to MersenneTwister)</td></font></tr> + * <tr><td>{@link org.apache.commons.math3.random.MersenneTwister MersenneTwister}</td><td>1</td></tr> + * <tr><td>{@link org.apache.commons.math3.random.JDKRandomGenerator JDKRandomGenerator}</td><td>between 0.96 and 1.16</td></tr> + * <tr><td>{@link org.apache.commons.math3.random.Well512a Well512a}</td><td>between 0.85 and 0.88</td></tr> + * <tr><td>{@link org.apache.commons.math3.random.Well1024a Well1024a}</td><td>between 0.63 and 0.73</td></tr> + * <tr><td>{@link org.apache.commons.math3.random.Well19937a Well19937a}</td><td>between 0.70 and 0.71</td></tr> + * <tr><td>{@link org.apache.commons.math3.random.Well19937c Well19937c}</td><td>between 0.57 and 0.71</td></tr> + * <tr><td>{@link org.apache.commons.math3.random.Well44497a Well44497a}</td><td>between 0.69 and 0.71</td></tr> + * <tr><td>{@link org.apache.commons.math3.random.Well44497b Well44497b}</td><td>between 0.65 and 0.71</td></tr> + * </table> + * + * <p>So for most simulation problems, the better generators like {@link + * org.apache.commons.math3.random.Well19937c Well19937c} and {@link + * org.apache.commons.math3.random.Well44497b Well44497b} are probably very good choices. + * + * <p>Note that <em>none</em> of these generators are suitable for cryptography. They are devoted to + * simulation, and to generate very long series with strong properties on the series as a whole + * (equidistribution, no correlation ...). They do not attempt to create small series but with very + * strong properties of unpredictability as needed in cryptography. + */ +package org.apache.commons.math3.random; |