#!/usr/bin/perl -w
#
# Count 3 letters words for M13 model
#--------------------------------------- 
use strict;
use warnings;
my $c = 0;
# (0) take at least one yass output file
( $#ARGV  >= 0) or die("\"give at least one filename on commandline\"");

# (1) clean the table
my %code       = ("|", 2, "." , 1 , ":" , 0 );
my %rcode      = (
		  0,":::", 
		  1,"::.",
		  2,"::|",
		  3,":.:", 
		  4,":..",
		  5,":.|",
		  6,":|:", 
		  7,":|.",
		  8,":||",
		  9,".::", 
		  10,".:.",
		  11,".:|",
		  12,"..:", 
		  13,"...",
		  14,"..|",
		  15,".|:", 
		  16,".|.",
		  17,".||",
		  18,"|::", 
		  19,"|:.",
		  20,"|:|",
		  21,"|.:", 
		  22,"|..",
		  23,"|.|",
		  24,"||:", 
		  25,"||.",
		  26,"|||"
		  );

# (1) clean the table
my @fullcount = ();
for my $code (0..4**3-1){
	$fullcount[$code] = 0;
}

# read the data on each file
foreach my $filename (@ARGV){
    my $alignment = "";
    open(IN, $filename) or die("Cant open \"$filename\"");
   
    while(<IN>){
	if ($_ =~ /^[|.: ]*$/ ){ 
            # (1) select only alignment fragments
	    chomp $_;
	    $alignment .= $_;
	}else{ 
	   
            # (2) new alignment found : process $alignment computed before
	    if ($_ =~ /^\*\(.*$/ ){ # new alignment found
		#remove end spaces of alignments
		$alignment =~ s/\s+$//;		
		#(temporary) : remove indels
		$alignment =~ s/\s+//g;

		#clear table of counts
		my @count = ();
		for my $code (0..3**3-1){
		    $count[$code][0] = 0;
		    $count[$code][1] = 0;
		    $count[$code][2] = 0;
		}
		
		my @bias = ();
		 for my $code (0..3-1){
		    $bias[$code] = 0;
		}
		my @trtv = ();
		 for my $code (0..2){
		    $trtv[$code] = 0;
		}

		
		# code computed here
		for my $j  (0..length($alignment)-3-1){
		    if (!(substr($alignment,$j,1) eq "|")){
			$bias[$j%3]++;
		    }    
		    
		    #trtv ratio
		    if ((substr($alignment,$j,1) eq ":")){
			$trtv[0]++;
		    }  
		    
		    if ((substr($alignment,$j,1) eq ".")){
			$trtv[1]++;
		    } 


		    my $cd = 0;
		    for my $i (0..3-1){
			$cd *= 3;
			$cd += $code{substr($alignment,$j+2-$i,1)};
		    }
		    $count[$cd][$j%3]+=1;
		}

		#take the maxbias
		my $maxbias = 0;
		my $sumbias = 0;
		my $delta   = 0;
		for my $i (0..3-1){
		    if($maxbias < $bias[$i]){
			$maxbias = $bias[$i];
			$delta   = $i;
		    }
		    $sumbias += $bias[$i];
		}

		#
		if ($maxbias >= 0.4 * $sumbias && $maxbias <= 0.55 * $sumbias && $maxbias >= 10 && $trtv[0] <=  $trtv[1]){
		    #print("".$bias[0]." ".$bias[1]." ".$bias[2]."\n");
		    #update fulltable
		    for my $code (0..3**3-1){
			$fullcount[$code] += $count[$code][$delta];
		    }
		}
		$alignment = "";
	    }
	}
    }
    close(IN);
}

# compute probs
my $sumcount = 0;
for my $cw (0..3**3-1){
    $sumcount +=  $fullcount[$cw];
}

if($sumcount == 0){
    $sumcount = 1;
}
if ($c){
    print("double W[27] = {\n");
}
for my $cw (0..3**3-1){
    if ($c){
	print("/* ".$rcode{$cw}." */\t".$fullcount[$cw]/$sumcount.",\n");
    } else {
	print($fullcount[$cw]/$sumcount.",");
    }
}
if ($c){
    print("};\n");
}
