/**********************************************************************/
/*                                                                    */
/*            (c) Copyright, 1997 by Professor Gabriel Robins         */
/*                                                                    */
/*      Department of Computer Science, University of Virginia        */
/*          Charlottesville, VA 22903-2442  (804) 982-2207            */
/*   robins@cs.virginia.edu    http://www.cs.virginia.edu/~robins/    */
/*                                                                    */
/*   This code may be freely used for all non-commercial purposes.    */
/*   All copies/portions of this code must contain this header.       */
/*                                                                    */
/**********************************************************************/

/**********************************************************************
**                     steiner5.c
**********************************************************************/

#include "geometry.h"

extern char **errorstrings;
extern void error(enum errors, char *);

void        find_neighbours(int);
void        closest_neighbour(void);
void	    update_tree();
void        generate_MST3(void);
void	    update_neighbour(FILE*,int quadrant, int dist,int j);

/*
 * debugging functions
 */
void 	     print_points(void);
void	     print_tree(FILE *, int);


#define	MAXINT	32767
#define FMT1	"number of points:\t%d\tmst_cost=\t%d\n"
#define DEBUGFMT1 "\n$$ beginning of linear_MST  %d $$\n"
#define DEBUGFMT2 "working on S.P. point X=%d, Y=%d\n"
#define DEBUGFMT3 "side %d: now on point %d, the max%d is %d, po%d is %d\n"



FILE 	*fp;
int 	num_candidates,id_num = 0;

/**********************************************************************
**	
**	This function finds the closest neighbour for a node in 	
**	its neighbour array
***********************************************************************/
void closest_neighbour()
{
	struct edge tmp;
	int i,
	    smallest_edge=0;


	for ( i = 0; i < 5; i++)
	        if (Neighbour[i].weight < Neighbour[smallest_edge].weight)
		        smallest_edge = i;
	if (smallest_edge != 0){
	        tmp = Neighbour[0];
	        Neighbour[0] = Neighbour[smallest_edge];
	        Neighbour[smallest_edge]= tmp;
	      }
}

/**********************************************************************
**	
**	This function updates a node's neighbour array 
**********************************************************************/

void	update_neighbour(FILE * fp,int quadrant, int dist,int j)
{
        Neighbour[quadrant].P2 = j;
	Neighbour[quadrant].weight = dist;

#ifdef DEBUG2
fprintf(fp,"\nupdate Neighbour[%d]= %d",quadrant, Neighbour[quadrant].weight);
#endif
}


/**********************************************************************
**	
**	This function finds up to four closest neighbours for a node 
**	in each quadrant
**********************************************************************/

void find_neighbours(int index)
{
        int i,j,
	    dx,dy,
	    dist,
	    quadrant;
	EDGE*	temp;
	   
	    

	for (i = 0; i < 5; i++)		/*initial the neighbour array */
	{
	   	Neighbour[i].P2 = NONE;
		Neighbour[i].weight= MAXINT;
	 }
  
	for (j=0;j<index;j++)
	{
	        dy= pointset[j].Y - pointset[index].Y;
		dx= pointset[j].X - pointset[index].X;

/*		if(dy>0) if(dx>0) if(dy<dx)   quadrant = 0; else quadrant = 1;
		         else     if(dy+dx>0) quadrant = 1; else quadrant = 2;
		else     if(dx<0) if(dy>dx)   quadrant = 2; else quadrant = 3;
		         else     if(dy+dx<0) quadrant = 3; else quadrant = 0;
*/
		
		if(dx > dy) if(dx + dy < 0)  quadrant = 3; else quadrant = 0;
		else     if(dx + dy <0)      quadrant = 2; else quadrant = 1;
         
		
		if (dy <0) dy = -dy;
		if (dx <0) dx = -dx;
		dist = dy + dx;
		
		temp = &Neighbour[quadrant];
		if (dist < temp ->weight){
	  	           temp ->P2 = j;
		           temp ->weight = dist;
		  }
				     
	 }
/*	closest_neighbour();*/
}



/***********************************************************************
**	
**	This function updates the tree by deleting the longest edge along
**	the cycle generated by adding the new point
**	called by  linear_MST
***********************************************************************/
void	update_Tree(int first, int last)
{
        int 	prev,curr,next;
	NODE	preNode,curNode;
  
#ifdef DEBUG
	FILE* fp;
	if (!(fp = fopen("c.update","a"))) {
	      fprintf(stderr,"***Error: cannot open file update for write.\n");
	      exit(-1);
	}
#endif

	if (first== last) return;

	if (Tree[first].parent == last){
              Tree[last]= Tree[first];
	      Tree[last].parent = first;
	      return;
	    }
	else {
	      prev = first;
	      curr= Tree[prev].parent;
	      next = Tree[curr].parent;
  
	      preNode = Tree[prev];
	      curNode = Tree[curr];
   
	      while (prev != last){
#ifdef DEBUG
fprintf(fp, "prev=%d, curr= %d, next= %d \n", prev, curr, next);
#endif
                      Tree[curr]= preNode;
                      Tree[curr].parent = prev;
                      preNode = curNode;
                      curNode = Tree[next];
                      prev= curr;
                      curr= next;
                      next= Tree[next].parent;

#ifdef DEBUG
fflush(fp);  
                      fclose(fp);
#endif
}
	    }
#ifdef DEBUG
fprintf(fp,"********end of function update_tree\n");
fflush(fp);  
fclose(fp);
#endif
}

/**********************************************************************
**
**	traverse the tree from a bottom node to a up level one
**********************************************************************/
void move_up(int *index, int *po)    
{
	if (Tree[*index].parent != NONE) {
		*index = Tree[*index].parent;
		Tree[*index].id = id_num;
		if (Tree[*index].weight > Tree[*po].weight)
			*po = *index;
	}
}

/**********************************************************************
**	trace back to the starting point
*
***********************************************************************/
int  reset_po(int first,int last)
{
        int i,tmp;
  
	tmp = first;
	for(i=first; i != last;i=Tree[i].parent)
	         if (Tree[i].weight > Tree[tmp].weight)
		        tmp = i;
	return tmp;
}

/**********************************************************************
**	
**	This function updates the MST by connecting the new node with
**	its up to four neighbours in the tree and deleting the longest 
**	edge along the induced cycle
**********************************************************************/
int linear_MST(int new_point)
{
	int cost=0, i = 1, j, index1, index2, dif, po1, po2, savings=0, PATH;
#ifdef DEBUG2
	FILE* fp;
	if (!(fp = fopen("c.out","a"))) error(EOPEN, "c.out");
fprintf(fp, DEBUGFMT1, new_point);
fprintf(fp,DEBUGFMT2,pointset[new_point].X,pointset[new_point].Y);
fprintf(fp,"the original cost is %d   \n", cost);
#endif

	Tree[new_point].parent = Neighbour[0].P2;
	if (new_point)		/* if not first point */	
  		Tree[new_point].weight = cost = Neighbour[0].weight;
	/*for (; i < 4 && Neighbour[i].P2 != NONE ; i++)*/
	for (; i < 4 ; i++) 
	  if ( Neighbour[i].P2 != NONE ) {
	        int parent1,
		    parent2,
		    weight1,
		    weight2,
		po2 = index2 = Neighbour[i].P2;
		Tree[index2].id = ++id_num;
		po1 = index1 = new_point;
		Tree[index1].id = id_num;
		dif = 0;
		PATH = 0;
		parent1 = Tree[index1].parent;
		parent2 = Tree[index2].parent;
		while (Tree[parent1].id != id_num
			&& Tree[parent2].id != id_num
			&& parent1 != parent2) {
			move_up(&index1, &po1);
			move_up(&index2, &po2);
#ifdef DEBUG2
fprintf(fp, DEBUGFMT3, 1, index1, 1, Tree[po1].weight, 1, po1);
fprintf(fp, DEBUGFMT3, 2, index2, 2, Tree[po2].weight, 2, po2);
#endif

			parent1 = Tree[index1].parent;	
			parent2 = Tree[index2].parent;
		}
		
		if (Tree[parent1].id == id_num) {
			if (parent1 == Neighbour[i].P2)
	    			PATH = 1;
			else po2 = reset_po(Neighbour[i].P2, parent1);
		}
		if (Tree[parent2].id == id_num) {
			if (parent2 == new_point)
				PATH = 2;
			else po1 = reset_po(new_point, parent2);
		}
#ifdef DEBUG2
fprintf(fp,"PATH %d, po1 is %d, po2 is %d\n", PATH, po1, po2);
#endif

		weight1 = Tree[po1].weight;
		weight2 = Tree[po2].weight;
		if (!PATH)
			PATH = (weight1 > weight2) ? 1 : 2;
		switch (PATH) {
			case 1:	if (weight1 > Neighbour[i].weight) {
					dif = weight1 - Neighbour[i].weight;
					update_Tree(new_point, po1);
					Tree[new_point].parent = Neighbour[i].P2;
					Tree[new_point].weight = Neighbour[i].weight;
				
    				}break;
			case 2:	if (weight2 > Neighbour[i].weight) {
					dif = weight2 - Neighbour[i].weight;
					update_Tree(Neighbour[i].P2, po2);
					Tree[Neighbour[i].P2].parent = new_point;
					Tree[Neighbour[i].P2].weight
						= Neighbour[i].weight;
				      }break;
			}
#ifdef DEBUG2
fprintf(fp,"****PATH %d, new_point %d, po1 %d, po2 %d, dif %d\n\n",
PATH,new_point,po1,po2,dif);
#endif
		savings += dif;
#ifdef DEBUG2
fprintf(fp, "------------Tree of this run on Neighbour %d----------\n",i);
print_tree(fp, new_point);
#endif

	}
#ifdef DEBUG2
fprintf(fp,"DONE checking neighbours on this S.P.  SAVE = %d \n", savings);
fclose(fp);
#endif

	return (savings - cost);
}	

/***************************************************************************/

void error(enum errors i, char *s)
{
	fprintf(stderr, errorstrings[i], s);
	exit(-1);
}

/**********************************************************************/

void print_points(void)
{
	int i;

	for(i = 0; i < number_of_points; i++) {
		printf("\n  i= %d pointset[i].X =%d, pointset[i].Y = %d\n",
			i,
			pointset[i].X,
			pointset[i].Y);
	}
}

/**********************************************************************/

void print_tree(FILE *fp, int new_point)
{
#ifdef DEBUG2
	int i;
	for (i = 0; i <= new_point; i++) {
		fprintf(fp,
			"tree[%d]  p=%d  weight= %d  \n",
             		i,
			Tree[i].parent,
			Tree[i].weight);
		fflush(fp);
	      }
#endif

}

/**********************************************************************/
/*	This function generates MST by adding the new nodes one at a time
**	calling linear_MST
**
***********************************************************************/

void	generate_MST3()
{
        int i,
	    j,
	    savings;

	Tree[0].weight = Tree[0].parent = Tree[0].id  = NONE;
	mst_cost = 0;
	for (j = 0; j < number_of_points; j++) {
	    find_neighbours(j);
            savings = linear_MST(j);
            mst_cost -= savings;
	    }
	/*printf("\nMST cost computed = %d\n", mst_cost);*/
	if(GRAPHICS == YES){
	    draw_point(pointset[0].X, pointset[0].Y, size);
	    for (i = 1; i < number_of_points; i++) {
		draw_point(pointset[i].X, pointset[i].Y, size);
		draw_edge(pointset[i], pointset[Tree[i].parent], 1);
	      }
	  }
}

/**********************************************************************
**	This function removes SPs of degree one and two
**
**********************************************************************/
void remove_invalid_SPs(int actual_number_of_points, int *SPs) {

int first_new_SP,i,j, k, prior_mst_cost;

number_of_points--;

    prior_mst_cost = mst_cost;

first_new_SP = number_of_points - *SPs;

for(j=actual_number_of_points;j<number_of_points;j++)
        degree_count[j] = 1;

for(j=1; j < number_of_points; j++)
        degree_count[Tree[j].parent]++;

for(i=number_of_points-1;i>=actual_number_of_points;i--)
if (degree_count[i] <= 2)
{
   if(i >= first_new_SP)
        *SPs = *SPs - 1;

   for(j=i;j<number_of_points-1;j++)
        pointset[j] = pointset[j+1];

   number_of_points--;
};

generate_MST3();
for(k=0;k<number_of_points;k++) old_Tree[k] = Tree[k];

number_of_points++;

};

/***********************************************************************
**	this function compute the saving  induced by the steiner points
**	using linear_MST updating scheme
***********************************************************************/

int generate_steiner5(void)
{
int   original_MST_cost,savings;
int     s, i, j, k, heuristic_num = 5, round = 0, nsp,
        actual_number_of_points, SPs_this_round;
int    num_candidates;

original_MST_cost = mst_cost;

for(i=0;i<number_of_points;i++)
    old_Tree[i] = Tree[i];
actual_number_of_points = number_of_points;
number_of_points++;

while(YES){
  num_candidates = 0;
  SPs_this_round = 0;
  for(i=0;i<actual_number_of_points;i++)
    for(j=0;j<actual_number_of_points;j++) if(i != j){
       pointset[number_of_points - 1].X = pointset[i].X;
       pointset[number_of_points - 1].Y = pointset[j].Y;

       find_neighbours(number_of_points - 1);
       if (Neighbour[0].weight != 0){
	 savings = linear_MST (number_of_points - 1);
	 if(savings > 0) {
           SP_candidates[num_candidates].X = pointset[number_of_points - 1].X;
	   SP_candidates[num_candidates].Y = pointset[number_of_points - 1].Y;
	   Steiner_savings[num_candidates] = savings;
	   num_candidates++;
         }
	 for(k=0;k<number_of_points;k++)
	   Tree[k] = old_Tree[k];
      }
    }	/* check all the Hanan points and set up the array for potential SPs*/

if(num_candidates == 0) break;

qsort_candidates(0, num_candidates-1);

for(k=0;k<num_candidates;k++)
   {
   pointset[number_of_points - 1] = SP_candidates[k];
  
   find_neighbours(number_of_points - 1);
   if (Neighbour[0].weight != 0){
	savings = linear_MST (number_of_points - 1);
        if(savings >= Steiner_savings[k])
	{
	  mst_cost -= savings;
	  number_of_points++;
	  SPs_this_round++;
	  for(i=0;i<number_of_points;i++)
	    old_Tree[i] = Tree[i];
	}
	else
	  for(i=0;i<number_of_points;i++)
	    Tree[i] = old_Tree[i];
      }
 }	/* add all the non-interferred SPs for this round */

round++;
remove_invalid_SPs(actual_number_of_points,&SPs_this_round);
if(SPs_this_round == 0) break;

}   

erase_window();
if(GRAPHICS == YES){
    draw_point(pointset[0].X, pointset[0].Y, size); 
    for (i = 1; i < number_of_points-1; i++){ 
      if (i < actual_number_of_points)	
	 draw_point(pointset[i].X, pointset[i].Y, size); 
      else
	draw_point(pointset[i].X, pointset[i].Y, size+2);

      draw_edge(pointset[i], pointset[Tree[i].parent],1);
    }
  }

steiner_cost[heuristic_num] = mst_cost;
nsp = number_of_points-1-actual_number_of_points;
printf("ST_cost  %d      MST_cost  %d \n", mst_cost, original_MST_cost);
printf("#SPs: %d\n", nsp);

if (OUT_TREE)
{
printf("#SPs: %d\n", nsp);

for (i = actual_number_of_points; i < number_of_points-1; i++) 
         printf("%d %d \n", pointset[i].X,pointset[i].Y);
		
}

improvement_percent[heuristic_num] = 
  (original_MST_cost - steiner_cost[heuristic_num]) * 100.0 / original_MST_cost;

num_steiner_pts[heuristic_num] = number_of_points - actual_number_of_points;
num_N_pts[heuristic_num][num_steiner_pts[heuristic_num]] ++;
ave_num_steiner_pts[heuristic_num] = ave_num_steiner_pts[heuristic_num] 
                                   + num_steiner_pts[heuristic_num];
min_num_steiner_pts[heuristic_num] = min(min_num_steiner_pts[heuristic_num],
                                         num_steiner_pts[heuristic_num]);
max_num_steiner_pts[heuristic_num] = max(max_num_steiner_pts[heuristic_num],
                                         num_steiner_pts[heuristic_num]);

min_round = min (round, min_round);
max_round = max (round, max_round);
total_round = total_round + round;

number_of_points = actual_number_of_points;

return(improvement_percent[heuristic_num]);
}

/**********************************************************************/