src/NNI.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 double wf(double lambda, double D_LR, double D_LU, double D_LD,
  10           double D_RU, double D_RD, double D_DU);*/
  11 /*NNI functions for unweighted OLS topological switches*/
  12
  13 /*fillTableUp fills all the entries in D associated with
  14   e->head,f->head and those edges g->head above e->head*/
  15 void fillTableUp(edge *e, edge *f, double **A, double **D, tree *T)
  16 {
  17   edge *g,*h;
  18   if (T->root == f->tail)
  19     {
  20       if (leaf(e->head))
  21         A[e->head->index][f->head->index] =
  22           A[f->head->index][e->head->index] =
  23           D[e->head->index2][f->tail->index2];
  24       else
  25         {
  26           g = e->head->leftEdge;
  27           h = e->head->rightEdge;
  28           A[e->head->index][f->head->index] =
  29             A[f->head->index][e->head->index] =
  30             (g->bottomsize*A[f->head->index][g->head->index]
  31              + h->bottomsize*A[f->head->index][h->head->index])
  32             /e->bottomsize;
  33         }
  34     }
  35   else
  36     {
  37       g = f->tail->parentEdge;
  38       fillTableUp(e,g,A,D,T); /*recursive call*/
  39       h = siblingEdge(f);
  40       A[e->head->index][f->head->index] =
  41         A[f->head->index][e->head->index] =
  42         (g->topsize*A[e->head->index][g->head->index]
  43          + h->bottomsize*A[e->head->index][h->head->index])/f->topsize;
  44     }
  45 }
  46
  47
  48 void makeOLSAveragesTable(tree *T, double **D, double **A);
  49
  50 double **buildAveragesTable(tree *T, double **D)
  51 {
  52   int i,j;
  53   double **A;
  54   A = (double **) malloc(T->size*sizeof(double *));
  55   for(i = 0; i < T->size;i++)
  56     {
  57       A[i] = (double *) malloc(T->size*sizeof(double));
  58       for(j=0;j<T->size;j++)
  59         A[i][j] = 0.0;
  60     }
  61   makeOLSAveragesTable(T,D,A);
  62   return(A);
  63 }
  64
  65 double wf2(double lambda, double D_AD, double D_BC, double D_AC, double D_BD,
  66            double D_AB, double D_CD)
  67 {
  68   double weight;
  69   weight = 0.5*(lambda*(D_AC + D_BD) + (1 - lambda)*(D_AD + D_BC)
  70                 + (D_AB + D_CD));
  71   return(weight);
  72 }
  73
  74 int NNIEdgeTest(edge *e, tree *T, double **A, double *weight)
  75 {
  76   int a,b,c,d;
  77   edge *f;
  78   double *lambda;
  79   double D_LR, D_LU, D_LD, D_RD, D_RU, D_DU;
  80   double w1,w2,w0;
  81
  82   if ((leaf(e->tail)) || (leaf(e->head)))
  83     return(NONE);
  84   lambda = (double *)malloc(3*sizeof(double));
  85   a = e->tail->parentEdge->topsize;
  86   f = siblingEdge(e);
  87   b = f->bottomsize;
  88   c = e->head->leftEdge->bottomsize;
  89   d = e->head->rightEdge->bottomsize;
  90
  91   lambda[0] = ((double) b*c + a*d)/((a + b)*(c+d));
  92   lambda[1] = ((double) b*c + a*d)/((a + c)*(b+d));
  93   lambda[2] = ((double) c*d + a*b)/((a + d)*(b+c));
  94
  95   D_LR = A[e->head->leftEdge->head->index][e->head->rightEdge->head->index];
  96   D_LU = A[e->head->leftEdge->head->index][e->tail->index];
  97   D_LD = A[e->head->leftEdge->head->index][f->head->index];
  98   D_RU = A[e->head->rightEdge->head->index][e->tail->index];
  99   D_RD = A[e->head->rightEdge->head->index][f->head->index];
 100   D_DU = A[e->tail->index][f->head->index];
 101
 102   w0 = wf2(lambda[0],D_RU,D_LD,D_LU,D_RD,D_DU,D_LR);
 103   w1 = wf2(lambda[1],D_RU,D_LD,D_DU,D_LR,D_LU,D_RD);
 104   w2 = wf2(lambda[2],D_DU,D_LR,D_LU,D_RD,D_RU,D_LD);
 105   free(lambda);
 106   if (w0 <= w1)
 107     {
 108       if (w0 <= w2) /*w0 <= w1,w2*/
 109         {
 110           *weight = 0.0;
 111           return(NONE);
 112         }
 113       else /*w2 < w0 <= w1 */
 114         {
 115           *weight = w2 - w0;
 116 /*        if (verbose)
 117             {
 118               printf("Swap across %s. ",e->label);
 119               printf("Weight dropping by %lf.\n",w0 - w2);
 120               printf("New weight should be %lf.\n",T->weight + w2 - w0);
 121             }*/
 122           return(RIGHT);
 123         }
 124     }
 125   else if (w2 <= w1) /*w2 <= w1 < w0*/
 126     {
 127       *weight = w2 - w0;
 128 /*      if (verbose)
 129         {
 130           printf("Swap across %s. ",e->label);
 131           printf("Weight dropping by %lf.\n",w0 - w2);
 132           printf("New weight should be %lf.\n",T->weight + w2 - w0);
 133         }*/
 134       return(RIGHT);
 135     }
 136   else /*w1 < w2, w0*/
 137     {
 138       *weight = w1 - w0;
 139 /*      if (verbose)
 140         {
 141           printf("Swap across %s. ",e->label);
 142           printf("Weight dropping by %lf.\n",w0 - w1);
 143           printf("New weight should be %lf.\n",T->weight + w1 - w0);
 144         }*/
 145       return(LEFT);
 146     }
 147 }
 148
 149 int *initPerm(int size);
 150
 151 void NNIupdateAverages(double **A, edge *e, edge *par, edge *skew,
 152                        edge *swap, edge *fixed, tree *T)
 153 {
 154   node *v;
 155   edge *elooper;
 156   v = e->head;
 157   /*first, v*/
 158   A[e->head->index][e->head->index] =
 159     (swap->bottomsize*
 160      ((skew->bottomsize*A[skew->head->index][swap->head->index]
 161        + fixed->bottomsize*A[fixed->head->index][swap->head->index])
 162       / e->bottomsize) +
 163      par->topsize*
 164      ((skew->bottomsize*A[skew->head->index][par->head->index]
 165        + fixed->bottomsize*A[fixed->head->index][par->head->index])
 166       / e->bottomsize)
 167      ) / e->topsize;
 168
 169   elooper = findBottomLeft(e); /*next, we loop over all the edges
 170                                  which are below e*/
 171   while (e != elooper)
 172     {
 173       A[e->head->index][elooper->head->index] =
 174         A[elooper->head->index][v->index]
 175         = (swap->bottomsize*A[elooper->head->index][swap->head->index] +
 176            par->topsize*A[elooper->head->index][par->head->index])
 177         / e->topsize;
 178       elooper = depthFirstTraverse(T,elooper);
 179     }
 180   elooper = findBottomLeft(swap); /*next we loop over all the edges below and
 181                                     including swap*/
 182   while (swap != elooper)
 183   {
 184     A[e->head->index][elooper->head->index] =
 185       A[elooper->head->index][e->head->index]
 186       = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
 187          fixed->bottomsize*A[elooper->head->index][fixed->head->index])
 188       / e->bottomsize;
 189     elooper = depthFirstTraverse(T,elooper);
 190   }
 191   /*now elooper = skew */
 192   A[e->head->index][elooper->head->index] =
 193     A[elooper->head->index][e->head->index]
 194     = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
 195        fixed->bottomsize* A[elooper->head->index][fixed->head->index])
 196     / e->bottomsize;
 197
 198   /*finally, we loop over all the edges in the tree
 199     on the far side of parEdge*/
 200   elooper = T->root->leftEdge;
 201   while ((elooper != swap) && (elooper != e)) /*start a top-first traversal*/
 202     {
 203       A[e->head->index][elooper->head->index] =
 204         A[elooper->head->index][e->head->index]
 205         = (skew->bottomsize * A[elooper->head->index][skew->head->index]
 206            + fixed->bottomsize* A[elooper->head->index][fixed->head->index])
 207         / e->bottomsize;
 208       elooper = topFirstTraverse(T,elooper);
 209     }
 210
 211   /*At this point, elooper = par.
 212     We finish the top-first traversal, excluding the subtree below par*/
 213   elooper = moveUpRight(par);
 214   while (NULL != elooper)
 215     {
 216       A[e->head->index][elooper->head->index]
 217         = A[elooper->head->index][e->head->index]
 218         = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
 219            fixed->bottomsize* A[elooper->head->index][fixed->head->index])
 220         / e->bottomsize;
 221       elooper = topFirstTraverse(T,elooper);
 222     }
 223 }
 224
 225
 226 void NNItopSwitch(tree *T, edge *e, int direction, double **A)
 227 {
 228   edge *par, *fixed;
 229   edge *skew, *swap;
 230
 231 /*  if (verbose)
 232     printf("Branch swap across edge %s.\n",e->label);*/
 233
 234   if (LEFT == direction)
 235     swap = e->head->leftEdge;
 236   else
 237     swap = e->head->rightEdge;
 238   skew = siblingEdge(e);
 239   fixed = siblingEdge(swap);
 240   par = e->tail->parentEdge;
 241
 242 /*  if (verbose)
 243     {
 244       printf("Branch swap: switching edges %s and %s.\n",skew->label,swap->label);
 245     }*/
 246   /*perform topological switch by changing f from (u,b) to (v,b)
 247     and g from (v,c) to (u,c), necessitatates also changing parent fields*/
 248
 249   swap->tail = e->tail;
 250   skew->tail = e->head;
 251
 252   if (LEFT == direction)
 253     e->head->leftEdge = skew;
 254   else
 255     e->head->rightEdge = skew;
 256   if (skew == e->tail->rightEdge)
 257     e->tail->rightEdge = swap;
 258   else
 259     e->tail->leftEdge = swap;
 260
 261   /*both topsize and bottomsize change for e, but nowhere else*/
 262
 263   e->topsize = par->topsize + swap->bottomsize;
 264   e->bottomsize = fixed->bottomsize + skew->bottomsize;
 265   NNIupdateAverages(A, e, par, skew, swap, fixed,T);
 266
 267 } /*end NNItopSwitch*/
 268
 269 void reHeapElement(int *p, int *q, double *v, int length, int i);
 270 void pushHeap(int *p, int *q, double *v, int length, int i);
 271 void popHeap(int *p, int *q, double *v, int length, int i);
 272
 273
 274 void NNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray,
 275                 double *weights, int *location, int *possibleSwaps)
 276 {
 277   int tloc;
 278   tloc = location[e->head->index+1];
 279   location[e->head->index+1] =
 280     NNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
 281   if (NONE == location[e->head->index+1])
 282     {
 283       if (NONE != tloc)
 284         popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
 285     }
 286   else
 287     {
 288       if (NONE == tloc)
 289         pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
 290       else
 291         reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
 292     }
 293 }
 294
 295 void permInverse(int *p, int *q, int length);
 296
 297 int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
 298
 299
 300 //void NNI(tree *T, double **avgDistArray, int *count)
 301 void NNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
 302 {
 303   edge *e, *centerEdge;
 304   edge **edgeArray;
 305   int *location;
 306   int *p,*q;
 307   int i,j;
 308   int possibleSwaps;
 309   double *weights;
 310   p = initPerm(T->size+1);
 311   q = initPerm(T->size+1);
 312   edgeArray = (edge **) malloc((T->size+1)*sizeof(double));
 313   weights = (double *) malloc((T->size+1)*sizeof(double));
 314   location = (int *) malloc((T->size+1)*sizeof(int));
 315
 316   double epsilon = 0.0;
 317   for (i=0; i<numSpecies; i++)
 318     for (j=0; j<numSpecies; j++)
 319       epsilon += D[i][j];
 320   epsilon = (epsilon / (numSpecies * numSpecies)) * EPSILON;
 321
 322   for (i=0;i<T->size+1;i++)
 323     {
 324       weights[i] = 0.0;
 325       location[i] = NONE;
 326     }
 327   e = findBottomLeft(T->root->leftEdge);
 328   /* *count = 0; */
 329   while (NULL != e)
 330     {
 331       edgeArray[e->head->index+1] = e;
 332       location[e->head->index+1] =
 333         NNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
 334       e = depthFirstTraverse(T,e);
 335     }
 336   possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
 337   permInverse(p,q,T->size+1);
 338   /*we put the negative values of weights into a heap, indexed by p
 339     with the minimum value pointed to by p[1]*/
 340   /*p[i] is index (in edgeArray) of edge with i-th position
 341     in the heap, q[j] is the position of edge j in the heap */
 342   while (weights[p[1]] + epsilon < 0)
 343     {
 344       centerEdge = edgeArray[p[1]];
 345       (*count)++;
 346       T->weight = T->weight + weights[p[1]];
 347       NNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
 348       location[p[1]] = NONE;
 349       weights[p[1]] = 0.0;  /*after the NNI, this edge is in optimal
 350                               configuration*/
 351       popHeap(p,q,weights,possibleSwaps--,1);
 352       /*but we must retest the other four edges*/
 353       e = centerEdge->head->leftEdge;
 354       NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
 355       e = centerEdge->head->rightEdge;
 356       NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
 357       e = siblingEdge(centerEdge);
 358       NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
 359       e = centerEdge->tail->parentEdge;
 360       NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
 361     }
 362   free(p);
 363   free(q);
 364   free(location);
 365   free(edgeArray);
 366 }
 367 /*
 368 void NNIwithoutMatrix(tree *T, double **D, int *count)
 369 {
 370   double **avgDistArray;
 371   avgDistArray = buildAveragesTable(T,D);
 372   NNI(T,avgDistArray,count);
 373 }
 374
 375 void NNIWithPartialMatrix(tree *T,double **D,double **A,int *count)
 376 {
 377   makeOLSAveragesTable(T,D,A);
 378   NNI(T,A,count);
 379 }
 380 */