Merge branch 'development' of https://github.com/revbayes/revbayes in…

…to development

src/core/analysis/MonteCarloAnalysis.cpp

@@ -475,12 +475,14 @@ void MonteCarloAnalysis::resetReplicates( void )
     }
 
 
+    size_t replicate_start = size_t(floor( (double(pid-active_PID)   / num_processes ) * replicates ) );
+    RandomNumberGenerator *rng = GLOBAL_RNG;
+    for (size_t j=0; j<(2*replicate_start); ++j) rng->uniform01();
+
+
     // redraw initial states for replicates
     for (size_t i = 0; i < replicates; ++i)
     {
-        RandomNumberGenerator *rng = GLOBAL_RNG;
-        for (size_t j=0; j<10; ++j) rng->uniform01();
-
 
         if ( i > 0 && runs[i] != NULL )
         {

src/core/dag/DeterministicNode.h

@@ -481,7 +481,7 @@ void RevBayesCore::DeterministicNode<valueType>::touchMe( DagNode *toucher, bool
     needs_update = true;
 
     // only if this function did not need an update we delegate the touch affected
-    if ( needed_update == false || was_touched == false || true )
+    if ( needed_update == false || was_touched == false )
     {
         // Dispatch the touch message to downstream nodes
         this->touchAffected( touchAll );

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

@@ -58,24 +58,37 @@ double HostSwitchRateModifier::computeRateMultiplier(std::vector<CharacterEvent*
     // loss event (independent of other hosts)
     if (from_state > to_state)
     {
+//        size_t num_one = sites_with_states[1].size();
+//        size_t num_two = sites_with_states[2].size();
 
         // repertoire must contain at least a single 2 (actual host)
-        // if the current repertoire contains one 2,
-        // and if our event causes the loss of state 2
-        if (sites_with_states[2].size() == 1 && from_state==2)
-        {
-            // cannot enter the null range (conditions on survival)
-            r = 0.0;
-            return r;
-        }
-        else
-        {
-            r = 1.0;
-            return r;
-        }
+//        if (num_two == 1 && from_state==2 && to_state==1)
+//        {
+//            // rate entering 0/1-repertoire is zero
+//            r = 0.0;
+//            return r;
+//        }
+//        else
+//        if (num_two == 0) {
+//            // rate leaving 0/1-repertoire is zero
+//            r = 0.0;
+//            return r;
+//        }
+//        else
+//        {
+//            r = 1.0;
+//            return r;
+//        }
+        return 1.0;
     }
     else
     {
+        // rate of leaving 0/1-repertoire equals zero
+        size_t num_two = sites_with_states[2].size();
+        if (num_two == 0) {
+            return 0.0;
+        }
+
         // gain event
         double scaler_value = scale[ to_state - 1 ];
 

src/core/distributions/phylogenetics/characterhistory/GeneralTreeHistoryCtmc.h

@@ -29,7 +29,7 @@ namespace RevBayesCore {
     class GeneralTreeHistoryCtmc : public TreeHistoryCtmc<charType> {
 
     public:
-        GeneralTreeHistoryCtmc(const TypedDagNode< Tree > *t, size_t nChars, size_t nSites, bool useAmbigChar=false);
+        GeneralTreeHistoryCtmc(const TypedDagNode< Tree > *t, size_t nChars, size_t nSites, bool useAmbigChar=false, bool internal=false);
         GeneralTreeHistoryCtmc(const GeneralTreeHistoryCtmc &n);                                                                         //!< Copy constructor
         virtual                                            ~GeneralTreeHistoryCtmc(void);                                                //!< Virtual destructor
 
@@ -104,7 +104,7 @@ namespace RevBayesCore {
 #include "RbConstants.h"
 
 template<class charType>
-RevBayesCore::GeneralTreeHistoryCtmc<charType>::GeneralTreeHistoryCtmc(const TypedDagNode<Tree> *tau, size_t nChars, size_t nSites, bool useAmbigChar) : TreeHistoryCtmc<charType>( tau, nChars, nSites, useAmbigChar )
+RevBayesCore::GeneralTreeHistoryCtmc<charType>::GeneralTreeHistoryCtmc(const TypedDagNode<Tree> *tau, size_t nChars, size_t nSites, bool useAmbigChar, bool internal) : TreeHistoryCtmc<charType>( tau, nChars, nSites, useAmbigChar, internal )
 {
 
     // initialize with default parameters

src/core/distributions/phylogenetics/characterhistory/TreeHistoryCtmc.h

@@ -29,7 +29,7 @@ namespace RevBayesCore {
 
     public:
         // Note, we need the size of the alignment in the constructor to correctly simulate an initial state
-        TreeHistoryCtmc(const TypedDagNode<Tree> *t, size_t nChars, size_t nSites, bool useAmbigChar=false);
+        TreeHistoryCtmc(const TypedDagNode<Tree> *t, size_t nChars, size_t nSites, bool useAmbigChar=false, bool internal=false);
         TreeHistoryCtmc(const TreeHistoryCtmc &n);                                                                              //!< Copy constructor
         virtual                                                            ~TreeHistoryCtmc(void);                                      //!< Virtual destructor
 
@@ -121,6 +121,7 @@ namespace RevBayesCore {
         bool                                                                tipsInitialized;
         bool                                                                branchHeterogeneousClockRates;
         bool                                                                useDirtyNodes;
+        bool                                                                store_internal_nodes;
 
         charType                                                            template_state;                                 //!< Template state used for ancestral state estimation. This makes sure that the state labels are preserved.
 
@@ -146,7 +147,7 @@ namespace RevBayesCore {
 
 
 template<class charType>
-RevBayesCore::TreeHistoryCtmc<charType>::TreeHistoryCtmc(const TypedDagNode<Tree> *t, size_t nChars, size_t nSites, bool useAmbigChar) : TypedDistribution< AbstractHomologousDiscreteCharacterData >(  new HomologousDiscreteCharacterData<charType>() ),
+RevBayesCore::TreeHistoryCtmc<charType>::TreeHistoryCtmc(const TypedDagNode<Tree> *t, size_t nChars, size_t nSites, bool useAmbigChar, bool internal) : TypedDistribution< AbstractHomologousDiscreteCharacterData >(  new HomologousDiscreteCharacterData<charType>() ),
     num_states( nChars ),
     num_sites( nSites ),
     num_site_rates( 1 ),
@@ -163,7 +164,8 @@ RevBayesCore::TreeHistoryCtmc<charType>::TreeHistoryCtmc(const TypedDagNode<Tree
     treatAmbiguousAsGaps( true ),
     tipsInitialized( false ),
     useDirtyNodes(false),
-    rootBranchLengthMultiplier(1)
+    rootBranchLengthMultiplier(1),
+    store_internal_nodes( internal )
 {
     // initialize with default parameters
     homogeneousClockRate        = new ConstantNode<double>("clockRate", new double(1.0) );
@@ -213,7 +215,8 @@ RevBayesCore::TreeHistoryCtmc<charType>::TreeHistoryCtmc(const TreeHistoryCtmc &
     treatAmbiguousAsGaps( n.treatAmbiguousAsGaps ),
     tipsInitialized( n.tipsInitialized ),
     template_state( n.template_state ),
-    rootBranchLengthMultiplier( n.rootBranchLengthMultiplier )
+    rootBranchLengthMultiplier( n.rootBranchLengthMultiplier ),
+    store_internal_nodes( n.store_internal_nodes )
 {
 
     homogeneousClockRate        = n.homogeneousClockRate;
@@ -829,12 +832,25 @@ void RevBayesCore::TreeHistoryCtmc<charType>::simulate(void)
 
     simulate(nd, bh, taxa);
 
-    // add the taxon data to the character data
+    // add the taxon data for tips to the character data
     for (size_t i = 0; i < tau->getValue().getNumberOfTips(); ++i)
     {
         this->value->addTaxonData( taxa[i] );
         //        this->value->getTaxonData( tau->getValue().getNodes()[i]->getName() );
     }
+    // add the taxon data for internal nodes to character data
+    if (store_internal_nodes) {
+        std::vector<TopologyNode*> nodes = tau->getValue().getNodes();
+        for (size_t i = tau->getValue().getNumberOfTips()+1; i < tau->getValue().getNumberOfNodes(); ++i)
+        {
+            std::stringstream ss;
+            size_t node_index = nodes[i]->getIndex();
+            ss << "Index_" << node_index + 1;
+            taxa[node_index].setTaxon( Taxon(ss.str()) );
+            this->value->addTaxonData( taxa[node_index] );
+            //        this->value->getTaxonData( tau->getValue().getNodes()[i]->getName() );
+        }
+    }
 
     TypedDistribution< AbstractHomologousDiscreteCharacterData >::setValue(this->value);
 }

src/core/math/RbStatisticsHelper.cpp

@@ -497,32 +497,32 @@ double RbStatistics::Helper::rndGamma(double s, RandomNumberGenerator& rng) {
 double RbStatistics::Helper::rndGamma1(double s, RandomNumberGenerator& rng)
 {
 
-	double			r, x = 0.0, small = 1e-37, w;
-	static double   a, p, uf, ss = 10.0, d;
-	
-	if (s != ss) {
-		a  = 1.0 - s;
-		p  = a / (a + s * exp(-a));
-		uf = p * pow(small / a, s);
-		d  = a * log(a);
-		ss = s;
+    double r, x = 0.0, small = 1e-37, w;
+    static double   a, p, uf, ss = 10.0, d;
+
+    if (s != ss) {
+        a  = 1.0 - s;
+        p  = a / (a + s * exp(-a));
+        uf = p * pow(small / a, s);
+        d  = a * log(a);
+        ss = s;
     }
-	for (;;) {
-		r = rng.uniform01();
-		if (r > p)        
-			x = a - log((1.0 - r) / (1.0 - p)), w = a * log(x) - d;
-		else if (r>uf)  
-			x = a * pow(r / p, 1.0 / s), w = x;
-		else            
-			return (0.0);
-		r = rng.uniform01();
-		if (1.0 - r <= w && r > 0.0)
+    for (;;) {
+        r = rng.uniform01();
+        if (r > p)
+            x = a - log((1.0 - r) / (1.0 - p)), w = a * log(x) - d;
+        else if (r>uf)
+            x = a * pow(r / p, 1.0 / s), w = x;
+        else
+            return (0.0);
+        r = rng.uniform01();
+        if (1.0 - r <= w && r > 0.0)
             if (r*(w + 1.0) >= 1.0 || -log(r) <= w)  
                 continue;
-		break;
+        break;
     }
 
-	return (x);
+    return (x);
 }
 
 /*!
@@ -533,30 +533,31 @@ double RbStatistics::Helper::rndGamma1(double s, RandomNumberGenerator& rng)
  * \return Returns a gamma-distributed random variable. 
  * \throws Does not throw an error.
  */
-double RbStatistics::Helper::rndGamma2(double s, RandomNumberGenerator& rng) {
-
-	double			r, d, f, g, x;
-	static double	b, h, ss = 0.0;
-
-	if (s != ss) {
-		b  = s - 1.0;
-		h  = sqrt(3.0 * s - 0.75);
-		ss = s;
+double RbStatistics::Helper::rndGamma2(double s, RandomNumberGenerator& rng)
+{
+
+    double              r, d, f, g, x;
+    static double       b, h, ss = 0.0;
+
+    if (s != ss) {
+        b  = s - 1.0;
+        h  = sqrt(3.0 * s - 0.75);
+        ss = s;
     }
-	for (;;) {
-		r = rng.uniform01();
-		g = r - r * r;
-		f = (r - 0.5) * h / sqrt(g);
-		x = b + f;
-		if (x <= 0.0) 
-			continue;
-		r = rng.uniform01();
-		d = 64.0 * r * r * g * g * g;
-		if (d * x < x - 2.0 * f * f || log(d) < 2.0 * (b * log(x / b) - f))  
-			break;
+    for (;;) {
+        r = rng.uniform01();
+        g = r - r * r;
+        f = (r - 0.5) * h / sqrt(g);
+        x = b + f;
+        if (x <= 0.0)
+            continue;
+        r = rng.uniform01();
+        d = 64.0 * r * r * g * g * g;
+        if (d * x < x - 2.0 * f * f || log(d) < 2.0 * (b * log(x / b) - f))  
+            break;
     }
 
-	return x;
+    return x;
 }
 
 /*!