forester/archive/RIO/others/phylip_mod/src/phylip.c

   1
   2 /* version 3.6. (c) Copyright 1993-2004 by the University of Washington.
   3    Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, Andrew Keeffe,
   4    and Dan Fineman.
   5    Permission is granted to copy and use this program provided no fee is
   6    charged for it and provided that this copyright notice is not removed. */
   7
   8 #ifdef OSX_CARBON
   9 #include <Carbon/Carbon.h>
  10 #endif
  11
  12 #include <stdio.h>
  13 #include <signal.h>
  14 #ifdef WIN32
  15 #include <windows.h>
  16 /* for console code (clear screen, text color settings) */
  17 CONSOLE_SCREEN_BUFFER_INFO savecsbi;
  18 HANDLE hConsoleOutput;
  19
  20 void phyClearScreen();
  21 void phySaveConsoleAttributes();
  22 void phySetConsoleAttributes();
  23 void phyRestoreConsoleAttributes();
  24 void phyFillScreenColor();
  25 #endif
  26
  27 #include "phylip.h"
  28
  29 #ifndef OLDC
  30 static void crash_handler(int signum);
  31
  32 #endif
  33 #if defined(OSX_CARBON) && defined(__MWERKS__)
  34 boolean fixedpath = false;
  35 #endif
  36 FILE *infile, *outfile, *intree, *intree2, *outtree, *weightfile, *catfile, *ancfile, *mixfile, *factfile;
  37 long spp, words, bits;
  38 boolean ibmpc, ansi, tranvsp;
  39 naym *nayme;                     /* names of species */
  40
  41 static void crash_handler(int sig_num)
  42 { /* when we crash, lets print out something usefull */
  43   printf("ERROR:  ");
  44   switch(sig_num) {
  45 #ifdef SIGSEGV
  46     case SIGSEGV:
  47       puts("This program has caused a Segmentation fault.");
  48       break;
  49 #endif /* SIGSEGV */
  50 #ifdef SIGFPE
  51     case SIGFPE:
  52       puts("This program has caused a Floating Point Exception");
  53       break;
  54 #endif  /* SIGFPE */
  55 #ifdef SIGILL
  56     case SIGILL:
  57       puts("This program has attempted an illegal instruction");
  58       break;
  59 #endif  /* SIGILL */
  60 #ifdef SIGPIPE
  61     case SIGPIPE:
  62       puts("This program tried to write to a broken pipe");
  63       break;
  64 #endif  /* SIGPIPE */
  65 #ifdef SIGBUS
  66     case SIGBUS:
  67       puts("This program had a bus error");
  68       break;
  69 #endif /* SIGBUS */
  70   }
  71   if (sig_num == SIGSEGV) {
  72     puts(
  73     "       This may have been caused by an incorrectly formatted input file");
  74     puts(
  75         "       or input tree file.  You should check those files carefully.");
  76     puts("       If this seems to be a bug, please mail joe@gs.washington.edu");
  77   }
  78   else {
  79     puts("       Most likely, you have encountered a bug in the program.");
  80     puts("       Since this seems to be a bug, please mail joe@gs.washington.edu");
  81   }
  82   puts("       with the name of the program, your computer system type,");
  83   puts("       a full description of the problem, and with the input data file.");
  84   puts("       (which should be in the body of the message, not as an Attachment).");
  85
  86 #ifdef WIN32
  87   puts ("Press Enter or Return to close program.");
  88   puts("  You may have to press Enter or Return twice.");
  89   getchar ();
  90   getchar ();
  91   phyRestoreConsoleAttributes();
  92 #endif
  93   abort();
  94 }
  95
  96
  97 void init(int argc, char** argv)
  98 { /* initialization routine for all programs
  99    * anything done at the beginig for every program should be done here */
 100
 101   /* set up signal handler for
 102    * segfault,floating point exception, illeagal instruction, bad pipe, bus error
 103    * there are more signals that can cause a crash, but these are the most common
 104    * even these aren't found on all machines.  */
 105 #ifdef SIGSEGV
 106   signal(SIGSEGV, crash_handler);
 107 #endif /* SIGSEGV */
 108 #ifdef SIGFPE
 109   signal(SIGFPE, crash_handler);
 110 #endif /* SIGFPE */
 111 #ifdef SIGILL
 112   signal(SIGILL, crash_handler);
 113 #endif /* SIGILL */
 114 #ifdef SIGPIPE
 115   signal(SIGPIPE, crash_handler);
 116 #endif /* SIGPIPE */
 117 #ifdef SIGBUS
 118   signal(SIGBUS, crash_handler);
 119 #endif /* SIGBUS */
 120
 121 #ifdef WIN32
 122   phySetConsoleAttributes();
 123   phyClearScreen();
 124 #endif
 125
 126 }
 127
 128 void scan_eoln(FILE *f)
 129 { /* eat everything to the end of line or eof*/
 130   char ch;
 131
 132   while (!eoff(f) && !eoln(f))
 133     gettc(f);
 134   if (!eoff(f))
 135     ch = gettc(f);
 136 }
 137
 138
 139 boolean eoff(FILE *f)
 140 { /* check for end of file */
 141     int ch;
 142
 143     if (feof(f))
 144       return true;
 145     ch = getc(f);
 146     if (ch == EOF) {
 147       ungetc(ch, f);
 148       return true;
 149     }
 150     ungetc(ch, f);
 151     return false;
 152 }  /*eoff*/
 153
 154
 155 boolean eoln(FILE *f)
 156 { /* check for end of line or eof*/
 157     register int ch;
 158
 159     ch = getc(f);
 160     if (ch == EOF)
 161       return true;
 162     ungetc(ch, f);
 163     return ((ch == '\n') || (ch == '\r'));
 164 }  /*eoln*/
 165
 166
 167 int filexists(char *filename)
 168 { /* check whether file already exists */
 169   FILE *fp;
 170   fp =fopen(filename,"rb");
 171   if (fp) {
 172     fclose(fp);
 173     return 1;
 174   } else
 175     return 0;
 176 }  /*filexists*/
 177
 178
 179 const char* get_command_name (const char *vektor)
 180 { /* returns the name of the program from vektor without the whole path */
 181   char *last_slash;
 182
 183   /* Point to the last slash... */
 184   last_slash = strrchr (vektor, DELIMITER);
 185
 186   if (last_slash)
 187     /* If there was a last slash, return the character after it */
 188     return last_slash + 1;
 189   else
 190     /* If not, return the vector */
 191     return vektor;
 192
 193 }  /*get_command_name*/
 194
 195
 196 void getstryng(char *fname)
 197 { /* read in a file name from stdin and take off newline if any */
 198
 199   fname = fgets(fname, 100, stdin);
 200   if (strchr(fname, '\n') != NULL)
 201     *strchr(fname, '\n') = '\0';
 202 } /* getstryng */
 203
 204
 205 void countup(long *loopcount, long maxcount)
 206 { /* count how many times this loop has tried to read data, bail out
 207      if exceeds maxcount */
 208
 209   (*loopcount)++;
 210   if ((*loopcount) >= maxcount) {
 211     printf("\nERROR: Made %ld attempts to read input in loop. Aborting run.\n",
 212             *loopcount);
 213     exxit(-1);
 214   }
 215 } /* countup */
 216
 217
 218 void openfile(FILE **fp,const char *filename,const char *filedesc,
 219               const char *mode,const char *application, char *perm)
 220 { /* open a file, testing whether it exists etc. */
 221   FILE *of;
 222   char file[FNMLNGTH];
 223   char filemode[3];
 224   char input[FNMLNGTH];
 225   char ch;
 226   const char *progname_without_path;
 227   long loopcount, loopcount2;
 228 #if defined(OSX_CARBON) && defined(__MWERKS__)
 229   ProcessSerialNumber myProcess;
 230   FSRef bundleLocation;
 231   unsigned char bundlePath[FNMLNGTH];
 232
 233   if(!fixedpath){
 234     /* change path to the bundle location instead of root directory */
 235     GetCurrentProcess(&myProcess);
 236     GetProcessBundleLocation(&myProcess, &bundleLocation);
 237     FSRefMakePath(&bundleLocation, bundlePath, FNMLNGTH);
 238     chdir((const char*)bundlePath);
 239     chdir(".."); /* get out of the .app directory */
 240
 241     fixedpath = true;
 242   }
 243 #endif
 244
 245   progname_without_path = get_command_name(application);
 246
 247   strcpy(file,filename);
 248   strcpy(filemode,mode);
 249   strcat(filemode,"b");
 250   loopcount = 0;
 251   while (1){
 252     if (filemode[0] == 'w' && filexists(file)){
 253       printf("\n%s: the file \"%s\" that you wanted to\n",
 254           progname_without_path, file);
 255       printf("     use as %s already exists.\n", filedesc);
 256       printf("     Do you want to Replace it, Append to it,\n");
 257       printf("     write to a new File, or Quit?\n");
 258       loopcount2 = 0;
 259       do {
 260         printf("     (please type R, A, F, or Q) \n");
 261 #ifdef WIN32
 262         phyFillScreenColor();
 263 #endif
 264         fgets(input, sizeof(input), stdin);
 265         ch  = input[0];
 266         uppercase(&ch);
 267         countup(&loopcount2, 10);
 268       } while (ch != 'A' && ch != 'R' && ch != 'F' && ch != 'Q');
 269       if (ch == 'Q')
 270         exxit(-1);
 271       if (ch == 'A') {
 272         strcpy(filemode,"ab");
 273         continue;
 274       }
 275       else if (ch == 'F') {
 276         file[0] = '\0';
 277         loopcount2 = 0;
 278         while (file[0] =='\0') {
 279           printf("Please enter a new file name> ");
 280           getstryng(file);
 281           countup(&loopcount2, 10);
 282         }
 283         strcpy(filemode,"wb");
 284         continue;
 285       }
 286     }
 287     of = fopen(file,filemode);
 288     if (of)
 289       break;
 290     else {
 291       switch (filemode[0]){
 292
 293       case 'r':
 294         printf("%s: can't find %s \"%s\"\n", progname_without_path,
 295             filedesc, file);
 296         file[0] = '\0';
 297         loopcount2 = 0;
 298         while (file[0] =='\0'){
 299           printf("Please enter a new file name> ");
 300           countup(&loopcount2, 10);
 301           getstryng(file);}
 302         break;
 303
 304       case 'w':
 305       case 'a':
 306         printf("%s: can't write %s \"%s\"\n", progname_without_path,
 307             filedesc, file);
 308         file[0] = '\0';
 309         loopcount2 = 0;
 310         while (file[0] =='\0'){
 311           printf("Please enter a new file name> ");
 312           countup(&loopcount2, 10);
 313           getstryng(file);}
 314         continue;
 315       default:
 316      printf("There is some error in the call of openfile. Unknown mode.\n");
 317         exxit(-1);
 318       }
 319     }
 320     countup(&loopcount, 20);
 321   }
 322   *fp = of;
 323   if (perm != NULL)
 324     strcpy(perm,file);
 325 } /* openfile */
 326
 327
 328 void cleerhome()
 329 { /* home cursor and clear screen, if possible */
 330 #ifdef WIN32
 331   if(ibmpc || ansi){
 332     phyClearScreen();
 333   } else {
 334     printf("\n\n");
 335   }
 336 #else
 337   printf("%s", ((ibmpc || ansi) ? ("\033[2J\033[H") : "\n\n"));
 338 #endif
 339 } /* cleerhome */
 340
 341
 342 double randum(longer seed)
 343 { /* random number generator -- slow but machine independent
 344   This is a multiplicative congruential 32-bit generator
 345   x(t+1) = 1664525 * x(t) mod 2^32, one that passes the
 346   Coveyou-Macpherson and Lehmer tests, see Knuth ACP vol. 2
 347   We here implement it representing each integer in base-64
 348   notation -- i.e. as an array of 6 six-bit chunks   */
 349   long i, j, k, sum;
 350   longer mult, newseed;
 351   double x;
 352
 353   mult[0] = 13;   /* these four statements set the multiplier */
 354   mult[1] = 24;   /* -- they are its "digits" in a base-64    */
 355   mult[2] = 22;   /*    notation: 1664525 = 13*64^3+24*64^2   */
 356   mult[3] = 6;    /*                         +22*64+6         */
 357   for (i = 0; i <= 5; i++)
 358     newseed[i] = 0;
 359   for (i = 0; i <= 5; i++) {
 360     sum = newseed[i];
 361     k = i;
 362     if (i > 3)
 363       k = 3;
 364     for (j = 0; j <= k; j++)
 365       sum += mult[j] * seed[i - j];
 366     newseed[i] = sum;
 367     for (j = i; j <= 4; j++) {
 368       newseed[j + 1] += newseed[j] / 64;
 369       newseed[j] &= 63;
 370     }
 371   }
 372   memcpy(seed, newseed, sizeof(longer));
 373   seed[5] &= 3;
 374   x = 0.0;
 375   for (i = 0; i <= 5; i++)
 376     x = x / 64.0 + seed[i];
 377   x /= 4.0;
 378   return x;
 379 }  /* randum */
 380
 381
 382 void randumize(longer seed, long *enterorder)
 383 { /* randomize input order of species */
 384   long i, j, k;
 385
 386   for (i = 0; i < spp; i++) {
 387     j = (long)(randum(seed) * (i+1));
 388     k = enterorder[j];
 389     enterorder[j] = enterorder[i];
 390     enterorder[i] = k;
 391   }
 392 } /* randumize */
 393
 394
 395 double normrand(longer seed)
 396 {/* standardized Normal random variate */
 397   double x;
 398
 399   x = randum(seed)+randum(seed)+randum(seed)+randum(seed)
 400        + randum(seed)+randum(seed)+randum(seed)+randum(seed)
 401        + randum(seed)+randum(seed)+randum(seed)+randum(seed)-6.0;
 402   return(x);
 403 } /* normrand */
 404
 405
 406 long readlong(const char *prompt)
 407 { /* read a long */
 408   long res, loopcount;
 409   char string[100];
 410
 411   loopcount = 0;
 412   do {
 413     printf("%s",prompt);
 414     getstryng(string);
 415     if (sscanf(string,"%ld",&res) == 1)
 416       break;
 417     countup(&loopcount, 10);
 418    } while (1);
 419   return res;
 420 }  /* readlong */
 421
 422
 423 void uppercase(Char *ch)
 424 { /* convert ch to upper case */
 425   *ch = (islower (*ch) ? toupper(*ch) : (*ch));
 426 }  /* uppercase */
 427
 428
 429 void initseed(long *inseed, long *inseed0, longer seed)
 430 { /* input random number seed */
 431   long i, loopcount;
 432
 433   loopcount = 0;
 434   do {
 435     printf("Random number seed (must be odd)?\n");
 436     scanf("%ld%*[^\n]", inseed);
 437     getchar();
 438     countup(&loopcount, 10);
 439   } while (((*inseed) < 0) || ((*inseed) & 1) == 0);
 440   *inseed0 = *inseed;
 441   for (i = 0; i <= 5; i++)
 442     seed[i] = 0;
 443   i = 0;
 444   do {
 445     seed[i] = *inseed & 63;
 446     *inseed /= 64;
 447     i++;
 448   } while (*inseed != 0);
 449 }  /*initseed*/
 450
 451
 452 void initjumble(long *inseed, long *inseed0, longer seed, long *njumble)
 453 { /* input number of jumblings for jumble option */
 454   long loopcount;
 455
 456   initseed(inseed, inseed0, seed);
 457   loopcount = 0;
 458   do {
 459     printf("Number of times to jumble?\n");
 460     scanf("%ld%*[^\n]", njumble);
 461     getchar();
 462     countup(&loopcount, 10);
 463   } while ((*njumble) < 1);
 464 }  /*initjumble*/
 465
 466
 467 void initoutgroup(long *outgrno, long spp)
 468 { /* input outgroup number */
 469   long loopcount;
 470   boolean done;
 471
 472   loopcount = 0;
 473   do {
 474     printf("Type number of the outgroup:\n");
 475     scanf("%ld%*[^\n]", outgrno);
 476     getchar();
 477     done = (*outgrno >= 1 && *outgrno <= spp);
 478     if (!done) {
 479       printf("BAD OUTGROUP NUMBER: %ld\n", *outgrno);
 480       printf("  Must be in range 1 - %ld\n", spp);
 481     }
 482     countup(&loopcount, 10);
 483   } while (done != true);
 484 }  /*initoutgroup*/
 485
 486
 487 void initthreshold(double *threshold)
 488 { /* input threshold for threshold parsimony option */
 489   long loopcount;
 490   boolean done;
 491
 492   loopcount = 0;
 493   do {
 494     printf("What will be the threshold value?\n");
 495     scanf("%lf%*[^\n]", threshold);
 496     getchar();
 497     done = (*threshold >= 1.0);
 498     if (!done)
 499       printf("BAD THRESHOLD VALUE:  it must be greater than 1\n");
 500     else
 501       *threshold = (long)(*threshold * 10.0 + 0.5) / 10.0;
 502     countup(&loopcount, 10);
 503   } while (done != true);
 504 }  /*initthreshold*/
 505
 506
 507 void initcatn(long *categs)
 508 { /* initialize category number for rate categories */
 509   long loopcount;
 510
 511   loopcount = 0;
 512   *categs = 0;
 513   do {
 514     printf("Number of categories (1-%d)?\n", maxcategs);
 515     scanf("%ld%*[^\n]", categs);
 516     getchar();
 517     countup(&loopcount, 10);
 518   } while (*categs > maxcategs || *categs < 1);
 519 }  /*initcatn*/
 520
 521
 522 void initcategs(long categs, double *rate)
 523 { /* initialize category rates for HMM rates */
 524   long i, loopcount, scanned;
 525   char line[100], rest[100];
 526   boolean done;
 527
 528   loopcount = 0;
 529   for (;;){
 530     printf("Rate for each category? (use a space to separate)\n");
 531     getstryng(line);
 532     done = true;
 533     for (i = 0; i < categs; i++){
 534       scanned = sscanf(line,"%lf %[^\n]", &rate[i],rest);
 535       if ((scanned < 2 && i < (categs - 1)) ||
 536        (scanned < 1 && i == (categs - 1))){
 537      printf("Please enter exactly %ld values.\n",categs);
 538      done = false;
 539      break;
 540       }
 541       strcpy(line,rest);
 542     }
 543     if (done)
 544       break;
 545     countup(&loopcount, 100);
 546   }
 547 }  /*initcategs*/
 548
 549
 550 void initprobcat(long categs, double *probsum, double *probcat)
 551 { /* input probabilities of rate categores for HMM rates */
 552   long i, loopcount, scanned;
 553   boolean done;
 554   char line[100], rest[100];
 555
 556   loopcount = 0;
 557   do {
 558     printf("Probability for each category?");
 559     printf(" (use a space to separate)\n");
 560     getstryng(line);
 561     done = true;
 562     for (i = 0; i < categs; i++){
 563       scanned = sscanf(line,"%lf %[^\n]",&probcat[i],rest);
 564       if ((scanned < 2 && i < (categs - 1)) ||
 565        (scanned < 1 && i == (categs - 1))){
 566      done = false;
 567      printf("Please enter exactly %ld values.\n",categs);
 568      break;}
 569       strcpy(line,rest);
 570     }
 571     if (!done)
 572       continue;
 573     *probsum = 0.0;
 574     for (i = 0; i < categs; i++)
 575       *probsum += probcat[i];
 576     if (fabs(1.0 - (*probsum)) > 0.001) {
 577       done = false;
 578       printf("Probabilities must add up to");
 579       printf(" 1.0, plus or minus 0.001.\n");
 580     }
 581     countup(&loopcount, 100);
 582   } while (!done);
 583 }  /*initprobcat*/
 584
 585
 586 void lgr(long m, double b, raterootarray lgroot)
 587 { /* For use by initgammacat.  Get roots of m-th Generalized Laguerre
 588      polynomial, given roots of (m-1)-th, these are to be
 589      stored in lgroot[m][] */
 590   long i;
 591   double upper, lower, x, y;
 592   boolean dwn;   /* is function declining in this interval? */
 593
 594   if (m == 1) {
 595     lgroot[1][1] = 1.0+b;
 596   } else {
 597     dwn = true;
 598     for (i=1; i<=m; i++) {
 599       if (i < m) {
 600         if (i == 1)
 601           lower = 0.0;
 602         else
 603           lower = lgroot[m-1][i-1];
 604         upper = lgroot[m-1][i];
 605       } else {                 /* i == m, must search above */
 606         lower = lgroot[m-1][i-1];
 607         x = lgroot[m-1][m-1];
 608         do {
 609           x = 2.0*x;
 610           y = glaguerre(m, b,x);
 611         } while ((dwn && (y > 0.0)) || ((!dwn) && (y < 0.0)));
 612         upper = x;
 613       }
 614       while (upper-lower > 0.000000001) {
 615         x = (upper+lower)/2.0;
 616         if (glaguerre(m, b, x) > 0.0) {
 617           if (dwn)
 618             lower = x;
 619           else
 620             upper = x;
 621         } else {
 622           if (dwn)
 623             upper = x;
 624           else
 625             lower = x;
 626         }
 627       }
 628       lgroot[m][i] = (lower+upper)/2.0;
 629       dwn = !dwn;                /* switch for next one */
 630     }
 631   }
 632 } /* lgr */
 633
 634
 635 double logfac (long n)
 636 { /* log(n!) values were calculated with Mathematica
 637      with a precision of 30 digits */
 638   long i;
 639   double x;
 640
 641   switch (n)
 642     {
 643     case 0:
 644       return 0.;
 645     case 1:
 646       return 0.;
 647     case 2:
 648       return 0.693147180559945309417232121458;
 649     case 3:
 650       return 1.791759469228055000812477358381;
 651     case 4:
 652       return 3.1780538303479456196469416013;
 653     case 5:
 654       return 4.78749174278204599424770093452;
 655     case 6:
 656       return 6.5792512120101009950601782929;
 657     case 7:
 658       return 8.52516136106541430016553103635;
 659     case 8:
 660       return 10.60460290274525022841722740072;
 661     case 9:
 662       return 12.80182748008146961120771787457;
 663     case 10:
 664       return 15.10441257307551529522570932925;
 665     case 11:
 666       return 17.50230784587388583928765290722;
 667     case 12:
 668       return 19.98721449566188614951736238706;
 669     default:
 670       x = 19.98721449566188614951736238706;
 671       for (i = 13; i <= n; i++)
 672         x += log(i);
 673       return x;
 674     }
 675 }
 676
 677
 678 double glaguerre(long m, double b, double x)
 679 { /* Generalized Laguerre polynomial computed recursively.
 680      For use by initgammacat */
 681   long i;
 682   double gln, glnm1, glnp1; /* L_n, L_(n-1), L_(n+1) */
 683
 684   if (m == 0)
 685     return 1.0;
 686   else {
 687     if (m == 1)
 688       return 1.0 + b - x;
 689     else {
 690       gln = 1.0+b-x;
 691       glnm1 = 1.0;
 692       for (i=2; i <= m; i++) {
 693         glnp1 = ((2*(i-1)+b+1.0-x)*gln - (i-1+b)*glnm1)/i;
 694         glnm1 = gln;
 695         gln = glnp1;
 696       }
 697       return gln;
 698     }
 699   }
 700 } /* glaguerre */
 701
 702
 703 void initlaguerrecat(long categs, double alpha, double *rate, double *probcat)
 704 { /* calculate rates and probabilities to approximate Gamma distribution
 705      of rates with "categs" categories and shape parameter "alpha" using
 706      rates and weights from Generalized Laguerre quadrature */
 707   long i;
 708   raterootarray lgroot; /* roots of GLaguerre polynomials */
 709   double f, x, xi, y;
 710
 711   alpha = alpha - 1.0;
 712   lgroot[1][1] = 1.0+alpha;
 713   for (i = 2; i <= categs; i++)
 714     lgr(i, alpha, lgroot);                   /* get roots for L^(a)_n */
 715   /* here get weights */
 716   /* Gamma weights are (1+a)(1+a/2) ... (1+a/n)*x_i/((n+1)^2 [L_{n+1}^a(x_i)]^2)  */
 717   f = 1;
 718   for (i = 1; i <= categs; i++)
 719     f *= (1.0+alpha/i);
 720   for (i = 1; i <= categs; i++) {
 721     xi = lgroot[categs][i];
 722     y = glaguerre(categs+1, alpha, xi);
 723     x = f*xi/((categs+1)*(categs+1)*y*y);
 724     rate[i-1] = xi/(1.0+alpha);
 725     probcat[i-1] = x;
 726   }
 727 } /* initlaguerrecat */
 728
 729
 730 double hermite(long n, double x)
 731 { /* calculates hermite polynomial with degree n and parameter x */
 732   /* seems to be unprecise for n>13 -> root finder does not converge*/
 733   double h1 = 1.;
 734   double h2 = 2. * x;
 735   double xx = 2. * x;
 736   long i;
 737
 738   for (i = 1; i < n; i++) {
 739     xx = 2. * x * h2 - 2. * (i) * h1;
 740     h1 = h2;
 741     h2 = xx;
 742   }
 743   return xx;
 744 } /* hermite */
 745
 746
 747 void root_hermite(long n, double *hroot)
 748 { /* find roots of Hermite polynmials */
 749   long z;
 750   long ii;
 751   long start;
 752
 753   if (n % 2 == 0) {
 754     start = n/2;
 755     z = 1;
 756   } else {
 757     start = n/2 + 1;
 758     z=2;
 759     hroot[start-1] = 0.0;
 760   }
 761   for (ii = start; ii < n; ii++) {         /* search only upwards*/
 762     hroot[ii] = halfroot(hermite,n,hroot[ii-1]+EPSILON, 1./n);
 763     hroot[start - z] = -hroot[ii];
 764     z++;
 765   }
 766 } /* root_hermite */
 767
 768
 769 double halfroot(double (*func)(long m, double x), long n, double startx,
 770                 double delta)
 771 { /* searches from the bound (startx) only in one direction
 772      (by positive or negative delta, which results in
 773      other-bound=startx+delta)
 774      delta should be small.
 775      (*func) is a function with two arguments  */
 776   double xl;
 777   double xu;
 778   double xm;
 779   double fu;
 780   double fl;
 781   double fm = 100000.;
 782   double gradient;
 783   boolean dwn;
 784
 785   /* decide if we search above or below startx and escapes to trace back
 786      to the starting point that most often will be
 787      the root from the previous calculation */
 788   if (delta < 0) {
 789     xu = startx;
 790     xl = xu + delta;
 791   } else {
 792     xl = startx;
 793     xu = xl + delta;
 794   }
 795   delta = fabs(delta);
 796   fu = (*func)(n, xu);
 797   fl = (*func)(n, xl);
 798   gradient = (fl-fu)/(xl-xu);
 799   while(fabs(fm) > EPSILON) {        /* is root outside of our bracket?*/
 800     if ((fu<0.0 && fl<0.0) || (fu>0.0 && fl > 0.0)) {
 801       xu += delta;
 802       fu = (*func)(n, xu);
 803       fl = (*func)(n, xl);
 804       gradient = (fl-fu)/(xl-xu);
 805       dwn = (gradient < 0.0) ? true : false;
 806     } else {
 807       xm = xl - fl / gradient;
 808       fm = (*func)(n, xm);
 809       if (dwn) {
 810         if (fm > 0.) {
 811           xl = xm;
 812           fl = fm;
 813         } else {
 814           xu = xm;
 815           fu = fm;
 816         }
 817       } else {
 818         if (fm > 0.) {
 819           xu = xm;
 820           fu = fm;
 821         } else {
 822           xl = xm;
 823           fl = fm;
 824         }
 825       }
 826       gradient = (fl-fu)/(xl-xu);
 827     }
 828   }
 829   return xm;
 830 } /* halfroot */
 831
 832
 833 void hermite_weight(long n, double * hroot, double * weights)
 834 {
 835   /* calculate the weights for the hermite polynomial at the roots
 836      using formula from Abramowitz and Stegun chapter 25.4.46 p.890 */
 837   long i;
 838   double hr2;
 839   double numerator;
 840
 841   numerator = exp(0.6931471805599 * ( n-1.) + logfac(n)) / (n*n);
 842   for (i = 0; i < n; i++) {
 843     hr2 = hermite(n-1, hroot[i]);
 844     weights[i] = numerator / (hr2*hr2);
 845   }
 846 } /* hermiteweight */
 847
 848
 849 void inithermitcat(long categs, double alpha, double *rate, double *probcat)
 850 { /* calculates rates and probabilities */
 851   long i;
 852   double *hroot;
 853   double std;
 854
 855   std = SQRT2 /sqrt(alpha);
 856   hroot = (double *) Malloc((categs+1) * sizeof(double));
 857   root_hermite(categs, hroot);         /* calculate roots */
 858   hermite_weight(categs, hroot, probcat);  /* set weights */
 859   for (i=0; i<categs; i++) {           /* set rates */
 860     rate[i] = 1.0 + std * hroot[i];
 861     probcat[i] = probcat[i];
 862   }
 863   free(hroot);
 864 } /* inithermitcat */
 865
 866
 867 void initgammacat (long categs, double alpha, double *rate, double *probcat)
 868 { /* calculate rates and probabilities to approximate Gamma distribution
 869    of rates with "categs" categories and shape parameter "alpha" using
 870    rates and weights from Generalized Laguerre quadrature or from
 871    Hermite quadrature */
 872
 873   if (alpha >= 100.0)
 874     inithermitcat(categs, alpha, rate, probcat);
 875   else
 876     initlaguerrecat(categs, alpha, rate, probcat);
 877 } /* initgammacat */
 878
 879
 880 void inithowmany(long *howmanny, long howoften)
 881 {/* input how many cycles */
 882   long loopcount;
 883
 884   loopcount = 0;
 885   do {
 886     printf("How many cycles of %4ld trees?\n", howoften);
 887     scanf("%ld%*[^\n]", howmanny);
 888     getchar();
 889     countup(&loopcount, 10);
 890   } while (*howmanny <= 0);
 891 }  /*inithowmany*/
 892
 893
 894
 895 void inithowoften(long *howoften)
 896 { /* input how many trees per cycle */
 897   long loopcount;
 898
 899   loopcount = 0;
 900   do {
 901     printf("How many trees per cycle?\n");
 902     scanf("%ld%*[^\n]", howoften);
 903     getchar();
 904     countup(&loopcount, 10);
 905   } while (*howoften <= 0);
 906 }  /*inithowoften*/
 907
 908
 909 void initlambda(double *lambda)
 910 { /* input patch length parameter for autocorrelated HMM rates */
 911   long loopcount;
 912
 913   loopcount = 0;
 914   do {
 915     printf("Mean block length of sites having the same rate (greater than 1)?\n");
 916     scanf("%lf%*[^\n]", lambda);
 917     getchar();
 918     countup(&loopcount, 10);
 919   } while (*lambda <= 1.0);
 920   *lambda = 1.0 / *lambda;
 921 }  /*initlambda*/
 922
 923
 924 void initfreqs(double *freqa, double *freqc, double *freqg, double *freqt)
 925 { /* input frequencies of the four bases */
 926   char input[100];
 927   long scanned, loopcount;
 928
 929   printf("Base frequencies for A, C, G, T/U (use blanks to separate)?\n");
 930   loopcount = 0;
 931   do {
 932     getstryng(input);
 933     scanned = sscanf(input,"%lf%lf%lf%lf%*[^\n]", freqa, freqc, freqg, freqt);
 934     if (scanned == 4)
 935       break;
 936     else
 937       printf("Please enter exactly 4 values.\n");
 938     countup(&loopcount, 100);
 939   } while (1);
 940 }  /* initfreqs */
 941
 942
 943 void initratio(double *ttratio)
 944 { /* input transition/transversion ratio */
 945   long loopcount;
 946
 947   loopcount = 0;
 948   do {
 949     printf("Transition/transversion ratio?\n");
 950     scanf("%lf%*[^\n]", ttratio);
 951     getchar();
 952     countup(&loopcount, 10);
 953   } while (*ttratio < 0.0);
 954 }  /* initratio */
 955
 956
 957 void initpower(double *power)
 958 {
 959   printf("New power?\n");
 960   scanf("%lf%*[^\n]", power);
 961   getchar();
 962 }  /*initpower*/
 963
 964
 965 void initdatasets(long *datasets)
 966 {
 967   /* handle multi-data set option */
 968   long loopcount;
 969   boolean done;
 970
 971   loopcount = 0;
 972   do {
 973     printf("How many data sets?\n");
 974     scanf("%ld%*[^\n]", datasets);
 975     getchar();
 976     done = (*datasets > 1);
 977       if (!done)
 978      printf("Bad data sets number:  it must be greater than 1\n");
 979     countup(&loopcount, 10);
 980   } while (!done);
 981 } /* initdatasets */
 982
 983
 984 void justweights(long *datasets)
 985 {
 986   /* handle multi-data set option by weights */
 987   long loopcount;
 988   boolean done;
 989
 990   loopcount = 0;
 991   do {
 992     printf("How many sets of weights?\n");
 993     scanf("%ld%*[^\n]", datasets);
 994     getchar();
 995     done = (*datasets >= 1);
 996       if (!done)
 997      printf("BAD NUMBER:  it must be greater than 1\n");
 998     countup(&loopcount, 10);
 999   } while (!done);
1000 } /* justweights */
1001
1002
1003 void initterminal(boolean *ibmpc, boolean *ansi)
1004 {
1005   /* handle terminal option */
1006   if (*ibmpc) {
1007     *ibmpc = false;
1008     *ansi = true;
1009   } else if (*ansi)
1010       *ansi = false;
1011     else
1012       *ibmpc = true;
1013 }  /*initterminal*/
1014
1015
1016 void initnumlines(long *screenlines)
1017 {
1018   long loopcount;
1019
1020   loopcount = 0;
1021   do {
1022     *screenlines = readlong("Number of lines on screen?\n");
1023     countup(&loopcount, 10);
1024   } while (*screenlines <= 12);
1025 }  /*initnumlines*/
1026
1027
1028 void initbestrees(bestelm *bestrees, long maxtrees, boolean glob)
1029 {
1030   /* initializes either global or local field of each array in bestrees */
1031   long i;
1032
1033   if (glob)
1034     for (i = 0; i < maxtrees; i++)
1035       bestrees[i].gloreange = false;
1036   else
1037     for (i = 0; i < maxtrees; i++)
1038       bestrees[i].locreange = false;
1039 } /* initbestrees */
1040
1041
1042 void newline(FILE *filename, long i, long j, long k)
1043 {
1044   /* go to new line if i is a multiple of j, indent k spaces */
1045   long m;
1046
1047   if ((i - 1) % j != 0 || i <= 1)
1048     return;
1049   putc('\n', filename);
1050   for (m = 1; m <= k; m++)
1051     putc(' ', filename);
1052 }  /* newline */
1053
1054
1055 void inputnumbersold(long *spp, long *chars, long *nonodes, long n)
1056 {
1057   /* input the numbers of species and of characters */
1058
1059   if (fscanf(infile, "%ld%ld", spp, chars) != 2 || *spp <= 0 || *chars <= 0) {
1060     printf(
1061     "ERROR: Unable to read the number of species or characters in data set\n");
1062     printf(
1063       "The input file is incorrect (perhaps it was not saved text only).\n");
1064   }
1065   *nonodes = *spp * 2 - n;
1066 }  /* inputnumbersold */
1067
1068
1069 void inputnumbers(long *spp, long *chars, long *nonodes, long n)
1070 {
1071   /* input the numbers of species and of characters */
1072
1073   if (fscanf(infile, "%ld%ld", spp, chars) != 2 || *spp <= 0 || *chars <= 0) {
1074     printf(
1075     "ERROR: Unable to read the number of species or characters in data set\n");
1076     printf(
1077       "The input file is incorrect (perhaps it was not saved text only).\n");
1078   }
1079   *nonodes = *spp * 2 - n;
1080 }  /* inputnumbers */
1081
1082
1083 void inputnumbers2(long *spp, long *nonodes, long n)
1084 {
1085   /* read species number */
1086
1087   if (fscanf(infile, "%ld", spp) != 1 || *spp <= 0) {
1088     printf("ERROR: Unable to read the number of species in data set\n");
1089     printf(
1090       "The input file is incorrect (perhaps it was not saved text only).\n");
1091   }
1092   fprintf(outfile, "\n%4ld Populations\n", *spp);
1093   *nonodes = *spp * 2 - n;
1094 }  /* inputnumbers2 */
1095
1096
1097 void inputnumbers3(long *spp, long *chars)
1098 {
1099   /* input the numbers of species and of characters */
1100
1101   if (fscanf(infile, "%ld%ld", spp, chars) != 2 || *spp <= 0 || *chars <= 0) {
1102     printf(
1103     "ERROR: Unable to read the number of species or characters in data set\n");
1104     printf(
1105        "The input file is incorrect (perhaps it was not saved text only).\n");
1106     exxit(-1);
1107   }
1108 }  /* inputnumbers3 */
1109
1110
1111 void samenumsp(long *chars, long ith)
1112 {
1113   /* check if spp is same as the first set in other data sets */
1114   long cursp, curchs;
1115
1116   if (eoln(infile))
1117     scan_eoln(infile);
1118   fscanf(infile, "%ld%ld", &cursp, &curchs);
1119   if (cursp != spp) {
1120     printf(
1121          "\n\nERROR: Inconsistent number of species in data set %ld\n\n", ith);
1122     exxit(-1);
1123   }
1124   *chars = curchs;
1125 } /* samenumsp */
1126
1127
1128 void samenumsp2(long ith)
1129 {
1130   /* check if spp is same as the first set in other data sets */
1131   long cursp;
1132
1133   if (eoln(infile))
1134     scan_eoln(infile);
1135   if (fscanf(infile, "%ld", &cursp) != 1) {
1136     printf("\n\nERROR: Unable to read number of species in data set %ld\n",
1137       ith);
1138     printf(
1139       "The input file is incorrect (perhaps it was not saved text only).\n");
1140     exxit(-1);
1141   }
1142   if (cursp != spp) {
1143     printf(
1144       "\n\nERROR: Inconsistent number of species in data set %ld\n\n", ith);
1145     exxit(-1);
1146   }
1147 } /* samenumsp2 */
1148
1149
1150 void readoptions(long *extranum, const char *options)
1151 { /* read option characters from input file */
1152   Char ch;
1153
1154   while (!(eoln(infile))) {
1155     ch = gettc(infile);
1156     uppercase(&ch);
1157     if (strchr(options, ch) != NULL)
1158      (* extranum)++;
1159     else if (!(ch == ' ' || ch == '\t')) {
1160       printf("BAD OPTION CHARACTER: %c\n", ch);
1161       exxit(-1);
1162     }
1163   }
1164   scan_eoln(infile);
1165 }  /* readoptions */
1166
1167
1168 void matchoptions(Char *ch, const char *options)
1169 {  /* match option characters to those in auxiliary options line */
1170
1171   *ch = gettc(infile);
1172   uppercase(ch);
1173   if (strchr(options, *ch) == NULL) {
1174     printf("ERROR: Incorrect auxiliary options line");
1175     printf(" which starts with %c\n", *ch);
1176     exxit(-1);
1177   }
1178 }  /* matchoptions */
1179
1180
1181 void inputweightsold(long chars, steptr weight, boolean *weights)
1182 {
1183   Char ch;
1184   int i;
1185
1186   for (i = 1; i < nmlngth ; i++)
1187     getc(infile);
1188
1189   for (i = 0; i < chars; i++) {
1190     do {
1191       if (eoln(infile))
1192         scan_eoln(infile);
1193       ch = gettc(infile);
1194       if (ch == '\n')
1195         ch = ' ';
1196     } while (ch == ' ');
1197     weight[i] = 1;
1198     if (isdigit(ch))
1199       weight[i] = ch - '0';
1200     else if (isalpha(ch)) {
1201       uppercase(&ch);
1202       weight[i] = ch - 'A' + 10;
1203     } else {
1204       printf("\n\nERROR: Bad weight character: %c\n\n", ch);
1205       exxit(-1);
1206     }
1207   }
1208   scan_eoln(infile);
1209   *weights = true;
1210 } /*inputweightsold*/
1211
1212
1213 void inputweights(long chars, steptr weight, boolean *weights)
1214 {
1215   /* input the character weights, 0-9 and A-Z for weights 0 - 35 */
1216   Char ch;
1217   long i;
1218
1219   for (i = 0; i < chars; i++) {
1220     do {
1221       if (eoln(weightfile))
1222         scan_eoln(weightfile);
1223       ch = gettc(weightfile);
1224       if (ch == '\n')
1225         ch = ' ';
1226     } while (ch == ' ');
1227     weight[i] = 1;
1228     if (isdigit(ch))
1229       weight[i] = ch - '0';
1230     else if (isalpha(ch)) {
1231       uppercase(&ch);
1232       weight[i] = ch - 'A' + 10;
1233     } else {
1234       printf("\n\nERROR: Bad weight character: %c\n\n", ch);
1235       exxit(-1);
1236     }
1237   }
1238   scan_eoln(weightfile);
1239   *weights = true;
1240 }  /* inputweights */
1241
1242
1243 void inputweights2(long a, long b, long *weightsum,
1244         steptr weight, boolean *weights, const char *prog)
1245 {
1246   /* input the character weights, 0 or 1 */
1247   Char ch;
1248   long i;
1249
1250   *weightsum = 0;
1251   for (i = a; i < b; i++) {
1252     do {
1253       if (eoln(weightfile))
1254         scan_eoln(weightfile);
1255       ch = gettc(weightfile);
1256     } while (ch == ' ');
1257     weight[i] = 1;
1258     if (ch == '0' || ch == '1')
1259       weight[i] = ch - '0';
1260     else {
1261       printf("\n\nERROR: Bad weight character: %c -- ", ch);
1262       printf("weights in %s must be 0 or 1\n", prog);
1263       exxit(-1);
1264     }
1265     *weightsum += weight[i];
1266   }
1267   *weights = true;
1268   scan_eoln(weightfile);
1269 }  /* inputweights2 */
1270
1271
1272 void printweights(FILE *filename, long inc, long chars,
1273         steptr weight, const char *letters)
1274 {
1275   /* print out the weights of sites */
1276   long i, j;
1277   boolean letterweights;
1278
1279   letterweights = false;
1280   for (i = 0; i < chars; i++)
1281     if (weight[i] > 9)
1282       letterweights = true;
1283   fprintf(filename, "\n    %s are weighted as follows:",letters);
1284   if (letterweights)
1285     fprintf(filename, " (A = 10, B = 11, etc.)\n");
1286   else
1287     putc('\n', filename);
1288   for (i = 0; i < chars; i++) {
1289     if (i % 60 == 0) {
1290       putc('\n', filename);
1291       for (j = 1; j <= nmlngth + 3; j++)
1292         putc(' ', filename);
1293     }
1294     if (weight[i+inc] < 10)
1295       fprintf(filename, "%ld", weight[i + inc]);
1296     else
1297       fprintf(filename, "%c", 'A'-10+(int)weight[i + inc]);
1298     if ((i+1) % 5 == 0 && (i+1) % 60 != 0)
1299       putc(' ', filename);
1300   }
1301   fprintf(filename, "\n\n");
1302 }  /* printweights */
1303
1304
1305 void inputcategs(long a, long b, steptr category, long categs,const char *prog)
1306 {
1307   /* input the categories, 1-9 */
1308   Char ch;
1309   long i;
1310
1311   for (i = a; i < b; i++) {
1312     do {
1313       if (eoln(catfile))
1314         scan_eoln(catfile);
1315       ch = gettc(catfile);
1316     } while (ch == ' ');
1317     if ((ch >= '1') && (ch <= ('0'+categs)))
1318       category[i] = ch - '0';
1319     else {
1320       printf("\n\nERROR: Bad category character: %c", ch);
1321       printf(" -- categories in %s are currently 1-%ld\n", prog, categs);
1322       exxit(-1);
1323     }
1324   }
1325   scan_eoln(catfile);
1326 }  /* inputcategs */
1327
1328
1329 void printcategs(FILE *filename, long chars, steptr category,
1330                  const char *letters)
1331 {
1332   /* print out the sitewise categories */
1333   long i, j;
1334
1335   fprintf(filename, "\n    %s are:\n",letters);
1336   for (i = 0; i < chars; i++) {
1337     if (i % 60 == 0) {
1338       putc('\n', filename);
1339       for (j = 1; j <= nmlngth + 3; j++)
1340         putc(' ', filename);
1341     }
1342     fprintf(filename, "%ld", category[i]);
1343     if ((i+1) % 10 == 0 && (i+1) % 60 != 0)
1344       putc(' ', filename);
1345   }
1346   fprintf(filename, "\n\n");
1347 }  /* printcategs */
1348
1349
1350 void inputfactors(long chars, Char *factor, boolean *factors)
1351 {
1352   /* reads the factor symbols */
1353   long i;
1354
1355   for (i = 0; i < (chars); i++) {
1356     if (eoln(factfile))
1357       scan_eoln(factfile);
1358     factor[i] = gettc(factfile);
1359     if (factor[i] == '\n')
1360       factor[i] = ' ';
1361   }
1362   scan_eoln(factfile);
1363   *factors = true;
1364 }  /* inputfactors */
1365
1366
1367 void printfactors(FILE *filename, long chars, Char *factor, const char *letters)
1368 {
1369   /* print out list of factor symbols */
1370   long i;
1371
1372   fprintf(filename, "Factors%s:\n\n", letters);
1373   for (i = 1; i <= nmlngth - 5; i++)
1374     putc(' ', filename);
1375   for (i = 1; i <= (chars); i++) {
1376     newline(filename, i, 55, nmlngth + 3);
1377     putc(factor[i - 1], filename);
1378     if (i % 5 == 0)
1379       putc(' ', filename);
1380   }
1381   putc('\n', filename);
1382 }  /* printfactors */
1383
1384
1385 void headings(long chars, const char *letters1, const char *letters2)
1386 {
1387   long i, j;
1388
1389   putc('\n', outfile);
1390   j = nmlngth + (chars + (chars - 1) / 10) / 2 - 5;
1391   if (j < nmlngth - 1)
1392     j = nmlngth - 1;
1393   if (j > 37)
1394     j = 37;
1395   fprintf(outfile, "Name");
1396   for (i = 1; i <= j; i++)
1397     putc(' ', outfile);
1398   fprintf(outfile, "%s\n", letters1);
1399   fprintf(outfile, "----");
1400   for (i = 1; i <= j; i++)
1401     putc(' ', outfile);
1402   fprintf(outfile, "%s\n\n", letters2);
1403 }  /* headings */
1404
1405
1406 void initname(long i)
1407 {
1408   /* read in species name */
1409   long j;
1410
1411   for (j = 0; j < nmlngth; j++) {
1412     if (eoff(infile) | eoln(infile)){
1413       printf("\n\nERROR: end-of-line or end-of-file");
1414       printf(" in the middle of species name for species %ld\n\n", i+1);
1415       exxit(-1);
1416     }
1417     nayme[i][j] = gettc(infile);
1418     if ((nayme[i][j] == '(') || (nayme[i][j] == ')') || (nayme[i][j] == ':')
1419         || (nayme[i][j] == ',') || (nayme[i][j] == ';') || (nayme[i][j] == '[')
1420         || (nayme[i][j] == ']')) {
1421       printf("\nERROR: Species name may not contain characters ( ) : ; , [ ] \n");
1422       printf("       In name of species number %ld there is character %c\n\n",
1423               i+1, nayme[i][j]);
1424       exxit(-1);
1425     }
1426   }
1427 } /* initname */
1428
1429
1430 void findtree(boolean *found,long *pos,long nextree,long *place,bestelm *bestrees)
1431 {
1432   /* finds tree given by array place in array bestrees by binary search */
1433   /* used by dnacomp, dnapars, dollop, mix, & protpars */
1434   long i, lower, upper;
1435   boolean below, done;
1436
1437   below = false;
1438   lower = 1;
1439   upper = nextree - 1;
1440   (*found) = false;
1441   while (!(*found) && lower <= upper) {
1442     (*pos) = (lower + upper) / 2;
1443     i = 3;
1444     done = false;
1445     while (!done) {
1446       done = (i > spp);
1447       if (!done)
1448         done = (place[i - 1] != bestrees[(*pos) - 1].btree[i - 1]);
1449       if (!done)
1450         i++;
1451     }
1452     (*found) = (i > spp);
1453     if (*found)
1454       break;
1455     below = (place[i - 1] <  bestrees[(*pos )- 1].btree[i - 1]);
1456     if (below)
1457       upper = (*pos) - 1;
1458     else
1459       lower = (*pos) + 1;
1460   }
1461   if (!(*found) && !below)
1462     (*pos)++;
1463 }  /* findtree */
1464
1465
1466 void addtree(long pos,long *nextree,boolean collapse,long *place,bestelm *bestrees)
1467 {
1468   /* puts tree from array place in its proper position in array bestrees */
1469   /* used by dnacomp, dnapars, dollop, mix, & protpars */
1470   long i;
1471
1472   for (i = *nextree - 1; i >= pos; i--){
1473     memcpy(bestrees[i].btree, bestrees[i - 1].btree, spp * sizeof(long));
1474     bestrees[i].gloreange = bestrees[i - 1].gloreange;
1475     bestrees[i - 1].gloreange = false;
1476     bestrees[i].locreange = bestrees[i - 1].locreange;
1477     bestrees[i - 1].locreange = false;
1478     bestrees[i].collapse = bestrees[i - 1].collapse;
1479   }
1480   for (i = 0; i < spp; i++)
1481     bestrees[pos - 1].btree[i] = place[i];
1482     bestrees[pos - 1].collapse = collapse;
1483   (*nextree)++;
1484 }  /* addtree */
1485
1486
1487 long findunrearranged(bestelm *bestrees, long nextree, boolean glob)
1488 {
1489   /* finds bestree with either global or local field false */
1490   long i;
1491
1492   if (glob) {
1493     for (i = 0; i < nextree - 1; i++)
1494       if (!bestrees[i].gloreange)
1495         return i;
1496   } else {
1497     for (i = 0; i < nextree - 1; i++)
1498       if (!bestrees[i].locreange)
1499         return i;
1500   }
1501   return -1;
1502 } /* findunrearranged */
1503
1504
1505 boolean torearrange(bestelm *bestrees, long nextree)
1506 { /* sees if any best tree is yet to be rearranged */
1507
1508   if (findunrearranged(bestrees, nextree, true) >= 0)
1509     return true;
1510   else if (findunrearranged(bestrees, nextree, false) >= 0)
1511     return true;
1512   else
1513     return false;
1514 } /* torearrange */
1515
1516
1517 void reducebestrees(bestelm *bestrees, long *nextree)
1518 {
1519   /* finds best trees with collapsible branches and deletes them */
1520   long i, j;
1521
1522   i = 0;
1523   j = *nextree - 2;
1524   do {
1525     while (!bestrees[i].collapse && i < *nextree - 1) i++;
1526     while (bestrees[j].collapse && j >= 0) j--;
1527     if (i < j) {
1528       memcpy(bestrees[i].btree, bestrees[j].btree, spp * sizeof(long));
1529       bestrees[i].gloreange = bestrees[j].gloreange;
1530       bestrees[i].locreange = bestrees[j].locreange;
1531       bestrees[i].collapse = false;
1532       bestrees[j].collapse = true;
1533     }
1534   } while (i < j);
1535   *nextree = i + 1;
1536 } /* reducebestrees */
1537
1538
1539 void shellsort(double *a, long *b, long n)
1540 { /* Shell sort keeping a, b in same order */
1541   /* used by dnapenny, dolpenny, & penny */
1542   long gap, i, j, itemp;
1543   double rtemp;
1544
1545   gap = n / 2;
1546   while (gap > 0) {
1547     for (i = gap + 1; i <= n; i++) {
1548       j = i - gap;
1549       while (j > 0) {
1550      if (a[j - 1] > a[j + gap - 1]) {
1551        rtemp = a[j - 1];
1552        a[j - 1] = a[j + gap - 1];
1553        a[j + gap - 1] = rtemp;
1554        itemp = b[j - 1];
1555        b[j - 1] = b[j + gap - 1];
1556        b[j + gap - 1] = itemp;
1557      }
1558      j -= gap;
1559       }
1560     }
1561     gap /= 2;
1562   }
1563 }  /* shellsort */
1564
1565
1566 void getch(Char *c, long *parens, FILE *treefile)
1567 { /* get next nonblank character */
1568
1569   do {
1570     if (eoln(treefile))
1571       scan_eoln(treefile);
1572     (*c) = gettc(treefile);
1573
1574     if ((*c) == '\n' || (*c) == '\t')
1575       (*c) = ' ';
1576   } while ( *c == ' ' && !eoff(treefile) );
1577   if ((*c) == '(')
1578     (*parens)++;
1579   if ((*c) == ')')
1580     (*parens)--;
1581 }  /* getch */
1582
1583
1584 void getch2(Char *c, long *parens)
1585 { /* get next nonblank character */
1586   do {
1587     if (eoln(intree))
1588       scan_eoln(intree);
1589     *c = gettc(intree);
1590     if (*c == '\n' || *c == '\t')
1591       *c = ' ';
1592   } while (!(*c != ' ' || eoff(intree)));
1593   if (*c == '(')
1594    (*parens)++;
1595   if (*c == ')')
1596     (*parens)--;
1597 }  /* getch2 */
1598
1599
1600 void findch(Char c, Char *ch, long which)
1601 { /* scan forward until find character c */
1602   boolean done;
1603   long dummy_parens;
1604   done = false;
1605   while (!done) {
1606     if (c == ',') {
1607       if (*ch == '(' || *ch == ')' || *ch == ';') {
1608         printf(
1609       "\n\nERROR in user tree %ld: unmatched parenthesis or missing comma\n\n",
1610           which);
1611         exxit(-1);
1612       } else if (*ch == ',')
1613         done = true;
1614     } else if (c == ')') {
1615       if (*ch == '(' || *ch == ',' || *ch == ';') {
1616         printf("\n\nERROR in user tree %ld: ", which);
1617         printf("unmatched parenthesis or non-bifurcated node\n\n");
1618         exxit(-1);
1619       } else {
1620         if (*ch == ')')
1621           done = true;
1622       }
1623     } else if (c == ';') {
1624       if (*ch != ';') {
1625         printf("\n\nERROR in user tree %ld: ", which);
1626         printf("unmatched parenthesis or missing semicolon\n\n");
1627         exxit(-1);
1628       } else
1629         done = true;
1630     }
1631     if (*ch != ')' && done)
1632       continue;
1633    getch(ch, &dummy_parens, intree);
1634   }
1635 }  /* findch */
1636
1637
1638 void findch2(Char c, long *lparens, long *rparens, Char *ch)
1639 { /* skip forward in user tree until find character c */
1640   boolean done;
1641   long dummy_parens;
1642   done = false;
1643   while (!done) {
1644     if (c == ',') {
1645       if (*ch == '(' || *ch == ')' || *ch == ':' || *ch == ';') {
1646         printf("\n\nERROR in user tree: ");
1647         printf("unmatched parenthesis, missing comma");
1648         printf(" or non-trifurcated base\n\n");
1649      exxit(-1);
1650       } else if (*ch == ',')
1651         done = true;
1652     } else if (c == ')') {
1653       if (*ch == '(' || *ch == ',' || *ch == ':' || *ch == ';') {
1654         printf(
1655    "\n\nERROR in user tree: unmatched parenthesis or non-bifurcated node\n\n");
1656         exxit(-1);
1657       } else if (*ch == ')') {
1658         (*rparens)++;
1659         if ((*lparens) > 0 && (*lparens) == (*rparens)) {
1660           if ((*lparens) == spp - 2) {
1661            getch(ch, &dummy_parens, intree);
1662             if (*ch != ';') {
1663               printf( "\n\nERROR in user tree: ");
1664               printf("unmatched parenthesis or missing semicolon\n\n");
1665               exxit(-1);
1666             }
1667           }
1668         }
1669      done = true;
1670       }
1671     }
1672     if (*ch != ')' && done)
1673       continue;
1674     if (*ch == ')')
1675      getch(ch, &dummy_parens, intree);
1676   }
1677 }  /* findch2 */
1678
1679
1680 void processlength(double *valyew, double *divisor, Char *ch,
1681         boolean *minusread, FILE *treefile, long *parens)
1682 { /* read a branch length from a treefile */
1683   long digit, ordzero;
1684   boolean pointread;
1685
1686   ordzero = '0';
1687   *minusread = false;
1688   pointread = false;
1689   *valyew = 0.0;
1690   *divisor = 1.0;
1691   getch(ch, parens, treefile);
1692   digit = (long)(*ch - ordzero);
1693   while ( ((digit <= 9) && (digit >= 0)) || *ch == '.' || *ch == '-') {
1694     if (*ch == '.' )
1695       pointread = true;
1696     else if (*ch == '-' )
1697       *minusread = true;
1698     else {
1699       *valyew = *valyew * 10.0 + digit;
1700       if (pointread)
1701         *divisor *= 10.0;
1702     }
1703     getch(ch, parens, treefile);
1704     digit = (long)(*ch - ordzero);
1705   }
1706   if (*minusread)
1707     *valyew = -(*valyew);
1708 }  /* processlength */
1709
1710
1711 void writename(long start, long n, long *enterorder)
1712 { /* write species name and number in entry order */
1713   long i, j;
1714
1715   for (i = start; i < start+n; i++) {
1716     printf(" %3ld. ", i+1);
1717     for (j = 0; j < nmlngth; j++)
1718       putchar(nayme[enterorder[i] - 1][j]);
1719     putchar('\n');
1720     fflush(stdout);
1721   }
1722 }  /* writename */
1723
1724
1725 void memerror()
1726 {
1727   printf("Error allocating memory\n");
1728   exxit(-1);
1729 }  /* memerror */
1730
1731
1732 void odd_malloc(long x)
1733 { /* error message if attempt to malloc too little or too much memory */
1734   printf ("ERROR: a function asked for an inappropriate amount of memory:");
1735   printf ("  %ld bytes\n", x);
1736   printf ("       This can mean one of two things:\n");
1737   printf ("       1.  The input file is incorrect");
1738   printf (" (perhaps it was not saved as Text Only),\n");
1739   printf ("       2.  There is a bug in the program.\n");
1740   printf ("       Please check your input file carefully.\n");
1741   printf ("       If it seems to be a bug, please mail joe@gs.washington.edu\n");
1742   printf ("       with the name of the program, your computer system type,\n");
1743   printf ("       a full description of the problem, and with the input data file.\n");
1744   printf ("       (which should be in the body of the message, not as an Attachment).\n");
1745
1746   /* abort() can be used to crash */
1747
1748   exxit(-1);
1749 }
1750
1751
1752 MALLOCRETURN *mymalloc(long x)
1753 { /* wrapper for malloc, allowing error message if too little, too much */
1754   MALLOCRETURN *new_block;
1755
1756   if ((x <= 0) ||
1757       (x > TOO_MUCH_MEMORY))
1758     odd_malloc(x);
1759
1760   new_block = (MALLOCRETURN *)calloc(1,x);
1761
1762   if (!new_block) {
1763     memerror();
1764     return (MALLOCRETURN *) new_block;
1765   } else
1766     return (MALLOCRETURN *) new_block;
1767 } /* mymalloc */
1768
1769
1770 void gnu(node **grbg, node **p)
1771 { /* this and the following are do-it-yourself garbage collectors.
1772      Make a new node or pull one off the garbage list */
1773
1774   if (*grbg != NULL) {
1775     *p = *grbg;
1776     *grbg = (*grbg)->next;
1777   } else
1778     *p = (node *)Malloc(sizeof(node));
1779
1780   (*p)->back       = NULL;
1781   (*p)->next       = NULL;
1782   (*p)->tip        = false;
1783   (*p)->times_in_tree = 0.0;
1784   (*p)->r          = 0.0;
1785   (*p)->theta      = 0.0;
1786   (*p)->x          = NULL;
1787   (*p)->protx           = NULL;        /* for the sake of proml     */
1788 }  /* gnu */
1789
1790
1791 void chuck(node **grbg, node *p)
1792 { /* collect garbage on p -- put it on front of garbage list */
1793   p->back = NULL;
1794   p->next = *grbg;
1795   *grbg = p;
1796 }  /* chuck */
1797
1798
1799 void zeronumnuc(node *p, long endsite)
1800 {
1801   long i,j;
1802
1803   for (i = 0; i < endsite; i++)
1804     for (j = (long)A; j <= (long)O; j++)
1805       p->numnuc[i][j] = 0;
1806 } /* zeronumnuc */
1807
1808
1809 void zerodiscnumnuc(node *p, long endsite)
1810 {
1811   long i,j;
1812
1813   for (i = 0; i < endsite; i++)
1814     for (j = (long)zero; j <= (long)seven; j++)
1815       p->discnumnuc[i][j] = 0;
1816 } /* zerodiscnumnuc */
1817
1818
1819 void allocnontip(node *p, long *zeros, long endsite)
1820 { /* allocate an interior node */
1821   /* used by dnacomp, dnapars, & dnapenny */
1822
1823   p->numsteps = (steptr)Malloc(endsite*sizeof(long));
1824   p->oldnumsteps = (steptr)Malloc(endsite*sizeof(long));
1825   p->base = (baseptr)Malloc(endsite*sizeof(long));
1826   p->oldbase = (baseptr)Malloc(endsite*sizeof(long));
1827   p->numnuc = (nucarray *)Malloc(endsite*sizeof(nucarray));
1828   memcpy(p->base, zeros, endsite*sizeof(long));
1829   memcpy(p->numsteps, zeros, endsite*sizeof(long));
1830   memcpy(p->oldbase, zeros, endsite*sizeof(long));
1831   memcpy(p->oldnumsteps, zeros, endsite*sizeof(long));
1832   zeronumnuc(p, endsite);
1833 }  /* allocnontip */
1834
1835
1836 void allocdiscnontip(node *p, long *zeros, unsigned char *zeros2, long endsite)
1837 { /* allocate an interior node */
1838   /* used by pars */
1839
1840   p->numsteps = (steptr)Malloc(endsite*sizeof(long));
1841   p->oldnumsteps = (steptr)Malloc(endsite*sizeof(long));
1842   p->discbase = (discbaseptr)Malloc(endsite*sizeof(unsigned char));
1843   p->olddiscbase = (discbaseptr)Malloc(endsite*sizeof(unsigned char));
1844   p->discnumnuc = (discnucarray *)Malloc(endsite*sizeof(discnucarray));
1845   memcpy(p->discbase, zeros2, endsite*sizeof(unsigned char));
1846   memcpy(p->numsteps, zeros, endsite*sizeof(long));
1847   memcpy(p->olddiscbase, zeros2, endsite*sizeof(unsigned char));
1848   memcpy(p->oldnumsteps, zeros, endsite*sizeof(long));
1849   zerodiscnumnuc(p, endsite);
1850 }  /* allocdiscnontip */
1851
1852
1853 void allocnode(node **anode, long *zeros, long endsite)
1854 { /* allocate a node */
1855   /* used by dnacomp, dnapars, & dnapenny */
1856
1857   *anode = (node *)Malloc(sizeof(node));
1858   allocnontip(*anode, zeros, endsite);
1859 }  /* allocnode */
1860
1861
1862 void allocdiscnode(node **anode, long *zeros, unsigned char *zeros2,
1863         long endsite)
1864 { /* allocate a node */
1865   /* used by pars */
1866
1867   *anode = (node *)Malloc(sizeof(node));
1868   allocdiscnontip(*anode, zeros, zeros2, endsite);
1869 }  /* allocdiscnontip */
1870
1871
1872 void gnutreenode(node **grbg, node **p, long i, long endsite, long *zeros)
1873 { /* this and the following are do-it-yourself garbage collectors.
1874      Make a new node or pull one off the garbage list */
1875
1876   if (*grbg != NULL) {
1877     *p = *grbg;
1878     *grbg = (*grbg)->next;
1879     memcpy((*p)->numsteps, zeros, endsite*sizeof(long));
1880     memcpy((*p)->oldnumsteps, zeros, endsite*sizeof(long));
1881     memcpy((*p)->base, zeros, endsite*sizeof(long));
1882     memcpy((*p)->oldbase, zeros, endsite*sizeof(long));
1883     zeronumnuc(*p, endsite);
1884   } else
1885     allocnode(p, zeros, endsite);
1886   (*p)->back = NULL;
1887   (*p)->next = NULL;
1888   (*p)->tip = false;
1889   (*p)->visited = false;
1890   (*p)->index = i;
1891   (*p)->numdesc = 0;
1892   (*p)->sumsteps = 0.0;
1893 }  /* gnutreenode */
1894
1895
1896 void gnudisctreenode(node **grbg, node **p, long i,
1897         long endsite, long *zeros, unsigned char *zeros2)
1898 { /* this and the following are do-it-yourself garbage collectors.
1899      Make a new node or pull one off the garbage list */
1900
1901   if (*grbg != NULL) {
1902     *p = *grbg;
1903     *grbg = (*grbg)->next;
1904     memcpy((*p)->numsteps, zeros, endsite*sizeof(long));
1905     memcpy((*p)->oldnumsteps, zeros, endsite*sizeof(long));
1906     memcpy((*p)->discbase, zeros2, endsite*sizeof(unsigned char));
1907     memcpy((*p)->olddiscbase, zeros2, endsite*sizeof(unsigned char));
1908     zerodiscnumnuc(*p, endsite);
1909   } else
1910     allocdiscnode(p, zeros, zeros2, endsite);
1911   (*p)->back = NULL;
1912   (*p)->next = NULL;
1913   (*p)->tip = false;
1914   (*p)->visited = false;
1915   (*p)->index = i;
1916   (*p)->numdesc = 0;
1917   (*p)->sumsteps = 0.0;
1918 }  /* gnudisctreenode */
1919
1920
1921 void chucktreenode(node **grbg, node *p)
1922 { /* collect garbage on p -- put it on front of garbage list */
1923
1924   p->back = NULL;
1925   p->next = *grbg;
1926   *grbg = p;
1927 }  /* chucktreenode */
1928
1929
1930 void setupnode(node *p, long i)
1931 { /* initialization of node pointers, variables */
1932
1933   p->next = NULL;
1934   p->back = NULL;
1935   p->times_in_tree = (double) i * 1.0;
1936   p->index = i;
1937   p->tip = false;
1938 }  /* setupnode */
1939
1940
1941 long count_sibs (node *p)
1942 { /* Count the number of nodes in a ring, return the total number of */
1943   /* nodes excluding the one passed into the function (siblings)     */
1944   node *q;
1945   long return_int = 0;
1946
1947   if (p->tip) {
1948    printf ("Error: the function count_sibs called on a tip.  This is a bug.\n");
1949     exxit (-1);
1950   }
1951
1952   q = p->next;
1953   while (q != p) {
1954     if (q == NULL) {
1955       printf ("Error: a loop of nodes was not closed.\n");
1956       exxit (-1);
1957     } else {
1958       return_int++;
1959       q = q->next;
1960     }
1961   }
1962
1963   return return_int;
1964 }  /* count_sibs */
1965
1966
1967 void inittrav (node *p)
1968 { /* traverse to set pointers uninitialized on inserting */
1969   long i, num_sibs;
1970   node *sib_ptr;
1971
1972   if (p == NULL)
1973     return;
1974   if (p->tip)
1975     return;
1976   num_sibs = count_sibs (p);
1977   sib_ptr  = p;
1978   for (i=0; i < num_sibs; i++) {
1979     sib_ptr              = sib_ptr->next;
1980     sib_ptr->initialized = false;
1981     inittrav(sib_ptr->back);
1982   }
1983 } /* inittrav */
1984
1985
1986 void commentskipper(FILE ***intree, long *bracket)
1987 { /* skip over comment bracket contents in reading tree */
1988   char c;
1989
1990   c = gettc(**intree);
1991
1992   while (c != ']') {
1993
1994     if(feof(**intree)) {
1995       printf("\n\nERROR: Unmatched comment brackets\n\n");
1996       exxit(-1);
1997     }
1998
1999     if(c == '[') {
2000       (*bracket)++;
2001       commentskipper(intree, bracket);
2002     }
2003     c = gettc(**intree);
2004   }
2005   (*bracket)--;
2006 }  /* commentskipper */
2007
2008
2009 long countcomma(FILE **treefile, long *comma)
2010 {
2011   /* Modified by Dan F. 11/10/96 */
2012
2013   /* The next line inserted so this function leaves the file pointing
2014      to where it found it, not just re-winding it. */
2015   long orig_position = ftell (*treefile);
2016
2017   Char c;
2018   long  lparen = 0;
2019   long bracket = 0;
2020   (*comma) = 0;
2021
2022
2023   for (;;){
2024     c = getc(*treefile);
2025     if (feof(*treefile))
2026       break;
2027     if (c == ';')
2028       break;
2029     if (c == ',')
2030       (*comma)++;
2031     if (c == '(')
2032          lparen++;
2033     if (c == '[') {
2034       bracket++;
2035       commentskipper(&treefile, &bracket);
2036     }
2037   }
2038
2039   /* Don't just rewind, */
2040   /* rewind (*treefile); */
2041   /* Re-set to where it pointed when the function was called */
2042
2043   fseek (*treefile, orig_position, SEEK_SET);
2044
2045   return lparen + (*comma);
2046 }  /*countcomma*/
2047 /* countcomma rewritten so it passes back both lparen+comma to allocate nodep
2048    and a pointer to the comma variable.  This allows the tree to know how many
2049    species exist, and the tips to be placed in the front of the nodep array */
2050
2051
2052 long countsemic(FILE **treefile)
2053 { /* Used to determine the number of user trees.  Return
2054      either a: the number of semicolons in the file outside comments
2055      or b: the first integer in the file */
2056   Char c;
2057   long return_val, semic = 0;
2058   long bracket = 0;
2059
2060   /* Eat all whitespace */
2061   c = gettc(*treefile);
2062   while ((c == ' ')  ||
2063       (c == '\t') ||
2064       (c == '\n')) {
2065     c = gettc(*treefile);
2066   }
2067
2068   /* Then figure out if the first non-white character is a digit; if
2069      so, return it */
2070   if (isdigit (c)) {
2071     ungetc(c, *treefile);
2072     fscanf((*treefile), "%ld", &return_val);
2073   } else {
2074
2075     /* Loop past all characters, count the number of semicolons
2076        outside of comments */
2077     for (;;){
2078       c = fgetc(*treefile);
2079       if (feof(*treefile))
2080      break;
2081       if (c == ';')
2082      semic++;
2083       if (c == '[') {
2084      bracket++;
2085      commentskipper(&treefile, &bracket);
2086       }
2087     }
2088     return_val = semic;
2089   }
2090
2091   rewind (*treefile);
2092   return return_val;
2093 }  /* countsemic */
2094
2095
2096 void hookup(node *p, node *q)
2097 { /* hook together two nodes */
2098
2099   p->back = q;
2100   q->back = p;
2101 }  /* hookup */
2102
2103
2104 void link_trees(long local_nextnum, long nodenum, long local_nodenum,
2105         pointarray nodep)
2106 {
2107   if(local_nextnum == 0)
2108     hookup(nodep[nodenum],nodep[local_nodenum]);
2109   else if(local_nextnum == 1)
2110       hookup(nodep[nodenum], nodep[local_nodenum]->next);
2111     else if(local_nextnum == 2)
2112         hookup(nodep[nodenum],nodep[local_nodenum]->next->next);
2113       else
2114         printf("Error in Link_Trees()");
2115 } /* link_trees() */
2116
2117
2118 void allocate_nodep(pointarray *nodep, FILE **treefile, long  *precalc_tips)
2119 { /* pre-compute space and allocate memory for nodep */
2120
2121   long numnodes;      /* returns number commas & (    */
2122   long numcom = 0;        /* returns number commas */
2123
2124   numnodes = countcomma(treefile, &numcom) + 1;
2125   *nodep      = (pointarray)Malloc(2*numnodes*sizeof(node *));
2126
2127   (*precalc_tips) = numcom + 1;        /* this will be used in placing the
2128                                           tip nodes in the front region of
2129                                           nodep.  Used for species check?  */
2130 } /* allocate_nodep -plc */
2131
2132
2133 void malloc_pheno (node *p, long endsite, long rcategs)
2134 { /* Allocate the phenotype arrays; used by dnaml */
2135   long i;
2136
2137   p->x  = (phenotype)Malloc(endsite*sizeof(ratelike));
2138   p->underflows = Malloc(endsite * sizeof(double));
2139   for (i = 0; i < endsite; i++)
2140     p->x[i]  = (ratelike)Malloc(rcategs*sizeof(sitelike));
2141 } /* malloc_pheno */
2142
2143
2144 void malloc_ppheno (node *p,long endsite, long rcategs)
2145 {
2146   /* Allocate the phenotype arrays; used by proml */
2147   long i;
2148
2149   p->protx  = (pphenotype)Malloc(endsite*sizeof(pratelike));
2150   p->underflows  = Malloc(endsite*sizeof(double));
2151
2152   for (i = 0; i < endsite; i++)
2153     p->protx[i]  = (pratelike)Malloc(rcategs*sizeof(psitelike));
2154 } /* malloc_ppheno */
2155
2156
2157 long take_name_from_tree (Char *ch, Char *str, FILE *treefile)
2158 {
2159   /* This loop takes in the name from the tree.
2160      Return the length of the name string.  */
2161
2162   long name_length = 0;
2163
2164   do {
2165     if ((*ch) == '_')
2166       (*ch) = ' ';
2167     str[name_length++] = (*ch);
2168     if (eoln(treefile))
2169       scan_eoln(treefile);
2170     (*ch) = gettc(treefile);
2171     if (*ch == '\n')
2172       *ch = ' ';
2173   } while ((*ch) != ':' && (*ch) != ',' && (*ch) != ')' &&
2174         (*ch) != '[' && (*ch) != ';' && name_length <= MAXNCH);
2175   return name_length;
2176 }  /* take_name_from_tree */
2177
2178
2179 void match_names_to_data (Char *str, pointarray treenode, node **p, long spp)
2180 {
2181   /* This loop matches names taken from treefile to indexed names in
2182      the data file */
2183
2184   boolean found;
2185   long i, n;
2186
2187   n = 1;
2188   do {
2189     found = true;
2190     for (i = 0; i < nmlngth; i++) {
2191       found = (found && ((str[i] == nayme[n - 1][i]) ||
2192         (((nayme[n - 1][i] == '_') && (str[i] == ' ')) ||
2193         ((nayme[n - 1][i] == ' ') && (str[i] == '\0')))));
2194     }
2195
2196     if (found)
2197       *p = treenode[n - 1];
2198     else
2199       n++;
2200
2201   } while (!(n > spp || found));
2202
2203   if (n > spp) {
2204     printf("\n\nERROR: Cannot find species: ");
2205     for (i = 0; (str[i] != '\0') && (i < MAXNCH); i++)
2206       putchar(str[i]);
2207     printf(" in data file\n\n");
2208     exxit(-1);
2209   }
2210 }  /* match_names_to_data */
2211
2212
2213 void addelement(node **p, node *q, Char *ch, long *parens, FILE *treefile,
2214         pointarray treenode, boolean *goteof, boolean *first, pointarray nodep,
2215         long *nextnode, long *ntips, boolean *haslengths, node **grbg,
2216         initptr initnode,boolean unifok, long maxnodes)
2217 {
2218   /* Recursive procedure adds nodes to user-defined tree
2219      This is the main (new) tree-reading procedure */
2220
2221   node *pfirst;
2222   long i, len = 0, nodei = 0;
2223   boolean notlast;
2224   Char str[MAXNCH];
2225   node *r;
2226   long furs = 0;
2227
2228   if ((*ch) == '(') {
2229     (*nextnode)++;          /* get ready to use new interior node */
2230     nodei = *nextnode;      /* do what needs to be done at bottom */
2231     if ( maxnodes != -1 && nodei > maxnodes) {
2232       printf("ERROR in input tree file: Attempting to allocate too\n");
2233       printf("many nodes. This is usually caused by a unifurcation.");
2234       printf("To use this tree with this program  use Retree to read\n");
2235       printf("and write this tree.\n");
2236       exxit(-1);
2237     }
2238     (*initnode)(p, grbg, q, len, nodei, ntips,
2239                   parens, bottom, treenode, nodep, str, ch, treefile);
2240     pfirst      = (*p);
2241     notlast = true;
2242     while (notlast) {          /* loop through immediate descendants */
2243       furs++;
2244       (*initnode)(&(*p)->next, grbg, q,
2245                    len, nodei, ntips, parens, nonbottom, treenode,
2246                    nodep, str, ch, treefile);
2247                        /* ... doing what is done before each */
2248       r = (*p)->next;
2249       getch(ch, parens, treefile);      /* look for next character */
2250
2251        /* handle blank names */
2252       if((*ch) == ',' || (*ch) == ':'){
2253         ungetc((*ch), treefile);
2254         *ch = 0;
2255       } else if((*ch)==')'){
2256         ungetc((*ch), treefile);
2257         (*parens)++;
2258         *ch = 0;
2259       }
2260
2261       addelement(&(*p)->next->back, (*p)->next, ch, parens, treefile,
2262         treenode, goteof, first, nodep, nextnode, ntips,
2263         haslengths, grbg, initnode,unifok,maxnodes);
2264
2265       (*initnode)(&r, grbg, q, len, nodei, ntips,
2266                     parens, hslength, treenode, nodep, str, ch, treefile);
2267                            /* do what is done after each about length */
2268       pfirst->numdesc++;               /* increment number of descendants */
2269       *p = r;                         /* make r point back to p */
2270
2271       if ((*ch) == ')') {
2272         notlast = false;
2273         do {
2274           getch(ch, parens, treefile);
2275         } while ((*ch) != ',' && (*ch) != ')' &&
2276            (*ch) != '[' && (*ch) != ';' && (*ch) != ':');
2277       }
2278     }
2279     if ( furs <= 1 && !unifok ) {
2280       printf("ERROR in input tree file: A Unifurcation was detetected.\n");
2281       printf("To use this tree with this program use retree to read and");
2282       printf(" write this tree\n");
2283       exxit(-1);
2284     }
2285
2286     (*p)->next = pfirst;
2287     (*p)       = pfirst;
2288
2289   } else if ((*ch) != ')') {       /* if it's a species name */
2290     for (i = 0; i < MAXNCH; i++)   /* fill string with nulls */
2291       str[i] = '\0';
2292
2293     len = take_name_from_tree (ch, str, treefile);  /* get the name */
2294
2295     if ((*ch) == ')')
2296       (*parens)--;         /* decrement count of open parentheses */
2297     (*initnode)(p, grbg, q, len, nodei, ntips,
2298                   parens, tip, treenode, nodep, str, ch, treefile);
2299                               /* do what needs to be done at a tip */
2300   } else
2301     getch(ch, parens, treefile);
2302   if (q != NULL)
2303     hookup(q, (*p));                    /* now hook up */
2304   (*initnode)(p, grbg, q, len, nodei, ntips,
2305                 parens, iter, treenode, nodep, str, ch, treefile);
2306                               /* do what needs to be done to variable iter */
2307   if ((*ch) == ':')
2308     (*initnode)(p, grbg, q, len, nodei, ntips,
2309                   parens, length, treenode, nodep, str, ch, treefile);
2310                                    /* do what needs to be done with length */
2311   else if ((*ch) != ';' && (*ch) != '[')
2312     (*initnode)(p, grbg, q, len, nodei, ntips,
2313                   parens, hsnolength, treenode, nodep, str, ch, treefile);
2314                              /* ... or what needs to be done when no length */
2315   if ((*ch) == '[')
2316     (*initnode)(p, grbg, q, len, nodei, ntips,
2317                   parens, treewt, treenode, nodep, str, ch, treefile);
2318                              /* ... for processing a tree weight */
2319   else if ((*ch) == ';')     /* ... and at end of tree */
2320     (*initnode)(p, grbg, q, len, nodei, ntips,
2321                   parens, unittrwt, treenode, nodep, str, ch, treefile);
2322 }  /* addelement */
2323
2324
2325 void treeread (FILE *treefile, node **root, pointarray treenode,
2326         boolean *goteof, boolean *first, pointarray nodep,
2327         long *nextnode, boolean *haslengths, node **grbg, initptr initnode,
2328         boolean unifok,long maxnodes)
2329 {
2330   /* read in user-defined tree and set it up */
2331   char  ch;
2332   long parens = 0;
2333   long ntips = 0;
2334
2335   (*goteof) = false;
2336   (*nextnode) = spp;
2337
2338   /* eat blank lines */
2339   while (eoln(treefile) && !eoff(treefile))
2340     scan_eoln(treefile);
2341
2342   if (eoff(treefile)) {
2343     (*goteof) = true;
2344     return;
2345   }
2346
2347   getch(&ch, &parens, treefile);
2348
2349   while (ch != '(') {
2350     /* Eat everything in the file (i.e. digits, tabs) until you
2351        encounter an open-paren */
2352     getch(&ch, &parens, treefile);
2353   }
2354   (*haslengths) = true;
2355   addelement(root, NULL, &ch, &parens, treefile,
2356          treenode, goteof, first, nodep, nextnode, &ntips,
2357          haslengths, grbg, initnode,unifok,maxnodes);
2358
2359   /* Eat blank lines and end of current line*/
2360   do {
2361     scan_eoln(treefile);
2362   }
2363   while (eoln(treefile) && !eoff(treefile));
2364
2365   (*first) = false;
2366   if (parens != 0) {
2367     printf("\n\nERROR in tree file: unmatched parentheses\n\n");
2368     exxit(-1);
2369   }
2370 }  /* treeread */
2371
2372
2373 void addelement2(node *q, Char *ch, long *parens, FILE *treefile,
2374         pointarray treenode, boolean lngths, double *trweight, boolean *goteof,
2375         long *nextnode, long *ntips, long no_species, boolean *haslengths,
2376         boolean unifok,long maxnodes)
2377 {
2378   /* recursive procedure adds nodes to user-defined tree
2379      -- old-style bifurcating-only version */
2380   node *pfirst = NULL, *p;
2381   long i, len, current_loop_index;
2382   boolean notlast, minusread;
2383   Char str[MAXNCH];
2384   double valyew, divisor;
2385   long furs = 0;
2386
2387   if ((*ch) == '(') {
2388
2389     current_loop_index = (*nextnode) + spp;
2390     (*nextnode)++;
2391
2392     if ( maxnodes != -1 && current_loop_index > maxnodes) {
2393       printf("ERROR in intree file: Attempting to allocate too many nodes\n");
2394       printf("This is usually caused by a unifurcation.  To use this\n");
2395       printf("intree with this program  use retree to read and write\n");
2396       printf("this tree.\n");
2397       exxit(-1);
2398     }
2399     /* This is an assignment of an interior node */
2400     p = treenode[current_loop_index];
2401     pfirst = p;
2402     notlast = true;
2403     while (notlast) {
2404       furs++;
2405       /* This while loop goes through a circle (triad for
2406       bifurcations) of nodes */
2407       p = p->next;
2408       /* added to ensure that non base nodes in loops have indices */
2409       p->index = current_loop_index + 1;
2410
2411       getch(ch, parens, treefile);
2412
2413       addelement2(p, ch, parens, treefile, treenode, lngths, trweight,
2414         goteof, nextnode, ntips, no_species, haslengths,unifok,maxnodes);
2415
2416       if ((*ch) == ')') {
2417         notlast = false;
2418         do {
2419           getch(ch, parens, treefile);
2420         } while ((*ch) != ',' && (*ch) != ')' &&
2421            (*ch) != '[' && (*ch) != ';' && (*ch) != ':');
2422       }
2423     }
2424     if ( furs <= 1 && !unifok ) {
2425       printf("ERROR in intree file: A Unifurcation was detetected.\n");
2426       printf("To use this intree with this program use retree to read and");
2427       printf(" write this tree\n");
2428       exxit(-1);
2429     }
2430
2431   } else if ((*ch) != ')') {
2432     for (i = 0; i < MAXNCH; i++)
2433       str[i] = '\0';
2434     len = take_name_from_tree (ch, str, treefile);
2435     match_names_to_data (str, treenode, &p, spp);
2436     pfirst = p;
2437     if ((*ch) == ')')
2438       (*parens)--;
2439     (*ntips)++;
2440     strncpy (p->nayme, str, len);
2441   } else
2442     getch(ch, parens, treefile);
2443
2444   if ((*ch) == '[') {    /* getting tree weight from last comment field */
2445     if (!eoln(treefile)) {
2446       fscanf(treefile, "%lf", trweight);
2447       getch(ch, parens, treefile);
2448       if (*ch != ']') {
2449         printf("\n\nERROR: Missing right square bracket\n\n");
2450         exxit(-1);
2451       }
2452       else {
2453         getch(ch, parens, treefile);
2454         if (*ch != ';') {
2455           printf("\n\nERROR: Missing semicolon after square brackets\n\n");
2456           exxit(-1);
2457         }
2458       }
2459     }
2460   }
2461   else if ((*ch) == ';') {
2462     (*trweight) = 1.0 ;
2463     if (!eoln(treefile))
2464       printf("WARNING: tree weight set to 1.0\n");
2465   }
2466   else
2467     (*haslengths) = ((*haslengths) && q == NULL);
2468
2469   if (q != NULL)
2470     hookup(q, pfirst);
2471
2472   if ((*ch) == ':') {
2473     processlength(&valyew, &divisor, ch,
2474        &minusread, treefile, parens);
2475     if (q != NULL) {
2476       if (!minusread)
2477         q->oldlen = valyew / divisor;
2478       else
2479         q->oldlen = 0.0;
2480       if (lngths) {
2481         q->v = valyew / divisor;
2482         q->back->v = q->v;
2483         q->iter = false;
2484         q->back->iter = false;
2485         q->back->iter = false;
2486     }
2487     }
2488   }
2489
2490 }  /* addelement2 */
2491
2492
2493 void treeread2 (FILE *treefile, node **root, pointarray treenode,
2494         boolean lngths, double *trweight, boolean *goteof,
2495         boolean *haslengths, long *no_species,boolean unifok,long maxnodes)
2496 {
2497   /* read in user-defined tree and set it up
2498      -- old-style bifurcating-only version */
2499   char  ch;
2500   long parens = 0;
2501   long ntips = 0;
2502   long nextnode;
2503
2504   (*goteof) = false;
2505   nextnode = 0;
2506
2507   /* Eats all blank lines at start of file */
2508   while (eoln(treefile) && !eoff(treefile))
2509     scan_eoln(treefile);
2510
2511   if (eoff(treefile)) {
2512     (*goteof) = true;
2513     return;
2514   }
2515
2516   getch(&ch, &parens, treefile);
2517
2518   while (ch != '(') {
2519     /* Eat everything in the file (i.e. digits, tabs) until you
2520        encounter an open-paren */
2521     getch(&ch, &parens, treefile);
2522   }
2523
2524   addelement2(NULL, &ch, &parens, treefile, treenode, lngths, trweight,
2525           goteof, &nextnode, &ntips, (*no_species), haslengths,unifok,maxnodes);
2526   (*root) = treenode[*no_species];
2527
2528   /*eat blank lines */
2529   while (eoln(treefile) && !eoff(treefile))
2530     scan_eoln(treefile);
2531
2532   (*root)->oldlen = 0.0;
2533
2534   if (parens != 0) {
2535     printf("\n\nERROR in tree file:  unmatched parentheses\n\n");
2536     exxit(-1);
2537   }
2538 }  /* treeread2 */
2539
2540
2541 void exxit(int exitcode)
2542 {
2543 #ifdef WIN32
2544   if (exitcode == 0)
2545 #endif
2546     exit (exitcode);
2547 #ifdef WIN32
2548   else {
2549     puts ("Hit Enter or Return to close program.");
2550     puts("  You may have to hit Enter or Return twice.");
2551     getchar ();
2552     getchar ();
2553     phyRestoreConsoleAttributes();
2554     exit (exitcode);
2555   }
2556 #endif
2557 } /* exxit */
2558
2559
2560 char gettc(FILE* file)
2561 { /* catch eof's so that other functions not expecting an eof
2562    * won't have to worry about it */
2563   int ch;
2564
2565   ch=getc(file);
2566
2567   if (ch == EOF )   {
2568     puts("Unexpected End of File");
2569     exxit(-1);
2570   }
2571
2572   if ( ch == '\r' ) {
2573     ch = getc(file);
2574     if ( ch != '\n' )
2575       ungetc(ch,file);
2576     ch = '\n';
2577   }
2578   return ch;
2579 } /* gettc */
2580
2581 void unroot(tree *t, long nonodes)
2582 {
2583   /* used by fitch, restml and contml */
2584   if (t->start->back == NULL) {
2585     if (t->start->next->back->tip)
2586       t->start = t->start->next->next->back;
2587     else  t->start = t->start->next->back;
2588   }
2589   if (t->start->next->back == NULL) {
2590     if (t->start->back->tip)
2591       t->start = t->start->next->next->back;
2592     else t->start = t->start->back;
2593   }
2594   if (t->start->next->next->back == NULL)  {
2595     if (t->start->back->tip)
2596       t->start = t->start->next->back;
2597     else t->start = t->start->back;
2598   }
2599
2600
2601   unroot_r(t->start,t->nodep,nonodes);
2602   unroot_r(t->start->back, t->nodep, nonodes);
2603 }
2604
2605 void unroot_here(node* root, node** nodep, long nonodes)
2606 {
2607   node* tmpnode;
2608   double newl;
2609   /* used by unroot */
2610   /* assumes bifurcation this is ok in the programs that use it */
2611
2612
2613   newl = root->next->oldlen + root->next->next->oldlen;
2614   root->next->back->oldlen = newl;
2615   root->next->next->back->oldlen = newl;
2616
2617   newl = root->next->v + root->next->next->v;
2618   root->next->back->v = newl;
2619   root->next->next->back->v = newl;
2620
2621   root->next->back->back=root->next->next->back;
2622   root->next->next->back->back = root->next->back;
2623
2624   while ( root->index != nonodes ) {
2625     tmpnode = nodep[ root->index ];
2626     nodep[root->index] = root;
2627     root->index++;
2628     root->next->index++;
2629     root->next->next->index++;
2630     nodep[root->index - 2] = tmpnode;
2631     tmpnode->index--;
2632     tmpnode->next->index--;
2633     tmpnode->next->next->index--;
2634   }
2635 }
2636
2637 void unroot_r(node* p, node** nodep, long nonodes)
2638 {
2639   /* used by unroot */
2640   node *q;
2641
2642   if ( p->tip) return;
2643
2644   q = p->next;
2645   while ( q != p ) {
2646     if (q->back == NULL)
2647       unroot_here(q,nodep,nonodes);
2648     else unroot_r(q->back,nodep,nonodes);
2649     q = q->next;
2650   }
2651 }
2652
2653 void clear_connections(tree *t, long nonodes)
2654 {
2655   long i;
2656   for ( i = 0 ; i < nonodes ; i++) {
2657     if ( i > spp) {
2658       t->nodep[i]->next->back = NULL;
2659       t->nodep[i]->next->v = 0;
2660       t->nodep[i]->next->next->back = NULL;
2661       t->nodep[i]->next->next->v = 0;
2662     }
2663     t->nodep[i]->back = NULL;
2664     t->nodep[i]->v = 0;
2665   }
2666 }
2667
2668 #ifdef WIN32
2669 void phySaveConsoleAttributes()
2670 {
2671   GetConsoleScreenBufferInfo( hConsoleOutput, &savecsbi );
2672 } /* PhySaveConsoleAttributes */
2673
2674
2675 void phySetConsoleAttributes()
2676 {
2677   hConsoleOutput = GetStdHandle(STD_OUTPUT_HANDLE);
2678
2679   phySaveConsoleAttributes();
2680
2681   SetConsoleTextAttribute(hConsoleOutput,
2682     BACKGROUND_GREEN | BACKGROUND_BLUE | BACKGROUND_INTENSITY);
2683 } /* phySetConsoleAttributes */
2684
2685
2686 void phyRestoreConsoleAttributes()
2687 {
2688   COORD coordScreen = { 0, 0 };
2689   DWORD cCharsWritten;
2690   DWORD dwConSize;
2691
2692   dwConSize = savecsbi.dwSize.X * savecsbi.dwSize.Y;
2693
2694   SetConsoleTextAttribute(hConsoleOutput, savecsbi.wAttributes);
2695
2696   FillConsoleOutputAttribute( hConsoleOutput, savecsbi.wAttributes,
2697          dwConSize, coordScreen, &cCharsWritten );
2698 } /* phyRestoreConsoleAttributes */
2699
2700
2701 void phyFillScreenColor()
2702 {
2703   COORD coordScreen = { 0, 0 };
2704   DWORD cCharsWritten;
2705   CONSOLE_SCREEN_BUFFER_INFO csbi; /* to get buffer info */
2706   DWORD dwConSize;
2707
2708   GetConsoleScreenBufferInfo( hConsoleOutput, &csbi );
2709   dwConSize = csbi.dwSize.X * csbi.dwSize.Y;
2710
2711   FillConsoleOutputAttribute( hConsoleOutput, csbi.wAttributes,
2712          dwConSize, coordScreen, &cCharsWritten );
2713 } /* PhyFillScreenColor */
2714
2715
2716 void phyClearScreen()
2717 {
2718    COORD coordScreen = { 0, 0 };    /* here's where we'll home the
2719                                        cursor */
2720    DWORD cCharsWritten;
2721    CONSOLE_SCREEN_BUFFER_INFO csbi; /* to get buffer info */
2722    DWORD dwConSize;                 /* number of character cells in
2723                                        the current buffer */
2724
2725    /* get the number of character cells in the current buffer */
2726
2727    GetConsoleScreenBufferInfo( hConsoleOutput, &csbi );
2728    dwConSize = csbi.dwSize.X * csbi.dwSize.Y;
2729
2730    /* fill the entire screen with blanks */
2731
2732    FillConsoleOutputCharacter( hConsoleOutput, (TCHAR) ' ',
2733       dwConSize, coordScreen, &cCharsWritten );
2734
2735    /* get the current text attribute */
2736
2737    GetConsoleScreenBufferInfo( hConsoleOutput, &csbi );
2738
2739    /* now set the buffer's attributes accordingly */
2740
2741    FillConsoleOutputAttribute( hConsoleOutput, csbi.wAttributes,
2742       dwConSize, coordScreen, &cCharsWritten );
2743
2744    /* put the cursor at (0, 0) */
2745
2746    SetConsoleCursorPosition( hConsoleOutput, coordScreen );
2747    return;
2748 } /* PhyClearScreen */
2749 #endif
2750