ic = dcov.inv; /* covkind == FULLC */      for (i=1;i<=vSize;i++)         vTmp[i] = v1[i] - v2[i];      sum = 0.0;      for (i=2;i<=vSize;i++) {         crow = ic[i];         for (j=1; j<i; j++)            sum += vTmp[i]*vTmp[j]*crow[j];      }      sum *= 2;      for (i=1;i<=vSize;i++)         sum += vTmp[i] * vTmp[i] * ic[i][i];      break;   default:      HError(7170,"Distance: bad cov kind %d",dck);   }   if (trace&T_CDI) {      ShowVector("   dvec 1",v1,10);      ShowVector("   dvec 2",v2,10);      printf("   distance = %f\n",sqrt(sum));   }   return sqrt(sum);}/* AllocateVectors: distribute all pool vectors amongst clusters   1..curNumCl.  The centres of these clusters have already set.  Each   vector is placed in cluster with nearest centre.  Also, sets   aveCost of each cluster and returns totalCost.  Index of any cluster   with < minClustSize elements is returned, 0 otherwise */static int AllocateVectors(float *totalCost){   int n, i, bestn, cs;   float d,min;   Vector v;   if (trace&T_CAL)      printf("  allocating pool amongst %d clusters\n",curNumCl);   /* zero all clusters */   *totalCost=0.0;   for (n=1; n<=curNumCl; n++){      ccs->cl[n].aveCost = 0.0;      ccs->cl[n].csize = 0;   }   /* scan pool of vectors */   for (i=1; i<=nItems; i++) {      v = (Vector)GetItem(cvp,i);      /* find centre nearest to i'th vector */       min=Distance(v,ccs->cl[1].vCtr); bestn = 1;      for (n=2; n<=curNumCl; n++) {         d = Distance(v,ccs->cl[n].vCtr);         if (d < min) {            min = d; bestn = n;         }      }      if (trace&T_CAL)         printf("   item %d -> cluster %d, cost = %f\n",i,bestn,min);      /* increment costs and allocate vector to bestn */      *totalCost += min;         ccs->cl[bestn].aveCost += min;      cmap[i] = bestn; ++ccs->cl[bestn].csize;   }   /* Check for any empty clusters and average costs */   for (n=1; n<=curNumCl; n++) {      cs = ccs->cl[n].csize;      if (cs < minClustSize) {         if (trace&T_CAL)            printf("   cluster %d empty\n",n);         return n;      }      ccs->cl[n].aveCost /= cs;   }   return 0;}/* SplitVectors: distribute all pool vectors in cluster n between   clusters n1 and n2.  The centres of these clusters have already   set.  Each vector is placed in cluster with nearest centre.   */static void SplitVectors(int n, int n1, int n2){   int i, bestn;   float d1,d2;   Vector v;   int c,c1,c2;   if (trace&T_CAL)      printf("  SplitVectors: ");   for (i=1; i<=nItems; i++)      if (cmap[i]==n) {         v = (Vector)GetItem(cvp,i);         /* find centre nearest to i'th vector */          d1=Distance(v,ccs->cl[n1].vCtr);          d2=Distance(v,ccs->cl[n2].vCtr);          bestn = (d1<d2)?n1:n2;         /* allocate vector to bestn */         cmap[i] = bestn; ++ccs->cl[bestn].csize;      }   /* Check for any empty clusters and average costs */   c = ccs->cl[n].csize;   c1 = ccs->cl[n1].csize;   c2 = ccs->cl[n2].csize;   if (trace&T_CAL)      printf("  clusters split %d[%d] ->%d[%d] + %d[%d]\n",             n,c,n1,c1,n2,c2);      if (c1 == 0 || c2 == 0)      HError(7120,"SplitVectors: empty cluster %d[%d] ->%d[%d] + %d[%d]",             n,c,n1,c1,n2,c2);}/* FindCentres: of clusters a..b */static void FindCentres(int a, int b){   int n,i,j,cidx;   Vector v,ctr;   float cs;      for (n=a; n<=b; n++)      ZeroVector(ccs->cl[n].vCtr);   for (i=1; i<=nItems; i++) {      cidx = cmap[i];      if (cidx>=a && cidx<=b){         v = (Vector)GetItem(cvp,i);         ctr = ccs->cl[cidx].vCtr;         for (j=1; j<=vSize; j++)            ctr[j] += v[j];      }   }   for (n=a; n<=b; n++){      ctr = ccs->cl[n].vCtr;      cs = ccs->cl[n].csize;      for (j=1; j<=vSize; j++) ctr[j] /= cs;   }}/* FindCovariance: of cluster n */static void FindCovariance(int n){   Vector v,mean;   TriMat t;   int i,j,k;   double nx,x,y;   DVector sqsum;   DMatrix xsum;   Matrix c;      nx = ccs->cl[n].csize;   mean = ccs->cl[n].vCtr;   switch(ccs->ck){   case NULLC:      ccs->cl[n].cov.var = NULL;      break;   case DIAGC:       /* diagonal covariance matrix */   case INVDIAGC:    /* inverse diag covariance matrix */      sqsum = CreateDVector(&gstack,vSize);      ZeroDVector(sqsum);      for (i=1; i<=nItems; i++) {         if (cmap[i] == n) {            v = (Vector) GetItem(cvp,i);            for (j=1; j<=vSize; j++){               x = v[j]-mean[j];               sqsum[j] += x*x;            }         }      }      v = ccs->cl[n].cov.var;      for (j=1; j<=vSize; j++)         v[j] = (ccs->ck==DIAGC)?sqsum[j]/nx:nx/sqsum[j];      FreeDVector(&gstack,sqsum);      break;   case FULLC:    /* inverse full covariance matrix */      xsum = CreateDMatrix(&gstack,vSize,vSize);      c = CreateMatrix(&gstack,vSize,vSize);      ZeroDMatrix(xsum);      for (i=1; i<=nItems; i++) {         if (cmap[i] == n) {            v = (Vector) GetItem(cvp,i);            for (j=1; j<=vSize; j++)               for (k=1; k<=j; k++){                  x = v[j]-mean[j];                  y = v[k]-mean[k];                  xsum[j][k] += x*y;               }         }      }      t = ccs->cl[n].cov.inv;      for (j=1; j<=vSize; j++)         for (k=1; k<=j; k++)            t[j][k] = xsum[j][k]/nx;      CovInvert(t,c); Mat2Tri(c,t);      FreeDMatrix(&gstack,xsum);      break;   default:      HError(7170,"FindCovariance: unsupported cov kind [%d]",ccs->ck);   }}/* BiggestCluster: return index of cluster with highest average cost */static int BiggestCluster(void){   float maxCost;   int n,biggest;      maxCost=ccs->cl[1].aveCost; biggest = 1;   for (n=2; n<=curNumCl; n++)      if (ccs->cl[n].aveCost > maxCost) {         maxCost=ccs->cl[n].aveCost;         if( ccs->cl[n].csize >= minClustSize )            biggest=n;      }   return biggest;}/* FullestCluster: return index of cluster with most elements */static int FullestCluster(void){   int max,n,fullest;      max=ccs->cl[1].csize; fullest = 1;   for (n=1; n<=curNumCl; n++)      if (ccs->cl[n].csize > max) {         max=ccs->cl[n].csize; fullest=n;      }   return fullest;}/* Perturb: copy perturbed versions of cluster n into n1 and n2. */static void Perturb(int n, int n1, int n2){   int i;   Vector v,v1,v2;   float x;      if (trace&T_CGE)      printf("  Perturb: cluster %d -> %d + %d\n",n,n1,n2);   v  = ccs->cl[n].vCtr;   v1 = ccs->cl[n1].vCtr;   v2 = ccs->cl[n2].vCtr;   if (n!=n1) CopyVector(v,v1);   if (n!=n2) CopyVector(v,v2);     for (i=1; i<=vSize; i++) {      x = fabs(v[i]*0.01);      if (x<0.0001) x=0.0001;      v1[i] += x; v2[i] -= x;   }}/* NumTreeNodes: check that nc is power of 2, return 2nc-1 */static int NumTreeNodes(int nc){   int n = nc;      while (n>2) {      if ((n % 2) == 1)         HError(7172,"NumTreeNodes: nc %d not a power of 2",nc);      n /= 2;   }   return 2*nc-1;}/* InitClustering: create the cluster set data structure */static void InitClustering(MemHeap *x, Sequence vpool, int nc,                           Boolean treeCluster, CovKind distck, CovKind clusck, Covariance distcov){   int i,numClust;   Vector v;   Covariance cov;   nItems = vpool->nItems;   if (nItems < nc)      HError(7120,"InitClustering: only %d items for %d clusters",nItems,nc);   cvp = vpool;      numClust = (treeCluster)?NumTreeNodes(nc):nc;   curNumCl = 1;   ccs = (ClusterSet *)New(x,sizeof(ClusterSet));   ccs->isTree = treeCluster;   ccs->numClust = numClust;   ccs->ck = clusck;   ccs->cl = (Cluster *)New(x,sizeof(Cluster)*numClust);   --ccs->cl;  /* index is 1..numClust */   cmap = CreateShortVec(&gstack,nItems);   for (i=1; i<=nItems; i++)  /* put all vecs in cluster 1 */      cmap[i] = 1;   v = (Vector)GetItem(cvp,1); vSize = VectorSize(v);   vTmp = CreateVector(&gstack,vSize);   dck = distck; dcov = distcov;   for (i=1; i<=numClust; i++){      ccs->cl[i].csize = 0;      ccs->cl[i].vCtr = CreateVector(x,vSize);      ccs->cl[i].aveCost = 0.0;      switch (clusck){      case NULLC:         cov.var = NULL;         break;      case DIAGC:      case INVDIAGC:         cov.var = CreateVector(x,vSize);         break;      case FULLC:         cov.inv = CreateTriMat(x,vSize);         break;      default:         HError(7170,"InitClustering: bad cluster ckind");      }      ccs->cl[i].cov = cov;   }   ccs->cl[1].csize = nItems;   ccs->cl[1].aveCost = 1.0;   FindCentres(1,1);}/* EXPORT->FlatCluster: apply linear clustering to vpool */ClusterSet *FlatCluster(MemHeap *x, Sequence vpool, int nc,                         CovKind dck, CovKind cck, Covariance dcov){   int i,c,cc,ce,iter,repairCount;   float oldCost,newCost;   Boolean converged;      InitClustering(x,vpool,nc,FALSE,dck,cck,dcov);   if (trace&T_CGE)      printf("FlatCluster: %d items -> %d clusters\n",             ccs->cl[1].csize,ccs->numClust);   for (c=2; c<=ccs->numClust; c++) {      if (trace&T_CGE)         printf(" increasing from %d -> %d clusters\n",c-1,c);      /* increase num clusters by splitting biggest */      cc = BiggestCluster();      Perturb(cc,cc,c); curNumCl = c;      oldCost = 1e10;      /* reallocate vectors until cost stabilises */      for (iter=0,converged=FALSE; !converged; iter++){         repairCount = 0; /* try to fill empty clusters by splitting fullest */         while ((ce=AllocateVectors(&newCost)) != 0  && ++repairCount <=c){            cc = FullestCluster();            Perturb(cc,cc,ce);      /* ce = empty cluster */         }         if (ce != 0)            HError(7120,"FlatCluster: Failed to make %d clusters at iter %d\n",c,iter);         if (trace & T_CLC){            printf("   c=%d, iter=%d, cost = %e\n",c,iter,newCost);            fflush(stdout);         }         converged = (iter>=maxIter) || ((oldCost-newCost) / oldCost < 0.001);         FindCentres(1,curNumCl); oldCost = newCost;      }      if (trace & T_DCM) DumpClusterMap();   }   for (i=1; i<=ccs->numClust; i++)      FindCovariance(i);   FreeShortVec(&gstack,cmap);   return ccs;}  /* EXPORT->TreeCluster: apply binary tree clustering to vpool */ClusterSet *TreeCluster(MemHeap *x, Sequence vpool, int nc,                         CovKind dck, CovKind cck, Covariance dcov){   int rowsize,c,i,j;      InitClustering(x,vpool,nc,TRUE,dck,cck,dcov);   if (trace&T_CGE)      printf("TreeCluster: %d items -> %d clusters[%d nodes]\n",             ccs->cl[1].csize,nc,ccs->numClust);   FindCovariance(1);   c = 1; rowsize = 1;   while (rowsize < nc){      for (i=c,j=c+rowsize; i<c+rowsize; i++, j+=2){         Perturb(i,j,j+1);         SplitVectors(i,j,j+1);         FindCentres(j,j+1);         FindCovariance(j);         FindCovariance(j+1);      }      c += rowsize;      rowsize *= 2;   }   FreeShortVec(&gstack,cmap);   return ccs;}/* EXPORT->FreeClusterSet: free storage used by cs */void FreeClusterSet(ClusterSet *cs){   Dispose(cs->x,cs);}/* EXPORT->ShowClusterSet: print the given clusterset */void ShowClusterSet(ClusterSet *cs){   int i;   Cluster *c;   printf("%s Cluster Set: %d nodes\n",          cs->isTree?"Tree":"Flat", cs->numClust);   for (i=1,c=cs->cl+1; i<=cs->numClust; i++,c++) {      printf("%d. size=%d",i,c->csize);      if (!cs->isTree) printf(" cost=%.3f",c->aveCost);      printf("\n");      ShowVector("  mean",c->vCtr,20);      switch(cs->ck){      case NULLC:          break;      case INVDIAGC:          ShowVector("  invdiagC",c->cov.var,20);         break;      case DIAGC:          ShowVector("  diagC",c->cov.var,20);         break;      case FULLC:         ShowTriMat("  fullC",c->cov.inv,20,20);         break;      }        }}/* ------------------------- Accumulators ---------------------- */static int muC,vaC,trC,wtC,prC;/* CreateMuAcc:  create an accumulator for means */static MuAcc *CreateMuAcc(MemHeap *x, int vSize, int nPara){   MuAcc *ma;   int count;   ma = (MuAcc *)New(x,sizeof(MuAcc)*nPara);   for(count=0;count<nPara;count++){     ma[count].mu = CreateVector(x,vSize);     ZeroVector(ma[count].mu);     ma[count].occ = 0.0;     ++muC;  }   return ma;}/* CreateVaAcc: create an accumulator for square sums */static VaAcc *CreateVaAcc(MemHeap *x, int vSize, CovKind ck, int nPara){   VaAcc *va;   int count;   va = (VaAcc *) New(x,sizeof(VaAcc)*nPara);   for(count=0;count<nPara;count++){     switch(ck){     case DIAGC:     case INVDIAGC:       va[count].cov.var = CreateVector(x,vSize);       ZeroVector(va[count].cov.var);       break;     case FULLC:       va[count].cov.inv = CreateTriMat(x,vSize);       ZeroTriMat(va[count].cov.inv);