sdmodel.c

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

  double sum;
  int *__silents;

  ARRAY_CALLOC (__silents, mo->N);

  for (i = 0; i < mo->N; i++) {
    for (m = 0, sum = 0; m < mo->M; m++) {
      sum += mo->s[i].b[m];
    }
    if (sum < __EPS) {
      __silents[i] = 1;
      nSilentStates++;
    }
    else
      __silents[i] = 0;
  }

  if (nSilentStates) {
    mo->model_type = GHMM_kSilentStates;
    mo->silent = __silents;
  }
  else {
    mo->model_type = GHMM_kNotSpecified;
    mo->silent = NULL;
    m_free (__silents);
  }
  return (nSilentStates);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return (0);
#undef CUR_PROC
}

/*============================================================================*/
/*
 * Before calling generate sequence, the global random number generator must
 * be initialized by calling ghmm_rng_init().
 *
 */
/*returns the extended sequence struct with state matrix*/
ghmm_dseq *ghmm_dsmodel_generate_sequences (ghmm_dsmodel * mo, int seed,
                                        int global_len, long seq_number,
                                        int Tmax)
{
# define CUR_PROC "ghmm_dsmodel_generate_sequences"

  /* An end state is characterized by not having out-going transition. */
  unsigned long tm;             /* Time seed */
  ghmm_dseq *sq = NULL;
  int state, n, i, j, m, reject_os, reject_tmax, badseq, trans_class;
  double p, sum;
  int len = global_len, up = 0, stillbadseq = 0, reject_os_tmp = 0;
  int obsLength = 0;
  int silent_len = 0, badSilentStates = 0;
  int lastStateSilent = 0;
  int matchcount = 0;

  sq = ghmm_dseq_calloc (seq_number);

  if (!sq) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  if (len <= 0)
    /* A specific length of the sequences isn't given. As a model should have
       an end state, the konstant MAX_SEQ_LEN is used. */
    len = (int) GHMM_MAX_SEQ_LEN;

  if (seed > 0) {
    GHMM_RNG_SET (RNG, seed);
  }
  else {
    tm = rand ();
    GHMM_RNG_SET (RNG, tm);
    fprintf (stderr, "# using rng '%s' seed=%ld\n", GHMM_RNG_NAME (RNG), tm);
  }

  n = 0;
  reject_os = reject_tmax = 0;

  while (n < seq_number) {
    /* Test: A new seed for each sequence */
    /*   ghmm_rng_timeseed(RNG); */
    stillbadseq = badseq = 0;
    matchcount = 0;
    ARRAY_CALLOC (sq->seq[n], len);
    ARRAY_CALLOC (sq->states[n], len);

    /* Get a random initial state i */
    p = GHMM_RNG_UNIFORM (RNG);
    sum = 0.0;
    for (i = 0; i < mo->N; i++) {
      sum += mo->s[i].pi;
      if (sum >= p)
        break;
    }
    /* assert( !mo->silent[i] ); */

    if (mo->model_type & GHMM_kSilentStates) {

      if (!mo->silent[i]) {     /* fDte emits */
        lastStateSilent = 0;
        silent_len = 0;
        /* Get a random initial output m */
        p = GHMM_RNG_UNIFORM (RNG);
        sum = 0.0;
        for (m = 0; m < mo->M; m++) {
          sum += mo->s[i].b[m];
          if (sum >= p)
            break;
        }
        sq->seq[n][0] = m;
        sq->states[n][0] = i;
        if (mo->s[i].countme) {
          matchcount++;
        }
        state = 1;
      }
      else {                    /* silent state: we do nothing, no output */
        state = 0;
        lastStateSilent = 1;
        silent_len = 1;
      }
    }
    else {
      /* Get a random initial output m */
      p = GHMM_RNG_UNIFORM (RNG);
      sum = 0.0;
      for (m = 0; m < mo->M; m++) {
        sum += mo->s[i].b[m];
        if (sum >= p)
          break;
      }
      sq->seq[n][0] = m;
      sq->states[n][0] = i;
      if (mo->s[i].countme) {
        matchcount++;
      }
      state = 1;
    }

    obsLength = 0;

    /* class NOT dependent on observable but on path!!!! */
    while (state < len && obsLength < Tmax) {
      /* Get a new state */

      if (mo->cos > 1)
        trans_class = mo->get_class (sq->seq[n], state);
      else
        trans_class = 0;
      p = GHMM_RNG_UNIFORM (RNG);
      sum = 0.0;
      for (j = 0; j < mo->s[i].out_states; j++) {
        sum += mo->s[i].out_a[trans_class][j];
        if (sum >= p)
          break;
      }

      if (sum == 0.0) {
        if (mo->s[i].out_states > 0) {
          /* Repudiate the sequence, if all smo->s[i].out_a[class][.] == 0,
             that is, class "class" isn't used in the original data:
             go out of the while-loop, n should not be counted. */
          /* printf("Zustand %d, class %d, len %d out_states %d \n", i, class,
             state, smo->s[i].out_states); */
          badseq = 1;
          /* break; */

          /* Try: If the class is "empty", try the neighbour class;
             first, sweep down to zero; if still no success, sweep up to
             COS - 1. If still no success --> Repudiate the sequence. */
          if (trans_class > 0 && up == 0) {
            /*trans_class--;*/
            continue;
          }
          else {
            if (trans_class < mo->cos - 1) {
              /*trans_class++;          // as it is not dependent on time!*/
              up = 1;
              continue;
            }
            else {
              stillbadseq = 1;
              break;
            }
          }
        }
        else {
          /* An end state is reached, get out of the while-loop */
          break;
        }
      }

      i = mo->s[i].out_id[j];

      if ((mo->model_type & GHMM_kSilentStates) && mo->silent[i]) {   /* Get a silent state i */
        silent_len++;
        if (silent_len >= Tmax) {
          badSilentStates = 1;
          break;                /* reject this sequence */
        }
        else {
          badSilentStates = 0;
          lastStateSilent = 1;
        }
      }
      else {
        lastStateSilent = 0;
        silent_len = 0;
        /* Get a random output m from state i */
        p = GHMM_RNG_UNIFORM (RNG);
        sum = 0.0;
        for (m = 0; m < mo->M; m++) {
          sum += mo->s[i].b[m];
          if (sum >= p)
            break;
        }
        sq->seq[n][state] = m;
        sq->states[n][state] = i;
        if (mo->s[i].countme) {
          matchcount++;
        }
        state++;
      }

      /* Decide the class for the next step */
      /*trans_class = mo->get_class(sq->seq[n],state);*/
      up = 0;
      obsLength++;
    }                           /* while (state < len) , global_len depends on the data */

    if (badseq) {
      reject_os_tmp++;
    }

    if (stillbadseq) {
      reject_os++;
      m_free (sq->seq[n]);
      m_free (sq->states[n]);
      /*  printf("cl %d, s %d, %d\n", class, i, n); */
    }
    else {
      if (badSilentStates) {
        reject_tmax++;
        m_free (sq->seq[n]);
        m_free (sq->states[n]);
      }
      else {
        if (obsLength > Tmax) {
          reject_tmax++;
          m_free (sq->seq[n]);
          m_free (sq->seq[n]);
        }
        else {
          if (state < len) {
            ARRAY_REALLOC (sq->seq[n], state);
            ARRAY_REALLOC (sq->states[n], state);
          }
          sq->seq_len[n] = state;
          /* ighmm_cvector_print(stdout, sq->seq[n], sq->seq_len[n]," "," ",""); */
          n++;
        }
      }
    }


    /*    printf("reject_os %d, reject_tmax %d\n", reject_os, reject_tmax); */
    if (reject_os > 10000) {
      GHMM_LOG(LCONVERTED, "Reached max. no. of rejections\n");
      break;
    }
    if (!(n % 1000))
      printf ("%d Seqs. generated\n", n);
  }                             /* n-loop, while (n < seq_number) */


  if (reject_os > 0)
    printf ("%d sequences rejected (os)!\n", reject_os);
  if (reject_os_tmp > 0)
    printf ("%d sequences changed class\n", reject_os_tmp - reject_os);
  if (reject_tmax > 0)
    printf ("%d sequences rejected (Tmax, %d)!\n", reject_tmax, Tmax);


  return (sq);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */

  ghmm_dseq_free (&sq);
  return (NULL);
# undef CUR_PROC
}                               /* data */

/*=======================================================================
generate sequences by calling generate_sequences_ext
========================================*/

/*
ghmm_dseq *ghmm_dsmodel_generate_sequences(ghmm_dsmodel* mo, int seed, int global_len,
				     long seq_number, int Tmax) {
  sequence_ext_t *sq = sdmodel_generate_sequences_ext(mo,seed,global_len,seq_number,Tmax);
  printf("jaja\n\n");
  return sequence_unext(&sq);
}
*/
/*============================================================================*/
/* Some outputs */
/*============================================================================*/

void ghmm_dsmodel_states_print (FILE * file, ghmm_dsmodel * mo)
{
  int i, j;
  fprintf (file, "Modelparameters: \n M = %d \t N = %d\n", mo->M, mo->N);
  for (i = 0; i < mo->N; i++) {
    fprintf (file,
             "\nState %d \n PI = %.3f \n out_states = %d \n in_states = %d \n",
             i, mo->s[i].pi, mo->s[i].out_states, mo->s[i].in_states);
    fprintf (file, " Output probability:\t");
    for (j = 0; j < mo->M; j++)
      fprintf (file, "%.3f \t", mo->s[i].b[j]);
    fprintf (file, "\n Transition probability \n");
    fprintf (file, "  Out states (Id, a):\t");
    for (j = 0; j < mo->s[i].out_states; j++)
      fprintf (file, "FIXME: out_a is a matrix"/*(%d, %.3f) \t", mo->s[i].out_id[j], mo->s[i].out_a[j]*/);
    fprintf (file, "\n");
    fprintf (file, "  In states (Id, a):\t");
    for (j = 0; j < mo->s[i].in_states; j++)
      fprintf (file, "FIXME: in_a is a matrix"/*(%d, %.3f) \t", mo->s[i].in_id[j], mo->s[i].in_a[j]*/);
    fprintf (file, "\n");
  }
}                               /* ghmm_dmodel_states_print */

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

void ghmm_dsmodel_Ak_print (FILE * file, ghmm_dsmodel * mo, int k, char *tab,
                       char *separator, char *ending)
{
  int i, j, out_state;
  for (i = 0; i < mo->N; i++) {
    out_state = 0;
    fprintf (file, "%s", tab);
    if (mo->s[i].out_states > 0 && mo->s[i].out_id[out_state] == 0) {
      fprintf (file, "%.2f", mo->s[i].out_a[k][out_state]);
      out_state++;
    }
    else
      fprintf (file, "0.00");
    for (j = 1; j < mo->N; j++) {
      if (mo->s[i].out_states > out_state && mo->s[i].out_id[out_state] == j) {
        fprintf (file, "%s %.2f", separator, mo->s[i].out_a[k][out_state]);
        out_state++;
      }
      else
        fprintf (file, "%s 0.00", separator);
    }
    fprintf (file, "%s\n", ending);
  }
}                               /* ghmm_dmodel_A_print */

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

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

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

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


void ghmm_dsmodel_from_dmodel (ghmm_dsmodel * smo, const ghmm_dmodel * mo, int klass)
{
#define CUR_PROC "ghmm_dsmodel_from_dmodel"
  int i, j, m, nachf, vorg;

  for (i = 0; i < mo->N; i++) {
    nachf = mo->s[i].out_states;
    vorg = mo->s[i].in_states;
    /* Copy the values */
    for (j = 0; j < nachf; j++) {
      smo->s[i].out_a[klass][j] = mo->s[i].out_a[j];
      smo->s[i].out_id[j] = mo->s[i].out_id[j];
    }
    for (j = 0; j < vorg; j++) {
      smo->s[i].in_a[klass][j] = mo->s[i].in_a[j];
      smo->s[i].in_id[j] = mo->s[i].in_id[j];
    }
    for (m = 0; m < mo->M; m++)
      smo->s[i].b[m] = mo->s[i].b[m];
    smo->s[i].pi = mo->s[i].pi;
    smo->s[i].out_states = nachf;
    smo->s[i].in_states = vorg;
  }
  smo->prior = mo->prior;
#undef CUR_PROC
}

/*===========================================================================*/
ghmm_dmodel * ghmm_dsmodel_to_dmodel(const ghmm_dsmodel * mo, int kclass)
{
#define CUR_PROC "ghmm_dsmodel_to_dmodel"
  /*
   * Set the pointer appropriately
   */
  int i, j, nachf, vorg, m;
  ghmm_dmodel *m2 = NULL;
  ARRAY_CALLOC (m2, 1);
  ARRAY_CALLOC (m2->s, mo->N);
  for (i = 0; i < mo->N; i++) {
    nachf = mo->s[i].out_states;
    vorg = mo->s[i].in_states;
    ARRAY_CALLOC (m2->s[i].out_id, nachf);
    ARRAY_CALLOC (m2->s[i].out_a, nachf);
    ARRAY_CALLOC (m2->s[i].in_id, vorg);
    ARRAY_CALLOC (m2->s[i].in_a, vorg);
    ARRAY_CALLOC (m2->s[i].b, mo->M);
    /* Copy the values */
    for (j = 0; j < nachf; j++) {
      m2->s[i].out_a[j] = mo->s[i].out_a[kclass][j];
      m2->s[i].out_id[j] = mo->s[i].out_id[j];
    }
    for (j = 0; j < vorg; j++) {
      m2->s[i].in_a[j] = mo->s[i].in_a[kclass][j];
      m2->s[i].in_id[j] = mo->s[i].in_id[j];
    }
    for (m = 0; m < mo->M; m++)
      m2->s[i].b[m] = mo->s[i].b[m];
    m2->s[i].pi = mo->s[i].pi;
    m2->s[i].out_states = nachf;
    m2->s[i].in_states = vorg;
  }
  m2->N = mo->N;
  m2->M = mo->M;
  m2->prior = mo->prior;

  m2->model_type = mo->model_type;
  if (mo->model_type & GHMM_kSilentStates) {
    assert (mo->silent != NULL);
    ARRAY_CALLOC (m2->silent, mo->N);
    for (i = 0; i < m2->N; i++) {
      m2->silent[i] = mo->silent[i];
    }
    m2->topo_order_length = mo->topo_order_length;
    ARRAY_CALLOC (m2->topo_order, mo->topo_order_length);
    for (i = 0; i < m2->topo_order_length; i++) {
      m2->topo_order[i] = mo->topo_order[i];
    }
  }
  return (m2);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  m_free (m2->silent);
  m_free (m2->topo_order);
  ghmm_dmodel_free (&m2);
  return (NULL);
#undef CUR_PROC
}

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

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

/*state* state_copy(ghmm_dstate *my_state) {
  ghmm_dstate* new_state = (ghmm_dstate*) malloc(sizeof(state));

  state_copy_to(my_state,new_state);

  return new_state;

  }*/ /* state_copy */

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

/*void state_copy_to(ghmm_dstate *source, ghmm_dstate* dest) {
  dest->pi         = source->pi;
  dest->out_states = source->out_states;
  dest->in_states  = source->in_states;
  dest->fix        = source->fix;

  dest->b          = malloc(xxx);
  memcpy(dest->b,source->b,xxx);

  dest->out_id     = malloc(xxx);
  memcpy(dest->out_id,source->out_id,xxx);

  dest->in_id      = malloc(xxx);
  memcpy(dest->in_id,source->in_id,xxx);

  dest->out_a      = malloc(xxx);
  memcpy(dest->out_a,source->out_a,xxx);

  dest->in_a       = malloc(xxx);
  memcpy(dest->in_a,source->in_a,xxx);
  }*/ /* state_copy_to */


/*===================== E n d   o f  f i l e  "model.c"       ===============*/