color_kit.c

C语言开发的微粒群优化算法源程序

		tot_min=total;
	}
}

// Find the nearest different color already used

new_color=(float)-NA; // Arbitrary too big value

for (i=0;i<color_list.size;i++)
{
	c=color_list.x[i];
	if (c==worst_color) continue;
		if (fabs(c-worst_color)<fabs(new_color-worst_color))
		{
			new_color=c;
		}
}

//printf("\n minimize_color, worst %f , new %f",worst_color,new_color);
// Replace the worst color

for (node=0;node<M.size;node++)
{

	c=partt.pos.x[node];

	if (c==worst_color) 
	{
		partt.pos.x[node]=new_color;
	}
}
goto check;

}


// ----------------------------------------------------------------------------- MINIMIZE_CONSTRAINT
struct particle	minimize_constraint(struct particle part)
{
/* 
For Coloring problem/Frequency assignment, try to adapt the position so that all constraints are respected, 
eventually by increasing the number of colors used

Note: matrix M and table dplus are global variables

while there is a non respected constraint
	choose the node which has the greatest number of  non respected constraints
	modify its color as few as possible so that all its constraints are respected


*/
float			c0,c,c1;
float			d;
int				i;
int				j;
int				k;

float			m;
float			total;
float			tot_max;

struct particle partt;
int				worst_node;



partt=part;

check:

// Find the worst node 

tot_max=0;
 
for (i=0;i<M.size;i++)
	{
		total=0;
		c=partt.pos.x[i]; // color of node i
		for (j=0;j<dplus.dplus[i].size;j++)
		{
			k=dplus.dplus[i].val[j]; // k = neighbour of node i
			c1=partt.pos.x[k]; // color of node k
			d=M.val[i][k]-(float)fabs(c-c1);
			if (d>0) 
			{
				total=total+d; // How much this constraint is not satisfied
			}
		}
		if (total>tot_max)
			{
			tot_max=total;
			worst_node=i;
			}
	}
	
if (tot_max==0) goto end;	
	
// If found, change its color as few as possible so that all
// its constraints are respected, if possible	
	
m=0;
c0=partt.pos.x[worst_node];

next_color:
m=m+1;	
c=c0+m; // color of worst_node

for (j=0;j<dplus.dplus[worst_node].size;j++)
{
	k=dplus.dplus[worst_node].val[j]; // k = neighbour of worst_node
	c1=partt.pos.x[k]; // color of node k
	d=M.val[worst_node][k]-(float)fabs(c-c1);
	if (d>0) 
	{
		goto minus_m;
	}
}
 goto assign_color;

minus_m:
c=c0-m;
for (j=0;j<dplus.dplus[worst_node].size;j++)
{
	k=dplus.dplus[worst_node].val[j]; // k = neighbour of worst_node
	c1=partt.pos.x[k]; // color of node k
	d=M.val[worst_node][k]-(float)fabs(c-c1);
	if (d>0) 
	{
		goto next_color;
	}
}

assign_color:
partt.pos.x[worst_node]=c;
goto check;



end:


return partt;

}

// ----------------------------------------------------------------------------- READ_MATRIX_COLOR
struct matrix read_matrix_color(FILE *f_data, float used_constr)
{
int				i;
int				j;
struct matrix	M;
int				nothing; // Just to skip column and row number by reading the file
float			z;

fscanf (f_data,"%i",&M.size); // WARNING ***: Must be equal to param.DD

for (i=0;i<M.size;i++) fscanf(f_data,"%i",&nothing); // Column numbers

for (i=0;i<M.size;i++)
	{
	fscanf(f_data,"%i",&nothing); // row number
	for (j=0;j<M.size;j++)
		fscanf(f_data,"%i",&M.val[i][j]);
	}

// Check matrix

M.nb_cont=0; // Number of constraints
M.max_cont=0; // Maximum constraint
for (i=0;i<M.size-1;i++)
{
	for (j=i+1;j<M.size;j++)
	{
		if (M.val[i][j]>0)
		{ 
			M.nb_cont=M.nb_cont+1; // Number of constraints
			if (M.max_cont<M.val[i][j]) M.max_cont=(float)M.val[i][j];
		}
		if (M.val[i][j]!=M.val[j][i])
		{
			printf("\nERROR, read_matrix. matrix not symmetrical (%i, %i)",i,j);

infinite:;goto infinite;
		}
	}
}

// Eventually, simplify it, relaxing some constraints

if (used_constr<1)
	{
	printf("\n WARNING! As required, I use only %f percent of the constraints",100*used_constr);
	for (i=0;i<M.size-1;i++)
		{
		for (j=i+1;j<M.size;j++)
			{
			if (M.val[j][i]>0)
				{
				z=alea_float(0,1);
				if (z>used_constr)
					{
					M.val[j][i]=0;
					M.val[i][j]=0;
					}	
				}
			}
		}
	}

// Analyze the matrix (list of neighbours for each node)
// Note: dplus is a global variable

for (i=0;i<M.size;i++)
{
	dplus.dplus[i].size=0;
	for (j=0;j<M.size;j++)
	{
		if (j==i) continue; // No loop
		if (M.val[i][j]>0)
		{
			dplus.dplus[i].val[dplus.dplus[i].size]=j; // Build the list of neighbours
			dplus.dplus[i].size=dplus.dplus[i].size+1;
		}
	}
}

dplus.dmax=dplus.dplus[0].size;
dplus.dtot=dplus.dmax;

for (i=1;i<M.size;i++) 
{
	if (dplus.dplus[i].size>dplus.dmax) dplus.dmax=dplus.dplus[i].size; // Max of neighbours
	dplus.dtot=dplus.dtot+dplus.dplus[i].size; // 2*(Total number of constraints)
}



//display_matrix(M);

return M;
}


// ----------------------------------------------------------------------------- TOT_COLOR
float tot_color(struct position pos,struct matrix M,int cmax,int save)
{
/* Objective function for Cloring/Frequency Assignment problem
M is the symmetrical constraint matrix.
cmax is the maximum number of colors wanted

We must have abs(pos(i)-pos(j))>=M(i,j)
So we define d(i,j)=0 if it's OK
                   =M(i,j)-abs(pos(i)-pos(j)) else
The first part of objective function value is the sum of all d(i,j), + 1

The second part is depending on the number of different colors/frequencies used, compared to the wanted value.
We define the function so that its maximum value is <1.

We assume the position is on a projective hyperplan, so that all coordinates have to be
divided by the dimension, in order to find the right pos(i)-pos(j)

*/

float	coeff1,coeff2;
float	DD;
float	d;
//float	fr[100];
int		i;
int		j;
int		k;
int		n_col;
int		option;
float 	total1,total2,total;

//printf("\n tot_color");
//display_position(pos);


option=0; // it is more important to have the minimum number of colors
option=1; // It is more important to have no interference


coeff1=1; // For number of unsatisfied constraints
coeff2=1;// For number of different colors


DD=(float)pos.size; // (Must be equal to M.size)
total1=0;
for (i=0;i<DD;i++) 
{
	for (k=0;k<dplus.dplus[i].size;k++)
	{
		j=dplus.dplus[i].val[k];
		
		//d=M.val[i][j]-fabs(pos.x[i]-pos.x[j])/DD;
		d=M.val[i][j]-(float)fabs(pos.x[i]-pos.x[j]);
		if (d>0) 
		{
			total1=total1+1;// Add to the number of unsatisfied constraints (* 2)
			//total1=total1+d; // How much this constraint is not satisfied (method 1)
			//total1=total1+d*d; // How much this constraint is not satisfied (method 2)
		}
	}
}

/*
if (total1>dplus.dtot)
{
	printf("\n ERROR. tot_color. total1 %f > dplus.dtot %i",total1,dplus.dtot);
	infinite:;goto infinite;
}
*/

total1=total1/dplus.dtot; // On [0,1]

//Take into account the number of colors used

n_col=diff_color(pos);

if (n_col>cmax)
	{
		//total=total+(n_col-cmax)/(M.size+1);

		total2=(float)(n_col-cmax);
	}
else
	total2=0;
	
if (n_col>=M.size)// On [0,1]
{
	total2=1;
}
else
{
	total2=total2/(M.size-cmax); 
}



if (save==12)
{
	fprintf (f_move_color,"\n %.0f %f %f",pos.label,total1,total2);
}

total=(coeff1*total1+coeff2*total2)/(coeff1+coeff2); // On [0,1]

//printf("\nlabel %.0f.  total1, total2, total %f %f %f",pos.label,total1, total2,total);

X1=total1; // For Leonardo visualization
X2=total2; 

return total;


}

// ----------------------------------------------------------------------------- UPDATE_VEL_COLOR
struct velocity update_vel_color(struct position pos,struct velocity v0,struct velocity v2,float alpha,float b)
{
/*
Note, dplus is a global variable for a edge list representation of M

Option 0: classical weighting
Option 1: "combinatorial" one
*/
float			ca,cab,cb;
struct vector	g;
int				i,j,k;
int				n;
int				n1,n2;
int				option;
struct vector	used={0};
struct velocity v;//{int size;float v[Max_vel];float ek;float v2[Max_vel];int combin;}

option=0;

if (option==0)
{
	v.size=v0.size;
	v.combin=v0.combin;

	for (i=0;i<v.size;i++)
	{
		v.v[i]=alpha*v0.v[i]+b*v2.v[i];
	}


	return v;
}

//  OLD 


v=v0; // First, just take the initial velocity
cab=alpha+b;
ca=alpha/cab;
cb=b/cab;
n=(int)(0.5+cb*v2.size);// Evaluate how many components of v2 has to be taken into account
if (n==0) return v;


//printf("\n update_vel_color, coeffs %f %f, n=%i",ca,cb,n);
/*
display_velocity(v0);
display_velocity(v2);
*/


/* --------
The graph has N nodes.
- choose a node at random
- from this node, build a connexe graph (as dense as possible) with n nodes
- take just the velocity components corresponding to the nodes of this subgraph

*/

//i=alea(0,v2.size-1); // Choose the first node
i=max_color_rank(pos);
g.x[0]=(float)i;
g.size=1;
used.x[i]=1;
//printf("\n   nodes %i",i);

if (n>1)
{
	for (j=0;j<g.size;j++) // For each node already used
	{
		n1=(int)g.x[j]; // n1 = node already used
		
		for (k=0;k<dplus.dplus[n1].size;k++)
		{
			n2=dplus.dplus[n1].val[k];
			if (used.x[n2]==1) continue; // Already used
			g.x[g.size]=(float)n2;		
			used.x[n2]=1;
			g.size=g.size+1;  // Note the trick: the upper bound of the loop is modified inside the loop
//printf(" %i",n2);
			if (g.size>=n) goto end;
		}
	}
}

end:

for (i=0;i<g.size;i++) 
{
	j=(int)g.x[i];
	v.v[j]=v2.v[j];
}

//display_velocity(v);
printf("\n %f",v.v[0]);

return v;
}

// ----------------------------------------------------------------------------- UPDATE_VEL_COLOR2
struct velocity update_vel_color2(struct position pos,struct velocity v0,struct velocity v2,float alpha,float b,int cmax)
{
/*
Note: the matrix M and the edge list dplus are global variables
		H is also a global variable (for projections).

x is a position/coloring C0
x+v0 is the next position C1 if the particle follows its own way
x+v2 is another position/coloring C2 (usually the best ever found the neighbourhood

The particle moves (virtually) to C1, and update this position according to C2 and the coeffs, 
but only if it improves it.
*/

int					i;
int					j;
int					k,k_tot;
float				lambda;
int					loop;
int					loop_max;
struct vector_int	nodes;
struct vector_int	not_used;
struct velocity		v;
struct position		pos1,pos2,pos3;


loop_max=dplus.dmax;

//printf("\n update_vel_color2. loop_max=%i. begin",loop_max);
//display_position(pos);
//display_velocity(v2);

pos1=pos;
pos2=pos;

for (i=0;i<pos.size;i++) 
{
	pos1.x[i]=pos.x[i]+v0.v[i];// Define position 1, if the particle follows its own way
	pos2.x[i]=pos.x[i]+v2.v[i]; // Position2, Usually the best known position in the neighbourhood
}

lambda=alpha/(alpha+b);
k_tot=(int)(1-lambda)*pos.size;

if (k_tot<1) k_tot=1;

//printf(" (k_tot=%i)",k_tot);

loop=0;
loop:

//printf(" %i",loop);

// Define a list of k nodes to modify

not_used.size=pos.size;

for (i=0;i<not_used.size;i++)
{
	not_used.val[i]=i;
}

nodes.size=0;

for (k=0;k<k_tot;k++)
{
	i=alea(0,not_used.size-1);
	nodes.val[nodes.size]=not_used.val[i]; // Node to use
	nodes.size=nodes.size+1;

	// Compact not_used

	if (i==not_used.size-1)
	{
		not_used.size=not_used.size-1;
		continue;
	}

	for (j=i;j<not_used.size-1;j++)
	{
		not_used.val[j]=not_used.val[j+1];
	}
	not_used.size=not_used.size-1;
}

// Modify the nodes


pos1.f=tot_color(pos1,M,cmax,0);
pos3=pos1;

for (k=0;k<nodes.size;k++) 
{
	i=nodes.val[k];
	pos3.x[i]=pos2.x[i];	
}

// Check if better

pos3.f=tot_color(pos3,M,cmax,0);

if (pos3.f<pos1.f) goto end;


	// Try again
	loop=loop+1;
	if (loop<loop_max) goto loop;

	
end:

//printf("\n update_vel_color2. %f %f",pos1.f,pos3.f);
//display_position(pos3);

v=v0;
for (i=0;i<v.size;i++)
{
	v.v[i]=pos3.x[i]-pos.x[i];
}


//printf("\n update_vel_color2 pos3(10), pos(10), v(10) %f %f %f",pos3.x[10],pos.x[10],v.v[10])
//printf("  end");

return v;
}