@counts = (
	   [191,309,488,838,370,167,247,722,809,328,376,412,362,131], # G
	   [149,359,370,391,365,827,666,266,180,398,863,471,302,192], # C
	   [764,369,487,327,237,159,152,506,483,704,116,297,618,577], # A
	   [556,623,315,104,688,507,595,166,188,230,305,480,378,760], # T
	   [  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0], # N
	  );

for($j = 0; $j < 14; $j++)
{
    $total = 0;
    for($i = 0; $i < 4; $i++)
    {
	$total += $counts[$i][$j];
    }
    for($i = 0; $i < 5; $i++)
    {
	$pn[$i][$j] = $counts[$i][$j]/$total;
    }
}

$pointer = 0;
@circular = (1,0,3,3,2,2,2,3,1,3,2,2,1,0);
$p_min = &probability();
@circular = (2,3,0,0,3,1,1,0,0,2,1,3,2,3);
$p_max = &probability();

until(<> =~ /^ORIGIN/)
{
    # eat lines
}

$Nbases = $nbases = &readbases();

for($i = 0; $i <= 100; $i++)
{
    $histogram[$i] = 0;
}

@circular = (4,4,4,4,4,4,4,4,4,4,4,4,4,4);
$pointer = 0;
while($nbases > 0)
{
    $circular[$pointer % 14] = $buffer[$Nbases - $nbases];
    $pointer++;
    $nbases--;
    $p = &probability();
    if($p != 0)
    {
	$l = -log($p/$p_max)/log(2);
	$bin = int(100*$l/24.795715);
	$histogram[$bin]++;
    }
    if($nbases == 0)
    {
	$Nbases = $nbases = &readbases();
    }
}

for($i = 0; $i <= 100; $i++)
{
    $bin = 24.795715*$i/100;
    $hist = $histogram[$i];
    write;
}

sub readbases
{
    local($incoming);
    $incoming = <>;
    chop($incoming);
    $incoming =~ s/ |\d//g;
    $incoming =~ tr/GCATgcat/4/c;
    $incoming =~ tr/GCATgcat/01230123/;
    @buffer = split(//, $incoming);
    length($incoming);
}

sub probability
{
    local($prob);
    $prob = 1;
    for (0,1,3,4,5,6,7,8,9,10,12,13)
    {
	$prob *= $pn[$circular[($pointer + $_) % 14]][$_];
    }
    $prob;
}

format STDOUT = 
@##.### @######.#
$bin $hist
.
