rand.c

Solaris环境下的数据挖掘算法：birch聚类算法。该算法适用于对大量数据的挖掘。

    sum += value;
    sumOfSqu += value*value;
    if (value > sampleMax) sampleMax = value;
    if (value < sampleMin) sampleMin = value;
    };

void FloatRandomStream::CollectStat(float value) {
    count++;
    sum += value;
    sumOfSqu += value*value;
    if (value > sampleMax) sampleMax = value;
    if (value < sampleMin) sampleMin = value;
    };

DiscreteUniform::DiscreteUniform(int n, int stream) {
    this->stream = stream;
    this->n   = n;
    statFlag  = 2; 	
    mean      = (n+1)*0.5;
    if (n == 1){
	std = 0;
    } else {
        std = sqrt((n*n-1)/(12.0));
    };
    };

int  DiscreteUniform::Draw() {
    int result;
    float sample = rint(n)*Uniform(this->stream); 
#if defined(linux ) || (defined(i386) && defined(sun))
    result = int(rint(floor(sample))) + 1;  // I hope that was the right idea
#else
    result = irint(floor(sample)) + 1;
#endif
    if (statFlag  == 2)  CollectStat(result);
    return result;
    };

void DiscreteUniform::PrintProfile() {
	printf("\n DiscreteUniform [1-%d] generator (%f, %f)",n,mean,std);
    };

Choose::Choose(float prob, int stream) {
    this->stream = stream;
    this->prob   = prob;
    statFlag     = 2;
    mean         = prob;
    std          = sqrt(prob*(1-prob));
    };

int Choose::Draw() {
    int	sample;
    sample = (Uniform(this->stream) <= prob) ? 1 : 0;
    if (statFlag  & 2) CollectStat(sample);
    return sample;
    };

void Choose::PrintProfile() {
    printf("\n Choose[%f] generator (%f, %f)",prob,mean,std);
    }

Uniform::Uniform(float low, float high, int stream) {
    this->stream = stream;
    this->low  = low;
    this->high = high;
    statFlag   = 2; 	
    mean       = (low+high)*0.5;
    std        = (high-low)/sqrt(12.0);
    };

float	Uniform::Draw() {
    double sample = low+(high-low)*Uniform(this->stream); 
    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };

void Uniform::PrintProfile() {
	printf("\n Uniform [%f,%f] generator (%f, %f)",low,high,mean,std);
    };

Uniform01::Uniform01(int stream) {
    this->stream = stream;
    statFlag     = 2; 	
    mean         = 0.5;
    std          = sqrt(1.0/12.0);
    };

float	Uniform01::Draw() {
    double sample = Uniform(this->stream); 
    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };

void Uniform01::PrintProfile() {
	printf("\n Uniform [0,1] generator (%f, %f)",mean,std);
    };

Exponential::Exponential(float mean, int stream) {
	this->stream = stream;
	this->mean = mean;
	this->std = mean;
	statFlag  = 2;
}

float	Exponential::Draw()
{
    float sample = -mean * ( log ( Uniform(this->stream)));
    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };

void Exponential::PrintProfile() {
	printf("\n Exponential(%f) generator (%f, %f)",mean,mean,std);
    };

HyperExp::HyperExp(float mean, float variance, int stream) 
{
	float	cx2;

	this->stream = stream;
	this->mean = mean;
	this->std = sqrt(variance);
	cx2  = variance/ (mean * mean);
	alpha = 1/(1 + cx2);
	mu1  =  mean * (1/( 1 + sqrt( 1/(alpha - 1.0)*(cx2 -1)/2)));
	mu2  =  mean * (1/( 1 - sqrt((1/(1 - alpha) - 1) * (cx2 -1)/2)));
}

HyperExp::HyperExp(float alpha, float mu1, float mu2, int stream) 
{
	this->stream = stream;
	this->alpha = alpha;
	this->mu1 = mu1;
	this->mu2 = mu2;

	mean = alpha*mu1 + (1 - alpha)*mu2;
	std = sqrt(alpha*mu1*mu1 + (1 - alpha)*mu2*mu2);
	std = sqrt(2*(alpha*mu1*mu1 + (1 - alpha)*mu2*mu2)
	        - (alpha*mu1+(1-alpha)*mu2)*(alpha*mu1+(1-alpha)*mu2));
}

float	HyperExp::Draw()
{
	float	sample;

	if (Uniform(stream) <= alpha)
	    sample = -log(Uniform(stream))*mu1;
	else	
	    sample = -log(Uniform(stream))*mu2;
	CollectStat(sample);
	return sample;
}

void	HyperExp::PrintProfile() {
	printf("\nHyperexponential[%f, %f,%f] generator (%f, %f)",
	    alpha,mu1,mu2,mean,std);
}

Normal::Normal(float mean, float sigma, int stream) {
	this->stream = stream;
	this->mean = mean;
	this->std = sigma;
	statFlag  = 2;
    };

float  Normal::Draw() {
    float	v1, v2;
    float	sample;

    do {
        v1 = - log(Uniform(this->stream));
        v2 = - log(Uniform(this->stream));
        } while( (2.0* v1) < ((v2 - 1.0)*(v2 - 1.0)));

    if (Uniform(this->stream) < 0.5)
        v2 *= -1;
    sample = std * v2 + mean;
	
    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };
void  Normal::PrintProfile() {
	printf("\n Normal[%f,%f] generator (%f, %f)",mean,std,mean,std);
    };

LogNormal::LogNormal(float mean, float sigma, int stream) {
	this->stream = stream;
	this->mean = mean;
	this->std = sigma;
	statFlag  = 2;
    };

float  LogNormal::Draw() {
    float	v1, v2;
    float	sample;

    do {
        v1 = - log(Uniform(this->stream));
        v2 = - log(Uniform(this->stream));
        } while( (2.0* v1) < ((v2 - 1.0)*(v2 - 1.0)));

    if (Uniform(this->stream) < 0.5)
        v2 *= -1;
    sample = exp(std * v2 + mean);
	
    if (statFlag  == 2) CollectStat(sample);
    return sample;
    };

void  LogNormal::PrintProfile() {
	printf("\n LogNormal[%f,%f] generator (%f, %f)",mean,std,mean,std);
    };

Gamma::Gamma(float alpha, float beta, int stream) {
	this->alpha = alpha;
	this->beta  = beta;
	mean        = alpha/beta;
	std         = sqrt(alpha)/beta;
	statFlag  = 2;
}

float	GammaDist(int stream, float alpha, float beta) {
	float a, v1, v2;

	a = alpha - 1.0;
	do {
	    v1 = -log( Uniform(stream));
	    v2 = -log( Uniform(stream));
	} while ( v2 < a* (v1 - log(v1) - 1.0));

	return alpha/beta*v1;
    };

float	Gamma::Draw() {
	float sample = GammaDist(this->stream,alpha,beta);
	if (statFlag  == 2) CollectStat(sample);
	return sample;
    };

void Gamma::PrintProfile()
{
	printf("\n Gamma [%f, %f]  generator (%f, %f)",alpha,beta,mean,std);
}

Beta::Beta(float alpha, float beta, int stream) 
{
	this->alpha  = alpha;
	this->beta  = beta;
	statFlag     = 2;
	mean = alpha/(alpha + beta);
	std = sqrt((alpha*beta / (alpha + beta + 1))/((alpha + beta)*(alpha + beta)));
}

float Beta::Draw()
{
	float	y1, y2, sample;

	y1 = GammaDist(stream,alpha,beta);
	y2 = GammaDist(stream,alpha,beta);
	sample  = y1/(y1+y2);
	if (statFlag  == 2) CollectStat(sample);
	return sample;
}

void Beta::PrintProfile() {
    printf("\n Beta[%f, %f] generator (%f, %f)",alpha,beta,mean,std);
    };

GeneralCDF::GeneralCDF(int n, float* v, float* p , int stream) {
    int	i;

    this->n = n;
    this->v = new float[n];
    this->p = new float[n];
    q       = new float[n];

    this->v[0] = v[0];
    this->p[0] = p[0];
    this->q[0] = p[0];
    for(i=1; i < n; i++){ 
	this->v[i] = v[i];
	this->p[i] = p[i];
        this->q[i] = p[i] + p[i-1];
        }
    }

float	GeneralCDF::Draw() {
	float	sample, result;

	sample = Uniform(stream);
	result = v[binary_search(0, n-1, q, sample)];
	CollectStat(result);
	return result;
}

GeneralCDF::~GeneralCDF()
{
	delete this->v;
	delete this->p;
	delete this->q;
}

void	GeneralCDF::PrintProfile()
{
	printf("\nGeneral CDF generator (%f, %f)",mean,std);
}

int	binary_search(int low, int high, float* q, float val)
{
	int	mid;

	while (low <= high) {
		mid = (low + high)/2;
		if (q[mid] == val) return mid;
		else if (q[mid] < val) {
			low = mid + 1;
			return binary_search(low, high, q, val);
		}
		else {
			high = mid - 1;
			return binary_search(low, high, q, val);
		}
	}
	return low;
}