[BUG] dnConstrainedTopology doesn't work with dnCDBDP/dnGLHBDSP

[willgearty]

Describe the bug
I'm currently trying to run an FBD BiSSE analysis using dnCDBDP/dnGLHBDSP. I'm following a combination of this BiSSE tutorial and this FBD tutorial. I've been able to run an FBD analysis with all of the clade constraints I want (based on an existing unresolved supertree), although I had to implement a few source code hacks (willgearty@3724075) that I found in various github issues related to having lots of constraints (thanks @bredelings and @wrightaprilm). Now I'm trying to convert my FBD analysis into a BiSSE analysis, but I always run into an issue when I try to enforce the constraints (using dnConstrainedTopology). I've tried this without success with both dnCDBDP (which results in the error cannot create std::vector larger than max_size()) and dnGLHBDSP (results in a segmentation fault). I initially was running this with an entire set of 109 constraints, but even just a single constraint causes the same errors. I thought it might be just constraints, but I've also tried just supplying an initialTree to dnConstrainedTopology and that also causes the same issues.

To Reproduce
The analysis can be reproduced using the files in this zip folder: croc_files.zip. Note that the analysis can be swapped from dnCDBDP to dnGLHBDSP by changing which lines are commented. I've also included a tree and commented code for supplying an initialTree instead of constraints.

Expected behavior
The tree distribution should be constrained and the analysis should continue.

Computer info
I'm using a slightly modified version of the current development branch (willgearty@3724075) on Ubuntu (via WSL2).

[brpetrucci]

Working on trying to solve it now, it seems like in the current revbayes dev branch I'm not getting a segfault with dnGLHBDSP, rather getting a weird error that literally just says vector on the dnConstrainedTopology line, both using clade28 and clade1.

[brpetrucci]

Tracked problem down to line 126 of Dist_GLHBDSP.cpp, and then to line 66 of GeneralizedLineageHeterogeneousBirthDeathSamplingProcess.cpp, i.e. to simulating the tree. Seems like the issue is specifically in the buildSerialSampledRandomBinaryTree member function (called in line 1073 of the GLHBDSP cpp file). The error seems to happen at exactly line 204 (TopologyNode* rightChild = active_nodes.at(right);), in iteration i = 119. That's the iteration where active_nodes becomes of size 0, so we ran out of nodes and therefore have no nodes to access for the right child (I assume this is where Will was getting a segfault as well, maybe some recent change made it so it's not a segfault anymore?).

[willgearty]

Thanks for looking into this @brpetrucci! Let me know if you need anything else from me.

[brpetrucci]

No worries!!

Solved the problem I mentioned above--it seems like the SimulateCoalescentAges method was simulating ages by drawing a uniform between the taxon's minimum age and the root age, which probably led most ages to be too high. Changed it to simulate ages between the min and max age of the taxon, and now am getting a segfault on line 397 of AbstractRootedTreeDistribution.cpp (a call to taxa.size()). Progress!(?)

[brpetrucci]

Gonna keep using this as a diary I guess, maybe someone has an epiphany I don't eventually.

Tracked down the call to taxa.size() to line 710 of TopologyConstrainedTreeDistribution.cpp. It's part of the simulateRootedTree member function, where we are just trying to get the number of taxa to start simulating. It seems like in like 709, when we take the member variable base_distribution (which is a Typed and typecast it to AbstractRootedTreeDistribution, the casting and/or the = operator are not working as intended, as the resulting distribution (tree_base_distribution) doesn't have a taxa member (or any member). So now just gotta figure out if the dynamic casting or the = operator is the problem. Since we haven't gotten this problem with other distributions, my guess is that the base_distribution variable coming from dnGLHBDSP has something wrong, but during debugging I can't see anything specifically wrong there.

I tried running just an FBD skyline model with constraints and it seems like the dynamic casting/= operator work fine, so I guess the problem when we try to make a GeneralizedLineageHeterogeneousBirthDeathSamplingProcess variable into a TypedDistribution or an AbstractRootedTreeDistribution. I think the problem is that we're trying to make a general distribution (TypedDistribution) into a more specific one (AbstractRootedTreeDistribution). If GLHBDSP was, like other BDSP classes in RevBayes, an AbstractBirthDeathProcess, I suspect that could solve the issue. Unfortunately, because GLHBDSP needs to use the TreeCharacterData class, it cannot be an AbstractBirthDeathProcess. I'm not sure what the solution is, outside of completely restructuring GLHBDSP to be able to inherit one of classes that inherit from AbstractRootedTreeDistribution.

[brpetrucci]

Set up a branch (dev_const_top_bug) to work on this. The only change right now is the fix to how GLHBDSP is calculating taxa ages during simulations (@mikeryanmay tagging you here so you can confirm this wasn't a mistake), which solved one problem. I'm currently not sure how to solve the segfault (i.e. how to get the GLHBDSP distribution to inherit something that can be cast as an AbstractRootedTreeDistribution) without completely restructuring the class inheritance scheme, but I'll keep thinking.

[bredelings]

@brpetrucci Thanks for looking into this.

I made a debug build using meson setup ../git gcc-13-asdf-debug --prefix=/home/bredelings/Devel/revbayes/local/gcc-13-asdf-debug --buildtype=debug. Basically, this turns off optimization (-O0) and turns on debugging (-g).

Then I used gdb to find the exception that was thrown. I usually use gdbtui.

 (gdb) catch throw

Finding the relevant exception is challenging because rb throws a lot of exceptions when determining what the argument types of functions are. These exceptions are not errors, and we need to skip them. Maybe there is a way to do that, but I don't know how.
So I did a binary search to try and find the last exception.

 (gdb) c 100
 (gdb) up
 (gdb) c 200
 (gdb) up
 (gdb) c 400
 (gdb) up
...

I finally found the exception we were looking for, and it was here:

│       46   */                                                                                                                            │
│       47  TopologyConstrainedTreeDistribution::TopologyConstrainedTreeDistribution(TypedDistribution<Tree>* base_dist, const std::vector<│
│       48  //    active_backbone_clades( base_dist->getValue().getNumberOfInteriorNodes(), RbBitSet() ),                                  │
│  >    49      active_clades( base_dist->getValue().getNumberOfInteriorNodes(), RbBitSet() ),                                             │
│       50      backbone_topology(NULL),                                                                                                   │

It turns out that base_dist->getValue().getNumberOfInteriorNodes() is returning -1, but casted to the unsigned type size_t, which yields 18446744073709551615.

The code for Tree::getNumberOfInteriorNodes() does:

/**
 * Calculate the number of interior nodes in the BranchLengthTree by deducing the number of
 * tips from number of nodes, and then subtract 1 more if the BranchLengthTree is rooted.
 */
size_t Tree::getNumberOfInteriorNodes( void ) const
{

    size_t preliminaryNumIntNodes = getNumberOfNodes() - getNumberOfTips();

    if ( isRooted() )
    {
	return preliminaryNumIntNodes - 1;
    }
    else
    {
        return preliminaryNumIntNodes;
    }

}

It turns out that getNumberOfNodes() == getNumberOfTips() == 1. Since the tree is rooted, this ends up returning 0-1 == -1, but casted to size_t.

If we change this to the following, then the crash doesn't happen.

	if (preliminaryNumIntNodes > 1)
	    return preliminaryNumIntNodes - 1;
	else
	    return 0;

Then we get a different error message:

> source("crocodylia_bisse.Rev")
   Processing file "crocodylia_bisse.Rev"
   Reading in the Data
   Successfully read one character matrix from file 'data/crocodylia_habitat2.nex'
   Setting up BiSSE
   Setting up FBD
   Error:	Could not find taxon with name 'Gavialis_bengawanicus'.
   Error:	Problem processing line 105 in file "crocodylia_bisse.Rev"
> 

[willgearty]

Thanks @bredelings for looking into this. I implemented the changes that you and @brpetrucci made on my forked branch, rebuilt rb, and I'm now running into a similar error (this is when using dnCDCladoBDP and the full set of constraints):

> source("scripts/crocodylia_bisse.Rev")
   Processing file "scripts/crocodylia_bisse.Rev"
   Reading in the Data
   Successfully read one character matrix from file 'data/crocodylia_habitat2.nex'
   Setting up BiSSE
   Setting up FBD
   Processing file "scripts/constraints_crocodylia.Rev"
   Processing of file "scripts/constraints_crocodylia.Rev" completed
   Error:       Could not find taxon with name 'Argochampsa_krebsi'.
   Error:       Problem processing line 105 in file "scripts/crocodylia_bisse.Rev"
>

When I use only the clade28 constraint (line 102), I get the exact same error as @bredelings.

I've double checked and both Argochampsa_krebsi and Gavialis_bengawanicus are in the taxon list, habitat character matrix, and the constraints file with the exact same names.

Further, when I switch to dnGLHBDSP, it still results in a segmentation fault:

> source("scripts/crocodylia_bisse.Rev")
   Processing file "scripts/crocodylia_bisse.Rev"
   Reading in the Data
   Successfully read one character matrix from file 'data/crocodylia_habitat2.nex'
   Setting up BiSSE
   Setting up FBD
   Processing file "scripts/constraints_crocodylia.Rev"
   Processing of file "scripts/constraints_crocodylia.Rev" completed
   Warning: simulating tree under uniform model.
   Warning: simulating tree under uniform model.
Segmentation fault

[willgearty]

I've played around with various random clade constraints and it just seems to complain about whichever taxon I give it (even if the clade constraint is a single taxon that is confirmed to be in the taxon list).

[willgearty]

Let me know if you need any other information from me @bredelings @brpetrucci.

[bredelings]

Hi @willgearty. The segmentation fault is fixed on the development branch. Now I'm seeing some problems with the script you provided. For example, if the clamping is commented out, then the sampled tree doesn't have any taxa associated with it.

Can you (or @brpetrucci) compile revbayes from the latest version of the development branch and provide a script that demonstrates the problem? Currently I'm having to guess what to uncomment.

[willgearty]

Hi @bredelings, thanks for getting back to me. I just pulled the revbayes/development branch, then added the changes suggested in this issue, and I'm still hitting a segmentation fault using dnGLHBDSP.

You should be able to replicate the issue with this set of files: crocodylia_bisse.zip. Note that the tensorPhylo plugin location will need to be changed on line 2.

If you comment line 105 and uncomment line 106, you should get the same error from our comments above using dnCDCladoBDP, which seems to reference taxa that are definitely included in the data files..

[bredelings]

Hmm... it looks like you forgot to include scripts/constraints_crocodylia.Rev

[willgearty]

I just realized I didn't put the files within the correct folders, I've updated the zip folder with the correct folder system here:
will_crocs (2).zip.

[bredelings]

Thanks!

It looks like the problem occurs here:

revbayes/src/core/distributions/phylogenetics/tree/TopologyConstrainedTreeDistribution.cpp

Line 709 in 3ad7874

AbstractRootedTreeDistribution* tree_base_distribution = dynamic_cast<AbstractRootedTreeDistribution*>( base_distribution );

The problem is that base_distribution has type (RevBayesCore::GeneralizedLineageHeterogeneousBirthDeathSamplingProcess *). But that type is not derived from AbstractRootedTreeDistribution* so it is not convertible to it. Non-convertibility is indicated by returning a NULL pointer. However the code doesn't check for a NULL pointer.

Two things should probably be fixed:

• first, GeneralizedLineageHeterogeneousBirthDeathSamplingProcess should probably be derived from AbstractRootedTreeDistribution instead of from TypedDistribution<Tree>.
• second, the code mentioned above should check if the returned pointer is NULL, and throw an exception with the message that says something like "TopologyConstrainedTreeDistribution::simulatedRootedTree( ): cannot simulate rooted tree from distribution "<<name<<" because it is not derived from AbstractRootedTreeDistribution".

Probably you should contact @mikeryanmay about fixing the first one of these.

[willgearty]

Thanks @bredelings, hopefully @mikeryanmay will be able to help with that.

And do any of you have any insight into the error with dnGLHBDSP where it says that taxa that clearly exist don't exist? #375 (comment)

[bredelings]

I think they actually don't exist at first. No taxa are given to dnGLHBDSP( ) here:

fbd_dist = dnGLHBDSP(rootAge = root_age,
                     lambda = lambda,
                     mu = mu,
                     eta = Q,
                     pi = pi,
                     phi = phi,
                     taxa = taxa,
                     rho = rho,
                     nStates = num_states)

The taxa only get added here:

# fbd_tree.clamp(tree)

So, I think that is actually correct for a version of the script where the clamp is commented out.

[willgearty]

I'm confused, there is indeed taxa = taxa in that dnGLHBDSP call? And I've checked the file that the taxa object is coming from, and those "error" taxa definitely exist.

[bredelings]

Oh, right, its that no taxa are given to dnCDCladoBDP( ). dnGLHBDSP crashes so we wouldn't get that far with it.

fbd_dist = dnCDCladoBDP(rootAge = root_age,
                        speciationRates = lambda,
                        extinctionRates = mu,
                        Q = Q,
                        pi = pi,
                        phi = phi,
                        rho = rho)

[willgearty]

Oh, yes, you're right, thanks for pointing that out!

[willgearty]

Wait, dnCDCladoBDP doesn't have a taxa argument. And I can't clamp the taxa object to the fbd_dist object. How am I supposed to provide taxa to the function? 😕

[willgearty]

@bredelings or @brpetrucci any ideas?

[bredelings]

I'm not sure. You might try asking on the revbayes-users mailing list. You could also open a new issue about the dnCDCladoBDP not having any taxa parameter.

My guess is that dnCDCladoBDP has previously only been used in situations where the tree is clamped, and in those cases the taxa can come from the observed phylogeny. For example:

classe ~ dnCDCladoBDP( rootAge         = observed_phylogeny.rootAge(),
                       cladoEventMap   = clado_matrix,
                       extinctionRates = extinction_rates,
                       Q               = ana_rate_matrix,
                       delta           = 1.0,
                       pi              = root_frequencies,
                       rho             = rho,
                       condition       = "time",
                       taxa            = taxa )

# clamp the model with the observed data
classe.clamp( observed_phylogeny )
classe.clampCharData( data_biogeo )

But I have never used it myself.