3 /* Copyright 2007 Vincent Lefort */
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
10 //boolean leaf(node *v);
11 /*edge *siblingEdge(edge *e);
12 edge *depthFirstTraverse(tree *T, edge *e);
13 edge *findBottomLeft(edge *e);
14 edge *topFirstTraverse(tree *T, edge *e);
15 edge *moveUpRight(edge *e);
16 double wf(double lambda, double D_LR, double D_LU, double D_LD,
17 double D_RU, double D_RD, double D_DU);*/
18 /*NNI functions for unweighted OLS topological switches*/
20 /*fillTableUp fills all the entries in D associated with
21 e->head,f->head and those edges g->head above e->head*/
22 void fillTableUp(edge *e, edge *f, double **A, double **D, tree *T)
25 if (T->root == f->tail)
28 A[e->head->index][f->head->index] =
29 A[f->head->index][e->head->index] =
30 D[e->head->index2][f->tail->index2];
33 g = e->head->leftEdge;
34 h = e->head->rightEdge;
35 A[e->head->index][f->head->index] =
36 A[f->head->index][e->head->index] =
37 (g->bottomsize*A[f->head->index][g->head->index]
38 + h->bottomsize*A[f->head->index][h->head->index])
44 g = f->tail->parentEdge;
45 fillTableUp(e,g,A,D,T); /*recursive call*/
47 A[e->head->index][f->head->index] =
48 A[f->head->index][e->head->index] =
49 (g->topsize*A[e->head->index][g->head->index]
50 + h->bottomsize*A[e->head->index][h->head->index])/f->topsize;
55 void makeOLSAveragesTable(tree *T, double **D, double **A);
57 double **buildAveragesTable(tree *T, double **D)
61 A = (double **) malloc(T->size*sizeof(double *));
62 for(i = 0; i < T->size;i++)
64 A[i] = (double *) malloc(T->size*sizeof(double));
65 for(j=0;j<T->size;j++)
68 makeOLSAveragesTable(T,D,A);
72 double wf2(double lambda, double D_AD, double D_BC, double D_AC, double D_BD,
73 double D_AB, double D_CD)
76 weight = 0.5*(lambda*(D_AC + D_BD) + (1 - lambda)*(D_AD + D_BC)
81 int NNIEdgeTest(edge *e, tree *T, double **A, double *weight)
86 double D_LR, D_LU, D_LD, D_RD, D_RU, D_DU;
89 if ((leaf(e->tail)) || (leaf(e->head)))
91 lambda = (double *)malloc(3*sizeof(double));
92 a = e->tail->parentEdge->topsize;
95 c = e->head->leftEdge->bottomsize;
96 d = e->head->rightEdge->bottomsize;
98 lambda[0] = ((double) b*c + a*d)/((a + b)*(c+d));
99 lambda[1] = ((double) b*c + a*d)/((a + c)*(b+d));
100 lambda[2] = ((double) c*d + a*b)/((a + d)*(b+c));
102 D_LR = A[e->head->leftEdge->head->index][e->head->rightEdge->head->index];
103 D_LU = A[e->head->leftEdge->head->index][e->tail->index];
104 D_LD = A[e->head->leftEdge->head->index][f->head->index];
105 D_RU = A[e->head->rightEdge->head->index][e->tail->index];
106 D_RD = A[e->head->rightEdge->head->index][f->head->index];
107 D_DU = A[e->tail->index][f->head->index];
109 w0 = wf2(lambda[0],D_RU,D_LD,D_LU,D_RD,D_DU,D_LR);
110 w1 = wf2(lambda[1],D_RU,D_LD,D_DU,D_LR,D_LU,D_RD);
111 w2 = wf2(lambda[2],D_DU,D_LR,D_LU,D_RD,D_RU,D_LD);
115 if (w0 <= w2) /*w0 <= w1,w2*/
120 else /*w2 < w0 <= w1 */
125 printf("Swap across %s. ",e->label);
126 printf("Weight dropping by %lf.\n",w0 - w2);
127 printf("New weight should be %lf.\n",T->weight + w2 - w0);
132 else if (w2 <= w1) /*w2 <= w1 < w0*/
137 printf("Swap across %s. ",e->label);
138 printf("Weight dropping by %lf.\n",w0 - w2);
139 printf("New weight should be %lf.\n",T->weight + w2 - w0);
148 printf("Swap across %s. ",e->label);
149 printf("Weight dropping by %lf.\n",w0 - w1);
150 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);
233 void NNItopSwitch(tree *T, edge *e, int direction, double **A)
239 printf("Branch swap across edge %s.\n",e->label);*/
241 if (LEFT == direction)
242 swap = e->head->leftEdge;
244 swap = e->head->rightEdge;
245 skew = siblingEdge(e);
246 fixed = siblingEdge(swap);
247 par = e->tail->parentEdge;
251 printf("Branch swap: switching edges %s and %s.\n",skew->label,swap->label);
253 /*perform topological switch by changing f from (u,b) to (v,b)
254 and g from (v,c) to (u,c), necessitatates also changing parent fields*/
256 swap->tail = e->tail;
257 skew->tail = e->head;
259 if (LEFT == direction)
260 e->head->leftEdge = skew;
262 e->head->rightEdge = skew;
263 if (skew == e->tail->rightEdge)
264 e->tail->rightEdge = swap;
266 e->tail->leftEdge = swap;
268 /*both topsize and bottomsize change for e, but nowhere else*/
270 e->topsize = par->topsize + swap->bottomsize;
271 e->bottomsize = fixed->bottomsize + skew->bottomsize;
272 NNIupdateAverages(A, e, par, skew, swap, fixed,T);
274 } /*end NNItopSwitch*/
276 void reHeapElement(int *p, int *q, double *v, int length, int i);
277 void pushHeap(int *p, int *q, double *v, int length, int i);
278 void popHeap(int *p, int *q, double *v, int length, int i);
281 void NNIRetestEdge(int *p, int *q, edge *e,tree *T, double **avgDistArray,
282 double *weights, int *location, int *possibleSwaps)
285 tloc = location[e->head->index+1];
286 location[e->head->index+1] =
287 NNIEdgeTest(e,T,avgDistArray,weights + e->head->index+1);
288 if (NONE == location[e->head->index+1])
291 popHeap(p,q,weights,(*possibleSwaps)--,q[e->head->index+1]);
296 pushHeap(p,q,weights,(*possibleSwaps)++,q[e->head->index+1]);
298 reHeapElement(p,q,weights,*possibleSwaps,q[e->head->index+1]);
302 void permInverse(int *p, int *q, int length);
304 int makeThreshHeap(int *p, int *q, double *v, int arraySize, double thresh);
307 //void NNI(tree *T, double **avgDistArray, int *count)
308 void NNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
310 edge *e, *centerEdge;
317 p = initPerm(T->size+1);
318 q = initPerm(T->size+1);
319 edgeArray = (edge **) malloc((T->size+1)*sizeof(double));
320 weights = (double *) malloc((T->size+1)*sizeof(double));
321 location = (int *) malloc((T->size+1)*sizeof(int));
323 double epsilon = 0.0;
324 for (i=0; i<numSpecies; i++)
325 for (j=0; j<numSpecies; j++)
327 epsilon = (epsilon / (numSpecies * numSpecies)) * EPSILON;
329 for (i=0;i<T->size+1;i++)
334 e = findBottomLeft(T->root->leftEdge);
338 edgeArray[e->head->index+1] = e;
339 location[e->head->index+1] =
340 NNIEdgeTest(e,T,avgDistArray,weights + e->head->index + 1);
341 e = depthFirstTraverse(T,e);
343 possibleSwaps = makeThreshHeap(p,q,weights,T->size+1,0.0);
344 permInverse(p,q,T->size+1);
345 /*we put the negative values of weights into a heap, indexed by p
346 with the minimum value pointed to by p[1]*/
347 /*p[i] is index (in edgeArray) of edge with i-th position
348 in the heap, q[j] is the position of edge j in the heap */
349 while (weights[p[1]] + epsilon < 0)
351 centerEdge = edgeArray[p[1]];
353 T->weight = T->weight + weights[p[1]];
354 NNItopSwitch(T,edgeArray[p[1]],location[p[1]],avgDistArray);
355 location[p[1]] = NONE;
356 weights[p[1]] = 0.0; /*after the NNI, this edge is in optimal
358 popHeap(p,q,weights,possibleSwaps--,1);
359 /*but we must retest the other four edges*/
360 e = centerEdge->head->leftEdge;
361 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
362 e = centerEdge->head->rightEdge;
363 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
364 e = siblingEdge(centerEdge);
365 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
366 e = centerEdge->tail->parentEdge;
367 NNIRetestEdge(p,q,e,T,avgDistArray,weights,location,&possibleSwaps);
375 void NNIwithoutMatrix(tree *T, double **D, int *count)
377 double **avgDistArray;
378 avgDistArray = buildAveragesTable(T,D);
379 NNI(T,avgDistArray,count);
382 void NNIWithPartialMatrix(tree *T,double **D,double **A,int *count)
384 makeOLSAveragesTable(T,D,A);