sreestimate.c

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

/*******************************************************************************
  author       : Bernhard Knab
  filename     : /zpr/bspk/src/hmm/ghmm/ghmm/sreestimate.c
  created      : TIME: 17:11:14     DATE: Mon 15. November 1999
  $Id: sreestimate.c,v 1.12 2003/12/19 15:06:17 cic99 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 <math.h>
#include <float.h>
#include "mprintf.h"
#include "mes.h"
#include "sreestimate.h"
#include "smodel.h"
#include "sfoba.h"
#include "matrix.h"
#include "vector.h"
#include "randvar.h"
#include "gauss_tail.h"
#include "const.h"
#include "root_finder.h"

/* switch: turn off (0)  MESCONTR and MESINFO */
#define MCI 1
/* set control output (MES_WIN, MES_FILE, MES_FILE_WIN, ...) */
#define MESCONTR MES_FILE
/* set info output  (logP, ...) */
#define MESINFO MES_FILE


/** needed for normaldensitypos (truncated normal density) */
#define ACC 1E-8
static double C_PHI;
static double CC_PHI;
/***/

static local_store_t *sreestimate_alloc(const smodel *smo);
static int sreestimate_free(local_store_t **r, int N);
static int sreestimate_init(local_store_t *r, const smodel *smo);
static int sreestimate_alloc_matvek(double ***alpha, double ***beta, 
				    double **scale, double ****b, 
				    int T, int N, int M);
static int sreestimate_precompute_b(smodel *smo, double *O, int T, double ***b);
static int sreestimate_free_matvec(double **alpha, double **beta, double *scale,
				   double ***b, int T, int N);
static int sreestimate_setlambda(local_store_t *r, smodel *smo);
static int sreestimate_one_step(smodel *smo, local_store_t *r, 
			       int seq_number, int *T, double **O,
			       double *log_p, double *seq_w);
/*----------------------------------------------------------------------------*/
/* various allocations */
static local_store_t *sreestimate_alloc(const smodel *smo) {
# define CUR_PROC "sreestimate_alloc"
  int i;
  local_store_t* r = NULL;
  if (!m_calloc(r, 1)) {mes_proc(); goto STOP;}
  r->cos=smo->cos;
  if (!m_calloc(r->pi_num, smo->N)) {mes_proc(); goto STOP;}
  if (!m_calloc(r->a_num, smo->N)) {mes_proc(); goto STOP;}
  for (i = 0; i < smo->N; i++) {
    r->a_num[i] = stat_matrix_d_alloc(smo->cos, smo->s[i].out_states);
    if (!r->a_num[i]) {mes_proc(); goto STOP;}
  }
  r->a_denom = stat_matrix_d_alloc(smo->N, smo->cos);
  if (!r->a_denom) {mes_proc(); goto STOP;}
  /***/
  if (!m_calloc(r->c_denom, smo->N)) {mes_proc(); goto STOP;}
  r->c_num = stat_matrix_d_alloc(smo->N, smo->M);
  if (!(r->c_num)) {mes_proc(); goto STOP;}
  r->mue_num = stat_matrix_d_alloc(smo->N, smo->M);
  if (!(r->mue_num)) {mes_proc(); goto STOP;}
  r->u_num = stat_matrix_d_alloc(smo->N, smo->M);
  if (!(r->u_num)) {mes_proc(); goto STOP;}
  r->mue_u_denom = stat_matrix_d_alloc(smo->N, smo->M);
  if (!(r->mue_u_denom)) {mes_proc(); goto STOP;}
  r->sum_gt_otot = stat_matrix_d_alloc(smo->N, smo->M);
  if (!(r->sum_gt_otot)) {mes_proc(); goto STOP;}
  r->sum_gt_logb = stat_matrix_d_alloc(smo->N, smo->M);
  if (!(r->sum_gt_logb)) {mes_proc(); goto STOP;}
  return(r);
STOP:
  sreestimate_free(&r, smo->N);
  return 0;
# undef CUR_PROC
} /* sreestimate_alloc */

/*----------------------------------------------------------------------------*/
static int sreestimate_free(local_store_t **r, int N) {
# define CUR_PROC "sreestimate_free"
  int i;
  mes_check_ptr(r, return(-1));
  if( !*r ) return(0);
  m_free((*r)->pi_num);  
  for (i = 0; i < N; i++)
    stat_matrix_d_free( &((*r)->a_num[i]));
  m_free((*r)->a_num);
  stat_matrix_d_free( &((*r)->a_denom));
  /***/
  m_free((*r)->c_denom);
  stat_matrix_d_free( &((*r)->c_num) );
  stat_matrix_d_free( &((*r)->mue_num) );
  stat_matrix_d_free( &((*r)->u_num) );
  stat_matrix_d_free( &((*r)->mue_u_denom) );
  stat_matrix_d_free( &((*r)->sum_gt_otot) );
  stat_matrix_d_free( &((*r)->sum_gt_logb) );
  m_free(*r);
  return(0);
# undef CUR_PROC
} /* sreestimate_free */

/*----------------------------------------------------------------------------*/
static int sreestimate_init(local_store_t *r, const smodel *smo) {
# define CUR_PROC "sreestimate_init"
  int i, j, m, osc;
  r->pi_denom = 0.0;
  for (i = 0; i < smo->N; i++) {
    r->pi_num[i] = 0.0;
    for (osc = 0; osc < smo->cos; osc++) {
      r->a_denom[i][osc] = 0.0;
      for (j = 0; j < smo->s[i].out_states; j++)
	r->a_num[i][osc][j] = 0.0;
    }
    r->c_denom[i] = 0.0;
    for (m = 0; m < smo->M; m++) {
      r->c_num[i][m] = 0.0;
      r->mue_num[i][m] = 0.0;
      r->u_num[i][m] = 0.0;
      r->mue_u_denom[i][m] = 0.0;
      r->sum_gt_otot[i][m] = 0.0;
      r->sum_gt_logb[i][m] = 0.0;
    }
  }
  return(0);
# undef CUR_PROC
} /* sreestimate_init */

/*----------------------------------------------------------------------------*/
static int sreestimate_alloc_matvek(double ***alpha, double ***beta, 
				    double **scale, double ****b, 
				    int T, int N, int M) {
# define CUR_PROC "sreestimate_alloc_matvek"
  int t, res = -1;
  *alpha = stat_matrix_d_alloc(T, N);
  if (!(*alpha)) {mes_proc(); goto STOP;}
  *beta = stat_matrix_d_alloc(T, N);
  if (!(*beta)) {mes_proc(); goto STOP;}
  if (!m_calloc(*scale, T)) {mes_proc(); goto STOP;}
  /* 3-dim. matrix for b[t][i][m] with m = 1..M(!): */
  if (!m_calloc(*b, T)) {mes_proc(); goto STOP;}
  for (t = 0; t < T; t++) {
    (*b)[t] = stat_matrix_d_alloc(N, M+1);
    if (!((*b)[t])) {mes_proc(); goto STOP;}
  }
  res = 0;
STOP:
  return(res);
# undef CUR_PROC
} /* sreestimate_alloc_matvek */

/*----------------------------------------------------------------------------*/
static int sreestimate_free_matvec(double **alpha,double **beta,double *scale,
				   double ***b, int T, int N) {
# define CUR_PROC "sreestimate_free_matvec"
  int t;
  stat_matrix_d_free(&alpha);
  stat_matrix_d_free(&beta);
  m_free(scale); 
  if (!b) return(0);
  for (t = 0; t < T; t++)
    stat_matrix_d_free(&b[t]);
  m_free(b);
  return(0);
# undef CUR_PROC
} /* sreestimate_free_matvec */   

/*----------------------------------------------------------------------------*/
static int sreestimate_precompute_b(smodel *smo, double *O, int T, double ***b){
# define CUR_PROC "sreestimate_precompute_b"
  int t, i, m;
  /* save sum (c_im * b_im(O_t))  in b[t][i][smo->M] */
  for (t = 0; t < T; t++)
    for (i = 0; i < smo->N; i++)
      b[t][i][smo->M] = 0.0;
  /* save c_im * b_im(O_t)  directly in  b[t][i][m] */
  for (t = 0; t < T; t++)
    for (i = 0; i < smo->N; i++)
      for (m = 0; m < smo->M; m++) {
	b[t][i][m] = smodel_calc_cmbm(smo, i, m, O[t]); 
	b[t][i][smo->M] += b[t][i][m];
      }
  return(0);
# undef CUR_PROC
} /* sreestimate_precompute_b */   


/*----------------------------------------------------------------------------*/
static int sreestimate_setlambda(local_store_t *r, smodel *smo) {
# define CUR_PROC "sreestimate_setlambda"
  int res = -1;
  int i, j, m, l, j_id, osc;
  double pi_factor, a_factor_i = 0.0, c_factor_i = 0.0, u_im, mue_im, 
    mue_left, mue_right, A, B, Atil, Btil; /* Q; */
  int a_num_pos, a_denom_pos, c_denom_pos, c_num_pos;
  char *str;

  double p_i;
  int g, h;

  if (r->pi_denom <= DBL_MIN)
    { mes_prot("pi: denominator == 0.0!\n");  goto STOP; } 
  else
    pi_factor = 1/r->pi_denom;

  for (i = 0; i < smo->N; i++) {
    /* Pi */
    smo->s[i].pi =  r->pi_num[i] * pi_factor;
    
    /* A */
    for (osc = 0; osc < smo->cos; osc++) {
      /* note: denom. might be 0; never reached state? */
      a_denom_pos = 1;
      if (r->a_denom[i][osc] <= DBL_MIN) {
	a_denom_pos = 0;
#if MCI
	if (smo->s[i].out_states > 0) { 
	  if (smo->s[i].in_states == 0)
	    mes(MESCONTR,"state %d: no in_states\n",i);
	  else {
	    /* Test: sum all ingoing == 0 ? */
	    p_i = smo->s[i].pi;
	    for (g = 0; g < smo->cos; g++)
	      for (h = 0; h < smo->s[i].in_states; h++)
		p_i += smo->s[i].in_a[g][h];
	    if (p_i == 0)
	      mes(MESCONTR,"state %d: P(in) = 0\n",i);
	    else
	      mes(MESCONTR,"state %d can be reached but a-denom. = 0.0\n",i);
	  }
	}
#endif
      }
      else
	a_factor_i = 1/r->a_denom[i][osc];
      
      a_num_pos = 0;
      for (j = 0; j < smo->s[i].out_states; j++) {
	j_id = smo->s[i].out_id[j];
	/* TEST: denom. < numerator */
	if ((r->a_denom[i][osc] - r->a_num[i][osc][j]) < -EPS_PREC) {
#if MCI
	  mes(MESCONTR,"a[%d][%d][%d]: numerator > denom.!\n", i, osc, j_id); 
#endif
	  smo->s[i].out_a[osc][j] = 1.0;
	}
	else if (a_denom_pos)
	  smo->s[i].out_a[osc][j] = r->a_num[i][osc][j] * a_factor_i;
	else 	  
	  smo->s[i].out_a[osc][j] = 0.0;
	
	if (r->a_num[i][osc][j] > 0.0) /* >= EPS_PREC ? */ 
	  a_num_pos = 1;  
      
	/* important: also update in_a  */
	l = 0;
	while (l < smo->s[j_id].in_states)
	  if  (smo->s[j_id].in_id[l] == i)
	    break;
	  else 
	    l++;
	if ( l == smo->s[j_id].in_states) { 
	  str = mprintf(NULL, 0, "no matching in_a for out_a(=a[%d][%d]) found!\n",i,j_id); 
	  mes_prot(str);
	  m_free(str);
	  goto STOP;
	}
	smo->s[j_id].in_a[osc][l] = smo->s[i].out_a[osc][j]; 
      } /* j-loop */
      
#if MCI
      if (!a_num_pos && smo->s[i].out_states > 0)
	mes(MESCONTR,"all numerators a[%d][%d][j]==0 (denom. = %.4f, P(in)=%.4f)!\n", i, osc, r->a_denom[i][osc], p_i);
#endif
    } /* osc-loop */
    
    
    /* C, Mue und U */

    /* if fix, continue to next state */
    if (smo->s[i].fix)
      continue;

    c_denom_pos = 1;
    if (r->c_denom[i] <= DBL_MIN) { /* < EPS_PREC ? */
#if MCI
      mes(MESCONTR,"c[%d][m]: denominator == 0.0!\n", i);
#endif
      c_denom_pos = 0;
    }
    else 
      c_factor_i = 1/r->c_denom[i];
    
    c_num_pos = 0;
    for (m = 0; m < smo->M; m++) {

      /* TEST: denom. < numerator */
      if ((r->c_denom[i] - r->c_num[i][m]) < 0.0) { /* < -EPS_PREC ? */
#if MCI
	mes(MESCONTR,"c[%d][%d]: numerator > denominator! (%.4f > %.4f)!\n", 
	    i, m, r->c_num[i][m], r->c_denom[i]);
#endif
	smo->s[i].c[m] = 1.0;
      }
      else if (c_denom_pos) 
	/* c_denom == 0: no change in c_im (?) */
	smo->s[i].c[m] = r->c_num[i][m] * c_factor_i;
      
      if (r->c_num[i][m] > 0.0)
	c_num_pos = 1;
      
      if ( fabs(r->mue_u_denom[i][m]) <= DBL_MIN) /* < EPS_PREC ? */
#if MCI
	mes(MESCONTR,"mue[%d][%d]: denominator == 0.0!\n", i, m);
#else
      ;
#endif