
#ifndef SHUTTLE_H
#define SHUTTLE_H

#define MAX_N_CLUSTERS 16

// ***** THIS FILE HANDLES THE READING OF THE DATA TO BE CLASSIFIED ******

// classification training data
Group RSGA[494021];
int nRSGA;


Group *SGA;
int nSGA;

// classification test data
Group TSGA[400000];
int nTSGA;

class Shuttle {

public:

	int mins[36];
	int maxs[36];
	float rangeinv[36];
	int catrange[36];
	float catrangeinv[36];

	char catmap[8][128][128]; // [axis][category][char]
	int numcatmap[8];
	int cattemp[8];
	char classmap[MAX_N_CLASSES][128];
	int numclasses;
	int classes[MAX_N_CLASSES];
	int numinclass[MAX_N_CLASSES];
	int numclustersinclass[MAX_N_CLASSES];
	// high-dimensional coordinates of the cluster centroids
	float *clusters[MAX_N_CLASSES][MAX_N_CLUSTERS];
	// for each projection, whether we consider each cluster centroid
	int clusterflags[MAX_N_PROJS][MAX_N_CLASSES][MAX_N_CLUSTERS];
	float centroidscreencoords[MAX_N_CLASSES][MAX_N_CLUSTERS][2];
	// for each projection, screen coordinates of cluster centroids
	float allcentroidscreencoords[MAX_N_PROJS][MAX_N_CLASSES][MAX_N_CLUSTERS][2];
	int selectedcentroidclass;
	int selectedcentroidnum;

	int nrows;


	float fmins[64];
	float fmaxs[64];

	// [projection][region][category]
	int stats[MAX_N_PROJS][32][32];
	int nonstats[MAX_N_PROJS][32]; // non-terminal projection,category
	int nontermclass[MAX_N_PROJS];
	Point *neighbors[MAX_N_PROJS][32][32];
	int no_neighbors[MAX_N_PROJS][32];

	int showall;

	Shuttle() {
		int i, j, p;
		for (p=0;p<MAX_N_PROJS;p++) {
			for (i=0;i<MAX_N_CLASSES;i++) {
				for (j=0;j<MAX_N_CLUSTERS;j++) {
					clusterflags[p][i][j] = 1;
				}
			}
		}
		showall = 0;
		//SGA = new Group[500000];
	};

	
	int Sfindredblue(int p, Group *gptr) {
		int r, a, anyhasred, hasred, hri;
		r = 1; // assume red
		anyhasred = 0;
		for (a=0;a<Group::ncatdims;a++) {
			// check that there's some red
			hasred = 0;
			for (hri=0;hri<catrange[a];hri++) {
				if (projParallel[p][a][hri]==1) { hasred = 1; anyhasred = 1; }
			}
			// if any axis is blue, then it is blue
			if (hasred&&(projParallel[p][a][gptr->catcoords[a] - mins[a]]==0)) r = 0;
		}
		if (anyhasred == 0) r = 0; // if there are no reds at all in all axes, then blue
		return r;
	}


	int readVeryStandard() {

		int i,j,j2,n;
		float range;
		int first = 1;
		nSGA = 0;
		FILE *fptr;
		int ncols, dummy;
		int foundclass;
		int ntotalcols;
		int *symbolics;
		float tempf;
		int tempi;
		int *maxcat;

		for (i=0;i<100000;i++) {
			randomshow[i] = (unsigned char)0;
		}

		for (j=0;j<MAX_N_CLASSES;j++) {
			numinclass[j] = 0;
		}

		fptr = fopen("../data/kddcup99/processed/header.txt","rt");
		fscanf(fptr,"%d %d %d", &nrows, &dummy, &ntotalcols);
		fclose(fptr);

		ncols = 0;
		symbolics = new int[ntotalcols];
		fptr = fopen("../data/kddcup99/processed/symbolics.txt","rt");
		for (i=0;i<ntotalcols;i++) {
			fscanf(fptr,"%d", &(symbolics[i]));
			if (symbolics[i]==0) ncols++;
		}
		fclose(fptr);

		maxcat = new int[ntotalcols-ncols];

		fptr = fopen("../data/kddcup99/processed/data.txt","rt");
		
		for (i=0;i<ntotalcols-ncols;i++) maxcat[i] = 0;

		SGA = RSGA;

		for (n=0;n<nrows;n++) {
			SGA[nSGA].size = 5.0;
			SGA[nSGA].SetDims(ncols);
			SGA[nSGA].SetCatDims(ntotalcols-ncols);

			j = 0; j2 = 0;
			for (i=0;i<ntotalcols;i++) {
				if (symbolics[i]==0) {
					fscanf(fptr,"%f", &tempf);
					SGA[nSGA].startingcoords[j] = tempf;
					j++;
				} else {
					fscanf(fptr,"%d", &tempi);
					SGA[nSGA].catcoords[j2] = tempi;
					if (tempi>maxcat[j2]) maxcat[j2] = tempi;
					j2++;
				}
			}

			fscanf(fptr,"%d", &(SGA[nSGA].category));

			// convert categories
			if (SGA[nSGA].category==11) {
				SGA[nSGA].category=0;
			} else if ((SGA[nSGA].category==5)||(SGA[nSGA].category==10)||(SGA[nSGA].category==15)||(SGA[nSGA].category==17)) {
				SGA[nSGA].category=1;
			} else if ((SGA[nSGA].category==0)||(SGA[nSGA].category==6)||(SGA[nSGA].category==9)||(SGA[nSGA].category==14)||(SGA[nSGA].category==18)||(SGA[nSGA].category==20)) {
				SGA[nSGA].category=2;
			} else if ((SGA[nSGA].category==1)||(SGA[nSGA].category==7)||(SGA[nSGA].category==12)||(SGA[nSGA].category==16)) {
				SGA[nSGA].category=3;
			} else if ((SGA[nSGA].category==2)||(SGA[nSGA].category==3)||(SGA[nSGA].category==4)||(SGA[nSGA].category==8)||(SGA[nSGA].category==13)||(SGA[nSGA].category==19)||(SGA[nSGA].category==21)||(SGA[nSGA].category==22)) {
				SGA[nSGA].category=4;
			}

			// set randomshow
			if (SGA[nSGA].category==0) {
				if (rand()%15 == 3) randomshow[nSGA/8] |= (1 << (nSGA%8));
			} else if (SGA[nSGA].category==1) {
				if (rand()%3 == 1) randomshow[nSGA/8] |= (1 << (nSGA%8));
			} else if (SGA[nSGA].category==2) {
				if (rand()%3 == 2) randomshow[nSGA/8] |= (1 << (nSGA%8));
			} else {
				randomshow[nSGA/8] |= (1 << (nSGA%8));
			}

			if (first) {
				for (i=0;i<ncols;i++) {
					fmins[i] = fmaxs[i] = SGA[nSGA].startingcoords[i];
				}
				first = 0;
			} else {
				for (i=0;i<ncols;i++) {
					if (SGA[nSGA].startingcoords[i] > fmaxs[i]) fmaxs[i] = SGA[nSGA].startingcoords[i];
					if (SGA[nSGA].startingcoords[i] < fmins[i]) fmins[i] = SGA[nSGA].startingcoords[i];
				}
			}

			foundclass = 0;
			for (j=0;j<numclasses;j++) {
				if (classes[j]==SGA[nSGA].category) {
					foundclass = 1;
					numinclass[j] = numinclass[j] + 1;
				}
			}
			if (foundclass==0) {
				classes[numclasses] = SGA[nSGA].category;
				numinclass[numclasses] = 1;
				numclasses++;
			}	




			nSGA++;

		}
		fclose(fptr);


		for (i=0;i<ncols;i++) {
			printf("dimension %d max %f min %f\n", i, fmaxs[i], fmins[i]);
			if ((fmaxs[i] - fmins[i])>0.01) {
				rangeinv[i] = 1.0/(fmaxs[i] - fmins[i]);
			} else {
				rangeinv[i] = 1.0;
			}
		}


		// normalize to the range 0.0 to 1.0

		for (i=0;i<nSGA;i++) {
			for (j=0;j<ncols;j++) {
				SGA[i].startingcoords[j] = rangeinv[j]*(SGA[i].startingcoords[j] - fmins[j]);
			}
		}

		for (i=0;i<ntotalcols-ncols;i++) {
			catrange[i] = maxcat[i] + 1;
			catrangeinv[i] = 1.0/(float)(catrange[i]);
		}
		for (i=0;i<nSGA;i++) {
			for (j=0;j<ntotalcols-ncols;j++) {
				SGA[i].parcoords[j] = catrangeinv[j]*(0.25+(0.5*(float)i/(float)nSGA)+(float)SGA[i].catcoords[j]);
			}
		}

		nRSGA = nSGA;		

		return ncols;

	};



	void readVeryStandardTest() {

		int i,j,j2,n;
		float range;
		int first = 1;
		
		FILE *fptr;
		int ncols, dummy;
		int foundclass;
		int ntotalcols;
		int *symbolics;
		float tempf;
		int tempi;
		int ntestrows;


		fptr = fopen("../data/kddcup99/processed/testheader.txt","rt");
		fscanf(fptr,"%d %d %d", &ntestrows, &dummy, &ntotalcols);
		fclose(fptr);

		ncols = 0;
		symbolics = new int[ntotalcols];
		fptr = fopen("../data/kddcup99/processed/symbolics.txt","rt");
		for (i=0;i<ntotalcols;i++) {
			fscanf(fptr,"%d", &(symbolics[i]));
			if (symbolics[i]==0) ncols++;
		}
		fclose(fptr);

		fptr = fopen("../data/kddcup99/processed/testdata.txt","rt");

		nTSGA = 0;
		
		for (n=0;n<ntestrows;n++) {
			TSGA[nTSGA].size = 5.0;
			TSGA[nTSGA].SetDims(ncols);
			TSGA[nTSGA].SetCatDims(ntotalcols-ncols);
			j = 0; j2 = 0;
			for (i=0;i<ntotalcols;i++) {
				if (symbolics[i]==0) {
					fscanf(fptr,"%f", &tempf);
					TSGA[nTSGA].startingcoords[j] = tempf;
					j++;
				} else {
					fscanf(fptr,"%d", &tempi);
					TSGA[nTSGA].catcoords[j2] = tempi;
					j2++;
				}
			}

			fscanf(fptr,"%d", &(TSGA[nTSGA].category));


			nTSGA++;

		}
		fclose(fptr);

		// normalize to the range 0.0 to 1.0

		for (i=0;i<nTSGA;i++) {
			for (j=0;j<ncols;j++) {
				TSGA[i].startingcoords[j] = rangeinv[j]*(TSGA[i].startingcoords[j] - fmins[j]);
			}
		}

		for (i=0;i<nTSGA;i++) {
			for (j=0;j<ntotalcols-ncols;j++) {
				TSGA[i].parcoords[j] = catrangeinv[j]*(0.25+(0.5*(float)i/(float)nTSGA)+(float)TSGA[i].catcoords[j]);
			}
		}

	};

	int ClassifyByClusters(Group G, int & finrule) {
		int c;
		c = COR.TestSavedOne(&G, finrule);
		return c;
	};


	
	int ClassifyByTree(Group G, int & finrule) {

		int c, cp, px, py, resolved, r, j, k;
		float x, y; 
		int cr, rulefound, crule;
			
		finrule = -1;
		c = -1; // invalidate predicted class
		for (j=0;((j<numProjs)&&(c<0));j++) { // for each projection
			if (projCoords[j]==0) { // star coordinates mode
				R.Project(&G,projAxes[j],expans[j],&x,&y,j,NULL); // project using radviz
				// compare radviz projection to projected regions
				// returned x,y value are on a -1.2 to 1.2 square
				// convert that to a 0 to 512 square
				px = (int)(((x+1.2)/2.4)*512.0);
				py = (int)(((y+1.2)/2.4)*512.0);
				// index into the regions projection to find the region
				if ((px>=0)&&(px<512)&&(py>=0)&&(py<512)) {
					r = (int)(Parray[j].mat[px][py]); // region
				} else {
					r = -1;
				}
			} else if (projCoords[j]==1) { // if parent projection is using parallel coordinates
				// find out if blue or red
				r = Sfindredblue(j,&G);
			}
			// look into rules to find which class corresponds to the region
			for (k=0;(k<Rooles.numRules)&&(c<0);k++) { // for each rule
				if ((Rooles.rarray[k].Flag&2)==0) { // if rule is an unconditional rule
					if ((Rooles.rarray[k].Projection==(unsigned char)j)&&
						(Rooles.rarray[k].Region==(unsigned char)r)) { // if this rule has the appropriate proj and reg
						if ((Rooles.rarray[k].Flag&1)!=0) { // if the rule definitively states a class
							c = Rooles.rarray[k].ProjClass; // predicted class
							finrule = k;
						} else { // predicted another projection: we need to project again
							resolved = 0;
							cp = Rooles.rarray[k].ProjClass; // predicted projection
							while (!resolved) {
								if (projCoords[cp]==0) { // star coordinates mode
									R.Project(&G,projAxes[cp],expans[cp],&x,&y,cp,NULL); // project using radviz
									// compare radviz projection to projected regions
									// returned x,y value are on a -1.2 to 1.2 square
									// convert that to a 0 to 512 square
									px = (int)(((x+1.2)/2.4)*512.0);
									py = (int)(((y+1.2)/2.4)*512.0);
									if ((px>=0)&&(px<512)&&(py>=0)&&(py<512)) {
										// index into the regions projection to find the region
										cr = (int)(Parray[cp].mat[px][py]); // region
									} else {
										cr = -1;
									}
								} else if (projCoords[cp]==1) { // if parent projection is using parallel coordinates
									// find out if blue or red
									cr = Sfindredblue(cp,&G);
								}
								rulefound = 0;
								for (crule=0;(crule<Rooles.numRules)&&(!rulefound);crule++) { // for each rule
								if ((Rooles.rarray[crule].Flag&2)!=0) { // if rule is a conditional projection rule
										if ((Rooles.rarray[crule].Projection==(unsigned char)cp)&&
											(Rooles.rarray[crule].Region==(unsigned char)cr)) { 
											// if proj and reg match
											if ((Rooles.rarray[crule].Flag&1)!=0) { // if the rule definitively states a class
												resolved = 1;
												c = Rooles.rarray[crule].ProjClass; // predicted class
											} else { // if rule results in projection
												cp = Rooles.rarray[crule].ProjClass;
											}
											// break out of for loop
											finrule = crule;
											rulefound = 1;
										}
									}
								}
								if (!rulefound) {
									resolved = 1;
									// here assign the majority class of the non-terminal projection because we cannot
									// find any further projection
									c = nontermclass[cp];
								}
							}
						}
					}
				}
			}
			if ((projConditions[j][0]<0)&&(finrule<0)) { // if unconditional condition
				c = nontermclass[j];
			}
		}

		return c;

	};


	void readVeryStandardTestAndClassify() {

		int i,j,j2,n;
		float range;
		int first = 1;
		nSGA = 0;
		FILE *fptr;
		int ncols, dummy;
		int foundclass;
		int ntotalcols;
		int *symbolics;
		float tempf;
		int tempi;
		int ntestrows;
		Group G;

		int k,r,c,cr,d,e;
		float x,y;
		int px,py;
		int resolved, rulefound, cp, crule;
		int numSuccess = 0;
		int numFailure = 0;
		int *rsucc;
		int *rfail;
		int *csucc;
		int *cfail;
		int finrule;
		int dist, mindist;
		int norule;

		int resultsMat[5][5];
		int totalrow;
		int costMat[5][5];
		int totalcost;

		FILE *rfile;

		norule = 0;
		rsucc = new int[Rooles.numRules];
		rfail = new int[Rooles.numRules];
		csucc = new int[COR.nsaved];
		cfail = new int[COR.nsaved];
		for (i=0;i<Rooles.numRules;i++) { rsucc[i] = rfail[i] = 0; }
		for (i=0;i<COR.nsaved;i++) { csucc[i] = cfail[i] = 0; }
		for (i=0;i<5;i++) {
			for (j=0;j<5;j++) {
				resultsMat[i][j] = 0;
			}
		}

		fptr = fopen("../data/kddcup99/processed/testheader.txt","rt");
		fscanf(fptr,"%d %d %d", &ntestrows, &dummy, &ntotalcols);
		fclose(fptr);

		ncols = 0;
		symbolics = new int[ntotalcols];
		fptr = fopen("../data/kddcup99/processed/symbolics.txt","rt");
		for (i=0;i<ntotalcols;i++) {
			fscanf(fptr,"%d", &(symbolics[i]));
			if (symbolics[i]==0) ncols++;
		}
		fclose(fptr);

		fptr = fopen("../data/kddcup99/processed/testdata.txt","rt");
		
		for (n=0;n<ntestrows;n++) {
			G.size = 5.0;
			G.SetDims(ncols);
			G.SetCatDims(ntotalcols-ncols);
			j = 0; j2 = 0;
			for (i=0;i<ntotalcols;i++) {
				if (symbolics[i]==0) {
					fscanf(fptr,"%f", &tempf);
					G.startingcoords[j] = tempf;
					j++;
				} else {
					fscanf(fptr,"%d", &tempi);
					G.catcoords[j2] = tempi;
					j2++;
				}
			}
			fscanf(fptr,"%d", &(G.category));


			// normalize to the range 0.0 to 1.0
			for (j=0;j<ncols;j++) {
				G.startingcoords[j] = rangeinv[j]*(G.startingcoords[j] - fmins[j]);
			}
			for (j=0;j<ntotalcols-ncols;j++) {
				G.parcoords[j] = catrangeinv[j]*(0.25+(0.5*(float)i/(float)nSGA)+(float)G.catcoords[j]);
			}

			// ***************** classify **********************************************************************

			c = ClassifyByClusters(G, finrule);
			if (c>=0) {
				if (c==G.category) {
					csucc[finrule]++;
				} else {
					cfail[finrule]++;
				}
			} else {
				c = ClassifyByTree(G, finrule);
			}

			// ***************** compare predicted class against known category ********************************

			if (c==G.category) {
				rsucc[finrule]++;
				numSuccess++;
			} else {
				if (finrule >= 0) {
					rfail[finrule]++;
				}
				numFailure++;
			}
			if ((c>=0) && (G.category>=0) && (c<5) && (G.category<6)) {
				if (G.category!=5) {
					resultsMat[G.category][c]++;
				} else {
					if (c==4) {
						resultsMat[4][4]++;
					} else if (c==3) {
						resultsMat[3][3]++;
					} else {
						resultsMat[4][c]++;
						resultsMat[3][c]++;
					}
				}
			}


		}
		fclose(fptr);


		printf("classification results: success %d failure %d\n", numSuccess, numFailure);
		printf("success rate %f\n", (float)numSuccess/(float)(numSuccess+numFailure));
		

		printf("Results matrix:\n");
		for (i=0;i<5;i++) {
			totalrow = 0;
			for (j=0;j<5;j++) {
				printf("%8d ", resultsMat[i][j]); totalrow += resultsMat[i][j];
			}
			printf("%6.3f\n", (float)resultsMat[i][i]/(float)totalrow);
		}
		for (i=0;i<5;i++) {
			totalrow = 0;
			for (j=0;j<5;j++) {
				totalrow += resultsMat[j][i];
			}
			printf("%8.3f ", (float)resultsMat[i][i]/(float)totalrow);
		}
		printf("\n");

		costMat[0][0] = 0;
		costMat[0][1] = 1;
		costMat[0][2] = 2;
		costMat[0][3] = 2;
		costMat[0][4] = 2;
		costMat[1][0] = 1;
		costMat[1][1] = 0;
		costMat[1][2] = 2;
		costMat[1][3] = 2;
		costMat[1][4] = 2;
		costMat[2][0] = 2;
		costMat[2][1] = 1;
		costMat[2][2] = 0;
		costMat[2][3] = 2;
		costMat[2][4] = 2;
		costMat[3][0] = 3;
		costMat[3][1] = 2;
		costMat[3][2] = 2;
		costMat[3][3] = 0;
		costMat[3][4] = 2;
		costMat[4][0] = 4;
		costMat[4][1] = 2;
		costMat[4][2] = 2;
		costMat[4][3] = 2;
		costMat[4][4] = 0;

		totalrow = totalcost = 0;
		for (i=0;i<5;i++) {
			for (j=0;j<5;j++) {
				totalrow += resultsMat[i][j];
				totalcost += costMat[i][j] * resultsMat[i][j];
			}
		}
		printf("average cost %8.5f\n", (float)totalcost/(float)totalrow);


		rfile = fopen("../saved/results.txt","wt");

		fprintf(rfile,"Results matrix:\n");
		for (i=0;i<5;i++) {
			totalrow = 0;
			for (j=0;j<5;j++) {
				fprintf(rfile,"%8d ", resultsMat[i][j]); totalrow += resultsMat[i][j];
			}
			fprintf(rfile,"%6.3f\n", (float)resultsMat[i][i]/(float)totalrow);
		}
		for (i=0;i<5;i++) {
			totalrow = 0;
			for (j=0;j<5;j++) {
				totalrow += resultsMat[j][i];
			}
			fprintf(rfile,"%8.3f ", (float)resultsMat[i][i]/(float)totalrow);
		}
		fprintf(rfile,"\n");
		totalrow = totalcost = 0;
		for (i=0;i<5;i++) {
			for (j=0;j<5;j++) {
				totalrow += resultsMat[i][j];
				totalcost += costMat[i][j] * resultsMat[i][j];
			}
		}
		fprintf(rfile,"average cost %8.5f\n", (float)totalcost/(float)totalrow);

		fclose(rfile);

	};


	

	void Display(int p) {
		int i;
		int ind1, ind2;
		for (i=0;i<nSGA;i++) {
			ind1 = i/8;
			ind2 = i%8;
			if ((obj_clusters[i]<(signed char)0)&&(showall||(((randomshow[ind1] >> ind2) & 1) == 1))) {
				if (SGA[i].category!=hide_category) R.Display(&(SGA[i]),1.0,p,0.1*(float)i/(float)nSGA);
			}
		}
		
	};
	
	void DisplayZoom(Axis *a, float *eptr, int p) {
		int i;
		int ind1, ind2;	
		for (i=0;i<nSGA;i++) {
			ind1 = i/8;
			ind2 = i%8;
			if ((obj_clusters[i]<(signed char)0)&&(showall||(((randomshow[ind1] >> ind2) & 1) == 1))) {
				if (SGA[i].category!=hide_category) R.DisplayZoom(&(SGA[i]),1.0,a,eptr,p,0.1*(float)i/(float)nSGA);
			}
		}
		
	};




	void BuildRules() {
		int i,j,k,r,np;
		float x,y;
		int px,py;
		int resolved, rulefound, cp;
		int maxclass, maxn;
		FILE *fptr;
		int m, n, totaln;
		Point *temp;
		int biggestnum, biggestclass;

		int negstats[MAX_N_PROJS][MAX_N_CLASSES];

		// ************************ build the non-terminal rules **********************************
		for (i=0;i<numProjs;i++) {
			if (projConditions[i][0]>=0) {
				// set to be an uncondtional rule leading to a projection
				Rooles.AddRule(0,projConditions[i][0],projConditions[i][1],0,i);
			}
		}
		// ************************ create statistics for terminal regions ************************
		for (i=0;i<numProjs;i++) {
			for (j=0;j<32;j++) {
				nonstats[i][j] = 0;
				negstats[i][j] = 0;
				for (k=0;k<32;k++) {
					stats[i][j][k] = 0;
				}
			}
		}
		// for each data point, follow non-terminal rules to final projection and final region
		for (i=0;i<nSGA;i++) { // for each training data point
			for (j=0;j<numProjs;j++) { // for each projection
				if (projConditions[j][0]<0) { // proceed only if the projection has no conditions
					// look into projections to find if any corresponds to the region
					resolved = 0;
					cp = j; // current projection
					while (!resolved) {
						if (projCoords[cp]==0) { // star coordinates mode
							R.Project(&(SGA[i]),projAxes[cp],expans[cp],&x,&y,cp,NULL); // project using radviz
							// compare radviz projection to projected regions
							// returned x,y value are on a -1.2 to 1.2 square
							// convert that to a 0 to 512 square
							px = (int)(((x+1.2)/2.4)*512.0);
							py = (int)(((y+1.2)/2.4)*512.0);
							// index into the regions projection to find the region
							if ((px>=0)&&(px<512)&&(py>=0)&&(py<512)) {
								r = (int)(Parray[cp].mat[px][py]); // region
							} else {
								r = -1;
							}
						} else if (projCoords[cp]==1) { // if parent projection is using parallel coordinates
							// find out if blue or red
							r = Sfindredblue(cp,&(SGA[i]));
						}
						rulefound = 0;
						for (np=0;((np<numProjs)&&(!rulefound));np++) { // for each projection
							if ((projConditions[np][0]==cp)&&(projConditions[np][1]==r)) {
								// if this proj and reg leads to another proj
								rulefound = 1;
								// here add statistics for non-terminal projection
								// nonstats[cp][SGA[i].category]++;
								cp = np;
							}
						}
						if (!rulefound) { // if this proj and reg are terminal
							if (r>=0) {
								stats[cp][r][SGA[i].category]++;
							} else {
								negstats[cp][SGA[i].category]++;
							}
							// add non terminal class for terminal projection as well
							// nonstats[cp][SGA[i].category]++;
							resolved = 1;
						}

					}
				}
			}
		}

		// **************** build terminal rules based on statistics ******************************
		for (i=0;i<numProjs;i++) { // for each projection
			for (j=0;j<32;j++) { // for each potential region
				totaln = 0;
				no_neighbors[i][j] = -1; // invalidate definite class for this projection-region pair
				maxclass = -1; maxn = 0;
				for (k=0;k<32;k++) { // for each potential class
					if (stats[i][j][k]>maxn) {
						maxclass = k;
						maxn = stats[i][j][k];
					}
					if (stats[i][j][k]>0) totaln += stats[i][j][k];
				}
				if (maxn>0) {
					Rooles.AddRule(0,i,j,1,maxclass);
					if ((((float)stats[i][j][maxclass]/(float)totaln)>0.8)&&(totaln>500)) {
						no_neighbors[i][j] = maxclass;
					}
				}
			}
			// for negative region
			maxclass = -1; maxn = 0;
			for (k=0;k<32;k++) { // for each potential class
				if (negstats[i][k]>maxn) {
					maxclass = k;
					maxn = negstats[i][k];
				}
			}
			if (maxn>0) Rooles.AddRule(0,i,-1,1,maxclass);
		}

		// **************** final pass to set the rule-cond flag for each rule ********************
		for (i=0;i<Rooles.numRules;i++) { // for each rule
			if (projConditions[Rooles.rarray[i].Projection][0]>=0) { // if projection is condtional
				Rooles.rarray[i].Flag |= 2; // set rule-cond flag
			}
		}

		// ***** for all rules, check if superceded by special rule	
		for (i=0;i<Rooles.numRules;i++) { // for each rule
			for (j=0;j<Rooles.numspecialrules;j++) {
				if ((Rooles.rarray[i].Projection==Rooles.specialrules[j][0]) &&
					(Rooles.rarray[i].Region==Rooles.specialrules[j][1])) {
					Rooles.rarray[i].ProjClass = Rooles.specialrules[j][2];
				}
			}
		}

		fptr = fopen("genrules.txt","wt");
		for (i=0;i<Rooles.numRules;i++) { // for each rule
			if (Rooles.rarray[i].Flag & 2) { 
				fprintf(fptr, "if condproj ");
			} else {
				fprintf(fptr, "if projection ");
			}
			fprintf(fptr, "%d region %d ", (int)Rooles.rarray[i].Projection, (int)Rooles.rarray[i].Region);
			if (Rooles.rarray[i].Flag & 1) { 
				fprintf(fptr, "then class ");
			} else {
				fprintf(fptr, "then projection ");
			}
			fprintf(fptr, "%d rulenum %d\n", (int)Rooles.rarray[i].ProjClass, i);
		}
		fclose(fptr);
		// ***************** set non terminal classes **********************************************
		for (i=0;i<numProjs;i++) {
			biggestnum = biggestclass = 0;
			// set to be an uncondtional rule leading to a projection
			for (j=0;j<Group::ndims;j++) {
				if (nonstats[i][j] > biggestnum) {
					biggestclass = j;
					biggestnum = nonstats[i][j];
				}
				}
			nontermclass[i] = biggestclass;
		}

	};



	void DisplayLegend() {
		int i,j;
		int c,r;
		char buffer[16];

		strcpy(buffer,"classes");
		glColor3f(0.0,0.0,0.0);
		for (i=0;i<strlen(buffer);i++) {
			glRasterPos2f(2.1+0.2*(float)i,0.3);
			glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10,buffer[i]);
		}
		
		for (i=0;i<numclasses;i++) {
			c = classes[i];
			for (j=0;j<numclustersinclass[i];j++) {
				
				glColor3f(0.5,0.0,0.0);
				for (r=0;r<strlen(buffer);r++) {
					glRasterPos2f(4.1+0.5*(float)c+0.2*(float)r,0.6+0.5*(float)j);
					glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10,buffer[r]);
				}

				glColor3f(RadVizColors[c][0],RadVizColors[c][1],RadVizColors[c][2]);
				glBegin(GL_QUADS);
				glVertex3f(4.1+0.5*(float)c,0.1+0.5*(float)j,0.0);
				glVertex3f(4.4+0.5*(float)c,0.1+0.5*(float)j,0.0);
				glVertex3f(4.4+0.5*(float)c,0.5+0.5*(float)j,0.0);
				glVertex3f(4.1+0.5*(float)c,0.5+0.5*(float)j,0.0);
				glEnd();
				if (clusterflags[currentProj][i][j]==1) {
					glColor3f(0.0,0.0,0.0);
					glLineWidth(1);
					glBegin(GL_LINES);
					glVertex3f(4.05+0.5*(float)c,0.05+0.5*(float)j,0.0);
					glVertex3f(4.45+0.5*(float)c,0.05+0.5*(float)j,0.0);
					glVertex3f(4.05+0.5*(float)c,0.05+0.5*(float)j,0.0);
					glVertex3f(4.05+0.5*(float)c,0.50+0.5*(float)j,0.0);
					glEnd();
				}
			}
		}

	};

};




#endif