#!/usr/bin/perl
#
# tad2csv.pl - unpack ICFF binary output files from Opticon Monitor 2
#
# Usage: perl tad2csv.pl [-s [N]] filename1 [filename2 [filename3 ...]]]
# where -s with no number performs default smoothing,
#       -s followed by a number N takes a running average of N points
#       filename(x) includes full path and the .tad2 extension
#
use warnings;
use Cwd;
my $debug = 0;
#
# Authors: 
#	Frank Zucker, Medical Genomics of Pathogenic Protozoa
# 	Thomas E. Kammeyer, former lead programmer for OM1 & 2
# ICFF format has:
# Last 2 4-byte integers (little endian) = 
#	tocoff = Table Of Contents (TOC) offset 
#	toclen = length of TOC entries = 264 for tad2 files
# At byte [tocoff]: 4 characters, "ICFF" 
# Next byte = 1st byte of 1st TOC entry
# Entries = 255 bytes name, 4 bytes offset, 4 length, 1 flag = 01.
# Name = null-terminated string, possibly followed by data
# Offset > 0 is position in file, offset < 0 is relative position in TOC entry
#    (going back from end of name or from start of name?)
# Data inside TOC may be big-endian integers; the smaller value is always correct
# Dye plate read names are "read%3d" with no 0 padding, e.g. read1 to read225
# Reads are 6161 bytes (b): 17 b time stamp header + 32 b x 96 wells x 2 dyes
# except for resaved files: 25 b time stamp + tag + 4 b avg x 96 wells x 2 dyes
# Reads are across rows, A1, A2, ... A12, ... H1, ... H12 for dye 1, then dye 2.
# The header in unresaved reads have:
#       4 b int tag = 1 ???
#	4 b int step number
#	4 b int cycle number
#	1 b = 0 ???
#	4 b float time stamp
# Resaved read headers have:
#	4 b int tag = 0x0202 = 514 decimal ???
#	1 b = 0 ???
#	4 b number of wells = 96 ???
#	4 b int step number
#	4 b int cycle number
#	4 b = 0 ???
#	4 b float time stamp
# point time in the last 4 bytes (see statraw).
# Each well (unresaved) has: 
#	4 b int N = number of raw reads
#	4 b pad (text string='Opticon Monitor 1 and 2 were brought to you by:Matthew A. Chaboud,Lucas S. Garsha,Earl L. Hathaway,Thomas L. Houser,Troy D. Jacobson,Jacob Johnson,Thomas E. Kammeyer,Timothy C. St. Clair,Jieping Xu."I have graven it within the hills, and my vengeance upon the dust within the rock."'
#	8 b double sum of observations
#	8 b double sumsq for calculating sd = sqrt(sumsq/N - (sum/N)^2)
#	4 b float min observation
#	4 b float max observation
# Resaved wells have just the 4b float average read.
# Entry statraw contains status reports including time stamp and temperature:
# 19-byte header, 36-byte records with:
#	4 b float procedure time
#	4 b float ???
#	4 b float hold time
#	4 b int step
#	4 b int cycle
#	4 b float temp1 (block? sample?)
#	4 b float temp2 (sample? block?)
#	4 b float lid_temp
#	4 b int max_cycle
# Dye calibration entries begin with "dye0:" or "dye1:" have copies of dye
# calibration files and include: Plate (name), Gains, Version, Notes, Dye (name),
# DyeUseName, FileType, Colors (number of channels), Wells (num wells/plate),
# and calibration reads with contents = dye or empty, Temp=20, 40, 60 and 80C, 
# orientation 0 or 180 as [contents][orientation]:[Temp]:PR e.g. empty180:80:PR
# The average of the empty calibration reads at 0 and 180 for each well are 
# subtracted from the average of the dye reads for that well ???
#
# Dye Calibration, for a single dye (dye0 = dye1 = HEX in sample case):
# For each channel, 
#    for each calibration temp 20, 40, 60 80C, 
#	average 0 and 180 degree empty reads = e (empty)
#       average 0 and 180 dye reads = d (dye)
#	subtract empty from dye, d - e = f (fluorescent signal)
#       take average dme of all wells for this temp = f_avg
#       divide f by f_avg to get c20, c40, c60, c80 = calibration
# For each read, normalized value n = raw read / calibration
# For each well, 
#	find minimum n for well = m
#	subtract m from each n to get reported value
#
# Test code for unpacking sample file follows

# Force flushing of STDOUT
$| = 1;
my $tadfile;
$tadfile = 'E:\My Documents\plans\Tm_rev_eng\20070227_Copro_Drug15.tad2' if $debug;
my $smooth = 0;
my $smooth_default = 3;
my $n_args = 1 + $#ARGV;
my $arg = "";
if ($n_args > 0)
{
	# Check for smoothing flag if any arguments are left
	if ($ARGV[0] eq "-s")
	{
		shift;
		$smooth = $smooth_default;
		# Check for smoothing argument = number of points to smooth
		if ($n_args > 1)
		{
			if ($ARGV[0] !~ /\D/)
			{
				$smooth = shift;
			}
		}
	}
}
# Show help message if no arguments given, or if only smoothing arguments given
if ($#ARGV < 0)
{
	printf "\nUSAGE:\n".
	  "    For no smoothing:\n\n\tperl $0 file1 [file2 ...]\n\n".
	  "    or for default smoothing:\n\n\tperl $0 -s file1 [file2 ...]\n\n".
	  "    or to smooth N points\n\n\tperl $0 -s N file1 [file2 ...]\n\n".
	  "where file1, file2 etc. include the file path, name and extension.\n";
	exit (1);
}	
while (@ARGV)
{
	$tadfile = shift;
	printf "Read file $tadfile\nsmooth=$smooth" if $debug;
	my $result = tad2csv($tadfile, $smooth, $debug);
	printf "Result: $result\n" if $debug;
}
printf "DONE\n";
exit (0);

#
# Convert tad to csv:
# 	Read binary ICFF output file from Opticon Monitor 2,
# 	Output comma-separated value file of temperatures and intensities
# Argument 1 = file path and name for tad2 file
# Argument 2 = debug flag
#
sub tad2csv 
{
	# Get arguments tad2 filename and debug flag
	local ($filename,$smooth_n,$dbg) = @_;
	my $dump_raw = 0;
	my $dump_reads = 0;
	# use_both_orientations: 0=use 0 only, 1=both forward, 2=180 reversed
	my $use_both_orientations = 1;
	my $subtract_empty_from_calibration = 1;
	my $divide_by_avg_calibration = 1;
	# sum_channels: 0=use chnl 2 only; 1=sum both; 2=root of sum of squares
	# 3=chnl 2 - chnl 1; 4=chnl2 * cal1 + chnl1 * cal2; 
	# 5=(chnl2/cal2 + chnl1*cal2/cal1^2)/2 = 
	# chnl2/cal2 + (chnl1/cal1)*(cal2/cal1)/2
	my $sum_channels = 0;
	# Divide calibrated read values by avg calibrated read value?
	my $divide_by_well_avg = 0;
	# subtract_baseline: 0=never subtract, 1=obey SubGlobalMin, 2=subtract
	my $subtract_baseline = 1;
	my $dump_calibrations = 0;
	# interpolation=(piece-wise) "linear", "cubic",
	#	"least square" (linear), "quadratic" (least square fit to parabola)
	my $interpolation = "linear";
	my $smooth_method = "avg";
	my $parameters_headers = "";
	my $parameters = "";
	if ($dbg)
	{
		$parameter_headers = "Parameters:,use_both_orientations,".
		 "subtract_empty_from_calibration,divide_by_avg_calibration,".
		 "sum_channels,divide_by_well_avg,subtract_baseline,".
		 "smooth_n,interpolation";
		$parameters = ",$use_both_orientations,".
		 "$subtract_empty_from_calibration,$divide_by_avg_calibration,".
		 "$sum_channels,$divide_by_well_avg,$subtract_baseline,".
		 "$smooth_n,$smooth_method,$interpolation";
		printf "\n$parameter_headers\n$parameters\n";
	}
	# Read binary file into $bin, number of bytes to $nin
	if ($filename !~ /[\\\/]/)
	{
		printf "Cwd: ".cwd()."\n" if $dbg;
		$filename = cwd()."/$filename";
		printf "File + path: $filename\n" if $dbg;
	}
	# $filename =~ s/\//\\/g if $dbg;
	printf "New filename: $filename\n" if $dbg;
	
	local($bin, $nin) = read_binary_file($filename, $dbg);
	if (! $nin)
	{
		printf "Cannot read file $filename\n";
		return ("CANNOT READ $filename");
	}
	my $tocend = $nin - 8;
	# Get output file name - append RE and replacing final .tad[2] with .csv
	my $csv = $filename;
	$csv =~ s/.tad2$/_RE.csv/g;
	$csv =~ s/.tad$/_RE.csv/g;
	$csv =~ s/.csv/_dbg.csv/g if $dbg;
	if ($csv eq $filename) {
		printf "Warning: wrong extension file found in $filename\n";
		$csv = $csv.".csv";
	}
	# Set up header for output file - open and write after processing
	# open (CSV, ">$csv");
	my $well_header = 
		"A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,".
		"B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,".
		"C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,".
		"D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,".
		"E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,".
		"F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,".
		"G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,".
		"H1,H2,H3,H4,H5,H6,H7,H8,H9,H10,H11,H12";
	my $csv_header = "Filename: $csv\n\n" .
			 "Wellset: \"Default\"\n\n".
			 "Step,Cycle,Dye,Temp,$well_header\n";
	$csv_header = "$parameter_headers\n$parameters\n$csv_header" if $dbg;
	printf "$csv\n";
	printf "Unpack $nin bytes from \n$filename to \n$csv\n" if $dbg;
	# Get TOC offset and entry length
	local ($tocoff, $toclen) = unpack("x$tocend  l2",$bin);
	printf "tocoff = $tocoff, toclen = $toclen\n" if $dbg;
	# Get 4 byte constant ICFF at TOC offset
	my $ICFF = unpack("x$tocoff a4",$bin);
	printf "Index start: $ICFF\n" if $dbg;
	if ($ICFF ne "ICFF") 
	{
		printf "ERROR: wrong format file found in $filename\n";
		return (0);
	}
	# Get hash of TOC names and locations
	$, = ',';
	$" = ',';
	local (%itemloc, %itemoff, %itemlen);
	my $loc = $tocoff + 4;
	my @names;
	$#names = ($tocend - $loc) / $toclen - 1;
	printf "$#names toc entries \n" if $dbg;
	my $i = 0;
	while ($loc < $tocend)
	{
		local ($name,$off,$len) = unpack("x$loc a255l2",$bin);
		my $namelen = index($name,"\0");
		my $nametext = substr($name,0,$namelen);
		# Make array to stuff into hash with location, offset, length
		$itemloc{$nametext} = $loc;
		$itemoff{$nametext} = $off;
		$itemlen{$nametext} = $len;
		$names[$i] = $nametext;
		# printf "$i=$nametext," if $dbg;
		$i += 1;
		$loc += $toclen;
	}
#	printf "\n$i names in\n" if $dbg;
	$#names = $i - 1;
	printf "\nGot $i = 0 to $#names names from toc:\n" if $dbg;
#	printf "\nIn string:\n @names\n" if $dbg;
#	printf "\nIn string:\n @names[0 .. ($nnames-1)] \n" if $dbg;
#	printf "\nOut of string:\n'", @names,"'\n" if $dbg;
	my $num_wells  = get_value($itemloc{"nwells" },
				   $itemoff{"nwells" },"N",$bin);
	my $num_colors = get_value($itemloc{"ncolors"},
				   $itemoff{"ncolors"},"N",$bin);
	my $num_reads = $num_wells * $num_colors;
	my $plate_read_len = 17 + 32 * $num_reads;
	my $plate_read_len2 = 25 + 4 * $num_reads;
	my $dye0name = get_value($itemloc{"dye0:Dye"},
				 $itemoff{"dye0:Dye"},
				"a$itemlen{'dye0:Dye'}",$bin);
	$dye0name =~ s/^\s*(.*)[\s*\0*]$/$1/g;
	my $dye1name = "";
	if (defined $itemloc{"dye1:Dye"})
	{
		$dye1name = get_value($itemloc{"dye1:Dye"},
				 $itemoff{"dye1:Dye"},
				"a$itemlen{'dye1:Dye'}",$bin);
		$dye1name =~ s/^\s*(.*)[\s*\0*]$/$1/g;
	}
	if ($dye1name and ($dye0name ne $dye1name))
	{
		printf "\n*** WARNING *** Two different dyes, ".
			"\"$dye0name\" and \"$dye1name\", were used.\n".
		    "This program's calibration assumes one dye, not two.\n";
	}
	printf "\n$num_wells wells, $num_colors dyes = $dye0name,$dye1name,".
		"plate read len=$plate_read_len or $plate_read_len2 bytes\n" if $dbg;
	# Get options such as subtract global minimum 
	my $off = $itemoff{"WellSets"};
	my $len = $itemlen{"WellSets"};
	my $limit = $off + $len;
	while ((substr($bin,$off,2) ne "MC") and ($off < $limit))
	{
		$off ++;
	}
	$i = 4;
	while ((substr($bin,$off+$i,1) ne "\000") and ($off + $i < $limit))
	{
		$i ++;
	}
	my $option_str = substr($bin, $off, $i);
	my @option_list = split "\r\n", $option_str;
	chomp(@option_list);
	printf "option string:\n$option_str\n" if $dbg;
	printf "\noption list:\n@option_list\n" if $dbg;
	local ($option, $tag, $val);
	my %option_hash;
	foreach $option (@option_list)
	{
		my $colon = index $option, ":";
		$tag = substr($option,0,$colon);
		$val = ($colon == length($option) - 1) ? "NONE"
						   : substr $option, $colon+1;
		printf "$tag = '$val'\n" if $dbg;
		$option_hash{$tag} = $val;
	}
	if (defined($option_hash{"SmoothWidth"}) and
		$option_hash{"SmoothWidth"} ne "NONE")
	{
		$smooth_n = $option_hash{"SmoothWidth"};
	}
#	printf "sub global min: ".$option_hash{"SubGlobalMin"}."\n" if $dbg;
	# Set up arrays for empty and dye values at 4 temps, 2 angles for raw
	# @cal = 96 wells * 2 channels * 4 temps
	local (@e_raw, @d_raw, @cal, @ecal);
	$#cal = $num_reads * 4 - 1;
#TEKCHANGE: add ecal
	$#ecal = $#cal;
	$#e_raw = $num_reads * 4 * 2 - 1;
	$#d_raw = $#e_raw;
	local ($content, $angle, $temp, $read_name, $j, $k, $e, $d);
	my @temps = ('20','40','60', '80');
	my @angles = ('0', '180');
	my @contents = ('dye', 'empty');
	
	# Read raw dye and empty calibration data for all temps and angles
	my $n_warnings = 0;
	foreach $content (@contents) 
	{
	    $i = 0;
	    foreach $temp (@temps) 
	    {
		printf "Get raw calibration data for $content $temp: " if $dbg;
		foreach $angle (@angles) 
		{
			$read_name = "dye0:$content$angle:$temp:PR";
			printf "$read_name  " if $dbg;
			$len = $itemlen{$read_name};
			$off = $itemoff{$read_name};
			if ($len != $plate_read_len and 
			    $len != $plate_read_len2 and
			    $n_warnings < 1)
			{
				printf "\n\n*** WARNING *** Calibration plate read".
				  " for $read_name has length $len instead of ".
			    	  "$plate_read_len or $plate_read_len2!\n\n";
			    	$n_warnings ++;
			}
			# Find last point of this read
			$j = $i + $num_reads - 1;
			if ($content eq "dye") 
			{
				# printf "Get plate read at $off\n" if $dbg;
				@d_raw[$i..$j] = get_plate_read($bin, $off, $len);
			}
			else 
			{
				@e_raw[$i..$j] = get_plate_read($bin, $off, $len);
			}
			# printf "Read $read_name $i to $j\n" if $dbg;
			# Increment starting point of next read to 1 past end
			$i = $j + 1;
		}
		printf "\n" if $dbg;
	    }
	}
	printf "Got all calibration data\n" if $dbg;
	if ($dump_reads)
	{
		# Dump read data to _read.csv
		my $read = $csv;
		$read =~ s/.csv$/_read.csv/g;
		open READS,">$read";
		# Print header
		printf READS "Item,Temp,Channel,$well_header\n";
		foreach $content (@contents) 
		{
		    $j = 0;
		    foreach $temp (@temps) 
		    {
			foreach $angle (@angles)
			{
			    printf READS "$content:$angle,$temp,1,";
			    if ($content eq "dye") 
			    {
				printf READS "@d_raw[$j..($j+$num_wells-1)]\n".
					",,2,@d_raw[($j+$num_wells)..
					($j+$num_reads-1)]\n";
			    }
			    else 
			    {
				printf READS "@e_raw[$j..($j+$num_wells-1)]\n".
					",,2,@e_raw[($j+$num_wells)..
					($j+$num_reads-1)]\n";
			    }
			    $j += $num_reads;
			}
		    }
		} 
				
	}
	# Average 0 and 180 orientation reads for dye and empty at each temp
	$i = 0;
	$j = 0;
	foreach $temp (@temps) 
	{
		my @t_sum = (0, 0);
		# Loop through each read = 96 wells * 2 channels
		foreach $k (0..$num_reads-1)
		{
		    $e = $e_raw[$j + $k];
		    $d = $d_raw[$j + $k];
		    if ($use_both_orientations == 1)
		    {
		    	# Average 0 and 180 reads in each well
			$e = ($e + $e_raw[$j + $k + $num_reads])/2;
			$d = ($d + $d_raw[$j + $k + $num_reads])/2;
		    } elsif ($use_both_orientations == 2)
		    {
		    	# Average A1 0 read with H12 180 read, etc.
		    	# First channel k=0 to num_wells-1, reverse=(n-1)-k
		    	# Second channel K=n to 2n-1, reverse    = (2n-1)-k
		    	my $r = ((int($k/$num_wells)+1)*$num_wells - 1) - $k;
#		    	my $r = ((($k<$num_wells)?1:2)*$num_wells - 1) - $k;
		    	
			$e = ($e + $e_raw[$j + $r + $num_reads])/2;
			$d = ($d + $d_raw[$j + $r + $num_reads])/2;
		    }
		    $d -= $e if $subtract_empty_from_calibration;
		    # Store and sum (avg dye - avg empty)
		    $cal[$i + $k] = $d;
#TEKCHANGE: store empty calibration
		    $ecal[$i + $k] = $e;
		    $t_sum[( $k < $num_wells ? 0 : 1 )] += $d;
		}
		# Divide each (dye - empty) by avg (dye - empty) for this temp
		$t_sum[0] /= $num_wells;	# Average for channel 0
		$t_sum[1] /= $num_wells;	# Average for channel 1
		if ($divide_by_avg_calibration)
		{
		    foreach $k (0..$num_reads-1)
		    {
			$cal[$i + $k] /= $t_sum[( $k < $num_wells ? 0 : 1 )];
			# printf "$temp:$k cal=$cal[$i+$k]\n" if $dbg;
		    }
		}
		# Increment starting point of avg by 96 * 2, raw by 96 * 2 * 2
		$i += $num_reads;
		$j += $num_reads * 2;
	}
	printf "Processed calibration data\n" if $dbg;

	# Report raw read / calibration - min(raw / cal) for each well
	# Set up arrays for data reads: 
	#	raw = 96 * 2 * 225 temps, 
	#	well_min = 96 * 2
	#	read_time = ~300 time points from reads
	#	stat_time = ~300 reads from statraw
	#	stat_temp = ~300 values from statraw
	
	local (@raw, @read_time, @read_temp, @read_cycle);
	local (@stat_time, @stat_cycle, @stat_temp);
	my $num_cycles = 300;
	$#raw = $num_reads * $num_cycles - 1;
	$#read_time = $num_cycles - 1;

	# Step 1: get temps so we can calibrate reads
	# 1a: Get times from reads
	$i = 0;
	$off = $itemoff{"read0"};
	$len = $itemlen{"read0"};
	do {
		$read_time[$i] = get_time($bin, $off, $len);
		$i += 1;
		$off = $itemoff{"read$i"};
		$len = $itemlen{"read$i"};
	} while ($off);
	$num_cycles = $i;
	$#read_time = $num_cycles - 1;
	$#read_temp = $num_cycles - 1;
	$#read_cycle = $num_cycles - 1;
	$#stat_time = $num_cycles;
	$#stat_temp = $num_cycles;
	$#raw = $num_reads * $num_cycles - 1 if ($num_cycles > 300);
	
	# 1b: Get times and temps for each cycle in statraw
	
	my $statloc = $itemoff{"statraw"} + 19;
	my $statend = $statloc + $itemlen{"statraw"};
	$i = 0;
	my $old_cycle = 1;
	my $old_temp = 0;
	my $old_time = 0;
	printf "Reading statraw records\n" if $dbg;
	while ($statloc + 36 < $statend)
	{
		# Each statraw record = time, step, cycle, block & sample temp
		# Temperature rises during cycles; reads occur between cycles
		# For each cycle store time + temp of first record of next cycle
		# starting with cycle 1 in step 2
		local($ptime,$step,$cycle,$t1,$t2) = 
					unpack("x$statloc fx8V2f2",$bin);
#TEKCHANGE ?: change if condition
#		if ($step == 2 and $cycle > $old_cycle)
		if ($step == 2)
		{
			$stat_time[$i] = $ptime;
#TEKCHANGE: don't round temperature to 1 decimal place, and use current cycle
#			$stat_temp[$i] = int(10*$t2 + 0.5) / 10;
#			$stat_cycle[$i] = $old_cycle;
			$stat_temp[$i] = $t2;
			$stat_cycle[$i] = $cycle;
			# printf "stat $i $old_cycle $stat_temp[$i] $ptime\n" if $dbg;
			$i += 1;
		}
		$old_cycle = $cycle;
		$statloc += 36;
	}
	my $max_stat = $i;
	# 1c: Take stat_temp as read_temp when last stat time <= read_time
	
	$j = 0;
	foreach $i (0..($num_cycles - 1))
	{
		# Find next stat time later than read time
#TEKCHANGE: Find next stat time at or later than read time
#		while (($j < $max_stat) and ($read_time[$i] > $stat_time[$j]))
		while (($j < $max_stat) and ($read_time[$i] >= $stat_time[$j]))
		{
			$j += 1;
		}
		# Get temp from that cycle, cycle # from previous cycle
		if ($j < $max_stat)
		{
			$read_temp[$i] = $stat_temp[$j];
			$read_cycle[$i] = $stat_cycle[$j];
		}
		else 
		{
			$read_temp[$i] = $read_temp[$i-1] + 0.2;
			$read_cycle[$i] = $read_cycle[$i-1] + 1;
		}
		# printf "t $i $j " if $dbg;
		# printf "$read_cycle[$i] $read_temp[$i] $read_time[$i]\n" if $dbg;
	}
	# printf "Got temps @read_temp\n" if $dbg;
	
	# 2. Loop through reads and normalize: divide by calibration, find min

	# 2a. Read raw values for each channel, well and temp
	
	$i = 0;
	$j = 0;
	printf "Reading raw data\n" if $dbg;
	while ($i < $num_cycles)
	{
		$off = $itemoff{"read$i"};
		$len = $itemlen{"read$i"};
		$k = $j + $num_reads - 1;
		# printf "$i $j $k $#raw $off\n" if $dbg;
		@raw[$j..$k] = get_plate_read($bin, $off, $len);
		printf "Read cycle $i\n" if $dbg and not ($i % 20);
		my $start = $i*$num_reads;
		my $stop = $start + $num_wells - 1;
		printf READS "read$i,$read_temp[$i],1,".
			"@raw[$start..$stop]\n,,2,".
			"@raw[($start+$num_wells)..($stop+$num_wells)]\n"
			 if $dump_reads;
		$i += 1;
		$j += $num_reads;
	}
	close READS if $dump_reads;

	# 2b. For each well, calibrate at each temp and find min
		
	# cta,b = calibration temps a and b to use in extrapolations
	# ca1 = calibration value for first  temp, first  channel 
	# ca2 = calibration value for first  temp, second channel 
	# cb1 = calibration value for second temp, first  channel 
	# cb2 = calibration value for second temp, second channel 
	# c1 = calibration value for read temp, first channel
	# c2 = calibration value for read temp, second channel
	# v = final value
	local ($rt, $cta, $ctb, $ca1, $ca2, $cb1, $cb2, $c1, $c2, $v1, $v2, $v);
	local ($eca1, $eca2, $ecb1, $ecb2, $ec1, $ec2);
	local ($well_min, $well_sum);
	# Values for cubic interpolation: cal at 3rd & 4th temps, 
	# coefficients a, b, c and d for channels 1 and 2
	# for c = a0nt^3 + a1nt^2 + a2nt + a3n
	local ($cc1, $cc2, $cd1, $cd2, $x, $x2);
	local ($a01, $a11, $a21, $a31, $a02, $a12, $a22, $a32);
	# Values for least squares fit
	local ($xav,$xx,$yav1,$yav2,$yy1,$yy2,$xy1,$xy1,$m1,$m2,$b1,$b2);
	# Values for smoothing:
	local ($half_width, $lo, $running_avg, $running_sum);
	$half_width = int($smooth_n / 2);
	my $ratio;
	my @value;
	my @smoothed;
	$#value = $num_wells * $num_cycles - 1;
	$#smoothed = $#value;
	if ($dump_calibrations)
	{
		open (CSV, ">$csv");
		printf CSV "CALIBRATIONS:\n\nWell,temp,ca1,cb1,ca2,cb2,c1,c2,".
				"raw1,raw2,v\n";
	}	
	foreach $i (0..($num_wells - 1))
	{
		printf "Process well $i\n" if $dbg and not ($i % 12);
		$running_sum = 0;
		$well_min = 100;
		$well_sum = 0;
		if ($interpolation ne "linear")
		{
			$ca1 = $cal[$i];
			$ca2 = $cal[$i + $num_wells];
			$cb1 = $cal[$i + $num_reads];
			$cb2 = $cal[$i + $num_reads + $num_wells];
			$cc1 = $cal[$i + $num_reads*2];
			$cc2 = $cal[$i + $num_reads*2 + $num_wells];
			$cd1 = $cal[$i + $num_reads*3];
			$cd2 = $cal[$i + $num_reads*3 + $num_wells];
			if ($interpolation eq "cubic") {
				# Get all 4 calibration temps for each channel
				$x = $rt/20 - 2;
				$x2 = $x*$x;
				# Find cubic interpolation for chn 1
				$a31 = $cb1;
				$a11 = ($cc1+$ca1)/2 - $cb1;
				$a01 = ($cd1 - 3*$cc1 + 3*$cb1 - $ca1) / 6;
				$a21 = $cc1 - ($a01 + $a11 + $a31);
				# Find cubic interpolation for chn 2
				$a32 = $cb2;
				$a12 = ($cc2+$ca2)/2 - $cb2;
				$a02 = ($cd2 - 3*$cc2 + 3*$cb2 - $ca2) / 6;
				$a22 = $cc2 - ($a02 + $a12 + $a32);
			} 
			elsif ($interpolation eq "least square")
			{
				# Calculate sum of dx^2, dy^2, and dx*dy
				# x=20,40,60,80; mean=50,dx=x-50=30,10,10,30
				# y=caN,...cdN; sum dy^2= sum y^2 - n(y avg)^2
				$xav = 50; 	# = (20+40+60+80)/4;
				$xx = 2000;	# = 900+100+100+900
				$yav1=($ca1 + $cb1 + $cc1 + $cd1)/4;
				$yav2=($ca2 + $cb2 + $cc2 + $cd2)/4;
				$yy1=$ca1*$ca1+$cb1*$cb1+$cc1*$cc1+$cd1*$cd1;
				$yy2=$ca2*$ca2+$cb2*$cb2+$cc2*$cc2+$cd2*$cd2;
				$yy1 -= 4*$yav1*$yav1;
				$yy2 -= 4*$yav2*$yav2;
				$xy1=$ca1*20+$cb1*40+$cc1*60+$cd1*80;
				$xy2=$ca2*20+$cb2*40+$cc2*60+$cd2*80;
				$xy1 -= 200*$yav1;				
				$xy2 -= 200*$yav2;
				$m1 = $xy1 / $xx;
				$m2 = $xy2 / $xx;
				$b1 = $yav1 - $m1 * $xav;
				$b2 = $yav2 - $m2 * $xav;
				printf "$ca2 $cb2 $cc2 $cd2\n\t$m2 $b2\n" if $dbg and not ($i % 20);
			} 
			elsif ($interpolation eq "quadratic")
			{
				$a21 = (3*$ca1+11*$cb1+9*$cc1-3*$cd1)/20;
				$a22 = (3*$ca2+11*$cb2+9*$cc2-3*$cd2)/20;
				$a11 = (-11*$ca1+3*$cb1+7*$cc1+$cd1)/20;
				$a12 = (-11*$ca2+3*$cb2+7*$cc2+$cd2)/20;
				$a01 = (5*$ca1-5*$cb1-5*$cc1+5*$cd1)/20;
				$a02 = (5*$ca2-5*$cb2-5*$cc2+5*$cd2)/20;
				printf "well $i cal $ca2 $cb2 $cc2 $cd2\n".
					"\t $a02 $a12 $a22\n"
					 if $dbg and not ($i%20);
				 
			}
		}
		foreach $j (0..($num_cycles-1))
		{
			$rt = $read_temp[$j];
			if ($interpolation eq "cubic") {
				# Get all 4 calibration temps for each channel
				$x = $rt/20 - 2;
				$x2 = $x*$x;
				# Do cubic interpolation for chn 1
				$c1 = $a01*$x*$x2+$a11*$x2+$a21*$x+$a31;
				# Do cubic interpolation for chn 2
				$c2 = $a02*$x*$x2+$a12*$x2+$a22*$x+$a32;
				printf "well $i cal pts $ca2 $cb2 $cc2 $cd2\n" 
					if $dbg and not ($j or ($i %20));
				printf "$rt $c2\n" 
					if $dbg and not (($j % 25) or ($i%20));
			}
			elsif ($interpolation eq "least square")
			{
				$c1 = $m1 * $rt + $b1;
				$c2 = $m2 * $rt + $b2;
				printf "$rt $c2\n" if $dbg and not (($j%25) or ($i%20));
			} 
			elsif ($interpolation eq "quadratic")
			{
				$x = $rt/20 - 2;
				$c1 = $a01*$x*$x+$a11*$x+$a21;
				$c2 = $a02*$x*$x+$a12*$x+$a22;
				printf "$rt $c2\n" if $dbg and not (($j%25) or ($i%20));
			}
			else
			{
				# First temp index = 0 for t<40, 1 for 40<t< 60, or 2
				$k = int($rt/20) - 1;
				$k = 2 if $k > 2;
				$cta = ($k + 1) * 20; # = 20, 40, 60 or 80
				$ctb = $cta + 20;
				# First temp, first and second channels
				$ca1 = $cal[$i + $num_reads * $k];
				$ca2 = $cal[$i + $num_reads * $k + $num_wells];
#TEKCHANGE: get empty calibrations at first temp, both channels
				$eca1 = $ecal[$i + $num_reads * $k];
				$eca2 = $ecal[$i + $num_reads * $k + $num_wells];
				# Second temp, first and second channels
				$cb1 = $cal[$i + $num_reads * ($k + 1)];
				$cb2 = $cal[$i + $num_reads * ($k + 1) + $num_wells];
#TEKCHANGE: get empty calibrations at second temps, both channels
				$ecb1 = $ecal[$i + $num_reads * ($k + 1)];
				$ecb2 = $ecal[$i + $num_reads * ($k + 1) + $num_wells];

				$ratio = ($rt - $cta) / ($ctb - $cta);

				$c1 = $ca1 + $ratio * ($cb1 - $ca1);
				$c2 = $ca2 + $ratio * ($cb2 - $ca2);
#TEKCHANGE: linear interpolation of empty calibration
				$ec1 = $eca1 + $ratio * ($ecb1 - $eca1);
				$ec2 = $eca2 + $ratio * ($ecb2 - $eca2);
	#			$c1 = $ca1 + ($rt - $cta)*($cb1 - $ca1)/($ctb - $cta);
	#			$c2 = $ca2 + ($rt - $cta)*($cb2 - $ca2)/($ctb - $cta);
			}
			# printf "$i $j max=$#raw ".($i+$j*$num_reads+$num_wells).
			#	" c1,2 =$c1 $c2\n" if $dbg;
			printf CSV "$i,$rt,$ca1,$cb1,$ca2,$cb2,$c1,$c2,"
				 if $dump_calibrations;
#TEKCHANGE: Allow subtraction of blanks from values
#			$v1 = $raw[$i + $j*$num_reads] / $c1;
#			$v2 = $raw[$i + $j*$num_reads + $num_wells] / $c2;
			$v1 = $raw[$i + $j*$num_reads];
			$v2 = $raw[$i + $j*$num_reads + $num_wells];
#FZ change to TEK change: Subtract blanks only if SubBlanks option is TRUE
			if (defined($option_hash{"SubBlanks"}) and
				    $option_hash{"SubBlanks"} eq "TRUE")
			{
				$v1 -= $ec1;
				$v2 -= $ec2;
			}
			$v1 /= $c1;
			$v2 /= $c2;

			$v = $v2;
			if ($sum_channels == 1)
			{
				$v += $v1;
			} elsif ($sum_channels == 2)
			{
				# Try square root of sum of squares
				$v = sqrt($v1*$v1 + $v2*$v2);
			} elsif ($sum_channels == 3)
			{
				# Subtract channel 1 from 2
				$v -= $v1;
			} elsif ($sum_channels == 4)
			{
				# Multiply channels  by cal. from other channel 
				$v = $v1*$c2 + $v2*$c1;
			} elsif ($sum_channels == 5)
			{
				$v = ($v1*$c2/$c1+$v2/$c2)/2
			}
			printf CSV "$raw[$i+$j*$num_reads],".
				"$raw[$i+$j*$num_reads+$num_wells],$v\n"
				if $dump_calibrations;
			$value[$i + $j * $num_wells] = $v;
			if ($smooth_n > 1)
			{
				# Smooth n consecutive points
				# Only smoothing method = average, for now

				# Use end points for points past the ends:
				# i.e. start with n * 1st point
				$running_sum = $v * $smooth_n if ($j == 0);

				# Add new point
				$running_sum += $v;

				# Subtract point before start # of current 
				# window, or subtract 1st point if at start
				$lo = $j - $smooth_n;
				$lo = 0 if $lo < 0;
				# Take off old point, if any
#				printf "j $j lo $lo\n" if $dbg;
				$running_sum -= $value[$i + $lo * $num_wells];
				
				# Assign value to temperature in the middle of 
				# the window
				if ($j >= $half_width)
				{
				    # Store running average value at midpoint
				    $running_avg = $running_sum / $smooth_n;
				    $smoothed[$i+($j-$half_width)*$num_wells] =
						$running_avg;
				    $well_min = $running_avg if 
						$running_avg < $well_min;
				    $well_sum += $running_avg;
				}
#				printf "$j,v=$v,avg=$running_avg,min=$well_min\n" 
#					if $dbg & $j > $num_cycles - 8;
			}
			else 
			{
				$well_min = $v if $v < $well_min;
				$well_sum += $v;
			}
		}
		if ($smooth_n > 1)
		{
			# Fill in the smoothed points at the end of the cycle,
			# using the last value for points past the end
			foreach $j (0..$half_width-1)
			{
				$lo = $num_cycles + $j - $smooth_n;
				# Add new point = last value, $v
				$running_sum += $v;
				# Take off old point
#				printf "j $j lo $lo\n" if $dbg;
				$running_sum -= $value[$i + $lo * $num_wells];
				$running_avg = $running_sum / $smooth_n;
				$well_min = $running_avg if 
						$running_avg < $well_min;
				$well_sum += $running_avg;
				$lo += $half_width + 1;
				$smoothed[$i + $lo * $num_wells] = $running_avg;
#				printf "j $j $lo,avg=$running_avg,min=$well_min\n" if $dbg;
			}
			foreach $j (0 .. $num_cycles)
			{
				$value[$i+$j*$num_wells] = 
					$smoothed[$i+$j*$num_wells];
			}
		}
		my $well_avg = $well_sum / $num_cycles;
		if (defined($option_hash{"SubAvg"}) and
			    $option_hash{"SubAvg"} eq "TRUE")
		{
			# Subtract well average value from each well
			$well_avg -= $well_min;
		} 
		foreach $j (0..($num_cycles-1))
		{
			$v = $value[$i + $j * $num_wells];
			if (defined($option_hash{"SubAvg"}) and
				    $option_hash{"SubAvg"} eq "TRUE")
			{
				# Subtract well average value from each well
				$v -= $well_avg;
			} 
			elsif ((defined($option_hash{"SubGlobalMin"}) and
				$option_hash{"SubGlobalMin"} eq "TRUE" and
				$subtract_baseline) or
				($subtract_baseline == 2))
			{
				# Subtract well minimum value from each well
				$v -= $well_min;
			}
			if ($divide_by_well_avg) 
			{
				$v /= $well_avg;
			}
			# Always round result to 4 decimal places
			$value[$i + $j * $num_wells] = int($v*10000+0.5)/10000;
		}
	}
	printf "Done processing data - starting CSV output\n" if $dbg;
	open (CSV, ">$csv") if not $dump_calibrations;
	printf CSV "$csv_header";
	foreach $i (0..($num_cycles-1))
	{
#TEKCHANGE: Use sprintf to round read temp to 1 decimal place
#		printf CSV "2,$read_cycle[$i], $dye0name,$read_temp[$i],".
		printf CSV "2, $read_cycle[$i], $dye0name, ".
			sprintf("%.1f,", $read_temp[$i]).
			"@value[($i * $num_wells) .. (($i+1)*$num_wells-1)]\n";
	}
	close (CSV);
	printf "Done writing CSV file\n" if $dbg;
	return ($nin) if not $dump_raw;

#############################################################################
###
###-----------------------  DUMP RAW DATA FOR DEVELOPMENT ----------------###
###
#############################################################################
	if ($dump_raw)
	{
		# Dump raw data to _raw.csv
		my $raw = $csv;
		$raw =~ s/.csv$/_raw.csv/g;
		open RAW,">$raw";
		# Print header
		printf RAW "Item,Time,Dye,$well_header\n";
		# Loop through all entries in TOC
		for (@names)
		{
			$len = $itemlen{$_};
			$off = $itemoff{$_};
			$loc = $itemloc{$_};
			printf "$_ ";
			if ($len == $plate_read_len)
			{ 
				my $time = get_time($bin, $off, $len);
				my @vals = get_plate_read($bin, $off, $len);
				# Dump both dye reads to raw data file
				printf RAW "$_,$time,1,@vals[0..95]\n".
					   "$_,$time,2,@vals[96..191]\n";
				# printf "$#vals vals, temp $temp\n";
			}
			elsif ($_ eq "statraw")
			{
				my $statloc = $off + 19;
				my $statend = $statloc + $len;
				my $i = 0;
				my @times;
				my @temps;
				$#times = 1000;
				$#temps = 1000;
				while ($statloc + 36 < $statend)
				{
				    # Dump each statraw record:
				    # time, step, cycle, block & sample temp
				    local($ptime,$step,$cycle,$t1,$t2) =
					unpack("x$statloc fx8V2f2",$bin);
				    $times[$i] = $ptime;
				    $temps[$i] = $t2;
				    $i += 1;
				    printf RAW "$_,,$ptime,$i,$step,".
						"$cycle,$t1,$t2\n";
				    $statloc += 36;
				}
			} 
			else
			{
				printf RAW "$_,,,$off,$len";
				my $oloc = $off;
				# Get data from rest of TOC entry if offset < 0
				$oloc = $loc+length($_)+1 if $off < 0
				;# Dump as both 
				if ($len == 4)
				{
					printf RAW ",".unpack("x$oloc l",$bin);
					printf RAW ",".unpack("x$oloc N",$bin);
					printf RAW ",".unpack("x$oloc f",$bin);
				} 
				elsif ($len == 8)
				{
					printf RAW ",".unpack("x$oloc d",$bin);
				}
				if ($off < 0 or $len < 255)
				{
				# printf RAW ",".unpack("x$oloc a$len",$bin)."\n";
				}
				printf RAW "\n";
			}				
		}
		close RAW;
	}
# test code done to here
#			printf CSV,"\n$nametext,$cycle,$dyename,$temp";
	return ($nin);
}

#	my $cycle = 0;
#	my $num_wells = 96;
#	my $temp;
#	my $dyename = "HEX";

#		if ($nametext =~ /read/ or ($dbg and $namelen = 6161)) 
#		{

sub get_value
{
	local ($itemloc,$itemoff,$value_type,$bin) = @_;
	my $off = $itemoff;
	$off += $itemloc + 255 	if ($off < 0);
	return (unpack("x$off $value_type",$bin));
}

# Get one plate read from tad2 file.
# Arguments: binary data from file, location of plate read in data
# Return: array of values for A1 .. A12 ... H1 .. H12 for dye0 then dye1
sub get_plate_read
{
	my $dbg = 0;
	local ($bin,$offset,$len) = @_;
#	my $timestamp = unpack("x$offset x13f",$bin);
	$read_end = $offset + $len;
	my $resaved = 0;
	if ($len == 6161)
	{
		$offset += 17;
	} elsif ($len == 793 and 
		substr($bin, $offset, 4) eq "\002\002\000\000")
	{
		$offset += 25;
		$resaved = 1;
	} else
	{
		return (0);
	}
	my $well = 0;
	my @values;
	$#values = 96 * 2 - 1;
#	printf "Plate read at time $timestamp offset $offset:\n" if $dbg;
	while ($offset < $read_end)
	{
		# Print text threaded through padding as test
		# printf unpack("x$offset x4 a4",$bin) if $dbg;
		if ($resaved)
		{
			$values[$well] = unpack("x$offset f",$bin);
			$offset += 4;
		} else
		{
			local($n, $sum) = unpack("x$offset lx4d",$bin);
			if ($n > 0) {
				$values[$well] = $sum / $n;
			} else {
				$values[$well] = $sum;
			}
			$offset += 32;
		}
		printf "$well:$sum/$n=$values[$well]\n" if $dbg and ($well%96) == 0;
		$well += 1;
	}
	$#values = $well - 1;
	printf "\nGot values 0 to $#values:@values\n" if $dbg;
	# exit (0) if $dbg;
	return (@values);
}

sub get_time
{
	local ($bin, $offset, $len) = @_;
	$offset += ($len == 793) ? 21 : 13;
	my $t = unpack("x$offset f",$bin);
	return ($t);
}

		
# Read a file of a given name into a string
sub read_binary_file 
{
	local ($filename, $dbg) = @_;
	printf "Trying to read $filename\n" if $dbg;
	if (! (-e $filename))
	{
		printf "ERROR: File $filename does not exist\n";
		return (0,0);
	} elsif (-z $filename)
	{
		printf "ERROR: File $filename is empty\n";
		return (0,0);
	}
	my $bin = "";
	my $buffer_size = 2**16;
	# $buffer_size = 16 if $dbg;
	open (TAD, $filename);
	binmode(TAD);
	# Read first buffer
	my $nin = sysread(TAD, $bin, $buffer_size);
	if (!(defined $nin) or $nin==0)
	{
		printf "ERROR: File $filename returned no data\n";
		close (TAD);
		return ($bin,0);
	}
	print "Reading $filename\n" if $dbg;
	my $offset = $nin;
	printf "First read: buffer = $buffer_size, got $nin bytes\n" if $dbg;
	# Read one buffer at a time to the end
	while ($nin == $buffer_size) 
	{
		$nin = sysread(TAD, $bin, $buffer_size, $offset);
		$offset += $nin;
	}
	close (TAD);
	return ($bin, $offset);
}