# FORWARD

States = c("1","2")
Letters = c("A","T","C","G")
initialProbabilities = c(.72,.28)
transitionMatrix = matrix(c(.57,.43,.56,.44),2)
emissionProbabilities=matrix(c(.37,.13,.21,.44,.22,.17,.20,.26),2)
names(initialProbabilities)   = States
dimnames(transitionMatrix)= list(from=States,to=States)
dimnames(emissionProbabilities)= list(states=States,symbols=Letters)


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

transitionMatrix[is.na(transitionMatrix)]       = 0
emissionProbabilities[is.na(emissionProbabilities)] = 0
sequenceLength  = length(sequence)
statesLength    = length(States)
f          = array(NA,c(statesLength,sequenceLength))
dimnames(f)= list(states=States,index=1:sequenceLength)
# Initialization
for(state in States)
{
  f[state,1] =initialProbabilities[state] * emissionProbabilities[state,sequence[1]]
}
# Iteration
for(k in 2:sequenceLength)
{
  for(state in States)
  {
    temp2 = 0
    for(previousState in States)
    {
      temp   = f[previousState,k-1] * transitionMatrix[previousState,state]
      temp2 = temp2 + temp
    }
    f[state,k] = (emissionProbabilities[state,sequence[k]]) * temp2
  }
}


# Calculate forward probablities
forwardProbabilities=f
finalAnswerF = sum(forwardProbabilities[,4])
print(forwardProbabilities)
print(finalAnswerF)

###############################################################################

# BACKWARD

States = c("1","2")
Letters = c("A","T","C","G")
initialProbabilities = c(.72,.28)
transitionMatrix = matrix(c(.57,.43,.56,.44),2)
emissionProbabilities=matrix(c(.37,.13,.21,.44,.22,.17,.20,.26),2)


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


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


transitionMatrix[is.na(transitionMatrix)]       = 0
emissionProbabilities[is.na(emissionProbabilities)] = 0
sequenceLength  = length(sequence)
statesLength    = length(States)
b          = array(NA,c(statesLength,sequenceLength))
dimnames(b)= list(states=States,index=1:sequenceLength)
# Initialization
for(state in States)
{
  b[state,sequenceLength] = initialProbabilities[state] * emissionProbabilities[state,"G"]
}
# Iteration
for(k in (sequenceLength-1):1)
{
  for(state in States)
  {

    temp2 = 0
    for(nextState in States)
    {
      temp   = b[nextState,k+1] * transitionMatrix[state,nextState]*emissionProbabilities[nextState,sequence[k+1]]
      temp2 = temp2 + temp

    }
    b[state,k] = temp2
  }
}

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


################### FORWARD-BACKWARD

FB = array(NA,c(statesLength,sequenceLength))

#using previously written functions
for(i in 1:sequenceLength)
{
  for(st in 1:statesLength){
    FB[st,i] = (forwardProbabilities[st,i]*backwardProbabilities[st,i])/finalAnswerF
  }
}
finalAnswerFB = sum(FB[,4])
print(finalAnswerFB)