# 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)


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

 transMatrix[is.na(transMatrix)]       = 0
 emissionMatrices[is.na(emissionMatrices)] = 0
 nObservations  = length(observation)
 nStates    = length(States)
 f          = array(NA,c(nStates,nObservations))
 dimnames(f)= list(states=States,index=1:nObservations)
 # Init
 for(state in States)
 {
   f[state,1] =initProb[state] * emissionMatrices[state,observation[1]]
 }
 # Iteration
 for(k in 2:nObservations)
 {
   for(state in States)
   {
     temp2 = 0
     for(previousState in States)
     {
       temp   = f[previousState,k-1] * transMatrix[previousState,state]
       temp2 = temp2 + temp
     }
     f[state,k] = (emissionMatrices[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)


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


 transMatrix[is.na(transMatrix)]       = 0
 emissionMatrices[is.na(emissionMatrices)] = 0
 nObservations  = length(observation)
 nStates    = length(States)
 b          = array(NA,c(nStates,nObservations))
 dimnames(b)= list(states=States,index=1:nObservations)
 # Init
 for(state in States)
 {
   b[state,nObservations] = initProb[state] * emissionMatrices[state,"G"]
 }
 # Iteration
 for(k in (nObservations-1):1)
 {
   for(state in States)
   {

     temp2 = 0
     for(nextState in States)
     {
       temp   = b[nextState,k+1] * transMatrix[state,nextState]*emissionMatrices[nextState,observation[k+1]]
       temp2 = temp2 + temp

     }
     b[state,k] = temp2
   }
 }

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