15.3 Modify Strom

(Linux version)

< 15.2 | 15.3 | 15.4 >

Changes to Strom

All that is needed now is to create an OutputManager object and call its functions at the appropriate places. The parts that need to be added (or modified) in the Strom class declaration in the file strom.hpp are marked with bold, blue text below.

#pragma once    

#include &lt;iostream&gt;
#include "data.hpp"
#include "likelihood.hpp"
#include "tree_summary.hpp"
#include "partition.hpp"
#include "lot.hpp"
#include "chain.hpp"
<span style="color:#0000ff"><strong>#include "output_manager.hpp"</strong></span>
#include &lt;boost/program_options.hpp&gt;
#include &lt;boost/filesystem.hpp&gt;
#include &lt;boost/algorithm/string/split.hpp&gt;
#include &lt;boost/algorithm/string/classification.hpp&gt;

namespace strom {

    class Strom {
        public:
                                                    Strom();
                                                    ~Strom();

            void                                    clear();
            void                                    processCommandLineOptions(int argc, const char * argv[]);
            void                                    run();
        
        private:
            bool                                    processAssignmentString(const std::string & which, const std::string & definition);
            void                                    handleAssignmentStrings(const boost::program_options::variables_map & vm, std::string label, const std::vector&lt;std::string&gt; & definitions, std::string default_definition);
            bool                                    splitAssignmentString(const std::string & definition, std::vector&lt;std::string&gt; & vector_of_subset_names, std::vector&lt;double&gt;  & vector_of_values);
            <span style="color:#0000ff"><strong>void                                    sample(unsigned iter, Chain & chain);</strong></span>

            double                                  _expected_log_likelihood;

            std::string                             _data_file_name;
            std::string                             _tree_file_name;
            Partition::SharedPtr                    _partition;

            Data::SharedPtr                         _data;
            Model::SharedPtr                        _model;
            Likelihood::SharedPtr                   _likelihood;
            TreeSummary::SharedPtr                  _tree_summary;
            Lot::SharedPtr                          _lot;

            unsigned                                _random_seed;
            unsigned                                _num_iter;
            unsigned                                _print_freq;
            unsigned                                _sample_freq;

            bool                                    _use_gpu;
            bool                                    _ambig_missing;
            bool                                    _use_underflow_scaling;

            static std::string                      _program_name;
            static unsigned                         _major_version;
            static unsigned                         _minor_version;
            
            <span style="color:#0000ff"><strong>OutputManager::SharedPtr                _output_manager;</strong></span>

    };  

Add the indicated line to the clear function body to initialize the _output_manager data member.

    inline void Strom::clear() {    
        _data_file_name             = "";
        _tree_file_name             = "";
        _tree_summary               = nullptr;
        _partition.reset(new Partition());
        _use_gpu                    = true;
        _ambig_missing              = true;
        _model.reset(new Model());
        _expected_log_likelihood    = 0.0;
        _data                       = nullptr;
        _likelihood                 = nullptr;
        _use_underflow_scaling      = false;
        _lot                        = nullptr;
        _random_seed                = 1;
        _num_iter                   = 1000;
        _print_freq                 = 1;
        _sample_freq                = 1;
        <span style="color:#0000ff"><strong>_output_manager             = nullptr;</strong></span>
    }   

You’ll note that I’ve added a private member function named sample. This function is not essential, but serves to make the run function more tidy by gathering all code lines related to saving a tree and parameter sample together so that they can be called by adding a single line to the run function whereever sampling needs to occur. Here is the body of the sample function.

    inline void Strom::sample(unsigned iteration, Chain & chain) {  
        if (chain.getHeatingPower() &lt; 1.0)
            return;
        
        bool time_to_sample = (bool)(iteration % _sample_freq == 0);
        bool time_to_report = (bool)(iteration % _print_freq == 0);
        if (time_to_sample || time_to_report) {
            double logLike = chain.getLogLikelihood();
            double logPrior = chain.calcLogJointPrior();
            double TL = chain.getTreeManip()-&gt;calcTreeLength();
            if (time_to_report) {
                if (logPrior == Updater::getLogZero())
                    _output_manager-&gt;outputConsole(boost::str(boost::format("%12d %12.5f %12s %12.5f") % iteration % logLike % "-infinity" % TL));
                else
                    _output_manager-&gt;outputConsole(boost::str(boost::format("%12d %12.5f %12.5f %12.5f") % iteration % logLike % logPrior % TL));
            }
            if (time_to_sample) {
                _output_manager-&gt;outputTree(iteration, chain.getTreeManip());
                _output_manager-&gt;outputParameters(iteration, logLike, logPrior, TL, chain.getModel());
            }
        }
    }   

Changes to the run function involve creating an OutputManager object and opening/closing tree and sample files. You will also note that some of the blue, bold text reflects sections of code that have been commented out because their functionality has been replaced by the new OutputManager object. The two calls to the sample function cause output to be sent to the console and/or sample and tree files at the very start (iteration 0) and at each samplefreq (or printfreq) iteration thereafter.

    inline void Strom::run() {  
        std::cout &lt;&lt; "Starting..." &lt;&lt; std::endl;
        std::cout &lt;&lt; "Current working directory: " &lt;&lt; boost::filesystem::current_path() &lt;&lt; std::endl;

        try {
            std::cout &lt;&lt; "\n*** Reading and storing the data in the file " &lt;&lt; _data_file_name &lt;&lt; std::endl;
            _data = Data::SharedPtr(new Data());
            _data-&gt;setPartition(_partition);
            _data-&gt;getDataFromFile(_data_file_name);
            
            _model-&gt;setSubsetNumPatterns(_data-&gt;calcNumPatternsVect());
            _model-&gt;setSubsetSizes(_partition-&gt;calcSubsetSizes());
            _model-&gt;activate();

            // Report information about data partition subsets
            unsigned nsubsets = _data-&gt;getNumSubsets();
            std::cout &lt;&lt; "\nNumber of taxa: " &lt;&lt; _data-&gt;getNumTaxa() &lt;&lt; std::endl;
            std::cout &lt;&lt; "Number of partition subsets: " &lt;&lt; nsubsets &lt;&lt; std::endl;
            for (unsigned subset = 0; subset &lt; nsubsets; subset++) {
                DataType dt = _partition-&gt;getDataTypeForSubset(subset);
                std::cout &lt;&lt; "  Subset " &lt;&lt; (subset+1) &lt;&lt; " (" &lt;&lt; _data-&gt;getSubsetName(subset) &lt;&lt; ")" &lt;&lt; std::endl;
                std::cout &lt;&lt; "    data type: " &lt;&lt; dt.getDataTypeAsString() &lt;&lt; std::endl;
                std::cout &lt;&lt; "    sites:     " &lt;&lt; _data-&gt;calcSeqLenInSubset(subset) &lt;&lt; std::endl;
                std::cout &lt;&lt; "    patterns:  " &lt;&lt; _data-&gt;getNumPatternsInSubset(subset) &lt;&lt; std::endl;
                std::cout &lt;&lt; "    ambiguity: " &lt;&lt; (_ambig_missing || dt.isCodon() ? "treated as missing data (faster)" : "handled appropriately (slower)") &lt;&lt; std::endl;
                }

            std::cout &lt;&lt; "\n*** Resources available to BeagleLib " &lt;&lt; _likelihood-&gt;beagleLibVersion() &lt;&lt; ":\n";
            std::cout &lt;&lt; _likelihood-&gt;availableResources() &lt;&lt; std::endl;

            std::cout &lt;&lt; "\n*** Creating the likelihood calculator" &lt;&lt; std::endl;
            _likelihood = Likelihood::SharedPtr(new Likelihood());
            _likelihood-&gt;setPreferGPU(_use_gpu);
            _likelihood-&gt;setAmbiguityEqualsMissing(_ambig_missing);
            _likelihood-&gt;setData(_data);
            _likelihood-&gt;useUnderflowScaling(_use_underflow_scaling);

            std::cout &lt;&lt; "\n*** Model description" &lt;&lt; std::endl;
            std::cout &lt;&lt; _model-&gt;describeModel() &lt;&lt; std::endl;
            _likelihood-&gt;setModel(_model);

            _likelihood-&gt;initBeagleLib();

            std::cout &lt;&lt; "\n*** Reading and storing the first tree in the file " &lt;&lt; _tree_file_name &lt;&lt; std::endl;
            _tree_summary = TreeSummary::SharedPtr(new TreeSummary());
            _tree_summary-&gt;readTreefile(_tree_file_name, 0);
            Tree::SharedPtr tree = _tree_summary-&gt;getTree(0);
            std::string newick = _tree_summary-&gt;getNewick(0);  

            if (tree-&gt;numLeaves() != _data-&gt;getNumTaxa())
                throw XStrom(boost::format("Number of taxa in tree (%d) does not equal the number of taxa in the data matrix (%d)") % tree-&gt;numLeaves() % _data-&gt;getNumTaxa());

            std::cout &lt;&lt; "\n*** Calculating the likelihood of the tree" &lt;&lt; std::endl;
            TreeManip tm(tree);
            tm.selectAllPartials();
            tm.selectAllTMatrices();
            double lnL = _likelihood-&gt;calcLogLikelihood(tree);
            std::cout &lt;&lt; boost::str(boost::format("log likelihood = %.5f") % lnL) &lt;&lt; std::endl;
            
            if (_expected_log_likelihood != 0.0) 
                std::cout &lt;&lt; boost::str(boost::format("      (expecting %.3f)") % _expected_log_likelihood) &lt;&lt; std::endl;
            
            // Create a Lot object that generates (pseudo)random numbers   
            _lot = Lot::SharedPtr(new Lot);
            _lot-&gt;setSeed(_random_seed);

            <span style="color:#0000ff"><strong>// Create an output manager and open output files</strong></span>
            <span style="color:#0000ff"><strong>_output_manager.reset(new OutputManager);</strong></span>

            // Create a Chain object and take _num_iter steps
            Chain chain;
            unsigned num_free_parameters = chain.createUpdaters(_model, _lot, _likelihood);
            if (num_free_parameters &gt; 0) {
                <span style="color:#0000ff"><strong>_output_manager-&gt;outputConsole(boost::str(boost::format("\n%12s %12s %12s %12s") % "iteration" % "logLike" % "logPrior" % "TL"));</strong></span>
                <span style="color:#0000ff"><strong>_output_manager-&gt;openTreeFile("trees.tre", _data);</strong></span>
                <span style="color:#0000ff"><strong>_output_manager-&gt;openParameterFile("params.txt", _model);</strong></span>
                chain.setTreeFromNewick(newick);
                chain.start();
                
                <span style="color:#0000ff"><strong>// // Output column headers and first line of output showing starting state (iteration 0)</strong></span>
                <span style="color:#0000ff"><strong>//std::cout &lt;&lt; boost::str(boost::format("\n%12s %12s %12s %12s %12s %12s\n")</strong></span>
                <span style="color:#0000ff"><strong>//    % "iteration"</strong></span>
                <span style="color:#0000ff"><strong>//    % "lnLike"</strong></span>
                <span style="color:#0000ff"><strong>//    % "lnPrior"</strong></span>
                <span style="color:#0000ff"><strong>//    % "ratevar"</strong></span>
                <span style="color:#0000ff"><strong>//    % "accept"</strong></span>
                <span style="color:#0000ff"><strong>//    % "samples");</strong></span>
                <span style="color:#0000ff"><strong>//GammaRateVarUpdater::SharedPtr ratevar_updater = std::dynamic_pointer_cast&lt;GammaRateVarUpdater&gt; (chain.findUpdaterByName("Gamma Rate Variance"));</strong></span>
                <span style="color:#0000ff"><strong>//std::cout &lt;&lt; boost::str(boost::format("%12d %12.5f %12.5f %12.5f %12.1f %12d\n")</strong></span>
                <span style="color:#0000ff"><strong>//    % 0</strong></span>
                <span style="color:#0000ff"><strong>//    % chain.getLogLikelihood()</strong></span>
                <span style="color:#0000ff"><strong>//    % chain.calcLogJointPrior()</strong></span>
                <span style="color:#0000ff"><strong>//    % ratevar_updater-&gt;getCurrentPoint()</strong></span>
                <span style="color:#0000ff"><strong>//    % 0</strong></span>
                <span style="color:#0000ff"><strong>//    % 0);</strong></span>
                
                <span style="color:#0000ff"><strong>sample(0, chain);</strong></span>
                for (unsigned iteration = 1; iteration &lt;= _num_iter; ++iteration) {
                    chain.nextStep(iteration);
                    <span style="color:#0000ff"><strong>sample(iteration, chain);</strong></span>
                    
                    <span style="color:#0000ff"><strong>//if (iteration % _sample_freq == 0) {</strong></span>
                    <span style="color:#0000ff"><strong>//    GammaRateVarUpdater::SharedPtr ratevar_updater = std::dynamic_pointer_cast&lt;GammaRateVarUpdater&gt; (chain.findUpdaterByName("Gamma Rate Variance"));</strong></span>
                    <span style="color:#0000ff"><strong>//    double log_prior = chain.calcLogJointPrior();</strong></span>
                    <span style="color:#0000ff"><strong>//    if (log_prior == Updater::getLogZero())</strong></span>
                    <span style="color:#0000ff"><strong>//        std::cout &lt;&lt; boost::str(boost::format("%12d %12.5f %12s %12.5f %12.1f %12d\n")</strong></span>
                    <span style="color:#0000ff"><strong>//            % iteration</strong></span>
                    <span style="color:#0000ff"><strong>//            % chain.getLogLikelihood()</strong></span>
                    <span style="color:#0000ff"><strong>//            % "-infinity"</strong></span>
                    <span style="color:#0000ff"><strong>//            % ratevar_updater-&gt;getCurrentPoint()</strong></span>
                    <span style="color:#0000ff"><strong>//            % ratevar_updater-&gt;getAcceptPct()</strong></span>
                    <span style="color:#0000ff"><strong>//            % ratevar_updater-&gt;getNumUpdates());</strong></span>
                    <span style="color:#0000ff"><strong>//    else</strong></span>
                    <span style="color:#0000ff"><strong>//        std::cout &lt;&lt; boost::str(boost::format("%12d %12.5f %12.5f %12.5f %12.1f %12d\n")</strong></span>
                    <span style="color:#0000ff"><strong>//            % iteration</strong></span>
                    <span style="color:#0000ff"><strong>//            % chain.getLogLikelihood()</strong></span>
                    <span style="color:#0000ff"><strong>//            % log_prior</strong></span>
                    <span style="color:#0000ff"><strong>//            % ratevar_updater-&gt;getCurrentPoint()</strong></span>
                    <span style="color:#0000ff"><strong>//            % ratevar_updater-&gt;getAcceptPct()</strong></span>
                    <span style="color:#0000ff"><strong>//            % ratevar_updater-&gt;getNumUpdates());</strong></span>
                    <span style="color:#0000ff"><strong>//}</strong></span>
                }
                chain.stop(); 
                
                <span style="color:#0000ff"><strong>// Close output files</strong></span>
                <span style="color:#0000ff"><strong>_output_manager-&gt;closeTreeFile();</strong></span>
                <span style="color:#0000ff"><strong>_output_manager-&gt;closeParameterFile();</strong></span>
                
                <span style="color:#0000ff"><strong>// Summarize updater efficiency</strong></span>
                <span style="color:#0000ff"><strong>_output_manager-&gt;outputConsole("\nChain summary:");</strong></span>
                <span style="color:#0000ff"><strong>std::vector&lt;std::string&gt; names = chain.getUpdaterNames();</strong></span>
                <span style="color:#0000ff"><strong>std::vector&lt;double&gt; lambdas    = chain.getLambdas();</strong></span>
                <span style="color:#0000ff"><strong>std::vector&lt;double&gt; accepts    = chain.getAcceptPercentages();</strong></span>
                <span style="color:#0000ff"><strong>std::vector&lt;unsigned&gt; nupdates = chain.getNumUpdates();</strong></span>
                <span style="color:#0000ff"><strong>unsigned n = (unsigned)names.size();</strong></span>
                <span style="color:#0000ff"><strong>_output_manager-&gt;outputConsole(boost::str(boost::format("%30s %15s %15s %15s") % "Updater" % "Tuning Param." % "Accept %" % "No. Updates"));</strong></span>
                <span style="color:#0000ff"><strong>for (unsigned i = 0; i &lt; n; ++i) {</strong></span>
                    <span style="color:#0000ff"><strong>_output_manager-&gt;outputConsole(boost::str(boost::format("%30s %15.3f %15.1f %15d") % names[i] % lambdas[i] % accepts[i] % nupdates[i]));</strong></span>
                <span style="color:#0000ff"><strong>}</strong></span>
            }
            else {
                <span style="color:#0000ff"><strong>_output_manager-&gt;outputConsole("MCMC skipped because there are no free parameters in the model");</strong></span>
            }   
        }
        catch (XStrom & x) {
            std::cerr &lt;&lt; "Strom encountered a problem:\n  " &lt;&lt; x.what() &lt;&lt; std::endl;
        }

        std::cout &lt;&lt; "\nFinished!" &lt;&lt; std::endl;
    }   

< 15.2 | 15.3 | 15.4 >