1,1347c
/* ----------------------------------------------------------- */
/*                                                             */
/*                          ___                                */
/*                       |_| | |_/   SPEECH                    */
/*                       | | | | \   RECOGNITION               */
/*                       =========   SOFTWARE                  */ 
/*                                                             */
/*                                                             */
/* ----------------------------------------------------------- */
/*         Copyright: Microsoft Corporation                    */
/*          1995-2000 Redmond, Washington USA                  */
/*                    http://www.microsoft.com                 */
/*                                                             */
/*   Use of this software is governed by a License Agreement   */
/*    ** See the file License for the Conditions of Use  **    */
/*    **     This banner notice must not be removed      **    */
/*                                                             */
/* ----------------------------------------------------------- */
/*         File: HDRest.c: HMM initialisation program          */
/* ----------------------------------------------------------- */

char *hrest_version = "!HVER!HRest:   3.1 [CUED 16/01/02]";
char *hrest_vc_id = "$Id: HRest.c,v 1.7 2002/01/16 18:11:29 ge204 Exp $";

/*
   This program is used to estimate the transition parameters,
   means, covariances and mixture weights of a
   hidden Markov model using Baum-Welch reestimation.
   It handles multiple streams and tying but it ignores
   stream weights by assuming that they are all unity.
*/

/* Trace Flags */
#define T_TOP    0001    /* Top level tracing */
#define T_LD0    0002    /* File Loading */
#define T_LD1    0004    /* + segments within each file */
#define T_OTP    0010    /* Observation Probabilities */
#define T_ALF    0020    /* Alpha matrices */
#define T_BET    0040    /* Beta matrices */
#define T_OCC    0100    /* Occupation Counters */
#define T_TAC    0200    /* Transition Counters */
#define T_MAC    0400    /* Mean Counters */
#define T_VAC   01000    /* Variance Counters */
#define T_WAC   02000    /* MixWeight Counters */
#define T_TRE   04000    /* Reestimated transition matrix */
#define T_WRE  010000    /* Reestimated mixture weights */
#define T_MRE  020000    /* Reestimated means */
#define T_VRE  040000    /* Reestimated variances */
#define T_LGP 0100000    /* Compare LogP via alpha and beta */


#include "HShell.h"     /* HMM ToolKit Modules */
#include "HMem.h"
#include "HMath.h"
#include "HSigP.h"
#include "HAudio.h"
#include "HWave.h"
#include "HVQ.h"
#include "HParm.h"
#include "HLabel.h"
#include "HModel.h"
#include "HTrain.h"
#include "HUtil.h"


/* Global Settings */
static char * segLab = NULL;     /* segment label if any */
static LabId  segId  = NULL;     /* and its index */
static char * labDir = NULL;     /* label file directory */
static char * labExt = "lab";    /* label file extension */
static char * outDir = NULL;     /* output macro file directory, if any */
static int  maxIter  = 20;       /* max iterations in parameter estimation */
static float epsilon = 1.0E-4;   /* convergence criterion */
static double DURFLOOR=1e-30;  /* an emperical(arbitrary) const used to set the floor values of duration PMF */
static int minSeg    = 3;        /* min segments to train a model */
static Boolean firstTime = TRUE; /* Flag used to enable InitSegStore */
static Boolean saveBinary = FALSE;  /* save output in binary  */
static FileFormat dff=UNDEFF;    /* data file format */
static FileFormat lff=UNDEFF;    /* label file format */
static float minVar  = 0.0;      /* minimum variance */
static float mixWeightFloor=0.0; /* Floor for mixture weights */
static float tMPruneThresh = 10.0;    /* tied mix prune threshold */
static char *hmmfn;              /* HMM definition file name */
static char *outfn=NULL;         /* output definition file name */
enum _UPDSet{UPMEANS=1,UPVARS=2,UPTRANS=4,UPMIXES=8, UPDURS=16};
typedef enum _UPDSet UPDSet;
static UPDSet uFlags = (UPDSet) (UPMEANS|UPVARS|UPTRANS|UPMIXES|UPDURS);     /* update flags */
static int  trace    = 0;        /* Trace level */
static ConfParam *cParm[MAXGLOBS];   /* configuration parameters */
static int nParm = 0;               /* total num params */
static Boolean segReject = TRUE; /* Enable short train segment rejection */

static int maxD = 0;		/* max duration length */ 

/* Global Data Structures */
static HMMSet hset;        /* The current unitary hmm set */
static HLink hmm;          /* link to the hmm itself, HLink=HMMDef, single HMM model */
static int nStates;        /* numStates of hmm */
static int nStreams;       /* numStreams of hmm */
static HSetKind hsKind;    /* kind of the HMM system: //PLAINHS, SHAREDHS, TIEDHS, DISCRETEHS*/
static int maxMixes;       /* max num mixtures across all streams */
static int maxMixInS[SMAX];/* array[1..swidth[0]] of max mixes */
static int nSeg;           /* num training segments */
static int nTokUsed;       /* actual number of tokens(segments) used */
static int maxT,minT,T;    /* max,min and current segment lengths */
static DMatrix alpha;      /* array[1..nStates][1..maxT] of forward prob */
static DMatrix beta;       /* array[1..nStates][1..maxT] of backward prob */

static DMatrix alphaStar;  /* array[1..nStates][1..maxT] of the second forward prob */
static DMatrix betaStar;   /* array[1..nStates][1..maxT] of the second backward prob */
static DMatrix exptimeStart;/* array[1..nStates][1..maxT] of the expected # times started at t or before | O */
static DMatrix exptimeEnd;  /* array[1..nStates][1..maxT] of the expected # times terminated before t| O */
static DMatrix gamma;	   /* array[1..nStates][1..maxT] of the expected # times terminated before t| O */


static Matrix outprob;     /* array[2..nStates-1][1..maxT] of output prob */

static Matrix durprob;	   /* array[2..nStates-1][1..maxD] of duration prob */

static Vector **stroutp;   /* array[1..maxT][2..nStates-1][1..nStreams] ...*/
                           /* ... of streamprob */
static Matrix **mixoutp;   /* array[2..nStates-1][1..maxT][1..nStreams]
                              [1..maxMixes] of mixprob */
/*static Matrix **gammaMix; */  /* array[2..nStates-1][1..maxT][1..nStreams]
                              //[1..maxMixes] of mixprob */

static Vector occr;        /* array[1..nStates-1] of occ count for cur time */
static Vector zot;         /* temp storage for zero mean obs vector */
static Vector numZeroEnd_alpha;	/*temp buffer used to restrict the special ending of alpha in duration model*/ 
static Vector numZeroBegin_betaS;/*temp buffer used to restrict the special beginning of beta in duration model*/ 

static Vector vFloor[SMAX];      /* variance floor - default is all zero */
static float vDefunct=0.0;       /* variance below which mixture defunct */

static SegStore segStore;        /* Storage for data segments, htrain */
static MemHeap segmentStack;     /* Used by segStore, hmem */
static MemHeap alphaBetaStack;   /* For storage of alpha and beta probs */
static MemHeap accsStack;        /* For storage of accumulators */
static MemHeap transStack;       /* For storage of transcription */
static MemHeap bufferStack;      /* For storage of buffer */
static ParmBuf pbuf;             /* Currently input parm buffer */

   
/* ------------------ Process Command Line ------------------------- */

/* SetConfParms: set conf parms relevant to HRest  */
void SetConfParms(void)
{
   int i;
   double d;
   Boolean b;

   nParm = GetConfig("HREST", TRUE, cParm, MAXGLOBS);  /* hshell */
   if (nParm>0) {
      if (GetConfInt(cParm,nParm,"TRACE",&i)) trace = i;
      if (GetConfBool(cParm,nParm,"SAVEBINARY",&b)) saveBinary = b;
      if (GetConfFlt(cParm,nParm,"VDEFUNCT",&d)) vDefunct = d;
   }
}

void ReportUsage(void)
{
   printf("\nUSAGE: HDRest [options] hmmFile trainFiles...\n\n");
   printf(" Option                                       Default\n\n");
   printf(" -e f    Set convergence factor epsilon       1.0E-4\n");
   printf(" -i N    Set max iterations to N              20\n");
   printf(" -l s    Set segment label to s               none\n");
   printf(" -m N    Set min segments needed              3\n");
   printf(" -t      Disable short segment rejection      on\n");
   printf(" -u tmvw Update t)rans m)eans v)ars w)ghts    tmvw\n");
   printf(" -v f    Set minimum variance to f            0.0\n");
   printf(" -c f    Tied Mixture pruning threshold       10.0\n");
   printf(" -w f    Set mix wt floor to f x MINMIX       0.0\n");
   PrintStdOpts("BFGHILMSTX");
   printf("\n\n");
}

void SetuFlags(void)
{
   char *s;
   
   s=GetStrArg();    /* hshell */
   uFlags=(UPDSet) 0;        
   while (*s != '\0')
      switch (*s++) {
      case 't': uFlags = (UPDSet) (uFlags+UPTRANS); break;
      case 'm': uFlags = (UPDSet) (uFlags+UPMEANS); break;
      case 'v': uFlags = (UPDSet) (uFlags+UPVARS); break;
      case 'w': uFlags = (UPDSet) (uFlags+UPMIXES); break;
      default: HError(2220,"SetuFlags: Unknown update flag %c",*s);
         break;
      }
}

int main(int argc, char *argv[])
{
   char *datafn, *s;
   void Initialise1(void);
   void Initialise2(void);
   void LoadFile(char *fn);
   void ReEstimateModel(void);
   void SaveModel(char *outfn);
 
   if(InitShell(argc,argv,hrest_version,hrest_vc_id)<SUCCESS)  /* hshell */
      HError(2200,"HRest: InitShell failed");

   InitMem();   InitLabel();  
   InitMath();  InitSigP();
   InitWave();  InitAudio();
   InitVQ();    InitModel();
   if(InitParm()<SUCCESS)  
      HError(2200,"HRest: InitParm failed");

   InitTrain(); InitUtil();

   if (!InfoPrinted() && NumArgs() == 0)   /* hshell */
      ReportUsage();
   if (NumArgs() == 0) Exit(0);   /* hshell */
   SetConfParms();		/* hrest */
   CreateHMMSet(&hset,&gstack,FALSE);    /* hmodel, gstack<-hmem */

   /*//parse the switches for HRest like HRest -e -i -l ....
   //NextArg --- Next Argument*/
   while (NextArg() == SWITCHARG) {
      s = GetSwtArg();		/* hshell */
      if (strlen(s)!=1) 
         HError(2219,"HRest: Bad switch %s; must be single letter",s);
      switch(s[0]){      
      case 'e':
         epsilon = GetChkedFlt(0.0,1.0,s); break;		/* hshell */
      case 'i':
         maxIter = GetChkedInt(1,100,s); break;		/* hshell */
      case 'l':
         if (NextArg() != STRINGARG)
            HError(2219,"HRest: Segment label expected");
         segLab = GetStrArg();		/* hshell */
         break;
      case 'm':
         minSeg = GetChkedInt(1,1000,s); break;		/* hshell */
      case 't':
         segReject = FALSE;
         break;
      case 'u':
         SetuFlags(); break;	/* hrest */
      case 'v':
         minVar = GetChkedFlt(0.0,100.0,s); break;		/* hshell */
      case 'c':
         tMPruneThresh = GetChkedFlt(0.0,1000.0,s); break;		/* hshell */
      case 'w':
         mixWeightFloor = MINMIX * GetChkedFlt(0.0,10000.0,s); break;	/* hshell */
      case 'B':
         saveBinary = TRUE;
         break;
      case 'F':
         if (NextArg() != STRINGARG)	/* hshell */
            HError(2219,"HRest: Data File format expected");
         if((dff = Str2Format(GetStrArg())) == ALIEN)	/*hwave, hshell*/
            HError(-2289,"HRest: Warning ALIEN Data file format set");
         break;
      case 'G':
         if (NextArg() != STRINGARG)
            HError(2219,"HRest: Label File format expected");
         if((lff = Str2Format(GetStrArg())) == ALIEN)
            HError(-2289,"HRest: Warning ALIEN Label file format set");
         break;
      case 'H':
         if (NextArg() != STRINGARG)
            HError(2219,"HRest: HMM macro file name expected");
         AddMMF(&hset,GetStrArg());		/*hmodel, hshell*/
         break;
      case 'I':
         if (NextArg() != STRINGARG)
            HError(2219,"HRest: MLF file name expected");
         LoadMasterFile(GetStrArg());	/*hlabel, hshell*/
         break;
      case 'L':
         if (NextArg()!=STRINGARG)
            HError(2219,"HRest: Label file directory expected");
         labDir = GetStrArg(); break;
      case 'M':
         if (NextArg()!=STRINGARG)
            HError(2219,"HRest: Output macro file directory expected");
         outDir = GetStrArg(); break;  
      case 'T':
         trace = GetChkedInt(0,0100000,s); break;
      case 'X':
         if (NextArg()!=STRINGARG)
            HError(2219,"HRest: Label file extension expected");
         labExt = GetStrArg(); break;
      default:
         HError(2219,"HRest: Unknown switch %s",s);
      }
   }
   if (NextArg()!=STRINGARG)
      HError(2219,"HRest: source HMM file name expected");
   hmmfn = GetStrArg();
   if (outfn==NULL) outfn = hmmfn;
   Initialise1();	/*hrest*/
   do {
      if (NextArg()!=STRINGARG)
         HError(2219,"HRest: training data file name expected");
      datafn = GetStrArg();
      LoadFile(datafn);		/*load in a single Seg each time; hrest*/
   } while (NumArgs()>0);
   nSeg = NumSegs(segStore);	/*htrain*/
   if (nSeg < minSeg)
      HError(2221,"HRest: Too Few Training Examples [%d]",nSeg);
   Initialise2();		/*hrest*/
   if (trace&T_TOP) {
      printf("%d Examples loaded, Max length = %d, Min length = %d\n",
             nSeg, maxT,minT);
      fflush(stdout);
   }
   ReEstimateModel();		/*hrest*/

   if(SaveHMMSet(&hset,outDir,NULL,saveBinary)<SUCCESS)		/*hmodel*/
      HError(2211,"HRest: SaveHMMSet failed");

   Exit(0);
   return (0);          /* never reached -- make compiler happy */
}

/* ------------------------ Initialisation ----------------------- */

/* PrintInitialInfo: print a header of program settings */
void PrintInitialInfo(void)
{   
   printf("Reestimating HMM %s . . . \n",hmmfn);
   PrintHMMProfile(stdout, hmm);
   if (segLab==NULL)
      printf(" SegLab   :  NONE\n");
   else
      printf(" SegLab   :  %s\n",segLab);
   printf(" MaxIter  :  %d\n",maxIter);
   printf(" Epsilon  :  %f\n",epsilon);
   printf(" Updating :  ");
   if (uFlags&UPTRANS) printf("Transitions "); 
   if (uFlags&UPMEANS) printf("Means "); 
   if (uFlags&UPVARS)  printf("Variances "); 
   if (uFlags&UPMIXES && maxMixes>1)  printf("MixWeights"); 
   printf("\n\n");
   printf(" - system is ");
   switch (hset.hsKind){
   case PLAINHS:  printf("PLAIN\n");  break;
   case SHAREDHS: printf("SHARED\n"); break;
   case TIEDHS:   printf("TIED\n");   break;
   case DISCRETEHS: printf("DISCRETE\n"); break;
   }   
   fflush(stdout);
}
   
/* Initialise1: 1st phase of init prior to loading dbase */
void Initialise1(void)
{
   MLink link; /*typedef struct _MacroDef *MLink; hmodel */
   LabId  hmmId;  /*typedef NameCell *LabId; NameCell is a Hash Table linked list item; hlabel */
   char base[MAXSTRLEN];
   char path[MAXSTRLEN];
   char ext[MAXSTRLEN];
   int s,i;

   /* Create and Load HMM def from hmmfn into hset*/
   if(MakeOneHMM( &hset,BaseOf(hmmfn,base))<SUCCESS)	/*hmodel,hshell*/
      HError(2128,"Initialise1: MakeOneHMM failed");
   if(LoadHMMSet( &hset,PathOf(hmmfn,path),ExtnOf(hmmfn,ext))<SUCCESS)		/*hmodel,hshell*/
      HError(2128,"Initialise1: LoadHMMSet failed");
   if (hset.hsKind!=PLAINHS)
      uFlags = (UPDSet) (uFlags & (~(UPMEANS|UPVARS)));	/*update mean and var only when it is PLAINHS*/


   /* Get a pointer to the physical HMM and set related globals */
   hmmId = GetLabId(base,FALSE);	/*hlabel*/
   link = FindMacroName(&hset,'h',hmmId);	/*hmodel*/
   hmm = (HLink)link->structure;  
   nStates = hmm->numStates;

   /*maxD = larget maximum allowed duration in this HMM*/
   for (i=2;i<nStates;i++)
	   if(maxD<hmm->svec[i].info->maxD) maxD=hmm->svec[i].info->maxD;

   nStreams = hset.swidth[0];
   hsKind = hset.hsKind;
   
   /* Stacks for global structures requiring memory allocation */
   CreateHeap(&segmentStack,"SegStore", MSTAK, 1, 0.0, 100000, LONG_MAX);	/*hmem*/
   CreateHeap(&alphaBetaStack,"AlphaBetaStore", MSTAK, 1, 0.0, 1000, 1000);
   /*void CreateHeap(MemHeap *x, char *name, HeapType type, size_t elemSize, 
                float growf, size_t numElem,  size_t maxElem);  */
   CreateHeap(&accsStack,"AccsStore", MSTAK, 1, 0.0, 1000, 1000);
   CreateHeap(&transStack,"TransStore", MSTAK, 1, 0.0, 1000, 1000);   /*transcription stack*/
   CreateHeap(&bufferStack,"BufferStore", MSTAK, 1, 0.0, 1000, 1000);
   AttachAccs(&hset, &accsStack);	/*htrain*/

   SetVFloor( &hset, vFloor, minVar);	/*hmodel*/

   if(segLab != NULL)
      segId = GetLabId(segLab,TRUE);
   if(trace&T_TOP)
      PrintInitialInfo();	/*hrest*/

   maxMixes = MaxMixtures(hmm);
   for(s=1; s<=nStreams; s++)
      maxMixInS[s] = MaxMixInS(hmm, s);	/*maximum number of mixtures in  stream*/
   T = maxT = 0; minT = 100000;
}

/* Initialise2: 2nd phase of init after loading dbase */
void Initialise2(void)
{
   int t,j,m,s;

   alpha = CreateDMatrix(&alphaBetaStack,nStates,maxT);	/*hmem*/
   beta = CreateDMatrix(&alphaBetaStack,nStates,maxT);
   alphaStar = CreateDMatrix(&alphaBetaStack,nStates,maxT);	/*hmem*/
   betaStar = CreateDMatrix(&alphaBetaStack,nStates,maxT);

   exptimeStart =  CreateDMatrix(&alphaBetaStack,nStates-1,maxT);	/*hmem*/
   exptimeEnd = CreateDMatrix(&alphaBetaStack,nStates-1,maxT);	/*hmem*/
   gamma = CreateDMatrix(&alphaBetaStack,nStates-1,maxT);

   durprob = CreateMatrix(&alphaBetaStack,nStates-1,maxD+1);

   outprob = CreateMatrix(&alphaBetaStack,nStates-1,maxT); /* row 1 not used */
   ZeroMatrix(outprob);		/*hmath*/
   if (maxMixes>1){
      mixoutp = (Matrix**)New(&alphaBetaStack, (nStates-2)*sizeof(Matrix*));	/*hmem*/
	  /*gammaMix = (Matrix**)New(&alphaBetaStack, (nStates-2)*sizeof(Matrix*));
	  gammaMix -= 2;*/
      mixoutp -= 2;
      for (j=2;j<nStates;j++){
         mixoutp[j] = (Matrix*)New(&alphaBetaStack, maxT*sizeof(Matrix));
		 /*gammaMix[j] = (Matrix*)New(&alphaBetaStack, maxT*sizeof(Matrix));
		 --gammaMix[j];*/
         --mixoutp[j];
         for (t=1;t<=maxT;t++){
            mixoutp[j][t] = CreateMatrix(&alphaBetaStack,nStreams,maxMixes);
			/*gammaMix[j][t] = CreateMatrix(&alphaBetaStack,nStreams,maxMixes);*/
            for (s=1;s<=nStreams;s++){
				for (m=1;m<=maxMixes;m++) {
                  mixoutp[j][t][s][m]=LZERO;
				  /*gammaMix[j][t][s][m]=LZERO;*/
				}
            }
         }
      }
   }
   if (nStreams>1){
      stroutp = (Vector**)New(&alphaBetaStack, maxT*sizeof(Vector*));
      --stroutp;
      for (t=1;t<=maxT;t++){
         stroutp[t] = (Vector*)New(&alphaBetaStack,(nStates-2)*sizeof(Vector));
         stroutp[t] -= 2;
         for (j=2;j<nStates;j++)
            stroutp[t][j] = CreateVector(&alphaBetaStack,nStreams);	/*hmem*/
      }
   }
   occr = CreateVector(&gstack,nStates-1);
   zot = CreateVector(&gstack,hset.vecSize);

   numZeroEnd_alpha=CreateVector(&gstack,nStates);
   numZeroBegin_betaS=CreateVector(&gstack,nStates);
}

/* ---------------------------- Load Data ------------------------- */


/* CheckData: check data file consistent with HMM definition */
void CheckData(char *fn, BufferInfo info) 
{
   char tpk[80];
   char mpk[80];
   
   if (info.tgtPK != hset.pkind)
      HError(2250,"CheckData: Parameterisation in %s[%s] is incompatible with hmm %s[%s]",
             fn,ParmKind2Str(info.tgtPK,tpk),hmmfn,ParmKind2Str(hset.pkind,mpk));	/*hparm*/

   if (info.tgtVecSize!=hset.vecSize)
      HError(2250,"CheckData: Vector size in %s[%d] is incompatible with hmm %s[%d]",
             fn,info.tgtVecSize,hmmfn,hset.vecSize);
}

/* InitSegStore : Initialise segStore for particular observation */
void InitSegStore(BufferInfo *info)
{
   Observation obs;
   Boolean eSep;

   SetStreamWidths(info->tgtPK,info->tgtVecSize,hset.swidth,&eSep);	/*hparm*/
   obs = MakeObservation(&gstack,hset.swidth,info->tgtPK,	/*hparm*/
                         hset.hsKind==DISCRETEHS,eSep);
   segStore = CreateSegStore(&segmentStack,obs,10);	/*htrain*/
   firstTime = FALSE;
}

/* LoadFile: load whole file or segments into segStore */
void LoadFile(char *fn)
{
   BufferInfo info;		/*hparm*/
   char labfn[80];
   Transcription *trans;	/*hlabel*/
   long segStIdx,segEnIdx;
   static int segIdx=1;  /* Between call handle on latest seg in segStore */  
   static int prevSegIdx=1;
   HTime tStart, tEnd;		/*typedef double HTime*/
   int k,i,s,ncas,nObs,segLen;
   LLink p;		/*hlabel*/
   Observation obs;		/*hparm*/

   if((pbuf=OpenBuffer(&bufferStack, fn, 10, dff, FALSE_dup, FALSE_dup))==NULL)	/*hparm*/
      HError(2250,"LoadFile: Config parameters invalid");
   GetBufferInfo(pbuf,&info);	/*hparm*/
   CheckData(fn,info);		/*hrest*/
   if (firstTime) InitSegStore(&info);		/*hrest*/

   if (segId == NULL)  {   /* load whole parameter file */
      nObs = ObsInBuffer(pbuf);		/*hparm*/
      tStart = 0.0;
      tEnd = (info.tgtSampRate * nObs);
      LoadSegment(segStore, tStart, tEnd, pbuf);	/*htrain*/
      if (nObs > maxT) 
         maxT=nObs; 
      if (nObs < minT)
         minT=nObs;      
      segIdx++;
   }
   else {                  /* load segment of parameter file */
      MakeFN(fn,labDir,labExt,labfn);	/*hshell*/
      trans = LOpen(&transStack,labfn,lff);		/*hlabel*/
      ncas = NumCases(trans->head,segId);      /*hlabel*/
      nObs = 0;
      if ( ncas > 0) {
         for (i=1,nObs=0; i<=ncas; i++) {
            p = GetCase(trans->head,segId,i);	/*hlabel*/
            segStIdx= (long) (p->start/info.tgtSampRate);
            segEnIdx  = (long) (p->end/info.tgtSampRate);
            if (segEnIdx >= ObsInBuffer(pbuf))		/*hparm*/
               segEnIdx = ObsInBuffer(pbuf)-1;		/*hparm*/
            if (((segEnIdx - segStIdx + 1 >= nStates-2) || !segReject) 
		&& (segStIdx <= segEnIdx)) {	/* skip short segments */
               LoadSegment(segStore, p->start, p->end, pbuf);	/*htrain*/
               if (trace&T_LD1)
                  printf("  loading seg %s %f[%ld]->%f[%ld]\n",segId->name,
                         p->start,segStIdx,p->end,segEnIdx);
               segLen = SegLength(segStore, segIdx);	/*htrain*/
               nObs += segLen;
               if (segLen > maxT) 
                  maxT=segLen; 
               if (segLen < minT)
                  minT=segLen;
               segIdx++;
            }else if (trace&T_LD1)
               printf("   seg %s %f->%f ignored\n",segId->name,
                      p->start,p->end);
         }        
      }   
   }
   if (hset.hsKind == DISCRETEHS){
      for (k=prevSegIdx; k<segIdx; k++){
         segLen = SegLength(segStore, k);	/*htrain*/
         for (i=1; i<=segLen; i++){
            obs = GetSegObs(segStore, k, i);	/*htrain*/
            for (s=1; s<=nStreams; s++){
               if( (obs.vq[s] < 1) || (obs.vq[s] > maxMixInS[s]))
                  HError(2250,"LoadFile: Discrete data value [ %d ] out of range in stream [ %d ] in file %s",obs.vq[s],s,fn);
            }
         }
      }
      prevSegIdx=segIdx;
   }

   if (trace&T_LD0)
      printf(" %d observations loaded from %s\n",nObs,fn);
   CloseBuffer(pbuf);		/*hparm*/
   ResetHeap(&transStack);		/*hmem*/
}


/* ------------------------ Trace Functions -------------------- */

/* ShowSegNum: if not already printed, print seg number */
void ShowSegNum(int seg)
{
   static int lastseg = -1;
   
   if (seg != lastseg){
      printf("---- Training Segment %d [%3d frames] ----\n",seg,T);
      lastseg = seg;
   }
}
   
/* ------------------------- Alpha-Beta ------------------------ */

/* SetOutP: Set the output and mix prob matrices */                        
void SetOutP(int seg)	/*compute all b_j(Ot) and N(Ot,mu_jk, U_jk)  */
{
   int i,t,d,m,mx,s,nMix;
   StreamElem *se;	/*hmodel*/
   MixtureElem *me;		/*hmodel*/
   StateInfo *si;	/*hmodel*/
   Matrix mixp;		/*typedef float **Matrix */
   LogFloat x,prob,streamP;		/*typedef float LogFloat */
   Vector strp = NULL;	/*typedef float *vector */
   Observation obs;		/*hparm*/
   TMixRec *tmRec = NULL;	/*hmodel*/
   float wght,tmp;
   MixPDF *mpdf;		/*hmodel*/
   PreComp *pMix;	/*htrain*/
   
   for (t=1;t<=T;t++) {
      obs = GetSegObs(segStore, seg, t);	/*htrain*/
      if (hsKind == TIEDHS)
         PrecomputeTMix(&hset,&obs,tMPruneThresh,0);	/*hmodel*/       
      if ((maxMixes>1) && (hsKind!=DISCRETEHS)){ /* Multiple Mix Case */
         for (i=2;i<nStates;i++) {
            prob = 0.0;
            si = hmm->svec[i].info;
            se = si->pdf+1; 
            mixp = mixoutp[i][t];
            if (nStreams>1) strp = stroutp[t][i];
            for (s=1;s<=nStreams;s++,se++){
               switch (hsKind){         /* Get nMix */
               case TIEDHS:
                  tmRec = &(hset.tmRecs[s]);
                  nMix = tmRec->nMix;
                  break;
               case PLAINHS:
               case SHAREDHS:
                  nMix = se->nMix;
                  break;
               } 
               streamP = LZERO;
               for (mx=1;mx<=nMix;mx++) {
                  m=(hsKind==TIEDHS)?tmRec->probs[mx].index:mx;
                  switch (hsKind){      /* Get wght and mpdf */
                  case TIEDHS:
                     wght=se->spdf.tpdf[m];
                     mpdf=tmRec->mixes[m];
                     break;
                  case PLAINHS:
                  case SHAREDHS:
                     me = se->spdf.cpdf+m;
                     wght=me->weight;
                     mpdf=me->mpdf;
                     break;
                  }
                  if (wght>MINMIX){
                     switch(hsKind) { /* Get mixture prob */
                     case TIEDHS:
                        tmp = tmRec->probs[mx].prob;
                        x = (tmp>=MINLARG)?log(tmp)+tmRec->maxP:LZERO;
                        break;
                     case SHAREDHS : 
                        pMix = (PreComp *)mpdf->hook;
                        if (pMix->time==t)
                           x = pMix->prob;
                        else {
                           x = MOutP(obs.fv[s],mpdf);	/*hmodel*/
                           pMix->prob = x; pMix->time = t;
                        }
                        break;
                     case PLAINHS : 
                        x=MOutP(obs.fv[s],mpdf);	/*hmodel*/
                        break;
                     default:
                        x=LZERO;
                        break;
                     }
                     mixp[s][m]=x;
                     streamP = LAdd(streamP,log(wght)+x);	/*hmath*/
                  } else
                     mixp[s][m]=LZERO;
               }               
               if (nStreams>1)
                  strp[s]=streamP;
               prob += streamP; /* note stream weights ignored */
            }   
            outprob[i][t]=prob;
         }
      } else 
         if (nStreams>1) {      /* Single Mixture multiple stream */
            for (i=2;i<nStates;i++) {
               prob = 0.0;
               si = hmm->svec[i].info;
               se = si->pdf+1;
               strp = stroutp[t][i];
               for (s=1;s<=nStreams;s++,se++){
                  streamP = SOutP(&hset,s,&obs,se);	/*hmodel*/
                  strp[s] = streamP;
                  prob += streamP; /* note stream weights ignored */
               }
               outprob[i][t]=prob;
            }
         } else                 /* Single Mixture - Single Stream */
            for (i=2;i<nStates;i++){
               si = hmm->svec[i].info;
               se = si->pdf+1;
               if (hsKind==DISCRETEHS)
                  outprob[i][t]=SOutP(&hset,1,&obs,se);		/*hmodel*/
               else
                  outprob[i][t]=OutP(&obs,hmm,i);		/*hmodel*/
			   /*outprob[i][t]=-80;*/
            }
   }

   /*ShowMatrix("OutProb",outprob,T,nStates-1);	*//*hmath*/
   if (trace  & T_OTP) {
      ShowSegNum(seg);
      ShowMatrix("OutProb",outprob,10,12);	/*hmath*/
   }
}

/*
void SetDurP()
{
   StreamElem *se;	//hmodel
   StateInfo *si;	//hmodel
   PoissonPDF *ppdf;
   int i,d;

//static Matrix durprob;	    array[2..nStates-1][1..maxD] of duration prob
   for (i=2;i<nStates;i++) {
        si = hmm->svec[i].info;
		ppdf=si->ppdf;
		printf("%f\n",ppdf->mu);
		for (d=1;d<=maxD;d++)
			durprob[i][d]=PoissonOut(d,ppdf);	//no normalization, hmodel
   }

}
*/

void SetDurP()
{

   StateInfo *si;	/*hmodel*/
   int i,d;

/*static Matrix durprob;	    //array[2..nStates-1][1..maxD] of duration prob*/

   for (i=2;i<nStates;i++) {
		si=hmm->svec[i].info;

		printf("state%d:\n",i);
		for (d=1;d<=si->maxD;d++){
			durprob[i][d]=si->dur[d];	/*no normalization, hmodel*/
			printf("%.3f ",durprob[i][d]);
		}
		printf("\n");
   }

}

/* SetAlpha: compute alpha matrix and return prob of given sequence */
/* the forward backward process */
LogDouble SetAlphaBeta(int seg)
{
   int i,j,k,t,d;
   LogDouble a;
   LogDouble probtran, outstay;
   LogDouble outAccprob;
   LogDouble POa, POb;
   Vector durprobj;
   StateInfo *si;

   /*Compute alpha and alphaStar*/
 
   alphaStar[1][0]=LZERO;
   alpha[1][0]=0.0;
   for(j=2;j<nStates;j++){
	   alphaStar[j][0]=hmm->transP[1][j];
	   alpha[j][0]=LZERO;
   }

   outAccprob=0.0;

   ZeroVector(numZeroEnd_alpha);

   for(i=2;i<nStates-1;i++)
	   if(hmm->transP[i][nStates]<LSMALL)
		   for(j=nStates-2;j>=i;j--)
			   numZeroEnd_alpha[i]+=(hmm->transP[j][nStates]<LSMALL);

   /* 1<=t<= min(T,maxD) */
   for(t=1;t<=T;t++) {
	   alphaStar[1][t]=LZERO; alphaStar[nStates][t]=LZERO;
	   alpha[1][t]=LZERO; alpha[nStates][t]=LZERO;
	   alphaStar[2][t]=LZERO;
	   for (j=2;j<nStates-1;j++) { 
		   si=hmm->svec[j].info;
		   durprobj=si->dur;
		   maxD=si->maxD;
 	  
		   if (j==2){	/*state 2	*/	   
			  outAccprob+=outprob[2][t];	
			  if ((t<=maxD)&&(t<=T-numZeroEnd_alpha[2])) {
				  if(durprobj!=NULL)
					alpha[2][t]=alphaStar[2][0]+durprobj[t]+outAccprob;
				  else
					alpha[2][t]=alphaStar[2][0]+0.0+outAccprob;
			  }
			  else
				  alpha[2][t]=LZERO;

			  a=hmm->transP[2][3];
			  alphaStar[3][t]=alpha[2][t]+a;
		   }
		   else {  /*for state 3 and above*/
			   if(t<=T-numZeroEnd_alpha[j]){
				  probtran=LZERO;
				  outstay=0.0;
				  for (d=1;d<=IntMin(t,maxD);d++){
					  outstay+=outprob[j][t-d+1];	/* b_j(O_t-d+1,...,O_t) */
					  if(durprob!=NULL)
						  probtran=LAdd(probtran,alphaStar[j][t-d]+durprobj[d]+outstay);
					  else
						  probtran=LAdd(probtran,alphaStar[j][t-d]+0.0+outstay);
				  }
				  alpha[j][t]=probtran;	
			   }
			   else
				  alpha[j][t]=LZERO;

			   alphaStar[j+1][t]=LZERO;
			   for(k=2;k<=j;k++){
				  a=hmm->transP[k][j+1];
				  if (a>LSMALL)
					  alphaStar[j+1][t]=LAdd(alphaStar[j+1][t],alpha[k][t]+a);
			   }
			   if (alphaStar[j+1][t]<LZERO)
				  alphaStar[j+1][t]=LZERO;
		   }

	   }
	   alpha[nStates-1][t]=LZERO;
   }
  
	/*Compute PO, at t=T*/

	probtran=LZERO;
	outstay=0.0;
	si=hmm->svec[nStates-1].info;
	durprobj=si->dur;
	maxD=si->maxD;
	for (d=1;d<=IntMin(T,maxD);d++){
		outstay+=outprob[nStates-1][T-d+1];	/* b_j(O_t-d+1,...,O_t) */
		if(durprobj!=NULL)
			probtran=LAdd(probtran,alphaStar[nStates-1][T-d]+durprobj[d]+outstay);
		else
			probtran=LAdd(probtran,alphaStar[nStates-1][T-d]+0.0+outstay);
	}

	alpha[nStates-1][T]=probtran;

	alphaStar[nStates-1][T]=LZERO;
	alpha[nStates][T]=0.0;
	alphaStar[nStates][T]=0.0;
/*
	for(i=2;i<nStates-1;i++)
		if(hmm->transP[i][nStates]<LSMALL)
			alpha[i][T]=LZERO; 
*/
	POa=LZERO;
	for(i=2;i<nStates;i++)
		if(alpha[i][T]>LSMALL)
			POa=LAdd(POa,alpha[i][T]);

/*
	printf("alpha:\n");
	for(i=2;i<nStates;i++){
		for(t=0;t<=T;t++)
			printf("%.3e ",alpha[i][t]);
		printf("\n");
	}
	printf("alpha*:\n");
	for(i=2;i<nStates;i++){
		for(t=0;t<=T;t++)
			printf("%.3e ",alphaStar[i][t]);
		printf("\n");
	}
*/

	/*Compute beta and betaStar*/
   beta[1][T]=LZERO; betaStar[nStates][T]=0.0; 
   
   for(j=2;j<nStates;j++){
	   beta[j][T]=hmm->transP[j][nStates];
	   betaStar[j][T]=LZERO;
   }


   ZeroVector(numZeroBegin_betaS);
   for(i=nStates-1;i>2;i--)
	   if(hmm->transP[1][i]<LSMALL)
		   for(j=2;j<=i;j++)
			   numZeroBegin_betaS[i]+=(hmm->transP[1][j]<LSMALL);

   outAccprob=0.0;
   /* 1<=t<= min(T,maxD)*/ 
   for(t=T-1;t>=0;t--) {
	   beta[1][t]=LZERO; beta[nStates][t]=LZERO;
	   betaStar[1][t]=LZERO; betaStar[nStates][t]=LZERO;
	   beta[nStates-1][t]=LZERO;
	   for (j=nStates-1;j>2;j--) { 
		   si=hmm->svec[j].info;
		   durprobj=si->dur;
		   maxD=si->maxD;

		   if(j==nStates-1){	/*state 4*/

			  outAccprob+=outprob[j][t+1];	
			  if((T-t<=maxD)&&(t>=numZeroBegin_betaS[nStates-1])) {
				  if(durprobj!=NULL)
					  betaStar[nStates-1][t]=durprobj[T-t]+outAccprob+beta[nStates-1][T];
				  else
					  betaStar[nStates-1][t]=0.0+outAccprob+beta[nStates-1][T];
			  }
			  else
				  betaStar[nStates-1][t]=LZERO;

			  a=hmm->transP[nStates-2][nStates-1];
			  beta[nStates-2][t]=a+betaStar[nStates-1][t];
		   }
		   else {  /*for state 3 and above*/
			   if(t>=numZeroBegin_betaS[j]) {
				  probtran=LZERO;
				  outstay=0.0;
				  for (d=1;d<=IntMin(T-t,maxD);d++){
					  outstay+=outprob[j][t+d];	/* b_j(O_t-d+1,...,O_t) */
					  if(durprobj!=NULL)
						  probtran=LAdd(probtran,durprobj[d]+outstay+beta[j][t+d]);
					  else
						  probtran=LAdd(probtran,0.0+outstay+beta[j][t+d]);
				  }
				  betaStar[j][t]=probtran;			  	
			   }
			   else
				   betaStar[j][t]=LZERO;

			  beta[j-1][t]=LZERO;
			  if(t>0){
				  for(k=j;k<=nStates-1;k++){
					  a=hmm->transP[j-1][k];
					  if (a>LSMALL)
						  beta[j-1][t]=LAdd(beta[j-1][t],a+betaStar[k][t]);
				  }
				  if(beta[j-1][t]<LZERO)
					  beta[j-1][t]=LZERO;
			  }
		   }
	   }
	   betaStar[2][t]=LZERO;	/*state 2 can only start at 0*/
   }

	/*Compute PO, at t=0*/

    probtran=LZERO;
    outstay=0.0;
	si=hmm->svec[2].info;
	durprobj=si->dur;
	maxD=si->maxD;
    for (d=1;d<=IntMin(T,maxD);d++){
  		outstay+=outprob[2][d];	/* b_j(O_t-d+1,...,O_t) */
		if(durprobj!=NULL)
			probtran=LAdd(probtran,durprobj[d]+outstay+beta[2][d]);
		else
			probtran=LAdd(probtran,0.0+outstay+beta[2][d]);
	}
	betaStar[2][0]=probtran;	
	beta[1][0]=LZERO;
	beta[nStates][0]=LZERO;

	/*betaStar[3][0]=LZERO;betaStar[4][0]=LZERO;beta[2][0]=LZERO;*/

	POb=LZERO;
	for(i=2;i<nStates;i++)
		POb=LAdd(POb,betaStar[i][0]);

/*	POb=betaStar[2][0];*/

/*
   if (trace  & T_ALF) {
      ShowSegNum(seg);
      ShowDMatrix("Alpha",alpha,10,12);		
      ShowDMatrix("AlphaStar",alphaStar,10,12);		//hmath
      printf("LogP= %10.3f\n\n",alpha[nStates][T]);
   }
   */
/*
	printf("beta:\n");
	for(i=2;i<nStates;i++){
		for(t=0;t<=T;t++)
			printf("%.3e ",beta[i][t]);
		printf("\n");
	}

	printf("beta*:\n");
	for(i=2;i<nStates;i++){
		for(t=0;t<=T;t++)
			printf("%.3e ",betaStar[i][t]);
		printf("\n");
	}
*/

    return (POa+POb)/2;
	/*return POa;*/
}



/*
LogDouble SetAlpha(int seg)
{
   int i,j,t,d;
   LogDouble a;
   LogDouble probstay,probtran, outstay;
   Vector outAccprob;

   //alpha initialization
   //alpha[i][t] = prob(O1...Ot & i ends at t); 
   //alphaStar[i][t] = prob(O1...Ot & i begins at t+1); 

   // t=0 
   alpha[1][0] = 0.0;	// the prob(state 1 terminate at t=0)=1 
   for(i=2;i<=nStates;i++)
	   alpha[i][0]=LZERO;

   //t=1
   alpha[1][1] = LZERO;   //prob{state 1 terminates at t=1}=0 
   outAccprob=CreateVector(&gstack,nStates);	// = b_j(O1,O2,...,Ot)
   ZeroVector(outAccprob);
   for (j=2;j<nStates;j++) {					// col 1 from entry state 

	  a=hmm->transP[1][j];
	  outAccprob[j]+=outprob[j][1];

      if (a<LSMALL)
         alpha[j][1] = LZERO;
      else {  
		 alpha[j][1]=a+durprob[j][1]+outAccprob[j];
	  }
   }
   alpha[nStates][1] = LZERO;  

   // 1<t< min(T,maxD) 
   for(t=2;t<=IntMin(T,maxD);t++) {
	   alpha[1][t]=alpha[nStates][t]=LZERO;	//prob{being in state[1], non-emitting at t=2}=0
	   for (j=2;j<nStates;j++) {            // col 1 from entry state 
		  a=hmm->transP[1][j];
		  outAccprob[j]+=outprob[j][t];
		  if (a>LSMALL)
			 probstay=a+durprob[j][t]+outAccprob[j];
		  else
			 probstay=LZERO;
		  probtran=LZERO;
		  outstay=0.0;
		  for (d=1;d<t;d++){
			  outstay+=outprob[j][t-d+1];	// b_j(O_t-d+1,...,O_t) 
			  for(i=2;i<nStates;i++) {
					a=hmm->transP[i][j];	
					if (a>LSMALL)
						probtran=LAdd(probtran,alpha[i][t-d]+a+durprob[j][d]+outstay);
			  }
		  }
		  alpha[j][t]=LAdd(probstay,probtran);
	   }
   }
  
   if(T>maxD)
   {
	//alpha induction
	   for (t=maxD+1;t<=T;t++) {            // cols 2 to T 
		  for (j=2;j<nStates;j++) {
			 probtran=LZERO;
			 outstay=0.0;
			 for(d=1;d<=maxD;d++) {
				outstay+=outprob[j][t-d+1];	// b_j(O_t-d+1,...,O_t) 
				for (i=2;i<nStates;i++) {
					a=hmm->transP[i][j];
					if (a>LSMALL)
						probtran=LAdd(probtran,alpha[i][t-d]+a+durprob[j][d]+outstay);
				}
			 }
		     alpha[j][t]=probtran;			// alpha inductive equation 
		  }
		  alpha[1][t] = alpha[nStates][t] = LZERO;
	   }
   }

   // finally calc seg prob, nStates terminate at T+1 
   probtran = LZERO ;                      
   for (i=2;i<nStates;i++) {
      a=hmm->transP[i][nStates];
      if (a>LSMALL)
         probtran=LAdd(probtran,alpha[i][T]+a);	
   }  

   //Calculate alphaStar from alpha
   for (t=0;t<T;t++) { 
 	   for(j=1;j<=nStates;j++) {
		   alphaStar[j][t]=LZERO;
		   for(i=1;i<=nStates;i++) {
			   a=hmm->transP[i][j];
			   if (a>LSMALL)
				  alphaStar[j][t]=LAdd(alphaStar[j][t],alpha[i][t]+a);
		   }
	   }
	   //alphaStar[1][t]=alphaStar[nStates][t]=LZERO;
   }

   for(j=1;j<=nStates;j++) 
	   alphaStar[j][T]=LZERO;


   for(i=1;i<=nStates;i++) {
	    a=hmm->transP[i][nStates];
		if (a>LSMALL)
			alphaStar[nStates][T]=LAdd(alphaStar[nStates][T],alpha[i][T]+a);
   }


   if (trace  & T_ALF) {
      ShowSegNum(seg);
      ShowDMatrix("Alpha",alpha,10,12);		
      ShowDMatrix("AlphaStar",alphaStar,10,12);		//hmath
      printf("LogP= %10.3f\n\n",alpha[nStates][T]);
   }
   FreeVector(&gstack,outAccprob);
   return probtran;
}

*/

/* SetBeta: compute beta matrix */
/* the backward process */
/*
LogDouble SetBeta(int seg)
{
   int i,j,t,d;
   LogDouble a;
   LogDouble probstay,probtran, outstay;
   Vector outAccprob;

	//betaStar initialization
   // betaStar[i][t]=Prob{Ot...O_T | i begins at t+1} 
   // beta=Prob{Ot...O_T | i terminates at t} 

   betaStar[nStates][T] = 0.0;   // prob(nStates begins at T+1)=1 
   for (i=1;i<nStates;i++)
	   betaStar[i][T]=LZERO;

   // t=T-1 
   outAccprob=CreateVector(&gstack,nStates);	// = b_j(O1,O2,...,Ot)
   ZeroVector(outAccprob);
   for (i=2;i<nStates;i++) {               
	   a=hmm->transP[i][nStates];
	   outAccprob[i]+=outprob[i][T];
	   if (a>LSMALL)
		   betaStar[i][T-1]=durprob[i][1]+outAccprob[i]+a;
	   else
		   betaStar[i][T-1]=LZERO;
   }
   betaStar[1][T-1] = LZERO;
   betaStar[nStates][T-1]=LZERO;

   // max(0,T-maxD)<t<T 
   for(t=T-2;t>=IntMax(0,T-maxD);t--) {
	   betaStar[1][t]=betaStar[nStates][t]=LZERO;	//prob{being in non-emitting states at t}=0
	   for (i=2;i<nStates;i++) {            //col 1 from entry state
		  a=hmm->transP[i][nStates];
		  outAccprob[i]+=outprob[i][t+1];
		  if (a>LSMALL)
			 probstay=durprob[i][T-t]+outAccprob[i]+a;
		  else
			 probstay=LZERO;
		  probtran=LZERO;
		  outstay=0.0;
		  for (d=1;d<=T-t-1;d++){
			  outstay+=outprob[i][t+d];	// b_j(O_t+1,...,O_t+d) 
			  for(j=2;j<nStates;j++) {
					a=hmm->transP[i][j];	
					if (a>LSMALL)
						probtran=LAdd(probtran,durprob[j][d]+outstay+a+betaStar[j][t+d]);
			  }
		  }
		  betaStar[i][t]=LAdd(probstay,probtran);
	   }
   }

   if(T-maxD>0)
   {
	//betaStar induction
	   for (t=T-maxD-1;t>=0;t--) {            // cols 2 to T 
		  for (i=2;i<nStates;i++) {
			 probtran=LZERO;
			 outstay=0.0;
			 for(d=1;d<=maxD;d++) {
				outstay+=outprob[i][t+d];	// b_j(O_t+1,...,O_t+d) 
				for (j=2;j<nStates;j++) {
					a=hmm->transP[i][j];
					if (a>LSMALL)
						probtran=LAdd(probtran,durprob[j][d]+outstay+a+betaStar[j][t+d]);
				}
			 }
		     betaStar[j][t]=probtran;			// alpha inductive equation 
		  }
		  betaStar[1][t] = betaStar[nStates][t] = LZERO;
	   }
   }

   // finally calc seg prob, t=-1 

   probtran = LZERO ;                      
   for (j=2;j<nStates;j++) {
      a=hmm->transP[1][j];
      if (a>LSMALL)
         probtran=LAdd(probtran,a+betaStar[j][0]);	
   }  


   //Calculate beta from betaStar
   for (t=1;t<=T;t++) { 
 	   for(i=1;i<=nStates;i++) {
		   beta[i][t]=LZERO;
		   for(j=1;j<=nStates;j++) {
			   a=hmm->transP[i][j];
			   if (a>LSMALL)
				  beta[i][t]=LAdd(beta[i][t],a+betaStar[j][t]);
		   }
	   }
   }


   for(i=1;i<=nStates;i++) 
	   beta[i][0]=LZERO;
   for(j=1;j<=nStates;j++) {
	   a=hmm->transP[1][j];
	   if (a>LSMALL)
		  beta[1][0]=LAdd(beta[1][0],a+betaStar[j][0]);
   }


   if (trace & T_BET) {
      ShowSegNum(seg);	//hrest
      ShowDMatrix("BetaStar",betaStar,10,12);	//hmath
      ShowDMatrix("Beta",beta,10,12);	//hmath
      printf("LogP=%10.3f\n\n",betaStar[1][0]);
   }
   FreeVector(&gstack,outAccprob);
   return probtran;
}
*/


/* --------------------- Record Statistics ---------------- */

/* SetOccr: set the global occupation counters occr for current seg */
/* compute occr[i]= expected number of times in state i, given seg */
/* occr[i] is equivalent to sum_t gamma_t(i) */

void SetOccr(LogDouble pr, int seg)
{
   int i,j,t;
   LogDouble x,y;
   /*Vector sum_i;

   //compute exptimeStart, S_t(i) and exptimeEnd E_t(i)

   //sum_i=CreateVector(&gstack,T);
   //ZeroVector(sum_i);*/
   for(i=2;i<nStates;i++) {
	   x=LZERO; y=LZERO;
	   occr[i]=0.0;
	   for(t=1; t<=T; t++) {
			x=LAdd(x,alphaStar[i][t-1]+betaStar[i][t-1]-pr);	 /* S_t(i) */
			y=LAdd(y,alpha[i][t-1]+beta[i][t-1]-pr);		/* E_t(i) */
			if(x>MINEARG)
				exptimeStart[i][t]=exp(x);	
			else
			    exptimeStart[i][t]=0.0;

			if(y>MINEARG)
			    exptimeEnd[i][t]=exp(y);		
			else
			    exptimeEnd[i][t]=0.0;
			
			gamma[i][t]=exptimeStart[i][t]-exptimeEnd[i][t];	
			
			occr[i]+=gamma[i][t];
			/*sum_i[t]+=gamma[i][t];*/
	   }

   }

   occr[1] = 1.0; occr[nStates]=1.0;

   /*FreeVector(&gstack,sum_i);*/

   if (trace & T_OCC){
      ShowSegNum(seg);
      ShowVector("OCC: ",occr,20);	/*hmath*/
   }
	/*ShowVector("occr: ",occr,20);	//hmath*/
}

/*
void SetOccr(LogDouble pr, int seg)
{
   int i,t;
   LogDouble x,y;

   //compute exptimeStart, S_t(i) and exptimeEnd E_t(i)
   for(i=2;i<nStates;i++) {
	   x=0.0; y=0.0;
	   occr[i]=0.0;
	   for(t=1; t<=T; t++) {

		   x+=exp(alphaStar[i][t-1]+betaStar[i][t-1]-pr);	 // S_t(i) 
		   y+=exp(alpha[i][t-1]+beta[i][t-1]-pr);		// E_t(i) 
		   exptimeStart[i][t]=x;	
		   exptimeEnd[i][t]=y;		
		   occr[i]+=exptimeStart[i][t]-exptimeEnd[i][t];
	   }
   }

   occr[1] = 1.0; occr[nStates]=1.0;

   if (trace & T_OCC){
      ShowSegNum(seg);
      ShowVector("OCC: ",occr,20);	//hmath
   }
}

*/

/* UpDurCounts: update duration PDF counts for each stream of each state */
void UpDurCounts(int j, StateInfo *si, LogDouble pr, DVector alphaStarj, DVector betaj)
{
   int d,t;

   DurAcc *da;	
   LogDouble Lr;
   LogDouble outstay;
   double y;
   short dsize;
   Vector durprobj;

   da = (DurAcc *)GetHook(si->dur);/*
   //printf("state%d:\n",j);
   //for(d=1;d<=15;d++)
   //printf("%.4f ",da->dur[d]);
   //printf("\n");
   //ShowVector("da->dur",da->dur,maxD);*/

   durprobj=si->dur;
   maxD=si->maxD;
   for (t=0; t<T; t++) {
	  outstay=0.0;
      Lr = LZERO;
	  for(d=1;d<=IntMin(T-t,maxD);d++) {
		  outstay+=outprob[j][t+d];
		  /*Lr=LAdd(Lr,alphaStarj[t]+durprob[j][d]+outstay+betaj[t+d]-pr);*/
		  Lr=alphaStarj[t]+durprobj[d]+outstay+betaj[t+d]-pr;
          if (Lr > MINEARG) {
             y = exp(Lr);
             if (uFlags&UPDURS) { 
				 da->dur[d]+=y;
				 da->occ+=y;
			 }
		  }

		  /*printf("%.4f ",da->dur[d]);*/
	  }
	  /*//printf("\n");
	  //ShowVector("da->dur",da->dur,maxD);*/
   } 
}

/*
void UpDurCounts(int j, StateInfo *si)	
{
   int tao,d,t;

   DurAcc *da;	
   LogDouble Lr;
   LogDouble outstay;
   double y;
//typedef struct {     
// float exptldur;       
//  float exptlnumbegin;        
//} DurAcc;

   da = (DurAcc *)si->hook;	
   da->exptldur+=occr[j];
   da->exptlnumbegin+=exptimeStart[j][T];
}

*/

/* UpTranCounts: update the transition counters in ta */
void UpTranCounts(LogDouble pr,int seg)
{
   int i,j,t;
   Matrix tran;
   Vector tran_i,outprob_j,a_i,occ;
   DVector alpha_i,betaStar_j;
   LogDouble x,a_ij;
   double y;
   TrAcc *ta;	/*htrain*/
   
   ta = (TrAcc *) GetHook(hmm->transP);	/*hmem*/
   tran = ta->tran; occ = ta->occ;
   for (i=2; i<nStates; i++)
      occ[i] += occr[i];	/*occ=occr accumulated_over_seg, occr[i]=expected duration in state i*/

   tran_i = tran[1];          /* transitions 1->j    1<j<nStates */
   a_i = hmm->transP[1];
   for (j=2;j<nStates;j++) {
      a_ij = a_i[j];
      if (a_ij>LSMALL) {
         x = alpha[1][0]+a_ij + betaStar[j][0] - pr;
         if (x>MINEARG) {
            y =  exp(x);
            tran_i[j] += y; occ[1] += y;	/*occ[1]=sum_j number of trans1->j, used to compute prior */
         }
      }
   }
   for (i=2;i<nStates;i++) {        /* transitions i->j    1<i,j<nStates */
      a_i = hmm->transP[i];
      alpha_i = alpha[i];
      tran_i = tran[i];
      for (j=2;j<nStates;j++) {
         a_ij=a_i[j];
         if (a_ij>LSMALL) {
            x=LZERO; betaStar_j=betaStar[j]; 
            for (t=1;t<=T-1;t++)
               x=LAdd(x,alpha_i[t]+a_ij+betaStar_j[t]-pr);	/*hmath*/
            if (x>MINEARG)
               tran_i[j] += exp(x);	/*expected number of times of transition i->j, C(i,j) in Ferguson */
         }
      }
   }
   for (i=2; i<nStates; i++) {    /* transitions i->nStates    1<i<nStates */
      a_ij = hmm->transP[i][nStates];
      if (a_ij>LSMALL) {
         x = alpha[i][T] + a_ij + betaStar[nStates][T]- pr;
         if (x>MINEARG)
            tran[i][nStates] += exp(x);
      }     
   }
   if (trace & T_TAC){
      ShowSegNum(seg);
      ShowMatrix("TRAN: ",tran,10,10);	/*hmath*/
      ShowVector("TOCC: ",occ,10);		/*hmath*/
      fflush(stdout);
   }
}


/* UpStreamCounts: update mean, cov & mixweight counts for given stream */
void UpStreamCounts(int j, int s, StreamElem *se, int vSize, LogDouble pr, int seg,
                    DVector alphaStarj, DVector betaj)
{
   int i,m,nMix,k,l,t,ss,idx;
   MixtureElem *me;	/*hmodel*/
   MixPDF *mpdf;	/*hmodel*/
   MuAcc *ma;	/*htrain*/
   WtAcc *wa;	/*hmodel*/
   VaAcc *va;	/*hmodel*/
   Matrix *mixp_j;
   Vector ot, strpt;
   LogFloat a_ij,w;
   LogDouble Lr;
   double y;
   Observation obs;		/*hparm*/
   TMixRec *tmRec = NULL;	/*hmodel*/
   float wght;
   
   wa = (WtAcc *)se->hook;
   switch (hsKind){       /* Get nMix */
   case TIEDHS:
      tmRec = &(hset.tmRecs[s]);
      nMix = tmRec->nMix;
      break;
   case PLAINHS:
   case SHAREDHS:
      nMix = se->nMix;
      break;      
   case DISCRETEHS:
      nMix = 1;                /* Only one code selected per observation */
      break;
   }
   mixp_j = (maxMixes>1) ? mixoutp[j] : NULL;
   for (m=1; m<=nMix; m++) {
      switch (hsKind){            /* Get mpdf, wght */
      case TIEDHS:               
         wght=se->spdf.tpdf[m];
         mpdf=tmRec->mixes[m];
         break;
      case DISCRETEHS:
         wght=1.0;                 /* weight for DISCRETEHS has to be 1 */
         mpdf=NULL;
         break;
      case PLAINHS:
      case SHAREDHS:
         me = se->spdf.cpdf+m;
         wght=me->weight;
         mpdf=me->mpdf;
         break;
      }
      if (hsKind!=DISCRETEHS){
         ma = (MuAcc *)GetHook(mpdf->mean);		/*hmem*/
         va = (VaAcc *)GetHook(mpdf->cov.var);	/*hmem*/
      }
      if (wght > MINMIX) {
         w = log(wght);
         for (t=1; t<=T; t++) {
         
            /* Get observation vec ot and zero mean zot */
            obs = GetSegObs(segStore, seg, t);	/*htrain*/
            ot = obs.fv[s];
            if (hsKind!=DISCRETEHS)
               for (k=1; k<=vSize; k++)
                  zot[k] = ot[k] - mpdf->mean[k];	/*zero mean observation*/
                                 
            /* Compute state/mix occupation log likelihood */
			 y=0;
            if (nMix==1 || (hsKind==DISCRETEHS)){
               /*Lr = log(gamma[j][t]);*/
				y=gamma[j][t];
            }
            else {
			   /*
               Lr = alphaStar[j][t-1];

               if (Lr>LSMALL) {
                  Lr += mixp_j[t][s][m] + w + betaj[t] - pr;
				//Lr, the residence prob in mixture m of state j stream s at time t, gamma_t[j][m] in stream s
                  if (nStreams>1) { // add contrib of parallel streams 
                     strpt = stroutp[t][j];
                     for (ss=1; ss<=nStreams; ss++)
                        if (ss!=s) Lr += strpt[ss];
                  }
               }*/
			   Lr = log(gamma[j][t]);
			   if(Lr>LSMALL) {
				   Lr += mixp_j[t][s][m] + w - outprob[j][t];
                   if (nStreams>1) { /* add contrib of parallel streams */
                       strpt = stroutp[t][j];
                       for (ss=1; ss<=nStreams; ss++)
                          if (ss!=s) Lr += strpt[ss];
                   }
			   }

			   if (Lr > MINEARG) 
				  y = exp(Lr);
            }
            
            if (y > 0) {
               /* Update Weight Counter */
               if (uFlags&UPMIXES) {   
                  idx = (hsKind==DISCRETEHS) ? obs.vq[s] : m;
                  wa->occ += y; wa->c[idx] += y;
               }
               
               /* Update Mean Counter */
               if (uFlags&UPMEANS){
                  ma->occ += y; 
                  for (k=1; k<=vSize; k++)
                     ma->mu[k] += zot[k]*y; /* sum zero mean */
               }
               
               /* Update Covariance Counter */
               if (uFlags&UPVARS){
                  va->occ += y;
                  if (mpdf->ckind==DIAGC){
                     for (k=1; k<=vSize; k++)
                        va->cov.var[k] += zot[k] * zot[k] * y;
                  } else{
                     for (k=1; k<=vSize; k++)
                        for (l=1; l<=k; l++)
                           va->cov.inv[k][l] += zot[k] * zot[l] * y;
                  }
               }  
            }
         }
      }
      if ((trace&(T_MAC|T_VAC))&&(hsKind!=DISCRETEHS)) {
         ShowSegNum(seg);
         printf("State %d, Stream %d, Mixture %d\n",j,s,m);
         if (trace&T_MAC){
            printf("MEAN OCC: %.2f\n",ma->occ);
            ShowVector("MEAN ACC: ",ma->mu,10);
         }
         if (trace&T_VAC) {
            if (mpdf->ckind==DIAGC){
               printf("VAR OCC: %.2f\n",va->occ);
               ShowVector("VAR ACC: ",va->cov.var,10);
            } else {
               printf("INV OCC: %.2f\n",va->occ);
               ShowMatrix("INV ACC: ",va->cov.inv,10,10);
            }
         }
      }
   }
   if (trace&T_WAC){
      ShowSegNum(seg);
      printf("State %d, Stream %d\n",j,s);
      printf("WT OCC: %.2f\n",wa->occ);
      ShowVector("WT ACC: ",wa->c,10);
   }     
}
   
/* UpPDFCounts: update output PDF counts for each stream of each state */
void UpPDFCounts(LogDouble pr, int seg)
{
   int j,s;
   StateInfo *si;	/*hmodel*/
   StreamElem *se;	/*hmodel*/
   DVector alStarj, betj;

   for (j=2; j<nStates; j++) {
      si = hmm->svec[j].info;
      alStarj=alphaStar[j]; betj=beta[j];
	  if(si->dtype>0)
		  UpDurCounts(j, si, pr, alStarj,betj);	

	  /*alj = alphaStar[j]; betj = beta[j];*/
      for (s=1,se = si->pdf+1; s<=nStreams; s++,se++)
          UpStreamCounts(j,s,se,hset.swidth[s],pr,seg,alStarj,betj);		/*hrest*/
   }
}


/* UpdateCounters: update the various counters */
void UpdateCounters(LogDouble pr, int seg)
{
   SetOccr(pr,seg);		/*hrest*/
   if (uFlags&UPTRANS) 
      UpTranCounts(pr,seg);	/*hrest*/
	
   if (uFlags&(UPMEANS|UPVARS|UPMIXES|UPDURS))
      UpPDFCounts(pr,seg);	/*hrest*/

}

/* ------------------------- Model Update ----------------------- */

/* RestTransP: reestimate transition probs */
void RestTransP(void)
{
   int i,j;
   float occi,x,sum;
   TrAcc *ta;	/*htrain*/

   ta = (TrAcc *) GetHook(hmm->transP);	/*read transition accumulators into hmm->transP, hmem*/
   
   for (i=1;i<nStates;i++) {		/*update row by row, state by state */
      hmm->transP[i][1] = LZERO;
      occi = ta->occ[i];		/*total expected number of times in state i */
      if (occi == 0.0)
         HError(2222,"RestTransP: Zero state %d occupation count",i);
      sum = 0.0;
      for (j=2;j<=nStates;j++) {
         x = ta->tran[i][j]/occi;		/* x = Prob{ i->j | i } = E{i->j}/E{i} */
         hmm->transP[i][j] = x; sum += x;
      }
      for (j=2;j<=nStates;j++) {
         x = hmm->transP[i][j]/sum;	/* normalize Prob{i->j|i}, guarantee sum to 1 */
         hmm->transP[i][j] = (x<MINLARG) ? LZERO : log(x);
      }
   }
   if (trace & T_TRE)
      ShowMatrix("NEW TRANS: ",hmm->transP,10,10);	/*hmath*/
}

/* FloorMixes: apply floor to given mix set */
void FloorMixes(MixtureElem *mixes, int M, float floor)
{
   float sum,fsum,scale;
   MixtureElem *me;		/*hmodel*/
   int m;
   
   sum = fsum = 0.0;
   for (m=1,me=mixes; m<=M; m++,me++) {
      if (me->weight>floor)
         sum += me->weight;
      else {
         fsum += floor; me->weight = floor;
      }
   }
   if (fsum>1.0)
      HError(2223,"FloorMixes: Floor sum too large");
   scale = (1.0-fsum)/sum;
   if (trace&T_WRE) printf("MIXW: ");
   for (m=1,me=mixes; m<=M; m++,me++){
      if (me->weight>floor)
         me->weight *= scale;
      if (trace&T_WRE) printf(" %.2f",me->weight);
   }
   if (trace&T_WRE) printf("\n");
}  

/* FloorTMMixes: apply floor to given tied mix set */
void FloorTMMixes(Vector mixes, int M, float floor)
{
   float sum,fsum,scale,fltWt;
   int m;
   
   sum = fsum = 0.0;
   for (m=1; m<=M; m++) {
      fltWt = mixes[m];
      if (fltWt>floor)
         sum += fltWt;
      else {
         fsum += floor;
         mixes[m] = floor;
      }
   }
   if (fsum>1.0) HError(2223,"FloorTMMixes: Floor sum too large");
   scale = (1.0-fsum)/sum;
   if (trace&T_WRE) printf("MIXW: ");
   for (m=1; m<=M; m++){
      fltWt = mixes[m];
      if (fltWt>floor)
         mixes[m] = fltWt*scale;
      if (trace&T_WRE) printf(" %.2f",fltWt);
   }
}

/* FloorDProbs: apply floor to given discrete prob set */
void FloorDProbs(ShortVec mixes, int M, float floor)
{
   float sum,fsum,scale,fltWt;
   int m;
   
   sum = fsum = 0.0;
   for (m=1; m<=M; m++) {
      fltWt = Short2DProb(mixes[m]);	/*hmodel*/
      if (fltWt>floor)
         sum += fltWt;
      else {
         fsum += floor;
         mixes[m] = DProb2Short(floor);		/*hmodel*/
      }
   }
   if (fsum>1.0) HError(2327,"FloorDProbs: Floor sum too large");
   if (fsum == 0.0) return;
   if (sum == 0.0) HError(2328,"FloorDProbs: No probabilities above floor");
   scale = (1.0-fsum)/sum;
   for (m=1; m<=M; m++){
      fltWt = Short2DProb(mixes[m]);	/*hmodel*/
      if (fltWt>floor)
         mixes[m] = DProb2Short(fltWt*scale);	/*hmodel*/
   }
}

/* RestMixWeights: reestimate the mixture weights */
void RestMixWeights(int state, int s, StreamElem *se)
{
   WtAcc *wa;
   int m,M;
   float x;
   MixtureElem *me;
   
   wa = (WtAcc *)se->hook;	/*WtAcc, Mixture weight accumulator*/
   if (wa->occ == 0.0)
      HError(2222,"RestMixWeights: Zero weight occupation count");
   switch (hsKind){
   case TIEDHS:
      M=hset.tmRecs[s].nMix;
      break;
   case PLAINHS:
   case SHAREDHS:
   case DISCRETEHS:
      M=se->nMix;
      break;
   }
   for (m=1; m<=M; m++){
      x = wa->c[m]/wa->occ;	/* x = Prob {frequency in Component m} */
      if (x>1.0) 
         HError(2290,"RestMixWeights: Mix wt>1 in %d.%d.%d",state,s,m);
      switch (hsKind){
      case DISCRETEHS:
         se->spdf.dpdf[m] = (x>MINMIX) ? DProb2Short(x) : DLOGZERO;
         break;
      case TIEDHS:
         se->spdf.tpdf[m] = (x>MINMIX) ? x : 0.0;
         break;
      case PLAINHS:
      case SHAREDHS:
         me=se->spdf.cpdf+m;
         me->weight = (x>MINMIX) ? x : 0.0;		/*output*/
         break;
      }      
   }
}

/* RestMean: reestimate the given mean vector */
void RestMean(Vector mean, int vSize)
{
   int k;
   MuAcc *ma;
   float x;
   
   ma = (MuAcc *)GetHook(mean);		
   if (ma->occ == 0.0)
      HError(2222,"RestMean: Zero mean occupation count");
   for (k=1; k<=vSize; k++){
      x = mean[k] + ma->mu[k]/ma->occ;	/*compute new mean*/
      ma->mu[k] = mean[k];  /* remember old mean */
      mean[k] = x;		/*update to new mean*/
   }
   if (trace&T_MRE)
      ShowVector("MEAN: ",mean,10);
}

/* RestCoVar: reestimate the given covariance and return FALSE
              if any diagonal component == 0.0 */
Boolean RestCoVar(MixPDF *mp, int vSize, Vector minV,
                  Vector oldMean, Vector newMean, Boolean shared)
{
   int k,l;
   VaAcc *va;
   float x,z;
   float muDiffk,muDiffl;
   
   va = (VaAcc *)GetHook(mp->cov.var);
   if (va->occ == 0.0)
      HError(2222,"RestCoVar: Zero variance occupation count");
   if (mp->ckind==DIAGC){
      for (k=1; k<=vSize; k++){
         muDiffk = (shared)?0.0:newMean[k]-oldMean[k];
         x = va->cov.var[k] / va->occ - muDiffk*muDiffk;
         if (x<minV[k]) x = minV[k];
         if (x<vDefunct) return FALSE;
         mp->cov.var[k] = x;
      }
      FixDiagGConst(mp);	/*hmodel*/
      if (trace&T_VRE)
         ShowVector("VARS: ",mp->cov.var,10);
   } else {
      for (k=1; k<=vSize; k++){
         muDiffk = (shared)?0.0:newMean[k]-oldMean[k];
         for (l=1; l<k; l++) {
            muDiffl = (shared)?0.0:newMean[l]-oldMean[l];
            x = va->cov.inv[k][l] / va->occ - muDiffk*muDiffl;
            mp->cov.inv[k][l] = x;
         }
         z = va->cov.inv[k][k]/va->occ - muDiffk*muDiffk;
         mp->cov.inv[k][k] = (z<minV[k])?minV[k]:z;
      }
      if (trace&T_VRE)
         ShowTriMat("COVM: ",mp->cov.inv,10,10);

      FixFullGConst(mp,CovInvert(mp->cov.inv,mp->cov.inv));	/*hmodel*/
   }
   return TRUE;
}


/* RestDur: reestimate Poisson PDF parameters */

void RestDur(Vector durvec)
{
   DurAcc *da;
   float x;
   int i;
   
   da = (DurAcc *)GetHook(durvec);
   maxD=VectorSize(durvec);
   for(i=1;i<=maxD;i++) {
	  /*//x=da->dur[i]/da->occ;		//unsmoothed estimation*/
      x=(da->dur[i]+DURFLOOR)/(da->occ+maxD*DURFLOOR);   /*//smoothing, CS346 notes*/

	  durvec[i] = (x<MINLARG) ? LZERO : log(x);
   }
}

/* RestStream: reestimate stream parameters */
void RestStream(int state, int s, StreamElem *se, int vSize)
{
   int m,M;
   MixtureElem *me;	/*hmodel*/
   MixPDF *mp;	/*hmodel*/
   MuAcc *ma;
   Boolean shared;
   float wght;

   if (trace&(T_WRE|T_MRE|T_VRE))
      printf("State %d, Stream %d\n",state,s);
   if (uFlags&UPMIXES)
      RestMixWeights(state,s,se);
   if ((hsKind != DISCRETEHS)&&(hsKind != TIEDHS)){ /*wts only DI'ETE & TIED*/
      M=se->nMix;
      for (m=1; m<=M; m++){
         me = se->spdf.cpdf+m;
         wght=me->weight;
         mp=me->mpdf;
         if (wght > MINMIX) {
            if (trace&(T_MRE|T_VRE) && M>1)
               printf("Mixture %d\n",m);
            if (uFlags&UPMEANS)
               RestMean(mp->mean,vSize);	/*reestimate mean in current mixture */
            /* NB old mean left in ma->mu */
            if (uFlags&UPVARS){
               shared = GetUse(mp->cov.var) > 1;	/*hmem*/
               ma = (MuAcc *)GetHook(mp->mean);
               if ( !RestCoVar(mp,vSize,vFloor[s],ma->mu,mp->mean,shared)) {
                  if (M > 1) {
                     HError(-2225,"RestStream: Defunct Mix %d.%d.%d",state,s,m);
                     me->weight = 0.0;
                  } else
                     HError(2222,"RestStream: Zero Covariance in %d.%d",state,s);
               }
            }
         }
      }
   }
   if (hsKind == TIEDHS)
      M=hset.tmRecs[s].nMix;
   else
      M=se->nMix;
   if (M>1){
      switch (hsKind){
      case DISCRETEHS:
         FloorDProbs(se->spdf.dpdf,M,mixWeightFloor);	/*hrest*/
         break;
      case TIEDHS:
         FloorTMMixes(se->spdf.tpdf,M,mixWeightFloor);	/*hrest*/
         break;
      case PLAINHS:
      case SHAREDHS:
         FloorMixes(se->spdf.cpdf+1,M,mixWeightFloor);	/*hrest*/
         break;
      }     
   }
}

/* UpdateTheModel: use accumulated statistics to update model */
void UpdateTheModel(void)
{
   int j,s;
   StateInfo *si;	/*hmodel*/
   StreamElem *se;	/*hmodel*/

   if (uFlags&UPTRANS)
      RestTransP();
   if (uFlags&(UPMEANS|UPVARS|UPMIXES|UPDURS))
      for (j=2; j<nStates; j++) {
         si = hmm->svec[j].info;
		 /*//if (uFlags&UPDURS)*/
		 if(si->dtype>0)
			RestDur(si->dur);
         for (s=1,se = si->pdf+1; s<=nStreams; s++,se++)
            RestStream(j,s,se,hset.swidth[s]);
      }
   if (uFlags&UPVARS)
      FixAllGConsts(&hset);	/*hmodel*/
}

/* ------------------------- Top Level Control ----------------------- */


/* ReEstimateModel: top level of algorithm */
void ReEstimateModel(void)
{
   LogFloat segProb,oldP,newP,delta;
   LogDouble ap,bp;
   int converged,iteration,seg;

   iteration=0; 
   oldP=LZERO;
   do {        /*main re-est loop*/   
      ZeroAccs(&hset); newP = 0.0; ++iteration;
      nTokUsed = 0;
	  /*//SetDurP();*/	/*procompute the Duration PDF and store them in durprob*/
      for (seg=1;seg<=nSeg;seg++) {
         T=SegLength(segStore,seg);	/*htrain*/
         SetOutP(seg);		/*compute the observation prob and individual mixture observation probability*/
		
         /*//if ((ap=SetAlpha(seg)) > LSMALL){	//Compute alpha matrix, ap=alpha[nStates][T]*/
         /*   bp = SetBeta(seg);	//compute beta matrix, bp=beta[1][1]*/
		 segProb=SetAlphaBeta(seg);
		 if(segProb>LSMALL){/*
         //   if (trace & T_LGP)
         //      printf("%d.  Pa = %e, Pb = %e, Diff = %e\n",seg,ap,bp,ap-bp);
         //   segProb = (ap + bp) / 2.0;  *//* reduce numeric error */
            newP += segProb; ++nTokUsed;
            UpdateCounters(segProb,seg);	/*update transition counter and PDF counter*/
         } else
            if (trace&T_TOP)		/*likelihood is too small*/
               printf("Example %d skipped\n",seg);
      }
      if (nTokUsed==0)
         HError(2226,"ReEstimateModel: No Usable Training Examples");
      UpdateTheModel();		/*Reestimate all the model parameters*/
      newP /= nTokUsed;
      delta=newP-oldP; oldP=newP;
      converged=(fabs(delta)<epsilon); 
      if (trace&T_TOP) {
         printf("Ave LogProb at iter %d = %10.5f using %d examples",
                iteration,oldP,nTokUsed);
         if (iteration > 1)
            printf("  change = %10.5f",delta);
         printf("\n");
         fflush(stdout);
      }
   } while ((iteration < maxIter) && !converged);
   if (trace&T_TOP) {
      if (converged)
         printf("Estimation converged at iteration %d\n",iteration);
      else
         printf("Estimation aborted at iteration %d\n",iteration);
      fflush(stdout);
   }
}
/* ----------------------------------------------------------- */
/*                      END:  HRest.c                          */
/* ----------------------------------------------------------- */
.