source: tags/initial/GDE/PHYLIP/protml.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: 118.3 KB
Line 
1/* Output from p2c, the Pascal-to-C translator */
2/* From input file "protml.pas" */
3
4
5#include "p2c.h"
6/* p2c: protml.pas, line 2:
7 * Note: Unexpected name "seqfile" in program header [262] */
8/* p2c: protml.pas, line 2:
9 * Note: Unexpected name "lklfile" in program header [262] */
10/* p2c: protml.pas, line 2:
11 * Note: Unexpected name "tpmfile" in program header [262] */
12
13
14/* Maximum Likelihood Inference of Protein Phylogeny */
15/* (seqfile, output); */
16
17/*            PROTML  Ver 1.01    January 6, 1993          */
18/*               J.Adachi  and  M.Hasegawa                 */
19/*        The Institute of Statistical Mathematics         */
20/*     4-6-7 Minami-Azabu, Minato-ku, Tokyo 106, Japan     */
21
22#define maxsp           20   /* maximum number of species                 */
23#define maxnode         37   /* maxsp * 2 - 3                             */
24#define maxpair         190   /* maxsp * (maxsp-1) / 2                     */
25#define maxsite         500   /* maximum number of sites                   */
26#define maxptrn         500   /* maximum number of different site patterns */
27#define maxtree         15   /* maximum number of user trees              */
28#define maxsmooth       30   /* number of smoothing algorithm             */
29#define maxiterat       10   /* number of iterates of Newton method       */
30
31#define epsilon         1.0e-5
32    /* stopping value of error                   */
33#define minarc          1.0e-3
34    /* lower limit on number of subsutitutions   */
35#define maxarc          3.0e+2
36    /* upper limit on number of subsutitutions   */
37
38#define prprtn          1   /* proportion of branch length               */
39#define maxboot         10000
40    /* number of bootstrap resamplings           */
41#define maxexe          1   /* number of jobs                            */
42#define maxline         60   /* length of sequence output per line        */
43#define maxname         10   /* max. number of characters in species name */
44#define maxami          20   /* number of amino acids                     */
45
46#define minreal         1.0e-55
47/* if job is in underflow error,
48                        then increase this value                  */
49
50#define seqfname        "seqfile"   /* input file of sequences data    */
51#define tpmfname        "default_not_use"
52    /* input file of trans probability */
53#define lklfname        "default_not_use"
54    /* output file of ln likelihood    */
55
56
57/* maxnode   = maxsp*2 -3;            */
58/* maxpair   = maxsp*(maxsp-1) DIV 2; */
59
60typedef enum {
61  ami
62} nacidty;
63typedef enum {
64  AA, RR, NN, DD, CC, QQ, EE, GG, HH, II, LL, KK, MM, FF, PP_, SS, TT, WW, YY,
65  VV
66} amity;
67typedef double aryamity[(long)VV - (long)AA + 1];
68typedef long iaryamity[(long)VV - (long)AA + 1];
69typedef aryamity daryamity[(long)VV - (long)AA + 1];
70typedef daryamity taryamity[(long)VV - (long)AA + 1];
71typedef long niaryamity[maxami];
72typedef double naryamity[maxami];
73typedef naryamity ndaryamity[maxami];
74typedef aryamity probamity[maxptrn];
75typedef Char charseqty[maxsp + 1][maxsite];
76typedef Char namety[maxname];
77typedef struct node *arynodepty[maxnode + 1];
78typedef double lengthty[maxnode];
79typedef long ilengthty[maxnode];
80typedef lengthty dlengthty[maxnode];
81typedef long sitety[maxsite];
82typedef long itreety[maxtree + 1];
83typedef double rtreety[maxtree + 1];
84typedef double lptrnty[maxptrn];
85typedef lptrnty ltrptty[maxtree + 1];
86typedef long spity[maxsp];
87typedef double rpairty[maxpair];
88typedef double rnodety[maxnode];
89typedef rnodety dpanoty[maxpair];
90typedef rpairty dnopaty[maxnode];
91typedef rnodety dnonoty[maxnode];
92typedef boolean umbty[maxsp + 1];
93typedef long lspty[maxsp + 1];
94typedef long longintty[6];
95
96typedef struct tree {
97  /* tree topology data */
98  arynodepty brnchp;   /* point to node                */
99  double lklhd;   /* ln likelihood of this tree   */
100  double vrilkl;   /* variance of likelihood       */
101  lengthty vrilnga;   /* variance of length(branch)   */
102  lengthty vrilnga2;   /* variance2 of length(branch)  */
103  lengthty blklhd;   /* ln likelihood(branch)        */
104  struct node *startp;   /* point to a basal node        */
105  double aic;   /* Akaike Information Criterion */
106  /* convergence of branch length */
107  boolean cnvrgnc;
108} tree;
109
110typedef struct node {
111  /* a node of tree */
112  struct node *isop;   /* pointer to next subnode           */
113  struct node *kinp;   /* pointer to an ascendant node      */
114  long diverg;   /* number of divergences             */
115  long number;   /* node number                       */
116  namety namesp;   /* species name. if this node is tip */
117  spity descen;   /* number of descenant nodes         */
118  nacidty nadata;
119  union {
120    struct {
121      probamity prba;   /* a partial likelihood */
122      double lnga;   /* branch length        */
123    } U0;
124  } UU;
125} node;
126
127
128Static FILE *seqfile;   /* data file of amino acid sequences    */
129Static FILE *lklfile;   /* output file of ln likelihood of site */
130/* transition probability matrix data   */
131Static FILE *tpmfile;
132Static long numsp;   /* number of species                   */
133Static long ibrnch1;   /* first numbering of internal branch  */
134Static long ibrnch2;   /* last numbering of internal branch   */
135Static long numpair;   /* number of species pairs             */
136Static long endsite;   /* number of sites                     */
137Static long endptrn;   /* max number of site patterns         */
138Static long numnw;   /* curent numbering of Newton method   */
139Static long numsm;   /* curent numbering of smoothing       */
140Static long numexe;   /* curent numbering of executes        */
141Static long numtrees;   /* total number of tree topologies     */
142Static long notree;   /* current numbering of tree           */
143Static long maxltr;   /* numbering of max ln likelihood tree */
144Static long minatr;   /* numbering of min AIC tree           */
145/* numbering of decomposable stage     */
146Static long stage;
147Static double maxlkl;   /* max ln likelihood of trees */
148/* min AIC of trees           */
149Static double minaic;
150/* current character of seqfile */
151Static Char chs;
152Static boolean normal;   /* break out error                   */
153Static boolean usertree;   /* U option  designate user trees    */
154Static boolean semiaoptn;   /* S option  semi-auto decomposition */
155Static boolean bootsoptn;   /* B option  bootstrap probability   */
156Static boolean writeoptn;   /* W option  print output data       */
157Static boolean debugoptn;   /* D option  print debug data        */
158Static boolean putlkoptn;   /* P option */
159Static boolean firstoptn;   /* F option */
160Static boolean lastoptn;   /* L option */
161Static boolean readtoptn;   /* R option */
162/*          */
163Static boolean smoothed;
164
165Static tree ctree;   /* current tree                         */
166Static aryamity freqa;   /* amino acid frequency(acid type)      */
167Static charseqty chsequen;   /* acid sequence data(species,site)     */
168Static sitety weight;   /* total number of a new site(new site) */
169Static sitety alias;   /* number of old site(new site)         */
170Static sitety weightw;   /* work area of weight                  */
171Static aryamity freqdyhf, eval, eval2;
172Static daryamity ev, iev;
173Static taryamity coefp;
174Static node *freenode;
175Static itreety paratree;   /* number of parameters         */
176Static rtreety lklhdtree;   /* ln likelihood                */
177Static rtreety lklboottree;   /* ln likelihood (bootstrap)    */
178Static rtreety aictree;   /* Akaike Information Criterion */
179Static rtreety boottree;   /* bootstrap probability        */
180Static ltrptty lklhdtrpt;   /* ln likelihood                */
181
182
183/*** print line ***/
184Static Void printline(num)
185long num;
186{
187  /* print '-' line to standard output */
188  long i;
189
190  putchar(' ');
191  for (i = 1; i <= num; i++)
192    putchar('-');
193  putchar('\n');
194}  /* printline */
195
196
197/*******************************************************/
198/*****    READ DATA, INITIALIZATION AND SET UP     *****/
199/*******************************************************/
200
201/*** GET NUMBERS OF SPECIES AND SITES ***/
202Static Void getnums()
203{
204  /* get species-number,sites-number from file*/
205  /* input number of species, number of sites */
206  fscanf(seqfile, "%ld%ld", &numsp, &endsite);
207  if (usertree) {
208    ibrnch1 = numsp + 1;
209    ibrnch2 = numsp * 2 - 3;
210  } else {
211    ibrnch1 = numsp + 1;
212    ibrnch2 = numsp;
213  }
214  /* numpair  := numsp*(numsp-1) DIV 2; */
215  numpair = (long)((numsp * numsp - numsp) / 2.0);
216  if (numsp > maxsp)
217    printf("TOO MANY SPECIES: adjust CONSTants\n");
218  if (endsite > maxsite)
219    printf("TOO MANY SITES:   adjust CONSTants\n");
220  normal = (numsp <= maxsp && endsite <= maxsite);
221}  /* getnums */
222
223
224Local boolean letter(chru)
225Char chru;
226{
227  /* tests whether chru is a letter, (ASCII and EBCDIC) */
228  return (chru >= 'A' && chru <= 'I' || chru >= 'J' && chru <= 'R' ||
229          chru >= 'S' && chru <= 'Z' || chru >= 'a' && chru <= 'i' ||
230          chru >= 'j' && chru <= 'r' || chru >= 's' && chru <= 'z');
231}  /* letter */
232
233
234/*** CONVERT CHARACTER TO UPPER CASE ***/
235Static Void uppercase(chru)
236Char *chru;
237{
238  /* convert character to upper case */
239  /* convert chru to upper case -- either ASCII or EBCDIC */
240  if (letter(*chru) &&
241      (strcmp("a", "A") > 0 && *chru >= 'a' ||
242       strcmp("a", "A") < 0 && *chru < 'A'))
243    *chru = _toupper(*chru);
244}  /* uppercase */
245
246
247/*** GET OPERATIONAL OPTIONS ***/
248Static Void getoptions()
249{
250  /* get option from sequences file */
251  Char chrop;
252
253  usertree = false;
254  semiaoptn = false;
255  firstoptn = false;
256  lastoptn = false;
257  putlkoptn = false;
258  bootsoptn = false;
259  writeoptn = false;
260  debugoptn = false;
261  readtoptn = false;
262  while (!P_eoln(seqfile)) {
263    chrop = getc(seqfile);
264    if (chrop == '\n')
265      chrop = ' ';
266    uppercase(&chrop);
267    if (chrop != 'U' && chrop != 'S' && chrop != 'F' && chrop != 'L' &&
268        chrop != 'B' && chrop != 'P' && chrop != 'W' && chrop != 'D' &&
269        chrop != 'R') {
270      if (chrop != ' ') {
271        printf(" BAD OPTION CHARACTER: %c\n", chrop);
272        normal = false;
273      }
274      continue;
275    }
276    switch (chrop) {
277
278    case 'U':
279      usertree = true;
280      break;
281
282    case 'S':
283      semiaoptn = true;
284      break;
285
286    case 'F':
287      firstoptn = true;
288      break;
289
290    case 'L':
291      lastoptn = true;
292      break;
293
294    case 'P':
295      putlkoptn = true;
296      break;
297
298    case 'B':
299      bootsoptn = true;
300      break;
301
302    case 'W':
303      writeoptn = true;
304      break;
305
306    case 'D':
307      debugoptn = true;
308      break;
309
310    case 'R':
311      readtoptn = true;
312      break;
313    }
314  }
315  if (!(numexe == 1 && writeoptn))
316    return;
317  printf("\n%15s", "Users Option");
318  printf("%14s%5s", "Usertree  : ", usertree ? "TRUE" : "FALSE");
319  printf("%14s%5s", "Semiaoptn : ", semiaoptn ? "TRUE" : "FALSE");
320  printf("%14s%5s\n", "Bootsoptn : ", bootsoptn ? "TRUE" : "FALSE");
321  printf("%15c", ' ');
322  printf("%14s%5s", "Putlkoptn : ", putlkoptn ? "TRUE" : "FALSE");
323  printf("%14s%5s", "Writeoptn : ", writeoptn ? "TRUE" : "FALSE");
324  printf("%14s%5s\n", "Debugoptn : ", debugoptn ? "TRUE" : "FALSE");
325  printf("%15c", ' ');
326  printf("%14s%5s", "Firstoptn : ", firstoptn ? "TRUE" : "FALSE");
327  printf("%14s%5s\n\n", "Lastoptn  : ", lastoptn ? "TRUE" : "FALSE");
328}  /* getoptions */
329
330
331/*** GET AMINO ACID FREQENCY FROM DATASET ***/
332Static Void readbasefreqs()
333{
334  /* get freqency from sequences file */
335  amity ba;
336
337  fscanf(seqfile, "%*[^\n]");
338  getc(seqfile);
339  for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
340    fscanf(seqfile, "%lg", &freqa[(long)ba - (long)AA]);
341}  /* readbasefreqs */
342
343
344/*** SETUP DATASTRUCTURE OF TREE ***/
345Static Void setuptree(tr)
346tree *tr;
347{
348  /* setup datastructure of tree */
349  long i, k;
350  node *p;
351  long FORLIM, FORLIM1;
352
353  FORLIM = numsp;
354  for (i = 1; i <= FORLIM; i++) {
355    p = (node *)Malloc(sizeof(node));
356    p->number = i;
357    p->UU.U0.lnga = 0.0;
358    p->isop = NULL;
359    p->diverg = 0;
360    FORLIM1 = numsp;
361    for (k = 0; k < FORLIM1; k++)
362      p->descen[k] = 0;
363    p->descen[i - 1] = 1;
364    tr->brnchp[i] = p;
365    p->nadata = ami;
366  }
367  tr->lklhd = -999999.0;
368  tr->aic = 0.0;
369  tr->startp = tr->brnchp[1];
370  tr->cnvrgnc = false;
371}  /* setuptree */
372
373
374/*** GET SEQUENCES AND PRINTOUT SEQUENCES ***/
375Static Void getsequence()
376{
377  /* get and print sequences from file */
378  long is, js, ks, ls, ms, ns, ichr, numchr, maxchr;
379  Char chrs;
380  Char chrcnt[30];
381  long ichrcnt[30];
382  boolean namechar, existenced;
383  long FORLIM, FORLIM1;
384
385  fscanf(seqfile, "%*[^\n]");
386  getc(seqfile);
387  if (debugoptn)
388    putchar('\n');
389  FORLIM = numsp;
390  for (is = 1; is <= FORLIM; is++) {
391    /* readln(seqfile); */
392    do {
393      if (P_eoln(seqfile)) {
394        fscanf(seqfile, "%*[^\n]");
395        getc(seqfile);
396      }
397      chrs = getc(seqfile);
398      if (chrs == '\n')
399        chrs = ' ';
400      if (debugoptn) {
401        if (chrs == ' ')
402          putchar(chrs);
403      }
404    } while (chrs == ' ');
405    if (debugoptn) {
406      printf("<SPACE>\n");
407      putchar(chrs);
408    }
409    namechar = true;
410    ns = 0;
411    do {
412      ns++;
413      if (ns <= maxname)
414        ctree.brnchp[is]->namesp[ns - 1] = chrs;
415      if (P_eoln(seqfile))
416        namechar = false;
417      else {
418        chrs = getc(seqfile);
419        if (chrs == '\n')
420          chrs = ' ';
421        if (chrs == ' ')
422          namechar = false;
423      }
424      if (debugoptn) {
425        if (namechar)
426          putchar(chrs);
427      }
428    } while (namechar);
429    if (ns < maxname) {
430      for (js = ns; js < maxname; js++)
431        ctree.brnchp[is]->namesp[js] = ' ';
432    }
433    /*  FOR js := 1 TO maxname DO
434        BEGIN
435           read(seqfile, chrs);
436           ctree.brnchp[is]^.namesp[js] := chrs;
437        END;  */
438    if (debugoptn) {
439      printf("<NAME>\n");
440      putchar(chrs);
441    }
442    FORLIM1 = endsite;
443    for (js = 1; js <= FORLIM1; js++) {
444      if (normal) {
445        do {
446          if (P_eoln(seqfile)) {
447            fscanf(seqfile, "%*[^\n]");
448            getc(seqfile);
449          }
450          chrs = getc(seqfile);
451          if (chrs == '\n')
452            chrs = ' ';
453        } while (chrs == ' ' || chrs >= '0' && chrs <= '9');
454        uppercase(&chrs);
455        if (debugoptn)
456          putchar(chrs);
457        if (chrs != 'A' && chrs != 'C' && chrs != 'D' && chrs != 'E' &&
458            chrs != 'F' && chrs != 'G' && chrs != 'H' && chrs != 'I' &&
459            chrs != 'K' && chrs != 'L' && chrs != 'M' && chrs != 'N' &&
460            chrs != 'P' && chrs != 'Q' && chrs != 'R' && chrs != 'S' &&
461            chrs != 'T' && chrs != 'V' && chrs != 'W' && chrs != 'Y' &&
462            chrs != 'B' && chrs != 'Z' && chrs != 'V' && chrs != 'X' &&
463            chrs != '*' && chrs != '-') {
464          printf(
465            "\n WARNING -- BAD AMINO ACID \"%c\" AT POSITION%5ld OF SPECIES%3ld\n",
466            chrs, js, is);
467          normal = false;
468        }
469        chsequen[is][js - 1] = chrs;
470      }
471    }
472    if (debugoptn)
473      putchar('\n');
474  }
475  fscanf(seqfile, "%*[^\n]");
476  getc(seqfile);
477
478  FORLIM = endsite;
479  for (ks = 0; ks < FORLIM; ks++) {
480    ichr = 0;
481    FORLIM1 = numsp;
482    for (js = 0; js < FORLIM1; js++)
483      ichrcnt[js] = 0;
484    FORLIM1 = numsp;
485    for (js = 0; js < FORLIM1; js++)
486      chrcnt[js] = ' ';
487    FORLIM1 = numsp;
488    for (js = 1; js <= FORLIM1; js++) {
489      chrs = chsequen[js][ks];
490      existenced = false;
491      for (ms = 0; ms < ichr; ms++) {
492        if (chrs == chrcnt[ms]) {
493          ichrcnt[ms]++;
494          existenced = true;
495        }
496      }
497      if (!existenced) {
498        ichr++;
499        chrcnt[ichr - 1] = chrs;
500        ichrcnt[ichr - 1] = 1;
501      }
502    }
503    maxchr = 0;
504    numchr = 0;
505    for (ms = 1; ms <= ichr; ms++) {
506      if (ichrcnt[ms - 1] > maxchr) {
507        maxchr = ichrcnt[ms - 1];
508        numchr = ms;
509      }
510    }
511    if (numchr != 0)
512      chsequen[0][ks] = chrcnt[numchr - 1];
513    else
514      chsequen[0][ks] = '?';
515  }
516
517  if (numexe != 1)
518    return;
519  printf("\n%15s", "Sequences Data");
520  printf("%5ld%9s%5ld%6s\n\n", numsp, "Species,", endsite, "Sites");
521  printf(" Species%5cAmino Acid Sequences\n", ' ');
522  printf(" -------%5c--------------------\n\n", ' ');
523  FORLIM = (endsite - 1) / maxline + 1;
524  for (is = 1; is <= FORLIM; is++) {
525    FORLIM1 = numsp;
526    for (js = 0; js <= FORLIM1; js++) {
527      if (js == 0) {
528        printf(" Consensus");   /* Consensus Consent */
529        for (ks = 1; ks <= maxname - 7; ks++)
530          putchar(' ');
531
532      } else {
533        putchar(' ');
534        for (ks = 0; ks < maxname; ks++)
535          putchar(ctree.brnchp[js]->namesp[ks]);
536        printf("  ");
537      }
538      ls = maxline * is;
539      if (ls > endsite)
540        ls = endsite;
541      for (ks = maxline * (is - 1) + 1; ks <= ls; ks++) {
542        if (normal) {
543          if (js > 0 && chsequen[js][ks - 1] == chsequen[0][ks - 1])
544            putchar('.');
545          else
546            putchar(chsequen[js][ks - 1]);
547          if (ks % 10 == 0 && ks % maxline != 0)
548            putchar(' ');
549/* p2c: protml.pas, line 466:
550 * Note: Using % for possibly-negative arguments [317] */
551/* p2c: protml.pas, line 466:
552 * Note: Using % for possibly-negative arguments [317] */
553        }
554      }
555      putchar('\n');
556    }
557    putchar('\n');
558  }
559}  /* getsequence */
560
561
562/*************************************/
563/***      READ SEQUENCES DATA      ***/
564/*************************************/
565
566Static Void getdata()
567{
568  /* read data from sequences file */
569  if (numexe == 1) {
570    printline(57L);
571    printf("              PROTML Ver.1.00b  in MOLPHY\n");
572    printf("    Maximum Likelihood Inference of Protein Phylogeny\n");
573    printf("                based on Dayhoff model\n");
574    printline(57L);
575  }
576  getnums();
577  if (normal)
578    getoptions();
579  /* IF NOT freqsfrom THEN readbasefreqs; */
580  if (numexe == 1) {
581    if (normal)
582      setuptree(&ctree);
583  }
584  if (normal)
585    getsequence();
586}  /* getdata */
587
588
589/*** SORT OF SEQUENCES ***/
590Static Void sitesort()
591{
592  /* sort sequences */
593  /* Shell sort keeping sites, weights in same order */
594  long gap, i, j, jj, jg, k, itemp;
595  boolean flip, tied;
596  long FORLIM;
597
598  FORLIM = endsite;
599  for (i = 1; i <= FORLIM; i++) {
600    alias[i - 1] = i;
601    weightw[i - 1] = 1;
602  }
603  gap = endsite / 2;
604  while (gap > 0) {
605    FORLIM = endsite;
606    for (i = gap + 1; i <= FORLIM; i++) {
607      j = i - gap;
608      flip = true;
609      while (j > 0 && flip) {
610        jj = alias[j - 1];
611        jg = alias[j + gap - 1];
612        flip = false;   /* fel */
613        tied = true;   /* fel */
614        k = 1;
615        while (k <= numsp && tied) {
616          flip = (chsequen[k][jj - 1] > chsequen[k][jg - 1]);
617          tied = (tied && chsequen[k][jj - 1] == chsequen[k][jg - 1]);
618          k++;
619        }
620        if (!flip)
621          break;
622        itemp = alias[j - 1];
623        alias[j - 1] = alias[j + gap - 1];
624        alias[j + gap - 1] = itemp;
625        itemp = weightw[j - 1];
626        weightw[j - 1] = weightw[j + gap - 1];
627        weightw[j + gap - 1] = itemp;
628        j -= gap;
629      }
630    }
631    gap /= 2;
632  }
633}  /* sitesort */
634
635
636/*** COMBINATION OF SITE ***/
637Static Void sitecombine()
638{
639  /* combine sites */
640  /* combine sites that have identical patterns */
641  long i, j, k;
642  boolean tied;
643
644  i = 1;
645  while (i < endsite) {
646    j = i + 1;
647    tied = true;
648    while (j <= endsite && tied) {
649      k = 1;   /* fel */
650      while (k <= numsp && tied) {
651        tied = (tied && chsequen[k][alias[i - 1] - 1] == chsequen[k]
652                        [alias[j - 1] - 1]);
653        k++;
654      }
655      if (tied && weightw[j - 1] > 0) {
656        weightw[i - 1] += weightw[j - 1];
657        weightw[j - 1] = 0;
658        alias[j - 1] = alias[i - 1];
659      }
660      j++;
661    }
662    i = j - 1;
663  }
664}  /* sitecombine */
665
666
667/*** SCRUNCH OF SITE ***/
668Static Void sitescrunch()
669{
670  /* move so positively weighted sites come first */
671  long i, j, itemp;
672  boolean done, found;
673
674  done = false;
675  i = 1;
676  j = 2;
677  while (!done) {
678    found = false;
679    if (weightw[i - 1] > 0)
680      i++;
681    else {
682      if (j <= i)
683        j = i + 1;
684      if (j <= endsite) {
685        found = false;
686        do {
687          found = (weightw[j - 1] > 0);
688          j++;
689        } while (!(found || j > endsite));
690        if (found) {
691          j--;
692          itemp = alias[i - 1];
693          alias[i - 1] = alias[j - 1];
694          alias[j - 1] = itemp;
695          itemp = weightw[i - 1];
696          weightw[i - 1] = weightw[j - 1];
697          weightw[j - 1] = itemp;
698        } else
699          done = true;
700      } else
701        done = true;
702    }
703    done = (done || i >= endsite);
704  }
705}  /* sitescrunch */
706
707
708/*** PRINT PATTERN ***/
709Static Void printpattern(maxorder)
710long maxorder;
711{
712  /* print patterned sequences */
713  /* print pattern */
714  long i, j, k, l, n, m, big, sml, kweight, FORLIM, FORLIM1;
715
716  /* */
717  if (debugoptn) {
718    printf("\nalias  :\n");
719    FORLIM = endptrn;
720    for (i = 0; i < FORLIM; i++)
721      printf("%4ld", alias[i]);
722    printf("\n\nweight :\n");
723    FORLIM = endptrn;
724    for (i = 0; i < FORLIM; i++)
725      printf("%4ld", weight[i]);
726    printf("\n\nendptrn =%5ld\n\n\n", endptrn);
727  }
728  if (debugoptn) {
729    printf("  num alias weight  pattern\n");
730    FORLIM = endptrn;
731    for (i = 0; i < FORLIM; i++) {
732      printf("%4ld%6ld%7ld   ", i + 1, alias[i], weight[i]);
733      FORLIM1 = numsp;
734      for (j = 1; j <= FORLIM1; j++)
735        printf("%2c", chsequen[j][alias[i] - 1]);
736      putchar('\n');
737    }
738    putchar('\n');
739  }
740  /* */
741  if (!writeoptn)
742    return;
743
744  printf("\n Species%5cPatternized Sequences\n", ' ');
745  printf(" -------%5c---------------------\n\n", ' ');
746  FORLIM = (endptrn - 1) / maxline;
747  for (i = 0; i <= FORLIM; i++) {
748    l = maxline * (i + 1);
749    if (l > endptrn)
750      l = endptrn;
751    FORLIM1 = numsp;
752    for (j = 1; j <= FORLIM1; j++) {
753      putchar(' ');
754      for (k = 0; k < maxname; k++)
755        putchar(ctree.brnchp[j]->namesp[k]);
756      printf("  ");
757      for (k = maxline * i + 1; k <= l; k++) {
758        if (normal) {
759          if (j > 1 &&
760              chsequen[j][alias[k - 1] - 1] == chsequen[1][alias[k - 1] - 1])
761            putchar('.');
762          else
763            putchar(chsequen[j][alias[k - 1] - 1]);
764          if (k % 10 == 0 && k % maxline != 0)
765            putchar(' ');
766/* p2c: protml.pas, line 685:
767 * Note: Using % for possibly-negative arguments [317] */
768/* p2c: protml.pas, line 685:
769 * Note: Using % for possibly-negative arguments [317] */
770        }
771      }
772      putchar('\n');
773    }
774    putchar('\n');
775
776    for (n = maxorder; n >= 1; n--) {
777      printf("%3c", ' ');
778      for (k = 1; k <= maxname; k++)
779        putchar(' ');
780      big = 1;
781      for (m = 1; m <= n; m++)
782        big *= 10;
783      sml = big / 10;
784      for (k = maxline * i + 1; k <= l; k++) {
785        if (normal) {
786          kweight = weight[k - 1] % big / sml;
787/* p2c: protml.pas, line 700:
788 * Note: Using % for possibly-negative arguments [317] */
789          if (kweight > 0 && kweight < 10)
790            printf("%ld", kweight);
791          else if (kweight == 0) {
792            if (weight[k - 1] % big == weight[k - 1])
793              putchar(' ');
794            else if (weight[k - 1] % big < weight[k - 1])
795              putchar('0');
796            else
797              putchar('*');
798/* p2c: protml.pas, line 705:
799 * Note: Using % for possibly-negative arguments [317] */
800          } else
801            putchar('?');
802          if (k % 10 == 0 && k % maxline != 0)
803            putchar(' ');
804/* p2c: protml.pas, line 714:
805 * Note: Using % for possibly-negative arguments [317] */
806/* p2c: protml.pas, line 714:
807 * Note: Using % for possibly-negative arguments [317] */
808        }
809      }
810      putchar('\n');
811    }
812    putchar('\n');
813  }
814
815/* p2c: protml.pas, line 707:
816 * Note: Using % for possibly-negative arguments [317] */
817}  /* printpattern */
818
819
820/***********************************/
821/***  ARRANGE SITES OF SEQUENCES  ***/
822/***********************************/
823
824Static Void makeweights()
825{
826  /* condense same site-pattern */
827  /* make up weights vector to avoid duplicate computations */
828  long iw, jw, kw, maxorder, maxweight, nweight, nw, FORLIM, FORLIM1;
829
830  sitesort();
831  sitecombine();
832  sitescrunch();
833  maxorder = 0;
834  maxweight = 0;
835  FORLIM = endsite;
836  for (iw = 0; iw < FORLIM; iw++) {
837    weight[iw] = weightw[iw];
838    if (weight[iw] > 0)
839      endptrn = iw + 1;
840    if (weight[iw] > maxweight) {
841      maxweight = weight[iw];
842      nweight = weight[iw];
843      nw = 0;
844      do {
845        nweight /= 10;
846        nw++;
847      } while (nweight != 0);
848      if (nw > maxorder)
849        maxorder = nw;
850    }
851  }
852  if (endptrn > maxptrn) {
853    printf(" TOO MANY PATTERNS: increase CONSTant maxptrn to at least %5ld\n",
854           endptrn);
855    normal = false;
856  }
857
858  kw = 0;
859  FORLIM = endptrn;
860  for (iw = 1; iw <= FORLIM; iw++) {
861    FORLIM1 = weight[iw - 1];
862    for (jw = 1; jw <= FORLIM1; jw++) {
863      kw++;
864      weightw[kw - 1] = iw;
865    }
866  }
867
868  printpattern(maxorder);
869
870}  /* makeweights */
871
872
873/*****************************************/
874/***  SET PARTIAL LIKELIHOODS AT TIPS  ***/
875/*****************************************/
876
877Static Void makevalues()
878{
879  /* set up fractional likelihoods */
880  /* set up fractional likelihoods at tips */
881  long i, j, k;
882  amity ba;
883  long FORLIM, FORLIM1;
884
885  FORLIM = endptrn;
886  for (k = 0; k < FORLIM; k++) {
887    j = alias[k];
888    FORLIM1 = numsp;
889    for (i = 1; i <= FORLIM1; i++) {
890      for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
891        ctree.brnchp[i]->UU.U0.prba[k][(long)ba - (long)AA] = 0.0;
892      switch (chsequen[i][j - 1]) {
893
894      case 'A':
895        ctree.brnchp[i]->UU.U0.prba[k][0] = 1.0;
896        break;
897
898      case 'R':
899        ctree.brnchp[i]->UU.U0.prba[k][(long)RR - (long)AA] = 1.0;
900        break;
901
902      case 'N':
903        ctree.brnchp[i]->UU.U0.prba[k][(long)NN - (long)AA] = 1.0;
904        break;
905
906      case 'D':
907        ctree.brnchp[i]->UU.U0.prba[k][(long)DD - (long)AA] = 1.0;
908        break;
909
910      case 'C':
911        ctree.brnchp[i]->UU.U0.prba[k][(long)CC - (long)AA] = 1.0;
912        break;
913
914      case 'Q':
915        ctree.brnchp[i]->UU.U0.prba[k][(long)QQ - (long)AA] = 1.0;
916        break;
917
918      case 'E':
919        ctree.brnchp[i]->UU.U0.prba[k][(long)EE - (long)AA] = 1.0;
920        break;
921
922      case 'G':
923        ctree.brnchp[i]->UU.U0.prba[k][(long)GG - (long)AA] = 1.0;
924        break;
925
926      case 'H':
927        ctree.brnchp[i]->UU.U0.prba[k][(long)HH - (long)AA] = 1.0;
928        break;
929
930      case 'I':
931        ctree.brnchp[i]->UU.U0.prba[k][(long)II - (long)AA] = 1.0;
932        break;
933
934      case 'L':
935        ctree.brnchp[i]->UU.U0.prba[k][(long)LL - (long)AA] = 1.0;
936        break;
937
938      case 'K':
939        ctree.brnchp[i]->UU.U0.prba[k][(long)KK - (long)AA] = 1.0;
940        break;
941
942      case 'M':
943        ctree.brnchp[i]->UU.U0.prba[k][(long)MM - (long)AA] = 1.0;
944        break;
945
946      case 'F':
947        ctree.brnchp[i]->UU.U0.prba[k][(long)FF - (long)AA] = 1.0;
948        break;
949
950      case 'P':
951        ctree.brnchp[i]->UU.U0.prba[k][(long)PP_ - (long)AA] = 1.0;
952        break;
953
954      case 'S':
955        ctree.brnchp[i]->UU.U0.prba[k][(long)SS - (long)AA] = 1.0;
956        break;
957
958      case 'T':
959        ctree.brnchp[i]->UU.U0.prba[k][(long)TT - (long)AA] = 1.0;
960        break;
961
962      case 'W':
963        ctree.brnchp[i]->UU.U0.prba[k][(long)WW - (long)AA] = 1.0;
964        break;
965
966      case 'Y':
967        ctree.brnchp[i]->UU.U0.prba[k][(long)YY - (long)AA] = 1.0;
968        break;
969
970      case 'V':
971        ctree.brnchp[i]->UU.U0.prba[k][(long)VV - (long)AA] = 1.0;
972        break;
973
974      case 'B':
975        ctree.brnchp[i]->UU.U0.prba[k][(long)DD - (long)AA] = 0.5;
976        ctree.brnchp[i]->UU.U0.prba[k][(long)NN - (long)AA] = 0.5;
977        break;
978
979      case 'Z':
980        ctree.brnchp[i]->UU.U0.prba[k][(long)EE - (long)AA] = 0.5;
981        ctree.brnchp[i]->UU.U0.prba[k][(long)QQ - (long)AA] = 0.5;
982        break;
983
984      case 'X':
985        for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
986          ctree.brnchp[i]->UU.U0.prba[k][(long)ba - (long)AA] = 1.0;
987        break;
988
989      case '?':
990        for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
991          ctree.brnchp[i]->UU.U0.prba[k][(long)ba - (long)AA] = 1.0;
992        break;
993
994      case '-':
995        for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
996          ctree.brnchp[i]->UU.U0.prba[k][(long)ba - (long)AA] = 1.0;
997        break;
998      }
999    }
1000  }
1001}  /* makevalues */
1002
1003
1004/*** EMPIRICAL FREQENCIES OF AMINO ACIDS ***/
1005Static Void empirifreqsA()
1006{
1007  /* calculate empirical frequencies */
1008  /* Get empirical frequencies of amino acids from the data */
1009  long ia, ja;
1010  amity ba;
1011  aryamity sfreqa;
1012  double sum;
1013  long FORLIM, FORLIM1;
1014
1015  for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
1016    sfreqa[(long)ba - (long)AA] = 0.0;
1017  sum = 0.0;
1018  FORLIM = numsp;
1019  for (ia = 1; ia <= FORLIM; ia++) {
1020    FORLIM1 = endptrn;
1021    for (ja = 0; ja < FORLIM1; ja++) {
1022      for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
1023        sfreqa[(long)ba - (long)AA] +=
1024          weight[ja] * ctree.brnchp[ia]->UU.U0.prba[ja][(long)ba - (long)AA];
1025    }
1026  }
1027  for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
1028    sum += sfreqa[(long)ba - (long)AA];
1029  for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
1030    freqa[(long)ba - (long)AA] = sfreqa[(long)ba - (long)AA] / sum;
1031  if (!(numexe == 1 && writeoptn))
1032    return;
1033  printf("\n%13s%7.1f\n", " Total acid :", sum);
1034  printf("%13s", "  A,R,N,D,C :");
1035  for (ba = AA; (long)ba <= (long)CC; ba = (amity)((long)ba + 1))
1036    printf("%7.1f", sfreqa[(long)ba - (long)AA]);
1037  printf("\n%13s", "  Q,E,G,H,I :");
1038  for (ba = QQ; (long)ba <= (long)II; ba = (amity)((long)ba + 1))
1039    printf("%7.1f", sfreqa[(long)ba - (long)AA]);
1040  printf("\n%13s", "  L,K,M,F,P :");
1041  for (ba = LL; (long)ba <= (long)PP_; ba = (amity)((long)ba + 1))
1042    printf("%7.1f", sfreqa[(long)ba - (long)AA]);
1043  printf("\n%13s", "  S,T,W,Y,V :");
1044  for (ba = SS; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
1045    printf("%7.1f", sfreqa[(long)ba - (long)AA]);
1046  putchar('\n');
1047}  /* empirifreqsA */
1048
1049
1050/*************************************/
1051/***      GET FREQUENCY            ***/
1052/*************************************/
1053
1054Static Void getbasefreqs()
1055{
1056  /* print frequencies */
1057  amity ba;
1058
1059  empirifreqsA();
1060  if (!(numexe == 1 && writeoptn))
1061    return;
1062  printf("\n Empirical");
1063  /* IF freqsfrom THEN */
1064  printf(" Amino Acid Frequencies:\n");
1065  printf("%13s", "A,R,N,D,C :");
1066  for (ba = AA; (long)ba <= (long)CC; ba = (amity)((long)ba + 1))
1067    printf("%7.3f", freqa[(long)ba - (long)AA]);
1068  printf("\n%13s", "Q,E,G,H,I :");
1069  for (ba = QQ; (long)ba <= (long)II; ba = (amity)((long)ba + 1))
1070    printf("%7.3f", freqa[(long)ba - (long)AA]);
1071  printf("\n%13s", "L,K,M,F,P :");
1072  for (ba = LL; (long)ba <= (long)PP_; ba = (amity)((long)ba + 1))
1073    printf("%7.3f", freqa[(long)ba - (long)AA]);
1074  printf("\n%13s", "S,T,W,Y,V :");
1075  for (ba = SS; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
1076    printf("%7.3f", freqa[(long)ba - (long)AA]);
1077  putchar('\n');
1078}  /* getbasefreqs */
1079
1080
1081Static Void getinput()
1082{
1083  /* read data and set up data */
1084  /* read input file */
1085  normal = true;
1086  if (normal)
1087    getdata();
1088  if (normal)
1089    makeweights();
1090  if (normal)
1091    makevalues();
1092  if (normal)
1093    getbasefreqs();
1094}  /* getinput */
1095
1096
1097Static Void printnmtrx(nmtrx)
1098naryamity *nmtrx;
1099{
1100  long i, j;
1101
1102  printf(" MATRIX(numtype)\n");
1103  for (i = 0; i <= 19; i++) {
1104    for (j = 1; j <= 20; j++) {
1105      printf("%7.3f", nmtrx[i][j - 1]);
1106      if (j == 10 || j == 20)
1107        putchar('\n');
1108    }
1109    putchar('\n');
1110  }
1111}  /* printnmtrx */
1112
1113
1114#define eps             1.0e-20
1115
1116
1117Static Void luinverse(omtrx, imtrx, nmtrx)
1118naryamity *omtrx, *imtrx;
1119long nmtrx;
1120{
1121  /* INVERSION OF MATRIX ON LU DECOMPOSITION */
1122  long i, j, k, l, maxi, idx, ix, jx;
1123  double sum, tmp, maxb, aw;
1124  niaryamity index;
1125  double *wk;
1126
1127  wk = (double *)Malloc(sizeof(naryamity));
1128  aw = 1.0;
1129  for (i = 0; i < nmtrx; i++) {
1130    maxb = 0.0;
1131    for (j = 0; j < nmtrx; j++) {
1132      if (fabs(omtrx[i][j]) > maxb)
1133        maxb = fabs(omtrx[i][j]);
1134    }
1135    if (maxb == 0.0)
1136      printf("PROC. LUINVERSE:  singular\n");
1137    wk[i] = 1.0 / maxb;
1138  }
1139  for (j = 1; j <= nmtrx; j++) {
1140    for (i = 1; i < j; i++) {
1141      sum = omtrx[i - 1][j - 1];
1142      for (k = 0; k <= i - 2; k++)
1143        sum -= omtrx[i - 1][k] * omtrx[k][j - 1];
1144      omtrx[i - 1][j - 1] = sum;
1145    }
1146    maxb = 0.0;
1147    for (i = j; i <= nmtrx; i++) {
1148      sum = omtrx[i - 1][j - 1];
1149      for (k = 0; k <= j - 2; k++)
1150        sum -= omtrx[i - 1][k] * omtrx[k][j - 1];
1151      omtrx[i - 1][j - 1] = sum;
1152      tmp = wk[i - 1] * fabs(sum);
1153      if (tmp >= maxb) {
1154        maxb = tmp;
1155        maxi = i;
1156      }
1157    }
1158    if (j != maxi) {
1159      for (k = 0; k < nmtrx; k++) {
1160        tmp = omtrx[maxi - 1][k];
1161        omtrx[maxi - 1][k] = omtrx[j - 1][k];
1162        omtrx[j - 1][k] = tmp;
1163      }
1164      aw = -aw;
1165      wk[maxi - 1] = wk[j - 1];
1166    }
1167    index[j - 1] = maxi;
1168    if (omtrx[j - 1][j - 1] == 0.0)
1169      omtrx[j - 1][j - 1] = eps;
1170    if (j != nmtrx) {
1171      tmp = 1.0 / omtrx[j - 1][j - 1];
1172      for (i = j; i < nmtrx; i++)
1173        omtrx[i][j - 1] *= tmp;
1174    }
1175  }
1176  for (jx = 0; jx < nmtrx; jx++) {
1177    for (ix = 0; ix < nmtrx; ix++)
1178      wk[ix] = 0.0;
1179    wk[jx] = 1.0;
1180    l = 0;
1181    for (i = 1; i <= nmtrx; i++) {
1182      idx = index[i - 1];
1183      sum = wk[idx - 1];
1184      wk[idx - 1] = wk[i - 1];
1185      if (l != 0) {
1186        for (j = l - 1; j <= i - 2; j++)
1187          sum -= omtrx[i - 1][j] * wk[j];
1188      } else if (sum != 0.0)
1189        l = i;
1190      wk[i - 1] = sum;
1191    }
1192    for (i = nmtrx - 1; i >= 0; i--) {
1193      sum = wk[i];
1194      for (j = i + 1; j < nmtrx; j++)
1195        sum -= omtrx[i][j] * wk[j];
1196      wk[i] = sum / omtrx[i][i];
1197    }
1198    for (ix = 0; ix < nmtrx; ix++)
1199      imtrx[ix][jx] = wk[ix];
1200  }
1201  Free(wk);
1202}  /* luinverse */
1203
1204#undef eps
1205
1206
1207Static Void mproduct(am, bm, cm, na, nb, nc)
1208naryamity *am, *bm, *cm;
1209long na, nb, nc;
1210{
1211  long ia, ib, ic;
1212  double sum;
1213
1214  for (ia = 0; ia < na; ia++) {
1215    for (ic = 0; ic < nc; ic++) {
1216      sum = 0.0;
1217      for (ib = 0; ib < nb; ib++)
1218        sum += am[ia][ib] * bm[ib][ic];
1219      cm[ia][ic] = sum;
1220    }
1221  }
1222}  /* mproduct */
1223
1224
1225Static Void convermtrx(amtrx, nmtrx)
1226aryamity *amtrx;
1227naryamity *nmtrx;
1228{
1229  amity ba1, ba2;
1230
1231  for (ba1 = AA; (long)ba1 <= (long)VV; ba1 = (amity)((long)ba1 + 1)) {
1232    for (ba2 = AA; (long)ba2 <= (long)VV; ba2 = (amity)((long)ba2 + 1))
1233      nmtrx[(int)ba1][(int)ba2] = amtrx[(long)ba1 - (long)AA]
1234        [(long)ba2 - (long)AA];
1235  }
1236}  /* convermtrx */
1237
1238
1239Static Void reconvermtrx(nmtrx, amtrx)
1240naryamity *nmtrx;
1241aryamity *amtrx;
1242{
1243  amity ba1, ba2;
1244
1245  for (ba1 = AA; (long)ba1 <= (long)VV; ba1 = (amity)((long)ba1 + 1)) {
1246    for (ba2 = AA; (long)ba2 <= (long)VV; ba2 = (amity)((long)ba2 + 1))
1247      amtrx[(long)ba1 - (long)AA][(long)ba2 - (long)AA] = nmtrx[(int)ba1]
1248        [(int)ba2];
1249  }
1250}  /* reconvermtrx */
1251
1252
1253Static Void printeigen()
1254{
1255  amity ba1, ba2;
1256
1257  printf(" EIGEN VECTOR\n");
1258  for (ba1 = AA; (long)ba1 <= (long)VV; ba1 = (amity)((long)ba1 + 1)) {
1259    for (ba2 = AA; (long)ba2 <= (long)VV; ba2 = (amity)((long)ba2 + 1)) {
1260      printf("%7.3f", ev[(long)ba1 - (long)AA][(long)ba2 - (long)AA]);
1261      if (ba2 == II || ba2 == VV)
1262        putchar('\n');
1263    }
1264    putchar('\n');
1265  }
1266  printf(" EIGEN VALUE\n");
1267  for (ba1 = AA; (long)ba1 <= (long)VV; ba1 = (amity)((long)ba1 + 1)) {
1268    printf("%7.3f", eval[(long)ba1 - (long)AA]);
1269    if (ba1 == II || ba1 == VV)
1270      putchar('\n');
1271  }
1272}  /* printeigen */
1273
1274
1275Static Void checkevector(imtrx, nn)
1276naryamity *imtrx;
1277long nn;
1278{
1279  long i, j;
1280
1281  for (i = 1; i <= nn; i++) {
1282    for (j = 1; j <= nn; j++) {
1283      if (i == j) {
1284        if (fabs(imtrx[i - 1][j - 1] - 1.0) > 1.0e-5)
1285          printf(" error1 eigen vector (checkevector)\n");
1286      } else {
1287        if (fabs(imtrx[i - 1][j - 1]) > 1.0e-5)
1288          printf(" error2 eigen vector (checkevector)\n");
1289      }
1290    }
1291  }
1292}  /* checkevector */
1293
1294
1295Static Void printamtrx(amtrx)
1296aryamity *amtrx;
1297{
1298  amity ba1, ba2;
1299
1300  printf(" MATRIX(amitype)\n");
1301  for (ba1 = AA; (long)ba1 <= (long)VV; ba1 = (amity)((long)ba1 + 1)) {
1302    for (ba2 = AA; (long)ba2 <= (long)VV; ba2 = (amity)((long)ba2 + 1)) {
1303      if (amtrx[(long)ba1 - (long)AA][(long)ba2 - (long)AA] > 0.1e-5)
1304        printf("%14.8f", amtrx[(long)ba1 - (long)AA][(long)ba2 - (long)AA]);
1305      else if (amtrx[(long)ba1 - (long)AA][(long)ba2 - (long)AA] <= 0.0)
1306        printf("%6c% .1E",
1307               ' ', amtrx[(long)ba1 - (long)AA][(long)ba2 - (long)AA]);
1308      else
1309        printf("%3c% .4E",
1310               ' ', amtrx[(long)ba1 - (long)AA][(long)ba2 - (long)AA]);
1311      if (ba2 == II || ba2 == VV || ba2 == CC || ba2 == PP_)
1312        putchar('\n');
1313    }
1314    putchar('\n');
1315  }
1316}  /* printamtrx */
1317
1318
1319Static Void printfreqdyhf()
1320{
1321  amity ba;
1322
1323  printf("\n Dayhoff's Amino Acid Frequencies:\n");
1324  printf("%13s", "A,R,N,D,C :");
1325  for (ba = AA; (long)ba <= (long)CC; ba = (amity)((long)ba + 1))
1326    printf("%7.3f", freqdyhf[(long)ba - (long)AA]);
1327  printf("\n%13s", "Q,E,G,H,I :");
1328  for (ba = QQ; (long)ba <= (long)II; ba = (amity)((long)ba + 1))
1329    printf("%7.3f", freqdyhf[(long)ba - (long)AA]);
1330  printf("\n%13s", "L,K,M,F,P :");
1331  for (ba = LL; (long)ba <= (long)PP_; ba = (amity)((long)ba + 1))
1332    printf("%7.3f", freqdyhf[(long)ba - (long)AA]);
1333  printf("\n%13s", "S,T,W,Y,V :");
1334  for (ba = SS; (long)ba <= (long)VV; ba = (amity)((long)ba + 1))
1335    printf("%7.3f", freqdyhf[(long)ba - (long)AA]);
1336  putchar('\n');
1337}  /* printfreqdyhf */
1338
1339
1340/*************************************/
1341/*** TRANSITION PROBABILITY MATRIX ***/
1342/*************************************/
1343/*** READ MATRIX ( EIGEN VALUE,VECTOR AND FREQUENCY ) ***/
1344Static Void readeigenv()
1345{
1346  amity ba1, ba2;
1347
1348  for (ba2 = AA; (long)ba2 <= (long)VV; ba2 = (amity)((long)ba2 + 1)) {
1349    for (ba1 = AA; (long)ba1 <= (long)VV; ba1 = (amity)((long)ba1 + 1)) {
1350      if (P_eoln(tpmfile)) {
1351        fscanf(tpmfile, "%*[^\n]");
1352        getc(tpmfile);
1353      }
1354      fscanf(tpmfile, "%lg", &ev[(long)ba1 - (long)AA][(long)ba2 - (long)AA]);
1355    }
1356  }
1357  fscanf(tpmfile, "%*[^\n]");
1358  getc(tpmfile);
1359  for (ba1 = AA; (long)ba1 <= (long)VV; ba1 = (amity)((long)ba1 + 1)) {
1360    if (P_eoln(tpmfile)) {
1361      fscanf(tpmfile, "%*[^\n]");
1362      getc(tpmfile);
1363    }
1364    fscanf(tpmfile, "%lg", &eval[(long)ba1 - (long)AA]);
1365  }
1366  fscanf(tpmfile, "%*[^\n]");
1367  getc(tpmfile);
1368  for (ba1 = AA; (long)ba1 <= (long)VV; ba1 = (amity)((long)ba1 + 1)) {
1369    if (P_eoln(tpmfile)) {
1370      fscanf(tpmfile, "%*[^\n]");
1371      getc(tpmfile);
1372    }
1373    fscanf(tpmfile, "%lg%*[^\n]", &freqdyhf[(long)ba1 - (long)AA]);
1374    getc(tpmfile);
1375  }
1376}  /* readeigenv */
1377
1378
1379/*** DATA MATRIX ( EIGEN VALUE,VECTOR AND FREQUENCY ) ***/
1380Static Void dataeigenv()
1381{
1382  ev[0][0] = -2.18920e-01;
1383  ev[(long)RR - (long)AA][0] = -4.01968e-02;
1384  ev[(long)NN - (long)AA][0] = -2.27758e-01;
1385  ev[(long)DD - (long)AA][0] = -4.04699e-01;
1386  ev[(long)CC - (long)AA][0] = -3.49083e-02;
1387  ev[(long)QQ - (long)AA][0] = -2.38091e-01;
1388  ev[(long)EE - (long)AA][0] = 5.24114e-01;
1389  ev[(long)GG - (long)AA][0] = -2.34117e-02;
1390  ev[(long)HH - (long)AA][0] = 9.97954e-02;
1391  ev[(long)II - (long)AA][0] = -8.45249e-03;
1392  ev[(long)LL - (long)AA][0] = 5.28066e-03;
1393  ev[(long)KK - (long)AA][0] = 2.05284e-02;
1394  ev[(long)MM - (long)AA][0] = -1.23853e-02;
1395  ev[(long)FF - (long)AA][0] = -9.50275e-03;
1396  ev[(long)PP_ - (long)AA][0] = -3.24017e-02;
1397  ev[(long)SS - (long)AA][0] = 5.89404e-01;
1398  ev[(long)TT - (long)AA][0] = -1.99416e-01;
1399  ev[(long)WW - (long)AA][0] = -1.52404e-02;
1400  ev[(long)YY - (long)AA][0] = -1.81586e-03;
1401  ev[(long)VV - (long)AA][0] = 4.98109e-02;
1402
1403  ev[0][(long)RR - (long)AA] = 4.00877e-02;
1404  ev[(long)RR - (long)AA][(long)RR - (long)AA] = 3.76303e-02;
1405  ev[(long)NN - (long)AA][(long)RR - (long)AA] = 7.25090e-01;
1406  ev[(long)DD - (long)AA][(long)RR - (long)AA] = -5.28042e-01;
1407  ev[(long)CC - (long)AA][(long)RR - (long)AA] = 1.46525e-02;
1408  ev[(long)QQ - (long)AA][(long)RR - (long)AA] = -8.66164e-02;
1409  ev[(long)EE - (long)AA][(long)RR - (long)AA] = 3.20299e-01;
1410  ev[(long)GG - (long)AA][(long)RR - (long)AA] = 7.74478e-03;
1411  ev[(long)HH - (long)AA][(long)RR - (long)AA] = -8.90059e-02;
1412  ev[(long)II - (long)AA][(long)RR - (long)AA] = -2.94900e-02;
1413  ev[(long)LL - (long)AA][(long)RR - (long)AA] = -3.50400e-03;
1414  ev[(long)KK - (long)AA][(long)RR - (long)AA] = -5.22374e-02;
1415  ev[(long)MM - (long)AA][(long)RR - (long)AA] = 1.94473e-02;
1416  ev[(long)FF - (long)AA][(long)RR - (long)AA] = 5.80870e-03;
1417  ev[(long)PP_ - (long)AA][(long)RR - (long)AA] = 1.62922e-02;
1418  ev[(long)SS - (long)AA][(long)RR - (long)AA] = -2.60397e-01;
1419  ev[(long)TT - (long)AA][(long)RR - (long)AA] = 4.22891e-02;
1420  ev[(long)WW - (long)AA][(long)RR - (long)AA] = 2.42879e-03;
1421  ev[(long)YY - (long)AA][(long)RR - (long)AA] = -1.46244e-02;
1422  ev[(long)VV - (long)AA][(long)RR - (long)AA] = -5.05405e-04;
1423
1424  ev[0][(long)NN - (long)AA] = -4.62992e-01;
1425  ev[(long)RR - (long)AA][(long)NN - (long)AA] = 2.14018e-02;
1426  ev[(long)NN - (long)AA][(long)NN - (long)AA] = 1.71750e-01;
1427  ev[(long)DD - (long)AA][(long)NN - (long)AA] = 1.02440e-01;
1428  ev[(long)CC - (long)AA][(long)NN - (long)AA] = 3.17009e-03;
1429  ev[(long)QQ - (long)AA][(long)NN - (long)AA] = 1.50618e-01;
1430  ev[(long)EE - (long)AA][(long)NN - (long)AA] = -1.07848e-01;
1431  ev[(long)GG - (long)AA][(long)NN - (long)AA] = 5.62329e-02;
1432  ev[(long)HH - (long)AA][(long)NN - (long)AA] = -1.09388e-01;
1433  ev[(long)II - (long)AA][(long)NN - (long)AA] = -6.52572e-01;
1434  ev[(long)LL - (long)AA][(long)NN - (long)AA] = 1.46263e-02;
1435  ev[(long)KK - (long)AA][(long)NN - (long)AA] = -5.22977e-02;
1436  ev[(long)MM - (long)AA][(long)NN - (long)AA] = 4.61043e-02;
1437  ev[(long)FF - (long)AA][(long)NN - (long)AA] = 4.16223e-02;
1438  ev[(long)PP_ - (long)AA][(long)NN - (long)AA] = 6.60296e-02;
1439  ev[(long)SS - (long)AA][(long)NN - (long)AA] = 4.86194e-02;
1440  ev[(long)TT - (long)AA][(long)NN - (long)AA] = 3.28837e-01;
1441  ev[(long)WW - (long)AA][(long)NN - (long)AA] = -4.61826e-03;
1442  ev[(long)YY - (long)AA][(long)NN - (long)AA] = -8.60905e-03;
1443  ev[(long)VV - (long)AA][(long)NN - (long)AA] = 3.84867e-01;
1444
1445  ev[0][(long)DD - (long)AA] = 2.47117e-02;
1446  ev[(long)RR - (long)AA][(long)DD - (long)AA] = -3.76030e-02;
1447  ev[(long)NN - (long)AA][(long)DD - (long)AA] = 5.56820e-01;
1448  ev[(long)DD - (long)AA][(long)DD - (long)AA] = 5.60598e-02;
1449  ev[(long)CC - (long)AA][(long)DD - (long)AA] = -2.51572e-02;
1450  ev[(long)QQ - (long)AA][(long)DD - (long)AA] = 3.40304e-01;
1451  ev[(long)EE - (long)AA][(long)DD - (long)AA] = -4.10919e-01;
1452  ev[(long)GG - (long)AA][(long)DD - (long)AA] = -7.16450e-02;
1453  ev[(long)HH - (long)AA][(long)DD - (long)AA] = -2.14248e-01;
1454  ev[(long)II - (long)AA][(long)DD - (long)AA] = 5.25788e-02;
1455  ev[(long)LL - (long)AA][(long)DD - (long)AA] = -1.22652e-02;
1456  ev[(long)KK - (long)AA][(long)DD - (long)AA] = -5.68214e-02;
1457  ev[(long)MM - (long)AA][(long)DD - (long)AA] = 1.75995e-02;
1458  ev[(long)FF - (long)AA][(long)DD - (long)AA] = -7.65321e-03;
1459  ev[(long)PP_ - (long)AA][(long)DD - (long)AA] = -5.79082e-02;
1460  ev[(long)SS - (long)AA][(long)DD - (long)AA] = 3.84128e-01;
1461  ev[(long)TT - (long)AA][(long)DD - (long)AA] = -4.36123e-01;
1462  ev[(long)WW - (long)AA][(long)DD - (long)AA] = -1.26346e-02;
1463  ev[(long)YY - (long)AA][(long)DD - (long)AA] = -3.19921e-03;
1464  ev[(long)VV - (long)AA][(long)DD - (long)AA] = 2.57284e-02;
1465
1466  ev[0][(long)CC - (long)AA] = -2.23607e-01;
1467  ev[(long)RR - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1468  ev[(long)NN - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1469  ev[(long)DD - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1470  ev[(long)CC - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1471  ev[(long)QQ - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1472  ev[(long)EE - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1473  ev[(long)GG - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1474  ev[(long)HH - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1475  ev[(long)II - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1476  ev[(long)LL - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1477  ev[(long)KK - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1478  ev[(long)MM - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1479  ev[(long)FF - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1480  ev[(long)PP_ - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1481  ev[(long)SS - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1482  ev[(long)TT - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1483  ev[(long)WW - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1484  ev[(long)YY - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1485  ev[(long)VV - (long)AA][(long)CC - (long)AA] = -2.23607e-01;
1486
1487  ev[0][(long)QQ - (long)AA] = 4.04873e-01;
1488  ev[(long)RR - (long)AA][(long)QQ - (long)AA] = 6.49103e-03;
1489  ev[(long)NN - (long)AA][(long)QQ - (long)AA] = -1.19891e-01;
1490  ev[(long)DD - (long)AA][(long)QQ - (long)AA] = -3.67766e-03;
1491  ev[(long)CC - (long)AA][(long)QQ - (long)AA] = -9.46397e-04;
1492  ev[(long)QQ - (long)AA][(long)QQ - (long)AA] = -1.89655e-01;
1493  ev[(long)EE - (long)AA][(long)QQ - (long)AA] = 5.91497e-02;
1494  ev[(long)GG - (long)AA][(long)QQ - (long)AA] = -8.40994e-02;
1495  ev[(long)HH - (long)AA][(long)QQ - (long)AA] = 8.85589e-02;
1496  ev[(long)II - (long)AA][(long)QQ - (long)AA] = -7.24798e-01;
1497  ev[(long)LL - (long)AA][(long)QQ - (long)AA] = 2.13078e-02;
1498  ev[(long)KK - (long)AA][(long)QQ - (long)AA] = 6.70215e-02;
1499  ev[(long)MM - (long)AA][(long)QQ - (long)AA] = 4.33730e-02;
1500  ev[(long)FF - (long)AA][(long)QQ - (long)AA] = 4.68020e-02;
1501  ev[(long)PP_ - (long)AA][(long)QQ - (long)AA] = -7.75734e-02;
1502  ev[(long)SS - (long)AA][(long)QQ - (long)AA] = -6.04936e-02;
1503  ev[(long)TT - (long)AA][(long)QQ - (long)AA] = -3.20011e-01;
1504  ev[(long)WW - (long)AA][(long)QQ - (long)AA] = 6.03313e-04;
1505  ev[(long)YY - (long)AA][(long)QQ - (long)AA] = -9.61127e-03;
1506  ev[(long)VV - (long)AA][(long)QQ - (long)AA] = 3.47473e-01;
1507
1508  ev[0][(long)EE - (long)AA] = -4.99671e-02;
1509  ev[(long)RR - (long)AA][(long)EE - (long)AA] = 4.11185e-02;
1510  ev[(long)NN - (long)AA][(long)EE - (long)AA] = 1.47778e-01;
1511  ev[(long)DD - (long)AA][(long)EE - (long)AA] = 1.64073e-01;
1512  ev[(long)CC - (long)AA][(long)EE - (long)AA] = -1.66608e-03;
1513  ev[(long)QQ - (long)AA][(long)EE - (long)AA] = -3.40600e-01;
1514  ev[(long)EE - (long)AA][(long)EE - (long)AA] = -2.75784e-02;
1515  ev[(long)GG - (long)AA][(long)EE - (long)AA] = -6.14738e-03;
1516  ev[(long)HH - (long)AA][(long)EE - (long)AA] = 7.45280e-02;
1517  ev[(long)II - (long)AA][(long)EE - (long)AA] = 6.19257e-02;
1518  ev[(long)LL - (long)AA][(long)EE - (long)AA] = 8.66441e-02;
1519  ev[(long)KK - (long)AA][(long)EE - (long)AA] = 3.88639e-02;
1520  ev[(long)MM - (long)AA][(long)EE - (long)AA] = -8.97628e-01;
1521  ev[(long)FF - (long)AA][(long)EE - (long)AA] = -3.67076e-03;
1522  ev[(long)PP_ - (long)AA][(long)EE - (long)AA] = 4.14761e-02;
1523  ev[(long)SS - (long)AA][(long)EE - (long)AA] = 1.90444e-03;
1524  ev[(long)TT - (long)AA][(long)EE - (long)AA] = -4.52208e-02;
1525  ev[(long)WW - (long)AA][(long)EE - (long)AA] = -8.08550e-03;
1526  ev[(long)YY - (long)AA][(long)EE - (long)AA] = -1.28114e-02;
1527  ev[(long)VV - (long)AA][(long)EE - (long)AA] = 4.57149e-02;
1528
1529  ev[0][(long)GG - (long)AA] = -6.76985e-02;
1530  ev[(long)RR - (long)AA][(long)GG - (long)AA] = 1.07693e-01;
1531  ev[(long)NN - (long)AA][(long)GG - (long)AA] = 1.83827e-01;
1532  ev[(long)DD - (long)AA][(long)GG - (long)AA] = 2.57799e-01;
1533  ev[(long)CC - (long)AA][(long)GG - (long)AA] = -8.80400e-04;
1534  ev[(long)QQ - (long)AA][(long)GG - (long)AA] = -5.07282e-01;
1535  ev[(long)EE - (long)AA][(long)GG - (long)AA] = -2.08116e-02;
1536  ev[(long)GG - (long)AA][(long)GG - (long)AA] = -9.91077e-03;
1537  ev[(long)HH - (long)AA][(long)GG - (long)AA] = 1.32778e-01;
1538  ev[(long)II - (long)AA][(long)GG - (long)AA] = 3.30619e-02;
1539  ev[(long)LL - (long)AA][(long)GG - (long)AA] = -5.54361e-02;
1540  ev[(long)KK - (long)AA][(long)GG - (long)AA] = -7.58113e-02;
1541  ev[(long)MM - (long)AA][(long)GG - (long)AA] = 7.67118e-01;
1542  ev[(long)FF - (long)AA][(long)GG - (long)AA] = -8.15430e-03;
1543  ev[(long)PP_ - (long)AA][(long)GG - (long)AA] = 5.85005e-02;
1544  ev[(long)SS - (long)AA][(long)GG - (long)AA] = 1.59381e-02;
1545  ev[(long)TT - (long)AA][(long)GG - (long)AA] = -6.67189e-02;
1546  ev[(long)WW - (long)AA][(long)GG - (long)AA] = -9.06354e-03;
1547  ev[(long)YY - (long)AA][(long)GG - (long)AA] = -6.80482e-03;
1548  ev[(long)VV - (long)AA][(long)GG - (long)AA] = -3.69299e-02;
1549
1550  ev[0][(long)HH - (long)AA] = 2.37414e-02;
1551  ev[(long)RR - (long)AA][(long)HH - (long)AA] = -1.31704e-02;
1552  ev[(long)NN - (long)AA][(long)HH - (long)AA] = 1.93358e-02;
1553  ev[(long)DD - (long)AA][(long)HH - (long)AA] = 2.78387e-02;
1554  ev[(long)CC - (long)AA][(long)HH - (long)AA] = 6.35183e-02;
1555  ev[(long)QQ - (long)AA][(long)HH - (long)AA] = 1.94232e-02;
1556  ev[(long)EE - (long)AA][(long)HH - (long)AA] = 2.72120e-02;
1557  ev[(long)GG - (long)AA][(long)HH - (long)AA] = 3.14653e-02;
1558  ev[(long)HH - (long)AA][(long)HH - (long)AA] = 8.67994e-03;
1559  ev[(long)II - (long)AA][(long)HH - (long)AA] = 3.64294e-03;
1560  ev[(long)LL - (long)AA][(long)HH - (long)AA] = -1.65127e-02;
1561  ev[(long)KK - (long)AA][(long)HH - (long)AA] = 1.45011e-02;
1562  ev[(long)MM - (long)AA][(long)HH - (long)AA] = 1.53344e-04;
1563  ev[(long)FF - (long)AA][(long)HH - (long)AA] = -7.14500e-02;
1564  ev[(long)PP_ - (long)AA][(long)HH - (long)AA] = 2.47183e-02;
1565  ev[(long)SS - (long)AA][(long)HH - (long)AA] = 1.83958e-02;
1566  ev[(long)TT - (long)AA][(long)HH - (long)AA] = 2.03333e-02;
1567  ev[(long)WW - (long)AA][(long)HH - (long)AA] = -9.90001e-01;
1568  ev[(long)YY - (long)AA][(long)HH - (long)AA] = -6.85327e-02;
1569  ev[(long)VV - (long)AA][(long)HH - (long)AA] = 1.15883e-02;
1570
1571  ev[0][(long)II - (long)AA] = 3.86884e-02;
1572  ev[(long)RR - (long)AA][(long)II - (long)AA] = 6.83612e-02;
1573  ev[(long)NN - (long)AA][(long)II - (long)AA] = 5.94583e-02;
1574  ev[(long)DD - (long)AA][(long)II - (long)AA] = 7.89097e-02;
1575  ev[(long)CC - (long)AA][(long)II - (long)AA] = -9.49564e-01;
1576  ev[(long)QQ - (long)AA][(long)II - (long)AA] = 7.75597e-02;
1577  ev[(long)EE - (long)AA][(long)II - (long)AA] = 7.93973e-02;
1578  ev[(long)GG - (long)AA][(long)II - (long)AA] = 5.96830e-02;
1579  ev[(long)HH - (long)AA][(long)II - (long)AA] = 5.15293e-02;
1580  ev[(long)II - (long)AA][(long)II - (long)AA] = 7.67383e-03;
1581  ev[(long)LL - (long)AA][(long)II - (long)AA] = 2.21331e-02;
1582  ev[(long)KK - (long)AA][(long)II - (long)AA] = 8.30887e-02;
1583  ev[(long)MM - (long)AA][(long)II - (long)AA] = 3.95955e-02;
1584  ev[(long)FF - (long)AA][(long)II - (long)AA] = -1.11612e-01;
1585  ev[(long)PP_ - (long)AA][(long)II - (long)AA] = 5.03976e-02;
1586  ev[(long)SS - (long)AA][(long)II - (long)AA] = 1.59569e-02;
1587  ev[(long)TT - (long)AA][(long)II - (long)AA] = 3.72414e-02;
1588  ev[(long)WW - (long)AA][(long)II - (long)AA] = -5.52876e-02;
1589  ev[(long)YY - (long)AA][(long)II - (long)AA] = -1.86929e-01;
1590  ev[(long)VV - (long)AA][(long)II - (long)AA] = 1.41872e-02;
1591
1592  ev[0][(long)LL - (long)AA] = 5.80295e-02;
1593  ev[(long)RR - (long)AA][(long)LL - (long)AA] = 1.02103e-01;
1594  ev[(long)NN - (long)AA][(long)LL - (long)AA] = 7.00404e-02;
1595  ev[(long)DD - (long)AA][(long)LL - (long)AA] = 9.76903e-02;
1596  ev[(long)CC - (long)AA][(long)LL - (long)AA] = 2.47574e-01;
1597  ev[(long)QQ - (long)AA][(long)LL - (long)AA] = 7.38092e-02;
1598  ev[(long)EE - (long)AA][(long)LL - (long)AA] = 9.10580e-02;
1599  ev[(long)GG - (long)AA][(long)LL - (long)AA] = 1.00958e-01;
1600  ev[(long)HH - (long)AA][(long)LL - (long)AA] = 3.49872e-02;
1601  ev[(long)II - (long)AA][(long)LL - (long)AA] = -1.13481e-01;
1602  ev[(long)LL - (long)AA][(long)LL - (long)AA] = -2.12223e-01;
1603  ev[(long)KK - (long)AA][(long)LL - (long)AA] = 9.20232e-02;
1604  ev[(long)MM - (long)AA][(long)LL - (long)AA] = -1.06117e-01;
1605  ev[(long)FF - (long)AA][(long)LL - (long)AA] = -5.51278e-01;
1606  ev[(long)PP_ - (long)AA][(long)LL - (long)AA] = 8.17221e-02;
1607  ev[(long)SS - (long)AA][(long)LL - (long)AA] = 7.48437e-02;
1608  ev[(long)TT - (long)AA][(long)LL - (long)AA] = 4.59234e-02;
1609  ev[(long)WW - (long)AA][(long)LL - (long)AA] = 4.34903e-01;
1610  ev[(long)YY - (long)AA][(long)LL - (long)AA] = -5.43823e-01;
1611  ev[(long)VV - (long)AA][(long)LL - (long)AA] = -6.69564e-02;
1612
1613  ev[0][(long)KK - (long)AA] = 1.68796e-01;
1614  ev[(long)RR - (long)AA][(long)KK - (long)AA] = 6.98733e-01;
1615  ev[(long)NN - (long)AA][(long)KK - (long)AA] = -3.64477e-02;
1616  ev[(long)DD - (long)AA][(long)KK - (long)AA] = 4.79840e-02;
1617  ev[(long)CC - (long)AA][(long)KK - (long)AA] = -2.83352e-02;
1618  ev[(long)QQ - (long)AA][(long)KK - (long)AA] = 2.07546e-03;
1619  ev[(long)EE - (long)AA][(long)KK - (long)AA] = 6.44875e-02;
1620  ev[(long)GG - (long)AA][(long)KK - (long)AA] = -1.24957e-01;
1621  ev[(long)HH - (long)AA][(long)KK - (long)AA] = -2.31511e-01;
1622  ev[(long)II - (long)AA][(long)KK - (long)AA] = -1.08720e-01;
1623  ev[(long)LL - (long)AA][(long)KK - (long)AA] = 6.27788e-02;
1624  ev[(long)KK - (long)AA][(long)KK - (long)AA] = -4.17790e-01;
1625  ev[(long)MM - (long)AA][(long)KK - (long)AA] = -1.98195e-01;
1626  ev[(long)FF - (long)AA][(long)KK - (long)AA] = -1.01464e-02;
1627  ev[(long)PP_ - (long)AA][(long)KK - (long)AA] = -2.34527e-01;
1628  ev[(long)SS - (long)AA][(long)KK - (long)AA] = 1.80480e-01;
1629  ev[(long)TT - (long)AA][(long)KK - (long)AA] = 2.62775e-01;
1630  ev[(long)WW - (long)AA][(long)KK - (long)AA] = -7.67023e-02;
1631  ev[(long)YY - (long)AA][(long)KK - (long)AA] = 8.53465e-03;
1632  ev[(long)VV - (long)AA][(long)KK - (long)AA] = -1.14869e-01;
1633
1634  ev[0][(long)MM - (long)AA] = 1.14030e-02;
1635  ev[(long)RR - (long)AA][(long)MM - (long)AA] = 1.03453e-01;
1636  ev[(long)NN - (long)AA][(long)MM - (long)AA] = 1.03058e-01;
1637  ev[(long)DD - (long)AA][(long)MM - (long)AA] = 1.25170e-01;
1638  ev[(long)CC - (long)AA][(long)MM - (long)AA] = -7.83690e-02;
1639  ev[(long)QQ - (long)AA][(long)MM - (long)AA] = 8.15278e-02;
1640  ev[(long)EE - (long)AA][(long)MM - (long)AA] = 1.09160e-01;
1641  ev[(long)GG - (long)AA][(long)MM - (long)AA] = 9.53752e-02;
1642  ev[(long)HH - (long)AA][(long)MM - (long)AA] = 1.40899e-01;
1643  ev[(long)II - (long)AA][(long)MM - (long)AA] = -2.81849e-01;
1644  ev[(long)LL - (long)AA][(long)MM - (long)AA] = -4.92274e-01;
1645  ev[(long)KK - (long)AA][(long)MM - (long)AA] = 9.83942e-02;
1646  ev[(long)MM - (long)AA][(long)MM - (long)AA] = -3.14326e-01;
1647  ev[(long)FF - (long)AA][(long)MM - (long)AA] = 2.70954e-01;
1648  ev[(long)PP_ - (long)AA][(long)MM - (long)AA] = 4.05413e-02;
1649  ev[(long)SS - (long)AA][(long)MM - (long)AA] = 5.49397e-02;
1650  ev[(long)TT - (long)AA][(long)MM - (long)AA] = -2.30001e-04;
1651  ev[(long)WW - (long)AA][(long)MM - (long)AA] = -6.60170e-02;
1652  ev[(long)YY - (long)AA][(long)MM - (long)AA] = 5.64557e-01;
1653  ev[(long)VV - (long)AA][(long)MM - (long)AA] = -2.78943e-01;
1654
1655  ev[0][(long)FF - (long)AA] = 1.68903e-01;
1656  ev[(long)RR - (long)AA][(long)FF - (long)AA] = -4.16455e-01;
1657  ev[(long)NN - (long)AA][(long)FF - (long)AA] = 1.75235e-01;
1658  ev[(long)DD - (long)AA][(long)FF - (long)AA] = -1.57281e-01;
1659  ev[(long)CC - (long)AA][(long)FF - (long)AA] = -3.02415e-02;
1660  ev[(long)QQ - (long)AA][(long)FF - (long)AA] = -1.99945e-01;
1661  ev[(long)EE - (long)AA][(long)FF - (long)AA] = -2.80839e-01;
1662  ev[(long)GG - (long)AA][(long)FF - (long)AA] = -1.65168e-01;
1663  ev[(long)HH - (long)AA][(long)FF - (long)AA] = 3.84179e-01;
1664  ev[(long)II - (long)AA][(long)FF - (long)AA] = -1.72794e-01;
1665  ev[(long)LL - (long)AA][(long)FF - (long)AA] = 4.08470e-02;
1666  ev[(long)KK - (long)AA][(long)FF - (long)AA] = 1.09632e-01;
1667  ev[(long)MM - (long)AA][(long)FF - (long)AA] = 3.10366e-02;
1668  ev[(long)FF - (long)AA][(long)FF - (long)AA] = 1.78854e-02;
1669  ev[(long)PP_ - (long)AA][(long)FF - (long)AA] = -2.43651e-01;
1670  ev[(long)SS - (long)AA][(long)FF - (long)AA] = 2.58272e-01;
1671  ev[(long)TT - (long)AA][(long)FF - (long)AA] = 4.85439e-01;
1672  ev[(long)WW - (long)AA][(long)FF - (long)AA] = 1.87008e-02;
1673  ev[(long)YY - (long)AA][(long)FF - (long)AA] = -8.19317e-02;
1674  ev[(long)VV - (long)AA][(long)FF - (long)AA] = -1.85319e-01;
1675
1676  ev[0][(long)PP_ - (long)AA] = 1.89407e-01;
1677  ev[(long)RR - (long)AA][(long)PP_ - (long)AA] = -5.67638e-01;
1678  ev[(long)NN - (long)AA][(long)PP_ - (long)AA] = -2.92067e-02;
1679  ev[(long)DD - (long)AA][(long)PP_ - (long)AA] = 4.63283e-02;
1680  ev[(long)CC - (long)AA][(long)PP_ - (long)AA] = -7.50547e-02;
1681  ev[(long)QQ - (long)AA][(long)PP_ - (long)AA] = -1.77709e-01;
1682  ev[(long)EE - (long)AA][(long)PP_ - (long)AA] = 2.11012e-02;
1683  ev[(long)GG - (long)AA][(long)PP_ - (long)AA] = 4.57586e-01;
1684  ev[(long)HH - (long)AA][(long)PP_ - (long)AA] = -2.73917e-01;
1685  ev[(long)II - (long)AA][(long)PP_ - (long)AA] = 3.02761e-02;
1686  ev[(long)LL - (long)AA][(long)PP_ - (long)AA] = -8.55620e-02;
1687  ev[(long)KK - (long)AA][(long)PP_ - (long)AA] = -4.68228e-01;
1688  ev[(long)MM - (long)AA][(long)PP_ - (long)AA] = -1.49034e-01;
1689  ev[(long)FF - (long)AA][(long)PP_ - (long)AA] = 4.32709e-02;
1690  ev[(long)PP_ - (long)AA][(long)PP_ - (long)AA] = 1.13488e-01;
1691  ev[(long)SS - (long)AA][(long)PP_ - (long)AA] = 1.12224e-01;
1692  ev[(long)TT - (long)AA][(long)PP_ - (long)AA] = 8.54433e-02;
1693  ev[(long)WW - (long)AA][(long)PP_ - (long)AA] = 1.49702e-01;
1694  ev[(long)YY - (long)AA][(long)PP_ - (long)AA] = 1.44889e-02;
1695  ev[(long)VV - (long)AA][(long)PP_ - (long)AA] = 9.94179e-02;
1696
1697  ev[0][(long)SS - (long)AA] = -4.51742e-02;
1698  ev[(long)RR - (long)AA][(long)SS - (long)AA] = 2.72843e-01;
1699  ev[(long)NN - (long)AA][(long)SS - (long)AA] = -6.34739e-02;
1700  ev[(long)DD - (long)AA][(long)SS - (long)AA] = -2.99613e-01;
1701  ev[(long)CC - (long)AA][(long)SS - (long)AA] = -1.26065e-02;
1702  ev[(long)QQ - (long)AA][(long)SS - (long)AA] = 6.62398e-02;
1703  ev[(long)EE - (long)AA][(long)SS - (long)AA] = -2.93927e-01;
1704  ev[(long)GG - (long)AA][(long)SS - (long)AA] = 2.07444e-01;
1705  ev[(long)HH - (long)AA][(long)SS - (long)AA] = 7.61271e-01;
1706  ev[(long)II - (long)AA][(long)SS - (long)AA] = 1.22763e-01;
1707  ev[(long)LL - (long)AA][(long)SS - (long)AA] = -9.12377e-02;
1708  ev[(long)KK - (long)AA][(long)SS - (long)AA] = -1.67370e-01;
1709  ev[(long)MM - (long)AA][(long)SS - (long)AA] = -7.69871e-02;
1710  ev[(long)FF - (long)AA][(long)SS - (long)AA] = 3.25168e-02;
1711  ev[(long)PP_ - (long)AA][(long)SS - (long)AA] = -8.69178e-02;
1712  ev[(long)SS - (long)AA][(long)SS - (long)AA] = -4.91352e-02;
1713  ev[(long)TT - (long)AA][(long)SS - (long)AA] = -9.28875e-02;
1714  ev[(long)WW - (long)AA][(long)SS - (long)AA] = -3.40158e-02;
1715  ev[(long)YY - (long)AA][(long)SS - (long)AA] = -9.70949e-02;
1716  ev[(long)VV - (long)AA][(long)SS - (long)AA] = 1.69233e-01;
1717
1718  ev[0][(long)TT - (long)AA] = -6.57152e-02;
1719  ev[(long)RR - (long)AA][(long)TT - (long)AA] = -2.96127e-01;
1720  ev[(long)NN - (long)AA][(long)TT - (long)AA] = 1.40073e-01;
1721  ev[(long)DD - (long)AA][(long)TT - (long)AA] = 3.67827e-01;
1722  ev[(long)CC - (long)AA][(long)TT - (long)AA] = 3.88665e-02;
1723  ev[(long)QQ - (long)AA][(long)TT - (long)AA] = 3.54579e-01;
1724  ev[(long)EE - (long)AA][(long)TT - (long)AA] = 3.95304e-01;
1725  ev[(long)GG - (long)AA][(long)TT - (long)AA] = -1.30872e-01;
1726  ev[(long)HH - (long)AA][(long)TT - (long)AA] = 5.22294e-01;
1727  ev[(long)II - (long)AA][(long)TT - (long)AA] = -9.20181e-02;
1728  ev[(long)LL - (long)AA][(long)TT - (long)AA] = 1.35098e-01;
1729  ev[(long)KK - (long)AA][(long)TT - (long)AA] = -2.61406e-01;
1730  ev[(long)MM - (long)AA][(long)TT - (long)AA] = -5.72654e-02;
1731  ev[(long)FF - (long)AA][(long)TT - (long)AA] = -1.06748e-01;
1732  ev[(long)PP_ - (long)AA][(long)TT - (long)AA] = -1.86932e-01;
1733  ev[(long)SS - (long)AA][(long)TT - (long)AA] = -8.60432e-02;
1734  ev[(long)TT - (long)AA][(long)TT - (long)AA] = -1.17435e-01;
1735  ev[(long)WW - (long)AA][(long)TT - (long)AA] = 4.21809e-02;
1736  ev[(long)YY - (long)AA][(long)TT - (long)AA] = 3.66170e-02;
1737  ev[(long)VV - (long)AA][(long)TT - (long)AA] = -1.03287e-01;
1738
1739  ev[0][(long)WW - (long)AA] = 8.12248e-02;
1740  ev[(long)RR - (long)AA][(long)WW - (long)AA] = -5.52327e-02;
1741  ev[(long)NN - (long)AA][(long)WW - (long)AA] = -5.60451e-02;
1742  ev[(long)DD - (long)AA][(long)WW - (long)AA] = -1.10967e-01;
1743  ev[(long)CC - (long)AA][(long)WW - (long)AA] = -4.36007e-02;
1744  ev[(long)QQ - (long)AA][(long)WW - (long)AA] = 7.49285e-02;
1745  ev[(long)EE - (long)AA][(long)WW - (long)AA] = -6.08685e-02;
1746  ev[(long)GG - (long)AA][(long)WW - (long)AA] = -3.69332e-01;
1747  ev[(long)HH - (long)AA][(long)WW - (long)AA] = 1.94411e-01;
1748  ev[(long)II - (long)AA][(long)WW - (long)AA] = 3.55141e-02;
1749  ev[(long)LL - (long)AA][(long)WW - (long)AA] = -8.63681e-02;
1750  ev[(long)KK - (long)AA][(long)WW - (long)AA] = -2.09830e-01;
1751  ev[(long)MM - (long)AA][(long)WW - (long)AA] = -9.78879e-02;
1752  ev[(long)FF - (long)AA][(long)WW - (long)AA] = -7.98604e-02;
1753  ev[(long)PP_ - (long)AA][(long)WW - (long)AA] = 8.35838e-01;
1754  ev[(long)SS - (long)AA][(long)WW - (long)AA] = 4.95931e-02;
1755  ev[(long)TT - (long)AA][(long)WW - (long)AA] = 7.86337e-02;
1756  ev[(long)WW - (long)AA][(long)WW - (long)AA] = 1.01576e-02;
1757  ev[(long)YY - (long)AA][(long)WW - (long)AA] = 8.79878e-02;
1758  ev[(long)VV - (long)AA][(long)WW - (long)AA] = 7.51898e-02;
1759
1760  ev[0][(long)YY - (long)AA] = -3.45102e-02;
1761  ev[(long)RR - (long)AA][(long)YY - (long)AA] = 6.45521e-02;
1762  ev[(long)NN - (long)AA][(long)YY - (long)AA] = -2.22780e-02;
1763  ev[(long)DD - (long)AA][(long)YY - (long)AA] = -3.19792e-02;
1764  ev[(long)CC - (long)AA][(long)YY - (long)AA] = 5.97900e-02;
1765  ev[(long)QQ - (long)AA][(long)YY - (long)AA] = 4.40846e-02;
1766  ev[(long)EE - (long)AA][(long)YY - (long)AA] = -2.90400e-02;
1767  ev[(long)GG - (long)AA][(long)YY - (long)AA] = 1.32178e-01;
1768  ev[(long)HH - (long)AA][(long)YY - (long)AA] = 4.83456e-02;
1769  ev[(long)II - (long)AA][(long)YY - (long)AA] = -3.33680e-01;
1770  ev[(long)LL - (long)AA][(long)YY - (long)AA] = 1.96787e-01;
1771  ev[(long)KK - (long)AA][(long)YY - (long)AA] = -1.27120e-02;
1772  ev[(long)MM - (long)AA][(long)YY - (long)AA] = -1.01360e-02;
1773  ev[(long)FF - (long)AA][(long)YY - (long)AA] = 5.00463e-01;
1774  ev[(long)PP_ - (long)AA][(long)YY - (long)AA] = 2.65762e-01;
1775  ev[(long)SS - (long)AA][(long)YY - (long)AA] = 8.18628e-03;
1776  ev[(long)TT - (long)AA][(long)YY - (long)AA] = -1.15726e-01;
1777  ev[(long)WW - (long)AA][(long)YY - (long)AA] = -3.46187e-02;
1778  ev[(long)YY - (long)AA][(long)YY - (long)AA] = -5.64673e-01;
1779  ev[(long)VV - (long)AA][(long)YY - (long)AA] = -4.02511e-01;
1780
1781  ev[0][(long)VV - (long)AA] = -2.37427e-02;
1782  ev[(long)RR - (long)AA][(long)VV - (long)AA] = 6.64165e-02;
1783  ev[(long)NN - (long)AA][(long)VV - (long)AA] = -3.07931e-02;
1784  ev[(long)DD - (long)AA][(long)VV - (long)AA] = -1.35114e-01;
1785  ev[(long)CC - (long)AA][(long)VV - (long)AA] = -1.26575e-02;
1786  ev[(long)QQ - (long)AA][(long)VV - (long)AA] = -4.34887e-02;
1787  ev[(long)EE - (long)AA][(long)VV - (long)AA] = -1.42949e-01;
1788  ev[(long)GG - (long)AA][(long)VV - (long)AA] = 1.85463e-01;
1789  ev[(long)HH - (long)AA][(long)VV - (long)AA] = 4.82054e-02;
1790  ev[(long)II - (long)AA][(long)VV - (long)AA] = -2.88402e-01;
1791  ev[(long)LL - (long)AA][(long)VV - (long)AA] = 2.93123e-01;
1792  ev[(long)KK - (long)AA][(long)VV - (long)AA] = 1.32034e-02;
1793  ev[(long)MM - (long)AA][(long)VV - (long)AA] = 6.77690e-02;
1794  ev[(long)FF - (long)AA][(long)VV - (long)AA] = -5.43584e-01;
1795  ev[(long)PP_ - (long)AA][(long)VV - (long)AA] = 1.46520e-01;
1796  ev[(long)SS - (long)AA][(long)VV - (long)AA] = 3.28990e-03;
1797  ev[(long)TT - (long)AA][(long)VV - (long)AA] = -7.67072e-02;
1798  ev[(long)WW - (long)AA][(long)VV - (long)AA] = -2.03106e-02;
1799  ev[(long)YY - (long)AA][(long)VV - (long)AA] = 5.89467e-01;
1800  ev[(long)VV - (long)AA][(long)VV - (long)AA] = -2.68367e-01;
1801
1802
1803  eval[0] = -2.03036e-02;
1804  eval[(long)RR - (long)AA] = -2.33761e-02;
1805  eval[(long)NN - (long)AA] = -1.71812e-02;
1806  eval[(long)DD - (long)AA] = -1.82705e-02;
1807  eval[(long)CC - (long)AA] = -1.55431e-15;
1808  eval[(long)QQ - (long)AA] = -1.60581e-02;
1809  eval[(long)EE - (long)AA] = -1.34008e-02;
1810  eval[(long)GG - (long)AA] = -1.38363e-02;
1811  eval[(long)HH - (long)AA] = -2.36636e-03;
1812  eval[(long)II - (long)AA] = -2.68394e-03;
1813  eval[(long)LL - (long)AA] = -2.91392e-03;
1814  eval[(long)KK - (long)AA] = -1.13524e-02;
1815  eval[(long)MM - (long)AA] = -4.34547e-03;
1816  eval[(long)FF - (long)AA] = -1.06999e-02;
1817  eval[(long)PP_ - (long)AA] = -5.48466e-03;
1818  eval[(long)SS - (long)AA] = -8.98371e-03;
1819  eval[(long)TT - (long)AA] = -7.23244e-03;
1820  eval[(long)WW - (long)AA] = -7.37540e-03;
1821  eval[(long)YY - (long)AA] = -7.91291e-03;
1822  eval[(long)VV - (long)AA] = -8.24441e-03;
1823
1824
1825  freqdyhf[0] = 0.8712669e-01;
1826  freqdyhf[(long)RR - (long)AA] = 0.4090396e-01;
1827  freqdyhf[(long)NN - (long)AA] = 0.4043229e-01;
1828  freqdyhf[(long)DD - (long)AA] = 0.4687196e-01;
1829  freqdyhf[(long)CC - (long)AA] = 0.3347348e-01;
1830  freqdyhf[(long)QQ - (long)AA] = 0.3825543e-01;
1831  freqdyhf[(long)EE - (long)AA] = 0.4953036e-01;
1832  freqdyhf[(long)GG - (long)AA] = 0.8861207e-01;
1833  freqdyhf[(long)HH - (long)AA] = 0.3361838e-01;
1834  freqdyhf[(long)II - (long)AA] = 0.3688570e-01;
1835  freqdyhf[(long)LL - (long)AA] = 0.8535736e-01;
1836  freqdyhf[(long)KK - (long)AA] = 0.8048168e-01;
1837  freqdyhf[(long)MM - (long)AA] = 0.1475275e-01;
1838  freqdyhf[(long)FF - (long)AA] = 0.3977166e-01;
1839  freqdyhf[(long)PP_ - (long)AA] = 0.5067984e-01;
1840  freqdyhf[(long)SS - (long)AA] = 0.6957710e-01;
1841  freqdyhf[(long)TT - (long)AA] = 0.5854172e-01;
1842  freqdyhf[(long)WW - (long)AA] = 0.1049367e-01;
1843  freqdyhf[(long)YY - (long)AA] = 0.2991627e-01;
1844  freqdyhf[(long)VV - (long)AA] = 0.6471762e-01;
1845}  /* dataeigenv */
1846
1847
1848/*** TRANSITION PROBABILITY MATRIX ***/
1849Static Void tprobmtrx(arc, tpr)
1850double arc;
1851aryamity *tpr;
1852{
1853  /* transition probability matrix */
1854  double sum;
1855  amity iba, jba, kba;
1856  aryamity exparc;
1857
1858  /* negative : BOOLEAN; */
1859  /* negative := FALSE; */
1860  for (kba = AA; (long)kba <= (long)VV; kba = (amity)((long)kba + 1))
1861    exparc[(long)kba - (long)AA] = exp(arc * eval[(long)kba - (long)AA]);
1862  for (iba = AA; (long)iba <= (long)VV; iba = (amity)((long)iba + 1)) {
1863    for (jba = AA; (long)jba <= (long)VV; jba = (amity)((long)jba + 1)) {
1864      sum = coefp[(long)iba - (long)AA][(long)jba - (long)AA][0] * exparc[0] +
1865          coefp[(long)iba - (long)AA][(long)jba - (long)AA][(long)RR - (long)
1866                  AA] * exparc[(long)RR - (long)AA] + coefp[(long)iba - (long)
1867                  AA][(long)jba - (long)AA][(long)NN - (long)AA] *
1868            exparc[(long)NN - (long)AA] + coefp[(long)iba - (long)AA][(long)
1869                  jba - (long)AA][(long)DD - (long)AA] * exparc[(long)DD -
1870                (long)AA] + coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1871            [(long)CC - (long)AA] * exparc[(long)CC - (long)AA] + coefp[(long)
1872                  iba - (long)AA][(long)jba - (long)AA][(long)QQ - (long)AA] *
1873            exparc[(long)QQ - (long)AA] + coefp[(long)iba - (long)AA][(long)
1874                  jba - (long)AA][(long)EE - (long)AA] * exparc[(long)EE -
1875                (long)AA] + coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1876            [(long)GG - (long)AA] * exparc[(long)GG - (long)AA] + coefp[(long)
1877                  iba - (long)AA][(long)jba - (long)AA][(long)HH - (long)AA] *
1878            exparc[(long)HH - (long)AA] + coefp[(long)iba - (long)AA][(long)
1879                  jba - (long)AA][(long)II - (long)AA] * exparc[(long)II -
1880                (long)AA] + coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1881            [(long)LL - (long)AA] * exparc[(long)LL - (long)AA] + coefp[(long)
1882                  iba - (long)AA][(long)jba - (long)AA][(long)KK - (long)AA] *
1883            exparc[(long)KK - (long)AA] + coefp[(long)iba - (long)AA][(long)
1884                  jba - (long)AA][(long)MM - (long)AA] * exparc[(long)MM -
1885                (long)AA] + coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1886            [(long)FF - (long)AA] * exparc[(long)FF - (long)AA] + coefp[(long)
1887                  iba - (long)AA][(long)jba - (long)AA]
1888            [(long)PP_ - (long)AA] * exparc[(long)PP_ - (long)AA] +
1889          coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1890            [(long)SS - (long)AA] * exparc[(long)SS - (long)AA] +
1891          coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1892            [(long)TT - (long)AA] * exparc[(long)TT - (long)AA] +
1893          coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1894            [(long)WW - (long)AA] * exparc[(long)WW - (long)AA] +
1895          coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1896            [(long)YY - (long)AA] * exparc[(long)YY - (long)AA] +
1897          coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1898            [(long)VV - (long)AA] * exparc[(long)VV - (long)AA];
1899/* p2c: protml.pas, line 1462: Note:
1900 * Line breaker spent 0.0+1.00 seconds, 5000 tries on line 1898 [251] */
1901      if (sum < minreal)   /* negative := TRUE; */
1902        tpr[(long)iba - (long)AA][(long)jba - (long)AA] = 0.0;
1903            /* attention */
1904      else
1905        tpr[(long)iba - (long)AA][(long)jba - (long)AA] = sum;
1906    }
1907  }
1908  /* IF negative THEN
1909     IF debugoptn THEN
1910     BEGIN
1911        write(output,' Trans. prob. 1 is negative !');
1912        writeln(output,'  arc =',arc:10:5);
1913     END; */
1914}  /* tprobmtrx */
1915
1916
1917/*** TRANSITION PROBABILITY AND DIFFRENCE MATRIX ***/
1918Static Void tdiffmtrx(arc, tpr, td1, td2)
1919double arc;
1920aryamity *tpr, *td1, *td2;
1921{
1922  /* transition probability matrix */
1923  double sum, sumd1, sumd2, aaa, rrr, nnn, ddd, ccc, qqq, eee, ggg, hhh, iii,
1924         lll, kkk, mmm, fff, ppp, sss, ttt, www, yyy, vvv;
1925  amity iba, jba, kba;
1926  aryamity exparc;
1927
1928  for (kba = AA; (long)kba <= (long)VV; kba = (amity)((long)kba + 1))
1929    exparc[(long)kba - (long)AA] = exp(arc * eval[(long)kba - (long)AA]);
1930  for (iba = AA; (long)iba <= (long)VV; iba = (amity)((long)iba + 1)) {
1931    for (jba = AA; (long)jba <= (long)VV; jba = (amity)((long)jba + 1)) {
1932      aaa = coefp[(long)iba - (long)AA][(long)jba - (long)AA][0] * exparc[0];
1933      rrr = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1934            [(long)RR - (long)AA] * exparc[(long)RR - (long)AA];
1935      nnn = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1936            [(long)NN - (long)AA] * exparc[(long)NN - (long)AA];
1937      ddd = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1938            [(long)DD - (long)AA] * exparc[(long)DD - (long)AA];
1939      ccc = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1940            [(long)CC - (long)AA] * exparc[(long)CC - (long)AA];
1941      qqq = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1942            [(long)QQ - (long)AA] * exparc[(long)QQ - (long)AA];
1943      eee = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1944            [(long)EE - (long)AA] * exparc[(long)EE - (long)AA];
1945      ggg = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1946            [(long)GG - (long)AA] * exparc[(long)GG - (long)AA];
1947      hhh = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1948            [(long)HH - (long)AA] * exparc[(long)HH - (long)AA];
1949      iii = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1950            [(long)II - (long)AA] * exparc[(long)II - (long)AA];
1951      lll = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1952            [(long)LL - (long)AA] * exparc[(long)LL - (long)AA];
1953      kkk = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1954            [(long)KK - (long)AA] * exparc[(long)KK - (long)AA];
1955      mmm = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1956            [(long)MM - (long)AA] * exparc[(long)MM - (long)AA];
1957      fff = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1958            [(long)FF - (long)AA] * exparc[(long)FF - (long)AA];
1959      ppp = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1960            [(long)PP_ - (long)AA] * exparc[(long)PP_ - (long)AA];
1961      sss = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1962            [(long)SS - (long)AA] * exparc[(long)SS - (long)AA];
1963      ttt = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1964            [(long)TT - (long)AA] * exparc[(long)TT - (long)AA];
1965      www = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1966            [(long)WW - (long)AA] * exparc[(long)WW - (long)AA];
1967      yyy = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1968            [(long)YY - (long)AA] * exparc[(long)YY - (long)AA];
1969      vvv = coefp[(long)iba - (long)AA][(long)jba - (long)AA]
1970            [(long)VV - (long)AA] * exparc[(long)VV - (long)AA];
1971
1972      sum = aaa + rrr + nnn + ddd + ccc + qqq + eee + ggg + hhh + iii + lll +
1973            kkk + mmm + fff + ppp + sss + ttt + www + yyy + vvv;
1974      sumd1 = aaa * eval[0] + rrr * eval[(long)RR - (long)AA] +
1975          nnn * eval[(long)NN - (long)AA] + ddd * eval[(long)DD - (long)AA] +
1976          ccc * eval[(long)CC - (long)AA] + qqq * eval[(long)QQ - (long)AA] +
1977          eee * eval[(long)EE - (long)AA] + ggg * eval[(long)GG - (long)AA] +
1978          hhh * eval[(long)HH - (long)AA] + iii * eval[(long)II - (long)AA] +
1979          lll * eval[(long)LL - (long)AA] + kkk * eval[(long)KK - (long)AA] +
1980          mmm * eval[(long)MM - (long)AA] + fff * eval[(long)FF - (long)AA] +
1981          ppp * eval[(long)PP_ - (long)AA] + sss * eval[(long)SS - (long)AA] +
1982          ttt * eval[(long)TT - (long)AA] + www * eval[(long)WW - (long)AA] +
1983          yyy * eval[(long)YY - (long)AA] + vvv * eval[(long)VV - (long)AA];
1984/* p2c: protml.pas, line 1542: Note:
1985 * Line breaker spent 0.0+1.00 seconds, 5000 tries on line 1983 [251] */
1986      sumd2 = aaa * eval2[0] + rrr * eval2[(long)RR - (long)AA] +
1987          nnn * eval2[(long)NN - (long)AA] +
1988          ddd * eval2[(long)DD - (long)AA] +
1989          ccc * eval2[(long)CC - (long)AA] +
1990          qqq * eval2[(long)QQ - (long)AA] +
1991          eee * eval2[(long)EE - (long)AA] +
1992          ggg * eval2[(long)GG - (long)AA] +
1993          hhh * eval2[(long)HH - (long)AA] +
1994          iii * eval2[(long)II - (long)AA] +
1995          lll * eval2[(long)LL - (long)AA] +
1996          kkk * eval2[(long)KK - (long)AA] +
1997          mmm * eval2[(long)MM - (long)AA] +
1998          fff * eval2[(long)FF - (long)AA] +
1999          ppp * eval2[(long)PP_ - (long)AA] +
2000          sss * eval2[(long)SS - (long)AA] +
2001          ttt * eval2[(long)TT - (long)AA] +
2002          www * eval2[(long)WW - (long)AA] +
2003          yyy * eval2[(long)YY - (long)AA] + vvv * eval2[(long)VV - (long)AA];
2004/* p2c: protml.pas, line 1542:
2005 * Note: Line breaker spent 0.0 seconds, 5000 tries on line 2003 [251] */
2006      if (sum < minreal) {   /* attention */
2007        tpr[(long)iba - (long)AA][(long)jba - (long)AA] = 0.0;
2008        td1[(long)iba - (long)AA][(long)jba - (long)AA] = 0.0;
2009        td2[(long)iba - (long)AA][(long)jba - (long)AA] = 0.0;
2010      } else {
2011        tpr[(long)iba - (long)AA][(long)jba - (long)AA] = sum;
2012        td1[(long)iba - (long)AA][(long)jba - (long)AA] = sumd1;
2013        td2[(long)iba - (long)AA][(long)jba - (long)AA] = sumd2;
2014      }
2015    }
2016  }
2017  /* IF negative THEN
2018     IF debugoptn THEN
2019     BEGIN
2020        write(output,' Trans. prob. 2 is negative !');
2021        writeln(output,'  arc =',arc:10:5);
2022     END; */
2023}  /* tdiffmtrx */
2024
2025
2026/*** MAKE TRANSITION PROBABILITY MATRIX ***/
2027Static Void maketransprob()
2028{
2029  /* make transition probability matrix */
2030  ndaryamity amatrix, bmatrix, cmatrix;
2031  amity iba, jba, kba;
2032
2033  if (readtoptn)
2034    readeigenv();
2035  else
2036    dataeigenv();
2037  if (writeoptn)
2038    printfreqdyhf();
2039
2040  if (readtoptn)
2041    printf("read trans. matrix\n");
2042  if (debugoptn)
2043    printeigen();
2044
2045  convermtrx(ev, amatrix);
2046  memcpy(cmatrix, amatrix, sizeof(ndaryamity));
2047  luinverse(cmatrix, bmatrix, 20L);
2048  /* IF debugoptn printnmtrx (bmatrix); */
2049  mproduct(amatrix, bmatrix, cmatrix, 20L, 20L, 20L);
2050  checkevector(cmatrix, 20L);
2051  /* IF debugoptn printnmtrx (cmatrix); */
2052  reconvermtrx(bmatrix, iev);
2053  for (iba = AA; (long)iba <= (long)VV; iba = (amity)((long)iba + 1)) {
2054    for (jba = AA; (long)jba <= (long)VV; jba = (amity)((long)jba + 1)) {
2055      for (kba = AA; (long)kba <= (long)VV; kba = (amity)((long)kba + 1))
2056        coefp[(long)iba - (long)AA][(long)jba - (long)AA]
2057          [(long)kba - (long)AA] = ev[(long)iba - (long)AA]
2058            [(long)kba - (long)AA] * iev[(long)kba - (long)AA]
2059            [(long)jba - (long)AA];
2060    }
2061  }
2062  for (kba = AA; (long)kba <= (long)VV; kba = (amity)((long)kba + 1))
2063    eval2[(long)kba - (long)AA] = eval[(long)kba - (long)AA] *
2064                                  eval[(long)kba - (long)AA];
2065  /*
2066     tprobmtrx ( 100.0, tprobt );
2067     printamtrx ( tprobt );
2068  */
2069}  /* maketransprob */
2070
2071
2072Static double randum(seed)
2073long *seed;
2074{
2075  /* random number generator -- slow but machine independent */
2076  long i, j, k, sum;
2077  longintty mult, newseed;
2078  double x;
2079
2080  mult[0] = 13;
2081  mult[1] = 24;
2082  mult[2] = 22;
2083  mult[3] = 6;
2084  for (i = 0; i <= 5; i++)
2085    newseed[i] = 0;
2086  for (i = 0; i <= 5; i++) {
2087    sum = newseed[i];
2088    k = i;
2089    if (i > 3)
2090      k = 3;
2091    for (j = 0; j <= k; j++)
2092      sum += mult[j] * seed[i - j];
2093    newseed[i] = sum;
2094    for (j = i; j <= 4; j++) {
2095      newseed[j + 1] += newseed[j] / 64;
2096      newseed[j] &= 63;
2097    }
2098  }
2099  memcpy(seed, newseed, sizeof(longintty));
2100  seed[5] &= 3;
2101  x = 0.0;
2102  for (i = 0; i <= 5; i++)
2103    x = x / 64.0 + seed[i];
2104  x /= 4.0;
2105  return x;
2106}  /* randum */
2107
2108
2109/*******************************************************/
2110/*****          ESTIMATE TREE TOPOLOGY             *****/
2111/*******************************************************/
2112
2113/**************************************/
2114/***       READ TREE STRUCTURE      ***/
2115/**************************************/
2116
2117/*** SERCH OF END CHARACTER ***/
2118Static Void serchchend()
2119{
2120  if (chs == ',' || chs == ')')
2121    return;
2122  do {
2123    if (P_eoln(seqfile)) {
2124      fscanf(seqfile, "%*[^\n]");
2125      getc(seqfile);
2126    }
2127    chs = getc(seqfile);
2128    if (chs == '\n')
2129      chs = ' ';
2130  } while (chs != ',' && chs != ')');
2131}  /* serchchend */
2132
2133
2134/*** NEXT CHARACTER ***/
2135Static Void nextchar(ch)
2136Char *ch;
2137{
2138  do {
2139    if (P_eoln(seqfile)) {
2140      fscanf(seqfile, "%*[^\n]");
2141      getc(seqfile);
2142    }
2143    *ch = getc(seqfile);
2144    if (*ch == '\n')
2145      *ch = ' ';
2146  } while (*ch == ' ');
2147}  /* nextchar */
2148
2149
2150/*** CREATE DATASTRUCTURE OF NODE ***/
2151Static Void clearnode(p)
2152node *p;
2153{
2154  long i, FORLIM;
2155
2156  p->isop = NULL;
2157  p->kinp = NULL;
2158  p->diverg = 1;
2159  p->number = 0;
2160  for (i = 0; i < maxname; i++)
2161    p->namesp[i] = ' ';
2162  FORLIM = numsp;
2163  for (i = 0; i < FORLIM; i++)
2164    p->descen[i] = 0;
2165  p->nadata = ami;
2166  p->UU.U0.lnga = 0.0;
2167}  /* clearnode */
2168
2169
2170/*** JUDGMENT SPEICIES NAME ***/
2171Static Void judgspname(tr, num)
2172tree *tr;
2173long *num;
2174{
2175  long ie;   /* number of name strings */
2176  namety namestr;   /* current species name of seqfile */
2177  boolean found;
2178
2179  for (ie = 0; ie < maxname; ie++)
2180    namestr[ie] = ' ';
2181  ie = 1;
2182  do {
2183    if (chs == '_')
2184      chs = ' ';
2185    namestr[ie - 1] = chs;
2186    if (P_eoln(seqfile)) {
2187      fscanf(seqfile, "%*[^\n]");
2188      getc(seqfile);
2189    }
2190    chs = getc(seqfile);
2191    if (chs == '\n')
2192      chs = ' ';
2193    ie++;
2194  } while (chs != ':' && chs != ',' && chs != ')' && ie <= maxname);
2195  *num = 1;
2196  do {
2197    found = true;
2198    for (ie = 0; ie < maxname; ie++)
2199      found = (found && namestr[ie] == tr->brnchp[*num]->namesp[ie]);
2200    if (!found)
2201      (*num)++;
2202  } while (!(*num > numsp || found));
2203  if (*num <= numsp)
2204    return;
2205  printf(" Cannot find species: ");
2206  for (ie = 0; ie < maxname; ie++)
2207    putchar(namestr[ie]);
2208  putchar('\n');
2209}  /* judgspname */
2210
2211
2212/*** ADD EXTERNAL NODE ***/
2213Static Void externalnode(tr, up)
2214tree *tr;
2215node *up;
2216{
2217  long num;   /* number of external nodes */
2218
2219  judgspname(tr, &num);
2220  tr->brnchp[num]->kinp = up;
2221  up->kinp = tr->brnchp[num];
2222  if (tr->startp->number > num)
2223    tr->startp = tr->brnchp[num];
2224  tr->brnchp[num]->UU.U0.lnga = 0.0;
2225  up->UU.U0.lnga = 0.0;
2226}  /* externalnode */
2227
2228
2229/*** ADD INTERNAL NODE ***/
2230Static Void internalnode(tr, up, ninode)
2231tree *tr;
2232node *up;
2233long *ninode;
2234{
2235  node *np;   /* new pointer to internal node */
2236  node *cp;   /* current pointer to internal node */
2237  long i, nd, dvg;
2238
2239  nextchar(&chs);
2240  if (chs != '(') {
2241    externalnode(tr, up);
2242    return;
2243  }
2244  (*ninode)++;
2245  nd = *ninode;
2246  if (freenode == NULL)
2247    np = (node *)Malloc(sizeof(node));
2248  else {
2249    np = freenode;
2250    freenode = np->isop;
2251    np->isop = NULL;
2252  }
2253  clearnode(np);
2254  np->isop = np;
2255  cp = np;
2256  np = NULL;
2257  cp->number = nd;
2258  tr->brnchp[*ninode] = cp;
2259  up->kinp = cp;
2260  cp->kinp = up;
2261  dvg = 0;
2262  while (chs != ')') {
2263    dvg++;
2264    if (freenode == NULL)
2265      np = (node *)Malloc(sizeof(node));
2266    else {
2267      np = freenode;
2268      freenode = np->isop;
2269      np->isop = NULL;
2270    }
2271    clearnode(np);
2272    np->isop = cp->isop;
2273    cp->isop = np;
2274    cp = np;
2275    np = NULL;
2276    cp->number = nd;
2277    internalnode(tr, cp, ninode);
2278    serchchend();
2279  }
2280  for (i = 0; i <= dvg; i++) {
2281    cp->diverg = dvg;
2282    cp = cp->isop;
2283  }
2284  nextchar(&chs);   /* */
2285}  /* internalnode */
2286
2287
2288/*** MAKE STRUCTURE OF TREE ***/
2289Static Void treestructure(tr)
2290tree *tr;
2291{
2292  long i, dvg, ninode;   /* number of internal nodes */
2293  node *np;   /* new pointer */
2294  node *cp;   /* current pointer to zero node(root) */
2295
2296  ninode = numsp;
2297  nextchar(&chs);
2298  if (chs == '(') {
2299    dvg = -1;
2300    while (chs != ')') {
2301      if (freenode == NULL)
2302        np = (node *)Malloc(sizeof(node));
2303      else {
2304        np = freenode;
2305        freenode = np->isop;
2306        np->isop = NULL;
2307      }
2308      clearnode(np);
2309      internalnode(tr, np, &ninode);
2310      dvg++;
2311      if (dvg == 0)
2312        np->isop = np;
2313      else {
2314        np->isop = cp->isop;
2315        cp->isop = np;
2316      }
2317      cp = np;
2318      np = NULL;
2319      serchchend();
2320    }
2321    tr->brnchp[0] = cp->isop;
2322    tr->startp = tr->brnchp[numsp];
2323    for (i = 0; i <= dvg; i++) {
2324      cp->diverg = dvg;
2325      cp = cp->isop;
2326    }
2327  }
2328  fscanf(seqfile, "%*[^\n]");
2329  getc(seqfile);
2330  ibrnch2 = ninode;
2331}  /* treestructure */
2332
2333
2334/*** MAKE STAR STRUCTURE OF TREE ***/
2335Static Void starstructure(tr)
2336tree *tr;
2337{
2338  long i;
2339  /* */
2340  long dvg;
2341  node *np;   /* new pointer */
2342  node *cp;   /* current pointer to zero node(root) */
2343  long FORLIM;
2344
2345  dvg = -1;
2346  FORLIM = numsp;
2347  for (i = 1; i <= FORLIM; i++) {
2348    if (freenode == NULL)
2349      np = (node *)Malloc(sizeof(node));
2350    else {
2351      np = freenode;
2352      freenode = np->isop;
2353      np->isop = NULL;
2354    }
2355    clearnode(np);
2356
2357    tr->brnchp[i]->kinp = np;
2358    np->kinp = tr->brnchp[i];
2359
2360    dvg++;
2361    if (dvg == 0)
2362      np->isop = np;
2363    else {
2364      np->isop = cp->isop;
2365      cp->isop = np;
2366    }
2367    cp = np;
2368    np = NULL;
2369  }
2370  tr->brnchp[0] = cp->isop;
2371  tr->startp = tr->brnchp[numsp];
2372  for (i = 0; i <= dvg; i++) {
2373    cp->diverg = dvg;
2374    cp = cp->isop;
2375  }
2376  ibrnch2 = numsp;
2377}  /* starstructure */
2378
2379
2380/* INSERT BRANCH IN TREE */
2381Static Void insertbranch(tr, cp1, cp2, bp1, bp2, np1, np2)
2382tree *tr;
2383node *cp1, *cp2, *bp1, *bp2, *np1, *np2;
2384{
2385  long i, num, dvg;
2386  node *ap;
2387
2388  i = cp1->number;
2389  if (cp1 == tr->brnchp[i]) {
2390    if (i != 0) {
2391      ap = np1;
2392      np1 = np2;
2393      np2 = ap;
2394    } else
2395      tr->brnchp[i] = np2;
2396  }
2397  /*  np2^.number := bp1^.number;
2398      cp1^.number := np1^.number;
2399      cp2^.number := np1^.number;  */
2400  bp1->isop = np2;
2401  if (cp1 == bp2)
2402    np2->isop = cp2->isop;
2403  else
2404    np2->isop = cp1->isop;
2405  if (cp1 != bp2)
2406    bp2->isop = cp2->isop;
2407  np1->isop = cp1;
2408  if (cp1 != bp2)
2409    cp1->isop = cp2;
2410  cp2->isop = np1;
2411  tr->brnchp[ibrnch2]->kinp->number = tr->brnchp[ibrnch2]->kinp->isop->number;
2412  ap = tr->brnchp[ibrnch2];
2413  num = tr->brnchp[ibrnch2]->number;
2414  do {
2415    ap = ap->isop;
2416    ap->number = num;
2417  } while (ap->isop != tr->brnchp[ibrnch2]);
2418  dvg = 0;
2419  ap = np1;
2420  do {
2421    ap = ap->isop;
2422    dvg++;
2423  } while (ap->isop != np1);
2424  ap = np1;
2425  do {
2426    ap = ap->isop;
2427    ap->diverg = dvg;
2428  } while (ap != np1);
2429  dvg = 0;
2430  ap = np2;
2431  do {
2432    ap = ap->isop;
2433    dvg++;
2434  } while (ap->isop != np2);
2435  ap = np2;
2436  do {
2437    ap = ap->isop;
2438    ap->diverg = dvg;
2439  } while (ap != np2);
2440}  /* insertbranch */
2441
2442
2443/* DELETE BRANCH IN TREE */
2444Static Void deletebranch(tr, cnode, cp1, cp2, bp1, bp2, np1, np2)
2445tree *tr;
2446long cnode;
2447node *cp1, *cp2, *bp1, *bp2, *np1, *np2;
2448{
2449  /* i, */
2450  long dvg;
2451  node *ap;
2452
2453  if (cnode != 0) {
2454    if (tr->brnchp[cnode] == cp1) {
2455      ap = np1;
2456      np1 = np2;
2457      np2 = ap;
2458    }
2459  } else {
2460    if (tr->brnchp[cnode] == np2)
2461      tr->brnchp[cnode] = cp1;
2462  }
2463  /* i := np2^.number;
2464     IF i = 0 THEN
2465        IF np2 = tr.brnchp[i] THEN tr.brnchp[i] := cp1
2466     ELSE
2467        IF cp1 = tr.brnchp[i] THEN
2468        BEGIN
2469           ap  := np1;
2470           np1 := np2;
2471           np2 := ap;
2472        END; */
2473  /* cp1^.number := np2^.number;
2474     cp2^.number := np2^.number; */
2475  if (cp1 == bp2)
2476    cp2->isop = np2->isop;
2477  else
2478    cp2->isop = bp2->isop;
2479  if (cp1 != bp2)
2480    cp1->isop = np2->isop;
2481  np1->isop = NULL;
2482  if (cp1 != bp2)
2483    bp2->isop = cp2;
2484  np2->isop = NULL;
2485  bp1->isop = cp1;
2486  dvg = 0;
2487  ap = cp1;
2488  do {
2489    ap = ap->isop;
2490    dvg++;
2491  } while (ap->isop != cp1);
2492  ap = cp1;
2493  do {
2494    ap = ap->isop;
2495    ap->diverg = dvg;
2496    ap->number = cnode;
2497  } while (ap != cp1);
2498  np1->diverg = 0;
2499  np2->diverg = 0;
2500}  /* deletebranch */
2501
2502
2503/*** PRINT STRUCTURE OF TREE ***/
2504Static Void printcurtree(tr)
2505tree *tr;
2506{
2507  node *ap;
2508  long i, j, k, num, FORLIM, FORLIM1;
2509
2510  printf("\nStructure of Tree\n");
2511  printf("%7s%5s%5s%8s%9s%11s%13s\n",
2512         "number", "kinp", "isop", "diverg", "namesp", "length",
2513         "descendant");
2514  FORLIM = numsp;
2515  for (i = 1; i <= FORLIM; i++) {
2516    ap = tr->brnchp[i];
2517    printf("%5ld", ap->number);
2518    printf("%6ld", ap->kinp->number);
2519    if (ap->isop == NULL)
2520      printf("%6s", "nil");
2521    else
2522      printf("%6ld", ap->isop->number);
2523    printf("%7ld", ap->diverg);
2524    printf("%3c", ' ');
2525    for (j = 0; j < maxname; j++)
2526      putchar(ap->namesp[j]);
2527    printf("%8.3f", ap->UU.U0.lnga);
2528    printf("%3c", ' ');
2529    FORLIM1 = numsp;
2530    for (j = 0; j < FORLIM1; j++)
2531      printf("%2ld", ap->descen[j]);
2532    putchar('\n');
2533  }
2534  FORLIM = ibrnch2 + 1;
2535  for (i = ibrnch1; i <= FORLIM; i++) {
2536    if (i == ibrnch2 + 1)
2537      num = 0;
2538    else
2539      num = i;
2540    k = 0;
2541    ap = tr->brnchp[num];
2542    do {
2543      k++;
2544      if (ap->number == tr->brnchp[num]->number && k > 1)
2545        printf("%5c", '.');
2546      else
2547        printf("%5ld", ap->number);
2548      printf("%6ld", ap->kinp->number);
2549      if (ap->isop->number == tr->brnchp[num]->number && k > 0)   /* 1 */
2550        printf("%6c", '.');
2551      else
2552        printf("%6ld", ap->isop->number);
2553      printf("%7ld", ap->diverg);
2554      printf("%3c", ' ');
2555      for (j = 1; j <= maxname; j++)
2556        putchar(' ');
2557      printf("%8.3f", ap->UU.U0.lnga);
2558      printf("%3c", ' ');
2559      FORLIM1 = numsp;
2560      for (j = 0; j < FORLIM1; j++)
2561        printf("%2ld", ap->descen[j]);
2562      putchar('\n');
2563      ap = ap->isop;
2564    } while (ap != tr->brnchp[num]);
2565  }
2566}  /* printcurtree */
2567
2568
2569/**************************************/
2570/***  ESTIMATE LIKELIHOOD OF TREE   ***/
2571/**************************************/
2572
2573Static Void initdescen(tr)
2574tree *tr;
2575{
2576  node *ap;
2577  long i, j, k, num, FORLIM, FORLIM1, FORLIM2;
2578
2579  FORLIM = numsp;
2580  for (i = 1; i <= FORLIM; i++) {
2581    ap = tr->brnchp[i];
2582    FORLIM1 = numsp;
2583    for (j = 0; j < FORLIM1; j++)
2584      ap->descen[j] = 0;
2585    ap->descen[i - 1] = 1;
2586  }
2587  FORLIM = ibrnch2 + 1;
2588  for (i = ibrnch1; i <= FORLIM; i++) {
2589    if (i == ibrnch2 + 1)
2590      num = 0;
2591    else
2592      num = i;
2593    ap = tr->brnchp[num];
2594    FORLIM1 = tr->brnchp[num]->diverg + 1;
2595    for (k = 1; k <= FORLIM1; k++) {
2596      FORLIM2 = numsp;
2597      for (j = 0; j < FORLIM2; j++)
2598        ap->descen[j] = 0;
2599      ap = ap->isop;
2600    }
2601  }
2602}  /* initdescen */
2603
2604
2605Static Void initnodeA(p)
2606node *p;
2607{
2608  node *cp;
2609  long i, n;
2610  double sumprb;
2611  amity ba;
2612  long FORLIM, FORLIM2;
2613
2614  if (p->isop == NULL)   /* TIP */
2615    return;
2616  cp = p->isop;
2617  FORLIM = p->diverg;
2618  for (n = 1; n <= FORLIM; n++) {
2619    initnodeA(cp->kinp);
2620    cp = cp->isop;
2621  }
2622  FORLIM = endptrn;
2623  for (i = 0; i < FORLIM; i++) {
2624    for (ba = AA; (long)ba <= (long)VV; ba = (amity)((long)ba + 1)) {
2625      sumprb = 0.0;
2626      cp = p->isop;
2627      FORLIM2 = p->diverg;
2628      for (n = 1; n <= FORLIM2; n++) {
2629        sumprb += cp->kinp->UU.U0.prba[i][(long)ba - (long)AA];
2630        cp = cp->isop;
2631      }
2632      if (sumprb != 0.0)
2633        p->UU.U0.prba[i][(long)ba - (long)AA] = sumprb / p->diverg;
2634      else
2635        p->UU.U0.prba[i][(long)ba - (long)AA] = 0.0;
2636    }
2637  }
2638  p->UU.U0.lnga = 10.0;   /* attention */
2639  p->kinp->UU.U0.lnga = 10.0;   /* attention */
2640  FORLIM = numsp;
2641  for (i = 0; i < FORLIM; i++) {
2642    cp = p->isop;
2643    FORLIM2 = p->diverg;
2644    for (n = 1; n <= FORLIM2; n++) {
2645      if (cp->kinp->descen[i] > 0)
2646        p->descen[i] = 1;
2647      cp = cp->isop;
2648    }
2649  }
2650  /*  IF debugoptn THEN
2651      BEGIN
2652         write(output,p^.kinp^.number:5,'-':2,p^.number:3,' ':3);
2653         '(':2,q^.number:3,r^.number:3,')':2,
2654         FOR i := 1 TO numsp DO  write(output,p^.descen[i]:2);
2655         writeln(output);
2656      END;  */
2657}  /* initnodeA */
2658
2659
2660Static Void initbranch(p)
2661node *p;
2662{
2663  long n;
2664  node *cp;
2665  long FORLIM;
2666
2667  if (p->isop == NULL) {   /* TIP */
2668    initnodeA(p->kinp);
2669    return;
2670  }
2671  cp = p->isop;
2672  FORLIM = p->diverg;
2673  for (n = 1; n <= FORLIM; n++) {
2674    initbranch(cp->kinp);
2675    cp = cp->isop;
2676  }
2677}  /* initbranch */
2678
2679
2680Static Void evaluateA(tr)
2681tree *tr;
2682{
2683  double arc, lkl, sum, prod, lnlkl;
2684  long i;
2685  node *p, *q;
2686  aryamity xa1, xa2;
2687  amity ai, aj;
2688  daryamity tprobe;
2689  long FORLIM;
2690
2691  p = tr->startp;
2692  q = p->kinp;
2693  arc = p->UU.U0.lnga;
2694  if (arc < minarc)
2695    arc = minarc;
2696  tprobmtrx(arc, tprobe);
2697  lkl = 0.0;
2698  FORLIM = endptrn;
2699  for (i = 0; i < FORLIM; i++) {
2700    memcpy(xa1, p->UU.U0.prba[i], sizeof(aryamity));
2701    memcpy(xa2, q->UU.U0.prba[i], sizeof(aryamity));
2702    sum = 0.0;
2703    for (ai = AA; (long)ai <= (long)VV; ai = (amity)((long)ai + 1)) {
2704      prod = freqdyhf[(long)ai - (long)AA] * xa1[(long)ai - (long)AA];
2705      for (aj = AA; (long)aj <= (long)VV; aj = (amity)((long)aj + 1))
2706        sum += prod * tprobe[(long)ai - (long)AA]
2707               [(long)aj - (long)AA] * xa2[(long)aj - (long)AA];
2708    }
2709    if (sum > 0.0)
2710      lnlkl = log(sum);
2711    else
2712      lnlkl = 0.0;
2713    lkl += lnlkl * weight[i];
2714    lklhdtrpt[notree][i] = lnlkl;
2715  }
2716  tr->lklhd = lkl;
2717  tr->aic = ibrnch2 * 2 - 2.0 * lkl;
2718}  /* evaluateA */
2719
2720
2721Static Void sublklhdA(p)
2722node *p;
2723{
2724  long i, n;
2725  double arc;
2726  node *cp, *sp;
2727  amity ai;
2728  daryamity tprob;
2729  long FORLIM, FORLIM1;
2730
2731  cp = p->isop;
2732  FORLIM = p->diverg;
2733  for (n = 1; n <= FORLIM; n++) {
2734    sp = cp->kinp;
2735    arc = sp->UU.U0.lnga;
2736    tprobmtrx(arc, tprob);
2737    FORLIM1 = endptrn;
2738    for (i = 0; i < FORLIM1; i++) {
2739      for (ai = AA; (long)ai <= (long)VV; ai = (amity)((long)ai + 1)) {
2740        if (n == 1)
2741          p->UU.U0.prba[i][(long)ai - (long)AA] = 1.0;
2742        if (p->UU.U0.prba[i][(long)ai - (long)AA] < minreal)
2743          p->UU.U0.prba[i][(long)ai - (long)AA] = 0.0;
2744        else  /* attention */
2745          p->UU.U0.prba[i][(long)ai - (long)AA] *= tprob[(long)ai - (long)AA]
2746                [0] * sp->UU.U0.prba[i][0] + tprob[(long)ai - (long)AA][(long)
2747                      RR - (long)AA] * sp->UU.U0.prba[i]
2748                [(long)RR - (long)AA] + tprob[(long)ai - (long)AA][(long)NN -
2749                    (long)AA] * sp->UU.U0.prba[i][(long)NN - (long)AA] +
2750              tprob[(long)ai - (long)AA][(long)DD - (long)AA] *
2751                sp->UU.U0.prba[i][(long)DD - (long)AA] + tprob[(long)ai -
2752                    (long)AA][(long)CC - (long)AA] * sp->UU.U0.prba[i][(long)
2753                      CC - (long)AA] + tprob[(long)ai - (long)AA][(long)QQ -
2754                    (long)AA] * sp->UU.U0.prba[i][(long)QQ - (long)AA] +
2755              tprob[(long)ai - (long)AA][(long)EE - (long)AA] *
2756                sp->UU.U0.prba[i][(long)EE - (long)AA] + tprob[(long)ai -
2757                    (long)AA][(long)GG - (long)AA] * sp->UU.U0.prba[i][(long)
2758                      GG - (long)AA] + tprob[(long)ai - (long)AA][(long)HH -
2759                    (long)AA] * sp->UU.U0.prba[i][(long)HH - (long)AA] +
2760              tprob[(long)ai - (long)AA][(long)II - (long)AA] *
2761                sp->UU.U0.prba[i][(long)II - (long)AA] + tprob[(long)ai -
2762                    (long)AA][(long)LL - (long)AA] * sp->UU.U0.prba[i][(long)
2763                      LL - (long)AA] + tprob[(long)ai - (long)AA][(long)KK -
2764                    (long)AA] * sp->UU.U0.prba[i][(long)KK - (long)AA] +
2765              tprob[(long)ai - (long)AA][(long)MM - (long)AA] *
2766                sp->UU.U0.prba[i][(long)MM - (long)AA] +
2767              tprob[(long)ai - (long)AA][(long)FF - (long)AA] *
2768                sp->UU.U0.prba[i]
2769                [(long)FF - (long)AA] + tprob[(long)ai - (long)AA]
2770                [(long)PP_ - (long)AA] * sp->UU.U0.prba[i]
2771                [(long)PP_ - (long)AA] + tprob[(long)ai - (long)AA]
2772                [(long)SS - (long)AA] * sp->UU.U0.prba[i]
2773                [(long)SS - (long)AA] + tprob[(long)ai - (long)AA]
2774                [(long)TT - (long)AA] * sp->UU.U0.prba[i]
2775                [(long)TT - (long)AA] + tprob[(long)ai - (long)AA]
2776                [(long)WW - (long)AA] * sp->UU.U0.prba[i]
2777                [(long)WW - (long)AA] + tprob[(long)ai - (long)AA]
2778                [(long)YY - (long)AA] * sp->UU.U0.prba[i]
2779                [(long)YY - (long)AA] + tprob[(long)ai - (long)AA]
2780                [(long)VV - (long)AA] * sp->UU.U0.prba[i]
2781                [(long)VV - (long)AA];
2782/* p2c: protml.pas, line 2226: Note:
2783 * Line breaker spent 0.0+1.00 seconds, 5000 tries on line 2781 [251] */
2784      }
2785    }
2786    cp = cp->isop;
2787  }
2788}  /* sublklhdA */
2789
2790
2791Static Void branchlengthA(p, it)
2792node *p;
2793long *it;
2794{
2795  long i, numloop;
2796  double arc, arcold, sum, sumd1, sumd2, lkl, lkld1, lkld2, prod1, prod2;
2797  boolean done;
2798  node *q;
2799  amity ai, aj;
2800  daryamity tprobl, tdiff1, tdiff2;
2801  long FORLIM;
2802
2803  q = p->kinp;
2804  arc = p->UU.U0.lnga;
2805  done = false;
2806  *it = 0;
2807  if (numsm < 3)
2808    numloop = 3;
2809  else
2810    numloop = maxiterat;
2811
2812  while (*it < numloop && !done) {
2813    (*it)++;
2814    if (arc < minarc)
2815      arc = minarc;
2816    if (arc > maxarc)   /* attention */
2817      arc = minarc;
2818    tdiffmtrx(arc, tprobl, tdiff1, tdiff2);
2819    lkl = 0.0;
2820    lkld1 = 0.0;
2821    lkld2 = 0.0;
2822    FORLIM = endptrn;
2823    for (i = 0; i < FORLIM; i++) {
2824      sum = 0.0;
2825      sumd1 = 0.0;
2826      sumd2 = 0.0;
2827      for (ai = AA; (long)ai <= (long)VV; ai = (amity)((long)ai + 1)) {
2828        prod1 = freqdyhf[(long)ai - (long)AA] * p->UU.U0.prba[i]
2829                [(long)ai - (long)AA];
2830        if (prod1 < minreal)   /* attention */
2831          prod1 = 0.0;
2832        for (aj = AA; (long)aj <= (long)VV; aj = (amity)((long)aj + 1)) {
2833          prod2 = prod1 * q->UU.U0.prba[i][(long)aj - (long)AA];
2834          if (prod2 < minreal)   /* attention */
2835            prod2 = 0.0;
2836          sum += prod2 * tprobl[(long)ai - (long)AA][(long)aj - (long)AA];
2837          sumd1 += prod2 * tdiff1[(long)ai - (long)AA][(long)aj - (long)AA];
2838          sumd2 += prod2 * tdiff2[(long)ai - (long)AA][(long)aj - (long)AA];
2839        }
2840      }
2841      if (sum > minreal) {  /* attention */
2842        lkl += log(sum) * weight[i];
2843        lkld1 += sumd1 / sum * weight[i];
2844        if (sum * sum > minreal)
2845          lkld2 += (sumd2 * sum - sumd1 * sumd1) / sum / sum * weight[i];
2846      } else {
2847        if (debugoptn)
2848          printf(" *check branchlength1*%4ld%4ld\n", p->number, i + 1);
2849      }  /* attention */
2850    }
2851    arcold = arc;
2852    if (lkld1 != lkld2)
2853      arc -= lkld1 / lkld2;
2854    else {
2855      arcold = arc + epsilon * 0.1;
2856      if (debugoptn)
2857        printf(" *check branchlength2*");
2858    }
2859    if (arc > maxarc)   /* attention */
2860      arc = minarc;
2861    done = (fabs(arcold - arc) < epsilon);
2862    /* IF debugoptn THEN writeln(output,' ':10,arc:10:5,
2863       arcold:10:5, -(lkld1/lkld2):10:5); */
2864  }
2865  smoothed = (smoothed && done);
2866  if (arc < minarc)
2867    arc = minarc;
2868  if (arc > maxarc)   /* attention */
2869    arc = minarc;
2870  p->UU.U0.lnga = arc;
2871  q->UU.U0.lnga = arc;
2872}  /* branchlengthA */
2873
2874
2875Static Void printupdate(tr, p, vold, it)
2876tree *tr;
2877node *p;
2878double vold;
2879long it;
2880{
2881  if (p == tr->startp->kinp)
2882    printf("%4ld", numsm);
2883  else
2884    printf("%4c", ' ');
2885  if (p->isop == NULL)   /* TIP */
2886    printf("%12c", ' ');
2887  else {
2888    printf("%4c%3ld", '(', p->isop->kinp->number);
2889    printf("%3ld )", p->isop->isop->kinp->number);
2890  }
2891  printf("%3ld -%3ld", p->number, p->kinp->number);
2892  if (p->kinp->isop == NULL)   /* TIP */
2893    printf("%10c", ' ');
2894  else {
2895    printf(" (%3ld", p->kinp->isop->kinp->number);
2896    printf("%3ld )", p->kinp->isop->isop->kinp->number);
2897  }
2898  printf("%2c%10.5f", ' ', p->UU.U0.lnga);
2899  printf(" (%10.5f )", p->UU.U0.lnga - vold);
2900  printf("%4ld\n", it);
2901}  /* printupdate */
2902
2903
2904Static Void updateA(tr, p, it)
2905tree *tr;
2906node *p;
2907long *it;
2908{
2909  double vold;
2910
2911  vold = p->UU.U0.lnga;
2912  if (p->isop != NULL)   /* TIP */
2913    sublklhdA(p);
2914  if (p->kinp->isop != NULL)   /* TIP */
2915    sublklhdA(p->kinp);
2916  branchlengthA(p, it);
2917  /* IF debugoptn AND ((numsm = 1) OR (numsm = maxsmooth)) THEN */
2918  if (debugoptn && numsm == 1)
2919    printupdate(tr, p, vold, *it);
2920}  /* updateA */
2921
2922
2923Static Void smooth2(tr, numit)
2924tree *tr;
2925long *numit;
2926{
2927  long n, it, FORLIM;
2928
2929  FORLIM = numsp;
2930  for (n = 1; n <= FORLIM; n++) {
2931    updateA(tr, tr->brnchp[n], &it);
2932    *numit += it;
2933  }
2934  FORLIM = ibrnch1;
2935  for (n = ibrnch2; n >= FORLIM; n--) {
2936    updateA(tr, tr->brnchp[n], &it);
2937    *numit += it;
2938  }
2939}  /* smooth */
2940
2941
2942Static Void smooth(tr, p, numit)
2943tree *tr;
2944node *p;
2945long *numit;
2946{
2947  long n, it;
2948  double vold;
2949  node *cp;
2950  long FORLIM;
2951
2952  vold = p->UU.U0.lnga;
2953  if (p->isop != NULL)   /* TIP */
2954    sublklhdA(p);
2955  if (p->kinp->isop != NULL)   /* TIP */
2956    sublklhdA(p->kinp);
2957  branchlengthA(p, &it);
2958  *numit += it;
2959  if (debugoptn && (numsm == 1 || numsm == maxsmooth))
2960    printupdate(tr, p, vold, it);
2961  if (p->isop == NULL)   /* TIP */
2962    return;
2963  /* smooth ( tr,p^.isop^.kinp,numit );
2964     smooth ( tr,p^.isop^.isop^.kinp,numit ); */
2965  cp = p->isop;
2966  FORLIM = p->diverg;
2967  for (n = 1; n <= FORLIM; n++) {
2968    smooth(tr, cp->kinp, numit);
2969    cp = cp->isop;
2970  }
2971}  /* smooth */
2972
2973
2974Static Void printsmooth(tr, numit)
2975tree *tr;
2976long numit;
2977{
2978  long i, FORLIM;
2979
2980  if (debugoptn)
2981    return;
2982  printf(" %3ld", numsm);
2983  printf("%4ld", numit);
2984  FORLIM = ibrnch2;
2985  for (i = 1; i <= FORLIM; i++)
2986    printf("%5.1f", tr->brnchp[i]->UU.U0.lnga);
2987  putchar('\n');
2988}  /* printsmooth */
2989
2990
2991Static Void leastsquares(am_, ym)
2992rnodety *am_;
2993double *ym;
2994{
2995  dnonoty am;
2996  long i, j, k;
2997  double pivot, element;
2998  dnonoty im;
2999  long FORLIM, FORLIM1, FORLIM2;
3000
3001  memcpy(am, am_, sizeof(dnonoty));
3002  FORLIM = ibrnch2;
3003  for (i = 1; i <= FORLIM; i++) {
3004    FORLIM1 = ibrnch2;
3005    for (j = 1; j <= FORLIM1; j++) {
3006      if (i == j)
3007        im[i - 1][j - 1] = 1.0;
3008      else
3009        im[i - 1][j - 1] = 0.0;
3010    }
3011  }
3012  FORLIM = ibrnch2;
3013  for (k = 0; k < FORLIM; k++) {
3014    pivot = am[k][k];
3015    ym[k] /= pivot;
3016    FORLIM1 = ibrnch2;
3017    for (j = 0; j < FORLIM1; j++) {
3018      am[k][j] /= pivot;
3019      im[k][j] /= pivot;
3020    }
3021    FORLIM1 = ibrnch2;
3022    for (i = 0; i < FORLIM1; i++) {
3023      if (k + 1 != i + 1) {
3024        element = am[i][k];
3025        ym[i] -= element * ym[k];
3026        FORLIM2 = ibrnch2;
3027        for (j = 0; j < FORLIM2; j++) {
3028          am[i][j] -= element * am[k][j];
3029          im[i][j] -= element * im[k][j];
3030        }
3031      }
3032    }
3033  }
3034}  /* leastsquares */
3035
3036
3037Static Void initlength(tr)
3038tree *tr;
3039{
3040  long ia, j1, j2, k, na, n1, n2;
3041  double suma, sumy;
3042  spity des;
3043  rpairty dfpair, ymt;
3044  rnodety atymt;
3045  dpanoty amt;
3046  dnopaty atmt;
3047  dnonoty atamt;
3048  long FORLIM, FORLIM1, FORLIM2, FORLIM3;
3049
3050  FORLIM = ibrnch2;
3051  for (ia = 1; ia <= FORLIM; ia++) {
3052    memcpy(des, tr->brnchp[ia]->descen, sizeof(spity));
3053    na = 0;
3054    FORLIM1 = numsp;
3055    for (j1 = 1; j1 < FORLIM1; j1++) {
3056      FORLIM2 = numsp;
3057      for (j2 = j1 + 1; j2 <= FORLIM2; j2++) {
3058        na++;
3059        /* writeln(output,' ',ia:3,j1:3,j2:3,na:3); */
3060        if (des[j1 - 1] != des[j2 - 1])
3061          amt[na - 1][ia - 1] = 1.0;
3062        else
3063          amt[na - 1][ia - 1] = 0.0;
3064        if (ia == 1) {
3065          dfpair[na - 1] = 0.0;
3066          FORLIM3 = endsite;
3067          for (k = 0; k < FORLIM3; k++) {
3068            if (chsequen[j1][k] != chsequen[j2][k])
3069              dfpair[na - 1] += 1.0;
3070          }
3071        }
3072        if (dfpair[na - 1] > 0.0)
3073          ymt[na - 1] = -log(1.0 - dfpair[na - 1] / endsite);
3074        else
3075          ymt[na - 1] = -log(1.0);
3076      }
3077    }
3078  }
3079
3080  if (debugoptn) {
3081    putchar('\n');
3082    FORLIM = numpair;
3083    for (na = 0; na < FORLIM; na++) {
3084      printf(" %3ld ", na + 1);
3085      FORLIM1 = ibrnch2;
3086      for (ia = 0; ia < FORLIM1; ia++)
3087        printf("%3ld", (long)amt[na][ia]);
3088      printf("%8.3f%5ld\n", ymt[na], (long)dfpair[na]);
3089    }
3090  }
3091
3092  FORLIM = numpair;
3093  for (ia = 0; ia < FORLIM; ia++) {
3094    FORLIM1 = ibrnch2;
3095    for (na = 0; na < FORLIM1; na++)
3096      atmt[na][ia] = amt[ia][na];
3097  }
3098  FORLIM = ibrnch2;
3099  for (n1 = 0; n1 < FORLIM; n1++) {
3100    FORLIM1 = ibrnch2;
3101    for (n2 = 0; n2 < FORLIM1; n2++) {
3102      suma = 0.0;
3103      FORLIM2 = numpair;
3104      for (ia = 0; ia < FORLIM2; ia++)
3105        suma += atmt[n1][ia] * amt[ia][n2];
3106      atamt[n1][n2] = suma;
3107    }
3108  }
3109  FORLIM = ibrnch2;
3110  for (n1 = 0; n1 < FORLIM; n1++) {
3111    sumy = 0.0;
3112    FORLIM1 = numpair;
3113    for (ia = 0; ia < FORLIM1; ia++)
3114      sumy += atmt[n1][ia] * ymt[ia];
3115    atymt[n1] = sumy;
3116  }
3117
3118  if (debugoptn) {
3119    putchar('\n');
3120    FORLIM = ibrnch2;
3121    for (n1 = 1; n1 <= FORLIM; n1++) {
3122      printf(" %3ld ", n1);
3123      FORLIM1 = ibrnch2;
3124      for (n2 = 0; n2 < FORLIM1; n2++)
3125        printf("%3ld", (long)atamt[n1 - 1][n2]);
3126      printf("%8.3f\n", atymt[n1 - 1]);
3127    }
3128  }
3129
3130  leastsquares(atamt, atymt);
3131  FORLIM = ibrnch2;
3132  for (na = 0; na < FORLIM; na++)
3133    atymt[na] = 100.0 * atymt[na];
3134
3135  if (!debugoptn) {
3136    if (writeoptn) {
3137      printf("\n%5s", " arc");
3138      printf("%4s", "  it");
3139      FORLIM = ibrnch2;
3140      for (na = 1; na <= FORLIM; na++)
3141        printf("%3ld%2c", na, ' ');
3142      printf("\n%4c", '0');
3143      printf("%4c", '0');
3144      FORLIM = ibrnch2;
3145      for (na = 0; na < FORLIM; na++)
3146        printf("%5.1f", atymt[na]);
3147      putchar('\n');
3148    }
3149  }
3150
3151  FORLIM = ibrnch2;
3152  for (na = 1; na <= FORLIM; na++) {
3153    tr->brnchp[na]->UU.U0.lnga = atymt[na - 1];
3154    tr->brnchp[na]->kinp->UU.U0.lnga = atymt[na - 1];
3155  }
3156
3157}  /* initlength */
3158
3159
3160Static Void judgconverg(oldarcs, newarcs, cnvrg)
3161double *oldarcs, *newarcs;
3162boolean *cnvrg;
3163{
3164  /* judgment of convergence */
3165  long i;
3166  boolean same;
3167  long FORLIM;
3168
3169  same = true;
3170  FORLIM = ibrnch2;
3171  for (i = 0; i < FORLIM; i++) {
3172    if (fabs(newarcs[i] - oldarcs[i]) > epsilon)
3173      same = false;
3174  }
3175  *cnvrg = same;
3176}  /* judgconverg */
3177
3178
3179Static Void variance(tr)
3180tree *tr;
3181{
3182  long k, m;
3183  aryamity xa1, xa2;
3184  amity ai, aj;
3185  double arcm, alkl, lkl, lkl2, sarc, sarc2, sblklhd, ld1, ld2, prod, vlkl;
3186  lengthty varc, varc2, blklhdm;
3187  daryamity tprobm, tdifm1, tdifm2;
3188  long FORLIM;
3189  double TEMP;
3190  long FORLIM1;
3191
3192  alkl = tr->lklhd / endsite;
3193  vlkl = 0.0;
3194  FORLIM = endptrn;
3195  for (k = 0; k < FORLIM; k++) {
3196    TEMP = lklhdtrpt[notree][k] - alkl;
3197    vlkl += TEMP * TEMP * weight[k];
3198  }
3199  tr->vrilkl = vlkl;
3200
3201  FORLIM = ibrnch2;
3202  for (m = 1; m <= FORLIM; m++) {
3203    arcm = tr->brnchp[m]->UU.U0.lnga;
3204    tdiffmtrx(arcm, tprobm, tdifm1, tdifm2);
3205    sarc = 0.0;
3206    sarc2 = 0.0;
3207    sblklhd = 0.0;
3208    FORLIM1 = endptrn;
3209    for (k = 0; k < FORLIM1; k++) {
3210      memcpy(xa1, tr->brnchp[m]->UU.U0.prba[k], sizeof(aryamity));
3211      memcpy(xa2, tr->brnchp[m]->kinp->UU.U0.prba[k], sizeof(aryamity));
3212      lkl2 = 0.0;
3213      ld1 = 0.0;
3214      ld2 = 0.0;
3215      for (ai = AA; (long)ai <= (long)VV; ai = (amity)((long)ai + 1)) {
3216        prod = freqdyhf[(long)ai - (long)AA] * xa1[(long)ai - (long)AA];
3217        if (prod < minreal)   /* attention */
3218          prod = 0.0;
3219        for (aj = AA; (long)aj <= (long)VV; aj = (amity)((long)aj + 1)) {
3220          if (xa2[(long)aj - (long)AA] < minreal)   /* attention */
3221            xa2[(long)aj - (long)AA] = 0.0;
3222          lkl2 += prod * tprobm[(long)ai - (long)AA]
3223                  [(long)aj - (long)AA] * xa2[(long)aj - (long)AA];
3224          ld1 += prod * tdifm1[(long)ai - (long)AA]
3225                 [(long)aj - (long)AA] * xa2[(long)aj - (long)AA];
3226          ld2 += prod * tdifm2[(long)ai - (long)AA]
3227                 [(long)aj - (long)AA] * xa2[(long)aj - (long)AA];
3228        }
3229      }
3230      lkl = exp(lklhdtrpt[notree][k]);
3231      if (lkl * lkl > minreal)   /* attention */
3232        sarc += (ld2 * lkl - ld1 * ld1) / lkl / lkl * weight[k];
3233      if (lkl2 * lkl2 > minreal)   /* attention */
3234        sarc2 += (ld2 * lkl2 - ld1 * ld1) / lkl2 / lkl2 * weight[k];
3235      if (lkl2 > minreal)   /* attention */
3236        sblklhd += log(lkl2) * weight[k];
3237    }
3238    varc[m - 1] = sarc;
3239    varc2[m - 1] = sarc2;
3240    blklhdm[m - 1] = sblklhd;
3241
3242    /*          IF debugoptn THEN BEGIN
3243                   write  (output,m:4);
3244                   FOR m := 1 TO ibrnch2 DO write(output,sarc[m]:9:6);
3245                   writeln(output);
3246                END; */
3247
3248  }
3249  FORLIM = ibrnch2;
3250  for (m = 0; m < FORLIM; m++)
3251    varc[m] = 1.0 / varc[m];
3252  FORLIM = ibrnch2;
3253  for (m = 0; m < FORLIM; m++)
3254    varc2[m] = 1.0 / varc2[m];
3255  memcpy(tr->vrilnga, varc, sizeof(lengthty));
3256  memcpy(tr->vrilnga2, varc2, sizeof(lengthty));
3257  memcpy(tr->blklhd, blklhdm, sizeof(lengthty));
3258}  /* variance */
3259
3260
3261Static Void manuvaluate(tr)
3262tree *tr;
3263{
3264  long n, it;
3265  lengthty newarcs, oldarcs;
3266  boolean cnvrg;
3267  long FORLIM;
3268
3269  if (writeoptn)
3270    putchar('\n');
3271  if (debugoptn)
3272    printf(" MANUALVALUATE\n");
3273  /*  IF debugoptn THEN printcurtree ( tr );  */
3274  if (debugoptn)
3275    putchar('\n');
3276  initdescen(tr);
3277  initbranch(tr->startp);
3278  initbranch(tr->startp->kinp);
3279  /*  IF debugoptn THEN printcurtree ( tr );  */
3280  initlength(tr);
3281  if (debugoptn)
3282    printcurtree(tr);
3283
3284  if (debugoptn)
3285    printf(" SMOOTH\n");
3286  cnvrg = false;
3287  FORLIM = ibrnch2;
3288  for (n = 1; n <= FORLIM; n++)
3289    newarcs[n - 1] = tr->brnchp[n]->UU.U0.lnga;
3290  numnw = 0;
3291  numsm = 0;
3292  do {
3293    numsm++;
3294    it = 0;
3295    smooth(tr, tr->startp->kinp, &it);
3296    memcpy(oldarcs, newarcs, sizeof(lengthty));
3297    FORLIM = ibrnch2;
3298    for (n = 1; n <= FORLIM; n++)
3299      newarcs[n - 1] = tr->brnchp[n]->UU.U0.lnga;
3300    judgconverg(oldarcs, newarcs, &cnvrg);
3301    numnw += it;
3302    if (writeoptn)
3303      printsmooth(tr, it);
3304  } while (!(numsm >= maxsmooth || cnvrg));
3305  tr->cnvrgnc = cnvrg;
3306
3307  /* IF debugoptn THEN printcurtree ( tr ); */
3308  evaluateA(tr);
3309  variance(tr);
3310  if (!usertree)
3311    return;
3312  lklhdtree[notree] = tr->lklhd;
3313  aictree[notree] = tr->aic;
3314  paratree[notree] = ibrnch2;
3315  if (notree == 1) {
3316    maxltr = 1;
3317    maxlkl = tr->lklhd;
3318    minatr = 1;
3319    minaic = tr->aic;
3320    return;
3321  }
3322  if (tr->lklhd > maxlkl) {
3323    maxltr = notree;
3324    maxlkl = tr->lklhd;
3325  }
3326  if (tr->aic < minaic) {
3327    minatr = notree;
3328    minaic = tr->aic;
3329  }
3330}  /* manuvaluate */
3331
3332
3333Static Void autovaluate(tr)
3334tree *tr;
3335{
3336  long n, it;
3337  lengthty newarcs, oldarcs;
3338  boolean cnvrg;
3339  long FORLIM;
3340
3341  if (writeoptn)
3342    putchar('\n');
3343  if (debugoptn)
3344    printf(" AUTOVALUATE\n");
3345  /*
3346     IF debugoptn THEN printcurtree( tr );
3347     IF debugoptn THEN writeln(output);
3348     initdescen( tr );
3349     initbranch (tr.startp);
3350     initbranch (tr.startp^.kinp);
3351     IF debugoptn THEN printcurtree ( tr );
3352     initlength ( tr );
3353     IF debugoptn THEN printcurtree ( tr );
3354  */
3355  it = 0;
3356  for (n = 1; n <= 3; n++) {
3357    sublklhdA(tr->brnchp[ibrnch2]);
3358    sublklhdA(tr->brnchp[ibrnch2]->kinp);
3359    branchlengthA(tr->brnchp[ibrnch2], &it);
3360  }
3361  if (debugoptn)
3362    printcurtree(tr);
3363
3364  if (debugoptn)
3365    printf(" SMOOTH\n");
3366  cnvrg = false;
3367  FORLIM = ibrnch2;
3368  for (n = 1; n <= FORLIM; n++)
3369    newarcs[n - 1] = tr->brnchp[n]->UU.U0.lnga;
3370  numnw = 0;
3371  numsm = 0;
3372  do {
3373    numsm++;
3374    it = 0;
3375    smooth(tr, tr->startp->kinp, &it);
3376    memcpy(oldarcs, newarcs, sizeof(lengthty));
3377    FORLIM = ibrnch2;
3378    for (n = 1; n <= FORLIM; n++)
3379      newarcs[n - 1] = tr->brnchp[n]->UU.U0.lnga;
3380    judgconverg(oldarcs, newarcs, &cnvrg);
3381    numnw += it;
3382    if (writeoptn)
3383      printsmooth(tr, it);
3384  } while (!(numsm >= maxsmooth || cnvrg));
3385  tr->cnvrgnc = cnvrg;
3386
3387  /* IF debugoptn THEN printcurtree ( tr ); */
3388  evaluateA(tr);
3389  variance(tr);
3390  if (usertree)
3391    return;
3392  lklhdtree[notree] = tr->lklhd;
3393  aictree[notree] = tr->aic;
3394  paratree[notree] = ibrnch2;
3395  if (notree == 1) {
3396    maxltr = 1;
3397    maxlkl = tr->lklhd;
3398    minatr = 1;
3399    minaic = tr->aic;
3400    return;
3401  }
3402  if (tr->lklhd > maxlkl) {
3403    maxltr = notree;
3404    maxlkl = tr->lklhd;
3405  }
3406  if (tr->aic < minaic) {
3407    minatr = notree;
3408    minaic = tr->aic;
3409  }
3410}  /* autovaluate */
3411
3412
3413#define maxover         50
3414#define maxleng         30
3415
3416
3417/**************************************/
3418/***    OUTPUT OF TREE TOPOLOGY     ***/
3419/**************************************/
3420
3421Static Void prbranch(up, depth, m, maxm, umbrella, length)
3422node *up;
3423long depth, m, maxm;
3424boolean *umbrella;
3425long *length;
3426{
3427  long i, j, n, d, maxn;
3428  node *cp;
3429  boolean over;
3430  long FORLIM, FORLIM1;
3431
3432  over = false;
3433  d = depth + 1;
3434  if ((long)up->UU.U0.lnga >= maxover) {   /* +4 */
3435    over = true;
3436    length[d] = maxleng;
3437  } else
3438    length[d] = (long)up->UU.U0.lnga + 3;
3439  if (up->isop == NULL) {   /* TIP */
3440    if (m == 1)
3441      umbrella[d - 1] = true;
3442    for (j = 0; j < d; j++) {
3443      if (umbrella[j])
3444        printf("%*c", (int)length[j], ':');
3445      else
3446        printf("%*c", (int)length[j], ' ');
3447    }
3448    if (m == maxm)
3449      umbrella[d - 1] = false;
3450    if (over) {
3451      FORLIM = length[d] - 3;
3452      for (i = 1; i <= FORLIM; i++)
3453        putchar('+');
3454    } else {
3455      FORLIM = length[d] - 3;
3456      for (i = 1; i <= FORLIM; i++)
3457        putchar('-');
3458    }
3459    if (up->number < 10)
3460      printf("--%ld.", up->number);
3461    else if (up->number < 100)
3462      printf("-%2ld.", up->number);
3463    else
3464      printf("%3ld.", up->number);
3465    for (i = 0; i < maxname; i++)
3466      putchar(up->namesp[i]);
3467    /* write(output,d:36); */
3468    putchar('\n');
3469    return;
3470  }
3471  cp = up->isop;
3472  maxn = up->diverg;
3473  for (n = 1; n <= maxn; n++) {
3474    prbranch(cp->kinp, d, n, maxn, umbrella, length);
3475    cp = cp->isop;
3476    if (n == maxn / 2) {
3477      if (m == 1)
3478        umbrella[d - 1] = true;
3479      if (n == 1)
3480        umbrella[d - 1] = true;
3481      for (j = 0; j < d; j++) {
3482        if (umbrella[j])
3483          printf("%*c", (int)length[j], ':');
3484        else
3485          printf("%*c", (int)length[j], ' ');
3486      }
3487      if (n == maxn)
3488        umbrella[d - 1] = false;
3489      if (m == maxm)
3490        umbrella[d - 1] = false;
3491      if (over) {
3492        FORLIM1 = length[d] - 3;
3493        for (i = 1; i <= FORLIM1; i++)
3494          putchar('+');
3495      } else {
3496        FORLIM1 = length[d] - 3;
3497        for (i = 1; i <= FORLIM1; i++)
3498          putchar('-');
3499      }
3500      if (up->number < 10)   /* ,':' */
3501        printf("--%ld", up->number);
3502      else if (up->number < 100)
3503        printf("-%2ld", up->number);
3504      else
3505        printf("%3ld", up->number);
3506      /* write(output,d:50); */
3507      putchar('\n');
3508    } else if (n != maxn) {
3509      for (j = 0; j <= d; j++) {
3510        if (umbrella[j])
3511          printf("%*c", (int)length[j], ':');
3512        else
3513          printf("%*c", (int)length[j], ' ');
3514      }
3515      putchar('\n');
3516    }
3517  }
3518
3519  /* ,':' */
3520  /* ,':' */
3521}  /* prbranch */
3522
3523#undef maxover
3524#undef maxleng
3525
3526
3527Static Void prtopology(tr)
3528tree *tr;
3529{
3530  long n, maxn, depth;
3531  node *cp;
3532  umbty umbrella;
3533  lspty length;
3534
3535  for (n = 0; n <= maxsp; n++) {
3536    umbrella[n] = false;
3537    if (n == 0)
3538      length[n] = 3;
3539    else
3540      length[n] = 6;
3541  }
3542  cp = tr->brnchp[0];
3543  maxn = cp->diverg + 1;
3544  putchar('\n');
3545  for (n = 1; n <= maxn; n++) {
3546    depth = 0;
3547    prbranch(cp->kinp, depth, n, maxn, umbrella, length);
3548    cp = cp->isop;
3549    if (n == maxn / 2)
3550      printf("%*s\n", (int)length[0], "0:");
3551    else if (n != maxn)
3552      printf("%*c\n", (int)length[0], ':');
3553  }
3554}  /* prtopology */
3555
3556
3557/*********************************/
3558/***  OUTPUT OF TREE TOPOLOGY  ***/
3559/*********************************/
3560
3561Static Void charasubtoporogy(up, nsp, nch)
3562node *up;
3563long *nsp, *nch;
3564{
3565  long n, maxn;
3566  node *cp;
3567
3568  if (up->isop == NULL) {   /* TIP */
3569    for (n = 0; n < maxname; n++) {
3570      if (up->namesp[n] != ' ') {
3571        putchar(up->namesp[n]);
3572        (*nch)++;
3573      }
3574    }
3575    (*nsp)++;
3576    return;
3577  }
3578  cp = up->isop;
3579  maxn = up->diverg;
3580  putchar('(');
3581  (*nch)++;
3582  for (n = 1; n <= maxn; n++) {
3583    charasubtoporogy(cp->kinp, nsp, nch);
3584    cp = cp->isop;
3585    if (n != maxn) {
3586      putchar(',');
3587      (*nch)++;
3588      if (*nch > maxline - 10) {
3589        printf("\n%3c", ' ');
3590        *nch = 3;
3591      }
3592    }
3593  }
3594  putchar(')');
3595  (*nch)++;
3596}  /* charasubtoporogy */
3597
3598
3599Static Void charatopology(tr)
3600tree *tr;
3601{
3602  long n, maxn, nsp, nch;
3603  node *cp;
3604
3605  nsp = 0;
3606  nch = 3;
3607  cp = tr->brnchp[0];
3608  maxn = cp->diverg + 1;
3609  printf("\n%3s", " ( ");
3610  for (n = 1; n <= maxn; n++) {
3611    charasubtoporogy(cp->kinp, &nsp, &nch);
3612    cp = cp->isop;
3613    if (n != maxn) {
3614      printf("%2s", ", ");
3615      nch += 2;
3616      if (nch > maxline - 15) {
3617        printf("\n%3c", ' ');
3618        nch = 3;
3619      }
3620    }
3621  }
3622  printf(" )\n");
3623}  /* charatopology */
3624
3625
3626/**************************************/
3627/***   OUTPUT OF THE FINAL RESULT   ***/
3628/**************************************/
3629
3630Static Void printbranch(tr)
3631tree *tr;
3632{
3633  long i, j;
3634  node *p, *q;
3635  long FORLIM;
3636
3637  FORLIM = ibrnch2;
3638  for (i = 0; i < FORLIM; i++) {
3639    p = tr->brnchp[i + 1];
3640    q = p->kinp;
3641    printf("%5c", ' ');
3642    if (p->isop == NULL) {   /* TIP */
3643      for (j = 0; j < maxname; j++)
3644        putchar(p->namesp[j]);
3645    } else {
3646      for (j = 1; j <= maxname; j++)
3647        putchar(' ');
3648    }
3649    printf("%5ld", p->number);
3650    if (p->UU.U0.lnga >= maxarc)
3651      printf("%12s", " infinity");
3652    else if (p->UU.U0.lnga <= minarc)
3653      printf("%12s", " zero    ");
3654    else
3655      printf("%12.5f", p->UU.U0.lnga);
3656    printf("%3s%9.5f%2s", " (", sqrt(fabs(tr->vrilnga[i])), " )");
3657    if (debugoptn) {
3658      printf("%12.5f", sqrt(fabs(tr->vrilnga2[i])));
3659      printf("%13.5f", tr->blklhd[i]);
3660    }
3661    putchar('\n');
3662  }
3663
3664}  /* printbranch */
3665
3666
3667Static Void summarize(tr)
3668tree *tr;
3669{
3670  putchar('\n');
3671  /* printline( 46 ); */
3672  if (usertree)
3673    printf("%4s%3ld%8c", " No.", notree, ' ');
3674  else
3675    printf("%4s%3ld%2s%3ld%3c", " No.", stage, " -", notree, ' ');
3676  printf("%7s%9s%11s", "number", "Length", "S.E.");
3677  if (!tr->cnvrgnc)
3678    printf("%21s", "non convergence");
3679  /* write  (output, 'convergence    ':21) */
3680  if (writeoptn)
3681    printf("%2c%3ld%2c%4ld", ':', numsm, ',', numnw);
3682  putchar('\n');
3683  printline(46L);
3684  printbranch(tr);
3685  printline(46L);
3686  printf("%8s", "ln L :");
3687  printf("%10.3f", tr->lklhd);
3688  printf("%2c%8.3f%2s", '(', sqrt(fabs(tr->vrilkl)), " )");
3689  printf("%7s%9.3f\n", "AIC :", tr->aic);
3690  printline(46L);
3691  /*       writeln(output); */
3692}  /* summarize */
3693
3694
3695Static Void cleartree(tr)
3696tree *tr;
3697{
3698  long i, n;
3699  node *cp, *sp, *dp;
3700  long FORLIM;
3701
3702  FORLIM = ibrnch2;
3703  for (i = 1; i <= FORLIM; i++) {
3704    tr->brnchp[i]->kinp->kinp = NULL;
3705    tr->brnchp[i]->kinp = NULL;
3706  }
3707  FORLIM = ibrnch2 + 1;
3708  for (i = ibrnch1; i <= FORLIM; i++) {
3709    if (i == ibrnch2 + 1)
3710      n = 0;
3711    else
3712      n = i;
3713    sp = tr->brnchp[n];
3714    cp = sp->isop;
3715    sp->isop = NULL;
3716    while (cp != sp) {
3717      dp = cp;
3718      cp = cp->isop;
3719      dp->isop = freenode;
3720      freenode = dp;
3721      dp = NULL;
3722      /* dp^.isop := NIL;
3723         DISPOSE( dp ); */
3724    }
3725    sp = NULL;
3726    tr->brnchp[n] = NULL;
3727    cp->isop = freenode;
3728    freenode = cp;
3729    cp = NULL;
3730    /* DISPOSE( cp ); */
3731  }
3732}  /* cleartree */
3733
3734
3735Static Void bootstrap()
3736{
3737  long ib, jb, nb, bsite, maxboottree, inseed;
3738  double maxlklboot;
3739  longintty seed;
3740  long FORLIM, FORLIM1, FORLIM2;
3741
3742  inseed = 12345;
3743  for (ib = 0; ib <= 5; ib++)
3744    seed[ib] = 0;
3745  ib = 0;
3746  do {
3747    seed[ib] = inseed & 63;
3748    inseed /= 64;
3749    ib++;
3750  } while (inseed != 0);
3751
3752  FORLIM = numtrees;
3753  for (jb = 1; jb <= FORLIM; jb++)
3754    boottree[jb] = 0.0;
3755  for (nb = 1; nb <= maxboot; nb++) {
3756    FORLIM1 = numtrees;
3757    for (jb = 1; jb <= FORLIM1; jb++)
3758      lklboottree[jb] = 0.0;
3759    FORLIM1 = endsite;
3760    for (ib = 1; ib <= FORLIM1; ib++) {
3761      bsite = weightw[(int)((long)(randum(seed) * endsite + 1)) - 1];
3762      FORLIM2 = numtrees;
3763      for (jb = 1; jb <= FORLIM2; jb++)
3764        lklboottree[jb] += lklhdtrpt[jb][bsite - 1];
3765    }
3766
3767    maxlklboot = lklboottree[1];
3768    maxboottree = 1;
3769    if (debugoptn)
3770      printf("%5ld", nb);
3771    FORLIM1 = numtrees;
3772    for (jb = 1; jb <= FORLIM1; jb++) {
3773      if (lklboottree[jb] > maxlklboot) {
3774        maxlklboot = lklboottree[jb];
3775        maxboottree = jb;
3776      }
3777      if (debugoptn)
3778        printf("%8.1f", lklboottree[jb]);
3779    }
3780    if (debugoptn)
3781      printf("%4ld\n", maxboottree);
3782    boottree[maxboottree] += 1.0;
3783  }
3784  if (maxboot == 0) {
3785    FORLIM = numtrees;
3786    for (jb = 1; jb <= FORLIM; jb++)
3787      boottree[jb] /= maxboot;
3788  }
3789}  /* bootstrap */
3790
3791
3792Static Void printlklhd()
3793{
3794  long i, j, cul;
3795  double suml;   /* difference of ln lklhd */
3796  double suml2;   /* absolute value of difference */
3797  double suma;   /* difference of AIC */
3798  double suma2;   /* absolute value of AIC */
3799  double lklsite;   /* likelihood of site */
3800  double aicsite;   /* AIC of site */
3801  double sdl;   /* standard error ln lklhd */
3802  double sda;   /* standard error of AIC */
3803  long FORLIM, FORLIM1;
3804  double TEMP;
3805
3806  cul = 71;
3807  putchar('\n');
3808  if (usertree)
3809    printf("%10c%4ld user trees\n", ' ', numtrees);
3810  else
3811    printf("%10s%3ld%6s%5ld%6s\n", " NO.", stage, "stage", notree, "trees");
3812  if (numtrees <= 0 || endsite <= 1)
3813    return;
3814  printline(cul);
3815  printf("%6s%6s%12s%12s%6s%12s%17s\n",
3816         " Tree", "ln L", "Diff ln L", "#Para", "AIC", "Diff AIC", "Boot P");
3817  printline(cul);
3818  FORLIM = numtrees;
3819  for (i = 1; i <= FORLIM; i++) {
3820    printf("%4ld%9.2f", i, lklhdtree[i]);
3821    if (maxltr == i)
3822      printf(" <-- best%9c", ' ');
3823    else {
3824      suml = 0.0;
3825      suml2 = 0.0;
3826      FORLIM1 = endptrn;
3827      for (j = 0; j < FORLIM1; j++) {
3828        lklsite = lklhdtrpt[maxltr][j] - lklhdtrpt[i][j];
3829        suml += weight[j] * lklsite;
3830        suml2 += weight[j] * lklsite * lklsite;
3831      }
3832      TEMP = suml / endsite;
3833      sdl = sqrt(endsite / (endsite - 1.0) * (suml2 - TEMP * TEMP));
3834      printf("%9.2f", lklhdtree[i] - maxlkl);
3835      printf("%3s%6.2f", "+-", sdl);
3836    }
3837    printf("%4ld", paratree[i]);
3838    printf("%9.2f", aictree[i]);
3839    if (minatr == i)
3840      printf(" <-- best%9c", ' ');
3841    else {
3842      suma = 0.0;
3843      suma2 = 0.0;
3844      FORLIM1 = endptrn;
3845      for (j = 0; j < FORLIM1; j++) {
3846        aicsite = -2.0 * (lklhdtrpt[maxltr][j] - lklhdtrpt[i][j]);
3847        suma += weight[j] * aicsite;
3848        suma2 += weight[j] * aicsite * aicsite;
3849      }
3850      TEMP = suma / endsite;
3851      sda = sqrt(endsite / (endsite - 1.0) * (suma2 - TEMP * TEMP));
3852      printf("%9.2f", aictree[i] - minaic);
3853      printf("%3s%6.2f", "+-", sda);
3854    }
3855    if (usertree)
3856      printf("%9.4f", boottree[i]);
3857    putchar('\n');
3858  }
3859  printline(cul);
3860}  /* printlklhd */
3861
3862
3863Static Void outlklhd()
3864{
3865  long i, j, k, nt, FORLIM, FORLIM1, FORLIM2;
3866
3867  fprintf(lklfile, "%5ld%5ld\n", numsp, endsite);
3868  i = 0;
3869  FORLIM = numtrees;
3870  for (nt = 1; nt <= FORLIM; nt++) {
3871    fprintf(lklfile, "%5ld\n", nt);
3872    FORLIM1 = endptrn;
3873    for (j = 0; j < FORLIM1; j++) {
3874      FORLIM2 = weight[j];
3875      for (k = 1; k <= FORLIM2; k++) {
3876        fprintf(lklfile, "% .5E", lklhdtrpt[nt][j]);
3877        i++;
3878        if (i == 6) {
3879          putc('\n', lklfile);
3880          i = 0;
3881        }
3882      }
3883    }
3884    if (i != 0)
3885      putc('\n', lklfile);
3886    i = 0;
3887  }
3888}  /* outlklhd */
3889
3890
3891Static Void manutree()
3892{
3893  long ntree, FORLIM;
3894
3895  /*PAGE(output);*/
3896  fscanf(seqfile, "%ld%*[^\n]", &numtrees);
3897  getc(seqfile);
3898  printf("\n %5ld user trees\n\n", numtrees);
3899  FORLIM = numtrees;
3900  for (ntree = 1; ntree <= FORLIM; ntree++) {
3901    notree = ntree;
3902    treestructure(&ctree);
3903    manuvaluate(&ctree);
3904    prtopology(&ctree);
3905    summarize(&ctree);
3906    cleartree(&ctree);
3907    /*PAGE(output);*/
3908  }
3909  if (bootsoptn)
3910    bootstrap();
3911  printlklhd();
3912  if (putlkoptn)
3913    outlklhd();
3914}  /* manutree */
3915
3916
3917Static Void newbranch(nbranch, np1, np2)
3918long nbranch;
3919node **np1, **np2;
3920{
3921  long j;
3922  node *np;
3923
3924  for (j = 1; j <= 2; j++) {
3925    if (freenode == NULL)
3926      np = (node *)Malloc(sizeof(node));
3927    else {
3928      np = freenode;
3929      freenode = np->isop;
3930      np->isop = NULL;
3931    }
3932    clearnode(np);
3933    if (j == 1)
3934      *np1 = np;
3935    else
3936      *np2 = np;
3937    np = NULL;
3938  }
3939  (*np1)->kinp = *np2;
3940  (*np2)->kinp = *np1;
3941  ctree.brnchp[nbranch] = *np1;
3942  ctree.brnchp[nbranch]->number = ibrnch2;
3943}  /* newbranch */
3944
3945
3946Static Void decomposition(cnode)
3947long cnode;
3948{
3949  long i1, i2, maxdvg;
3950  node *cp1, *cp2, *bp1, *bp2, *lp1, *lp2, *np1, *np2;
3951  double maxlkls;
3952
3953  if (ctree.brnchp[cnode]->diverg <= 2) {
3954    return;
3955  }  /* IF */
3956  /*
3957     PAGE(output);
3958  */
3959  stage++;
3960  notree = 0;
3961  ibrnch2++;
3962  newbranch(ibrnch2, &np1, &np2);
3963  cp1 = ctree.brnchp[cnode];
3964  cp2 = cp1;
3965  bp1 = cp1;
3966  do {
3967    bp1 = bp1->isop;
3968  } while (bp1->isop != cp1);
3969  bp2 = bp1;
3970  maxlkls = -999999.0;
3971  maxdvg = ctree.brnchp[cnode]->diverg;
3972  if (maxdvg == 3) {
3973    maxdvg--;
3974    cp1 = cp1->isop;
3975    cp2 = cp1;
3976    bp1 = bp1->isop;
3977    bp2 = bp1;
3978  }
3979  for (i1 = 1; i1 <= maxdvg; i1++) {
3980    for (i2 = i1 + 1; i2 <= maxdvg + 1; i2++) {
3981      if (debugoptn) {
3982        ibrnch2--;
3983        printcurtree(&ctree);
3984        ibrnch2++;
3985      }
3986      notree++;
3987      bp2 = cp2;
3988      cp2 = cp2->isop;
3989      if (debugoptn)
3990        printf(" AUTO-INS%3ld%3ld%4c%3ld%3ld%4c%3ld%3ld%4c%3ld%3ld\n",
3991               i1, i2, 'c', cp1->kinp->number, cp2->kinp->number, 'b',
3992               bp1->kinp->number, bp2->kinp->number, 'n', np1->number,
3993               np2->number);
3994      insertbranch(&ctree, cp1, cp2, bp1, bp2, np1, np2);
3995      autovaluate(&ctree);
3996      prtopology(&ctree);
3997      charatopology(&ctree);
3998      summarize(&ctree);
3999      if (ctree.lklhd > maxlkls) {
4000        maxlkls = ctree.lklhd;
4001        lp1 = bp1;
4002        lp2 = bp2;
4003      }
4004      if (debugoptn)
4005        printf(" AUTO-DEL%3ld%3ld%4c%3ld%3ld%4c%3ld%3ld%4c%3ld%3ld\n",
4006               i1, i2, 'c', cp1->kinp->number, cp2->kinp->number, 'b',
4007               bp1->kinp->number, bp2->kinp->number, 'n', np1->number,
4008               np2->number);
4009      deletebranch(&ctree, cnode, cp1, cp2, bp1, bp2, np1, np2);
4010    }
4011    bp1 = cp1;
4012    cp1 = bp1->isop;
4013    bp2 = bp1;
4014    cp2 = bp1->isop;
4015  }  /* FOR */
4016  /*
4017     PAGE(output);
4018  */
4019  numtrees = notree;
4020  printlklhd();
4021  notree = 0;
4022  cp1 = lp1->isop;
4023  cp2 = lp2->isop;
4024  bp1 = lp1;
4025  bp2 = lp2;
4026  if (debugoptn)
4027    printf(" AUTO-MAX%4c%3ld%3ld%4c%3ld%3ld%4c%3ld%3ld\n",
4028           'c', cp1->kinp->number, cp2->kinp->number, 'b', bp1->kinp->number,
4029           bp2->kinp->number, 'n', np1->number, np2->number);
4030  insertbranch(&ctree, cp1, cp2, bp1, bp2, np1, np2);
4031  autovaluate(&ctree);
4032  prtopology(&ctree);
4033  summarize(&ctree);
4034}  /* decomposition */
4035
4036
4037Static Void autotree()
4038{
4039  long i, j, dvg, cnode, FORLIM;
4040
4041  stage = 0;
4042  notree = 1;
4043  if (semiaoptn)
4044    treestructure(&ctree);
4045  else
4046    starstructure(&ctree);
4047  manuvaluate(&ctree);
4048  prtopology(&ctree);
4049  summarize(&ctree);
4050
4051  if (!semiaoptn) {
4052    do {
4053      decomposition(0L);
4054    } while (ctree.brnchp[0]->diverg >= 3);
4055    return;
4056  }
4057
4058  /* printlklhd; */
4059  /* IF putlkoptn THEN outlklhd; */
4060  if (firstoptn) {
4061    do {
4062      decomposition(0L);
4063    } while (ctree.brnchp[0]->diverg >= 3);
4064    do {
4065      dvg = 2;   /* ctree.brnchp[0]^.diverg */
4066      cnode = 0;
4067      FORLIM = ibrnch2;
4068      for (i = ibrnch1; i <= FORLIM; i++) {
4069        if (ctree.brnchp[i]->diverg > dvg) {
4070          dvg = ctree.brnchp[i]->diverg;
4071          cnode = i;
4072        }
4073      }
4074      decomposition(cnode);
4075    } while (dvg >= 3);
4076    return;
4077  }
4078  if (lastoptn) {
4079    do {
4080      dvg = 2;   /* ctree.brnchp[0]^.diverg */
4081      cnode = 0;
4082      FORLIM = ibrnch2;
4083      for (i = ibrnch1; i <= FORLIM; i++) {
4084        if (ctree.brnchp[i]->diverg > dvg) {
4085          dvg = ctree.brnchp[i]->diverg;
4086          cnode = i;
4087        }
4088      }
4089      decomposition(cnode);
4090    } while (dvg >= 3);
4091    do {
4092      decomposition(0L);
4093    } while (ctree.brnchp[0]->diverg >= 3);
4094    return;
4095  }
4096  do {
4097    dvg = 2;   /* ctree.brnchp[0]^.diverg */
4098    cnode = 0;
4099    FORLIM = ibrnch2 + 1;
4100    for (i = ibrnch1; i <= FORLIM; i++) {
4101      if (i == ibrnch2 + 1)
4102        j = 0;
4103      else
4104        j = i;
4105      if (ctree.brnchp[j]->diverg > dvg) {
4106        dvg = ctree.brnchp[j]->diverg;
4107        cnode = j;
4108      }
4109    }
4110    decomposition(cnode);
4111  } while (dvg >= 3);
4112
4113  /* auto mode */
4114}  /* autotree */
4115
4116
4117Static Void mainio()
4118{
4119  long nexe;
4120
4121  freenode = NULL;
4122  for (nexe = 1; nexe <= maxexe; nexe++) {
4123    numexe = nexe;
4124    if (maxexe > 1)
4125      printf(" * JOB :%4ld *\n", numexe);
4126    getinput();
4127    if (normal && numexe == 1)
4128      maketransprob();
4129    if (normal) {
4130      if (usertree)
4131        manutree();
4132      else
4133        autotree();
4134    }
4135    /* IF maxexe > 1 THEN PAGE(output); */
4136  }
4137}  /* mainio */
4138
4139
4140main(argc, argv)
4141int argc;
4142Char *argv[];
4143{
4144  /* ASSIGN(seqfile,seqfname);      If TURBO Pascal, use */
4145  /* ASSIGN(tpmfile,tpmfname); */
4146  /* IF putlkoptn THEN
4147     ASSIGN (lklfile,lklfname); */
4148
4149  PASCAL_MAIN(argc, argv);
4150  tpmfile = NULL;
4151  lklfile = NULL;
4152  seqfile = NULL;
4153  rewind(seqfile);
4154  /* If SUN Pascal, don't use */
4155  /* RESET (tpmfile); */
4156  /* IF putlkoptn THEN
4157     REWRITE(lklfile); */
4158
4159  /* RESET (seqfile,seqfname);      If SUN Pascal, use */
4160  /* RESET (tpmfile,tpmfname); */
4161  /* IF putlkoptn THEN
4162     REWRITE(lklfile,lklfname); */
4163
4164  mainio();
4165
4166  /* CLOSE (seqfile);               If TURBO Pascal, use */
4167  /* CLOSE (tpmfile); */
4168  /* IF putlkoptn THEN
4169     CLOSE (lklfile); */
4170
4171  if (seqfile != NULL)
4172    fclose(seqfile);
4173  if (lklfile != NULL)
4174    fclose(lklfile);
4175  if (tpmfile != NULL)
4176    fclose(tpmfile);
4177  exit(EXIT_SUCCESS);
4178}  /* PROTML */
4179
4180
4181
4182
4183/* End. */
Note: See TracBrowser for help on using the repository browser.