oblong/oblong.c

/*************************************************************************************
 *  OBLONG
 *      a program for finding most parsimonious 
 *      trees with minimum RAM requirements
 *
 *  Copyright April 2013 by Pablo A. Goloboff
 *
 *  This program is free software; you may redistribute it and/or modify its
 *  under the terms of the GNU General Public License as published by the Free
 *  Software Foundation; either version 2 of the License, or (at your option)
 *  any later version.
 *
 *  This program is distributed in the hope that it will be useful, but
 *  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 *  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 *  for more details.
 *
 *  If publishing results based on this program, please cite Goloboff (forthcoming)
 *
 ***********************************************************************************
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#include <ctype.h>
#include <math.h>
#define MAXNAMELEN      100
#define CHUNKSZ         1000
#define MAXMULTITAXA    1000 

FILE * inp = NULL , * outp = NULL , * tmpp , * partsfile , * treeinp = NULL , * multinp ;

#ifdef WATCOM
long long * fpos ;
#else 
fpos_t * fpos ;
#endif

unsigned long long bestgloballen , lenlag = 0 , rearrangs = 0 ; 
unsigned long int ** mat , * wtchg , * movelist , lenlimit ;
unsigned long int gup , gmix , gdn , nitdn ;
long int nc2 = 0 , nc3 = 0 , nc4 = 0 , trashed = 0 , totnc , maxnc  , piecesz , totnc4 ;
double * gbrem , * tbrem ;
long int * abrem , fasta = 0 ; 
int nt , nc , nodz , nstord , together = 1 , use_spr = 0 , reportreplication = 0 ,
    numreplications = 0 , reportperiod = 1 , ratcycles = 0 , ratcyclesdone = 0 ,
    dojak = 0 , saveall = 0 , swapcollapse = 0 , recordgoodmoves = 0 , usenexus = 0 ;
int ranseed = 1 , do_branch_swapping = 1 , random_starts = 0 , numtrees = 0 , numhits ,
    ratcheting = 0 , changewts = 0 , numsubs , jakprob = 36 , numcolmoves , numcolroots , 
    savejusttrees = 0 , felipe = 0 , multifiles = 0 , nummultifiles , recorddn = 0 ,   
    equal_taxon_sequence = 0 , numparts , quickbremer = 0 , doingbremer = 0 ,
    bremertype = 2 , shownegbremers = 0 , verbose = 0 ;  
short int * nodlist , * anx , * lef , * rig , * sis , * rlist ,
          * tgtlist , * lassis , * nodfork , nodtoskip = 0 , * ladd , * nodmatch , 
          * blist , ** memtrees , * unsuprt , * distoroot , * colmoves , * colroots ;
char strsp[ MAXNAMELEN ] , maskin[ 256 ] , bitct [ 65536 ] , reportheader ,
     reporttoscreen = 1 ;
char legindic [ ] = { '=' , 0 } ; 
unsigned char ** matbuf , * tmpstates , * isabs = NULL ;
double maxram = -1 ; 

typedef struct { unsigned long int up , dn ; int root , tgt ; } BremTyp ;
BremTyp * bremdif ; 

void transfer_bremer ( int , int , BremTyp * ) ;

unsigned long long downpass_all ( void ) ;
unsigned long long downpass_split ( void ) ;
unsigned long long ( * downpass ) ( void ) ; 

unsigned long long wagner_all ( void ) ;
unsigned long long wagner_split ( void ) ;
unsigned long long ( * wagner ) ( void ) ; 

unsigned long long bswap_all ( void ) ;
unsigned long long bswap_split ( void ) ;
unsigned long long ( * bswap ) ( void ) ; 

#define VERYBIG 10000000000000000 

#ifdef BITS64
#define FLDS2 32 
#define FLDS3 21 
#define FLDS4 16
#define HI4   9838263505978427528UL
#define LO4   8608480567731124087UL
#define HI3   5270498306774157604UL
#define LO3   3952873730080618203UL
#define HI2   12297829382473034410UL
#define LO2   6148914691236517205UL
#else 
#define FLDS2 16 
#define FLDS3 10 
#define FLDS4 8
#define HI4  2290649224
#define LO4  2004318071
#define HI3  613566756
#define LO3  460175067
#define HI2  2863311530
#define LO2  1431655765
#endif

/*********** TRIVIAL BUT USEFUL FUNCTIONS  ******************/

void display_help ( int kind )
{
    if ( kind == 0 ) 
    fprintf ( stdout ,
    "\n"
    "  OBLONG (Merriam-Webster): \"a rectangle with one of its\n"
    "  sides much longer than the other\"\n"
    "\n"
    "                 by P.A. Goloboff (2013)\n"
    "\n" 
    "     For phylogenetic analysis of wide matrices, with low memory\n"
    "  requirements.  Recommended for data sets with many characters\n"
    "  and few taxa (not advised for hundreds of taxa or more, as\n"
    "  it is optimized for memory, not speed)\n" 
    "     Uses temporary disk buffers, \"oblong.tmp\" and/or \"oblong.pts\"\n"
    "  to avoid holding all data in RAM. As much disk space as the size\n"
    "  of the data set may be required.\n" 
    "     Screen output can be redirected to a file with \"2> file\"\n"
    "\n"
    "  Use \"oblong -h\" for help on program options\n"
    "\n"
    ) ; 
    else
    fprintf ( stdout , 
    "\n"
    " Input/Output:\n"
    "\n"
    "    -ifilename   use \"filename\" as data file.\n"
    "    -p           data files in Phylip/RAxML instead of TNT format.\n"
    "    -f           data files in Fasta instead of TNT format.\n"
    "    -@list       get from file \"list\" the list of data files, must\n"
    "                 be in either Phylip (no \"-p\" switch required), or\n"
    "                 Fasta format (\"-f\" switch required). By default, \n"
    "                 taxon names are matched with identical strings; a \n"
    "                 taxon absent from a data set gets all \"?\".\n"
    "    -e           for file lists, instead of matching taxa based on name,\n"
    "                 assume all data files have the same taxon sequence.\n" 
    "    -ooutfile    output results to \"outfile\"; this is a TNT-ready\n"
    "                 data set with no characters (just tree topologies).\n"
    "\n"
    " Searches:\n" 
    "\n"
    "    -w           use random trees as starting points, not Wagner trees\n"
    "    -rN          use N starting points.\n"
    "    -spr         use SPR branch swapping instead of the default TBR.\n"
    "    -b           no branch swapping.\n"
    "    -a           save the tree for each replication (even if suboptimal)\n"
    "    -jN          jackknife the data, with N = P(del) (no N = 36). A \n"
    "                 single (possibly collapsed) tree saved per replicate.\n"
    "                 Note: uninformative characters/states are never \n" 
    "                 deactivated (they don't influence supports).\n" 
    "    -mX          use no more than X MBytes of RAM. If not using \"-j\" or\n"
    "                 \"-x\", results are identical to the default but slower;\n"
    "                 otherwise, randomization sequences may differ.\n"
    "    -RN          use N as random seed.\n"
    "    -xN          use N iterations of ratchet (no N = 10). Note that for\n"
    "                 C informative characters, \"-x\" requires C/2 additional\n"
    "                 bytes of RAM (even under \"-m\").\n"
    "    -k           if jackknifing, do TBR/SPR tree collapsing (recommended).\n"
    "\n"
    " Other:\n" 
    "\n"
    "    -N           output Nexus format (i.e. tree-file is TreeView-ready)\n" 
    "    -t           save only the trees (no \"xread\").\n"
    "    -Ttreefile   read trees from \"treefile\" and use them as starting\n"
    "                 points (if \"-b\" and no \"-x\", it just calculates\n"
    "                 tree-scores).  Treefile must be the output of Oblong\n"
    "                 (i.e. what is in \"outfile\"), with number of trees in\n"
    "                 tree-title.\n"
    "    -q           use a quick approximmation for group supports. Default is\n"
    "                 S = 0.64 * R ^ 1/A (where R = relative bremer support,\n"
    "                 and A = absolute bremer support); this is also indicated\n" 
    "                 with \"-q2\".  You can also indicate \"-q0\", absolute bremer\n" 
    "                 support, or \"-q1\", relative. Supports are obtained via\n" 
    "                 branch-swapping, recording score difference for each move.\n"
    "                 You can use this to measure strength of contradiction, starting\n"
    "                 from suboptimal trees (\"-T\") and skipping swapping (\"-b\").\n"
    "    -v           verbose mode: report rearrangements examined as they are done\n"
    "                 (note: repeatedly writing to screen slows down program).\n" 
    "\n"
    " You can also view this info with \"oblong -h | less\" or \"oblong -h | more\"\n"
    "\n"  
    ) ; 
}    

void errout ( int cond , char * msg )
{
    if ( !cond ) return ;
    if ( msg != NULL ) 
       fprintf ( stderr , "\n%s\n" , msg ) ; //output to the screen, usually error
    exit ( 0 ) ; //exit(0)正常运行程序并退出程序；exit(1)非正常运行导致退出程序
}

void * myalloc ( unsigned long int size )
{
    void * vd = malloc ( size ) ;
    errout ( vd == NULL , "Sorry, not enough memory available --try with the \"-m\" switch" ) ;
    return vd ;
}

int stringis ( char * a , char * b )
{
    while ( * a && tolower ( * a ) == tolower ( * b ) ) { ++a ; ++ b ; } //tolower: conver letters from uppercase to lowercase
    if ( * a || * b ) return 0 ; // or || ,return fake
    return 1 ; //return true
}
    
void procargs ( int nargs , char ** arg ) // int nargs =2 for [oblong -f] **arg = *arg [] 字符串数组
{
    int i ;
    char * cp ;
    if ( nargs == 1 ) { // no arguments if nargs == 1
        display_help ( 0 ) ; //not show the help ouput just no argument output
        exit ( 0 ) ; }
    for ( i = 1 ; i < nargs ; ++ i ) {
        cp = arg [ i ] ; 
        if ( * cp  != '-' ) { // if arguments not start with -
            fprintf ( stderr , "\nWrong argument: %s\n" , arg [ i ] ) ; //show which argument is wrong
            exit ( 0 ) ; }
        switch ( * ++ cp ) {
            case 'a': saveall = 1 ; break ; //save all trees
            case 'b': do_branch_swapping = 0 ; break ; //no branch swapping
            case 'e': equal_taxon_sequence = 1 ; break ; 
            case 'f': fasta = 1 ; break ; //fasta format
            case 'h': display_help ( 1 ) ; exit ( 0 ) ; //help
            case 'i': //use the filename -i${filename}
            case '@':
                if ( * cp == '@' ) multifiles = 1 ; //there are multi files
                ++ cp ; 
                errout ( stringis ( cp , "oblong.tmp" ) , "Cannot use \"oblong.tmp\" as input file; rename it and try again" ) ; 
                errout ( stringis ( cp , "oblong.pts" ) , "Cannot use \"oblong.pts\" as input file; rename it and try again" ) ; 
                errout ( inp != NULL , "Can't open TWO input files" ) ;
                errout ( ( inp = fopen ( cp , "rb" ) ) == NULL , "Can't open input file" ) ; //input file, if the input list can't be readable in the bianry mode, on windows, binary is different
                break ;
            case 'j': //j${NumberofJakknife}
               ++ cp ;
               if ( isdigit ( * cp ) ) //to find if the character is the numer
                    jakprob = atoi ( cp ) ; //atoi: convert character to the int
               dojak = changewts = 1 ;
               break ; 
            case 'k': swapcollapse = 1 ; break ; //do TBR/SPR collapsing
	    case 'm': errout ( ( ( maxram = atof ( ++ cp ) ) <= 0 ) , "You must indicate RAM to use" ) ; together = 0 ; break ; //atof: transfer chracter into double
            case 'N': usenexus = 1 ; break ; //output the nexus file
            case 'o': //-o${outfile}
                ++ cp ; 
                errout ( stringis ( cp , "oblong.tmp" ) , "Cannot use \"oblong.tmp\" as output file; use a different name" ) ; 
                errout ( stringis ( cp , "oblong.pts" ) , "Cannot use \"oblong.pts\" as output file; use a different name" ) ; 
                errout ( outp != NULL , "Can't open TWO output files" ) ;
                errout ( ( outp = fopen ( cp , "wb" ) ) == NULL , "Can't open output file" ) ;
                break ;
            case 'p': felipe = 1 ; break ; //file in Phylip style
            case 'q': //a quick approximmation
                ++ cp ; 
                quickbremer = 1 ;
                if ( * cp == '0' ) bremertype = 0 ; 
                if ( * cp == '1' ) bremertype = 1 ; 
                if ( * cp == '2' ) bremertype = 2 ; 
                break ; 
            case 'r': numreplications = atoi ( ++ cp ) ; break ; //use N starting points. -r${Number}
            case 'R': ranseed = atoi ( ++ cp ) ; break ; //random seed
            case 's': //spr use spr instead of tbr
                if ( !strcmp ( cp , "spr" ) ) { use_spr = 1 ; break ; }
            case 't': savejusttrees = 1 ; break ; //output standrad tre format
            case 'T': //set start tree -T${treefile}
                ++ cp ; 
                errout ( treeinp != NULL , "Can't open TWO tree files" ) ;
                errout ( ( treeinp = fopen ( cp , "rb" ) ) == NULL , "Can't open tree file" ) ;
                break ;
            case 'v': verbose = 1 ; break ; //report rearrangements
            case 'w': random_starts = 1 ; break ; //use random trees as start tree
            case 'x': //-x${Numer of ratcher}
                ++ cp ;
                if ( isdigit ( * cp ) ) ratcycles = atoi ( cp ) ;//if the char is number, make ratecycles=
                else ratcycles = 10 ;
                break ;
            default: fprintf ( stderr , "\nUnrecognized option: %s\n" , arg [ i ] ) ; exit ( 0 ) ; }}
    errout ( ( fasta && felipe ) , "Cannot use both Phylip and Fasta formats" ) ; //fasta = 1 and felipe = 1, it can't be identified as different type 
    errout ( ( inp == NULL ) , "No data set specified" ) ; 
    errout ( ( outp == NULL ) , "No output file specified" ) ;
    if ( quickbremer && dojak ) {
        dojak = changewts = 0 ; //if quick approximmation and jack both are set, the jack won't be done 
        fprintf ( stderr , "\nNOTE: Cannot do both jackknifing and bremer at the same time!\n" ) ;
        if ( bremertype == 2 )
             fprintf ( stderr , "      deletion prob P=%.2f will be used to rescale supports (1-P)\n" , ( double ) jakprob / 100 ) ;
        else
          if ( bremertype == 0 ) fprintf ( stderr , "    will do bremer supports instead\n" ) ;
          else                   fprintf ( stderr , "    will do relative bremer supports instead\n" ) ;
        errout ( ( !jakprob && bremertype == 2 ) , "    ...but cannot use a deletion prob. P=0 to rescale!" ) ;  
        fprintf ( stderr , "\n" ) ; }
    if ( quickbremer && usenexus ) legindic [ 0 ] = 0 ; // *** left here temp 
    if ( !do_branch_swapping && !ratcycles ) shownegbremers = 1 ; // *** left here temp  
    if ( treeinp != NULL && numreplications )
       fprintf ( stderr , "\nNote: \"-T\" overrides \"-r\": number of trees in file determines starting points\n" ) ; 
    if ( !numreplications ) numreplications = 1 ; 
    if ( maxram > 0 ) {
      fprintf ( stderr , "\nUsing %.2f MBytes of RAM for matrix buffers\n" , maxram ) ;
        fflush ( stderr ) ; //清空缓存区 
        maxram *= 1024 * 1024 ; } // maxram = maxram * 1024 * 1024;
    reportheader = '\r' ;
    if ( !isatty( fileno( stderr ) ) ) {
        errout ( verbose , "Verbose allowed only when STDERR goes to screen" ) ; 
        reporttoscreen = 0 ; 
        reportheader = '\n' ;
        reportperiod = 600 ; 
        fprintf ( stderr , "\n\nProgram options:\n\n   oblong " ) ;
        for ( i = 1 ; i < nargs ; ++ i )
            fprintf ( stderr , "%s " , arg [ i ] ) ;
        fprintf ( stderr , "\n\n" ) ; }
    return ; 
}

void init_bit_counters ( void )
{
    int a , b , c ; // bitct[65536]
    for ( a = 0 ; a < 65536 ; ++ a ) { // a 0 -> 65535 
      for ( b = 1 , c = 0 ; b <= a ; b <<= 1 ) // <= 是 < 或 =; 
        if ( ( b & a ) ) ++ c ; // c++; 
      bitct [ a ] = c ; }
}        

void maskinit ( void )
{
    int a ; 
    for ( a = 0 ; a < 256 ; ++ a ) maskin [ a ] = 0 ; 
    maskin [ 'a' ] = maskin [ 'A' ] = maskin [ '0' ] = 1 ; 
    maskin [ 'g' ] = maskin [ 'G' ] = maskin [ '1' ] = 2 ; 
    maskin [ 'c' ] = maskin [ 'C' ] = maskin [ '2' ] = 4 ; 
    maskin [ 't' ] = maskin [ 'T' ] = maskin [ '3' ] = 8 ; 
    maskin [ 'R' ] = maskin [ 'r' ] = maskin [ 'a' ] | maskin [ 'g' ] ; 
    maskin [ 'Y' ] = maskin [ 'y' ] = maskin [ 't' ] | maskin [ 'c' ] ; 
    maskin [ 'W' ] = maskin [ 'w' ] = maskin [ 'a' ] | maskin [ 't' ] ; 
    maskin [ 'S' ] = maskin [ 's' ] = maskin [ 'c' ] | maskin [ 'g' ] ; 
    maskin [ 'M' ] = maskin [ 'm' ] = maskin [ 'a' ] | maskin [ 'c' ] ; 
    maskin [ 'K' ] = maskin [ 'k' ] = maskin [ 'g' ] | maskin [ 't' ] ; 
    maskin [ 'B' ] = maskin [ 'b' ] = maskin [ 'c' ] | maskin [ 'g' ] | maskin [ 't' ] ; 
    maskin [ 'D' ] = maskin [ 'd' ] = maskin [ 'a' ] | maskin [ 'g' ] | maskin [ 't' ] ; 
    maskin [ 'H' ] = maskin [ 'h' ] = maskin [ 'a' ] | maskin [ 'c' ] | maskin [ 't' ] ; 
    maskin [ 'V' ] = maskin [ 'v' ] = maskin [ 'a' ] | maskin [ 'c' ] | maskin [ 'g' ] ; 
    maskin [ 'N' ] = maskin [ 'n' ] = maskin [ '-' ] = maskin [ '?' ] = 1 | 2 | 4 | 8 ;
    maskin [ '.' ] = 127 ; 
    return ; 
}

/********* CHECKING AND SORTING INPUT *************************/

void prepare ( int off , unsigned char intwo )
{
    int a , b , c , totbytes ;
    char msk , shif ;
    unsigned short int * pt ;
    unsigned long int * lp , written ;
    void * vd ; 
    int trsf ; 
    char mybits = bitct [ intwo ] ; 
    char convert[ 16 ] , goesto[ 4 ] ;
    if ( mybits == 2 ) { msk = 3 ; ++ nc2 ; }
    if ( mybits == 3 ) { msk = 7 ; ++ nc3 ; }
    if ( mybits == 4 ) { msk = 15 ; ++ nc4 ; }
    for ( a = b = 1 , c = 0 ; a < 16 ; a <<= 1 , ++ c ) 
        if ( ( a & intwo ) ) {
            goesto [ c ] = b ;
            b <<= 1 ; }
        else goesto [ c ] = 0 ; 
    convert [ 0 ] = msk ; 
    for ( a = 1 ; a < 16 ; ++ a ) {
       convert [ a ] = 0 ; 
       for ( b = 1 , c = 0 ; b <= a ; b <<= 1 , ++ c )
          if ( ( a & b ) ) convert [ a ] |= goesto [ c ] ;
       if ( !convert [ a ] ) convert [ a ] = msk ; }
    for ( a = 0 ; a < nt ; ++ a ) 
       matbuf [ a ] [ off ] = convert [ matbuf [ a ] [ off ] ] ;
    fprintf ( tmpp , "%c" , mybits ) ;
    pt = ( unsigned short int * ) tmpstates ;
    for ( a = ( nt + 1 ) / 2 ; a -- > 0 ; ) 
        pt [ a ] = 0 ; 
    vd = ( void * ) tmpstates ; 
    for ( a = shif = 0 ; a < nt ; ++ a ) {
        lp = ( unsigned long int * ) pt ;
        trsf = matbuf [ a ] [ off ] ;
        trsf <<= shif ; 
        * lp |= trsf ;
        shif += mybits ;
        if ( shif >= 16 ) {
            ++ pt ;
            shif -= 16 ; }}
    totbytes = ( ( mybits * ( nt + 1 ) ) - 1 ) / 8 ;
    written = fwrite ( vd , 1 , totbytes + 1 , tmpp ) ;
    errout ( written != totbytes + 1 , "Sorry, cannot write to temporary file (disk full?)" ) ; 
}

int domin ( int a )
{
    int best = 10 ;
    int b , c , misinsome ; 
    for ( c = 16 ; c -- ; ) {
        for ( b = misinsome = 0; b < nt && !misinsome ; ++ b ) 
            if ( ! ( matbuf [ b ] [ a ] & c ) ) misinsome = 1 ;
        if ( !misinsome )
           if ( best > bitct [ c ] ) best = bitct [ c ] ; }
   return best - 1 ;          
}

int dotrueextra ( long int a , unsigned long int intwo )
{
    int st , ret = 0 , b ; 
    for ( b = 0 ; b < nt ; ++ b ) {
        st = matbuf [ b ] [ a ] ;
        if ( st == 15 ) continue ;
        if ( ( st & ~intwo ) && !( st & intwo ) ) ret ++ ; }
    return ret ; 
}

void process_chunk ( int start )
{
    long int a , b , d ;
    unsigned char c ; 
    unsigned char * cp ;
    unsigned char refval ;
    unsigned char intwo , inone ;
    nstord = CHUNKSZ ;
    if ( start + nstord > nc ) nstord = nc - start ;
    for ( d = 0 ; d < nt ; ++ d ) 
      for ( b = nstord ; b -- ; )
         matbuf [ d ] [ b ] = 0 ; 
    for ( d = 0 ; d < nt ; ++ d ) {
       cp = matbuf [ d ] ;
       if ( isabs != NULL ) 
          if ( isabs [ d ] ) {  // This can only happen when reading multiple files and taxon is absent from this matrix... 
             for ( b = 0 ; b < nstord ; ++ b ) 
                 * cp ++ = 15 ;
             continue ; }
#ifdef WATCOM
       _fseeki64 ( inp , fpos [ d ] , SEEK_SET ) ; 
#else
       fsetpos ( inp , &fpos [ d ] ) ;
#endif       
       for ( b = 0 ; b < nstord ; ++ b ) {
           fscanf ( inp , " %c" , &c ) ;
           if ( c == '[' ) {
             refval = 0 ;  
             while ( 1 ) {
               fscanf ( inp , " %c" , &c ) ; 
               if ( c == ']' ) break ; 
               refval |= maskin [ c ] ; }}
           else 
           if ( c == '.' ) {
             errout ( !d , "Cannot use \".\" for the first taxon" ) ;
             refval = matbuf [ d ] [ b ] ; }
           else refval = maskin [ c ] ; 
           * cp ++ = refval ; }
#ifdef WATCOM
       fpos [ d ] = _ftelli64 ( inp ) ; 
#else
       fgetpos ( inp , &fpos [ d ] ) ;
#endif
       }
    for ( a = 0 ; a < nstord ; ++ a ) {
        inone = intwo = 0 ;
        for ( b = 0 ; b < nt ; ++ b ) {
            if ( matbuf [ b ] [ a ] == 15 ) continue ; 
            intwo |= inone & matbuf [ b ] [ a ] ;
            inone |= matbuf [ b ] [ a ] ; }
        if ( bitct [ intwo ] > 1 ) {
            if ( inone != intwo ) 
                if ( bitct [ inone ] == 1 ) 
                    lenlag += bitct [ ( inone & ~intwo ) ] ;
                else lenlag += dotrueextra ( a , intwo ) ; 
            prepare ( a , intwo ) ; }
        else {
            lenlag += domin ( a ) ; 
            ++ trashed ; }}
}

void popmat ( int part )   // populate matrix, for split optimization
{
    unsigned long int a , b , c , shif , mybits , totbytes ;
    unsigned long int inmsk , msk , ent , inshif , val , numelems ;
    unsigned long int * lp ;
    unsigned short int * pt ;
    unsigned long int sz = piecesz * FLDS4 , numnc = 0 ;
    void * vd ; 
    for ( b = 0 ; b < piecesz ; ++ b )
       for ( a = 0 ; a < nt ; ++ a ) 
          mat [ a ] [ b ] = -1 ; 
    vd = ( void * ) tmpstates ; 
    for ( a = 0 , c = part * FLDS4 ; a < sz && c < totnc ; ++ a , ++ c ) {
       mybits = getc ( tmpp ) ;
       if ( mybits == 2 ) msk = 3 ; 
       if ( mybits == 3 ) msk = 7 ; 
       if ( mybits == 4 ) msk = 15 ; 
       ent = numnc / FLDS4 ;
       inshif = 4 * ( numnc % FLDS4 ) ;
       inmsk = ( unsigned long int ) 15 << inshif ;
       ++ numnc ;
       totbytes = ( ( mybits * ( nt + 1 ) ) - 1 ) / 8 ;
       numelems = fread ( vd , 1 , totbytes + 1 , tmpp ) ;
       errout ( numelems != totbytes + 1 , "OOOPS, unexpected error loading matrix!.\nSorry, abandoning run..." ) ;
       pt = ( unsigned short int * ) tmpstates ;
       for ( b = shif = 0 ; b < nt ; ++ b ) {
           lp = ( unsigned long int * ) pt ; 
           val = ( * lp >> shif ) & msk ;
           mat [ b ] [ ent ] &= ~inmsk ;
           mat [ b ] [ ent ] |= val << inshif ; 
           shif += mybits ;
           if ( shif >= 16 ) {
               ++ pt ;
               shif -= 16 ; }}}
}

void create_parts_file ( void )
{
    unsigned long int atnc ;
    unsigned long int written ;
    void * vd ;
    unsigned long int * at ; 
    errout ( ( partsfile = fopen ( "oblong.pts" , "wb" ) ) == NULL , "Cannot open temporary file \"oblong.pts\"" ) ;  
    for ( atnc = 0 ; atnc < nc4 ; atnc += piecesz ) {
        popmat ( atnc ) ;
        at = mat [ 0 ] ; 
        vd = ( void * ) at ; 
        written = fwrite ( vd , 1 , sizeof ( unsigned long int ) * nt * piecesz , partsfile ) ;
        errout ( written != sizeof ( unsigned long int ) * nt * piecesz , "Sorry, cannot write to temporary file (disk full?)" )  ; }
    fputc ( 32 , partsfile ) ; 
    fclose ( partsfile ) ;
    partsfile = fopen ( "oblong.pts" , "rb" ) ; 
    return ;
}

void tmpparse ( void )
{
    int atnc2 = 0 , atnc3 = 0 , atnc4 = 0 , totbytes ;
    unsigned char mybits ;
    int a , b , shif , inshif ;
    unsigned short int * pt ;
    unsigned long int * lp , inmsk , ent ;
    signed long int sizis , toalloc , subtract ;
    unsigned long int j , * segposit , msk , val ; 
    double prop , prevprop ;
    totnc = nc2 + nc3 + nc4 ; 
    nc2 = ( nc2 + ( FLDS2 - 1 ) ) / FLDS2 ;
    nc3 = ( nc3 + ( FLDS3 - 1 ) ) / FLDS3 ;  
    nc4 = ( nc4 + ( FLDS4 - 1 ) ) / FLDS4 ;  
    nc3 += nc2 ;
    nc4 += nc3 ;
    nodz = 2 * nt - 1 ;
    nodlist = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ; 
    anx = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ; 
    lef = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
    rig = ( short int * ) myalloc ( nt * sizeof ( short int ) ) ;
    rig -= nt ; 
    sis = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
    rlist = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
    tgtlist = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
    lassis = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
    nodfork = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
    ladd = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
    blist = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
    mat = ( unsigned long int ** ) myalloc ( nodz * sizeof ( long int * ) ) ;
    if ( swapcollapse || quickbremer ) {
       unsuprt = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
       distoroot = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
       colmoves = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
       colroots = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
       if ( quickbremer ) {
          nodmatch = ( short int * ) myalloc ( nt * sizeof ( short int ) ) ;
          nodmatch -= nt ; 
          abrem = ( long int * ) myalloc ( nt * sizeof ( long int ) ) ;
          abrem -= nt ; 
          gbrem = ( double * ) myalloc ( nt * sizeof ( double ) ) ;
          gbrem -= nt ; 
          tbrem = ( double * ) myalloc ( nt * sizeof ( double ) ) ;
          tbrem -= nt ; 
          bremdif = ( BremTyp * ) myalloc ( nt * nt * sizeof ( BremTyp ) ) ; }}
    b = numreplications ;
    if ( ratcycles ) b = ( numreplications * ratcycles ) + 1 ; 
    memtrees = ( short int ** ) myalloc ( b * sizeof ( short int * ) ) ;
    for ( a = 0 ; a < b ; ++ a )
       memtrees [ a ] = ( short int * ) myalloc ( nodz * sizeof ( short int ) ) ;
    toalloc = nc4 ;
    downpass = downpass_all ;
    bswap = bswap_all ;
    wagner = wagner_all ; 
    if ( !together ) {
        sizis = ( signed long int ) maxram ;
        subtract = ( 12 * nodz * sizeof ( short int ) ) +   // tree-buffers... 
               ( nt * nt * sizeof ( long int ) ) +      // buffer to store moves...
               ( b * sizeof ( short int ** ) ) + ( b * nodz * sizeof ( short int ) ) ; // trees
        subtract += nodz * sizeof ( long int * ) ;
        if ( subtract >= sizis ) {
            fprintf ( stderr , "\nSorry: %.0f Bytes are insufficient to run this data set. Increase \"-m\"\n" , maxram ) ;
            exit ( 0 ) ; }
        toalloc = ( sizis - subtract ) / ( nodz * sizeof ( long int ) ) ;
        j = ( totnc + ( FLDS4 - 1 ) ) / FLDS4 ;  // needed total 
        numparts = 1 + ( j / toalloc ) ;
        toalloc = j / numparts ;
        while ( toalloc * numparts * FLDS4 < totnc )
            toalloc ++ ; 
        if ( toalloc < 1 ) {
            fprintf ( stderr , "\nSorry: %.0f Bytes are insufficient to run this data set. Increase -m\n" , maxram ) ;
            exit ( 0 ) ; }
        if ( toalloc * FLDS4 < totnc ) {
            piecesz = toalloc ;
            nc2 = nc3 = 0 ; 
            nc4 = totnc4 = numparts * piecesz ; 
            fprintf ( stderr , "%cWill process the chars in %i chunks of %li\n" , reportheader , numparts , toalloc * FLDS4 ) ; fflush ( stderr ) ; 
            movelist = ( unsigned long int * ) myalloc ( nt * nt * sizeof ( long int ) ) ; 
            downpass = downpass_split ;
            wagner = wagner_split ;
            bswap = bswap_split ; }
        else {
            fprintf ( stderr , "%cLucky you! All data can be held in RAM! No disk swaps needed...\n" , reportheader ) ; fflush ( stderr ) ; 
            toalloc = nc4 ;
            together = 1 ; }}
    if ( ratcycles || ( ( dojak || quickbremer ) && together ) )
         wtchg = ( unsigned long int * ) myalloc ( nc4 * sizeof ( long int ) ) ;
    else if ( ( dojak || quickbremer ) && !together )
           wtchg = ( unsigned long int * ) myalloc ( piecesz * sizeof ( long int ) ) ;
    tmpp = fopen ( "oblong.tmp" , "rb" ) ;
    if ( !together ) {
        /** Alloc matrix for terms in a contiguous memory segment, for faster loading *******/
        segposit = ( unsigned long int * ) myalloc ( nt * toalloc * sizeof ( long int ) ) ;
        for ( a = 0 ; a < nt ; ++ a ) {
             mat [ a ] = segposit ;
             segposit += toalloc ; }
        for ( a = nt ; a < nodz ; ++ a )
            mat [ a ] = ( unsigned long int * ) myalloc ( toalloc * sizeof ( long int ) ) ; 
        create_parts_file( ) ; 
        fclose ( tmpp ) ; 
        remove ( "oblong.tmp" ) ;
        return ; }
    /**  Just alloc matrix as normally done ****/ 
    for ( a = 0 ; a < nodz ; ++ a ) {
        mat [ a ] = ( unsigned long int * ) myalloc ( nc4 * sizeof ( long int ) ) ;
        if ( a < nt ) 
           for ( b = 0 ; b < nc4 ; ++ b ) mat [ a ] [ b ] = -1 ; }
    prevprop = 0 ; 
    for ( a = 0 ; a < totnc ; ++ a ) {
       if ( reporttoscreen && ( a % 327 ) == 0 ) {
           prop =  ( double ) a / totnc ;
           if ( prop - prevprop > 0.015 ) {
               fprintf ( stderr , "\rReading matrix... %.1f%%" , 50 + ( prop * 100 / 2 ) ) ;
               prevprop = prop ; }}
       mybits = getc ( tmpp ) ;
       if ( mybits == 2 ) { 
          msk = 3 ; ent = atnc2 / FLDS2 ; inshif = 2 * ( atnc2 % FLDS2 ) ; 
          inmsk = ( unsigned long int ) 3 << inshif ; ++ atnc2 ; }
       if ( mybits == 3 ) { 
          msk = 7 ; ent = nc2 + ( atnc3 / FLDS3 ) ; inshif = 3 * ( atnc3 % FLDS3 ) ; 
          inmsk = ( unsigned long int ) 7 << inshif ; ++ atnc3 ; }
       if ( mybits == 4 ) { 
          msk = 15 ; ent = nc3 + ( atnc4 / FLDS4 ) ; inshif = 4 * ( atnc4 % FLDS4 ) ; 
          inmsk = ( unsigned long int ) 15 << inshif ; ++ atnc4 ; }
       totbytes = ( ( mybits * ( nt + 1 ) ) - 1 ) / 8 ; 
       for ( b = 0 ; b <= totbytes ; ++ b )
          tmpstates [ b ] = getc ( tmpp ) ;
       pt = ( unsigned short int * ) tmpstates ;
       for ( b = shif = 0 ; b < nt ; ++ b ) {
           lp = ( unsigned long int * ) pt ; 
           val = ( * lp >> shif ) & msk ;
           mat [ b ] [ ent ] &= ~inmsk ;
           mat [ b ] [ ent ] |= val << inshif ; 
           shif += mybits ;
           if ( shif >= 16 ) {
               ++ pt ;
               shif -= 16 ; }}}
    fclose ( tmpp ) ; 
    remove ( "oblong.tmp" ) ;
}

void get_fasta_dims ( int * mync , int * mynt )
{
    int i ;
    char c ; 
    int curnt = 0 , curnc = 0 , atnc ; 
    char * nampt = strsp ;
    while ( !feof ( inp ) && i != '>' ) i = getc ( inp ) ; //没读到 EOF 并且不是 > 继续，直到读到 >
    if ( i != '>' ) { //读到最后没有 > 报错
      fprintf ( stderr , "\nError reading fasta file --expected taxon name\n" ) ;
      errout ( 1 , NULL ) ; }
    while ( !feof ( inp ) ) { //没读到 EOF 不停
        fscanf ( inp , " %s" , nampt ) ; //从 inp 读到 nampt
        i = 0 ;
        ++ curnt ; 
        if ( curnt == 1) //第一次
            while ( 1 ) {
                fscanf ( inp , " %c" , &c ) ;
                if ( c == '>' || feof ( inp ) ) break ; //读到 > 或者 EOF
                ++ curnc ; } // 读到 > 或者 EOF curnc ++
        else { //第二次及之后
            atnc = 0 ; 
            while ( 1 ) {
                fscanf ( inp , " %c" , &c ) ;
                if ( c == '>' || feof ( inp ) ) break ;
                ++ atnc ; }
            if ( atnc != curnc ) { //到最后比较
                fprintf ( stderr , "\nError reading fasta file --it seems not to be aligned\n" ) ;
                errout ( 1 , NULL ) ; }}}
    * mync = curnc ;
    * mynt = curnt ; 
    rewind ( inp ) ; 
}

void read_multifiles ( void )
{
    char * nampt = strsp ;
    unsigned char state ;
    long int a , b , i , curnt = 0 , curnc = 0 , matched , prvnt ;
    int atfile = 0 , numfiles = 0 ; 
    char ** taxnam ; 
    double prop , secs , prevprop ;
    clock_t beg , end ; // 定义起始时间
    beg = clock () ; //计算时间
    errout ( ( tmpp = fopen ( "oblong.tmp" , "wb" ) ) == NULL , "Cannot open temporary file, \"oblong.tmp\"" ) ; 
    // Find out total NT and total NC
    multinp = inp ; 
    maskinit ( ) ; //maskinit 函数
    taxnam = ( char ** ) myalloc ( MAXMULTITAXA * sizeof ( char * ) ) ; //default MAXMULTITAXA 1000 L32
    fscanf ( inp , " %s" , nampt ) ; //FILE * inp = NULL 读入文件到 nampt
    while ( !feof ( multinp ) ) { //feof 如果 EOF 返回非 0，没有EOF 返回 0，
        ++ numfiles ; //numfiles 初始为 0
        if ( ( inp = fopen ( nampt , "rb" ) ) == NULL ) { // r readbale b binary 二进制可以防止 Windows下文件问题
            fprintf ( stderr , "\nCannot open file %s (from list of multiple files)\n" , nampt ) ;
            exit ( 0 ) ; }
        if ( reporttoscreen ) // L60: default 1 L309:0 
           fprintf ( stderr , "%cChecking input file %3i %-40s" , reportheader , numfiles , nampt ) ; 
        if ( fasta ) get_fasta_dims ( &nc , &nt ) ; //get nc and nt
        else {
           i = fscanf ( inp , " %i %i" , &nt , &nc ) ;
           if ( i < 2 ) {
               fprintf ( stderr , "\nError reading number of taxa/characters in data set %s\nMake sure it is in proper Phylip format\n" , bitct ) ;
               exit ( 0 ) ; }}
        if ( numfiles == 1 ) prvnt = nt ; 
        errout ( ( equal_taxon_sequence && nt != prvnt ) , "Under \"-e\" all data files must have the same number of taxa" ) ;
        curnc += nc ; 
        for ( a = 0 ; a < nt ; ++ a ) {
            fscanf ( inp , " %s" , nampt ) ;
            if ( fasta ) ++ nampt ; 
            if ( equal_taxon_sequence ) 
                if ( numfiles == 1 ) matched = -1 ;
                else matched = a ; 
            else 
               for ( i = 0 , matched = -1 ; i < curnt && matched < 0 ; ++ i )
                  if ( !strcmp ( nampt , taxnam[ i ]  ) ) matched = i ; 
            if ( matched < 0 ) {
                errout( curnt >= MAXMULTITAXA , "Sorry, cannot deal with so many distinct taxon names (change MAXMULTITAXA and recompile)" ) ;  
                taxnam [ curnt ] = ( char * ) myalloc ( MAXNAMELEN * sizeof ( char ) ) ; //MAXNAMELEN = 1000 default
                strcpy ( taxnam [ curnt ++ ] , nampt ) ; }
            if ( fasta ) -- nampt ; 
            for ( b = 0 ; b < nc ; ++ b ) {
               fscanf ( inp , " %c" , &state ) ;
               if ( !maskin [ state ] ) {
                   fprintf ( stderr , "\nUnrecognized symbol: %c\n" , state ) ;
                   exit ( 0 ) ; }}}
        fclose ( inp ) ;
        fscanf ( multinp , " %s" , nampt ) ; }
    // Now do some initializations and process
    rewind ( multinp ) ;
#ifdef WATCOM
    fpos = ( long long * ) myalloc ( curnt * sizeof ( long long ) ) ;
#else 
    fpos = ( fpos_t * ) myalloc ( curnt * sizeof ( fpos_t ) ) ;
#endif
    matbuf = ( unsigned char ** ) myalloc ( curnt * sizeof ( char * ) ) ;
    isabs = ( unsigned char * ) myalloc ( curnt * sizeof ( char ) ) ; 
    for ( a = 0 ; a < curnt ; ++ a ) 
        if ( curnc < CHUNKSZ )
             matbuf [ a ] = ( char * ) myalloc ( curnc * sizeof ( char ) ) ; 
        else matbuf [ a ] = ( char * ) myalloc ( CHUNKSZ * sizeof ( char ) ) ; 
    a = curnt ;
    if ( ( a & 1 ) ) ++ a ; 
    tmpstates = ( unsigned char * ) myalloc ( a * sizeof ( char ) ) ;
    fscanf ( multinp , " %s" , nampt ) ;
    while ( !feof ( multinp ) ) {
        if ( ( inp = fopen ( nampt , "rb" ) ) == NULL ) {  // it was just opened ok, but re-check just in case... 
            fprintf ( stderr , "\nCannot open file %s (from list of multiple files)\n" , nampt ) ;
            exit ( 0 ) ; }
        if ( !fasta ) fscanf ( inp , " %i %i" , &nt , &nc ) ;
        for ( a = 0 ; a < curnt ; ++ a ) isabs [ a ] = 1 ;
        for ( a = 0 ; a < nt ; ++ a ) {
           fscanf ( inp , " %s" , nampt ) ;
           if ( fasta ) ++ nampt ; 
           if ( equal_taxon_sequence )
                  i = a ;
           else 
             for ( i = 0 ; i < curnt ; ++ i )
                if ( !strcmp ( nampt , taxnam[ i ] ) ) break ;  // there has to be a match...
           if ( fasta ) -- nampt ;  
#ifdef WATCOM
           fpos [ i ] = _ftelli64 ( inp ) ; 
#else
           fgetpos ( inp , &fpos [ i ] ) ;
#endif  
           isabs [ i ] = 0 ;
           for ( b = 0 ; b < nc ; ++ b ) fscanf ( inp , " %c" , &state ) ; }
        nt = curnt ;
        atfile ++ ;
        prevprop = 0 ; 
        for ( a = 0 ; a < nc ; a += CHUNKSZ ) {
            if ( reporttoscreen ) {
                prop = ( double ) a / nc ;
                if ( ( prop - prevprop ) > 0.009 ) {
                    fprintf ( stderr , "%cReading matrix %i of %i... %.1f%%               " , reportheader , atfile , numfiles , prop * 100 ) ;
                    prevprop = prop ; }}
            process_chunk ( a ) ; }
        fclose ( inp ) ;
        fscanf ( multinp , " %s" , nampt ) ; }
    free ( matbuf ) ; 
    fclose ( tmpp ) ;
    // Assign bufers and (if appropriate), read data in, from tmp file...
    nt = curnt ; 
    tmpparse ( ) ; 
    fprintf ( stderr , "%cRead all data (%i files), %i taxa, %li chars. (%li informative)  \n" , reportheader , numfiles , nt , curnc , curnc - trashed ) ; 
    end = clock () ;
    secs = ( double ) ( end - beg ) / CLOCKS_PER_SEC ; 
    fprintf ( stderr , "Time to read data: %.2f secs.\n\n" , secs ) ; 
    // Begin writing to output file... 
    if ( !usenexus ) {
      if ( !savejusttrees ) 
          fprintf ( outp , "xread\n'Results from oblong'\n0 %i\n" , nt ) ; }
    else fprintf ( outp , "#NEXUS\nbegin trees;\nTranslate\n" ) ;
    for ( a = 0 ; a < nt ; ++ a ) 
        if ( !usenexus ) {
          if ( !savejusttrees ) 
            fprintf ( outp , "%s " , taxnam [ a ] ) ; }
        else {
            fprintf ( outp , "    %li %s" , a + 1 , taxnam [ a ] ) ;
            if ( a < nt - 1 ) fprintf ( outp , "," ) ;
            fprintf ( outp , "\n" ) ; }
    if ( !usenexus ) {
      if ( !savejusttrees ) 
        fprintf ( outp , ";\ncollapse none;\n" ) ; }
    else fprintf ( outp , "    ;\n" ) ; 
    return ;
}

void read_input ( void )
{
    char * nampt = strsp ;
    unsigned char state ;
    int a , b ;
    char c ; 
    double prop , secs , prevprop ;
    clock_t beg , end ;
    strcpy ( nampt , " " ) ; // nampt = " "
    if ( treeinp != NULL ) { // int is NULL (L34), but can specific the tree file (L273)
        while ( strcmp ( nampt , "tread" ) && !feof ( treeinp ) ) //feof 如果 EOF 返回非 0，没有EOF 返回 0， strcmp 0 相同，非0 则执行，即 非 tread 并且 没有 EOF 且执行
           fscanf ( treeinp , " %7s" , nampt ) ; // 输入长度为 7 的 string
        errout ( feof ( treeinp ) , "String \"tread\" not found in tree file" ) ;
        fscanf ( treeinp , " %c" , &c ) ;
        errout ( c != 39 , "Number of trees should be indicated in the title of tree-file" ) ;
        errout ( !fscanf ( treeinp , " %i" , &numreplications ) , "Number of trees should be indicated in the title of tree-file" ) ;
        while ( getc ( treeinp ) != 39 ) ; }
    if ( multifiles ) {
        read_multifiles ( ) ;
        return ; }
    beg = clock () ;
    if ( !felipe && !fasta ) {
      while ( strcmp ( nampt , "xread" ) ) {
          fscanf ( inp , " %7s" , nampt ) ;
          if ( feof ( inp ) ) break ; }
      errout ( feof ( inp ) , "String \"xread\" not found in data set\n\nMaybe it's in Phylip (\"-p\") or Fasta (\"-f\") format?" ) ;
      fscanf ( inp , " %c" , &c ) ;
      if ( c == 39 )
         while ( getc ( inp ) != 39 && !feof ( inp ) ) ;
      else ungetc ( c , inp ) ; }
    if ( fasta ) get_fasta_dims ( &nc , &nt ) ; 
    else fscanf ( inp , " %i %i" , &nc , &nt ) ;
    if ( felipe ) {
        a = nc ;
        nc = nt ;
        nt = a ; }
    errout ( !nc , "Illegal number of characters" ) ; 
    errout ( nt < 4 , "Illegal number of taxa" ) ; 
#ifdef WATCOM
    fpos = ( long long * ) myalloc ( nt * sizeof ( long long ) ) ;
#else 
    fpos = ( fpos_t * ) myalloc ( nt * sizeof ( fpos_t ) ) ;
#endif
    matbuf = ( unsigned char ** ) myalloc ( nt * sizeof ( char * ) ) ;
    a = nt ;
    if ( ( a & 1 ) ) ++ a ; 
    tmpstates = ( unsigned char * ) myalloc ( a * sizeof ( char ) ) ;
    maskinit ( ) ; 
    // Get file positions for first character for each taxon...