diff --git a/CUTAL.c b/CUTAL.c index 33b1afb..e42cebb 100644 --- a/CUTAL.c +++ b/CUTAL.c @@ -1,11 +1,11 @@ -/* CUTAL (version 0.1) a program for cutting alignement into non-recombinant blocks under - * a parsimony framework +/* CUTAL (version 0.1) a program for cutting alignement into non-recombinant blocks under + * a parsimony framework * Copyright May 2018 by Celine Scornavacca and Mark Jones - * - * The program is heavily based on Parsimonator-1.0.2 by Alexandros Stamatakis, which was partially + * + * The program is heavily based on Parsimonator-1.0.2 by Alexandros Stamatakis, which was partially * derived from fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen - * - * and + * + * and * * programs of the PHYLIP package by Joe Felsenstein. * @@ -18,9 +18,9 @@ * 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. - * * - * For any other enquiries send an email to celine.scornavacca@umontpellier.fr and + * + * For any other enquiries send an email to celine.scornavacca@umontpellier.fr and * markelliotlloyd@gmail.com * */ @@ -46,7 +46,7 @@ -#if ! (defined(__ppc) || defined(__powerpc__) || defined(PPC)) +#if !(defined(__ppc) || defined(__powerpc__) || defined(PPC)) #include #endif @@ -61,12 +61,12 @@ FILE *INFILE; -char run_id[128] = "", - seq_file[1024] = "", - tree_file[1024]="", - resultFileName[1024] = "", - infoFileName[1024] = "", - randomFileName[1024] = ""; +char run_id[128] = "", + seq_file[1024] = "", + tree_file[1024]="", + resultFileName[1024] = "", + infoFileName[1024] = "", + randomFileName[1024] = ""; static boolean whitechar (int ch); @@ -74,382 +74,382 @@ static void treeEchoContext (FILE *fp1, FILE *fp2, int n); static stringHashtable *initStringHashTable(hashNumberType n) { - /* - init with primes - */ - - static const hashNumberType initTable[] = {53, 97, 193, 389, 769, 1543, 3079, 6151, 12289, 24593, 49157, 98317, - 196613, 393241, 786433, 1572869, 3145739, 6291469, 12582917, 25165843, - 50331653, 100663319, 201326611, 402653189, 805306457, 1610612741}; - - - /* init with powers of two - - static const hashNumberType initTable[] = {64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, - 32768, 65536, 131072, 262144, 524288, 1048576, 2097152, - 4194304, 8388608, 16777216, 33554432, 67108864, 134217728, - 268435456, 536870912, 1073741824, 2147483648U}; - */ - - stringHashtable *h = (stringHashtable*)malloc(sizeof(stringHashtable)); - - hashNumberType - tableSize, - i, - primeTableLength = sizeof(initTable)/sizeof(initTable[0]), - maxSize = (hashNumberType)-1; - - assert(n <= maxSize); - - i = 0; - - while(initTable[i] < n && i < primeTableLength) - i++; - - assert(i < primeTableLength); - - tableSize = initTable[i]; - - h->table = (stringEntry**)calloc(tableSize, sizeof(stringEntry*)); - h->tableSize = tableSize; - - return h; + /* + init with primes + */ + + static const hashNumberType initTable[] = {53, 97, 193, 389, 769, 1543, 3079, 6151, 12289, 24593, 49157, 98317, + 196613, 393241, 786433, 1572869, 3145739, 6291469, 12582917, 25165843, + 50331653, 100663319, 201326611, 402653189, 805306457, 1610612741}; + + + /* init with powers of two + + static const hashNumberType initTable[] = {64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, + 32768, 65536, 131072, 262144, 524288, 1048576, 2097152, + 4194304, 8388608, 16777216, 33554432, 67108864, 134217728, + 268435456, 536870912, 1073741824, 2147483648U}; + */ + + stringHashtable *h = (stringHashtable*)malloc(sizeof(stringHashtable)); + + hashNumberType + tableSize, + i, + primeTableLength = sizeof(initTable)/sizeof(initTable[0]), + maxSize = (hashNumberType)-1; + + assert(n <= maxSize); + + i = 0; + + while(initTable[i] < n && i < primeTableLength) + i++; + + assert(i < primeTableLength); + + tableSize = initTable[i]; + + h->table = (stringEntry**)calloc(tableSize, sizeof(stringEntry*)); + h->tableSize = tableSize; + + return h; } static hashNumberType hashString(char *p, hashNumberType tableSize) { - hashNumberType h = 0; - - for(; *p; p++) - h = 31 * h + *p; - - return (h % tableSize); + hashNumberType h = 0; + + for(; *p; p++) + h = 31 * h + *p; + + return (h % tableSize); } - + static void addword(char *s, stringHashtable *h, int nodeNumber) { - hashNumberType position = hashString(s, h->tableSize); - stringEntry *p = h->table[position]; - - for(; p!= NULL; p = p->next) - { - if(strcmp(s, p->word) == 0) - return; - } - - p = (stringEntry *)malloc(sizeof(stringEntry)); - - assert(p); - - p->nodeNumber = nodeNumber; - p->word = (char *)malloc((size_t)(strlen(s) + 1) * sizeof(char)); - - strcpy(p->word, s); - - p->next = h->table[position]; - - h->table[position] = p; + hashNumberType position = hashString(s, h->tableSize); + stringEntry *p = h->table[position]; + + for(; p!= NULL; p = p->next) + { + if(strcmp(s, p->word) == 0) + return; + } + + p = (stringEntry *)malloc(sizeof(stringEntry)); + + assert(p); + + p->nodeNumber = nodeNumber; + p->word = (char *)malloc((size_t)(strlen(s) + 1) * sizeof(char)); + + strcpy(p->word, s); + + p->next = h->table[position]; + + h->table[position] = p; } static int lookupWord(char *s, stringHashtable *h) { - hashNumberType position = hashString(s, h->tableSize); - stringEntry *p = h->table[position]; - - for(; p!= NULL; p = p->next) - { - if(strcmp(s, p->word) == 0) - return p->nodeNumber; - } - - return -1; + hashNumberType position = hashString(s, h->tableSize); + stringEntry *p = h->table[position]; + + for(; p!= NULL; p = p->next) + { + if(strcmp(s, p->word) == 0) + return p->nodeNumber; + } + + return -1; } static int treeFinishCom (FILE *fp, char **strp) { - int ch; - - while ((ch = getc(fp)) != EOF && ch != ']') { - if (strp != NULL) *(*strp)++ = ch; /* save character */ - if (ch == '[') { /* nested comment; find its end */ - if ((ch = treeFinishCom(fp, strp)) == EOF) break; - if (strp != NULL) *(*strp)++ = ch; /* save closing ] */ - } - } - - if (strp != NULL) **strp = '\0'; /* terminate string */ - return ch; + int ch; + + while ((ch = getc(fp)) != EOF && ch != ']') { + if (strp != NULL) *(*strp)++ = ch; /* save character */ + if (ch == '[') { /* nested comment; find its end */ + if ((ch = treeFinishCom(fp, strp)) == EOF) break; + if (strp != NULL) *(*strp)++ = ch; /* save closing ] */ + } + } + + if (strp != NULL) **strp = '\0'; /* terminate string */ + return ch; } /* treeFinishCom */ static int treeGetCh (FILE *fp) /* get next nonblank, noncomment character */ { /* treeGetCh */ - int ch; - - while ((ch = getc(fp)) != EOF) { - if (whitechar(ch)) ; - else if (ch == '[') { /* comment; find its end */ - if ((ch = treeFinishCom(fp, (char **) NULL)) == EOF) break; - } - else break; - } - - return ch; + int ch; + + while ((ch = getc(fp)) != EOF) { + if (whitechar(ch)); + else if (ch == '[') { /* comment; find its end */ + if ((ch = treeFinishCom(fp, (char **) NULL)) == EOF) break; + } + else break; + } + + return ch; } /* treeGetCh */ static boolean treeLabelEnd (int ch) { - switch (ch) - { - case EOF: - case '\0': - case '\t': - case '\n': - case '\r': - case ' ': - case ':': - case ',': - case '(': - case ')': - case ';': - return TRUE; - default: - break; - } - return FALSE; -} + switch (ch) + { + case EOF: + case '\0': + case '\t': + case '\n': + case '\r': + case ' ': + case ':': + case ',': + case '(': + case ')': + case ';': + return TRUE; + default: + break; + } + return FALSE; +} static boolean treeGetLabel (FILE *fp, char *lblPtr, int maxlen) { - int ch; - boolean done, quoted, lblfound; - - if (--maxlen < 0) - lblPtr = (char *) NULL; - else - if (lblPtr == NULL) - maxlen = 0; - - ch = getc(fp); - done = treeLabelEnd(ch); - - lblfound = ! done; - quoted = (ch == '\''); - if (quoted && ! done) - { - ch = getc(fp); - done = (ch == EOF); - } - - while (! done) - { - if (quoted) - { - if (ch == '\'') - { - ch = getc(fp); - if (ch != '\'') - break; - } + int ch; + boolean done, quoted, lblfound; + + if (--maxlen < 0) + lblPtr = (char *) NULL; + else + if (lblPtr == NULL) + maxlen = 0; + + ch = getc(fp); + done = treeLabelEnd(ch); + + lblfound = !done; + quoted = (ch == '\''); + if (quoted && !done) + { + ch = getc(fp); + done = (ch == EOF); } - else - if (treeLabelEnd(ch)) break; - if (--maxlen >= 0) *lblPtr++ = ch; - ch = getc(fp); - if (ch == EOF) break; - } + while (!done) + { + if (quoted) + { + if (ch == '\'') + { + ch = getc(fp); + if (ch != '\'') + break; + } + } + else + if (treeLabelEnd(ch)) break; + + if (--maxlen >= 0) *lblPtr++ = ch; + ch = getc(fp); + if (ch == EOF) break; + } - if (ch != EOF) (void) ungetc(ch, fp); + if (ch != EOF) (void) ungetc(ch, fp); - if (lblPtr != NULL) *lblPtr = '\0'; + if (lblPtr != NULL) *lblPtr = '\0'; - return lblfound; + return lblfound; } static boolean treeFlushLabel (FILE *fp) -{ - return treeGetLabel(fp, (char *) NULL, (int) 0); -} +{ + return treeGetLabel(fp, (char *) NULL, (int) 0); +} -static int treeFindTipByLabelString(char *str, tree *tr) +static int treeFindTipByLabelString(char *str, tree *tr) { - int lookup = lookupWord(str, tr->nameHash); - - if(lookup > 0) - { - assert(! tr->nodep[lookup]->back); - return lookup; - } - else - { - printf("ERROR: Cannot find tree species: %s\n", str); - return 0; - } + int lookup = lookupWord(str, tr->nameHash); + + if(lookup > 0) + { + assert(!tr->nodep[lookup]->back); + return lookup; + } + else + { + printf("ERROR: Cannot find tree species: %s\n", str); + return 0; + } } static int treeFindTipName(FILE *fp, tree *tr) { - char str[nmlngth+2]; - int n; + char str[nmlngth+2]; + int n; - if(treeGetLabel(fp, str, nmlngth+2)) - n = treeFindTipByLabelString(str, tr); - else - n = 0; - + if(treeGetLabel(fp, str, nmlngth+2)) + n = treeFindTipByLabelString(str, tr); + else + n = 0; - return n; -} + + return n; +} static boolean treeProcessLength (FILE *fp, double *dptr) { - int ch; - - if ((ch = treeGetCh(fp)) == EOF) return FALSE; /* Skip comments */ - (void) ungetc(ch, fp); - - if (fscanf(fp, "%lf", dptr) != 1) { - printf("ERROR: treeProcessLength: Problem reading branch length\n"); - treeEchoContext(fp, stdout, 40); - printf("\n"); - return FALSE; - } - - return TRUE; + int ch; + + if ((ch = treeGetCh(fp)) == EOF) return FALSE; /* Skip comments */ + (void) ungetc(ch, fp); + + if (fscanf(fp, "%lf", dptr) != 1) { + printf("ERROR: treeProcessLength: Problem reading branch length\n"); + treeEchoContext(fp, stdout, 40); + printf("\n"); + return FALSE; + } + + return TRUE; } static int treeFlushLen (FILE *fp) { - double dummy; - int ch; - - ch = treeGetCh(fp); - - if (ch == ':') - { - ch = treeGetCh(fp); - - ungetc(ch, fp); - if(!treeProcessLength(fp, & dummy)) return 0; - return 1; - } - - - - if (ch != EOF) (void) ungetc(ch, fp); - return 1; -} + double dummy; + int ch; + + ch = treeGetCh(fp); + + if (ch == ':') + { + ch = treeGetCh(fp); + + ungetc(ch, fp); + if(!treeProcessLength(fp, &dummy)) return 0; + return 1; + } + + + + if (ch != EOF) (void) ungetc(ch, fp); + return 1; +} static void treeEchoContext (FILE *fp1, FILE *fp2, int n) { /* treeEchoContext */ - int ch; - boolean waswhite; - - waswhite = TRUE; - - while (n > 0 && ((ch = getc(fp1)) != EOF)) { - if (whitechar(ch)) { - ch = waswhite ? '\0' : ' '; - waswhite = TRUE; - } - else { - waswhite = FALSE; - } - - if (ch > '\0') {putc(ch, fp2); n--;} - } + int ch; + boolean waswhite; + + waswhite = TRUE; + + while (n > 0 && ((ch = getc(fp1)) != EOF)) { + if (whitechar(ch)) { + ch = waswhite ? '\0' : ' '; + waswhite = TRUE; + } + else { + waswhite = FALSE; + } + + if (ch > '\0') {putc(ch, fp2); n--;} + } } /* treeEchoContext */ static boolean treeNeedCh (FILE *fp, int c1, char *where) { - int c2; - - if ((c2 = treeGetCh(fp)) == c1) return TRUE; - - printf("ERROR: Expecting '%c' %s tree; found:", c1, where); - if (c2 == EOF) - { - printf("End-of-File"); - } - else - { - ungetc(c2, fp); - treeEchoContext(fp, stdout, 40); - } - putchar('\n'); - - if(c1 == ':') - printf("RAxML may be expecting to read a tree that contains branch lengths\n"); - - return FALSE; -} + int c2; + + if ((c2 = treeGetCh(fp)) == c1) return TRUE; + + printf("ERROR: Expecting '%c' %s tree; found:", c1, where); + if (c2 == EOF) + { + printf("End-of-File"); + } + else + { + ungetc(c2, fp); + treeEchoContext(fp, stdout, 40); + } + putchar('\n'); + + if(c1 == ':') + printf("RAxML may be expecting to read a tree that contains branch lengths\n"); + + return FALSE; +} static void addElementLen (FILE *fp, tree *tr, nodeptr p) -{ - nodeptr q; - int n, ch, fres; - - if ((ch = treeGetCh(fp)) == '(') - { - n = (tr->nextnode)++; - if (n > 2*(tr->mxtips) - 2) - { - if (n > 2*(tr->mxtips) - 1) - { - printf("ERROR: Too many internal nodes. Is tree rooted?\n"); - printf(" Deepest splitting should be a trifurcation.\n"); - assert(0); - } - else - assert(0); - } - - q = tr->nodep[n]; - - addElementLen(fp, tr, q->next); - if (! treeNeedCh(fp, ',', "in")) - assert(0); - addElementLen(fp, tr, q->next->next); - if (! treeNeedCh(fp, ')', "in")) - assert(0); - - (void) treeFlushLabel(fp); - } - else - { - ungetc(ch, fp); - if ((n = treeFindTipName(fp, tr)) <= 0) - assert(0); - q = tr->nodep[n]; - tr->nodesInTree[n] = 1; - - if (tr->start->number > n) - tr->start = q; - (tr->ntips)++; - } - - - fres = treeFlushLen(fp); - if(!fres) - assert(0); - - - p->back = q; - q->back = p; -} +{ + nodeptr q; + int n, ch, fres; + + if ((ch = treeGetCh(fp)) == '(') + { + n = (tr->nextnode)++; + if (n > 2*(tr->mxtips) - 2) + { + if (n > 2*(tr->mxtips) - 1) + { + printf("ERROR: Too many internal nodes. Is tree rooted?\n"); + printf(" Deepest splitting should be a trifurcation.\n"); + assert(0); + } + else + assert(0); + } + + q = tr->nodep[n]; + + addElementLen(fp, tr, q->next); + if (!treeNeedCh(fp, ',', "in")) + assert(0); + addElementLen(fp, tr, q->next->next); + if (!treeNeedCh(fp, ')', "in")) + assert(0); + + (void) treeFlushLabel(fp); + } + else + { + ungetc(ch, fp); + if ((n = treeFindTipName(fp, tr)) <= 0) + assert(0); + q = tr->nodep[n]; + tr->nodesInTree[n] = 1; + + if (tr->start->number > n) + tr->start = q; + (tr->ntips)++; + } + + + fres = treeFlushLen(fp); + if(!fres) + assert(0); + + + p->back = q; + q->back = p; +} @@ -458,61 +458,61 @@ static void addElementLen (FILE *fp, tree *tr, nodeptr p) void treeReadLen (FILE *fp, tree *tr) { - nodeptr - p; - - int - i, - ch; - - for (i = 1; i <= tr->mxtips; i++) - tr->nodep[i]->back = (node *) NULL; - - for(i = tr->mxtips + 1; i < 2 * tr->mxtips; i++) - { - tr->nodep[i]->back = (nodeptr)NULL; - tr->nodep[i]->next->back = (nodeptr)NULL; - tr->nodep[i]->next->next->back = (nodeptr)NULL; - tr->nodep[i]->number = i; - tr->nodep[i]->next->number = i; - tr->nodep[i]->next->next->number = i; - } - - tr->start = tr->nodep[tr->mxtips + 1]; - - tr->ntips = 0; - tr->nextnode = tr->mxtips + 1; - - p = tr->nodep[(tr->nextnode)++]; - - while((ch = treeGetCh(fp)) != '('); - - addElementLen(fp, tr, p); - - if (!treeNeedCh(fp, ',', "in")) - assert(0); - - addElementLen(fp, tr, p->next); - - - if ((ch = treeGetCh(fp)) == ',') - addElementLen(fp, tr, p->next->next); - else - assert(0); - - if (! treeNeedCh(fp, ')', "in")) - assert(0); - - - treeFlushLabel(fp); - - if (! treeFlushLen(fp)) - assert(0); - - if (! treeNeedCh(fp, ';', "at end of")) - assert(0); - - + nodeptr + p; + + int + i, + ch; + + for (i = 1; i <= tr->mxtips; i++) + tr->nodep[i]->back = (node *) NULL; + + for(i = tr->mxtips + 1; i < 2 * tr->mxtips; i++) + { + tr->nodep[i]->back = (nodeptr)NULL; + tr->nodep[i]->next->back = (nodeptr)NULL; + tr->nodep[i]->next->next->back = (nodeptr)NULL; + tr->nodep[i]->number = i; + tr->nodep[i]->next->number = i; + tr->nodep[i]->next->next->number = i; + } + + tr->start = tr->nodep[tr->mxtips + 1]; + + tr->ntips = 0; + tr->nextnode = tr->mxtips + 1; + + p = tr->nodep[(tr->nextnode)++]; + + while((ch = treeGetCh(fp)) != '('); + + addElementLen(fp, tr, p); + + if (!treeNeedCh(fp, ',', "in")) + assert(0); + + addElementLen(fp, tr, p->next); + + + if ((ch = treeGetCh(fp)) == ',') + addElementLen(fp, tr, p->next->next); + else + assert(0); + + if (!treeNeedCh(fp, ')', "in")) + assert(0); + + + treeFlushLabel(fp); + + if (!treeFlushLen(fp)) + assert(0); + + if (!treeNeedCh(fp, ';', "at end of")) + assert(0); + + } @@ -526,33 +526,33 @@ void treeReadLen (FILE *fp, tree *tr) extern FILE *myfopen(const char *path, const char *mode) { - FILE *fp = fopen(path, mode); - - if(strcmp(mode,"r") == 0 || strcmp(mode,"rb") == 0) - { - if(fp) - return fp; - else - { - - printf("The file %s you want to open for reading does not exist, exiting ...\n", path); - exit(-1); - return (FILE *)NULL; - } - } - else - { - if(fp) - return fp; - else - { - - printf("The file %s RAxML wants to open for writing or appending can not be opened [mode: %s], exiting ...\n", - path, mode); - exit(-1); - return (FILE *)NULL; - } - } + FILE *fp = fopen(path, mode); + + if(strcmp(mode,"r") == 0 || strcmp(mode,"rb") == 0) + { + if(fp) + return fp; + else + { + + printf("The file %s you want to open for reading does not exist, exiting ...\n", path); + exit(-1); + return (FILE *)NULL; + } + } + else + { + if(fp) + return fp; + else + { + + printf("The file %s RAxML wants to open for writing or appending can not be opened [mode: %s], exiting ...\n", + path, mode); + exit(-1); + return (FILE *)NULL; + } + } } @@ -560,18 +560,18 @@ extern FILE *myfopen(const char *path, const char *mode) void printBothOpen(const char* format, ... ) { - FILE *f = myfopen(infoFileName, "ab"); + FILE *f = myfopen(infoFileName, "ab"); - va_list args; - va_start(args, format); - vfprintf(f, format, args ); - va_end(args); + va_list args; + va_start(args, format); + vfprintf(f, format, args ); + va_end(args); - va_start(args, format); - vprintf(format, args ); - va_end(args); + va_start(args, format); + vprintf(format, args ); + va_end(args); - fclose(f); + fclose(f); } @@ -591,15 +591,15 @@ void printBothOpen(const char* format, ... ) double gettime(void) { #ifdef WIN32 - time_t tp; - struct tm localtm; - tp = time(NULL); - localtm = *localtime(&tp); - return 60.0*localtm.tm_min + localtm.tm_sec; + time_t tp; + struct tm localtm; + tp = time(NULL); + localtm = *localtime(&tp); + return 60.0*localtm.tm_min + localtm.tm_sec; #else - struct timeval ttime; - gettimeofday(&ttime , NULL); - return ttime.tv_sec + ttime.tv_usec * 0.000001; + struct timeval ttime; + gettimeofday(&ttime, NULL); + return ttime.tv_sec + ttime.tv_usec * 0.000001; #endif } @@ -608,19 +608,19 @@ double gettime(void) static int filexists(char *filename) { - FILE *fp; - int res; - fp = fopen(filename,"rb"); - - if(fp) - { - res = 1; - fclose(fp); - } - else - res = 0; - - return res; + FILE *fp; + int res; + fp = fopen(filename,"rb"); + + if(fp) + { + res = 1; + fclose(fp); + } + else + res = 0; + + return res; } @@ -634,27 +634,27 @@ static int filexists(char *filename) static void getnums (rawdata *rdta) { - if (fscanf(INFILE, "%d %d", & rdta->numsp, & rdta->sites) != 2) - { - - printf("ERROR: Problem reading number of species and sites\n"); - exit(-1); - } - - if (rdta->numsp < 4) - { - - printf("TOO FEW SPECIES\n"); - exit(-1); - } - - if (rdta->sites < 1) - { - printf("TOO FEW SITES\n"); - exit(-1); - } - - return; + if (fscanf(INFILE, "%d %d", &rdta->numsp, &rdta->sites) != 2) + { + + printf("ERROR: Problem reading number of species and sites\n"); + exit(-1); + } + + if (rdta->numsp < 4) + { + + printf("TOO FEW SPECIES\n"); + exit(-1); + } + + if (rdta->sites < 1) + { + printf("TOO FEW SITES\n"); + exit(-1); + } + + return; } @@ -663,18 +663,18 @@ static void getnums (rawdata *rdta) static boolean whitechar (int ch) { - return (ch == ' ' || ch == '\n' || ch == '\t' || ch == '\r'); + return (ch == ' ' || ch == '\n' || ch == '\t' || ch == '\r'); } static void uppercase (int *chptr) { - int ch; + int ch; - ch = *chptr; - if ((ch >= 'a' && ch <= 'i') || (ch >= 'j' && ch <= 'r') - || (ch >= 's' && ch <= 'z')) - *chptr = ch + 'A' - 'a'; + ch = *chptr; + if ((ch >= 'a' && ch <= 'i') || (ch >= 'j' && ch <= 'r') + || (ch >= 's' && ch <= 'z')) + *chptr = ch + 'A' - 'a'; } @@ -682,97 +682,97 @@ static void uppercase (int *chptr) // Comments are by Mark Jones trying to figure out what is happening static void getyspace (rawdata *rdta) { - size_t size = 4 * ((size_t)(rdta->sites / 4 + 1)); // memory size of a 'row' in the 2d array . somehow 4 things are enough to represent 4 elements I guess? - int i; - unsigned char *y0; + size_t size = 4 * ((size_t)(rdta->sites / 4 + 1)); // memory size of a 'row' in the 2d array . somehow 4 things are enough to represent 4 elements I guess? + int i; + unsigned char *y0; - rdta->y = (unsigned char **) malloc((size_t)(rdta->numsp + 1) * sizeof(unsigned char *)); // create space for the top-level array - assert(rdta->y); //we found space right? + rdta->y = (unsigned char **) malloc((size_t)(rdta->numsp + 1) * sizeof(unsigned char *)); // create space for the top-level array + assert(rdta->y); //we found space right? - y0 = (unsigned char *) malloc(((size_t)(rdta->numsp + 1)) * size * sizeof(unsigned char)); //this a pointer address, after which is the requested amount of free space... - assert(y0); + y0 = (unsigned char *) malloc(((size_t)(rdta->numsp + 1)) * size * sizeof(unsigned char)); //this a pointer address, after which is the requested amount of free space... + assert(y0); - // the location of y and its array may be different from the location of the y_0 and the associated actual 2d array. - // each element of y points to the corresponding point in the 2d array space corresponding to its row. - // this isn't going to make sense to anyone but me I guess but that's ok, I'm going to delete this all before pushing. - // SO... what I'm going to want to do is something like this, but for 3 dimensions??? + // the location of y and its array may be different from the location of the y_0 and the associated actual 2d array. + // each element of y points to the corresponding point in the 2d array space corresponding to its row. + // this isn't going to make sense to anyone but me I guess but that's ok, I'm going to delete this all before pushing. + // SO... what I'm going to want to do is something like this, but for 3 dimensions??? - rdta->y0 = y0; // pointer of first element of 2d array + rdta->y0 = y0; // pointer of first element of 2d array - for (i = 0; i <= rdta->numsp; i++) - { - rdta->y[i] = y0; // row i of rdta is assigned to pointer (original y_0) + i*(row size) - y0 += size; - } + for (i = 0; i <= rdta->numsp; i++) + { + rdta->y[i] = y0; // row i of rdta is assigned to pointer (original y_0) + i*(row size) + y0 += size; + } - return; + return; } // Returns a 2d int array of dimensions [d1][d2], and allocates memory for it static int **get2dIntArray(int d1, int d2) - { - //printf("creating 2d int array... "); - int **arr; - arr = (int **)malloc(sizeof(int *) * d1); - arr[0] = (int *)malloc(sizeof(int) * d2 * d1); - int i; - for(i = 0; i < d1; i++) - arr[i] = (*arr + d2 * i); - //printf("done \n"); - return arr; - - } +{ + //printf("creating 2d int array... "); + int **arr; + arr = (int **)malloc(sizeof(int *) * d1); + arr[0] = (int *)malloc(sizeof(int) * d2 * d1); + int i; + for(i = 0; i < d1; i++) + arr[i] = (*arr + d2 * i); + //printf("done \n"); + return arr; + +} // Returns a 3d int array of dimensions [d1][d2][d3], and allocates memory for it static int ***get3dIntArray(int d1, int d2, int d3) - { - //printf("creating 3d int array... "); - int b,i; - int ***arr; - arr = (int ***)malloc(sizeof(int **) * d1); - arr[0] = (int **)malloc(sizeof(int *) * d2 * d1); - arr[0][0] = (int *)malloc(sizeof(int) * d3 * d2 * d1); - - for(b = 0; b < d1; b++) - { - arr[b] = (*arr + d2 * b); - for(i = 0; i < d2; i++) - { - arr[b][i] = (**arr + d3 * i); - } - } - //printf("done \n"); - return arr; - } +{ + //printf("creating 3d int array... "); + int b,i; + int ***arr; + arr = (int ***)malloc(sizeof(int **) * d1); + arr[0] = (int **)malloc(sizeof(int *) * d2 * d1); + arr[0][0] = (int *)malloc(sizeof(int) * d3 * d2 * d1); + + for(b = 0; b < d1; b++) + { + arr[b] = (*arr + d2 * b); + for(i = 0; i < d2; i++) + { + arr[b][i] = (**arr + d3 * i); + } + } + //printf("done \n"); + return arr; +} // Returns a 3d int array of dimensions [d1][d2][d3], and allocates memory for it static float ***get3dFloatArray(int d1, int d2, int d3) - { - //printf("creating 3d float array... "); - int b,i; - float ***arr; - arr = (float ***)malloc(sizeof(float **) * d1); - arr[0] = (float **)malloc(sizeof(float *) * d2 * d1); - arr[0][0] = (float *)malloc(sizeof(float) * d3 * d2 * d1); - - for(b = 0; b < d1; b++) - { - arr[b] = (*arr + d2 * b); - for(i = 0; i < d2; i++) - { - arr[b][i] = (**arr + d3 * i); - for(int j = 0; j < d3; j++){ - arr[b][i][j]=FLT_MAX; - } - - } - } - //printf("done \n"); - return arr; - } +{ + //printf("creating 3d float array... "); + int b,i; + float ***arr; + arr = (float ***)malloc(sizeof(float **) * d1); + arr[0] = (float **)malloc(sizeof(float *) * d2 * d1); + arr[0][0] = (float *)malloc(sizeof(float) * d3 * d2 * d1); + + for(b = 0; b < d1; b++) + { + arr[b] = (*arr + d2 * b); + for(i = 0; i < d2; i++) + { + arr[b][i] = (**arr + d3 * i); + for(int j = 0; j < d3; j++) { + arr[b][i][j]=FLT_MAX; + } + + } + } + //printf("done \n"); + return arr; +} @@ -780,293 +780,293 @@ static float ***get3dFloatArray(int d1, int d2, int d3) static boolean setupTree (tree *tr) { - nodeptr p0, p, q; - int - i, - j, - tips, - inter; + nodeptr p0, p, q; + int + i, + j, + tips, + inter; - - tips = tr->mxtips; - inter = tr->mxtips - 1; - - tr->yVector = (unsigned char **) malloc((size_t)(tr->mxtips + 1) * sizeof(unsigned char *)); + tips = tr->mxtips; + inter = tr->mxtips - 1; - - tr->treeStringLength = tr->mxtips * (nmlngth+128) + 256 + tr->mxtips * 2; + tr->yVector = (unsigned char **) malloc((size_t)(tr->mxtips + 1) * sizeof(unsigned char *)); - tr->tree_string = (char*)calloc(tr->treeStringLength, sizeof(char)); - /*TODO, must that be so long ?*/ - + tr->treeStringLength = tr->mxtips * (nmlngth+128) + 256 + tr->mxtips * 2; - tr->nameList = (char **)malloc(sizeof(char *) * (size_t)(tips + 1)); + tr->tree_string = (char*)calloc(tr->treeStringLength, sizeof(char)); - if (!(p0 = (nodeptr) malloc((size_t)(tips + 3*inter) * sizeof(node)))) - { - printf("ERROR: Unable to obtain sufficient tree memory\n"); - return FALSE; - } + /*TODO, must that be so long ?*/ - if (!(tr->nodep = (nodeptr *) malloc((size_t)(2*tr->mxtips) * sizeof(nodeptr)))) - { - printf("ERROR: Unable to obtain sufficient tree memory, too\n"); - return FALSE; - } - tr->nodep[0] = (node *) NULL; /* Use as 1-based array */ - for (i = 1; i <= tips; i++) - { - p = p0++; + tr->nameList = (char **)malloc(sizeof(char *) * (size_t)(tips + 1)); - - p->x = 0; - p->number = i; - p->next = p; - p->back = (node *)NULL; - + if (!(p0 = (nodeptr) malloc((size_t)(tips + 3*inter) * sizeof(node)))) + { + printf("ERROR: Unable to obtain sufficient tree memory\n"); + return FALSE; + } - tr->nodep[i] = p; - } + if (!(tr->nodep = (nodeptr *) malloc((size_t)(2*tr->mxtips) * sizeof(nodeptr)))) + { + printf("ERROR: Unable to obtain sufficient tree memory, too\n"); + return FALSE; + } - for (i = tips + 1; i <= tips + inter; i++) - { - q = (node *) NULL; - for (j = 1; j <= 3; j++) - { - p = p0++; - if(j == 1) - p->x = 1; - else - p->x = 0; - p->number = i; - p->next = q; - - p->back = (node *) NULL; - + tr->nodep[0] = (node *) NULL; /* Use as 1-based array */ + + for (i = 1; i <= tips; i++) + { + p = p0++; - q = p; - } - p->next->next->next = p; - tr->nodep[i] = p; - } + p->x = 0; + p->number = i; + p->next = p; + p->back = (node *)NULL; - - tr->start = (node *) NULL; - + tr->nodep[i] = p; + } - tr->ntips = 0; - tr->nextnode = 0; + for (i = tips + 1; i <= tips + inter; i++) + { + q = (node *) NULL; + for (j = 1; j <= 3; j++) + { + p = p0++; + if(j == 1) + p->x = 1; + else + p->x = 0; + p->number = i; + p->next = q; - + p->back = (node *) NULL; - return TRUE; + + + q = p; + } + p->next->next->next = p; + tr->nodep[i] = p; + } + + + tr->start = (node *) NULL; + + + + tr->ntips = 0; + tr->nextnode = 0; + + + + return TRUE; } static void checkTaxonName(char *buffer, int len) { - int i; - - for(i = 0; i < len - 1; i++) - { - boolean valid; - - switch(buffer[i]) - { - case '\0': - case '\t': - case '\n': - case '\r': - case ' ': - case ':': - case ',': - case '(': - case ')': - case ';': - case '[': - case ']': - case '\'': - valid = FALSE; - break; - default: - valid = TRUE; - } - - if(!valid) - { - printf("ERROR: Taxon Name \"%s\" is invalid at position %d, it contains illegal character %c\n", buffer, i, buffer[i]); - printf("Illegal characters in taxon-names are: tabulators, carriage returns, spaces, \":\", \",\", \")\", \"(\", \";\", \"]\", \"[\", \"\'\"\n"); - printf("Exiting\n"); - exit(-1); - } - - } - assert(buffer[len - 1] == '\0'); + int i; + + for(i = 0; i < len - 1; i++) + { + boolean valid; + + switch(buffer[i]) + { + case '\0': + case '\t': + case '\n': + case '\r': + case ' ': + case ':': + case ',': + case '(': + case ')': + case ';': + case '[': + case ']': + case '\'': + valid = FALSE; + break; + default: + valid = TRUE; + } + + if(!valid) + { + printf("ERROR: Taxon Name \"%s\" is invalid at position %d, it contains illegal character %c\n", buffer, i, buffer[i]); + printf("Illegal characters in taxon-names are: tabulators, carriage returns, spaces, \":\", \",\", \")\", \"(\", \";\", \"]\", \"[\", \"\'\"\n"); + printf("Exiting\n"); + exit(-1); + } + + } + assert(buffer[len - 1] == '\0'); } static boolean getdata(rawdata *rdta, tree *tr) { - size_t - i, - j, - basesread, - basesnew; - - int - ch, my_i, meaning, - len, - meaningDNA[256]; - - boolean - allread, - firstpass; - - char - buffer[nmlngth + 2]; - - for(i = 0; i < 256; i++) - meaningDNA[i] = -1; - - meaningDNA['A'] = 1; - meaningDNA['B'] = 14; - meaningDNA['C'] = 2; - meaningDNA['D'] = 13; - meaningDNA['G'] = 4; - meaningDNA['H'] = 11; - meaningDNA['K'] = 12; - meaningDNA['M'] = 3; - meaningDNA['R'] = 5; - meaningDNA['S'] = 6; - meaningDNA['T'] = 8; - meaningDNA['U'] = 8; - meaningDNA['V'] = 7; - meaningDNA['W'] = 9; - meaningDNA['Y'] = 10; - - meaningDNA['N'] = - meaningDNA['O'] = - meaningDNA['X'] = - meaningDNA['-'] = - meaningDNA['?'] = - 15; - - - - /*******************************************************************/ - - basesread = basesnew = 0; - - allread = FALSE; - firstpass = TRUE; - ch = ' '; - - while (! allread) - { - for (i = 1; i <= (size_t)tr->mxtips; i++) - { - if (firstpass) - { - ch = getc(INFILE); - while(ch == ' ' || ch == '\n' || ch == '\t' || ch == '\r') - ch = getc(INFILE); - - my_i = 0; - - do - { - buffer[my_i] = ch; - ch = getc(INFILE); - my_i++; - if(my_i >= nmlngth) - { - printf("Taxon Name to long at taxon %zd, adapt constant nmlngth in\n", i); - printf("axml.h, current setting %d\n", nmlngth); - exit(0); - } - } - while(ch != ' ' && ch != '\n' && ch != '\t' && ch != '\r'); - - while(ch == ' ' || ch == '\n' || ch == '\t' || ch == '\r') - ch = getc(INFILE); - - ungetc(ch, INFILE); - - buffer[my_i] = '\0'; - len = strlen(buffer) + 1; - checkTaxonName(buffer, len); - tr->nameList[i] = (char *)malloc(sizeof(char) * (size_t)len); - strcpy(tr->nameList[i], buffer); - } - - j = basesread; - - while ((j < (size_t)rdta->sites) && ((ch = getc(INFILE)) != EOF) && (ch != '\n') && (ch != '\r')) - { - uppercase(& ch); - meaning = meaningDNA[ch]; - - if (meaning != -1) - { - j++; - rdta->y[i][j - 1] = meaning; - } - else - { - if(!whitechar(ch)) - { - printf("ERROR: Bad base (%c) at site %zd of sequence %zd\n", - ch, j + 1, i); - return FALSE; - } - } - } - - - - if (ch == EOF) - { - printf("ERROR: End-of-file at site %zd of sequence %zd\n", j + 1, i); - return FALSE; - } - - if (! firstpass && (j == basesread)) - i--; - else - { - if (i == 1) - basesnew = j; - else - if (j != basesnew) - { - printf("ERROR: Sequences out of alignment\n"); - printf("%zd (instead of %zd) residues read in sequence %zd %s\n", - j - basesread, basesnew - basesread, i, tr->nameList[i]); - return FALSE; - } - } - while (ch != '\n' && ch != EOF && ch != '\r') ch = getc(INFILE); /* flush line *//* PC-LINEBREAK*/ - } - - firstpass = FALSE; - basesread = basesnew; - allread = (basesread >= (size_t)rdta->sites); - } - - - - - return TRUE; + size_t + i, + j, + basesread, + basesnew; + + int + ch, my_i, meaning, + len, + meaningDNA[256]; + + boolean + allread, + firstpass; + + char + buffer[nmlngth + 2]; + + for(i = 0; i < 256; i++) + meaningDNA[i] = -1; + + meaningDNA['A'] = 1; + meaningDNA['B'] = 14; + meaningDNA['C'] = 2; + meaningDNA['D'] = 13; + meaningDNA['G'] = 4; + meaningDNA['H'] = 11; + meaningDNA['K'] = 12; + meaningDNA['M'] = 3; + meaningDNA['R'] = 5; + meaningDNA['S'] = 6; + meaningDNA['T'] = 8; + meaningDNA['U'] = 8; + meaningDNA['V'] = 7; + meaningDNA['W'] = 9; + meaningDNA['Y'] = 10; + + meaningDNA['N'] = + meaningDNA['O'] = + meaningDNA['X'] = + meaningDNA['-'] = + meaningDNA['?'] = + 15; + + + + /*******************************************************************/ + + basesread = basesnew = 0; + + allread = FALSE; + firstpass = TRUE; + ch = ' '; + + while (!allread) + { + for (i = 1; i <= (size_t)tr->mxtips; i++) + { + if (firstpass) + { + ch = getc(INFILE); + while(ch == ' ' || ch == '\n' || ch == '\t' || ch == '\r') + ch = getc(INFILE); + + my_i = 0; + + do + { + buffer[my_i] = ch; + ch = getc(INFILE); + my_i++; + if(my_i >= nmlngth) + { + printf("Taxon Name to long at taxon %zd, adapt constant nmlngth in\n", i); + printf("axml.h, current setting %d\n", nmlngth); + exit(0); + } + } + while(ch != ' ' && ch != '\n' && ch != '\t' && ch != '\r'); + + while(ch == ' ' || ch == '\n' || ch == '\t' || ch == '\r') + ch = getc(INFILE); + + ungetc(ch, INFILE); + + buffer[my_i] = '\0'; + len = strlen(buffer) + 1; + checkTaxonName(buffer, len); + tr->nameList[i] = (char *)malloc(sizeof(char) * (size_t)len); + strcpy(tr->nameList[i], buffer); + } + + j = basesread; + + while ((j < (size_t)rdta->sites) && ((ch = getc(INFILE)) != EOF) && (ch != '\n') && (ch != '\r')) + { + uppercase(&ch); + meaning = meaningDNA[ch]; + + if (meaning != -1) + { + j++; + rdta->y[i][j - 1] = meaning; + } + else + { + if(!whitechar(ch)) + { + printf("ERROR: Bad base (%c) at site %zd of sequence %zd\n", + ch, j + 1, i); + return FALSE; + } + } + } + + + + if (ch == EOF) + { + printf("ERROR: End-of-file at site %zd of sequence %zd\n", j + 1, i); + return FALSE; + } + + if (!firstpass && (j == basesread)) + i--; + else + { + if (i == 1) + basesnew = j; + else + if (j != basesnew) + { + printf("ERROR: Sequences out of alignment\n"); + printf("%zd (instead of %zd) residues read in sequence %zd %s\n", + j - basesread, basesnew - basesread, i, tr->nameList[i]); + return FALSE; + } + } + while (ch != '\n' && ch != EOF && ch != '\r') ch = getc(INFILE); /* flush line *//* PC-LINEBREAK*/ + } + + firstpass = FALSE; + basesread = basesnew; + allread = (basesread >= (size_t)rdta->sites); + } + + + + + return TRUE; } @@ -1081,38 +1081,38 @@ static boolean getdata(rawdata *rdta, tree *tr) static void getinput(analdef *adef, rawdata *rdta, tree *tr) //, cruncheddata *cdta //not needed anymore { - int i; - - - INFILE = myfopen(seq_file, "rb"); - - getnums(rdta); - - tr->mxtips = rdta->numsp; - - - - getyspace(rdta); - - - setupTree(tr); - - - if(!getdata(rdta, tr)) - { - printf("Problem reading alignment file \n"); - exit(1); - } - // cdta->endsite = rdta->sites; // not needed anymore - - if(adef->restart) - { - tr->nameHash = initStringHashTable(10 * tr->mxtips); - for(i = 1; i <= tr->mxtips; i++) - addword(tr->nameList[i], tr->nameHash, i); - } - - fclose(INFILE); + int i; + + + INFILE = myfopen(seq_file, "rb"); + + getnums(rdta); + + tr->mxtips = rdta->numsp; + + + + getyspace(rdta); + + + setupTree(tr); + + + if(!getdata(rdta, tr)) + { + printf("Problem reading alignment file \n"); + exit(1); + } + // cdta->endsite = rdta->sites; // not needed anymore + + if(adef->restart) + { + tr->nameHash = initStringHashTable(10 * tr->mxtips); + for(i = 1; i <= tr->mxtips; i++) + addword(tr->nameList[i], tr->nameHash, i); + } + + fclose(INFILE); } @@ -1134,16 +1134,16 @@ static void getinput(analdef *adef, rawdata *rdta, tree *tr) //, cruncheddata *c static boolean makevalues(rawdata *rdta, tree *tr) //, cruncheddata *cdta not needed anymore { - int - i; + int + i; + + tr->rdta = rdta; + //tr->cdta = cdta; //not needed anymore - tr->rdta = rdta; - //tr->cdta = cdta; //not needed anymore - - for(i = 0; i <= rdta->numsp; i++) - tr->yVector[i] = rdta->y[i]; + for(i = 0; i <= rdta->numsp; i++) + tr->yVector[i] = rdta->y[i]; - return TRUE; + return TRUE; } @@ -1160,10 +1160,10 @@ static boolean makevalues(rawdata *rdta, tree *tr) //, cruncheddata *cdta not ne static void initAdef(analdef *adef) -{ - adef->restart = FALSE; - adef->parsimonySeed = 0; - adef->numberOfTrees = 1; +{ + adef->restart = FALSE; + adef->parsimonySeed = 0; + adef->numberOfTrees = 1; } @@ -1173,70 +1173,70 @@ static void initAdef(analdef *adef) static int mygetopt(int argc, char **argv, char *opts, int *optind, char **optarg) { - static int sp = 1; - register int c; - register char *cp; - - if(sp == 1) - { - if(*optind >= argc || argv[*optind][0] != '-' || argv[*optind][1] == '\0') - return -1; - } - else - { - if(strcmp(argv[*optind], "--") == 0) - { - *optind = *optind + 1; - return -1; - } - } - - c = argv[*optind][sp]; - if(c == ':' || (cp=strchr(opts, c)) == 0) - { - printf(": illegal option -- %c \n", c); - if(argv[*optind][++sp] == '\0') - { - *optind = *optind + 1; - sp = 1; - } - return('?'); - } - if(*++cp == ':') - { - if(argv[*optind][sp+1] != '\0') - { - *optarg = &argv[*optind][sp+1]; - *optind = *optind + 1; - } - else - { - *optind = *optind + 1; - if(*optind >= argc) - { - printf(": option requires an argument -- %c\n", c); - sp = 1; - return('?'); - } - else - { - *optarg = argv[*optind]; - *optind = *optind + 1; - } - } - sp = 1; - } - else - { - if(argv[*optind][++sp] == '\0') - { - sp = 1; - *optind = *optind + 1; - } - *optarg = 0; - } - return(c); - } + static int sp = 1; + register int c; + register char *cp; + + if(sp == 1) + { + if(*optind >= argc || argv[*optind][0] != '-' || argv[*optind][1] == '\0') + return -1; + } + else + { + if(strcmp(argv[*optind], "--") == 0) + { + *optind = *optind + 1; + return -1; + } + } + + c = argv[*optind][sp]; + if(c == ':' || (cp=strchr(opts, c)) == 0) + { + printf(": illegal option -- %c \n", c); + if(argv[*optind][++sp] == '\0') + { + *optind = *optind + 1; + sp = 1; + } + return('?'); + } + if(*++cp == ':') + { + if(argv[*optind][sp+1] != '\0') + { + *optarg = &argv[*optind][sp+1]; + *optind = *optind + 1; + } + else + { + *optind = *optind + 1; + if(*optind >= argc) + { + printf(": option requires an argument -- %c\n", c); + sp = 1; + return('?'); + } + else + { + *optarg = argv[*optind]; + *optind = *optind + 1; + } + } + sp = 1; + } + else + { + if(argv[*optind][++sp] == '\0') + { + sp = 1; + *optind = *optind + 1; + } + *optarg = 0; + } + return(c); +} @@ -1244,92 +1244,92 @@ static int mygetopt(int argc, char **argv, char *opts, int *optind, char **optar static void analyzeRunId(char id[128]) { - int i = 0; - - while(id[i] != '\0') - { - if(i >= 128) - { - printf("Error: run id after \"-n\" is too long, it has %d characters please use a shorter one\n", i); - assert(0); - } - - if(id[i] == '/') - { - printf("Error character %c not allowed in run ID\n", id[i]); - assert(0); - } - - - i++; - } - - if(i == 0) - { - printf("Error: please provide a string for the run id after \"-n\" \n"); - assert(0); - } + int i = 0; + + while(id[i] != '\0') + { + if(i >= 128) + { + printf("Error: run id after \"-n\" is too long, it has %d characters please use a shorter one\n", i); + assert(0); + } + + if(id[i] == '/') + { + printf("Error character %c not allowed in run ID\n", id[i]); + assert(0); + } + + + i++; + } + + if(i == 0) + { + printf("Error: please provide a string for the run id after \"-n\" \n"); + assert(0); + } } static void get_args(int argc, char *argv[], analdef *adef) { - - int - optind = 1, - c; - - boolean - bad_opt =FALSE; - - char *optarg; - - run_id[0] = 0; - seq_file[0] = 0; - - adef->parsimonySeed=1; - adef->numberOfBlocks=2; - - while(!bad_opt && ((c = mygetopt(argc,argv,"p:n:s:t:N:b:", &optind, &optarg))!=-1)) - { - switch(c) - { - case 'N': - sscanf(optarg,"%d", &(adef->numberOfTrees)); - if(adef->numberOfTrees <= 0) - { - printf("number of trees can't be smaller than 1\n"); - exit(-1); - } - break; - case 'b': - sscanf(optarg,"%d", &(adef->numberOfBlocks)); - break; - case 't': - strcpy(tree_file, optarg); - adef->restart = TRUE; - break; - case 's': - strcpy(seq_file, optarg); - break; - case 'p': - sscanf(optarg,"%ld", &(adef->parsimonySeed)); - if(adef->parsimonySeed <= 0) - { - printf("Parsimony seed specified via -p must be greater than zero\n"); - exit(-1); - } - break; - case 'n': - strcpy(run_id, optarg); - analyzeRunId(run_id); - break; - default: - assert(0); - } - } - - return; + + int + optind = 1, + c; + + boolean + bad_opt =FALSE; + + char *optarg; + + run_id[0] = 0; + seq_file[0] = 0; + + adef->parsimonySeed=1; + adef->numberOfBlocks=2; + + while(!bad_opt && ((c = mygetopt(argc,argv,"p:n:s:t:N:b:", &optind, &optarg))!=-1)) + { + switch(c) + { + case 'N': + sscanf(optarg,"%d", &(adef->numberOfTrees)); + if(adef->numberOfTrees <= 0) + { + printf("number of trees can't be smaller than 1\n"); + exit(-1); + } + break; + case 'b': + sscanf(optarg,"%d", &(adef->numberOfBlocks)); + break; + case 't': + strcpy(tree_file, optarg); + adef->restart = TRUE; + break; + case 's': + strcpy(seq_file, optarg); + break; + case 'p': + sscanf(optarg,"%ld", &(adef->parsimonySeed)); + if(adef->parsimonySeed <= 0) + { + printf("Parsimony seed specified via -p must be greater than zero\n"); + exit(-1); + } + break; + case 'n': + strcpy(run_id, optarg); + analyzeRunId(run_id); + break; + default: + assert(0); + } + } + + return; } @@ -1339,25 +1339,25 @@ static void get_args(int argc, char *argv[], analdef *adef) static void makeFileNames(void) -{ - strcpy(resultFileName, "RAxML_parsimonyTree."); - strcpy(infoFileName, "RAxML_info."); - - strcat(resultFileName, run_id); - strcat(infoFileName, run_id); - - if(filexists(infoFileName)) - { - printf("RAxML output files with the run ID <%s> already exist \n", run_id); - - exit(-1); - } +{ + strcpy(resultFileName, "RAxML_parsimonyTree."); + strcpy(infoFileName, "RAxML_info."); + + strcat(resultFileName, run_id); + strcat(infoFileName, run_id); + + if(filexists(infoFileName)) + { + printf("RAxML output files with the run ID <%s> already exist \n", run_id); + + exit(-1); + } } - + @@ -1375,48 +1375,48 @@ static void makeFileNames(void) static int iterated_bitcount(unsigned int n) { - int - count=0; - - while(n) - { - count += n & 0x1u ; - n >>= 1 ; - } - - return count; + int + count=0; + + while(n) + { + count += n & 0x1u; + n >>= 1; + } + + return count; } static char bits_in_16bits [0x1u << 16]; static void compute_bits_in_16bits(void) { - unsigned int i; - - assert(sizeof(unsigned int) == 4); - - for (i = 0; i < (0x1u<<16); i++) - bits_in_16bits[i] = iterated_bitcount(i); - - return ; + unsigned int i; + + assert(sizeof(unsigned int) == 4); + + for (i = 0; i < (0x1u<<16); i++) + bits_in_16bits[i] = iterated_bitcount(i); + + return; } unsigned int precomputed16_bitcount (unsigned int n) { - /* works only for 32-bit int*/ - - return bits_in_16bits [n & 0xffffu] - + bits_in_16bits [(n >> 16) & 0xffffu] ; + /* works only for 32-bit int*/ + + return bits_in_16bits [n & 0xffffu] + + bits_in_16bits [(n >> 16) & 0xffffu]; } -static float getHomoplasyRatio(int i, int j, int homoplasyValue) +static float getHomoplasyRatio(int i, int j, int homoplasyValue) { - assert (i <= j); // program breaks if called with i greater than j - return ((float) homoplasyValue) / ((float)(j-i+1)); + assert (i <= j);// program breaks if called with i greater than j + return ((float) homoplasyValue) / ((float)(j-i+1)); } @@ -1428,277 +1428,277 @@ static float getHomoplasyRatio(int i, int j, int homoplasyValue) static int getOptHomoplasyOfContiguousBlock(tree *tr, analdef *adef, int startCharacter, int endCharacter) { - assert (startCharacter <= endCharacter); - - int parsimonyScore = makeParsimonyTreeFastDNA(tr, adef, startCharacter, endCharacter); - // TODO: makeParsimonyTreeFastDNA changes the tree, right? - // Which means that after the first time we call makeParsimonyTreeFastDNA, subsequent calls won't be starting with a 'clean' tree. - // This might be useful later on, but to begin with we'd probably like to use a fresh tree (with the same ->rdata ) - - //int minPossibleParsimonyScore = 3*(endCharacter - startCharacter + 1); - // TODO calculate actual minimum homoplasy score - //(the above calculation assumes that A,C,G,T appear in each character, and so the min parsimony score of each character is 3) - int minPossibleParsimonyScore = 0; // simplest way to ensure no negative homoplasy ratios while the min parsimony score is unknown - - int returnVal = parsimonyScore - minPossibleParsimonyScore; - return returnVal; - //return 0; + assert (startCharacter <= endCharacter); + + int parsimonyScore = makeParsimonyTreeFastDNA(tr, adef, startCharacter, endCharacter); + // TODO: makeParsimonyTreeFastDNA changes the tree, right? + // Which means that after the first time we call makeParsimonyTreeFastDNA, subsequent calls won't be starting with a 'clean' tree. + // This might be useful later on, but to begin with we'd probably like to use a fresh tree (with the same ->rdata ) + + //int minPossibleParsimonyScore = 3*(endCharacter - startCharacter + 1); + // TODO calculate actual minimum homoplasy score + //(the above calculation assumes that A,C,G,T appear in each character, and so the min parsimony score of each character is 3) + int minPossibleParsimonyScore = 0; // simplest way to ensure no negative homoplasy ratios while the min parsimony score is unknown + + int returnVal = parsimonyScore - minPossibleParsimonyScore; + return returnVal; + //return 0; } -static void getOptBlockPartition(tree *tr, analdef *adef) +static void getOptBlockPartition(tree *tr, analdef *adef) { - int debugOutput = 0; + int debugOutput = 0; - int numsp = tr->rdta->numsp; // number of species - int sites = tr->rdta->sites; // number of sites - int desiredBlockCount = adef->numberOfBlocks; + int numsp = tr->rdta->numsp; // number of species + int sites = tr->rdta->sites; // number of sites + int desiredBlockCount = adef->numberOfBlocks; - // BP_blockScoreArray[i][j] stores the optimum homoplasy score of the block on characters i to j (inclusive) - //int **BP_blockScoreArray = get2dIntArray(sites, sites); + // BP_blockScoreArray[i][j] stores the optimum homoplasy score of the block on characters i to j (inclusive) + //int **BP_blockScoreArray = get2dIntArray(sites, sites); - int BP_blockScoreArray[sites][sites];// = get2dIntArray(sites, sites); + int BP_blockScoreArray[sites][sites];// = get2dIntArray(sites, sites); - int impossiblyHighHomoplasyScore = numsp*sites; + int impossiblyHighHomoplasyScore = numsp*sites; -// Entry BP_DPHomoplasyRatioLookupTable[b][i][j] for i<= j is the minimum value of -// the homoplasy ratio of any block partition P of character 0 to j, +// Entry BP_DPHomoplasyRatioLookupTable[b][i][j] for i<= j is the minimum value of +// the homoplasy ratio of any block partition P of character 0 to j, // for which the last block is from i to j and there are b+1 blocks in total; -// where "homoplasy ratio of P" means the maximum homoplasy ratio of -// any block in P. A value of Float.POSITIVE_INFINITY means there is no solution +// where "homoplasy ratio of P" means the maximum homoplasy ratio of +// any block in P. A value of Float.POSITIVE_INFINITY means there is no solution // for this set of values (e.g. if i > j) - //float ***BP_DPHomoplasyRatioLookupTable = get3dFloatArray(desiredBlockCount, sites, sites); +//float ***BP_DPHomoplasyRatioLookupTable = get3dFloatArray(desiredBlockCount, sites, sites); - float BP_DPHomoplasyRatioLookupTable[desiredBlockCount][sites][sites]; + float BP_DPHomoplasyRatioLookupTable[desiredBlockCount][sites][sites]; // Used to reconstruct an optimal solution. -// BP_backtrackTable[b][i][j] is the value i' such that an +// BP_backtrackTable[b][i][j] is the value i' such that an // optimal block partition with b+1 blocks and last block [i,j] // has [i',j-1] as its second-to-last block. - //int ***BP_backtrackTable = get3dIntArray(desiredBlockCount, sites, sites); +//int ***BP_backtrackTable = get3dIntArray(desiredBlockCount, sites, sites); - int BP_backtrackTable[desiredBlockCount][sites][sites]; + int BP_backtrackTable[desiredBlockCount][sites][sites]; - //printf("populating BP_blockScoreArray\n"); - // Populate BP_blockScoreArray, which stores the optimum homoplasy score every contiguous block of characters - // Todo there is stuff to do with declaring data structures and assigning memory space and stuff like that that I am super not doing at the moment. + //printf("populating BP_blockScoreArray\n"); + // Populate BP_blockScoreArray, which stores the optimum homoplasy score every contiguous block of characters + // Todo there is stuff to do with declaring data structures and assigning memory space and stuff like that that I am super not doing at the moment. - for (int i = 0; i < sites; i++) - { - for (int j = 0; j < sites; j++) + for (int i = 0; i < sites; i++) { - if (i > j) - BP_blockScoreArray[i][j] = INT_MAX; - else - BP_blockScoreArray[i][j] = getOptHomoplasyOfContiguousBlock(tr, adef, i, j); - - } - if (debugOutput > 0) - { - printf("Calculated parsimony score for blocks [%d -- j] for all j\n",i); - } - - } - - //printf("populating lookup table.\n"); - - for (int i = 0; i < sites; i++) - { - - for (int j = 0; j < sites; j++) - { - for (int b = 0; b < desiredBlockCount; b++) - { - //printf("%f\n ",BP_DPHomoplasyRatioLookupTable[0][0][0]); - float val = FLT_MAX; - int bestPrevI = -1; - - if (b == 0 && i > 0) + for (int j = 0; j < sites; j++) { - val = FLT_MAX; // if first block, i should be 0 - } - if (b > 0 && i == 0) - { - val = FLT_MAX; // if not first block, i should be bigger than 0 - } + if (i > j) + BP_blockScoreArray[i][j] = INT_MAX; + else + BP_blockScoreArray[i][j] = getOptHomoplasyOfContiguousBlock(tr, adef, i, j); - if (b == 0 && i == 0) // if first block, take homoplasy ratio of that block. - { - val = getHomoplasyRatio(i,j, BP_blockScoreArray[i][j]); - //printf("%d %d %d = %f\n ",b, i,j,val); } - if (b > 0 && i > 0) + if (debugOutput > 0) { - if (i > j) // i should be at most j. - { - val = FLT_MAX; - } - else - { - - //find the minimum possible homoplasy score for any b-1 block partition up to this point - float currentMinimum = FLT_MAX; - int lastj = i-1; - for(int lasti = 0; lasti <= lastj; lasti++) - { - float tempVal = BP_DPHomoplasyRatioLookupTable[b-1][lasti][lastj]; - - if (tempVal < currentMinimum ) - { - currentMinimum = tempVal; - // update the bestPrevI value because we have a new best block partition - bestPrevI = lasti; - } - } - // currentMinimum is now the optimum homoplasy ratio for any block partition of 0 to i-1 - // using b-1 blocks - float valCurrentBlock = getHomoplasyRatio(i,j, BP_blockScoreArray[i][j]); - // If the current block has worse homoplasy ratio than what we could - // get up to this point, that ratio is now the best we can do for a block partition of this type. - if (currentMinimum < valCurrentBlock ) - { - val = valCurrentBlock; - } - // Otherwise, the best we could do up to this point is still the best we can do - else - { - val = currentMinimum; - } - } + printf("Calculated parsimony score for blocks [%d -- j] for all j\n",i); } - // populate the lookup table with the calculated value - BP_DPHomoplasyRatioLookupTable[b][i][j] = val; + } + + //printf("populating lookup table.\n"); - // also populate the backtrack table so we can reconstruct a solution later - BP_backtrackTable[b][i][j] = bestPrevI; + for (int i = 0; i < sites; i++) + { - if (debugOutput > 0) + for (int j = 0; j < sites; j++) { - printf("BP_DPHomoplasyRatioLookupTable[%d]",b); - printf("[%d]",i); - printf("[%d]",j); - printf("= %f ",val); - printf("BP_backtrackTable[%d]",b); - printf("[%d]",i); - printf("[%d]",j); - printf("= %d\n",bestPrevI); + for (int b = 0; b < desiredBlockCount; b++) + { + //printf("%f\n ",BP_DPHomoplasyRatioLookupTable[0][0][0]); + float val = FLT_MAX; + int bestPrevI = -1; + + if (b == 0 && i > 0) + { + val = FLT_MAX; // if first block, i should be 0 + } + if (b > 0 && i == 0) + { + val = FLT_MAX; // if not first block, i should be bigger than 0 + } + + if (b == 0 && i == 0) // if first block, take homoplasy ratio of that block. + { + val = getHomoplasyRatio(i,j, BP_blockScoreArray[i][j]); + //printf("%d %d %d = %f\n ",b, i,j,val); + } + if (b > 0 && i > 0) + { + if (i > j) // i should be at most j. + { + val = FLT_MAX; + } + else + { + + //find the minimum possible homoplasy score for any b-1 block partition up to this point + float currentMinimum = FLT_MAX; + int lastj = i-1; + for(int lasti = 0; lasti <= lastj; lasti++) + { + float tempVal = BP_DPHomoplasyRatioLookupTable[b-1][lasti][lastj]; + + if (tempVal < currentMinimum ) + { + currentMinimum = tempVal; + // update the bestPrevI value because we have a new best block partition + bestPrevI = lasti; + } + } + // currentMinimum is now the optimum homoplasy ratio for any block partition of 0 to i-1 + // using b-1 blocks + float valCurrentBlock = getHomoplasyRatio(i,j, BP_blockScoreArray[i][j]); + // If the current block has worse homoplasy ratio than what we could + // get up to this point, that ratio is now the best we can do for a block partition of this type. + if (currentMinimum < valCurrentBlock ) + { + val = valCurrentBlock; + } + // Otherwise, the best we could do up to this point is still the best we can do + else + { + val = currentMinimum; + } + } + } + + // populate the lookup table with the calculated value + BP_DPHomoplasyRatioLookupTable[b][i][j] = val; + + // also populate the backtrack table so we can reconstruct a solution later + BP_backtrackTable[b][i][j] = bestPrevI; + + if (debugOutput > 0) + { + printf("BP_DPHomoplasyRatioLookupTable[%d]",b); + printf("[%d]",i); + printf("[%d]",j); + printf("= %f ",val); + printf("BP_backtrackTable[%d]",b); + printf("[%d]",i); + printf("[%d]",j); + printf("= %d\n",bestPrevI); + } + } } - } } - } - // Okay, now we have the DP homoplasy lookup table and also the backtrack table. it remains to just find the best block partition using these tabes! + // Okay, now we have the DP homoplasy lookup table and also the backtrack table. it remains to just find the best block partition using these tabes! - printf("determining optimal block partition.\n"); + printf("determining optimal block partition.\n"); - // Given DPHomoplasyRatioLookupTable, it remains to find the value corresponding a block partition of the full set of characters (i.e. j = characterCount-1) that has minimum homoplasy ratio. - float currentOptimum = FLT_MAX; - // Find the minimum value over all entries DPHomoplasyRatioLookupTable[b][i][j] with j = characterCount-1. - // (Note that b may be smaller than desiredBlockCount-1, if splitting into fewer than desiredBlockcount blocks - // actually gives better homoplasy ratio. - int optFinalI = -1; - int optFinalB = -1; - for(int b = 0; b < desiredBlockCount; b++) - { - for(int i = 0; i < sites; i++) + // Given DPHomoplasyRatioLookupTable, it remains to find the value corresponding a block partition of the full set of characters (i.e. j = characterCount-1) that has minimum homoplasy ratio. + float currentOptimum = FLT_MAX; + // Find the minimum value over all entries DPHomoplasyRatioLookupTable[b][i][j] with j = characterCount-1. + // (Note that b may be smaller than desiredBlockCount-1, if splitting into fewer than desiredBlockcount blocks + // actually gives better homoplasy ratio. + int optFinalI = -1; + int optFinalB = -1; + for(int b = 0; b < desiredBlockCount; b++) { - float tempVal = BP_DPHomoplasyRatioLookupTable[b][i][sites-1]; - if (tempVal < currentOptimum ) - { - currentOptimum = tempVal; - optFinalI = i; - optFinalB = b; - } + for(int i = 0; i < sites; i++) + { + float tempVal = BP_DPHomoplasyRatioLookupTable[b][i][sites-1]; + if (tempVal < currentOptimum ) + { + currentOptimum = tempVal; + optFinalI = i; + optFinalB = b; + } + } } - } - // Store the minimum value found - float BP_optHomoplasyRatio = currentOptimum; - int optBlockCount = optFinalB + 1; - // Now use backtracking to reconstruct the details of an optimal solution + // Store the minimum value found + float BP_optHomoplasyRatio = currentOptimum; + int optBlockCount = optFinalB + 1; + // Now use backtracking to reconstruct the details of an optimal solution - int BP_optBlockStarts[optBlockCount]; // List of the first elements of each block in an optimal block partition - int BP_optBlockEnds[optBlockCount]; // List of the last elements of each block in an optimal block partition - //int BP_optTrees[optFinalB + 1]; // List of (the indices of?????) optimal trees for each block in an optimal block partition + int BP_optBlockStarts[optBlockCount]; // List of the first elements of each block in an optimal block partition + int BP_optBlockEnds[optBlockCount]; // List of the last elements of each block in an optimal block partition + //int BP_optTrees[optFinalB + 1]; // List of (the indices of?????) optimal trees for each block in an optimal block partition - // The first block we store information on is actually the last block, which we know ends with the last character. - int currentBlockStart = optFinalI; - int currentBlockEnd = sites - 1; - //int currentBlockCount = optFinalB; + // The first block we store information on is actually the last block, which we know ends with the last character. + int currentBlockStart = optFinalI; + int currentBlockEnd = sites - 1; + //int currentBlockCount = optFinalB; - // TEMP COMMENTING OUT TO GET AROUND ANOTHER BUG - - // Work backwards, finding the previous block and recording its details, until we reach the first block - for (int b = optBlockCount-1; b >= 0; b--) - { - // printf("heeey%d\n", b); - if (b >0) - { - assert (currentBlockStart <= currentBlockEnd); - assert (currentBlockStart >= 0); - assert (currentBlockEnd >= 0); - } - // append the details of the current block to their respective lists. - BP_optBlockStarts[b] = currentBlockStart; - BP_optBlockEnds[b] = currentBlockEnd; - //printf("hooo %d\n", currentBlockStart); - //printf("hrrrrrro %d\n", currentBlockEnd); - //BP_optTrees.add(BP_optTreePerBlock[currentBlockStart][currentBlockEnd]); + // TEMP COMMENTING OUT TO GET AROUND ANOTHER BUG - int previousBlockStart = BP_backtrackTable[b][currentBlockStart][currentBlockEnd]; - int previousBlockEnd = currentBlockStart - 1; - //int previousBlockCount= currentBlockCount - 1; - //printf("haaaaa\n"); + // Work backwards, finding the previous block and recording its details, until we reach the first block + for (int b = optBlockCount-1; b >= 0; b--) + { + // printf("heeey%d\n", b); + if (b >0) + { + assert (currentBlockStart <= currentBlockEnd); + assert (currentBlockStart >= 0); + assert (currentBlockEnd >= 0); + } + // append the details of the current block to their respective lists. + BP_optBlockStarts[b] = currentBlockStart; + BP_optBlockEnds[b] = currentBlockEnd; + //printf("hooo %d\n", currentBlockStart); + //printf("hrrrrrro %d\n", currentBlockEnd); + //BP_optTrees.add(BP_optTreePerBlock[currentBlockStart][currentBlockEnd]); + int previousBlockStart = BP_backtrackTable[b][currentBlockStart][currentBlockEnd]; + int previousBlockEnd = currentBlockStart - 1; + //int previousBlockCount= currentBlockCount - 1; + //printf("haaaaa\n"); - // previous block becomes the new current block - currentBlockStart = previousBlockStart; - currentBlockEnd = previousBlockEnd; - //currentBlockCount = previousBlockCount; - } + // previous block becomes the new current block + currentBlockStart = previousBlockStart; + currentBlockEnd = previousBlockEnd; + //currentBlockCount = previousBlockCount; + } - // Tada! We have found the optimal block partition. Hopefully. Time to report it! - printf("Optimal block partition: %d blocks", optFinalB+1); + // Tada! We have found the optimal block partition. Hopefully. Time to report it! + printf("Optimal block partition: %d blocks", optFinalB+1); - for (int h = 0; h <= optFinalB; h++) - { - int tempBlockStart = BP_optBlockStarts[h]; - int tempBlockEnd = BP_optBlockEnds[h]; - printf("["); - printf("%d", tempBlockStart); - printf("---"); - printf("%d", tempBlockEnd); - printf("]"); - if (h < optFinalB) - { - printf(", "); - } - else + for (int h = 0; h <= optFinalB; h++) { - printf("\n"); + + int tempBlockStart = BP_optBlockStarts[h]; + int tempBlockEnd = BP_optBlockEnds[h]; + printf("["); + printf("%d", tempBlockStart); + printf("---"); + printf("%d", tempBlockEnd); + printf("]"); + if (h < optFinalB) + { + printf(", "); + } + else + { + printf("\n"); + } } - } - printf("Maximum homoplasy ratio: %f\n", BP_optHomoplasyRatio); - + printf("Maximum homoplasy ratio: %f\n", BP_optHomoplasyRatio); + - // TODO: report back more details, like the homoplasy score / optimal tree for each block? + // TODO: report back more details, like the homoplasy score / optimal tree for each block? } @@ -1707,59 +1707,57 @@ static void getOptBlockPartition(tree *tr, analdef *adef) int main (int argc, char *argv[]) { - rawdata *rdta; - //cruncheddata *cdta; //not needed anymore - tree *tr; - analdef *adef; + rawdata *rdta; + //cruncheddata *cdta; //not needed anymore + tree *tr; + analdef *adef; - - - adef = (analdef *)malloc(sizeof(analdef)); - rdta = (rawdata *)malloc(sizeof(rawdata)); - // cdta = (cruncheddata *)malloc(sizeof(cruncheddata)); //not needed anymore - tr = (tree *)malloc(sizeof(tree)); - /* initialize lookup table for fast bit counter */ - compute_bits_in_16bits(); - - initAdef(adef); + adef = (analdef *)malloc(sizeof(analdef)); + rdta = (rawdata *)malloc(sizeof(rawdata)); + // cdta = (cruncheddata *)malloc(sizeof(cruncheddata)); //not needed anymore + tr = (tree *)malloc(sizeof(tree)); - - get_args(argc,argv, adef); - + /* initialize lookup table for fast bit counter */ - getinput(adef, rdta, tr); // , cdta //not needed anymore - + compute_bits_in_16bits(); - //makeFileNames(); //we do not want to print for now - + initAdef(adef); - makevalues(rdta, tr); //, cdta //not needed anymore - - //we do not want to print for now - //printBothOpen("\n%s version %s a fast open-source code for builiding DNA parsimony start trees\n\n", programName, programVersion); - //printBothOpen("Released under GNU GPL in %s by Alexandros Stamatakis\n\n", programDate); - //printBothOpen("Alignment has %d sites and %d taxa\n\n", tr->cdta->endsite, tr->mxtips); - //printBothOpen("%d randomized stepwise addition order parsimony trees with a couple of SPR moves will be computed\n\n", adef->numberOfTrees); + get_args(argc,argv, adef); + getinput(adef, rdta, tr); // , cdta //not needed anymore - //int minParsimonyScoreOverRuns = makeParsimonyTreeFastDNA(tr, adef, 0, rdta->sites); //rdta->sites); - //printf ("doing a test\n"); - //int minParsimonyScoreOverRuns = makeParsimonyTreeFastDNA(tr, adef, 0, rdta->sites - 1); //rdta->sites); - //printf ("test done\n"); - //minParsimonyScoreOverRuns = makeParsimonyTreeFastDNA(tr, adef, 0, 0); //rdta->sites); + //makeFileNames(); //we do not want to print for now - getOptBlockPartition(tr, adef); + makevalues(rdta, tr); //, cdta //not needed anymore - //printf("Best parsimony score %u\n", minParsimonyScoreOverRuns); + //we do not want to print for now + //printBothOpen("\n%s version %s a fast open-source code for builiding DNA parsimony start trees\n\n", programName, programVersion); + //printBothOpen("Released under GNU GPL in %s by Alexandros Stamatakis\n\n", programDate); + //printBothOpen("Alignment has %d sites and %d taxa\n\n", tr->cdta->endsite, tr->mxtips); + //printBothOpen("%d randomized stepwise addition order parsimony trees with a couple of SPR moves will be computed\n\n", adef->numberOfTrees); - return 0; -} + //int minParsimonyScoreOverRuns = makeParsimonyTreeFastDNA(tr, adef, 0, rdta->sites); //rdta->sites); + //printf ("doing a test\n"); + //int minParsimonyScoreOverRuns = makeParsimonyTreeFastDNA(tr, adef, 0, rdta->sites - 1); //rdta->sites); + //printf ("test done\n"); + //minParsimonyScoreOverRuns = makeParsimonyTreeFastDNA(tr, adef, 0, 0); //rdta->sites); + + + getOptBlockPartition(tr, adef); + + + //printf("Best parsimony score %u\n", minParsimonyScoreOverRuns); + + + return 0; +} diff --git a/CUTAL.h b/CUTAL.h index da5031a..bb476ce 100644 --- a/CUTAL.h +++ b/CUTAL.h @@ -1,11 +1,11 @@ -/* CUTAL (version 0.1) a program for cutting alignement into non-recombinant blocks under - * a parsimony framework +/* CUTAL (version 0.1) a program for cutting alignement into non-recombinant blocks under + * a parsimony framework * Copyright May 2018 by Celine Scornavacca and Mark Jones - * - * The program is heavily based on Parsimonator-1.0.2 by Alexandros Stamatakis, which was partially + * + * The program is heavily based on Parsimonator-1.0.2 by Alexandros Stamatakis, which was partially * derived from fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen - * - * and + * + * and * * programs of the PHYLIP package by Joe Felsenstein. * @@ -18,9 +18,9 @@ * 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. - * * - * For any other enquiries send an email to celine.scornavacca@umontpellier.fr and + * + * For any other enquiries send an email to celine.scornavacca@umontpellier.fr and * markelliotlloyd@gmail.com * */ @@ -90,10 +90,10 @@ typedef unsigned int parsimonyNumber; typedef struct ratec { - double accumulatedSiteLikelihood; - double rate; + double accumulatedSiteLikelihood; + double rate; } - rateCategorize; +rateCategorize; @@ -114,13 +114,13 @@ struct noderec; typedef struct noderec { - - struct noderec *next; - struct noderec *back; - int number; - char x; + + struct noderec *next; + struct noderec *back; + int number; + char x; } - node, *nodeptr; +node, *nodeptr; @@ -128,21 +128,21 @@ typedef struct noderec typedef struct { - int numsp; - int sites; - unsigned char **y; - unsigned char *y0; - unsigned char *yBUF; - int *wgt; + int numsp; + int sites; + unsigned char **y; + unsigned char *y0; + unsigned char *yBUF; + int *wgt; } rawdata; typedef struct { - int *alias; /* site representing a pattern */ - int *aliaswgt; /* weight by pattern */ - int *rateCategory; - int endsite; /* # of sequence patterns */ - double *patrat; /* rates per pattern */ - double *patratStored; + int *alias; /* site representing a pattern */ + int *aliaswgt;/* weight by pattern */ + int *rateCategory; + int endsite; /* # of sequence patterns */ + double *patrat;/* rates per pattern */ + double *patratStored; } cruncheddata; @@ -152,32 +152,32 @@ typedef struct { -typedef struct +typedef struct { - int left; - int right; - double likelihood; + int left; + int right; + double likelihood; } lhEntry; -typedef struct +typedef struct { - int count; - int size; - lhEntry *entries; + int count; + int size; + lhEntry *entries; } lhList; -typedef struct List_{ - void *value; - struct List_ *next; +typedef struct List_ { + void *value; + struct List_ *next; } List; struct stringEnt { - int nodeNumber; - char *word; - struct stringEnt *next; + int nodeNumber; + char *word; + struct stringEnt *next; }; typedef struct stringEnt stringEntry; @@ -186,51 +186,51 @@ typedef unsigned int hashNumberType; typedef struct { - hashNumberType tableSize; - stringEntry **table; + hashNumberType tableSize; + stringEntry **table; } - stringHashtable; +stringHashtable; typedef struct { - - - - parsimonyNumber **parsimonyState_A; - parsimonyNumber **parsimonyState_C; - parsimonyNumber **parsimonyState_G; - parsimonyNumber **parsimonyState_T; - unsigned int *parsimonyScore; - int *ti; - unsigned int compressedWidth; - - unsigned char **yVector; - - - - - stringHashtable *nameHash; - - - node **nodep; - node *start; - int mxtips; - - int ntips; - int nextnode; - - - rawdata *rdta; - //cruncheddata *cdta; //not needed anymore - - char **nameList; - char *tree_string; - int *nodesInTree; - - int treeStringLength; - unsigned int bestParsimony; - - nodeptr removeNode; - nodeptr insertNode; + + + + parsimonyNumber **parsimonyState_A; + parsimonyNumber **parsimonyState_C; + parsimonyNumber **parsimonyState_G; + parsimonyNumber **parsimonyState_T; + unsigned int *parsimonyScore; + int *ti; + unsigned int compressedWidth; + + unsigned char **yVector; + + + + + stringHashtable *nameHash; + + + node **nodep; + node *start; + int mxtips; + + int ntips; + int nextnode; + + + rawdata *rdta; + //cruncheddata *cdta; //not needed anymore + + char **nameList; + char *tree_string; + int *nodesInTree; + + int treeStringLength; + unsigned int bestParsimony; + + nodeptr removeNode; + nodeptr insertNode; } tree; @@ -249,10 +249,10 @@ typedef struct { typedef struct { - long parsimonySeed; - boolean restart; - int numberOfTrees; - int numberOfBlocks; + long parsimonySeed; + boolean restart; + int numberOfTrees; + int numberOfBlocks; } analdef; diff --git a/fastDNAparsimony.c b/fastDNAparsimony.c index 7c16daf..f91cb1c 100644 --- a/fastDNAparsimony.c +++ b/fastDNAparsimony.c @@ -1,11 +1,11 @@ -/* CUTAL (version 0.1) a program for cutting alignement into non-recombinant blocks under - * a parsimony framework +/* CUTAL (version 0.1) a program for cutting alignement into non-recombinant blocks under + * a parsimony framework * Copyright May 2018 by Celine Scornavacca and Mark Jones - * - * The program is heavily based on Parsimonator-1.0.2 by Alexandros Stamatakis, which was partially + * + * The program is heavily based on Parsimonator-1.0.2 by Alexandros Stamatakis, which was partially * derived from fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen - * - * and + * + * and * * programs of the PHYLIP package by Joe Felsenstein. * @@ -18,9 +18,9 @@ * 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. - * * - * For any other enquiries send an email to celine.scornavacca@umontpellier.fr and + * + * For any other enquiries send an email to celine.scornavacca@umontpellier.fr and * markelliotlloyd@gmail.com * */ @@ -30,12 +30,12 @@ #include #include #include -#include +#include #endif #include #include -#include +#include #include #include #include @@ -46,7 +46,7 @@ #include #include - + #endif #ifdef __AVX @@ -60,15 +60,15 @@ #include "CUTAL.h" extern char run_id[128]; -extern char seq_file[1024]; -extern char tree_file[1024]; -extern char resultFileName[1024]; -extern char infoFileName[1024]; -extern char randomFileName[1024]; +extern char seq_file[1024]; +extern char tree_file[1024]; +extern char resultFileName[1024]; +extern char infoFileName[1024]; +extern char randomFileName[1024]; -const unsigned int mask32[32] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536, 131072, - 262144, 524288, 1048576, 2097152, 4194304, 8388608, 16777216, 33554432, 67108864, 134217728, - 268435456, 536870912, 1073741824, 2147483648U}; +const unsigned int mask32[32] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536, 131072, + 262144, 524288, 1048576, 2097152, 4194304, 8388608, 16777216, 33554432, 67108864, 134217728, + 268435456, 536870912, 1073741824, 2147483648U}; /* vector-specific stuff */ @@ -108,32 +108,32 @@ const unsigned int mask32[32] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 20 #define BYTE_ALIGNMENT 16 #endif -static void *malloc_aligned(size_t size, size_t align) +static void *malloc_aligned(size_t size, size_t align) { - void *ptr = (void *)NULL; - int res; - + void *ptr = (void *)NULL; + int res; + #if defined (__APPLE__) - /* - presumably malloc on MACs always returns - a 16-byte aligned pointer - */ + /* + presumably malloc on MACs always returns + a 16-byte aligned pointer + */ + + ptr = malloc(size); - ptr = malloc(size); - - if(ptr == (void*)NULL) - assert(0); + if(ptr == (void*)NULL) + assert(0); #else - res = posix_memalign( &ptr, align, size ); + res = posix_memalign( &ptr, align, size ); + + if(res != 0) + assert(0); +#endif - if(res != 0) - assert(0); -#endif - - return ptr; + return ptr; } @@ -143,52 +143,52 @@ static void *malloc_aligned(size_t size, size_t align) static void hookupParsimony(nodeptr p, nodeptr q) { - p->back = q; - q->back = p; + p->back = q; + q->back = p; } static void getxnodeLocal (nodeptr p) { - nodeptr s; + nodeptr s; - if((s = p->next)->x || (s = s->next)->x) - { - p->x = s->x; - s->x = 0; - } + if((s = p->next)->x || (s = s->next)->x) + { + p->x = s->x; + s->x = 0; + } } static void computeTraversalInfoParsimony(nodeptr p, int *ti, int *counter, int maxTips, boolean full) -{ - nodeptr - q = p->next->back, - r = p->next->next->back; - - if(! p->x) - getxnodeLocal(p); - - if(full) - { - if(q->number > maxTips) - computeTraversalInfoParsimony(q, ti, counter, maxTips, full); - - if(r->number > maxTips) - computeTraversalInfoParsimony(r, ti, counter, maxTips, full); - } - else - { - if(q->number > maxTips && !q->x) - computeTraversalInfoParsimony(q, ti, counter, maxTips, full); - - if(r->number > maxTips && !r->x) - computeTraversalInfoParsimony(r, ti, counter, maxTips, full); - } - - - ti[*counter] = p->number; - ti[*counter + 1] = q->number; - ti[*counter + 2] = r->number; - *counter = *counter + 4; +{ + nodeptr + q = p->next->back, + r = p->next->next->back; + + if(!p->x) + getxnodeLocal(p); + + if(full) + { + if(q->number > maxTips) + computeTraversalInfoParsimony(q, ti, counter, maxTips, full); + + if(r->number > maxTips) + computeTraversalInfoParsimony(r, ti, counter, maxTips, full); + } + else + { + if(q->number > maxTips && !q->x) + computeTraversalInfoParsimony(q, ti, counter, maxTips, full); + + if(r->number > maxTips && !r->x) + computeTraversalInfoParsimony(r, ti, counter, maxTips, full); + } + + + ti[*counter] = p->number; + ti[*counter + 1] = q->number; + ti[*counter + 2] = r->number; + *counter = *counter + 4; } @@ -201,26 +201,26 @@ static void computeTraversalInfoParsimony(nodeptr p, int *ti, int *counter, int static inline unsigned int populationCount(INT_TYPE v_N) { - unsigned int - res[INTS_PER_VECTOR] __attribute__ ((aligned (BYTE_ALIGNMENT))); - - unsigned int - i, - a = 0; - - VECTOR_STORE((CAST)res, v_N); - - for(i = 0; i < INTS_PER_VECTOR; i++) - a += BIT_COUNT(res[i]); - - return a; + unsigned int + res[INTS_PER_VECTOR] __attribute__ ((aligned (BYTE_ALIGNMENT))); + + unsigned int + i, + a = 0; + + VECTOR_STORE((CAST)res, v_N); + + for(i = 0; i < INTS_PER_VECTOR; i++) + a += BIT_COUNT(res[i]); + + return a; } #else static inline unsigned int populationCount(unsigned int n) { - return BIT_COUNT(n); + return BIT_COUNT(n); } #endif @@ -231,152 +231,152 @@ static inline unsigned int populationCount(unsigned int n) static void newviewParsimonyIterativeFast(tree *tr) -{ - INT_TYPE - allOne = SET_ALL_BITS_ONE; - - int - *ti = tr->ti, - count = ti[0], - index; - - unsigned int - width = tr->compressedWidth; - - - - for(index = 4; index < count; index += 4) - { - int - pNumber = ti[index], - qNumber = ti[index + 1], - rNumber = ti[index + 2]; - - parsimonyNumber - *leftState_A = tr->parsimonyState_A[qNumber], - *rightState_A = tr->parsimonyState_A[rNumber], - *thisState_A = tr->parsimonyState_A[pNumber], - - *leftState_C = tr->parsimonyState_C[qNumber], - *rightState_C = tr->parsimonyState_C[rNumber], - *thisState_C = tr->parsimonyState_C[pNumber], - - *leftState_G = tr->parsimonyState_G[qNumber], - *rightState_G = tr->parsimonyState_G[rNumber], - *thisState_G = tr->parsimonyState_G[pNumber], - - *leftState_T = tr->parsimonyState_T[qNumber], - *rightState_T = tr->parsimonyState_T[rNumber], - *thisState_T = tr->parsimonyState_T[pNumber]; - - - - unsigned int - i, - ts = 0; - +{ + INT_TYPE + allOne = SET_ALL_BITS_ONE; + + int + *ti = tr->ti, + count = ti[0], + index; + + unsigned int + width = tr->compressedWidth; + + + + for(index = 4; index < count; index += 4) + { + int + pNumber = ti[index], + qNumber = ti[index + 1], + rNumber = ti[index + 2]; + + parsimonyNumber + *leftState_A = tr->parsimonyState_A[qNumber], + *rightState_A = tr->parsimonyState_A[rNumber], + *thisState_A = tr->parsimonyState_A[pNumber], + + *leftState_C = tr->parsimonyState_C[qNumber], + *rightState_C = tr->parsimonyState_C[rNumber], + *thisState_C = tr->parsimonyState_C[pNumber], + + *leftState_G = tr->parsimonyState_G[qNumber], + *rightState_G = tr->parsimonyState_G[rNumber], + *thisState_G = tr->parsimonyState_G[pNumber], + + *leftState_T = tr->parsimonyState_T[qNumber], + *rightState_T = tr->parsimonyState_T[rNumber], + *thisState_T = tr->parsimonyState_T[pNumber]; + + + + unsigned int + i, + ts = 0; + #pragma omp parallel for reduction(+:ts) - for(i = 0; i < width; i += INTS_PER_VECTOR) - { - INT_TYPE - s_r, s_l, v_N, - l_A, l_C, l_G, l_T, - v_A, v_C, v_G, v_T; - - s_l = VECTOR_LOAD((CAST)(&leftState_A[i])); - s_r = VECTOR_LOAD((CAST)(&rightState_A[i])); - l_A = VECTOR_BIT_AND(s_l, s_r); - v_A = VECTOR_BIT_OR(s_l, s_r); - - s_l = VECTOR_LOAD((CAST)(&leftState_C[i])); - s_r = VECTOR_LOAD((CAST)(&rightState_C[i])); - l_C = VECTOR_BIT_AND(s_l, s_r); - v_C = VECTOR_BIT_OR(s_l, s_r); - - s_l = VECTOR_LOAD((CAST)(&leftState_G[i])); - s_r = VECTOR_LOAD((CAST)(&rightState_G[i])); - l_G = VECTOR_BIT_AND(s_l, s_r); - v_G = VECTOR_BIT_OR(s_l, s_r); - - s_l = VECTOR_LOAD((CAST)(&leftState_T[i])); - s_r = VECTOR_LOAD((CAST)(&rightState_T[i])); - l_T = VECTOR_BIT_AND(s_l, s_r); - v_T = VECTOR_BIT_OR(s_l, s_r); - - v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T)); - - VECTOR_STORE((CAST)(&thisState_A[i]), VECTOR_BIT_OR(l_A, VECTOR_AND_NOT(v_N, v_A))); - VECTOR_STORE((CAST)(&thisState_C[i]), VECTOR_BIT_OR(l_C, VECTOR_AND_NOT(v_N, v_C))); - VECTOR_STORE((CAST)(&thisState_G[i]), VECTOR_BIT_OR(l_G, VECTOR_AND_NOT(v_N, v_G))); - VECTOR_STORE((CAST)(&thisState_T[i]), VECTOR_BIT_OR(l_T, VECTOR_AND_NOT(v_N, v_T))); - - v_N = VECTOR_AND_NOT(v_N, allOne); - - ts += populationCount(v_N); - } - - tr->parsimonyScore[pNumber] = ts + tr->parsimonyScore[rNumber] + tr->parsimonyScore[qNumber]; - } + for(i = 0; i < width; i += INTS_PER_VECTOR) + { + INT_TYPE + s_r, s_l, v_N, + l_A, l_C, l_G, l_T, + v_A, v_C, v_G, v_T; + + s_l = VECTOR_LOAD((CAST)(&leftState_A[i])); + s_r = VECTOR_LOAD((CAST)(&rightState_A[i])); + l_A = VECTOR_BIT_AND(s_l, s_r); + v_A = VECTOR_BIT_OR(s_l, s_r); + + s_l = VECTOR_LOAD((CAST)(&leftState_C[i])); + s_r = VECTOR_LOAD((CAST)(&rightState_C[i])); + l_C = VECTOR_BIT_AND(s_l, s_r); + v_C = VECTOR_BIT_OR(s_l, s_r); + + s_l = VECTOR_LOAD((CAST)(&leftState_G[i])); + s_r = VECTOR_LOAD((CAST)(&rightState_G[i])); + l_G = VECTOR_BIT_AND(s_l, s_r); + v_G = VECTOR_BIT_OR(s_l, s_r); + + s_l = VECTOR_LOAD((CAST)(&leftState_T[i])); + s_r = VECTOR_LOAD((CAST)(&rightState_T[i])); + l_T = VECTOR_BIT_AND(s_l, s_r); + v_T = VECTOR_BIT_OR(s_l, s_r); + + v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T)); + + VECTOR_STORE((CAST)(&thisState_A[i]), VECTOR_BIT_OR(l_A, VECTOR_AND_NOT(v_N, v_A))); + VECTOR_STORE((CAST)(&thisState_C[i]), VECTOR_BIT_OR(l_C, VECTOR_AND_NOT(v_N, v_C))); + VECTOR_STORE((CAST)(&thisState_G[i]), VECTOR_BIT_OR(l_G, VECTOR_AND_NOT(v_N, v_G))); + VECTOR_STORE((CAST)(&thisState_T[i]), VECTOR_BIT_OR(l_T, VECTOR_AND_NOT(v_N, v_T))); + + v_N = VECTOR_AND_NOT(v_N, allOne); + + ts += populationCount(v_N); + } + + tr->parsimonyScore[pNumber] = ts + tr->parsimonyScore[rNumber] + tr->parsimonyScore[qNumber]; + } } static unsigned int evaluateParsimonyIterativeFast(tree *tr) { - INT_TYPE - allOne = SET_ALL_BITS_ONE; - - int - pNumber = tr->ti[1], - qNumber = tr->ti[2]; - - unsigned int - bestScore = tr->bestParsimony, - width = tr->compressedWidth, - i, - sum; - - parsimonyNumber - *rightState_A = tr->parsimonyState_A[pNumber], - *leftState_A = tr->parsimonyState_A[qNumber], - - *rightState_C = tr->parsimonyState_C[pNumber], - *leftState_C = tr->parsimonyState_C[qNumber], - - *rightState_G = tr->parsimonyState_G[pNumber], - *leftState_G = tr->parsimonyState_G[qNumber], - - *rightState_T = tr->parsimonyState_T[pNumber], - *leftState_T = tr->parsimonyState_T[qNumber]; - - - if(tr->ti[0] > 4) - newviewParsimonyIterativeFast(tr); - - sum = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber]; - -#pragma omp parallel for reduction(+:sum) - for(i = 0; i < width; i += INTS_PER_VECTOR) - { - INT_TYPE - l_A = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_A[i])), VECTOR_LOAD((CAST)(&rightState_A[i]))), - l_C = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_C[i])), VECTOR_LOAD((CAST)(&rightState_C[i]))), - l_G = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_G[i])), VECTOR_LOAD((CAST)(&rightState_G[i]))), - l_T = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_T[i])), VECTOR_LOAD((CAST)(&rightState_T[i]))), - v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T)); - - - v_N = VECTOR_AND_NOT(v_N, allOne); - - sum += populationCount(v_N); + INT_TYPE + allOne = SET_ALL_BITS_ONE; + + int + pNumber = tr->ti[1], + qNumber = tr->ti[2]; + + unsigned int + bestScore = tr->bestParsimony, + width = tr->compressedWidth, + i, + sum; + + parsimonyNumber + *rightState_A = tr->parsimonyState_A[pNumber], + *leftState_A = tr->parsimonyState_A[qNumber], + + *rightState_C = tr->parsimonyState_C[pNumber], + *leftState_C = tr->parsimonyState_C[qNumber], + + *rightState_G = tr->parsimonyState_G[pNumber], + *leftState_G = tr->parsimonyState_G[qNumber], + + *rightState_T = tr->parsimonyState_T[pNumber], + *leftState_T = tr->parsimonyState_T[qNumber]; + + + if(tr->ti[0] > 4) + newviewParsimonyIterativeFast(tr); + + sum = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber]; + +#pragma omp parallel for reduction(+:sum) + for(i = 0; i < width; i += INTS_PER_VECTOR) + { + INT_TYPE + l_A = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_A[i])), VECTOR_LOAD((CAST)(&rightState_A[i]))), + l_C = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_C[i])), VECTOR_LOAD((CAST)(&rightState_C[i]))), + l_G = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_G[i])), VECTOR_LOAD((CAST)(&rightState_G[i]))), + l_T = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_T[i])), VECTOR_LOAD((CAST)(&rightState_T[i]))), + v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T)); + + + v_N = VECTOR_AND_NOT(v_N, allOne); + + sum += populationCount(v_N); #ifndef _OPENMP - if(sum >= bestScore) - break; + if(sum >= bestScore) + break; #endif - } - - return sum; + } + + return sum; } @@ -385,131 +385,131 @@ static unsigned int evaluateParsimonyIterativeFast(tree *tr) static void newviewParsimonyIterativeFast(tree *tr) -{ - int - count = tr->ti[0], - *ti = tr->ti, - index; - - unsigned int - width = tr->compressedWidth; - - for(index = 4; index < count; index += 4) - { - int - pNumber = ti[index], - qNumber = ti[index + 1], - rNumber = ti[index + 2]; - - parsimonyNumber - *leftState_A = tr->parsimonyState_A[qNumber], - *rightState_A = tr->parsimonyState_A[rNumber], - *thisState_A = tr->parsimonyState_A[pNumber], - - *leftState_C = tr->parsimonyState_C[qNumber], - *rightState_C = tr->parsimonyState_C[rNumber], - *thisState_C = tr->parsimonyState_C[pNumber], - - *leftState_G = tr->parsimonyState_G[qNumber], - *rightState_G = tr->parsimonyState_G[rNumber], - *thisState_G = tr->parsimonyState_G[pNumber], - - *leftState_T = tr->parsimonyState_T[qNumber], - *rightState_T = tr->parsimonyState_T[rNumber], - *thisState_T = tr->parsimonyState_T[pNumber]; - - register parsimonyNumber - o_A, - o_C, - o_G, - o_T, - t_A, - t_C, - t_G, - t_T, - t_N; - - unsigned int - i, - ts = 0; - - for(i = 0; i < width; i++) - { - t_A = leftState_A[i] & rightState_A[i]; - t_C = leftState_C[i] & rightState_C[i]; - t_G = leftState_G[i] & rightState_G[i]; - t_T = leftState_T[i] & rightState_T[i]; - - o_A = leftState_A[i] | rightState_A[i]; - o_C = leftState_C[i] | rightState_C[i]; - o_G = leftState_G[i] | rightState_G[i]; - o_T = leftState_T[i] | rightState_T[i]; - - t_N = ~(t_A | t_C | t_G | t_T); - - thisState_A[i] = t_A | (t_N & o_A); - thisState_C[i] = t_C | (t_N & o_C); - thisState_G[i] = t_G | (t_N & o_G); - thisState_T[i] = t_T | (t_N & o_T); - - ts += populationCount(t_N); - } - - tr->parsimonyScore[pNumber] = ts + tr->parsimonyScore[rNumber] + tr->parsimonyScore[qNumber]; - } +{ + int + count = tr->ti[0], + *ti = tr->ti, + index; + + unsigned int + width = tr->compressedWidth; + + for(index = 4; index < count; index += 4) + { + int + pNumber = ti[index], + qNumber = ti[index + 1], + rNumber = ti[index + 2]; + + parsimonyNumber + *leftState_A = tr->parsimonyState_A[qNumber], + *rightState_A = tr->parsimonyState_A[rNumber], + *thisState_A = tr->parsimonyState_A[pNumber], + + *leftState_C = tr->parsimonyState_C[qNumber], + *rightState_C = tr->parsimonyState_C[rNumber], + *thisState_C = tr->parsimonyState_C[pNumber], + + *leftState_G = tr->parsimonyState_G[qNumber], + *rightState_G = tr->parsimonyState_G[rNumber], + *thisState_G = tr->parsimonyState_G[pNumber], + + *leftState_T = tr->parsimonyState_T[qNumber], + *rightState_T = tr->parsimonyState_T[rNumber], + *thisState_T = tr->parsimonyState_T[pNumber]; + + register parsimonyNumber + o_A, + o_C, + o_G, + o_T, + t_A, + t_C, + t_G, + t_T, + t_N; + + unsigned int + i, + ts = 0; + + for(i = 0; i < width; i++) + { + t_A = leftState_A[i] & rightState_A[i]; + t_C = leftState_C[i] & rightState_C[i]; + t_G = leftState_G[i] & rightState_G[i]; + t_T = leftState_T[i] & rightState_T[i]; + + o_A = leftState_A[i] | rightState_A[i]; + o_C = leftState_C[i] | rightState_C[i]; + o_G = leftState_G[i] | rightState_G[i]; + o_T = leftState_T[i] | rightState_T[i]; + + t_N = ~(t_A | t_C | t_G | t_T); + + thisState_A[i] = t_A | (t_N & o_A); + thisState_C[i] = t_C | (t_N & o_C); + thisState_G[i] = t_G | (t_N & o_G); + thisState_T[i] = t_T | (t_N & o_T); + + ts += populationCount(t_N); + } + + tr->parsimonyScore[pNumber] = ts + tr->parsimonyScore[rNumber] + tr->parsimonyScore[qNumber]; + } } static unsigned int evaluateParsimonyIterativeFast(tree *tr) { - int - pNumber = tr->ti[1], - qNumber = tr->ti[2]; - - unsigned int - bestScore = tr->bestParsimony, - width = tr->compressedWidth, - i, - sum, - t_A, - t_C, - t_G, - t_T, - t_N; - - parsimonyNumber - *rightState_A = tr->parsimonyState_A[pNumber], - *leftState_A = tr->parsimonyState_A[qNumber], - - *rightState_C = tr->parsimonyState_C[pNumber], - *leftState_C = tr->parsimonyState_C[qNumber], - - *rightState_G = tr->parsimonyState_G[pNumber], - *leftState_G = tr->parsimonyState_G[qNumber], - - *rightState_T = tr->parsimonyState_T[pNumber], - *leftState_T = tr->parsimonyState_T[qNumber]; - - if(tr->ti[0] > 4) - newviewParsimonyIterativeFast(tr); - - sum = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber]; - - for(i = 0; i < width; i++) - { - t_A = leftState_A[i] & rightState_A[i]; - t_C = leftState_C[i] & rightState_C[i]; - t_G = leftState_G[i] & rightState_G[i]; - t_T = leftState_T[i] & rightState_T[i]; - - t_N = ~(t_A | t_C | t_G | t_T); - - sum += populationCount(t_N); - if(sum >= bestScore) - break; - } - - return sum; + int + pNumber = tr->ti[1], + qNumber = tr->ti[2]; + + unsigned int + bestScore = tr->bestParsimony, + width = tr->compressedWidth, + i, + sum, + t_A, + t_C, + t_G, + t_T, + t_N; + + parsimonyNumber + *rightState_A = tr->parsimonyState_A[pNumber], + *leftState_A = tr->parsimonyState_A[qNumber], + + *rightState_C = tr->parsimonyState_C[pNumber], + *leftState_C = tr->parsimonyState_C[qNumber], + + *rightState_G = tr->parsimonyState_G[pNumber], + *leftState_G = tr->parsimonyState_G[qNumber], + + *rightState_T = tr->parsimonyState_T[pNumber], + *leftState_T = tr->parsimonyState_T[qNumber]; + + if(tr->ti[0] > 4) + newviewParsimonyIterativeFast(tr); + + sum = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber]; + + for(i = 0; i < width; i++) + { + t_A = leftState_A[i] & rightState_A[i]; + t_C = leftState_C[i] & rightState_C[i]; + t_G = leftState_G[i] & rightState_G[i]; + t_T = leftState_T[i] & rightState_T[i]; + + t_N = ~(t_A | t_C | t_G | t_T); + + sum += populationCount(t_N); + if(sum >= bestScore) + break; + } + + return sum; } #endif @@ -521,58 +521,58 @@ static unsigned int evaluateParsimonyIterativeFast(tree *tr) static unsigned int evaluateParsimony(tree *tr, nodeptr p, boolean full) { - volatile unsigned int result; - nodeptr q = p->back; - int - *ti = tr->ti, - counter = 4; - - ti[1] = p->number; - ti[2] = q->number; - - if(full) - { - if(p->number > tr->mxtips) - computeTraversalInfoParsimony(p, ti, &counter, tr->mxtips, full); - if(q->number > tr->mxtips) - computeTraversalInfoParsimony(q, ti, &counter, tr->mxtips, full); - } - else - { - if(p->number > tr->mxtips && !p->x) - computeTraversalInfoParsimony(p, ti, &counter, tr->mxtips, full); - if(q->number > tr->mxtips && !q->x) - computeTraversalInfoParsimony(q, ti, &counter, tr->mxtips, full); - } - - - - ti[0] = counter; - - - - result = evaluateParsimonyIterativeFast(tr); - - - - return result; + volatile unsigned int result; + nodeptr q = p->back; + int + *ti = tr->ti, + counter = 4; + + ti[1] = p->number; + ti[2] = q->number; + + if(full) + { + if(p->number > tr->mxtips) + computeTraversalInfoParsimony(p, ti, &counter, tr->mxtips, full); + if(q->number > tr->mxtips) + computeTraversalInfoParsimony(q, ti, &counter, tr->mxtips, full); + } + else + { + if(p->number > tr->mxtips && !p->x) + computeTraversalInfoParsimony(p, ti, &counter, tr->mxtips, full); + if(q->number > tr->mxtips && !q->x) + computeTraversalInfoParsimony(q, ti, &counter, tr->mxtips, full); + } + + + + ti[0] = counter; + + + + result = evaluateParsimonyIterativeFast(tr); + + + + return result; } -static void newviewParsimony(tree *tr, nodeptr p) -{ - if(p->number <= tr->mxtips) - return; - - { - int - counter = 4; - - computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE); - tr->ti[0] = counter; - - newviewParsimonyIterativeFast(tr); - } +static void newviewParsimony(tree *tr, nodeptr p) +{ + if(p->number <= tr->mxtips) + return; + + { + int + counter = 4; + + computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE); + tr->ti[0] = counter; + + newviewParsimonyIterativeFast(tr); + } } @@ -583,724 +583,724 @@ static void newviewParsimony(tree *tr, nodeptr p) static void insertParsimony (tree *tr, nodeptr p, nodeptr q) { - nodeptr r; - - r = q->back; - - hookupParsimony(p->next, q); - hookupParsimony(p->next->next, r); - - newviewParsimony(tr, p); -} + nodeptr r; + + r = q->back; + + hookupParsimony(p->next, q); + hookupParsimony(p->next->next, r); + + newviewParsimony(tr, p); +} static nodeptr buildNewTip (tree *tr, nodeptr p) -{ - nodeptr q; +{ + nodeptr q; + + q = tr->nodep[(tr->nextnode)++]; + hookupParsimony(p, q); + q->next->back = (nodeptr)NULL; + q->next->next->back = (nodeptr)NULL; - q = tr->nodep[(tr->nextnode)++]; - hookupParsimony(p, q); - q->next->back = (nodeptr)NULL; - q->next->next->back = (nodeptr)NULL; - - return q; -} + return q; +} static void buildSimpleTree (tree *tr, int ip, int iq, int ir) -{ - nodeptr p, s; - int i; - - i = MIN(ip, iq); - if (ir < i) i = ir; - tr->start = tr->nodep[i]; - tr->ntips = 3; - p = tr->nodep[ip]; - hookupParsimony(p, tr->nodep[iq]); - s = buildNewTip(tr, tr->nodep[ir]); - insertParsimony(tr, s, p); +{ + nodeptr p, s; + int i; + + i = MIN(ip, iq); + if (ir < i) i = ir; + tr->start = tr->nodep[i]; + tr->ntips = 3; + p = tr->nodep[ip]; + hookupParsimony(p, tr->nodep[iq]); + s = buildNewTip(tr, tr->nodep[ir]); + insertParsimony(tr, s, p); } static void testInsertParsimony (tree *tr, nodeptr p, nodeptr q) -{ - unsigned int mp; - - nodeptr r = q->back; - - insertParsimony(tr, p, q); - - mp = evaluateParsimony(tr, p->next->next, FALSE); - - if(mp < tr->bestParsimony) - { - tr->bestParsimony = mp; - tr->insertNode = q; - tr->removeNode = p; - } - - hookupParsimony(q, r); - p->next->next->back = p->next->back = (nodeptr) NULL; - - return; -} +{ + unsigned int mp; + + nodeptr r = q->back; + + insertParsimony(tr, p, q); + + mp = evaluateParsimony(tr, p->next->next, FALSE); + + if(mp < tr->bestParsimony) + { + tr->bestParsimony = mp; + tr->insertNode = q; + tr->removeNode = p; + } + + hookupParsimony(q, r); + p->next->next->back = p->next->back = (nodeptr) NULL; + + return; +} static void restoreTreeParsimony(tree *tr, nodeptr p, nodeptr q) -{ - nodeptr - r = q->back; - - int counter = 4; - - hookupParsimony(p->next, q); - hookupParsimony(p->next->next, r); - - computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE); - tr->ti[0] = counter; - - newviewParsimonyIterativeFast(tr); +{ + nodeptr + r = q->back; + + int counter = 4; + + hookupParsimony(p->next, q); + hookupParsimony(p->next->next, r); + + computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE); + tr->ti[0] = counter; + + newviewParsimonyIterativeFast(tr); } static void addTraverseParsimony (tree *tr, nodeptr p, nodeptr q, int mintrav, int maxtrav, boolean doAll) -{ - if (doAll || (--mintrav <= 0)) - testInsertParsimony(tr, p, q); - - if (((q->number > tr->mxtips)) && ((--maxtrav > 0) || doAll)) - { - addTraverseParsimony(tr, p, q->next->back, mintrav, maxtrav, doAll); - addTraverseParsimony(tr, p, q->next->next->back, mintrav, maxtrav, doAll); - } +{ + if (doAll || (--mintrav <= 0)) + testInsertParsimony(tr, p, q); + + if (((q->number > tr->mxtips)) && ((--maxtrav > 0) || doAll)) + { + addTraverseParsimony(tr, p, q->next->back, mintrav, maxtrav, doAll); + addTraverseParsimony(tr, p, q->next->next->back, mintrav, maxtrav, doAll); + } } static nodeptr findAnyTipFast(nodeptr p, int numsp) -{ - return (p->number <= numsp)? p : findAnyTipFast(p->next->back, numsp); -} +{ + return (p->number <= numsp) ? p : findAnyTipFast(p->next->back, numsp); +} static double randum(long *seed) { - long sum, mult0, mult1, seed0, seed1, seed2, newseed0, newseed1, newseed2; - double res; - - mult0 = 1549; - seed0 = *seed & 4095; - sum = mult0 * seed0; - newseed0 = sum & 4095; - sum >>= 12; - seed1 = (*seed >> 12) & 4095; - mult1 = 406; - sum += mult0 * seed1 + mult1 * seed0; - newseed1 = sum & 4095; - sum >>= 12; - seed2 = (*seed >> 24) & 255; - sum += mult0 * seed2 + mult1 * seed1; - newseed2 = sum & 255; - - *seed = newseed2 << 24 | newseed1 << 12 | newseed0; - res = 0.00390625 * (newseed2 + 0.000244140625 * (newseed1 + 0.000244140625 * newseed0)); - - return res; + long sum, mult0, mult1, seed0, seed1, seed2, newseed0, newseed1, newseed2; + double res; + + mult0 = 1549; + seed0 = *seed & 4095; + sum = mult0 * seed0; + newseed0 = sum & 4095; + sum >>= 12; + seed1 = (*seed >> 12) & 4095; + mult1 = 406; + sum += mult0 * seed1 + mult1 * seed0; + newseed1 = sum & 4095; + sum >>= 12; + seed2 = (*seed >> 24) & 255; + sum += mult0 * seed2 + mult1 * seed1; + newseed2 = sum & 255; + + *seed = newseed2 << 24 | newseed1 << 12 | newseed0; + res = 0.00390625 * (newseed2 + 0.000244140625 * (newseed1 + 0.000244140625 * newseed0)); + + return res; } static void makePermutationFast(int *perm, int n, analdef *adef) -{ - int i, j, k; +{ + int i, j, k; + + assert(adef->parsimonySeed != 0); - assert(adef->parsimonySeed != 0); - - - for (i = 1; i <= n; i++) - perm[i] = i; - for (i = 1; i <= n; i++) - { - double d = randum(&adef->parsimonySeed); + for (i = 1; i <= n; i++) + perm[i] = i; - k = (int)((double)(n + 1 - i) * d); - - j = perm[i]; + for (i = 1; i <= n; i++) + { + double d = randum(&adef->parsimonySeed); - perm[i] = perm[i + k]; - perm[i + k] = j; - } + k = (int)((double)(n + 1 - i) * d); + + j = perm[i]; + + perm[i] = perm[i + k]; + perm[i + k] = j; + } } static nodeptr removeNodeParsimony (nodeptr p) -{ - nodeptr q, r; +{ + nodeptr q, r; - q = p->next->back; - r = p->next->next->back; - - hookupParsimony(q, r); + q = p->next->back; + r = p->next->next->back; - p->next->next->back = p->next->back = (node *) NULL; - - return q; + hookupParsimony(q, r); + + p->next->next->back = p->next->back = (node *) NULL; + + return q; } -static int rearrangeParsimony(tree *tr, nodeptr p, int mintrav, int maxtrav, boolean doAll) -{ - nodeptr p1, p2, q, q1, q2; - int mintrav2; - - if (maxtrav > tr->ntips - 3) - maxtrav = tr->ntips - 3; - - - if(maxtrav < mintrav) - return 0; - - q = p->back; - - if(p->number > tr->mxtips) - { - p1 = p->next->back; - p2 = p->next->next->back; - - if ((p1->number > tr->mxtips) || (p2->number > tr->mxtips)) - { - removeNodeParsimony(p); - - if ((p1->number > tr->mxtips)) - { - addTraverseParsimony(tr, p, p1->next->back, mintrav, maxtrav, doAll); - addTraverseParsimony(tr, p, p1->next->next->back, mintrav, maxtrav, doAll); - } - - if ((p2->number > tr->mxtips)) - { - addTraverseParsimony(tr, p, p2->next->back, mintrav, maxtrav, doAll); - addTraverseParsimony(tr, p, p2->next->next->back, mintrav, maxtrav, doAll); - } - - - hookupParsimony(p->next, p1); - hookupParsimony(p->next->next, p2); - - newviewParsimony(tr, p); - } - } - - if ((q->number > tr->mxtips) && maxtrav > 0) - { - q1 = q->next->back; - q2 = q->next->next->back; - - if ( - ( - (q1->number > tr->mxtips) && - ((q1->next->back->number > tr->mxtips) || (q1->next->next->back->number > tr->mxtips)) - ) - || - ( - (q2->number > tr->mxtips) && - ((q2->next->back->number > tr->mxtips) || (q2->next->next->back->number > tr->mxtips)) - ) - ) - { - - removeNodeParsimony(q); - - mintrav2 = mintrav > 2 ? mintrav : 2; - - if ((q1->number > tr->mxtips)) - { - addTraverseParsimony(tr, q, q1->next->back, mintrav2 , maxtrav, doAll); - addTraverseParsimony(tr, q, q1->next->next->back, mintrav2 , maxtrav, doAll); - } - - if ((q2->number > tr->mxtips)) - { - addTraverseParsimony(tr, q, q2->next->back, mintrav2 , maxtrav, doAll); - addTraverseParsimony(tr, q, q2->next->next->back, mintrav2 , maxtrav, doAll); - } - - hookupParsimony(q->next, q1); - hookupParsimony(q->next->next, q2); - - newviewParsimony(tr, q); - } - } - - return 1; -} +static int rearrangeParsimony(tree *tr, nodeptr p, int mintrav, int maxtrav, boolean doAll) +{ + nodeptr p1, p2, q, q1, q2; + int mintrav2; + + if (maxtrav > tr->ntips - 3) + maxtrav = tr->ntips - 3; + + + if(maxtrav < mintrav) + return 0; + + q = p->back; + + if(p->number > tr->mxtips) + { + p1 = p->next->back; + p2 = p->next->next->back; + + if ((p1->number > tr->mxtips) || (p2->number > tr->mxtips)) + { + removeNodeParsimony(p); + + if ((p1->number > tr->mxtips)) + { + addTraverseParsimony(tr, p, p1->next->back, mintrav, maxtrav, doAll); + addTraverseParsimony(tr, p, p1->next->next->back, mintrav, maxtrav, doAll); + } + + if ((p2->number > tr->mxtips)) + { + addTraverseParsimony(tr, p, p2->next->back, mintrav, maxtrav, doAll); + addTraverseParsimony(tr, p, p2->next->next->back, mintrav, maxtrav, doAll); + } + + + hookupParsimony(p->next, p1); + hookupParsimony(p->next->next, p2); + + newviewParsimony(tr, p); + } + } + + if ((q->number > tr->mxtips) && maxtrav > 0) + { + q1 = q->next->back; + q2 = q->next->next->back; + + if ( + ( + (q1->number > tr->mxtips) && + ((q1->next->back->number > tr->mxtips) || (q1->next->next->back->number > tr->mxtips)) + ) + || + ( + (q2->number > tr->mxtips) && + ((q2->next->back->number > tr->mxtips) || (q2->next->next->back->number > tr->mxtips)) + ) + ) + { + + removeNodeParsimony(q); + + mintrav2 = mintrav > 2 ? mintrav : 2; + + if ((q1->number > tr->mxtips)) + { + addTraverseParsimony(tr, q, q1->next->back, mintrav2, maxtrav, doAll); + addTraverseParsimony(tr, q, q1->next->next->back, mintrav2, maxtrav, doAll); + } + + if ((q2->number > tr->mxtips)) + { + addTraverseParsimony(tr, q, q2->next->back, mintrav2, maxtrav, doAll); + addTraverseParsimony(tr, q, q2->next->next->back, mintrav2, maxtrav, doAll); + } + + hookupParsimony(q->next, q1); + hookupParsimony(q->next->next, q2); + + newviewParsimony(tr, q); + } + } + + return 1; +} static void restoreTreeRearrangeParsimony(tree *tr) -{ - removeNodeParsimony(tr->removeNode); - restoreTreeParsimony(tr, tr->removeNode, tr->insertNode); +{ + removeNodeParsimony(tr->removeNode); + restoreTreeParsimony(tr, tr->removeNode, tr->insertNode); } static boolean isInformative(tree *tr, int site) { - int - informativeCounter = 0, - check[256], - j, - undetermined = 15; - - unsigned char - nucleotide, - target = 0; - - for(j = 0; j < 256; j++) - check[j] = 0; - - for(j = 1; j <= tr->mxtips; j++) - { - nucleotide = tr->yVector[j][site]; - check[nucleotide] = check[nucleotide] + 1; - } - - - if(check[1] > 1) - { - informativeCounter++; - target = target | 1; - } - if(check[2] > 1) - { - informativeCounter++; - target = target | 2; - } - if(check[4] > 1) - { - informativeCounter++; - target = target | 4; - } - if(check[8] > 1) - { - informativeCounter++; - target = target | 8; - } - - if(informativeCounter >= 2) - return TRUE; - else - { - for(j = 0; j < undetermined; j++) - { - if(j == 3 || j == 5 || j == 6 || j == 7 || j == 9 || j == 10 || j == 11 || - j == 12 || j == 13 || j == 14) - { - if(check[j] > 1) - { - if(!(target & j)) - return TRUE; - } - } - } - } - - return FALSE; + int + informativeCounter = 0, + check[256], + j, + undetermined = 15; + + unsigned char + nucleotide, + target = 0; + + for(j = 0; j < 256; j++) + check[j] = 0; + + for(j = 1; j <= tr->mxtips; j++) + { + nucleotide = tr->yVector[j][site]; + check[nucleotide] = check[nucleotide] + 1; + } + + + if(check[1] > 1) + { + informativeCounter++; + target = target | 1; + } + if(check[2] > 1) + { + informativeCounter++; + target = target | 2; + } + if(check[4] > 1) + { + informativeCounter++; + target = target | 4; + } + if(check[8] > 1) + { + informativeCounter++; + target = target | 8; + } + + if(informativeCounter >= 2) + return TRUE; + else + { + for(j = 0; j < undetermined; j++) + { + if(j == 3 || j == 5 || j == 6 || j == 7 || j == 9 || j == 10 || j == 11 || + j == 12 || j == 13 || j == 14) + { + if(check[j] > 1) + { + if(!(target & j)) + return TRUE; + } + } + } + } + + return FALSE; } static void determineUninformativeSites(tree *tr, int *informative) { - int - i; - //number = 0; - - /* - Not all characters are useful in constructing a parsimony tree. - Invariant characters, those that have the same state in all taxa, - are obviously useless and are ignored by the method. Characters in - which a state occurs in only one taxon are also ignored. - All these characters are called parsimony uninformative. - - Alternative definition: informative columns contain at least two types - of nucleotides, and each nucleotide must appear at least twice in each - column. Kind of a pain if we intend to check for this when using, e.g., - amibiguous DNA encoding. - */ - - for(i = 0; i < tr->rdta->sites; i++) - { - if(isInformative(tr, i)) - informative[i] = 1; - else - { - informative[i] = 0; - //number++; - } - } - - - /* printf("Uninformative Patterns: %d\n", number); */ + int + i; + //number = 0; + + /* + Not all characters are useful in constructing a parsimony tree. + Invariant characters, those that have the same state in all taxa, + are obviously useless and are ignored by the method. Characters in + which a state occurs in only one taxon are also ignored. + All these characters are called parsimony uninformative. + + Alternative definition: informative columns contain at least two types + of nucleotides, and each nucleotide must appear at least twice in each + column. Kind of a pain if we intend to check for this when using, e.g., + amibiguous DNA encoding. + */ + + for(i = 0; i < tr->rdta->sites; i++) + { + if(isInformative(tr, i)) + informative[i] = 1; + else + { + informative[i] = 0; + //number++; + } + } + + + /* printf("Uninformative Patterns: %d\n", number); */ } static void reorderNodes(tree *tr, nodeptr *np, nodeptr p, int *count) { - int i, found = 0; - - if((p->number <= tr->mxtips)) - return; - else - { - for(i = tr->mxtips + 1; (i <= (tr->mxtips + tr->mxtips - 1)) && (found == 0); i++) - { - if (p == np[i] || p == np[i]->next || p == np[i]->next->next) - { - if(p == np[i]) - tr->nodep[*count + tr->mxtips + 1] = np[i]; - else - { - if(p == np[i]->next) - tr->nodep[*count + tr->mxtips + 1] = np[i]->next; - else - tr->nodep[*count + tr->mxtips + 1] = np[i]->next->next; - } - - found = 1; - *count = *count + 1; - } - } - - reorderNodes(tr, np, p->next->back, count); - reorderNodes(tr, np, p->next->next->back, count); - } + int i, found = 0; + + if((p->number <= tr->mxtips)) + return; + else + { + for(i = tr->mxtips + 1; (i <= (tr->mxtips + tr->mxtips - 1)) && (found == 0); i++) + { + if (p == np[i] || p == np[i]->next || p == np[i]->next->next) + { + if(p == np[i]) + tr->nodep[*count + tr->mxtips + 1] = np[i]; + else + { + if(p == np[i]->next) + tr->nodep[*count + tr->mxtips + 1] = np[i]->next; + else + tr->nodep[*count + tr->mxtips + 1] = np[i]->next->next; + } + + found = 1; + *count = *count + 1; + } + } + + reorderNodes(tr, np, p->next->back, count); + reorderNodes(tr, np, p->next->next->back, count); + } } static void nodeRectifierFast(tree *tr) { - nodeptr *np = (nodeptr *)malloc(2 * (size_t)tr->mxtips * sizeof(nodeptr)); - int i; - int count = 0; - - tr->start = tr->nodep[1]; - - for(i = tr->mxtips + 1; i <= (tr->mxtips + tr->mxtips - 1); i++) - np[i] = tr->nodep[i]; - - reorderNodes(tr, np, tr->start->back, &count); - - free(np); + nodeptr *np = (nodeptr *)malloc(2 * (size_t)tr->mxtips * sizeof(nodeptr)); + int i; + int count = 0; + + tr->start = tr->nodep[1]; + + for(i = tr->mxtips + 1; i <= (tr->mxtips + tr->mxtips - 1); i++) + np[i] = tr->nodep[i]; + + reorderNodes(tr, np, tr->start->back, &count); + + free(np); } - + static parsimonyNumber *compressDNA(tree *tr, int *informative) { - size_t - i, - entries = 0, - compressedEntries, - compressedEntriesPadded; - - parsimonyNumber - *compressedScratch; - - for(i = 0; i < (size_t)tr->rdta->sites; i++) - if(informative[i]) - entries += 1; - - compressedEntries = entries / PCF; - - if(entries % PCF != 0) - compressedEntries++; - - /* printf("compression %d -> %d\n", entries, compressedEntries); */ + size_t + i, + entries = 0, + compressedEntries, + compressedEntriesPadded; + + parsimonyNumber + *compressedScratch; + + for(i = 0; i < (size_t)tr->rdta->sites; i++) + if(informative[i]) + entries += 1; + + compressedEntries = entries / PCF; + + if(entries % PCF != 0) + compressedEntries++; + + /* printf("compression %d -> %d\n", entries, compressedEntries); */ #if (defined(__SIM_SSE3) || defined(__AVX)) - if(compressedEntries % INTS_PER_VECTOR != 0) - compressedEntriesPadded = compressedEntries + (INTS_PER_VECTOR - (compressedEntries % INTS_PER_VECTOR)); - else - compressedEntriesPadded = compressedEntries; + if(compressedEntries % INTS_PER_VECTOR != 0) + compressedEntriesPadded = compressedEntries + (INTS_PER_VECTOR - (compressedEntries % INTS_PER_VECTOR)); + else + compressedEntriesPadded = compressedEntries; #else - compressedEntriesPadded = compressedEntries; + compressedEntriesPadded = compressedEntries; #endif - /* printf("padded %d\n", compressedEntriesPadded); */ - - compressedScratch = (parsimonyNumber *)malloc_aligned((size_t)compressedEntriesPadded * 8 * (size_t)tr->mxtips * sizeof(parsimonyNumber), BYTE_ALIGNMENT); - - - for(i = 0; i < compressedEntriesPadded * 8 * tr->mxtips; i++) - compressedScratch[i] = 0; - - - - for(i = 0; i < (size_t)tr->mxtips; i++) - { - parsimonyNumber - *compressedTip_A = &compressedScratch[(compressedEntriesPadded * 4 * (i + 1))], - *compressedTip_C = &compressedScratch[(compressedEntriesPadded * 4 * (i + 1)) + compressedEntriesPadded * 1], - *compressedTip_G = &compressedScratch[(compressedEntriesPadded * 4 * (i + 1)) + compressedEntriesPadded * 2], - *compressedTip_T = &compressedScratch[(compressedEntriesPadded * 4 * (i + 1)) + compressedEntriesPadded * 3], - compressedValue_A = 0, - compressedValue_C = 0, - compressedValue_G = 0, - compressedValue_T = 0; - - size_t - w = 0, - compressedIndex = 0, - compressedCounter = 0, - index = 0; - - for(index = 0; index < (size_t)tr->rdta->sites; index++) - { - if(informative[index]) - { - parsimonyNumber value = (parsimonyNumber)(tr->yVector[i + 1][index]); - - for(w = 0; w < 1; w++) - { - if(value & 1) - compressedValue_A |= mask32[compressedCounter]; - if(value & 2) - compressedValue_C |= mask32[compressedCounter]; - if(value & 4) - compressedValue_G |= mask32[compressedCounter]; - if(value & 8) - compressedValue_T |= mask32[compressedCounter]; - - compressedCounter++; - - if(compressedCounter == PCF) - { - compressedTip_A[compressedIndex] = compressedValue_A; - compressedTip_C[compressedIndex] = compressedValue_C; - compressedTip_G[compressedIndex] = compressedValue_G; - compressedTip_T[compressedIndex] = compressedValue_T; - - compressedValue_A = 0; - compressedValue_C = 0; - compressedValue_G = 0; - compressedValue_T = 0; - - compressedCounter = 0; - compressedIndex++; - } - } - } - } - - for(;compressedIndex < compressedEntriesPadded; compressedIndex++) - { - for(;compressedCounter < PCF; compressedCounter++) - { - compressedValue_A |= mask32[compressedCounter]; - compressedValue_C |= mask32[compressedCounter]; - compressedValue_G |= mask32[compressedCounter]; - compressedValue_T |= mask32[compressedCounter]; - } - - compressedTip_A[compressedIndex] = compressedValue_A; - compressedTip_C[compressedIndex] = compressedValue_C; - compressedTip_G[compressedIndex] = compressedValue_G; - compressedTip_T[compressedIndex] = compressedValue_T; - - compressedValue_A = 0; - compressedValue_C = 0; - compressedValue_G = 0; - compressedValue_T = 0; - - compressedCounter = 0; - } - } - - - tr->parsimonyState_A = (parsimonyNumber**)malloc_aligned(sizeof(parsimonyNumber*) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); - tr->parsimonyState_C = (parsimonyNumber**)malloc_aligned(sizeof(parsimonyNumber*) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); - tr->parsimonyState_G = (parsimonyNumber**)malloc_aligned(sizeof(parsimonyNumber*) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); - tr->parsimonyState_T = (parsimonyNumber**)malloc_aligned(sizeof(parsimonyNumber*) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); - - tr->parsimonyScore = (unsigned int*)malloc_aligned(sizeof(unsigned int) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); - - for(i = 0; i < 2 * (size_t)tr->mxtips; i++) - { - tr->parsimonyState_A[i] = &compressedScratch[i * 4 * compressedEntriesPadded]; - tr->parsimonyState_C[i] = &compressedScratch[i * 4 * compressedEntriesPadded + compressedEntriesPadded * 1]; - tr->parsimonyState_G[i] = &compressedScratch[i * 4 * compressedEntriesPadded + compressedEntriesPadded * 2]; - tr->parsimonyState_T[i] = &compressedScratch[i * 4 * compressedEntriesPadded + compressedEntriesPadded * 3]; - - tr->parsimonyScore[i] = 0; - } - - tr->compressedWidth = (unsigned int)compressedEntriesPadded; - - return compressedScratch; + /* printf("padded %d\n", compressedEntriesPadded); */ + + compressedScratch = (parsimonyNumber *)malloc_aligned((size_t)compressedEntriesPadded * 8 * (size_t)tr->mxtips * sizeof(parsimonyNumber), BYTE_ALIGNMENT); + + + for(i = 0; i < compressedEntriesPadded * 8 * tr->mxtips; i++) + compressedScratch[i] = 0; + + + + for(i = 0; i < (size_t)tr->mxtips; i++) + { + parsimonyNumber + *compressedTip_A = &compressedScratch[(compressedEntriesPadded * 4 * (i + 1))], + *compressedTip_C = &compressedScratch[(compressedEntriesPadded * 4 * (i + 1)) + compressedEntriesPadded * 1], + *compressedTip_G = &compressedScratch[(compressedEntriesPadded * 4 * (i + 1)) + compressedEntriesPadded * 2], + *compressedTip_T = &compressedScratch[(compressedEntriesPadded * 4 * (i + 1)) + compressedEntriesPadded * 3], + compressedValue_A = 0, + compressedValue_C = 0, + compressedValue_G = 0, + compressedValue_T = 0; + + size_t + w = 0, + compressedIndex = 0, + compressedCounter = 0, + index = 0; + + for(index = 0; index < (size_t)tr->rdta->sites; index++) + { + if(informative[index]) + { + parsimonyNumber value = (parsimonyNumber)(tr->yVector[i + 1][index]); + + for(w = 0; w < 1; w++) + { + if(value & 1) + compressedValue_A |= mask32[compressedCounter]; + if(value & 2) + compressedValue_C |= mask32[compressedCounter]; + if(value & 4) + compressedValue_G |= mask32[compressedCounter]; + if(value & 8) + compressedValue_T |= mask32[compressedCounter]; + + compressedCounter++; + + if(compressedCounter == PCF) + { + compressedTip_A[compressedIndex] = compressedValue_A; + compressedTip_C[compressedIndex] = compressedValue_C; + compressedTip_G[compressedIndex] = compressedValue_G; + compressedTip_T[compressedIndex] = compressedValue_T; + + compressedValue_A = 0; + compressedValue_C = 0; + compressedValue_G = 0; + compressedValue_T = 0; + + compressedCounter = 0; + compressedIndex++; + } + } + } + } + + for(; compressedIndex < compressedEntriesPadded; compressedIndex++) + { + for(; compressedCounter < PCF; compressedCounter++) + { + compressedValue_A |= mask32[compressedCounter]; + compressedValue_C |= mask32[compressedCounter]; + compressedValue_G |= mask32[compressedCounter]; + compressedValue_T |= mask32[compressedCounter]; + } + + compressedTip_A[compressedIndex] = compressedValue_A; + compressedTip_C[compressedIndex] = compressedValue_C; + compressedTip_G[compressedIndex] = compressedValue_G; + compressedTip_T[compressedIndex] = compressedValue_T; + + compressedValue_A = 0; + compressedValue_C = 0; + compressedValue_G = 0; + compressedValue_T = 0; + + compressedCounter = 0; + } + } + + + tr->parsimonyState_A = (parsimonyNumber**)malloc_aligned(sizeof(parsimonyNumber*) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); + tr->parsimonyState_C = (parsimonyNumber**)malloc_aligned(sizeof(parsimonyNumber*) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); + tr->parsimonyState_G = (parsimonyNumber**)malloc_aligned(sizeof(parsimonyNumber*) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); + tr->parsimonyState_T = (parsimonyNumber**)malloc_aligned(sizeof(parsimonyNumber*) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); + + tr->parsimonyScore = (unsigned int*)malloc_aligned(sizeof(unsigned int) * 2 * (size_t)tr->mxtips, BYTE_ALIGNMENT); + + for(i = 0; i < 2 * (size_t)tr->mxtips; i++) + { + tr->parsimonyState_A[i] = &compressedScratch[i * 4 * compressedEntriesPadded]; + tr->parsimonyState_C[i] = &compressedScratch[i * 4 * compressedEntriesPadded + compressedEntriesPadded * 1]; + tr->parsimonyState_G[i] = &compressedScratch[i * 4 * compressedEntriesPadded + compressedEntriesPadded * 2]; + tr->parsimonyState_T[i] = &compressedScratch[i * 4 * compressedEntriesPadded + compressedEntriesPadded * 3]; + + tr->parsimonyScore[i] = 0; + } + + tr->compressedWidth = (unsigned int)compressedEntriesPadded; + + return compressedScratch; } static void stepwiseAddition(tree *tr, nodeptr p, nodeptr q) -{ - nodeptr - r = q->back; - - unsigned int - mp, - bestParsimony = tr->bestParsimony; - - int - counter = 4; - - p->next->back = q; - q->back = p->next; - - p->next->next->back = r; - r->back = p->next->next; - - computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE); - tr->ti[0] = counter; - - newviewParsimonyIterativeFast(tr); - +{ + nodeptr + r = q->back; + + unsigned int + mp, + bestParsimony = tr->bestParsimony; + + int + counter = 4; + + p->next->back = q; + q->back = p->next; + + p->next->next->back = r; + r->back = p->next->next; + + computeTraversalInfoParsimony(p, tr->ti, &counter, tr->mxtips, FALSE); + tr->ti[0] = counter; + + newviewParsimonyIterativeFast(tr); + #if (defined(__SIM_SSE3) || defined(__AVX)) - { - INT_TYPE - allOne = SET_ALL_BITS_ONE; - - int - pNumber = p->number, - qNumber = p->back->number; - - unsigned int - width = tr->compressedWidth, - i; - - parsimonyNumber - *rightState_A = tr->parsimonyState_A[pNumber], - *leftState_A = tr->parsimonyState_A[qNumber], - - *rightState_C = tr->parsimonyState_C[pNumber], - *leftState_C = tr->parsimonyState_C[qNumber], - - *rightState_G = tr->parsimonyState_G[pNumber], - *leftState_G = tr->parsimonyState_G[qNumber], - - *rightState_T = tr->parsimonyState_T[pNumber], - *leftState_T = tr->parsimonyState_T[qNumber]; - - mp = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber]; + { + INT_TYPE + allOne = SET_ALL_BITS_ONE; + + int + pNumber = p->number, + qNumber = p->back->number; + + unsigned int + width = tr->compressedWidth, + i; + + parsimonyNumber + *rightState_A = tr->parsimonyState_A[pNumber], + *leftState_A = tr->parsimonyState_A[qNumber], + + *rightState_C = tr->parsimonyState_C[pNumber], + *leftState_C = tr->parsimonyState_C[qNumber], + + *rightState_G = tr->parsimonyState_G[pNumber], + *leftState_G = tr->parsimonyState_G[qNumber], + + *rightState_T = tr->parsimonyState_T[pNumber], + *leftState_T = tr->parsimonyState_T[qNumber]; + + mp = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber]; #pragma omp parallel for reduction(+:mp) - for(i = 0; i < width; i += INTS_PER_VECTOR) - { - INT_TYPE - l_A = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_A[i])), VECTOR_LOAD((CAST)(&rightState_A[i]))), - l_C = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_C[i])), VECTOR_LOAD((CAST)(&rightState_C[i]))), - l_G = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_G[i])), VECTOR_LOAD((CAST)(&rightState_G[i]))), - l_T = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_T[i])), VECTOR_LOAD((CAST)(&rightState_T[i]))), - v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T)); + for(i = 0; i < width; i += INTS_PER_VECTOR) + { + INT_TYPE + l_A = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_A[i])), VECTOR_LOAD((CAST)(&rightState_A[i]))), + l_C = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_C[i])), VECTOR_LOAD((CAST)(&rightState_C[i]))), + l_G = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_G[i])), VECTOR_LOAD((CAST)(&rightState_G[i]))), + l_T = VECTOR_BIT_AND(VECTOR_LOAD((CAST)(&leftState_T[i])), VECTOR_LOAD((CAST)(&rightState_T[i]))), + v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T)); - v_N = VECTOR_AND_NOT(v_N, allOne); + v_N = VECTOR_AND_NOT(v_N, allOne); - mp += populationCount(v_N); + mp += populationCount(v_N); #ifndef _OPENMP - if(mp >= bestParsimony) - goto SKIP; + if(mp >= bestParsimony) + goto SKIP; #endif - } - } + } + } #else - { - int - pNumber = p->number, - qNumber = p->back->number; - - unsigned int - width = tr->compressedWidth, - i, - t_A, - t_C, - t_G, - t_T, - t_N; - - parsimonyNumber - *rightState_A = tr->parsimonyState_A[pNumber], - *leftState_A = tr->parsimonyState_A[qNumber], - - *rightState_C = tr->parsimonyState_C[pNumber], - *leftState_C = tr->parsimonyState_C[qNumber], - - *rightState_G = tr->parsimonyState_G[pNumber], - *leftState_G = tr->parsimonyState_G[qNumber], - - *rightState_T = tr->parsimonyState_T[pNumber], - *leftState_T = tr->parsimonyState_T[qNumber]; - - mp = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber]; - - for(i = 0; i < width; i++) - { - t_A = leftState_A[i] & rightState_A[i]; - t_C = leftState_C[i] & rightState_C[i]; - t_G = leftState_G[i] & rightState_G[i]; - t_T = leftState_T[i] & rightState_T[i]; - - t_N = ~(t_A | t_C | t_G | t_T); - - mp += populationCount(t_N); - - if(mp >= bestParsimony) - goto SKIP; - } - } + { + int + pNumber = p->number, + qNumber = p->back->number; + + unsigned int + width = tr->compressedWidth, + i, + t_A, + t_C, + t_G, + t_T, + t_N; + + parsimonyNumber + *rightState_A = tr->parsimonyState_A[pNumber], + *leftState_A = tr->parsimonyState_A[qNumber], + + *rightState_C = tr->parsimonyState_C[pNumber], + *leftState_C = tr->parsimonyState_C[qNumber], + + *rightState_G = tr->parsimonyState_G[pNumber], + *leftState_G = tr->parsimonyState_G[qNumber], + + *rightState_T = tr->parsimonyState_T[pNumber], + *leftState_T = tr->parsimonyState_T[qNumber]; + + mp = tr->parsimonyScore[pNumber] + tr->parsimonyScore[qNumber]; + + for(i = 0; i < width; i++) + { + t_A = leftState_A[i] & rightState_A[i]; + t_C = leftState_C[i] & rightState_C[i]; + t_G = leftState_G[i] & rightState_G[i]; + t_T = leftState_T[i] & rightState_T[i]; + + t_N = ~(t_A | t_C | t_G | t_T); + + mp += populationCount(t_N); + + if(mp >= bestParsimony) + goto SKIP; + } + } #endif #ifdef _OPENMP - if(mp < tr->bestParsimony) - { -#endif - tr->bestParsimony = mp; - tr->insertNode = q; + if(mp < tr->bestParsimony) + { +#endif + tr->bestParsimony = mp; + tr->insertNode = q; #ifdef _OPENMP - } +} #endif - - SKIP: - q->back = r; - r->back = q; - - if(q->number > tr->mxtips && tr->parsimonyScore[q->number] > 0) - { - stepwiseAddition(tr, p, q->next->back); - stepwiseAddition(tr, p, q->next->next->back); - } + +SKIP: + q->back = r; + r->back = q; + + if(q->number > tr->mxtips && tr->parsimonyScore[q->number] > 0) + { + stepwiseAddition(tr, p, q->next->back); + stepwiseAddition(tr, p, q->next->next->back); + } } static char *Tree2StringREC(char *treestr, tree *tr, nodeptr p) -{ - if(p->number > 0 && p->number <= tr->mxtips) - { - sprintf(treestr, "%s", tr->nameList[p->number]); - while (*treestr) treestr++; - } - else - { - *treestr++ = '('; - treestr = Tree2StringREC(treestr, tr, p->next->back); - *treestr++ = ','; - treestr = Tree2StringREC(treestr, tr, p->next->next->back); - if(p == tr->start->back) - { - *treestr++ = ','; - treestr = Tree2StringREC(treestr, tr, p->back); - } - *treestr++ = ')'; - } - - if(p == tr->start->back) - sprintf(treestr, ";\n"); - else - sprintf(treestr, "%s", "\0"); - - while (*treestr) treestr++; - return treestr; +{ + if(p->number > 0 && p->number <= tr->mxtips) + { + sprintf(treestr, "%s", tr->nameList[p->number]); + while (*treestr) treestr++; + } + else + { + *treestr++ = '('; + treestr = Tree2StringREC(treestr, tr, p->next->back); + *treestr++ = ','; + treestr = Tree2StringREC(treestr, tr, p->next->next->back); + if(p == tr->start->back) + { + *treestr++ = ','; + treestr = Tree2StringREC(treestr, tr, p->back); + } + *treestr++ = ')'; + } + + if(p == tr->start->back) + sprintf(treestr, ";\n"); + else + sprintf(treestr, "%s", "\0"); + + while (*treestr) treestr++; + return treestr; } @@ -1308,209 +1308,205 @@ static char *Tree2StringREC(char *treestr, tree *tr, nodeptr p) int makeParsimonyTreeFastDNA(tree *tr, analdef *adef, int startSite, int endSite) -{ - nodeptr - p, - f; - - int - treeCounter = 0, - i, - nextsp, - *perm = (int *)malloc((size_t)(tr->mxtips + 1) * sizeof(int)), - *informative = (int *)malloc(sizeof(int) * (size_t)tr->rdta->sites); - - unsigned int - randomMP, - startMP; - - parsimonyNumber - *parsimonySpace; - - determineUninformativeSites(tr, informative); - - for(int i = 0; i <= startSite-1; i++){ - informative[i]=0; - } - - for(int i = endSite+1; i < tr->rdta->sites; i++){ - informative[i]=0; - } - - parsimonySpace = compressDNA(tr, informative); //TODO: Could we avoid to create compressDNA from scratch each time - - - free(informative); - - tr->ti = (int*)malloc(sizeof(int) * 4 * (size_t)tr->mxtips); - - if(adef->restart) - tr->nodesInTree = (int*)calloc((size_t)(tr->mxtips + 1), sizeof(int)); - - - int minParsimonyScoreOverRuns = INT_MAX; - - - while(treeCounter < adef->numberOfTrees) - { - double t = gettime(); - - if(adef->restart) - { - int j = 0; - - FILE *treeFile = fopen(tree_file, "rb"); - - treeReadLen(treeFile, tr); - - //we do not want to print for now - //printBothOpen("Read Starting tree with %d tips, total %d\n", tr->ntips, tr->mxtips); - - fclose(treeFile); - - tr->start = findAnyTipFast(tr->start, tr->rdta->numsp); - - tr->bestParsimony = INT_MAX; - - evaluateParsimony(tr, tr->start->back, TRUE); - - - assert(tr->start); - assert(adef->parsimonySeed != 0); - - j = tr->ntips + 1; - - for(i = 1; i <= tr->mxtips; i++) - if(tr->nodesInTree[i] == 0) - perm[j++] = i; - - for(i = tr->ntips + 1; i <= tr->mxtips; i++) - { - int k, j; - - k = (int)((double)(tr->mxtips + 1 - i) * randum(&adef->parsimonySeed)); - - assert(i + k <= tr->mxtips); - j = perm[i]; - perm[i] = perm[i + k]; - perm[i + k] = j; - } - - f = tr->start; - } - else - { - makePermutationFast(perm, tr->mxtips, adef); - - tr->ntips = 0; - - tr->nextnode = tr->mxtips + 1; - - buildSimpleTree(tr, perm[1], perm[2], perm[3]); - - f = tr->start; - } - - while(tr->ntips < tr->mxtips) - { - nodeptr q; - - tr->bestParsimony = INT_MAX; - nextsp = ++(tr->ntips); - p = tr->nodep[perm[nextsp]]; - q = tr->nodep[(tr->nextnode)++]; - p->back = q; - q->back = p; - - - stepwiseAddition(tr, q, f->back); - - { - nodeptr - r = tr->insertNode->back; - - int counter = 4; - - hookupParsimony(q->next, tr->insertNode); - hookupParsimony(q->next->next, r); - - computeTraversalInfoParsimony(q, tr->ti, &counter, tr->mxtips, FALSE); - tr->ti[0] = counter; - - newviewParsimonyIterativeFast(tr); - } - } - - - /*printf("ADD: %d\n", tr->bestParsimony); */ - - - nodeRectifierFast(tr); - - randomMP = tr->bestParsimony; - - do - { - startMP = randomMP; - nodeRectifierFast(tr); - for(i = 1; i <= tr->mxtips + tr->mxtips - 2; i++) - { - rearrangeParsimony(tr, tr->nodep[i], 1, 20, FALSE); - if(tr->bestParsimony < randomMP) - { - restoreTreeRearrangeParsimony(tr); - randomMP = tr->bestParsimony; - } - } - } - while(randomMP < startMP); - - - /* printf("OPT: %d\n", tr->bestParsimony); */ - - { - FILE - *f; - - //we do not want to print for now -/* char - buf[64], - fileName[2048]; - - sprintf(buf, "%d", treeCounter); - - strcpy(fileName, resultFileName); - strcat(fileName, "."); - strcat(fileName, buf); - - f = myfopen(fileName, "wb"); - - Tree2StringREC(tr->tree_string, tr, tr->start->back); - - fprintf(f, "%s", tr->tree_string); - - fclose(f); - - printBothOpen("Parsimony tree [%d] with length %u computed in %f seconds written to file: %s\n", treeCounter, tr->bestParsimony, (gettime() - t), fileName); - */ - //int parsScore =tr->bestParsimony; - //printf("Parsimony score %u\n", tr->bestParsimony); - minParsimonyScoreOverRuns = MIN(minParsimonyScoreOverRuns,tr->bestParsimony); - //printf("Best parsimony score %u\n", minParsimonyScoreOverRuns); - } - - treeCounter++; - } - - free(perm); - free(parsimonySpace); - - - //fprintf("Parsimony score %d", tr->bestParsimony); - return minParsimonyScoreOverRuns ;//tr->bestParsimony; - -} - - - - +{ + nodeptr + p, + f; + + int + treeCounter = 0, + i, + nextsp, + *perm = (int *)malloc((size_t)(tr->mxtips + 1) * sizeof(int)), + *informative = (int *)malloc(sizeof(int) * (size_t)tr->rdta->sites); + + unsigned int + randomMP, + startMP; + + parsimonyNumber + *parsimonySpace; + + determineUninformativeSites(tr, informative); + + for(int i = 0; i <= startSite-1; i++) { + informative[i]=0; + } + + for(int i = endSite+1; i < tr->rdta->sites; i++) { + informative[i]=0; + } + + parsimonySpace = compressDNA(tr, informative); //TODO: Could we avoid to create compressDNA from scratch each time + + + free(informative); + + tr->ti = (int*)malloc(sizeof(int) * 4 * (size_t)tr->mxtips); + + if(adef->restart) + tr->nodesInTree = (int*)calloc((size_t)(tr->mxtips + 1), sizeof(int)); + + + int minParsimonyScoreOverRuns = INT_MAX; + + + while(treeCounter < adef->numberOfTrees) + { + double t = gettime(); + + if(adef->restart) + { + int j = 0; + + FILE *treeFile = fopen(tree_file, "rb"); + + treeReadLen(treeFile, tr); + + //we do not want to print for now + //printBothOpen("Read Starting tree with %d tips, total %d\n", tr->ntips, tr->mxtips); + + fclose(treeFile); + + tr->start = findAnyTipFast(tr->start, tr->rdta->numsp); + + tr->bestParsimony = INT_MAX; + + evaluateParsimony(tr, tr->start->back, TRUE); + + + assert(tr->start); + assert(adef->parsimonySeed != 0); + + j = tr->ntips + 1; + + for(i = 1; i <= tr->mxtips; i++) + if(tr->nodesInTree[i] == 0) + perm[j++] = i; + + for(i = tr->ntips + 1; i <= tr->mxtips; i++) + { + int k, j; + + k = (int)((double)(tr->mxtips + 1 - i) * randum(&adef->parsimonySeed)); + + assert(i + k <= tr->mxtips); + j = perm[i]; + perm[i] = perm[i + k]; + perm[i + k] = j; + } + + f = tr->start; + } + else + { + makePermutationFast(perm, tr->mxtips, adef); + + tr->ntips = 0; + + tr->nextnode = tr->mxtips + 1; + + buildSimpleTree(tr, perm[1], perm[2], perm[3]); + + f = tr->start; + } + + while(tr->ntips < tr->mxtips) + { + nodeptr q; + + tr->bestParsimony = INT_MAX; + nextsp = ++(tr->ntips); + p = tr->nodep[perm[nextsp]]; + q = tr->nodep[(tr->nextnode)++]; + p->back = q; + q->back = p; + + + stepwiseAddition(tr, q, f->back); + + { + nodeptr + r = tr->insertNode->back; + + int counter = 4; + + hookupParsimony(q->next, tr->insertNode); + hookupParsimony(q->next->next, r); + + computeTraversalInfoParsimony(q, tr->ti, &counter, tr->mxtips, FALSE); + tr->ti[0] = counter; + + newviewParsimonyIterativeFast(tr); + } + } + + + /*printf("ADD: %d\n", tr->bestParsimony); */ + + + nodeRectifierFast(tr); + + randomMP = tr->bestParsimony; + + do + { + startMP = randomMP; + nodeRectifierFast(tr); + for(i = 1; i <= tr->mxtips + tr->mxtips - 2; i++) + { + rearrangeParsimony(tr, tr->nodep[i], 1, 20, FALSE); + if(tr->bestParsimony < randomMP) + { + restoreTreeRearrangeParsimony(tr); + randomMP = tr->bestParsimony; + } + } + } + while(randomMP < startMP); + + + /* printf("OPT: %d\n", tr->bestParsimony); */ + + { + FILE + *f; + + //we do not want to print for now +/* char + buf[64], + fileName[2048]; + + sprintf(buf, "%d", treeCounter); + + strcpy(fileName, resultFileName); + strcat(fileName, "."); + strcat(fileName, buf); + + f = myfopen(fileName, "wb"); + + Tree2StringREC(tr->tree_string, tr, tr->start->back); + + fprintf(f, "%s", tr->tree_string); + + fclose(f); + + printBothOpen("Parsimony tree [%d] with length %u computed in %f seconds written to file: %s\n", treeCounter, tr->bestParsimony, (gettime() - t), fileName); + */ + //int parsScore =tr->bestParsimony; + //printf("Parsimony score %u\n", tr->bestParsimony); + minParsimonyScoreOverRuns = MIN(minParsimonyScoreOverRuns,tr->bestParsimony); + //printf("Best parsimony score %u\n", minParsimonyScoreOverRuns); + } + + treeCounter++; + } + + free(perm); + free(parsimonySpace); + + + //fprintf("Parsimony score %d", tr->bestParsimony); + return minParsimonyScoreOverRuns;//tr->bestParsimony; + +}