model.c

来自「General Hidden Markov Model Library 一个通用」· C语言 代码 · 共 2,232 行 · 第 1/5 页

{
# define CUR_PROC "ghmm_dmodel_A_print_transp"
  int i, j;
  int *out_state;

  ARRAY_CALLOC (out_state, mo->N);
  for (i = 0; i < mo->N; i++)
    out_state[i] = 0;

  for (j = 0; j < mo->N; j++) {
    fprintf (file, "%s", tab);
    if (mo->s[0].out_states != 0 && mo->s[0].out_id[out_state[0]] == j) {
      fprintf (file, "%.2f", mo->s[0].out_a[out_state[0]]);
      (out_state[0])++;
    }
    else
      fprintf (file, "0.00");
    for (i = 1; i < mo->N; i++) {
      if (mo->s[i].out_states != 0 && mo->s[i].out_id[out_state[i]] == j) {
        fprintf (file, "%s %.2f", separator, mo->s[i].out_a[out_state[i]]);
        (out_state[i])++;
      }
      else
        fprintf (file, "%s 0.00", separator);
    }
    fprintf (file, "%s\n", ending);
  }
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  m_free (out_state);
  return;
# undef CUR_PROC
}                               /* ghmm_dmodel_A_print_transp */

/*============================================================================*/

void ghmm_dmodel_B_print_transp (FILE * file, ghmm_dmodel * mo, char *tab,
                           char *separator, char *ending)
{
  int i, j;
  for (j = 0; j < mo->M; j++) {
    fprintf (file, "%s", tab);
    fprintf (file, "%.2f", mo->s[0].b[j]);
    for (i = 1; i < mo->N; i++)
      fprintf (file, "%s %.2f", separator, mo->s[i].b[j]);
    fprintf (file, "%s\n", ending);
  }
}                               /* ghmm_dmodel_B_print_transp */

/*============================================================================*/

void ghmm_dmodel_Pi_print_transp (FILE * file, ghmm_dmodel * mo, char *tab, char *ending)
{
  int i;
  for (i = 0; i < mo->N; i++)
    fprintf (file, "%s%.2f%s\n", tab, mo->s[i].pi, ending);
}                               /* ghmm_dmodel_Pi_print_transp */

/*============================================================================*/

void ghmm_dl_print (FILE * file, ghmm_dmodel * mo, char *tab, char *separator,
                        char *ending)
{
  int i;
  fprintf (file, "%s%d", tab, mo->label[0]);
  for (i = 1; i < mo->N; i++)
    fprintf (file, "%s %d", separator, mo->label[i]);
  fprintf (file, "%s\n", ending);
}                               /* ghmm_dl_print */

/*============================================================================*/
void ghmm_dmodel_print (FILE * file, ghmm_dmodel * mo)
{
  fprintf (file, "HMM = {\n\tM = %d;\n\tN = %d;\n", mo->M, mo->N);
  fprintf (file, "\tprior = %.3f;\n", mo->prior);
  fprintf (file, "\tModelType = %d;\n", mo->model_type);
  fprintf (file, "\tA = matrix {\n");
  ghmm_dmodel_A_print (file, mo, "\t", ",", ";");
  fprintf (file, "\t};\n\tB = matrix {\n");
  ghmm_dmodel_B_print (file, mo, "\t", ",", ";");
  fprintf (file, "\t};\n\tPi = vector {\n");
  ghmm_dmodel_Pi_print (file, mo, "\t", ",", ";");
  fprintf (file, "\t};\n\tfix_state = vector {\n");
  ghmm_dmodel_fix_print (file, mo, "\t", ",", ";");
  if (mo->model_type & GHMM_kLabeledStates) {
    fprintf (file, "\t};\n\tlabel_state = vector {\n");
    ghmm_dl_print (file, mo, "\t", ",", ";");
  }
  fprintf (file, "\t};\n};\n\n");
}                               /* ghmm_dmodel_print */

/*============================================================================*/

#ifdef GHMM_OBSOLETE
void ghmm_dmodel_direct_print (FILE * file, model_direct * mo_d, int multip)
{
  int i, j;
  for (i = 0; i < multip; i++) {
    fprintf (file, "HMM = {\n\tM = %d;\n\tN = %d;\n", mo_d->M, mo_d->N);
    fprintf (file, "\tprior = %.3f;\n", mo_d->prior);
    fprintf (file, "\tA = matrix {\n");
    ighmm_cmatrix_print (file, mo_d->A, mo_d->N, mo_d->N, "\t", ",", ";");
    fprintf (file, "\t};\n\tB = matrix {\n");
    ighmm_cmatrix_print (file, mo_d->B, mo_d->N, mo_d->M, "\t", ",", ";");
    fprintf (file, "\t};\n\tPi = vector {\n");
    fprintf (file, "\t%.4f", mo_d->Pi[0]);
    for (j = 1; j < mo_d->N; j++)
      fprintf (file, ", %.4f", mo_d->Pi[j]);
    fprintf (file, ";\n\t};\n");
    fprintf (file, "\tfix_state = vector {\n");
    fprintf (file, "\t%d", mo_d->fix_state[0]);
    for (j = 1; j < mo_d->N; j++)
      fprintf (file, ", %d", mo_d->fix_state[j]);
    fprintf (file, ";\n\t};\n");
    fprintf (file, "};\n\n");
  }
}                               /* ghmm_dmodel_direct_print */

/*============================================================================*/

void ghmm_dmodel_direct_clean (model_direct * mo_d, hmm_check_t * check)
{
#define CUR_PROC "ghmm_dmodel_direct_clean"
  int i;
  if (!mo_d)
    return;
  mo_d->M = mo_d->N = 0;
  mo_d->prior = -1;
  if (mo_d->A) {
    for (i = 0; i < check->r_a; i++)
      m_free (mo_d->A[i]);
    m_free (mo_d->A);
  }
  if (mo_d->B) {
    for (i = 0; i < check->r_b; i++)
      m_free (mo_d->B[i]);
    m_free (mo_d->B);
  }
  if (mo_d->Pi){
    m_free (mo_d->Pi);
  }
  if (mo_d->fix_state){
    m_free (mo_d->fix_state);
  }
  
  mo_d->A = mo_d->B = NULL;
  mo_d->Pi = NULL;
  mo_d->fix_state = NULL;
#undef CUR_PROC
}                               /* ghmm_dmodel_direct_clean */

/*============================================================================*/

int ghmm_dmodel_direct_check_data (model_direct * mo_d, hmm_check_t * check)
{
#define CUR_PROC "ghmm_dmodel_direct_check_data"
  char *str;
  if (check->r_a != mo_d->N || check->c_a != mo_d->N) {
    str = ighmm_mprintf (NULL, 0, "Incompatible dim. A (%d X %d) and N (%d)\n",
                   check->r_a, check->c_a, mo_d->N);
    GHMM_LOG(LCONVERTED, str);
    m_free (str);
    return (-1);
  }
  if (check->r_b != mo_d->N || check->c_b != mo_d->M) {
    str = ighmm_mprintf (NULL, 0,
            "Incompatible dim. B (%d X %d) and N X M (%d X %d)\n",
            check->r_b, check->c_b, mo_d->N, mo_d->M);
    GHMM_LOG(LCONVERTED, str);
    m_free (str);
    return (-1);
  }
  if (check->len_pi != mo_d->N) {
    str = ighmm_mprintf (NULL, 0, "Incompatible dim. Pi (%d) and N (%d)\n",
                   check->len_pi, mo_d->N);
    GHMM_LOG(LCONVERTED, str);
    m_free (str);
    return (-1);
  }
  if (check->len_fix != mo_d->N) {
    str = ighmm_mprintf (NULL, 0, "Incompatible dim. fix_state (%d) and N (%d)\n",
                   check->len_fix, mo_d->N);
    GHMM_LOG(LCONVERTED, str);
    m_free (str);
    return (-1);
  }

  return 0;
#undef CUR_PROC
}                               /* ghmm_dmodel_direct_check_data */
#endif /* GHMM_OBSOLETE */


/*============================================================================*/
/* XXX symmetric not implemented yet */
double ghmm_dmodel_prob_distance (ghmm_dmodel * m0, ghmm_dmodel * m, int maxT, int symmetric,
                            int verbose)
{
#define CUR_PROC "ghmm_dmodel_prob_distance"

#define STEPS 40

  double p0, p;
  double d = 0.0;
  double *d1;
  ghmm_dseq *seq0 = NULL;
  ghmm_dseq *tmp = NULL;
  ghmm_dmodel *mo1, *mo2;
  int i, t, a, k;
  int true_len;
  int true_number;
  int left_to_right = 0;
  long total, index;
  int step_width = 0;
  int steps = 1;

  /* printf("***  model_prob_distance:\n"); */

  if (verbose) {                /* If we are doing it verbosely we want to have 40 steps */
    step_width = maxT / 40;
    steps = STEPS;
  }
  else                          /* else just one */
    step_width = maxT;

  ARRAY_CALLOC (d1, steps);

  mo1 = m0;
  mo2 = m;

  for (k = 0; k < 2; k++) {     /* Two passes for the symmetric case */

    /* seed = 0 -> no reseeding. Call  ghmm_rng_timeseed(RNG) externally */
    seq0 = ghmm_dmodel_generate_sequences (mo1, 0, maxT + 1, 1, maxT + 1);



    if (seq0 == NULL) {
      GHMM_LOG(LCONVERTED, " generate_sequences failed !");
      goto STOP;
    }

    if (seq0->seq_len[0] < maxT) {      /* There is an absorbing state */

      /* NOTA BENE: Assumpting the model delivers an explicit end state, 
         the condition of a fix initial state is removed. */

      /* For now check that Pi puts all weight on state */
      /*
         t = 0;
         for (i = 0; i < mo1->N; i++) {
         if (mo1->s[i].pi > 0.001)
         t++;
         }    
         if (t > 1) {
         GHMM_LOG(LCONVERTED, "ERROR: No proper left-to-right model. Multiple start states");
         goto STOP;
         } */

      left_to_right = 1;
      total = seq0->seq_len[0];

      while (total <= maxT) {

        /* create a additional sequences at once */
        a = (maxT - total) / (total / seq0->seq_number) + 1;
        /* printf("total=%d generating %d", total, a); */
        tmp = ghmm_dmodel_generate_sequences (mo1, 0, 0, a, a);
        if (tmp == NULL) {
          GHMM_LOG(LCONVERTED, " generate_sequences failed !");
          goto STOP;
        }
        ghmm_dseq_free (&tmp);
        ghmm_dseq_add (seq0, tmp);

        total = 0;
        for (i = 0; i < seq0->seq_number; i++)
          total += seq0->seq_len[i];
      }
    }

    if (left_to_right) {

      for (t = step_width, i = 0; t <= maxT; t += step_width, i++) {

        index = 0;
        total = seq0->seq_len[0];

        /* Determine how many sequences we need to get a total of t
           and adjust length of last sequence to obtain total of 
           exactly t */

        while (total < t) {
          index++;
          total += seq0->seq_len[index];
        }

        true_len = seq0->seq_len[index];
        true_number = seq0->seq_number;

        if ((total - t) > 0)
          seq0->seq_len[index] = total - t;
        seq0->seq_number = index;

        p0 = ghmm_dmodel_likelihood (mo1, seq0);
        if (p0 == +1 || p0 == -1) {     /* error! */
          GHMM_LOG(LCONVERTED, "problem: ghmm_dmodel_likelihood failed !");
          goto STOP;
        }
        p = ghmm_dmodel_likelihood (mo2, seq0);
        if (p == +1 || p == -1) {       /* what shall we do now? */
          GHMM_LOG(LCONVERTED, "problem: ghmm_dmodel_likelihood failed !");
          goto STOP;
        }

        d = 1.0 / t * (p0 - p);

        if (symmetric) {
          if (k == 0)
            /* save d */
            d1[i] = d;
          else {
            /* calculate d */
            d = 0.5 * (d1[i] + d);
          }
        }

        if (verbose && (!symmetric || k == 1))
          printf ("%d\t%f\t%f\t%f\n", t, p0, p, d);

        seq0->seq_len[index] = true_len;
        seq0->seq_number = true_number;
      }
    }

    else {

      true_len = seq0->seq_len[0];

      for (t = step_width, i = 0; t <= maxT; t += step_width, i++) {
        seq0->seq_len[0] = t;

        p0 = ghmm_dmodel_likelihood (mo1, seq0);
        /* printf("   P(O|m1) = %f\n",p0); */
        if (p0 == +1) {         /* error! */
          GHMM_LOG(LCONVERTED, "seq0 can't be build from mo1!");
          goto STOP;
        }
        p = ghmm_dmodel_likelihood (mo2, seq0);
        /* printf("   P(O|m2) = %f\n",p); */
        if (p == +1) {          /* what shall we do now? */
          GHMM_LOG(LCONVERTED, "problem: seq0 can't be build from mo2!");
          goto STOP;
        }

        d = (1.0 / t) * (p0 - p);

        if (symmetric) {
          if (k == 0)
            /* save d */
            d1[i] = d;
          else {
            /* calculate d */
            d = 0.5 * (d1[i] + d);
          }
        }

        if (verbose && (!symmetric || k == 1))
          printf ("%d\t%f\t%f\t%f\n", t, p0, p, d);

      }
      seq0->seq_len[0] = true_len;
    }

    if (symmetric) {
      ghmm_dseq_free (&seq0);
      mo1 = m;
      mo2 = m0;
    }
    else
      break;

  }                             /* k = 1,2 */

  ghmm_dseq_free (&seq0);
  free (d1);
  return d;

STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ghmm_dseq_free (&seq0);
  free (d1);
  return (0.0);
#undef CUR_PROC
}


/*============================================================================*/

void ghmm_dstate_clean (ghmm_dstate * my_state)
{
#define CUR_PROC "ghmm_dstate_clean"

  mes_check_ptr (my_state, return;);

  if (my_state->b){
    m_free (my_state->b);
  }
  if (my_state->out_id){
    m_free (my_state->out_id);
  }
  if (my_state->in_id){
    m_free (my_state->in_id);
  }
  if (my_state->out_a){
    m_free (my_state->out_a);
  }
  if (my_state->in_a){
    m_free (my_state->in_a);
  }
  my_state->pi = 0;
  my_state->b = NULL;
  my_state->out_id = NULL;
  my_state->in_id = NULL;
  my_state->out_a = NULL;
  my_state->in_a = NULL;
  my_state->out_states = 0;
  my_state->in_states = 0;
  my_state->fix = 0;

#undef CUR_PROC
}                               /* ghmm_dstate_clean */



/*==========================  Labeled HMMs  ================================*/

ghmm_dseq *ghmm_dmodel_label_generate_sequences(
    ghmm_dmodel * mo, int seed, int global_len, long seq_number, int Tmax)
{
#define CUR_PROC "ghmm_dmodel_label_generate_sequences"

  ghmm_dseq *sq = NULL;
  int state;
  int j, m, j_id;
  double p, sum, max_sum;
  int len = global_len;
  int n = 0;
  int pos, label_pos;