#include "SeeGraphHelper.h"
#include "mst.h"

/* Windows */
int mainWindow;

/* Data */
int viewMode=VIEW_SPLAT;
int viewOption=OPT_SINGLEGVLSTEP;
graph* g;
#define MAX_SUBGRAPHS 16
int numSGs = 0;
graph* sg[MAX_SUBGRAPHS];

/* Render info */
viewData view;
unsigned char* voxels;
int* dimSize = NULL;
float** tlut   = NULL;
GLubyte*** colorBar = NULL;
//float** btdtlut = NULL;
float degreeRot =0.f;
float degreeStep=0.f;
float textColor[3];

/* Viewing data structs */
ArcBall Spaceball;
GLUquadricObj* quadratic;
//int resX=3,resY=3;//default
//int resX=8,resY=4;//better
int resX=12,resY=12;//for pictures

/* Viewing Options */
int  GL_WIN_WIDTH        = 800;
int  GL_WIN_HEIGHT       = 600;
bool blackBackground     = 0;
bool fullScreen          = 0;
bool showText			    = 1;
bool calcFrameRate       = 1;
bool showNormals         = 0;
bool updateNormals       = 0;
bool showDisconnected    = 0;
bool showUnselected      = 1;
bool showVertHalos       = 0;
bool showEdges           = 1;
bool showEdgesVisibleVerts = 1;
bool normalizeEdgeColors = 1;
bool separateEdgeColors  = 1;
bool showOnlyLCC         = 0;
bool showFullEdges       = 0;
bool showFullVerts       = 0;
bool activeQuery         = 1;
bool allowPicking        = 0;
bool loadBTDMatrixPGM    = 0;
bool showBTDMatrix       = 1;
bool showBTDMatrixBelt   = 1;
bool showBTDMatrixStretch= 0;
bool showBTDSelections   = 1;
bool showBTDasLOD        = 0;
bool showSubgraphs       = 1;
bool showSeeEdit         = 1;
bool showVertLabels      = 0;
bool saveAndExit         = 0;
bool useTwoDLayout       = 0;
bool useEBLayout         = 0;
bool useCustomLUT        = 1;
bool useBipartiteLayout  = 0;

/* Command line parameters */
float threshold;
int absoluteValue;

/* Dynamic threshold handling */
int numFiles = 0;
int displayThresh=0;//current-thresh index
float startThresh;
float stopThresh;
float stepThresh;

/* Hypothesis definition */
int currentQuery = 0;
int numFeatures  = 1;
float textX = 0.2f;
float textY = 0.985f;
float textH = 0.02f;

/* MST Value */
float MSTExtraEdgeValue = 1.05f;

/* Neural Net data */
NeuralNet* n = NULL;
Data *learn = NULL, *test = NULL;

/* BTD Matrix data */
BTD* d2;
graph* lod;
graph* glod; //LOD switch pointer

/* Useful things (voxel hacks) */
int numDisplayVerts;
int numDisplayEdges;
int geneNodeColor = 255;
int geneNodeColorQ1 = 220;
int geneNodeColorQ2 = 50;
int geneTransparentColor=0,geneTransparentAlpha=0;
int geneNodeSelectColor = 128;
int geneNodeUnSelectColor = 192;
clock_t oldTime, newTime;

unsigned char* pixels = NULL;
int pixCount=0;
int karCount=0;
char* pixFile = NULL;
char* karFile = NULL;

extern bool RecordActions;

//------------------------------------------------------------------------
// Prints the correct usage and exits upon incorrect parameters
//------------------------------------------------------------------------
void PrintUsage() {
	printf("Usage: Seegraph viewOption[1] filename.wel threshold[0.0-1.0] absValFlag[0,1] [r]\n");
	printf("   Ex: Seegraph 1 rma_pearsons.85.wel 0.85 1\n");
	printf("Usage: Seegraph viewOption[2] filename.wel startThresh endThresh stepThresh absValFlag [r]\n");
	printf("   Ex: Seegraph 2 rma_pearsons.85.wel 0.85 0.95 0.01 1\n");
	printf("Usage: Seegraph viewOption[3] filename.wel startThresh stepThresh absValFlag [r]\n");
	printf("   Ex: Seegraph 3 rma_pearsons.85.wel 0.85 0.01 1\n");
	printf("Usage: Seegraph filename.sgs\n");
	printf("   Ex: Seegraph rma_script.sgs\n");
	exit(-1);
}

//------------------------------------------------------------------------
// Function that specifies how voxels should be colored during splatting
//------------------------------------------------------------------------
void voxtransf(float* rgba, void* vox)
{
  unsigned char* voxel = (unsigned char*) vox;
  int id = voxel[0];
  rgba[0] = tlut[id][0];
  rgba[1] = tlut[id][1];
  rgba[2] = tlut[id][2];
  rgba[3] = tlut[id][3];
}

//------------------------------------------------------------------------
// Color table used throughout SeeGraph
//------------------------------------------------------------------------
void initColor(int index) {
  int i;
  float* tlut1d = NULL;

  // Voxel splat color table
  if(tlut!=NULL) {
	  for(i=0;i<256;i++) if(tlut[i]!=NULL) free(tlut[i]);
  }
  free(tlut);
  if(tlut1d!=NULL) free(tlut1d);
  tlut=(float**)malloc(256*sizeof(float*));
  if(tlut==NULL) {
	  perror("  initColor(): tlut"); exit(-1);
  }
  for (i=0; i<256; i++) {
	  tlut[i] = (float*)malloc(4*sizeof(float));
	  if(tlut[i]==NULL) {
		  perror("  initColor(): tlut"); exit(-1);
	  }
	  tlut[i][3] = 1.f;
  }

  tlut1d = (float *)malloc(256*4*sizeof(float));
  if(tlut==NULL) {
	  perror("  initColor(): tlut"); exit(-1);
  }
  if(index==HOTCOLD) vcbCustomizeColorTableColdHot(tlut1d, 0, 255);
  else if (index==FULLCOLOR) vcbCustomizeColorTableFullColor(tlut1d, 0, 255);
  else if (index==BLACKYELLOW) vcbCustomizeColorTable(tlut1d, 0.f,0.f,0.f, 1.f,1.f,0.f, 0,255);
  else if (index==BLACKREDYELLOW) {
	vcbCustomizeColorTable(tlut1d, 0.f,0.f,0.f, 0.5686274f,0.2f,0.160784f, 0,128);
	vcbCustomizeColorTable(tlut1d, 0.5686274f,0.2f,0.160784f, 1.f,1.f,0.0549019f,128,255);
	//vcb error? the below line puts the last color on 158=128+20 and doesn't initialize the rest
	//vcbCustomizeColorTable(&(tlut1d[128]), 0.5686274f,0.2f,0.160784f, 1.f,1.f,0.0549019f,0,127);
  }
  else { printf("Unknown index\n"); exit(1); }
  for(i=0;i<256;i++) tlut[i][0]=tlut1d[i*4], tlut[i][1]=tlut1d[i*4+1], tlut[i][2]=tlut1d[i*4+2];
  tlut[1][3]=0.05f;
  tlut[0][3]=0;
  tlut[geneTransparentColor][3]=(float)geneTransparentAlpha;
  free(tlut1d);
}

//------------------------------------------------------------------------
// Color table used throughout SeeGraph
//------------------------------------------------------------------------
int initColorFromFile(char* fn) {
   int i,r,g,b;
   FILE* fin;
   fin=fopen(fn,"r");
   if(fin==NULL) return -1;

   // Voxel splat color table
   if(tlut!=NULL) {
     for(i=0;i<256;i++) if(tlut[i]!=NULL) free(tlut[i]);
   }
   free(tlut);
   tlut=(float**)malloc(256*sizeof(float*));
   if(tlut==NULL) {
     perror("  initColor(): tlut"); return -1;
   }
   for (i=0; i<256; i++) {
     tlut[i] = (float*)malloc(4*sizeof(float));
     if(tlut[i]==NULL) {
	     perror("  initColor(): tlut"); return -1;
     }
     tlut[i][3] = 1.f;
   }

   for(i=0;i<256;i++) {
      fscanf(fin,"%d %d %d\n",&r,&g,&b);
      tlut[i][0]=r/255.f;
      tlut[i][1]=g/255.f;
      tlut[i][2]=b/255.f;
   }

   tlut[1][3]=0.05f;
   tlut[0][3]=0;
   tlut[geneTransparentColor][3]=(float)geneTransparentAlpha;
   return 0;
}

//------------------------------------------------------------------------
/// Allocate and initialize all queries
//------------------------------------------------------------------------
int initQueriesSmart(int numQueries,int numFeatures,graph* g) {
   char* filenameData;

   filenameData=(char*)malloc( (strlen(g->name)+1)*sizeof(char) );
   if(filenameData==NULL) {
      perror("initQueriesSmart(): filenameData");
   }
   strcpy(filenameData,g->name);
   strcpy(&filenameData[strlen(g->name)-4],".txt");
   /*if(!fexists(filenameData)) {
      printf("Not cleared bc %s doesn't exist.\n",filenameData);
      free(filenameData);
      return -1;
   }*/
initQueries(g->numQueries,numFeatures,g);
   if(showSeeEdit && strstr(filenameData,"demo3")!=NULL) initHypothesesSmartSE_D3(g);
   else if(showSeeEdit) initHypothesesSmartSE(g);
   else initHypothesesSmart(g);

   free(filenameData);
   return 0;
}

//------------------------------------------------------------------------
/// The code here initializes data from our dataset(s)
//------------------------------------------------------------------------
void initData(int argc, char** argv) {
	int i, errFlag, numVertLocs, layoutArgc;
	float numSteps;
	char *filename, *extension, *filenameGVL, *filenameBTD, *filenameData;
	char** layoutArgv;

	int x = 2;

	//Check for Recording script
	if(!strcmp(argv[argc-1],"r"))
	{
		InitRecordScript();
		argc--;
	}
	//Read the input file
	if(argc<2 || argc>7) {
		PrintUsage();
	}
	else if (argc==2 && strcmp(argv[1]+(strlen(argv[1])-4),".sgs")==0)
	{
		InitSeeGScript(argv[1]); //Script should load data with proper parameters
		return;
	}
	else if (argc==2 || argc==3) {
		filename=argv[1];
		if(strcmp(filename+strlen(filename)-4,".MST")==0)
		{
			createMST(filename,MSTExtraEdgeValue);
			strcpy(filename+strlen(filename)-3,"wel");
		}
		x = 1;
		if(argc==3) {
			threshold = (float)atof(argv[2]);
		}
		else {
			printf("Reading the edge list for threshold...\n");
			if( getLowestWeight(filename,&threshold)==-1 ) {
				printf("  Read %s\n",filename);
				PrintUsage();
			}
			printf("  Read %s\n",filename);
		}

		viewMode     = VIEW_QUADRICS;
		viewOption   = OPT_SINGLEGVLSTEP;
		absoluteValue= 1;
		startThresh  = (float)(threshold-(1e-6));
		stopThresh   = 1.f;
		stepThresh   = stopThresh-startThresh;

		//How many steps of threshold do we want?
		numSteps = 1.f;
		if( ((int)(stepThresh*100))%15==0 ) numSteps=15.f;
		else if( ((int)(stepThresh*100))%10==0 ) numSteps=10.f;
		stepThresh/=numSteps;
	}
	
	else if(atoi(argv[1])==1) {
		if(argc==3) {
			filename=argv[2];
			threshold     = -3.0;
			absoluteValue = 1;
		} 
		else if(argc==4) {
			filename  = argv[2];
			threshold = (float)atof(argv[3]);
			absoluteValue = 1;
		}
		else {
			filename      = argv[2];
			threshold     = (float)atof(argv[3]);
			absoluteValue = atoi(argv[4]);
		}
	}
	else if(atoi(argv[1])==2) {
		if(argc==7) {
			viewMode   = VIEW_SPLAT;
			viewOption = OPT_MULTIGVL;
			filename      = argv[2];
			threshold     = (float)atof(argv[3]);
			startThresh   = (float)atof(argv[3]);
			stopThresh    = (float)atof(argv[4]);
			stepThresh    = (float)atof(argv[5]);
			absoluteValue = atoi(argv[6]);
		}
		else PrintUsage();
	}
	else if(atoi(argv[1])==3) {
		if(argc==6) {
			viewMode   = VIEW_QUADRICS;
			viewOption = OPT_SINGLEGVLSTEP;
			filename      = argv[2];
			threshold     = (float)atof(argv[3]);
			startThresh   = (float)atof(argv[3]);
			stopThresh = 0.9999f;
			stepThresh    = (float)atof(argv[4]);
			absoluteValue = atoi(argv[5]);
		}
		else if(argc==7) {
			viewMode   = VIEW_QUADRICS;
			viewOption = OPT_SINGLEGVLSTEP;
			filename      = argv[2];
			threshold     = (float)atof(argv[3]);
			startThresh   = (float)atof(argv[3]);
			stepThresh    = (float)atof(argv[4]);
			stopThresh    = (float)atof(argv[5]);
			absoluteValue = atoi(argv[6]);
		}
		else PrintUsage();
	}
	else {
		printf("Unknown view option selected\n");
		PrintUsage();
	}

	
	if(strcmp(filename+strlen(filename)-4,".MST")==0)
	{
		createMST(filename,MSTExtraEdgeValue);
		strcpy(filename+strlen(filename)-3,"wel");
	}

	if(RecordActions)
		RecordScript_Load(argv,argc,x);

	//Main graph data structure
	for(i=0;i<MAX_SUBGRAPHS;i++) sg[i]=NULL;
	g=createGraph();
	if (g->name!=NULL) free(g->name);
	g->name=(char*)malloc( (strlen(filename)+1)*sizeof(char) );
	if(g->name==NULL) {
		perror("initData() g->name");
	}
	strcpy(g->name,filename);

	//Assign colors for color-coding query results
	g->numQueries = 3;
	if((g->queryColors=(int*)malloc(g->numQueries*sizeof(int))) == NULL) {
		perror("initData() queryColors"); exit(-1);
	}
	if(g->numQueries>0) g->queryColors[0]=geneNodeColor;
	if(g->numQueries>1) g->queryColors[1]=geneNodeColorQ1;
	if(g->numQueries>2) g->queryColors[2]=geneNodeColorQ2;
	if(g->numQueries>3) {
		printf(" WARNING: InitData() needs new query color-codes\n");
		for(i=3;i<g->numQueries;i++) g->queryColors[i]=0;
	}

	//Read the input file
	printf("Reading the edge list...\n");
	
	//Assumes 3-letter extension convention
	extension=filename+(strlen(filename)-4);
	if(strcmp(extension,".ltm")==0)
		g->numVerts = readLTMDefaultWeighted2D(filename, threshold,absoluteValue, g);
	else
		g->numVerts = readDefaultWeighted2D(filename, threshold, absoluteValue, g);

	if(g->numVerts<0) exit(-1);
	else printf("  Read %s\n",filename);

	//printMinMaxWeight(g,absoluteValue);
	numDisplayEdges=initNumEdges(g);

	//Write data to another format
	//write_graph("DefaultOutMatrix.txt", g, "row_major");

	//Show all edges and vertices
	if(g->vertColorList!=NULL) free(g->vertColorList);
	if((g->vertColorList=(int*)malloc(g->numVerts*sizeof(int))) == NULL) {
		perror("showAll() vertColorList"); exit(-1);
	}
	showAll(g,geneNodeColor);

	VCB_DEBUG_CODE {
		printf("  Num Verts = %d\n",g->numVerts);
		printf("  Weight 1 0 = %d\n",g->mat2d[1][0]);
		printf("  Weight 2 0 = %d\n",g->mat2d[2][0]);
		printf("  Weight 2 1 = %d\n",g->mat2d[2][1]);
		printf("  Weight %d %d = %d\n",g->numVerts-1,g->numVerts-3,g->mat2d[g->numVerts-1][g->numVerts-3]);
		printf("  Weight %d %d = %d\n",g->numVerts-1,g->numVerts-2,g->mat2d[g->numVerts-1][g->numVerts-2]);
	}

	//Initialize view parameters
	view.eye[0]=0.f;
	view.eye[1]=0.f;
	view.eye[2]=230.f;//32800.0f}; (for SHRT_MIN/MAX)
	view.zoom=0.f;
	// Light Position
	view.lightPos[0]=view.eye[0];
	view.lightPos[1]=view.eye[1];
	view.lightPos[2]=view.eye[2];
	view.lightPos[3]=1.f;
	// Ambient Light Values
	view.lightSpec[0]=1.0f;
	view.lightSpec[1]=1.0f;
	view.lightSpec[2]=1.0f;
	view.lightSpec[3]=1.0f;
	// Ambient Light Values
	view.lightAmb[0]=0.01f;
	view.lightAmb[1]=0.01f;
	view.lightAmb[2]=0.01f;
	view.lightAmb[3]=1.0f;
	// Diffuse Light Values
	view.lightDif[0]=0.8f;
	view.lightDif[1]=0.8f;
	view.lightDif[2]=0.8f;
	view.lightDif[3]=1.0f;

   if(blackBackground) {
      textColor[0] = 1.f;
      textColor[1] = 1.f;
      textColor[2] = 1.f;
   } else {
      textColor[0] = 0.f;
      textColor[1] = 0.f;
      textColor[2] = 0.f;
   }

	//Generate 3D layout if it doesn't exist yet
	printf("Getting vertex locations...\n");
	filenameGVL=(char*)malloc( (strlen(filename)+4)*sizeof(char) );
	if(filenameGVL==NULL) {perror("initData(): filenameGVL");}
	strcpy(filenameGVL,filename);
	if(useTwoDLayout) strcpy(&filenameGVL[strlen(filename)-4],"_2d.gvl");
	else strcpy(&filenameGVL[strlen(filename)-4],".gvl");

	if( !fexists(filenameGVL) ) {
		if(useBipartiteLayout == 1) {
			if(createBipartiteLayout(filenameGVL,g)==1)
				printf("  This graph is bipartite, created special layout\n");
			else
			{
				layoutConverser(&layoutArgc,&layoutArgv,filename,threshold,absoluteValue);  
				createLayout(layoutArgc,layoutArgv,useTwoDLayout,useEBLayout);
			}
		}
		else {
			layoutConverser(&layoutArgc,&layoutArgv,filename,threshold,absoluteValue); 
			createLayout(layoutArgc,layoutArgv,useTwoDLayout,useEBLayout);
		}
	}

	if (viewOption==OPT_MULTIGVL) {
		for(float t=startThresh;t<=(stopThresh+(1e-6));t+=stepThresh) numFiles++;
		printf("NumFiles = %d\n",numFiles);
		view.locList=(short int**)malloc(numFiles*sizeof(short int*));
		float currentThresh=startThresh;
		for(i=0;i<numFiles;i++) {
			numVertLocs=loadGVL(filenameGVL,&(view.locList[i]));

			if(numVertLocs < 0) {
				fprintf(stderr,"%s does not exist"); return;
			}
			printf("  Read %s\n",filenameGVL);
			free(filenameGVL);
			if( (currentThresh==startThresh) && (g->numVerts!=numVertLocs) ) {
				printf("GVL does not correspond with input file\n");
				exit(1);
			}
			Normalize(g->numVerts,0,DIM_SIZE-1,&(view.locList[i]));
			VCB_DEBUG_CODE {
				printf("  After normalization\n");
				printf("  Vert 1 = (%d,%d,%d)\n",view.locList[i][0],view.locList[i][1],view.locList[i][2]);
				printf("  Vert 2 = (%d,%d,%d)\n",view.locList[i][3],view.locList[i][4],view.locList[i][5]);
				printf("  Vert %d = (%d,%d,%d)\n",g->numVerts-1,view.locList[i][3*g->numVerts-3],view.locList[i][3*g->numVerts-2],view.locList[i][3*g->numVerts-1]);
			}
			currentThresh+=stepThresh;
		}
	}
	else if (viewOption==OPT_SINGLEGVLSTEP) {
		for(float t=startThresh;t<=(stopThresh+(1e-6));t+=stepThresh) numFiles++;
		printf("  Num Steps Available = %d\n",numFiles);
		errFlag = initOneLocation(filenameGVL, g);
		if(errFlag != 0) exit(-1);
	}
	else {
		errFlag = initOneLocation(filenameGVL, g);
		if(errFlag != 0) exit(-1);
	}
	free(filenameGVL);
	//Initialize vertex locations randomly
	//initRandom(3,0,DIM_SIZE-1,g->numVerts,&(view.locList[0]));

	if(useCustomLUT) errFlag=initColorFromFile((char*)"./default.lut");
	else initColor(HOTCOLD);
	if(errFlag==-1) {
		printf("  Custom color table not found, using default\n");
		useCustomLUT=0;
		initColor(HOTCOLD);
	}
	//colorBar = initColorBar(tlut);

	//Load voxel data
	printf("Initializing voxels...\n");
	dimSize=(int*)malloc(3*sizeof(int));
	for(i=0;i<3;i++) dimSize[i]=DIM_SIZE;
	initVoxelsWithRandNormals(g->numVerts,&voxels,geneNodeColor,view);
	for(i=0;i<g->numVerts;i++) g->vertColorList[i]=(int)voxels[i*10];
	//initVoxelsWithZNormals(g->numVerts,&voxels,geneNodeColor,view);

	//Init graph data structure
	printf("Counting vertice edges...\n");
	g->vertPosn = view.locList[0];
	g->degree=NULL;
	//initDegreeBase(g);
	initDegreeAll(g,startThresh,stopThresh,stepThresh,absoluteValue);
	numDisplayVerts=showAllVerts(g);
	updateVoxelPosition(displayThresh, g);

	//Initialize our seeEdit query viewport
	if(showSeeEdit) {
		printf("Installing seeEdit viewport...\n");
		//if(seinit(mainWindow, GL_WIN_WIDTH/4,0, 3*GL_WIN_WIDTH/4,100)==-1) exit(-1);
		if(seinit(mainWindow, 10,10, GL_WIN_WIDTH-2*10,70)==INVALID) exit(-1);
	}
	secolor(SE_BACKGROUND,0.93f,0.93f,0.93f);
	secolor(SE_INFO,0.8f,0.2f,0.2f);
	secolor(SE_TEXT,0.2f,0.5f,0.2f);
	secolor(SE_BBOX,0.f,0.f,0.f);
	registerCallback(queryChanged);

	//Build b-tree for searching
	printf("Building b-tree...\n");
	filenameData=(char*)malloc( (strlen(filename)+1)*sizeof(char) );
	if(filenameData==NULL) {
		 perror("initData(): filenameData");
	}
	strcpy(filenameData,filename);
	strcpy(&filenameData[strlen(filename)-4],".txt");
	if(fexists(filenameData)) {
		loadVertData(filenameData,g);
	}
	else BuildTree(numFeatures, g);

	initQueriesSmart(g->numQueries,numFeatures,g);
	BuildLookup(g);

	//Load block-tridiagonalized data for this file
	if(showBTDMatrix) {
		printf("Getting block tridiagonalization matrix...\n");
		filenameBTD=(char*)malloc( (strlen(filename)+1)*sizeof(char) );
		strcpy(filenameBTD,filename);


		if(!useCustomLUT) {
			initColor(BLACKREDYELLOW); //Change LUT for BTD 
			//initColor(BLACKYELLOW); //Change LUT for BTD 
			//initColor(FULLCOLOR); //Change LUT for BTD 
			//Just for testing custom color tables
			//vcbCustomizeColorTable(tlut1d, 0.f,0.f,0.f, 0.5686274f,0.2f,0.160784f, 0,128);
			//vcbCustomizeColorTable(tlut1d, 0.5686274f,0.2f,0.160784f, 1.f,1.f,0.0549019f,128,255);
		}
		
		if(loadBTDMatrixPGM) {
			strcpy(&filenameBTD[strlen(filename)-4],".pgm");
			d2=getBTDMatrixPGM(filenameBTD,g,threshold,tlut,GL_TEXTURE2);
		} else {
			strcpy(&filenameBTD[strlen(filename)-4],".btv");
			d2=getBTDMatrixPV(filenameBTD,g,threshold,tlut,GL_TEXTURE2);
		}
		if(!useCustomLUT) initColor(HOTCOLD);  //Then change back
		if(d2==NULL) {
			loadBTDMatrixPGM    = 0;
			showBTDMatrix       = 0;
			showBTDMatrixBelt   = 0;
			showBTDMatrixStretch= 0;
			showBTDSelections   = 0;
		} 
		else d2->drawingBelt = showBTDMatrixBelt;
		free(filenameBTD);
	}

	//--Neural Network stuff
	//Allocate and store all datasets (train,validate,test)
	if(! (n=(NeuralNet*)malloc(sizeof(NeuralNet))) ) { 
		perror("neural net"); exit(-1);
	}
	createNN(n,1,INPUT_NODES,HIDDEN_NODES,OUTPUT_NODES,LEARNING_RATE);
	assignRandomWeights(n);
	if(! (learn=(Data*)malloc(sizeof(Data))) ) {
		perror("learn data"); exit(-1);
	}
	if(! (test=(Data*)malloc(sizeof(Data))) ) {
		perror("test data"); exit(-1);
	}
	learn->item=NULL; learn->in=learn->sup=NULL;
	initLearnData(learn, g, 0);
	test->item=NULL; test->in=test->sup=NULL;
	initLearnData(test, g, 1);

   //For creating graphs in Excel of a diagonal of a BTD belt
   //int test[3] = {6148, 6000, 1268};
   //graphBTD(d2,g,"output-graph.txt",3,test);
   //printf("  wrote %s\n","output-graph.txt");

	//Spring Model Calculation
	//printf("Calculating spring model...\n");
	//g = init_graph(filename, g->numVerts, g->numVerts*g->numVerts/2, mat2d, label);

	//springModel3D(g->numVerts,(const short**)mat2d,&location,DIM_SIZE,1);

	//printf("Calculating edge-length heuristic...\n");
	//edgeLength(g->numVerts,mat,&location);

	//printf("Printing the edge list...\n");
	//printDefaultWeightedEdgeList("rma_pearsons.85.out.wel", g->numVerts, (void*)mat, label->list, -3.0);
}

//------------------------------------------------------------------------
/// The code here is mainly designed to set OpenGl's initial configuration.
//------------------------------------------------------------------------
void initOGL()
{
	glShadeModel(GL_SMOOTH);							// Enable Smooth Shading
	glEnable(GL_DEPTH_TEST);							// Enables Depth  Testing
	glDepthFunc(GL_LEQUAL);								// The Type Of Depth Testing To Do
	glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_FASTEST);	// Fast Perspective Calculations
	if(blackBackground) glClearColor(0.0f, 0.0f, 0.0f, 0.5f); // Black Background
	else glClearColor(1.0f, 1.0f, 1.0f, 0.5f);	// White Background
	glClearDepth(1.0f);									// Depth Buffer Setup

	Spaceball.SetWidth(GL_WIN_WIDTH);
	Spaceball.SetHeight(GL_WIN_HEIGHT);

	quadratic=gluNewQuadric();	// Create A Pointer To The Quadric Object
	gluQuadricOrientation(quadratic, GLU_OUTSIDE);
	gluQuadricNormals(quadratic, GLU_SMOOTH);	// Create Smooth Normals (or GLU_FLAT/NONE)
	gluQuadricTexture(quadratic, GL_TRUE);		// Create Texture Coords

	glLightfv(GL_LIGHT1, GL_SPECULAR, view.lightSpec);	// Setup The Ambient Light
	glLightfv(GL_LIGHT1, GL_AMBIENT, view.lightAmb);	// Setup The Ambient Light
	glLightfv(GL_LIGHT1, GL_DIFFUSE, view.lightDif);	// Setup The Diffuse Light
	glLightfv(GL_LIGHT1, GL_POSITION,view.lightPos);	// Position The Light
	glEnable(GL_LIGHT1);								// Enable Light One
	glEnable(GL_LIGHTING);

	// Blending is currently only used for splatting of vertex halos
	glEnable(GL_ALPHA_TEST);
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);

	glEnable(GL_COLOR_MATERIAL);
	if(fullScreen) glutFullScreen();

	RunVertQuery(g);
	g->edgeList1d=NULL;
	//RunEdgeQuery(getCurrentThresh(),absoluteValue,g);
	updateDisplayVars(g);

	return;
}

//------------------------------------------------------------------------
// Prints the current query data
//------------------------------------------------------------------------
void printQuery(int queryID, graph* g) {
  int i;

  if(queryID>=g->numQueries || g->fltQuery==NULL) return;
  for(i=0;i<g->numVertFeatures;i++) {
    printf("[%.1f,%.1f] ",g->fltQuery[queryID][i*2],g->fltQuery[queryID][i*2+1]);
  }
  printf("\n");
}

//------------------------------------------------------------------------
// Initialize user-interactive query data
//------------------------------------------------------------------------
void initQueries(int numQueries, int numFeatures, graph* g) {
	int i;
	g->numQueries=numQueries;

  if(g->fltQuery!=NULL) {
    for(i=0;i<g->numQueries;i++) {
       if(g->fltQuery[i]!=NULL) {
          free(g->fltQuery[i]);
          g->fltQuery[i]=NULL;
       }
    }
    free(g->fltQuery);
  }
	g->fltQuery=(float**)malloc(g->numQueries*sizeof(float*));
  if(g->fltQuery==NULL) {
    perror("initQueries() fltQuery");
  }

	for(i=0;i<g->numQueries;i++) g->fltQuery[i]=NULL;
	if( !showSeeEdit) {
    if(g->strQuery!=NULL) {
      free(g->strQuery);
    }
		g->strQuery=(char**)malloc(g->numQueries*sizeof(char*));
    if(g->strQuery==NULL) {
      perror("initQueries() g->strQuery");
    }
		for(i=0;i<g->numQueries;i++) g->strQuery[i]=NULL;
	}
}

//------------------------------------------------------------------------
// Hypothesis to be used for querying
//------------------------------------------------------------------------
void initHypotheses(graph* g) {
	int i, j;
	if(g->numVertFeatures>0) {
      if(g->fltQuery==NULL) {
         perror("initHypotheses() fltQuery");
      }
		for(i=0;i<g->numQueries;i++) {
         if(g->fltQuery[i]!=NULL) free(g->fltQuery[i]);
			g->fltQuery[i]=(float*)malloc(2*g->numVertFeatures*sizeof(float));
         if(g->fltQuery[i]==NULL) {
            perror("initHypotheses() fltQuery[]");
         }
			for(j=0;j<g->numVertFeatures;j++) {
				g->fltQuery[i][2*j]=0.f;
				g->fltQuery[i][2*j+1]=255.f;
			}
			getDisplayEquation(g->numVertFeatures,g->fltQuery[i],&(g->strQuery[i]));
		}
	}
	else VCB_WARNING("No hypothesis to initialize: build tree data first\n");
}

//------------------------------------------------------------------------
// Hypothesis to be used for querying
//------------------------------------------------------------------------
void initHypothesesSmart(graph* g) {
	int i, j;
	float min,max;
	if(g->numVertFeatures>0) {
      if(g->fltQuery==NULL) {
         perror("initHypothesesSmart() fltQuery");
      }
		for(i=0;i<g->numQueries;i++) {
         if(g->fltQuery[i]!=NULL) free(g->fltQuery[i]);
			g->fltQuery[i]=(float*)malloc(2*g->numVertFeatures*sizeof(float));
         if(g->fltQuery[i]==NULL) {
            perror("initHypothesesSmart() fltQuery[]");
         }
			for(j=0;j<g->numVertFeatures;j++) {
				getMinMax(g,j,&min,&max);
				g->fltQuery[i][2*j]=min;
				g->fltQuery[i][2*j+1]=max;
			}
			getDisplayEquation(g->numVertFeatures,g->fltQuery[i],&(g->strQuery[i]));
		}
	}
	else VCB_WARNING("No hypothesis to initialize: build tree data first\n");
	//printf("%s\n",g->strQuery[0]);
}

//------------------------------------------------------------------------
// Hypothesis to be used for querying
//------------------------------------------------------------------------
void initHypothesesSmartSE(graph* g) {
	int i, j, tlutid, queryID;
	float min,max;
	char queryString[1024];

	//Add Queries to viewport
	seformatdefault("%.1f");
	if(g->numVertFeatures>0) {
      if(g->fltQuery==NULL) {
         perror("initHypothesesSmartSE() fltQuery");
      }
		for(i=0;i<g->numQueries;i++) {

         if(g->fltQuery[i]!=NULL) free(g->fltQuery[i]);
			g->fltQuery[i]=(float*)malloc(2*g->numVertFeatures*sizeof(float));
         if(g->fltQuery[i]==NULL) {
            perror("initHypothesesSmartSE() fltQuery[]");
         }
			for(j=0;j<g->numVertFeatures;j++) {
				getMinMax(g,j,&min,&max);
				g->fltQuery[i][2*j]=min;
				g->fltQuery[i][2*j+1]=max;
			}
			
			sprintf(queryString,"VQ%d =",i+1);
			for(j=0;j<g->numVertFeatures;j++) strcat(queryString," [%se,%se]");
			strcat(queryString,"\n");
			//queryID=seadd("Query = [%se,%se]\n",2*g->numVertFeatures,g->fltQuery[i]);
			queryID=seadd(queryString,2*g->numVertFeatures,g->fltQuery[i]);

			//Change seeEdit format
			tlutid=g->queryColors[i];
			secolor(queryID,tlut[tlutid][0],tlut[tlutid][1],tlut[tlutid][2]);
			for(j=0;j<g->numVertFeatures;j++) {
				min=g->fltQuery[i][2*j];
				max=g->fltQuery[i][2*j+1];
				seminmax(queryID,2*j,  min,max);
				seminmax(queryID,2*j+1,min,max);
				sestep(queryID,2*j,  1.f+(1e-6f)); //for float precision
				sestep(queryID,2*j+1,1.f+(1e-6f)); //for float precision
				seinfo(queryID,2*j,"This changes LOWERBOUND of DEGREE");
				seinfo(queryID,2*j+1,"This changes UPPERBOUND of DEGREE");
			}
			//getDisplayEquation(g->numVertFeatures,g->fltQuery[i],&(g->strQuery[i]));
		}
	}
	else VCB_WARNING("No hypothesis to initialize: build tree data first\n");

	// Customized layout - will change with each new query
   if( strcmp(g->name,"demo5.wel")==0 ) {
      seloc(0, 0, 51, 18);
      seloc(0, 1, 71, 18);    seloc(0, 2, 101, 18);   seloc(0, 3, 120, 18);
      seloc(0, 4, 147, 18);   seloc(0, 5, 166, 18);   seloc(0, 6, 192, 18);
      seloc(0, 7, 212, 18);   seloc(0, 8, 238, 18);   seloc(0, 9, 259, 18);
      seloc(0, 10, 284, 18);  seloc(0, 11, 305, 18);  seloc(0, 12, 329, 18);
      seloc(0, 13, 350, 18);  seloc(0, 14, 378, 18);  seloc(0, 15, 397, 18);
      seloc(0, 16, 424, 18);  seloc(0, 17, 442, 18);  seloc(0, 18, 468, 18);
      seloc(0, 19, 489, 18);  seloc(0, 20, 516, 18);  seloc(0, 21, 534, 18);
      seloc(0, 22, 561, 18);  seloc(0, 23, 580, 18);  seloc(0, 24, 606, 18);
      seloc(0, 25, 625, 18);

      seloc(1, 0, 51, 33);
      seloc(1, 1, 71, 33);    seloc(1, 2, 101, 33);   seloc(1, 3, 120, 33);
      seloc(1, 4, 147, 33);   seloc(1, 5, 166, 33);   seloc(1, 6, 192, 33);
      seloc(1, 7, 212, 33);   seloc(1, 8, 238, 33);   seloc(1, 9, 259, 33);
      seloc(1, 10, 284, 33);  seloc(1, 11, 305, 33);  seloc(1, 12, 329, 33);
      seloc(1, 13, 350, 33);  seloc(1, 14, 378, 33);  seloc(1, 15, 397, 33);
      seloc(1, 16, 424, 33);  seloc(1, 17, 442, 33);  seloc(1, 18, 468, 33);
      seloc(1, 19, 489, 33);  seloc(1, 20, 516, 33);  seloc(1, 21, 534, 33);
      seloc(1, 22, 561, 33);  seloc(1, 23, 580, 33);  seloc(1, 24, 606, 33);
      seloc(1, 25, 625, 33);

      seloc(2, 0, 51, 48);
      seloc(2, 1, 71, 48);    seloc(2, 2, 101, 48);   seloc(2, 3, 120, 48);
      seloc(2, 4, 147, 48);   seloc(2, 5, 166, 48);   seloc(2, 6, 192, 48);
      seloc(2, 7, 212, 48);   seloc(2, 8, 238, 48);   seloc(2, 9, 259, 48);
      seloc(2, 10, 284, 48);  seloc(2, 11, 305, 48);  seloc(2, 12, 329, 48);
      seloc(2, 13, 350, 48);  seloc(2, 14, 378, 48);  seloc(2, 15, 397, 48);
      seloc(2, 16, 424, 48);  seloc(2, 17, 442, 48);  seloc(2, 18, 468, 48);
      seloc(2, 19, 489, 48);  seloc(2, 20, 516, 48);  seloc(2, 21, 534, 48);
      seloc(2, 22, 561, 48);  seloc(2, 23, 580, 48);  seloc(2, 24, 606, 48);
      seloc(2, 25, 625, 48);
      for(queryID=0;queryID<g->numQueries;queryID++) {
			for(j=0;j<g->numVertFeatures;j++) {
            if(j==0) {
			      seinfo(queryID,2*j,"This changes LOWERBOUND of COMPONENT ID");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of COMPONENT ID");
            } else if (j==1) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of BLUE GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of BLUE GROUP");
            }
            else if (j==2) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of GREEN GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of GREEN GROUP");
            }
            else if (j==3) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of PURPLE GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of PURPLE GROUP");
            }
            else if (j==4) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of RED GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of RED GROUP");
            }
            else if (j==5) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of FIG4a GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of FIG4a GROUP");
            }
            else if (j==6) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of FIG4b GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of FIG4b GROUPD");
            }
            else if (j==7) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of FIG4c GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of FIG4c GROUP");
            }
            else if (j==8) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of FIG4d GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of FIG4d GROUP");
            }
            else if (j==9) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of FIG5a GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of FIG5a GROUP");
            }
            else if (j==10) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of FIG5b GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of FIG5b GROUP");
            }
            else if (j==11) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of FIG5c GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of FIG5c GROUP");
            }
            else if (j==12) {
               seinfo(queryID,2*j,"This changes LOWERBOUND of FIG5d GROUP");
			      seinfo(queryID,2*j+1,"This changes UPPERBOUND of FIG5d GROUP");
            }
         }
      }
   }
   else if( strcmp(g->name,"demo7.wel")==0 ) {
      //Voxel correlation values so no features
		/*for(i=0;i<g->numQueries;i++) {
			for(j=0;j<g->numVertFeatures;j++) {
            if(j==0) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of Vx");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of Vx");
            }
            if(j==1) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of Vy");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of Vy");
            }
            if(j==2) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of Vz");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of Vz");
            }
            if(j==3) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of ENERGY");
				   seinfo(queryID,2*j+1,"This changes LOWERBOUND of ENERGY");
            }
			}
      }*/
   }
   else if( strcmp(g->name,"demo7.wel")==0 ) {
      seloc(0, 0, 50, 18);	   seloc(0, 1, 91, 18);
      seloc(0, 2, 127, 18);	seloc(0, 3, 171, 18);
      seloc(0, 4, 208, 18);	seloc(0, 5, 250, 18);
      seloc(0, 6, 284, 18);	seloc(0, 7, 327, 18);
	   seformat(0,6,"%.2f");
      seformat(0,7,"%.2f");

      seloc(1, 0, 51, 33);    seloc(1, 1, 91, 33);
      seloc(1, 2, 127, 33);	seloc(1, 3, 171, 33);
      seloc(1, 4, 208, 33);	seloc(1, 5, 250, 33);
      seloc(1, 6, 284, 33);	seloc(1, 7, 327, 33);
	   seformat(1,6,"%.2f");
      seformat(1,7,"%.2f");
  
      seloc(2, 0, 51, 48);    seloc(2, 1, 91, 48);
      seloc(2, 2, 127, 48);	seloc(2, 3, 171, 48);
      seloc(2, 4, 208, 48);	seloc(2, 5, 250, 48);
      seloc(2, 6, 284, 48);	seloc(2, 7, 327, 48);
	   seformat(2,6,"%.2f");
      seformat(2,7,"%.2f");
		for(i=0;i<g->numQueries;i++) {
			for(j=0;j<g->numVertFeatures;j++) {
            if(j==0) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of X");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of X");
            }
            if(j==1) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of Y");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of Y");
            }
            if(j==2) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of Z");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of Z");
            }
            if(j==3) {
				   seinfo(queryID,2*j,"This changes ATOM TYPE [0=F,1=H,2=O]");
				   seinfo(queryID,2*j+1,"This changes ATOM TYPE [0=F,1=H,2=O]");
            }
			}
      }
   }
   else if( strcmp(g->name,"demo8.wel")==0 ) {
		for(i=0;i<g->numQueries;i++) {
			for(j=0;j<g->numVertFeatures;j++) {
            if(j==0) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of X");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of X");
            }
            if(j==1) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of Y");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of Y");
            }
            if(j==2) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of X");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of Z");
            }
            if(j==3) {
				   seinfo(queryID,2*j,"This changes LOWERBOUND of ENERGY");
				   seinfo(queryID,2*j+1,"This changes UPPERBOUND of ENERGY");
            }
			}
      }
      seloc(0, 0, 51, 18);	   seloc(0, 1, 87, 18);
      seloc(0, 2, 112, 18);	seloc(0, 3, 149, 18);
      seloc(0, 4, 174, 18);	seloc(0, 5, 209, 18);
      seloc(0, 6, 243, 18);	seloc(0, 7, 271, 18);
	   
      seloc(1, 0, 51, 33);    seloc(1, 1, 87, 33);
      seloc(1, 2, 112, 33);	seloc(1, 3, 149, 33);
      seloc(1, 4, 174, 33);	seloc(1, 5, 209, 33);
      seloc(1, 6, 243, 33);	seloc(1, 7, 271, 33);
	   
      seloc(2, 0, 51, 48);    seloc(2, 1, 87, 48);
      seloc(2, 2, 112, 48);	seloc(2, 3, 149, 48);
      seloc(2, 4, 174, 48);	seloc(2, 5, 209, 48);
      seloc(2, 6, 243, 48);	seloc(2, 7, 271, 48);
	}
}

//------------------------------------------------------------------------
// Hypothesis to be used for querying
//------------------------------------------------------------------------
void initHypothesesSmartSE_D3(graph* g) {
	int i, j, tlutid, queryID;
	float min,max;
	char queryString[1024], infoString[1024];

	//Add Queries to viewport
	if(g->numVertFeatures>0) {
		seformatdefault("%.0f");
		for(i=0;i<g->numQueries;i++) {

         if(g->fltQuery[i]!=NULL) free(g->fltQuery[i]);
			g->fltQuery[i]=(float*)malloc(2*g->numVertFeatures*sizeof(float));
			if(g->fltQuery[i]==NULL) { perror("initHypSmartSE()"); exit(-1);}
			for(j=0;j<g->numVertFeatures;j++) {
				getMinMax(g,j,&min,&max);
				g->fltQuery[i][2*j]=min;
				g->fltQuery[i][2*j+1]=max;
			}
			
			sprintf(queryString,"VQ%d =",i+1);
			for(j=0;j<g->numVertFeatures;j++) strcat(queryString," [%se,%se]");
			strcat(queryString,"\n");
			//queryID=seadd("Query = [%se,%se]\n",2*g->numVertFeatures,g->fltQuery[i]);
			queryID=seadd(queryString,2*g->numVertFeatures,g->fltQuery[i]);

			//Change seeEdit format
			tlutid=g->queryColors[i];
			secolor(queryID,tlut[tlutid][0],tlut[tlutid][1],tlut[tlutid][2]);
			for(j=0;j<g->numVertFeatures;j++) {
				min=g->fltQuery[i][2*j];
				max=g->fltQuery[i][2*j+1];
				seminmax(queryID,2*j,  min,max);
				seminmax(queryID,2*j+1,min,max);
				sestep(queryID,2*j,  (max-min)/100.f);
				sestep(queryID,2*j+1,(max-min)/100.f);
			}

			// Customized info for D3 - will change with each new query
			//Should be F0[1,778], F1[-1,43], F2[-111.4327,29475018],
			//          F3[-1,29476031.19], F4[-1,36], F5[0,.779], F6[0,0.074309373], F7[0,1]
			for(j=0;j<g->numVertFeatures;j++) {
				if(j>=5 && j<=7) {
					seformat(queryID,2*j,"%.2f");
					seformat(queryID,2*j+1,"%.2f");
				}

				if(j==0) {
               //TEST - VIDEO
               //strcpy(infoString,"The DEGREE of the vertex");
					strcpy(infoString,"The QTL LOCUS the gene is most closely related to");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==1) {
					strcpy(infoString,"The CHROMOSOME on which the gene resides");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==2) {
					strcpy(infoString,"MEGABASE location");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==3) {
					strcpy(infoString,"GENE MEGABASE location");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==4) {
					strcpy(infoString,"PARACLIQUE this gene is a member of");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==5) {
					strcpy(infoString,"HERITABILITY computed using RMA");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==6) {
					strcpy(infoString,"SQUARED ERROR of RMA");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==7) {
					strcpy(infoString,"P-VALUE of RMA");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}

			}

			//getDisplayEquation(g->numVertFeatures,g->fltQuery[i],&(g->strQuery[i]));
		}
	}
	else VCB_WARNING("No hypothesis to initialize: build tree data first\n");

	// Customized layout - will change with each new query
	seloc(0, 0, 59, 20);
	seloc(0, 1, 91, 20);    seloc(0, 2, 137, 20);   seloc(0, 3, 172, 20);
	seloc(0, 4, 207, 20);   seloc(0, 5, 246, 20);   seloc(0, 6, 287, 20);
	seloc(0, 7, 321, 20);   seloc(0, 8, 369, 20);   seloc(0, 9, 406, 20);
	seloc(0, 10, 444, 20);  seloc(0, 11, 479, 20);  seloc(0, 12, 522, 20);
	seloc(0, 13, 557, 20);  seloc(0, 14, 599, 20);  seloc(0, 15, 633, 20);

	seloc(1, 0, 59, 35);
	seloc(1, 1, 91, 35);    seloc(1, 2, 137, 35);   seloc(1, 3, 172, 35);
	seloc(1, 4, 207, 35);   seloc(1, 5, 246, 35);   seloc(1, 6, 287, 35);
	seloc(1, 7, 321, 35);   seloc(1, 8, 369, 35);   seloc(1, 9, 406, 35);
	seloc(1, 10, 444, 35);  seloc(1, 11, 479, 35);  seloc(1, 12, 522, 35);
	seloc(1, 13, 557, 35);  seloc(1, 14, 599, 35);  seloc(1, 15, 633, 35);

	seloc(2, 0, 59, 50);
	seloc(2, 1, 91, 50);    seloc(2, 2, 137, 50);   seloc(2, 3, 172, 50);
	seloc(2, 4, 207, 50);   seloc(2, 5, 246, 50);   seloc(2, 6, 287, 50);
	seloc(2, 7, 321, 50);   seloc(2, 8, 369, 50);   seloc(2, 9, 406, 50);
	seloc(2, 10, 444, 50);  seloc(2, 11, 479, 50);  seloc(2, 12, 522, 50);
	seloc(2, 13, 557, 50);  seloc(2, 14, 599, 50);  seloc(2, 15, 633, 50);
}

//------------------------------------------------------------------------
// Hypothesis to be used for querying
//------------------------------------------------------------------------
void initHypothesesSmartSE_D3c(graph* g) {
	int i, j, tlutid, queryID;
	float min,max;
	char queryString[1024], infoString[1024];

   //Change value to unique for all
   for(i=0;i<g->numVerts;i++) {
      g->vertData[i*2]=(float)i;
   }

	//Add Queries to viewport
	if(g->numVertFeatures>0) {
		seformatdefault("%.0f");
		for(i=0;i<g->numQueries;i++) {

         if(g->fltQuery[i]!=NULL) free(g->fltQuery[i]);
			g->fltQuery[i]=(float*)malloc(2*g->numVertFeatures*sizeof(float));
			if(g->fltQuery[i]==NULL) { perror("initHypSmartSE()"); exit(-1);}

			for(j=0;j<g->numVertFeatures;j++) {
				getMinMax(g,j,&min,&max);
				g->fltQuery[i][2*j]=min;
				g->fltQuery[i][2*j+1]=max;
			}
			
			sprintf(queryString,"VQ%d =",i+1);
			for(j=0;j<g->numVertFeatures;j++) strcat(queryString," [%se,%se]");
			strcat(queryString,"\n");
			//queryID=seadd("Query = [%se,%se]\n",2*g->numVertFeatures,g->fltQuery[i]);
			queryID=seadd(queryString,2*g->numVertFeatures,g->fltQuery[i]);

			//Change seeEdit format
			tlutid=g->queryColors[i];
			secolor(queryID,tlut[tlutid][0],tlut[tlutid][1],tlut[tlutid][2]);
			for(j=0;j<g->numVertFeatures;j++) {
				min=g->fltQuery[i][2*j];
				max=g->fltQuery[i][2*j+1];
				seminmax(queryID,2*j,  min,max);
				seminmax(queryID,2*j+1,min,max);
				sestep(queryID,2*j,  (max-min)/100.f);
				sestep(queryID,2*j+1,(max-min)/100.f);
			}

			// Customized info for D3 - will change with each new query
			//Should be F0[1,778], F1[-1,43], F2[-111.4327,29475018],
			//          F3[-1,29476031.19], F4[-1,36], F5[0,.779], F6[0,0.074309373], F7[0,1]
			for(j=0;j<g->numVertFeatures;j++) {
				if(j>=5 && j<=7) {
					seformat(queryID,2*j,"%.2f");
					seformat(queryID,2*j+1,"%.2f");
				}

				if(j==0) {
					strcpy(infoString,"The QTL LOCUS the gene is most closely related to");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==1) {
					strcpy(infoString,"The CHROMOSOME on which the gene resides");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==2) {
					strcpy(infoString,"MEGABASE location");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==3) {
					strcpy(infoString,"GENE MEGABASE location");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==4) {
					strcpy(infoString,"PARACLIQUE this gene is a member of");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==5) {
					strcpy(infoString,"HERITABILITY computed using RMA");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==6) {
					strcpy(infoString,"SQUARED ERROR of RMA");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}
				if(j==7) {
					strcpy(infoString,"P-VALUE of RMA");
					seinfo(queryID,2*j,infoString);
					seinfo(queryID,2*j+1,infoString);
				}

			}

			//getDisplayEquation(g->numVertFeatures,g->fltQuery[i],&(g->strQuery[i]));
		}
	}
	else VCB_WARNING("No hypothesis to initialize: build tree data first\n");

	// Customized layout - will change with each new query
	seloc(0, 0, 59, 20);
	seloc(0, 1, 91, 20);    seloc(0, 2, 137, 20);   seloc(0, 3, 172, 20);
	seloc(0, 4, 207, 20);   seloc(0, 5, 246, 20);   seloc(0, 6, 287, 20);
	seloc(0, 7, 321, 20);   seloc(0, 8, 369, 20);   seloc(0, 9, 406, 20);
	seloc(0, 10, 444, 20);  seloc(0, 11, 479, 20);  seloc(0, 12, 522, 20);
	seloc(0, 13, 557, 20);  seloc(0, 14, 599, 20);  seloc(0, 15, 633, 20);

	seloc(1, 0, 59, 35);
	seloc(1, 1, 91, 35);    seloc(1, 2, 137, 35);   seloc(1, 3, 172, 35);
	seloc(1, 4, 207, 35);   seloc(1, 5, 246, 35);   seloc(1, 6, 287, 35);
	seloc(1, 7, 321, 35);   seloc(1, 8, 369, 35);   seloc(1, 9, 406, 35);
	seloc(1, 10, 444, 35);  seloc(1, 11, 479, 35);  seloc(1, 12, 522, 35);
	seloc(1, 13, 557, 35);  seloc(1, 14, 599, 35);  seloc(1, 15, 633, 35);

	seloc(2, 0, 59, 50);
	seloc(2, 1, 91, 50);    seloc(2, 2, 137, 50);   seloc(2, 3, 172, 50);
	seloc(2, 4, 207, 50);   seloc(2, 5, 246, 50);   seloc(2, 6, 287, 50);
	seloc(2, 7, 321, 50);   seloc(2, 8, 369, 50);   seloc(2, 9, 406, 50);
	seloc(2, 10, 444, 50);  seloc(2, 11, 479, 50);  seloc(2, 12, 522, 50);
	seloc(2, 13, 557, 50);  seloc(2, 14, 599, 50);  seloc(2, 15, 633, 50);
}


//------------------------------------------------------------------------
// Initializes the vertex locations from the initial threshold's *.gvl
//------------------------------------------------------------------------
int initOneLocation(char* filenameGVL, graph* g) {
	int numVertLocs;
	view.locList=(short int**)malloc(1*sizeof(short int*));
   if(view.locList==NULL) {
      return -1;
   }

   numVertLocs=loadGVL(filenameGVL,&(view.locList[0]));
   if(numVertLocs < 0) {
      return -1;
   }

	printf("  Read %s\n",filenameGVL);
	if(g->numVerts != numVertLocs) {
		printf("GVL does not correspond with input file\n");
		exit(1);
	}
	VCB_DEBUG_CODE {
		printf("  Before normalization:\n");
		printf("    Vert 1 = (%d,%d,%d)\n",view.locList[0][0],view.locList[0][1],view.locList[0][2]);
		printf("    Vert 2 = (%d,%d,%d)\n",view.locList[0][3],view.locList[0][4],view.locList[0][5]);
		printf("    Vert %d = (%d,%d,%d)\n",g->numVerts-1,view.locList[0][3*g->numVerts-3],view.locList[0][3*g->numVerts-2],view.locList[0][3*g->numVerts-1]);
	}
	Normalize(g->numVerts,0,DIM_SIZE-1,&(view.locList[0]));
	VCB_DEBUG_CODE {
		printf("  After normalization:\n");
		printf("    Vert 1 = (%d,%d,%d)\n",view.locList[0][0],view.locList[0][1],view.locList[0][2]);
		printf("    Vert 2 = (%d,%d,%d)\n",view.locList[0][3],view.locList[0][4],view.locList[0][5]);
		printf("    Vert %d = (%d,%d,%d)\n",g->numVerts-1,view.locList[0][3*g->numVerts-3],view.locList[0][3*g->numVerts-2],view.locList[0][3*g->numVerts-1]);
	}
   return 0;
}

//------------------------------------------------------------------------
// Returns the current threshold since this is calculated and not stored
//------------------------------------------------------------------------
float getCurrentThresh() {return (startThresh+(float)displayThresh*stepThresh);}

//------------------------------------------------------------------------
// Cleanup function to deallocate all unused memory after initialization
//------------------------------------------------------------------------
void cleanup (void)
{
	vcbVolrCloseS();
	fflush(stdout);
}

//------------------------------------------------------------------------
// Calculates the current frame rate
//------------------------------------------------------------------------
void calculateFrameRate()
{
    static int framesPerSecond   = 0; // This will store our fps
    static long lastTime         = 0;    // This will hold the time from the last frame
    static char strFrameRate[50] = {0};  // We will store the string here for the window title
    
    // Increase the frame counter
    ++framesPerSecond;
	newTime=clock();
    if( newTime-oldTime >= CLOCKS_PER_SEC )
    {
        oldTime=newTime;
        // Copy the frames per second into a string to display in the window
        sprintf(strFrameRate, "Current Frames Per Second: %d", framesPerSecond);
        // Reset the frames per second
        framesPerSecond = 0;
    }

    // draw frame rate on screen
    drawText(strFrameRate, GLUT_BITMAP_HELVETICA_10,
			    textColor[0],textColor[1],textColor[2], 0.4f,0.003f );
}

//------------------------------------------------------------------------
// Draws the text to the screen
//------------------------------------------------------------------------
void drawText(char *str, void *font, GLclampf r, GLclampf g, GLclampf b,
              GLfloat x, GLfloat y) {
  /* font: font to use, e.g., GLUT_BITMAP_HELVETICA_10
     r, g, b: text colour
     x, y: text position in window: vcbRange [0,0] (bottom left of window)
           to [1,1] (top right of window). */

  char *ch;
  GLint matrixMode;
  GLboolean lightingOn;

 lightingOn= glIsEnabled(GL_LIGHTING);        /* lighting on? */
 if (lightingOn) glDisable(GL_LIGHTING);
 glGetIntegerv(GL_MATRIX_MODE, &matrixMode);  /* matrix mode? */
 glMatrixMode(GL_PROJECTION);
 glPushMatrix();
   glLoadIdentity();
   gluOrtho2D(0.0, 1.0, 0.0, 1.0);
   glMatrixMode(GL_MODELVIEW);
   glPushMatrix();
     glLoadIdentity();
     glPushAttrib(GL_COLOR_BUFFER_BIT);       /* save current colour */
       glColor3f(r, g, b);
       glRasterPos3f(x, y, 0.0);
       for(ch= str; *ch; ch++) {
          glutBitmapCharacter(font, (int)*ch);
       }
     glPopAttrib();
   glPopMatrix();
   glMatrixMode(GL_PROJECTION);
 glPopMatrix();
 glMatrixMode(matrixMode);
 if (lightingOn) glEnable(GL_LIGHTING);
}

//------------------------------------------------------------------------
// Converts world space into 2d screen location
//------------------------------------------------------------------------
void world2screen(float objX, float objY, float objZ, double* sx, double* sy) {
	double sz[1];
	double modelview[16];
	double projection[16];
	int viewport[4];

	glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
	glGetDoublev(GL_PROJECTION_MATRIX, projection);
	glGetIntegerv(GL_VIEWPORT, viewport);
	gluProject(objX,objY,objZ, modelview,projection,viewport, sx,sy,sz);
}

//------------------------------------------------------------------------
// Draws the labels for the given graph to the appropriate location
//------------------------------------------------------------------------
void drawVertLabels(graph* g) {
	int i;
	float red=0.9f,green=0.9f,blue=0.9f;
   //TEST - CLUT font color
   //float red=0.f,green=0.f,blue=0.9f;
	double x,y;
   char* label;

	for(i=0;i<g->numDisplayVerts;i++) {
		/*world2screen(g->offset[0]+g->vertPosn[g->vertList[i]*3  ],
					 g->offset[1]+g->vertPosn[g->vertList[i]*3+1],
					 g->offset[2]+g->vertPosn[g->vertList[i]*3+2],&x,&y);*/
      world2screen(g->vertPosn[g->vertList[i]*3  ],
					    g->vertPosn[g->vertList[i]*3+1],
					    g->vertPosn[g->vertList[i]*3+2],&x,&y);
		x = x/GL_WIN_WIDTH;
		y = y/GL_WIN_HEIGHT;

      if(g->vertLabelAlias==0) label=g->edgeLabels->list[g->vertList[i]];
      else label=g->vertLabels[g->vertList[i]];

		//drawText(label, GLUT_BITMAP_HELVETICA_10, red,green,blue, x,y);
		drawText(label, GLUT_BITMAP_HELVETICA_12,	red,green,blue, (GLfloat)x,(GLfloat)y);
      //TEST - big font
      //drawText(label, GLUT_BITMAP_HELVETICA_18,	red,green,blue, x,y);
	}
}

//------------------------------------------------------------------------
// Draw the current threshold to the screen
//------------------------------------------------------------------------
void drawCurrentThresh() {
    static char strThresh[50] = {0};  // We will store the string here for the window title
	if(viewOption==OPT_MULTIGVL || viewOption==OPT_SINGLEGVLSTEP) {
		if(!absoluteValue)	sprintf(strThresh, "Current Threshold: %1.4f", getCurrentThresh());
		else sprintf(strThresh, "Current Threshold: |%1.4f|", getCurrentThresh());
	}
	else {
		if(!absoluteValue)	sprintf(strThresh, "Current Threshold: %1.4f", threshold);
		else sprintf(strThresh, "Current Threshold: |%1.4f|", threshold);
	}
    drawText(strThresh, GLUT_BITMAP_HELVETICA_10,
			textColor[0],textColor[1],textColor[2], textX,textY );
}

//------------------------------------------------------------------------
// Draw the number of current vertices to the screen
//------------------------------------------------------------------------
void drawNumVerts() {
    static char strNumVerts[50] = {0};  // We will store the string here for the window title
	sprintf(strNumVerts, "Number of Vertices: %d", numDisplayVerts);
    //drawText(strNumVerts, GLUT_BITMAP_HELVETICA_10,
	//		textColor[0],textColor[1],textColor[2], textX,textY-textH ); //draw under
	drawText(strNumVerts, GLUT_BITMAP_HELVETICA_10,
			textColor[0],textColor[1],textColor[2], (GLfloat)(textX+0.2),(GLfloat)textY); //draw to right
}

//------------------------------------------------------------------------
// Draw the number of current edges to the screen
//------------------------------------------------------------------------
void drawNumEdges() {
    static char strNumEdges[50] = {0};  // We will store the string here for the window title
	sprintf(strNumEdges, "Number of Edges: %d", numDisplayEdges);
    //drawText(strNumEdges, GLUT_BITMAP_HELVETICA_10, 
			//textColor[0],textColor[1],textColor[2], textX,textY-2*textH ); //draw under
	drawText(strNumEdges, GLUT_BITMAP_HELVETICA_10, 
			textColor[0],textColor[1],textColor[2], (GLfloat)(textX+2*0.2),(GLfloat)textY ); //draw to right
}

//------------------------------------------------------------------------
// Draw the current vertex hypothesis to the screen
//------------------------------------------------------------------------
void drawVertHyp() {
	int i;
	float charHeight=0.02f;
	for(i=0;i<g->numQueries;i++) {
		drawText(g->strQuery[i], GLUT_BITMAP_HELVETICA_12, 
				tlut[g->queryColors[i]][0],tlut[g->queryColors[i]][1],tlut[g->queryColors[i]][2],
				(GLfloat)(textX+0.2),(GLfloat)(textY-i*1.5*charHeight));
	}
}

//------------------------------------------------------------------------
// Draw busy message to inform the user that we are computing
//------------------------------------------------------------------------
void drawBusy() {
    static char strWarning[50] = {0};  // We will store the string here for the window title
	sprintf(strWarning, "Please Wait...");
    drawText(strWarning, GLUT_BITMAP_HELVETICA_18, 1.0f, 1.0f, 0.0f, 0.43f, 0.5f );
    glFlush();
    glutSwapBuffers();
}

//------------------------------------------------------------------------
// Draws edges with colors based upon table lookup
//------------------------------------------------------------------------
void drawEdges(graph* g, float** tlut, float dataThresh, int normalize, int separate) {
	int i, colorIndex, v1, v2;
	
	glDisable(GL_LIGHTING);//Removes abrupt lighting change artifacts

	glLoadIdentity();
    glTranslatef (-view.eye[0], -view.eye[1], -(view.eye[2]-view.zoom*view.eye[2]));
	glTranslatef (g->offset[0],g->offset[1],g->offset[2]);
	glMultMatrixf(Spaceball.GetMatrix());
	glScalef(g->scale[0],g->scale[1],g->scale[2]);
	glRotatef(degreeRot,0.f,1.f,0.f);
	glTranslated(-DIM_SIZE/2,-DIM_SIZE/2,-DIM_SIZE/2);

	glBegin(GL_LINES);
	for(i=0;i<g->numDisplayEdges;i++) {
		v1=g->edgeList[2*i+0];
		v2=g->edgeList[2*i+1];
		colorIndex = getColorIndex(g->mat2d[v1][v2],dataThresh,normalize,separate);

		//Need to find v1 and v2 from edge index
		VCB_DEBUG_CODE {
			if(colorIndex>255 || colorIndex<0 || colorIndex<(int)(threshold*255.f))
				printf("Error on color index\n");
		}
		glColor4fv(tlut[colorIndex]);
		glVertex3f(g->vertPosn[v1*3],g->vertPosn[v1*3+1],g->vertPosn[v1*3+2]);
		glVertex3f(g->vertPosn[v2*3],g->vertPosn[v2*3+1],g->vertPosn[v2*3+2]);
	}
	glEnd();

	glEnable(GL_LIGHTING);//Removes abrupt lighting change artifacts
}

//------------------------------------------------------------------------
// Draw normals for our splats to see if they're behaving correctly
//------------------------------------------------------------------------
void drawHairySplats() {
	int i;
	unsigned short *pos;

	glLoadIdentity();
    glTranslatef (-view.eye[0], -view.eye[1], -(view.eye[2]-view.zoom*view.eye[2]));
	glMultMatrixf(Spaceball.GetMatrix());
	Spaceball.Update();
	glRotatef(degreeRot,0.f,1.f,0.f);
	glTranslated(-DIM_SIZE/2,-DIM_SIZE/2,-DIM_SIZE/2);

	for(i=0;i<g->numVerts;i++) {
		glBegin(GL_LINES);
			pos = (unsigned short*)&voxels[i*10+4];
			glVertex3f(pos[0],pos[1],pos[2]);
			glVertex3f(pos[0]+5.f*(float)voxels[i*10+1]/255.f,pos[1]+5.f*(float)voxels[i*10+2]/255.f,pos[2]+5.f*(float)voxels[i*10+3]/255.f);
		glEnd();
	}
}

//------------------------------------------------------------------------
// Draws all vertices as quadrics
//------------------------------------------------------------------------
void drawVertQuadrics(graph* g, float** tlut) {
	int i,vid;
	float mat[4] = {0.66f, 0.66f, 1.f, 1.f};


	glLoadIdentity();
   glTranslatef (-view.eye[0], -view.eye[1], -(view.eye[2]-view.zoom*view.eye[2]));
	glTranslatef (g->offset[0],g->offset[1],g->offset[2]);
	glMultMatrixf(Spaceball.GetMatrix());
	Spaceball.Update();
	glScalef(g->scale[0],g->scale[1],g->scale[2]);
	glRotatef(degreeRot,0.f,1.f,0.f);
	glTranslated(-DIM_SIZE/2,-DIM_SIZE/2,-DIM_SIZE/2);

	vcbSetOGLMaterial(GL_FRONT, MATERIAL_PEWTER, 1.0);
	glMaterialfv(GL_FRONT,GL_SPECULAR,mat);
	glMaterialf(GL_FRONT,GL_SHININESS, 2.5f);
	glShadeModel(GL_SMOOTH);

	//glEnable(GL_DEPTH_TEST);
	//glDepthFunc(GL_LEQUAL);

	for(i=0;i<g->numDisplayVerts;i++) {
		vid = g->vertList[i];
		glLoadName(vid);					// Assign Object A Name (ID)
		if( strcmp("LOD layout",g->name)==0) { //LOD graph, just pass through
			glColor4fv(tlut[g->vertColorList[vid]]);
			glTranslated(g->vertPosn[vid*3],g->vertPosn[vid*3+1],g->vertPosn[vid*3+2]);
			gluSphere(quadratic,1.f,resX,resY);
			glTranslated(-g->vertPosn[vid*3],-g->vertPosn[vid*3+1],-g->vertPosn[vid*3+2]);
		}
		else if( strcmp("subgraph",g->name)==0) { //subgraph, just pass through
			glColor4fv(tlut[g->vertColorList[vid]]);
			glTranslated(g->vertPosn[vid*3],g->vertPosn[vid*3+1],g->vertPosn[vid*3+2]);
			gluSphere(quadratic,1.f,resX,resY);
			glTranslated(-g->vertPosn[vid*3],-g->vertPosn[vid*3+1],-g->vertPosn[vid*3+2]);
		}
		else if( showFullVerts ) { //showing all, just pass through
			glColor4fv(tlut[g->vertColorList[vid]]);
			glTranslated(g->vertPosn[vid*3],g->vertPosn[vid*3+1],g->vertPosn[vid*3+2]);
			gluSphere(quadratic,1.f,resX,resY);
			glTranslated(-g->vertPosn[vid*3],-g->vertPosn[vid*3+1],-g->vertPosn[vid*3+2]);
		}
    else if( g->vertColorList[vid]==geneNodeColor || 
			g->vertColorList[vid]==geneNodeColorQ1 || 
			g->vertColorList[vid]==geneNodeColorQ2) {
			glColor4fv(tlut[g->vertColorList[vid]]);
			glTranslated(g->vertPosn[vid*3],g->vertPosn[vid*3+1],g->vertPosn[vid*3+2]);
			gluSphere(quadratic,1.f,resX,resY);
			glTranslated(-g->vertPosn[vid*3],-g->vertPosn[vid*3+1],-g->vertPosn[vid*3+2]);
		}
		else if (showDisconnected && g->vertColorList[vid]==geneTransparentColor) {
			glColor4fv(tlut[g->vertColorList[vid]]);
			glTranslated(g->vertPosn[vid*3],g->vertPosn[vid*3+1],g->vertPosn[vid*3+2]);
			gluSphere(quadratic,1.f,resX,resY);
			glTranslated(-g->vertPosn[vid*3],-g->vertPosn[vid*3+1],-g->vertPosn[vid*3+2]);
		}
		else if (g->vertColorList[vid]==geneNodeSelectColor) {
			glColor4fv(tlut[g->vertColorList[vid]]);
			glTranslated(g->vertPosn[vid*3],g->vertPosn[vid*3+1],g->vertPosn[vid*3+2]);
			gluSphere(quadratic,1.f,resX,resY);
			glTranslated(-g->vertPosn[vid*3],-g->vertPosn[vid*3+1],-g->vertPosn[vid*3+2]);
		}
		else if (showUnselected && !showDisconnected && g->vertColorList[vid]==geneNodeUnSelectColor) {
			glColor4fv(tlut[g->vertColorList[vid]]);
			//glColor4fv(tlut[geneTransparentColor]);
			glTranslated(g->vertPosn[vid*3],g->vertPosn[vid*3+1],g->vertPosn[vid*3+2]);
			gluSphere(quadratic,1.f,resX,resY);
			glTranslated(-g->vertPosn[vid*3],-g->vertPosn[vid*3+1],-g->vertPosn[vid*3+2]);
		}
		//else don't draw
		//can't remove the translations to below if statements
	}
	//glDisable(GL_DEPTH_TEST);
}

//------------------------------------------------------------------------
// Draw a single large sphere
//------------------------------------------------------------------------
void drawSphere() {
	glLoadIdentity();
    glTranslatef (-view.eye[0], -view.eye[1], -(view.eye[2]-view.zoom*view.eye[2]));
	glMultMatrixf(Spaceball.GetMatrix());
	Spaceball.Update();
	glRotatef(degreeRot,0.f,1.f,0.f);
	glTranslated(-DIM_SIZE/2,-DIM_SIZE/2,-DIM_SIZE/2);

	glTranslated(50,50,50);
	glColor4f(0.5f,0.5f,0.5f,1.f);
	gluSphere(quadratic,100.f,32,32);
}

//------------------------------------------------------------------------
// Draw our vertices as splats with hardware acceleration
//------------------------------------------------------------------------
void drawVertexSplats(unsigned char* voxels, graph* g, viewData* view) {
	glLoadIdentity();
    glTranslatef (-view->eye[0], -view->eye[1], -(view->eye[2]-view->zoom*view->eye[2]));
	glMultMatrixf(Spaceball.GetMatrix());
	Spaceball.Update();
	glScalef(g->scale[0],g->scale[1],g->scale[2]);
	glRotatef(degreeRot,0.f,1.f,0.f);
	glTranslated(-DIM_SIZE/2,-DIM_SIZE/2,-DIM_SIZE/2);

	vcbVolrSsplats(voxels, g->numVerts, dimSize);
}

//------------------------------------------------------------------------
// Draw halos for all the given vertices
//------------------------------------------------------------------------
void drawVertexHalos(unsigned char* voxels, graph* g, viewData* view) {
	int i;
	unsigned short* pos;

	glLoadIdentity();
    glTranslatef (-view->eye[0], -view->eye[1], -(view->eye[2]-view->zoom*view->eye[2]));
	glTranslatef (g->offset[0],g->offset[1],g->offset[2]);
	glMultMatrixf(Spaceball.GetMatrix());
	//Spaceball.Update();
	glScalef(g->scale[0],g->scale[1],g->scale[2]);
	glRotatef(degreeRot,0.f,1.f,0.f);
	glTranslated(-DIM_SIZE/2,-DIM_SIZE/2,-DIM_SIZE/2);

	//HACK: Draw custom halo
	for(i=0;i<g->numVerts;i++) voxels[i*10]=(unsigned char)0;
	for(i=0;i<g->numDisplayVerts;i++) {
		voxels[g->vertList[i]*10]=(unsigned char)g->haloColorIdx;
	}

	for(i=0;i<g->numVerts;i++) {
		pos = (unsigned short*)&(voxels[i*10+4]);
		(*pos++) = g->vertPosn[i*3  ]; // x-coordinate
		(*pos++) = g->vertPosn[i*3+1]; // y-coordinate
		(*pos)   = g->vertPosn[i*3+2]; // z-coordinate
		/*printf("(%d,%d,%d)->(%d,%d,%d)\n",g->offset[0]+g->vertPosn[i*3  ],
			g->offset[0]+g->vertPosn[i*3+1],g->offset[0]+g->vertPosn[i*3+2],
			pos[-2],pos[-1],pos[0]);*/
	}

	vcbVolrSsplats(voxels, g->numVerts, dimSize);

	for(i=0;i<g->numVerts;i++)
		voxels[i*10]=(unsigned char)g->vertColorList[i];
	for(i=0;i<g->numVerts;i++) {
		pos = (unsigned short*)&(voxels[i*10+4]);
		(*pos++) = view->locList[0][i*3+ 0]; // x-coordinate
		(*pos++) = view->locList[0][i*3+ 1]; // y-coordinate
		(*pos)   = view->locList[0][i*3+ 2]; // z-coordinate
	}
}

//------------------------------------------------------------------------
// Sets the OpenGL state machine so that it can draw normal stuff
//------------------------------------------------------------------------
void disableSplatting() {
	glColor3f(1.f,1.f,1.f);
	glActiveTexture(GL_TEXTURE0);
	glDisable(GL_TEXTURE_1D);
	glActiveTexture(GL_TEXTURE1);
	glDisable(GL_TEXTURE_2D);
	glEnable(GL_LIGHTING);
}

//------------------------------------------------------------------------
// Sets the OpenGL state machine so that texture-based splatting can work
//------------------------------------------------------------------------
void enableSplatting() {
	glColor3f(1.f,1.f,1.f);
	glDisable(GL_LIGHTING);
	glActiveTexture(GL_TEXTURE1);
	glEnable(GL_TEXTURE_2D);
	glActiveTexture(GL_TEXTURE0);
	glEnable(GL_TEXTURE_1D);
}

//------------------------------------------------------------------------
// Move last vertex in the specified direction
//------------------------------------------------------------------------
void moveVertex(graph* g, unsigned char key) {
	int vid = g->numDisplayVerts-1;
	if(key=='8') {				//up
		g->vertPosn[vid*3+1]++;
	} else if (key=='2') {		//down
		g->vertPosn[vid*3+1]--;
	} else if (key=='6') {		//right
		g->vertPosn[vid*3]++;
	} else if (key=='4') {      //left
		g->vertPosn[vid*3]--;
	} else if (key=='/') {      //back
		g->vertPosn[vid*3+2]--;
	} else if (key=='*') {		//forward
		g->vertPosn[vid*3+2]++;
	}
	/*printf(" Vertex %d = (%d,%d,%d)\n",vid,g->vertPosn[vid*3],
											g->vertPosn[vid*3+1],
											g->vertPosn[vid*3+2]);*/
}

//------------------------------------------------------------------------
// Move the graph in the specified direction
//------------------------------------------------------------------------
void moveGraph(graph* g, unsigned char key) {
	int vid = g->numDisplayVerts-1;
	int moveSpeed = 7;
	if(key=='8') {				//up
		g->offset[1]+=moveSpeed;
	} else if (key=='2') {		//down
		g->offset[1]-=moveSpeed;
	} else if (key=='6') {		//right
		g->offset[0]+=moveSpeed;
	} else if (key=='4') {      //left
		g->offset[0]-=moveSpeed;
	} else if (key=='/') {      //back
		g->offset[2]-=moveSpeed;
	} else if (key=='*') {		//forward
		g->offset[2]+=moveSpeed;
	}
	/*printf(" Graph Offset = (%d,%d,%d)\n",g->offset[0],
										   g->offset[1],
										   g->offset[2]);*/
}

//------------------------------------------------------------------------
// The desired threshold was changed, alter vertex locations appropriately
//------------------------------------------------------------------------
void updateVoxelPosition(int index, graph* g) {
	int i;
	unsigned short* pos;
	for(i=0;i<g->numVerts;i++) {
		pos = (unsigned short*)&voxels[i*10+4];
		pos[0]=view.locList[index][3*i+0];
		pos[1]=view.locList[index][3*i+1];
		pos[2]=view.locList[index][3*i+2];
	}
	updateDisplayVars(g);
}

//------------------------------------------------------------------------
// Update the display variables (#verts&edges)
//------------------------------------------------------------------------
void updateDisplayVars(graph* g) {
	clock_t startTime, endTime;
	startTime=clock();
	drawBusy();
	if (showOnlyLCC) {
		printf("Current threshold = %f\n",getCurrentThresh());
		numDisplayVerts=showLCC(getCurrentThresh(),absoluteValue,g);
		numDisplayEdges=inducedSubgraph(getCurrentThresh(),absoluteValue,g);
		updateVoxelVisibility(&voxels, geneTransparentColor, g);
		return;
	}
	else
		numDisplayVerts=filterDisconnectedVertsUsingDegree(g,displayThresh);

	updateVoxelVisibility(&voxels, geneTransparentColor, g);
	endTime=clock();
	VCB_DEBUG("  Initializing Vertices took %d ms\n",endTime-startTime);

	startTime=clock();
	
	//Visibility Option - show edges for only certain-colored vertices or all visible vertices
	if(!showEdgesVisibleVerts) numDisplayEdges=initVisibleEdgeList(geneNodeColor, g);
	else numDisplayEdges=initVisibleEdgeListAll(g);

	if(!showFullEdges) {
		numDisplayEdges=filterEdgesFailThresh(getCurrentThresh(),absoluteValue,g);
		//numDisplayEdges=RunEdgeQuery(getCurrentThresh(),absoluteValue,g);
	}
	//numDisplayEdges=filterInducedSubgraphShow(g, geneNodeColor);

	endTime=clock();
	VCB_INFO("  Initializing Edges took %d ms\n",endTime-startTime);
}

//------------------------------------------------------------------------
// Update splat normals every frame if desired
//------------------------------------------------------------------------
void updateSplatNormals() {
	int i;
	unsigned short* pos;
	float ModelView[16];
	float ModelViewInvT[16];
	float normalDir[4];

	glLoadIdentity();
    glTranslatef (-view.eye[0], -view.eye[1], -(view.eye[2]-view.zoom*view.eye[2]));
	glMultMatrixf(Spaceball.GetMatrix());
	Spaceball.Update();
	glTranslated(-DIM_SIZE/2,-DIM_SIZE/2,-DIM_SIZE/2);

	glGetFloatv (GL_MODELVIEW_MATRIX, ModelView);
	for(i=0;i<g->numVerts;i++) {
		pos = (unsigned short*)&voxels[i*10+4];
		vcbInvMat(ModelViewInvT,ModelView);
		vcbTransposeMatr(ModelView);
		vcbMatMult4f(normalDir, ModelViewInvT, (float)pos[0], (float)pos[1], (float)pos[2], 1.f);
		vcbNormalize3(&normalDir[0]);
		//printf("Normals %d = %f %f %f\n",i,normalDir[0],normalDir[1],normalDir[2]);
		voxels[i*10+1] = (unsigned char)normalDir[0];	// x-component of normal
		voxels[i*10+2] = (unsigned char)normalDir[1];	// y-component of normal
		voxels[i*10+3] = (unsigned char)normalDir[2];	// z-component of normal
	}
}

//------------------------------------------------------------------------
// Allows users to pick vertices
//------------------------------------------------------------------------
void pickVertices(int x, int y, int button) {
	int choose,depth,loop;
	GLuint selectB[128];		//storage buffer for hits
	GLint hits, viewport[4];	//num hits & viewport data
	glGetIntegerv(GL_VIEWPORT, viewport); //viewport x,y,length,width

	glSelectBuffer(128, selectB);	//tell OpenGL to use this as buffer
	glRenderMode(GL_SELECT); //store results in buffer, don't draw to screen
	glInitNames();					//initializes the name stack
	glPushName(0);					//pushes one entry onto the stack

	glMatrixMode(GL_PROJECTION);								// Selects The Projection Matrix
	glPushMatrix();												// Push The Projection Matrix
	glLoadIdentity();											// Resets The Matrix
	// This Creates A Matrix That Will Zoom Up To A Small Portion Of The Screen, Where The Mouse Is.
	gluPickMatrix((GLdouble) x, (GLdouble) (viewport[3]-y), 1.0f, 1.0f, viewport);
	// Apply The Perspective Matrix
	gluPerspective(45.0f, (GLfloat) (viewport[2]-viewport[0])/(GLfloat) (viewport[3]-viewport[1]), 0.1f, 1000.0f);
	glMatrixMode(GL_MODELVIEW);		// Select The Modelview Matrix
	// Render The Targets To The Selection Buffer
	if(viewMode==VIEW_QUADRICS) drawVertQuadrics(g,tlut);
	else if (viewMode==VIEW_SPLAT) drawVertexSplats(voxels,g,&view);
	else {printf("Unrecognized picking mode\n");}
	glMatrixMode(GL_PROJECTION);	// Select The Projection Matrix
	glPopMatrix();					// Pop The Projection Matrix
	glMatrixMode(GL_MODELVIEW);		// Select The Modelview Matrix
	hits=glRenderMode(GL_RENDER);	// Switch To Render Mode, Find Out How Many

	if (hits > 0) {									// If There Were Some Hits
		choose = selectB[3];					// Make Our Selection The First Object
		depth = selectB[1];						// Store How Far Away It Is 
		for (loop = 1; loop < hits; loop++) {	// Loop Through All The Detected Hits
			// If This Object Is Closer To Us Than The One We Have Selected
			if (selectB[loop*4+1] < GLuint(depth)) {
				choose = selectB[loop*4+3];			// Select The Closer Object
				depth = selectB[loop*4+1];			// Store How Far Away It Is
			}       
		}
		if(button==GLUT_LEFT_BUTTON) {
			if (g->vertColorList[choose]==geneNodeSelectColor) g->vertColorList[choose]=geneNodeColor;
			else g->vertColorList[choose]=geneNodeSelectColor;
		} else {
			if (g->vertColorList[choose]==geneNodeUnSelectColor) g->vertColorList[choose]=geneNodeColor;
			else g->vertColorList[choose]=geneNodeUnSelectColor;
		}
    }
}

//------------------------------------------------------------------------
// Unselect/label all vertices
//------------------------------------------------------------------------
void clearVerts() {
	int i;
	for(i=0;i<g->numVerts;i++) {
		g->vertColorList[i]=geneNodeColor;
	}
}

//------------------------------------------------------------------------
// Train neural net on selected vertices and find similar
//------------------------------------------------------------------------
int classify(NeuralNet* nn, Data* learn, Data* test) {
	int sum,numCorrect,minIdx=0;
	float minErr;
	clock_t start;

	//---Learning
	//printf("Learning ...\n");
	numCorrect=nn->nepochs=sum=0; minErr=RAND_MAX;
	//printf("  Epoch: %d, MSE: %f\n",nn->nepochs,nn->a.mse);
	start=clock();
	while ( numCorrect!=learn->nitems && nn->nepochs<MAX_EPOCHS) {
		//Go through all input data once (one epoch) and learn each time
		numCorrect = onlineNN(nn,learn,1);
		//printf("  Epoch: %d, MSE: %f (%d correct)\n",nn->nepochs,nn->a.mse,numCorrect);
		/*if(nn->nepochs%NUM_EPOCHS==0) {
			saveAgent(nn);
		}*/
	}

	//---Testing
	//printf("Testing ...\n");
	numCorrect=onlineNN(nn,test,0);
	//printf("  Training and classifiction took %0.6f seconds\n",(clock()-start)/(float)CLOCKS_PER_SEC);
	return numCorrect;
}

//------------------------------------------------------------------------
// Initialize data to be used by neural network (works only for degree)
//------------------------------------------------------------------------
int initLearnData(Data* d, graph* g, int allData) {
	int i,count;

	if(d->item != NULL) free(d->item);
	if(d->in != NULL) free(d->in);
	if(d->sup != NULL) free(d->sup);
			
	
	count=0;
	for(i=0;i<g->numVerts;i++) {
		if (allData==1) count++;
		else if(g->vertColorList[i]==geneNodeSelectColor || 
			    g->vertColorList[i]==geneNodeUnSelectColor)
			count++;
	}
	if(count==0) return -1;

	//Number of data items
	d->nitems=count;
	d->nin=numFeatures+1;	//Extra bias input
	d->nsup=2;				//Show or no show

	//Allocate space safely
	if(! (d->item=(int*)malloc(d->nitems*sizeof(int))) ) {
		perror("init item"); exit(-1);
	}
	if(! (d->in=(float**)malloc(d->nitems*sizeof(float*))) ) {
		perror("init in"); exit(-1);
	}
	if(! (d->sup=(float**)malloc(d->nitems*sizeof(float*))) ) {
		perror("init sup"); exit(-1);
	}
	//for each record, allocate space for input/output layers
	for(i=0;i<d->nitems;i++) {
		if(! (d->in[i]=(float*)malloc(d->nin*sizeof(float))) ) {
			perror("init in[]"); exit(-1);
		}
		if(! (d->sup[i]=(float*)malloc(d->nsup*sizeof(float))) ) {
			perror("init sup[]"); exit(-1);
		}
	}

	count=0;
	if(allData==1) {
		for(i=0;i<g->numVerts;i++) {
			d->item[count]=i;
			d->in[count][0]=(float)(g->degree[displayThresh][i]);
			d->in[count][d->nin-1]=1.f;
			//Don't care about supervisor because we'll classify these
			count++;
		}
	} else {
		for(i=0;i<g->numVerts;i++) {
			if(g->vertColorList[i]==geneNodeSelectColor ||
			   g->vertColorList[i]==geneNodeUnSelectColor) {
				d->item[count]=i;
				d->in[count][0]=(float)(g->degree[displayThresh][i]);
				d->in[count][d->nin-1]=1.f;
				if(g->vertColorList[i]==geneNodeSelectColor) {//Selected
					d->sup[count][0]=0.f;	d->sup[count][1]=1.f;
				} else if (g->vertColorList[i]==geneNodeUnSelectColor) {//UnSelected
					d->sup[count][0]=1.f;	d->sup[count][1]=0.f;
				}
				//printf(" Degree node %d=%.0f\n",i,d->in[count][0]);
				count++;
			}
		}
	}
	return 0;
}

//------------------------------------------------------------------------
// Initialize LOD data to be used by NN (works only for degree)
//------------------------------------------------------------------------
int initLearnDataLOD(Data* d, graph* g, int allData) {
	int i,j,v1,count;

	count=0;
	for(i=0;i<g->numVerts;i++) {
		if (allData==1)
			count+=g->numBin[i];
		else if(g->vertColorList[i]==geneNodeSelectColor || 
			    g->vertColorList[i]==geneNodeUnSelectColor)
			count+=g->numBin[i];
	}
	if(count==0) return -1;

	//Number of data items
	d->nitems=count;
	d->nin=numFeatures+1;	//Extra bias input
	d->nsup=2;				//Show or no show

	if(d->item != NULL) free(d->item);
	if(d->in != NULL) free(d->in);
	if(d->sup != NULL) free(d->sup);

	//Allocate space safely
	if(! (d->item=(int*)malloc(d->nitems*sizeof(int))) ) {
		perror("init item"); exit(-1);
	}
	if(! (d->in=(float**)malloc(d->nitems*sizeof(float*))) ) {
		perror("init in"); exit(-1);
	}
	if(! (d->sup=(float**)malloc(d->nitems*sizeof(float*))) ) {
		perror("init sup"); exit(-1);
	}
	//for each record, allocate space for input/output layers
	for(i=0;i<d->nitems;i++) {
		if(! (d->in[i]=(float*)malloc(d->nin*sizeof(float))) ) {
			perror("init in[]"); exit(-1);
		}
		if(! (d->sup[i]=(float*)malloc(d->nsup*sizeof(float))) ) {
			perror("init sup[]"); exit(-1);
		}
	}

	count=0;
	if(allData==1) {
		for(i=0;i<g->numVerts;i++) {
			for(j=0;j<g->numBin[i];j++) {
				v1=g->lodBin[i][j];
				d->item[count]=v1;
				d->in[count][0]=(float)g->degree[displayThresh][v1];
				d->in[count][d->nin-1]=1.f;
				//Don't care about supervisor because we'll classify these
				count++;
			}
		}
	} else {
		for(i=0;i<g->numVerts;i++) {
			if(g->vertColorList[i]==geneNodeSelectColor ||
			   g->vertColorList[i]==geneNodeUnSelectColor) {
				for(j=0;j<g->numBin[i];j++) {
					v1=g->lodBin[i][j];
					d->item[count]=v1;
					d->in[count][0]=(float)g->degree[displayThresh][v1];
					d->in[count][d->nin-1]=1.f;
					if(g->vertColorList[i]==geneNodeSelectColor) {//Selected
						//printf(" Ex: %d in bin %d\n",v1,i);
						d->sup[count][0]=0.f;	d->sup[count][1]=1.f;
					} else if (g->vertColorList[i]==geneNodeUnSelectColor) {//UnSelected
						//printf(" Cex: %d in bin %d\n",v1,i);
						d->sup[count][0]=1.f;	d->sup[count][1]=0.f;
					}
					//printf(" Degree node %d=%.0f\n",v1,d->in[count][0]);
					count++;
				}
			}
		}
	}
	return 0;
}

//------------------------------------------------------------------------
// Change voxel color based upon classification
//------------------------------------------------------------------------
void updateClassification(NeuralNet* n, graph* g) {
	int i, count, targetOnColor, targetOffColor;
	int remove = 1;

	if(remove==0) {
		targetOnColor=geneNodeSelectColor;
		targetOffColor=geneNodeUnSelectColor;
	}
	else {
		targetOnColor=geneNodeColor;
		targetOffColor=geneTransparentColor;
	}

	count=0;
	for(i=0;i<n->a.nguess;i++) {
		//This turns all vertices red but does not preserve query color
		// KI: NN classification of Subgraph vs. query color
		if(n->a.guess[i]==1) g->vertColorList[i]=targetOnColor;
		else if (n->a.guess[i]==0) {
			g->vertColorList[i]=-1;
			count++;
		}
	}

	if(remove==1) {
		numDisplayEdges = filterInducedSubgraphAllButColor(g, -1);
	}
	g->numDisplayVerts = g->numVerts-count;
	numDisplayVerts = g->numDisplayVerts;
	//printf("  Example = %d, CounterExample = %d\n",g->numVerts-count,count);
}

//------------------------------------------------------------------------
// The query has changed, so do what is necessary to update the data
//------------------------------------------------------------------------
int queryChanged() {
	int i;
	if( !showSeeEdit)
		for(i=0;i<g->numQueries;i++)
			getDisplayEquation(g->numVertFeatures,g->fltQuery[i],&(g->strQuery[i]));
	//getAllDisplayEquations(g->numQueries,g->numVertFeatures,g->fltQuery,&(g->strQuery));
	if(activeQuery) {
		drawBusy();
		numDisplayVerts=RunVertQuery(g);//numDisplayEdges=RunEdgeQuery(getCurrentThresh(),absoluteValue,g);
		updateDisplayVars(g);
	}
	return 0;
}

//------------------------------------------------------------------------
// Draw pixels over corner (better to use triangles) of BTD belt (testing)
//------------------------------------------------------------------------
void drawBTDCorners(BTD* btd) {
	int i,j,x,y;
	GLboolean lightingOn,texturingOn;

	if(btd->numSArea==0) return;

	lightingOn  = glIsEnabled(GL_LIGHTING);			// lighting on?
	texturingOn = glIsEnabled(GL_TEXTURE_2D);		// texturing on?
	if (lightingOn) glDisable(GL_LIGHTING);
	if (texturingOn) glDisable(GL_TEXTURE_2D);

	glMatrixMode(GL_PROJECTION);
	glPushMatrix(); 
	glLoadIdentity();
	glOrtho(0,GL_WIN_WIDTH,0,GL_WIN_HEIGHT,-1,1);
	glMatrixMode(GL_MODELVIEW);
	glPushMatrix(); 
	glLoadIdentity(); // Reset The Current Modelview Matrix

	// Now you are in 2D...
	glColor3f(1.f, 0.f, 0.f);
	glBegin(GL_POINTS);
	for(i=btd->drawX;i<btd->drawX+btd->drawWidth;i++) {
		for(j=btd->drawY;j<btd->drawX+btd->drawHeight;j++) {
			x=i;	y=j;
			screen2belt(btd, &x,&y);
			if( x<abs(y-btd->drawHeight/2) ||
				btd->drawWidth-x<abs(y-btd->drawHeight/2) ) {
				glVertex2f((GLfloat)i,(GLfloat)j);
			}
		}
	}
	glEnd();

	// Ending ortho mode...
	glPopMatrix();
	glMatrixMode( GL_PROJECTION );
	glPopMatrix();
	glMatrixMode( GL_MODELVIEW );
	if (texturingOn) glEnable(GL_TEXTURE_2D);
	if (lightingOn) glEnable(GL_LIGHTING);
}

//------------------------------------------------------------------------
// Show block tridiagonalized matrix as thin belt
//------------------------------------------------------------------------
void drawBTDBelt(BTD* btd) {
	GLboolean lightingOn;

	lightingOn  = glIsEnabled(GL_LIGHTING);			// lighting on?
	if (lightingOn) glDisable(GL_LIGHTING);

	glMatrixMode(GL_PROJECTION);
	glPushMatrix(); 
	glLoadIdentity();
	glOrtho(0,GL_WIN_WIDTH,0,GL_WIN_HEIGHT,-1,1);
	glMatrixMode(GL_MODELVIEW);
	glPushMatrix(); 
	glLoadIdentity(); // Reset The Current Modelview Matrix

	// Now you are in 2D...
	glTranslated(btd->drawX,btd->drawY,0);
	glScaled(btd->drawWidth, btd->drawHeight,1);

	glActiveTexture(btd->texUnit);
	glEnable(GL_TEXTURE_2D);

	glBegin(GL_QUADS);
		glMultiTexCoord2f(btd->texUnit, -btd->beltOffset, btd->beltOffset);
			glVertex3f(0.0f, 0.0f, 0.0f);
		glMultiTexCoord2f(btd->texUnit, 1.f-btd->beltOffset, 1.f+btd->beltOffset);
			glVertex3f(1.0f, 0.0f, 0.0f);
		glMultiTexCoord2f(btd->texUnit, 1.f+btd->beltOffset, 1.f-btd->beltOffset);
			glVertex3f(1.0f, 1.0f, 0.0f);
		glMultiTexCoord2f(btd->texUnit, btd->beltOffset, -btd->beltOffset);
			glVertex3f(0.0f, 1.0f, 0.0f);
	glEnd();

	glDisable(GL_TEXTURE_2D);

	// Ending ortho mode...
	glPopMatrix();
	glMatrixMode( GL_PROJECTION );
	glPopMatrix();
	glMatrixMode( GL_MODELVIEW );
	if (lightingOn) glEnable(GL_LIGHTING);
}

//------------------------------------------------------------------------
// Draw block tridiagonalized matrix selections over the texture
//------------------------------------------------------------------------
void drawBTDSelections(BTD* btd) {
	int i,offset;
	GLboolean lightingOn,texturingOn;

	if(btd->numSArea==0) return;

	lightingOn  = glIsEnabled(GL_LIGHTING);			// lighting on?
	texturingOn = glIsEnabled(GL_TEXTURE_2D);		// texturing on?
	if (lightingOn) glDisable(GL_LIGHTING);
	if (texturingOn) glDisable(GL_TEXTURE_2D);

	glMatrixMode(GL_PROJECTION);
	glPushMatrix(); 
	glLoadIdentity();
	glOrtho(0,GL_WIN_WIDTH,0,GL_WIN_HEIGHT,-1,1);
	glMatrixMode(GL_MODELVIEW);
	glPushMatrix(); 
	glLoadIdentity(); // Reset The Current Modelview Matrix

	// Now you are in 2D...
	glColor3f(0.f, 0.f, 0.f);
	for(i=0;i<btd->numSArea;i++) {
      if(btd->sArea[4*i+3]!=-1) {
		   glBegin(GL_LINE_LOOP);
			   glVertex2f((GLfloat)(btd->sArea[4*i]),  (GLfloat)(btd->sArea[4*i+1]));
			   glVertex2f((GLfloat)(btd->sArea[4*i+2]),(GLfloat)(btd->sArea[4*i+1]));
			   glVertex2f((GLfloat)(btd->sArea[4*i+2]),(GLfloat)(btd->sArea[4*i+3]));
			   glVertex2f((GLfloat)(btd->sArea[4*i]),  (GLfloat)(btd->sArea[4*i+3]));
		   glEnd();
      }
      else {
         offset=btd->sArea[4*i+2];
		   glBegin(GL_LINE_LOOP);
			   glVertex2f((GLfloat)(btd->sArea[4*i]-offset), (GLfloat)(btd->sArea[4*i+1]       ));
            glVertex2f((GLfloat)(btd->sArea[4*i]       ), (GLfloat)(btd->sArea[4*i+1]-offset));
            glVertex2f((GLfloat)(btd->sArea[4*i]+offset), (GLfloat)(btd->sArea[4*i+1]       ));
            glVertex2f((GLfloat)(btd->sArea[4*i]       ), (GLfloat)(btd->sArea[4*i+1]+offset));
		   glEnd();
      }
	}

	// Ending ortho mode...
	glPopMatrix();
	glMatrixMode( GL_PROJECTION );
	glPopMatrix();
	glMatrixMode( GL_MODELVIEW );
	if (texturingOn) glEnable(GL_TEXTURE_2D);
	if (lightingOn) glEnable(GL_LIGHTING);
}

//------------------------------------------------------------------------
// Filter to only show subgraph selected in BTD
//------------------------------------------------------------------------
void showOnlyBTD(BTD* btd, graph* g) {
	int i, index, rmFlag=-1, onFlag=255;

	if(btd->numSArea==0) return;

	// Free and allocate new memory for displayable constructs
	if(g->vertList!=NULL) free(g->vertList);
	if((g->vertList=(int*)malloc(g->numVerts*sizeof(int))) == NULL) {
		perror("showOnlyBTD() vertList"); exit(-1);
	}

	// Mark all vertices to not be included
	for(i=0;i<g->numVerts;i++) g->vertColorList[i]=rmFlag;

	// Add all edges to eventually be filtered
	showAllEdges(g);

	// Assign all BTD selections the given color index
	if(showBTDMatrixBelt) {
		btd->drawingBelt = 1;
		markBTDMatrixBelt(btd,g,onFlag);
	} else {
		btd->drawingBelt = 0;
		markBTDMatrix(btd,g,onFlag);
	}

	// Add flagged vertices to the new display list
	index=0;
	for(i=0;i<g->numVerts;i++) {
		if(g->vertColorList[i]==onFlag) {
			g->vertList[index]=i;
			index++;
		}
	}
	numDisplayVerts = g->numDisplayVerts = index;
	numDisplayEdges = filterInducedSubgraphAllButColor(g,rmFlag);
}
//------------------------------------------------------------------------
// Show block tridiagonalized matrix as quad texture
//------------------------------------------------------------------------
void drawTexture(GLuint tu, int x, int y, int width, int height) {
	GLboolean lightingOn;

	lightingOn  = glIsEnabled(GL_LIGHTING);			// lighting on?
	if (lightingOn) glDisable(GL_LIGHTING);

	glMatrixMode(GL_PROJECTION);
	glPushMatrix(); 
	glLoadIdentity();
	glOrtho(0,GL_WIN_WIDTH,0,GL_WIN_HEIGHT,-1,1);
	glMatrixMode(GL_MODELVIEW);
	glPushMatrix(); 
	glLoadIdentity(); // Reset The Current Modelview Matrix

	// Now you are in 2D...
	glTranslated(x,y,0);
	glScaled(width,height,1);

	glActiveTexture(tu);
	glEnable(GL_TEXTURE_2D);

	glBegin(GL_QUADS);
		glMultiTexCoord2f(tu, 1.0f, 0.0f); glVertex3f(0.0f, 0.0f, 0.0f);
		glMultiTexCoord2f(tu, 1.0f, 1.0f); glVertex3f(1.0f, 0.0f, 0.0f);
		glMultiTexCoord2f(tu, 0.0f, 1.0f); glVertex3f(1.0f, 1.0f, 0.0f);
		glMultiTexCoord2f(tu, 0.0f, 0.0f); glVertex3f(0.0f, 1.0f, 0.0f);
	glEnd();

	glDisable(GL_TEXTURE_2D);

	// Ending ortho mode...
	glPopMatrix();
	glMatrixMode( GL_PROJECTION );
	glPopMatrix();
	glMatrixMode( GL_MODELVIEW );
	if (lightingOn) glEnable(GL_LIGHTING);
}

//------------------------------------------------------------------------
// Draw Image directly to frame buffer
//------------------------------------------------------------------------
void drawImagef(float* pixels, float x, float y, int width, int height) {
     GLint matrixMode;
     GLboolean lightingOn, blendOn, atestOn;

     // setup opengl states for text, save modified states 
     if (lightingOn = glIsEnabled(GL_LIGHTING))
             glDisable(GL_LIGHTING);
     if (blendOn    = glIsEnabled(GL_BLEND))
             glDisable(GL_BLEND);
     if (atestOn = glIsEnabled(GL_ALPHA_TEST))
             glDisable(GL_ALPHA_TEST);
     glGetIntegerv(GL_MATRIX_MODE, &matrixMode);

     glPushAttrib(GL_COLOR_BUFFER_BIT);

     glMatrixMode(GL_MODELVIEW);
     glPushMatrix();
     glLoadIdentity();

     glMatrixMode(GL_PROJECTION);
     glPushMatrix();
     glLoadIdentity();
     glOrtho(0., 1., 0., 1., -1., 1.);
     /* end set up */

     glRasterPos3f(x, y, 0.f);

     glDrawPixels(width,height,GL_RGB,GL_FLOAT,pixels);

     /* restore opengl state machine */
     glPopAttrib();
     glPopMatrix();
     glMatrixMode(GL_MODELVIEW);
     glPopMatrix();
     glMatrixMode(matrixMode);
     if (lightingOn) glEnable(GL_LIGHTING);
     if (blendOn)    glEnable(GL_BLEND);
     if (atestOn)    glEnable(GL_ALPHA_TEST);

     glEnable(GL_BLEND);
     glEnable(GL_ALPHA_TEST);
}

//------------------------------------------------------------------------
// Draw the block tridiagonalized matrix
//------------------------------------------------------------------------
void drawBTDMatrix(BTD* btd, int x, int y, int width, int height) {
	float beltHeight;
   int maxDrawHeight = 150;
   //TEST - for marking
   //int maxDrawHeight = 250;

	btd->drawX=x; btd->drawY=y;
	btd->drawWidth=width; btd->drawHeight=height;

	//Draw the BTD texture
	if(showBTDMatrix && !showBTDMatrixBelt && showBTDMatrixStretch) {
		btd->drawWidth=btd->width; btd->drawHeight=btd->height;
		drawTexture(btd->texUnit, x,y, btd->drawWidth,btd->drawHeight);
		//drawImagef(btd.texData, x/(float)GL_WIN_WIDTH,y/(float)GL_WIN_HEIGHT, btd.width,btd.height);
	}
	else if(showBTDMatrix && showBTDMatrixBelt) {
		beltHeight = SQRT2*btd->maxBlockSize; //pythagorean theorem reduced
		if(beltHeight>maxDrawHeight/2) beltHeight=(float)maxDrawHeight;

      if(showBTDMatrixStretch) { //btd->height>GL_WIN_WIDTH-2*x || 
			btd->drawWidth  = GL_WIN_WIDTH-2*x;
			btd->drawHeight = (int)beltHeight;
		}
		else {
         //btd->drawWidth  = width;
			btd->drawWidth  = (int)(beltHeight*btd->maxBlockSize/btd->drawHeight); //match aspect ratio
			btd->drawHeight = (int)beltHeight;
		}
		drawBTDBelt(btd);
	}
	else if(showBTDMatrix) {
		drawTexture(btd->texUnit, btd->drawX, btd->drawY,
								  btd->drawWidth, btd->drawHeight);
		//drawImagef(btd.texData, x/(float)GL_WIN_WIDTH,y/(float)GL_WIN_HEIGHT, width,height);
	}

	//Draw the BTD selections
	if(showBTDMatrix && showBTDSelections) {
		drawBTDSelections(btd);
		//drawBTDRowLines(btd);
	}
}

//------------------------------------------------------------------------
// Store the given graph as a new subgraph in the array
//------------------------------------------------------------------------
void createSubgraph(graph* g) {
	graph* h;
	
	h = copy(g);

	if(h->name!=NULL) free(h->name);
	h->name=(char*)malloc( (strlen("subgraph")+1)*sizeof(char) );
	strcpy(h->name,"subgraph");

	//Hack: Hard-code position for the first 4 graphs
	if(numSGs==0) {
		h->offset[0]=285; h->offset[1]=128; h->offset[2]=-473;
	}
	else if(numSGs==1) {
		h->offset[0]=299; h->offset[1]=-114; h->offset[2]=-473;
	}
	else if(numSGs==2) {
		h->offset[0]=-320; h->offset[1]=-114; h->offset[2]=-473;
	}
	else if(numSGs==3) {
		h->offset[0]=-313; h->offset[1]=130; h->offset[2]=-473;
	}

	/*h->scale[0]=0.25f; h->scale[1]=0.25f; h->scale[2]=0.25f;*/
   h->haloColorIdx = (h->haloColorIdx+(numSGs+1)*58)%255; //HACK: same eqn createLOD
	
	/*for(i=0;i<h->numVerts;i++) {
		h->vertColorList[i]=0;
	}
	for(i=0;i<h->numDisplayVerts;i++) {
		h->vertColorList[h->vertList[i]] = h->haloColorIdx;
	}*/

	//Store the new graph
	numSGs++;
	if(numSGs==MAX_SUBGRAPHS) {perror("  Too many subgraphs"); exit(-1);}
	sg[numSGs-1]=h;
}

void takeScreenshot(char* filename) {
	if(pixels==NULL) {
		if( (pixels=(unsigned char*)malloc(GL_WIN_WIDTH*GL_WIN_HEIGHT*4*sizeof(unsigned char))) == NULL) {
			perror("takeScreenshot() pixels"); exit(-1);
		}
	}
	glReadPixels(0,0,GL_WIN_WIDTH,GL_WIN_HEIGHT,GL_RGBA,GL_UNSIGNED_BYTE,pixels);
	if(pixFile==NULL) {
		if( (pixFile=(char*)malloc( 10*sizeof(char))) == NULL) {
			perror("pixfile"); exit(-1);
		}
	}
	vcbImgWriteBMP(filename,pixels,4, GL_WIN_WIDTH, GL_WIN_HEIGHT);
	free(pixFile); pixFile=NULL;
	free(pixels);  pixels =NULL;
}

int copyFile(char* fromFile, char* toFile) {
  FILE *from, *to;
  char ch;

  /* open source file */
  if((from = fopen(fromFile, "rb"))==NULL) {
    printf("Cannot open source file.\n");
    exit(1);
  }

  /* open destination file */
  if((to = fopen(toFile, "wb"))==NULL) {
    printf("Cannot open destination file.\n");
    exit(1);
  }

  /* copy the file */
  while(!feof(from)) {
    ch = fgetc(from);
    if(ferror(from)) {
      printf("Error reading source file.\n");
      exit(1);
    }
    if(!feof(from)) fputc(ch, to);
    if(ferror(to)) {
      printf("Error writing destination file.\n");
      exit(1);
    }
  }

  if(fclose(from)==EOF) {
    printf("Error closing source file.\n");
    exit(1);
  }

  if(fclose(to)==EOF) {
    printf("Error closing destination file.\n");
    exit(1);
  }

  return 0;
}

void callR(graph* g) {

   //printDefaultWeightedEdgeList("R-input.txt", 2, (const void*)g->mat2d, (const char**)g->edgeLabels, (double)getCurrentThresh());
   write_graph ((char*)"R-input.txt", g, (char*)"veg"); //vertList, edgeList, graph LTM
   //copyFile("R-input.txt","R-output.txt");
   printf("  Write data to R-output.txt to change the current graph\n");
   system((char*)"R");
   read_graph ((char*)"R-output.txt",g,(char*)"veg");
   printf("  Data read into current working graph, continue interaction...\n");
   //when close, check for R-ouput.txt, give error notifications if not usable graph
}

//------------------------------------------------------------------------
// Initialize thermometer-like colorbar (need to add shader support to see)
//------------------------------------------------------------------------
GLubyte*** initColorBar(float** tlut) {
  int i,j,k;
  GLubyte*** tmpColorBar;

  tmpColorBar=(GLubyte***)malloc(512*sizeof(GLubyte**));
  if(tmpColorBar==NULL) {perror("initColorBar()"); return NULL;}
  for(i=0;i<512;i++) {
     tmpColorBar[i]=(GLubyte**)malloc(20*sizeof(GLubyte*));
     if(tmpColorBar[i]==NULL) {perror("initColorBar()"); return NULL;}
     for(j=0;j<20;j++) {
        tmpColorBar[i][j]=(GLubyte*)malloc(4*sizeof(GLubyte));
        if(tmpColorBar[i][j]==NULL) {perror("initColorBar()"); return NULL;}
     }
  }

  for(i=0;i<512;i++) for(j=0;j<20;j++) for(k=0;k<4;k++) {
	  tmpColorBar[i][j][k] = (GLubyte)(tlut[i/2][k]*255.f);
	  printf("%d %d= %u %u %u %u\n",i,j,tmpColorBar[i][j][0],
        tmpColorBar[i][j][1],tmpColorBar[i][j][2],tmpColorBar[i][j][3]);
  }
  return tmpColorBar;
}

//------------------------------------------------------------------------
//This function draws the LUT indices on the screen using a 2D colorbar
//------------------------------------------------------------------------
void oglDrawImage(float x, float y, GLubyte ***image, int width, int height, int blending) {
	GLint matrixMode;
	GLboolean lightingOn, blendOn, atestOn;

	/* setup opengl states for image output, saving modified states too */
	if (lightingOn = glIsEnabled(GL_LIGHTING))
		glDisable(GL_LIGHTING);
	if (blendOn  = glIsEnabled(GL_BLEND))
		if (blending) glDisable(GL_BLEND);
	if (atestOn = glIsEnabled(GL_ALPHA_TEST))
		glDisable(GL_ALPHA_TEST);
	glGetIntegerv(GL_MATRIX_MODE, &matrixMode);

	glPushAttrib(GL_COLOR_BUFFER_BIT);

	glMatrixMode(GL_MODELVIEW);
	glPushMatrix();
	glLoadIdentity();

	glMatrixMode(GL_PROJECTION);
	glPushMatrix();
	glLoadIdentity();
	glOrtho(0., 1., 0., 1., -1., 1.);
	/* end set up */

	glRasterPos3f(x, y, 0.f);
	glDrawPixels(width, height, GL_RGBA, GL_UNSIGNED_BYTE, image);

	/* restore opengl state machine */
	glPopAttrib();
	glPopMatrix();
	glMatrixMode(GL_MODELVIEW);
	glPopMatrix();
	glMatrixMode(matrixMode);
	if (lightingOn) glEnable(GL_LIGHTING);
	if (blendOn)    glEnable(GL_BLEND);
	if (atestOn)    glEnable(GL_ALPHA_TEST);

	glEnable(GL_BLEND);
	glEnable(GL_ALPHA_TEST);
}
