3 //boolean leaf(node *v);
4 /*edge *siblingEdge(edge *e);
5 edge *depthFirstTraverse(tree *T, edge *e);
6 edge *findBottomLeft(edge *e);
7 edge *topFirstTraverse(tree *T, edge *e);
8 edge *moveUpRight(edge *e);
9 double wf(double lambda, double D_LR, double D_LU, double D_LD,
10 double D_RU, double D_RD, double D_DU);*/
11 /*NNI functions for unweighted OLS topological switches*/
13 /*fillTableUp fills all the entries in D associated with
14 e->head,f->head and those edges g->head above e->head*/
15 void fillTableUp(edge *e, edge *f, double **A, double **D, tree *T)
18 if (T->root == f->tail)
21 A[e->head->index][f->head->index] =
22 A[f->head->index][e->head->index] =
23 D[e->head->index2][f->tail->index2];
26 g = e->head->leftEdge;
27 h = e->head->rightEdge;
28 A[e->head->index][f->head->index] =
29 A[f->head->index][e->head->index] =
30 (g->bottomsize*A[f->head->index][g->head->index]
31 + h->bottomsize*A[f->head->index][h->head->index])
37 g = f->tail->parentEdge;
38 fillTableUp(e,g,A,D,T); /*recursive call*/
40 A[e->head->index][f->head->index] =
41 A[f->head->index][e->head->index] =
42 (g->topsize*A[e->head->index][g->head->index]
43 + h->bottomsize*A[e->head->index][h->head->index])/f->topsize;
48 void makeOLSAveragesTable(tree *T, double **D, double **A);
50 double **buildAveragesTable(tree *T, double **D)
54 A = (double **) malloc(T->size*sizeof(double *));
55 for(i = 0; i < T->size;i++)
57 A[i] = (double *) malloc(T->size*sizeof(double));
58 for(j=0;j<T->size;j++)
61 makeOLSAveragesTable(T,D,A);
65 double wf2(double lambda, double D_AD, double D_BC, double D_AC, double D_BD,
66 double D_AB, double D_CD)
69 weight = 0.5*(lambda*(D_AC + D_BD) + (1 - lambda)*(D_AD + D_BC)
74 int NNIEdgeTest(edge *e, tree *T, double **A, double *weight)
79 double D_LR, D_LU, D_LD, D_RD, D_RU, D_DU;
82 if ((leaf(e->tail)) || (leaf(e->head)))
84 lambda = (double *)malloc(3*sizeof(double));
85 a = e->tail->parentEdge->topsize;
88 c = e->head->leftEdge->bottomsize;
89 d = e->head->rightEdge->bottomsize;
91 lambda[0] = ((double) b*c + a*d)/((a + b)*(c+d));
92 lambda[1] = ((double) b*c + a*d)/((a + c)*(b+d));
93 lambda[2] = ((double) c*d + a*b)/((a + d)*(b+c));
95 D_LR = A[e->head->leftEdge->head->index][e->head->rightEdge->head->index];
96 D_LU = A[e->head->leftEdge->head->index][e->tail->index];
97 D_LD = A[e->head->leftEdge->head->index][f->head->index];
98 D_RU = A[e->head->rightEdge->head->index][e->tail->index];
99 D_RD = A[e->head->rightEdge->head->index][f->head->index];
100 D_DU = A[e->tail->index][f->head->index];
102 w0 = wf2(lambda[0],D_RU,D_LD,D_LU,D_RD,D_DU,D_LR);
103 w1 = wf2(lambda[1],D_RU,D_LD,D_DU,D_LR,D_LU,D_RD);
104 w2 = wf2(lambda[2],D_DU,D_LR,D_LU,D_RD,D_RU,D_LD);
108 if (w0 <= w2) /*w0 <= w1,w2*/
113 else /*w2 < w0 <= w1 */
118 printf("Swap across %s. ",e->label);
119 printf("Weight dropping by %lf.\n",w0 - w2);
120 printf("New weight should be %lf.\n",T->weight + w2 - w0);
125 else if (w2 <= w1) /*w2 <= w1 < w0*/
130 printf("Swap across %s. ",e->label);
131 printf("Weight dropping by %lf.\n",w0 - w2);
132 printf("New weight should be %lf.\n",T->weight + w2 - w0);
141 printf("Swap across %s. ",e->label);
142 printf("Weight dropping by %lf.\n",w0 - w1);
143 printf("New weight should be %lf.\n",T->weight + w1 - w0);
149 int *initPerm(int size);
151 void NNIupdateAverages(double **A, edge *e, edge *par, edge *skew,
152 edge *swap, edge *fixed, tree *T)
158 A[e->head->index][e->head->index] =
160 ((skew->bottomsize*A[skew->head->index][swap->head->index]
161 + fixed->bottomsize*A[fixed->head->index][swap->head->index])
164 ((skew->bottomsize*A[skew->head->index][par->head->index]
165 + fixed->bottomsize*A[fixed->head->index][par->head->index])
169 elooper = findBottomLeft(e); /*next, we loop over all the edges
173 A[e->head->index][elooper->head->index] =
174 A[elooper->head->index][v->index]
175 = (swap->bottomsize*A[elooper->head->index][swap->head->index] +
176 par->topsize*A[elooper->head->index][par->head->index])
178 elooper = depthFirstTraverse(T,elooper);
180 elooper = findBottomLeft(swap); /*next we loop over all the edges below and
182 while (swap != elooper)
184 A[e->head->index][elooper->head->index] =
185 A[elooper->head->index][e->head->index]
186 = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
187 fixed->bottomsize*A[elooper->head->index][fixed->head->index])
189 elooper = depthFirstTraverse(T,elooper);
191 /*now elooper = skew */
192 A[e->head->index][elooper->head->index] =
193 A[elooper->head->index][e->head->index]
194 = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
195 fixed->bottomsize* A[elooper->head->index][fixed->head->index])
198 /*finally, we loop over all the edges in the tree
199 on the far side of parEdge*/
200 elooper = T->root->leftEdge;
201 while ((elooper != swap) && (elooper != e)) /*start a top-first traversal*/
203 A[e->head->index][elooper->head->index] =
204 A[elooper->head->index][e->head->index]
205 = (skew->bottomsize * A[elooper->head->index][skew->head->index]
206 + fixed->bottomsize* A[elooper->head->index][fixed->head->index])
208 elooper = topFirstTraverse(T,elooper);
211 /*At this point, elooper = par.
212 We finish the top-first traversal, excluding the subtree below par*/
213 elooper = moveUpRight(par);
214 while (NULL != elooper)
216 A[e->head->index][elooper->head->index]
217 = A[elooper->head->index][e->head->index]
218 = (skew->bottomsize * A[elooper->head->index][skew->head->index] +
219 fixed->bottomsize* A[elooper->head->index][fixed->head->index])
221 elooper = topFirstTraverse(T,elooper);
226 void NNItopSwitch(tree *T, edge *e, int direction, double **A)
232 printf("Branch swap across edge %s.\n",e->label);*/
234 if (LEFT == direction)
235 swap = e->head->leftEdge;
237 swap = e->head->rightEdge;
238 skew = siblingEdge(e);
239 fixed = siblingEdge(swap);
240 par = e->tail->parentEdge;
244 printf("Branch swap: switching edges %s and %s.\n",skew->label,swap->label);
246 /*perform topological switch by changing f from (u,b) to (v,b)
247 and g from (v,c) to (u,c), necessitatates also changing parent fields*/
249 swap->tail = e->tail;
250 skew->tail = e->head;
252 if (LEFT == direction)
253 e->head->leftEdge = skew;
255 e->head->rightEdge = skew;
256 if (skew == e->tail->rightEdge)
257 e->tail->rightEdge = swap;
259 e->tail->leftEdge = swap;
261 /*both topsize and bottomsize change for e, but nowhere else*/
263 e->topsize = par->topsize + swap->bottomsize;
264 e->bottomsize = fixed->bottomsize + skew->bottomsize;
265 NNIupdateAverages(A, e, par, skew, swap, fixed,T);
267 } /*end NNItopSwitch*/
269 void reHeapElement(int *p, int *q, double *v, int length, int i);
270 void pushHeap(int *p, int *q, double *v, int length, int i);
271 void popHeap(int *p, int *q, double *v, int length, int i);
274 void NNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray,
275 double *weights, int *location, int *possibleSwaps)
278 tloc = location[e->head->index+1];
279 location[e->head->index+1] =
280 NNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
281 if (NONE == location[e->head->index+1])
284 popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
289 pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
291 reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
295 void permInverse(int *p, int *q, int length);
297 int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
300 //void NNI(tree *T, double **avgDistArray, int *count)
301 void NNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
303 edge *e, *centerEdge;
310 p = initPerm(T->size+1);
311 q = initPerm(T->size+1);
312 edgeArray = (edge **) malloc((T->size+1)*sizeof(double));
313 weights = (double *) malloc((T->size+1)*sizeof(double));
314 location = (int *) malloc((T->size+1)*sizeof(int));
316 double epsilon = 0.0;
317 for (i=0; i<numSpecies; i++)
318 for (j=0; j<numSpecies; j++)
320 epsilon = (epsilon / (numSpecies * numSpecies)) * EPSILON;
322 for (i=0;i<T->size+1;i++)
327 e = findBottomLeft(T->root->leftEdge);
331 edgeArray[e->head->index+1] = e;
332 location[e->head->index+1] =
333 NNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
334 e = depthFirstTraverse(T,e);
336 possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
337 permInverse(p,q,T->size+1);
338 /*we put the negative values of weights into a heap, indexed by p
339 with the minimum value pointed to by p[1]*/
340 /*p[i] is index (in edgeArray) of edge with i-th position
341 in the heap, q[j] is the position of edge j in the heap */
342 while (weights[p[1]] + epsilon < 0)
344 centerEdge = edgeArray[p[1]];
346 T->weight = T->weight + weights[p[1]];
347 NNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
348 location[p[1]] = NONE;
349 weights[p[1]] = 0.0; /*after the NNI, this edge is in optimal
351 popHeap(p,q,weights,possibleSwaps--,1);
352 /*but we must retest the other four edges*/
353 e = centerEdge->head->leftEdge;
354 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
355 e = centerEdge->head->rightEdge;
356 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
357 e = siblingEdge(centerEdge);
358 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
359 e = centerEdge->tail->parentEdge;
360 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
368 void NNIwithoutMatrix(tree *T, double **D, int *count)
370 double **avgDistArray;
371 avgDistArray = buildAveragesTable(T,D);
372 NNI(T,avgDistArray,count);
375 void NNIWithPartialMatrix(tree *T,double **D,double **A,int *count)
377 makeOLSAveragesTable(T,D,A);