/* ----------------------------------------------------------- *//*                                                             *//*                          ___                                *//*                       |_| | |_/   SPEECH                    *//*                       | | | | \   RECOGNITION               *//*                       =========   SOFTWARE                  */ /*                                                             *//*                                                             *//* ----------------------------------------------------------- *//* developed at:                                               *//*                                                             *//*      Speech Vision and Robotics group                       *//*      Cambridge University Engineering Department            *//*      http://svr-www.eng.cam.ac.uk/                          *//*                                                             *//*      Entropic Cambridge Research Laboratory                 *//*      (now part of Microsoft)                                *//*                                                             *//* ----------------------------------------------------------- *//*         Copyright: Microsoft Corporation                    *//*          1995-2000 Redmond, Washington USA                  *//*                    http://www.microsoft.com                 *//*                                                             *//*              2002  Cambridge University                     *//*                    Engineering Department                   *//*                                                             *//*   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: HTrain.c   HMM Training Support Routines      *//* ----------------------------------------------------------- */char *htrain_version = "!HVER!HTrain:   3.2 [CUED 09/12/02]";char *htrain_vc_id = "$Id: HTrain.c,v 1.8 2002/12/19 16:37:11 ge204 Exp $";#include "HShell.h"#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 "HUtil.h"#include "HTrain.h"/* --------------------------- Trace Flags ------------------------- */#define T_TOP  00001       /* Top Level tracing */#define T_SEQ  00002       /* Trace sequence operations */#define T_CGE  00004       /* General cluster tracing */#define T_CLC  00010       /* Trace cluster convergence */#define T_DCM  00020       /* Dump cluster map after each cycle */#define T_CAL  00040       /* Trace item to cluster allocation */#define T_CDI  00100       /* Trace cluster distance calc */#define T_NAC  00200       /* Count accs created */#define T_ALD  00400       /* Trace acc load/dump */static int trace = 0;/* --------------------------- Initialisation ---------------------- */static ConfParam *cParm[MAXGLOBS];      /* config parameters */static int nParm = 0;static int maxIter = 10;               /* max num cluster iterations */static int minClustSize = 3;           /* min num vectors in cluster */static Boolean ldBinary = TRUE;        /* load/dump in binary */#define DoPreComps(hsKind) (hsKind==SHAREDHS||hsKind==PLAINHS)/* EXPORT->InitTrain: initialise configuration parameters */void InitTrain(void){   int i;   Boolean b;   Register(htrain_version,htrain_vc_id);   nParm = GetConfig("HTRAIN", TRUE, cParm, MAXGLOBS);   if (nParm>0){      if (GetConfInt(cParm,nParm,"TRACE",&i)) trace = i;      if (GetConfInt(cParm,nParm,"MAXCLUSTITER",&i)) maxIter = i;      if (GetConfInt(cParm,nParm,"MINCLUSTSIZE",&i)) minClustSize = i;      if (GetConfBool(cParm,nParm,"BINARYACCFORMAT",&b)) ldBinary = b;   }}/* -------------------- Generic Sequence Type ------------------- *//*  Items are stored within chained lists of blocks.  Items are   normally added to the end of the last block.  If the last  block is full, then the item is inserted into the first block  with a free item slot.  If all blocks are full, then a new  block is created and chained to the end of the list.*/   /* NewItemBlock: create and return an emptry ItemBlock */static IBLink NewItemBlock(MemHeap *x, int blkSize){   IBLink b;      if (trace&T_SEQ)      printf("HTrain: creating item block, %d items on stack %s\n",             blkSize,x->name);   b = (IBLink)New(x,sizeof(ItemBlock));   b->used = 0; b->next = NULL;   b->items = (Ptr *)New(x,sizeof(Ptr)*blkSize);   return b;}/* EXPORT->CreateSequence: create and return an empty sequence */Sequence CreateSequence(MemHeap *x, int blkSize){   Sequence s;      s = (Sequence)New(x,sizeof(SequenceInfo));   s->mem = x;    s->nItems = 0;   s->nFree = s->blkSize = blkSize;   s->hd = s->tl = NewItemBlock(x,blkSize);   return s;}/* EXPORT->StoreItem: store item in seq, extend if needed */void StoreItem(Sequence seq, Ptr item){   IBLink b;   if (seq->nFree==0){     /* add a block */      b = NewItemBlock(seq->mem,seq->blkSize);      seq->tl->next = b; seq->tl = b;      seq->nFree = seq->blkSize;   }   if (seq->tl->used < seq->blkSize) { /* try last block first */      b = seq->tl;      b->items[b->used++] = item;   }else{   /* otherwise search blocks for empty slot */      for (b = seq->hd; b->used == seq->blkSize; b = b->next)         if (b==NULL)            HError(7190,"StoreItem: free slot missing");      b->items[b->used++] = item;   }   --seq->nFree; ++seq->nItems;}/* FindItemBlock: find block containing i'th item, j is index in block */static IBLink FindItemBlock(Sequence seq, int i, int *j){   IBLink b;      if (i > seq->nItems)      HError(7171,"FindItemBlock: %d'th item from %d requested",             i,seq->nItems);   b = seq->hd;   while (i>b->used){      i -= b->used; b = b->next;   }   *j = i-1;   return b;}/* EXPORT->DeleteItem: delete i'th item */void DeleteItem(Sequence seq, int i){   IBLink b;   int j,k;      b = FindItemBlock(seq,i,&j);   --b->used;   for (k=j; k<b->used; k++)      b->items[k] = b->items[k+1];   ++seq->nFree; --seq->nItems;}/* EXPORT->GetItem: return ptr to the i'th item in sequence */Ptr GetItem(Sequence seq, int i){   IBLink b;   int j;      b = FindItemBlock(seq,i,&j);   return b->items[j];  }/* ------------ Segment (Observation Sequence) Storage ------------ *//* EXPORT->CreateSegStore: Create and return an empty segment store   suitable for observations of form obs. */SegStore CreateSegStore(MemHeap *x, Observation obs, int segLen){   SegStore ss;      ss = (SegStore)New(x,sizeof(SegStoreRec));   ss->mem = x;   ss->o = obs; ss->segLen = segLen;   ss->hasfv = (obs.pk&BASEMASK) != DISCRETE;   ss->hasvq = (obs.pk&HASVQ)  || (obs.pk&BASEMASK) == DISCRETE;   if (ss->hasfv) ss->fvSegs = CreateSequence(x,100);   if (ss->hasvq) ss->vqSegs = CreateSequence(x,100);   return ss;}/* EXPORT->LoadSegment: of obs in pbuf from start to end */void LoadSegment(SegStore ss, HTime start, HTime end, ParmBuf pbuf){   Sequence vq;   Sequence fv;   BufferInfo info;   long i,st,en,len;   int s,S = ss->o.swidth[0];   short *vqItem;   Vector *fvItem;      GetBufferInfo(pbuf,&info);   st = (long)(start/info.tgtSampRate);   en = (long)(end/info.tgtSampRate);   len = ObsInBuffer(pbuf);   if (en >= len) en = len-1;   if (st > en) {      HError(-7173,"LoadSegment: empty segment");      return;   }   if (ss->hasvq) vq = CreateSequence(ss->mem,ss->segLen);   if (ss->hasfv) fv = CreateSequence(ss->mem,ss->segLen);   for (i=st; i<=en; i++) {      if (ss->hasfv) { /* put new vectors in ss->o */         for (s=1; s<=S; s++)             ss->o.fv[s] = CreateVector(ss->mem,ss->o.swidth[s]);              }      ReadAsTable(pbuf,i,&(ss->o));      if (ss->hasvq) {         vqItem = (short *)New(ss->mem,sizeof(short)*(S+1));         for (s=1; s<=S; s++)             vqItem[s] = ss->o.vq[s];         StoreItem(vq, (Ptr) vqItem);      }      if (ss->hasfv) {         fvItem = (Vector *)New(ss->mem,sizeof(Vector)*(S+1));         for (s=1; s<=S; s++)             fvItem[s] = ss->o.fv[s];         StoreItem(fv, (Ptr) fvItem);      }   }   if (ss->hasvq) StoreItem(ss->vqSegs, (Ptr)vq);   if (ss->hasfv) StoreItem(ss->fvSegs, (Ptr)fv);}/* EXPORT->SegLength: Return num obs in i'th segment */int SegLength(SegStore ss, int i){   Sequence s;      if (ss->hasfv)       s = (Sequence) GetItem(ss->fvSegs,i);   else if (ss->hasvq)      s = (Sequence) GetItem(ss->vqSegs,i);   else      HError(7191,"SegLength: no fv or vq seg stored");   return s->nItems;}/* EXPORT->NumSegs: Return num segments in ss */int NumSegs(SegStore ss){   if (ss->hasfv)       return ss->fvSegs->nItems;   else if (ss->hasvq)      return ss->vqSegs->nItems;   else      HError(7191,"NumSegs: no fv or vq segs stored");   return 0;  /* to keep compiler happy */}/* EXPORT->GetSegObs: Return j'th observation from i'th segment */Observation GetSegObs(SegStore ss, int i, int j){   Sequence vq;   Sequence fv;   int s,S = ss->o.swidth[0];   short *vqItem;   Vector *fvItem;      if (ss->hasvq) {      vq = (Sequence) GetItem(ss->vqSegs,i);      vqItem = (short *) GetItem(vq,j);      for (s=1; s<=S; s++)          ss->o.vq[s] = vqItem[s];   }        if (ss->hasfv) {      fv = (Sequence) GetItem(ss->fvSegs,i);      fvItem = (Vector *) GetItem(fv,j);      for (s=1; s<=S; s++)          ss->o.fv[s] = fvItem[s];   }        return ss->o;}/* EXPORT->SequenceMean: compute mean of vectors in ss */void SequenceMean(Sequence ss, Vector mean){   Vector v;   int i,j,size;   double n;   DVector sum;   size = VectorSize(mean);   sum = CreateDVector(&gstack,size);   ZeroDVector(sum);   n = ss->nItems;   for (i=1; i<=ss->nItems; i++) {      v = (Vector) GetItem(ss,i);      for (j=1; j<=size; j++)         sum[j] += v[j];   }   for (j=1; j<=size; j++)      mean[j] = sum[j]/n;   FreeDVector(&gstack,sum);}/* EXPORT->SequenceCov: compute covariance of vectors in ss */void SequenceCov(Sequence ss, CovKind ck, Covariance cov, Vector mean){   Vector v;   int i,j,k,size;   double n,x,y;   DVector sqsum;   DMatrix xsum;      switch(ck){   case DIAGC:       /* diagonal covariance matrix */   case INVDIAGC:    /* inverse diag covariance matrix */      size = VectorSize(cov.var);      sqsum = CreateDVector(&gstack,size);      ZeroDVector(sqsum);      n = ss->nItems;      for (i=1; i<=ss->nItems; i++) {         v = (Vector) GetItem(ss,i);         for (j=1; j<=size; j++){            x = v[j]-mean[j];            sqsum[j] += x*x;         }      }      for (j=1; j<=size; j++)         cov.var[j] = (ck==DIAGC)?sqsum[j]/n:n/sqsum[j];      FreeDVector(&gstack,sqsum);      break;   case FULLC:    /* inverse full covariance matrix */      size = TriMatSize(cov.inv);      xsum = CreateDMatrix(&gstack,size,size);      ZeroDMatrix(xsum);      n = ss->nItems;      for (i=1; i<=ss->nItems; i++) {         v = (Vector) GetItem(ss,i);         for (j=1; j<=size; j++)            for (k=1; k<=j; k++){               x = v[j]-mean[j];               y = v[k]-mean[k];               xsum[j][k] += x*y;            }      }      for (j=1; j<=size; j++)         for (k=1; k<=j; k++)            cov.inv[j][k] = xsum[j][k]/n;      CovInvert(cov.inv,cov.inv);      FreeDMatrix(&gstack,xsum);      break;   default:      HError(7170,"SequenceCov: unsupported cov kind [%d]",ck);   }}/* --------------------- Vector Clustering -------------------- *//*    The routines in this section implement a top-down clustering   algorithm.  There are two clustering modes: linear and tree   based.  Both are top-down and start from a single cluster.      In the default linear mode, the largest cluster is   successively split until the required number of clusters are   formed.  Each cluster centre is the average of all vectors in the   cluster.  Splitting is performed by perturbing the cluster centre to   form two cluster centres.  The largest cluster is defined as that   which has the largest average construction cost.      In tree clustering mode, the current leaves are split regardless   of their size.      Linear mode is best for clustering small data sets (eg as in HInit).   Either mode can be used for large data sets (eg as in HQuant).*/static int curNumCl;    /* num clust currently created */static int nItems;      /* num items in pool */static ClusterSet *ccs; /* current cluster set */static Sequence cvp;    /* current vector pool being clustered */static ShortVec cmap;   /* array[1..cvp->nItems]of cluster index */static CovKind dck;     /* defines type of distance calc */static Covariance dcov; /* covariance to use in distance calc */static int vSize;       /* size of vectors */static Vector vTmp;     /* temp vector *//* DumpClusterMap: dump the vector pool -> cluster map */static void DumpClusterMap(void){   int i,lc=0;      for (i=1; i<=nItems; i++) {      printf("%3d",cmap[i]);      if (++lc == 20) { printf("\n"); lc = 0; }   }   if (lc != 0) printf("\n");   fflush(stdout);}/* Distance: compute distance between v1 and v2 */static float Distance(Vector v1, Vector v2){   Vector iv,crow;   TriMat ic;   double sum,x;   int i,j;   switch(dck){   case NULLC:      sum = 0.0;      for (i=1; i<=vSize; i++){         x = v1[i]-v2[i]; sum += x*x;      }      break;   case DIAGC:      iv = dcov.var;  /* covkind == DIAGC */      sum = 0.0;      for (i=1; i<=vSize; i++){         x = v1[i]-v2[i]; sum += x*x/iv[i];      }      break;   case INVDIAGC:      iv = dcov.var;  /* covkind == INVDIAGC */      sum = 0.0;      for (i=1; i<=vSize; i++){         x = v1[i]-v2[i]; sum += x*x*iv[i];      }      break;   case FULLC: