#!/usr/bin/perl -w
#use strict;
##chiamata: perl HMM_es_1.pl



#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 = '3151166666624';

my @obs_list = split(//,$dice_seq);

my %P_emission = (
	'F' => { 1 => 0, 2 => 0, 3 => 0, 4 => 0, 5 => 0, 6 => 0 },
	'L' => { 1 => -0.51, 2 => -0.51, 3 => -0.51, 4 => -0.51 , 5 => -0.51 , 6 => 1.10 } );
my %P_transition = ( 'FF' => 0.64, 'LL' => 0.59, 'FL' => -2.30, 'LF' => -1.61 );

my %path;


########################ALGORITMO DI VITERBI#######################################################


#INIZIALIZZAZIONE TRAMITE UN HASH DI ARRAY:nello stato non truccato ammetto che non ci sia preferenza tra stato F o L e dunque la probabilità di
#andare verso F o L è uguale quindi diventa uno zero; nello stato truccato invece l'emissione è diversa perchè il dado è truccato
my %v = ( 'F' => [0], 'L' => [-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(F L) {
		my ($max_v, $max_prev_state, $max_path) = (-1,'X','XX');

		foreach my $state1 qw(F L) {
			my $tmp_trans = $state1.$state2;


			####### FILL IN #### VITERBI ITERATION STEP ######   
			# my $tmp_v =  

			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 (F o L)
my $last_state = 'X';
if( $v{'F'}->[$idx-1] > $v{'L'}->[$idx-1] ) {
  $last_state = 'F';
} else {
  $last_state = 'L';
}


#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{"F"};
my @Vf = @{ $arrayref };   #estrazione di array da un hash di array
my $Fvit=join("\t", @Vf);
print "F\t$Fvit\n";
$arrayref = $v{"L"};
my @Vl = @{ $arrayref };   #estrazione di array da un hash di array
my $Lvit=join("\t", @Vl);
print "L\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:                                                                                 #
# 3151166666624                                                                           #
# FFFFFFLLLLLLL										  #
#											  #
#Stato finale L (v=10,61)								  #
###########################################################################################
