20.3 Updating the Updater

(Linux version)

< 20.2 | 20.3 | 20.4 >

To finish up, we need to add some functionality to the Updater class. First, add a data member and a member function to the class declaration to inform an Updater object that the heating power should only be applied to the likelihood, not both prior and likelihood.

    class Updater { 
    
        friend class Chain;

        public:
            typedef std::shared_ptr&lt;Updater&gt;        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);
            <span style="color:#0000ff"><strong>void                                    setHeatLikelihoodOnly(bool yes);</strong></span>
            void                                    setTuning(bool on);
            void                                    setTargetAcceptanceRate(double target);
            void                                    setPriorParameters(const std::vector&lt;double&gt; & c);
            void                                    setTopologyPriorOptions(bool resclass, double C);
            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&lt;double&gt;                     _prior_parameters;

            <span style="color:#0000ff"><strong>bool                                    _heat_likelihood_only;</strong></span>
            double                                  _heating_power;
            mutable PolytomyTopoPriorCalculator     _topo_prior_calculator;
            
            static const double                     _log_zero;
    }; 

Next, initialize the data member in the clear function.

    inline void Updater::clear() { 
        _name                   = "updater";
        _tuning                 = true;
        _lambda                 = 0.0001;
        _weight                 = 1.0;
        _prob                   = 0.0;
        _target_acceptance      = 0.3;
        _naccepts               = 0;
        _nattempts              = 0;
        _heating_power          = 1.0;
        <span style="color:#0000ff"><strong>_heat_likelihood_only   = false;</strong></span>
        _prior_parameters.clear();
        reset();
    } 

We use the value stored in the _heat_likelihood_only data member inside the update function:

    inline double Updater::update(double prev_lnL) { 
        double prev_log_prior = calcLogPrior();
        
        // Clear any nodes previously selected so that we can detect those nodes
        // whose partials and/or transition probabilities need to be recalculated
        _tree_manipulator-&gt;deselectAllPartials();
        _tree_manipulator-&gt;deselectAllTMatrices();
        
        // Set model to proposed state and calculate _log_hastings_ratio
        proposeNewState();
        
        // Use alternative partials and transition probability buffer for any selected nodes
        // This allows us to easily revert to the previous values if the move is rejected
        _tree_manipulator-&gt;flipPartialsAndTMatrices();

        // Calculate the log-likelihood and log-prior for the proposed state
        double log_likelihood = calcLogLikelihood();
        double log_prior = calcLogPrior();
        
        // Decide whether to accept or reject the proposed state
        bool accept = true;
        if (log_prior &gt; _log_zero) {
            double log_R = 0.0;
            log_R += _heating_power*(log_likelihood - prev_lnL);
            <span style="color:#0000ff"><strong>//log_R += _heating_power*(log_prior - prev_log_prior);</strong></span>
            <span style="color:#0000ff"><strong>log_R += (_heat_likelihood_only ? 1.0 : _heating_power)*(log_prior - prev_log_prior);</strong></span>
            log_R += _log_hastings_ratio;
            log_R += _log_jacobian;
            
            double logu = _lot-&gt;logUniform();
            if (logu &gt; log_R)
                accept = false;
        }
        else
            accept = false;

        if (accept) {
            _naccepts++;
        }
        else {
            revert();
            _tree_manipulator-&gt;flipPartialsAndTMatrices();
            log_likelihood = prev_lnL;
        }

        tune(accept);
        reset();

        return log_likelihood;
    } 

Finally, add a setting function for the _heat_likelihood_only data member.

    inline void Updater::setHeatLikelihoodOnly(bool yes) {  
        _heat_likelihood_only = yes;
    } 