distributions.cc

经网络提出了一种基于蚁群聚类算法的径向基神经网络. 利用蚁群算法的并行寻优特征和挥发系数方法的自适应更改信息量的能力,并以球面聚类的方式确定了径向基神经网络中基函数的位置, 同时通过比较隐层神经元的相

/* $Id: distributions.cc,v 1.3 2006-08-12 20:50:06 jonathan Exp $
* Jonathan Ledlie, Harvard University.
* Copyright 2005.  All rights reserved.
*/

#include <stdio.h>
#include <stdlib.h>
#include "error.h"
#include "distributions.h"

/*
* Parts of this are borrowed (with thanks) from:
*
* Perl Math::Random
* Math::Random was put together by John Venier and Barry W. Brown
* with help from SWIG.  For version 0.61, Geoffrey Rommel made various
* cosmetic changes.
*
* Kenneth J. Christensen
* Associate Professor
* Department of Computer Science and Engineering
* University of South Florida
*
*/

float INVALID = -10099.99104;

double ParetoDistribution::next ()
{
	return (scale/(pow(randPct(),(1./shape))));
};

double ParetoDistribution::getMean() 
{
	if (shape > 1.)
		return (scale * shape) / (shape - 1.);
	else
		return 0;
}

double NormalDistribution::next ()
{
	double normal = stddev*snorm()+mean;
	return normal;
};

double NormalDistribution::getMean() 
{
	return mean;
}

double UniformDistribution::next () 
{
	return randPct();
};

double UniformDistribution::getMean() 
{
	return 0.;
}

int UniformIntegerDistribution::next () 
{
	return unifRand(0,maxElement);
};

double UniformIntegerDistribution::getMean()
{
	return maxElement/2.;
}

double GnutellaBWDistribution::next () 
{
	double z = randPct();

	if (z < .2)
	{
		return 1.0672;
	} 
	else if (z < .65) 
	{
		return 10.672;
	}
	else if (z < .95)
	{
		return 106.72;
	}
	else if (z < .999) 
	{
		return 1067.2;
	} 
	else
	{
		return 10672;
	}
};

double GnutellaBWDistribution::getMean() 
{
	ASSERT (0);
	return 0.;
}

double computeZipfNormalization (double alpha, int elements) 
{
	double c = 0.;
	for (int i=1; i<=elements; i++)
	{
		c = c + (1.0 / pow((double) i, alpha));
	}
	c = 1.0 / c;
	return c;
}

ZipfDistribution::ZipfDistribution (double a, int n) : alpha(a), elements(n)
{
	c = computeZipfNormalization (alpha, elements);

	// Map z to the value
	double sum_prob = 0.;
	for (int i=1; i<=elements; i++) 
	{
		sum_prob = sum_prob + c / pow((double) i, alpha);
		zipfKeys.insert(pair<double,double>(sum_prob,randPct()));
	}
}

ZipfIntegerDistribution::ZipfIntegerDistribution (double a, int n) 
			: alpha(a), elements(n) 
{
	c = computeZipfNormalization (alpha, elements);

	vector<int> elementList;
	for (int i=0; i<elements; i++)
	{
		elementList.push_back(i);
	}
	//把elementList内的数据随机队顺一次
	random_shuffle (elementList.begin(),elementList.end());

	// Map z to the value
	double sum_prob = 0.;
	for (int i=1; i<=elements; i++) 
	{
		sum_prob = sum_prob + c / pow((double) i, alpha);
		zipfKeys.insert(pair<double,int>(sum_prob,elementList[i-1]));
	}
}


double ZipfDistribution::next () 
{
	double z;                     // Uniform random number (0 < z < 1)
	z = randPct();

	map<double,double>::iterator p;
	p = zipfKeys.upper_bound (z);//在zipfkeys中返回一个比Z大的值的位置
	if (p == zipfKeys.end())
		p = zipfKeys.begin();

	return p->second;
};

int ZipfIntegerDistribution::next ()
{
	double z;                     // Uniform random number (0 < z < 1)
	z = randPct();

	map<double,int>::iterator p;
	p = zipfKeys.upper_bound (z);
	if (p == zipfKeys.end())
		p = zipfKeys.begin();

	return p->second;
};


double ZipfDistribution::getMean() 
{
	return alpha;
}

double ZipfIntegerDistribution::getMean() 
{
	return alpha;
}

double PoissonDistribution::next () 
{

	double u;
	while ((u = randPct()) == 0.)
		;
	return -mean * log(u);
};

double PoissonDistribution::getMean() 
{
	return mean;
}

double ConstantDistribution::next ()
{
	return value;
};

double ConstantDistribution::getMean() 
{
	return value;
}

IntegerDistribution* DistributionFactory::newIntegerDistribution (char *str)
{
	char distType;
	float para1, para2;
	getParameters (str, distType, para1, para2);
	//printf ("%c p1 %f p2 %f\n", distType, para1, para2);

	switch (distType)
	{
	case 'Z':
		if (para1 == INVALID || para2 == INVALID)
			return NULL;
		return new ZipfIntegerDistribution ((double)para1, (int)para2);
		break;
	case 'U':
		if (para1 == INVALID || para2 != INVALID)
			return NULL;
		return new UniformIntegerDistribution ((int)para1);
		break;
	default:
		return NULL;    
	}

	return NULL;
}


// str[in]  distType[out] para1[out]  para[out]
void getParameters (char *str, char &distType, float &para1, float &para2) 
{
	para1 = INVALID;
	para2 = INVALID;
	int ret;

	if ( (sscanf (str, "%c/%f/%f", &distType, &para1, &para2) != 3) &&
		(sscanf (str, "%c/%f", &distType, &para1) != 2) &&
		(sscanf (str, "%c", &distType) != 1))
	{
		return;
	}
}

Distribution* DistributionFactory::newDistribution (char *str)
{
	char distType;
	float para1, para2;
	getParameters (str, distType, para1, para2);

	switch (distType) 
	{
	case 'P':
		if (para1 == INVALID || para2 == INVALID)
			return NULL;
		return new ParetoDistribution (para1, para2);
		break;
	case 'N':
		if (para1 == INVALID || para2 == INVALID)
			return NULL;
		return new NormalDistribution (para1, para2);
		break;
	case 'Z':
		if (para1 == INVALID || para2 == INVALID) 
			return NULL;
		return new ZipfDistribution ((double)para1, (int)para2);
		break;
	case 'U':
		if (para1 != INVALID || para2 != INVALID) 
			return NULL;
		return new UniformDistribution ();
		break;
	case 'F':
		if (para1 == INVALID || para2 != INVALID)
			return NULL;
		return new PoissonDistribution (para1);
		break;
	case 'C':
		if (para1 == INVALID || para2 != INVALID)
			return NULL;
		return new ConstantDistribution (para1);
		break;
	case 'G':
		if (para1 != INVALID || para2 != INVALID) 
			return NULL;
		return new GnutellaBWDistribution ();
		break;
	default:
		return NULL;    
	}

	return NULL;
}

double snorm(void)
/*
**********************************************************************


(STANDARD-)  N O R M A L  DISTRIBUTION                           


**********************************************************************
**********************************************************************

FOR DETAILS SEE:                                                 

AHRENS, J.H. AND DIETER, U.                            
EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM             
SAMPLING FROM THE NORMAL DISTRIBUTION.                 
MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.          

ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  
(M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  

Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   
SUNIF.  The argument IR thus goes away.                          

**********************************************************************
THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
*/
{
	static double a[32] = {
		0.0,              0.03917608550309, 0.07841241273311, 0.11776987457909,
		0.15731068461017, 0.19709908429430, 0.23720210932878, 0.27769043982157,
		0.31863936396437, 0.36012989178957, 0.40225006532172, 0.44509652498551,
		0.48877641111466, 0.53340970624127, 0.57913216225555, 0.62609901234641,
		0.67448975019607, 0.72451438349236, 0.77642176114792, 0.83051087820539,
		0.88714655901887, 0.94678175630104, 1.00999016924958, 1.07751556704027,
		1.15034938037600, 1.22985875921658, 1.31801089730353, 1.41779713799625,
		1.53412054435253, 1.67593972277344, 1.86273186742164, 2.15387469406144
	};
	static double d[31] = {
		0.0, 0.0, 0.0, 0.0,
		0.0,              0.26368432217502, 0.24250845238097, 0.22556744380930,
		0.21163416577204, 0.19992426749317, 0.18991075842246, 0.18122518100691,
		0.17360140038056, 0.16684190866667, 0.16079672918053, 0.15534971747692,
		0.15040938382813, 0.14590257684509, 0.14177003276856, 0.13796317369537,
		0.13444176150074, 0.13117215026483, 0.12812596512583, 0.12527909006226,
		0.12261088288608, 0.12010355965651, 0.11774170701949, 0.11551189226063,
		0.11340234879117, 0.11140272044119, 0.10950385201710
	};
	static double t[31] = {
		7.6738283767E-4,   2.30687039764E-3,  3.86061844387E-3,  5.43845406707E-3,
		7.05069876857E-3,  8.70839582019E-3,  1.042356984914E-2, 1.220953194966E-2,
		1.408124734637E-2, 1.605578804548E-2, 1.815290075142E-2, 2.039573175398E-2,
		2.281176732513E-2, 2.543407332319E-2, 2.830295595118E-2, 3.146822492920E-2,
		3.499233438388E-2, 3.895482964836E-2, 4.345878381672E-2, 4.864034918076E-2,
		5.468333844273E-2, 6.184222395816E-2, 7.047982761667E-2, 8.113194985866E-2,
		9.462443534514E-2, 0.11230007889456,  0.13649799954975,  0.17168856004707,
		0.22762405488269,  0.33049802776911,  0.58470309390507
	};
	static double h[31] = {
		3.920617164634E-2, 3.932704963665E-2, 3.950999486086E-2, 3.975702679515E-2,
		4.007092772490E-2, 4.045532602655E-2, 4.091480886081E-2, 4.145507115859E-2,
		4.208311051344E-2, 4.280748137995E-2, 4.363862733472E-2, 4.458931789605E-2,
		4.567522779560E-2, 4.691571371696E-2, 4.833486978119E-2, 4.996298427702E-2,
		5.183858644724E-2, 5.401138183398E-2, 5.654656186515E-2, 5.953130423884E-2,
		6.308488965373E-2, 6.737503494905E-2, 7.264543556657E-2, 7.926471414968E-2,
		8.781922325338E-2, 9.930398323927E-2, 0.11555994154118,  0.14043438342816,
		0.18361418337460,  0.27900163464163,  0.70104742502766
	};
	static long i;
	static double snormboth,u,s,ustar,aa,w,y,tt;
	//    u = ranf();

	u = randPct();//返回0-1之间的随机数
	s = 0.0;
	if(u > 0.5)
		s = 1.0;
	u += (u-s);
	u = 32.0*u;
	i = (long) (u);
	if(i == 32)
		i = 31;
	if(i == 0) 
		goto S100;
	/*
	START CENTER
	*/
	ustar = u-(double)i;
	aa = *(a+i-1);
S40:
	if(ustar <= *(t+i-1)) goto S60;
	w = (ustar-*(t+i-1))**(h+i-1);
S50:
	/*
	EXIT   (BOTH CASES)
	*/
	y = aa+w;
	snormboth = y;
	if(s == 1.0) 
		snormboth = -y;
	return snormboth;
S60:
	/*
	CENTER CONTINUED
	*/
	//    u = ranf();
	u = randPct();
	w = u*(*(a+i)-aa);
	tt = (0.5*w+aa)*w;
	goto S80;
S70:
	tt = u;
	//    ustar = ranf();
	ustar = randPct();
S80:
	if(ustar > tt) goto S50;
	//    u = ranf();
	u = randPct();
	if(ustar >= u) goto S70;
	//    ustar = ranf();
	ustar = randPct();
	goto S40;
S100:
	/*
	START TAIL
	*/
	i = 6;
	aa = *(a+31);
	goto S120;
S110:
	aa += *(d+i-1);
	i += 1;
S120:
	u += u;
	if(u < 1.0)
		goto S110;
	u -= 1.0;
S140:
	w = u**(d+i-1);
	tt = (0.5*w+aa)*w;
	goto S160;
S150:
	tt = u;
S160:
	//    ustar = ranf();
	ustar = randPct();
	if(ustar > tt) 
		goto S50;
	//    u = ranf();
	u = randPct();
	if(ustar >= u)
		goto S150;
	//    u = ranf();
	u = randPct();
	goto S140;
}