20.1 Adding settings for steppingstone

(Linux version)

< 20.0 | 20.1 | 20.2 >

We will need to add a function (calcMarginalLikelihood) and a couple of variables (_steppingstone and _ss_alpha) to strom.hpp to hold settings that allow users to specify that they would like to carry out a steppingstone analysis rather than an ordinary MCMC analysis.

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

            void                                    clear();
            void                                    processCommandLineOptions(int argc, const char * argv[]);
            void                                    run();
                    
        private:
            bool                                    processAssignmentString(Model::SharedPtr m, const std::string & which, const std::string & definition);
            void                                    handleAssignmentStrings(Model::SharedPtr m, 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);
            void                                    sample(unsigned iter, Chain & chain);

            void                                    readData();
            void                                    readTrees();
            void                                    showPartitionInfo();
            void                                    showBeagleInfo();
            void                                    showMCMCInfo();
            void                                    calcHeatingPowers();
            <span style="color:#0000ff"><strong>void                                    calcMarginalLikelihood() const;</strong></span>
            void                                    initChains();
            void                                    startTuningChains();
            void                                    stopTuningChains();
            void                                    stepChains(unsigned iteration, bool sampling);
            void                                    swapChains();
            void                                    stopChains();
            void                                    swapSummary() const;
            void                                    showChainTuningInfo() const;

            double                                  _expected_log_likelihood;
            
            double                                  _topo_prior_C;
            bool                                    _allow_polytomies;
            bool                                    _resolution_class_prior;

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

            Data::SharedPtr                         _data;
            std::vector&lt;Likelihood::SharedPtr&gt;      _likelihoods;
            TreeSummary::SharedPtr                  _tree_summary;
            Lot::SharedPtr                          _lot;

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

            unsigned                                _num_burnin_iter; 
            bool                                    _using_stored_data;
            bool                                    _use_gpu;
            <span style="color:#0000ff"><strong>bool                                    _steppingstone;</strong></span>
            <span style="color:#0000ff"><strong>double                                  _ss_alpha;</strong></span>

            bool                                    _ambig_missing;
            unsigned                                _num_chains;
            double                                  _heating_lambda;
            std::vector&lt;Chain&gt;                      _chains;
            std::vector&lt;double&gt;                     _heating_powers;
            std::vector&lt;unsigned&gt;                   _swaps;

            bool                                    _use_underflow_scaling;

            static std::string                      _program_name;
            static unsigned                         _major_version;
            static unsigned                         _minor_version;
            
            OutputManager::SharedPtr                _output_manager;

    };  

Initialize the new variables in the clear function:

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

        _using_stored_data          = true;
        _likelihoods.clear();
        _num_burnin_iter            = 1000;
        _heating_lambda             = 0.5;
        _num_chains                 = 1;
        _chains.resize(0);
        _heating_powers.resize(0);
        _swaps.resize(0);
    }   

Adding program options

Create program options that the user can access in the processCommandLineOptions function:

    inline void Strom::processCommandLineOptions(int argc, const char * argv[]) {   
        std::vector&lt;std::string&gt; partition_statefreq;
        std::vector&lt;std::string&gt; partition_rmatrix;
        std::vector&lt;std::string&gt; partition_omega;
        std::vector&lt;std::string&gt; partition_ratevar;
        std::vector&lt;std::string&gt; partition_pinvar;
        std::vector&lt;std::string&gt; partition_ncateg;
        std::vector&lt;std::string&gt; partition_subsets;
        std::vector&lt;std::string&gt; partition_relrates;
        std::vector&lt;std::string&gt; partition_tree;
        boost::program_options::variables_map vm;
        boost::program_options::options_description desc("Allowed options");
        desc.add_options()
            ("help,h", "produce help message")
            ("version,v", "show program version")
            ("seed,z",        boost::program_options::value(&_random_seed)-&gt;default_value(1),   "pseudorandom number seed")
            ("niter,n",       boost::program_options::value(&_num_iter)-&gt;default_value(1000),   "number of MCMC iterations")
            ("printfreq",  boost::program_options::value(&_print_freq)-&gt;default_value(1),   "skip this many iterations before reporting progress")
            ("samplefreq",  boost::program_options::value(&_sample_freq)-&gt;default_value(1),   "skip this many iterations before sampling next")
            ("datafile,d",  boost::program_options::value(&_data_file_name)-&gt;required(), "name of a data file in NEXUS format")
            ("treefile,t",  boost::program_options::value(&_tree_file_name)-&gt;required(), "name of a tree file in NEXUS format")
            ("subset",  boost::program_options::value(&partition_subsets), "a string defining a partition subset, e.g. 'first:1-1234\3' or 'default[codon:standard]:1-3702'")
            ("ncateg,c", boost::program_options::value(&partition_ncateg), "number of categories in the discrete Gamma rate heterogeneity model")
            ("statefreq", boost::program_options::value(&partition_statefreq), "a string defining state frequencies for one or more data subsets, e.g. 'first,second:0.1,0.2,0.3,0.4'")
            ("omega", boost::program_options::value(&partition_omega), "a string defining the nonsynonymous/synonymous rate ratio omega for one or more data subsets, e.g. 'first,second:0.1'")
            ("rmatrix", boost::program_options::value(&partition_rmatrix), "a string defining the rmatrix for one or more data subsets, e.g. 'first,second:1,2,1,1,2,1'")
            ("ratevar", boost::program_options::value(&partition_ratevar), "a string defining the among-site rate variance for one or more data subsets, e.g. 'first,second:2.5'")
            ("pinvar", boost::program_options::value(&partition_pinvar), "a string defining the proportion of invariable sites for one or more data subsets, e.g. 'first,second:0.2'")
            ("relrate", boost::program_options::value(&partition_relrates), "a string defining the (unnormalized) relative rates for all data subsets (e.g. 'default:3,1,6').")
            ("tree", boost::program_options::value(&partition_tree), "the index of the tree in the tree file (first tree has index = 1)")
            ("topopriorC", boost::program_options::value(&_topo_prior_C)-&gt;default_value(1.0), "topology prior C: tree (or resolution class) with m internal nodes has probability C time greater than tree (or resolution class) with m+1 internal nodes.")
            ("allowpolytomies", boost::program_options::value(&_allow_polytomies)-&gt;default_value(true), "yes or no; if yes, then topopriorC and polytomyprior are used, otherwise topopriorC and polytomyprior are ignored")
            ("resclassprior", boost::program_options::value(&_resolution_class_prior)-&gt;default_value(true), "if yes, topologypriorC will apply to resolution classes; if no, topologypriorC will apply to individual tree topologies")
            ("expectedLnL", boost::program_options::value(&_expected_log_likelihood)-&gt;default_value(0.0), "log likelihood expected")
            ("nchains",       boost::program_options::value(&_num_chains)-&gt;default_value(1),                "number of chains")
            ("heatfactor",    boost::program_options::value(&_heating_lambda)-&gt;default_value(0.5),          "determines how hot the heated chains are")
            ("burnin",        boost::program_options::value(&_num_burnin_iter)-&gt;default_value(100),         "number of iterations used to burn in chains")
            ("usedata",       boost::program_options::value(&_using_stored_data)-&gt;default_value(true),      "use the stored data in calculating likelihoods (specify no to explore the prior)")
            ("gpu",           boost::program_options::value(&_use_gpu)-&gt;default_value(true),                "use GPU if available")
            ("ambigmissing",  boost::program_options::value(&_ambig_missing)-&gt;default_value(true),          "treat all ambiguities as missing data")
            ("underflowscaling",  boost::program_options::value(&_use_underflow_scaling)-&gt;default_value(true),          "scale site-likelihoods to prevent underflow (slower but safer)")
            <span style="color:#0000ff"><strong>("steppingstone", boost::program_options::value(&_steppingstone)-&gt;default_value(false),                "use heated chains to compute marginal likelihood with the steppingstone method")</strong></span>
            <span style="color:#0000ff"><strong>("ssalpha", boost::program_options::value(&_ss_alpha)-&gt;default_value(0.25),                "determines how bunched steppingstone chain powers are toward the prior: chain k of K total chains has power (k/K)^{1/ssalpha}")</strong></span>
        ;   

Modifying the sample function

Add a conditional to the sample function to account for the fact that the chains are doing something quite differently under steppingstone than under normal MCMC sampling:

    inline void Strom::sample(unsigned iteration, Chain & chain) {  
        <span style="color:#0000ff"><strong>if (_steppingstone) {</strong></span>
            <span style="color:#0000ff"><strong>bool time_to_sample = (bool)(iteration % _sample_freq == 0);</strong></span>
            <span style="color:#0000ff"><strong>if (time_to_sample && iteration &gt; 0) {</strong></span>
                <span style="color:#0000ff"><strong>chain.storeLogLikelihood();</strong></span>
            <span style="color:#0000ff"><strong>}</strong></span>

<span style="color:#0000ff"><strong></strong></span>
            <span style="color:#0000ff"><strong>assert(_heating_powers.size() &gt; 0);</strong></span>
            <span style="color:#0000ff"><strong>double largest_power = *(_heating_powers.rbegin());</strong></span>
            <span style="color:#0000ff"><strong>if (chain.getHeatingPower() != largest_power)</strong></span>
                <span style="color:#0000ff"><strong>return;</strong></span>
                
<span style="color:#0000ff"><strong></strong></span>
            <span style="color:#0000ff"><strong>bool time_to_report = (bool)(iteration % _print_freq == 0);</strong></span>
            <span style="color:#0000ff"><strong>if (time_to_report) {</strong></span>
                <span style="color:#0000ff"><strong>double logLike = chain.getLogLikelihood();</strong></span>
                <span style="color:#0000ff"><strong>double logPrior = chain.calcLogJointPrior();</strong></span>
                <span style="color:#0000ff"><strong>double TL = chain.getTreeManip()-&gt;calcTreeLength();</strong></span>
                <span style="color:#0000ff"><strong>unsigned m = chain.getTreeManip()-&gt;calcResolutionClass();</strong></span>
                <span style="color:#0000ff"><strong>if (time_to_report) {</strong></span>
                    <span style="color:#0000ff"><strong>if (logPrior == Updater::getLogZero())</strong></span>
                        <span style="color:#0000ff"><strong>_output_manager-&gt;outputConsole(boost::str(boost::format("%12d %12d %12.5f %12s %12.5f") % iteration % m % logLike % "-infinity" % TL));</strong></span>
                    <span style="color:#0000ff"><strong>else</strong></span>
                        <span style="color:#0000ff"><strong>_output_manager-&gt;outputConsole(boost::str(boost::format("%12d %12d %12.5f %12.5f %12.5f") % iteration % m % logLike % logPrior % TL));</strong></span>
                <span style="color:#0000ff"><strong>}</strong></span>
            <span style="color:#0000ff"><strong>}</strong></span>
        <span style="color:#0000ff"><strong>}</strong></span>
        <span style="color:#0000ff"><strong>else {</strong></span>
            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();
                unsigned m = chain.getTreeManip()-&gt;calcResolutionClass();
                if (time_to_report) {
                    if (logPrior == Updater::getLogZero())
                        _output_manager-&gt;outputConsole(boost::str(boost::format("%12d %12d %12.5f %12s %12.5f") % iteration % m % logLike % "-infinity" % TL));
                    else
                        _output_manager-&gt;outputConsole(boost::str(boost::format("%12d %12d %12.5f %12.5f %12.5f") % iteration % m % logLike % logPrior % TL));
                }
                if (time_to_sample) {
                    _output_manager-&gt;outputTree(iteration, chain.getTreeManip());
                    _output_manager-&gt;outputParameters(iteration, logLike, logPrior, TL, m, chain.getModel());
                }
            }
        <span style="color:#0000ff"><strong>}</strong></span>
    }   

You will note that the steppingstone section is almost the same as the non-steppingstone section, but there are some key differences. Ordinarily (i.e. when we are not doing steppingstone), we only produce a progress report for the cold chain. That explains why we return from the function early if the heating power is less than 1. Under steppingstone, it is arbitrary which chain we use when reporting progress, so I’ve chosen here to report progress for the chain with the largest power, which is the final entry in the _heating_powers vector.

You’ll notice that, in the steppingstone section, I’ve also left out code to save parameter values to the sampled parameters and sampled trees files. Samples from these chains have no use except to estimate the marginal likelihood.

Finally, at the beginning of the steppingstone section, a new function storeLogLikelihood is called for each chain if it is time to sample. This function, as you’ll see when you create it, just stores the current log-likelihood in a vector named _ss_loglikes, which is a data member of the Chain class.

Modifying the calcHeatingPowers function

In the calcHeatingPowers function, we need to setup the chain powers differently for steppingstone than for ordinary Metropolis-coupled MCMC. In MCMCMC, it is not wise to let the chain powers get too small. For example, if the chain power were set to zero, the chain would explore the prior and would almost never serve as a useful scout for the cold chain. On the other hand, the chain powers used in steppingstone are concentrated toward the prior end of the spectrum. The reason for this is best explained graphically. Returning to our coin-flipping example, compare the power posterior distributions for evenly-spaced powers (left) to the power posterior distributions for powers concentrated toward the prior end of the spectrum (right). Comparison of even versus bunched steppingstone heating powers As you can see, evenly-spaced powers do not result in homogeneous ratios between adjacent distributions. It is as if the steppingstones placed to help cross a stream involved a really big first jump and then really tiny jumps thereafter. Bunching powers close to the prior results in ratios that are more even and thus more easily and accurately estimated.

The degree to which the heating powers are pushed towards 0.0 is determined by the _ss_alpha setting. The code highlighted below allocates heating powers for steppingstone based on the number of chains specified and the value of _ss_alpha. Note that if _ss_alpha equals 1, then the powers will be evenly spaced between 0 and 1 (but note that power equal to 1 is never used). If _ss_alpha is less than 1, then the nominal (equally spaced) power is raised to the value 1/_ss_alpha. An example is provided in the comments for K=5 chains and _ss_alpha = 0.25.

    inline void Strom::calcHeatingPowers() { 
        <span style="color:#0000ff"><strong>if (_steppingstone) {</strong></span>
            <span style="color:#0000ff"><strong>// Specify chain heating power for steppingstone</strong></span>
            <span style="color:#0000ff"><strong>// For K = 5 chains and alpha = 0.25 (1/alpha = 4):</strong></span>
            <span style="color:#0000ff"><strong>//   k   chain power</strong></span>
            <span style="color:#0000ff"><strong>// ---------------------</strong></span>
            <span style="color:#0000ff"><strong>//   0   (0/5)^4 = 0.0000 &lt;-- prior</strong></span>
            <span style="color:#0000ff"><strong>//   1   (1/5)^4 = 0.0016</strong></span>
            <span style="color:#0000ff"><strong>//   2   (2/5)^4 = 0.0256</strong></span>
            <span style="color:#0000ff"><strong>//   3   (3/5)^4 = 0.1296</strong></span>
            <span style="color:#0000ff"><strong>//   4   (4/5)^4 = 0.4096</strong></span>
            <span style="color:#0000ff"><strong>//   5   (5/5)^4 = 1.0000 &lt;-- posterior not used</strong></span>
            <span style="color:#0000ff"><strong>double inv_alpha = 1.0/_ss_alpha;</strong></span>
            <span style="color:#0000ff"><strong>double k = 0.0;</strong></span>
            <span style="color:#0000ff"><strong>double K = (double)_heating_powers.size();</strong></span>
            <span style="color:#0000ff"><strong>for (auto & h : _heating_powers) {</strong></span>
                <span style="color:#0000ff"><strong>h = pow(k++/K, inv_alpha);</strong></span>
            <span style="color:#0000ff"><strong>}</strong></span>
        <span style="color:#0000ff"><strong>}</strong></span>
        <span style="color:#0000ff"><strong>else {</strong></span>
            // Specify chain heating power (e.g. _heating_lambda = 0.2)
            // chain_index  power
            //      0       1.000 = 1/(1 + 0.2*0)
            //      1       0.833 = 1/(1 + 0.2*1)
            //      2       0.714 = 1/(1 + 0.2*2)
            //      3       0.625 = 1/(1 + 0.2*3)
            unsigned i = 0;
            for (auto & h : _heating_powers) {
                h = 1.0/(1.0 + _heating_lambda*i++);
            }
        <span style="color:#0000ff"><strong>}</strong></span>
    } 

Adding the calcMarginalLikelihood function

Once all heated chains have finished, the calcMarginalLikelihood function can be called to do the final computation.

    inline void Strom::calcMarginalLikelihood() const { 
        if (_steppingstone) {
            // Calculate the log ratio for each steppingstone
            std::vector&lt;std::pair&lt;double, double&gt; &gt; log_ratio;
            for (auto & c : _chains) {
                log_ratio.push_back(std::make_pair(c.getHeatingPower(), c.calcLogSteppingstoneRatio()));
            }
            
            // Sort log_ratio vector from lowest to highest power
            std::sort(log_ratio.begin(), log_ratio.end());
            
            _output_manager-&gt;outputConsole("\nSteppingstone results:");
            _output_manager-&gt;outputConsole(boost::str(boost::format("%20s %20s %20s") % "beta" % "log(ratio)" % "cumulative"));
            double log_marginal_likelihood = 0.0;
            for (auto p : log_ratio) {
                double beta = p.first;
                double logratio = p.second;
                log_marginal_likelihood += logratio;
                _output_manager-&gt;outputConsole(boost::str(boost::format("%20.5f %20.5f %20.5f") % beta % logratio % log_marginal_likelihood));
            }
            _output_manager-&gt;outputConsole(boost::str(boost::format("\nlog(marginal likelihood) = %.5f") % log_marginal_likelihood));
        }
    }

This function first creates a vector of pairs in which the first element is the heating power of the chain and the second element is the steppingstone ratio calculated for that chain. The next step is to sort this vector so that elements are in order of heating power. This sorting is purely cosmetic: sorting causes the output to have increasing chain powers and not a random ordering.

I mentioned in the introduction that the estimated marginal likelihood is simply the product of the ratios calculated from each steppingstone power posterior distribution. On the log scale, this amounts to adding up the log ratios computed by each chain. This simple summation is performed in the last section of this function while outputting the heating powers and the estimated log ratios from each chain.

Miscellany

Before moving on to making necessary modifications and additions to the Chain class, I want to finish up work on the Strom class.

Because it makes no sense to swap chains during steppingstone, we need to be sure to turn off chain swapping if steppingstone is being performed. In the swapChains function, this is as simple as changing if (_num_chains == 1) return; to if (_num_chains == 1 || _steppingstone) return; in the first statement inside the function body:

    inline void Strom::swapChains() {   
        <span style="color:#0000ff"><strong>//if (_num_chains == 1)</strong></span>
        <span style="color:#0000ff"><strong>if (_num_chains == 1 || _steppingstone)</strong></span>
            return;
            
        // Select two chains at random to swap
        // If _num_chains = 3...
        //  i  j  = (i + 1 + randint(0,1)) % _num_chains
        // ---------------------------------------------
        //  0  1  = (0 + 1 +      0      ) %     3
        //     2  = (0 + 1 +      1      ) %     3
        // ---------------------------------------------
        //  1  2  = (1 + 1 +      0      ) %     3
        //     0  = (1 + 1 +      1      ) %     3
        // ---------------------------------------------
        //  2  0  = (2 + 1 +      0      ) %     3
        //     1  = (2 + 1 +      1      ) %     3
        // ---------------------------------------------
        unsigned i = (unsigned)_lot-&gt;randint(0, _num_chains-1);
        unsigned j = i + 1 + (unsigned)_lot-&gt;randint(0, _num_chains-2);
        j %= _num_chains;

        assert(i != j && i &gt;=0 && i &lt; _num_chains && j &gt;= 0 && j &lt; _num_chains);

        // Determine upper and lower triangle cells in _swaps vector
        unsigned smaller = _num_chains;
        unsigned larger  = _num_chains;
        double index_i   = _chains[i].getChainIndex();
        double index_j   = _chains[j].getChainIndex();
        if (index_i &lt; index_j) {
            smaller = index_i;
            larger  = index_j;
        }
        else {
            smaller = index_j;
            larger  = index_i;
        }
        unsigned upper = smaller*_num_chains + larger;
        unsigned lower = larger*_num_chains  + smaller;
        _swaps[upper]++;

        // Propose swap of chains i and j
        // Proposed state swap will be successful if a uniform random deviate is less than R, where
        //    R = Ri * Rj = (Pi(j) / Pi(i)) * (Pj(i) / Pj(j))
        // Chain i: power = a, kernel = pi
        // Chain j: power = b, kernel = pj
        //      pj^a         pi^b
        // Ri = ----    Rj = ----
        //      pi^a         pj^b
        // log R = (a-b) [log(pj) - log(pi)]

        double heat_i       = _chains[i].getHeatingPower();
        double log_kernel_i = _chains[i].calcLogLikelihood() + _chains[i].calcLogJointPrior();

        double heat_j       = _chains[j].getHeatingPower();
        double log_kernel_j = _chains[j].calcLogLikelihood() + _chains[j].calcLogJointPrior();

        double logR = (heat_i - heat_j)*(log_kernel_j - log_kernel_i);

        double logu = _lot-&gt;logUniform();
        if (logu &lt; logR) {
            // accept swap
            _swaps[lower]++;
            _chains[j].setHeatingPower(heat_i);
            _chains[i].setHeatingPower(heat_j);
            _chains[j].setChainIndex(index_i);
            _chains[i].setChainIndex(index_j);
            std::vector&lt;double&gt; lambdas_i = _chains[i].getLambdas();
            std::vector&lt;double&gt; lambdas_j = _chains[j].getLambdas();
            _chains[i].setLambdas(lambdas_j);
            _chains[j].setLambdas(lambdas_i);
        }
    }   

Make the same kind of modification to the swapSummary function:

    inline void Strom::swapSummary() const { 
        <span style="color:#0000ff"><strong>//if (_num_chains &gt; 1) {</strong></span>
        <span style="color:#0000ff"><strong>if (_num_chains &gt; 1 && !_steppingstone) {</strong></span>
            unsigned i, j;
            std::cout &lt;&lt; "\nSwap summary (upper triangle = no. attempted swaps; lower triangle = no. successful swaps):" &lt;&lt; std::endl;

            // column headers
            std::cout &lt;&lt; boost::str(boost::format("%12s") % " ");
            for (i = 0; i &lt; _num_chains; ++i)
                std::cout &lt;&lt; boost::str(boost::format(" %12d") % i);
            std::cout &lt;&lt; std::endl;

            // top line
            std::cout &lt;&lt; boost::str(boost::format("%12s") % "------------");
            for (i = 0; i &lt; _num_chains; ++i)
                std::cout &lt;&lt; boost::str(boost::format("-%12s") % "------------");
            std::cout &lt;&lt; std::endl;

            // table proper
            for (i = 0; i &lt; _num_chains; ++i) {
                std::cout &lt;&lt; boost::str(boost::format("%12d") % i);
                for (j = 0; j &lt; _num_chains; ++j) {
                    if (i == j)
                        std::cout &lt;&lt; boost::str(boost::format(" %12s") % "---");
                    else
                        std::cout &lt;&lt; boost::str(boost::format(" %12.5f") % _swaps[i*_num_chains + j]);
                }
                std::cout &lt;&lt; std::endl;
            }

            // bottom line
            std::cout &lt;&lt; boost::str(boost::format("%12s") % "------------");
            for (i = 0; i &lt; _num_chains; ++i)
                std::cout &lt;&lt; boost::str(boost::format("-%12s") % "------------");
            std::cout &lt;&lt; std::endl;
        }
    }  

In order to calculate their particular steppingstone ratio, each chain must know not only the heating power that it is using but also the heating power of the power posterior distribution that is nested just inside it. The chains are told not only their heating power but this “next heating power” in the initChains function:

    inline void Strom::initChains() {   
        //...  
        // Create _num_chains chains
        _chains.resize(_num_chains);
        
        // Create _num_chains by _num_chains swap matrix
        _swaps.assign(_num_chains*_num_chains, 0);

        // Create heating power vector
        _heating_powers.assign(_num_chains, 1.0);
        calcHeatingPowers();
        
        // Initialize chains
        for (unsigned chain_index = 0; chain_index &lt; _num_chains; ++chain_index) {
            auto & c        = _chains[chain_index];
            auto likelihood = _likelihoods[chain_index];
            auto m          = likelihood-&gt;getModel();
            
            // Finish setting up models
            m-&gt;setTopologyPriorOptions(_allow_polytomies, _resolution_class_prior, _topo_prior_C);
            m-&gt;setSubsetNumPatterns(_data-&gt;calcNumPatternsVect());
            m-&gt;setSubsetSizes(_partition-&gt;calcSubsetSizes());
            m-&gt;activate();
            if (chain_index == 0)
                std::cout &lt;&lt; "\n" &lt;&lt; m-&gt;describeModel() &lt;&lt; std::endl;
            else
                m-&gt;describeModel();
                
            // Finish setting up likelihoods
            likelihood-&gt;setData(_data);
            likelihood-&gt;useUnderflowScaling(_use_underflow_scaling);
            likelihood-&gt;initBeagleLib();
            likelihood-&gt;useStoredData(_using_stored_data);
            
            // Build list of updaters, one for each free parameter in the model
            unsigned num_free_parameters = c.createUpdaters(m, _lot, likelihood);
            if (num_free_parameters == 0)
                throw XStrom("MCMC skipped because there are no free parameters in the model");

            // Tell the chain that it should adapt its updators (at least initially)
            c.startTuning();

            // Set heating power to precalculated value
            c.setChainIndex(chain_index);
            c.setHeatingPower(_heating_powers[chain_index]);
            <span style="color:#0000ff"><strong>if (_steppingstone) {</strong></span>
                <span style="color:#0000ff"><strong>if (chain_index == _num_chains - 1)</strong></span>
                    <span style="color:#0000ff"><strong>c.setNextHeatingPower(1.0);</strong></span>
                <span style="color:#0000ff"><strong>else</strong></span>
                    <span style="color:#0000ff"><strong>c.setNextHeatingPower(_heating_powers[chain_index + 1]);</strong></span>
            <span style="color:#0000ff"><strong>}</strong></span>
                        
            // Give the chain a starting tree
            std::string newick = _tree_summary-&gt;getNewick(m-&gt;getTreeIndex());
            c.setTreeFromNewick(newick);

            // Print headers in output files and make sure each updator has its starting value
            c.start();
        }
    }   

The last modifications we need to make in the Strom class are to the run function: (1) add the calcMarginalLikelihood call and (2) avoid opening the parameter or tree output files if steppingstone is in progress. Note that calcMarginalLikelihood is a no-op (i.e. no-operation, nothing done) if we are not currently carrying out steppingstone.

    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 {
            readData();
            readTrees();
            showPartitionInfo();

            // Create a Lot object that generates (pseudo)random numbers
            _lot = Lot::SharedPtr(new Lot);
            _lot-&gt;setSeed(_random_seed);

            // Create  Chain objects
            initChains();
            
            showBeagleInfo();
            showMCMCInfo();

            // Create an output manager and open output files
            _output_manager.reset(new OutputManager);
            _output_manager-&gt;outputConsole(boost::str(boost::format("\n%12s %12s %12s %12s %12s") % "iteration" % "m" % "logLike" % "logPrior" % "TL"));
            <span style="color:#0000ff"><strong>if (!_steppingstone) {</strong></span>
                _output_manager-&gt;openTreeFile("trees.tre", _data);
                _output_manager-&gt;openParameterFile("params.txt", _chains[0].getModel());
            <span style="color:#0000ff"><strong>}</strong></span>
            sample(0, _chains[0]);
            
            // Burn-in the chains
            startTuningChains();
            for (unsigned iteration = 1; iteration &lt;= _num_burnin_iter; ++iteration) {
                stepChains(iteration, false);
                swapChains();
            }
            stopTuningChains();

            // Sample the chains
            for (unsigned iteration = 1; iteration &lt;= _num_iter; ++iteration) {
                stepChains(iteration, true);
                swapChains();
            }
            showChainTuningInfo();
            stopChains();
            
            // Create swap summary
            swapSummary();
            
            <span style="color:#0000ff"><strong>// Estimate the marginal likelihood if doing steppingstone</strong></span>
            <span style="color:#0000ff"><strong>calcMarginalLikelihood();</strong></span>
            
            // Close output files
            <span style="color:#0000ff"><strong>if (!_steppingstone) {</strong></span>
                _output_manager-&gt;closeTreeFile();
                _output_manager-&gt;closeParameterFile();
            <span style="color:#0000ff"><strong>}</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;
    } 

< 20.0 | 20.1 | 20.2 >