taginference.cpp

文本倾向性分析程序

#include "TagInference.h"

using namespace std;

CTagInference::CTagInference()
{
}

CTagInference::~CTagInference()
{
}

void CTagInference::Load()
{
	// Read setting
	FILE* fileptr;
    fileptr = fopen("settings.txt", "r");
	fscanf(fileptr, "class number %d\n", &nclass);
    fscanf(fileptr, "var max iter %d\n", &emmax);
    fscanf(fileptr, "var convergence %lf\n", &epsilon);
    fclose(fileptr);

	// Read words
	ifstream inWord(WORD_LIST);
	string strBuffer;
	while( getline(inWord, strBuffer) )
	{
		m_vec_word.push_back(strBuffer);
	}
	nlex = m_vec_word.size();
	inWord.close();

	// Read tags
	/*ifstream inTag(TAG_LIST);
	while( getline(inTag, strBuffer) )
	{
		m_vec_tag.push_back(strBuffer);
	}
	ntag = m_vec_tag.size();
	inTag.close();*/

	// Read alpha
	ifstream inAlpha(FINAL_ALPHA);
	alpha = (double *)calloc(nclass, sizeof(double));
	getline(inAlpha, strBuffer);
	strtokenizer strToken(strBuffer, "\t");
	int i;
	for( i = 0; i < nclass; i++)
	{
		alpha[i] = atof(strToken.token(i).c_str());
	}
	inAlpha.close();

	// Read beta
	ifstream inBeta(FINAL_BETA);
	beta = dmatrix(nlex, nclass);
	if (beta == NULL)
		return;
	for( i = 0; i < nlex; i++)
	{
		getline(inBeta, strBuffer);
		strtokenizer strToken(strBuffer, "\t");
		for(int j = 0; j < nclass; j++)
		{
			beta[i][j] = atof(strToken.token(j).c_str());
		}
	}
	inBeta.close();

	// Read eta
/*	ifstream inEta(FINAL_ETA);
	eta = dmatrix(ntag, nclass);
	if (eta == NULL)
		return;
	for( i = 0; i < ntag; i++)
	{
		getline(inEta, strBuffer);
		strtokenizer strToken(strBuffer, "\t");
		for(int j = 0; j < nclass; j++)
		{
			eta[i][j] = atof(strToken.token(j).c_str());
		}
	}
	inEta.close();*/

	gamma = (double *)calloc(nclass, sizeof(double));
	ap = (double *)calloc(nclass, sizeof(double));
	nt = (double *)calloc(nclass, sizeof(double));
}


void CTagInference::Evaluate()
{
	Load();

	ifstream inData(TEST_DATA);
	string strBuffer;
	double accuracy = 0.0;
	int count = 0;
	while( getline(inData, strBuffer) )
	{
		count++;
		vector<rank_t> vecModel;
		vector<rank_t> vecTrue;
		vector<rank_t> vecWord;

		// Read current document
		strtokenizer strText(strBuffer," ");
		document cur_doc;
		cur_doc.len = strText.count_tokens();
		cur_doc.size = 0;
		cur_doc.id  = (int *)calloc(cur_doc.len, sizeof(int));
		cur_doc.cnt = (double *)calloc(cur_doc.len, sizeof(double));
		for(int i=0; i<strText.count_tokens(); i++)
		{
			strtokenizer strTerm(strText.token(i),":");
			cur_doc.id[i] = atoi(strTerm.token(0).c_str());
			cur_doc.cnt[i] = atoi(strTerm.token(1).c_str());
			cur_doc.size += atoi(strTerm.token(1).c_str());
			rank_t temp_rank;
			temp_rank.r_str = cur_doc.id[i];
			temp_rank.r_val = cur_doc.cnt[i];
			vecWord.push_back(temp_rank);
		}
		sort(vecWord.begin(), vecWord.end(), Comp);

		getline(inData, strBuffer);
		strtokenizer strTag(strBuffer," ");
		for( i=0; i<strTag.count_tokens(); i++)
		{
			rank_t temp_rank;
			strtokenizer strTerm(strTag.token(i),":");
			temp_rank.r_str = atoi(strTerm.token(0).c_str());
			temp_rank.r_val = atoi(strTerm.token(1).c_str());
			vecTrue.push_back(temp_rank);
		}
		sort(vecTrue.begin(), vecTrue.end(), Comp);

		// Inference
		q = dmatrix(cur_doc.len, nclass);
		if (q == NULL)
			return;

		EM(cur_doc, gamma, q, nt, ap, alpha, (const double **)beta, cur_doc.len, nclass, emmax, epsilon);
	
		double z = 0;
		int k;
		for( k = 0; k < nclass; k++)
			z += gamma[k];
		for( k = 0; k < nclass; k++)
		{
			gamma[k] = gamma[k]/z;
			cout << k << ":" << gamma[k] << " "<<endl;// nt[k] << "\t"; hzm
		}
		cout << endl << endl;

		free_dmatrix(q, cur_doc.len);

		// Evaluate
		for( i = 0; i < ntag; i++)
		{
			double cur_prob = 0.0;
			for(int k = 0; k < nclass; k++)
			{
				cur_prob += gamma[k]*eta[i][k];
			}
			rank_t temp_rank;
			temp_rank.r_str = i;
			temp_rank.r_val = cur_prob;
			vecModel.push_back(temp_rank);
		}
		sort(vecModel.begin(), vecModel.end(), Comp);

		// output
		/*int size = TARGET_SIZE < vecWord.size() ? TARGET_SIZE : vecWord.size();
		for(int i = 0; i < size; i++)
		{
			cout << m_vec_word[vecWord[i].r_str] << "\t";
		}
		cout << endl << endl;
		size = TARGET_SIZE < vecTrue.size() ? TARGET_SIZE : vecTrue.size();
		for(int i = 0; i < size; i++)
		{
			cout << m_vec_tag[vecTrue[i].r_str] << "\t";
		}
		cout << endl << endl;
		for(int i = 0; i < size; i++)
		{
			cout << m_vec_tag[vecModel[i].r_str] << "\t";
		}
		cout << endl << endl;*/

		map<int, int> mapTrue;
		map<int, int> mapModel;
		map<int, int>::iterator iter;
/*		int size = TARGET_SIZE < vecTrue.size() ? TARGET_SIZE : vecTrue.size();
		for( i = 0; i < size; i++)
		{
			mapTrue.insert(map<int, int>::value_type(vecTrue[i].r_str,0));
			mapModel.insert(map<int, int>::value_type(vecModel[i].r_str,0));
		}
		double precision = 0.0;
		double recall = 0.0;
		for(iter = mapTrue.begin(); iter != mapTrue.end(); iter++)
		{
			//cout << m_vec_tag[iter->first] << "\t";
			if(mapModel.find(iter->first) != mapModel.end())
				recall += 1;
		}
		//cout << endl << endl;
		for(iter = mapModel.begin(); iter != mapModel.end(); iter++)
		{
			//cout << m_vec_tag[iter->first] << "\t";
			if(mapTrue.find(iter->first) != mapTrue.end())
				precision += 1;
		}
		//cout << endl << endl;
		precision = precision/size;
		recall = recall/size;
		accuracy += precision;*/
	}

	cout << "average accuracy: " << accuracy/count << endl;
	
	free(gamma);
	free(ap);
	free(nt);
	free_dmatrix(beta, nlex);
	free(alpha);
	
	inData.close();
}

double CTagInference::EM(document &d, double *gamma, double **q, double *nt, double *ap, const double *alpha, const double **beta, int L, int K, int emmax, double epsilon)
{
	int j, k, l;
	double z;
	double lik = 0, plik = 0;

	for (k = 0; k < K; k++)
	{
		nt[k] = (double) d.size / K;
		gamma[k] = alpha[k] + nt[k];
	}
	
	for (j = 0; j < emmax; j++)
	{
		/* estep */
		for (k = 0; k < K; k++)
			ap[k] = exp(psi(gamma[k]));//guo
		/* accumulate q */
		for (l = 0; l < L; l++)
			for (k = 0; k < K; k++)
				q[l][k] = beta[d.id[l]][k] * ap[k];
		/* normalize q */
		for (l = 0; l < L; l++) {
			z = 0;
			for (k = 0; k < K; k++)
				z += q[l][k];
			for (k = 0; k < K; k++)
				q[l][k] /= z;
		}
		/* mstep */
		for (k = 0; k < K; k++) {
			z = 0;
			for (l = 0; l < L; l++)
				z += q[l][k] * d.cnt[l];
			nt[k] = z;
			gamma[k] = alpha[k] + nt[k];
		}
		/* modified by guo */
		lik = Likelihood(d, alpha, beta, gamma, q, L, K); 
		if ((j > 0) && (fabs((plik-lik)/plik) < epsilon))
		{
			//cout << "Converged!" << endl;
			break;
		}
		plik = lik;
	}
	return lik;
}

double CTagInference::Likelihood(document &d, const double *alpha, const double **beta, double *gamma, double** phi, int L, int K)
{
	double lik;
	int l, k;
	double sum_alpha, sum_gamma;
	lik = 0;

	for(k = 0, sum_alpha = 0, sum_gamma = 0; k < K; k++)
	{
		sum_alpha += alpha[k];
		sum_gamma += gamma[k];
	}
	
	lik += lgamma(sum_alpha) - lgamma(sum_gamma);

	for(k = 0; k < K; k++)
	{
		lik += lgamma(gamma[k]) - lgamma(alpha[k]);
	}
	
	for (l = 0; l < L; l++) {
		for (k = 0; k < K; k++) {
			if(phi[l][k] > 0.0)
			{
				/*if (beta[d.id[l]][k] == 0.0)
					lik += d.cnt[l] * phi[l][k] * ( -100.0 - log(phi[l][k]));
				else*/
				lik += d.cnt[l] * phi[l][k] * ( log(beta[d.id[l]][k]) - log(phi[l][k]));
			}
		}
	}

	return lik;
}

double CTagInference::digamma(double x)
{
  double result;
  static const double
	  neginf = -HUGE_VAL,//,-1.0/0.0,guo
	  c = 12,
	  s = 1e-6,
	  d1 = -0.57721566490153286,
	  d2 = 1.6449340668482264365, /* pi^2/6 */
	  s3 = 1./12,
	  s4 = 1./120,
	  s5 = 1./252,
	  s6 = 1./240,
	  s7 = 1./132,
	  s8 = 691/32760,
	  s9 = 1/12,
	  s10 = 3617/8160;
  /* Illegal arguments */
  if((x == neginf) || _isnan(x)) {
    return HUGE_VAL;//0.0/0.0;//guo
  }
  /* Singularities */
  if((x <= 0) && (floor(x) == x)) {
    return neginf;
  }
  /* Negative values */
  /* Use the reflection formula (Jeffrey 11.1.6):
   * digamma(-x) = digamma(x+1) + pi*cot(pi*x)
   *
   * This is related to the identity
   * digamma(-x) = digamma(x+1) - digamma(z) + digamma(1-z)
   * where z is the fractional part of x
   * For example:
   * digamma(-3.1) = 1/3.1 + 1/2.1 + 1/1.1 + 1/0.1 + digamma(1-0.1)
   *               = digamma(4.1) - digamma(0.1) + digamma(1-0.1)
   * Then we use
   * digamma(1-z) - digamma(z) = pi*cot(pi*z)
   */
  if(x < 0) {
    return digamma(1-x) + M_PI/tan(-M_PI*x);
  }
  /* Use Taylor series if argument <= S */
  if(x <= s) return d1 - 1/x + d2*x;
  /* Reduce to digamma(X + N) where (X + N) >= C */
  result = 0;
  while(x < c) {
    result -= 1/x;
    x++;
  }
  /* Use de Moivre's expansion if argument >= C */
  /* This expansion can be computed in Maple via asympt(Psi(x),x) */
  if(x >= c) {
    double r = 1/x;
    result += log(x) - 0.5*r;
    r *= r;
    result -= r * (s3 - r * (s4 - r * (s5 - r * (s6 - r * s7))));
  }
  return result;
}

double CTagInference::lgamma(double x)
{
    double x0,x2,xp,gl,gl0;
    int n,k;
    static double a[] = {
        8.333333333333333e-02,
       -2.777777777777778e-03,
        7.936507936507937e-04,
       -5.952380952380952e-04,
        8.417508417508418e-04,
       -1.917526917526918e-03,
        6.410256410256410e-03,
       -2.955065359477124e-02,
        1.796443723688307e-01,
       -1.39243221690590
    };
    
    x0 = x;
    if (x <= 0.0) return 1e308;
    else if ((x == 1.0) || (x == 2.0)) return 0.0;
    else if (x <= 7.0) {
        n = (int)(7-x);
        x0 = x+n;
    }
    x2 = 1.0/(x0*x0);
    xp = 2.0*M_PI;
    gl0 = a[9];
    for (k=8;k>=0;k--) {
        gl0 = gl0*x2 + a[k];
    }
    gl = gl0/x0+0.5*log(xp)+(x0-0.5)*log(x0)-x0;
    if (x <= 7.0) {
        for (k=1;k<=n;k++) {
            gl -= log(x0-1.0);
            x0 -= 1.0;
        }
    }
    return gl;
}

double** CTagInference::dmatrix(int rows, int cols)
{
	double **matrix;
	int i;

	matrix = (double **)calloc(rows, sizeof(double));
	if (matrix == NULL)
		return NULL;
	for (i = 0; i < rows; i++) {
		matrix[i] = (double *)calloc(cols, sizeof(double));
		if (matrix[i] == NULL)
			return NULL;
	}
	
	return matrix;
}

void CTagInference::free_dmatrix(double **matrix, int rows)
{
	int i;
	for (i = 0; i < rows; i++)
		free(matrix[i]);
	free(matrix);
}

bool CTagInference::Comp( rank_t &e1, rank_t &e2 )
{
	return e1.r_val > e2.r_val;
}