/*--------------------------------------------------------------------
 *    $Id: filter1d.c,v 1.17 2004/08/14 00:38:19 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
 *--------------------------------------------------------------------*/
/*
 * filter1d will read N columns of data from file/GMT_stdin and return
 * filtered output at user-selected positions.  The time variable
 * can be in any specified column of the input.  Several filters
 * are available, with robustness as an option.
 *
 * Filters:
 * Convolutions: Boxcar, Cosine Arch, Gaussian, Median, or Mode.
 * Geospatial:   Median, Mode, Extreme values.
 * Robust:	Option replaces outliers with medians in filter.
 * Output:	At input times, or from t_start to t_stop by t_int.
 * Lack:	Option checks for data gaps in input series.
 * Symmetry:	Option checks for asymmetry in filter window.
 * Quality:	Option checks for low mean weight in window.
 *
 * Author:	Walter H. F. Smith
 * Date:	13 January, 1989
 * Version:	4
 */

#include "gmt.h"

#define BOXCAR 0
#define COS_ARCH 1
#define GAUSSIAN 2
#define MEDIAN 3
#define MODE 4
#define LOWER_ALL 5
#define LOWER_POS 6
#define UPPER_ALL 7
#define UPPER_NEG 8
#define N_FILTERS 9
#define CONVOLVE 2		/* If filter_type > CONVOLVE then a MEDIAN, MODE, or EXTREME filter is selected  */

PFD get_weight[3];		/* Selects desired weight function.  */

double	boxcar_weight (double radius, double half_width);	/* Returns the weight for a given delta_time  */
double	cosine_weight (double radius, double half_width);
double	gaussian_weight (double radius, double half_width);

int	load_data_and_check_extrema (void);		/* Does just what it says  */
int	set_up_filter (void);	/* Creates filter weights or reads them from file and sets start & stop  */
int	do_the_filter (void);	/* Does the job  */
int	allocate_space (void);	/* Does just what it says  */
void	free_space (void);		/* Does just what it says  */
void	get_robust_estimates (int j, int n, int both);	/* Does just what it says  */
void	allocate_more_data_space (void);	/* Does just what it says  */
void	allocate_more_work_space (void);	/* Does just what it says  */
BOOLEAN	lack_check (int i_col, int left, int right);		/* Tests for lack of data condition  (optional)  */

double	*f_wt;			/* Pointer for array of filter coefficients  */
double	*wt_sum;		/* Pointer for array of weight sums [each column]  */
double	*data_sum;		/* Pointer for array of data * weight sums [columns]  */
double	*min_loc;		/* Pointer for array of values, one per [column]  */
double	*max_loc;
double	*last_loc;
double	*this_loc;
double	*min_scl;
double	*max_scl;
double	*last_scl;
double	*this_scl;
double	**work;			/* Pointer to array of pointers to doubles for work  */	
double	**data;			/* Pointer to array of pointers to doubles for data  */	
double	dt;			/* Delta time resolution for filter computation  */
double	q_factor;		/* Quality level for mean weights or n in median  */
double filter_width;		/* Full width of filter in user's units */
double half_width;		
double ignore_val;		/* Value to skip over in input stream  */	
double t_start;			/* x-value of first output point */
double t_stop;			/* x-value of last output point */
double t_int;			/* Output interval */
double sym_coeff = 0.0;		/* Symmetry coefficient  */
double lack_width = 0.0;	/* Lack of data width  */
double extreme = 0.0;		/* Extreme value [for the l|L|u|U filter] */

int	*n_this_col;		/* Pointer to array of counters [one per column]  */
int	*n_left;		/* Pointer to array of counters [one per column]  */
int	*n_right;		/* Pointer to array of counters [one per column]  */

int	n_rows;			/* Number of rows of input  */
int	n_cols = 2;		/* Number of columns of input  */
int	t_col = 0;		/* Column of time abscissae (independent variable)  */
int	n_f_wts;		/* Number of filter weights  */
int	half_n_f_wts;		/* Half the number of filter weights  */
int	n_row_alloc = GMT_CHUNK;		/* Number of rows of data to allocate  */
int	n_work_alloc = GMT_CHUNK;	/* Number of rows of workspace to allocate  */
int	filter_type = 0;	/* Flag indicating desired filter type  */
int	kind = 0;			/* -1 skip +ve, +1 skip -ve, else use all  [for the l|L|u|U filter] */
int	way;			/* -1 find minimum, +1 find maximum  [for the l|L|u|U filter] */

BOOLEAN	*good_one;		/* Pointer to array of logicals [one per column]  */

BOOLEAN	use_ends = FALSE;	/* True to start/stop at ends of series instead of 1/2 width inside  */
BOOLEAN check_asym = FALSE;	/* TRUE to test whether the data are asymmetric about the output time  */
BOOLEAN check_lack = FALSE;	/* TRUE to test for lack of data (gap) in the filter window */
BOOLEAN check_q = FALSE;	/* TRUE to test average weight or N in median */
BOOLEAN robust = FALSE;		/* Look for outliers in data when TRUE */
BOOLEAN equidist = TRUE;	/* Data is evenly sampled in t */
BOOLEAN out_at_time = FALSE;	/* TRUE when output is required at evenly spaced intervals */
BOOLEAN ignore_test = FALSE;	/* TRUE when input values that eqaul ignore_val should be discarded */
BOOLEAN error = FALSE;
BOOLEAN f_operator = FALSE;	/* TRUE if custom weights coefficients sum to zero */
BOOLEAN n_cols_set = FALSE;

char	buffer[BUFSIZ];		/* Used to scan input lines  */

FILE	*fp_in = NULL;		/* File pointer to data file or GMT_stdin */
FILE	*fp_wt = NULL;		/* File pointer to custom weight coefficients (optional) */

int GMT_mode_selection = 0, GMT_n_multiples = 0;

main (int argc, char **argv)
{

	int	i;
	char	c;
	
	argc = GMT_begin (argc, argv);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
		
				/* Common parameters */
			
				case 'H':
				case 'V':
				case ':':
				case 'b':
				case 'f':
				case '\0':
					error += GMT_get_common_args (argv[i], 0, 0, 0, 0);
					break;
				
				/* Supplemental parameters */
				
				case 'F':	/* Filter selection  */
					switch (argv[i][2]) {
						case 'B':
							robust = TRUE;
						case 'b':
							filter_type = BOXCAR;
							filter_width = atof (&argv[i][3]);
							break;
						case 'C':
							robust = TRUE;
						case 'c':
							filter_type = COS_ARCH;
							filter_width = atof (&argv[i][3]);
							break;
						case 'G':
							robust = TRUE;
						case 'g':
							filter_type = GAUSSIAN;
							filter_width = atof (&argv[i][3]);
							break;
						case 'M':
							robust = TRUE;
						case 'm':
							filter_type = MEDIAN;
							filter_width = atof (&argv[i][3]);
							break;
						case 'P':
							robust = TRUE;
						case 'p':
							filter_type = MODE;
							filter_width = atof (&argv[i][3]);
							c = argv[i][strlen(argv[i]-1)];
							if (c == '-') GMT_mode_selection = -1;
							if (c == '+') GMT_mode_selection = +1;
							break;
						case 'l':
							filter_type = LOWER_ALL;
							filter_width = atof (&argv[i][3]);
							way = -1;
							extreme = DBL_MAX;
							break;
						case 'L':
							filter_type = LOWER_POS;
							filter_width = atof (&argv[i][3]);
							way = -1;
							kind = +1;
							break;
						case 'u':
							filter_type = UPPER_ALL;
							filter_width = atof (&argv[i][3]);
							way = +1;
							extreme = -DBL_MAX;
							break;
						case 'U':
							filter_type = UPPER_NEG;
							filter_width = atof (&argv[i][3]);
							way = +1;
							kind = -1;
							break;
						case 'F':
							robust = TRUE;
						case 'f':
							if ((fp_wt = fopen(&argv[i][3],"r")) == NULL) {
								fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, &argv[i][3]);
								exit (EXIT_FAILURE);
							}
							break;
						default:
							error = TRUE;
							break;
					}
					break;
					
				case 'D':
					dt = atof (&argv[i][2]);
					equidist = FALSE;
					break;
				case 'E':
					use_ends = TRUE;
					break;
				case 'I':
					ignore_val = atof (&argv[i][2]);
					ignore_test = TRUE;
					break;
				case 'L':
					check_lack = TRUE;
					lack_width = atof (&argv[i][2]);
					break;
				case 'N':
					sscanf (&argv[i][2], "%d/%d", &n_cols, &t_col);
					n_cols_set = TRUE;
					break;
				case 'Q':
					q_factor = atof(&argv[i][2]);
					check_q = TRUE;
					break;
				case 'S':
					check_asym = TRUE;
					sym_coeff = atof (&argv[i][2]);
					break;
				case 'T':
					sscanf (&argv[i][2], "%lf/%lf/%lf", &t_start, &t_stop, &t_int);
					out_at_time = TRUE;
					if (t_int == 0.0) error = TRUE;
					break;
				default:
					error = TRUE;
                                        GMT_default_error (argv[i][1]);
					break;
			}
		}
		else {
			if ((fp_in = GMT_fopen(argv[i], GMT_io.r_mode)) == NULL) {
				fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[i]);
				exit (EXIT_FAILURE);
			}
		}
	}

	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "filter1d %s - Time domain filtering of 1-D time series\n\n", GMT_VERSION);
		fprintf (stderr, "usage: filter1d [infile] -F<type><width>[<mode>] [-D<increment>] [-E] [-H[<nrec>]]\n");
		fprintf (stderr, "\t[-I<ignore_val>] [-L<lack_width>] [-N<n_cols>/<t_col>] -Q<q_factor>]\n");
		fprintf (stderr, "\t[[-S<symmetry_factor>] [-T<start>/<stop>/<int>] [-V] [-bi[s][<n>]]\n\t[-bo[s][<n>]] [-f[i|o]<colinfo>]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t-F sets Filtertype.  Choose from convolution and non-convolution filters.\n");
		fprintf (stderr, "\t   and append full filter <width> in same units as time column,\n");
		fprintf (stderr, "\t   Convolution filters:\n");
		fprintf (stderr, "\t     b: Boxcar : Weights are equal.\n");
		fprintf (stderr, "\t     c: Cosine arch : Weights given by cosine arch.\n");
		fprintf (stderr, "\t     g: Gaussian : Weights given gy Gaussian function.\n");
		fprintf (stderr, "\t     f<name>: Custom : Weights given in one-column file <name>.\n");
		fprintf (stderr, "\t   Non-convolution filters:\n");
		fprintf (stderr, "\t     m: Median : Return the median value.\n");
		fprintf (stderr, "\t     p: Maximum likelihood probability (mode) estimator : Return the mode.\n");
		fprintf (stderr, "\t        By default, we return the average if more than one mode is found.\n");
		fprintf (stderr, "\t        Append - or + to the width to instead return the smallest or largest mode.\n");
		fprintf (stderr, "\t     l: Lower : Return minimum of all points.\n");
		fprintf (stderr, "\t     L: Lower+ : Return minimum of all +ve points.\n");
		fprintf (stderr, "\t     u: Upper : Return maximum of all points.\n");
		fprintf (stderr, "\t     U: Upper- : Return maximum of all -ve points.\n");
		fprintf (stderr, "\t   Upper case type B, C, G, M, P, F will use robust filter versions:\n");
		fprintf (stderr, "\t   i.e., replace outliers (2.5 L1 scale off median) with median during filtering.\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		
		fprintf (stderr, "\t-D used when series is NOT equidistantly sampled.\n");
		fprintf (stderr, "\t   Then <increment> will be the abscissae resolution, i.e. all abscissae\n");
		fprintf (stderr, "\t   will be rounded off to a multiple of <increment>.\n");
		fprintf (stderr, "\t-E include Ends of time series in output.  Default loses half_width at each end.\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I to ignore values; If an input value == <ignore_val> it will be set to NaN.\n");
		fprintf (stderr, "\t-L checks for Lack of data condition.  If input data has a gap exceeding\n");
		fprintf (stderr, "\t   <width> then no output will be given at that point.  Default does not check Lack.\n");
		fprintf (stderr, "\t-N sets # columns in input and which column contains the independent\n");
		fprintf (stderr, "\t   variable (time). The left-most column is # 0, the right-most is # (<n_cols> - 1).\n");
		fprintf (stderr, "\t   Default is <n_cols> = 2, <t_col> = 0; i.e. file has t, f(t) pairs.\n");
		fprintf (stderr, "\t-Q assess Quality of output value by checking mean weight in convolution.\n");
		fprintf (stderr, "\t   Enter <q_factor> between 0 and 1.  If mean weight < q_factor, output is suppressed\n");
		fprintf (stderr, "\t   at this point.  Default does not check Quality.\n");
		fprintf (stderr, "\t-S checks symmetry of data about window center.  Enter a factor\n");
		fprintf (stderr, "\t   between 0 and 1.  If ( (abs(n_left - n_right)) / (n_left + n_right) ) > factor,\n");
		fprintf (stderr, "\t   then no output will be given at this point.  Default does not check Symmetry.\n");
		fprintf (stderr, "\t-T make evenly spaced timesteps from <start> to <stop> by <int>. Default uses input times.\n");
		GMT_explain_option ('V');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		GMT_explain_option ('o');
		GMT_explain_option ('n');
		GMT_explain_option ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}

	/* Check arguments */

	if (out_at_time && (t_stop - t_start) < filter_width) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -T option:  Output interval < filterwidth\n", GMT_program);
		error++;
	}
	if (check_lack && (lack_width < 0.0 || lack_width > filter_width) ) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -L option:  Unreasonable lack-of-data interval\n", GMT_program);
		error++;
	}
	if (check_asym && (sym_coeff < 0.0 || sym_coeff > 1.0) ) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Enter a factor between 0 and 1\n", GMT_program);
		error++;
	}
	if (filter_type < 0 || filter_type >= N_FILTERS) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -F option.  Correct syntax:\n", GMT_program);
		fprintf (stderr, "-FX<width>, X one of BbCcGgMmPpFflLuU\n");
		error++;
	}
	if (t_col >= n_cols) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -N option.  time column exceeds number of columns\n", GMT_program);
		error++;
	}
	if (check_q && q_factor < 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Q option:  Enter a factor between 0 and 1\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] < n_cols) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data must have at least %d fields\n", GMT_program, n_cols);
		error++;
	}
	if (GMT_io.binary[GMT_IN] && !n_cols_set && GMT_io.ncol[GMT_IN] == 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Must specify number of columns in binary input data (-bi)\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]]);
	}


#ifdef SET_IO_MODE
	GMT_setmode (GMT_OUT);
#endif
	if (filter_type > CONVOLVE) robust = FALSE;

	if (allocate_space () ) {
		fprintf(stderr,"%s: fatal error memory allocation.\n", GMT_program);
		free_space ();
		exit (EXIT_FAILURE);
	}

	if (load_data_and_check_extrema () ) {
		fprintf(stderr,"%s: fatal error during data read.\n", GMT_program);
		free_space ();
		exit (EXIT_FAILURE);
	}

	if (set_up_filter () ) {
		fprintf(stderr,"%s: fatal error during coefficient setup.\n", GMT_program);
		free_space ();
		exit (EXIT_FAILURE);
	}

	if (do_the_filter () ) {
		fprintf(stderr,"%s: fatal error in filtering routine.\n", GMT_program);
		free_space ();
		exit (EXIT_FAILURE);
	}

	free_space ();
	
	if (GMT_n_multiples > 0 && gmtdefs.verbose) fprintf (stderr, "%s: WARNING: %d multiple modes found\n", GMT_program, GMT_n_multiples);
	
	GMT_end (argc, argv);
}
	
int load_data_and_check_extrema (void)
{	

	int	i, n_fields, n_expected_fields;
	double	last_time, new_time, *in;

	if (fp_in == NULL) {
		fp_in = GMT_stdin;
#ifdef SET_IO_MODE
		GMT_setmode (GMT_IN);
#endif
	}
	
	if (gmtdefs.io_header[GMT_IN]) {
		for (i = 0; i < gmtdefs.n_header_recs; i++) {
			GMT_fgets (buffer, BUFSIZ, fp_in);
			if (gmtdefs.io_header[GMT_OUT]) fputs (buffer, GMT_stdout);
		}
	}

	n_rows = 0;
	if (robust || (filter_type == MEDIAN) ) {
		for (i = 0; i < n_cols; i++) {
			min_loc[i] = DBL_MAX;
			max_loc[i] = -DBL_MAX;
		}
	}
	last_time = -DBL_MAX;

	if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] > 0)
		n_expected_fields = GMT_io.ncol[GMT_IN];
	else
		n_expected_fields = n_cols;
	
	while ((n_fields = GMT_input (fp_in, &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_rows);
			exit (EXIT_FAILURE);
		}

		for (i = 0; i < n_cols; i++) {
			if (ignore_test && in[i] == ignore_val)
				data[i][n_rows] = GMT_d_NaN;
			else {
				data[i][n_rows] = in[i];
				if (robust || (filter_type == MEDIAN) ) {
					if (in[i] > max_loc[i]) max_loc[i] = in[i];
					if (in[i] < min_loc[i]) min_loc[i] = in[i];
				}
			}
		}
		
		if (!GMT_is_dnan (data[t_col][n_rows]) )
		{	/* Don't n_rows++ if time is bad this row  */
		
			new_time = data[t_col][n_rows];
			n_rows++;
			if (new_time < last_time) {
				fprintf(stderr,"%s: Error! Time decreases at line # %d\n", GMT_program, n_rows);
				fprintf(stderr,"\tUse UNIX utility sort and then try again.\n");
				return(-1);
			}
			last_time = new_time;
		}
		
		if (n_rows == n_row_alloc) {	/* Need more memory */
			n_row_alloc += GMT_CHUNK;
			allocate_more_data_space ();
		}
	}
	
	if (fp_in != GMT_stdin) GMT_fclose (fp_in);
	
	for (i = 0; i < n_cols; i++) data[i] = (double *) GMT_memory ((void *)data[i], (size_t)n_rows, sizeof (double), GMT_program);
	if (gmtdefs.verbose) fprintf (stderr, "%s: Read %d data lines\n", GMT_program, n_rows);

	/* Initialize scale parameters and last_loc based on min and max of data  */
	
	if (robust || (filter_type == MEDIAN) ) {
		for (i = 0; i < n_cols; i++) {
			min_scl[i] = 0.0;
			max_scl[i] = 0.5 * (max_loc[i] - min_loc[i]);
			last_scl[i] = 0.5 * max_scl[i];
			last_loc[i] = 0.5 * (max_loc[i] + min_loc[i]);
		}
	}

	return (0);
}

int	set_up_filter(void)
{
	int	i, i1, i2;
	double	t_0, t_1, time, w_sum;
	BOOLEAN	normalize = FALSE;
	
	
	
	t_0 = data[t_col][0];
	t_1 = data[t_col][n_rows-1];
	if (equidist) dt = (t_1 - t_0) / (n_rows - 1);
	

	if (fp_wt != NULL) {
		f_wt = (double *) GMT_memory (VNULL, (size_t)n_work_alloc, sizeof(double), GMT_program);
		n_f_wts = 0;
		w_sum = 0.0;
		while (fscanf (fp_wt,"%lf\n", &f_wt[n_f_wts]) != EOF) {
			w_sum += f_wt[n_f_wts];
			n_f_wts++;
			if (n_f_wts == n_work_alloc) {	/* Need more memory */
				n_work_alloc += GMT_CHUNK;
				allocate_more_work_space ();
			}
		}
		fclose (fp_wt);
		f_operator = (w_sum == 0.0);	/* If weights sum to zero it is an operator like {-1 1] or [1 -2 1] */
		half_n_f_wts = n_f_wts / 2;
		half_width = half_n_f_wts * dt;
		f_wt = (double *) GMT_memory ((void *)f_wt, (size_t)n_f_wts, sizeof (double), GMT_program);
		if (gmtdefs.verbose) fprintf (stderr, "%s: Read %d filter weights from file.\n", GMT_program, n_f_wts);
	}
	else if (filter_type <= CONVOLVE) {

		get_weight[BOXCAR] = (PFD)boxcar_weight;
		get_weight[COS_ARCH] = (PFD)cosine_weight;
		get_weight[GAUSSIAN] = (PFD)gaussian_weight;
		half_width = 0.5 * filter_width;
		half_n_f_wts = (int)floor (half_width / dt);
		n_f_wts = 2 * half_n_f_wts + 1;
		
		f_wt = (double *) GMT_memory (VNULL, (size_t)n_f_wts, sizeof(double), GMT_program);
		for (i = 0; i <= half_n_f_wts; i++) {
			time = i * dt;
			i1 = half_n_f_wts - i;
			i2 = half_n_f_wts + i;
			f_wt[i1] = f_wt[i2] = ( *get_weight[filter_type]) (time, half_width);
		}
		if (normalize) {
			w_sum = 0.0;
			for (i = 0; i < n_f_wts; i++) w_sum += f_wt[i];
			for (i = 0; i < n_f_wts; i++) f_wt[i] /= w_sum;
		}
	}
	else {
		half_width = 0.5 * filter_width;
	}

	/* Initialize start/stop time */

	if (out_at_time) {
		
		if (use_ends) {
			while (t_start < t_0) t_start += t_int;
			while (t_stop > t_1) t_stop -= t_int;
		}
		else {
			while ( (t_start - half_width) < t_0) t_start += t_int;
			while ( (t_stop + half_width) > t_1) t_stop -= t_int;
		}
	}
	else {
		if (use_ends) {
			t_start = t_0;
			t_stop = t_1;
		}
		else {
			for (i = 0; (data[t_col][i] - t_0) < half_width; i++);
			t_start = data[t_col][i];
			for (i = n_rows - 1; (t_1 - data[t_col][i]) < half_width; i--);
			t_stop = data[t_col][i];
		}
	}

	if (gmtdefs.verbose) fprintf (stderr, "F width: %g Resolution: %g Start: %g Stop: %g\n",
		filter_width, dt, t_start, t_stop);
	
	return(0);
}

int	do_the_filter (void)
{

	int	i_row, i_col, left, right, n_l, n_r, iq, i_f_wt;
	int	i_t_output, n_in_filter, n_for_call, n_good_ones;
	double	time, delta_time, *outval, wt, val, med, scl, small;

	outval = (double *)GMT_memory (VNULL, (size_t)n_cols, sizeof (double), GMT_program);
	
	if(!out_at_time) {	/* Position i_t_output at first output time  */
		for(i_t_output = 0; data[t_col][i_t_output] < t_start; i_t_output++);
		small = (data[t_col][1] - data[t_col][0]);
	}
	else
		small = t_int;
	
	small *= GMT_CONV_LIMIT;
	time = t_start;
	left = right = 0;		/* Left/right end of filter window */

	iq = irint (q_factor);

	/* while (time <= t_stop) { */
	while (time <= (t_stop + small)) {
		while ((time - data[t_col][left]) > half_width) left++;
		while (right < n_rows && (data[t_col][right] - time) <= half_width) right++;
		n_in_filter = right - left;
		if ( (!(n_in_filter)) || (check_lack && ( (filter_width / n_in_filter) > lack_width) ) ) {
			if (out_at_time)
				time += t_int;
			else {
				i_t_output++;
				time = (i_t_output < n_rows) ? data[t_col][i_t_output] : t_stop + 1.0;
			}
			continue;
		}
			
		for (i_col = 0; i_col < n_cols; i_col++) {
			n_this_col[i_col] = 0;
			wt_sum[i_col] = 0.0;
			data_sum[i_col] = 0.0;
			if (i_col == t_col) {
				good_one[i_col] = FALSE;
			}
			else if (check_lack) {
				good_one[i_col] = !(lack_check(i_col, left, right));
			}
			else {
				good_one[i_col] = TRUE;
			}
			if (check_asym) {
				n_left[i_col] = 0;
				n_right[i_col] = 0;
			}
		}
			
		if (robust || filter_type > CONVOLVE) {
			if (n_in_filter > n_work_alloc) {
				n_work_alloc = n_in_filter;
				allocate_more_work_space ();
			}
			for (i_row = left; i_row < right; i_row++) {
				for (i_col = 0; i_col < n_cols; i_col++) {
					if (!(good_one[i_col])) continue;
					if (!GMT_is_dnan (data[i_col][i_row])) {
						work[i_col][n_this_col[i_col]] = data[i_col][i_row];
						n_this_col[i_col]++;
						if (check_asym) {
							if (data[t_col][i_row] < time) n_left[i_col]++;
							if (data[t_col][i_row] > time) n_right[i_col]++;
						}
					}
				}
			}
			if (check_asym) {
				for (i_col = 0; i_col < n_cols; i_col++) {
					if (!(good_one[i_col])) continue;
					n_l = n_left[i_col];
					n_r = n_right[i_col];
					if ((((double)abs(n_l - n_r))/(n_l + n_r)) > sym_coeff) good_one[i_col] = FALSE;
				}
			}
			if ( (filter_type > CONVOLVE) && check_q) {
				for (i_col = 0; i_col < n_cols; i_col++) {
					if (n_this_col[i_col] < iq) good_one[i_col] = FALSE;
				}
			}
				
			for (i_col = 0; i_col < n_cols; i_col++) {
				if (good_one[i_col]) {
					n_for_call = n_this_col[i_col];
					get_robust_estimates (i_col, n_for_call, robust);
				}
			}
				
		}	/* That's it for the robust work  */

		if (filter_type > CONVOLVE) {
			
			/* Need to count how many good ones; use data_sum area  */
				
			n_good_ones = 0;
			for (i_col = 0; i_col < n_cols; i_col++) {
				if (i_col == t_col) {
					data_sum[i_col] = time;
				}
				else if (good_one[i_col]) {
					data_sum[i_col] = this_loc[i_col];
					n_good_ones++;
				}
				else {
					data_sum[i_col] = GMT_d_NaN;
				}
			}
			if (n_good_ones) GMT_output (GMT_stdout, n_cols, data_sum);
		}
			
		else {
			
			if (robust) for (i_col = 0; i_col < n_cols; i_col++) n_this_col[i_col] = 0;
						
			for (i_row = left; i_row < right; i_row++) {
				delta_time = time - data[t_col][i_row];
				i_f_wt = half_n_f_wts + (int)floor(0.5 + delta_time/dt);
				if ( (i_f_wt < 0) || (i_f_wt >= n_f_wts) ) continue;
					
				for(i_col = 0; i_col < n_cols; i_col++) {
					if ( !(good_one[i_col]) ) continue;
					if (!GMT_is_dnan (data[i_col][i_row])) {
						wt = f_wt[i_f_wt];
						val = data[i_col][i_row];
						if (robust) {
							med = this_loc[i_col];
							scl = this_scl[i_col];
							val = ((fabs(val-med)) > (2.5 * scl)) ? med : val;
						}
						else if (check_asym) {	/* This wasn't already done  */
							if (data[t_col][i_row] < time) n_left[i_col]++;
							if (data[t_col][i_row] > time) n_right[i_col]++;
						}
						wt_sum[i_col] += wt;
						data_sum[i_col] += (wt * val);
						n_this_col[i_col]++;
					}
				}
			}
			n_good_ones = 0;
			for (i_col = 0; i_col < n_cols; i_col++) {
				if ( !(good_one[i_col]) ) continue;
				if ( !(n_this_col[i_col]) ) {
					good_one[i_col] = FALSE;
					continue;
				}
				if (check_asym && !(robust) ) {
					n_l = n_left[i_col];
					n_r = n_right[i_col];
					if ((((double)abs(n_l - n_r))/(n_l + n_r)) > sym_coeff) {
						good_one[i_col] = FALSE;
						continue;
					}
				}
				if (check_q && ((wt_sum[i_col] / n_this_col[i_col]) < q_factor)) {
					good_one[i_col] = FALSE;
					continue;
				}
				n_good_ones++;
			}
			if (n_good_ones) {
				for (i_col = 0; i_col < n_cols; i_col++) {
					if (i_col == t_col)
						outval[i_col] = time;
					else if (good_one[i_col])
						outval[i_col] = (f_operator) ? data_sum[i_col] : data_sum[i_col] / wt_sum[i_col];
					else
						outval[i_col] = GMT_d_NaN;
				}
				GMT_output (GMT_stdout, n_cols, outval);
			}
		}

		/* Go to next output time */

		if (out_at_time)
			time += t_int;
		else {
			i_t_output++;
			time = (i_t_output < n_rows) ? data[t_col][i_t_output] : t_stop + 1.0;
		}
	}
	
	GMT_free ((void *)outval);
	
	return(0);
}

int	allocate_space (void)
{

	int	i;

	data = (double **) GMT_memory (VNULL, (size_t)n_cols, sizeof(double *), GMT_program);
	for (i = 0; i < n_cols; i++) data[i] = (double *) GMT_memory (VNULL, (size_t)n_row_alloc, sizeof(double), GMT_program);
	
	wt_sum = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
	data_sum = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
	n_this_col = (int *) GMT_memory (VNULL, (size_t)n_cols, sizeof(int), GMT_program);
	good_one = (int *) GMT_memory (VNULL, (size_t)n_cols, sizeof(int), GMT_program);

	if (check_asym) n_left = (int *) GMT_memory (VNULL, (size_t)n_cols, sizeof(int), GMT_program);
	if (check_asym) n_right = (int *) GMT_memory (VNULL, (size_t)n_cols, sizeof(int), GMT_program);
	
	if (robust || (filter_type > CONVOLVE) ) {	/* Then we need workspace  */

		work = (double **) GMT_memory (VNULL, (size_t)n_cols, sizeof(double *), GMT_program);
		for (i = 0; i < n_cols; i++) work[i] = (double *) GMT_memory (VNULL, (size_t)n_work_alloc, sizeof(double), GMT_program);
		min_loc = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
		max_loc = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
		last_loc = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
		this_loc = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
		min_scl = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
		max_scl = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
		this_scl = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
		last_scl = (double *) GMT_memory (VNULL, (size_t)n_cols, sizeof(double), GMT_program);
	}
	return (0);
}

void	allocate_more_data_space (void)
{
	int	i;

	for (i = 0; i < n_cols; i++) data[i] = (double *) GMT_memory ((void *)data[i], (size_t)n_row_alloc, sizeof(double), GMT_program);
}

void	allocate_more_work_space (void)
{
	int	i;

	for (i = 0; i < n_cols; i++) work[i] = (double *) GMT_memory ((void *)work[i], (size_t)n_work_alloc, sizeof(double), GMT_program);
}


void	free_space (void)
{
	GMT_free ((void *)data);
	GMT_free ((void *)wt_sum);
	GMT_free ((void *)data_sum);
	GMT_free ((void *)n_this_col);
	GMT_free ((void *)good_one);
	GMT_free ((void *)n_left);
	GMT_free ((void *)n_right);
	GMT_free ((void *)work);
	GMT_free ((void *)min_loc);
	GMT_free ((void *)max_loc);
	GMT_free ((void *)last_loc);
	GMT_free ((void *)this_loc);
	GMT_free ((void *)min_scl);
	GMT_free ((void *)max_scl);
	GMT_free ((void *)last_scl);
	GMT_free ((void *)this_scl);
}

BOOLEAN	lack_check (int i_col, int left, int right)
{			

	int	last_row, this_row;
	BOOLEAN	lacking = FALSE;
	double	last_t;
	
	last_row = left;
	while (!(GMT_is_dnan (data[i_col][last_row])) && last_row < (right - 1)) last_row++;
	
	last_t = data[t_col][last_row];
	this_row = last_row + 1;
	while ( !(lacking) && this_row < (right - 1)) {
		while (!(GMT_is_dnan (data[i_col][this_row])) && this_row < (right - 1)) this_row++;
		
		if ( (data[t_col][this_row] - last_t) > lack_width)
			lacking = TRUE;
		else {
			last_t = data[t_col][this_row];
			last_row = this_row;
			this_row++;
		}
	}
	return (lacking);
}

void	get_robust_estimates (int j, int n, int both)
{
	int i, n_smooth, sort_me = TRUE;
	double	low, high, last, temp;
	
	if (filter_type == MODE) {

		n_smooth = n / 2;
	
		GMT_mode (work[j], n, n_smooth, sort_me, GMT_mode_selection, &GMT_n_multiples, &temp);
	}
	else if (filter_type == MEDIAN) {

		low = min_loc[j];
		high = max_loc[j];
		last = last_loc[j];
	
		GMT_median (work[j], n, low, high, last, &temp);
	}
	else {	/* EXTREME filter **/
		temp = GMT_extreme (work[j], n, extreme, kind, way);
	}

	last_loc[j] = this_loc[j] = temp;
	
	if (both) {
		for (i = 0; i < n; i++)
			work[j][i] = fabs(work[j][i] - this_loc[j]);
		low = min_scl[j];
		high = max_scl[j];
		last = last_scl[j];
		GMT_median(work[j], n, low, high, last, &temp);
		last_scl[j] = this_scl[j] = temp;
	} 
}
	
double	boxcar_weight (double radius, double half_width)
{
	double weight;
	
	if (radius > half_width)
		weight = 0.0;
	else
		weight = 1.0;
		
	return (weight);
}

double	cosine_weight (double radius, double half_width)
{
	double weight, cosine_constant;
	cosine_constant = M_PI / half_width;
	if (radius > half_width)
		weight = 0.0;
	else
		weight = 1.0 + cos (radius * cosine_constant);
		
	return (weight);
}

double	gaussian_weight (double radius, double half_width)
{
	double weight, gauss_constant;
	gauss_constant = -4.5 / (half_width * half_width);
	if (radius > half_width)
		weight = 0.0;
	else
		weight = exp (radius * radius * gauss_constant);
		
	return (weight);
}
