X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=src%2FbNNI.c;h=9363ea5c5d0308d7b150dbfecc45ac2068bc4409;hb=f3426364b40c7c0e6aadf6ea2690716425abdfc9;hp=fe7fcce3623fa3a914069d3997ccbf4e02a6e944;hpb=3f528bd2c47301b64ed03bd28039bbf22a7510ad;p=ape.git diff --git a/src/bNNI.c b/src/bNNI.c index fe7fcce..9363ea5 100644 --- a/src/bNNI.c +++ b/src/bNNI.c @@ -1,9 +1,3 @@ -/*#include -#include -#include -#include "graph.h" -#include "main.h" -*/ #include "me.h" /*boolean leaf(node *v); @@ -27,17 +21,17 @@ void pushHeap(int *p, int *q, double *v, int length, int i); void popHeap(int *p, int *q, double *v, int length, int i); -void bNNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray, +void bNNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray, double *weights, int *location, int *possibleSwaps) { int tloc; tloc = location[e->head->index+1]; - location[e->head->index+1] = + location[e->head->index+1] = bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1); if (NONE == location[e->head->index+1]) { if (NONE != tloc) - popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]); + popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]); } else { @@ -74,7 +68,7 @@ void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies edgeArray = (edge **) malloc((T->size+1)*sizeof(double)); weights = (double *) malloc((T->size+1)*sizeof(double)); location = (int *) malloc((T->size+1)*sizeof(int)); - + double epsilon = 0.0; for (i=0; iroot->leftEdge); + e = findBottomLeft(T->root->leftEdge); while (NULL != e) { edgeArray[e->head->index+1] = e; - location[e->head->index+1] = + location[e->head->index+1] = bNNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1); e = depthFirstTraverse(T,e); - } + } possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0); permInverse(p,q,T->size+1); /*we put the negative values of weights into a heap, indexed by p with the minimum value pointed to by p[1]*/ - /*p[i] is index (in edgeArray) of edge with i-th position + /*p[i] is index (in edgeArray) of edge with i-th position in the heap, q[j] is the position of edge j in the heap */ while (weights[p[1]] + epsilon < 0) { @@ -156,10 +150,10 @@ void updateSubTreeAfterNNI(double **A, node *v, edge *rootEdge, node *closer, no updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,UP); updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP); sib = siblingEdge(v->parentEdge); - A[rootEdge->head->index][v->index] = - A[v->index][rootEdge->head->index] = + A[rootEdge->head->index][v->index] = + A[v->index][rootEdge->head->index] = 0.5*A[rootEdge->head->index][sib->head->index] + - 0.5*A[rootEdge->head->index][v->parentEdge->tail->index]; + 0.5*A[rootEdge->head->index][v->parentEdge->tail->index]; break; case DOWN: /*rootEdge is above the center edge of the NNI*/ sib = siblingEdge(rootEdge); @@ -168,8 +162,8 @@ void updateSubTreeAfterNNI(double **A, node *v, edge *rootEdge, node *closer, no if (NULL != rootEdge->tail->parentEdge) updateSubTreeAfterNNI(A, v, rootEdge->tail->parentEdge, closer, further, 0.5*dcoeff, DOWN); updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, DOWN); - A[rootEdge->head->index][v->index] = - A[v->index][rootEdge->head->index] = + A[rootEdge->head->index][v->index] = + A[v->index][rootEdge->head->index] = 0.5*A[rootEdge->head->index][v->leftEdge->head->index] + 0.5*A[rootEdge->head->index][v->rightEdge->head->index]; break; @@ -179,8 +173,8 @@ void updateSubTreeAfterNNI(double **A, node *v, edge *rootEdge, node *closer, no if (NULL != rootEdge->head->rightEdge) updateSubTreeAfterNNI(A, v, rootEdge->head->rightEdge, closer, further, 0.5*dcoeff,SKEW); updatePair(A, rootEdge, rootEdge, closer, further, dcoeff, UP); - A[rootEdge->head->index][v->index] = - A[v->index][rootEdge->head->index] = + A[rootEdge->head->index][v->index] = + A[v->index][rootEdge->head->index] = 0.5*A[rootEdge->head->index][v->leftEdge->head->index] + 0.5*A[rootEdge->head->index][v->rightEdge->head->index]; break; @@ -188,18 +182,17 @@ void updateSubTreeAfterNNI(double **A, node *v, edge *rootEdge, node *closer, no } /*swapping across edge whose head is v*/ -void bNNIupdateAverages(double **A, node *v, edge *par, edge *skew, +void bNNIupdateAverages(double **A, node *v, edge *par, edge *skew, edge *swap, edge *fixed) -{ - A[v->index][v->index] = 0.25*(A[fixed->head->index][par->head->index] + - A[fixed->head->index][swap->head->index] + - A[skew->head->index][par->head->index] + +{ + A[v->index][v->index] = 0.25*(A[fixed->head->index][par->head->index] + + A[fixed->head->index][swap->head->index] + + A[skew->head->index][par->head->index] + A[skew->head->index][swap->head->index]); updateSubTreeAfterNNI(A, v, fixed, skew->head, swap->head, 0.25, UP); updateSubTreeAfterNNI(A, v, par, swap->head, skew->head, 0.25, DOWN); - updateSubTreeAfterNNI(A, v, skew, fixed->head, par->head, 0.25, UP); + updateSubTreeAfterNNI(A, v, skew, fixed->head, par->head, 0.25, UP); updateSubTreeAfterNNI(A, v, swap, par->head, fixed->head, 0.25, SKEW); - } @@ -309,11 +302,10 @@ int bNNIEdgeTest(edge *e, tree *T, double **A, double *weight) printf("Weight dropping by %lf.\n",w0 - w1); printf("New weight should be %lf.\n",T->weight + w1 - w0); }*/ - return(LEFT); + return(LEFT); } } - /*limitedFillTableUp fills all the entries in D associated with e->head,f->head and those edges g->head above e->head, working recursively and stopping when trigger is reached*/ @@ -324,8 +316,7 @@ void limitedFillTableUp(edge *e, edge *f, double **A, edge *trigger) if (f != trigger) limitedFillTableUp(e,g,A,trigger); h = siblingEdge(f); - A[e->head->index][f->head->index] = - A[f->head->index][e->head->index] = - 0.5*(A[e->head->index][g->head->index] + A[e->head->index][h->head->index]); + A[e->head->index][f->head->index] = + A[f->head->index][e->head->index] = + 0.5*(A[e->head->index][g->head->index] + A[e->head->index][h->head->index]); } -