/*--------------------------------------------------------------------
 *    $Id: blockmean.c,v 1.5.4.3 2002/02/27 17:58:55 pwessel Exp $
 *
 *	Copyright (c) 1991-2002 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
 *--------------------------------------------------------------------*/

/*  
   blockmean.c
   reads x, y, 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.
   
   Latest method uses a hash table and linked lists if the region will
   not simply fit in core memory.
      
   Author:      Walter H. F. Smith
   Version:     3.0, testing hash tables
   Date:        4 October, 1988
   Modified:	26 April 1991 by WHFS for gmt v2.0
   Modified:	3 Jan 1995 by PW for gmt 3.0
   Modified:	3 May 1998 by PW for gmt 3.1
   Modified:	18 Oct 1999 by PW to add -S
   Modified:	3.3.5: 10 Jul 2000 by PW to add -L
   Version:	3.4: 01-MAR-2001 by PW, Use -F instead of -N, and add -C
   Version:	3.4.1
*/

#include "gmt.h"

struct	WXYZ {
	double w;
	double	x;
	double	y;
	double	z;
}	*wxyz;	

struct LIST {
	int		key;
	struct WXYZ	sum;
	struct LIST	*p_right;
	struct LIST	*p_left;
}	**hash, *this_box;

int	set_up_arrays(int n_x, int n_y, int *method);
void	write_output_and_free_space (int n_x, int n_y, int method, int report_weight, BOOLEAN mean_xy, double west, double delta_x, double south, double delta_y, double xy_off, int *n_cells_filled);
struct	LIST *find_box(int i, int j, int method);
BOOLEAN sum_only = FALSE;

main (int argc, char **argv)
{

	BOOLEAN	error, weighted, offset, report_weight, nofile = TRUE, done = FALSE, first = TRUE, mean_xy = TRUE;
	
	FILE *fp = NULL;
	
	double	west, east, south, north, delta_x, delta_y, del_off, xy_off, *in;
	double	idx, idy, weight;
	
	int	i, j, method, n_x, n_y, ij, n_expected_fields, n_fields, n_req;
	int	n_lost, n_read, n_pitched, n_cells_filled, n_files = 0, fno, n_args;
	
	char	modifier, line[BUFSIZ], format[BUFSIZ];

	argc = GMT_begin (argc, argv);

	west = east = south = north = delta_x = delta_y = 0.0;
	del_off = 0.5;
	error = weighted = offset = report_weight = 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 '\0':
                                      error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
                                      break;
                              
				/* Supplemental parameters */
                              
				case 'b':
					error += GMT_io_selection (&argv[i][2]);
					break;
					
				case 'C':
					mean_xy = FALSE;	/* Report center of block instead */
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &delta_x, &delta_y);
					break;
				case 'L':
					GMT_geographic_in = TRUE;
					break;
				case 'N':	/* Backward compatible with 3.3.6 */
				case 'F':
					offset = TRUE;
					break;
				case 'S':
					sum_only = TRUE;
					break;
				case 'W':
					if ( (modifier = argv[i][2]) == 'i' || modifier == 'I') {
						weighted = TRUE;
						report_weight = FALSE;
					}
					else if (modifier == 'O' || modifier == 'o') {
						report_weight = TRUE;
						weighted = FALSE;
					}
					else
						weighted = report_weight = TRUE;
					break;
					
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}

	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "blockmean %s - Block averaging by L2 norm\n\n", GMT_VERSION);
		fprintf (stderr, "usage: blockmean [infile(s)] -I<xinc[m]>[/<yinc>[m]] -R<west/east/south/north>\n");
		fprintf (stderr, "\t[-C] [-F] [-H[<nrec>]] [-L] [-S] [-V] [-W[i][o] ] [-:] [-bi[s][<n>]] [-bo[s][<n>]]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t-I sets Increment of the grid; enter xinc, optionally xinc/yinc.\n");
		fprintf (stderr, "\t   Default is yinc = xinc.  Append an m [or c] to xinc or yinc to indicate minutes [or seconds],\n");
		fprintf (stderr, "\t   e.g., -I10m/5m grids longitude every 10 minutes, latitude every 5 minutes.\n");
		GMT_explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-C Output center of block and mean z-value.  [Default outputs mean x-y location]\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-L means that x is longitude, i.e. assumed to be periodic in 360\n");
		fprintf (stderr, "\t-S report block sums rather than mean values [Default is mean values].\n");
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W sets Weight options.  -WI reads Weighted Input (4 cols: x,y,z,w) but writes only (x,y,z) Output.\n");
		fprintf (stderr, "\t   -WO reads unWeighted Input (3 cols: x,y,z) but reports sum (x,y,z,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 columns (4 if -W is set).\n");
		GMT_explain_option ('o');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	
	if (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
		error++;
	}
	if (delta_x <= 0.0 || delta_y <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", GMT_program);
		error++;
	}
	
	if (GMT_io.binary[0] && gmtdefs.io_header) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
	n_req = (weighted) ? 4 : 3;
	if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = n_req;
	if (GMT_io.binary[0] && n_req > GMT_io.ncol[0]) {
		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 .gmtcommands */
	
	if (GMT_io.binary[0] && gmtdefs.verbose) {
		char *type[2] = {"double", "single"};
		fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]);
	}

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif
	idx = 1.0 / delta_x;
	idy = 1.0 / delta_y;
	n_x = irint ((east - west) * idx) + 1;
	n_y = irint ((north - south) * idy) + 1;
	
	if (offset) {
		n_x--;
		n_y--;
		del_off = 0.0;
	}
	xy_off = 0.5 - del_off;	/* Use to calculate mean location of block */
	
	if ( (set_up_arrays(n_x, n_y, &method) ) ) exit (EXIT_FAILURE);
	
	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, west, east, south, north, n_x, n_y);
		if (method) {
			fprintf (stderr, "%s: Direct method requires %d bytes for array.\n", GMT_program, (n_x*n_y*sizeof(struct WXYZ)));
			fprintf (stderr, "%s: Constructing hash table and linked lists.\n", GMT_program);
		}
		else {
			fprintf (stderr, "%s: Region fits inside core memory.\n", GMT_program);
			fprintf (stderr, "%s: Using %d bytes for array.\n", GMT_program, (n_x*n_y*sizeof(struct WXYZ)));
		}
	}

	n_read = n_pitched = 0;

	n_expected_fields = (GMT_io.binary[0]) ? GMT_io.ncol[0] : 3 + weighted;
	
	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 (0);
#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 (gmtdefs.io_header) {
			for (i = 0; i < gmtdefs.n_header_recs; i++) {
				GMT_fgets (line, BUFSIZ, fp);
				line[strlen(line)-1] = 0;
				if (first) (report_weight && !(weighted)) ? fprintf (GMT_stdout, "%s weights\n", line) : fprintf (GMT_stdout, "%s\n", line);
			}
			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 %d\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++;
			
			i = (int)floor(((in[0] - west) * idx) + del_off);
			if ( i < 0 || i >= n_x ) continue;
			j = (int)floor(((in[1] - south) * idy) + del_off);
			if ( j < 0 || j >= n_y ) continue;
			weight = (weighted) ? in[3] : 1.0;
			if (method) {
				this_box = find_box(i, j, method);
				this_box->key = (n_x > n_y) ? j : i;
				if (mean_xy) {
					this_box->sum.x += (in[0]*weight);
					this_box->sum.y += (in[1]*weight);
				}
				else if (this_box->sum.w == 0.0) {	/* So we only do this once */
					this_box->sum.x = west  + (i + xy_off) * delta_x;
					this_box->sum.y = south + (j + xy_off) * delta_y;
				}
				this_box->sum.w += weight;
				this_box->sum.z += (in[2]*weight);
			}
			else {
				ij = i * n_y + j;
				wxyz[ij].w += weight;
				if (mean_xy) {
					wxyz[ij].x += (in[0]*weight);
					wxyz[ij].y += (in[1]*weight);
				}
				wxyz[ij].z += (in[2]*weight);
			}
			n_pitched++;
		}
		if (fp != GMT_stdin) GMT_fclose(fp);
		
	}

	write_output_and_free_space (n_x, n_y, method, report_weight, mean_xy, west, delta_x, south, delta_y, xy_off, &n_cells_filled);

	n_lost = n_read - n_pitched;
	if (gmtdefs.verbose) {
		fprintf (stderr, "%s: N read: %d\tN used: %d\tN outside: %d\tN cells filled: %d\n",
		GMT_program, n_read, n_pitched, n_lost, n_cells_filled);
	}

        GMT_end (argc, argv);
}

int	set_up_arrays(int n_x, int n_y, int *method)
{

	if ( (wxyz = (struct WXYZ *) calloc ((size_t)(n_x * n_y), sizeof (struct WXYZ) ) ) == NULL) {

		/* Full matrix did not fit in core, or failed for other reason.  Try hash table:  */

		if (n_x > n_y) {	/* Hash on x (i); linked list on y (j)  */
			*method = 1;
			hash = (struct LIST **) GMT_memory (VNULL, (size_t)n_x, sizeof (struct LIST *), GMT_program);
		}
		else {	/* Hash on y (j); linked list on x (i)  */
			*method = -1;
			hash = (struct LIST **) GMT_memory (VNULL, (size_t)n_y, sizeof (struct LIST *), GMT_program);
		}
	}
	else {
		/* Easy method fits in memory.  */
		*method = 0;
	}
	return(0);
}

void	write_output_and_free_space (int n_x, int n_y, int method, int report_weight, BOOLEAN mean_xy, double west, double delta_x, double south, double delta_y, double xy_off, int *n_cells_filled)
{
	int	loop, limit, n_out, i, j;
	double	out[4], iw;
	
	n_out = (report_weight) ? 4 : 3;

	*n_cells_filled = 0;
	
	if (method) {	/* We have to loop over the linked lists  */
	
		if (method == 1)
			limit = n_x;
		else
			limit = n_y;
		
		for (loop = 0; loop < limit; loop++) {

			/* If this list is empty, continue:  */
			
			if (hash[loop] == NULL) continue;
			
			/* Go to the leftmost box in the list  */
			
			this_box = hash[loop];
			while (this_box->p_left) this_box = this_box->p_left;
			
			/* While we are at a box, write it, and move right  */
			
			while (this_box) {
			
				(*n_cells_filled)++;
				out[3] = this_box->sum.w;
				iw = 1.0 / out[3];
				if (mean_xy) {
					out[0] = this_box->sum.x * iw;
					out[1] = this_box->sum.y * iw;
				}
				else {	/* This was only set once to the center location */
					out[0] = this_box->sum.x;
					out[1] = this_box->sum.y;
				}
				out[2] = (sum_only) ? this_box->sum.z : this_box->sum.z * iw;

				GMT_output (GMT_stdout, n_out, out);
				
				if (this_box->p_right) {
					this_box = this_box->p_right;
					GMT_free ((void *)this_box->p_left);
				}
				else {
					GMT_free ((void *)this_box);
					this_box = NULL;
				}
			}
		}
		GMT_free ((void *)hash);
	}
	
	else {	/* We have a simple array in memory  */
	
		limit = n_x * n_y;
		
		for (loop = 0; loop < limit; loop++) {
		
			if (wxyz[loop].w == 0.0) continue;
			
			(*n_cells_filled)++;
			out[3] = wxyz[loop].w;
			iw = 1.0 / out[3];
			if (mean_xy) {
				out[0] = wxyz[loop].x * iw;
				out[1] = wxyz[loop].y * iw;
			}
			else {	/* Use block center */
				j = loop % n_y;
				i = loop / n_y;
				out[0] = west  + (i + xy_off) * delta_x;
				out[1] = south + (j + xy_off) * delta_y;
			}
			out[2] = (sum_only) ? wxyz[loop].z : wxyz[loop].z * iw;

			GMT_output (GMT_stdout, n_out, out);
		}
		
		free ((void *)wxyz);	/* Free because allocated with calloc instead of GMT_memory */
	}
}
		
struct LIST *find_box(int i, int j, int method)
{
	int	hash_key, list_key;
	struct LIST *current_box, *temp_box;

	static int last_i = -1;
	static int last_j = -1;
	static struct LIST *last_box = NULL;
	
	if ( (i == last_i) && (j == last_j) ) return (last_box);
	
	/* Get here when we have to search  */
	
	if (method > 0) {
		hash_key = i;
		list_key = j;
	}
	else {
		hash_key = j;
		list_key = i;
	}
	
	current_box = hash[hash_key];
	
	if (current_box) {
	
		/* Hash table has an entry into a linked list;
			if it doesn't match the list_key, search the list */
		
		if (current_box->key < list_key) {
		
			/* Current's key too low; move right while necessary  */
		
			while ( (current_box->p_right) && (current_box->p_right->key <= list_key) )
				current_box = current_box->p_right;
			
			if (current_box->key < list_key) {
			
				/* Current's key still too low  */
			
				if (current_box->p_right) {
				
					/* Next's key too high; insert a box in between  */
					
					temp_box = (struct LIST *) GMT_memory (VNULL, (size_t)1, sizeof (struct LIST), GMT_program);
					temp_box->p_right = current_box->p_right;
					current_box->p_right->p_left = temp_box;
					current_box->p_right = temp_box;
					temp_box->p_left = current_box;
					current_box = temp_box;
				}
				else {
					/* There is no next; make new box at end of list  */
					
					temp_box = (struct LIST *) GMT_memory (VNULL, (size_t)1, sizeof (struct LIST), GMT_program);
					temp_box->p_left = current_box;
					current_box->p_right = temp_box;
					current_box = temp_box;
				}
			}
		}
		
		if (current_box->key > list_key) {
		
			/* Current's key too high; move left while necessary  */
		
			while ( (current_box->p_left) && (current_box->p_left->key >= list_key) )
				current_box = current_box->p_left;
				
			if (current_box->key > list_key) {
			
				/* Current's key still too high  */
			
				if (current_box->p_left) {
				
					/* Next's key too low; insert a box in between  */
					
					temp_box = (struct LIST *) GMT_memory (VNULL, (size_t)1, sizeof (struct LIST), GMT_program);
					temp_box->p_left = current_box->p_left;
					current_box->p_left->p_right = temp_box;
					current_box->p_left = temp_box;
					temp_box->p_right = current_box;
					current_box = temp_box;
				}
				else {
					/* There is no next; make new box at end of list  */
					
					temp_box = (struct LIST *) GMT_memory (VNULL, (size_t)1, sizeof (struct LIST), GMT_program);
					temp_box->p_right = current_box;
					current_box->p_left = temp_box;
					current_box = temp_box;
				}
			}
		}
			

			
		
		
	}
	else  {	/* Hash table is NULL; create first box in linked list  */
	
		current_box = (struct LIST *) GMT_memory (VNULL, (size_t)1, sizeof (struct LIST), GMT_program);
	}
		
	/* Now we set the state of the static variables and return  */
	
	last_i = i;
	last_j = j;
	last_box = current_box;
	hash[hash_key] = current_box;	/* So we enter lists where we were last time  */
	return (current_box);
}