sgenerate.c

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

      ARRAY_REALLOC (sq->seq[n], t);
    sq->seq_len[n] = t;

  }                             /* for n .. < seq_number */

  ighmm_cmatrix_free (&alpha, max_short_len);
  m_free (scale);
  return sq;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ighmm_cmatrix_free (&alpha, max_short_len);
  ghmm_cseq_free (&sq);
  return (NULL);
# undef CUR_PROC
}                               /* ghmm_sgenerate_extensions */

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

double *ghmm_sgenerate_single_ext (ghmm_cmodel * smo, double *O, const int len,
                              int *new_len, double **alpha,
                              sgeneration_mode_t mode)
{
# define CUR_PROC "ghmm_sgenerate_single_ext"
  int i, j, m, t, class=0, up = 0;
  double *new_O = NULL, *scale = NULL, *initial_distribution = NULL;
  double log_p, sum, p;
  int k;
  char * str;

  /* TEMP */
  if (mode == all_viterbi || mode == viterbi_viterbi || mode == viterbi_all) {
    GHMM_LOG(LCONVERTED, "Error: mode not implemented yet\n");
    goto STOP;
  }
  if (len <= 0) {
    GHMM_LOG(LCONVERTED, "Error: sequence with zero or negativ length\n");
    goto STOP;
  }
  ARRAY_CALLOC (new_O, (int) GHMM_MAX_SEQ_LEN);
  ARRAY_CALLOC (scale, len);
  ARRAY_CALLOC (initial_distribution, smo->N);
  ghmm_cseq_copy (new_O, O, len);
  *new_len = len;
  /* Initial Distribution ???
     Pi(i) = alpha_t(i)/P(O|lambda) */
  if (mode == all_all || mode == all_viterbi) {
    if (ghmm_cmodel_forward (smo, O, len, NULL /* ?? */ , alpha, scale, &log_p)) {
      GHMM_LOG(LCONVERTED, "error from sfoba_forward, unable to extend\n");
      ARRAY_REALLOC (new_O, *new_len);
      return new_O;
    }
    sum = 0.0;


    for (i = 0; i < smo->N; i++) {
      /* alpha is scaled! */
      initial_distribution[i] = alpha[len - 1][i];
      sum += initial_distribution[i];
    }
    /* nicht ok.? scale to one? */
    for (i = 0; i < smo->N; i++) {
      initial_distribution[i] /= sum;
    }
  }

  p = GHMM_RNG_UNIFORM (RNG);
  sum = 0.0;
  for (i = 0; i < smo->N; i++) {
    sum += initial_distribution[i];
    if (sum >= p)
      break;
  }
  /* error due to incorrect normalization ?? */
  if (i == smo->N) {
    i--;
    while (i > 0 && initial_distribution[i] == 0.0)
      i--;
  }
  /* TEST */
  /* Already at the beginning in an end state? How is that possible? */
  if (smo->s[i].out_states == 0) {
    printf ("Beginn: Endzustand, State %d\n", i);
    for (k = 0; k < len; k++)
      printf ("%.2f ", O[k]);
    printf ("\n");
  }
  /* End Test */

  for (t = 0; t < len; t++)
    if (smo->cos == 1) {
      class = 0;
    }
    else {
      if (!smo->class_change->get_class) {
        printf ("ERROR: get_class not initialized\n");
        goto STOP;
      }
      /*printf("1: cos = %d, k = %d, t = %d\n",smo->cos,smo->class_change->k,t);*/
      class = smo->class_change->get_class (smo, O, 0, t);      /*XXX No sequence number */
    }

  t = len;
  while (t < (int) GHMM_MAX_SEQ_LEN) {
    if (smo->s[i].out_states == 0)
      /* reached final state, exit while loop */
      break;
    p = GHMM_RNG_UNIFORM (RNG);
    sum = 0.0;
    for (j = 0; j < smo->s[i].out_states; j++) {
      sum += smo->s[i].out_a[class][j];
      if (sum >= p)
        break;
    }
    /* error due to incorrect normalization ?? */
    if (j == smo->s[i].out_states) {
      j--;
      while (j > 0 && smo->s[i].out_a[class][j] == 0.0)
        j--;
    }
    if (sum == 0.0) {
      /* Test: If an "empty" class, try the neighbour class;
         first, sweep down to zero, if still no success sweep up 
         to COs - 1. If still no success --> discard the sequence.
       */
      if (class > 0 && up == 0) {
        class--;
        continue;
      }
      else if (class < smo->cos - 1) {
        class++;
        up = 1;
        continue;
      }
      else {
        str = ighmm_mprintf (NULL, 0, "unable to extend seq (all osc empty)\n");
        GHMM_LOG(LCONVERTED, str);
        m_free (str);
        goto STOP;
      }
    }                           /* sum == 0 */

    /* new state */
    i = smo->s[i].out_id[j];

    if (smo->M == 1)
      m = 0;
    else {
      p = GHMM_RNG_UNIFORM (RNG);
      sum = 0.0;
      for (m = 0; m < smo->M; m++) {
        sum += smo->s[i].c[m];
        if (sum >= p)
          break;
      }
      if (m == smo->M) {
        m--;
        while (m > 0 && smo->s[i].c[m] == 0.0)
          m--;
      }
    }
    /* Output in state i, komp. m */
    /* random variable from density function */
    new_O[t] = ghmm_cmodel_get_random_var (smo, i, m);

    if (smo->cos == 1) {
      class = 0;
    }
    else {
      if (!smo->class_change->get_class) {
        printf ("ERROR: get_class not initialized\n");
        goto STOP;
      }
      /*printf("1: cos = %d, k = %d, t = %d\n",smo->cos,smo->class_change->k,t);*/
      class = smo->class_change->get_class (smo, new_O, 0, t);  /* XXX sequence number ? */
    }

    t++;
    up = 0;
  }                             /* while (t < MAX_SEQ_LEN) */
  if (t < (int) GHMM_MAX_SEQ_LEN)
    ARRAY_REALLOC (new_O, t);

  *new_len = t;


  m_free (scale);
  m_free (initial_distribution);

  return new_O;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  m_free (new_O);
  m_free (scale);
  m_free (initial_distribution);
  return NULL;
# undef CUR_PROC
}                               /* ghmm_sgenerate_single_ext */


/* generate a single next value bases on a trained model and on a seq und
   to length "len". Use the most prob. state given the seq as an initial state
   and determin the next state und the symbol with the RNG.

   If this function is called for one O successive "len", then it's possible
   to optimize some more.
*/

double ghmm_sgenerate_next_value (ghmm_cmodel * smo, double *O, const int len)
{
# define CUR_PROC "ghmm_sgenerate_next_value"
  double **alpha = NULL;
  double res = -1.0, sum, p;
  double *scale = NULL, log_p, max_val = -1000000;
  int i, j, m, init_state = -1;

  if (smo->cos > 1) {
    GHMM_LOG(LCONVERTED, "ghmm_sgenerate_next_value only for COS == 1\n");
    goto STOP;
  }

  alpha = ighmm_cmatrix_alloc (len, smo->N);
  if (!alpha) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  ARRAY_CALLOC (scale, len);
  if (ghmm_cmodel_forward (smo, O, len, NULL /* ?? */ , alpha, scale, &log_p)) {
    GHMM_LOG(LCONVERTED, "error from sfoba_forward\n");
    goto STOP;
  }

  /* find inititial state */
  sum = 0.0;
  for (i = 0; i < smo->N; i++)
    sum += alpha[len - 1][i];
  if (sum < 0.9 || sum > 1.1) {
    printf ("Error sum = %.4f (!= 1)\n", sum);
    goto STOP;
  }
  /* max state */
  for (i = 0; i < smo->N; i++) {
    if (alpha[len - 1][i] > max_val) {
      init_state = i;
      max_val = alpha[len - 1][i];
    }
  }

  /* random state */
  /*
     p = GHMM_RNG_UNIFORM(RNG);
     sum = 0.0;
     for (i = 0; i < smo->N; i++) {
     sum += alpha[len - 1][i];
     if (sum >= p)
     break;    
     }   
     if (i == smo->N) {
     i--;
     while (i > 0 && alpha[len - 1][i] == 0.0) i--;
     }  
     init_state = i;
   */


  if (init_state == -1 || smo->s[init_state].out_states == 0)
    goto STOP;

  p = GHMM_RNG_UNIFORM (RNG);
  sum = 0.0;
  for (j = 0; j < smo->s[init_state].out_states; j++) {
    sum += smo->s[init_state].out_a[0][j];
    if (sum >= p)
      break;
  }
  /* error due to incorrect normalization ?? */
  if (j == smo->s[init_state].out_states) {
    j--;
    while (j > 0 && smo->s[init_state].out_a[0][j] == 0.0)
      j--;
  }

  /* new state */
  i = smo->s[init_state].out_id[j];

  if (smo->M == 1)
    m = 0;
  else {
    p = GHMM_RNG_UNIFORM (RNG);
    sum = 0.0;
    for (m = 0; m < smo->M; m++) {
      sum += smo->s[i].c[m];
      if (sum >= p)
        break;
    }
    if (m == smo->M) {
      m--;
      while (m > 0 && smo->s[i].c[m] == 0.0)
        m--;
    }
  }
  /* Output in state i, komp. m */
  /* random variable from density function */
  res = ghmm_cmodel_get_random_var (smo, i, m);

STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ighmm_cmatrix_free (&alpha, len);
  m_free (scale);
  return res;
# undef CUR_PROC
}