#!/usr/bin/perl -w
#
# Count 3letters words for HMM model
#--------------------------------------- 
use strict;
use warnings;
my $c = 0;
# number of hidden states
my $H          = 4;
# table % of similarity separators between states
my @pH         = (0.0,0.60,0.70,0.80,1.00);
# window to evaluate % of similarity
my $winsize    = 32;


# (0) take at least one yass output file
( $#ARGV  >= 0) or die("\"give at least one filename on commandline\"");

# global variables
my %code       = ("|", 2, "." , 1 , ":" , 0 );
my %error      = ("|", 0, "." , 1 , ":" , 1 );
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 tables
# p table
my @fullcount = ();
for my $h (0..$H-1) {
    for my $code   (0..4**3) {
	$fullcount[$h][$code] = 0;
    }
}

# q table
my @fullswitch = ();
for my $h1 (0..$H-1) {
    for my $h2 (0..$H-1) {
	$fullswitch[$h1][$h2] = 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;

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

		# (b) count and bias done here
		my $winerr = 0;	
		my $ch = $H - 1;
		for my $j (0..length($alignment)-3-1) {

		    # compute number of errors inside the window
		    if (substr($alignment,$j,1) eq "|"){
			$winerr++;
		    }
		    if ($j-$winsize >= 0) {
			if (substr($alignment,$j-$winsize,1) eq "|"){
			    $winerr--;
			}

			# if first pos
			if (!($j%3)) {
			    my $oldch = $ch;
			    # reset current hidden state H
			    for my $h (0..$H-1){
				#print $winerr / $winsize."\n"; 
				if ( $pH[$h] <= ($winerr / $winsize) ){
				    $ch = $h;
				}
			    }
			    $swicth[$oldch][$ch]++;
			}
		    }
		    
		    #compute the bias
		    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]++;
		    } 


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

		

		#(c) alignment proceeded : compute full sums

		# take the maxbias
		my $maxbias = 0;
		my $sumbias = 0;
		my $delta   = 0;
		for my $i (0..3-1){
		    if( $bias[$i] > $maxbias ){
			$maxbias = $bias[$i];
			$delta   = $i;
		    }
		    $sumbias += $bias[$i];
		}
		
		# (d) sum all the counts 
		if ($maxbias >= 0.4 * $sumbias && $maxbias <= 0.55 * $sumbias && $maxbias >= 10 && $trtv[0] <=  $trtv[1]){
		    
		    #update fulltable
		    for  my $code (0..3**3-1){
			$fullcount[$ch][$code] += $count[$ch][$code][$delta];
		    }
		    #update full switch
		    for my $h1 (0..$H-1) {
			for my $h2 (0..$H-1) {
			    $fullswitch[$h1][$h2] += $swicth[$h1][$h2];
			}
		    }
		}
		$alignment = "";
	    }
	}
    }
    close(IN);
}


# init
my @sumcount = ();
for my $h (0..$H-1) {
    $sumcount[$h] = 0;
}
my @sumswitch = ();
for my $h (0..$H-1) {
    $sumswitch[$h] = 0;
}

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


# output HMM (H part) 
if ($c) {
    print("double H[".($H*$H)."] =  {\n");
}
for my $h1 (0..$H-1) {
    if ($sumswitch[$h1] == 0) {
	$sumswitch[$h1] = 1;
    }
    if ($c) {
	print "\t /* q".$h1." */\t";
    }
    for my $h2 (0..$H-1) { 
	print "".($fullswitch[$h1][$h2]/($sumswitch[$h1]))."," ;
    }
    if ($c) {
	print "\n";
    }
} 
if ($c) {
    print("};\n");
}
# output HMM (M13 part)
for my $h (0..$H-1) {   
    if ($sumcount[$h] == 0) {
	$sumcount[$h] = 1;
    } 
    if ($c) {
	print("double WM".$h."[27] =  {\n");
    }
    for my $cw (0..3**3-1){
	if ($c) {
	    print("\t /* ".$rcode{$cw}." */\t".(($fullcount[$h][$cw])/($sumcount[$h])).",\n");
	} else {
	    print("".(($fullcount[$h][$cw])/($sumcount[$h])).",");
	}
    }
    if ($c) {
	print("};\n");
    }
}

if ($c) {
    print "double * WMM[".$H."] = {\n";
    for my $h (0..$H-1) {
	print "\t WM".$h.",\n";
    }
    print "};\n";
}

