src/bNNI.c

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