main.c

这是Clerc最新的文章balanced PSO,包括论文和与其配套的源程序

	switch(param.strat)
	{
		default: break;
		case 1: case 2:
		printf("\n w0Rate %f",param.w0Rate);
		 break;

	}

	printf("\n c1 = %f, c2 = %f",param.c1,param.c2);
if(param.R==1) printf("\nStandard deviation %f",param.sigma);
	printf("\nLocal search %i",param.localSearch);


	switch(param.localSearch)
	{
		case 1:
		case 2:
		printf("\n Initial exploitation rate reference %f",param.exploit1);
	 break;

		case 3:
		printf("\n Exploitation area relative side length: either %f or %f",1./param.S,param.exploit2);
		case 4:
		printf("\n Initial exploitation rate reference %f",param.exploit1);
				break;
	}

	switch(param.explOption)
	{
		default:
		break;
		case 3:
		printf("\n Exploitation area relative side length: either %f or %f",1./param.S,param.exploit2);		
		break;
		case 4:
		printf("\n Exploitation area relative side length: uniform in [0,%f]",param.exploit2);
		break;
		case 5:
		printf("\n Exploitation area relative side length: Gaussian (%f,%f)",1./param.S,param.exploit2);
		break;
	}
	//---------------
	errorMean = 0;	    
	evalMean = 0;	    
	nFailure = 0;	
	//------------------------------------- RUNS	
	for (run = 0; run < runMax; run++)  
	{	
		FEclass=0;
		//srand (clock () / 100);	// May improve pseudo-randomness 
       
		result = PSO (param, pb, 0);
		error = result.error;
	
		if (error > pb.epsilon) // Failure
		{
			nFailure = nFailure + 1;
		}
		
		// Result display
		printf ("\nRun %i. Eval %f. Error %f ", run+1, result.nEval, error);
		printf("  x(0)= %f",(double)result.SW.P[result.SW.best].x[0]);
			// Save result
			 fprintf( f_run, "\n%i %f %0.10f ", run, result.nEval,  error );
			 for ( d = 0; d < pb.SS.D; d++ ) fprintf( f_run, " %0.10f",  (double)result.SW.P[result.SW.best].x[d] );
			
	
		// Compute/store some statistical information
		if (run == 0)
		{
			errorMin = error;
			bestBest=result.SW.P[result.SW.best];
		}
		else if (error < errorMin)
		{
			errorMin = error;
			bestBest=result.SW.P[result.SW.best];
		}

		evalMean = evalMean + result.nEval;	
		errorMean = errorMean + error;	
		errorMeanBest[run] = error;
		logProgressMean  = logProgressMean - log(error);		
	}		// End loop on "run"


	// Display statistical information
	// Mean
	evalMean = evalMean / (double) runMax;   
	errorMean = errorMean / (double) runMax;
	logProgressMean = logProgressMean/(double) runMax;
			
	printf ("\n Eval. (mean)= %f", evalMean);	
	printf ("\n Error (mean) = %.12f", errorMean);
		
	// Variance
	variance = 0;
			
	for (run = 0; run < runMax; run++)
				variance = variance + pow (errorMeanBest[run] - errorMean, 2);
			
	variance = sqrt (variance / runMax);	    
	printf ("\n Std. dev. %f", variance); 
	printf("\n Log_progress (mean) = %f", logProgressMean);	
	// Success rate and minimum value
	printf("\n Failure(s) %i",nFailure);
	successRate = 100 * (1 - nFailure / (double) runMax);			
	printf ("\n Success rate = %.2f%%", successRate);
			
	if (run > 1)
		printf ("\n Best min value = %f", errorMin);
			
	// Save
	/*
		fprintf(f_synth,"\n"); for (d=0;d<SS.D;d++) fprintf(f_synth,"%f ",
			pb.offset[d]);
		fprintf(f_synth,"    %f %f %f %.0f%% %f",errorMean,variance,errorMin,
			successRate,evalMean); 

	 fprintf( f_synth, "\n%f %f %f %f %.0f%% %f ", shift,
			errorMean, variance, errorMin, successRate, evalMean );
	*/
	fprintf (f_synth, "\n");		
	fprintf (f_synth, "%f %f %.0f%% %f   ",
					errorMean, variance, successRate, evalMean);		   
	
	// Save the best result
	fprintf( f_synth, "\n%f ", errorMin );
	for ( d = 0; d < pb.SS.D; d++ ) fprintf( f_synth, " %f", (double) bestBest.x[d] );

	// Save the position in the "init" format, if functionCode<0

	if(functionCode<0)
	{
		dim=pb.SS.D/pb.nb;
		fprintf(f_init_save,"%i\n",pb.nb);
		for(n=0;n<pb.nb;n++)
		{
			for(d=0;d<dim;d++) 	fprintf(f_init_save,"%f ",(double) bestBest.x[d+n*dim]); 
			fprintf(f_init_save,"\n");
		}
	}
		
	// Save the curve "success rate" vs FEmax
	for(j=0;j<FEmaxSize-1;j++)
	{
		zz=0;
		for(run=0;run<runMax;run++)
		{
			zz=zz+successFlag[run][j];
		}
		zz=zz/runMax;
		fprintf(fCurve,"%i %f\n",FEmax[j],zz);
	}


	return 0; // End of main program
}
// ======================================================== PSO
struct result PSO (struct param param, struct problem pb, int level) 
{  
	struct velocity aleaV; 
	double c1,c2; 
	int d; 
	int deb;
	double error;   
	double errorPrev;
	struct exploit ex;
	double explBestSize;
	int explForced;
	double explRate0, explRate;
	int g[S_max]; // Best informants
	int sBest; // Rank in g of the best of the bests
	struct velocity GX;   
	int index[S_max], indexTemp[S_max];     
	int initLinks;	// Flag to (re)init or not the information links
	int iter; 		// Iteration number (time step)
	int iterBegin;
	int length;
	int LINKS[S_max][S_max];	// Information links
	struct position local;
	int m; 
	int noEval; 	
	int noStop;
	int outside;
	double p;
	struct velocity PX;	
	struct result R;
	int rank;
	double rho;
	int s0=0;
	int s,s1; 
	int t;
	float z;
	double zz,zzz;

	aleaV.size=pb.SS.D;

	// -----------------------------------------------------
	// INITIALISATION

	p=param.p; // Probability threshold for random topology
	R.SW.S = param.S; // Size of the current swarm
	local.size=pb.SS.D;
	
	// Position
		
	switch(param.init)
	{
		case -1: // Read from file 
		f_init=fopen("f_initS20D30.txt","r");
		fscanf(f_init,"%i",&s0);
		if(s0!=R.SW.S)
		{
			printf("\nThe number of points (%i) in the file",s0);
			printf("\nis not consistent with the swarm size S (%i)",R.SW.S);
			ERROR(" ");
		}
		for (s = 0; s < R.SW.S; s++)   
		{	
			for (d = 0; d < pb.SS.D; d++)  
			{ 
				fscanf(f_init,"%f",&z);
				R.SW.X[s].x[d]=z;		
			}
		}	
		break;

		default: // 0. Uniform
		for (s = 0; s < R.SW.S; s++)   
		{	
			for (d = 0; d < pb.SS.D; d++)  
			{  
				R.SW.X[s].x[d] = alea (pb.SS.min[d], pb.SS.max[d]);
			}
		}
		break;

		case 1: // Centre biased
		for (s = 0; s < R.SW.S; s++)   
		{

			for (d = 0; d < pb.SS.D; d++)  
			{  
				zz=0.5*(pb.SS.max[d]+pb.SS.min[d]);
				if(alea(0,1)<0.5) R.SW.X[s].x[d]=zz*pow(alea(0,1),1./pb.SS.D);
				else R.SW.X[s].x[d]=zz*(2-pow(alea(0,1),1./pb.SS.D));
			}

		}
		break;
		
		case 2: // Improved Hammersley
		deb=alea_integer(0,pb.SS.D-1);
		for (s = 0; s < R.SW.S; s++)   
		{
				R.SW.X[s]=hammersley(pb.SS,R.SW.S, deb,s);
		}
		break;

		case 3: // Random or Improved Hammersley	
		deb=alea_integer(0,pb.SS.D-1);
		for (s = 0; s < R.SW.S; s++)   
		{
			if(alea(0,1)<0.5)
				R.SW.X[s]=hammersley(pb.SS,R.SW.S, deb,s);
			else
				for (d = 0; d < pb.SS.D; d++)  
				{  
					R.SW.X[s].x[d] = alea (pb.SS.min[d], pb.SS.max[d]);
				}
		}
		break;

		case 4: // 	Random/Regular	
		for (s = 0; s < R.SW.S; s++)   
		{	
			for (d = 0; d < pb.SS.D; d++)  
			{  
				R.SW.X[s].x[d] = floor(1+(pb.nb+1)*alea (pb.SS.min[d], pb.SS.max[d]))/(pb.nb+1);
			}
		}
		break;

		case 5: // Centre biased 2
		for (s = 0; s < R.SW.S; s++)   
		{
			for (d = 0; d < pb.SS.D; d++)  
			{  
				z=0.5*(pb.SS.min[d]+pb.SS.max[d]); // centre
				zz=alea (0, 1)-0.5; // Normalised distance to the centre
				zzz=(pb.SS.max[d]-pb.SS.min[d])*pow(fabs(zz),1+1./R.SW.S); // decrease it

				if(zz>0)
				 R.SW.X[s].x[d]=z+zzz; // re-add it
				else
				 R.SW.X[s].x[d]=z-zzz;
			}
		}
		break;
	}


	// Velocity
// 0 => on each dimension U(xmin,xmax) -U(xminxmax)
		// 1 => on each dimension U(xmin,xmax) - X(d)
		// 2 =>  Y - X, where Y is the initial position of another particle, chosen at random
		// 3 => set to zero
	switch(param.initVel)
	{
		case 0:
			for (s = 0; s < R.SW.S; s++) // 1, 0.5 
			{	
				for (d = 0; d < pb.SS.D; d++)  
				{				
					R.SW.V[s].v[d] =
					alea (pb.SS.min[d], pb.SS.max[d])-alea (pb.SS.min[d], pb.SS.max[d]);
				}
			}
		break;

		case 1:
		case1:
			for (s = 0; s < R.SW.S; s++)   // xi - xj
			{	
				s1=alea_integer(0,R.SW.S-1);
				for (d = 0; d < pb.SS.D; d++)  
				{  
					R.SW.V[s].v[d] =alea(pb.SS.min[d], pb.SS.max[d])-R.SW.X[s].x[d];
				}
			}
		break;

		case 2:
		case2:
			for (s = 0; s < R.SW.S; s++)   // xi - xj
			{	
				s1=alea_integer(0,R.SW.S-1);
				for (d = 0; d < pb.SS.D; d++)  
				{  
					R.SW.V[s].v[d] =R.SW.X[s1].x[d]-R.SW.X[s].x[d];
				}
			}
		break;

		case 3:
			for (s = 0; s < R.SW.S; s++)  // 0  
			{	
				for (d = 0; d < pb.SS.D; d++)  
				{				
					R.SW.V[s].v[d] = 0;
				}
			}
		break;

		case 4:
		if(alea(0,1)<0.5) goto case1; else goto case2;

		case 5: // Half-diff
		for (s = 0; s < R.SW.S; s++) // 1, 0.5 
			{	
				for (d = 0; d < pb.SS.D; d++)  
				{				
					R.SW.V[s].v[d] =
					0.5*(alea (pb.SS.min[d], pb.SS.max[d])-alea (pb.SS.min[d], pb.SS.max[d]));
				}
			}
		break;
	}


	// Sizes
	for (s = 0; s < R.SW.S; s++)   
	{
		R.SW.X[s].size = pb.SS.D;
		R.SW.V[s].size = pb.SS.D;
//printf("\n--");
//for(d=0;d<pb.SS.D;d++) printf("\n%f ",R.SW.X[s].x[d]);
	}
		
	if(param.init==-1) fclose (f_init);

	// Distribution method(s)
	switch(param.R)
	{
		default:
		for (s = 0; s < R.SW.S; s++)  R.SW.R[s] =param.R;
		break;

		case 3:
 		for (s = 0; s < R.SW.S; s++)
		{
			if(alea(0,1)<0.5 ) R.SW.R[s]=0; // Uniform
			else R.SW.R[s]=2; //  Decreasing function of the 1D distance
		}
		break;
	}

	// Strategies (different w)
	switch(param.strat)
	{
		default:
		for (s = 0; s < R.SW.S; s++)   
		{
				R.SW.w[s]=param.w[0];
				R.SW.c[s]=param.c1;
		}
		break;
		case 1: // Random w (uniform) for some particles
		for (s = 0; s < R.SW.S; s++)   
		{
			if(alea(0,1)<param.w0Rate)
				 R.SW.w[s]=param.w[0]; 
			else R.SW.w[s]=param.w[1];

			R.SW.c[s]=param.c1;
		}
		break;

		case 2: // Random w (non-uniform) for some particles