1
2 #include "phylip.h"
3 #include "seq.h"
4
5 /* version 3.696.
6 Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
7
8 Copyright (c) 1993-2014, Joseph Felsenstein
9 All rights reserved.
10
11 Redistribution and use in source and binary forms, with or without
12 modification, are permitted provided that the following conditions are met:
13
14 1. Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16
17 2. Redistributions in binary form must reproduce the above copyright notice,
18 this list of conditions and the following disclaimer in the documentation
19 and/or other materials provided with the distribution.
20
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
25 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 long nonodes, endsite, outgrno, nextree, which;
35 boolean interleaved, printdata, outgropt, treeprint, dotdiff, transvp;
36 steptr weight, category, alias, location, ally;
37 sequence y;
38
39
fix_x(node * p,long site,double maxx,long rcategs)40 void fix_x(node* p,long site, double maxx, long rcategs)
41 { /* dnaml dnamlk */
42 long i,j;
43 p->underflows[site] += log(maxx);
44
45 for ( i = 0 ; i < rcategs ; i++ ) {
46 for ( j = 0 ; j < ((long)T - (long)A + 1) ; j++)
47 p->x[site][i][j] /= maxx;
48 }
49 } /* fix_x */
50
51
fix_protx(node * p,long site,double maxx,long rcategs)52 void fix_protx(node* p,long site, double maxx, long rcategs)
53 { /* proml promlk */
54 long i,m;
55
56 p->underflows[site] += log(maxx);
57
58 for ( i = 0 ; i < rcategs ; i++ )
59 for (m = 0; m <= 19; m++)
60 p->protx[site][i][m] /= maxx;
61 } /* fix_protx */
62
63
alloctemp(node ** temp,long * zeros,long endsite)64 void alloctemp(node **temp, long *zeros, long endsite)
65 {
66 /*used in dnacomp and dnapenny */
67 *temp = (node *)Malloc(sizeof(node));
68 (*temp)->numsteps = (steptr)Malloc(endsite*sizeof(long));
69 (*temp)->base = (baseptr)Malloc(endsite*sizeof(long));
70 (*temp)->numnuc = (nucarray *)Malloc(endsite*sizeof(nucarray));
71 memcpy((*temp)->base, zeros, endsite*sizeof(long));
72 memcpy((*temp)->numsteps, zeros, endsite*sizeof(long));
73 zeronumnuc(*temp, endsite);
74 } /* alloctemp */
75
76
freetemp(node ** temp)77 void freetemp(node **temp)
78 {
79 /* used in dnacomp, dnapars, & dnapenny */
80 free((*temp)->numsteps);
81 free((*temp)->base);
82 free((*temp)->numnuc);
83 free(*temp);
84 } /* freetemp */
85
86
freetree2(pointarray treenode,long nonodes)87 void freetree2 (pointarray treenode, long nonodes)
88 {
89 /* The natural complement to alloctree2. Free all elements of all
90 the rings (normally triads) in treenode */
91 long i;
92 node *p, *q;
93
94 /* The first spp elements are just nodes, not rings */
95 for (i = 0; i < spp; i++)
96 free (treenode[i]);
97
98 /* The rest are rings */
99 for (i = spp; i < nonodes; i++) {
100 p = treenode[i]->next;
101 while (p != treenode[i]) {
102 q = p->next;
103 free (p);
104 p = q;
105 }
106 /* p should now point to treenode[i], which has yet to be freed */
107 free (p);
108 }
109 free (treenode);
110 } /* freetree2 */
111
112
113 /*void inputdata(long chars)
114 {
115 * input the names and sequences for each species *
116 * used by dnacomp, dnadist, dnainvar, dnaml, dnamlk, dnapars, & dnapenny *
117 long i, j, k, l, basesread, basesnew=0;
118 Char charstate;
119 boolean allread, done;
120
121 if (printdata)
122 headings(chars, "Sequences", "---------");
123 basesread = 0;
124 allread = false;
125 while (!(allread)) {
126 * eat white space -- if the separator line has spaces on it*
127 do {
128 charstate = gettc(infile);
129 } while (charstate == ' ' || charstate == '\t');
130 ungetc(charstate, infile);
131 if (eoln(infile))
132 scan_eoln(infile);
133 i = 1;
134 while (i <= spp) {
135 if ((interleaved && basesread == 0) || !interleaved)
136 initname(i-1);
137 j = (interleaved) ? basesread : 0;
138 done = false;
139 while (!done && !eoff(infile)) {
140 if (interleaved)
141 done = true;
142 while (j < chars && !(eoln(infile) || eoff(infile))) {
143 charstate = gettc(infile);
144 if (charstate == '\n' || charstate == '\t')
145 charstate = ' ';
146 if (charstate == ' ' || (charstate >= '0' && charstate <= '9'))
147 continue;
148 uppercase(&charstate);
149 if ((strchr("ABCDGHKMNRSTUVWXY?O-",charstate)) == NULL){
150 printf("ERROR: bad base: %c at site %5ld of species %3ld\n",
151 charstate, j+1, i);
152 if (charstate == '.') {
153 printf(" Periods (.) may not be used as gap characters.\n");
154 printf(" The correct gap character is (-)\n");
155 }
156 exxit(-1);
157 }
158 j++;
159 y[i - 1][j - 1] = charstate;
160 }
161 if (interleaved)
162 continue;
163 if (j < chars)
164 scan_eoln(infile);
165 else if (j == chars)
166 done = true;
167 }
168 if (interleaved && i == 1)
169 basesnew = j;
170
171 scan_eoln(infile);
172
173 if ((interleaved && j != basesnew) ||
174 (!interleaved && j != chars)) {
175 printf("\nERROR: sequences out of alignment at position %ld", j+1);
176 printf(" of species %ld\n\n", i);
177 exxit(-1);
178 }
179 i++;
180 }
181
182 if (interleaved) {
183 basesread = basesnew;
184 allread = (basesread == chars);
185 } else
186 allread = (i > spp);
187 }
188 if (!printdata)
189 return;
190 for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
191 for (j = 1; j <= spp; j++) {
192 for (k = 0; k < nmlngth; k++)
193 putc(nayme[j - 1][k], outfile);
194 fprintf(outfile, " ");
195 l = i * 60;
196 if (l > chars)
197 l = chars;
198 for (k = (i - 1) * 60 + 1; k <= l; k++) {
199 if (dotdiff && (j > 1 && y[j - 1][k - 1] == y[0][k - 1]))
200 charstate = '.';
201 else
202 charstate = y[j - 1][k - 1];
203 putc(charstate, outfile);
204 if (k % 10 == 0 && k % 60 != 0)
205 putc(' ', outfile);
206 }
207 putc('\n', outfile);
208 }
209 putc('\n', outfile);
210 }
211 putc('\n', outfile);
212 } * inputdata */
213
214
alloctree(pointarray * treenode,long nonodes,boolean usertree)215 void alloctree(pointarray *treenode, long nonodes, boolean usertree)
216 {
217 /* allocate treenode dynamically */
218 /* used in dnapars, dnacomp, dnapenny & dnamove */
219 long i, j;
220 node *p, *q;
221
222 *treenode = (pointarray)Malloc(nonodes*sizeof(node *));
223 for (i = 0; i < spp; i++) {
224 (*treenode)[i] = (node *)Malloc(sizeof(node));
225 (*treenode)[i]->tip = true;
226 (*treenode)[i]->index = i+1;
227 (*treenode)[i]->iter = true;
228 (*treenode)[i]->branchnum = 0;
229 (*treenode)[i]->initialized = true;
230 }
231 if (!usertree)
232 for (i = spp; i < nonodes; i++) {
233 q = NULL;
234 for (j = 1; j <= 3; j++) {
235 p = (node *)Malloc(sizeof(node));
236 p->tip = false;
237 p->index = i+1;
238 p->iter = true;
239 p->branchnum = 0;
240 p->initialized = false;
241 p->next = q;
242 q = p;
243 }
244 p->next->next->next = p;
245 (*treenode)[i] = p;
246 }
247 } /* alloctree */
248
249
allocx(long nonodes,long rcategs,pointarray treenode,boolean usertree)250 void allocx(long nonodes, long rcategs, pointarray treenode, boolean usertree)
251 {
252 /* allocate x dynamically */
253 /* used in dnaml & dnamlk */
254 long i, j, k;
255 node *p;
256
257 for (i = 0; i < spp; i++){
258 treenode[i]->x = (phenotype)Malloc(endsite*sizeof(ratelike));
259 treenode[i]->underflows = (double *)Malloc(endsite * sizeof (double));
260 for (j = 0; j < endsite; j++)
261 treenode[i]->x[j] = (ratelike)Malloc(rcategs*sizeof(sitelike));
262 }
263 if (!usertree) {
264 for (i = spp; i < nonodes; i++) {
265 p = treenode[i];
266 for (j = 1; j <= 3; j++) {
267 p->underflows = (double *)Malloc (endsite * sizeof (double));
268 p->x = (phenotype)Malloc(endsite*sizeof(ratelike));
269 for (k = 0; k < endsite; k++)
270 p->x[k] = (ratelike)Malloc(rcategs*sizeof(sitelike));
271 p = p->next;
272 }
273 }
274 }
275 } /* allocx */
276
277
prot_allocx(long nonodes,long rcategs,pointarray treenode,boolean usertree)278 void prot_allocx(long nonodes, long rcategs, pointarray treenode,
279 boolean usertree)
280 {
281 /* allocate x dynamically */
282 /* used in proml */
283 long i, j, k;
284 node *p;
285
286 for (i = 0; i < spp; i++){
287 treenode[i]->protx = (pphenotype)Malloc(endsite*sizeof(pratelike));
288 treenode[i]->underflows = (double *)Malloc(endsite*sizeof(double));
289 for (j = 0; j < endsite; j++)
290 treenode[i]->protx[j] = (pratelike)Malloc(rcategs*sizeof(psitelike));
291 }
292 if (!usertree) {
293 for (i = spp; i < nonodes; i++) {
294 p = treenode[i];
295 for (j = 1; j <= 3; j++) {
296 p->protx = (pphenotype)Malloc(endsite*sizeof(pratelike));
297 p->underflows = (double *)Malloc(endsite*sizeof(double));
298 for (k = 0; k < endsite; k++)
299 p->protx[k] = (pratelike)Malloc(rcategs*sizeof(psitelike));
300 p = p->next;
301 }
302 }
303 }
304 } /* prot_allocx */
305
306
307
308
setuptree(pointarray treenode,long nonodes,boolean usertree)309 void setuptree(pointarray treenode, long nonodes, boolean usertree)
310 {
311 /* initialize treenodes */
312 long i;
313 node *p;
314
315 for (i = 1; i <= nonodes; i++) {
316 if (i <= spp || !usertree) {
317 treenode[i-1]->back = NULL;
318 treenode[i-1]->tip = (i <= spp);
319 treenode[i-1]->index = i;
320 treenode[i-1]->numdesc = 0;
321 treenode[i-1]->iter = true;
322 treenode[i-1]->initialized = true;
323 treenode[i-1]->tyme = 0.0;
324 }
325 }
326 if (!usertree) {
327 for (i = spp + 1; i <= nonodes; i++) {
328 p = treenode[i-1]->next;
329 while (p != treenode[i-1]) {
330 p->back = NULL;
331 p->tip = false;
332 p->index = i;
333 p->numdesc = 0;
334 p->iter = true;
335 p->initialized = false;
336 p->tyme = 0.0;
337 p = p->next;
338 }
339 }
340 }
341 } /* setuptree */
342
343
setuptree2(tree * a)344 void setuptree2(tree *a)
345 {
346 /* initialize a tree */
347 /* used in dnaml, dnamlk, & restml */
348
349 a->likelihood = -999999.0;
350 a->start = a->nodep[0]->back;
351 a->root = NULL;
352 } /* setuptree2 */
353
354
alloctip(node * p,long * zeros)355 void alloctip(node *p, long *zeros)
356 { /* allocate a tip node */
357 /* used by dnacomp, dnapars, & dnapenny */
358
359 p->numsteps = (steptr)Malloc(endsite*sizeof(long));
360 p->oldnumsteps = (steptr)Malloc(endsite*sizeof(long));
361 p->base = (baseptr)Malloc(endsite*sizeof(long));
362 p->oldbase = (baseptr)Malloc(endsite*sizeof(long));
363 memcpy(p->base, zeros, endsite*sizeof(long));
364 memcpy(p->numsteps, zeros, endsite*sizeof(long));
365 memcpy(p->oldbase, zeros, endsite*sizeof(long));
366 memcpy(p->oldnumsteps, zeros, endsite*sizeof(long));
367 } /* alloctip */
368
369
370
371
getbasefreqs(double freqa,double freqc,double freqg,double freqt,double * freqr,double * freqy,double * freqar,double * freqcy,double * freqgr,double * freqty,double * ttratio,double * xi,double * xv,double * fracchange,boolean freqsfrom,boolean printdata)372 void getbasefreqs(double freqa, double freqc, double freqg, double freqt,
373 double *freqr, double *freqy, double *freqar, double *freqcy,
374 double *freqgr, double *freqty, double *ttratio, double *xi,
375 double *xv, double *fracchange, boolean freqsfrom,
376 boolean printdata)
377 {
378 /* used by dnadist, dnaml, & dnamlk */
379 double aa, bb;
380
381 if (printdata) {
382 putc('\n', outfile);
383 if (freqsfrom)
384 fprintf(outfile, "Empirical ");
385 fprintf(outfile, "Base Frequencies:\n\n");
386 fprintf(outfile, " A %10.5f\n", freqa);
387 fprintf(outfile, " C %10.5f\n", freqc);
388 fprintf(outfile, " G %10.5f\n", freqg);
389 fprintf(outfile, " T(U) %10.5f\n", freqt);
390 fprintf(outfile, "\n");
391 }
392 *freqr = freqa + freqg;
393 *freqy = freqc + freqt;
394 *freqar = freqa / *freqr;
395 *freqcy = freqc / *freqy;
396 *freqgr = freqg / *freqr;
397 *freqty = freqt / *freqy;
398 aa = *ttratio * (*freqr) * (*freqy) - freqa * freqg - freqc * freqt;
399 bb = freqa * (*freqgr) + freqc * (*freqty);
400 *xi = aa / (aa + bb);
401 *xv = 1.0 - *xi;
402 if (*xi < 0.0) {
403 printf("\n WARNING: This transition/transversion ratio\n");
404 printf(" is impossible with these base frequencies!\n");
405 *xi = 0.0;
406 *xv = 1.0;
407 (*ttratio) = (freqa*freqg+freqc*freqt)/((*freqr)*(*freqy));
408 printf(" Transition/transversion parameter reset\n");
409 printf(" so transition/transversion ratio is %10.6f\n\n", (*ttratio));
410 }
411 if (freqa <= 0.0)
412 freqa = 0.000001;
413 if (freqc <= 0.0)
414 freqc = 0.000001;
415 if (freqg <= 0.0)
416 freqg = 0.000001;
417 if (freqt <= 0.0)
418 freqt = 0.000001;
419 *fracchange = (*xi) * (2 * freqa * (*freqgr) + 2 * freqc * (*freqty)) +
420 (*xv) * (1.0 - freqa * freqa - freqc * freqc - freqg * freqg
421 - freqt * freqt);
422 } /* getbasefreqs */
423
424
empiricalfreqs(double * freqa,double * freqc,double * freqg,double * freqt,steptr weight,pointarray treenode)425 void empiricalfreqs(double *freqa, double *freqc, double *freqg,
426 double *freqt, steptr weight, pointarray treenode)
427 {
428 /* Get empirical base frequencies from the data */
429 /* used in dnaml & dnamlk */
430 long i, j, k;
431 double sum, suma, sumc, sumg, sumt, w;
432
433 *freqa = 0.25;
434 *freqc = 0.25;
435 *freqg = 0.25;
436 *freqt = 0.25;
437 for (k = 1; k <= 8; k++) {
438 suma = 0.0;
439 sumc = 0.0;
440 sumg = 0.0;
441 sumt = 0.0;
442 for (i = 0; i < spp; i++) {
443 for (j = 0; j < endsite; j++) {
444 w = weight[j];
445 sum = (*freqa) * treenode[i]->x[j][0][0];
446 sum += (*freqc) * treenode[i]->x[j][0][(long)C - (long)A];
447 sum += (*freqg) * treenode[i]->x[j][0][(long)G - (long)A];
448 sum += (*freqt) * treenode[i]->x[j][0][(long)T - (long)A];
449 suma += w * (*freqa) * treenode[i]->x[j][0][0] / sum;
450 sumc += w * (*freqc) * treenode[i]->x[j][0][(long)C - (long)A] / sum;
451 sumg += w * (*freqg) * treenode[i]->x[j][0][(long)G - (long)A] / sum;
452 sumt += w * (*freqt) * treenode[i]->x[j][0][(long)T - (long)A] / sum;
453 }
454 }
455 sum = suma + sumc + sumg + sumt;
456 *freqa = suma / sum;
457 *freqc = sumc / sum;
458 *freqg = sumg / sum;
459 *freqt = sumt / sum;
460 }
461 if (*freqa <= 0.0)
462 *freqa = 0.000001;
463 if (*freqc <= 0.0)
464 *freqc = 0.000001;
465 if (*freqg <= 0.0)
466 *freqg = 0.000001;
467 if (*freqt <= 0.0)
468 *freqt = 0.000001;
469 } /* empiricalfreqs */
470
471
sitesort(long chars,steptr weight)472 void sitesort(long chars, steptr weight)
473 {
474 /* Shell sort keeping sites, weights in same order */
475 /* used in dnainvar, dnapars, dnacomp & dnapenny */
476 long gap, i, j, jj, jg, k, itemp;
477 boolean flip, tied;
478
479 gap = chars / 2;
480 while (gap > 0) {
481 for (i = gap + 1; i <= chars; i++) {
482 j = i - gap;
483 flip = true;
484 while (j > 0 && flip) {
485 jj = alias[j - 1];
486 jg = alias[j + gap - 1];
487 tied = true;
488 k = 1;
489 while (k <= spp && tied) {
490 flip = (y[k - 1][jj - 1] > y[k - 1][jg - 1]);
491 tied = (tied && y[k - 1][jj - 1] == y[k - 1][jg - 1]);
492 k++;
493 }
494 if (!flip)
495 break;
496 itemp = alias[j - 1];
497 alias[j - 1] = alias[j + gap - 1];
498 alias[j + gap - 1] = itemp;
499 itemp = weight[j - 1];
500 weight[j - 1] = weight[j + gap - 1];
501 weight[j + gap - 1] = itemp;
502 j -= gap;
503 }
504 }
505 gap /= 2;
506 }
507 } /* sitesort */
508
509
sitecombine(long chars)510 void sitecombine(long chars)
511 {
512 /* combine sites that have identical patterns */
513 /* used in dnapars, dnapenny, & dnacomp */
514 long i, j, k;
515 boolean tied;
516
517 i = 1;
518 while (i < chars) {
519 j = i + 1;
520 tied = true;
521 while (j <= chars && tied) {
522 k = 1;
523 while (k <= spp && tied) {
524 tied = (tied &&
525 y[k - 1][alias[i - 1] - 1] == y[k - 1][alias[j - 1] - 1]);
526 k++;
527 }
528 if (tied) {
529 weight[i - 1] += weight[j - 1];
530 weight[j - 1] = 0;
531 ally[alias[j - 1] - 1] = alias[i - 1];
532 }
533 j++;
534 }
535 i = j - 1;
536 }
537 } /* sitecombine */
538
539
sitescrunch(long chars)540 void sitescrunch(long chars)
541 {
542 /* move so one representative of each pattern of
543 sites comes first */
544 /* used in dnapars & dnacomp */
545 long i, j, itemp;
546 boolean done, found;
547
548 done = false;
549 i = 1;
550 j = 2;
551 while (!done) {
552 if (ally[alias[i - 1] - 1] != alias[i - 1]) {
553 if (j <= i)
554 j = i + 1;
555 if (j <= chars) {
556 do {
557 found = (ally[alias[j - 1] - 1] == alias[j - 1]);
558 j++;
559 } while (!(found || j > chars));
560 if (found) {
561 j--;
562 itemp = alias[i - 1];
563 alias[i - 1] = alias[j - 1];
564 alias[j - 1] = itemp;
565 itemp = weight[i - 1];
566 weight[i - 1] = weight[j - 1];
567 weight[j - 1] = itemp;
568 } else
569 done = true;
570 } else
571 done = true;
572 }
573 i++;
574 done = (done || i >= chars);
575 }
576 } /* sitescrunch */
577
578
sitesort2(long sites,steptr aliasweight)579 void sitesort2(long sites, steptr aliasweight)
580 {
581 /* Shell sort keeping sites, weights in same order */
582 /* used in dnaml & dnamnlk */
583 long gap, i, j, jj, jg, k, itemp;
584 boolean flip, tied;
585
586 gap = sites / 2;
587 while (gap > 0) {
588 for (i = gap + 1; i <= sites; i++) {
589 j = i - gap;
590 flip = true;
591 while (j > 0 && flip) {
592 jj = alias[j - 1];
593 jg = alias[j + gap - 1];
594 tied = (category[jj - 1] == category[jg - 1]);
595 flip = (category[jj - 1] > category[jg - 1]);
596 k = 1;
597 while (k <= spp && tied) {
598 flip = (y[k - 1][jj - 1] > y[k - 1][jg - 1]);
599 tied = (tied && y[k - 1][jj - 1] == y[k - 1][jg - 1]);
600 k++;
601 }
602 if (!flip)
603 break;
604 itemp = alias[j - 1];
605 alias[j - 1] = alias[j + gap - 1];
606 alias[j + gap - 1] = itemp;
607 itemp = aliasweight[j - 1];
608 aliasweight[j - 1] = aliasweight[j + gap - 1];
609 aliasweight[j + gap - 1] = itemp;
610 j -= gap;
611 }
612 }
613 gap /= 2;
614 }
615 } /* sitesort2 */
616
617
sitecombine2(long sites,steptr aliasweight)618 void sitecombine2(long sites, steptr aliasweight)
619 {
620 /* combine sites that have identical patterns */
621 /* used in dnaml & dnamlk */
622 long i, j, k;
623 boolean tied;
624
625 i = 1;
626 while (i < sites) {
627 j = i + 1;
628 tied = true;
629 while (j <= sites && tied) {
630 tied = (category[alias[i - 1] - 1] == category[alias[j - 1] - 1]);
631 k = 1;
632 while (k <= spp && tied) {
633 tied = (tied &&
634 y[k - 1][alias[i - 1] - 1] == y[k - 1][alias[j - 1] - 1]);
635 k++;
636 }
637 if (!tied)
638 break;
639 aliasweight[i - 1] += aliasweight[j - 1];
640 aliasweight[j - 1] = 0;
641 ally[alias[j - 1] - 1] = alias[i - 1];
642 j++;
643 }
644 i = j;
645 }
646 } /* sitecombine2 */
647
648
sitescrunch2(long sites,long i,long j,steptr aliasweight)649 void sitescrunch2(long sites, long i, long j, steptr aliasweight)
650 {
651 /* move so positively weighted sites come first */
652 /* used by dnainvar, dnaml, dnamlk, & restml */
653 long itemp;
654 boolean done, found;
655
656 done = false;
657 while (!done) {
658 if (aliasweight[i - 1] > 0)
659 i++;
660 else {
661 if (j <= i)
662 j = i + 1;
663 if (j <= sites) {
664 do {
665 found = (aliasweight[j - 1] > 0);
666 j++;
667 } while (!(found || j > sites));
668 if (found) {
669 j--;
670 itemp = alias[i - 1];
671 alias[i - 1] = alias[j - 1];
672 alias[j - 1] = itemp;
673 itemp = aliasweight[i - 1];
674 aliasweight[i - 1] = aliasweight[j - 1];
675 aliasweight[j - 1] = itemp;
676 } else
677 done = true;
678 } else
679 done = true;
680 }
681 done = (done || i >= sites);
682 }
683 } /* sitescrunch2 */
684
685
makevalues(pointarray treenode,long * zeros,boolean usertree)686 void makevalues(pointarray treenode, long *zeros, boolean usertree)
687 {
688 /* set up fractional likelihoods at tips */
689 /* used by dnacomp, dnapars, & dnapenny */
690 long i, j;
691 char ns = 0;
692 node *p;
693
694 setuptree(treenode, nonodes, usertree);
695 for (i = 0; i < spp; i++)
696 alloctip(treenode[i], zeros);
697 if (!usertree) {
698 for (i = spp; i < nonodes; i++) {
699 p = treenode[i];
700 do {
701 allocnontip(p, zeros, endsite);
702 p = p->next;
703 } while (p != treenode[i]);
704 }
705 }
706 for (j = 0; j < endsite; j++) {
707 for (i = 0; i < spp; i++) {
708 switch (y[i][alias[j] - 1]) {
709
710 case 'A':
711 ns = 1 << A;
712 break;
713
714 case 'C':
715 ns = 1 << C;
716 break;
717
718 case 'G':
719 ns = 1 << G;
720 break;
721
722 case 'U':
723 ns = 1 << T;
724 break;
725
726 case 'T':
727 ns = 1 << T;
728 break;
729
730 case 'M':
731 ns = (1 << A) | (1 << C);
732 break;
733
734 case 'R':
735 ns = (1 << A) | (1 << G);
736 break;
737
738 case 'W':
739 ns = (1 << A) | (1 << T);
740 break;
741
742 case 'S':
743 ns = (1 << C) | (1 << G);
744 break;
745
746 case 'Y':
747 ns = (1 << C) | (1 << T);
748 break;
749
750 case 'K':
751 ns = (1 << G) | (1 << T);
752 break;
753
754 case 'B':
755 ns = (1 << C) | (1 << G) | (1 << T);
756 break;
757
758 case 'D':
759 ns = (1 << A) | (1 << G) | (1 << T);
760 break;
761
762 case 'H':
763 ns = (1 << A) | (1 << C) | (1 << T);
764 break;
765
766 case 'V':
767 ns = (1 << A) | (1 << C) | (1 << G);
768 break;
769
770 case 'N':
771 ns = (1 << A) | (1 << C) | (1 << G) | (1 << T);
772 break;
773
774 case 'X':
775 ns = (1 << A) | (1 << C) | (1 << G) | (1 << T);
776 break;
777
778 case '?':
779 ns = (1 << A) | (1 << C) | (1 << G) | (1 << T) | (1 << O);
780 break;
781
782 case 'O':
783 ns = 1 << O;
784 break;
785
786 case '-':
787 ns = 1 << O;
788 break;
789 }
790 treenode[i]->base[j] = ns;
791 treenode[i]->numsteps[j] = 0;
792 }
793 }
794 } /* makevalues */
795
796
makevalues2(long categs,pointarray treenode,long endsite,long spp,sequence y,steptr alias)797 void makevalues2(long categs, pointarray treenode, long endsite,
798 long spp, sequence y, steptr alias)
799 {
800 /* set up fractional likelihoods at tips */
801 /* used by dnaml & dnamlk */
802 long i, j, k, l;
803 bases b;
804
805 for (k = 0; k < endsite; k++) {
806 j = alias[k];
807 for (i = 0; i < spp; i++) {
808 for (l = 0; l < categs; l++) {
809 for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
810 treenode[i]->x[k][l][(long)b - (long)A] = 0.0;
811 switch (y[i][j - 1]) {
812
813 case 'A':
814 treenode[i]->x[k][l][0] = 1.0;
815 break;
816
817 case 'C':
818 treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
819 break;
820
821 case 'G':
822 treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
823 break;
824
825 case 'T':
826 treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
827 break;
828
829 case 'U':
830 treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
831 break;
832
833 case 'M':
834 treenode[i]->x[k][l][0] = 1.0;
835 treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
836 break;
837
838 case 'R':
839 treenode[i]->x[k][l][0] = 1.0;
840 treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
841 break;
842
843 case 'W':
844 treenode[i]->x[k][l][0] = 1.0;
845 treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
846 break;
847
848 case 'S':
849 treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
850 treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
851 break;
852
853 case 'Y':
854 treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
855 treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
856 break;
857
858 case 'K':
859 treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
860 treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
861 break;
862
863 case 'B':
864 treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
865 treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
866 treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
867 break;
868
869 case 'D':
870 treenode[i]->x[k][l][0] = 1.0;
871 treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
872 treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
873 break;
874
875 case 'H':
876 treenode[i]->x[k][l][0] = 1.0;
877 treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
878 treenode[i]->x[k][l][(long)T - (long)A] = 1.0;
879 break;
880
881 case 'V':
882 treenode[i]->x[k][l][0] = 1.0;
883 treenode[i]->x[k][l][(long)C - (long)A] = 1.0;
884 treenode[i]->x[k][l][(long)G - (long)A] = 1.0;
885 break;
886
887 case 'N':
888 for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
889 treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
890 break;
891
892 case 'X':
893 for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
894 treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
895 break;
896
897 case '?':
898 for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
899 treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
900 break;
901
902 case 'O':
903 for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
904 treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
905 break;
906
907 case '-':
908 for (b = A; (long)b <= (long)T; b = (bases)((long)b + 1))
909 treenode[i]->x[k][l][(long)b - (long)A] = 1.0;
910 break;
911 }
912 }
913 }
914 }
915 } /* makevalues2 */
916
917
fillin(node * p,node * left,node * rt)918 void fillin(node *p, node *left, node *rt)
919 {
920 /* sets up for each node in the tree the base sequence
921 at that point and counts the changes. */
922 long i, j, k, n, purset, pyrset;
923 node *q;
924
925 purset = (1 << (long)A) + (1 << (long)G);
926 pyrset = (1 << (long)C) + (1 << (long)T);
927 if (!left) {
928 memcpy(p->base, rt->base, endsite*sizeof(long));
929 memcpy(p->numsteps, rt->numsteps, endsite*sizeof(long));
930 q = rt;
931 } else if (!rt) {
932 memcpy(p->base, left->base, endsite*sizeof(long));
933 memcpy(p->numsteps, left->numsteps, endsite*sizeof(long));
934 q = left;
935 } else {
936 for (i = 0; i < endsite; i++) {
937 p->base[i] = left->base[i] & rt->base[i];
938 p->numsteps[i] = left->numsteps[i] + rt->numsteps[i];
939 if (p->base[i] == 0) {
940 p->base[i] = left->base[i] | rt->base[i];
941 if (transvp) {
942 if (!((p->base[i] == purset) || (p->base[i] == pyrset)))
943 p->numsteps[i] += weight[i];
944 }
945 else p->numsteps[i] += weight[i];
946 }
947 }
948 q = rt;
949 }
950 if (left && rt) n = 2;
951 else n = 1;
952 for (i = 0; i < endsite; i++)
953 for (j = (long)A; j <= (long)O; j++)
954 p->numnuc[i][j] = 0;
955 for (k = 1; k <= n; k++) {
956 if (k == 2) q = left;
957 for (i = 0; i < endsite; i++) {
958 for (j = (long)A; j <= (long)O; j++) {
959 if (q->base[i] & (1 << j))
960 p->numnuc[i][j]++;
961 }
962 }
963 }
964 } /* fillin */
965
966
getlargest(long * numnuc)967 long getlargest(long *numnuc)
968 {
969 /* find the largest in array numnuc */
970 long i, largest;
971
972 largest = 0;
973 for (i = (long)A; i <= (long)O; i++)
974 if (numnuc[i] > largest)
975 largest = numnuc[i];
976 return largest;
977 } /* getlargest */
978
979
multifillin(node * p,node * q,long dnumdesc)980 void multifillin(node *p, node *q, long dnumdesc)
981 {
982 /* sets up for each node in the tree the base sequence
983 at that point and counts the changes according to the
984 changes in q's base */
985 long i, j, b, largest, descsteps, purset, pyrset;
986
987 memcpy(p->oldbase, p->base, endsite*sizeof(long));
988 memcpy(p->oldnumsteps, p->numsteps, endsite*sizeof(long));
989 purset = (1 << (long)A) + (1 << (long)G);
990 pyrset = (1 << (long)C) + (1 << (long)T);
991 for (i = 0; i < endsite; i++) {
992 descsteps = 0;
993 for (j = (long)A; j <= (long)O; j++) {
994 b = 1 << j;
995 if ((descsteps == 0) && (p->base[i] & b))
996 descsteps = p->numsteps[i]
997 - (p->numdesc - dnumdesc - p->numnuc[i][j]) * weight[i];
998 }
999 if (dnumdesc == -1)
1000 descsteps -= q->oldnumsteps[i];
1001 else if (dnumdesc == 0)
1002 descsteps += (q->numsteps[i] - q->oldnumsteps[i]);
1003 else
1004 descsteps += q->numsteps[i];
1005 if (q->oldbase[i] != q->base[i]) {
1006 for (j = (long)A; j <= (long)O; j++) {
1007 b = 1 << j;
1008 if (transvp) {
1009 if (b & purset) b = purset;
1010 if (b & pyrset) b = pyrset;
1011 }
1012 if ((q->oldbase[i] & b) && !(q->base[i] & b))
1013 p->numnuc[i][j]--;
1014 else if (!(q->oldbase[i] & b) && (q->base[i] & b))
1015 p->numnuc[i][j]++;
1016 }
1017 }
1018 largest = getlargest(p->numnuc[i]);
1019 if (q->oldbase[i] != q->base[i]) {
1020 p->base[i] = 0;
1021 for (j = (long)A; j <= (long)O; j++) {
1022 if (p->numnuc[i][j] == largest)
1023 p->base[i] |= (1 << j);
1024 }
1025 }
1026 p->numsteps[i] = (p->numdesc - largest) * weight[i] + descsteps;
1027 }
1028 } /* multifillin */
1029
1030
sumnsteps(node * p,node * left,node * rt,long a,long b)1031 void sumnsteps(node *p, node *left, node *rt, long a, long b)
1032 {
1033 /* sets up for each node in the tree the base sequence
1034 at that point and counts the changes. */
1035 long i;
1036 long ns, rs, ls, purset, pyrset;
1037
1038 if (!left) {
1039 memcpy(p->numsteps, rt->numsteps, endsite*sizeof(long));
1040 memcpy(p->base, rt->base, endsite*sizeof(long));
1041 } else if (!rt) {
1042 memcpy(p->numsteps, left->numsteps, endsite*sizeof(long));
1043 memcpy(p->base, left->base, endsite*sizeof(long));
1044 } else {
1045 purset = (1 << (long)A) + (1 << (long)G);
1046 pyrset = (1 << (long)C) + (1 << (long)T);
1047 for (i = a; i < b; i++) {
1048 ls = left->base[i];
1049 rs = rt->base[i];
1050 ns = ls & rs;
1051 p->numsteps[i] = left->numsteps[i] + rt->numsteps[i];
1052 if (ns == 0) {
1053 ns = ls | rs;
1054 if (transvp) {
1055 if (!((ns == purset) || (ns == pyrset)))
1056 p->numsteps[i] += weight[i];
1057 }
1058 else p->numsteps[i] += weight[i];
1059 }
1060 p->base[i] = ns;
1061 }
1062 }
1063 } /* sumnsteps */
1064
1065
sumnsteps2(node * p,node * left,node * rt,long a,long b,long * threshwt)1066 void sumnsteps2(node *p,node *left,node *rt,long a,long b,long *threshwt)
1067 {
1068 /* counts the changes at each node. */
1069 long i, steps;
1070 long ns, rs, ls, purset, pyrset;
1071 long term;
1072
1073 if (a == 0) p->sumsteps = 0.0;
1074 if (!left)
1075 memcpy(p->numsteps, rt->numsteps, endsite*sizeof(long));
1076 else if (!rt)
1077 memcpy(p->numsteps, left->numsteps, endsite*sizeof(long));
1078 else {
1079 purset = (1 << (long)A) + (1 << (long)G);
1080 pyrset = (1 << (long)C) + (1 << (long)T);
1081 for (i = a; i < b; i++) {
1082 ls = left->base[i];
1083 rs = rt->base[i];
1084 ns = ls & rs;
1085 p->numsteps[i] = left->numsteps[i] + rt->numsteps[i];
1086 if (ns == 0) {
1087 ns = ls | rs;
1088 if (transvp) {
1089 if (!((ns == purset) || (ns == pyrset)))
1090 p->numsteps[i] += weight[i];
1091 }
1092 else p->numsteps[i] += weight[i];
1093 }
1094 }
1095 }
1096 for (i = a; i < b; i++) {
1097 steps = p->numsteps[i];
1098 if ((long)steps <= threshwt[i])
1099 term = steps;
1100 else
1101 term = threshwt[i];
1102 p->sumsteps += (double)term;
1103 }
1104 } /* sumnsteps2 */
1105
1106
multisumnsteps(node * p,node * q,long a,long b,long * threshwt)1107 void multisumnsteps(node *p, node *q, long a, long b, long *threshwt)
1108 {
1109 /* computes the number of steps between p and q */
1110 long i, j, steps, largest, descsteps, purset, pyrset, b1;
1111 long term;
1112
1113 if (a == 0) p->sumsteps = 0.0;
1114 purset = (1 << (long)A) + (1 << (long)G);
1115 pyrset = (1 << (long)C) + (1 << (long)T);
1116 for (i = a; i < b; i++) {
1117 descsteps = 0;
1118 for (j = (long)A; j <= (long)O; j++) {
1119 if ((descsteps == 0) && (p->base[i] & (1 << j)))
1120 descsteps = p->numsteps[i] -
1121 (p->numdesc - 1 - p->numnuc[i][j]) * weight[i];
1122 }
1123 descsteps += q->numsteps[i];
1124 largest = 0;
1125 for (j = (long)A; j <= (long)O; j++) {
1126 b1 = (1 << j);
1127 if (transvp) {
1128 if (b1 & purset) b1 = purset;
1129 if (b1 & pyrset) b1 = pyrset;
1130 }
1131 if (q->base[i] & b1)
1132 p->numnuc[i][j]++;
1133 if (p->numnuc[i][j] > largest)
1134 largest = p->numnuc[i][j];
1135 }
1136 steps = (p->numdesc - largest) * weight[i] + descsteps;
1137 if ((long)steps <= threshwt[i])
1138 term = steps;
1139 else
1140 term = threshwt[i];
1141 p->sumsteps += (double)term;
1142 }
1143 } /* multisumnsteps */
1144
1145
multisumnsteps2(node * p)1146 void multisumnsteps2(node *p)
1147 {
1148 /* counts the changes at each multi-way node. Sums up
1149 steps of all descendants */
1150 long i, j, largest, purset, pyrset, b1;
1151 node *q;
1152 baseptr b;
1153
1154 purset = (1 << (long)A) + (1 << (long)G);
1155 pyrset = (1 << (long)C) + (1 << (long)T);
1156 for (i = 0; i < endsite; i++) {
1157 p->numsteps[i] = 0;
1158 q = p->next;
1159 while (q != p) {
1160 if (q->back) {
1161 p->numsteps[i] += q->back->numsteps[i];
1162 b = q->back->base;
1163 for (j = (long)A; j <= (long)O; j++) {
1164 b1 = (1 << j);
1165 if (transvp) {
1166 if (b1 & purset) b1 = purset;
1167 if (b1 & pyrset) b1 = pyrset;
1168 }
1169 if (b[i] & b1) p->numnuc[i][j]++;
1170 }
1171 }
1172 q = q->next;
1173 }
1174 largest = getlargest(p->numnuc[i]);
1175 p->base[i] = 0;
1176 for (j = (long)A; j <= (long)O; j++) {
1177 if (p->numnuc[i][j] == largest)
1178 p->base[i] |= (1 << j);
1179 }
1180 p->numsteps[i] += ((p->numdesc - largest) * weight[i]);
1181 }
1182 } /* multisumnsteps2 */
1183
alltips(node * forknode,node * p)1184 boolean alltips(node *forknode, node *p)
1185 {
1186 /* returns true if all descendants of forknode except p are tips;
1187 false otherwise. */
1188 node *q, *r;
1189 boolean tips;
1190
1191 tips = true;
1192 r = forknode;
1193 q = forknode->next;
1194 do {
1195 if (q->back && q->back != p && !q->back->tip)
1196 tips = false;
1197 q = q->next;
1198 } while (tips && q != r);
1199 return tips;
1200 } /* alltips */
1201
1202
gdispose(node * p,node ** grbg,pointarray treenode)1203 void gdispose(node *p, node **grbg, pointarray treenode)
1204 {
1205 /* go through tree throwing away nodes */
1206 node *q, *r;
1207
1208 p->back = NULL;
1209 if (p->tip)
1210 return;
1211 treenode[p->index - 1] = NULL;
1212 q = p->next;
1213 while (q != p) {
1214 gdispose(q->back, grbg, treenode);
1215 q->back = NULL;
1216 r = q;
1217 q = q->next;
1218 chuck(grbg, r);
1219 }
1220 chuck(grbg, q);
1221 } /* gdispose */
1222
1223
preorder(node * p,node * r,node * root,node * removing,node * adding,node * changing,long dnumdesc)1224 void preorder(node *p, node *r, node *root, node *removing, node *adding,
1225 node *changing, long dnumdesc)
1226 {
1227 /* recompute number of steps in preorder taking both ancestoral and
1228 descendent steps into account. removing points to a node being
1229 removed, if any */
1230 node *q, *p1, *p2;
1231
1232 if (p && !p->tip && p != adding) {
1233 q = p;
1234 do {
1235 if (p->back != r) {
1236 if (p->numdesc > 2) {
1237 if (changing)
1238 multifillin (p, r, dnumdesc);
1239 else
1240 multifillin (p, r, 0);
1241 } else {
1242 p1 = p->next;
1243 if (!removing)
1244 while (!p1->back)
1245 p1 = p1->next;
1246 else
1247 while (!p1->back || p1->back == removing)
1248 p1 = p1->next;
1249 p2 = p1->next;
1250 if (!removing)
1251 while (!p2->back)
1252 p2 = p2->next;
1253 else
1254 while (!p2->back || p2->back == removing)
1255 p2 = p2->next;
1256 p1 = p1->back;
1257 p2 = p2->back;
1258 if (p->back == p1) p1 = NULL;
1259 else if (p->back == p2) p2 = NULL;
1260 memcpy(p->oldbase, p->base, endsite*sizeof(long));
1261 memcpy(p->oldnumsteps, p->numsteps, endsite*sizeof(long));
1262 fillin(p, p1, p2);
1263 }
1264 }
1265 p = p->next;
1266 } while (p != q);
1267 q = p;
1268 do {
1269 preorder(p->next->back, p->next, root, removing, adding, NULL, 0);
1270 p = p->next;
1271 } while (p->next != q);
1272 }
1273 } /* preorder */
1274
1275
updatenumdesc(node * p,node * root,long n)1276 void updatenumdesc(node *p, node *root, long n)
1277 {
1278 /* set p's numdesc to n. If p is the root, numdesc of p's
1279 descendants are set to n-1. */
1280 node *q;
1281
1282 q = p;
1283 if (p == root && n > 0) {
1284 p->numdesc = n;
1285 n--;
1286 q = q->next;
1287 }
1288 do {
1289 q->numdesc = n;
1290 q = q->next;
1291 } while (q != p);
1292 } /* updatenumdesc */
1293
1294
add(node * below,node * newtip,node * newfork,node ** root,boolean recompute,pointarray treenode,node ** grbg,long * zeros)1295 void add(node *below,node *newtip,node *newfork,node **root,
1296 boolean recompute,pointarray treenode,node **grbg,long *zeros)
1297 {
1298 /* inserts the nodes newfork and its left descendant, newtip,
1299 to the tree. below becomes newfork's right descendant.
1300 if newfork is NULL, newtip is added as below's sibling */
1301 /* used in dnacomp & dnapars */
1302 node *p;
1303
1304 if (below != treenode[below->index - 1])
1305 below = treenode[below->index - 1];
1306 if (newfork) {
1307 if (below->back != NULL)
1308 below->back->back = newfork;
1309 newfork->back = below->back;
1310 below->back = newfork->next->next;
1311 newfork->next->next->back = below;
1312 newfork->next->back = newtip;
1313 newtip->back = newfork->next;
1314 if (*root == below)
1315 *root = newfork;
1316 updatenumdesc(newfork, *root, 2);
1317 } else {
1318 gnutreenode(grbg, &p, below->index, endsite, zeros);
1319 p->back = newtip;
1320 newtip->back = p;
1321 p->next = below->next;
1322 below->next = p;
1323 updatenumdesc(below, *root, below->numdesc + 1);
1324 }
1325 if (!newtip->tip)
1326 updatenumdesc(newtip, *root, newtip->numdesc);
1327 (*root)->back = NULL;
1328 if (!recompute)
1329 return;
1330 if (!newfork) {
1331 memcpy(newtip->back->base, below->base, endsite*sizeof(long));
1332 memcpy(newtip->back->numsteps, below->numsteps, endsite*sizeof(long));
1333 memcpy(newtip->back->numnuc, below->numnuc, endsite*sizeof(nucarray));
1334 if (below != *root) {
1335 memcpy(below->back->oldbase, zeros, endsite*sizeof(long));
1336 memcpy(below->back->oldnumsteps, zeros, endsite*sizeof(long));
1337 multifillin(newtip->back, below->back, 1);
1338 }
1339 if (!newtip->tip) {
1340 memcpy(newtip->back->oldbase, zeros, endsite*sizeof(long));
1341 memcpy(newtip->back->oldnumsteps, zeros, endsite*sizeof(long));
1342 preorder(newtip, newtip->back, *root, NULL, NULL, below, 1);
1343 }
1344 memcpy(newtip->oldbase, zeros, endsite*sizeof(long));
1345 memcpy(newtip->oldnumsteps, zeros, endsite*sizeof(long));
1346 preorder(below, newtip, *root, NULL, newtip, below, 1);
1347 if (below != *root)
1348 preorder(below->back, below, *root, NULL, NULL, NULL, 0);
1349 } else {
1350 fillin(newtip->back, newtip->back->next->back,
1351 newtip->back->next->next->back);
1352 if (!newtip->tip) {
1353 memcpy(newtip->back->oldbase, zeros, endsite*sizeof(long));
1354 memcpy(newtip->back->oldnumsteps, zeros, endsite*sizeof(long));
1355 preorder(newtip, newtip->back, *root, NULL, NULL, newfork, 1);
1356 }
1357 if (newfork != *root) {
1358 memcpy(below->back->base, newfork->back->base, endsite*sizeof(long));
1359 memcpy(below->back->numsteps, newfork->back->numsteps, endsite*sizeof(long));
1360 preorder(newfork, newtip, *root, NULL, newtip, NULL, 0);
1361 } else {
1362 fillin(below->back, newtip, NULL);
1363 fillin(newfork, newtip, below);
1364 memcpy(below->back->oldbase, zeros, endsite*sizeof(long));
1365 memcpy(below->back->oldnumsteps, zeros, endsite*sizeof(long));
1366 preorder(below, below->back, *root, NULL, NULL, newfork, 1);
1367 }
1368 if (newfork != *root) {
1369 memcpy(newfork->oldbase, below->base, endsite*sizeof(long));
1370 memcpy(newfork->oldnumsteps, below->numsteps, endsite*sizeof(long));
1371 preorder(newfork->back, newfork, *root, NULL, NULL, NULL, 0);
1372 }
1373 }
1374 } /* add */
1375
1376
findbelow(node ** below,node * item,node * fork)1377 void findbelow(node **below, node *item, node *fork)
1378 {
1379 /* decide which of fork's binary children is below */
1380
1381 if (fork->next->back == item)
1382 *below = fork->next->next->back;
1383 else
1384 *below = fork->next->back;
1385 } /* findbelow */
1386
1387
re_move(node * item,node ** fork,node ** root,boolean recompute,pointarray treenode,node ** grbg,long * zeros)1388 void re_move(node *item, node **fork, node **root, boolean recompute,
1389 pointarray treenode, node **grbg, long *zeros)
1390 {
1391 /* removes nodes item and its ancestor, fork, from the tree.
1392 the new descendant of fork's ancestor is made to be
1393 fork's second descendant (other than item). Also
1394 returns pointers to the deleted nodes, item and fork.
1395 If item belongs to a node with more than 2 descendants,
1396 fork will not be deleted */
1397 /* used in dnacomp & dnapars */
1398 node *p, *q, *other = NULL, *otherback = NULL;
1399
1400 if (item->back == NULL) {
1401 *fork = NULL;
1402 return;
1403 }
1404 *fork = treenode[item->back->index - 1];
1405 if ((*fork)->numdesc == 2) {
1406 updatenumdesc(*fork, *root, 0);
1407 findbelow(&other, item, *fork);
1408 otherback = other->back;
1409 if (*root == *fork) {
1410 *root = other;
1411 if (!other->tip)
1412 updatenumdesc(other, *root, other->numdesc);
1413 }
1414 p = item->back->next->back;
1415 q = item->back->next->next->back;
1416 if (p != NULL)
1417 p->back = q;
1418 if (q != NULL)
1419 q->back = p;
1420 (*fork)->back = NULL;
1421 p = (*fork)->next;
1422 while (p != *fork) {
1423 p->back = NULL;
1424 p = p->next;
1425 }
1426 } else {
1427 updatenumdesc(*fork, *root, (*fork)->numdesc - 1);
1428 p = *fork;
1429 while (p->next != item->back)
1430 p = p->next;
1431 p->next = item->back->next;
1432 }
1433 if (!item->tip) {
1434 updatenumdesc(item, item, item->numdesc);
1435 if (recompute) {
1436 memcpy(item->back->oldbase, item->back->base, endsite*sizeof(long));
1437 memcpy(item->back->oldnumsteps, item->back->numsteps, endsite*sizeof(long));
1438 memcpy(item->back->base, zeros, endsite*sizeof(long));
1439 memcpy(item->back->numsteps, zeros, endsite*sizeof(long));
1440 preorder(item, item->back, *root, item->back, NULL, item, -1);
1441 }
1442 }
1443 if ((*fork)->numdesc >= 2)
1444 chuck(grbg, item->back);
1445 item->back = NULL;
1446 if (!recompute)
1447 return;
1448 if ((*fork)->numdesc == 0) {
1449 memcpy(otherback->oldbase, otherback->base, endsite*sizeof(long));
1450 memcpy(otherback->oldnumsteps, otherback->numsteps, endsite*sizeof(long));
1451 if (other == *root) {
1452 memcpy(otherback->base, zeros, endsite*sizeof(long));
1453 memcpy(otherback->numsteps, zeros, endsite*sizeof(long));
1454 } else {
1455 memcpy(otherback->base, other->back->base, endsite*sizeof(long));
1456 memcpy(otherback->numsteps, other->back->numsteps, endsite*sizeof(long));
1457 }
1458 p = other->back;
1459 other->back = otherback;
1460 if (other == *root)
1461 preorder(other, otherback, *root, otherback, NULL, other, -1);
1462 else
1463 preorder(other, otherback, *root, NULL, NULL, NULL, 0);
1464 other->back = p;
1465 if (other != *root) {
1466 memcpy(other->oldbase,(*fork)->base, endsite*sizeof(long));
1467 memcpy(other->oldnumsteps,(*fork)->numsteps, endsite*sizeof(long));
1468 preorder(other->back, other, *root, NULL, NULL, NULL, 0);
1469 }
1470 } else {
1471 memcpy(item->oldbase, item->base, endsite*sizeof(long));
1472 memcpy(item->oldnumsteps, item->numsteps, endsite*sizeof(long));
1473 memcpy(item->base, zeros, endsite*sizeof(long));
1474 memcpy(item->numsteps, zeros, endsite*sizeof(long));
1475 preorder(*fork, item, *root, NULL, NULL, *fork, -1);
1476 if (*fork != *root)
1477 preorder((*fork)->back, *fork, *root, NULL, NULL, NULL, 0);
1478 memcpy(item->base, item->oldbase, endsite*sizeof(long));
1479 memcpy(item->numsteps, item->oldnumsteps, endsite*sizeof(long));
1480 }
1481 } /* remove */
1482
1483
postorder(node * p)1484 void postorder(node *p)
1485 {
1486 /* traverses an n-ary tree, suming up steps at a node's descendants */
1487 /* used in dnacomp, dnapars, & dnapenny */
1488 node *q;
1489
1490 if (p->tip)
1491 return;
1492 q = p->next;
1493 while (q != p) {
1494 postorder(q->back);
1495 q = q->next;
1496 }
1497 zeronumnuc(p, endsite);
1498 if (p->numdesc > 2)
1499 multisumnsteps2(p);
1500 else
1501 fillin(p, p->next->back, p->next->next->back);
1502 } /* postorder */
1503
1504
getnufork(node ** nufork,node ** grbg,pointarray treenode,long * zeros)1505 void getnufork(node **nufork,node **grbg,pointarray treenode,long *zeros)
1506 {
1507 /* find a fork not used currently */
1508 long i;
1509
1510 i = spp;
1511 while (treenode[i] && treenode[i]->numdesc > 0) i++;
1512 if (!treenode[i])
1513 gnutreenode(grbg, &treenode[i], i, endsite, zeros);
1514 *nufork = treenode[i];
1515 } /* getnufork */
1516
1517
reroot(node * outgroup,node * root)1518 void reroot(node *outgroup, node *root)
1519 {
1520 /* reorients tree, putting outgroup in desired position. used if
1521 the root is binary. */
1522 /* used in dnacomp & dnapars */
1523 node *p, *q;
1524
1525 if (outgroup->back->index == root->index)
1526 return;
1527 p = root->next;
1528 q = root->next->next;
1529 p->back->back = q->back;
1530 q->back->back = p->back;
1531 p->back = outgroup;
1532 q->back = outgroup->back;
1533 outgroup->back->back = q;
1534 outgroup->back = p;
1535 } /* reroot */
1536
1537
reroot2(node * outgroup,node * root)1538 void reroot2(node *outgroup, node *root)
1539 {
1540 /* reorients tree, putting outgroup in desired position. */
1541 /* used in dnacomp & dnapars */
1542 node *p;
1543
1544 p = outgroup->back->next;
1545 while (p->next != outgroup->back)
1546 p = p->next;
1547 root->next = outgroup->back;
1548 p->next = root;
1549 } /* reroot2 */
1550
1551
reroot3(node * outgroup,node * root,node * root2,node * lastdesc,node ** grbg)1552 void reroot3(node *outgroup, node *root, node *root2, node *lastdesc,
1553 node **grbg)
1554 {
1555 /* reorients tree, putting back outgroup in original position. */
1556 /* used in dnacomp & dnapars */
1557 node *p;
1558
1559 p = root->next;
1560 while (p->next != root)
1561 p = p->next;
1562 chuck(grbg, root);
1563 p->next = outgroup->back;
1564 root2->next = lastdesc->next;
1565 lastdesc->next = root2;
1566 } /* reroot3 */
1567
1568
savetraverse(node * p)1569 void savetraverse(node *p)
1570 {
1571 /* sets BOOLEANs that indicate which way is down */
1572 node *q;
1573
1574 p->bottom = true;
1575 if (p->tip)
1576 return;
1577 q = p->next;
1578 while (q != p) {
1579 q->bottom = false;
1580 savetraverse(q->back);
1581 q = q->next;
1582 }
1583 } /* savetraverse */
1584
1585
newindex(long i,node * p)1586 void newindex(long i, node *p)
1587 {
1588 /* assigns index i to node p */
1589
1590 while (p->index != i) {
1591 p->index = i;
1592 p = p->next;
1593 }
1594 } /* newindex */
1595
1596
flipindexes(long nextnode,pointarray treenode)1597 void flipindexes(long nextnode, pointarray treenode)
1598 {
1599 /* flips index of nodes between nextnode and last node. */
1600 long last;
1601 node *temp;
1602
1603 last = nonodes;
1604 while (treenode[last - 1]->numdesc == 0)
1605 last--;
1606 if (last > nextnode) {
1607 temp = treenode[nextnode - 1];
1608 treenode[nextnode - 1] = treenode[last - 1];
1609 treenode[last - 1] = temp;
1610 newindex(nextnode, treenode[nextnode - 1]);
1611 newindex(last, treenode[last - 1]);
1612 }
1613 } /* flipindexes */
1614
1615
parentinmulti(node * anode)1616 boolean parentinmulti(node *anode)
1617 {
1618 /* sees if anode's parent has more than 2 children */
1619 node *p;
1620
1621 while (!anode->bottom) anode = anode->next;
1622 p = anode->back;
1623 while (!p->bottom)
1624 p = p->next;
1625 return (p->numdesc > 2);
1626 } /* parentinmulti */
1627
1628
sibsvisited(node * anode,long * place)1629 long sibsvisited(node *anode, long *place)
1630 {
1631 /* computes the number of nodes which are visited earlier than anode among
1632 its siblings */
1633 node *p;
1634 long nvisited;
1635
1636 while (!anode->bottom) anode = anode->next;
1637 p = anode->back->next;
1638 nvisited = 0;
1639 do {
1640 if (!p->bottom && place[p->back->index - 1] != 0)
1641 nvisited++;
1642 p = p->next;
1643 } while (p != anode->back);
1644 return nvisited;
1645 } /* sibsvisited */
1646
1647
smallest(node * anode,long * place)1648 long smallest(node *anode, long *place)
1649 {
1650 /* finds the smallest index of sibling of anode */
1651 node *p;
1652 long min;
1653
1654 while (!anode->bottom) anode = anode->next;
1655 p = anode->back->next;
1656 if (p->bottom) p = p->next;
1657 min = nonodes;
1658 do {
1659 if (p->back && place[p->back->index - 1] != 0) {
1660 if (p->back->index <= spp) {
1661 if (p->back->index < min)
1662 min = p->back->index;
1663 } else {
1664 if (place[p->back->index - 1] < min)
1665 min = place[p->back->index - 1];
1666 }
1667 }
1668 p = p->next;
1669 if (p->bottom) p = p->next;
1670 } while (p != anode->back);
1671 return min;
1672 } /* smallest */
1673
1674
bintomulti(node ** root,node ** binroot,node ** grbg,long * zeros)1675 void bintomulti(node **root, node **binroot, node **grbg, long *zeros)
1676 { /* attaches root's left child to its right child and makes
1677 the right child new root */
1678 node *left, *right, *newnode, *temp;
1679
1680 right = (*root)->next->next->back;
1681 left = (*root)->next->back;
1682 if (right->tip) {
1683 (*root)->next = right->back;
1684 (*root)->next->next = left->back;
1685 temp = left;
1686 left = right;
1687 right = temp;
1688 right->back->next = *root;
1689 }
1690 gnutreenode(grbg, &newnode, right->index, endsite, zeros);
1691 newnode->next = right->next;
1692 newnode->back = left;
1693 left->back = newnode;
1694 right->next = newnode;
1695 (*root)->next->back = (*root)->next->next->back = NULL;
1696 *binroot = *root;
1697 (*binroot)->numdesc = 0;
1698 *root = right;
1699 (*root)->numdesc++;
1700 (*root)->back = NULL;
1701 } /* bintomulti */
1702
1703
backtobinary(node ** root,node * binroot,node ** grbg)1704 void backtobinary(node **root, node *binroot, node **grbg)
1705 { /* restores binary root */
1706 node *p;
1707
1708 binroot->next->back = (*root)->next->back;
1709 (*root)->next->back->back = binroot->next;
1710 p = (*root)->next;
1711 (*root)->next = p->next;
1712 binroot->next->next->back = *root;
1713 (*root)->back = binroot->next->next;
1714 chuck(grbg, p);
1715 (*root)->numdesc--;
1716 *root = binroot;
1717 (*root)->numdesc = 2;
1718 } /* backtobinary */
1719
1720
outgrin(node * root,node * outgrnode)1721 boolean outgrin(node *root, node *outgrnode)
1722 { /* checks if outgroup node is a child of root */
1723 node *p;
1724
1725 p = root->next;
1726 while (p != root) {
1727 if (p->back == outgrnode)
1728 return true;
1729 p = p->next;
1730 }
1731 return false;
1732 } /* outgrin */
1733
1734
flipnodes(node * nodea,node * nodeb)1735 void flipnodes(node *nodea, node *nodeb)
1736 { /* flip nodes */
1737 node *backa, *backb;
1738
1739 backa = nodea->back;
1740 backb = nodeb->back;
1741 backa->back = nodeb;
1742 backb->back = nodea;
1743 nodea->back = backb;
1744 nodeb->back = backa;
1745 } /* flipnodes */
1746
1747
moveleft(node * root,node * outgrnode,node ** flipback)1748 void moveleft(node *root, node *outgrnode, node **flipback)
1749 { /* makes outgroup node to leftmost child of root */
1750 node *p;
1751 boolean done;
1752
1753 p = root->next;
1754 done = false;
1755 while (p != root && !done) {
1756 if (p->back == outgrnode) {
1757 *flipback = p;
1758 flipnodes(root->next->back, p->back);
1759 done = true;
1760 }
1761 p = p->next;
1762 }
1763 } /* moveleft */
1764
1765
savetree(node * root,long * place,pointarray treenode,node ** grbg,long * zeros)1766 void savetree(node *root, long *place, pointarray treenode,
1767 node **grbg, long *zeros)
1768 { /* record in place where each species has to be
1769 added to reconstruct this tree */
1770 /* used by dnacomp & dnapars */
1771 long i, j, nextnode, nvisited;
1772 node *p, *q, *r = NULL, *root2, *lastdesc,
1773 *outgrnode, *binroot, *flipback;
1774 boolean done, newfork;
1775
1776 binroot = NULL;
1777 lastdesc = NULL;
1778 root2 = NULL;
1779 flipback = NULL;
1780 outgrnode = treenode[outgrno - 1];
1781 if (root->numdesc == 2)
1782 bintomulti(&root, &binroot, grbg, zeros);
1783 if (outgrin(root, outgrnode)) {
1784 if (outgrnode != root->next->back)
1785 moveleft(root, outgrnode, &flipback);
1786 } else {
1787 root2 = root;
1788 lastdesc = root->next;
1789 while (lastdesc->next != root)
1790 lastdesc = lastdesc->next;
1791 lastdesc->next = root->next;
1792 gnutreenode(grbg, &root, outgrnode->back->index, endsite, zeros);
1793 root->numdesc = root2->numdesc;
1794 reroot2(outgrnode, root);
1795 }
1796 savetraverse(root);
1797 nextnode = spp + 1;
1798 for (i = nextnode; i <= nonodes; i++)
1799 if (treenode[i - 1]->numdesc == 0)
1800 flipindexes(i, treenode);
1801 for (i = 0; i < nonodes; i++)
1802 place[i] = 0;
1803 place[root->index - 1] = 1;
1804 for (i = 1; i <= spp; i++) {
1805 p = treenode[i - 1];
1806 while (place[p->index - 1] == 0) {
1807 place[p->index - 1] = i;
1808 while (!p->bottom)
1809 p = p->next;
1810 r = p;
1811 p = p->back;
1812 }
1813 if (i > 1) {
1814 q = treenode[i - 1];
1815 newfork = true;
1816 nvisited = sibsvisited(q, place);
1817 if (nvisited == 0) {
1818 if (parentinmulti(r)) {
1819 nvisited = sibsvisited(r, place);
1820 if (nvisited == 0)
1821 place[i - 1] = place[p->index - 1];
1822 else if (nvisited == 1)
1823 place[i - 1] = smallest(r, place);
1824 else {
1825 place[i - 1] = -smallest(r, place);
1826 newfork = false;
1827 }
1828 } else
1829 place[i - 1] = place[p->index - 1];
1830 } else if (nvisited == 1) {
1831 place[i - 1] = place[p->index - 1];
1832 } else {
1833 place[i - 1] = -smallest(q, place);
1834 newfork = false;
1835 }
1836 if (newfork) {
1837 j = place[p->index - 1];
1838 done = false;
1839 while (!done) {
1840 place[p->index - 1] = nextnode;
1841 while (!p->bottom)
1842 p = p->next;
1843 p = p->back;
1844 done = (p == NULL);
1845 if (!done)
1846 done = (place[p->index - 1] != j);
1847 if (done) {
1848 nextnode++;
1849 }
1850 }
1851 }
1852 }
1853 }
1854 if (flipback)
1855 flipnodes(outgrnode, flipback->back);
1856 else {
1857 if (root2) {
1858 reroot3(outgrnode, root, root2, lastdesc, grbg);
1859 root = root2;
1860 }
1861 }
1862 if (binroot)
1863 backtobinary(&root, binroot, grbg);
1864 } /* savetree */
1865
1866
addnsave(node * p,node * item,node * nufork,node ** root,node ** grbg,boolean multf,pointarray treenode,long * place,long * zeros)1867 void addnsave(node *p, node *item, node *nufork, node **root, node **grbg,
1868 boolean multf, pointarray treenode, long *place, long *zeros)
1869 { /* adds item to tree and save it. Then removes item. */
1870 node *dummy;
1871
1872 if (!multf)
1873 add(p, item, nufork, root, false, treenode, grbg, zeros);
1874 else
1875 add(p, item, NULL, root, false, treenode, grbg, zeros);
1876 savetree(*root, place, treenode, grbg, zeros);
1877 if (!multf)
1878 re_move(item, &nufork, root, false, treenode, grbg, zeros);
1879 else
1880 re_move(item, &dummy, root, false, treenode, grbg, zeros);
1881 } /* addnsave */
1882
1883
addbestever(long * pos,long * nextree,long maxtrees,boolean collapse,long * place,bestelm * bestrees)1884 void addbestever(long *pos, long *nextree, long maxtrees, boolean collapse,
1885 long *place, bestelm *bestrees)
1886 { /* adds first best tree */
1887
1888 *pos = 1;
1889 *nextree = 1;
1890 initbestrees(bestrees, maxtrees, true);
1891 initbestrees(bestrees, maxtrees, false);
1892 addtree(*pos, nextree, collapse, place, bestrees);
1893 } /* addbestever */
1894
1895
addtiedtree(long pos,long * nextree,long maxtrees,boolean collapse,long * place,bestelm * bestrees)1896 void addtiedtree(long pos, long *nextree, long maxtrees, boolean collapse,
1897 long *place, bestelm *bestrees)
1898 { /* add tied tree */
1899
1900 if (*nextree <= maxtrees)
1901 addtree(pos, nextree, collapse, place, bestrees);
1902 } /* addtiedtree */
1903
1904
clearcollapse(pointarray treenode)1905 void clearcollapse(pointarray treenode)
1906 {
1907 /* clears collapse status at a node */
1908 long i;
1909 node *p;
1910
1911 for (i = 0; i < nonodes; i++) {
1912 treenode[i]->collapse = undefined;
1913 if (!treenode[i]->tip) {
1914 p = treenode[i]->next;
1915 while (p != treenode[i]) {
1916 p->collapse = undefined;
1917 p = p->next;
1918 }
1919 }
1920 }
1921 } /* clearcollapse */
1922
1923
clearbottom(pointarray treenode)1924 void clearbottom(pointarray treenode)
1925 {
1926 /* clears boolean bottom at a node */
1927 long i;
1928 node *p;
1929
1930 for (i = 0; i < nonodes; i++) {
1931 treenode[i]->bottom = false;
1932 if (!treenode[i]->tip) {
1933 p = treenode[i]->next;
1934 while (p != treenode[i]) {
1935 p->bottom = false;
1936 p = p->next;
1937 }
1938 }
1939 }
1940 } /* clearbottom */
1941
1942
collabranch(node * collapfrom,node * tempfrom,node * tempto)1943 void collabranch(node *collapfrom, node *tempfrom, node *tempto)
1944 { /* collapse branch from collapfrom */
1945 long i, j, b, largest, descsteps;
1946 boolean done;
1947
1948 for (i = 0; i < endsite; i++) {
1949 descsteps = 0;
1950 for (j = (long)A; j <= (long)O; j++) {
1951 b = 1 << j;
1952 if ((descsteps == 0) && (collapfrom->base[i] & b))
1953 descsteps = tempfrom->oldnumsteps[i]
1954 - (collapfrom->numdesc - collapfrom->numnuc[i][j])
1955 * weight[i];
1956 }
1957 done = false;
1958 for (j = (long)A; j <= (long)O; j++) {
1959 b = 1 << j;
1960 if (!done && (tempto->base[i] & b)) {
1961 descsteps += (tempto->numsteps[i]
1962 - (tempto->numdesc - collapfrom->numdesc
1963 - tempto->numnuc[i][j]) * weight[i]);
1964 done = true;
1965 }
1966 }
1967 for (j = (long)A; j <= (long)O; j++)
1968 tempto->numnuc[i][j] += collapfrom->numnuc[i][j];
1969 largest = getlargest(tempto->numnuc[i]);
1970 tempto->base[i] = 0;
1971 for (j = (long)A; j <= (long)O; j++) {
1972 if (tempto->numnuc[i][j] == largest)
1973 tempto->base[i] |= (1 << j);
1974 }
1975 tempto->numsteps[i] = (tempto->numdesc - largest) * weight[i] + descsteps;
1976 }
1977 } /* collabranch */
1978
1979
allcommonbases(node * a,node * b,boolean * allsame)1980 boolean allcommonbases(node *a, node *b, boolean *allsame)
1981 { /* see if bases are common at all sites for nodes a and b */
1982 long i;
1983 boolean allcommon;
1984
1985 allcommon = true;
1986 *allsame = true;
1987 for (i = 0; i < endsite; i++) {
1988 if ((a->base[i] & b->base[i]) == 0)
1989 allcommon = false;
1990 else if (a->base[i] != b->base[i])
1991 *allsame = false;
1992 }
1993 return allcommon;
1994 } /* allcommonbases */
1995
1996
findbottom(node * p,node ** bottom)1997 void findbottom(node *p, node **bottom)
1998 { /* find a node with field bottom set at node p */
1999 node *q;
2000
2001 if (p->bottom)
2002 *bottom = p;
2003 else {
2004 q = p->next;
2005 while(!q->bottom && q != p)
2006 q = q->next;
2007 *bottom = q;
2008 }
2009 } /* findbottom */
2010
2011
moresteps(node * a,node * b)2012 boolean moresteps(node *a, node *b)
2013 { /* see if numsteps of node a exceeds those of node b */
2014 long i;
2015
2016 for (i = 0; i < endsite; i++)
2017 if (a->numsteps[i] > b->numsteps[i])
2018 return true;
2019 return false;
2020 } /* moresteps */
2021
2022
passdown(node * desc,node * parent,node * start,node * below,node * item,node * added,node * total,node * tempdsc,node * tempprt,boolean multf)2023 boolean passdown(node *desc, node *parent, node *start, node *below,
2024 node *item, node *added, node *total, node *tempdsc,
2025 node *tempprt, boolean multf)
2026 { /* track down to node start to see if an ancestor branch can be collapsed */
2027 node *temp;
2028 boolean done, allsame;
2029
2030 done = (parent == start);
2031 while (!done) {
2032 desc = parent;
2033 findbottom(parent->back, &parent);
2034 if (multf && start == below && parent == below)
2035 parent = added;
2036 memcpy(tempdsc->base, tempprt->base, endsite*sizeof(long));
2037 memcpy(tempdsc->numsteps, tempprt->numsteps, endsite*sizeof(long));
2038 memcpy(tempdsc->oldbase, desc->base, endsite*sizeof(long));
2039 memcpy(tempdsc->oldnumsteps, desc->numsteps, endsite*sizeof(long));
2040 memcpy(tempprt->base, parent->base, endsite*sizeof(long));
2041 memcpy(tempprt->numsteps, parent->numsteps, endsite*sizeof(long));
2042 memcpy(tempprt->numnuc, parent->numnuc, endsite*sizeof(nucarray));
2043 tempprt->numdesc = parent->numdesc;
2044 multifillin(tempprt, tempdsc, 0);
2045 if (!allcommonbases(tempprt, parent, &allsame))
2046 return false;
2047 else if (moresteps(tempprt, parent))
2048 return false;
2049 else if (allsame)
2050 return true;
2051 if (parent == added)
2052 parent = below;
2053 done = (parent == start);
2054 if (done && ((start == item) || (!multf && start == below))) {
2055 memcpy(tempdsc->base, tempprt->base, endsite*sizeof(long));
2056 memcpy(tempdsc->numsteps, tempprt->numsteps, endsite*sizeof(long));
2057 memcpy(tempdsc->oldbase, start->base, endsite*sizeof(long));
2058 memcpy(tempdsc->oldnumsteps, start->numsteps, endsite*sizeof(long));
2059 multifillin(added, tempdsc, 0);
2060 tempprt = added;
2061 }
2062 }
2063 temp = tempdsc;
2064 if (start == below || start == item)
2065 fillin(temp, tempprt, below->back);
2066 else
2067 fillin(temp, tempprt, added);
2068 return !moresteps(temp, total);
2069 } /* passdown */
2070
2071
trycollapdesc(node * desc,node * parent,node * start,node * below,node * item,node * added,node * total,node * tempdsc,node * tempprt,boolean multf,long * zeros)2072 boolean trycollapdesc(node *desc, node *parent, node *start,
2073 node *below, node *item, node *added, node *total,
2074 node *tempdsc, node *tempprt, boolean multf, long *zeros)
2075 { /* see if branch between nodes desc and parent can be collapsed */
2076 boolean allsame;
2077
2078 if (desc->numdesc == 1)
2079 return true;
2080 if (multf && start == below && parent == below)
2081 parent = added;
2082 memcpy(tempdsc->base, zeros, endsite*sizeof(long));
2083 memcpy(tempdsc->numsteps, zeros, endsite*sizeof(long));
2084 memcpy(tempdsc->oldbase, desc->base, endsite*sizeof(long));
2085 memcpy(tempdsc->oldnumsteps, desc->numsteps, endsite*sizeof(long));
2086 memcpy(tempprt->base, parent->base, endsite*sizeof(long));
2087 memcpy(tempprt->numsteps, parent->numsteps, endsite*sizeof(long));
2088 memcpy(tempprt->numnuc, parent->numnuc, endsite*sizeof(nucarray));
2089 tempprt->numdesc = parent->numdesc - 1;
2090 multifillin(tempprt, tempdsc, -1);
2091 tempprt->numdesc += desc->numdesc;
2092 collabranch(desc, tempdsc, tempprt);
2093 if (!allcommonbases(tempprt, parent, &allsame) ||
2094 moresteps(tempprt, parent)) {
2095 if (parent != added) {
2096 desc->collapse = nocollap;
2097 parent->collapse = nocollap;
2098 }
2099 return false;
2100 } else if (allsame) {
2101 if (parent != added) {
2102 desc->collapse = tocollap;
2103 parent->collapse = tocollap;
2104 }
2105 return true;
2106 }
2107 if (parent == added)
2108 parent = below;
2109 if ((start == item && parent == item) ||
2110 (!multf && start == below && parent == below)) {
2111 memcpy(tempdsc->base, tempprt->base, endsite*sizeof(long));
2112 memcpy(tempdsc->numsteps, tempprt->numsteps, endsite*sizeof(long));
2113 memcpy(tempdsc->oldbase, start->base, endsite*sizeof(long));
2114 memcpy(tempdsc->oldnumsteps, start->numsteps, endsite*sizeof(long));
2115 memcpy(tempprt->base, added->base, endsite*sizeof(long));
2116 memcpy(tempprt->numsteps, added->numsteps, endsite*sizeof(long));
2117 memcpy(tempprt->numnuc, added->numnuc, endsite*sizeof(nucarray));
2118 tempprt->numdesc = added->numdesc;
2119 multifillin(tempprt, tempdsc, 0);
2120 if (!allcommonbases(tempprt, added, &allsame))
2121 return false;
2122 else if (moresteps(tempprt, added))
2123 return false;
2124 else if (allsame)
2125 return true;
2126 }
2127 return passdown(desc, parent, start, below, item, added, total, tempdsc,
2128 tempprt, multf);
2129 } /* trycollapdesc */
2130
2131
setbottom(node * p)2132 void setbottom(node *p)
2133 { /* set field bottom at node p */
2134 node *q;
2135
2136 p->bottom = true;
2137 q = p->next;
2138 do {
2139 q->bottom = false;
2140 q = q->next;
2141 } while (q != p);
2142 } /* setbottom */
2143
zeroinsubtree(node * subtree,node * start,node * below,node * item,node * added,node * total,node * tempdsc,node * tempprt,boolean multf,node * root,long * zeros)2144 boolean zeroinsubtree(node *subtree, node *start, node *below, node *item,
2145 node *added, node *total, node *tempdsc, node *tempprt,
2146 boolean multf, node* root, long *zeros)
2147 { /* sees if subtree contains a zero length branch */
2148 node *p;
2149
2150 if (!subtree->tip) {
2151 setbottom(subtree);
2152 p = subtree->next;
2153 do {
2154 if (p->back && !p->back->tip &&
2155 !((p->back->collapse == nocollap) && (subtree->collapse == nocollap))
2156 && (subtree->numdesc != 1)) {
2157 if ((p->back->collapse == tocollap) && (subtree->collapse == tocollap)
2158 && multf && (subtree != below))
2159 return true;
2160 /* when root->numdesc == 2
2161 * there is no mandatory step at the root,
2162 * instead of checking at the root we check around it
2163 * we only need to check p because the first if
2164 * statement already gets rid of it for the subtree */
2165 else if ((p->back->index != root->index || root->numdesc > 2) &&
2166 trycollapdesc(p->back, subtree, start, below, item, added, total,
2167 tempdsc, tempprt, multf, zeros))
2168 return true;
2169 else if ((p->back->index == root->index && root->numdesc == 2) &&
2170 !(root->next->back->tip) && !(root->next->next->back->tip) &&
2171 trycollapdesc(root->next->back, root->next->next->back, start,
2172 below, item,added, total, tempdsc, tempprt, multf, zeros))
2173 return true;
2174 }
2175 p = p->next;
2176 } while (p != subtree);
2177 p = subtree->next;
2178 do {
2179 if (p->back && !p->back->tip) {
2180 if (zeroinsubtree(p->back, start, below, item, added, total,
2181 tempdsc, tempprt, multf, root, zeros))
2182 return true;
2183 }
2184 p = p->next;
2185 } while (p != subtree);
2186 }
2187 return false;
2188 } /* zeroinsubtree */
2189
2190
collapsible(node * item,node * below,node * temp,node * temp1,node * tempdsc,node * tempprt,node * added,node * total,boolean multf,node * root,long * zeros,pointarray treenode)2191 boolean collapsible(node *item, node *below, node *temp, node *temp1,
2192 node *tempdsc, node *tempprt, node *added, node *total,
2193 boolean multf, node *root, long *zeros, pointarray treenode)
2194 {
2195 /* sees if any branch can be collapsed */
2196 node *belowbk;
2197 boolean allsame;
2198
2199 if (multf) {
2200 memcpy(tempdsc->base, item->base, endsite*sizeof(long));
2201 memcpy(tempdsc->numsteps, item->numsteps, endsite*sizeof(long));
2202 memcpy(tempdsc->oldbase, zeros, endsite*sizeof(long));
2203 memcpy(tempdsc->oldnumsteps, zeros, endsite*sizeof(long));
2204 memcpy(added->base, below->base, endsite*sizeof(long));
2205 memcpy(added->numsteps, below->numsteps, endsite*sizeof(long));
2206 memcpy(added->numnuc, below->numnuc, endsite*sizeof(nucarray));
2207 added->numdesc = below->numdesc + 1;
2208 multifillin(added, tempdsc, 1);
2209 } else {
2210 fillin(added, item, below);
2211 added->numdesc = 2;
2212 }
2213 fillin(total, added, below->back);
2214 clearbottom(treenode);
2215 if (below->back) {
2216 if (zeroinsubtree(below->back, below->back, below, item, added, total,
2217 tempdsc, tempprt, multf, root, zeros))
2218 return true;
2219 }
2220 if (multf) {
2221 if (zeroinsubtree(below, below, below, item, added, total,
2222 tempdsc, tempprt, multf, root, zeros))
2223 return true;
2224 } else if (!below->tip) {
2225 if (zeroinsubtree(below, below, below, item, added, total,
2226 tempdsc, tempprt, multf, root, zeros))
2227 return true;
2228 }
2229 if (!item->tip) {
2230 if (zeroinsubtree(item, item, below, item, added, total,
2231 tempdsc, tempprt, multf, root, zeros))
2232 return true;
2233 }
2234 if (multf && below->back && !below->back->tip) {
2235 memcpy(tempdsc->base, zeros, endsite*sizeof(long));
2236 memcpy(tempdsc->numsteps, zeros, endsite*sizeof(long));
2237 memcpy(tempdsc->oldbase, added->base, endsite*sizeof(long));
2238 memcpy(tempdsc->oldnumsteps, added->numsteps, endsite*sizeof(long));
2239 if (below->back == treenode[below->back->index - 1])
2240 belowbk = below->back->next;
2241 else
2242 belowbk = treenode[below->back->index - 1];
2243 memcpy(tempprt->base, belowbk->base, endsite*sizeof(long));
2244 memcpy(tempprt->numsteps, belowbk->numsteps, endsite*sizeof(long));
2245 memcpy(tempprt->numnuc, belowbk->numnuc, endsite*sizeof(nucarray));
2246 tempprt->numdesc = belowbk->numdesc - 1;
2247 multifillin(tempprt, tempdsc, -1);
2248 tempprt->numdesc += added->numdesc;
2249 collabranch(added, tempdsc, tempprt);
2250 if (!allcommonbases(tempprt, belowbk, &allsame))
2251 return false;
2252 else if (allsame && !moresteps(tempprt, belowbk))
2253 return true;
2254 else if (belowbk->back) {
2255 fillin(temp, tempprt, belowbk->back);
2256 fillin(temp1, belowbk, belowbk->back);
2257 return !moresteps(temp, temp1);
2258 }
2259 }
2260 return false;
2261 } /* collapsible */
2262
2263
replaceback(node ** oldback,node * item,node * forknode,node ** grbg,long * zeros)2264 void replaceback(node **oldback, node *item, node *forknode,
2265 node **grbg, long *zeros)
2266 { /* replaces back node of item with another */
2267 node *p;
2268
2269 p = forknode;
2270 while (p->next->back != item)
2271 p = p->next;
2272 *oldback = p->next;
2273 gnutreenode(grbg, &p->next, forknode->index, endsite, zeros);
2274 p->next->next = (*oldback)->next;
2275 p->next->back = (*oldback)->back;
2276 p->next->back->back = p->next;
2277 (*oldback)->next = (*oldback)->back = NULL;
2278 } /* replaceback */
2279
2280
putback(node * oldback,node * item,node * forknode,node ** grbg)2281 void putback(node *oldback, node *item, node *forknode, node **grbg)
2282 { /* restores node to back of item */
2283 node *p, *q;
2284
2285 p = forknode;
2286 while (p->next != item->back)
2287 p = p->next;
2288 q = p->next;
2289 oldback->next = p->next->next;
2290 p->next = oldback;
2291 oldback->back = item;
2292 item->back = oldback;
2293 oldback->index = forknode->index;
2294 chuck(grbg, q);
2295 } /* putback */
2296
2297
savelocrearr(node * item,node * forknode,node * below,node * tmp,node * tmp1,node * tmp2,node * tmp3,node * tmprm,node * tmpadd,node ** root,long maxtrees,long * nextree,boolean multf,boolean bestever,boolean * saved,long * place,bestelm * bestrees,pointarray treenode,node ** grbg,long * zeros)2298 void savelocrearr(node *item, node *forknode, node *below, node *tmp,
2299 node *tmp1, node *tmp2, node *tmp3, node *tmprm, node *tmpadd,
2300 node **root, long maxtrees, long *nextree, boolean multf,
2301 boolean bestever, boolean *saved, long *place,
2302 bestelm *bestrees, pointarray treenode, node **grbg,
2303 long *zeros)
2304 { /* saves tied or better trees during local rearrangements by removing
2305 item from forknode and adding to below */
2306 node *other, *otherback = NULL, *oldfork, *nufork, *oldback;
2307 long pos;
2308 boolean found, collapse;
2309
2310 if (forknode->numdesc == 2) {
2311 findbelow(&other, item, forknode);
2312 otherback = other->back;
2313 oldback = NULL;
2314 } else {
2315 other = NULL;
2316 replaceback(&oldback, item, forknode, grbg, zeros);
2317 }
2318 re_move(item, &oldfork, root, false, treenode, grbg, zeros);
2319 if (!multf)
2320 getnufork(&nufork, grbg, treenode, zeros);
2321 else
2322 nufork = NULL;
2323 addnsave(below, item, nufork, root, grbg, multf, treenode, place, zeros);
2324 pos = 0;
2325 findtree(&found, &pos, *nextree, place, bestrees);
2326 if (other) {
2327 add(other, item, oldfork, root, false, treenode, grbg, zeros);
2328 if (otherback->back != other)
2329 flipnodes(item, other);
2330 } else
2331 add(forknode, item, NULL, root, false, treenode, grbg, zeros);
2332 *saved = false;
2333 if (found) {
2334 if (oldback)
2335 putback(oldback, item, forknode, grbg);
2336 } else {
2337 if (oldback)
2338 chuck(grbg, oldback);
2339 re_move(item, &oldfork, root, true, treenode, grbg, zeros);
2340 collapse = collapsible(item, below, tmp, tmp1, tmp2, tmp3, tmprm,
2341 tmpadd, multf, *root, zeros, treenode);
2342 if (!collapse) {
2343 if (bestever)
2344 addbestever(&pos, nextree, maxtrees, collapse, place, bestrees);
2345 else
2346 addtiedtree(pos, nextree, maxtrees, collapse, place, bestrees);
2347 }
2348 if (other)
2349 add(other, item, oldfork, root, true, treenode, grbg, zeros);
2350 else
2351 add(forknode, item, NULL, root, true, treenode, grbg, zeros);
2352 *saved = !collapse;
2353 }
2354 } /* savelocrearr */
2355
2356
clearvisited(pointarray treenode)2357 void clearvisited(pointarray treenode)
2358 {
2359 /* clears boolean visited at a node */
2360 long i;
2361 node *p;
2362
2363 for (i = 0; i < nonodes; i++) {
2364 treenode[i]->visited = false;
2365 if (!treenode[i]->tip) {
2366 p = treenode[i]->next;
2367 while (p != treenode[i]) {
2368 p->visited = false;
2369 p = p->next;
2370 }
2371 }
2372 }
2373 } /* clearvisited */
2374
2375
2376 /*void hyprint(long b1, long b2, struct LOC_hyptrav *htrav,
2377 pointarray treenode, Char *basechar)
2378 {
2379 * print out states in sites b1 through b2 at node *
2380 long i, j, k, n;
2381 boolean dot;
2382 bases b;
2383
2384 if (htrav->bottom) {
2385 if (!outgropt)
2386 fprintf(outfile, " ");
2387 else
2388 fprintf(outfile, "root ");
2389 } else
2390 fprintf(outfile, "%4ld ", htrav->r->back->index - spp);
2391 if (htrav->r->tip) {
2392 for (i = 0; i < nmlngth; i++)
2393 putc(nayme[htrav->r->index - 1][i], outfile);
2394 } else
2395 fprintf(outfile, "%4ld ", htrav->r->index - spp);
2396 if (htrav->bottom)
2397 fprintf(outfile, " ");
2398 else if (htrav->nonzero)
2399 fprintf(outfile, " yes ");
2400 else if (htrav->maybe)
2401 fprintf(outfile, " maybe ");
2402 else
2403 fprintf(outfile, " no ");
2404 for (i = b1; i <= b2; i++) {
2405 j = location[ally[i - 1] - 1];
2406 htrav->tempset = htrav->r->base[j - 1];
2407 htrav->anc = htrav->hypset[j - 1];
2408 if (!htrav->bottom)
2409 htrav->anc = treenode[htrav->r->back->index - 1]->base[j - 1];
2410 dot = dotdiff && (htrav->tempset == htrav->anc && !htrav->bottom);
2411 if (dot)
2412 putc('.', outfile);
2413 else if (htrav->tempset == (1 << A))
2414 putc('A', outfile);
2415 else if (htrav->tempset == (1 << C))
2416 putc('C', outfile);
2417 else if (htrav->tempset == (1 << G))
2418 putc('G', outfile);
2419 else if (htrav->tempset == (1 << T))
2420 putc('T', outfile);
2421 else if (htrav->tempset == (1 << O))
2422 putc('-', outfile);
2423 else {
2424 k = 1;
2425 n = 0;
2426 for (b = A; b <= O; b = b + 1) {
2427 if (((1 << b) & htrav->tempset) != 0)
2428 n += k;
2429 k += k;
2430 }
2431 putc(basechar[n - 1], outfile);
2432 }
2433 if (i % 10 == 0)
2434 putc(' ', outfile);
2435 }
2436 putc('\n', outfile);
2437 } * hyprint */
2438
2439
gnubase(gbases ** p,gbases ** garbage,long endsite)2440 void gnubase(gbases **p, gbases **garbage, long endsite)
2441 {
2442 /* this and the following are do-it-yourself garbage collectors.
2443 Make a new node or pull one off the garbage list */
2444 if (*garbage != NULL) {
2445 *p = *garbage;
2446 *garbage = (*garbage)->next;
2447 } else {
2448 *p = (gbases *)Malloc(sizeof(gbases));
2449 (*p)->base = (baseptr)Malloc(endsite*sizeof(long));
2450 }
2451 (*p)->next = NULL;
2452 } /* gnubase */
2453
2454
chuckbase(gbases * p,gbases ** garbage)2455 void chuckbase(gbases *p, gbases **garbage)
2456 {
2457 /* collect garbage on p -- put it on front of garbage list */
2458 p->next = *garbage;
2459 *garbage = p;
2460 } /* chuckbase */
2461
2462
2463 /*void hyptrav(node *r_, long *hypset_, long b1, long b2, boolean bottom_,
2464 pointarray treenode, gbases **garbage, Char *basechar)
2465 {
2466 * compute, print out states at one interior node *
2467 struct LOC_hyptrav Vars;
2468 long i, j, k;
2469 long largest;
2470 gbases *ancset;
2471 nucarray *tempnuc;
2472 node *p, *q;
2473
2474 Vars.bottom = bottom_;
2475 Vars.r = r_;
2476 Vars.hypset = hypset_;
2477 gnubase(&ancset, garbage, endsite);
2478 tempnuc = (nucarray *)Malloc(endsite*sizeof(nucarray));
2479 Vars.maybe = false;
2480 Vars.nonzero = false;
2481 if (!Vars.r->tip)
2482 zeronumnuc(Vars.r, endsite);
2483 for (i = b1 - 1; i < b2; i++) {
2484 j = location[ally[i] - 1];
2485 Vars.anc = Vars.hypset[j - 1];
2486 if (!Vars.r->tip) {
2487 p = Vars.r->next;
2488 for (k = (long)A; k <= (long)O; k++)
2489 if (Vars.anc & (1 << k))
2490 Vars.r->numnuc[j - 1][k]++;
2491 do {
2492 for (k = (long)A; k <= (long)O; k++)
2493 if (p->back->base[j - 1] & (1 << k))
2494 Vars.r->numnuc[j - 1][k]++;
2495 p = p->next;
2496 } while (p != Vars.r);
2497 largest = getlargest(Vars.r->numnuc[j - 1]);
2498 Vars.tempset = 0;
2499 for (k = (long)A; k <= (long)O; k++) {
2500 if (Vars.r->numnuc[j - 1][k] == largest)
2501 Vars.tempset |= (1 << k);
2502 }
2503 Vars.r->base[j - 1] = Vars.tempset;
2504 }
2505 if (!Vars.bottom)
2506 Vars.anc = treenode[Vars.r->back->index - 1]->base[j - 1];
2507 Vars.nonzero = (Vars.nonzero || (Vars.r->base[j - 1] & Vars.anc) == 0);
2508 Vars.maybe = (Vars.maybe || Vars.r->base[j - 1] != Vars.anc);
2509 }
2510 hyprint(b1, b2, &Vars, treenode, basechar);
2511 Vars.bottom = false;
2512 if (!Vars.r->tip) {
2513 memcpy(tempnuc, Vars.r->numnuc, endsite*sizeof(nucarray));
2514 q = Vars.r->next;
2515 do {
2516 memcpy(Vars.r->numnuc, tempnuc, endsite*sizeof(nucarray));
2517 for (i = b1 - 1; i < b2; i++) {
2518 j = location[ally[i] - 1];
2519 for (k = (long)A; k <= (long)O; k++)
2520 if (q->back->base[j - 1] & (1 << k))
2521 Vars.r->numnuc[j - 1][k]--;
2522 largest = getlargest(Vars.r->numnuc[j - 1]);
2523 ancset->base[j - 1] = 0;
2524 for (k = (long)A; k <= (long)O; k++)
2525 if (Vars.r->numnuc[j - 1][k] == largest)
2526 ancset->base[j - 1] |= (1 << k);
2527 if (!Vars.bottom)
2528 Vars.anc = ancset->base[j - 1];
2529 }
2530 hyptrav(q->back, ancset->base, b1, b2, Vars.bottom,
2531 treenode, garbage, basechar);
2532 q = q->next;
2533 } while (q != Vars.r);
2534 }
2535 chuckbase(ancset, garbage);
2536 } * hyptrav *
2537
2538
2539 void hypstates(long chars, node *root, pointarray treenode,
2540 gbases **garbage, Char *basechar)
2541 {
2542 * fill in and describe states at interior nodes *
2543 * used in dnacomp, dnapars, & dnapenny *
2544 long i, n;
2545 baseptr nothing;
2546
2547 fprintf(outfile, "\nFrom To Any Steps? State at upper node\n");
2548 fprintf(outfile, " ");
2549 if (dotdiff)
2550 fprintf(outfile, " ( . means same as in the node below it on tree)\n");
2551 nothing = (baseptr)Malloc(endsite*sizeof(long));
2552 for (i = 0; i < endsite; i++)
2553 nothing[i] = 0;
2554 for (i = 1; i <= ((chars - 1) / 40 + 1); i++) {
2555 putc('\n', outfile);
2556 n = i * 40;
2557 if (n > chars)
2558 n = chars;
2559 hyptrav(root, nothing, i * 40 - 39, n, true, treenode, garbage, basechar);
2560 }
2561 free(nothing);
2562 } * hypstates */
2563
2564
initbranchlen(node * p)2565 void initbranchlen(node *p)
2566 {
2567 node *q;
2568
2569 p->v = 0.0;
2570 if (p->back)
2571 p->back->v = 0.0;
2572 if (p->tip)
2573 return;
2574 q = p->next;
2575 while (q != p) {
2576 initbranchlen(q->back);
2577 q = q->next;
2578 }
2579 q = p->next;
2580 while (q != p) {
2581 q->v = 0.0;
2582 q = q->next;
2583 }
2584 } /* initbranchlen */
2585
2586
initmin(node * p,long sitei,boolean internal)2587 void initmin(node *p, long sitei, boolean internal)
2588 {
2589 long i;
2590
2591 if (internal) {
2592 for (i = (long)A; i <= (long)O; i++) {
2593 p->cumlengths[i] = 0;
2594 p->numreconst[i] = 1;
2595 }
2596 } else {
2597 for (i = (long)A; i <= (long)O; i++) {
2598 if (p->base[sitei - 1] & (1 << i)) {
2599 p->cumlengths[i] = 0;
2600 p->numreconst[i] = 1;
2601 } else {
2602 p->cumlengths[i] = -1;
2603 p->numreconst[i] = 0;
2604 }
2605 }
2606 }
2607 } /* initmin */
2608
2609
initbase(node * p,long sitei)2610 void initbase(node *p, long sitei)
2611 {
2612 /* traverse tree to initialize base at internal nodes */
2613 node *q;
2614 long i, largest;
2615
2616 if (p->tip)
2617 return;
2618 q = p->next;
2619 while (q != p) {
2620 if (q->back) {
2621 memcpy(q->numnuc, p->numnuc, endsite*sizeof(nucarray));
2622 for (i = (long)A; i <= (long)O; i++) {
2623 if (q->back->base[sitei - 1] & (1 << i))
2624 q->numnuc[sitei - 1][i]--;
2625 }
2626 if (p->back) {
2627 for (i = (long)A; i <= (long)O; i++) {
2628 if (p->back->base[sitei - 1] & (1 << i))
2629 q->numnuc[sitei - 1][i]++;
2630 }
2631 }
2632 largest = getlargest(q->numnuc[sitei - 1]);
2633 q->base[sitei - 1] = 0;
2634 for (i = (long)A; i <= (long)O; i++) {
2635 if (q->numnuc[sitei - 1][i] == largest)
2636 q->base[sitei - 1] |= (1 << i);
2637 }
2638 }
2639 q = q->next;
2640 }
2641 q = p->next;
2642 while (q != p) {
2643 initbase(q->back, sitei);
2644 q = q->next;
2645 }
2646 } /* initbase */
2647
2648
inittreetrav(node * p,long sitei)2649 void inittreetrav(node *p, long sitei)
2650 {
2651 /* traverse tree to clear boolean initialized and set up base */
2652 node *q;
2653
2654 if (p->tip) {
2655 initmin(p, sitei, false);
2656 p->initialized = true;
2657 return;
2658 }
2659 q = p->next;
2660 while (q != p) {
2661 inittreetrav(q->back, sitei);
2662 q = q->next;
2663 }
2664 initmin(p, sitei, true);
2665 p->initialized = false;
2666 q = p->next;
2667 while (q != p) {
2668 initmin(q, sitei, true);
2669 q->initialized = false;
2670 q = q->next;
2671 }
2672 } /* inittreetrav */
2673
2674
compmin(node * p,node * desc)2675 void compmin(node *p, node *desc)
2676 {
2677 /* computes minimum lengths up to p */
2678 long i, j, minn, cost, desclen, descrecon=0, maxx;
2679
2680 maxx = 10 * spp;
2681 for (i = (long)A; i <= (long)O; i++) {
2682 minn = maxx;
2683 for (j = (long)A; j <= (long)O; j++) {
2684 if (transvp) {
2685 if (
2686 (
2687 ((i == (long)A) || (i == (long)G))
2688 && ((j == (long)A) || (j == (long)G))
2689 )
2690 || (
2691 ((j == (long)C) || (j == (long)T))
2692 && ((i == (long)C) || (i == (long)T))
2693 )
2694 )
2695 cost = 0;
2696 else
2697 cost = 1;
2698 } else {
2699 if (i == j)
2700 cost = 0;
2701 else
2702 cost = 1;
2703 }
2704 if (desc->cumlengths[j] == -1) {
2705 desclen = maxx;
2706 } else {
2707 desclen = desc->cumlengths[j];
2708 }
2709 if (minn > cost + desclen) {
2710 minn = cost + desclen;
2711 descrecon = 0;
2712 }
2713 if (minn == cost + desclen) {
2714 descrecon += desc->numreconst[j];
2715 }
2716 }
2717 p->cumlengths[i] += minn;
2718 p->numreconst[i] *= descrecon;
2719 }
2720 p->initialized = true;
2721 } /* compmin */
2722
2723
minpostorder(node * p,pointarray treenode)2724 void minpostorder(node *p, pointarray treenode)
2725 {
2726 /* traverses an n-ary tree, computing minimum steps at each node */
2727 node *q;
2728
2729 if (p->tip) {
2730 return;
2731 }
2732 q = p->next;
2733 while (q != p) {
2734 if (q->back)
2735 minpostorder(q->back, treenode);
2736 q = q->next;
2737 }
2738 if (!p->initialized) {
2739 q = p->next;
2740 while (q != p) {
2741 if (q->back)
2742 compmin(p, q->back);
2743 q = q->next;
2744 }
2745 }
2746 } /* minpostorder */
2747
2748
branchlength(node * subtr1,node * subtr2,double * brlen,pointarray treenode)2749 void branchlength(node *subtr1, node *subtr2, double *brlen,
2750 pointarray treenode)
2751 {
2752 /* computes a branch length between two subtrees for a given site */
2753 long i, j, minn, cost, nom, denom;
2754 node *temp;
2755
2756 if (subtr1->tip) {
2757 temp = subtr1;
2758 subtr1 = subtr2;
2759 subtr2 = temp;
2760 }
2761 if (subtr1->index == outgrno) {
2762 temp = subtr1;
2763 subtr1 = subtr2;
2764 subtr2 = temp;
2765 }
2766 minpostorder(subtr1, treenode);
2767 minpostorder(subtr2, treenode);
2768 minn = 10 * spp;
2769 nom = 0;
2770 denom = 0;
2771 for (i = (long)A; i <= (long)O; i++) {
2772 for (j = (long)A; j <= (long)O; j++) {
2773 if (transvp) {
2774 if (
2775 (
2776 ((i == (long)A) || (i == (long)G))
2777 && ((j == (long)A) || (j == (long)G))
2778 )
2779 || (
2780 ((j == (long)C) || (j == (long)T))
2781 && ((i == (long)C) || (i == (long)T))
2782 )
2783 )
2784 cost = 0;
2785 else
2786 cost = 1;
2787 } else {
2788 if (i == j)
2789 cost = 0;
2790 else
2791 cost = 1;
2792 }
2793 if (subtr1->cumlengths[i] != -1 && (subtr2->cumlengths[j] != -1)) {
2794 if (subtr1->cumlengths[i] + cost + subtr2->cumlengths[j] < minn) {
2795 minn = subtr1->cumlengths[i] + cost + subtr2->cumlengths[j];
2796 nom = 0;
2797 denom = 0;
2798 }
2799 if (subtr1->cumlengths[i] + cost + subtr2->cumlengths[j] == minn) {
2800 nom += subtr1->numreconst[i] * subtr2->numreconst[j] * cost;
2801 denom += subtr1->numreconst[i] * subtr2->numreconst[j];
2802 }
2803 }
2804 }
2805 }
2806 *brlen = (double)nom/(double)denom;
2807 } /* branchlength */
2808
2809
2810 /*void printbranchlengths(node *p)
2811 {
2812 node *q;
2813 long i;
2814
2815 if (p->tip)
2816 return;
2817 q = p->next;
2818 do {
2819 fprintf(outfile, "%6ld ",q->index - spp);
2820 if (q->back->tip) {
2821 for (i = 0; i < nmlngth; i++)
2822 putc(nayme[q->back->index - 1][i], outfile);
2823 } else
2824 fprintf(outfile, "%6ld ", q->back->index - spp);
2825 fprintf(outfile, " %f\n",q->v);
2826 if (q->back)
2827 printbranchlengths(q->back);
2828 q = q->next;
2829 } while (q != p);
2830 } * printbranchlengths */
2831
2832
branchlentrav(node * p,node * root,long sitei,long chars,double * brlen,pointarray treenode)2833 void branchlentrav(node *p, node *root, long sitei, long chars,
2834 double *brlen, pointarray treenode)
2835 {
2836 /* traverses the tree computing tree length at each branch */
2837 node *q;
2838
2839 if (p->tip)
2840 return;
2841 if (p->index == outgrno)
2842 p = p->back;
2843 q = p->next;
2844 do {
2845 if (q->back) {
2846 branchlength(q, q->back, brlen, treenode);
2847 q->v += ((weight[sitei - 1] / 10.0) * (*brlen)/chars);
2848 q->back->v += ((weight[sitei - 1] / 10.0) * (*brlen)/chars);
2849 if (!q->back->tip)
2850 branchlentrav(q->back, root, sitei, chars, brlen, treenode);
2851 }
2852 q = q->next;
2853 } while (q != p);
2854 } /* branchlentrav */
2855
2856
treelength(node * root,long chars,pointarray treenode)2857 void treelength(node *root, long chars, pointarray treenode)
2858 {
2859 /* calls branchlentrav at each site */
2860 long sitei;
2861 double trlen;
2862
2863 initbranchlen(root);
2864 for (sitei = 1; sitei <= endsite; sitei++) {
2865 trlen = 0.0;
2866 initbase(root, sitei);
2867 inittreetrav(root, sitei);
2868 branchlentrav(root, root, sitei, chars, &trlen, treenode);
2869 }
2870 } /* treelength */
2871
2872
coordinates(node * p,long * tipy,double f,long * fartemp)2873 void coordinates(node *p, long *tipy, double f, long *fartemp)
2874 {
2875 /* establishes coordinates of nodes for display without lengths */
2876 node *q, *first, *last;
2877 node *mid1 = NULL, *mid2 = NULL;
2878 long numbranches, numb2;
2879
2880 if (p->tip) {
2881 p->xcoord = 0;
2882 p->ycoord = *tipy;
2883 p->ymin = *tipy;
2884 p->ymax = *tipy;
2885 (*tipy) += down;
2886 return;
2887 }
2888 numbranches = 0;
2889 q = p->next;
2890 do {
2891 coordinates(q->back, tipy, f, fartemp);
2892 numbranches += 1;
2893 q = q->next;
2894 } while (p != q);
2895 first = p->next->back;
2896 q = p->next;
2897 while (q->next != p)
2898 q = q->next;
2899 last = q->back;
2900 numb2 = 1;
2901 q = p->next;
2902 while (q != p) {
2903 if (numb2 == (long)(numbranches + 1)/2)
2904 mid1 = q->back;
2905 if (numb2 == (long)(numbranches/2 + 1))
2906 mid2 = q->back;
2907 numb2 += 1;
2908 q = q->next;
2909 }
2910 p->xcoord = (long)((double)(last->ymax - first->ymin) * f);
2911 p->ycoord = (long)((mid1->ycoord + mid2->ycoord) / 2);
2912 p->ymin = first->ymin;
2913 p->ymax = last->ymax;
2914 if (p->xcoord > *fartemp)
2915 *fartemp = p->xcoord;
2916 } /* coordinates */
2917
2918
2919 /*void drawline(long i, double scale, node *root)
2920 {
2921 * draws one row of the tree diagram by moving up tree *
2922 node *p, *q, *r, *first =NULL, *last =NULL;
2923 long n, j;
2924 boolean extra, done, noplus;
2925
2926 p = root;
2927 q = root;
2928 extra = false;
2929 noplus = false;
2930 if (i == (long)p->ycoord && p == root) {
2931 if (p->index - spp >= 10)
2932 fprintf(outfile, " %2ld", p->index - spp);
2933 else
2934 fprintf(outfile, " %ld", p->index - spp);
2935 extra = true;
2936 noplus = true;
2937 } else
2938 fprintf(outfile, " ");
2939 do {
2940 if (!p->tip) {
2941 r = p->next;
2942 done = false;
2943 do {
2944 if (i >= r->back->ymin && i <= r->back->ymax) {
2945 q = r->back;
2946 done = true;
2947 }
2948 r = r->next;
2949 } while (!(done || r == p));
2950 first = p->next->back;
2951 r = p->next;
2952 while (r->next != p)
2953 r = r->next;
2954 last = r->back;
2955 }
2956 done = (p == q);
2957 n = (long)(scale * (p->xcoord - q->xcoord) + 0.5);
2958 if (n < 3 && !q->tip)
2959 n = 3;
2960 if (extra) {
2961 n--;
2962 extra = false;
2963 }
2964 if ((long)q->ycoord == i && !done) {
2965 if (noplus) {
2966 putc('-', outfile);
2967 noplus = false;
2968 }
2969 else
2970 putc('+', outfile);
2971 if (!q->tip) {
2972 for (j = 1; j <= n - 2; j++)
2973 putc('-', outfile);
2974 if (q->index - spp >= 10)
2975 fprintf(outfile, "%2ld", q->index - spp);
2976 else
2977 fprintf(outfile, "-%ld", q->index - spp);
2978 extra = true;
2979 noplus = true;
2980 } else {
2981 for (j = 1; j < n; j++)
2982 putc('-', outfile);
2983 }
2984 } else if (!p->tip) {
2985 if ((long)last->ycoord > i && (long)first->ycoord < i
2986 && i != (long)p->ycoord) {
2987 putc('!', outfile);
2988 for (j = 1; j < n; j++)
2989 putc(' ', outfile);
2990 } else {
2991 for (j = 1; j <= n; j++)
2992 putc(' ', outfile);
2993 }
2994 noplus = false;
2995 } else {
2996 for (j = 1; j <= n; j++)
2997 putc(' ', outfile);
2998 noplus = false;
2999 }
3000 if (p != q)
3001 p = q;
3002 } while (!done);
3003 if ((long)p->ycoord == i && p->tip) {
3004 for (j = 0; j < nmlngth; j++)
3005 putc(nayme[p->index - 1][j], outfile);
3006 }
3007 putc('\n', outfile);
3008 } * drawline *
3009
3010
3011 void printree(node *root, double f)
3012 {
3013 * prints out diagram of the tree *
3014 * used in dnacomp, dnapars, & dnapenny *
3015 long i, tipy, dummy;
3016 double scale;
3017
3018 putc('\n', outfile);
3019 if (!treeprint)
3020 return;
3021 putc('\n', outfile);
3022 tipy = 1;
3023 dummy = 0;
3024 coordinates(root, &tipy, f, &dummy);
3025 scale = 1.5;
3026 putc('\n', outfile);
3027 for (i = 1; i <= (tipy - down); i++)
3028 drawline(i, scale, root);
3029 fprintf(outfile, "\n remember:");
3030 if (outgropt)
3031 fprintf(outfile, " (although rooted by outgroup)");
3032 fprintf(outfile, " this is an unrooted tree!\n\n");
3033 } * printree */
3034
3035
writesteps(long chars,boolean weights,steptr oldweight,node * root)3036 void writesteps(long chars, boolean weights, steptr oldweight, node *root)
3037 {
3038 /* used in dnacomp, dnapars, & dnapenny */
3039 long i, j, k, l;
3040
3041 putc('\n', outfile);
3042 if (weights)
3043 fprintf(outfile, "weighted ");
3044 fprintf(outfile, "steps in each site:\n");
3045 fprintf(outfile, " ");
3046 for (i = 0; i <= 9; i++)
3047 fprintf(outfile, "%4ld", i);
3048 fprintf(outfile, "\n *------------------------------------");
3049 fprintf(outfile, "-----\n");
3050 for (i = 0; i <= (chars / 10); i++) {
3051 fprintf(outfile, "%5ld", i * 10);
3052 putc('|', outfile);
3053 for (j = 0; j <= 9; j++) {
3054 k = i * 10 + j;
3055 if (k == 0 || k > chars)
3056 fprintf(outfile, " ");
3057 else {
3058 l = location[ally[k - 1] - 1];
3059 if (oldweight[k - 1] > 0)
3060 fprintf(outfile, "%4ld",
3061 oldweight[k - 1] *
3062 (root->numsteps[l - 1] / weight[l - 1]));
3063 else
3064 fprintf(outfile, " 0");
3065 }
3066 }
3067 putc('\n', outfile);
3068 }
3069 } /* writesteps */
3070
3071
treeout(node * p,long nextree,long * col,node * root)3072 void treeout(node *p, long nextree, long *col, node *root)
3073 {
3074 /* write out file with representation of final tree */
3075 /* used in dnacomp, dnamove, dnapars, & dnapenny */
3076 node *q;
3077 long i, n;
3078 Char c;
3079
3080 if (p->tip) {
3081 n = strlen(nayme[p->index - 1]);
3082 /*for (i = 1; i <= nmlngth; i++) {
3083 if (nayme[p->index - 1][i - 1] != ' ')
3084 n = i;
3085 }*/
3086 for (i = 0; i < n; i++) {
3087 c = nayme[p->index - 1][i];
3088 /*if (c == ' ')
3089 c = '_';*/
3090 putc(c, outtree);
3091 }
3092 //*col += n;
3093 } else {
3094 putc('(', outtree);
3095 //(*col)++;
3096 q = p->next;
3097 while (q != p) {
3098 treeout(q->back, nextree, col, root);
3099 q = q->next;
3100 if (q == p)
3101 break;
3102 putc(',', outtree);
3103 /*(*col)++;
3104 if (*col > 60) {
3105 putc('\n', outtree);
3106 *col = 0;
3107 }*/
3108 }
3109 putc(')', outtree);
3110 //(*col)++;
3111 }
3112 if (p != root)
3113 return;
3114 if (nextree > 2)
3115 fprintf(outtree, "[%6.4f];\n", 1.0 / (nextree - 1));
3116 else
3117 fprintf(outtree, ";\n");
3118 } /* treeout */
3119
3120
treeout3(node * p,long nextree,long * col,node * root)3121 void treeout3(node *p, long nextree, long *col, node *root)
3122 {
3123 /* write out file with representation of final tree */
3124 /* used in dnapars -- writes branch lengths */
3125 node *q;
3126 long i, n, w;
3127 double x;
3128 Char c;
3129
3130 if (p->tip) {
3131 n = 0;
3132 /*for (i = 1; i <= nmlngth; i++) {
3133 if (nayme[p->index - 1][i - 1] != ' ')
3134 n = i;
3135 }*/
3136 n = strlen(nayme[p->index - 1]);
3137 for (i = 0; i < n; i++) {
3138 c = nayme[p->index - 1][i];
3139 /*if (c == ' ')
3140 c = '_';*/
3141 putc(c, outtree);
3142 }
3143 //*col += n;
3144 } else {
3145 putc('(', outtree);
3146 //(*col)++;
3147 q = p->next;
3148 while (q != p) {
3149 treeout3(q->back, nextree, col, root);
3150 q = q->next;
3151 if (q == p)
3152 break;
3153 putc(',', outtree);
3154 /*(*col)++;
3155 if (*col > 60) {
3156 putc('\n', outtree);
3157 *col = 0;
3158 }*/
3159 }
3160 putc(')', outtree);
3161 //(*col)++;
3162 }
3163 x = p->v;
3164 if (x > 0.0)
3165 w = (long)(0.43429448222 * log(x));
3166 else if (x == 0.0)
3167 w = 0;
3168 else
3169 w = (long)(0.43429448222 * log(-x)) + 1;
3170 if (w < 0)
3171 w = 0;
3172 if (p != root) {
3173 fprintf(outtree, ":%*.5f", (int)(w + 7), x);
3174 //*col += w + 8;
3175 }
3176 if (p != root)
3177 return;
3178 if (nextree > 2)
3179 fprintf(outtree, "[%6.4f];\n", 1.0 / (nextree - 1));
3180 else
3181 fprintf(outtree, ";\n");
3182 } /* treeout3 */
3183
3184
3185 /* FIXME curtree should probably be passed by reference *
3186 void drawline2(long i, double scale, tree curtree)
3187 {
3188 fdrawline2(outfile, i, scale, &curtree);
3189 }*
3190
3191 void fdrawline2(FILE *fp, long i, double scale, tree *curtree)
3192 {
3193 * draws one row of the tree diagram by moving up tree *
3194 * used in dnaml, proml, & restml *
3195 node *p, *q;
3196 long n, j;
3197 boolean extra;
3198 node *r, *first =NULL, *last =NULL;
3199 boolean done;
3200
3201 p = curtree->start;
3202 q = curtree->start;
3203 extra = false;
3204 if (i == (long)p->ycoord && p == curtree->start) {
3205 if (p->index - spp >= 10)
3206 fprintf(fp, " %2ld", p->index - spp);
3207 else
3208 fprintf(fp, " %ld", p->index - spp);
3209 extra = true;
3210 } else
3211 fprintf(fp, " ");
3212 do {
3213 if (!p->tip) {
3214 r = p->next;
3215 done = false;
3216 do {
3217 if (i >= r->back->ymin && i <= r->back->ymax) {
3218 q = r->back;
3219 done = true;
3220 }
3221 r = r->next;
3222 } while (!(done || (p != curtree->start && r == p) ||
3223 (p == curtree->start && r == p->next)));
3224 first = p->next->back;
3225 r = p;
3226 while (r->next != p)
3227 r = r->next;
3228 last = r->back;
3229 if (p == curtree->start)
3230 last = p->back;
3231 }
3232 done = (p->tip || p == q);
3233 n = (long)(scale * (q->xcoord - p->xcoord) + 0.5);
3234 if (n < 3 && !q->tip)
3235 n = 3;
3236 if (extra) {
3237 n--;
3238 extra = false;
3239 }
3240 if ((long)q->ycoord == i && !done) {
3241 if ((long)p->ycoord != (long)q->ycoord)
3242 putc('+', fp);
3243 else
3244 putc('-', fp);
3245 if (!q->tip) {
3246 for (j = 1; j <= n - 2; j++)
3247 putc('-', fp);
3248 if (q->index - spp >= 10)
3249 fprintf(fp, "%2ld", q->index - spp);
3250 else
3251 fprintf(fp, "-%ld", q->index - spp);
3252 extra = true;
3253 } else {
3254 for (j = 1; j < n; j++)
3255 putc('-', fp);
3256 }
3257 } else if (!p->tip) {
3258 if ((long)last->ycoord > i && (long)first->ycoord < i &&
3259 (i != (long)p->ycoord || p == curtree->start)) {
3260 putc('|', fp);
3261 for (j = 1; j < n; j++)
3262 putc(' ', fp);
3263 } else {
3264 for (j = 1; j <= n; j++)
3265 putc(' ', fp);
3266 }
3267 } else {
3268 for (j = 1; j <= n; j++)
3269 putc(' ', fp);
3270 }
3271 if (q != p)
3272 p = q;
3273 } while (!done);
3274 if ((long)p->ycoord == i && p->tip) {
3275 for (j = 0; j < nmlngth; j++)
3276 putc(nayme[p->index-1][j], fp);
3277 }
3278 putc('\n', fp);
3279 } * drawline2 *
3280
3281
3282 void drawline3(long i, double scale, node *start)
3283 {
3284 * draws one row of the tree diagram by moving up tree *
3285 * used in dnapars *
3286 node *p, *q;
3287 long n, j;
3288 boolean extra;
3289 node *r, *first =NULL, *last =NULL;
3290 boolean done;
3291
3292 p = start;
3293 q = start;
3294 extra = false;
3295 if (i == (long)p->ycoord) {
3296 if (p->index - spp >= 10)
3297 fprintf(outfile, " %2ld", p->index - spp);
3298 else
3299 fprintf(outfile, " %ld", p->index - spp);
3300 extra = true;
3301 } else
3302 fprintf(outfile, " ");
3303 do {
3304 if (!p->tip) {
3305 r = p->next;
3306 done = false;
3307 do {
3308 if (i >= r->back->ymin && i <= r->back->ymax) {
3309 q = r->back;
3310 done = true;
3311 }
3312 r = r->next;
3313 } while (!(done || (r == p)));
3314 first = p->next->back;
3315 r = p;
3316 while (r->next != p)
3317 r = r->next;
3318 last = r->back;
3319 }
3320 done = (p->tip || p == q);
3321 n = (long)(scale * (q->xcoord - p->xcoord) + 0.5);
3322 if (n < 3 && !q->tip)
3323 n = 3;
3324 if (extra) {
3325 n--;
3326 extra = false;
3327 }
3328 if ((long)q->ycoord == i && !done) {
3329 if ((long)p->ycoord != (long)q->ycoord)
3330 putc('+', outfile);
3331 else
3332 putc('-', outfile);
3333 if (!q->tip) {
3334 for (j = 1; j <= n - 2; j++)
3335 putc('-', outfile);
3336 if (q->index - spp >= 10)
3337 fprintf(outfile, "%2ld", q->index - spp);
3338 else
3339 fprintf(outfile, "-%ld", q->index - spp);
3340 extra = true;
3341 } else {
3342 for (j = 1; j < n; j++)
3343 putc('-', outfile);
3344 }
3345 } else if (!p->tip) {
3346 if ((long)last->ycoord > i && (long)first->ycoord < i &&
3347 (i != (long)p->ycoord || p == start)) {
3348 putc('|', outfile);
3349 for (j = 1; j < n; j++)
3350 putc(' ', outfile);
3351 } else {
3352 for (j = 1; j <= n; j++)
3353 putc(' ', outfile);
3354 }
3355 } else {
3356 for (j = 1; j <= n; j++)
3357 putc(' ', outfile);
3358 }
3359 if (q != p)
3360 p = q;
3361 } while (!done);
3362 if ((long)p->ycoord == i && p->tip) {
3363 for (j = 0; j < nmlngth; j++)
3364 putc(nayme[p->index-1][j], outfile);
3365 }
3366 putc('\n', outfile);
3367 } * drawline3 */
3368
3369
copynode(node * c,node * d,long categs)3370 void copynode(node *c, node *d, long categs)
3371 {
3372 long i, j;
3373
3374 for (i = 0; i < endsite; i++)
3375 for (j = 0; j < categs; j++)
3376 memcpy(d->x[i][j], c->x[i][j], sizeof(sitelike));
3377 memcpy(d->underflows,c->underflows,sizeof(double) * endsite);
3378 d->tyme = c->tyme;
3379 d->v = c->v;
3380 d->xcoord = c->xcoord;
3381 d->ycoord = c->ycoord;
3382 d->ymin = c->ymin;
3383 d->ymax = c->ymax;
3384 d->iter = c->iter; /* iter used in dnaml only */
3385 d->haslength = c->haslength; /* haslength used in dnamlk only */
3386 d->initialized = c->initialized; /* initialized used in dnamlk only */
3387 } /* copynode */
3388
3389
prot_copynode(node * c,node * d,long categs)3390 void prot_copynode(node *c, node *d, long categs)
3391 {
3392 /* a version of copynode for proml */
3393 long i, j;
3394
3395 for (i = 0; i < endsite; i++)
3396 for (j = 0; j < categs; j++)
3397 memcpy(d->protx[i][j], c->protx[i][j], sizeof(psitelike));
3398 memcpy(d->underflows,c->underflows,sizeof(double) * endsite);
3399 d->tyme = c->tyme;
3400 d->v = c->v;
3401 d->xcoord = c->xcoord;
3402 d->ycoord = c->ycoord;
3403 d->ymin = c->ymin;
3404 d->ymax = c->ymax;
3405 d->iter = c->iter; /* iter used in dnaml only */
3406 d->haslength = c->haslength; /* haslength used in dnamlk only */
3407 d->initialized = c->initialized; /* initialized used in dnamlk only */
3408 } /* prot_copynode */
3409
3410
copy_(tree * a,tree * b,long nonodes,long categs)3411 void copy_(tree *a, tree *b, long nonodes, long categs)
3412 {
3413 /* used in dnamlk */
3414 long i;
3415 node *p, *q, *r, *s, *t;
3416
3417 for (i = 0; i < spp; i++) {
3418 copynode(a->nodep[i], b->nodep[i], categs);
3419 if (a->nodep[i]->back) {
3420 if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1])
3421 b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1];
3422 else if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1]->next)
3423 b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next;
3424 else
3425 b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next->next;
3426 }
3427 else b->nodep[i]->back = NULL;
3428 }
3429 for (i = spp; i < nonodes; i++) {
3430 if (a->nodep[i]) {
3431 p = a->nodep[i];
3432 q = b->nodep[i];
3433 r = p;
3434 do {
3435 copynode(p, q, categs);
3436 if (p->back) {
3437 s = a->nodep[p->back->index - 1];
3438 t = b->nodep[p->back->index - 1];
3439 if (s->tip) {
3440 if(p->back == s)
3441 q->back = t;
3442 } else {
3443 do {
3444 if (p->back == s)
3445 q->back = t;
3446 s = s->next;
3447 t = t->next;
3448 } while (s != a->nodep[p->back->index - 1]);
3449 }
3450 }
3451 else
3452 q->back = NULL;
3453 p = p->next;
3454 q = q->next;
3455 } while (p != r);
3456 }
3457 }
3458 b->likelihood = a->likelihood;
3459 b->start = a->start; /* start used in dnaml only */
3460 b->root = a->root; /* root used in dnamlk only */
3461 } /* copy_ */
3462
3463
prot_copy_(tree * a,tree * b,long nonodes,long categs)3464 void prot_copy_(tree *a, tree *b, long nonodes, long categs)
3465 {
3466 /* used in promlk */
3467 /* identical to copy_() except for calls to prot_copynode rather */
3468 /* than copynode. */
3469 long i;
3470 node *p, *q, *r, *s, *t;
3471
3472 for (i = 0; i < spp; i++) {
3473 prot_copynode(a->nodep[i], b->nodep[i], categs);
3474 if (a->nodep[i]->back) {
3475 if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1])
3476 b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1];
3477 else if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1]->next
3478 )
3479 b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next;
3480 else
3481 b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next->next;
3482 }
3483 else b->nodep[i]->back = NULL;
3484 }
3485 for (i = spp; i < nonodes; i++) {
3486 if (a->nodep[i]) {
3487 p = a->nodep[i];
3488 q = b->nodep[i];
3489 r = p;
3490 do {
3491 prot_copynode(p, q, categs);
3492 if (p->back) {
3493 s = a->nodep[p->back->index - 1];
3494 t = b->nodep[p->back->index - 1];
3495 if (s->tip)
3496 {
3497 if(p->back == s)
3498 q->back = t;
3499 } else {
3500 do {
3501 if (p->back == s)
3502 q->back = t;
3503 s = s->next;
3504 t = t->next;
3505 } while (s != a->nodep[p->back->index - 1]);
3506 }
3507 }
3508 else
3509 q->back = NULL;
3510 p = p->next;
3511 q = q->next;
3512 } while (p != r);
3513 }
3514 }
3515 b->likelihood = a->likelihood;
3516 b->start = a->start; /* start used in dnaml only */
3517 b->root = a->root; /* root used in dnamlk only */
3518 } /* prot_copy_ */
3519
3520
standev(long chars,long numtrees,long minwhich,double minsteps,double * nsteps,long ** fsteps,longer seed)3521 void standev(long chars, long numtrees, long minwhich, double minsteps,
3522 double *nsteps, long **fsteps, longer seed)
3523 { /* do paired sites test (KHT or SH test) on user-defined trees */
3524 /* used in dnapars & protpars */
3525 long i, j, k;
3526 double wt, sumw, sum, sum2, sd;
3527 double temp;
3528 double **covar, *P, *f, *r;
3529
3530 #define SAMPLES 1000
3531 if (numtrees == 2) {
3532 fprintf(outfile, "Kishino-Hasegawa-Templeton test\n\n");
3533 fprintf(outfile, "Tree Steps Diff Steps Its S.D.");
3534 fprintf(outfile, " Significantly worse?\n\n");
3535 which = 1;
3536 while (which <= numtrees) {
3537 fprintf(outfile, "%3ld%10.1f", which, nsteps[which - 1] / 10);
3538 if (minwhich == which)
3539 fprintf(outfile, " <------ best\n");
3540 else {
3541 sumw = 0.0;
3542 sum = 0.0;
3543 sum2 = 0.0;
3544 for (i = 0; i < endsite; i++) {
3545 if (weight[i] > 0) {
3546 wt = weight[i] / 10.0;
3547 sumw += wt;
3548 temp = (fsteps[which - 1][i] - fsteps[minwhich - 1][i]) / 10.0;
3549 sum += wt * temp;
3550 sum2 += wt * temp * temp;
3551 }
3552 }
3553 sd = sqrt(sumw / (sumw - 1.0) * (sum2 - sum * sum / sumw));
3554 fprintf(outfile, "%10.1f%12.4f",
3555 (nsteps[which - 1] - minsteps) / 10, sd);
3556 if ((sum > 0.0) && (sum > 1.95996 * sd))
3557 fprintf(outfile, " Yes\n");
3558 else
3559 fprintf(outfile, " No\n");
3560 }
3561 which++;
3562 }
3563 fprintf(outfile, "\n\n");
3564 } else { /* Shimodaira-Hasegawa test using normal approximation */
3565 if(numtrees > MAXSHIMOTREES){
3566 fprintf(outfile, "Shimodaira-Hasegawa test on first %d of %ld trees\n\n"
3567 , MAXSHIMOTREES, numtrees);
3568 numtrees = MAXSHIMOTREES;
3569 } else {
3570 fprintf(outfile, "Shimodaira-Hasegawa test\n\n");
3571 }
3572 covar = (double **)Malloc(numtrees*sizeof(double *));
3573 sumw = 0.0;
3574 for (i = 0; i < endsite; i++)
3575 sumw += weight[i] / 10.0;
3576 for (i = 0; i < numtrees; i++)
3577 covar[i] = (double *)Malloc(numtrees*sizeof(double));
3578 for (i = 0; i < numtrees; i++) { /* compute covariances of trees */
3579 sum = nsteps[i]/(10.0*sumw);
3580 for (j = 0; j <=i; j++) {
3581 sum2 = nsteps[j]/(10.0*sumw);
3582 temp = 0.0;
3583 for (k = 0; k < endsite; k++) {
3584 if (weight[k] > 0) {
3585 wt = weight[k]/10.0;
3586 temp = temp + wt*(fsteps[i][k]/10.0-sum)
3587 *(fsteps[j][k]/10.0-sum2);
3588 }
3589 }
3590 covar[i][j] = temp;
3591 if (i != j)
3592 covar[j][i] = temp;
3593 }
3594 }
3595 for (i = 0; i < numtrees; i++) { /* in-place Cholesky decomposition
3596 of trees x trees covariance matrix */
3597 sum = 0.0;
3598 for (j = 0; j <= i-1; j++)
3599 sum = sum + covar[i][j] * covar[i][j];
3600 if (covar[i][i] <= sum)
3601 temp = 0.0;
3602 else
3603 temp = sqrt(covar[i][i] - sum);
3604 covar[i][i] = temp;
3605 for (j = i+1; j < numtrees; j++) {
3606 sum = 0.0;
3607 for (k = 0; k < i; k++)
3608 sum = sum + covar[i][k] * covar[j][k];
3609 if (fabs(temp) < 1.0E-12)
3610 covar[j][i] = 0.0;
3611 else
3612 covar[j][i] = (covar[j][i] - sum)/temp;
3613 }
3614 }
3615 f = (double *)Malloc(numtrees*sizeof(double)); /* resampled sums */
3616 P = (double *)Malloc(numtrees*sizeof(double)); /* vector of P's of trees */
3617 r = (double *)Malloc(numtrees*sizeof(double)); /* store Normal variates */
3618 for (i = 0; i < numtrees; i++)
3619 P[i] = 0.0;
3620 sum2 = nsteps[0]/10.0; /* sum2 will be smallest # of steps */
3621 for (i = 1; i < numtrees; i++)
3622 if (sum2 > nsteps[i]/10.0)
3623 sum2 = nsteps[i]/10.0;
3624 for (i = 1; i <= SAMPLES; i++) { /* loop over resampled trees */
3625 for (j = 0; j < numtrees; j++) /* draw Normal variates */
3626 r[j] = normrand(seed);
3627 for (j = 0; j < numtrees; j++) { /* compute vectors */
3628 sum = 0.0;
3629 for (k = 0; k <= j; k++)
3630 sum += covar[j][k]*r[k];
3631 f[j] = sum;
3632 }
3633 sum = f[1];
3634 for (j = 1; j < numtrees; j++) /* get min of vector */
3635 if (f[j] < sum)
3636 sum = f[j];
3637 for (j = 0; j < numtrees; j++) /* accumulate P's */
3638 if (nsteps[j]/10.0-sum2 <= f[j] - sum)
3639 P[j] += 1.0/SAMPLES;
3640 }
3641 fprintf(outfile, "Tree Steps Diff Steps P value");
3642 fprintf(outfile, " Significantly worse?\n\n");
3643 for (i = 0; i < numtrees; i++) {
3644 fprintf(outfile, "%3ld%10.1f", i+1, nsteps[i]/10);
3645 if ((minwhich-1) == i)
3646 fprintf(outfile, " <------ best\n");
3647 else {
3648 fprintf(outfile, " %9.1f %10.3f", nsteps[i]/10.0-sum2, P[i]);
3649 if (P[i] < 0.05)
3650 fprintf(outfile, " Yes\n");
3651 else
3652 fprintf(outfile, " No\n");
3653 }
3654 }
3655 fprintf(outfile, "\n");
3656 free(P); /* free the variables we Malloc'ed */
3657 free(f);
3658 free(r);
3659 for (i = 0; i < numtrees; i++)
3660 free(covar[i]);
3661 free(covar);
3662 }
3663 } /* standev */
3664
3665
standev2(long numtrees,long maxwhich,long a,long b,double maxlogl,double * l0gl,double ** l0gf,steptr aliasweight,longer seed)3666 void standev2(long numtrees, long maxwhich, long a, long b, double maxlogl,
3667 double *l0gl, double **l0gf, steptr aliasweight, longer seed)
3668 { /* do paired sites test (KHT or SH) for user-defined trees */
3669 /* used in dnaml, dnamlk, proml, promlk, and restml */
3670 double **covar, *P, *f, *r;
3671 long i, j, k;
3672 double wt, sumw, sum, sum2, sd;
3673 double temp;
3674
3675 #define SAMPLES 1000
3676 if (numtrees == 2) {
3677 fprintf(outfile, "Kishino-Hasegawa-Templeton test\n\n");
3678 fprintf(outfile, "Tree logL Diff logL Its S.D.");
3679 fprintf(outfile, " Significantly worse?\n\n");
3680 which = 1;
3681 while (which <= numtrees) {
3682 fprintf(outfile, "%3ld %9.1f", which, l0gl[which - 1]);
3683 if (maxwhich == which)
3684 fprintf(outfile, " <------ best\n");
3685 else {
3686 sumw = 0.0;
3687 sum = 0.0;
3688 sum2 = 0.0;
3689 for (i = a; i <= b; i++) {
3690 if (aliasweight[i] > 0) {
3691 wt = aliasweight[i];
3692 sumw += wt;
3693 temp = l0gf[which - 1][i] - l0gf[maxwhich - 1][i];
3694 sum += temp * wt;
3695 sum2 += wt * temp * temp;
3696 }
3697 }
3698 temp = sum / sumw;
3699 sd = sqrt(sumw / (sumw - 1.0) * (sum2 - sum * sum / sumw ));
3700 fprintf(outfile, "%10.1f %11.4f", (l0gl[which - 1])-maxlogl, sd);
3701 if ((sum < 0.0) && ((-sum) > 1.95996 * sd))
3702 fprintf(outfile, " Yes\n");
3703 else
3704 fprintf(outfile, " No\n");
3705 }
3706 which++;
3707 }
3708 fprintf(outfile, "\n\n");
3709 } else { /* Shimodaira-Hasegawa test using normal approximation */
3710 if(numtrees > MAXSHIMOTREES){
3711 fprintf(outfile, "Shimodaira-Hasegawa test on first %d of %ld trees\n\n"
3712 , MAXSHIMOTREES, numtrees);
3713 numtrees = MAXSHIMOTREES;
3714 } else {
3715 fprintf(outfile, "Shimodaira-Hasegawa test\n\n");
3716 }
3717 covar = (double **)Malloc(numtrees*sizeof(double *));
3718 sumw = 0.0;
3719 for (i = a; i <= b; i++)
3720 sumw += aliasweight[i];
3721 for (i = 0; i < numtrees; i++)
3722 covar[i] = (double *)Malloc(numtrees*sizeof(double));
3723 for (i = 0; i < numtrees; i++) { /* compute covariances of trees */
3724 sum = l0gl[i]/sumw;
3725 for (j = 0; j <=i; j++) {
3726 sum2 = l0gl[j]/sumw;
3727 temp = 0.0;
3728 for (k = a; k <= b ; k++) {
3729 if (aliasweight[k] > 0) {
3730 wt = aliasweight[k];
3731 temp = temp + wt*(l0gf[i][k]-sum)
3732 *(l0gf[j][k]-sum2);
3733 }
3734 }
3735 covar[i][j] = temp;
3736 if (i != j)
3737 covar[j][i] = temp;
3738 }
3739 }
3740 for (i = 0; i < numtrees; i++) { /* in-place Cholesky decomposition
3741 of trees x trees covariance matrix */
3742 sum = 0.0;
3743 for (j = 0; j <= i-1; j++)
3744 sum = sum + covar[i][j] * covar[i][j];
3745 if (covar[i][i] <= sum)
3746 temp = 0.0;
3747 else
3748 temp = sqrt(covar[i][i] - sum);
3749 covar[i][i] = temp;
3750 for (j = i+1; j < numtrees; j++) {
3751 sum = 0.0;
3752 for (k = 0; k < i; k++)
3753 sum = sum + covar[i][k] * covar[j][k];
3754 if (fabs(temp) < 1.0E-12)
3755 covar[j][i] = 0.0;
3756 else
3757 covar[j][i] = (covar[j][i] - sum)/temp;
3758 }
3759 }
3760 f = (double *)Malloc(numtrees*sizeof(double)); /* resampled likelihoods */
3761 P = (double *)Malloc(numtrees*sizeof(double)); /* vector of P's of trees */
3762 r = (double *)Malloc(numtrees*sizeof(double)); /* store Normal variates */
3763 for (i = 0; i < numtrees; i++)
3764 P[i] = 0.0;
3765 for (i = 1; i <= SAMPLES; i++) { /* loop over resampled trees */
3766 for (j = 0; j < numtrees; j++) /* draw Normal variates */
3767 r[j] = normrand(seed);
3768 for (j = 0; j < numtrees; j++) { /* compute vectors */
3769 sum = 0.0;
3770 for (k = 0; k <= j; k++)
3771 sum += covar[j][k]*r[k];
3772 f[j] = sum;
3773 }
3774 sum = f[1];
3775 for (j = 1; j < numtrees; j++) /* get max of vector */
3776 if (f[j] > sum)
3777 sum = f[j];
3778 for (j = 0; j < numtrees; j++) /* accumulate P's */
3779 if (maxlogl-l0gl[j] <= sum-f[j])
3780 P[j] += 1.0/SAMPLES;
3781 }
3782 fprintf(outfile, "Tree logL Diff logL P value");
3783 fprintf(outfile, " Significantly worse?\n\n");
3784 for (i = 0; i < numtrees; i++) {
3785 fprintf(outfile, "%3ld%10.1f", i+1, l0gl[i]);
3786 if ((maxwhich-1) == i)
3787 fprintf(outfile, " <------ best\n");
3788 else {
3789 fprintf(outfile, " %9.1f %10.3f", l0gl[i]-maxlogl, P[i]);
3790 if (P[i] < 0.05)
3791 fprintf(outfile, " Yes\n");
3792 else
3793 fprintf(outfile, " No\n");
3794 }
3795 }
3796 fprintf(outfile, "\n");
3797 free(P); /* free the variables we Malloc'ed */
3798 free(f);
3799 free(r);
3800 for (i = 0; i < numtrees; i++)
3801 free(covar[i]);
3802 free(covar);
3803 }
3804 } /* standev */
3805
3806
freetip(node * anode)3807 void freetip(node *anode)
3808 {
3809 /* used in dnacomp, dnapars, & dnapenny */
3810
3811 free(anode->numsteps);
3812 free(anode->oldnumsteps);
3813 free(anode->base);
3814 free(anode->oldbase);
3815 } /* freetip */
3816
3817
freenontip(node * anode)3818 void freenontip(node *anode)
3819 {
3820 /* used in dnacomp, dnapars, & dnapenny */
3821
3822 free(anode->numsteps);
3823 free(anode->oldnumsteps);
3824 free(anode->base);
3825 free(anode->oldbase);
3826 free(anode->numnuc);
3827 } /* freenontip */
3828
3829
freenodes(long nonodes,pointarray treenode)3830 void freenodes(long nonodes, pointarray treenode)
3831 {
3832 /* used in dnacomp, dnapars, & dnapenny */
3833 long i;
3834 node *p;
3835
3836 for (i = 0; i < spp; i++)
3837 freetip(treenode[i]);
3838 for (i = spp; i < nonodes; i++) {
3839 if (treenode[i] != NULL) {
3840 p = treenode[i]->next;
3841 do {
3842 freenontip(p);
3843 p = p->next;
3844 } while (p != treenode[i]);
3845 freenontip(p);
3846 }
3847 }
3848 } /* freenodes */
3849
3850
freenode(node ** anode)3851 void freenode(node **anode)
3852 {
3853 /* used in dnacomp, dnapars, & dnapenny */
3854
3855 freenontip(*anode);
3856 free(*anode);
3857 } /* freenode */
3858
3859
freetree(long nonodes,pointarray treenode)3860 void freetree(long nonodes, pointarray treenode)
3861 {
3862 /* used in dnacomp, dnapars, & dnapenny */
3863 long i;
3864 node *p, *q;
3865
3866 for (i = 0; i < spp; i++)
3867 free(treenode[i]);
3868 for (i = spp; i < nonodes; i++) {
3869 if (treenode[i] != NULL) {
3870 p = treenode[i]->next;
3871 do {
3872 q = p->next;
3873 free(p);
3874 p = q;
3875 } while (p != treenode[i]);
3876 free(p);
3877 }
3878 }
3879 free(treenode);
3880 } /* freetree */
3881
3882
prot_freex_notip(long nonodes,pointarray treenode)3883 void prot_freex_notip(long nonodes, pointarray treenode)
3884 {
3885 /* used in proml */
3886 long i, j;
3887 node *p;
3888
3889 for (i = spp; i < nonodes; i++) {
3890 p = treenode[i];
3891 if ( p == NULL ) continue;
3892 do {
3893 for (j = 0; j < endsite; j++){
3894 free(p->protx[j]);
3895 p->protx[j] = NULL;
3896 }
3897 free(p->underflows);
3898 p->underflows = NULL;
3899 free(p->protx);
3900 p->protx = NULL;
3901 p = p->next;
3902 } while (p != treenode[i]);
3903 }
3904 } /* prot_freex_notip */
3905
3906
prot_freex(long nonodes,pointarray treenode)3907 void prot_freex(long nonodes, pointarray treenode)
3908 {
3909 /* used in proml */
3910 long i, j;
3911 node *p;
3912
3913 for (i = 0; i < spp; i++) {
3914 for (j = 0; j < endsite; j++)
3915 free(treenode[i]->protx[j]);
3916 free(treenode[i]->protx);
3917 free(treenode[i]->underflows);
3918 }
3919 for (i = spp; i < nonodes; i++) {
3920 p = treenode[i];
3921 do {
3922 for (j = 0; j < endsite; j++)
3923 free(p->protx[j]);
3924 free(p->protx);
3925 free(p->underflows);
3926 p = p->next;
3927 } while (p != treenode[i]);
3928 }
3929 } /* prot_freex */
3930
3931
freex_notip(long nonodes,pointarray treenode)3932 void freex_notip(long nonodes, pointarray treenode)
3933 {
3934 /* used in dnaml & dnamlk */
3935 long i, j;
3936 node *p;
3937
3938 for (i = spp; i < nonodes; i++) {
3939 p = treenode[i];
3940 if ( p == NULL ) continue;
3941 do {
3942 for (j = 0; j < endsite; j++)
3943 free(p->x[j]);
3944 free(p->underflows);
3945 free(p->x);
3946 p = p->next;
3947 } while (p != treenode[i]);
3948 }
3949 } /* freex_notip */
3950
3951
freex(long nonodes,pointarray treenode)3952 void freex(long nonodes, pointarray treenode)
3953 {
3954 /* used in dnaml & dnamlk */
3955 long i, j;
3956 node *p;
3957
3958 for (i = 0; i < spp; i++) {
3959 for (j = 0; j < endsite; j++)
3960 free(treenode[i]->x[j]);
3961 free(treenode[i]->x);
3962 free(treenode[i]->underflows);
3963 }
3964 for (i = spp; i < nonodes; i++) {
3965 if(treenode[i]){
3966 p = treenode[i];
3967 do {
3968 for (j = 0; j < endsite; j++)
3969 free(p->x[j]);
3970 free(p->x);
3971 free(p->underflows);
3972 p = p->next;
3973 } while (p != treenode[i]);
3974 }
3975 }
3976 } /* freex */
3977
3978
3979
freegarbage(gbases ** garbage)3980 void freegarbage(gbases **garbage)
3981 {
3982 /* used in dnacomp, dnapars, & dnapenny */
3983 gbases *p;
3984
3985 while (*garbage) {
3986 p = *garbage;
3987 *garbage = (*garbage)->next;
3988 free(p->base);
3989 free(p);
3990 }
3991 } /*freegarbage */
3992
3993
freegrbg(node ** grbg)3994 void freegrbg(node **grbg)
3995 {
3996 /* used in dnacomp, dnapars, & dnapenny */
3997 node *p;
3998
3999 while (*grbg) {
4000 p = *grbg;
4001 *grbg = (*grbg)->next;
4002 freenontip(p);
4003 free(p);
4004 }
4005 } /*freegrbg */
4006
4007
collapsetree(node * p,node * root,node ** grbg,pointarray treenode,long * zeros)4008 void collapsetree(node *p, node *root, node **grbg, pointarray treenode,
4009 long *zeros)
4010 {
4011 /* Recurse through tree searching for zero length brances between */
4012 /* nodes (not to tips). If one exists, collapse the nodes together, */
4013 /* removing the branch. */
4014 node *q, *x1, *y1, *x2, *y2;
4015 long i, j, index, index2, numd;
4016 if (p->tip)
4017 return;
4018 q = p->next;
4019 do {
4020 if (!q->back->tip && q->v == 0.000000) {
4021 /* merge the two nodes. */
4022 x1 = y2 = q->next;
4023 x2 = y1 = q->back->next;
4024 while(x1->next != q)
4025 x1 = x1-> next;
4026 while(y1->next != q->back)
4027 y1 = y1-> next;
4028 x1->next = x2;
4029 y1->next = y2;
4030
4031 index = q->index;
4032 index2 = q->back->index;
4033 numd = treenode[index-1]->numdesc + q->back->numdesc -1;
4034 chuck(grbg, q->back);
4035 chuck(grbg, q);
4036 q = x2;
4037
4038 /* update the indicies around the node circle */
4039 do{
4040 if(q->index != index){
4041 q->index = index;
4042 }
4043 q = q-> next;
4044 }while(x2 != q);
4045 updatenumdesc(treenode[index-1], root, numd);
4046
4047 /* Alter treenode to point to real nodes, and update indicies */
4048 /* acordingly. */
4049 j = 0; i=0;
4050 for(i = (index2-1); i < nonodes-1 && treenode[i+1]; i++){
4051 treenode[i]=treenode[i+1];
4052 treenode[i+1] = NULL;
4053 x1=x2=treenode[i];
4054 do{
4055 x1->index = i+1;
4056 x1 = x1 -> next;
4057 } while(x1 != x2);
4058 }
4059
4060 /* Create a new empty fork in the blank spot of treenode */
4061 x1=NULL;
4062 for(i=1; i <=3 ; i++){
4063 gnutreenode(grbg, &x2, index2, endsite, zeros);
4064 x2->next = x1;
4065 x1 = x2;
4066 }
4067 x2->next->next->next = x2;
4068 treenode[nonodes-1]=x2;
4069 if (q->back)
4070 collapsetree(q->back, root, grbg, treenode, zeros);
4071 } else {
4072 if (q->back)
4073 collapsetree(q->back, root, grbg, treenode, zeros);
4074 q = q->next;
4075 }
4076 } while (q != p);
4077 } /* collapsetree */
4078
4079
collapsebestrees(node ** root,node ** grbg,pointarray treenode,bestelm * bestrees,long * place,long * zeros,long chars,boolean recompute,boolean progress)4080 void collapsebestrees(node **root, node **grbg, pointarray treenode,
4081 bestelm *bestrees, long *place, long *zeros,
4082 long chars, boolean recompute, boolean progress)
4083
4084 {
4085 /* Goes through all best trees, collapsing trees where possible, and */
4086 /* deleting trees that are not unique. */
4087 long i,j, k, pos, nextnode, oldnextree;
4088 boolean found;
4089 node *dummy;
4090
4091 oldnextree = nextree;
4092 for(i = 0 ; i < (oldnextree - 1) ; i++){
4093 bestrees[i].collapse = true;
4094 }
4095
4096 if(progress)
4097 printf("Collapsing best trees\n ");
4098 k = 0;
4099 for(i = 0 ; i < (oldnextree - 1) ; i++){
4100 if(progress){
4101 if(i % (((oldnextree-1) / 72) + 1) == 0)
4102 putchar('.');
4103 fflush(stdout);
4104 }
4105 while(!bestrees[k].collapse)
4106 k++;
4107 /* Reconstruct tree. */
4108 *root = treenode[0];
4109 add(treenode[0], treenode[1], treenode[spp], root, recompute,
4110 treenode, grbg, zeros);
4111 nextnode = spp + 2;
4112 for (j = 3; j <= spp; j++) {
4113 if (bestrees[k].btree[j - 1] > 0)
4114 add(treenode[bestrees[k].btree[j - 1] - 1], treenode[j - 1],
4115 treenode[nextnode++ - 1], root, recompute, treenode, grbg,
4116 zeros);
4117 else
4118 add(treenode[treenode[-bestrees[k].btree[j - 1]-1]->back->index-1],
4119 treenode[j - 1], NULL, root, recompute, treenode, grbg, zeros);
4120 }
4121 reroot(treenode[outgrno - 1], *root);
4122
4123 treelength(*root, chars, treenode);
4124 collapsetree(*root, *root, grbg, treenode, zeros);
4125 savetree(*root, place, treenode, grbg, zeros);
4126 /* move everything down in the bestree list */
4127 for(j = k ; j < (nextree - 2) ; j++){
4128 memcpy(bestrees[j].btree, bestrees[j + 1].btree, spp * sizeof(long));
4129 bestrees[j].gloreange = bestrees[j + 1].gloreange;
4130 bestrees[j + 1].gloreange = false;
4131 bestrees[j].locreange = bestrees[j + 1].locreange;
4132 bestrees[j + 1].locreange = false;
4133 bestrees[j].collapse = bestrees[j + 1].collapse;
4134 }
4135 pos=0;
4136 findtree(&found, &pos, nextree-1, place, bestrees);
4137
4138 /* put the new tree at the end of the list if it wasn't found */
4139 nextree--;
4140 if(!found)
4141 addtree(pos, &nextree, false, place, bestrees);
4142
4143 /* Deconstruct the tree */
4144 for (j = 1; j < spp; j++){
4145 re_move(treenode[j], &dummy, root, recompute, treenode,
4146 grbg, zeros);
4147 }
4148 }
4149 if (progress) {
4150 putchar('\n');
4151 #ifdef WIN32
4152 phyFillScreenColor();
4153 #endif
4154 }
4155 }
4156