bwfw.hmm<-function(x, pii, A, pc, f0, f1)
{

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

## USAGE
 # bwfw.hmm(x, pii, A, pc, f0, f1)

## ARGUMENTS
 # x=(x[1], ..., x[m]): the observed data
 # pii=(pii[0], pii[1]): the initial state distribution
 # A=(A[0,0], A[0,1]\\ A[1,0], A[1,1]): the transition matrix
 # f0=(mu, sigma): the parameters for null distribution
 # pc=(c[1], ..., c[L])
 # --the probability weights in the mixture for each component
 # f1=(mu[1], sigma[1]\\...\\mu[L], sigma[L])
 # --the parameter set for the non-null distribution

## DETAILS
 # bwfw.hmm calculates values for backward, forward variables, probabilities of hidden states, 
 # --the lfdr variables and etc. 
 # --using the forward-backward procedure (Baum et al.) 
 # --based on a sequence of observations for a given hidden markov model M=(pii, A, f)
 # the underflow problem was fixed by using the rescaled forward and backward variables
 # -- see Stamp (2004) for a detailed instruction on the coding of this algorithm

## VALUES
 # alpha: rescaled backward variables
 # beta: rescaled forward variables
 # lfdr: lfdr variables
 # gamma: probabilities of hidden states
 # dgamma: rescaled transition variables
 # omega: rescaled weight variables

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

## Initialize

NUM<-length(x)
L<-length(pc)

## Densities
f0x<-dnorm(x, f0[1], f0[2])
f1x<-rep(0, NUM)
for (c in 1:L)
{
  f1x<-f1x+pc[c]*dnorm(x, f1[c, 1], f1[c, 2])
}

## the backward-forward procedure

# a. the backward variables
# --rescaled 

alpha<-matrix(rep(0, NUM*2), NUM, 2, byrow=T)
# scaling variable c_0
c0<-rep(0, NUM)

alpha[1, 1]<-pii[1]*f0x[1]
alpha[1, 2]<-pii[2]*f1x[1]
# rescaling alpha
c0[1]<-1/sum(alpha[1, ])
alpha[1, ]<-c0[1]*alpha[1, ]

for (k in 1:(NUM-1))
{ 
  alpha[k+1, 1]<-(alpha[k, 1]*A[1, 1]+alpha[k, 2]*A[2, 1])*f0x[k+1]
  alpha[k+1, 2]<-(alpha[k, 1]*A[1, 2]+alpha[k, 2]*A[2, 2])*f1x[k+1]
  # rescaling alpha
  c0[k+1]<-1/sum(alpha[k+1, ])
  alpha[k+1, ]<-c0[k+1]*alpha[k+1, ]
}

# b. the forward variables
# --rescaled

beta<-matrix(rep(0, NUM*2), NUM, 2, byrow=T)

beta[NUM, 1]<-c0[NUM]
beta[NUM, 2]<-c0[NUM]

for (k in (NUM-1):1)
{ 
  beta[k, 1]<-A[1, 1]*f0x[k+1]*beta[k+1, 1]+A[1, 2]*f1x[k+1]*beta[k+1, 2]
  beta[k, 2]<-A[2, 1]*f0x[k+1]*beta[k+1, 1]+A[2, 2]*f1x[k+1]*beta[k+1, 2]
  # rescaling beta
  # using the same scaling factors as alpha 
  beta[k, ]<-c0[k]*beta[k, ]
}

# c. lfdr variables
# --original
# --the same formulae hold for the rescaled alpha and beta

lfdr<-rep(0, NUM)

for (k in 1:NUM)
{ 
  q1<-alpha[k, 1]*beta[k, 1]
  q2<-alpha[k, 2]*beta[k, 2]
  lfdr[k]<-q1/(q1+q2)
}

# d. probabilities of hidden states
# -- and transition variables
# -- both are rescaled
 
gamma<-matrix(1:(NUM*2), NUM, 2, byrow=T)
# initialize gamma[NUM]
gamma[NUM, ]<-c(lfdr[NUM], 1-lfdr[NUM])
dgamma<-array(rep(0, (NUM-1)*4), c(2, 2, (NUM-1)))

for (k in 1:(NUM-1))
{
  denom<-0
  for (i in 0:1)
  {
    for (j in 0:1)
    { 
      fx<-(1-j)*f0x[k+1]+j*f1x[k+1]
      denom<-denom+alpha[k, i+1]*A[i+1, j+1]*fx*beta[k+1, j+1]
    }
  }
  for (i in 0:1)
  {
    gamma[k, i+1]<-0
    for (j in 0:1)
    { 
      fx<-(1-j)*f0x[k+1]+j*f1x[k+1]
      dgamma[i+1, j+1, k]<-alpha[k, i+1]*A[i+1, j+1]*fx*beta[k+1, j+1]/denom
      gamma[k, i+1]<-gamma[k, i+1]+dgamma[i+1, j+1, k]  
    }
  }
}

# e. weight variables
# --rescaled
# --the same formula holds using the rescaled gamma

omega<-matrix(rep(0, NUM*L), NUM, L, byrow=T)

for (k in 1:NUM)
{
  for (c in 1:L)
  { 
    f1c<-dnorm(x[k], f1[c, 1], f1[c, 2])
    omega[k, c]<-gamma[k, 2]*pc[c]*f1c/f1x[k]
  }
}

# f. return the results of the bwfw proc.

bwfw.var<-list(bw=alpha, fw=beta, lsi=lfdr, pr=gamma, ts=dgamma, wt=omega)

return(bwfw.var)
  
}