#!/usr/bin/env perl
use strict;
use warnings;


# This program performs a block jacknife on genomic data to obtain confidence intervals

my $inputfile = $ARGV[0];

unless (open DATAINPUT, $inputfile) {
	print "Can not find the input file.\n";
	exit;
}

my @sumsites;
my $fDM_numsites;
my @fDM_numsites;
my @fdm_values;
my @weighted_fdm_values;
my $counter=0;
my $y;
my $x;
my @temp; 
my $weighted_fdm_values;
my $non_weighted_fdm_value;
my $ABBA=0;
my $BABA=0;
my $BBAA=0;

while ( my $line = <DATAINPUT>) {
	@temp=split('\t',$line);
	if($temp[0] ne 'chromosome'){
		# if fDM is not numeric, it shouldn't contribute to the weighted average
		if($temp[9] !~ /NAN/){
			# load the number of ABBA and BABA sites for each window
			$sumsites[$counter]=$temp[3]+$temp[4];
			# load the number of fDM sites for each window
			$fDM_numsites[$counter]=$temp[10];
			# add them up to eventually get the average
			$fDM_numsites+=$temp[10];
			# also load the fDM stat for each window
			$fdm_values[$counter]=$temp[9];
			# also calculate the non-weighted value
			$non_weighted_fdm_value+=$temp[9];
			# keep track of how many windows are loaded
			$counter+=1;
			$ABBA+=$temp[3];
			$BABA+=$temp[4];
			$BBAA+=$temp[5];
		}
	}	
}		

# Calculate real stat
# make $fDM_numsites the average per site
$fDM_numsites=$fDM_numsites/($counter);
my $weighted_average=0;
# calculate the weighted average of fDM
for ($y = 0 ; $y <= $#fdm_values; $y++ ) {
	$weighted_average+=($fDM_numsites[$y]*$fdm_values[$y]/$fDM_numsites);
}	
print "The number of ABBA sites is ",$ABBA,"\n";
print "The number of BABA sites is ",$BABA,"\n";
print "The number of BBAA sites is ",$BBAA,"\n";
print "The average number of fDM sites per window is ",$fDM_numsites,"\n";
#print "The non-weighted average is ",sprintf("%.6f",$non_weighted_fdm_value/($counter)),"\n";
$weighted_average=$weighted_average/($counter);
print "The weighted average is ",sprintf("%.6f",$weighted_average),"\n";

# now calculate the standard error.
my @jack_sumsites;
my $jack_averagesites;
my @jack_fdm_values;
my $jack_weighted_average=0;
my @jack_weighted_fdm_values;
my @jackarray;
my $counter2=0;

for ($y = 0 ; $y < $counter; $y++ ) {
	for ($x = 0 ; $x < $counter; $x++ ) {
		# leave out one row for each jackknfe replicates
		if($y != $x){
			# load the number of fDM sites for each window
			$jack_sumsites[$counter2]=$fDM_numsites[$x];
			# add them up to eventually get the average
			$jack_averagesites+=$fDM_numsites[$x];
			# also load the fDM stat for each window
			$jack_fdm_values[$counter2]=$fdm_values[$x];
			# keep track of the number of windows
			$counter2+=1;			
		}
	}
	# make the $jack_averagesites an average
	$jack_averagesites=$jack_averagesites/($counter2);
	# calculate the weighted average
	for ($x = 0 ; $x < $counter2; $x++ ) {
		$jack_weighted_average+=($jack_sumsites[$x]*$jack_fdm_values[$x]/$jack_averagesites);
	}
	push(@jackarray,($jack_weighted_average/$counter2));
	# reset variables
	$jack_weighted_average=0;
	$jack_averagesites=0;
	@jack_fdm_values=();
	$counter2=0;
}

# now calculate the variance of the jackknife replicates

# first we need the mean
my $jack_mean=0;
for ($x = 0 ; $x < $counter; $x++ ) {
	$jack_mean+=$jackarray[$x];
}
$jack_mean=$jack_mean/($counter);

my $jack_var=0;
for ($x = 0 ; $x < $counter; $x++ ) {
	$jack_var+=($jack_mean-$jackarray[$x])**2;
}

# for the sample variance, divide by (n-1)
print "jackvar ",sprintf("%.9f",$jack_var/($counter-1)),"\n";
my $sterr = sqrt($counter*($jack_var/($counter-1)));
print "The standard error of the weighted fDM is ",sprintf("%.5f",$sterr),"\n";
print "The 95\%CI of the weighted fDM is ",
sprintf("%.6f",($weighted_average-1.96*$sterr))," - ",sprintf("%.6f",($weighted_average+1.96*$sterr)),"\n";

close DATAINPUT;



