binaries/src/ViennaRNA/Cluster/split.c

   1 /* Last changed Time-stamp: <95/07/12 19:49:53 ivo>*/
   2 /*
   3                Split Decomposition of Distance Matrices
   4              as described by H.J.Bandelt and A.W.M.Dress
   5                        Adv Math, 92 (1992) p47
   6                  c Peter Stadler and Ivo L. Hofacker
   7 */
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <string.h>
  11 #include "utils.h"
  12
  13 #define  PUBLIC
  14 #define  PRIVATE      static
  15 #define  DEBUG        0
  16 #define  SHORT_OUTPUT 1
  17
  18 #define  DINFTY       1.e32
  19 #define  ZERO         1.e-10
  20
  21
  22 typedef struct {
  23    short   *splitlist[2];
  24    int      splitsize;
  25    double   isolation_index; } Split;
  26
  27 PUBLIC Split *split_decomposition(float **dist);
  28 PUBLIC void free_Split(Split *x);
  29 PUBLIC void print_Split(Split *x);
  30 PUBLIC void sort_Split(Split *x);
  31
  32 PUBLIC Split *split_decomposition(float **dist)
  33 {
  34
  35    int      elm, n_of_splits;
  36    short   *full_slots;
  37    int      i,j,sp,new_sp,spp;
  38    int      i2,j1,j2,maxlen;
  39    int      x,y,z;
  40    double   alpha,beta,tmp;
  41    double   test1,test2;
  42    Split   *SD, *S;
  43    int number_of_points;
  44
  45    number_of_points = (int) dist[0][0];
  46
  47    /* Initialize */
  48    elm = 2;
  49    n_of_splits = 1;
  50    maxlen = (number_of_points+1)*(number_of_points+1);
  51    full_slots = (short *) space(2*(maxlen+1)*sizeof(short));
  52    full_slots[1]=1;
  53
  54    SD = space((maxlen+1)*sizeof(Split));
  55
  56    SD[0].splitlist[0]    = NULL;
  57    SD[0].splitlist[1]    = NULL;
  58    SD[0].splitsize       = 1;
  59    SD[0].isolation_index = 0.0;
  60    /**/
  61    SD[1].splitlist[0]    = (short *) space((number_of_points+1)*sizeof(short));
  62    SD[1].splitlist[0][0] = number_of_points;
  63    SD[1].splitlist[0][1] = 1;
  64    SD[1].splitlist[1]    = (short *) space((number_of_points+1)*sizeof(short));
  65    SD[1].splitlist[1][1] = 2;
  66    SD[1].splitsize       = 1;
  67    SD[1].isolation_index = dist[1][2];
  68
  69
  70    /* Iteration */
  71
  72    for( elm=3; elm <= number_of_points; elm++){
  73
  74       for (sp=1, spp=n_of_splits;sp<=n_of_splits; sp++){
  75
  76          SD[sp].splitlist[0][0] = elm;
  77
  78          /* 1.  Split = {old_A | old_B+{elm}}  */
  79
  80          alpha = 2.0*SD[sp].isolation_index;
  81          for(i2=1; i2<=(elm - SD[sp].splitsize); i2++){
  82             if(i2==(elm-SD[sp].splitsize)) {
  83                x = elm;
  84             }
  85             else{
  86                x= SD[sp].splitlist[1][i2];
  87             }
  88             for(j1=1; j1 <= SD[sp].splitsize; j1++){
  89                y= SD[sp].splitlist[0][j1];
  90                for(j2=1; j2 <= (SD[sp].splitsize); j2++){
  91                   z= SD[sp].splitlist[0][j2];
  92
  93                   /* calculate the value beta = beta(elm,x; y,z) */
  94
  95                   beta = dist[elm][y] + dist[x][z];
  96                   tmp  = dist[elm][z] + dist[x][y];
  97                   if(tmp>beta) beta=tmp;
  98                   tmp  = dist[elm][x] + dist[y][z];
  99                   if(tmp>beta) beta=tmp;
 100                   beta -= ( dist[elm][x] + dist[y][z] );
 101
 102                   if(beta<alpha) alpha=beta;
 103                }
 104             }
 105          }
 106          alpha/=2.;
 107
 108          if (alpha > ZERO){
 109             /* add {old_A | old_B + {elm}} to the split-list */
 110             spp++;
 111             full_slots[spp] =1;
 112             SD[spp].splitsize       = SD[sp].splitsize;
 113             SD[spp].isolation_index = alpha;
 114             SD[spp].splitlist[0]    = (short *) space((number_of_points+1)*sizeof(short));
 115             SD[spp].splitlist[0][0] = elm;
 116             SD[spp].splitlist[1]    = (short *) space((number_of_points+1)*sizeof(short));
 117             for(i=1; i<=SD[spp].splitsize; i++)
 118                SD[spp].splitlist[0][i] = SD[sp].splitlist[0][i];
 119             for(i=1; i<=(elm-1-SD[spp].splitsize); i++)
 120                SD[spp].splitlist[1][i] = SD[sp].splitlist[1][i];
 121             SD[spp].splitlist[1][elm-SD[spp].splitsize] = elm;
 122          }
 123
 124          /* 2. Split = {old_A+{elm} | old_B}   */
 125
 126          alpha = 2.0*SD[sp].isolation_index;
 127          for(i2=1; i2<= (SD[sp].splitsize+1); i2++){
 128             if(i2==(SD[sp].splitsize+1)){
 129                x= elm;
 130             }
 131             else
 132                {
 133                   x= SD[sp].splitlist[0][i2];
 134                }
 135             for(j1=1; j1 <= (elm-1-SD[sp].splitsize); j1++){
 136                y= SD[sp].splitlist[1][j1];
 137                for(j2=1; j2 <= (elm-1-SD[sp].splitsize); j2++){
 138                   z= SD[sp].splitlist[1][j2];
 139
 140                   /* calculate the value beta = beta(elm,x; y,z) */
 141
 142                   beta = dist[elm][y] + dist[x][z];
 143                   tmp  = dist[elm][z] + dist[x][y];
 144                   if(tmp>beta) beta=tmp;
 145                   tmp  = dist[elm][x] + dist[y][z];
 146                   if(tmp>beta) beta=tmp;
 147                   beta -= ( dist[elm][x] + dist[y][z] );
 148
 149                   if(beta<alpha) alpha=beta;
 150                }
 151             }
 152          }
 153          alpha/=2.;
 154
 155          if (alpha > ZERO){
 156             /* replace {old_A | old_B} by {old_A+{elm} | old_B} in the splitlist */
 157             SD[sp].splitsize++;
 158             SD[sp].splitlist[0][SD[sp].splitsize] = elm;
 159             SD[sp].isolation_index = alpha;
 160          }
 161          else{
 162             /* remove {old_A | old_B} from the split list */
 163             full_slots[sp] = 0;
 164             free(SD[sp].splitlist[0]);
 165             SD[sp].splitlist[0]    = NULL;
 166             free(SD[sp].splitlist[1]);
 167             SD[sp].splitlist[1]    = NULL;
 168             SD[sp].splitsize       = 0;
 169             SD[sp].isolation_index = 0.0;
 170          }
 171       }
 172
 173       /* 3.  Split =  { {elm} | {1,...,elm-1} } */
 174
 175       alpha=DINFTY;
 176       for(i=1;i<=elm-1;i++){
 177          for(j=1; j<=elm-1;j++){
 178             tmp = dist[elm][i]+dist[elm][j] - dist[i][j];
 179             if( tmp < alpha) alpha = tmp;
 180          }
 181       }
 182       alpha/=2.;
 183       if (alpha > ZERO){
 184          spp++;
 185          full_slots[spp] = 1;
 186          SD[spp].splitsize         = 1;
 187          SD[spp].isolation_index   = alpha;
 188          SD[spp].splitlist[0]      = (short *) space((number_of_points+1)*sizeof(short));
 189          SD[spp].splitlist[0][0]   = elm;
 190          SD[spp].splitlist[1]      = (short *) space((number_of_points+1)*sizeof(short));
 191          SD[spp].splitlist[0][1]   = elm;
 192          for (i=1; i<= elm-1; i++)
 193             SD[spp].splitlist[1][i] = i;
 194       }
 195       /* note that spp now points to the last entry in the splitlist !!! */
 196
 197       /* garbage collection .... unfortunately a bit boring */
 198
 199       for ( sp=1, new_sp=0; sp<= spp; sp++ ){
 200          if (full_slots[sp]){
 201             new_sp++;
 202             if(sp != new_sp){
 203                /* copy all the junk from sp to new_sp */
 204                full_slots[new_sp]=1; full_slots[sp]=0;
 205                SD[new_sp].splitsize        = SD[sp].splitsize;
 206                SD[new_sp].isolation_index  = SD[sp].isolation_index;
 207                if(SD[new_sp].splitlist[0]==NULL)
 208                   SD[new_sp].splitlist[0] = (short *) space((number_of_points+1)*sizeof(short));
 209                if(SD[new_sp].splitlist[1]==NULL)
 210                   SD[new_sp].splitlist[1] = (short *) space((number_of_points+1)*sizeof(short));
 211                for(i=0;i<=elm;i++){
 212                   SD[new_sp].splitlist[0][i]  = SD[sp].splitlist[0][i];
 213                   SD[new_sp].splitlist[1][i]  = SD[sp].splitlist[1][i];
 214                }
 215                free(SD[sp].splitlist[0]);
 216                SD[sp].splitlist[0] = NULL;
 217                free(SD[sp].splitlist[1]);
 218                SD[sp].splitlist[1] = NULL;
 219                SD[sp].splitsize       = 0;
 220                SD[sp].isolation_index = 0.0;
 221             }
 222          }
 223       }
 224       n_of_splits = new_sp;
 225       SD[0].splitsize = n_of_splits;
 226 #if DEBUG
 227       for(test1=0, i=2; i<=elm; i++) for( j=1; j<i; j++) test1+=dist[i][j];
 228       for(test2=0, i=1; i<= n_of_splits; i++)
 229          test2 += (elm-SD[i].splitsize)*SD[i].splitsize*SD[i].isolation_index;
 230       SD[0].isolation_index = (test1 - test2)/test1;
 231       print_Split(SD);
 232 #endif
 233    } /* End of iteration */
 234
 235
 236    /* Calculate fraction of split-prime part for the full matrix
 237       if not already done */
 238
 239    for(test1=0.0, i=2; i<=number_of_points; i++)
 240       for( j=1; j<i; j++) test1+=dist[i][j];
 241    for(test2=0.0, i=1; i<= n_of_splits; i++)
 242       test2 += (number_of_points-SD[i].splitsize)*
 243          SD[i].splitsize*SD[i].isolation_index;
 244    SD[0].isolation_index = (test1 - test2)/test1;
 245
 246 #if DEBUG
 247    print_Split(SD);
 248 #endif
 249    /* free superfluous space */
 250
 251    S = space((n_of_splits+1)*sizeof(Split));
 252    for (i=0; i<= n_of_splits; i++) S[i] = SD[i];
 253    free(SD);
 254
 255    if(DEBUG)
 256       print_Split(S);
 257
 258    return(S);
 259 }
 260
 261 /* -------------------------------------------------------------------------- */
 262
 263 PUBLIC void free_Split(Split *x)
 264 {
 265    int i;
 266    for (i=0; i<=x[0].splitsize; i++) {
 267       if( x[i].splitlist[0] != NULL ) free( x[i].splitlist[0] );
 268       if( x[i].splitlist[1] != NULL ) free( x[i].splitlist[1] );
 269    }
 270    free(x);
 271 }
 272
 273 /* -------------------------------------------------------------------------- */
 274
 275 PUBLIC void print_Split(Split *x)
 276 {
 277    int j,k;
 278    printf("> %d Split Decomposition",x[0].splitsize);
 279    for (k=1; k<=x[0].splitsize; k++) {
 280       if( (x[k].splitlist[0])&&(x[k].splitlist[1]) ) {
 281          printf("\n%3d   %8.4f  : {",k,x[k].isolation_index);
 282          for(j=1;j<=x[k].splitsize;j++)
 283             printf(" %3d",x[k].splitlist[0][j]);
 284          printf("   |");
 285 #if SHORT_OUTPUT
 286          printf(" ...");
 287 #else
 288          for(j=1;j<=(x[k].splitlist[0][0]-x[k].splitsize);j++)
 289             printf(" %3d",x[k].splitlist[1][j]);
 290          printf(" } ");
 291 #endif
 292       }
 293    }
 294    printf("\n      %8.4f  : { [Split prime fraction] }\n", x[0].isolation_index);
 295 }
 296
 297
 298 /* -------------------------------------------------------------------------- */
 299
 300 int CompareSplit(const void *v1 , const void *v2)
 301 {
 302    Split *S1, *S2;
 303
 304    S1 = (Split *) v1; S2 = (Split *) v2;
 305    if(S1->isolation_index < S2->isolation_index) return +1;
 306    if(S2->isolation_index < S1->isolation_index) return -1;
 307    return 0;
 308 }
 309
 310 /* -------------------------------------------------------------------------- */
 311
 312 PUBLIC void sort_Split(Split *x)
 313 {
 314    int i,j;
 315    int t1;
 316
 317    qsort((void *) &(x[1]), x[0].splitsize ,sizeof(Split),CompareSplit );
 318
 319    for(i=1; i<=x[0].splitsize; i++) {
 320       if(x[i].splitsize > x[i].splitlist[0][0]-x[i].splitsize) {
 321          for(j=1; j<=x[i].splitsize; j++) {
 322             t1 = x[i].splitlist[0][j];
 323             x[i].splitlist[0][j] = x[i].splitlist[1][j];
 324             x[i].splitlist[1][j] = t1;
 325          }
 326          x[i].splitsize = x[i].splitlist[0][0]-x[i].splitsize;
 327       }
 328    }
 329 }
 330
 331
 332 /* -------------------------------------------------------------------------- */