hfb.c

实现HMM算法

}

/* StepAlpha: calculate alphat column for time t and return
   forward beam limits in startq and endq */
static void StepAlpha(AlphaBeta *ab, int t, int *start, int *end, 
                      int Q, int T, LogDouble pr, int skipstart, int skipend)
{
   DVector aq,laq,*tmp, *alphat,*alphat1;
   PruneInfo *p;
   float ***outprob;
   int sq,eq,i,j,q,Nq,lNq;
   LogDouble x=0.0,y,a,a1N=0.0;
   HLink hmm;
   
   alphat  = ab->alphat;
   alphat1 = ab->alphat1;

   /* First prune beta beam further to get alpha beam */
   p = ab->pInfo;
   sq = p->qLo[t-1];    /* start start-point at bottom of beta beam at t-1 */

   while (pr-MaxModelProb(ab,sq,t-1,sq)>pruneSetting.minFrwdP){
      ++sq;                /* raise start point */
      if (sq>p->qHi[t]) 
         HError(7390,"StepAlpha: Alpha prune failed sq(%d) > qHi(%d)",sq,p->qHi[t]);
   }
   if (sq<p->qLo[t])       /* start-point below beta beam so pull it back */
      sq = p->qLo[t];
   
   eq = p->qHi[t-1]<Q?p->qHi[t-1]+1:p->qHi[t-1];
   /* start end-point at top of beta beam at t-1  */
   /* JJO : + 1 to allow for state q-1[N] -> q[1] */
   /*       + 1 for each tee model following eq.  */
   while (pr-MaxModelProb(ab,eq,t-1,sq)>pruneSetting.minFrwdP){
      --eq;             /* lower end-point */
      if (eq<sq) 
         HError(7390,"StepAlpha: Alpha prune failed eq(%d) < sq(%d)",eq,sq);
   }
   while (eq<Q && ab->qDms[eq]==0) eq++;
   if (eq>p->qHi[t])  /* end point above beta beam so pull it back */
      eq = p->qHi[t]; 
      
   if (trace&T_PRU && NonSkipRegion(skipstart,skipend,t)){
      printf("%d: Alpha Beam %d->%d \n",t,sq,eq);
      fflush(stdout);
   } 
   
   /* Now compute current alpha column */
   tmp = ab->alphat1; ab->alphat1 = ab->alphat; ab->alphat = tmp;
   alphat  = ab->alphat;
   alphat1 = ab->alphat1;
   
   if (sq>1) ZeroAlpha(ab,1,sq-1);
   
   Nq = (sq == 1) ? 0:ab->al_qList[sq-1]->numStates;/*此时Nq中保存的是前一个模型的状态数*/
   
   for (q = sq; q <= eq; q++) {
      lNq = Nq; hmm = ab->al_qList[q]; Nq = hmm->numStates; /*此时lNq中保存的是前一个模型的状态数*/
      aq = alphat[q];/*第q个模型当前时刻t的alpha向量*/ 
      laq = alphat1[q];/*第q个模型t-1时刻的alpha向量*/ 
      if (laq == NULL)
         HError(7322,"StepAlpha: laq gone wrong!");
      if((outprob = ab->otprob[t][q]) == NULL)
         HError(7322,"StepAlpha: Outprob NULL at time %d model %d in StepAlpha",t,q);
      if (q==1)/*计算HTKBOOK第129页Alpha计算公式的第(1)个式子的第一部分*/
         aq[1] = LZERO;
      else{/*计算HTKBOOK第129页Alpha计算公式的第(1)个式子的第二部分*/
         aq[1] = alphat1[q-1][lNq];
         if (q>sq && a1N>LSMALL) /* tee Model */
            aq[1] = LAdd(aq[1],alphat[q-1][1]+a1N);
      }
      for (j=2;j<Nq;j++) {/*计算HTKBOOK第129页Alpha计算公式的第(2)个式子*/
         a = hmm->transP[1][j];
         x = (a>LSMALL)?a+aq[1]:LZERO;/*Alpha计算公式的第(2)个式子的第一个求和公式*/
         for (i=2;i<Nq;i++){/*Alpha计算公式的第(2)个式子的第二个求和公式*/
            a = hmm->transP[i][j]; y = laq[i];
            if (a>LSMALL && y>LSMALL)
               x = LAdd(x,y+a);
         }
         aq[j] = x + outprob[j][0][0];/*Alpha计算公式的第(2)个式子的乘积部分*/
      }
      x = LZERO;/*下面计算HTKBOOK第129页Alpha计算公式的第(3)个式子*/
      for (i=2;i<Nq;i++){
         a = hmm->transP[i][Nq]; y = aq[i];
         if (a>LSMALL && y>LSMALL)
            x = LAdd(x,y+a);
      }
      aq[Nq] = x; a1N = hmm->transP[1][Nq];/*第(3)个式子计算完毕*/
   }/*q循环结束*/
   if (eq<Q) ZeroAlpha(ab,eq+1,Q);

   if (trace&T_PRU && p->pruneThresh < NOPRUNE)
	  CheckPruning(ab,t,skipstart,skipend);
   if (t==T){
      if (fabs((x-pr)/T) > 0.001)/*当前不执行此处*/
         HError(-7391,"StepAlpha: Forward/Backward Disagree %f/%f",x,pr);
      if (trace&T_PRU && p->pruneThresh < NOPRUNE) /*当前不执行此处*/
	     SummarisePruning(p, Q, T);	  
   }

   *start=sq; *end=eq;

}


/* CreateBeta: create Q and T pointer arrays for beta */
static void CreateBeta(AlphaBeta *ab, int T)
{
   int t;
   PruneInfo *p;
   DVector **beta;

   p = ab->pInfo;
   p->qHi = CreateShortVec(&ab->abMem, T); /* storage for min and max q vals */
   p->qLo = CreateShortVec(&ab->abMem, T);
   beta=(DVector **)New(&ab->abMem, T*sizeof(DVector *));
   --beta;
   for (t=1;t<=T;t++){
      beta[t] = NULL;
   }

   ab->beta = beta;
}

/* CreateBetaQ: column of DVectors covering current beam */
static DVector *CreateBetaQ(MemHeap *x, int qLo,int qHi,int Q)
{
   int q;
   DVector *v;
   /*当前qLo=1*/
   qLo--; qLo--;  if (qLo<1) qLo=1;/*每次qHi都在变化，依次减1*/
   qHi++;  if (qHi>Q) qHi=Q;
   v = (DVector *)New(x, (qHi-qLo+1)*sizeof(DVector));
   v -= qLo;
   for (q=qLo;q<=qHi;q++) v[q] = NULL;
   return(v);
}
  

/* CreateOtprob: create T pointer arrays for Otprob */
static void CreateOtprob(AlphaBeta *ab, int T)
{
   int t;
   float *****otprob;
   
   otprob=(float *****)New(&ab->abMem, T*sizeof(float ****));
   --otprob;
   for (t=1;t<=T;t++){
      otprob[t] = NULL;
   }
   
   ab->otprob = otprob;

}

/* CreateOqprob: create Q pointer arrays for Otprob */
static float ****CreateOqprob(MemHeap *x, int qLo,int qHi)
{
   int q;
   float ****v;
   
   v=(float ****)New(x, (qHi-qLo+1)*sizeof(float ***));
   v-=qLo;
   for (q=qLo;q<=qHi;q++) v[q] = NULL;
   return(v);
}

/* NewBetaVec: create prob vector size 1..N */
static DVector NewBetaVec(MemHeap *x, int N)
{
   DVector v;
   
   v=(DVector)New(x, N*sizeof(double));
   --v;
   return v;
}

/* CreateOjsprob: create [2..N-1][0..S] arrays for Otprob */
static float *** CreateOjsprob(MemHeap *x, int N, int S)
{
   float ***v;
   int SS,j,s;
   
   SS=(S==1)?1:S+1;
   v=(float ***)New(x, (N-2)*sizeof(float **));
   v -= 2;
   for (j=2;j<N;j++) {
      v[j]=(float **)New(x, SS*sizeof(float *));
      if (S==1)
	 v[j][0] = NULL;
      else {
	 v[j][0] = (float *)New(x, sizeof(float));
	 v[j][0][0] = LZERO;
	 for (s=1;s<SS;s++)
	    v[j][s] = NULL;
      }
   }
   return v;
}

/* NewOtprobVec: create prob vector size [0..M] */
static float * NewOtprobVec(MemHeap *x, int M)
{
   float *v;
   int MM,m;

   MM=(M==1)?1:M+1;
   v=(float *)New(x, MM*sizeof(float));
   v[0]=LZERO;
   if (M>1)
      for (m=1;m<=M;m++)
	 v[m] = LZERO;
   return v;
}

/* ShStrP: Stream Outp calculation exploiting sharing */
static float * ShStrP(HMMSet *hset, StreamElem *ste, Vector v, int t,
		       AdaptXForm *xform, MemHeap *abmem)
{
   WtAcc *wa;
   MixtureElem *me;
   MixPDF *mp;
   float *outprobjs;
   int m,M;
   PreComp *pMix;
   LogFloat det,x,mixp,wt;
   Vector otvs;
   
   wa = (WtAcc *)ste->hook;
   if (wa->time==t) /*此判断条件可为真也可为假*/         /* seen this state before */
	  outprobjs = wa->prob;
   else {
      M = ste->nMix;
      outprobjs = NewOtprobVec(abmem,M);
      me = ste->spdf.cpdf+1;
      if (M==1){/*当前此判断条件为真*//* Single Mix Case */
         mp = me->mpdf;
         pMix = (PreComp *)mp->hook;
         if ((pMix != NULL) && (pMix->time == t))/*当前此判断条件为假*/
            x = pMix->prob;            
         else {/*当前执行此分支*/
            x = MOutP(ApplyCompFXForm(mp,v,xform,&det,t),mp);/* returns log prob of vector x for given mixture */
			/*当前ApplyCompFXForm的执行结果是返回Vector v,且det=0*/	
            x += det;
            if (pMix != NULL) {/*当前此判断条件为真*/
               pMix->prob = x; pMix->time = t;
            }
         }
      } else if (sharedMix) { /*当前以下部分不执行此处*//* Multiple Mixture Case - general case */
         x = LZERO;
         for (m=1;m<=M;m++,me++) {
            wt = MixLogWeight(hset,me->weight);
            if (wt>LMINMIX){
               mp = me->mpdf;
               pMix = (PreComp *)mp->hook;
               if ((pMix != NULL) && (pMix->time == t))
                  mixp = pMix->prob;
               else {
                  mixp = MOutP(ApplyCompFXForm(mp,v,xform,&det,t),mp);
		  mixp += det;
                  if (pMix != NULL) {
                     pMix->prob = mixp; pMix->time = t;
                  }
               }
               x = LAdd(x,wt+mixp);
	       outprobjs[m] = mixp;
            }
         }
      } else if (!pde) { /* Multiple Mixture Case - no shared mix case */
         x = LZERO;
         for (m=1;m<=M;m++,me++) {
            wt = MixLogWeight(hset,me->weight);
            if (wt>LMINMIX){
               mp = me->mpdf;
	       mixp = MOutP(ApplyCompFXForm(mp,v,xform,&det,t),mp);
	       mixp += det;
               x = LAdd(x,wt+mixp);
	       outprobjs[m] = mixp;
            }
         }
      } else {    /* Partial distance elimination */
	 /* first Gaussian computed exactly in PDE */
	 wt = MixLogWeight(hset,me->weight);
	 mp = me->mpdf;
	 otvs = ApplyCompFXForm(mp,v,xform,&det,t);
	 mixp = IDOutP(otvs,VectorSize(otvs),mp); /* INVDIAGC assumed */
	 mixp += det;
	 x = wt+mixp;
	 outprobjs[1] = mixp;
	 for (m=2,me=ste->spdf.cpdf+2;m<=M;m++,me++) {
            wt = MixLogWeight(hset,me->weight);
	    if (wt>LMINMIX){
	       mp = me->mpdf;
	       otvs = ApplyCompFXForm(mp,v,xform,&det,t);
	       if (PDEMOutP(otvs,mp,&mixp,x-wt-det) == TRUE) {
		  mixp += det;
		  x = LAdd(x,wt+mixp);
	       }
	       outprobjs[m] = mixp; /* LZERO if PDEMOutP returns FALSE */
	    }
	 }
      }/*当前从下面开始执行*/
      outprobjs[0] = x;
      wa->prob = outprobjs;
      wa->time = t;
   }
   return outprobjs;
}
   
/* Setotprob: allocate and calculate otprob matrix at time t */
static void Setotprob(AlphaBeta *ab, FBInfo *fbInfo, ParmBuf pbuf, 
                      Observation ot, int t, int S, int qHi, int qLo)
{
   int q,j,Nq,s;
   float ***outprob, **outprobj, *****otprob;
   StreamElem *ste;
   HLink hmm;
   LogFloat sum;
   PruneInfo *p;
   int skipstart, skipend;
   HMMSet *hset;
   
   hset = fbInfo->al_hset;
   skipstart = fbInfo->skipstart;
   skipend = fbInfo->skipend;
   p = ab->pInfo;
   otprob = ab->otprob;
   /*otprob是array[1..T][1..Q][2..Nq-1][0..S][0..M]of prob,T是帧数，
   Q是模型数（多个模型一起估计），Nq是状态数，S是流数，M是混合数*/
   ReadAsTable(pbuf,t-1,&ot);/*将pbuf中的数据读入ot中*/
   if (hset->hsKind == TIEDHS)/*当前不此判断条件为假*/
	  PrecomputeTMix(hset,&ot,pruneSetting.minFrwdP,0);	  
   if (trace&T_OUT && NonSkipRegion(skipstart,skipend,t)) 
      printf(" Output Probs at time %d\n",t);
   if (qLo>1) --qLo;
   otprob[t] = CreateOqprob(&ab->abMem,qLo,qHi);
   /*otprob[t]是array[1..Q][2..Nq-1][0..S][0..M]of prob*/
   for (q=qHi;q>=qLo;q--) {
      if (trace&T_OUT && NonSkipRegion(skipstart,skipend,t)) 
         printf(" Q%2d: ",q);