You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mahout.apache.org by sr...@apache.org on 2009/11/25 04:31:49 UTC
svn commit: r883972 [5/10] - in /lucene/mahout/trunk:
core/src/main/java/org/apache/mahout/fpm/pfpgrowth/
matrix/src/main/java/org/apache/mahout/jet/math/
matrix/src/main/java/org/apache/mahout/jet/random/
matrix/src/main/java/org/apache/mahout/jet/ran...
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/NegativeBinomial.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/NegativeBinomial.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/NegativeBinomial.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/NegativeBinomial.java Wed Nov 25 03:31:47 2009
@@ -24,22 +24,20 @@
* <p>
* J.H. Ahrens, U. Dieter (1974): Computer methods for sampling from gamma, beta, Poisson and binomial distributions, Computing 12, 223--246.
*
- * @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 NegativeBinomial extends AbstractDiscreteDistribution {
- protected int n;
- protected double p;
+ protected int n;
+ protected double p;
- protected Gamma gamma;
- protected Poisson poisson;
-
- // The uniform random number generated shared by all <b>static</b> methods.
- protected static NegativeBinomial shared = new NegativeBinomial(1,0.5,makeDefaultGenerator());
+ protected Gamma gamma;
+ protected Poisson poisson;
+
+ // The uniform random number generated shared by all <b>static</b> methods.
+ protected static NegativeBinomial shared = new NegativeBinomial(1,0.5,makeDefaultGenerator());
/**
* Constructs a Negative Binomial distribution.
* Example: n=1, p=0.5.
@@ -48,16 +46,16 @@
* @param randomGenerator a uniform random number generator.
*/
public NegativeBinomial(int n, double p, RandomEngine randomGenerator) {
- setRandomGenerator(randomGenerator);
- setNandP(n,p);
- this.gamma = new Gamma(n,1.0,randomGenerator);
- this.poisson = new Poisson(0.0,randomGenerator);
+ setRandomGenerator(randomGenerator);
+ setNandP(n,p);
+ this.gamma = new Gamma(n,1.0,randomGenerator);
+ this.poisson = new Poisson(0.0,randomGenerator);
}
/**
* Returns the cumulative distribution function.
*/
public double cdf(int k) {
- return Probability.negativeBinomial(k,n,p);
+ return Probability.negativeBinomial(k,n,p);
}
/**
* Returns a deep copy of the receiver; the copy will produce identical sequences.
@@ -66,18 +64,18 @@
* @return a copy of the receiver.
*/
public Object clone() {
- NegativeBinomial copy = (NegativeBinomial) super.clone();
- if (this.poisson != null) copy.poisson = (Poisson) this.poisson.clone();
- copy.poisson.setRandomGenerator(copy.getRandomGenerator());
- if (this.gamma != null) copy.gamma = (Gamma) this.gamma.clone();
- copy.gamma.setRandomGenerator(copy.getRandomGenerator());
- return copy;
+ NegativeBinomial copy = (NegativeBinomial) super.clone();
+ if (this.poisson != null) copy.poisson = (Poisson) this.poisson.clone();
+ copy.poisson.setRandomGenerator(copy.getRandomGenerator());
+ if (this.gamma != null) copy.gamma = (Gamma) this.gamma.clone();
+ copy.gamma.setRandomGenerator(copy.getRandomGenerator());
+ return copy;
}
/**
* Returns a random number from the distribution.
*/
public int nextInt() {
- return nextInt(n,p);
+ return nextInt(n,p);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
@@ -107,17 +105,17 @@
* *
******************************************************************/
- double x = p /(1.0 - p);
- double p1 = p;
- double y = x * this.gamma.nextDouble(n,1.0);
- return this.poisson.nextInt(y);
+ double x = p /(1.0 - p);
+ double p1 = p;
+ double y = x * this.gamma.nextDouble(n,1.0);
+ return this.poisson.nextInt(y);
}
/**
* Returns the probability distribution function.
*/
public double pdf(int k) {
- if (k > n) throw new IllegalArgumentException();
- return org.apache.mahout.jet.math.Arithmetic.binomial(n,k) * Math.pow(p,k) * Math.pow(1.0-p,n-k);
+ if (k > n) throw new IllegalArgumentException();
+ return org.apache.mahout.jet.math.Arithmetic.binomial(n,k) * Math.pow(p,k) * Math.pow(1.0-p,n-k);
}
/**
* Sets the parameters number of trials and the probability of success.
@@ -125,8 +123,8 @@
* @param p the probability of success.
*/
public void setNandP(int n, double p) {
- this.n = n;
- this.p = p;
+ this.n = n;
+ this.p = p;
}
/**
* Returns a random number from the distribution with the given parameters n and p.
@@ -134,23 +132,23 @@
* @param p the probability of success.
*/
public static int staticNextInt(int n, double p) {
- synchronized (shared) {
- return shared.nextInt(n,p);
- }
+ synchronized (shared) {
+ return shared.nextInt(n,p);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+n+","+p+")";
+ return this.getClass().getName()+"("+n+","+p+")";
}
/**
* 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);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Normal.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Normal.java Wed Nov 25 03:31:47 2009
@@ -14,17 +14,17 @@
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.
+ 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>
@@ -43,107 +43,107 @@
*/
@Deprecated
public class Normal extends AbstractContinousDistribution {
- protected double mean;
- protected double variance;
- protected double standardDeviation;
+ protected double mean;
+ protected double variance;
+ protected double standardDeviation;
- protected double cache; // cache for Box-Mueller algorithm
- protected boolean cacheFilled; // Box-Mueller
+ protected double cache; // cache for Box-Mueller algorithm
+ protected boolean cacheFilled; // Box-Mueller
- protected double SQRT_INV; // performance cache
+ 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());
+ // 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);
+ setRandomGenerator(randomGenerator);
+ setState(mean,standardDeviation);
}
/**
* Returns the cumulative distribution function.
*/
public double cdf(double x) {
- return Probability.normal(mean,variance,x);
+ return Probability.normal(mean,variance,x);
}
/**
* Returns a random number from the distribution.
*/
public double nextDouble() {
- return nextDouble(this.mean,this.standardDeviation);
+ 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;
+ // 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));
+ 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;
+ 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);
- }
+ 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);
- }
+ synchronized (shared) {
+ return shared.nextDouble(mean,standardDeviation);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+mean+","+standardDeviation+")";
+ 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);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Poisson.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Poisson.java Wed Nov 25 03:31:47 2009
@@ -37,53 +37,51 @@
* 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;
+ 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
-
+ // 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());
+ // 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);
+ setRandomGenerator(randomGenerator);
+ setMean(mean);
}
/**
* Returns the cumulative distribution function.
*/
public double cdf(int k) {
- return Probability.poisson(k,this.mean);
+ return Probability.poisson(k,this.mean);
}
/**
* Returns a deep copy of the receiver; the copy will produce identical sequences.
@@ -92,18 +90,18 @@
* @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;
+ 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);
+ 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);
+ return nextInt(this.mean);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
@@ -125,215 +123,215 @@
* 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;
+ 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;
- }
+ 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);
+ 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;
+ 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();
- }
+ synchronized (shared) {
+ shared.setMean(mean);
+ return shared.nextInt();
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+mean+")";
+ 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);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/PoissonSlow.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/PoissonSlow.java Wed Nov 25 03:31:47 2009
@@ -26,38 +26,36 @@
* 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;
+ 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};
+ // 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());
+ // 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);
+ setRandomGenerator(randomGenerator);
+ setMean(mean);
}
/**
* Returns the value ln(Gamma(xx) for xx > 0. Full accuracy is obtained for
@@ -65,119 +63,119 @@
* (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);
+ 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);
+ 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;
- }
+ /*
+ * 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;
+ 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);
- }
- }
+ 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();
- }
+ synchronized (shared) {
+ shared.setMean(mean);
+ return shared.nextInt();
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+mean+")";
+ 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);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Stack.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Stack.java Wed Nov 25 03:31:47 2009
@@ -12,76 +12,76 @@
* 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 */
+ 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];
+ 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;
+ 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;
+ 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]);
+ 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;
+ return s->i + 1;
}
static void free_stack(stack_t *s)
{
- free((char *)(s->v));
- free((char *)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];
+ 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;
+ 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;
+ return i+1;
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/StudentT.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/StudentT.java Wed Nov 25 03:31:47 2009
@@ -31,19 +31,17 @@
* 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 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());
+ 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.
@@ -51,20 +49,20 @@
* @throws IllegalArgumentException if <tt>freedom <= 0.0</tt>.
*/
public StudentT(double freedom, RandomEngine randomGenerator) {
- setRandomGenerator(randomGenerator);
- setState(freedom);
+ setRandomGenerator(randomGenerator);
+ setState(freedom);
}
/**
* Returns the cumulative distribution function.
*/
public double cdf(double x) {
- return Probability.studentT(freedom,x);
+ return Probability.studentT(freedom,x);
}
/**
* Returns a random number from the distribution.
*/
public double nextDouble() {
- return nextDouble(this.freedom);
+ return nextDouble(this.freedom);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
@@ -72,32 +70,32 @@
* @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;
+ /*
+ * 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);
+ 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));
+ 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);
+ return this.TERM * Math.pow((1+ x*x/freedom), -(freedom+1) * 0.5);
}
/**
* Sets the distribution parameter.
@@ -105,11 +103,11 @@
* @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);
+ 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.
@@ -117,23 +115,23 @@
* @throws IllegalArgumentException if <tt>freedom <= 0.0</tt>.
*/
public static double staticNextDouble(double freedom) {
- synchronized (shared) {
- return shared.nextDouble(freedom);
- }
+ synchronized (shared) {
+ return shared.nextDouble(freedom);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+freedom+")";
+ 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);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Uniform.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Uniform.java Wed Nov 25 03:31:47 2009
@@ -17,210 +17,208 @@
* <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());
+ 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));
+ 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);
+ 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);
+ 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);
+ 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;
+ 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();
+ 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();
+ 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);
+ 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));
+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()));
+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));
- }
+ /* 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));
- }
+ // 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;
+ // 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);
+ 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;
+ 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();
- }
+ 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();
- }
+ 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);
- }
+ 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);
- }
+ 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);
- }
+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);
- }
+ 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);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+min+","+max+")";
+ return this.getClass().getName()+"("+min+","+max+")";
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/VonMises.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/VonMises.java Wed Nov 25 03:31:47 2009
@@ -27,36 +27,34 @@
* <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;
+ protected double my_k;
- // cached vars for method nextDouble(a) (for performance only)
- private double k_set = -1.0;
- private double tau,rho,r;
+ // 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());
+ // 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);
+ setRandomGenerator(randomGenerator);
+ setState(freedom);
}
/**
* Returns a random number from the distribution.
*/
public double nextDouble() {
- return nextDouble(this.my_k);
+ return nextDouble(this.my_k);
}
/**
* Returns a random number from the distribution; bypasses the internal state.
@@ -81,59 +79,59 @@
* *
* Implemented by F. Niederl, August 1992 *
******************************************************************/
- double u,v,w,c,z;
+ 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;
- }
+ 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 //
+ // 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;
+ 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);
- }
+ synchronized (shared) {
+ return shared.nextDouble(freedom);
+ }
}
/**
* Returns a String representation of the receiver.
*/
public String toString() {
- return this.getClass().getName()+"("+my_k+")";
+ 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);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Zeta.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Zeta.java Wed Nov 25 03:31:47 2009
@@ -30,29 +30,27 @@
* <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;
+ 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;
+ // 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());
+ // 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);
+ setRandomGenerator(randomGenerator);
+ setState(ro,pk);
}
/**
* Returns a zeta distributed random number.
@@ -96,69 +94,69 @@
* Variate Generation, Clarendon Press, Oxford. *
* *
******************************************************************/
- double u,v,e,x;
- long k;
+ 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;
+ 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);
+ return (int) generateZeta(ro, pk, randomGenerator);
}
/**
* Sets the parameters.
*/
public void setState(double ro, double pk) {
- this.ro = ro;
- this.pk = 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();
- }
+ 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+")";
+ 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);
- }
+ synchronized (shared) {
+ shared.setRandomGenerator(randomGenerator);
+ }
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/Benchmark.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/Benchmark.java Wed Nov 25 03:31:47 2009
@@ -75,8 +75,6 @@
*
*
* @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.
@@ -87,167 +85,167 @@
* Makes this class non instantiable, but still let's others inherit from it.
*/
protected Benchmark() {
- throw new RuntimeException("Non instantiable");
+ throw new RuntimeException("Non instantiable");
}
/**
* Benchmarks <tt>raw()</tt> for various uniform generation engines.
*/
public static void benchmark(int times) {
- org.apache.mahout.matrix.Timer timer = new org.apache.mahout.matrix.Timer();
- RandomEngine gen;
+ org.apache.mahout.matrix.Timer timer = new org.apache.mahout.matrix.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.");
- */
+ 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");
-
+ 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);
- }
+ 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;
+ RandomEngine random;
- /*
- System.out.println("raw():");
- random = (RandomEngine) randomEngine.clone();
- //org.apache.mahout.matrix.Timer timer = new org.apache.mahout.matrix.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();
- }
+ /*
+ System.out.println("raw():");
+ random = (RandomEngine) randomEngine.clone();
+ //org.apache.mahout.matrix.Timer timer = new org.apache.mahout.matrix.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");
+ //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.matrix.set.OpenMultiFloatHashSet multiset = new org.apache.mahout.matrix.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.matrix.Timer timer = new org.apache.mahout.matrix.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");
+ System.out.println("from="+from+", to="+to);
+
+ //org.apache.mahout.matrix.set.OpenMultiFloatHashSet multiset = new org.apache.mahout.matrix.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.matrix.Timer timer = new org.apache.mahout.matrix.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");
}
}
Modified: 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=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/DRand.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/DRand.java Wed Nov 25 03:31:47 2009
@@ -33,8 +33,6 @@
* 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
*/
@@ -43,20 +41,20 @@
*/
@Deprecated
public class DRand extends RandomEngine {
- private int current;
- public static final int DEFAULT_SEED = 1;
+ 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);
+ 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);
+ setSeed(seed);
}
/**
* Constructs and returns a random number generator seeded with the given date.
@@ -64,16 +62,16 @@
* @param d typically <tt>new java.util.Date()</tt>
*/
public DRand(Date d) {
- this((int)d.getTime());
+ 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;
+ current *= 0x278DDE6D; /* z(i+1)=a*z(i) (mod 2**32) */
+ // a == 0x278DDE6D == 663608941
+
+ return current;
}
/**
* Sets the receiver's seed.
@@ -82,10 +80,10 @@
* @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;
+ 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;
+ this.current = 4*seed+1;
}
}