# FORWARD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 # Symbols = c("A","T","C","G")
 # initProb = c(.65,.35)
 # transMatrix = matrix(c(.68,.55,.32,.45),2)
 # emissionMatrices=matrix(c(.25,.45,.11,.19,.14,.24,.47,.15),2)
 #
States = c("1","2")
Symbols = c("A","T","C","G")
initProb = c(.72,.28)
transMatrix = matrix(c(.57,.43,.56,.44),2)
emissionMatrices=matrix(c(.37,.13,.21,.44,.22,.17,.20,.26),2)
 names(initProb)   = States
 dimnames(transMatrix)= list(from=States,to=States)
 dimnames(emissionMatrices)= list(states=States,symbols=Symbols)


 hmm = list(States=States,Symbols=symbols,startProbs=initProb,transProbs=transMatrix,
            emissionProbs=emissionMatrices)

 observation = c("G","T","T","C","G","A","T","G")

 hmm$transProbs[is.na(hmm$transProbs)]       = 0
 hmm$emissionProbs[is.na(hmm$emissionProbs)] = 0
 nObservations  = length(observation)
 nStates    = length(hmm$States)
 f          = array(NA,c(nStates,nObservations))
 dimnames(f)= list(states=hmm$States,index=1:nObservations)
 # Init
 for(state in hmm$States)
 {
   f[state,1] =hmm$startProbs[state] * hmm$emissionProbs[state,observation[1]]
 }
 # Iteration
 for(k in 2:nObservations)
 {
   for(state in hmm$States)
   {
     temp2 = 0
     for(previousState in hmm$States)
     {
       temp   = f[previousState,k-1] * hmm$transProbs[previousState,state]
       temp2 = temp2 + temp
     }
     f[state,k] = (hmm$emissionProbs[state,observation[k]]) * temp2
   }
 }


 observations = c("A","T","G","C","C","A","T","G")
 # Calculate forward probablities
 forwardProbabilities=f
 finalAnswer = sum(forwardProbabilities[,4])
 print(finalAnswer)

 print(f)
 # final = sum(f[,4])
 # print(final)
 ###############################################################################

 # BACKWARD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 States = c("1","2")
 Symbols = c("A","T","C","G")
 initProb = c(.72,.28)
 transMatrix = matrix(c(.57,.43,.56,.44),2)
 emissionMatrices=matrix(c(.37,.13,.21,.44,.22,.17,.20,.26),2)


 names(initProb)   = States
 dimnames(transMatrix)= list(from=States,to=States)
 dimnames(emissionMatrices)= list(states=States,symbols=Symbols)


 hmm = list(States=States,Symbols=symbols,startProbs=initProb,transProbs=transMatrix,
            emissionProbs=emissionMatrices)

 observation = c("G","T","T","C","G","A","T","G")


hmm$transProbs[is.na(hmm$transProbs)]       = 0
hmm$emissionProbs[is.na(hmm$emissionProbs)] = 0
nObservations  = length(observation)
nStates    = length(hmm$States)
b          = array(NA,c(nStates,nObservations))
dimnames(b)= list(states=hmm$States,index=1:nObservations)
# Init
for(state in hmm$States)
{
 b[state,nObservations] = hmm$startProbs[state] * hmm$emissionProbs[state,observation[1]]
}
# Iteration
for(k in (nObservations-1):1)
{
 for(state in hmm$States)
 {

   temp2 = 0
   for(nextState in hmm$States)
   {
     temp   = b[nextState,k+1] * hmm$transProbs[state,nextState]*hmm$emissionProbs[nextState,observation[k+1]]
     temp2 = temp2 + temp

   }
   b[state,k] = temp2
 }
}

 print(b)
 backwardProbabilities=b
 finalAnswer = sum(backwardProbabilities[,4])
 print(finalAnswer)