src/me_balanced.c

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