/*  CLUSTER.SRC              copyright by R.Schlittgen, 4.April.1999
**                           ( e-mail: schlitt@hermes1.econ.uni-hamburg.de )
**
**           Clustering and dendrogram
**
**  AGGLOM     agglomerative clustering                               15
**  DENDRO     plotting a dendrogram                                 123
**  DENDROCO   determination of cooordinates fr DENDRO               186
**  GROUP      presents the clusters built up to a given distance    235
**  KMEANS     k-means clustering                                    283
*/

/*
**  AGGLOM
**  Purpose: agglomerative clustering
**
**  Format:  {o,lev} = AGGLOM(D,m)
**  INPUT :  D = (n,n) distance matrix
**           m = scalar, gives the clustering method:
**               1 : single lincage
**               2 : complete lincage
**               3 : average lincage
**               4 : centroid method
**
**  OUTPUT:  o = (n-1,n)-matrix
**                each row gives the newly built cluster
**           lev = (n-1,1) vector, the distances at which the clusters are
**                 glued together
**
*/


proc (2)= AGGLOM(D,method);

  local ap,aq,be,c,c0,dmin,d1,ga,i,index,j,n,ncd,np,nq,p,q,r,rs,s2,s1,w,w1,w2;


  if not D==D'; print "Distance matrix must be symmeric."; end; endif;
  if method==1;                           @ Single Linkage @
    ap = 0.5; aq = 0.5; be=0; ga=-0.5;
   elseif method == 2;                    @ Complete Linkage @
    ap = 0.5; aq = 0.5; be=0; ga=0.5;
   elseif method == 3;                    @ Average Linkage @
    ap = 0.5; aq = 0.5; be=0; ga=0;
   elseif method == 4;                    @ Centroid method @
    ap = 0.5; aq = 0.5; be=-0.25; ga=0;
   else;
    print "Only 4 methods are implemented (method=1,2,3,4).";stop;
  endif;
  n = rows(D);
  index = seqa(1,1,n)';
  ncd = {};
  c0 = zeros(n,1);
  c = zeros(n-1,n);
  i = 0;
  do while rows(D) >= 2;
    i = i+1;
    r = rows(D);
    rs=seqa(1,1,r);
    dmin = minc(minc(D + 1E32*eye(r)));
    D1 = (D+1E32*eye(r)).==dmin;
    p = selif(rs,maxc(D1));
    p = p[1];
    q = selif(rs,D1[.,p].==1);
    q = q[1];
    ncd = ncd|(index[p q]~dmin);
    if i ==1;
       c[i,1 2] = index[p q];
       np = 1;
       nq = 1;
    else;
     s1 = ceil(maxc(indexcat(vecr(c[1:i-1,.]),ncd[i,1]))/n);
     s2 = ceil(maxc(indexcat(vecr(c[1:i-1,.]),ncd[i,2]))/n);
     if scalmiss(s1) and scalmiss(s2);
        np = 1;
        nq = 1;
        c[i,1 2] = ncd[i,1 2];
      elseif scalmiss(s1) and (not scalmiss(s2));
        np = 1;
        w  = delif(c[s2,.]',c[s2,.]'.==c0);
        nq = rows(w);
        w = w|ncd[i,1]';
        c[i,1:rows(w)]=w';
      elseif scalmiss(s2) and (not scalmiss(s1));
        nq = 1;
        w = delif(c[s1,.]',c[s1,.]'.== c0);
        np = rows(w);
        w = w|ncd[i,2]';
        c[i,1:rows(w)]=w';
      else;
        w1 = delif(c[minc(s1|s2),.]',c[minc(s1|s2),.]'.==c0);
        np = rows(w1);
        w2 = delif(c[maxc(s1|s2),.]',c[maxc(s1|s2),.]'.==c0);
        nq = rows(w2);
        w = w1|w2;
        j = rows(w)+1;
        do while j > 2 ;
           j=j-1;
             if not scalmiss(indexcat(w[1:(j-1)],w[j]));
                w[j]=0;
             endif;
           endo;
         w = delif(w,w.==zeros(rows(w),1))';
         c[i,1:cols(w)]=w;
      endif;
    endif;

/* updating  ap,aq, be for centroid method */

    if method == 4;
      ap = np/(np+nq);
      aq = nq/(np+nq);
      be = -ap*aq;
    endif;

    D[.,p] = ap*D[.,p]+aq*D[.,q]+be*dmin+ga*abs(D[.,p]-D[.,q]);
    D[p,.] = D[.,p]';
    D[p,p] = 0;
    if q < r;
      D = D[1:(q-1) (q+1):r,1:(q-1) (q+1):r];
      index = index[1:(q-1) (q+1):r];
    else;
      D = D[1:(q-1),1:(q-1)];
      index = index[1:(q-1)];
    endif;
 endo;
 retp(c,ncd[.,3]);
 endp;


/*
**  DENDRO
**  Purpose: plots dendrogram
**
**  Format:  DENDRO(o,lev,lab,space)
**  INPUT :  o = (n-1,n) matrix, output of AGGLOM
**           lev = (n-1,1) vector, the distances at which the clusters
**                 are glued together (output of AGGLOM)
**           lab = (n,1) string array, the labels of objects, or " "
**           space, scalar used for positioning the labels for the objects
**
**  REMARK: this proc uses DENDROCO
**
*/


proc (0) = DENDRO(o,lev,lab,space);

  local i,t,x,y,n;

  graphset;
  n = rows(o)+1;
  {x,y} = dendroco(o,lev);
  t = rows(x);

  _pline =  ones(t,9);
  _pline[.,2] = _pline[.,2]*6;
  _pline[., 3:6] = x[.,1]~y[.,1]~x[.,2]~y[.,2];
  _pline[.,9] = _pline[.,9]*2;

  if rows(lab) /= n;
     lab = "Obj" $+ ftocv(o[n-1,.]',1,0);
    else;
     lab = lab[o[n-1,.]];
  endif;
   _pmsgstr = lab[1];
    i=1;
  do while i < n;
     i=i+1;
    _pmsgstr =_pmsgstr$+"\000"$+lab[i];
   endo;
   _pmsgctl = ones(n,1)*(-space~0~(3/n)~0~1~1~2);
   _pmsgctl[.,2] = seqa(n-.25,-1,n);

   _pdate="";
   _pframe={0,1};
   _pmcolor = zeros(8,1)|15;
   scale((-space)|maxc(lev),0|n);
   _paxes={1,0};
  xy(0|0,0|n);
 graphset;
 endp;



/*
**  DENDROCO
**  Purpose: determines the coordinates for DENDRO
**
**  Format:  {x,y} = DENDROCO(O,LEV)
**  INPUT :  o = (n-1,n) matrix, output of AGGLOM
**           LEV = (n-1,1) vector, output of AGGLOM
**  OUTPUT:  x = (t,2) matrix of x coordinates
**           y = (t,2) matrix of y coordinates
**
*/


 proc (2) = DENDROCO(o,lev);

 local  i,j,lab,m1,m2,n,tr,xcoord,x0,ycoord,ym,ind,y0;

 n = cols(o);
 tr = o;
 lab = tr[n-1,.]';
  x0 = zeros(n,1);
  y0 = seqa(n,-1,n);
  xcoord = {};
  ycoord = {};
  i=0;
    do while i<n-1;
      i = i+1;
      ind = indexcat(lab,tr[i,1]);
         j = 1;
          do while (tr[i,j] ne 0) and (j < n) ;
             j = j+1;
             ind = ind|indexcat(lab,tr[i,j]);
          endo;
          ind = packr(ind);
          m1 = minc(y0[ind]);
          m1 = meanc(selif(ind, y0[ind].==(m1*ones(rows(ind),1))));
          m2 = maxc(y0[ind]);
          m2 = meanc(selif(ind,y0[ind].==(m2*ones(rows(ind),1))));
          ym = m1|m2;
          xcoord = xcoord|((x0[ym]|lev[i])~(lev[i]*ones(rows(ym)+1,1)));
          ycoord = ycoord|((y0[ym]~y0[ym])|y0[minc(ym)|maxc(ym)]');
          x0[ind] = lev[i]*ones(rows(ind),1);
          y0[ind] = meanc(y0[ym])*ones(rows(ind),1);
      endo;
   retp(xcoord,ycoord);
 endp;



/*
**  GROUP
**  Purpose: gives clusters at a given level of distance
**
**  Format:  G = GROUP(o,lev,d)
**  INPUT :  o = (n-1,n) matrix, output of agglom
**           lev = (n-1,1) vector, level at which clusters are glued,
**                 output of agglom
**           d = scalar, the distance level
**
**  OUTPUT:  G = (k,n) matrix
**              each row represents one cluster;
**              rows are filled up with missings
**
*/

 proc (1) = GROUP(o,lev,d);

  local g,g1,i,index,j,m,n;

  n = cols(o);
  g = delif(o,lev.>d);
  if ismiss(g); print "All elements in atoms"; end; endif;

  g = rev(miss(g,0));
  j = 0;
  do while j < rows(g);
     j = j+1;
     g1 = sortc(packr(g[j,.]'),1);
     g[j,1:rows(g1)] = g1';
      m = rows(g);
      index = seqa(1,1,m);
      i = j;
      do while i < m;
        i = i+1;
        if scalmiss(setdif(packr(g[i,.]'),g1,1)); index[i]=0; endif;
      endo;
      g = delif(g,index.==zeros(m,1));
   endo;
   g1 = setdif(seqa(1,1,n),packr(vecr(g)),1);
   if not scalmiss(g1);
     g = g|(g1~miss(zeros(rows(g1),n-1),0));
   endif;
   g = sortc(g,1);
 retp(g);
 endp;


/*
**  KMEANS
**  PURPOSE: finds a given number of clusters at a local minimum
**           of within group variance
**
**  FORMAT:  {ind,c} = KMEANS(x,g)
**  INPUT :  x = (n,p) matrix, data
**           g = scalar, number of clusters
**
**  OUTPUT:  ind = (n,1) vector. each row indicats the cluster number
**                 the corresponding row of x belongs to
**           c   = (g,p) matrix, centers of clusters
**
*/

 proc (2) = KMEANS(x,g);

  local c,d,dist,h,i,ind,indold,iter,n;

   n = rows(x);
   indold=zeros(n,1);

   /* initial cluster centers */

   d = sumc((x'-meanc(x))^2);
   d = sortc(seqa(1,1,n)~d,2);
   d = d[1+seqa(0,1,g)*n/g,1];
   c = x[d,.]';

   /* main loop */

   iter = 0;
   do while iter < 100;
      iter = iter+1;
      dist = {};
      i = 0;
      do while i < g;
         i = i+1;
         dist = dist|sumc((x'-c[.,i])^2)';
      endo;
      ind = minindc(dist);
      if ind==indold; goto ende; endif;
      indold=ind;

      /* new centers */

      i = 0;
      do while i < g;
         i = i+1;
         h = selif(x,ind.==i*ones(n,1));
         if scalmiss(h); print "One cluster got emty"; end; endif;
         c[.,i] = meanc(h);
      endo;
   endo;
   ende:
   retp(ind,c');
   endp;



