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

void ReadHMM(FILE *fp, HMM *phmm){
	int i, j, k;

	fscanf(fp, "M= %d\n", &(phmm->M));
	fscanf(fp, "N= %d\n", &(phmm->N));
	fscanf(fp, "A:\n");
	phmm->A = (double **) dmatrix(1, phmm->N, 1, phmm->N);
	for (i = 1; i <= phmm->N; i++) {
		for (j = 1; j <= phmm->N; j++) {
			fscanf(fp, "%lf", &(phmm->A[i][j]));
		}
		fscanf(fp,"\n");
	}

	fscanf(fp, "B:\n");
	phmm->B = (double **) dmatrix(1, phmm->N, 1, phmm->M);
	for (j = 1; j <= phmm->N; j++) {
		for (k = 1; k <= phmm->M; k++) {
			fscanf(fp, "%lf", &(phmm->B[j][k])); 
		}
		fscanf(fp,"\n");
	}
	fscanf(fp, "pi:\n");
	phmm->pi = (double *) dvector(1, phmm->N);
	for (i = 1; i <= phmm->N; i++) 
		fscanf(fp, "%lf", &(phmm->pi[i])); 
}

void FreeHMM(HMM *phmm)
{
	free_dmatrix(phmm->A, 1, phmm->N, 1, phmm->N);
	free_dmatrix(phmm->B, 1, phmm->N, 1, phmm->M);
	free_dvector(phmm->pi, 1, phmm->N);
}

void InitHMM(HMM *phmm, int N, int M)
{
	int i, j, k;
	double sum;
        srand(1000);
       	phmm->M = M;
        phmm->N = N;
        phmm->A = (double **) dmatrix(1, phmm->N, 1, phmm->N);
        for (i = 1; i <= phmm->N; i++) {
		sum = 0.0;
                for (j = 1; j <= phmm->N; j++) {
                        phmm->A[i][j] = rand(); 
			sum += phmm->A[i][j];
		}
                for (j = 1; j <= phmm->N; j++) 
			 phmm->A[i][j] /= sum;
	}

        phmm->B = (double **) dmatrix(1, phmm->N, 1, phmm->M);

        for (j = 1; j <= phmm->N; j++) {
		sum = 0.0;	
                for (k = 1; k <= phmm->M; k++) {
                        phmm->B[j][k] = rand();
			sum += phmm->B[j][k];
		}
                for (k = 1; k <= phmm->M; k++) 
			phmm->B[j][k] /= sum;
	}
 
        phmm->pi = (double *) dvector(1, phmm->N);
	sum = 0.0;
        for (i = 1; i <= phmm->N; i++) {
                phmm->pi[i] = rand(); 
		sum += phmm->pi[i];
	}
        for (i = 1; i <= phmm->N; i++) 
		phmm->pi[i] /= sum;
}

void PrintHMM(FILE *fp, HMM *phmm)
{
        int i, j, k;
	fprintf(fp, "M= %d\n", phmm->M);
	fprintf(fp, "N= %d\n", phmm->N);
	fprintf(fp, "A:\n");
        for (i = 1; i <= phmm->N; i++) {
                for (j = 1; j <= phmm->N; j++) {
                        fprintf(fp, "%f ", phmm->A[i][j] );
		}
		fprintf(fp, "\n");
	}

	fprintf(fp, "B:\n");
        for (j = 1; j <= phmm->N; j++) {
                for (k = 1; k <= phmm->M; k++){
                        fprintf(fp, "%f ", phmm->B[j][k]);
		}
		fprintf(fp, "\n");
	}

	fprintf(fp, "pi:\n");
        for (i = 1; i <= phmm->N; i++) {
		fprintf(fp, "%f ", phmm->pi[i]);
	}

	fprintf(fp, "\n\n");
}