smodel.c

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

        goto STOP;
      }
      if (!n_read) {
        ighmm_scanner_error (s, "need N as a range for A");
        goto STOP;
      }
      if (a_read) {
        ighmm_scanner_error (s, "identifier A twice");
        goto STOP;
      }
      ARRAY_CALLOC (a_3dmatrix, smo->cos);
      for (osc = 0; osc < smo->cos; osc++) {
        ARRAY_CALLOC (a_3dmatrix[osc], smo->N);
        ighmm_scanner_get_name (s);
        if (!strcmp (s->id, "matrix")) {
          if (ighmm_cmatrix_read (s, a_3dmatrix[osc], smo->N, smo->N)) {
            ighmm_scanner_error (s, "unable to read matrix Ak");
            goto STOP;
          }
        }
        else {
          ighmm_scanner_error (s, "unknown identifier");
          goto STOP;
        }
        if (osc < smo->cos - 1) {
          ighmm_scanner_consume (s, ';');
          if (s->err)
            goto STOP;
          /* read next matrix */
          ighmm_scanner_get_name (s);
          if (strncmp (s->id, "Ak_", 3)) {
            ighmm_scanner_error (s, "not enough matrices Ak");
            goto STOP;
          }
          ighmm_scanner_consume (s, '=');
          if (s->err)
            goto STOP;
        }
      }
      a_read = 1;
    }
    else if (!strcmp (s->id, "C")) {    /* weight for output components */
      if ((!n_read) || (!m_read)) {
        ighmm_scanner_error (s, "need M and N as a range for C");
        goto STOP;
      }
      if (c_read) {
        ighmm_scanner_error (s, "identifier C twice");
        goto STOP;
      }
      ARRAY_CALLOC (c_matrix, smo->N);
      ighmm_scanner_get_name (s);
      if (!strcmp (s->id, "matrix")) {
        if (ighmm_cmatrix_read (s, c_matrix, smo->N, M)) {
          ighmm_scanner_error (s, "unable to read matrix C");
          goto STOP;
        }
        c_read = 1;
      }
      else {
        ighmm_scanner_error (s, "unknown identifier");
        goto STOP;
      }
    }
    else if (!strcmp (s->id, "Mue")) {  /* mean of normal density */
      if ((!n_read) || (!m_read)) {
        ighmm_scanner_error (s, "need M and N as a range for Mue");
        goto STOP;
      }
      if (mue_read) {
        ighmm_scanner_error (s, "identifier Mue twice");
        goto STOP;
      }
      ARRAY_CALLOC (mue_matrix, smo->N);
      ighmm_scanner_get_name (s);
      if (!strcmp (s->id, "matrix")) {
        if (ighmm_cmatrix_read (s, mue_matrix, smo->N, M)) {
          ighmm_scanner_error (s, "unable to read matrix Mue");
          goto STOP;
        }
        mue_read = 1;
      }
      else {
        ighmm_scanner_error (s, "unknown identifier");
        goto STOP;
      }
    }
    else if (!strcmp (s->id, "U")) {    /* variances of normal density */
      if ((!n_read) || (!m_read)) {
        ighmm_scanner_error (s, "need M and N as a range for U");
        goto STOP;
      }
      if (u_read) {
        ighmm_scanner_error (s, "identifier U twice");
        goto STOP;
      }
      ARRAY_CALLOC (u_matrix, smo->N);
      ighmm_scanner_get_name (s);
      if (!strcmp (s->id, "matrix")) {
        if (ighmm_cmatrix_read (s, u_matrix, smo->N, M)) {
          ighmm_scanner_error (s, "unable to read matrix U");
          goto STOP;
        }
        u_read = 1;
      }
      else {
        ighmm_scanner_error (s, "unknown identifier");
        goto STOP;
      }
    }
    else {
      ighmm_scanner_error (s, "unknown identifier");
      goto STOP;
    }
    ighmm_scanner_consume (s, ';');
    if (s->err)
      goto STOP;
  }                             /* while(..s->c-'}') */
  ighmm_scanner_consume (s, '}');
  if (s->err)
    goto STOP;

  if (m_read == 0 || n_read == 0 || pi_read == 0 || a_read == 0 ||
      c_read == 0 || mue_read == 0 || u_read == 0 || density_read == 0) {
    ighmm_scanner_error (s,
                   "some identifier not specified (N, M, Pi, A, C, Mue, U or density)");
    goto STOP;
  }
  /* set prior to -1 it none was specified */
  if (prior_read == 0)
    smo->prior = -1;
  /* default for fix is 0 */
  if (fix_read == 0) {
    ARRAY_CALLOC (fix_vektor, smo->N);
    for (i = 0; i < smo->N; i++)
      fix_vektor[i] = 0;
  }
  /* memory alloc for all transition matrices. If a transition is possible in one
     class --> alloc memory for all classes */
  for (i = 0; i < smo->N; i++) {
    /*smo->density[i] = density;*/
    for (j = 0; j < smo->N; j++) {
      out = in = 0;
      for (osc = 0; osc < smo->cos; osc++) {
        if (a_3dmatrix[osc][i][j] > 0.0)
          out = 1;
        if (a_3dmatrix[osc][j][i] > 0.0)
          in = 1;
      }
      smo->s[i].out_states += out;
      smo->s[i].in_states += in;
    }
    if (ghmm_cstate_alloc (smo->s + i, M, smo->s[i].in_states,
                            smo->s[i].out_states, smo->cos)) {
      GHMM_LOG_QUEUED(LCONVERTED);
      goto STOP;
    }
    /* given no support of the smo file, make all densities of same type and
       the cutoff value for the tail gaussian equal to the previously 
       used constant */
    smo->M = M;
    for (j=0; j < smo->M; j++){
      smo->s[i].density[j] = (ghmm_density_t)density;
      smo->s[i].a[j] = -GHMM_EPS_NDT;
      smo->s[i].M = M;
    }
    /* copy values read to smodel */
    if (smodel_copy_vectors
        (smo, i, pi_vektor, fix_vektor, a_3dmatrix, c_matrix, mue_matrix,
         u_matrix)) {
      GHMM_LOG_QUEUED(LCONVERTED);
      goto STOP;
    }

  }
  if (a_3dmatrix)
    for (i = 0; i < smo->cos; i++)
      ighmm_cmatrix_free (&(a_3dmatrix[i]), smo->N);
  m_free (a_3dmatrix);
  ighmm_cmatrix_free (&c_matrix, smo->N);
  ighmm_cmatrix_free (&mue_matrix, smo->N);
  ighmm_cmatrix_free (&u_matrix, smo->N);
  m_free (pi_vektor);
  m_free (fix_vektor);
  return (smo);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  if (a_3dmatrix)
    for (i = 0; i < smo->cos; i++)
      ighmm_cmatrix_free (&(a_3dmatrix[i]), smo->N);
  m_free (a_3dmatrix);
  ighmm_cmatrix_free (&c_matrix, smo->N);
  ighmm_cmatrix_free (&mue_matrix, smo->N);
  ighmm_cmatrix_free (&u_matrix, smo->N);
  m_free (pi_vektor);
  m_free (fix_vektor);
  ghmm_cmodel_free (&smo);
  return NULL;
#undef CUR_PROC
}                               /* ghmm_cmodel_read_block */

#endif /* GHMM_OBSOLETE */


/*============================================================================*/
ghmm_cmodel **ghmm_cmodel_xml_read(const char *filename, int *smo_number)
{
#define CUR_PROC "ghmm_cmodel_xml_read"
  ghmm_xmlfile* f;
  ghmm_cmodel** smo;

  f = ghmm_xmlfile_parse(filename);
  if (f) {
    assert(f->modelType &  GHMM_kContinuousHMM);
    *smo_number = (f->noModels);
    smo = f->model.c;
    
    free(f); /* XXX - by now, we free f */  
    return smo;
  } else
    return NULL;
#undef CUR_PROC
}


/*============================================================================*/
int ghmm_cmodel_xml_write(ghmm_cmodel** smo, const char* file, int smo_number)
{
#define CUR_PROC "ghmm_cmodel_xml_write"
  ghmm_xmlfile* f;

  ARRAY_MALLOC(f,1);
  f->noModels = smo_number;
  f->modelType = GHMM_kContinuousHMM;
  f->model.c = smo;
  ghmm_xmlfile_write(f, file);
  free(f);  
  return 0;
STOP:  
  return -1;
#undef CUR_PROC
}

/*============================================================================*/
int ghmm_cmodel_free (ghmm_cmodel ** smo)
{
#define CUR_PROC "ghmm_cmodel_free"
  int i;
  mes_check_ptr (smo, return (-1));

  for (i = 0; i < (*smo)->N && (*smo)->s; i++) {
    
    /* if there are no out_states field was never allocated */ 
    if ((*smo)->s[i].out_states > 0){
      m_free ((*smo)->s[i].out_id);
    }  
    /* if there are no in_states field was never allocated */ 
    if ((*smo)->s[i].in_states > 0){
      m_free ((*smo)->s[i].in_id);
    }  
    ighmm_cmatrix_free (&((*smo)->s[i].out_a), (*smo)->cos);
    ighmm_cmatrix_free (&((*smo)->s[i].in_a), (*smo)->cos);
    m_free ((*smo)->s[i].c);
    m_free ((*smo)->s[i].mue);
    m_free ((*smo)->s[i].u);
    m_free ((*smo)->s[i].a);
    m_free ((*smo)->s[i].mixture_fix);
    m_free ((*smo)->s[i].density);
    
  }
  if ((*smo)->s) m_free ((*smo)->s);

  if ((*smo)->class_change) {
    if ((*smo)->class_change->user_data) {
      m_free ((*smo)->class_change->user_data);
    }
    m_free ((*smo)->class_change);
  }
  /*m_free ((*smo)->density);*/
  m_free (*smo);
  return (0);
#undef CUR_PROC
}                               /* ghmm_cmodel_free */


/*============================================================================*/
ghmm_cmodel *ghmm_cmodel_copy (const ghmm_cmodel * smo)
{
# define CUR_PROC "ghmm_cmodel_copy"
  int i, k, j, nachf, vorg, m;
  ghmm_cmodel *sm2 = NULL;
  ARRAY_CALLOC (sm2, 1);
  ARRAY_CALLOC (sm2->s, smo->N);
  for (i = 0; i < smo->N; i++) {
    nachf = smo->s[i].out_states;
    vorg = smo->s[i].in_states;
    ARRAY_CALLOC (sm2->s[i].out_id, nachf);
    sm2->s[i].out_a = ighmm_cmatrix_alloc (smo->cos, nachf);
    if (!sm2->s[i].out_a) {
      GHMM_LOG_QUEUED(LCONVERTED);
      goto STOP;
    }
    ARRAY_CALLOC (sm2->s[i].in_id, vorg);
    sm2->s[i].in_a = ighmm_cmatrix_alloc (smo->cos, vorg);
    if (!sm2->s[i].in_a) {
      GHMM_LOG_QUEUED(LCONVERTED);
      goto STOP;
    }
    ARRAY_CALLOC (sm2->s[i].c, smo->s[i].M);
    ARRAY_CALLOC (sm2->s[i].mue, smo->s[i].M);
    ARRAY_CALLOC (sm2->s[i].u, smo->s[i].M);
    ARRAY_CALLOC (sm2->s[i].mixture_fix, smo->s[i].M);
    ARRAY_CALLOC (sm2->s[i].a, smo->s[i].M);
    ARRAY_CALLOC (sm2->s[i].density, smo->s[i].M);
    
    /* copy values */
    for (j = 0; j < nachf; j++) {
      for (k = 0; k < smo->cos; k++)
        sm2->s[i].out_a[k][j] = smo->s[i].out_a[k][j];
      sm2->s[i].out_id[j] = smo->s[i].out_id[j];
    }
    for (j = 0; j < vorg; j++) {
      for (k = 0; k < smo->cos; k++)
        sm2->s[i].in_a[k][j] = smo->s[i].in_a[k][j];
      sm2->s[i].in_id[j] = smo->s[i].in_id[j];
    }
    for (m = 0; m < smo->s[i].M; m++) {
      sm2->s[i].c[m] = smo->s[i].c[m];
      sm2->s[i].mue[m] = smo->s[i].mue[m];
      sm2->s[i].u[m] = smo->s[i].u[m];
      sm2->s[i].mixture_fix[m] = smo->s[i].mixture_fix[m];
      sm2->s[i].a[m] = smo->s[i].a[m];
      sm2->s[i].density[m] = smo->s[i].density[m];
    }
    sm2->s[i].M = smo->s[i].M;
    sm2->s[i].pi = smo->s[i].pi;
    sm2->s[i].fix = smo->s[i].fix;
    sm2->s[i].out_states = nachf;
    sm2->s[i].in_states = vorg;
  }
  sm2->cos = smo->cos;
  sm2->N = smo->N;
  sm2->M = smo->M;
  sm2->prior = smo->prior;
  return (sm2);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ghmm_cmodel_free (&sm2);
  return (NULL);
# undef CUR_PROC
}                               /* ghmm_cmodel_copy */


/*============================================================================*/
int ghmm_cmodel_check (const ghmm_cmodel * smo)
{
# define CUR_PROC "ghmm_cmodel_check"
  int valid = 0;
  double sum;
  int i, k, j;
  char * str;
  /* sum  Pi[i] == 1 ? */
  sum = 0.0;
  
  for (i = 0; i < smo->N; i++) {
    sum += smo->s[i].pi;
  }
  if (fabs (sum - 1.0) >= GHMM_EPS_PREC) {
    GHMM_LOG(LCONVERTED, "sum Pi[i] != 1.0\n");
    valid = -1;
    /*goto STOP; */
  }
  /* only 0/1 in fix? */
  for (i = 0; i < smo->N; i++) {
    if (smo->s[i].fix != 0 && smo->s[i].fix != 1) {
      GHMM_LOG(LCONVERTED, "in vector fix_state only 0/1 values!\n");
      valid = -1;
      /*goto STOP;*/
    }
  }
  for (i = 0; i < smo->N; i++) {
    if (smo->s[i].out_states == 0) {
      str = ighmm_mprintf (NULL, 0, "out_states = 0 (state %d -> final state!)\n", i);
      GHMM_LOG(LCONVERTED, str);
      m_free (str);
    }
    /* sum  a[i][k][j] */
    for (k = 0; k < smo->cos; k++) {
      sum = 0.0;
      for (j = 0; j < smo->s[i].out_states; j++) {
        sum += smo->s[i].out_a[k][j];
      }
      if (fabs (sum - 1.0) >= GHMM_EPS_PREC) {
        str = ighmm_mprintf (NULL, 0, "sum out_a[j] = %.4f != 1.0 (state %d, class %d)\n", sum, i,k);
        GHMM_LOG(LCONVERTED, str);
        m_free (str);
        valid = -1;
        /*goto STOP; */
      }
    }
    /* sum c[j] */
    sum = 0.0;
    for (j = 0; j < smo->s[i].M; j++)
      sum += smo->s[i].c[j];
    if (fabs (sum - 1.0) >= GHMM_EPS_PREC) {
      str = ighmm_mprintf (NULL, 0, "sum c[j] = %.2f != 1.0 (state %d)\n", sum, i);
      GHMM_LOG(LCONVERTED, str);
      m_free (str);
      valid = -1;            
      /* goto STOP; */
    }
    /* check mue, u ? */
  }

/* for lazy-evaluation-like model checking
   uncomment 'goto STOP' statements and 'STOP:' line  */  
/* STOP: */
  return (valid);
# undef CUR_PROC
}                               /* ghmm_cmodel_check */


/*============================================================================*/
int ghmm_cmodel_check_compatibility (ghmm_cmodel ** smo, int smodel_number)
{
#define CUR_PROC "ghmm_cmodel_check_compatibility"
/* XXX - old function not used any more !!! */
  int i, j;
  char * str;

  for (i = 0; i < smodel_number; i++)
    for (j = i + 1; j < smodel_number; j++) {
      if (smo[i]->N != smo[j]->N) {
        str = ighmm_mprintf (NULL, 0,
                "ERROR: different number of states in smodel %d (%d) and smodel %d (%d)",
                i, smo[i]->N, j, smo[j]->N);
        GHMM_LOG(LCONVERTED, str);
        m_free (str);
        return (-1);
      }
      if (smo[i]->s[0].M != smo[j]->s[0].M) {
        str = ighmm_mprintf (NULL, 0,
                "ERROR: different number of possible outputs in smodel  %d (%d) and smodel %d (%d)",
                i, smo[i]->s[0].M, j, smo[j]->s[0].M);
        GHMM_LOG(LCONVERTED, str);
        m_free (str);
        return (-1);
      }
    }
  return 0;
#undef CUR_PROC
}                               /* ghmm_cmodel_check_compatibility */


/*============================================================================*/
double ghmm_cmodel_get_random_var (ghmm_cmodel * smo, int state, int m)
{
# define CUR_PROC "ghmm_cmodel_get_random_var"
  switch (smo->s[state].density[m]) {
  case normal_approx:
  case normal:
    return (ighmm_rand_normal (smo->s[state].mue[m], smo->s[state].u[m], 0));
  case normal_right:
    return (ighmm_rand_normal_right (smo->s[state].a[m],smo->s[state].mue[m], smo->s[state].u[m], 0));
  case normal_left:
    return -(ighmm_rand_normal_right (-smo->s[state].a[m],-smo->s[state].mue[m], smo->s[state].u[m], 0));
  case uniform:
	return (ighmm_rand_uniform_cont (0, smo->s[state].mue[m] ,smo->s[state].u[m]));
  default:
    ighmm_mes (MES_WIN, "Warning: density function not specified!\n");
    return (-1);
  }
# undef CUR_PROC
}                               /* ghmm_cmodel_get_random_var */