source: tags/initial/GDE/PHYLIP/dnamlk2.c

Last change on this file was 2, checked in by oldcode, 24 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.8 KB
Line 
1#include "phylip.h"
2
3#define maxcategories   5    /* maximum number of site types                 */
4#define smoothings      2    /* number of passes through smoothing algorithm */
5#define iterations      10   /* number of iterates for each branch           */
6#define nmlngth         10   /* number of characters max. in species name    */
7
8#define epsilon         0.0001   /* used in makenewv, getthree, update */
9
10#define point           '.'
11
12#define ibmpc0          false
13#define ansi0           true
14#define vt520           false
15
16
17typedef enum {
18  A, C, G, T
19} base;
20typedef double sitelike[(short)T - (short)A + 1];
21typedef sitelike **phenotype;
22typedef double *contribarr;
23typedef char **sequence;
24
25
26typedef short longer[6];
27
28typedef struct node {
29  struct node *next, *back;
30  phenotype x;
31  double tyme, v;
32  short number, xcoord, ycoord, ymin, ymax;
33  boolean tip, haslength, initialized;
34} node;
35
36
37typedef struct tree {
38  node **nodep;
39  double likelihood;
40  node *root;
41} tree;
42
43
44extern short      categs,spp,endsite,sites,nonodes,weightsum,datasets,
45                  njumble,jumb;
46extern double     rate[maxcategories];
47extern double     xi,xv,freqa,freqc,freqg,freqt,freqr,freqy,freqar,freqcy,
48                  freqgr,freqty,lambda,lambda1,fracchange;
49extern boolean    usertree,globle,jumble,lengths,trout,weights,ctgry,auto_,
50                  printdata,progress,treeprint;
51extern tree       curtree,bestree,bestree2;
52extern contribarr probcat;
53extern double      **contribution;
54extern short       *alias,*ally,*location,*aliasweight;
55extern FILE       *infile,*outfile,*treefile;
56extern            Char **naym;
57extern            short *enterorder;
58extern            longer seed;
59
60/* end pseudo-include file */
61
62#define down            2
63#define over            60
64
65
66/* Local variables for maketree, propagated globally for C version: */
67short    k, which, maxwhich,col;
68double  like, bestyet, tdelta, lnlike, slope, curv, maxlogl;
69boolean lastrearr, smoothed;
70double  *l0gl;
71double  x[3], lnl[3];
72double  expon1i[maxcategories], expon1v[maxcategories],
73        expon2i[maxcategories], expon2v[maxcategories];
74node    *there;
75double  **l0gf;
76Char ch;
77
78double randum(seed)
79short *seed;
80{
81  /* random number generator -- slow but machine independent */
82  short i, j, k, sum;
83  longer mult, newseed;
84  double x;
85
86  mult[0] = 13;
87  mult[1] = 24;
88  mult[2] = 22;
89  mult[3] = 6;
90  for (i = 0; i <= 5; i++)
91    newseed[i] = 0;
92  for (i = 0; i <= 5; i++) {
93    sum = newseed[i];
94    k = i;
95    if (i > 3)
96      k = 3;
97    for (j = 0; j <= k; j++)
98      sum += mult[j] * seed[i - j];
99    newseed[i] = sum;
100    for (j = i; j <= 4; j++) {
101      newseed[j + 1] += newseed[j] / 64;
102      newseed[j] &= 63;
103    }
104  }
105  memcpy(seed, newseed, sizeof(longer));
106  seed[5] &= 3;
107  x = 0.0;
108  for (i = 0; i <= 5; i++)
109    x = x / 64.0 + seed[i];
110  x /= 4.0;
111  return x;
112}  /* randum */
113
114
115void exmake(lz, n)
116double lz;
117short n;
118{
119  /* pretabulate tables of exponentials so need not do for each site */
120  short i;
121  double rat;
122
123  for (i = 0; i < categs; i++) {
124    rat = rate[i];
125    switch (n) {
126
127    case 1:
128      expon1i[i] = exp(rat * xi * lz);
129      expon1v[i] = exp(rat * xv * lz);
130      break;
131
132    case 2:
133      expon2i[i] = exp(rat * xi * lz);
134      expon2v[i] = exp(rat * xv * lz);
135      break;
136    }
137  }
138}  /* exmake */
139
140double evaluate(p)
141node *p;
142{
143  double sum, sum2, sumc, y, lz, y1, z1zz, z1yy, prod12, prod1, prod2, prod3,
144         sumterm, lterm;
145  static contribarr like,term,clai;
146  short i, j, lai;
147  boolean recorduser;
148  node *q, *r;
149  sitelike x1, x2;
150
151  if (!like){
152    like   = (double *)Malloc(categs * sizeof(double));
153    term   = (double *)Malloc(categs * sizeof(double));
154    clai   = (double *)Malloc(categs * sizeof(double));
155  }
156  recorduser = (!auto_ && usertree);
157  sum = 0.0;
158  if (p != curtree.root) {
159    r = p;
160    q = p->back;
161    y = fabs(r->tyme - q->tyme);
162  } else {
163    r = p->next->back;
164    q = p->next->next->back;
165    y = r->tyme + q->tyme - 2 * p->tyme;
166  }
167  lz = -y;
168  exmake(lz, 1L);
169  for (i = 0; i < endsite; i++) {
170    for (j = 0; j < categs; j++) {
171      y1 = 1.0 - expon1v[j];
172      z1zz = expon1i[j] * expon1v[j];
173      z1yy = (1.0 - expon1i[j]) * expon1v[j];
174      memcpy(x1, r->x[i][j], sizeof(sitelike));
175      prod1 = freqa * x1[0] + freqc * x1[(short)C - (short)A] +
176              freqg * x1[(short)G - (short)A] + freqt * x1[(short)T - (short)A];
177      memcpy(x2, q->x[i][j], sizeof(sitelike));
178      prod2 = freqa * x2[0] + freqc * x2[(short)C - (short)A] +
179              freqg * x2[(short)G - (short)A] + freqt * x2[(short)T - (short)A];
180      prod3 = (x1[0] * freqa + x1[(short)G - (short)A] * freqg) *
181              (x2[0] * freqar + x2[(short)G - (short)A] * freqgr) +
182          (x1[(short)C - (short)A] * freqc + x1[(short)T - (short)A] * freqt) *
183          (x2[(short)C - (short)A] * freqcy + x2[(short)T - (short)A] * freqty);
184      prod12 = freqa * x1[0] * x2[0] +
185               freqc * x1[(short)C - (short)A] * x2[(short)C - (short)A] +
186               freqg * x1[(short)G - (short)A] * x2[(short)G - (short)A] +
187               freqt * x1[(short)T - (short)A] * x2[(short)T - (short)A];
188      term[j] = z1zz * prod12 + z1yy * prod3 + y1 * prod1 * prod2;
189    }
190    sumterm = 0.0;
191    for (j = 0; j < categs; j++)
192      sumterm += probcat[j] * term[j];
193    lterm = log(sumterm);
194    for (j = 0; j < categs; j++)
195      clai[j] = term[j] / sumterm;
196    memcpy(contribution[i], clai, categs * sizeof(double ));
197    if (recorduser)
198      l0gf[which - 1][i] = lterm;
199    sum += aliasweight[i] * lterm;
200  }
201  if (auto_) {
202    for (j = 0; j < categs; j++)
203      like[j] = 1.0;
204    for (i = 0; i < sites; i++) {
205      if ((ally[i] > 0) && (location[ally[i]-1] > 0)) {
206        sumc = 0.0;
207        for (k = 1; k <= categs; k++)
208          sumc += probcat[k - 1] * like[k - 1];
209        sumc *= lambda;
210        lai = location[ally[i] - 1];
211        memcpy(clai, contribution[lai - 1], categs * sizeof(double));
212        for (j = 0; j < categs; j++)
213          like[j] = (lambda1 * like[j] + sumc) * clai[j];
214      } else {
215        for (j = 0; j < categs; j++)
216          like[j] = lambda1 * like[j] + lambda;
217      }
218    }
219    sum2 = 0.0;
220    for (i = 0; i < categs; i++)
221      sum2 += probcat[i] * like[i];
222    sum += log(sum2);
223  }
224  curtree.likelihood = sum;
225  if (!recorduser)
226    return sum;
227  l0gl[which - 1] = sum;
228  if (which == 1) {
229    maxwhich = 1;
230    maxlogl = sum;
231    return sum;
232  }
233  if (sum > maxlogl) {
234    maxwhich = which;
235    maxlogl = sum;
236  }
237  return sum;
238}  /* evaluate */
239
240void nuview(p)
241node *p;
242{
243  short i, j;
244  double lw, ww1, ww2, zz1, zz2, vv1, vv2, yy1, yy2, sum1, sum2, sumr1, sumr2,
245         sumy1, sumy2;
246  node *q, *r;
247  sitelike xx1, xx2, xx3;
248
249  q = p->next->back;
250  r = p->next->next->back;
251  if (q != NULL)
252    lw = -fabs(p->tyme - q->tyme);
253  else
254    lw = 0.0;
255  exmake(lw, 1L);
256  if (r != NULL)
257    lw = -fabs(p->tyme - r->tyme);
258  else
259    lw = 0.0;
260  exmake(lw, 2L);
261  for (i = 0; i < endsite; i++) {
262    for (j = 0; j < categs; j++) {
263      ww1 = expon1i[j];
264      zz1 = expon1v[j];
265      vv1 = 1.0 - ww1;
266      yy1 = 1.0 - zz1;
267      ww2 = expon2i[j];
268      zz2 = expon2v[j];
269      vv2 = 1.0 - ww2;
270      yy2 = 1.0 - zz2;
271      if (q != NULL)
272        memcpy(xx1, q->x[i][j], sizeof(sitelike));
273      if (r != NULL)
274        memcpy(xx2, r->x[i][j], sizeof(sitelike));
275      if (q == NULL) {
276        xx1[0] = 1.0;
277        xx1[(short)C - (short)A] = 1.0;
278        xx1[(short)G - (short)A] = 1.0;
279        xx1[(short)T - (short)A] = 1.0;
280      }
281      if (r == NULL) {
282        xx2[0] = 1.0;
283        xx2[(short)C - (short)A] = 1.0;
284        xx2[(short)G - (short)A] = 1.0;
285        xx2[(short)T - (short)A] = 1.0;
286      }
287      sum1 = yy1 * (freqa * xx1[0] + freqc * xx1[(short)C - (short)A] +
288            freqg * xx1[(short)G - (short)A] + freqt * xx1[(short)T - (short)A]);
289      sum2 = yy2 * (freqa * xx2[0] + freqc * xx2[(short)C - (short)A] +
290            freqg * xx2[(short)G - (short)A] + freqt * xx2[(short)T - (short)A]);
291      sumr1 = vv1 * (freqar * xx1[0] + freqgr * xx1[(short)G - (short)A]);
292      sumr2 = vv2 * (freqar * xx2[0] + freqgr * xx2[(short)G - (short)A]);
293      sumy1 = vv1 * (freqcy * xx1[(short)C - (short)A] +
294                     freqty * xx1[(short)T - (short)A]);
295      sumy2 = vv2 * (freqcy * xx2[(short)C - (short)A] +
296                     freqty * xx2[(short)T - (short)A]);
297      xx3[0] = (sum1 + zz1 * (ww1 * xx1[0] + sumr1)) *
298               (sum2 + zz2 * (ww2 * xx2[0] + sumr2));
299      xx3[(short)C - (short)A] =
300        (sum1 + zz1 * (ww1 * xx1[(short)C - (short)A] + sumy1)) *
301        (sum2 + zz2 * (ww2 * xx2[(short)C - (short)A] + sumy2));
302      xx3[(short)G - (short)A] =
303        (sum1 + zz1 * (ww1 * xx1[(short)G - (short)A] + sumr1)) *
304        (sum2 + zz2 * (ww2 * xx2[(short)G - (short)A] + sumr2));
305      xx3[(short)T - (short)A] =
306        (sum1 + zz1 * (ww1 * xx1[(short)T - (short)A] + sumy1)) *
307        (sum2 + zz2 * (ww2 * xx2[(short)T - (short)A] + sumy2));
308      memcpy(p->x[i][j], xx3, sizeof(sitelike));
309    }
310  }
311}  /* nuview */
312
313void getthree(p)
314node *p;
315{
316  /* compute likelihood, slope, curvature at a new triple of points */
317  double tt, td;
318  node *q, *r, *sdown;
319
320  sdown = curtree.nodep[p->number - 1];
321  q = sdown->next->back;
322  r = sdown->next->next->back;
323  sdown = sdown->back;
324  td = fabs(tdelta);
325  if (td < epsilon)
326    td = epsilon;
327  tt = p->tyme;
328  if (q->tyme - tt < td)
329    td = (q->tyme - tt) / 2.0;
330  if (r->tyme - tt < td)
331    td = (r->tyme - tt) / 2.0;
332  if (sdown != NULL) {
333    if (tt - sdown->tyme < td)
334      td = (tt - sdown->tyme) / 2.0;
335  }
336  if (-tt > epsilon && td < epsilon)
337    td = epsilon;
338  if (tt == -td)
339    td += epsilon;
340  p->tyme = tt + td;
341  x[0] = tt + td;
342  nuview(p);
343  lnl[0] = evaluate(p);
344  p->tyme = tt - td;
345  x[2] = tt - td;
346  nuview(p);
347  lnl[2] = evaluate(p);
348  p->tyme = tt;
349  x[1] = tt;
350  nuview(p);
351  lnl[1] = evaluate(p);
352}  /* getthree */
353
354
355void makenewv(p)
356node *p;
357{
358  /* improve a node time */
359  short it, imin, imax, i;
360  double tt, tfactor, tlow, thigh, oldlike, ymin, ymax, s32, s21, yold;
361  boolean done, already;
362  node *q, *r, *s, *sdown;
363
364  s = curtree.nodep[p->number - 1];
365  q = s->next->back;
366  r = s->next->next->back;
367  sdown = s->back;
368  if (s == curtree.root)
369    tlow = -10.0;
370  else
371    tlow = sdown->tyme;
372  thigh = q->tyme;
373  if (r->tyme < thigh)
374    thigh = r->tyme;
375  done = (thigh - tlow < epsilon || thigh - s->tyme < epsilon
376           || s->tyme - tlow < epsilon);
377  it = 1;
378  tdelta = (thigh - tlow) / 10.0;
379  getthree(s);
380  tfactor = 1.0;
381  while (it < iterations && !done) {
382    tt = s->tyme;
383    ymax = lnl[0];
384    imax = 1;
385    for (i = 2; i <= 3; i++) {
386      if (lnl[i - 1] > ymax) {
387        ymax = lnl[i - 1];
388        imax = i;
389      }
390    }
391    if (imax != 2) {
392      ymax = x[1];
393      x[1] = x[imax - 1];
394      x[imax - 1] = ymax;
395      ymax = lnl[1];
396      lnl[1] = lnl[imax - 1];
397      lnl[imax - 1] = ymax;
398    }
399    tt = x[1];
400    oldlike = lnl[1];
401    yold = tt;
402    s32 = (lnl[2] - lnl[1]) / (x[2] - x[1]);
403    s21 = (lnl[1] - lnl[0]) / (x[1] - x[0]);
404    if (fabs(x[2] - x[0]) > epsilon)
405      curv = (s32 - s21) / ((x[2] - x[0]) / 2);
406    else
407      curv = 0.0;
408    slope = (s32 + s21) / 2 - curv *      (x[2] - 2 * x[1] + x[0]) / 4;
409    if (curv >= 0.0) {
410      if (slope < 0)
411        tdelta = -fabs(tdelta);
412      else
413        tdelta = fabs(tdelta);
414    } else
415      tdelta = -(tfactor * slope / curv);
416    if (tt + tdelta <= tlow + epsilon)
417      tdelta = (tlow - tt) / 2;
418    if (tt + tdelta >= thigh - epsilon)
419      tdelta = (thigh - tt) / 2;
420    tt += tdelta;
421    done = (fabs(yold - tt) < epsilon || fabs(tdelta) < epsilon);
422    s->tyme = tt;
423    nuview(s);
424    lnlike = evaluate(s);
425    ymin = lnl[0];
426    imin = 1;
427    for (i = 2; i <= 3; i++) {
428      if (lnl[i - 1] < ymin) {
429        ymin = lnl[i - 1];
430        imin = i;
431      }
432    }
433    already = (tt == x[0]) || (tt == x[1]) || (tt == x[2]);
434    if (!already && ymin < lnlike) {
435      x[imin - 1] = tt;
436      lnl[imin - 1] = lnlike;
437    }
438    if (already || lnlike < oldlike) {
439      tt = x[1];
440      tfactor /= 2;
441      tdelta /= 2;
442      curtree.likelihood = oldlike;
443      lnlike = oldlike;
444    } else
445      tfactor = 1.0;
446    p->tyme = tt;
447    p->next->tyme = tt;
448    p->next->next->tyme = tt;
449    nuview(p);
450    nuview(p->next);
451    nuview(p->next->next);
452    it++;
453  }
454  smoothed = smoothed && done;
455}  /* makenewv */
456
457void update(p)
458node *p;
459{
460  /* improve time and recompute views at a node */
461  if (p == NULL)
462    return;
463  if (p->back != NULL) {
464    if (!p->back->tip)
465      nuview(p->back);
466  }
467  if (p->next->back != NULL) {
468    if (!p->next->back->tip)
469      nuview(p->next->back);
470  }
471  if (p->next->next->back != NULL) {
472    if (!p->next->next->back->tip)
473      nuview(p->next->next->back);
474  }
475  if (!lengths) {
476    makenewv(p);
477    return;
478  }
479  nuview(p);
480  nuview(p->next);
481  nuview(p->next->next);
482}  /* update */
483
484void smooth(p)
485node *p;
486{
487  if (p == NULL)
488    return;
489  if (p->tip)
490    return;
491  update(p);
492  smooth(p->next->back);
493  update(p);
494  smooth(p->next->next->back);
495  update(p);
496}  /* smooth */
497
498void add(below, newtip, newfork)
499node *below, *newtip, *newfork;
500{
501  /* inserts the nodes newfork and its descendant, newtip, in the tree. */
502  short i;
503  boolean done;
504  node *p;
505
506  below = curtree.nodep[below->number - 1];
507  if (below->back != NULL)
508    below->back->back = newfork;
509  newfork->back = below->back;
510  below->back = newfork->next->next;
511  newfork->next->next->back = below;
512  newfork->next->back = newtip;
513  newtip->back = newfork->next;
514  if (newtip->tyme < below->tyme)
515    p = newtip;
516  else p = below;
517  newfork->tyme = p->tyme;
518  if (curtree.root == below)
519    curtree.root = newfork;
520  if (newfork->back != NULL) {
521    if (p->tyme > newfork->back->tyme)
522      newfork->tyme = (p->tyme + newfork->back->tyme) / 2.0;
523    else newfork->tyme = p->tyme - epsilon;
524    newfork->next->tyme = newfork->tyme;
525    newfork->next->next->tyme = newfork->tyme;
526    do {
527      p = curtree.nodep[p->back->number - 1];
528      done = (p == curtree.root);
529      if (!done)
530        done = (curtree.nodep[p->back->number - 1]->tyme < p->tyme - epsilon);
531      if (!done) {
532        curtree.nodep[p->back->number - 1]->tyme = p->tyme - epsilon;
533        curtree.nodep[p->back->number - 1]->next->tyme = p->tyme - epsilon;
534        curtree.nodep[p->back->number - 1]->next->next->tyme = p->tyme - epsilon;
535      }
536    } while (!done);
537  } else {
538      newfork->tyme = newfork->tyme - 0.1;
539      newfork->next->tyme = newfork->tyme;
540      newfork->next->next->tyme = newfork->tyme;
541    }
542  nuview(newfork);
543  nuview(newfork->next);
544  nuview(newfork->next->next);
545  smoothed = false;
546  i = 1;
547  while (i < smoothings && !smoothed) {
548    smoothed = true;
549    smooth(newfork);
550    smooth(newfork->back);
551    i++;
552  }
553}  /* add */
554
555void re_move(item, fork)
556node **item, **fork;
557{
558  /* removes nodes item and its ancestor, fork, from the tree.
559    the new descendant of fork's ancestor is made to be
560    fork's second descendant (other than item).  Also
561    returns pointers to the deleted nodes, item and fork */
562  node *p, *q;
563
564  if ((*item)->back == NULL) {
565    *fork = NULL;
566    return;
567  }
568  *fork = curtree.nodep[(*item)->back->number - 1];
569  if (curtree.root == *fork) {
570    if (*item == (*fork)->next->back)
571      curtree.root = (*fork)->next->next->back;
572    else
573      curtree.root = (*fork)->next->back;
574  }
575  p = (*item)->back->next->back;
576  q = (*item)->back->next->next->back;
577  if (p != NULL)
578    p->back = q;
579  if (q != NULL)
580    q->back = p;
581  (*fork)->back = NULL;
582  p = (*fork)->next;
583  while (p != *fork) {
584    p->back = NULL;
585    p = p->next;
586  }
587  (*item)->back = NULL;
588}  /* remove */
589
590void copynode(c, d)
591node *c, *d;
592{
593  short i, j;
594
595  for (i=0;i<endsite;i++)
596    for (j = 0; j < categs; j++)
597      memcpy(d->x[i][j],c->x[i][j], sizeof(sitelike));
598  d->tyme = c->tyme;
599  d->v = c->v;
600  d->xcoord = c->xcoord;
601  d->ycoord = c->ycoord;
602  d->ymin = c->ymin;
603  d->ymax = c->ymax;
604  d->haslength = c->haslength;
605  d->initialized = c->initialized;
606}  /* copynode */
607
608void copy_(a, b)
609tree *a, *b;
610{
611  short i, j=0;
612  node *p, *q;
613
614  for (i = 0; i < spp; i++) {
615    copynode(a->nodep[i], b->nodep[i]);
616    if (a->nodep[i]->back) {
617      if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->number - 1])
618        b->nodep[i]->back = b->nodep[a->nodep[i]->back->number - 1];
619      else if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->number - 1]->next)
620        b->nodep[i]->back = b->nodep[a->nodep[i]->back->number - 1]->next;
621      else
622        b->nodep[i]->back = b->nodep[a->nodep[i]->back->number - 1]->next->next;
623    }
624    else b->nodep[i]->back = NULL;
625  }
626  for (i = spp; i < nonodes; i++) {
627    p = a->nodep[i];
628    q = b->nodep[i];
629    for (j = 1; j <= 3; j++) {
630      copynode(p, q);
631      if (p->back) {
632        if (p->back == a->nodep[p->back->number - 1])
633          q->back = b->nodep[p->back->number - 1];
634        else if (p->back == a->nodep[p->back->number - 1]->next)
635          q->back = b->nodep[p->back->number - 1]->next;
636        else
637          q->back = b->nodep[p->back->number - 1]->next->next;
638      }
639      else
640        q->back = NULL;
641      p = p->next;
642      q = q->next;
643    }
644  }
645  b->likelihood = a->likelihood;
646  b->root = a->root;
647}  /* copy */
648
649void tryadd(p, item,nufork)
650node *p,**item,**nufork;
651
652{
653  /* temporarily adds one fork and one tip to the tree.
654    if the location where they are added yields greater
655    likelihood than other locations tested up to that
656    time, then keeps that location as there */
657  add(p, *item, *nufork);
658  like = evaluate(p);
659  if (lastrearr) {
660    if (like >= bestree.likelihood)
661      copy_(&curtree, &bestree);
662  }
663  if (like > bestyet) {
664    bestyet = like;
665    there = p;
666  }
667  re_move(item,nufork);
668}  /* tryadd */
669
670void addpreorder(p, item_, nufork_)
671node *p, *item_, *nufork_;
672{
673  /* traverses a binary tree, calling PROCEDURE tryadd
674    at a node before calling tryadd at its descendants */
675  node *item, *nufork;
676
677  item = item_;
678  nufork = nufork_;
679  if (p == NULL)
680    return;
681  tryadd(p, &item,&nufork);
682  if (!p->tip) {
683    addpreorder(p->next->back, item, nufork);
684    addpreorder(p->next->next->back,item,nufork);
685  }
686}  /* addpreorder */
687
688
689void tryrearr(p, success)
690node *p;
691boolean *success;
692{
693  /* evaluates one rearrangement of the tree.
694    if the new tree has greater likelihood than the old
695    one sets success = TRUE and keeps the new tree.
696    otherwise, restores the old tree */
697  node *frombelow, *whereto, *forknode;
698  double oldlike;
699
700  if (p == curtree.root)
701    return;
702  forknode = curtree.nodep[p->back->number - 1];
703  if (forknode == curtree.root)
704    return;
705  oldlike = bestyet;
706  if (forknode->next->back == p)
707    frombelow = forknode->next->next->back;
708  else
709    frombelow = forknode->next->back;
710  whereto = curtree.nodep[forknode->back->number - 1];
711  re_move(&p, &forknode);
712  nuview(whereto);
713  add(whereto, p, forknode);
714  like = evaluate(p);
715  if (like <= oldlike) {
716    re_move(&p, &forknode);
717    add(frombelow, p, forknode);
718  } else {
719    (*success) = true;
720    bestyet = like;
721  }
722}  /* tryrearr */
723
724void repreorder(p, success)
725node *p;
726boolean *success;
727{
728  /* traverses a binary tree, calling PROCEDURE tryrearr
729    at a node before calling tryrearr at its descendants */
730  if (p == NULL)
731    return;
732  tryrearr(p, success);
733  if (p->tip)
734    return;
735  if (!(*success))
736    repreorder(p->next->back, success);
737  if (!(*success))
738    repreorder(p->next->next->back,success);
739}  /* repreorder */
740
741Local Void rearrange(r)
742node **r;
743{
744  /* traverses the tree (preorder), finding any local
745    rearrangement which decreases the number of steps.
746    if traversal succeeds in increasing the tree's
747    likelihood, PROCEDURE rearrange runs traversal again */
748  boolean success;
749  success = true;
750  while (success) {
751    success = false;
752    repreorder(*r, &success);
753  }
754}  /* rearrange */
755
756void getch(c)
757Char *c;
758{
759  /* get next nonblank character */
760  do {
761    if (eoln(infile)) {
762      fscanf(infile, "%*[^\n]");
763      getc(infile);
764    }
765    *c = getc(infile);
766    if (*c == '\n')
767      *c = ' ';
768  } while (*c == ' ');
769}  /* getch */
770
771void findch(c,rparens)
772Char c;
773short *rparens;
774{
775  /* scan forward until find character c */
776  boolean done;
777
778  done = false;
779  while (!(done)) {
780    if (c == ',') {
781      if (ch == '(' || ch == ')' || ch == ';') {
782        printf("\nERROR IN USER TREE: UNMATCHED PARENTHESIS OR MISSING COMMA\n");
783        exit(-1);
784      } else if (ch == ',')
785        done = true;
786    } else if (c == ')') {
787      if (ch == '(' || ch == ',' || ch == ';') {
788        printf("\nERROR IN USER TREE: UNMATCHED PARENTHESIS OR NOT BIFURCATED NODE\n");
789        exit(-1);
790      } else {
791        if (ch == ')') {
792          done = true;
793          (*rparens)++;
794        }
795      }
796    } else if (c == ';') {
797      if (ch != ';') {
798        printf("\nERROR IN USER TREE: UNMATCHED PARENTHESIS OR MISSING SEMICOLON\n");
799        exit(-1);
800      } else
801        done = true;
802    }
803    if ((done && ch == ')') || !(done))
804      getch(&ch);
805  }
806}  /* findch */
807
808
809void processlength(p)
810node *p;
811{
812  short digit, ordzero;
813  double valyew, divisor;
814  boolean pointread;
815
816  ordzero = '0';
817  pointread = false;
818  valyew = 0.0;
819  divisor = 1.0;
820  getch(&ch);
821  digit = ch - ordzero;
822  while (((unsigned short)digit <= 9) ||         ch == point){
823    if (ch == point)
824      pointread = true;
825    else {
826      valyew = valyew * 10.0 + digit;
827      if (pointread)
828        divisor *= 10.0;
829    }
830    getch(&ch);
831    digit = ch - ordzero;
832  }
833  p->v = valyew / divisor / fracchange;
834}  /* processlength */
835
836void addelement(p, nextnode,lparens,rparens,names)
837node **p;
838short *nextnode,*lparens,*rparens;
839boolean *names;
840{
841  /* recursive procedure adds nodes to user-defined tree */
842  node *q;
843  short i, n;
844  boolean found;
845  Char str[nmlngth];
846
847  getch(&ch);
848  if (ch == '(') {
849    if ((*lparens) >= spp - 1) {
850      printf("\nERROR IN USER TREE: TOO MANY LEFT PARENTHESES\n");
851      exit(-1);
852    } else {
853      (*lparens)++;
854      (*nextnode)++;
855      q = curtree.nodep[(*nextnode) - 1];
856      addelement(&q->next->back, nextnode,lparens,rparens,names);
857      q->next->back->back = q->next;
858      findch(',',rparens);
859      addelement(&q->next->next->back, nextnode,lparens,rparens,names);
860      q->next->next->back->back = q->next->next;
861      findch(')',rparens);
862      *p = q;
863    }
864  } else  {
865    for (i = 0; i < nmlngth; i++)
866      str[i] = ' ';
867    n = 1;
868    do {
869      if (ch == '_')
870        ch = ' ';
871      str[n - 1] = ch;
872      if (eoln(infile)) {
873        fscanf(infile, "%*[^\n]");
874        getc(infile);
875      }
876      ch = getc(infile);
877      if (ch == '\n')
878        ch = ' ';
879      n++;
880    } while (ch != ':' && ch != ',' && ch != ')' && n <= nmlngth);
881    if (ch == ')')
882      (*rparens)++;
883    n = 1;
884    do {
885      found = true;
886      for (i = 0; i < nmlngth; i++)
887        found = (found && str[i] == naym[n - 1][i]);
888      if (found) {
889        if (names[n - 1] == false) {
890          *p = curtree.nodep[n - 1];
891          names[n - 1] = true;
892        } else {
893          printf("\nERROR IN USER TREE: DUPLICATE NAME FOUND -- ");
894          for (i = 0; i < nmlngth; i++)
895            putchar(naym[n - 1][i]);
896          putchar('\n');
897          exit(-1);
898        }
899      } else
900        n++;
901    } while (!(n > spp || found));
902    if (n > spp) {
903      printf("Cannot find species: ");
904      for (i = 0; i < nmlngth; i++)
905        putchar(str[i]);
906      putchar('\n');
907      exit(-1);
908    }
909  }
910  if (!lengths)
911    return;
912  if (ch == ':' || ch == ';') {
913    if ((*rparens) < spp && ch == ':')
914      processlength(*p);
915    return;
916  }
917  printf("ERROR IN USER TREE: UNMATCHED PARENTHESES\n");
918  printf(" OR BRANCH LENGTHS OPTION CHOSEN IN MENU\n");
919  printf(" WITH NO INPUT BRANCH LENGTHS\n");
920  exit(-1);
921}  /* addelement */
922
923void tymetrav(p, x)
924node *p;
925double *x;
926{
927  /* set up times of nodes */
928  if (!p->tip) {
929    tymetrav(p->next->back, x);
930    tymetrav(p->next->next->back,x);
931  } else
932    (*x)     = 0.0;
933    p->tyme  = (*x);
934    (*x)    -= p->v;
935}  /* tymetrav */
936
937void treeread()
938{
939  /* read in user-defined tree and set it up */
940  short nextnode, lparens, rparens;
941  boolean *names;
942  double x;
943  short i;
944
945  nextnode = spp;
946  names = (boolean *)Malloc(spp*sizeof(boolean));
947  for (i = 0; i < spp; i++)
948    names[i] = false;
949  lparens = 0;
950  rparens = 0;
951  addelement(&curtree.root, &nextnode,&lparens,&rparens,names);
952  curtree.root->back = NULL;
953  findch(';',&rparens);
954  if (progress)
955    printf("\n\n");
956  if (lengths)
957    tymetrav(curtree.root, &x);
958  fscanf(infile, "%*[^\n]");
959  getc(infile);
960  free(names);
961}  /* treeread */
962
963void coordinates(p, tipy)
964node *p;
965short *tipy;
966{
967  /* establishes coordinates of nodes */
968  node *q, *first, *last;
969
970  if (p->tip) {
971    p->xcoord = 0;
972    p->ycoord = (*tipy);
973    p->ymin   = (*tipy);
974    p->ymax   = (*tipy);
975    (*tipy)  += down;
976    return;
977  }
978  q = p->next;
979  do {
980    coordinates(q->back, tipy);
981    q = q->next;
982  } while (p != q);
983  first = p->next->back;
984  q = p->next;
985  while (q->next != p)
986    q = q->next;
987  last = q->back;
988  p->xcoord = (short)(0.5 - over * p->tyme);
989  p->ycoord = (first->ycoord + last->ycoord) / 2;
990  p->ymin = first->ymin;
991  p->ymax = last->ymax;
992}  /* coordinates */
993
994void drawline(i, scale)
995short i;
996double scale;
997{
998  /* draws one row of the tree diagram by moving up tree */
999  node *p, *q, *r, *first, *last;
1000  short n, j;
1001  boolean extra, done;
1002
1003  p = curtree.root;
1004  q = curtree.root;
1005  extra = false;
1006  if (p->ycoord == i) {
1007    if (p->number - spp >= 10)
1008      fprintf(outfile, "-%2hd", p->number - spp);
1009    else
1010      fprintf(outfile, "--%hd", p->number - spp);
1011    extra = true;
1012  } else
1013    fprintf(outfile, "  ");
1014  do {
1015    if (!p->tip) {
1016      r = p->next;
1017      done = false;
1018      do {
1019        if (i >= r->back->ymin && i <= r->back->ymax) {
1020          q = r->back;
1021          done = true;
1022        }
1023        r = r->next;
1024      } while (!(done || r == p));
1025      first = p->next->back;
1026      r = p->next;
1027      while (r->next != p)
1028        r = r->next;
1029      last = r->back;
1030    }
1031    done = (p == q);
1032    n = (short)(scale * (p->xcoord - q->xcoord) + 0.5);
1033    if (n < 3 && !q->tip)
1034      n = 3;
1035    if (extra) {
1036      n--;
1037      extra = false;
1038    }
1039    if (q->ycoord == i && !done) {
1040      if (p->ycoord != q->ycoord)
1041        putc('+', outfile);
1042      else
1043        putc('-', outfile);
1044      if (!q->tip) {
1045        for (j = 1; j <= n - 2; j++)
1046          putc('-', outfile);
1047        if (q->number - spp >= 10)
1048          fprintf(outfile, "%2hd", q->number - spp);
1049        else
1050          fprintf(outfile, "-%hd", q->number - spp);
1051        extra = true;
1052      } else {
1053        for (j = 1; j < n; j++)
1054          putc('-', outfile);
1055      }
1056    } else if (!p->tip) {
1057      if (last->ycoord > i && first->ycoord < i && i != p->ycoord) {
1058        putc('!', outfile);
1059        for (j = 1; j < n; j++)
1060          putc(' ', outfile);
1061      } else {
1062        for (j = 1; j <= n; j++)
1063          putc(' ', outfile);
1064      }
1065    } else {
1066      for (j = 1; j <= n; j++)
1067        putc(' ', outfile);
1068    }
1069    if (p != q)
1070      p = q;
1071  } while (!done);
1072  if (p->ycoord == i && p->tip) {
1073    for (j = 0; j < nmlngth; j++)
1074      putc(naym[p->number - 1][j], outfile);
1075  }
1076  putc('\n', outfile);
1077}  /* drawline */
1078
1079void printree()
1080{
1081  /* prints out diagram of the tree */
1082 /* Local variables for printree: */
1083  short tipy;
1084  double scale;
1085  short i;
1086  node *p;
1087
1088  if (!treeprint)
1089    return;
1090  putc('\n', outfile);
1091  tipy = 1;
1092  coordinates(curtree.root, &tipy);
1093  p = curtree.root;
1094  while (!p->tip)
1095    p = p->next->back;
1096  scale = 1.0 / (short)(p->tyme - curtree.root->tyme + 1.000);
1097  putc('\n', outfile);
1098  for (i = 1; i <= tipy - down; i++)
1099    drawline(i, scale);
1100  putc('\n', outfile);
1101}  /* printree */
1102
1103#undef down
1104#undef over
1105
1106void describe(p)
1107node *p;
1108{
1109  short i;
1110  double v;
1111
1112  if (p == curtree.root)
1113    fprintf(outfile, " root         ");
1114  else
1115    fprintf(outfile, "%4hd          ", p->back->number - spp);
1116  if (p->tip) {
1117    for (i = 0; i < nmlngth; i++)
1118      putc(naym[p->number - 1][i], outfile);
1119  } else
1120    fprintf(outfile, "%4hd      ", p->number - spp);
1121  if (p != curtree.root) {
1122    fprintf(outfile, "%10.5f", fracchange * (p->tyme - curtree.root->tyme));
1123    v = fracchange * (p->tyme - curtree.nodep[p->back->number - 1]->tyme);
1124    fprintf(outfile, "%12.5f", v);
1125  }
1126  putc('\n', outfile);
1127  if (!p->tip) {
1128    describe(p->next->back);
1129    describe(p->next->next->back);
1130  }
1131}  /* describe */
1132
1133void summarize()
1134{
1135  short i, j, mx;
1136  double mode, sum;
1137  double *like,*nulike;
1138  short **mp;
1139
1140  like    = (double *)Malloc(categs * sizeof(double));
1141  nulike  = (double *)Malloc(categs * sizeof(double));
1142  mp      = (short **)Malloc(sites * sizeof(short *));
1143  for (i=0; i<=sites-1; i++)
1144     mp[i] = (short *)Malloc(sizeof(short)*categs);
1145  fprintf(outfile, "\nLn Likelihood = %11.5f\n\n", curtree.likelihood);
1146  fprintf(outfile, " Ancestor      Node      Node Time     Length\n");
1147  fprintf(outfile, " --------      ----      ---- ----     ------\n");
1148  describe(curtree.root);
1149  putc('\n', outfile);
1150  if (ctgry && categs > 1) {
1151    for (i = 0; i < categs; i++)
1152      like[i] = 1.0;
1153    for (i = sites - 1; i >= 0; i--) {
1154      sum = 0.0;
1155      for (j = 0; j < categs; j++) {
1156        nulike[j] = (lambda1 + lambda * probcat[j]) * like[j]; 
1157        mp[i][j] = j + 1;
1158        for (k = 1; k <= categs; k++) {
1159          if (k != j + 1) {
1160            if (lambda * probcat[k - 1] * like[k - 1] > nulike[j]) {
1161              nulike[j] = lambda * probcat[k - 1] * like[k - 1];
1162              mp[i][j] = k;
1163            }
1164          }
1165        }
1166        if ((ally[i] > 0) && (location[ally[i]-1] > 0))
1167          nulike[j] *= contribution[location[ally[i] - 1] - 1][j];
1168        sum += nulike[j];
1169      }
1170      for (j = 0; j < categs; j++)
1171        nulike[j] /= sum;
1172      memcpy(like, nulike, categs * sizeof(double));
1173    }
1174    mode = 0.0;
1175    mx = 1;
1176    for (i = 1; i <= categs; i++) {
1177      if (probcat[i - 1] * like[i - 1] > mode) {
1178        mx = i;
1179        mode = probcat[i - 1] * like[i - 1];
1180      }
1181    }
1182    fprintf(outfile,
1183 "\nCombination of categories that contributes the most to the likelihood:\n");
1184    for (i = 1; i <= nmlngth + 3; i++)
1185      putc(' ', outfile);
1186    for (i = 1; i <= sites; i++) {
1187      fprintf(outfile, "%hd", mx);
1188      if (i % 10 == 0)
1189        putc(' ', outfile);
1190      if (i % 60 == 0 && i != sites) {
1191        putc('\n', outfile);
1192        for (j = 1; j <= nmlngth + 3; j++)
1193          putc(' ', outfile);
1194      }
1195    mx = mp[i - 1][mx - 1];
1196    }
1197  }
1198  putc('\n', outfile);
1199  free(like);
1200  free(nulike);
1201  for (i=0;i<sites;++i)
1202    free(mp[i]);
1203  free(mp);
1204}  /* summarize */
1205
1206void treeout(p)
1207node *p;
1208{
1209  /* write out file with representation of final tree */
1210  short i, n, w;
1211  Char c;
1212  double x;
1213
1214  if (p->tip) {
1215    n = 0;
1216    for (i = 1; i <= nmlngth; i++) {
1217      if (naym[p->number - 1][i - 1] != ' ')
1218        n = i;
1219    }
1220    for (i = 0; i < n; i++) {
1221      c = naym[p->number - 1][i];
1222      if (c == ' ')
1223        c = '_';
1224      putc(c, treefile);
1225    }
1226    col += n;
1227  } else {
1228    putc('(', treefile);
1229    col++;
1230    treeout(p->next->back);
1231    putc(',', treefile);
1232    col++;
1233    if (col > 55) {
1234      putc('\n', treefile);
1235      col = 0;
1236    }
1237    treeout(p->next->next->back);
1238    putc(')', treefile);
1239    col++;
1240  }
1241  if (p == curtree.root) {
1242    fprintf(treefile, ";\n");
1243    return;
1244  }
1245  x = fracchange * (p->tyme - curtree.nodep[p->back->number - 1]->tyme);
1246  if (x > 0.0)
1247    w = (short)(0.4342944822 * log(x));
1248  else if (x == 0.0)
1249    w = 0;
1250  else
1251    w = (short)(0.4342944822 * log(-x)) + 1;
1252  if (w < 0)
1253    w = 0;
1254  fprintf(treefile, ":%*.5f", (int)(w + 7), x);
1255  col += w + 8;
1256}  /* treeout */
1257
1258
1259void nodeinit(p)
1260node *p;
1261{
1262  /* set up times at one node */
1263  node *q, *r;
1264  double lowertyme;
1265
1266  q = p->next->back;
1267  r = p->next->next->back;
1268  lowertyme = q->tyme;
1269  if (r->tyme < q->tyme)
1270    lowertyme = r->tyme;
1271  p->tyme = lowertyme - 0.1;
1272  p->next->tyme = p->tyme;
1273  p->next->next->tyme = p->tyme;
1274  q->v = q->tyme - p->tyme;
1275  p->next->v = q->v;
1276  r->v = r->tyme - p->tyme;
1277  p->next->next->v = r->v;
1278}  /* nodeinit */
1279
1280void initrav(p)
1281node *p;
1282{
1283  /* traverse to set up times throughout tree */
1284  if (p->tip)
1285    return;
1286  initrav(p->next->back);
1287  initrav(p->next->next->back);
1288  nodeinit(p);
1289}  /* initrav */
1290
1291
1292void travinit(p)
1293node *p;
1294{
1295  /* traverse to set up initial values */
1296  if (p == NULL)
1297    return;
1298  if (p->tip)
1299    return;
1300  if (p->initialized)
1301    return;
1302  travinit(p->next->back);
1303  travinit(p->next->next->back);
1304  nuview(p);
1305  p->initialized = true;
1306}  /* travinit */
1307
1308void travsp(p)
1309node *p;
1310{
1311  /* traverse to find tips */
1312  if (p == curtree.root)
1313    travinit(p);
1314  if (p->tip)
1315    travinit(p->back);
1316  else {
1317    travsp(p->next->back);
1318    travsp(p->next->next->back);
1319  }
1320}  /* travsp */
1321
1322void treevaluate()
1323{
1324  /* evaluate likelihood of tree, after iterating branch lengths */
1325  short i;
1326  double dummy;
1327
1328  for (i = 0; i < spp; i++)
1329    curtree.nodep[i]->initialized = false;
1330  for (i = spp; i < nonodes; i++) {
1331    curtree.nodep[i]->initialized = false;
1332    curtree.nodep[i]->next->initialized = false;
1333    curtree.nodep[i]->next->next->initialized = false;
1334  }
1335  if (!lengths)
1336    initrav(curtree.root);
1337  travsp(curtree.root);
1338  for (i = 1; i <= smoothings * 4; i++)
1339    smooth(curtree.root);
1340  dummy = evaluate(curtree.root);
1341}  /* treevaluate */
1342
1343
1344void maketree()
1345{
1346  /* constructs a binary tree from the pointers in curtree.nodep,
1347    adds each node at location which yields highest likelihood
1348    then rearranges the tree for greatest likelihood */
1349
1350  short i, j, numtrees, num;
1351  double bestlike, gotlike, sum, sum2, sd;
1352  node *item, *nufork, *dummy, *p1, *p2, *p3;
1353  double TEMP;
1354
1355  if (!usertree) {
1356    for (i = 1; i <= spp; i++)
1357      enterorder[i - 1] = i;
1358    if (jumble) {
1359      for (i = 0; i < spp; i++) {
1360        j = (short)(randum(seed) * spp) + 1;
1361        k = enterorder[j - 1];
1362        enterorder[j - 1] = enterorder[i];
1363        enterorder[i] = k;
1364      }
1365    }
1366    curtree.root = curtree.nodep[spp];
1367    if (progress) {
1368      printf("\nAdding species:\n");
1369      printf("   ");
1370    }
1371    add(curtree.nodep[enterorder[0] - 1], curtree.nodep[enterorder[1] - 1],
1372        curtree.nodep[spp]);
1373    if (progress) {
1374      for (i = 0; i < nmlngth; i++)
1375        putchar(naym[enterorder[0] - 1][i]);
1376      printf("\n   ");
1377      for (i = 0; i < nmlngth; i++)
1378        putchar(naym[enterorder[1] - 1][i]);
1379      putchar('\n');
1380    }
1381    lastrearr = false;
1382    for (i = 3; i <= spp; i++) {
1383      bestree.likelihood = -999999.0;
1384      bestyet = -999999.0;
1385      there = curtree.root;
1386      item = curtree.nodep[enterorder[i - 1] - 1];
1387      nufork = curtree.nodep[spp + i - 2];
1388      addpreorder(curtree.root, item, nufork);
1389      add(there, item, nufork);
1390      like = bestyet;
1391      rearrange(&curtree.root);
1392      if (curtree.likelihood > bestree.likelihood)
1393        copy_(&curtree, &bestree);
1394      if (progress) {
1395        printf("   ");
1396        for (j = 0; j < nmlngth; j++)
1397          putchar(naym[enterorder[i - 1] - 1][j]);
1398        putchar('\n');
1399      }
1400      lastrearr = (i == spp);
1401      if (lastrearr && globle) {
1402        if (progress) {
1403          printf("Doing global rearrangements\n");
1404          printf("  !");
1405          for (j = 1; j <= nonodes; j++)
1406            putchar('-');
1407          printf("!\n");
1408        }
1409        bestlike = bestyet;
1410        do {
1411          if (progress)
1412            printf("   ");
1413          gotlike = bestlike;
1414          for (j = 0; j < nonodes; j++) {
1415            bestyet = -999999.00;
1416            item = curtree.nodep[j];
1417            if (item != curtree.root) {
1418              nufork = curtree.nodep[curtree.nodep[j]->back->number - 1];
1419              re_move(&item, &nufork);
1420              there = curtree.root;
1421              addpreorder(curtree.root, item, nufork);
1422              add(there, item, nufork);
1423            }
1424            if (progress)
1425              putchar('.');
1426          }
1427          if (progress)
1428            putchar('\n');
1429        } while (bestlike < gotlike);
1430      }
1431    }
1432    if (njumble > 1 && lastrearr) {
1433      for (i = 0; i < spp; i++)
1434        re_move(&curtree.nodep[i], &dummy);
1435      if (jumb == 1 || bestree2.likelihood < bestree.likelihood)
1436        copy_(&bestree, &bestree2);
1437    }
1438    if (jumb == njumble) {
1439      if (njumble > 1) copy_(&bestree2, &curtree);
1440      else copy_(&bestree, &curtree);
1441      fprintf(outfile, "\n\n");
1442      curtree.likelihood = evaluate(curtree.root);
1443      printree();
1444      summarize();
1445      if (trout) {
1446        col = 0;
1447        treeout(curtree.root);
1448      }
1449    }
1450  } else {
1451    fscanf(infile, "%hd%*[^\n]", &numtrees);
1452    l0gl = (double *)Malloc(numtrees * sizeof(double));
1453    l0gf = (double **)Malloc(numtrees * sizeof(double *));
1454    for (i=0;i<numtrees;++i)
1455      l0gf[i] = (double *)Malloc(endsite * sizeof(double));
1456    getc(infile);
1457    if (treeprint) {
1458      fprintf(outfile, "User-defined tree");
1459      if (numtrees > 1)
1460        putc('s', outfile);
1461      fprintf(outfile, ":\n\n");
1462    }
1463    fprintf(outfile, "\n\n");
1464    which = 1;
1465    while (which <= numtrees) {
1466      setuptree(&curtree);
1467      treeread();
1468      treevaluate();
1469      printree();
1470      summarize();
1471      if (trout) {
1472        col = 0;
1473        treeout(curtree.root);
1474      }
1475      which++;
1476    }
1477    putc('\n', outfile);
1478    if (!auto_ && numtrees > 1 && weightsum > 1 ) {
1479      fprintf(outfile, "Tree    Ln L      Diff Ln L     Its S.D.");
1480      fprintf(outfile, "   Significantly worse?\n\n");
1481      num = numtrees;
1482      for (which = 1; which <= num; which++) {
1483        fprintf(outfile, "%3hd%12.5f", which, l0gl[which - 1]);
1484        if (maxwhich == which)
1485          fprintf(outfile, "  <------ best\n");
1486        else {
1487          sum = 0.0;
1488          sum2 = 0.0;
1489          for (j = 0; j < endsite; j++) {
1490            sum += aliasweight[j] *
1491                   (l0gf[maxwhich - 1][j] - l0gf[which - 1][j]);
1492            TEMP = l0gf[maxwhich - 1][j] - l0gf[which - 1][j];
1493            sum2 += aliasweight[j] * (TEMP * TEMP);
1494          }
1495          sd = sqrt(weightsum / (weightsum - 1.0) * (sum2 - sum * sum / weightsum));
1496          fprintf(outfile, "%12.5f%12.4f",
1497                  l0gl[which - 1] - maxlogl, sd);
1498          if (sum > 1.95996 * sd)
1499            fprintf(outfile, "           Yes\n");
1500          else
1501            fprintf(outfile, "            No\n");
1502        }
1503      }
1504      fprintf(outfile, "\n\n");
1505    }
1506    free(l0gl);
1507    for (i=0;i<numtrees;++i)
1508      free(l0gf[i]);
1509    free(l0gf);
1510  }
1511  if (jumb == njumble) {
1512    if (progress) {
1513      printf("\nOutput written to output file\n\n");
1514      if (trout)
1515        printf("Tree(s) also written onto file\n\n");
1516    }
1517    for (i=0; i<spp;i++){
1518      for (j=0;j<endsite;++j){
1519        free(curtree.nodep[i]->x[j]);
1520        if (!usertree) {
1521          free(bestree.nodep[i]->x[j]);
1522          if (njumble > 1)
1523            free(bestree2.nodep[i]->x[j]);
1524        }
1525      }
1526      free(curtree.nodep[i]->x);
1527      if (!usertree) {
1528        free(bestree.nodep[i]->x);
1529        if (njumble > 1)
1530          free(bestree2.nodep[i]->x);
1531      }
1532    }
1533    for (i=spp; i < nonodes ; i++){
1534      p1 = curtree.nodep[i];
1535      if (!usertree) {
1536        p2 = bestree.nodep[i];
1537        if (njumble > 1)
1538          p3 = bestree2.nodep[i];
1539      }
1540      for (k=0;k<3;++k) {
1541        for (j=0;j<endsite;++j){
1542          free(p1->x[j]);
1543          if (!usertree) {
1544            free(p2->x[j]);
1545            if (njumble > 1)
1546              free(p3->x[j]);
1547          }
1548        }
1549        free(p1->x);
1550        if (!usertree) {
1551          free(p2->x);
1552          if (njumble > 1)
1553            free(p3->x);
1554        }
1555        p1=p1->next;
1556        if (!usertree) {
1557          p2=p2->next;
1558          if (njumble > 1)
1559          p3=p3->next;
1560        }
1561      }
1562    }
1563  }
1564}  /* maketree */
1565
1566
Note: See TracBrowser for help on using the repository browser.