/*--------------------------------------------------------------------
 *    $Id: minmax.c,v 1.21 2004/08/11 18:40:30 pwessel Exp $
 *
 *	Copyright (c) 1991-2004 by P. Wessel and W. H. F. Smith
 *	See COPYING file for copying and redistribution conditions.
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * minmax.c will read ascii or binary tables and report the
 * extreme values for all columns
 *
 * Author:	Paul Wessel
 * Date:	20-JUN-2000
 * Revised:	20-FEB-2001 BCH-J: Added -D option for dateline discontinuity
 * Version:	4.0
 */

#include "gmt.h"

int strip_blanks_and_output (double x, int col);

main (int argc, char **argv)
{
	int n, i, j, ncol, n_files = 0, fno, n_args, n_fields, n_delta = 0, n_expected_fields;
	int find_min = 0, find_max = 0, ecol = -1;

	BOOLEAN  error = FALSE, nofile = TRUE, done = FALSE, got_stuff = FALSE, first, step = FALSE, columns = FALSE, give_r_string = FALSE;
	BOOLEAN  dateline = FALSE, brackets = FALSE, zstep = FALSE, extrema = FALSE, work_on_abs_value, special = FALSE;

	double delta[BUFSIZ], *xyzmin, *xyzmax, west, east, south, north, low, high, value, e_min = DBL_MAX, e_max = -DBL_MAX, *in;

	char line[BUFSIZ], file[BUFSIZ], chosen[BUFSIZ];

	FILE *fp = NULL;
	
	argc = GMT_begin (argc, argv);
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				/* Common parameters */
                      
				case 'H':
				case ':':
  				case 'b':
  				case 'f':
				case '\0':
					error += GMT_get_common_args (argv[i], 0, 0, 0, 0);
					break;

				/* Supplemental parameters */

				case 'C':
					columns = TRUE;
					break;
                                case 'D':
                                        dateline = TRUE;
                                        break;
                                case 'E':
                                        extrema = TRUE;
					switch (argv[i][2]) {
						case 'L':
							find_min = 2;
							break;
						case 'l':
							find_min = 1;
							break;
						case 'H':
							find_max = 2;
							break;
						case 'h':
							find_max = 1;
							break;
						default:
							error ++;
							fprintf (stderr, "%s: GMT SYNTAX ERROR -E. Flags are L|l|H|h\n", GMT_program);
							break;
					}
					ecol = atoi (&argv[i][3]);
                                        break;
				case 'M':
					GMT_multisegment (&argv[i][2]);
					break;
				case 'I':
					if (argv[i][2] == 'p') special = TRUE;
					j = (special) ? 3 : 2;
					n_delta = GMT_getincn (&argv[i][j], delta, BUFSIZ);
					step = TRUE;
					break;
				case 'L':	/* Obsolete, but backward compatibility prevails [use -f instead] */
					GMT_io.in_col_type[0] = GMT_io.out_col_type[0] = GMT_IS_LON;
					GMT_io.in_col_type[1] = GMT_io.out_col_type[1] = GMT_IS_LAT;
					fprintf (stderr, "%s: Option -L is obsolete (but is processed correctly).  Please use -fg instead\n", GMT_program);
					break;
				case 'T':
					delta[0] = atof (&argv[i][2]);
					zstep = TRUE;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}
	
	if (error || GMT_quick) {	/* Because it's ok to give no arguments */
		fprintf (stderr, "minmax %s - Find extreme values in ASCII tables\n\n", GMT_VERSION);
		fprintf (stderr, "usage: minmax [files] [-C] [-D] [-E<L|l|H|h><col>] [-H[<nrec>]] [-Idx[/dy[/<dz>..]] [-M[<flag>]]\n");
		fprintf (stderr, "\t[-T<dz>] [-:] [-bi[s][<n>] [-f[i|o]<colinfo>]\n");
             
              	if (GMT_quick) exit (EXIT_FAILURE);
 
		fprintf (stderr, "\t-C formats the min and max into separate columns\n");
		fprintf (stderr, "\t-D sets longitude discontinuity to dateline (requires -fi0x,1y)\n");
		fprintf (stderr, "\t-E Return the record with extreme value in specified column <col>\n");
		fprintf (stderr, "\t   Specify l or h for min or max value, respectively.  Upper case L or H\n");
		fprintf (stderr, "\t   means we operate instead on the absolute values of the data.\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I returns textstring -Rw/e/s/n to nearest multiple of dx/dy (assumes 2+ col data)\n");
		fprintf (stderr, "\t   If -C is set then no -R string is issued.  Instead, the number of increments\n");
		fprintf (stderr, "\t   given determines how many columns are rounded off to the nearest multiple.\n");
		GMT_explain_option ('M');
		fprintf (stderr, "\t-T returns textstring -Tzmin/zmax/dz to nearest multiple of the given dz\n");
		fprintf (stderr, "\t   Calculations are based on the first column only\n");
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 2 input columns\n");
		GMT_explain_option ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	if (step && !special && n_delta == 1) {		/* Special case of dy = dx if not given */
		delta[1] = delta[0];
		n_delta = 2;
	}
	if (step && !columns && n_delta < 2) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  -Ip requires -C\n", GMT_program);
		error++;
	}
	if (step && zstep) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Only one of -I and -T can be specified\n", GMT_program);
		error++;
	}
	if (zstep && delta[0] <= 0.0 ) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -T option.  Must specify a positive increment\n", GMT_program);
		error++;
	}
	if (step) {
		for (i = 0; i < n_delta; i++) {
			if (delta[i] <= 0.0) {
				fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", GMT_program);
				error++;
			}
		}
	}
	if (GMT_io.binary[GMT_IN] && gmtdefs.io_header[GMT_IN]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] == 0) GMT_io.ncol[GMT_IN] = 2;
	if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] < 1) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data (-bi) must have at least 1 column\n", GMT_program);
		error++;
	}
	if (GMT_io.in_col_type[0] != GMT_IS_LON && dateline) {
	  	fprintf (stderr, "%s: GMT SYNTAX ERROR -D option: requires -fi0x,1y\n", GMT_program);
		error++;
	}
	if (extrema && ecol < 0) {
	  	fprintf (stderr, "%s: GMT SYNTAX ERROR -E option: requires a positive column\n", GMT_program);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands4 */

	if (GMT_io.binary[GMT_IN] && gmtdefs.verbose) {
		char *type[2] = {"double", "single"};
		fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[GMT_IN], type[GMT_io.single_precision[GMT_IN]]);
	}

	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	
	n_args = (argc > 1) ? argc : 2;
	west = south = DBL_MAX;	east = north = -DBL_MAX;
	
	xyzmin = (double *) GMT_memory (VNULL, (size_t)1, sizeof (double), GMT_program);
	xyzmax = (double *) GMT_memory (VNULL, (size_t)1, sizeof (double), GMT_program);
	
	give_r_string = (step && !columns);
	brackets = !columns;
	work_on_abs_value = (extrema && (find_min == 2 || find_max == 2));

	n_expected_fields = (GMT_io.ncol[GMT_IN]) ? GMT_io.ncol[GMT_IN] : BUFSIZ;

	for (fno = 1; !done && fno < n_args; fno++) {     /* Loop over input files, if any */
		if (!nofile && argv[fno][0] == '-') continue;

		
		if (nofile) {   /* Just read standard input */
			fp = GMT_stdin;
			done = TRUE;
			strcpy (file, "<stdin>");
			if (gmtdefs.verbose) fprintf (stderr, "%s: Reading from standard input\n", GMT_program);
#ifdef SET_IO_MODE
			GMT_setmode (GMT_IN);
#endif
		}
		else {
			strcpy (file, argv[fno]);
			if ((fp = GMT_fopen (file, GMT_io.r_mode)) == NULL) {
				fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, file);
				continue;
			}
		}
		
		if (gmtdefs.io_header[GMT_IN]) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (line, BUFSIZ, fp);
		
		n = ncol = 0;
		n_expected_fields = (GMT_io.ncol[GMT_IN]) ? GMT_io.ncol[GMT_IN] : BUFSIZ;
		first = TRUE;

		n_fields = GMT_input (fp, &n_expected_fields, &in);

		while (! (GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

			while ((GMT_io.status & GMT_IO_SEGMENT_HEADER) && !(GMT_io.status & GMT_IO_EOF)) {
					n_fields = GMT_input (fp, &n_expected_fields, &in);
			}
			if ((GMT_io.status & GMT_IO_EOF)) continue;	/* At EOF */

			if (GMT_io.status & GMT_IO_MISMATCH) {
				fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected_fields, n);
				exit (EXIT_FAILURE);
			}

			if (first) {	/* First time, allocate # of columns */
			
				ncol = n_expected_fields;
				if (extrema && ecol >= ncol) {
	  				fprintf (stderr, "%s: GMT SYNTAX ERROR -E option: Chosen column exceeds column range\n", GMT_program);
					exit (EXIT_FAILURE);
				}

				/* Now we know # of columns, so allocate memory */
			
				xyzmin = (double *) GMT_memory ((void *)xyzmin, (size_t)ncol, sizeof (double), GMT_program);
				xyzmax = (double *) GMT_memory ((void *)xyzmax, (size_t)ncol, sizeof (double), GMT_program);
			
				for (i = 0; i < ncol; i++) {	/* Initialize */
					xyzmin[i] = +DBL_MAX;
					xyzmax[i] = -DBL_MAX;
				}
				if (step && ncol < 2 && !columns) step = FALSE;
				first = FALSE;
			}
		
			/* Decode all fields and update minmax arrays */

			while (dateline && in[0] > 180.0) in[0] -= 360.0;

			if (extrema && !GMT_is_dnan (in[ecol])) {
				value = (work_on_abs_value) ? fabs (in[ecol]) : in[ecol];
				if (find_min && value < e_min) {
					e_min = value;
					strcpy (chosen, GMT_io.current_record);
				}
				else if (find_max && value > e_max) {
					e_max = value;
					strcpy (chosen, GMT_io.current_record);
				}
			}
			else if (!extrema) {
				for (i = 0; i < ncol; i++) {
					if (GMT_is_dnan (in[i])) continue;
					if (in[i] < xyzmin[i]) xyzmin[i] = in[i];
					if (in[i] > xyzmax[i]) xyzmax[i] = in[i];
				}
			}
			
			n++;

			n_fields = GMT_input (fp, &n_expected_fields, &in);
		}
		if (fp != GMT_stdin) GMT_fclose (fp);
		
		if (!got_stuff) got_stuff = (n > 0);	/* We were able to open and read at least 1 record from a file */

		if (give_r_string && got_stuff) {
			west  = MIN (west, xyzmin[0]);		east  = MAX (east, xyzmax[0]);
			south = MIN (south, xyzmin[1]);		north = MAX (north, xyzmax[1]);
		}
		else if (zstep && got_stuff) {
			west  = MIN (west, xyzmin[0]);		east  = MAX (east, xyzmax[0]);
		}
		else if (!extrema && n > 0) {
			if (!columns) fprintf (GMT_stdout, "%s: N = %d\t", file, n);
			for (i = 0; i < ncol; i++) {
				if (xyzmin[i] == DBL_MAX) {
					low = high = GMT_d_NaN;
				}
				else if (i < n_delta) {	/* Special treatment for x and y if selected */
					low  = (step) ? floor (xyzmin[i] / delta[i]) * delta[i] : xyzmin[i];
					high = (step) ? ceil (xyzmax[i] / delta[i]) * delta[i] : xyzmax[i];
				}
				else {
					low = xyzmin[i];
					high = xyzmax[i];
				}
				if (brackets) fputc ('<', GMT_stdout);
				GMT_ascii_output_one (GMT_stdout, low, i);
				(columns) ? fputc ('\t', GMT_stdout) : fputc ('/', GMT_stdout);
				GMT_ascii_output_one (GMT_stdout, high, i);
				if (brackets) fputc ('>', GMT_stdout);
				if (i < (ncol - 1)) fputc ('\t', GMT_stdout);
			}
			fprintf (GMT_stdout, "\n");
		}
	}
	if (got_stuff) {
		if (give_r_string) {
			west  = floor (west / delta[0]) * delta[0];	east  = ceil (east / delta[0]) * delta[0];
			south = floor (south / delta[1]) * delta[1];	north = ceil (north / delta[1]) * delta[1];
			if (east < west) east += 360.0;
			fprintf (GMT_stdout, "-R");
			strip_blanks_and_output (west, 0);	fputc ('/', GMT_stdout);
			strip_blanks_and_output (east, 0);	fputc ('/', GMT_stdout);
			strip_blanks_and_output (south, 1);	fputc ('/', GMT_stdout);
			strip_blanks_and_output (north, 1);	fputc ('\n', GMT_stdout);
		}
		else if (zstep) {	/* -T option */
			west  = floor (west / delta[0]) * delta[0];	east  = ceil (east / delta[0]) * delta[0];
			fprintf (GMT_stdout, "-T");
			strip_blanks_and_output (west, 0);	fputc ('/', GMT_stdout);
			strip_blanks_and_output (east, 0);	fputc ('/', GMT_stdout);
			strip_blanks_and_output (delta[0], 0);	fputc ('\n', GMT_stdout);
		}
		else if (extrema)
			fprintf (GMT_stdout, "%s", chosen);
	}
	else if (!got_stuff)
		fprintf (stderr, "%s: No input data found!\n", GMT_program);
	
	GMT_free ((void *)xyzmin);
	GMT_free ((void *)xyzmax);
	
	GMT_end (argc, argv);
}

int strip_blanks_and_output (double x, int col)
{
	/* Alternative to GMT_ascii_output_one that strips off leading blanks first */
	
	int k;
	char text[64];
	
	GMT_ascii_format_one (text, x, GMT_io.out_col_type[col]);
	for (k = 0; text[k] && text[k] == ' '; k++);
	return (fprintf (GMT_stdout, "%s", &text[k]));
}
