src/SPR.c

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