/**********************************************************************/
/*                                                                    */
/*            (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.       */
/*                                                                    */
/**********************************************************************/

/**********************************************************************/
/*                      Calculations.c                                */
/**********************************************************************/

#include "geometry.h"


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

/*char *errorstrings[] = {
	"***** Error: Storage allocation for %s\n",
	"***** Error: Cannot open file %s\n",
};*/

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

void	generate_random_points();
void	print_separator();
void	do_randomize();

void multi_loop_stats(void)
{
int i, j, old_gridsize, old_nop;

old_gridsize = grid_size;
old_nop = number_of_points;

for(i=3;i<=50;i++)
   {
   for(j=0;j<9;j++)
      {
      number_of_points = i;
      grid_size = grid_sizes[j];
      loop_stats();
      }
   print_separator();
   }
   
grid_size = old_gridsize;
number_of_points = old_nop;

}

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

void loop_stats(void)
{
int i, span_cost, stei_cost, mst_total, steiner_total;
float total_perf;

/*printf("Seed = %u\n", (int) next); fflush(stdout);*/

mst_total =  steiner_total = 0;
total_perf = 0.0;
for(i=0; i< number_of_iterations;i++)    {
   erase_window();
   if(choice !=2) do_randomize();	
   generate_MST3(number_of_points); /* getchar();*/
   span_cost = mst_cost;
   mst_total = mst_total + span_cost;
   generate_steiner5();  getchar();
   stei_cost = mst_cost; /*steiner_cost[4];*/
   steiner_total = steiner_total + stei_cost;

   total_perf = total_perf +
     (float)(span_cost-stei_cost)/(float)span_cost;
 }
}


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

int max(int i, int j) { if (i > j) return(i); else return(j); }
int min(int i, int j) { if (i < j) return(i); else return(j); }
double fmax(double i, double j) { if (i > j) return(i); else return(j); }
double fmin(double i, double j) { if (i < j) return(i); else return(j); }
int sig(int i){ if (i>0) return(1); else if (i<0) return(-1); else return(0); }

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

int generate_unique_identifier(void)
{
return(unique_identifier++);
}

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

void draw_left_half_L(struct point p1, struct point p2, int thickness)
{
struct point p3, p4;

  p3.X = min(p1.X, p2.X); p3.Y = min(p1.Y, p2.Y);
  p4.X = min(p1.X, p2.X); p4.Y = max(p1.Y, p2.Y);
  draw_line(p3, p4, thickness);
}

void draw_right_half_L(struct point p1, struct point p2, int thickness)
{
struct point p3, p4;

  p3.X = max(p1.X, p2.X); p3.Y = min(p1.Y, p2.Y);
  p4.X = max(p1.X, p2.X); p4.Y = max(p1.Y, p2.Y);
  draw_line(p3, p4, thickness);
}

void draw_top_half_L(struct point p1, struct point p2, int thickness)
{
struct point p3, p4;

  p3.X = min(p1.X, p2.X); p3.Y = max(p1.Y, p2.Y);
  p4.X = max(p1.X, p2.X); p4.Y = max(p1.Y, p2.Y);
  draw_line(p3, p4, thickness);
}

void draw_bottom_half_L(struct point p1, struct point p2, int thickness)
{
struct point p3, p4;

  p3.X = min(p1.X, p2.X); p3.Y = min(p1.Y, p2.Y);
  p4.X = max(p1.X, p2.X); p4.Y = min(p1.Y, p2.Y);
  draw_line(p3, p4, thickness);
}

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

void draw_rectilinear_edge(struct point p1, struct point p2,
       int thickness, int L_orientation)
{
/*#ifdef MAC*/
int rnd_orientation[4] = { high_L, low_L, left_L, right_L };

switch (L_orientation)
  {
  case left_L:
  draw_left_half_L(p1, p2, thickness);
  if( sig(p1.X - p2.X) == sig(p1.Y - p2.Y) )
    draw_top_half_L(p1, p2, thickness);
  else draw_bottom_half_L(p1, p2, thickness);
  break;
 case right_L:
  draw_right_half_L(p1, p2, thickness);
  if( sig(p1.X - p2.X) == sig(p1.Y - p2.Y) )
    draw_bottom_half_L(p1, p2, thickness);
  else draw_top_half_L(p1, p2, thickness);
  break;
 case low_L:
  draw_bottom_half_L(p1, p2, thickness);
  if( sig(p1.X - p2.X) == sig(p1.Y - p2.Y) )
    draw_right_half_L(p1, p2, thickness);
  else draw_left_half_L(p1, p2, thickness);
  break;
 case high_L:
  draw_top_half_L(p1, p2, thickness);
  if( sig(p1.X - p2.X) == sig(p1.Y - p2.Y) )
    draw_left_half_L(p1, p2, thickness);
  else draw_right_half_L(p1, p2, thickness);
  break;
 case random_L:
  draw_rectilinear_edge(p1, p2, thickness, rnd_orientation[random(0,3)]);
  break;
}

/*#endif*/

}

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

void draw_edge(struct point p1, struct point p2, int thickness)
{

/*#ifdef MAC*/

switch (metric_type)
  {
  case Euclidean :
    draw_line(p1,p2,thickness);
    break;

  case Manhattan:
    draw_rectilinear_edge(p1, p2, thickness, L_edge_orientation);
    break;

  case Linfinity:
    draw_rectilinear_edge(p1, p2, thickness, L_edge_orientation);
    break;
  }

/*#endif*/

}

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

void draw_points(void)
{
/*#ifdef MAC*/
int i;

for(i=0;i<number_of_points;i++) 
   draw_point(pointset[i].X, pointset[i].Y, size);
   
/*#endif*/

}

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

int geom_dist(int i, int j)
{
int dist;
double dx, dy, rdist;

switch (metric_type)
  {
  case Manhattan:
         dist = abs(pointset[i].X - pointset[j].X)
                    + abs(pointset[i].Y - pointset[j].Y);
     break;

  case Euclidean:
     dx = pointset[i].X - pointset[j].X;
     dy = pointset[i].Y - pointset[j].Y;
     rdist = 10.0 * sqrt(dx*dx + dy*dy) + 0.5;
     dist = (int) rdist;
  break;

  case Linfinity:
           dist = max(abs(pointset[i].X - pointset[j].X),
                      abs(pointset[i].Y - pointset[j].Y));
     break;
  }
return(dist);
}

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

void generate_edges(void)
{
int i, j;

number_of_edges = 0;
for(i=0;i<number_of_points-1;i++)
  for(j=i+1;j<number_of_points;j++)
      {
      edges[number_of_edges].P1 = i;
      edges[number_of_edges].P2 = j;
      number_of_edges ++;
      }
}

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

void sort_edges(void)
{
qsort_edges(0, number_of_edges - 1);
}

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

void qsort_edges(int left, int right)
{
struct edge tmp;
int i, j, last, mid, smallest_edge, smallest_dist;

if(left>= right) return;

if(left>= (right-5))
   {
   for(i=left; i < right; i++)
      {
      smallest_dist = edges[i].weight; smallest_edge = i;
      for(j=i + 1; j <= right; j++)
         {
         if(edges[j].weight < smallest_dist)
            { smallest_dist = edges[j].weight; smallest_edge = j; }
         }
      tmp = edges[i]; edges[i] = edges[smallest_edge]; edges[smallest_edge] = tmp;
      }
   return;
   }
mid = (left + right) / 2;
tmp = edges[left]; edges[left] = edges[mid]; edges[mid] = tmp;
last = left;
for(i=left+1; i <= right; i++)
  if (edges[i].weight  < edges[left].weight)
     {
     last++;
     tmp = edges[last]; edges[last] = edges[i]; edges[i] = tmp;
     }
tmp = edges[last]; edges[last] = edges[left]; edges[left] = tmp;
qsort_edges(left, last - 1);
qsort_edges(last + 1, right);
}

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

void draw_edges(void)
{
/*#ifdef MAC*/

int  j;

for(j=0;j<number_of_edges;j++)
if(edges[j].flag == MARKED)
   draw_edge(pointset[edges[j].P1], pointset[edges[j].P2], 1);
 
/*#endif*/

}

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

void do_randomize(void)
{
generate_random_points();
/*redraw_points();*/
}

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

void print_separator(void)
{
printf("==============================================================\n");
}

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


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

 The following routine initializes a couple of global variables.
 GRAPHICS is used to display graphically the execution results of the code.
 
 **********************************************************************/

void init_vars(void)
{
   int i;

   metric_type = Manhattan;

   /*printf("Allocating memory for arrays.\n");*/

   pointset = (struct point *) calloc((max_number_of_points*4),
				      sizeof(struct point));
   SP_candidates = (struct point *) calloc((max_number_of_points *
					    max_number_of_points), 
					   sizeof(struct point));
   edges = (struct edge *) calloc((max_number_of_points *
				   max_number_of_points / 2), 
				  sizeof(struct edge));
        
   if((pointset == NULL) || (edges == NULL) ||
      (num_N_pts == NULL) || (SP_candidates == NULL))
      printf("*** Error: A:Storage allocation problem.\n");

   for (i = 0; i<HEURISTICS;i++)
   {
      num_N_pts[i] = (int *) calloc(max_number_of_points, sizeof(int));
      if(num_N_pts[i] == NULL)
	 printf("*** Error: B:Storage allocation problem.\n");
   }

   hit_list = (int *) calloc((max_number_of_points), sizeof(int));
   degree_count = (int *) calloc((max_number_of_points), sizeof(int));
   Steiner_savings = (int *) calloc((max_number_of_points *
				     max_number_of_points), sizeof(int));

   if((hit_list == NULL) || (degree_count == NULL) ||
      (Steiner_savings == NULL))
      printf("*** Error: C:Storage allocation problem.\n");

   if (!(Tree = (NODE*) calloc(1 + max_number_of_points,sizeof(NODE))))
      error(ECALLOC, "Tree");
Tree++;
Tree[NONE].id = Tree[NONE].parent = Tree[NONE].weight = NONE;
Tree[0].weight = Tree[0].parent = Tree[0].id  = NONE;
if (!(old_Tree = (NODE*) calloc(1 + max_number_of_points,sizeof(NODE))))
	error(ECALLOC, "old_Tree");
old_Tree++;
old_Tree[NONE].id = old_Tree[NONE].parent = old_Tree[NONE].weight = NONE;
old_Tree[0].weight = old_Tree[0].parent = old_Tree[0].id  = NONE;
if (!(Neighbour = (EDGE*) calloc(5, sizeof(EDGE))))
	error(ECALLOC, "Neighbour");
for(i=0;i<HEURISTICS;i++)
   {
   num_steiner_pts[i]=0;
   min_num_steiner_pts[i]=MAXINT;
   max_num_steiner_pts[i]=0;
   }
}

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

void generate_MST2(void)
{
int i, k, tmp, e1, e2;
struct edge e;

generate_edges();

/* for(i=0;i<number_of_edges;i++) edges[i].weight = geom_dist(edges[i].P1, edges[i].P2);
*/

for(i=0;i<number_of_edges;i++) 
if(metric_type == Manhattan)
  { e = edges[i]; e1 = e.P1; e2 = e.P2;
    tmp = pointset[e1].X - pointset[e2].X;
    if(tmp < 0) tmp = 0 - tmp;
    edges[i].weight = tmp;
    tmp = pointset[e1].Y - pointset[e2].Y;
    if(tmp < 0) tmp = 0 - tmp;
    edges[i].weight = edges[i].weight + tmp;
  }
else edges[i].weight = geom_dist(edges[i].P1, edges[i].P2);

sort_edges();

for(i=0;i<number_of_edges;i++) edges[i].flag = NOT_MARKED;
for(i=0;i<number_of_points;i++) hit_list[i] = NOT_MARKED;

hit_list[edges[0].P1] = MARKED;
for(i=1;i<number_of_points;i++)
  {
  k=0;
  while (hit_list[edges[k].P1] == hit_list[edges[k].P2]) k++;
  edges[k].flag = MARKED;
  hit_list[edges[k].P1] = MARKED;
  hit_list[edges[k].P2] = MARKED;
  }

mst_cost = 0;
for(i=0;i<number_of_edges;i++) 
   if(edges[i].flag == MARKED)
      mst_cost = mst_cost + edges[i].weight;

}

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

     This function actually generates the random points.

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


void generate_random_points(void)
{
int i;

for(i=0;i<number_of_points;i++)
  {
  pointset[i].X = random(1, grid_size);
  pointset[i].Y = random(1, grid_size);

  pointset[i].flag = NO;
  pointset[i].id = generate_unique_identifier();
  }

}

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

     This function "randomizes" the random number seed by setting it
     to the some function of the CPU clock and the real time.

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

void randomize_seed(void)
{
time_t current_time;

initial_seed = (unsigned long int) clock() + 
               (unsigned long int) time(&current_time);
next = initial_seed;
}

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

     This function generates a random integer within the given range,
     using a congruential multiplier algorithm.

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

int random(int i, int j)
{
int k;

next = next * 1103515245 + 54321;
k = (next/65536 % 32768) * (j-i) / 32768 + i;

return(min(max(k,i),j));
}