pviterbi_propagate.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/pviterbi_propagate.c
*       Authors:  Matthias Heinig, Janne Grunau
*
*       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 $.
*
*******************************************************************************/

/* Possible sources of errors: initialisation, pushback (the loop after) */

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

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

#include "ghmm.h"
#include "mprintf.h"
#include "mes.h"
#include "matrix.h"
#include "pmodel.h"
#include "psequence.h"
#include "pviterbi.h"
#include "pviterbi_propagate.h"
#include "ghmm_internals.h"

#define PROP_EPS 10e-6


/* --- typedefs --- */
/** Matrix cell **/
typedef struct cell {
  /** x coordinate of the cell **/
  int x;
  /** y coordinate of the cell **/
  int y;
  /** state coordinate of the cell **/
  int state;
  /** state that was inhabited before **/
  int previous_state;
  /** probability up to the previous cell **/
  double log_p;
  /** transition from the provious state to this state **/
  double log_a;
  int ref_count;
} cell;

typedef struct plocal_propagate_store_t {
  /** precomputed log probabilities for transitions into the states 
      for each transition class of the source state **/
  double *** log_in_a;
  /** precomputed log probabilities for each state for the emissions **/
  double **log_b;
  /** lookback matrix for the last offset steps **/
  double *** phi;
  /** log probabilities for the current u, v and every state **/
  double *phi_new;
  /** traceback matrix of cell pointers **/
  cell **** end_of_first;
  /** for convinience store a pointer to the model **/
  ghmm_dpmodel * mo;
  /** for the debug mode store information of matrix sizes **/
  /** length of sequence X determines size of psi **/
  int len_x;
  /** length of sequence Y determines size of phi and psi **/
  int len_y;
  /** for the recursion reuse memory here is the start index **/
  int start_x;
  int start_y;
  /** non functional stuff **/
  int    *topo_order;
  int    topo_order_length;
} plocal_propagate_store_t;


/*  --==--  local declarations --==-- */

static plocal_propagate_store_t * pviterbi_propagate_alloc (ghmm_dpmodel *mo, int len_y);
  
static int pviterbi_propagate_free (plocal_propagate_store_t **v, int n, int max_offset_x,
				    int max_offset_y, int len_y);

static cell * init_cell(int x, int y, int ghmm_dstate, int previous_state, double log_p,
			double log_a);

static void pviterbi_prop_precompute (ghmm_dpmodel *mo, plocal_propagate_store_t *pv);


static cell * pviterbi_propagate_step (ghmm_dpmodel *mo, ghmm_dpseq * X, ghmm_dpseq * Y,
				       cell * start, cell * stop, double * log_p,
				       plocal_propagate_store_t * pv);

static int * pviterbi_propagate_recursion(ghmm_dpmodel *mo, ghmm_dpseq * X, ghmm_dpseq * Y,
				   double *log_p, int *path_length, cell *start,
				   cell *stop, double max_size,
				   plocal_propagate_store_t * pv);


/*------------        Here comes the Propagate stuff          ------------- */
/*============================================================================*/
static cell * init_cell (int x, int y, int state, int previous_state,
			 double log_p, double log_a) {
#define CUR_PROC "init_cell"
  cell * mcell;
  ARRAY_CALLOC (mcell, 1);
  /* printf("Alloc cell: %i\n", sizeof(*mcell)); */
  mcell->x = x;
  mcell->y = y;
  mcell->state = state;
  mcell->previous_state = previous_state;
  mcell->log_p = log_p;
  mcell->log_a = log_a;
  return mcell;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return(NULL);
#undef CUR_PROC
}

/*============================================================================*/
cell * ghmm_dpmodel_copy_cell(cell * c) {
  if (c)
    return init_cell(c->x, c->y, c->state, c->previous_state, c->log_p, c->log_a);
  else
    return NULL;
}

/*============================================================================*/
void ghmm_dpmodel_print_cell(cell * c) {
  if (c)
    printf("(%i, %i) state: %i previous state: %i log p: %f log a: %f\n", c->x, c->y, c->state, c->previous_state, c->log_p, c->log_a);
  else
    printf("No cell\n");
}

/*============================================================================*/
static plocal_propagate_store_t * pviterbi_propagate_alloc (ghmm_dpmodel *mo, int len_y) {
#define CUR_PROC "pviterbi_propagate_alloc"
  plocal_propagate_store_t* v = NULL;
  int i, j, k;
  ARRAY_CALLOC (v, 1);

  v->mo = mo;
  v->len_y = len_y;
  /* Allocate the log_in_a's -> individal lenghts */
  ARRAY_CALLOC (v->log_in_a, mo->N);
  /* first index of log_in_a: target state */
  for (j = 0; j < mo->N; j++){ 
    /* second index: source state */
    ARRAY_CALLOC (v->log_in_a[j], mo->s[j].in_states);
    for (i=0; i<mo->s[j].in_states; i++) {
      /* third index: transition classes of source state */
      ARRAY_CALLOC (v->log_in_a[j][i], mo->s[mo->s[j].in_id[i]].kclasses);
    }
  }

  ARRAY_CALLOC (v->log_b, mo->N);
  for (j=0; j<mo->N; j++) {
    ARRAY_CALLOC (v->log_b[j], ghmm_dpmodel_emission_table_size(mo, j) + 1);
  }
  if (!(v->log_b)) {GHMM_LOG_QUEUED(LCONVERTED); goto STOP;}
  v->phi = ighmm_cmatrix_3d_alloc(mo->max_offset_x + 1, len_y + mo->max_offset_y + 1, mo->N);
  if (!(v->phi)) {GHMM_LOG_QUEUED(LCONVERTED); goto STOP;}
  ARRAY_CALLOC (v->phi_new, mo->N);
  ARRAY_CALLOC (v->end_of_first, mo->max_offset_x + 1);
  for (j=0; j<mo->max_offset_x + 1; j++) {
    ARRAY_CALLOC (v->end_of_first[j], len_y + mo->max_offset_y + 1);
    for (i=0; i<len_y + mo->max_offset_y + 1; i++) {
      ARRAY_CALLOC (v->end_of_first[j][i], mo->N);
      for (k=0; k<mo->N; k++)
	v->end_of_first[j][i][k] = NULL;
	/*ARRAY_CALLOC (v->end_of_first[j][i][k], 1);*/
    }
  }
  v->topo_order_length = 0;
  ARRAY_CALLOC (v->topo_order, mo->N);
  return(v);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  pviterbi_propagate_free(&v, mo->N, mo->max_offset_x, mo->max_offset_y, len_y);
  return(NULL);
#undef CUR_PROC
}

/*============================================================================*/
static void pviterbi_prop_precompute (ghmm_dpmodel *mo, plocal_propagate_store_t *pv)
{
#define CUR_PROC "pviterbi_precompute"
  int i, j, emission, t_class;
  
  /* Precomputing the log(a_ij) */  
  for (j = 0; j < mo->N; j++)
    for (i = 0; i < mo->s[j].in_states; i++) 
      for (t_class = 0; t_class < mo->s[mo->s[j].in_id[i]].kclasses; t_class++){
	if ( mo->s[j].in_a[i][t_class] == 0.0 )  
	  pv->log_in_a[j][i][t_class] = +1; 
	else
	  pv->log_in_a[j][i][t_class] = log( mo->s[j].in_a[i][t_class] );
      } 

  /* Precomputing the log emission probabilities for each state*/
  for (j = 0; j < mo->N; j++) {
    for (emission = 0; emission < ghmm_dpmodel_emission_table_size(mo,j); emission++) {
      if (1) {
	if ( mo->s[j].b[emission] == 0.0 )   /* DBL_EPSILON ? */ 
	  pv->log_b[j][emission] = +1; 
	else
	  pv->log_b[j][emission] = log( mo->s[j].b[emission] );
      }
      else{
	pv->log_b[j][emission] = 0.0; 
      }
    }
    pv->log_b[j][emission] = 1; /* the last field is for invalid emissions */
  }
#undef CUR_PROC
}/* viterbi_precompute */


/*============================================================================*/
static int pviterbi_propagate_free (plocal_propagate_store_t **v, int n,
				    int max_offset_x, int max_offset_y,
				    int len_y) {
#define CUR_PROC "pviterbi_propagate_free"
  int j, i;
  mes_check_ptr(v, return(-1));
  if( !*v ) return(0);
  for (j = 0; j < n; j++) {
    for (i=0; i<(*v)->mo->s[j].in_states; i++)
      m_free((*v)->log_in_a[j][i]);
    m_free((*v)->log_in_a[j]);
  }
  m_free((*v)->log_in_a);
  for (j=0; j<n; j++)
    m_free((*v)->log_b[j]);
  m_free((*v)->log_b);
  ighmm_cmatrix_3d_free(&((*v)->phi), max_offset_x + 1, len_y + max_offset_y + 1);
  m_free((*v)->phi_new);
  if ((*v)->end_of_first) {
    for (j=0; j<max_offset_x + 1; j++) {
      if ((*v)->end_of_first[j]) {
	for (i=0; i<len_y + max_offset_y + 1; i++)
	  if ((*v)->end_of_first[j][i])
	    m_free((*v)->end_of_first[j][i]);
	m_free((*v)->end_of_first[j]);
      }
    }
    m_free((*v)->end_of_first);
  }
  m_free((*v)->topo_order);
  m_free(*v);
  return(0);
#undef CUR_PROC
}

/*============================================================================*/
static double sget_log_in_a_prop(plocal_propagate_store_t * pv, int i, int j,
				 ghmm_dpseq * X, ghmm_dpseq * Y, int index_x,
				 int index_y) {
  /* determine the transition class for the source state */
  int id = pv->mo->s[i].in_id[j];
  int cl = pv->mo->s[id].class_change->get_class(pv->mo, X, Y, index_x,index_y,
					pv->mo->s[id].class_change->user_data);
  return pv->log_in_a[i][j][cl];
}

/*============================================================================*/
static double log_b_prop(plocal_propagate_store_t * pv, int state, int emission) {
#ifdef DEBUG
  if (state > pv->mo->N) 
    fprintf(stderr, "log_b_prop: State index out of bounds %i > %i\n", state, pv->mo->N);
  if (emission > ghmm_dpmodel_emission_table_size(pv->mo, state))
    fprintf(stderr, "log_b_prop: Emission index out of bounds %i > %i for state %i\n",
	    emission, ghmm_dpmodel_emission_table_size(pv->mo, state), state);
#endif
  return pv->log_b[state][emission];
}

/*============================================================================*/
static double get_phi_prop(plocal_propagate_store_t * pv, int x, int y,
			   int offset_x, int offset_y, int state){
#ifdef DEBUG
  if (y > pv->len_y || state > pv->mo->N || y < - pv->mo->max_offset_y)
    fprintf(stderr, "get_phi_prop: out of bounds %i %i %i\n", 
	    x + pv->mo->max_offset_x, y + pv->mo->max_offset_y, state);
#endif

  if (y - offset_y + pv->mo->max_offset_y >= pv->start_y)
    return pv->phi[offset_x][y - offset_y + pv->mo->max_offset_y][state];
  else
    return 1;
}

/*============================================================================*/
/* set the value of this matrix cell */
static void set_phi_prop (plocal_propagate_store_t * pv, int x, int y,
			 int state, double prob){
#ifdef DEBUG
  if (y > pv->len_y || state > pv->mo->N || y < - pv->mo->max_offset_y)
    fprintf(stderr, "set_phi_prop: out of bounds %i %i %i %f\n", 
	    x + pv->mo->max_offset_x, y + pv->mo->max_offset_y, 
	    state, prob);
#endif
  pv->phi[0][y + pv->mo->max_offset_y][state] = prob;
}

/*============================================================================*/
static void set_end_of_first (plocal_propagate_store_t * pv, int x, int y,
			     int state, cell * end_of_first){
#ifdef DEBUG
  if (y > pv->len_y || state > pv->mo->N || y < - pv->mo->max_offset_y) {
    fprintf(stderr, "set_end_of_first: out of bounds %i %i %i\n", 
	    x + pv->mo->max_offset_x, y + pv->mo->max_offset_y, 
	    state);
    ghmm_dpmodel_print_cell(end_of_first);
  }
#endif
  pv->end_of_first[0][y + pv->mo->max_offset_y][state] = end_of_first;
}

/*============================================================================*/
static cell * get_end_of_first (plocal_propagate_store_t * pv, int x, int y,
				int offset_x, int offset_y, int state) {