smodel.c

一个通用的隐性马尔可夫C代码库 开发环境:C语言 简要说明:这是一个通用的隐性马尔可夫C代码库

/*******************************************************************************
  author       : Bernhard Knab
  filename     : /zpr/bspk/src/hmm/ghmm/ghmm/smodel.c
  created      : TIME: 21:54:32     DATE: Sun 14. November 1999
  $Id: smodel.c,v 1.21 2004/01/05 19:23:12 schliep Exp $

Copyright (C) 1998-2001, ZAIK/ZPR, Universit鋞 zu K鰈n

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


*******************************************************************************/

#include <float.h>
#include <math.h>
#include "mprintf.h"
#include "mes.h"
#include "smodel.h"
#include "sfoba.h"
#include "sreestimate.h"
#include "matrix.h"
#include "vector.h"
#include "sequence.h"
#include "const.h"
#include "rng.h"
#include "randvar.h"
#include "string.h"

/*----------------------------------------------------------------------------*/
static int smodel_state_alloc(sstate *state,
			      int M,
			      int in_states,
			      int out_states,
			      int cos) {
# define CUR_PROC "smodel_state_alloc"
  int res = -1;
  if(!m_calloc(state->c, M)) {mes_proc(); goto STOP;}
  if(!m_calloc(state->mue, M)) {mes_proc(); goto STOP;}
  if(!m_calloc(state->u, M)) {mes_proc(); goto STOP;}
  if (out_states > 0) {
    if(!m_calloc(state->out_id, out_states)) {mes_proc(); goto STOP;}
    state->out_a = matrix_d_alloc(cos, out_states);
    if(!state->out_a) {mes_proc(); goto STOP;}
  }
  if (in_states > 0) {
    if(!m_calloc(state->in_id, in_states)) {mes_proc(); goto STOP;}
    state->in_a = matrix_d_alloc(cos, in_states);
    if(!state->in_a) {mes_proc(); goto STOP;}
  }
  res = 0;
STOP:
  return(res);
# undef CUR_PROC
} /* smodel_state_alloc */


/*----------------------------------------------------------------------------*/
static int smodel_copy_vectors(smodel *smo, int index, double *pi, int *fix,
			       double ***a_matrix, double **c_matrix, 
			       double **mue_matrix, double **u_matrix) {
#define CUR_PROC "smodel_alloc_vectors"
  int i, c, exist, count_out = 0, count_in = 0;
  smo->s[index].pi = pi[index];
  smo->s[index].fix = fix[index];
  for (i = 0; i < smo->M; i++){
    smo->s[index].c[i] = c_matrix[index][i];
    smo->s[index].mue[i] = mue_matrix[index][i];
    smo->s[index].u[i] = u_matrix[index][i];
  }

  for (i = 0; i < smo->N; i++) {
    exist = 0;
    for (c = 0; c < smo->cos; c++) {
      if (a_matrix[c][index][i]) { 
	exist = 1;
	break;
      }
    }
    /* transition to successor possible at least in one transition class */
    if (exist) { 
      if (count_out >= smo->s[index].out_states) {mes_proc(); return(-1);}
      smo->s[index].out_id[count_out] = i;
      for (c = 0; c < smo->cos; c++) 
	smo->s[index].out_a[c][count_out] = a_matrix[c][index][i];
      count_out++;
    }
    exist = 0;
    for (c = 0; c < smo->cos; c++) {
      if (a_matrix[c][i][index]) { 
	exist = 1;
	break;
      }
    }
    /* transition to predecessor possible at least in one transition class */
    if (exist) {
      if (count_in >= smo->s[index].in_states) {mes_proc(); return(-1);}
      smo->s[index].in_id[count_in] = i;
      for (c = 0; c < smo->cos; c++)
	smo->s[index].in_a[c][count_in] = a_matrix[c][i][index];
      count_in++;
    }
  }
  return(0);
#undef CUR_PROC
} /* smodel_alloc_vectors */


/*============================================================================*/
smodel **smodel_read(const char *filename, int *smo_number) {
#define CUR_PROC "smodel_read"
  int j;
  long new_models = 0;
  scanner_t *s = NULL;
  smodel **smo = NULL;
  *smo_number = 0;
  s = scanner_alloc(filename);  if(!s) {mes_proc(); goto STOP;}
  while(!s->err && !s->eof) {
    scanner_get_name(s);
    scanner_consume(s, '='); if(s->err) goto STOP;
    if (!strcmp(s->id, "SHMM") || !strcmp(s->id, "shmm")) {
      (*smo_number)++;
      /* more mem */
      if (m_realloc(smo, *smo_number)) { mes_proc(); goto STOP; }
      smo[*smo_number - 1] = smodel_read_block(s, (int *) &new_models); 
      if (!smo[*smo_number - 1]) { mes_proc(); goto STOP; }
      /* copy smodel */
      if (new_models > 1) { 
	/* "-1" due to  (*smo_number)++ from above  */
	if (m_realloc(smo, *smo_number - 1 + new_models)) { 
	  mes_proc(); goto STOP; 
	}
	for (j = 1; j < new_models; j++) {
	  smo[*smo_number] = smodel_copy(smo[*smo_number - 1]);
	  if (!smo[*smo_number]) { mes_proc(); goto STOP; }
	  (*smo_number)++;
	}
      }
    }  
    else {
      scanner_error(s, "unknown identifier");
      goto STOP;
    }
    scanner_consume(s, ';'); if(s->err) goto STOP; 
  } /* while(!s->err && !s->eof) */
  /*if (smodel_check_compatibility(smo, *smo_number) == -1) {
    mes_proc(); goto STOP;
  }*/
  return smo;
STOP:
  return NULL;
#undef CUR_PROC
} /* smodel_read */


/*============================================================================*/
smodel *smodel_read_block(scanner_t *s, int *multip){
#define CUR_PROC "smodel_read_block"
  int i,j,osc, m_read, n_read, pi_read, a_read, c_read, mue_read, cos_read, 
    u_read, len, density_read, out, in, prior_read, fix_read;
  smodel *smo      = NULL;
  double *pi_vektor = NULL, **a_matrix = NULL, ***a_3dmatrix = NULL;
  double **c_matrix = NULL, **mue_matrix = NULL, **u_matrix = NULL;
  int *fix_vektor = NULL;
  
  m_read = n_read = pi_read = a_read = c_read = mue_read = u_read 
    = density_read = prior_read = fix_read = cos_read = 0;
  *multip = 1; /* default */
  if (!(m_calloc(smo, 1))) { mes_proc(); goto STOP; }
  scanner_consume( s, '{' ); if(s->err) goto STOP; 
  while(!s->err && !s->eof && s->c - '}') {    
    scanner_get_name(s);
    if (strcmp(s->id, "M") && strcmp(s->id, "N") && strcmp(s->id, "Pi") &&
	strcmp(s->id, "A") && strcmp(s->id, "C") && strcmp(s->id, "Mue") &&
	strcmp(s->id, "U") && strcmp(s->id, "multip") && strcmp(s->id, "cos") &&
	strcmp(s->id, "prior") && strcmp(s->id, "fix_state") &&
	strcmp(s->id, "density") && strncmp(s->id, "Ak_", 3) ) {
      scanner_error(s, "unknown identifier"); goto STOP;
    }    
    scanner_consume(s, '='); if(s->err) goto STOP;
    if (!strcmp(s->id, "multip")) {
      *multip = scanner_get_int(s);
      if (*multip < 1) {/* ignore: makes no sense */
	*multip = 1;
	mes_prot("Multip < 1 ignored\n");
      }
    }
    else if (!strcmp(s->id, "M")) {/* number of output components */
      if (m_read) {scanner_error(s, "identifier M twice"); goto STOP;}
      smo->M = scanner_get_int(s);
      m_read = 1;
    }
    else if (!strcmp(s->id, "N")) {/* number of states */
      if (n_read) {scanner_error(s, "identifier N twice"); goto STOP;}
      smo->N = scanner_get_int(s);
      if (!m_calloc(smo->s, smo->N)) {mes_proc(); goto STOP;}
      n_read = 1;
    }   
    else if (!strcmp(s->id, "density")) {/* which density function? */
      if (density_read) {scanner_error(s,"identifier density twice");goto STOP;}
      smo->density = (density_t)scanner_get_int(s);
      if ((int)smo->density < 0 || smo->density >= density_number)
	{ scanner_error(s, "unknown typ of density function"); goto STOP; }
      density_read = 1;
    }   
    else if (!strcmp(s->id, "prior")) {/* modelprior */
      if (prior_read) {scanner_error(s,"identifier prior twice");goto STOP;}
      smo->prior = scanner_get_edouble(s);
      if ((smo->prior < 0 || smo->prior > 1) && smo->prior != -1)
	{ scanner_error(s, "invalid model prior"); goto STOP; }
      prior_read = 1;
    }
    else if (!strcmp(s->id, "cos")) {/* number of transition classes */
      if (cos_read) {scanner_error(s,"identifier cos twice");goto STOP;}
      smo->cos = scanner_get_int(s);
      if (smo->cos <= 0)
	{ scanner_error(s, "invalid model cos"); goto STOP; }
      cos_read = 1;
    }   
    else if (!strcmp(s->id, "Pi")) {/* Initial State Prob. */
      if (!n_read) {scanner_error(s, "need N as a range for Pi"); goto STOP;}
      if (pi_read) {scanner_error(s, "identifier Pi twice"); goto STOP;}
      scanner_get_name(s);
      if (!strcmp(s->id, "vector")) {
	scanner_consume( s, '{' ); if(s->err) goto STOP; 
	pi_vektor = scanner_get_double_earray(s, &len);
	if (len != smo->N) 
	  {scanner_error(s, "wrong number of elements in PI"); goto STOP;}
	scanner_consume( s, ';' ); if(s->err) goto STOP;
	scanner_consume( s, '}' ); if(s->err) goto STOP; 
	pi_read = 1;
      }
      else {
	scanner_error(s, "unknown identifier"); goto STOP;
      }
    }
    else if (!strcmp(s->id, "fix_state")) {
      if (!n_read) {scanner_error(s, "need N as a range for fix_state"); goto STOP;}
      if (fix_read) {scanner_error(s, "identifier fix_state twice"); goto STOP;}
      scanner_get_name(s);
      if (!strcmp(s->id, "vector")) {
	scanner_consume( s, '{' ); if(s->err) goto STOP; 
	fix_vektor = scanner_get_int_array(s, &len);
	if (len != smo->N) 
	  {scanner_error(s, "wrong number of elements in fix_state"); goto STOP;}
	scanner_consume( s, ';' ); if(s->err) goto STOP;
	scanner_consume( s, '}' ); if(s->err) goto STOP; 
	fix_read = 1;
      }
      else {
	scanner_error(s, "unknown identifier"); goto STOP;
      }
    }
    /* 1. Possibility: one matrix for all transition classes */
    else if (!strcmp(s->id, "A")) {
      if (!cos_read) {scanner_error(s, "cos unknown (needed for dim(A))"); goto STOP;}
      if (!n_read) {scanner_error(s, "need N as a range for A"); goto STOP;}
      if (a_read) {scanner_error(s, "Identifier A twice"); goto STOP;}
      if (!m_calloc(a_matrix, smo->N)) {mes_proc(); goto STOP;}
      scanner_get_name(s);
      if (!strcmp(s->id, "matrix")) {
	if (matrix_d_read(s, a_matrix, smo->N, smo->N)){
	  scanner_error(s, "unable to read matrix A"); goto STOP;
	}
	a_read = 1;
      }
      else {
	scanner_error(s, "unknown identifier"); goto STOP;
      }
      /* copy transition matrix to all transition classes */
      if (!m_calloc(a_3dmatrix, smo->cos)) {mes_proc(); goto STOP;}
      a_3dmatrix[0] = a_matrix;
      for (i = 1; i < smo->cos; i++) {
	a_3dmatrix[i] = matrix_d_alloc_copy(smo->N, smo->N, a_matrix);
	if (!a_3dmatrix[i]) {mes_proc(); goto STOP;}
      }
    }
    /* 2. Possibility: one matrix for each transition class specified */
    else if (!strncmp(s->id, "Ak_", 3)) {
      if (!cos_read)  {
	scanner_error(s, "cos unknown (needed for dim(A))"); 
	goto STOP;
      }
      if (!n_read) {scanner_error(s, "need N as a range for A"); goto STOP;}
      if (a_read) {scanner_error(s, "identifier A twice"); goto STOP;}
      if (!m_calloc(a_3dmatrix, smo->cos)) {mes_proc(); goto STOP;}
      for (osc = 0; osc < smo->cos; osc++) {
	if (!m_calloc(a_3dmatrix[osc], smo->N)) {mes_proc(); goto STOP;}
	scanner_get_name(s);
	if (!strcmp(s->id, "matrix")) {
	  if (matrix_d_read(s, a_3dmatrix[osc], smo->N, smo->N)){
	    scanner_error(s, "unable to read matrix Ak"); goto STOP;
	  }
	}
	else {
	  scanner_error(s, "unknown identifier"); goto STOP;
	}
	if (osc < smo->cos-1) {
	  scanner_consume( s, ';' ); if(s->err) goto STOP;
	  /* read next matrix */
	  scanner_get_name(s);
	  if (strncmp(s->id, "Ak_", 3)) {
	    scanner_error(s, "not enough matrices Ak"); goto STOP;
	  }
	  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)){
	scanner_error(s, "need M and N as a range for C"); goto STOP;
      }
      if (c_read) {scanner_error(s, "identifier C twice"); goto STOP;}
      if (!m_calloc(c_matrix, smo->N)) {mes_proc(); goto STOP;}
      scanner_get_name(s);
      if (!strcmp(s->id, "matrix")) {
	if(matrix_d_read(s, c_matrix, smo->N, smo->M)) {
	  scanner_error(s, "unable to read matrix C"); goto STOP;
	}
	c_read = 1;
      }
      else {
	scanner_error(s, "unknown identifier"); goto STOP;
      }
    }
    else if (!strcmp(s->id, "Mue")) {/* mean of normal density */
      if ((!n_read)||(!m_read)){
	scanner_error(s, "need M and N as a range for Mue"); goto STOP;
      }
      if (mue_read) {scanner_error(s, "identifier Mue twice"); goto STOP;}
      if (!m_calloc(mue_matrix, smo->N)) {mes_proc(); goto STOP;}
      scanner_get_name(s);
      if (!strcmp(s->id, "matrix")) {
	if(matrix_d_read(s, mue_matrix, smo->N, smo->M)) {
	  scanner_error(s, "unable to read matrix Mue"); goto STOP;
	}
	mue_read = 1;
      }
      else {
	scanner_error(s, "unknown identifier"); goto STOP;
      }
    }
    else if (!strcmp(s->id, "U")) {/* variances of normal density */
      if ((!n_read)||(!m_read)){
	scanner_error(s, "need M and N as a range for U"); goto STOP;
      }
      if (u_read) {scanner_error(s, "identifier U twice"); goto STOP;}
      if (!m_calloc(u_matrix, smo->N)) {mes_proc(); goto STOP;}
      scanner_get_name(s);
      if (!strcmp(s->id, "matrix")) {
	if(matrix_d_read(s, u_matrix, smo->N, smo->M)) {
	  scanner_error(s, "unable to read matrix U"); goto STOP;
	}
	u_read = 1;
      }
      else {
	scanner_error(s, "unknown identifier"); goto STOP;
      }
    }
    else {
      scanner_error(s, "unknown identifier");
      goto STOP;
    }
    scanner_consume( s, ';' ); if(s->err) goto STOP;
  } /* while(..s->c-'}') */
  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) {
    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) {
    if (!m_calloc(fix_vektor, smo->N)) {mes_proc(); goto STOP;}
    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++) {
    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 (smodel_state_alloc(smo->s + i, smo->M, smo->s[i].in_states,
			   smo->s[i].out_states, smo->cos)) { mes_proc(); goto STOP; }
    /* copy values read to smodel */
    if(smodel_copy_vectors(smo, i, pi_vektor, fix_vektor, a_3dmatrix, c_matrix, mue_matrix,
			   u_matrix)) {mes_proc(); goto STOP;}
  }
  if (a_3dmatrix)
    for (i = 0; i < smo->cos; i++) matrix_d_free(&(a_3dmatrix[i]), smo->N);
  m_free(a_3dmatrix);
  matrix_d_free(&c_matrix, smo->N);
  matrix_d_free(&mue_matrix, smo->N);
  matrix_d_free(&u_matrix, smo->N);
  m_free(pi_vektor); m_free(fix_vektor);
  return(smo);
STOP:
  if (a_3dmatrix) 
    for (i = 0; i < smo->cos; i++) matrix_d_free(&(a_3dmatrix[i]), smo->N);
  m_free(a_3dmatrix);
  matrix_d_free(&c_matrix, smo->N);
  matrix_d_free(&mue_matrix, smo->N);
  matrix_d_free(&u_matrix, smo->N);
  m_free(pi_vektor); m_free(fix_vektor);
  smodel_free(&smo);
  return NULL;
#undef CUR_PROC
} /* smodel_read_block */


/*============================================================================*/
int smodel_free(smodel **smo) {
#define CUR_PROC "smodel_free"
  int i;
  mes_check_ptr(smo, return(-1));
  if( !*smo ) return(0);
  for (i = 0; i < (*smo)->N; i++) {
    m_free((*smo)->s[i].out_id);
    m_free((*smo)->s[i].in_id);