#####################################################################
#
# (c) 2001-2003 Bruz Marzolf
#
# quantitation.pl - functions to load quantitation files into a
# standard format and checking grid alignment.
#
# Quantition file loading - works with QuantArray, Dapple
# and AnalyzerDG output files. Checks first line in quantitation
# file to determine spotfinder used, should ID QuantArray and
# Analyzer correctly, assumes anything else is dapple
#
# Grid alignment is based on correlation between slides,
# determines if there is a better alignment by shifting the grid
# up to 1/8 of total # of spots in each dimension x and y 
#
# Last Update November 5, 2003
#
#####################################################################

use strict;
use Benchmark;
use PDL;
use PDL::Fit::Polynomial;
use PDL::Dbg;

my $DEBUG = 0;

#METAROW     - number of rows of grids
#METACOL     - number of rows columns of grids
#MAXROW      - rows of spots per grid
#MAXCOL      - columns of spots per grid
#MAXSLIDEROW - rows of spots if entire array is considered one big grid
#MAXSLIDEROW - columns of spots if entire array is considered one big grid

my $MAXMETAROW;
my $MAXMETACOL;
my $MAXROW;
my $MAXCOL;
my $MAXSLIDEROW;
my $MAXSLIDECOL;

#default is not to perform alignment test
my $alignChoice = "none";   

#immediately flush ouput
$|=1;                       

#Initialize all size measures to 0
$MAXSLIDEROW = $MAXSLIDECOL = $MAXMETAROW = $MAXMETACOL = $MAXROW = $MAXCOL = 0;

#####################################################################

# read quantitation file
sub readSpotFile {
    my $spotfile = shift @_;     #quantitation file name
    my $spots = shift @_;        #spot information for a slide
   
    open(INFILE, "$spotfile");   
    my $header = <INFILE>;       #toss header
    $header =~ s/\c@//g;         #remove unicode for QA

    #split up header contents
    my @headerItems = split /\t|,| /, $header;
    
    #Look for identification:
    # All quantarray output files start with 'User Name'
    if( $header =~ /User Name/ ) {       
	readQuantArrayFile( \*INFILE, $spots );
    }
    # AnalyzerDG first line  must contain 'Spot ...'
    elsif( $header =~ /Spot/ ) {	 
	readAnalyzerDGFile( \*INFILE, $spots, $spotfile, $header );
    }
    # If we don't know what it is and the first line has two #'s,
    # it's probably dapple format...
    elsif( $headerItems[0] =~ /\d/ && $headerItems[1] =~ /\d/ ) {	
	readDappleFile( \*INFILE, $spots );   
    }
#####################################################################
# Place check for other spotfinder here
#
# elsif( some characteristic of particular quantitation header) {
#       readSomeKindOfFile( \*INFILE, $spots );
# }
#####################################################################
    else {
	die "Unrecognized file format encountered - $spotfile\n"; 	
    }
    
    close(INFILE);
}

#####################################################################

sub readQuantArrayFile {
    my $infile = shift @_;
    my $spots = shift @_;
    my (%cols, @entry, $n, $header);
    my ($metarow, $metacol, $row, $col);

    $spots->{spotfinder} = "QuantArray";

    #first half of quantarray file is composed of a header and a bunch of
    #ratio data we don't want, so we skip all lines until 'Begin Data'
    while(<INFILE>) {
	s/\c@//g;                    
	last if $_=~/Begin Data/;	 
    }
    $header = <INFILE>;	                 #grab header line,
    $header =~ s/\c@//g;                 #remove unicode

    #find number of channels in header by looking at the highest
    #ch# header column
    if( $header =~ /ch4/ )    { $spots->{chan} = 4; }
    elsif( $header =~ /ch3/ ) { $spots->{chan} = 3; }
    elsif( $header =~ /ch2/ ) { $spots->{chan} = 2; }
    elsif( $header =~ /ch1/ ) { $spots->{chan} = 1; }

    #assign which column intensity and background are in for each
    #channel in the QuantArray output file
    for my $i ( 1..$spots->{chan} ) {
	push @{$cols{int}}, (4 + 4 * $i);
	push @{$cols{bck}}, (5 + 4 * $i);
    }

    #Places to look for spot coordinates, these are constant for quantarray
    $cols{metarow} = 1;             
    $cols{metacol} = 2;
    $cols{row} = 3;
    $cols{col} = 4;

    $cols{offset} = 1;  #Subtract 1 from each coordinate to get 0->n numbering
    while(<INFILE>) {
	s/\c@//g;
	@entry = split /\t/;        #split on tab, comma, space

	#For each coordinate, Subtract 1 to get 0->n numbering
	$metarow = $entry[$cols{metarow}] - 1;     
	$metacol = $entry[$cols{metacol}] - 1;
	$row = $entry[$cols{row}] - 1;
	$col = $entry[$cols{col}] - 1;

	#Grab intensity and background for each channel
	for $n ( 0..$#{$cols{int}} ) {
	    $spots->{int}[$n][$metarow][$metacol][$row][$col] = 
		$entry[$cols{int}[$n]];
	    $spots->{bck}[$n][$metarow][$metacol][$row][$col] = 
		$entry[$cols{bck}[$n]];
	}	
    }
}

#####################################################################

sub readAnalyzerDGFile {
    my $infh = shift @_;
    my $spots = shift @_;
    my $spotfile = shift @_;
    my $header = shift @_;
    my (%cols, @entry, $n, @dyes);
    my ($metarow, $metacol, $row, $col);

    $spots->{spotfinder} = "AnalyzerDG";
    
    if( $header =~ /Spot Mean/ ) { $spots->{intType} = "Mean"; }
    elsif ( $header =~ /Spot Median/ ) { $spots->{intType} = "Median"; }
    else { die "Neither Spot mean nor median intensity found!\n"; }

    #Look for all channel names, and count them
    my $matches = 0;
    while( $header =~ /Spot $spots->{intType}.*?\((.*?)\)/g ) {   
	push @dyes, $1;
	$matches++;
    }
    $spots->{chan} = $matches;

    #Sub returns a hash of column #'s for each piece of info needed
    %cols = getDGcols($header, $spots->{intType}, @dyes);
    
    # Find maximum # of meta rows and cols if appropriate columns
    # exist in ADG file.  Assumes last line of file will have highest
    # number possible for super rows and cols.
    my $maxMetaRow = 0;
    my $maxMetaCol = 0;
    my $lastline; 
    if( $cols{superrow} > 0 || $cols{supercol} > 0 ) {
	while(<$infh>) { $lastline = $_; }  
	@entry = split /\t|,| /, $lastline;    #split on tab, comma, space
	if( $entry[$cols{metarow}] > $maxMetaRow ) {
	    $maxMetaRow = $entry[$cols{metarow}];
	}
	if( $entry[$cols{metacol}] > $maxMetaCol ) {
	    $maxMetaCol = $entry[$cols{metacol}]; 
	}
    }    
    close($infh);

    open(INFILE, "$spotfile") or die "Couldn't open $spotfile to read\n";
    $header = <INFILE>;
    my $superrow = 0;
    my $supercol = 0;
    while(<INFILE>) {
	chomp( $_ );
	@entry = split /\t|,| /;        #split on tab, comma, space
	if( $cols{superrow} > 0 || $cols{supercol} > 0 ) {
	    $superrow = $entry[$cols{superrow}]-1;
	    $supercol = $entry[$cols{supercol}]-1;
	}
        #change layout info to acct. for supergrids, also subtract
	#1 from each coordinate to get 0->n numbering
	$metarow = $entry[$cols{metarow}] - 1 + $superrow*$maxMetaRow;  
	$metacol = $entry[$cols{metacol}] - 1 + $supercol*$maxMetaCol;  
	$row = $entry[$cols{row}] - 1;
	$col = $entry[$cols{col}] - 1;
	for $n ( 0..$#{$cols{int}} ) {
	    $spots->{int}[$n][$metarow][$metacol][$row][$col] = 
		   $entry[$cols{int}[$n]];
	    $spots->{bck}[$n][$metarow][$metacol][$row][$col] = 
		   $entry[$cols{bck}[$n]];
	    $spots->{conf}[$n][$metarow][$metacol][$row][$col] = 
		   $entry[$cols{conf}[$n]];
	}       	
    }
}

#Determine where appropriate columns are in MolecularWare files
sub getDGcols {
    my $header = shift @_;
    my $intType = shift @_;
    my $dye;
    print "dyes: $_[0] $_[1] $_[2] $_[3]\n" if $DEBUG==2;
    my %cols;
    
    my @header = split /\t|,/, $header;
    for my $i (0.. $#header) {
	print $header[$i], "\n" if $DEBUG > 2;
	if( $header[$i] =~ /Multi-Set Row/ ) {
	    $cols{superrow} = $i;
	}
	elsif( $header[$i] =~ /Multi-Set Col/ ) {
	    $cols{supercol} = $i;
	}
	elsif( $header[$i] =~ /Set Row/ ) {
	    $cols{metarow} = $i;
	}
	elsif( $header[$i] =~ /Set Column/ ) {
	    $cols{metacol} = $i;
	}
	elsif( $header[$i] =~ /Row/ ) {
	    $cols{row} = $i;
	}
	elsif( $header[$i] =~ /Column/ ) {
	    $cols{col} = $i;
	}
	foreach $dye (@_) {
	    if( $header[$i] =~ /Spot $intType Intensity \($dye\)/ ){ 
		push @{$cols{int}},$i; 
	    }
	    elsif( $header[$i] =~ /Background $intType Intensity \($dye\)/ ){
		push @{$cols{bck}},$i; 
	    }
	    elsif( $header[$i] =~ /Spot Confidence \($dye\)/ ){ 
		push @{$cols{conf}},$i; 
	    }
	}
    }
    
    return %cols;
}

#####################################################################

sub readDappleFile {
    my $infile = shift @_;
    my $spots = shift @_;
    my (%cols, @entry, $n);
    my ($metarow, $metacol, $row, $col);
    
    $spots->{spotfinder} = "Dapple";
    $spots->{chan} = 2;
    $cols{int} = [5,9];
    $cols{bck} = [6,10];
    $cols{metarow} = 1;             #these are constant for dapple
    $cols{metacol} = 0;
    $cols{row} = 3;
    $cols{col} = 2;
    while(<INFILE>) {
	@entry = split /\t|,| /;        #split on space
	for $n ( 0..$#{$cols{int}} ) {
	    $metarow = $entry[$cols{metarow}];     
	    $metacol = $entry[$cols{metacol}];
	    $row = $entry[$cols{row}];
	    $col = $entry[$cols{col}];
	    $spots->{int}[$n][$metarow][$metacol][$row][$col] = 
		$entry[$cols{int}[$n]];
	    $spots->{bck}[$n][$metarow][$metacol][$row][$col] = 
		$entry[$cols{bck}[$n]];	    
	}
    }
}

#####################################################################

#####################################################################
# Place code to read other spotfinder outputs here
#
# sub readSomeOtherFile {
#       ...
# }
#####################################################################

#####################################################################

# find dimensions for array in $spotfile
sub findMaxRowCol {
    my $quant = shift @_;
    my $numSlides = shift @_;
    my ($MetaRow ,$MetaCol, $Row, $Col);
    my $spots;
    
    #Set maximum dimensions as dimensions of first quantitation
    $spots = pdl $quant->[0]{int}[0];
    ($MAXCOL, $MAXROW, $MAXMETACOL, $MAXMETAROW) = $spots->dims;
    
    $MAXMETAROW -= 1;  
    $MAXMETACOL -= 1;
    $MAXROW -= 1;
    $MAXCOL -= 1;

    #Check dimensions on first quantitation versus all others,
    #complain and return -1 if they don't all match
    for my $n (1..$numSlides) {
	$spots = pdl $quant->[$n]{int};
	($Col, $Row, $MetaCol, $MetaRow) = $spots->dims;
	if($MAXMETACOL != $MetaCol-1 || $MAXMETAROW != $MetaRow-1 || 
	   $MAXCOL != $Col-1 || $MAXROW != $Row-1) {
	    print "Array Sizes aren't matching between slides!!!\n" .
		  "ref: $MAXCOL, $MAXROW, $MAXMETACOL, $MAXMETAROW vs " .
	          "test: $Col, $Row, $MetaCol, $MetaRow";
	    return -1;    #sizes not matching
	}
    }
    
    return 1;             #1=no discrepancies in array size
}

#####################################################################

#Performs an alignment check between first slide in a bunch and the
#rest of the slides, displays output as an html table
sub checkAlign {
    #Remnants of benchmarking to optimize this code,
    #which is still not lightning fast
    my ($t0, $t1) if $DEBUG == 1;
    $t0 = new Benchmark if $DEBUG == 1;
    
    my ($quant) = shift @_;      #quantitation data
    my $numSlides = shift @_;    #number of slides in quant
    my ($quantFiles) = shift @_; #quantation file names
    my $alignChoice = shift @_;  #choice of alignment:text or html

    #This function can't work without 2 slides
    if( $numSlides < 1 ) { 
	print "# Need at least two slides for alignment!\n"; return -1; 
    }

    #These are used to build an HTML table
    my @singlerow;
    my @rrows;
    my $comment;

    #current correlation coefficient
    my $cor;

    #Can only look for fits with 1/4 of a grid shift (up to 4 rows
    #if there are 16 rows in each grid
    my $max_x_shift = int( $MAXCOL / 4 );       
    my $max_y_shift = int( $MAXROW / 4 );       
    
    #Not sure if this is even necessary now
    my ($refpdl, $testpdl) = 0;

    #Sum all channels on reference array
    for my $c ( 0..($quant->[0]{chan}-1) ) {
				$refpdl += pdl $quant->[0]{int}[$c];
    }

    #Optimization 
    $t1 = new Benchmark if $DEBUG == 1;
    showTimeDiff( $t0, $t1, "sum channels" ) if $DEBUG == 1;

    #Iterate through all slides
    for my $n ( 1..$numSlides ) {

	#Sum all channels on test array
	$testpdl = 0;
	for my $c ( 0..($quant->[$n]{chan}-1) ) {
	    $testpdl += pdl $quant->[$n]{int}[$c];
	}

  # temporarily reshape PDLs to allow correlation assessment
	my @refdims = $refpdl->dims;
	reshape $refpdl, $refdims[0]*$refdims[1]*$refdims[2]*$refdims[3];
	my @testdims = $testpdl->dims;
	reshape $testpdl, $testdims[0]*$testdims[1]*$testdims[2]*$testdims[3];


	#Print out name of slides being compared,
	#either in html or text
	if( $alignChoice eq "html" ) {
	    print "<p><b>$quantFiles->[0] vs. $quantFiles->[$n]</b></p>";
			if( correl_PDL( $refpdl, $testpdl ) < 0.7 ) {
					print "<p><b><font color=\"#ff0000\" >WARNING: Overall correlation under 0.7, results may not be accurate</font></b></p>";
			}
	    for my $metacol (0..$MAXMETACOL) { 
					$singlerow[$metacol] = $metacol + 1; 
	    }
	    unshift(@singlerow, "");
	    push( @rrows, td( \@singlerow ) );
	}
	elsif( $alignChoice eq "text" ) {
	    print "\n$quantFiles->[0] vs. $quantFiles->[$n]\n\t";
	    for my $metacol (0..$MAXMETACOL) { print $metacol + 1, "\t"; }
	    print "\n";
	}
	
	# return PDLs to original shape
	reshape $refpdl, $refdims[0],$refdims[1],$refdims[2],$refdims[3];
	reshape $testpdl, $testdims[0],$testdims[1],$testdims[2],$testdims[3];

	my $position = null;
	my @positionName;

	#PDL slices for testing different amounts of shift
	my ( $testShift, $refShift );
	
	#The following are for the correlation routine
	#done inline here, but also located in correl_PDL
	#with more explanation
	my $points;
	my ( $sumsqrx, $sumxsqr, $ssx );
	my ( $sumsqry, $sumysqr, $ssy );
	my ( $sumprodxy, $sumxyprod, $ssxy );

	#Optimization
	$t0 = new Benchmark if $DEBUG == 1;

	#Iterate through all allowed shifts in y
	for my $y_shift (-$max_y_shift..$max_y_shift) {
	    #Iterate through all allowed shifts in x
	    for my $x_shift (-$max_x_shift..$max_x_shift) {
		#Logic for choosing frame of comparison on reference.
		#Basically, keeps slice within boundaries of pdl,
		#while changing what data is framed
		my $y_start = 0;
		my $y_stop = $MAXROW;
		my $x_start = 0;
		my $x_stop = $MAXCOL;
		if( $y_shift < 0 ) {
		    $y_start = -$y_shift;
		}
		elsif($y_shift > 0) { 
		    $y_stop = $MAXROW-$y_shift; 
		}
		if( $x_shift < 0 ) {
		    $x_start = -$x_shift; 			
		}
		elsif($x_shift > 0) { 
		    $x_stop = $MAXCOL-$x_shift;
		}
		
		#actually make the slice for reference array
		$refShift = 
		    $refpdl->slice("$x_start:$x_stop,$y_start:$y_stop,:,:");
		
		#logic for test
		$y_start = 0;
		$y_stop = $MAXROW;
		$x_start = 0;
		$x_stop = $MAXCOL;
		if( $y_shift < 0 ) {
		    $y_stop = $MAXROW+$y_shift; 
		}
		elsif($y_shift > 0) { 
		    $y_start = $y_shift;
		}
		if( $x_shift < 0 ) {
		    $x_stop = $MAXCOL+$x_shift; 
		}
		elsif($x_shift > 0) { 
		    $x_start = $x_shift;
		}
		
		#make the slice for test array
		$testShift = 
		    $testpdl->slice("$x_start:$x_stop,$y_start:$y_stop,:,:");

		#Same as correl_PDL, but put inline here to increase
		#speed by not having overhead of passing args. See
		#that function for comments
		$points = ($x_stop - $x_start + 1)*($y_stop - $y_start + 1);
		$sumsqrx = sumover( sumover( $refShift * $refShift ) );
		$sumxsqr = sumover( sumover( $refShift ) ) ** 2;
		$ssx = $sumsqrx - $sumxsqr / $points;

		$sumsqry = sumover( sumover( $testShift * $testShift ) );
		$sumysqr = sumover( sumover( $testShift ) ) ** 2;
		$ssy = $sumsqry - $sumysqr / $points;

		$sumprodxy = sumover( sumover( $refShift * $testShift ) );
		$sumxyprod = sumover( sumover( $refShift ) ) * 
		    sumover( sumover( $testShift ) );
		$ssxy = $sumprodxy - $sumxyprod / $points;

	        #resulting correlation coefficient
		$cor = $ssxy / sqrt( $ssx * $ssy );
		
		#put correlations into a pdl of all correlations
		#in order of posiition for this array<->array comparison 
		if( $position->isempty ) {
		    $position = $cor->cat;
		}
		else {
		    $position = cat( $position->dog, $cor );
		}
		#record the names of the positions that corresponds
		#to the pdl $position
		push @positionName, fit_direction( $y_shift, $x_shift );
	    }
	}

	#Optimization
	$t1 = new Benchmark if $DEBUG == 1;
	showTimeDiff( $t0, $t1, "correlation code" ) if $DEBUG == 1;
	
	#Reorder positional correlations to index sort each grid
	$position = $position->reorder( 2,1,0 );
	my $positionSort = $position->qsorti;
	#index in positionSort for each grid's highest correlation 
	my $high_position = ( $max_x_shift * 2 + 1 ) * 
	                    ( $max_y_shift * 2 + 1 ) - 1;
	#calculate highest - lowest correlation for each grid
	my $cor_diffs = 
	    $position->index( $positionSort->slice("($high_position),:,:") ) -
	    $position->index( $positionSort->slice("(0),:,:") );

	#find average of these correlation differences (highest - lowest)
	my $avg_diff;
	( $avg_diff ) = $cor_diffs->stats;

	#Set the cutoff for a good vs. bad grid at half of the avgerage
	#difference.
	#***This method of doing things will fail if too many grids are 'bad'
	#and don't have any decent correlation
	my $cutoff = $avg_diff / 2;

	#AnalyzerDG numbering system, used for display of output
	my $dgNum = 1;                     
	     
	#direction recommended to move grid
	my $direction;
	
	#HTML table option
	if( $alignChoice eq "html" ) {
	    #clear out table contents
	    @rrows = ();
	
	    for my $metarow (0..$MAXMETAROW) {
		#clear out individual row
		@singlerow = ();
		for my $metacol (0..$MAXMETACOL) {
		    #get direction name based on index sort of correlations
		    $direction = $positionName[ $positionSort->at( 
				 $high_position, $metarow, $metacol ) ];
		    
		    $comment = "<font color=\"#000000\" ><b>$dgNum:</b>",
		               "</font>";
		    #if correlation is below cutoff, it's a 'bad grid'
		    if( $cor_diffs->at( $metarow, $metacol ) < $cutoff ) {
			$comment .= "<font color=\"#ff0000\" >bad grid",
			            "</font>";
		    }	      
		    #if it's not a bad grid, it's either good or needs to be
		    #moved in order to make it good
		    else {
			if( $direction eq "good" ){
			    $comment .= "<font color=\"#007700\" >good",
				        "</font>";
			}
			else {
			    $comment .= "<font color=\"#0000ff\" >$direction",
				        "</font>";
			}
		    }
		    $singlerow[$metacol] = $comment;

		    #up analyzerDG's numbering scheme
		    $dgNum += 1;
		}

		#add a row number at the beginning of each row
		unshift( @singlerow, $metarow + 1 );

		#push the single row onto the row of rows (table)
		push( @rrows,td( \@singlerow ) );
	    }	
	    
	    #print the result as an HTML table
	    print table({-border=>undef,-width=>'25%',-bgcolor=>'#FFFFFF'}, 
			Tr(\@rrows) ); 
	}
	
	#Text table option
	elsif( $alignChoice eq "text" ) {
	    for my $metarow (0..$MAXMETAROW) {

		#add a row number at the beginning of each row
		print "$metarow\t";

		for my $metacol (0..$MAXMETACOL) {
		    #get direction name based on index sort of correlations
		    $direction = $positionName[ $positionSort->at( 
				 $high_position, $metarow, $metacol ) ];
		    print $dgNum;

		    #if correlation is below cutoff, it's a 'bad grid'
		    if( $cor_diffs->at( $metarow, $metacol ) < $cutoff ) {
			print "$dgNum:bad grid\t";
		    }	

		    #if it's not a bad grid, it's either good or needs to be
		    #moved in order to make it good
		    else {
			if( $direction eq "good" ){
			    print "good\t";
			}
			else {
			    print "$direction\t";
			}
		    }
		    #up analyzerDG's numbering scheme
		    $dgNum += 1;
		}
		print "\n";
	    }	
	    print "\n";
	}

    }
}

#This function is used exclusively by checkAlign,
#generates a text description of how to move grid based on best fit values
sub fit_direction {
    #Parameters are the number of rows and columns that need to be shifted
    my $y_shift = shift @_;
    my $x_shift = shift @_;
    
    #direction starts off blank
    my $direction = "";
    
    if( $y_shift < 0 ) {           #negative y_shift means it must move up
	$y_shift = -$y_shift; 	   #make positive
	$direction .= "$y_shift"."up";  
    }
    elsif( $y_shift > 0 ) {        #negative y_shift means it must move down
	$direction .= "$y_shift"."down"; 
    }		

    if( $x_shift < 0 ) {           #negative x_shift means it must move left
	$x_shift = -$x_shift; 	   #make positive
	$direction .= "$x_shift"."left"; 
    }
    elsif( $x_shift > 0 ) {        #negative x_shift means it must move right
	$direction .= "$x_shift"."right";
    }

    #if we can't find anything bad to say about the grid, call it good.
    #this may be 
    if($direction eq "") { $direction = "good"; }

    return $direction;
}

#####################################################################

# Show time of benchmark
sub showTimeDiff {
    #get two times to compare, and a name for the function we're timing
    my $timeA = shift @_;
    my $timeB = shift @_;
    my $funcName = shift @_;
    
    #find and print the amount of time between points A and B if we're
    #in 'debug' mode
    my $timeDiff = timediff($timeB, $timeA);
    print "Time to $funcName: ",timestr($timeDiff),"<br>" if $DEBUG == 1;
}

#####################################################################

#Takes 2 lists and finds the correlation coefficient (R) between the
#two, BUT the arrays must be the same size and have numbers in all
#their positions (no empty slots)
sub correl_PDL {
    #Two lists to be compared
    my $x = shift @_;
    my $y = shift @_;

    #Turn them into PDLs
    my @x_dims = $x->dims;
    my @y_dims = $y->dims;

    #Number of points being compared
    my $points;
    #Sum of squares of x, sum of x's squared, sum for squared deviates for x
    my ($sumsqrx, $sumxsqr, $ssx);
    #Sum of squares of y, sum of y's squared, sum for squared deviates for y
    my ($sumsqry, $sumysqr, $ssy);
    #Sum of products of x's and y's, product of the sum of x's and sum of y's,
    #sum of co-deviates for paired x's and y's
    my ($sumprodxy, $sumxyprod, $ssxy);

    #Freak out if the lists are not equal size
    if( $x_dims[0] != $y_dims[0] ) {
	die "Attempt to correlate two unequal pdls...\n";
    }

    #Find number of points, and return 0 if there aren't any
    $points = $x_dims[0];
    return 0 if $points == 0;

    #Some stats
    $sumsqrx = sum( $x * $x );
    $sumxsqr = sum( $x ) ** 2;
    $ssx = $sumsqrx - $sumxsqr / $points;
    
    $sumsqry = sum( $y * $y );
    $sumysqr = sum( $y ) ** 2;
    $ssy = $sumsqry - $sumysqr / $points;

    $sumprodxy = sum( $x * $y );
    $sumxyprod = sum( $x ) * sum( $y );
    $ssxy = $sumprodxy - $sumxyprod / $points;
    
    #return 0 if the result would require taking the
    #square root of a negative number
    if( $ssx * $ssy > 0 ) {
	return $ssxy / sqrt( $ssx * $ssy );
    }
    else {
	return 0;
    }
}

#####################################################################

1;
