# Viterbi algorithm ---------------------------------------------------------------


# initializations ---------------------------------------------------------


States = c("1","2","3","4")
Letters = c("Giraffes","Live","In","The","Savannah")
iter = 0
transitionMatrix = matrix(c(.22,0,.41,.32,.38,.21,.09,0,.15,.42,.3,.29,.25,.37,.2,.39),4)
#emissionProbabilities=matrix(c(.13,.27,.42,.11,.09,.31,.09,.24,.06,.3,.15,.21,.3,.12,.22,.08,.13,.16,.47,.16),5)
emissionProbabilities=matrix(c(.13,.31,.15,.08,.27,.09,.21,.13,.42,.24,.3,.16,.11,.06,.12,.47,.09,.3,.22,.16),4)
sequence = c("Giraffes","Live","In","The","Savannah")
initialProbabilities = rep(1/2,4)
initialProbabilities[1]=0.23
initialProbabilities[2]=0.31
initialProbabilities[3]=0.27
initialProbabilities[4]=0.19

names(initialProbabilities)   = States
dimnames(transitionMatrix)= list(from=States,to=States)
dimnames(emissionProbabilities)= list(states=States,symbols=Letters)


transitionMatrix[is.na(transitionMatrix)]       = 0
emissionProbabilities[is.na(emissionProbabilities)] = 0
diff = c()
#################################################################
transitionMatrix[is.na(transitionMatrix)]       = 0
emissionProbabilities[is.na(emissionProbabilities)] = 0
nSequence=length(sequence)
nStates=length(States)
v= array(NA,c(nStates,nSequence))
dimnames(v)= list(states=States,index=1:nSequence)
# Init
for(state in States)
{
  v[state,1] =initialProbabilities[state] * emissionProbabilities[state,sequence[1]]
}
# Iteration
for(k in 2:nSequence)
{
  for(state in States)
  {
    maxi = NULL
    for(previousState in States)
    {
      temp = v[previousState,k-1] * transitionMatrix[previousState,state] * emissionProbabilities[state,sequence[k]]
      maxi = max(maxi, temp)
    }
    v[state,k] = maxi
  }
}
# Traceback
finalVitPath = rep(NA,nSequence)
for(state in States)
{
  if(max(v[,nSequence])==v[state,nSequence])
  {
    finalVitPath[nSequence] = state
    break
  }
}
for(k in (nSequence-1):1)
{
  for(state in States)
  {
    if(max(v[,k]*transitionMatrix[,finalVitPath[k+1]])==v[state,k]*transitionMatrix[state,finalVitPath[k+1]])
    {
      finalVitPath[k] = state
      break
    }
  }
}
######################################################################

# final answer ------------------------------------------------------------

print(finalVitPath)