/*
 * The author of this software is Michael Trick.  Copyright (c) 1994 by 
 * Michael Trick.
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose without fee is hereby granted, provided that this entire notice
 * is included in all copies of any software which is or includes a copy
 * or modification of this software and in all copies of the supporting
 * documentation for such software.
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR DOES NOT MAKE ANY
 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
 */

/*
   COLOR.C: Easy code for graph coloring
   Author: Michael A. Trick, Carnegie Mellon University, trick+@cmu.edu
   Last Modified: November 2, 1994


Graph is input in a file.  First line contains the number of nodes and
edges.  All following contain the node numbers (from 1 to n) incident to 
each edge.  Sample:

4 4
1 2
2 3
3 4
1 4

represents a four node cycle graph.

Code is probably insufficiently debugged, but may be useful to some people.

For more information on this code, see Anuj Mehrotra and Michael A. Trick,
"A column generation approach to graph coloring", GSIA Technical report series.

Modified by Gang Chen <chen.gang1983@gmail.com> at Jun 25, 2010
  * use dynamic memory allocation to handle large networks
  * print coloring result to screen

*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/times.h>

#define MAX_RAND (2.0*(1 << 30))
#define MAX_NODE 40000
#define TRUE 1
#define FALSE 0
#define INF 100000.0

struct tms buffer;			/* structure for timing              */
int current_time,start_time;
double utime;

int **adj;
int BestColoring;
int num_node;
int *ColorClass;
int prob_count;
int *Order;
int *Handled;
int **ColorAdj;
int *ColorCount;
int lb;
int num_prob,max_prob;


int best_clique;


int greedy_clique(valid,clique)
int *valid,*clique;

{
  int i,j,k;
  int max;
  int place,done;
  int *order;
  int *weight;
  weight = (int *)malloc(sizeof(int)*num_node);
  
  for (i=0;i<num_node;i++) clique[i] = 0;
  order = (int *)calloc(num_node+1,sizeof(int));
  place = 0;
  for (i=0;i<num_node;i++) {
    if (valid[i]) {
      order[place] = i;
      place++;
    }
  }
  for (i=0;i<num_node;i++)
    weight[i] = 0;
  for (i=0;i<num_node;i++) 
    {
      if (!valid[i]) continue;
      for (j=0;j<num_node;j++) 
	{
	  if (!valid[j]) continue;
	  if (adj[i][j]) weight[i]++;
	}
    }
  

  done = FALSE;
  while (!done) {
    done = TRUE;
    for (i=0;i<place-1;i++) {
      j = order[i];
      k = order[i+1];
      if (weight[j] < weight[k]) {
	order[i] = k;
	order[i+1] = j;
	done = FALSE;
      }
    }
  }


  clique[order[0]] = TRUE;
  for (i=1;i<place;i++) {
    j = order[i];
    for (k=0;k<i;k++) {
      if (clique[order[k]] && !adj[j][order[k]]) break;
    }
    if (k==i) {
      clique[j] = TRUE;
    }
    else clique[j] = FALSE;
    
  }
  max = 0;
  for (i=0;i<place;i++) 
    if (clique[order[i]]) max ++;
  
  free(order);
/*  printf("Clique found of size %d\n",max);*/
  free(weight);
  
  return max;
}

/* Target is a goal value:  once a clique is found with value target
   it is possible to return
   
   Lower is a bound representing an already found clique:  once it is
   determined that no clique exists with value better than lower, it
   is permitted to return with a suboptimal clique.

   Note, to find a clique of value 1, it is not permitted to just set
   the lower to 1:  the recursion will not work.  Lower represents a
   value that is the goal for the recursion.
*/

int max_w_clique(valid,clique,lower,target)
int *valid,*clique;
int lower,target;


{
  int start,j,k;
  int incumb,new_weight;
  int *valid1,*clique1;
  int *order;
  int *value;
  int i,place,finish,done,place1,place2;
  int total_left;
  
/*  printf("entered with lower %d target %d\n",lower,target);*/
  num_prob++;
  if (num_prob > max_prob) return -1;
  for (j=0;j<num_node;j++) clique[j] = 0;
  total_left = 0;
  for (i=0;i<num_node;i++)
    if (valid[i]) total_left ++;
  if (total_left < lower) {
    return 0.0;
  }

  order = (int *)calloc(num_node+1,sizeof(int));
  value = (int *) calloc(num_node,sizeof(int));
  incumb = greedy_clique(valid,clique);
  if (incumb >=target) return incumb;
  if (incumb > best_clique) {
    best_clique=incumb;
/*    printf("Clique of size %5d found.\n",best_clique);*/
  }
/*  printf("Greedy gave %f\n",incumb);*/
  
  place = 0;
  for (i=0;i<num_node;i++) {
    if (clique[i]) {
      order[place] = i;
      total_left --;
      place++;
    }
  }
  start = place;
  for (i=0;i<num_node;i++) {
    if (!clique[i]&&valid[i]) {
      order[place] = i;
      place++;
    }
  }
  finish = place;
  for (place=start;place<finish;place++) {
    i = order[place];
    value[i] = 0;
    for (j=0;j<num_node;j++) {
      if (valid[j] && adj[i][j]) value[i]++;
    }
  }

  done = FALSE;
  while (!done) {
    done = TRUE;
    for (place=start;place<finish-1;place++) {
      i = order[place];
      j = order[place+1];
      if (value[i] < value[j] ) {
	order[place] = j;
	order[place+1] = i;
	done = FALSE;
      }
    }
  }
  free(value);
  for (place=start;place<finish;place++) {
    if (incumb + total_left < lower) {
      return 0;
    }
    
    j = order[place];
    total_left --;
    
    if (clique[j]) continue;
    
    valid1 = (int *)calloc(num_node,sizeof(int));
    clique1 = (int *)calloc(num_node,sizeof(int));
    for (place1=0;place1 < num_node;place1++) valid1[place1] = FALSE;
    for (place1=0;place1<place;place1++) {
      k = order[place1];
      if (valid[k] && (adj[j][k])){
	valid1[k] = TRUE;
      }
      else
	valid1[k] = FALSE;
    }
    new_weight = max_w_clique(valid1,clique1,incumb-1,target-1);
    if (new_weight+1 > incumb)  {
/*      printf("Taking new\n");*/
      incumb = new_weight+1;
      for (k=0;k<num_node;k++) clique[k] = clique1[k];
      clique[j] = TRUE;
      if (incumb > best_clique) {
	best_clique=incumb;
/*	printf("Clique of size %5d found.\n",best_clique);*/
      }
    }
    
  /*    else printf("Taking incumb\n");*/
    free(valid1);
    free(clique1);
    if (incumb >=target) break;
  }
  free(order);
  return(incumb);
}



AssignColor(node,color)
int node,color;

{
  int node1;

/*  printf("  %d color +%d\n",node,color);*/
  ColorClass[node] = color;
  for (node1=0;node1<num_node;node1++) 
    {
      if (node==node1) continue;
      if (adj[node][node1]) 
	{
	  if (ColorAdj[node1][color]==0) ColorCount[node1]++;
	  ColorAdj[node1][color]++;
	  ColorAdj[node1][0]--;
	  if (ColorAdj[node1][0] < 0) printf("ERROR on assign\n");	
	}
    }
  
}


RemoveColor(node,color)
int node,color;

{
  int node1;
/*  printf("  %d color -%d\n",node,color);  */
  ColorClass[node] = 0;
  for (node1=0;node1<num_node;node1++) 
    {
      if (node==node1) continue;
      if (adj[node][node1]) 
	{
	  ColorAdj[node1][color]--;
	  if (ColorAdj[node1][color]==0) ColorCount[node1]--;
	  if (ColorAdj[node1][color] < 0) printf("ERROR on assign\n");
	  ColorAdj[node1][0]++;
	}
    }
  
}

/*int lower_bound(current_color) 
     int current_color;

{
  int i,j,k;
  int max;
  int done;
  int clique[MAX_NODE];
  int order[MAX_NODE];
  int weight[MAX_NODE];
  
  for (i=0;i<num_node;i++) clique[i] = 0;
  for (i=0;i<num_node;i++) {
    order[i] = i;
    weight[i] = ColorAdj[i][0];
    for (k=1;k<=current_color;k++)
      if (ColorClass[i] != k) weight[i] +=ColorAdj[i][k];
  }
  
  done = FALSE;
  while (!done) {
    done = TRUE;
    for (i=0;i<num_node-1;i++) {
      j = order[i];
      k = order[i+1];
      if (weight[j] < weight[k]) {
	order[i] = k;
	order[i+1] = j;
	done = FALSE;
      }
    }
  }
  clique[order[0]] = TRUE;
  printf("Node %d in clique\n",order[0]);
  for (i=1;i<num_node;i++) {
    j = order[i];
    for (k=0;k<i;k++) {
      if (clique[order[k]]) {
	if (!adj[j][order[k]]) 
	  {
	    if ((ColorClass[order[k]]==0) || (ColorClass[j]==0) ||(ColorClass[order[k]]==ColorClass[j]))
	      break;
	  }
      }
    }
    if (k==i) {
      clique[j] = TRUE;
      printf("Node %d in clique\n",j);
    }
  }
  max = 0;
  for (i=0;i<num_node;i++) 
    if (clique[order[i]]) max++;
  printf("Lower bound of %d found with Best Coloring %d\n",max,BestColoring);
  return max;
}
*/

/*int lower_bound(current_color)
int current_color;
{
  return 0;
  
}
*/

int color(i,current_color)
     int i;
{
  int j,new_val;
  int k,max,count,place;
  

  prob_count++;
  if (current_color >=BestColoring) return(current_color);
  if (BestColoring <=lb) return(BestColoring);
  
  if (i >= num_node) return(current_color);
/*  printf("Node %d, num_color %d\n",i,current_color);*/
  
/* Find node with maximum color_adj */
  max = -1;
  place = -1;
  for(k=0;k<num_node;k++) 
    {
/*      count = 0;*/
      if (Handled[k]) continue;
/*      for (j=1;j<=current_color;j++)
	if (ColorAdj[k][j] > 0) count++;
      if (count!=ColorCount[k]) printf("Trouble with color count\n");
*/      
/*      printf("ColorCount[%3d] = %d\n",k,ColorCount[k]);*/
      if ((ColorCount[k] > max) || ((ColorCount[k]==max)&&(ColorAdj[k][0]>ColorAdj[place][0])))
	{
/*	  printf("Best now at %d\n",k);*/
	  max = ColorCount[k];
	  place = k;
	}
    }
  if (place==-1) 
    {
      printf("Graph is disconnected.  This code needs to be updated for that case.\n");
      exit(1);
    }

  
  Order[i] = place;
  Handled[place] = TRUE;
/*  printf("Using node %d at level %d\n",place,i);*/
  for (j=1;j<=current_color;j++) 
    {
      if (!ColorAdj[place][j]) 
	{
	  ColorClass[place] = j;
	  AssignColor(place,j);
	  new_val = color(i+1,current_color);
	  if (new_val < BestColoring){
	    BestColoring = new_val;
	    print_colors();
	  }
	  RemoveColor(place,j);
	  if (BestColoring<=current_color) {
	    Handled[place] = FALSE;
	    return(BestColoring);
	  }
	}
    }
  if (current_color+1 < BestColoring) 
    {
      ColorClass[place] = current_color+1;
      AssignColor(place,current_color+1);
      new_val = color(i+1,current_color+1);
      if (new_val < BestColoring) {
	BestColoring = new_val;
	print_colors();
      }
      
      RemoveColor(place,current_color+1);
    }
  Handled[place] = FALSE;
  return(BestColoring);
}

print_colors() 
{
  int i,j;

  times(&buffer);
  current_time = buffer.tms_utime;
  
  printf("Best coloring is %d at time %7.1f\n",BestColoring,(current_time-start_time)/60.0);
  
  for (i=0;i<num_node;i++)
    printf("%d = %d\n",i,ColorClass[i]);
  for (i=0;i<num_node;i++)
    for (j=0;j<num_node;j++) 
      {
	if (i==j) continue;
	if (adj[i][j] && (ColorClass[i]==ColorClass[j]))
	  printf("Error with nodes %d and %d and color %d\n",i,j,ColorClass[i]);
      }
}

main(argc,argv)
  int argc; 
  char *argv[];
{
  FILE *fp;
  int m,i,j,k,val;
  int *valid,*clique;
  int place;



  if (argc < 2) 
    {
      printf("Provide data file in command.\n");
      exit(1);
    }
  
  if ((fp = fopen(argv[1],"r"))==NULL) 
    {
      printf("No such graph as %s\n",argv[1]);
      exit(1);
    }
  
  fscanf(fp,"%d %d",&num_node,&m);

  valid = (int *)malloc(sizeof(int)*num_node);
  clique = (int *)malloc(sizeof(int)*num_node);

  adj = (int **)malloc(sizeof(int *)*num_node);
  if(adj == NULL) return(1);
  for (i=0; i<num_node; i++){
  	adj[i] = (int *)malloc(sizeof(int)*num_node);
  	if(adj[i] == NULL) return(1);
  }

  ColorClass = (int *)malloc(sizeof(int)*num_node);
  Order = (int *)malloc(sizeof(int)*num_node);
  Handled = (int *)malloc(sizeof(int)*num_node);
  
  ColorAdj = (int **)malloc(sizeof(int *)*num_node);
  for (i=0; i<num_node; i++){
  	ColorAdj[i] = (int *)malloc(sizeof(int)*num_node);
  }

  ColorCount = (int *)malloc(sizeof(int)*num_node);

  for (i=0;i<num_node;i++)
    for (j=0;j<num_node;j++) 
		adj[i][j] = FALSE;

  for (k=0;k<m;k++) {
    fscanf(fp,"%d %d",&i,&j);
    adj[i-1][j-1] = TRUE;
    adj[j-1][i-1] = TRUE;
  }
  fclose(fp);


  prob_count = 0;
  for (i=0;i<num_node;i++)
    for (j=0;j<=num_node;j++)
      ColorAdj[i][j] = 0;
  for (i=0;i<num_node;i++)
    for (j=0;j<num_node;j++)
      if (adj[i][j]) ColorAdj[i][0]++;
  for (i=0;i<num_node;i++)
    ColorCount[i]=0;
  times(&buffer);
  start_time=buffer.tms_utime;
  printf("Graph %s read with %d nodes and %d edges\n",argv[1],num_node,m);
  for (i=0;i<num_node;i++)
    Handled[i] = FALSE;
  BestColoring = num_node+1;
/*  ColorClass[0] = 1;
  AssignColor(0,1);
  Handled[0] = TRUE;*/
  for (i=0;i<num_node;i++) valid[i] = TRUE;
  best_clique = 0;
  num_prob = 0;
  max_prob = 10000;
  
  lb = max_w_clique(valid,clique,0,num_node);
  place = 0;
  for (i=0;i<num_node;i++) 
    {
      if (clique[i]) 
	{
	  Order[place] = i;
	  Handled[i] = TRUE;
	  place++;
	  AssignColor(i,place);
	  for (j=0;j<num_node;j++)
	    if ((i!=j)&&clique[j] && (!adj[i][j])) printf("Result is not a clique!\n");
	  
	}
    }
  
  printf("Lower bound is %d",lb);
  if (num_prob >=max_prob) printf(" (not confirmed)\n");
  else printf("\n");
  val = color(place,place);
  times(&buffer);
  current_time=buffer.tms_utime;
  
  printf("Best coloring has value %d, subproblems: %d time:%7.1f\n",val,prob_count,(current_time-start_time)/60.0);

}