modelutil.c

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

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, see http://ghmm.org
*
*       Filename: ghmm/ghmm/modelutil.c
*       Authors:  Wasinee Rungsarityotin
*
*       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: 1713 $
*                       from $Date: 2006-10-16 16:06:28 +0200 (Mon, 16 Oct 2006) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

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

#include "mes.h"
#include "ghmm.h"
#include "model.h"
#include "matrix.h"
#include "ghmm_internals.h"


typedef enum eDFSCOLORS { GRAY = 0, BLACK = 1, WHITE = 2, NONE =
    -1 } DFSCOLORS;

typedef struct local_store_topo {
  int *topo_order;
  int topo_order_length;
  int *queue;
  int head, tail;
} local_store_topo;

static local_store_topo *topo_alloc (ghmm_dmodel * mo, int len);
static int topo_free (local_store_topo ** v, int n, int len);


/*----------------------------------------------------------------------------*/
static local_store_topo *topo_alloc (ghmm_dmodel * mo, int len)
{
#define CUR_PROC "sdtopo_alloc"
  local_store_topo *v = NULL;

  ARRAY_CALLOC (v, 1);
  ARRAY_CALLOC (v->queue, mo->N);

  v->topo_order_length = 0;
  v->head = 0;                  /* initialize static queue (array implementation) */
  v->tail = 0;
  ARRAY_CALLOC (v->topo_order, mo->N);

  return (v);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  topo_free (&v, mo->N, len);
  return (NULL);
#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
static int topo_free (local_store_topo ** v, int n, int len)
{
#define CUR_PROC "topo_free"

  mes_check_ptr (v, return (-1));
  if (!*v)
    return (0);
  (*v)->head = 0;
  (*v)->tail = 0;
  m_free ((*v)->topo_order);
  m_free ((*v)->queue);
  m_free (*v);
  return (0);
#undef CUR_PROC
}

/*---------------------------------------------------------------------------*/
/*
 * Implementation of DFSVisit with recursion (WS)
 *
 */
static void model_DFSVisit (ghmm_dmodel * c_model, int nodev, int *timevisit,
                            int *parents, int *colors, int **edge_classes)
{
#define CUR_PROC "DFSVisit"
  int i, w;
  colors[nodev] = GRAY;
  ++(*timevisit);
  for (i = 0; i < c_model->s[nodev].out_states; i++) {
    w = c_model->s[nodev].out_id[i];    /* Explore edge (v,w)*/
    if (edge_classes[nodev][w] == NONE) {       /* First exploration*/
      edge_classes[nodev][w] = colors[w];
      /*fprintf(stderr, " %d edge (%s, %s)\n", (int) colors[w], c_model->s[nodev].label, c_model->s[w].label); */
    }
    if (colors[w] == WHITE) {
      parents[w] = nodev;
      model_DFSVisit (c_model, w, timevisit, parents, colors, edge_classes);
    }
  }
  colors[nodev] = BLACK;        /* finished  */
  ++(*timevisit);
#undef CUR_PROC
}


/** Return classification of edges in the model graph 
    WHITE EDGE => edges in the DFS search tree
    GRAY EDGE  => edges that form cycles which must be removed 
                  before running topological sort
*/
int **ghmm_dmodel_DFS (ghmm_dmodel * c_model)
{
#define CUR_PROC "ghmm_dmodel_DFS"
  int i, j, initials=0;
  int timevisit = 0;
  int *colors;
  int *parents=NULL;
  int **edge_classes=NULL;

  ARRAY_CALLOC (colors, c_model->N);
  ARRAY_CALLOC (parents, c_model->N);

  /*edge_classes=(int**)calloc(c_model->N,sizeof(int*));*/
  /*for(i=0; i < c_model->N; i++) {*/
  /*  edge_classes[i]=(int*)calloc(c_model->N,sizeof(int));*/
  /*}*/


  edge_classes = ighmm_dmatrix_stat_alloc (c_model->N, c_model->N);


  for (i = 0; i < c_model->N; i++) {
    if (c_model->s[i].pi == 1.0)
      initials = i;             /* assuming only one initial state */
    colors[i] = WHITE;
    parents[i] = -1;
    for (j = 0; j < c_model->N; j++) {
      edge_classes[i][j] = NONE;
    }
  }

  model_DFSVisit (c_model, initials, &timevisit, parents, colors,
                  edge_classes);
  for (i = 0; i < c_model->N; i++) {    /* handle subtrees */
    if (colors[i] == WHITE) {
      model_DFSVisit (c_model, i, &timevisit, parents, colors, edge_classes);
    }
  }

STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  m_free (colors);
  m_free (parents);
  return edge_classes;
#undef CUR_PROC
}

static void __topological_sort (ghmm_dmodel * c_model, local_store_topo * v,
                                int **edge_classes)
{
  int i, j;
  int nodeu, nodew, dels_cnt;
  int *indegrees = (int *) calloc (c_model->N, sizeof (int));

  for (i = 0; i < c_model->N; i++) {
    indegrees[i] = c_model->s[i].in_states;
  }
  for (i = 0; i < c_model->N; i++) {    /* don't consider back edges in top sort */
    for (j = 0; j < c_model->N; j++) {
      if (edge_classes[i][j] == GRAY) {
        indegrees[j]--;
        /* fprintf(stderr, "BACK edge (%s, %s)\n", c_model->s[i].label, c_model->s[j].label); */
      }
    }
  }
  /* Create a queue q and enqueue all nodes of indegree 0. */
  v->head = 0;
  v->tail = 0;
  for (i = 0; i < c_model->N; i++) {
    if (indegrees[i] == 0)
      v->queue[v->tail++] = i;
  }
  dels_cnt = 0;
  while (v->head != v->tail) {
    nodeu = v->queue[v->head++];        /* dequeue */
    if (c_model->silent[nodeu]) {
      v->topo_order[dels_cnt++] = nodeu;        /* append it to the list */
      /* printf("Silent state: %s\n", c_model->s[nodeu].label); */
    }
    for (i = 0; i < c_model->s[nodeu].out_states; i++) {
      nodew = c_model->s[nodeu].out_id[i];
      if (edge_classes[nodeu][nodew] != GRAY) {
        indegrees[nodew]--;
        if (nodew != nodeu && indegrees[nodew] == 0) {
          v->queue[v->tail++] = nodew;  /* enqueue */
        }
      }
    }
  }
  v->topo_order_length = dels_cnt;
  free (indegrees);
}

/*----------------------------------------------------------------------------*/
void ghmm_dmodel_order_topological (ghmm_dmodel * mo)
{
#define CUR_PROC "ghmm_dmodel_topo_order"
  int i;
  local_store_topo *v;
  int **edge_cls;

  v = topo_alloc (mo, 1);
  if (!v) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }

  edge_cls = ghmm_dmodel_DFS (mo);
  __topological_sort (mo, v, edge_cls);
  mo->topo_order_length = v->topo_order_length;

  /* free topo_order if existing. XXX maybe safe to skip sorting if topo_order
     exists? Can a model change topological? */
  if (mo->topo_order)
    m_free(mo->topo_order);
  ARRAY_CALLOC (mo->topo_order, mo->topo_order_length);

  for (i = 0; i < v->topo_order_length; i++) {
    mo->topo_order[i] = v->topo_order[i];
  }

  /*fprintf(stderr,"Ordering silent states....\n\t");
     for(i=0; i < mo->topo_order_length; i++) {
     fprintf(stderr, "%d, ", mo->topo_order[i]);
     }
     fprintf(stderr,"\n\n");  

     for(i=0; i < mo->N; i++) { 
     for(j=0; j < mo->N; j++) {
     fprintf(stderr,"edge_cls[%d][%d]=%d\n",i,j,edge_cls[i][j]);           
     }
     } */

  ighmm_dmatrix_stat_free (&edge_cls);
  topo_free (&v, mo->N, 1);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  i = 0;
#undef CUR_PROC
}