Skip to content

Commit

Permalink
-added an automatic option to compute the SH-like support values (-f …
Browse files Browse the repository at this point in the history 
…J option) on a given tree

 for each partition individually. If you invoke -f J with a partition file, the program will
 automatically also print out an additional tree file where the SH-like support values
 are displayed as branch labels in the form [s1,s2,...,sN] where s1...sN are the SH-like support values
 for each partition.
 At first glance this seems to be quite interesting because we already observed huge differences
 between the per-partition SH-like support values on a small test dataset.
  • Loading branch information
@stamatak
stamatak committed Jul 16, 2013
1 parent 7df7678 commit b355115
Show file tree
Hide file tree
Showing 9 changed files with 140 additions and 56 deletions.
2 changes: 1 addition & 1 deletion README
View file
@@ -1,4 +1,4 @@
Standard RAxML version 7.6.9
Standard RAxML version 7.7.0

Compiling under Linux:

Expand Down
2 changes: 1 addition & 1 deletion Release-Notes
View file
@@ -1,4 +1,4 @@
RAxML v7.6.9
RAxML v7.7.0

We view RAxML code development as permanent work in progress.

Expand Down
28 changes: 12 additions & 16 deletions axml.c
View file
Expand Up @@ -5964,7 +5964,7 @@ void printResult(tree *tr, analdef *adef, boolean finalPrint)
case TREE_EVALUATION:


Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, finalPrint, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, finalPrint, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE, FALSE);

logFile = myfopen(temporaryFileName, "wb");
fprintf(logFile, "%s", tr->tree_string);
Expand Down Expand Up @@ -5992,8 +5992,7 @@ void printResult(tree *tr, analdef *adef, boolean finalPrint)
{
case GAMMA:
case GAMMA_I:
Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, finalPrint, adef,
SUMMARIZE_LH, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, finalPrint, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE, FALSE);

logFile = myfopen(temporaryFileName, "wb");
fprintf(logFile, "%s", tr->tree_string);
Expand All @@ -6003,8 +6002,7 @@ void printResult(tree *tr, analdef *adef, boolean finalPrint)
printTreePerGene(tr, adef, temporaryFileName, "wb");
break;
case CAT:
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef,
NO_BRANCHES, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, FALSE, FALSE);

logFile = myfopen(temporaryFileName, "wb");
fprintf(logFile, "%s", tr->tree_string);
Expand All @@ -6017,8 +6015,7 @@ void printResult(tree *tr, analdef *adef, boolean finalPrint)
}
else
{
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef,
NO_BRANCHES, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, FALSE, FALSE);
logFile = myfopen(temporaryFileName, "wb");
fprintf(logFile, "%s", tr->tree_string);
fclose(logFile);
Expand Down Expand Up @@ -6048,7 +6045,7 @@ void printBootstrapResult(tree *tr, analdef *adef, boolean finalPrint)
{
if(adef->bootstrapBranchLengths)
{
Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, finalPrint, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, finalPrint, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE, FALSE);

logFile = myfopen(fileName, "ab");
fprintf(logFile, "%s", tr->tree_string);
Expand All @@ -6059,7 +6056,7 @@ void printBootstrapResult(tree *tr, analdef *adef, boolean finalPrint)
}
else
{
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, FALSE, FALSE);

logFile = myfopen(fileName, "ab");
fprintf(logFile, "%s", tr->tree_string);
Expand All @@ -6084,15 +6081,15 @@ void printBipartitionResult(tree *tr, analdef *adef, boolean finalPrint, boolean

if(!printIC)
{
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, TRUE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, printIC);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, TRUE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, printIC, FALSE);

logFile = myfopen(bipartitionsFileName, "ab");

fprintf(logFile, "%s", tr->tree_string);
fclose(logFile);
}

Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, TRUE, FALSE, printIC);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, TRUE, FALSE, printIC, FALSE);

if(printIC)
logFile = myfopen(icFileNameBranchLabels, "ab");
Expand Down Expand Up @@ -6168,7 +6165,7 @@ void printLog(tree *tr, analdef *adef, boolean finalPrint)
sprintf(treeID, "%d", tr->checkPointCounter);
strcat(checkPoints, treeID);

Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, FALSE, FALSE);

logFile = myfopen(checkPoints, "ab");
fprintf(logFile, "%s", tr->tree_string);
Expand Down Expand Up @@ -6198,7 +6195,7 @@ void printStartingTree(tree *tr, analdef *adef, boolean finalPrint)
FILE *treeFile;
char temporaryFileName[1024] = "", treeID[64] = "";

Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, finalPrint, adef, NO_BRANCHES, FALSE, FALSE, FALSE, FALSE);

if(adef->randomStartingTree)
strcpy(temporaryFileName, randomFileName);
Expand Down Expand Up @@ -8543,8 +8540,7 @@ static void computeAllLHs(tree *tr, analdef *adef, char *bootStrapFileName)
list[i].tree = i;
list[i].lh = tr->likelihood;

Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE,
TRUE, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, TRUE, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE, FALSE);

fprintf(result, "%s", tr->tree_string);

Expand Down Expand Up @@ -10004,7 +10000,7 @@ static void thoroughTreeOptimization(tree *tr, analdef *adef, rawdata *rdta, cru
strcat(bestTreeFileName, "RAxML_bestTree.");
strcat(bestTreeFileName, run_id);

Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, TRUE, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, TRUE, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE, FALSE);
f = myfopen(bestTreeFileName, "wb");
fprintf(f, "%s", tr->tree_string);
fclose(f);
Expand Down
7 changes: 4 additions & 3 deletions axml.h
View file
Expand Up @@ -163,7 +163,7 @@
#define PointGamma(prob,alpha,beta) PointChi2(prob,2.0*(alpha))/(2.0*(beta))

#define programName "RAxML"
#define programVersion "7.6.9"
#define programVersion "7.7.0"
#define programDate "July 16 2013"


Expand Down Expand Up @@ -453,7 +453,8 @@ typedef struct
epaBranchData *epa;

unsigned int *vector;
int support;
int support;
int *supports;
double ic;
double icAll;
struct noderec *oP;
Expand Down Expand Up @@ -1195,7 +1196,7 @@ extern boolean freeBestTree ( bestlist *bt );


extern char *Tree2String ( char *treestr, tree *tr, nodeptr p, boolean printBranchLengths, boolean printNames, boolean printLikelihood,
boolean rellTree, boolean finalPrint, analdef *adef, int perGene, boolean branchLabelSupport, boolean printSHSupport, boolean printIC);
boolean rellTree, boolean finalPrint, analdef *adef, int perGene, boolean branchLabelSupport, boolean printSHSupport, boolean printIC, boolean printSHSupports);
extern void printTreePerGene(tree *tr, analdef *adef, char *fileName, char *permission);


Expand Down
102 changes: 85 additions & 17 deletions fastSearch.c
View file
Expand Up @@ -749,7 +749,7 @@ static void storeBranches(tree *tr, nodeptr p, double *pqz, double *pz1, double
}


static int SHSupport(int nPos, int nBootstrap, int *col, double loglk[3], double *siteloglk[3])
static int SHSupport(int nPos, int nBootstrap, int *col, double loglk[3], double *siteloglk[3], int lower, int upper, boolean perPartition)
{
double
resampleDelta,
Expand Down Expand Up @@ -777,15 +777,17 @@ static int SHSupport(int nPos, int nBootstrap, int *col, double loglk[3], double
diff = ABS(loglk[1] - loglk[0]);
/*printf("%1.80f %1.80f\n", loglk[0], loglk[1]);*/
shittySplit = TRUE;
assert(diff < 0.1);
if(!perPartition)
assert(diff < 0.1);
}

if(loglk[2] >= loglk[0])
{
diff = ABS(loglk[2] - loglk[0]);
/*printf("%1.80f %1.80f\n", loglk[0], loglk[2]);*/
shittySplit = TRUE;
assert(diff < 0.1);
shittySplit = TRUE;
if(!perPartition)
assert(diff < 0.1);
}

if(loglk[0] > loglk[2] && loglk[0] > loglk[1])
Expand All @@ -800,6 +802,7 @@ static int SHSupport(int nPos, int nBootstrap, int *col, double loglk[3], double

if(shittySplit)
return 0;


/*
printf("%f %f %f\n", loglk[0], loglk[1], loglk[2]);
Expand All @@ -813,9 +816,11 @@ static int SHSupport(int nPos, int nBootstrap, int *col, double loglk[3], double
for (i = 0; i < 3; i++)
resampled[i] = -loglk[i];

for (j = 0; j < nPos; j++)
for (j = lower; j < upper; j++)
{
int pos = col[iBoot * nPos + j];
int
pos = col[iBoot * nPos + j];

for (i = 0; i < 3; i++)
resampled[i] += pos * siteloglk[i][j];
}
Expand Down Expand Up @@ -886,14 +891,20 @@ static void doNNIs(tree *tr, nodeptr p, double *lhVectors[3], boolean shSupport,
if(!isTip(q->number, tr->mxtips))
{
int
model,
whichNNI = 0;

nodeptr
qb1 = q->next->back,
qb2 = q->next->next->back;

double
lh[3];
lh[3],
*lhv[3];

lhv[0] = (double*)rax_malloc(sizeof(double) * tr->NumberOfModels);
lhv[1] = (double*)rax_malloc(sizeof(double) * tr->NumberOfModels);
lhv[2] = (double*)rax_malloc(sizeof(double) * tr->NumberOfModels);

*innerBranches = *innerBranches + 1;

Expand All @@ -908,6 +919,9 @@ static void doNNIs(tree *tr, nodeptr p, double *lhVectors[3], boolean shSupport,
evaluateGeneric(tr, p);

lh[0] = tr->likelihood;

for(model = 0; model < tr->NumberOfModels; model++)
lhv[0][model] = tr->perPartitionLH[model];

storeBranches(tr, p, pqz_0, pz1_0, pz2_0, qz1_0, qz2_0);

Expand All @@ -934,7 +948,10 @@ static void doNNIs(tree *tr, nodeptr p, double *lhVectors[3], boolean shSupport,
else
evaluateGeneric(tr, p);

lh[1] = tr->likelihood;
lh[1] = tr->likelihood;

for(model = 0; model < tr->NumberOfModels; model++)
lhv[1][model] = tr->perPartitionLH[model];

storeBranches(tr, p, pqz_1, pz1_1, pz2_1, qz1_1, qz2_1);

Expand Down Expand Up @@ -966,6 +983,9 @@ static void doNNIs(tree *tr, nodeptr p, double *lhVectors[3], boolean shSupport,

lh[2] = tr->likelihood;

for(model = 0; model < tr->NumberOfModels; model++)
lhv[2][model] = tr->perPartitionLH[model];

storeBranches(tr, p, pqz_2, pz1_2, pz2_2, qz1_2, qz2_2);

if(lh[2] > lh[0] && lh[2] > lh[1])
Expand Down Expand Up @@ -1016,7 +1036,31 @@ static void doNNIs(tree *tr, nodeptr p, double *lhVectors[3], boolean shSupport,
*interchanges = *interchanges + 1;

if(shSupport)
p->bInf->support = SHSupport(tr->cdta->endsite, 1000, tr->resample, lh, lhVectors);
{
p->bInf->support = SHSupport(tr->cdta->endsite, 1000, tr->resample, lh, lhVectors, 0, tr->cdta->endsite, FALSE);

//printf("ALL: %f %f %f ", lh[0], lh[1], lh[2]);

for(model = 0; model < tr->NumberOfModels; model++)
{
double
perPartitionLH[3];

perPartitionLH[0] = lhv[0][model];
perPartitionLH[1] = lhv[1][model];
perPartitionLH[2] = lhv[2][model];

p->bInf->supports[model] = SHSupport(tr->cdta->endsite, 1000, tr->resample, perPartitionLH, lhVectors, tr->partitionData[model].lower, tr->partitionData[model].upper, TRUE);

//printf(" p%d %f %f %f ", model, perPartitionLH[0], perPartitionLH[1], perPartitionLH[2]);
}

//printf("\n");
}

rax_free(lhv[0]);
rax_free(lhv[1]);
rax_free(lhv[2]);
}


Expand Down Expand Up @@ -1258,7 +1302,7 @@ void fastSearch(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta)



Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, FALSE, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, FALSE, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE, FALSE);

f = myfopen(bestTreeFileName, "wb");
fprintf(f, "%s", tr->tree_string);
Expand Down Expand Up @@ -1290,17 +1334,23 @@ void shSupports(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta)
*f;

int
i,
interchanges = 0,
counter = 0;

assert(adef->restart);

tr->resample = permutationSH(tr, 1000, 12345);


//a tiny bit dirty here, all allocs should be cleaned up!

lhVectors[0] = (double *)rax_malloc(sizeof(double) * tr->cdta->endsite);
lhVectors[1] = (double *)rax_malloc(sizeof(double) * tr->cdta->endsite);
lhVectors[2] = (double *)rax_malloc(sizeof(double) * tr->cdta->endsite);
tr->bInf = (branchInfo*)rax_malloc(sizeof(branchInfo) * (tr->mxtips - 3));
tr->bInf = (branchInfo*)rax_malloc(sizeof(branchInfo) * (tr->mxtips - 3));

for(i = 0; i < tr->mxtips - 3; i++)
tr->bInf[i].supports = (int*)rax_malloc(sizeof(int) * tr->NumberOfModels);

initModel(tr, rdta, cdta, adef);

Expand All @@ -1316,7 +1366,7 @@ void shSupports(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta)
else
{
evaluateGenericInitrav(tr, tr->start);
modOpt(tr, adef, FALSE, 10.0);
modOpt(tr, adef, FALSE, 1.0);
}

printBothOpen("Time after model optimization: %f\n", gettime() - masterTime);
Expand Down Expand Up @@ -1354,7 +1404,7 @@ void shSupports(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta)
strcat(bestTreeFileName, "RAxML_fastTree.");
strcat(bestTreeFileName, run_id);

Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, FALSE, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, FALSE, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE, FALSE);

f = myfopen(bestTreeFileName, "wb");
fprintf(f, "%s", tr->tree_string);
Expand All @@ -1365,17 +1415,35 @@ void shSupports(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta)
strcat(shSupportFileName, "RAxML_fastTreeSH_Support.");
strcat(shSupportFileName, run_id);

Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, FALSE, adef, SUMMARIZE_LH, FALSE, TRUE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, FALSE, adef, SUMMARIZE_LH, FALSE, TRUE, FALSE, FALSE);

f = myfopen(shSupportFileName, "wb");
fprintf(f, "%s", tr->tree_string);
fclose(f);

printBothOpen("RAxML NNI-optimized tree written to file: %s\n", bestTreeFileName);

printBothOpen("Same tree with SH-like supports written to file: %s\n", shSupportFileName);
printBothOpen("\nSame tree with SH-like supports written to file: %s\n", shSupportFileName);

if(tr->NumberOfModels > 1)
{
char
shSupportPerPartitionFileName[1024];

strcpy(shSupportPerPartitionFileName, workdir);
strcat(shSupportPerPartitionFileName, "RAxML_fastTree_perPartition_SH_Support.");
strcat(shSupportPerPartitionFileName, run_id);

printBothOpen("Total execution time: %f\n", gettime() - masterTime);
Tree2String(tr->tree_string, tr, tr->start->back, TRUE, TRUE, FALSE, FALSE, FALSE, adef, SUMMARIZE_LH, FALSE, FALSE, FALSE, TRUE);

f = myfopen(shSupportPerPartitionFileName, "wb");
fprintf(f, "%s", tr->tree_string);
fclose(f);

printBothOpen("\nSame tree with SH-like support for each partition written to file: %s\n", shSupportPerPartitionFileName);
}

printBothOpen("\nTotal execution time: %f\n", gettime() - masterTime);

exit(0);
}
2 changes: 1 addition & 1 deletion leaveDropping.c
View file
Expand Up @@ -1642,7 +1642,7 @@ void computeRogueTaxa(tree *tr, char* treeSetFileName, analdef *adef)

assert((unsigned)tips == ((unsigned)tr->mxtips - droppedTaxaNum));

Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, FALSE, adef, NO_BRANCHES, FALSE, FALSE, FALSE);
Tree2String(tr->tree_string, tr, tr->start->back, FALSE, TRUE, FALSE, FALSE, FALSE, adef, NO_BRANCHES, FALSE, FALSE, FALSE, FALSE);
fprintf(outf, "%s", tr->tree_string);

/*printf("%u %u\n", (unsigned)tips, ((unsigned)tr->mxtips - droppedTaxaNum));*/
Expand Down

0 comments on commit b355115

Please sign in to comment.