From dee0849a9704d532af0b550146cbafbaa6ee1d19 Mon Sep 17 00:00:00 2001 From: Raymond Date: Thu, 2 Apr 2015 10:43:13 -0700 Subject: third party library: apache-commons-math Change-Id: I52a325624a7f0dd652b362a9840626d6d9f3c42b --- .../math/random/AbstractRandomGenerator.java | 272 ++++++ .../apache/commons/math/random/AbstractWell.java | 186 ++++ .../commons/math/random/BitsStreamGenerator.java | 153 ++++ .../random/CorrelatedRandomVectorGenerator.java | 304 +++++++ .../commons/math/random/EmpiricalDistribution.java | 133 +++ .../math/random/EmpiricalDistributionImpl.java | 479 ++++++++++ .../math/random/GaussianRandomGenerator.java | 47 + .../commons/math/random/JDKRandomGenerator.java | 50 ++ .../commons/math/random/MersenneTwister.java | 259 ++++++ .../math/random/NormalizedRandomGenerator.java | 38 + .../apache/commons/math/random/RandomAdaptor.java | 198 +++++ .../org/apache/commons/math/random/RandomData.java | 272 ++++++ .../apache/commons/math/random/RandomDataImpl.java | 966 +++++++++++++++++++++ .../commons/math/random/RandomGenerator.java | 148 ++++ .../commons/math/random/RandomVectorGenerator.java | 35 + .../random/UncorrelatedRandomVectorGenerator.java | 93 ++ .../math/random/UniformRandomGenerator.java | 61 ++ .../random/UnitSphereRandomVectorGenerator.java | 84 ++ .../apache/commons/math/random/ValueServer.java | 385 ++++++++ .../org/apache/commons/math/random/Well1024a.java | 106 +++ .../org/apache/commons/math/random/Well19937a.java | 108 +++ .../org/apache/commons/math/random/Well19937c.java | 115 +++ .../org/apache/commons/math/random/Well44497a.java | 111 +++ .../org/apache/commons/math/random/Well44497b.java | 119 +++ .../org/apache/commons/math/random/Well512a.java | 107 +++ .../org/apache/commons/math/random/package.html | 132 +++ 26 files changed, 4961 insertions(+) create mode 100644 src/main/java/org/apache/commons/math/random/AbstractRandomGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/AbstractWell.java create mode 100644 src/main/java/org/apache/commons/math/random/BitsStreamGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/EmpiricalDistribution.java create mode 100644 src/main/java/org/apache/commons/math/random/EmpiricalDistributionImpl.java create mode 100644 src/main/java/org/apache/commons/math/random/GaussianRandomGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/JDKRandomGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/MersenneTwister.java create mode 100644 src/main/java/org/apache/commons/math/random/NormalizedRandomGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/RandomAdaptor.java create mode 100644 src/main/java/org/apache/commons/math/random/RandomData.java create mode 100644 src/main/java/org/apache/commons/math/random/RandomDataImpl.java create mode 100644 src/main/java/org/apache/commons/math/random/RandomGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/RandomVectorGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/UncorrelatedRandomVectorGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/UniformRandomGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/UnitSphereRandomVectorGenerator.java create mode 100644 src/main/java/org/apache/commons/math/random/ValueServer.java create mode 100644 src/main/java/org/apache/commons/math/random/Well1024a.java create mode 100644 src/main/java/org/apache/commons/math/random/Well19937a.java create mode 100644 src/main/java/org/apache/commons/math/random/Well19937c.java create mode 100644 src/main/java/org/apache/commons/math/random/Well44497a.java create mode 100644 src/main/java/org/apache/commons/math/random/Well44497b.java create mode 100644 src/main/java/org/apache/commons/math/random/Well512a.java create mode 100644 src/main/java/org/apache/commons/math/random/package.html (limited to 'src/main/java/org/apache/commons/math/random') diff --git a/src/main/java/org/apache/commons/math/random/AbstractRandomGenerator.java b/src/main/java/org/apache/commons/math/random/AbstractRandomGenerator.java new file mode 100644 index 0000000..afec63e --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/AbstractRandomGenerator.java @@ -0,0 +1,272 @@ +/* + * 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.math.random; + +import org.apache.commons.math.exception.NotStrictlyPositiveException; +import org.apache.commons.math.util.FastMath; + +/** + * Abstract class implementing the {@link RandomGenerator} interface. + * Default implementations for all methods other than {@link #nextDouble()} and + * {@link #setSeed(long)} are provided. + *

+ * All data generation methods are based on {@code code nextDouble()}. + * Concrete implementations must override + * this method and should provide better / more + * performant implementations of the other methods if the underlying PRNG + * supplies them.

+ * + * @since 1.1 + * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ + */ +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}. Implemementations 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 underyling random number generator using a + * {@code long} seed. Sequences of values generated starting with the + * same seeds should be identical. + *

+ * 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. + *

+ * 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 = 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 232 possible {@code int} values + * should be produced with (approximately) equal probability. + *

+ * The default implementation provided here returns + * 

+     * (int) (nextDouble() * Integer.MAX_VALUE)
+     *

+ * + * @return the next pseudorandom, uniformly distributed {@code int} + * value from this random number generator's sequence + */ + public int nextInt() { + return (int) (nextDouble() * 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. + *

+ * The default implementation returns + * 

+     * (int) (nextDouble() * n
+     *

+ * + * @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 + * 264 possible {@code long} values + * should be produced with (approximately) equal probability. + *

+ * The default implementation returns + * 

+     * (long) (nextDouble() * Long.MAX_VALUE)
+     *

+ * + * @return the next pseudorandom, uniformly distributed {@code long} + *value from this random number generator's sequence + */ + public long nextLong() { + return (long) (nextDouble() * Long.MAX_VALUE); + } + + /** + * Returns the next pseudorandom, uniformly distributed + * {@code boolean} value from this random number generator's + * sequence. + *

+ * The default implementation returns + * 

+     * nextDouble() <= 0.5
+     *

+ * + * @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. + *

+ * The default implementation returns + * 

+     * (float) nextDouble() 
+     *

+ * + * @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. + *

+ * 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. + *

+ * The default implementation uses the Polar Method + * due to G.E.P. Box, M.E. Muller and G. Marsaglia, as described in + * D. Knuth, The Art of Computer Programming, 3.4.1C.

+ *

+ * 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/math/random/AbstractWell.java b/src/main/java/org/apache/commons/math/random/AbstractWell.java new file mode 100644 index 0000000..96e18a2 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/AbstractWell.java @@ -0,0 +1,186 @@ +/* + * 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.math.random; + +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. + + *

This generator is described in a paper by François Panneton, + * Pierre L'Ecuyer and Makoto Matsumoto Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2 ACM + * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper + * are in wellrng-errata.txt.

+ + * @see WELL Random number generator + * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $ + * @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. + *

The instance is initialized using the current time 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, System.currentTimeMillis()); + } + + /** 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. + *

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. + *

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 related to the current time + */ + @Override + public void setSeed(final int[] seed) { + + if (seed == null) { + setSeed(System.currentTimeMillis()); + return; + } + + System.arraycopy(seed, 0, v, 0, Math.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; + + } + + /** Reinitialize the generator as if just built with the given long seed. + *

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/math/random/BitsStreamGenerator.java b/src/main/java/org/apache/commons/math/random/BitsStreamGenerator.java new file mode 100644 index 0000000..8979473 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/BitsStreamGenerator.java @@ -0,0 +1,153 @@ +/* + * 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.math.random; + +import org.apache.commons.math.exception.NotStrictlyPositiveException; +import org.apache.commons.math.util.FastMath; + +/** Base class for random number generators that generates bits streams. + + * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ + * @since 2.0 + + */ +public abstract class BitsStreamGenerator implements RandomGenerator { + + /** 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. + *

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 void nextBytes(byte[] bytes) { + int i = 0; + final int iEnd = bytes.length - 3; + while (i < iEnd) { + final int random = next(32); + bytes[i] = (byte) (random & 0xff); + bytes[i + 1] = (byte) ((random >> 8) & 0xff); + bytes[i + 2] = (byte) ((random >> 16) & 0xff); + bytes[i + 3] = (byte) ((random >> 24) & 0xff); + i += 4; + } + int random = next(32); + while (i < bytes.length) { + bytes[i++] = (byte) (random & 0xff); + random = random >> 8; + } + } + + /** {@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} */ + public int nextInt(int n) throws IllegalArgumentException { + + if (n < 1) { + throw new NotStrictlyPositiveException(n); + } + + // find bit mask for n + int mask = n; + mask |= mask >> 1; + mask |= mask >> 2; + mask |= mask >> 4; + mask |= mask >> 8; + mask |= mask >> 16; + + while (true) { + final int random = next(32) & mask; + if (random < n) { + return random; + } + } + + } + + /** {@inheritDoc} */ + public long nextLong() { + final long high = ((long) next(32)) << 32; + final long low = ((long) next(32)) & 0xffffffffL; + return high | low; + } + +} diff --git a/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java b/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java new file mode 100644 index 0000000..f1df51d --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java @@ -0,0 +1,304 @@ +/* + * 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.math.random; + +import org.apache.commons.math.DimensionMismatchException; +import org.apache.commons.math.linear.MatrixUtils; +import org.apache.commons.math.linear.NotPositiveDefiniteMatrixException; +import org.apache.commons.math.linear.RealMatrix; +import org.apache.commons.math.util.FastMath; + +/** + * A {@link RandomVectorGenerator} that generates vectors with with + * correlated components. + *

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.

+ *

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 + * Multivariate Normal Distribution. 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}.

+ *

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 C = UT.U + * where C is the covariance matrix and U + * is an upper-triangular matrix, we compute C = B.BT + * where B is a rectangular matrix having + * more rows than columns. The number of columns of B 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.

+ * + * @version $Revision: 1043908 $ $Date: 2010-12-09 12:53:14 +0100 (jeu. 09 déc. 2010) $ + * @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; + + /** Permutated Cholesky root of the covariance matrix. */ + private RealMatrix root; + + /** Rank of the covariance matrix. */ + private int rank; + + /** Simple constructor. + *

Build 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 + * @exception IllegalArgumentException if there is a dimension + * mismatch between the mean vector and the covariance matrix + * @exception NotPositiveDefiniteMatrixException if the + * covariance matrix is not strictly positive definite + * @exception DimensionMismatchException if the mean and covariance + * arrays dimensions don't match + */ + public CorrelatedRandomVectorGenerator(double[] mean, + RealMatrix covariance, double small, + NormalizedRandomGenerator generator) + throws NotPositiveDefiniteMatrixException, DimensionMismatchException { + + int order = covariance.getRowDimension(); + if (mean.length != order) { + throw new DimensionMismatchException(mean.length, order); + } + this.mean = mean.clone(); + + decompose(covariance, small); + + this.generator = generator; + normalized = new double[rank]; + + } + + /** Simple constructor. + *

Build 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 + * @exception NotPositiveDefiniteMatrixException if the + * covariance matrix is not strictly positive definite + */ + public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small, + NormalizedRandomGenerator generator) + throws NotPositiveDefiniteMatrixException { + + int order = covariance.getRowDimension(); + mean = new double[order]; + for (int i = 0; i < order; ++i) { + mean[i] = 0; + } + + decompose(covariance, small); + + this.generator = generator; + normalized = new double[rank]; + + } + + /** Get the underlying normalized components generator. + * @return underlying uncorrelated components generator + */ + public NormalizedRandomGenerator getGenerator() { + return generator; + } + + /** Get the root of the covariance matrix. + * The root is the rectangular matrix B such that + * the covariance matrix is equal to B.BT + * @return root of the square matrix + * @see #getRank() + */ + public RealMatrix getRootMatrix() { + return root; + } + + /** 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 rectangular + * matrix of the decomposition. + * @return rank of the square matrix. + * @see #getRootMatrix() + */ + public int getRank() { + return rank; + } + + /** Decompose the original square matrix. + *

The decomposition is based on a Choleski decomposition + * where additional transforms are performed: + *

+ * This means that rather than computing M = UT.U where U + * is an upper triangular matrix, this method computed M=B.BT + * where B is a rectangular matrix. + * @param covariance covariance matrix + * @param small diagonal elements threshold under which column are + * considered to be dependent on previous ones and are discarded + * @exception NotPositiveDefiniteMatrixException if the + * covariance matrix is not strictly positive definite + */ + private void decompose(RealMatrix covariance, double small) + throws NotPositiveDefiniteMatrixException { + + int order = covariance.getRowDimension(); + double[][] c = covariance.getData(); + double[][] b = new double[order][order]; + + int[] swap = new int[order]; + int[] index = new int[order]; + for (int i = 0; i < order; ++i) { + index[i] = i; + } + + rank = 0; + for (boolean loop = true; loop;) { + + // find maximal diagonal element + swap[rank] = rank; + for (int i = rank + 1; i < order; ++i) { + int ii = index[i]; + int isi = index[swap[i]]; + if (c[ii][ii] > c[isi][isi]) { + swap[rank] = i; + } + } + + + // swap elements + if (swap[rank] != rank) { + int tmp = index[rank]; + index[rank] = index[swap[rank]]; + index[swap[rank]] = tmp; + } + + // check diagonal element + int ir = index[rank]; + if (c[ir][ir] < small) { + + if (rank == 0) { + throw new NotPositiveDefiniteMatrixException(); + } + + // check remaining diagonal elements + for (int i = rank; i < order; ++i) { + if (c[index[i]][index[i]] < -small) { + // there is at least one sufficiently negative diagonal element, + // the covariance matrix is wrong + throw new NotPositiveDefiniteMatrixException(); + } + } + + // all remaining diagonal elements are close to zero, + // we consider we have found the rank of the covariance matrix + ++rank; + loop = false; + + } else { + + // transform the matrix + double sqrt = FastMath.sqrt(c[ir][ir]); + b[rank][rank] = sqrt; + double inverse = 1 / sqrt; + for (int i = rank + 1; i < order; ++i) { + int ii = index[i]; + double e = inverse * c[ii][ir]; + b[i][rank] = e; + c[ii][ii] -= e * e; + for (int j = rank + 1; j < i; ++j) { + int ij = index[j]; + double f = c[ii][ij] - e * b[j][rank]; + c[ii][ij] = f; + c[ij][ii] = f; + } + } + + // prepare next iteration + loop = ++rank < order; + + } + + } + + // build the root matrix + root = MatrixUtils.createRealMatrix(order, rank); + for (int i = 0; i < order; ++i) { + for (int j = 0; j < rank; ++j) { + root.setEntry(index[i], j, b[i][j]); + } + } + + } + + /** 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 < rank; ++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 < rank; ++j) { + correlated[i] += root.getEntry(i, j) * normalized[j]; + } + } + + return correlated; + + } + +} diff --git a/src/main/java/org/apache/commons/math/random/EmpiricalDistribution.java b/src/main/java/org/apache/commons/math/random/EmpiricalDistribution.java new file mode 100644 index 0000000..7f08a06 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/EmpiricalDistribution.java @@ -0,0 +1,133 @@ +/* + * 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.math.random; + +import java.io.IOException; +import java.io.File; +import java.net.URL; +import java.util.List; + +import org.apache.commons.math.stat.descriptive.StatisticalSummary; +import org.apache.commons.math.stat.descriptive.SummaryStatistics; + +/** + * Represents an + * empirical probability distribution -- a probability distribution derived + * from observed data without making any assumptions about the functional form + * of the population distribution that the data come from.

+ * Implementations of this interface maintain data structures, called + * distribution digests, that describe empirical distributions and + * support the following operations:

+ * Applications can use EmpiricalDistribution implementations 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.

+ * + * @version $Revision: 817128 $ $Date: 2009-09-21 03:30:53 +0200 (lun. 21 sept. 2009) $ + */ +public interface EmpiricalDistribution { + + /** + * Computes the empirical distribution from the provided + * array of numbers. + * + * @param dataArray the data array + */ + void load(double[] dataArray); + + /** + * Computes the empirical distribution from the input file. + * + * @param file the input file + * @throws IOException if an IO error occurs + */ + void load(File file) throws IOException; + + /** + * Computes the empirical distribution using data read from a URL. + * + * @param url url of the input file + * @throws IOException if an IO error occurs + */ + void load(URL url) throws IOException; + + /** + * Generates a random value from this distribution. + * Preconditions: + * @return the random value. + * + * @throws IllegalStateException if the distribution has not been loaded + */ + double getNextValue() throws IllegalStateException; + + + /** + * Returns a + * {@link org.apache.commons.math.stat.descriptive.StatisticalSummary} + * describing this distribution. + * Preconditions: + * + * @return the sample statistics + * @throws IllegalStateException if the distribution has not been loaded + */ + StatisticalSummary getSampleStats() throws IllegalStateException; + + /** + * Property indicating whether or not the distribution has been loaded. + * + * @return true if the distribution has been loaded + */ + boolean isLoaded(); + + /** + * Returns the number of bins. + * + * @return the number of bins + */ + int getBinCount(); + + /** + * Returns a list of + * {@link org.apache.commons.math.stat.descriptive.SummaryStatistics} + * containing statistics describing the values in each of the bins. The + * List is indexed on the bin number. + * + * @return List of bin statistics + */ + List getBinStats(); + + /** + * Returns the array of upper bounds for the bins. Bins are:
+ * [min,upperBounds[0]],(upperBounds[0],upperBounds[1]],..., + * (upperBounds[binCount-2], upperBounds[binCount-1] = max]. + * + * @return array of bin upper bounds + */ + double[] getUpperBounds(); + +} diff --git a/src/main/java/org/apache/commons/math/random/EmpiricalDistributionImpl.java b/src/main/java/org/apache/commons/math/random/EmpiricalDistributionImpl.java new file mode 100644 index 0000000..a3ae8a7 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/EmpiricalDistributionImpl.java @@ -0,0 +1,479 @@ +/* + * 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.math.random; + +import java.io.BufferedReader; +import java.io.File; +import java.io.FileReader; +import java.io.IOException; +import java.io.InputStreamReader; +import java.io.Serializable; +import java.net.URL; +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math.MathRuntimeException; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.stat.descriptive.StatisticalSummary; +import org.apache.commons.math.stat.descriptive.SummaryStatistics; +import org.apache.commons.math.util.FastMath; + +/** + * Implements EmpiricalDistribution interface. This implementation + * uses what amounts to the + * + * Variable Kernel Method with Gaussian smoothing:

+ * Digesting the input file + *

  1. Pass the file once to compute min and max.
    2. + *
  2. Divide the range from min-max into binCount "bins."
    3. + *
  3. Pass the data file again, computing bin counts and univariate + * statistics (mean, std dev.) for each of the bins
    4. + *
  4. Divide the interval (0,1) into subintervals associated with the bins, + * with the length of a bin's subinterval proportional to its count.
+ * Generating random values from the distribution
    + *
  1. Generate a uniformly distributed value in (0,1)
    2. + *
  2. Select the subinterval to which the value belongs. + *
  3. Generate a random Gaussian value with mean = mean of the associated + * bin and std dev = std dev of associated bin.

+ *USAGE NOTES:

+ * + * @version $Revision: 1003886 $ $Date: 2010-10-02 23:04:44 +0200 (sam. 02 oct. 2010) $ + */ +public class EmpiricalDistributionImpl implements Serializable, EmpiricalDistribution { + + /** Serializable version identifier */ + private static final long serialVersionUID = 5729073523949762654L; + + /** List of SummaryStatistics objects characterizing the bins */ + private final List 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; + + /** RandomData instance to use in repeated calls to getNext() */ + private final RandomData randomData = new RandomDataImpl(); + + /** + * Creates a new EmpiricalDistribution with the default bin count. + */ + public EmpiricalDistributionImpl() { + binCount = 1000; + binStats = new ArrayList(); + } + + /** + * Creates a new EmpiricalDistribution with the specified bin count. + * + * @param binCount number of bins + */ + public EmpiricalDistributionImpl(int binCount) { + this.binCount = binCount; + binStats = new ArrayList(); + } + + /** + * Computes the empirical distribution from the provided + * array of numbers. + * + * @param in the input data array + */ + public void load(double[] in) { + DataAdapter da = new ArrayDataAdapter(in); + try { + da.computeStats(); + fillBinStats(in); + } catch (IOException e) { + throw new MathRuntimeException(e); + } + loaded = true; + + } + + /** + * Computes the empirical distribution using data read from a URL. + * @param url url of the input file + * + * @throws IOException if an IO error occurs + */ + public void load(URL url) throws IOException { + BufferedReader in = + new BufferedReader(new InputStreamReader(url.openStream())); + try { + DataAdapter da = new StreamDataAdapter(in); + da.computeStats(); + if (sampleStats.getN() == 0) { + throw MathRuntimeException.createEOFException(LocalizedFormats.URL_CONTAINS_NO_DATA, + url); + } + in = new BufferedReader(new InputStreamReader(url.openStream())); + fillBinStats(in); + loaded = true; + } finally { + try { + in.close(); + } catch (IOException ex) { + // ignore + } + } + } + + /** + * Computes the empirical distribution from the input file. + * + * @param file the input file + * @throws IOException if an IO error occurs + */ + public void load(File file) throws IOException { + BufferedReader in = new BufferedReader(new FileReader(file)); + try { + DataAdapter da = new StreamDataAdapter(in); + da.computeStats(); + in = new BufferedReader(new FileReader(file)); + fillBinStats(in); + loaded = true; + } finally { + try { + in.close(); + } catch (IOException ex) { + // ignore + } + } + } + + /** + * Provides methods for computing sampleStats and + * beanStats 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; + + } + + /** + * Factory of DataAdapter objects. For every supported source + * of data (array of doubles, file, etc.) an instance of the proper object + * is returned. + */ + private class DataAdapterFactory{ + /** + * Creates a DataAdapter from a data object + * + * @param in object providing access to the data + * @return DataAdapter instance + */ + public DataAdapter getAdapter(Object in) { + if (in instanceof BufferedReader) { + BufferedReader inputStream = (BufferedReader) in; + return new StreamDataAdapter(inputStream); + } else if (in instanceof double[]) { + double[] inputArray = (double[]) in; + return new ArrayDataAdapter(inputArray); + } else { + throw MathRuntimeException.createIllegalArgumentException( + LocalizedFormats.INPUT_DATA_FROM_UNSUPPORTED_DATASOURCE, + in.getClass().getName(), + BufferedReader.class.getName(), double[].class.getName()); + } + } + } + /** + * DataAdapter 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 + */ + public 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.valueOf(str).doubleValue(); + sampleStats.addValue(val); + } + inputStream.close(); + inputStream = null; + } + } + + /** + * DataAdapter 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 + */ + public ArrayDataAdapter(double[] in){ + super(); + 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 in object providing access to the data + * @throws IOException if an IO error occurs + */ + private void fillBinStats(Object in) throws IOException { + // Set up grid + min = sampleStats.getMin(); + max = sampleStats.getMax(); + delta = (max - min)/(Double.valueOf(binCount)).doubleValue(); + + // 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 + DataAdapterFactory aFactory = new DataAdapterFactory(); + DataAdapter da = aFactory.getAdapter(in); + 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. + * + * @return the random value. + * @throws IllegalStateException if the distribution has not been loaded + */ + public double getNextValue() throws IllegalStateException { + + if (!loaded) { + throw MathRuntimeException.createIllegalStateException(LocalizedFormats.DISTRIBUTION_NOT_LOADED); + } + + // Start with a uniformly distributed random number in (0,1) + double x = FastMath.random(); + + // Use this to select the bin and generate a Gaussian within the bin + for (int i = 0; i < binCount; i++) { + if (x <= upperBounds[i]) { + SummaryStatistics stats = binStats.get(i); + if (stats.getN() > 0) { + if (stats.getStandardDeviation() > 0) { // more than one obs + return randomData.nextGaussian + (stats.getMean(),stats.getStandardDeviation()); + } else { + return stats.getMean(); // only one obs in bin + } + } + } + } + throw new MathRuntimeException(LocalizedFormats.NO_BIN_SELECTED); + } + + /** + * Returns a {@link StatisticalSummary} describing this distribution. + * Preconditions:
    + *
  • the distribution must be loaded before invoking this method
+ * + * @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 getBinStats() { + return binStats; + } + + /** + *

Returns a fresh copy of the array of upper bounds for the bins. + * Bins are:
+ * [min,upperBounds[0]],(upperBounds[0],upperBounds[1]],..., + * (upperBounds[binCount-2], upperBounds[binCount-1] = max].

+ * + *

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]; + binUpperBounds[0] = min + delta; + for (int i = 1; i < binCount - 1; i++) { + binUpperBounds[i] = binUpperBounds[i-1] + delta; + } + 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.

+ * + *

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 + */ + 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; + } +} diff --git a/src/main/java/org/apache/commons/math/random/GaussianRandomGenerator.java b/src/main/java/org/apache/commons/math/random/GaussianRandomGenerator.java new file mode 100644 index 0000000..6bd9775 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/GaussianRandomGenerator.java @@ -0,0 +1,47 @@ +/* + * 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.math.random; + +/** + * This class is a gaussian normalized random generator for scalars. + *

This class is a simple wrapper around the {@link + * RandomGenerator#nextGaussian} method.

+ * @version $Revision: 811827 $ $Date: 2009-09-06 17:32:50 +0200 (dim. 06 sept. 2009) $ + * @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/math/random/JDKRandomGenerator.java b/src/main/java/org/apache/commons/math/random/JDKRandomGenerator.java new file mode 100644 index 0000000..573ef61 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/JDKRandomGenerator.java @@ -0,0 +1,50 @@ +/* + * 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.math.random; + +import java.util.Random; + +/** + * Extension of java.util.Random to implement + * {@link RandomGenerator}. + * + * @since 1.1 + * @version $Revision: 811685 $ $Date: 2009-09-05 19:36:48 +0200 (sam. 05 sept. 2009) $ + */ +public class JDKRandomGenerator extends Random implements RandomGenerator { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -7745277476784028798L; + + /** {@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); + } + +} diff --git a/src/main/java/org/apache/commons/math/random/MersenneTwister.java b/src/main/java/org/apache/commons/math/random/MersenneTwister.java new file mode 100644 index 0000000..6a6fadd --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/MersenneTwister.java @@ -0,0 +1,259 @@ +/* + * 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.math.random; + +import java.io.Serializable; + +import org.apache.commons.math.util.FastMath; + + +/** This class implements a powerful pseudo-random number generator + * developed by Makoto Matsumoto and Takuji Nishimura during + * 1996-1997. + + *

This generator features an extremely long period + * (219937-1) and 623-dimensional equidistribution up to 32 + * bits accuracy. The home page for this generator is located at + * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html.

+ + *

This generator is described in a paper by Makoto Matsumoto and + * Takuji Nishimura in 1998: Mersenne + * Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random + * Number Generator, ACM Transactions on Modeling and Computer + * Simulation, Vol. 8, No. 1, January 1998, pp 3--30

+ + *

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:

+ + * + * + + * + + * + *
Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + * All rights reserved.
Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + *
    + *
  1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer.
    2. + *
  2. 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.
    3. + *
  3. The names of its contributors may not be used to endorse or promote + * products derived from this software without specific prior written + * permission.
    4. + *
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.
+ + * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ + * @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. + *

The instance is initialized using the current time as the + * seed.

+ */ + public MersenneTwister() { + mt = new int[N]; + setSeed(System.currentTimeMillis()); + } + + /** 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. + *

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; + 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; + } + } + + /** Reinitialize the generator as if just built with the given int array seed. + *

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 related to the current time + */ + @Override + public void setSeed(int[] seed) { + + if (seed == null) { + setSeed(System.currentTimeMillis()); + 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 + + } + + /** Reinitialize the generator as if just built with the given long seed. + *

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. + *

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/math/random/NormalizedRandomGenerator.java b/src/main/java/org/apache/commons/math/random/NormalizedRandomGenerator.java new file mode 100644 index 0000000..0f13b81 --- /dev/null +++ b/src/main/java/org/apache/commons/math/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.math.random; + +/** + * This interface represent a normalized random generator for + * scalars. + * Normalized generator provide null mean and unit standard deviation scalars. + * @version $Revision: 811786 $ $Date: 2009-09-06 11:36:08 +0200 (dim. 06 sept. 2009) $ + * @since 1.2 + */ +public interface NormalizedRandomGenerator { + + /** Generate a random scalar with null mean and unit standard deviation. + *

This method does not 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/math/random/RandomAdaptor.java b/src/main/java/org/apache/commons/math/random/RandomAdaptor.java new file mode 100644 index 0000000..8091a48 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/RandomAdaptor.java @@ -0,0 +1,198 @@ +/* + * 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.math.random; + +import java.util.Random; + +/** + * Extension of java.util.Random wrapping a + * {@link RandomGenerator}. + * + * @since 1.1 + * @version $Revision: 1003886 $ $Date: 2010-10-02 23:04:44 +0200 (sam. 02 oct. 2010) $ + */ +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 Random using the supplied + * RandomGenerator. + * + * @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 + * boolean value from this random number generator's + * sequence. + * + * @return the next pseudorandom, uniformly distributed + * boolean 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 + * double value between 0.0 and + * 1.0 from this random number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed + * double value between 0.0 and + * 1.0 from this random number generator's sequence + */ + @Override + public double nextDouble() { + return randomGenerator.nextDouble(); + } + + /** + * Returns the next pseudorandom, uniformly distributed float + * value between 0.0 and 1.0 from this random + * number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed float + * value between 0.0 and 1.0 from this + * random number generator's sequence + */ + @Override + public float nextFloat() { + return randomGenerator.nextFloat(); + } + + /** + * Returns the next pseudorandom, Gaussian ("normally") distributed + * double value with mean 0.0 and standard + * deviation 1.0 from this random number generator's sequence. + * + * @return the next pseudorandom, Gaussian ("normally") distributed + * double value with mean 0.0 and + * standard deviation 1.0 from this random number + * generator's sequence + */ + @Override + public double nextGaussian() { + return randomGenerator.nextGaussian(); + } + + /** + * Returns the next pseudorandom, uniformly distributed int + * value from this random number generator's sequence. + * All 232 possible int values + * should be produced with (approximately) equal probability. + * + * @return the next pseudorandom, uniformly distributed int + * value from this random number generator's sequence + */ + @Override + public int nextInt() { + return randomGenerator.nextInt(); + } + + /** + * Returns a pseudorandom, uniformly distributed 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 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 long + * value from this random number generator's sequence. All + * 264 possible long values + * should be produced with (approximately) equal probability. + * + * @return the next pseudorandom, uniformly distributed long + *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/math/random/RandomData.java b/src/main/java/org/apache/commons/math/random/RandomData.java new file mode 100644 index 0000000..0fc5136 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/RandomData.java @@ -0,0 +1,272 @@ +/* + * 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.math.random; +import java.util.Collection; + +/** + * Random data generation utilities. + * @version $Revision: 780975 $ $Date: 2009-06-02 11:05:37 +0200 (mar. 02 juin 2009) $ + */ +public interface RandomData { + /** + * Generates a random string of hex characters of length + * len. + *

+ * The generated string will be random, but not cryptographically + * secure. To generate cryptographically secure strings, use + * nextSecureHexString

+ *

+ * Preconditions:

    + *
  • len > 0 (otherwise an IllegalArgumentException + * is thrown.)
    • + *

+ * + * @param len the length of the string to be generated + * @return random string of hex characters of length len + */ + String nextHexString(int len); + + /** + * Generates a uniformly distributed random integer between + * lower and upper (endpoints included). + *

+ * The generated integer will be random, but not cryptographically secure. + * To generate cryptographically secure integer sequences, use + * nextSecureInt.

+ *

+ * Preconditions:

    + *
  • lower < upper (otherwise an IllegalArgumentException + * is thrown.)
    • + *

+ * + * @param lower lower bound for generated integer + * @param upper upper bound for generated integer + * @return a random integer greater than or equal to lower + * and less than or equal to upper. + */ + int nextInt(int lower, int upper); + + /** + * Generates a uniformly distributed random long integer between + * lower and upper (endpoints included). + *

+ * The generated long integer values will be random, but not + * cryptographically secure. + * To generate cryptographically secure sequences of longs, use + * nextSecureLong

+ *

+ * Preconditions:

    + *
  • lower < upper (otherwise an IllegalArgumentException + * is thrown.)
    • + *

+ * + * @param lower lower bound for generated integer + * @param upper upper bound for generated integer + * @return a random integer greater than or equal to lower + * and less than or equal to upper. + */ + long nextLong(long lower, long upper); + + /** + * Generates a random string of hex characters from a secure random + * sequence. + *

+ * If cryptographic security is not required, + * use nextHexString().

+ *

+ * Preconditions:

    + *
  • len > 0 (otherwise an IllegalArgumentException + * is thrown.)
    • + *

+ * @param len length of return string + * @return the random hex string + */ + String nextSecureHexString(int len); + + /** + * Generates a uniformly distributed random integer between + * lower and upper (endpoints included) + * from a secure random sequence. + *

+ * Sequences of integers generated using this method will be + * cryptographically secure. If cryptographic security is not required, + * nextInt should be used instead of this method.

+ *

+ * Definition: + * + * Secure Random Sequence

+ *

+ * Preconditions:

    + *
  • lower < upper (otherwise an IllegalArgumentException + * is thrown.)
    • + *

+ * + * @param lower lower bound for generated integer + * @param upper upper bound for generated integer + * @return a random integer greater than or equal to lower + * and less than or equal to upper. + */ + int nextSecureInt(int lower, int upper); + + /** + * Generates a random long integer between lower + * and upper (endpoints included). + *

+ * Sequences of long values generated using this method will be + * cryptographically secure. If cryptographic security is not required, + * nextLong should be used instead of this method.

+ *

+ * Definition: + * + * Secure Random Sequence

+ *

+ * Preconditions:

    + *
  • lower < upper (otherwise an IllegalArgumentException + * is thrown.)
    • + *

+ * + * @param lower lower bound for generated integer + * @param upper upper bound for generated integer + * @return a long integer greater than or equal to lower + * and less than or equal to upper. + */ + long nextSecureLong(long lower, long upper); + + /** + * Generates a random value from the Poisson distribution with + * the given mean. + *

+ * Definition: + * + * Poisson Distribution

+ *

+ * Preconditions:

    + *
  • The specified mean must be positive (otherwise an + * IllegalArgumentException is thrown.)
    • + *

+ * @param mean Mean of the distribution + * @return poisson deviate with the specified mean + */ + long nextPoisson(double mean); + + /** + * Generates a random value from the + * Normal (or Gaussian) distribution with the given mean + * and standard deviation. + *

+ * Definition: + * + * Normal Distribution

+ *

+ * Preconditions:

    + *
  • sigma > 0 (otherwise an IllegalArgumentException + * is thrown.)
    • + *

+ * @param mu Mean of the distribution + * @param sigma Standard deviation of the distribution + * @return random value from Gaussian distribution with mean = mu, + * standard deviation = sigma + */ + double nextGaussian(double mu, double sigma); + + /** + * Generates a random value from the exponential distribution + * with expected value = mean. + *

+ * Definition: + * + * Exponential Distribution

+ *

+ * Preconditions:

    + *
  • mu >= 0 (otherwise an IllegalArgumentException + * is thrown.)
    • + *

+ * @param mean Mean of the distribution + * @return random value from exponential distribution + */ + double nextExponential(double mean); + + /** + * Generates a uniformly distributed random value from the open interval + * (lower,upper) (i.e., endpoints excluded). + *

+ * Definition: + * + * Uniform Distribution lower and + * upper - lower are the + * + * location and scale parameters, respectively.

+ *

+ * Preconditions:

    + *
  • lower < upper (otherwise an IllegalArgumentException + * is thrown.)
    • + *

+ * + * @param lower lower endpoint of the interval of support + * @param upper upper endpoint of the interval of support + * @return uniformly distributed random value between lower + * and upper (exclusive) + */ + double nextUniform(double lower, double upper); + + /** + * Generates an integer array of length k whose entries + * are selected randomly, without repetition, from the integers + * 0 through n-1 (inclusive). + *

+ * Generated arrays represent permutations + * of n taken k at a time.

+ *

+ * Preconditions:

    + *
  • k <= n
    • + *
  • n > 0
    • + *
+ * If the preconditions are not met, an IllegalArgumentException is + * thrown.

+ * + * @param n domain of the permutation + * @param k size of the permutation + * @return random k-permutation of n + */ + int[] nextPermutation(int n, int k); + + /** + * Returns an array of k objects selected randomly + * from the Collection c. + *

+ * Sampling from c + * is without replacement; but if c contains identical + * objects, the sample may include repeats. If all elements of + * c are distinct, the resulting object array represents a + * + * Simple Random Sample of size + * k from the elements of c.

+ *

+ * Preconditions:

    + *
  • k must be less than or equal to the size of c
    • + *
  • c must not be empty
    • + *
+ * If the preconditions are not met, an IllegalArgumentException is + * thrown.

+ * + * @param c collection to be sampled + * @param k size of the sample + * @return random sample of k elements from c + */ + Object[] nextSample(Collection c, int k); +} diff --git a/src/main/java/org/apache/commons/math/random/RandomDataImpl.java b/src/main/java/org/apache/commons/math/random/RandomDataImpl.java new file mode 100644 index 0000000..e9ccab7 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/RandomDataImpl.java @@ -0,0 +1,966 @@ +/* + * 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.math.random; + +import java.io.Serializable; +import java.security.MessageDigest; +import java.security.NoSuchAlgorithmException; +import java.security.NoSuchProviderException; +import java.security.SecureRandom; +import java.util.Collection; + +import org.apache.commons.math.MathException; +import org.apache.commons.math.distribution.BetaDistributionImpl; +import org.apache.commons.math.distribution.BinomialDistributionImpl; +import org.apache.commons.math.distribution.CauchyDistributionImpl; +import org.apache.commons.math.distribution.ChiSquaredDistributionImpl; +import org.apache.commons.math.distribution.ContinuousDistribution; +import org.apache.commons.math.distribution.FDistributionImpl; +import org.apache.commons.math.distribution.GammaDistributionImpl; +import org.apache.commons.math.distribution.HypergeometricDistributionImpl; +import org.apache.commons.math.distribution.IntegerDistribution; +import org.apache.commons.math.distribution.PascalDistributionImpl; +import org.apache.commons.math.distribution.TDistributionImpl; +import org.apache.commons.math.distribution.WeibullDistributionImpl; +import org.apache.commons.math.distribution.ZipfDistributionImpl; +import org.apache.commons.math.exception.MathInternalError; +import org.apache.commons.math.exception.NotStrictlyPositiveException; +import org.apache.commons.math.exception.NumberIsTooLargeException; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.util.FastMath; +import org.apache.commons.math.util.MathUtils; + +/** + * 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 nextSecureXxx methods. If no + * RandomGenerator is provided in the constructor, the default is + * to use a generator based on {@link java.util.Random}. To plug in a different + * implementation, either implement RandomGenerator directly or + * extend {@link AbstractRandomGenerator}. + *

+ * Supports reseeding the underlying pseudo-random number generator (PRNG). The + * SecurityProvider and Algorithm used by the + * SecureRandom instance can also be reset. + *

+ *

+ * For details on the default PRNGs, see {@link java.util.Random} and + * {@link java.security.SecureRandom}. + *

+ *

+ * Usage Notes: + *

    + *
  • + * Instance variables are used to maintain RandomGenerator and + * SecureRandom instances used in data generation. Therefore, to + * generate a random sequence of values or strings, you should use just + * one RandomDataImpl instance repeatedly.
    • + *
  • + * 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.
    • + *
  • + * When a new RandomDataImpl is created, the underlying random + * number generators are not initialized. If you do not + * explicitly seed the default non-secure generator, it is seeded with the + * current time in milliseconds on first use. The same holds for the secure + * generator. If you provide a RandomGenerator to the constructor, + * however, this generator is not reseeded by the constructor nor is it reseeded + * on first use.
    • + *
  • + * The reSeed and reSeedSecure methods delegate to the + * corresponding methods on the underlying RandomGenerator and + * SecureRandom instances. Therefore, reSeed(long) + * 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 not + * reinitialize the secure random number generator (so secure sequences started + * with calls to reseedSecure(long) won't be identical).
    • + *
  • + * This implementation is not synchronized. + *
+ *

+ * + * @version $Revision: 1061496 $ $Date: 2011-01-20 21:32:16 +0100 (jeu. 20 janv. 2011) $ + */ +public class RandomDataImpl 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 SecureRandom secRand = null; + + /** + * Construct a RandomDataImpl. + */ + public RandomDataImpl() { + } + + /** + * 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 + * @since 1.1 + */ + public RandomDataImpl(RandomGenerator rand) { + super(); + this.rand = rand; + } + + /** + * {@inheritDoc} + *

+ * Algorithm Description: hex strings are generated using a + * 2-step process. + *

    + *
  1. + * len/2+1 binary bytes are generated using the underlying Random
    2. + *
  2. + * Each binary byte is translated into 2 hex digits
    3. + *
+ *

+ * + * @param len + * the desired string length. + * @return the random string. + * @throws NotStrictlyPositiveException if {@code len <= 0}. + */ + public String nextHexString(int len) { + if (len <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.LENGTH, len); + } + + // Get a random number generator + RandomGenerator ran = getRan(); + + // 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); + } + + /** + * Generate a random int value uniformly distributed between + * lower and upper, inclusive. + * + * @param lower + * the lower bound. + * @param upper + * the upper bound. + * @return the random integer. + * @throws NumberIsTooLargeException if {@code lower >= upper}. + */ + public int nextInt(int lower, int upper) { + if (lower >= upper) { + throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND, + lower, upper, false); + } + double r = getRan().nextDouble(); + return (int) ((r * upper) + ((1.0 - r) * lower) + r); + } + + /** + * Generate a random long value uniformly distributed between + * lower and upper, inclusive. + * + * @param lower + * the lower bound. + * @param upper + * the upper bound. + * @return the random integer. + * @throws NumberIsTooLargeException if {@code lower >= upper}. + */ + public long nextLong(long lower, long upper) { + if (lower >= upper) { + throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND, + lower, upper, false); + } + double r = getRan().nextDouble(); + return (long) ((r * upper) + ((1.0 - r) * lower) + r); + } + + /** + * {@inheritDoc} + *

+ * Algorithm Description: hex strings are generated in + * 40-byte segments using a 3-step process. + *

    + *
  1. + * 20 random bytes are generated using the underlying + * SecureRandom.
    2. + *
  2. + * SHA-1 hash is applied to yield a 20-byte binary digest.
    3. + *
  3. + * Each byte of the binary digest is converted to 2 hex digits.
    4. + *
+ *

+ * + * @param len + * the length of the generated string + * @return the random string + * @throws NotStrictlyPositiveException if {@code len <= 0}. + */ + public String nextSecureHexString(int len) { + if (len <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.LENGTH, len); + } + + // Get SecureRandom and setup Digest provider + SecureRandom 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); + } + + /** + * Generate a random int value uniformly distributed between + * lower and upper, inclusive. This algorithm uses + * a secure random number generator. + * + * @param lower + * the lower bound. + * @param upper + * the upper bound. + * @return the random integer. + * @throws NumberIsTooLargeException if {@code lower >= upper}. + */ + public int nextSecureInt(int lower, int upper) { + if (lower >= upper) { + throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND, + lower, upper, false); + } + SecureRandom sec = getSecRan(); + return lower + (int) (sec.nextDouble() * (upper - lower + 1)); + } + + /** + * Generate a random long value uniformly distributed between + * lower and upper, inclusive. This algorithm uses + * a secure random number generator. + * + * @param lower + * the lower bound. + * @param upper + * the upper bound. + * @return the random integer. + * @throws NumberIsTooLargeException if {@code lower >= upper}. + */ + public long nextSecureLong(long lower, long upper) { + if (lower >= upper) { + throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND, + lower, upper, false); + } + SecureRandom sec = getSecRan(); + return lower + (long) (sec.nextDouble() * (upper - lower + 1)); + } + + /** + * {@inheritDoc} + *

+ * Algorithm Description: + *

  • For small means, uses simulation of a Poisson process + * using Uniform deviates, as described + * here. + * The Poisson process (and hence value returned) is bounded by 1000 * mean.
    • + * + *
  • For large means, uses the rejection algorithm described in
    + * Devroye, Luc. (1981).The Computer Generation of Poisson Random Variables + * Computing vol. 26 pp. 197-207.

+ * + * @param mean mean of the Poisson distribution. + * @return the random Poisson value. + * @throws NotStrictlyPositiveException if {@code mean <= 0}. + */ + public long nextPoisson(double mean) { + if (mean <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean); + } + + final RandomGenerator generator = getRan(); + + final double pivot = 40.0d; + if (mean < pivot) { + double p = FastMath.exp(-mean); + long n = 0; + double r = 1.0d; + double rnd = 1.0d; + + while (n < 1000 * mean) { + rnd = generator.nextDouble(); + r = r * rnd; + if (r >= p) { + n++; + } else { + return n; + } + } + return n; + } else { + final double lambda = FastMath.floor(mean); + final double lambdaFractional = mean - lambda; + final double logLambda = FastMath.log(lambda); + final double logLambdaFactorial = MathUtils.factorialLog((int) lambda); + final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional); + final double delta = FastMath.sqrt(lambda * FastMath.log(32 * lambda / FastMath.PI + 1)); + final double halfDelta = delta / 2; + final double twolpd = 2 * lambda + delta; + final double a1 = FastMath.sqrt(FastMath.PI * twolpd) * FastMath.exp(1 / 8 * lambda); + final double a2 = (twolpd / delta) * FastMath.exp(-delta * (1 + delta) / twolpd); + final double aSum = a1 + a2 + 1; + final double p1 = a1 / aSum; + final double p2 = a2 / aSum; + final double c1 = 1 / (8 * lambda); + + double x = 0; + double y = 0; + double v = 0; + int a = 0; + double t = 0; + double qr = 0; + double qa = 0; + for (;;) { + final double u = nextUniform(0.0, 1); + if (u <= p1) { + final double n = nextGaussian(0d, 1d); + x = n * FastMath.sqrt(lambda + halfDelta) - 0.5d; + if (x > delta || x < -lambda) { + continue; + } + y = x < 0 ? FastMath.floor(x) : FastMath.ceil(x); + final double e = nextExponential(1d); + v = -e - (n * n / 2) + c1; + } else { + if (u > p1 + p2) { + y = lambda; + break; + } else { + x = delta + (twolpd / delta) * nextExponential(1d); + y = FastMath.ceil(x); + v = -nextExponential(1d) - delta * (x + 1) / twolpd; + } + } + a = x < 0 ? 1 : 0; + t = y * (y + 1) / (2 * lambda); + if (v < -t && a == 0) { + y = lambda + y; + break; + } + qr = t * ((2 * y + 1) / (6 * lambda) - 1); + qa = qr - (t * t) / (3 * (lambda + a * (y + 1))); + if (v < qa) { + y = lambda + y; + break; + } + if (v > qr) { + continue; + } + if (v < y * logLambda - MathUtils.factorialLog((int) (y + lambda)) + logLambdaFactorial) { + y = lambda + y; + break; + } + } + return y2 + (long) y; + } + } + + /** + * Generate a random value from a Normal (a.k.a. Gaussian) distribution with + * the given mean, mu and the given standard deviation, + * sigma. + * + * @param mu + * the mean of the distribution + * @param sigma + * the standard deviation of the distribution + * @return the random Normal value + * @throws NotStrictlyPositiveException if {@code sigma <= 0}. + */ + public double nextGaussian(double mu, double sigma) { + if (sigma <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.STANDARD_DEVIATION, sigma); + } + return sigma * getRan().nextGaussian() + mu; + } + + /** + * Returns a random value from an Exponential distribution with the given + * mean. + *

+ * Algorithm Description: Uses the Inversion + * Method to generate exponentially distributed random values from + * uniform deviates. + *

+ * + * @param mean the mean of the distribution + * @return the random Exponential value + * @throws NotStrictlyPositiveException if {@code mean <= 0}. + */ + public double nextExponential(double mean) { + if (mean <= 0.0) { + throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean); + } + final RandomGenerator generator = getRan(); + double unif = generator.nextDouble(); + while (unif == 0.0d) { + unif = generator.nextDouble(); + } + return -mean * FastMath.log(unif); + } + + /** + * {@inheritDoc} + *

+ * Algorithm Description: 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). + *

+ * + * @param lower + * the lower bound. + * @param upper + * the upper bound. + * @return a uniformly distributed random value from the interval (lower, + * upper) + * @throws NumberIsTooLargeException if {@code lower >= upper}. + */ + public double nextUniform(double lower, double upper) { + if (lower >= upper) { + throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND, + lower, upper, false); + } + final RandomGenerator generator = getRan(); + + // ensure nextDouble() isn't 0.0 + double u = generator.nextDouble(); + while (u <= 0.0) { + u = generator.nextDouble(); + } + + return lower + u * (upper - lower); + } + + /** + * Generates a random value from the {@link BetaDistributionImpl Beta Distribution}. + * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) 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 + * @throws MathException if an error occurs generating the random value + * @since 2.2 + */ + public double nextBeta(double alpha, double beta) throws MathException { + return nextInversionDeviate(new BetaDistributionImpl(alpha, beta)); + } + + /** + * Generates a random value from the {@link BinomialDistributionImpl Binomial Distribution}. + * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) 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 + * @throws MathException if an error occurs generating the random value + * @since 2.2 + */ + public int nextBinomial(int numberOfTrials, double probabilityOfSuccess) throws MathException { + return nextInversionDeviate(new BinomialDistributionImpl(numberOfTrials, probabilityOfSuccess)); + } + + /** + * Generates a random value from the {@link CauchyDistributionImpl Cauchy Distribution}. + * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) 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 + * @throws MathException if an error occurs generating the random value + * @since 2.2 + */ + public double nextCauchy(double median, double scale) throws MathException { + return nextInversionDeviate(new CauchyDistributionImpl(median, scale)); + } + + /** + * Generates a random value from the {@link ChiSquaredDistributionImpl ChiSquare Distribution}. + * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion} + * to generate random values. + * + * @param df the degrees of freedom of the ChiSquare distribution + * @return random value sampled from the ChiSquare(df) distribution + * @throws MathException if an error occurs generating the random value + * @since 2.2 + */ + public double nextChiSquare(double df) throws MathException { + return nextInversionDeviate(new ChiSquaredDistributionImpl(df)); + } + + /** + * Generates a random value from the {@link FDistributionImpl F Distribution}. + * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) 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 MathException if an error occurs generating the random value + * @since 2.2 + */ + public double nextF(double numeratorDf, double denominatorDf) throws MathException { + return nextInversionDeviate(new FDistributionImpl(numeratorDf, denominatorDf)); + } + + /** + * Generates a random value from the {@link GammaDistributionImpl Gamma Distribution}. + * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion} + * to generate random values. + * + * @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 MathException if an error occurs generating the random value + * @since 2.2 + */ + public double nextGamma(double shape, double scale) throws MathException { + return nextInversionDeviate(new GammaDistributionImpl(shape, scale)); + } + + /** + * Generates a random value from the {@link HypergeometricDistributionImpl 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 MathException if an error occurs generating the random value + * @since 2.2 + */ + public int nextHypergeometric(int populationSize, int numberOfSuccesses, int sampleSize) throws MathException { + return nextInversionDeviate(new HypergeometricDistributionImpl(populationSize, numberOfSuccesses, sampleSize)); + } + + /** + * Generates a random value from the {@link PascalDistributionImpl 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 + * @throws MathException if an error occurs generating the random value + * @since 2.2 + */ + public int nextPascal(int r, double p) throws MathException { + return nextInversionDeviate(new PascalDistributionImpl(r, p)); + } + + /** + * Generates a random value from the {@link TDistributionImpl T Distribution}. + * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion} + * to generate random values. + * + * @param df the degrees of freedom of the T distribution + * @return random value from the T(df) distribution + * @throws MathException if an error occurs generating the random value + * @since 2.2 + */ + public double nextT(double df) throws MathException { + return nextInversionDeviate(new TDistributionImpl(df)); + } + + /** + * Generates a random value from the {@link WeibullDistributionImpl Weibull Distribution}. + * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) 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 + * @throws MathException if an error occurs generating the random value + * @since 2.2 + */ + public double nextWeibull(double shape, double scale) throws MathException { + return nextInversionDeviate(new WeibullDistributionImpl(shape, scale)); + } + + /** + * Generates a random value from the {@link ZipfDistributionImpl 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 + * @throws MathException if an error occurs generating the random value + * @since 2.2 + */ + public int nextZipf(int numberOfElements, double exponent) throws MathException { + return nextInversionDeviate(new ZipfDistributionImpl(numberOfElements, exponent)); + } + + /** + * Returns the RandomGenerator used to generate non-secure random data. + *

+ * Creates and initializes a default generator if null. + *

+ * + * @return the Random used to generate random data + * @since 1.1 + */ + private RandomGenerator getRan() { + if (rand == null) { + rand = new JDKRandomGenerator(); + rand.setSeed(System.currentTimeMillis()); + } + return rand; + } + + /** + * Returns the SecureRandom used to generate secure random data. + *

+ * Creates and initializes if null. + *

+ * + * @return the SecureRandom used to generate secure random data + */ + private SecureRandom getSecRan() { + if (secRand == null) { + secRand = new SecureRandom(); + secRand.setSeed(System.currentTimeMillis()); + } + return secRand; + } + + /** + * Reseeds the random number generator with the supplied seed. + *

+ * Will create and initialize if null. + *

+ * + * @param seed + * the seed value to use + */ + public void reSeed(long seed) { + if (rand == null) { + rand = new JDKRandomGenerator(); + } + rand.setSeed(seed); + } + + /** + * Reseeds the secure random number generator with the current time in + * milliseconds. + *

+ * Will create and initialize if null. + *

+ */ + public void reSeedSecure() { + if (secRand == null) { + secRand = new SecureRandom(); + } + secRand.setSeed(System.currentTimeMillis()); + } + + /** + * Reseeds the secure random number generator with the supplied seed. + *

+ * Will create and initialize if null. + *

+ * + * @param seed + * the seed value to use + */ + public void reSeedSecure(long seed) { + if (secRand == null) { + secRand = new SecureRandom(); + } + secRand.setSeed(seed); + } + + /** + * Reseeds the random number generator with the current time in + * milliseconds. + */ + public void reSeed() { + if (rand == null) { + rand = new JDKRandomGenerator(); + } + rand.setSeed(System.currentTimeMillis()); + } + + /** + * Sets the PRNG algorithm for the underlying SecureRandom instance using + * the Security Provider API. The Security Provider API is defined in + * Java Cryptography Architecture API Specification & Reference. + *

+ * USAGE NOTE: This method carries significant + * 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 = SecureRandom.getInstance(algorithm, provider); + } + + /** + * Generates an integer array of length k whose entries are + * selected randomly, without repetition, from the integers + * 0 through n-1 (inclusive). + *

+ * Generated arrays represent permutations of n taken + * k at a time. + *

+ *

+ * Preconditions: + *

    + *
  • k <= n
    • + *
  • n > 0
    • + *
+ * If the preconditions are not met, an IllegalArgumentException is thrown. + *

+ *

+ * Uses a 2-cycle permutation shuffle. The shuffling process is described + * here. + *

+ * + * @param n + * domain of the permutation (must be positive) + * @param k + * size of the permutation (must satisfy 0 < k <= n). + * @return the random permutation as an int array + * @throws NumberIsTooLargeException if {@code k > n}. + * @throws NotStrictlyPositiveException if {@code k <= 0}. + */ + public int[] nextPermutation(int n, int k) { + 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 = getNatural(n); + shuffle(index, n - k); + int[] result = new int[k]; + for (int i = 0; i < k; i++) { + result[i] = index[n - i - 1]; + } + + return result; + } + + /** + * Uses a 2-cycle permutation shuffle to generate a random permutation. + * Algorithm Description: Uses a 2-cycle permutation + * shuffle to generate a random permutation of c.size() 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, + * here + * + * @param c + * Collection to sample from. + * @param k + * sample size. + * @return the random sample. + * @throws NumberIsTooLargeException if {@code k > c.size()}. + * @throws NotStrictlyPositiveException if {@code k <= 0}. + */ + public Object[] nextSample(Collection c, int k) { + 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; + } + + /** + * Generate a random deviate from the given distribution using the + * inversion method. + * + * @param distribution Continuous distribution to generate a random value from + * @return a random value sampled from the given distribution + * @throws MathException if an error occurs computing the inverse cumulative distribution function + * @since 2.2 + */ + public double nextInversionDeviate(ContinuousDistribution distribution) throws MathException { + return distribution.inverseCumulativeProbability(nextUniform(0, 1)); + + } + + /** + * Generate a random deviate from the given distribution using the + * inversion method. + * + * @param distribution Integer distribution to generate a random value from + * @return a random value sampled from the given distribution + * @throws MathException if an error occurs computing the inverse cumulative distribution function + * @since 2.2 + */ + public int nextInversionDeviate(IntegerDistribution distribution) throws MathException { + final double target = nextUniform(0, 1); + final int glb = distribution.inverseCumulativeProbability(target); + if (distribution.cumulativeProbability(glb) == 1.0d) { // No mass above + return glb; + } else { + return glb + 1; + } + } + + // ------------------------Private methods---------------------------------- + + /** + * Uses a 2-cycle permutation shuffle to randomly re-order the last elements + * of list. + * + * @param list + * list to be shuffled + * @param end + * element past which shuffling begins + */ + private void shuffle(int[] list, int end) { + int target = 0; + for (int i = list.length - 1; i >= end; i--) { + if (i == 0) { + target = 0; + } else { + target = nextInt(0, i); + } + int temp = list[target]; + list[target] = list[i]; + list[i] = temp; + } + } + + /** + * Returns an array representing n. + * + * @param n + * the natural number to represent + * @return array with entries = elements of n + */ + private int[] getNatural(int n) { + int[] natural = new int[n]; + for (int i = 0; i < n; i++) { + natural[i] = i; + } + return natural; + } + +} diff --git a/src/main/java/org/apache/commons/math/random/RandomGenerator.java b/src/main/java/org/apache/commons/math/random/RandomGenerator.java new file mode 100644 index 0000000..0b86c2a --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/RandomGenerator.java @@ -0,0 +1,148 @@ +/* + * 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.math.random; + + +/** + * Interface extracted from java.util.Random. This interface is + * implemented by {@link AbstractRandomGenerator}. + * + * @since 1.1 + * @version $Revision: 949750 $ $Date: 2010-05-31 16:06:04 +0200 (lun. 31 mai 2010) $ + */ +public interface RandomGenerator { + + /** + * Sets the seed of the underlying random number generator using an + * int seed. + *

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 + * int array seed. + *

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 + * long seed. + *

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 int + * value from this random number generator's sequence. + * All 232 possible int values + * should be produced with (approximately) equal probability. + * + * @return the next pseudorandom, uniformly distributed int + * value from this random number generator's sequence + */ + int nextInt(); + + /** + * Returns a pseudorandom, uniformly distributed 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 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 long + * value from this random number generator's sequence. All + * 264 possible long values + * should be produced with (approximately) equal probability. + * + * @return the next pseudorandom, uniformly distributed long + *value from this random number generator's sequence + */ + long nextLong(); + + /** + * Returns the next pseudorandom, uniformly distributed + * boolean value from this random number generator's + * sequence. + * + * @return the next pseudorandom, uniformly distributed + * boolean value from this random number generator's + * sequence + */ + boolean nextBoolean(); + + /** + * Returns the next pseudorandom, uniformly distributed float + * value between 0.0 and 1.0 from this random + * number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed float + * value between 0.0 and 1.0 from this + * random number generator's sequence + */ + float nextFloat(); + + /** + * Returns the next pseudorandom, uniformly distributed + * double value between 0.0 and + * 1.0 from this random number generator's sequence. + * + * @return the next pseudorandom, uniformly distributed + * double value between 0.0 and + * 1.0 from this random number generator's sequence + */ + double nextDouble(); + + /** + * Returns the next pseudorandom, Gaussian ("normally") distributed + * double value with mean 0.0 and standard + * deviation 1.0 from this random number generator's sequence. + * + * @return the next pseudorandom, Gaussian ("normally") distributed + * double value with mean 0.0 and + * standard deviation 1.0 from this random number + * generator's sequence + */ + double nextGaussian(); +} diff --git a/src/main/java/org/apache/commons/math/random/RandomVectorGenerator.java b/src/main/java/org/apache/commons/math/random/RandomVectorGenerator.java new file mode 100644 index 0000000..15abbd7 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/RandomVectorGenerator.java @@ -0,0 +1,35 @@ +/* + * 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.math.random; + + +/** This interface represents a random generator for whole vectors. + * + * @since 1.2 + * @version $Revision: 811786 $ $Date: 2009-09-06 11:36:08 +0200 (dim. 06 sept. 2009) $ + * + */ + +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/math/random/UncorrelatedRandomVectorGenerator.java b/src/main/java/org/apache/commons/math/random/UncorrelatedRandomVectorGenerator.java new file mode 100644 index 0000000..d365f9b --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/UncorrelatedRandomVectorGenerator.java @@ -0,0 +1,93 @@ +/* + * 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.math.random; + +import java.util.Arrays; + +import org.apache.commons.math.exception.DimensionMismatchException; + +/** + * A {@link RandomVectorGenerator} that generates vectors with uncorrelated + * components. Components of generated vectors follow (independent) Gaussian + * distributions, with parameters supplied in the constructor. + * + * @version $Revision: 962515 $ $Date: 2010-07-09 15:15:28 +0200 (ven. 09 juil. 2010) $ + * @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. + *

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. + *

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/math/random/UniformRandomGenerator.java b/src/main/java/org/apache/commons/math/random/UniformRandomGenerator.java new file mode 100644 index 0000000..54492d0 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/UniformRandomGenerator.java @@ -0,0 +1,61 @@ +/* + * 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.math.random; + +import org.apache.commons.math.util.FastMath; + +/** + * This class implements a normalized uniform random generator. + *

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 + * + * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ + */ + +public class UniformRandomGenerator implements NormalizedRandomGenerator { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 1569292426375546027L; + + /** 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. + *

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/math/random/UnitSphereRandomVectorGenerator.java b/src/main/java/org/apache/commons/math/random/UnitSphereRandomVectorGenerator.java new file mode 100644 index 0000000..cac8f18 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/UnitSphereRandomVectorGenerator.java @@ -0,0 +1,84 @@ +/* + * 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.math.random; + +import org.apache.commons.math.util.FastMath; + + +/** + * Generate random vectors isotropically located on the surface of a sphere. + * + * @since 2.1 + * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ + */ + +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]; + + double normSq; + do { + normSq = 0; + for (int i = 0; i < dimension; i++) { + final double comp = 2 * rand.nextDouble() - 1; + v[i] = comp; + normSq += comp * comp; + } + } while (normSq > 1); + + 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/math/random/ValueServer.java b/src/main/java/org/apache/commons/math/random/ValueServer.java new file mode 100644 index 0000000..9146e69 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/ValueServer.java @@ -0,0 +1,385 @@ +/* + * 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.math.random; +import java.io.BufferedReader; +import java.io.IOException; +import java.io.InputStreamReader; +import java.net.MalformedURLException; +import java.net.URL; + +import org.apache.commons.math.MathRuntimeException; +import org.apache.commons.math.exception.util.LocalizedFormats; + +/** + * Generates values for use in simulation applications. + *

+ * How values are generated is determined by the mode + * property.

+ *

+ * Supported mode values are:

    + *
  • DIGEST_MODE -- uses an empirical distribution
    • + *
  • REPLAY_MODE -- replays data from valuesFileURL
    • + *
  • UNIFORM_MODE -- generates uniformly distributed random values with + * mean = mu
    • + *
  • EXPONENTIAL_MODE -- generates exponentially distributed random values + * with mean = mu
    • + *
  • GAUSSIAN_MODE -- generates Gaussian distributed random values with + * mean = mu and + * standard deviation = sigma
    • + *
  • CONSTANT_MODE -- returns mu every time.

+ * + * @version $Revision: 1003886 $ $Date: 2010-10-02 23:04:44 +0200 (sam. 02 oct. 2010) $ + * + */ +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 RandomData randomData; + + // Data generation modes ====================================== + + /** Creates new ValueServer */ + public ValueServer() { + randomData = new RandomDataImpl(); + } + + /** + * Construct a ValueServer instance using a RandomData as its source + * of random data. + * + * @param randomData the RandomData instance used to source random data + * @since 1.1 + */ + public ValueServer(RandomData randomData) { + this.randomData = randomData; + } + + /** + * 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 + */ + public double getNext() throws IOException { + 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 MathRuntimeException.createIllegalStateException( + 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 + */ + public void fill(double[] values) throws IOException { + for (int i = 0; i < values.length; i++) { + values[i] = getNext(); + } + } + + /** + * Returns an array of length length 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 + */ + public double[] fill(int length) throws IOException { + 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 valuesFileURL, using the default number of bins. + *

+ * valuesFileURL must exist and be + * readable by *this at runtime.

+ *

+ * This method must be called before using getNext() + * with mode = DIGEST_MODE

+ * + * @throws IOException if an I/O error occurs reading the input file + */ + public void computeDistribution() throws IOException { + empiricalDistribution = new EmpiricalDistributionImpl(); + empiricalDistribution.load(valuesFileURL); + } + + /** + * Computes the empirical distribution using values from the file + * in valuesFileURL and binCount bins. + *

+ * valuesFileURL must exist and be readable by this process + * at runtime.

+ *

+ * This method must be called before using getNext() + * with mode = DIGEST_MODE

+ * + * @param binCount the number of bins used in computing the empirical + * distribution + * @throws IOException if an error occurs reading the input file + */ + public void computeDistribution(int binCount) + throws IOException { + empiricalDistribution = new EmpiricalDistributionImpl(binCount); + empiricalDistribution.load(valuesFileURL); + mu = empiricalDistribution.getSampleStats().getMean(); + sigma = empiricalDistribution.getSampleStats().getStandardDeviation(); + } + + /** Getter for property mode. + * @return Value of property mode. + */ + public int getMode() { + return mode; + } + + /** Setter for property mode. + * @param mode New value of property mode. + */ + public void setMode(int mode) { + this.mode = mode; + } + + /** + * Getter for valuesFileURL + * @return Value of property valuesFileURL. + */ + public URL getValuesFileURL() { + return valuesFileURL; + } + + /** + * Sets the valuesFileURL 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 valuesFileURL + * @param url New value of property valuesFileURL. + */ + public void setValuesFileURL(URL url) { + this.valuesFileURL = url; + } + + /** Getter for property empiricalDistribution. + * @return Value of property empiricalDistribution. + */ + public EmpiricalDistribution getEmpiricalDistribution() { + return empiricalDistribution; + } + + /** + * Resets REPLAY_MODE file pointer to the beginning of the valuesFileURL. + * + * @throws IOException if an error occurs opening the file + */ + public void resetReplayFile() throws IOException { + if (filePointer != null) { + try { + filePointer.close(); + filePointer = null; + } catch (IOException ex) { + // ignore + } + } + filePointer = new BufferedReader(new InputStreamReader(valuesFileURL.openStream())); + } + + /** + * Closes 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; + } + } + + /** Getter for property mu. + * @return Value of property mu. + */ + public double getMu() { + return mu; + } + + /** Setter for property mu. + * @param mu New value of property mu. + */ + public void setMu(double mu) { + this.mu = mu; + } + + /** Getter for property sigma. + * @return Value of property sigma. + */ + public double getSigma() { + return sigma; + } + + /** Setter for property sigma. + * @param sigma New value of property sigma. + */ + public void setSigma(double sigma) { + this.sigma = sigma; + } + + //------------- private methods --------------------------------- + + /** + * Gets a random value in DIGEST_MODE. + *

+ * Preconditions:

    + *
  • Before this method is called, computeDistribution() + * must have completed successfully; otherwise an + * IllegalStateException will be thrown

+ * + * @return next random value from the empirical distribution digest + */ + private double getNextDigest() { + if ((empiricalDistribution == null) || + (empiricalDistribution.getBinStats().size() == 0)) { + throw MathRuntimeException.createIllegalStateException(LocalizedFormats.DIGEST_NOT_INITIALIZED); + } + return empiricalDistribution.getNextValue(); + } + + /** + * Gets next sequential value from the valuesFileURL. + *

+ * Throws an IOException if the read fails.

+ *

+ * This method will open the valuesFileURL if there is no + * replay file open.

+ *

+ * The valuesFileURL 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 valuesFileURL is + * empty.

+ * + * @return next value from the replay file + * @throws IOException if there is a problem reading from the file + * @throws NumberFormatException if an invalid numeric string is + * encountered in the file + */ + private double getNextReplay() throws IOException { + 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 MathRuntimeException.createEOFException(LocalizedFormats.URL_CONTAINS_NO_DATA, + valuesFileURL); + } + } + return Double.valueOf(str).doubleValue(); + } + + /** + * Gets a uniformly distributed random value with mean = mu. + * + * @return random uniform value + */ + private double getNextUniform() { + return randomData.nextUniform(0, 2 * mu); + } + + /** + * Gets an exponentially distributed random value with mean = mu. + * + * @return random exponential value + */ + private double getNextExponential() { + return randomData.nextExponential(mu); + } + + /** + * Gets a Gaussian distributed random value with mean = mu + * and standard deviation = sigma. + * + * @return random Gaussian value + */ + private double getNextGaussian() { + return randomData.nextGaussian(mu, sigma); + } + +} diff --git a/src/main/java/org/apache/commons/math/random/Well1024a.java b/src/main/java/org/apache/commons/math/random/Well1024a.java new file mode 100644 index 0000000..8406ed5 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/Well1024a.java @@ -0,0 +1,106 @@ +/* + * 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.math.random; + + +/** This class implements the WELL1024a pseudo-random number generator + * from François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + + *

This generator is described in a paper by François Panneton, + * Pierre L'Ecuyer and Makoto Matsumoto Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2 ACM + * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper + * are in wellrng-errata.txt.

+ + * @see WELL Random number generator + * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $ + * @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. + *

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/math/random/Well19937a.java b/src/main/java/org/apache/commons/math/random/Well19937a.java new file mode 100644 index 0000000..97f9bfd --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/Well19937a.java @@ -0,0 +1,108 @@ +/* + * 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.math.random; + + +/** This class implements the WELL19937a pseudo-random number generator + * from François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + + *

This generator is described in a paper by François Panneton, + * Pierre L'Ecuyer and Makoto Matsumoto Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2 ACM + * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper + * are in wellrng-errata.txt.

+ + * @see WELL Random number generator + * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $ + * @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. + *

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/math/random/Well19937c.java b/src/main/java/org/apache/commons/math/random/Well19937c.java new file mode 100644 index 0000000..53029c1 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/Well19937c.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.math.random; + + +/** This class implements the WELL19937c pseudo-random number generator + * from François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + + *

This generator is described in a paper by François Panneton, + * Pierre L'Ecuyer and Makoto Matsumoto Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2 ACM + * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper + * are in wellrng-errata.txt.

+ + * @see WELL Random number generator + * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $ + * @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. + *

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 ^ ((z4 << 7) & 0xe46e1700); + z4 = z4 ^ ((z4 << 15) & 0x9b868000); + + return z4 >>> (32 - bits); + + } + +} diff --git a/src/main/java/org/apache/commons/math/random/Well44497a.java b/src/main/java/org/apache/commons/math/random/Well44497a.java new file mode 100644 index 0000000..70d672c --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/Well44497a.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.math.random; + + +/** This class implements the WELL44497a pseudo-random number generator + * from François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + + *

This generator is described in a paper by François Panneton, + * Pierre L'Ecuyer and Makoto Matsumoto Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2 ACM + * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper + * are in wellrng-errata.txt.

+ + * @see WELL Random number generator + * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $ + * @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. + *

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/math/random/Well44497b.java b/src/main/java/org/apache/commons/math/random/Well44497b.java new file mode 100644 index 0000000..28b5dbb --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/Well44497b.java @@ -0,0 +1,119 @@ +/* + * 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.math.random; + + +/** This class implements the WELL44497b pseudo-random number generator + * from François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + + *

This generator is described in a paper by François Panneton, + * Pierre L'Ecuyer and Makoto Matsumoto Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2 ACM + * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper + * are in wellrng-errata.txt.

+ + * @see WELL Random number generator + * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $ + * @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. + *

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 ^ ((z4 << 7) & 0x93dd1400); + z4 = z4 ^ ((z4 << 15) & 0xfa118000); + + return z4 >>> (32 - bits); + + } + +} diff --git a/src/main/java/org/apache/commons/math/random/Well512a.java b/src/main/java/org/apache/commons/math/random/Well512a.java new file mode 100644 index 0000000..c822c9d --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/Well512a.java @@ -0,0 +1,107 @@ +/* + * 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.math.random; + + +/** This class implements the WELL512a pseudo-random number generator + * from François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. + + *

This generator is described in a paper by François Panneton, + * Pierre L'Ecuyer and Makoto Matsumoto Improved + * Long-Period Generators Based on Linear Recurrences Modulo 2 ACM + * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper + * are in wellrng-errata.txt.

+ + * @see WELL Random number generator + * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $ + * @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. + *

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/math/random/package.html b/src/main/java/org/apache/commons/math/random/package.html new file mode 100644 index 0000000..9978ca0 --- /dev/null +++ b/src/main/java/org/apache/commons/math/random/package.html @@ -0,0 +1,132 @@ + + + + +

Random number and random data generators.

+

Commons-math provides a few pseudo random number generators. The top level interface is RandomGenerator. + It is implemented by three classes: +

    +
  • {@link org.apache.commons.math.random.JDKRandomGenerator JDKRandomGenerator} + that extends the JDK provided generator
    • +
  • AbstractRandomGenerator as a helper for users generators
    • +
  • BitStreamGenerator which is an abstract class for several generators and + which in turn is extended by: +
      +
    • {@link org.apache.commons.math.random.MersenneTwister MersenneTwister}
      • +
    • {@link org.apache.commons.math.random.Well512a Well512a}
      • +
    • {@link org.apache.commons.math.random.Well1024a Well1024a}
      • +
    • {@link org.apache.commons.math.random.Well19937a Well19937a}
      • +
    • {@link org.apache.commons.math.random.Well19937c Well19937c}
      • +
    • {@link org.apache.commons.math.random.Well44497a Well44497a}
      • +
    • {@link org.apache.commons.math.random.Well44497b Well44497b}
      • +
    +
    • +
+

+ +

+ 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: 219937 - 1 (which is a Mersenne prime). This generator + is described in a paper by Makoto Matsumoto and Takuji Nishimura in 1998: Mersenne Twister: + A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator, 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 2512 - 1 + to 244497 - 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 Improved Long-Period + Generators Based on Linear Recurrences Modulo 2 ACM Transactions on Mathematical Software, + 32, 1 (2006). The errata for the paper are in wellrng-errata.txt. +

+ +

+ 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 2n-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. +

+ +

+ 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 escape zeroland 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). +

+ +

+ 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. +

+ +

+ 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: +

+ +

+ + + + + + + + + + + +
Example of performances
Namegeneration rate (relative to MersenneTwister)
{@link org.apache.commons.math.random.MersenneTwister MersenneTwister}1
{@link org.apache.commons.math.random.JDKRandomGenerator JDKRandomGenerator}between 0.96 and 1.16
{@link org.apache.commons.math.random.Well512a Well512a}between 0.85 and 0.88
{@link org.apache.commons.math.random.Well1024a Well1024a}between 0.63 and 0.73
{@link org.apache.commons.math.random.Well19937a Well19937a}between 0.70 and 0.71
{@link org.apache.commons.math.random.Well19937c Well19937c}between 0.57 and 0.71
{@link org.apache.commons.math.random.Well44497a Well44497a}between 0.69 and 0.71
{@link org.apache.commons.math.random.Well44497b Well44497b}between 0.65 and 0.71
+

+ +

+ So for most simulation problems, the better generators like {@link + org.apache.commons.math.random.Well19937c Well19937c} and {@link + org.apache.commons.math.random.Well44497b Well44497b} are probably very good choices. +

+ +

+ Note that none 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. +

+ + + -- cgit v1.2.3