All Downloads are FREE. Search and download functionalities are using the official Maven repository.

cern.jet.stat.tfloat.quantile.FloatQuantileFinderFactory Maven / Gradle / Ivy

Go to download

Parallel Colt is a multithreaded version of Colt - a library for high performance scientific computing in Java. It contains efficient algorithms for data analysis, linear algebra, multi-dimensional arrays, Fourier transforms, statistics and histogramming.

The newest version!
/*
Copyright (C) 1999 CERN - European Organization for Nuclear Research.
Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose 
is hereby granted without fee, provided that the above copyright notice appear in all copies and 
that both that copyright notice and this permission notice appear in supporting documentation. 
CERN makes no representations about the suitability of this software for any purpose. 
It is provided "as is" without expressed or implied warranty.
 */
package cern.jet.stat.tfloat.quantile;

// import cern.it.util.Floats;
import cern.jet.math.tfloat.FloatArithmetic;
import cern.jet.random.tfloat.engine.FloatRandomEngine;

/**
 * Factory constructing exact and approximate quantile finders for both known
 * and unknown N. Also see {@link hep.aida.tfloat.bin.QuantileFloatBin1D},
 * demonstrating how this package can be used.
 * 
 * The approx. algorithms compute approximate quantiles of large data sequences
 * in a single pass. The approximation guarantees are explicit, and apply for
 * arbitrary value distributions and arrival distributions of the data sequence.
 * The main memory requirements are smaller than for any other known technique
 * by an order of magnitude.
 * 
 * 

* The approx. algorithms are primarily intended to help applications scale. * When faced with a large data sequences, traditional methods either need very * large memories or time consuming disk based sorting. In constrast, the * approx. algorithms can deal with > 10^10 values without disk based sorting. * *

* All classes can be seen from various angles, for example as *

1. Algorithm to compute quantiles. *
2. 1-dim-equi-depth histogram. *
3. 1-dim-histogram arbitrarily rebinnable in real-time. *
4. A space efficient MultiSet data structure using lossy compression. *
5. A space efficient value preserving bin of a 2-dim or d-dim histogram. *
(All subject to an accuracy specified by the user.) * *

* Use methods newXXX(...) to get new instances of one of the * following quantile finders. * *

* 1. Exact quantile finding algorithm for known and unknown N * requiring large main memory. *

* The folkore algorithm: Keeps all elements in main memory, sorts the list, * then picks the quantiles. * * * * *

*

* 2. Approximate quantile finding algorithm for known N requiring * only one pass and little main memory. *

* *

* Needs as input the following parameters: *

*

1. N - the number of values of the data sequence over which * quantiles are to be determined. *
2. quantiles - the number of quantiles to be computed. If * unknown in advance, set this number large, e.g. * quantiles >= 10000. *
3. epsilon - the allowed approximation error on quantiles. * The approximation guarantee of this algorithm is explicit. * *

* It is also possible to couple the approximation algorithm with random * sampling to further reduce memory requirements. With sampling, the * approximation guarantees are explicit but probabilistic, i.e. they apply with * respect to a (user controlled) confidence parameter "delta". * *

4. delta - the probability allowed that the approximation * error fails to be smaller than epsilon. Set delta to zero for * explicit non probabilistic guarantees. * *

* After Gurmeet Singh Manku, Sridhar Rajagopalan and Bruce G. Lindsay, * Approximate Medians and other Quantiles in One Pass and with Limited Memory, * Proc. of the 1998 ACM SIGMOD Int. Conf. on Management of Data, Paper * available here. * * * * *

*

* 3. Approximate quantile finding algorithm for unknown N * requiring only one pass and little main memory. *

* This algorithm requires at most two times the memory of a corresponding * approx. quantile finder knowing N. * *

* Needs as input the following parameters: *

*

2. quantiles - the number of quantiles to be computed. If * unknown in advance, set this number large, e.g. quantiles >= 1000. *
2. epsilon - the allowed approximation error on quantiles. * The approximation guarantee of this algorithm is explicit. * *

* It is also possible to couple the approximation algorithm with random * sampling to further reduce memory requirements. With sampling, the * approximation guarantees are explicit but probabilistic, i.e. they apply with * respect to a (user controlled) confidence parameter "delta". * *

3. delta - the probability allowed that the approximation * error fails to be smaller than epsilon. Set delta to zero for * explicit non probabilistic guarantees. * *

* After Gurmeet Singh Manku, Sridhar Rajagopalan and Bruce G. Lindsay, Random * Sampling Techniques for Space Efficient Online Computation of Order * Statistics of Large Datasets. Proc. of the 1999 ACM SIGMOD Int. Conf. on * Management of Data, Paper available here. * *

* Example usage: * *

 * _TODO_
 * 
* *

* * * @author [email protected] * @version 1.0, 09/24/99 * @see KnownFloatQuantileEstimator * @see UnknownFloatQuantileEstimator */ public class FloatQuantileFinderFactory extends Object { /** * Make this class non instantiable. Let still allow others to inherit. */ protected FloatQuantileFinderFactory() { } /** * 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 * returnSamplingRate array, which, therefore must be at least of * length 1. * * @param N * the 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. 0.001) (0 <= epsilon <= 1). To * get exact result, set epsilon=0.0; * @param delta * the probability that the approximation error is more than than * epsilon (e.g. 0.0001) (0 <= delta <= 1 * ). To avoid probabilistic answers, set delta=0.0. * @param quantiles * the number of quantiles to be computed (e.g. 100) ( * quantiles >= 1). If unknown in advance, set this * number large, e.g. quantiles >= 10000. * @param returnSamplingRate * output parameter, a float[1] where the sampling rate * is to be filled in. * @return long[2] - long[0]=the number of buffers, * long[1]=the number of elements per buffer, * returnSamplingRate[0]=the required sampling rate. */ public static long[] known_N_compute_B_and_K(long N, float epsilon, float delta, int quantiles, float[] returnSamplingRate) { returnSamplingRate[0] = 1.0f; if (epsilon <= 0.0) { // no way around exact quantile search long[] result = new long[2]; result[0] = 1; result[1] = N; return result; } if (epsilon >= 1.0 || delta >= 1.0) { // can make any error we wish long[] result = new long[2]; result[0] = 2; result[1] = 1; return result; } if (delta > 0.0) { return known_N_compute_B_and_K_slow(N, epsilon, delta, quantiles, returnSamplingRate); } return known_N_compute_B_and_K_quick(N, epsilon); } /** * Computes the number of buffers and number of values per buffer such that * quantiles can be determined with a guaranteed approximation error * no more than epsilon. Assumes that quantiles are to be computed over N * values. * * @return long[2] - long[0]=the number of buffers, * long[1]=the number of elements per buffer. * @param N * the anticipated number of values over which quantiles shall be * determined. * @param epsilon * the approximation error which is guaranteed not to be exceeded * (e.g. 0.001) (0 <= epsilon <= 1). To * get exact result, set epsilon=0.0; */ protected static long[] known_N_compute_B_and_K_quick(long N, float epsilon) { final int maxBuffers = 50; final int maxHeight = 50; final float N_float = N; final float c = N_float * epsilon * 2.0f; int[] heightMaximums = new int[maxBuffers - 1]; // for each b, determine maximum height, i.e. the height for which x<=0 // and x is a maximum // with x = binomial(b+h-2, h-1) - binomial(b+h-3, h-3) + // binomial(b+h-3, h-2) - N * epsilon * 2.0 for (int b = 2; b <= maxBuffers; b++) { int h = 3; while (h <= maxHeight && // skip heights until x<=0 (h - 2) * (FloatArithmetic.binomial(b + h - 2, h - 1)) - (FloatArithmetic.binomial(b + h - 3, h - 3)) + (FloatArithmetic.binomial(b + h - 3, h - 2)) - c > 0.0) { h++; } // from now on x is monotonically growing... while (h <= maxHeight && // skip heights until x>0 (h - 2) * (FloatArithmetic.binomial(b + h - 2, h - 1)) - (FloatArithmetic.binomial(b + h - 3, h - 3)) + (FloatArithmetic.binomial(b + h - 3, h - 2)) - c <= 0.0) { h++; } h--; // go back to last height // was x>0 or did we loop without finding anything? int hMax; if (h >= maxHeight && (h - 2) * (FloatArithmetic.binomial(b + h - 2, h - 1)) - (FloatArithmetic.binomial(b + h - 3, h - 3)) + (FloatArithmetic.binomial(b + h - 3, h - 2)) - c > 0.0) { hMax = Integer.MIN_VALUE; } else { hMax = h; } heightMaximums[b - 2] = hMax; // safe some space } // end for // for each b, determine the smallest k satisfying the constraints, i.e. // for each b, determine kMin, with kMin = N/binomial(b+hMax-2,hMax-1) long[] kMinimums = new long[maxBuffers - 1]; for (int b = 2; b <= maxBuffers; b++) { int h = heightMaximums[b - 2]; long kMin = Long.MAX_VALUE; if (h > Integer.MIN_VALUE) { float value = (FloatArithmetic.binomial(b + h - 2, h - 1)); long tmpK = (long) (Math.ceil(N_float / value)); if (tmpK <= Long.MAX_VALUE) { kMin = tmpK; } } kMinimums[b - 2] = kMin; } // from all b's, determine b that minimizes b*kMin long multMin = Long.MAX_VALUE; int minB = -1; for (int b = 2; b <= maxBuffers; b++) { if (kMinimums[b - 2] < Long.MAX_VALUE) { long mult = (b) * (kMinimums[b - 2]); if (mult < multMin) { multMin = mult; minB = b; } } } long b, k; if (minB != -1) { // epsilon large enough? b = minB; k = kMinimums[minB - 2]; } else { // epsilon is very small or zero. b = 1; // the only possible solution without violating the k = N; // approximation guarantees is exact quantile search. } long[] result = new long[2]; result[0] = b; result[1] = k; return result; } /** * 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 returnSamplingRate 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. 0.001) (0 <= epsilon <= 1). To * get exact result, set epsilon=0.0; * @param delta * the probability that the approximation error is more than than * epsilon (e.g. 0.0001) (0 <= delta <= 1 * ). To avoid probabilistic answers, set delta=0.0. * @param quantiles * the number of quantiles to be computed (e.g. 100) ( * quantiles >= 1). If unknown in advance, set this * number large, e.g. quantiles >= 10000. * @param samplingRate * a float[1] where the sampling rate is to be filled * in. * @return long[2] - long[0]=the number of buffers, * long[1]=the number of elements per buffer, * returnSamplingRate[0]=the required sampling rate. */ protected static long[] known_N_compute_B_and_K_slow(long N, float epsilon, float delta, int quantiles, float[] returnSamplingRate) { final int maxBuffers = 50; final int maxHeight = 50; final float N_float = N; // One possibility is to use one buffer of size N // long ret_b = 1; long ret_k = N; float sampling_rate = 1.0f; 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 float logarithm = (float) Math.log(2.0 * quantiles / delta); final float c = 2.0f * epsilon * N_float; for (long b = 2; b < maxBuffers; b++) for (long h = 3; h < maxHeight; h++) { float binomial = FloatArithmetic.binomial(b + h - 2, h - 1); long tmp = (long) Math.ceil(N_float / binomial); if ((b * tmp < memory) && ((h - 2) * binomial - FloatArithmetic.binomial(b + h - 3, h - 3) + FloatArithmetic.binomial(b + h - 3, h - 2) <= c)) { ret_k = tmp; ret_b = b; memory = ret_k * b; sampling_rate = 1.0f; } if (delta > 0.0) { float t = (h - 2) * FloatArithmetic.binomial(b + h - 2, h - 1) - FloatArithmetic.binomial(b + h - 3, h - 3) + FloatArithmetic.binomial(b + h - 3, h - 2); float u = logarithm / epsilon; float v = FloatArithmetic.binomial(b + h - 2, h - 1); float w = logarithm / (2.0f * 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 float x = (float) (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_float * 2.0f * epsilon * epsilon / logarithm; } } } long[] result = new long[2]; result[0] = ret_b; result[1] = ret_k; returnSamplingRate[0] = sampling_rate; return result; } /** * Returns a quantile finder that minimizes the amount of memory needed * under the user provided constraints. * * Many applications don't know in advance over how many elements quantiles * are to be computed. However, some of them can give an upper limit, which * will assist the factory in choosing quantile finders with minimal memory * requirements. For example if you select values from a database and fill * them into histograms, then you probably don't know how many values you * will fill, but you probably do know that you will fill at most S * elements, the size of your database. * * @param known_N * specifies whether the number of elements over which quantiles * are to be computed is known or not. * @param N * if known_N==true, the number of elements over which * quantiles are to be computed. if known_N==false, the * upper limit on the number of elements over which quantiles are * to be computed. If such an upper limit is a-priori unknown, * then set N = Long.MAX_VALUE. * @param epsilon * the approximation error which is guaranteed not to be exceeded * (e.g. 0.001) (0 <= epsilon <= 1). To * get exact result, set epsilon=0.0; * @param delta * the probability that the approximation error is more than than * epsilon (e.g. 0.0001) (0 <= delta <= 1). To avoid * probabilistic answers, set delta=0.0. * @param quantiles * the number of quantiles to be computed (e.g. 100) ( * quantiles >= 1). If unknown in advance, set this * number large, e.g. quantiles >= 10000. * @param generator * a uniform random number generator. Set this parameter to * null to use a default generator. * @return the quantile finder minimizing memory requirements under the * given constraints. */ public static FloatQuantileFinder newFloatQuantileFinder(boolean known_N, long N, float epsilon, float delta, int quantiles, FloatRandomEngine generator) { // boolean known_N = true; // if (N==Long.MAX_VALUE) known_N = false; // check parameters. // if they are illegal, keep quite and return an exact finder. if (epsilon <= 0.0 || N < 1000) return new ExactFloatQuantileFinder(); if (epsilon > 1) epsilon = 1; if (delta < 0) delta = 0; if (delta > 1) delta = 1; if (quantiles < 1) quantiles = 1; if (quantiles > N) N = quantiles; KnownFloatQuantileEstimator finder; if (known_N) { float[] samplingRate = new float[1]; long[] resultKnown = known_N_compute_B_and_K(N, epsilon, delta, quantiles, samplingRate); long b = resultKnown[0]; long k = resultKnown[1]; if (b == 1) return new ExactFloatQuantileFinder(); return new KnownFloatQuantileEstimator((int) b, (int) k, N, samplingRate[0], generator); } else { long[] resultUnknown = unknown_N_compute_B_and_K(epsilon, delta, quantiles); long b1 = resultUnknown[0]; long k1 = resultUnknown[1]; long h1 = resultUnknown[2]; float preComputeEpsilon = -1.0f; if (resultUnknown[3] == 1) preComputeEpsilon = epsilon; // if (N==Long.MAX_VALUE) { // no maximum N provided by user. // if (true) fixes bug reported by [email protected] if (true) { // no maximum N provided by user. if (b1 == 1) return new ExactFloatQuantileFinder(); return new UnknownFloatQuantileEstimator((int) b1, (int) k1, (int) h1, preComputeEpsilon, generator); } // determine whether UnknownFinder or KnownFinder with maximum N // requires less memory. float[] samplingRate = new float[1]; // IMPORTANT: for known finder, switch sampling off (delta == 0) !!! // with knownN-sampling we can only guarantee the errors if the // input sequence has EXACTLY N elements. // with knownN-no sampling we can also guarantee the errors for // sequences SMALLER than N elements. long[] resultKnown = known_N_compute_B_and_K(N, epsilon, 0, quantiles, samplingRate); long b2 = resultKnown[0]; long k2 = resultKnown[1]; if (b2 * k2 < b1 * k1) { // the KnownFinder is smaller if (b2 == 1) return new ExactFloatQuantileFinder(); return new KnownFloatQuantileEstimator((int) b2, (int) k2, N, samplingRate[0], generator); } // the UnknownFinder is smaller if (b1 == 1) return new ExactFloatQuantileFinder(); return new UnknownFloatQuantileEstimator((int) b1, (int) k1, (int) h1, preComputeEpsilon, generator); } } /** * Convenience method that computes phi's for equi-depth histograms. This is * simply a list of numbers with i / (float)quantiles for * i={1,2,...,quantiles-1}. * * @return the equi-depth phi's */ public static cern.colt.list.tfloat.FloatArrayList newEquiDepthPhis(int quantiles) { cern.colt.list.tfloat.FloatArrayList phis = new cern.colt.list.tfloat.FloatArrayList(quantiles - 1); for (int i = 1; i <= quantiles - 1; i++) phis.add(i / (float) quantiles); return phis; } /** * 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. * * @param epsilon * the approximation error which is guaranteed not to be exceeded * (e.g. 0.001) (0 <= epsilon <= 1). To * get exact results, set epsilon=0.0; * @param delta * the probability that the approximation error is more than than * epsilon (e.g. 0.0001) (0 <= delta <= 1 * ). To get exact results, set delta=0.0. * @param quantiles * the number of quantiles to be computed (e.g. 100) ( * quantiles >= 1). If unknown in advance, set this * number large, e.g. quantiles >= 10000. * @return long[4] - long[0]=the number of buffers, * long[1]=the number of elements per buffer, * long[2]=the tree height where sampling shall start, * long[3]==1 if precomputing is better, otherwise 0; */ public static long[] unknown_N_compute_B_and_K(float epsilon, float delta, int quantiles) { return unknown_N_compute_B_and_K_raw(epsilon, delta, quantiles); // move stuff from _raw(..) here and delete _raw(...) /* * long[] result_1 = * unknown_N_compute_B_and_K_raw(epsilon,delta,quantiles); long b1 = * result_1[0]; long k1 = result_1[1]; * * * int quantilesToPrecompute = (int) Floats.ceiling(1.0 / epsilon); * * if (quantiles>quantilesToPrecompute) { // try if precomputing * quantiles requires less memory. long[] result_2 = * unknown_N_compute_B_and_K_raw(epsilon/2.0,delta,quantilesToPrecompute); * * long b2 = result_2[0]; long k2 = result_2[1]; if (b2*k2 < b1*k1) { * result_2[3] = 1; //precomputation is better result_1 = result_2; } } * return result_1; */ } /** * 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. You never need to call this * method. It is only for curious users wanting to gain some insight * into the workings of the algorithms. * * @param epsilon * the approximation error which is guaranteed not to be exceeded * (e.g. 0.001) (0 <= epsilon <= 1). To * get exact result, set epsilon=0.0; * @param delta * the probability that the approximation error is more than than * epsilon (e.g. 0.0001) (0 <= delta <= 1 * ). To get exact results, set delta=0.0. * @param quantiles * the number of quantiles to be computed (e.g. 100) ( * quantiles >= 1). If unknown in advance, set this * number large, e.g. quantiles >= 10000. * @return long[4] - long[0]=the number of buffers, * long[1]=the number of elements per buffer, * long[2]=the tree height where sampling shall start, * long[3]==1 if precomputing is better, otherwise 0; */ protected static long[] unknown_N_compute_B_and_K_raw(float epsilon, float delta, int quantiles) { // delta can be set to zero, i.e., all quantiles should be approximate // with probability 1 if (epsilon <= 0.0) { long[] result = new long[4]; result[0] = 1; result[1] = Long.MAX_VALUE; result[2] = Long.MAX_VALUE; result[3] = 0; return result; } if (epsilon >= 1.0 || delta >= 1.0) { // can make any error we wish long[] result = new long[4]; result[0] = 2; result[1] = 1; result[2] = 3; result[3] = 0; return result; } if (delta <= 0.0) { // no way around exact quantile search long[] result = new long[4]; result[0] = 1; result[1] = Long.MAX_VALUE; result[2] = Long.MAX_VALUE; result[3] = 0; return result; } int max_b = 50; int max_h = 50; int max_H = 50; int max_Iterations = 2; long best_b = Long.MAX_VALUE; long best_k = Long.MAX_VALUE; long best_h = Long.MAX_VALUE; long best_memory = Long.MAX_VALUE; float pow = (float) Math.pow(2.0, max_H); float logDelta = (float) (Math.log(2.0 / (delta / quantiles)) / (2.0 * epsilon * epsilon)); // float logDelta = Math.log(2.0/(quantiles*delta)) / // (2.0*epsilon*epsilon); while (best_b == Long.MAX_VALUE && max_Iterations-- > 0) { // until we // find a // solution // identify that combination of b and h that minimizes b*k. // exhaustive search. for (int b = 2; b <= max_b; b++) { for (int h = 2; h <= max_h; h++) { float Ld = FloatArithmetic.binomial(b + h - 2, h - 1); float Ls = FloatArithmetic.binomial(b + h - 3, h - 1); // now we have k>=c*(1-alpha)^-2. // let's compute c. // float c = Math.log(2.0/(delta/quantiles)) / // (2.0*epsilon*epsilon*Math.min(Ld, 8.0*Ls/3.0)); float c = (float) (logDelta / Math.min(Ld, 8.0 * Ls / 3.0)); // now we have k>=d/alpha. // let's compute d. float beta = Ld / Ls; float cc = (float) ((beta - 2.0) * (max_H - 2.0) / (beta + pow - 2.0)); float d = (float) ((h + 3 + cc) / (2.0 * epsilon)); /* * float d = (Ld*(h+max_H-1.0) + Ls*((h+1)*pow - * 2.0*(h+max_H))) / (Ld + Ls*(pow-2.0)); d = (d + 2.0) / * (2.0*epsilon); */ // now we have c*(1-alpha)^-2 == d/alpha. // we solve this equation for alpha yielding two solutions // alpha_1,2 = (c + 2*d +- Sqrt(c*c + 4*c*d))/(2*d) float f = c * c + 4.0f * c * d; if (f < 0.0) continue; // non real solution to equation float root = (float) Math.sqrt(f); float alpha_one = (float) ((c + 2.0 * d + root) / (2.0 * d)); float alpha_two = (float) ((c + 2.0 * d - root) / (2.0 * d)); // any alpha must satisfy 0 0) { // valid solution? long memory = b * k; if (memory < best_memory) { // found a solution requiring less memory best_k = k; best_b = b; best_h = h; best_memory = memory; } } } } // end for h } // end for b if (best_b == Long.MAX_VALUE) { System.out.println("Warning: Computing b and k looks like a lot of work!"); // no solution found so far. very unlikely. Anyway, try again. max_b *= 2; max_h *= 2; max_H *= 2; } } // end while long[] result = new long[4]; result[3] = 0; if (best_b == Long.MAX_VALUE) { // no solution found. // no way around exact quantile search. result[0] = 1; result[1] = Long.MAX_VALUE; result[2] = Long.MAX_VALUE; } else { result[0] = best_b; result[1] = best_k; result[2] = best_h; } return result; } }





© 2015 - 2025 Weber Informatics LLC | Privacy Policy