📄 pso_binary.c
字号:
break;
case 17: // Like 16, but with condition "the particle has just improved its position"
Gp=P_m[g];
dist=(int)log(D); // We suppose here D>=2
r=alea(1,dist);
error1=total_error(P[s].f);
error2=total_error(P_m[s].f);
i=error1<=error2; // Improvement test
for (k=0;k<r;k++)
{
d=alea_integer(0,D-1);
if((i && B_mod[s].b[d]<1) || !i)
{
Gp.x.b[d]=1- Gp.x.b[d]; // Around g
}
}
for (d=0;d<D;d++)
{
if (Gp.x.b[d]!=P[s].x.b[d]) B_mod[s].b[d]=1; // Has been modified
else B_mod[s].b[d]=0; // Has not been modified
P[s].x.b[d] =Gp.x.b[d];
}
break;
case 19: // Like 11, but with "taboo"
//dist=log(D); // We suppose here D>=2
dist=(int)(1+Diam/(double)nb_pivots); if (dist<3) dist=3;
r=alea(1,dist);
P[s]=P_m[g];
error1=total_error(P[g].f);
error2=total_error(P_m[g].f);
i=error1<=error2; // g has just improved its position
for (k=0;k<r;k++) // Around g
{
d=alea_integer(0,D-1);
if(!i || (i && B_mod[s].b[d]<1)) // "Taboo"
{
P[s].x.b[d]=1- P[s].x.b[d];
B_mod[s].b[d]=1;
}
else
{
B_mod[s].b[d]=0;
}
}
break;
case 99: // TEST . 7 with taboo
if (error<error_prev) goto skip99;
z7=z7+(double)K/(double)S;
z=log(z7);
cp1=cp0/z;
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; //Decreasing. The more you trust g, the less you update it
skip99:
Pp=P_m[s];Gp=P_m[g];
for (k=0;k<n2;k++) // Around p. At most n2 modified components
{
d=alea_integer(0,D-1); Pp.x.b[d]=1- Pp.x.b[d];
}
for (k=0;k<n3;k++) // Around g. At most n3 modified components
{
d=alea_integer(0,D-1); Gp.x.b[d]=1- Gp.x.b[d];
}
for (d=0;d<D;d++) // Choose each component
{
if(B_mod[s].b[d]<1) // Taboo
{
bit1=P[s].x.b[d];
bit2=Pp.x.b[d]; bit3=Gp.x.b[d];
//Majority choice
if(bit2==bit3) P[s].x.b[d]=bit2;
else P[s].x.b[d]=alea_integer(0,1);
if (P[s].x.b[d]!=bit1) B_mod[s].b[d]=1;
}
else B_mod[s].b[d]=0;
}
break;
//=========== Methods a priori designed for constant circular neighbourhood
// But, of course you may try them with random link (cf fixed_link parameter)
case 100: // KE-C Method (Jim Kennedy and Russ Eberhart method, improved by Maurice Clerc)
cp0=1;
cp1=2.0;
cp2=cp1;
for (d=0;d<D;d++)
{
cp3=cp0*V[s].b[d] + alea(0,cp1)*(P_m[s].x.b[d]-P[s].x.b[d]); // cp0*v + rand(0,cp1)*(p-x)
cp3=cp3+alea(0,cp2)*(P_m[g].x.b[d]-P[s].x.b[d]); // + rand(0,cp2)*(g-x)
//Note that g is either the local best or the global best
// depending on the neighbourhood size
V[s].b[d]=(int)cp3; // New integer velocity
cp3= P[s].x.b[d]+ V[s].b[d]; // Add velocity
cp3=1/(1+exp(-cp3)); // Keep position in ]0,1[
pr=alea(0,1); // Probabilistic binary choice
if (pr<cp3) P[s].x.b[d]=1; else P[s].x.b[d]=0;
}
break;
case 101: // TL Method (Tasgetiren and Liang. 2004)
// Suggested S=2*D, K=S
cp0=4.0; // Max velocity (absolute value)
cp1=2.0;
cp2=cp1;
for (d=0;d<D;d++)
{
cp3=V[s].b[d]+ alea(0,cp1)*(P_m[s].x.b[d]-P[s].x.b[d]); // cp0*v + rand(0,cp1)*(p-x)
cp3=cp3+alea(0,cp2)*(P_m[g].x.b[d]-P[s].x.b[d]); // + rand(0,cp2)*(g-x)
//Note that g is either the local best or the global best
// depending on the neighbourhood size
if(cp3>cp0) cp3=cp0; else if(cp3<-cp0) cp3=-cp0; // Take vmax into account
V[s].b[d]=(int)cp3 ; // New velocity
// Probabilistic binary choice for position
cp3=1/(1+exp(-cp3)); // Keep in ]0,1[
pr=alea(0,1);
if (pr<cp3) P[s].x.b[d]=1; else P[s].x.b[d]=0;
}
break;
case 111: // Pivot method -----------------------------------------
// Works pretty well on some problems .. and pretty bad on some others
Gp=P_m[g]; // Best previous position g of the best informer
dist=(int)log(D); // We suppose here D>=2
r=alea(1,dist);
//r=1+alea_normal(0,dist);// Test
for (k=0;k<r;k++)
{
d=alea_integer(0,D-1);
Gp.x.b[d]=1- Gp.x.b[d]; // Around g
}
P[s]=Gp;
break;
default:
printf("\n SORRY, no method %i",method);
break;
} // End of switch (method)
// ... evaluate the new position
for (i=0;i<n_f;i++) P[s].f.fi[i]=fabs(perf(s,function[i])-fmin.fi[i]);
switch (tribe)
{
default: // ... update the best position
if (best_(P[s].f,P_m[s].f)==1) P_m[s]=P[s];
break;
case 1: // ... or update the best position of the whole tribe
for(s1=0;s1<S;s1++)
{
if (best_(P[s].f,P_m[s1].f)==1) P_m[s1]=P[s];
}
break;
case 2: // update the worst
i=s; // rank of the worst
for(s1=0;s1<S;s1++)
{
if(s1==s) continue;
if(LINKS[s1][s]==0) continue;
if (best_(P_m[s1].f,P_m[i].f)==1) continue;
i=s1;
}
// if(i!=s) printf("\n %i modify %i",s,i);
if (best_(P[s].f,P_m[i].f)==1) P_m[i]=P[s];
break;
} // end switch (tribe)
// ... update the best of the bests
if (best_(P_m[s].f,best.f)) best=P_m[s];
// End of evaluation
if(print_level>0)
{
printf("\ncurrent ");
for (d=0;d<D;d++) printf("%1i",(int)P[s].x.b[d]);
for (i=0;i<n_f;i++) printf(" %f",P[s].f.fi[i]);
printf("\n best ");
for (d=0;d<D;d++) printf("%1i",(int)P_m[s].x.b[d]);
for (i=0;i<n_f;i++) printf(" %f",P_m[s].f.fi[i]);
printf("\n best best ");
for (d=0;d<D;d++) printf("%1i",(int)best.x.b[d]);
for (i=0;i<n_f;i++) printf(" %f",best.f.fi[i]);
}
}
// Check if finished
error_prev2=error_prev;
error_prev=error;
error=total_error(best.f);
//printf("\n error %f",error);
/*
// If no improvement, information links will be reinitialized
if(error<error_prev)
{
init_links=0;
}
else
{
init_links=1;
}
*/
init_links=1;
//fprintf(f_trace,"\n%i %f",nb_eval,save_time);
if (error>eps && nb_eval<eval_max) goto loop;
if (error>eps) n_failure=n_failure+1;
// Result display-save
printf("\n\nExec %i Eval= %i ",n_exec,nb_eval);
printf("\n");
for (d=0;d<D;d++) printf("%1i",(int)best.x.b[d]);
for (i=0;i<n_f;i++) printf(" %.4f",best.f.fi[i]);
if (n_exec==1) min=error; else if(error<min) min=error;
if (function[0]==11) // Special save for this case
{
DD=D/Sp;
for (i=0;i<Sp;i++)
{
printf("\n s%i ",i+1);
for (d=0;d<DD;d++) printf("%i",best.x.b[d+DD*i]);
for (d=0;d<DD;d++) fprintf(f_trace,"%i ",best.x.b[d+DD*i]);
fprintf(f_trace,"\n");
}
for (i=0;i<Sp-1;i++) // "particle" i
{
printf("\n");
for (j=i+1;j<Sp;j++) // "particle j
{
// Compute the distance
dist=0;
for (d=0;d<DD;d++)
{
if (best.x.b[d+DD*i]!=best.x.b[d+DD*j]) dist=dist+1;
}
printf(" %i",dist);
}
}
fprintf(f_trace,"\n");
}
// Save error for stats
error_[n_exec]=error;
printf("\n %i %f",n_exec,error_[n_exec]);
// Multiobjective result save
if (n_f>1)
{
for (i=0;i<n_f;i++)
fprintf(f_multi,"%f ",best.f.fi[i]);
for (d=0;d<D;d++) fprintf(f_multi," %.4i",(int)best.x.b[d]);
fprintf(f_multi,"\n");
}
if (adapt_K==1) printf("\n %f %f",K_choice[0][0],K_choice[2][0]);
// Compute and display some statistical information
eval_mean=eval_mean+nb_eval;
eps_mean=eps_mean+error;
if (n_exec<n_exec_max) goto init;
t2=clock();
printf("\n\n Time (clocks) %.0f",t2-t1);
save_time= save_time/eval_mean;
printf("\n Saved time rate %f",save_time);
eval_mean=eval_mean/(double)n_exec;
eps_mean=eps_mean/(double)n_exec;
printf("\n Eval. (mean)= %f",eval_mean);
printf("\n Error (mean) = %f",eps_mean);
success_rate=1- n_failure/(double)n_exec;
printf("\n Success rate = %.6f%%",100*success_rate);
if (n_exec>1) printf("\n Best min value = %f",min);
r=0;
for(k=1;k<=n_exec;k++)
{r=r+error_[k];
//printf("\n k %i r %f",k,r);
}
r=r/n_exec; // Mean
pr=0 ;
for(k=1;k<=n_exec;k++) pr=pr+(r-error_[k])*(r-error_[k]);
pr=sqrt(pr/n_exec);
printf("\n Mean %f StDev %f",r,pr);
fprintf(f_synth,"%f ",eval_mean);
fprintf(f_synth,"%f ",eps_mean);
fprintf(f_synth,"%.6f%% ",100*success_rate);
fprintf(f_synth,"%f ",min);
fprintf(f_synth,"%f ",save_time);
fprintf(f_synth,"\n");
if (eval_max==eval_max1) quality=success_rate*(alpha-pow(alpha,eval_max1-1));
else
if (eval_max<eval_max2)
quality=quality+success_rate*(pow(alpha,eval_max-delta_eval_max-1)-pow(alpha,eval_max-1));
else
quality=quality+success_rate*(pow(alpha,eval_max2-delta_eval_max-1));
//fprintf(f_synth,"\n %i %f ",eval_max,quality);
} // End loop on eval_max
z=quality/alpha;
printf("\n Quality: %f",z);
fprintf(f_synth,"\nQuality %f ",z);
printf("\n k-distribution stat."); for (k=1;k<=dist+1;k++) printf(" %.0f",k_tot[k]);
end:;
}
//===========================================================
double alea(double a,double b)
{ // rand number (uniform distribution) in [a b]
double r;
r=(double)rand_kiss()/RAND_MAX_KISS;
// r=(double)rand(); r=r/RAND_MAX;
return a+r*(b-a);
}
//===========================================================
int alea_integer(int a,int b)
{// Integer rand number in [a b]
int ir;
double r;
r=alea(0,1); ir=(int)(a+r*(b+1-a)); if (ir>b) ir=b;
return ir;
}
//===========================================================
double alea_normal(double mean, double std_dev)
{
/*
Use the polar form of the Box-Muller transformation to obtain
a pseudo random number from a Gaussian distribution
Note that it is always >= mean
*/
double x1, x2, w, y1;
//double y2;
do {
x1 = 2.0 * alea(0,1) - 1.0;
x2 = 2.0 * alea(0,1) - 1.0;
w = x1 * x1 + x2 * x2;
} while ( w >= 1.0);
w = sqrt( -2.0 * log( w ) / w );
y1 = x1 * w;
// y2 = x2 * w;
y1=y1*std_dev+mean;
return y1;
}
//===========================================================
int best_(struct f fx,struct f fy)
{
/*
Compare two position evaluations. In case of multiobjective
it is the Pareto dominance
Return 1 is fx better than fy, 0 else
*/
int i;
for(i=0;i<fx.size;i++)
{
if (fx.fi[i]>fy.fi[i]) return 0;
}
// All fx values are <=
for(i=0;i<fx.size;i++)
{
if (fx.fi[i]<fy.fi[i]) return 1; // At least one is <
}
return 0; //You may also return 2, if this particular case
// is interesting
}
//===========================================================
int best_info(int s,int option)
{
/*
Return the rank of the best informer
option = 0 => classical (really the best)
= 1 => using pseudo-gradient
*/
double delta1,delta2;
double dist;
double error_1,error_2;
int g;
double grad;
int m;
struct f min_f;
int s2;
switch (option)
{
case 1: // Pseudo gradient
g=s; // If nothing better, the best informer is
// the particle itself
grad=0;
error_1=total_error(P_m[s].f);
for (s2=0;s2<S;s2++)
{
if (LINKS[s2][s]==0) continue; // No link
error_2=total_error(P_m[s2].f);
delta1=error_1-error_2;
dist=distance(s,s2,0);
if (dist<=0) continue;
delta2=delta1/dist;
if (delta2>grad)
{
g=s2; grad=delta2;
}
}
// if (g==s) goto case0; // Possible variant
break;
case 0: // Classical option: the "real" best
for (g=0;g<S;g++) {if (LINKS[g][s]==0) continue; goto next;}
next:
min_f=P_m[g].f;
for (m=g+1;m<S;m++)
{
if (LINKS[m][s]==0) continue;
if (best_(P_m[m].f,min_f)==1) {g=m;min_f=P_m[m].f;}
}
break;
}
return g;
}
//===========================================================
struct bits binary(unsigned int a, int size)
{
/* Compute the bit string (truncated at "size")
corresponding to the integer a
*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -