lib.c

#include "lib.h"
#include <cstring>
/*
	new --> calloc and delete --> free:
	May 20, 2004: 
		int** get_cols

	free --> delete:
	May 20, 2004:
		zfree

	ttd:
	1.	in fgets, check for overrun.

--------------------------------------------------------------------------------

	cat2		- cat together two strings.
	choose(k,n)	- n choose k.
	drand		- random number generator.
	entropy		- naive estimate of entropy.
	entropy2	- naive estimate of second moment of log(p).
	get_chars	- parse file and return list of character
			  strings pointing to each element on each line.
	get_cols	- get columns from file (including stdin).
	header_info	- get info from headers in ".dat" files.
	isnumber	- turns string into a number with lots of checking.
	kl_divergence	- Kullback-Leibler distance between two distributions.
	length_string	- returns length of \0 terminated string, including \0.
	max		- finds maximum of array.
	min		- finds minimum of array.
	mmult		- multiply matrix times matrix or matrix times vector.
	cint		- 2-D integer array, using calloc.
	cfloat		- 2-D float array, using calloc.
	cdouble		- 2-D double array, using calloc.
	newchar		- 2-D character array.
	newdouble	- 2-D double array, using new.
	newfloat	- 2-D float array.
	newint		- 2-D integer array.
	init_constant	- initizlize 1 and 2 D arrays.
	mi		- computes mutual information
	rabs		- absolute value.
	read_binary	- read binary file
	read_n		- read n lines of a file and return data in kth column.
	read_plexon	- read plexon file.
	readline	- read line of file.
	parse_c		- parse data file -- return pointer to string.
	parse_comline	- extract info from command line.
	s_parse		- simple string parse.
	unit_no		- return cell number in unit format given
			  input in _either_ plexon or unit format.
	write_err	- write an error message and exit(1).
*/

char* cat2(char* prefix, char* suffix)
{
/*
	cat together two strings and return the result
*/
	char* dum = new char[80];
	strcpy(dum, prefix);
	strcat(dum, suffix);
	return dum;
}

int choose(int k, int n)
{
/*
	compute n!/(k! (n-k)!)

	will get overflow if n and k are too big.
*/
	if (k > n) return 0;

 	int c = 1;
	k = k > n/2 ? n-k: k;
	for (int i=0; i < k; i++) c *= (n-i);
	for (int i=0; i < k; i++) c /= (i+1);

	return c;
}

double drand()
{
	// return a random number between 0 and 1
	// use srand(seed) to start seed.
	// 4/4/00:  RAND_MAX --> RAND_MAX-1, so 1 cannot be returned.
	// 8/4/04:  RAND_MAX --> RAND_MAX. why I thought RAND_MAX-1
	//                       would fix this is a mystery. would
	//                       like to go to RAND_MAX+1, but this
	//                       causes overflow and returns negative
	//                       numbers.
	return (double) rand()/RAND_MAX;
}

double drand(float x)
{
	// return a random number between 0 and x
	return x*rand()/RAND_MAX;
}

double drand(float xmin, float xmax)
{
	// return a random number between xmin and xmax
	return xmin + (xmax-xmin)*rand()/RAND_MAX;
}

double entropy(int n, int* p)
{
/*
	overloaded function (int* p --> double s)

	compute entropy of a probability distribution, p.
*/
	// ---find ntot
	float ntot=0;
	for (int i=0; i < n; i++) ntot += p[i] > 0 ? p[i] : 0;

	double s=0;
	for (int i=0; i < n; i++)
		s -= p[i] > 0 ? (p[i]/ntot)*log(p[i]/ntot) : 0;

	return s/log(2);
}

double entropy(int n, float* p, int norm)
{
/*
	overloaded function (float* p --> double s)

	compute entropy of a probability distribution, p.
*/
	double s=0;
	for (int i=0; i < n; i++) s -= p[i] > 0 ? p[i]*log(p[i]) : 0;

	if (norm)
	{
		// normalize probability distribution
		double z=0;
		for (int i=0; i < n; i++) z += p[i] > 0 ? p[i] : 0;
		if (z > 0) s = s/z + log(z);
	}

	return s/log(2);
}

double entropy(int n, double* p, int norm)
{
/*
	overloaded function (double* p --> double s)

	compute entropy of a probability distribution, p.
*/
	double s=0;
	for (int i=0; i < n; i++) s -= p[i] > 0 ? p[i]*log(p[i]) : 0;

	if (norm)
	{
		// normalize probability distribution
		double z=0;
		for (int i=0; i < n; i++) z += p[i] > 0 ? p[i] : 0;
		if (z > 0) s = s/z + log(z);
	}

	return s/log(2);
}

double entropy2(int n, int* p)
{
/*
	overloaded function (int* p --> double s)

	compute second moment of log p.
*/
	// ---find ntot
	float ntot=0;
	for (int i=0; i < n; i++) ntot += p[i] > 0 ? p[i] : 0;

	double s=0;
	for (int i=0; i < n; i++)
		s += p[i] > 0 ? (p[i]/ntot)*pow(log(p[i]/ntot), 2) : 0;

	return s/(log(2)*log(2));
}

double entropy2(int n, float* p, int norm=0)
{
/*
	overloaded function (float* p --> double s)

	compute second moment of log p.
*/
	double z=1;
	if (norm)
	{
		// normalize probability distribution
		double z=0;
		for (int i=0; i < n; i++) z += p[i] > 0 ? p[i] : 0;
	}

	double s=0;
	for (int i=0; i < n; i++)
		s += p[i] > 0 ? (p[i]/z)*pow(log(p[i]/z), 2) : 0;

	return s/(log(2)*log(2));
}

double entropy2(int n, double* p, int norm=0)
{
/*
	overloaded function (double* p --> double s)

	compute second moment of log p.
*/
	double z=1;
	if (norm)
	{
		// normalize probability distribution
		double z=0;
		for (int i=0; i < n; i++) z += p[i] > 0 ? p[i] : 0;
	}

	double s=0;
	for (int i=0; i < n; i++)
		s += p[i] > 0 ? (p[i]/z)*pow(log(p[i]/z), 2) : 0;

	return s/(log(2)*log(2));
}

charplus get_chars(FILE* f)
{
/*
	read from file f, except skips lines whose first element is #.
	delimiter is white space or tab.

	RETURNS:
	charplus z:
	z.x[i][j]	- jth element in line i, in char* format
	z.m[i]		- number of elements on line i
	z.n		- number of lines.  returns 0 if the file
			  does not exist.
*/
	// ---declarations
	struct charplus z;

	// ---check that file exists.  if not, return with z.n = 0.
	//    reserve space just in case.
	if (!f)
	{
		z.n=0;
		z.m = new int[1];
		z.x = new char**[1];
		z.x[0] = new char*[1];
		z.x[0][0] = new char[1];
		return z;
	}

	// reserve space and make a few declarations
	char* line = new char[MLL_lib];
	int n_elements;
	ostrstream m_ptr;

	strstream inout;
	int line_no=0;
	while (fgets(line, MLL_lib, f))
	{
		inout << line;
		line_no++;
	}

	// ---reserve slightly too much memory, since some lines are
	//    comments.
	z.x = new char**[line_no];
	z.n = 0;
	for (int i=0; i < line_no; i++)
	{
		inout.getline(line, MLL_lib);

		// parse
		int* nw = parse_c(line, n_elements);

		// check that there is at least one element on the line
		if (n_elements == 0)
		{
			z.x[z.n] = new char*[1];
			z.x[z.n][0] = new char[1];
			z.n++;
			m_ptr.write((char*) &n_elements, sizeof(n_elements));
		}

		// check for comment
		else if (line[nw[0]] != '#')
		{
			// ---data exists

			// reserve space
			z.x[z.n] = new char*[n_elements];
			for (int i=0; i < n_elements; i++)
			{
				z.x[z.n][i] =
					new char[length_string(&line[nw[i]])];
				strcpy(z.x[z.n][i], &line[nw[i]]);
			}

			m_ptr.write((char*) &n_elements, sizeof(n_elements));
			z.n++;
		}

		delete [] nw;
	}

	delete [] line;

	z.m = (int*) m_ptr.str();
	return z;
}

float** get_cols(int n_col, int* col, int& n, FILE* f, int silent)
{
/*
	read from file f.  skips lines whose first element is #

	INPUT:
	n_col	- see col below.
	col	- reads columns col[0], col[1], ..., col[n_col-1].
		  skips lines with fewer than max(col) elements.
	f	- file.  in calling program:
		  FILE *f;
		  f = fopen("filename", "r");
		  or
		  f = stdin;
	silent	- if nonzero, do not announce empty lines


	OUTPUT:
	n	- number of columns

	RETURNS:
	x	- x[n][j] = jth element in column col[n]
*/

	// some checking
	int err=0;
	if (n_col < 1)
	{
		cerr << "function get_cols:  first argument is "
			<< "less than 1.  code will probably bomb." << endl;
		err=1;
	}
	for (int i=0; i < n_col; i++)
	{
		if (col[i] < 0)
		{
			cerr << "function get_cols:  col[" << i << "]"
				<< " is less than 0.  code will probably bomb."
				<< endl;
			err=1;
		}
	}
	if (err)
	{
		float** x = new float*[1];
		x[0] = new float[1];
		return x;
	}

	// reserve space
	float** x = new float*[n_col];

	// find maximum column
	int colmax = col[0];
	for (int i=1; i < n_col; i++)
		colmax = col[i] > colmax ?  col[i] : colmax;
	colmax += 1;

	// reserve space and make a few declarations
	char* line = new char[MLL_lib];
	int line_no=0, n_elements;
	ostrstream* x_ptr = new ostrstream[n_col];

	n=0;
	while (fgets(line, MLL_lib, f))
	{
		line_no++;

		// parse
		int* nw = parse_c(line, n_elements);

		// check for empty line
		if (n_elements == 0) {}

		// check for comment
		else if (line[nw[0]] == '#') {}

		// check for enough elements.
		else if (n_elements < colmax)
		{
			if (!silent)
				cerr << "line " << line_no << ":  fewer than "
					<< colmax << " elements; skipping"
					<< endl;
		}
		else
		{
			n++;

			for (int i=0; i < n_col; i++)
			{
				float y = atof(&line[nw[col[i]]]);
				x_ptr[i].write((char*) &y, sizeof(y));
			}
		}

		delete [] nw;
	}

	for (int i=0; i < n_col; i++) x[i] = (float*) x_ptr[i].str();

	delete [] line;
	fclose(f);

	return x;
}

int** get_cols(
	FILE* f, int& lines, int& n_elements, int flag, int min_elements)
{
/*
	read from file f, except skips lines whose first element is #.

	INPUT:
	f		- file to be read from.
	flag		- determines orientation of array x; see below.

	OPTIONAL INPUT:
	min_elements	- only lines with at least min_elements will
			  be saved. default=0

	OUTPUT:
	lines		- number of lines
	n_elements	- number of elements on each line.

	RETURNS:
	x[i][j]		- flag=1:  element i, line j.
			  flag=0:  line i, element j.
*/
	int** x;

	// ---check that file exists.  if not, return with lines =
	//    n_elements=0;
	if (!f)
	{
		lines=n_elements=0;
		return x;
	}

	// ---dummy variable; not used.
	char** dat;

	ostrstream x_ptr;
	ostrstream m_ptr;
	parse_file(f, x_ptr, m_ptr, lines, 0, dat, -1.0, ' ', min_elements);

	// ---check
	if (lines == 0)
	{
		n_elements=0;
		return x;
	}

	// ---number of elements on each line
	int* m = (int*) m_ptr.str();

	// ---check that m[i] is the same for all lines.  if it isn't,
	//    return with lines=-1.
	n_elements = m[0];
	for (int i=1; i < lines; i++) if (m[i] != n_elements)
	{
		lines=-1;
		return x;
	}

	// ---temporarily store the whole file
	double* xtmp = (double*) x_ptr.str();

	// ---transfer to x

	int cnt=0;
	if (flag)
	{
		// ---x[i][j] is element i, line j
		x = cint(n_elements, lines);

		for (int i=0; i < lines; i++) for (int j=0; j < n_elements; j++)
			x[j][i] = (int) floor(xtmp[cnt++] + 0.5);
	}
	else
	{
		// ---x[i][j] is line i, element j
		x = cint(lines, n_elements);

		for (int i=0; i < lines; i++) for (int j=0; j < n_elements; j++)
			x[i][j] = (int) floor(xtmp[cnt++] + 0.5);
	}

	delete [] xtmp;
	fclose(f);

	return x;
}

float** get_cols(
	FILE* f, int& lines, int& n_elements, float flag, int min_elements)
{
/*
	read from file f, except skips lines whose first element is #.

	INPUT:
	f		- file to be read from.
	flag		- determines orientation of array x; see below.

	OPTIONAL INPUT:
	min_elements	- only lines with at least min_elements will
			  be saved. default=0

	OUTPUT:
	lines		- number of lines
	n_elements	- number of elements on each line.

	RETURNS:
	x[i][j]		- flag=1:  element i, line j.
			  flag=0:  line i, element j.
*/
	float** x;

	// ---check that file exists.  if not, return with lines =
	//    n_elements=0;
	if (!f)
	{
		lines=n_elements=0;
		return x;
	}

	// ---dummy variable; not used.
	char** dat;

	ostrstream x_ptr;
	ostrstream m_ptr;
	parse_file(f, x_ptr, m_ptr, lines, 0, dat, -1.0, ' ', min_elements);

	// ---check
	if (lines == 0)
	{
		n_elements=0;
		return x;
	}

	// ---number of elements on each line
	int* m = (int*) m_ptr.str();

	// ---check that m[i] is the same for all lines.  if it isn't,
	//    return with lines=-1.
	n_elements = m[0];
	for (int i=1; i < lines; i++) if (m[i] != n_elements)
	{
		lines=-1;
		return x;
	}

	// ---temporarily store the whole file
	double* xtmp = (double*) x_ptr.str();

	// ---transfer to x

	int cnt=0;
	if (flag)
	{
		// ---x[i][j] is element i, line j
		x = newfloat(n_elements, lines);

		for (int i=0; i < lines; i++) for (int j=0; j < n_elements; j++)
			x[j][i] = xtmp[cnt++];
	}
	else
	{
		// ---x[i][j] is line i, element j
		x = newfloat(lines, n_elements);

		for (int i=0; i < lines; i++) for (int j=0; j < n_elements; j++)
			x[i][j] = xtmp[cnt++];
	}

	delete [] xtmp;
	fclose(f);

	return x;
}

double** get_cols(
	FILE* f, int& lines, int& n_elements, double flag, int min_elements)
{
/*
	read from file f, except skips lines whose first element is #.

	INPUT:
	f		- file to be read from.
	flag		- determines orientation of array x; see below.

	OPTIONAL INPUT:
	min_elements	- only lines with at least min_elements will
			  be saved. default=0

	OUTPUT:
	lines		- number of lines
	n_elements	- number of elements on each line.

	RETURNS:
	x[i][j]		- flag=1:  element i, line j.
			  flag=0:  line i, element j.
*/
	double** x;

	// ---check that file exists.  if not, return with lines =
	//    n_elements=0;
	if (!f)
	{
		lines=n_elements=0;
		x = (double**) calloc(1, sizeof(double*));
		x[0] = (double*) calloc(1, sizeof(double));
		return x;
	}

	// ---dummy variable; not used.
	char** dat;

	ostrstream x_ptr;
	ostrstream m_ptr;
	parse_file(f, x_ptr, m_ptr, lines, 0, dat, -1.0, ' ', min_elements);

	// ---check
	if (lines == 0)
	{
		n_elements=0;
		return x;
	}

	// ---number of elements on each line
	int* m = (int*) m_ptr.str();

	// ---check that m[i] is the same for all lines.  if it isn't,
	//    return with lines=-1.
	n_elements = m[0];
	for (int i=1; i < lines; i++) if (m[i] != n_elements)
	{
		lines=-1;
		return x;
	}

	// ---temporarily store the whole file
	double* xtmp = (double*) x_ptr.str();

	// ---transfer to x

	int cnt=0;
	if (flag)
	{
		// ---x[i][j] is element i, line j
		x = newdouble(n_elements, lines);

		for (int i=0; i < lines; i++) for (int j=0; j < n_elements; j++)
			x[j][i] = xtmp[cnt++];
	}
	else
	{
		// ---x[i][j] is line i, element j
		x = newdouble(lines, n_elements);

		for (int i=0; i < lines; i++) for (int j=0; j < n_elements; j++)
			x[i][j] = xtmp[cnt++];
	}

	delete [] xtmp;
	fclose(f);

	return x;
}

floatplus get_cols(FILE* f)
{
/*
	read from file f, except skips lines whose first element is #.
	call:
		floatplus z = get_cols(stdin);
	or (I believe)
		floatplus z = get_cols(fopen("filename", "r"));

	INPUT:
	f	- input file.

	RETURNS:
	floatplus z:
	z.x[i][j]	- jth element in line i
	z.m[i]		- number of elements on line i
	z.n		- number of lines.  returns 0 if the file
			  does not exist.
*/
	// ---declarations
	struct floatplus z;

	// ---check that file exists.  if not, return with z.n = 0.
	//    reserve space just in case.
	if (!f)
	{
		z.n=0;
		z.m = new int[1];
		z.x = new double*[1];
		z.x[0] = new double[1];
		z.ndat=1;
		z.ddat = new double[1];
		return z;
	}

	// reserve space and make a few declarations
	ostrstream x_ptr;
	ostrstream m_ptr;
	char** dat;

	z.ddat = parse_file(f, x_ptr, m_ptr, z.n, 0, dat, 0.0, ' ', 0);

	z.m = (int*) m_ptr.str();

	// ---temporarily store the whole file
	double* x = (double*) x_ptr.str();

	// ---transfer to z.x
	z.x = newdouble(z.n, z.m);

	int cnt=0;
	for (int i=0; i < z.n; i++) for (int j=0; j < z.m[i]; j++)
		z.x[i][j] = x[cnt++];

	delete [] x;
	fclose(f);

	return z;
}

floatplus get_cols(
	FILE* f, int ndat, char** dat, double dflt=-1.0, char sep=' ',
	int min_elements)
{
/*
	read from file f, except skips lines whose first element is #.
	call:
		floatplus z = get_cols(stdin);
	or (I believe)
		floatplus z = get_cols(fopen("filename", "r"));

	INPUT:
	f		- input file.
	ndat		- number of elements in dat.
	dat[i]		- element to be matched to comments; see parse_file
			  and xdat below.
	dflt		- value returned in z.ddat[i] (see below) if
			  there is no match with dat[i].
	sep		- separator; see xdat below.

			  basically, if, say, the input file contains
			  the line "# noise_sheaf_length=17",
			  dat[3]="noise_sheaf_length" and sep = '=',
			  then z.ddat[3] will equal 17.  if the input
			  file contains the line "# noise_sheaf_length
			  17", again dat[3]="noise_sheaf_length", but
			  this time sep = ' ' (the default), then
			  again z.ddat[3] will equal 17.
	min_elements	- only lines with at least min_elements will
			  be saved.

	RETURNS:
	floatplus z:
	z.x[i][j]	- jth element in line i
	z.m[i]		- number of elements on line i
	z.n		- number of lines.  returns 0 if the file
			  does not exist.
	z.ddat		- data from comment; see above.  if dat[i] is
			  not found, z.ddat[i] comes back equal to dflt.
	ndat		- number of elements in z.dat.
*/
	// ---declarations
	struct floatplus z;

	// ---check that file exists.  if not, return with z.n = 0.
	//    reserve space just in case.
	if (!f)
	{
		z.n=0;
		z.m = new int[1];
		z.x = new double*[1];
		z.x[0] = new double[1];
		z.ndat=1;
		z.ddat = new double[1];
		return z;
	}

	// reserve space and make a few declarations
	ostrstream x_ptr;
	ostrstream m_ptr;

	z.ddat = parse_file(f, x_ptr, m_ptr, z.n, ndat, dat, dflt,
		sep, min_elements);

	z.m = (int*) m_ptr.str();

	// ---temporarily store the whole file
	double* x = (double*) x_ptr.str();

	// ---transfer to z.x
	z.x = newdouble(z.n, z.m);

	int cnt=0;
	for (int i=0; i < z.n; i++) for (int j=0; j < z.m[i]; j++)
		z.x[i][j] = x[cnt++];

	delete [] x;
	fclose(f);

	return z;
}

int* header_info(FILE* f, int n, char** s, int dflt=-1)
{
/*
	get info from header in ".dat" files.  assume header has the
	form

	# ...
	...
	# string=n
	...

	where string is a string and n is an integer.  stops searching
	after all strings are found, or if at the first line _not_
	beginning with #.

	INPUT:
	n	- number of strings.
	s[i]	- a list of strings, i = 0, n-1.
	dflt	- if string is not found, then dflt is returned.

	RETURN:
	info[i]	- the info pointed to by string.  e.g., if
		  s[3]="noise_random_seed" and the .dat file had the
		  line "noise_random_seed=5", then info[3] would come
		  back equal to 5.  if noise_random_seed was _not_ in
		  the .dat file, then info[3] would come back set to
		  dflt.
*/
	int* info = init_constant(n, dflt);
	char* line = new char[MLL_lib];

	while (fgets(line, MLL_lib, f))
	{
		// parse
		int n_elements;
		int* nw = parse_c(line, n_elements);

		// ---if there is at least 1 element on the line, check
		//    that first element is '#'.  if not, return.
		if (n_elements > 0) if (line[nw[0]] != '#') return info;

		// only consider lines with 2 or more elements.
		if (n_elements > 1)
		{
			// ---search for equal sign in second element.
			for (int i=0; ; i++)
			{
				if (line[nw[1]+i] == '\0') break;
				else if (line[nw[1]+i] == '=')
				{
					line[nw[1]+i] = '\0';
		
					// ---look for match
					for (int j=0; j < n; j++)
						if (!strcmp(s[j], &line[nw[1]]))
							info[j]=atoi(&line[nw[1]+i+1]);
		
					break;
				}
			}
		}
		
		delete [] nw;

		// ---see if we are all done
		int k;
		for (k=0; k < n; k++) if (info[k] == dflt) break;
		if (k == n) break;
	}

	delete [] line;

	return info;
}

double* parse_file(
	FILE* f, ostrstream& x_ptr, ostrstream& m_ptr, int& lines,
	int ndat, char** dat, double dflt, char sep=' ', int min_elements=0)
{
/*
	parse file

	INPUT:
	f		- input file.
	ndat		- number of elements in dat.
	dat[i]		- element to be matched to comments; see parse_file
			  and xdat below.
	dflt		- value returned in z.ddat[i] (see below) if
			  there is no match with dat[i].
	sep		- separator; see xdat below.

			  basically, if, say, the input file contains
			  the line "# noise_sheaf_length=17",
			  dat[3]="noise_sheaf_length" and sep = '=',
			  then z.ddat[3] will equal 17.  if the input
			  file contains the line "# noise_sheaf_length
			  17", again dat[3]="noise_sheaf_length", but
			  this time sep = ' ' (the default), then
			  again z.ddat[3] will equal 17.
	min_elements	- only lines with at least min_elements will
			  be saved.

	OUTPUT:
	x_ptr	- pointer to data.
	m_ptr	- point to number of elements on each line.

	RETURNS:
	ddat	- see comments in get_cols above.
*/
	char* line = new char[MLL_lib];
	int line_no=0, n_elements;
	lines=0;
	int cnt=0;
	int* cdat = init_constant((ndat < 1 ? 1 : ndat), 0);
	double* ddat = init_constant((ndat < 1 ? 1 : ndat), dflt);

	while (fgets(line, MLL_lib, f))
	{
		line_no++;

		// parse
		int* nw = parse_c(line, n_elements);

		// check that there is at least one element on the line
		if (n_elements == 0)
		{
			if (n_elements >= min_elements)
			{
				lines++;
				m_ptr.write((char*) &n_elements,
					sizeof(n_elements));
			}
		}
		else if (line[nw[0]] != '#')
		{
			// data exists and is not a comment
			for (int i=0; i < n_elements; i++)
			{
				double y = atof(&line[nw[i]]);
				x_ptr.write((char*) &y, sizeof(y));
			}

			lines++;
			m_ptr.write((char*) &n_elements, sizeof(n_elements));
		}
		else if (line[nw[0]] == '#' && cnt < ndat)
		{
			cnt += xdat(ndat, cdat, dat, sep, n_elements, nw, line,
				ddat);
		}

		delete [] nw;
	}

	delete [] line;

	return ddat;
}

int xdat(int ndat, int* cdat, char** dat, char sep,
	int n_line, int* nw, char* line, double* ddat)
{
/*
	see if any of the elements in line appear in dat.

	INPUT:
	ndat	- number of elements in dat.
	cdat[i]	- if 1, i^th element has been found; if 1, it has not.
		  once an element is found, it is not searched for again.
	dat[i]	- i^th elements.
	sep	- separater.  two cases:
		  sep = ' ':	compare each element in line to each
		  		element in dat until a match is found
				with dat[i].  set ddat[i] to next
				element in line.  if no next element,
				set ddat[i] to 0.
		  sep != ' ':	say sep = '='.  compare each element
		  		in line up to '=' until a match is
				found with dat[i].  set ddat[i] to the
				element after '='.  if nothing is
				there, set ddat[i] to zero.
	n_line	- number of elements in line.
	nw	- list of pointers to elements in line.
	line	- contains elements to be compared to dat.  see sep above.

	OUTPUT:
	ddat	- see sep above.

	RETURNS:
	cnt	- 1 if a match is found, 0 if no match.
*/
	int cnt=0;

	// ---loop over elments in line
	for (int i=0; i < n_line; i++)
	{
		int j, found=1;
		if (sep != ' ')
		{
			// ---parse to sep and put in '\0'.
			for (j=0;; j++)
			{
				if (line[nw[i]+j] == '\0')
				{
					found=0;
					break;
				}

				if (line[nw[i]+j] == sep)
				{
					line[nw[i]+j] = '\0';
					found=1;
					break;
				}
			}

		}

		// ---loop over lements in dat
		if (found) for (int n=0; n < ndat; n++)
		{
			if (cdat[n] == 0) if (!strcmp(dat[n], &line[nw[i]]))
			{
				cdat[n]=1;
				cnt=1;
				if (sep == ' ')
				{
					if (i+1 < n_line)
						ddat[n]=atof(&line[nw[i+1]]);
					else ddat[n]=0.0;
				}
				else
				{
					if (line[nw[i]+j+1] != '\0')
						ddat[n]=atof(&line[nw[i]+j+1]);
					else ddat[n]=0.0;
				}
			}
		}
	}
	
	return cnt;
}

int isnumber(char* s, double& x)
{
/*
	If s is a legal number, the number is transfered to x and a 1 is
	returne. Otherwise a 0 is returned.
*/
	// ---the strategy is to cycle through the elements of s and
	//    determine whether each element is legal. what is legal
	//    depends on what has come before, since numbers can
	//    includ things like 1.6e-2.
	int flag=1;
	for (int i=0; i < strlen(s); i++)
	{
		if (flag == 1)
		{
			// ---allowed values are +/-, 0-9, .
			if (s[i] == '+' || s[i] == '-')		flag=2;
			else if (isdigit(s[i]))			flag=3;
			else if (s[i] == '.')			flag=6;
			else					flag=0;
		}
		else if (flag == 2)
		{
			// ---allowed values are 0-9, .
			if (isdigit(s[i]))			flag=3;
			else if (s[i] == '.')			flag=6;
			else					flag=0;
		}
		else if (flag == 3)
		{
			// ---allowed values are e/E, 0-9, .
			if (s[i] == 'e' || s[i] == 'E')		flag=4;
			else if (isdigit(s[i]))			flag=3;
			else if (s[i] == '.')			flag=7;
			else					flag=0;
		}
		else if (flag == 4)
		{
			// ---allowed values are +/-, 0-9
			if (s[i] == '+' || s[i] == '-')		flag=5;
			else if (isdigit(s[i]))			flag=5;
			else					flag=0;
		}
		else if (flag == 5)
		{
			// ---allowed values are 0-9
			if (isdigit(s[i]))			flag=5;
			else					flag=0;
		}
		else if (flag == 6)
		{
			// ---allowed values are e/E, 0-9, .
			if (isdigit(s[i]))			flag=7;
			else					flag=0;
		}
		else if (flag == 7)
		{
			// ---allowed values are e/E, 0-9
			if (s[i] == 'e' || s[i] == 'E')		flag=4;
			else if (isdigit(s[i]))			flag=7;
			else					flag=0;
		}

		if (flag == 0) return 0;
	}

	// --- can't end on flag=2, since this indicates a + or - all alone.
	if (flag == 2) return 0;

	// ---if we've made it this far, s must be a number.
	x=atof(s);
	return 1;
}

double kl_divergence(int n, float* p, float* q)
{
/*
	computes sum_i p[i] log p[i]/q[i].  assumes p and q are
	normalized.  skips terms with q[i]=0.

	compute entropy of a probability distribution, p.
*/
	double kl=0;
	for (int i=0; i < n; i++)
		kl += p[i] > 0 && q[i] > 0 ? p[i]*log(p[i]/q[i]) : 0;

	return kl/log(2);
}

int length_string(char* s)
{
/*
	return length of \0 terminated string
*/
	int i=0;
	while (s[i++] != '\0');
	return i;
}

int max(int n, int* x)
{
/*
	return maximum of x
*/
	if (n == 0) return 0;

	int xmax = x[0];

	for (int i=1; i < n; i++) xmax = x[i] > xmax ? x[i] : xmax;

	return xmax;
}

float max(int n, float* x)
{
/*
	return maximum of x
*/
	if (n == 0) return 0;

	float xmax = x[0];

	for (int i=1; i < n; i++) xmax = x[i] > xmax ? x[i] : xmax;

	return xmax;
}

double max(int n, double* x)
{
/*
	return maximum of x
*/
	if (n == 0) return 0;

	double xmax = x[0];

	for (int i=1; i < n; i++) xmax = x[i] > xmax ? x[i] : xmax;

	return xmax;
}

int min(int n, int* x)
{
/*
	return minimum of x
*/
	if (n == 0) return 0;

	int xmin = x[0];

	for (int i=1; i < n; i++) xmin = x[i] < xmin ? x[i] : xmin;

	return xmin;
}

float min(int n, float* x)
{
/*
	return minimum of x
*/
	if (n == 0) return 0;

	float xmin = x[0];

	for (int i=1; i < n; i++) xmin = x[i] < xmin ? x[i] : xmin;

	return xmin;
}

double min(int n, double* x)
{
/*
	return minimum of x
*/
	if (n == 0) return 0;

	double xmin = x[0];

	for (int i=1; i < n; i++) xmin = x[i] < xmin ? x[i] : xmin;

	return xmin;
}

double* mmult(int n, int m, double** a, double* x)
{
/*
	return a*x.  a is n by m matrix, x is an m-dim vector
*/
	double* y = init_constant(n, (double) 0);
	for (int i=0; i < n; i++) for (int j=0; j < m; j++)
		y[i] += a[i][j]*x[j];

	return y;
}

double** mmult(int n, int m, int l, double** a, double** x)
{
/*
	return a*x.  a is n by m matrix, x is an m by l matrix
*/
	double** y = init_constant(n, l, (double) 0);
	for (int i=0; i < n; i++) for (int j=0; j < l; j++)
		for (int k=0; k < m; k++) y[i][j] += a[i][j]*x[j][k];

	return y;
}

char** cchar(int n1, int n2)
{
	char** x = (char**) calloc(n1, sizeof(char*));
	for (int n=0; n < n1; n++) x[n] = (char*) calloc(n2, sizeof(char));
	return x;
}

int** cint(int n1, int n2)
{
	int** x = (int**) calloc(n1, sizeof(int*));
	for (int n=0; n < n1; n++) x[n] = (int*) calloc(n2, sizeof(int));
	return x;
}

float** cfloat(int n1, int n2)
{
	float** x = (float**) calloc(n1, sizeof(float*));
	for (int n=0; n < n1; n++) x[n] = (float*) calloc(n2, sizeof(float));
	return x;
}

double** cdouble(int n1, int n2)
{
	double** x = (double**) calloc(n1, sizeof(double*));
	for (int n=0; n < n1; n++) x[n] = (double*) calloc(n2, sizeof(double));
	return x;
}

double** cdouble(int n1, int* n2)
{
	double** x = (double**) calloc(n1, sizeof(double*));
	for (int n=0; n < n1; n++) x[n] =
		(double*) calloc(n2[n], sizeof(double));
	return x;
}

void cfree(int n, char** x)
{
	for (int i=0; i < n; i++) free(x[i]);
	free(x);
}

void cfree(int n, int** x)
{
	for (int i=0; i < n; i++) free(x[i]);
	free(x);
}

void cfree(int n, float** x)
{
	for (int i=0; i < n; i++) free(x[i]);
	free(x);
}

void cfree(int n, double** x)
{
	for (int i=0; i < n; i++) free(x[i]);
	free(x);
}

char** newchar(int n1, int n2)
{
	char** x = new char*[n1];
	for (int n=0; n < n1; n++) x[n] = new char[n2];
	return x;
}

int** newint(int n1, int n2)
{
	int** x = new int*[n1];
	for (int n=0; n < n1; n++) x[n] = new int[n2];
	return x;
}

int** newint(int n1, int* n2)
{
	int** x = new int*[n1];
	for (int n=0; n < n1; n++) x[n] = new int[n2[n]];
	return x;
}

float** newfloat(int n1, int n2)
{
	float** x = new float*[n1];
	for (int n=0; n < n1; n++) x[n] = new float[n2];
	return x;
}

float** newfloat(int n1, int* n2)
{
	float** x = new float*[n1];
	for (int n=0; n < n1; n++) x[n] = new float[n2[n]];
	return x;
}

double** newdouble(int n1, int n2)
{
	double** x = new double*[n1];
	for (int n=0; n < n1; n++) x[n] = new double[n2];
	return x;
}

double** newdouble(int n1, int* n2)
{
	double** x = new double*[n1];
	for (int n=0; n < n1; n++) x[n] = new double[n2[n]];
	return x;
}

int* init_constant(int n, int constant)
{
/*
	overloaded function:  initialize n dimensional vector to value constant.
*/
	int* x = new int[n];
	for (int i=0; i < n; i++) x[i] = constant;
	return x;
}

int** init_constant(int n, int m, int constant)
{
/*
	overloaded function:  initialize n by m array to value constant.
*/
	int** x = newint(n, m);
	for (int i=0; i < n; i++) for (int j=0; j < m; j++) x[i][j] = constant;
	return x;
}

int** init_constant(int n, int* m, int constant)
{
/*
	overloaded function:  initialize n by m[i] array to value constant.
*/
	int** x = new int*[n];
	for (int j=0; j < n; j++)
	{
		x[j] = new int[m[j]];
		for (int i=0; i < m[j]; i++) x[j][i] = constant;
	}
	return x;
}

float* init_constant(int n, float constant)
{
/*
	overloaded function:  initialize n dimensional vector to value constant.
*/
	float* x = new float[n];
	for (int i=0; i < n; i++) x[i] = constant;
	return x;
}

float** init_constant(int n, int m, float constant)
{
/*
	overloaded function:  initialize n by m array to value constant.
*/
	float** x = newfloat(n, m);
	for (int i=0; i < n; i++) for (int j=0; j < m; j++) x[i][j] = constant;
	return x;
}

float** init_constant(int n, int* m, float constant)
{
/*
	overloaded function:  initialize n by m[i] array to value constant.
*/
	float** x = new float*[n];
	for (int j=0; j < n; j++)
	{
		x[j] = new float[m[j]];
		for (int i=0; i < m[j]; i++) x[j][i] = constant;
	}
	return x;
}

double* init_constant(int n, double constant)
{
/*
	overloaded function:  initialize n dimensional vector to value constant.
*/
	double* x = new double[n];
	for (int i=0; i < n; i++) x[i] = constant;
	return x;
}

double** init_constant(int n, int m, double constant)
{
/*
	overloaded function:  initialize n by m array to value constant.
*/
	double** x = newdouble(n, m);
	for (int i=0; i < n; i++) for (int j=0; j < m; j++) x[i][j] = constant;
	return x;
}

double** init_constant(int n, int* m, double constant)
{
/*
	overloaded function:  initialize n by m[i] array to value constant.
*/
	double** x = new double*[n];
	for (int j=0; j < n; j++)
	{
		x[j] = new double[m[j]];
		for (int i=0; i < m[j]; i++) x[j][i] = constant;
	}
	return x;
}

float mi(float** p, int imax, int jmax)
{
/*
	compute mutual information:

		M = sum_ij p_ij log(p_ij/px_i*py_j)
		px_i = sum_j p_ij
		py_j = sum_i p_ij
		

	i runs from 0 to imax - 1
	j runs from 0 to jmax - 1
*/

	// compute px and py
	float* px = new float[imax];
	for (int i=0; i < imax; i++)
	{
		px[i] = 0;
		for (int j=0; j < jmax; j++) px[i] += p[i][j];
	}
	float* py = new float[jmax];
	for (int j=0; j < jmax; j++)
	{
		py[j] = 0;
		for (int i=0; i < imax; i++) py[j] += p[i][j];
	}

	// ---mutual information
	float mi=0;
	// H(X)
	for (int i=0; i < imax; i++) mi -= px[i] > 0 ? px[i]*log(px[i]) : 0;

	// H(Y)
	for (int j=0; j < jmax; j++) mi -= py[j] > 0 ? py[j]*log(py[j]) : 0;

	// H(X,Y)
	for (int i=0; i < imax; i++) for (int j=0; j < jmax; j++)
		mi += p[i][j] > 0 ? p[i][j]*log(p[i][j]) : 0;

	mi /= log(2);

	delete [] px;
	delete [] py;

	return mi;
}

int rabs(int x)
{
	return x > 0 ? x : -x;
}

float rabs(float x)
{
	return x > 0 ? x : -x;
}

double rabs(double x)
{
	return x > 0 ? x : -x;
}

unsigned short* read_binary(FILE* f, int& n, short dum)
{
/*
	read binary file of chars.

	INPUT:
	f	- input file.
	dum	- specifies type (overloaded file).

	OUTPUT:
	n	- number of elements in file.  set to 0 if file does
		  not exist.

	RETURNS:
		- elements in file
*/
	// ---check that file exists.
	unsigned short* x;
	if (!f)
	{
		n=0;
		return x;
	}

	// ---find length
	fseek(f, 0, SEEK_END);
	int nbytes = ftell(f);
	n = nbytes/sizeof(short);
	rewind(f);

	// ---get data
	x = new unsigned short[n];
	fread(x, sizeof(short), n, f);

	// ---if dum < 0, reverse order of bytes
	if (dum < 0) for (int i=0; i < n; i++) x[i] = x[i]/256 + 256*(x[i]%256);

	// ---return
	return x;
}

int* read_binary(FILE* f, int& n, int dum)
{
/*
	read binary file of ints.

	INPUT:
	f	- input file.
	dum	- specifies type (overloaded file).

	OUTPUT:
	n	- number of elements in file.  set to 0 if file does
		  not exist.

	RETURNS:
		- elements in file
*/
	// ---check that file exists.
	int* x;
	if (!f)
	{
		n=0;
		return x;
	}

	// ---find length
	fseek(f, 0, SEEK_END);
	int nbytes = ftell(f);
	n = nbytes/sizeof((int) nbytes);
	rewind(f);

	// ---get data
	x = new int[n];
	fread(x, sizeof(x[0]), n, f);

	// ---return
	return x;
}

float* read_binary(FILE* f, int& n, float dum)
{
/*
	read binary file of floats.

	INPUT:
	f	- input file.
	dum	- specifies type (overloaded file).

	OUTPUT:
	n	- number of elements in file

	RETURNS:
		- elements in file
*/
	// ---check that file exists.
	float* x;
	if (!f)
	{
		n=0;
		return x;
	}

	// ---find length
	fseek(f, 0, SEEK_END);
	int nbytes = ftell(f);
	n = nbytes/sizeof((float) nbytes);
	rewind(f);

	// ---get data
	x = new float[n];
	fread(x, sizeof(x[0]), n, f);

	// ---return
	return x;
}

double* read_binary(FILE* f, int& n, double dum)
{
/*
	read binary file of doubles.

	INPUT:
	f	- input file.
	dum	- specifies type (overloaded file).

	OUTPUT:
	n	- number of elements in file

	RETURNS:
		- elements in file
*/
	// ---check that file exists.
	double* x;
	if (!f)
	{
		n=0;
		return x;
	}

	// ---find length
	fseek(f, 0, SEEK_END);
	int nbytes = ftell(f);
	n = nbytes/sizeof((double) nbytes);
	rewind(f);

	// ---get data
	x = new double[n];
	fread(x, sizeof(x[0]), n, f);

	// ---return
	return x;
}

int* read_binary(char* infile, int& n, int dum)
{
/*
	read binary file of ints.

	INPUT:
	infile	- file name
	dum	- specifies type (overloaded file).

	OUTPUT:
	n	- number of elements in file.  set to 0 if file does
		  not exist.

	RETURNS:
		- elements in file
*/
	// ---open file
	FILE *f = fopen(infile, "r");

	// ---check that file exists.
	int* x;
	if (!f)
	{
		n=0;
		return x;
	}

	// ---find length
	fseek(f, 0, SEEK_END);
	int nbytes = ftell(f);
	n = nbytes/sizeof((int) nbytes);
	rewind(f);

	// ---get data
	x = new int[n];
	fread(x, sizeof(x[0]), n, f);

	// ---return
	return x;
}

float* read_binary(char* infile, int& n, float dum)
{
/*
	read binary file of floats.

	INPUT:
	infile	- file name
	dum	- specifies type (overloaded file).

	OUTPUT:
	n	- number of elements in file

	RETURNS:
		- elements in file
*/
	// ---open file
	FILE *f = fopen(infile, "r");

	// ---check that file exists.
	float* x;
	if (!f)
	{
		n=0;
		return x;
	}

	// ---find length
	fseek(f, 0, SEEK_END);
	int nbytes = ftell(f);
	n = nbytes/sizeof((float) nbytes);
	rewind(f);

	// ---get data
	x = new float[n];
	fread(x, sizeof(x[0]), n, f);

	// ---return
	return x;
}

double* read_binary(char* infile, int& n, double dum)
{
/*
	read binary file of doubles.

	INPUT:
	infile	- file name
	dum	- specifies type (overloaded file).

	OUTPUT:
	n	- number of elements in file

	RETURNS:
		- elements in file
*/
	// ---open file
	FILE *f = fopen(infile, "r");

	// ---check that file exists.
	double* x;
	if (!f)
	{
		n=0;
		return x;
	}

	// ---find length
	fseek(f, 0, SEEK_END);
	int nbytes = ftell(f);
	n = nbytes/sizeof((double) nbytes);
	rewind(f);

	// ---get data
	x = new double[n];
	fread(x, sizeof(x[0]), n, f);

	// ---return
	return x;
}

double* read_n(FILE* f, int lines, int col, int& n, int silent=1)
{
/*
	read n lines of a file and return data in kth column.
	read from file f.  skips lines whose first element is #.

	NOTE:  does not check to see that file exists, or that col and
	       lines > 0!!!  do that in the calling program.

	INPUT:
	f	- file.
	lines	- number of lines to read.
	col	- column to extract.
	silent	- if on, don't list missing elements.

	OUTPUT:
	n	- number of lines actually used.

	RETURNS:
	x	- array containing n numbers.
	
*/

	// ---reserve space
	char s[MLL_lib];
	double* x = new double[lines];

	// ---main loop
	n=0;
	while (fgets(s, MLL_lib, f))
	{
		// ---get element from column "col".  ugly, but it works.
		int i=0, isav=0, n_col=0;
		while (s[i])
		{
			if (s[i] == '#') break;
			while (isspace(s[i])) i++;
			if (s[i] == '\0') break;
		
			while (!isspace(s[i]) && s[i] != '\0') i++;
			if (s[i] == '\0') break;
			s[i++] = '\0';
			if (++n_col == col) break;
			isav=i;
		}

		// ---store data
		if (n_col == col) x[n++] = atof(&s[isav]);
		if (n == lines) break;
	}

	return x;
}

plexon read_plexon(FILE* f, FILE* g=NULL, int sort=1)
{
/*
	read file assumed to be in plexon format.  if g != NULL, write
	header to file g.

	sort if sort=1.

*/
	// reserve space and make a few declarations
	char* line = new char[MLL_lib];
	ostrstream t_ptr, c_ptr;

	// ---write and write (if g != NULL) header
	while (fgets(line, MLL_lib, f))
	{
		if (!isdigit(line[0])) {if (g != NULL) fprintf(g, "%s", line);}
		else break;
	}

	// ---now get actual data.  break as soon as a non-digit or
	//    eof is encountered.  don't forget that "line" now contains
	//    the first line of data.
	char* s;
	struct plexon d;
	d.n=0;
	for (;;)
	{
		// ---parse line
		s=line;
		while (isdigit(*s)) s++;
		*s = '\0';
		s++;
		while (isspace(*s)) s++;

		// ---extract code and time
		int code = atoi(line);
		double time = atof(s);

		// ---store in osttstream
		c_ptr.write((char*) &code, sizeof(code));
		t_ptr.write((char*) &time, sizeof(time));
		d.n++;

		// ---exit at eof or non-digit (files can end with comments).
		if (!fgets(line, MLL_lib, f)) break;
		if (!isdigit(line[0])) break;
	}

	// ---save data
	d.c = (int*) c_ptr.str();
	d.t = (double*) t_ptr.str();

	// ---sort data.  use an unsophisticate method, since data
	//    should be almost sorted anyway.
	if (sort)
	{
		// ---sort data.  use an unsophisticate method, since
		//    data should be almost sorted anyway.

		int i=1;
		for (;;)
		{
			if (d.t[i] < d.t[i-1])
			{
				// ---switch and decrement i
				double tmp = d.t[i];
				d.t[i] = d.t[i-1];
				d.t[i-1] = tmp;
				int cmp = d.c[i];
				d.c[i] = d.c[i-1];
				d.c[i-1] = cmp;
				i -= i == 1 ? 0 : 1;
			}
			else i++;
			if (i == d.n) break;
		}
	}

	return d;
}

int readline(ifstream& infile, int nmax, int line, int& nchars, char* buf)
{
/*
	read a line of a file.  a couple things to note:
		1.  buf must be dimensioned to at least nmax+1.  this
		    is so a \0 can be stuck on the end.
		2.  it is important to read the whole line even if the
		    maximum number of characters has been reached.
		    otherwise, on the next call the rest of the
		    current line will be read.

	INPUT:
	infile	- input file
	nmax	- maximum number of characters allowed on line.  see
		  note (1) above.
	line	- line number; used to tell which line(s) have too
		  many characters.  if negative, no error message will
		  be written.

	OUTPUT:
	buf	- character array containing contents of line.
	nchars	- number of characters on a line.

	RETURNS:
	data	- 1 if this line has data one it.
		  0 if there are no more lines in the file.
*/
	char ch;
	int data=1, exceed=0;
	nchars=0;

	for (;;)
	{
		if (!infile.get(ch))
		{
			data = 0;
			break;
		}
		if (ch == '\n')
		{
			buf[nchars] = '\0';
			data = 1;
			break;
		}
		else if (nchars >= nmax) exceed = 1;
		else buf[nchars++] = ch;
	}
	if (exceed && line > 0)
		cerr << "too many characters on line " << line << endl;
//	cout << "readline:   " << buf << endl;
	return data;
}

int* parse_c(char* buf, int& n)
{
/*
	find all instances of whitespace delimited strings.

	INPUT:
	buf	- character array containing line to be parsed.
		  assumed to end in '\0'.

	OUTPUT:
	buf	- '\0' placed at end of every string.
	n	- number of strings in buf.

	RETURNS:
	nw	- strings are contained in buf[nw[i]], i=0, n-1.
*/
	// initialize
	n = 0;

	// skip past any blank spaces
	int i=0;
	while (isspace(buf[i])) i++;

	// variable length beast
	ostrstream i_ptr;

	// parse
	while (buf[i])
	{
		while (isspace(buf[i])) i++;
		if (buf[i] == '\0') break;

		i_ptr.write((char*) &i, sizeof(i));
		n++;

		while (!isspace(buf[i]) && buf[i] != '\0') i++;
		if (buf[i] == '\0') break;
		buf[i++] = '\0';
	}

	int* nw = (int*) i_ptr.str();

	return nw;
}

int parse_comline(int n, char** t, char* s, int inc, int dflt)
{
/*
	parse command line.  looks for match betwen t[i] and s.  two cases:
	inc > 0:  returns t[i+inc] if t[i]=s
	inc = 0:  returns 1 if t[i] = s
	in either case, if there is no match, returns dflt

	INPUT:
	n	- number of elements on command line (argc)
	t	- list of strings on command line (argv)
	s	- string to be matched with t[i]
	inc	- inc > 0:  returns t[i+inc] if t[i]=s
		  inc = 0:  returns 1 if t[i] = s
		  inc < 0:  returns t[i+|inc|] if t[i]=s
		            no error if there is nothing after flag.
	dflt	- returns dflt if no match.

	RETURNS:
		- see inc and dflt above
*/
	int iflag = inc < 0 ? 0 : 1;
	inc = inc < 0 ? -inc : inc;

	int flag;
	for (flag=0; flag < n; flag++) {if (!strcmp(s, t[flag])) break;}
	if (flag == n) return dflt;

	// special case:  return 1 if this is a stand-alone flag.
	if (inc == 0) return 1;

	if (flag+inc < n)
	{
		// ---check for -flag
		int dash=0;
		for (int i=1; i <= inc; i++)
			if (t[flag+i][0] == '-' && isalpha(t[flag+i][1]))
				dash++;
		if (!dash) return atoi(t[flag+inc]);
	}

	// ---return default if inc > 1
	if (inc > 1) return dflt;

	// ---default:  error
	if (iflag) write_err("missing value after ", s, "; inc = ", inc, ".");
	else return dflt;
}

float parse_comline(int n, char** t, char* s, int inc, float dflt)
{
/*
	parse command line.  looks for match betwen t[i] and s.  two cases:
	inc > 0:  returns t[i+inc] if t[i]=s
	inc = 0:  returns 1 if t[i] = s
	in either case, if there is no match, returns dflt

	INPUT:
	n	- number of elements on command line (argc)
	t	- list of strings on command line (argv)
	s	- string to be matched with t[i]
	inc	- inc > 0:  returns t[i+inc] if t[i]=s
		  inc = 0:  returns 1 if t[i] = s
	dflt	- returns dflt if no match.

	RETURNS:
		- see inc and dflt above
*/
	int flag;
	for (flag=0; flag < n; flag++) {if (!strcmp(s, t[flag])) break;}
	if (flag == n) return dflt;

	if (flag < 0) return dflt;

	// special case:  return 1 if this is a stand-alone flag.
	if (inc == 0) return 1;

	if (flag+inc < n)
	{
		// ---check for -flag
		int dash=0;
		for (int i=1; i <= inc; i++)
			if (t[flag+i][0] == '-' && isalpha(t[flag+i][1]))
				dash++;
		if (!dash) return atof(t[flag+inc]);
	}

	// ---return default if inc > 1
	if (inc > 1) return dflt;

	// ---default:  error
	cerr << "missing value after " << s << "; inc = " << inc << "." << endl;
	exit(1);
}

double parse_comline(int n, char** t, char* s, int inc, double dflt)
{
/*
	parse command line.  looks for match betwen t[i] and s.  two cases:
	inc > 0:  returns t[i+inc] if t[i]=s
	inc = 0:  returns 1 if t[i] = s
	in either case, if there is no match, returns dflt

	INPUT:
	n	- number of elements on command line (argc)
	t	- list of strings on command line (argv)
	s	- string to be matched with t[i]
	inc	- inc > 0:  returns t[i+inc] if t[i]=s
		  inc = 0:  returns 1 if t[i] = s
	dflt	- returns dflt if no match.

	RETURNS:
		- see inc and dflt above
*/
	int flag;
	for (flag=0; flag < n; flag++) {if (!strcmp(s, t[flag])) break;}
	if (flag == n) return dflt;

	if (flag < 0) return dflt;

	// special case:  return 1 if this is a stand-alone flag.
	if (inc == 0) return 1;

	if (flag+inc < n)
	{
		// ---check for -flag
		int dash=0;
		for (int i=1; i <= inc; i++)
			if (t[flag+i][0] == '-' && isalpha(t[flag+i][1]))
				dash++;
		if (!dash) return atof(t[flag+inc]);
	}

	// ---return default if inc > 1
	if (inc > 1) return dflt;

	// ---default:  error
	cerr << "missing value after " << s << "; inc = " << inc << "." << endl;
	exit(1);
}

int* parse_comline(int n, char** t, char* s, int& nmatch)
{
/*
	NOTE:  this should be phased out!!!

	parse command line.  looks for match betwen t[i] and s.
	returns all elements after match until next "-" flag.
	n is number of elements.

	INPUT:
	n	- number of elements on command line (argc)
	t	- list of strings on command line (argv)
	s	- string to be matched with t[i]

	OUTPUT:
	nmatch	- number of elements after flag.
		  if nmatch < 0, exit with error if number of elements
		  is not equal to -nmatch.

	RETURNS:
		- array containing elements.
*/
	int target = nmatch < 0 ? -nmatch : 0;

	// ---find match to flag
	int i1;
	for (i1=0; i1 < n; i1++) {if (!strcmp(s, t[i1])) break;}

	// ---no match
	if (i1 == n)
	{
		if (target) write_err(s, " needs ", target, " arguments.");
		nmatch = 0;
		int* x = new int[1];
		return x;
	}

	// ---find next occurrance of "-" followed by non-digit
	int i2;
	for (i2=i1+1; i2 < n; i2++)
	{
		if (t[i2][0] == '-') if (!isdigit(t[i2][1])) break;
	}

	nmatch = i2 - i1 - 1;
	int* x = new int[nmatch > 0 ? nmatch : 1];
	for (int i=0; i < nmatch; i++) x[i] = atoi(t[i+i1+1]);

	if (target && target != nmatch)
		write_err(s, " needs ", target, " arguments.");

	return x;
}

int* parse_comline(int n, char** t, char* s, int dum, int& nmatch)
{
/*
	a bit of ugliness to distinguish this overloaded function from
	float*, etc.:  need an int in second to last argument.

	parse command line.  looks for match betwen t[i] and s.
	returns all elements after match until next "-" flag.
	n is number of elements.

	INPUT:
	n	- number of elements on command line (argc)
	t	- list of strings on command line (argv)
	s	- string to be matched with t[i]
	dum	- see below

	OUTPUT:
	nmatch	- number of elements after flag.
		  if dum < 0, exit with error if number of elements
		  is not equal to -dum.

	RETURNS:
		- array containing elements.
*/
	int target = dum < 0 ? -dum : 0;

	// ---find match to flag
	int i1;
	for (i1=0; i1 < n; i1++) {if (!strcmp(s, t[i1])) break;}

	// ---no match
	if (i1 == n)
	{
		if (target) write_err(s, " needs ", target, " arguments.");
		nmatch = 0;
		int* x = new int[1];
		return x;
	}

	// ---find next occurrance of "-" followed by non-digit
	int i2;
	for (i2=i1+1; i2 < n; i2++)
	{
		if (t[i2][0] == '-') if (!isdigit(t[i2][1])) break;
	}

	nmatch = i2 - i1 - 1;
	int* x = new int[nmatch > 0 ? nmatch : 1];
	for (int i=0; i < nmatch; i++) x[i] = atoi(t[i+i1+1]);

	if (target && target != nmatch)
		write_err(s, " needs ", target, " arguments.");

	return x;
}

int* parse_comline(int n, char** t, char* s, int dum, int& nmatch, char* dflt)
{
/*
	a bit of ugliness to distinguish this overloaded function:
	need an int in third to last argument.

	parse command line.  looks for match betwen t[i] and s.
	returns all elements after match until next "-" flag.  n is
	number of elements; i.e., the size of t.

	INPUT:
	n	- number of elements on command line (argc)
	t	- list of strings on command line (argv)
	s	- string to be matched with t[i]
	dum	- < 0:  exit with error if number of elements is not
			equal to -dum.  this feature turned off if
			dflt is present.
		- > 0:  use default string (dflt) to fill in any
		        missing values, up to dum elements.
	dflt	- string containing defaults. if empty, none will be applied.

	OUTPUT:
	nmatch	- number of elements after flag.

	RETURNS:
		- array containing elements.
*/
	int ndum = (int) dum;
	int target = ndum < 0 ? -ndum : 0;

	// ---find match to flag
	int i1;
	for (i1=0; i1 < n; i1++) {if (!strcmp(s, t[i1])) break;}

	// ---no match
	if (ndum <= 0 && i1 == n)
	{
		if (target) write_err(s, " needs ", target, " arguments.");
		nmatch = 0;
		int* x = new int[1];
		return x;
	}

	// ---find next occurrance of "-" followed by non-digit and
	//    non-decimal point.
	int i2;
	for (i2=i1+1; i2 < n; i2++)
	{
		if (t[i2][0] == '-')
			if (!isdigit(t[i2][1]) && t[i2][1] != '.') break;
	}

	// ---compute length of array.
	nmatch = i2 - i1 - 1;
	int* x = new int[1 + (nmatch > ndum ? nmatch : ndum)];

	// ---set elements
	for (int i=0; i < nmatch; i++) x[i] = atoi(t[i+i1+1]);

	// ---defaults
	if (nmatch < ndum)
	{
		double* tmp = s_parse(ndum, dflt);
		for (int i=nmatch; i < ndum; i++) x[i] = (int) (0.5+tmp[i]);
		free(tmp);
	}

	// ---check that we we have hit the target number of elements,
	//    target=0 if dum > 0.
	if (target && target != nmatch)
		write_err(s, " needs ", target, " arguments.");

	return x;
}

float* parse_comline(
	int n, char** t, char* s, float dum, int& nmatch, char* dflt)
{
/*
	a bit of ugliness to distinguish this overloaded function from
	int*:  need a (float) 1 (or some such double) in third to last argument.

	parse command line.  looks for match betwen t[i] and s.
	returns all elements after match until next "-" flag.  n is
	number of elements; i.e., the size of t.

	INPUT:
	n	- number of elements on command line (argc)
	t	- list of strings on command line (argv)
	s	- string to be matched with t[i]
	dum	- < 0:  exit with error if number of elements is not
			equal to -dum.  this feature turned off if
			dflt is present.
		- > 0:  use default string (dflt) to fill in any
		        missing values, up to dum elements.
	dflt	- string containing defaults. if empty, none will be applied.

	OUTPUT:
	nmatch	- number of elements after flag.

	RETURNS:
		- array containing elements.
*/
	int ndum = (int) dum;
	int target = ndum < 0 ? -ndum : 0;

	// ---find match to flag
	int i1;
	for (i1=0; i1 < n; i1++) {if (!strcmp(s, t[i1])) break;}

	// ---no match
	if (ndum <= 0 && i1 == n)
	{
		if (target) write_err(s, " needs ", target, " arguments.");
		nmatch = 0;
		float* x = new float[1];
		return x;
	}

	// ---find next occurrance of "-" followed by non-digit and
	//    non-decimal point.
	int i2;
	for (i2=i1+1; i2 < n; i2++)
	{
		if (t[i2][0] == '-')
			if (!isdigit(t[i2][1]) && t[i2][1] != '.') break;
	}

	// ---compute length of array.
	nmatch = i2 - i1 - 1;
	float* x = new float[1 + (nmatch > ndum ? nmatch : ndum)];

	// ---set elements
	for (int i=0; i < nmatch; i++) x[i] = atof(t[i+i1+1]);

	// ---defaults
	if (nmatch < ndum && strcmp(dflt, ""))
	{
		double* tmp = s_parse(ndum, dflt);
		for (int i=nmatch; i < ndum; i++) x[i] = tmp[i];
		free(tmp);
		nmatch=ndum;
	}

	// ---check that we we have hit the target number of elements,
	//    target=0 if dum > 0.
	if (target && target != nmatch)
		write_err(s, " needs ", target, " arguments.");

	return x;
}

double* parse_comline(
	int n, char** t, char* s, double dum, int& nmatch, char* dflt)
{
/*
	a bit of ugliness to distinguish this overloaded function from
	int*:  need a 1.0 (or some such double) in second to last argument.

	parse command line.  looks for match betwen t[i] and s.
	returns all elements after match until next "-" flag.
	n is number of elements.

	INPUT:
	n	- number of elements on command line (argc)
	t	- list of strings on command line (argv)
	s	- string to be matched with t[i]
	dum	- < 0:  exit with error if number of elements is not
			equal to -dum.  this feature turned off if
			dflt is present.
		  > 0:  use default string (dflt) to fill in any
		  	missing values, up to dum elements.
	dflt	- string containing defaults.  if empty, none will be applied.

	OUTPUT:
	nmatch	- number of elements after flag.

	RETURNS:
		- array containing elements.
*/
	int ndum = (int) dum;
	int target = ndum < 0 ? -ndum : 0;

	// ---find match to flag
	int i1;
	for (i1=0; i1 < n; i1++) {if (!strcmp(s, t[i1])) break;}

	// ---no match
	if (i1 == n)
	{
		if (target) write_err(s, " needs ", target, " arguments.");
		nmatch = 0;
		double* x = new double[1];
		return x;
	}

	// ---find next occurrance of "-" followed by non-digit and
	//    non-decimal point.
	int i2;
	for (i2=i1+1; i2 < n; i2++)
	{
		if (t[i2][0] == '-')
			if (!isdigit(t[i2][1]) && t[i2][1] != '.') break;
	}

	// ---compute length of array.
	nmatch = i2 - i1 - 1;
	double* x = new double[nmatch > 0 ? nmatch : 1];

	// ---set elements
	for (int i=0; i < nmatch; i++) x[i] = atof(t[i+i1+1]);

	// ---defaults
	if (nmatch < ndum)
	{
		double* tmp = s_parse(ndum, dflt);
		for (int i=nmatch; i < ndum; i++) x[i] = tmp[i];
		free(tmp);
	}

	// ---check that we we have hit the target number of elements
	//    target=0 if dum > 0.
	if (target && target != nmatch)
		write_err(s, " needs ", target, " arguments.");

	return x;
}

char* parse_comline(int n, char** t, char* s, int inc, char* dflt)
{
/*
	parse command line.  looks for match betwen t[i] and s.  two cases:
	inc > 0:  returns t[i+inc] if t[i]=s
	inc = 0:  returns 1 if t[i] = s
	in either case, if there is no match, returns dflt

	INPUT:
	n	- number of elements on command line (argc)
	t	- list of strings on command line (argv)
	s	- string to be matched with t[i]
	inc	- inc > 0:  returns t[i+inc] if t[i]=s
		  inc = 0:  returns 1 if t[i] = s
	dflt	- returns dflt if no match.

	RETURNS:
		- see inc and dflt above
*/
	int flag;
	for (flag=0; flag < n; flag++) {if (!strcmp(s, t[flag])) break;}
	if (flag == n) return dflt;

	// special case:  return 1 if this is a stand-alone flag.
	if (inc == 0) return "1";

	if (flag+inc < n)
	{
		// ---check for -flag
		int dash=0;
		for (int i=1; i <= inc; i++) if (t[flag+i][0] == '-') dash++;
		if (!dash) return t[flag+inc];
	}

	// ---return default if inc > 1
	if (inc > 1) return dflt;

	// ---default:  error
	cerr << "missing value after " << s << "; inc = " << inc << "." << endl;
	exit(1);
}

double* s_parse(int n, char* v)
{
/*
	simple string parse:  turn the first n elements of v into
	doubles and store in x.  if there are fewer than n elements,
	repeat the last element that occurs.  if there are no
	elements, fill x with zeros.
*/
	int cnt=0, i=0, j=0;
	char *s, *t = (char*) calloc(1024, sizeof(char));
	double* x = (double*) calloc(n, sizeof(double));

	strcpy(t, v);
	s=t;

	for (;;)
	{
		while (*s && isspace(*s)) {s++; j++; i++;}
		while (*s && !isspace(*s)) {s++; j++;}
		if (i == j) break;

		*s='\0';
		x[cnt] = atof(&t[i]);
		cnt++;

		j++;
		s++;
		i=j;

		if (cnt == n) break;
	}

	// ---fill rest of array
	if (cnt == 0) {x[0]=0.0; cnt++;}
	if (cnt < n) for (int i=cnt; i < n; i++) x[i] = x[cnt-1];

	free(t);
	return x;
}

int unit_no(char* s)
{
/*
	get unit number.  three possibilities:
	1.  s is an integer.  return atoi(s).
	2.  s is an integer followed by a, b, c or d.  return
	    4*(n-1)+f(alpha) where n is the integer part and
	    f(a)=1, f(b)=2, f(c)=3, f(d)=4.
	3.  s is none of the above.  return 0.
*/
	int n = strlen(s);

	// ---check that first part of s is OK.
	for (int i=0; i < n-1; i++) if (!isdigit(s[i])) return 0;

	// ---if last character in s is digit, return atoi(s).
	if (isdigit(s[n-1])) return atoi(s);
	else
	{
		// ---translate last character via f(alpha) above; if
		//    last character is not a-d, return 0.
		int offset;
		if (s[n-1] == 'a') offset=1;
		else if (s[n-1] == 'b') offset=2;
		else if (s[n-1] == 'c') offset=3;
		else if (s[n-1] == 'd') offset=4;
		else return 0;

		// ---translate from plexon to unit format and return result.
		char* dum = new char[n];
		for (int i=0; i < n-1; i++) dum[i] = s[i];
		dum[n-1] = '\0';
		int cell_unit = 4*(atoi(dum)-1)+offset;

		delete [] dum;
		return cell_unit;
	}
}

void write_err(char* s)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s\n", s);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s, int n)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%d\n", s, n);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s1, int n, char* s2)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%d%s\n", s1, n, s2);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s, double z)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%f\n", s, z);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s1, char* s2)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%s\n", s1, s2);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s1, char* s2, char* s3)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%s%s\n", s1, s2, s3);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s1, char* s2, int n, char* s3)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%s%d%s\n", s1, s2, n, s3);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s1, char* s2, int n1, char* s3, int n2)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%s%d%s%d\n", s1, s2, n1, s3, n2);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s1, char* s2, char* s3, int n, char* s4)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%s%s%d%s\n", s1, s2, s3, n, s4);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s1, char* s2, char* s3, int n1, char* s4, int n2)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%s%s%d%s%d\n", s1, s2, s3, n1, s4, n2);
	fprintf(stderr, "");
	exit(1);
}

void write_err(char* s1, char* s2, char* s3, char* s4, char* s5)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%s%s%s%s\n", s1, s2, s3, s4, s5);
	fprintf(stderr, "");
	exit(1);
}

void write_warning(char* s)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s\n", s);
	fprintf(stderr, "");
}

void write_warning(char* s1, char* s2, char* s3)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%s%s\n", s1, s2, s3);
	fprintf(stderr, "");
}

void write_warning(char* s1, int n, char* s2)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%d%s\n", s1, n, s2);
	fprintf(stderr, "");
}

void write_warning(char* s1, int n1, char* s2, int n2, char* s3)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%d%s%d%s\n", s1, n1, s2, n2, s3);
	fprintf(stderr, "");
}

void write_warning(char* s1, char* s2, char* s3, int n, char* s4)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%s%s%d%s\n", s1, s2, s3, n, s4);
	fprintf(stderr, "");
}

void write_warning(char* s1, int n1, char* s2, char* s3, int n2, char* s4)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%d%s%d%s\n", s1, n1, s2, n2, s3);
	fprintf(stderr, "");
}

void write_warning(char* s1, double x, char* s2)
{
	fprintf(stderr, "");
	fprintf(stderr, "%s%f%s\n", s1, x, s2);
	fprintf(stderr, "");
}

void write_message(char* s, char* t)
{
FILE* f = fopen(s, "w");

fprintf(f, t);
fclose(f);

char sys[80];
sprintf(sys, "more %s", s);
system(sys);

sprintf(sys, "rm %s", s);
system(sys);
}

void zfree(floatplus& z)
{
	for (int i=0; i < z.n; i++) delete [] z.x[i];
	delete [] z.x;
	delete [] z.m;
}