Add the 4 new function prototypes and 3 new data members highlighted below to the Model class declaration.
#pragma once
#include <algorithm>
#include <vector>
#include "datatype.hpp"
#include "qmatrix.hpp"
#include "asrv.hpp"
#include "libhmsbeagle/beagle.h"
#include <boost/format.hpp>
#include <boost/math/distributions/gamma.hpp>
#include <Eigen/Dense>
namespace strom {
class Likelihood;
class Model {
friend class Likelihood;
public:
typedef std::vector<ASRV::SharedPtr> asrv_vect_t;
typedef std::vector<QMatrix::SharedPtr> qmat_vect_t;
typedef std::vector<unsigned> subset_sizes_t;
typedef std::vector<DataType> subset_datatype_t;
typedef std::vector<double> subset_relrate_vect_t;
typedef std::vector<QMatrix::SharedPtr> state_freq_params_t;
typedef std::vector<QMatrix::SharedPtr> exchangeability_params_t;
typedef std::vector<QMatrix::SharedPtr> omega_params_t;
typedef std::vector<ASRV::SharedPtr> ratevar_params_t;
typedef std::vector<ASRV::SharedPtr> pinvar_params_t;
typedef boost::shared_ptr<Model> SharedPtr;
Model();
~Model();
void activate();
void inactivate();
std::string describeModel();
void setSubsetDataTypes(const subset_datatype_t & datatype_vect);
void setSubsetRateVar(ASRV::ratevar_ptr_t ratevar, unsigned subset, bool fixed);
void setSubsetPinvar(ASRV::pinvar_ptr_t pinvar, unsigned subset, bool fixed);
void setSubsetExchangeabilities(QMatrix::freq_xchg_ptr_t exchangeabilities, unsigned subset, bool fixed);
void setSubsetStateFreqs(QMatrix::freq_xchg_ptr_t state_frequencies, unsigned subset, bool fixed);
void setSubsetOmega(QMatrix::omega_ptr_t omega, unsigned subset, bool fixed);
void setSubsetRelRates(subset_relrate_vect_t & relrates, bool fixed);
subset_relrate_vect_t & getSubsetRelRates();
bool isFixedSubsetRelRates() const;
double calcNormalizingConstantForSubsetRelRates() const;
void setTreeIndex(unsigned i, bool fixed);
unsigned getTreeIndex() const;
bool isFixedTree() const;
<span style="color:#0000ff"><strong>void setTopologyPriorOptions(bool allow_polytomies, bool resclass, double C);</strong></span>
<span style="color:#0000ff"><strong>bool isResolutionClassTopologyPrior() const;</strong></span>
<span style="color:#0000ff"><strong>double getTopologyPriorC() const;</strong></span>
<span style="color:#0000ff"><strong>bool isAllowPolytomies() const;</strong></span>
unsigned getNumSubsets() const;
unsigned getNumSites() const;
unsigned getSubsetNumSites(unsigned subset) const;
const QMatrix & getQMatrix(unsigned subset) const;
const ASRV & getASRV(unsigned subset) const;
void setSubsetIsInvarModel(bool is_invar, unsigned subset);
bool getSubsetIsInvarModel(unsigned subset) const;
void setSubsetSizes(const subset_sizes_t nsites_vect);
subset_sizes_t & getSubsetSizes();
void setSubsetNumPatterns(const subset_sizes_t npatterns_vect);
unsigned getSubsetNumPatterns(unsigned subset) const;
void setSubsetNumCateg(unsigned ncateg, unsigned subset);
unsigned getSubsetNumCateg(unsigned subset) const;
state_freq_params_t & getStateFreqParams();
exchangeability_params_t & getExchangeabilityParams();
omega_params_t & getOmegaParams();
ratevar_params_t & getRateVarParams();
pinvar_params_t & getPinvarParams();
int setBeagleEigenDecomposition(int beagle_instance, unsigned subset, unsigned instance_subset);
int setBeagleStateFrequencies(int beagle_instance, unsigned subset, unsigned instance_subset);
int setBeagleAmongSiteRateVariationRates(int beagle_instance, unsigned subset, unsigned instance_subset);
int setBeagleAmongSiteRateVariationProbs(int beagle_instance, unsigned subset, unsigned instance_subset);
std::string paramNamesAsString(std::string sep) const;
std::string paramValuesAsString(std::string sep) const;
private:
void clear();
unsigned _num_subsets;
unsigned _num_sites;
subset_sizes_t _subset_sizes;
subset_sizes_t _subset_npatterns;
subset_datatype_t _subset_datatypes;
qmat_vect_t _qmatrix;
asrv_vect_t _asrv;
bool _tree_index;
bool _tree_fixed;
<span style="color:#0000ff"><strong>bool _allow_polytomies;</strong></span>
<span style="color:#0000ff"><strong>bool _resolution_class_prior;</strong></span>
<span style="color:#0000ff"><strong>double _topo_prior_C;</strong></span>
bool _subset_relrates_fixed;
subset_relrate_vect_t _subset_relrates;
state_freq_params_t _state_freq_params;
exchangeability_params_t _exchangeability_params;
omega_params_t _omega_params;
ratevar_params_t _ratevar_params;
pinvar_params_t _pinvar_params;
};
Initialize the 2 new data members (_resolution_class_prior and _topo_prior_C) in the clear function.
inline void Model::clear() {
_num_subsets = 0;
_num_sites = 0;
_tree_index = 0;
_tree_fixed = false;
<span style="color:#0000ff"><strong>_allow_polytomies = true;</strong></span>
<span style="color:#0000ff"><strong>_resolution_class_prior = true;</strong></span>
<span style="color:#0000ff"><strong>_topo_prior_C = 1.0;</strong></span>
_subset_relrates_fixed = false;
_subset_relrates.clear();
_subset_sizes.clear();
_subset_npatterns.clear();
_subset_datatypes.clear();
_qmatrix.clear();
_asrv.clear();
}
Add the body of the setTopologyPriorOptions function to model.hpp somehwere before the bracket that closes the strom namespace. This function provides a way to inform the Model object of all 3 parameters of the topology prior.
inline void Model::setTopologyPriorOptions(bool allow_polytomies, bool resclass, double C) {
_allow_polytomies = allow_polytomies;
_resolution_class_prior = resclass;
_topo_prior_C = C;
}
Finally, add bodies of the 3 accessor functions isAllowPolytomies, isResolutionClassTopologyPrior, and getTopologyPriorC.
inline bool Model::isAllowPolytomies() const {
return _allow_polytomies;
}
inline bool Model::isResolutionClassTopologyPrior() const {
return _resolution_class_prior;
}
inline double Model::getTopologyPriorC() const {
return _topo_prior_C;
}
The only tweak we need to make to OutputManager is to add the resolution class of the current tree to the quantities output to the parameter file.
In the class declaration in the file output_manager.hpp, add unsigned m to the parameter list of the outputParameters function.
#pragma once
#include "data.hpp"
#include "tree_manip.hpp"
#include "model.hpp"
#include "xstrom.hpp"
#include <fstream>
namespace strom {
class OutputManager {
public:
typedef std::shared_ptr< OutputManager > SharedPtr;
OutputManager();
~OutputManager();
void setModel(Model::SharedPtr model) {_model = model;}
void setTreeManip(TreeManip::SharedPtr tm) {_tree_manip = tm;}
void openTreeFile(std::string filename, Data::SharedPtr data);
void openParameterFile(std::string filename, Model::SharedPtr model);
void closeTreeFile();
void closeParameterFile();
void outputConsole(std::string s);
void outputTree(unsigned iter, TreeManip::SharedPtr tm);
<span style="color:#0000ff"><strong>void outputParameters(unsigned iter, double lnL, double lnP, double TL, unsigned m, Model::SharedPtr model);</strong></span>
private:
TreeManip::SharedPtr _tree_manip;
Model::SharedPtr _model;
std::ofstream _treefile;
std::ofstream _parameterfile;
std::string _tree_file_name;
std::string _param_file_name;
};
In the function openParameterFile, replace the line highlighted below. The difference is that a column heading "m" has been added for the column showing the resolution class.
inline void OutputManager::openParameterFile(std::string filename, Model::SharedPtr model) {
assert(model);
assert(!_parameterfile.is_open());
_param_file_name = filename;
_parameterfile.open(_param_file_name.c_str());
if (!_parameterfile.is_open())
throw XStrom(boost::str(boost::format("Could not open parameter file \"%s\"") % _param_file_name));
<span style="color:#0000ff"><strong>_parameterfile << boost::str(boost::format("%s\t%s\t%s\t%s\t%s\t%s") % "iter" % "lnL" % "lnPr" % "TL" % "m" % model->paramNamesAsString("\t")) << std::endl;</strong></span>
}
Finally, in the outputParameters function, replace the lines highlighted below. The difference is that a parameter m has been added to the function and the value of m is now output in the new resolution class column.
<span style="color:#0000ff"><strong>inline void OutputManager::outputParameters(unsigned iter, double lnL, double lnP, double TL, unsigned m, Model::SharedPtr model) {</strong></span>
assert(model);
assert(_parameterfile.is_open());
<span style="color:#0000ff"><strong>_parameterfile << boost::str(boost::format("%d\t%.5f\t%.5f\t%.5f\t%d\t%s") % iter % lnL % lnP % TL % m % model->paramValuesAsString("\t")) << std::endl;</strong></span>
}