PageRenderTime 32ms CodeModel.GetById 24ms app.highlight 6ms RepoModel.GetById 1ms app.codeStats 0ms

/ghmm-0.7.0/tests/shmm_viterbi_test.c

http://github.com/jcnossen/hmmgenefinder
C | 77 lines | 49 code | 16 blank | 12 comment | 12 complexity | b91159309cd8e2c2e3e9b8ea6b1390e1 MD5 | raw file
Possible License(s): GPL-2.0
 1/*******************************************************************************
 2  author       : Bernd Wichern
 3  filename     : /zpr/bspk/src/hmm/ghmm/tests/shmm_viterbi_test.c
 4  created      : TIME: 12:52:37     DATE: Tue 05. June 2001
 5  last-modified: TIME: 15:33:06     DATE: Tue 05. June 2001
 6  $Id: shmm_viterbi_test.c 201 2002-02-22 23:48:44Z achim $
 7*******************************************************************************/
 8
 9#include <string.h>
10#include <ghmm/smodel.h>
11#include <ghmm/sequence.h>
12#include <ghmm/sviterbi.h>
13#include <ghmm/mes.h>
14
15static int viterbi_test(char* argv[]);
16
17/*============================================================================*/
18
19static int viterbi_test(char* argv[]) {
20#define CUR_PROC "viterbi_test"
21  FILE *outfile = NULL;
22  sequence_d_t **sqd = NULL;
23  smodel** smo = NULL;
24  char outfilename[256];
25  int smo_number, sqd_number, *state_seq = NULL, model, t, cnt = 0;
26  long i, j;
27  double log_p;
28
29
30  /* read array of models */
31  smo = smodel_read(argv[1], &smo_number);
32  if (!(smo)) {mes_proc(); goto STOP;}
33  
34  /* read sequences */
35  sqd = sequence_d_read(argv[2], &sqd_number);
36  if (!sqd) {mes_proc(); goto STOP;}
37  
38  strcpy(outfilename, argv[3]);
39  if(!(outfile = mes_fopen(outfilename, "wt"))) {mes_proc(); goto STOP;}
40
41  /* calculate viterbi path for every possible sequence-model combination */
42  for (model = 0; model < smo_number; model++) {
43    for (i = 0;  i < sqd_number; i++) {
44      for (j = 0; j < sqd[i]->seq_number; j++) {
45	cnt++;
46	state_seq = sviterbi(smo[model], sqd[i]->seq[j], sqd[i]->seq_len[j],
47			     &log_p);
48	if (state_seq == NULL) {mes_proc(); goto STOP;}
49	fprintf(outfile, "%d %ld (Seq.ID %d): logp %.4f\t", model, i, 
50		(int) sqd[i]->seq_id[j], log_p);
51	for (t = 0; t < sqd[i]->seq_len[j]; t++)
52	  fprintf(outfile, "%2d ", state_seq[t]);
53	fprintf(outfile, "\n");
54	m_free(state_seq);
55	if (!cnt%100) printf("%d\n", cnt);
56      }
57    }
58  }
59
60
61STOP:
62  return -1;
63# undef CUR_PROC
64}
65
66/*============================================================================*/
67
68int main(int argc, char* argv[]) {
69
70  if (argc != 4) {
71    printf("Usage: shmm_viterbi_test <Model File> <Sequence File> <Output File> \n");
72    exit(1);
73  }
74
75  return viterbi_test(argv);
76}
77