src/bNNI.c

   1 /* bNNI.c    2007-09-05 */
   2
   3 /* Copyright 2007 Vincent Lefort */
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 #include "me.h"
   9
  10 /*boolean leaf(node *v);
  11 edge *siblingEdge(edge *e);
  12 edge *depthFirstTraverse(tree *T, edge *e);
  13 edge *findBottomLeft(edge *e);
  14 edge *topFirstTraverse(tree *T, edge *e);
  15 edge *moveUpRight(edge *e);*/
  16
  17 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger);
  18 void assignBMEWeights(tree *T, double **A);
  19 //void updateAveragesMatrix(tree *T, double **A, node *v,int direction);
  20 void bNNItopSwitch(tree *T, edge *e, int direction, double **A);
  21 int bNNIEdgeTest(edge *e, tree *T, double **A, double *weight);
  22 void updatePair(double **A, edge *nearEdge, edge *farEdge, node *closer, node *further, double dcoeff, int direction);
  23
  24 int *initPerm(int size);
  25
  26 void reHeapElement(int *p, int *q, double *v, int length, int i);
  27 void pushHeap(int *p, int *q, double *v, int length, int i);
  28 void popHeap(int *p, int *q, double *v, int length, int i);
  29
  30
  31 void bNNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray,
  32                 double *weights, int *location, int *possibleSwaps)
  33 {
  34   int tloc;
  35   tloc = location[e->head->index+1];
  36   location[e->head->index+1] =
  37     bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
  38   if (NONE == location[e->head->index+1])
  39     {
  40       if (NONE != tloc)
  41         popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
  42     }
  43   else
  44     {
  45       if (NONE == tloc)
  46         pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
  47       else
  48         reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
  49     }
  50 }
  51
  52 int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
  53
  54 void permInverse(int *p, int *q, int length);
  55
  56 void weighTree(tree *T)
  57 {
  58   edge *e;
  59   T->weight = 0;
  60   for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
  61     T->weight += e->distance;
  62 }
  63
  64 //void bNNI(tree *T, double **avgDistArray, int *count)
  65 void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
  66 {
  67   edge *e, *centerEdge;
  68   edge **edgeArray;
  69   int *p, *location, *q;
  70   int i,j;
  71   int possibleSwaps;
  72   double *weights;
  73   p = initPerm(T->size+1);
  74   q = initPerm(T->size+1);
  75   edgeArray = (edge **) malloc((T->size+1)*sizeof(double));
  76   weights = (double *) malloc((T->size+1)*sizeof(double));
  77   location = (int *) malloc((T->size+1)*sizeof(int));
  78
  79   double epsilon = 0.0;
  80   for (i=0; i<numSpecies; i++)
  81     for (j=0; j<numSpecies; j++)
  82       epsilon += D[i][j];
  83   epsilon = (epsilon / (numSpecies * numSpecies)) * EPSILON;
  84
  85   for (i=0;i<T->size+1;i++)
  86     {
  87       weights[i] = 0.0;
  88       location[i] = NONE;
  89     }
  90 /*  if (verbose)
  91     {
  92       assignBMEWeights(T,avgDistArray);
  93       weighTree(T);
  94     }*/
  95   e = findBottomLeft(T->root->leftEdge);
  96   while (NULL != e)
  97     {
  98       edgeArray[e->head->index+1] = e;
  99       location[e->head->index+1] =
 100         bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
 101       e = depthFirstTraverse(T,e);
 102     }
 103   possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
 104   permInverse(p,q,T->size+1);
 105   /*we put the negative values of weights into a heap, indexed by p
 106     with the minimum value pointed to by p[1]*/
 107   /*p[i] is index (in edgeArray) of edge with i-th position
 108     in the heap, q[j] is the position of edge j in the heap */
 109   while (weights[p[1]] + epsilon < 0)
 110     {
 111       centerEdge = edgeArray[p[1]];
 112       (*count)++;
 113 /*      if (verbose)
 114         {
 115           T->weight = T->weight + weights[p[1]];
 116           printf("New tree weight is %lf.\n",T->weight);
 117         }*/
 118       bNNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
 119       location[p[1]] = NONE;
 120       weights[p[1]] = 0.0;  /*after the bNNI, this edge is in optimal
 121                               configuration*/
 122       popHeap(p,q,weights,possibleSwaps--,1);
 123       /*but we must retest the other edges of T*/
 124       /*CHANGE 2/28/2003 expanding retesting to _all_ edges of T*/
 125       e = depthFirstTraverse(T,NULL);
 126       while (NULL != e)
 127         {
 128           bNNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
 129           e = depthFirstTraverse(T,e);
 130         }
 131     }
 132   free(p);
 133   free(q);
 134   free(location);
 135   free(edgeArray);
 136   free(weights);
 137   assignBMEWeights(T,avgDistArray);
 138 }
 139
 140
 141 /*This function is the meat of the average distance matrix recalculation*/
 142 /*Idea is: we are looking at the subtree rooted at rootEdge.  The subtree
 143 rooted at closer is closer to rootEdge after the NNI, while the subtree
 144 rooted at further is further to rootEdge after the NNI.  direction tells
 145 the direction of the NNI with respect to rootEdge*/
 146 void updateSubTreeAfterNNI(double **A, node *v, edge *rootEdge, node *closer, node *further,
 147                            double dcoeff, int direction)
 148 {
 149   edge *sib;
 150   switch(direction)
 151     {
 152     case UP: /*rootEdge is below the center edge of the NNI*/
 153       /*recursive calls to subtrees, if necessary*/
 154       if (NULL != rootEdge->head->leftEdge)
 155         updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,UP);
 156       if (NULL != rootEdge->head->rightEdge)
 157         updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,UP);
 158       updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
 159       sib = siblingEdge(v->parentEdge);
 160       A[rootEdge->head->index][v->index] =
 161         A[v->index][rootEdge->head->index] =
 162         0.5*A[rootEdge->head->index][sib->head->index] +
 163         0.5*A[rootEdge->head->index][v->parentEdge->tail->index];
 164       break;
 165     case DOWN: /*rootEdge is above the center edge of the NNI*/
 166       sib = siblingEdge(rootEdge);
 167       if (NULL != sib)
 168         updateSubTreeAfterNNI(A, v, sib, closer, further, 0.5*dcoeff, SKEW);
 169       if (NULL != rootEdge->tail->parentEdge)
 170         updateSubTreeAfterNNI(A, v, rootEdge->tail->parentEdge, closer, further, 0.5*dcoeff, DOWN);
 171       updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, DOWN);
 172       A[rootEdge->head->index][v->index] =
 173         A[v->index][rootEdge->head->index] =
 174         0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
 175         0.5*A[rootEdge->head->index][v->rightEdge->head->index];
 176       break;
 177     case SKEW: /*rootEdge is in subtree skew to v*/
 178       if (NULL != rootEdge->head->leftEdge)
 179         updateSubTreeAfterNNI(A, v, rootEdge->head->leftEdge, closer, further, 0.5*dcoeff,SKEW);
 180       if (NULL != rootEdge->head->rightEdge)
 181         updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,SKEW);
 182       updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP);
 183       A[rootEdge->head->index][v->index] =
 184         A[v->index][rootEdge->head->index] =
 185         0.5*A[rootEdge->head->index][v->leftEdge->head->index] +
 186         0.5*A[rootEdge->head->index][v->rightEdge->head->index];
 187       break;
 188     }
 189 }
 190
 191 /*swapping across edge whose head is v*/
 192 void bNNIupdateAverages(double **A, node *v, edge *par, edge *skew,
 193                         edge *swap, edge *fixed)
 194 {
 195   A[v->index][v->index] = 0.25*(A[fixed->head->index][par->head->index] +
 196                                 A[fixed->head->index][swap->head->index] +
 197                                 A[skew->head->index][par->head->index] +
 198                                 A[skew->head->index][swap->head->index]);
 199   updateSubTreeAfterNNI(A, v, fixed, skew->head, swap->head, 0.25, UP);
 200   updateSubTreeAfterNNI(A, v, par, swap->head, skew->head, 0.25, DOWN);
 201   updateSubTreeAfterNNI(A, v, skew, fixed->head, par->head, 0.25, UP);
 202   updateSubTreeAfterNNI(A, v, swap, par->head, fixed->head, 0.25, SKEW);
 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 /*limitedFillTableUp fills all the entries in D associated with
 317   e->head,f->head and those edges g->head above e->head, working
 318   recursively and stopping when trigger is reached*/
 319 void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger)
 320 {
 321   edge *g,*h;
 322   g = f->tail->parentEdge;
 323   if (f != trigger)
 324     limitedFillTableUp(e,g,A,trigger);
 325   h = siblingEdge(f);
 326   A[e->head->index][f->head->index] =
 327     A[f->head->index][e->head->index] =
 328     0.5*(A[e->head->index][g->head->index] + A[e->head->index][h->head->index]);
 329 }