#!/usr/bin/perl
################################################################################
#
# sec_calc.pl = Input and calculate scores for size exclusion chromatograms 
#		for PROSPERO
#
# Revision date: 21 May 2012
# Authors:	Frank Zucker (Jul. 2010), 
#		Hae Young Kim (Sept. 2010)
#		Bio-Molecular Structure Center
#		Dept. of Biochemistry, 
#		University of Washington, Seattle
#
# Input module and gaussian curve fitter for SEC curves for PROSPERO,
# Prediction of outcome from Sequence and Protein Experimental Results Online
#
# Fit Gaussians to .asc curve and doc files from PrimeView, write to XML.
#
################################################################################
#
# Input: Absorbance vs time or volume curve (e.g. file XXXXXcurve.asc),
#	 Metadata = flow rate, fraction size, start and end times (e.g. Xdoc.asc),
#		start and end of pooled fractions (e.g. from pooled_fractions.csv)
#
# Output: XML, <SEC_DATA> around <DERIVED_VALUES>, <RAW_DATA> & <MODEL_VALUES>:
#	<SEC_DATA> has attributes:
#		BATCH (sample identification)
#		CALCULATED,	 VERSION (sec_calc time stamp and version #)
#			and SEC run metadata including: 
#		FLOWRATE,	FRACTION_VOLUME, 
#		STARTTIME,	ENDTIME,
#		POOL_FROM_FRACTION,	POOL_TO_FRACTION, 
#		STARTDATE,	ENDDATE
#	<DERIVED_VALUES> has attributes:
#			AREA	= total of observed - background
#			R1	= Rabs for 1 Gaussian
#			PP1	= Pool purity for 1 Gaussian
#			DIRECTORY = gnuplot intermediate file folder
#		as well as a csv table of npeaks, Rabs, RMSD, and max PP for
#		1 peak, for the "best" model, and in verbose mode for all models.
#		It also contains elements:
#			<BEST_MODEL> and, for verbose mode, 
#			<MODEL> for each model tried.
#		Both <BEST_MODEL> and <MODEL> elements have attribute:
#			NPEAKS  = number of Gaussians in model
#			RABS	= sum of |obs-calc| / sum of (obs - bg)
#			RMSD	= sqrt(sum of (obs-calc)^2 / N points)
#			PP	= Purity of pooled fractions
#		<BEST_MODEL> contains attribute PEAKMIN and a csv table with headers:
#			Peak Number = order by peak start position,
#			Start, End = center -/+ 2*sigma, covering 95% of peak area
#			Fraction of Total (observed area), and 
#			Fraction of Pooled (area).
#		with one row for each peak plus totals for areas (should be 1.000)
#
#		The "best" model is the one with the lowest BIC value of those
#		which have a change in Rabs >= peakmin (default = 5% i.e. 0.05).
#		The program tries up to peak_max (now = 7) peaks, stopping when
#		the residual change falls below peakmin.
#		Other criteria which could be (but are not currently) applied:
#	 		maximum residual < peakmin * maximum observed absorbance,
#			or change in max. resid. < peakmin * max. resid, 
#			or minimum peak area < peakmin * total area,
#
#		<BEST_MODEL> & <MODEL> elements contain <BACKGROUND> and <PEAK>s
#			<BACKGROUND> has attrubutes INTERCEPT and SLOPE
#				where background = INTERCEPT + x*SLOPE
#				Later versions may include attribute SQUARE,
#				where bg = INTERCEPT + x*SLOPE + x^2*SQUARE
#			<PEAK> has attributes PEAKNUM, CENTER, HEIGHT, SIGMA, AREA,
#				where y = HEIGHT*e^(-(x-CENTER)^2/2*SIGMA^2)
#	<RAW_DATA> = csv with min, mAu (observed)
#	<MODEL_VALUES> = csv with fraction, absorbance (observed), pool, models
#		where frxn = (time - start time)*flow/fraction size,
#		pool	= 0 or max. abs. in pooled fractions
#		and models = sets of columns for G1, G2, ... G<max>; for each Gn,
#		Gn	= total for n-Gaussian model (bg + n peaks), 
#		Gn_res	= residual = obs. abs - Gn total, 
#		Gn_bg	= background values
#		Gn_1, Gn_2, ... Gn_n = Gaussian curves 1 to n for n-peak model
#		
#		Note 1: each Gn starts at column (5n + n^2)/2 counting from 1
#		Note 2: to plot an individual peak m you'll want Gn_bg + Gn_m,
#			i.e. for c=(5n+n*n)/2, you want #(c+2) + #(c+2+m)
#
################################################################################
#
# Main wrapper: accept batch number or file in, file out and fractions pooled
# Pass to sec_parse to pull out data and pass it to sec_derive, which uses sec_cif


use strict;
use warnings;
use Params::Validate qw(:all);

use PDL;
use IO::File;
use IO::CaptureOutput qw(capture);
use File::Copy;
use Data::Dumper;
use File::Path qw(make_path remove_tree);
use DateTime;
my $GNUPLOT = "gnuplot";
my $GNUPLOT_SERVER = "/usr/local/bin/gnuplot";
$GNUPLOT = $GNUPLOT_SERVER if -e $GNUPLOT_SERVER;

{
my $version_ID = "sec_calc_1.09";
my $script = substr($0, 1 + CORE::rindex($0,($0 =~ /\\/) ? "\\" : "/"));

my $help_text = "\n\tSEC Calculator $version_ID: fit Gaussians to SEC curves.\n".
"USAGE:\n".
"$script -b <batch> [<path>]   |   <in_file>   |   - \n".
"    [-o <out_file>]   [-p peakmin]   [-n npeaks]\n".
"    [-f <from_frxn>] [-t <to_frxn>] [-r|--raw] [-m|--model] [-g|gnudir tmpdir]\n".
"    [-s|--slope] [-v|--verbose] [-d|--debug] [--version] [-h | --help]\n".
"where:\n".
	"\tItems in <angle brackets> should be replaced by the real thing,\n".
	"\tItems in [square brackets] are optional. \"a | b\" means a or b.\n".
	"\t-h, --help, or no input specifier prints this text.\n".
"****  You must supply a batch, input file or - meaning STANDARD INPUT  ****\n".
"If you give <in_file> as [<path>/]<batch>curve.asc, or give -b <batch> [<path>],\n".
"then it looks for 3 files from <path> if given, or from the current directory:\n".
"  - <batch>curve.asc     = absorbance vs time or volume from PrimeView),\n".
"  - <batch>doc.asc       = fraction volume and other metadata from PrimeView\n".
"  - pooled_fractions.csv = text file with comma-separted values you create with\n".
"                           a header line like \"batch,from,to\"\n".
"Peak purity is calculated if <from_frxn> and <to_frxn> are given on the command\n".
"line or read from pooled_fractions.csv. If <batch>doc.asc file is not found,\n".
"then metadata defaults are fraction volume=5.0 mL, flow rate=0.5 mL/min, \n".
"start time = 160.33 min.  Other defaults are:\n".
	"\t<out_file> ... \t[<batch>|<in_file>].sec.csv if either is given,\n".
		"\t\t\tor STANDARD OUTPUT (e.g. screen) if not\n".
	"\t<peakmin> ... \t5, peak finder stops if next peak <5% of total area\n".
	"\t<npeaks> ... \t7 peaks maximum\n".
	"   NO\t-r or --raw\tDo NOT write raw data as standardized csv\n".
	"   NO\t-m or --model\tDo NOT write model values as standardized csv\n".
	"   NO\t-g or --gnudir\tDo NOT read or write gnuplot intermediate files\n".
	"   NO\t-v or --verbose\tDo NOT write peaks and residuals for all models\n".
	"   NO\t-s or --slope\tDo NOT limit the background to a constant value\n".
"Default output: peak table, residual values & pool purity for best model only.\n";

# Parse command line options

my ($batch, $infile, $outfile, $peakmin, $fromfrxn, $tofrxn, $raw, $model, 
	$verbose, $debug, $gnudir,$slopeflag,$width, $npeaks, $plot) = 
	get_command_line_args($version_ID, $help_text);

# $verbose = $debug if ($debug && !$verbose);
$peakmin = 5 if ! $peakmin;
print  "Read from ", $infile || "STD IN",
	" to ", $outfile || "STD OUT","\n" if ($verbose);

# Pass infile, outfile etc. to parsing routine

my $result = sec_parse(	{	batch	=> $batch,
				infile	=> $infile, 
				outfile	=> $outfile, 
				peakmin	=> $peakmin,
				fromfrxn=> $fromfrxn,
				tofrxn	=> $tofrxn,
				raw	=> $raw,
				model	=> $model,
				verbose	=> $verbose, 
				debug	=> $debug,
				version	=> $version_ID, 
				gnudir	=> $gnudir,
				slopeflag => ! $slopeflag,
				width	=> $width,
				npeaks	=> $npeaks,
				plot	=> $plot,
			} );

print "Result $result" if ($verbose);
print STDERR "ERROR IN PARSING SEC" if (!$result);
exit(!$result);

}

################################################################################
#
# get_command_line_args = generic command to get input module command lines
#
# In: non-standard option(s) i.e. other than file in, out, verbose, debug,
#						version and help.
#	 version string, help string
# Out: input and output file names, non-standard parameter, verbose, debug
#
# Exits after printing version, help, or command-line parsing error.
#
################################################################################

sub get_command_line_args {
	my ($version_ID, $help_text) = @_;
	use Getopt::Long;
	Getopt::Long::Configure("bundling", "ignorecase_always");
	my $batch	= 0;
	my $infile	= 0;
	my $outfile	= 0;
	my $peakmin	= 0;
	my $fromfrxn	= 0;
	my $tofrxn	= 0;
	my $raw		= 0;
	my $model	= 0;
	my $verbose	= 0;
	my $debug	= 0;
	my $gnudir	= 0;
	my $version	= 0;
	my $help	= 0;
	my $slopeflag	= 0;
	my $width 	= 0;
	my $npeaks	= 0;
	my $plot	= 0;
	my $result = GetOptions("in|i=s"	=> \$infile,
				"out|o=s"	=> \$outfile,
				"batch|b=s"	=> \$batch,
				"peakmin|p=f"	=> \$peakmin,
				"fromfrxn|f=i"	=> \$fromfrxn,
				"tofrxn|t=i"	=> \$tofrxn,
				"slope|s"	=> \$slopeflag,
				"width|w=f"	=> \$width,
				"raw|r"		=> \$raw,
				"model|m"	=> \$model,
				"verbose|v+"	=> \$verbose,
				"debug|d+"	=> \$debug,
				"gnudir|g=s"	=> \$gnudir,
				"npeaks|n=i"	=> \$npeaks,
				"chart|c=s"	=> \$plot,
				"version"	=> \$version,
				"help|h"	=> \$help);
#	print STDERR "ARGV: |@ARGV|" if $debug;
#	print STDERR "Infile: |$infile|" if $debug;
	$infile = shift @ARGV if (@ARGV && !$infile);
	if ($version) 
	{
		print("Version $version_ID\n");
		exit(0);
	}
	if ($help || !$result || !$infile)
	{
		print $help_text;
		exit(!$result);
	}
	return ($batch, $infile, $outfile, $peakmin, $fromfrxn, $tofrxn, $raw, 
		$model, $verbose, $debug, $gnudir, $slopeflag, $width, $npeaks, $plot);
}

################################################################################
#
# parse_sec = parser for SEC files for PROSPERO SEC input/fit module
#
# In:	batch and/or input filename (default: standard in), 
#	output filename (default: standard out),
#	peakmin = cutoff for finding residual peaks (default: 5%)
#	from and to fractions
#	raw and model flags
#	verbose and debug flags (default none)
# Out: write derived, raw and model data to output file
# Return: 1 on success, 0 on failure
#
################################################################################

sub sec_parse {
	
	my %params = validate( @_, {	
					batch	=> { type => SCALAR, default=>0},
					infile	=> { type => SCALAR, default=>'-'},
					outfile => { type => SCALAR, default=>'-'},
					peakmin	=> { type => SCALAR, default=>5},
					fromfrxn=> { type => SCALAR, default=>0},
					tofrxn	=> { type => SCALAR, default=>0},
					raw	=> { type => SCALAR, default=>0},
					model	=> { type => SCALAR, default=>0},
					verbose	=> { type => SCALAR, default=>0},
					debug	=> { type => SCALAR, default=>0},
					version => { type => SCALAR, 
						default => "sec_calc_v1.00" },
					gnudir	=> { type => SCALAR, default=>0},
					slopeflag=> { type => SCALAR, default=>0},
					width	=> { type => SCALAR, default=>0},
					npeaks	=> { type => SCALAR, default=>7},
					plot	=> { type => SCALAR, default=>0},
				   } );
	my $batch	= $params{'batch'};
	my $infile	= $params{'infile'};
	my $outfile	= $params{'outfile'};
	my $peakmin	= $params{'peakmin'};
	my $fromfrxn	= $params{'fromfrxn'};
	my $tofrxn	= $params{'tofrxn'};
	my $print_raw	= $params{'raw'};
	my $print_model	= $params{'model'};
	my $verbose	= $params{'verbose'};
	my $debug	= $params{'debug'};
	my $version_ID	= $params{'version'};
	my $gnudir	= $params{'gnudir'};	
	my $slopeflag	= $params{'slopeflag'};	
	my $width	= $params{'width'};	
	$width = 2.0 if ! $width;
	my $npeaks	= $params{'npeaks'};	
	$npeaks = 7 if ! $npeaks;
	my $plot	= $params{'plot'};	
	my $gnuin	= 0;
	my $gnuout	= 0;
	# Read from existing gnu directory, or write to non-existant directory
	if ($gnudir) {
		if (-e $gnudir) {
			$gnuin = $gnudir;
			print "Read from gnuplot intermediate directory $gnuin\n"
				if $verbose;
		} else {
			$gnuout = $gnudir;
			print "Write to gnuplot intermediate directory $gnuout\n"
				if $verbose;
		}
	}
	# print "**** gnuin: $gnuin  gnuout: $gnuout\n";
	
	my $sourcefile = '';


	print  "Batch $batch, In $infile, Out $outfile\n" if $debug;

	# Constants:	c2f = c (sigma) to FWHM = 2*sqrt(2 ln 2)
	#		sq2p = sqrt(2 pi), 
	#	where pi = 4*arctan(angle from 0,0 to 1,1), and a*c*sq2p = area
	my $c2f = 2.0 * sqrt( 2 * log(2.0) );
	my $sq2p = sqrt( 8 * atan2(1, 1) );

	# Open input file if needed 

	my $delim = 0;
	my $ncols = 0;
	my $headline = 0;
	my $header1  = 'min';
	my $fixedcol = -1;
	my $datacol  = -1;
	my $flowrate = 0.5;	# default, in mL / min
	my $frxnvol  = 5.0;	# default, in mL per fraction
	my $starttime = 160.33; # default, in min
	my $endtime  = 1000.0; # default, in min
	my $startdate = 0;
	my $enddate  = 0;
	my $docunit = 'min';
	
	#
	# Find input files = if <batch>, or <infile> = <inpath>/<batch>curve.asc,
	# then	<inpath>/<batch>curve.asc and 
	#	<inpath>/<batch>doc.asc and 
	#	<inpath>/pooled_fractions.csv 
	# else if <infile> = - and batch = 0, then STANDARD INPUT, else 
	#	<infile> only.
	#
	$infile .= ".asc" if $infile =~ /curve$/;
	my $curvefile	= $infile;
	my $docfile	= 0;
	my $inpath	= "$infile";
	my $poolfile	= "pooled_fractions.csv";
	my $warnings	= 0;
	my $flowlog	= 0;
	if ($infile =~ /curve\.asc$/) {
		# curve -> doc
		($docfile = $curvefile) =~ s/curve\.asc$/doc.asc/g;
		my $slash	= CORE::rindex($infile,"/") + 1;
		# prefix pool file with everything up to last slash in infile
		$slash = CORE::rindex($infile,"\\") + 1 if ! $slash;
		$poolfile = substr($infile, 0, $slash).$poolfile if $slash;
	} elsif ($batch) {
		$curvefile = $batch."curve.asc";
		$docfile = $batch."doc.asc";
		print  "Before: curve $curvefile, doc $docfile\n" if $debug;
		if ($infile) {
			$infile = substr($infile,0,-1) if $infile =~ /[\/\\]$/;
			$curvefile = "$infile/$curvefile";
			$docfile   = "$infile/$docfile";
			$poolfile  = "$infile/$poolfile";
			print  "After:  curve $curvefile, doc $docfile\n" if $debug;
			print  "		in $infile, pool $poolfile\n" if $debug;
		}
	} else {
		print  "WARNING: no doc file name specified:\n" if $debug;
#			"	   - peak purity calculation may not be correct.\n";
	}
	
	# Read doc file into array to be parsed 
	my @doclines = ();
	if ($docfile && $docfile ne '-') 
	{
		print  "Looking for doc file $docfile\n" if $verbose;
		if (-e $docfile) {
			print  "Opening doc file $docfile\n" if $verbose;
			if (open INDOC, $docfile) {
				print  "Reading doc file $docfile\n" if $verbose;
				while (<INDOC>) {
					chomp;
					print  "Doc file line $_\n" if $verbose > 1;
					push @doclines, $_;
				}
				print  "Got doc lines 0 to $#doclines\n" if $debug > 1;
				close INDOC;
			} else {
				print STDERR "WARNING: Cannot open doc file ".
				  "$docfile!\n$!";
			}
		} else {
			print STDERR "WARNING: Cannot find doc file $docfile!\n";
		}
	}
	
	# Read curve into array to be parsed after delimiter is detected
	my @lines = ();
	if ($curvefile ne '-') 
	{
		if ($curvefile =~ /\.asc$/) {
			$delim = "\t";
			$ncols = 2;
			$headline = 2;
			$fixedcol = 0;
			$datacol  = 1;
			print "Assuming data in line 3 columns 1&2\n" if $verbose;
		}
		open INFILE, $curvefile or die "ERROR OPENING INFILE $infile\n$!";
		while (<INFILE>) {
			chomp;
			push @lines, $_;
		}
		close INFILE;
	} else {
		# Accept concatenated doc and curve files
		my $docflag = (@doclines);
		while (<>) {
			chomp;
			$docflag = 2 if $_ =~ /Log/;
			$docflag = 0 if $_ eq "10";	# start of curve file?
			if ($docflag) {
				push @doclines, $_;
			} else {
				push @lines, $_;
			}
			# 2 blank lines at end of doc file
			$docflag -- if ($docfile && !$_);
		}
	}
	if (@doclines) {
#		print  "Call parse doc: verbose $verbose\n" if $debug;
		my $docvalues = parse_doc( {	
						doclines => \@doclines,
						debug	 => $debug,
						verbose	 => $verbose
					  } );
		$flowrate  = $$docvalues{'flowrate'} if exists $$docvalues{'flowrate'};;
		$frxnvol   = $$docvalues{'frxnvol'} if exists $$docvalues{'frxnvol'};
		$starttime = $$docvalues{'starttime'};
		$endtime   = $$docvalues{'endtime'};
		$startdate = $$docvalues{'startdate'};
		$enddate   = $$docvalues{'enddate'};
		$sourcefile= $$docvalues{'sourcefile'};
		$docunit   = $$docvalues{'unit'} if exists $$docvalues{'unit'};
		print "flow rate $flowrate, frxn vol $frxnvol, time $starttime".
			" to $endtime\n" if $debug;
		$warnings = $$docvalues{'warnings'} if exists $$docvalues{'warnings'};
		$flowlog  = $$docvalues{'flowlog'} if exists $$docvalues{'flowlog'};
		print "Flow change warnings:\n",join("\n",@$warnings),"\n" if
			$debug && $warnings;
		# Report flow change warnings if curve X axis or docunit is min:
		if (($docunit eq 'min' || $header1 eq 'min') && $warnings) {
			my $nwarn = @$warnings;
			$nwarn = 3 if $nwarn > 3;
			for (my $i=0; $i<$nwarn; $i++) {
				print STDERR "$$warnings[$i]\n";
			}
			print STDERR 
				"         Export curve file with X axix in ml\n"
				if $header1 ne 'ml';
			print STDERR "         Export documentation in ml\n" 
				if $docunit ne 'ml';
		} else {
			$warnings = 0;
		}
		
	} else {
		print STDERR "NOTE: Assuming\tflow rate\t= $flowrate mL/min,\n".
		"\t\tfraction volume\t= $frxnvol mL,\n".
		"\t\tstart time\t= $starttime min\n".
		"\t\tend time\t= $endtime min\n".
		"If flow, volume or start are wrong, peak purity calculations may also be wrong.\n";
	}

	# Read pooled_fractions.csv if needed to get from and to fraction
	# File pooled_fractions.csv has columns batch, from_frxn, to_frxn
	# If it has other columns or is in another order, it must includes
	# a header line with headers including the text "batch","from" and "to".
	
	if ($poolfile && $batch && ! ($fromfrxn && $tofrxn)) {
		if (-e $poolfile) {
			if (open POOL, $poolfile) {
				my $batch_col = 0;
				my $from_col = 1;
				my $to_col = 2;
				while (<POOL>) {
					my @vals = split /,/;
					# Find batch, from and to columns
					if ($_ =~ /batch/i && !$batch_col) {
						for my $i (0..$#vals) {
							if ($vals[$i] =~ /batch/i) {
								$batch_col = $i;
							} elsif ($vals[$i]=~/from/i) {
							$from_col = $i;
							} elsif ($vals[$i]=~/to/i) {
							$to_col = $i;
							}
						}
					} elsif($vals[$batch_col]=~/\s*$batch\s*/){
						$fromfrxn = $vals[$from_col];
						$fromfrxn =~ s/^ *//g;
						$fromfrxn =~ s/ *$//g;
						$tofrxn	= $vals[$to_col];
						$tofrxn =~ s/^ *//g;
						$tofrxn =~ s/ *$//g;
						last;
					}
				}
				close POOL;
			} else {
				print STDERR 
					"WARNING: Cannot open file $poolfile!\n$!\n";
			}
		} else {
			print STDERR "WARNING: Cannot find file $poolfile!\n";
		}
	}
	if (! ($fromfrxn && $tofrxn)) {
		print STDERR "No pooled fractions given; peak purity will not".
			" be estimated.\n";
	}		
	my $poolstart= $fromfrxn;
	my $poolend  = $tofrxn;
	if ($poolstart =~ /[A-H][0-9]/) {
		# Convert from A01-H12 to 1-48, 1-48 zig-zagging up A, down B etc.
		my $c = ord($fromfrxn) - ord 'A';	# A = 0, B = 1 ... H = 7
		my $r = substr $fromfrxn, 1;		# Cnn -> nn
		$r = 13 - $r if ($c % 2);		# swap 1 to 12 -> 12 to 1 for odd c1
		$poolstart = (($r + 12*$c - 1)%48) + 1;	# A=1-12, B=13-24, C=25-36, D=37-48
							# E=1-12, F=13-24, G=25-36, H=37-48 
		$c = ord($tofrxn) - ord 'A';		# Same for pool end
		$r = substr $tofrxn, 1;
		$r = 13 - $r if ($c % 2);
		$poolend = (($r + 12 * $c - 1) % 48) + 1;
	}

	print "Before: delimiter |$delim| ncols |$ncols| last line $#lines\n" if $debug;
	if (! $delim || !$ncols) {
		# Get a hash reference from subroutine
		my $delim_values =find_delimiter( {	lines	=> \@lines,
							debug	=> $debug,
							verbose	=>$verbose } );
		# print ":",join(":",@delim_values), "\n" if $debug;
		$delim = $$delim_values{'delim'};
		$ncols = $$delim_values{'ncols'};
	}
	print "After:  delimiter |$delim| ncols |$ncols|\n" if $debug;
	print "Before find_headline: headline $headline, fixed $fixedcol\n" if $debug;
	# if (! $headline || $fixedcol < 0) {
		my $head_values = find_headline( {	lines => \@lines, 
							delim => $delim,
							head => "min|ml",
							ncols => $ncols,
							debug	=> $debug,
							verbose	=>$verbose } );
		# If find_headline fails, bail out
		die "ERROR: no header line found!\n" if (! $head_values);
		$headline = $$head_values{'header_line'};
		$fixedcol = $$head_values{'fixed_col'};
		$header1  = $$head_values{'header_value'};
	# }
	$datacol = $fixedcol + 1 if ($ncols > 1 && $datacol < 0 || $datacol == $fixedcol);
	print "Header is line ",1+$headline,", $header1 is column ",$fixedcol+1,
		"\n" if $debug;
	# Is outfile a directory? I.e. does it end in / or \?
	my $last_char = $outfile ? substr $outfile, -1, 1 : '';
	my $outdir = $last_char eq '/' || $last_char eq '\\';
	if ($outfile && $outfile ne '-' && ! $outdir) {
		$outfile = $outfile . ".sec" if ($outfile !~ /.sec/);
		$outfile = $outfile . ".csv" if ($outfile !~ /.csv$/);
	} elsif ($curvefile && $outfile ne '-') {
		if ($outdir) {
			my $last_slash = 1 + rindex $curvefile, '/';
			$last_slash = rindex $curvefile, '\\' if ! $last_slash;
			$outfile .= substr $curvefile, $last_slash;
		} else {
			$outfile = $curvefile;
		}
		$outfile =~ s/curve//g;
		$outfile =~ s/\.asc$//g;
		$outfile =~ s/.csv$//g;
		$outfile = $outfile.".sec.csv";
	}

	my $out = 0;
	if ($outfile && $outfile ne '-') {
		if (-e $outfile) {
			my $outfile_orig = $outfile;
			my $outbase = $outfile;
			my $outext  = '';
			# Output file names should end in .sec.csv.
			# Replace xxx.sec.csv    with xxx.v2.sec.csv
			# Replace xxx.vM.sec.csv with xxx.vN.sec.csv
			# Replace xxx.yyy.zzz    with xxx.v2.yyy.zzz
			# Replace xxx.vM.yyy.zzz with xxx.vN.yyy.zzz
			# Replace xxx    with xxx.v2
			# Replace xxx.VM with xxx.vN
			my $dot	= CORE::index $outfile, '.sec.csv';
			$dot	= CORE::index $outfile, '.v'	if $dot < 0;
			$dot	= CORE::index $outfile, '.' 	if $dot < 0;
			if ($dot > 0) {
				$outbase = substr $outfile, 0, $dot;
				$outext  = substr $outfile, $dot;
				$outext =~ s/^v[0-9]*//g;
				$outext =~ s/^.//g;
				print "-- outbase $outbase outext $outext\n" if $debug;
			}
			my $out_version = 1;	# start numbering with version v2
			while (-e "$outfile") {
				$out_version ++;
				$outfile = "$outbase.v$out_version.$outext";
				print "Trying output file version $outfile\n" if $verbose;
			}
			print STDERR "WARNING: Output file  $outfile_orig already exists!\n".
				     "         Writing to   $outfile instead.\n";
		}
		open OUTFILE, ">$outfile" or die 
			"ERROR: Cannot open output file $outfile\n$!";
		$out = \*OUTFILE;
	} else {
		$out = \*STDOUT;
	} 

	# Print wrapper tag with all available metadata attributes
	# my $timestamp = localtime;	# local time in human-readable format
	# DateTime->now->set_time_zone( 'America/Los_Angeles' );
	my $timestamp = DateTime->now->ymd('-')." ".DateTime->now->time(':');
	print $out "<SEC_DATA";
	print $out "\n\tBATCH=\"$batch\""		 if $batch;
	print $out "\n\tCALCULATED=\"$timestamp\"";
	print $out   "\tVERSION=\"$version_ID\""	 if $version_ID;
	print $out "\n\tNPEAKS=\"$npeaks\"";
	print $out " WIDTH=\"$width\"";
	print $out " PEAKMIN=\"$peakmin\""	 	if $peakmin;
	print $out "\n\tSOURCEFILE=\"$sourcefile\""	 if $sourcefile;
	print $out   "\tCURVEFILE=\"$curvefile\""	if $curvefile;
	print $out "\n\tFLOWRATE=\"$flowrate ml/min\""	 if $flowrate;
	print $out   "\tFRACTION_VOLUME=\"$frxnvol ml\"" if $frxnvol;
	print $out "\n\tSTARTFRXNS=\"$starttime $docunit\""	 if $starttime;
	print $out   "\tENDFRXNS=\"$endtime $docunit\""	 if $endtime;
	print $out "\n\tPOOL_FROM_FRACTION=\"$fromfrxn\"" if $fromfrxn;
	print $out   " POOL_TO_FRACTION=\"$tofrxn\""	 if $tofrxn;
	print $out "\n\tPOOL_FROM_FRACTION_NUMBER=\"$poolstart\"" if $poolstart;
	print $out   " POOL_TO_FRACTION_NUMBER=\"$poolend\""	 if $poolend;
	print $out "\n\tSTARTDATE=\"$startdate\""	 if $startdate;
	print $out   "\tENDDATE=\"$enddate\""		 if $enddate;
#	print $out ">\n";
	
	my $column_values = get_column( {  lines	=> \@lines,
					   headline	=> $headline,
					   fixedcol	=> $fixedcol,
					   datacol	=> $datacol,
					   delim	=> $delim,
					   debug	=> $debug,
					   verbose	=> $verbose,
					} );
	# Skip if no data was found for this column
	if (! $column_values) {
		print STDERR "WARNING: No data for column $datacol!\n";
		print $out "</SEC_DATA>";
		close $out;		
		exit(1);
	}
	my $xin		= $$column_values{'xin'};
	my $yin		= $$column_values{'abs'};
	my $header2	= $$column_values{'header'};
	if ($debug > 1) {
		print "Column ",1+$fixedcol," $header1: ",join(' ',@$xin[0..5]),"\n";
		print "Column ",1+$datacol," $header2: ",join(' ',@$yin[0..5]),"\n";
	}
	my $i;
	my $j;
	my @min = ();
	my @frxn = ();
	my @abs  = ();
	my $max_abs = -1000;
	my $ifromfrxn	= -1;
	my $itofrxn	= -1;
	# Convert to time domain first, assuming a constant flow rate and 
	# assuming both the curve and doc file used the same starting point
	# (both start of run or both injection time, but not one of each).
	# The first assumption just makes computing easier; the 2nd is crucial
	# to analyzing the right segment of the data.
	# 
	# After converting to time, test for the right time range and 
	# convert to fractions.
	#
	# Convert to time by applying linear constants timemult and timeadd
	# derived from the doc file to xin, start and end times as needed:
	#   For xin (column 1) in min, docunit in min: t=xin = xin*1 + 0,
	#		no change to start and end times
	#	For xin in min, docunit in ml, t = xin*1 + 0,
	#		divide start and end (volumes) by flow to get times
	#   For xin in ml, docunit in ml: t = xin/flow (if flow !=0)
	#				(or = xin/flow + start/flow ??),
	#		divide start and end volumes by flow
	#	For xin in ml, docunit in min: t = xin/flow (+ start ??),
	#		no change to start and end times.
	#
	# In other words, timemult = 1 for xin in min, 1/flow for ml; timeadd=0;
	# start and end = start and end / flow for docunit = ml.
	#
	# In all cases, frxn = ((xin * timemult + timeadd) - starttime) / size
	# 	for all xin * timeult + timeadd > starttime and < endtime

	# TO BE DONE: recover flow log and change flow rate based on X axis
	
#	my $iflow = 0;
#	my $next_change = 0;
#	my $nflow = @$flowlog;
#	if ($nflow $$ $header1 =~ /min/i) {
#		my $change = $$flowlog[0];
#		$flowrate = $$change{'flow'};
#		$next_change = $change{'time'};
#	}
	my $flowfactor = ($flowrate > 0.0) ? $flowrate : 1.0;
	my $timemult = 1.0;
	my $timeadd  = 0;
	if ($header1 =~ /ml/i) {
		$timemult = 1/ $flowfactor;
		# $timeadd  = $starttime;
	}
	if ($docunit =~ /ml/i) {
		$starttime /= $flowfactor;
		$endtime /= $flowfactor;
	}
	my $time2frxn = $flowfactor / $frxnvol ;
	my @pooled = ();
	print  "Flow $flowrate, Size $frxnvol, time=x*$timemult + $timeadd\n" if $debug;
	print  "Pool = $poolstart to $poolend\n" if $debug;

	my $firstin = -1;
	for ($i=0; $i<@$xin; $i++ ) {
		my $x = $$xin[$i]; 
		my $time = $x * $timemult + $timeadd;
		# Start at start time, stop at end time to simplify indeces
		next if ($time < $starttime);
		$firstin = $i if $firstin < 0;
		last if ($time > $endtime);
		my $f = 1 + ($time - $starttime) * $time2frxn;
		my $this_abs = $$yin[$i];
		push @frxn, $f;
		push @min, $time;
		push @abs, $this_abs;
		# Find first and last rows in pooled fractions
		push @pooled, ($f >= $poolstart && $f < $poolend + 1.0);
		# Find max. absorbance in pooled fractions
		$max_abs = $this_abs if $this_abs > $max_abs && 
				$f >= $poolstart && $f < $poolend + 1.0;
		my $j=$#frxn;
		$ifromfrxn = $j if ($f >= $poolstart && $ifromfrxn < 0);
		$itofrxn   = $j if ($f < $poolend + 1);
		print "pt $i:$j= $x, $min[$j], $f, $abs[$j]\n" if $debug>1;
	}
	print " frxn index = $ifromfrxn to $itofrxn\n" if $debug;

	# Only use points from start time to end time
	my $npoints = @abs;
	print "Found $npoints points, frxns $frxn[0] to $frxn[$#frxn]\n" if $debug;	
	print $out "\n\tOBSERVED_POINTS=\"$npoints\"";
	print $out ">\n";

	print $out "\t<WARNINGS>\n",join("\n",@$warnings),"\n</WARNINGS>" 
		if ($warnings);

	# Pass data to sec_fit for recursion, including pool but not guess -
	# let sec_fit do its own guessing.

	my $maxpeaks = 0;
	my $model = 0;
	my $variance = 1;
	my $temp_dir = 0;
	if (1) {	# always get derived value
		# Accept peakmin as percent; ALWAYS convert to a fraction
		$peakmin = $peakmin/100;	# if $peakmin && $peakmin >= 1;
		print "Calling sec_fit, peak min $peakmin\n" if $debug;
		my $all_models = sec_fit( {
						fraction	=> \@frxn,
						absorbance	=> \@abs,
						name		=> $batch,
						peak_count	=> $npeaks,
						peakmin		=> $peakmin,
						pooled		=> \@pooled,
						verbose		=> $verbose,
						debug		=> $debug,
						gnuin		=> $gnuin,
						gnuout		=> $gnuout,
						slopeflag	=> $slopeflag,
						width		=> $width,
						plot		=> $plot,
					} );
		print Dumper( $all_models ) if $debug > 1;
		
		$maxpeaks = $$all_models{'max_peaks'};
		$model = $$all_models{'model'};
		$variance = $$all_models{'variance'};
		$temp_dir = $$all_models{'directory'};
		
		# Returned hash is: my %all_models = 
		# { max_peaks=> N,   # this is the max. iteration
		#	model	 => { 
		#			1	=> {  
		#				solution	=> {
		#					slope		=> 1
		#					intercept	=> 1
		#					a1			=> 1
		#					b1			=> 1
		#					c1			=> 1
		#				},
		#				y_fitted	=> [],
		#				y_residual	=> [],
		#				R_abs		=> 1,
		#				RMSD		=> 1,
		#				mean		=> 1,
		#			},
		#			.....
		#			N	=> {  
		#				solution	=> {
		#					slope		=> 1
		#					intercept	=> 1
		#					a1			=> 1
		#					b1			=> 1
		#					c1			=> 1
		#					....
		#					an			=> 1
		#					bn			=> 1
		#					cn			=> 1
		#				},
		#				y_fitted	=> [],
		#				y_residual	=> [],
		#				R_abs		=> 1,
		#				RMSD		=> 1,
		#				mean		=> 1,
		#			},
		#		}
		# }
	}
					
	# Print derived values, if any
	my @write_list = ();

	if ($maxpeaks <= 0) {
		print STDERR "WARNING: NO PEAKS FOUND in $curvefile\n";
	} else {
		print  "Sec_fit found $maxpeaks peaks\n" if $debug;
		my @all_rabs	= (1);		# store Rabs for i peaks in rabs[i]
		my @all_rms	= (0);
		my @all_pp	= (0);
		my @best_pp	= ();
		my @all_bic	= (0);
		my @all_area	= (0);
		my @all_resid	= (0);
		my @all_pp_resid= (0);

		# Calculate summary values for each peak first, and
		# find best model = peak with lowest BIC = n*ln(R) + (2+3g)*ln(n)
		# where n = number of points, g=number of gaussians, assuming
		# each model has 2 parameters for bg, 3 for each gaussian.
		#
		# Then write summary values for best peak, then for all peaks

		my $bicmin	= 1e10;
		my $best_npk	= 0;
		my %best_list	= ();
		my $frxn_factor = (@frxn ? ($frxn[-1]-$frxn[0]) / @frxn : 1);
		my $rabs_old	= 1;
		my $d_rabs_stop = 0;
		foreach my $g (1..$maxpeaks) {
			my $npeak    = $$model{"$g"};
			my $solution = $$npeak{'solution'};
			my $slope	= $$solution{'slope'};
			my $interc	= $$solution{'intercept'};
			my $rabs  = $$npeak{'R_abs'};

			# Old logic: 
			# Stop after first peak with delta Rabs < peakmin,
			# 	except in verbose mode
			
			# New logic: find last peak with delta Rabs < peakmin

			# $d_rabs_stop = $d_rabs_stop || ($g > 1 && $rabs_old-$rabs < $peakmin);
			# last if ($d_rabs_stop && ! $verbose);
			$d_rabs_stop = $rabs_old-$rabs < $peakmin;

			my $rmsd  = $$npeak{'RMSD'};
##			my $resid = $$npeak{'y_residual'};
##			my $backg = $$npeak{'background'};
##			my $gauss = $$npeak{'gaussian'};
			print  "Model $g: m $slope b $interc \n" if $debug;
			
##			my $ssq  = 0;		# sum of squared residuals
##			my $rabs = 0;		# sum of abs. value of resid.
			my @gpp  = (0); 	# pool abs. in each gaussian
			my $pp   = 0;		# best pooled purity
			my $tot = 0;		# total area
			my $ptot = 0;		# area of pool
			my $i;
			
			# Add up residuals, total obs - bg, pool obs - bg
			
			foreach $i (0..$#abs) {
				print "Point $i: $frxn[$i], $abs[$i]\n" if $debug>1;
##				my $r = $$resid[$i];
##				$ssq += $r * $r; 
##				$rabs += abs($r);
				my $f = $frxn[$i];
				# Calculate obs - bg, bg = m*x + b
				my $y = $abs[$i] - ($slope*$f+$interc);
				$tot += $y;
				$ptot += $y if $pooled[$i];
			}
			
			$tot *= $frxn_factor;
			my $tot_left = $tot;
			my $pp_left = $ptot;
			# Add up each Gaussian's contribution to pooled frxns
			
			my $i_gpp_max = 1;
			foreach my $p (1..$g) {
##				my $gp = $$solution{"$p"};
				my $a = abs($$solution{"a$p"});
				my $b =     $$solution{"b$p"};
				my $c = abs($$solution{"c$p"});
				$tot_left -= $a * $c * $sq2p;
				push @gpp, 0;	# add gpp[$p]
				# $gpp[$p] = 0;
				foreach $i ($ifromfrxn .. $itofrxn) {
					my $f = $frxn[$i] - $b;
					# y(p) = a * e^(-(x-b)^2/(2c^2))
					$gpp[-1] += $a * exp(-$f*$f/(2*$c*$c));
					print "     Peak $p: pool part $i, ".
						"frxn $frxn[$i], dx $f =",
						$a * exp(-$f*$f/(2*$c*$c)),"\n" 
						if $debug > 1
				}
				$pp_left -= $gpp[-1];
				$gpp[-1] /= ($ptot > 0 ? $ptot : 1);
				$i_gpp_max = $p if $gpp[-1] > $gpp[$i_gpp_max];
				print "  Peak $p: pool=$gpp[-1]\n" if $debug;
			}

			# Deal with summary statistics
			
#			my $rms = sqrt($ssq/$npoints);	# root mean sqr resid.
#			my $bic = $npoints * log($rmsd)+(2+3*$g)*log($npoints);
			# Use error variance = RMSD ^ 2; log(R^2) = 2*logR
			my $bic = $npoints*2*log($rmsd)+(2+3*$g)*log($npoints);
			# Find best Bayesian Information Criterion (bic) 
			# among peaks with > peakmin change in Rabs
			if ($g == 1 || ( $bic < $bicmin && ! $d_rabs_stop) ) {
				$bicmin = $bic;
				$best_npk = $g;
				@best_pp = @gpp;
			}
			push @all_rms,  $rmsd;
			push @all_rabs, $rabs;	#  / ($tot ? $tot : 1);
			push @all_pp,   $gpp[$i_gpp_max];
			push @all_bic,	$bic;
			# Area by fraction = total for each point * range/npoints
			push @all_area,	$tot;
			push @all_resid, $tot_left / ($tot ? $tot : 1);
			push @all_pp_resid,	$pp_left / ($ptot ? $ptot : 1); 
			$rabs_old = $rabs;
		}
		#########################################################################
		# <DERIVED_VALUES> has attributes:
		#	R1	= Rabs for 1 Gaussian
		#	PP1	= Pool purity for 1 Gaussian
		#	DIRECTORY = gnuplot intermediate file folder
		# and contains a table with NPEAKS,RABS,RMSD,PP,BIC for 1 peak
		# and for the best model, or in verbose mode all models,
		# plus elements:
		#	<BEST_MODEL> and, for verbose mode, 
		#	<MODEL> for each model tried.
		# Both <BEST_MODEL> and <MODEL> elements have attribute:
		#	NPEAKS  = number of Gaussians in model
		#	RABS	= sum of |obs-calc| / sum of (obs - bg)
		#	RMSD	= sqrt(sum of (obs-calc)^2 / N points)
		#	PP	= Purity of pooled fractions
		#	AREA	= total of observed - background
		# <BEST_MODEL> also contains the attribute:
		#	PEAKMIN = minimum change in Rabs as a fraction of total area
		#	and a csv table with headers:
		#	Peak Number = order by peak start position,
		#	Start, End = center -/+ 2*sigma, covering 95% of peak area
		#	Fraction of Total (observed area), and 
		#	Fraction of Pooled (area).
		# with one row for each peak plus totals for areas (should be 1.000)
		############################################################################
		print  "Start printing derived values\n" if $debug;
		print $out "\n\t<DERIVED_VALUES";
		my $ws = " ";
		printf $out "%sR1=\"%.4f\"", $ws,$all_rabs[1];
		printf $out "%sPP1=\"%.4f\"", $ws,$all_pp[1];
		printf $out "%sSMOOTHING_WIDTH=\"%.2f\"", $ws,$width;
		printf $out "%sDIRECTORY=\"%s\"", $ws,$temp_dir if $temp_dir;
		print $out ">\n";
		$best_npk=1 if ! $best_npk;

		# Write CSV with best or all model stats
		@write_list = $verbose ? (1..$maxpeaks) : ($best_npk);
		print $out "\nNPEAKS,RABS,dRABS,RMSD,AREA,PP,BIC\n";
		foreach my $m (@write_list) {		
			printf $out "%i,%.4f,%.4f,%.4f,%.1f,%.3f,%.3f\n", 
				$m, $all_rabs[$m], $all_rabs[$m-1]-$all_rabs[$m],
				$all_rms[$m], $all_area[$m], $all_pp[$m], $all_bic[$m];
		}
		
		# Print values for min BIC
		print $out "\n\t\t<BEST_MODEL";
		printf $out "%sPEAKMIN=\"%.4f\"", $ws,$peakmin;
		print $out $ws,"NPEAKS=\"$best_npk\"";
		printf $out "%sRABS=\"%.4f\"", $ws,$all_rabs[$best_npk];
		printf $out "%sRMSD=\"%.4f\"", $ws,$all_rms[$best_npk];
		printf $out "%sPP=\"%.4f\"",   $ws,$all_pp[$best_npk];
		printf $out "%sBIC=\"%.3f\"",  $ws,$all_bic[$best_npk];
		printf $out "%sAREA=\"%.3f\"", $ws,$all_area[$best_npk];
		print $out ">\n";
		my $k = 0;
				
		# Write best model peak table
		
		my $best_npeak = $$model{$best_npk};
		my $best_sol   = $$best_npeak{'solution'};
		print "best_npk $best_npk\n" if $debug;
		
		# Write csv peak table in ndx order
		print $out "Peak Number,Start,End,Fraction of Total Area,Fraction of Pool Area\n";
		foreach my $p (1..$best_npk) {
			my $a = abs($$best_sol{"a$p"});
			my $b =     $$best_sol{"b$p"};
			my $c = abs($$best_sol{"c$p"});
			printf $out "%i,%.2f,%.2f,%.3f,%.3f\n", 
				$p, 		# peak num
				$b - 2 * $c,	# start
				$b + 2 * $c,	# end
				$a * $c * $sq2p / $all_area[$best_npk],	# frxn of total area
				$best_pp[$p];	# frxn of pool
		}
		printf $out "Total,,,%.3f,%.3f\n",
			1 - $all_resid[$best_npk],	# frxn of total
			1 - $all_pp_resid[$best_npk];	# frxn of pool
		print $out "\t\t</BEST_MODEL>\n";
		
		# Models

#		<MODEL> contains <BG> and <PEAK>s and has attributes 
#			NPEAKS, RABS (total absolute residual), PP (peak purity).
#			<BACKGROUND> has attrubutes INTERCEPT, SLOPE, SQUARE 
#				where background=INTERCEPT+x*SLOPE+x^2*SQUARE
#			<PEAK> has attributes PEAKNUM, CENTER, HEIGHT, SIGMA

		print $out "<!--\n".
			"#\tBackground = INTERCEPT + SLOPE*x\n".	# + SQUARE*x^2\n".
			"#\tEach Gaussian peak is ".
			" HEIGHT * exp(-(x-CENTER)^2/(2*SIGMA^2))\n".
			"#\t'Start' and 'End' of peaks are fractions at 2*SIGMA from CENTER,\n".
			"#\tenclosing 95% of the peak's area.\n".
			"-->\n";
			#	"#\twhere SIGMA*2*sqrt(2ln(2)) = FWHM = ".
			#	" full width at half max.\n#\n";

		foreach my $g (@write_list) {
			my $npeak = $$model{$g};
			my $solution = $$npeak{'solution'};
			my $slope  = $$solution{'slope'};
			my $interc = $$solution{'intercept'};
			my $sqr    = 0;
			$sqr	  = $$solution{'square'} if exists $$solution{'square'} ;
#			my $gauss = $$npeak{'gaussian'};

			# print <Model npeaks="n" rabs="r" pp="p">

			print  $out "\t\t<MODEL NPEAKS=\"$g\"";
			printf $out " RABS=\"%.4f\"",	$all_rabs[$g];
			printf $out " RMSD=\"%.4f\"",	$all_rms[$g];
			printf $out " PP=\"%.4f\"",	$all_pp[$g];
			printf $out " AREA=\"%.3f\"",	$all_area[$g];
			printf $out " BIC=\"%.3f\">\n", $all_bic[$g];

			# print <background intercept= slope= square=/>

			print  $out "\t\t\t<BACKGROUND ";
			printf $out " INTERCEPT=\"%.3f\"", $interc;
			printf $out " SLOPE=\"%.3f\"",	   $slope if $slopeflag;
			printf $out " SQUARE=\"%.3f$\"",   $$sqr if 
							exists $$solution{'square'};
			print $out "/>\n";

			foreach (1..$g) {
				# print <Peak peaknum ctr ht fwhm/>
#				my $pk = $$gauss{"$_"};
				my $a = abs($$solution{"a$_"});
				my $b =     $$solution{"b$_"};
				my $c = abs($$solution{"c$_"});
				print  $out "\t\t\t<PEAK PEAKNUM=\"$_\"";
				printf $out " HEIGHT=\"%.3f\"",	$a;
				printf $out " CENTER=\"%.3f\"",	$b;
				printf $out " SIGMA=\"%.3f\"",	$c;	# *$c2f;
				printf $out " AREA=\"%.3f\"/>\n",$a*$c*$sq2p;	
			}
			print $out "\t\t</MODEL>\n";
		}
		print $out "\t</DERIVED_VALUES>\n";

	}

	if ($print_raw) {
		# Print raw data if requested
		# <RAW_DATA npoints="n"> = csv with min, mAu (observed)
		# Print wrapper and headers
		if (! $header1) {
			print STDERR "WARNING: no header for first column!\n".
					 "		 - using default of 'min'\n";
			$header1 = 'min';
		}
		if (! $header2) {
			print STDERR "WARNING: no header for second column!\n".
					 "		 - using default of 'mAu'\n";
			$header2 = 'mAu';
		}
		print $out	"\n\t<RAW_DATA NPOINTS=\"$npoints\">\n".
				"$header1,$header2\n";
		foreach $i (0..($npoints-1)) {
			printf $out "%.3f,%.3f\n", $$xin[$firstin + $i], $abs[$i];
		}
		print $out "\t</RAW_DATA>\n";
	}	

	if ($print_model && $maxpeaks > 0) {

		# Print model values if requested
		# <MODEL_VALUES> = csv with headers: 
		# fraction, absorbance, pool, Gn, Gn_res, Gn_bg, Gn_1, ... Gn_n
		# for each model (repeating the 1st 3 for plotting convenience)
		# where fraction = (time - start time)*flow rate/fraction size
		# For each model Gn = G1, G2, ... G<max>,
		#	Gn	= total for n-Gaussian model (bg + n peaks), 
		#	Gn_res	= residual = obs. abs - Gn total, 
		#	Gn_bg	= background values
		#	Gn_1, Gn_2, ... Gn_n = Gaussian curves 1 to n
		
		# Wrapper tag
		print $out "\n\t<MODEL_VALUES MAXPEAKS=\"$write_list[-1]\">\n";

		# Old: do not repeat 1st 3 columns in verbose mode:		
		# Headers = fraction,absorbance,pool,
		#			G1, G1_res, G1_bg, G1_1,
		#			G2, G2_res, G2_bg, G2_1, G2_2,
		#			G3, G3_res, G3_bg, G3_1, G3_2, G3_3
		# This means data for a model with n Gaussians starts in column
		# 3 + 3(n-1) + n*(n-1)/2 = (n^2 + 5n)/2 (numbering from 1 up)
		#		
		# New: repeat absorbance and pool in verbose mode: 
		# Headers =fraction, 
		#	   absorbance, pool, G1, G1_res, G1_bg, G1_1,
		# 	   absorbance, pool, G2, G2_res, G2_bg, G2_1, G2_2,
		# 	   absorbance, pool, G3, G3_res, G3_bg, G3_1, G3_2, G3_3
		# This means data for a model with n Gaussians starts in column
		# 1 + 5(n-1) + n*(n-1)/2 = (n^2 + 9n)/2 - 4 and goes to
		#     5n     + n*(n+1)/2 = (n^2 + 11n)/2 (numbering from 1 up)

		# Old 1st 3 columns
		# print $out "fraction,absorbance,pool";
		print $out "fraction";
		foreach $i (@write_list) {
			# New 1st 2 columns
			print $out ",absorbance,pool";
			# Add ,gn,gn_res,gn_bg
			print $out ",G$i,G$i","_res,G$i","_bg";
			# Add ,Gi_j
			for ($j=1; $j<=$i; $j++) {
				print $out ",G$i","_$j";
			}
		}
		print $out "\n";
		
		# Print data, formatting as you go for simplicity
		foreach my $f (0..($npoints - 1)) {
			my $frx = $frxn[$f];
			print  "Fraction $frx at $f\n" if $debug>2 && $f<5;

			# OLD Print fraction, observerd absorbance
			# printf $out "%.3f,%.3f,%i", $frx,$abs[$f], $pooled[$f]*$max_abs;
			printf $out "%.3f", $frx;
			foreach my $g (@write_list) {
				# NEW Print fraction, observerd absorbance
				printf $out ",%.3f,%i", $abs[$f],
							$pooled[$f]*$max_abs;
				my $npeak = $$model{$g};
				my $yfit  = $$npeak{'y_fitted'};
				my $resid = $$npeak{'y_residual'};

				print  "Npeak $g keys: ",
					join(",",keys %$npeak),"\n" if $debug>2 && $f<5;

				# Print total fit, residual for g-peak model
				printf $out ",%.3f,%.3f", $$yfit[$f], $$resid[$f];
				my $solution = $$npeak{'solution'};
				my $slope  = $$solution{'slope'};
				my $intercept  = $$solution{'intercept'};

				# Print background for g-peak model
				printf $out ",%.3f", $frx*$slope + $intercept;
#				my $gauss = $$npeak{'gaussian'};
##				print  "Gauss $g keys: ",
##					join(",",keys %$gauss),"\n" if $debug>2 && $f<5;
				foreach my $p (1..$g) {
##					my $gp = $$gauss{"$p"};
##					print  "Peak $p keys: ",
##						join(",",keys %$gp),"\n" if $debug>2 && $f<5;
					my $a = abs($$solution{"a$p"});
					my $b = $$solution{"b$p"};
					my $c = $$solution{"c$p"};
					my $fbc = ($frx - $b) / $c;

					# Print each peak for g-peak model
					printf $out ",%.3f",$a*exp(-$fbc*$fbc/2);
				}
			}
			print $out "\n";
		}
		print $out "\t</MODEL_VALUES>";
	}		
	print $out "\n</SEC_DATA>\n";
	close OUTFILE;
	print STDERR "Done\n" if $verbose;
	exit(0);
}

####################################################################
#
# find_headline - find header line of table = 1st with N delimiters 
#
# Pass in array of lines, delimiter, and number of columns
# Pass out index of first line with the right number of non-blank columns,
# and the column of the fixed text e.g. "Tem" for temperature.
####################################################################

sub find_headline {
	my %params = validate( @_, {	lines => { type => ARRAYREF },
						delim => { type=>SCALAR, default=>',' },
						head  => { type=>SCALAR,default=>'Tem'},
					ncols => { type=>SCALAR, default=>2 }, 
						debug => { type => SCALAR, default => 0 },
						verbose =>{type => SCALAR, default => 0 }
				  } );
	my @lines = @{$params{'lines'}};
	my $delim =   $params{'delim'};
	my $head  =   $params{'head'};
	my $ncols =   $params{'ncols'};
	my $debug =   $params{'debug'};
	my $verbose = $params{'verbose'};
	my $i = 0;
	my $n = $#lines;
	my $old_count = -1;
	my $fallbackline = -1;
	my $fallbackcol = -1;
	my @heads = split /\|/, $head;
	my $foundhead = 0;
	print "Find headline with header ",join(" or ",@heads),"\n" if $debug;
	foreach my $line (@lines) {
		print "Line |$line|\n" if $debug > 1;	
		next if (! $line =~ /$delim/);
		# Count number of non-blank fields split by delimiter
		my $count = 0;
		my @fields = split /$delim/, $line;
		print "Found ",$#fields+1," fields\n" if $debug;
		my $headcol = -1;
		my $j = 0;
		foreach (@fields) { 
			s/\s*//g;
			if ($_ ne '') { 
				$count ++; 
				foreach my $h (@heads) {
					print "Look for header $h\n" if $debug>1;
					if ($_ =~ /^$h/) {
						$headcol = $j;
						$foundhead = $h;
						print "FOUND $h at ",$j+1,"\n" if
								$debug;
					}
				}
			}
			$j ++;
		}
		return ( { header_line => $i, fixed_col => $headcol, 
			  header_value => $foundhead } )
			  if ($count > 1 && $headcol >= 0);
		$fallbackline = $i if ($count == $ncols);
		$fallbackcol = $headcol if ($headcol >= 0);
		# if ($count) {
			# # Stop on finding 2 in a row with almost the right number
			# return ($i-1) if ($count == $old_count);
			# $old_count = $count;
		# }
		$i ++;

		# Give up half-way through the file.
		last if ($i > $n / 2);
	}
	my @next_line_fields = split /$delim/, $lines[$fallbackline + 1];
	for (my $i=0; $i<@next_line_fields; $i++) {
		return ( { header_line => $fallbackline, fixed_col => $i } )
			if ( $next_line_fields[$i] =~ /^\s*[0-9][0-9]\.[0-9]$/);
	}
	return (0);
}

####################################################################
#
# find_delimiter = find comma, tab or other consistent delimiter 
#
# Returns delimiter character if any, number of columns if fixed, or 0,
#	   and number of columns found.
# Pass in reference to file handle and (optionally) number of digits expected:
#	digits_left_1  = Number to left of decimal place in first real number
#	digits_right_1 = Number of right of decimal place in first real number
#	digits_left_2  = Number to left of decimal place in other real numbers
#	digits_right_2 = Number of right of decimal place in other real numbers
#
# #################################################################

sub find_delimiter {
	my %params = validate( @_,{	
			lines => { type	=> ARRAYREF },
			digits_left_1	=> { type => SCALAR, default => 2 },
			digits_right_1	=> { type => SCALAR, default => 1 },
			digits_left_2	=> { type => SCALAR, default => 1 },
			digits_right_2	=> { type => SCALAR, default => 4 },
				debug		=> { type => SCALAR, default => 0 },
				verbose		=> { type => SCALAR, default => 0 }
				  } );

	my @lines = @{$params{'lines'}};
	my $nlines = @lines;
	my $dig_l1 = $params{'digits_left_1'};
	my $dig_r1 = $params{'digits_right_1'};
	my $dig_l2 = $params{'digits_left_2'};
	my $dig_r2 = $params{'digits_right_2'};
	my $debug  = $params{'debug'};
	my $verbose = $params{'verbose'};
	print "find_delimiter in $nlines, digits = $dig_l1,$dig_r1,$dig_l2,$dig_r2\n" if ($debug > 1);

	# Scan lines of input file for 10 lines with consistent delimiter 
	# Count number of each possible delimiter character on each line
	#
	# Store number of times the count is the same for each character
	# Stop when any character's "same" count reaches 10, or the number 
	# of lines read reaches 100 (no delimiter found)

	my $delim_str = ",\tsw:;|~!@#\$%^&*_=i/?()[]{}<>";
	my @delim_arr = split //, $delim_str;
	print join(',',@delim_arr),"\n" if $debug > 1;
	my @count = ();
	my $n = 0;
	my $i = 0;
	my $j = 0;
	foreach $i (@delim_arr) {
		push @count, 0;
	}
	my $m1 = 50;
	my $m2 = 150;
	if ($nlines < 150) {
		$m1 = $nlines / 3;
		$m2 = $nlines;
	}
	my $m = $m2 + 1 - $m1;
	print "m=$m m1=$m1 m2=$m2\n" if $debug > 1;
	# my $maxlines = 100;
	# my $start_test = $nlines/2;

	# Test the lines from $m1 to $m2 for delimiter

	# Initiallize hash of delim and linenum-delim
	my %occurs;
	my %totals;
	foreach $j (@delim_arr) {
		$j = 't' if ($j eq "\t");
		$totals{$j} = 0;
		foreach $i ($m1..$m2) {
			$occurs{ "$i-$j" } = 0;
		}
	}
	# $totals{'white'} = 0;
	# foreach $i ($m1..$m2) { $occurs{ "$i-white" } = 0; }

	foreach $i ($m1..$m2) {
		my $white = 1;
		print "$i: $lines[$i]\n" if $debug >2;
		my @chars = split //, $lines[$i];
		print join(':',@chars),"\n" if $debug>3;
		foreach $j (@chars) {
			print "'$j'\n" if $debug >3;
			if ($j eq ' ' || ($j ne 'w' && CORE::index($delim_str, $j) >= 0))
				{
				print CORE::index($delim_str,$j),"\n" if $debug >2;
				$j = 't' if ($j eq "\t");
				$j = 's' if ($j eq ' ');
				$occurs{ "$i-$j" } = $occurs{ "$i-$j" } + 1;
				$totals{ $j } = $totals{ $j } + 1;
				print "'$j' '$i-$j'\n" if $debug >2;
				# print "occurs $occurs{"$i-$j"}\n" if $debug >3;
			}
			if ( $j =~ /\s/ ) {
				if (! $white ) {
					$occurs{"$i-w"} =$occurs{"$i-w"}+1;
					$totals{"w"} =$totals{"w"}+1;
					$white = 1;
				}
			} else  { 
				$white = 0 
			};
		}
	}
	my @any = ();
	my $nperline = 0;
	foreach $j (@delim_arr) {
		$j = 't' if ($j eq "\t");
		$j = 's' if ($j eq ' ');
		print "'$j' total: $totals{$j}\n" if $debug > 1;
		if ($totals{$j}) {
			$nperline = $occurs{"$m1-$j"};
			print "$nperline\n" if $debug > 1;
			foreach $i (($m1+1)..$m2) {
				print $occurs{"$i-$j"},"=? $nperline\n" if $debug > 2;
				if ($occurs{"$i-$j"} < $nperline) {
					$nperline = -1;
					last;
				}
			}
			print "$j $nperline\n" if $debug > 1;
			if ($nperline > 0) {
				$j = "\t" if ($j eq 't');
				$j = " " if ($j eq 's');
				print "return '$j',",($nperline+1),"\n" if $debug > 1;
				# return ( ( "$j", $nperline+1 ) );
				return( {delim => "$j",ncols =>($nperline+1)} );
			} else {
				push @any, $j;
			}
		}
	}
	# Check for fixed length - to be done
	# Return first delimiter that occurs at least once on each line
	my $max = 0; 
	my $most = 0;
	foreach $j (@any) {
		$i = $totals{$j};
		return ( {delim => $j, ncols => int($i/$m) }) if ($i >= $m);
		if ($i > $max) {
			$max = $i;
			$most = $j;
		}
	}
	return ( { delim => "$most", ncols => (int($max/$m)+1) } );
}



###################################################################
#
# get_column - pull out a pair of columns from data lines
#
# Pass in data lines, header line number, fixed and data column numbers
# and delimiter.
# Pass out header value and arrays for data and fixed columns.
# datacol is column number after fixed col, counting from 1 = 1st to the right
#
###################################################################

sub get_column {
	my %params = validate ( @_, {	lines	=> { type => ARRAYREF },
					headline=> { type=>SCALAR, default=>0 },
					fixedcol=> { type=>SCALAR, default=>0 },
					datacol	=> { type=>SCALAR, default=>1 },
					delim	=> {type=>SCALAR,default=>"\t"},
					debug	=> { type=>SCALAR, default=>0},
					verbose	=> { type=>SCALAR, default=>0}
					} );
# start:
	my $i;
	my @fixed = ();
	my @data  = ();
	my $dcol = 0;
	my $head = 0;
	my @lines	= @{$params{'lines'}};
	my $delim	= $params{'delim'};
	my $headline	= $params{'headline'};
	my @headers	= split /$delim/, $lines[$headline];
	my $datacol	= $params{'datacol'};
	my $fixedcol	= $params{'fixedcol'};
	my $debug	= $params{'debug'};
	my $verbose	= $params{'verbose'};
	print  "Get columns after line ",1+$headline,": |$lines[$headline]|\n" if $debug;
	print  "Delim: |$delim|, columns $fixedcol, $datacol\n" if $debug;
#	print  "Headers: |",join("|,|",@headers),"|\n" if $debug;
	
	if ($datacol =~ /^[0-9]+$/) {
 		$dcol = $datacol;		# For Tm: + $fixedcol;
		return (0) if ($#headers < $dcol);
		$head = $headers[$dcol];
	} else {
		$head = $datacol;
		foreach $i (($fixedcol+1)..@headers) {
			if ($headers[$i] =~ /^\s*datacol\s*$/) {
				$dcol = $i;
				last;
			}
		}
		return (0) if (! $dcol);
	}
	$head =~ s/^\s+//g;
	$head =~ s/\s+$//g;
#	print  "Data column head: |$head|\n" if $debug;
	$delim = ' ' if $delim eq 's';
	$delim = '\s' if $delim eq 'w';
	my $ilim = $headline + 5;
	foreach $i (($headline+1)..@lines) {
#		print STDERR "Line $i: |$lines[$i]|\n" if $debug && $i < $ilim;
		next if (! $lines[$i]);
		my @fields = split /$delim/, $lines[$i];
		foreach my $c ($fixedcol, $dcol) {
			$fields[$c] =~ s/^\s+//g;
			$fields[$c] =~ s/\s+$//g;
		}
#		print STDERR "Fields: |",join("|,|",@fields),"|\n" if $debug && $i < $ilim;
		# Stop if time or volume is 0 and we're not at the start
		last if (@fields < $fixedcol || 
				($i > $ilim && ! $fields[$fixedcol]));
		# Skip if not enough fields in record
		next if (@fields < $dcol);
		push @fixed, $fields[$fixedcol];
		push @data,  $fields[$dcol];
#		print STDERR "fixed: |",join("|,|",@fixed),"|\n" if $debug && $i < $ilim;
	}
	return ( { header => $head,
		   xin	  => \@fixed,
		   abs	  => \@data } );
}
		

##############################################################################
#
# sec_fit: Pass an initial guess to LM fitting routine and return new params
#
#   In: 
#	fraction	= from time or volume 
#	absorbance	= mAu
#	npeaks		= number of peaks to fit
#
#   Out:
#	Pointer to hash of same shape as guess with fit values, plus version
#		and reason for stopping fit; or 0 if fit fails to converge.
#
##############################################################################

	

sub sec_fit {
	my %params = validate( @_, {	
					absorbance	=> { type => ARRAYREF },
					fraction	=> { type => ARRAYREF },
					peak_count	=> { type => SCALAR, default => 1},
					peakmin		=> { type => SCALAR, default => .05},
					name		=> { 
										type	=> SCALAR, 
										default	=> 'batch_unamed'
					},
					pooled		=> { type => ARRAYREF },
					verbose		=> { type => SCALAR, default => 0 },
					debug		=> { type => SCALAR, default => 0 },
					gnuin		=> { type => SCALAR, default => 0 },
					gnuout		=> { type => SCALAR, default => 0 },
					slopeflag	=> { type => SCALAR, default => 0 },
					width		=> { type => SCALAR, default => 0 },
					plot		=> { type => SCALAR, default => 0 },
				} );



	# Pull arguments from hash
	my $fraction	= $params{'fraction'};
	my $absorbance	= $params{'absorbance'};
	my $pooled	= $params{'pooled'};
	my $peak_count	= $params{'peak_count'};
	my $peakmin	= $params{'peakmin'};
	my $verbose	= $params{'verbose'};
	my $debug	= $params{'debug'};
	my $name	= $params{'name'};
	my $gnuin	= $params{'gnuin'};
	my $gnuout	= $params{'gnuout'};
	my $slopeflag	= $params{'slopeflag'};
	my $width	= $params{'width'};
	my $plot	= $params{'plot'};

	# Make (and maybe remove a preexisting) directory to store gnuplot's files
	DateTime->now->set_time_zone( 'America/Los_Angeles' );
	my $ymd = DateTime->now->ymd('_');
	my $hms = DateTime->now->time('_');

	my $working_directory = "/tmp/$name" . "_" . $ymd . "_" . $hms;
	if ($gnuin && -e $gnuin && -d $gnuin) {
		$working_directory = $gnuin;
		print "Reading from gnuplot intermediate directory $gnuin\n" 
			if $verbose;
	} else {
		$gnuin = 0;
		if ($gnuout) {
			$working_directory = $gnuout;
			print "Writing to gnuplot intermediate directory $gnuout\n" 
				if $verbose;
			if (-e $working_directory) {
				print "OOPS! gnuplot directory $gnuout exists!\n" 
					if $verbose;
				if (!remove_tree($working_directory)) {
					print "Failed to remove $gnuout.\n" 
						if $verbose;
					my $suffix = 1;
					while (-e $working_directory."_$suffix") {
						print $gnuout."_$suffix exists!\n" 
							if $verbose;
						$suffix ++;
					}
					$working_directory .= "_$suffix";
					print "Writing to gnuplot directory $working_directory\n" 
						if $verbose;
				}
			}
		}
	}
	
	if (! $gnuin) {
		my @existing_directory = make_path( $working_directory );
		
		if (scalar( @existing_directory) == 0 ) {
			my $new_directory_suffix = int( rand(10) );
			$working_directory .= $new_directory_suffix;
			if ($verbose ) {
				print "You already had an existing directory.   \t\nWe are going to try another directory  $working_directory..\n";
			}
			# we tried to remove the pre-existing directory to remake it, but was getting problems in windows.  
			#remove_tree( $working_directory );
			#make_path( $working_directory );
		}
	}
	my $file_stem = $working_directory .  "/" .  $name  ;

	if ($verbose) {
		print "writing to path at $working_directory\n"
	}
	
	# sanity test that all inputs have the same length
	my $observed_points = _sanity_check( { 
				fraction	=> $fraction,  
				absorbance	=> $absorbance, 
				pooled		=> $pooled, 
				verbose		=> $verbose,
				debug		=> $debug,
		}  
	);

	# Create two PDL arrays and fill with the x and y values
	my $x = PDL->new($fraction);		
	my $y = PDL->new($absorbance);		
	my $p = PDL->new($pooled);


	my $sum_y = sum( $y );

	if ( $debug ) {
		print "\t******* Beginning the SEC fit process ******\n";
		print "\tSum of y is $sum_y\n";
		print "\t$observed_points data points are observed\n";
	}


	# setup the input data files with the observed values into name.dat
	# and the initial parameters file
	my $gp_data_file = _generate_gnuplot_datafilename( { stem => $file_stem  } );
	wcols $x, $y, $gp_data_file;

	my $iteration = 0;
	my $rabs_old = 1;

	# my $a_init_1 = 0;
	# my $a_init_n = 0;
	
	my %all_models	= ( 	directory	=> $working_directory,
#				model		=> (),
			  );

	while ( $iteration < $peak_count ) {

		$iteration ++;
		my $gp_params_file = _generate_gnuplot_datafilename( { stem => $file_stem, version =>  $iteration, type => "parameter" } );

		if ($verbose) {
		}

		my $init_params =
		_setup_parameters(	{
						name		=> $name,
						iteration	=> $iteration,
						x		=> $x,
						y		=> $y,
						verbose		=> $verbose,
						debug		=> $debug,
						sum_y		=> $sum_y,
						file_stem	=> $file_stem,
						gnuin		=> $gnuin,
						slopeflag	=> $slopeflag,
						width		=> $width,
					}
		);
		print "Initial params for iteration $iteration: ",Dumper($init_params) if $debug;
		# DO NOT USE THE DELTA-A TEST FOR ENDING ITERATION
		# Stop iterating when maximum residual is below peakmin * max y,
		# or when max residual is within peakmin of previous max residual
		# my $a_init = $$init_params{'a_init'};
		# if ($iteration == 1) { 
		# 	$a_init_1 = $a_init;
		# 	$a_init_n = $a_init; 
		# } elsif ($peakmin > 0 && 
		# 		($a_init < $peakmin * $a_init_1 ||
		# 		abs($a_init - $a_init_n) < $peakmin * $a_init_n)) {
		# 	$peak_count = $iteration - 1;
		# 	last;
		# } else { 
		# 	$a_init_n = $a_init;
		# }
		# printf STDERR "Iteration %i: initial height = %.3f\n",$iteration,$a_init;
		
		# capture makes gnuplot shut up.  otherwise, we'd get a ton of output from gnuplot to STDERR
		if (! $gnuin) {
		    print "Gaussian fit for $iteration peaks\n" if $debug;
		    capture {
			_gaussian(	{	
						input_file	=> $gp_data_file,
						params_file	=> $gp_params_file,
						peak_count	=> $iteration,
						peakmin		=> $peakmin,
						name		=> $name,
						verbose		=> $verbose,
						debug		=> $debug,
						file_stem	=> $file_stem,
						slopeflag	=> $slopeflag,
						plot		=> $plot,
					}
			)
		    } ; 

		}

		# Append solutions to hash one at a time, to make it easier
		# to stop when delta Rabs falls below peakmin. 
		# Use _sec_fit_append_return_hash, NOT _sec_fit_build_return_hash

		my $new_model = _sec_fit_append_return_hash ( {
						x		=> $x,
						y		=> $y,
						iteration	=> $iteration,
						name		=> $name,
						verbose		=> $verbose,
#						pooled		=> $p,
						sum_y		=> $sum_y,
						file_stem	=> $file_stem,
						debug		=> $debug,
						directory	=> $working_directory,
						peakmin		=> $peakmin,
			} );

		if (! $new_model) {
			$peak_count = $iteration - 1;
			print "No model found for iteration $iteration!\n" if $verbose;
			last;
		}
		print "New model returned: ",Dumper($new_model) if $debug;
		# $rabs_old = $all_models{models}{$iteration-1}{R_abs} if ($iteration>1);
		my $rabs = $$new_model{R_abs};
		my $d_rabs = $rabs_old - $rabs;
		# Do NOT stop iterating if Rabs change is less than peakmin,
		# since we're looking for the LAST model with delta Rabs > peakmin
		# Stop only if the change is less than peakmin/10, if peakmin<>0;
		# don't stop even then if in verbose mode.
		
		# Stop and ignore the last model if it includes ridiculously wide peaks
		# e.g. sigma > fraction range/2
		my $c1 = $$new_model{'solution'}->{'c1'};
		
		if ($c1 > $$fraction[-1] / 2 ) {
			printf "STOPPING ITERATION: Peak in round $iteration gave has wide sigma %.0f\n",$c1;
			$peak_count = $iteration -1;
			last;
		}
		
#		if ($d_rabs < $peakmin && $iteration > 1) 
		if ($peakmin && $d_rabs < $peakmin / 10 && $iteration > 1) 
		
		{
			printf "DELTA R_ABS < PEAK MIN / 10:\n".
				"\tR_abs from %i to %i peaks = %.4f - %.4f ".
				"= %.4f < peakmin/10 = %.4f\n",
				$iteration-1, $iteration, $rabs_old, $rabs,
				$rabs_old - $rabs, $peakmin/10 if $verbose;

			# Stop iterating if delta Rabs < peakmin,
			# or if delta Rabs < peakmin/10 in verbose mode

			$peak_count = $iteration - 1 if (! $verbose);
						# || $d_rabs < $peakmin/10;

			# Do NOT append the new model if delta Rabs < peakmin/10,
			# unless we're in verbose mode
			
			last if ! $verbose;

		} 
		if ($rabs < $peakmin) {
			printf "R_ABS < PEAK MIN\n".
				"\tR_abs for %i peaks = %.4f < peakmin = %.4f\n",
				$iteration, $rabs, $peakmin if $verbose;
			# Stop iterating if Rabs < peakmin, or
			# Rabs < peakmin/10 in verbose mode.
			# Since Rabs > 0, further delta Rabs must be < peakmin.
			$peak_count = $iteration if (! $verbose) ||
							$rabs < $peakmin/10;
		}
		$rabs_old = $rabs;
		$$new_model{solution}{slope} = 0 if ! $slopeflag;
		$all_models{model}{$iteration} = $new_model;
	}
	$all_models{max_peaks} = $peak_count;
	return \%all_models;
}

# For iteration n, calculate files for residual to be saved to "name_n_residual.dat" 
sub _sec_fit_append_return_hash  {
	my %params = validate( @_, {	
					x			=> { type => SCALARREF },
					y			=> { type => SCALARREF },
#					pooled		=> { type => SCALARREF },
					sum_y		=> { type => SCALAR },
					name		=> { type => SCALAR },
					file_stem	=> { type => SCALAR },
					iteration	=> { type => SCALAR },
					directory	=> { type => SCALAR },
					debug		=> { type => SCALAR, 
								default => 0 },
					verbose		=> { type => SCALAR, 
								default => 0 },
					peakmin		=> { type => SCALAR, 
								default => .05 },
				}
			);
	my $x		= $params{'x'};
	my $y		= $params{'y'};
#	my $pooled	= $params{'pooled'};
	my $sum_y	= $params{'sum_y'};
	my $name	= $params{'name'};
	my $file_stem	= $params{'file_stem'};
	my $iteration	= $params{'iteration'};
	my $verbose	= $params{'verbose'};
	my $debug	= $params{'debug'};
	my $directory	= $params{'directory'};
	my $peakmin	= $params{'peakmin'};

	my $ny		= nelem($y);
	my $y_mean	= sum($y)/$ny;
#	my $y_minus_mean= $y - ($y_mean / $ny);
#	my $variance	= sum( $y_minus_mean * $y_minus_mean ) / $ny;
#	$all_models{'variance'} = $variance;
	return 0 if ! nelem($x);
	my $frxn_factor = (max($x)-min($x)) / nelem($x);
	my $sq2p = sqrt( 8 * atan2(1, 1) );

	print "Iteration $iteration \n" if $verbose;
	my $residual_file	= _generate_gnuplot_datafilename( { stem => $file_stem , version =>  ($iteration) , type => "residual" });
	my $background_file = _generate_gnuplot_datafilename( { stem => $file_stem , version => ($iteration) , type => "background" });
	my $fitted_file = _generate_gnuplot_datafilename( { stem => $file_stem , version => ($iteration) , type => "fitted" });
	my $parameter_file = _generate_gnuplot_datafilename( { stem => $file_stem , version => ($iteration) , type => "parameter" });

	return 0 if (! -e $residual_file);
	my $y_residual	= rcols( $residual_file,  1 , { EXCLUDE => '/^#/' }, { INCLUDE => '/s*\d/' } );
	my $background	= rcols( $background_file, 1, { EXCLUDE => '/^#/' }, { INCLUDE => '/s*\d/' });
	my $y_fitted	= rcols( $fitted_file, 1, { EXCLUDE => '/^#/' }, { INCLUDE => '/s*\d/' });

	my $observed_points	= nelem( $y_residual );
	my $sum_abs_residual	= sum( abs( $y_residual ) );
	my $sum_background	= sum( $background );
	my $R_abs		= $sum_abs_residual / ($sum_y - $sum_background);
	my $total_area		= ($sum_y - $sum_background) * $frxn_factor;
#		my $R_abs_pooled	= 1;

	my $RMSD	= sqrt( sum( $y_residual * $y_residual ) / $observed_points );
	my $mean	= ( $sum_y - $sum_background) / $observed_points;
	my $RMSD_div_mean = $RMSD / $mean;


	# Using the $params and $values lists, create a hash with the keys in $params and values in $values.
	# We had some oddness with rcols reading in the wrong columns of text
	my ($values, $params) = rcols($parameter_file, {PERLCOLS => [0] }, 2);

	# Sort solution peaks in order by "start" = b - 2c
	
	my %solution_unsorted;
	@solution_unsorted{@$params} = list( $values );
	my $sort_result = sort_peaks( { solution => \%solution_unsorted, 
					npeaks	 => $iteration } );
	my $ndx_ref = $$sort_result{'index'};
	my @ndx = @$ndx_ref;
	print "Peak order: ",join(",",@ndx),"\n" if $debug;
	# Keep original background values
	my %solution = (
			slope	  => $solution_unsorted{'slope'},
			intercept => $solution_unsorted{'intercept'}, 
			);
	my $min_area = 1;
	# Add peaks in order by start, and find min. peak area
	foreach my $ip (1..$iteration) {
		my $i = $ndx[$ip];
		foreach my $j ('a', 'b', 'c') {
			$solution{"$j$ip"} = $solution_unsorted{"$j$i"};
		}
		my $area = abs($solution{"a$ip"} * $solution{"c$ip"})*$sq2p
					/ ($total_area ? $total_area : 1);
		$min_area = $area if $area < $min_area;
	}

	print "Solution: m=$solution{'slope'}, b=$solution{'intercept'}\n".
		"\ta1=$solution{'a1'},b1=$solution{'b1'},c1=$solution{'c1'}\n"
		if $debug;

	if ( $verbose ) {
		print "\tSum( abs( residual) )		$sum_abs_residual	\n";
		print "\tSum( background )		$sum_background	\n";
		print "\ttotal_area			$total_area\n";
		print "\tR_abs				$R_abs\n";
		print "\tRMSD				$RMSD\n";
		print "\tmean				$mean\n";
		print "\tRMSD/mean			$RMSD_div_mean\n";
		print "\tSolution:",Dumper(\%solution);
	}
	
	# DO NOT USE THE FOLLOWING MIN. PEAK AREA TEST
	# Stop appending to results hash if min. peak < pean_min * total area
	# if ($min_area < $peakmin) {
	# 	$all_models{'max_peaks'} = $iteration;
	# 	printf "Stopping iteration: min peak area for model %i ".
	# 		"= %.4f < peakmin = %.4f\n", $iteration,
	# 		$min_area, $peakmin if $verbose; 
	# 	$peak_count = $iteration;
	# }

	# Return model for the iteration to be appended

	my $model = 	{ 
				solution	=> { %solution },
				y_fitted	=> [ list( $y_fitted ) ] ,
				y_residual	=> [ list( $y_residual ) ],
				R_abs		=> $R_abs,
				RMSD		=> $RMSD,
				mean		=> $mean,
			};
	print "New model to be returned: ",Dumper($model) if $debug > 1;
	return $model;
}


# USE _sec_fit_append_return_hash ABOVE, NOT THE FOLLOWING ROUTINE
# For iteration n, calculate files for residual to be saved to "name_n_residual.dat" 
sub _sec_fit_build_return_hash  {
	my %params = validate( @_, {	
					x			=> { type => SCALARREF },
					y			=> { type => SCALARREF },
#					pooled		=> { type => SCALARREF },
					sum_y		=> { type => SCALAR },
					name		=> { type => SCALAR },
					file_stem	=> { type => SCALAR },
					peak_count	=> { type => SCALAR },
					directory	=> { type => SCALAR },
					debug		=> { type => SCALAR, default => 0 },
					verbose		=> { type	=> SCALAR, default => 0 },
					peakmin	=> { type => SCALAR, default => .05 },			
				}
			);
	my $x		= $params{'x'};
	my $y		= $params{'y'};
#	my $pooled	= $params{'pooled'};
	my $sum_y	= $params{'sum_y'};
	my $name	= $params{'name'};
	my $file_stem	= $params{'file_stem'};
	my $peak_count	= $params{'peak_count'};
	my $verbose	= $params{'verbose'};
	my $debug	= $params{'debug'};
	my $directory	= $params{'directory'};
	my $peakmin	= $params{'peakmin'};

	my %all_models	= ( 	max_peaks 	=> $peak_count,
				directory	=> $directory,

				);
	my $ny		= nelem($y);
	my $y_mean	= sum($y)/$ny;
#	my $y_minus_mean= $y - ($y_mean / $ny);
#	my $variance	= sum( $y_minus_mean * $y_minus_mean ) / $ny;
#	$all_models{'variance'} = $variance;
	return 0 if ! nelem($x);
	my $frxn_factor = (max($x)-min($x)) / nelem($x);
	my $sq2p = sqrt( 8 * atan2(1, 1) );

	my $iteration = 0;
	while ( $iteration < $peak_count ) {
		$iteration++;
		print "Iteration $iteration \n" if $verbose;
		my $residual_file	= _generate_gnuplot_datafilename( { stem => $file_stem , version =>  ($iteration) , type => "residual" });
		my $background_file = _generate_gnuplot_datafilename( { stem => $file_stem , version => ($iteration) , type => "background" });
		my $fitted_file = _generate_gnuplot_datafilename( { stem => $file_stem , version => ($iteration) , type => "fitted" });
		my $parameter_file = _generate_gnuplot_datafilename( { stem => $file_stem , version => ($iteration) , type => "parameter" });

		my $y_residual	= rcols( $residual_file,  1 , { EXCLUDE => '/^#/' }, { INCLUDE => '/s*\d/' } );
		my $background	= rcols( $background_file, 1, { EXCLUDE => '/^#/' }, { INCLUDE => '/s*\d/' });
		my $y_fitted	= rcols( $fitted_file, 1, { EXCLUDE => '/^#/' }, { INCLUDE => '/s*\d/' });

		my $observed_points	= nelem( $y_residual );
		my $sum_abs_residual	= sum( abs( $y_residual ) );
		my $sum_background	= sum( $background );
		my $R_abs		= $sum_abs_residual / ($sum_y - $sum_background);
		my $total_area		= ($sum_y - $sum_background) * $frxn_factor;
#		my $R_abs_pooled	= 1;

		my $RMSD	= sqrt( sum( $y_residual * $y_residual ) / $observed_points );
		my $mean	= ( $sum_y - $sum_background) / $observed_points;
		my $RMSD_div_mean = $RMSD / $mean;


		# Using the $params and $values lists, create a hash with the keys in $params and values in $values.
		# We had some oddness with rcols reading in the wrong columns of text
		my ($values, $params) = rcols($parameter_file, {PERLCOLS => [0] }, 2);

		# Sort solution peaks in order by "start" = b - 2c
		
		my %solution_unsorted;
		@solution_unsorted{@$params} = list( $values );
		my $sort_result = sort_peaks( { solution => \%solution_unsorted, 
						npeaks	 => $iteration } );
		my $ndx_ref = $$sort_result{'index'};
		my @ndx = @$ndx_ref;
		print "Peak order: ",join(",",@ndx),"\n" if $debug;
		# Keep original background values
		my %solution = (
				slope	  => $solution_unsorted{'slope'},
				intercept => $solution_unsorted{'intercept'}, 
				);
		my $min_area = 1;
		# Add peaks in order by start, and find min. peak area
		foreach my $ip (1..$iteration) {
			my $i = $ndx[$ip];
			foreach my $j ('a', 'b', 'c') {
				$solution{"$j$ip"} = $solution_unsorted{"$j$i"};
			}
			my $area = abs($solution{"a$ip"} * $solution{"c$ip"})*$sq2p
						/ $total_area;
			$min_area = $area if $area < $min_area;
		}

		print "Solution: m=$solution{'slope'}, b=$solution{'intercept'}\n".
			"\ta1=$solution{'a1'},b1=$solution{'b1'},c1=$solution{'c1'}\n"
			if $debug;

		if ( $verbose ) {
			print "\tSum( abs( residual) )		$sum_abs_residual	\n";
			print "\tSum( background )		$sum_background	\n";
			print "\ttotal_area			$total_area\n";
			print "\tR_abs				$R_abs\n";
			print "\tRMSD				$RMSD\n";
			print "\tmean				$mean\n";
			print "\tRMSD/mean			$RMSD_div_mean\n";
			print "\tSolution:",Dumper(\%solution);
		}
		
		# DO NOT USE THE FOLLOWING MIN. PEAK AREA TEST
		# Stop compiling results if min. peak < pean_min * total area
		if (0) {
			if ($min_area < $peakmin) {
				$all_models{'max_peaks'} = $iteration;
				printf "Stopping iteration: min peak area for model %i ".
					"= %.4f < peakmin = %.4f\n", $iteration,
					$min_area, $peakmin if $verbose; 
				$peak_count = $iteration;
			}
		}		
		$all_models{model}{$iteration} = { 
					solution	=> { %solution },
					y_fitted	=> [ list( $y_fitted ) ] ,
					y_residual	=> [ list( $y_residual ) ],
					R_abs		=> $R_abs,
					RMSD		=> $RMSD,
					mean		=> $mean,
		};
	}
	return \%all_models;

	#print Dumper(\%all_models);
}


# For iteration n, write out to file "name_n_parameter.dat"
sub _setup_parameters {
	my %params = validate( @_, {	
					x		=> { type => SCALARREF },
					y		=> { type => SCALARREF },
					sum_y		=> { type => SCALAR },
					file_stem	=> { type => SCALAR },
					name		=> { type => SCALAR},
					iteration	=> { type => SCALAR },
					verbose		=> { type => SCALAR, default => 0 },
					debug		=> { type => SCALAR, default => 0 },
					gnuin		=> { type => SCALAR, default => 0 },
					slopeflag	=> { type => SCALAR, default => 0 },
					width	=> { type => SCALAR, default => 0 },
				}
			);
	my $x		= $params{'x'};
	my $y		= $params{'y'};
	my $sum_y	= $params{'sum_y'};
	my $file_stem	= $params{'file_stem'};
	my $name	= $params{'name'};
	my $iteration	= $params{'iteration'};
	my $verbose	= $params{'verbose'};
	my $debug	= $params{'debug'};
	my $gnuin	= $params{'gnuin'};
	my $slopeflag	= $params{'slopeflag'};
	my $width	= $params{'width'};
	$width = 4.0 if ! $width;
	my $parameters	= '';
	my %init_params = ();

	my $filename = _generate_gnuplot_datafilename( { stem => $file_stem, version => $iteration, type => "parameter" } );



	if ($debug) {
		print "\n***************** $iteration Gaussian Initial Parameters ***********\n";
	}

	# for the first iteration, it is rather simple!  we just take the max peak point
	if ( $iteration == 1 ) {
		# initial parameters
		
		my $a_init = max $y;

		# the the mean, or mu, where the max peak is.  The center of the peak.
		my $tmp_index = which( $y == $a_init );
		my $tmp_b = $x->at( $tmp_index->at(0)  );
		my $b_init = $tmp_b; 

		# estimate the the minimum width at half max.  

		# my $c_init = sqrt( _width_of_greatest_peak(	
		#				x		=> $x,
		#				y		=> $y,
		#				y_min	=> $a_init / 2,
		#				verbose	=> $verbose,
		#			));
		
		my $c_init = 1;
		my $m = 0;
		my $intercept = 0;
		$parameters = <<"EOL";
intercept	= $intercept
a1		= $a_init
b1		= $b_init
c1		= $c_init
EOL
		$parameters = $parameters."slope = $m \n" if $slopeflag;
		if ($debug) {
			print "\tFirst iteration, setting up file $filename with max peak\n";
			#print "index is $tmp_index\n";
			#print "index_2 is $tmp_index_2\n";
			print "\ta = $a_init\n";
			print "\tb = $b_init\n";
			print "\tc = $c_init\n";
		}
		if (! $gnuin) {
			my $fh = new IO::File;
			$fh->open($filename, "w");
			print $fh  $parameters;
			$fh->close;
		}
		%init_params = ( a_init => $a_init, b_init => $b_init, c_init => $c_init );

	} else {
		# we have an iteration step beyond 1, so we must build on the previous solved parameters
		# determine the files with the previous solution, and load the data into PDLs
		my $previous_iteration = $iteration - 1;
		my $residual_file	= _generate_gnuplot_datafilename( { stem => $file_stem , version =>  ($previous_iteration) , type => "residual" });
		my $background_file	= _generate_gnuplot_datafilename( { stem => $file_stem , version => ($previous_iteration) , type => "background" });
		my $solved_params_file	= _generate_gnuplot_datafilename( { stem => $file_stem , version => ($previous_iteration) , type => "parameter", extension => ".dat" });

		# Load the values, stored in the 2nd column of the files into PDLs
		#print "residual $residual_file\n";
		my $residual	= rcols( $residual_file,  1 , { EXCLUDE => '/^#/' }, { INCLUDE => '/s*\d/' } );
		#print "background\n";
		my $background	= rcols( $background_file, 1, { EXCLUDE => '/^#/' }, { INCLUDE => '/s*\d/' });
		#print "solved_parameters\n";
		my $solved_parameters	= rcols( $solved_params_file, 2 );
		#print "observed_points\n";
		my $observed_points		= nelem( $residual );

		my $sum_abs_residual	= sum( abs( $residual ) );
		my $sum_background		= sum( $background );
		my $R_abs	= $sum_abs_residual / ($sum_y - $sum_background);
		my $RMSD	= sqrt( sum( $residual * $residual ) / $observed_points );
		my $mean	= ( $sum_y - $sum_background) / $observed_points;

		# Find max of running average of $width fraction's worth of points
		my $points_to_smooth = int($width * $observed_points / 
				($x->at($observed_points - 1) - $x->at(0)));
		# Always include center point by making number of points odd
		$points_to_smooth += 1 if ! $points_to_smooth % 2;
		my $smooth_residual = _smooth({	x => $residual,
						n => $points_to_smooth,
						# width => $width,
						type => "Gaussian",
						debug => $debug });
		# my $smooth_residual = $$smooth_result{'smoothed'};
		print "returned from _smooth\n" if $debug;
		print "smoothed resid pdl: ",Dumper($smooth_residual) if $debug;
		if ($verbose > 1) {
			print "Smoothed resid\n";
			for (my $i=0; $i<$observed_points; $i++) {
				printf "%.4f\n", at($smooth_residual,$i);
			}
		}

		# print "smoothed resid ref: ",Dumper($$smooth_residual) if $debug;
		my $a_init = max $smooth_residual;
		my $tmp_index = which( $smooth_residual == $a_init );
		my $b_init = at( $x, at($tmp_index,0) );
		
		print "Max smoothed residual = $a_init at ",at($tmp_index,0),
			" = frxn $b_init\n" if $debug;
		print "Initial peak estimate: frxn $b_init, height $a_init\n" if $verbose;
#		my $residual_greater_half = $residual * ($residual > $a_init / 2);

#		my $top = sum( $residual_greater_half * $x);
#		my $bottom  = sum( $residual_greater_half );

#		my $b_init = $top / $bottom;
		
		my $c_init = 1;
		# my $c_init = sqrt(_width_of_greatest_peak(
		#					x		=> $x,
		#					y		=> $y,
		#					y_min	=> $a_init / 2,
		#					verbose	=> $verbose,
		#				));
		# $c_init = 1;

		if ($debug) {
			print "\tBuild on the previous solved parameters in $solved_params_file\n";
			print "\tsum of the abs of the residual is $sum_abs_residual\n";
			print "\tsum of the backgroundis $sum_background\n";
			#print "\n\tPrevious params are: $solved_parameters with " .  nelem( $solved_parameters ) . " elements\n";

			#print "\tResiduals have: " . nelem( $residual ) . " elements\n";

			print "\tR_abs	= $R_abs\n";
			print "\tRMSD	= $RMSD\n";
			print "\tmean	= $mean\n";
			print "\tRMSD/mean = " . $RMSD / $mean . "\n";

			# initial parameters
			#print "\t $residual_greater_half\n";

			#print "\ttop is $top\n";
			#print "\tbottom is $bottom\n";

			#print "\tcopy $solved_params_file to file $filename and append the following initial params\n";
			print "\n";
			print "\ta$iteration = $a_init\n";
			print "\tb$iteration = $b_init\n";
			print "\tc$iteration = $c_init\n";
		}
		if (! $gnuin) {
			copy($solved_params_file, $filename) or die "Copy failed: $!";
			my $fh = new IO::File $filename, O_WRONLY | O_APPEND;
			$parameters = <<"EOL";
a$iteration		= $a_init
b$iteration		= $b_init
c$iteration		= $c_init
EOL
			print $fh  $parameters;
			$fh->close;
		}
		%init_params = ( a_init => $a_init, b_init => $b_init, c_init => $c_init );
	}
	return \%init_params;
}

# Take running average of N point of X
sub _smooth {
	my %params = validate( @_, {	
					x => { type => SCALARREF },
					n => { 	type => SCALAR, 
							DEFAULT => 20},
				#	width => { type => SCALAR, DEFAULT => 1.0},
					type => { type => SCALAR, 
							DEFAULT => "Gaussian"},
					debug => { type => SCALAR, 
							DEFAULT => 0},
				} );
	my $x = $params{'x'};
	my $n = $params{'n'};
	my $stype = $params{'type'};
	# my $width = $params{'width'};
	my $gflag = lc($stype) =~ /^g/;
	my $debug = $params{'debug'};

	# Copy positive X's into Y from n/2 to len + n/2 - 1 = len places:
	# E.g. len = 10, n=5, n2 = 2, copy x(0 to 9) to y(2 to 11)

	my $len = nelem($x);
	my $n2 = int($n/2);
	my $ln2m1 = $n2 + $len - 1;
	my $y = zeroes($len + $n - 1);
	my $ycenter = slice($y, "$n2:$ln2m1");

	# Use only positive X's since negative X's are either background
	# fluctuations or over-estimates from another peak.
	
	my $w0 = zeroes($ycenter);
	# give values < 0 weight w0
	$ycenter .= $x * ($x > 0.0 - $w0*(1-$x>0.0));
	
	my $z = zeroes($len);
	my $w = zeroes($n);
	$w .= 1.0 / $n;
	if ($gflag) {
		# Use Gaussian weighting: exp(-x^2/2), x=-2 to 2 = 95% of peak
		# gx = -1 to 1 in n steps
		my $gx = (sequence($n)*2.0/($n-1.0) - 1.0)*2.0;
		my $gy = exp(-($gx*$gx)/2.0);
		$w = $gy / sum $gy;
		print "Gaussian weights:",$w,"\n" if $debug > 1;
	}
	
	my $j = 0;
	for (my $i=0; $i<$n; $i++) {
		$j = $i + $len - 1;
#		print "_smooth 3: i,j=$i,$j\n";
		$z = $z + slice( $y, "$i:$j") * $w->at($i);
		print "_smooth $i: z(20)=",at($z,20),"\n" if $debug > 2;
	}
	print "_smooth: ",Dumper($z) if $debug > 2;
	print "_smooth: ",slice($z,'0:20'),"\n" if $debug > 1;
#	my %return_hash = ( smoothed => $z );
	return $z;
}

# call gnuplot to fit, calculate residual, and perhaps even plot the fit vs the observed
sub _gaussian {

	my %params = validate( @_, {	
#					x			=> { type => SCALARREF },
#					y			=> { type => SCALARREF },
					peak_count	=> { type => SCALAR },
					peakmin	=> { type => SCALAR },
					file_stem	=> { type => SCALAR },
					name		=> { type => SCALAR	},
					input_file	=> { type => SCALAR	},
					params_file	=> { type => SCALAR	},
#					pooled		=> { type => ARRAYREF },
					verbose		=> { type	=> SCALAR, default => 0 },
					debug		=> { type => SCALAR, default => 0 },
					fit_limit	=> { type => SCALAR, default => "5e-6"},
					fit_maxiter	=> { type => SCALAR, default => "500"},
					slopeflag	=> { type => SCALAR, default => 0 },
					plot		=> { type => SCALAR, default => 0 },
				} );

	# Pull arguments from hash
#	my $x			= $params{'x'};
#	my $y			= $params{'y'};
#	my $pooled		= $params{'pooled'};
	my $verbose		= $params{'verbose'};
	my $debug		= $params{'debug'};
	my $file_stem	= $params{'file_stem'};
	my $name		= $params{'name'};
	my $peak_count	= $params{'peak_count'};
	my $peakmin	= $params{'peakmin'};
	my $input_file	= $params{'input_file'};
	my $fit_limit	= $params{'fit_limit'};
	my $fit_maxiter	= $params{'fit_maxiter'};
	my $params_file	= $params{'params_file'};
	my $slopeflag	= $params{'slopeflag'};
	my $plot	= $params{'plot'};

	# Overwrite any existing files
	
	my $gp_fit_log		= _generate_gnuplot_datafilename( { stem => $file_stem , type => 'fit' , extension => ".log" } );
	my $gp_graph_ps		= _generate_gnuplot_datafilename( { stem => $file_stem , version => $peak_count, extension => ".ps" } );
	my $gp_graph_2_ps	= _generate_gnuplot_datafilename( { stem => $file_stem , version => $peak_count, type =>  "solution", extension =>  ".ps" } );
	my $gp_params_file	= $params_file;
	my $gp_fitted		= _generate_gnuplot_datafilename( { stem => $file_stem , version => $peak_count, type => "fitted" } );
	my $gp_residual		= _generate_gnuplot_datafilename( { stem => $file_stem , version => $peak_count, type => "residual" } );
	my $gp_residual_initial		= _generate_gnuplot_datafilename( { stem => $file_stem , version => $peak_count, type => "residual_initial" } );
	my $gp_background	= _generate_gnuplot_datafilename( { stem => $file_stem , version => $peak_count, type => "background" } );

	print "gp params $gp_params_file  gp fitted $gp_fitted\n" if $debug;

	# gnuplot parameters
	my $gp_title	= "Batch $name fit with $peak_count Gaussians";
	my $gp_xlabel	= "Fraction";
	my $gp_ylabel	= "absorption (mAu)";
	my $gp_image_width = 640;
	my $gp_image_height = 400;
	my $gp_function	= "g$peak_count";

	open ( GNUPLOT, "|$GNUPLOT") or die "\tCan't open $GNUPLOT, please make sure that gnuplot is in your path! $1";

#fit g1(x) "y.dat" via "a.dat"
#plot "t_1.dat" using 1:(g1($1)-$2), f1(x), "t_1.dat" to get the residuals....


	# commands sent to gnuplot
	# set up the gaussians in fit command does the actual fit
	my $bg_command = $slopeflag	? "background(x) = slope*x + intercept\n" 
					: "background(x) = intercept\n";
	my $fit_command = "";
	my $g_string = "g0(x) = background(x)\n";
	foreach my $i ( 1..$peak_count ) {
		$fit_command .= "gaussian_$i(x) = abs(a$i) * exp(- (x-b$i)*(x-b$i)/ ( 2*c$i*c$i ) )\n";

# gaussian_1(x) = abs(a1) * exp(- (x-b1)*(x-b1)/ ( 2*c1*c1 ) )
# gaussian_2(x) = abs(a2) * exp(- (x-b2)*(x-b2)/ ( 2*c2*c2 ) )
# gaussian_3(x) = abs(a3) * exp(- (x-b3)*(x-b3)/ ( 2*c3*c3 ) )
# gaussian_4(x) = abs(a4) * exp(- (x-b4)*(x-b4)/ ( 2*c4*c4 ) )
# gaussian_5(x) = abs(a5) * exp(- (x-b5)*(x-b5)/ ( 2*c5*c5 ) )
# gaussian_6(x) = abs(a6) * exp(- (x-b6)*(x-b6)/ ( 2*c6*c6 ) )
# gaussian_7(x) = abs(a7) * exp(- (x-b7)*(x-b7)/ ( 2*c7*c7 ) )

# df(x) = -((x-b ) / (c*c) )* abs(a) * exp( - (x-b)*(x-b)/ ( 2*c*c ) )

		$g_string .= "g$i(x) = g".($i-1)."(x) + gaussian_$i(x)\n";

# g1(x) = background(x) + gaussian_1(x);
# g2(x) = background(x) + gaussian_1(x) + gaussian_2(x);
# g3(x) = background(x) + gaussian_1(x) + gaussian_2(x) + gaussian_3(x);
# g4(x) = background(x) + gaussian_1(x) + gaussian_2(x) + gaussian_3(x) + gaussian_4(x);
# g5(x) = background(x) + gaussian_1(x) + gaussian_2(x) + gaussian_3(x) + gaussian_4(x) + gaussian_5(x);
# g6(x) = background(x) + gaussian_1(x) + gaussian_2(x) + gaussian_3(x) + gaussian_4(x) + gaussian_5(x) + gaussian_6(x);
# g7(x) = background(x) + gaussian_1(x) + gaussian_2(x) + gaussian_3(x) + gaussian_4(x) + gaussian_5(x) + gaussian_6(x) + gaussian_7(x);

	}
	$fit_command .= $g_string;
	$fit_command .= << "EOL";
FIT_LIMIT = $fit_limit
FIT_MAXITER = $fit_maxiter
set table "$gp_residual_initial"
load "$gp_params_file"
plot "$input_file" using 1:(\$2 - $gp_function(\$1)) 
unset table

set fit logfile "$gp_fit_log"
fit $gp_function(x) "$input_file" via "$gp_params_file"
update "$gp_params_file"
EOL

	$fit_command = $bg_command.$fit_command;

	# This plot command will draw the raw inputs, the fitted curve, and the residual
	my $terminal = 'postscript color portrait';
	if ($plot && $plot !~ /^(1|ps)$/i) {
		$terminal       = $plot;
		if ($plot !~ /^postscript/ && $plot !~ / /) {
			$gp_graph_ps    =~ s/\.ps$/.$plot/g;
			$gp_graph_2_ps  =~ s/\.ps$/.$plot/g;
		}
	}

	my $plot_command = << "EOL";
set title	"$gp_title"
set xlabel	"$gp_xlabel"
set ylabel	"$gp_ylabel"
set border
set noxzeroaxis
set noyzeroaxis
set grid
set terminal $terminal
set output "$gp_graph_ps"
plot  "$input_file" title "Observed", $gp_function(x) title "Fitted", "$input_file" using 1:(\$2 - $gp_function(\$1)) title "Residual"

EOL
	# plot the observed, background, and gaussians that make up the fitted curve
	my $gaussians_only = '';
	my $count = 0;
	while ( $count < $peak_count ) {
		$count++;
		$gaussians_only = $gaussians_only . ", gaussian_$count(x) title \"g$count\"";

	}
	$gp_title =  "Batch $name, solution with background and $peak_count Gaussians";
	my $plot_gaussians = <<"EOL";
set title	"$gp_title"
set xlabel	"$gp_xlabel"
set ylabel	"$gp_ylabel"
set border
set noxzeroaxis
set noyzeroaxis
set grid
set terminal $terminal
set output "$gp_graph_2_ps"
load "$gp_params_file"
plot "$input_file" title "Observed", $gp_function(x) title "Fitted", background(x) title "Background" $gaussians_only
EOL

	# save the new y_fitted to name_n_fitted.dat
	my $save_fitted = << "EOL";
set table "$gp_fitted"
plot "$input_file" using 1:($gp_function(\$1)) title "Y Fitted"

EOL

	# save the residual y_residual to name_n_residual.dat
	my $save_residual = << "EOL";
set table "$gp_residual"
plot "$input_file" using 1:((\$2 - $gp_function(\$1))) title "Residual"
EOL

	


	# save the background to name_n_background.dat
	my $save_background = << "EOL";
set table "$gp_background"
plot "$input_file" using 1:(background(\$1)) title "background"
EOL




	print GNUPLOT $fit_command;

	if ($plot) {
		print GNUPLOT $plot_command;
		print GNUPLOT $plot_gaussians;
	#print $plot_command;
	#print $fit_command;
	}
	print GNUPLOT $save_fitted;
	print GNUPLOT $save_residual;
	print GNUPLOT $save_background;


	close (GNUPLOT);
}



#  make sure that the input arrays are all of the same length.  Also enforce that
#  pooled is only filled booleans, of value 0 or 1
sub _sanity_check {

	my %params = validate( @_ , {
									fraction	=>
										{	type	=> ARRAYREF },
									absorbance =>
										{	type	=> ARRAYREF },
									pooled		=>
										{	type	=> ARRAYREF },
									verbose		=> 
										{	type	=> SCALAR,
											default => 0 },
									debug		=> 
										{	type	=> SCALAR,
											default => 0 },
								} );

	my $fraction	= $params{'fraction'};
	my $absorbance	= $params{'absorbance'};
	my $pooled		= $params{'pooled'};

	my $fraction_length	= @$fraction;
	my $absorbance_length = @$absorbance;
	my $pooled_length = @$pooled;

	if ($params{'debug'} ) {
		print "Sanity Checking that x y pooled are all the same length\n";
		print "\tfraction has $fraction_length elements and absorbance has $absorbance_length elements\n";
	}

	if ( ($fraction_length == $absorbance_length) && ($pooled_length == $fraction_length)) {
		if ( $fraction_length > 0 ) {
			return $fraction_length;
		} 
		warn "Unacceptable inputs.  Data has 0 acceptable inputs";

		exit();
	}

	warn "Arrays are of differing lengths!  Please correct the input columns of fraction and absorbance so that are both the same count.\n";
	warn "fraction\t$fraction_length\n";
	warn "absorbance\t$absorbance_length\n";
	warn "pooled\t$pooled_length\n";
	exit();
}


sub _generate_gnuplot_datafilename {
	my %params = validate( @_ , {
									stem		=> { type => SCALAR },
									type		=> { type => SCALAR, optional => 1, default => ''},
									version		=> { type => SCALAR, optional => 1, default => ''},
									extension	=> { type => SCALAR, optional => 1, default => '.dat' },
									verbose		=> { type => SCALAR, optional => 1, default => 0 },
								} );

	my $stem		= $params{'stem'};
	my $type		= $params{'type'};
	my $version		= $params{'version'};
	my $extension	= $params{'extension'};
	my $verbose		= $params{'verbose'};

	if (grep {defined && length} $version ) {
		$version = "_" . $version;
	}

	if (grep {defined && length} $type ) {
		$type = "_" . $type;
	}


	my $filename = $stem . $version . $type . $extension;
	if ($verbose) {
		# print "### _generate_gnuplot_datafilename() has generated filename is $filename\n";
	}
	return $filename;
}


# return the width of the greatest peak > y_min 
sub _width_of_greatest_peak {
	my %params = validate( @_ , {
									x			=> { type => SCALARREF },
									y			=> { type => SCALARREF },
									y_min		=> { type => SCALAR }, 
									verbose		=> { type => SCALAR, optional => 1, default => 0 },
									debug		=> { type => SCALAR, optional => 1, default => 0 },
								} );

	my $x			= $params{'x'};
	my $y			= $params{'y'};
	my $y_min		= $params{'y_min'};
	my $verbose		= $params{'verbose'};
	my $debug		= $params{'debug'};

	my $y_greater_half = $y * ($y > $y_min);


	my $max			= max( $y_greater_half );
	my $tmp_index	= which( $max == $y_greater_half);

	
	my $max_index	= $tmp_index->at(0);
	my $left_index	= $max_index;
	my $right_index	= $max_index;

	while ($left_index  > $y_min) {
		if( $y_greater_half->at($left_index) <=  0) {
			$left_index++;
			last;
		}
		$left_index--;
	}

	while ($right_index  > $y_min) {
		if( $y_greater_half->at($right_index) <=  0) {
			$right_index--;
			last;
		}
		$right_index++;
	}

	my $width = $x->at($right_index) - $x->at($left_index);

	if ($debug) {
		print "** in _width_of_peak **\n";
		print "\tmax\t$max\n";
		print "\tmax_index\t$max_index\n";
		print "\tleft_index\t$left_index\n";
		print "\tright_index\t$right_index\n";
		print "\ty_greater_half->slice = " .  $y_greater_half->slice("$left_index:$right_index") . "\n";
		print "\tx->slice = " .  $x->slice("$left_index:$right_index") . "\n";
		print "\tx[$left_index] = " . $x->at($left_index) . "\n";
		print "\tx[$right_index] = " . $x->at($right_index) . "\n";
		print "\twidth\t$width\n";
	}
	return $width;
}



#################################################################################
# parse_doc = parse doc.asc file to get flow rate, fraction volume, start & end
# In: lines from <batch>doc.asc file = log of SEC run
# Out: hash with flowrate, frxnvol, starttime, endtime, etc +flow changes
#
# Log should record at least fraction collection start & end, and fraction size
# If units are time, flow rate is also required.
# Since curve file usually records numbers from injection time/volume,
# injection time/volume should be determined if possible and subtracted from
# start and end time/volume.
#################################################################################

sub parse_doc {
	my %params = validate( @_, {	doclines => { type => ARRAYREF },
					verbose	 => { type => SCALAR, default=>0},
					debug	 => { type => SCALAR, default=>0}
				   } );
#	my $lineref	= $params{'doclines'};
	my @lines	= @{$params{'doclines'}};
#	my @lines	= @$lineref;
	my $verbose	= $params{'verbose'};
	my $debug	= $params{'debug'};
	
#	print STDERR "Doc lines:\n|",join("\n",@lines),"|\n" if $verbose;
	
# Sample doc file to be parsed (ignore #):
# Log Book
#	0.00 min  Method Run 1/6/2010, 12:30:00 PM Pacific Standard Time, Method : , Result : c:\...\prime\2010Jan06no001.RES
#	0.00 min  Batch ID: 24E3F234-BE5D-458F-8EAE-097AF690EA3F
#	0.00 min  Base Time {}
#	0.00 min  Concentration 0 %B
#	0.00 min  Autozero
#	0.01 min  Flow 0.5 ml/min
#  200.12 min  Break point 2
#  200.12 min  Autozero
#  200.13 min  Injection Valve Inj
#  360.46 min  Break point 3
#  360.46 min  Fraction size 5.0 ml
#  818.61 min  Break point 4
#  818.61 min  Fraction Off
#  818.61 min  Injection Valve Load
# 1037.36 min  End 1/7/2010, 5:47:09 AM Pacific Standard Time

# Sample of version 2 doc file:
# Logbook
#     0.00 min  Method Run 1/20/2010, 2:26:44 PM, Method : SEC X2 v3, Result : c:\...\Default\SSGCID Results\SEC Purifications\SEC X2 v3055\SEC X2 v3001.res
#     0.00 min  Base CV, 318.56 {ml}, HiLoad_26/60_Superdex_75_prep_grade
#     0.13 min  Flow     1.00 ml/min
#     0.13 min  SampleFlow     10.0 ml/min
#     2.20 min  SampleFlow      1.0 ml/min
#     2.74 min  SampleFlow      0.0 ml/min
#     2.74 min  Flow     1.00 ml/min
#   106.31 min  Flow     1.00 ml/min
#   106.31 min  MethodBase SamplePump
#   106.33 min  Hold 1/20/2010, 4:13:01 PM
#   106.33 min  DirectLoad      1.0 ml/min,     50.0 ml, No ,      0.1 ml/min
#   119.54 min  Block Complete_Sample_Load
#   119.54 min  Base Volume {ml}
#   119.54 min  InjectionValve Inject 
#   119.54 min  SampleFlow      1.0 ml/min
#   129.54 min  InjectionValve Load 
#   129.54 min  SampleFlow      0.0 ml/min
#   129.54 min  Flow     1.00 ml/min
#   129.54 min  Block Volume_Fractionation
#   189.54 min  Fractionation 18mm , 5.000  ml , NextTube , Volume 
#   428.48 min  Block Fractionation_Stop
#   428.48 min  Base SameAsMain
#   428.48 min  FractionationStop 
#   428.50 min  Peak_FracStop 
#   428.52 min  End Block
#   428.52 min  End 1/20/2010, 9:35:20 PM

	my %results = ();
	my $start = 0;
	my $now = 0;
	my $maxtime = 0;
	my $lastdate;

	# Status: 1=initial, 2=inject, 3=collecting fractions, 4=stop collecting,
	#	  5=end after injection, 6=end without injection or fraction
	# Logic: Change status to 1 on injection, to 2 on fraction collection,
	#	from 3 to 4 or from 2 to 5 or 1 to 6 on end i.e. to 7-status.
	# Record flow, time and status of all flow changes before end; 
	# analyze after parsing (see below for that logic).

	my $status = 1;
	my @flowchange = ();	# hash of flow, time, status
	my $flow = 0;
	my $nflows = 0;
	my $unit = 0;
	my $frxnstart = 0;
	my $time_per_frxn = 0;
	my $frxn = 0;
	my $old_frxn = 0;
	
	foreach (@lines) {
		print "Parse doc line \"$_\"\n" if $debug > 2;
		chomp;
		s/\s+$//g;
		my $uc = $_;
		$_ = lc;
		if ($_ =~ /^ *[-0-9]+\.[0-9]* *m[li]/) {
			($now = $_) =~ s/ m[li].*//g;
			$now =~ s/ //g;
			if (! $unit) {
				($unit=$_)=~s/^ *[0-9]+\.[0-9]* *//g;
				$unit =~ s/ .*$//g;
				$unit =~ s/ //g;
				$results{'unit'} = $unit;
				print "Doc unit: $unit\n" if $debug;
			}
			$frxn = $old_frxn + ($now - $frxnstart)/$time_per_frxn 
				if $status == 3 && $time_per_frxn > 0.0;
			print "Now = $now, status = $status, frxn $frxn\n" if $debug > 1;
		}
		if ($_ =~ / method run /) {
			print "Run line: |$uc|\n" if $debug > 1;
			my $sourcefile = $uc;
			$sourcefile =~ s/.*, result[ :]*//gi;
			$results{'sourcefile'} = $sourcefile;
			print "Source file: |$sourcefile|\n" if $debug;
			s/.*run //gi;
			s/m.*/m/gi;
			$results{'startdate'} = $_;
			print "Start date = |$_|\n" if $debug;
		} elsif ($_ =~ / end *[1-9]/ && $_ !~ / end block/) {
			print "End line, frxn $frxn: $uc\n" if $debug > 1;
			s/.* end //g;
			s/m .*/m/g;
			$results{'enddate'} = $_;
			print "End date = $_\n" if $debug;
			$lastdate = $now;
			$status = 7 - $status if ($status < 4);
			if (!grep{ defined && length} $results{'endtime'}) {
				$results{'endtime'} = $now - $start;
				print "End time: $now\n" if $debug;
			}
		} elsif ($_ =~ / flow +[0-9.]+ ml\/min/ &&
			 ! grep{ defined && length} $results{'endtime'}) {
			print "Flow rate line, frxn $frxn: $uc\n" if $debug > 1;
			s/.* flow +//g;
			s/ ml\/min.*//g;
			# Store all flow rate changes 
			push @flowchange, {	oldflow => $flow,
						flow => $_ , 
						time => $now, 
						frxn => $frxn,
						status => $status } 
				if ($nflows && $_ != $flow);
			# Return last non-zero flow up to frxn 5 i.e.
			# status < 3 (no fractions yet) or < 5 frxns collected:
			$results{'flowrate'} = $_ if $_ > 0.00 && 
				(  $frxn < 5.0 || ! $nflows );
			if ($frxn < 5.0 && $status == 3 && $unit !~ /ml/i) {
				# Set frxnstart and time_per_frxn to count frxns 
				# starting from now with the new flow rate:
				$old_frxn = $frxn;
				$frxnstart = $now;
				if ($_ > 0.0 && grep{defined && length} 
				    $results{'frxnvol'}) {
					$time_per_frxn = $results{'frxnvol'}/$_;
				} else {
					$time_per_frxn = 0.0;
				}
				print "Frxn start $frxnstart, time per frxn $time_per_frxn\n" if $debug;
			}
			$flow = $_;
			$nflows++;
			print "Flow rate $nflows = $_\n" if $debug;
		} elsif ($_ =~ / fraction size [0-9.]+ ml/ ||
			 $_ =~ / fractionation .* [0-9.]+ *ml/ ||
			 $_ =~ / fractionation_[0-9a-z]* *[0-9.]+ *ml/) {
			print "Fraction size line: $uc\n" if $debug > 1;
			if ($_ =~ / fractionation .* [0-9.]+ *ml/) {
				s/.* fractionation +//g;
				# print "Fraction size 2a: $_\n" if $debug;
				s/^[0-9.]+ *mm *, *//gi;
				# print "Fraction size 2b: $_\n" if $debug;
			} elsif ($_ =~ / fractionation_.* [0-9.]+ *ml/) {
				# s/.* \([0-9.]+ ml\)/$1/g;
				print "Fractionation_ found\n" if $debug;
				s/.* fractionation_[0-9a-z]* *//g;
				print "Change to '$_'\n" if $debug;
			} else {
				s/.* fraction size //g;
				print "Fraction size 3: $_\n" if $debug;
			}
			s/ ml.*//g;
			s/ //g;
			$results{'frxnvol'} = $_;
			print "Fraction size  = $_ at $now\n" if $debug;
			$results{'starttime'} = $now - $start 
				if ! exists($results{'starttime'});
			$status = 3;
			# time_per_frxn = time (or vol) for 1 fraction
			# For unit min: = size ml/frxn / flow ml/min = frxn/min
			# For unit ml:  = size ml/frxn / 1.
			# Use defaults of 0.5 ml/min, 5 ml/frxn if missing data
			$frxnstart = $now;
			$frxn = 0;
			$time_per_frxn = (($_ > 0.0) ? $_ : 5.0) / 
				(  ($unit =~ /ml/i) ? 1.0 : 
				  (($flow > 0.0 ? $flow : 0.5))  );
			print "Frxn start $frxnstart, time per frxn $time_per_frxn flow $flow\n" if $debug;
		} elsif ($_ =~ / Injection Valve Inj/i || $_ =~ / DirectLoad /i) {
			print "Injection line: $uc\n" if $debug > 1;
			$start = $now;
			print "Injection time = $now\n" if $debug;
			$status = 2 if $status == 1;
		} elsif ($_ =~ / Break point 3/i) {
			print "Break point 3 line: $uc\n" if $debug > 1;
			$results{'starttime'} = $now - $start;
			print "BP 3 time $now, start time $start\n" if $debug;
		} elsif ($_ =~ / Break point 4/i || 
			 $_ =~ / FractionationStop/i ||
			 $_ =~ / Fraction Off/i) {
			print "Break point 4 line, frxn $frxn: $uc\n" if $debug;
			$results{'endtime'} = $now - $start;
			print "BP 4 time $now, end time $results{'endtime'}\n" if
									$debug;
			$status = 7 - $status if ($status < 4);
		}
		$maxtime = $now if ($now =~ /^[0-9.]+$/) && $now > $maxtime;

	}
	$status = 7 - $status if ($status < 4);
	print "Start time: |$results{'starttime'}| or |$start| or 0\n" if $debug;
	$results{'starttime'} = (grep{defined && length} $start) 
		? $start : 0
		if ! grep{defined && length} $results{'starttime'};

	print "End time: $results{'endtime'} = ",$results{'endtime'}+$start,
			" - $start, or (".
			"|$results{'lastdate'}| or ".
			"|$now| or |$maxtime|) - $start\n" if $debug;
	$results{'endtime'} = 
	(
		(grep{defined && length} $results{'lastdate'}) ? $lastdate : 
		(grep{defined && length} $now) ? $now : $maxtime
	) - $start if ! grep{defined && length} $results{'endtime'};
	print "End time-> $results{'endtime'}\n" if $debug;

	# Status: 1=initial, 2=inject, 3=collecting fractions, 4=stop collecting,
	#	  5=end after injection, 6=end without injection or fraction
	# Return all warnings with:
	# 	Final status 4, warning status 3 or (status 2 and unit ml)
	#	Final status 5, warning status 2 (after injection if no frxn)
	#	Final status 6, warning status 1 (all changes if no inj.)
	# I.e. final status + warning status = 7
	# Let the calling program decide whether or not to print these warnings 
	# based on the units used in the curve file - if it's ml, don't print.
	
	if (@flowchange) {
		my @flow_warnings = ();
		my @flowlog = ();
		my @wstatus_text = ("DURING RUN", "AFTER INJECTION", 
					"DURING FRACTION COLLECTION");
		for (my $i=0; $i<$#flowchange; $i++) {
			my $change = $flowchange[$i];
			print "Change $i: ",Dumper($change) if $debug;
			my $oldflow = $$change{'oldflow'};
			my $flow = $$change{'flow'};
			my $time = $$change{'time'};
			my $wstatus = $$change{'status'};
			my $frxn = $$change{'frxn'};
			print "Flow changed to $flow at $time, status $wstatus, frxn $frxn\n"
				if $debug;
			my $rfrxn = int(100.*$frxn+0.5)/100;
			if ($status + $wstatus == 7 || 
			    $status == 4 && $unit =~ /ml/i) {
				push @flow_warnings, "WARNING: FLOW RATE CHANGED ".
				    $wstatus_text[$wstatus-1]." from $oldflow ".
				    "to $flow at $time $unit = frxn $rfrxn";
				push @flowlog, $change;
			}
		}
		$results{'warnings'} = \@flow_warnings if (@flow_warnings);
		$results{'flowlog'} = \@flow_warnings if (@flowlog);
	}
	
	print "Parsing results: ",Dumper(\%results) if $debug;
	
	return (\%results);
}


###############################################################################
#
# sort_peaks - sort peaks by start of peak = center - 2 sigma
#
# In = ref. to solution hash returned by sec_fit: a1,b1,c1 ... an,bn,cn, and N
# Out = array of indeces in order by start of peak = b-2c
#
################################################################################

sub sort_peaks {
	my %params = validate( @_, { solution	=> HASHREF,
				     npeaks	=> SCALAR } );
	my %solution	= %{$params{'solution'}};
	my $n		= $params{'npeaks'};
#	print "sort peaks n = $n\n";
	my @ndx = (0, 1);
	my %results = ( index => \@ndx, min => 100000 );
	if ($n > 1) {
		my @start = (-100000000);
		foreach my $i (1..$n) {
			push @start, $solution{"b$i"} - 2*abs($solution{"c$i"});
		}
		# Sort numbers from 1 to n in order by start position of peak 1 to n
#		print "starts: ",join(', ',@start),"\n";
		@ndx = sort {$start[$a] <=> $start[$b]} 0..$n;
#		print "index: ",join(', ',@ndx),"\n";
		$results{'index'} = \@ndx;
	}
#	print "sort peak results: ",Dumper(%results);
	return (\%results);
}
