Merge commit 'c2f106d768ce47f72aa60c93439da3f57f574f5f'

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<CharacterEvent*> currState, CharacterEventDiscrete* newState, std::vector<std::set<size_t> > sites_with_states, double age=0.0)
-double HostSwitchRateModifier::computeRateMultiplier(std::vector<CharacterEvent *> currState, CharacterEventDiscrete* newState, double age)
+
+double HostSwitchRateModifier::computeRateMultiplier(std::vector<CharacterEvent*> currState, CharacterEventDiscrete* newState, std::vector<std::set<size_t> > 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<CharacterEventDiscrete*>(currState[index])->getState();
+    size_t from_state = static_cast<CharacterEventDiscrete*>(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<CharacterEventDiscrete*>(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<CharacterEvent*> currState, CharacterEventDiscrete* newState, double age)
+{
+    std::vector<std::set<size_t> > sites_with_states(num_states);
+    for (size_t i = 0; i < currState.size(); i++)
+    {
+        sites_with_states[ static_cast<CharacterEventDiscrete*>(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<CharacterEvent*> currState, CharacterEventDiscrete* newState, std::vector<size_t> counts, double age)
+{
+    std::vector<std::set<size_t> > sites_with_states(num_states);
+    for (size_t i = 0; i < currState.size(); i++)
     {
-        size_t s = static_cast<CharacterEventDiscrete*>(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<CharacterEventDiscrete*>(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<double>& s)

src/core/datatypes/phylogenetics/ratemap/HostSwitchRateModifier.h

@@ -33,6 +33,8 @@ namespace RevBayesCore
 
         HostSwitchRateModifier&   assign(const Assignable &m);
         double                              computeRateMultiplier(std::vector<CharacterEvent*> currState, CharacterEventDiscrete* newState, double age=0.0);
+        double                              computeRateMultiplier(std::vector<CharacterEvent*> currState, CharacterEventDiscrete* newState, std::vector<size_t> counts, double age=0.0);
+        double                              computeRateMultiplier(std::vector<CharacterEvent*> currState, CharacterEventDiscrete* newState, std::vector<std::set<size_t> > 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);
 

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];
 

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;
 

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

src/core/functions/phylogenetics/ratemap/HostSwitchRateModifierFunction.cpp

@@ -17,7 +17,7 @@ using namespace RevBayesCore;
 
 HostSwitchRateModifierFunction::HostSwitchRateModifierFunction(const TypedDagNode<Tree>* t,
                                                                const TypedDagNode<RbVector<double> >* s) :
-TypedFunction<CharacterHistoryRateModifier>( new HostSwitchRateModifier(0, t->getValue().getNumberOfTips() ) ),
+TypedFunction<CharacterHistoryRateModifier>( new HostSwitchRateModifier(3, t->getValue().getNumberOfTips() ) ),
 tau(t),
 scale(s)
 {

src/revlanguage/workspace/RbRegister_Type.cpp

@@ -128,23 +128,6 @@ void RevLanguage::Workspace::initializeTypeGlobalWorkspace(void)
     try
     {
 
-        AddWorkspaceVectorType<Taxon,4>::addTypeToWorkspace( *this, new Taxon() );
-        AddWorkspaceVectorType<RateGenerator,3>::addTypeToWorkspace( *this, new RateGenerator() );
-        AddWorkspaceVectorType<CladogeneticProbabilityMatrix,3>::addTypeToWorkspace( *this, new CladogeneticProbabilityMatrix() );
-        AddWorkspaceVectorType<CladogeneticSpeciationRateMatrix,3>::addTypeToWorkspace( *this, new CladogeneticSpeciationRateMatrix() );
-        AddWorkspaceVectorType<DistanceMatrix,3>::addTypeToWorkspace( *this, new DistanceMatrix() );
-        AddWorkspaceVectorType<MatrixReal,3>::addTypeToWorkspace( *this, new MatrixReal() );
-        AddWorkspaceVectorType<MatrixRealPos,3>::addTypeToWorkspace( *this, new MatrixRealPos() );
-        AddWorkspaceVectorType<MatrixRealSymmetric,3>::addTypeToWorkspace( *this, new MatrixRealSymmetric() );
-        AddWorkspaceVectorType<AbstractHomologousDiscreteCharacterData,3>::addTypeToWorkspace( *this, new AbstractHomologousDiscreteCharacterData() );
-        AddWorkspaceVectorType<ContinuousCharacterData,3>::addTypeToWorkspace( *this, new ContinuousCharacterData() );
-        AddWorkspaceVectorType<CharacterHistoryRateModifier,3>::addTypeToWorkspace( *this, new CharacterHistoryRateModifier() );
-        AddWorkspaceVectorType<TimeTree,3>::addTypeToWorkspace( *this, new TimeTree() );
-		AddWorkspaceVectorType<BranchLengthTree,3>::addTypeToWorkspace( *this, new BranchLengthTree() );
-        AddWorkspaceVectorType<Tree,3>::addTypeToWorkspace( *this, new Tree() );
-        AddWorkspaceVectorType<Clade,3>::addTypeToWorkspace( *this, new Clade() );
-//        AddWorkspaceVectorType<Dist_bdp,3>::addTypeToWorkspace( *this, new Dist_bdp() );
-
         addTypeWithConstructor( new Clade() );
         addTypeWithConstructor( new Taxon() );
 

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 <sstream>
+#include <vector>
+#include <set>
+#include <cstdlib>
+
+/* 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<Taxon,4>::addTypeToWorkspace( *this, new Taxon() );
+        AddWorkspaceVectorType<RateGenerator,3>::addTypeToWorkspace( *this, new RateGenerator() );
+        AddWorkspaceVectorType<CladogeneticProbabilityMatrix,3>::addTypeToWorkspace( *this, new CladogeneticProbabilityMatrix() );
+        AddWorkspaceVectorType<CladogeneticSpeciationRateMatrix,3>::addTypeToWorkspace( *this, new CladogeneticSpeciationRateMatrix() );
+        AddWorkspaceVectorType<DistanceMatrix,3>::addTypeToWorkspace( *this, new DistanceMatrix() );
+        AddWorkspaceVectorType<MatrixReal,3>::addTypeToWorkspace( *this, new MatrixReal() );
+        AddWorkspaceVectorType<MatrixRealPos,3>::addTypeToWorkspace( *this, new MatrixRealPos() );
+        AddWorkspaceVectorType<MatrixRealSymmetric,3>::addTypeToWorkspace( *this, new MatrixRealSymmetric() );
+        AddWorkspaceVectorType<AbstractHomologousDiscreteCharacterData,3>::addTypeToWorkspace( *this, new AbstractHomologousDiscreteCharacterData() );
+        AddWorkspaceVectorType<ContinuousCharacterData,3>::addTypeToWorkspace( *this, new ContinuousCharacterData() );
+        AddWorkspaceVectorType<CharacterHistoryRateModifier,3>::addTypeToWorkspace( *this, new CharacterHistoryRateModifier() );
+        AddWorkspaceVectorType<TimeTree,3>::addTypeToWorkspace( *this, new TimeTree() );
+        AddWorkspaceVectorType<BranchLengthTree,3>::addTypeToWorkspace( *this, new BranchLengthTree() );
+        AddWorkspaceVectorType<Tree,3>::addTypeToWorkspace( *this, new Tree() );
+        AddWorkspaceVectorType<Clade,3>::addTypeToWorkspace( *this, new Clade() );
+        //        AddWorkspaceVectorType<Dist_bdp,3>::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);
+    }
+}