src/me_balanced.c

   1 #include "me.h"
   2
   3 void BalWFext(edge *e, double **A) /*works except when e is the one edge
   4                                   inserted to new vertex v by firstInsert*/
   5 {
   6   edge *f, *g;
   7   if ((leaf(e->head)) && (leaf(e->tail)))
   8     e->distance = A[e->head->index][e->head->index];
   9   else if (leaf(e->head))
  10     {
  11       f = e->tail->parentEdge;
  12       g = siblingEdge(e);
  13       e->distance = 0.5*(A[e->head->index][g->head->index]
  14                          + A[e->head->index][f->head->index]
  15                          - A[g->head->index][f->head->index]);
  16     }
  17   else
  18     {
  19       f = e->head->leftEdge;
  20       g = e->head->rightEdge;
  21       e->distance = 0.5*(A[g->head->index][e->head->index]
  22                          + A[f->head->index][e->head->index]
  23                          - A[f->head->index][g->head->index]);
  24     }
  25 }
  26
  27 void BalWFint(edge *e, double **A)
  28 {
  29   int up, down, left, right;
  30   up = e->tail->index;
  31   down = (siblingEdge(e))->head->index;
  32   left = e->head->leftEdge->head->index;
  33   right = e->head->rightEdge->head->index;
  34   e->distance = 0.25*(A[up][left] + A[up][right] + A[left][down] + A[right][down]) - 0.5*(A[down][up] + A[left][right]);
  35 }
  36
  37 void assignBMEWeights(tree *T, double **A)
  38 {
  39   edge *e;
  40   e = depthFirstTraverse(T,NULL);
  41   while (NULL != e) {
  42     if ((leaf(e->head)) || (leaf(e->tail)))
  43       BalWFext(e,A);
  44     else
  45       BalWFint(e,A);
  46     e = depthFirstTraverse(T,e);
  47   }
  48 }
  49
  50 void BMEcalcDownAverage(tree *T, node *v, edge *e, double **D, double **A)
  51 {
  52   edge  *left, *right;
  53   if (leaf(e->head))
  54     A[e->head->index][v->index] = D[v->index2][e->head->index2];
  55   else
  56     {
  57       left = e->head->leftEdge;
  58       right = e->head->rightEdge;
  59       A[e->head->index][v->index] = 0.5 * A[left->head->index][v->index]
  60         + 0.5 * A[right->head->index][v->index];
  61     }
  62 }
  63
  64 void BMEcalcUpAverage(tree *T, node *v, edge *e, double **D, double **A)
  65 {
  66   edge *up,*down;
  67   if (T->root == e->tail)
  68     A[v->index][e->head->index] = D[v->index2][e->tail->index2];
  69   /*for now, use convention
  70     v->index first => looking up
  71     v->index second => looking down */
  72   else
  73     {
  74       up = e->tail->parentEdge;
  75       down = siblingEdge(e);
  76       A[v->index][e->head->index] = 0.5 * A[v->index][up->head->index]
  77         +0.5  * A[down->head->index][v->index];
  78     }
  79 }
  80
  81
  82 void BMEcalcNewvAverages(tree *T, node *v, double **D, double **A)
  83 {
  84   /*loop over edges*/
  85   /*depth-first search*/
  86   edge *e;
  87   e = NULL;
  88   e = depthFirstTraverse(T,e);  /*the downward averages need to be
  89                                   calculated from bottom to top */
  90   while(NULL != e)
  91     {
  92       BMEcalcDownAverage(T,v,e,D,A);
  93       e = depthFirstTraverse(T,e);
  94     }
  95
  96   e = topFirstTraverse(T,e);   /*the upward averages need to be calculated
  97                                  from top to bottom */
  98   while(NULL != e)
  99     {
 100       BMEcalcUpAverage(T,v,e,D,A);
 101       e = topFirstTraverse(T,e);
 102     }
 103 }
 104
 105
 106 /*update Pair updates A[nearEdge][farEdge] and makes recursive call to subtree
 107   beyond farEdge*/
 108 /*root is head or tail of edge being split, depending on direction toward
 109   v*/
 110 void updatePair(double **A, edge *nearEdge, edge *farEdge, node *v,
 111                 node *root, double dcoeff, int direction)
 112 {
 113   edge *sib;
 114   switch(direction) /*the various cases refer to where the new vertex has
 115                       been inserted, in relation to the edge nearEdge*/
 116     {
 117     case UP: /*this case is called when v has been inserted above
 118                or skew to farEdge*/
 119       /*do recursive calls first!*/
 120       if (NULL != farEdge->head->leftEdge)
 121         updatePair(A,nearEdge,farEdge->head->leftEdge,v,root,dcoeff,UP);
 122       if (NULL != farEdge->head->rightEdge)
 123         updatePair(A,nearEdge,farEdge->head->rightEdge,v,root,dcoeff,UP);
 124       A[farEdge->head->index][nearEdge->head->index] =
 125         A[nearEdge->head->index][farEdge->head->index]
 126         = A[farEdge->head->index][nearEdge->head->index]
 127         + dcoeff*A[farEdge->head->index][v->index]
 128         - dcoeff*A[farEdge->head->index][root->index];
 129       break;
 130     case DOWN: /*called when v has been inserted below farEdge*/
 131       if (NULL != farEdge->tail->parentEdge)
 132         updatePair(A,nearEdge,farEdge->tail->parentEdge,v,root,dcoeff,DOWN);
 133       sib = siblingEdge(farEdge);
 134       if (NULL != sib)
 135         updatePair(A,nearEdge,sib,v,root,dcoeff,UP);
 136       A[farEdge->head->index][nearEdge->head->index] =
 137         A[nearEdge->head->index][farEdge->head->index]
 138         = A[farEdge->head->index][nearEdge->head->index]
 139         + dcoeff*A[v->index][farEdge->head->index]
 140         - dcoeff*A[farEdge->head->index][root->index];
 141     }
 142 }
 143
 144 void updateSubTree(double **A, edge *nearEdge, node *v, node *root,
 145                    node *newNode, double dcoeff, int direction)
 146 {
 147   edge *sib;
 148   switch(direction)
 149     {
 150     case UP: /*newNode is above the edge nearEdge*/
 151       A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
 152       A[newNode->index][nearEdge->head->index] =
 153         A[nearEdge->head->index][newNode->index] =
 154         A[nearEdge->head->index][root->index];
 155       if (NULL != nearEdge->head->leftEdge)
 156         updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff, UP);
 157       if (NULL != nearEdge->head->rightEdge)
 158         updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff, UP);
 159       updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
 160       break;
 161     case DOWN: /*newNode is below the edge nearEdge*/
 162       A[nearEdge->head->index][v->index] = A[v->index][nearEdge->head->index];
 163       A[newNode->index][nearEdge->head->index] =
 164         A[nearEdge->head->index][newNode->index] =
 165         0.5*(A[nearEdge->head->index][root->index]
 166              + A[v->index][nearEdge->head->index]);
 167       sib = siblingEdge(nearEdge);
 168       if (NULL != sib)
 169         updateSubTree(A, sib, v, root, newNode, 0.5*dcoeff, SKEW);
 170       if (NULL != nearEdge->tail->parentEdge)
 171         updateSubTree(A, nearEdge->tail->parentEdge, v, root, newNode, 0.5*dcoeff, DOWN);
 172       updatePair(A, nearEdge, nearEdge, v, root, dcoeff, DOWN);
 173       break;
 174     case SKEW: /*newNode is neither above nor below nearEdge*/
 175       A[v->index][nearEdge->head->index] = A[nearEdge->head->index][v->index];
 176       A[newNode->index][nearEdge->head->index] =
 177         A[nearEdge->head->index][newNode->index] =
 178         0.5*(A[nearEdge->head->index][root->index] +
 179              A[nearEdge->head->index][v->index]);
 180       if (NULL != nearEdge->head->leftEdge)
 181         updateSubTree(A, nearEdge->head->leftEdge, v, root, newNode, 0.5*dcoeff,SKEW);
 182       if (NULL != nearEdge->head->rightEdge)
 183         updateSubTree(A, nearEdge->head->rightEdge, v, root, newNode, 0.5*dcoeff,SKEW);
 184       updatePair(A, nearEdge, nearEdge, v, root, dcoeff, UP);
 185     }
 186 }
 187
 188
 189 /*we update all the averages for nodes (u1,u2), where the insertion point of
 190   v is in "direction" from both u1 and u2 */
 191 /*The general idea is to proceed in a direction from those edges already corrected
 192  */
 193
 194 /*r is the root of the tree relative to the inserted node*/
 195
 196 void BMEupdateAveragesMatrix(double **A, edge *e, node *v,node *newNode)
 197 {
 198   edge *sib, *par, *left, *right;
 199   /*first, update the v,newNode entries*/
 200   A[newNode->index][newNode->index] = 0.5*(A[e->head->index][e->head->index]
 201                                            + A[v->index][e->head->index]);
 202   A[v->index][newNode->index] = A[newNode->index][v->index] =
 203     A[v->index][e->head->index];
 204   A[v->index][v->index] =
 205     0.5*(A[e->head->index][v->index] + A[v->index][e->head->index]);
 206   left = e->head->leftEdge;
 207   right = e->head->rightEdge;
 208   if (NULL != left)
 209     updateSubTree(A,left,v,e->head,newNode,0.25,UP); /*updates left and below*/
 210   if (NULL != right)
 211     updateSubTree(A,right,v,e->head,newNode,0.25,UP); /*updates right and below*/
 212   sib = siblingEdge(e);
 213   if (NULL != sib)
 214     updateSubTree(A,sib,v,e->head,newNode,0.25,SKEW); /*updates sib and below*/
 215   par = e->tail->parentEdge;
 216   if (NULL != par)
 217     updateSubTree(A,par,v,e->head,newNode,0.25,DOWN); /*updates par and above*/
 218
 219   /*must change values A[e->head][*] last, as they are used to update
 220     the rest of the matrix*/
 221   A[newNode->index][e->head->index] = A[e->head->index][newNode->index]
 222     = A[e->head->index][e->head->index];
 223   A[v->index][e->head->index] = A[e->head->index][v->index];
 224
 225   updatePair(A,e,e,v,e->head,0.5,UP); /*updates e->head fields only*/
 226 }
 227
 228 /*A is tree below sibling, B is tree below edge, C is tree above edge*/
 229 double wf3(double D_AB, double D_AC, double D_kB, double D_kC)
 230 {
 231   return(D_AC + D_kB - D_AB - D_kC);
 232 }
 233
 234 void BMEtestEdge(edge *e, node *v, double **A)
 235 {
 236   edge *up, *down;
 237   down = siblingEdge(e);
 238   up = e->tail->parentEdge;
 239   e->totalweight = wf3(A[e->head->index][down->head->index],
 240                       A[down->head->index][e->tail->index],
 241                       A[e->head->index][v->index],
 242                       A[v->index][e->tail->index])
 243     + up->totalweight;
 244 }
 245
 246 void BMEsplitEdge(tree *T, node *v, edge *e, double **A)
 247 {
 248   edge *newPendantEdge;
 249   edge *newInternalEdge;
 250   node *newNode;
 251   char nodeLabel[NODE_LABEL_LENGTH];
 252   char edgeLabel1[EDGE_LABEL_LENGTH];
 253   char edgeLabel2[EDGE_LABEL_LENGTH];
 254   snprintf(nodeLabel,1,"");
 255   //sprintf(edgeLabel1,"E%d",T->size);
 256   //sprintf(edgeLabel2,"E%d",T->size+1);
 257   snprintf(edgeLabel1,EDGE_LABEL_LENGTH,"E%d",T->size);
 258   snprintf(edgeLabel2,EDGE_LABEL_LENGTH,"E%d",T->size+1);
 259
 260   /*make the new node and edges*/
 261   newNode = makeNewNode(nodeLabel,T->size+1);
 262   newPendantEdge = makeEdge(edgeLabel1,newNode,v,0.0);
 263   newInternalEdge = makeEdge(edgeLabel2,newNode,e->head,0.0);
 264
 265   /*update the matrix of average distances*/
 266   BMEupdateAveragesMatrix(A,e,v,newNode);
 267
 268   /*put them in the correct topology*/
 269   newNode->parentEdge = e;
 270   e->head->parentEdge = newInternalEdge;
 271   v->parentEdge = newPendantEdge;
 272   e->head = newNode;
 273
 274   T->size = T->size + 2;
 275
 276   if (e->tail->leftEdge == e)
 277     /*actually this is totally arbitrary and probably unnecessary*/
 278     {
 279       newNode->leftEdge = newInternalEdge;
 280       newNode->rightEdge = newPendantEdge;
 281     }
 282   else
 283     {
 284       newNode->leftEdge = newInternalEdge;
 285       newNode->rightEdge = newPendantEdge;
 286     }
 287 }
 288
 289 tree *BMEaddSpecies(tree *T,node *v, double **D, double **A)
 290      /*the key function of the program addSpeices inserts
 291        the node v to the tree T.  It uses testEdge to see what the relative
 292        weight would be if v split a particular edge.  Once insertion point
 293        is found, v is added to T, and A is updated.  Edge weights
 294        are not assigned until entire tree is build*/
 295
 296 {
 297   tree *T_e;
 298   edge *e; /*loop variable*/
 299   edge *e_min; /*points to best edge seen thus far*/
 300   double w_min = 0.0;   /*used to keep track of tree weights*/
 301
 302   /*initialize variables as necessary*/
 303
 304   /*CASE 1: T is empty, v is the first node*/
 305   if (NULL == T)  /*create a tree with v as only vertex, no edges*/
 306     {
 307       T_e = newTree();
 308       T_e->root = v;
 309       /*note that we are rooting T arbitrarily at a leaf.
 310         T->root is not the phylogenetic root*/
 311       v->index = 0;
 312       T_e->size = 1;
 313       return(T_e);
 314     }
 315   /*CASE 2: T is a single-vertex tree*/
 316   if (1 == T->size)
 317         {
 318           v->index = 1;
 319           e = makeEdge("",T->root,v,0.0);
 320           //sprintf(e->label,"E1");
 321           snprintf(e->label,EDGE_LABEL_LENGTH,"E1");
 322           A[v->index][v->index] = D[v->index2][T->root->index2];
 323           T->root->leftEdge = v->parentEdge = e;
 324           T->size = 2;
 325           return(T);
 326         }
 327   /*CASE 3: T has at least two nodes and an edge.  Insert new node
 328     by breaking one of the edges*/
 329
 330   v->index = T->size;
 331   BMEcalcNewvAverages(T,v,D,A);
 332   /*calcNewvAverages will update A for the row and column
 333     include the node v.  Will do so using pre-existing averages in T and
 334     information from A,D*/
 335   e_min = T->root->leftEdge;
 336   e = e_min->head->leftEdge;
 337   while (NULL != e)
 338     {
 339       BMEtestEdge(e,v,A);
 340       /*testEdge tests weight of tree if loop variable
 341         e is the edge split, places this value in the e->totalweight field */
 342       if (e->totalweight < w_min)
 343         {
 344           e_min = e;
 345           w_min = e->totalweight;
 346         }
 347       e = topFirstTraverse(T,e);
 348     }
 349   /*e_min now points at the edge we want to split*/
 350 /*  if (verbose)
 351     printf("Inserting %s between %s and %s on %s\n",v->label,e_min->tail->label,
 352            e_min->head->label,e_min->label);*/
 353   BMEsplitEdge(T,v,e_min,A);
 354   return(T);
 355 }
 356
 357 /*calcUpAverages will ensure that A[e->head->index][f->head->index] is
 358   filled for any f >= g.  Works recursively*/
 359 void calcUpAverages(double **D, double **A, edge *e, edge *g)
 360 {
 361   node *u,*v;
 362   edge *s;
 363   if (!(leaf(g->tail)))
 364     {
 365       calcUpAverages(D,A,e,g->tail->parentEdge);
 366       s = siblingEdge(g);
 367       u = g->tail;
 368       v = s->head;
 369       A[e->head->index][g->head->index] = A[g->head->index][e->head->index]
 370         = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
 371     }
 372 }
 373
 374 void makeBMEAveragesTable(tree *T, double **D, double **A)
 375 {
 376   edge *e, *f, *exclude;
 377   node *u,*v;
 378   /*first, let's deal with the averages involving the root of T*/
 379   e = T->root->leftEdge;
 380   f = depthFirstTraverse(T,NULL);
 381   while (NULL != f) {
 382     if (leaf(f->head)) {
 383       A[e->head->index][f->head->index] = A[f->head->index][e->head->index]
 384         = D[e->tail->index2][f->head->index2];
 385         }
 386     else
 387       {
 388         u = f->head->leftEdge->head;
 389         v = f->head->rightEdge->head;
 390         A[e->head->index][f->head->index] = A[f->head->index][e->head->index]
 391           = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
 392       }
 393     f = depthFirstTraverse(T,f);
 394   }
 395  e = depthFirstTraverse(T,NULL);
 396   while (T->root->leftEdge != e) {
 397     f = exclude = e;
 398     while (T->root->leftEdge != f) {
 399       if (f == exclude)
 400         exclude = exclude->tail->parentEdge;
 401       else if (leaf(e->head))
 402         {
 403           if (leaf(f->head))
 404             A[e->head->index][f->head->index] =
 405               A[f->head->index][e->head->index]
 406               = D[e->head->index2][f->head->index2];
 407           else
 408             {
 409               u = f->head->leftEdge->head; /*since f is chosen using a
 410                                              depth-first search, other values
 411                                              have been calculated*/
 412               v = f->head->rightEdge->head;
 413               A[e->head->index][f->head->index]
 414                 = A[f->head->index][e->head->index]
 415                 = 0.5*(A[e->head->index][u->index] + A[e->head->index][v->index]);
 416             }
 417         }
 418       else
 419         {
 420           u = e->head->leftEdge->head;
 421           v = e->head->rightEdge->head;
 422           A[e->head->index][f->head->index] = A[f->head->index][e->head->index] = 0.5*(A[f->head->index][u->index] + A[f->head->index][v->index]);
 423         }
 424       f = depthFirstTraverse(T,f);
 425     }
 426     e = depthFirstTraverse(T,e);
 427   }
 428   e = depthFirstTraverse(T,NULL);
 429   while (T->root->leftEdge != e)
 430     {
 431       calcUpAverages(D,A,e,e); /*calculates averages for
 432                                  A[e->head->index][g->head->index] for
 433                                  any edge g in path from e to root of tree*/
 434       e = depthFirstTraverse(T,e);
 435     }
 436 } /*makeAveragesMatrix*/