From 18d1e745aeb72325fa2aa911bf33e1797ac0b163 Mon Sep 17 00:00:00 2001 From: Michael Landis Date: Wed, 1 Aug 2018 14:20:49 -0400 Subject: [PATCH 1/3] update host switch rate mod --- .../ratemap/HostSwitchRateModifier.cpp | 100 +++++++++++++----- .../ratemap/HostSwitchRateModifier.h | 2 + .../HostSwitchRateModifierFunction.cpp | 2 +- 3 files changed, 75 insertions(+), 29 deletions(-) diff --git a/src/core/datatypes/phylogenetics/ratemap/HostSwitchRateModifier.cpp b/src/core/datatypes/phylogenetics/ratemap/HostSwitchRateModifier.cpp index 2bbcde554..c11d2eb9a 100644 --- a/src/core/datatypes/phylogenetics/ratemap/HostSwitchRateModifier.cpp +++ b/src/core/datatypes/phylogenetics/ratemap/HostSwitchRateModifier.cpp @@ -8,9 +8,10 @@ using namespace RevBayesCore; -HostSwitchRateModifier::HostSwitchRateModifier(size_t ns, size_t nc) : CharacterHistoryRateModifier(ns, nc), +HostSwitchRateModifier::HostSwitchRateModifier(size_t ns, size_t nc) : CharacterHistoryRateModifier(3, nc), scale( 1.0 ), num_branches( 0 ) + { ; } @@ -41,51 +42,88 @@ HostSwitchRateModifier& HostSwitchRateModifier::assign(const Assignable &m) } } -//computeRateMultiplier(std::vector currState, CharacterEventDiscrete* newState, std::vector > sites_with_states, double age=0.0) -double HostSwitchRateModifier::computeRateMultiplier(std::vector currState, CharacterEventDiscrete* newState, double age) + +double HostSwitchRateModifier::computeRateMultiplier(std::vector currState, CharacterEventDiscrete* newState, std::vector > sites_with_states, double age) { + // which character will change? - size_t index = newState->getSiteIndex(); + size_t to_index = newState->getSiteIndex(); // what is the state change? - size_t from_state = static_cast(currState[index])->getState(); + size_t from_state = static_cast(currState[to_index])->getState(); size_t to_state = newState->getState(); + double r = 1.0; + // loss event (independent of other hosts) if (from_state > to_state) { - return 1.0; + if (sites_with_states[1].size() == 1 && sites_with_states[2].size() == 0) + { + // cannot enter the null range (conditions on survival) + r = 0.0; + return r; + } + else + { + r = 1.0; + return r; + } } - - // gain event - double scaler_value = scale[ to_state - 1 ]; - - // if the gain event level's scaling factor equals zero, then there's no effect - if ( scaler_value == 0.0 ) + else { - return 1.0; + // gain event + double scaler_value = scale[ to_state - 1 ]; + + // if the gain event level's scaling factor equals zero, then there's no effect + if ( scaler_value == 0.0 ) + { + return 1.0; + } + + // sum of phylo.distance-scaled rates + double delta = 0.0; + size_t n_on = 0; + for (size_t from_index = 0; from_index < this->num_characters; from_index++) + { + size_t s = static_cast(currState[from_index])->getState(); + if (s != 0) { + delta += distance[from_index][to_index]; + n_on += 1; + } else { + ; // do nothing + } + } + + double delta_mean = delta / n_on; + r = std::exp( -scaler_value * delta_mean ); + } + return r; +} + +double HostSwitchRateModifier::computeRateMultiplier(std::vector currState, CharacterEventDiscrete* newState, double age) +{ + std::vector > sites_with_states(num_states); + for (size_t i = 0; i < currState.size(); i++) + { + sites_with_states[ static_cast(currState[i])->getState() ].insert(i); } - - // distance-scaled gain event - double r = 0.0; - - // sum of phylo.distance-scaled rates - for (size_t i = 0; i < this->num_characters; i++) + return computeRateMultiplier(currState, newState, sites_with_states, age); +} + +double HostSwitchRateModifier::computeRateMultiplier(std::vector currState, CharacterEventDiscrete* newState, std::vector counts, double age) +{ + std::vector > sites_with_states(num_states); + for (size_t i = 0; i < currState.size(); i++) { - size_t s = static_cast(currState[i])->getState(); - if (s != 0) { - double v = std::exp( -scaler_value * distance[i][index] ); - r += v; -// std::cout << i << " -> " << index << " , " << to_state << " = " << v << "\n"; - } else { - ; - } + sites_with_states[ static_cast(currState[i])->getState() ].insert(i); } - return r; + return computeRateMultiplier(currState, newState, sites_with_states, age); } + double HostSwitchRateModifier::computeSiteRateMultiplier(const TopologyNode& node, CharacterEvent* currState, CharacterEvent* newState, double age) { return 1.0; @@ -113,6 +151,12 @@ void HostSwitchRateModifier::setTree(const RevBayesCore::Tree &t) num_branches = tau.getNumberOfNodes() - 1; distance = *RevBayesCore::TreeUtilities::getDistanceMatrix ( tau ); + double max_distance = 2 * tau.getRoot().getAge(); + for (size_t i = 0; i < distance.size(); i++) { + for (size_t j = 0; j < distance[i].size(); j++) { + distance[i][j] /= max_distance; + } + } } void HostSwitchRateModifier::setScale(const std::vector& s) diff --git a/src/core/datatypes/phylogenetics/ratemap/HostSwitchRateModifier.h b/src/core/datatypes/phylogenetics/ratemap/HostSwitchRateModifier.h index 11ed64fec..1e6183e4d 100644 --- a/src/core/datatypes/phylogenetics/ratemap/HostSwitchRateModifier.h +++ b/src/core/datatypes/phylogenetics/ratemap/HostSwitchRateModifier.h @@ -33,6 +33,8 @@ namespace RevBayesCore HostSwitchRateModifier& assign(const Assignable &m); double computeRateMultiplier(std::vector currState, CharacterEventDiscrete* newState, double age=0.0); + double computeRateMultiplier(std::vector currState, CharacterEventDiscrete* newState, std::vector counts, double age=0.0); + double computeRateMultiplier(std::vector currState, CharacterEventDiscrete* newState, std::vector > sites_with_states, double age=0.0); double computeSiteRateMultiplier(const TopologyNode& node, CharacterEvent* currState, CharacterEvent* newState, double age=0.0); double computeSiteRateMultiplier(const TopologyNode& node, unsigned currState, unsigned newState, unsigned charIdx=0, double age=0.0); diff --git a/src/core/functions/phylogenetics/ratemap/HostSwitchRateModifierFunction.cpp b/src/core/functions/phylogenetics/ratemap/HostSwitchRateModifierFunction.cpp index b7f2a2b41..806ce2c27 100644 --- a/src/core/functions/phylogenetics/ratemap/HostSwitchRateModifierFunction.cpp +++ b/src/core/functions/phylogenetics/ratemap/HostSwitchRateModifierFunction.cpp @@ -17,7 +17,7 @@ using namespace RevBayesCore; HostSwitchRateModifierFunction::HostSwitchRateModifierFunction(const TypedDagNode* t, const TypedDagNode >* s) : -TypedFunction( new HostSwitchRateModifier(0, t->getValue().getNumberOfTips() ) ), +TypedFunction( new HostSwitchRateModifier(3, t->getValue().getNumberOfTips() ) ), tau(t), scale(s) { From 01979ab4f1c7baa903ac2155120d7ae7dcf33849 Mon Sep 17 00:00:00 2001 From: Jiansi Gao Date: Wed, 1 Aug 2018 18:53:43 -0700 Subject: [PATCH 2/3] allowing the probability (density) of node ages (unlabeled history) to be correctly computed (rather than returning nan) under a ssbdp when the sampling probability at the present (rho) is 0 and there is no extant sample (in which scenario the probability should still be positive; this is reasonable for fast-evolving pathogens; see Stadler et al. 2012 MBE) --- .../tree/AbstractRootedTreeDistribution.cpp | 11 ++++++++++- .../ConstantRateSerialSampledBirthDeathProcess.cpp | 7 +++++-- ...iecewiseConstantSerialSampledBirthDeathProcess.cpp | 5 ++++- 3 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/core/distributions/phylogenetics/tree/AbstractRootedTreeDistribution.cpp b/src/core/distributions/phylogenetics/tree/AbstractRootedTreeDistribution.cpp index 3186503cc..5f61622a6 100644 --- a/src/core/distributions/phylogenetics/tree/AbstractRootedTreeDistribution.cpp +++ b/src/core/distributions/phylogenetics/tree/AbstractRootedTreeDistribution.cpp @@ -660,8 +660,17 @@ void AbstractRootedTreeDistribution::simulateTree( void ) ra = rng->uniform01() * ( max_age - max_node_age ) + max_node_age; } + + double min_node_age = max_node_age; + for (size_t j = 0; j < nodes.size(); ++j) + { + if ( nodes[j]->getAge() < min_node_age ) + { + min_node_age = nodes[j]->getAge(); + } + } - simulateClade(nodes, ra, 0.0); + simulateClade(nodes, ra, min_node_age); TopologyNode *root = nodes[0]; diff --git a/src/core/distributions/phylogenetics/tree/birthdeath/ConstantRateSerialSampledBirthDeathProcess.cpp b/src/core/distributions/phylogenetics/tree/birthdeath/ConstantRateSerialSampledBirthDeathProcess.cpp index 8b36867c6..bb39d1083 100644 --- a/src/core/distributions/phylogenetics/tree/birthdeath/ConstantRateSerialSampledBirthDeathProcess.cpp +++ b/src/core/distributions/phylogenetics/tree/birthdeath/ConstantRateSerialSampledBirthDeathProcess.cpp @@ -166,8 +166,11 @@ double ConstantRateSerialSampledBirthDeathProcess::computeLnProbabilityTimes( vo } // add the log probability for sampling the extant taxa - lnProbTimes += num_extant_taxa * log( 4.0 * sampling_prob ); - + if (num_extant_taxa > 0) + { + lnProbTimes += num_extant_taxa * log( 4.0 * sampling_prob ); + } + // add the log probability of the initial sequences lnProbTimes += -lnQ(process_time, c1, c2) * num_initial_lineages; diff --git a/src/core/distributions/phylogenetics/tree/birthdeath/PiecewiseConstantSerialSampledBirthDeathProcess.cpp b/src/core/distributions/phylogenetics/tree/birthdeath/PiecewiseConstantSerialSampledBirthDeathProcess.cpp index 037cdfe40..e77c050bc 100644 --- a/src/core/distributions/phylogenetics/tree/birthdeath/PiecewiseConstantSerialSampledBirthDeathProcess.cpp +++ b/src/core/distributions/phylogenetics/tree/birthdeath/PiecewiseConstantSerialSampledBirthDeathProcess.cpp @@ -445,7 +445,10 @@ double PiecewiseConstantSerialSampledBirthDeathProcess::computeLnProbabilityTime else { // add the extant tip age term - lnProbTimes += num_extant_taxa * log( homogeneous_rho->getValue() ); + if (num_extant_taxa > 0) + { + lnProbTimes += num_extant_taxa * log( homogeneous_rho->getValue() ); + } } // add the sampled ancestor age terms From c2f106d768ce47f72aa60c93439da3f57f574f5f Mon Sep 17 00:00:00 2001 From: Sebastian Hoehna Date: Thu, 2 Aug 2018 09:30:07 +0200 Subject: [PATCH 3/3] Breaking up type registration for windows compilation --- src/revlanguage/workspace/RbRegister_Type.cpp | 17 -- .../workspace/RbRegister_VectorType.cpp | 164 ++++++++++++++++++ src/revlanguage/workspace/Workspace.cpp | 1 + src/revlanguage/workspace/Workspace.h | 3 +- 4 files changed, 167 insertions(+), 18 deletions(-) create mode 100644 src/revlanguage/workspace/RbRegister_VectorType.cpp diff --git a/src/revlanguage/workspace/RbRegister_Type.cpp b/src/revlanguage/workspace/RbRegister_Type.cpp index 2b287ba55..028ca7c60 100644 --- a/src/revlanguage/workspace/RbRegister_Type.cpp +++ b/src/revlanguage/workspace/RbRegister_Type.cpp @@ -128,23 +128,6 @@ void RevLanguage::Workspace::initializeTypeGlobalWorkspace(void) try { - AddWorkspaceVectorType::addTypeToWorkspace( *this, new Taxon() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new RateGenerator() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new CladogeneticProbabilityMatrix() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new CladogeneticSpeciationRateMatrix() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new DistanceMatrix() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new MatrixReal() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new MatrixRealPos() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new MatrixRealSymmetric() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new AbstractHomologousDiscreteCharacterData() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new ContinuousCharacterData() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new CharacterHistoryRateModifier() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new TimeTree() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new BranchLengthTree() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new Tree() ); - AddWorkspaceVectorType::addTypeToWorkspace( *this, new Clade() ); -// AddWorkspaceVectorType::addTypeToWorkspace( *this, new Dist_bdp() ); - addTypeWithConstructor( new Clade() ); addTypeWithConstructor( new Taxon() ); diff --git a/src/revlanguage/workspace/RbRegister_VectorType.cpp b/src/revlanguage/workspace/RbRegister_VectorType.cpp new file mode 100644 index 000000000..260fd71e7 --- /dev/null +++ b/src/revlanguage/workspace/RbRegister_VectorType.cpp @@ -0,0 +1,164 @@ +/** + * @file + * This file contains the Workspace function that adds types and functions + * to the global workspace, registering them with the interpreter/compiler + * during the process. + * + * @brief Function registering language objects + * + * Instructions + * + * This is the central registry of Rev objects. It is a large file and needs + * to be properly organized to facilitate maintenance. Follow these simple + * guidelines to ensure that your additions follow the existing structure. + * + * 1. All headers are added in groups corresponding to directories in the + * revlanguage code base. + * 2. All objects (types, distributions, and functions) are registered in + * groups corresponding to directories in the revlanguage code base. + * 3. All entries in each group are listed in alphabetical order. + * + * Some explanation of the directory structure is provided in the comments + * in this file. Consult these comments if you are uncertain about where + * to add your objects in the code. + */ + + +#include +#include +#include +#include + +/* Files including helper classes */ +#include "AddContinuousDistribution.h" +#include "AddDistribution.h" +#include "AddWorkspaceVectorType.h" +#include "AddVectorizedWorkspaceType.h" +#include "RbException.h" +#include "RevAbstractType.h" +#include "RlUserInterface.h" +#include "Workspace.h" + +/// Miscellaneous types /// + +/* Base types (in folder "datatypes") */ +#include "RevObject.h" +#include "AbstractModelObject.h" + +/* Container types (in folder "datatypes/container") */ +#include "RlCorrespondenceAnalysis.h" +#include "RlMatrixReal.h" +#include "RlMatrixRealPos.h" +#include "RlMatrixRealSymmetric.h" +#include "RlRateGeneratorSequence.h" +#include "RlRateMatrix.h" +#include "RlSimplex.h" + +/* Container types (in folder "datatypes/math") */ +#include "ModelVector.h" +#include "WorkspaceVector.h" + +/* Container types (in folder "distributions/phylogenetics") */ +#include "Dist_bdp.h" + +/* Evolution types (in folder "datatypes/phylogenetics") */ +#include "RlDistanceMatrix.h" + +/* Character state types (in folder "datatypes/phylogenetics/character") */ +#include "RlAminoAcidState.h" +#include "RlDnaState.h" +#include "RlRnaState.h" +#include "RlStandardState.h" + +/* Character data types (in folder "datatypes/phylogenetics/datamatrix") */ +#include "RlAbstractCharacterData.h" +#include "RlAbstractHomologousDiscreteCharacterData.h" +#include "RlContinuousCharacterData.h" + +/* Tree types (in folder "datatypes/phylogenetics/trees") */ +#include "RlClade.h" +#include "RlRootedTripletDistribution.h" + + +/* Taxon types (in folder "datatypes/phylogenetics") */ +#include "RlTaxon.h" + +/* Inference types (in folder "analysis") */ +#include "RlBootstrapAnalysis.h" +#include "RlBurninEstimationConvergenceAssessment.h" +#include "RlHillClimber.h" +#include "RlMcmc.h" +#include "RlMcmcmc.h" +#include "RlModel.h" +#include "RlPathSampler.h" +#include "RlPosteriorPredictiveAnalysis.h" +#include "RlPosteriorPredictiveSimulation.h" +#include "RlPowerPosteriorAnalysis.h" +#include "RlSteppingStoneSampler.h" +#include "RlValidationAnalysis.h" +#include "RlAncestralStateTrace.h" + +/// Stopping Rules /// +#include "RlMaxIterationStoppingRule.h" +#include "RlMaxTimeStoppingRule.h" +#include "RlMinEssStoppingRule.h" +#include "RlGelmanRubinStoppingRule.h" +#include "RlGewekeStoppingRule.h" +#include "RlStationarityStoppingRule.h" + + +/// Types /// + +/* These types are needed as template types for the moves */ +#include "RlBranchLengthTree.h" +#include "RlCharacterHistoryRateModifier.h" +#include "RlMonitor.h" +#include "RlMove.h" +#include "RlRateGenerator.h" +#include "RlCladogeneticProbabilityMatrix.h" +#include "RlCladogeneticSpeciationRateMatrix.h" +#include "RlTimeTree.h" + + + +/** Initialize global workspace */ +void RevLanguage::Workspace::initializeVectorTypeGlobalWorkspace(void) +{ + + try + { + + AddWorkspaceVectorType::addTypeToWorkspace( *this, new Taxon() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new RateGenerator() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new CladogeneticProbabilityMatrix() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new CladogeneticSpeciationRateMatrix() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new DistanceMatrix() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new MatrixReal() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new MatrixRealPos() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new MatrixRealSymmetric() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new AbstractHomologousDiscreteCharacterData() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new ContinuousCharacterData() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new CharacterHistoryRateModifier() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new TimeTree() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new BranchLengthTree() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new Tree() ); + AddWorkspaceVectorType::addTypeToWorkspace( *this, new Clade() ); + // AddWorkspaceVectorType::addTypeToWorkspace( *this, new Dist_bdp() ); + + } + catch(RbException& rbException) + { + + RBOUT("Caught an exception while initializing types in the workspace\n"); + std::ostringstream msg; + rbException.print(msg); + msg << std::endl; + RBOUT(msg.str()); + + RBOUT("Please report this bug to the RevBayes Development Core Team"); + + RBOUT("Press any character to exit the program."); + getchar(); + exit(1); + } +} diff --git a/src/revlanguage/workspace/Workspace.cpp b/src/revlanguage/workspace/Workspace.cpp index d8a48f4bb..56d5bd16f 100755 --- a/src/revlanguage/workspace/Workspace.cpp +++ b/src/revlanguage/workspace/Workspace.cpp @@ -259,6 +259,7 @@ void Workspace::initializeGlobalWorkspace( void ) { initializeBasicTypeGlobalWorkspace(); + initializeVectorTypeGlobalWorkspace(); initializeTypeGlobalWorkspace(); initializeMonitorGlobalWorkspace(); initializeMoveGlobalWorkspace(); diff --git a/src/revlanguage/workspace/Workspace.h b/src/revlanguage/workspace/Workspace.h index 10614b411..9319ede00 100755 --- a/src/revlanguage/workspace/Workspace.h +++ b/src/revlanguage/workspace/Workspace.h @@ -100,7 +100,8 @@ namespace RevLanguage { void initializeMonitorGlobalWorkspace(void); //!< Initialize global workspace for monitors void initializeMoveGlobalWorkspace(void); //!< Initialize global workspace for moves void initializeTypeGlobalWorkspace(void); //!< Initialize global workspace for types - + void initializeVectorTypeGlobalWorkspace(void); //!< Initialize global workspace for types + TypeTable typeTable; //!< Type table bool typesInitialized; //!< Are types initialized?