hfb.c

隐马尔科夫模型工具箱

{
   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);
      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);

      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;

   return(qt);
}

/* 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));
   ab->alphat  = alphat;
   ab->alphat1 = alphat1;
}


/* 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;
      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]:LZERO;
      }
      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;
      a1N = hmm->transP[1][Nq];
   }
   ZeroAlpha(ab,eq+1,Q);
   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;
}

/* 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;

   for (q = sq; q <= eq; q++) {
      lNq = Nq; hmm = ab->al_qList[q]; Nq = hmm->numStates; 
      aq = alphat[q]; 
      laq = alphat1[q];
      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)
         aq[1] = LZERO;
      else{
         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++) {
         a = hmm->transP[1][j];
         x = (a>LSMALL)?a+aq[1]:LZERO;
         for (i=2;i<Nq;i++){
            a = hmm->transP[i][j]; y = laq[i];
            if (a>LSMALL && y>LSMALL)
               x = LAdd(x,y+a);
         }
         aq[j] = x + outprob[j][0];
      }
      x = LZERO;
      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];
   }
   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--; qLo--;  if (qLo<1) qLo=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;
}

/* NewOtprobVec: create prob matrix size [2..N-1][0..S] */
static float ** NewOtprobVec(MemHeap *x, int N, int S)
{
   float **v;
   int SS,i;
   
   SS=(S==1)?1:S+1;
   v=(float **)New(x, (N-2)*sizeof(float *));
   v -= 2;
   for (i=2;i<N;i++)
      v[i]=(float *)New(x, SS*sizeof(float));
   return v;
}

/* ShStrP: Stream Outp calculation exploiting sharing */
static LogFloat ShStrP(HMMSet *hset, Vector v, int t, HLink hmm, 
                       int state, int stream)
{
   WtAcc *wa;
   StreamElem *ste;
   MixtureElem *me;
   MixPDF *mp;
   int m,M;
   PreComp *pMix;
   LogFloat x,mixp,wt;
   
   ste = hmm->svec[state].info->pdf+stream;
   wa = (WtAcc *)ste->hook;
   if (wa->time==t)           /* seen this state before */
      x = wa->prob;
   else {
      M = ste->nMix;
      me = ste->spdf.cpdf+1;
      if (M==1){                 /* Single Mix Case */
         mp = me->mpdf;
         pMix = (PreComp *)mp->hook;
         if (pMix->time == t)
            x = pMix->prob;
         else {
            x = MOutP(v,mp);