This is an automated email from the git hooks/post-receive script. tille pushed a commit to branch master in repository raxml.
commit d7986931af8a5b2ec08fc36bd9ec92a685f0adbb Author: Andreas Tille <[email protected]> Date: Thu May 19 18:08:34 2016 +0200 Imported Upstream version 8.2.8+dfsg --- .gitignore | 2 + Makefile.AVX.gcc | 2 +- README | 2 +- axml.c | 574 ++++++++++++++++++++++++++++++++++++++++++++--------- axml.h | 7 +- fastDNAparsimony.c | 64 +++--- treeIO.c | 13 +- 7 files changed, 533 insertions(+), 131 deletions(-) diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..2cd51c6 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +*.o +/raxmlHPC* diff --git a/Makefile.AVX.gcc b/Makefile.AVX.gcc index 5c5354b..e9d7dbe 100644 --- a/Makefile.AVX.gcc +++ b/Makefile.AVX.gcc @@ -5,7 +5,7 @@ CC = gcc CFLAGS = -D__SIM_SSE3 -msse3 -D_GNU_SOURCE -O2 -fomit-frame-pointer -funroll-loops -D__AVX #-Wall -Wunused-parameter -Wredundant-decls -Wreturn-type -Wswitch-default -Wunused-value -Wimplicit -Wimplicit-function-declaration -Wimplicit-int -Wimport -Wunused -Wunused-function -Wunused-label -Wno-int-to-pointer-cast -Wbad-function-cast -Wmissing-declarations -Wmissing-prototypes -Wnested-externs -Wold-style-definition -Wstrict-prototypes -Wpointer-sign -Wextra -Wredundant-decls -W [...] -LIBRARIES = -lm +LIBRARIES = -lm RM = rm -f diff --git a/README b/README index 1132b20..1befc7e 100644 --- a/README +++ b/README @@ -1,4 +1,4 @@ -Standard RAxML version 8.2.7 +Standard RAxML version 8.2.8 ============================================================================================================ diff --git a/axml.c b/axml.c index e10fc12..1927e39 100644 --- a/axml.c +++ b/axml.c @@ -29,6 +29,8 @@ */ #ifdef WIN32 +#define WIN32_LEAN_AND_MEAN // skips unwanted headers like socket etc. +#include <windows.h> #include <direct.h> #endif @@ -49,6 +51,8 @@ #include <limits.h> #include <inttypes.h> #include <getopt.h> +//#include <stdbool.h> + #if (defined(_WAYNE_MPI) || defined (_QUARTET_MPI)) #include <mpi.h> @@ -396,13 +400,17 @@ static void setRateHetAndDataIncrement(tree *tr, analdef *adef) 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; -#else +#ifdef WIN32 // WINDOWS build + FILETIME tm; + ULONGLONG t; +#if defined(NTDDI_WIN8) && NTDDI_VERSION >= NTDDI_WIN8 // >= WIN8 + GetSystemTimePreciseAsFileTime( &tm ); +#else // < WIN8 + GetSystemTimeAsFileTime( &tm ); +#endif + t = ((ULONGLONG)tm.dwHighDateTime << 32) | (ULONGLONG)tm.dwLowDateTime; + return (double)t / 10000000.0; +#else // Unixoid build struct timeval ttime; gettimeofday(&ttime , NULL); return ttime.tv_sec + ttime.tv_usec * 0.000001; @@ -1975,14 +1983,14 @@ static void getinput(analdef *adef, rawdata *rdta, cruncheddata *cdta, tree *tr) static unsigned char buildStates(int secModel, unsigned char v1, unsigned char v2) { - unsigned char new = 0; + unsigned char newChar = 0; switch(secModel) { case SECONDARY_DATA: - new = v1; - new = new << 4; - new = new | v2; + newChar = v1; + newChar = newChar << 4; + newChar = newChar | v2; break; case SECONDARY_DATA_6: { @@ -2030,38 +2038,38 @@ static unsigned char buildStates(int secModel, unsigned char v1, unsigned char v unsigned char n1 = meaningDNA[allowedStates[i][0]]; unsigned char n2 = meaningDNA[allowedStates[i][1]]; - new = n1; - new = new << 4; - new = new | n2; + newChar = n1; + newChar = newChar << 4; + newChar = newChar | n2; - intermediateBinaryStates[i] = new; + intermediateBinaryStates[i] = newChar; } - new = v1; - new = new << 4; - new = new | v2; + newChar = v1; + newChar = newChar << 4; + newChar = newChar | v2; for(i = 0; i < length; i++) { - if(new == intermediateBinaryStates[i]) + if(newChar == intermediateBinaryStates[i]) break; } if(i < length) - new = finalBinaryStates[i]; + newChar = finalBinaryStates[i]; else { - new = 0; + newChar = 0; for(i = 0; i < length; i++) { if(v1 & meaningDNA[allowedStates[i][0]]) { /*printf("Adding %c%c\n", allowedStates[i][0], allowedStates[i][1]);*/ - new |= finalBinaryStates[i]; + newChar |= finalBinaryStates[i]; } if(v2 & meaningDNA[allowedStates[i][1]]) { /*printf("Adding %c%c\n", allowedStates[i][0], allowedStates[i][1]);*/ - new |= finalBinaryStates[i]; + newChar |= finalBinaryStates[i]; } } } @@ -2112,25 +2120,25 @@ static unsigned char buildStates(int secModel, unsigned char v1, unsigned char v unsigned char n1 = meaningDNA[allowedStates[i][0]]; unsigned char n2 = meaningDNA[allowedStates[i][1]]; - new = n1; - new = new << 4; - new = new | n2; + newChar = n1; + newChar = newChar << 4; + newChar = newChar | n2; - intermediateBinaryStates[i] = new; + intermediateBinaryStates[i] = newChar; } - new = v1; - new = new << 4; - new = new | v2; + newChar = v1; + newChar = newChar << 4; + newChar = newChar | v2; for(i = 0; i < 6; i++) { /* exact match */ - if(new == intermediateBinaryStates[i]) + if(newChar == intermediateBinaryStates[i]) break; } if(i < 6) - new = finalBinaryStates[i]; + newChar = finalBinaryStates[i]; else { /* distinguish between exact mismatches and partial mismatches */ @@ -2142,20 +2150,20 @@ static unsigned char buildStates(int secModel, unsigned char v1, unsigned char v { /* printf("partial mismatch\n"); */ - new = 0; + newChar = 0; for(i = 0; i < 6; i++) { if((v1 & meaningDNA[allowedStates[i][0]]) && (v2 & meaningDNA[allowedStates[i][1]])) { /*printf("Adding %c%c\n", allowedStates[i][0], allowedStates[i][1]);*/ - new |= finalBinaryStates[i]; + newChar |= finalBinaryStates[i]; } else - new |= finalBinaryStates[6]; + newChar |= finalBinaryStates[6]; } } else - new = finalBinaryStates[6]; + newChar = finalBinaryStates[6]; } } break; @@ -2163,7 +2171,7 @@ static unsigned char buildStates(int secModel, unsigned char v1, unsigned char v assert(0); } - return new; + return newChar; } @@ -3704,6 +3712,7 @@ static void initAdef(analdef *adef) adef->setThreadAffinity = FALSE; adef->bootstopPermutations = 100; adef->fcThreshold = 99; + adef->sampleQuartetsWithoutReplacement = FALSE; } static int modelExists(char *model, analdef *adef) @@ -4773,7 +4782,7 @@ static void parseOutgroups(char *outgr, tree *tr) static void printVersionInfo(boolean terminal, FILE *infoFile) { char - text[11][1024]; + text[12][1024]; int i; @@ -4785,13 +4794,14 @@ static void printVersionInfo(boolean terminal, FILE *infoFile) sprintf(text[4], "Alexey Kozlov (HITS)\n"); sprintf(text[5], "Kassian Kobert (HITS)\n"); sprintf(text[6], "David Dao (KIT and HITS)\n"); - sprintf(text[7], "Nick Pattengale (Sandia)\n"); - sprintf(text[8], "Wayne Pfeiffer (SDSC)\n"); - sprintf(text[9], "Akifumi S. Tanabe (NRIFS)\n"); - sprintf(text[10], "Charlie Taylor (UF)\n\n"); + sprintf(text[7], "Sarah Lutteropp (KIT and HITS)\n"); + sprintf(text[8], "Nick Pattengale (Sandia)\n"); + sprintf(text[9], "Wayne Pfeiffer (SDSC)\n"); + sprintf(text[10], "Akifumi S. Tanabe (NRIFS)\n"); + sprintf(text[11], "Charlie Taylor (UF)\n\n"); - for(i = 0; i < 10; i++) + for(i = 0; i < 12; i++) { if(terminal) printf("%s", text[i]); @@ -5430,6 +5440,11 @@ static void printREADME(void) printf(" The allowed minimum number is 100!\n"); printf("\n"); printf(" DEFAULT: 100\n"); + printf("\n"); + printf(" --quartets-without-replacement specify that quartets are randomly subsampled, but without replacement.\n"); + printf("\n"); + printf(" DEFAULT: random sampling with replacements\n"); + printf("\n"); printf("\n\n\n\n"); } @@ -5579,7 +5594,7 @@ static void get_args(int argc, char *argv[], analdef *adef, tree *tr) while(1) { static struct - option long_options[16] = + option long_options[17] = { {"mesquite", no_argument, &flag, 1}, {"silent", no_argument, &flag, 1}, @@ -5596,6 +5611,7 @@ static void get_args(int argc, char *argv[], analdef *adef, tree *tr) {"HKY85", no_argument, &flag, 1}, {"set-thread-affinity", no_argument, &flag, 1}, {"bootstop-perms", required_argument, &flag, 1}, + {"quartets-without-replacement", no_argument, &flag, 1}, {0, 0, 0, 0} }; @@ -5800,6 +5816,9 @@ static void get_args(int argc, char *argv[], analdef *adef, tree *tr) adef->fcThreshold = perms - round((double)perms / 100.0); } break; + case 15: + adef->sampleQuartetsWithoutReplacement = TRUE; + break; default: if(flagCheck) { @@ -8355,6 +8374,7 @@ static void finalizeInfoFile(tree *tr, analdef *adef) default: assert(0); } + } break; case GENERIC_64: assert(0); @@ -8442,13 +8462,12 @@ static void finalizeInfoFile(tree *tr, analdef *adef) } } break; - case BINARY_DATA: - params += 1; - break; - default: - assert(0); - } - } + case BINARY_DATA: + params += 1; + break; + default: + assert(0); + } if(adef->useInvariant) params += 2; @@ -8743,7 +8762,7 @@ static void threadFixModelIndices(tree *tr, tree *localTree, int tid, int n) { if(i % (size_t)n == (size_t)tid) { - localTree->partitionData[model].wgt[localCounter] = tr->cdta->aliaswgt[globalCounter]; + localTree->partitionData[model].wgt[localCounter] = tr->cdta->aliaswgt[globalCounter]; localTree->partitionData[model].invariant[localCounter] = tr->invariant[globalCounter]; localTree->partitionData[model].rateCategory[localCounter] = tr->cdta->rateCategory[globalCounter]; @@ -11490,6 +11509,127 @@ unsigned int precomputed16_bitcount (unsigned int n) /* functions to compute likelihoods on quartets */ +/*** functions by Sarah for drawing quartets without replacement ***/ + +static uint64_t f2(int n, int a) +{ + return ((n - a) * (n - 1 - a) * (n - 2 - a ) / 6); +}; + +static uint64_t f3(int n, int b) +{ + return ((n - b) * (n - 1 - b) / 2); +}; + +static uint64_t f4(int n, int c) +{ + return (n-c); +}; + +static void preprocessQuartetPrefix(int numberOfTaxa, uint64_t *prefixSumF2, uint64_t *prefixSumF3, uint64_t *prefixSumF4) +{ + int + i, + n = numberOfTaxa; + + + + prefixSumF2[0] = 1; + prefixSumF3[0] = 1; + prefixSumF4[0] = 1; + + for (i = 1; i < n - 3; ++i) + { + prefixSumF2[i] = prefixSumF2[i - 1] + f2(n, i); + prefixSumF3[i] = prefixSumF3[i - 1] + f3(n, i+1); + prefixSumF4[i] = prefixSumF4[i - 1] + f4(n, i+2); + } +} + +static unsigned int binarySearch(uint64_t* array, uint64_t z, int n) +{ + unsigned int + first = 0, + last = n-3, + middle = (first + last) / 2, + lastSmaller = 0; + + while(first <= last) + { + if(array[middle] < z) + { + first = middle + 1; + lastSmaller = middle; + } + else + { + if (array[middle] > z) + last = middle-1; + else + { + // array[middle] == z + lastSmaller = middle; + break; + } + } + + middle = (first + last)/2; + } + + return lastSmaller; +} + +static void mapNumberToQuartet(int numberOfTaxa, uint64_t z, int *t1, int *t2, int *t3, int *t4, uint64_t *prefixSumF2, uint64_t *prefixSumF3, uint64_t *prefixSumF4) +{ + uint64_t + wantedT1 = z; + + *t1 = binarySearch(prefixSumF2, z, numberOfTaxa) + 1; + + uint64_t + foundT1 = prefixSumF2[*t1 - 1]; + + if(wantedT1 == foundT1) + { + *t2 = *t1+1; + *t3 = *t1+2; + *t4 = *t1+3; + return; + } + + uint64_t + wantedT2 = (prefixSumF3[*t1 - 1]) + (wantedT1 - foundT1); + + *t2 = binarySearch(prefixSumF3, wantedT2, numberOfTaxa) + 2; + + uint64_t + foundT2 = prefixSumF3[*t2 - 2]; + + if(wantedT2 == foundT2) + { + *t3 = *t2 + 1; + *t4 = *t2 + 2; + return; + } + + uint64_t + wantedT3 = (prefixSumF4[*t2 - 2]) + (wantedT2 - foundT2); + + *t3 = binarySearch(prefixSumF4, wantedT3, numberOfTaxa) + 3; + + uint64_t + foundT3 = prefixSumF4[*t3 - 3]; + + if (wantedT3 == foundT3) + { + *t4 = *t3 + 1; + return; + } + + *t4 = wantedT3 - foundT3 + *t3 + 1; +} + + /* a parser error function */ static void parseError(int c) @@ -11755,6 +11895,8 @@ static void startQuartetMaster(tree *tr, FILE *f) #endif + + static void computeAllThreeQuartets(tree *tr, nodeptr q1, nodeptr q2, int t1, int t2, int t3, int t4, FILE *f, analdef *adef) { /* set the tip nodes to different sequences @@ -11767,7 +11909,7 @@ static void computeAllThreeQuartets(tree *tr, nodeptr q1, nodeptr q2, int t1, in p4 = tr->nodep[t4]; double - l; + l; #ifdef _QUARTET_MPI quartetResult @@ -11833,32 +11975,210 @@ static void computeAllThreeQuartets(tree *tr, nodeptr q1, nodeptr q2, int t1, in #define RANDOM_QUARTETS 1 #define GROUPED_QUARTETS 2 +static void sampleQuartetsWithoutReplacementA(tree *tr, int numberOfTaxa, int64_t seed, uint64_t numberOfQuartets, uint64_t randomQuartets, nodeptr q1, nodeptr q2, uint64_t *prefixSumF2, uint64_t *prefixSumF3, uint64_t *prefixSumF4, FILE *f, analdef *adef, uint64_t actVal) +{ + int64_t + myseed = seed; + + uint64_t + sampleSize = randomQuartets, + quartetCounter = 0, + top = numberOfQuartets - sampleSize, + s; + + int + t1, + t2, + t3, + t4; + + double + NReal = (double)numberOfQuartets, + v, + quot; + + while(sampleSize >= 2) + { + v = randum(&myseed); + s = 0; + quot = top / NReal; + + while (quot > v) + { + s++; + top--; + NReal--; + quot = (quot * top) / NReal; + } + // Skip over the next s records and select the following one for the sample + actVal += s+1; + mapNumberToQuartet(numberOfTaxa, actVal, &t1, &t2, &t3, &t4, prefixSumF2, prefixSumF3, prefixSumF4); + computeAllThreeQuartets(tr, q1, q2, t1, t2, t3, t4, f, adef); + quartetCounter++; + + NReal--; + sampleSize--; + } + + // Special case sampleSize == 1 + s = trunc(round(NReal) * randum(&myseed)); + // Skip over the next s records and select the following one for the sample + actVal += s+1; + + mapNumberToQuartet(numberOfTaxa, actVal, &t1, &t2, &t3, &t4, prefixSumF2, prefixSumF3, prefixSumF4); + computeAllThreeQuartets(tr, q1, q2, t1, t2, t3, t4, f, adef); + quartetCounter++; + assert(quartetCounter == randomQuartets); +} + +static void sampleQuartetsWithoutReplacementD(tree *tr, int numberOfTaxa, int64_t seed, uint64_t numberOfQuartets, uint64_t randomQuartets, nodeptr q1, nodeptr q2, + uint64_t *prefixSumF2, uint64_t *prefixSumF3, uint64_t *prefixSumF4, FILE *f, analdef *adef, uint64_t actVal) +{ + int64_t + myseed = seed; + + uint64_t + sampleSize = randomQuartets, + quartetCounter = 0, + s, + qu1, + threshold, + t, + limit; + + int + t1, + t2, + t3, + t4; + + double + negalphainv = -1.0/13, + nreal = sampleSize, + ninv = 1.0 / nreal, + Nreal = numberOfQuartets, + vprime = exp(log(randum(&myseed)) * ninv), + qu1real, + nmin1inv, + x, + u, + negSreal, + y1, + y2, + top, + bottom; + + qu1 = -sampleSize + 1 + numberOfQuartets; + qu1real = -nreal + 1.0 + Nreal; + threshold = -negalphainv * sampleSize; + + while((sampleSize > 1) && (threshold < numberOfQuartets)) + { + nmin1inv = 1.0 / (-1.0 + nreal); + while(TRUE) + { + while (TRUE) + // step D2: Generate U and X + { + x = Nreal * (-vprime + 1.0); + s = trunc(x); + if (s < qu1) break; + vprime = exp(log(randum(&myseed)) * ninv); + } + u = randum(&myseed); + negSreal = (double) s * (-1); + // step D3: Accept? + y1 = exp(log(u * Nreal / qu1real) * nmin1inv); + vprime = y1 * (-x / Nreal + 1.0) * (qu1real / (negSreal + qu1real)); + if (vprime <= 1.0) break; // Accept! test (2.8) is true + // step D4: Accept? + y2 = 1.0; + top = -1.0 + Nreal; + if(-1 + sampleSize > s) + { + bottom = -nreal + Nreal; + limit = -s + numberOfQuartets; + } + else + { + bottom = -1.0 + negSreal + Nreal; + limit = qu1; + } + for (t = -1 + numberOfQuartets; t >= limit; t--) + { + y2 = (y2 * top)/bottom; + top = -1.0 + top; + bottom = -1.0 + bottom; + } + + if(Nreal / (-x + Nreal) >= y1 * exp(log(y2) * nmin1inv)) + { + // Accept! + vprime = exp(log(randum(&myseed)) * nmin1inv); + break; + } + vprime = exp(log(randum(&myseed)) * ninv); + } + // Step D5: Select the (s+1)st record + // Skip over the next s records and select the following one for the sample + actVal += s+1; + mapNumberToQuartet(numberOfTaxa, actVal, &t1, &t2, &t3, &t4, prefixSumF2, prefixSumF3, prefixSumF4); + computeAllThreeQuartets(tr, q1, q2, t1, t2, t3, t4, f, adef); + quartetCounter++; + + numberOfQuartets = -s + (-1 + numberOfQuartets); + Nreal = negSreal + (-1.0 + Nreal); + sampleSize--; + nreal = nreal - 1.0; + ninv = nmin1inv; + qu1 = qu1 - s; + qu1real += negSreal; + threshold += negalphainv; + } + if (sampleSize > 1) + { + // Use Method A to finish the sampling + assert(quartetCounter == randomQuartets - sampleSize); + sampleQuartetsWithoutReplacementA(tr, numberOfTaxa, seed, numberOfQuartets, sampleSize, q1, q2, prefixSumF2, prefixSumF3, prefixSumF4, f, adef, actVal); + } + else // Special case sampleSize == 1 + { + s = trunc(numberOfQuartets * vprime); + // Skip over the next s records and select the following one for the sample + actVal += s+1; + mapNumberToQuartet(numberOfTaxa, actVal, &t1, &t2, &t3, &t4, prefixSumF2, prefixSumF3, prefixSumF4); + computeAllThreeQuartets(tr, q1, q2, t1, t2, t3, t4, f, adef); + quartetCounter++; + assert(quartetCounter == randomQuartets); + } +} static void computeQuartets(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta) { /* some indices for generating quartets in an arbitrary way */ int - flavor = ALL_QUARTETS, - i, + flavor = ALL_QUARTETS, //type of quartet calculation + i, t1, t2, t3, t4, *groups[4], - groupSize[4]; - + groupSize[4]; + double - fraction = 0.0, + fraction = 0.0, //fraction of random quartets to compute t; uint64_t - randomQuartets = (uint64_t)(adef->multipleRuns), - quartetCounter = 0, - numberOfQuartets = ((uint64_t)tr->mxtips * ((uint64_t)tr->mxtips - 1) * ((uint64_t)tr->mxtips - 2) * ((uint64_t)tr->mxtips - 3)) / 24; - - /* use two inner nodes for building quartet trees */ + randomQuartets = (uint64_t)(adef->multipleRuns), //number of random quartets to compute + quartetCounter = 0, + //total number of possible quartets, note that we count the following ((A,B),(C,D)), ((A,C),(B,D)), ((A,D),(B,C)) as one quartet here + numberOfQuartets = ((uint64_t)tr->mxtips * ((uint64_t)tr->mxtips - 1) * ((uint64_t)tr->mxtips - 2) * ((uint64_t)tr->mxtips - 3)) / 24; + + /* use two inner tree nodes for building quartet trees */ nodeptr q1 = tr->nodep[tr->mxtips + 1], @@ -11870,6 +12190,11 @@ static void computeQuartets(tree *tr, analdef *adef, rawdata *rdta, cruncheddata FILE *f; + /***********************************/ + + + + /* build output file name */ strcpy(quartetFileName, workdir); @@ -11878,8 +12203,6 @@ static void computeQuartets(tree *tr, analdef *adef, rawdata *rdta, cruncheddata /* open output file */ - - #ifdef _QUARTET_MPI if(processID == 0) #endif @@ -11894,22 +12217,29 @@ static void computeQuartets(tree *tr, analdef *adef, rawdata *rdta, cruncheddata if(!adef->useBinaryModelFile) { #ifdef _QUARTET_MPI + //the parallel version requires a pre-computed model parameter file as input! assert(0); #endif - /* get a starting tree: either reads in a tree or computes a randomized stepwise addition parsimony tree */ + /* get a starting tree on which we optimize the likelihood model parameters: either reads in a tree or computes a randomized stepwise addition parsimony tree */ getStartingTree(tr, adef); - /* optimize model parameters on that comprehensive tree that can subsequently be used for qyartet building */ -#ifndef __BLACKRIM + /* optimize model parameters on that comprehensive tree that can subsequently be used for evaluation of quartet likelihoods */ + +#ifndef __BLACKRIM //if BLACKRIM is defined, the model parameters will be optimized for each quartet individually modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); #endif printBothOpen("Time for parsing input tree or building parsimony tree and optimizing model parameters: %f\n\n", gettime() - masterTime); + +#ifndef __BLACKRIM + printBothOpen("Tree likelihood: %f\n\n", tr->likelihood); +#endif } else { + //if a binary model parameter file has been specified, we just read the model parameters from there readBinaryModel(tr, adef); printBothOpen("Time for reading model parameters: %f\n\n", gettime() - masterTime); @@ -11920,10 +12250,17 @@ static void computeQuartets(tree *tr, analdef *adef, rawdata *rdta, cruncheddata if(adef->useQuartetGrouping) { + //quartet grouping evaluates all possible quartets from four disjoint + //sets of user-specified taxon names + flavor = GROUPED_QUARTETS; + + //parse the four disjoint sets of taxon names specified by the user from file groupingParser(quartetGroupingFileName, groups, groupSize, tr); #ifdef __BLACKRIM + //special implementation where we only sub-sample randomly from the quartets + //defined by the four user-specified taxon sets numberOfQuartets = (uint64_t)groupSize[0] * (uint64_t)groupSize[1] * (uint64_t)groupSize[2] * (uint64_t)groupSize[3]; if(randomQuartets > numberOfQuartets) @@ -11936,13 +12273,18 @@ static void computeQuartets(tree *tr, analdef *adef, rawdata *rdta, cruncheddata } else { + //if the user specified more random quartets to sample than there actually + //exist for the number of taxa, then fix this. if(randomQuartets > numberOfQuartets) randomQuartets = 1; if(randomQuartets == 1) + //change flavor if randomQuartets > possibleQuartets flavor = ALL_QUARTETS; else { + //compute the fraction of random quartets to sample + //there may be an issue here with the unit64_t <-> double cast fraction = (double)randomQuartets / (double)numberOfQuartets; flavor = RANDOM_QUARTETS; } @@ -11999,6 +12341,9 @@ static void computeQuartets(tree *tr, analdef *adef, rawdata *rdta, cruncheddata tr->mxtips is the maximum number of tips in the alignment/tree */ + + //now do the respective quartet evaluations by switching over the three distinct flavors + #ifdef _QUARTET_MPI if(processID > 0) #endif @@ -12027,32 +12372,77 @@ static void computeQuartets(tree *tr, analdef *adef, rawdata *rdta, cruncheddata } break; case RANDOM_QUARTETS: - { - /* randomly sub-sample a fraction of all quartets */ + { + //code contributed by Sarah for drawing quartets without replacement :-) - for(t1 = 1; t1 <= tr->mxtips; t1++) - for(t2 = t1 + 1; t2 <= tr->mxtips; t2++) - for(t3 = t2 + 1; t3 <= tr->mxtips; t3++) - for(t4 = t3 + 1; t4 <= tr->mxtips; t4++) - { - double - r = randum(&adef->parsimonySeed); - - if(r < fraction) - { + if(adef->sampleQuartetsWithoutReplacement) + { + uint64_t + *prefixSumF2 = (uint64_t*)rax_malloc(sizeof(uint64_t) * (size_t)(tr->mxtips - 2)), + *prefixSumF3 = (uint64_t*)rax_malloc(sizeof(uint64_t) * (size_t)(tr->mxtips - 2)), + *prefixSumF4 = (uint64_t*)rax_malloc(sizeof(uint64_t) * (size_t)(tr->mxtips - 2)); + + preprocessQuartetPrefix(tr->mxtips, prefixSumF2, prefixSumF3, prefixSumF4); + + if(randomQuartets >= numberOfQuartets / 13) + sampleQuartetsWithoutReplacementA(tr, tr->mxtips, adef->parsimonySeed, numberOfQuartets, randomQuartets, q1, q2, prefixSumF2, prefixSumF3, prefixSumF4, f, adef, 0); + else + sampleQuartetsWithoutReplacementD(tr, tr->mxtips, adef->parsimonySeed, numberOfQuartets, randomQuartets, q1, q2, prefixSumF2, prefixSumF3, prefixSumF4, f, adef, 0); + + rax_free(prefixSumF2); + rax_free(prefixSumF3); + rax_free(prefixSumF4); + } + else + { + //endless loop ta make sure we randomly sub-sample exactly as many quartets as the user specified + + //This is not very elegant, but it works, note however, that especially when the number of + //random quartets to be sampled is large, that is, close to the total number of quartets + //some quartets may be sampled twice by pure chance. To randomly sample unique quartets + //using hashes or bitmaps to store which quartets have already been sampled is not memory efficient. + //Insetad, we need to use a random number generator that can generate a unique series of random numbers + //and then have a function f() that maps those random numbers to the corresponding index quartet (t1, t2, t3, t4), + //see above + + do + { + //loop over all quartets + for(t1 = 1; t1 <= tr->mxtips; t1++) + for(t2 = t1 + 1; t2 <= tr->mxtips; t2++) + for(t3 = t2 + 1; t3 <= tr->mxtips; t3++) + for(t4 = t3 + 1; t4 <= tr->mxtips; t4++) + { + //chose a random number + double + r = randum(&adef->parsimonySeed); + + //if the random number is smaller than the fraction of quartets to subsample + //evaluate the likelihood of the current quartet + if(r < fraction) + { #ifdef _QUARTET_MPI - if((quartetCounter % (uint64_t)(processes - 1)) == (uint64_t)(processID - 1)) + //MPI version very simple and naive way to determine which processor + //is goingt to do the likelihood calculations for this quartet + if((quartetCounter % (uint64_t)(processes - 1)) == (uint64_t)(processID - 1)) #endif - computeAllThreeQuartets(tr, q1, q2, t1, t2, t3, t4, f, adef); - quartetCounter++; - } - - if(quartetCounter == randomQuartets) - goto DONE; - } - - DONE: - assert(quartetCounter == randomQuartets); + //function that computes the likelihood for all three possible unrooted trees + //defined by the given quartet of taxa + computeAllThreeQuartets(tr, q1, q2, t1, t2, t3, t4, f, adef); + //increment quartet counter that counts how many quartets we have evaluated + quartetCounter++; + } + + //exit endless loop if we have randomly sub-sampled as many quartets as the user specified + if(quartetCounter == randomQuartets) + goto DONE; + } + } + while(1); + + DONE: + assert(quartetCounter == randomQuartets); + } } break; case GROUPED_QUARTETS: @@ -12115,6 +12505,8 @@ static void computeQuartets(tree *tr, analdef *adef, rawdata *rdta, cruncheddata } #endif + + t = gettime() - t; printBothOpen("\nPure quartet computation time: %f secs\n", t); diff --git a/axml.h b/axml.h index ee7933a..58aeb8a 100644 --- a/axml.h +++ b/axml.h @@ -168,9 +168,9 @@ #define PointGamma(prob,alpha,beta) PointChi2(prob,2.0*(alpha))/(2.0*(beta)) #define programName "RAxML" -#define programVersion "8.2.7" -#define programVersionInt 8270 -#define programDate "March 10 2016" +#define programVersion "8.2.8" +#define programVersionInt 8280 +#define programDate "March 23 2016" #define TREE_EVALUATION 0 @@ -1175,6 +1175,7 @@ typedef struct { boolean setThreadAffinity; int bootstopPermutations; int fcThreshold; + boolean sampleQuartetsWithoutReplacement; } analdef; diff --git a/fastDNAparsimony.c b/fastDNAparsimony.c index 9d81560..dda14f4 100644 --- a/fastDNAparsimony.c +++ b/fastDNAparsimony.c @@ -253,13 +253,13 @@ void newviewParsimonyIterativeFast(tree *tr) parsimonyNumber *left[2], *right[2], - *this[2]; + *thisOne[2]; for(k = 0; k < 2; k++) { left[k] = &(tr->partitionData[model].parsVect[(width * 2 * qNumber) + width * k]); right[k] = &(tr->partitionData[model].parsVect[(width * 2 * rNumber) + width * k]); - this[k] = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]); + thisOne[k] = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]); } for(i = 0; i < width; i += INTS_PER_VECTOR) @@ -281,8 +281,8 @@ void newviewParsimonyIterativeFast(tree *tr) v_N = VECTOR_BIT_OR(l_A, l_C); - VECTOR_STORE((CAST)(&this[0][i]), VECTOR_BIT_OR(l_A, VECTOR_AND_NOT(v_N, v_A))); - VECTOR_STORE((CAST)(&this[1][i]), VECTOR_BIT_OR(l_C, VECTOR_AND_NOT(v_N, v_C))); + VECTOR_STORE((CAST)(&thisOne[0][i]), VECTOR_BIT_OR(l_A, VECTOR_AND_NOT(v_N, v_A))); + VECTOR_STORE((CAST)(&thisOne[1][i]), VECTOR_BIT_OR(l_C, VECTOR_AND_NOT(v_N, v_C))); v_N = VECTOR_AND_NOT(v_N, allOne); @@ -295,13 +295,13 @@ void newviewParsimonyIterativeFast(tree *tr) parsimonyNumber *left[4], *right[4], - *this[4]; + *thisOne[4]; for(k = 0; k < 4; k++) { left[k] = &(tr->partitionData[model].parsVect[(width * 4 * qNumber) + width * k]); right[k] = &(tr->partitionData[model].parsVect[(width * 4 * rNumber) + width * k]); - this[k] = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]); + thisOne[k] = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]); } for(i = 0; i < width; i += INTS_PER_VECTOR) @@ -333,10 +333,10 @@ void newviewParsimonyIterativeFast(tree *tr) v_N = VECTOR_BIT_OR(VECTOR_BIT_OR(l_A, l_C), VECTOR_BIT_OR(l_G, l_T)); - VECTOR_STORE((CAST)(&this[0][i]), VECTOR_BIT_OR(l_A, VECTOR_AND_NOT(v_N, v_A))); - VECTOR_STORE((CAST)(&this[1][i]), VECTOR_BIT_OR(l_C, VECTOR_AND_NOT(v_N, v_C))); - VECTOR_STORE((CAST)(&this[2][i]), VECTOR_BIT_OR(l_G, VECTOR_AND_NOT(v_N, v_G))); - VECTOR_STORE((CAST)(&this[3][i]), VECTOR_BIT_OR(l_T, VECTOR_AND_NOT(v_N, v_T))); + VECTOR_STORE((CAST)(&thisOne[0][i]), VECTOR_BIT_OR(l_A, VECTOR_AND_NOT(v_N, v_A))); + VECTOR_STORE((CAST)(&thisOne[1][i]), VECTOR_BIT_OR(l_C, VECTOR_AND_NOT(v_N, v_C))); + VECTOR_STORE((CAST)(&thisOne[2][i]), VECTOR_BIT_OR(l_G, VECTOR_AND_NOT(v_N, v_G))); + VECTOR_STORE((CAST)(&thisOne[3][i]), VECTOR_BIT_OR(l_T, VECTOR_AND_NOT(v_N, v_T))); v_N = VECTOR_AND_NOT(v_N, allOne); @@ -349,13 +349,13 @@ void newviewParsimonyIterativeFast(tree *tr) parsimonyNumber *left[20], *right[20], - *this[20]; + *thisOne[20]; for(k = 0; k < 20; k++) { left[k] = &(tr->partitionData[model].parsVect[(width * 20 * qNumber) + width * k]); right[k] = &(tr->partitionData[model].parsVect[(width * 20 * rNumber) + width * k]); - this[k] = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]); + thisOne[k] = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]); } for(i = 0; i < width; i += INTS_PER_VECTOR) @@ -379,7 +379,7 @@ void newviewParsimonyIterativeFast(tree *tr) } for(j = 0; j < 20; j++) - VECTOR_STORE((CAST)(&this[j][i]), VECTOR_BIT_OR(l_A[j], VECTOR_AND_NOT(v_N, v_A[j]))); + VECTOR_STORE((CAST)(&thisOne[j][i]), VECTOR_BIT_OR(l_A[j], VECTOR_AND_NOT(v_N, v_A[j]))); v_N = VECTOR_AND_NOT(v_N, allOne); @@ -392,7 +392,7 @@ void newviewParsimonyIterativeFast(tree *tr) parsimonyNumber *left[32], *right[32], - *this[32]; + *thisOne[32]; assert(states <= 32); @@ -400,7 +400,7 @@ void newviewParsimonyIterativeFast(tree *tr) { left[k] = &(tr->partitionData[model].parsVect[(width * states * qNumber) + width * k]); right[k] = &(tr->partitionData[model].parsVect[(width * states * rNumber) + width * k]); - this[k] = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]); + thisOne[k] = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]); } for(i = 0; i < width; i += INTS_PER_VECTOR) @@ -424,7 +424,7 @@ void newviewParsimonyIterativeFast(tree *tr) } for(j = 0; j < states; j++) - VECTOR_STORE((CAST)(&this[j][i]), VECTOR_BIT_OR(l_A[j], VECTOR_AND_NOT(v_N, v_A[j]))); + VECTOR_STORE((CAST)(&thisOne[j][i]), VECTOR_BIT_OR(l_A[j], VECTOR_AND_NOT(v_N, v_A[j]))); v_N = VECTOR_AND_NOT(v_N, allOne); @@ -646,7 +646,7 @@ void newviewParsimonyIterativeFast(tree *tr) parsimonyNumber *left[2], *right[2], - *this[2]; + *thisOne[2]; parsimonyNumber o_A, @@ -659,7 +659,7 @@ void newviewParsimonyIterativeFast(tree *tr) { left[k] = &(tr->partitionData[model].parsVect[(width * 2 * qNumber) + width * k]); right[k] = &(tr->partitionData[model].parsVect[(width * 2 * rNumber) + width * k]); - this[k] = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]); + thisOne[k] = &(tr->partitionData[model].parsVect[(width * 2 * pNumber) + width * k]); } for(i = 0; i < width; i++) @@ -672,8 +672,8 @@ void newviewParsimonyIterativeFast(tree *tr) t_N = ~(t_A | t_C); - this[0][i] = t_A | (t_N & o_A); - this[1][i] = t_C | (t_N & o_C); + thisOne[0][i] = t_A | (t_N & o_A); + thisOne[1][i] = t_C | (t_N & o_C); totalScore += populationCount(t_N); } @@ -684,13 +684,13 @@ void newviewParsimonyIterativeFast(tree *tr) parsimonyNumber *left[4], *right[4], - *this[4]; + *thisOne[4]; for(k = 0; k < 4; k++) { left[k] = &(tr->partitionData[model].parsVect[(width * 4 * qNumber) + width * k]); right[k] = &(tr->partitionData[model].parsVect[(width * 4 * rNumber) + width * k]); - this[k] = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]); + thisOne[k] = &(tr->partitionData[model].parsVect[(width * 4 * pNumber) + width * k]); } parsimonyNumber @@ -718,10 +718,10 @@ void newviewParsimonyIterativeFast(tree *tr) t_N = ~(t_A | t_C | t_G | t_T); - this[0][i] = t_A | (t_N & o_A); - this[1][i] = t_C | (t_N & o_C); - this[2][i] = t_G | (t_N & o_G); - this[3][i] = t_T | (t_N & o_T); + thisOne[0][i] = t_A | (t_N & o_A); + thisOne[1][i] = t_C | (t_N & o_C); + thisOne[2][i] = t_G | (t_N & o_G); + thisOne[3][i] = t_T | (t_N & o_T); totalScore += populationCount(t_N); } @@ -732,7 +732,7 @@ void newviewParsimonyIterativeFast(tree *tr) parsimonyNumber *left[20], *right[20], - *this[20]; + *thisOne[20]; parsimonyNumber o_A[20], @@ -743,7 +743,7 @@ void newviewParsimonyIterativeFast(tree *tr) { left[k] = &(tr->partitionData[model].parsVect[(width * 20 * qNumber) + width * k]); right[k] = &(tr->partitionData[model].parsVect[(width * 20 * rNumber) + width * k]); - this[k] = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]); + thisOne[k] = &(tr->partitionData[model].parsVect[(width * 20 * pNumber) + width * k]); } for(i = 0; i < width; i++) @@ -762,7 +762,7 @@ void newviewParsimonyIterativeFast(tree *tr) t_N = ~t_N; for(k = 0; k < 20; k++) - this[k][i] = t_A[k] | (t_N & o_A[k]); + thisOne[k][i] = t_A[k] | (t_N & o_A[k]); totalScore += populationCount(t_N); } @@ -773,7 +773,7 @@ void newviewParsimonyIterativeFast(tree *tr) parsimonyNumber *left[32], *right[32], - *this[32]; + *thisOne[32]; parsimonyNumber o_A[32], @@ -786,7 +786,7 @@ void newviewParsimonyIterativeFast(tree *tr) { left[k] = &(tr->partitionData[model].parsVect[(width * states * qNumber) + width * k]); right[k] = &(tr->partitionData[model].parsVect[(width * states * rNumber) + width * k]); - this[k] = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]); + thisOne[k] = &(tr->partitionData[model].parsVect[(width * states * pNumber) + width * k]); } for(i = 0; i < width; i++) @@ -803,7 +803,7 @@ void newviewParsimonyIterativeFast(tree *tr) t_N = ~t_N; for(k = 0; k < states; k++) - this[k][i] = t_A[k] | (t_N & o_A[k]); + thisOne[k][i] = t_A[k] | (t_N & o_A[k]); totalScore += populationCount(t_N); } diff --git a/treeIO.c b/treeIO.c index f323d8b..a9ded68 100644 --- a/treeIO.c +++ b/treeIO.c @@ -1120,9 +1120,16 @@ static boolean addElementLen (FILE *fp, tree *tr, nodeptr p, boolean readBranchL endCounter, branchLabel = -1; - if (! treeNeedCh(fp, ':', "in")) return FALSE; - if (! treeProcessLength(fp, &branch, &branchLabel, storeBranchLabels, tr)) - return FALSE; + if (! treeNeedCh(fp, ':', "in")) + { + printf("ERROR: problem reading branch length ... RAxML will abort with a failing assertion\n\n"); + return FALSE; + } + if (! treeProcessLength(fp, &branch, &branchLabel, storeBranchLabels, tr)) + { + printf("ERROR: problem reading branch length ... RAxML will abort with a failing assertion\n\n"); + return FALSE; + } endCounter = tr->branchLabelCounter; -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/raxml.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
