phylip.c

一个神经网络工具箱

/* version 3.6. (c) Copyright 1993-2002 by the University of Washington.
   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, Andrew Keeffe,
   and Dan Fineman.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */

#include <stdio.h>
#include <signal.h>
#ifdef WIN32
#include <windows.h>
/* for console code (clear screen, text color settings) */
CONSOLE_SCREEN_BUFFER_INFO savecsbi;
HANDLE hConsoleOutput;

void phyClearScreen();
void phySaveConsoleAttributes();
void phySetConsoleAttributes();
void phyRestoreConsoleAttributes();
void phyFillScreenColor();
#endif

#include "phylip.h"

#ifndef OLDC
static void crash_handler(int signum);

#endif

FILE *infile, *outfile, *intree, *intree2, *outtree, *weightfile, *catfile, *ancfile, *mixfile, *factfile;
long spp, words, bits;
boolean ibmpc, ansi, tranvsp;
naym *nayme;                     /* names of species */

static void crash_handler(int sig_num)
{ /* when we crash, lets print out something usefull */
  printf("ERROR:  ");
  switch(sig_num) {
#ifdef SIGSEGV
    case SIGSEGV:
      puts("This program has caused a Segmentation fault.");
      break;
#endif /* SIGSEGV */
#ifdef SIGFPE
    case SIGFPE:
      puts("This program has caused a Floating Point Exception");
      break;
#endif  /* SIGFPE */
#ifdef SIGILL
    case SIGILL:
      puts("This program has attempted an illegal instruction");
      break;
#endif  /* SIGILL */
#ifdef SIGPIPE 
    case SIGPIPE:
      puts("This program tried to write to a broken pipe");
      break;
#endif  /* SIGPIPE */
#ifdef SIGBUS
    case SIGBUS:
      puts("This program had a bus error");
      break;
#endif /* SIGBUS */
  }   
  if (sig_num == SIGSEGV) {
    puts("       This may have been caused by an incorrectly formatted");
    puts("       input tree file or input file.  You should check those");
    puts("       files carefully.");
    puts("       If this seems to be a bug, please mail [e-mail]");
  }
  else {
    puts("       Most likely, you have encountered a bug in the program.");
    puts("       Since this seems to be a bug, please mail [e-mail]");
  }
  puts("       with the name of the program, your computer system type,");
  puts("       a full description of the problem, and with the input data file.");
  puts("       (which should be in the body of the message, not as an Attachment).");

#ifdef WIN32
  puts ("Hit Enter or Return to close program.");
  puts("  You may have to hit Enter or Return twice.");
  getchar ();
  getchar ();
  phyRestoreConsoleAttributes();
#endif
  abort();
}


void init(int argc, char** argv) 
{ /* initialization routine for all programs 
   * anything done at the beginig for every program should be done here */ 
 
  /* set up signal handler for 
   * segfault,floating point exception, illeagal instruction, bad pipe, bus error
   * there are more signals that can cause a crash, but these are the most common
   * even these aren't found on all machines.  */
#ifdef SIGSEGV
  signal(SIGSEGV, crash_handler);
#endif /* SIGSEGV */
#ifdef SIGFPE
  signal(SIGFPE, crash_handler);
#endif /* SIGFPE */
#ifdef SIGILL
  signal(SIGILL, crash_handler);
#endif /* SIGILL */
#ifdef SIGPIPE
  signal(SIGPIPE, crash_handler);
#endif /* SIGPIPE */
#ifdef SIGBUS
  signal(SIGBUS, crash_handler);
#endif /* SIGBUS */

#ifdef WIN32
  phySetConsoleAttributes();
  phyClearScreen();
#endif

}

void scan_eoln(FILE *f) 
{ /* eat everything to the end of line or eof*/
  char ch;

  while (!eoff(f) && !eoln(f)) 
    gettc(f);
  if (!eoff(f)) 
    ch = gettc(f);
  if (ch == '\r' && !eoff(f) && eoln(f))
    gettc(f);
}


boolean eoff(FILE *f)
{ /* check for end of file */
    int ch;

    if (feof(f)) 
      return true;
    ch = getc(f);
    if (ch == EOF) {
      ungetc(ch, f);
      return true;
    }
    ungetc(ch, f);
    return false;
}  /*eoff*/


boolean eoln(FILE *f)
{ /* check for end of line or eof*/
    register int ch;

    ch = getc(f);
    if (ch == EOF)
      return true;
    ungetc(ch, f);
    return ((ch == '\n') || (ch == '\r'));
}  /*eoln*/


int filexists(char *filename)
{ /* check whether file already exists */
  FILE *fp;
  fp =fopen(filename,"r");
  if (fp) {
    fclose(fp);
    return 1;
  } else
    return 0;
}  /*filexists*/


const char* get_command_name (const char *vektor)
{ /* returns the name of the program from vektor without the whole path */
  char *last_slash;

  /* Point to the last slash... */
  last_slash = strrchr (vektor, DELIMITER);

  if (last_slash)
    /* If there was a last slash, return the character after it */
    return last_slash + 1;
  else
    /* If not, return the vector */
    return vektor;

}  /*get_command_name*/


void getstryng(char *fname)
{ /* read in a file name from stdin and take off newline if any */

  fname = fgets(fname, 100, stdin);
  if (strchr(fname, '\n') != NULL)
    *strchr(fname, '\n') = '\0';
} /* getstryng */


void countup(long *loopcount, long maxcount)
{ /* count how many times this loop has tried to read data, bail out
     if exceeds maxcount */

  (*loopcount)++;
  if ((*loopcount) >= maxcount) {
    printf("\nERROR: Made %ld attempts to read input in loop. Aborting run.\n",
            *loopcount);
    exxit(-1);
  }
} /* countup */


void openfile(FILE **fp,const char *filename,const char *filedesc,
              const char *mode,const char *application, char *perm)
{ /* open a file, testing whether it exists etc. */
  FILE *of;
  char file[FNMLNGTH];
  char filemode[2];
  char input[FNMLNGTH];
  char ch;
  const char *progname_without_path;
  long loopcount, loopcount2;

  progname_without_path = get_command_name(application);

  strcpy(file,filename);
  strcpy(filemode,mode);
  loopcount = 0;



	  while (1){
		if (filemode[0] == 'w' && filexists(file)){

/*
		  printf("\n%s: the file \"%s\" that you wanted to\n",
			  progname_without_path, file);
		  printf("     use as %s already exists.\n", filedesc);
		  printf("     Do you want to Replace it, Append to it,\n");
		  printf("     write to a new File, or Quit?\n");
		  loopcount2 = 0;
*/

	/*
		  do {
			printf("     (please type R, A, F, or Q) \n");
	#ifdef WIN32
			phyFillScreenColor();
	#endif
			fgets(input, sizeof(input), stdin);
			ch  = input[0];
			uppercase(&ch);
			countup(&loopcount2, 10);
		  } while (ch != 'A' && ch != 'R' && ch != 'F' && ch != 'Q');

	*/

ch = 'R';	// always replace file

	  if (ch == 'Q')
        exxit(-1);
      if (ch == 'A') {
        strcpy(filemode,"a");
        continue;
      }
      else if (ch == 'F') {
        file[0] = '\0';
        loopcount2 = 0;
        while (file[0] =='\0') {
          printf("Please enter a new file name> ");
          getstryng(file);
          countup(&loopcount2, 10);
        }
        strcpy(filemode,"w");
        continue;
      }
    }
    of = fopen(file,filemode);
    if (of)
      break;
    else {
      switch (filemode[0]){

      case 'r':
        printf("%s: can't find %s \"%s\"\n", progname_without_path,
            filedesc, file);
        file[0] = '\0';
        loopcount2 = 0;
        while (file[0] =='\0'){
          printf("Please enter a new file name> ");
          countup(&loopcount2, 10);
          getstryng(file);}
        break;

      case 'w':
      case 'a':
        printf("%s: can't write %s file %s\n", progname_without_path,
            filedesc, file);
        file[0] = '\0';
        loopcount2 = 0;
        while (file[0] =='\0'){
          printf("Please enter a new file name> ");
          countup(&loopcount2, 10);
          getstryng(file);}
        continue;
      default:
     printf("There is some error in the call of openfile. Unknown mode.\n");
        exxit(-1);
      }
    }
    countup(&loopcount, 20);
  }
  *fp = of;
  if (perm != NULL)
    strcpy(perm,file);
} /* openfile */


void cleerhome()
{ /* home cursor and clear screen, if possible */
  printf("%s", ((ibmpc || ansi) ? ("\033[2J\033[H") : "\n\n"));
} /* cleerhome */


double randum(longer seed)
{ /* random number generator -- slow but machine independent
  This is a multiplicative congruential 32-bit generator
  x(t+1) = 1664525 * x(t) mod 2^32, one that passes the
  Coveyou-Macpherson and Lehmer tests, see Knuth ACP vol. 2   */
  long i, j, k, sum;
  longer mult, newseed;
  double x;

  mult[0] = 13;   /* these four statements set the multiplier */
  mult[1] = 24;   /* -- they are its "digits" in a base-64    */
  mult[2] = 22;   /*    notation: 1664525 = 13*64^3+24*64^2   */
  mult[3] = 6;    /*                         +22*64+6         */
  for (i = 0; i <= 5; i++)
    newseed[i] = 0;
  for (i = 0; i <= 5; i++) {
    sum = newseed[i];
    k = i;
    if (i > 3)
      k = 3;
    for (j = 0; j <= k; j++)
      sum += mult[j] * seed[i - j];
    newseed[i] = sum;
    for (j = i; j <= 4; j++) {
      newseed[j + 1] += newseed[j] / 64;
      newseed[j] &= 63;
    }
  }
  memcpy(seed, newseed, sizeof(longer));
  seed[5] &= 3;
  x = 0.0;
  for (i = 0; i <= 5; i++)
    x = x / 64.0 + seed[i];
  x /= 4.0;
  return x;
}  /* randum */


void randumize(longer seed, long *enterorder)
{ /* randomize input order of species */
  long i, j, k;

  for (i = 0; i < spp; i++) {
    j = (long)(randum(seed) * (i+1));
    k = enterorder[j];
    enterorder[j] = enterorder[i];
    enterorder[i] = k;
  }
} /* randumize */


double normrand(longer seed)
{/* standardized Normal random variate */
  double x;

  x = randum(seed)+randum(seed)+randum(seed)+randum(seed)
       + randum(seed)+randum(seed)+randum(seed)+randum(seed)
       + randum(seed)+randum(seed)+randum(seed)+randum(seed)-6.0;
  return(x);
} /* normrand */ 


long readlong(const char *prompt)
{ /* read a long */
  long res, loopcount;
  char string[100];

  loopcount = 0;
  do {
    printf("%s",prompt);
    getstryng(string);
    if (sscanf(string,"%ld",&res) == 1)
      break;
    countup(&loopcount, 10);
   } while (1);
  return res;
}  /* readlong */


void uppercase(Char *ch)
{ /* convert ch to upper case */
  *ch = (islower (*ch) ? toupper(*ch) : (*ch));
}  /* uppercase */


void initseed(long *inseed, long *inseed0, longer seed)
{ /* input random number seed */
  long i, loopcount;

  loopcount = 0;
  do {
    printf("Random number seed (must be odd)?\n");
    scanf("%ld%*[^\n]", inseed);
    getchar();
    countup(&loopcount, 10);
  } while (((*inseed) < 0) || ((*inseed) & 1) == 0);
  *inseed0 = *inseed;
  for (i = 0; i <= 5; i++)
    seed[i] = 0;
  i = 0;
  do {
    seed[i] = *inseed & 63;
    *inseed /= 64;
    i++;
  } while (*inseed != 0);
}  /*initseed*/


void initjumble(long *inseed, long *inseed0, longer seed, long *njumble)
{ /* input number of jumblings for jumble option */
  long loopcount;

  initseed(inseed, inseed0, seed);
  loopcount = 0;
  do {
    printf("Number of times to jumble?\n");
    scanf("%ld%*[^\n]", njumble);
    getchar();
    countup(&loopcount, 10);
  } while ((*njumble) < 1);
}  /*initjumble*/


void initoutgroup(long *outgrno, long spp)
{ /* input outgroup number */
  long loopcount;
  boolean done;

  loopcount = 0;
  do {
    printf("Type number of the outgroup:\n");
    scanf("%ld%*[^\n]", outgrno);
    getchar();
    done = (*outgrno >= 1 && *outgrno <= spp);
    if (!done) {
      printf("BAD OUTGROUP NUMBER: %ld\n", *outgrno);
      printf("  Must be in range 1 - %ld\n", spp);
    }
    countup(&loopcount, 10);
  } while (done != true);
}  /*initoutgroup*/


void initthreshold(double *threshold)
{ /* input threshold for threshold parsimony option */
  long loopcount;
  boolean done;

  loopcount = 0;
  do {
    printf("What will be the threshold value?\n");
    scanf("%lf%*[^\n]", threshold);
    getchar();
    done = (*threshold >= 1.0);
    if (!done)
      printf("BAD THRESHOLD VALUE:  it must be greater than 1\n");
    else
      *threshold = (long)(*threshold * 10.0 + 0.5) / 10.0;
    countup(&loopcount, 10);
  } while (done != true);
}  /*initthreshold*/


void initcatn(long *categs)
{ /* initialize category number for rate categories */
  long loopcount;

  loopcount = 0;
  do {
    printf("Number of categories (1-%d)?\n", maxcategs);
    scanf("%ld%*[^\n]", categs);
    getchar();
    countup(&loopcount, 10);
  } while (*categs > maxcategs || *categs < 1);
}  /*initcatn*/


void initcategs(long categs, double *rate)
{ /* initialize category rates for HMM rates */
  long i, loopcount, scanned;
  char line[100], rest[100];
  boolean done;

  loopcount = 0;
  for (;;){
    printf("Rate for each category? (use a space to separate)\n");
    getstryng(line);