pso_binary.c

采用二进制对PSO进行编码

 int i,max,rank,z1,z2;
 struct bits bit;

 for (i=0;i<size;i++) bit.b[i]=0;

 bit.size=size;
 if(a==0) return bit;
 if(a==1) {bit.b[0]=1; return bit; }
 
 z1=a;
 rank=0;
 max=1;
 loop:
 z2=z1/2;
 if(z2==0) {bit.b[rank]=1; goto end;}

if (rank<=size-1)
  {
    if(2*z2<z1) bit.b[rank]=1;
    z1=z2; rank=rank+1; goto loop;
  }

end:
// for (i=0;i<size;i++) printf("%1i",bit.b[size-i-1]);

return bit;
}

//===========================================================
double combin(int Dim,int k)
{
   double c;
   int d;
  c=Dim;
  for (d=1;d<=k-1;d++) c=c*(Dim-d)/(double)(d+1);
  return c;
}
//===========================================================
int diameter()   // Swarm diameter
{
 int diam,dist;
 int s1,s2;
 diam=0;
 for (s1=0;s1<S-1; s1++)
 for (s2=s1+1;s2<S;s2++)
 {
   // dist=distance(s1,s2,-1); // Memory swarm
   dist=(int)distance(s1,s2,-2); // Current swarm

   if(dist>diam) diam=dist;
 }
  return diam;
}
//===========================================================
double	distance(int s1,int s2, int option)
{
// Distance between two positions
// option = 0 => s1 and s2 are ranks of best memorized positions
// option = 1 => s1 is the rank of a best memorized position
//               and s2 the rank of a current position
// option <0 => Hamming distance
// WARNING. Depending on your computer, it may not work
// for high dimension
  int	d;
  double dist;
  int n1;
  
  if(option<0) goto hamming;

  dist=0;
  for (d=0;d<P[s1].x.size;d++)
  {
    if(option==0) n1=abs(P_m[s1].x.b[d]- P_m[s2].x.b[d]);
    else n1=abs(P_m[s1].x.b[d]- P[s2].x.b[d]);
    if (n1>0)
    {
      dist=dist+pow(2.0,d);

    }
  }
  if (dist>mean_val) dist=2*mean_val-dist; // Cyclic distance
  return dist;

hamming:
 switch(option)

 {
case -3:
  dist=0;
  for (d=0;d<P_m[s1].x.size;d++)
  {

    n1=abs(P[s1].x.b[d]- P_m[s2].x.b[d]);
	dist=dist+n1;
  }

  return dist;


 default: // Between  current positions
  dist=0;
  for (d=0;d<P[s1].x.size;d++)
  {
    n1=abs(P[s1].x.b[d]- P[s2].x.b[d]);
	dist=dist+n1;
  }
  return dist;

 case -1:
  dist=0;
  for (d=0;d<P_m[s1].x.size;d++)
  {
    n1=abs(P_m[s1].x.b[d]- P_m[s2].x.b[d]);
	dist=dist+n1;
  }
  return dist;
 }

}
//===========================================================
void  fundamental_hyp(int function)
{  // Statistically checks the fundamental hypothesis
 // "Nearer is better"
 // NOTE: this program is in fact "independent". No optimisation

   int  bit;
   int d;
  double fx,fy;
  double f[3];
  int m,n;
//  double dis[3];
  double dist,distx,disty;
//  int loop;
 double  n_sample;
 int nn;
 int option;
// int r;
 int radius;
 float rate_sample;
 double rho;
 int s0,s1;
 double  true;
 double v;

//option: 0 => use distance to the solution
           // (check global improvement)
   // 1 => check local improvement
 

 printf("\n Global 0, local 1, local 2, local 3 ?: ");
 scanf("%i", &option);
 
 printf("\n Sampling rate? ");
 scanf("%f", &rate_sample);

  v=pow(pow(2,D),2);
  n_sample=v*rate_sample;
 //printf("\nDimension is %i. Radius? ",D);
 //scanf("%i",&radius);
  P[0].x.size=D;  P[1].x.size=D;P[2].x.size=D;
    true=0;
switch(option)
{

    case 0:  // For global search
  switch(function) // Solutions are supposed to be known
  {
    case 1:
    for (d=0;d<D;d++) P[0].x.b[d]=1;
    break;

    case 2:
    // Note that there are in fact  pow(D/12,2) solutions

    for(n=0;n<=D/12;n=n+1)
    {
      bit=alea_integer(0,1);

      m=12*n;
      for (d=m;d<m+6;d++) P[0].x.b[d]=bit;
      for (d=m+6;d<m+12;d++) P[0].x.b[d]=1-bit;
    }
    break;

  case 3:
  for (d=0;d<D;d++) P[0].x.b[d]=0;
  break;

  case 6:
  for (d=0;d<D;d++) P[0].x.b[d]=0;
  for(m=3;m<=D;m=m+6)
  {
    for (d=m;d<m+3;d++) P[0].x.b[d]=1;
  }
  break;
  

  case 7:
  for (d=0;d<D-1;d++) P[0].x.b[d]=0;
   P[0].x.b[D-1]=1;
  break;

  case 9:
  for (d=0;d<D;d++) P[0].x.b[d]=0;
  break;

  case -1:
  // Solution for Q12: 000010101001
  for (d=0;d<D;d++) P[0].x.b[d]=0;
   P[0].x.b[4]=1;
   P[0].x.b[6]=1;
   P[0].x.b[8]=1;
   P[0].x.b[11]=1; 
  break;

  case 8:
  // Solution for f_landscape_12_20: 011111111110
    for (d=0;d<D;d++) P[0].x.b[d]=1;
   P[0].x.b[0]=0;
   P[0].x.b[11]=0;
  break;
  
  case 99:
  for (d=0;d<D;d++) P[0].x.b[d]=0;
  break;
  
  default:

  printf("\n NOT YET for function %i, SORRY",function) ;
  goto end;
 }  // end switch function
 //---------------
   printf("\n Solution\n");
 for(d=0;d<D;d++) printf("%i", P[0].x.b[d]); printf("\n");
       // Compute "solution" fitness
      f[0]=perf(0,function);

 for (radius=2;radius<=D;radius++)
 {
    nn=0;
 for (n=0;n<n_sample;n++)
 { 

      //Define at random point x
      for (d=0;d<D;d++) P[1].x.b[d]=alea_integer(0,1);
 
      // Define at random y at distance <=radius

 // for (d=0;d<D;d++) P[2].x.b[d]=alea_integer(0,1);
 
      P[2]=P[1];
      dist=alea_integer(1,radius); 
      for (m=0;m<dist;m++)
      {
          d=alea_integer(0,D-1);
          P[2].x.b[d]=1- P[1].x.b[d];
      }

      // Check P2#P1
      for (d=0;d<D;d++) if(P[2].x.b[d]!=  P[1].x.b[d]) goto global3;
      continue; // They are the same
   global3:

   nn=nn+1;  
      // Compute dist(x,solution)
       distx=distance(1,0,-2);
      // Compute dist(y,solution)
      disty=distance(2,0,-2);
          
      // Compute x fitness
      f[1]=perf(1,function);
      // Compute y fitness
      f[2]=perf(2,function);

      // Check the hypothesis
 
      fx=fabs(f[1]-f[0]); fy=fabs(f[2]-f[0]);

  if((distx-disty)*(fx-fy)>0)
  {
 //printf("\n %f %f %f %f",distx,disty,fx,fy);
    true=true+1;continue;
    }
 // if(distx==disty) true=true+0.5;
  }
   true=0.5*true/nn;

 printf("\n%i %f",radius,true);
 fprintf(f_trace,"%i %i %f %i %f \n", function, D, rate_sample, radius, true);
  }
  break;
  
 case 1: // For local search
 for (radius=2;radius<=D;radius++)
 {
    nn=0;
 for (n=0;n<n_sample;n++)
 {
      for (d=0;d<D;d++) P[0].x.b[d]=alea_integer(0,1); // Random point 0
      // Random point 1
      P[1]=P[0];
      dist=alea_integer(1,radius);
      for (m=0;m<dist;m++)

      {
          d=alea_integer(0,D-1);
          P[1].x.b[d]=1- P[0].x.b[d];
      }
      // Random point 2
      P[2]=P[0];
      dist=alea_integer(1,radius);
      for (m=0;m<dist;m++)
      {

          d=alea_integer(0,D-1);
          P[2].x.b[d]=1- P[0].x.b[d];
      }

      // Check P2#P1
      for (d=0;d<D;d++) if(P[2].x.b[d]!=  P[1].x.b[d]) goto local3;
      continue; // They are the same

      local3:
      nn=nn+1;
      for(m=0;m<3;m++) // Evaluations
        f[m]=perf(m,function);

      distx=distance(1,0,-2);
      disty=distance(2,0,-2);

      fx=fabs(f[1]-f[0]); fy=fabs(f[2]-f[0]);

    if((distx-disty)*(fx-fy)>0) {true=true+1;continue;}
   // if(distx==disty) true=true+0.5;
  }
   true=true/nn;
 printf("\n%i %f",radius,true);
 fprintf(f_trace,"%i %i %f %i %f \n", function, D, rate_sample, radius, true);
  }

 break;
 case 2:   // For local search. Derivation 11 simulation
  for (radius=2;radius<=D;radius++)

 {
    nn=0;
 for (n=0;n<n_sample;n++)
 {
      for (d=0;d<D;d++) P[0].x.b[d]=alea_integer(0,1); // Random point 0
      // Random point 1
      P[1]=P[0];
    //  dist=alea_integer(1,radius);
    //  for (m=0;m<dist;m++)

    rho=alea(1,radius);
    for(m=0;m<rho;m++)
      
      {
          d=alea_integer(0,D-1);
          P[1].x.b[d]=1- P[0].x.b[d];
      }

      for (m=0;m<2;m++) f[m]=perf(m,function); // Evaluations
       s0=0;s1=1;
       if (f[1]>f[0]) {s0=1;s1=0;}   // s1 is the best  (seen as "g")
  
      // Random point 2 around the best
      P[2]=P[s1];
      dist=alea_integer(1,radius);
      for (m=0;m<dist;m++)
      {
          d=alea_integer(0,D-1);
          P[2].x.b[d]=1- P[s1].x.b[d];
      }

      nn=nn+1;

    f[2]=perf(2,function);
    // If the new point is better ...
    if(f[2]<f[s0]) {true=true+1; continue;}
  //  if(f[2]==f[s0]) {true=true+0.5; continue;}
  }
   true=true/nn;
 printf("\n%i %f",radius,true);
 fprintf(f_trace,"%i %i %f %i %f \n", function, D, rate_sample, radius, true);
  }

 break;

 case 3:  // Local search just around the current position
   for (radius=2;radius<=D;radius++)
 {
    nn=0;
 for (n=0;n<n_sample;n++)
 {
      for (d=0;d<D;d++) P[0].x.b[d]=alea_integer(0,1); // Random point 0
      // Random point 1
      P[1]=P[0];
      dist=alea_integer(1,radius);
      for (m=0;m<dist;m++)
      {
          d=alea_integer(0,D-1);
          P[1].x.b[d]=1- P[0].x.b[d];
      }

      for (m=0;m<2;m++) f[m]=perf(m,function); // Evaluations
 
      nn=nn+1;

      
    // If the new point is better ...
    if(f[1]<f[0]) {true=true+1; continue;}
    //if(f[2]==f[s0]) {true=true+0.5; continue;}
  }
   true=true/nn;
 printf("\n%i %f",radius,true);
 fprintf(f_trace,"%i %i %f %i %f \n", function, D, rate_sample, radius, true);
  }

 break;
 } // end switch option
 
  end:;
}
//===========================================================
void init_swarm(int option)
{  /*
 option 0 => positions completely at random
 option 1 => (removed. Was not good)
 option <0 => read from file init.txt
*/
 double CDk[S_max];
int check;
int d;
double dd;
int dist0[S_max],dist1[S_max];
int Dist[S_max][S_max];
//int dmin[1000];// Just for information
int H[D_max]; // Just for information about distances
int k,kk,kkk;
int k1,k2;
double mean0,mean1;
double max;
int min;
double min_r;
int Nk[D_max];
//int norm;
double r;
double R[S_max][D_max];
int St;
int s,s1,s2;
struct position sample[2000];
int sample_s[2000];
int sample_n[S_max];
double sample_d[2000];
double sample_R[2000][D_max];
int space;
double tot;
 double var0,var1;
 
 for (s=0;s<S;s++) P[s].x.size=D;

switch(option)
{
  default: // Completely at random
for (s=0;s<S;s++)
{
	for (d=0;d<D;d++) P[s].x.b[d]=alea_integer(0,1);
}
break;
 


case 2: // Regular distribution

for(d=0;d<D;d++) P[0].x.b[d]=0;  // First position 000...000
//for(d=0;d<D;d++) P[1].x.b[d]=1;  // Second position 111...111
St=0;
loop2:
St=St+1;   // New position rank
for (d=0;d<D;d++)
{
   // Try 0
	P[St].x.b[d]=0;
  for (s=0;s<St;s++)
  {
    dist0[s]=0;
	for (k=0;k<=d;k++)
	{
	  if(P[St].x.b[k] != P[s].x.b[k]) dist0[s]=dist0[s]+1;
	}
   }
   // Mean
    mean0=0;
     for (s=0;s<St;s++) mean0=mean0+dist0[s];
    mean0=mean0/St;
   
   // Variance
    var0=0;
      for (s=0;s<St;s++) var0=var0+ (dist0[s]-mean0)*(dist0[s]-mean0);
   
   // Try 1
	P[St].x.b[d]=1;
   for (s=0;s<St;s++)
  {
 	dist1[s]=0;
	for (k=0;k<=d;k++)
	{
	  if(P[St].x.b[k]!= P[s].x.b[k]) dist1[s]=dist1[s]+1;
	}
   }

   // Mean
   mean1=0;
     for (s=0;s<St;s++) mean1=mean1+dist1[s];
    mean1=mean1/St;
   // Variance
   var1=0;
      for (s=0;s<St;s++) var1=var1+ (dist1[s]-mean1)*(dist1[s]-mean1);


   // Choose the best
	  if(mean1>mean0) {P[St].x.b[d]=1; continue;}
	  if(mean1<mean0) {P[St].x.b[d]=0; continue;}

		if(var1<var0) P[St].x.b[d]=1;
		else P[St].x.b[d]=0;

//printf("\nSt %i, %f %f, %i",St,var0,var1,P[St].x.b[