19.8 Modifying Model and OutputManager

(Linux version)

< 19.7 | 19.8 | 19.9 >

Update the Model class

Add the 4 new function prototypes and 3 new data members highlighted below to the Model class declaration.

#pragma once    

#include &lt;algorithm&gt;
#include &lt;vector&gt;
#include "datatype.hpp"
#include "qmatrix.hpp"
#include "asrv.hpp"
#include "libhmsbeagle/beagle.h"
#include &lt;boost/format.hpp&gt;
#include &lt;boost/math/distributions/gamma.hpp&gt;
#include &lt;Eigen/Dense&gt;

namespace strom {
    
    class Likelihood;

    class Model { 
        
        friend class Likelihood;

        public:
            typedef std::vector&lt;ASRV::SharedPtr&gt;      asrv_vect_t;
            typedef std::vector&lt;QMatrix::SharedPtr&gt;   qmat_vect_t;
            typedef std::vector&lt;unsigned&gt;             subset_sizes_t;
            typedef std::vector&lt;DataType&gt;             subset_datatype_t;
            typedef std::vector&lt;double&gt;               subset_relrate_vect_t;
            typedef std::vector&lt;QMatrix::SharedPtr&gt;   state_freq_params_t;
            typedef std::vector&lt;QMatrix::SharedPtr&gt;   exchangeability_params_t;
            typedef std::vector&lt;QMatrix::SharedPtr&gt;   omega_params_t;
            typedef std::vector&lt;ASRV::SharedPtr&gt;      ratevar_params_t;
            typedef std::vector&lt;ASRV::SharedPtr&gt;      pinvar_params_t;
            typedef boost::shared_ptr&lt;Model&gt;          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;
        };
    

Modify the clear function

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 setTopologyPriorOptions member function

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

Add the accessors isAllowPolytomies, isResolutionClassTopologyPrior, and getTopologyPriorC

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

Update the OutputManager class

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 &lt;fstream&gt;

namespace strom {

    class OutputManager {
        public:
            typedef std::shared_ptr&lt; OutputManager &gt;            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 &lt;&lt; boost::str(boost::format("%s\t%s\t%s\t%s\t%s\t%s") % "iter" % "lnL" % "lnPr" % "TL" % "m" % model-&gt;paramNamesAsString("\t")) &lt;&lt; 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 &lt;&lt; boost::str(boost::format("%d\t%.5f\t%.5f\t%.5f\t%d\t%s") % iter % lnL % lnP % TL % m % model-&gt;paramValuesAsString("\t")) &lt;&lt; std::endl;</strong></span>
    }

< 19.7 | 19.8 | 19.9 >