/************************************************************************
* editcm reads a file in cm format and edits outliers, add     		*
* uncertainties, and optionally applies a regression model to depth     *
* 		                                                        *
************************************************************************/
/************************************************************************
* Creator: David Sandwell                                               *
* Date   : March 26, 2008					        *
************************************************************************/
/************************************************************************
* Modification History: Seung Hee Kim					* 
* Date	:	June 27, 2008				                *
*                                                                       * 
************************************************************************/ 

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

float torben_median(float m[], int n) {

	int i, less, greater, equal;
	float min, max, guess, maxltguess, mingtguess;

	min = max = m[0];

	for (i=1;i<n;i++) {
		if (m[i]<min) min=m[i];
		if (m[i]>max) max=m[i];
	}

	while(1) {
		guess = (min+max)/2;
		less = 0;
		greater = 0;
		equal = 0;
		maxltguess = min;
		mingtguess = max;
		for (i=0;i<n;i++) {
			if (m[i]<guess) {
				less++;
				if (m[i]>maxltguess) maxltguess = m[i];
			} else if (m[i]>guess) {
				greater++;
				if (m[i]<mingtguess) mingtguess = m[i];
			} else equal++;
		}
		if (less <= (n+1)/2 && greater <= (n+1)/2) break;
		else if (less>greater) max = maxltguess;
		else min = mingtguess;
	}
	if (less >= (n+1)/2) return maxltguess;
	else if (less+equal >= (n+1)/2) return guess;
	else return mingtguess;
}

int main(int argc, char **argv) {
	
	#define NUM_REC (1000*10000)
	int i, imode, nrec, ndata, nflagged, sid;
	float *dd, *sigma;
	float d, MAD_torben, pflagged;
	
	struct cmRecType {
		int    time;
		double lon;
		double lat;
		int    depth;
		int    sig_h;
		int    sig_d;
		int    id;
		int    pred;
	} *jjPtr, *scan;
	
	scan = jjPtr = malloc(sizeof(struct cmRecType)*NUM_REC);
	dd = malloc(sizeof(float)*NUM_REC);
	
	if(argc < 2){
		fprintf(stderr,"usage: %s imode < file.cm > file.cm\n\n",argv[0]);
		fprintf(stderr,"          imode:\n");
		fprintf(stderr,"          0 - dump measured depth, pred depth, and absolute difference \n");
		fprintf(stderr,"          1 - calculate Median Absolute Deviation (MAD) - returns SID, MAD, Number of Records, Number of Unflagged Records \n");
		fprintf(stderr,"          2 - edit outliers and add depth uncertainty \n");
	exit(-1);
	}
	imode = atoi(argv[1]);
	if(imode < 0 || imode > 2) {
		printf(" imode out of range \n");
		exit(1);
	}
	
	/* read the data */
	
	nrec=0;
	scan = jjPtr;
	while (fscanf(stdin,"%d %lf %lf %d %d %d %d %d", 
	  &scan->time, &scan->lat, &scan->lon, &scan-> depth, &scan->sig_h, &scan->sig_d, &scan->id, &scan->pred)
		   != EOF) {
		
		nrec++;
		if(nrec > NUM_REC) {
			printf(" increase NUM_REC \n");
			exit(1);
		}
		scan++;
	}
	
	/* dump the measured and pred depth */
	
	scan = jjPtr;
	
	if (imode == 0) {
		for(i=0;i<nrec;i++) {
			if(scan->sig_d < 9999) 
				d=fabsf(scan->pred - scan->depth);
/*				d=fminf(d,fabsf(scan->pred));		*/
				printf("%d %d %f\n", scan->pred, scan->depth, d);
			scan++;
		}
		exit(1);
	}
	
	if (imode > 0) {
		
		ndata = 0;
		scan = jjPtr;
		for(i=0;i<nrec;i++) {
			sid = scan->id;
			if(scan->sig_d != 9999) {
				d=fabsf(scan->pred - scan->depth);
				dd[ndata] = d;
				ndata++;
			}
			scan++;
		}


/* Use qsort to determine median and compare to torben method */

		MAD_torben = torben_median(dd, ndata);

		if (imode == 1) {
			fprintf(stderr,"%d %f %d %d\n", sid, MAD_torben, nrec, ndata);
		}

	} 
	
	if(imode > 1){
		
		/* now flag the outliers where ... */

		nflagged = 0;
		scan = jjPtr;
		for(i=0;i<nrec;i++) {
			if(scan->sig_d != 9999) {
				d=fabsf(scan->pred - scan->depth);
				if( d > 8*MAD_torben && d >=200) {
					scan->sig_d = 9999;
					nflagged++;
				}
			}
			fprintf(stderr," %d %d %f \n",i,scan,d);
			scan++;
		}
		pflagged = (float)100 * ((float)nflagged / (float)ndata);
		fprintf(stderr,"%d %f %d %d %d %f\n", sid, MAD_torben, nrec, ndata, nflagged, pflagged);


		/* write all the data records */
		
		scan = jjPtr;
		for(i=0;i<nrec;i++) {
			printf("%d %.5f %.5f %d %d %d %d %d \n", scan->time, scan->lat, 
				   scan->lon, scan-> depth, scan->sig_h, scan->sig_d, scan->id, scan->pred);
			scan++;
		}
	}

	return 0;

}