Skip to content

Commit

Permalink
minor updates on species tree proposals
Browse files Browse the repository at this point in the history
  • Loading branch information
@hoehna
hoehna committed Jul 13, 2015
1 parent df45ff2 commit c016d2d
Show file tree
Hide file tree
Showing 6 changed files with 70 additions and 40 deletions.
2 changes: 1 addition & 1 deletion examples/JukesCantor.Rev
View file
Expand Up @@ -22,7 +22,7 @@ data <- readDiscreteCharacterData("data/primates_cytb.nex")
# Get some useful variables from the data. We need these later on.
n_species <- data.ntaxa()
names <- data.names()

# set my move index
mi = 0

Expand Down
1 change: 1 addition & 0 deletions projects/cmake/regenerate.sh
View file
Expand Up @@ -141,6 +141,7 @@ set(PROJECT_SOURCE_DIR ${CMAKE_SOURCE_DIR}/../../src)
SET(BOOST_ROOT ../../boost_1_55_0)
SET(Boost_USE_STATIC_RUNTIME true)
SET(Boost_USE_STATIC_LIBS ON)
#find_package(Boost 1.55.0 COMPONENTS filesystem regex signals context system thread date_time program_options iostreams serialization math_c99 math_c99f math_tr1f math_tr1l REQUIRED)
find_package(Boost
1.55.0
Expand Down
18 changes: 8 additions & 10 deletions src/core/distributions/phylogenetics/tree/MultispeciesCoalescent.cpp
View file
Expand Up @@ -19,10 +19,15 @@ MultispeciesCoalescent::MultispeciesCoalescent(const TypedDagNode<TimeTree> *sp,
taxa(t),
speciesTree( sp ),
Nes( NULL ),
Ne( NULL ),
Ne( new ConstantNode<double>("Ne", new double(1.0) ) ),
numTaxa( taxa.size() ),
logTreeTopologyProb (0.0)
{
// add the parameters to our set (in the base class)
// in that way other class can easily access the set of our parameters
// this will also ensure that the parameters are not getting deleted before we do
addParameter( speciesTree );
addParameter( Ne );

std::set<std::string> speciesNames;
for (std::vector<Taxon>::const_iterator it=taxa.begin(); it!=taxa.end(); ++it)
Expand All @@ -39,16 +44,8 @@ MultispeciesCoalescent::MultispeciesCoalescent(const TypedDagNode<TimeTree> *sp,

logTreeTopologyProb = (numTaxa - 1) * RbConstants::LN2 - 2.0 * lnFact - std::log( numTaxa ) ;

//Default value for Ne
Ne = new ConstantNode<double>("Ne", new double(1.0) );

redrawValue();

// add the parameters to our set (in the base class)
// in that way other class can easily access the set of our parameters
// this will also ensure that the parameters are not getting deleted before we do
addParameter( speciesTree );
addParameter( Ne );


}
Expand Down Expand Up @@ -602,7 +599,8 @@ double MultispeciesCoalescent::getNe(size_t index) const



void MultispeciesCoalescent::redrawValue( void ) {
void MultispeciesCoalescent::redrawValue( void )
{

simulateTree();

Expand Down
39 changes: 33 additions & 6 deletions src/core/moves/MetropolisHastingsMove.cpp
View file
Expand Up @@ -212,6 +212,21 @@ void MetropolisHastingsMove::performMove( double lHeat, double pHeat )
proposal->prepareProposal();
double lnHastingsRatio = proposal->doProposal();


affectedNodes.clear();
for (std::set<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it)
{
(*it)->getAffectedNodes( affectedNodes );
}
for (std::set<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it)
{
DagNode *the_node = *it;
if ( affectedNodes.find(the_node) != affectedNodes.end() )
{
affectedNodes.erase(the_node);
}
}

// first we touch all the nodes
// that will set the flags for recomputation
for (std::set<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it)
Expand All @@ -222,34 +237,41 @@ void MetropolisHastingsMove::performMove( double lHeat, double pHeat )
double lnPriorRatio = 0.0;
double lnLikelihoodRatio = 0.0;

std::cerr << "\nCore Node:\n";

// compute the probability of the current value for each node
for (std::set<DagNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it)
{
DagNode *the_node = *it;
std::cerr << the_node->getName() << " <" << the_node << ">" << std::endl;
if ( RbMath::isAComputableNumber(lnPriorRatio) && RbMath::isAComputableNumber(lnLikelihoodRatio) && RbMath::isAComputableNumber(lnHastingsRatio) )
{
if ( (*it)->isClamped() )
if ( the_node->isClamped() )
{
lnLikelihoodRatio += (*it)->getLnProbabilityRatio();
lnLikelihoodRatio += the_node->getLnProbabilityRatio();
}
else
{
lnPriorRatio += (*it)->getLnProbabilityRatio();
lnPriorRatio += the_node->getLnProbabilityRatio();
}
}
}

std::cerr << "\nAffected Node:\n";
// then we recompute the probability for all the affected nodes
for (std::set<DagNode*>::iterator it = affectedNodes.begin(); it != affectedNodes.end(); ++it)
{
DagNode *the_node = *it;
std::cerr << the_node->getName() << " <" << the_node << ">" << std::endl;
if ( RbMath::isAComputableNumber(lnPriorRatio) && RbMath::isAComputableNumber(lnLikelihoodRatio) && RbMath::isAComputableNumber(lnHastingsRatio) )
{
if ( (*it)->isClamped() )
if ( the_node->isClamped() )
{
lnLikelihoodRatio += (*it)->getLnProbabilityRatio();
lnLikelihoodRatio += the_node->getLnProbabilityRatio();
}
else
{
lnPriorRatio += (*it)->getLnProbabilityRatio();
lnPriorRatio += the_node->getLnProbabilityRatio();
}
}
}
Expand All @@ -260,6 +282,7 @@ void MetropolisHastingsMove::performMove( double lHeat, double pHeat )

if ( RbMath::isAComputableNumber(lnPosteriorRatio) == false )
{
std::cerr << "Reject\n";

proposal->undoProposal();

Expand All @@ -277,6 +300,7 @@ void MetropolisHastingsMove::performMove( double lHeat, double pHeat )

if (lnAcceptanceRatio >= 0.0)
{
std::cerr << "Accept\n";
numAccepted++;

// call accept for each node
Expand All @@ -288,6 +312,7 @@ void MetropolisHastingsMove::performMove( double lHeat, double pHeat )
}
else if (lnAcceptanceRatio < -300.0)
{
std::cerr << "Reject\n";
proposal->undoProposal();

// call restore for each node
Expand All @@ -303,6 +328,7 @@ void MetropolisHastingsMove::performMove( double lHeat, double pHeat )
double u = GLOBAL_RNG->uniform01();
if (u < r)
{
std::cerr << "Accept\n";
numAccepted++;

// call accept for each node
Expand All @@ -315,6 +341,7 @@ void MetropolisHastingsMove::performMove( double lHeat, double pHeat )
}
else
{
std::cerr << "Reject\n";
proposal->undoProposal();

// call restore for each node
Expand Down
24 changes: 13 additions & 11 deletions src/core/moves/proposal/tree/SpeciesSubtreeScaleBetaProposal.cpp
View file
Expand Up @@ -104,28 +104,30 @@ double SpeciesSubtreeScaleBetaProposal::doProposal( void )
TreeUtilities::getOldestTip(&tau, node, min_age);

// draw new ages
double currentValue = my_age / (parent_age - min_age);
double current_value = my_age / (parent_age - min_age);
double a = alpha + 1.0;
double b = (a-1.0) / currentValue - a + 2.0;
double b = (a-1.0) / current_value - a + 2.0;
double new_value = RbStatistics::Beta::rv(a, b, *rng);

new_value = current_value;

double my_new_age = new_value * (parent_age - min_age);

double scalingFactor = my_new_age / my_age;
double scaling_factor = my_new_age / my_age;

size_t nNodes = node->getNumberOfNodesInSubtree( false );
size_t num_nodes = node->getNumberOfNodesInSubtree( false );

for ( size_t i=0; i<geneTrees.size(); ++i )
{
// get the i-th gene tree
TimeTree& geneTree = geneTrees[i]->getValue();
TimeTree& gene_tree = geneTrees[i]->getValue();

std::vector<TopologyNode*> nodes = getOldestNodesInPopulation(geneTree, *node );
std::vector<TopologyNode*> nodes = getOldestNodesInPopulation(gene_tree, *node );

for (size_t j=0; j<nodes.size(); ++j)
{
// add the number of nodes that we are going to scale in the subtree
nNodes += nodes[j]->getNumberOfNodesInSubtree( false );
num_nodes += nodes[j]->getNumberOfNodesInSubtree( false );

if ( nodes[j]->isTip() == true )
{
Expand All @@ -138,7 +140,7 @@ double SpeciesSubtreeScaleBetaProposal::doProposal( void )
}

// rescale the subtree of this gene tree
TreeUtilities::rescaleSubtree(&geneTree, nodes[j], scalingFactor );
TreeUtilities::rescaleSubtree(&gene_tree, nodes[j], scaling_factor );

}

Expand All @@ -157,14 +159,14 @@ double SpeciesSubtreeScaleBetaProposal::doProposal( void )
// }

// rescale the subtree of the species tree
TreeUtilities::rescaleSubtree(&tau, node, scalingFactor );
TreeUtilities::rescaleSubtree(&tau, node, scaling_factor );

// compute the Hastings ratio
double forward = RbStatistics::Beta::lnPdf(a, b, new_value);
double new_a = alpha + 1.0;
double new_b = (a-1.0) / new_value - a + 2.0;
double backward = RbStatistics::Beta::lnPdf(new_a, new_b, currentValue);
double lnHastingsratio = (backward - forward) * (nNodes-1);
double backward = RbStatistics::Beta::lnPdf(new_a, new_b, current_value);
double lnHastingsratio = (backward - forward) * (num_nodes-1);

return lnHastingsratio;
}
Expand Down
26 changes: 14 additions & 12 deletions src/core/moves/proposal/tree/SpeciesSubtreeScaleProposal.cpp
View file
Expand Up @@ -104,21 +104,23 @@ double SpeciesSubtreeScaleProposal::doProposal( void )
// draw new ages and compute the hastings ratio at the same time
double my_new_age = min_age + (parent_age - min_age) * rng->uniform01();

double scalingFactor = my_new_age / my_age;
my_new_age = my_age;

size_t nNodes = node->getNumberOfNodesInSubtree( false );
double scaling_factor = my_new_age / my_age;

size_t num_nodes = node->getNumberOfNodesInSubtree( false );

for ( size_t i=0; i<geneTrees.size(); ++i )
{
// get the i-th gene tree
TimeTree& geneTree = geneTrees[i]->getValue();
TimeTree& gene_tree = geneTrees[i]->getValue();

std::vector<TopologyNode*> nodes = getOldestNodesInPopulation(geneTree, *node );
std::vector<TopologyNode*> nodes = getOldestNodesInPopulation(gene_tree, *node );

for (size_t j=0; j<nodes.size(); ++j)
{
// add the number of nodes that we are going to scale in the subtree
nNodes += nodes[j]->getNumberOfNodesInSubtree( false );
num_nodes += nodes[j]->getNumberOfNodesInSubtree( false );

if ( nodes[j]->isTip() == true )
{
Expand All @@ -131,7 +133,7 @@ double SpeciesSubtreeScaleProposal::doProposal( void )
}

// rescale the subtree of this gene tree
TreeUtilities::rescaleSubtree(&geneTree, nodes[j], scalingFactor );
TreeUtilities::rescaleSubtree(&gene_tree, nodes[j], scaling_factor );

}

Expand All @@ -151,10 +153,10 @@ double SpeciesSubtreeScaleProposal::doProposal( void )


// rescale the subtree of the species tree
TreeUtilities::rescaleSubtree(&tau, node, scalingFactor );
TreeUtilities::rescaleSubtree(&tau, node, scaling_factor );

// compute the Hastings ratio
double lnHastingsratio = (nNodes > 1 ? log( scalingFactor ) * (nNodes-1) : 0.0 );
double lnHastingsratio = (num_nodes > 1 ? log( scaling_factor ) * (num_nodes-1) : 0.0 );

return lnHastingsratio;
}
Expand All @@ -174,14 +176,14 @@ std::vector<TopologyNode*> SpeciesSubtreeScaleProposal::getOldestNodesInPopulati
}

// get all the taxa from the species tree that are descendants of node i
std::vector<TopologyNode*> speciesTaxa;
TreeUtilities::getTaxaInSubtree( &n, speciesTaxa );
std::vector<TopologyNode*> species_taxa;
TreeUtilities::getTaxaInSubtree( &n, species_taxa );

// get all the individuals
std::set<TopologyNode*> individualTaxa;
for (size_t i = 0; i < speciesTaxa.size(); ++i)
for (size_t i = 0; i < species_taxa.size(); ++i)
{
const std::string &name = speciesTaxa[i]->getName();
const std::string &name = species_taxa[i]->getName();
std::vector<TopologyNode*> ind = tau.getTipNodesWithSpeciesName( name );
for (size_t j = 0; j < ind.size(); ++j)
{
Expand Down

0 comments on commit c016d2d

Please sign in to comment.