From aca18599c14a3f87d0910ab3c611a739d1d44bec Mon Sep 17 00:00:00 2001 From: boussau Date: Fri, 4 Aug 2017 16:18:52 +0200 Subject: [PATCH 1/2] Better help for the species tree moves in the MSC framework. --- .../tree/Move_SpeciesNodeTimeSlideUniform.cpp | 231 ++++++++++++++--- .../tree/Move_SpeciesNodeTimeSlideUniform.h | 29 ++- .../moves/tree/Move_SpeciesSubtreeScale.cpp | 230 ++++++++++++++--- .../moves/tree/Move_SpeciesSubtreeScale.h | 31 ++- .../tree/Move_SpeciesSubtreeScaleBeta.cpp | 232 +++++++++++++++--- .../moves/tree/Move_SpeciesSubtreeScaleBeta.h | 37 +-- .../moves/tree/Move_SpeciesTreeScale.cpp | 229 ++++++++++++++--- .../moves/tree/Move_SpeciesTreeScale.h | 34 ++- 8 files changed, 885 insertions(+), 168 deletions(-) diff --git a/src/revlanguage/moves/tree/Move_SpeciesNodeTimeSlideUniform.cpp b/src/revlanguage/moves/tree/Move_SpeciesNodeTimeSlideUniform.cpp index 3cd1c8c7b..0f0b26a72 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesNodeTimeSlideUniform.cpp +++ b/src/revlanguage/moves/tree/Move_SpeciesNodeTimeSlideUniform.cpp @@ -20,12 +20,12 @@ using namespace RevLanguage; */ Move_SpeciesNodeTimeSlideUniform::Move_SpeciesNodeTimeSlideUniform() : Move() { - + // add method for call "addGeneTreeVariable" as a function ArgumentRules* addGeneTreeArgRules = new ArgumentRules(); addGeneTreeArgRules->push_back( new ArgumentRule( "geneTree" , TimeTree::getClassTypeSpec(), "A gene tree.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); methods.addFunction( new MemberProcedure( "addGeneTreeVariable", RlUtils::Void, addGeneTreeArgRules) ); - + } @@ -37,7 +37,7 @@ Move_SpeciesNodeTimeSlideUniform::Move_SpeciesNodeTimeSlideUniform() : Move() */ Move_SpeciesNodeTimeSlideUniform* Move_SpeciesNodeTimeSlideUniform::clone(void) const { - + return new Move_SpeciesNodeTimeSlideUniform(*this); } @@ -56,15 +56,15 @@ void Move_SpeciesNodeTimeSlideUniform::constructInternalObject( void ) { // we free the memory first delete value; - + // now allocate a new move double w = static_cast( weight->getRevObject() ).getValue(); RevBayesCore::TypedDagNode* tmp = static_cast( speciesTree->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *st = static_cast *>( tmp ); - + RevBayesCore::Proposal *p = new RevBayesCore::TreeNodeAgeUpdateProposal(st); value = new RevBayesCore::MetropolisHastingsMove(p,w); - + } @@ -75,9 +75,9 @@ void Move_SpeciesNodeTimeSlideUniform::constructInternalObject( void ) */ const std::string& Move_SpeciesNodeTimeSlideUniform::getClassType(void) { - + static std::string rev_type = "Move_SpeciesNodeTimeSlideUniform"; - + return rev_type; } @@ -89,9 +89,9 @@ const std::string& Move_SpeciesNodeTimeSlideUniform::getClassType(void) */ const TypeSpec& Move_SpeciesNodeTimeSlideUniform::getClassTypeSpec(void) { - + static TypeSpec rev_type_spec = TypeSpec( getClassType(), new TypeSpec( Move::getClassTypeSpec() ) ); - + return rev_type_spec; } @@ -105,7 +105,7 @@ std::string Move_SpeciesNodeTimeSlideUniform::getMoveName( void ) const { // create a constructor function name variable that is the same for all instance of this class std::string c_name = "SpeciesNodeTimeSlideUniform"; - + return c_name; } @@ -120,21 +120,21 @@ std::string Move_SpeciesNodeTimeSlideUniform::getMoveName( void ) const */ const MemberRules& Move_SpeciesNodeTimeSlideUniform::getParameterRules(void) const { - + static MemberRules memberRules; static bool rules_set = false; - + if ( !rules_set ) { - memberRules.push_back( new ArgumentRule( "speciesTree", TimeTree::getClassTypeSpec() , "The species tree on which this move operates.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); - + memberRules.push_back( new ArgumentRule( "speciesTree", TimeTree::getClassTypeSpec() , "The ultrametric species tree on which this move operates.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); + /* Inherit weight from Move, put it after variable */ const MemberRules& inheritedRules = Move::getParameterRules(); memberRules.insert( memberRules.end(), inheritedRules.begin(), inheritedRules.end() ); - + rules_set = true; } - + return memberRules; } @@ -146,43 +146,43 @@ const MemberRules& Move_SpeciesNodeTimeSlideUniform::getParameterRules(void) con */ const TypeSpec& Move_SpeciesNodeTimeSlideUniform::getTypeSpec( void ) const { - + static TypeSpec type_spec = getClassTypeSpec(); - + return type_spec; } RevPtr Move_SpeciesNodeTimeSlideUniform::executeMethod(const std::string& name, const std::vector& args, bool &found) { - + if ( name == "addGeneTreeVariable" ) { found = true; - + RevBayesCore::TypedDagNode* tmp = static_cast( args[0].getVariable()->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *gt = static_cast *>( tmp ); - + RevBayesCore::MetropolisHastingsMove *m = static_cast(this->value); RevBayesCore::TreeNodeAgeUpdateProposal &p = static_cast( m->getProposal() ); p.addGeneTree( gt ); - + return NULL; } else if ( name == "removeGeneTreeVariable" ) { found = true; - + RevBayesCore::TypedDagNode* tmp = static_cast( args[0].getVariable()->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *gt = static_cast *>( tmp ); - + RevBayesCore::MetropolisHastingsMove *m = static_cast(this->value); RevBayesCore::TreeNodeAgeUpdateProposal &p = static_cast( m->getProposal() ); p.removeGeneTree( gt ); - + return NULL; } - + return Move::executeMethod( name, args, found ); } @@ -193,7 +193,7 @@ RevPtr Move_SpeciesNodeTimeSlideUniform::executeMethod(const std::s */ void Move_SpeciesNodeTimeSlideUniform::printValue(std::ostream &o) const { - + o << "SpeciesNodeTimeSlideUniform("; if (speciesTree != NULL) { @@ -204,7 +204,7 @@ void Move_SpeciesNodeTimeSlideUniform::printValue(std::ostream &o) const o << "?"; } o << ")"; - + } @@ -220,7 +220,7 @@ void Move_SpeciesNodeTimeSlideUniform::printValue(std::ostream &o) const */ void Move_SpeciesNodeTimeSlideUniform::setConstParameter(const std::string& name, const RevPtr &var) { - + if ( name == "speciesTree" ) { speciesTree = var; @@ -229,9 +229,178 @@ void Move_SpeciesNodeTimeSlideUniform::setConstParameter(const std::string& name { Move::setConstParameter(name, var); } - + } + + /** + * Get the author(s) of this function so they can receive credit (and blame) for it. + */ + std::vector Move_SpeciesNodeTimeSlideUniform::getHelpAuthor(void) const + { + // create a vector of authors for this function + std::vector authors; + authors.push_back( "Sebastian Hoehna, Bastien Boussau" ); + + return authors; + } + + + /** + * Get the (brief) description for this function + */ + std::vector Move_SpeciesNodeTimeSlideUniform::getHelpDescription(void) const + { + // create a variable for the description of the function + std::vector descriptions; + descriptions.push_back( "Makes a node time slide move both in the species tree and in the gene trees that contain nodes of the relevant populations. Tree topologies are not altered." ); + + return descriptions; + } + + + /** + * Get the more detailed description of the function + */ + std::vector Move_SpeciesNodeTimeSlideUniform::getHelpDetails(void) const + { + // create a variable for the description of the function + std::vector details; + + std::string details_1 = ""; + details_1 += "The species tree must be ultrametric."; + + details.push_back( details_1 ); + + std::string details_2 = ""; + details_2 += "All the gene trees that evolved along the species tree according to some form of multispecies coalescent must be added to the move using the addGeneTreeVariable method. "; + + details.push_back( details_2 ); + + std::string details_3 = ""; + details_3 += "This move jointly performs node time slides (branch length alterations, keeping the topologies fixed) on the species tree and on gene trees, all of which must be ultrametric. "; + + details.push_back( details_3 ); + + return details; + } + + + /** + * Get an executable and instructive example. + * These examples should help the users to show how this function works but + * are also used to test if this function still works. + */ + std::string Move_SpeciesNodeTimeSlideUniform::getHelpExample(void) const + { + // create an example as a single string variable. + std::string example = ""; + example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; + example += "dataFolder = \"simulatedTrees\" \n"; + example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; + example += "n_species <- 10\n"; + example += "n_genes <- 2\n"; + example += "n_alleles <- 3\n"; + example += "# we simulate an ultrametric species tree:\n"; + example += "# Species names:\n"; + example += "for (i in 1:n_species) {\n"; + example += " species[i] <- taxon(taxonName=\"Species_\"+i, speciesName=\"Species_\"+i)\n"; + example += "}\n"; + example += "spTree ~ dnBirthDeath(lambda=0.3, mu=0.2, rootAge=10, rho=1, samplingStrategy=\"uniform\", condition=\"nTaxa\", taxa=species)\n"; + example += "print(spTree)\n"; + example += "# let's pick a constant effective population size of 50:\n"; + example += "popSize <- 50\n"; + example += "# let's simulate gene trees now:\n"; + example += "# taxa names:\n"; + example += "for (g in 1:n_genes) {\n"; + example += " for (i in 1:n_species) {\n"; + example += " for (j in 1:n_alleles) {\n"; + example += " taxons[g][(i-1)*n_alleles+j] <- taxon(taxonName=\"Species_\"+i+\"_\"+j, speciesName=\"Species_\"+i)\n"; + example += " }\n"; + example += " }\n"; + example += " geneTrees[g] ~ dnMultiSpeciesCoalescent(speciesTree=spTree, Ne=popSize, taxa=taxons[g])\n"; + example += " print(geneTrees[g])\n"; + example += "}\n"; + example += "# We can save the species tree and the gene trees: \n"; + example += "write(speciesTree, filename=dataFolder+\"speciesTree\")\n"; + example += "# Saving the gene trees\n"; + example += "for (i in 1:(n_genes)) {\n"; + example += " write(geneTrees[i], filename=dataFolder+\"geneTree_\"+i+\".tree\")\n"; + example += "}\n"; + example += "# set my move index\n"; + example += "mi = 0\n"; + example += "move_species_node_time_slide = mvSpeciesNodeTimeSlideUniform( speciesTree=psi, weight=5 )\n"; + example += "for (i in 1:n_genes) {\n"; + example += " move_species_node_time_slide.addGeneTreeVariable( geneTree[i] )\n"; + example += "}\n"; + example += "moves[++mi] = move_species_node_time_slide\n"; + example += "# We get a handle on our model.\n"; + example += "# We can use any node of our model as a handle, here we choose to use the topology.\n"; + example += "mymodel = model(spTree)\n"; + example += "# Monitors to check the progression of the program\n"; + example += "monitors[1] = mnScreen(printgen=10, sptree)\n"; + example += "# Here we use a plain MCMC. You could also set nruns=2 for a replicated analysis\n"; + example += "# or use mcmcmc with heated chains.\n"; + example += "mymcmc = mcmc(mymodel, monitors, moves, nruns=4)\n"; + example += "mymcmc.run(generations=1000)\n"; + example += "mymcmc.operatorSummary()\n"; + + return example; + } + + + /** + * Get some references/citations for this function + * + */ + std::vector Move_SpeciesNodeTimeSlideUniform::getHelpReferences(void) const + { + // create an entry for each reference + std::vector references; + RevBayesCore::RbHelpReference ref = RevBayesCore::RbHelpReference(); + ref.setCitation("Guided Tree Topology Proposals for Bayesian Phylogenetic Inference. Sebastian Höhna, Alexei J. Drummond. Syst Biol (2012) 61 (1): 1-11."); + ref.setDoi("https://doi.org/10.1093/sysbio/syr074"); + ref.setUrl("https://academic.oup.com/sysbio/article-lookup/doi/10.1093/sysbio/syr074"); + + references.push_back(ref); + + RevBayesCore::RbHelpReference ref2 = RevBayesCore::RbHelpReference(); + ref2.setCitation("Algorithmic improvements to species delimitation and phylogeny estimation under the multispecies coalescent. Graham Jones. Journal of Mathematical Biology, 2016."); + ref2.setDoi("https://doi.org/10.1007/s00285-016-1034-0"); + ref2.setUrl("http://link.springer.com/article/10.1007/s00285-016-1034-0"); + + references.push_back(ref2); + + return references; + } + + + /** + * Get the names of similar and suggested other functions + */ + std::vector Move_SpeciesNodeTimeSlideUniform::getHelpSeeAlso(void) const + { + // create an entry for each suggested function + std::vector see_also; + see_also.push_back( "mvSpeciesSubtreeScale" ); + see_also.push_back( "mvSpeciesSubtreeScaleBeta" ); + see_also.push_back( "mvSpeciesNarrow" ); + see_also.push_back( "mvSpeciesTreeScale" ); + + return see_also; + } + + + /** + * Get the title of this help entry + */ + std::string Move_SpeciesNodeTimeSlideUniform::getHelpTitle(void) const + { + // create a title variable + std::string title = "Node time slide joint move on species tree and gene trees for multispecies coalescent models. "; + + return title; + } diff --git a/src/revlanguage/moves/tree/Move_SpeciesNodeTimeSlideUniform.h b/src/revlanguage/moves/tree/Move_SpeciesNodeTimeSlideUniform.h index 186190dfd..6d23f00c9 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesNodeTimeSlideUniform.h +++ b/src/revlanguage/moves/tree/Move_SpeciesNodeTimeSlideUniform.h @@ -8,8 +8,8 @@ #include namespace RevLanguage { - - + + /** * The RevLanguage wrapper of the species-node-time-slide (uniform) move. * @@ -26,11 +26,11 @@ namespace RevLanguage { * */ class Move_SpeciesNodeTimeSlideUniform : public Move { - + public: - + Move_SpeciesNodeTimeSlideUniform(void); //!< Default constructor - + // Basic utility functions virtual Move_SpeciesNodeTimeSlideUniform* clone(void) const; //!< Clone object void constructInternalObject(void); //!< We construct the a new internal SlidingMove. @@ -40,19 +40,28 @@ namespace RevLanguage { const MemberRules& getParameterRules(void) const; //!< Get member rules (const) virtual const TypeSpec& getTypeSpec(void) const; //!< Get language type of the object virtual void printValue(std::ostream& o) const; //!< Print value (for user) - + // Member method functions virtual RevPtr executeMethod(const std::string& name, const std::vector& args, bool &f); //!< Map member methods to internal functions protected: - + + + std::vector getHelpAuthor(void) const; //!< Get the author(s) of this function + std::vector getHelpDescription(void) const; //!< Get the description for this function + std::vector getHelpDetails(void) const; //!< Get the more detailed description of the function + std::string getHelpExample(void) const; //!< Get an executable and instructive example + std::vector getHelpReferences(void) const; //!< Get some references/citations for this function + std::vector getHelpSeeAlso(void) const; //!< Get suggested other functions + std::string getHelpTitle(void) const; //!< Get the title of this help entry + void setConstParameter(const std::string& name, const RevPtr &var); //!< Set member variable - + RevPtr speciesTree; RevPtr geneTrees; - + }; - + } #endif diff --git a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.cpp b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.cpp index 140dc35b6..6c5b85864 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.cpp +++ b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.cpp @@ -20,7 +20,7 @@ using namespace RevLanguage; */ Move_SpeciesSubtreeScale::Move_SpeciesSubtreeScale() : Move() { - + // add method for call "addGeneTreeVariable" as a function ArgumentRules* addGeneTreeArgRules = new ArgumentRules(); addGeneTreeArgRules->push_back( new ArgumentRule( "geneTree" , TimeTree::getClassTypeSpec(), "A gene tree to scale.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); @@ -36,7 +36,7 @@ Move_SpeciesSubtreeScale::Move_SpeciesSubtreeScale() : Move() */ Move_SpeciesSubtreeScale* Move_SpeciesSubtreeScale::clone(void) const { - + return new Move_SpeciesSubtreeScale(*this); } @@ -55,49 +55,49 @@ void Move_SpeciesSubtreeScale::constructInternalObject( void ) { // we free the memory first delete value; - + // now allocate a new move double w = static_cast( weight->getRevObject() ).getValue(); RevBayesCore::TypedDagNode* tmp = static_cast( speciesTree->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *st = static_cast *>( tmp ); - + RevBayesCore::Proposal *p = new RevBayesCore::SpeciesSubtreeScaleProposal(st); value = new RevBayesCore::MetropolisHastingsMove(p,w); - + } RevPtr Move_SpeciesSubtreeScale::executeMethod(const std::string& name, const std::vector& args, bool &found) { - + if ( name == "addGeneTreeVariable" ) { found = true; - + RevBayesCore::TypedDagNode* tmp = static_cast( args[0].getVariable()->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *gt = static_cast *>( tmp ); - + RevBayesCore::MetropolisHastingsMove *m = static_cast(this->value); RevBayesCore::SpeciesSubtreeScaleProposal &p = static_cast( m->getProposal() ); p.addGeneTree( gt ); - + return NULL; } else if ( name == "removeGeneTreeVariable" ) { found = true; - + RevBayesCore::TypedDagNode* tmp = static_cast( args[0].getVariable()->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *gt = static_cast *>( tmp ); - + RevBayesCore::MetropolisHastingsMove *m = static_cast(this->value); RevBayesCore::SpeciesSubtreeScaleProposal &p = static_cast( m->getProposal() ); p.removeGeneTree( gt ); - + return NULL; } - + return Move::executeMethod( name, args, found ); } @@ -109,9 +109,9 @@ RevPtr Move_SpeciesSubtreeScale::executeMethod(const std::string& n */ const std::string& Move_SpeciesSubtreeScale::getClassType(void) { - + static std::string rev_type = "Move_SpeciesSubtreeScale"; - + return rev_type; } @@ -123,9 +123,9 @@ const std::string& Move_SpeciesSubtreeScale::getClassType(void) */ const TypeSpec& Move_SpeciesSubtreeScale::getClassTypeSpec(void) { - + static TypeSpec rev_type_spec = TypeSpec( getClassType(), new TypeSpec( Move::getClassTypeSpec() ) ); - + return rev_type_spec; } @@ -139,7 +139,7 @@ std::string Move_SpeciesSubtreeScale::getMoveName( void ) const { // create a constructor function name variable that is the same for all instance of this class std::string c_name = "SpeciesSubtreeScale"; - + return c_name; } @@ -154,21 +154,21 @@ std::string Move_SpeciesSubtreeScale::getMoveName( void ) const */ const MemberRules& Move_SpeciesSubtreeScale::getParameterRules(void) const { - + static MemberRules memberRules; static bool rules_set = false; - + if ( !rules_set ) { memberRules.push_back( new ArgumentRule( "speciesTree", TimeTree::getClassTypeSpec(), "The species variable on which this move operates.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); - + /* Inherit weight from Move, put it after variable */ const MemberRules& inheritedRules = Move::getParameterRules(); memberRules.insert( memberRules.end(), inheritedRules.begin(), inheritedRules.end() ); - + rules_set = true; } - + return memberRules; } @@ -180,9 +180,9 @@ const MemberRules& Move_SpeciesSubtreeScale::getParameterRules(void) const */ const TypeSpec& Move_SpeciesSubtreeScale::getTypeSpec( void ) const { - + static TypeSpec type_spec = getClassTypeSpec(); - + return type_spec; } @@ -192,7 +192,7 @@ const TypeSpec& Move_SpeciesSubtreeScale::getTypeSpec( void ) const */ void Move_SpeciesSubtreeScale::printValue(std::ostream &o) const { - + o << "SpeciesSubtreeScale("; if (speciesTree != NULL) { @@ -203,7 +203,7 @@ void Move_SpeciesSubtreeScale::printValue(std::ostream &o) const o << "?"; } o << ")"; - + } @@ -219,7 +219,7 @@ void Move_SpeciesSubtreeScale::printValue(std::ostream &o) const */ void Move_SpeciesSubtreeScale::setConstParameter(const std::string& name, const RevPtr &var) { - + if ( name == "speciesTree" ) { speciesTree = var; @@ -228,9 +228,181 @@ void Move_SpeciesSubtreeScale::setConstParameter(const std::string& name, const { Move::setConstParameter(name, var); } - + } + /** + * Get the author(s) of this function so they can receive credit (and blame) for it. + */ + std::vector Move_SpeciesSubtreeScale::getHelpAuthor(void) const + { + // create a vector of authors for this function + std::vector authors; + authors.push_back( "Sebastian Hoehna, Bastien Boussau" ); + + return authors; + } + + + /** + * Get the (brief) description for this function + */ + std::vector Move_SpeciesSubtreeScale::getHelpDescription(void) const + { + // create a variable for the description of the function + std::vector descriptions; + descriptions.push_back( "Makes a subtree scale move both in the species tree and in the gene trees that contain nodes of the relevant populations. Tree topologies are not altered." ); + + return descriptions; + } + + + /** + * Get the more detailed description of the function + */ + std::vector Move_SpeciesSubtreeScale::getHelpDetails(void) const + { + // create a variable for the description of the function + std::vector details; + + std::string details_1 = ""; + details_1 += "The species tree must be ultrametric."; + + details.push_back( details_1 ); + + std::string details_2 = ""; + details_2 += "All the gene trees that evolved along the species tree according to some form of multispecies coalescent must be added to the move using the addGeneTreeVariable method. "; + + details.push_back( details_2 ); + + std::string details_3 = ""; + details_3 += "This move jointly performs a subtree scale move (a whole subtree is scaled up or down, keeping the topology fixed) on the species tree and on gene trees, all of which must be ultrametric. "; + + details.push_back( details_3 ); + + std::string details_4 = ""; + details_4 += "How this works: we pick a random node which is not the root.\nThen, we uniformly pick an age between the parent and the oldest sampled descendant.\nThe picked subtree is then scaled to this new age.\nAll gene-trees that are present in the population will be scaled accordingly."; + + details.push_back( details_4 ); + + return details; + } + + + /** + * Get an executable and instructive example. + * These examples should help the users to show how this function works but + * are also used to test if this function still works. + */ + std::string Move_SpeciesSubtreeScale::getHelpExample(void) const + { + // create an example as a single string variable. + std::string example = ""; + example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; + example += "dataFolder = \"simulatedTrees\" \n"; + example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; + example += "n_species <- 10\n"; + example += "n_genes <- 2\n"; + example += "n_alleles <- 3\n"; + example += "# we simulate an ultrametric species tree:\n"; + example += "# Species names:\n"; + example += "for (i in 1:n_species) {\n"; + example += " species[i] <- taxon(taxonName=\"Species_\"+i, speciesName=\"Species_\"+i)\n"; + example += "}\n"; + example += "spTree ~ dnBirthDeath(lambda=0.3, mu=0.2, rootAge=10, rho=1, samplingStrategy=\"uniform\", condition=\"nTaxa\", taxa=species)\n"; + example += "print(spTree)\n"; + example += "# let's pick a constant effective population size of 50:\n"; + example += "popSize <- 50\n"; + example += "# let's simulate gene trees now:\n"; + example += "# taxa names:\n"; + example += "for (g in 1:n_genes) {\n"; + example += " for (i in 1:n_species) {\n"; + example += " for (j in 1:n_alleles) {\n"; + example += " taxons[g][(i-1)*n_alleles+j] <- taxon(taxonName=\"Species_\"+i+\"_\"+j, speciesName=\"Species_\"+i)\n"; + example += " }\n"; + example += " }\n"; + example += " geneTrees[g] ~ dnMultiSpeciesCoalescent(speciesTree=spTree, Ne=popSize, taxa=taxons[g])\n"; + example += " print(geneTrees[g])\n"; + example += "}\n"; + example += "# We can save the species tree and the gene trees: \n"; + example += "write(speciesTree, filename=dataFolder+\"speciesTree\")\n"; + example += "# Saving the gene trees\n"; + example += "for (i in 1:(n_genes)) {\n"; + example += " write(geneTrees[i], filename=dataFolder+\"geneTree_\"+i+\".tree\")\n"; + example += "}\n"; + example += "# set my move index\n"; + example += "mi = 0\n"; + example += "move_species_subtree_scale = mvSpeciesSubtreeScale( speciesTree=psi, weight=5 )\n"; + example += "for (i in 1:n_genes) {\n"; + example += " move_species_subtree_scale.addGeneTreeVariable( geneTree[i] )\n"; + example += "}\n"; + example += "moves[++mi] = move_species_subtree_scale\n"; + example += "# We get a handle on our model.\n"; + example += "# We can use any node of our model as a handle, here we choose to use the topology.\n"; + example += "mymodel = model(spTree)\n"; + example += "# Monitors to check the progression of the program\n"; + example += "monitors[1] = mnScreen(printgen=10, sptree)\n"; + example += "# Here we use a plain MCMC. You could also set nruns=2 for a replicated analysis\n"; + example += "# or use mcmcmc with heated chains.\n"; + example += "mymcmc = mcmc(mymodel, monitors, moves, nruns=4)\n"; + example += "mymcmc.run(generations=1000)\n"; + example += "mymcmc.operatorSummary()\n"; + + return example; + } + + + /** + * Get some references/citations for this function + * + */ + std::vector Move_SpeciesSubtreeScale::getHelpReferences(void) const + { + // create an entry for each reference + std::vector references; + RevBayesCore::RbHelpReference ref = RevBayesCore::RbHelpReference(); + ref.setCitation("Guided Tree Topology Proposals for Bayesian Phylogenetic Inference. Sebastian Höhna, Alexei J. Drummond. Syst Biol (2012) 61 (1): 1-11."); + ref.setDoi("https://doi.org/10.1093/sysbio/syr074"); + ref.setUrl("https://academic.oup.com/sysbio/article-lookup/doi/10.1093/sysbio/syr074"); + + references.push_back(ref); + + RevBayesCore::RbHelpReference ref2 = RevBayesCore::RbHelpReference(); + ref2.setCitation("Algorithmic improvements to species delimitation and phylogeny estimation under the multispecies coalescent. Graham Jones. Journal of Mathematical Biology, 2016."); + ref2.setDoi("https://doi.org/10.1007/s00285-016-1034-0"); + ref2.setUrl("http://link.springer.com/article/10.1007/s00285-016-1034-0"); + + references.push_back(ref2); + + return references; + } + + + /** + * Get the names of similar and suggested other functions + */ + std::vector Move_SpeciesSubtreeScale::getHelpSeeAlso(void) const + { + // create an entry for each suggested function + std::vector see_also; + see_also.push_back( "mvSpeciesNodeTimeSlideUniform" ); + see_also.push_back( "mvSpeciesSubtreeScaleBeta" ); + see_also.push_back( "mvSpeciesNarrow" ); + see_also.push_back( "mvSpeciesTreeScale" ); + + return see_also; + } + + + /** + * Get the title of this help entry + */ + std::string Move_SpeciesSubtreeScale::getHelpTitle(void) const + { + // create a title variable + std::string title = "Subtree scale move on species tree and gene trees for multispecies coalescent models. "; + return title; + } diff --git a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.h b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.h index 3aa389111..c610eea9b 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.h +++ b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.h @@ -8,8 +8,8 @@ #include namespace RevLanguage { - - + + /** * The RevLanguage wrapper of the species-subtree-scale move. * @@ -26,11 +26,11 @@ namespace RevLanguage { * */ class Move_SpeciesSubtreeScale : public Move { - + public: - + Move_SpeciesSubtreeScale(void); //!< Default constructor - + // Basic utility functions virtual Move_SpeciesSubtreeScale* clone(void) const; //!< Clone object void constructInternalObject(void); //!< We construct the a new internal move. @@ -40,20 +40,29 @@ namespace RevLanguage { const MemberRules& getParameterRules(void) const; //!< Get member rules (const) virtual const TypeSpec& getTypeSpec(void) const; //!< Get language type of the object virtual void printValue(std::ostream& o) const; //!< Print value (for user) - + // Member method functions virtual RevPtr executeMethod(const std::string& name, const std::vector& args, bool &f); //!< Map member methods to internal functions - + protected: - + + + std::vector getHelpAuthor(void) const; //!< Get the author(s) of this function + std::vector getHelpDescription(void) const; //!< Get the description for this function + std::vector getHelpDetails(void) const; //!< Get the more detailed description of the function + std::string getHelpExample(void) const; //!< Get an executable and instructive example + std::vector getHelpReferences(void) const; //!< Get some references/citations for this function + std::vector getHelpSeeAlso(void) const; //!< Get suggested other functions + std::string getHelpTitle(void) const; //!< Get the title of this help entry + void setConstParameter(const std::string& name, const RevPtr &var); //!< Set member variable - + RevPtr speciesTree; RevPtr geneTrees; - + }; - + } #endif diff --git a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.cpp b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.cpp index a9db55c6c..f251fed24 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.cpp +++ b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.cpp @@ -21,7 +21,7 @@ using namespace RevLanguage; */ Move_SpeciesSubtreeScaleBeta::Move_SpeciesSubtreeScaleBeta() : Move() { - + // add method for call "addGeneTreeVariable" as a function ArgumentRules* addGeneTreeArgRules = new ArgumentRules(); addGeneTreeArgRules->push_back( new ArgumentRule( "geneTree" , TimeTree::getClassTypeSpec(), "A gene tree.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); @@ -38,7 +38,7 @@ Move_SpeciesSubtreeScaleBeta::Move_SpeciesSubtreeScaleBeta() : Move() */ Move_SpeciesSubtreeScaleBeta* Move_SpeciesSubtreeScaleBeta::clone(void) const { - + return new Move_SpeciesSubtreeScaleBeta(*this); } @@ -57,51 +57,51 @@ void Move_SpeciesSubtreeScaleBeta::constructInternalObject( void ) { // we free the memory first delete value; - + // now allocate a new sliding move double w = static_cast( weight->getRevObject() ).getValue(); double a = static_cast( alpha->getRevObject() ).getValue(); RevBayesCore::TypedDagNode* tmp = static_cast( speciesTree->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *st = static_cast *>( tmp ); - + RevBayesCore::Proposal *p = new RevBayesCore::SpeciesSubtreeScaleBetaProposal(st,a); - + bool t = static_cast( tune->getRevObject() ).getValue(); value = new RevBayesCore::MetropolisHastingsMove(p,w,t); - + } RevPtr Move_SpeciesSubtreeScaleBeta::executeMethod(const std::string& name, const std::vector& args, bool &found) { - + if ( name == "addGeneTreeVariable" ) { found = true; - + RevBayesCore::TypedDagNode* tmp = static_cast( args[0].getVariable()->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *gt = static_cast *>( tmp ); - + RevBayesCore::MetropolisHastingsMove *m = static_cast(this->value); RevBayesCore::SpeciesSubtreeScaleBetaProposal &p = static_cast( m->getProposal() ); p.addGeneTree( gt ); - + return NULL; } else if ( name == "removeGeneTreeVariable" ) { found = true; - + RevBayesCore::TypedDagNode* tmp = static_cast( args[0].getVariable()->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *gt = static_cast *>( tmp ); - + RevBayesCore::MetropolisHastingsMove *m = static_cast(this->value); RevBayesCore::SpeciesSubtreeScaleBetaProposal &p = static_cast( m->getProposal() ); p.removeGeneTree( gt ); - + return NULL; } - + return Move::executeMethod( name, args, found ); } @@ -113,9 +113,9 @@ RevPtr Move_SpeciesSubtreeScaleBeta::executeMethod(const std::strin */ const std::string& Move_SpeciesSubtreeScaleBeta::getClassType(void) { - + static std::string rev_type = "Move_SpeciesSubtreeScaleBeta"; - + return rev_type; } @@ -127,9 +127,9 @@ const std::string& Move_SpeciesSubtreeScaleBeta::getClassType(void) */ const TypeSpec& Move_SpeciesSubtreeScaleBeta::getClassTypeSpec(void) { - + static TypeSpec rev_type_spec = TypeSpec( getClassType(), new TypeSpec( Move::getClassTypeSpec() ) ); - + return rev_type_spec; } @@ -143,7 +143,7 @@ std::string Move_SpeciesSubtreeScaleBeta::getMoveName( void ) const { // create a constructor function name variable that is the same for all instance of this class std::string c_name = "SpeciesSubtreeScaleBeta"; - + return c_name; } @@ -158,23 +158,23 @@ std::string Move_SpeciesSubtreeScaleBeta::getMoveName( void ) const */ const MemberRules& Move_SpeciesSubtreeScaleBeta::getParameterRules(void) const { - + static MemberRules memberRules; static bool rules_set = false; - + if ( !rules_set ) { memberRules.push_back( new ArgumentRule( "speciesTree", TimeTree::getClassTypeSpec() , "The species tree on which this move operates on.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); memberRules.push_back( new ArgumentRule( "alpha" , RealPos::getClassTypeSpec() , "The concentration parameter.", ArgumentRule::BY_VALUE , ArgumentRule::ANY , new RealPos(10.0) ) ); memberRules.push_back( new ArgumentRule( "tune" , RlBoolean::getClassTypeSpec() , "Should we tune the concentration parameter during burnin?", ArgumentRule::BY_VALUE , ArgumentRule::ANY , new RlBoolean( true ) ) ); - + /* Inherit weight from Move, put it after variable */ const MemberRules& inheritedRules = Move::getParameterRules(); memberRules.insert( memberRules.end(), inheritedRules.begin(), inheritedRules.end() ); - + rules_set = true; } - + return memberRules; } @@ -186,9 +186,9 @@ const MemberRules& Move_SpeciesSubtreeScaleBeta::getParameterRules(void) const */ const TypeSpec& Move_SpeciesSubtreeScaleBeta::getTypeSpec( void ) const { - + static TypeSpec type_spec = getClassTypeSpec(); - + return type_spec; } @@ -198,7 +198,7 @@ const TypeSpec& Move_SpeciesSubtreeScaleBeta::getTypeSpec( void ) const */ void Move_SpeciesSubtreeScaleBeta::printValue(std::ostream &o) const { - + o << "SpeciesSubtreeScaleBeta("; if (speciesTree != NULL) { @@ -209,7 +209,7 @@ void Move_SpeciesSubtreeScaleBeta::printValue(std::ostream &o) const o << "?"; } o << ")"; - + } @@ -225,7 +225,7 @@ void Move_SpeciesSubtreeScaleBeta::printValue(std::ostream &o) const */ void Move_SpeciesSubtreeScaleBeta::setConstParameter(const std::string& name, const RevPtr &var) { - + if ( name == "speciesTree" ) { speciesTree = var; @@ -246,9 +246,181 @@ void Move_SpeciesSubtreeScaleBeta::setConstParameter(const std::string& name, co { Move::setConstParameter(name, var); } - + } + /** + * Get the author(s) of this function so they can receive credit (and blame) for it. + */ + std::vector Move_SpeciesSubtreeScaleBeta::getHelpAuthor(void) const + { + // create a vector of authors for this function + std::vector authors; + authors.push_back( "Sebastian Hoehna, Bastien Boussau" ); + + return authors; + } + + + /** + * Get the (brief) description for this function + */ + std::vector Move_SpeciesSubtreeScaleBeta::getHelpDescription(void) const + { + // create a variable for the description of the function + std::vector descriptions; + descriptions.push_back( "Makes a subtree scale move both in the species tree and in the gene trees that contain nodes of the relevant populations. Tree topologies are not altered. Uses a beta distribution to propose a new age value." ); + + return descriptions; + } + + + /** + * Get the more detailed description of the function + */ + std::vector Move_SpeciesSubtreeScaleBeta::getHelpDetails(void) const + { + // create a variable for the description of the function + std::vector details; + + std::string details_1 = ""; + details_1 += "The species tree must be ultrametric."; + + details.push_back( details_1 ); + + std::string details_2 = ""; + details_2 += "All the gene trees that evolved along the species tree according to some form of multispecies coalescent must be added to the move using the addGeneTreeVariable method. "; + + details.push_back( details_2 ); + + std::string details_3 = ""; + details_3 += "This move jointly performs a subtree scale move (a whole subtree is scaled up or down, keeping the topology fixed) on the species tree and on gene trees, all of which must be ultrametric. "; + + details.push_back( details_3 ); + + std::string details_4 = ""; + details_4 += "How this works: we pick a random node which is not the root.\nThen, we pick a new age between the parent and the oldest sampled descendant according to a beta distribution.\nThe picked subtree is then scaled to this new age.\nAll gene-trees that are present in the population will be scaled accordingly."; + + details.push_back( details_4 ); + + return details; + } + + + /** + * Get an executable and instructive example. + * These examples should help the users to show how this function works but + * are also used to test if this function still works. + */ + std::string Move_SpeciesSubtreeScaleBeta::getHelpExample(void) const + { + // create an example as a single string variable. + std::string example = ""; + example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; + example += "dataFolder = \"simulatedTrees\" \n"; + example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; + example += "n_species <- 10\n"; + example += "n_genes <- 2\n"; + example += "n_alleles <- 3\n"; + example += "# we simulate an ultrametric species tree:\n"; + example += "# Species names:\n"; + example += "for (i in 1:n_species) {\n"; + example += " species[i] <- taxon(taxonName=\"Species_\"+i, speciesName=\"Species_\"+i)\n"; + example += "}\n"; + example += "spTree ~ dnBirthDeath(lambda=0.3, mu=0.2, rootAge=10, rho=1, samplingStrategy=\"uniform\", condition=\"nTaxa\", taxa=species)\n"; + example += "print(spTree)\n"; + example += "# let's pick a constant effective population size of 50:\n"; + example += "popSize <- 50\n"; + example += "# let's simulate gene trees now:\n"; + example += "# taxa names:\n"; + example += "for (g in 1:n_genes) {\n"; + example += " for (i in 1:n_species) {\n"; + example += " for (j in 1:n_alleles) {\n"; + example += " taxons[g][(i-1)*n_alleles+j] <- taxon(taxonName=\"Species_\"+i+\"_\"+j, speciesName=\"Species_\"+i)\n"; + example += " }\n"; + example += " }\n"; + example += " geneTrees[g] ~ dnMultiSpeciesCoalescent(speciesTree=spTree, Ne=popSize, taxa=taxons[g])\n"; + example += " print(geneTrees[g])\n"; + example += "}\n"; + example += "# We can save the species tree and the gene trees: \n"; + example += "write(speciesTree, filename=dataFolder+\"speciesTree\")\n"; + example += "# Saving the gene trees\n"; + example += "for (i in 1:(n_genes)) {\n"; + example += " write(geneTrees[i], filename=dataFolder+\"geneTree_\"+i+\".tree\")\n"; + example += "}\n"; + example += "# set my move index\n"; + example += "mi = 0\n"; + example += "move_species_subtree_scale_beta = mvSpeciesSubtreeScaleBeta( speciesTree=psi, weight=5 )\n"; + example += "for (i in 1:n_genes) {\n"; + example += " move_species_subtree_scale_beta.addGeneTreeVariable( geneTree[i] )\n"; + example += "}\n"; + example += "moves[++mi] = move_species_subtree_scale_beta\n"; + example += "# We get a handle on our model.\n"; + example += "# We can use any node of our model as a handle, here we choose to use the topology.\n"; + example += "mymodel = model(spTree)\n"; + example += "# Monitors to check the progression of the program\n"; + example += "monitors[1] = mnScreen(printgen=10, sptree)\n"; + example += "# Here we use a plain MCMC. You could also set nruns=2 for a replicated analysis\n"; + example += "# or use mcmcmc with heated chains.\n"; + example += "mymcmc = mcmc(mymodel, monitors, moves, nruns=4)\n"; + example += "mymcmc.run(generations=1000)\n"; + example += "mymcmc.operatorSummary()\n"; + + return example; + } + + + /** + * Get some references/citations for this function + * + */ + std::vector Move_SpeciesSubtreeScaleBeta::getHelpReferences(void) const + { + // create an entry for each reference + std::vector references; + RevBayesCore::RbHelpReference ref = RevBayesCore::RbHelpReference(); + ref.setCitation("Guided Tree Topology Proposals for Bayesian Phylogenetic Inference. Sebastian Höhna, Alexei J. Drummond. Syst Biol (2012) 61 (1): 1-11."); + ref.setDoi("https://doi.org/10.1093/sysbio/syr074"); + ref.setUrl("https://academic.oup.com/sysbio/article-lookup/doi/10.1093/sysbio/syr074"); + + references.push_back(ref); + + RevBayesCore::RbHelpReference ref2 = RevBayesCore::RbHelpReference(); + ref2.setCitation("Algorithmic improvements to species delimitation and phylogeny estimation under the multispecies coalescent. Graham Jones. Journal of Mathematical Biology, 2016."); + ref2.setDoi("https://doi.org/10.1007/s00285-016-1034-0"); + ref2.setUrl("http://link.springer.com/article/10.1007/s00285-016-1034-0"); + + references.push_back(ref2); + + return references; + } + + + /** + * Get the names of similar and suggested other functions + */ + std::vector Move_SpeciesSubtreeScaleBeta::getHelpSeeAlso(void) const + { + // create an entry for each suggested function + std::vector see_also; + see_also.push_back( "mvSpeciesNodeTimeSlideUniform" ); + see_also.push_back( "mvSpeciesSubtreeScale" ); + see_also.push_back( "mvSpeciesNarrow" ); + see_also.push_back( "mvSpeciesTreeScale" ); + + return see_also; + } + + + /** + * Get the title of this help entry + */ + std::string Move_SpeciesSubtreeScaleBeta::getHelpTitle(void) const + { + // create a title variable + std::string title = "Subtree scale move on species tree and gene trees for multispecies coalescent models. "; + return title; + } diff --git a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.h b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.h index a6e455d32..2583c815f 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.h +++ b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.h @@ -8,12 +8,12 @@ #include namespace RevLanguage { - - + + /** - * The RevLanguage wrapper of the species-subtree-scale move. + * The RevLanguage wrapper of the species-subtree-scale-beta move. * - * The RevLanguage wrapper of the narrow-exchange move simply + * The RevLanguage wrapper of the species-subtree-scale-beta move simply * manages the interactions through the Rev with our core. * That is, the internal move object can be constructed and hooked up * in a DAG-nove (variable) that it works on. @@ -26,36 +26,45 @@ namespace RevLanguage { * */ class Move_SpeciesSubtreeScaleBeta : public Move { - + public: - + Move_SpeciesSubtreeScaleBeta(void); //!< Default constructor - + // Basic utility functions virtual Move_SpeciesSubtreeScaleBeta* clone(void) const; //!< Clone object - void constructInternalObject(void); //!< We construct the a new internal SlidingMove. + void constructInternalObject(void); //!< We construct a new internal SlidingMove. static const std::string& getClassType(void); //!< Get Rev type static const TypeSpec& getClassTypeSpec(void); //!< Get class type spec std::string getMoveName(void) const; //!< Get the name used for the constructor function in Rev. const MemberRules& getParameterRules(void) const; //!< Get member rules (const) virtual const TypeSpec& getTypeSpec(void) const; //!< Get language type of the object virtual void printValue(std::ostream& o) const; //!< Print value (for user) - + // Member method functions virtual RevPtr executeMethod(const std::string& name, const std::vector& args, bool &f); //!< Map member methods to internal functions - + protected: - + + + std::vector getHelpAuthor(void) const; //!< Get the author(s) of this function + std::vector getHelpDescription(void) const; //!< Get the description for this function + std::vector getHelpDetails(void) const; //!< Get the more detailed description of the function + std::string getHelpExample(void) const; //!< Get an executable and instructive example + std::vector getHelpReferences(void) const; //!< Get some references/citations for this function + std::vector getHelpSeeAlso(void) const; //!< Get suggested other functions + std::string getHelpTitle(void) const; //!< Get the title of this help entry + void setConstParameter(const std::string& name, const RevPtr &var); //!< Set member variable - + RevPtr speciesTree; RevPtr geneTrees; RevPtr alpha; RevPtr tune; //!< If autotuning should be used. - + }; - + } #endif diff --git a/src/revlanguage/moves/tree/Move_SpeciesTreeScale.cpp b/src/revlanguage/moves/tree/Move_SpeciesTreeScale.cpp index 6b5ba3a62..2d8bf2ba5 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesTreeScale.cpp +++ b/src/revlanguage/moves/tree/Move_SpeciesTreeScale.cpp @@ -21,7 +21,7 @@ using namespace RevLanguage; */ Move_SpeciesTreeScale::Move_SpeciesTreeScale() : Move() { - + // add method for call "addGeneTreeVariable" as a function ArgumentRules* addGeneTreeArgRules = new ArgumentRules(); addGeneTreeArgRules->push_back( new ArgumentRule( "geneTree" , TimeTree::getClassTypeSpec(), "A gene tree variable.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); @@ -37,7 +37,7 @@ Move_SpeciesTreeScale::Move_SpeciesTreeScale() : Move() */ Move_SpeciesTreeScale* Move_SpeciesTreeScale::clone(void) const { - + return new Move_SpeciesTreeScale(*this); } @@ -56,53 +56,53 @@ void Move_SpeciesTreeScale::constructInternalObject( void ) { // we free the memory first delete value; - + // now allocate a new move double w = static_cast( weight->getRevObject() ).getValue(); RevBayesCore::TypedDagNode* tmp = static_cast( speciesTree->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *st = static_cast *>( tmp ); - + RevBayesCore::TypedDagNode* tmp_a = static_cast( rootAge->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *ra = static_cast *>( tmp_a ); double d = static_cast( delta->getRevObject() ).getValue(); bool tune = static_cast( tuning->getRevObject() ).getValue(); - + RevBayesCore::Proposal *p = new RevBayesCore::SpeciesTreeScaleProposal(st,ra,d); value = new RevBayesCore::MetropolisHastingsMove(p, w, tune); - + } RevPtr Move_SpeciesTreeScale::executeMethod(const std::string& name, const std::vector& args, bool &found) { - + if ( name == "addGeneTreeVariable" ) { found = true; - + RevBayesCore::TypedDagNode* tmp = static_cast( args[0].getVariable()->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *gt = static_cast *>( tmp ); - + RevBayesCore::MetropolisHastingsMove *m = static_cast(this->value); RevBayesCore::SpeciesTreeScaleProposal &p = static_cast( m->getProposal() ); p.addGeneTree( gt ); - + return NULL; } else if ( name == "removeGeneTreeVariable" ) { found = true; - + RevBayesCore::TypedDagNode* tmp = static_cast( args[0].getVariable()->getRevObject() ).getDagNode(); RevBayesCore::StochasticNode *gt = static_cast *>( tmp ); - + RevBayesCore::MetropolisHastingsMove *m = static_cast(this->value); RevBayesCore::SpeciesTreeScaleProposal &p = static_cast( m->getProposal() ); p.removeGeneTree( gt ); - + return NULL; } - + return Move::executeMethod( name, args, found ); } @@ -114,9 +114,9 @@ RevPtr Move_SpeciesTreeScale::executeMethod(const std::string& name */ const std::string& Move_SpeciesTreeScale::getClassType(void) { - + static std::string rev_type = "Move_SpeciesTreeScale"; - + return rev_type; } @@ -128,9 +128,9 @@ const std::string& Move_SpeciesTreeScale::getClassType(void) */ const TypeSpec& Move_SpeciesTreeScale::getClassTypeSpec(void) { - + static TypeSpec rev_type_spec = TypeSpec( getClassType(), new TypeSpec( Move::getClassTypeSpec() ) ); - + return rev_type_spec; } @@ -144,7 +144,7 @@ std::string Move_SpeciesTreeScale::getMoveName( void ) const { // create a constructor function name variable that is the same for all instance of this class std::string c_name = "SpeciesTreeScale"; - + return c_name; } @@ -159,24 +159,24 @@ std::string Move_SpeciesTreeScale::getMoveName( void ) const */ const MemberRules& Move_SpeciesTreeScale::getParameterRules(void) const { - + static MemberRules memberRules; static bool rules_set = false; - + if ( !rules_set ) { memberRules.push_back( new ArgumentRule( "speciesTree", TimeTree::getClassTypeSpec() , "The species tree on which this move operates.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); memberRules.push_back( new ArgumentRule( "rootAge" , RealPos::getClassTypeSpec() , "The root age variable.", ArgumentRule::BY_REFERENCE, ArgumentRule::STOCHASTIC ) ); memberRules.push_back( new ArgumentRule( "delta" , RealPos::getClassTypeSpec() , "The strength of the proposal", ArgumentRule::BY_VALUE , ArgumentRule::ANY , new RealPos( 1.0 ) ) ); memberRules.push_back( new ArgumentRule( "tune" , RlBoolean::getClassTypeSpec() , "Should we tune the strength during burnin?", ArgumentRule::BY_VALUE , ArgumentRule::ANY , new RlBoolean( true ) ) ); - + /* Inherit weight from Move, put it after variable */ const MemberRules& inheritedRules = Move::getParameterRules(); memberRules.insert( memberRules.end(), inheritedRules.begin(), inheritedRules.end() ); - + rules_set = true; } - + return memberRules; } @@ -188,9 +188,9 @@ const MemberRules& Move_SpeciesTreeScale::getParameterRules(void) const */ const TypeSpec& Move_SpeciesTreeScale::getTypeSpec( void ) const { - + static TypeSpec type_spec = getClassTypeSpec(); - + return type_spec; } @@ -200,7 +200,7 @@ const TypeSpec& Move_SpeciesTreeScale::getTypeSpec( void ) const */ void Move_SpeciesTreeScale::printValue(std::ostream &o) const { - + o << "SpeciesTreeScale("; if (speciesTree != NULL) { @@ -211,7 +211,7 @@ void Move_SpeciesTreeScale::printValue(std::ostream &o) const o << "?"; } o << ")"; - + } @@ -227,7 +227,7 @@ void Move_SpeciesTreeScale::printValue(std::ostream &o) const */ void Move_SpeciesTreeScale::setConstParameter(const std::string& name, const RevPtr &var) { - + if ( name == "speciesTree" ) { speciesTree = var; @@ -252,9 +252,178 @@ void Move_SpeciesTreeScale::setConstParameter(const std::string& name, const Rev { Move::setConstParameter(name, var); } - + } + /** + * Get the author(s) of this function so they can receive credit (and blame) for it. + */ + std::vector Move_SpeciesTreeScale::getHelpAuthor(void) const + { + // create a vector of authors for this function + std::vector authors; + authors.push_back( "Sebastian Hoehna, Bastien Boussau" ); + + return authors; + } + + + /** + * Get the (brief) description for this function + */ + std::vector Move_SpeciesTreeScale::getHelpDescription(void) const + { + // create a variable for the description of the function + std::vector descriptions; + descriptions.push_back( "Makes a tree scale move both in the species tree and in the gene trees. Tree topologies are not altered." ); + + return descriptions; + } + + + /** + * Get the more detailed description of the function + */ + std::vector Move_SpeciesTreeScale::getHelpDetails(void) const + { + // create a variable for the description of the function + std::vector details; + + std::string details_1 = ""; + details_1 += "The species tree must be ultrametric."; + + details.push_back( details_1 ); + + std::string details_2 = ""; + details_2 += "All the gene trees that evolved along the species tree according to some form of multispecies coalescent must be added to the move using the addGeneTreeVariable method. "; + + details.push_back( details_2 ); + + std::string details_3 = ""; + details_3 += "This move jointly performs a tree scale move (the entire tree is scaled up or down, keeping the topology fixed) on the species tree and on gene trees, all of which must be ultrametric. "; + + details.push_back( details_3 ); + + return details; + } + + + /** + * Get an executable and instructive example. + * These examples should help the users to show how this function works but + * are also used to test if this function still works. + */ + std::string Move_SpeciesTreeScale::getHelpExample(void) const + { + // create an example as a single string variable. + std::string example = ""; + example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; + example += "dataFolder = \"simulatedTrees\" \n"; + example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; + example += "n_species <- 10\n"; + example += "n_genes <- 2\n"; + example += "n_alleles <- 3\n"; + example += "# we simulate an ultrametric species tree:\n"; + example += "# Species names:\n"; + example += "for (i in 1:n_species) {\n"; + example += " species[i] <- taxon(taxonName=\"Species_\"+i, speciesName=\"Species_\"+i)\n"; + example += "}\n"; + example += "root ~ dnNormal(mean=75,sd=2.5,min=0.0, max=Inf)\n"; + example += "spTree ~ dnBirthDeath(lambda=0.3, mu=0.2, rootAge=root, rho=1, samplingStrategy=\"uniform\", condition=\"nTaxa\", taxa=species)\n"; + example += "print(spTree)\n"; + example += "# let's pick a constant effective population size of 50:\n"; + example += "popSize <- 50\n"; + example += "# let's simulate gene trees now:\n"; + example += "# taxa names:\n"; + example += "for (g in 1:n_genes) {\n"; + example += " for (i in 1:n_species) {\n"; + example += " for (j in 1:n_alleles) {\n"; + example += " taxons[g][(i-1)*n_alleles+j] <- taxon(taxonName=\"Species_\"+i+\"_\"+j, speciesName=\"Species_\"+i)\n"; + example += " }\n"; + example += " }\n"; + example += " geneTrees[g] ~ dnMultiSpeciesCoalescent(speciesTree=spTree, Ne=popSize, taxa=taxons[g])\n"; + example += " print(geneTrees[g])\n"; + example += "}\n"; + example += "# We can save the species tree and the gene trees: \n"; + example += "write(speciesTree, filename=dataFolder+\"speciesTree\")\n"; + example += "# Saving the gene trees\n"; + example += "for (i in 1:(n_genes)) {\n"; + example += " write(geneTrees[i], filename=dataFolder+\"geneTree_\"+i+\".tree\")\n"; + example += "}\n"; + example += "# set my move index\n"; + example += "mi = 0\n"; + example += "move_species_tree_scale = mvSpeciesTreeScale( speciesTree=psi, root=root, weight=5 )\n"; + example += "for (i in 1:n_genes) {\n"; + example += " move_species_tree_scale.addGeneTreeVariable( geneTree[i] )\n"; + example += "}\n"; + example += "moves[++mi] = move_species_tree_scale\n"; + example += "# We get a handle on our model.\n"; + example += "# We can use any node of our model as a handle, here we choose to use the topology.\n"; + example += "mymodel = model(spTree)\n"; + example += "# Monitors to check the progression of the program\n"; + example += "monitors[1] = mnScreen(printgen=10, sptree)\n"; + example += "# Here we use a plain MCMC. You could also set nruns=2 for a replicated analysis\n"; + example += "# or use mcmcmc with heated chains.\n"; + example += "mymcmc = mcmc(mymodel, monitors, moves, nruns=4)\n"; + example += "mymcmc.run(generations=1000)\n"; + example += "mymcmc.operatorSummary()\n"; + + return example; + } + + + /** + * Get some references/citations for this function + * + */ + std::vector Move_SpeciesTreeScale::getHelpReferences(void) const + { + // create an entry for each reference + std::vector references; + RevBayesCore::RbHelpReference ref = RevBayesCore::RbHelpReference(); + ref.setCitation("Guided Tree Topology Proposals for Bayesian Phylogenetic Inference. Sebastian Höhna, Alexei J. Drummond. Syst Biol (2012) 61 (1): 1-11."); + ref.setDoi("https://doi.org/10.1093/sysbio/syr074"); + ref.setUrl("https://academic.oup.com/sysbio/article-lookup/doi/10.1093/sysbio/syr074"); + + references.push_back(ref); + + RevBayesCore::RbHelpReference ref2 = RevBayesCore::RbHelpReference(); + ref2.setCitation("Algorithmic improvements to species delimitation and phylogeny estimation under the multispecies coalescent. Graham Jones. Journal of Mathematical Biology, 2016."); + ref2.setDoi("https://doi.org/10.1007/s00285-016-1034-0"); + ref2.setUrl("http://link.springer.com/article/10.1007/s00285-016-1034-0"); + + references.push_back(ref2); + + return references; + } + + + /** + * Get the names of similar and suggested other functions + */ + std::vector Move_SpeciesTreeScale::getHelpSeeAlso(void) const + { + // create an entry for each suggested function + std::vector see_also; + see_also.push_back( "mvSpeciesNodeTimeSlideUniform" ); + see_also.push_back( "mvSpeciesSubtreeScaleBeta" ); + see_also.push_back( "mvSpeciesNarrow" ); + see_also.push_back( "mvSpeciesSubtreeScale" ); + + return see_also; + } + + + /** + * Get the title of this help entry + */ + std::string Move_SpeciesTreeScale::getHelpTitle(void) const + { + // create a title variable + std::string title = "Tree scale move on species tree and gene trees for multispecies coalescent models. "; + + return title; + } diff --git a/src/revlanguage/moves/tree/Move_SpeciesTreeScale.h b/src/revlanguage/moves/tree/Move_SpeciesTreeScale.h index c06cbc134..110cdedbb 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesTreeScale.h +++ b/src/revlanguage/moves/tree/Move_SpeciesTreeScale.h @@ -8,12 +8,12 @@ #include namespace RevLanguage { - - + + /** - * The RevLanguage wrapper of the species-subtree-scale move. + * The RevLanguage wrapper of the species-tree-scale move. * - * The RevLanguage wrapper of the narrow-exchange move simply + * The RevLanguage wrapper of the species-tree-scale move simply * manages the interactions through the Rev with our core. * That is, the internal move object can be constructed and hooked up * in a DAG-nove (variable) that it works on. @@ -26,11 +26,11 @@ namespace RevLanguage { * */ class Move_SpeciesTreeScale : public Move { - + public: - + Move_SpeciesTreeScale(void); //!< Default constructor - + // Basic utility functions virtual Move_SpeciesTreeScale* clone(void) const; //!< Clone object void constructInternalObject(void); //!< We construct the a new internal SlidingMove. @@ -40,23 +40,31 @@ namespace RevLanguage { const MemberRules& getParameterRules(void) const; //!< Get member rules (const) virtual const TypeSpec& getTypeSpec(void) const; //!< Get language type of the object virtual void printValue(std::ostream& o) const; //!< Print value (for user) - + // Member method functions virtual RevPtr executeMethod(const std::string& name, const std::vector& args, bool &f); //!< Map member methods to internal functions - + protected: - + + std::vector getHelpAuthor(void) const; //!< Get the author(s) of this function + std::vector getHelpDescription(void) const; //!< Get the description for this function + std::vector getHelpDetails(void) const; //!< Get the more detailed description of the function + std::string getHelpExample(void) const; //!< Get an executable and instructive example + std::vector getHelpReferences(void) const; //!< Get some references/citations for this function + std::vector getHelpSeeAlso(void) const; //!< Get suggested other functions + std::string getHelpTitle(void) const; //!< Get the title of this help entry + void setConstParameter(const std::string& name, const RevPtr &var); //!< Set member variable - + RevPtr speciesTree; RevPtr geneTrees; RevPtr rootAge; RevPtr delta; RevPtr tuning; - + }; - + } #endif From 613e000f4689c777fc5ccc21fd7b311abfa62832 Mon Sep 17 00:00:00 2001 From: boussau Date: Fri, 4 Aug 2017 16:28:17 +0200 Subject: [PATCH 2/2] Corrected a small typo in the example code for MSC-related moves and distributions. --- .../phylogenetics/tree/Dist_constPopMultispCoal.cpp | 2 +- .../tree/Dist_multispeciesCoalescentInverseGammaPrior.cpp | 2 +- .../tree/Dist_multispeciesCoalescentUniformPrior.cpp | 2 +- src/revlanguage/moves/tree/Move_SpeciesNarrowExchange.cpp | 2 +- src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.cpp | 2 +- src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.cpp | 2 +- src/revlanguage/moves/tree/Move_SpeciesTreeScale.cpp | 2 +- 7 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/revlanguage/distributions/phylogenetics/tree/Dist_constPopMultispCoal.cpp b/src/revlanguage/distributions/phylogenetics/tree/Dist_constPopMultispCoal.cpp index 9ced1eac6..d30af7c79 100644 --- a/src/revlanguage/distributions/phylogenetics/tree/Dist_constPopMultispCoal.cpp +++ b/src/revlanguage/distributions/phylogenetics/tree/Dist_constPopMultispCoal.cpp @@ -269,7 +269,7 @@ const TypeSpec& Dist_constPopMultispCoal::getTypeSpec( void ) const // create an example as a single string variable. std::string example = ""; example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; - example += "dataFolder = \"simulatedTrees\" \n"; + example += "dataFolder = \"simulatedTrees/\" \n"; example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; example += "n_species <- 10\n"; example += "n_genes <- 2\n"; diff --git a/src/revlanguage/distributions/phylogenetics/tree/Dist_multispeciesCoalescentInverseGammaPrior.cpp b/src/revlanguage/distributions/phylogenetics/tree/Dist_multispeciesCoalescentInverseGammaPrior.cpp index 51578c1c6..6f6edce4f 100644 --- a/src/revlanguage/distributions/phylogenetics/tree/Dist_multispeciesCoalescentInverseGammaPrior.cpp +++ b/src/revlanguage/distributions/phylogenetics/tree/Dist_multispeciesCoalescentInverseGammaPrior.cpp @@ -271,7 +271,7 @@ std::string Dist_multispeciesCoalescentInverseGammaPrior::getHelpExample(void) c // create an example as a single string variable. std::string example = ""; example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; - example += "dataFolder = \"simulatedTrees\" \n"; + example += "dataFolder = \"simulatedTrees/\" \n"; example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; example += "n_species <- 10\n"; example += "n_genes <- 2\n"; diff --git a/src/revlanguage/distributions/phylogenetics/tree/Dist_multispeciesCoalescentUniformPrior.cpp b/src/revlanguage/distributions/phylogenetics/tree/Dist_multispeciesCoalescentUniformPrior.cpp index eb2817fc8..c8e934f49 100644 --- a/src/revlanguage/distributions/phylogenetics/tree/Dist_multispeciesCoalescentUniformPrior.cpp +++ b/src/revlanguage/distributions/phylogenetics/tree/Dist_multispeciesCoalescentUniformPrior.cpp @@ -250,7 +250,7 @@ const TypeSpec& Dist_multispeciesCoalescentUniformPrior::getTypeSpec( void ) con // create an example as a single string variable. std::string example = ""; example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; - example += "dataFolder = \"simulatedTrees\" \n"; + example += "dataFolder = \"simulatedTrees/\" \n"; example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; example += "n_species <- 10\n"; example += "n_genes <- 2\n"; diff --git a/src/revlanguage/moves/tree/Move_SpeciesNarrowExchange.cpp b/src/revlanguage/moves/tree/Move_SpeciesNarrowExchange.cpp index 5456810bc..ad9cba26f 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesNarrowExchange.cpp +++ b/src/revlanguage/moves/tree/Move_SpeciesNarrowExchange.cpp @@ -299,7 +299,7 @@ void Move_SpeciesNarrowExchange::setConstParameter(const std::string& name, cons // create an example as a single string variable. std::string example = ""; example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; - example += "dataFolder = \"simulatedTrees\" \n"; + example += "dataFolder = \"simulatedTrees/\" \n"; example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; example += "n_species <- 10\n"; example += "n_genes <- 2\n"; diff --git a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.cpp b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.cpp index 6c5b85864..ef18ade0b 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.cpp +++ b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScale.cpp @@ -301,7 +301,7 @@ void Move_SpeciesSubtreeScale::setConstParameter(const std::string& name, const // create an example as a single string variable. std::string example = ""; example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; - example += "dataFolder = \"simulatedTrees\" \n"; + example += "dataFolder = \"simulatedTrees/\" \n"; example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; example += "n_species <- 10\n"; example += "n_genes <- 2\n"; diff --git a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.cpp b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.cpp index f251fed24..9d9f5059a 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.cpp +++ b/src/revlanguage/moves/tree/Move_SpeciesSubtreeScaleBeta.cpp @@ -319,7 +319,7 @@ void Move_SpeciesSubtreeScaleBeta::setConstParameter(const std::string& name, co // create an example as a single string variable. std::string example = ""; example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; - example += "dataFolder = \"simulatedTrees\" \n"; + example += "dataFolder = \"simulatedTrees/\" \n"; example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; example += "n_species <- 10\n"; example += "n_genes <- 2\n"; diff --git a/src/revlanguage/moves/tree/Move_SpeciesTreeScale.cpp b/src/revlanguage/moves/tree/Move_SpeciesTreeScale.cpp index 2d8bf2ba5..b3786d27f 100644 --- a/src/revlanguage/moves/tree/Move_SpeciesTreeScale.cpp +++ b/src/revlanguage/moves/tree/Move_SpeciesTreeScale.cpp @@ -321,7 +321,7 @@ void Move_SpeciesTreeScale::setConstParameter(const std::string& name, const Rev // create an example as a single string variable. std::string example = ""; example += "# We are going to save the trees we simulate in the folder simulatedTrees:\n"; - example += "dataFolder = \"simulatedTrees\" \n"; + example += "dataFolder = \"simulatedTrees/\" \n"; example += "# Let’s simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:\n"; example += "n_species <- 10\n"; example += "n_genes <- 2\n";