gradescent.c

来自「General Hidden Markov Model Library 一个通用」· C语言 代码 · 共 556 行 · 第 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/gradescent.c
*       Authors:  Janne Grunau, Alexander Riemer
*
*       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 <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "ghmm.h"
#include "mes.h"
#include "mprintf.h"
#include "matrix.h"
#include "model.h"
#include "foba.h"
#include "gradescent.h"
#include "ghmm_internals.h"


static double compute_performance (ghmm_dmodel * mo, ghmm_dseq * sq);

/*----------------------------------------------------------------------------*/
static void gradient_descent_gfree (double **matrix_b, double *matrix_a,
                             double *matrix_pi, int N)
{
#define CUR_PROC "gradient_descent_gfree"

  int i;

  if (matrix_b)
    for (i = 0; i < N; i++) {
      m_free (matrix_b[i]);
    }
  m_free (matrix_b);

  m_free (matrix_a);
  m_free (matrix_pi);

#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
/** allocates memory for m and n matrices: */
static int gradient_descent_galloc (double ***matrix_b, double **matrix_a,
                             double **matrix_pi, ghmm_dmodel * mo)
{
#define CUR_PROC "gradient_descent_galloc"

  int i;

  /* first allocate memory for matrix_b */
  ARRAY_MALLOC (*matrix_b, mo->N);
  for (i = 0; i < mo->N; i++)
    ARRAY_CALLOC ((*matrix_b)[i], ghmm_ipow (mo, mo->M, mo->order[i] + 1));

  /* matrix_a(i,j) = matrix_a[i*mo->N+j] */
  ARRAY_CALLOC (*matrix_a, mo->N * mo->N);

  /* allocate memory for matrix_pi */
  ARRAY_CALLOC (*matrix_pi, mo->N);

  return 0;

STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  gradient_descent_gfree (*matrix_b, *matrix_a, *matrix_pi, mo->N);
  return -1;

#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
/**
   computes matrices of n and m variables (expected values for how often a
   certain parameter from A or B is used)
   computes Baum-Welch variables implicit 
   @return                 0/-1 success/error
   @param mo:              pointer to a ghmm_dmodel
   @param alpha:           matrix of forward variables
   @param backward:        matrix of backward variables
   @param scale:           scaling vector from forward-backward-algorithm
   @param seq:             sequence in internal representation
   @param seq_len:         length of sequence
   @param matrix_b:        matrix for parameters from B (n_b or m_b)
   @param matrix_a:        matrix for parameters from A (n_a or m_a)
   @param vec_pi:          vector for parameters in PI (n_pi or m_pi)
 */
int ghmm_dmodel_label_gradient_expectations (ghmm_dmodel* mo, double **alpha,
                                     double **beta, double *scale, int *seq,
                                     int seq_len, double **matrix_b,
                                     double *matrix_a, double *vec_pi)
{
#define CUR_PROC "ghmm_dmodel_label_gradient_expectations"

  int h, i, j, t;

  int size, j_id, e_index;

  double gamma, xi;
  double foba_sum;

  /* initialise matrices with zeros */
  for (i = 0; i < mo->N; i++) {
    for (j = 0; j < mo->N; j++)
      matrix_a[i * mo->N + j] = 0;
    size = ghmm_ipow (mo, mo->M, mo->order[i] + 1);
    for (h = 0; h < size; h++)
      matrix_b[i][h] = 0;
  }

  for (t = 0; t < seq_len; t++) {

    /* sum products of forward and backward variables over all states: */
    foba_sum = 0.0;
    for (i = 0; i < mo->N; i++)
      foba_sum += alpha[t][i] * beta[t][i];
    if (GHMM_EPS_PREC > fabs (foba_sum)) {
      printf
        ("gradescent_compute_expect: foba_sum (%g) smaller as EPS_PREC (%g). t = %d.\n",
         foba_sum, GHMM_EPS_PREC, t);
      return -1;
    }

    for (i = 0; i < mo->N; i++) {

      /* compute gamma implicit */
      gamma = alpha[t][i] * beta[t][i] / foba_sum;

      /* n_pi is easiest: n_pi(i) = gamma(0,i) */
      if (0 == t)
        vec_pi[i] = gamma;

      /* n_b(i,c) = sum[t, hist(t)=c | gamma(t,i)] / sum[t | gamma(t,i)] */
      e_index = get_emission_index (mo, i, seq[t], t);
      if (-1 != e_index)
        matrix_b[i][e_index] += gamma;
    }

    /* updating history, xi needs the right e_index for the next state */
    update_emission_history (mo, seq[t]);

    for (i = 0; i < mo->N; i++) {
      /* n_a(i,j) = sum[t=0..T-2 | xi(t,i,j)] / sum[t=0..T-2 | gamma(t,i)] */
      /* compute xi only till the state before the last */
      for (j = 0; (j < mo->s[i].out_states) && (t < seq_len - 1); j++) {
        j_id = mo->s[i].out_id[j];

        /* compute xi implicit */
	xi = 0;
        e_index = get_emission_index (mo, j_id, seq[t + 1], t + 1);
        if (e_index != -1)
          xi = alpha[t][i] * beta[t + 1][j_id] * mo->s[i].out_a[j]
            * mo->s[j_id].b[e_index] / (scale[t + 1] * foba_sum);

        matrix_a[i * mo->N + j_id] += xi;
      }
    }
  }

  return 0;
#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
static double compute_performance (ghmm_dmodel * mo, ghmm_dseq * sq)
{
#define CUR_PROC "compute_performance"

  int k, seq_len, success = 1;
  /* log P[O | lambda, labeling] as computed by the forward algorithm */
  double log_p;
  /* sum over log P (calculated by forward_label) for all sequences
     used to compute the performance of the training */
  double log_p_sum = 0.0;

  /* loop over all sequences */
  for (k = 0; k < sq->seq_number && success; k++) {
    success = 0;
    seq_len = sq->seq_len[k];

    if (-1 !=
        ghmm_dmodel_label_logp (mo, sq->seq[k], sq->state_labels[k], seq_len,
                         &log_p)) {
      success = 1;
      log_p_sum += log_p;

      if (-1 != ghmm_dmodel_logp (mo, sq->seq[k], seq_len, &log_p))
        log_p_sum -= log_p;
      else {
        printf ("ghmm_dmodel_logp error in sequence %d, length: %d\n", k,
                seq_len);
        success = 0;
      }
    }
    else
      printf ("ghmm_dl_logp error in sequence %d, length: %d\n", k,
              seq_len);
  }

  /* return log_p_sum in success or +1.0 a probality of 0.0 on error */
  if (success)
    return log_p_sum;
  else
    return 1.0;
#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
/**
   Trains the model with a set of annotated sequences using gradient descent.
   Model must not have silent states. (one iteration)
   @return            0/-1 success/error
   @param mo:         pointer to a ghmm_dmodel
   @param sq:         struct of annotated sequences
   @param eta:        training parameter for gradient descent
 */
static int gradient_descent_onestep (ghmm_dmodel * mo, ghmm_dseq * sq, double eta)
{
#define CUR_PROC "gradient_descent_onestep"

  int g, h, i, j, k;
  int seq_len, j_id, hist, size;

  /* expected usage of model parameters (A and B entries) */
  double *m_pi, *n_pi, *m_a, *n_a;
  double **m_b, **n_b;

  /* forward & backward variables w/ their scaling vector */
  double **alpha, **beta, *scale;

  /* variables to store values associated with the algorithm */
  double pi_sum, a_row_sum, b_block_sum, gradient;

  /* log P[O | lambda, labeling] as computed by the forward algorithm */
  double log_p;

  /* allocate memory for the parameters used for reestimation */
  if (-1 == gradient_descent_galloc (&m_b, &m_a, &m_pi, mo))
    return -1;
  if (-1 == gradient_descent_galloc (&n_b, &n_a, &n_pi, mo))