/*******************************************************************************
*
* This file is part of the General Hidden Markov Model Library,
* GHMM version __VERSION__, 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: schliep@ghmm.org
*
* 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 10:06:28 -0400 (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
}