Skip to content

Commit

Permalink
finishing the command compareref
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhangchicool
zhangchicool committed Nov 22, 2015
1 parent a67c55e commit bc4cf96
Show file tree
Hide file tree
Showing 3 changed files with 315 additions and 47 deletions.
6 changes: 3 additions & 3 deletions src/command.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12947,9 +12947,9 @@ else if (!strcmp(helpTkn, "Set"))
MrBayesPrint (" tree distance values. \n");
MrBayesPrint (" \n");
MrBayesPrint (" Note that the \"Sumt\" command provides a different set of convergence diag- \n");
MrBayesPrint (" nostics tools that you may also want to explore. Unlike \"Comparetree\", \"Sumt\"\n");
MrBayesPrint (" can compare more than two tree samples and will calculate consensus trees and \n");
MrBayesPrint (" split frequencies from the pooled samples. \n");
MrBayesPrint (" nostics tools that you may also want to explore. Unlike \"Comparetree\", \n");
MrBayesPrint (" \"Sumt\" can compare more than two tree samples and will calculate consensus \n");
MrBayesPrint (" trees and split frequencies from the pooled samples. \n");
MrBayesPrint (" \n");
MrBayesPrint (" Options: \n");
MrBayesPrint (" \n");
Expand Down
90 changes: 70 additions & 20 deletions src/mcmc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -507,12 +507,12 @@ int AddTreeSamples (int from, int to, int saveToList)
{
do {
if (fgets (lineBuf, longestLine, fp) == NULL)
{
SafeFclose (&fp);
free (lineBuf);
free (tempStr);
{
SafeFclose (&fp);
free (lineBuf);
free (tempStr);
return ERROR;
}
}
word = strtok (lineBuf, " ");
} while (strcmp (word, "tree") != 0);
if (k>=from)
Expand Down Expand Up @@ -1708,7 +1708,7 @@ void CalcPartFreqStats (PFNODE *p, STATS *stat)

n = chainParams.numRuns;
min = (int)(chainParams.minPartFreq * stat->numSamples);
if (((MrBFlt)min) != chainParams.minPartFreq * stat->numSamples)
if ((MrBFlt)min != chainParams.minPartFreq * stat->numSamples)
min++;

/* recursively compute partition frequencies for all subpartitions */
Expand Down Expand Up @@ -1817,6 +1817,56 @@ void CalcTopoConvDiagn (int numSamples)
}


/* used in DoCompRefTree */
void PartFreq (PFNODE *p, STATS *stat, int *ntrees)
{
int i, n = chainParams.numRuns;
MrBFlt f, sum, sumsq, stdev;

/* recursively compute partition frequencies for all subpartitions */
if (p->left != NULL)
PartFreq (p->left, stat, ntrees);
if (p->right != NULL)
PartFreq (p->right, stat, ntrees);

sum = sumsq = 0.0;
for (i=0; i<chainParams.numRuns; i++)
{
f = (MrBFlt)(p->count[i]) / (MrBFlt)ntrees[i];
sum += f;
sumsq += f * f;
}

f = (sumsq - sum * sum / n) / (n - 1);
if (f < 0.0)
stdev = 0.0;
else
stdev = sqrt (f);

stat->sum += stdev;
if (stat->max < stdev)
stat->max = stdev;

stat->numPartitions++;
}
void CalcTopoConvDiagn2 (int *nTrees)
{
int n;
STATS *stat;

for (n=0; n<numTopologies; n++)
{
stat = &chainParams.stat[n];
stat->numPartitions = 0.0;
stat->sum = stat->max = 0.0;

PartFreq (partFreqTreeRoot[n], stat, nTrees);

stat->avgStdDev = stat->sum / stat->numPartitions;
}
}


int CheckTemperature (void)
{
if (chainParams.userDefinedTemps == YES)
Expand Down Expand Up @@ -16102,16 +16152,15 @@ int RunChain (RandLong *seed)
for (i=0; i<numTopologies; i++)
chainParams.stat[i].pair = NULL;
}
if (chainParams.relativeBurnin == YES)

if ((chainParams.dtree = AllocateTree (numLocalTaxa)) == NULL)
{
if ((chainParams.dtree = AllocateTree (numLocalTaxa)) == NULL)
{
MrBayesPrint ("%s Could not allocate chainParams.dtree in RunChain\n", spacer);
return ERROR;
}
else
memAllocs[ALLOC_DIAGNTREE] = YES;
MrBayesPrint ("%s Could not allocate chainParams.dtree in RunChain\n", spacer);
return ERROR;
}
else
memAllocs[ALLOC_DIAGNTREE] = YES;

if (chainParams.allComps == YES)
{
for (i=0; i<numTopologies; i++)
Expand Down Expand Up @@ -16871,12 +16920,13 @@ int RunChain (RandLong *seed)
}

/* print mcmc diagnostics. Blocking for MPI */
if (chainParams.mcmcDiagn == YES && (n % chainParams.diagnFreq == 0 || n == chainParams.numGen || (chainParams.isSS == YES && (n-lastStepEndSS) % numGenInStepSS == 0)))
if (chainParams.mcmcDiagn == YES && (n % chainParams.diagnFreq == 0
|| n == chainParams.numGen
|| (chainParams.isSS == YES && (n-lastStepEndSS) % numGenInStepSS == 0)))
{
if (chainParams.numRuns > 1
&& ((n > 0 && chainParams.relativeBurnin == YES &&
(chainParams.isSS == NO || (n > chainParams.burninSS * chainParams.sampleFreq && (n-lastStepEndSS) > numGenInStepBurninSS)))
|| (n >= chainParams.chainBurnIn * chainParams.sampleFreq && chainParams.relativeBurnin == NO)))
if (chainParams.numRuns > 1 &&
((n > 0 && chainParams.relativeBurnin == YES && (chainParams.isSS == NO || (n > chainParams.burninSS * chainParams.sampleFreq && (n-lastStepEndSS) > numGenInStepBurninSS)))
|| (n >= chainParams.chainBurnIn * chainParams.sampleFreq && chainParams.relativeBurnin == NO)))
{
/* we need some space for coming output */
MrBayesPrint ("\n");
Expand Down Expand Up @@ -16992,7 +17042,7 @@ int RunChain (RandLong *seed)
}
if (chainParams.isSS == YES)
splitfreqSS[i*chainParams.numStepsSS+chainParams.numStepsSS-stepIndexSS-1] = f;
if (chainParams.stat[i].numPartitions == 0 && f > chainParams.stopVal)
if (chainParams.stat[i].numPartitions == 0 || f > chainParams.stopVal)
stopChain = NO;
}
if (n < chainParams.numGen - chainParams.printFreq && (chainParams.stopRule == NO || stopChain == NO))
Expand Down
Loading

0 comments on commit bc4cf96

Please sign in to comment.