hmm.cpp

一个用于 词性标注的 HMM程序。 包含 训练和测试功能。

//   file : hmm.cc
// version: 1.03 [August 21, 1995]
/* Copyright (C) 1994 Richard Myers and James Whitson
   Permission is granted to any individual or institution to use, copy, or
   redistribute this software so long as all of the original files are 
   included unmodified, that it is not sold for profit, and that this 
   copyright notice is retained.
*/

/*  *************************************************************
    Updated by Emdad Khan (Computer & Information Science dept,
    U.C Santa Cruz and National Semiconductor, Santa Clara, CA).
    Date 5/12/95. Brief description of the changes follows:

    (1)
    The original model written by Richard & James seemed to works for equal
    length string only. For unequal length strings (i.e vq coded words),
    i.e when different people say the same words, the program does not 
    converge (basically oscillates). The authors used "instanteneous"
    updates for each strings for the unequal case. Such scheme nullifies
    the learning a string by another etc. The commonly used approach
    is the "batch" update: a good reference is Rabiner's book or paper
    (as alos referenced by the authors). 
     
    So, I updated the code (training routine and other associated routines)
    to use batch update for unequal length strings. More detailed descrip-
    tion of the changes are given in corresponding routines i.e
    batch_train, beta_I, compute_gamma, alpha_F and reestimate routines.

    (2) 
    Some of equations related to training (namely, gamma, beta, alpha)
    showed discrepencies compared to the equations given in HMM literatures.
    Those equations are corrected and details of the changes are shown in
    corresponding routines. These routines are basically same as the 
    routines mentioned above.

    (3) 
    Some of the parameters were initialised both in the main routines as 
    well as in the function declarations in the .h file. Some of the 
    compilers does not like this. So, initialization of parameters are
    changed so that it is present only in one place. 
    
    *******************************************************************
*/ 
 
/*  *************************************************************
    Updated by Arthur Stephens ([e-mail])
    Date 8/30/95.

    Purified code, plugged memory leaks, corrected some out of bounds
    array conditions and got it to compile on Solaris. (Thanks Art --rm)
*/


#include <iostream.h>
#include <fstream.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "hmm.h"
//#include "random.h"

/************************* HMM DECLARATION ****************************/

// create hmm, initialize model parameters from file.
//
HMM::HMM(char* filename)
{
  // set max_states, max_symbols and model parameters from file
  load_model(filename);
  
  cout << "\nUsing predefined model...\n";
  // show_probs();
}

// create hmm, initialize model parameters with random values.
//
HMM::HMM(int symbols, int states, int new_seed)
{
  if (new_seed != -1)
    set_seed(new_seed);

  // set globals
  max_symbols = symbols;
  max_states = states;

  // allocate space for model parameters
  alloc_model_matrices();

  // rnd initialize the trans and emit probs
  rnd_init_hmm_probs();

  cout << "\nUsing random model...\n";
  // show_probs();
}

// delete hmm, delete arrays and matrixes.
//
HMM::~HMM()
{
  // delete state and some numerator matrixes
  int i;
  for (i=0; i < max_states; i++) {
   
    delete States[i];
  }
  
  delete [] States;

  // delete strings matrixes
  for (i=0; i<num_strings; i++) {
    delete [] strings[i];
  }
  delete [] strings;
  delete [] string_len;

  // delete alpha,beta and gamma matrixes
  for (i=0; i < trellis_width; i++) {
    delete [] alpha[i];
    delete [] beta[i];
    delete [] gamma[i];
    delete [] sum_gamma[i];
  }
  delete [] alpha;
  delete [] beta;
  delete [] gamma;
  delete [] sum_gamma;

  // delete numerator, denominator and scaling_factors arrays
 
  delete [] scaling_factors;
}

// train hmm from sequence file until the sum of the probability
// changes is less than min_delta_psum.
// *** Changes by Emdad Khan: the code for the unequal length strings
// is changed to batch mode. It does not work in nonbatch mode.
// The equations given in Rabiner and elsewhere also indicates batch 
// mode training. Yi[] and yi[][] are used to add the batch mode
// update. ***

void HMM::batch_train(char *filename, double min_delta_psum) 
{
	int symbol_count,str;

	// local arrays and integers used in the case with strings of unequal lengths

	int i,j,k;

	// dynamic allocation of arrays (ANSI)


	// load string_matrix and set globals
	symbol_count = load_string_matrix(filename);

	// define trellis dimensions
	trellis_width = symbol_count+1;

	// allocate space for all the training arrays
	alloc_training_matrices();

	// loop until local min in probs
	double prob_sum=999999.0; // magic number,ick
	double alpha_tot,last_alpha_tot;
	double temp1,temp2,temp3;

	for (str=0,alpha_tot=0.0; str<num_strings; str++)//根据模型计算样本的的初始概率
		alpha_tot+=test(strings[str],string_len[str]);
	cout << "\nlog(alpha)\tprob_sum\n";
	cout << "----------\t--------\n";
	cout << alpha_tot <<"\n";

	// the following code does batch update for each case:
	// diff_len == true and != true 
	while (prob_sum > min_delta_psum) 
	{
	 
		
		last_alpha_tot = alpha_tot;

		

        //将计数清空 准备训练
		for (i=0 ; i < max_states; i++)  
		{
							
			for (k=0; k < max_symbols; k++) 
				sum_b[i][k] =0;

			for (j=0; j<max_states; j++)
				sum_aij[i][j]=0;
			
			for(k=0; k<trellis_width; k++ ) 
				sum_gamma[k][i]=0;
	
		}


         //针对每一个实例进行计算，最后统一进行参数估计
		for (str=0, alpha_tot=0.0; str<num_strings; str++) 
		{
			//计算alpha,beta,gamma，epsum的值
			set_ab_counts(strings[str],string_len[str]);

			//累计gamma
			for(k=0; k<trellis_width; k++)
				for(i=0; i<max_states; i++)
					sum_gamma[k][i]+=gamma[k][i];

			//累计aij
			for (i=0;i< max_states; i++) 
				for(j=0; j<max_states; j++)
				{
					
					temp1 =0;
					temp2 =0;
					for(k=0; k<trellis_width; k++)
					{
						temp1 += epsum[k][i][j];
						temp2 +=gamma[k][i]; 
					}
					sum_aij[i][j] +=temp1/temp2;
				}
            //累计bj
			for (k=0;k< string_len[str]; k++) 
				for(i=0; i<max_states; i++)
				{
					temp3=0;
					for (int l=0; l<trellis_width ;l++)
						temp3 +=gamma[l][i];

					sum_b[i][strings[str][k]] +=gamma[k][i]/temp3;
					
				
				}
		
		
		}
		// end of top for loop


			// 重新估计模型的参数
		prob_sum = reestimate();

		for (str=0, alpha_tot=0.0; str<num_strings; str++) 
			alpha_tot+=test(strings[str],string_len[str]);     
		
	
		cout << alpha_tot << "\t";
		cout << prob_sum << "\n";
		if (alpha_tot < last_alpha_tot) 
		{
			cerr << "\tWARNING: alpha values reversed!\n"; // should not occur!
			// break;
		}
	}
	
	cout << "\nDone training model " << filename << "...\n";
	// show_probs();


}  

// test hmm from sequence file.
//
double HMM::batch_test(char *filename)
{
  int symbol_count;
  FILE* ofp; // output file to write alpha values

  ofp = fopen("alphaout", "w");
  // load string_matrix and set globals
  symbol_count = load_string_matrix(filename);

  // define trellis dimensions
  trellis_width = symbol_count+1;

  // create space for testing strs
  alloc_testing_matrices();

  // loop through all test strings
  cout << "\ntest on strings...\n";
  double total=0.0,val;
  for (int i=0; i<num_strings; i++) {
    if (diff_len==TRUE) 
      val = test(strings[i],string_len[i]);
    else
      val = test(strings[i]);
    cout << "alpha for string[" << i << "] = " << val << " ";
    if (diff_len==TRUE)
      cout << " len= " << string_len[i];
    cout << "\n";
    fprintf(ofp, " %f", val); // writing alpha values to file
    if (i == (num_strings - 1))
    fprintf(ofp, "\n");
    total+=val;
  }
  fclose(ofp);
  cout << "------------ \nbatch alpha = "<<total<<"\n";

  cout << "\nDone testing model on" << filename << "...\n"; 
return(total);
}

// output the model parameters.
//
void HMM::show_probs()
{
  cout.precision(6);
  cout.setf(ios::fixed,ios::floatfield);

  cout << "\n\nstates: " << max_states;
  cout << "\nsymbols: " << max_symbols << "\n";
  for (int s=0; s < max_states; s++) 
  {
    cout << States[s]->start << " ";
	cout << "\n";
    int k;

	for (k=0; k < max_states; k++)
      cout << States[s]->next_trans[k] << " ";

	cout << "\n";

    for (k=0; k < max_symbols; k++)
      cout << States[s]->recur_out[k] << " ";
    cout << "\n";
   
  }
  cout << "\n";
}

// output the model parameters to a file.
//输出模型
void HMM::dump_model(char* filename)
{
  ofstream model(filename);
  model.precision(6);
  model.setf(ios::fixed,ios::floatfield);

  if (!model) {
    cerr << "ERROR: Couldn't create file for model dump.";
    exit(-1);
  }

  model << "states: " << max_states << "\n";
  model << "symbols: " << max_symbols << "\n";
  for (int s=0; s < max_states; s++) 
  {
	model << States[s]->start << "\t";
	model << "\n";
    int k;
	for(k=0; k<max_states; k++)
      model << States[s]->next_trans[k] << "\t";

	 model << "\n";
    
    for (k=0; k < max_symbols; k++)
      model << States[s]->recur_out[k] << "\t";
   
    
    model << "\n\n";
  }
  model << "\n";
  model.close();
  cout << "\nDumped model to file ==> " << filename << "\n";
}



/************************* SPACE ALLOCATION ****************************/

// allocate space for state probs.
//
void HMM::alloc_model_matrices() {

  // check params
  if (max_states <= 0 || max_symbols <= 0) {
    cerr <<"ERROR: alloc_model_matrices(), must define global sizes\n";
    exit(-1);
  }

  States = new STATE* [max_states];

  for (int i=0; i < max_states; i++) {
    States[i] = new STATE(max_states,max_symbols);
  }
}

// allocate space for training trellises.
//训练的时候需要的参数
void HMM::alloc_training_matrices()
{
  // check parms
  if (max_states <= 0 || max_symbols <= 0 || trellis_width <=0 ) {
    cerr <<"ERROR: alloc_training_matrices(), must define global sizes\n";
    exit(-1);
  }
  

  alpha = new double* [trellis_width];
  beta = new double*  [trellis_width];
  gamma = new double* [trellis_width];
  
  epsum =new double **[trellis_width];


  ///////////
  sum_aij= new double *[max_states];
  sum_b= new double *[max_states];
  sum_gamma = new double* [trellis_width];

  /////////////////////

  int i; 
  int j;
  for (i=0; i < trellis_width; i++) 
  {
    alpha[i] = new double [max_states];
    beta[i] = new double [max_states];
    gamma[i] = new double [max_states];
	sum_gamma[i] = new double [max_states];
    epsum[i] = new double* [max_states];
	for(j=0; j < max_states; j++)
		epsum[i][j]= new double [max_states];
  }

  for(i =0; i<max_states; i++)
  {
	  sum_aij[i]= new double [max_states];
	  sum_b[i] =new double [max_symbols];
  }

  // scaling factors 
  // fixed [trellis_width-1] (rem) 8/17/95
  scaling_factors = new double[trellis_width];

 
}

// allocate space for testing trellises.
//
void HMM::alloc_testing_matrices()
{
  // check params
  if (max_states <= 0 || max_symbols <= 0 || trellis_width <= 0 ) {
    cerr <<"ERROR: alloc_testing_matrices(), must define global sizes\n";
    exit(-1);
  }

  // space for calculation tables
  // two dimensional array decls
  alpha = new double* [trellis_width];
 
  for (int i=0; i < trellis_width; i++) {
    alpha[i] = new double [max_states];
  }

  // scaling factors
  scaling_factors = new double[trellis_width];
}

/************************* MODEL INITIALIZATION ****************************/

// randomly initialize the model parameters
//随机初始化发射概率
void HMM::rnd_init_hmm_probs()
{
  double prob_trans, *prob_syms,*prob_states;

  // init rnd gen
  //srandom(seed);
  srand(seed);

  // dynamic allocation of arrays (ANSI)
  prob_syms = new double [max_symbols];
  prob_states = new double [max_states];

  for (int i=0; i < max_states; i++)
    {
      //prob_trans = double(random() % 100) / 100.0;    
      prob_trans = double(rand() % 100) / 100.0;    
      rnd_sym_probs(prob_syms);
	  rnd_state_probs(prob_states);
      
      int j;

      for (j=0; j < max_symbols; j++)
      	States[i]->set_recur_out(j,prob_syms[j]);
      rnd_sym_probs(prob_syms);
      

      for (j=0; j < max_states; j++)
		States[i]->set_next_trans(j,prob_states[j]);
    }
 

  delete [] prob_syms;
  delete [] prob_states;
}

// set output symbol probabilities.
//
void HMM::rnd_sym_probs(double prob_syms[])
{
  // determine the trans and out probs
  double prob_sum = 0.0;
  int s;
  for (s=0; s < max_symbols; s++) {
    //prob_syms[s] = double(random() % 100) / 100.0;
    prob_syms[s] = double(rand() % 100) / 100.0;
    prob_sum += prob_syms[s];
  }
  for (s=0; s < max_symbols; s++)  // normalize so add to one
    prob_syms[s] = prob_syms[s] / prob_sum;
}

//随机初始化状态转移概率
void HMM::rnd_state_probs(double prob_states[])
{
  // determine the trans and out probs
  double prob_sum = 0.0;
  int s;
  for (s=0; s < max_states; s++) {
    
    prob_states[s] = double(rand() % 100) / 100.0;
    prob_sum += prob_states[s];