
#ifndef CORRELATION_H
#define CORRELATION_H

#define MAX_G_ARRAY 1024
#define BRUSH_MASK 1024
#define MAX_AUX 65536
#define MAX_COR_DIMS 64
#define MAX_N_SAVED 128
#define MAX_INTERDIM 16
#define MAX_CONTROLS 8

class SavedInterdim {

public:

	int dim1;
	int dim2;
	int trendflag; // linear, non-linear, const, erratic
	float first[2];
	float last[2];
	float error;

	SavedInterdim() { };

	void Copy(const SavedInterdim &other) {
		dim1 = other.dim1;
		dim2 = other.dim2;
		trendflag = other.trendflag;
		first[0] = other.first[0];
		first[1] = other.first[1];
		last[0] = other.last[0];
		last[1] = other.last[1];
		error = other.error;
	};

};


class SavedCorrelation {

public:

	int category;
	int dimflags[MAX_COR_DIMS]; // inc, dcr, const, erratic
	float dimpoints[MAX_COR_DIMS][8][2]; // for trending
	int num_dimpoints[MAX_COR_DIMS]; // for trending
	float dimpointerrors[MAX_COR_DIMS][8]; // for trending/constant
	SavedInterdim idims[MAX_INTERDIM];
	int nidims;

	SavedCorrelation() : nidims(0) { };

	void Print(FILE *fptr) {
		int i;
		fprintf(fptr,"class %d\n", category);
		for (i=0;i<Group::ndims;i++) {
			if (dimflags[i]==0) {
				fprintf(fptr,"%s linear, start %f, gradient %f\n",
					R.dimensionnames[i], dimpoints[i][0][0],
					dimpoints[i][1][0] - dimpoints[i][0][0]);
			} else if (dimflags[i]==1) {
				fprintf(fptr,"%s non-linear\n", R.dimensionnames[i]);
			} else if (dimflags[i]==2) {
				fprintf(fptr,"%s constant = %f\n", R.dimensionnames[i], dimpoints[i][0][0]);
			} else {
				fprintf(fptr,"%s erratic\n", R.dimensionnames[i]);
			}
		}
		for (i=0;i<nidims;i++) {
			fprintf(fptr,"%s trends with %s\n", R.dimensionnames[idims[i].dim1], R.dimensionnames[idims[i].dim2]);
		}

	};


	void Save(FILE *fptr) {
		int i, j;
		fprintf(fptr,"class %d\n", category);
		for (i=0;i<Group::ndims;i++) {
			fprintf(fptr,"%d\n", dimflags[i]);
			if (dimflags[i]==0) { // linear
				fprintf(fptr,"%e %e %e %e %e\n", dimpoints[i][0][0], dimpoints[i][1][0], 
					dimpoints[i][0][1], dimpoints[i][1][1], dimpointerrors[i][0]);
			} else if (dimflags[i]==1) { // non-linear
				fprintf(fptr,"%d\n", num_dimpoints[i]);
				for (j=0;j<num_dimpoints[i];j++) {
					fprintf(fptr,"%e %e %e\n", dimpoints[i][j][0], dimpoints[i][j][1], 
						dimpointerrors[i][j]);
				}
			} else if (dimflags[i]==2) { // constant
				fprintf(fptr,"%e %e\n", dimpoints[i][0][0], dimpointerrors[i][0]);
			}
		}
		fprintf(fptr,"%d\n", nidims); // inter-dimensional trends
		for (i=0;i<nidims;i++) {
			fprintf(fptr,"%d %d %d %e %e %e %e %e\n", idims[i].dim1, idims[i].dim2, idims[i].trendflag, 
				idims[i].first[0], idims[i].first[1], idims[i].last[0], idims[i].last[1], 
				idims[i].error); 
		}
	};


	void Read(FILE *fptr) {
		int i, j;
		char dum[32];
		fscanf(fptr,"%s %d", dum, &category);
		for (i=0;i<Group::ndims;i++) {
			fscanf(fptr,"%d",&(dimflags[i]));
			if (dimflags[i]==0) { // linear
				fscanf(fptr,"%e %e %e %e %e", &(dimpoints[i][0][0]), &(dimpoints[i][1][0]), 
					&(dimpoints[i][0][1]), &(dimpoints[i][1][1]), &(dimpointerrors[i][0]));
				if (dimpointerrors[i][0] < 0.00000001) dimpointerrors[i][0] = 0.00000001;
			} else if (dimflags[i]==1) { // non-linear
				fscanf(fptr,"%d", &(num_dimpoints[i]));
				for (j=0;j<num_dimpoints[i];j++) {
					fscanf(fptr,"%e %e %e", &(dimpoints[i][j][0]), &(dimpoints[i][j][1]), 
						&(dimpointerrors[i][j]));
					if (dimpointerrors[i][j] < 0.00000001) dimpointerrors[i][j] = 0.00000001;
				}
			} else if (dimflags[i]==2) { // constant
				fscanf(fptr,"%e %e", &(dimpoints[i][0][0]), &(dimpointerrors[i][0]));
				if (dimpointerrors[i][0] < 0.00000001) dimpointerrors[i][0] = 0.00000001;
			}
		}
		fscanf(fptr,"%d", &(nidims)); // inter-dimensional trends
		for (i=0;i<nidims;i++) {
			fscanf(fptr,"%d %d %d %e %e %e %e %e", &(idims[i].dim1), &(idims[i].dim2), &(idims[i].trendflag), 
				&(idims[i].first[0]), &(idims[i].first[1]), &(idims[i].last[0]), &(idims[i].last[1]), 
				&(idims[i].error)); 
		}

	};


};


class Correlation {

	Group *garray[MAX_G_ARRAY];
	Group *auxiliary[MAX_AUX];
	int SGAinds[MAX_G_ARRAY];
	int auxSGAinds[MAX_AUX];
	float coords[MAX_G_ARRAY][2];
	unsigned char brushed[MAX_G_ARRAY];
	unsigned char removed[MAX_G_ARRAY];
	float auxcoords[MAX_AUX][2];
	unsigned char auxbrushed[MAX_AUX];
	unsigned char auxremoved[MAX_AUX];
	float minx, maxx, miny, maxy;
	int dimflags[MAX_COR_DIMS]; // inc, dcr, const, erratic
	int numgraphrows;
	float dimextremes[MAX_COR_DIMS][2];
	int dimgraphcoords[MAX_COR_DIMS][2];
	unsigned char brushmask[BRUSH_MASK][BRUSH_MASK];
	int screenw, screenh;
	float dimpoints[MAX_COR_DIMS][MAX_CONTROLS][2]; // for trending
	int num_dimpoints[MAX_COR_DIMS]; // for trending
	float dimpointerrors[MAX_COR_DIMS][MAX_CONTROLS]; // for trending/constant

	int interdim[MAX_COR_DIMS]; // interdimensional correlation
	int interdim_trend; // inc, dcr, const, erratic
	

public:

	int ng;
	int naux;

	int correlation_mode;
	int interdim_mode;
	float margin;
	int negative;

	int nsaved;
	SavedCorrelation saved[MAX_N_SAVED];

	int nidims;
	SavedInterdim idims[MAX_INTERDIM];

	Correlation() : ng(0) , naux(0) , correlation_mode(0) , interdim_mode(0) ,
		margin(0.001) , negative(0) , nsaved(0) , interdim_trend(0) { };

	void Print(FILE *fptr) {
		int i;
		fprintf(fptr, "num correlations %d\n", nsaved);
		for (i=0;i<nsaved;i++) {
			fprintf(fptr,"Saved correlation %d\n", i);
			saved[i].Print(fptr);
		}
	};

	
	void SaveAll(FILE *fptr) {
		int i;
		fprintf(fptr, "num correlations %d\n", nsaved);
		for (i=0;i<nsaved;i++) {
			fprintf(fptr,"Saved correlation %d\n", i);
			saved[i].Save(fptr);
		}
	};

	void ReadAll(FILE *fptr, Group *gpp, int ng) {
		int i;
		char dum[32];
		// nullify flags
		for (i=0;i<ng;i++) obj_clusters[i] = (signed char)-1;
		fscanf(fptr,"%s %s %d", dum, dum, &nsaved);
		for (i=0;i<nsaved;i++) {
			fscanf(fptr,"%s %s %s", dum, dum, dum);
			saved[i].Read(fptr);
		}
		TestSaved(gpp, ng); // to set the flags
	};

	void Insert(Group *gptr, int ind, float x, float y) {
		int i;
		int found;
		if ((ng==0)&&(naux==0)) {
			minx = maxx = x;
			miny = maxy = y;
		} else {
			if (x<minx) minx = x;
			if (x>maxx) maxx = x;
			if (y<miny) miny = y;
			if (y>maxy) maxy = y;
		}
		found = 0;
		for (i=0;i<ng;i++) {
			if ((coords[i][0]>x-0.002)&&(coords[i][0]<x+0.002)&&
				(coords[i][1]>y-0.002)&&(coords[i][1]<y+0.002)) {
				found=1;
			}
		}
		if (found==1) {
			auxcoords[naux][0] = x;
			auxcoords[naux][1] = y;
			auxiliary[naux] = gptr;
			auxSGAinds[naux] = ind;
			naux++;
			if (naux >= MAX_AUX-1) exit(1);
		} else {
			coords[ng][0] = x;
			coords[ng][1] = y;
			garray[ng] = gptr;
			SGAinds[ng] = ind;
			ng++;
			if (ng >= MAX_G_ARRAY-1) exit(1);
		}
	};

	void Sort() {
		int ind;
		float tx, ty;
		int i,j;
		int tSGAind;
		Group *tgptr;

		ind = 1;
		if ((maxx - minx) > (maxy - miny)) ind = 0;
		for (i=0;i<ng;i++) {
			for (j=i+1;j<ng;j++) {
				if (coords[j][ind] < coords[i][ind]) {
					// swap
					tgptr = garray[j];
					tx = coords[j][0];
					ty = coords[j][1];
					tSGAind = SGAinds[j];
					garray[j] = garray[i];
					coords[j][0] = coords[i][0];
					coords[j][1] = coords[i][1];
					SGAinds[j] = SGAinds[i];
					garray[i] = tgptr;
					coords[i][0] = tx;
					coords[i][1] = ty;
					SGAinds[i] = tSGAind;
				}
			}
		}
	};

	void Correlate() {
		int i,j,n;
		int ind;
		int numicr, numdcr, numzro;
		int delta;
		int cont;

		ind = 1;
		if ((maxx - minx) > (maxy - miny)) ind = 0;
		
		for (i=0;i<ng;i++) {
			if (removed[i]==(unsigned char)0) {
				printf("%f ", coords[i][ind]);
			}
		}
		printf("\n");
		for (n=0;n<Group::ndims;n++) {
			for (i=0;i<ng;i++) {
				if (removed[i]==(unsigned char)0) {
					printf("%f ",garray[i]->startingcoords[n]);
				}
			}
			printf("\n");
		}
		
		for (n=0;n<Group::ndims;n++) {
			numicr = numdcr = numzro = 0;
			for (i=0;i<ng-1;i++) {
				delta = 1;
				cont = 1;
				// find the next index whose screen coords is significantly different
				if (removed[i]==(unsigned char)0) {
					for (j=i+1;(j<ng)&&(cont);j++) {
						if ((coords[j][ind] > coords[i][ind] + margin) && 
							(removed[j]==(unsigned char)0)) {
							cont = 0;	
						} else {
							//delta = delta * 2;
							if (j+delta >= ng) {
								cont = 0;
								j = -1;
							}
						}
					}
					j -= delta;
				
				
					if (j>0) {
						if (garray[j]->startingcoords[n] - garray[i]->startingcoords[n] > 0.000001) {
							numicr++;
						} else if (garray[j]->startingcoords[n] - garray[i]->startingcoords[n] < -0.000001) {
							numdcr++;
						} else {
							numzro++;
						}
					}
				}
			}
			printf("%d %d %d ", numicr, numdcr, numzro);
			if ((numdcr==0) && (numzro==0)) {
				printf("dimension %d strictly increasing\n", n);
				dimflags[n] = 0;
			} else if ((numicr==0) && (numzro==0)) {
				printf("dimension %d strictly decreasing\n", n);
				dimflags[n] = 0;
			} else if ((numicr==0) && (numdcr==0)) {
				printf("dimension %d all constant\n", n);
				dimflags[n] = 2;
			} else {
				printf("dimension %d erratic\n", n);
				dimflags[n] = 3;
			}
		}

		for (i=0;i<Group::ndims;i++) {
			if (dimflags[i]<2) {
				for (j=i+1;j<Group::ndims;j++) {
					if (dimflags[j]==0) {
						// found dimensions i and j that show trend
						// now lets see if their ratios are constant
						printf("dim %d and %d: ", i, j); 
						for (n=0;n<ng;n++) {
							printf("%f ", garray[n]->startingcoords[i]/garray[n]->startingcoords[j]);
						}
						printf("\n");
					}
				}
			}
		}

		// determine num rows and cols
		for (numgraphrows=0;numgraphrows*numgraphrows<Group::ndims;numgraphrows++);
		// for each dimension, determine i,j coord of graph
		for (i=0;i<Group::ndims;i++) {
			dimgraphcoords[i][0] = i%numgraphrows;
			dimgraphcoords[i][1] = i/numgraphrows;
		}
		// for each dimension, determine max,min value
		for (i=0;i<Group::ndims;i++) {
			dimextremes[i][0] = dimextremes[i][1] = garray[0]->startingcoords[i];
			for (j=0;j<ng;j++) {
				if (garray[j]->startingcoords[i]<dimextremes[i][0]) {
					dimextremes[i][0] = garray[j]->startingcoords[i];
				} else if (garray[j]->startingcoords[i]>dimextremes[i][1]) {
					dimextremes[i][1] = garray[j]->startingcoords[i];
				}
			}
			for (j=0;j<naux;j++) {
				if (auxiliary[j]->startingcoords[i]<dimextremes[i][0]) {
					dimextremes[i][0] = auxiliary[j]->startingcoords[i];
				} else if (auxiliary[j]->startingcoords[i]>dimextremes[i][1]) {
					dimextremes[i][1] = auxiliary[j]->startingcoords[i];
				}
			}
		}

		for (i=0;i<BRUSH_MASK;i++) {
			for (j=0;j<BRUSH_MASK;j++) {
				brushmask[i][j] = (unsigned char)0;
			}
		}

		for (i=0;i<MAX_G_ARRAY;i++) {
			brushed[i] = (unsigned char)0;
		}

		for (i=0;i<MAX_AUX;i++) {
			auxbrushed[i] = (unsigned char)0;
		}
	};

	void Clear() {
		int i;
		ng = 0; naux = 0;
		for (i=0;i<Group::ndims;i++) {
			interdim[i] = 0;
		}
	};

	void PreRemove() {
		int i;
		for (i=0;i<ng;i++) removed[i] = (unsigned char)0;
		for (i=0;i<naux;i++) auxremoved[i] = (unsigned char)0;
	};

	void Remove() {
		int i,j;
		for (i=0;i<ng;i++) {
			if (brushed[i]==(unsigned char)1) {
				removed[i] = (unsigned char)1;
				brushed[i] = (unsigned char)0;
			}
		}
		for (i=0;i<naux;i++) {
			if (auxbrushed[i]==(unsigned char)1) {
				auxremoved[i] = (unsigned char)1;
				auxbrushed[i] = (unsigned char)0;
			}
		}
		for (i=0;i<BRUSH_MASK;i++) {
			for (j=0;j<BRUSH_MASK;j++) {
				brushmask[i][j] = (unsigned char)0;
			}
		}
		Correlate();
	};

	void GetFirstLast(Group **first, Group **last, 
		float &firstcoord, float &lastcoord, int ind) {

		int i;
		int validfound;

		validfound = 0;
		for (i=0;i<ng;i++) {
			if (removed[i]==0) {
				if (validfound==0) {
					*first = *last = garray[i];
					firstcoord = lastcoord = coords[i][ind];
					validfound = 1;
				} else {
					if (coords[i][ind] < firstcoord) {
						firstcoord = coords[i][ind];
						*first = garray[i];
					}
					if (coords[i][ind] > lastcoord) {
						lastcoord = coords[i][ind];
						*last = garray[i];
					}
				}
			}
		}
		for (i=0;i<naux;i++) {
			if (auxremoved[i]==0) {
				if (validfound==0) {
					*first = *last = auxiliary[i];
					firstcoord = lastcoord = auxcoords[i][ind];
					validfound = 1;
				} else {
					if (auxcoords[i][ind] < firstcoord) {
						firstcoord = auxcoords[i][ind];
						*first = auxiliary[i];
					}
					if (auxcoords[i][ind] > lastcoord) {
						lastcoord = auxcoords[i][ind];
						*last = auxiliary[i];
					}
				}
			}
		}

	};

	void Finalize() {
		int i;
		int ind;
		Group *first;
		Group *last;
		float firstcoord, lastcoord;
		int n;
		float diff, error, bestvalue;

		ind = 1;
		if ((maxx - minx) > (maxy - miny)) ind = 0;
		
		GetFirstLast(&first, &last, firstcoord, lastcoord, ind);
		
		for (i=0;i<Group::ndims;i++) {
			if (dimflags[i]==0) {
				dimpoints[i][0][0] = first->startingcoords[i];
				dimpoints[i][0][1] = firstcoord;
				dimpoints[i][1][0] = last->startingcoords[i];
				dimpoints[i][1][1] = lastcoord;
				num_dimpoints[i] = 2;
				// set the error
				error = 0.0;
				for (n=0;n<ng;n++) {
					if (removed[n]==0) {
						bestvalue = first->startingcoords[i] + 
							((coords[n][ind] - firstcoord) / (lastcoord - firstcoord)) *
							(last->startingcoords[i] - first->startingcoords[i]);
						diff = garray[n]->startingcoords[i] - bestvalue;
						if (diff < 0.0) diff = -diff;
						if (diff > error) error = diff;
					}
				}
				for (n=0;n<naux;n++) {
					if (auxremoved[n]==0) {
						bestvalue = first->startingcoords[i] + 
							((auxcoords[n][ind] - firstcoord) / (lastcoord - firstcoord)) *
							(last->startingcoords[i] - first->startingcoords[i]);
						diff = auxiliary[n]->startingcoords[i] - bestvalue;
						if (diff < 0.0) diff = -diff;
						if (diff > error) error = diff;
					}
				}
				dimpointerrors[i][0] = error * (lastcoord - firstcoord) /
					(last->startingcoords[i] - first->startingcoords[i]);
				if (dimpointerrors[i][0] < 0.0) dimpointerrors[i][0] = -dimpointerrors[i][0];
				dimpointerrors[i][0] = 2.0 * dimpointerrors[i][0];
				printf("dimension %d trend, first %f last %f error %f\n", i, dimpoints[i][0][1], dimpoints[i][1][1], dimpointerrors[i][0]);
			} else if (dimflags[i]==2) {
				dimpoints[i][0][0] = first->startingcoords[i];
				dimpointerrors[i][0] = 0.0;
				for (n=0;n<ng;n++) {
					if (removed[n]==0) {
						error = garray[n]->startingcoords[i] - dimpoints[i][0][0];
						if (error < 0.0) error = -error;
						if (error > dimpointerrors[i][0]) dimpointerrors[i][0] = error;
					}
				}
				for (n=0;n<naux;n++) {
					if (auxremoved[n]==0) {
						error = auxiliary[n]->startingcoords[i] - dimpoints[i][0][0];
						if (error < 0.0) error = -error;
						if (error > dimpointerrors[i][0]) dimpointerrors[i][0] = error;
					}
				}
				dimpointerrors[i][0] = 2.0 * dimpointerrors[i][0];
				printf("dimension %d constant, error %f\n", i, dimpointerrors[i][0]);
			}
		}

	};

	void NonLinearTrend(int d) {
		int i, j, n;
		int ind;
		Group *first;
		Group *last;
		float firstcoord, lastcoord;
		int worst_ind;
		float worst_error;
		int worst_interval;
		int found;
		int interval;
		float interpolated, diff;
		int t;

		ind = 1;
		if ((maxx - minx) > (maxy - miny)) ind = 0;

		// insert the first and the last
		GetFirstLast(&first, &last, firstcoord, lastcoord, ind);
		dimpoints[d][0][0] = first->startingcoords[d];
		dimpoints[d][0][1] = firstcoord;
		dimpoints[d][1][0] = last->startingcoords[d];
		dimpoints[d][1][1] = lastcoord;
		num_dimpoints[d] = 2;
		printf("points ");
		for (t=0;t<num_dimpoints[d];t++) {
			printf("%f %f ,", dimpoints[d][t][0], dimpoints[d][t][1]);
		}
		printf("\n");
		// insert 6 others
		for (n=0;n<MAX_CONTROLS-2;n++) {		
			worst_ind = -1;
			worst_error = 0.0;
			for (i=0;i<ng;i++) {
				if (removed[i]==0) {
					// find the interval this point belongs to
					found = 0;
					for (j=1;((j<num_dimpoints[d]) && (!found));j++) {
						if ((dimpoints[d][j-1][1] < coords[i][ind]) &&
							(dimpoints[d][j][1] > coords[i][ind])) {
							found = 1;
							interval = j;
						}
					}
					if (found) {
						// interpolate the interval
						interpolated = dimpoints[d][interval-1][0] + 
							((coords[i][ind]-dimpoints[d][interval-1][1])/
							 (dimpoints[d][interval][1]-dimpoints[d][interval-1][1]))*
							(dimpoints[d][interval][0]-dimpoints[d][interval-1][0]);
						// compare with actual coords
						diff = interpolated - garray[i]->startingcoords[d];
						if (diff < 0.0) diff = -diff;
						// update worst_ind and worst_error if necessary
						if (diff > worst_error) {
							worst_error = diff;
							worst_ind = i;
							worst_interval = interval;
						}
					}
				}
			}
			// add the worst ind and shift rest up
			if (worst_ind >= 0) {
				for (i=num_dimpoints[d];i>worst_interval;i--) {
					dimpoints[d][i][0] = dimpoints[d][i-1][0];
					dimpoints[d][i][1] = dimpoints[d][i-1][1];
				}
				dimpoints[d][worst_interval][0] = garray[worst_ind]->startingcoords[d];
				dimpoints[d][worst_interval][1] = coords[worst_ind][ind];
				num_dimpoints[d]++;
			}
		}
		// set error margin
		dimpointerrors[d][0] = 0.1 * (lastcoord - firstcoord);
		dimpointerrors[d][1] = 0.1 * (dimextremes[d][1] - dimextremes[d][0]);
	};

	void DrawNonLinear(int d, float startingx, float startingy, 
						float graphw, float graphh, int ind) {

		int i;
		float px, py;

		glPointSize(3);
		glLineWidth(1);
		glColor3f(0.0,1.0,1.0);
		for (i=0;i<num_dimpoints[d];i++) {
			
			if (ind==0) {
				px = 0.1 + 0.8 * (dimpoints[d][i][1] - minx) / (maxx - minx);
			} else {
				px = 0.1 + 0.8 * (dimpoints[d][i][1] - miny) / (maxy - miny);
			}
			py = 0.1 + 0.8 * (dimpoints[d][i][0] - dimextremes[d][0]) /
						(dimextremes[d][1] - dimextremes[d][0]);
			px = startingx + px * graphw;
			py = startingy + py * graphh;
			glBegin(GL_POINTS);
			glVertex3f(px,py,0.0);
			glEnd();
			if (i<num_dimpoints[d]-1) {
				glBegin(GL_LINES);
				glVertex3f(px,py,0.0);
				if (ind==0) {
					px = 0.1 + 0.8 * (dimpoints[d][i+1][1] - minx) / (maxx - minx);
				} else {
					px = 0.1 + 0.8 * (dimpoints[d][i+1][1] - miny) / (maxy - miny);
				}
				py = 0.1 + 0.8 * (dimpoints[d][i+1][0] - dimextremes[d][0]) /
							(dimextremes[d][1] - dimextremes[d][0]);
				px = startingx + px * graphw;
				py = startingy + py * graphh;
				glVertex3f(px,py,0.0);
				glEnd();
			}
		}
	};

	void SelfTest() {
		int i;
		int ntested;
		int n;
		n = ntested = 0;
		for (i=0;i<ng;i++) {
			if (removed[i]==(unsigned char)0) {
				ntested++;
				n += Test(garray[i]);
			}
		}
		printf("success rate %d out of %d\n", n, ntested);
		n = ntested = 0;
		for (i=0;i<naux;i++) {
			if (auxremoved[i]==(unsigned char)0) {
				ntested++;
				n += Test(auxiliary[i]);
			}
		}
		printf("success rate %d out of %d\n", n, ntested);
	};


	int Test(Group *gptr) {
		int i, j;
		int d1, d2;
		float predictedscreen, firstpredictedscreen;
		int found, interval;
		float upperd, lowerd, upperscreen, lowerscreen, diff1, diff2;
		int okay = 1;
		int firstfound = 0;
		for (i=0;(i<Group::ndims)&&okay;i++) {
			if (dimflags[i]==2) { // const dim
				if ((gptr->startingcoords[i]>dimpoints[i][0][0]+dimpointerrors[i][0])||
					(gptr->startingcoords[i]<dimpoints[i][0][0]-dimpointerrors[i][0])) {
					okay = 0;
				}
			} else if (dimflags[i]==0) { // linear trend
				predictedscreen = dimpoints[i][0][1] + 
					((gptr->startingcoords[i]-dimpoints[i][0][0])/(dimpoints[i][1][0]-dimpoints[i][0][0]))*
					(dimpoints[i][1][1]-dimpoints[i][0][1]);
				if (firstfound==0) {
					firstfound = 1;
					firstpredictedscreen = predictedscreen;
					if (predictedscreen > dimpoints[i][1][1] + 0.01 * (dimpoints[i][1][1]-dimpoints[i][0][1])) okay = 0;
					if (predictedscreen < dimpoints[i][0][1] - 0.01 * (dimpoints[i][1][1]-dimpoints[i][0][1])) okay = 0;
				} else {
					if ((predictedscreen > firstpredictedscreen + dimpointerrors[i][0]) ||
						(predictedscreen < firstpredictedscreen - dimpointerrors[i][0])) {
						okay = 0;
					}
				}
			}
		}
		
		if (firstfound!=0) {
			for (i=0;(i<Group::ndims)&&okay;i++) {
				if (dimflags[i]==1) { // non-linear trend
					// add and subtract error margin to get upper and lower y bounds
					upperd = gptr->startingcoords[i] + dimpointerrors[i][1];
					lowerd = gptr->startingcoords[i] - dimpointerrors[i][1];

					// translate to upper and lower x bounds
					
					found = 0;
					for (interval=1;((interval<num_dimpoints[i]) && (!found));interval++) {
						if (((dimpoints[i][interval][0] >= upperd)&&
							 (dimpoints[i][interval-1][0] <= upperd)) ||
							 ((dimpoints[i][interval][0] <= upperd)&&
							 (dimpoints[i][interval-1][0] >= upperd))) found = 1;
					}
					if (found) {
						// find predicted screen coord by linear interpolation of interval
						upperscreen = dimpoints[i][interval-1][1] + 
								((upperd-dimpoints[i][interval-1][0])/(dimpoints[i][interval][0]-dimpoints[i][interval-1][0]))*
							(dimpoints[i][interval][1]-dimpoints[i][interval-1][1]);
					} else {
						// find predicted screen coord by comparing whether first or last point is closer
						diff1 = upperd-dimpoints[i][0][0];
						if (diff1 < 0.0) diff1 = -diff1;
						diff2 = upperd-dimpoints[i][num_dimpoints[i]-1][0];
						if (diff2 < 0.0) diff2 = -diff2;
						if (diff1 < diff2) {
							upperscreen = dimpoints[i][0][1];
						} else {
							upperscreen = dimpoints[i][num_dimpoints[i]-1][1];
						}
					}

					found = 0;
					for (interval=1;((interval<num_dimpoints[i]) && (!found));interval++) {
						if (((dimpoints[i][interval][0] >= lowerd)&&
							 (dimpoints[i][interval-1][0] <= lowerd)) ||
							 ((dimpoints[i][interval][0] <= lowerd)&&
							 (dimpoints[i][interval-1][0] >= lowerd))) found = 1;
					}
					if (found) {
						// find predicted screen coord by linear interpolation of interval
						lowerscreen = dimpoints[i][interval-1][1] + 
								((lowerd-dimpoints[i][interval-1][0])/(dimpoints[i][interval][0]-dimpoints[i][interval-1][0]))*
							(dimpoints[i][interval][1]-dimpoints[i][interval-1][1]);
					} else {
						// find predicted screen coord by comparing whether first or last point is closer
						diff1 = upperd-dimpoints[i][0][0];
						if (diff1 < 0.0) diff1 = -diff1;
						diff2 = upperd-dimpoints[i][num_dimpoints[i]-1][0];
						if (diff2 < 0.0) diff2 = -diff2;
						if (diff1 < diff2) {
							lowerscreen = dimpoints[i][0][1];
						} else {
							lowerscreen = dimpoints[i][num_dimpoints[i]-1][1];
						}
					}

					// check if predictedscreen is within these bounds
					if ((predictedscreen > upperscreen)||(predictedscreen < lowerscreen)) okay = 0;

				}
			}
		}
		// test inter-dimensional trends
		for (i=0;i<nidims;i++) {
			d1 = idims[i].dim1;
			d2 = idims[i].dim2;
			if (idims[i].trendflag==0) { // linear trend
				predictedscreen = idims[i].first[0] + ((gptr->startingcoords[d2]-idims[i].first[1])/(idims[i].last[1]-idims[i].first[1]))*
					(idims[i].last[0]-idims[i].first[0]);
				if ((predictedscreen > gptr->startingcoords[d2] + idims[i].error) ||
					(predictedscreen < gptr->startingcoords[d2] - idims[i].error)) {
					okay = 0;
				}
			} else if (idims[i].trendflag==1) { // non-linear trend
			}
		}
		return okay;
	};


	void DisplayPoint(int g, unsigned char *removedptr, unsigned char *brushedptr, Group **ggptr,
						float startingx, float startingy, float graphw, float graphh, int ind, int n, float coord) {
		float px, py;
		if (removedptr[g]==(unsigned char)0) {
			if (ind==0) {
				px = 0.1 + 0.8 * (coord - minx) / (maxx - minx);
			} else {
				px = 0.1 + 0.8 * (coord - miny) / (maxy - miny);
			}
			if (dimflags[n]==0) {
				glColor3f(1.0,1.0,0.0);
				py = 0.1 + 0.8 * (ggptr[g]->startingcoords[n] - dimextremes[n][0]) /
						(dimextremes[n][1] - dimextremes[n][0]);
			} else if ((dimflags[n]==0)||(dimflags[n]==1)) {
				glColor3f(1.0,1.0,1.0);
				py = 0.1 + 0.8 * (ggptr[g]->startingcoords[n] - dimextremes[n][0]) /
						(dimextremes[n][1] - dimextremes[n][0]);
			} else if (dimflags[n]==3) {
				glColor3f(0.0,1.0,0.0);
				py = 0.1 + 0.8 * (ggptr[g]->startingcoords[n] - dimextremes[n][0]) /
						(dimextremes[n][1] - dimextremes[n][0]);
			} else {
				glColor3f(0.0,0.0,1.0);
				py = 0.1 + ggptr[g]->startingcoords[n] * 0.8;
			}
			px = startingx + px * graphw;
			py = startingy + py * graphh;
			if (brushmask[(int)px][(int)py]==(unsigned char)1) {
				brushedptr[g] = (unsigned char)1;
			}
			glPointSize(1);
			glBegin(GL_POINTS);
			glVertex3f(px,py,0.0);
			glEnd();
		}

	};

	
	void DisplayHighPoint(int g, unsigned char *removedptr, unsigned char *brushedptr, Group **ggptr,
						float startingx, float startingy, float graphw, float graphh, int ind, int n, float coord) {
		float px, py;
		if ((brushedptr[g]==(unsigned char)1)&&
			(removedptr[g]==(unsigned char)0)) {
			glColor3f(1.0,0.0,0.0);
			if (ind==0) {
				px = 0.1 + 0.8 * (coord - minx) / (maxx - minx);
			} else {
				px = 0.1 + 0.8 * (coord - miny) / (maxy - miny);
			}
			if ((dimflags[n]==0)||(dimflags[n]==1)) {
				py = 0.1 + 0.8 * (ggptr[g]->startingcoords[n] - dimextremes[n][0]) /
						(dimextremes[n][1] - dimextremes[n][0]);
			} else if (dimflags[n]==3) {
				py = 0.1 + 0.8 * (ggptr[g]->startingcoords[n] - dimextremes[n][0]) /
						(dimextremes[n][1] - dimextremes[n][0]);
			} else {
				py = 0.1 + ggptr[g]->startingcoords[n] * 0.8;
			}
			px = startingx + px * graphw;
			py = startingy + py * graphh;
			glPointSize(3);
			glBegin(GL_POINTS);
			glVertex3f(px,py,0.0);
			glEnd();
		}
	};


	void Display(int sw, int sh) {
		if (interdim_mode==0) {
			NormalDisplay(sw,sh);
		} else {
			InterDimensionalDisplay(sw, sh);
		}
	};


	void NormalDisplay(int sw, int sh) {

		int i,j,n,g,c;
		int ind;
		float startingx, startingy;
		float graphw, graphh;
		float px, py;
		char trendname[32];


		screenw = sw;
		screenh = sh;

		ind = 1;
		if ((maxx - minx) > (maxy - miny)) ind = 0;
		graphw = screenw / (float) numgraphrows;
		graphh = screenh / (float) numgraphrows;
		

		glMatrixMode(GL_PROJECTION);
		glLoadIdentity();
		gluOrtho2D(0.0,screenw,0.0,screenh);
				// left, right, bottom, top
		glViewport(0.0,0.0,screenw,screenh);
				// startx, starty, xsize, ysize
				// coordinates begin from lower left corner of window	
		glMatrixMode(GL_MODELVIEW);

		for (i=0;i<screenw;i++) {
			for (j=0;j<screenh;j++) {
				if (brushmask[i][j]==(unsigned char)1) {
					glPointSize(1);
					glColor3f(0.5,0.5,0.5);
					glBegin(GL_POINTS);
					glVertex3f(i,j,-0.1);
					glEnd();
				}
			}
		}

		for (i=0;i<ng;i++) {
			brushed[i] = (unsigned char)0;
		}
		for (i=0;i<naux;i++) {
			auxbrushed[i] = (unsigned char)0;
		}

		for (n=0;n<Group::ndims;n++) {
			i = dimgraphcoords[n][0];
			j = dimgraphcoords[n][1];
			startingx = graphw * (float)i;
			startingy = graphh * (float)j;
			
			glLineWidth(1);
			glColor3f(0.5,0.5,0.5);
			glBegin(GL_LINES);
			glVertex3f(startingx+0.1*graphw,startingy+0.1*graphh,-0.1);
			glVertex3f(startingx+0.9*graphw,startingy+0.1*graphh,-0.1);
			glVertex3f(startingx+0.1*graphw,startingy+0.1*graphh,-0.1);
			glVertex3f(startingx+0.1*graphw,startingy+0.9*graphh,-0.1);
			glEnd();

			if (dimflags[n]==0) {
				glColor3f(1.0,1.0,0.0);
				strcpy(trendname,"linear trend");
			} else if (dimflags[n]==1) {
				glColor3f(1.0,1.0,1.0);
				strcpy(trendname,"non-linear trend");
			} else if (dimflags[n]==3) {
				glColor3f(0.0,1.0,0.0);
				strcpy(trendname,"erratic");
			} else {
				glColor3f(0.0,0.0,1.0);
				strcpy(trendname,"constant");
			}
			glBegin(GL_QUADS);
			glVertex3f(startingx+0.02*graphw,startingy+0.02*graphh,0.0);
			glVertex3f(startingx+0.08*graphw,startingy+0.02*graphh,0.0);
			glVertex3f(startingx+0.08*graphw,startingy+0.08*graphh,0.0);
			glVertex3f(startingx+0.02*graphw,startingy+0.08*graphh,0.0);
			glEnd();
	
			for (c=0;c<strlen(trendname);c++) {
				glRasterPos2f(startingx+0.1*graphw+5.0*(float)c,startingy);
				glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10,trendname[c]);
			}		
			
			for (g=0;g<ng;g++) {
				DisplayPoint(g, removed, brushed, garray,
						startingx, startingy, graphw, graphh, ind, n, coords[g][ind]);
			}
			for (g=0;g<naux;g++) {
				DisplayPoint(g, auxremoved, auxbrushed, auxiliary,
						startingx, startingy, graphw, graphh, ind, n, auxcoords[g][ind]);
			}

			if (interdim[n]==0) {
				glColor3f(1.0,1.0,1.0);
			} else {
				glColor3f(1.0,0.0,0.0);
			}
			glBegin(GL_QUADS);
			glVertex3f(startingx+0.52*graphw,startingy+0.02*graphh,0.0);
			glVertex3f(startingx+0.58*graphw,startingy+0.02*graphh,0.0);
			glVertex3f(startingx+0.58*graphw,startingy+0.08*graphh,0.0);
			glVertex3f(startingx+0.52*graphw,startingy+0.08*graphh,0.0);
			glEnd();

		}
		// second pass to highlight brushed points
		for (n=0;n<Group::ndims;n++) {
			i = dimgraphcoords[n][0];
			j = dimgraphcoords[n][1];
			startingx = graphw * (float)i;
			startingy = graphh * (float)j;
			for (g=0;g<ng;g++) {
				DisplayHighPoint(g, removed, brushed, garray,
						startingx, startingy, graphw, graphh, ind, n, coords[g][ind]);
			}
			for (g=0;g<naux;g++) {
				DisplayHighPoint(g, auxremoved, auxbrushed, auxiliary,
						startingx, startingy, graphw, graphh, ind, n, auxcoords[g][ind]);
			}
					
			// draw the non-linear trends
			if (dimflags[n]==1) DrawNonLinear(n, startingx, startingy, graphw, graphh, ind);
	
		}

		glutSwapBuffers(); 
	};

	void InterDimensionalDisplay(int sw, int sh) {

		int i;
		int d1, d2;
		int c;
		char trendname[32];

		d1 = d2 = -1;
		for (i=0;i<Group::ndims;i++) {
			if (interdim[i]==1) {
				if (d1 < 0) {
					d1 = i;
				} else {
					d2 = i;
				}
			}
		}

		glMatrixMode(GL_PROJECTION);
		glLoadIdentity();
		gluOrtho2D(-0.2,1.2,-0.2,1.2);
				// left, right, bottom, top
		glViewport(0.0,0.0,sw,sh);
				// startx, starty, xsize, ysize
				// coordinates begin from lower left corner of window	
		glMatrixMode(GL_MODELVIEW);


		
		if (interdim_trend==0) {
			glColor3f(1.0,1.0,0.0);
			strcpy(trendname,"linear trend");
		} else if (interdim_trend==1) {
			glColor3f(1.0,1.0,1.0);
			strcpy(trendname,"non-linear trend");
		} else if (interdim_trend==3) {
			glColor3f(0.0,1.0,0.0);
			strcpy(trendname,"erratic");
		} else {
			glColor3f(0.0,0.0,1.0);
			strcpy(trendname,"constant");
		}
		glBegin(GL_QUADS);
		glVertex3f(-0.02,-0.02,0.0);
		glVertex3f(-0.02,-0.01,0.0);
		glVertex3f(-0.01,-0.01,0.0);
		glVertex3f(-0.01,-0.02,0.0);
		glEnd();
	
		for (c=0;c<strlen(trendname);c++) {
			glRasterPos2f(0.02*(float)c,-0.02);
			glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_10,trendname[c]);
		}		


		glColor3f(1.0,1.0,0.0);
		glPointSize(2);
		glBegin(GL_POINTS);

		for (i=0;i<ng;i++) {
			if (removed[i]==0) {
				glVertex3f(garray[i]->startingcoords[d1],
					garray[i]->startingcoords[d2],0.0);
			}
		}

		for (i=0;i<naux;i++) {
			if (auxremoved[i]==0) {
				glVertex3f(auxiliary[i]->startingcoords[d1],
					auxiliary[i]->startingcoords[d2],0.0);
			}
		}

		glEnd();

		glutSwapBuffers(); 
	};

	void MouseMotion(int brushsize, int x, int y) {
		
		int bsx, bsy, bex, bey;
		int i, j;
		unsigned char b;

		bsx = x - brushsize; if (bsx<0) bsx = 0;
		bsy = y - brushsize; if (bsy<0) bsy = 0;
		bex = x + brushsize; if (bex>=screenw) bex = screenw-1;
		bey = y + brushsize; if (bey>=screenh) bey = screenh-1;
		b = (unsigned char)1;
		if (negative) b = (unsigned char)0;
		for (i=bsx;i<=bex;i++) {
			for (j=bsy;j<=bey;j++) {
				brushmask[i][j] = (unsigned char)b;
			}
		}
		
	};


	void MouseButton(int brushsize, int state, int x, int y) {

		int bsx, bsy, bex, bey;
		int i, j, n;
		unsigned char b;
		int ind;
		float startingx, startingy;
		float graphw, graphh;
		int changetrend = 0;

		ind = 1;
		if ((maxx - minx) > (maxy - miny)) ind = 0;
		graphw = screenw / (float) numgraphrows;
		graphh = screenh / (float) numgraphrows;

		if (state==GLUT_DOWN) {

			if (interdim_mode==0) {

				for (n=0;n<Group::ndims;n++) {
					i = dimgraphcoords[n][0];
					j = dimgraphcoords[n][1];
					startingx = graphw * (float)i;
					startingy = graphh * (float)j;
					if (((float)x>startingx+0.01*graphw)&&
						((float)x<startingx+0.09*graphw)&&
						((float)y>startingy+0.01*graphh)&&
						((float)y<startingy+0.09*graphh)) {
						changetrend = 1;
						if (dimflags[n]<3) {
							dimflags[n]++;
						} else {
							dimflags[n] = 0;
						}
						if (dimflags[n]==1) NonLinearTrend(n);
					}
					
					if (((float)x>startingx+0.51*graphw)&&
						((float)x<startingx+0.59*graphw)&&
						((float)y>startingy+0.01*graphh)&&
						((float)y<startingy+0.09*graphh)) {
						changetrend = 1;
						interdim[n]= 1 - interdim[n];
					}

				}
				if (!changetrend) {
					bsx = x - brushsize; if (bsx<0) bsx = 0;
					bsy = y - brushsize; if (bsy<0) bsy = 0;
					bex = x + brushsize; if (bex>=screenw) bex = screenw-1;
					bey = y + brushsize; if (bey>=screenh) bey = screenh-1;
					b = (unsigned char)1;
					if (negative) b = (unsigned char)0;
					for (i=bsx;i<=bex;i++) {
						for (j=bsy;j<=bey;j++) {
							brushmask[i][j] = (unsigned char)b;
						}
					}
				}

			} else {

				if (interdim_trend<3) {
					interdim_trend++;
				} else {
					interdim_trend = 0;
				}

			}
		}
	};

	int TestSavedOne(Group *gpp, int & finrule) {

		int i, j, k;

		for (j=0;j<nsaved;j++) {
			for (i=0;i<Group::ndims;i++) {
				dimflags[i] = saved[j].dimflags[i];
				num_dimpoints[i] = saved[j].num_dimpoints[i];
				for (k=0;k<8;k++) {
					dimpoints[i][k][0] = saved[j].dimpoints[i][k][0];
					dimpoints[i][k][1] = saved[j].dimpoints[i][k][1];
					dimpointerrors[i][k] = saved[j].dimpointerrors[i][k];
				}
			}
			nidims = saved[j].nidims;
			for (i=0;i<nidims;i++) {
				idims[i].Copy(saved[j].idims[nidims]);
			}
			if (Test(&(gpp[i]))) {
				finrule = i;
				return saved[j].category;
			}
		}
		return -1;

	};


	
	void TestSaved(Group *gpp, int ng) {

		int i, j, k, tp, tn, fp, fn, c;

		for (j=0;j<nsaved;j++) {
			c = saved[j].category;
			for (i=0;i<Group::ndims;i++) {
				dimflags[i] = saved[j].dimflags[i];
				num_dimpoints[i] = saved[j].num_dimpoints[i];
				for (k=0;k<8;k++) {
					dimpoints[i][k][0] = saved[j].dimpoints[i][k][0];
					dimpoints[i][k][1] = saved[j].dimpoints[i][k][1];
					dimpointerrors[i][k] = saved[j].dimpointerrors[i][k];
				}
			}
			nidims = saved[j].nidims;
			for (i=0;i<nidims;i++) {
				idims[i].Copy(saved[j].idims[nidims]);
			}
			tp = tn = fp = fn = 0;
			for (i=0;i<ng;i++) {
				if (gpp[i].category == c) {	
					if (Test(&(gpp[i]))) {
						tp++;
						if (obj_clusters[i]<(signed char)0) {
							obj_clusters[i] = (signed char)j;
						}
					} else {
						fn++;
					}
				} else {	
					if (Test(&(gpp[i]))) {
						fp++;
					} else {
						tn++;
					}
				}
			}
			printf("cluster %d true positive %d negative %d false positive %d negative %d\n", j, tp, tn, fp, fn);
		}

	};


	void Commit(int c) {
		int i, j;
		saved[nsaved].category = c;
		for (i=0;i<Group::ndims;i++) {
			saved[nsaved].dimflags[i] = dimflags[i];
			saved[nsaved].num_dimpoints[i] = num_dimpoints[i];
			for (j=0;j<8;j++) {
				saved[nsaved].dimpoints[i][j][0] = dimpoints[i][j][0];
				saved[nsaved].dimpoints[i][j][1] = dimpoints[i][j][1];
				saved[nsaved].dimpointerrors[i][j] = dimpointerrors[i][j];
			}
		}
		saved[nsaved].nidims = nidims;
		for (i=0;i<nidims;i++) {
			saved[nsaved].idims[nidims].Copy(idims[i]);
		}
		for (i=0;i<ng;i++) {
			if (removed[i]==0) {
				obj_clusters[SGAinds[i]] = (signed char)nsaved;
			}
		}
		for (i=0;i<naux;i++) {
			if (auxremoved[i]==0) {
				obj_clusters[auxSGAinds[i]] = (signed char)nsaved;
			}
		}
		nsaved++;
	};

	void SaveInterdim() {

		int i, n;
		int d1, d2;
		float mini, maxi;
		float minj, maxj;
		float error, bestvalue, diff;

		d1 = d2 = -1;
		for (i=0;i<Group::ndims;i++) {
			if (interdim[i]==1) {
				if (d1 < 0) {
					d1 = i;
				} else {
					d2 = i;
				}
			}
		}
		interdim[d1] = interdim[d2] = 0;

		// save the 2 dimensions involved
		idims[nidims].dim1 = d1;
		idims[nidims].dim2 = d2;
		// save the trend
		idims[nidims].trendflag = interdim_trend;
		// find the minimum and maximum points in dimension i
		mini = 2.0; maxi = -1.0;
		for (i=0;i<ng;i++) {
			if (removed[i]==0) {
				if (garray[i]->startingcoords[d1] < mini) {
					mini = garray[i]->startingcoords[d1];
					minj = garray[i]->startingcoords[d2];
				}
				if (garray[i]->startingcoords[d1] > maxi) {
					maxi = garray[i]->startingcoords[d1];
					maxj = garray[i]->startingcoords[d2];
				}
			}
		}
		for (i=0;i<naux;i++) {
			if (auxremoved[i]==0) {
				if (auxiliary[i]->startingcoords[d1] < mini) {
					mini = auxiliary[i]->startingcoords[d1];
					minj = auxiliary[i]->startingcoords[d2];
				}
				if (garray[i]->startingcoords[d1] > maxi) {
					maxi = auxiliary[i]->startingcoords[d1];
					maxj = auxiliary[i]->startingcoords[d2];
				}
			}
		}
		// store their (i,j) coordinates
		idims[nidims].first[0] = mini;
		idims[nidims].first[1] = minj;
		idims[nidims].last[0] = maxi;
		idims[nidims].last[1] = maxj;
		// find the error

		if (interdim_trend==0) { // linear trend
			error = 0.0;
			for (n=0;n<ng;n++) {
				if (removed[n]==0) {
					bestvalue = idims[nidims].first[1] + 
						((garray[n]->startingcoords[d1] - idims[nidims].first[0]) / (idims[nidims].last[0] - idims[nidims].first[0])) *
						(idims[nidims].last[1] - idims[nidims].first[1]);
					diff = garray[n]->startingcoords[d2] - bestvalue;
					if (diff < 0.0) diff = -diff;
					if (diff > error) error = diff;
				}
			}
			for (n=0;n<naux;n++) {
				if (auxremoved[n]==0) {
					bestvalue = idims[nidims].first[1] + 
						((auxiliary[n]->startingcoords[d1] - idims[nidims].first[0]) / (idims[nidims].last[0] - idims[nidims].first[0])) *
						(idims[nidims].last[1] - idims[nidims].first[1]);
					diff = auxiliary[n]->startingcoords[d2] - bestvalue;
					if (diff < 0.0) diff = -diff;
					if (diff > error) error = diff;
				}
			}
			idims[nidims].error = error * (idims[nidims].last[0] - idims[nidims].first[0]) /
											(idims[nidims].last[1] - idims[nidims].first[1]);
			if (idims[nidims].error < 0.0) idims[nidims].error = idims[nidims].error;
			idims[nidims].error = 2.0 * idims[nidims].error;
			printf("dimension %d %d trend, first %f %f last %f %f error %f\n", d1, d2, idims[nidims].first[0], idims[nidims].first[1],
				idims[nidims].last[0], idims[nidims].last[1], idims[nidims].error);
		} else if (interdim_trend==1) { // non-linear trend
			
			
		} else if (interdim_trend==2) { // constant

			idims[nidims].error = 0.0;
			for (n=0;n<ng;n++) {
				if (removed[n]==0) {
					error = garray[n]->startingcoords[d2] - idims[nidims].first[1];
					if (error < 0.0) error = -error;
					if (error > idims[nidims].error) idims[nidims].error = error;
				}
			}
			for (n=0;n<naux;n++) {
				if (auxremoved[n]==0) {
					error = auxiliary[n]->startingcoords[d2] - idims[nidims].first[1];
					if (error < 0.0) error = -error;
					if (error > idims[nidims].error) idims[nidims].error = error;
				}
			}
			idims[nidims].error = 2.0 * idims[nidims].error;
			printf("dimension %d %d constant, error %f\n", d1, d2, idims[nidims].error);
		}

		// increment number of interdims saved
		nidims++;
	};

};

#endif

