text: #ifndef STAN_ANALYZE_MCMC_COMPUTE_POTENTIAL_SCALE_REDUCTION_HPP
  #define STAN_ANALYZE_MCMC_COMPUTE_POTENTIAL_SCALE_REDUCTION_HPP
  
  #include <stan/math/prim.hpp>
  #include <stan/analyze/mcmc/autocovariance.hpp>
  #include <stan/analyze/mcmc/split_chains.hpp>
  #include <boost/accumulators/accumulators.hpp>
  #include <boost/accumulators/statistics/stats.hpp>
  #include <boost/accumulators/statistics/mean.hpp>
  #include <boost/accumulators/statistics/variance.hpp>
  #include <boost/math/special_functions/fpclassify.hpp>
  #include <algorithm>
  #include <cmath>
  #include <vector>
  #include <limits>
  
  namespace stan {
  namespace analyze {
  
  /**
   * Computes the potential scale reduction (Rhat) for the specified
   * parameter across all kept samples.
   *
   * See more details in Stan reference manual section "Potential
   * Scale Reduction". http://mc-stan.org/users/documentation
   *
   * Current implementation assumes draws are stored in contiguous
   * blocks of memory.  Chains are trimmed from the back to match the
   * length of the shortest chain.
   *
   * @param draws stores pointers to arrays of chains
   * @param sizes stores sizes of chains
   * @return potential scale reduction for the specified parameter
   */
  inline double compute_potential_scale_reduction(
      std::vector<const double*> draws, std::vector<size_t> sizes) {
    int num_chains = sizes.size();
    size_t num_draws = sizes[0];
    for (int chain = 1; chain < num_chains; ++chain) {
      num_draws = std::min(num_draws, sizes[chain]);
    }
  
    // check if chains are constant; all equal to first draw's value
    bool are_all_const = false;
    Eigen::VectorXd init_draw = Eigen::VectorXd::Zero(num_chains);
  
    for (int chain = 0; chain < num_chains; chain++) {
      Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, 1>> draw(
          draws[chain], sizes[chain]);
  
      for (int n = 0; n < num_draws; n++) {
        if (!boost::math::isfinite(draw(n))) {
          return std::numeric_limits<double>::quiet_NaN();
        }
      }
  
      init_draw(chain) = draw(0);
  
      if (draw.isApproxToConstant(draw(0))) {
        are_all_const |= true;
      }
    }
  
    if (are_all_const) {
      // If all chains are constant then return NaN
      // if they all equal the same constant value
      if (init_draw.isApproxToConstant(init_draw(0))) {
        return std::numeric_limits<double>::quiet_NaN();
      }
    }
  
    using boost::accumulators::accumulator_set;
    using boost::accumulators::stats;
    using boost::accumulators::tag::mean;
    using boost::accumulators::tag::variance;
  
    Eigen::VectorXd chain_mean(num_chains);
    accumulator_set<double, stats<variance>> acc_chain_mean;
    Eigen::VectorXd chain_var(num_chains);
    double unbiased_var_scale = num_draws / (num_draws - 1.0);
  
    for (int chain = 0; chain < num_chains; ++chain) {
      accumulator_set<double, stats<mean, variance>> acc_draw;
      for (int n = 0; n < num_draws; ++n) {
        acc_draw(draws[chain][n]);
      }
  
      chain_mean(chain) = boost::accumulators::mean(acc_draw);
      acc_chain_mean(chain_mean(chain));
      chain_var(chain)
          = boost::accumulators::variance(acc_draw) * unbiased_var_scale;
    }
  
    double var_between = num_draws * boost::accumulators::variance(acc_chain_mean)
                         * num_chains / (num_chains - 1);
    double var_within = chain_var.mean();
  
    // rewrote [(n-1)*W/n + B/n]/W as (n-1+ B/W)/n
    return sqrt((var_between / var_within + num_draws - 1) / num_draws);
  }
  
  /**
   * Computes the potential scale reduction (Rhat) for the specified
   * parameter across all kept samples.
   *
   * See more details in Stan reference manual section "Potential
   * Scale Reduction". http://mc-stan.org/users/documentation
   *
   * Current implementation assumes draws are stored in contiguous
   * blocks of memory.  Chains are trimmed from the back to match the
   * length of the shortest chain.  Argument size will be broadcast to
   * same length as draws.
   *
   * @param draws stores pointers to arrays of chains
   * @param sizes stores sizes of chains
   * @return potential scale reduction for the specified parameter
   */
  inline double compute_potential_scale_reduction(
      std::vector<const double*> draws, size_t size) {
    int num_chains = draws.size();
    std::vector<size_t> sizes(num_chains, size);
    return compute_potential_scale_reduction(draws, sizes);
  }
  
  /**
   * Computes the split potential scale reduction (Rhat) for the
   * specified parameter across all kept samples.  When the number of
   * total draws N is odd, the (N+1)/2th draw is ignored.
   *
   * See more details in Stan reference manual section "Potential
   * Scale Reduction". http://mc-stan.org/users/documentation
   *
   * Current implementation assumes draws are stored in contiguous
   * blocks of memory.  Chains are trimmed from the back to match the
   * length of the shortest chain.
   *
   * @param draws stores pointers to arrays of chains
   * @param sizes stores sizes of chains
   * @return potential scale reduction for the specified parameter
   */
  inline double compute_split_potential_scale_reduction(
      std::vector<const double*> draws, std::vector<size_t> sizes) {
    int num_chains = sizes.size();
    size_t num_draws = sizes[0];
    for (int chain = 1; chain < num_chains; ++chain) {
      num_draws = std::min(num_draws, sizes[chain]);
    }
  
    std::vector<const double*> split_draws = split_chains(draws, sizes);
  
    double half = num_draws / 2.0;
    std::vector<size_t> half_sizes(2 * num_chains, std::floor(half));
  
    return compute_potential_scale_reduction(split_draws, half_sizes);
  }
  
  /**
   * Computes the split potential scale reduction (Rhat) for the
   * specified parameter across all kept samples.  When the number of
   * total draws N is odd, the (N+1)/2th draw is ignored.
   *
   * See more details in Stan reference manual section "Potential
   * Scale Reduction". http://mc-stan.org/users/documentation
   *
   * Current implementation assumes draws are stored in contiguous
   * blocks of memory.  Chains are trimmed from the back to match the
   * length of the shortest chain.  Argument size will be broadcast to
   * same length as draws.
   *
   * @param draws stores pointers to arrays of chains
   * @param sizes stores sizes of chains
   * @return potential scale reduction for the specified parameter
   */
  inline double compute_split_potential_scale_reduction(
      std::vector<const double*> draws, size_t size) {
    int num_chains = draws.size();
    std::vector<size_t> sizes(num_chains, size);
    return compute_split_potential_scale_reduction(draws, sizes);
  }
  
  }  // namespace analyze
  }  // namespace stan
  
  #endif
  