src/me_balanced.c

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