From 1fdcf365eddbbd433de685df4cb037e2f6b2d402 Mon Sep 17 00:00:00 2001 From: Kuoi Date: Mon, 7 Mar 2022 20:38:24 +0000 Subject: [PATCH] init --- oblong.c | 2456 +++++++++++++++++++++++++++++++++++++++++++++++++++++ readme.md | 32 + 2 files changed, 2488 insertions(+) create mode 100644 oblong.c create mode 100644 readme.md diff --git a/oblong.c b/oblong.c new file mode 100644 index 0000000..25bbe15 --- /dev/null +++ b/oblong.c @@ -0,0 +1,2456 @@ +/************************************************************************************* + * 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 +#include +#include +#include +#include +#include +#include +#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 ) ; + exit ( 0 ) ; +} + +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 ; } + if ( * a || * b ) return 0 ; + return 1 ; +} + +void procargs ( int nargs , char ** arg ) +{ + int i ; + char * cp ; + if ( nargs == 1 ) { + display_help ( 0 ) ; + exit ( 0 ) ; } + for ( i = 1 ; i < nargs ; ++ i ) { + cp = arg [ i ] ; + if ( * cp != '-' ) { + fprintf ( stderr , "\nWrong argument: %s\n" , arg [ i ] ) ; + exit ( 0 ) ; } + switch ( * ++ cp ) { + case 'a': saveall = 1 ; break ; + case 'b': do_branch_swapping = 0 ; break ; + case 'e': equal_taxon_sequence = 1 ; break ; + case 'f': fasta = 1 ; break ; + case 'h': display_help ( 1 ) ; exit ( 0 ) ; + case 'i': + case '@': + if ( * cp == '@' ) multifiles = 1 ; + ++ 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" ) ; + break ; + case 'j': + ++ cp ; + if ( isdigit ( * cp ) ) + jakprob = atoi ( cp ) ; + dojak = changewts = 1 ; + break ; + case 'k': swapcollapse = 1 ; break ; + case 'm': errout ( ( ( maxram = atof ( ++ cp ) ) <= 0 ) , "You must indicate RAM to use" ) ; together = 0 ; break ; + case 'N': usenexus = 1 ; break ; + case 'o': + ++ 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 ; + case 'q': + ++ 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 ; + case 'R': ranseed = atoi ( ++ cp ) ; break ; + case 's': + if ( !strcmp ( cp , "spr" ) ) { use_spr = 1 ; break ; } + case 't': savejusttrees = 1 ; break ; + case 'T': + ++ 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 ; + case 'w': random_starts = 1 ; break ; + case 'x': + ++ cp ; + if ( isdigit ( * cp ) ) ratcycles = atoi ( cp ) ; + 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" ) ; + errout ( ( inp == NULL ) , "No data set specified" ) ; + errout ( ( outp == NULL ) , "No output file specified" ) ; + if ( quickbremer && dojak ) { + dojak = changewts = 0 ; + 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 ; + if ( !do_branch_swapping && !ratcycles ) shownegbremers = 1 ; + 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 ; } + 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 ; + for ( a = 0 ; a < 65536 ; ++ a ) { + for ( b = 1 , c = 0 ; b <= a ; b <<= 1 ) + if ( ( b & a ) ) ++ 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 ) ; + if ( i != '>' ) { + fprintf ( stderr , "\nError reading fasta file --expected taxon name\n" ) ; + errout ( 1 , NULL ) ; } + while ( !feof ( inp ) ) { + fscanf ( inp , " %s" , nampt ) ; + i = 0 ; + ++ curnt ; + if ( curnt == 1) + while ( 1 ) { + fscanf ( inp , " %c" , &c ) ; + if ( c == '>' || feof ( inp ) ) break ; + ++ 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 ( ) ; + taxnam = ( char ** ) myalloc ( MAXMULTITAXA * sizeof ( char * ) ) ; + fscanf ( inp , " %s" , nampt ) ; + while ( !feof ( multinp ) ) { + ++ numfiles ; + if ( ( inp = fopen ( nampt , "rb" ) ) == NULL ) { + fprintf ( stderr , "\nCannot open file %s (from list of multiple files)\n" , nampt ) ; + exit ( 0 ) ; } + if ( reporttoscreen ) + fprintf ( stderr , "%cChecking input file %3i %-40s" , reportheader , numfiles , nampt ) ; + if ( fasta ) get_fasta_dims ( &nc , &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 ) ) ; + 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 , " " ) ; + if ( treeinp != NULL ) { + while ( strcmp ( nampt , "tread" ) && !feof ( treeinp ) ) + fscanf ( treeinp , " %7s" , nampt ) ; + 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 ; +} + + diff --git a/readme.md b/readme.md new file mode 100644 index 0000000..aa9c4a0 --- /dev/null +++ b/readme.md @@ -0,0 +1,32 @@ + +Oblong is self-explanatory; just call it without arguments. + +For Mac or Linux, you may have to make the file executable; +open a command shell, go to the directory where oblong-linux +or oblong-mac is, and type "chmod a+x oblong ". + +You can rename binaries; the program does not keep track of +its own name. + + + +-------------------------------------- +Notes for compilation + + +For Linux or Mac: + + gcc -O3 -o oblong -DBITS64 oblong.c + + For linux, you may have to add "/usr/lib/libm.a" + after "oblong.c" (the library containing calls to + some math functions). + +For Windows (Watcom): + + create a project, and define (under Options/C compiler switches) + a macro "WATCOM" under "Source switches" --this is so that + the program will use functions in the Watcom library, enabling + it to handle larger files. + +