/* TableOfReal_extensions.c
 *
 * Copyright (C) 1993-2003 David Weenink
 *
 * 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; either version 2 of the License, or (at
 * your option) any later version.
 *
 * 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.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 djmw 20000202 17 typos in F1-3,L1-3 table corrected 
 djmw 20030512 Latest modification
*/

#include <ctype.h>
#include "NUMclapack.h"
#include "NUM2.h"
#include "SVD.h"
#include "Sequence.h"
#include "TableOfReal_extensions.h"
#include "regularExp.h"
#include "Formula.h"

#define EMPTY_STRING(s) ((s) == NULL || s[0] == '\0')
#define MAX(m,n) ((m) > (n) ? (m) : (n))
#define MIN(m,n) ((m) < (n) ? (m) : (n))

#define Graphics_ARROW 1
#define Graphics_TWOWAYARROW 2
#define Graphics_LINE 3


int TableOfReal_hasRowLabels (I)
{
	iam (TableOfReal);
	long i;
	
	if (my rowLabels == NULL) return 0;
	for (i = 1; i <= my numberOfRows; i++)
	{
		if (EMPTY_STRING(my rowLabels[i])) return 0;
	}
	return 1;
}

int TableOfReal_hasColumnLabels (I)
{
	iam (TableOfReal);
	long i;
	
	if (my columnLabels == NULL) return 0;
	for (i = 1; i <= my numberOfColumns; i++)
	{
		if (EMPTY_STRING (my columnLabels[i])) return 0;
	}
	return 1;
}

TableOfReal TableOfReal_createIrisDataset (void)
{
	TableOfReal me = NULL;
	long i, j; 
	float iris[150][4] = {
	5.1,3.5,1.4,0.2,4.9,3.0,1.4,0.2,4.7,3.2,1.3,0.2,4.6,3.1,1.5,0.2,5.0,3.6,1.4,0.2,
	5.4,3.9,1.7,0.4,4.6,3.4,1.4,0.3,5.0,3.4,1.5,0.2,4.4,2.9,1.4,0.2,4.9,3.1,1.5,0.1,
	5.4,3.7,1.5,0.2,4.8,3.4,1.6,0.2,4.8,3.0,1.4,0.1,4.3,3.0,1.1,0.1,5.8,4.0,1.2,0.2,
	5.7,4.4,1.5,0.4,5.4,3.9,1.3,0.4,5.1,3.5,1.4,0.3,5.7,3.8,1.7,0.3,5.1,3.8,1.5,0.3,
	5.4,3.4,1.7,0.2,5.1,3.7,1.5,0.4,4.6,3.6,1.0,0.2,5.1,3.3,1.7,0.5,4.8,3.4,1.9,0.2,
	5.0,3.0,1.6,0.2,5.0,3.4,1.6,0.4,5.2,3.5,1.5,0.2,5.2,3.4,1.4,0.2,4.7,3.2,1.6,0.2,
	4.8,3.1,1.6,0.2,5.4,3.4,1.5,0.4,5.2,4.1,1.5,0.1,5.5,4.2,1.4,0.2,4.9,3.1,1.5,0.2,
	5.0,3.2,1.2,0.2,5.5,3.5,1.3,0.2,4.9,3.6,1.4,0.1,4.4,3.0,1.3,0.2,5.1,3.4,1.5,0.2,
	5.0,3.5,1.3,0.3,4.5,2.3,1.3,0.3,4.4,3.2,1.3,0.2,5.0,3.5,1.6,0.6,5.1,3.8,1.9,0.4,
	4.8,3.0,1.4,0.3,5.1,3.8,1.6,0.2,4.6,3.2,1.4,0.2,5.3,3.7,1.5,0.2,5.0,3.3,1.4,0.2,
	7.0,3.2,4.7,1.4,6.4,3.2,4.5,1.5,6.9,3.1,4.9,1.5,5.5,2.3,4.0,1.3,6.5,2.8,4.6,1.5,
	5.7,2.8,4.5,1.3,6.3,3.3,4.7,1.6,4.9,2.4,3.3,1.0,6.6,2.9,4.6,1.3,5.2,2.7,3.9,1.4,
	5.0,2.0,3.5,1.0,5.9,3.0,4.2,1.5,6.0,2.2,4.0,1.0,6.1,2.9,4.7,1.4,5.6,2.9,3.6,1.3,
	6.7,3.1,4.4,1.4,5.6,3.0,4.5,1.5,5.8,2.7,4.1,1.0,6.2,2.2,4.5,1.5,5.6,2.5,3.9,1.1,
	5.9,3.2,4.8,1.8,6.1,2.8,4.0,1.3,6.3,2.5,4.9,1.5,6.1,2.8,4.7,1.2,6.4,2.9,4.3,1.3,
	6.6,3.0,4.4,1.4,6.8,2.8,4.8,1.4,6.7,3.0,5.0,1.7,6.0,2.9,4.5,1.5,5.7,2.6,3.5,1.0,
	5.5,2.4,3.8,1.1,5.5,2.4,3.7,1.0,5.8,2.7,3.9,1.2,6.0,2.7,5.1,1.6,5.4,3.0,4.5,1.5,
	6.0,3.4,4.5,1.6,6.7,3.1,4.7,1.5,6.3,2.3,4.4,1.3,5.6,3.0,4.1,1.3,5.5,2.5,4.0,1.3,
	5.5,2.6,4.4,1.2,6.1,3.0,4.6,1.4,5.8,2.6,4.0,1.2,5.0,2.3,3.3,1.0,5.6,2.7,4.2,1.3,
	5.7,3.0,4.2,1.2,5.7,2.9,4.2,1.3,6.2,2.9,4.3,1.3,5.1,2.5,3.0,1.1,5.7,2.8,4.1,1.3,
	6.3,3.3,6.0,2.5,5.8,2.7,5.1,1.9,7.1,3.0,5.9,2.1,6.3,2.9,5.6,1.8,6.5,3.0,5.8,2.2,
	7.6,3.0,6.6,2.1,4.9,2.5,4.5,1.7,7.3,2.9,6.3,1.8,6.7,2.5,5.8,1.8,7.2,3.6,6.1,2.5,
	6.5,3.2,5.1,2.0,6.4,2.7,5.3,1.9,6.8,3.0,5.5,2.1,5.7,2.5,5.0,2.0,5.8,2.8,5.1,2.4,
	6.4,3.2,5.3,2.3,6.5,3.0,5.5,1.8,7.7,3.8,6.7,2.2,7.7,2.6,6.9,2.3,6.0,2.2,5.0,1.5,
	6.9,3.2,5.7,2.3,5.6,2.8,4.9,2.0,7.7,2.8,6.7,2.0,6.3,2.7,4.9,1.8,6.7,3.3,5.7,2.1,
	7.2,3.2,6.0,1.8,6.2,2.8,4.8,1.8,6.1,3.0,4.9,1.8,6.4,2.8,5.6,2.1,7.2,3.0,5.8,1.6,
	7.4,2.8,6.1,1.9,7.9,3.8,6.4,2.0,6.4,2.8,5.6,2.2,6.3,2.8,5.1,1.5,6.1,2.6,5.6,1.4,
	7.7,3.0,6.1,2.3,6.3,3.4,5.6,2.4,6.4,3.1,5.5,1.8,6.0,3.0,4.8,1.8,6.9,3.1,5.4,2.1,
	6.7,3.1,5.6,2.4,6.9,3.1,5.1,2.3,5.8,2.7,5.1,1.9,6.8,3.2,5.9,2.3,6.7,3.3,5.7,2.5,
	6.7,3.0,5.2,2.3,6.3,2.5,5.0,1.9,6.5,3.0,5.2,2.0,6.2,3.4,5.4,2.3,5.9,3.0,5.1,1.8	
	};

	if (! (me = TableOfReal_create (150, 4))) return NULL;
	
	TableOfReal_setColumnLabel (me, 1, "sl");
	TableOfReal_setColumnLabel (me, 2, "sw");
	TableOfReal_setColumnLabel (me, 3, "pl");
	TableOfReal_setColumnLabel (me, 4, "pw");
	for (i=1; i <= 150; i++)
	{
		int kind = (i - 1) / 50 + 1;
		char *label = kind == 1 ? "1" : kind == 2 ? "2" : "3";
		for (j=1; j <= 4; j++) my data[i][j] = iris[i-1][j-1];
		TableOfReal_setRowLabel (me, i, label);
	}
	Thing_setName (me, "iris");
	return me;
}

Strings TableOfReal_extractRowLabels (I)
{
	iam (TableOfReal);
	Strings thee = new (Strings);
	long i, n = my numberOfRows;
	
	if (thee == NULL) return NULL;
	if (n < 1) return thee;
	
	thy strings = NUMpvector (1, n);
	if (thy strings == NULL) goto end;
	thy numberOfStrings = n;
	
	for (i = 1; i <= n; i++)
	{
		char *label = my rowLabels[i] ? my rowLabels[i] : "?";
		thy strings[i] = Melder_strdup (label);
		if (thy strings[i] == NULL) goto end; 
	}

end:

	if (Melder_hasError()) forget (thee);
	return thee;	
}

Strings TableOfReal_extractColumnLabels (I)
{
	iam (TableOfReal);
	Strings thee = new (Strings);
	long i, n = my numberOfColumns;
		
	if (thee == NULL) return NULL;
	if (n < 1) return thee;
	
	thy strings = NUMpvector (1, n);
	if (thy strings == NULL) goto end;
	thy numberOfStrings = n;
	
	for (i = 1; i <= n; i++)
	{
		char *label = my columnLabels[i] ? my columnLabels[i] : "?";
		thy strings[i] = Melder_strdup (label);
		if (thy strings[i] == NULL) goto end; 
	}

end:

	if (Melder_hasError()) forget (thee);
	return thee;	
}

TableOfReal TableOfReal_transpose (I)
{
	iam (TableOfReal);
	TableOfReal thee;
	long i, j;
	
	thee = TableOfReal_create (my numberOfColumns, my numberOfRows);
	if (thee == NULL) return NULL;

	for (i = 1; i <= my numberOfRows; i++)
	{
		for (j = 1; j <= my numberOfColumns; j++)
		{
			thy data[j][i] = my data[i][j];
		}
	}
		
	if (! NUMstrings_copyElements (my rowLabels, thy columnLabels, 
		1, my numberOfRows) ||
		! NUMstrings_copyElements (my columnLabels, thy rowLabels, 
		1, my numberOfColumns)) forget (thee);
	return thee;
}

int TableOfReal_split (I, long fromrow, long torow, long fromcol, long tocol,
	Pattern *p, Categories *c)
{
	iam (TableOfReal);
	long i, j, ncol = my numberOfColumns, nrow = my numberOfRows, row, col;
	int status = 1;

	if (fromrow == torow && fromrow == 0)
	{
		fromrow = 1; torow = nrow;
	}
	else if (fromrow > 0 && fromrow <= nrow && torow == 0)
	{
		torow = nrow;
	}
	else if (! (fromrow > 0 && torow <= nrow && fromrow <= torow))
	{
		return Melder_error ("TableOfReal_split: illegal row selection.");
	}
	if (fromcol == tocol && fromcol == 0)
	{
		fromcol = 1; tocol = ncol;
	}
	else if (fromcol > 0 && fromcol <= ncol && tocol == 0)
	{
		tocol = ncol;
	}
	else if (! (fromcol > 0 && tocol <= ncol && fromcol <= tocol))
	{
		return Melder_error ("TableOfReal_split: illegal col selection.");
	}
	nrow = torow - fromrow + 1;
	ncol = tocol - fromcol + 1;
	
	*c = NULL;
	if (! (*p = Pattern_create (nrow, ncol)) ||
		! (*c = Categories_create ())) goto end;
	
	for (row=1, i=fromrow; i <= torow; i++, row++)
	{
		char *s = my rowLabels[i] ? my rowLabels[i] : "?";
		SimpleString item = SimpleString_create (s);
		if (! item || ! Collection_addItem (*c, item)) goto end;
		for (col=1, j=fromcol; j <= tocol; j++, col++)
		{
			(*p)->z[row][col] = my data[i][j];
		}
	}
	
end:

	if (Melder_hasError ())
	{
		forget (*p);
		forget (*c);
		status = 0;
	}
	return status;
}

void TableOfReal_getColumnExtrema (I, long col, double *min, double *max)
{
	iam (TableOfReal); long i;
	if (col < 1 || col > my numberOfColumns)
	{
		(void) Melder_error ("TableOfReal_getColumnExtrema: not a valid column.");
		*min = NUMundefined; *max = NUMundefined; return;
	}
	*min = *max = my data[1][col];
	for (i=2; i <= my numberOfRows; i++)
	{
		if (my data[i][col] > *max) *max = my data[i][col];
		else if (my data[i][col] < *min) *min = my data[i][col];
	}
}

void TableOfReal_drawRowsAsHistogram (I, Graphics g, char *rows, long colb, long cole,
	double ymin, double ymax, double xoffsetFraction, double interbarFraction, double grey,
	double interbarsFraction, int labels, int garnish)	
{
	iam (TableOfReal);
	char *proc = "TableOfReal_drawRowsAsHistogram";
	long i, j, irow, nrows, ncols;
	double *irows, bar_width = 1, xb, dx, x1, x2, y1, y2;
	
	if (cole >= colb)
	{
		colb = 1; cole = my numberOfColumns;
	}
	if (colb <= cole && colb < 1 || cole > my numberOfColumns)
	{
		    Melder_warning ("%s: Invalid columns", proc);
			return;
	}
	
	irows = NUMstring_to_numbers (rows, &nrows);
	if (irows == NULL)
	{
		 Melder_warning ("%s: No rows!", proc);
		 return;
	}
	
	for (i = 1; i <= nrows; i++)
	{
		irow = irows[i];
		if (irow < 0 || irow > my numberOfRows)
		{
			 Melder_warning ("%s: Invalid row (%d)", proc, irow);
			 goto end;
		}
		if (ymin >= ymax)
		{
			double min, max;
			NUMdvector_extrema (my data[irow], colb, cole, &min, &max);
			if (i > 1)
			{
				if (min < ymin) ymin = min;
				if (max > ymax) ymax = max;
			}
			else
			{
				ymin = min; ymax = max;
			}
		}
	}
    
    Graphics_setWindow (g, 0, 1, ymin, ymax);
    Graphics_setInner (g);
	
	ncols = cole - colb + 1;	
	bar_width /= ncols * nrows + 2 * xoffsetFraction + (ncols - 1) * interbarsFraction + 
		ncols * (nrows -1) * interbarFraction;
    dx = (interbarsFraction + nrows + (nrows -1) * interbarFraction) * bar_width;

	for (i = 1; i <= nrows; i++)
	{
		irow = irows[i];
		xb = xoffsetFraction * bar_width + (i - 1) * (1 + interbarFraction) * bar_width;  
	
		x1 = xb;
    	for (j = colb; j <= cole; j++)
    	{
			x2 = x1 + bar_width;
			y1 = ymin; y2 = my data[irow][j];
			if (y2 > ymin)
			{
				if (y2 > ymax) y2 = ymax;
				Graphics_setGrey (g, grey);
				Graphics_fillRectangle (g, x1, x2, y1, y2);
				Graphics_setGrey (g, 0); /* Black */
				Graphics_rectangle (g, x1, x2, y1, y2);
			}
			x1 += dx;
    	}
	}
	
    Graphics_unsetInner (g);
	
	if (labels)
	{
		xb = (xoffsetFraction + 0.5 * (nrows + (nrows - 1) * interbarFraction)) * bar_width;
		for (j = colb; j <= cole; j++)
		{
			if (my columnLabels[j]) Graphics_markBottom (g, xb, 0, 0, 0, my columnLabels[j]);
			xb += dx;
		}
	}
	if (garnish)
	{
		Graphics_drawInnerBox (g);
		Graphics_marksLeft (g, 2, 1, 1, 0);
	}

end:
	
	NUMdvector_free (irows, 1);
}

void TableOfReal_drawRowAsHistogram (I, Graphics g, long irow, long colb, long cole,
	double ymin, double ymax, double xboffsetFraction, double xeoffsetFraction,
	double intercolumnFraction, double grey, double labeloffsetFraction, int garnish)
{
	char *proc = "TableOfReal_drawRowAsHistogram";
	iam (TableOfReal);
 
	long j, ncols, nx = my numberOfColumns;
	double bar_width = 1, xlab, xb, dx, x1, x2, y1, y2;

	if (cole >= colb) {colb = 1; cole = nx; }
	if (colb <= cole && colb < 1 || cole > nx)
	{
		    Melder_warning ("%s: Invalid columns", proc);
			return;
	}
	if (irow < 0 || irow > my numberOfRows)
	{ 
		Melder_warning ("%s: Invalid row", proc);
		return;
	}
	if (ymin >= ymax) NUMdvector_extrema (my data[irow], colb, cole, &ymin, &ymax);
	
    ncols = cole - colb + 1;

    bar_width /= ncols + xboffsetFraction + xeoffsetFraction + (ncols - 1) * intercolumnFraction;

    xb = xboffsetFraction * bar_width;
    dx = (intercolumnFraction + 1) * bar_width;
    
    Graphics_setWindow (g, 0, 1, ymin, ymax);
    Graphics_setInner (g);

    x1 = xb;
    for (j = colb; j <= cole; j++)
    {
		x2 = x1 + bar_width;
		y1 = ymin; y2 = my data[irow][j];
		if (y2 > ymin)
		{
			if (y2 > ymax) y2 = ymax;
			Graphics_setGrey (g, grey);
			Graphics_fillRectangle (g, x1, x2, y1, y2);
			Graphics_setGrey (g, 0); /* Black */
			Graphics_rectangle (g, x1, x2, y1, y2);
		}
        x1 += dx;
    }
    
    Graphics_unsetInner (g);

	if (garnish)
	{
		Graphics_drawInnerBox (g);
		xlab = labeloffsetFraction * bar_width;
	    for (j = colb; j <= cole; j++)
	    {
			if (my columnLabels[j]) Graphics_markBottom (g, xlab, 0, 0, 0, my columnLabels[j]);
			xlab += dx;
	    }
	}
}


void TableOfReal_drawBiplot (I, Graphics g, double xmin, double xmax, 
	double ymin, double ymax, double sv_splitfactor, int labelsize, int garnish)
{
	iam (TableOfReal);
	SVD svd;
	double lambda1, lambda2;
	float *x = NULL, *y = NULL;
	long i, numberOfZeroed;
	long nr = my numberOfRows, nc = my numberOfColumns, nmin;
	long nPoints = nr + nc; 
	int fontsize = Graphics_inqFontSize (g);

	svd = SVD_create (nr, nc);
	if (svd == NULL) goto end;
	
	NUMdmatrix_copyElements (my data, svd -> u, 1, nr, 1, nc);
	NUMcentreColumns_d (svd -> u, 1, nr, 1, nc, NULL);
	
	if (! SVD_compute (svd)) goto end; 	
	numberOfZeroed = SVD_zeroSmallSingularValues (svd, 0);
	
	nmin = MIN (nr, nc) - numberOfZeroed;
	if (nmin < 2)
	{
		Melder_warning ("TableOfReal_drawBiplot: There must be at least two "
			"(independent) columns in the table.");
		goto end;
	}

	x = NUMfvector (1, nPoints);
	if (x == NULL) goto end;
	y = NUMfvector (1, nPoints);
	if (y == NULL) goto end;

	lambda1 = pow (svd -> d[1], sv_splitfactor);
	lambda2 = pow (svd -> d[2], sv_splitfactor);
	for (i = 1; i <= nr; i++)
	{
		x[i] = svd -> u[i][1] * lambda1;
		y[i] = svd -> u[i][2] * lambda2;
	}
	lambda1 = svd -> d[1] / lambda1;
	lambda2 = svd -> d[2] / lambda2;
	for (i = 1; i <= nc; i++)
	{
			x[nr + i] = svd -> v[i][1] * lambda1;
			y[nr + i] = svd -> v[i][2] * lambda2;
	}
		
	if (xmax <= xmin) NUMfvector_extrema (x, 1, nPoints, &xmin, &xmax);
	if (xmax <= xmin) { xmax += 1; xmin -= 1; }
	if (ymax <= ymin) NUMfvector_extrema (y, 1, nPoints, &ymin, &ymax);
	if (ymax <= ymin) { ymax += 1; ymin -= 1; }
	
    Graphics_setWindow (g, xmin, xmax, ymin, ymax);	
    Graphics_setInner (g);
	if (labelsize > 0) Graphics_setFontSize (g, labelsize);			
    Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
	
	for (i = 1; i <= nPoints; i++)
	{
		char *label;
		if (i <= nr)
		{
			label = my rowLabels[i];
			if (label == NULL) label = "?__r_";
		}
		else
		{
			label = my columnLabels[i - nr];
			if (label == NULL) label = "?__c_";
		}
		Graphics_text (g, x[i], y[i], label);
	}
	
    Graphics_unsetInner (g);
	
	if (garnish)
	{
		Graphics_drawInnerBox (g);
		Graphics_marksLeft (g, 2, 1, 1, 0);
		Graphics_marksBottom (g, 2, 1, 1, 0);
	}
	
	if (labelsize > 0) Graphics_setFontSize (g, fontsize);
				
end:
	forget (svd);
	NUMfvector_free (x, 1);
	NUMfvector_free (y, 1);
}

/*
	Draw a box plot of data[1..ndata]. The vertical center line of the plot
	is at position 'x'. The rectangle box is 2*w wide, the whisker 2*r.
	All drawing outside [ymin, ymax] is clipped. 
*/
static void Graphics_drawBoxPlot (Graphics g, double data[], long ndata, 
	double x, double r, double w, double ymin, double ymax)
{
	double lowerOuterFence, lowerInnerFence, mean;
	double q75, q25, q50, upperInnerFence, upperOuterFence, hspread;
	double lowerWhisker, upperWhisker, y1, y2;
	int lineType = Graphics_inqLineType (g);
	long i, ie;

	Melder_assert (r > 0 && w > 0);
	if (ndata < 3) return;
	
	/*
		Sort the data (increasing: data[1] <= ... <= data[ndata]).
		Get the median (q50) and the upper and lower quartile points 
		(q25 and q75).
		Now q25 and q75 are the lower and upper hinges, respectively.
		The fances can be calcultaed from q25 and q75.
		The spread is defined as the interquartile range or midrange 
		|q75 - q25|.
		The fences are defined as:
		(lower/upper) innerfence = (lower/upper) hinge +/- 1.5 hspread
		(lower/upper) outerfence = (lower/upper) hinge +/- 3.0 hspread
	*/
		
	NUMsort_d (ndata, data);

	if (ymax <= ymin)
	{
		ymin = data[1]; ymax = data[ndata];
	}		
	if (data[1] > ymax || data[ndata] < ymin) return;

	for (mean=0, i=1; i <= ndata; i++) mean += data[i];
	mean /= ndata;
			
	q25 = NUMquantile_d (ndata, data, 0.25);
	q50 = NUMquantile_d (ndata, data, 0.5);
	q75 = NUMquantile_d (ndata, data, 0.75);
		
	hspread = fabs (q75 - q25);
	lowerOuterFence = q25 - 3.0 * hspread;
	lowerInnerFence = q25 - 1.5 * hspread;
	upperInnerFence = q75 + 1.5 * hspread;
	upperOuterFence = q75 + 3.0 * hspread;
		
	/*
		Decide whether there are outliers that have to be drawn.
		First process data from below (data are sorted).
	*/
	
	i = 1; ie = ndata;
	while (i <= ie && data[i] < ymin) i++;			
    Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
	while (i <= ie && data[i] < lowerOuterFence)
	{
		Graphics_text (g, x, data[i], "o"); i++;
	}
	while (i <= ie && data[i] < lowerInnerFence)
	{
		Graphics_text (g, x, data[i], "*"); i++;
	}
	lowerWhisker = data[i] < q25 ? data[i] : lowerInnerFence;
	if (lowerWhisker > ymax) return;
		
	/*
		Next process data from above.
	*/
		
	i = ndata; ie = i;				
	while (i >= ie && data[i] > ymax) i--;
	while (i >= ie && data[i] > upperOuterFence)
	{
		Graphics_text (g, x, data[i], "o"); i--;
	}
	while (i >= ie && data[i] > upperInnerFence)
	{
		Graphics_text (g, x, data[i], "*"); i--;
	}
	upperWhisker = data[i] > q75 ? data[i] : upperInnerFence;
	if (upperWhisker < ymin) return;
		
	/*
		Determine what parts of the "box" have to be drawn within the 
		range [ymin, ymax].
		Horizontal lines first.
	*/
		
	y1 = lowerWhisker;
	if (ymax > y1 && ymin < y1)
		Graphics_line (g, x - r, y1, x + r, y1);
	y1 = q25;
	if (ymax > y1 && ymin < y1)
		Graphics_line (g, x - w, y1, x + w, y1);
	y1 = q50;
	if (ymax > y1 && ymin < y1)
		Graphics_line (g, x - w, y1, x + w, y1);
	y1 = q75;
	if (ymax > y1 && ymin < y1)
		Graphics_line (g, x - w, y1, x + w, y1);
	y1 = upperWhisker;
	if (ymax > y1 && ymin < y1)
		Graphics_line (g, x - r, y1, x + r, y1);
		
	/*
		Extension: draw the mean too.
	*/
		
	y1 = mean;
	if (ymax > y1 && ymin < y1)
	{
		Graphics_setLineType (g, Graphics_DOTTED);
		Graphics_line (g, x - w, y1, x + w, y1);
		Graphics_setLineType (g, lineType); 
	}
		
	/*
		Now process the vertical lines.
	*/

	y1 = lowerWhisker; y2 = q25; 	
	if (ymax > y1 && ymin < y2)
	{
		y1 = MAX (y1, ymin);
		y2 = MIN (y2, ymax);
		Graphics_line (g, x, y1, x, y2);
	}
	y1 = q25; y2 = q75; 	
	if (ymax > y1 && ymin < y2)
	{
		y1 = MAX (y1, ymin);
		y2 = MIN (y2, ymax);
		Graphics_line (g, x - w, y1, x - w, y2);
		Graphics_line (g, x + w, y1, x + w, y2);
	}
	y1 = q75; y2 = upperWhisker; 	
	if (ymax > y1 && ymin < y2)
	{
		y1 = MAX (y1, ymin);
		y2 = MIN (y2, ymax);
		Graphics_line (g, x, y1, x, y2);
	}
}

void TableOfReal_drawBoxPlots (I, Graphics g, long rowmin, long rowmax, 
	long colmin, long colmax, double ymin, double ymax, int garnish)
{
	iam (TableOfReal); 
	double *data = NULL; 
	long i, j, numberOfRows;

	if (rowmax < rowmin || rowmax < 1)
	{
		rowmin = 1; rowmax = my numberOfRows;
	}
	if (rowmin < 1) rowmin = 1;
	if (rowmax > my numberOfRows) rowmax = my numberOfRows;
	numberOfRows = rowmax - rowmin + 1;
	if (colmax < colmin || colmax < 1)
	{
		colmin = 1; colmax = my numberOfColumns;
	}
	if (colmin < 1) colmin = 1; 
	if (colmax > my numberOfColumns) colmax = my numberOfColumns;
	if (ymax <= ymin) NUMdmatrix_extrema (my data, rowmin, rowmax, 
		colmin, colmax, &ymin, &ymax);
	if ((data = NUMdvector (1, numberOfRows)) == NULL) return;
	
    Graphics_setWindow (g, colmin - 0.5, colmax + 0.5, ymin, ymax);	
    Graphics_setInner (g);
	
	for (j=colmin; j <= colmax; j++)
	{
		double x = j, r = 0.05, w = 0.2, t; long ndata = 0;
		
		for (i=1; i <= numberOfRows; i++)
		{
			if ((t = my data[rowmin+i-1][j]) != NUMundefined) data[++ndata] = t;
		}
		Graphics_drawBoxPlot (g, data, ndata, x, r, w, ymin, ymax);
	}
	Graphics_unsetInner (g);
	if (garnish)
	{
		Graphics_drawInnerBox (g);
		for (j=colmin; j <= colmax; j++)
		{
			if (my columnLabels && my columnLabels[j] && my columnLabels[j][0])
				Graphics_markBottom (g, j, 0, 1, 0, my columnLabels [j]);
		}
    	Graphics_marksLeft (g, 2, 1, 1, 0);
	}
}

int TableOfReal_equalLabels (I, thou, int rowLabels, int columnLabels)
{
	iam (TableOfReal); thouart (TableOfReal); long i;
	Melder_assert (rowLabels || columnLabels);
	if (rowLabels)
	{
		if (my numberOfRows != thy numberOfRows) return 0;
		if (my rowLabels == thy rowLabels) return 1;
		for (i=1; i <= my numberOfRows; i++)
		{
			if (NUMstrcmp (my rowLabels[i], thy rowLabels[i])) return 0;
		}
	}
	if (columnLabels)
	{
		if (my numberOfColumns != thy numberOfColumns) return 0;
		if (my columnLabels == thy columnLabels) return 1;
		for (i=1; i <= my numberOfColumns; i++) 
		{
			if (NUMstrcmp (my columnLabels[i], thy columnLabels[i])) return 0;
		}
	}
	return 1;
}

int TableOfReal_copyLabels (I, thou, int rowOrigin, int columnOrigin)
{
	iam (TableOfReal); thouart (TableOfReal); long i;
	if (rowOrigin == 1)
	{
		if (my numberOfRows != thy numberOfRows) return 0;
		for (i=1; i <= thy numberOfRows; i++) if (my rowLabels[i])
			TableOfReal_setRowLabel (thee, i, my rowLabels[i]);
	}
	else if (rowOrigin == -1)
	{
		if (my numberOfColumns != thy numberOfRows) return 0;
		for (i=1; i <= thy numberOfRows; i++) if (my columnLabels[i]) 
			TableOfReal_setRowLabel (thee, i, my columnLabels[i]);
	}
	if (columnOrigin == 1)
	{
		if (my numberOfColumns != thy numberOfColumns) return 0;
		for (i=1; i <= thy numberOfColumns; i++) if (my columnLabels[i])
			TableOfReal_setColumnLabel (thee, i, my columnLabels[i]);
	}
	else if (columnOrigin == -1)
	{
		if (my numberOfRows != thy numberOfColumns) return 0;
		for (i=1; i <= thy numberOfColumns; i++) if (my rowLabels[i])
			TableOfReal_setColumnLabel (thee, i, my rowLabels[i]);
	}
	return 1;
}

void TableOfReal_labelsFromCollectionItemNames (I, thou, int row, int column)
{
	iam (TableOfReal); 
	thouart (Collection); 
	long i; 
	char *name;
	
	if (row)
	{
		Melder_assert (my numberOfRows == thy size);
		for (i = 1; i <= my numberOfRows; i++)
		{
			name = Thing_getName (thy item[i]);
			if (name != NULL) TableOfReal_setRowLabel (me, i, name);
		}
	}
	if (column)
	{
		Melder_assert (my numberOfColumns == thy size);
		for (i=1; i <= my numberOfColumns; i++)
		{
			name = Thing_getName (thy item[i]);
			if (name != NULL) TableOfReal_setColumnLabel (me, i, name);
		}
	}
}

void TableOfReal_centreColumns (I)
{
	iam (TableOfReal);
	NUMcentreColumns_d (my data, 1, my numberOfRows, 1, my numberOfColumns,
		NULL);
}

int TableOfReal_and_Categories_setRowLabels (I, Categories thee)
{
	iam (TableOfReal);
	Categories c = NULL;
	long i;
	if (my numberOfRows != thy size) return Melder_error 
		("TableOfReal_and_Categories_swap: unequal number of items.");
		
	/* 
		If anything goes wrong we must leave the Table intact.
		We first copy the Categories, swap the labels 
		and then delete the newly created categories.
	*/
	
	c = Data_copy (thee);
	if (c == NULL) return 0;
	
	for (i=1; i <= my numberOfRows; i++)
	{
		SimpleString s = c -> item[i];
		char *t = s -> string;
		s -> string = my rowLabels[i];
		my rowLabels[i] = t;
	}
	
	forget (c);
	return 1;
}

void TableOfReal_centreColumns_byRowLabel (I)
{
	iam (TableOfReal); char *label = my rowLabels[1];
	long i, index = 1;
		
	for (i=2; i <= my numberOfRows; i++)
	{
		char *li = my rowLabels[i];
		if (li != NULL && li != label && strcmp (li, label))
		{
			NUMcentreColumns_d (my data, index, i - 1, 1, my numberOfColumns, NULL);
			label = li; index = i;
		}
	}
	NUMcentreColumns_d (my data, index, my numberOfRows, 1, my numberOfColumns, NULL);
}

double TableOfReal_getRowSum (I, long index)
{
	iam (TableOfReal);
	double sum = 0;
	long j;
	if (index < 1 || index > my numberOfRows)
	{
		return NUMundefined;
	}
	for (j = 1; j <= my numberOfColumns; j++)
	{
		sum += my data[index][j];
	}
	return sum;
}

double TableOfReal_getColumnSum (I, long index)
{
	iam (TableOfReal);
	double sum = 0;
	long i;
	if (index < 1 || index > my numberOfColumns)
	{
		return NUMundefined;
	}
	for (i = 1; i <= my numberOfRows; i++)
	{
		sum += my data[i][index];
	}
	return sum;
}

double TableOfReal_getGrandSum (I)
{
	iam (TableOfReal);
	double sum = 0;
	long i, j;
	for (i = 1; i <= my numberOfRows; i++)
	{
		for (j = 1; j <= my numberOfColumns; j++)
		{
			sum += my data[i][j];
		}
	}
	return sum;
}

void TableOfReal_centreRows (I)
{
	iam (TableOfReal);
	NUMcentreRows_d (my data, 1, my numberOfRows, 1, my numberOfColumns);
}

void TableOfReal_doubleCentre (I)
{
	iam (TableOfReal);
	NUMdoubleCentre_d (my data, 1, my numberOfRows, 1, my numberOfColumns);
}

void TableOfReal_normalizeColumns (I, double norm)
{
	iam (TableOfReal);
	NUMnormalizeColumns_d (my data, my numberOfRows, my numberOfColumns, norm);
}

void TableOfReal_normalizeRows (I, double norm)
{
	iam (TableOfReal);
	NUMnormalizeRows_d (my data, my numberOfRows, my numberOfColumns, norm);
}

void TableOfReal_standardizeColumns (I)
{
	iam (TableOfReal);
	NUMstandardizeColumns (my data, 1, my numberOfRows, 1, my numberOfColumns);
}

void TableOfReal_standardizeRows (I)
{
	iam (TableOfReal);
	NUMstandardizeRows (my data, 1, my numberOfRows, 1, my numberOfColumns);
}

void TableOfReal_normalizeTable (I, double norm)
{
	iam (TableOfReal);
	NUMnormalize_d (my data, my numberOfRows, my numberOfColumns, norm);
}

double TableOfReal_getTableNorm (I)
{
	iam (TableOfReal); 
	double sumsq = 0; 
	long i, j;
	for (i = 1; i <= my numberOfRows; i++)
	{
		for (j = 1; j <= my numberOfColumns; j++)
		{
			sumsq += my data[i][j] * my data[i][j];
		}
	}
	return sqrt (sumsq);
}

int TableOfReal_checkPositive (I)
{
	iam (TableOfReal);
	long i, j, negative = 0;
	
	for (i = 1; i <= my numberOfRows; i++)
	{
		for (j = 1; j <= my numberOfColumns; j++)
		{
			if (my data[i][j] < 0) { negative ++; break; }
		}
	}
	return negative == 0 ? 1 : 
		Melder_error ("All matrix entries should be positive!");
}

/* NUMundefined ??? */
void NUMdmatrix_getColumnExtrema (double **a, long rowb, long rowe, long icol, double *min, double *max);
void NUMdmatrix_getColumnExtrema (double **a, long rowb, long rowe, long icol, double *min, double *max)
{
	long i;
	*min = *max = a[rowb][icol];
	for (i=rowb+1; i <= rowe; i++)
	{
		double t = a[i][icol];
		if (t > *max) *max = t;
		else if (t < *min) *min = t;
	} 
}

void TableOfReal_drawScatterPlotMatrix (I, Graphics g, long colb, long cole, 
	double fractionWhite)
{
	iam (TableOfReal); 
	double *xmin = NULL, *xmax = NULL;
	long i, j, k, m = my numberOfRows, n;
	
	if (colb == 0 && cole == 0)
	{
		colb = 1; cole = my numberOfColumns;
	}
	else if (cole < colb || colb < 1 || cole > my numberOfColumns) return;
	
	n = cole - colb + 1;
	if (n == 1) return;
		
	if (! (xmin = NUMdvector (colb, cole)) ||
		! (xmax = NUMdvector (colb, cole))) goto end;
	
	for (j=colb; j <= cole; j++)
	{
		xmin[j] = xmax[j] = my data[1][j];
	}
	for (i=2; i <= m; i++)
	{
		for (j=colb; j <= cole; j++)
		{
			if (my data[i][j] > xmax[j]) xmax[j] = my data[i][j];
			else if (my data[i][j] < xmin[j]) xmin[j] = my data[i][j];
		}
	}
	for (j=colb; j <= cole; j++)
	{
		double extra = fractionWhite * fabs (xmax[j] - xmin[j]);
		if (extra == 0) extra = 0.5;
		xmin[j] -= extra; xmax[j] += extra;
	}
	
	Graphics_setWindow (g, 0, n, 0, n);
	Graphics_line (g, 0, n, n, n);
	Graphics_line (g, 0, 0, 0, n);
	Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
	
	for (i=1; i <= n; i++)
	{
		long xcol, ycol = colb + i - 1; char *mark, label[20];
		Graphics_line (g, 0, n - i, n, n - i);
		Graphics_line (g, i, n, i, 0);
		for (j=1; j <= n; j++)
		{
			xcol = colb + j -1;
			if (i == j)
			{
				mark = my columnLabels[xcol];
				if (! mark)
				{
					sprintf (label, "Column %d", xcol); mark = label;
				}
				Graphics_text (g, j - 0.5, n - i + 0.5, mark);
			}
			else
			{
				for (k=1; k <= m; k++)
				{
					double x = j - 1 + (my data[k][xcol] - xmin[xcol]) / 
						(xmax[xcol] - xmin[xcol]);
					double y = n - i + (my data[k][ycol] - xmin[ycol]) / 
						(xmax[ycol] - xmin[ycol]); 
					mark = EMPTY_STRING (my rowLabels[k]) ? "+" : my rowLabels[k];
					Graphics_text (g, x, y, mark);
				}
			}
		}
	}
end:
	NUMdvector_free (xmin, colb); NUMdvector_free (xmax, colb);
}

void TableOfReal_drawScatterPlot (I, Graphics g, long icx, long icy, long rowb, 
	long rowe, double xmin, double xmax, double ymin, double ymax, int noRowLabels, 
	int garnish)
{
    iam (TableOfReal); 
	double tmp, m = my numberOfRows, n = my numberOfColumns;
    long i, ix = labs (icx), iy = labs (icy);
    
    if (ix < 1 || ix > n || iy < 1 || iy > n) return;
    if (rowb < 1) rowb = 1;
    if (rowe > m) rowe = m;
    if (rowe <= rowb)
    {
    	rowb = 1; rowe = m;
    }
    
    if (xmax <= xmin)
    {
		NUMdmatrix_getColumnExtrema (my data, rowb, rowe, ix, & xmin, & xmax);
		tmp = xmax - xmin == 0 ? 0.5 : 0.0;
		xmin -= tmp; xmax += tmp;
    }
    if (ymax <= ymin)
    {
		NUMdmatrix_getColumnExtrema (my data, rowb, rowe, iy, & ymin, & ymax);
		tmp = ymax - ymin == 0 ? 0.5 : 0.0;
		ymin -= tmp; ymax += tmp;
    }
    if (icx < 0)
    {
    	double t = xmin; xmin = xmax; xmax = t;
    }
    if (icy < 0)
    {
    	double t = ymin; ymin = ymax; ymax = t;
    }
	
    Graphics_setInner (g);
    Graphics_setWindow (g, xmin, xmax, ymin, ymax);
	Graphics_setTextAlignment (g, Graphics_CENTRE, Graphics_HALF);
	
    for (i=rowb; i <= rowe; i++)
    {
    	double x = my data[i][ix], y = my data[i][iy];
    	char *mark = noRowLabels || EMPTY_STRING (my rowLabels[i]) ? "+" : my rowLabels[i];
    	if (x >= xmin && x <= xmax && y >= ymin && y <= ymax)
    		Graphics_text (g, x, y, mark);
    }
	
    Graphics_unsetInner (g);
	
    if (garnish)
    {
		Graphics_drawInnerBox (g);
		if (my columnLabels[ix]) Graphics_textBottom (g, 1, my columnLabels[ix]);
		if (my columnLabels[iy]) Graphics_textLeft (g, 1, my columnLabels[iy]);
		Graphics_marksLeft (g, 2, 1, 1, 0);
		Graphics_marksBottom (g, 2, 1, 1, 0);
		if (ymax * ymin < 0 && xmax * xmin < 0)
		{
			Graphics_markLeft (g, 0, 0, 1, 1, "");
			Graphics_markBottom (g, 0, 0, 1, 1, "");
		}
	}
}

/****************  TABLESOFREAL **************************************/

class_methods (TablesOfReal, Ordered)
class_methods_end

int TablesOfReal_init (I, void *klas)
{
	iam (TablesOfReal);
	if (! me || ! Ordered_init (me, klas, 10)) return 0;
	return 1;
}

TablesOfReal TablesOfReal_create (void)
{
	TablesOfReal me = new (TablesOfReal);
	if (! me || ! TablesOfReal_init (me, classTableOfReal)) forget (me);
	return me;
}

TableOfReal TablesOfReal_sum (I)
{
	iam (TablesOfReal); TableOfReal thee;
	long i, j, k;
	if (my size <= 0) return NULL;
	if (! (thee = Data_copy (my item[1]))) { forget (thee); return NULL; }
	for (i=2; i <= my size; i++)
	{
		TableOfReal him = my item[i];
		if (thy numberOfRows != his numberOfRows || thy numberOfColumns != his numberOfColumns ||
			! TableOfReal_equalLabels (thee, him, 1, 1))
		{
			forget (thee);
			return Melder_errorp ("TablesOfReal_sum: dimensions or labels of items 1 and %d differ", i);
		}
		for (j=1; j <= thy numberOfRows; j++)
			for (k=1; k <= thy numberOfColumns; k++) thy data[j][k] += his data[j][k]; 
	}
	return thee;
}

int TablesOfReal_checkDimensions (I)
{
	iam (TablesOfReal); TableOfReal t1; long i;
	if (my size < 2) return 1;
	t1 = my item[1];
	for (i=2; i <= my size; i++)
	{
		TableOfReal t = my item[i];
		if (t -> numberOfColumns != t1 -> numberOfColumns ||
			t -> numberOfRows != t1 -> numberOfRows) return 0;
	}
	return 1;
}

double TableOfReal_getColumnQuantile (I, long col, double quantile)
{
	iam (TableOfReal); long i, m = my numberOfRows;
	double *values, r;
	
	if (col < 1 || col > my numberOfColumns) return NUMundefined;
	
	if (! (values = NUMdvector (1, m))) return NUMundefined;
	
	for (i=1; i <= m; i++)
	{
		values[i] = my data[i][col];
	}
	
	NUMsort_d (m, values);
	r = NUMquantile_d (m, values, quantile);
	
	NUMdvector_free (values, 1);
	
	return r;
}

/*
	Recognize a string of integers separated by whitespace or ranges like 1:3
	or 3:1 or a combination of both. 
	Only accepted characters are '0123456789', colon ':' and space ' '.
	All input numbers must be >= min and <= max.
	selection = "1 5 14 6:8"
	Returns: numberOfIntegers = 6, integers[]=1,5,14,6,7,8
	selection = "6:2"
	Returns: numberOfIntegers = 5, integers[]=6,5,4,3,2
*/
static int parse (const char *selection, long min, long max, 
	long *integers, long *numberOfIntegers)
{
	int digit = 0, colon = 0;
	const char *p, *pc;
	
	*numberOfIntegers = 0;

	p = pc = selection;
	while (1)
	{
		if (isdigit (*p))
		{
			if (! digit)
			{
				pc = p;
				digit = 1;
			}
		} 
		else if (*p == ' ' || *p == '\0' || (*p == ':' && ! colon))
		{
			if (digit)
			{
				long i, number = atol (pc), nplus = 1;
				
				if (number < min || number > max) return Melder_error 
					("parse: Index \"%d\" out of allowed range [%d,%d]",
						number, min, max);
				
				if (colon)
				{
					long number1 = integers[*numberOfIntegers];
					int ascending = number1 <= number;
					
					if (! ascending)
					{
						long itmp = number1; number1 = number; number = itmp;
					}
					nplus = number - number1;
					for (i = 1; i <= *numberOfIntegers - 1; i++)
					{
						if (integers[i] >= number1 && integers[i] <= number)
							return Melder_error ("parse: Range [%d:%d] "
								"overlaps previous selected column %d.",
								number1, number, integers[i]);
					}
					
					if (ascending)
					{
						for (i = 1; i <= nplus; i++)
						{
							integers[*numberOfIntegers + i] = number1 + i;
						}		
					}
					else
					{
						for (i = 1; i <= nplus; i++)
						{
							integers[*numberOfIntegers + i] = number - i;
						}		
					}
					*numberOfIntegers += nplus;
					colon = 0;
				}
				else
				{
					for (i = 1; i <= *numberOfIntegers; i++)
					{
						if (integers[i] == number) return Melder_error ("parse: "
							"You cannot select a column (%d) more than once.", number);
					}
					*numberOfIntegers += 1;
					integers[*numberOfIntegers] = number;
				}
				
				digit = 0;
			}
			if (*p == ':') colon = 1;
			else if (*p == '\0') break;
		}
		else
		{
			return Melder_error ("parse: illegal character at position %d in "
				"string \"%s\".", p - selection + 1, selection);
		}
		p++;
	}
	return 1;
}

TableOfReal TableOfReal_selectColumnsWhereRow (I, const char *columnSelection, 
	const char *rowConditionFormula)
{
	char *proc = "TableOfReal_selectColumnsWhereRow";
	iam (TableOfReal);
	TableOfReal thee = NULL;
	long i, j, k, irow, ncols, nrows, *colSelect = NULL;
	long nx = my numberOfColumns, ny = my numberOfRows;
	double result;
	
	colSelect = NUMlvector (1, nx);
	if (colSelect == NULL) return NULL;
	if (! parse (columnSelection, 1, nx, colSelect, &ncols)) goto end;

	if (ncols == 0)
	{
		(void) Melder_error ("%s: No columns selected.", proc);
		goto end;
	}
		
	/*
		Determine the number of rows in the new table:
		select a row if a number <> 0 occurs in it.
	*/

	if (! Formula_compile (NULL, me, rowConditionFormula, FALSE, TRUE)) return NULL;
	
	for (nrows = 0, i = 1; i <= ny; i++)
	{
		for (j = 1; j <= nx; j++)
		{
			if (! Formula_run (i, j, &result, NULL)) goto end;
			if (result != 0.0)
			{
				nrows++; break;
			}
		}
	}
	
	if (nrows < 1)
	{
		(void) Melder_error ("%s: No rows selected", proc);
		goto end;
	}
		
	thee = TableOfReal_create (nrows, ncols);
	if (thee == NULL) goto end;
	
	for (k = 1; k <= ncols; k++)
	{
		char *label = my columnLabels[colSelect[k]];
		if (label)
		{
			TableOfReal_setColumnLabel (thee, k, label);
		} 
	}
	
	for (irow = 0, i = 1; i <= ny; i++)
	{
		for (j = 1; j <= nx; j++)
		{
			if (! Formula_run (i, j, &result, NULL)) goto end;
			if (result != 0.0)
			{
				irow++;
				for (k = 1; k <= ncols; k++)
				{
					thy data[irow][k] = my data[i][colSelect[k]];
				}
				if (my rowLabels[i])
				{
					TableOfReal_setRowLabel (thee, irow, my rowLabels[i]);
				}
				break;
			}
		}
	}
	
end:

	NUMlvector_free (colSelect, 1);
	if (Melder_hasError()) forget (thee);
	return thee;
}

/*TableOfReal TableOfReal_selectColumnsWhereRowoud (I, const char *columnSelection, 
	const char *rowConditionFormula)
{
	iam (TableOfReal);
	TableOfReal thee = NULL;
	long i, j, k, irow, ncols, nrows;
	long nx = my numberOfColumns, ny = my numberOfRows;
	Matrix original = NULL, conditions = NULL;
	long *colSelect = NULL;
	char *name = "TableOfReal_selectColumnsWhereRow";

	colSelect = NUMlvector (1, nx);
	if (colSelect == NULL) return NULL;
	if (! parse (columnSelection, 1, nx, colSelect, &ncols)) goto end;

	if (ncols == 0)
	{
		(void) Melder_error ("%s: No columns selected.");
		goto end;
	}
	
	original = TableOfReal_to_Matrix (me);
	if (original == NULL) goto end;
	conditions = Data_copy (original);
	if (conditions == NULL) goto end;
	
	if (! Matrix_formula (original, rowConditionFormula, 0, 0, 0, 
			conditions)) goto end;
	
	for (nrows = 0, i = 1; i <= ny; i++)
	{
		for (j = 1; j <= nx; j++)
		{
			if (conditions -> z[i][j] != 0.0)
			{
				nrows++; break;
			}
		}
	}
	
	if (nrows < 1)
	{
		(void) Melder_error ("%s: No rows selected", name);
		goto end;
	}
		
	thee = TableOfReal_create (nrows, ncols);
	if (thee == NULL) goto end;
	
	for (k = 1; k <= ncols; k++)
	{
		char *label = my columnLabels[colSelect[k]];
		if (label)
		{
			TableOfReal_setColumnLabel (thee, k, label);
		} 
	}
	
	for (irow = 0, i = 1; i <= ny; i++)
	{
		for (j = 1; j <= nx; j++)
		{
			if (conditions -> z[i][j] != 0.0)
			{
				irow++;
				for (k = 1; k <= ncols; k++)
				{
					thy data[irow][k] = my data[i][colSelect[k]];
				}
				if (my rowLabels[i])
				{
					TableOfReal_setRowLabel (thee, irow, my rowLabels[i]);
				}
				break;
			}
		}
	}
	
end:

	forget (original);
	forget (conditions);
	NUMlvector_free (colSelect, 1);
	if (Melder_hasError()) forget (thee);
	return thee;
}*/

TableOfReal TableOfReal_createFromPolsData_50males (int include_levels)
{
	TableOfReal me = new (TableOfReal);
	long i, j, nrows = 600, ncols = include_levels ? 6 : 3;
	char *columnLabels[6] = {"F__1_", "F__2_", "F__3_", 
		"L__1_", "L__2_", "L__3_"};
	char *vowelLabels[12] = {"u", "a", "o", "\\as", "\\o/", "i",
		"y", "e", "\\yc", "\\ep", "\\ct", "\\ic"};
		
	short data[600][6] = {
    /* speaker 1 */
    320,  630,  2560,  6,  13,  48,
    780, 1300,  2460,  6,   8,  30,
    500,  940,  2420,  3,  12,  35,
    720, 1060,  2420,  3,   8,  27,
    430, 1580,  2260,  2,  24,  36,
    280, 2300,  2780, 14,  22,  27,
    320, 1680,  2140,  6,  23,  30,
    420, 2000,  2620,  5,  20,  23,
    420, 1540,  2380,  4,  19,  24,
    600, 1720,  2700,  3,  17,  29,
    520, 1000,  2520,  4,  13,  31,
    350, 2000,  2520,  7,  19,  18,
    /* speaker 2 */
    440,  780,  2600,  7,  20,  35, 
    940, 1300,  2780,  5,  13,  26,  
    500,  740,  2700,  7,  11,  32,
    800, 1000,  2480, 11,  11,  31,
    460, 1500,  2300, 10,  20,  20,
    320, 2400,  3040,  9,  28,  27,
    340, 1600,  2050,  7,  30,  32,
    420, 2200,  2650,  4,  10,  12,
    540, 1370,  2320,  6,  23,  31,
    760, 1660,  2600,  8,  13,  18,
    560,  800,  2800,  5,  17,  37,
    430, 2030,  2660,  7,  18,  17,
    /* speaker 3 */
    280,  740,  2160,  8,  32,  46,
    860, 1500,  2580,  7,  13,  19,
    480,  820,  2280,  4,  11,  33,
    680, 1020,  2460,  9,  14,  25,
    360, 1520,  2080,  5,  16,  21,
    270, 2040,  2860,  5,  18,  18,
    280, 1600,  1900,  5,  16,  21,
    380, 1940,  2580,  5,  15,  17,
    400, 1520,  2120,  4,  21,  24,
    560, 1840,  2520,  7,  18,  22,
    500,  820,  2520,  3,  11,  30,
    350, 2000,  2660,  4,  15,  19,
    /* speaker 4 */
    360,  820,  2220,  6,  11,  30,
    840, 1300,  2280,  9,   7,  22,
    500,  900,  2320,  7,  13,  35,
    680, 1000,  2480, 12,  12,  27,
    440, 1320,  2060,  5,  12,  27,
    240, 2060,  2580,  5,  20,  23,
    300, 1540,  2020,  4,  18,  24,
    460, 1920,  2460,  7,  18,  17,
    480, 1320,  2200,  6,  15,  22,
    660, 1660,  2340,  7,  16,  21,
    500,  800,  2520,  7,  14,  34,
    440, 1920,  2380,  3,  19,  18,
    /* speaker 5 */
    440,  880,  2300,  6,  22,  41,
    820, 1420,  2180,  6,   7,  28,
    540,  960,  2460,  4,  10,  28,
    780, 1040,  2620,  8,  10,  32,
    540, 1540,  2160,  2,  17,  19,
    300, 2300,  2900,  5,  23,  30,
    360, 1860,  2200,  2,  15,  15,
    520, 1960,  2400,  4,  21,  21,
    560, 1440,  2280,  2,  17,  27,
    700, 1720,  2300,  7,  13,  22,
    600,  880,  3000,  5,   5,  25,
    560, 1920,  2400,  5,  22,  25,
    /* speaker 6 */
    260,  700,  2550,  3,  24,  45,
    820, 1460,  2760, 10,  15,  27,
    450,  900,  2460, 19,  20,  45,
    700, 1080,  2660, 13,  22,  32,
    460, 1750,  2300,  7,  32,  40,
    240, 2500,  3000,  3,  33,  37,
    260, 2100,  2500, 10,  37,  43,
    300, 2320,  2860,  8,  26,  26,
    440, 1700,  2660,  7,  27,  32,
    560, 2080,  2840, 12,  15,  16,
    550,  900,  2740,  7,  25,  30,
    340, 2340,  3000,  8,  31,  31,
    /* speaker 7 */
    280,  860,  2340,  5,  15,  33,
    800, 1320,  2540,  7,  14,  26,
    520,  920,  2600,  8,  15,  33,
    600, 1000,  2760,  7,   7,  25,
    450, 1660,  2260,  7,  20,  24,
    260, 2340,  2640,  3,  17,  20,
    280, 1780,  2160,  3,  25,  29,
    400, 2040,  2400,  9,  19,  21,
    460, 1560,  2400,  4,  19,  22,
    620, 1760,  2560,  6,  18,  25,
    560,  960,  2760,  9,  16,  33,
    340, 2000,  2600,  6,  14,  18,
    /* speaker 8 */
    320,  880,  2200,  6,  16,  40,
    800, 1160,  2600,  6,  13,  26,
    560,  980,  2360,  5,   8,  35,
    700, 1080,  2540, 10,  14,  34,
    500, 1480,  2300,  8,  17,  29,
    280, 2080,  2620,  3,  27,  27,
    320, 1760,  2060,  4,  25,  32,
    400, 1940,  2540,  9,  17,  22,
    400, 1560,  2280,  6,  18,  28,
    540, 1860,  2540,  4,  14,  19,
    560,  920,  2320,  8,  13,  32,
    340, 1960,  2480,  6,  13,  17,
    /* speaker 9 */
    300,  680,  2400,  3,  19,  32,
    860, 1300,  2660, 12,  18,  26,  
    500,  940,  2500,  3,  13,  33,
    700, 1120,  2620,  9,  17,  24,
    500, 1500,  2280,  3,  17,  22,
    300, 2380,  2960,  2,  20,  22,
    300, 1760,  2160,  1,  19,  23,
    480, 2100,  2580,  4,  15,  15,
    500, 1580,  2400,  5,  12,  22,
    640, 1700,  2620,  3,  15,  19,
    560,  900,  2940,  4,   8,  31,
    400, 2040,  2600,  7,  17,  19,
    /* speaker 10 */
    360,  900,  2220, 11,  21,  46,
    880, 1400,  2660,  9,  17,  29,
    460,  940,  2400,  3,  13,  37,
    660, 1040,  2660,  5,   5,  31,
    460, 1580,  2360,  4,  22,  26,
    340, 2200,  2920, 14,  30,  26,
    400, 1880,  2800, 12,  19,  33,
    460, 2080,  2800,  4,  18,  17,
    460, 1480,  2260,  5,  27,  31,
    600, 1860,  2640,  7,  18,  23,
    520,  860,  2880,  2,  24,  32,
    460, 2100,  2800,  6,  19,  22,
    /* speaker 11 */
    320,  830,  2060,  5,  16,  45,
    820, 1340,  2200,  7,   9,  26,
    520,  840,  2040,  7,  13,  31,
    660, 1060,  2300,  8,   9,  28,
    440, 1520,  2040,  7,  16,  22,
    300, 2100,  2600,  5,  23,  30,
    300, 1740,  2040,  2,  17,  19,
    340, 2040,  2460,  5,  20,  24,
    500, 1440,  2200,  8,  17,  21,
    600, 1760,  2380, 10,  15,  19,
    560,  900,  2300,  9,  10,  28,
    300, 1960,  2400,  7,  22,  24,
    /* speaker 12 */
    400,  860,  2700,  7,  15,  46,
    940, 1520,  3040, 10,  19,  33,
    580, 1040,  2960,  5,  17,  42,
    860, 1280,  3000,  8,  19,  34,
    480, 1600,  2620,  7,  15,  32,
    300, 2500,  2880,  3,  29,  29,
    300, 1670,  2350,  3,  31,  39,
    480, 2220,  2640,  7,  30,  31,
    380, 1600,  2540,  4,  22,  33,
    630, 2140,  2880,  7,  27,  27,
    640,  900,  3000,  8,  13,  39,
    360, 2220,  2780,  5,  23,  31,
    /* speaker 13 */
    260,  780,  2460,  5,  20,  44,
    900, 1560,  2860,  9,  13,  26,
    440,  850,  2600,  3,  17,  42,
    860, 1140,  2820, 10,  12,  30,
    460, 1580,  2400,  3,  19,  27,
    260, 2560,  3240,  3,  28,  34,
    300, 1960,  2500,  2,  31,  29,
    460, 2320,  2960,  6,  25,  27,
    460, 1600,  2460,  6,  22,  30,
    680, 2100,  2940, 14,  19,  25,
    540,  800,  2740,  7,  13,  36,
    380, 2500,  2980,  7,  20,  23,
    /* speaker 14 */
    340,  720,  2500,  6,  16,  47,
    900, 1500,  3020, 10,  13,  33,
    450,  900,  2700,  4,  16,  42,
    600, 1000,  2720, 12,  10,  31,
    420, 1740,  2560,  8,  21,  25,
    360, 2500,  3000,  8,  10,  13,
    360, 1900,  2420,  6,  16,  12,
    380, 1780,  2420,  4,  36,  18,
    500, 1640,  2620,  4,  23,  29, 
    600, 1940,  2700,  3,  16,  19,
    500,  800,  2800,  3,  10,  38,
    400, 2360,  2740,  7,  13,  27,
    /* speaker 15 */
    360,  780,  2320,  5,  23,  50,
    860, 1420,  2420,  8,  15,  35,
    440,  840,  2480,  9,  20,  38,
    660,  980,  2500,  6,  13,  33,
    460, 1660,  2200,  2,  21,  31,
    300, 2360,  2960,  7,  29,  27,
    320, 1740,  2220,  2,  28,  30,
    400, 2240,  2560,  5,  23,  24,
    480, 1420,  2220,  6,  23,  37,
    680, 1640,  2340,  9,  21,  28,
    560,  860,  2780,  4,  15,  44,
    440, 2120,  2500,  2,  22,  23,
    /* speaker 16 */
    360,  760,  2300,  7,  23,  52,
    660, 1000,  2500,  9,  15,  35,
    500,  920,  2520,  7,  15,  37,
    780, 1060,  2380,  5,   9,  38,
    440, 1560,  2260,  7,  25,  31,
    280, 2200,  2880,  7,  43,  38,
    380, 1720,  2200,  7,  29,  39,
    360, 2140,  2620,  3,  26,  28,
    360, 1600,  2400,  4,  26,  33,
    520, 1800,  2480,  8,  32,  35,
    540,  920,  2640,  6,  20,  44,
    340, 2080,  2680,  3,  28,  27,
    /* speaker 17 */
    400,  820,  2200,  5,  15,  48,
    1100,1480,  2260, 10,  14,  30,
    520,  940,  2560,  5,  15,  43,
    660,  940,  2820, 11,  11,  35,
    500, 1720,  2400,  5,  19,  23,
    360, 2300,  3260, 11,  25,  17,
    360, 2100,  2420, 10,  19,  20,
    440, 2360,  2860,  6,  20,  25,
    500, 1760,  2600,  6,  15,  25,
    660, 1840,  2620,  7,  17,  24,
    540,  860,  2860,  3,   9,  41,
    400, 2440,  3000,  5,  28,  30,
    /* speaker 18 */
    360,  860,  2520,  6,  15,  35,
    740, 1300,  2660,  7,   9,  22,
    460,  800,  2620,  3,  20,  35,
    740, 1040,  2800, 15,  13,  30,
    440, 1400,  2200,  4,  21,  25,
    340, 2040,  2500,  6,  18,  20,
    340, 1340,  2040,  3,  15,  25,
    420, 1760,  2420,  6,  18,  21,
    460, 1380,  2200,  6,  20,  22,
    560, 1640,  2400,  9,  13,  19,
    540,  920,  2520,  9,  10,  23,
    400, 1800,  2400,  8,  19,  21,
    /* speaker 19 */
    320,  840,  2360,  5,  20,  45,
    780, 1140,  2740,  8,  15,  36,
    460, 1020,  2700,  6,  11,  43,
    800, 1100,  2720, 15,  17,  42,
    440, 1500,  2500,  5,  23,  37,
    260, 2000,  2680,  4,  39,  40,
    300, 1540,  2100,  5,  33,  42,
    400, 1900,  2680,  6,  25,  34,
    440, 1500,  2400,  4,  26,  42,
    600, 1700,  2640, 18,  33,  40,
    500,  800,  3000,  5,  20,  47,
    400, 2000,  2500,  5,  27,  37,
    /* speaker 20 */
    400,  960,  2400, 10,  16,  37,
    800, 1220,  2380,  7,  13,  25,
    500, 1080,  2500, 10,  14,  31,
    780, 1100,  2600, 10,  11,  30,
    400, 1480,  2380,  8,  20,  25,
    280, 2380,  2720,  6,  23,  24,
    300, 1760,  2220,  5,  14,  22,
    400, 2000,  2600,  4,  20,  24,
    440, 1500,  2440,  4,  15,  20,
    440, 1800,  2620, 12,  18,  22,
    460,  860,  2600, 12,  13,  33,
    400, 2040,  2640,  3,  18,  23,
    /* speaker 21 */
    300,  700,  2100,  3,  15,  37,
    800, 1100,  2300, 11,  15,  30,
    420,  700,  2440,  5,  12,  33,
    660,  920,  2520,  6,   9,  32,
    400, 1300,  2000,  2,  22,  30,
    300, 1940,  2620,  5,  25,  24,
    260, 1920,  2900,  4,  21,  31,
    300, 1900,  2340,  4,  20,  24,
    400, 1200,  2000,  4,  20,  29,
    540, 1500,  2280,  4,  15,  22,
    400,  740,  2580,  3,  11,  34,
    300, 1900,  2380,  2,  20,  20,
    /* speaker 22 */
    400,  780,  2500,  6,  14,  40,
    760, 1260,  2620,  6,  14,  30,
    540,  860,  2600,  4,   7,  34,
    660, 1100,  2460,  9,  12,  30,
    540, 1460,  2260,  5,  12,  20,
    300, 2300,  2800,  7,  20,  17,
    300, 1980,  2900,  6,  14,  32,
    420, 2100,  2600,  4,  17,  21,
    500, 1440,  2300,  4,  18,  26,
    540, 1900,  2600, 10,  20,  22,
    550,  900,  3000,  5,  11,  34,
    300, 2200,  2700,  6,  14,  20,
    /* speaker 23 */
    360,  900,  2140,  3,  14,  37,
    800, 1250,  2650, 12,  13,  26,
    520,  960,  2200,  5,   7,  28,
    760, 1120,  2700, 12,  10,  29,
    400, 1500,  2160,  7,  20,  30,
    300, 2260,  3000,  4,  18,  17,
    320, 1800,  2500,  6,  17,  32,
    460, 2020,  2800,  3,  10,  16,
    500, 1500,  2300,  7,  15,  23,
    640, 1600,  2500,  8,  22,  32,
    550,  940,  2420,  7,  12,  40,
    460, 2100,  2880,  2,  21,  24,
    /* speaker 24 */
    360,  860,  2460,  4,  15,  33,
    840, 1400,  2500,  5,   8,  20,
    460,  900,  2520,  1,   8,  31,
    620, 1020,  2770,  7,   6,  26,
    410, 1460,  2360,  2,  11,  18,
    270, 2140,  2580,  2,  15,  21,
    300, 1870,  2300,  1,  11,  16,
    360, 2000,  2520,  2,  17,  19,
    400, 1520,  2400,  2,  12,  17,
    600, 1600,  2580,  3,  12,  20,
    500,  900,  2700,  2,  11,  33,
    360, 1940,  2550,  1,  19,  21,
    /* speaker 25 */
    360,  860,  2200,  8,  22,  46,
    880, 1240,  2400,  7,  14,  25,
    460,  920,  3300,  2,  11,  40,
    600, 1000,  2600,  4,   7,  32,
    400, 1440,  2160,  2,  31,  37,
    360, 2240,  2760,  6,  23,  21,
    380, 1660,  2000, 13,  16,  26,
    460, 2000,  2520,  5,  23,  32,
    460, 1500,  2300,  7,  17,  25,
    540, 1700,  2460, 10,  22,  35,
    540, 1000,  2600,  8,  13,  40,
    340, 2040,  2580,  9,  27,  30,
    /* speaker 26 */
    400,  800,  2500,  6,  23,  42,
    960, 1300,  2640,  9,   8,  29,
    460,  860,  2460,  7,  16,  39,
    740, 1140,  2400,  9,  10,  30,
    400, 1600,  2400,  6,  21,  25,
    360, 2500,  2840, 10,  19,  20,
    360, 1800,  2400,  8,  21,  27,
    360, 2080,  2680,  5,  14,  16,
    400, 1620,  2440,  5,  15,  21,
    600, 1940,  2600,  5,  13,  23,
    560,  980,  2900,  5,  13,  32,
    400, 2060,  2540,  5,  21,  21,
    /* speaker 27 */
    300,  900,  2300,  3,  15,  38, 
    780, 1300,  2400,  9,  18,  32,
    550, 1000,  2480,  5,  10,  35,
    680, 1050,  2550,  5,  10,  35,
    520, 1480,  2400,  5,  16,  27,
    260, 2180,  2560,  1,  30,  30,
    250, 1720,  2220,  1,  26,  32,
    360, 2100,  2650,  4,  31,  25,
    440, 1440,  2440,  4,  21,  26,
    600, 1600,  2500,  5,  15,  27,
    560,  950,  2700,  3,  10,  40,
    360, 1900,  2600,  3,  26,  31,
    /* speaker 28 */
    280,  740,  2500,  3,  20,  45,
    780, 1300,  2840,  4,  12,  25,
    440,  860,  2860,  5,  13,  41,
    440,  700,  3040,  9,  10,  37,
    450, 1520,  2320,  2,  25,  30,
    220, 2340,  2960,  1,  35,  36,
    240, 1800,  2140,  3,  36,  37,
    300, 2200,  2600,  4,  28,  33,
    440, 1500,  2480,  3,  26,  37,
    500, 1660,  2620, 10,  30,  38,
    420,  700,  3000,  4,  12,  35,
    300, 2140,  2760,  6,  32,  31,
    /* speaker 29 */
    340,  660,  2320,  4,  13,  37,
    640, 1250,  2480, 10,  16,  29, 
    560, 1000,  2480, 10,  14,  31,
    720, 1150,  2600, 10,  13,  26,
    480, 1400,  2160,  9,  22,  28,
    300, 2040,  2640,  5,  20,  22,
    280, 1540,  1960,  3,  21,  23,
    460, 1760,  2320,  7,  19,  20,
    440, 1550,  2200,  8,  16,  23,
    480, 1660,  1960, 10,  16,  23,
    480,  840,  2840,  9,  12,  28,
    400, 1780,  2360,  7,  20,  23, 
    /* speaker 30 */
    360,  800,  2540,  1,  11,  40,
    600, 1300,  2600,  6,   8,  27, 
    500,  860,  2440,  2,   7,  36,
    750, 1140,  2640,  4,   9,  30,
    460, 1400,  2340,  1,  23,  28,
    340, 2300,  2620,  2,  26,  25,
    300, 1540,  2300,  5,  25,  35,
    440, 2000,  2540,  1,  14,  19,
    440, 1360,  2360,  1,  19,  26,
    620, 1840,  2560,  3,  18,  23,
    520,  820,  2680,  2,   7,  34,
    420, 2000,  2640,  1,  20,  25, /*L2 (10) corrected 20021211 */
    /* speaker 31 */
    340,  740,  2240,  5,  15,  44,
    820, 1200,  2250,  5,  17,  40,
    440,  820,  2540,  3,  13,  37,
    760, 1060,  2340,  8,  15,  42,
    460, 1540,  2380,  3,  21,  25,
    280, 2260,  2620,  8,  31,  32,
    300, 1800,  2220,  6,  36,  37,
    460, 1900,  2260,  3,  34,  31,
    380, 1540,  2400, 15,  36,  40,
    500, 1740,  2400,  4,  26,  39,
    460,  840,  2580,  7,  28,  35,
    320, 2100,  2460,  7,  30,  27,
    /* speaker 32 */
    360,  900,  2200,  5,  18,  39,
    640, 1280,  2340,  7,  15,  26,
    460,  920,  2360,  7,  14,  33,
    720, 1200,  2580, 11,  14,  26,
    420, 1520,  2260,  6,  17,  23,
    400, 2000,  2560,  5,  15,  20,
    380, 1700,  2100,  6,  18,  21,
    440, 1740,  2420,  4,  14,  17,
    500, 1520,  2440,  5,  14,  21,
    580, 1540,  2460,  5,  12,  24,
    580, 1020,  2700,  9,  13,  30,
    460, 1720,  2400,  7,  22,  20,
    /* speaker 33 */
    400,  700,  2600,  4,  17,  45,
    900, 1440,  2600, 10,  17,  36, 
    460,  860,  2600,  7,  18,  45,
    680, 1000,  2200,  7,  13,  39,
    460, 1600,  2540,  5,  28,  37,
    300, 2260,  2880,  7,  28,  26,
    320, 1860,  2200,  7,  22,  28,
    440, 2180,  2660,  3,  24,  28,
    380, 1560,  2360,  8,  26,  33,
    620, 1720,  2060,  6,  21,  28,
    600,  860,  2900, 12,  18,  37,
    440, 2040,  2600,  4,  26,  30,
    /* speaker 34 */
    370,  900,  2230,  3,  17,  35,
    700, 1200,  2580, 12,  17,  30,
    500,  840,  2460,  4,  13,  37,
    720, 1080,  2640,  7,  13,  37,
    440, 1300,  2220,  4,  20,  31,
    300, 2040,  2580,  7,  25,  25,
    320, 1540,  2080,  8,  21,  23,
    380, 1860,  2450,  7,  24,  33,
    460, 1200,  2360,  3,  20,  30,
    580, 1500,  2380,  9,  22,  25,
    480,  820,  2580,  8,  15,  32,
    400, 1800,  2360,  6,  23,  26,
    /* speaker 35 */
    280, 1040,  2340,  7,  26,  41,
    820, 1300,  2760, 10,  15,  32,
    440, 1220,  2580,  6,  18,  29,
    600, 1040,  2540,  8,  13,  27,
    420, 1560,  2480,  6,  22,  26,
    300, 2160,  2700,  5,  29,  30,
    250, 1760,  2320,  9,  30,  38,
    440, 1940,  2550,  5,  25,  28,
    400, 1600,  2460,  8,  26,  29,
    580, 1820,  2460,  5,  23,  30,
    460,  860,  2660,  5,  21,  37,
    400, 2100,  2640,  8,  27,  27,
    /* speaker 36 */
    360,  740,  2160,  2,  21,  40,
    660, 1260,  2540, 10,  18,  21, 
    500,  900,  2600,  9,  20,  30,
    640, 1000,  2880, 11,  17,  29,
    460, 1300,  2140,  8,  19,  25,
    300, 1900,  2580, 11,  18,  22,
    320, 1660,  2060,  9,  17,  20,
    400, 1780,  2320,  8,  20,  21,
    380, 1360,  2200,  6,  20,  25,
    540, 1600,  2260,  8,  21,  22,
    540,  860,  2720,  7,  20,  32,
    400, 1740,  2340,  5,  22,  23,
    /* speaker 37 */
    300,  900,  2140,  5,  21,  39,  
    700, 1240,  2460, 10,  18,  28,
    480,  960,  2140,  4,  12,  32,
    640, 1120,  2480,  9,  14,  32,
    460, 1520,  2160,  4,  18,  25,
    320, 2120,  2600, 10,  20,  26,
    320, 1800,  2200,  8,  25,  27,
    320, 1920,  2460,  8,  21,  27,
    480, 1460,  2260,  5,  22,  27,
    600, 1600,  2480,  7,  17,  24,
    500,  950,  2450,  4,  14,  38,
    460, 1820,  2480,  6,  18,  26,
    /* speaker 38 */
    320,  760,  2080, 11,  23,  41,
    840, 1180,  2700, 13,  20,  37, 
    500,  920,  2400, 13,  17,  37,
    660, 1060,  2700, 13,  17,  36,
    440, 1400,  2220,  5,  32,  37,
    280, 2240,  2700, 14,  29,  35,
    300, 1640,  2080, 12,  31,  31,
    440, 2040,  2600, 11,  19,  24,
    400, 1460,  2160,  8,  32,  38,
    580, 1700,  1900, 13,  26,  26,
    500,  840,  2920, 10,  18,  40,
    360, 2060,  2440,  7,  21,  27,
    /* speaker 39 */
    320,  760,  2480,  9,  21,  46,
    700, 1420,  2680,  9,  16,  31,
    500,  940,  2500,  4,  16,  41,
    700, 1060,  2720,  8,  10,  38,
    440, 1580,  2260, 11,  34,  39,
    260, 2200,  2700,  5,  29,  32,
    200, 1600,  2060,  7,  33,  34,
    400, 2200,  2600, 13,  29,  31,
    380, 1500,  2220,  8,  25,  37,
    540, 1750,  2420,  8,  23,  32,
    520,  820,  2560, 10,  19,  43,
    400, 1700,  2320,  3,  38,  23,
    /* speaker 40 */
    300,  680,  1920,  7,  28,  48,
    740, 1200,  2550,  8,  10,  24, 
    420,  860,  2420,  7,  17,  37,
    640, 1120,  2500, 12,  17,  37,
    360, 1500,  2180,  3,  27,  35,
    280, 2160,  2920,  4,  27,  31,
    260, 1560,  2050,  2,  26,  27,
    360, 2020,  2500,  4,  26,  28,
    440, 1400,  2320,  4,  21,  32,
    460, 1660,  2460,  5,  19,  27,
    500,  840,  2580,  6,  14,  35,
    360, 1920,  2560,  3,  31,  31,
    /* speaker 41 */
    360,  880,  2320,  2,  12,  43,  
    840, 1200,  2500, 12,  17,  37,
    580, 1060,  2300, 11,  10,  32,
    580, 1100,  2680,  8,  12,  33,
    560, 1600,  2200,  9,  22,  27,
    300, 2260,  2800,  5,  25,  26,
    320, 1760,  2100,  5,  23,  25,
    500, 2020,  2660,  5,  17,  21,
    420, 1520,  2320,  1,  20,  28,
    700, 1800,  2620,  8,  17,  22,
    540,  860,  2720,  7,  13,  35,
    420, 2080,  2600,  4,  21,  25,
    /* speaker 42 */
    420,  800,  2400,  5,  15,  40,
    800, 1400,  2900,  9,  16,  41,
    420,  820,  2480,  5,   8,  40,
    600, 1200,  2760,  6,  12,  34,
    400, 1560,  2120,  2,  20,  31,
    320, 2360,  2820,  8,  25,  27,
    340, 1680,  2240,  9,  19,  35,
    400, 2180,  2760,  2,  19,  20,
    400, 1440,  2360,  3,  15,  26,
    700, 1700,  2340, 11,  18,  29,
    500,  780,  2840,  7,  14,  38,
    380, 2120,  2720,  2,  21,  25,
    /* speaker 43 */
    300,  760,  2020,  3,  16,  38, 
    740, 1200,  2360,  8,  15,  29,
    460,  860,  2200,  8,  12,  39,
    620,  900,  2500,  8,  12,  27,
    400, 1340,  2100,  7,  20,  31,
    240, 2000,  2340,  2,  22,  28,
    240, 1580,  1860,  3,  16,  24,
    360, 1640,  2080,  5,  19,  26,
    400, 1340,  2060,  3,  16,  28,
    580, 1400,  2120,  6,  13,  24,
    500,  800,  2460,  6,   6,  31,
    440, 1720,  2100,  7,  19,  24,
    /* speaker 44 */
    260,  800,  2400,  3,  16,  48, 
    780, 1300,  2700,  6,  14,  28,
    480,  900,  2500,  5,   8,  35,
    620, 1000,  2820,  5,   9,  28,
    420, 1400,  2300,  3,  18,  29,
    240, 2040,  2680,  1,  31,  28,
    260, 1580,  2260,  3,  31,  28,
    380, 2000,  2600,  5,  29,  26,
    420, 1420,  2400,  2,  21,  26,
    540, 1640,  2440,  5,  19,  26,
    480,  840,  2800,  6,  13,  32,
    280, 1960,  2560,  5,  27,  28,
    /* speaker 45 */
    300,  840,  3060,  3,  10,  38,
    800, 1220,  2280,  6,  10,  26,
    500,  920,  2120,  6,   8,  31,
    700, 1020,  2600,  3,  11,  28,
    400, 1260,  2020,  6,  17,  24,
    260, 1960,  2440,  1,  22,  22,
    300, 1480,  1940,  2,  18,  22,
    440, 1880,  2380,  6,  17,  17,
    320, 1400,  2140,  5,  18,  27,
    500, 1560,  2300,  7,  18,  22,
    540,  780,  2400,  8,  13,  34,
    360, 1860,  2300,  4,  20,  21,
    /* speaker 46 */
    320,  860,  2380,  3,  19,  41,
    660, 1400,  2540, 11,  20,  27,
    520,  940,  2580,  7,  11,  34,
    700, 1040,  2720,  4,   8,  23,
    400, 1600,  2280,  2,  27,  29,
    320, 2340,  3140,  3,  33,  29,
    300, 1860,  2160,  2,  25,  25,
    420, 2200,  2760,  1,  17,  23,
    460, 2320,  3360,  2,  33,  37,
    500, 2100,  2760,  6,  23,  38,
    600,  920,  2700,  8,  17,  29,
    420, 2200,  2740,  3,  30,  32,
    /* speaker 47 */
    360,  800,  2120,  3,  18,  33,
    700, 1220,  2760,  6,  13,  27,
    540,  940,  2640,  2,   8,  33,
    620, 1080,  2800,  4,  10,  33,
    500, 1400,  2200,  6,  18,  25,
    320, 2240,  2940,  2,  22,  27,
    320, 1800,  2100,  1,  26,  27,
    420, 2040,  2400,  3,  19,  24,
    460, 1440,  2140,  1,  13,  25,
    600, 1600,  2520,  6,  13,  27,
    560,  700,  2780,  5,  11,  31,
    440, 1920,  2560,  3,  22,  24,
    /* speaker 48 */
    300,  760,  1900,  3,  17,  42,
    800, 1260,  2740,  7,  10,  28,
    460,  840,  1840,  4,  13,  38,
    540,  900,  2400, 10,  14,  28,
    420, 1380,  2100,  2,  16,  29,
    220, 2080,  2900,  2,  28,  21,
    220, 1760,  2120,  1,  22,  25,
    440, 2060,  2780,  1,  19,  21,
    440, 1440,  2560,  3,  19,  31,
    580, 1400,  2100,  5,  15,  22,
    520,  900,  2300,  3,  10,  32,
    420, 1720,  2720,  6,  22,  21,
    /* speaker 49 */
    320, 1000,  2220,  3,  24,  43,
    700, 1280,  2500,  3,  13,  30,
    460, 1060,  2380,  4,  13,  31,
    620, 1100,  2840, 10,  18,  33,
    340, 1440,  2260,  4,  21,  30,
    280, 2140,  2580,  3,  31,  32,
    280, 1820,  2220,  2,  36,  35,
    340, 2100,  2500,  2,  29,  31,
    380, 1460,  2400,  2,  26,  38,
    500, 1640,  2500,  7,  27,  31,
    500,  960,  2720,  4,  17,  31,
    420, 1960,  2700,  1,  32,  32,
    /* speaker 50 */
    340,  780,  2020, 11,  22,  36,
    660, 1220,  2500,  9,  14,  22,
    420,  760,  2440,  2,  17,  33,
    560, 1000,  2600,  6,  13,  25,
    400, 1320,  2120,  8,  18,  21,
    300, 1860,  2440,  6,  22,  22,
    280, 1600,  1900,  6,  16,  19,
    340, 1740,  2260,  3,  12,  17,
    400, 1360,  2160,  4,  16,  22,
    520, 1580,  2240,  2,  12,  16,
    380,  800,  2560,  7,  11,  25,
    360, 1740,  2260,  5,  14,  17};
		
	if (me == NULL || ! TableOfReal_init (me, nrows, ncols))
	{
		forget (me); return NULL;
	}
	for (j = 1; j <= ncols; j++)
	{
		TableOfReal_setColumnLabel (me, j, columnLabels[j-1]);
	}
	for (i = 1; i <= nrows; i++)
	{
		for (j = 1; j <= ncols; j++)
		{
			my data [i][j] = data[i-1][j-1];
		}
		TableOfReal_setRowLabel (me, i, vowelLabels[(i-1) % 12]);
	}	 	
	return me;
}

TableOfReal TableOfReal_randomizeRows (TableOfReal me)
{
	TableOfReal thee = NULL;
	Sequence s = NULL;
	long i;

	thee = TableOfReal_create (my numberOfRows, my numberOfColumns);
	if (thee == NULL) return NULL;

	/*
		Copy column labels.
	*/

	if (! NUMstrings_copyElements (my columnLabels, thy columnLabels, 1, 
		my numberOfColumns)) goto end;

	s = Sequence_create (my numberOfRows, SEQUENCE_PERMUTE, 0);
	if (s == NULL) goto end;
	
	for (i = 1; i <= my numberOfRows; i++)
	{
		long p = Sequence_ith (s, i);
		NUMdvector_copyElements (my data[p], thy data[i], 
			1, my numberOfColumns);
		if (my rowLabels[p])
		{
			TableOfReal_setRowLabel (thee, i, my rowLabels[p]);
		}
	}
	
end:
	forget (s);
	if (Melder_hasError ()) forget (thee);
	return thee;
}

TableOfReal TableOfReal_bootstrap (TableOfReal me)
{
	TableOfReal thee = NULL;
	long i;

	thee = TableOfReal_create (my numberOfRows, my numberOfColumns);
	if (thee == NULL) return NULL;

	/*
		Copy column labels.
	*/

	for (i = 1; i <= my numberOfColumns; i++)
	{
		if (my columnLabels[i])
		{
			TableOfReal_setColumnLabel (thee, i, my columnLabels[i]);
		}
	}

	/*
		Select randomly from table with replacement. Because of replacement
		you do not get back the original data set. A random fraction, 
		typically 1/e (37%) are replaced by duplicates.
	*/
	
	for (i = 1; i <= my numberOfRows; i++)
	{
		long p = NUMrandomInteger (1, my numberOfRows);
		NUMdvector_copyElements (my data[p], thy data[i], 
			1, my numberOfColumns);
		if (my rowLabels[p])
		{
			TableOfReal_setRowLabel (thee, i, my rowLabels[p]);
		}
	}	

	if (Melder_hasError ()) forget (thee);
	return thee;
}

int TableOfReal_changeRowLabels (I, char *search, char *replace, 
	int maximumNumberOfReplaces, long *nmatches, long *nstringmatches, 
	int use_regexp)
{
	iam (TableOfReal);
	char ** rowLabels = strs_replace (my rowLabels, 1, my numberOfRows, 
		search, replace, maximumNumberOfReplaces, nmatches, 
		nstringmatches, use_regexp);
	if (rowLabels == NULL) return 0;
	NUMstrings_free (my rowLabels, 1, my numberOfRows);
	my rowLabels = rowLabels;
	return 1;
}

int TableOfReal_changeColumnLabels (I, char *search, char *replace, 
	int maximumNumberOfReplaces, long *nmatches, long *nstringmatches, 
	int use_regexp)
{
	iam (TableOfReal);
	char ** columnLabels = strs_replace (my columnLabels, 1, my numberOfColumns, 
		search, replace, maximumNumberOfReplaces, nmatches, 
		nstringmatches, use_regexp);
	if (columnLabels == NULL) return 0;
	NUMstrings_free (my columnLabels, 1, my numberOfColumns);
	my columnLabels = columnLabels;
	return 1;
}

long TableOfReal_getNumberOfLabelMatches (I, char *search, int columnLabels, 
	int use_regexp)
{
	iam (TableOfReal);
	long i, nmatches = 0, numberOfLabels = my numberOfRows;
	char **labels = my rowLabels;
	regexp *compiled_regexp = NULL;
	
	if (search == NULL || strlen (search) == 0) return 0;
	if (columnLabels)
	{
		numberOfLabels = my numberOfColumns;
		labels = my columnLabels;
	}	
	if (use_regexp)
	{
		char *compileMsg;
		compiled_regexp = CompileRE (search, &compileMsg, 0);
		if (compiled_regexp == NULL) return Melder_error (compileMsg);
	}
	for (i = 1; i <= numberOfLabels; i++)
	{
		if (labels[i] == NULL) continue;
		if (use_regexp)
		{
			if (ExecRE (compiled_regexp, NULL, labels[i], NULL, 0, 
				'\0', '\0', NULL)) nmatches++;
		}
		else if (strequ (labels[i], search)) nmatches++;
	}
	if (use_regexp) free (compiled_regexp);
	return nmatches;
}

void TableOfReal_drawVectors (I, Graphics g, long col1, long col2, double xmin,
	double xmax, double ymin, double ymax, int vectype, int labelsize,
	int garnish)
{
	iam (TableOfReal);
	long i, nx = my numberOfColumns, ny = my numberOfRows;
	double min, max;
	int fontsize = Graphics_inqFontSize (g);

	if (col1 < 1 || col1 >= nx || col2 < 1 || col2 >= nx) return;
	if (col1 == col2) return;
 
	if (xmin >= xmax)
	{
		NUMdmatrix_extrema (my data, 1, ny, col1, col1, &min, &max);
		NUMdmatrix_extrema (my data, 1, ny, col2, col2, &xmin, &xmax);
		if (min < xmin) xmin = min;
		if (max > xmax) xmax = max;
	}
	if (ymin >= ymax)
	{
		NUMdmatrix_extrema (my data, 1, ny, col1+1, col1+1, &min, &max);
		NUMdmatrix_extrema (my data, 1, ny, col2+1, col2+1, &ymin, &ymax);
		if (min < ymin) ymin = min;
		if (max > ymax) ymax = max;
	}
	if (xmin == xmax)
	{
		if (ymin == ymax) return;
		xmin -= 0.5;
		xmax += 0.5;
	}
	if (ymin == ymax)
	{
		ymin -= 0.5;
		ymax += 0.5;
	}

	Graphics_setWindow (g, xmin, xmax, ymin, ymax);
	Graphics_setInner (g);
	if (labelsize > 0) Graphics_setFontSize (g, labelsize);			
	for (i = 1; i <= ny; i++)
	{
		float x1 = my data[i][col1];
		float y1 = my data[i][col1 + 1];	 
		float x2 = my data[i][col2];
		float y2 = my data[i][col2 + 1];	 
		if (vectype == Graphics_LINE)
			Graphics_line (g, x1, y1, x2, y2);
		else if (vectype == Graphics_TWOWAYARROW)
		{
			Graphics_arrow (g, x1, y1, x2, y2);
			Graphics_arrow (g, x2, y2, x1, y1);
		}
		else /*if (vectype == Graphics_ARROW) */	
			Graphics_arrow (g, x1, y1, x2, y2);
		if (labelsize <= 0) continue;
		
	}
	if (labelsize > 0) Graphics_setFontSize (g, fontsize);
	Graphics_unsetInner (g);
	if (garnish)
	{
	    Graphics_drawInnerBox (g);
    	Graphics_marksLeft (g, 2, 1, 1, 0);
    	Graphics_marksBottom (g, 2, 1, 1, 0);
	}
}

TableOfReal TableOfReal_sortRowsByIndex (I, long *index, int reverse)
{
	iam (TableOfReal);
	TableOfReal thee = NULL;
	double min, max;
	long i, j;

	if (my rowLabels == NULL) return NULL;
	
	NUMlvector_extrema (index, 1, my numberOfRows, &min, &max);
	if (min < 1 || max > my numberOfRows) return Melder_errorp
		("TableOfReal_sortRowsByIndex: one or more indices out of range [1, %d].",
		my numberOfRows);
	
	thee = TableOfReal_create (my numberOfRows, my numberOfColumns);
	if (thee == NULL) return NULL;
	
	for (i = 1; i <= my numberOfRows; i++)
	{
		long    myindex = reverse ? i : index[i];
		long   thyindex = reverse ? index[i] : i;
		char   *mylabel = my rowLabels[myindex];
		double  *mydata = my data[myindex];
		double *thydata = thy data[thyindex];
		
		/*
			Copy the row label
		*/
		
		if (mylabel != NULL)
		{
			thy rowLabels[i] = Melder_strdup (mylabel);
			if (thy rowLabels[i] == NULL) goto end;
		}
		
		/*
			Copy the row values
		*/
		
		for (j = 1; j <= my numberOfColumns; j++)
		{
			thydata[j] = mydata[j];
		}
	}
	
	/*
		Copy column labels.
	*/
	
	(void) NUMstrings_copyElements (my columnLabels, thy columnLabels, 
		1, my numberOfColumns);
end:

	if (Melder_hasError()) forget (thee);
	return thee;
}

long *TableOfReal_rowLabelsIndex (I)
{
	iam (TableOfReal); 
	long *index = NUMlvector (1, my numberOfRows);
	
	if (index == NULL) return NULL;
	NUMindexx_s (my rowLabels, my numberOfRows, index);
	return index;
}

TableOfReal TableOfReal_sortByRowLabels_simple (I)
{
	iam (TableOfReal); 
	TableOfReal thee = NULL;
	long *index = NULL;
	
	index = TableOfReal_rowLabelsIndex (me);
	if (index == NULL) return NULL;
	
	thee = TableOfReal_sortRowsByIndex (me, index, 0);

	NUMlvector_free (index, 1);
	return thee;
}

TableOfReal TableOfReal_averageColumnsByRowLabels (I)
{
	iam (TableOfReal);
	TableOfReal thee = NULL, sorted = NULL;
	char *label, **tmp;
	long *index = NULL, indexi = 1, i;

	index = TableOfReal_rowLabelsIndex (me);
	if (index == NULL) return NULL;
	
	sorted = TableOfReal_sortRowsByIndex (me, index, 0);
	if (sorted == NULL) goto end;

	label = sorted -> rowLabels[1];
	for (i = 2; i <= my numberOfRows; i++)
	{
		char *li = sorted -> rowLabels[i];
		if (li != NULL && li != label && strcmp (li, label))
		{
			NUMaverageColumns (sorted -> data, indexi, i - 1, 1, 
				my numberOfColumns);
			label = li; indexi = i;
		}
	}
		
	NUMaverageColumns (sorted -> data, indexi, my numberOfRows, 1, 
		my numberOfColumns);

	/*
		Now inverse the table.
	*/
	tmp = sorted -> rowLabels; sorted -> rowLabels = my rowLabels;
	thee = TableOfReal_sortRowsByIndex (sorted, index, 1);
	sorted -> rowLabels = tmp;
	
end:
	forget (sorted);
	NUMlvector_free (index, 1);
	if (Melder_hasError()) forget (thee);
	return thee;
}

TableOfReal TableOfReal_rankColumns (I)
{
	iam (TableOfReal);
	TableOfReal thee = Data_copy (me);
	
	if (thee == NULL) return NULL;
	if (! NUMrankColumns_d (thy data, 1, thy numberOfRows, 
		1, thy numberOfColumns)) forget (thee);
	return thee;	
}
	
int TableOfReal_setSequentialColumnLabels (I, long from, long to, 
	char *precursor, long number, long increment)
{
	iam (TableOfReal);
	
	if (from == 0) from = 1;
	if (to == 0) to = my numberOfColumns;
	if (from < 1 || from > my numberOfColumns || to < from ||
		to > my numberOfColumns) return Melder_error 
			("TableOfReal_setSequentialColumnLabels: wrong column indices.");
	return NUMstrings_setSequentialNumbering (my columnLabels, from, to, 
		precursor, number, increment);
}
	
int TableOfReal_setSequentialRowLabels (I, long from, long to, 
	char *precursor, long number, long increment)
{
	iam (TableOfReal);
	
	if (from == 0) from = 1;
	if (to == 0) to = my numberOfRows;
	if (from < 1 || from > my numberOfRows || to < from ||
		to > my numberOfRows) return Melder_error 
			("TableOfReal_setSequentialRowLabels: wrong row indices.");
	return NUMstrings_setSequentialNumbering (my rowLabels, from, to, 
		precursor, number, increment);
}

/* For the inheritors */
TableOfReal TableOfReal_to_TableOfReal (I)
{
	iam (TableOfReal);
	TableOfReal thee = TableOfReal_create (my numberOfRows, my numberOfColumns);
	if (thee == NULL) return NULL;
	
	NUMdmatrix_copyElements (my data, thy data, 1, my numberOfRows, 
		1, my numberOfColumns);
	if (! TableOfReal_copyLabels (me, thee, 1, 1)) forget (thee);
	return thee;
}

TableOfReal TableOfReal_choleskyDecomposition (TableOfReal me, int upper, int inverse)
{
	char *proc = "TableOfReal_choleskyDecomposition", uplo = 'U', diag = 'N';
	long i, j, n = my numberOfColumns, lda = my numberOfRows, info;
	TableOfReal thee;

	if (n != lda) return Melder_errorp("%s: Must be a square symmetric matrix.", proc);
	if ((thee = Data_copy (me)) == NULL) return NULL;

	if (upper)
	{
		uplo = 'L'; /* Fortran storage */
		for (i = 2; i <= n; i++) for (j = 1; j < i; j++) thy data[i][j] = 0;
	}
	else
	{
		for (i = 1; i < n; i++) for (j = i+1; j <= n; j++) thy data[i][j] = 0;
	}
	(void) NUMlapack_dpotf2 (&uplo, &n, &thy data[1][1], &lda, &info);
	if (info != 0) goto end;

	if (inverse)
	{
		(void) NUMlapack_dtrtri (&uplo, &diag, &n, &thy data[1][1], &lda, &info);   
	}

end:
	if (Melder_hasError()) forget (thee);
	return thee;
}

#undef EMPTY_STRING
#undef MAX
#undef MIN

/* End of file TableOfReal_extensions.c */
