src/SPR.c

   1 /* SPR.c    2009-01-12 */
   2
   3 /* Copyright 2009 Richard Desper */
   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 /*functions from bNNI.c*/
  11 void makeBMEAveragesTable(tree *T, double **D, double **A);
  12 void assignBMEWeights(tree *T, double **A);
  13
  14 /*from me.c*/
  15 edge *siblingEdge(edge *e);
  16 void weighTree(tree *T);
  17 void freeMatrix(double **D, int size);
  18 edge *depthFirstTraverse(tree *T, edge *e);
  19 double **initDoubleMatrix(int d);
  20
  21 /*from below*/
  22 node *indexedNode(tree *T, int i);
  23 edge *indexedEdge(tree *T, int i);
  24 void assignSPRWeights(node *v, double **A, double ***swapWeights);
  25 void SPRTopShift(tree *T, node *vmove, edge *esplit, int UpOrDown);
  26 void assignDownWeightsUp(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights);
  27 void assignDownWeightsSkew(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights);
  28 void assignDownWeightsDown(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights);
  29 void assignUpWeights(edge *etest, node *vtest, node *va, edge *back, node *cprve, double oldD_AB, double coeff, double **A, double ***swapWeights);
  30
  31 void zero3DMatrix(double ***X, int h, int l, int w)
  32 {
  33         int i,j,k;
  34         for(i=0;i<h;i++)
  35                 for(j=0;j<l;j++)
  36                         for(k=0;k<w;k++)
  37                                 X[i][j][k] = 0.0;
  38 }
  39
  40
  41 void findTableMin(int *imin, int *jmin, int *kmin, int n, double ***X, double *min)
  42 {
  43   int i,j,k;
  44   for(i=0;i<2;i++)
  45     for(j=0;j<n;j++)
  46       for(k=0;k<n;k++)
  47         {
  48           if (X[i][j][k] < *min)
  49             {
  50               *min = X[i][j][k];
  51               *imin = i;
  52               *jmin = j;
  53               *kmin = k;
  54             }
  55         }
  56 }
  57
  58
  59 void SPR(tree *T, double **D, double **A, int *count)
  60 {
  61   int i,j,k;
  62   node *v;
  63   /*FILE *treefile;*/
  64   edge *e,*f;
  65   /* char filename[MAX_LABEL_LENGTH];*/
  66   double ***swapWeights;
  67   double swapValue = 0.0;
  68   swapWeights = (double ***)malloc(2*sizeof(double **));
  69   makeBMEAveragesTable(T,D,A);
  70   assignBMEWeights(T,A);
  71   weighTree(T);
  72   /*if (verbose)
  73     printf("Before SPRs: tree length is %lf.\n",T->weight);*/
  74   for(i=0;i<2;i++)
  75     swapWeights[i] = initDoubleMatrix(T->size);
  76   do
  77     {
  78       swapValue=0.0;
  79       zero3DMatrix(swapWeights,2,T->size,T->size);
  80       i = j = k = 0;
  81       for(e=depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
  82         assignSPRWeights(e->head,A,swapWeights);
  83       findTableMin(&i,&j,&k,T->size,swapWeights,&swapValue);
  84       swapValue = swapWeights[i][j][k];
  85       if (swapValue < -EPSILON)
  86         {
  87 //        if (verbose)
  88 //          printf("New tree weight should be %lf.\n",T->weight + 0.25*swapValue);
  89           v = indexedNode(T,j);
  90           f = indexedEdge(T,k);
  91 //        if (verbose)
  92 //          printf("Swapping tree below %s to split edge %s with head %s and tail %s\n",
  93 //                         v->parentEdge->label,f->label,f->head->label,f->tail->label);
  94           SPRTopShift(T,v,f,2-i);
  95           makeBMEAveragesTable(T,D,A);
  96           assignBMEWeights(T,A);
  97           weighTree(T);
  98           (*count)++;
  99           /*sprintf(filename,"tree%d.new",*count);*/
 100 //        if (verbose)
 101 //          printf("After %d SPRs, tree weight is %lf.\n\n",*count,T->weight);
 102           /*treefile = fopen(filename,"w");
 103           NewickPrintTree(T,treefile);
 104           fclose(treefile);*/
 105           }
 106     } while (swapValue < -EPSILON);
 107   for(i=0;i<2;i++)
 108     freeMatrix(swapWeights[i],T->size);
 109   free(swapWeights);
 110   /*if (verbose)
 111     readOffTree(T);*/
 112 }
 113
 114 /*assigns values to array swapWeights*/
 115 /*swapWeights[0][j][k] will be the value of removing the tree below the edge whose head node has index j
 116 and reattaching it to split the edge whose head has the index k*/
 117 /*swapWeights[1][j][k] will be the value of removing the tree above the edge whose head node has index j
 118 and reattaching it to split the edge whose head has the index k*/
 119 void assignSPRWeights(node *vtest, double **A, double ***swapWeights)
 120 {
 121   edge *etest, *left, *right, *sib, *par;
 122   etest = vtest->parentEdge;
 123   left = vtest->leftEdge;
 124   right = vtest->rightEdge;
 125   par = etest->tail->parentEdge;
 126   sib = siblingEdge(etest);
 127   if (NULL != par)
 128     assignDownWeightsUp(par,vtest,sib->head,NULL,NULL,0.0,1.0,A,swapWeights);
 129   if (NULL != sib)
 130     assignDownWeightsSkew(sib,vtest,sib->tail,NULL,NULL,0.0,1.0,A,swapWeights);
 131   /*assigns values for moving subtree rooted at vtest, starting with edge
 132     parental to tail of edge parental to vtest*/
 133   if (NULL != left)
 134     {
 135       assignUpWeights(left,vtest,right->head,NULL,NULL,0.0,1.0,A,swapWeights);
 136       assignUpWeights(right,vtest,left->head,NULL,NULL,0.0,1.0,A,swapWeights);
 137     }
 138 }
 139
 140
 141 /*recall NNI formula: change in tree length from AB|CD split to AC|BD split is
 142 proportional to D_AC + D_BD - D_AB - D_CD*/
 143 /*in our case B is the tree being moved (below vtest), A is the tree backwards below back, but
 144   with the vtest subtree removed, C is the sibling tree of back and D is the tree above etest*/
 145 /*use va to denote the root of the sibling tree to B in the original tree*/
 146 /*please excuse the multiple uses of the same letters: A,D, etc.*/
 147 void assignDownWeightsUp(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights)
 148 {
 149   edge *par, *sib, *skew;
 150   double D_AC, D_BD, D_AB, D_CD;
 151   par = etest->tail->parentEdge;
 152   skew = siblingEdge(etest);
 153   if (NULL == back) /*first recursive call*/
 154     {
 155       if (NULL == par)
 156         return;
 157       else /*start the process of assigning weights recursively*/
 158         {
 159           assignDownWeightsUp(par,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
 160           assignDownWeightsSkew(skew,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
 161         }
 162     }
 163   else /*second or later recursive call*/
 164     {
 165       sib = siblingEdge(back);
 166       D_BD = A[vtest->index][etest->head->index]; /*straightforward*/
 167       D_CD = A[sib->head->index][etest->head->index]; /*this one too*/
 168       D_AC = A[sib->head->index][back->head->index] + coeff*(A[sib->head->index][va->index]
 169                                                              - A[sib->head->index][vtest->index]);
 170       D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
 171       swapWeights[0][vtest->index][etest->head->index] = swapWeights[0][vtest->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
 172       if (NULL != par)
 173         {
 174           assignDownWeightsUp(par,vtest,va,etest,sib->head,D_AB,0.5*coeff,A,swapWeights);
 175           assignDownWeightsSkew(skew,vtest,va,etest,sib->head,D_AB,0.5*coeff,A,swapWeights);
 176         }
 177     }
 178 }
 179
 180 void assignDownWeightsSkew(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights)
 181 {
 182   /*same general idea as assignDownWeights, except needing to keep track of things a bit differently*/
 183   edge *par, *left, *right;
 184   /*par here = sib before
 185     left, right here = par, skew before*/
 186   double D_AB, D_CD, D_AC, D_BD;
 187   /*B is subtree being moved - below vtest
 188     A is subtree remaining fixed - below va, unioned with all trees already passed by B*/
 189   /*C is subtree being passed by B, in this case above par
 190     D is subtree below etest, fixed on other side*/
 191   par = etest->tail->parentEdge;
 192   left = etest->head->leftEdge;
 193   right = etest->head->rightEdge;
 194   if (NULL == back)
 195     {
 196       if (NULL == left)
 197         return;
 198       else
 199         {
 200           assignDownWeightsDown(left,vtest,va,etest,etest->tail,A[vtest->index][etest->tail->index],0.5,A,swapWeights);
 201           assignDownWeightsDown(right,vtest,va,etest,etest->tail,A[vtest->index][etest->tail->index],0.5,A,swapWeights);
 202         }
 203     }
 204   else
 205     {
 206       D_BD = A[vtest->index][etest->head->index];
 207       D_CD = A[par->head->index][etest->head->index];
 208       D_AC = A[back->head->index][par->head->index] + coeff*(A[va->index][par->head->index] - A[vtest->index][par->head->index]);
 209       D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
 210       swapWeights[0][vtest->index][etest->head->index] = swapWeights[0][vtest->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
 211       if (NULL != left)
 212         {
 213           assignDownWeightsDown(left,vtest, va, etest, etest->tail, D_AB, 0.5*coeff, A, swapWeights);
 214           assignDownWeightsDown(right,vtest, va, etest, etest->tail, D_AB, 0.5*coeff, A, swapWeights);
 215         }
 216     }
 217 }
 218
 219 void assignDownWeightsDown(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A, double ***swapWeights)
 220 {
 221   /*again the same general idea*/
 222   edge *sib, *left, *right;
 223   /*sib here = par in assignDownWeightsSkew
 224     rest is the same as assignDownWeightsSkew*/
 225   double D_AB, D_CD, D_AC, D_BD;
 226   /*B is below vtest, A is below va unioned with all trees already passed by B*/
 227   /*C is subtree being passed - below sib*/
 228   /*D is tree below etest*/
 229   sib = siblingEdge(etest);
 230   left = etest->head->leftEdge;
 231   right = etest->head->rightEdge;
 232   D_BD = A[vtest->index][etest->head->index];
 233   D_CD = A[sib->head->index][etest->head->index];
 234   D_AC = A[sib->head->index][back->head->index] + coeff*(A[sib->head->index][va->index] - A[sib->head->index][vtest->index]);
 235   D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
 236   swapWeights[0][vtest->index][etest->head->index] = swapWeights[0][vtest->index][back->head->index] + ( D_AC + D_BD - D_AB - D_CD);
 237   if (NULL != left)
 238     {
 239       assignDownWeightsDown(left,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
 240       assignDownWeightsDown(right,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
 241     }
 242 }
 243
 244 void assignUpWeights(edge *etest, node *vtest, node *va, edge *back, node *cprev, double oldD_AB, double coeff, double **A,
 245                      double ***swapWeights)
 246 {
 247         /*SPR performed on tree above vtest...*/
 248         /*same idea as above, with appropriate selections of edges and nodes*/
 249   edge *sib, *left, *right;
 250   /*B is above vtest, A is other tree below vtest unioned with trees in path to vtest*/
 251   /*sib is tree C being passed by B*/
 252   /*D is tree below etest*/
 253   double D_AB, D_CD, D_AC, D_BD;
 254   double thisWeight;
 255   sib = siblingEdge(etest);
 256   left = etest->head->leftEdge;
 257   right = etest->head->rightEdge;
 258   if (NULL == back) /*first recursive call*/
 259     {
 260       if (NULL == left)
 261         return;
 262       else /*start the process of assigning weights recursively*/
 263         {
 264           assignUpWeights(left,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
 265           assignUpWeights(right,vtest,va,etest,va,A[va->index][vtest->index],0.5,A,swapWeights);
 266         }
 267     }
 268   else /*second or later recursive call*/
 269     {
 270       D_BD = A[vtest->index][etest->head->index];
 271       D_CD = A[sib->head->index][etest->head->index];
 272       D_AC = A[back->head->index][sib->head->index] + coeff*(A[va->index][sib->head->index] - A[vtest->index][sib->head->index]);
 273       D_AB = 0.5*(oldD_AB + A[vtest->index][cprev->index]);
 274       thisWeight = swapWeights[1][vtest->index][etest->head->index] = swapWeights[1][vtest->index][back->head->index] + (D_AC + D_BD - D_AB - D_CD);
 275       if (NULL != left)
 276         {
 277           assignUpWeights(left,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
 278           assignUpWeights(right,vtest, va, etest, sib->head, D_AB, 0.5*coeff, A, swapWeights);
 279         }
 280     }
 281 }
 282
 283 void pruneSubtree(edge *p, edge *u, edge *d)
 284 /*starting with edge u above edges p, d*/
 285 /*removes p, d from tree, u connects to d->head to compensate*/
 286 {
 287   p->tail->parentEdge = NULL;/*remove p subtree*/
 288   u->head = d->head;
 289   d->head->parentEdge = u;      /*u connected to d->head*/
 290   d->head = NULL; /*d removed from tree*/
 291 }
 292
 293 void SPRsplitEdge(edge *e, edge *p, edge *d)
 294 /*splits edge e to make it parental to p,d.  d is parental to what
 295   previously was below e*/
 296 {
 297   d->head = e->head;
 298   e->head = p->tail;
 299   p->tail->parentEdge = e;
 300   d->head->parentEdge = d;
 301 }
 302
 303
 304 /*topological shift function*/
 305 /*removes subtree rooted at v and re-inserts to spilt e*/
 306 void SPRDownShift(tree *T, node *v, edge *e)
 307 {
 308   edge *vup, *vdown, *vpar;
 309   vpar = v->parentEdge;
 310   vdown = siblingEdge(vpar);
 311   vup = vpar->tail->parentEdge;
 312   /*topological shift*/
 313   pruneSubtree(vpar,vup,vdown);
 314   /*removes v subtree and vdown, extends vup*/
 315   SPRsplitEdge(e,vpar,vdown); /*splits e to make e sibling edge to vpar,
 316                                 both below vup*/
 317 }
 318
 319 void SPRUpShift(tree *T, node *vmove, edge *esplit)
 320 /*an inelegant iterative version*/
 321 {
 322   edge *f;
 323   edge **EPath;
 324   edge **sib;
 325   node **v;
 326   int i;
 327   int pathLength;
 328   for(f=esplit->tail->parentEdge,pathLength=1;f->tail != vmove;f=f->tail->parentEdge)
 329     pathLength++;
 330   /*count number of edges to vmove*/
 331   /*note that pathLength > 0 will hold*/
 332
 333   EPath = (edge **)malloc(pathLength*sizeof(edge *));
 334   v = (node **)malloc(pathLength*sizeof(edge *));
 335   sib = (edge **)malloc((pathLength+1)*sizeof(edge *));
 336   /*there are pathLength + 1 side trees, one at the head and tail of each edge in the path*/
 337
 338   sib[pathLength] = siblingEdge(esplit);
 339   i = pathLength;
 340   f = esplit->tail->parentEdge;
 341   while (i > 0)
 342     {
 343       i--;
 344       EPath[i] = f;
 345       sib[i] = siblingEdge(f);
 346       v[i] = f->head;
 347       f = f->tail->parentEdge;
 348     }
 349   /*indexed so head of Epath is v[i], tail is v[i-1] and sibling edge is sib[i]*/
 350   /*need to assign head, tail of each edge in path
 351     as well as have proper values for the left and right fields*/
 352
 353   if (esplit == esplit->tail->leftEdge)
 354     {
 355       vmove->leftEdge = esplit;
 356       vmove->rightEdge = EPath[pathLength-1];
 357     }
 358   else
 359     {
 360       vmove->rightEdge = esplit;
 361       vmove->leftEdge = EPath[pathLength-1];
 362     }
 363   esplit->tail = vmove;
 364   /*espilt->head remains unchanged*/
 365   /*vmove has the proper fields for left, right, and parentEdge*/
 366
 367   for(i=0;i<(pathLength-1);i++)
 368     EPath[i]->tail = v[i+1];
 369
 370   /*this bit flips the orientation along the path
 371     the tail of Epath[i] is now v[i+1] instead of v[i-1]*/
 372
 373   EPath[pathLength-1]->tail = vmove;
 374
 375   for(i=1;i<pathLength;i++)
 376     {
 377       if (sib[i+1] == v[i]->leftEdge)
 378         v[i]->rightEdge = EPath[i-1];
 379       else
 380         v[i]->leftEdge = EPath[i-1];
 381     }
 382   if (sib[1] == v[0]->leftEdge)
 383     v[0]->rightEdge = sib[0];
 384   else
 385     v[0]->leftEdge = sib[0];
 386   sib[0]->tail = v[0];
 387   free(EPath);
 388   free(v);
 389   free(sib);
 390 }
 391
 392
 393 void SPRTopShift(tree *T, node *vmove, edge *esplit, int UpOrDown)
 394 {
 395   if (DOWN == UpOrDown)
 396     SPRDownShift(T,vmove,esplit);
 397   else
 398     SPRUpShift(T,vmove,esplit);
 399 }
 400
 401 node *indexedNode(tree *T, int i)
 402 {
 403   edge *e;
 404   for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
 405     if (i == e->head->index)
 406       return(e->head);
 407   if (i == T->root->index)
 408     return(T->root);
 409   return(NULL);
 410 }
 411
 412 edge *indexedEdge(tree *T, int i)
 413 {
 414   edge *e;
 415   for(e = depthFirstTraverse(T,NULL);NULL!=e;e=depthFirstTraverse(T,e))
 416     if (i == e->head->index)
 417       return(e);
 418   return(NULL);
 419 }