1 | |
---|
2 | #include "phylip.h" |
---|
3 | #include "disc.h" |
---|
4 | #include "wagner.h" |
---|
5 | |
---|
6 | /* version 3.6. (c) Copyright 1993-2002 by the University of Washington. |
---|
7 | Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe. |
---|
8 | Permission is granted to copy and use this program provided no fee is |
---|
9 | charged for it and provided that this copyright notice is not removed. */ |
---|
10 | |
---|
11 | #define maxtrees 100 /* maximum number of trees to be printed out */ |
---|
12 | #define often 100 /* how often to notify how many trees examined */ |
---|
13 | #define many 1000 /* how many multiples of howoften before stop */ |
---|
14 | |
---|
15 | typedef long *treenumbers; |
---|
16 | typedef double *valptr; |
---|
17 | typedef long *placeptr; |
---|
18 | |
---|
19 | #ifndef OLDC |
---|
20 | /* function prototypes */ |
---|
21 | void getoptions(void); |
---|
22 | void allocrest(void); |
---|
23 | void doinit(void); |
---|
24 | void inputoptions(void); |
---|
25 | void doinput(void); |
---|
26 | void supplement(bitptr); |
---|
27 | void evaluate(node2 *); |
---|
28 | void addtraverse(node2 *,node2 *,node2 *,long *,long *,valptr,placeptr); |
---|
29 | void addit(long); |
---|
30 | void reroot(node2 *); |
---|
31 | void describe(void); |
---|
32 | void maketree(void); |
---|
33 | /* function prototypes */ |
---|
34 | #endif |
---|
35 | |
---|
36 | Char infilename[FNMLNGTH], outfilename[FNMLNGTH], outtreename[FNMLNGTH], weightfilename[FNMLNGTH], ancfilename[FNMLNGTH], mixfilename[FNMLNGTH]; |
---|
37 | node2 *root; |
---|
38 | long outgrno, rno, howmany, howoften, col, msets, ith; |
---|
39 | /* outgrno indicates outgroup */ |
---|
40 | |
---|
41 | boolean weights, thresh, ancvar, questions, allsokal, allwagner, |
---|
42 | mixture, simple, trout, noroot, didreroot, outgropt, |
---|
43 | progress, treeprint, stepbox, ancseq, mulsets, firstset; |
---|
44 | boolean *ancone, *anczero, *ancone0, *anczero0, justwts; |
---|
45 | pointptr2 treenode; /* pointers to all nodes in tree */ |
---|
46 | double fracdone, fracinc; |
---|
47 | double threshold; |
---|
48 | double *threshwt; |
---|
49 | bitptr wagner, wagner0; |
---|
50 | boolean *added; |
---|
51 | Char *guess; |
---|
52 | steptr numsteps; |
---|
53 | long **bestorders, **bestrees; |
---|
54 | steptr numsone, numszero; |
---|
55 | gbit *garbage; |
---|
56 | |
---|
57 | long examined, mults; |
---|
58 | boolean firsttime, done, full; |
---|
59 | double like, bestyet; |
---|
60 | treenumbers current, order; |
---|
61 | long fullset; |
---|
62 | bitptr zeroanc, oneanc; |
---|
63 | bitptr suppsteps; |
---|
64 | |
---|
65 | |
---|
66 | void getoptions() |
---|
67 | { |
---|
68 | /* interactively set options */ |
---|
69 | long loopcount, loopcount2; |
---|
70 | Char ch, ch2; |
---|
71 | |
---|
72 | fprintf(outfile, "\nPenny algorithm, version %s\n",VERSION); |
---|
73 | fprintf(outfile, " branch-and-bound to find all"); |
---|
74 | fprintf(outfile, " most parsimonious trees\n\n"); |
---|
75 | howoften = often; |
---|
76 | howmany = many; |
---|
77 | outgrno = 1; |
---|
78 | outgropt = false; |
---|
79 | simple = true; |
---|
80 | thresh = false; |
---|
81 | threshold = spp; |
---|
82 | trout = true; |
---|
83 | weights = false; |
---|
84 | justwts = false; |
---|
85 | ancvar = false; |
---|
86 | allsokal = false; |
---|
87 | allwagner = true; |
---|
88 | mixture = false; |
---|
89 | printdata = false; |
---|
90 | progress = true; |
---|
91 | treeprint = true; |
---|
92 | stepbox = false; |
---|
93 | ancseq = false; |
---|
94 | loopcount = 0; |
---|
95 | for(;;) { |
---|
96 | cleerhome(); |
---|
97 | printf("\nPenny algorithm, version %s\n",VERSION); |
---|
98 | printf(" branch-and-bound to find all most parsimonious trees\n\n"); |
---|
99 | printf("Settings for this run:\n"); |
---|
100 | printf(" X Use Mixed method? %s\n", |
---|
101 | mixture ? "Yes" : "No"); |
---|
102 | printf(" P Parsimony method? %s\n", |
---|
103 | (allwagner && !mixture) ? "Wagner" : |
---|
104 | (!(allwagner || mixture)) ? "Camin-Sokal" : "(methods in mixture)"); |
---|
105 | printf(" F How often to report, in trees:%5ld\n",howoften); |
---|
106 | printf(" H How many groups of%5ld trees:%6ld\n",howoften,howmany); |
---|
107 | printf(" O Outgroup root?"); |
---|
108 | if (outgropt) |
---|
109 | printf(" Yes, at species number%3ld\n", outgrno); |
---|
110 | else |
---|
111 | printf(" No, use as outgroup species%3ld\n", outgrno); |
---|
112 | printf(" S Branch and bound is simple? %s\n", |
---|
113 | simple ? "Yes" : "No. reconsiders order of species"); |
---|
114 | printf(" T Use Threshold parsimony?"); |
---|
115 | if (thresh) |
---|
116 | printf(" Yes, count steps up to%4.1f per char.\n", threshold); |
---|
117 | else |
---|
118 | printf(" No, use ordinary parsimony\n"); |
---|
119 | printf(" A Use ancestral states in input file? %s\n", |
---|
120 | ancvar ? "Yes" : "No"); |
---|
121 | printf(" W Sites weighted? %s\n", |
---|
122 | (weights ? "Yes" : "No")); |
---|
123 | printf(" M Analyze multiple data sets?"); |
---|
124 | if (mulsets) |
---|
125 | printf(" Yes, %2ld %s\n", msets, |
---|
126 | (justwts ? "sets of weights" : "data sets")); |
---|
127 | else |
---|
128 | printf(" No\n"); |
---|
129 | printf(" 0 Terminal type (IBM PC, ANSI, none)? %s\n", |
---|
130 | ibmpc ? "IBM PC" : ansi ? "ANSI" : "(none)"); |
---|
131 | printf(" 1 Print out the data at start of run %s\n", |
---|
132 | printdata ? "Yes" : "No"); |
---|
133 | printf(" 2 Print indications of progress of run %s\n", |
---|
134 | progress ? "Yes" : "No"); |
---|
135 | printf(" 3 Print out tree %s\n", |
---|
136 | treeprint ? "Yes" : "No"); |
---|
137 | printf(" 4 Print out steps in each character %s\n", |
---|
138 | stepbox ? "Yes" : "No"); |
---|
139 | printf(" 5 Print states at all nodes of tree %s\n", |
---|
140 | ancseq ? "Yes" : "No"); |
---|
141 | printf(" 6 Write out trees onto tree file? %s\n", |
---|
142 | trout ? "Yes" : "No"); |
---|
143 | if(weights && justwts){ |
---|
144 | printf( |
---|
145 | "WARNING: W option and Multiple Weights options are both on. "); |
---|
146 | printf( |
---|
147 | "The W menu option is unnecessary and has no additional effect. \n"); |
---|
148 | } |
---|
149 | printf("\nAre these settings correct?"); |
---|
150 | printf(" (type Y or the letter for one to change)\n"); |
---|
151 | #ifdef WIN32 |
---|
152 | phyFillScreenColor(); |
---|
153 | #endif |
---|
154 | scanf("%c%*[^\n]", &ch); |
---|
155 | getchar(); |
---|
156 | uppercase(&ch); |
---|
157 | if (ch == 'Y') |
---|
158 | break; |
---|
159 | if (strchr("WHFSOMPATX1234560",ch) != NULL){ |
---|
160 | switch (ch) { |
---|
161 | |
---|
162 | case 'X': |
---|
163 | mixture = !mixture; |
---|
164 | break; |
---|
165 | |
---|
166 | case 'P': |
---|
167 | allwagner = !allwagner; |
---|
168 | break; |
---|
169 | |
---|
170 | case 'A': |
---|
171 | ancvar = !ancvar; |
---|
172 | break; |
---|
173 | |
---|
174 | case 'H': |
---|
175 | inithowmany(&howmany, howoften); |
---|
176 | break; |
---|
177 | |
---|
178 | case 'F': |
---|
179 | inithowoften(&howoften); |
---|
180 | break; |
---|
181 | |
---|
182 | case 'S': |
---|
183 | simple = !simple; |
---|
184 | break; |
---|
185 | |
---|
186 | case 'O': |
---|
187 | outgropt = !outgropt; |
---|
188 | if (outgropt) |
---|
189 | initoutgroup(&outgrno, spp); |
---|
190 | else outgrno = 1; |
---|
191 | break; |
---|
192 | |
---|
193 | case 'T': |
---|
194 | thresh = !thresh; |
---|
195 | if (thresh) |
---|
196 | initthreshold(&threshold); |
---|
197 | break; |
---|
198 | |
---|
199 | case 'W': |
---|
200 | weights = !weights; |
---|
201 | break; |
---|
202 | |
---|
203 | case 'M': |
---|
204 | mulsets = !mulsets; |
---|
205 | if (mulsets){ |
---|
206 | printf("Multiple data sets or multiple weights?"); |
---|
207 | loopcount2 = 0; |
---|
208 | do { |
---|
209 | printf(" (type D or W)\n"); |
---|
210 | #ifdef WIN32 |
---|
211 | phyFillScreenColor(); |
---|
212 | #endif |
---|
213 | scanf("%c%*[^\n]", &ch2); |
---|
214 | getchar(); |
---|
215 | if (ch2 == '\n') |
---|
216 | ch2 = ' '; |
---|
217 | uppercase(&ch2); |
---|
218 | countup(&loopcount2, 10); |
---|
219 | } while ((ch2 != 'W') && (ch2 != 'D')); |
---|
220 | justwts = (ch2 == 'W'); |
---|
221 | if (justwts) |
---|
222 | justweights(&msets); |
---|
223 | else |
---|
224 | initdatasets(&msets); |
---|
225 | } |
---|
226 | break; |
---|
227 | |
---|
228 | case '0': |
---|
229 | initterminal(&ibmpc, &ansi); |
---|
230 | break; |
---|
231 | |
---|
232 | case '1': |
---|
233 | printdata = !printdata; |
---|
234 | break; |
---|
235 | |
---|
236 | case '2': |
---|
237 | progress = !progress; |
---|
238 | break; |
---|
239 | |
---|
240 | case '3': |
---|
241 | treeprint = !treeprint; |
---|
242 | break; |
---|
243 | |
---|
244 | case '4': |
---|
245 | stepbox = !stepbox; |
---|
246 | break; |
---|
247 | |
---|
248 | case '5': |
---|
249 | ancseq = !ancseq; |
---|
250 | break; |
---|
251 | |
---|
252 | case '6': |
---|
253 | trout = !trout; |
---|
254 | break; |
---|
255 | } |
---|
256 | } else |
---|
257 | printf("Not a possible option!\n"); |
---|
258 | countup(&loopcount, 100); |
---|
259 | } |
---|
260 | allsokal = (!allwagner && !mixture); |
---|
261 | } /* getoptions */ |
---|
262 | |
---|
263 | |
---|
264 | void allocrest() |
---|
265 | { |
---|
266 | long i; |
---|
267 | |
---|
268 | weight = (steptr)Malloc(chars*sizeof(steptr)); |
---|
269 | threshwt = (double *)Malloc(chars*sizeof(double)); |
---|
270 | bestorders = (long **)Malloc(maxtrees*sizeof(long *)); |
---|
271 | bestrees = (long **)Malloc(maxtrees*sizeof(long *)); |
---|
272 | for (i = 1; i <= maxtrees; i++) { |
---|
273 | bestorders[i - 1] = (long *)Malloc(spp*sizeof(long)); |
---|
274 | bestrees[i - 1] = (long *)Malloc(spp*sizeof(long)); |
---|
275 | } |
---|
276 | numsteps = (steptr)Malloc(chars*sizeof(steptr)); |
---|
277 | guess = (Char *)Malloc(chars*sizeof(Char)); |
---|
278 | numszero = (steptr)Malloc(chars*sizeof(steptr)); |
---|
279 | numsone = (steptr)Malloc(chars*sizeof(steptr)); |
---|
280 | current = (treenumbers)Malloc(spp*sizeof(long)); |
---|
281 | order = (treenumbers)Malloc(spp*sizeof(long)); |
---|
282 | nayme = (naym *)Malloc(spp*sizeof(naym)); |
---|
283 | added = (boolean *)Malloc(nonodes*sizeof(boolean)); |
---|
284 | ancone = (boolean *)Malloc(chars*sizeof(boolean)); |
---|
285 | anczero = (boolean *)Malloc(chars*sizeof(boolean)); |
---|
286 | ancone0 = (boolean *)Malloc(chars*sizeof(boolean)); |
---|
287 | anczero0 = (boolean *)Malloc(chars*sizeof(boolean)); |
---|
288 | wagner = (bitptr)Malloc(words*sizeof(long)); |
---|
289 | wagner0 = (bitptr)Malloc(words*sizeof(long)); |
---|
290 | zeroanc = (bitptr)Malloc(words*sizeof(long)); |
---|
291 | oneanc = (bitptr)Malloc(words*sizeof(long)); |
---|
292 | suppsteps = (bitptr)Malloc(words*sizeof(long)); |
---|
293 | extras = (steptr)Malloc(chars*sizeof(steptr)); |
---|
294 | } /* allocrest */ |
---|
295 | |
---|
296 | |
---|
297 | void doinit() |
---|
298 | { |
---|
299 | /* initializes variables */ |
---|
300 | |
---|
301 | inputnumbers(&spp, &chars, &nonodes, 1); |
---|
302 | words = chars / bits + 1; |
---|
303 | getoptions(); |
---|
304 | if (printdata) |
---|
305 | fprintf(outfile, "%2ld species, %3ld characters\n", spp, chars); |
---|
306 | alloctree2(&treenode); |
---|
307 | setuptree2(treenode); |
---|
308 | allocrest(); |
---|
309 | } /* doinit */ |
---|
310 | |
---|
311 | void inputoptions() |
---|
312 | { |
---|
313 | /* input the information on the options */ |
---|
314 | long i; |
---|
315 | if(justwts){ |
---|
316 | if(firstset){ |
---|
317 | scan_eoln(infile); |
---|
318 | if (ancvar) { |
---|
319 | inputancestorsnew(anczero0, ancone0); |
---|
320 | } |
---|
321 | if (mixture) { |
---|
322 | inputmixturenew(wagner0); |
---|
323 | } |
---|
324 | } |
---|
325 | for (i = 0; i < (chars); i++) |
---|
326 | weight[i] = 1; |
---|
327 | inputweights(chars, weight, &weights); |
---|
328 | for (i = 0; i < (words); i++) { |
---|
329 | if (mixture) |
---|
330 | wagner[i] = wagner0[i]; |
---|
331 | else if (allsokal) |
---|
332 | wagner[i] = 0; |
---|
333 | else |
---|
334 | wagner[i] = (1L << (bits + 1)) - (1L << 1); |
---|
335 | } |
---|
336 | } |
---|
337 | else { |
---|
338 | if (!firstset) { |
---|
339 | samenumsp(&chars, ith); |
---|
340 | } |
---|
341 | scan_eoln(infile); |
---|
342 | for (i = 0; i < (chars); i++) |
---|
343 | weight[i] = 1; |
---|
344 | if (ancvar) { |
---|
345 | inputancestorsnew(anczero0, ancone0); |
---|
346 | } |
---|
347 | if (mixture) { |
---|
348 | inputmixturenew(wagner0); |
---|
349 | } |
---|
350 | if (weights) |
---|
351 | inputweights(chars, weight, &weights); |
---|
352 | for (i = 0; i < (words); i++) { |
---|
353 | if (mixture) |
---|
354 | wagner[i] = wagner0[i]; |
---|
355 | else if (allsokal) |
---|
356 | wagner[i] = 0; |
---|
357 | else |
---|
358 | wagner[i] = (1L << (bits + 1)) - (1L << 1); |
---|
359 | } |
---|
360 | } |
---|
361 | for (i = 0; i < (chars); i++) { |
---|
362 | if (!ancvar) { |
---|
363 | anczero[i] = true; |
---|
364 | ancone[i] = (((1L << (i % bits + 1)) & wagner[i / bits]) != 0); |
---|
365 | } else { |
---|
366 | anczero[i] = anczero0[i]; |
---|
367 | ancone[i] = ancone0[i]; |
---|
368 | } |
---|
369 | } |
---|
370 | noroot = true; |
---|
371 | questions = false; |
---|
372 | for (i = 0; i < (chars); i++) { |
---|
373 | if (weight[i] > 0) { |
---|
374 | noroot = (noroot && ancone[i] && anczero[i] && |
---|
375 | ((((1L << (i % bits + 1)) & wagner[i / bits]) != 0) |
---|
376 | || threshold <= 2.0)); |
---|
377 | } |
---|
378 | questions = (questions || (ancone[i] && anczero[i])); |
---|
379 | threshwt[i] = threshold * weight[i]; |
---|
380 | } |
---|
381 | } /* inputoptions */ |
---|
382 | |
---|
383 | |
---|
384 | void doinput() |
---|
385 | { |
---|
386 | /* reads the input data */ |
---|
387 | inputoptions(); |
---|
388 | if(!justwts || firstset) |
---|
389 | inputdata2(treenode); |
---|
390 | } /* doinput */ |
---|
391 | |
---|
392 | |
---|
393 | void supplement(bitptr local_suppsteps) |
---|
394 | { |
---|
395 | /* determine minimum number of steps more which will |
---|
396 | be added when rest of species are put in tree */ |
---|
397 | long i, j, k, l; |
---|
398 | long defone, defzero, a; |
---|
399 | |
---|
400 | k = 0; |
---|
401 | for (i = 0; i < (words); i++) { |
---|
402 | defone = 0; |
---|
403 | defzero = 0; |
---|
404 | a = 0; |
---|
405 | for (l = 1; l <= bits; l++) { |
---|
406 | k++; |
---|
407 | if (k <= chars) { |
---|
408 | if (!ancone[k - 1]) |
---|
409 | defzero = ((long)defzero) | (1L << l); |
---|
410 | if (!anczero[k - 1]) |
---|
411 | defone = ((long)defone) | (1L << l); |
---|
412 | } |
---|
413 | } |
---|
414 | for (j = 0; j < (spp); j++) { |
---|
415 | defone |= treenode[j]->empstte1[i] & (~treenode[j]->empstte0[i]); |
---|
416 | defzero |= treenode[j]->empstte0[i] & (~treenode[j]->empstte1[i]); |
---|
417 | if (added[j]) |
---|
418 | a |= defone & defzero; |
---|
419 | } |
---|
420 | local_suppsteps[i] = defone & defzero & (~a); |
---|
421 | } |
---|
422 | } /* supplement */ |
---|
423 | |
---|
424 | |
---|
425 | void evaluate(node2 *r) |
---|
426 | { |
---|
427 | /* Determines the number of steps needed for a tree. |
---|
428 | This is the minimum number needed to evolve chars on |
---|
429 | this tree */ |
---|
430 | long i, stepnum, smaller; |
---|
431 | double sum; |
---|
432 | |
---|
433 | sum = 0.0; |
---|
434 | for (i = 0; i < (chars); i++) { |
---|
435 | numszero[i] = 0; |
---|
436 | numsone[i] = 0; |
---|
437 | } |
---|
438 | supplement(suppsteps); |
---|
439 | for (i = 0; i < (words); i++) |
---|
440 | zeroanc[i] =fullset; |
---|
441 | full = true; |
---|
442 | postorder(r, fullset, full, wagner, zeroanc); |
---|
443 | cpostorder(r, full, zeroanc, numszero, numsone); |
---|
444 | count(r->fulstte1, zeroanc, numszero, numsone); |
---|
445 | count(suppsteps, zeroanc, numszero, numsone); |
---|
446 | for (i = 0; i < (words); i++) |
---|
447 | zeroanc[i] = 0; |
---|
448 | full = false; |
---|
449 | postorder(r, fullset, full, wagner, zeroanc); |
---|
450 | cpostorder(r, full, zeroanc, numszero, numsone); |
---|
451 | count(r->empstte0, zeroanc, numszero, numsone); |
---|
452 | count(suppsteps, zeroanc, numszero, numsone); |
---|
453 | for (i = 0; i < (chars); i++) { |
---|
454 | smaller = spp * weight[i]; |
---|
455 | numsteps[i] = smaller; |
---|
456 | if (anczero[i]) { |
---|
457 | numsteps[i] = numszero[i]; |
---|
458 | smaller = numszero[i]; |
---|
459 | } |
---|
460 | if (ancone[i] && numsone[i] < smaller) |
---|
461 | numsteps[i] = numsone[i]; |
---|
462 | stepnum = numsteps[i] + extras[i]; |
---|
463 | if (stepnum <= threshwt[i]) |
---|
464 | sum += stepnum; |
---|
465 | else |
---|
466 | sum += threshwt[i]; |
---|
467 | guess[i] = '?'; |
---|
468 | if (!ancone[i] || (anczero[i] && numszero[i] < numsone[i])) |
---|
469 | guess[i] = '0'; |
---|
470 | else if (!anczero[i] || (ancone[i] && numsone[i] < numszero[i])) |
---|
471 | guess[i] = '1'; |
---|
472 | } |
---|
473 | if (examined == 0 && mults == 0) |
---|
474 | bestyet = -1.0; |
---|
475 | like = sum; |
---|
476 | } /* evaluate */ |
---|
477 | |
---|
478 | |
---|
479 | void addtraverse(node2 *a, node2 *b, node2 *c, long *m, long *n, |
---|
480 | valptr valyew, placeptr place) |
---|
481 | { |
---|
482 | /* traverse all places to add b */ |
---|
483 | if (done) |
---|
484 | return; |
---|
485 | if ((*m) <= 2 || !(noroot && (a == root || a == root->next->back))) { |
---|
486 | add3(a, b, c, &root, treenode); |
---|
487 | (*n)++; |
---|
488 | evaluate(root); |
---|
489 | examined++; |
---|
490 | if (examined == howoften) { |
---|
491 | examined = 0; |
---|
492 | mults++; |
---|
493 | if (mults == howmany) |
---|
494 | done = true; |
---|
495 | if (progress) { |
---|
496 | printf("%6ld", mults); |
---|
497 | if (bestyet >= 0) |
---|
498 | printf("%18.5f", bestyet); |
---|
499 | else |
---|
500 | printf(" - "); |
---|
501 | printf("%17ld%20.2f\n", nextree - 1, fracdone * 100); |
---|
502 | #ifdef WIN32 |
---|
503 | phyFillScreenColor(); |
---|
504 | #endif |
---|
505 | } |
---|
506 | } |
---|
507 | valyew[(*n) - 1] = like; |
---|
508 | place[(*n) - 1] = a->index; |
---|
509 | re_move3(&b, &c, &root, treenode); |
---|
510 | } |
---|
511 | if (!a->tip) { |
---|
512 | addtraverse(a->next->back, b, c, m,n,valyew,place); |
---|
513 | addtraverse(a->next->next->back, b, c, m,n,valyew,place); |
---|
514 | } |
---|
515 | } /* addtraverse */ |
---|
516 | |
---|
517 | |
---|
518 | void addit(long m) |
---|
519 | { |
---|
520 | /* adds the species one by one, recursively */ |
---|
521 | long n; |
---|
522 | valptr valyew; |
---|
523 | placeptr place; |
---|
524 | long i, j, n1, besttoadd = 0; |
---|
525 | valptr bestval; |
---|
526 | placeptr bestplace; |
---|
527 | double oldfrac, oldfdone, sum, bestsum; |
---|
528 | |
---|
529 | valyew = (valptr)Malloc(nonodes*sizeof(double)); |
---|
530 | bestval = (valptr)Malloc(nonodes*sizeof(double)); |
---|
531 | place = (placeptr)Malloc(nonodes*sizeof(long)); |
---|
532 | bestplace = (placeptr)Malloc(nonodes*sizeof(long)); |
---|
533 | if (simple && !firsttime) { |
---|
534 | n = 0; |
---|
535 | added[order[m - 1] - 1] = true; |
---|
536 | addtraverse(root, treenode[order[m - 1] - 1], |
---|
537 | treenode[spp + m - 2], &m,&n,valyew,place); |
---|
538 | besttoadd = order[m - 1]; |
---|
539 | memcpy(bestplace, place, nonodes*sizeof(long)); |
---|
540 | memcpy(bestval, valyew, nonodes*sizeof(double)); |
---|
541 | } else { |
---|
542 | bestsum = -1.0; |
---|
543 | for (i = 1; i <= (spp); i++) { |
---|
544 | if (!added[i - 1]) { |
---|
545 | n = 0; |
---|
546 | added[i - 1] = true; |
---|
547 | addtraverse(root, treenode[i - 1], treenode[spp + m - 2], &m, |
---|
548 | &n,valyew,place); |
---|
549 | added[i - 1] = false; |
---|
550 | sum = 0.0; |
---|
551 | for (j = 0; j < (n); j++) |
---|
552 | sum += valyew[j]; |
---|
553 | if (sum > bestsum) { |
---|
554 | bestsum = sum; |
---|
555 | besttoadd = i; |
---|
556 | memcpy(bestplace, place, nonodes*sizeof(long)); |
---|
557 | memcpy(bestval, valyew, nonodes*sizeof(double)); |
---|
558 | } |
---|
559 | } |
---|
560 | } |
---|
561 | } |
---|
562 | order[m - 1] = besttoadd; |
---|
563 | memcpy(place, bestplace, nonodes*sizeof(long)); |
---|
564 | memcpy(valyew, bestval, nonodes*sizeof(double)); |
---|
565 | shellsort(valyew, place, n); |
---|
566 | oldfrac = fracinc; |
---|
567 | oldfdone = fracdone; |
---|
568 | n1 = 0; |
---|
569 | for (i = 0; i < (n); i++) { |
---|
570 | if (valyew[i] <= bestyet || bestyet < 0.0) |
---|
571 | n1++; |
---|
572 | } |
---|
573 | if (n1 > 0) |
---|
574 | fracinc /= n1; |
---|
575 | for (i = 0; i < (n); i++) { |
---|
576 | if (valyew[i] <= bestyet || bestyet < 0.0) { |
---|
577 | current[m - 1] = place[i]; |
---|
578 | add3(treenode[place[i] - 1], treenode[besttoadd - 1], |
---|
579 | treenode[spp + m - 2], &root, treenode); |
---|
580 | added[besttoadd - 1] = true; |
---|
581 | if (m < spp) |
---|
582 | addit(m + 1); |
---|
583 | else { |
---|
584 | if (valyew[i] < bestyet || bestyet < 0.0) { |
---|
585 | nextree = 1; |
---|
586 | bestyet = valyew[i]; |
---|
587 | } |
---|
588 | if (nextree <= maxtrees) { |
---|
589 | memcpy(bestorders[nextree - 1], order, |
---|
590 | spp*sizeof(long)); |
---|
591 | memcpy(bestrees[nextree - 1], current, |
---|
592 | spp*sizeof(long)); |
---|
593 | } |
---|
594 | nextree++; |
---|
595 | firsttime = false; |
---|
596 | } |
---|
597 | re_move3(&treenode[besttoadd - 1], &treenode[spp + m - 2], &root, treenode); |
---|
598 | added[besttoadd - 1] = false; |
---|
599 | } |
---|
600 | fracdone += fracinc; |
---|
601 | } |
---|
602 | fracinc = oldfrac; |
---|
603 | fracdone = oldfdone; |
---|
604 | free(valyew); |
---|
605 | free(bestval); |
---|
606 | free(place); |
---|
607 | free(bestplace); |
---|
608 | } /* addit */ |
---|
609 | |
---|
610 | |
---|
611 | void reroot(node2 *outgroup) |
---|
612 | { |
---|
613 | /* reorients tree, putting outgroup in desired position. */ |
---|
614 | node2 *p, *q, *newbottom, *oldbottom; |
---|
615 | |
---|
616 | if (outgroup->back->index == root->index) |
---|
617 | return; |
---|
618 | newbottom = outgroup->back; |
---|
619 | p = treenode[newbottom->index - 1]->back; |
---|
620 | while (p->index != root->index) { |
---|
621 | oldbottom = treenode[p->index - 1]; |
---|
622 | treenode[p->index - 1] = p; |
---|
623 | p = oldbottom->back; |
---|
624 | } |
---|
625 | p = root->next; |
---|
626 | q = root->next->next; |
---|
627 | p->back->back = q->back; |
---|
628 | q->back->back = p->back; |
---|
629 | p->back = outgroup; |
---|
630 | q->back = outgroup->back; |
---|
631 | outgroup->back->back = root->next->next; |
---|
632 | outgroup->back = root->next; |
---|
633 | treenode[newbottom->index - 1] = newbottom; |
---|
634 | } /* reroot */ |
---|
635 | |
---|
636 | |
---|
637 | void describe() |
---|
638 | { |
---|
639 | /* prints ancestors, steps and table of numbers of steps in |
---|
640 | each character */ |
---|
641 | |
---|
642 | if (stepbox) { |
---|
643 | putc('\n', outfile); |
---|
644 | writesteps(weights, numsteps); |
---|
645 | } |
---|
646 | if (questions && (!noroot || didreroot)) |
---|
647 | guesstates(guess); |
---|
648 | if (ancseq) { |
---|
649 | hypstates(fullset, full, noroot, didreroot, root, wagner, |
---|
650 | zeroanc, oneanc, treenode, guess, garbage); |
---|
651 | putc('\n', outfile); |
---|
652 | } |
---|
653 | if (trout) { |
---|
654 | col = 0; |
---|
655 | treeout2(root, &col, root); |
---|
656 | } |
---|
657 | } /* describe */ |
---|
658 | |
---|
659 | |
---|
660 | void maketree() |
---|
661 | { |
---|
662 | /* tree construction recursively by branch and bound */ |
---|
663 | long i, j, k; |
---|
664 | node2 *dummy; |
---|
665 | |
---|
666 | fullset = (1L << (bits + 1)) - (1L << 1); |
---|
667 | if (progress) { |
---|
668 | printf("\nHow many\n"); |
---|
669 | printf("trees looked Approximate\n"); |
---|
670 | printf("at so far Length of How many percentage\n"); |
---|
671 | printf("(multiples shortest tree trees this long searched\n"); |
---|
672 | printf("of %4ld): found so far found so far so far\n", |
---|
673 | howoften); |
---|
674 | printf("---------- ------------ ------------ ------------\n"); |
---|
675 | #ifdef WIN32 |
---|
676 | phyFillScreenColor(); |
---|
677 | #endif |
---|
678 | } |
---|
679 | done = false; |
---|
680 | mults = 0; |
---|
681 | examined = 0; |
---|
682 | nextree = 1; |
---|
683 | root = treenode[0]; |
---|
684 | firsttime = true; |
---|
685 | for (i = 0; i < (spp); i++) |
---|
686 | added[i] = false; |
---|
687 | added[0] = true; |
---|
688 | order[0] = 1; |
---|
689 | k = 2; |
---|
690 | fracdone = 0.0; |
---|
691 | fracinc = 1.0; |
---|
692 | bestyet = -1.0; |
---|
693 | addit(k); |
---|
694 | if (done) { |
---|
695 | if (progress) { |
---|
696 | printf("Search broken off! Not guaranteed to\n"); |
---|
697 | printf(" have found the most parsimonious trees.\n"); |
---|
698 | } |
---|
699 | if (treeprint) { |
---|
700 | fprintf(outfile, "Search broken off! Not guaranteed to\n"); |
---|
701 | fprintf(outfile, " have found the most parsimonious\n"); |
---|
702 | fprintf(outfile, " trees, but here is what we found:\n"); |
---|
703 | } |
---|
704 | } |
---|
705 | if (treeprint) { |
---|
706 | fprintf(outfile, "\nrequires a total of %18.3f\n\n", bestyet); |
---|
707 | if (nextree == 2) |
---|
708 | fprintf(outfile, "One most parsimonious tree found:\n"); |
---|
709 | else |
---|
710 | fprintf(outfile, "%5ld trees in all found\n", nextree - 1); |
---|
711 | } |
---|
712 | if (nextree > maxtrees + 1) { |
---|
713 | if (treeprint) |
---|
714 | fprintf(outfile, "here are the first%4ld of them\n", (long)maxtrees); |
---|
715 | nextree = maxtrees + 1; |
---|
716 | } |
---|
717 | if (treeprint) |
---|
718 | putc('\n', outfile); |
---|
719 | for (i = 0; i < (spp); i++) |
---|
720 | added[i] = true; |
---|
721 | for (i = 0; i <= (nextree - 2); i++) { |
---|
722 | for (j = k; j <= (spp); j++) |
---|
723 | add3(treenode[bestrees[i][j - 1] - 1], |
---|
724 | treenode[bestorders[i][j - 1] - 1], treenode[spp + j - 2], |
---|
725 | &root, treenode); |
---|
726 | if (noroot) |
---|
727 | reroot(treenode[outgrno - 1]); |
---|
728 | didreroot = (outgropt && noroot); |
---|
729 | evaluate(root); |
---|
730 | printree(treeprint, noroot, didreroot, root); |
---|
731 | describe(); |
---|
732 | for (j = k - 1; j < (spp); j++) |
---|
733 | re_move3(&treenode[bestorders[i][j] - 1], &dummy, &root, treenode); |
---|
734 | } |
---|
735 | if (progress) { |
---|
736 | printf("\nOutput written to file \"%s\"\n\n", outfilename); |
---|
737 | if (trout) |
---|
738 | printf("Trees also written onto file \"%s\"\n\n", outtreename); |
---|
739 | } |
---|
740 | if (ancseq) |
---|
741 | freegarbage(&garbage); |
---|
742 | } /* maketree */ |
---|
743 | |
---|
744 | |
---|
745 | int main(int argc, Char *argv[]) |
---|
746 | { /* Penny's branch-and-bound method */ |
---|
747 | /* Reads in the number of species, number of characters, |
---|
748 | options and data. Then finds all most parsimonious trees */ |
---|
749 | #ifdef MAC |
---|
750 | argc = 1; /* macsetup("Penny",""); */ |
---|
751 | argv[0] = "Penny"; |
---|
752 | #endif |
---|
753 | init(argc,argv); |
---|
754 | openfile(&infile,INFILE,"input file", "r",argv[0],infilename); |
---|
755 | openfile(&outfile,OUTFILE,"output file", "w",argv[0],outfilename); |
---|
756 | ibmpc = IBMCRT; |
---|
757 | ansi = ANSICRT; |
---|
758 | mulsets = false; |
---|
759 | msets = 1; |
---|
760 | firstset = true; |
---|
761 | garbage = NULL; |
---|
762 | bits = 8*sizeof(long) - 1; |
---|
763 | doinit(); |
---|
764 | if (weights || justwts) |
---|
765 | openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename); |
---|
766 | if (trout) |
---|
767 | openfile(&outtree,OUTTREE,"output tree file", "w",argv[0],outtreename); |
---|
768 | if(ancvar) |
---|
769 | openfile(&ancfile,ANCFILE,"ancestors file", "r",argv[0],ancfilename); |
---|
770 | if(mixture) |
---|
771 | openfile(&mixfile,MIXFILE,"mixture file", "r",argv[0],mixfilename); |
---|
772 | |
---|
773 | for (ith = 1; ith <= msets; ith++) { |
---|
774 | if(firstset){ |
---|
775 | if (allsokal && !mixture) |
---|
776 | fprintf(outfile, "Camin-Sokal parsimony method\n\n"); |
---|
777 | if (allwagner && !mixture) |
---|
778 | fprintf(outfile, "Wagner parsimony method\n\n"); |
---|
779 | } |
---|
780 | doinput(); |
---|
781 | if (msets > 1 && !justwts) { |
---|
782 | fprintf(outfile, "Data set # %ld:\n\n",ith); |
---|
783 | if (progress) |
---|
784 | printf("\nData set # %ld:\n",ith); |
---|
785 | } |
---|
786 | if (justwts){ |
---|
787 | if(firstset && mixture && printdata) |
---|
788 | printmixture(outfile, wagner); |
---|
789 | fprintf(outfile, "Weights set # %ld:\n\n", ith); |
---|
790 | if (progress) |
---|
791 | printf("\nWeights set # %ld:\n\n", ith); |
---|
792 | } |
---|
793 | else if (mixture && printdata) |
---|
794 | printmixture(outfile, wagner); |
---|
795 | if (printdata){ |
---|
796 | if (weights || justwts) |
---|
797 | printweights(outfile, 0, chars, weight, "Characters"); |
---|
798 | if (ancvar) |
---|
799 | printancestors(outfile, anczero, ancone); |
---|
800 | } |
---|
801 | if (ith == 1) |
---|
802 | firstset = false; |
---|
803 | maketree(); |
---|
804 | } |
---|
805 | FClose(infile); |
---|
806 | FClose(outfile); |
---|
807 | FClose(outtree); |
---|
808 | #ifdef MAC |
---|
809 | fixmacfile(outfilename); |
---|
810 | fixmacfile(outtreename); |
---|
811 | #endif |
---|
812 | #ifdef WIN32 |
---|
813 | phyRestoreConsoleAttributes(); |
---|
814 | #endif |
---|
815 | return 0; |
---|
816 | } /* Penny's branch-and-bound method */ |
---|
817 | |
---|