/*
Combinatorial kit for Particle Swarm Optimization
Last update: 2001-10-25 Maurice.Clerc@WriteMe.com

A "velocity" is essentially a list of pairs of values (v,v2), meaning
"assign value v according to value v2"

Special v2 values, when velocity applied to a position:
NA means "any value at random"
NO means "no modification"
AS means "assign v"

The "norm" of a velocity is the number of v2 not equal to NO

*/

/*------------------------------------------------------------------- ALEA_VELOCITY */
struct velocity	alea_velocity(int size,struct param param)
{ /* Just for combinatorial kit

NA value for about "size" components, NO either

*/
int				i;
int				num;
struct velocity vt;


vt.size=param.DD;
vt.combin=param.combin;

for (i=0;i<vt.size;i++)
	{
	num=alea(1,vt.size);
	if (num>size) vt.v2[i]=NO;
		else vt.v2[i]=NA;
	}	
return vt;
}

/*------------------------------------------------------------------- AUTO_MOVE */
struct particle auto_move(struct particle part,struct param param, struct model model)
/* The particle moves slowly, looking for a better position */
{
int 			i;
struct particle	partt,partt1;
int				type;
struct velocity v;
float			vmax;

type=model.local_search;
if (type==0) return part;
partt=part;


if (type==1) // Lazy Descent Method
	{
	for (i=0;i<part.pos.size;i++) // Move at random. 
		{
		vmax=(param.H.max[i]-param.H.min[i])/param.H.dx[i];
		v=alea_velocity(alea(1,(int)vmax),param);
		partt.pos=pos_plus_vel(partt.pos,v,param);

		partt.pos.f=tot(partt.pos,param,model);
		partt.pos.f=eval(partt.pos.f,param.target);

		if(partt.pos.f<part.pos.f) // Stop as soon as a better position is found
			{
			return partt;
			}
		}
	return part;	
	}

	
if (type==2) // Energetic
	{

	partt1=partt;
	energetic:
	for (i=0;i<part.pos.size;i++) // Move at random
		{
		vmax=(param.H.max[i]-param.H.min[i])/param.H.dx[i];
		v=alea_velocity(alea(1,(int)vmax),param);
		partt1.pos=pos_plus_vel(partt.pos,v,param);
		
		partt1.pos.f=tot(partt1.pos,param,model);
		partt1.pos.f=eval(partt1.pos.f,param.target);

		if(partt1.pos.f<partt.pos.f) // Move as long as a better position is found
			{
			partt=partt1;
			goto energetic;
			}
		}
	return partt;	
	}	
	
	
if (type==3) /* Move as long as a better position is found
				Check all immediate physical neighbours. If one is better, move, and check again
			*/

	{
	v.size=1;
printf("\n WARNING (auto_move). type=2 is not implemented !");
printf("\n    I use type=1 (energetic) instead"); goto energetic;

	return partt;
	}



printf("\n ERROR. auto_move: unknown type");

exit(1);		
}




/*------------------------------------------------------------------- COEFF_TIMES_VEL */	
struct velocity	coeff_times_vel(struct velocity v,float coeff, struct param param)
{
float			coefft;
int				i,j,k,l;
int				n;
int				rank[Max_vel];
int				rank_size;
struct velocity vt;

coefft=coeff;
if (param.printlevel>2)
{
	printf("\n coeff_times_vel %f",coeff);
	display_velocity(v);
}

if (coeff<0)
{
	//printf("\n *** WARNING. coeff_times_vel, coeff <0 for combinatorial process: %f",coeff);
	//return v;
	coefft=-coeff;
}


if (coeff>1)
{
	//printf("\n *** WARNING. coeff_times_vel, coeff >1 for combinatorial process: %f",coeff);
	//return v;
	//coefft=1;
	coefft=coeff-(int)coeff;

}
 if (coeff==1) 
{
	return v;;
}

vt=v;

if (coeff==0) 
{
	for (i=0;i<vt.size;i++)
	{
		vt.v2[i]=NO; // "null" velocity
	}
	return vt;
}


	n=(int)MAX(1,(1-coefft)*v.size);
//printf("\n v.size %i,  n  %i",v.size,n);
	
	rank_size=vt.size;
	for (i=0;i<rank_size;i++) rank[i]=i;
	
	// Choose at random n different components and set them to "zero"
	for (i=0;i<n;i++)
	{
		j=alea(0,rank_size-1);
		l=rank[j];
		//vt.v2[l]=NA;
		
		vt.v2[l]=NO; 
		
		// Compact rank[]
		if (j<rank_size-1)
		{
			for (k=j;k<rank_size-1;k++)
			{
				rank[k]=rank[k+1];
			}
			rank_size=rank_size-1;
		}
	
	}
if (param.printlevel>2)
{
	display_velocity(vt);
}
return vt;
}

//------------------------------------------------------------------- COEFF_TIMES_VEL2 	
struct velocity	coeff_times_vel2(struct velocity v1,struct velocity v2,float c1,float c2, int rank_d)
{
/*
c1*v1 + c2*v2

c1 >=0
c2 >=0

*/
int		i;
int		j;
int		k;
int		l;
int		n;
int		rank_size;
int		rank[Max_vel];
struct 	velocity vt;
float	c;

/*
printf("\n coeff_times_vel2: c1 %f, c2 %f => %f*v1 + %f*v2",c1,c2,c1/(c1+c2),c2/(c1+c2));
display_velocity(v1);
display_velocity(v2);
*/

if (rank_d>=0) // Special case dimensions independent
{
	vt=v1;
	vt.v[rank_d]=v2.v[rank_d];
	vt.v2[rank_d]=v2.v2[rank_d];
	return vt;
}

if (c1==0) return v2;
if (c2==0) return v1;


vt=v1;

if (c1+c2==0)
{
printf("\n ERROR coeff_times_vel2");
printf("c1,c2 %f %f",c1,c2);
infinite:; goto infinite;
}

n=(int)((c2/(c1+c2))*vt.size);
	
	rank_size=vt.size;
	for (i=0;i<rank_size;i++) rank[i]=i;
	
	// Choose at random n different components and set them to v2 values
	for (i=0;i<n;i++)
	{
		j=alea(0,rank_size-1);
		l=rank[j];
		
		vt.v[l]=v2.v[i];
		vt.v2[l]=v2.v2[i]; 
		
		// Compact rank[]
		if (j<rank_size-1)
		{
			for (k=j;k<rank_size-1;k++)
			{
				rank[k]=rank[k+1];
			}
			rank_size=rank_size-1;
		}
	}

return vt;
}


/*------------------------------------------------------------- POS_MINUS_POS */
struct velocity	pos_minus_pos(struct position p1,struct position p2,int trace)
	/* Applying v to p2 gives p1, or, in short, p2+v=p1, or v=p1-p2 
	
	*/
{
int				i;
struct velocity vt;

if (trace>2) 
{
	printf("\npos_minus_pos");
	display_position(p1);
	display_position(p2);
	
}

vt.size=p1.size;
vt.ek=0; // No meaning ...


for (i=0;i<vt.size;i++)
{

	vt.v[i]=p1.x[i];
	vt.v2[i]=AS;	
/*
	vt.v[i]=p2.x[i];
	vt.v2[i]=p1.x[i];
*/
/*	
	if (p1.x[i]==p2.x[i])
	{
		vt.v2[i]=NO;
		continue;
	}

	
	vt.v[i]=p1.x[i];
	vt.v2[i]=AS;
*/
	
}
vt.combin=1;
if (trace>2) display_velocity(vt);
return vt;

}
/*-------------------------------------------------------------------POS_PLUS_VEL */
struct position	pos_plus_vel (struct position p,struct velocity v,struct param param)
{
/* Modify a position according to a velocity, for combinatorial kit
	position + velocity => new position
WARNING: just _position_ is modified. Not kinetic energy nor f value.
They MUST be modified after that.
*/

int				i;
int				max,min;
struct position	pt;

pt=p;

if (param.printlevel>2)
{
	printf("\n pos_plus_vel");
	display_position(p);
	display_velocity(v);
}

for (i=0;i<v.size;i++)
{	
	if (v.v2[i]==NO) continue; // No modification
	if (v.v2[i]==AS) {pt.x[i]=v.v[i];continue;} // Assign
		
	/*
	if (v.v[i]==pt.x[i]) // Can be applied
	{
		pt.x[i]=v.v2[i];
		continue;
	}
	*/
	
	
	
	// All possible values are equivalent
	//if (v.v2[i]==NA)
	{ 
		min=(int)param.H.min[i]; // WARNING. We DO suppose coding values are all integers in [min, max]
		max=(int)param.H.max[i];
		pt.x[i]=(float)alea(min,max);
		continue;
	}
		
	
}
	
if (param.printlevel>2)
{
	display_position(pt);
}	
	
return pt;
}


/*------------------------------------------------------------------- 	VEL_PLUS_VEL */
struct velocity	vel_plus_vel(struct velocity v1,struct velocity v2,struct param param) 
{
int				i;
int				max;
int				min;
struct velocity vt;
float			z;

if (param.printlevel>2)
{
	printf("\n vel_plus_vel");
	display_velocity(v1);
	display_velocity(v2);
}
vt=v1;

for (i=0;i<v1.size;i++)
{
	if (vt.v2[i]==NA) continue; // random + v2 = random
	if (v2.v2[i]==NO) continue; // v + null = v
	if (v2.v2[i]==NA) {vt.v2[i]=NA;continue;} // v + random = random
	if (vt.v2[i]==NO) {vt.v[i]=v2.v[i];vt.v2[i]=v2.v2[i];continue;} // null + v2 = v2
}

/*
for (i=0;i<v1.size;i++)
{
	if (v2.v[i]==v1.v2[i]) // Can be applied (v1,v1_2) + (v2,v2_2) => (v1,v2_2)
	{
		vt.v2[i]=v2.v2[i];
		continue;
	}
	
	if (v2.v2[i]==NO) continue; // No modification

		// else, all values are equivalent
		min=(int)param.H.min[i]; // WARNING. We DO suppose coding values are all integers in [min, max]
		max=(int)param.H.max[i];
		vt.v2[i]=(float)alea(min,max);
}
*/
/*
for (i=0;i<v1.size;i++)
{
	z=alea_float(0,1);
	if (z<0.5) continue;
		vt.v[i]=v2.v[i];
		vt.v2[i]=v2.v2[i];
	
}

*/
if(param.printlevel>2)
{
	display_velocity(vt);
	printf("\n");
}

return vt;

}