#!/usr/bin/perl -w
#use strict;
#eseguo: perl HMM_es_4.pl state_sequences.txt



#DICHIARAZIONE DELLE VARIABILI: $dice_seq è la sequenza, @obs_list è un array che contiene ogni singolo valore che compone $dice_seq, my $P_emission
#è un hash di hash (hash a due chiavi) e contiene le  probabilità di emissione, my %P_transition è un hash a chiave singola e contiene le probabilità
# di transizione 


my $dice_seq = 'KKIIFFFFL';

my @obs_list = split(//,$dice_seq);

my %P_emission = (
	'S' => { 'F' => 0.047 , 'S' => 0.096 , 'T' => 0.061 , 'N' => 0.059 ,  'K' => 0.071 ,  'Y' => 0.039 , 'E' => 0.063, 'V' => 0.053 , 'Q' => 0.037 , 'M' => 0.021 , 'C' => 0.012 , 'L' => 0.089 , 'A' => 0.054 , 'W' => 0.014 , 'P' => 0.043 , 'H' => 0.023 , 'D' => 0.058 , 'I' => 0.059 , 'R' => 0.050 , 'G' => 0.051},
	'T' => {'S' => 0.068 , 'F' => 0.107 , 'T' => 0.052 , 'N' => 0.021 , 'K' => 0.009 , 'E' => 0.010 , 'Y' => 0.041 , 'V' => 0.098 , 'Q' => 0.014 , 'M' =>0.032 , 'C' => 0.027 , 'L' => 0.158 , 'A' => 0.084 , 'W' => 0.023 , 'P' =>0.029 , 'H' => 0.008 , 'D' => 0.009 , 'I' => 0.129 , 'R' => 0.006, 'G' => 0.076 } );
my %P_transition = ( 'SS' => 0.209 , 'TT' => 0.005, 'ST' => 0.029 , 'TS' => 0.028 );

my %path;


########################ALGORITMO DI VITERBI#######################################################


#INIZIALIZZAZIONE TRAMITE UN HASH DI ARRAY
my %v = ( 'S' => [0], 'T' => [-0.51] );

#Siamo nella prima posizione:
my $idx = 1;

#ITERAZIONE
foreach my $tmp_obs (@obs_list) {
  if($idx > $#obs_list){$obsidx = "END";}else{$obsidx=$obs_list[$idx];}
  print "i =  ".($idx)." obs_i = $obsidx i-1 = ".($idx-1)." obs_i-1 = $obs_list[$idx-1]\n";
  my $max_state = '';   #con il primo ciclo foreach per ogni osservazione settiamo a zero il massimo stato di viterbi
  my $max_state_v = 0;
	foreach my $state2 qw(S T) {
		my ($max_v, $max_prev_state, $max_path) = (-1,'X','XX');

		foreach my $state1 qw(S T) {
			my $tmp_trans = $state1.$state2;


			  
			my $tmp_v =  $P_emission{$state2}->{$tmp_obs} + $P_transition{$tmp_trans} + $v{$state1}->[$idx-1];

			print sprintf("[%d]%s -> %s: %.2f * %.2f * %.2f = %.2f\n", 
											$idx, $state1, $state2, $P_emission{$state2}->{$tmp_obs}, $P_transition{$tmp_trans},
											$v{$state1}->[$idx-1], $tmp_v);											
#nel ciclo if testo da dove arriva il valore state2
			if( $tmp_v > $max_v ) {
				$max_v = $tmp_v; 
        $max_prev_state = $state1;
				$max_path = $state1.$state2;
			}
		}
		print sprintf("MAX PATH to %s: %s, %s%d = %.2f\n", $state2, $max_path, 
										$state2,$idx+1,$max_v);
	      
    $path{$idx}->{$state2} = $max_prev_state;  #salvataggio delle variabili per lo step successivo

		$v{$state2}->[$idx] = $max_v;

	}
	print "\n";
	$idx += 1;    #EQUIVALE A DIRE $idx++
}

#DETERMINAZIONE DELLO STATO FINALE: col ciclo if..else controllo i valori dell'ultima colonna e il valore maggiore determina lo stato finale (S o T)
my $last_state = 'X';
if( $v{'S'}->[$idx-1] > $v{'T'}->[$idx-1] ) {
  $last_state = 'S';
} else {
  $last_state = 'T';
}


#Stampa lo stato finale e il relativo valore
print "Final state: $last_state (v= $v{$last_state}->[$idx-1])\n";

# TRACEBACK: abbiamo il valore massimo dell'ultima colonna; al primo giro usa il valore dell'ultima colonna, prende lo stato precedente e lo mette
#nel vettore degli stati; poi fa un update così che al prossimo giro parte dal valore precedente. L'ordine degli stati emessi dal traceback
#devono poi essere invertiti con un reverse

my @state_seq = ();
for(my $i=$idx-1;$i>0;$i--) {
  my $prev_state = $path{$i}->{$last_state};
  push(@state_seq,$prev_state);
  $last_state = $prev_state;
}

# Stampa la sequenza di stati che ha prodotto la sequenza di oservazioni $dice_seq con maggiore probabilità
print "\n[Output]\n";
print join("",@obs_list),"\n";
print join("",reverse @state_seq),"\n";


##### OUTPUT ADDITIONAL INFO ###########################################

##### STAMPA LA MATRICE DI VITERBI

print "\n\nVITERBI matrix:\n\n";
my $jseq = join("   \t", @obs_list);
print "\tB\t$jseq\n";
my $arrayref = $v{"S"};
my @Vf = @{ $arrayref };   #estrazione di array da un hash di array
my $Fvit=join("\t", @Vf);
print "S\t$Fvit\n";
$arrayref = $v{"T"};
my @Vl = @{ $arrayref };   #estrazione di array da un hash di array
my $Lvit=join("\t", @Vl);
print "T\t$Lvit\n\n";

##### STAMPA LA MATRICE DI TRACEBACK ###############
print "\n\nTRACEBACK pointers:\n\n";
    my $rHoH = \%path;
    my $curpos=0;
    for my $k1 ( sort {$a <=> $b} (keys %$rHoH) ) {
        print "Step: $k1  \tEmission: $obs_list[$curpos]\n";
        for my $k2 ( keys %{$rHoH->{ $k1 }} ) {
            print "\tmax path to $k2 :  $rHoH->{ $k1 }{ $k2 }\t";
			if($k2 ne $rHoH->{ $k1 }{ $k2 }){print "*\n";}else{print "\n";}
        }
		$curpos++;
		print "\n";
    }


#####################################################################################
# OUTPUT									    #
# KKIIFFFFL									    #	
# SSSSSSSSS									    #	
# stato finale: S(v=2.418)							    #
#####################################################################################