discrime.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/discrime.c
*       Authors:  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: 1754 $
*                       from $Date: 2006-11-03 17:29:44 +0100 (Fri, 03 Nov 2006) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

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

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

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

#define LAMBDA 0.14
#define TRIM(o, n) ((1-LAMBDA)*(o) + LAMBDA*(n))

#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
#else
#define logl(A) log(A)
#define expl(A) exp(A)
#endif

/* forward declaration */
static int discrime_galloc (ghmm_dmodel ** mo, ghmm_dseq ** sqs, int noC,
                     double ******matrix_b, double *****matrix_a,
                     double *****matrix_pi, long double ***omega,
                     long double ****omegati, double ****log_p);

static void discrime_gfree (ghmm_dmodel ** mo, ghmm_dseq ** sqs, int noC,
                     double *****matrix_b, double ****matrix_a,
                     double ****matrix_pi, long double **omega,
                     long double ***omegati, double ***log_p);

static void discrime_trim_gradient (double *new, int length);

static double discrime_lambda = 0.0;
static double discrime_alpha = 1.0;


/*----------------------------------------------------------------------------*/
/** allocates memory for m and n matrices: */
static int discrime_galloc (ghmm_dmodel ** mo, ghmm_dseq ** sqs, int noC,
                     double ******matrix_b, double *****matrix_a,
                     double *****matrix_pi, long double ***omega,
                     long double ****omegati, double ****log_p)
{
#define CUR_PROC "discrime_galloc"

  int i, k, l, m;

  /* first allocate memory for matrix_b: */
  ARRAY_CALLOC (*matrix_b, noC);
  for (k = 0; k < noC; k++) {
    ARRAY_CALLOC ((*matrix_b)[k], sqs[k]->seq_number);
    for (l = 0; l < sqs[k]->seq_number; l++) {
      ARRAY_CALLOC ((*matrix_b)[k][l], noC);
      for (m = 0; m < noC; m++) {
        ARRAY_CALLOC ((*matrix_b)[k][l][m], mo[m]->N);
        for (i = 0; i < mo[m]->N; i++)
          ARRAY_CALLOC ((*matrix_b)[k][l][m][i], ghmm_ipow (mo[m], mo[m]->M, mo[m]->order[i] + 1));
      }
    }
  }

  /* matrix_a(k,l,i,j) = matrix_a[k][l][i*mo->N + j] */
  ARRAY_CALLOC (*matrix_a, noC);
  for (k = 0; k < noC; k++) {
    ARRAY_CALLOC ((*matrix_a)[k], sqs[k]->seq_number);
    for (l = 0; l < sqs[k]->seq_number; l++) {
      ARRAY_CALLOC ((*matrix_a)[k][l], noC);
      for (m = 0; m < noC; m++)
        ARRAY_CALLOC ((*matrix_a)[k][l][m], mo[m]->N * mo[m]->N);
    }
  }

  /* allocate memory for matrix_pi */
  ARRAY_CALLOC (*matrix_pi, noC);
  for (k = 0; k < noC; k++) {
    ARRAY_CALLOC ((*matrix_pi)[k], sqs[k]->seq_number);
    for (l = 0; l < sqs[k]->seq_number; l++) {
      ARRAY_CALLOC ((*matrix_pi)[k][l], noC);
      for (m = 0; m < noC; m++)
        ARRAY_CALLOC ((*matrix_pi)[k][l][m], mo[m]->N);
    }
  }

  /* allocate memory for matrices of likelihoods 
     log_p[k][l][m] =
     log_prob of l-th sequence of k-th class under the m-th ghmm_dmodel */
  ARRAY_CALLOC (*log_p, noC);
  for (k = 0; k < noC; k++) {
    ARRAY_CALLOC ((*log_p)[k], sqs[k]->seq_number);
    for (l = 0; l < sqs[k]->seq_number; l++)
      ARRAY_CALLOC ((*log_p)[k][l], noC);
  }

  /* allocate memory for outer derivatives */
  ARRAY_CALLOC (*omega, noC);
  for (k = 0; k < noC; k++)
    ARRAY_CALLOC ((*omega)[k], sqs[k]->seq_number);

  /* allocate memory for omega tilde. NB: size(omega)*noC == size(omegati) */
  ARRAY_CALLOC (*omegati, noC);
  for (k = 0; k < noC; k++) {
    ARRAY_CALLOC ((*omegati)[k], sqs[k]->seq_number);
    for (l = 0; l < sqs[k]->seq_number; l++)
      ARRAY_CALLOC ((*omegati)[k][l], noC);
  }

  return 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  discrime_gfree (mo, sqs, noC, *matrix_b, *matrix_a, *matrix_pi,
                  *omega, *omegati, *log_p);
  return -1;
#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
/** frees memory for expectation matrices for all classes & training sequences*/
static void discrime_gfree (ghmm_dmodel ** mo, ghmm_dseq ** sqs, int noC,
                     double *****matrix_b, double ****matrix_a,
                     double ****matrix_pi, long double **omega,
                     long double ***omegati, double ***log_p)
{
#define CUR_PROC "discrime_gfree"

  int i, k, l, m;

  for (k = 0; k < noC; k++) {
    for (l = 0; l < sqs[k]->seq_number; l++) {
      for (m = 0; m < noC; m++) {
        for (i = 0; i < mo[m]->N; i++)
          m_free (matrix_b[k][l][m][i]);
        m_free (matrix_b[k][l][m]);
      }
      m_free (matrix_b[k][l]);
    }
    m_free (matrix_b[k]);
  }
  m_free (matrix_b);

  for (k = 0; k < noC; k++) {
    for (l = 0; l < sqs[k]->seq_number; l++) {
      for (m = 0; m < noC; m++)
        m_free (matrix_a[k][l][m]);
      m_free (matrix_a[k][l]);
    }
    m_free (matrix_a[k]);
  }
  m_free (matrix_a);

  for (k = 0; k < noC; k++) {
    for (l = 0; l < sqs[k]->seq_number; l++) {
      for (m = 0; m < noC; m++)
        m_free (matrix_pi[k][l][m]);
      m_free (matrix_pi[k][l]);
    }
    m_free (matrix_pi[k]);
  }
  m_free (matrix_pi);

  for (k = 0; k < noC; k++) {
    for (l = 0; l < sqs[k]->seq_number; l++)
      m_free (log_p[k][l]);
    m_free (log_p[k]);
  }
  m_free (log_p);

  for (k = 0; k < noC; k++)
    m_free (omega[k]);
  m_free (omega);

  for (k = 0; k < noC; k++) {
    for (l = 0; l < sqs[k]->seq_number; l++)
      m_free (omegati[k][l]);
    m_free (omegati[k]);
  }
  m_free (omegati);

  return;
#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
static int discrime_calculate_omega (ghmm_dmodel ** mo, ghmm_dseq ** sqs, int noC,
                              long double **omega, long double ***omegati,
                              double ***log_p)
{
#define CUR_PROC "discrime_calculate_omega"
  int k, l, m;
  int seq_no, argmax = 0;
  double sum, max;
  double exponent;
  long double sigmoid;

  /* compute outer derivatives */
  /* iterate over all classes */
  for (k = 0; k < noC; k++) {
    seq_no = sqs[k]->seq_number;
    /*iterate over all training sequences */
    for (l = 0; l < seq_no; l++) {

      max = 1.0;
      /* calculate log(sum[prob]) for logarithmic probabilities
         first determine the maximum */
      for (m = 0; m < noC; m++) {
        if (m != k
            && (1.0 == max || max < (log_p[k][l][m] + log (mo[m]->prior)))) {
          max = log_p[k][l][m] + log (mo[m]->prior);
          argmax = m;
        }
      }

      /* sum */
      sum = 1.0;
      for (m = 0; m < noC; m++) {
        if (m != k && m != argmax)
          sum += exp (log_p[k][l][m] + log (mo[m]->prior) - max);
      }
      sum = log (sum) + max;

      exponent = log_p[k][l][k] + log (mo[k]->prior) - sum;
      if (exponent < logl (LDBL_MIN)) {
        printf ("exponent was too large (%g) cut it down!\n", exponent);
        exponent = (double) logl (LDBL_MIN);
      }

      sigmoid = 1.0 / (1.0 + expl ((-discrime_alpha) * exponent));

      omega[k][l] = discrime_alpha * sigmoid * (1.0 - sigmoid);
/*       printf("omega[%d][%d] = %Lg\n",k ,l, omega[k][l]); */
      /* omegati is not useful for m==k equation (2.9) */
      for (m = 0; m < noC; m++) {
        if (m == k)
          omegati[k][l][m] = 0;
        else
          omegati[k][l][m] =
            omega[k][l] * expl (log_p[k][l][m] -
                                sum) * (long double) mo[m]->prior;
/* 	printf("%Lg, ", expl(log_p[k][l][m] - sum) * (long double)mo[m]->prior); */
/* 	printf("omegati[%d][%d][%d] = %Lg\n",k ,l, m, omegati[k][l][m]); */
      }
/*       printf("\n"); */
    }
  }
  return 0;
#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
static int discrime_precompute (ghmm_dmodel ** mo, ghmm_dseq ** sqs, int noC,
                         double *****expect_b, double ****expect_a,
                         double ****expect_pi, double ***log_p)
{
#define CUR_PROC "discrime_precompute"

  int k, l, m;
  int seq_len, success = 1;

  /* forward, backward variables and scaler */
  double **alpha, **beta, *scale;

  ghmm_dseq *sq;

  /* iterate over all classes */
  for (k = 0; k < noC; k++) {

    sq = sqs[k];

    /*iterate over all training sequences */
    for (l = 0; l < sq->seq_number; l++) {

      seq_len = sq->seq_len[l];

      /* iterate over all classes */
      for (m = 0; m < noC; m++) {

        if (-1 ==
            ighmm_reestimate_alloc_matvek (&alpha, &beta, &scale, seq_len,
                                     mo[m]->N))
          return -1;

        /* calculate forward and backward variables without labels: */
        if (-1 == ghmm_dmodel_forward (mo[m], sq->seq[l], seq_len, alpha, scale,
                                &(log_p[k][l][m]))) {
          success = 0;
          printf ("forward\n");
          goto FREE;
        }
        if (-1 == ghmm_dmodel_backward (mo[m], sq->seq[l], seq_len, beta, scale)) {
          success = 0;
          printf ("backward\n");
          goto FREE;
        }

        /* compute expectation matrices
           expect_*[k][l][m] holds the expected number how ofter a particular
           parameter of model m is used by l-th training sequence of class k */
        ghmm_dmodel_label_gradient_expectations (mo[m], alpha, beta, scale,
                                         sq->seq[l], seq_len,
                                         expect_b[k][l][m], expect_a[k][l][m],
                                         expect_pi[k][l][m]);

      FREE:
        ighmm_reestimate_free_matvek (alpha, beta, scale, seq_len);
        if (!success)
          return -1;
      }
    }
  }
  return 0;
#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
double ghmm_dmodel_label_discrim_perf(ghmm_dmodel** mo, ghmm_dseq** sqs, int noC)
{
#define CUR_PROC "ghmm_dmodel_label_discrim_perf"

  int k, l, m, temp;
  int argmax = 0;
  double sum, max;
  double *logp;
  double exponent;
  long double sigmoid;
  double performance = 0.0;

  ghmm_dseq *sq;

  ARRAY_CALLOC (logp, noC);

  /* iterate over all classes */
  for (k = 0; k < noC; k++) {
    sq = sqs[k];
    /*iterate over all training sequences */
    for (l = 0; l < sq->seq_number; l++) {

      /* iterate over all classes */
      for (m = 0; m < noC; m++) {
        temp = ghmm_dmodel_logp (mo[m], sq->seq[l], sq->seq_len[l], &(logp[m]));
        if (0 != temp)
          printf ("ghmm_dmodel_logp error in sequence[%d][%d] under model %d (%g)\n",
                  k, l, m, logp[m]);
        /*printf("ghmm_dmodel_logp sequence[%d][%d] under model %d (%g)\n", k, l, m, logp[m]);*/
      }

      max = 1.0;
      for (m = 0; m < noC; m++) {
        if (m != k && (1.0 == max || max < (logp[m] + log (mo[m]->prior)))) {
          max = logp[m] + log (mo[m]->prior);
          argmax = m;
        }
      }

      /* sum */
      sum = 1.0;
      for (m = 0; m < noC; m++) {
        if (m != k && m != argmax)
          sum += exp (logp[m] + log (mo[m]->prior) - max);
      }
      sum = log (sum) + max;

      exponent = logp[k] + log (mo[k]->prior) - sum;

      if (exponent < logl (LDBL_MIN)) {
        printf ("exponent was too large (%g) cut it down!\n", exponent);
        exponent = (double) logl (LDBL_MIN);
      }

      sigmoid = 1.0 / (1.0 + expl ((-discrime_alpha) * exponent));

      performance += (double) sigmoid;
    }
  }
  m_free (logp);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  return performance;
#undef CUR_PROC
}


/*----------------------------------------------------------------------------*/
static void discrime_print_statistics(ghmm_dmodel** mo, ghmm_dseq** sqs, int noC,
                                int* falseP, int* falseN)