[med-svn] [raxml] 01/06: Imported Upstream version 8.2.8+dfsg
Andreas Tille
tille at debian.org
Thu May 19 16:36:23 UTC 2016
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 <tille at debian.org>
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
More information about the debian-med-commit
mailing list