/*******************************************************************

    Program:  minmax -- started by Tom Solomon  3/11/93

web pointer to this program:
   http://www.physics.emory.edu/~weeks/software/maxmin.html
limited support for this program available via email:
   weeks@physics.emory.edu

This program is public domain, but please leave the email &
web information above intact.  This program originally
written by Tom Solomon, with features added by Eric Weeks.
Please contact Eric if you have questions, don't bug Tom.
Note this program has several undocumented features, which
you can perhaps figure out by staring at the source code.
I make no guarantees as to the usefulness of this code, but
I believe it to be bug-free.



9/06/94:	(Eric Weeks) added windowing, -n -w -d -s options
9/10/94:  Try to add in intelligent double-windowing
9/11/94:  rewrite double-windowing; read in all data before starting; add in
		the scanit function to check against two successive minima (or maxima)
9/11/94:	add in filtering (-f option)
9/12/94:	add in interpolation (-i option)
		Note that currently, program is fully backwards compatible.  ie:
		minmax			: works exactly like Tom's original version
		minmax -w 10		: works with a single window
		minmax -w 10 -W 20	: works with double-window scheme
		minmax -nif 4		: turn on some new options
6/22/97:  clean up a little bit for web distribution
5/18/15:  fix a bug found by Minwoong Joe; remove -i option



 *******************************************************************/
#include <stdio.h>
#include <math.h>

double tdat[100000], ddat[100000];
int info[100000];

/* ---------------------------------------------------------------- */
void output(int flag,double t,double datum)  {
   if (flag == 1) { printf("%g %g\n",t,datum); } else { printf("%g\n",datum); }
}

/* ---------------------------------------------------------------- */
void output2(int value,int minfl,int maxfl,int debug,int pairfl,int i, int w) {

	/* value=0 for nothing, 1 for min, 2 for max, 3 for both */

	if (debug==1)  printf("w%1d",w);
	/* print if it's a minimum */
	if (((value==1)||(value==3)) && (minfl==1))  {
		if (debug==1)  printf("min ");
		output(pairfl,tdat[i],ddat[i]);
	}
	/* print if a max, & we haven't printed it out */
	if ((value>1)&&(maxfl==1))  {
		if ((value==2)||(minfl==0)) {
				if (debug==1)  printf("MAX ");
				output(pairfl,tdat[i],ddat[i]);
		}
	}
}

/* ---------------------------------------------------------------- */
double addnoise(double d, double *yy)  {
	double xx,mult;
	xx = (*yy);
	xx = 4.0*xx*(1.0-xx);
	while ((xx<.01)||(xx>.99))  xx = 4.0*xx*(1.0-xx);
	mult = .995 + .010*xx;
	if ((d>.2)||(d<(-.2)))  { d*= mult; } else { d+= 0.0001*xx; }
	(*yy) = xx;
	return d;
}

/* ---------------------------------------------------------------- */
/* function: check												*
 *	This function scans the specified window, to see if the middle data	*
 *	point is a maximum or a minimum within the window.				*/
int check(int window,int middle,int strict)  {
	int minmin,maxmax,value,j;
	double x;

	minmin=1;maxmax=1;
	if (strict==1)  {
		for (j=middle-window+1;((j<middle+window)&&((maxmax+minmin)>0));j++) {
			if (j != middle)  {
				if (ddat[middle]>=ddat[j])  minmin=0;
				if (ddat[middle]<=ddat[j])  maxmax=0;
			}
		}
	} else {
		for (j=middle-window+1;((j<middle+window)&&((maxmax+minmin)>0));j++) {
			if (ddat[middle]>ddat[j])  minmin=0;
			if (ddat[middle]<ddat[j])  maxmax=0;
		}
	}
	value= minmin+2*maxmax;
	return value;
}

/* ---------------------------------------------------------------- */
double gettol(int window,int middle)  {
	double x,y,z;

	x = ddat[middle] - ddat[middle-window+1]; if (x<0)  x=-x;
	y = ddat[middle] - ddat[middle+window-1]; if (y<0)  y=-y;
	z = (x>y) ? y : x; return z;
	/* returns absolute value of WORST case */
}
/* ---------------------------------------------------------------- */

/* Function scanit:												*
 *   two successive extrema have been found that are both the same (either	*
 *	maxima or minima.)  So, we know that in between these two we MUST be	*
 *	able to find another extrema (of the opposite type).  So...  find it! */
int scanit(int value, int prevptr, int numptr)  {
	int t,ptr;
	double max,min;

	if (value==1)  {
		max=ddat[prevptr];
		for (t=prevptr+1;t<numptr;t++)
			if (max<ddat[t])  { max=ddat[t]; ptr=t; }
	} else {
		min=ddat[prevptr];
		for (t=prevptr+1;t<numptr;t++)
			if (min>ddat[t])  { min=ddat[t]; ptr=t; }
	}
	return ptr;
}
/* ---------------------------------------------------------------------- */

main(argc,argv)
int argc;
char **argv;
{
	int c,window1,minfl,maxfl,atoi(),ptr,r;
	int window2,twowin,previous,prevptr;
	int pairfl,increase,endflag,strict,flag,minmin,maxmax,i,j;
	int noise,debug,value,num,numptr,filter;
	double datum,olddatum,lastpeak,t,oldt,mult,xx,tol,thetol,ave;
	FILE *fp, *fopen();
	extern int optind;
	extern char *optarg;

	window1 = 0;  window2 = 0; twowin=0;
	minfl = 1; maxfl = 1; pairfl = 0; endflag = 0;
	strict = 1;  noise = 0;  debug = 0; tol = 0.02; filter = 0;
	/******************************************************************/
	/************ PROCESS COMMAND LINE ARGUMENTS **********************/
	/******************************************************************/
	while ((c = getopt(argc, argv, "f:dnehmMpsW:w:t:i")) != EOF)
        switch (c) {
case 'h': fprintf(stderr,"This program looks for local peaks and troughs of a time series.\n");
fprintf(stderr,"Usage:  minmax [options] [<] (filename)");
fprintf(stderr,"\n\nOptions:\n\n");
fprintf(stderr," -h  :  this help message\n");
fprintf(stderr," -w #:  small window size (number of points) [default: no window]\n");
fprintf(stderr," -m  :  Look for maxima only  [default: print both]\n");
fprintf(stderr," -M  :  Look for minima only  [default: print both]\n");
fprintf(stderr," -p  :  Data in x-y pairs. (look for extrema of y)\n");
fprintf(stderr," -e  :  Print out first and last point.\n");
fprintf(stderr," -s  :  strictness off -- use >= for comparison [ > ]\n");
/* UNDOCUMENTED FEATURES ------------------------------
fprintf(stderr," -W #:  large window size (number of points)\n");
fprintf(stderr,"        (if small window is specified but not large, then\n");
fprintf(stderr,"        the program uses only single-windowing)\n");
fprintf(stderr," -n  :  add a small amount of noise (1%%)\n");
fprintf(stderr," -d  :  debug on\n");
fprintf(stderr," -t #:  tolerance for double-windowing [%.2f]\n",tol);
fprintf(stderr,"        (if maximum-minimum<tol within the small window,\n");
fprintf(stderr,"        check the large window.)\n");
fprintf(stderr," -f #:  turn on boxcar filtering over # boxes [default:no]\n");
*/
          exit(1);
          break;
case 'w': window1 = atoi(optarg);
		if (window1<2) {
			fprintf(stderr,"error! window must be bigger than 1\n");
			exit(1);
		}
		break;
case 'W': window2 = atoi(optarg);
		if (window2<2) {
			fprintf(stderr,"error! window must be bigger than 1\n");
			exit(1);
		}
		break;
case 't': tol=atof(optarg);
		if (tol<=0)  {
			fprintf(stderr,"error! tolerance must be positive.\n");
			exit(1);
		}
		/* fprintf(stderr,"tolerance set ---->%f\n",tol); */
		break;
case 'f': filter=atoi(optarg);
		break;
case 'm': minfl = 0;
		break;
case 'M': maxfl = 0;
		break;
case 'p': pairfl = 1;
		break;
case 'e': endflag = 1;
		break;
case 's': strict = 0;
		break;
case 'n': noise = 1;
		break;
case 'd': debug = 1;
		break;
        }

	if ((window1>0)&&(window2<=window1)&&(window2>0))  {
		fprintf(stderr,"large window isn't larger than small window\n");
		fprintf(stderr,"-->setting large window to twice small window\n");
		window2=2*window1;
	}
	if ((window1>0)&&(window2>0))  twowin=1;
	if ((window1==0)&&(window2>0))  {
		fprintf(stderr,"error! large window specified but not small.\n");
		fprintf(stderr,"I'm assuming you want to use one window only.\n");
		window1=window2;
	}
	if ((window1>0)&&(window2==0))  window2=window1;
	/*********************** DONE WITH COMMAND LINE ARGUMENTS ************/

	/******** INPUT **********/
	argc -= (optind-1) ; argv += (optind-1) ;
	fp = (argc > 1) ? fopen(*++argv, "r") : stdin;

	num=0;
	if (pairfl==1)  {
		oldt=99999.9;
		while (fscanf(fp,"%lf %lf",&t,&datum)==2)  {
			num++; tdat[num]=t; ddat[num]=datum;
/* printf("%f %f\n",tdat[num],ddat[num]); */
			oldt=t; olddatum=datum;
		}
	} else {
		while (fscanf(fp,"%lf",&datum)==1)  {
			num++; tdat[num]=0.0; ddat[num]=datum;
		}
	}
	/****** DONE WITH INPUT ******/
	filter--;
	if (filter>0)  {
		ddat[0]=0.0;
		mult = 1.0/((double)(filter+1.0));
		for (i=1;i<=num-filter;i++)  {
			tdat[i]=(tdat[i]+tdat[i+filter])*0.5;
			ave = 0.0;
			for (j=i;j<=i+filter;j++)  ave += ddat[j];
			ddat[i]=ave*mult;
		}
		num -= filter;
	}
	/********* DONE WITH FILTERING *******/


	
	/***** print out first point, if we want to *****/
	if (endflag == 1) output(pairfl,tdat[1],ddat[1]);

	t=tdat[1]; datum=ddat[1];
	olddatum = datum; oldt = t;
	/***** GET FIRST LINE *****/
	t=tdat[2]; numptr=2;

	if (window1==0)  {
		/***** ARE WE GOING UP OR DOWN??  *****/
		t=tdat[numptr];  datum=ddat[numptr];
		if (olddatum < datum) { increase = 1; } else { increase = 0; }
		while (numptr<=num)  {
			t=tdat[numptr];  datum=ddat[numptr];
			if ((datum<olddatum)&&(increase==1)) {
				if (maxfl == 1)  output(pairfl,oldt,olddatum);
				increase = 0;
			}
			if ((datum>olddatum)&&(increase==0)) {
				if (minfl == 1)  output(pairfl,oldt,olddatum);
				increase = 1;
			}
			olddatum=datum; oldt = t;
			numptr++;
		}
	} else {
		/*********************** WINDOWING *******************/
		if (noise==1)  {
			xx = .01234;
			for (i=1;i<=num;i++)  ddat[i] = addnoise(ddat[i],&xx);
		}
		info[1] = 0; info[num] = 0;
		for (i=2;i<num;i++)  {
			minmin=1; maxmax=1;
			if ( (ddat[i]>ddat[i-1]) || (ddat[i]>ddat[i+1]) )  minmin=0;
			if ( (ddat[i]<ddat[i-1]) || (ddat[i]<ddat[i+1]) )  maxmax=0;
			info[i] = minmin+2*maxmax;
		}
		previous=-1;
		for (numptr=window2;numptr<=(num-window2);numptr++)  {
			value = 0;
			if (info[numptr]>0)  value=check(window1,numptr,strict);
			if (value>0)  {
				if ((twowin==0) || (gettol(window1,numptr)>tol))  {
					if (previous==value)  {
						ptr=scanit(value,prevptr,numptr);
						output2(3-value,minfl,maxfl,debug,pairfl,ptr,4);
					}
					output2(value,minfl,maxfl,debug,pairfl,numptr,1);
					previous=value;  prevptr=numptr;
				} else {     /* DOUBLE WINDOWING */
					value=check(window2,numptr,strict);
					if (previous==value)  {
						ptr=scanit(value,prevptr,numptr);
						output2(3-value,minfl,maxfl,debug,pairfl,ptr,4);
					}
					if (value>0)  {
						output2(value,minfl,maxfl,debug,pairfl,numptr,2);
						previous=value; prevptr=numptr;
					}
				}	
			}
		} /* end of for loop */
		/******************** END OF WINDOWING ***************/
	}

	/***** PRINT OUT LAST POINT *****/
	if (endflag == 1)  output(pairfl,tdat[num],ddat[num]);
}