The following highlighted changes need to be made to the class declaration in updater.hpp:
#pragma once
#include "tree.hpp"
#include "tree_manip.hpp"
#include "lot.hpp"
#include "xstrom.hpp"
#include "likelihood.hpp"
<span style="color:#0000ff"><strong>#include "topo_prior_calculator.hpp"</strong></span>
namespace strom {
class Chain;
class Updater {
friend class Chain;
public:
typedef std::shared_ptr<Updater> SharedPtr;
TreeManip::SharedPtr getTreeManip() const;
Updater();
virtual ~Updater();
void setLikelihood(Likelihood::SharedPtr likelihood);
void setTreeManip(TreeManip::SharedPtr treemanip);
void setLot(Lot::SharedPtr lot);
void setLambda(double lambda);
void setHeatingPower(double p);
void setTuning(bool on);
void setTargetAcceptanceRate(double target);
void setPriorParameters(const std::vector<double> & c);
<span style="color:#0000ff"><strong>void setTopologyPriorOptions(bool resclass, double C);</strong></span>
void setWeight(double w);
void calcProb(double wsum);
double getLambda() const;
double getWeight() const;
double getProb() const;
double getAcceptPct() const;
double getNumUpdates() const;
std::string getUpdaterName() const;
virtual void clear();
virtual double calcLogPrior() = 0;
double calcLogTopologyPrior() const;
double calcLogEdgeLengthPrior() const;
double calcLogLikelihood() const;
virtual double update(double prev_lnL);
static double getLogZero();
protected:
virtual void reset();
virtual void tune(bool accepted);
virtual void revert() = 0;
virtual void proposeNewState() = 0;
Lot::SharedPtr _lot;
Likelihood::SharedPtr _likelihood;
TreeManip::SharedPtr _tree_manipulator;
std::string _name;
double _weight;
double _prob;
double _lambda;
double _log_hastings_ratio;
double _log_jacobian;
double _target_acceptance;
unsigned _naccepts;
unsigned _nattempts;
bool _tuning;
std::vector<double> _prior_parameters;
double _heating_power;
<span style="color:#0000ff"><strong>mutable PolytomyTopoPriorCalculator _topo_prior_calculator;</strong></span>
static const double _log_zero;
};
// member function bodies go here
This function must be modified to take account of polytomous trees. It uses the TopoPriorCalculator class to compute the prior for a tree topology with a specified number of internal nodes.
inline double Updater::calcLogTopologyPrior() const {
<span style="color:#0000ff"><strong>Tree::SharedPtr tree = _tree_manipulator->getTree();</strong></span>
<span style="color:#0000ff"><strong>assert(tree);</strong></span>
<span style="color:#0000ff"><strong>if (tree->isRooted())</strong></span>
<span style="color:#0000ff"><strong>_topo_prior_calculator.chooseRooted();</strong></span>
<span style="color:#0000ff"><strong>else</strong></span>
<span style="color:#0000ff"><strong>_topo_prior_calculator.chooseUnrooted();</strong></span>
<span style="color:#0000ff"><strong>_topo_prior_calculator.setNTax(tree->numLeaves());</strong></span>
<span style="color:#0000ff"><strong>unsigned m = tree->numInternals();</strong></span>
<span style="color:#0000ff"><strong></strong></span>
<span style="color:#0000ff"><strong>double log_topology_prior = _topo_prior_calculator.getLogNormalizedTopologyPrior(m);</strong></span>
<span style="color:#0000ff"><strong></strong></span>
<span style="color:#0000ff"><strong>return log_topology_prior;</strong></span>
}
This function calls the appropriate functions in the PolytomyTopoPriorCalculator object to set the C parameter for the topology prior as well as the type of polytomy prior (resolution class or standard polytomy).
inline void Updater::setTopologyPriorOptions(bool resclass, double C) {
_topo_prior_calculator.setC(C);
if (resclass)
_topo_prior_calculator.chooseResolutionClassPrior();
else
_topo_prior_calculator.choosePolytomyPrior();
}
The calcLogEdgeLengthPrior function must be modified because, if polytomies may be present, we must now count the number of edges in the tree and can can no longer use a formula that only works for binary trees.
inline double Updater::calcLogEdgeLengthPrior() const {
double log_prior = 0.0;
Tree::SharedPtr tree = _tree_manipulator->getTree();
assert(tree);
double TL = _tree_manipulator->calcTreeLength();
<span style="color:#0000ff"><strong>//double n = tree->numLeaves();</strong></span>
<span style="color:#0000ff"><strong>double num_edges = _tree_manipulator->countEdges();</strong></span>
double a = _prior_parameters[0]; // shape of Gamma prior on TL
double b = _prior_parameters[1]; // scale of Gamma prior on TL
double c = _prior_parameters[2]; // parameter of Dirichlet prior on edge length proportions
// Calculate Gamma prior on tree length (TL)
double log_gamma_prior_on_TL = (a - 1.0)*log(TL) - TL/b - a*log(b) - std::lgamma(a);
// Calculate Dirichlet prior on edge length proportions
//
// Note that, for n edges, the Dirichlet prior density is
//
// p1^{c-1} p2^{c-1} ... pn^{c-1}
// ------------------------------
// Gamma(c)^n / Gamma(n*c)
//
// where n = num_edges, pk = edge length k / TL and Gamma is the Gamma function.
// If c == 1, then both numerator and denominator equal 1, so it is pointless
// do loop over edge lengths.
double log_edge_length_proportions_prior = std::lgamma(num_edges*c);
if (c != 1.0) {
for (auto nd : tree->_preorder) {
double edge_length_proportion = nd->_edge_length/TL;
log_edge_length_proportions_prior += (c - 1.0)*log(edge_length_proportion);
}
log_edge_length_proportions_prior -= std::lgamma(c)*num_edges;
}
log_prior = log_gamma_prior_on_TL + log_edge_length_proportions_prior;
return log_prior;
}