model.c

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

  sq = ghmm_dseq_calloc(seq_number);
  if (!sq) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }

  /* allocating additional fields for the state sequence in the ghmm_dseq struct */
  ARRAY_CALLOC(sq->states, seq_number);
  ARRAY_CALLOC(sq->states_len, seq_number);

  /* allocating additional fields for the labels in the ghmm_dseq struct */
  ARRAY_CALLOC(sq->state_labels, seq_number);
  ARRAY_CALLOC(sq->state_labels_len, seq_number);

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

  if (seed > 0) {
    GHMM_RNG_SET(RNG, seed);
  }

  /* initialize the emission history */
  mo->emission_history = 0;

  while (n < seq_number) {
    ARRAY_CALLOC(sq->seq[n], len);
    
    /* for silent models we have to allocate for the maximal possible number
       of lables and states */
    if (mo->model_type & GHMM_kSilentStates) {
      ARRAY_CALLOC(sq->states[n], len * mo->N);
      ARRAY_CALLOC(sq->state_labels[n], len * mo->N);
    }
     else {
      ARRAY_CALLOC(sq->states[n], len);
      ARRAY_CALLOC(sq->state_labels[n], len);
    }

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

    while (pos < len) {
      /* save the state path and label */
      sq->states[n][label_pos] = state;
      sq->state_labels[n][label_pos] = mo->label[state];
      label_pos++;

      /* Get a random output m if the state is not a silent state */
      if (!(mo->model_type & GHMM_kSilentStates) || !(mo->silent[state])) {
        m = get_random_output(mo, state, pos);
        update_emission_history(mo, m);
        sq->seq[n][pos] = m;
        pos++;
      }

      /* get next state */
      p = GHMM_RNG_UNIFORM(RNG);
      if (pos < mo->maxorder) {
        max_sum = 0.0;
        for (j = 0; j < mo->s[state].out_states; j++) {
          j_id = mo->s[state].out_id[j];
          if (!(mo->model_type & GHMM_kHigherOrderEmissions) || mo->order[j_id] < pos)
            max_sum += mo->s[state].out_a[j];
        }
        if (j && fabs(max_sum) < GHMM_EPS_PREC) {
          GHMM_LOG_PRINTF(LERROR, LOC, "No possible transition from state %d "
                          "since all successor states require more history "
                          "than seen up to position: %d.",
                          state, pos);
          break;
        }
        if (j)
          p *= max_sum;
      }

      sum = 0.0;
      for (j = 0; j < mo->s[state].out_states; j++) {
        j_id = mo->s[state].out_id[j];
        if (!(mo->model_type & GHMM_kHigherOrderEmissions) || mo->order[j_id] < pos) {
          sum += mo->s[state].out_a[j];
          if (sum >= p)
            break;
        }
      }

      if (sum == 0.0) {
        GHMM_LOG_PRINTF(LINFO, LOC, "final state (%d) reached at position %d "
                        "of sequence %d", state, pos, n);
	break;
      }

      state = j_id;
    }                           /* while (pos < len) */

    /* realocate state path and label sequence to actual size */ 
    if (mo->model_type & GHMM_kSilentStates) {
      ARRAY_REALLOC(sq->states[n], label_pos);
      ARRAY_REALLOC(sq->state_labels[n], label_pos);
    }

    sq->seq_len[n] = pos;
    sq->states_len[n] = label_pos;
    sq->state_labels_len[n] = label_pos;
    n++;
  }                             /* while( n < seq_number ) */

  return (sq);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ghmm_dseq_free(&sq);
  return NULL;
#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
/* Scales the output and transitions probs of all states in a given model */
int ghmm_dmodel_normalize (ghmm_dmodel * mo)
{
#define CUR_PROC "ghmm_dmodel_normalize"
  double pi_sum=0.0;
  int i, j, m, j_id, i_id=0, res=0;
  int size = 1;
  char *str;

  for (i = 0; i < mo->N; i++) {
    if (mo->s[i].pi >= 0.0)
      pi_sum += mo->s[i].pi;
    else
      mo->s[i].pi = 0.0;

    /* check model_type before using state order */
    if (mo->model_type & GHMM_kHigherOrderEmissions)
      size = ghmm_ipow (mo, mo->M, mo->order[i]);

    /* normalize transition probabilities */
    if (ighmm_cvector_normalize (mo->s[i].out_a, mo->s[i].out_states) == -1) {
      res = -1;
    }
    /* for every outgoing probability update the corrosponding incoming probability */
    for (j = 0; j < mo->s[i].out_states; j++) {
      j_id = mo->s[i].out_id[j];
      for (m = 0; m < mo->s[j_id].in_states; m++) {
        if (i == mo->s[j_id].in_id[m]) {
          i_id = m;
          break;
        }
      }
      if (i_id == mo->s[j_id].in_states) {
        str = ighmm_mprintf (NULL, 0, "Outgoing transition from state %d to \
           state %d has no corresponding incoming transition.\n", i, j_id);
        GHMM_LOG(LCONVERTED, str);
        return -1;
      }
      mo->s[j_id].in_a[i_id] = mo->s[i].out_a[j];
    }
    /* normalize emission probabilities, but not for silent states */
    if (!((mo->model_type & GHMM_kSilentStates) && mo->silent[i])) {
      for (m = 0; m < size; m++) {
        if (ighmm_cvector_normalize (&(mo->s[i].b[m * mo->M]), mo->M) == -1)
          res = -1;
      }
    }
  }
  for (i = 0; i < mo->N; i++)
    mo->s[i].pi /= pi_sum;

  return res;
#undef CUR_PROC
}                               /* ghmm_dmodel_normalize */


/*----------------------------------------------------------------------------*/
int ghmm_dmodel_add_noise (ghmm_dmodel * mo, double level, int seed)
{
#define CUR_PROC "model_add_noise_A"

  int h, i, j, hist;
  int size = 1;

  if (level > 1.0)
    level = 1.0;

  for (i = 0; i < mo->N; i++) {
    for (j = 0; j < mo->s[i].out_states; j++)
      /* add noise only to out_a, in_a is updated on normalisation */
      mo->s[i].out_a[j] *= (1 - level) + (GHMM_RNG_UNIFORM (RNG) * 2 * level);

    if (mo->model_type & GHMM_kHigherOrderEmissions)
      size = ghmm_ipow (mo, mo->M, mo->order[i]);
    for (hist = 0; hist < size; hist++)
      for (h = hist * mo->M; h < hist * mo->M + mo->M; h++)
        mo->s[i].b[h] *= (1 - level) + (GHMM_RNG_UNIFORM (RNG) * 2 * level);
  }

  return ghmm_dmodel_normalize (mo);

#undef CUR_PROC
}

/*----------------------------------------------------------------------------*/
static int ghmm_dstate_transition_add(ghmm_dstate * s, int start, int dest, double prob)
{
#define CUR_PROC "ghmm_dmodel_transition_add"

  int i;

  /* resize the arrays */
  ARRAY_REALLOC (s[dest].in_id, s[dest].in_states + 1);
  ARRAY_REALLOC (s[dest].in_a, s[dest].in_states + 1);
  ARRAY_REALLOC (s[start].out_id, s[start].out_states + 1);
  ARRAY_REALLOC (s[start].out_a, s[start].out_states + 1);

  s[dest].in_states += 1;
  s[start].out_states += 1;

  /* search the right place to insert while moving greater entrys one field back */
  for (i = s[start].out_states - 1; i >= 0; i--) {
    if (i == 0 || dest > s[start].out_id[i - 1]) {
      s[start].out_id[i] = dest;
      s[start].out_a[i] = prob;
      break;
    }
    else {
      s[start].out_id[i] = s[start].out_id[i - 1];
      s[start].out_a[i] = s[start].out_a[i - 1];
    }
  }

  /* search the right place to insert while moving greater entrys one field back */
  for (i = s[dest].in_states - 1; i >= 0; i--)
    if (i == 0 || start > s[dest].in_id[i - 1]) {
      s[dest].in_id[i] = start;
      s[dest].in_a[i] = prob;
      break;
    }
    else {
      s[dest].in_id[i] = s[dest].in_id[i - 1];
      s[dest].in_a[i] = s[dest].in_a[i - 1];
    }

  return 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return -1;
#undef CUR_PROC
}

/*----------------------------------------------------------------------------*/
static int ghmm_dstate_transition_del(ghmm_dstate * s, int start, int dest)
{
#define CUR_PROC "ghmm_dmodel_transition_del"

  int i, j;

  /* search ... */
  for (j = 0; dest != s[start].out_id[j]; j++)
    if (j == s[start].out_states) {
      GHMM_LOG(LCONVERTED, "No such transition");
      return -1;
    }
  /* ... and replace outgoing */
  for (i = j + 1; i < s[start].out_states; i++) {
    s[start].out_id[i - 1] = s[start].out_id[i];
    s[start].out_a[i - 1] = s[start].out_a[i];
  }

  /* search ... */
  for (j = 0; start != s[dest].in_id[j]; j++)
    if (j == s[dest].in_states) {
      GHMM_LOG(LCONVERTED, "No such transition");
      return -1;
    }
  /* ... and replace incoming */
  for (i = j + 1; i < s[dest].in_states; i++) {
    s[dest].in_id[i - 1] = s[dest].in_id[i];
    s[dest].in_a[i - 1] = s[dest].in_a[i];
  }

  /* reset number */
  s[dest].in_states -= 1;
  s[start].out_states -= 1;

  /* free memory */
  ARRAY_REALLOC (s[dest].in_id, s[dest].in_states);
  ARRAY_REALLOC (s[dest].in_a, s[dest].in_states);
  ARRAY_REALLOC (s[start].out_id, s[start].out_states);
  ARRAY_REALLOC (s[start].out_a, s[start].out_states);

  return 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return -1;
#undef CUR_PROC
}

/*----------------------------------------------------------------------------*/
/** 
   Allocates a new ghmm_dbackground struct and assigs the arguments to
   the respective fields. Note: The arguments need allocation outside of this
   function.
   
   @return     :               0 on success, -1 on error
   @param mo   :               one model
   @param cur  :               a id of a state
   @param times:               number of times the state cur is at least evaluated
*/
int ghmm_dmodel_duration_apply (ghmm_dmodel * mo, int cur, int times)
{
#define CUR_PROC "ghmm_dmodel_duration_apply"

  int i, j, last, size, failed=0;

  if (mo->model_type & GHMM_kSilentStates) {
    GHMM_LOG(LCONVERTED, "Sorry, apply_duration doesn't support silent states yet\n");
    return -1;
  }

  last = mo->N;
  mo->N += times - 1;

  ARRAY_REALLOC (mo->s, mo->N);
  if (mo->model_type & GHMM_kSilentStates) {
    ARRAY_REALLOC (mo->silent, mo->N);
    ARRAY_REALLOC (mo->topo_order, mo->N);
  }
  if (mo->model_type & GHMM_kTiedEmissions)
    ARRAY_REALLOC (mo->tied_to, mo->N);
  if (mo->model_type & GHMM_kBackgroundDistributions)
    ARRAY_REALLOC (mo->background_id, mo->N);

  size = ghmm_ipow (mo, mo->M, mo->order[cur] + 1);
  for (i = last; i < mo->N; i++) {
    /* set the new state */
    mo->s[i].pi = 0.0;
    mo->order[i] = mo->order[cur];
    mo->s[i].fix = mo->s[cur].fix;
    mo->label[i] = mo->label[cur];
    mo->s[i].in_a = NULL;
    mo->s[i].in_id = NULL;
    mo->s[i].in_states = 0;
    mo->s[i].out_a = NULL;
    mo->s[i].out_id = NULL;
    mo->s[i].out_states = 0;

    ARRAY_MALLOC (mo->s[i].b, size);
    for (j = 0; j < size; j++)
      mo->s[i].b[j] = mo->s[cur].b[j];

    if (mo->model_type & GHMM_kSilentStates) {
      mo->silent[i] = mo->silent[cur];
      /* XXX what to do with topo_order
         mo->topo_order[i] = ????????????; */
    }
    if (mo->model_type & GHMM_kTiedEmissions)
      /* XXX is there a clean solution for tied states?
         what if the current state is a tie group leader?
         the last added state should probably become
         the new tie group leader */
      mo->tied_to[i] = GHMM_kUntied;
    if (mo->model_type & GHMM_kBackgroundDistributions)
      mo->background_id[i] = mo->background_id[cur];
  }

  /* move the outgoing transitions to the last state */
  while (mo->s[cur].out_states > 0) {
    if (mo->s[cur].out_id[0] == cur) {
      ghmm_dstate_transition_add (mo->s, mo->N - 1, mo->N - 1, mo->s[cur].out_a[0]);
      ghmm_dstate_transition_del (mo->s, cur, mo->s[cur].out_id[0]);
    }
    else {
      ghmm_dstate_transition_add (mo->s, mo->N - 1, mo->s[cur].out_id[0],
                            mo->s[cur].out_a[0]);
      ghmm_dstate_transition_del (mo->s, cur, mo->s[cur].out_id[0]);
    }
  }

  /* set the linear transitions through all added states */
  ghmm_dstate_transition_add (mo->s, cur, last, 1.0);
  for (i = last + 1; i < mo->N; i++) {
    ghmm_dstate_transition_add (mo->s, i - 1, i, 1.0);
  }

  if (ghmm_dmodel_normalize (mo))
    goto STOP;

  return 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  /* Fail hard if these realloc fail. They shouldn't because we have the memory
     and try only to clean up! */
  if (failed++)
    exit (1);
  
  ARRAY_REALLOC (mo->s, last);
  ARRAY_REALLOC (mo->tied_to, last);
  ARRAY_REALLOC (mo->background_id, last);

  mo->N = last;
  return -1;
#undef CUR_PROC
}

/*----------------------------------------------------------------------------*/
/** 
   Allocates a new ghmm_dbackground struct and assigns
   the arguments to the respective fields.
   Note: The arguments need allocation outside of this function.
   
   @return:      new pointer to a ghmm_dbackground struct or NULL
   @param n:     number of distributions
   @param order: orders of the distribtions (optional)
   @param B:     matrix of distribution parameters (optional)
*/
ghmm_dbackground *ghmm_dbackground_alloc (int n, int m, int *orders, double **B)
{
#define CUR_PROC "ghmm_dbackground_alloc"
  ghmm_dbackground *ptbackground;

  ARRAY_CALLOC(ptbackground, 1);

  ptbackground->n = n;
  ptbackground->m = m;
  if (orders)
    ptbackground->order = orders;
  if (B)
    ptbackground->b = B;

  ARRAY_CALLOC(ptbackground->name, n);

  return ptbackground;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return NULL;
#undef CUR_PROC
}

/*----------