binaries/src/ViennaRNA/Cluster/cluster.c

   1 /* Last changed Time-stamp: <95/07/12 19:52:41 ivo> */
   2 /*
   3                  Cluster Analysis using Ward's Method
   4                 Ward J Amer Stat Ass, 58 (1963), p236
   5                    c Peter Stadler and Ivo Hofacker
   6 */
   7 #include <stdio.h>
   8 #include <stdlib.h>
   9 #include <string.h>
  10 #include "utils.h"
  11
  12 #define PUBLIC
  13 #define PRIVATE static
  14
  15 #define INFINITY   1000000
  16
  17 typedef struct{
  18         int   set1;
  19         int   set2;
  20         float distance;
  21         float distance2;
  22         } Union;
  23
  24 typedef struct {
  25                  int  type;
  26                  int  weight;
  27                  int  father;
  28                  int  sons;
  29                  int  leftmostleaf;
  30                } Postorder_list;
  31
  32 PUBLIC Union *wards_cluster(float **clmat);
  33 PUBLIC Union *neighbour_joining(float **clmat);
  34 PUBLIC void   printf_phylogeny(Union *tree, char *type);
  35
  36
  37 /*--------------------------------------------------------------------*/
  38
  39 PUBLIC Union *wards_cluster(float **clmat)
  40 {
  41    float  **d;
  42    int      *indic;
  43    int      *size;
  44    float    *help;
  45    Union    *tree;
  46
  47    float    min,deno,xa,xb,x;
  48    int      i,j,step,s=0,t=0,n;
  49
  50    n= (int)(clmat[0][0]);
  51
  52    size  = (int *)     space((n+1)*sizeof(int));
  53    d     = (float **)  space((n+1)*sizeof(float *));
  54    for(i=0;i<=n;i++)
  55       d[i] = (float *) space((n+1)*sizeof(float));
  56    indic = (int *)     space((n+1)*sizeof(int));
  57    help  = (float *)   space((n+1)*sizeof(float));
  58    tree  = (Union *)   space((n+1)*sizeof(Union));
  59
  60    tree[0].set1      = n;
  61    tree[0].set2      = 0;
  62    tree[0].distance  = 0.0;
  63    tree[0].distance2 = 0.0;
  64
  65    for (i=1;i<=n;i++) size[i]=1;
  66    for (i=1; i<=n; i++){
  67       for(j=1; j<=n; j++){
  68          d[i][j] = clmat[i][j];
  69       }
  70    }
  71
  72    /* look for the indices [s,t]  with minimum  d[s][t]*/
  73    for(step=1;step<n; step++){
  74       min = INFINITY;
  75       for (i=1; i<=n; i++){
  76          if (indic[i]==0){
  77             for (j=1; j<=n; j++){
  78                if(j!=i){
  79                   if (indic[j]==0){
  80                      if(d[i][j] < min) {
  81                         min = d[i][j];
  82                         s = i;
  83                         t = j;
  84                      }
  85                   }
  86                }
  87             }
  88          }
  89       }
  90
  91       /* now we have to join the clusters s and t and update the working array*/
  92
  93       tree[step].set1     = s;
  94       tree[step].set2     = t;
  95       tree[step].distance = min;
  96       indic[t] =1;
  97
  98       for (i=1; i<=n; i++){
  99          if (indic[i]==0){
 100             deno = (float) (size[i]+size[s]+size[t]);
 101             xa = ((float) (size[i]+size[s]))/deno;
 102             xb = ((float) (size[i]+size[t]))/deno;
 103              x = ((float) size[i])/deno;
 104             help[i] = xa*d[i][s] + xb*d[i][t] - x*d[s][t];
 105          }
 106       }
 107       for (i=1; i<=n; i++){
 108          if (indic[i]==0){
 109             d[s][i] = help[i];
 110             d[i][s] = help[i];
 111          }
 112       }
 113       size[s] += size[t];
 114    }
 115
 116    free(help);
 117    free(indic);
 118    for(i=0;i<=n;i++) free(d[i]);
 119    free(d);
 120    free(size);
 121
 122    return tree;
 123 }
 124
 125 /*--------------------------------------------------------------------*/
 126
 127 PUBLIC Union *neighbour_joining(float **clmat)
 128 {
 129   int n,i,j,k,l,step,ll[3];
 130   float b1,b2,b3,nn,tot,tmin,d1,d2;
 131   int mini=0, minj=0;
 132   int    *indic;
 133   float  *av, *temp;
 134   float **d;
 135   Union   *tree;
 136
 137   n = (int) (clmat[0][0]);
 138
 139   tree = (Union *) space((n+1)*sizeof(Union));
 140   indic = (int   *) space((n+1)*sizeof(int)   );
 141   av   = (float *) space((n+1)*sizeof(float) );
 142   temp = (float *) space((n+1)*sizeof(float) );
 143   d    = (float**) space((n+1)*sizeof(float*));
 144   for(i=0;i<=n;i++) d[i] = (float*) space((n+1)*sizeof(float));
 145
 146
 147   for (i=1; i<=n; i++){
 148      for(j=1; j<=n; j++){
 149         d[i][j] = clmat[i][j];
 150      }
 151      d[i][i]=0.0;
 152   }
 153
 154   tree[0].set1      = n;
 155   tree[0].set2      = 0;
 156   tree[0].distance  = 0.0;
 157   tree[0].distance2 = 0.0;
 158
 159   nn = (float) n;
 160
 161   for(step=1;step<=n-3;step++) {
 162 /*    for(j=2;j<=n;j++) for(i=1;i<j;i++) d[j][i]=d[i][j]; */
 163      tmin=99999.9;
 164      for(l=2; l<=n; l++){
 165         if(!indic[l]) {
 166            for(k=1;k<l;k++){
 167               if(!indic[k]) {
 168                  d1=0.0; d2=0.0;
 169                  for(i=1; i<=n;i++) {
 170                     d1 += d[i][k];
 171                     d2 += d[i][l];
 172                  }
 173                  tot=(nn-2.0)*d[k][l]-d1-d2;
 174                  if(tot<tmin){
 175                     tmin=tot;
 176                     mini=k;
 177                     minj=l;
 178                  }
 179               }
 180            }
 181         }
 182      }
 183
 184      d1=0.0; d2=0.0;
 185      for(i=1;i<=n;i++) {
 186         d1 +=d[i][mini];
 187         d2 +=d[i][minj];
 188      }
 189      d1 = (d1-d[mini][minj])/(nn-2.0);
 190      d2 = (d2-d[mini][minj])/(nn-2.0);
 191
 192      tree[step].set1      = mini;
 193      tree[step].distance  = (d[mini][minj]+d1-d2)*0.5-av[mini];
 194      tree[step].set2      = minj;
 195      tree[step].distance2 = d[mini][minj]-(d[mini][minj]+d1-d2)*0.5-av[minj];
 196
 197      av[mini]=d[mini][minj]*0.5;
 198
 199      nn=nn-1.0;
 200      indic[minj]=1;
 201      for(j=1;j<=n;j++) {
 202         if(!indic[j])
 203            temp[j]=(d[mini][j]+d[minj][j])*0.5;
 204      }
 205      for(j=1;j<=n;j++) {
 206         if(!indic[j]) {
 207            d[mini][j] = temp[j];
 208            d[j][mini] = temp[j];
 209         }
 210      }
 211      for(j=1;j<=n;j++){
 212          d[minj][j] = 0.0;
 213          d[j][minj] = 0.0;
 214          d[j][j]    = 0.0;
 215      }
 216   }
 217
 218   j=0;
 219   for(i=1;i<=n;i++) {
 220      if(!indic[i]){
 221         ll[j]=i;
 222         j++;
 223      }
 224   }
 225   b1=(d[ll[0]][ll[1]]+d[ll[0]][ll[2]]-d[ll[1]][ll[2]])*0.5;
 226   b2=d[ll[0]][ll[1]]-b1;
 227   b3=d[ll[0]][ll[2]]-b1;
 228   b1 -= av[ll[0]];
 229   b2 -= av[ll[1]];
 230   b3 -= av[ll[2]];
 231   tree[step].set1      = ll[1];
 232   tree[step].distance  = b2;
 233   tree[step].set2      = ll[2];
 234   tree[step].distance2 = b3;
 235   step++;
 236   tree[step].set1      = ll[0];
 237   tree[step].distance  = 0.0;
 238   tree[step].set2      = ll[1];
 239   tree[step].distance2 = b1;
 240
 241   for(i=0;i<=n;i++) free(d[i]); free(d);
 242   free(temp);
 243   free(av);
 244   free(indic);
 245
 246   return tree;
 247 }
 248
 249 /*--------------------------------------------------------------------*/
 250
 251
 252 PUBLIC void   printf_phylogeny(Union *tree, char *type)
 253 {
 254    int i,n;
 255    n = tree[0].set1;
 256
 257    printf("> %d %s ( Phylogeny using ",n, type);
 258    switch(type[0]){
 259     case 'W'  :
 260       printf("Ward's Method )\n");
 261       printf("> Nodes      Variance\n");
 262       for(i=1; i<n; i++)
 263          printf("%3d %3d    %9.4f \n",
 264          tree[i].set1, tree[i].set2, tree[i].distance);
 265       break;
 266     case 'N' :
 267       printf("Saitou's Neighbour Joining Method )\n");
 268       printf("> Nodes      Branch Length in Tree\n");
 269       for(i=1; i<n; i++)
 270          printf("%3d %3d    %9.4f   %9.4f\n",
 271          tree[i].set1, tree[i].set2, tree[i].distance, tree[i].distance2);
 272       break;
 273     default:
 274       printf(" non-identified Method )\n");
 275       for(i=1; i<n; i++)
 276          printf("%3d %3d    %9.4f   %9.4f\n",
 277          tree[i].set1, tree[i].set2, tree[i].distance, tree[i].distance2);
 278    }
 279 }
 280
 281 /*--------------------------------------------------------------------*/