#include <stdio.h>
#include "hmm.h"
#include "nrutil.h"

void Forward(HMM *phmm, int T, int *O, double **alpha, double *pprob)
{
    int     i, j;
    int     t;
    double sum;

    for(i= 1; i <= phmm->N; i++)
	alpha[1][i] = phmm->pi[i]* phmm->B[i][O[1]];

    for (t = 1; t < T; t++)
    {
	for (j = 1; j <= phmm->N; j++)
	{
	     sum = 0.0;
	     for (i = 1; i <= phmm->N; i++)
		sum += alpha[t][i]* (phmm->A[i][j]);
		alpha[t+1][j] = sum*(phmm->B[j][O[t+1]]);
	}
     }

     *pprob = 0.0;
     for (i = 1; i <= phmm->N; i++)
         *pprob += alpha[T][i];
}

void ForwardWithScale(HMM *phmm, int T, int *O, double **alpha,	double *scale, double *pprob)
{
     int i, j;

     int t;
     double sum;

     scale[1] = 0.0;
     for (i = 1; i <= phmm->N; i++)
     {
	alpha[1][i] = phmm->pi[i]* (phmm->B[i][O[1]]);
	scale[1] += alpha[1][i];
     }

     for (i = 1; i <= phmm->N; i++)
	alpha[1][i] /= scale[1];

     for (t = 1; t <= T - 1; t++)
     {
	scale[t+1] = 0.0;
	for (j = 1; j <= phmm->N; j++)
	{
	   sum = 0.0;
	   for (i = 1; i <= phmm->N; i++)
	      sum += alpha[t][i]* (phmm->A[i][j]);

	   alpha[t+1][j] = sum*(phmm->B[j][O[t+1]]);
	   scale[t+1] += alpha[t+1][j];
	}

	for (j = 1; j <= phmm->N; j++)
	   alpha[t+1][j] /= scale[t+1];
      }

      *pprob = 0.0;
      for (t = 1; t <= T; t++)
	 *pprob += log(scale[t]);
}