smodel.c

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

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, see http://ghmm.org
*
*       Filename: ghmm/ghmm/smodel.c
*       Authors:  Bernhard Knab, Benjamin Georgi
*
*       Copyright (C) 1998-2004 Alexander Schliep
*       Copyright (C) 1998-2001 ZAIK/ZPR, Universitaet zu Koeln
*       Copyright (C) 2002-2004 Max-Planck-Institut fuer Molekulare Genetik,
*                               Berlin
*
*       Contact: [e-mail]
*
*       This library is free software; you can redistribute it and/or
*       modify it under the terms of the GNU Library General Public
*       License as published by the Free Software Foundation; either
*       version 2 of the License, or (at your option) any later version.
*
*       This library 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
*       Library General Public License for more details.
*
*       You should have received a copy of the GNU Library General Public
*       License along with this library; if not, write to the Free
*       Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*
*       This file is version $Revision: 1831 $
*                       from $Date: 2007-05-23 14:18:45 +0200 (Wed, 23 May 2007) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

#ifdef WIN32
#  include "win_config.h"
#endif

#ifdef HAVE_CONFIG_H
#  include "../config.h"
#endif

#include <float.h>
#include <math.h>
#include <assert.h>

#include "ghmm.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 "rng.h"
#include "randvar.h"
#include "string.h"
#include "ghmm_internals.h"
#include "xmlreader.h"
#include "xmlwriter.h"

#include "obsolete.h"

/*----------------------------------------------------------------------------*/
int ghmm_cstate_alloc (ghmm_cstate * s,
                        int M, int in_states, int out_states, int cos)
{
# define CUR_PROC "ghmm_cstate_alloc"
  int res = -1;
  int i;
  ARRAY_CALLOC (s->c, M);
  ARRAY_CALLOC (s->mue, M);
  ARRAY_CALLOC (s->u, M);
  ARRAY_CALLOC (s->a,M);

  ARRAY_CALLOC (s->mixture_fix, M);
  ARRAY_CALLOC (s->density,M);
  s->xPosition = 0.0;
  s->yPosition = 0.0;
   

  /* mixture component fixing deactivated by default */
  for (i = 0; i < M; i++) {
    s->mixture_fix[i] = 0;
  }

  if (out_states > 0) {
    ARRAY_CALLOC (s->out_id, out_states);
    s->out_a = ighmm_cmatrix_alloc (cos, out_states);
    if (!s->out_a) {
      GHMM_LOG_QUEUED(LCONVERTED);
      goto STOP;
    }
  }
  if (in_states > 0) {
    ARRAY_CALLOC (s->in_id, in_states);
    s->in_a = ighmm_cmatrix_alloc (cos, in_states);
    if (!s->in_a) {
      GHMM_LOG_QUEUED(LCONVERTED);
      goto STOP;
    }
  }
  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return (res);
# undef CUR_PROC
}                               /* ghmm_cstate_alloc */

/*----------------------------------------------------------------------------*/
ghmm_cmodel * ghmm_cmodel_calloc(int N, int modeltype) {
#define CUR_PROC "ghmm_cmodel_calloc"
  ghmm_cmodel * mo;

  assert(modeltype & GHMM_kContinuousHMM);
  ARRAY_CALLOC(mo, 1);
  mo->N = N;
  mo->M = 0;
  mo->model_type = modeltype;
  ARRAY_CALLOC(mo->s, N);
  return mo;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ghmm_cmodel_free(&mo);
  return NULL;
#undef CUR_PROC
}

int ghmm_cmodel_class_change_alloc (ghmm_cmodel * smo)
{
#define CUR_PROC "ghmm_cmodel_class_change_alloc"
  ghmm_cmodel_class_change_context *c = NULL;
  ARRAY_CALLOC (c, 1);

  c->k = -1;
  c->get_class = NULL;

  smo->class_change = c;
  return (0);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return (-1);
# undef CUR_PROC
}

#ifdef GHMM_OBSOLETE
/*----------------------------------------------------------------------------*/
static int smodel_copy_vectors (ghmm_cmodel * 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"
/* XXX - this function is only used in the smo file reading functions. It should be removed 
since it does not take care of het .mixtures densities  */
  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->s[index].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) {
        GHMM_LOG_QUEUED(LCONVERTED);
        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) {
        GHMM_LOG_QUEUED(LCONVERTED);
        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 */


/*============================================================================*/
ghmm_cmodel **ghmm_cmodel_read (const char *filename, int *smo_number)
{
#define CUR_PROC "ghmm_cmodel_read"
/* XXX - old function. No support to heterogeneous densities  */
  int j;
  long new_models = 0;
  scanner_t *s = NULL;
  ghmm_cmodel **smo = NULL;
  *smo_number = 0;
  s = ighmm_scanner_alloc (filename);
  if (!s) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  while (!s->err && !s->eof) {
    ighmm_scanner_get_name (s);
    ighmm_scanner_consume (s, '=');
    if (s->err)
      goto STOP;
    if (!strcmp (s->id, "SHMM") || !strcmp (s->id, "shmm")) {
      (*smo_number)++;
      /* more mem */
      ARRAY_REALLOC (smo, *smo_number);
      smo[*smo_number - 1] = ghmm_cmodel_read_block (s, (int *) &new_models);
      if (!smo[*smo_number - 1]) {
        GHMM_LOG_QUEUED(LCONVERTED);
        goto STOP;
      }
      /* copy ghmm_cmodel */
      if (new_models > 1) {
        /* "-1" due to  (*smo_number)++ from above  */
        ARRAY_REALLOC (smo, *smo_number - 1 + new_models);
        for (j = 1; j < new_models; j++) {
          smo[*smo_number] = ghmm_cmodel_copy (smo[*smo_number - 1]);
          if (!smo[*smo_number]) {
            GHMM_LOG_QUEUED(LCONVERTED);
            goto STOP;
          }
          (*smo_number)++;
        }
      }
    }
    else {
      ighmm_scanner_error (s, "unknown identifier");
      goto STOP;
    }
    ighmm_scanner_consume (s, ';');
    if (s->err)
      goto STOP;
  }                             /* while(!s->err && !s->eof) */
  /*if (ghmm_cmodel_check_compatibility(smo, *smo_number) == -1) {
     GHMM_LOG_QUEUED(LCONVERTED); goto STOP;
     } */

  ighmm_scanner_free (&s);             
  return smo;
  
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
    ighmm_scanner_free (&s);
    return NULL;
#undef CUR_PROC
}    /* ghmm_cmodel_read */


/*============================================================================*/
ghmm_cmodel *ghmm_cmodel_read_block (scanner_t * s, int *multip)
{
#define CUR_PROC "ghmm_cmodel_read_block"
/* XXX - old function. No support to heterogeneous densities  */
  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, M;
  ghmm_cmodel *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;
  ghmm_density_t  density = 0;

  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 */
  ARRAY_CALLOC (smo, 1);

  ighmm_scanner_consume (s, '{');
  if (s->err)
    goto STOP;
  while (!s->err && !s->eof && s->c - '}') {
    ighmm_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)) {
      ighmm_scanner_error (s, "unknown identifier");
      goto STOP;
    }
    ighmm_scanner_consume (s, '=');
    if (s->err)
      goto STOP;
    if (!strcmp (s->id, "multip")) {
      *multip = ighmm_scanner_get_int (s);
      if (*multip < 1) {        /* ignore: makes no sense */
        *multip = 1;
        GHMM_LOG(LCONVERTED, "Multip < 1 ignored\n");
      }
    }
    else if (!strcmp (s->id, "M")) {    /* number of output components */
      if (m_read) {
        ighmm_scanner_error (s, "identifier M twice");
        goto STOP;
      }
      M = ighmm_scanner_get_int (s);
      m_read = 1;
    }
    else if (!strcmp (s->id, "N")) {    /* number of states */
      if (n_read) {
        ighmm_scanner_error (s, "identifier N twice");
        goto STOP;
      }
      smo->N = ighmm_scanner_get_int (s);
      ARRAY_CALLOC (smo->s, smo->N);
      /*ARRAY_CALLOC (smo->density, smo->N);*/
      n_read = 1;
    }
    else if (!strcmp (s->id, "density")) {      /* which density function? */
      if (density_read) {
        ighmm_scanner_error (s, "identifier density twice");
        goto STOP;
      }
      /*current smo format only suport models with same densities*/      
      density = (ghmm_density_t) ighmm_scanner_get_int (s);
      if ((int) density < 0 || density >= density_number) {
        ighmm_scanner_error (s, "unknown typ of density function");
        goto STOP;
      }
      density_read = 1;
    }
    else if (!strcmp (s->id, "prior")) {        /* modelprior */
      if (prior_read) {
        ighmm_scanner_error (s, "identifier prior twice");
        goto STOP;
      }
      smo->prior = ighmm_scanner_get_edouble (s);
      if ((smo->prior < 0 || smo->prior > 1) && smo->prior != -1) {
        ighmm_scanner_error (s, "invalid model prior");
        goto STOP;
      }
      prior_read = 1;
    }
    else if (!strcmp (s->id, "cos")) {  /* number of transition classes */
      if (cos_read) {
        ighmm_scanner_error (s, "identifier cos twice");
        goto STOP;
      }
      smo->cos = ighmm_scanner_get_int (s);
      if (smo->cos <= 0) {
        ighmm_scanner_error (s, "invalid model cos");
        goto STOP;
      }
      cos_read = 1;
    }
    else if (!strcmp (s->id, "Pi")) {   /* Initial State Prob. */
      if (!n_read) {
        ighmm_scanner_error (s, "need N as a range for Pi");
        goto STOP;
      }
      if (pi_read) {
        ighmm_scanner_error (s, "identifier Pi twice");
        goto STOP;
      }
      ighmm_scanner_get_name (s);
      if (!strcmp (s->id, "vector")) {
        ighmm_scanner_consume (s, '{');
        if (s->err)
          goto STOP;
        pi_vektor = scanner_get_double_earray (s, &len);
        if (len != smo->N) {
          ighmm_scanner_error (s, "wrong number of elements in PI");
          goto STOP;
        }
        ighmm_scanner_consume (s, ';');
        if (s->err)
          goto STOP;
        ighmm_scanner_consume (s, '}');
        if (s->err)
          goto STOP;
        pi_read = 1;
      }
      else {
        ighmm_scanner_error (s, "unknown identifier");
        goto STOP;
      }
    }
    else if (!strcmp (s->id, "fix_state")) {
      if (!n_read) {
        ighmm_scanner_error (s, "need N as a range for fix_state");
        goto STOP;
      }
      if (fix_read) {
        ighmm_scanner_error (s, "identifier fix_state twice");
        goto STOP;
      }
      ighmm_scanner_get_name (s);
      if (!strcmp (s->id, "vector")) {
        ighmm_scanner_consume (s, '{');
        if (s->err)
          goto STOP;
        fix_vektor = scanner_get_int_array (s, &len);
        if (len != smo->N) {
          ighmm_scanner_error (s, "wrong number of elements in fix_state");
          goto STOP;
        }
        ighmm_scanner_consume (s, ';');
        if (s->err)
          goto STOP;
        ighmm_scanner_consume (s, '}');
        if (s->err)
          goto STOP;
        fix_read = 1;
      }
      else {
        ighmm_scanner_error (s, "unknown identifier");
        goto STOP;
      }
    }
    /* 1. Possibility: one matrix for all transition classes */
    else if (!strcmp (s->id, "A")) {
      if (!cos_read) {
        ighmm_scanner_error (s, "cos unknown (needed for dim(A))");
        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_matrix, smo->N);
      ighmm_scanner_get_name (s);
      if (!strcmp (s->id, "matrix")) {
        if (ighmm_cmatrix_read (s, a_matrix, smo->N, smo->N)) {
          ighmm_scanner_error (s, "unable to read matrix A");
          goto STOP;
        }
        a_read = 1;
      }
      else {
        ighmm_scanner_error (s, "unknown identifier");
        goto STOP;
      }
      /* copy transition matrix to all transition classes */
      ARRAY_CALLOC (a_3dmatrix, smo->cos);
      a_3dmatrix[0] = a_matrix;
      for (i = 1; i < smo->cos; i++) {
        a_3dmatrix[i] = ighmm_cmatrix_alloc_copy (smo->N, smo->N, a_matrix);
        if (!a_3dmatrix[i]) {
          GHMM_LOG_QUEUED(LCONVERTED);
          goto STOP;
        }
      }
    }
    /* 2. Possibility: one matrix for each transition class specified */
    else if (!strncmp (s->id, "Ak_", 3)) {
      if (!cos_read) {
        ighmm_scanner_error (s, "cos unknown (needed for dim(A))");