#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "nrutil.h"
#include "hmm.h"
#include "forward.c"
#include "nrutil.c"
#include "hmmutils.c"
#include "hmmrand.c"
#include "sequence.c"
#include "viterbi.c"

int main (int argc, char **argv)
{
	int 	t, T;
	HMM  	hmm;
	int	*O;
	int	*q;
	double **delta;
	int	**psi;
	double 	proba, logproba; 
	FILE	*fp;
	if (argc != 3) {
		printf("Usage error \n");
		printf("Usage: testvit <model.hmm> <obs.seq> \n");
		exit (1);
	}

	fp = fopen(argv[1], "r");
	if (fp == NULL) {
		fprintf(stderr, "Error: File %s not found\n", argv[1]);
		exit (1);
	}
	ReadHMM(fp, &hmm);
	fclose(fp);

	fp = fopen(argv[2], "r");
	if (fp == NULL) {
		fprintf(stderr, "Error: File %s not found\n", argv[2]);
		exit (1);
	}
	ReadSequence(fp, &T, &O);
	fclose(fp);
	q = ivector(1,T);
	delta = dmatrix(1, T, 1, hmm.N);
	psi = imatrix(1, T, 1, hmm.N);

	printf("------------------------------------\n");
	printf("Viterbi using direct probabilities\n");
	Viterbi(&hmm, T, O, delta, psi, q, &proba);
	fprintf(stdout, "Viterbi  MLE log prob = %E\n", log(proba));
	fprintf(stdout, "Optimal state sequence:\n");
	PrintSequence(stdout, T, q);

	printf("------------------------------------\n");
	printf("Viterbi using log probabilities\n");

	ViterbiLog(&hmm, T, O, delta, psi, q, &logproba); 

	fprintf(stdout, "Viterbi  MLE log prob = %E\n", logproba);
	fprintf(stdout, "Optimal state sequence:\n");
	PrintSequence(stdout, T, q);
	printf("------------------------------------\n");
	printf("The two log probabilites and optimal state sequences\n");
	printf("should identical (within numerical precision). \n");

	free_ivector(q, 1, T);
	free_ivector(O, 1, T);
	free_imatrix(psi, 1, T, 1, hmm.N);
	free_dmatrix(delta, 1, T, 1, hmm.N);
	FreeHMM(&hmm);
}