MCData.cc

Go to the documentation of this file.

#include "casm/monte_carlo/MCData.hh"

#include "casm/external/boost.hh"

namespace CASM {

  MCDataEquilibration::MCDataEquilibration(const Eigen::VectorXd &observations, double prec) {

    size_type start1, start2, N;
    double sum1, sum2;
    double eps = (observations(0) == 0.0) ? 1e-8 : std::abs(observations(0)) * 1e-8;

    N = observations.size();
    bool is_even = ((N % 2) == 0);

    // -------------------------------------------------
    // if the values are all the same to observations(0)*1e-8, set m_is_equilibrated = true; m_equil_samples = 0;
    bool all_same = true;
    for(size_type i = 0; i < N; i++)
      if(std::abs(observations(i) - observations(0)) > eps) {
        all_same = false;
        break;
      }
    if(all_same) {
      m_is_equilibrated = true;
      m_equil_samples = 0;
      return;
    }

    // find partitions
    start1 = 0;
    start2 = (is_even) ? N / 2 : (N / 2) + 1;

    // find sums for each partition
    sum1 = observations.head(start2).sum();
    // mean1 = sum1 / (start2 - start1)
    sum2 = observations.segment(start2, N - start2).sum();
    // mean2 = sum2 / (N - start2)

    //std::cout << observations << std::endl;

    // increment start1 (and update start2, sum1, and sum2) until abs(mean1 - mean2) < prec
    while(std::abs((sum1 / (start2 - start1)) - (sum2 / (N - start2))) > prec && start1 < N - 2) {
      /*std::cout << "  start1: " << start1 << "  start2: " << start2 << "  N: " << N
                << "  mean1: " << (sum1 / (start2 - start1)) << "  mean2: " << (sum2 / (N - start2))
                << "  diff: " << std::abs((sum1 / (start2 - start1)) - (sum2 / (N - start2)))
                << "  prec: " << prec << std::endl;*/
      if(is_even) {
        sum1 -= observations(start1);
        sum1 += observations(start2);
        sum2 -= observations(start2);
        start2++;
      }
      else {
        sum1 -= observations(start1);
      }

      start1++;
      is_even = !is_even;
    }

    // ensure that there are values on either side of the total mean
    double mean_tot = (sum1 + sum2) / (N - start1);

    if(observations(start1) < mean_tot) {
      while(observations(start1) < mean_tot && start1 < N - 1)
        start1++;
    }
    else {
      while(observations(start1) > mean_tot && start1 < N - 1)
        start1++;
    }

    m_is_equilibrated = (start1 < N - 1);
    m_equil_samples = start1;
  }


  MCDataConvergence::MCDataConvergence(const Eigen::VectorXd &observations, double conf) :
    m_is_converged(false),
    m_mean(observations.mean()),
    m_squared_norm(observations.squaredNorm()) {

    // will check if Var <= criteria
    double z_alpha = sqrt(2.0) * boost::math::erf_inv(conf);

    bool found_rho;
    double rho;
    double CoVar0;
    size_type N = observations.size();

    // try to calculate variance taking into account correlations
    std::tie(found_rho, rho, CoVar0) = _calc_rho(observations);

    if(!found_rho) {
      m_is_converged = false;
      m_calculated_prec = 1.0 / 0.0;
      return;
    }


    double var_of_mean = (CoVar0 / N) * (1.0 + rho) / (1.0 - rho);

    m_calculated_prec = z_alpha * sqrt(var_of_mean);
  }

  double covariance(const Eigen::VectorXd &X, const Eigen::VectorXd &Y, double mean) {
    Eigen::VectorXd vmean = Eigen::VectorXd::Ones(X.size()) * mean;
    return (X - vmean).dot(Y - vmean) / X.size();
  }


  std::tuple<bool, double, double> MCDataConvergence::_calc_rho(const Eigen::VectorXd &observations) {

    size_type N = observations.size();
    double CoVar0 = m_squared_norm / N - m_mean * m_mean;

    // if there is essentially no variation, return rho(l==1)
    if(std::abs(CoVar0 / m_mean) < 1e-8 || CoVar0 == 0.0) {
      CoVar0 = 0.0;
      return std::make_tuple(true, pow(2.0, -1.0 / 1), CoVar0);
    }

    // simple incremental search for now, could try bracketing / bisection
    for(size_type i = 1; i < N; ++i) {

      size_type range_size = N - i;

      double cov = covariance(observations.segment(0, range_size), observations.segment(i, range_size), m_mean);

      if(std::abs(cov / CoVar0) <= 0.5) {
        return std::make_tuple(true, pow(2.0, (-1.0 / i)), CoVar0);
      }
    }

    // if could not find:
    return std::make_tuple(false, 0.0, CoVar0);
  }

}