/*--------------------------------------------------------------------
 *	$Id: pshistogram.c,v 1.4.4.2 2002/02/27 17:41:10 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
 *--------------------------------------------------------------------*/
/*
 * pshistogram.c -- a program for plotting histograms
 *
 *
 * Author:	Walter H. F. Smith
 * Date:	15 February, 1988
 *
 * Updated to v2.0 5-May-1991 Paul Wessel
 * Updated to v3.0 1-Jan-1995 Paul Wessel
 * Updated to v3.1 13-Jun-1998 Paul Wessel
 * Updated to v3.3 23-Apr-1999 Paul Wessel
 *		     Now uses two passes to get accurate mean and stdev
 * Version:	3.4.1
 */

#include "gmt.h"

#define	COUNTS 0
#define FREQ_PCT 1
#define LOG_COUNTS 2
#define LOG_FREQ_PCT 3
#define LOG10_COUNTS 4
#define LOG10_FREQ_PCT 5

double	yy0, yy1, xmin, xmax;

float	*x;
int	*boxh;

int	n = 0;
int	n_boxes = 0;
int	n_counted = 0;
double	box_width = 0.0;
double	west = 0.0, east = 0.0, south = 0.0, north = 0.0;
float	stats[6];

FILE	*fp_in = NULL;

char	buffer[BUFSIZ], format[BUFSIZ];

int	center_box = FALSE, cumulative = FALSE, draw_outline = FALSE;
int	second = FALSE, hist_type = COUNTS;

struct GMT_PEN pen;
struct GMT_FILL fill;
	
int read_data(void);
int fill_boxes(void);
int plot_boxes(BOOLEAN stairs, BOOLEAN flip_to_y);
int get_loc_scl(float *x, int n, float *stats);

main (int argc, char **argv)
{
	int i;
	BOOLEAN inquire = FALSE, error = FALSE, automatic = FALSE;
	int fill_bar = FALSE, stairs = FALSE, flip_to_y = FALSE, dump_boxes = FALSE;
	double tmp;
	char string[BUFSIZ], *text;
	
	argc = GMT_begin (argc, argv);

	GMT_init_pen (&pen, GMT_PENWIDTH);
	GMT_init_fill (&fill, -1, 0, 0);	/* Do not fill is default */
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
		
				/* Common parameters */
			
				case 'B':
				case 'H':
				case 'J':
				case 'K':
				case 'O':
				case 'P':
				case 'U':
				case 'V':
				case 'X':
				case 'x':
				case 'Y':
				case 'y':
				case 'c':
				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 '2':
					second = TRUE;
					break;
				case 'A':
					flip_to_y = TRUE;
					break;
				case 'C':
					center_box = TRUE;
					break;
				case 'E':
					sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
					break;
				case 'G':
					if (GMT_getfill (&argv[i][2], &fill)) {
						GMT_fill_syntax ('G');
						error++;
					}
					fill_bar = TRUE;
					break;
				case 'L':		/* Set line attributes */
					if (GMT_getpen (&argv[i][2], &pen)) {
						GMT_pen_syntax ('L');
						error++;
					}
					draw_outline = TRUE;
					break;
				case 'I':
					inquire = TRUE;
					if (argv[i][2] == 'o') dump_boxes = TRUE;
					break;
				case 'Q':
					cumulative = TRUE;
					break;
				case 'S':
					stairs = TRUE;
					break;
				case 'R':
					project_info.region_supplied = TRUE;

					strcpy (string, &argv[i][2]);
					text = strtok (string, "/");
					west = GMT_ddmmss_to_degree (text);
					text = strtok (CNULL, "/");
					east = GMT_ddmmss_to_degree (text);
					text = strtok (CNULL, "/");
					south = GMT_ddmmss_to_degree (text);
					text = strtok (CNULL, "/");
					north = GMT_ddmmss_to_degree (text);
					break;
				case 'W':
					box_width = atof (&argv[i][2]);
					break;
				case 'Z':
					hist_type = atoi (&argv[i][2]);
					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,"pshistogram %s - Calculate and plot histograms\n\n", GMT_VERSION);
		fprintf (stderr, "usage: pshistogram [file] -Jx|X<scale> -W<width> [-2] [-B<tickinfo>] [-C] [-Eaz/el] [-G<fill>]\n");
		fprintf (stderr, "\t[-H[<nrec>]] [-I[o]] [-K] [-L<pen>] [-O] [-P] [-Q] [-Rw/e/s/n] [-S] [-U[text]] [-V]\n");
		fprintf (stderr, "\t[-X<x_shift>] [-Y<y_shift>] [-Z[type]] [-c<ncopies>] [-bi[s][<n>]]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		fprintf (stderr, "\t-Jx|X for linear projection.  Scale in %s/units (or width in %s)\n", GMT_unit_names[gmtdefs.measure_unit], GMT_unit_names[gmtdefs.measure_unit]);
		fprintf (stderr, "\t    Use / to specify separate x/y scaling.\n");
		fprintf (stderr, "\t    If -JX is used then give axes lengths in %s rather than scales\n", GMT_unit_names[gmtdefs.measure_unit]);

		fprintf (stderr, "\t-W sets the bin width\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-2 use second rather than first column\n");
		GMT_explain_option ('b');
		fprintf (stderr, "\t-A will plot horizontal bars [Default is vertical]\n");
		fprintf (stderr, "\t-C will center the bins\n");
		fprintf (stderr, "\t-E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
		fprintf (stderr, "\t-G Specify color or pattern for columns. fill can be either\n");
		fprintf (stderr, "\t   1) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0]\n");
		fprintf (stderr, "\t   2) p[or P]<dpi>/<pattern> for predefined patterns (1-90).\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I will inquire about min/max x and y.  No plotting is done\n");
		fprintf (stderr, "\t   Append o to output the resulting x, y data\n");
		GMT_explain_option ('K');
		fprintf (stderr, "\t-L <pen> is penwidth. Specify r/g/b for colored line.\n");
		GMT_explain_option ('O');
		GMT_explain_option ('P');
		fprintf (stderr, "\t-Q plot a cumulative histogram\n");
		GMT_explain_option ('r');
		fprintf (stderr, "\t   If neither -R nor -I are set, w/e/s/n will be based on input data\n");
		fprintf (stderr, "\t-S Draws a stairs-step diagram instead of histogram.\n");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		GMT_explain_option ('X');
		fprintf (stderr, "\t-Z to choose type of vertical axis.  Select from\n");
		fprintf (stderr, "\t   0 - Counts [Default]\n");
		fprintf (stderr, "\t   1 - Frequency percent\n");
		fprintf (stderr, "\t   2 - Log (1+counts)\n");
		fprintf (stderr, "\t   3 - Log (1+frequency percent)\n");
		fprintf (stderr, "\t   4 - Log10 (1+counts)\n");
		fprintf (stderr, "\t   5 - Log10 (1+frequency percent)\n");
		GMT_explain_option ('c');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 2 input columns\n");
		GMT_explain_option ('.');

		exit (EXIT_FAILURE);
	}
	
	if (!inquire && (project_info.projection < 0 || project_info.projection > 9)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -J option:  Only linear projection supported.\n", GMT_program);
		error++;
	}
	if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", GMT_program);
		error++;
	}
	if (hist_type < COUNTS || hist_type > LOG10_FREQ_PCT) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  histogram type must be in 0-5 range\n", GMT_program);
		error++;
	}
	if (box_width <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -W option:  bin width must be nonzero\n", GMT_program);
		error++;
	}
	/* Now must specify either fill color with -G or outline pen with -L; JLL */
 
	if (!(inquire || fill_bar || draw_outline)) {
		fprintf (stderr, "%s: Must specify either fill color (-G), outline pen attributes (-L), or both.\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++;
	}
        if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 2;
        if (GMT_io.binary[0] && GMT_io.ncol[0] < 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 .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]]);
	}

	if (inquire) gmtdefs.verbose = TRUE;
	if (!inquire && !project_info.region_supplied) automatic = TRUE;
	
	if (fp_in == NULL) {
		fp_in = GMT_stdin;
#ifdef SET_IO_MODE
		GMT_setmode (0);
#endif
	}

	if ( read_data() ) {
		fprintf (stderr, "%s: Fatal error, read only 0 points.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	
	if (gmtdefs.verbose) {
		sprintf (format, "%%s: Extreme values of the data :\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, format, GMT_program, x[0], x[n-1]);
		sprintf (format, "%%s: Locations: L2, L1, LMS; Scales: L2, L1, LMS\t%s\t%s\t%s\t%s\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, format, GMT_program, stats[0], stats[1], stats[2], stats[3], stats[4], stats[5]);
	}
	
	if (west == east) {	/* Set automatic x range [ and tickmarks] */
		if (frame_info.anot_int[0] == 0.0) {
			tmp = pow (10.0, floor (d_log10 (xmax-xmin)));
			if (((xmax-xmin) / tmp) < 3.0) tmp *= 0.5;
		}
		else
			tmp = frame_info.anot_int[0];
		west = floor (xmin / tmp) * tmp;
		east = ceil (xmax / tmp) * tmp;
		if (frame_info.anot_int[0] == 0.0) {
			frame_info.anot_int[0] = frame_info.frame_int[0] = tmp;
			frame_info.plot = TRUE;
		}
	}
	
	if ( fill_boxes () ) {
		fprintf (stderr, "%s: Fatal error during box fill.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	
	if (gmtdefs.verbose) {
		sprintf (format, "\n%%s: min/max values are :\t%s\t%s\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, format, GMT_program, xmin, xmax, yy0, yy1);
	}
	
	if (inquire) {	/* Only info requested, quit before plotting */
		if (dump_boxes) {
			int ibox;
			double xx, yy;
			sprintf (format, "%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format);
			for (ibox = 0; ibox < n_boxes; ibox++) {
				if (boxh[ibox] == 0) continue;
				xx = west + ibox * box_width;
				if (center_box) xx -= (0.5 * box_width);
				if (hist_type == LOG_COUNTS)
					yy = d_log1p( (double)boxh[ibox]);
				else if (hist_type == LOG10_COUNTS)
					yy = d_log101p( (double)boxh[ibox]);
				else if (hist_type == FREQ_PCT)
					yy = (100.0 * boxh[ibox]) / n_counted;
				else if (hist_type == LOG_FREQ_PCT)
					yy = d_log1p( 100.0 * boxh[ibox] / n_counted );
				else if (hist_type == LOG10_FREQ_PCT)
					yy = d_log101p( 100.0 * boxh[ibox] / n_counted );
				else
					yy = boxh[ibox];
				fprintf (GMT_stdout, format, xx, yy);
			}
		}
		GMT_free ((void *) x);
		GMT_free ((void *) boxh);
		exit (EXIT_SUCCESS);
	}
	
	if (automatic) {	/* Set up s/n based on 'clever' rounding up of the minmax values */
		project_info.region = TRUE;
		south = 0.0;
		if (frame_info.anot_int[1] == 0.0) {
			tmp = pow (10.0, floor (d_log10 (yy1)));
			if ((yy1 / tmp) < 3.0) tmp *= 0.5;
		}
		else
			tmp = frame_info.anot_int[1];
		north = ceil (yy1 / tmp) * tmp;
		if (frame_info.anot_int[1] == 0.0) {	/* Tickmarks not set */
			frame_info.anot_int[1] = frame_info.frame_int[1] = tmp;
			frame_info.plot = TRUE;
		}
		if (project_info.pars[0] == 0.0) {	/* Must give default xscale */
			project_info.pars[0] = gmtdefs.x_axis_length / (east - west);
			project_info.projection = LINEAR;
		}
		if (project_info.pars[1] == 0.0) {	/* Must give default yscale */
			project_info.pars[1] = gmtdefs.y_axis_length / north;
			project_info.projection = LINEAR;
		}
	}
	
	if (automatic && gmtdefs.verbose) {
		sprintf (format, "%%s: Use w/e/s/n = %s/%s/%s/%s and x-tick/y-tick = %s/%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, format, GMT_program, west, east, south, north, frame_info.anot_int[0], frame_info.anot_int[1]);
	}

	if (flip_to_y) {
		char buffer[128];
		d_swap (frame_info.anot_int[0], frame_info.anot_int[1]);
		d_swap (frame_info.frame_int[0], frame_info.frame_int[1]);
		d_swap (frame_info.grid_int[0], frame_info.grid_int[1]);
		text = frame_info.label[0];
		strcpy (buffer, frame_info.label[0]);
		strcpy(frame_info.label[0], frame_info.label[1]);
		strcpy (frame_info.label[1], buffer);
		GMT_map_setup (south, north, west, east);
	}
	else
		GMT_map_setup (west, east, south, north);

	ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
		gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
		gmtdefs.dpi, GMT_INCH, gmtdefs.paper_width, gmtdefs.page_rgb, GMT_epsinfo (argv[0]));
	GMT_echo_command (argc, argv);
	if (gmtdefs.unix_time) GMT_timestamp (argc, argv);
	
	if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);

	GMT_map_clip_on (GMT_no_rgb, 3);
	if ( plot_boxes (stairs, flip_to_y) ) {
		fprintf (stderr, "%s: Fatal error during box plotting.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	
	GMT_map_clip_off();
	
	if (frame_info.plot) GMT_map_basemap ();

	if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
	ps_plotend (gmtdefs.last_page);
	
	GMT_free ((void *) x);
	GMT_free ((void *) boxh);
	
	GMT_end (argc, argv);
}

int	read_data (void) {
	int i, n_alloc = GMT_CHUNK, n_expected_fields, n_fields;
	double *in;
	
	x = (float *) GMT_memory (VNULL, (size_t)n_alloc , sizeof (float), GMT_program);
	
	n = 0;
	xmin = DBL_MAX;	xmax = -DBL_MAX;
	n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 1 + second;
	
	if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (buffer, BUFSIZ, fp_in);
	
	while ((n_fields = GMT_input (fp_in, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_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, 2, n);
			exit (EXIT_FAILURE);
		}

		x[n] = (float)((second) ? in[1] : in[0]);
		if (!GMT_is_fnan (x[n])) n++;

		if (n == n_alloc) {
			n_alloc += GMT_CHUNK;
			x = (float *) GMT_memory ((void *)x, (size_t) n_alloc, sizeof (float), GMT_program);
		}
		xmin = MIN (xmin, x[n-1]);
		xmax = MAX (xmax, x[n-1]);
	}
	
	if (gmtdefs.verbose) {
		fprintf (stderr,"%s: %d points read\n", GMT_program, n);
	}

	x = (float *) GMT_memory ((void *)x, (size_t)n , sizeof (float), GMT_program);
	if (fp_in != GMT_stdin) GMT_fclose (fp_in);

	get_loc_scl (x, n, stats);
	
	return( !(n) );
}

int	fill_boxes (void) {

	double	add_half = 0.0;
	
	int	b0, b1, i, ibox, count_sum;

	n_boxes = (int)ceil( (east - west) / box_width);
	
	if (center_box) {
		n_boxes++;
		add_half = 0.5;
	}
	
	if (n_boxes <= 0) return (-1);

	boxh = (int *) GMT_memory (VNULL, (size_t)n_boxes , sizeof (int), GMT_program);

	n_counted = 0;
	
	/* First fill boxes with counts  */
	
	for (i = 0; i < n; i++) {
		ibox = (int)floor( ( (x[i] - west) / box_width) + add_half);
		if (ibox < 0 || ibox >= n_boxes) continue;
		boxh[ibox] ++;
		n_counted++;
	}

	if (cumulative) {
		count_sum = 0;
		for (ibox = 0; ibox < n_boxes; ibox++) {
			count_sum += boxh[ibox];
			boxh[ibox] = count_sum;
		}
		b0 = 0;
		b1 = count_sum;
	}
	else {
		b0 = n_counted;
		b1 = 0;
		for (ibox = 0; ibox < n_boxes; ibox++) {
			if (b0 > boxh[ibox]) b0 = boxh[ibox];
			if (b1 < boxh[ibox]) b1 = boxh[ibox];
		}
	}

	/* Now find out what the min max y will be  */
	
	if (b0) {
		if (hist_type == LOG_COUNTS)
			yy0 = d_log1p((double)b0);
		else if (hist_type == LOG10_COUNTS)
			yy0 = d_log101p( (double)b0);
		else if (hist_type == FREQ_PCT)
			yy0 = (100.0 * b0) / n_counted;
		else if (hist_type == LOG_FREQ_PCT)
			yy0 = d_log1p( 100.0 * b0 / n_counted );
		else if (hist_type == LOG10_FREQ_PCT)
			yy0 = d_log101p( 100.0 * b0 / n_counted );
		else
			yy0 = b0;
	}
	else {
		yy0 = 0;
	}
	if (b1) {
		if (hist_type == LOG_COUNTS)
			yy1 = d_log1p( (double)b1);
		else if (hist_type == LOG10_COUNTS)
			yy1 = d_log101p( (double)b1);
		else if (hist_type == FREQ_PCT)
			yy1 = (100.0 * b1) / n_counted;
		else if (hist_type == LOG_FREQ_PCT)
			yy1 = d_log1p( 100.0 * b1 / n_counted );
		else if (hist_type == LOG10_FREQ_PCT)
			yy1 = d_log101p( 100.0 * b1 / n_counted );
		else
			yy1 = b1;
	}
	else {
		yy1 = 0;
	}
	return (0);
}

int	plot_boxes (BOOLEAN stairs, BOOLEAN flip_to_y)
{

	int	i, ibox, first = TRUE;
	
	double	x[4], y[4], xx, yy, *px, *py;
	
	if (draw_outline) GMT_setpen (&pen);
	
	if (flip_to_y) {
		px = y;
		py = x;
	}
	else {
		px = x;
		py = y;
	}

	for (ibox = 0; ibox < n_boxes; ibox++) {
		if (stairs || boxh[ibox]) {
			x[0] = west + ibox * box_width;
			if (center_box) x[0] -= (0.5 * box_width);
			x[1] = x[0] + box_width;
			if (x[0] < west) x[0] = west;
			if (x[1] > east) x[1] = east;
			x[2] = x[1];
			x[3] = x[0];
			y[0] = y[1] = south;
			if (hist_type == LOG_COUNTS)
				y[2] = d_log1p( (double)boxh[ibox]);
			else if (hist_type == LOG10_COUNTS)
				y[2] = d_log101p( (double)boxh[ibox]);
			else if (hist_type == FREQ_PCT)
				y[2] = (100.0 * boxh[ibox]) / n_counted;
			else if (hist_type == LOG_FREQ_PCT)
				y[2] = d_log1p( 100.0 * boxh[ibox] / n_counted );
			else if (hist_type == LOG10_FREQ_PCT)
				y[2] = d_log101p( 100.0 * boxh[ibox] / n_counted );
			else
				y[2] = boxh[ibox];
			y[3] = y[2];
			
			for (i = 0; i < 4; i++) {
				GMT_geo_to_xy (px[i], py[i], &xx, &yy);
				if (project_info.three_D) GMT_xyz_to_xy (xx, yy, project_info.z_level, &xx, &yy);
				px[i] = xx;	py[i] = yy;
			}
			
			if (stairs) {
				if (first) {
					first = FALSE;
					ps_plot (px[0], py[0], 3);
				}
				ps_plot (px[3], py[3], 2);
				ps_plot (px[2], py[2], 2);
			}
			else
				GMT_fill (px, py, 4, &fill, draw_outline);
		}
	}

	if (stairs) ps_plot (px[1], py[1], 2);
	
	return(0);
}

int	get_loc_scl (float *x, int n, float *stats)
{
	/* Returns stats[] = L2, L1, LMS location, L2, L1, LMS scale  */

	int	i, j;
	double	xsum, x2sum, dx, median, mode, mad, lmsscl;

	if (n < 3) return(-1);
	
	qsort ((void *)x, (size_t)n, sizeof(float), GMT_comp_float_asc);

	/* Get median */
	j = n/2;
	median = (n%2) ? x[j] : (0.5 * (x[j] + x[j-1]));
	stats[1] = (float)median;

	/* Get mode */
	GMT_mode_f (x, n, j, 0, &mode);
	stats[2] = (float)mode;
	
	/* Get MAD for L1 */
	GMT_getmad_f (x, n, median, &mad);
	stats[4] = (float)mad;
	
	/* Get LMSscale for mode */
	GMT_getmad_f (x, n, mode, &lmsscl);
	stats[5] = (float)lmsscl;

	/* Calculate mean and stdev in two passes to minimize risk of overflow */

	xsum = x2sum = 0.0;
	for (i = 0; i < n; i++) xsum += x[i];
	xsum /= n;	/* xsum is now the mean value */
	for (i = 0; i < n; i++) {
		dx = x[i] - xsum;
		x2sum += (dx * dx);
	}
	stats[0] = (float)xsum;
	stats[3] = (float)sqrt (x2sum / (n - 1));

	return (0);
}