#!/usr/bin/perl -w
# Cross_correlation.pl
#    VERSION: 7/22/2004
#    PURPOSE: 
#        Calculates normalized product between two data sets at each point
#          given by:
#        
#                         (point from group A) * (point from group B)
#          correlation =  -------------------------------------------
#                                   normalization factor
#
#          normalization factor = 
#                Square root [(sum of squares of A) * (sum of squares of B)]
#
#        Group B point may be offset by a specified amount
#        For example, an offset of -3 causes a group A point to be compared
#          to a group B point 3 rows up
#        
#        A range of offsets may be specified as +/- offset width
#
#    INPUT FILES:  Tab delimited file
#        1st line: headers for columns of data
#        Subsequent lines: columns of data
#
#    OUTPUT FILES: Tab delimited file
#        1st line: headers for columns of data
#        Subsequent lines: offset, cross_correlation

############## LIBRARIES AND PRAGMAS ################

   use strict;

#################### CONSTANTS ######################

   my $LF = "\n";            # Line feed
   my $tab = "\t";           # Tab
   my $offset_width = 20;    # Range of offsets (+/- value)
   my $group1 = 0;           # Column of data to be considered
   my $group2 = 1;           # Column of data to be considered
 
#################### VARIABLES ######################

   my @points;               # Holds all data
   my @header;               # Holds column headers
   my $number_of_points = 0; # Number of data rows in file

###################### FILES ########################

   my $data_file = 'deltoid.txt';

   my $output_file_name = 'deltoid_correlation.txt';
   open OUTPUT, ">$output_file_name" 
             or die "Can't open $output_file_name: $!\n";

################### MAIN PROGRAM ####################
### Read data from file into array @points
### Correlate columns of data specified by group1 and group2

   Read_data_set($data_file);
   print "Number of items read: ", $number_of_points, $LF;
 
   Correlate_data();

#################### SUBROUTINES ####################

#### READ_DATA_SET (path to input file)
#    Reads in the data set from given input file
#    Stores columns of data in array, @points
#    @points is 2-dimensional array
#        First dimension is group (column of data)
#        Second dimension is time (row of data)
  sub Read_data_set
  {
     my ($path) = @_;  	# sets file path location with variable 
                        # passed to Read_data_set 
     open DATA_FILE, "<$path" 
 	             or die "Can't open $path: $!\n";
     my $line = <DATA_FILE>; 
     @header = (split (" ", $line));
     $line = <DATA_FILE>;
		 
     while ($line) {   # while there is data left in the input file
        chomp $line;   # strip off carriage return from input line
        push @points, [split (" ", $line)];
        $number_of_points = $number_of_points + 1;
        $line = <DATA_FILE>;
     }
  }

#### CORRELATE_DATA
#    Calculates correlations and prints values file
# 1. Prints header
# 2. Calculates group1 part of normalization factor 
# 3. Loops through all offsets
#    a. Prints offset to file
#    b. Loops through data
#       i. Calculates crossproduct
#      ii. Calculates group2 part of normalization factor
#    c. Calculates normalization factor
#    d. Prints correlation to file

  sub Correlate_data {
     my $first_point = $offset_width;
     my $last_point = $number_of_points - $offset_width - 1;
     my $sum_of_squares_group1 = 0;
     my $sum_of_squares_group2;
     my $sum_of_products12;
     my $normalization_factor12;

     print OUTPUT "", $tab, $header[$group1], " vs ", $header[$group2], $LF;

     foreach my $i ($first_point .. $last_point) {
        $sum_of_squares_group1 +=   $points[$i][$group1]
                                  * $points[$i][$group1];
     }

     foreach my $offset (-$offset_width .. $offset_width) {
        $sum_of_products12 = 0;
        $sum_of_squares_group2 = 0;

        print OUTPUT $offset, $tab;
 
        foreach my $i ($first_point .. $last_point) {
           $sum_of_products12 +=   $points[$i][$group1] 
                                 * $points[$i + $offset][$group2];
           $sum_of_squares_group2 +=   $points[$i + $offset][$group2]
                                     * $points[$i + $offset][$group2];
        }

        $normalization_factor12 = sqrt(  $sum_of_squares_group1
                                       * $sum_of_squares_group2);  
 
        printf OUTPUT "%8.5f", $sum_of_products12 / $normalization_factor12;
        print OUTPUT $LF;
     }
  }
