You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mahout.apache.org by gs...@apache.org on 2009/11/23 16:14:38 UTC
svn commit: r883365 [5/47] - in /lucene/mahout/trunk: ./ examples/ matrix/
matrix/src/ matrix/src/main/ matrix/src/main/java/
matrix/src/main/java/org/ matrix/src/main/java/org/apache/
matrix/src/main/java/org/apache/mahout/ matrix/src/main/java/org/ap...
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Normal.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Normal.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Normal.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Normal.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,149 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random;
+
+import org.apache.mahout.jet.random.engine.RandomEngine;
+import org.apache.mahout.jet.stat.Probability;
+/**
+Normal (aka Gaussian) distribution; See the <A HREF="http://www.cern.ch/RD11/rkb/AN16pp/node188.html#SECTION0001880000000000000000"> math definition</A>
+and <A HREF="http://www.statsoft.com/textbook/glosn.html#Normal Distribution"> animated definition</A>.
+<pre>
+ 1 2
+ pdf(x) = --------- exp( - (x-mean) / 2v )
+ sqrt(2pi*v)
+
+ x
+ -
+ 1 | | 2
+ cdf(x) = --------- | exp( - (t-mean) / 2v ) dt
+ sqrt(2pi*v)| |
+ -
+ -inf.
+</pre>
+where <tt>v = variance = standardDeviation^2</tt>.
+<p>
+Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
+<dt>
+Static methods operate on a default uniform random number generator; they are synchronized.
+<p>
+<b>Implementation:</b> Polar Box-Muller transformation. See
+G.E.P. Box, M.E. Muller (1958): A note on the generation of random normal deviates, Annals Math. Statist. 29, 610-611.
+<p>
+@author wolfgang.hoschek@cern.ch
+@version 1.0, 09/24/99
+*/
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class Normal extends AbstractContinousDistribution {
+ protected double mean;
+ protected double variance;
+ protected double standardDeviation;
+
+ protected double cache; // cache for Box-Mueller algorithm
+ protected boolean cacheFilled; // Box-Mueller
+
+ protected double SQRT_INV; // performance cache
+
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static Normal shared = new Normal(0.0,1.0,makeDefaultGenerator());
+/**
+ * Constructs a normal (gauss) distribution.
+ * Example: mean=0.0, standardDeviation=1.0.
+ */
+public Normal(double mean, double standardDeviation, RandomEngine randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ setState(mean,standardDeviation);
+}
+/**
+ * Returns the cumulative distribution function.
+ */
+public double cdf(double x) {
+ return Probability.normal(mean,variance,x);
+}
+/**
+ * Returns a random number from the distribution.
+ */
+public double nextDouble() {
+ return nextDouble(this.mean,this.standardDeviation);
+}
+/**
+ * Returns a random number from the distribution; bypasses the internal state.
+ */
+public double nextDouble(double mean, double standardDeviation) {
+ // Uses polar Box-Muller transformation.
+ if (cacheFilled && this.mean == mean && this.standardDeviation == standardDeviation) {
+ cacheFilled = false;
+ return cache;
+ };
+
+ double x,y,r,z;
+ do {
+ x = 2.0*randomGenerator.raw() - 1.0;
+ y = 2.0*randomGenerator.raw() - 1.0;
+ r = x*x+y*y;
+ } while (r >= 1.0);
+
+ z = Math.sqrt(-2.0*Math.log(r)/r);
+ cache = mean + standardDeviation*x*z;
+ cacheFilled = true;
+ return mean + standardDeviation*y*z;
+}
+/**
+ * Returns the probability distribution function.
+ */
+public double pdf(double x) {
+ double diff = x-mean;
+ return SQRT_INV * Math.exp(-(diff*diff) / (2.0*variance));
+}
+/**
+ * Sets the uniform random generator internally used.
+ */
+protected void setRandomGenerator(RandomEngine randomGenerator) {
+ super.setRandomGenerator(randomGenerator);
+ this.cacheFilled = false;
+}
+/**
+ * Sets the mean and variance.
+ */
+public void setState(double mean, double standardDeviation) {
+ if (mean!=this.mean || standardDeviation!=this.standardDeviation) {
+ this.mean = mean;
+ this.standardDeviation = standardDeviation;
+ this.variance = standardDeviation*standardDeviation;
+ this.cacheFilled = false;
+
+ this.SQRT_INV = 1.0 / Math.sqrt(2.0*Math.PI*variance);
+ }
+}
+/**
+ * Returns a random number from the distribution with the given mean and standard deviation.
+ */
+public static double staticNextDouble(double mean, double standardDeviation) {
+ synchronized (shared) {
+ return shared.nextDouble(mean,standardDeviation);
+ }
+}
+/**
+ * Returns a String representation of the receiver.
+ */
+public String toString() {
+ return this.getClass().getName()+"("+mean+","+standardDeviation+")";
+}
+/**
+ * Sets the uniform random number generated shared by all <b>static</b> methods.
+ * @param randomGenerator the new uniform random number generator to be shared.
+ */
+private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Normal.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Poisson.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Poisson.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Poisson.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Poisson.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,339 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random;
+
+import org.apache.mahout.jet.math.Arithmetic;
+import org.apache.mahout.jet.random.engine.RandomEngine;
+import org.apache.mahout.jet.stat.Probability;
+/**
+ * Poisson distribution (quick); See the <A HREF="http://www.cern.ch/RD11/rkb/AN16pp/node208.html#SECTION0002080000000000000000"> math definition</A>
+ * and <A HREF="http://www.statsoft.com/textbook/glosp.html#Poisson Distribution"> animated definition</A>.
+ * <p>
+ * <tt>p(k) = (mean^k / k!) * exp(-mean)</tt> for <tt>k >= 0</tt>.
+ * <p>
+ * Valid parameter ranges: <tt>mean > 0</tt>.
+ * Note: if <tt>mean <= 0.0</tt> then always returns zero.
+ * <p>
+ * Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
+ * <dt>
+ * Static methods operate on a default uniform random number generator; they are synchronized.
+ * <p>
+ * <b>Implementation:</b> High performance implementation.
+ * Patchwork Rejection/Inversion method.
+ * <dt>This is a port of <tt>pprsc.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
+ * C-RAND's implementation, in turn, is based upon
+ * <p>
+ * H. Zechner (1994): Efficient sampling from continuous and discrete unimodal distributions,
+ * Doctoral Dissertation, 156 pp., Technical University Graz, Austria.
+ * <p>
+ * Also see
+ * <p>
+ * Stadlober E., H. Zechner (1999), <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">The patchwork rejection method for sampling from unimodal distributions</A>,
+ * to appear in ACM Transactions on Modelling and Simulation.
+ *
+ * @author wolfgang.hoschek@cern.ch
+ * @version 1.0, 09/24/99
+ */
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class Poisson extends AbstractDiscreteDistribution {
+ protected double mean;
+
+ // precomputed and cached values (for performance only)
+ // cache for < SWITCH_MEAN
+ protected double my_old = -1.0;
+ protected double p,q,p0;
+ protected double[] pp = new double[36];
+ protected int llll;
+
+ // cache for >= SWITCH_MEAN
+ protected double my_last = -1.0;
+ protected double ll;
+ protected int k2, k4, k1, k5;
+ protected double dl, dr, r1, r2, r4, r5, lr, l_my, c_pm;
+ protected double f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
+
+ // cache for both;
+ protected int m;
+
+
+ protected static final double MEAN_MAX = Integer.MAX_VALUE; // for all means larger than that, we don't try to compute a poisson deviation, but return the mean.
+ protected static final double SWITCH_MEAN = 10.0; // switch from method A to method B
+
+
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static Poisson shared = new Poisson(0.0,makeDefaultGenerator());
+/**
+ * Constructs a poisson distribution.
+ * Example: mean=1.0.
+ */
+public Poisson(double mean, RandomEngine randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ setMean(mean);
+}
+/**
+ * Returns the cumulative distribution function.
+ */
+public double cdf(int k) {
+ return Probability.poisson(k,this.mean);
+}
+/**
+ * Returns a deep copy of the receiver; the copy will produce identical sequences.
+ * After this call has returned, the copy and the receiver have equal but separate state.
+ *
+ * @return a copy of the receiver.
+ */
+public Object clone() {
+ Poisson copy = (Poisson) super.clone();
+ if (this.pp != null) copy.pp = (double[]) this.pp.clone();
+ return copy;
+}
+private static double f(int k, double l_nu, double c_pm) {
+ return Math.exp(k * l_nu - Arithmetic.logFactorial(k) - c_pm);
+}
+/**
+ * Returns a random number from the distribution.
+ */
+public int nextInt() {
+ return nextInt(this.mean);
+}
+/**
+ * Returns a random number from the distribution; bypasses the internal state.
+ */
+public int nextInt(double theMean) {
+/******************************************************************
+ * *
+ * Poisson Distribution - Patchwork Rejection/Inversion *
+ * *
+ ******************************************************************
+ * *
+ * For parameter my < 10 Tabulated Inversion is applied. *
+ * For my >= 10 Patchwork Rejection is employed: *
+ * The area below the histogram function f(x) is rearranged in *
+ * its body by certain point reflections. Within a large center *
+ * interval variates are sampled efficiently by rejection from *
+ * uniform hats. Rectangular immediate acceptance regions speed *
+ * up the generation. The remaining tails are covered by *
+ * exponential functions. *
+ * *
+ *****************************************************************/
+ RandomEngine gen = this.randomGenerator;
+ double my = theMean;
+
+ double t,g,my_k;
+
+ double gx,gy,px,py,e,x,xx,delta,v;
+ int sign;
+
+ //static double p,q,p0,pp[36];
+ //static long ll,m;
+ double u;
+ int k,i;
+
+ if (my < SWITCH_MEAN) { // CASE B: Inversion- start new table and calculate p0
+ if (my != my_old) {
+ my_old = my;
+ llll = 0;
+ p = Math.exp(-my);
+ q = p;
+ p0 = p;
+ //for (k=pp.length; --k >=0; ) pp[k] = 0;
+ }
+ m = (my > 1.0) ? (int)my : 1;
+ for(;;) {
+ u = gen.raw(); // Step U. Uniform sample
+ k = 0;
+ if (u <= p0) return(k);
+ if (llll != 0) { // Step T. Table comparison
+ i = (u > 0.458) ? Math.min(llll,m) : 1;
+ for (k = i; k <=llll; k++) if (u <= pp[k]) return(k);
+ if (llll == 35) continue;
+ }
+ for (k = llll +1; k <= 35; k++) { // Step C. Creation of new prob.
+ p *= my/(double)k;
+ q += p;
+ pp[k] = q;
+ if (u <= q) {
+ llll = k;
+ return(k);
+ }
+ }
+ llll = 35;
+ }
+ } // end my < SWITCH_MEAN
+ else if (my < MEAN_MAX ) { // CASE A: acceptance complement
+ //static double my_last = -1.0;
+ //static long int m, k2, k4, k1, k5;
+ //static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
+ // f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
+ int Dk, X, Y;
+ double Ds, U, V, W;
+
+ m = (int) my;
+ if (my != my_last) { // set-up
+ my_last = my;
+
+ // approximate deviation of reflection points k2, k4 from my - 1/2
+ Ds = Math.sqrt(my + 0.25);
+
+ // mode m, reflection points k2 and k4, and points k1 and k5, which
+ // delimit the centre region of h(x)
+ k2 = (int) Math.ceil(my - 0.5 - Ds);
+ k4 = (int) (my - 0.5 + Ds);
+ k1 = k2 + k2 - m + 1;
+ k5 = k4 + k4 - m;
+
+ // range width of the critical left and right centre region
+ dl = (double) (k2 - k1);
+ dr = (double) (k5 - k4);
+
+ // recurrence constants r(k) = p(k)/p(k-1) at k = k1, k2, k4+1, k5+1
+ r1 = my / (double) k1;
+ r2 = my / (double) k2;
+ r4 = my / (double)(k4 + 1);
+ r5 = my / (double)(k5 + 1);
+
+ // reciprocal values of the scale parameters of expon. tail envelopes
+ ll = Math.log(r1); // expon. tail left
+ lr = -Math.log(r5); // expon. tail right
+
+ // Poisson constants, necessary for computing function values f(k)
+ l_my = Math.log(my);
+ c_pm = m * l_my - Arithmetic.logFactorial(m);
+
+ // function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5
+ f2 = f(k2, l_my, c_pm);
+ f4 = f(k4, l_my, c_pm);
+ f1 = f(k1, l_my, c_pm);
+ f5 = f(k5, l_my, c_pm);
+
+ // area of the two centre and the two exponential tail regions
+ // area of the two immediate acceptance regions between k2, k4
+ p1 = f2 * (dl + 1.0); // immed. left
+ p2 = f2 * dl + p1; // centre left
+ p3 = f4 * (dr + 1.0) + p2; // immed. right
+ p4 = f4 * dr + p3; // centre right
+ p5 = f1 / ll + p4; // expon. tail left
+ p6 = f5 / lr + p5; // expon. tail right
+ } // end set-up
+
+ for (;;) {
+ // generate uniform number U -- U(0, p6)
+ // case distinction corresponding to U
+ if ((U = gen.raw() * p6) < p2) { // centre left
+
+ // immediate acceptance region R2 = [k2, m) *[0, f2), X = k2, ... m -1
+ if ((V = U - p1) < 0.0) return(k2 + (int)(U/f2));
+ // immediate acceptance region R1 = [k1, k2)*[0, f1), X = k1, ... k2-1
+ if ((W = V / dl) < f1 ) return(k1 + (int)(V/f1));
+
+ // computation of candidate X < k2, and its counterpart Y > k2
+ // either squeeze-acceptance of X or acceptance-rejection of Y
+ Dk = (int)(dl * gen.raw()) + 1;
+ if (W <= f2 - Dk * (f2 - f2/r2)) { // quick accept of
+ return(k2 - Dk); // X = k2 - Dk
+ }
+ if ((V = f2 + f2 - W) < 1.0) { // quick reject of Y
+ Y = k2 + Dk;
+ if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0)) {// quick accept of
+ return(Y); // Y = k2 + Dk
+ }
+ if (V <= f(Y, l_my, c_pm)) return(Y); // final accept of Y
+ }
+ X = k2 - Dk;
+ }
+ else if (U < p4) { // centre right
+ // immediate acceptance region R3 = [m, k4+1)*[0, f4), X = m, ... k4
+ if ((V = U - p3) < 0.0) return(k4 - (int)((U - p2)/f4));
+ // immediate acceptance region R4 = [k4+1, k5+1)*[0, f5)
+ if ((W = V / dr) < f5 ) return(k5 - (int)(V/f5));
+
+ // computation of candidate X > k4, and its counterpart Y < k4
+ // either squeeze-acceptance of X or acceptance-rejection of Y
+ Dk = (int)(dr * gen.raw()) + 1;
+ if (W <= f4 - Dk * (f4 - f4*r4)) { // quick accept of
+ return(k4 + Dk); // X = k4 + Dk
+ }
+ if ((V = f4 + f4 - W) < 1.0) { // quick reject of Y
+ Y = k4 - Dk;
+ if (V <= f4 + Dk * (1.0 - f4)/ dr) { // quick accept of
+ return(Y); // Y = k4 - Dk
+ }
+ if (V <= f(Y, l_my, c_pm)) return(Y); // final accept of Y
+ }
+ X = k4 + Dk;
+ }
+ else {
+ W = gen.raw();
+ if (U < p5) { // expon. tail left
+ Dk = (int)(1.0 - Math.log(W)/ll);
+ if ((X = k1 - Dk) < 0) continue; // 0 <= X <= k1 - 1
+ W *= (U - p4) * ll; // W -- U(0, h(x))
+ if (W <= f1 - Dk * (f1 - f1/r1)) return(X); // quick accept of X
+ }
+ else { // expon. tail right
+ Dk = (int)(1.0 - Math.log(W)/lr);
+ X = k5 + Dk; // X >= k5 + 1
+ W *= (U - p5) * lr; // W -- U(0, h(x))
+ if (W <= f5 - Dk * (f5 - f5*r5)) return(X); // quick accept of X
+ }
+ }
+
+ // acceptance-rejection test of candidate X from the original area
+ // test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)
+ // log f(X) = (X - m)*log(my) - log X! + log m!
+ if (Math.log(W) <= X * l_my - Arithmetic.logFactorial(X) - c_pm) return(X);
+ }
+ }
+ else { // mean is too large
+ return (int) my;
+ }
+}
+/**
+ * Returns the probability distribution function.
+ */
+public double pdf(int k) {
+ return Math.exp(k*Math.log(this.mean) - Arithmetic.logFactorial(k) - this.mean);
+
+ // Overflow sensitive:
+ // return (Math.pow(mean,k) / cephes.Arithmetic.factorial(k)) * Math.exp(-this.mean);
+}
+/**
+ * Sets the mean.
+ */
+public void setMean(double mean) {
+ this.mean = mean;
+}
+/**
+ * Returns a random number from the distribution with the given mean.
+ */
+public static int staticNextInt(double mean) {
+ synchronized (shared) {
+ shared.setMean(mean);
+ return shared.nextInt();
+ }
+}
+/**
+ * Returns a String representation of the receiver.
+ */
+public String toString() {
+ return this.getClass().getName()+"("+mean+")";
+}
+/**
+ * Sets the uniform random number generated shared by all <b>static</b> methods.
+ * @param randomGenerator the new uniform random number generator to be shared.
+ */
+private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Poisson.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/PoissonSlow.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/PoissonSlow.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/PoissonSlow.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/PoissonSlow.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,183 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random;
+
+import org.apache.mahout.jet.random.engine.RandomEngine;
+/**
+ * Poisson distribution; See the <A HREF="http://www.cern.ch/RD11/rkb/AN16pp/node208.html#SECTION0002080000000000000000"> math definition</A>
+ * and <A HREF="http://www.statsoft.com/textbook/glosp.html#Poisson Distribution"> animated definition</A>.
+ * <p>
+ * <tt>p(k) = (mean^k / k!) * exp(-mean)</tt> for <tt>k >= 0</tt>.
+ * <p>
+ * Valid parameter ranges: <tt>mean > 0</tt>.
+ * Note: if <tt>mean <= 0.0</tt> then always returns zero.
+ * <p>
+ * Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
+ * <dt>
+ * Static methods operate on a default uniform random number generator; they are synchronized.
+ * <p>
+ * <b>Implementation:</b>
+ * This is a port of <A HREF="http://wwwinfo.cern.ch/asd/lhc++/clhep/manual/RefGuide/Random/RandPoisson.html">RandPoisson</A> used in <A HREF="http://wwwinfo.cern.ch/asd/lhc++/clhep">CLHEP 1.4.0</A> (C++).
+ * CLHEP's implementation, in turn, is based upon "W.H.Press et al., Numerical Recipes in C, Second Edition".
+ *
+ * @author wolfgang.hoschek@cern.ch
+ * @version 1.0, 09/24/99
+ */
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class PoissonSlow extends AbstractDiscreteDistribution {
+ protected double mean;
+
+ // precomputed and cached values (for performance only)
+ protected double cached_sq;
+ protected double cached_alxm;
+ protected double cached_g;
+
+ protected static final double MEAN_MAX = Integer.MAX_VALUE; // for all means larger than that, we don't try to compute a poisson deviation, but return the mean.
+ protected static final double SWITCH_MEAN = 12.0; // switch from method A to method B
+
+ protected static final double[] cof = { // for method logGamma()
+ 76.18009172947146,-86.50532032941677,
+ 24.01409824083091, -1.231739572450155,
+ 0.1208650973866179e-2, -0.5395239384953e-5};
+
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static PoissonSlow shared = new PoissonSlow(0.0,makeDefaultGenerator());
+/**
+ * Constructs a poisson distribution.
+ * Example: mean=1.0.
+ */
+public PoissonSlow(double mean, RandomEngine randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ setMean(mean);
+}
+/**
+ * Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
+ * xx > 1. For 0 < xx < 1. the reflection formula (6.1.4) can be used first.
+ * (Adapted from Numerical Recipes in C)
+ */
+public static double logGamma(double xx) {
+ double x = xx - 1.0;
+ double tmp = x + 5.5;
+ tmp -= (x + 0.5) * Math.log(tmp);
+ double ser = 1.000000000190015;
+
+ double[] coeff = cof;
+ for (int j = 0; j <= 5; j++ ) {
+ x++;
+ ser += coeff[j]/x;
+ }
+ return -tmp + Math.log(2.5066282746310005*ser);
+}
+/**
+ * Returns a random number from the distribution.
+ */
+public int nextInt() {
+ return nextInt(this.mean);
+}
+/**
+ * Returns a random number from the distribution; bypasses the internal state.
+ */
+private int nextInt(double theMean) {
+ /*
+ * Adapted from "Numerical Recipes in C".
+ */
+ double xm = theMean;
+ double g = this.cached_g;
+
+ if (xm == -1.0 ) return 0; // not defined
+ if (xm < SWITCH_MEAN ) {
+ int poisson = -1;
+ double product = 1;
+ do {
+ poisson++;
+ product *= randomGenerator.raw();
+ } while ( product >= g );
+ // bug in CLHEP 1.4.0: was "} while ( product > g );"
+ return poisson;
+ }
+ else if (xm < MEAN_MAX ) {
+ double t;
+ double em;
+ double sq = this.cached_sq;
+ double alxm = this.cached_alxm;
+
+ RandomEngine rand = this.randomGenerator;
+ do {
+ double y;
+ do {
+ y = Math.tan(Math.PI*rand.raw());
+ em = sq*y + xm;
+ } while (em < 0.0);
+ em = (double) (int)(em); // faster than em = Math.floor(em); (em>=0.0)
+ t = 0.9*(1.0 + y*y)* Math.exp(em*alxm - logGamma(em + 1.0) - g);
+ } while (rand.raw() > t);
+ return (int) em;
+ }
+ else { // mean is too large
+ return (int) xm;
+ }
+}
+/**
+ * Returns a random number from the distribution.
+ */
+protected int nextIntSlow() {
+ final double bound = Math.exp( - mean);
+ int count = 0;
+ double product;
+ for (product = 1.0; product >= bound && product > 0.0; count++) {
+ product *= randomGenerator.raw();
+ }
+ if (product<=0.0 && bound>0.0) return (int) Math.round(mean); // detected endless loop due to rounding errors
+ return count-1;
+}
+/**
+ * Sets the mean.
+ */
+public void setMean(double mean) {
+ if (mean != this.mean) {
+ this.mean = mean;
+ if (mean == -1.0) return; // not defined
+ if (mean < SWITCH_MEAN) {
+ this.cached_g = Math.exp(-mean);
+ }
+ else {
+ this.cached_sq = Math.sqrt(2.0*mean);
+ this.cached_alxm = Math.log(mean);
+ this.cached_g = mean*cached_alxm - logGamma(mean + 1.0);
+ }
+ }
+}
+/**
+ * Returns a random number from the distribution with the given mean.
+ */
+public static int staticNextInt(double mean) {
+ synchronized (shared) {
+ shared.setMean(mean);
+ return shared.nextInt();
+ }
+}
+/**
+ * Returns a String representation of the receiver.
+ */
+public String toString() {
+ return this.getClass().getName()+"("+mean+")";
+}
+/**
+ * Sets the uniform random number generated shared by all <b>static</b> methods.
+ * @param randomGenerator the new uniform random number generator to be shared.
+ */
+private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/PoissonSlow.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Stack.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Stack.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Stack.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Stack.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,87 @@
+/*
+Copyright 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random;
+
+/**
+ * Not yet commented.
+ */
+class Stack {
+ int N; /* max number of elts on stack */
+ int[]v; /* array of values on the stack */
+ int i; /* index of top of stack */
+/**
+ * Constructs a new stack with the given capacity.
+ */
+public Stack(int capacity) {
+ this.N = capacity;
+ this.i = -1; // indicates stack is empty
+ this.v = new int[N];
+/*
+static stack_t *
+new_stack(int N) {
+ stack_t *s;
+ s = (stack_t *)malloc(sizeof(stack_t));
+ s->N = N;
+ s->i = -1; // indicates stack is empty
+ s->v = (int *)malloc(sizeof(int)*N);
+ return s;
+}
+static void
+push_stack(stack_t *s, int v)
+{
+ s->i += 1;
+ if ((s->i) >= (s->N)) {
+ fprintf(stderr,"Cannot push stack!\n");
+ exit(0); // fatal!!
+ }
+ (s->v)[s->i] = v;
+}
+static int pop_stack(stack_t *s)
+{
+ if ((s->i) < 0) {
+ fprintf(stderr,"Cannot pop stack!\n");
+ exit(0);
+ }
+ s->i -= 1;
+ return ((s->v)[s->i + 1]);
+}
+static inline int size_stack(const stack_t *s)
+{
+ return s->i + 1;
+}
+static void free_stack(stack_t *s)
+{
+ free((char *)(s->v));
+ free((char *)s);
+}
+*/
+}
+/**
+ * Returns the topmost element.
+ */
+public int pop() {
+ if (this.i < 0) throw new InternalError("Cannot pop stack!");
+ this.i--;
+ return this.v[this.i+1];
+}
+/**
+ * Places the given value on top of the stack.
+ */
+public void push(int value) {
+ this.i++;
+ if (this.i >= this.N) throw new InternalError("Cannot push stack!");
+ this.v[this.i] = value;
+}
+/**
+ * Returns the number of elements contained.
+ */
+public int size() {
+ return i+1;
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Stack.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/StudentT.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/StudentT.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/StudentT.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/StudentT.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,139 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random;
+
+import org.apache.mahout.jet.random.engine.RandomEngine;
+import org.apache.mahout.jet.stat.Probability;
+/**
+ * StudentT distribution (aka T-distribution); See the <A HREF="http://www.cern.ch/RD11/rkb/AN16pp/node279.html#SECTION0002790000000000000000"> math definition</A>
+ * and <A HREF="http://www.statsoft.com/textbook/gloss.html#Student's t Distribution"> animated definition</A>.
+ * <p>
+ * <tt>p(x) = k * (1+x^2/f) ^ -(f+1)/2</tt> where <tt>k = g((f+1)/2) / (sqrt(pi*f) * g(f/2))</tt> and <tt>g(a)</tt> being the gamma function and <tt>f</tt> being the degrees of freedom.
+ * <p>
+ * Valid parameter ranges: <tt>freedom > 0</tt>.
+ * <p>
+ * Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
+ * <dt>
+ * Static methods operate on a default uniform random number generator; they are synchronized.
+ * <p>
+ * <b>Implementation:</b>
+ * <dt>
+ * Method: Adapted Polar Box-Muller transformation.
+ * <dt>
+ * This is a port of <A HREF="http://wwwinfo.cern.ch/asd/lhc++/clhep/manual/RefGuide/Random/RandStudentT.html">RandStudentT</A> used in <A HREF="http://wwwinfo.cern.ch/asd/lhc++/clhep">CLHEP 1.4.0</A> (C++).
+ * CLHEP's implementation, in turn, is based on <tt>tpol.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
+ * C-RAND's implementation, in turn, is based upon
+ * <p>R.W. Bailey (1994): Polar generation of random variates with the t-distribution, Mathematics of Computation 62, 779-781.
+ *
+ * @author wolfgang.hoschek@cern.ch
+ * @version 1.0, 09/24/99
+ */
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class StudentT extends AbstractContinousDistribution {
+ protected double freedom;
+
+ protected double TERM; // performance cache for pdf()
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static StudentT shared = new StudentT(1.0,makeDefaultGenerator());
+/**
+ * Constructs a StudentT distribution.
+ * Example: freedom=1.0.
+ * @param freedom degrees of freedom.
+ * @throws IllegalArgumentException if <tt>freedom <= 0.0</tt>.
+ */
+public StudentT(double freedom, RandomEngine randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ setState(freedom);
+}
+/**
+ * Returns the cumulative distribution function.
+ */
+public double cdf(double x) {
+ return Probability.studentT(freedom,x);
+}
+/**
+ * Returns a random number from the distribution.
+ */
+public double nextDouble() {
+ return nextDouble(this.freedom);
+}
+/**
+ * Returns a random number from the distribution; bypasses the internal state.
+ * @param a degrees of freedom.
+ * @throws IllegalArgumentException if <tt>a <= 0.0</tt>.
+ */
+public double nextDouble(double degreesOfFreedom) {
+ /*
+ * The polar method of Box/Muller for generating Normal variates
+ * is adapted to the Student-t distribution. The two generated
+ * variates are not independent and the expected no. of uniforms
+ * per variate is 2.5464.
+ *
+ * REFERENCE : - R.W. Bailey (1994): Polar generation of random
+ * variates with the t-distribution, Mathematics
+ * of Computation 62, 779-781.
+ */
+ if (degreesOfFreedom<=0.0) throw new IllegalArgumentException();
+ double u,v,w;
+
+ do {
+ u = 2.0 * randomGenerator.raw() - 1.0;
+ v = 2.0 * randomGenerator.raw() - 1.0;
+ }
+ while ((w = u * u + v * v) > 1.0);
+
+ return(u * Math.sqrt( degreesOfFreedom * ( Math.exp(- 2.0 / degreesOfFreedom * Math.log(w)) - 1.0) / w));
+}
+/**
+ * Returns the probability distribution function.
+ */
+public double pdf(double x) {
+ return this.TERM * Math.pow((1+ x*x/freedom), -(freedom+1) * 0.5);
+}
+/**
+ * Sets the distribution parameter.
+ * @param freedom degrees of freedom.
+ * @throws IllegalArgumentException if <tt>freedom <= 0.0</tt>.
+ */
+public void setState(double freedom) {
+ if (freedom<=0.0) throw new IllegalArgumentException();
+ this.freedom = freedom;
+
+ double val = Fun.logGamma((freedom+1)/2) - Fun.logGamma(freedom/2);
+ this.TERM = Math.exp(val)/Math.sqrt(Math.PI*freedom);
+}
+/**
+ * Returns a random number from the distribution.
+ * @param freedom degrees of freedom.
+ * @throws IllegalArgumentException if <tt>freedom <= 0.0</tt>.
+ */
+public static double staticNextDouble(double freedom) {
+ synchronized (shared) {
+ return shared.nextDouble(freedom);
+ }
+}
+/**
+ * Returns a String representation of the receiver.
+ */
+public String toString() {
+ return this.getClass().getName()+"("+freedom+")";
+}
+/**
+ * Sets the uniform random number generated shared by all <b>static</b> methods.
+ * @param randomGenerator the new uniform random number generator to be shared.
+ */
+private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/StudentT.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Uniform.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Uniform.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Uniform.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Uniform.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,226 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random;
+
+import org.apache.mahout.jet.random.engine.RandomEngine;
+/**
+ * Uniform distribution; <A HREF="http://www.cern.ch/RD11/rkb/AN16pp/node292.html#SECTION0002920000000000000000"> Math definition</A>
+ * and <A HREF="http://www.statsoft.com/textbook/glosu.html#Uniform Distribution"> animated definition</A>.
+ * <p>
+ * Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
+ * <dt>
+ * Static methods operate on a default uniform random number generator; they are synchronized.
+ * <p>
+ * @author wolfgang.hoschek@cern.ch
+ * @version 1.0, 09/24/99
+ */
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class Uniform extends AbstractContinousDistribution {
+ protected double min;
+ protected double max;
+
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static Uniform shared = new Uniform(makeDefaultGenerator());
+/**
+ * Constructs a uniform distribution with the given minimum and maximum, using a {@link org.apache.mahout.jet.random.engine.MersenneTwister} seeded with the given seed.
+ */
+public Uniform(double min, double max, int seed) {
+ this(min,max, new org.apache.mahout.jet.random.engine.MersenneTwister(seed));
+}
+/**
+ * Constructs a uniform distribution with the given minimum and maximum.
+ */
+public Uniform(double min, double max, RandomEngine randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ setState(min,max);
+}
+/**
+ * Constructs a uniform distribution with <tt>min=0.0</tt> and <tt>max=1.0</tt>.
+ */
+public Uniform(RandomEngine randomGenerator) {
+ this(0,1,randomGenerator);
+}
+/**
+ * Returns the cumulative distribution function (assuming a continous uniform distribution).
+ */
+public double cdf(double x) {
+ if (x <= min) return 0.0;
+ if (x >= max) return 1.0;
+ return (x-min) / (max-min);
+}
+/**
+ * Returns a uniformly distributed random <tt>boolean</tt>.
+ */
+public boolean nextBoolean() {
+ return randomGenerator.raw() > 0.5;
+}
+/**
+ * Returns a uniformly distributed random number in the open interval <tt>(min,max)</tt> (excluding <tt>min</tt> and <tt>max</tt>).
+ */
+public double nextDouble() {
+ return min+(max-min)*randomGenerator.raw();
+}
+/**
+ * Returns a uniformly distributed random number in the open interval <tt>(from,to)</tt> (excluding <tt>from</tt> and <tt>to</tt>).
+ * Pre conditions: <tt>from <= to</tt>.
+ */
+public double nextDoubleFromTo(double from, double to) {
+ return from+(to-from)*randomGenerator.raw();
+}
+/**
+ * Returns a uniformly distributed random number in the open interval <tt>(from,to)</tt> (excluding <tt>from</tt> and <tt>to</tt>).
+ * Pre conditions: <tt>from <= to</tt>.
+ */
+public float nextFloatFromTo(float from, float to) {
+ return (float) nextDoubleFromTo(from,to);
+}
+/**
+ * Returns a uniformly distributed random number in the closed interval <tt>[min,max]</tt> (including <tt>min</tt> and <tt>max</tt>).
+ */
+public int nextInt() {
+ return nextIntFromTo((int)Math.round(min), (int)Math.round(max));
+}
+/**
+ * Returns a uniformly distributed random number in the closed interval <tt>[from,to]</tt> (including <tt>from</tt> and <tt>to</tt>).
+ * Pre conditions: <tt>from <= to</tt>.
+ */
+public int nextIntFromTo(int from, int to) {
+ return (int) ((long)from + (long)((1L + (long)to - (long)from)*randomGenerator.raw()));
+}
+/**
+ * Returns a uniformly distributed random number in the closed interval <tt>[from,to]</tt> (including <tt>from</tt> and <tt>to</tt>).
+ * Pre conditions: <tt>from <= to</tt>.
+ */
+public long nextLongFromTo(long from, long to) {
+ /* Doing the thing turns out to be more tricky than expected.
+ avoids overflows and underflows.
+ treats cases like from=-1, to=1 and the like right.
+ the following code would NOT solve the problem: return (long) (Doubles.randomFromTo(from,to));
+
+ rounding avoids the unsymmetric behaviour of casts from double to long: (long) -0.7 = 0, (long) 0.7 = 0.
+ checking for overflows and underflows is also necessary.
+ */
+
+ // first the most likely and also the fastest case.
+ if (from>=0 && to<Long.MAX_VALUE) {
+ return from + (long) (nextDoubleFromTo(0.0,to-from+1));
+ }
+
+ // would we get a numeric overflow?
+ // if not, we can still handle the case rather efficient.
+ double diff = ((double)to) - (double)from + 1.0;
+ if (diff <= Long.MAX_VALUE) {
+ return from + (long) (nextDoubleFromTo(0.0,diff));
+ }
+
+ // now the pathologic boundary cases.
+ // they are handled rather slow.
+ long random;
+ if (from==Long.MIN_VALUE) {
+ if (to==Long.MAX_VALUE) {
+ //return Math.round(nextDoubleFromTo(from,to));
+ int i1 = nextIntFromTo(Integer.MIN_VALUE,Integer.MAX_VALUE);
+ int i2 = nextIntFromTo(Integer.MIN_VALUE,Integer.MAX_VALUE);
+ return ((i1 & 0xFFFFFFFFL) << 32) | (i2 & 0xFFFFFFFFL);
+ }
+ random = Math.round(nextDoubleFromTo(from,to+1));
+ if (random > to) random = from;
+ }
+ else {
+ random = Math.round(nextDoubleFromTo(from-1,to));
+ if (random < from) random = to;
+ }
+ return random;
+}
+/**
+ * Returns the probability distribution function (assuming a continous uniform distribution).
+ */
+public double pdf(double x) {
+ if (x <= min || x >= max) return 0.0;
+ return 1.0 / (max-min);
+}
+/**
+ * Sets the internal state.
+ */
+public void setState(double min, double max) {
+ if (max<min) { setState(max,min); return; }
+ this.min=min;
+ this.max=max;
+}
+/**
+ * Returns a uniformly distributed random <tt>boolean</tt>.
+ */
+public static boolean staticNextBoolean() {
+ synchronized (shared) {
+ return shared.nextBoolean();
+ }
+}
+/**
+ * Returns a uniformly distributed random number in the open interval <tt>(0,1)</tt> (excluding <tt>0</tt> and <tt>1</tt>).
+ */
+public static double staticNextDouble() {
+ synchronized (shared) {
+ return shared.nextDouble();
+ }
+}
+/**
+ * Returns a uniformly distributed random number in the open interval <tt>(from,to)</tt> (excluding <tt>from</tt> and <tt>to</tt>).
+ * Pre conditions: <tt>from <= to</tt>.
+ */
+public static double staticNextDoubleFromTo(double from, double to) {
+ synchronized (shared) {
+ return shared.nextDoubleFromTo(from,to);
+ }
+}
+/**
+ * Returns a uniformly distributed random number in the open interval <tt>(from,to)</tt> (excluding <tt>from</tt> and <tt>to</tt>).
+ * Pre conditions: <tt>from <= to</tt>.
+ */
+public static float staticNextFloatFromTo(float from, float to) {
+ synchronized (shared) {
+ return shared.nextFloatFromTo(from,to);
+ }
+}
+/**
+ * Returns a uniformly distributed random number in the closed interval <tt>[from,to]</tt> (including <tt>from</tt> and <tt>to</tt>).
+ * Pre conditions: <tt>from <= to</tt>.
+ */
+public static int staticNextIntFromTo(int from, int to) {
+ synchronized (shared) {
+ return shared.nextIntFromTo(from,to);
+ }
+}
+/**
+ * Returns a uniformly distributed random number in the closed interval <tt>[from,to]</tt> (including <tt>from</tt> and <tt>to</tt>).
+ * Pre conditions: <tt>from <= to</tt>.
+ */
+public static long staticNextLongFromTo(long from, long to) {
+ synchronized (shared) {
+ return shared.nextLongFromTo(from,to);
+ }
+}
+/**
+ * Sets the uniform random number generation engine shared by all <b>static</b> methods.
+ * @param randomGenerator the new uniform random number generation engine to be shared.
+ */
+public static void staticSetRandomEngine(RandomEngine randomGenerator) {
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
+}
+/**
+ * Returns a String representation of the receiver.
+ */
+public String toString() {
+ return this.getClass().getName()+"("+min+","+max+")";
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Uniform.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/VonMises.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/VonMises.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/VonMises.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/VonMises.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,139 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random;
+
+import org.apache.mahout.jet.random.engine.RandomEngine;
+/**
+ * Von Mises distribution.
+ * <p>
+ * Valid parameter ranges: <tt>k > 0</tt>.
+ * <p>
+ * Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
+ * <dt>
+ * Static methods operate on a default uniform random number generator; they are synchronized.
+ * <p>
+ * <b>Implementation:</b>
+ * <dt>
+ * Method: Acceptance Rejection.
+ * <dt>
+ * This is a port of <tt>mwc.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
+ * C-RAND's implementation, in turn, is based upon
+ * <p>
+ * D.J. Best, N.I. Fisher (1979): Efficient simulation of the von Mises distribution, Appl. Statist. 28, 152-157.
+ *
+ * @author wolfgang.hoschek@cern.ch
+ * @version 1.0, 09/24/99
+ */
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class VonMises extends AbstractContinousDistribution {
+ protected double my_k;
+
+ // cached vars for method nextDouble(a) (for performance only)
+ private double k_set = -1.0;
+ private double tau,rho,r;
+
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static VonMises shared = new VonMises(1.0,makeDefaultGenerator());
+/**
+ * Constructs a Von Mises distribution.
+ * Example: k=1.0.
+ * @throws IllegalArgumentException if <tt>k <= 0.0</tt>.
+ */
+public VonMises(double freedom, RandomEngine randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ setState(freedom);
+}
+/**
+ * Returns a random number from the distribution.
+ */
+public double nextDouble() {
+ return nextDouble(this.my_k);
+}
+/**
+ * Returns a random number from the distribution; bypasses the internal state.
+ * @throws IllegalArgumentException if <tt>k <= 0.0</tt>.
+ */
+public double nextDouble(double k) {
+/******************************************************************
+ * *
+ * Von Mises Distribution - Acceptance Rejection *
+ * *
+ ******************************************************************
+ * *
+ * FUNCTION : - mwc samples a random number from the von Mises *
+ * distribution ( -Pi <= x <= Pi) with parameter *
+ * k > 0 via rejection from the wrapped Cauchy *
+ * distibution. *
+ * REFERENCE: - D.J. Best, N.I. Fisher (1979): Efficient *
+ * simulation of the von Mises distribution, *
+ * Appl. Statist. 28, 152-157. *
+ * SUBPROGRAM: - drand(seed) ... (0,1)-Uniform generator with *
+ * unsigned long integer *seed. *
+ * *
+ * Implemented by F. Niederl, August 1992 *
+ ******************************************************************/
+ double u,v,w,c,z;
+
+ if (k<=0.0) throw new IllegalArgumentException();
+
+ if (k_set!=k) { // SET-UP
+ tau = 1.0 + Math.sqrt(1.0 + 4.0*k*k);
+ rho = (tau-Math.sqrt(2.0*tau)) / (2.0*k);
+ r = (1.0+rho*rho) / (2.0*rho);
+ k_set = k;
+ }
+
+ // GENERATOR
+ do {
+ u = randomGenerator.raw(); // U(0/1)
+ v = randomGenerator.raw(); // U(0/1)
+ z = Math.cos(Math.PI * u);
+ w = (1.0+r*z) / (r+z);
+ c = k*(r-w);
+ } while ((c*(2.0-c) < v) && (Math.log(c/v)+1.0 < c)); // Acceptance/Rejection
+
+ return (randomGenerator.raw() > 0.5)? Math.acos(w): -Math.acos(w); // Random sign //
+ // 0 <= x <= Pi : -Pi <= x <= 0 //
+}
+/**
+ * Sets the distribution parameter.
+ * @throws IllegalArgumentException if <tt>k <= 0.0</tt>.
+ */
+public void setState(double k) {
+ if (k<=0.0) throw new IllegalArgumentException();
+ this.my_k = k;
+}
+/**
+ * Returns a random number from the distribution.
+ * @throws IllegalArgumentException if <tt>k <= 0.0</tt>.
+ */
+public static double staticNextDouble(double freedom) {
+ synchronized (shared) {
+ return shared.nextDouble(freedom);
+ }
+}
+/**
+ * Returns a String representation of the receiver.
+ */
+public String toString() {
+ return this.getClass().getName()+"("+my_k+")";
+}
+/**
+ * Sets the uniform random number generated shared by all <b>static</b> methods.
+ * @param randomGenerator the new uniform random number generator to be shared.
+ */
+private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/VonMises.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Zeta.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Zeta.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Zeta.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Zeta.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,164 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random;
+
+import org.apache.mahout.jet.random.engine.RandomEngine;
+/**
+ * Zeta distribution.
+ * <p>
+ * Valid parameter ranges: <tt>ro > 0</tt> and <tt>pk >= 0</tt>.
+ * <dt>
+ * If either <tt>ro > 100</tt> or <tt>k > 10000</tt> numerical problems in
+ * computing the theoretical moments arise, therefore <tt>ro <= 100</tt> and
+ * <tt>k <= 10000</tt> are recommended.
+ * <p>
+ * Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
+ * <dt>
+ * Static methods operate on a default uniform random number generator; they are synchronized.
+ * <p>
+ * <b>Implementation:</b>
+ * <dt>Method: Acceptance/Rejection.
+ * High performance implementation.
+ * <dt>This is a port and adaption of <tt>Zeta.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
+ * C-RAND's implementation, in turn, is based upon
+ * <p>
+ * J. Dagpunar (1988): Principles of Random Variate Generation, Clarendon Press, Oxford.
+ *
+ * @author wolfgang.hoschek@cern.ch
+ * @version 1.0, 09/24/99
+ */
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class Zeta extends AbstractDiscreteDistribution {
+ protected double ro;
+ protected double pk;
+
+ // cached values (for performance)
+ protected double c,d,ro_prev = -1.0,pk_prev = -1.0;
+ protected double maxlongint = Long.MAX_VALUE - 1.5;
+
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static Zeta shared = new Zeta(1.0,1.0,makeDefaultGenerator());
+/**
+ * Constructs a Zeta distribution.
+ */
+public Zeta(double ro, double pk, RandomEngine randomGenerator) {
+ setRandomGenerator(randomGenerator);
+ setState(ro,pk);
+}
+/**
+ * Returns a zeta distributed random number.
+ */
+protected long generateZeta(double ro, double pk, RandomEngine randomGenerator) {
+/******************************************************************
+ * *
+ * Zeta Distribution - Acceptance Rejection *
+ * *
+ ******************************************************************
+ * *
+ * To sample from the Zeta distribution with parameters ro and pk *
+ * it suffices to sample variates x from the distribution with *
+ * density function f(x)=B*{[x+0.5]+pk}^(-(1+ro)) ( x > .5 ) *
+ * and then deliver k=[x+0.5]. *
+ * 1/B=Sum[(j+pk)^-(ro+1)] (j=1,2,...) converges for ro >= .5 . *
+ * It is not necessary to compute B, because variates x are *
+ * generated by acceptance rejection using density function *
+ * g(x)=ro*(c+0.5)^ro*(c+x)^-(ro+1). *
+ * * *
+ * Integer overflow is possible, when ro is small (ro <= .5) and *
+ * pk large. In this case a new sample is generated. If ro and pk *
+ * satisfy the inequality ro > .14 + pk*1.85e-8 + .02*ln(pk) *
+ * the percentage of overflow is less than 1%, so that the *
+ * result is reliable. *
+ * NOTE: The comment above is likely to be nomore valid since *
+ * the C-version operated on 32-bit integers, while this Java *
+ * version operates on 64-bit integers. However, the following is *
+ * still valid: * *
+ * * *
+ * If either ro > 100 or k > 10000 numerical problems in *
+ * computing the theoretical moments arise, therefore ro<=100 and *
+ * k<=10000 are recommended. *
+ * *
+ ******************************************************************
+ * *
+ * FUNCTION: - zeta samples a random number from the *
+ * Zeta distribution with parameters ro > 0 and *
+ * pk >= 0. *
+ * REFERENCE: - J. Dagpunar (1988): Principles of Random *
+ * Variate Generation, Clarendon Press, Oxford. *
+ * *
+ ******************************************************************/
+ double u,v,e,x;
+ long k;
+
+ if (ro != ro_prev || pk != pk_prev) { // Set-up
+ ro_prev = ro;
+ pk_prev = pk;
+ if (ro<pk) {
+ c = pk-0.5;
+ d = 0;
+ }
+ else {
+ c = ro-0.5;
+ d = (1.0+ro)*Math.log((1.0+pk)/(1.0+ro));
+ }
+ }
+ do {
+ do {
+ u=randomGenerator.raw();
+ v=randomGenerator.raw();
+ x = (c+0.5)*Math.exp(-Math.log(u)/ro) - c;
+ } while (x<=0.5 || x>=maxlongint);
+
+ k = (int) (x+0.5);
+ e = -Math.log(v);
+ } while ( e < (1.0+ro)*Math.log((k+pk)/(x+c)) - d );
+
+ return k;
+}
+/**
+ * Returns a random number from the distribution.
+ */
+public int nextInt() {
+ return (int) generateZeta(ro, pk, randomGenerator);
+}
+/**
+ * Sets the parameters.
+ */
+public void setState(double ro, double pk) {
+ this.ro = ro;
+ this.pk = pk;
+}
+/**
+ * Returns a random number from the distribution.
+ */
+public static int staticNextInt(double ro, double pk) {
+ synchronized (shared) {
+ shared.setState(ro,pk);
+ return shared.nextInt();
+ }
+}
+/**
+ * Returns a String representation of the receiver.
+ */
+public String toString() {
+ return this.getClass().getName()+"("+ro+","+pk+")";
+}
+/**
+ * Sets the uniform random number generated shared by all <b>static</b> methods.
+ * @param randomGenerator the new uniform random number generator to be shared.
+ */
+private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Zeta.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/Benchmark.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/Benchmark.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/Benchmark.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/Benchmark.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,253 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random.engine;
+
+/**
+ * Benchmarks the performance of the currently provided uniform pseudo-random number generation engines.
+ * <p>
+ * All distributions are obtained by using a <b>uniform</b> pseudo-random number generation engine.
+ * followed by a transformation to the desired distribution.
+ * Therefore, the performance of the uniform engines is crucial.
+ * <p>
+ * <h2 align=center>Comparison of uniform generation engines</h2>
+ * <center>
+ * <table border>
+ * <tr>
+ * <td align="center" width="40%">Name</td>
+ * <td align="center" width="20%">Period</td>
+ * <td align="center" width="40%">
+ * <p>Speed<br>
+ * [# million uniform random numbers generated/sec]<br>
+ * Pentium Pro 200 Mhz, JDK 1.2, NT</p>
+ * </td>
+ * </tr>
+ * <tr>
+ * <td align="center" width="40%"> <tt>MersenneTwister</tt></td>
+ * <td align="center" width="20%">2<sup>19937</sup>-1 (=10<sup>6001</sup>)</td>
+ * <td align="center" width="40">2.5</td>
+ * </tr>
+ * <tr>
+ * <td align="center" width="40%"> <tt>Ranlux</tt> (default luxury level 3) </td>
+ * <td align="center" width="20%">10<sup>171</sup></td>
+ * <td align="center" width="40">0.4</td>
+ * </tr>
+ * <tr>
+ * <td align="center" width="40"> <tt>Ranmar</tt></td>
+ * <td align="center" width="20">10<sup>43</sup></td>
+ * <td align="center" width="40%">1.6</td>
+ * </tr>
+ * <tr>
+ * <td align="center" width="40%"> <tt>Ranecu</tt> </td>
+ * <td align="center" width="20">10<sup>18</sup></td>
+ * <td align="center" width="40%">1.5</td>
+ * </tr>
+ * <tr>
+ * <td align="center"> <tt>java.util.Random.nextFloat() </tt><tt>
+ * </tt></td>
+ * <td align="center"><font size=+3>?</font></td>
+ * <td align="center">2.4</td>
+ * </tr>
+ * </table>
+ * </center>
+ * <p>
+ * <b>Note:</b> Methods working on the default uniform random generator are <b>synchronized</b> and therefore in current VM's <b>slow</b> (as of June '99).
+ * Methods taking as argument a uniform random generator are <b>not synchronized</b> and therefore much <b>quicker</b>.
+ * Thus, if you need a lot of random numbers, you should use the unsynchronized approach:
+ * <p>
+ * <b>Example usage:</b><pre>
+ * edu.cornell.lassp.houle.RngPack.RandomElement generator;
+ * generator = new org.apache.mahout.jet.random.engine.MersenneTwister(new java.util.Date());
+ * //generator = new edu.cornell.lassp.houle.RngPack.Ranecu(new java.util.Date());
+ * //generator = new edu.cornell.lassp.houle.RngPack.Ranmar(new java.util.Date());
+ * //generator = new edu.cornell.lassp.houle.RngPack.Ranlux(new java.util.Date());
+ * //generator = makeDefaultGenerator();
+ * for (int i=1000000; --i >=0; ) {
+ * double uniform = generator.raw();
+ * ...
+ * }
+ * </pre>
+ *
+ *
+ * @see org.apache.mahout.jet.random
+ * @author wolfgang.hoschek@cern.ch
+ * @version 1.0, 09/24/99
+ */
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class Benchmark {
+/**
+ * Makes this class non instantiable, but still let's others inherit from it.
+ */
+protected Benchmark() {
+ throw new RuntimeException("Non instantiable");
+}
+/**
+ * Benchmarks <tt>raw()</tt> for various uniform generation engines.
+ */
+public static void benchmark(int times) {
+ org.apache.mahout.colt.Timer timer = new org.apache.mahout.colt.Timer();
+ RandomEngine gen;
+
+ timer.reset().start();
+ for (int i=times; --i>=0; ) ; // no operation
+ timer.stop().display();
+ float emptyLoop = timer.elapsedTime();
+ System.out.println("empty loop timing done.");
+
+ gen = new MersenneTwister();
+ System.out.println("\n MersenneTwister:");
+ timer.reset().start();
+ for (int i=times; --i>=0; ) gen.raw();
+ timer.stop().display();
+ System.out.println(times/(timer.elapsedTime()-emptyLoop)+ " numbers per second.");
+
+
+ gen = new MersenneTwister64();
+ System.out.println("\n MersenneTwister64:");
+ timer.reset().start();
+ for (int i=times; --i>=0; ) gen.raw();
+ timer.stop().display();
+ System.out.println(times/(timer.elapsedTime()-emptyLoop)+ " numbers per second.");
+
+ /*
+ gen = new edu.stanford.mt.MersenneTwister();
+ System.out.println("\n edu.stanford.mt.MersenneTwister:");
+ timer.reset().start();
+ for (int i=times; --i>=0; ) gen.raw();
+ timer.stop().display();
+ System.out.println(times/(timer.elapsedTime()-emptyLoop)+ " numbers per second.");
+ */
+
+
+ gen = new DRand();
+ System.out.println("\nDRand:");
+ timer.reset().start();
+ for (int i=times; --i>=0; ) gen.raw();
+ timer.stop().display();
+ System.out.println(times/(timer.elapsedTime()-emptyLoop)+ " numbers per second.");
+
+
+ java.util.Random javaGen = new java.util.Random();
+ System.out.println("\njava.util.Random.nextFloat():");
+ timer.reset().start();
+ for (int i=times; --i>=0; ) javaGen.nextFloat(); // nextDouble() is slower
+ timer.stop().display();
+ System.out.println(times/(timer.elapsedTime()-emptyLoop)+ " numbers per second.");
+
+ /*
+ gen = new edu.cornell.lassp.houle.RngPack.Ranecu();
+ System.out.println("\nRanecu:");
+ timer.reset().start();
+ for (int i=times; --i>=0; ) gen.raw();
+ timer.stop().display();
+ System.out.println(times/(timer.elapsedTime()-emptyLoop)+ " numbers per second.");
+
+ gen = new edu.cornell.lassp.houle.RngPack.Ranmar();
+ System.out.println("\nRanmar:");
+ timer.reset().start();
+ for (int i=times; --i>=0; ) gen.raw();
+ timer.stop().display();
+ System.out.println(times/(timer.elapsedTime()-emptyLoop)+ " numbers per second.");
+
+ gen = new edu.cornell.lassp.houle.RngPack.Ranlux();
+ System.out.println("\nRanlux:");
+ timer.reset().start();
+ for (int i=times; --i>=0; ) gen.raw();
+ timer.stop().display();
+ System.out.println(times/(timer.elapsedTime()-emptyLoop)+ " numbers per second.");
+ */
+
+ System.out.println("\nGood bye.\n");
+
+}
+/**
+ * Tests various methods of this class.
+ */
+public static void main(String args[]) {
+ long from = Long.parseLong(args[0]);
+ long to = Long.parseLong(args[1]);
+ int times = Integer.parseInt(args[2]);
+ int runs = Integer.parseInt(args[3]);
+ //testRandomFromTo(from,to,times);
+ //benchmark(1000000);
+ //benchmark(1000000);
+ for (int i=0; i<runs; i++) {
+ benchmark(times);
+ //benchmarkSync(times);
+ }
+}
+/**
+ * Prints the first <tt>size</tt> random numbers generated by the given engine.
+ */
+public static void test(int size, RandomEngine randomEngine) {
+ RandomEngine random;
+
+ /*
+ System.out.println("raw():");
+ random = (RandomEngine) randomEngine.clone();
+ //org.apache.mahout.colt.Timer timer = new org.apache.mahout.colt.Timer().start();
+ for (int j=0, i=size; --i>=0; j++) {
+ System.out.print(" "+random.raw());
+ if (j%8==7) System.out.println();
+ }
+
+ System.out.println("\n\nfloat():");
+ random = (RandomEngine) randomEngine.clone();
+ for (int j=0, i=size; --i>=0; j++) {
+ System.out.print(" "+random.nextFloat());
+ if (j%8==7) System.out.println();
+ }
+
+ System.out.println("\n\ndouble():");
+ random = (RandomEngine) randomEngine.clone();
+ for (int j=0, i=size; --i>=0; j++) {
+ System.out.print(" "+random.nextDouble());
+ if (j%8==7) System.out.println();
+ }
+ */
+ System.out.println("\n\nint():");
+ random = (RandomEngine) randomEngine.clone();
+ for (int j=0, i=size; --i>=0; j++) {
+ System.out.print(" "+random.nextInt());
+ if (j%8==7) System.out.println();
+ }
+
+ //timer.stop().display();
+ System.out.println("\n\nGood bye.\n");
+}
+/**
+ * Tests various methods of this class.
+ */
+private static void xtestRandomFromTo(long from, long to, int times) {
+ System.out.println("from="+from+", to="+to);
+
+ //org.apache.mahout.colt.set.OpenMultiFloatHashSet multiset = new org.apache.mahout.colt.set.OpenMultiFloatHashSet();
+
+ java.util.Random randomJava = new java.util.Random();
+ //edu.cornell.lassp.houle.RngPack.RandomElement random = new edu.cornell.lassp.houle.RngPack.Ranecu();
+ //edu.cornell.lassp.houle.RngPack.RandomElement random = new edu.cornell.lassp.houle.RngPack.MT19937B();
+ //edu.cornell.lassp.houle.RngPack.RandomElement random = new edu.stanford.mt.MersenneTwister();
+ RandomEngine random = new MersenneTwister();
+ int _from=(int)from, _to=(int)to;
+ org.apache.mahout.colt.Timer timer = new org.apache.mahout.colt.Timer().start();
+ for (int j=0, i=times; --i>=0; j++) {
+ //randomJava.nextInt(10000);
+ //Integers.randomFromTo(_from,_to);
+ System.out.print(" "+random.raw());
+ if (j%8==7) System.out.println();
+ //multiset.add(nextIntFromTo(_from,_to));
+ }
+
+ timer.stop().display();
+ //System.out.println(multiset); //check the distribution
+ System.out.println("Good bye.\n");
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/Benchmark.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/DRand.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/DRand.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/DRand.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/DRand.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,91 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random.engine;
+
+import java.util.Date;
+/**
+ * Quick medium quality uniform pseudo-random number generator.
+ *
+ * Produces uniformly distributed <tt>int</tt>'s and <tt>long</tt>'s in the closed intervals <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt> and <tt>[Long.MIN_VALUE,Long.MAX_VALUE]</tt>, respectively,
+ * as well as <tt>float</tt>'s and <tt>double</tt>'s in the open unit intervals <tt>(0.0f,1.0f)</tt> and <tt>(0.0,1.0)</tt>, respectively.
+ * <p>
+ * The seed can be any integer satisfying <tt>0 < 4*seed+1 < 2<sup>32</sup></tt>.
+ * In other words, there must hold <tt>seed >= 0 && seed < 1073741823</tt>.
+ * <p>
+ * <b>Quality:</b> This generator follows the multiplicative congruential method of the form
+ * <dt>
+ * <tt>z(i+1) = a * z(i) (mod m)</tt> with
+ * <tt>a=663608941 (=0X278DDE6DL), m=2<sup>32</sup></tt>.
+ * <dt>
+ * <tt>z(i)</tt> assumes all different values <tt>0 < 4*seed+1 < m</tt> during a full period of 2<sup>30</sup>.
+ *
+ * <p>
+ * <b>Performance:</b> TO_DO
+ * <p>
+ * <b>Implementation:</b> TO_DO
+ * <p>
+ * Note that this implementation is <b>not synchronized</b>.
+ * <p>
+ *
+ * @author wolfgang.hoschek@cern.ch
+ * @version 1.0, 09/24/99
+ * @see MersenneTwister
+ * @see java.util.Random
+ */
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class DRand extends RandomEngine {
+ private int current;
+ public static final int DEFAULT_SEED = 1;
+/**
+ * Constructs and returns a random number generator with a default seed, which is a <b>constant</b>.
+ */
+public DRand() {
+ this(DEFAULT_SEED);
+}
+/**
+ * Constructs and returns a random number generator with the given seed.
+ * @param seed should not be 0, in such a case <tt>DRand.DEFAULT_SEED</tt> is substituted.
+ */
+public DRand(int seed) {
+ setSeed(seed);
+}
+/**
+ * Constructs and returns a random number generator seeded with the given date.
+ *
+ * @param d typically <tt>new java.util.Date()</tt>
+ */
+public DRand(Date d) {
+ this((int)d.getTime());
+}
+/**
+ * Returns a 32 bit uniformly distributed random number in the closed interval <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt> (including <tt>Integer.MIN_VALUE</tt> and <tt>Integer.MAX_VALUE</tt>).
+ */
+public int nextInt() {
+ current *= 0x278DDE6D; /* z(i+1)=a*z(i) (mod 2**32) */
+ // a == 0x278DDE6D == 663608941
+
+ return current;
+}
+/**
+ * Sets the receiver's seed.
+ * This method resets the receiver's entire internal state.
+ * The following condition must hold: <tt>seed >= 0 && seed < (2<sup>32</sup>-1) / 4</tt>.
+ * @param seed if the above condition does not hold, a modified seed that meets the condition is silently substituted.
+ */
+protected void setSeed(int seed) {
+ if (seed<0) seed = -seed;
+ int limit = (int)((Math.pow(2,32)-1) /4); // --> 536870911
+ if (seed >= limit) seed = seed >> 3;
+
+ this.current = 4*seed+1;
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/DRand.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,267 @@
+package org.apache.mahout.jet.random.engine;
+
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+import java.util.Date;
+/**
+MersenneTwister (MT19937) is one of the strongest uniform pseudo-random number generators known so far; at the same time it is quick.
+Produces uniformly distributed <tt>int</tt>'s and <tt>long</tt>'s in the closed intervals <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt> and <tt>[Long.MIN_VALUE,Long.MAX_VALUE]</tt>, respectively,
+as well as <tt>float</tt>'s and <tt>double</tt>'s in the open unit intervals <tt>(0.0f,1.0f)</tt> and <tt>(0.0,1.0)</tt>, respectively.
+The seed can be any 32-bit integer except <tt>0</tt>. Shawn J. Cokus commented that perhaps the seed should preferably be odd.
+<p>
+<b>Quality:</b> MersenneTwister is designed to pass the k-distribution test. It has an astronomically large period of 2<sup>19937</sup>-1 (=10<sup>6001</sup>) and 623-dimensional equidistribution up to 32-bit accuracy.
+It passes many stringent statistical tests, including the <A HREF="http://stat.fsu.edu/~geo/diehard.html">diehard</A> test of G. Marsaglia and the load test of P. Hellekalek and S. Wegenkittl.
+<p>
+<b>Performance:</b> Its speed is comparable to other modern generators (in particular, as fast as <tt>java.util.Random.nextFloat()</tt>).
+2.5 million calls to <tt>raw()</tt> per second (Pentium Pro 200 Mhz, JDK 1.2, NT).
+Be aware, however, that there is a non-negligible amount of overhead required to initialize the data
+structures used by a MersenneTwister. Code like
+<pre>
+ double sum = 0.0;
+ for (int i=0; i<100000; ++i) {
+ RandomElement twister = new MersenneTwister(new java.util.Date());
+ sum += twister.raw();
+ }
+</pre>
+will be wildly inefficient. Consider using
+<pre>
+ double sum = 0.0;
+ RandomElement twister = new MersenneTwister(new java.util.Date());
+ for (int i=0; i<100000; ++i) {
+ sum += twister.raw();
+ }
+</pre>
+instead. This allows the cost of constructing the MersenneTwister object
+to be borne only once, rather than once for each iteration in the loop.
+<p>
+<b>Implementation:</b> After M. Matsumoto and T. Nishimura,
+"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.
+<dt>More info on <A HREF="http://www.math.keio.ac.jp/~matumoto/eindex.html"> Masumoto's homepage</A>.
+<dt>More info on <A HREF="http://www.ncsa.uiuc.edu/Apps/CMP/RNG/www-rng.html"> Pseudo-random number generators is on the Web</A>.
+<dt>Yet <A HREF="http://nhse.npac.syr.edu/random"> some more info</A>.
+<p>
+The correctness of this implementation has been verified against the published output sequence
+<a href="http://www.math.keio.ac.jp/~nisimura/random/real2/mt19937-2.out">mt19937-2.out</a> of the C-implementation
+<a href="http://www.math.keio.ac.jp/~nisimura/random/real2/mt19937-2.c">mt19937-2.c</a>.
+(Call <tt>test(1000)</tt> to print the sequence).
+<dt>
+Note that this implementation is <b>not synchronized</b>.
+<p>
+<b>Details:</b> MersenneTwister is designed with consideration of the flaws of various existing generators in mind.
+It is an improved version of TT800, a very successful generator.
+MersenneTwister is based on linear recurrences modulo 2.
+Such generators are very fast, have extremely long periods, and appear quite robust.
+MersenneTwister produces 32-bit numbers, and every <tt>k</tt>-dimensional vector of such numbers appears the same number of times as <tt>k</tt> successive values over the
+period length, for each <tt>k <= 623</tt> (except for the zero vector, which appears one time less).
+If one looks at only the first <tt>n <= 16</tt> bits of each number, then the property holds for even larger <tt>k</tt>, as shown in the following table (taken from the publication cited above):
+<div align="center">
+<table width="75%" border="1" cellspacing="0" cellpadding="0">
+ <tr>
+ <td width="2%"> <div align="center">n</div> </td>
+ <td width="6%"> <div align="center">1</div> </td>
+ <td width="5%"> <div align="center">2</div> </td>
+ <td width="5%"> <div align="center">3</div> </td>
+ <td width="5%"> <div align="center">4</div> </td>
+ <td width="5%"> <div align="center">5</div> </td>
+ <td width="5%"> <div align="center">6</div> </td>
+ <td width="5%"> <div align="center">7</div> </td>
+ <td width="5%"> <div align="center">8</div> </td>
+ <td width="5%"> <div align="center">9</div> </td>
+ <td width="5%"> <div align="center">10</div> </td>
+ <td width="5%"> <div align="center">11</div> </td>
+ <td width="10%"> <div align="center">12 .. 16</div> </td>
+ <td width="10%"> <div align="center">17 .. 32</div> </td>
+ </tr>
+ <tr>
+ <td width="2%"> <div align="center">k</div> </td>
+ <td width="6%"> <div align="center">19937</div> </td>
+ <td width="5%"> <div align="center">9968</div> </td>
+ <td width="5%"> <div align="center">6240</div> </td>
+ <td width="5%"> <div align="center">4984</div> </td>
+ <td width="5%"> <div align="center">3738</div> </td>
+ <td width="5%"> <div align="center">3115</div> </td>
+ <td width="5%"> <div align="center">2493</div> </td>
+ <td width="5%"> <div align="center">2492</div> </td>
+ <td width="5%"> <div align="center">1869</div> </td>
+ <td width="5%"> <div align="center">1869</div> </td>
+ <td width="5%"> <div align="center">1248</div> </td>
+ <td width="10%"> <div align="center">1246</div> </td>
+ <td width="10%"> <div align="center">623</div> </td>
+ </tr>
+</table>
+</div>
+<p>
+MersenneTwister generates random numbers in batches of 624 numbers at a time, so the caching and pipelining of modern systems is exploited.
+The generator is implemented to generate the output by using the fastest arithmetic operations only: 32-bit additions and bit operations (no division, no multiplication, no mod).
+These operations generate sequences of 32 random bits (<tt>int</tt>'s).
+<tt>long</tt>'s are formed by concatenating two 32 bit <tt>int</tt>'s.
+<tt>float</tt>'s are formed by dividing the interval <tt>[0.0,1.0]</tt> into 2<sup>32</sup> sub intervals, then randomly choosing one subinterval.
+<tt>double</tt>'s are formed by dividing the interval <tt>[0.0,1.0]</tt> into 2<sup>64</sup> sub intervals, then randomly choosing one subinterval.
+<p>
+@author wolfgang.hoschek@cern.ch
+@version 1.0, 09/24/99
+@see java.util.Random
+*/
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class MersenneTwister extends RandomEngine {
+ private int mti;
+ private int[] mt = new int[N]; /* set initial seeds: N = 624 words */
+
+ /* Period parameters */
+ private static final int N=624;
+ private static final int M=397;
+ private static final int MATRIX_A=0x9908b0df; /* constant vector a */
+ private static final int UPPER_MASK=0x80000000; /* most significant w-r bits */
+ private static final int LOWER_MASK=0x7fffffff; /* least significant r bits */
+
+ /* for tempering */
+ private static final int TEMPERING_MASK_B=0x9d2c5680;
+ private static final int TEMPERING_MASK_C=0xefc60000;
+
+ private static final int mag0 = 0x0;
+ private static final int mag1 = MATRIX_A;
+ //private static final int[] mag01=new int[] {0x0, MATRIX_A};
+ /* mag01[x] = x * MATRIX_A for x=0,1 */
+
+ public static final int DEFAULT_SEED = 4357;
+/**
+ * Constructs and returns a random number generator with a default seed, which is a <b>constant</b>.
+ * Thus using this constructor will yield generators that always produce exactly the same sequence.
+ * This method is mainly intended to ease testing and debugging.
+ */
+public MersenneTwister() {
+ this(DEFAULT_SEED);
+}
+/**
+ * Constructs and returns a random number generator with the given seed.
+ */
+public MersenneTwister(int seed) {
+ setSeed(seed);
+}
+/**
+ * Constructs and returns a random number generator seeded with the given date.
+ *
+ * @param d typically <tt>new java.util.Date()</tt>
+ */
+public MersenneTwister(Date d) {
+ this((int)d.getTime());
+}
+/**
+ * Returns a copy of the receiver; the copy will produce identical sequences.
+ * After this call has returned, the copy and the receiver have equal but separate state.
+ *
+ * @return a copy of the receiver.
+ */
+public Object clone() {
+ MersenneTwister clone = (MersenneTwister) super.clone();
+ clone.mt = (int[]) this.mt.clone();
+ return clone;
+}
+/**
+ * Generates N words at one time.
+ */
+protected void nextBlock() {
+ /*
+ // ******************** OPTIMIZED **********************
+ // only 5-10% faster ?
+ int y;
+
+ int kk;
+ int[] cache = mt; // cached for speed
+ int kkM;
+ int limit = N-M;
+ for (kk=0,kkM=kk+M; kk<limit; kk++,kkM++) {
+ y = (cache[kk]&UPPER_MASK)|(cache[kk+1]&LOWER_MASK);
+ cache[kk] = cache[kkM] ^ (y >>> 1) ^ ((y & 0x1) == 0 ? mag0 : mag1);
+ }
+ limit = N-1;
+ for (kkM=kk+(M-N); kk<limit; kk++,kkM++) {
+ y = (cache[kk]&UPPER_MASK)|(cache[kk+1]&LOWER_MASK);
+ cache[kk] = cache[kkM] ^ (y >>> 1) ^ ((y & 0x1) == 0 ? mag0 : mag1);
+ }
+ y = (cache[N-1]&UPPER_MASK)|(cache[0]&LOWER_MASK);
+ cache[N-1] = cache[M-1] ^ (y >>> 1) ^ ((y & 0x1) == 0 ? mag0 : mag1);
+
+ this.mt = cache;
+ this.mti = 0;
+ */
+
+
+
+ // ******************** UNOPTIMIZED **********************
+ int y;
+
+ int kk;
+
+ for (kk=0;kk<N-M;kk++) {
+ y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+ mt[kk] = mt[kk+M] ^ (y >>> 1) ^ ((y & 0x1) == 0 ? mag0 : mag1);
+ }
+ for (;kk<N-1;kk++) {
+ y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+ mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ ((y & 0x1) == 0 ? mag0 : mag1);
+ }
+ y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
+ mt[N-1] = mt[M-1] ^ (y >>> 1) ^ ((y & 0x1) == 0 ? mag0 : mag1);
+
+ this.mti = 0;
+
+}
+/**
+ * Returns a 32 bit uniformly distributed random number in the closed interval <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt> (including <tt>Integer.MIN_VALUE</tt> and <tt>Integer.MAX_VALUE</tt>).
+ */
+public int nextInt() {
+ /* Each single bit including the sign bit will be random */
+ if (mti == N) nextBlock(); // generate N ints at one time
+
+ int y = mt[mti++];
+ y ^= y >>> 11; // y ^= TEMPERING_SHIFT_U(y );
+ y ^= (y << 7) & TEMPERING_MASK_B; // y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
+ y ^= (y << 15) & TEMPERING_MASK_C; // y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
+ // y &= 0xffffffff; //you may delete this line if word size = 32
+ y ^= y >>> 18; // y ^= TEMPERING_SHIFT_L(y);
+
+ return y;
+}
+/**
+ * Sets the receiver's seed.
+ * This method resets the receiver's entire internal state.
+ */
+protected void setSeed(int seed) {
+ mt[0] = seed & 0xffffffff;
+ for (int i = 1; i < N; i++) {
+ mt[i] = (1812433253 * (mt[i-1] ^ (mt[i-1] >> 30)) + i);
+ /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
+ /* In the previous versions, MSBs of the seed affect */
+ /* only MSBs of the array mt[]. */
+ /* 2002/01/09 modified by Makoto Matsumoto */
+ mt[i] &= 0xffffffff;
+ /* for >32 bit machines */
+ }
+ //System.out.println("init done");
+ mti = N;
+
+ /*
+ old version was:
+ for (int i = 0; i < N; i++) {
+ mt[i] = seed & 0xffff0000;
+ seed = 69069 * seed + 1;
+ mt[i] |= (seed & 0xffff0000) >>> 16;
+ seed = 69069 * seed + 1;
+ }
+ //System.out.println("init done");
+ mti = N;
+ */
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister.java
------------------------------------------------------------------------------
svn:eol-style = native
Added: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister64.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister64.java?rev=883365&view=auto
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister64.java (added)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister64.java Mon Nov 23 15:14:26 2009
@@ -0,0 +1,51 @@
+/*
+Copyright � 1999 CERN - European Organization for Nuclear Research.
+Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
+is hereby granted without fee, provided that the above copyright notice appear in all copies and
+that both that copyright notice and this permission notice appear in supporting documentation.
+CERN makes no representations about the suitability of this software for any purpose.
+It is provided "as is" without expressed or implied warranty.
+*/
+package org.apache.mahout.jet.random.engine;
+
+import java.util.Date;
+/**
+ * Same as <tt>MersenneTwister</tt> except that method <tt>raw()</tt> returns 64 bit random numbers instead of 32 bit random numbers.
+ *
+ * @author wolfgang.hoschek@cern.ch
+ * @version 1.0, 09/24/99
+ * @see MersenneTwister
+ */
+/**
+ * @deprecated until unit tests are in place. Until this time, this class/interface is unsupported.
+ */
+@Deprecated
+public class MersenneTwister64 extends MersenneTwister {
+/**
+ * Constructs and returns a random number generator with a default seed, which is a <b>constant</b>.
+ */
+public MersenneTwister64() {
+ super();
+}
+/**
+ * Constructs and returns a random number generator with the given seed.
+ * @param seed should not be 0, in such a case <tt>MersenneTwister64.DEFAULT_SEED</tt> is silently substituted.
+ */
+public MersenneTwister64(int seed) {
+ super(seed);
+}
+/**
+ * Constructs and returns a random number generator seeded with the given date.
+ *
+ * @param d typically <tt>new java.util.Date()</tt>
+ */
+public MersenneTwister64(Date d) {
+ super(d);
+}
+/**
+ * Returns a 64 bit uniformly distributed random number in the open unit interval <code>(0.0,1.0)</code> (excluding 0.0 and 1.0).
+ */
+public double raw() {
+ return nextDouble();
+}
+}
Propchange: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister64.java
------------------------------------------------------------------------------
svn:eol-style = native