sreestimate.c

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

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, see http://ghmm.org
*
*       Filename: ghmm/ghmm/sreestimate.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: 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 <math.h>
#include <float.h>
#include "ghmm.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 "ghmm_internals.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

typedef struct local_store_t {
  int cos;
  double *pi_num;
  double pi_denom;
  double ***a_num;
  double **a_denom;
  double **c_num;
  double *c_denom;
  double **mue_num;
  double **u_num;
  double **mue_u_denom;       /* mue-denom. = u-denom. for sym. normal density */
  double **sum_gt_otot;       /* for truncated normal density */
  double **sum_gt_logb;       /* Control according to Q-function */
} local_store_t;

/** 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 ghmm_cmodel * smo);
static int sreestimate_free (local_store_t ** r, int N);
static int sreestimate_init (local_store_t * r, const ghmm_cmodel * 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 (ghmm_cmodel * 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, ghmm_cmodel * smo);
static int sreestimate_one_step (ghmm_cmodel * 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 ghmm_cmodel * smo)
{
# define CUR_PROC "sreestimate_alloc"
  int i;
  local_store_t *r = NULL;
  ARRAY_CALLOC (r, 1);
  r->cos = smo->cos;
  ARRAY_CALLOC (r->pi_num, smo->N);
  ARRAY_CALLOC (r->a_num, smo->N);
  for (i = 0; i < smo->N; i++) {
    r->a_num[i] = ighmm_cmatrix_stat_alloc (smo->cos, smo->s[i].out_states);
    if (!r->a_num[i]) {
      GHMM_LOG_QUEUED(LCONVERTED);
      goto STOP;
    }
  }
  r->a_denom = ighmm_cmatrix_stat_alloc (smo->N, smo->cos);
  if (!r->a_denom) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  /* For sacke of simplicity, a NXmax(M) is being allocated,
  even though not all emissins components are of size max(M) */
  ARRAY_CALLOC (r->c_denom, smo->N);  
  r->c_num = ighmm_cmatrix_stat_alloc (smo->N, smo->M);
  if (!(r->c_num)) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  r->mue_num = ighmm_cmatrix_stat_alloc (smo->N, smo->M);
  if (!(r->mue_num)) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  r->u_num = ighmm_cmatrix_stat_alloc (smo->N, smo->M);
  if (!(r->u_num)) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  r->mue_u_denom = ighmm_cmatrix_stat_alloc (smo->N, smo->M);
  if (!(r->mue_u_denom)) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  r->sum_gt_otot = ighmm_cmatrix_stat_alloc (smo->N, smo->M);
  if (!(r->sum_gt_otot)) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  r->sum_gt_logb = ighmm_cmatrix_stat_alloc (smo->N, smo->M);
  if (!(r->sum_gt_logb)) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  return (r);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  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++)
    ighmm_cmatrix_stat_free (&((*r)->a_num[i]));
  m_free ((*r)->a_num);
  ighmm_cmatrix_stat_free (&((*r)->a_denom));
  /***/
  m_free ((*r)->c_denom);
  ighmm_cmatrix_stat_free (&((*r)->c_num));
  ighmm_cmatrix_stat_free (&((*r)->mue_num));
  ighmm_cmatrix_stat_free (&((*r)->u_num));
  ighmm_cmatrix_stat_free (&((*r)->mue_u_denom));
  ighmm_cmatrix_stat_free (&((*r)->sum_gt_otot));
  ighmm_cmatrix_stat_free (&((*r)->sum_gt_logb));
  m_free (*r);
  return (0);
# undef CUR_PROC
}                               /* sreestimate_free */

/*----------------------------------------------------------------------------*/
static int sreestimate_init (local_store_t * r, const ghmm_cmodel * 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 = ighmm_cmatrix_stat_alloc (T, N);
  if (!(*alpha)) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  *beta = ighmm_cmatrix_stat_alloc (T, N);
  if (!(*beta)) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  ARRAY_CALLOC (*scale, T);
  /* 3-dim. matrix for b[t][i][m] with m = 1..M(!): */
  ARRAY_CALLOC (*b, T);
  for (t = 0; t < T; t++) {
    (*b)[t] = ighmm_cmatrix_stat_alloc (N, M + 1);
    if (!((*b)[t])) {
      GHMM_LOG_QUEUED(LCONVERTED);
      goto STOP;
    }
  }
  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  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;
  ighmm_cmatrix_stat_free (&alpha);
  ighmm_cmatrix_stat_free (&beta);
  m_free (scale);
  if (!b)
    return (0);
  for (t = 0; t < T; t++)
    ighmm_cmatrix_stat_free (&b[t]);
  m_free (b);
  return (0);
# undef CUR_PROC
}                               /* sreestimate_free_matvec */

/*----------------------------------------------------------------------------*/
static int sreestimate_precompute_b (ghmm_cmodel * 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->s[i].M; m++) {
        b[t][i][m] = ghmm_cmodel_calc_cmbm (smo, i, m, O[t]);
        b[t][i][smo->s[i].M] += b[t][i][m];
      }
  return (0);
# undef CUR_PROC
}                               /* sreestimate_precompute_b */


/*----------------------------------------------------------------------------*/
static int sreestimate_setlambda (local_store_t * r, ghmm_cmodel * smo)
{
# define CUR_PROC "sreestimate_setlambda"
  int res = -1;
  int i, j, m, l, j_id, osc, fix_flag;
  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, fix_w, unfix_w;    /* Q; */
  int a_num_pos, a_denom_pos, c_denom_pos, c_num_pos;
  char *str;

  double p_i=0.0;
  int g, h;

  if (r->pi_denom <= DBL_MIN) {
    GHMM_LOG(LCONVERTED, "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)
            ighmm_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)
              ighmm_mes (MESCONTR, "state %d: P(in) = 0\n", i);
            else
              ighmm_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]) < -GHMM_EPS_PREC) {
#if MCI
          ighmm_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 {
          /* if(smo->s[i].out_states == 1 &&  smo->s[i].out_id[0]==i) {
            printf("%d is an absorbing end state -> don't change transitions !\n",i);
          }   */          
#if MCI
          ighmm_mes (MESCONTR, "state %d has no observed outgoing transitions, parameters will not be reestimated. \n", i );   
#endif
          continue;                    
          /* 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 = ighmm_mprintf (NULL, 0,
                  "no matching in_a for out_a(=a[%d][%d]) found!\n", i, j_id);
          GHMM_LOG(LCONVERTED, 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)
        ighmm_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 */


    /* if fix, continue to next state */
    if (smo->s[i].fix){
      continue;
    }
    /* C, Mue und U */


    c_denom_pos = 1;
    if (r->c_denom[i] <= DBL_MIN) {     /* < EPS_PREC ? */
#if MCI
      ighmm_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;
    fix_w = 1.0;
    unfix_w = 0.0;
    fix_flag = 0;
    for (m = 0; m < smo->s[i].M; m++) {

      /* if fixed continue to next component */
      if (smo->s[i].mixture_fix[m] == 1) {
        /*printf("state %d, component %d is fixed !\n",i,m);*/
        fix_w = fix_w - smo->s[i].c[m];
        fix_flag = 1;           /* we have to normalize weights -> fix flag is set to one */
        continue;
      }

      /* TEST: denom. < numerator */
      if ((r->c_denom[i] - r->c_num[i][m]) < 0.0) {     /* < -EPS_PREC ? */
#if MCI
        ighmm_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;

      unfix_w = unfix_w + smo->s[i].c[m];

      if (fabs (r->mue_u_denom[i][m]) <= DBL_MIN)       /* < EPS_PREC ? */
#if MCI
        ighmm_mes (MESCONTR, "mue[%d][%d]: denominator == 0.0!\n", i, m);
#else
        ;
#endif
      else {
        /* set mue_im */
        smo->s[i].mue[m] = r->mue_num[i][m] / r->mue_u_denom[i][m];
      }

      /* TEST: u_denom == 0.0 ? */