9 //boolean leaf(node *v);
10 /*edge *siblingEdge(edge *e);
11 edge *depthFirstTraverse(tree *T, edge *e);
12 edge *findBottomLeft(edge *e);
13 edge *topFirstTraverse(tree *T, edge *e);
14 edge *moveUpRight(edge *e);
15 double wf(double lambda, double D_LR, double D_LU, double D_LD,
16 double D_RU, double D_RD, double D_DU);*/
17 /*NNI functions for unweighted OLS topological switches*/
19 /*fillTableUp fills all the entries in D associated with
20 e->head,f->head and those edges g->head above e->head*/
21 void fillTableUp(edge *e, edge *f, double **A, double **D, tree *T)
24 if (T->root == f->tail)
27 A[e->head->index][f->head->index] =
28 A[f->head->index][e->head->index] =
29 D[e->head->index2][f->tail->index2];
32 g = e->head->leftEdge;
33 h = e->head->rightEdge;
34 A[e->head->index][f->head->index] =
35 A[f->head->index][e->head->index] =
36 (g->bottomsize*A[f->head->index][g->head->index]
37 + h->bottomsize*A[f->head->index][h->head->index])
43 g = f->tail->parentEdge;
44 fillTableUp(e,g,A,D,T); /*recursive call*/
46 A[e->head->index][f->head->index] =
47 A[f->head->index][e->head->index] =
48 (g->topsize*A[e->head->index][g->head->index]
49 + h->bottomsize*A[e->head->index][h->head->index])/f->topsize;
54 void makeOLSAveragesTable(tree *T, double **D, double **A);
56 double **buildAveragesTable(tree *T, double **D)
60 A = (double **) malloc(T->size*sizeof(double *));
61 for(i = 0; i < T->size;i++)
63 A[i] = (double *) malloc(T->size*sizeof(double));
64 for(j=0;j<T->size;j++)
67 makeOLSAveragesTable(T,D,A);
71 double wf2(double lambda, double D_AD, double D_BC, double D_AC, double D_BD,
72 double D_AB, double D_CD)
75 weight = 0.5*(lambda*(D_AC + D_BD) + (1 - lambda)*(D_AD + D_BC)
80 int NNIEdgeTest(edge *e, tree *T, double **A, double *weight)
85 double D_LR, D_LU, D_LD, D_RD, D_RU, D_DU;
88 if ((leaf(e->tail)) || (leaf(e->head)))
90 lambda = (double *)malloc(3*sizeof(double));
91 a = e->tail->parentEdge->topsize;
94 c = e->head->leftEdge->bottomsize;
95 d = e->head->rightEdge->bottomsize;
97 lambda[0] = ((double) b*c + a*d)/((a + b)*(c+d));
98 lambda[1] = ((double) b*c + a*d)/((a + c)*(b+d));
99 lambda[2] = ((double) c*d + a*b)/((a + d)*(b+c));
101 D_LR = A[e->head->leftEdge->head->index][e->head->rightEdge->head->index];
102 D_LU = A[e->head->leftEdge->head->index][e->tail->index];
103 D_LD = A[e->head->leftEdge->head->index][f->head->index];
104 D_RU = A[e->head->rightEdge->head->index][e->tail->index];
105 D_RD = A[e->head->rightEdge->head->index][f->head->index];
106 D_DU = A[e->tail->index][f->head->index];
108 w0 = wf2(lambda[0],D_RU,D_LD,D_LU,D_RD,D_DU,D_LR);
109 w1 = wf2(lambda[1],D_RU,D_LD,D_DU,D_LR,D_LU,D_RD);
110 w2 = wf2(lambda[2],D_DU,D_LR,D_LU,D_RD,D_RU,D_LD);
114 if (w0 <= w2) /*w0 <= w1,w2*/
119 else /*w2 < w0 <= w1 */
124 printf("Swap across %s. ",e->label);
125 printf("Weight dropping by %lf.\n",w0 - w2);
126 printf("New weight should be %lf.\n",T->weight + w2 - w0);
131 else if (w2 <= w1) /*w2 <= w1 < w0*/
136 printf("Swap across %s. ",e->label);
137 printf("Weight dropping by %lf.\n",w0 - w2);
138 printf("New weight should be %lf.\n",T->weight + w2 - w0);
147 printf("Swap across %s. ",e->label);
148 printf("Weight dropping by %lf.\n",w0 - w1);
149 printf("New weight should be %lf.\n",T->weight + w1 - w0);
156 int *initPerm(int size);
158 void NNIupdateAverages(double **A, edge *e, edge *par, edge *skew,
159 edge *swap, edge *fixed, tree *T)
165 A[e->head->index][e->head->index] =
167 ((skew->bottomsize*A[skew->head->index][swap->head->index]
168 + fixed->bottomsize*A[fixed->head->index][swap->head->index])
171 ((skew->bottomsize*A[skew->head->index][par->head->index]
172 + fixed->bottomsize*A[fixed->head->index][par->head->index])
176 elooper = findBottomLeft(e); /*next, we loop over all the edges
180 A[e->head->index][elooper->head->index] =
181 A[elooper->head->index][v->index]
182 = (swap->bottomsize*A[elooper->head->index][swap->head->index] +
183 par->topsize*A[elooper->head->index][par->head->index])
185 elooper = depthFirstTraverse(T,elooper);
187 elooper = findBottomLeft(swap); /*next we loop over all the edges below and
189 while (swap != elooper)
191 A[e->head->index][elooper->head->index] =
192 A[elooper->head->index][e->head->index]
193 = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
194 fixed->bottomsize*A[elooper->head->index][fixed->head->index])
196 elooper = depthFirstTraverse(T,elooper);
198 /*now elooper = skew */
199 A[e->head->index][elooper->head->index] =
200 A[elooper->head->index][e->head->index]
201 = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
202 fixed->bottomsize* A[elooper->head->index][fixed->head->index])
205 /*finally, we loop over all the edges in the tree
206 on the far side of parEdge*/
207 elooper = T->root->leftEdge;
208 while ((elooper != swap) && (elooper != e)) /*start a top-first traversal*/
210 A[e->head->index][elooper->head->index] =
211 A[elooper->head->index][e->head->index]
212 = (skew->bottomsize * A[elooper->head->index][skew->head->index]
213 + fixed->bottomsize* A[elooper->head->index][fixed->head->index])
215 elooper = topFirstTraverse(T,elooper);
218 /*At this point, elooper = par.
219 We finish the top-first traversal, excluding the subtree below par*/
220 elooper = moveUpRight(par);
221 while (NULL != elooper)
223 A[e->head->index][elooper->head->index]
224 = A[elooper->head->index][e->head->index]
225 = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
226 fixed->bottomsize* A[elooper->head->index][fixed->head->index])
228 elooper = topFirstTraverse(T,elooper);
234 void NNItopSwitch(tree *T, edge *e, int direction, double **A)
240 printf("Branch swap across edge %s.\n",e->label);*/
242 if (LEFT == direction)
243 swap = e->head->leftEdge;
245 swap = e->head->rightEdge;
246 skew = siblingEdge(e);
247 fixed = siblingEdge(swap);
248 par = e->tail->parentEdge;
252 printf("Branch swap: switching edges %s and %s.\n",skew->label,swap->label);
254 /*perform topological switch by changing f from (u,b) to (v,b)
255 and g from (v,c) to (u,c), necessitatates also changing parent fields*/
257 swap->tail = e->tail;
258 skew->tail = e->head;
260 if (LEFT == direction)
261 e->head->leftEdge = skew;
263 e->head->rightEdge = skew;
264 if (skew == e->tail->rightEdge)
265 e->tail->rightEdge = swap;
267 e->tail->leftEdge = swap;
269 /*both topsize and bottomsize change for e, but nowhere else*/
271 e->topsize = par->topsize + swap->bottomsize;
272 e->bottomsize = fixed->bottomsize + skew->bottomsize;
273 NNIupdateAverages(A, e, par, skew, swap, fixed,T);
275 } /*end NNItopSwitch*/
277 void reHeapElement(int *p, int *q, double *v, int length, int i);
278 void pushHeap(int *p, int *q, double *v, int length, int i);
279 void popHeap(int *p, int *q, double *v, int length, int i);
282 void NNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray,
283 double *weights, int *location, int *possibleSwaps)
286 tloc = location[e->head->index+1];
287 location[e->head->index+1] =
288 NNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
289 if (NONE == location[e->head->index+1])
292 popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
297 pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
299 reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
303 void permInverse(int *p, int *q, int length);
305 int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
308 //void NNI(tree *T, double **avgDistArray, int *count)
309 void NNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
311 edge *e, *centerEdge;
318 p = initPerm(T->size+1);
319 q = initPerm(T->size+1);
320 edgeArray = (edge **) malloc((T->size+1)*sizeof(double));
321 weights = (double *) malloc((T->size+1)*sizeof(double));
322 location = (int *) malloc((T->size+1)*sizeof(int));
324 double epsilon = 0.0;
325 for (i=0; i<numSpecies; i++)
326 for (j=0; j<numSpecies; j++)
328 epsilon = (epsilon / (numSpecies * numSpecies)) * EPSILON;
330 for (i=0;i<T->size+1;i++)
335 e = findBottomLeft(T->root->leftEdge);
339 edgeArray[e->head->index+1] = e;
340 location[e->head->index+1] =
341 NNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
342 e = depthFirstTraverse(T,e);
344 possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
345 permInverse(p,q,T->size+1);
346 /*we put the negative values of weights into a heap, indexed by p
347 with the minimum value pointed to by p[1]*/
348 /*p[i] is index (in edgeArray) of edge with i-th position
349 in the heap, q[j] is the position of edge j in the heap */
350 while (weights[p[1]] + epsilon < 0)
352 centerEdge = edgeArray[p[1]];
354 T->weight = T->weight + weights[p[1]];
355 NNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
356 location[p[1]] = NONE;
357 weights[p[1]] = 0.0; /*after the NNI, this edge is in optimal
359 popHeap(p,q,weights,possibleSwaps--,1);
360 /*but we must retest the other four edges*/
361 e = centerEdge->head->leftEdge;
362 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
363 e = centerEdge->head->rightEdge;
364 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
365 e = siblingEdge(centerEdge);
366 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
367 e = centerEdge->tail->parentEdge;
368 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
376 void NNIwithoutMatrix(tree *T, double **D, int *count)
378 double **avgDistArray;
379 avgDistArray = buildAveragesTable(T,D);
380 NNI(T,avgDistArray,count);
383 void NNIWithPartialMatrix(tree *T,double **D,double **A,int *count)
385 makeOLSAveragesTable(T,D,A);