Created
March 15, 2013 10:23
-
-
Save tomonari-masada/5168843 to your computer and use it in GitHub Desktop.
A Revised Inference for Correlated Topic Model
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #include <stdio.h> | |
| #include <stdlib.h> | |
| #include <string.h> | |
| #include <math.h> | |
| #include <time.h> | |
| #include <omp.h> | |
| #include <xmmintrin.h> | |
| #include <pmmintrin.h> | |
| #include <lbfgs.h> | |
| #define DIAGONAL_RATIO 0.0 | |
| #define REG_CONST 1.0 | |
| #define BUFFSIZE 100000 | |
| #define TOKENLEN 256 | |
| #define HASH_SIZE 65537 | |
| #define DTYPE double | |
| #define MIN_HYP 1e-3 | |
| #define HYP_CONV 1e-3 | |
| #define NUM_OF_CGS_ITERATIONS 1000 | |
| #define NUM_OF_THREADS 16 | |
| #define NUM_OF_VB_ITERATIONS 200 | |
| FILE *fOut; | |
| int *nWordNums; | |
| int **nWordIDs; | |
| int nNumOfTopics; | |
| int nNumOfTopics2; | |
| int nNumOfWords; | |
| char **sWords; | |
| int nNumOfDocs; | |
| int *nDocHeadCell; | |
| int *nDocNumOfCells; | |
| int *nDocLen; | |
| char **sDocName; | |
| int nMaxDocNumOfCells; | |
| int nNumOfTokens; | |
| int nNumOfTestTokens; | |
| int nNumOfCells; | |
| int *nWordCell; | |
| int *nTFCell; | |
| int *nTestTFCell; | |
| int **nTopicArray; | |
| DTYPE *dNu; | |
| DTYPE *dMeanR; | |
| DTYPE *dStdevS; | |
| DTYPE *dInvSigma; | |
| DTYPE dLogDetSigma; | |
| DTYPE *dMu; | |
| DTYPE *dLambda; | |
| DTYPE *dLambdaSum; | |
| DTYPE *dZeta; | |
| DTYPE *dZetaSum; | |
| DTYPE *dZeta2; | |
| DTYPE *dZetaSum2; | |
| DTYPE *dAlpha; | |
| DTYPE dAlphaSum; | |
| DTYPE *dBeta; | |
| DTYPE dBetaSum; | |
| DTYPE dDiagonal; | |
| //-------- | |
| typedef struct { int nID; DTYPE dData; } sortData; | |
| int sort_compdata(const void *a, const void *b) | |
| { | |
| if (((sortData *) a)->dData < ((sortData *) b)->dData) { return 1; } | |
| else if (((sortData *) a)->dData > ((sortData *) b)->dData) { return -1; } | |
| else { return 0; } | |
| } | |
| //-------- | |
| time_t start_time; | |
| void put_time(FILE *f) | |
| { | |
| time_t current_time; | |
| time(¤t_time); | |
| fprintf(f, "%.0lf", difftime(current_time, start_time)); | |
| fflush(f); | |
| } | |
| //-------- | |
| DTYPE dDigammaOne; | |
| DTYPE digamma(DTYPE dX) | |
| { | |
| DTYPE dXX, dXX2, dTemp, dSum; | |
| dXX = dX; dTemp = 0.0; dSum = 0.0; | |
| while (dXX < 8.0) { dTemp += 1.0 / dXX; dXX = dXX + 1.0; } | |
| dXX2 = dXX * dXX; | |
| dSum = 1.0 / (120.0 * dXX2) - 1.0 / 12.0; | |
| dSum = dSum / dXX2 - 1.0 / (2.0 * dXX) + logl(dXX) - dTemp; | |
| return dSum; | |
| } | |
| DTYPE trigamma(DTYPE dX) | |
| { | |
| DTYPE dXX, dXX2, dTemp, dSum; | |
| dXX = dX; dTemp = 0.0; dSum = 0.0; | |
| while (dXX < 8.0) { dTemp -= 1.0 / (dXX * dXX); dXX = dXX + 1.0; } | |
| dXX2 = dXX * dXX; | |
| dSum = - 1.0 / (30.0 * dXX2) + 1.0 / 6.0; | |
| dSum = (dSum / dXX2 + 1.0 / (2.0 * dXX) + 1.0) / dXX - dTemp; | |
| return dSum; | |
| } | |
| DTYPE tetragamma(DTYPE dX) | |
| { | |
| DTYPE dXX, dXX2, dTemp, dSum; | |
| dXX = dX; dTemp = 0.0; dSum = 0.0; | |
| while (dXX < 8.0) { dTemp += 2.0 / (dXX * dXX * dXX); dXX = dXX + 1.0; } | |
| dXX2 = dXX * dXX; | |
| dSum = (1.0 / 6.0) / dXX2 - 1.0 / 2.0; | |
| dSum = (dSum / dXX2 - 1.0 / dXX - 1.0) / dXX2 - dTemp; | |
| return dSum; | |
| } | |
| DTYPE invDigamma(DTYPE dY) | |
| { | |
| DTYPE dX; | |
| int nCount; | |
| if (dY >= -2.22) dX = exp(dY) + 0.5; | |
| else dX = - 1.0 / (dY - dDigammaOne); | |
| for (nCount = 0; nCount < 5; nCount ++) | |
| dX = dX - (digamma(dX) - dY) / trigamma(dX); | |
| return dX; | |
| } | |
| //-------- MT | |
| /* Period parameters */ | |
| #define N 624 | |
| #define M 397 | |
| #define MATRIX_A 0x9908b0dfUL /* constant vector a */ | |
| #define UPPER_MASK 0x80000000UL /* most significant w-r bits */ | |
| #define LOWER_MASK 0x7fffffffUL /* least significant r bits */ | |
| unsigned long *mt; /* the array for the state vector */ | |
| int *mti; /* mti==N+1 means mt[N] is not initialized */ | |
| void init_genrand(unsigned long s, int nThreadID) | |
| { | |
| int nTemp; | |
| mt[nThreadID * N + 0]= s & 0xffffffffUL; | |
| for (mti[nThreadID] = 1; mti[nThreadID] < N; mti[nThreadID] ++) { | |
| nTemp = nThreadID * N + mti[nThreadID]; | |
| mt[nTemp] = (1812433253UL * (mt[nTemp - 1] ^ (mt[nTemp - 1] >> 30)) + mti[nThreadID]); | |
| mt[nTemp] &= 0xffffffffUL; | |
| } | |
| } | |
| unsigned long genrand_int32(int nThreadID) | |
| { | |
| unsigned long y; | |
| unsigned long mag01[2]={0x0UL, MATRIX_A}; | |
| int nTemp; | |
| if (mti[nThreadID] >= N) { /* generate N words at one time */ | |
| int kk; | |
| if (mti[nThreadID] == N+1) /* if init_genrand() has not been called, */ | |
| init_genrand(5489UL, nThreadID); /* a default initial seed is used */ | |
| for (kk = 0; kk < N - M; kk ++) { | |
| nTemp = nThreadID * N + kk; | |
| y = (mt[nTemp] & UPPER_MASK) | (mt[nTemp + 1] & LOWER_MASK); | |
| mt[nTemp] = mt[nTemp + M] ^ (y >> 1) ^ mag01[y & 0x1UL]; | |
| } | |
| for (; kk < N - 1; kk ++) { | |
| nTemp = nThreadID * N + kk; | |
| y = (mt[nTemp] & UPPER_MASK) | (mt[nTemp + 1] & LOWER_MASK); | |
| mt[nTemp] = mt[nTemp + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1UL]; | |
| } | |
| nTemp = nThreadID * N; | |
| y = (mt[nTemp + N - 1] & UPPER_MASK) | (mt[nTemp] & LOWER_MASK); | |
| mt[nTemp + N - 1] = mt[nTemp + M - 1] ^ (y >> 1) ^ mag01[y & 0x1UL]; | |
| mti[nThreadID] = 0; | |
| } | |
| y = mt[nThreadID * N + mti[nThreadID]]; | |
| mti[nThreadID] ++; | |
| y ^= (y >> 11); | |
| y ^= (y << 7) & 0x9d2c5680UL; | |
| y ^= (y << 15) & 0xefc60000UL; | |
| y ^= (y >> 18); | |
| return y; | |
| } | |
| long genrand_int31(int nT) { return (long)(genrand_int32(nT)>>1); } | |
| double genrand_real1(int nT) { return genrand_int32(nT)*(1.0/4294967295.0); /* divided by 2^32-1 */ } | |
| double genrand_real2(int nT) { return genrand_int32(nT)*(1.0/4294967296.0); /* divided by 2^32 */ } | |
| double genrand_real3(int nT) { return (((double)genrand_int32(nT)) + 0.5)*(1.0/4294967296.0); /* divided by 2^32 */ } | |
| double genrand_res53(int nT) { unsigned long a=genrand_int32(nT)>>5, b=genrand_int32(nT)>>6; return(a*67108864.0+b)*(1.0/9007199254740992.0); } | |
| void initRand() | |
| { | |
| int nCount; | |
| mt = (unsigned long *) malloc(sizeof(unsigned long) * NUM_OF_THREADS * N); /* the array for the state vector */ | |
| mti = (int *) malloc(sizeof(int) * NUM_OF_THREADS); | |
| init_genrand((unsigned) time(NULL), 0); | |
| for (nCount = 1; nCount < NUM_OF_THREADS; nCount ++) | |
| init_genrand(genrand_int31(nCount - 1), nCount); | |
| return; | |
| } | |
| //-------- | |
| int hash_func(const char *sStr) | |
| { | |
| int n1, n2, n3, n4, nL; | |
| nL = (int) strlen(sStr); if (nL < 0) nL = - nL; n1 = (int) sStr[0]; | |
| if (strlen(sStr) > 1) { n2 = (int) sStr[1]; n3 = (int) sStr[nL - 2]; n4 = (int) sStr[nL - 1]; } | |
| else { n2 = (int) sStr[0]; n3 = (int) sStr[0]; n4 = (int) sStr[0]; } | |
| if (n1 < 0) n1 = - n1; if (n2 < 0) n2 = - n2; if (n3 < 0) n3 = - n3; if (n4 < 0) n4 = - n4; | |
| return ((n2 + n1 * 1025 + nL * 4097 + n3 * 32769) % HASH_SIZE + (n4 * 65537) % HASH_SIZE) % HASH_SIZE; | |
| } | |
| int initHash() | |
| { | |
| int nCount; | |
| nNumOfWords = 0; | |
| nWordNums = (int *) malloc(sizeof(int) * HASH_SIZE); | |
| nWordIDs = (int **) malloc(sizeof(int *) * HASH_SIZE); | |
| if ((nWordNums == NULL) || (nWordIDs == NULL)) return 1; | |
| for (nCount = 0; nCount < HASH_SIZE; nCount ++) nWordNums[nCount] = 0; | |
| return 0; | |
| } | |
| int addWord(const char *s) | |
| { | |
| unsigned nCount; | |
| int nID; | |
| int nHashKey; | |
| nID = 0; | |
| nHashKey = hash_func(s); | |
| for (nCount = 0; nCount < nWordNums[nHashKey]; nCount ++) | |
| if (strcmp(s, sWords[nWordIDs[nHashKey][nCount] - 1]) == 0) { | |
| nID = nWordIDs[nHashKey][nCount]; | |
| break; | |
| } | |
| if (nID == 0) { | |
| nNumOfWords ++; | |
| nWordNums[nHashKey] ++; | |
| if (nWordNums[nHashKey] == 1) | |
| nWordIDs[nHashKey] = (int *) malloc(sizeof(int)); | |
| else | |
| nWordIDs[nHashKey] = (int *) realloc(nWordIDs[nHashKey], sizeof(int) * nWordNums[nHashKey]); | |
| nWordIDs[nHashKey][nWordNums[nHashKey] - 1] = nNumOfWords; | |
| if (nNumOfWords == 1) | |
| sWords = (char **) malloc(sizeof(char *)); | |
| else | |
| sWords = (char **) realloc(sWords, sizeof(char *) * nNumOfWords); | |
| sWords[nNumOfWords - 1] = (char *) malloc(sizeof(char) * (strlen(s) + 1)); | |
| memset(sWords[nNumOfWords - 1], 0, strlen(s) + 1); | |
| strcpy(sWords[nNumOfWords - 1], s); | |
| nID = nNumOfWords; | |
| } | |
| return nID; | |
| } | |
| //-------- matrix calculations | |
| DTYPE inverse(DTYPE *dMat, int nDim) | |
| { | |
| int nC, nOrder, nLdim, *nPivot, nLwork, nInfo; | |
| DTYPE dRet, *dWork; | |
| nOrder = nDim; nLdim = nDim; | |
| nLwork = nDim * nDim; | |
| nPivot = (int *) malloc(sizeof(int) * (nDim + 1)); | |
| dWork = (DTYPE *) malloc(sizeof(DTYPE) * nLwork); | |
| dgetrf_(&nOrder, &nOrder, dMat, &nLdim, nPivot, &nInfo); | |
| if (nInfo != 0) { fprintf(stderr, "fatal error in inverse %d (%d)\n",nInfo,nDim); exit(1); } | |
| dRet = 0.0; | |
| for (nC = 0; nC < nDim; nC ++) dRet += log(fabs(dMat[nC + nDim * nC])); | |
| dgetri_(&nOrder, dMat, &nLdim, nPivot, dWork, &nLwork, &nInfo); | |
| free(nPivot); | |
| free(dWork); | |
| return dRet; | |
| } | |
| //-------- | |
| //DTYPE dConst = REG_CONST * DIAGONAL_RATIO; | |
| DTYPE dConst = REG_CONST; | |
| static lbfgsfloatval_t LBFGS_MeanR( | |
| void *instance, | |
| const lbfgsfloatval_t *dX, | |
| lbfgsfloatval_t *dG, | |
| const int nDim, | |
| const lbfgsfloatval_t dStep | |
| ) | |
| { | |
| int nCD, nCK, nCK2; | |
| int nTemp, nTemp2; | |
| lbfgsfloatval_t dFX; | |
| DTYPE dTemp; | |
| nCD = (*((int *) instance)); | |
| nTemp = nCD * nNumOfTopics; | |
| dFX = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| nTemp2 = nCK * nNumOfTopics; | |
| dG[nCK] = 0.0; | |
| for (nCK2 = 0; nCK2 < nNumOfTopics; nCK2 ++) | |
| dG[nCK] += (dX[nCK2] - dMu[nCK2]) * dInvSigma[nTemp2 + nCK2]; | |
| dTemp = exp(dX[nCK] | |
| + 0.5 * dX[nNumOfTopics + nCK] * dX[nNumOfTopics + nCK]); | |
| dTemp *= ((DTYPE) nDocLen[nCD]) / dNu[nCD]; | |
| dFX -= dX[nCK] * dG[nCK]; | |
| dFX += 0.5 * dX[nCK] * dX[nCK] * dInvSigma[nTemp2 + nCK]; | |
| dFX += dX[nCK] * dLambda[nTemp + nCK]; | |
| dFX -= dTemp; | |
| dFX -= 0.5 * dX[nNumOfTopics + nCK] * dX[nNumOfTopics + nCK] | |
| * dInvSigma[nTemp2 + nCK]; | |
| dFX += log(fabs(dX[nNumOfTopics + nCK])); | |
| dG[nCK] = dLambda[nTemp + nCK] - dTemp - dG[nCK]; | |
| dG[nNumOfTopics + nCK] = - dTemp | |
| - fabs(dX[nNumOfTopics + nCK]) * dInvSigma[nTemp2 + nCK] | |
| + 1.0 / fabs(dX[nNumOfTopics + nCK]); | |
| dTemp = dX[nCK] - log(dLambda[nTemp + nCK] + dAlpha[nCK]); | |
| dFX -= dConst * dTemp * dTemp; | |
| dG[nCK] -= 2.0 * dConst * dTemp; | |
| } | |
| dFX = - dFX; | |
| for (nCK = 0; nCK < nNumOfTopics * 2; nCK ++) dG[nCK] = - dG[nCK]; | |
| return dFX; | |
| } | |
| static int LBFGS_progress( | |
| void *instance, | |
| const lbfgsfloatval_t *x, | |
| const lbfgsfloatval_t *g, | |
| const lbfgsfloatval_t fx, | |
| const lbfgsfloatval_t xnorm, | |
| const lbfgsfloatval_t gnorm, | |
| const lbfgsfloatval_t step, | |
| int n, | |
| int k, | |
| int ls | |
| ) | |
| { | |
| /* | |
| printf("Iteration %d:\n", k); | |
| printf(" fx = %f, x[0] = %f, x[1] = %f\n", fx, x[0], x[1]); | |
| printf(" xnorm = %f, gnorm = %f, step = %f\n", xnorm, gnorm, step); | |
| printf("\n"); | |
| */ | |
| if (k > 1000) return 1; | |
| return 0; | |
| } | |
| //-------- | |
| int statistics() | |
| { | |
| int nCD, nCK, nCW, nCC, nCToken; | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) { | |
| dLambdaSum[nCD] = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dLambda[nCD * nNumOfTopics + nCK] = 0.0; | |
| } | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| dZetaSum[nCK] = 0.0; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) | |
| dZeta[nCK * nNumOfWords + nCW] = 0.0; | |
| } | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) | |
| for (nCC = nDocHeadCell[nCD]; nCC < nDocHeadCell[nCD] + nDocNumOfCells[nCD]; nCC ++) | |
| for (nCToken = 0; nCToken < nTFCell[nCC]; nCToken ++) { | |
| dLambda[nCD * nNumOfTopics + nTopicArray[nCC][nCToken]] += 1.0; | |
| dZeta[nTopicArray[nCC][nCToken] * nNumOfWords + nWordCell[nCC]] += 1.0; | |
| } | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dLambdaSum[nCD] += dLambda[nCD * nNumOfTopics + nCK]; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) | |
| dZetaSum[nCK] += dZeta[nCK * nNumOfWords + nCW]; | |
| return 0; | |
| } | |
| //-------- updating hyper parameters | |
| int optHypA2(FILE *fOut, int nI) | |
| { | |
| int nCK, nIter; | |
| DTYPE dPrevAlphaSum; | |
| DTYPE dConv; | |
| fprintf(fOut, "#hA\t%d", nI); | |
| dAlphaSum = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dAlphaSum += dAlpha[nCK]; | |
| nIter = 0; | |
| do { | |
| #pragma omp parallel | |
| { | |
| int nCD, nCK; | |
| DTYPE dNSum, dDSum; | |
| #ifdef _OPENMP | |
| int nThreadID = omp_get_thread_num(); | |
| int nNumOfThreads = omp_get_num_threads(); | |
| #else | |
| int nThreadID = 0; | |
| int nNumOfThreads = 1; | |
| #endif | |
| for (nCK = nThreadID; nCK < nNumOfTopics; nCK += nNumOfThreads) { | |
| dNSum = 0.0; | |
| dDSum = 0.0; | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) { | |
| dNSum += digamma(dLambda[nCD * nNumOfTopics + nCK] + dAlpha[nCK]); | |
| dDSum += digamma(dLambdaSum[nCD] + dAlphaSum); | |
| } | |
| dNSum = dAlpha[nCK] * (dNSum - digamma(dAlpha[nCK]) * ((DTYPE) nNumOfDocs)); | |
| dDSum = 1.0 + dDSum - digamma(dAlphaSum) * ((DTYPE) nNumOfDocs); | |
| dAlpha[nCK] = dNSum / dDSum; | |
| if (dAlpha[nCK] < MIN_HYP) dAlpha[nCK] = MIN_HYP; | |
| } | |
| } | |
| nIter ++; | |
| dPrevAlphaSum = dAlphaSum; | |
| dAlphaSum = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dAlphaSum += dAlpha[nCK]; | |
| dConv = fabs((dAlphaSum - dPrevAlphaSum) / dAlphaSum); | |
| } while (dConv > HYP_CONV); | |
| fprintf(fOut, "\t%.6lf\t%d\n", dAlphaSum / ((DTYPE) nNumOfTopics), nIter); | |
| return 0; | |
| } | |
| int optHypB2(FILE *fOut, int nI) | |
| { | |
| int nCW, nIter; | |
| DTYPE dPrevBetaSum; | |
| DTYPE dConv; | |
| fprintf(fOut, "#hB\t%d", nI); | |
| dBetaSum = 0.0; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) dBetaSum += dBeta[nCW]; | |
| nIter = 0; | |
| do { | |
| #pragma omp parallel | |
| { | |
| int nCK, nCW; | |
| DTYPE dNSum, dDSum; | |
| #ifdef _OPENMP | |
| int nThreadID = omp_get_thread_num(); | |
| int nNumOfThreads = omp_get_num_threads(); | |
| #else | |
| int nThreadID = 0; | |
| int nNumOfThreads = 1; | |
| #endif | |
| for (nCW = nThreadID; nCW < nNumOfWords; nCW += nNumOfThreads) { | |
| dNSum = 0.0; | |
| dDSum = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| dNSum += digamma(dZeta[nCK * nNumOfWords + nCW] + dBeta[nCW]); | |
| dDSum += digamma(dZetaSum[nCK] + dBetaSum); | |
| } | |
| dNSum = dBeta[nCW] * (dNSum - digamma(dBeta[nCW]) * ((DTYPE) nNumOfTopics)); | |
| dDSum = 1.0 + dDSum - digamma(dBetaSum) * ((DTYPE) nNumOfTopics); | |
| dBeta[nCW] = dNSum / dDSum; | |
| if (dBeta[nCW] < MIN_HYP) dBeta[nCW] = MIN_HYP; | |
| } | |
| } | |
| nIter ++; | |
| dPrevBetaSum = dBetaSum; | |
| dBetaSum = 0.0; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) dBetaSum += dBeta[nCW]; | |
| dConv = fabs((dBetaSum - dPrevBetaSum) / dBetaSum); | |
| } while (dConv > HYP_CONV); | |
| fprintf(fOut, "\t%.6lf\t%d\n", dBetaSum / ((DTYPE) nNumOfWords), nIter); | |
| return 0; | |
| } | |
| //-------- for evaluation | |
| double logEviLB() | |
| { | |
| double dRet; | |
| dRet = 0.0; | |
| #pragma omp parallel | |
| { | |
| double dLocalRet; | |
| double dLR; | |
| int nCD, nCK, nCK2, nCW, nCC; | |
| int nTemp, nTemp3; | |
| //int nCD2; | |
| double *dPi; | |
| double dTemp; | |
| #ifdef _OPENMP | |
| int nThreadID = omp_get_thread_num(); | |
| int nNumOfThreads = omp_get_num_threads(); | |
| #else | |
| int nThreadID = 0; | |
| int nNumOfThreads = 1; | |
| #endif | |
| dPi = (double *) malloc(sizeof(double) * nNumOfTopics); | |
| #pragma omp barrier | |
| dLocalRet = 0.0; | |
| for (nCD = nThreadID; nCD < nNumOfDocs; nCD += nNumOfThreads) { | |
| nTemp = nCD * nNumOfTopics; | |
| for (nCC = nDocHeadCell[nCD]; nCC < nDocHeadCell[nCD] + nDocNumOfCells[nCD]; nCC ++) { | |
| nCW = nWordCell[nCC]; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dPi[nCK] = dMeanR[nTemp + nCK] | |
| + dZeta2[nCK * nNumOfWords + nCW] - dZetaSum2[nCK]; | |
| dTemp = dPi[0]; | |
| for (nCK = 1; nCK < nNumOfTopics; nCK ++) if (dTemp < dPi[nCK]) dTemp = dPi[nCK]; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dPi[nCK] = exp(dPi[nCK] - dTemp); | |
| dTemp = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dTemp += dPi[nCK]; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dPi[nCK] /= dTemp; | |
| dLR = dLocalRet; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dLocalRet += ((double) nTFCell[nCC]) * dPi[nCK] | |
| * (dMeanR[nTemp + nCK] - log(dNu[nCD]) | |
| + dZeta2[nCK * nNumOfWords + nCW] - dZetaSum2[nCK]); | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| if (dPi[nCK] > 1e-5) | |
| dLocalRet -= ((double) nTFCell[nCC]) * dPi[nCK] * log(dPi[nCK]); | |
| {if((!isnan(dLR))&&(isnan(dLocalRet))){printf("####00\t%d\t%d\t%.8lf\t%.8lf\n",nCD,nCC,dLR,dLocalRet);fflush(stdout);}} | |
| } | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| nTemp = nCD * nNumOfTopics + nCK; | |
| nTemp3 = nCK * nNumOfTopics + nCK; | |
| dLR = dLocalRet; | |
| dLocalRet -= 0.5 * (dStdevS[nTemp] * dStdevS[nTemp] * dInvSigma[nTemp3]); | |
| {if((!isnan(dLR))&&(isnan(dLocalRet))){printf("####01a\t%d\t%d\t%.8lf\t%.8lf\n",nCD,nCK,dLR,dLocalRet);fflush(stdout);}} | |
| dLR = dLocalRet; | |
| dLocalRet += log(dStdevS[nTemp]); | |
| {if((!isnan(dLR))&&(isnan(dLocalRet))){printf("####01b\t%d\t%d\t%.8lf\t%.8lf\t%.8lf\n",nCD,nCK,dLR,dLocalRet,dStdevS[nTemp]);fflush(stdout);}} | |
| dLR = dLocalRet; | |
| dLocalRet -= 0.5 * (dMeanR[nTemp] - dMu[nCK]) * (dMeanR[nTemp] - dMu[nCK]) * dInvSigma[nTemp3]; | |
| {if((!isnan(dLR))&&(isnan(dLocalRet))){printf("####01d\t%d\t%d\t%.8lf\t%.8lf\n",nCD,nCK,dLR,dLocalRet);fflush(stdout);}} | |
| nTemp = nCD * nNumOfTopics; | |
| nTemp3 = nCK * nNumOfTopics; | |
| dLR = dLocalRet; | |
| for (nCK2 = nCK + 1; nCK2 < nNumOfTopics; nCK2 ++) | |
| dLocalRet -= (dMeanR[nTemp + nCK] - dMu[nCK]) * (dMeanR[nTemp + nCK2] - dMu[nCK2]) * dInvSigma[nTemp3 + nCK2]; | |
| {if((!isnan(dLR))&&(isnan(dLocalRet))){printf("####02\t%d\t%d\t%.8lf\t%.8lf\n",nCD,nCK,dLR,dLocalRet);fflush(stdout);}} | |
| } | |
| } | |
| #pragma omp barrier | |
| for (nCK = nThreadID; nCK < nNumOfTopics; nCK += nNumOfThreads) { | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) { | |
| nTemp = nCK * nNumOfWords + nCW; | |
| dLR = dLocalRet; | |
| dLocalRet += (dBeta[nCW] - dZeta[nTemp]) * (dZeta2[nTemp] - dZetaSum2[nCK]) + lgamma(dZeta[nTemp]); | |
| {if((!isnan(dLR))&&(isnan(dLocalRet))){printf("####04\t%d\t%s\t%.8lf\t%.8lf\n",nCK,sWords[nCW],dLR,dLocalRet);fflush(stdout);}} | |
| {if((!isinf(dLR))&&(isinf(dLocalRet))){printf("#---04\t%d\t%s\t%.8lf\t%.8lf\t%.8lf\n",nCK,sWords[nCW],dLR,dLocalRet,dZeta[nTemp]);fflush(stdout);}} | |
| } | |
| dLR = dLocalRet; | |
| dLocalRet -= lgamma(dZetaSum[nCK]); | |
| {if((!isnan(dLR))&&(isnan(dLocalRet))){printf("####05\t%d\t%.8lf\t%.8lf\n",nCK,dLR,dLocalRet);fflush(stdout);}} | |
| {if((!isinf(dLR))&&(isinf(dLocalRet))){printf("#---05\t%d\t%.8lf\t%.8lf\n",nCK,dLR,dLocalRet);fflush(stdout);}} | |
| } | |
| #pragma omp barrier | |
| for (nCW = nThreadID; nCW < nNumOfWords; nCW += nNumOfThreads) { | |
| dLR = dLocalRet; | |
| dLocalRet -= ((DTYPE) nNumOfTopics) * lgamma(dBeta[nCW]); | |
| {if((!isnan(dLR))&&(isnan(dLocalRet))){printf("####06\t%s\t%.8lf\t%.8lf\n",sWords[nCW],dLR,dLocalRet);fflush(stdout);}} | |
| {if((!isinf(dLR))&&(isinf(dLocalRet))){printf("#---06\t%s\t%.8lf\t%.8lf\n",sWords[nCW],dLR,dLocalRet);fflush(stdout);}} | |
| } | |
| #pragma omp critical | |
| { | |
| dRet += dLocalRet; | |
| } | |
| free(dPi); | |
| } | |
| dRet += ((DTYPE) nNumOfTopics) * lgamma(dBetaSum); | |
| dRet += ((DTYPE) nNumOfDocs) * ((DTYPE) nNumOfTopics); | |
| // Determinants cannot be calculated correctly. | |
| //dPrevRet = dRet; | |
| dRet -= 0.5 * ((DTYPE) nNumOfDocs) * dLogDetSigma; | |
| //{if((!isnan(dPrevRet))&&(isnan(dRet))){printf("####07\t%.8lf\t%.8lf\n",dPrevRet,dRet);fflush(stdout);}} | |
| //dPrevRet = dRet; | |
| //dRet -= 0.5 * ((DTYPE) nNumOfTopics) * dLogDetGram; | |
| //{if((!isnan(dPrevRet))&&(isnan(dRet))){printf("####08\t%.8lf\t%.8lf\n",dPrevRet,dRet);fflush(stdout);}} | |
| return dRet; | |
| } | |
| double perplexityLDA(int nFlag) // call setTempDocTopic before calling this function when nFlag != 0 | |
| { | |
| double dRet = 0.0; | |
| #pragma omp parallel | |
| { | |
| double dLocalRet; | |
| double dProb; | |
| double dTempSum; | |
| int nCD, nCK, nCC, nCW; | |
| #ifdef _OPENMP | |
| int nThreadID = omp_get_thread_num(); | |
| int nNumOfThreads = omp_get_num_threads(); | |
| #else | |
| int nThreadID = 0; | |
| int nNumOfThreads = 1; | |
| #endif | |
| dLocalRet = 0.0; | |
| for (nCD = nThreadID; nCD < nNumOfDocs; nCD += nNumOfThreads) { | |
| dTempSum = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dTempSum += exp(dMeanR[nCD * nNumOfTopics + nCK]); | |
| for (nCC = nDocHeadCell[nCD]; nCC < nDocHeadCell[nCD] + nDocNumOfCells[nCD]; nCC ++) | |
| if (nTestTFCell[nCC] > 0) { | |
| nCW = nWordCell[nCC]; | |
| dProb = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| if (nFlag == 0) | |
| dProb += (dLambda[nCD * nNumOfTopics + nCK] + dAlpha[nCK]) | |
| * (dZeta[nCK * nNumOfWords + nCW] + dBeta[nCW]) / (dZetaSum[nCK] + dBetaSum); | |
| else | |
| dProb += exp(dMeanR[nCD * nNumOfTopics + nCK]) | |
| * dZeta[nCK * nNumOfWords + nCW] / dZetaSum[nCK]; | |
| if (nFlag == 0) | |
| dProb /= (dLambdaSum[nCD] + dAlphaSum); | |
| else | |
| dProb /= dTempSum; | |
| dLocalRet += ((double) nTestTFCell[nCC]) * log(dProb); | |
| } | |
| } | |
| #pragma omp critical | |
| { | |
| dRet -= dLocalRet; | |
| } | |
| } | |
| return exp(dRet / ((double) nNumOfTestTokens)); | |
| } | |
| int setTestData() | |
| { | |
| int nCD, nCC; | |
| int nCToken; | |
| nNumOfTestTokens = 0; | |
| nTestTFCell = (int *) malloc(sizeof(int) * nNumOfCells); | |
| nTopicArray = (int **) malloc(sizeof(int *) * nNumOfCells); | |
| for (nCC = 0; nCC < nNumOfCells; nCC ++) nTestTFCell[nCC] = 0; | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) | |
| for (nCC = nDocHeadCell[nCD]; nCC < nDocHeadCell[nCD] + nDocNumOfCells[nCD]; nCC ++) { | |
| nTestTFCell[nCC] = 0; | |
| for (nCToken = 0; nCToken < nTFCell[nCC]; nCToken ++) | |
| if (genrand_int32(nCC % NUM_OF_THREADS) % 10 == 0) nTestTFCell[nCC] ++; | |
| nTFCell[nCC] -= nTestTFCell[nCC]; | |
| nNumOfTestTokens += nTestTFCell[nCC]; | |
| nNumOfTokens -= nTestTFCell[nCC]; | |
| if (nTFCell[nCC] > 0) nTopicArray[nCC] = (int *) malloc(sizeof(int) * nTFCell[nCC]); | |
| } | |
| return 0; | |
| } | |
| int printLog(int nI) | |
| { | |
| sortData *stData; | |
| int nCD, nCK, nCW; | |
| stData = (sortData *) malloc(sizeof(sortData) * nNumOfWords); | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) { | |
| stData[nCW].nID = nCW; | |
| stData[nCW].dData = dZeta[nCK * nNumOfWords + nCW]; | |
| } | |
| qsort(stData, nNumOfWords, sizeof(sortData), sort_compdata); | |
| printf("#KW %d %d", nI, nCK); | |
| for (nCW = 0; nCW < 100; nCW ++) | |
| if (stData[nCW].dData > 0.0) | |
| printf(" %s %.4le", sWords[stData[nCW].nID], stData[nCW].dData); | |
| printf("\n"); | |
| fflush(stdout); | |
| } | |
| free(stData); | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) { | |
| printf("#JK %d %s", nI, sDocName[nCD]); | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| printf(" %.4le", dLambda[nCD * nNumOfTopics + nCK]); | |
| printf("\n"); | |
| } | |
| if (nI > 0) | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) { | |
| printf("#Me %d %s", nI, sDocName[nCD]); | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| printf(" %.4le", exp(dMeanR[nCD * nNumOfTopics + nCK])); | |
| printf("\n"); | |
| } | |
| return 0; | |
| } | |
| //-------- | |
| int main(int argc, char **argv) | |
| { | |
| char sBuff[BUFFSIZE]; | |
| int nDocID, nWordID; | |
| char sTempDoc[TOKENLEN]; | |
| char sPrevDoc[TOKENLEN]; | |
| char sTempWord[TOKENLEN]; | |
| int nTF; | |
| int nIteration; | |
| if (argc < 2) { | |
| fprintf(stderr, "usage : %s num_of_topics\n", argv[0]); | |
| exit(1); | |
| } | |
| nNumOfTopics = atoi(argv[1]); | |
| if (nNumOfTopics < 2) { fprintf(stderr, "%s : illegal argument : %s\n", argv[0], argv[1]); exit(1); } | |
| // denormal values | |
| _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON); | |
| _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); | |
| fOut = stdout; | |
| initHash(); | |
| initRand(); | |
| nNumOfDocs = 0; | |
| nDocHeadCell = (int *) malloc(sizeof(int)); | |
| nDocNumOfCells = (int *) malloc(sizeof(int)); | |
| nDocLen = (int *) malloc(sizeof(int)); | |
| sDocName= (char **) malloc(sizeof(char *)); | |
| nNumOfCells = 0; | |
| nWordCell = (int *) malloc(sizeof(int)); | |
| nTFCell = (int *) malloc(sizeof(int)); | |
| memset(sPrevDoc, 0, TOKENLEN); | |
| while (! feof(stdin)) { | |
| memset(sBuff, 0, BUFFSIZE); | |
| fgets(sBuff, BUFFSIZE - 1, stdin); | |
| memset(sTempWord, 0, TOKENLEN); | |
| memset(sTempDoc, 0, TOKENLEN); | |
| if (sscanf(sBuff, "W %s %s %d\n", sTempDoc, sTempWord, &nTF) == 3) { | |
| if (strcmp(sPrevDoc, sTempDoc) != 0) { | |
| nNumOfDocs ++; | |
| nDocHeadCell = (int *) realloc(nDocHeadCell, sizeof(int) * nNumOfDocs); | |
| nDocNumOfCells = (int *) realloc(nDocNumOfCells, sizeof(int) * nNumOfDocs); | |
| nDocLen = (int *) realloc(nDocLen, sizeof(int) * nNumOfDocs); | |
| sDocName = (char **) realloc(sDocName, sizeof(char *) * nNumOfDocs); | |
| nDocHeadCell[nNumOfDocs - 1] = nNumOfCells; | |
| nDocNumOfCells[nNumOfDocs - 1] = 0; | |
| nDocLen[nNumOfDocs - 1] = 0; | |
| sDocName[nNumOfDocs - 1] = (char *) malloc(sizeof(char) * (strlen(sTempDoc) + 1)); | |
| memset(sDocName[nNumOfDocs - 1], 0, strlen(sTempDoc) + 1); | |
| strcpy(sDocName[nNumOfDocs - 1], sTempDoc); | |
| } | |
| nDocID = nNumOfDocs - 1; | |
| nWordID = addWord(sTempWord); | |
| if (nWordID > 0) { | |
| nNumOfCells ++; | |
| nDocNumOfCells[nDocID] ++; | |
| nWordCell = (int *) realloc(nWordCell, sizeof(int) * nNumOfCells); | |
| nWordCell[nNumOfCells - 1] = nWordID - 1; | |
| nTFCell = (int *) realloc(nTFCell, sizeof(int) * nNumOfCells); | |
| nTFCell[nNumOfCells - 1] = nTF; | |
| nDocLen[nDocID] += nTF; | |
| nNumOfTokens += nTF; | |
| } else { fprintf(stderr, "Invalid data : %s\n", sTempWord); exit(1); } | |
| memset(sPrevDoc, 0, TOKENLEN); | |
| strcpy(sPrevDoc, sTempDoc); | |
| } | |
| } | |
| if (nNumOfDocs == 0) { fprintf(stderr, "no valid data.\n"); exit(1); } | |
| { | |
| int nCD; | |
| nMaxDocNumOfCells = 0; | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) | |
| if (nDocNumOfCells[nCD] > nMaxDocNumOfCells) nMaxDocNumOfCells = nDocNumOfCells[nCD]; | |
| } | |
| setTestData(); | |
| // print statistics | |
| { | |
| fprintf(fOut, "### J : %d; ", nNumOfDocs); | |
| fprintf(fOut, "W : %d; ", nNumOfWords); | |
| fprintf(fOut, "K : %d; ", nNumOfTopics); | |
| fprintf(fOut, "JW pairs : %d; ", nNumOfCells); | |
| fprintf(fOut, "M : %d + %d; ", nNumOfTokens, nNumOfTestTokens); | |
| fprintf(fOut, "X : %d; ", nMaxDocNumOfCells); | |
| fprintf(fOut, "\n"); | |
| fflush(fOut); | |
| } | |
| // initialization | |
| { | |
| int nCD, nCC, nCToken; | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) | |
| for (nCC = nDocHeadCell[nCD]; nCC < nDocHeadCell[nCD] + nDocNumOfCells[nCD]; nCC ++) | |
| for (nCToken = 0; nCToken < nTFCell[nCC]; nCToken ++) | |
| nTopicArray[nCC][nCToken] = genrand_int32(genrand_int32(nCD % NUM_OF_THREADS) % NUM_OF_THREADS) % nNumOfTopics; | |
| } | |
| dDigammaOne = digamma(1.0); | |
| { | |
| int nCK, nCW; | |
| dAlpha = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfTopics); | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dAlpha[nCK] = 50.0 / ((DTYPE) nNumOfTopics); | |
| dAlphaSum = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dAlphaSum += dAlpha[nCK]; | |
| dBeta = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfWords); | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) dBeta[nCW] = 0.01; | |
| dBetaSum = 0.0; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) dBetaSum += dBeta[nCW]; | |
| dLambda = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfDocs * nNumOfTopics); | |
| dLambdaSum = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfDocs); | |
| dZeta = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfTopics * nNumOfWords); | |
| dZetaSum = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfTopics); | |
| dZeta2 = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfTopics * nNumOfWords); | |
| dZetaSum2 = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfTopics); | |
| dNu = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfDocs); | |
| dMeanR = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfDocs * nNumOfTopics); | |
| dStdevS = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfDocs * nNumOfTopics); | |
| dMu = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfTopics); | |
| dInvSigma = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfTopics * nNumOfTopics); | |
| } | |
| statistics(); | |
| // initialization by CGS for LDA | |
| nIteration = 0; | |
| time(&start_time); | |
| do { | |
| #pragma omp parallel | |
| { | |
| DTYPE *dTempDist; | |
| int nCD, nCC, nCToken, nCK, nCW; | |
| int nOldTopic, nNewTopic; | |
| DTYPE dTemp; | |
| #ifdef _OPENMP | |
| int nThreadID = omp_get_thread_num(); | |
| int nNumOfThreads = omp_get_num_threads(); | |
| #else | |
| int nThreadID = 0; | |
| int nNumOfThreads = 1; | |
| #endif | |
| dTempDist = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfTopics); | |
| for (nCD = nThreadID; nCD < nNumOfDocs; nCD += nNumOfThreads) { | |
| for (nCC = nDocHeadCell[nCD]; nCC < nDocHeadCell[nCD] + nDocNumOfCells[nCD]; nCC ++) { | |
| nCW = nWordCell[nCC]; | |
| for (nCToken = 0; nCToken < nTFCell[nCC]; nCToken ++) { | |
| nOldTopic = nTopicArray[nCC][nCToken]; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| if (nCK == nOldTopic) | |
| dTempDist[nCK] = (dLambda[nCD * nNumOfTopics + nCK] + dAlpha[nCK] - 1.0) | |
| * (dZeta[nCK * nNumOfWords + nCW] + dBeta[nCW] - 1.0) / (dZetaSum[nCK] + dBetaSum - 1.0); | |
| else | |
| dTempDist[nCK] = (dLambda[nCD * nNumOfTopics + nCK] + dAlpha[nCK]) | |
| * (dZeta[nCK * nNumOfWords + nCW] + dBeta[nCW]) / (dZetaSum[nCK] + dBetaSum); | |
| if (dTempDist[nCK] < 0.0) dTempDist[nCK] = 0.0; | |
| } | |
| dTemp = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dTemp += dTempDist[nCK]; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dTempDist[nCK] /= dTemp; | |
| dTemp = genrand_res53(nThreadID); | |
| nCK = 0; | |
| while ((dTemp > 0.0) && (nCK < nNumOfTopics)) dTemp -= dTempDist[nCK ++]; | |
| nNewTopic = nCK - 1; | |
| if (nNewTopic == -1) nNewTopic = 0; | |
| if (nOldTopic != nNewTopic) { | |
| nTopicArray[nCC][nCToken] = nNewTopic; | |
| dLambda[nCD * nNumOfTopics + nOldTopic] -= 1.0; | |
| dLambda[nCD * nNumOfTopics + nNewTopic] += 1.0; | |
| #pragma omp critical | |
| { | |
| dZeta[nOldTopic * nNumOfWords + nCW] -= 1.0; | |
| dZeta[nNewTopic * nNumOfWords + nCW] += 1.0; | |
| dZetaSum[nOldTopic] -= 1.0; | |
| dZetaSum[nNewTopic] += 1.0; | |
| } | |
| } | |
| } | |
| } | |
| } | |
| free(dTempDist); | |
| } | |
| nIteration ++; | |
| fprintf(fOut, "#i\t%d\t", nIteration); | |
| fprintf(fOut, "\t%.3lf\t", perplexityLDA(0)); | |
| put_time(fOut); | |
| fprintf(fOut, "\n"); | |
| fflush(fOut); | |
| if ((nIteration > 10) | |
| && ((nIteration < 30) | |
| || ((nIteration < 100) && (nIteration % 10 == 0)) | |
| || (nIteration % 20 == 0))) { | |
| optHypA2(fOut, nIteration); | |
| optHypB2(fOut, nIteration); | |
| } | |
| } while (nIteration < NUM_OF_CGS_ITERATIONS); | |
| // | |
| { | |
| int nCD, nCK, nCW; | |
| int nTemp; | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) { | |
| nTemp = nCD * nNumOfTopics; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dMeanR[nTemp + nCK] = log(dLambda[nTemp + nCK] + dAlpha[nCK]); | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dStdevS[nTemp + nCK] = 1.0; | |
| } | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| nTemp = nCK * nNumOfWords; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) dZeta[nTemp + nCW] += dBeta[nCW]; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) dZeta2[nTemp + nCW] += digamma(dZeta[nTemp + nCW]); | |
| dZetaSum[nCK] = 0.0; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) dZetaSum[nCK] += dZeta[nTemp + nCW]; | |
| dZetaSum2[nCK] = digamma(dZetaSum[nCK]); | |
| } | |
| } | |
| // VB for LOCCTM | |
| nIteration = 0; | |
| do { | |
| #pragma omp parallel | |
| { | |
| int nCD, nCK, nCK2, nCW, nCC; | |
| DTYPE dTemp; | |
| int nTemp; | |
| DTYPE *dPi; | |
| lbfgsfloatval_t dFX; | |
| lbfgsfloatval_t *dX; | |
| lbfgs_parameter_t param; | |
| #ifdef _OPENMP | |
| int nThreadID = omp_get_thread_num(); | |
| int nNumOfThreads = omp_get_num_threads(); | |
| #else | |
| int nThreadID = 0; | |
| int nNumOfThreads = 1; | |
| #endif | |
| dPi = (DTYPE *) malloc(sizeof(DTYPE) * nNumOfTopics); | |
| dX = lbfgs_malloc(nNumOfDocs * 2); | |
| #pragma omp barrier | |
| for (nCK = nThreadID; nCK < nNumOfTopics; nCK += nNumOfThreads) { | |
| dMu[nCK] = 0.0; | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) | |
| dMu[nCK] += dMeanR[nCD * nNumOfTopics + nCK]; | |
| dMu[nCK] /= ((DTYPE) nNumOfDocs); | |
| } | |
| #pragma omp barrier | |
| for (nCK = nThreadID; nCK < nNumOfTopics; nCK += nNumOfThreads) { | |
| dPi[nCK] = 0.0; | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) { | |
| dTemp = dStdevS[nCD * nNumOfTopics + nCK]; | |
| dPi[nCK] += dTemp * dTemp; | |
| } | |
| dPi[nCK] /= ((DTYPE) nNumOfDocs); | |
| nTemp = nCK * nNumOfTopics; | |
| for (nCK2 = nCK; nCK2 < nNumOfTopics; nCK2 ++) { | |
| dInvSigma[nTemp + nCK2] = 0.0; | |
| for (nCD = 0; nCD < nNumOfDocs; nCD ++) | |
| dInvSigma[nTemp + nCK2] += (dMeanR[nCD * nNumOfTopics + nCK] - dMu[nCK]) | |
| * (dMeanR[nCD * nNumOfTopics + nCK2] - dMu[nCK2]); | |
| dInvSigma[nTemp + nCK2] /= ((DTYPE) nNumOfDocs); | |
| if (nCK != nCK2) | |
| dInvSigma[nCK2 * nNumOfTopics + nCK] = dInvSigma[nTemp + nCK2]; | |
| } | |
| nTemp = nCK * nNumOfWords; | |
| dZetaSum[nCK] = dBetaSum; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) dZeta[nTemp + nCW] = dBeta[nCW]; | |
| } | |
| #pragma omp barrier | |
| if (nThreadID == 0) { | |
| dDiagonal = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dDiagonal += dInvSigma[nCK * nNumOfTopics + nCK]; | |
| dDiagonal /= ((DTYPE) nNumOfTopics); | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dInvSigma[nCK * nNumOfTopics + nCK] | |
| += (1.0 - DIAGONAL_RATIO) * dPi[nCK] + DIAGONAL_RATIO * dDiagonal; | |
| dLogDetSigma = inverse(dInvSigma, nNumOfTopics); | |
| } | |
| #pragma omp barrier | |
| for (nCD = nThreadID; nCD < nNumOfDocs; nCD += nNumOfThreads) { | |
| nTemp = nCD * nNumOfTopics; | |
| dNu[nCD] = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| dNu[nCD] += exp(dMeanR[nTemp + nCK] | |
| + 0.5 * dStdevS[nTemp + nCK] * dStdevS[nTemp + nCK]); | |
| } | |
| dLambdaSum[nCD] = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dLambda[nTemp + nCK] = 0.0; | |
| for (nCC = nDocHeadCell[nCD]; nCC < nDocHeadCell[nCD] + nDocNumOfCells[nCD]; nCC ++) { | |
| nCW = nWordCell[nCC]; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dPi[nCK] = dMeanR[nTemp + nCK] | |
| + dZeta2[nCK * nNumOfWords + nCW] - dZetaSum2[nCK]; | |
| dTemp = dPi[0]; | |
| for (nCK = 1; nCK < nNumOfTopics; nCK ++) if (dTemp < dPi[nCK]) dTemp = dPi[nCK]; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dPi[nCK] = exp(dPi[nCK] - dTemp); | |
| dTemp = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) dTemp += dPi[nCK]; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dPi[nCK] = ((DTYPE) nTFCell[nCC]) * dPi[nCK] / dTemp; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| dLambda[nTemp + nCK] += dPi[nCK]; | |
| dLambdaSum[nCD] += dPi[nCK]; | |
| } | |
| #pragma omp critical | |
| { | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| dZeta[nCK * nNumOfWords + nCW] += dPi[nCK]; | |
| dZetaSum[nCK] += dPi[nCK]; | |
| } | |
| } | |
| } | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| dX[nCK] = dMeanR[nTemp + nCK]; | |
| dX[nNumOfTopics + nCK] = dStdevS[nTemp + nCK]; | |
| } | |
| lbfgs_parameter_init(¶m); | |
| param.linesearch = LBFGS_LINESEARCH_BACKTRACKING_ARMIJO; | |
| lbfgs(nNumOfTopics * 2, dX, &dFX, LBFGS_MeanR, LBFGS_progress, &nCD, ¶m); | |
| //printf("LBFGS(%d)\t%d\t", nRet, nCD); | |
| //printf("fx = %.6le, x[0] = %.6le, x[1] = %.6le\n", dFX, dX[nNumOfTopics], dX[nNumOfTopics + 1]); | |
| if ((! isnan(dFX)) && (! isinf(dFX))) { | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) { | |
| dMeanR[nTemp + nCK] = dX[nCK]; | |
| if (dX[nNumOfTopics + nCK] > 0.0) | |
| dStdevS[nTemp + nCK] = dX[nNumOfTopics + nCK]; | |
| } | |
| nTemp = nCD * nNumOfTopics; | |
| dNu[nCD] = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dNu[nCD] += exp(dMeanR[nTemp + nCK] | |
| + 0.5 * dStdevS[nTemp + nCK] * dStdevS[nTemp + nCK]); | |
| } | |
| } | |
| #pragma omp barrier | |
| for (nCK = nThreadID; nCK < nNumOfTopics; nCK += nNumOfThreads) { | |
| dZetaSum2[nCK] = digamma(dZetaSum[nCK]); | |
| nTemp = nCK * nNumOfWords; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) | |
| dZeta2[nTemp + nCW] = digamma(dZeta[nTemp + nCW]); | |
| } | |
| #pragma omp barrier | |
| for (nCW = nThreadID; nCW < nNumOfWords; nCW += nNumOfThreads) { | |
| dTemp = dBeta[nCW]; | |
| dBeta[nCW] = 0.0; | |
| for (nCK = 0; nCK < nNumOfTopics; nCK ++) | |
| dBeta[nCW] += (dZeta2[nCK * nNumOfWords + nCW] - dZetaSum2[nCK]); | |
| dBeta[nCW] = invDigamma(digamma(dBetaSum) + dBeta[nCW] / ((DTYPE) nNumOfTopics)); | |
| if ((dBeta[nCW] < MIN_HYP) || isnan(dBeta[nCW]) || isinf(dBeta[nCW])) dBeta[nCW] = dTemp; | |
| } | |
| #pragma omp barrier | |
| free(dPi); | |
| lbfgs_free(dX); | |
| } | |
| { | |
| int nCW; | |
| dBetaSum = 0.0; | |
| for (nCW = 0; nCW < nNumOfWords; nCW ++) dBetaSum += dBeta[nCW]; | |
| } | |
| fprintf(fOut, "#I\t%d\t", nIteration); | |
| fprintf(fOut, "%.3lf\t", logEviLB()); | |
| fprintf(fOut, "%.3lf\t", perplexityLDA(1)); | |
| put_time(fOut); | |
| fprintf(fOut, "\t%d\t%.6lf", nIteration, dBetaSum / ((DTYPE) nNumOfWords)); | |
| fprintf(fOut, "\n"); | |
| fflush(fOut); | |
| nIteration ++; | |
| if (nIteration == 1) printLog(nIteration); | |
| } while (nIteration < NUM_OF_VB_ITERATIONS); | |
| printLog(nIteration); | |
| exit(0); | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment