MrBayes
MrBayes copied to clipboard
Minimal example does not run
I had observed that MrBayes behaves vastly different from RevBayes and Beast2 for a fossilized birth death prior with sampled ancestors when set to nominally the same parameter values, so I tried to construct a minimal model using that code in the hopes to find the issue.
#NEXUS
BEGIN DATA;
DIMENSIONS NTAX=2 NCHAR=1;
FORMAT DATATYPE = restriction GAP = - MISSING = ? ;
MATRIX
A ?
B ?
;
END;
begin trees;
[ Set starting tree ]
Tree starting = (A:100,B:100);
end;
begin mrbayes;
set autoclose=no nowarn=yes;
prset brlenspr = clock:fossilization;
prset treeagepr = uniform(0, 500000000000) ;
prset speciationpr = fixed(0.1) ;
prset extinctionpr = fixed(0.0);
prset fossilizationpr = fixed(1.0);
prset sampleprob = 1.0;
prset samplestrat = random ;
prset nodeagepr = calibrated;
prset clockvarpr = igr;
prset igrvarpr = fixed(1.0);
prset statefreqpr = dirichlet(1, 1);
lset coding = NoAbsenceSites;
mcmcp ngen=12500000 Starttree=Current relburnin=yes burninfrac=0.2 printfreq=1000 swapfreq=50 samplefreq=1000 nchains=4 savebrlens=yes;
startvals tau=starting;
showmodel;
mcmc;
end;
What is the current observed behaviour?
Running the current version of MrBayes dies with the message
There must be at least two included taxa, now there is only one
This message is nonsense; Un-commenting (#if 0
→ #if 1
) the debug output clearly shows that there are two taxa, and that the issue is that in line 21238, if (numLocalTaxa <= 2)
includes equality.
Changing line 21238 to if (numLocalTaxa < 2)
leads to a SEGFAULT:
Program received signal SIGSEGV, Segmentation fault. 0x00005555556edc6b in GetNodeDownPass (t=0x5555557de800, p=0x5555557de968, i=0x7ffffffccc80, j=0x7ffffffccc84) at utils.c:3924 3924 t->allDownPass[(*j)++] = p;
The traceback is
#1703 0x00005555556edaf5 in GetDownPass (t=0x5555557de800) at utils.c:3899 #1704 0x00005555556ea314 in BuildRandomUTopology (t=0x5555557de800, seed=0x555555788de8 <globalSeed>) at utils.c:2663 #1705 0x000055555563cdeb in FillTreeParams (seed=0x555555788de8 <globalSeed>, fromChain=0, toChain=8) at model.c:12801 #1706 0x000055555566d8a9 in SetUpAnalysis (seed=0x555555788de8 <globalSeed>) at model.c:21300 #1707 0x000055555556dd03 in DoMatrix () at command.c:5223
#1708 0x000055555558c094 in ParseCommand (s=0x5555557d6550 ";\n") at command.c:14018
Note the call to BuildRandomUTopology, even though the model specifies a clock and thus a rooted tree.
I put a breakpoint on line 1036, which should lead to 1037 isRooted = YES;
. The breakpoint is not encountered by MrBayes before the segfault happens.
Dear @Anaphory ,
thank you so much for your report. Understanding if and why different software reports different results given the same implementation of a method is important.
However, and this is just my conjecture, your minimal example may not be feasible to analyse in MrBayes. And this because the fundamental design was aimed for analyzing biological data sets (with more than two samples, and with some actual data for those samples).
Good software design would obviously contain tests for exceptions and corner-cases, and the code did indeed try to accommodate a minimalistic data set such as yours (but did run in to additional problems). We may try to see if we can catch further exceptions (the code should not crash but end safely), but I fear that making sure that data such as yours will indeed give some output can be hard. The main developers of the coda have to decide if this is a priority.
Cheers Johan