/* 
 * Search for witnesses to R(5,5) > k
 * M. Pleszkoch & H. Motteler, 10/25/92
 *
 * command line arguments
 *	S <n>	search <n> graphs
 *	M <n>	explore mutants to a depth <n>
 *	G <c>	select initial graph type <c>
 *	d <f>	set debugging flags to <f>
 *	s <n>	set random number generator's seed to <n>
 *
 * initial graph types
 *	r	random initital values for G
 *	i	read adjacency representation of G from stdin
 *	1	41 node symmetric G
 *
 * debug flags
 *      1	misc. initialization printing
 *	2	print initial classification vector
 *	4	print first part of initial 'badness' array
 *	8	print state info at outer iterations
 *	16	print state info at inner iterations
 *	32	perform some consistency checks
 */

#include <mpl.h>
#include <mp_libc.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <reduce.h> 

#define N 41		/* number of vertices */

/* the following types are for the MasPar model MP-1
 */
#define BOOL	char
#define VERT	char
#define INT8	char
#define INT16	short
#define INT32	int

#define maxval(x) ((1<<(sizeof(x)*8-1))-1)
#define acopy(A,B,L) {int i; for(i=0; i<L; i++) B[i] = A[i];}

/* Singular Data (variables allocated on ACU)
 */
INT32 T[34];		/* G's current classification vector */
INT32 Tnew[34];		/* temporary classification vector */
INT32 debug=0;		/* debugging flags */
INT32 gfnx='r'; 	/* code passed to graph generator */
INT32 genmat=0;		/* flag to generate NN training data */
INT32 search=0;		/* flag to enter search loop */
INT32 limit=1000;	/* default number of graphs to check */
INT32 num5sets;		/* number of distinct 5-groups */

/* Plural Data (variables allocated on each PE)
 */
plural BOOL G[N][N];	/* adjacency rep. of current graph */
plural BOOL A[5][5];	/* adj. mat. for local 5-group */
plural VERT v1, v2, v3;	/* vertices of the local 3-clusters */
plural BOOL in3set;	/* flag processors in 3-clusters */
plural VERT u, v;  	/* current edge of interest */
plural INT8 class;	/* classification of a 5-group */
plural INT8 class2;
plural BOOL in5set;
plural INT16 B;
plural BOOL inB;	/* NxN subarray; 'bad' edge count */

/* Classification Function
 *
 * Apply Mark's classification function to the set of 5-clusters
 * defined by the 3-clusters v1, v2, v3, and the edge (u,v).
 * The (N choose 3) 3-clusters are distributed over the PE array.
 * Classify() assumes that v1, v2, v3, u, and v are all distinct.
 * Classify() returns a classification code in the range 0-33
 * for each of the (N choose 3) 3-clusters for current edge (u, v).
 */
plural INT8 classify(u, v, edge) 
   plural VERT u, v; plural BOOL edge;
{
   extern void exit();
   INT32 i, j;
   plural INT8 t;

   plural VERT V[5];		/* current 5-group in sorted order */
   plural INT8 D[5];		/* out degrees of current 5-group */
   plural INT8 flip, ecnt, class;

   /* update local adjacency matrix */

   A[0][0] = A[1][1] = A[2][2] = A[3][3] = A[4][4] = 0;
   A[0][1] = G[v1][v2];
   A[0][2] = G[v1][v3];
   A[1][0] = G[v2][v1];
   A[1][2] = G[v2][v3];
   A[2][0] = G[v3][v1];
   A[2][1] = G[v3][v2];

   A[0][3] = G[v1][u];
   A[0][4] = G[v1][v];
  
   A[1][3] = G[v2][u];
   A[1][4] = G[v2][v];
 
   A[2][3] = G[v3][u];
   A[2][4] = G[v3][v];

   A[3][0] = G[v1][u];
   A[3][1] = G[v2][u];
   A[3][2] = G[v3][u];
   A[3][4] = edge;
 
   A[4][0] = G[v1][v];
   A[4][1] = G[v2][v];
   A[4][2] = G[v3][v];
   A[4][3] = edge;

   /* find local adjacencies and local edge count */
   ecnt = 0;
   for (i=0; i<5; i++) {	/* 15% time this loop */
      V[i] = i;
      D[i] = 0;
      for (j=0; j<5; j++)
	 D[i] += A[i][j];
      ecnt += D[i];
      }
   ecnt = ecnt >> 1;

   /* use G or ~G, whichever has fewer edges */
   if (ecnt <= 5)
      flip = 1;
   else {
      ecnt = 10-ecnt;
      for (i=0; i<5; i++) 
	 D[i] = 4 - D[i];
      flip = 0;
      }

   /* sort by out-degree and permute V[] same way */
   for (i=0; i<4; i++)
      for (j=i+1; j<5; j++)			/* 10% time this loop */
	 if (D[j] > D[i]) {
	    t = D[i]; D[i] = D[j]; D[j] = t;
	    t = V[i]; V[i] = V[j]; V[j] = t;
	    }

   /* Mark's Classification Formula: */
   class = D[0] * (D[0]-2) + 2 * (D[1]+D[3]) + D[4];

   if (ecnt == 4 && class == 7) 
      class += 4 + (A[V[3]][V[4]] == flip);

   if (ecnt == 5) {
      class += 4;
      if (class == 16 && D[1] == 2) 
	 class +=  3 * (A[V[0]][V[4]] == flip) - 1;
      }

   if (flip == 0)
      class = 33 - class;

   if (class < 0 || 33 < class) {
      /* let's just look at one of the errors... */
      if (enumerate() == 0) {
	 p_fprintf(stderr, "Classify(): illegal class code %d\n", class);
	 p_fprintf(stderr, "v1=%d v2=%d v3=%d u=%d v=%d\n",v1,v2,v3,u,v);
	 { int i, j;
	   for (i=0; i<5; i++) {
	      for (j=0; j<5; j++) 
		 p_fprintf(stderr, "%d ", A[i][j]);
	      p_fprintf(stderr, "\n");
	      }
	   }
	 }
	exit(127);
	}
   return class;
   }


/* Incremental Classification 
 *
 * Compute the classification vector Tnew[] that would result
 * if the edge (u,v) was flipped.  G is not modified.
 * Assumes v1 < v2 < v3 and u < v.
 */
void newclass(u, v)
   plural VERT u, v;
{
   static INT32 i;
   plural static INT8 T_p[34];
   extern void check_cv();

   for (i=0; i<34; i++) T_p[i] = 0;		/* 5% time here !? */

   /* compute change in classification on each processor */
   in5set = 
	(u!= v1 && u != v2 && u != v3 && v != v1 && v != v2 && v != v3);
   if (in5set) {

      class = classify(u, v,  G[u][v]);
      T_p[class]--;

      class2 = classify(u, v, ! G[u][v]);
      T_p[class2]++;
      
      /* compute new classification vector from individual changes */
      for (i=0; i<34; i++) {			/* 16% time this loop */
	 Tnew[i] = T[i];
	 Tnew[i] += reduceAdd16(T_p[i]);  /* (reduceAdd8 won't work) */
	 }
      }
   if (debug & 32) check_cv(Tnew);

   }


/* Maintain 'badness count' B for edges of G
 *
 * UpdateB() increments or decrementsedge counters in B for vertex 
 * pairs from vertices v1, v2, v3, u, v.  (A 5-clique in G or ~G.)
 * UpdateB() assumes that u < v and v1 < v2 < v3.
 * B is kept as a submatrix of the processor array, in the upper
 * LHS corner, and only the upper triangular part of B is used.
 */
void updateB(u, v, d)
   plural VERT u, v; plural INT8 d;
{
   plural INT8 i, j, k;
   plural INT8 a[4], b[3], c[5];
   plural INT16 dest;
   plural INT16 Btmp;
   plural INT16 t;

   /* first merge u,v with v1, v2, v3 to get a
    * sorted list of these vertices in c[0]-c[4]
    */
   a[0] = v1; a[1] = v2; a[2] = v3; a[3] = maxval(a[3]);
   b[0] = u; b[1] = v; b[2] = maxval(b[2]);
   for(i=0,j=0,k=0; k<5; k++) 
      if (a[i] < b[j]) 
	 c[k] = a[i++];
      else
	 c[k] = b[j++];
/*
   if (enumerate() == 0) {
	p_printf("v1=%d, v2=%d, v3=%d, u=,%d, v=%d\n", v1,v2,v3,u,v);
	p_printf("%d %d %d %d %d\n", c[0],c[1],c[2],c[3],c[4]);
	}
*/
   /* update B[vi][vj] for vi < vj 
    */
   for (i=0; i<4; i++)
      for (j=i+1; j<5; j++) {
	 dest = nxproc*c[i]+c[j];		/* row order */
	 t = router[dest].B;			/* get B[vi][vj] */
	 all if (inB) Btmp = 0;
	 p_sendwithAdd16( d, dest, &Btmp );	/* update B[vi][vj] */
	 all if (inB) B += Btmp;
/*	 router[dest].B = t+d;	/* TEMP: won't handle collisions!!! */
	 }
   }


/* Return the edge with the k-th largest 'bad' count
 * k=0 => return baddest, k=1 => return 2nd baddest, etc.
 * A sort is only done when k=0.
 */
void worstB(w1, w2, k)
   VERT *w1, *w2; INT16 k;
{
   plural static INT16 t;
   if (inB) {
      if (k == 0) t = rank16(B);
      if (t == (N*(N-1)/2-1) - k)  {
	 *w1 = globalor(iyproc);	/* row*/
	 *w2 = globalor(ixproc); 	/* column */
         }
      }
   }

/* pseudo random numbers 
 */
int seed = 19;
#define P1 127
#define P2 8191
int myrand(max) int max; {
   seed = (seed * P1) % P2;
   return seed % max;
   }

/* print a classification vector in 2 columns
 */
void print_cv(T) 
   int T[];
{
   int i;
   printf("Classification Vector:\n");
   for (i=0; i<17; i++) {
      printf("%2d=%5d, %2d=%5d\n", i, T[i], 33-i, T[33-i]);
      }
   printf("\n");
   }

/* perform a consistency check on a classsification vector
 */
void check_cv(T)
   INT32 T[];
{
   INT32 t; INT8 i;
   extern void exit();
   for (t=0,i=0; i<34; i++) t += T[i];
   if (t != num5sets) {
      printf("newclass(): bad classification vector\n");
      print_cv(T);
      exit(127);
      }
   }

/* Print mutation distance + classification vector to channel 'fp'
 * in a matlab readable form (i.e., a 'flat file').
 */
void fprint_cv2(fp, d) FILE *fp; int d; {
   int i;
   fprintf(fp, "%d ", d);
   for (i=0; i<34; i++) fprintf(fp, "%d ", T[i]);
   fprintf(fp, "\n");
   }

/* Initialize the graph G
 *   fnx code	function
 *	i	read G from stdin
 *	r	random
 *	1	41 node symmetric
 *	2	not implemented
 */
void initial_G(fnx) int fnx; {
  extern void exit();
  int i, j, t;
  INT8 g1[41] = {0,2,1,1,2,1,2,1,2,1,1,2,2,1,2,1,1,1,2,2,2, 
		    2,2,2,1,1,1,2,1,2,2,1,1,2,1,2,1,2,1,1,2};
   if (fnx == 'i') {
      if (debug & 1) printf("reading initial graph from stdin\n");
      for (i=0; i<N; i++)
	 for (j=0; j<N; j++) {
	    scanf("%d", &t);
	    G[i][j] = (plural BOOL) t;
	    }
      }
   else if (fnx == 'r') {
      if (debug & 1) printf("generating random initial graph\n");
      for (i=0; i<(N-1); i++)
	 for (j=i+1; j<N; j++)
	       G[j][i] = G[i][j] = (plural BOOL) myrand(2);
      }
   else if (fnx == '1') {
      if (debug & 1) printf("41 node symmetric initial graph\n");
      for (i=0; i<(N-1); i++)
	 for (j=i+1; j<N; j++)
	       G[j][i] = G[i][j] = (plural BOOL) g1[j-i]-1;
      }
   else {
      fprintf(stderr, "initial_G(): bad function code\n");
      exit(127);
      }
  }

/* print adjacency representation of G 
 */
void write_G(s)
   char *s;
{
   int i, j;
   FILE * fp;
   fp = fopen(s, "w");
   for (i=0; i<N; i++) {
      for (j=0; j<N; j++) 
	 fprintf(fp, " %d", proc[0].G[i][j]);
      fprintf(fp, "\n");
      }
   fclose(fp);
   } 

/* return 'cost' of a classification vector T[]
 */
float costfnx(T) 
   int T[];
{
   int i;
   float t;

   if (T[0] == 0 && T[33] == 0) return (float) 0;

   return
      (float)(T[0]+T[33]) / 10.0
    + (float)abs(T[1]+T[32]-10086) / 2000.0
    + (float)(T[2]+T[3]+T[30]+T[31]-51250) / 2000.0
    + (float)(T[4]+T[5]+T[6]+T[7]+T[26]+T[27]+T[28]+T[29]-160146) / 2000.0;

/*

   return (float) (T[0]+T[33]) / 10.0;

   return f_sqrt((float) T[0]*T[0] + T[33]*T[33]) / 10.0
          +  (float) abs((T[1]+T[32]) - 10000) / 500.0
	  +  (float) (T[2]+T[3]+T[30]+T[31]) / 51000.0;
	 
   return   (T[0]+T[33])/10.0  +  (T[1]+T[32])/10000.0  +
	  (T[2] + T[3] + T[30] + T[31])/50000.0;
*/
   }


void flipedge(i,j) 
   INT8 i, j;
{
   G[i][j] = ! G[i][j];
   G[j][i] = ! G[j][i];
   }

/* print current state in a compact form
 * k = graph num, (v1, v2) = edge, c = cost, T = class vec
 */
void print_stat(k, v1, v2, c, T)
int k, v1, v2; 
float c;
int T[];
{
/* 
   printf("%2d, (%2d,%2d), cost = %6.3f, class = [%3d, %3d, %3d, %3d, %3d]\n",
	  k, v1, v2, c, 
	  T[0], T[33], T[1], T[32], 
	  T[2] + T[3] + T[30] + T[31]);
*/
   printf("%2d, (%2d,%2d), cost = %6.3f, class = [%3d, %3d, %3d, %3d]\n",
	  k, v1, v2, c, 
	  T[0] + T[33], 
	  T[1] + T[32], 
	  T[2]+T[3] + T[30]+T[31],
	  T[4]+T[5]+T[6]+T[7] + T[26]+T[27]+T[28]+T[29] );
   }


/*	Action starts here
 */
main (argc, argv) 
   int argc; char *argv[];
{
   extern void exit();
   extern int getopt();
   extern void dump_weights();
   int i, j, k, t;
   int limit, mdepth;
   int gen = 0;
   plural INT32 classcall=0;
   plural INT8  class;
   plural INT32 T_p[34];

   /* Get action-to-perform from command line argument
    */
   int c;
   extern int atoi();
   extern int optind;
   extern char *optarg;
   while ((c=getopt(argc, argv, "S:M:d:s:G:")) != EOF)
      switch(c) {
       case 'S':
	 search = 1;
	 limit = atoi(optarg);
	 if (debug & 1) printf("search %d graphs\n", limit);
	 break;
       case 'M':
	 genmat = 1;
	 mdepth = atoi(optarg);
	 if (debug & 1) printf("mutant depth %d\n", mdepth);
	 break;
       case 'd':
	 debug = atoi(optarg);
	 if (debug & 1) printf("debug flags %d\n", debug);
	 break;
       case 's':
         seed = atoi(optarg);
	 if (debug & 1) printf("seed is %d\n", seed);
	 break;
       case 'G':
	 gfnx = *optarg;
	 break;
       case '?':
	 exit(127);
	 }


   /* initialize the graph G */
   initial_G(gfnx);

   /* Enumerate all distinct 3-sets; each 3-set has a unique 
    * assignment of vertices to v1-v3, in increasing order.
    */
   t = 0;
   for (i=0; i<(N-2); i++)
      for (j=i+1; j<(N-1); j++)
	 for (k=j+1; k<N; k++) {
	    proc[t].v1 = i;
	    proc[t].v2 = j;
	    proc[t].v3 = k;
	    proc[t].in3set = 1;
	    t++;
	    }

   if (debug & 1) printf("%d 3-clusters initialized\n", t);
   if (t > nproc) {
      fprintf(stderr, "%d processors < %d 3-clusters\n", nproc, t);
      exit(127);
      }

   for(; t<nproc; t++) proc[t].in3set = 0;

   /* initialize the 'bad count' sub-array B 
    */
   for (i=0; i<nproc; i++)
      proc[i].inB = 0;
   for (k=0,i=0; i<N-1; i++) 
	 for (j=i+1; j<N; j++) {
	    proc[i][j].B = 0;
	    proc[i][j].inB = 1;
	    }
   if (debug & 1) printf("'bad count' array initialized\n");

   if (in3set) {

      /* initialize PE and ACU cost vectors */
      for (i=0; i<34; i++) {
	 T[i] = 0;
	 T_p[i] = 0;
	 }

      /* step through pairs (u,v) to generate all 5-clusters 
       */
      if (debug & 1) printf("classifying all 5-clusters...\n");
      for (u=v3+1; u<N-1; u++)
	 for (v=u+1; v<N; v++) {
	    class = classify(u, v, G[u][v]);
	    classcall++;
	    T_p[class]++;
	    if (class == 0 || class == 33)
	       updateB(u, v, (plural) 1);
	    }

      num5sets = reduceAdd32(classcall);

      for (i=0; i<34; i++) {
	 T[i] = reduceAdd32(T_p[i]);
	 }

      /* check operation of classify */
      if (debug & 2) {
         printf("%d initial calls to classify()\n", num5sets);
         print_cv(T);
	 }
      check_cv(T);

      /* display the first few 'worst' vert's in B */
      if (debug & 4) {
        int i, j;
	VERT w1, w2;
	printf("First values in 'badness' array\n");
	for (i=0; i<16; i++) {
	   worstB(&w1, &w2, i);
	   j = proc[w1][w2].B;
	   printf("edge=(%2d,%2d), badness=%d, k=%d\n", w1, w2, j, i);
	   }
        }

#define MAXDEPTH 1000	/* genmat path depth */

      /* option to generate training data for NN cost function
       */
      if (genmat) {
	 FILE *fopen(); FILE *fp; 
	 int i, j, k, t, Tinit[34];
	 float temp;
	 VERT m[MAXDEPTH], n[MAXDEPTH], w1, w2;
	 if (debug & 1) printf("Generating training data...\n");
/*	 if (debug & 8) fwd_init();  */
	 fp = fopen("train.mat", "w");
	 fprint_cv2(fp, 0, T);
	 if (debug & 8) print_stat(0, -1, -1, costfnx(T), T);

	 /* root has 20 descendants */
	 acopy(T, Tinit, 34);
	 for (i=1; i<=20; i++) {
	    newclass((plural VERT) 0, (plural VERT) i);
	    flipedge(0,i);
	    acopy(Tnew, T, 34);
	    fprint_cv2(fp, 1, T);
	    if (debug & 8) print_stat(1, 0, i, costfnx(T), T);

	    /* root-child heads a walk 'mdepth' steps long */
	    for (j=0; j<mdepth; j++) {
	       m[j] = myrand(N);
	       while ((n[j]=myrand(N)) == m[j]);
	       if (n[j] < m[j]) { t=n[j]; n[j]=m[j]; m[j]=t; }
	       newclass((plural VERT) m[j], (plural VERT) n[j]);
	       flipedge(m[j],n[j]);
	       acopy(Tnew, T, 34);
	       fprint_cv2(fp, j+2, T);
	       if (debug & 8) print_stat(j+2, m[j], n[j], costfnx(T), T);
	       }
	    /* unroll changes in G */
	    for (j=mdepth-1; j>=0; j--)
	       flipedge(m[j],n[j]);
	    flipedge(0,i);
	    acopy(Tinit, T, 34);
	    }
	 fclose(fp);
	 exit(0);
	 }

      /* search loop
       */
      if (search) {
	 VERT w1, w2, w1a, w2a;
	 float cost, newcost;
	 int i, k;
	 gen = 0;
/*	 fwd_init();  		*/
	 cost = costfnx(T);
	 if (debug & 1) {
	    printf("starting search loop\n");
	    printf("initial cost=%g\n", cost);
	    }
	 while (cost != 0.0 && gen < limit) {

	    /* check 'worst' edge , second 'worst' edge, etc 
	     */
	    worstB(&w1, &w2, 0);
	    w1a = w1; 
	    w2a = w2;
	    newclass((plural) w1, (plural) w2);
	    newcost = costfnx(Tnew);
	    gen++; 
	    if (debug & 16) print_stat(gen, w1, w2, newcost, Tnew);

/*	    cost = cost * (float) 1.05;  */
/*	    cost = cost +  0.2 / (cost + 1.0);	*/
	    for (k=1; k < N*(N-1)/2 && newcost >= cost; k++, gen++) {
	       worstB(&w1, &w2, k);
	       newclass((plural) w1, (plural) w2);
	       newcost = costfnx(Tnew);
	       if (debug & 16) print_stat(gen, w1, w2, newcost, Tnew);
	       }

	    if (k == N*(N-1)/2) {
	       /* no improved cost from bad array, so take 'baddest' */
	       w1 = w1a; 
 	       w2 = w2a;
	       newclass((plural) w1, (plural) w2);
	       newcost = costfnx(Tnew);
	       }

	    /* keep the change we have found */
	    if (in5set) {
	       if (class == 0 || class == 33) 
		  updateB((plural) w1, (plural) w2, (plural) -1);
	       if (class2 == 0 || class2 == 33) 
		  updateB((plural) w1, (plural) w2, (plural) 1);
	       }
	    flipedge(w1, w2);
	    acopy(Tnew, T, 34);
	    cost = newcost;
	    if (debug & 8) print_stat(gen, w1, w2, newcost, Tnew);
	    }

	 /* save whatever we've got
	  */
	 write_G("lastG.mat");
	 if (debug & 1) {
	    printf("checked %d graphs\n", gen);
	    printf("final cost %g\n", cost);
	    }
	 }
      }	/* end of 'if in3set' block */
   }
