src/NNI.c

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