/* RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees * Copyright August 2006 by Alexandros Stamatakis * * Partially derived from * fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen * * and * * Programs of the PHYLIP package by Joe Felsenstein. * * This program is free software; you may redistribute it and/or modify its * under the terms of the GNU General Public License as published by the Free * Software Foundation; either version 2 of the License, or (at your option) * any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License * for more details. * * * For any other enquiries send an Email to Alexandros Stamatakis * Alexandros.Stamatakis@epfl.ch * * When publishing work that is based on the results from RAxML-VI-HPC please cite: * * Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models". * Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446 */ #ifndef WIN32 #include #endif #include #include #include #include #include #include #include "axml.h" #ifdef _WAYNE_MPI #include extern int processID; extern int processes; #endif extern int optimizeRatesInvocations; extern int optimizeRateCategoryInvocations; extern int optimizeAlphaInvocations; extern int optimizeInvarInvocations; extern int checkPointCounter; extern int Thorough; extern int partCount; extern char tree_file[1024]; extern const unsigned int mask32[32]; extern double masterTime; extern FILE *INFILE, *permutationFile, *logFile, *infoFile; extern char seq_file[1024]; extern char permFileName[1024], resultFileName[1024], logFileName[1024], checkpointFileName[1024], infoFileName[1024], run_id[128], workdir[1024], bootStrapFile[1024], bootstrapFileName[1024], bipartitionsFileName[1024],bipartitionsFileNameBranchLabels[1024]; void catToGamma(tree *tr, analdef *adef) { assert(tr->rateHetModel == CAT); if(adef->useInvariant) tr->rateHetModel = GAMMA_I; else tr->rateHetModel = GAMMA; #ifdef _USE_PTHREADS masterBarrier(THREAD_CAT_TO_GAMMA, tr); #endif } void gammaToCat(tree *tr) { assert(tr->rateHetModel == GAMMA || tr->rateHetModel == GAMMA_I); tr->rateHetModel = CAT; #ifdef _USE_PTHREADS masterBarrier(THREAD_GAMMA_TO_CAT, tr); #endif } static void singleBootstrap(tree *tr, int i, analdef *adef, rawdata *rdta, cruncheddata *cdta) { tr->treeID = i; tr->checkPointCounter = 0; computeNextReplicate(tr, &adef->boot, (int*)NULL, (int*)NULL, FALSE, FALSE); initModel(tr, rdta, cdta, adef); getStartingTree(tr, adef); computeBIGRAPID(tr, adef, TRUE); if(adef->bootstrapBranchLengths) { switch(tr->rateHetModel) { case GAMMA: case GAMMA_I: modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); break; case CAT: tr->likelihood = unlikely; catToGamma(tr, adef); initModel(tr, rdta, cdta, adef); if(i == 0) modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); else modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); gammaToCat(tr); break; default: assert(0); } } printBootstrapResult(tr, adef, TRUE); } /***************************** EXPERIMENTAL FUNCTIONS ********************************************************/ static int compareTopolRell(const void *p1, const void *p2) { topolRELL **rc1 = (topolRELL **)p1; topolRELL **rc2 = (topolRELL **)p2; double i = (*rc1)->likelihood; double j = (*rc2)->likelihood; if (i > j) return (-1); if (i < j) return (1); return (0); } void fixModelIndices(tree *tr, int endsite, boolean fixRates) { int model, i; assert(tr->NumberOfModels > 0); tr->partitionData[0].lower = 0; model = tr->model[0]; i = 1; while(i < endsite) { if(tr->model[i] != model) { tr->partitionData[model].upper = i; tr->partitionData[model + 1].lower = i; model = tr->model[i]; } i++; } tr->partitionData[tr->NumberOfModels - 1].upper = endsite; for(model = 0; model < tr->NumberOfModels; model++) tr->partitionData[model].width = tr->partitionData[model].upper - tr->partitionData[model].lower; #ifndef _USE_PTHREADS for(model = 0; model < tr->NumberOfModels; model++) { int j, lower = tr->partitionData[model].lower; /* SOS what about sumBuffer? */ /* tr->partitionData[model].sumBuffer = &tr->sumBuffer[offset]; */ tr->partitionData[model].perSiteLL = &tr->perSiteLL[lower]; tr->partitionData[model].wgt = &tr->cdta->aliaswgt[lower]; tr->partitionData[model].invariant = &tr->invariant[lower]; tr->partitionData[model].rateCategory = &tr->cdta->rateCategory[lower]; for(j = 1; j <= tr->mxtips; j++) tr->partitionData[model].yVector[j] = &(tr->yVector[j][tr->partitionData[model].lower]); { int width = tr->partitionData[model].width, undetermined = getUndetermined(tr->partitionData[model].dataType), j; tr->partitionData[model].gapVectorLength = ((int)width / 32) + 1; memset(tr->partitionData[model].gapVector, 0, tr->partitionData[model].initialGapVectorSize); for(j = 1; j <= tr->mxtips; j++) for(i = 0; i < width; i++) if(tr->partitionData[model].yVector[j][i] == undetermined) tr->partitionData[model].gapVector[tr->partitionData[model].gapVectorLength * j + i / 32] |= mask32[i % 32]; } } #else masterBarrier(THREAD_FIX_MODEL_INDICES, tr); #endif if(fixRates) updatePerSiteRates(tr, TRUE); } void reductionCleanup(tree *tr, int *originalRateCategories, int *originalInvariant) { tr->cdta->endsite = tr->originalCrunchedLength; memcpy(tr->cdta->aliaswgt, tr->originalWeights, sizeof(int) * tr->cdta->endsite); memcpy(tr->model, tr->originalModel, sizeof(int) * tr->cdta->endsite); memcpy(tr->dataVector, tr->originalDataVector, sizeof(int) * tr->cdta->endsite); memcpy(tr->cdta->rateCategory, originalRateCategories, sizeof(int) * tr->cdta->endsite); memcpy(tr->invariant, originalInvariant, sizeof(int) * tr->cdta->endsite); memcpy(tr->rdta->y0, tr->rdta->yBUF, ((size_t)tr->rdta->numsp) * ((size_t)tr->cdta->endsite) * sizeof(char)); tr->cdta->endsite = tr->originalCrunchedLength; fixModelIndices(tr, tr->originalCrunchedLength, TRUE); } /* old problematic code by Andre void computeNextReplicate(tree *tr, long *randomSeed, int *originalRateCategories, int *originalInvariant, boolean isRapid, boolean fixRates) { int j, count = 0, model, *weightBuffer, endsite, *weights, i, l; for(j = 0; j < tr->originalCrunchedLength; j++) tr->cdta->aliaswgt[j] = 0; for(model = 0; model < tr->NumberOfModels; ++model) { int maxWeight = 0, ctr = 0, // totalWeight = 0, realNumSites = tr->origNumSitePerModel[model], *wgtVirtualAln = (int*) rax_calloc(realNumSites, sizeof(int)); for(j = 0; j < tr->originalCrunchedLength; ++j) { if(tr->originalModel[j] == model) { wgtVirtualAln[ctr++] = tr->originalWeights[j]; //totalWeight += tr->originalWeights[j] ; if(maxWeight < tr->originalWeights[j]) maxWeight = tr->originalWeights[j]; } } weightBuffer = (int *)rax_calloc(realNumSites, sizeof(int)); j = 0; while( j < realNumSites ) { int pos = (int) (realNumSites * randum(randomSeed)); //if((int) (totalWeight * randum(randomSeed)) < wgtVirtualAln[pos] ) if((int) (maxWeight * randum(randomSeed)) < wgtVirtualAln[pos] ) { weightBuffer[pos]++; j++; } } ctr = 0; for(j = 0; j < tr->originalCrunchedLength; ++j) { if(model == tr->originalModel[j]) { tr->cdta->aliaswgt[j] = weightBuffer[ctr]; ctr++; } } rax_free(weightBuffer); rax_free(wgtVirtualAln); } endsite = 0; for (j = 0; j < tr->originalCrunchedLength; j++) if(tr->cdta->aliaswgt[j] > 0) endsite++; weights = tr->cdta->aliaswgt; for(i = 0; i < tr->rdta->numsp; i++) { unsigned char *yPos = &(tr->rdta->y0[((size_t)tr->originalCrunchedLength) * ((size_t)i)]), *origSeq = &(tr->rdta->yBUF[((size_t)tr->originalCrunchedLength) * ((size_t)i)]); int l, j; for(j = 0, l = 0; j < tr->originalCrunchedLength; j++) if(tr->cdta->aliaswgt[j] > 0) yPos[l++] = origSeq[j]; } for(j = 0, l = 0; j < tr->originalCrunchedLength; j++) { if(weights[j]) { tr->cdta->aliaswgt[l] = tr->cdta->aliaswgt[j]; tr->dataVector[l] = tr->originalDataVector[j]; tr->model[l] = tr->originalModel[j]; if(isRapid) { tr->cdta->rateCategory[l] = originalRateCategories[j]; tr->invariant[l] = originalInvariant[j]; } l++; } } tr->cdta->endsite = endsite; fixModelIndices(tr, endsite, fixRates); count = 0; for(j = 0; j < tr->cdta->endsite; j++) count += tr->cdta->aliaswgt[j]; if(count != tr->rdta->sites) printf("count=%d\ttr->rdta->sites=%d\n",count, tr->rdta->sites ); assert(count == tr->rdta->sites); } */ void computeNextReplicate(tree *tr, long *randomSeed, int *originalRateCategories, int *originalInvariant, boolean isRapid, boolean fixRates) { int j, model, w, *weightBuffer, endsite, *weights, i, l; for(j = 0; j < tr->originalCrunchedLength; j++) tr->cdta->aliaswgt[j] = 0; for(model = 0; model < tr->NumberOfModels; model++) { int nonzero = 0, pos = 0; for (j = 0; j < tr->originalCrunchedLength; j++) { if(tr->originalModel[j] == model) nonzero += tr->originalWeights[j]; } weightBuffer = (int *)rax_calloc(nonzero, sizeof(int)); for (j = 0; j < nonzero; j++) weightBuffer[(int) (nonzero*randum(randomSeed))]++; for(j = 0; j < tr->originalCrunchedLength; j++) { if(model == tr->originalModel[j]) { for(w = 0; w < tr->originalWeights[j]; w++) { tr->cdta->aliaswgt[j] += weightBuffer[pos]; pos++; } } } rax_free(weightBuffer); } endsite = 0; for (j = 0; j < tr->originalCrunchedLength; j++) { if(tr->cdta->aliaswgt[j] > 0) endsite++; } weights = tr->cdta->aliaswgt; for(i = 0; i < tr->rdta->numsp; i++) { unsigned char *yPos = &(tr->rdta->y0[((size_t)tr->originalCrunchedLength) * ((size_t)i)]), *origSeq = &(tr->rdta->yBUF[((size_t)tr->originalCrunchedLength) * ((size_t)i)]); for(j = 0, l = 0; j < tr->originalCrunchedLength; j++) if(tr->cdta->aliaswgt[j] > 0) yPos[l++] = origSeq[j]; } for(j = 0, l = 0; j < tr->originalCrunchedLength; j++) { if(weights[j]) { tr->cdta->aliaswgt[l] = tr->cdta->aliaswgt[j]; tr->dataVector[l] = tr->originalDataVector[j]; tr->model[l] = tr->originalModel[j]; if(isRapid) { tr->cdta->rateCategory[l] = originalRateCategories[j]; tr->invariant[l] = originalInvariant[j]; } l++; } } tr->cdta->endsite = endsite; fixModelIndices(tr, endsite, fixRates); { int count = 0; for(j = 0; j < tr->cdta->endsite; j++) count += tr->cdta->aliaswgt[j]; if(count != tr->rdta->sites) printf("count=%d\ttr->rdta->sites=%d\n",count, tr->rdta->sites ); assert(count == tr->rdta->sites); } } static pInfo *allocParams(tree *tr) { int i; pInfo *partBuffer = (pInfo*)rax_malloc(sizeof(pInfo) * tr->NumberOfModels); for(i = 0; i < tr->NumberOfModels; i++) { const partitionLengths *pl = getPartitionLengths(&(tr->partitionData[i])); partBuffer[i].EIGN = (double*)rax_malloc(pl->eignLength * sizeof(double)); partBuffer[i].EV = (double*)rax_malloc(pl->evLength * sizeof(double)); partBuffer[i].EI = (double*)rax_malloc(pl->eiLength * sizeof(double)); partBuffer[i].substRates = (double *)rax_malloc(pl->substRatesLength * sizeof(double)); partBuffer[i].frequencies = (double*)rax_malloc(pl->frequenciesLength * sizeof(double)); partBuffer[i].tipVector = (double *)rax_malloc(pl->tipVectorLength * sizeof(double)); } return partBuffer; } static void rax_freeParams(int numberOfModels, pInfo *partBuffer) { int i; for(i = 0; i < numberOfModels; i++) { rax_free(partBuffer[i].EIGN); rax_free(partBuffer[i].EV); rax_free(partBuffer[i].EI); rax_free(partBuffer[i].substRates); rax_free(partBuffer[i].frequencies); rax_free(partBuffer[i].tipVector); } } static void copyParams(int numberOfModels, pInfo *dst, pInfo *src, tree *tr) { int i; assert(src != dst); for(i = 0; i < numberOfModels; i++) { const partitionLengths *pl = getPartitionLengths(&src[i]); dst[i].dataType = src[i].dataType; memcpy(dst[i].EIGN, src[i].EIGN, pl->eignLength * sizeof(double)); memcpy(dst[i].EV, src[i].EV, pl->evLength * sizeof(double)); memcpy(dst[i].EI, src[i].EI, pl->eiLength * sizeof(double)); memcpy(dst[i].substRates, src[i].substRates, pl->substRatesLength * sizeof(double)); memcpy(dst[i].frequencies, src[i].frequencies, pl->frequenciesLength * sizeof(double)); memcpy(dst[i].tipVector, src[i].tipVector, pl->tipVectorLength * sizeof(double)); } #ifdef _USE_PTHREADS masterBarrier(THREAD_COPY_PARAMS, tr); #endif } #ifdef _WAYNE_MPI static void copyFiles(FILE *source, FILE *destination) { size_t c; char buf[1024]; assert(processID == 0); while((c = fread(buf, sizeof(char), 1024, source)) > 0) fwrite(buf, sizeof(char), c, destination); } static void concatenateBSFiles(int processes) { if(processID == 0) { int i; FILE *destination = myfopen(bootstrapFileName, "w"), *source = (FILE*)NULL; char sourceName[1024]; strcpy(sourceName, bootstrapFileName); strcat(sourceName, ".PID."); for(i = 0; i < processes; i++) { char buf[64], temporary[1024]; sprintf(buf, "%d", i); strcpy(temporary, sourceName); strcat(temporary, buf); source = myfopen(temporary, "r"); copyFiles(source, destination); fclose(source); } fclose(destination); } } static void removeBSFiles(int processes) { if(processID == 0) { int i; char sourceName[1024]; strcpy(sourceName, bootstrapFileName); strcat(sourceName, ".PID."); for(i = 0; i < processes; i++) { char buf[64], temporary[1024]; sprintf(buf, "%d", i); strcpy(temporary, sourceName); strcat(temporary, buf); remove(temporary); } } } #endif void doAllInOne(tree *tr, analdef *adef) { int i, n, bestIndex, bootstrapsPerformed; #ifdef _WAYNE_MPI int bootStopTests = 1, j, bootStrapsPerProcess = 0; #endif double loopTime; int *originalRateCategories; int *originalInvariant; #ifdef _WAYNE_MPI int slowSearches, fastEvery; #else int slowSearches, fastEvery = 5; #endif int treeVectorLength = -1; topolRELL_LIST *rl; double bestLH, mlTime, overallTime; long radiusSeed = adef->rapidBoot; FILE *f; char bestTreeFileName[1024]; hashtable *h = (hashtable*)NULL; unsigned int **bitVectors = (unsigned int**)NULL; boolean bootStopIt = FALSE; double pearsonAverage = 0.0; pInfo *catParams = allocParams(tr); pInfo *gammaParams = allocParams(tr); unsigned int vLength; n = adef->multipleRuns; #ifdef _WAYNE_MPI if(n % processes != 0) n = processes * ((n / processes) + 1); #endif if(adef->bootStopping) { h = initHashTable(tr->mxtips * 100); treeVectorLength = adef->multipleRuns; bitVectors = initBitVector(tr, &vLength); } rl = (topolRELL_LIST *)rax_malloc(sizeof(topolRELL_LIST)); initTL(rl, tr, n); originalRateCategories = (int*)rax_malloc(tr->cdta->endsite * sizeof(int)); originalInvariant = (int*)rax_malloc(tr->cdta->endsite * sizeof(int)); initModel(tr, tr->rdta, tr->cdta, adef); if(adef->grouping) printBothOpen("\n\nThe topologies of all Bootstrap and ML trees will adhere to the constraint tree specified in %s\n", tree_file); if(adef->constraint) printBothOpen("\n\nThe topologies of all Bootstrap and ML trees will adhere to the bifurcating backbone constraint tree specified in %s\n", tree_file); #ifdef _WAYNE_MPI long parsimonySeed0 = adef->parsimonySeed; long replicateSeed0 = adef->rapidBoot; n = n / processes; #endif for(i = 0; i < n && !bootStopIt; i++) { #ifdef _WAYNE_MPI j = i + n * processID; tr->treeID = j; #else tr->treeID = i; #endif tr->checkPointCounter = 0; loopTime = gettime(); #ifdef _WAYNE_MPI if(i == 0) { if(parsimonySeed0 != 0) adef->parsimonySeed = parsimonySeed0 + 10000 * processID; adef->rapidBoot = replicateSeed0 + 10000 * processID; radiusSeed = adef->rapidBoot; } #endif if(i % 10 == 0) { if(i > 0) reductionCleanup(tr, originalRateCategories, originalInvariant); if(adef->grouping || adef->constraint) { FILE *f = myfopen(tree_file, "rb"); assert(adef->restart); partCount = 0; if (! treeReadLenMULT(f, tr, adef)) exit(-1); fclose(f); } else makeParsimonyTree(tr, adef); tr->likelihood = unlikely; if(i == 0) { double t; onlyInitrav(tr, tr->start); treeEvaluate(tr, 1); t = gettime(); modOpt(tr, adef, FALSE, 5.0); #ifdef _WAYNE_MPI printBothOpen("\nTime for BS model parameter optimization on Process %d: %f seconds\n", processID, gettime() - t); #else printBothOpen("\nTime for BS model parameter optimization %f\n", gettime() - t); #endif memcpy(originalRateCategories, tr->cdta->rateCategory, sizeof(int) * tr->cdta->endsite); memcpy(originalInvariant, tr->invariant, sizeof(int) * tr->cdta->endsite); if(adef->bootstrapBranchLengths) { if(tr->rateHetModel == CAT) { copyParams(tr->NumberOfModels, catParams, tr->partitionData, tr); assert(tr->cdta->endsite == tr->originalCrunchedLength); catToGamma(tr, adef); modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); copyParams(tr->NumberOfModels, gammaParams, tr->partitionData, tr); gammaToCat(tr); copyParams(tr->NumberOfModels, tr->partitionData, catParams, tr); } else { assert(tr->cdta->endsite == tr->originalCrunchedLength); } } } } computeNextReplicate(tr, &adef->rapidBoot, originalRateCategories, originalInvariant, TRUE, TRUE); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 1); computeBOOTRAPID(tr, adef, &radiusSeed); #ifdef _WAYNE_MPI saveTL(rl, tr, j); #else saveTL(rl, tr, i); #endif if(adef->bootstrapBranchLengths) { double lh = tr->likelihood; if(tr->rateHetModel == CAT) { copyParams(tr->NumberOfModels, tr->partitionData, gammaParams, tr); catToGamma(tr, adef); resetBranches(tr); onlyInitrav(tr, tr->start); treeEvaluate(tr, 2.0); gammaToCat(tr); copyParams(tr->NumberOfModels, tr->partitionData, catParams, tr); tr->likelihood = lh; } else { treeEvaluate(tr, 2.0); tr->likelihood = lh; } } printBootstrapResult(tr, adef, TRUE); loopTime = gettime() - loopTime; writeInfoFile(adef, tr, loopTime); if(adef->bootStopping) #ifdef _WAYNE_MPI { int nn = (i + 1) * processes; if((nn > START_BSTOP_TEST) && (i * processes < FC_SPACING * bootStopTests) && ((i + 1) * processes >= FC_SPACING * bootStopTests) ) { MPI_Barrier(MPI_COMM_WORLD); concatenateBSFiles(processes); MPI_Barrier(MPI_COMM_WORLD); bootStopIt = computeBootStopMPI(tr, bootstrapFileName, adef, &pearsonAverage); bootStopTests++; } } #else bootStopIt = bootStop(tr, h, i, &pearsonAverage, bitVectors, treeVectorLength, vLength); #endif } #ifdef _WAYNE_MPI MPI_Barrier(MPI_COMM_WORLD); bootstrapsPerformed = i * processes; bootStrapsPerProcess = i; concatenateBSFiles(processes); removeBSFiles(processes); MPI_Barrier(MPI_COMM_WORLD); #else bootstrapsPerformed = i; #endif rax_freeParams(tr->NumberOfModels, catParams); rax_free(catParams); rax_freeParams(tr->NumberOfModels, gammaParams); rax_free(gammaParams); if(adef->bootStopping) { freeBitVectors(bitVectors, 2 * tr->mxtips); rax_free(bitVectors); freeHashTable(h); rax_free(h); } { double t; printBothOpenMPI("\n\n"); if(adef->bootStopping) { if(bootStopIt) { switch(tr->bootStopCriterion) { case FREQUENCY_STOP: printBothOpenMPI("Stopped Rapid BS search after %d replicates with FC Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("Pearson Average of %d random splits: %f\n",BOOTSTOP_PERMUTATIONS , pearsonAverage); break; case MR_STOP: printBothOpenMPI("Stopped Rapid BS search after %d replicates with MR-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; case MRE_STOP: printBothOpenMPI("Stopped Rapid BS search after %d replicates with MRE-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; case MRE_IGN_STOP: printBothOpenMPI("Stopped Rapid BS search after %d replicates with MRE_IGN-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; default: assert(0); } } else { switch(tr->bootStopCriterion) { case FREQUENCY_STOP: printBothOpenMPI("Rapid BS search did not converge after %d replicates with FC Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("Pearson Average of %d random splits: %f\n",BOOTSTOP_PERMUTATIONS , pearsonAverage); break; case MR_STOP: printBothOpenMPI("Rapid BS search did not converge after %d replicates with MR-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; case MRE_STOP: printBothOpenMPI("Rapid BS search did not converge after %d replicates with MRE-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; case MRE_IGN_STOP: printBothOpenMPI("Rapid BS search did not converge after %d replicates with MR_IGN-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; default: assert(0); } } } t = gettime() - masterTime; printBothOpenMPI("Overall Time for %d Rapid Bootstraps %f seconds\n", bootstrapsPerformed, t); printBothOpenMPI("Average Time per Rapid Bootstrap %f seconds\n", (double)(t/((double)bootstrapsPerformed))); if(!adef->allInOne) { printBothOpenMPI("All %d bootstrapped trees written to: %s\n", bootstrapsPerformed, bootstrapFileName); #ifdef _WAYNE_MPI MPI_Finalize(); #endif exit(0); } } /* ML-search */ mlTime = gettime(); double t = mlTime; printBothOpenMPI("\nStarting ML Search ...\n\n"); /***CLEAN UP reduction stuff */ reductionCleanup(tr, originalRateCategories, originalInvariant); /****/ #ifdef _WAYNE_MPI restoreTL(rl, tr, n * processID); #else restoreTL(rl, tr, 0); #endif resetBranches(tr); evaluateGenericInitrav(tr, tr->start); modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); #ifdef _WAYNE_MPI if(bootstrapsPerformed <= 100) fastEvery = 5; else fastEvery = bootstrapsPerformed / 20; for(i = 0; i < bootstrapsPerformed; i++) rl->t[i]->likelihood = unlikely; for(i = 0; i < bootStrapsPerProcess; i++) { j = i + n * processID; if(i % fastEvery == 0) { restoreTL(rl, tr, j); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 1); optimizeRAPID(tr, adef); saveTL(rl, tr, j); } } #else for(i = 0; i < bootstrapsPerformed; i++) { rl->t[i]->likelihood = unlikely; if(i % fastEvery == 0) { restoreTL(rl, tr, i); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 1); optimizeRAPID(tr, adef); saveTL(rl, tr, i); } } #endif printBothOpenMPI("Fast ML optimization finished\n\n"); t = gettime() - t; #ifdef _WAYNE_MPI printBothOpen("Fast ML search on Process %d: Time %f seconds\n\n", processID, t); j = n * processID; qsort(&(rl->t[j]), n, sizeof(topolRELL*), compareTopolRell); restoreTL(rl, tr, j); #else printBothOpen("Fast ML search Time: %f seconds\n\n", t); qsort(&(rl->t[0]), bootstrapsPerformed, sizeof(topolRELL*), compareTopolRell); restoreTL(rl, tr, 0); #endif t = gettime(); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); slowSearches = bootstrapsPerformed / 5; if(bootstrapsPerformed % 5 != 0) slowSearches++; slowSearches = MIN(slowSearches, 10); #ifdef _WAYNE_MPI if(processes > 1) { if(slowSearches % processes == 0) slowSearches = slowSearches / processes; else slowSearches = (slowSearches / processes) + 1; } for(i = 0; i < slowSearches; i++) { j = i + n * processID; restoreTL(rl, tr, j); rl->t[j]->likelihood = unlikely; evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 1.0); thoroughOptimization(tr, adef, rl, j); } #else for(i = 0; i < slowSearches; i++) { restoreTL(rl, tr, i); rl->t[i]->likelihood = unlikely; evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 1.0); thoroughOptimization(tr, adef, rl, i); } #endif /*************************************************************************************************************/ if(tr->rateHetModel == CAT) { catToGamma(tr, adef); modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); } bestIndex = -1; bestLH = unlikely; #ifdef _WAYNE_MPI for(i = 0; i < slowSearches; i++) { j = i + n * processID; restoreTL(rl, tr, j); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 2); printBothOpen("Slow ML Search %d Likelihood: %f\n", j, tr->likelihood); if(tr->likelihood > bestLH) { bestLH = tr->likelihood; bestIndex = j; } } /*printf("processID = %d, bestIndex = %d; bestLH = %f\n", processID, bestIndex, bestLH);*/ #else for(i = 0; i < slowSearches; i++) { restoreTL(rl, tr, i); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 2); printBothOpen("Slow ML Search %d Likelihood: %f\n", i, tr->likelihood); if(tr->likelihood > bestLH) { bestLH = tr->likelihood; bestIndex = i; } } #endif printBothOpenMPI("Slow ML optimization finished\n\n"); t = gettime() - t; #ifdef _WAYNE_MPI printBothOpen("Slow ML search on Process %d: Time %f seconds\n", processID, t); #else printBothOpen("Slow ML search Time: %f seconds\n", t); #endif t = gettime(); restoreTL(rl, tr, bestIndex); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 2); Thorough = 1; tr->doCutoff = FALSE; treeOptimizeThorough(tr, 1, 10); evaluateGenericInitrav(tr, tr->start); modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); t = gettime() - t; #ifdef _WAYNE_MPI printBothOpen("Thorough ML search on Process %d: Time %f seconds\n", processID, t); #else printBothOpen("Thorough ML search Time: %f seconds\n", t); #endif #ifdef _WAYNE_MPI bestLH = tr->likelihood; printf("\nprocessID = %d, bestLH = %f\n", processID, bestLH); if(processes > 1) { double *buffer; int bestProcess; buffer = (double *)rax_malloc(sizeof(double) * processes); for(i = 0; i < processes; i++) buffer[i] = unlikely; buffer[processID] = bestLH; for(i = 0; i < processes; i++) MPI_Bcast(&buffer[i], 1, MPI_DOUBLE, i, MPI_COMM_WORLD); bestLH = buffer[0]; bestProcess = 0; for(i = 1; i < processes; i++) if(buffer[i] > bestLH) { bestLH = buffer[i]; bestProcess = i; } rax_free(buffer); if(processID != bestProcess) { MPI_Finalize(); exit(0); } } #endif printBothOpen("\nFinal ML Optimization Likelihood: %f\n", tr->likelihood); printBothOpen("\nModel Information:\n\n"); printModelParams(tr, adef); strcpy(bestTreeFileName, workdir); 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, FALSE); f = myfopen(bestTreeFileName, "wb"); fprintf(f, "%s", tr->tree_string); fclose(f); if(adef->perGeneBranchLengths) printTreePerGene(tr, adef, bestTreeFileName, "w"); overallTime = gettime() - masterTime; mlTime = gettime() - mlTime; printBothOpen("\nML search took %f secs or %f hours\n", mlTime, mlTime / 3600.0); printBothOpen("\nCombined Bootstrap and ML search took %f secs or %f hours\n", overallTime, overallTime / 3600.0); printBothOpen("\nDrawing Bootstrap Support Values on best-scoring ML tree ...\n\n"); freeTL(rl); rax_free(rl); calcBipartitions(tr, adef, bestTreeFileName, bootstrapFileName); overallTime = gettime() - masterTime; printBothOpen("Program execution info written to %s\n", infoFileName); printBothOpen("All %d bootstrapped trees written to: %s\n\n", bootstrapsPerformed, bootstrapFileName); printBothOpen("Best-scoring ML tree written to: %s\n\n", bestTreeFileName); if(adef->perGeneBranchLengths && tr->NumberOfModels > 1) printBothOpen("Per-Partition branch lengths of best-scoring ML tree written to %s.PARTITION.0 to %s.PARTITION.%d\n\n", bestTreeFileName, bestTreeFileName, tr->NumberOfModels - 1); printBothOpen("Best-scoring ML tree with support values written to: %s\n\n", bipartitionsFileName); printBothOpen("Best-scoring ML tree with support values as branch labels written to: %s\n\n", bipartitionsFileNameBranchLabels); printBothOpen("Overall execution time for full ML analysis: %f secs or %f hours or %f days\n\n", overallTime, overallTime/3600.0, overallTime/86400.0); #ifdef _WAYNE_MPI MPI_Finalize(); #endif exit(0); } /*******************************************EXPERIMENTAL FUNCTIONS END *****************************************************/ void doBootstrap(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta) { int bootstrapsPerformed, i, n, treeVectorLength = -1; unsigned int vLength = 0; #ifdef _WAYNE_MPI int j, bootStopTests = 1; #endif double loopTime, pearsonAverage; hashtable *h = (hashtable*)NULL; unsigned int **bitVectors = (unsigned int **)NULL; boolean bootStopIt = FALSE; n = adef->multipleRuns; #ifdef _WAYNE_MPI if(n % processes != 0) n = processes * ((n / processes) + 1); adef->multipleRuns = n; #endif if(adef->bootStopping) { h = initHashTable(tr->mxtips * 100); bitVectors = initBitVector(tr, &vLength); treeVectorLength = adef->multipleRuns; } #ifdef _WAYNE_MPI long parsimonySeed0 = adef->parsimonySeed; long replicateSeed0 = adef->rapidBoot; long bootstrapSeed0 = adef->boot; n = n / processes; #endif for(i = 0; i < n && !bootStopIt; i++) { loopTime = gettime(); #ifdef _WAYNE_MPI if(i == 0) { if(parsimonySeed0 != 0) adef->parsimonySeed = parsimonySeed0 + 10000 * processID; adef->rapidBoot = replicateSeed0 + 10000 * processID; adef->boot = bootstrapSeed0 + 10000 * processID; } j = i + n*processID; singleBootstrap(tr, j, adef, rdta, cdta); #else singleBootstrap(tr, i, adef, rdta, cdta); #endif loopTime = gettime() - loopTime; writeInfoFile(adef, tr, loopTime); if(adef->bootStopping) #ifdef _WAYNE_MPI { int nn = (i + 1) * processes; if((nn > START_BSTOP_TEST) && (i * processes < FC_SPACING * bootStopTests) && ((i + 1) * processes >= FC_SPACING * bootStopTests) ) { MPI_Barrier(MPI_COMM_WORLD); concatenateBSFiles(processes); MPI_Barrier(MPI_COMM_WORLD); bootStopIt = computeBootStopMPI(tr, bootstrapFileName, adef, &pearsonAverage); bootStopTests++; } } #else bootStopIt = bootStop(tr, h, i, &pearsonAverage, bitVectors, treeVectorLength, vLength); #endif } #ifdef _WAYNE_MPI MPI_Barrier(MPI_COMM_WORLD); bootstrapsPerformed = i * processes; if(processID == 0) { if(!adef->bootStopping) concatenateBSFiles(processes); removeBSFiles(processes); } MPI_Barrier(MPI_COMM_WORLD); #else bootstrapsPerformed = i; #endif adef->multipleRuns = bootstrapsPerformed; if(adef->bootStopping) { freeBitVectors(bitVectors, 2 * tr->mxtips); rax_free(bitVectors); freeHashTable(h); rax_free(h); if(bootStopIt) { switch(tr->bootStopCriterion) { case FREQUENCY_STOP: printBothOpenMPI("Stopped Standard BS search after %d replicates with FC Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("Pearson Average of %d random splits: %f\n",BOOTSTOP_PERMUTATIONS , pearsonAverage); break; case MR_STOP: printBothOpenMPI("Stopped Standard BS search after %d replicates with MR-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; case MRE_STOP: printBothOpenMPI("Stopped Standard BS search after %d replicates with MRE-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; case MRE_IGN_STOP: printBothOpenMPI("Stopped Standard BS search after %d replicates with MRE_IGN-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; default: assert(0); } } else { switch(tr->bootStopCriterion) { case FREQUENCY_STOP: printBothOpenMPI("Standard BS search did not converge after %d replicates with FC Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("Pearson Average of %d random splits: %f\n",BOOTSTOP_PERMUTATIONS , pearsonAverage); break; case MR_STOP: printBothOpenMPI("Standard BS search did not converge after %d replicates with MR-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; case MRE_STOP: printBothOpenMPI("Standard BS search did not converge after %d replicates with MRE-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; case MRE_IGN_STOP: printBothOpenMPI("Standard BS search did not converge after %d replicates with MR_IGN-based Bootstopping criterion\n", bootstrapsPerformed); printBothOpenMPI("WRF Average of %d random splits: %f\n", BOOTSTOP_PERMUTATIONS, pearsonAverage); break; default: assert(0); } } } } void doInference(tree *tr, analdef *adef, rawdata *rdta, cruncheddata *cdta) { int i, n; #ifdef _WAYNE_MPI int j; #endif double loopTime; topolRELL_LIST *rl = (topolRELL_LIST *)NULL; int best = -1, newBest = -1; double bestLH = unlikely; FILE *f; char bestTreeFileName[1024]; double overallTime; n = adef->multipleRuns; #ifdef _WAYNE_MPI if(n % processes != 0) n = processes * ((n / processes) + 1); #endif if(!tr->catOnly) { rl = (topolRELL_LIST *)rax_malloc(sizeof(topolRELL_LIST)); initTL(rl, tr, n); } #ifdef _WAYNE_MPI long parsimonySeed0 = adef->parsimonySeed; n = n / processes; #endif for(i = 0; i < n; i++) { #ifdef _WAYNE_MPI if(i == 0) { if(parsimonySeed0 != 0) adef->parsimonySeed = parsimonySeed0 + 10000 * processID; } j = i + n * processID; tr->treeID = j; #else tr->treeID = i; #endif tr->checkPointCounter = 0; loopTime = gettime(); initModel(tr, rdta, cdta, adef); if(i == 0) printBaseFrequencies(tr); getStartingTree(tr, adef); computeBIGRAPID(tr, adef, TRUE); #ifdef _WAYNE_MPI if(tr->likelihood > bestLH) { best = j; bestLH = tr->likelihood; } if(!tr->catOnly) saveTL(rl, tr, j); #else if(tr->likelihood > bestLH) { best = i; bestLH = tr->likelihood; } if(!tr->catOnly) saveTL(rl, tr, i); #endif loopTime = gettime() - loopTime; writeInfoFile(adef, tr, loopTime); } assert(best >= 0); #ifdef _WAYNE_MPI MPI_Barrier(MPI_COMM_WORLD); n = n * processes; #endif if(tr->catOnly) { printBothOpenMPI("\n\nNOT conducting any final model optimizations on all %d trees under CAT-based model ....\n", n); printBothOpenMPI("\nREMEMBER that CAT-based likelihood scores are meaningless!\n\n", n); #ifdef _WAYNE_MPI if(processID != 0) { MPI_Finalize(); exit(0); } #endif } else { printBothOpenMPI("\n\nConducting final model optimizations on all %d trees under GAMMA-based models ....\n\n", n); #ifdef _WAYNE_MPI n = n / processes; #endif if(tr->rateHetModel == GAMMA || tr->rateHetModel == GAMMA_I) { restoreTL(rl, tr, best); evaluateGenericInitrav(tr, tr->start); if(!adef->useBinaryModelFile) modOpt(tr, adef, FALSE, adef->likelihoodEpsilon); else { readBinaryModel(tr); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 2); } bestLH = tr->likelihood; tr->likelihoods[best] = tr->likelihood; saveTL(rl, tr, best); tr->treeID = best; printResult(tr, adef, TRUE); newBest = best; for(i = 0; i < n; i++) { #ifdef _WAYNE_MPI j = i + n * processID; if(j != best) { restoreTL(rl, tr, j); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 1); tr->likelihoods[j] = tr->likelihood; if(tr->likelihood > bestLH) { newBest = j; bestLH = tr->likelihood; saveTL(rl, tr, j); } tr->treeID = j; printResult(tr, adef, TRUE); } if(n == 1 && processes == 1) printBothOpen("Inference[%d] final GAMMA-based Likelihood: %f tree written to file %s\n", i, tr->likelihoods[i], resultFileName); else printBothOpen("Inference[%d] final GAMMA-based Likelihood: %f tree written to file %s.RUN.%d\n", j, tr->likelihoods[j], resultFileName, j); #else if(i != best) { restoreTL(rl, tr, i); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 1); tr->likelihoods[i] = tr->likelihood; if(tr->likelihood > bestLH) { newBest = i; bestLH = tr->likelihood; saveTL(rl, tr, i); } tr->treeID = i; printResult(tr, adef, TRUE); } if(n == 1) printBothOpen("Inference[%d] final GAMMA-based Likelihood: %f tree written to file %s\n", i, tr->likelihoods[i], resultFileName); else printBothOpen("Inference[%d] final GAMMA-based Likelihood: %f tree written to file %s.RUN.%d\n", i, tr->likelihoods[i], resultFileName, i); #endif } } else { catToGamma(tr, adef); #ifdef _WAYNE_MPI for(i = 0; i < n; i++) { j = i + n*processID; rl->t[j]->likelihood = unlikely; } #else for(i = 0; i < n; i++) rl->t[i]->likelihood = unlikely; #endif initModel(tr, rdta, cdta, adef); restoreTL(rl, tr, best); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); modOpt(tr, adef, TRUE, adef->likelihoodEpsilon); tr->likelihoods[best] = tr->likelihood; bestLH = tr->likelihood; saveTL(rl, tr, best); tr->treeID = best; printResult(tr, adef, TRUE); newBest = best; for(i = 0; i < n; i++) { #ifdef _WAYNE_MPI j = i + n*processID; if(j != best) { restoreTL(rl, tr, j); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 2); tr->likelihoods[j] = tr->likelihood; if(tr->likelihood > bestLH) { newBest = j; bestLH = tr->likelihood; saveTL(rl, tr, j); } tr->treeID = j; printResult(tr, adef, TRUE); } if(n == 1 && processes == 1) printBothOpen("Inference[%d] final GAMMA-based Likelihood: %f tree written to file %s\n", i, tr->likelihoods[i], resultFileName); else printBothOpen("Inference[%d] final GAMMA-based Likelihood: %f tree written to file %s.RUN.%d\n", j, tr->likelihoods[j], resultFileName, j); #else if(i != best) { restoreTL(rl, tr, i); resetBranches(tr); evaluateGenericInitrav(tr, tr->start); treeEvaluate(tr, 2); tr->likelihoods[i] = tr->likelihood; if(tr->likelihood > bestLH) { newBest = i; bestLH = tr->likelihood; saveTL(rl, tr, i); } tr->treeID = i; printResult(tr, adef, TRUE); } if(n == 1) printBothOpen("Inference[%d] final GAMMA-based Likelihood: %f tree written to file %s\n", i, tr->likelihoods[i], resultFileName); else printBothOpen("Inference[%d] final GAMMA-based Likelihood: %f tree written to file %s.RUN.%d\n", i, tr->likelihoods[i], resultFileName, i); #endif } } assert(newBest >= 0); #ifdef _WAYNE_MPI if(processes > 1) { double *buffer; int bestProcess; buffer = (double *)rax_malloc(sizeof(double) * processes); for(i = 0; i < processes; i++) buffer[i] = unlikely; buffer[processID] = bestLH; for(i = 0; i < processes; i++) MPI_Bcast(&buffer[i], 1, MPI_DOUBLE, i, MPI_COMM_WORLD); bestLH = buffer[0]; bestProcess = 0; for(i = 1; i < processes; i++) if(buffer[i] > bestLH) { bestLH = buffer[i]; bestProcess = i; } rax_free(buffer); if(processID != bestProcess) { MPI_Finalize(); exit(0); } } #endif restoreTL(rl, tr, newBest); evaluateGenericInitrav(tr, tr->start); printBothOpen("\n\nStarting final GAMMA-based thorough Optimization on tree %d likelihood %f .... \n\n", newBest, tr->likelihoods[newBest]); Thorough = 1; tr->doCutoff = FALSE; treeOptimizeThorough(tr, 1, 10); evaluateGenericInitrav(tr, tr->start); printBothOpen("Final GAMMA-based Score of best tree %f\n\n", tr->likelihood); strcpy(bestTreeFileName, workdir); 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, FALSE); f = myfopen(bestTreeFileName, "wb"); fprintf(f, "%s", tr->tree_string); fclose(f); if(adef->perGeneBranchLengths) printTreePerGene(tr, adef, bestTreeFileName, "w"); } overallTime = gettime() - masterTime; printBothOpen("Program execution info written to %s\n", infoFileName); if(!tr->catOnly) { printBothOpen("Best-scoring ML tree written to: %s\n\n", bestTreeFileName); if(adef->perGeneBranchLengths && tr->NumberOfModels > 1) printBothOpen("Per-Partition branch lengths of best-scoring ML tree written to %s.PARTITION.0 to %s.PARTITION.%d\n\n", bestTreeFileName, bestTreeFileName, tr->NumberOfModels - 1); } printBothOpen("Overall execution time: %f secs or %f hours or %f days\n\n", overallTime, overallTime/3600.0, overallTime/86400.0); if(!tr->catOnly) { freeTL(rl); rax_free(rl); } #ifdef _WAYNE_MPI MPI_Finalize(); #endif exit(0); }