ghmmunittests.c

General Hidden Markov Model Library 一个通用的隐马尔科夫模型的C代码库

#ifdef HAVE_CONFIG_H
#include "../config.h"
#endif /* HAVE_CONFIG_H */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <ghmm/smodel.h>
#include <check.h>

#include <ghmm/psequence.h>
#include <ghmm/pmodel.h>
#include <ghmm/pviterbi.h>
#include <ghmm/pviterbi_propagate.h>
#include "testpair.h"

/*******************************************************************************

  Constants for all test cases.

*******************************************************************************/

#define kSMODEL_SAMPLE_SMO_FILE "sample.smo"


/****** smodel_suite *********************************************************/
START_TEST (smodel_read_free) /* Add smodel_read_free in  smodel_suite() below */
{
  /* Test whether we can read an ghmm_cmodel, iterate over it, and free it again */

  ghmm_cmodel **model_array;
  int model_counter, result;
  char *inFileName = kSMODEL_SAMPLE_SMO_FILE;
  model_array = ghmm_c_read(inFileName, &model_counter);

  fail_unless(model_counter == 1, 
	      "Read %d number of models instead of one", model_counter);

  fail_unless(model_array != NULL,
	      "ghmm_c_read returned null pointer");

  result = ghmm_c_free(model_array);  
  
  fail_unless(result == NULL,
	      "ghmm_c_free failed");  
}
END_TEST







Suite *smodel_suite(void)
{
  Suite *s = suite_create("smodel");
  TCase *tc_core = tcase_create("Core");

  suite_add_tcase (s, tc_core); 
  
  tcase_add_test(tc_core, smodel_read_free);

  return s;
}
/****** end of smodel_suite ************************************************/


/****** PairHMM suite ******************************************************/

START_TEST (mysequence_test) /* Add mysequence_test in  pmodel_suite() below */
{
  /* check if init, access to fields and free work */

  mysequence * x;
  int length = 100;
  int i, j;
  int discrete = 1;
  int continuous = 1;
  /* check the initialisation */
  x = init_mysequence(length, discrete, continuous);
  fail_unless(x->length == length, "Init sequence length");
  fail_unless(x->number_of_alphabets == discrete, "Init sequence discrete");
  fail_unless(x->number_of_d_seqs == continuous, "Init sequence continous");
  /* check the field access */
  for (i=0; i<length; i++) {
    for (j=0; j<discrete; j++) {
      x->seq[j][i] = i % 4;
      fail_unless(get_char_mysequence(x, j, i) == i % 4, "get_char_mysequence returns %i instead of %i at position %i", get_char_mysequence(x, j, i), i % 4, i);
    }
    for (j=0; j<continuous; j++) {
      x->d_value[j][i] = 1.5 * i;
      fail_unless(get_double_mysequence(x, j, i) == 1.5 * i, "get_double_mysequence returns %i instead of %i at position %i", get_double_mysequence(x, j, i), 1.5 * i, i);
    }
  }
  /* check the free */
  int result = free_mysequence(x, discrete, continuous);
  fail_unless(result == 0, "Failure in free_mysequence");
}
END_TEST

START_TEST (pviterbi_path_test) /* Add ghmm_dp_viterbi_test in  pmodel_suite() below */
{
  /* test whether the computed viterbi path correctly aligns the two sequences */
  ghmm_dpmodel * mo = get_test_pmodel();

  /* first get two complex sequences */
  mysequence * x;
  mysequence * y;
  int length = 100;
  int i, j;
  x = init_mysequence(length, 1, 0);
  for (i=0; i<length; i++)
    x->seq[0][i] = i % 4;
  y = init_mysequence(length, 1, 0);
  for (i=0; i<length; i++)
    y->seq[0][i] = i % 4;

  /* introduce some mismatches */  
  x->seq[0][3] = 2;
  y->seq[0][4] = 2;

  /* run the viterbi algorithm */
  double logp = 1.0;
  length = 0;
  int * path;
  path = ghmm_dp_viterbi(mo, x, y, &logp, &length);

  /* check if logp is the expected log p */
  double exp_logp = -168.239906;
  fail_unless(abs(logp - exp_logp) < PROP_EPS, "Log p of the path should be %f instead of %f", exp_logp, logp);

  /* double check with the logp function */
  double check_logp = ghmm_dp_viterbi_logp(mo, x, y, path, length);
  fail_unless(abs(check_logp - exp_logp) < PROP_EPS, "Log p computed by the ghmm_dp_viterbi_logp function of the path should be %f instead of %f", exp_logp, check_logp);
  fail_unless(abs(check_logp - logp) < PROP_EPS, "Log p should be equal (ghmm_dp_viterbi vs pviterbi_logp: %f != %f", logp, check_logp);

  /* the begin of the path has to be: 0,0,0,1,0,2 because of the mismatches */
  fail_unless(path[3] == 1 && path[4] == 0 && path[5] == 2, "Alignment impossible!!!");

  /* free everything */
  m_free(path);
  free_mysequence(x, 1, 0);
  free_mysequence(y, 1, 0);
  ghmm_dp_free(mo);
}
END_TEST

START_TEST (pviterbi_propagate_test) /* Add ghmm_dp_viterbi_test in  pmodel_suite() below */
{
  /* test whether the computed viterbi path correctly aligns the two sequences */
  ghmm_dpmodel * mo = get_test_pmodel();

  /* first get two complex sequences */
  mysequence * x;
  mysequence * y;
  int length = 100;
  int i, j;
  x = init_mysequence(length, 1, 0);
  for (i=0; i<length; i++)
    x->seq[0][i] = i % 4;
  y = init_mysequence(length, 1, 0);
  for (i=0; i<length; i++)
    y->seq[0][i] = i % 4;

  /* introduce some mismatches */  
  x->seq[0][3] = 2;
  y->seq[0][4] = 2;

  /* run the viterbi algorithm */
  double naive_logp = 1.0;
  int naive_length = 0;
  int * naive_path;
  naive_path = ghmm_dp_viterbi_propagate(mo, x, y, &naive_logp, &naive_length, 500);

  /* run the linear memory version */
  double logp = 1.0;
  length = 0;
  int * path;
  path = ghmm_dp_viterbi_propagate(mo, x, y, &logp, &length, 500);

  /* check if the lengths of the paths are equal */
  fail_unless(length == naive_length, "Path length of the linear verion (%i) differs from naive implementation (%i)", length, naive_length);

  /* check if logp is the expected log p */
  double exp_logp = -168.239906;
  fail_unless(abs(logp - exp_logp) < PROP_EPS, "Log p of the path should be %f instead of %f", exp_logp, logp);

  /* double check with the logp function */
  double check_logp = ghmm_dp_viterbi_logp(mo, x, y, path, length);
  fail_unless(abs(check_logp - exp_logp) < PROP_EPS, "Log p computed by the ghmm_dp_viterbi_logp function of the path should be %f instead of %f", exp_logp, check_logp);
  fail_unless(abs(check_logp - logp) < PROP_EPS, "Log p should be equal (ghmm_dp_viterbi vs pviterbi_logp: %f != %f", logp, check_logp);

  /* the begin of the path has to be: 0,0,0,1,0,2 because of the mismatches */
  fail_unless(path[3] == 1 && path[4] == 0 && path[5] == 2, "Alignment impossible!!!");

  /* do a state by state comparison of the paths */
  for (i=0; i<length; i++) {
    fail_unless(path[i] == naive_path[i], "States differ at position %i: naive %i vs linear %i", i, naive_path[i], path[i]);
  }

  /* free everything */
  m_free(path);
  m_free(naive_path);
  free_mysequence(x, 1, 0);
  free_mysequence(y, 1, 0);
  ghmm_dp_free(mo);
}
END_TEST

Suite *pmodel_suite(void)
{
  Suite *s = suite_create("pmodel");
  TCase *tc_core = tcase_create("Core");

  suite_add_tcase (s, tc_core); 
  
  tcase_add_test(tc_core, mysequence_test);
  tcase_add_test(tc_core, pviterbi_path_test);
  tcase_add_test(tc_core, pviterbi_propagate_test);

  return s;
}


/****** end of PairHMM suite ***********************************************/


int main(void)
{
  int nf;
  //Suite *s = smodel_suite();
  Suite *s = pmodel_suite();
  SRunner *sr = srunner_create(s);

  ghmm_rng_init();  /* Important! initialise rng  XXX UNTESTED*/

  srunner_run_all(sr, CK_NORMAL);
  nf = srunner_ntests_failed(sr);
  srunner_free(sr);
  return (nf == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
}