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