pso_binary.c

C语言编写的基本微粒群算法的应用介绍源码

  goto end;

//----------------------------------------------- OPTIMISATION

optim:
t1=clock();
mean_val=pow(2.0,D-1);

// Save if multiobjective
if (n_f>1)
	f_multi=fopen("f_multi.txt","w");

alpha=pow(0.99,1.0/(3*D-1)); // This is just to estimate the "quality" of the algorithm.
							// Better to run on several eval_max.
							// For the "ideal" algo, we want a succes rate of 0.99
							// after at most 3 evaluations for each dimension

quality_ref=0;

//------------------------------LOOP on maximum number of evaluations
for (eval_max=eval_max1;eval_max<=eval_max2 ;eval_max=eval_max+delta_eval_max) 
{

// Initialisation of information variables
n_exec=0; eval_mean=0; eps_mean=0; n_failure=0;
save_time=0;

// x_p = best previous position
// x_g = best previous position of the best informer
cp1=1.;
if(deter==1) {cp2=1; cp3=1;} // Deterministic: always towards x_p and x_g

init:          // Initialisations of positions and velocities
init=1;
n_iter=-1;
n2=1; n3=1; // For method 7

if (adapt_method>0) // If adaptive method
{
	n_adapt=0; method_rank=0;
	method=strategies[0];
	adapt_latency=1;
	n_iter0=n_adapt_max;
	n_iter1=0;
	n_iter1=n_adapt_max;
// if(method==100) K=35; else K=alea_integer(2,3);   //TEST
}


if (adapt_K==1) // If adaptive K
{
	K=1+2*alea_integer(0,1); // 1 or 3
	K_choice[0][0]=0; K_choice[2][0]=0; // Cumulated improvement
	K_choice[0][1]=0; K_choice[2][1]=0; // Number of times K is used
	n_adaptK=0;
}

// Some confidence coefficients for some methods
pr0=0.5+0.5/log((1+log((double)D))); // a
//pr0=1/pow(2,D); // b
c_pr0=0.5*log((double)D)/log((double)S);

n_exec=n_exec+1;


  // Initialise particles and velocities
  init_swarm(init_option);

// First evaluations
nb_eval=0;

for (s=0;s<S;s++)
{
  P[s].f.size=n_f;
	for (i=0;i<n_f;i++)
		P[s].f.fi[i]=fabs(perf(s,function[i])-fmin.fi[i]);
	P_m[s]=P[s]; // Best position = current one
}



// Save the best of the bests
best=P_m[0];
for (s=0;s<S;s++) if (best_(P_m[s].f,best.f)==1) best=P_m[s];


// Current min error
error=total_error(best.f);
error_prev=error;
error_prev2=error_prev;
init_links=1; // So that information links will be initialized

// For method 7
cp0=1-0.5*log((double)D)/log((double)S);
if (K>1) cp0=cp0*log((double)K);
 //z7=(double)nb_eval*K/(double)S;
  z7=(double)K/(double)S;
 if (z7>1) z7=log(z7);

 cp1=cp0/z7;
 cp2=cp1;

 if(K>2)
  cp3=cp1/log(K-1);// Decreasing. 
 else cp3=cp1;

   n2=(int)(cp2*D);// Decreasing. The more you trust p, the less you update it
   n3=(int)(cp3*D);

   //Save the state
		for (s=0;s<S;s++) {P_prev[s]=P[s]; P_m_prev[s]=P_m[s];V_prev[s]=V[s];}
		best_prev=best;

 // Modified bits. Useful for some methods
for (s=0;s<S;s++)
{
  B_mod[s].size=D;
  for (d=0;d<D;d++) B_mod[s].b[d]=0; // Not yet modified
}
   
//-------------------------------------- ITERATIONS
init=0;
loop:
Diam=diameter();
//fprintf(f_trace,"\n%i %i",nb_eval,Diam);
//fprintf(f_trace,"\n%i %i",nb_eval,nb_pivots);
n_iter=n_iter+1;

 if(method==50)
 {



    n_adapt=n_adapt+1;
    if(n_adapt>n_adapt_max)
    {
      n_adapt=0;
      if (error>=error0) k0=k0+1; else {k0=k0-1;if(k0<1) k0=1;}
 //printf("\n %i %f %f",k0, error0,error);     
      error0=error;
    } 
  }


 if (adapt_method==1) 
 {
  if (error>=error_prev)// If there has been no improvement ...
  {

	 n_adapt=n_adapt+1;
	if(n_adapt>n_adapt_max) // ... and if 
	{
    method_rank=(method_rank+1)%strat_number;
	//method_rank=alea_integer(0,strat_number-1); // Test
		method=strategies[method_rank];
		n_adapt=0;
	}
  }
  else n_adapt=0; // TEST
 }


if (adapt_method==2) // Try TWO strategies starting from the same state
					// and keep the best one
{
//fprintf(f_trace,"\n%i %i %f",n_iter,method, error);	
	if(n_iter==n_iter1+adapt_latency*n_adapt_max)
	{
		n_iter0=n_iter;
		error1=error;
   //Save the state
		for (s=0;s<S;s++) {P_prev[s]=P[s]; P_m_prev[s]=P_m[s];V_prev[s]=V[s];}
    //*** (TO DO: add saving for B_mod if "taboo" is use, and restore below)
		best_prev=best;
//fprintf(f_trace,"\n save the state");
	}
	if(n_iter==n_iter0+n_adapt_max)
	{
		error2=error; // Error after iter n_iter-1
//fprintf(f_trace,"\niter %i method %i error2 %f",n_iter,method,error2);

     // Save the final state after having used method rank 0
		for (s=0;s<S;s++) {P_temp[s]=P[s]; P_m_temp[s]=P_m[s];V_temp[s]=V[s];}
		best_temp=best;

		for (s=0;s<S;s++) {P[s]=P_prev[s];P_m[s]=P_m_prev[s];V[s]=V_prev[s];} // Back

		best=best_prev;

		method_rank=(method_rank+1)%2; // Switch method
		method=strategies[method_rank];
		init_links=0; // Keep the same links
	}
	if(n_iter==n_iter0+2*n_adapt_max) 
	{
		error3=error; // Error after back and having tried method rank 0
//fprintf(f_trace,"\niter %i method %i error3 %f",n_iter,method,error3);
		if(error3>error2) // If no improvement
		{			
			method_rank=(method_rank+1)%2; // Switch (i.e. back to the previous method)
			method=strategies[method_rank];
			init_links=0; // Keep the same links

			for (s=0;s<S;s++) {P[s]=P_temp[s];P_m[s]=P_m_temp[s];V[s]=V_temp[s];} // Back
			best=best_temp;
		}

//printf("\n error2 %f error3 %f n_iter %i method %i",error2,error3,n_iter,method);
		n_iter1=n_iter;
		adapt_latency=adapt_latency+1;
	}
 // if(method==100) K=35; else K=alea_integer(2,3);  //TEST
}

//printf("\n method %i for iter %i, error %f",method,n_iter, error);

if(adapt_K==1)

{
	K_choice[K-1][1]=K_choice[K-1][1]+1; // Number of times K has been used

	//cp1=(error_prev2-error_prev)/error_prev2;



	cp2=(error_prev-error) /error_prev;
	//cp2=cp2*cp2;
	cp2=log(1+cp2);
	K_choice[K-1][0]=K_choice[K-1][0]+cp2; // Cumulated improvement

//if (cp2>0) printf("\n %f", cp2);
	n_adaptK=n_adaptK+1; // Let some time to the current K

	if(n_adaptK>n_adapt_max)
	{
		if (K_choice[0][1]*K_choice[2][1]==0)
		{
			K=4-K;// Switch. 
		}
		else
		{
			cp0=K_choice[0][0]/K_choice[0][1];
			cp1=K_choice[2][0]/K_choice[2][1];
			pr=alea(0, cp0+cp1);
			if (pr<cp0) K=1; else K=3;
//printf("\n %.0f %.0f %f %f pr %f =>K %i", K_choice[0][1],K_choice[2][1],cp0,cp1,pr,K);
//printf("\n %i",K);
		}
		n_adaptK=0;
	}
}



if(init_links==1 && !fixed_link) // Who informs who, at random
{
  for (s=0;s<S;s++) //
	{
		for (m=0;m<S;m++) LINKS[m][s]=0;
		LINKS[s][s]=1; // However, each particle informs itself
	}

  // Other links
	  for (m=0;m<S;m++)
	  {
		s1=0;s2=S-1;

		if (K>1)
			for (k=0;k<K;k++)
			{
			s=alea_integer(s1,s2);
			LINKS[m][s]=1;
		//	LINKS[s][m]=1;
			}
	  }
}
if (fixed_link && nb_eval<2*S) // Constant circular neighbourhood
{
	k1=(K-1)/2;
	if(K%2==0) k2=K/2; else k2=k1;
	for (s=0;s<S;s++) //
	{
		for (m=0;m<S;m++) LINKS[m][s]=0;
		for (m=s;m<=s+k2;m++) 
		{i=m%S; LINKS[i][s]=1;}
		for (m=s;m>=s-k1;m=m-1) 
		{i=(S+m)%S; LINKS[i][s]=1;}
	}
}

if(print_level>1) // Display links

{
   printf("\nLINKS");
	for (s=0;s<S;s++) //
	{ printf("\n");
		for (m=0;m<S;m=m++) printf("%i ",LINKS[m][s]);
	}
}


// The swarm MOVES

for (s=0;s<S;s++) pivots[s]=0;

for (s=0;s<S;s++)  // Best informants
{
	g=best_info(s,option_best);
  pivots[g]=1;
    // Note. TO DO. Memorise g in a table and use it below
 }
 // Number of pivots . The mean number is S/(K-1)
 nb_pivots=0; for (s=0;s<S;s++) nb_pivots=nb_pivots+ pivots[s];
 // fprintf(f_trace,"\n %i",nb_pivots);
  
for (s=0;s<S;s++)  // For each particle ...
{  
	g=best_info(s,option_best);
 
switch (method)
{
	
	case 0: // Magic PSO
 /*
	if (S==1) 
		cp2=alea(-1,1);// cp2=alea(0,1); // Test for S=1
	else
	{  
		if(deter==0) // Random choices in pure binary algo
		{  
			i=alea_integer(0,2);
			//i=alea_integer(0,3);
		if (i==0){cp2=1;cp3=1;}   // Towards x_p and x_g
		if (i==1) {cp2=1;cp3=-1;} // Towards x_p away from x_g
		if (i==2) {cp2=-1;cp3=1;} // Away from x_p and towards x_g
		if (i==3) {cp2=-1;cp3=-1;} // Away from x_p and away from x_g
		}

	} 
 */
  cp2=2*alea_integer(0,1)-1;   cp3=2*alea_integer(0,1)-1;
	for (d=0;d<D;d++)
	{
		bit1=P[s].x.b[d]; bit2=P_m[s].x.b[d];

		if (S==1) 
			V[s].b[d]=(int)(V[s].b[d] + 2*cp2*(bit2-bit1)); // Test for S=1
		else
		{
			bit3=P_m[g].x.b[d];

			if(deter==-1) // more randomness in pure binary algo

			{
			i=alea_integer(0,2);
			// i=alea_integer(0,3);
			if (i==0){cp2=1;cp3=1;}   // Towards x_p and x_g
			if (i==1) {cp2=1;cp3=-1;} // Towards x_p away from x_g
			if (i==2) {cp2=-1;cp3=1;} // Away from x_p and towards x_g
			if (i==3) {cp2=-1;cp3=-1;} // Away from x_p and away from x_g
			}
 /*
   if ( V[s].b[d]<-1) // Useful only when using several strategies
   {
    V[s].b[d]=-1; 
    }
    else if (V[s].b[d]>1) V[s].b[d]=1; else V[s].b[d]=0;
 */   
		V[s].b[d]=(int)(V[s].b[d] + cp2*(bit2-bit1) +cp3*(bit3-bit1)); // in [-3...3]
		P[s].x.b[d]=P[s].x.b[d]+V[s].b[d];  // in {-3...4}     
		P[s].x.b[d]=(4+(int)P[s].x.b[d])%2;  // new position in {0,1}
//		V[s].b[d]= V[s].b[d]-bit1+  P[s].x.b[d];
 //printf("\n %f",   V[s].b[d]);
		V[s].b[d]=(int)(3+V[s].b[d])%3-1; // new velocity in {-1,0,1}
 /*       // TEST
    if (V[s].b[d]<-1 || V[s].b[d] >1)
    {printf("\nERROR 0: %f",  V[s].b[d]);
      printf("\n %i %i %i",bit1,bit2,bit3);
      printf("\n %f %f",cp2,cp3);
  scanf("%i",&i);
    }
 */
		}
	}
	break;

case 1:

   // Method 1 -----------------------------------
cp0=0.5;  // Random move, relative probability
if (deter!=-1)
{  
  cp1=cp0;

  cp2=(double)nb_eval*K/(double)S;
  cp2=1-cp1/log(cp2);
  cp3=cp2;
  cp=cp1+cp2+cp3;
}

	for (d=0;d<D;d++)

	{
    if (deter==-1)
    {

      cp1=alea(cp0/2,cp0);
      cp2=(double)nb_eval*K/(double)S;
      cp2=1-cp1/log(cp2);
      cp3=cp2;
      cp=cp1+cp2+cp3;
    }
    bit1=alea_integer(0,1);
    bit2=P_m[s].x.b[d];   bit3=P_m[g].x.b[d];
    
    pr=alea(0,cp);
    if (pr<=cp1) bit= bit1;

    else 
      if (pr<=cp1+cp2) bit=bit2; else bit=bit3;


     P[s].x.b[d]=bit; 
    //V[s].b[d]=  P[s].x.b[d]-bit0; // -1 or 0 or 1. 
	}
break;


case 2:
// Method 2 ---------------------------
cp0=2/3.;

cp0=0.8;
 cp1=cp0;

  cp2=(double)nb_eval*K/(double)S;
  cp2=1-cp1/log(cp2);

  cp3=cp2;

	for (d=0;d<D;d++)
	{
    if (deter==-1)
    {
      cp1=alea(cp0/2,cp0);
      cp2=(double)nb_eval*K/(double)S;
    cp2=1-cp1/log(cp2);
    cp3=cp2;
    }
    
    bit1=alea_integer(0,1);  bit2=P_m[s].x.b[d];   bit3=P_m[g].x.b[d];
    pr=alea(0,1); if (pr<cp1) choice[0]= bit1; else choice[0]=1-bit1;
    pr=alea(0,1); if (pr<cp2) choice[1]=bit2; else choice[1]=1-bit2;
    pr=alea(0,1); if (pr<cp3) choice[2]=bit3; else choice[2]=1-bit3;

    //Majority
    bit=choice[0]+choice[1]+choice[2];
    if(bit<2)  P[s].x.b[d]=0; else P[s].x.b[d]=1;

	}
break;

case 3:
// Method 3 -----------------------
	cp1=0.5;
	cp=(double)nb_eval*K/(double)S;
	cp=cp1/log(cp);


	for (d=0;d<D;d++)
	{
		bit0=0; bit1=0;

		for (s1=0;s1<S;s1++) // Ask all informers
		{
			if(LINKS[s1][s]==0) continue;
				bit=P_m[s1].x.b[d];

				if (bit==0) bit0=bit0+1; else bit1=bit1+1;		
		}

		pr=alea(0,bit0+bit1);
		if(pr<=bit0) bit=0; else bit=1; // Probabilistic majority...

		pr=alea(0,1); if(pr<=cp) bit=1-bit;// ... not sure 
		P[s].x.b[d]=bit;
	}

break;

case 4:
// Method 4 ------------------------------------------
  cp1=0.5;
  cp2=(double)nb_eval*K/(double)S;
  cp2=1-cp1/log(cp2);
  cp=cp1+K*cp2;

	for (d=0;d<D;d++)
	{
		bit0=0; bit1=0;
		for (s1=0;s1<S;s1++) // Ask all informers
		{
			if(LINKS[s1][s]==0) continue;

				bit=P_m[s1].x.b[d];
				if (bit==0) bit0=bit0+1; else bit1=bit1+1;
		}
     // Weighted majority
		if(bit0>bit1) bit=0;else if (bit1>bit0) bit=1; else bit=alea_integer(0,1);

    pr=alea(0,cp); if(pr<cp1) bit=alea_integer(0,1);
		P[s].x.b[d]=bit;
	}
break;


case 5:
// Method 5---------------------------------

  //cp1: Confidence coefficient for velocity
  //cp2: Confidence coefficient for p
  //cp3: Confidence coefficient for g

cp0=1-0.5*log((double)D)/log((double)S);
 if (K>1) cp0=cp0*log((double)K);


 z=(double)nb_eval*K/(double)S;
 if (z>1) z=log(z);