#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "nrutil.h"
#include "hmm.h"

void GenSequenceArray(HMM *phmm, int T, int *O, int *q)
{
        int     t = 1;
        int     q_t, o_t;

    	srand(1000);

        q[1] = GenInitalState(phmm);
        O[1] = GenSymbol(phmm, q[1]);

        for (t = 2; t <= T; t++) {
                q[t] = GenNextState(phmm, q[t-1]);
                O[t] = GenSymbol(phmm, q[t]);
        }
}

int GenInitalState(HMM *phmm)
{
        double val, accum;
        int i, q_t;

        val = rand();
        accum = 0.0;
        q_t = phmm->N;
        for (i = 1; i <= phmm->N; i++) {
                if (val < phmm->pi[i] + accum) {
                        q_t = i;
                        break;
                }
                else {
                                accum += phmm->pi[i];
                }
        }

        return q_t;
}

int GenNextState(HMM *phmm, int q_t)
{
        double val, accum;
        int j, q_next;

        val = rand();
        accum = 0.0;
        q_next = phmm->N;
        for (j = 1; j <= phmm->N; j++) {
                if ( val < phmm->A[q_t][j] + accum ) {
                        q_next = j;
                        break;
                }
                else
                        accum += phmm->A[q_t][j];
        }

        return q_next;
}

int GenSymbol(HMM *phmm, int q_t)
{
        double val, accum;
        int j, o_t;

        val = rand();
        accum = 0.0;
        o_t = phmm->M;
        for (j = 1; j <= phmm->M; j++) {
                if ( val < phmm->B[q_t][j] + accum ) {
                       o_t = j;
                       break;
                }
                else
                        accum += phmm->B[q_t][j];
        }
        return o_t;
}

void ReadSequence(FILE *fp, int *pT, int **pO)
{
        int *O;
        int i;

        fscanf(fp, "T= %d\n", pT);
        O = ivector(1,*pT);
        for (i=1; i <= *pT; i++)
                fscanf(fp,"%d", &O[i]);
        *pO = O;
}

void PrintSequence(FILE *fp, int T, int *O)
{
        int i;

        fprintf(fp, "T= %d\n", T);
        for (i=1; i <= T; i++) 
                fprintf(fp,"%d ", O[i]);
	printf("\n");
}