/* UpWeights: update given mixture weights */
void UpWeights(int i, int s, int M, WtAcc *wa, StreamElem *ste)
{
   int m;
   float sum=0.0;
   
   if (wa->occ == 0.0)
      HError(2127,"UpWeights: zero occ i=%d/s=%d",i,s);
   for (m=1; m<=M; m++){
      sum += wa->c[m];
      switch(hset.hsKind){
      case PLAINHS:
      case SHAREDHS:
         ste->spdf.cpdf[m].weight = wa->c[m] / wa->occ;
         break;
      case TIEDHS:
         ste->spdf.tpdf[m] = wa->c[m] / wa->occ;
         break;
      }
   }
   if (fabs(sum-wa->occ)/sum > 0.001)
      HError(2190,"UpWeights: mix weight sum error");
}

/* UpMeans: update mean, leave old mean in acc */
void UpMeans(int i, int s, int m, int size, MuAcc *ma, Vector mean)
{
   int k;
   float x;
   
   if (ma->occ == 0.0)
      HError(2127,"UpMeans: zero occ i=%d/s=%d/m=%d",i,s,m);
   for (k=1; k<=size; k++){
      x = mean[k] + ma->mu[k]/ma->occ;
      ma->mu[k] = mean[k];  /* remember old mean */
      if (uFlags&UPMEANS) mean[k] = x;
   }
}

/* FloorMixes: apply floor to given mix set */
void FloorMixes(HMMSet *hset, MixtureElem *mixes, int M, float floor)
{
   float sum,fsum,scale;
   MixtureElem *me;
   int m;
   
   if (hset->logWt == TRUE) HError(999,"FloorMixes requires linear weights");
   sum = fsum = 0.0;
   for (m=1,me=mixes; m<=M; m++,me++) {
      if (MixWeight(hset,me->weight)>floor)
         sum += me->weight;
      else {
         fsum += floor; me->weight = floor;
      }
   }
   if (fsum>1.0) HError(2327,"FloorMixes: Floor sum too large");
   if (fsum == 0.0) return;
   if (sum == 0.0) HError(2328,"FloorMixes: No mixture weights above floor");
   scale = (1.0-fsum)/sum;
   for (m=1,me=mixes; m<=M; m++,me++)
      if (me->weight>floor) me->weight *= scale;
}

void CheckUpdateSetUp()
{
  AdaptXForm *xf;

  xf = xfInfo.paXForm;
  if ((xfInfo.paXForm != NULL) && !(uFlags&UPXFORM)) {
    while (xf != NULL) {
       if ((xf->xformSet->xkind != CMLLR) && (xf->xformSet->xkind != SEMIT))
	HError(999,"SAT only supported with SEMIT/CMLLR transforms");
      xf = xf->parentXForm;
    }
  }
}

/* GenSentences: top level control of the sentence generator */
void  GenSentences(char * latfn, char * dicfn)
{
   int i,min,max,len;
   double e,p;
   MemHeap lheap;
   FILE *f;
   Boolean isPipe;

   InitVocab(&voc);
   if(ReadDict(dicfn,&voc)<SUCCESS)
      HError(3413,"GenSententces:ReadDict failed" );
   CreateHeap(&lheap,"Lattice Heap",MSTAK,1,0.4,1000,5000);
   if ((f=FOpen(latfn,NetFilter,&isPipe)) == NULL)
      HError(3410,"GenSentences: Can't open lattice file %s",latfn);
   if((lat = ReadLattice(f, &lheap, &voc, TRUE, FALSE))==NULL)
      HError(3410,"GenSentences: ReadLattice failed");
   FClose(f,isPipe);

   if (trace&T_TOP)
      printf("HSGen %d sents from lattice %s/dictionary %s\n",
             ngen,latfn,dicfn);
   psSum = 0.0; lenSum = 0; min = 100000; max = 0;
   if (trace&T_DET) quiet = TRUE;  /* kill output if detailed trace */
   for (i=1; i<=ngen; i++){
      len = GenSent(i);
      lenSum += len;
      if (len>max) max = len;
      if (len<min) min = len;
   }
   if (stats)  {
      ComputeVSize();
      e = psSum / lenSum;
      p = exp(e);
      e = e / log(2.0);
      printf("Entropy = %f,  Perplexity = %f\n",e,p);
      printf("%d Sentences: average len = %.1f, min=%d, max=%d\n",
             ngen,(float)lenSum/ngen,min,max);
   }
}

/* GetNode: read a node definition and create node */
static VQNode GetNode(Source *src, CovKind ck, short width)
{
   char buf[MAXSTRLEN];
   VQNode n;
   short vqidx,nid,lid,rid;
   Vector mean;
   Covariance cov;

   vqidx = GetVal(src,0,0,"VQ Index");
   nid = GetVal(src,0,0,"Node Id");
   lid = GetVal(src,0,0,"Left Id");
   rid = GetVal(src,0,0,"Right Id");
   mean = CreateVector(&vqHeap,width);
   if (!ReadVector(src, mean, FALSE))
      HError(6150,"GetNode: cannot read mean vector at %s",
             SrcPosition(*src, buf));
   switch(ck){
   case NULLC:
      cov.var = NULL;
      n = CreateVQNode(vqidx,nid,lid,rid,mean,ck,cov);
      break;
   case INVDIAGC:
      cov.var = CreateVector(&vqHeap,width);
      if (!ReadVector(src, cov.var, FALSE))
         HError(6150,"GetNode: cannot read variance vector at %s",
                SrcPosition(*src, buf));
      n = CreateVQNode(vqidx,nid,lid,rid,mean,ck,cov);
      break;
   case FULLC:
      cov.inv = CreateTriMat(&vqHeap,width);
      if (!ReadTriMat(src, cov.inv, FALSE))
         HError(6150,"GetNode: cannot read covariance matrix at %s",
                SrcPosition(*src, buf));
      n = CreateVQNode(vqidx,nid,lid,rid,mean,ck,cov);
      break;
   }
   return n;
}

/* EXPORT->LoadMasterFile: Load the Master Label File stored in fname 
                           and append the entries to the MLF table */
void LoadMasterFile(char *fname)
{
   char buf[MAXFNAMELEN];
   char *men;        /* end of mode indicator */
   char *pst,*pen;   /* start/end of pattern (inc quotes) */
   char *dst=NULL,*den=NULL;   /* start/end of subdirectory (inc quotes) */
   Boolean inEntry = FALSE;   /* ignore ".." within an entry */
   MLFEntry *e;
   FILE *f;
   
   if (numMLFs == MAXMLFS)
      HError(6520,"LoadMasterFile: MLF file limit reached [%d]",MAXMLFS);
   if ((f = fopen(fname,"rb")) == NULL)
      HError(6510,"LoadMasterFile: cannot open MLF %s",fname);
   if (fgets(buf,MAXFNAMELEN,f) == NULL)
      HError(6513,"LoadMasterFile: MLF file is empty");
   if (NoMLFHeader(buf))
      HError(6551,"LoadMasterFile: MLF file header is missing"); 
   incSpaces=FALSE;
   while (fgets(buf,MAXFNAMELEN,f) != NULL){
      if (!inEntry && FindMLFStr(buf,&pst,&pen)) {
         e = (MLFEntry *)New(&mlfHeap,sizeof(MLFEntry));
         e->type = FindMLFType(pen+1,&men);
         if (e->type == MLF_IMMEDIATE) {
            e->def.immed.fidx = numMLFs;
            e->def.immed.offset = ftell(f);
            if (e->def.immed.offset < 0)
               HError(6521,"LoadMasterFile: cant ftell on MLF file");
            inEntry = TRUE;
         } else {
            if (!FindMLFStr(men+1,&dst,&den))
               HError(6551,"LoadMasterFile: Missing subdir in MLF\n(%s)",buf);
            *den = '\0';
            e->def.subdir = NewString(&mlfHeap,den-dst-1);
            strcpy(e->def.subdir,dst+1);
         }
         *pen = '\0';         /* overwrite trailing pattern quote */
         ++pst;               /* skipover leading pattern quote */
         e->patType = ClassifyMLFPattern(pst);
         if (e->patType == PAT_ANYPATH) 
            pst += 2;         /* skipover leading "* /" */
         e->pattern = NewString(&mlfHeap,pen-pst);
         strcpy(e->pattern,pst);
         e->patHash = (e->patType==PAT_GENERAL)?0:MLFHash(e->pattern);
         StoreMLFEntry(e);
      } else
         if (inEntry && IsDotLine(buf)) inEntry = FALSE;
   }
   if (compatMode && incSpaces)
      HError(-6551,"LoadMasterFile: . or %s on line with spaces in %s",
             LEVELSEP,fname);
   mlfile[numMLFs++] = f;
}

/* EXPORT->ASpec2LPCep: Perform IDFT to get autocorrelation values then 
           produce autoregressive coeffs. and cepstral transform them */
void ASpec2LPCep (Vector as, Vector ac, Vector lp, Vector c, DMatrix cm)
{
   float lpcGain, E;

   /* Do IDFT to get autocorrelation values */
   E = MatrixIDFT(as, ac, cm);
   lp[VectorSize(lp)] = 0.0;    /* init to make Purify et al. happy */
   /* do Durbin recursion to get predictor coefficients */
   lpcGain = Durbin(NULL,lp,ac,E,VectorSize(ac)-1);
   if (lpcGain<=0) 
      HError(-5323,"ASpec2LPCep: Negative lpcgain");
   LPC2Cepstrum(lp,c);
   c[VectorSize(c)] = (float) -log((double) 1.0/lpcGain); /* value forms C0 */
}

/* StoreOOV: store OOV wdid/count into ps */
static void StoreOOV(PStats *ps, LabId wdid, int count)
{
   int n;

   ps->oov[ps->uniqOOV].wdid  = wdid;
   ps->oov[ps->uniqOOV].count = count;
   ps->uniqOOV++; ps->nOOV++;
   if (ps->uniqOOV==MAX_OOV) {
     n = SortOOV(ps);
     printf("StoreOOV: sorting OOVs, compacting %d -> %d\n",MAX_OOV,n);
     if (n==MAX_OOV)
       HError(16630,"Maximum number of unique OOV's [%d] reached\n",MAX_OOV);
  }
}

/* UpTrans: update transition parameters */
void UpTrans(TrAcc *ta, Matrix tr)
{
   int i,j;
   float occi,x,sum;

   for (i=1; i<nStates; i++){
      occi = ta->occ[i];
      if (occi == 0.0)
         HError(2127,"UpTrans: zero occ in state %d",i);
      sum = 0.0;
      tr[i][1] = LZERO;
      for (j=2;j<=nStates;j++) {
         x = ta->tran[i][j]/occi;
         tr[i][j] = x; sum += x;
      }
      if (fabs(sum-1.0) > 0.001)
         HError(2190,"UpTrans: row %d, sum=%f",i,sum,occi);
      for (j=2;j<=nStates;j++) {
         x = tr[i][j]/sum;
         tr[i][j] = (x<MINLARG) ? LZERO : log(x);
      }
   }
}

/* FindVQTable: find VQ table with given name and unless magic is 0, check
                that it is same, return NULL if not found */
static VQTable FindVQTable(char * tabFN, short magic)
{
   Boolean found=FALSE;
   VQTable p;

   for (p=vqList; p!=NULL && !found; p=p->next)
      if (strcmp(tabFN,p->tabFN)==0){
         if (magic != 0 && magic != p->magic)
            HError(6170,"FindVQTable: %s has magic=%d but new magic=%d",
                   tabFN,p->magic,magic);
         return p;
      }
   return NULL;
}

/* FindBestMixes: for each state/obs pair find most likely mix component */
void FindBestMixes(int segNum, int segLen, IntVec states, IntVec *mixes)
{
   int i,s,m,bestm,M=0;
   StreamElem *ste;
   IntVec smix;
   Observation obs;
   Vector v;
   LogFloat bestP,p;
   MixtureElem *me;
   MixPDF *mp;

   if (trace&T_MIX)
      printf(" Mixture component alignment\n");
   for (i=1; i<=segLen; i++){
      ste = hmmLink->svec[states[i]].info->pdf+1;
      obs = GetSegObs(segStore, segNum, i);
      if (hset.hsKind == TIEDHS)
         PrecomputeTMix(&hset, &obs, 0.0, 1);
      for (s=1; s<=nStreams; s++,ste++){
         if (hset.hsKind != TIEDHS)
            M = ste->nMix;
         smix = mixes[s];
         if (hset.hsKind==TIEDHS) /* PrecomputeTMix has already sorted probs */
            bestm = hset.tmRecs[s].probs[1].index;
         else if (M==1)
            bestm = 1;   
         else{
            v = obs.fv[s];
            bestP = LZERO; bestm=0;
            if (trace&T_MIX)
               printf("  seg %d, stream %d: ",i,s);
            for (m=1; m<=M; m++){
               me =  ste->spdf.cpdf+m;
               mp = me->mpdf;
               p = MOutP(v,mp);
               if (p>bestP){
                  bestP=p; bestm=m;
               }
               if (trace&T_MIX)
                  printf(" P(mix[%d])=%.1f",m,p);
            }
            if (bestm==0)
               HError(2125,"FindBestMixes: no best mix");
            if (trace&T_MIX)
               printf(" [best=%d]\n",bestm);
         }
         smix[i] = bestm;
      }
   }
}

void LoadNetwork()
{
   FILE *nf;
   Boolean isPipe;
   int n=0;

   CreateHeap(&wdNetHeap,"Lattice heap",MSTAK,1,0.0,4000,4000);
   if ( (nf = FOpen(wdNetFn,NetFilter,&isPipe)) == NULL)
      HError(3210,"LoadNetwork: Cannot open Word Net file %s",wdNetFn);
   if((wdNet = ReadLattice(nf,&wdNetHeap,&vocab,TRUE,FALSE))==NULL)
      HError(3210,"LoadNetwork: ReadLattice failed");
   FClose(nf,isPipe);

   printf("Read Word Network with %d nodes / %d arcs\n",wdNet->nn,wdNet->na);

   CreateHeap(&netHeap,"Net heap",MSTAK,1,0,
      wdNet->na*sizeof(NetLink),wdNet->na*sizeof(NetLink));

   net = ExpandWordNet(&netHeap,wdNet,&vocab,&hset);
   printf("Created network with %d nodes / %d links\n",
      net->numNode,net->numLink);

}

/* EXPORT->GetVQ: get vq indices for vectors in fv */
void GetVQ(VQTable vqTab, int numS, Vector *fv, short *vq)
{
   short s,idx,size;
   float bestx, x,xl,xr;
   VQNode n,bestn;
   Vector v;

   for (s=1; s<=numS; s++) {
      n = vqTab->tree[s]; v = fv[s];
      size = VectorSize(v);
      if (n==NULL)
         HError(6174,"GetVQ: null tree in stream %d",s);
      if (size != vqTab->swidth[s])
         HError(6174,"GetVQ: stream %d width incompatible",s);
      switch(vqTab->type){
      case linTree:
         bestn = n; bestx = VQNodeScore(n,v,size,vqTab->ckind);
         for(n=bestn->right; n != NULL; n=n->right){
            x = VQNodeScore(n,v,size,vqTab->ckind);
            if (x<bestx) {
               bestx = x; bestn = n;
            }
         }
         idx = bestn->vqidx;
         break;
      case binTree:
         while (n->right != NULL){
            xr = VQNodeScore(n->right,v,size,vqTab->ckind);
            xl = VQNodeScore(n->left,v,size,vqTab->ckind);
            n =  (xr<xl)?n->right:n->left;
         }
         idx = n->vqidx;
         break;
      }
      vq[s] = idx;
   }
}

/* LoadHTKLabels: load a HTK transcription */
static void LoadHTKLabels(MemHeap *x, Transcription *t, Source *src)
{
   LabList *ll;
   int alt = 0;
   
   InitTrScan();
   GetTrSym(src,TRUE);
   if (trSym==TRNUM || trSym==TRSTR){
      ll = LoadHTKList(x,src,++alt);
      AddLabelList(ll,t);
      while (trSym == TRLEV){
         GetTrSym(src,TRUE);
         if (trSym != TREOL)
            HError(6550,"LoadHTKList: End of Line after /// Expected");       
         GetTrSym(src,TRUE);
         ll = LoadHTKList(x,src,++alt);
         AddLabelList(ll,t);
      }
   }
   while (trSym==TREOL) 
      GetTrSym(src,TRUE);
   if (trSym != TREOF)
      HError(6550,"LoadHTKLabels: Junk at end of HTK transcription");
}

/* CreateDataFileList: */
static void CreateDataFileList()
{
   char *datafn;
   g_nDataFileNum = 0;
   do {
      if (NextArg() != STRINGARG)
         HError(2319, "HMGenS: data file name expected");
      datafn = GetStrArg();

      g_pDataFileList[g_nDataFileNum] = (TDataFile *) New(&gstack, sizeof(TDataFile));
      strcpy(g_pDataFileList[g_nDataFileNum]->datafn, datafn);
      g_pDataFileList[g_nDataFileNum]->bValid = TRUE;
      g_nDataFileNum++;
   } while (NumArgs() > 0);
}

/* CheckGenSetUp: Check & setup GenInfo structure */
static void CheckGenSetUp(void)
{
   int d, p, r, s, stream;
   PdfStream *pst = NULL;

   /* # of PdfStreams */
   genInfo->nPdfStream[0] = (nPdfStr == NULL) ? hset.swidth[0] : IntVecSize(nPdfStr);
   if (genInfo->nPdfStream[0] > hset.swidth[0])
      HError(6604, "CheckGenSetUp: # of PdfStreams (%d) is too large (should be less than %d)", genInfo->nPdfStream[0], hset.swidth[0]);

   /* size of each PdfStreams */
   r = hset.swidth[0];
   for (p = stream = 1; p <= genInfo->nPdfStream[0]; stream += genInfo->nPdfStream[p++]) {
      pst = &(genInfo->pst[p]);

      /* # of streams in this PdfStream */
      genInfo->nPdfStream[p] = (nPdfStr == NULL) ? 1 : nPdfStr[p];
      r -= genInfo->nPdfStream[p];

      /* calculate vector size for this PdfStream */
      for (s = stream, pst->vSize = 0; s < stream + genInfo->nPdfStream[p]; s++)
         pst->vSize += hset.swidth[s];

      /* order (vecSize of static feature vector) of this PdfStream */
      pst->order = (pdfStrOrder == NULL) ? 1 : pdfStrOrder[p];
      if (pst->order < 1 || pst->order > pst->vSize)
         HError(6604, "CheckGenSetUp: Order of PdfStream %d should be within 1--%d", p, pst->vSize);

      /* window coefficients */
      if (winFn[p] == NULL)
         HError(6604, "CheckGenSetUp: window file names are not specified");
      pst->win.num = (int) winFn[p][0][0];
      if (pst->win.num > MAXWINNUM)
         HError(6604, "CheckGenSetUp: # of window out of range");
      if (pst->win.num * pst->order != pst->vSize)
         HError(6604, "CheckGenSetUp: # of window (%d) times order (%d) should be equal to vSize (%d)", pst->win.num, pst->order, pst->vSize);
      for (d = 0; d < pst->win.num; d++)
         MakeFN(winFn[p][d + 1], winDir, winExt, pst->win.fn[d]);
   }
   if (r != 0)
      HError(6604, "CheckGenSetUp: # of streams in HMMSet (%d) and PdfStreams (%d) are inconsistent", hset.swidth[0], genInfo->nPdfStream[0]);

   /* output trace information */
   if (trace & T_TOP) {
      for (p = 1; p <= genInfo->nPdfStream[0]; p++) {
         printf("PdfStream [%d]:\n", p);
         printf("  #streams: %d (vSize=%d)\n", genInfo->nPdfStream[p], genInfo->pst[p].vSize);
         printf("  #order:   %d\n", genInfo->pst[p].order);
         printf("  file ext: %s\n", genInfo->pst[p].ext);
         for (d = 0; d < genInfo->pst[p].win.num; d++)
            printf("  %d-th window: %s\n", d, genInfo->pst[p].win.fn[d]);
      }
      printf("\n");
      fflush(stdout);
   }
}

/* CalcCovs: calculate covariance of speech data */
void CalcCovs(void)
{
   int x,y,s,V;
   float meanx,meany,varxy,n;
   Matrix fullMat;
   
   if (totalCount<2)
      HError(2021,"CalcCovs: Only %d speech frames accumulated",totalCount);
   if (trace&T_TOP)
      printf("%ld speech frames accumulated\n", totalCount);
   n = (float)totalCount;     /* to prevent rounding to integer below */
   for (s=1; s<=hset.swidth[0]; s++){  /* For each stream   */
      V = hset.swidth[s];
      for (x=1; x<=V; x++)            /* For each coefficient ... */
         accs[s].meanSum[x] /= n;         /* ... calculate mean */
      for (x=1;x<=V;x++) {
         meanx = accs[s].meanSum[x];      /* ... and [co]variance */
         if (fullcNeeded[s]) {
            for (y=1; y<=x; y++) {
               meany = accs[s].meanSum[y];
               varxy = accs[s].squareSum.inv[x][y]/n - meanx*meany;
               accs[s].squareSum.inv[x][y] =
                  (x != y || varxy > minVar) ? varxy : minVar;    
            }
         }
         else {
            varxy = accs[s].squareSum.var[x]/n - meanx*meanx;
            accs[s].fixed.var[x] = (varxy > minVar) ? varxy :minVar;
         }
      }
      if (fullcNeeded[s]) { /* invert covariance matrix */
         fullMat=CreateMatrix(&gstack,V,V);
         ZeroMatrix(fullMat); 
         CovInvert(accs[s].squareSum.inv,fullMat);
         Mat2Tri(fullMat,accs[s].fixed.inv);
         FreeMatrix(&gstack,fullMat);
      }
      if (trace&T_COVS) {
         printf("Stream %d\n",s);
         if (meanUpdate)
            ShowVector(" Mean Vector ", accs[s].meanSum,12);
         if (fullcNeeded[s]) {
            ShowTriMat(" Covariance Matrix ",accs[s].squareSum.inv,12,12);
         } else
            ShowVector(" Variance Vector ", accs[s].fixed.var,12);
      }
   }
}

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

/* ClusterVecs: apply required clustering to vectors in stream s */
void ClusterVecs(Sequence *seq, int s)
{
   switch(tType) {
   case binTree:
      cs[s] = TreeCluster(&cStack,seq[s],cbSizes[s],ck,ck,cov[s]);
      break;
   
   case linTree:
      cs[s] = FlatCluster(&cStack,seq[s],cbSizes[s],ck,ck,cov[s]);
      break;
   
   default:
      HError(2590,"ClusterVecs: Unsupported tree type %d",tType);
   }
   if(trace&T_CLUST) ShowClusterSet(cs[s]);
}