cern.jet.math.Arithmetic Java Examples

/**
 * Constructs an approximate quantile finder with b buffers, each having k elements.
 * @param b the number of buffers
 * @param k the number of elements per buffer
 * @param N the total number of elements over which quantiles are to be computed.
 * @param samplingRate 1.0 --> all elements are consumed. 10.0 --> Consumes one random element from successive blocks of 10 elements each. Etc.
 * @param generator a uniform random number generator.
 */
public KnownDoubleQuantileEstimator(int b, int k, long N, double samplingRate, RandomEngine generator) {
	this.samplingRate = samplingRate;
	this.N = N;

	if (this.samplingRate <= 1.0) {
		this.samplingAssistant = null;
	}
	else {
		this.samplingAssistant = new RandomSamplingAssistant(Arithmetic.floor(N/samplingRate), N, generator);
	}
	
	setUp(b,k);
	this.clear();
}
 

/**
 * 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 parameters number of trials and the probability of success.
 * @param n the number of trials
 * @param p the probability of success.
 * @throws IllegalArgumentException if <tt>n*Math.min(p,1-p) &lt;= 0.0</tt>
 */
public void setNandP(int n, double p) {
	if (n*Math.min(p,1-p) <= 0.0) throw new IllegalArgumentException();
	this.n = n;
	this.p = p;
	
	this.log_p = Math.log(p);
	this.log_q = Math.log(1.0-p);
	this.log_n = Arithmetic.logFactorial(n);
}
 

/**
 * Removes all elements from the receiver.  The receiver will
 * be empty after this call returns, and its memory requirements will be close to zero.
 */
public void clear() {
	super.clear();
	this.beta=1.0;
	this.weHadMoreThanOneEmptyBuffer = false;
	//this.setSamplingRate(samplingRate,N);

	RandomSamplingAssistant assist = this.samplingAssistant;
	if (assist != null) {
		this.samplingAssistant = new RandomSamplingAssistant(Arithmetic.floor(N/samplingRate), N, assist.getRandomGenerator());
	}
}
 

/**
 * Constructs an approximate quantile finder with b buffers, each having k elements.
 * @param b the number of buffers
 * @param k the number of elements per buffer
 * @param N the total number of elements over which quantiles are to be computed.
 * @param samplingRate 1.0 --> all elements are consumed. 10.0 --> Consumes one random element from successive blocks of 10 elements each. Etc.
 * @param generator a uniform random number generator.
 */
public KnownDoubleQuantileEstimator(int b, int k, long N, double samplingRate, RandomEngine generator) {
	this.samplingRate = samplingRate;
	this.N = N;

	if (this.samplingRate <= 1.0) {
		this.samplingAssistant = null;
	}
	else {
		this.samplingAssistant = new RandomSamplingAssistant(Arithmetic.floor(N/samplingRate), N, generator);
	}
	
	setUp(b,k);
	this.clear();
}
 

/**
 * 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 parameters number of trials and the probability of success.
 * @param n the number of trials
 * @param p the probability of success.
 * @throws IllegalArgumentException if <tt>n*Math.min(p,1-p) &lt;= 0.0</tt>
 */
public void setNandP(int n, double p) {
	if (n*Math.min(p,1-p) <= 0.0) throw new IllegalArgumentException();
	this.n = n;
	this.p = p;
	
	this.log_p = Math.log(p);
	this.log_q = Math.log(1.0-p);
	this.log_n = Arithmetic.logFactorial(n);
}
 

/**
 * Removes all elements from the receiver.  The receiver will
 * be empty after this call returns, and its memory requirements will be close to zero.
 */
public void clear() {
	super.clear();
	this.beta=1.0;
	this.weHadMoreThanOneEmptyBuffer = false;
	//this.setSamplingRate(samplingRate,N);

	RandomSamplingAssistant assist = this.samplingAssistant;
	if (assist != null) {
		this.samplingAssistant = new RandomSamplingAssistant(Arithmetic.floor(N/samplingRate), N, assist.getRandomGenerator());
	}
}
 

/**
 * Returns the probability distribution function.
 */
public double pdf(int k) {
	if (k < 0) throw new IllegalArgumentException();
	int r = this.n - k;
	return Math.exp(this.log_n - Arithmetic.logFactorial(k) - Arithmetic.logFactorial(r) + this.log_p * k + this.log_q * r);
}
 

private static double f(int k, double l_nu, double c_pm) {
	return  Math.exp(k * l_nu - Arithmetic.logFactorial(k) - c_pm);
}
 

/**
 * Returns the probability distribution function.
 */
public double pdf(int k) {
	return Arithmetic.binomial(my_s, k) * Arithmetic.binomial(my_N - my_s, my_n - k) 
		/ Arithmetic.binomial(my_N, my_n);
}
 

/**
 * Computes the number of buffers and number of values per buffer such that
 * quantiles can be determined with an approximation error no more than epsilon with a certain probability.
 * Assumes that quantiles are to be computed over N values.
 * The required sampling rate is computed and stored in the first element of the provided <tt>returnSamplingRate</tt> array, which, therefore must be at least of length 1.
 * @param N the anticipated number of values over which quantiles shall be computed (e.g 10^6).
 * @param epsilon the approximation error which is guaranteed not to be exceeded (e.g. <tt>0.001</tt>) (<tt>0 &lt;= epsilon &lt;= 1</tt>). To get exact result, set <tt>epsilon=0.0</tt>;
 * @param delta the probability that the approximation error is more than than epsilon (e.g. <tt>0.0001</tt>) (<tt>0 &lt;= delta &lt;= 1</tt>). To avoid probabilistic answers, set <tt>delta=0.0</tt>.
 * @param quantiles the number of quantiles to be computed (e.g. <tt>100</tt>) (<tt>quantiles &gt;= 1</tt>). If unknown in advance, set this number large, e.g. <tt>quantiles &gt;= 10000</tt>.
 * @param samplingRate a <tt>double[1]</tt> where the sampling rate is to be filled in.
 * @return <tt>long[2]</tt> - <tt>long[0]</tt>=the number of buffers, <tt>long[1]</tt>=the number of elements per buffer, <tt>returnSamplingRate[0]</tt>=the required sampling rate.
 */
protected static long[] known_N_compute_B_and_K_slow(long N, double epsilon, double delta, int quantiles, double[] returnSamplingRate) {
	final int maxBuffers = 50;
	final int maxHeight = 50;
	final double N_double = N;

	// One possibility is to use one buffer of size N
	//
	long ret_b = 1;
	long ret_k = N;
	double sampling_rate = 1.0;
	long memory = N;


	// Otherwise, there are at least two buffers (b >= 2)
	// and the height of the tree is at least three (h >= 3)
	//
	// We restrict the search for b and h to MAX_BINOM, a large enough value for
	// practical values of    epsilon >= 0.001   and    delta >= 0.00001
	//
	final double logarithm = Math.log(2.0*quantiles/delta);
	final double c = 2.0 * epsilon * N_double;
	for (long b=2 ; b<maxBuffers ; b++)
		for (long h=3 ; h<maxHeight ; h++) {
			double binomial = Arithmetic.binomial(b+h-2, h-1);
			long tmp = (long) Math.ceil(N_double / binomial);
			if ((b * tmp < memory) && 
					((h-2) * binomial - Arithmetic.binomial(b+h-3, h-3) + Arithmetic.binomial(b+h-3, h-2)
					<= c) ) {
				ret_k = tmp ;
				ret_b = b ;
				memory = ret_k * b;
				sampling_rate = 1.0 ;
			}
			if (delta > 0.0) {
				double t = (h-2) * Arithmetic.binomial(b+h-2, h-1) - Arithmetic.binomial(b+h-3, h-3) + Arithmetic.binomial(b+h-3, h-2) ;
				double u = logarithm / epsilon ;
				double v = Arithmetic.binomial (b+h-2, h-1) ;
				double w = logarithm / (2.0*epsilon*epsilon) ;

				// From our SIGMOD 98 paper, we have two equantions to satisfy:
				// t  <= u * alpha/(1-alpha)^2
				// kv >= w/(1-alpha)^2
				//
				// Denoting 1/(1-alpha)    by x,
				// we see that the first inequality is equivalent to
				// t/u <= x^2 - x
				// which is satisfied by x >= 0.5 + 0.5 * sqrt (1 + 4t/u)
				// Plugging in this value into second equation yields
				// k >= wx^2/v

				double x = 0.5 + 0.5 * Math.sqrt(1.0 + 4.0*t/u) ;
				long k = (long) Math.ceil(w*x*x/v) ;
				if (b * k < memory) {
					ret_k = k ;
					ret_b = b ;
					memory = b * k ;
					sampling_rate = N_double*2.0*epsilon*epsilon / logarithm ;
				}
			}
		}
		
	long[] result = new long[2];
	result[0]=ret_b;
	result[1]=ret_k;
	returnSamplingRate[0]=sampling_rate;
	return result;
}
 

/**
 * Returns a random number from the distribution.
 */
protected int hmdu(int N, int M, int n, RandomEngine randomGenerator) {

  int            I, K;
  double              p, nu, c, d, U;

  if (N != N_last || M != M_last || n != n_last) {   // set-up           */
		N_last = N;
	 M_last = M;
	 n_last = n;

	 Mp = (double) (M + 1);
	 np = (double) (n + 1);  N_Mn = N - M - n;

	 p  = Mp / (N + 2.0);
	 nu = np * p;                             /* mode, real       */
	 if ((m = (int) nu) == nu && p == 0.5) {     /* mode, integer    */
		mp = m--;
	 }
	 else {
		mp = m + 1;                           /* mp = m + 1       */
		}

 /* mode probability, using the external function flogfak(k) = ln(k!)    */
	 fm = Math.exp(Arithmetic.logFactorial(N - M) - Arithmetic.logFactorial(N_Mn + m) - Arithmetic.logFactorial(n - m)
		 + Arithmetic.logFactorial(M)     - Arithmetic.logFactorial(M - m)    - Arithmetic.logFactorial(m)
		 - Arithmetic.logFactorial(N)     + Arithmetic.logFactorial(N - n)    + Arithmetic.logFactorial(n)   );

 /* safety bound  -  guarantees at least 17 significant decimal digits   */
 /*                  b = min(n, (long int)(nu + k*c')) */
		b = (int) (nu + 11.0 * Math.sqrt(nu * (1.0 - p) * (1.0 - n/(double)N) + 1.0));	
		if (b > n) b = n;
	}

	for (;;) {
		if ((U = randomGenerator.raw() - fm) <= 0.0)  return(m);
		c = d = fm;

 /* down- and upward search from the mode                                */
		for (I = 1; I <= m; I++) {
			K  = mp - I;                                   /* downward search  */
			c *= (double)K/(np - K) * ((double)(N_Mn + K)/(Mp - K));
			if ((U -= c) <= 0.0)  return(K - 1);

			K  = m + I;                                    /* upward search    */
			d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
			if ((U -= d) <= 0.0)  return(K);
		}

 /* upward search from K = 2m + 1 to K = b                               */
		for (K = mp + m; K <= b; K++) {
			d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
			if ((U -= d) <= 0.0)  return(K);
		}
	}
}
 

private static double fc_lnpk(int k, int N_Mn, int M, int n) {
	return(Arithmetic.logFactorial(k) + Arithmetic.logFactorial(M - k) + Arithmetic.logFactorial(n - k) + Arithmetic.logFactorial(N_Mn + k));
}
 

private static double fc_lnpk(int k, int N_Mn, int M, int n) {
	return(Arithmetic.logFactorial(k) + Arithmetic.logFactorial(M - k) + Arithmetic.logFactorial(n - k) + Arithmetic.logFactorial(N_Mn + k));
}
 

/**
 * Returns a random number from the distribution.
 */
protected int hmdu(int N, int M, int n, RandomEngine randomGenerator) {

  int            I, K;
  double              p, nu, c, d, U;

  if (N != N_last || M != M_last || n != n_last) {   // set-up           */
		N_last = N;
	 M_last = M;
	 n_last = n;

	 Mp = (double) (M + 1);
	 np = (double) (n + 1);  N_Mn = N - M - n;

	 p  = Mp / (N + 2.0);
	 nu = np * p;                             /* mode, real       */
	 if ((m = (int) nu) == nu && p == 0.5) {     /* mode, integer    */
		mp = m--;
	 }
	 else {
		mp = m + 1;                           /* mp = m + 1       */
		}

 /* mode probability, using the external function flogfak(k) = ln(k!)    */
	 fm = Math.exp(Arithmetic.logFactorial(N - M) - Arithmetic.logFactorial(N_Mn + m) - Arithmetic.logFactorial(n - m)
		 + Arithmetic.logFactorial(M)     - Arithmetic.logFactorial(M - m)    - Arithmetic.logFactorial(m)
		 - Arithmetic.logFactorial(N)     + Arithmetic.logFactorial(N - n)    + Arithmetic.logFactorial(n)   );

 /* safety bound  -  guarantees at least 17 significant decimal digits   */
 /*                  b = min(n, (long int)(nu + k*c')) */
		b = (int) (nu + 11.0 * Math.sqrt(nu * (1.0 - p) * (1.0 - n/(double)N) + 1.0));	
		if (b > n) b = n;
	}

	for (;;) {
		if ((U = randomGenerator.raw() - fm) <= 0.0)  return(m);
		c = d = fm;

 /* down- and upward search from the mode                                */
		for (I = 1; I <= m; I++) {
			K  = mp - I;                                   /* downward search  */
			c *= (double)K/(np - K) * ((double)(N_Mn + K)/(Mp - K));
			if ((U -= c) <= 0.0)  return(K - 1);

			K  = m + I;                                    /* upward search    */
			d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
			if ((U -= d) <= 0.0)  return(K);
		}

 /* upward search from K = 2m + 1 to K = b                               */
		for (K = mp + m; K <= b; K++) {
			d *= (np - K)/(double)K * ((Mp - K)/(double)(N_Mn + K));
			if ((U -= d) <= 0.0)  return(K);
		}
	}
}
 

/**
 * Returns the probability distribution function.
 */
public double pdf(int k) {
	return Arithmetic.binomial(my_s, k) * Arithmetic.binomial(my_N - my_s, my_n - k) 
		/ Arithmetic.binomial(my_N, my_n);
}
 

/**
 * Returns the probability distribution function.
 */
public double pdf(int k) {
	if (k < 0) throw new IllegalArgumentException();
	int r = this.n - k;
	return Math.exp(this.log_n - Arithmetic.logFactorial(k) - Arithmetic.logFactorial(r) + this.log_p * k + this.log_q * r);
}
 

private static double f(int k, double l_nu, double c_pm) {
	return  Math.exp(k * l_nu - Arithmetic.logFactorial(k) - c_pm);
}
 

/**
 * Computes the number of buffers and number of values per buffer such that
 * quantiles can be determined with an approximation error no more than epsilon with a certain probability.
 * Assumes that quantiles are to be computed over N values.
 * The required sampling rate is computed and stored in the first element of the provided <tt>returnSamplingRate</tt> array, which, therefore must be at least of length 1.
 * @param N the anticipated number of values over which quantiles shall be computed (e.g 10^6).
 * @param epsilon the approximation error which is guaranteed not to be exceeded (e.g. <tt>0.001</tt>) (<tt>0 &lt;= epsilon &lt;= 1</tt>). To get exact result, set <tt>epsilon=0.0</tt>;
 * @param delta the probability that the approximation error is more than than epsilon (e.g. <tt>0.0001</tt>) (<tt>0 &lt;= delta &lt;= 1</tt>). To avoid probabilistic answers, set <tt>delta=0.0</tt>.
 * @param quantiles the number of quantiles to be computed (e.g. <tt>100</tt>) (<tt>quantiles &gt;= 1</tt>). If unknown in advance, set this number large, e.g. <tt>quantiles &gt;= 10000</tt>.
 * @param samplingRate a <tt>double[1]</tt> where the sampling rate is to be filled in.
 * @return <tt>long[2]</tt> - <tt>long[0]</tt>=the number of buffers, <tt>long[1]</tt>=the number of elements per buffer, <tt>returnSamplingRate[0]</tt>=the required sampling rate.
 */
protected static long[] known_N_compute_B_and_K_slow(long N, double epsilon, double delta, int quantiles, double[] returnSamplingRate) {
	final int maxBuffers = 50;
	final int maxHeight = 50;
	final double N_double = N;

	// One possibility is to use one buffer of size N
	//
	long ret_b = 1;
	long ret_k = N;
	double sampling_rate = 1.0;
	long memory = N;


	// Otherwise, there are at least two buffers (b >= 2)
	// and the height of the tree is at least three (h >= 3)
	//
	// We restrict the search for b and h to MAX_BINOM, a large enough value for
	// practical values of    epsilon >= 0.001   and    delta >= 0.00001
	//
	final double logarithm = Math.log(2.0*quantiles/delta);
	final double c = 2.0 * epsilon * N_double;
	for (long b=2 ; b<maxBuffers ; b++)
		for (long h=3 ; h<maxHeight ; h++) {
			double binomial = Arithmetic.binomial(b+h-2, h-1);
			long tmp = (long) Math.ceil(N_double / binomial);
			if ((b * tmp < memory) && 
					((h-2) * binomial - Arithmetic.binomial(b+h-3, h-3) + Arithmetic.binomial(b+h-3, h-2)
					<= c) ) {
				ret_k = tmp ;
				ret_b = b ;
				memory = ret_k * b;
				sampling_rate = 1.0 ;
			}
			if (delta > 0.0) {
				double t = (h-2) * Arithmetic.binomial(b+h-2, h-1) - Arithmetic.binomial(b+h-3, h-3) + Arithmetic.binomial(b+h-3, h-2) ;
				double u = logarithm / epsilon ;
				double v = Arithmetic.binomial (b+h-2, h-1) ;
				double w = logarithm / (2.0*epsilon*epsilon) ;

				// From our SIGMOD 98 paper, we have two equantions to satisfy:
				// t  <= u * alpha/(1-alpha)^2
				// kv >= w/(1-alpha)^2
				//
				// Denoting 1/(1-alpha)    by x,
				// we see that the first inequality is equivalent to
				// t/u <= x^2 - x
				// which is satisfied by x >= 0.5 + 0.5 * sqrt (1 + 4t/u)
				// Plugging in this value into second equation yields
				// k >= wx^2/v

				double x = 0.5 + 0.5 * Math.sqrt(1.0 + 4.0*t/u) ;
				long k = (long) Math.ceil(w*x*x/v) ;
				if (b * k < memory) {
					ret_k = k ;
					ret_b = b ;
					memory = b * k ;
					sampling_rate = N_double*2.0*epsilon*epsilon / logarithm ;
				}
			}
		}
		
	long[] result = new long[2];
	result[0]=ret_b;
	result[1]=ret_k;
	returnSamplingRate[0]=sampling_rate;
	return result;
}