hfb.c

实现HMM算法

         }
         if (trace&T_TOP) {
            printf("Retrying Beta pass at %5.1f\n",pruneThresh);
         }
      }
   while(pruneThresh<=pruneSetting.pruneLim);
   
   if (lbeta<LSMALL)
      HError(7323,"StepBack: Beta prune error");
   return TRUE;
}

/* ---------------------- Statistics Accumulation -------------------- */

/* UpTranParms: update the transition counters of given hmm */
static void UpTranParms(FBInfo *fbInfo, HLink hmm, int t, int q,
                        DVector aqt, DVector bqt, DVector bqt1, DVector bq1t, 
                        LogDouble pr)
{
   int i,j,N;
   Vector ti,ai;
   float ***outprob,***outprob1;
   double sum,x;
   TrAcc *ta;
   AlphaBeta *ab;

   N = hmm->numStates;
   ab = fbInfo->ab;
   ta = (TrAcc *) GetHook(hmm->transP);/*转移矩阵Accumulators*/
   outprob = ab->otprob[t][q]; /*t时刻第q个模型的观察概率*/
   if (bqt1!=NULL) outprob1 = ab->otprob[t+1][q];  /* Bug fix */
   else outprob1 = NULL;
   for (i=1;i<N;i++)
      ta->occ[i] += ab->occt[i];/* 状态占有累加器，每个状态i的状态占有等于每个时刻t该状态i的状态占有的累加state occupation */
   for (i=1;i<N;i++) {
      ti = ta->tran[i]; ai = hmm->transP[i];/*ti为第i个状态的转移向量的累加器，ai为当前时刻t第i个状态的转移向量*/
      for (j=2;j<=N;j++) {
         if (i==1 && j<N) {                  /* 入口状态1到其它状态j(j<N)的转移entry transition */
            x = aqt[1]+ai[j]+outprob[j][0][0]+bqt[j]-pr;/*计算HTKBOOK第131页的公式(1)分子部分的时间t的累加单项，此为对数概率形式*/
            if (x>MINEARG) ti[j] += exp(x);/*状态i和状态j之间转移的时间t的累加单项，恢复为原始概率*/
         } else
            if (i>1 && j<N && bqt1!=NULL) {     /* 内部状态之间的转移internal transition */
               x = aqt[i]+ai[j]+outprob1[j][0][0]+bqt1[j]-pr;/*计算HTKBOOK第130页的公式(1)分子部分的时间t的累加单项，此为对数概率形式*/
               if (x>MINEARG) ti[j] += exp(x);/*状态i和状态j之间转移的时间t的累加单项，恢复为原始概率*/
            } else
               if (i>1 && j==N) {/*计算HTKBOOK第131页的公式(2)分子部分的时间t的累加单项，此为对数概率形式*/ /* 状态i(i>1)到模型出口状态Nq的转移exit transition */
                  x = aqt[i]+ai[N]+bqt[N]-pr;
                  if (x>MINEARG) ti[N] += exp(x);/*状态i和状态N之间转移的时间t的累加单项，恢复为原始概率*/
               }
         if (i==1 && j==N && ai[N]>LSMALL && bq1t != NULL){ /*计算HTKBOOK第131页的公式(3)分子部分的时间t的累加单项，此为对数概率形式*//* 模型的入口状态1到出口状态Nq的转移tee transition */
            x = aqt[1]+ai[N]+bq1t[1]-pr;
            if (x>MINEARG) ti[N] += exp(x);/*状态1和状态N之间转移的时间t的累加单项，恢复为原始概率*/
         }
      }
   }
   if (trace&T_TRA && NonSkipRegion(fbInfo->skipstart,fbInfo->skipend,t)) {/*当前不执行此处*/
      printf("Tran Counts at time %d, Model Q%d %s\n",t,q,ab->qIds[q]->name);
      for (i=1;i<=N;i++) {
         printf("  %d. Occ %8.2f: Trans ",i,ta->occ[i]);
         sum = 0.0;
         for (j=2; j<=N; j++) {
            x = ta->tran[i][j]; sum += x;
            printf("%7.2f ",x);
         }
         printf(" [%8.2f]\n",sum);
      }
   }
}

/* UpMixParms: update mu/va accs of given hmm在当前时刻t，更新第q个模型的相关累加器 */
static void UpMixParms(FBInfo *fbInfo, int q, HLink hmm, HLink al_hmm,
                       Observation ot, Observation ot2, 
                       int t, DVector aqt, DVector aqt1, DVector bqt, int S,
                       Boolean twoDataFiles, LogDouble pr)
{
   int i,s,j,k,kk,m=0,mx,M=0,N,vSize;
   Vector mu_jm,var,mean=NULL,invk,otvs;
   TMixRec *tmRec = NULL;
   float **outprob;
   Matrix inv;
   LogFloat c_jm,a,prob=0.0;
   LogDouble x,initx = LZERO;
   float zmean,zmeanlr,zmean2,tmp;
   double Lr,steSumLr;
   HMMSet *hset;
   HSetKind hsKind;
   AlphaBeta *ab;
   StreamElem *ste, *al_ste=NULL;
   MixtureElem *me;
   MixPDF *mp=NULL;
   MuAcc *ma;
   VaAcc *va;
   WtAcc *wa = NULL;
   Boolean mmix=FALSE;  /* TRUE if multiple mixture */
   float wght=0.0;
   /* variables for 2-model reestimation */
   Vector comp_prob=NULL;        /* array[1..M] of Component probability */
   float norm=0.0;               /* total mixture prob */
   LogFloat det;
   AdaptXForm *inxform;

   ab     = fbInfo->ab;
   hset   = fbInfo->up_hset;
   hsKind = fbInfo->hsKind;

   if (trace&T_MIX && fbInfo->uFlags&UPMIXES && 
       NonSkipRegion(fbInfo->skipstart,fbInfo->skipend,t)){/*当前不执行此处*/
      printf("Mixture Weights at time %d, model Q%d %s\n",
             t,q,ab->qIds[q]->name);
   }

   if (fbInfo->twoModels)/*当前不执行此处*/
       comp_prob = CreateVector(&gstack,fbInfo->maxM);

   N = hmm->numStates;
   for (j=2;j<N;j++) {/*第1个状态和最后一个状态都没有观察*/
      if (fbInfo->maxM>1){/*当前此判断条件为假，不执行此处*/
         initx = hmm->transP[1][j] + aqt[1];
         if (t>1)
            for (i=2;i<N;i++){
               a = hmm->transP[i][j];
               if (a>LSMALL)
                  initx = LAdd(initx,aqt1[i]+a);
            }
         initx += bqt[j] - pr;
      }/*以上不被执行*/
      if (trace&T_MIX && fbInfo->uFlags&UPMIXES && 
          NonSkipRegion(fbInfo->skipstart,fbInfo->skipend,t))/*当前此判断条件为假，不执行此处*/
         printf("  State %d: ",j);/*当前从下面开始执行*/
      ste = hmm->svec[j].info->pdf+1;/*当前要更新模型的第j个状态的流信息指针*/
      outprob = ab->otprob[t][q][j];/*当前时刻t第q个模型的第j个状态的观察概率*/
      for (s=1;s<=S;s++,ste++){
         /* Get observation vector for this state/stream */
         vSize = hset->swidth[s];/*当前特征向量的长度vSize=39*/
      
         switch (hsKind){/*当前ksKind=1,即等于SHAREDHS*/
         case TIEDHS:             /* if tied mixtures then we only */
            tmRec = &(hset->tmRecs[s]); /* want to process the non-pruned */
            M = tmRec->topM;            /* components */
            mmix = TRUE;
            break;
         case DISCRETEHS:
            M = 1;
            mmix = FALSE;
            break;
         case PLAINHS:
         case SHAREDHS:/*当前执行此分支*/
            M = ste->nMix;/*当前流的混合数*/
            mmix = (M>1);
            break;
         }
         /* update weight occupation count 更新权占有计数*/
         wa = (WtAcc *) ste->hook; steSumLr = 0.0;/*状态元素Lr加和*/

         if (fbInfo->twoModels) { /* component probs of update hmm *//*当前不执行这里*/
            norm = LZERO;
            for (mx=1; mx<=M; mx++) {
               if (alCompLevel) {
                  al_ste = al_hmm->svec[j].info->pdf+1;
                  if (al_ste->nMix != M)
                     HError(999,"Cannot align at the component level if number of components is different!");
               }
               if (alCompLevel) {
                  me = al_ste->spdf.cpdf+mx;
               } else {
                  me = ste->spdf.cpdf+mx;
               }
	           inxform = fbInfo->inXForm;
               mp=me->mpdf;
               if (twoDataFiles){  
                  if (alCompLevel) {
                     otvs = ApplyCompFXForm(mp,ot.fv[s],fbInfo->al_inXForm,&det,t);
                  }
                  else {
                     otvs = ApplyCompFXForm(mp,ot2.fv[s],inxform,&det,t);
                  }
               } else {
                  otvs = ApplyCompFXForm(mp,ot.fv[s],inxform,&det,t);
               }
               wght = MixLogWeight(hset,me->weight);
               comp_prob[mx]=wght+MOutP(otvs,mp)+det;
               norm = LAdd(norm,comp_prob[mx]);
            }
         }/*从下面开始执行*/

         for (mx=1;mx<=M;mx++) { 
            /* process mixtures */
            switch (hsKind){    /* Get wght and mpdf */
            case TIEDHS:
               m=tmRec->probs[mx].index;
               wght=MixLogWeight(hset,ste->spdf.tpdf[m]);
               mp=tmRec->mixes[m];
               break;
            case DISCRETEHS:
               if (twoDataFiles)
                  m=ot2.vq[s];
               else
                  m=ot.vq[s];
               wght = 0.0; /* This is the log-weight, just for consistency! */
               mp=NULL;
               break;
            case PLAINHS:
            case SHAREDHS:/*当前执行这里*/
               m = mx;/*当前混合索引*/
               me = ste->spdf.cpdf+m;/*当前MixtureElem */
               wght = MixLogWeight(hset,me->weight);/*混合权的对数值*/
               mp=me->mpdf;/* 当前混合元素的 mixture pdf */
               break;
            }
            if (wght>LMINMIX){/*当前执行这里*/
              /* compute mixture likelihood  */
              if (!mmix || (hsKind==DISCRETEHS))/*当前mmix=F，所以此判断条件为真*/ /* For DISCRETEHS calcs are*/
                 x = aqt[j]+bqt[j]-pr;/* same as single mix*//*计算了HTKBOOK第8页的公式(1.26)*/
              else if (fbInfo->twoModels) { /*下面部分不被执行*/     /* note: only SHAREDHS or PLAINHS */
                 x = comp_prob[m]+aqt[j]+bqt[j]-pr-norm;
              }
              else {
                 c_jm=wght;
                 x = initx+c_jm;
                 switch(hsKind) {
                 case TIEDHS :
                    tmp = tmRec->probs[mx].prob;
                    prob = (tmp>=MINLARG)?log(tmp)+tmRec->maxP:LZERO;
                    break;
                 case PLAINHS : 
                 case SHAREDHS:
		            if (S==1)
		              prob = outprob[0][mx];
		            else
		              prob = outprob[s][mx];
                    break;
                 default:
                    x=LZERO;
                    break;
                 }
                 x += prob;
                 if (S>1)      /* adjust for parallel streams */
                    x += outprob[s][0];
              }/*以上部分不被执行， c_jm=wght;之上的else分支结束*/
              if (twoDataFiles){ /*当前此判断条件为假*/ /* switch to new data for mu & var est */
                 otvs = ot2.fv[s];
              }
              if (-x<pruneSetting.minFrwdP) {
		         if (twoDataFiles){  /* switch to new data for mu & var est */
		            otvs = ApplyCompFXForm(mp,ot2.fv[s],fbInfo->paXForm,&det,t);
				 } else {/*执行这里*/
		         otvs = ApplyCompFXForm(mp,ot.fv[s],fbInfo->paXForm,&det,t);/*det=0; return fv[s]，otvs为从调用参数传过来的观察的特征向量*/
				 }
                 Lr = exp(x);/*将对数概率转换成了原始概率，计算了HTKBOOK第8页的公式(1.26)*/
                  /* More diagnostics */
                  /* if (Lr>0.000001 && ab->occt[j]>0.000001 &&
                     (Lr/ab->occt[j])>1.00001)
                     printf("Too big %d %d %s : %5.3f %10.2f %8.2f (%4.2f)\n",t,q,
                     ab->qIds[q]->name,Lr/ab->occt[j],Lr,ab->occt[j],prob); */
            
                  /* update occupation counts */
                  steSumLr += Lr;
                  /* update the adaptation statistic counts */
                  if (fbInfo->uFlags&UPXFORM)/*当前不执行此处*/
                     AccAdaptFrame(Lr, otvs, mp, t);
                  /* update mean counts */
                  if ((fbInfo->uFlags&UPMEANS) || (fbInfo->uFlags&UPVARS))/*当前执行此处*/
                     mean = mp->mean; 
                  if ((fbInfo->uFlags&UPMEANS) && (fbInfo->uFlags&UPVARS)) {/*当前执行此处*/
                     ma = (MuAcc *) GetHook(mean);
                     va = (VaAcc *) GetHook(mp->cov.var);
                     ma->occ += Lr;/* HTKbook第7页第(1.13)公式的分母累加单项，共享这个混合密度函数的状态占有数occ for states sharing this mpdf */
                     va->occ += Lr; /* HTKbook第7页第(1.14)公式的分母累加单项，共享这个混合密度函数的状态占有数occ for states sharing this mpdf */
                     mu_jm = ma->mu;/* mean vector counts */
                     if ((mp->ckind==DIAGC)||(mp->ckind==INVDIAGC)){/*当前执行此处*/
                       /*
						ste为从调用参数中传递过来的HMM的第j个状态元素的状态信息的流元素
                        otvs为从调用参数传过来的观察的特征向量
                        mp为为从调用参数中传递过来的HMM的第j个状态元素的状态信息的流元素的第m个混合的pdf 
						 用新的观察更新方差
						*/ 
						var = va->cov.var;
                        for (k=1;k<=vSize;k++) {/*更新方差，计算HTKBOOK第7页公式(1.14)中分子的求和单项*/
                           zmean=otvs[k]-mean[k];/*当前观察-均值*/
                           zmeanlr=zmean*Lr;/*将上面的差乘以权值L*//* 当前时刻t，当前状态j的均值累加，HTKbook第7页第(1.13)公式的分子累加单项*/
                           mu_jm[k] += zmeanlr;
                           var[k] += zmean*zmeanlr;;/*计算HTKBOOK第7页公式(1.14)中分子的求和单项，即当前时刻t，当前状态j的方差累加*/
                        }/*k循环结束*/
                     } else {/*当前不执行此分之*/
                        inv = va->cov.inv;
                        for (k=1;k<=vSize;k++) {
                           invk = inv[k];
                           zmean=otvs[k]-mean[k];
                           zmeanlr=zmean*Lr;
                           mu_jm[k] += zmeanlr;
                           for (kk=1;kk<=k;kk++) {
                              zmean2 = otvs[kk]-mean[kk];
                              invk[kk] += zmean2*zmeanlr;
                           }/*kk循环结束*/
                        }/*k循环结束*/
                     }/*inv = va->cov.inv;之上的else分支结束*/
                  }/*ma = (MuAcc *) GetHook(mean);之上的if分支结束*/
                  else if (fbInfo->uFlags&UPMEANS){/*当前不执行这里*/
                     ma = (MuAcc *) GetHook(mean);
                     mu_jm = ma->mu;
                     ma->occ += Lr;
                     for (k=1;k<=vSize;k++)     /* sum zero mean */
                        mu_jm[k] += (otvs[k]-mean[k])*Lr;
                  }
                  else if (fbInfo->uFlags&UPVARS){/*当前不执行这里*/
                     /* update covariance counts */
                     va = (VaAcc *) GetHook(mp->cov.var);
                     va->occ += Lr;
                     if ((mp->ckind==DIAGC)||(mp->ckind==INVDIAGC)){
                        var = va->cov.var;
                        for (k=1;k<=vSize;k++) {
                           zmean=otvs[k]-mean[k];
                           var[k] += zmean*zmean*Lr;
                        }
                     } else {
                        inv = va->cov.inv;
                        for (k=1;k<=vSize;k++) {
                           invk = inv[k];
                           zmean=otvs[k]-mean[k];               
                           for (kk=1;kk<=k;kk++) {
                              zmean2 = otvs[kk]-mean[kk];
                              invk[kk] += zmean*zmean2*Lr;
                           }
                        }
                     }
                  }
            
                  /* update mixture weight counts */
                  if (fbInfo->uFlags&UPMIXES) {/*当前执行这里*/
                     wa->c[m] +=Lr;
                     if (trace&T_MIX && NonSkipRegion(fbInfo->skipstart,fbInfo->skipend,t))
                        printf("%3d. %7.2f",m,wa->c[m]);
                  }
               }
            }
            if (twoDataFiles){ /* Switch back to old data for prob calc */
               otvs = ot.fv[s]