/*--------------------------------------------------------------------
 *    $Id: blockmode.c,v 1.36 2006/05/14 01:43:04 pwessel Exp $
 *
 *	Copyright (c) 1991-2006 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
 *--------------------------------------------------------------------*/

/*
   blockmode.c
   Takes lon, lat, data, [weight] on GMT_stdin or file and writes out one value
   per cell, where cellular region is bounded by West East South North and
   cell dimensions are delta_x, delta_y.
      
   Author: 	Walter H. F. Smith
   Date:	28 October, 1998
   Version:	1st version based on GMT3.1b blockmedian.c, but with changes
   		to struct BLK_DATA and to sort comparison routines to facilitate
   		mode coding.
   Modified:	3.3.5: 10-JUL-2000 PW: Added -L option
   Version:	3.4 01-MAR-2001 by PW, Replace -N with -F, and added -C
   Version:	4 01-MAR-2003 by PW
   Version:	4.1: 14-SEP-2005 by PW, Added enhanced -I
   		4-APR-2006 by PW: Added -E for LMS scale, low, and high value
			Also implemented size_t counters to be 64-bit compatible.
*/

#define BLOCKMODE

#include "gmt.h"
#include "block_subs.h"

int main (int argc, char **argv)
{
	BOOLEAN	error = FALSE, nofile = TRUE, done = FALSE, first = TRUE, mode_xy;

	/* Default value for go_quickly = FALSE for backward compatibility with 3.0  */

	FILE *fp = NULL;

	double	*in, out[4], wesn[4], i_n_in_cell, weight, *z_tmp = NULL;

	int	i, j, ix, iy, fno, n_files = 0, n_args, n_req;
	int	n_expected_fields, n_fields, n_out, w_col;
	
	size_t index, first_in_cell, first_in_new_cell, n_lost, n_read;
	size_t n_pitched, n_cells_filled, n_alloc, n_in_cell, nz, nz_alloc;

	char	modifier, buffer[BUFSIZ], format[BUFSIZ];

	struct GRD_HEADER h;

	struct BLK_DATA *data;
	struct BLOCKMODE_CTRL *Ctrl;
	
	double weighted_mode (struct BLK_DATA *d, double wsum, size_t n, int k);
	
	argc = GMT_begin (argc, argv);

	Ctrl = (struct BLOCKMODE_CTRL *) New_Blockmode_Ctrl ();	/* Allocate and initialize a new control structure */

	GMT_grd_init (&h, argc, argv, FALSE);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
              
				/* Common parameters */
                      
				case 'H':
				case 'R':
				case 'V':
				case ':':
				case 'b':
				case 'f':
				case '\0':
                                      error += GMT_parse_common_options (argv[i], &h.x_min, &h.x_max, &h.y_min, &h.y_max);
                                      break;
                              
				/* Supplemental parameters */
                              
				case 'C':
					Ctrl->C.active = TRUE;
					break;
				case 'E':
					Ctrl->E.active = TRUE;		/* Extended report with standard deviation, min, and max in cols 4-6 */
					break;
				case 'I':
					Ctrl->I.active = TRUE;
					if (GMT_getinc (&argv[i][2], &Ctrl->I.xinc, &Ctrl->I.yinc)) {
						GMT_inc_syntax ('I', 1);
						error = 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 -f instead\n", GMT_program);
					break;
				case 'N':	/* Backward compatible with 3.3.6 */
				case 'F':
					Ctrl->F.active = TRUE;
					break;
				case 'Q':
					Ctrl->Q.active = TRUE;
					break;
				case 'W':
					Ctrl->W.active = TRUE;
					if ( (modifier = argv[i][2]) == 'i' || modifier == 'I')
						Ctrl->W.weighted[GMT_IN] = TRUE;
					else if (modifier == 'O' || modifier == 'o')
						Ctrl->W.weighted[GMT_OUT] = TRUE;
					else
						Ctrl->W.weighted[GMT_IN] = Ctrl->W.weighted[GMT_OUT] = TRUE;
					break;

				default:
					error = TRUE;
                                        GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}

	if (argc == 1 || GMT_give_synopsis_and_exit) {
		fprintf (stderr, "blockmode %s - Block averaging by mode estimates\n\n", GMT_VERSION);
		fprintf (stderr, "usage: blockmode [infile(s)] %s %s\n", GMT_I_OPT, GMT_Rgeo_OPT);
		fprintf (stderr, "\t[-C] [-E] [-F] [%s] [-Q] [-V] [-W[i][o] ] [%s] [%s]\n", GMT_H_OPT, GMT_t_OPT, GMT_b_OPT);
		fprintf (stderr, "\t[%s]\n\n", GMT_f_OPT);

		if (GMT_give_synopsis_and_exit) exit (EXIT_FAILURE);

		GMT_inc_syntax ('I', 0);
		GMT_explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-C Output center of block and mode z-value  [Default is mode location (but see -Q)]\n");
		fprintf (stderr, "\t-E Extend output with LMS scale (s), low (l), and high (h) value per block, i,e,\n");
		fprintf (stderr, "\t   output (x,y,z,s,l,h[,w]) [Default outputs (x,y,z[,w]); see -W regarding w.\n");
		fprintf (stderr, "\t-F Offsets registration so block edges are on gridlines (pixel reg.).  [Default: grid reg.]\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-Q Quicker; get mode z and mean x,y.  [Default gets mode x, mode y, mode z.]\n");
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W sets Weight options.\n");
		fprintf (stderr, "\t-  -Wi reads Weighted Input (4 cols: x,y,z,w) but writes only (x,y,z[,s,l,h]) Output.\n");
		fprintf (stderr, "\t   -Wo reads unWeighted Input (3 cols: x,y,z) but reports sum (x,y,z[,s,l,h],w) Output.\n");
		fprintf (stderr, "\t   -W with no modifier has both weighted Input and Output; Default is no weights used.\n");
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 3 [6 with -E] columns (or 4 [7] if -W is set).\n");
		GMT_explain_option ('o');
		GMT_explain_option ('n');
		GMT_explain_option ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}

	if (Ctrl->C.active && Ctrl->Q.active) {
		fprintf (stderr, "%s: GMT WARNING:  -C overrides -Q\n", GMT_program);
		Ctrl->Q.active = FALSE;
	}
	if (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
		error++;
	}
	if (Ctrl->I.xinc <= 0.0 || Ctrl->I.yinc <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", GMT_program);
		error = TRUE;
	}
	if (GMT_io.binary[GMT_IN] && GMT_io.io_header[GMT_IN]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
	n_req = (Ctrl->W.weighted[GMT_IN]) ? 4 : 3;
	if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] == 0) GMT_io.ncol[GMT_IN] = n_req;
	if (GMT_io.binary[GMT_IN] && n_req > GMT_io.ncol[GMT_IN]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  binary input data must have at least %d columns\n", GMT_program, n_req);
		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]]);
	}

#ifdef SET_IO_MODE
	GMT_setmode (GMT_OUT);
#endif
	
	mode_xy = !Ctrl->C.active;
	h.x_inc = Ctrl->I.xinc;
	h.y_inc = Ctrl->I.yinc;
	h.node_offset = Ctrl->F.active;
	GMT_RI_prepare (&h);	/* Ensure -R -I consistency and set nx, ny */
	h.xy_off = 0.5 * h.node_offset;	/* Use to calculate mean location of block */

	if (gmtdefs.verbose) {
		sprintf (format, "%%s: W: %s E: %s S: %s N: %s nx: %%d ny: %%d\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, format, GMT_program, h.x_min, h.x_max, h.y_min, h.y_max, h.nx, h.ny);
	}

	n_read = n_pitched = 0;
        n_alloc = GMT_CHUNK;
        data = (struct BLK_DATA *) GMT_memory (VNULL, n_alloc, sizeof(struct BLK_DATA), GMT_program);

	GMT_set_xy_domain (wesn, &h);	/* May include some padding if gridline-registered */
		
	/* Read the input data  */

	n_expected_fields = (GMT_io.binary[GMT_IN]) ? GMT_io.ncol[GMT_IN] : 3 + Ctrl->W.weighted[GMT_IN];

	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	n_args = (argc > 1) ? argc : 2;

	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;
#ifdef SET_IO_MODE
			GMT_setmode (GMT_IN);
#endif
		}
		else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) {
			fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]);
			continue;
		}

		if (!nofile && gmtdefs.verbose) fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[fno]);

		if (GMT_io.io_header[GMT_IN]) {
			for (i = 0; i < GMT_io.n_header_recs; i++) {
				GMT_fgets (buffer, BUFSIZ, fp);
				buffer[strlen(buffer)-1] = 0;
				if (first && GMT_io.io_header[GMT_OUT]) (Ctrl->W.weighted[GMT_OUT] && !(Ctrl->W.weighted[GMT_IN])) ? fprintf (GMT_stdout, "%s weights\n", buffer) : fprintf (GMT_stdout, "%s\n", buffer);
			}
			first = FALSE;
		}

		while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

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

			if (GMT_is_dnan (in[2])) continue;	/* Skip when z = NaN */

			n_read++;

			if (GMT_y_is_outside (in[1],  wesn[2], wesn[3])) continue;	/* Outside y-range */
			if (GMT_x_is_outside (&in[0], wesn[0], wesn[1])) continue;	/* Outside x-range */

			ix = GMT_x_to_i (in[0], h.x_min, h.x_inc, h.xy_off, h.nx);
			if ( ix < 0 || ix >= h.nx ) continue;
			iy = GMT_y_to_j (in[1], h.y_min, h.y_inc, h.xy_off, h.ny);
			if ( iy < 0 || iy >= h.ny ) continue;

			index = ((size_t)iy) * ((size_t)h.nx) + ((size_t)ix);		/* 64-bit safe 1-D index */

			data[n_pitched].i = index;
			if (mode_xy) {
				data[n_pitched].a[BLK_X] = in[0];
				data[n_pitched].a[BLK_Y] = in[1];
			}
			data[n_pitched].a[BLK_Z] = in[2];
			data[n_pitched].a[BLK_W] = (Ctrl->W.weighted[GMT_IN]) ? in[3] : 1.0;

			n_pitched++;
			if (n_pitched == n_alloc) {
				n_alloc += GMT_CHUNK;
				data = (struct BLK_DATA *) GMT_memory ((void *)data, n_alloc, sizeof(struct BLK_DATA), GMT_program);
			}
		}
		if (fp != GMT_stdin) GMT_fclose(fp);

	}

	data = (struct BLK_DATA *) GMT_memory ((void *)data, n_pitched, sizeof(struct BLK_DATA), GMT_program);
	n_lost = n_read - n_pitched;
	if (gmtdefs.verbose) fprintf(stderr,"%s: N read: %u N used: %u N outside_area: %u\n", GMT_program, n_read, n_pitched, n_lost);

	/* Ready to go. */

	n_out = (Ctrl->W.weighted[GMT_OUT]) ? 4 : 3;
	if (Ctrl->E.active) n_out += 3;
	w_col = n_out - 1;	/* Weights always reported in last output column */
	nz_alloc = GMT_CHUNK;
	if (Ctrl->E.active) z_tmp = (double *) GMT_memory (VNULL, nz_alloc, sizeof (double), GMT_program);	/* tmp array for L1 scale calculation */

	/* Sort on index and Z value */

	qsort((void *)data, n_pitched, sizeof (struct BLK_DATA), BLK_compare_index_z);

	/* Find n_in_cell and write appropriate output  */

	first_in_cell = 0;
	n_cells_filled = 0;
	while (first_in_cell < n_pitched) {
		weight = data[first_in_cell].a[BLK_W];
		if (Ctrl->E.active) z_tmp[0] = data[first_in_cell].a[BLK_Z];
		nz = 1;
		if (Ctrl->C.active) {
			j = data[first_in_cell].i / h.nx;
			i = data[first_in_cell].i % h.nx;
			out[0] = GMT_i_to_x (i, h.x_min, h.x_max, h.x_inc, h.xy_off, h.nx);
			out[1] = GMT_j_to_y (j, h.y_min, h.y_max, h.y_inc, h.xy_off, h.ny);
		}
		else {
			out[0] = data[first_in_cell].a[BLK_X];
			out[1] = data[first_in_cell].a[BLK_Y];
		}
		first_in_new_cell = first_in_cell + 1;
		while ( (first_in_new_cell < n_pitched) && (data[first_in_new_cell].i == data[first_in_cell].i) ) {
			weight += data[first_in_new_cell].a[BLK_W];	/* Summing up weights */
			if (mode_xy) {
				out[0] += data[first_in_new_cell].a[BLK_X];
				out[1] += data[first_in_new_cell].a[BLK_Y];
			}
			if (Ctrl->E.active) {	/* Must get a temporary copy of the sorted z array */
				z_tmp[nz] = data[first_in_new_cell].a[BLK_Z];
				nz++;
				if (nz == nz_alloc) {
					nz_alloc += GMT_CHUNK;
					z_tmp = (double *) GMT_memory ((void *)z_tmp, nz_alloc, sizeof (double), GMT_program);
				}
			}
			first_in_new_cell++;
		}
		n_in_cell = first_in_new_cell - first_in_cell;
		if (n_in_cell > 2) {	/* data are already sorted on z; get z mode  */
			out[2] = weighted_mode (&data[first_in_cell], weight, n_in_cell, 2);
			if (Ctrl->Q.active) {
				i_n_in_cell = 1.0 / n_in_cell;
				out[0] *= i_n_in_cell;
				out[1] *= i_n_in_cell;
			}
			else if (mode_xy) {
				qsort((void *)&data[first_in_cell], n_in_cell, sizeof (struct BLK_DATA), BLK_compare_x);
				out[0] = weighted_mode (&data[first_in_cell], weight, n_in_cell, 0);

				qsort((void *)&data[first_in_cell], n_in_cell, sizeof (struct BLK_DATA), BLK_compare_y);
				out[1] = weighted_mode (&data[first_in_cell], weight, n_in_cell, 1);
			}
		}
		else if (n_in_cell == 2) {
			if (data[first_in_cell].a[BLK_W] > data[first_in_cell+1].a[BLK_W]) {
				out[2] = data[first_in_cell].a[BLK_Z];
				if (Ctrl->Q.active) {
					out[0] *= 0.5;
					out[1] *= 0.5;
				}
				else if (mode_xy) {
					out[0] = data[first_in_cell].a[BLK_X];
					out[1] = data[first_in_cell].a[BLK_Y];
				}
			}
			else if (data[first_in_cell].a[BLK_W] < data[first_in_cell+1].a[BLK_W]) {
				out[2] = data[first_in_cell+1].a[BLK_Z];
				if (Ctrl->Q.active) {
					out[0] *= 0.5;
					out[1] *= 0.5;
				}
				else if (mode_xy) {
					out[0] = data[first_in_cell+1].a[BLK_X];
					out[1] = data[first_in_cell+1].a[BLK_Y];
				}
			}
			else {
				if (mode_xy) {	/* Need average location */
					out[0] *= 0.5;
					out[1] *= 0.5;
				}
				out[2] = 0.5 * (data[first_in_cell].a[BLK_Z] + data[first_in_cell+1].a[BLK_Z]);
			}
		}
		else
			out[2] = data[first_in_cell].a[BLK_Z];

		if (Ctrl->E.active) {
			out[4] = z_tmp[0];	/* Low value */
			out[5] = z_tmp[nz-1];	/* High value */
			/* Turn z_tmp into absolute deviations from the mode (out[2]) */
			if (nz > 1) {
				for (index = 0; index < nz; index++) z_tmp[index] = fabs (z_tmp[index] - out[2]);
				qsort ((void *)z_tmp, nz, sizeof (double), GMT_comp_double_asc);
				out[3] = (nz%2) ? z_tmp[nz/2] : 0.5 * (z_tmp[(nz-1)/2] + z_tmp[nz/2]);
				out[3] *= 1.4826;	/* This will be LMS MAD-based scale */
			}
			else
				out[3] = GMT_d_NaN;
		}
		if (Ctrl->W.weighted[GMT_OUT]) out[w_col] = weight;
		
		GMT_output (GMT_stdout, n_out, out);

		n_cells_filled++;
		first_in_cell = first_in_new_cell;
	}

	if (gmtdefs.verbose) fprintf(stderr,"%s: N_cells_filled: %u\n", GMT_program, n_cells_filled);

	GMT_free ((void *)data);
	if (Ctrl->E.active) GMT_free ((void *)z_tmp);

        GMT_end (argc, argv);

	Free_Blockmode_Ctrl (Ctrl);	/* Deallocate control structure */

	exit (EXIT_SUCCESS);
}

double	weighted_mode (struct BLK_DATA *d, double wsum, size_t n, int k)
{
	/* Estimate mode by finding a maximum in the estimated
	   pdf of weighted data.  Estimate the pdf as the finite 
	   difference of the cumulative frequency distribution 
	   over points from i to j.  This has the form top/bottom,
	   where top is the sum of the weights from i to j, and
	   bottom is (data[j] - data[i]).  Strategy is to start
	   with i=0, j=n-1, and then move i or j toward middle
	   while j-i > n/2 and bottom > 0.  At end while, midpoint
	   of range from i to j is the mode estimate.  Choose 
	   to move either i or j depending on which one will
	   cause greatest increase in pdf estimate.  If a tie,
	   move both.
	   
	   Strictly, the pdf estimated this way would need to be
	   scaled by (1/wsum), but this is constant so we don't
	   use it here, as we are seeking a relative minimum.  
	   
	   I assumed n > 2 when I wrote this.  */
	   
	double	top, topj, topi, bottomj, bottomi, pj, pi;
	size_t	i, j, nh;

	i = 0;
	j = n - 1;
	nh = n / 2;
	top = wsum;

	while (j-i > nh) {
		topi = top - d[i].a[BLK_W];
		topj = top - d[j].a[BLK_W];
		bottomi = d[j].a[k] - d[i+1].a[k];
		bottomj = d[j-1].a[k] - d[i].a[k];

		if (bottomj == 0.0) {
			return (d[j-1].a[k]);
		}
		else if (bottomi == 0.0) {
			return (d[i+1].a[k]);
		}
		else {
			pi = topi/bottomi;
			pj = topj/bottomj;
			if (pi > pj) {
				i++;
				top = topi;
			}
			else if (pi < pj) {
				j--;
				top = topj;
			}
			else {
				top -= (d[i].a[BLK_W] + d[j].a[BLK_W]);
                                i++;
                                j--;
			}
		}
	}
	return (0.5*(d[j].a[k] + d[i].a[k]));
}

#include "block_subs.c"
