myfunction.c

用C编写的部落寻优的优化算法

// ----------------------------------------------------------------------------- MYFUNCTION
/*
   You can add any function you want.
   Don't forget to also modify functions.txt
*/
struct f MyFunction(struct position pos,int funct,double target,int level)
{

// Position evaluation

double	a_1,a_2,a_3;
double	alpha,beta,gamma,delta,p,q;
double	better;
double	b;
struct vector_i bit,bit_x;
int		bit_sum;
double   c,c1,c2;
int			choice;
int		d,dmax;
double	discr;
int		D;
int 			DD;
double	f_model;
double   f;
double      f1,f2;
double	gen[Max_DD]; // For Moving Peaks
int		i;
int		ix,ix2;
int     j;
int     k;
int		m;
int		n1,n2;
int		num;
//struct position post={0};
struct position post;
double      pr;
double			product;
double      r;
double			rho;
double			s;
double  sigma;
double			sum1,sum2;
double			theta;
struct f 			total;
int		used[Max_DD]={0};
double	x1,x2,x3,x4,x5,x6,x7,x8,x9;
double			x;
double      xd;
double			xid;
double	y,y1,y2,y3,y4;
double	z1,z2,z3;
struct roots z;

// Data for model tuning example (see case 22)

static float 	data[5][2] =
{
{	10	,	1480	},
{	20	,	2200	},
{	30	,	5000	},
{	40	,	5880	},
{	50	,	11590	}
};


static double	Fmax=1000.0;
static double	S=189000.0;
static double	lmax=14.0;
static double	dmin=0.2;
static double	Dmax=3.0;
static double	Fp=300;
static double	spm=6.0;
static double	sw=1.25;
static double	G=11500000;

// Variables specific to Coil compressing spring
// See also constrain.c
double Cf,K,sp,lf;

//----------------
int d2=5;

// Data for Master Mind  (hard coded solution)
static int data_M[4] =
{ 1,2,3,4};

/*
// Specific Informational table for 4 different colors
static int r_w[5][5]=
{

{16	,192,	120,	20,	1},
{152,	192,	24,	0,	0},
{312,	108,	6,	0,	0},
{136,	8,	0,	0	,0},
{9	,0,	0	,0,	0}
};
 */
// For Foxholes problem
static int a[2][25] =
{
   {-32, -16, 0, 16, 32, -32, -16, 0, 16, 32, -32, -16, 0, 16, 32,
      -32, -16, 0, 16, 32, -32, -16, 0, 16, 32  },
   {-32, -32, -32, -32, -32, -16, -16, -16, -16, -16,
      16, 16, 16, 16, 16, 32, 32, 32, 32, 32 }
};

// For Jeannet-Messine problem (cf. ROADEF 2003, p 273)
static double a1[6]=
{ 0.5, 0.3, 0.8, 0.1, 0.9, 0.12};
static double a2[6]=
{-0.5, 0.6, 0.1, 1.5, -1, 0.8};

// For polynom fittint
   int const M=60;
   double py, dx=(double)M;
 //-----------------------------------------
 
eval_f[level]=eval_f[level]+1;
eval_f_tot[level]=eval_f_tot[level]+1;

if (problem[level].printlevel>2)
	printf("\n my_function. tot. eval_f_tot %6.0f",eval_f_tot[level]);

DD=pos.p.size;
total.size=pos.f.size;

switch (funct)
{

case 1:  //DD-sum
	total.f[0] = 0;
	for( d=0;d<DD;d++)                       //for each dimension
	{
		total.f[0] = total.f[0] + pos.p.x[ d];               //add DD if i's coordinates
	}

	break;

case 2:  // DD-product.
	total.f[0] = 1;
	for( d=0;d<DD;d++)
		total.f[0] = total.f[0]*pos.p.x[ d];

	break;

case 3: //Parabole () Sphere
	//      Global min f(x)=0, x(i)=0
	total.f[0] = 0;
	for( d=0;d<DD;d++)
	{
		total.f[0] = total.f[0] + pos.p.x[ d]*pos.p.x[ d];
	}
	break;

case 4: // Rosenbrock's valley, Banana function
//  Global min f(x)=0, x(i)=1
	total.f[0]=0;

	for (d=0;d<DD-1;d++)
		{
		xid=1-pos.p.x[d];
		total.f[0]=total.f[0]+xid*xid;
		xid= pos.p.x[d]*pos.p.x[d]-pos.p.x[d+1];
		total.f[0]=total.f[0]+100*xid*xid;
		}
	break;


case 5: // Clerc's f1, Alpine function, min 0 in (0,...,0)
	total.f[0]=0;
	/*
	for( d=0;d<DD;d++)
		{
		xid=sin(pos.p.x[d]);
		total.f[0]=total.f[0]+sqrt(fabs(pos.p.x[d]*xid));
		}
		*/
	for( d=0;d<DD;d++)
	{
		xid=pos.p.x[d];
		total.f[0]=total.f[0]+fabs(xid*sin(xid)+0.1*xid);
	}
	break;

case 6: // Griewank's function. Min =0 at (100,100 ... 100)
	if (DD==2) goto Gr2;
	total.f[0]=0;
	product=1;
	for (d=0;d<DD;d++)
		{
		xid=pos.p.x[ d]-100;
		total.f[0]=total.f[0]+xid*xid;
		product=product*cos (xid/sqrt(d+1));
		}
	total.f[0]=total.f[0]/4000 -product +1;
	break;

Gr2:	// -------- 2D only 
	x1=pos.p.x[0];
	x2=pos.p.x[1];
	total.f[0]=(x1*x1+x2*x2)/2 - cos(two_pi*x1)*cos(two_pi*x2)+1;
	break;


case 7:  // Rastrigin. Minimum value 0. Solution (0,0 ...0)
	k=10;

	total.f[0]=0;
	for (d=0;d<DD;d++)
		{
		xid=pos.p.x[ d];
		total.f[0]=total.f[0]+xid*xid - k*cos(two_pi*xid);
		}
	total.f[0]=total.f[0]+DD*k;

	break;

case 8:  // Fifty-fifty problem
	total.f[0]=tot_fifty_fifty(pos);
	break;


case 9:  // Ackley
	sum1=0;
	sum2=0;
	x=DD;
	for (d=0;d<DD;d++)
		{
		xid=pos.p.x[ d];
		sum1=sum1+xid*xid;
		sum2=sum2+cos(two_pi*xid);
		}
	total.f[0]=(-20*exp(-0.2*sqrt(sum1/x))-exp(sum2/x)+20+E);
	break;


case 10: //Foxholes 2D
	total.f[0]=0;
	for (j=0; j<25; j++)
   		{
		sum1=0;
      	for (d=0; d<2; d++)
      		{
         	sum1 =sum1+ pow (pos.p.x[d] - a[d][j],6);
      		}
      	total.f[0]=total.f[0]+1/(j+1+sum1);
   		}
	total.f[0]=1/(0.002+total.f[0]);
	break;

case 11: // ======================= Apple trees
	total.f[0]=apple_trees(pos);
	post= homogen_to_carte(pos); // For Leonardo visualization (for example)
	break;


case 12: // ======================== El Farol
/*
1 dimension = 1 Irishman
Just two values/dim.: 0 <=> stay at home
                      1 <=> go to the pub

Note: the current swarm sw is a global variable
*/
	sum1=60; // Maximum "tolerance": maximum number of people each Irishman does accept in the pub
	total.f[0]=0;
	for (d=0;d<DD;d++)
	{
		total.f[0]=total.f[0]+pos.p.x[d];
	}

	total.f[0]=total.f[0]*(1-total.f[0]/sum1);
	break;

 case 13: //  ======================Fermat

/* Suggested parameters:
- dimension 3
- xmin 2
- xmax 100
- target 0
- granul 1
- eps 0





You will then find three integer numbers x,y,z
so that x^2 + y^2 = z^2
*/

total.f[0] = 0;
dmax=pos.p.size-1;
for( d=0;d<dmax;d++)  //for each dimension

	{
	total.f[0]=total.f[0]+pos.p.x[d]*pos.p.x[d];

	}
total.f[0]=fabs(total.f[0]-pos.p.x[dmax]*pos.p.x[dmax]);
break;

case 14: //  ===========================================Knapsack
/* Suggested parameters:
- dimension 10
- xmin 1

- xmax 100
- target 100
- granul 1
- eps 0
*/
total.f[0] = 0;
	for( d=0;d<pos.p.size;d++)
		total.f[0] = total.f[0] + pos.p.x[ d];
break;


case 15: //  ===================================== (x*x+y*y-1)*(x*x+y*y-1)+(sin(10*x)-y)*(sin(10*x)-y)
// dimension =2
x=pos.p.x[0];

y=pos.p.x[1];
total.f[0]=(x*x+y*y-1)*(x*x+y*y-1)+(sin(10*x)-y)*(sin(10*x)-y);
break;

case 16: //  =====================================  Intersection of two circles
//dimension =2
x=pos.p.x[0];
y=pos.p.x[1];
total.f[0]=fabs(x*x+y*y-1)+fabs((x-2)*(x-2)+y*y-4);
break;

case 17: //  =====================================  fabs(sin(x)-y-sqrt(3)/2+1)+fabs(x*x-log(y)-(pi/3)*(pi/3))
// dimension =2 , in [0,10]
x=pos.p.x[0];
y=pos.p.x[1];
total.f[0]=fabs(sin(x)-y-sqrt(3)/2+1)+fabs(x*x-log(y)-(pi/3)*(pi/3));
break;

case 18: //  ===================================== 2D Linear system
// dimension =2, solution (5,2)
x=pos.p.x[0];
y=pos.p.x[1];
total.f[0]=fabs(3*x+2*y-19)+fabs(2*x-y-8);
break;

case 19: //  ======================================== sin wave
// dimension 1
x=pos.p.x[0];
total.f[0]=sin(x);
break;


case 20: //  =================== 3D linear system. Chinese problem of the 100 fowls
// Several integer solutions
//printf("\n\n");
	total.f[0]=fabs (pos.p.x[0]+pos.p.x[1]+pos.p.x[2]-100);
	total.f[0]=total.f[0]+fabs (5*pos.p.x[0]+3*pos.p.x[1]+(1/3)*pos.p.x[2]-100);
break;


case 21: // ===================================== Magic square
d2=(int)(sqrt(pos.p.size)); // Nb of lines = Nb of columns=sqrt(dimension)
total.f[0]=0;

for (i=0;i<d2-1;i++)  // Rows
					// For each pair of rows, compute the difference x
					// of the sums. "distance" to minimize =x*x
{

	for (j=i+1;j<d2;j++) 
	{
		x=0;	
		for (d=0;d<d2;d++)
			x=x+pos.p.x[i*d2+d]-pos.p.x[j*d2+d];
	}
	total.f[0]=total.f[0]+x*x;
}

for (i=0;i<d2-1;i++)  // Columns
{
	for (j=i+1;j<d2;j++)
	{
		x=0;
		for (d=0;d<d2;d++)
		x=x+pos.p.x[i+d*d2]-pos.p.x[j+d*d2];
		total.f[0]=total.f[0]+x*x;
	}
}

break;

case 22: // ======================== Model tuning ================
/* We have some data and a 2 parameters model
Find the "best" set of parameters
*/

/*
Model lambda*(D^mu)
column 1: D
column 2: function value
*/

total.f[0]=0;
for (d=0;d<d2;d++)
	{
	f_model=pos.p.x[0]*pow(data[d][0],pos.p.x[1]);
	total.f[0]=total.f[0]+(f_model-data[d][1])*(f_model-data[d][1]);
	}
 total.f[0]=sqrt(total.f[0]);  
break;

case 23: // =================== 4 positions 6 colors Master Mind

/*
Just a test. PSO alone is quite bad for this problem.
If you are really interested in Master Mind, you should have a look at my "ultimate" program,
on my math web site: you can't beat it.
*/

// "Normal" estimation
n1=0; // right color, right position
n2=0; // right color, wrong position


for (d=0;d<4;d++) //Right position
{
	ix=(int)pos.p.x[d];
	if (data_M[d]!=ix) continue;
		 n1=n1+1;// Right position
		used[ix]=1;
}


for (d=0;d<4;d++) // Wrong position
{
	ix=(int)pos.p.x[d];
	if (used[ix]>0) continue;

	for (d2=0;d2<4;d2++)
	{
		if (d2==d) continue;
			ix2=data_M[d2];
			if (ix2!=ix) continue;
				if (used[ix2]>0) continue;
					n2=n2+1;
					used[ix]=1;
	}
}

total.f[0]=400-100*n1 - n2;
break;


case 24: //======================= Catalan's conjecture  x^m - y^n = 1 has just one integer solution (3^2 - 2^3 = 1)
/* dimension = 2x2 =4
  a) the solution (3,2,2,3) is easily found
  b) you may try as many times as you want, with any initial position, you'll never find another one

    => the conjecture is _probably_ true
*/
total.f[0]=pow(pos.p.x[0],pos.p.x[1])-pow(pos.p.x[2],pos.p.x[3])-1;
if (pos.p.x[0]==pos.p.x[2]) total.f[0]=total.f[0]+pow(pos.p.x[0],pos.p.x[1]); // Just to avoid the local minimum x0=x2 and x1=x3
break;

case 25:
break;

case 26: //================================ Chinese-MIT problem

a_1=100;

a_2=200;
a_3=300;
total.f[0] = 0;
	for( d=0;d<pos.p.size;d++)
		total.f[0] = total.f[0]+ pos.p.x[d]; // Any positive function


n1=(int)(float)pos.p.size/3;
n2=2*n1;

x1=0;
for (d=0;d<n1;d++) x1=x1+pos.p.x[d];
x2=0;
for (d=n1-1;d<n2;d++) x2=x2+pos.p.x[d];


x3=0;
for (d=n2-1;d<pos.p.size;d++) x3=x3+pos.p.x[d];

x1=(a_1-x1);
x2=(a_2-x2);
x3=(a_3-x3);
total.f[0]=total.f[0]+x1*x1+x2*x2+x3*x3;
//total.f[0]=x1*x1+x2*x2+x3*x3;

break;

case 28: // ================================== Test surface 2D
	x1=pos.p.x[0];
	x2=pos.p.x[1];
	total.f[0]=20*x1*sin(3*x1*x2)+5*x2*sin(2*x1);

	break;


case 29: // =================          Cognitive harmony
/*
P  = Weight (square) matrix. Weights in [-1,1]
x = Activation vector. Activations in [0,1]
=> harmony = transp(x)Px, to maximize

In fact, we compute here a dissonance to minimize
*/

// Evaluate the position
MM=pos.p.size;
total.f[0]=0;
for (d=0;d<MM;d++) // Compute harmony
{
	for(m=0;m<MM;m++) {total.f[0]=total.f[0]+pos.p.x[d]*problem[level].P.val[d][m]*pos.p.x[m];}
}
total.f[0]=MM*(MM-1)-total.f[0]; // Dissonance
 //total.f[0]=MM-total.f[0];    //**** WARNING, THIS MAY GIVE NEGATIVE RESULT
break;

 case 30: // Non specific TSP
 /*
   Just a test, can't be very good. See specific PSO_for_TSP for better result.
   For this problem, you must use
   granularity=1 (normal) or 0 (continuous)
   all_different=1
 */
// Evaluate the position
total.f[0]=f_tsp(pos,1,1,-1,-1,level);
 break;

 case 31: // Sum of absolute values
   // Note: you have a classical integer problem by setting granularity to 1, in the problem description
   	total.f[0] = 0;
	for( d=0;d<DD;d++)                       //for each dimension
	{
		total.f[0] = total.f[0] + fabs(pos.p.x[ d]);               //add DD if i's coordinates
	}
 break;

  case 32: // Non specific QAP
 /*
   Just a test, can't be very good.
   For this problem, you must use
   dimension = MM
   granularity=1 (normal)
   all_different=1
 */

// Evaluate the position

total.f[0]=f_qap(pos,1,1,-1,-1,level);
//printf("\nmyfunction. total %f",total.f[0]);

 break;



 case 33: //  Multiobjective Lis-Eiben 1
 /*
    Each run with a different random initialization
    gives a point whose (f1,f2) is near of the Pareto front
*/
x1=pos.p.x[ 0];
x2=pos.p.x[ 1];
x1=x1*x1+x2*x2;
f1=pow(x1,1.0/8);

x1= pos.p.x[ 0]-0.5;
x2= pos.p.x[ 1]-0.5;
x2=x1*x1+x2*x2;

f2=pow(x2,1.0/4);

total.f[0]=f1;
total.f[1]=f2;
break;

 case 34:   // Multiobjective Schaffer
 // x in (-5,10]
    x=pos.p.x[ 0];
    f2=(x-5)*(x-5);
    if (x<=1) {f1=-x;  goto end_34;}
    if (x<=3) {f1=-2+x;  goto end_34; }
    if (x<=4) {f1=4-x;  goto end_34;}
    f1=-4+x;
    end_34:
 total.f[0]=f1+1;   // Just to have a sure positive value
total.f[1]=f2;
   break;

      case 35:   // Multiobjective F1
 // x1 in [0,1]
 //   x2 in [0,1]

    x1=pos.p.x[ 0];
    x2=pos.p.x[ 1];
    f1=(x1*x1+x2*x2)/2;
    x1=x1-2;
    x2=x2-2;
    f2=(x1*x1+x2*x2)/2;

 total.f[0]=f1;
total.f[1]=f2;
   break;

         case 36:   // Multiobjective F2
 // x1 in [0,1]
 //   x2 in [0,1]


    x1=pos.p.x[ 0];
    x2=pos.p.x[ 1];
    f1=x1;
    x=1+9*x2;
    f2= x-sqrt(f1*x);

 total.f[0]=f1;
total.f[1]=f2;
   break;
   
    case 37:   // Multiobjective F3
 // x1 in (0,1]
 //   x2 in (0,1]

    f1=pos.p.x[ 0];
    x2=pos.p.x[ 1];
    x=1+9*x2;
    f2=x*(1-f1*f1/(x*x));