src/bNNI.c

   1 /*#include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include "graph.h"
   5 #include "main.h"
   6 */
   7 #include "me.h"
   8
   9 /*boolean leaf(node *v);
  10 edge *siblingEdge(edge *e);
  11 edge *depthFirstTraverse(tree *T, edge *e);
  12 edge *findBottomLeft(edge *e);
  13 edge *topFirstTraverse(tree *T, edge *e);
  14 edge *moveUpRight(edge *e);*/
  15
  16 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger);
  17 void assignBMEWeights(tree *T, double **A);
  18 //void updateAveragesMatrix(tree *T, double **A, node *v,int direction);
  19 void bNNItopSwitch(tree *T, edge *e, int direction, double **A);
  20 int bNNIEdgeTest(edge *e, tree *T, double **A, double *weight);
  21 void updatePair(double **A, edge *nearEdge, edge *farEdge, node *closer, node *further, double dcoeff, int direction);
  22
  23 int *initPerm(int size);
  24
  25 void reHeapElement(int *p, int *q, double *v, int length, int i);
  26 void pushHeap(int *p, int *q, double *v, int length, int i);
  27 void popHeap(int *p, int *q, double *v, int length, int i);
  28
  29
  30 void bNNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray,
  31                 double *weights, int *location, int *possibleSwaps)
  32 {
  33   int tloc;
  34   tloc = location[e->head->index+1];
  35   location[e->head->index+1] =
  36     bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
  37   if (NONE == location[e->head->index+1])
  38     {
  39       if (NONE != tloc)
  40         popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
  41     }
  42   else
  43     {
  44       if (NONE == tloc)
  45         pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
  46       else
  47         reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
  48     }
  49 }
  50
  51 int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
  52
  53 void permInverse(int *p, int *q, int length);
  54
  55 void weighTree(tree *T)
  56 {
  57   edge *e;
  58   T->weight = 0;
  59   for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
  60     T->weight += e->distance;
  61 }
  62
  63 //void bNNI(tree *T, double **avgDistArray, int *count)
  64 void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
  65 {
  66   edge *e, *centerEdge;
  67   edge **edgeArray;
  68   int *p, *location, *q;
  69   int i,j;
  70   int possibleSwaps;
  71   double *weights;
  72   p = initPerm(T->size+1);
  73   q = initPerm(T->size+1);
  74   edgeArray = (edge **) malloc((T->size+1)*sizeof(double));
  75   weights = (double *) malloc((T->size+1)*sizeof(double));
  76   location = (int *) malloc((T->size+1)*sizeof(int));
  77
  78   double epsilon = 0.0;
  79   for (i=0; i<numSpecies; i++)
  80     for (j=0; j<numSpecies; j++)
  81       epsilon += D[i][j];
  82   epsilon = (epsilon / (numSpecies * numSpecies)) * EPSILON;
  83
  84   for (i=0;i<T->size+1;i++)
  85     {
  86       weights[i] = 0.0;
  87       location[i] = NONE;
  88     }
  89 /*  if (verbose)
  90     {
  91       assignBMEWeights(T,avgDistArray);
  92       weighTree(T);
  93     }*/
  94   e = findBottomLeft(T->root->leftEdge);
  95   while (NULL != e)
  96     {
  97       edgeArray[e->head->index+1] = e;
  98       location[e->head->index+1] =
  99         bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
 100       e = depthFirstTraverse(T,e);
 101     }
 102   possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
 103   permInverse(p,q,T->size+1);
 104   /*we put the negative values of weights into a heap, indexed by p
 105     with the minimum value pointed to by p[1]*/
 106   /*p[i] is index (in edgeArray) of edge with i-th position
 107     in the heap, q[j] is the position of edge j in the heap */
 108   while (weights[p[1]] + epsilon < 0)
 109     {
 110       centerEdge = edgeArray[p[1]];
 111       (*count)++;
 112 /*      if (verbose)
 113         {
 114           T->weight = T->weight + weights[p[1]];
 115           printf("New tree weight is %lf.\n",T->weight);
 116         }*/
 117       bNNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
 118       location[p[1]] = NONE;
 119       weights[p[1]] = 0.0;  /*after the bNNI, this edge is in optimal
 120                               configuration*/
 121       popHeap(p,q,weights,possibleSwaps--,1);
 122       /*but we must retest the other edges of T*/
 123       /*CHANGE 2/28/2003 expanding retesting to _all_ edges of T*/
 124       e = depthFirstTraverse(T,NULL);
 125       while (NULL != e)
 126         {
 127           bNNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
 128           e = depthFirstTraverse(T,e);
 129         }
 130     }
 131   free(p);
 132   free(q);
 133   free(location);
 134   free(edgeArray);
 135   free(weights);
 136   assignBMEWeights(T,avgDistArray);
 137 }
 138
 139
 140 /*This function is the meat of the average distance matrix recalculation*/
 141 /*Idea is: we are looking at the subtree rooted at rootEdge.  The subtree
 142 rooted at closer is closer to rootEdge after the NNI, while the subtree
 143 rooted at further is further to rootEdge after the NNI.  direction tells
 144 the direction of the NNI with respect to rootEdge*/
 145 void updateSubTreeAfterNNI(double **A, node *v, edge *rootEdge, node *closer, node *further,
 146                            double dcoeff, int direction)
 147 {
 148   edge *sib;
 149   switch(direction)
 150     {
 151     case UP: /*rootEdge is below the center edge of the NNI*/
 152       /*recursive calls to subtrees, if necessary*/
 153       if (NULL != rootEdge->head->leftEdge)
 154         updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,UP);
 155       if (NULL != rootEdge->head->rightEdge)
 156         updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,UP);
 157       updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
 158       sib = siblingEdge(v->parentEdge);
 159       A[rootEdge->head->index][v->index] =
 160         A[v->index][rootEdge->head->index] =
 161         0.5*A[rootEdge->head->index][sib->head->index] +
 162         0.5*A[rootEdge->head->index][v->parentEdge->tail->index];
 163       break;
 164     case DOWN: /*rootEdge is above the center edge of the NNI*/
 165       sib = siblingEdge(rootEdge);
 166       if (NULL != sib)
 167         updateSubTreeAfterNNI(A, v, sib, closer, further, 0.5*dcoeff, SKEW);
 168       if (NULL != rootEdge->tail->parentEdge)
 169         updateSubTreeAfterNNI(A, v, rootEdge->tail->parentEdge, closer, further, 0.5*dcoeff, DOWN);
 170       updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, DOWN);
 171       A[rootEdge->head->index][v->index] =
 172         A[v->index][rootEdge->head->index] =
 173         0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
 174         0.5*A[rootEdge->head->index][v->rightEdge->head->index];
 175       break;
 176     case SKEW: /*rootEdge is in subtree skew to v*/
 177       if (NULL != rootEdge->head->leftEdge)
 178         updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,SKEW);
 179       if (NULL != rootEdge->head->rightEdge)
 180         updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,SKEW);
 181       updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
 182       A[rootEdge->head->index][v->index] =
 183         A[v->index][rootEdge->head->index] =
 184         0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
 185         0.5*A[rootEdge->head->index][v->rightEdge->head->index];
 186       break;
 187     }
 188 }
 189
 190 /*swapping across edge whose head is v*/
 191 void bNNIupdateAverages(double **A, node *v, edge *par, edge *skew,
 192                         edge *swap, edge *fixed)
 193 {
 194   A[v->index][v->index] = 0.25*(A[fixed->head->index][par->head->index] +
 195                                 A[fixed->head->index][swap->head->index] +
 196                                 A[skew->head->index][par->head->index] +
 197                                 A[skew->head->index][swap->head->index]);
 198   updateSubTreeAfterNNI(A, v, fixed, skew->head, swap->head, 0.25, UP);
 199   updateSubTreeAfterNNI(A, v, par, swap->head, skew->head, 0.25, DOWN);
 200   updateSubTreeAfterNNI(A, v, skew, fixed->head, par->head, 0.25, UP);
 201   updateSubTreeAfterNNI(A, v, swap, par->head, fixed->head, 0.25, SKEW);
 202
 203 }
 204
 205
 206 void bNNItopSwitch(tree *T, edge *e, int direction, double **A)
 207 {
 208   edge *down, *swap, *fixed;
 209   node *u, *v;
 210 /*  if (verbose)
 211     {
 212       printf("Performing branch swap across edge %s ",e->label);
 213       printf("with ");
 214       if (LEFT == direction)
 215         printf("left ");
 216       else printf("right ");
 217       printf("subtree.\n");
 218     }*/
 219   down = siblingEdge(e);
 220   u = e->tail;
 221   v = e->head;
 222   if (LEFT == direction)
 223     {
 224       swap = e->head->leftEdge;
 225       fixed = e->head->rightEdge;
 226       v->leftEdge = down;
 227     }
 228   else
 229     {
 230       swap = e->head->rightEdge;
 231       fixed = e->head->leftEdge;
 232       v->rightEdge = down;
 233     }
 234   swap->tail = u;
 235   down->tail = v;
 236   if(e->tail->leftEdge == e)
 237     u->rightEdge = swap;
 238   else
 239     u->leftEdge = swap;
 240   bNNIupdateAverages(A, v, e->tail->parentEdge, down, swap, fixed);
 241 }
 242
 243 double wf5(double D_AD, double D_BC, double D_AC, double D_BD,
 244            double D_AB, double D_CD)
 245 {
 246   double weight;
 247   weight = 0.25*(D_AC + D_BD + D_AD + D_BC) + 0.5*(D_AB + D_CD);
 248   return(weight);
 249 }
 250
 251 int bNNIEdgeTest(edge *e, tree *T, double **A, double *weight)
 252 {
 253   edge *f;
 254   double D_LR, D_LU, D_LD, D_RD, D_RU, D_DU;
 255   double w1,w2,w0;
 256 /*  if (verbose)
 257     printf("Branch swap: testing edge %s.\n",e->label);*/
 258   if ((leaf(e->tail)) || (leaf(e->head)))
 259     return(NONE);
 260
 261   f = siblingEdge(e);
 262
 263   D_LR = A[e->head->leftEdge->head->index][e->head->rightEdge->head->index];
 264   D_LU = A[e->head->leftEdge->head->index][e->tail->index];
 265   D_LD = A[e->head->leftEdge->head->index][f->head->index];
 266   D_RU = A[e->head->rightEdge->head->index][e->tail->index];
 267   D_RD = A[e->head->rightEdge->head->index][f->head->index];
 268   D_DU = A[e->tail->index][f->head->index];
 269
 270   w0 = wf5(D_RU,D_LD,D_LU,D_RD,D_DU,D_LR); /*weight of current config*/
 271   w1 = wf5(D_RU,D_LD,D_DU,D_LR,D_LU,D_RD); /*weight with L<->D switch*/
 272   w2 = wf5(D_DU,D_LR,D_LU,D_RD,D_RU,D_LD); /*weight with R<->D switch*/
 273   if (w0 <= w1)
 274     {
 275       if (w0 <= w2) /*w0 <= w1,w2*/
 276         {
 277           *weight = 0.0;
 278           return(NONE);
 279         }
 280       else /*w2 < w0 <= w1 */
 281         {
 282           *weight = w2 - w0;
 283 /*        if (verbose)
 284             {
 285               printf("Possible swap across %s. ",e->label);
 286               printf("Weight dropping by %lf.\n",w0 - w2);
 287               printf("New weight would be %lf.\n",T->weight + w2 - w0);
 288             }*/
 289           return(RIGHT);
 290         }
 291     }
 292   else if (w2 <= w1) /*w2 <= w1 < w0*/
 293     {
 294       *weight = w2 - w0;
 295 /*      if (verbose)
 296         {
 297           printf("Possible swap across %s. ",e->label);
 298           printf("Weight dropping by %lf.\n",w0 - w2);
 299           printf("New weight should be %lf.\n",T->weight + w2 - w0);
 300         }*/
 301       return(RIGHT);
 302     }
 303   else /*w1 < w2, w0*/
 304     {
 305       *weight = w1 - w0;
 306 /*      if (verbose)
 307         {
 308           printf("Possible swap across %s. ",e->label);
 309           printf("Weight dropping by %lf.\n",w0 - w1);
 310           printf("New weight should be %lf.\n",T->weight + w1 - w0);
 311         }*/
 312       return(LEFT);
 313     }
 314 }
 315
 316
 317 /*limitedFillTableUp fills all the entries in D associated with
 318   e->head,f->head and those edges g->head above e->head, working
 319   recursively and stopping when trigger is reached*/
 320 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger)
 321 {
 322   edge *g,*h;
 323   g = f->tail->parentEdge;
 324   if (f != trigger)
 325     limitedFillTableUp(e,g,A,trigger);
 326   h = siblingEdge(f);
 327   A[e->head->index][f->head->index] =
 328     A[f->head->index][e->head->index] =
 329     0.5*(A[e->head->index][g->head->index] + A[e->head->index][h->head->index]);
 330 }
 331