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...
    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 ) {
        fscanf ( inp , " %s" , nampt ) ;
        if ( fasta ) ++ nampt ; 
        if ( !usenexus ) {
          if ( !savejusttrees ) 
            fprintf ( outp , "%s " , nampt ) ; }
        else {
            fprintf ( outp , "    %i %s" , a + 1 , nampt ) ;
            if ( a < nt - 1 ) fprintf ( outp , "," ) ;
            fprintf ( outp , "\n" ) ; }
        if ( fasta ) -- nampt ; 
#ifdef WATCOM
        fpos [ a ] = _ftelli64 ( inp ) ; 
#else
        fgetpos ( inp , &fpos [ a ] ) ;
#endif        
        if ( nc < CHUNKSZ )
             matbuf [ a ] = ( char * ) myalloc ( nc * sizeof ( char ) ) ; 
        else matbuf [ a ] = ( char * ) myalloc ( CHUNKSZ * sizeof ( char ) ) ; 
        for ( b = 0 ; b < nc ; ++ b ) {
            fscanf ( inp , " %c" , &state ) ;
            if ( state == '[' ) {
              while ( 1 ) {
               fscanf ( inp , " %c" , &state ) ;
               if ( state == ']' ) break ; 
               if ( !maskin [ state ] ) {
                  fprintf ( stderr , "\nUnrecognized symbol in matrix: %c\n" , state ) ;
                  exit ( 0 ) ; }}}
            else 
            if ( !maskin [ state ] ) {
                fprintf ( stderr , "\nUnrecognized symbol in matrix: %c\n" , state ) ;
                exit ( 0 ) ; }}}
    if ( !usenexus ) {
      if ( !savejusttrees ) 
        fprintf ( outp , ";\ncollapse none;\n" ) ; }
    else fprintf ( outp , "    ;\n" ) ; 
    // Read pieces, and output temporary file as you go...
    errout ( ( tmpp = fopen ( "oblong.tmp" , "wb" ) ) == NULL , "Cannot open temporary file, \"oblong.tmp\"" ) ; 
    prevprop = 0 ; 
    for ( a = 0 ; a < nc ; a += CHUNKSZ ) {
        if ( reporttoscreen ) {
            prop = ( double ) a / nc ;
            if ( together ) prop /= 2 ;
            if ( ( prop - prevprop ) > 0.009 ) {
                fprintf ( stderr , "\rReading matrix... %.1f%%" , prop * 100 ) ;
                prevprop = prop ; }}
        process_chunk ( a ) ; }
    free ( matbuf ) ; 
    fclose ( tmpp ) ;
    // Assign bufers and (if appropriate), read data in, from tmp file...
    tmpparse ( ) ; 
    fprintf ( stderr , "%cReading matrix, %i taxa, %i chars. (%li informative)  \n" , reportheader , nt , nc , nc - trashed ) ; 
    fclose ( inp ) ; 
    end = clock () ;
    secs = ( double ) ( end - beg ) / CLOCKS_PER_SEC ; 
    fprintf ( stderr , "Time to read data: %.2f secs.\n\n" , secs ) ; 
    return ;
}        

/************  GENERAL TREE-HANDLING STUFF ************/

void reset_nodz ( short int * anpt )
{
  short int * bu = rlist ; 
  short int * tbu = tgtlist ; 
  int a , b , atno , lg ; 
  for ( a = nodz ; a -- > nt ; ) tbu [ a ] = -1;
  for ( a = nt ; a -- ; ) tbu [ a ] = a ; 
  atno = nt ;
  lg = anpt [ nt ] ; 
  tbu [ nt ] = nt ;
  tbu [ lg ] = lg ; 
  for ( a = 0 ; a < nt ; ++a) {
      b = anpt [ a ] ;
      while ( tbu [ b ] < 0 ) {
          tbu [ b ] = ++ atno ; 
          b = anpt [ b ] ; }}
  for ( a = lg ; a-- ; ) 
      bu [ tbu [ a ] ] = tbu [ anpt [ a ] ] ;
  for ( a = lg ; a-- ; ) anpt [ a ] = bu [ a ] ;
  anpt [ nt ] = ++ atno ; 
  return ; 
}

int equal_trees ( short int * a , short int * b )
{
    int c ; 
    for ( c = 0 ; c < nodz ; ++ c )
       if ( * a ++ != * b ++ ) return 0 ;
    return 1 ;
}    
            
void storetree ( void )
{
    int b , isdup = 0 ;
    for ( b = 0 ; b < nodz ; ++ b ) memtrees [ numtrees ] [ b ] = anx [ b ] ;
    reset_nodz ( memtrees [ numtrees ] ) ; 
    if ( !dojak && !saveall ) 
       for ( b = 0 ; b < numtrees && !isdup ; ++ b )
          if ( equal_trees ( memtrees [ b ] , memtrees [ numtrees ] ) ) isdup = 1 ; 
    if ( !isdup ) 
        ++ numtrees ;
}

void randomlist ( int num , short int * dali , int firstis )
{ short int * bakrand ;
  int a , todo , i , n ;  
  if ( ranseed ) srand ( ranseed ) ; 
  for ( a = 0 ; a < num ; ++ a ) dali [ a ] = a ; 
  if ( firstis == 0 ) {
     todo = num - 1 ; 
     ++ dali ;  // skip 0
     bakrand = dali + todo ; 
     for ( a = 1 ; a < num ; ++ a ) {
        i = rand () % todo -- ; 
        n = * -- bakrand ; 
        * bakrand = dali [ i ] ; 
        dali [ i ] = n ; }}
  else {
     todo = num ; 
     bakrand = dali + todo ; 
     for ( a = 0 ; a < num ; ++ a ) {
        i = rand () % todo -- ; 
        n = * -- bakrand ; 
        * bakrand = dali [ i ] ; 
        dali [ i ] = n ; }}
}

void nodzert ( int what, int where )
{ 
    int a, b, c ;
    a = anx [ where ] ;
    c = sis [ where ] ; 
    b = anx [ what ] = anx[ where ] = anx[nt]++;
    anx [ b ] = a;
    sis [ what ] = where;
    sis [ where ] = what;
    sis [ c ] = b;
    sis [ b ] = c;
    lef [ a ] = c;
    rig [ a ] = b;
    lef [ b ] = where ;
    rig [ b ] = what ; 
}     


int savebremervals = 0 ;

void savant ( int which )
{
    int a , didsome = 0 ;
    if ( savebremervals && which != nt ) {
       if ( gbrem [ which ] > 0 || shownegbremers ) fprintf ( outp , "(" ) ; }
    else 
      fprintf ( outp , "(" ) ;
    for ( a = 0 ; a < anx [ nt ] ; ++ a )
      if ( anx [ a ] == which ) {
         if ( usenexus && didsome ) 
             fprintf ( outp , "," ) ;
         didsome = 1 ; 
         if ( a < nt ) 
             if ( usenexus ) fprintf ( outp , "%i " , a + 1 ) ;
             else fprintf ( outp , "%i " , a ) ;
         else savant ( a ) ; }
    if ( savebremervals && which != nt && anx [ which ] != nt ) {
       if ( gbrem [ which ] > 0 )
          if ( bremertype == 0 ) fprintf ( outp , ")%s%.0f " , legindic , gbrem [ which ] ) ; 
          else fprintf ( outp , ")%s%.0f " , legindic , gbrem [ which ] * 100 ) ;
       else if ( shownegbremers ) 
               if ( bremertype == 0 ) fprintf ( outp , ")%s%.0f " , legindic , gbrem [ which ] ) ; 
               else fprintf ( outp , ")%s(%.0f) " , legindic , gbrem [ which ] * 100 ) ; }
    else 
      fprintf ( outp , ")" ) ;
}

void savetrees ( void )
{
   short int * anxwas = anx ;
   int i ; 
   if ( quickbremer )
     fprintf ( stderr , "\nSaving supports and %i trees to output file...\n\n" , numtrees ) ; 
   else 
     fprintf ( stderr , "\nSaving %i trees to output file...\n\n" , numtrees ) ; 
   if ( !usenexus ) 
       fprintf ( outp , "tread '%i Oblong trees'\n" , numtrees ) ;
   for ( i = 0 ; i < numtrees ; ++ i ) {
      if ( quickbremer && !i ) savebremervals = 1 ;
      if ( !usenexus ) {
           if ( i ) fprintf ( outp , "*\n" ) ; }
      else fprintf ( outp , "tree OBLONG_%i = [&U] " , i + 1 ) ; 
      anx = memtrees [ i ] ;
      savant ( nt ) ;
      if ( usenexus ) 
         fprintf ( outp , ";\n" ) ;
      if ( quickbremer && !i ) {
          savebremervals = quickbremer = 0 ;
          if ( !usenexus ) fprintf ( outp , "*\n" ) ; 
          i = -1 ; }}
   anx = anxwas ;
   if ( !usenexus ) 
       fprintf ( outp , ";\nproc/;\n" ) ;
   else fprintf ( outp , "End ; \n" ) ; 
}       

int binoptlist ( short int * list , short int * tt )
{
    short int * x , * y ;
    x = y = list ; 
    * y ++ = nt ;
    while ( x < y ) {
        if ( lef [ * x ] >= nt ) * y ++ = lef [ * x ] ;
        if ( rig [ * x ] >= nt ) * y ++ = rig [ * x ] ;
        ++ x ; }
    return y - list ;
}

void initnet ( int un, int dos, int tres )
{
    anx [ un ] = anx [ nt + 1] = nt;
    anx [ dos ] = anx [ tres ] = nt + 1;
    anx [ nt ] = nt + 2 ; 
    lef [ nt ] = un;
    rig [ nt ] = sis [ un ] = nt + 1;
    lef [ nt + 1 ] = dos;
    rig [ nt + 1 ] = sis [ dos ] = tres;
    sis [ nt + 1 ] = un; 
    sis [ tres ] = dos;
} 

/************ OPTIMIZATION STUFF *********************/

void dowtfilter ( int prob , int keepprevious , int dalim )
{
    int i , j ;
    unsigned long int msk ; 
    for ( i = 0 ; i < nc2 ; ++ i ) {
        if ( keepprevious )
             wtchg [ i ] |= wtchg [ i ] << 1 ; 
        else wtchg [ i ] = LO2 ; 
        for ( j = 0 , msk = 1 ; j < FLDS2 ; ++ j , msk <<= 2 ) 
            if ( ( rand () % 100 ) < prob ) 
                wtchg [ i ] &= ~msk ; }
    for (       ; i < nc3 ; ++ i ) {
        if ( keepprevious )
             wtchg [ i ] |= wtchg [ i ] >> 1 ; 
        else wtchg [ i ] = HI3 ; 
        for ( j = 0 , msk = 4 ; j < FLDS3 ; ++ j , msk <<= 3 ) 
            if ( ( rand () % 100 ) < prob ) 
                wtchg [ i ] &= ~msk ; }
    for (       ; i < dalim ; ++ i ) {
        if ( keepprevious )
             wtchg [ i ] |= wtchg [ i ] >> 1 ; 
        else wtchg [ i ] = HI4 ; 
        for ( j = 0 , msk = 8 ; j < FLDS4 ; ++ j , msk <<= 4 ) 
            if ( ( rand () % 100 ) < prob ) 
                wtchg [ i ] &= ~msk ; }
}       

void resetwtfilter ( void )
{
    int i ;
    for ( i = 0 ; i < nc2 ; ++ i ) 
        wtchg [ i ] = ( wtchg [ i ] >> 1 ) & LO2 ; 
    for (       ; i < nc3 ; ++ i ) 
        wtchg [ i ] = ( wtchg [ i ] << 1 ) & HI3 ; 
    for (       ; i < nc4 ; ++ i ) 
        wtchg [ i ] = ( wtchg [ i ] << 1 ) & HI4 ; 
}       

void storewtfilter ( void )
{
    int i ;
    for ( i = 0 ; i < nc2 ; ++ i ) 
        wtchg [ i ] |= ( wtchg [ i ] << 1 ) ; 
    for (       ; i < nc3 ; ++ i ) 
        wtchg [ i ] |= ( wtchg [ i ] >> 1 ) ; 
    for (       ; i < nc4 ; ++ i ) 
        wtchg [ i ] |= ( wtchg [ i ] >> 1 ) ; 
}       


unsigned long long downpass_all ( void )
{
    int nunos , len , a , nod , c ;
    unsigned long long totlen ; 
    unsigned long int * lp , * rp , * np , x , y , z ;
    unsigned long int * wts ; 
    nunos = binoptlist ( nodlist , anx ) ;
    totlen = lenlag ;
    for ( a = nunos ; a -- ; ) {
        nod = nodlist [ a ] ;
        lp = mat [ lef [ nod ] ] ; 
        rp = mat [ rig [ nod ] ] ; 
        np = mat [ nod ] ; 
        wts = wtchg - 1 ;
        len = 0 ; 
        for ( c = 0 ; c < nc2 ; ++ c ) {
            x = * lp ++ & * rp ++ ;
            y = LO2 ^ ( ( x & LO2 ) | ( ( x & HI2 ) >> 1 ) ) ;
            if ( changewts ) y &= * ++ wts ; 
#ifdef BITS64
            len += bitct [ 65535 & ( y | ( y >> 15 ) ) ] + bitct [ 65535 & ( ( y >> 32 ) | ( y >> 47 ) ) ] ; 
#else 
            len += bitct [ 65535 & ( y | ( y >> 15 ) ) ] ; 
#endif
            y |= y << 1 ;
            * np ++ = x | y ; }
        for (   ; c < nc3 ; ++ c ) {
            x = * lp & * rp ;
            y = * lp ++ | * rp ++ ;
            z = HI3 & ~( x | ( ( x & LO3 ) + LO3 ) ) ;
            if ( changewts ) z &= * ++ wts ; 
#ifdef BITS64
            len += bitct [ 65535 & ( z | ( z >> 16 ) ) ] + bitct [ 65535 & ( z >> 32 | ( z >> 49 ) ) ] ;  
#else 
            len += bitct [ 65535 & ( z | ( z >> 17 ) ) ] ;
#endif
            z |= z - ( z >> 2 ) ;
            * np ++ = x | ( z & y ) ; }
        for (   ; c < nc4 ; ++ c ) {
            x = * lp & * rp ;
            y = * lp ++ | * rp ++ ;
            z = HI4 & ~( x | ( ( x & LO4 ) + LO4 ) ) ;
            if ( changewts ) z &= * ++ wts ; 
#ifdef BITS64
            len += bitct [ 65535 & ( z | ( z >> 17 ) | ( z >> 34 ) | ( z >> 51 ) ) ] ;
#else 
            len += bitct [ 65535 & ( z | ( z >> 17 ) ) ] ;
#endif
            z |= z - ( z >> 3 ) ;
            * np ++ = x | ( z & y ) ; }
        totlen += len ; }
    return totlen ; 
}    

void loadmat ( int part ) 
{
    void * vd ;
    unsigned long int * seg ;
    unsigned long int numelems ; 
    if ( !part ) rewind ( partsfile ) ;
    seg = mat [ 0 ] ;
    vd = ( void * ) seg ; 
    numelems = fread ( vd , 1 , sizeof ( unsigned long int ) * nt * piecesz , partsfile ) ;
    errout ( numelems != sizeof ( unsigned long int ) * nt * piecesz , "OOOPS, unexpected error loading matrix!.\nSorry, abandoning run..." ) ;
    if ( dojak && !ratcycles && !together ) {
        if ( !part ) srand ( ranseed ) ;
        dowtfilter ( jakprob , 0 , piecesz ) ; }
}    

unsigned long long downpass_split ( void )
{
    int nunos , len , a , nod , c ;
    unsigned long long totlen ; 
    unsigned long int * lp , * rp , * np , x , y , z ;
    unsigned long int atnc ;
    unsigned long int * wts ; 
    nunos = binoptlist ( nodlist , anx ) ;
    totlen = lenlag ;
    for ( atnc = 0 ; atnc < nc4 ; atnc += piecesz ) {
        loadmat ( atnc ) ; 
        for ( a = nunos ; a -- ; ) {
            nod = nodlist [ a ] ;
            lp = mat [ lef [ nod ] ] ; 
            rp = mat [ rig [ nod ] ] ; 
            np = mat [ nod ] ;
            if ( ratcycles ) wts = wtchg + atnc - 1 ;
            else wts = wtchg - 1 ;
            len = 0 ; 
            for ( c = 0 ; c < piecesz ; ++ c ) {
                x = * lp & * rp ;
                y = * lp ++ | * rp ++ ;
                z = HI4 & ~( x | ( ( x & LO4 ) + LO4 ) ) ;
                if ( changewts ) z &= * ++ wts ; 
#ifdef BITS64
                len += bitct [ 65535 & ( z | ( z >> 17 ) | ( z >> 34 ) | ( z >> 51 ) ) ] ;
#else 
                len += bitct [ 65535 & ( z | ( z >> 17 ) ) ] ;
#endif
                z |= z - ( z >> 3 ) ;
                * np ++ = x | ( z & y ) ; }
            totlen += len ; 
                }}
    return totlen ; 
}    

/*********  SUPPORT FUNCTIONS FOR WAGNER AND BRANCH-SWAPPING ************************/

void zert ( int whatp, int wherep )
{
    int a, b, c ;
    if ( wherep == 0 ) {
        anx [ lef [ nt ] = sis [ whatp ] = 0 ] = nt ; 
        anx [ rig [ nt ] = sis [ 0 ] = whatp ] = nt ;
        return ; } 
    a = anx [ wherep ] ;
    c = sis [ wherep ] ; 
    b = anx [ wherep ] = anx [ whatp ] ; 
    anx [ b ] = a;
    sis [ sis [ wherep ] = rig [ b ] = whatp ] = lef [ b ] = wherep;
    sis [ lef [ a ] = c ] = b ;
    sis [ rig [ a ] = b ] = c ;
}     

void unzert ( int whatp )
{  
    int a, b, c;
    a = anx [ whatp ] ;
    if ( a == nt ) {
        lef [ anx [ 0 ] = nt ] = rig [ nt ] = sis [ 0 ] = 0 ;
        anx [ whatp ] = whatp ;
        return ; } 
    b = sis [ whatp ] ;
    c = sis [ a ] ; 
    a = anx [ a ] ; 
    lef [ a ] = b ;
    rig [ anx [ b ] = anx [ c ] = a ] = sis [ sis [ c ] = b ] = c ;
}

void tbreroot ( int new, int old )
{
   short int * an = anx ; 
   int i , j , k , m , n , nx , first ;
   if ( an [ new ] == old || old < nt ) return ; 
   j = nx = first = an [ i = new ] ;
   n = sis [ new ] ;
   while ( nx != old ) {
       k = i ; 
       i = j ;
       j = nx ;
       m = n ;
       n = sis [ j ] ; 
       if ( k == lef [ i ] ) { lef [ i ] = j ; sis [ sis [ rig [ i ] ] = j ] = rig [ i ] ; } 
       else { rig [ i ] = j ; sis [ sis [ lef [ i ] ] = j ] = lef [ i ] ; }
       nx = an [ j ] ; 
       an [ j ] = i ; }
   an [ lef [ j ] = m ] = an [ rig [ j ] = n ] = j ;
   sis [ sis [ n ] = m ] = n ;
   an [ lef [ old ] = new ] = an [ rig [ old ] = first ] = old ;
   sis [ sis [ new ] = first ] = new ; 
   return ;
}

void binmksis ( short int * tt ) 
{
    int a , b , c , lownod = tt [ nt ] ;
    for ( a = nodz ; a -- ; ) nodfork [ a ] = 0 ;
    for ( a = nt ; a -- ; ) {
        b = tt [ c = a ] ;
        while ( b != lownod ) {
                 if ( !nodfork [ b ] ) 
                     lef [ b ] = lassis [ b ] = c ;
                 else 
                     sis [ sis [ c ] = lassis [ b ] ] = rig [ b ] = c ;
                 if ( ++ nodfork [ b ] > 1 ) break ;
                 c = b ;
                 b = tt [ b ] ; }}
    return ; 
}         

int divoptlist ( short int * list , int clip , int sisis )
{
    short int * x , * y ;
    x = y = list ; 
    if ( sisis ) * y ++ = nt ;
    if ( clip > nt ) * y ++ = clip ; 
    while ( x < y ) {
        if ( lef [ * x ] >= nt ) * y ++ = lef [ * x ] ;
        if ( rig [ * x ] >= nt ) * y ++ = rig [ * x ] ;
        ++ x ; }
    return y - list ;
}

/***********  HANDLE SWAP-COLLAPSING FOR JACKKNIFING *************/

void delnodz ( void )
{ short int  * tmp , *cs ;
  int a , b , at , last ;
  char stopped ;
  short int * an = anx ;
  short int * dellist = unsuprt ; 
  tmp = nodlist ;
  cs =  lef ; 
  dellist[ nt ] = 0 ;
  at = nt + 1 ;
  for ( a = nt ; a -- ; ) cs [ a ] = a ;
  for ( a = nt ; a < nodz ; ++ a ) cs [ a ] = 0 ; 
  cs [ nt ] = nt ; 
  for ( a = nodz ; a -- ; ) tmp[ a ] = 0 ; 
  for ( a = nt ; a -- ; ) {
      b = last = a ;
      stopped = 0 ; 
      while ( b != nt ) {
          b = an [ b ] ;
          if ( cs [ b ] ) { tmp[ cs [ last ] ] = cs [ b ] ; stopped = 1 ; break ; } 
          else
          if ( ! dellist [ b ] ) { 
             tmp [ cs [ last ] ] = at ;
             cs [ last = b ] = at ++ ; } 
         } 
      if ( !stopped ) tmp [ cs [ last ] ] = nt ; }
  for ( a = nodz ; a -- ; ) an [ a ] = tmp [ a ] ;
  an [ nt ] = at ;
} 

void filldistoroot ( void )
{
    short int * x = nodlist , * y ;
    binmksis ( anx ) ; 
    distoroot [ nt ] = 0 ;
    nodlist [ 0 ] = lef [ nt ] ;
    nodlist [ 1 ] = rig [ nt ] ;
    x = nodlist - 1 ; 
    y = nodlist + 2 ; 
    while ( ++ x < y ) {
        distoroot [ * x ] = distoroot [ anx [ * x ] ] + 1 ;
        if ( * x > nt ) {
            * y ++ = lef [ * x ] ;
            * y ++ = rig [ * x ] ; }}
    return ; 
}         

void collapse_with_swaps ( int clip , int sisis )
{
    int a , b , c , comnod ;
    int da , db , dest , rootis ;
    zert ( clip , sisis ) ;
    for ( c = 0 ; c < numcolmoves ; ++ c ) {
        dest = colmoves [ c ] ; 
        if ( dest == sisis ) continue ; 
        da = distoroot [ a = clip ] ; 
        db = distoroot [ b = dest ] ;
        while ( distoroot [ a ] > db ) a = anx [ a ] ;
        while ( distoroot [ b ] > da ) b = anx [ b ] ;
        while ( a != b ) { a = anx [ a ] ; b = anx [ b ] ; }
        comnod = a ;
        if ( comnod != dest ) 
            for ( a = anx [ dest ] ; a != comnod ; a = anx [ a ] ) unsuprt [ a ] = 1 ; 
        for ( b = anx [ clip ] ; b != comnod ; b = anx [ b ] ) unsuprt [ b ] = 1 ; }
    if ( !use_spr && clip > nt ) 
      for ( c = 0 ; c < numcolroots ; ++ c ) {
        rootis = colroots [ c ] ; 
        if ( anx [ rootis ] != clip ) {
            a = anx [ rootis ] ;
            while ( a != clip ) {
               unsuprt [ a ] = 1 ;
               a = anx [ a ] ; }}}
    unzert ( clip ) ; 
    return ; 
}

void collapse_the_tree ( void )
{
    int i ;
    if ( !swapcollapse ) return ;
    for ( i = 0 ; i < nodz ; ++ i ) unsuprt [ i ] = 0 ;
    filldistoroot ( ) ;
    recordgoodmoves = 1 ;
    bswap ( ) ; 
    recordgoodmoves = 0 ;
    delnodz ( ) ; 
}

/***************  WAGNER AND BRANCH-SWAPPING FOR SINGLE BUFFER ****************/

void fullpass_all ( int clip , int sisis )
{
    int nunos , a , nod , c ;
    unsigned long int * lp , * rp , * np , * ap , x , y , z , bak ; 
    nunos = divoptlist ( nodlist , clip , sisis ) ;
    for ( a = nunos ; a -- ; ) {
        nod = nodlist [ a ] ;
        if ( nodtoskip == nod ) continue ; 
        lp = mat [ lef [ nod ] ] ; 
        rp = mat [ rig [ nod ] ] ; 
        np = mat [ nod ] ; 
        for ( c = 0 ; c < nc2 ; ++ c ) {
            x = * lp ++ & * rp ++ ;
            y = LO2 ^ ( ( x & LO2 ) | ( ( x & HI2 ) >> 1 ) ) ;
            y |= y << 1 ;
            * np ++ = x | y ; }
        for (   ; c < nc3 ; ++ c ) {
            x = * lp & * rp ;
            y = * lp ++ | * rp ++ ;
            z = HI3 & ~( x | ( ( x & LO3 ) + LO3 ) ) ;
            z |= z - ( z >> 2 ) ;
            * np ++ = x | ( z & y ) ; }
        for (   ; c < nc4 ; ++ c ) {
            x = * lp & * rp ;
            y = * lp ++ | * rp ++ ;
            z = HI4 & ~( x | ( ( x & LO4 ) + LO4 ) ) ;
            z |= z - ( z >> 3 ) ;
            * np ++ = x | ( z & y ) ; }}
    for ( a = 0 ; ++ a < nunos ; ) {
        nod = nodlist [ a ] ;
        if ( nod == clip || anx [ nod ] == nodtoskip || nod == nodtoskip ) continue ; 
        lp = mat [ lef [ nod ] ] - 1 ; 
        rp = mat [ rig [ nod ] ] - 1 ; 
        np = mat [ nod ] - 1 ;
        ap = mat [ anx [ nod ] ] - 1 ;
        for ( c = 0 ; c < nc2 ; ++ c ) {
            bak = * ++ lp | * ++ rp ;
            x = * ++ ap & ~ * ++ np ;
            x = ( x & LO2 ) | ( ( x & HI2 ) >> 1 ) ;
            x |= ( x << 1 ) ;
            * np = ( * ap & ~ x ) | ( x & ( * np | * ap & bak ) ) ; } 
        for (   ; c < nc3 ; ++ c ) {
            x = * ++ lp & * ++ rp ;
            z = * lp | * rp ;
            y = HI3 & ~( x | ( ( x & LO3 ) + LO3 ) ) ;
            y |= y - ( y >> 2 ) ;
            bak = z | y ; 
            x = * ++ ap & ~ * ++ np ;
            x = ( x | ( ( x & LO3 ) + LO3 ) ) & HI3 ;
            x |= x - ( x >> 2 ) ;
            * np = ( * ap & ~ x ) | ( x & ( * np | * ap & bak ) ) ; } 
        for (   ; c < nc4 ; ++ c ) {
            x = * ++ lp & * ++ rp ;
            z = * lp | * rp ; 
            y =  HI4 & ~( x | ( ( x & LO4 ) + LO4 ) ) ; 
            y |= y - ( y >> 3 ) ; 
            bak = z | y  ;
            x = * ++ ap & ~ * ++ np ;
            x = ( x | ( ( x & LO4 ) + LO4 ) ) & HI4 ;
            x |= x - ( x >> 3 ) ;
            * np = ( * ap & ~ x ) | ( x & ( * np | * ap & bak ) ) ; }}
    nodtoskip = 0 ; 
}    

unsigned long int scorediff_all ( short int clipis , short int srcis , short int tgtis )
{
    unsigned long int len = 0 ;
    unsigned long int * tgt , * tgta , * src , * srca ;
    unsigned long int untgt , unsrc ;
    unsigned long int x , y , z ;
    unsigned long int * wts ; 
    long int a ;
    if ( clipis < nt ) srca = mat [ clipis ] - 1 ;
    else srca = mat [ anx [ srcis ] ] - 1 ;
    src = mat [ srcis ] - 1 ;
    if ( !tgtis ) tgta = mat [ 0 ] - 1 ;
    else tgta = mat [ anx [ tgtis ] ] - 1 ;
    tgt = mat [ tgtis ] - 1 ;
    wts = wtchg - 1 ; 
    for ( a = 0 ; a < nc2 ; ++ a ) {
        unsrc = * ++ src | * ++ srca ; 
        if ( srcis < nt && clipis > nt ) {
           x = * srca & ~ * src ;
           x = ( x & LO2 ) | ( ( x & HI2 ) >> 1 ) ;
           x |= ( x << 1 ) ;
           unsrc = * srca | ( * src & x ) ; }
        untgt = * ++ tgt | * ++ tgta ;
        if ( tgtis && tgtis < nt ) {
           x = * tgta & ~ * tgt ;
           x = ( x & LO2 ) | ( ( x & HI2 ) >> 1 ) ;
           x |= ( x << 1 ) ;
           untgt = * tgta | ( * tgt & x ) ; }
        x = unsrc & untgt ;
        y = LO2 ^ ( ( x & LO2 ) | ( ( x & HI2 ) >> 1 ) ) ;
        if ( changewts ) 
           if ( recorddn ) {
              ++ wts ; 
              if ( recorddn == 1 ) * wts = y ;
              x = * wts & y ;
              z = x ^ y ;
#ifdef BITS64
              gmix += bitct [ 65535 & ( x | ( x >> 15 ) ) ] + bitct [ 65535 & ( ( x >> 32 ) | ( x >> 47 ) ) ] ; 
              gup += bitct [ 65535 & ( z | ( z >> 15 ) ) ] + bitct [ 65535 & ( ( z >> 32 ) | ( z >> 47 ) ) ] ; 
#else 
              gmix += bitct [ 65535 & ( x | ( x >> 15 ) ) ] ;
              gup += bitct [ 65535 & ( z | ( z >> 15 ) ) ] ;
#endif
              continue ; }
          else 
            y &= * ++ wts ; 
#ifdef BITS64
        len += bitct [ 65535 & ( y | ( y >> 15 ) ) ] + bitct [ 65535 & ( ( y >> 32 ) | ( y >> 47 ) ) ] ; 
#else 
        len += bitct [ 65535 & ( y | ( y >> 15 ) ) ] ; 
#endif
        if ( len > lenlimit ) return len ; }
    for (       ; a < nc3 ; ++ a ) {
        unsrc = * ++ src | * ++ srca ; 
        if ( srcis < nt && clipis > nt ) {
           x = * srca & ~ * src ;
           x = ( x | ( ( x & LO3 ) + LO3 ) ) & HI3 ;
           x |= x - ( x >> 2 ) ;
           unsrc = * srca | ( * src & x ) ; }
        untgt = * ++ tgt | * ++ tgta ;
        if ( tgtis && tgtis < nt ) {
           x = * tgta & ~ * tgt ;
           x = ( x | ( ( x & LO3 ) + LO3 ) ) & HI3 ;
           x |= x - ( x >> 2 ) ;
           untgt = * tgta | ( * tgt & x ) ; } 
        x = unsrc & untgt ;
        z = HI3 & ~( x | ( ( x & LO3 ) + LO3 ) ) ;
        if ( changewts ) 
          if ( recorddn ) {
             ++ wts ; 
             if ( recorddn == 1 ) * wts = z ;
             y = * wts & z ; 
             x = y ^ z ; 
#ifdef BITS64
             gmix += bitct [ 65535 & ( y | ( y >> 16 ) ) ] + bitct [ 65535 & ( y >> 32 | ( y >> 49 ) ) ] ;  
             gup += bitct [ 65535 & ( x | ( x >> 16 ) ) ] + bitct [ 65535 & ( x >> 32 | ( x >> 49 ) ) ] ;  
#else 
             gmix += bitct [ 65535 & ( y | ( y >> 17 ) ) ] ;
             gup += bitct [ 65535 & ( x | ( x >> 17 ) ) ] ;
#endif
             continue ; }
          else 
             z &= * ++ wts ; 
#ifdef BITS64
        len += bitct [ 65535 & ( z | ( z >> 16 ) ) ] + bitct [ 65535 & ( z >> 32 | ( z >> 49 ) ) ] ;  
#else 
        len += bitct [ 65535 & ( z | ( z >> 17 ) ) ] ;
#endif
        if ( len > lenlimit ) return len ; }
    for (       ; a < nc4 ; ++ a ) {
        unsrc = * ++ src | * ++ srca ; 
        if ( srcis < nt && clipis > nt ) {
           x = * srca & ~ * src ;
           x = ( x | ( ( x & LO4 ) + LO4 ) ) & HI4 ;
           x |= x - ( x >> 3 ) ;
           unsrc = * srca | ( * src & x ) ; }
        untgt = * ++ tgt | * ++ tgta ;
        if ( tgtis && tgtis < nt ) {
           x = * tgta & ~ * tgt ;
           x = ( x | ( ( x & LO4 ) + LO4 ) ) & HI4 ;
           x |= x - ( x >> 3 ) ;
           untgt = * tgta | ( * tgt & x ) ; } 
        x = unsrc & untgt ;
        z = HI4 & ~( x | ( ( x & LO4 ) + LO4 ) ) ;
        if ( changewts )
          if ( recorddn ) {
             ++ wts ; 
             if ( recorddn == 1 ) * wts = z ;
             y = * wts & z ; 
             x = y ^ z ; 
#ifdef BITS64
             gmix += bitct [ 65535 & ( y | ( y >> 17 ) | ( y >> 34 ) | ( y >> 51 ) ) ] ;
             gup += bitct [ 65535 & ( x | ( x >> 17 ) | ( x >> 34 ) | ( x >> 51 ) ) ] ;
#else 
             gmix += bitct [ 65535 & ( y | ( y >> 17 ) ) ] ;
             gup += bitct [ 65535 & ( x | ( x >> 17 ) ) ] ;
#endif
             continue ; }
          else 
           z &= * ++ wts ; 
#ifdef BITS64
        len += bitct [ 65535 & ( z | ( z >> 17 ) | ( z >> 34 ) | ( z >> 51 ) ) ] ;
#else 
        len += bitct [ 65535 & ( z | ( z >> 17 ) ) ] ;
#endif
        if ( len > lenlimit ) return len ; }
    return len ; 
}

unsigned long long wagner_all ( void )
{
    int a , added, fac, basnod , bestpos , clip , tgt ;
    unsigned long int lim , diff ;
    unsigned long long len ; 
    short int * inlist = nodfork ; 
    unsigned long int seedwas = ranseed ; 
    randomlist( nt , ladd , 0 ) ;
    initnet ( ladd [ 0 ] , ladd [ 1 ] , ladd [ 2 ] ) ; 
    blist [ 0 ] = ladd [ 1 ] ; 
    blist [ 1 ] = ladd [ 2 ] ; 
    blist [ 2 ] = nt + 1;
    added = 2;
    fac = 3 ;
    basnod = nt + 2;
    ranseed = 0 ; 
    while ( ++ added < nt ) {
        if ( random_starts ) {
           nodzert ( ladd [ added ] , blist [ rand() % fac ] ) ;
           blist [ fac ++ ] = ladd [ added ] ;
           blist [ fac ++ ] = basnod ++ ;
           continue ; } 
        randomlist ( fac , inlist , 1 ) ;
        clip = ladd [ added ] ; 
        fullpass_all ( clip , clip ) ; 
        lim = lenlimit = VERYBIG ;
        if ( verbose )
           fprintf ( stderr , "\rRep %3i        %4i taxa added" , reportreplication , added ) ; 
        for ( a = 0 ; a < fac ; ++ a ) {
             tgt = blist [ inlist [ a ] ] ; 
             diff = scorediff_all ( clip , clip , tgt ) ;
             if ( diff < lim ) {
                 lim = diff ;
                 lenlimit = lim - 1 ; 
                 bestpos = tgt ; }}
        nodzert ( clip , bestpos ) ;
        blist [ fac ++ ] = clip ;
        blist [ fac ++ ] = basnod ++ ; } 
   len = downpass_all () ;
   ranseed = seedwas ; 
   return len ; 
}

unsigned long long bswap_all ( void )
{
    unsigned long int diff , lim , nod ;
    unsigned long long len , achk ; 
    int clip , bestpos , begpos , sisis , nunos , bestroot , begroot , i , isdone ;
    int curroot , tgt , rounds , maxrounds = 0 , subsdone = 0 ; 
    short int * rdo , * rdone , * tgtdo , * tgtdone ;
    double prop = 0 , secs ;
    static int clockinit = 0 ;    
    static clock_t last , now ;
    BremTyp * brempt ; 
    anx [ nt ] = nodz ;
    binmksis ( anx ) ;
    if ( !rig [ nt ] ) {
        rig [ nt ] = lef [ nt ] ;
        lef [ nt ] = 0 ; }
    nunos = binoptlist ( nodlist , anx ) ;
    len = downpass_all ( ) ;
    clip = 0 ;
    if ( !clockinit ) {
        clockinit = 1 ; 
        last = clock () ; }
    if ( ratcycles && verbose ) 
         prop = ( double ) ratcyclesdone / ratcycles ; 
    for ( rounds = 1 ; rounds < nodz ; ++ rounds ) {
        if ( ++ clip >= nodz ) clip = 1 ; 
        if ( clip == nt ) continue ; 
        if ( rounds > maxrounds ) maxrounds = rounds ;
        now = clock () ;
        secs = ( double ) ( now  - last ) / CLOCKS_PER_SEC ; 
        if ( secs > reportperiod && !ratcheting && !doingbremer ) {
             last = clock () ; 
             prop = ( double ) ( maxrounds + 1 ) / nodz ;
             if ( ratcycles ) 
                 prop = ( double ) ratcyclesdone / ratcycles ; 
             if ( recordgoodmoves )
                fprintf ( stderr , "%cRep %3i  collapse LEN %-10lli  %12lli rearrangs" , reportheader , reportreplication , len , rearrangs ) ;
             else 
                fprintf ( stderr , "%cRep %3i  %6.2f%%  LEN %-10lli  %12lli rearrangs" , reportheader , reportreplication , prop * 100 , len , rearrangs ) ; }
        unzert ( clip ) ; 
        begpos = bestpos = sisis = sis [ clip ] ;
        if ( sisis ) nodtoskip = anx [ clip ] ;
        fullpass_all ( clip , sisis ) ; 
        lenlimit = VERYBIG ;
        if ( doingbremer ) {
            recorddn = changewts = 1 ;
            gup = gmix = 0 ; }
        if ( clip < nt ) 
             lim = scorediff_all ( clip , clip , sisis ) ;
        else 
           if ( lef [ clip ] > nt ) 
                lim = scorediff_all ( clip , lef [ clip ] , sisis ) ;
           else lim = scorediff_all ( clip , rig [ clip ] , sisis ) ;
        if ( doingbremer ) {
            brempt = bremdif ; 
            recorddn = 2 ; 
            nitdn = gmix ; }
        numcolmoves = numcolroots = 0 ;
        if ( !doingbremer ) 
           if ( ratcheting || recordgoodmoves )
                lenlimit = lim ; 
           else lenlimit = lim - 1 ; 
        if ( lim == 0 && !( ratcheting || recordgoodmoves || doingbremer ) ) {
            zert ( clip , bestpos ) ;
            continue ; }
        tgtdo = tgtdone = tgtlist ;
        if ( !sisis ) * tgtdo ++ = 0 ;
        else {
            * tgtdo ++ = sis [ 0 ] ;
            while ( tgtdone < tgtdo ) {
                nod = * tgtdone ++ ;
                if ( nod > nt ) {
                    * tgtdo ++ = lef [ nod ] ; 
                    * tgtdo ++ = rig [ nod ] ; }}}
        rdo = rdone = rlist ;
        if ( clip < nt ) * rdo ++ = begroot = bestroot = clip ;
        else {
            * rdo ++ = begroot = bestroot = rig [ clip ] ;
            if ( lef [ clip ] > nt ) {
                * rdo ++ = lef [ lef [ clip ] ] ;  
                * rdo ++ = rig [ lef [ clip ] ] ; }}
        while ( rdone < rdo ) {
             curroot = * rdone ;
             if ( curroot > nt ) {
                 * rdo ++ = lef [ curroot ] ; 
                 * rdo ++ = rig [ curroot ] ; }
             for ( tgtdone = tgtlist ; tgtdone < tgtdo ; ++ tgtdone ) {
                  tgt = * tgtdone ; 
                  ++ rearrangs ;
                  if ( verbose )
                     if ( !ratcheting && !recordgoodmoves && !doingbremer ) 
                          fprintf ( stderr , "\rRep %3i  %6.2f%%  LEN %-10lli  %12lli rearrangs" , reportreplication , prop * 100 , len , rearrangs ) ; 
                  gmix = gup = 0 ; 
                  diff = scorediff_all ( clip , curroot , tgt ) ;
                  if ( doingbremer ) {
                     gdn = nitdn - gmix ;
                     brempt -> up = gup ;
                     brempt -> dn = gdn ;
                     brempt -> root = curroot ;
                     brempt -> tgt = tgt ; 
                     ++ brempt ; 
                     continue ; }
                  if ( recordgoodmoves && diff == lim ) {
                      if ( rdone > rlist ) {
                          for ( i = isdone = 0 ; i < numcolroots && !isdone ; ++ i )
                                 if ( colroots [ i ] == curroot ) isdone = 1 ;
                          if ( !isdone ) colroots [ numcolroots ++ ] = curroot ; }
                      for ( i = isdone = 0 ; i < numcolmoves && !isdone ; ++ i )
                                 if ( colmoves [ i ] == tgt ) isdone = 1 ;
                      if ( !isdone ) colmoves [ numcolmoves ++ ] = tgt ; }
                  if ( ratcheting && diff == lim ) {
                       bestpos = tgt ;
                       bestroot = curroot ; }
                  if ( diff < lim ) {
                        len -= lim - diff ; 
                        lim = diff ;
                        lenlimit = lim - 1 ;
                        if ( ratcheting ) ++ lenlimit ; 
                        bestpos = tgt ;
                        bestroot = curroot ; }
                }  // targets 
             if ( use_spr ) break ; 
             ++ rdone ; }   // rootings 
        if ( doingbremer )
            transfer_bremer ( clip , sisis , brempt ) ; 
        if ( recordgoodmoves )
            collapse_with_swaps ( clip , sisis ) ; 
        if ( clip > nt ) 
            if ( !use_spr && bestroot != begroot ) tbreroot ( bestroot , clip ) ;
        zert ( clip , bestpos ) ;
        if ( bestpos != begpos || bestroot != begroot ) {
            rounds = 0 ;
            if ( ratcheting && ++ subsdone == numsubs ) break ; }} // clippings 
    achk = len ; 
    len = downpass_all ( ) ;
    if ( len != achk && !doingbremer ) {
        fprintf ( stderr , "\nOOPS!!!\nAn internal error occured!!!\nlength is %lli, not %lli\n" , len , achk ) ;
        exit ( 0 ) ; }
    return len ;
}

/**********  WAGNER AND BRANCH-SWAPPING WITH DISK BUFFER ****************/

void fullpass_split ( int clip , int sisis )
{
    int nunos , a , nod , c ;
    unsigned long int * lp , * rp , * np , * ap , x , y , z , bak ; 
    nunos = divoptlist ( nodlist , clip , sisis ) ;
    for ( a = nunos ; a -- ; ) {
        nod = nodlist [ a ] ;
        if ( nodtoskip == nod ) continue ; 
        lp = mat [ lef [ nod ] ] ; 
        rp = mat [ rig [ nod ] ] ; 
        np = mat [ nod ] ; 
        for ( c = 0 ; c < piecesz ; ++ c ) {
            x = * lp & * rp ;
            y = * lp ++ | * rp ++ ;
            z = HI4 & ~( x | ( ( x & LO4 ) + LO4 ) ) ;
            z |= z - ( z >> 3 ) ;
            * np ++ = x | ( z & y ) ; }}
    for ( a = 0 ; ++ a < nunos ; ) {
        nod = nodlist [ a ] ;
        if ( nod == clip || anx [ nod ] == nodtoskip || nod == nodtoskip ) continue ; 
        lp = mat [ lef [ nod ] ] - 1 ; 
        rp = mat [ rig [ nod ] ] - 1 ; 
        np = mat [ nod ] - 1 ;
        ap = mat [ anx [ nod ] ] - 1 ;
        for ( c = 0 ; c < piecesz ; ++ c ) {
            x = * ++ lp & * ++ rp ;
            z = * lp | * rp ; 
            y =  HI4 & ~( x | ( ( x & LO4 ) + LO4 ) ) ; 
            y |= y - ( y >> 3 ) ; 
            bak = z | y  ;
            x = * ++ ap & ~ * ++ np ;
            x = ( x | ( ( x & LO4 ) + LO4 ) ) & HI4 ;
            x |= x - ( x >> 3 ) ;
            * np = ( * ap & ~ x ) | ( x & ( * np | * ap & bak ) ) ; }}
    nodtoskip = 0 ; 
} 

unsigned long int scorediff_split ( short int clipis , short int srcis , short int tgtis , unsigned long int atnc )
{
    int len = 0 ;
    unsigned long int * tgt , * tgta , * src , * srca ;
    unsigned long int untgt , unsrc ;
    unsigned long int x , z , y ; 
    unsigned long int * wts ; 
    long int a ;
    if ( clipis < nt ) srca = mat [ clipis ] - 1 ;
    else srca = mat [ anx [ srcis ] ] - 1 ;
    src = mat [ srcis ] - 1 ;
    if ( !tgtis ) tgta = mat [ 0 ] - 1 ;
    else tgta = mat [ anx [ tgtis ] ] - 1 ;
    tgt = mat [ tgtis ] - 1 ;
    if ( ratcycles ) wts = wtchg + atnc - 1 ;
    else wts = wtchg - 1 ;
    for ( a = 0 ; a < piecesz ; ++ a ) {
        unsrc = * ++ src | * ++ srca ; 
        if ( srcis < nt && clipis > nt ) {
           x = * srca & ~ * src ;
           x = ( x | ( ( x & LO4 ) + LO4 ) ) & HI4 ;
           x |= x - ( x >> 3 ) ;
           unsrc = * srca | ( * src & x ) ; }
        untgt = * ++ tgt | * ++ tgta ;
        if ( tgtis && tgtis < nt ) {
           x = * tgta & ~ * tgt ;
           x = ( x | ( ( x & LO4 ) + LO4 ) ) & HI4 ;
           x |= x - ( x >> 3 ) ;
           untgt = * tgta | ( * tgt & x ) ; } 
        x = unsrc & untgt ;
        z = HI4 & ~( x | ( ( x & LO4 ) + LO4 ) ) ;
        if ( changewts ) 
            if ( recorddn ) {
               ++ wts ;
               if ( recorddn == 1 ) * wts = z ;
               x = z & * wts ; 
               y = x ^ z ; 
#ifdef BITS64
               gmix += bitct [ 65535 & ( x | ( x >> 17 ) | ( x >> 34 ) | ( x >> 51 ) ) ] ;
               gup += bitct [ 65535 & ( y | ( y >> 17 ) | ( y >> 34 ) | ( y >> 51 ) ) ] ;
#else 
               gmix += bitct [ 65535 & ( x | ( x >> 17 ) ) ] ;
               gup += bitct [ 65535 & ( y | ( y >> 17 ) ) ] ;
#endif
               continue ; }
            else 
               z &= * ++ wts ; 
#ifdef BITS64
        len += bitct [ 65535 & ( z | ( z >> 17 ) | ( z >> 34 ) | ( z >> 51 ) ) ] ;
#else 
        len += bitct [ 65535 & ( z | ( z >> 17 ) ) ] ;
#endif
        }
    return len ; 
}

unsigned long long wagner_split ( void )
{
    int a , added , fac , basnod , bestpos , clip , tgt ;
    unsigned long int lim , diff , atnc , dapos ;
    unsigned long long len ; 
    short int * inlist = nodfork ;
    unsigned long int * mp ;
    unsigned long int seedwas = ranseed ; 
    int lastpartdone = -1 , atpart = 0 ; 
    randomlist( nt , ladd , 0 ) ;
    initnet ( ladd [ 0 ] , ladd [ 1 ] , ladd [ 2 ] ) ; 
    blist [ 0 ] = ladd [ 1 ] ; 
    blist [ 1 ] = ladd [ 2 ] ; 
    blist [ 2 ] = nt + 1;
    added = 2;
    fac = 3 ;
    basnod = nt + 2;
    ranseed = 0 ; 
    while ( ++ added < nt ) {
        if ( random_starts ) {
           nodzert ( ladd [ added ] , blist [ rand() % fac ] ) ;
           blist [ fac ++ ] = ladd [ added ] ;
           blist [ fac ++ ] = basnod ++ ;
           continue ; } 
        if ( verbose )
           fprintf ( stderr , "\rRep %3i        %4i taxa added" , reportreplication , added ) ; 
        randomlist ( fac , inlist , 1 ) ;
        clip = ladd [ added ] ; 
        for ( a = 0 , mp = movelist ; a < fac ; ++ a , ++ mp ) * mp = 0 ; 
        for ( atnc = 0 ; atnc < nc4 ; atnc += piecesz ) {
            if ( atpart != lastpartdone )  // if loaded already, don't load again... 
                loadmat ( atpart ) ;
            lastpartdone = dapos = atpart ;
            atpart += piecesz ; 
            if ( atpart >= nc4 ) atpart = 0 ; 
            fullpass_split ( clip , clip ) ; 
            for ( a = 0 , mp = movelist ; a < fac ; ++ a , ++ mp ) {
                 tgt = blist [ inlist [ a ] ] ; 
                 * mp += scorediff_split ( clip , clip , tgt , dapos ) ; }}
        atpart = lastpartdone ; 
        lim = VERYBIG ;
        for ( a = 0 , mp = movelist ; a < fac ; ++ a , ++ mp ) {
             diff = * mp ; 
             if ( diff < lim ) {
                 lim = diff ;
                 bestpos = blist [ inlist [ a ] ] ; }}
        nodzert ( clip , bestpos ) ;
        blist [ fac ++ ] = clip ;
        blist [ fac ++ ] = basnod ++ ; } 
   len = downpass_split () ;
   ranseed = seedwas ; 
   return len ; 
}

unsigned long long bswap_split ( void )
{
    unsigned long int diff , lim , nod , atnc , j , initrearrangs , dapos ;
    unsigned long long len , achk ; 
    int clip , bestpos , begpos , sisis , nunos , bestroot , begroot , i , isdone ;
    int curroot , tgt , rounds , maxrounds = 0 , subsdone = 0 ; 
    short int * rdo , * rdone , * tgtdo , * tgtdone ;
    double prop = 0 , secs , curprop ;
    static int clockinit = 0 ;
    unsigned long int * mp ; 
    static clock_t last , now ;
    int lastpartdone = -1 , atpart = 0 ; 
    BremTyp * brempt ; 
    anx [ nt ] = nodz ;
    binmksis ( anx ) ;
    if ( !rig [ nt ] ) {
        rig [ nt ] = lef [ nt ] ;
        lef [ nt ] = 0 ; }
    nunos = binoptlist ( nodlist , anx ) ;
    len = downpass_split ( ) ;
    clip = 0 ;
    if ( !clockinit ) {
        clockinit = 1 ; 
        last = clock () ; }
    if ( ratcycles && verbose ) 
         prop = ( double ) ratcyclesdone / ratcycles ; 
    for ( rounds = 1 ; rounds < nodz ; ++ rounds ) {
        if ( ++ clip >= nodz ) clip = 1 ; 
        if ( clip == nt ) continue ; 
        if ( rounds > maxrounds ) maxrounds = rounds ;
        now = clock () ;
        secs = ( double ) ( now  - last ) / CLOCKS_PER_SEC ; 
        if ( secs > reportperiod && !ratcheting && !doingbremer ) {
             last = clock () ; 
             prop = ( double ) ( maxrounds + 1 ) / nodz ;
             if ( ratcycles ) 
                 prop = ( double ) ratcyclesdone / ratcycles ; 
             if ( recordgoodmoves )
                fprintf ( stderr , "%cRep %3i  collapse LEN %-10lli  %12lli rearrangs" , reportheader , reportreplication , len , rearrangs ) ;
             else 
                fprintf ( stderr , "%cRep %3i  %6.2f%%  LEN %-10lli  %12lli rearrangs" , reportheader , reportreplication , prop * 100 , len , rearrangs ) ; }
        unzert ( clip ) ; 
        begpos = bestpos = sisis = sis [ clip ] ;
        if ( sisis ) nodtoskip = anx [ clip ] ;
        initrearrangs = rearrangs ; 
        /***  Make lists of rootings and targets *********/
        tgtdo = tgtdone = tgtlist ;
        if ( !sisis ) * tgtdo ++ = 0 ;
        else {
            * tgtdo ++ = sis [ 0 ] ;
            while ( tgtdone < tgtdo ) {
                nod = * tgtdone ++ ;
                if ( nod > nt ) {
                    * tgtdo ++ = lef [ nod ] ; 
                    * tgtdo ++ = rig [ nod ] ; }}}
        rdo = rdone = rlist ;
        if ( clip < nt ) * rdo ++ = begroot = bestroot = clip ;
        else {
            * rdo ++ = begroot = bestroot = rig [ clip ] ;
            if ( lef [ clip ] > nt ) {
                * rdo ++ = lef [ lef [ clip ] ] ;  
                * rdo ++ = rig [ lef [ clip ] ] ; }}
        while ( rdone < rdo ) {
             curroot = * rdone ++ ;
             if ( curroot > nt ) {
                 * rdo ++ = lef [ curroot ] ; 
                 * rdo ++ = rig [ curroot ] ; }}
        /***   Calculate scores for each move, one matrix piece at a time *****/
        lim = 0 ;
        for ( atnc = 0 ; atnc < nc4 ; atnc += piecesz ) {
            if ( atpart != lastpartdone )  // if loaded already, don't load again... 
                loadmat ( atpart ) ;  
            lastpartdone = dapos = atpart ;
            atpart += piecesz ; 
            if ( atpart >= nc4 ) atpart = 0 ; 
            fullpass_split ( clip , sisis ) ; 
            if ( doingbremer ) {
               recorddn = changewts = 1 ;
               gup = gmix = 0 ; }
            lim += scorediff_split ( clip , * rlist , sisis , dapos ) ;
            if ( doingbremer ) {
                brempt = bremdif ; 
                recorddn = 2 ; 
                nitdn = gmix ; }
            mp = movelist ; 
            rdone = rlist ;
            if ( verbose )
              if ( !ratcheting && !recordgoodmoves && !doingbremer ) {
                 i = ( rdo - rlist ) * ( tgtdo - tgtlist ) ; 
                 curprop = ( double ) atnc / nc4 ; 
                 j = ( i * curprop ) + initrearrangs ; 
                 fprintf ( stderr , "\rRep %3i  %6.2f%%  LEN %-10lli  %12li rearrangs" , reportreplication , prop * 100 , len , j ) ; }
            while ( rdone < rdo ) {
                 curroot = * rdone ++ ;
                 for ( tgtdone = tgtlist ; tgtdone < tgtdo ; ++ tgtdone ) {
                      tgt = * tgtdone ; 
                      if ( !atnc ) {
                          ++ rearrangs ;
                          * mp = 0 ;
                          if ( doingbremer ) {
                             brempt -> up = 0 ;
                             brempt -> dn = 0 ; }}
                      gmix = gup = 0 ;
                      * mp += scorediff_split ( clip , curroot , tgt , dapos ) ;
                      if ( doingbremer ) {
                          gdn = nitdn - gmix ;
                          brempt -> up += gup ;
                          brempt -> dn += gdn ;
                          brempt -> root = curroot ;
                          brempt -> tgt = tgt ; 
                          ++ brempt ; }
                      ++ mp ; }  // targets 
                 if ( use_spr ) break ; }  // rootings 
             }  // matrix pieces
        atpart = lastpartdone ; 
        /*** now repeat rootings & targets, using total scores... *********/
        numcolmoves = numcolroots = 0 ; 
        mp = movelist ; 
        rdone = rlist ;
        if ( doingbremer ) rdone = rdo ; 
        while ( rdone < rdo ) {
             curroot = * rdone ++ ;
             for ( tgtdone = tgtlist ; tgtdone < tgtdo ; ++ tgtdone ) {
                      tgt = * tgtdone ; 
                      diff = * mp ++ ;
                      if ( recordgoodmoves && diff == lim ) {
                          if ( rdone > rlist ) {
                              for ( i = isdone = 0 ; i < numcolroots && !isdone ; ++ i )
                                     if ( colroots [ i ] == curroot ) isdone = 1 ;
                              if ( !isdone ) colroots [ numcolroots ++ ] = curroot ; }
                          for ( i = isdone = 0 ; i < numcolmoves && !isdone ; ++ i )
                                     if ( colmoves [ i ] == tgt ) isdone = 1 ;
                          if ( !isdone ) colmoves [ numcolmoves ++ ] = tgt ; }
                      if ( ratcheting && diff == lim ) {
                           bestpos = tgt ;
                           bestroot = curroot ; }
                      if ( diff < lim ) {
                            len -= lim - diff ; 
                            lim = diff ;
                            bestpos = tgt ;
                            bestroot = curroot ; }}  // targets
             if ( use_spr ) break ; }  // rootings 
        if ( doingbremer )
            transfer_bremer ( clip , sisis , brempt ) ; 
        if ( recordgoodmoves )
            collapse_with_swaps ( clip , sisis ) ; 
        if ( clip > nt ) 
            if ( !use_spr && bestroot != begroot ) tbreroot ( bestroot , clip ) ;
        zert ( clip , bestpos ) ;
        if ( bestpos != begpos || bestroot != begroot ) {
            rounds = 0 ;
            if ( ratcheting && ++ subsdone == numsubs ) break ; }} // clippings 
    achk = len ;
    if ( !doingbremer ) {
       len = downpass_split ( ) ;
       if ( len != achk ) {
           fprintf ( stderr , "\nOOPS!!!\nAn internal error occured!!!\nlength is %lli, not %lli\n" , len , achk ) ;
           exit ( 0 ) ; }}
    return len ;
}

/*******  CALCULATE (AND POSSIBLY COMBINE) BREMER SUPPORTS *****************/

void calcusizes ( short int * to , short int * daprvtree , short int * inthtulist )
{
  int a , b , nodn , i ;
  short int * lp ; 
  for ( a = nt ; a -- ; ) {
     inthtulist [ a ] = a ; 
     daprvtree [ a ] = 1 ; }
  binmksis ( to ) ;
  binoptlist ( nodlist , to ) ;
  nodn = nodz - nt ;
  lp = nodlist + nodn;
  for( i = nodn ; i -- ; ) {
      a = * ( -- lp ) ;
      daprvtree [ a ] = 0 ;
      inthtulist [ a ] = -1 ; 
      b = lef [ a ] ; 
      daprvtree [ a ] += daprvtree [ b ] ;
      if ( inthtulist [ b ] >= 0 )
           inthtulist [ a ] = inthtulist [ b ] ;
      b = rig [ a ] ; 
      daprvtree [ a ] += daprvtree [ b ] ;
      if ( inthtulist [ b ] >= 0 )
           inthtulist [ a ] = inthtulist [ b ] ; }
}

void do_node_equivalences ( short int * from , short int * to , short int * save )
{
  int b , c , d , e , i ;
  int mysz ;
  short int * dolist = ladd , * done , * todo , * lp , * nodok = lassis ;
  int isitok ;
  short int * conslist = tgtlist , consnodz , * tmpsize = rlist ;
  short int * bufprvtree = blist , * thet = from ;
  short int * waglist = unsuprt , * htulist = colmoves  ; 
  calcusizes ( to , bufprvtree , htulist ) ;
  consnodz = nodz - nt ; 
  for ( i = 0 ; i < consnodz ; ++ i ) conslist [ i ] = nodlist [ i ] ; 
  calcusizes ( thet , tmpsize , waglist ) ;
  for ( i = nt ; i < nodz ; ++ i ) 
     nodok [ i ] = save [ i ] = -1 ; 
  lp = conslist + consnodz ;
  for( i = consnodz ; i -- ; ) {
       b = * ( -- lp ) ;
       c = htulist [ b ] ;
       if ( c < 0 ) continue ;
       if ( ! thet [ c ] ) continue ; 
       mysz = bufprvtree [ b ] ;
       while ( tmpsize [ c ] < mysz && c != nt ) c = thet [ c ] ;
       if ( tmpsize [ c ] != mysz ) continue ; 
       * ( todo = dolist ) = c ;
       done = todo ++ ; 
       isitok = 1 ; 
       while ( done < todo ) {
           d = * done ++ ;
           if ( d >= nt ) {
             if ( nodok [ d ] >= 0 ) {  
               for ( e = nodok [ d ] ; e != nt && bufprvtree [ e ] <= mysz && e != b ; e = to [ e ] ) ;
               if ( e != b ) { 
                 isitok = 0 ; 
                 break ; }}
             else {
                 * todo ++ = lef [ d ] ;
                 * todo ++ = rig [ d ] ; }}
           else {
              for ( e = d ; e != nt && bufprvtree [ e ] <= mysz && e != b ; e = to [ e ] ) ;
              if ( e != b ) {
                  isitok = 0 ;
                  break ; }}}
       if ( isitok ) {
           save [ b ] = c ; 
           nodok [ c ] = htulist [ b ] ; }}
}

void transfer_bremer ( int clip , int sisis , BremTyp * hidest )
{
    BremTyp * at = bremdif - 1 ;
    int a , b , comnod , rootis ;
    int da , db , valisneg ;
    unsigned long int mid ; 
    double abs , rel , eval , jprob ;
    eval = ( double ) jakprob / 100 ;
    jprob = 1 - eval ;
    zert ( clip , sisis ) ;
    while ( ++ at < hidest ) {
        valisneg = 0 ; 
        if  ( shownegbremers ) 
           if ( at -> up < at -> dn ) {
               mid = at -> up ;
               at -> up = at -> dn ;
               at -> dn = mid ;
               valisneg = 1 ; }
        abs = at -> up - at -> dn ;
        if ( abs > 0 ) {
           rel = abs / at -> up ;
           if ( bremertype == 2 ) 
                eval = pow ( jprob * rel , 1 / abs ) ; 
           else if ( bremertype == 1 ) eval = rel ;
                else eval = abs ; }
        else eval = 0 ;
        if ( valisneg ) eval = -eval ; 
        rootis = at -> root ; 
        if ( !use_spr && clip > nt && anx [ rootis ] != clip ) {
            a = anx [ rootis ] ;
            while ( a != clip ) {
               if ( tbrem [ a ] > eval ) tbrem [ a ] = eval ; 
               a = anx [ a ] ; }}
        if ( at -> tgt == sisis ) continue ; 
        da = distoroot [ a = clip ] ; 
        db = distoroot [ b = at -> tgt ] ;
        while ( distoroot [ a ] > db ) a = anx [ a ] ;
        while ( distoroot [ b ] > da ) b = anx [ b ] ;
        while ( a != b ) { a = anx [ a ] ; b = anx [ b ] ; }
        comnod = a ;
        if ( comnod != at -> tgt ) 
            for ( a = anx [ at -> tgt ] ; a != comnod ; a = anx [ a ] ) {
                if ( abrem [ a ] > abs ) abrem [ a ] = abs ;
                if ( abrem [ a ] < abs ) continue ; // only consider moves within absolute BS
                if ( tbrem [ a ] > eval ) tbrem [ a ] = eval ;  }
        for ( b = anx [ clip ] ; b != comnod ; b = anx [ b ] ) {
           if ( abrem [ a ] > abs ) abrem [ a ] = abs ;
           if ( abrem [ a ] < abs ) continue ; // ditto
           if ( tbrem [ b ] > eval ) tbrem [ b ] = eval ; }}
    unzert ( clip ) ; 
}

void do_quick_bremers ( void )
{
    int i , j , is ;
    short int * anxwas = anx ;
    long long thisdiff ; 
    doingbremer = 1 ;
    rearrangs = 0 ;
    if ( reporttoscreen ) 
       fprintf ( stderr , "\nCalculating supports from %i trees..." , numtrees ) ; 
    for ( i = 0 ; i < numtrees ; ++ i ) {
      for ( j = nt ; j < nodz ; ++ j ) {
          tbrem [ j ] = VERYBIG ; 
          abrem [ j ] = 2000000000 ; }
        anx = memtrees [ i ] ; 
        filldistoroot ( ) ;
        bswap ( ) ;
        if ( bremertype == 0 ) 
           thisdiff = downpass () - bestgloballen ; 
        if ( !i ) {
           for ( j = nt + 1 ; j < nodz ; ++ j ) 
              gbrem [ j ] = tbrem [ j ] ; }
        else {
            do_node_equivalences ( memtrees [ i ] , memtrees [ 0 ] , nodmatch ) ;
            for ( j = nt + 1 ; j < nodz ; ++ j ) {
                is = nodmatch [ j ] ;
                if ( is < 0 ) {
                    if ( bremertype == 0 ) {
                        if ( gbrem [ j ] > thisdiff ) gbrem [ j ] = thisdiff ; }
                    else 
                       gbrem [ j ] = 0 ;
                    continue ; }
                if ( gbrem [ j ] > tbrem [ is ] )
                    gbrem [ j ] = tbrem [ is ] ; }}}
    fprintf ( stderr , "%cCalculated supports, %lli additional rearrangs.\n" , reportheader , rearrangs ) ; 
    anx = anxwas ; 
}

/**************** HANDLING TREE SEARCHES ***********************/

char getatree ( void  )
{
    int htu , atnode , i ;
    char car = '(' ;
    for( htu = nt ; htu -- ; ) anx [ htu ] = 0 ; 
    htu = atnode = nt - 1 ; 
    while ( car != '*' && car != ';' ) {
       if ( car == '(' )
           { errout ( htu > nodz , "Tree contains too many open parenths ... can't read it" ) ; 
             anx [ htu + 1 ] = atnode ;
             atnode = ++ htu ; } 
       else if ( car == ')' ) {
                errout ( atnode == nodz , "Tree contains too many close parenths ... can't read it" ) ; 
                atnode = anx [ atnode ] ; }
             else {
                 ungetc ( car , treeinp ) ; 
                 fscanf ( treeinp , " %i" , &i ) ; 
                 errout ( anx [ i ] , "Duplicate taxa in tree" ) ; 
                 anx [ i ] = atnode ; }
       while ( isspace ( car = getc ( treeinp ) ) ) ; }
    for( i = nt ; i -- ; )
        errout ( anx [ i ] == 0 , "Tree is incomplete ... can't read it\nAre you sure it corresponds to the same matrix?" ) ; 
    anx [ nt ] = ++ htu ;
    errout ( htu != nodz , "Tree is not binary ... can't process it" ) ;
    binmksis ( anx ) ; 
    return car ;
} 

unsigned long long doratchet ( void )
{
    int i , j , numtreeswas = numtrees ;
    unsigned long long len , beslen = downpass ( ) ;
    numsubs = nt / 2 ;
    for ( i = 0 ; i < nodz ; ++ i )
       memtrees [ numtrees ] [ i ] = anx [ i ] ; 
    ratcyclesdone = 0 ;
    if ( dojak ) storewtfilter () ;
    else
      if ( beslen <= bestgloballen || saveall ) {
           bestgloballen = beslen ; 
           storetree () ; }
    for ( j = 0 ; j < ratcycles ; ++ j ) {
        dowtfilter ( 5 , dojak , nc4 ) ;
        ratcheting = changewts = 1 ;
        ratcyclesdone = j + 1 ; 
        bswap () ;
        if ( dojak ) resetwtfilter ( ) ; 
        ratcheting = 0 ; 
        changewts = dojak ;
        len = bswap () ;
        if ( !dojak ) {
           if ( len < beslen ) beslen = len ; 
           if ( len < bestgloballen ) {
              if ( !saveall ) numtrees = 0 ;
              else numtrees = numtreeswas ; 
              numhits = 0 ; 
              bestgloballen = len ; }
           if ( len == bestgloballen ) {
               ++ numhits ; 
               storetree () ; }}
        else 
          if ( len < beslen ) {
             beslen = len ;
             for ( i = 0 ; i < nodz ; ++ i )
                memtrees [ numtrees ] [ i ] = anx [ i ] ; }}
    ratcyclesdone = 0 ; 
    if ( dojak ) 
       for ( i = 0 ; i < nodz ; ++ i )
          anx [ i ] = memtrees [ numtrees ] [ i ] ; 
    return beslen ; 
}

void dosearches ( void )
{
    int a , i , exists , nxt ;
    unsigned long long len ;
    clock_t beg , end ;
    double secs ;
    beg = clock () ; 
    bestgloballen = VERYBIG ;
    for ( a = 0 ; a ++ < numreplications ; ) {
        srand ( ranseed ) ; 
        ranseed = rand () ;
        reportreplication = a ;
        if ( dojak && ( together || ratcycles ) ) 
            dowtfilter ( jakprob , 0 , nc4 ) ; 
        if ( treeinp == NULL ) 
              len = wagner () ;
        else {
            while ( isspace ( i = getc ( treeinp ) ) ) ;
            errout ( i != '(' , "Expected an opening parent in tree-file" ) ;
            nxt = getatree ( ) ;
            len = downpass ( ) ; }
        if ( do_branch_swapping ) {
            reset_nodz ( anx ) ;
            exists = 0 ; 
            if ( !dojak ) 
              for ( i = 0 ; i < numtrees && !exists ; ++ i )
                  if ( equal_trees ( anx , memtrees [ i ] ) ) exists = 1 ; 
            if ( !exists ) 
                len = bswap ( ) ; }
        if ( ratcycles )
           len = doratchet ( ) ; 
        fprintf ( stderr , "%cRep %3i           LEN %-10lli  %12lli rearrangs" , reportheader , reportreplication, len , rearrangs ) ;
        if ( reporttoscreen ) fprintf ( stderr , "\n" ) ; 
        if ( !dojak ) 
          if ( len < bestgloballen ) {
              if ( !saveall ) numtrees = 0 ;
              numhits = 0 ; 
              bestgloballen = len ; }
        if ( dojak ) collapse_the_tree ( ) ;
        if ( !ratcycles || dojak )  // jackknifing saves a single tree per ratchet set (sorry!)
           if ( len == bestgloballen || dojak || saveall ) {
               ++ numhits ; 
               storetree () ; }
        if ( treeinp != NULL && nxt == ';' ) break ; }
    if ( numreplications > 1 && !dojak )
         fprintf ( stderr , "\n             BEST LEN %lli, HIT %i TIMES \n" , bestgloballen , numhits ) ; 
    else fprintf ( stderr , "\n" ) ; 
    if ( quickbremer )
        do_quick_bremers ( ) ; 
    if ( !together ) {
        fclose ( partsfile ) ; 
        remove ( "oblong.pts" ) ; }
    end = clock () ;
    secs = ( double ) ( end - beg ) / CLOCKS_PER_SEC ; 
    fprintf ( stderr , "\nTime used in tree calculations: %.2f secs.\n" , secs ) ;
}     

/*******  MAIN --where it all begins... *****************/

int main ( int argc , char ** argv )
{
  procargs ( argc , argv ) ;
  init_bit_counters ( ) ; 
  read_input ( ) ;
  dosearches () ;
  savetrees () ;
  fclose ( outp ) ;
  return 1 ; 
}