#!/usr/common/usg/languages/perl/5.16.0/bin/perl -w
#
# insertSizeCorrection2.pl by Jarrod Chapman <jchapman@lbl.gov> Tue Jun 16 19:56:11 PDT 2009
# Copyright 2009 Jarrod Chapman. All rights reserved.
#

use Getopt::Std;

my %opts = ();
my $validLine = getopts('m:i:G', \%opts);
my @required = ("m","i");
my $nRequired = 0;
map {$nRequired += exists($opts{$_})} @required;
$validLine &= ($nRequired == @required);
if ($validLine != 1) {
    print "Usage: ./insertSizeCorrection2.pl <-m merSize> <-i insertSizeDistributionFile> <<-G(aussianApproximation?)>>\n";
    exit;
}

my $merSize = $opts{"m"};
my $insertSizeFile = $opts{"i"};

my $gaussian = 0;
if (exists($opts{"G"})) {
    $gaussian = 1;
}

my $minSize = undef;
my $maxSize = undef;
my %sizeDist = ();

my $totalCounts = 0;

open (I,$insertSizeFile) || die "Couldn't open $insertSizeFile\n";
while (my $line = <I>) {
    chomp $line;
    if ($line =~ /^#/) {
	next;
    }
    my ($size,$freq) = $line =~ /(\d+)\s+(\d+)/;
    unless (defined($minSize)) {
	$minSize = $size;
    } else {
	if ($size < $minSize) {
	    $minSize = $size;
	}
    }
    unless (defined($maxSize)) {
	$maxSize = $size;
    } else {
	if ($size > $maxSize) {
	    $maxSize = $size;
	}
    }

    $sizeDist{$size} = $freq;
    $totalCounts += $freq;

}
close I;

my $meanSize = 0;
my $meanSquareSize = 0;

my @observedSizes = sort {$a <=> $b} (keys(%sizeDist));

foreach my $size (@observedSizes) {
    $sizeDist{$size} /= $totalCounts;
    $meanSize += $size*$sizeDist{$size};
    $meanSquareSize += $size*$size*$sizeDist{$size};
}

my $stdDev = sqrt($meanSquareSize-$meanSize*$meanSize);

print STDERR "#Mean end separation: $meanSize +/- $stdDev\n";

if ($gaussian) {
    for (my $s = $minSize; $s <= $maxSize; $s++) {
	$sizeDist{$s} = exp(-(($s-$meanSize)*($s-$meanSize))/(2*($stdDev*$stdDev)));
    }
    @observedSizes = sort {$a <=> $b} (keys(%sizeDist));
}

my $twoM = 2*$merSize;

for (my $g = -$merSize+1; $g < $maxSize-$twoM; $g++) {

    my $num = 0;
    my $den = 0;

    foreach my $s (@observedSizes) {
	my $p = ($g <= $s-$twoM) ? ($s-$twoM-$g+1)*$sizeDist{$s} : 0;
	$num += $s*$p;
	$den += $p;
    }

    my $meanSpan = $num/$den;
    my $meanAnchor = $meanSpan - $g;
    print "$g\t$meanAnchor\n";

}

	
