hfb.c

实现HMM算法

      CreateHeap(&utt->dataStack2,"dataStore2",MSTAK, 1, 0.5, 1000,  10000);
   utt->pbuf = NULL; utt->pbuf2 = NULL; utt->tr = NULL;
}

/* InitPruneStats: initialise pruning stats */
static void InitPruneStats(AlphaBeta *ab)
{
   PruneInfo *p;

   ab->pInfo = (PruneInfo *) New(&ab->abMem, sizeof(PruneInfo));
   p = ab->pInfo;
   p->maxBeamWidth = 0;
   p->maxAlphaBeta = LZERO;
   p->minAlphaBeta = 1.0;
}


/* -------------------------- Trace Support ----------------------- */

/* CreateTraceOcc: create the array of acc occ counts */
static void CreateTraceOcc(AlphaBeta *ab, UttInfo *utt)
{
   int q;
   Vector *occa;

   printf("\n");
   ab->occa=(Vector *)New(&ab->abMem, utt->Q*sizeof(Vector));
   occa = ab->occa;
   --occa;
   for (q=1;q<=utt->Q;q++){
      occa[q] = CreateVector(&ab->abMem, ab->al_qList[q]->numStates);
      ZeroVector(occa[q]);
   }
}

/* TraceOcc: print current accumulated occ counts for all models */
static void TraceOcc(AlphaBeta *ab, UttInfo *utt, int t)
{
   int Nq, q, i;
   Vector occaq;
   HLink hmm;
   float max;

   printf("Accumulated Occ Counts at time %d\n",t);
   for (q=1; q<=utt->Q; q++){
      occaq = ab->occa[q]; hmm = ab->al_qList[q]; Nq = hmm->numStates;
      max = 0.0;        /* ignore zero vectors */
      for (i=1;i<=Nq;i++)
         if (occaq[i]>max) max = occaq[i];
      if (max>0.0) {    /* not zero so print it */
         printf("  Q%2d: %5s", q,ab->qIds[q]->name);
         for (i=1;i<=Nq;i++) printf("%7.2f",occaq[i]);
         printf("\n");
      }
   }
}

/* SetOcct: set the global occupation count for given hmm */
static void SetOcct(HLink hmm, int q, Vector occt, Vector *occa,
                    DVector aqt, DVector bqt, DVector bq1t, LogDouble pr)
{
   int i,N;
   double x;
   Vector occaq;
   
   N=hmm->numStates;
   for (i=1;i<=N;i++) {
      x = aqt[i]+bqt[i];/*计算了HTKBOOK第8页中第(1.25)个公式*/
      if (i==1 && bq1t != NULL && hmm->transP[1][N] > LSMALL)/*当前此判断条件为假*/
         x = LAdd(x,aqt[1]+bq1t[1]+hmm->transP[1][N]);
      x -= pr;/*计算了HTKBOOK第8页中第(1.26)个公式*/
      occt[i] = (x>MINEARG) ? exp(x) : 0.0;/*将对数概率恢复原始概率*/
   }
   if (trace&T_OCC) {/*当前不执行这里*/
      occaq = occa[q];
      for (i=1;i<=N;i++) occaq[i] += occt[i];
   }
}


/* NonSkipRegion: returns true if t is not in the skip region */
static Boolean NonSkipRegion(int skipstart, int skipend, int t)
{
   return skipstart<1 || t<skipstart || t>skipend;
}

/* PrLog: print a log value */
void PrLog(LogDouble x)
{
   if (x<LSMALL)
      printf("       LZERO");
   else
      printf("%12.5f",x);
}



/* -------------------------------------------------------------------*/

/* GetInputObs: Get input Observations for t */
void GetInputObs( UttInfo *utt, int t, HSetKind hsKind )
{

   if (utt->twoDataFiles)
      ReadAsTable(utt->pbuf2,t-1,&(utt->ot2));
   ReadAsTable(utt->pbuf,t-1,&(utt->ot));

   if (hsKind == TIEDHS)
      if (utt->twoDataFiles)
         ReadAsTable(utt->pbuf,t-1,&(utt->ot2));

}

/* --------------------------- Forward-Backward --------------------- */



/* CheckPruning: record peak alpha.beta product and position */
static void CheckPruning(AlphaBeta *ab, int t, int skipstart, int skipend)
{
   int i,q,Nq,bestq,besti,margin;
   PruneInfo *p;
   DVector aq,bq;
   HLink hmm;
   LogDouble lx,maxL;

   bestq = besti = 0;
   maxL = LZERO;
   p = ab->pInfo;
   for (q=p->qLo[t];q<=p->qHi[t];q++){
      hmm = ab->al_qList[q]; Nq = hmm->numStates;   
      aq = ab->alphat[q]; bq=ab->beta[t][q];
      for (i=2;i<Nq;i++){
         if ((lx=aq[i]+bq[i])>maxL){
            bestq = q; besti = i; maxL=lx;
         }
      }
   }
   if (maxL > p->maxAlphaBeta) p->maxAlphaBeta = maxL;
   if (maxL < p->minAlphaBeta) p->minAlphaBeta = maxL;
   margin = p->qHi[t] - p->qLo[t]+1;
   if (margin>p->maxBeamWidth) p->maxBeamWidth = margin;
   if (NonSkipRegion(skipstart, skipend, t)){
      if (bestq == 0) 
         printf("%3d. No max found in alpha.beta\n",t);
      else
         printf("%3d. Max Alpha.Beta = %9.4f at q=%d i=%d [%s]\n",
                t,maxL,bestq,besti,ab->qIds[bestq]->name);
   }
}

/* SummarisePruning: print out pruning stats */
static void SummarisePruning(PruneInfo *p, int Q, int T)
{
   long totalQ=0;
   float e;
   int t;
   
   for (t=1;t<=T;t++)
      totalQ += p->qHi[t]-p->qLo[t]+1;
   e = (1.0 - (float) totalQ / ((float) T*Q)) * 100.0;
   printf(" Pruning %.1f%%; MaxBeamWidth %d; PeakShortFall %.2f\n",
          e,p->maxBeamWidth,p->maxAlphaBeta - p->minAlphaBeta);
   fflush(stdout);
}

/* CreateInsts: create array of hmm instances for current transcription */
static int CreateInsts(FBInfo *fbInfo, AlphaBeta *ab, int Q, Transcription *tr)
{
   int q,qt;
   LLink lab;
   MLink macroName;
   TrAcc *ta;
   LabId  *qIds;
   short *qDms;
   HLink *al_qList, *up_qList;
   HMMSet *al_hset, *up_hset;

   al_hset=fbInfo->al_hset; up_hset=fbInfo->up_hset;
   /* init logical hmm list */
   up_qList=(HLink *)New(&ab->abMem, Q*sizeof(HLink)); 
   --up_qList;

   /* 2-model re-estimation update models */
   if (fbInfo->twoModels) {/*当前这里不被执行*/
      al_qList=(HLink *)New(&ab->abMem, Q*sizeof(HLink));
      --al_qList;
   }
   else /* use same list for update and align */
      al_qList = up_qList;

   qIds = (LabId *)New(&ab->abMem, Q*sizeof(LabId));
   --qIds;
   qDms = (short *)New(&ab->abMem, Q*sizeof(short));
   --qDms;

   qt=0;
   for (lab=tr->head->head->succ,q=1; lab->succ!= NULL; lab=lab->succ,q++){
      /* align models */
      if((macroName=FindMacroName(al_hset,'l',lab->labid))==NULL)
         HError(7321,"CreateInsts: Unknown label %s",lab->labid->name);/*执行后macroName->id->name的值类似为OO，即单音素名称*/
      al_qList[q] = (HLink)macroName->structure;
      /* 2-model re-estimation update models */
      if (fbInfo->twoModels){/*当前不执行这里*/
         if((macroName=FindMacroName(up_hset,'l',lab->labid))==NULL)
            HError(2321,"CreateInsts: Unknown update label %s",lab->labid->name);
         up_qList[q] = (HLink)macroName->structure;
         /* need equal num states */
         if ((al_qList[q])->numStates != (up_qList[q])->numStates)
            HError(999,"Num states differ in align and update models (%d %d)",
                   (al_qList[q])->numStates,(up_qList[q])->numStates);      
      }
      qIds[q] = macroName->id;
      ta = (TrAcc *)GetHook(al_qList[q]->transP);
      qt += (qDms[q] = ta->minDur);
	  /*ta->minDur中保存的是Min no of frames to get through trans mat，当前开始和结束时取值为2，中间取值为3 */
      
      if (q>1 && qDms[q]==0 && qDms[q-1]==0)
         HError(7332,"CreateInsts: Cannot have successive Tee models");
      if (al_hset->hsKind==SHAREDHS)/*当前此判断条件为真*/
		 ResetHMMPreComps(al_qList[q],al_hset->swidth[0]);
   }
   if ((qDms[1]==0)||(qDms[Q]==0))
      HError(7332,"CreateInsts: Cannot have Tee models at start or end of transcription");

   /* store in struct*/
   ab->al_qList = al_qList; 
   ab->up_qList = up_qList;
   ab->qIds  = qIds;
   ab->qDms  = qDms;/*qDms是一个矩阵，长度为Q，qDms[q]表示第q个模型的最小驻留*/

   return(qt);/*qt中保存了总的Min no of frames to get through trans mat*/
}

/* CreateAlpha: allocate alpha columns */
static void CreateAlpha(AlphaBeta *ab, HMMSet *hset, int Q)
{

   int q;
   DVector *alphat, *alphat1;
 
   /* Create Storage Space - two columns */
   alphat = (DVector *)New(&ab->abMem, Q*sizeof(DVector));
   --alphat;
   for (q=1;q<=Q;q++)
       alphat[q] = CreateDVector(&ab->abMem, (ab->al_qList[q])->numStates);
   alphat1=(DVector *)New(&ab->abMem, Q*sizeof(DVector));
   --alphat1;
   for (q=1;q<=Q;q++)
       alphat1[q] = CreateDVector(&ab->abMem, (ab->al_qList[q])->numStates);

   ab->occt = CreateVector(&ab->abMem,MaxStatesInSet(hset));/*occ probs for current time t*/
   ab->alphat  = alphat;/*array[1..Q][1..Nq] of prob*/
   ab->alphat1 = alphat1;/*alpha[t-1]*/
}


/* ZeroAlpha: zero alpha's of given models */
static void ZeroAlpha(AlphaBeta *ab, int qlo, int qhi)
{
   HLink hmm;
   int Nq,j,q;
   DVector aq;
   
   for (q=qlo;q<=qhi;q++) {   
      hmm = ab->al_qList[q]; 
      Nq = hmm->numStates; 
      aq = ab->alphat[q];
      for (j=1;j<=Nq;j++)
         aq[j] = LZERO;
   }
}

/* InitAlpha: initialise alpha columns for time t=1 */
static void InitAlpha(AlphaBeta *ab, int *start, int *end, 
                      int Q, int skipstart, int skipend)
{
   int i,j,Nq,eq,q;
   PruneInfo *p;
   HLink hmm;
   DVector aq;
   float ***outprob;
   LogDouble x,a,a1N=0.0;
   
   p = ab->pInfo;
   eq = p->qHi[1];
   for (q=1; q<=eq; q++){
      hmm = ab->al_qList[q]; Nq = hmm->numStates;
      aq = ab->alphat[q];
      aq[1] = (q==1)?0.0:ab->alphat[q-1][1]+a1N;/*初始化HTKbook第129页Alpha初始化的第(1)个公式*/
      if((outprob = ab->otprob[1][q]) == NULL)
         HError(7322,"InitAlpha: Outprob NULL in model %d in InitAlpha",q);
      for (j=2;j<Nq;j++) {
         a = hmm->transP[1][j];
         aq[j] = (a>LSMALL)?aq[1]+a+outprob[j][0][0]:LZERO;/*初始化HTKbook第129页Alpha初始化的第(2)个公式*/
      }
      x = LZERO;
      for (i=2;i<Nq;i++) {
         a = hmm->transP[i][Nq];
         if (a>LSMALL)
            x = LAdd(x,aq[i]+a);
      }
      aq[Nq] = x;/*初始化HTKbook第129页Alpha初始化的第(1)个公式*/
      a1N = hmm->transP[1][Nq];
   }
   ZeroAlpha(ab,eq+1,Q);/*其它值都为0*/
   if (trace&T_PRU && p->pruneThresh < NOPRUNE)
      CheckPruning(ab,1,skipstart,skipend);
   *start = 1; *end = eq;
}

/* MaxModelProb: Calc max probability of being in model q at
   time t, return LZERO if cannot do so */
static LogDouble MaxModelProb(AlphaBeta *ab, int q, int t, int minq)
{
   DVector aq,bq,bq1;
   LogDouble maxP,x;
   int Nq1,Nq,i,qx,qx1;
   HLink hmm;
   
   if (q==1)
      maxP = LZERO;
   else {
      bq1 = ab->beta[t][q-1]; Nq1 = ab->al_qList[q-1]->numStates;
      maxP = (bq1==NULL)?LZERO:ab->alphat[q-1][Nq1] + bq1[Nq1];
      for (qx=q-1;qx>minq && ab->al_qList[qx]->transP[1][Nq1] > LSMALL;qx--){
         qx1 = qx-1;
         bq1 = ab->beta[t][qx1]; Nq1 = ab->al_qList[qx1]->numStates;
         x=(bq1==NULL)?LZERO:ab->alphat[qx1][Nq1]+bq1[Nq1];
         if (x > maxP) maxP = x;
      }
   }
   hmm = ab->al_qList[q]; Nq = hmm->numStates;   
   bq=ab->beta[t][q];
   if (bq != NULL) {
      aq = ab->alphat[q]; 
      for (i=1;i<Nq;i++)
         if ((x=aq[i]+bq[i]) > maxP) maxP = x;
   }
   return maxP;