Metropolis.hpp

Go to the documentation of this file.

/*******************************************************************************
 Causal Dynamical Triangulations in C++ using CGAL
 
 Copyright © 2015 Adam Getchell
 ******************************************************************************/
 
 
#ifndef INCLUDE_METROPOLIS_HPP_
#define INCLUDE_METROPOLIS_HPP_
 
// CDT headers
#include "Move_command.hpp"
#include "Move_strategy.hpp"
#include "S3Action.hpp"
 
using Gmpzf = CGAL::Gmpzf;
 
template <typename ManifoldType>
class MoveStrategy<Strategies::METROPOLIS, ManifoldType>
{
  using Counter = move_tracker::MoveTracker<ManifoldType>;
 
  long double m_Alpha{};
 
  long double m_K{};
 
  long double m_Lambda{};
 
  Int_precision m_passes{1};
 
  Int_precision m_checkpoint{1};
 
  Geometry<ManifoldType::dimension> m_geometry;
 
  Counter m_proposed_moves;
 
  Counter m_accepted_moves;
 
  Counter m_rejected_moves;
 
  Counter m_attempted_moves;
 
  Counter m_succeeded_moves;
 
  Counter m_failed_moves;
 
 public:
  MoveStrategy() = default;
 
  [[maybe_unused]] MoveStrategy(long double const Alpha, long double const K,
                                long double const   Lambda,
                                Int_precision const passes,
                                Int_precision const checkpoint)
      : m_Alpha(Alpha)
      , m_K(K)
      , m_Lambda(Lambda)
      , m_passes(passes)
      , m_checkpoint{checkpoint}
  {
#ifndef NDEBUG
    spdlog::debug("{} called.\n", __PRETTY_FUNCTION__);
#endif
  }
 
  [[nodiscard]] auto Alpha() const noexcept { return m_Alpha; }
 
  [[nodiscard]] auto K() const noexcept { return m_K; }
 
  [[nodiscard]] auto Lambda() const noexcept { return m_Lambda; }
 
  [[nodiscard]] auto passes() const noexcept { return m_passes; }
 
  [[nodiscard]] auto checkpoint() const noexcept { return m_checkpoint; }
 
  auto get_proposed() const { return m_proposed_moves; }
 
  auto get_accepted() const { return m_accepted_moves; }
 
  auto get_rejected() const { return m_rejected_moves; }
 
  auto get_attempted() const { return m_attempted_moves; }
 
  auto get_succeeded() const { return m_succeeded_moves; }
 
  auto get_failed() const { return m_failed_moves; }
 
  auto CalculateA1(move_tracker::move_type move) const noexcept
  {
    auto all_moves = m_proposed_moves.total();
    auto this_move = m_proposed_moves[move];
    // Set precision for initialization and assignment functions
    mpfr_set_default_prec(PRECISION);
 
    // Initialize for MPFR
    mpfr_t r_1;
    mpfr_t r_2;
    mpfr_t a_1;
    mpfr_inits2(PRECISION, a_1, nullptr);  // NOLINT
 
    mpfr_init_set_si(r_1, this_move, MPFR_RNDD);  // r_1 = this_move NOLINT
    mpfr_init_set_si(r_2, all_moves, MPFR_RNDD);  // r_2 = total_moves NOLINT
 
    // The result
    mpfr_div(a_1, r_1, r_2, MPFR_RNDD);  // a_1 = r_1/r_2 NOLINT
 
    // Convert mpfr_t total to Gmpzf result by using Gmpzf(double d)
    auto result = mpfr_get_ld(a_1, MPFR_RNDD);  // NOLINT
 
    // Free memory
    mpfr_clears(r_1, r_2, a_1, nullptr);  // NOLINT
 
#ifndef NDEBUG
    spdlog::debug("{} called.\n", __PRETTY_FUNCTION__);
    spdlog::trace("Total proposed moves = {}\n", all_moves);
    spdlog::trace("A1 is {}\n", result);
#endif
 
    return result;
  }  // CalculateA1()
 
  template <int dimension>
  auto CalculateA2(move_tracker::move_type const move) const noexcept
  {
    if (dimension == 3)
    {
      auto currentS3Action =
          S3_bulk_action(m_geometry.N1_TL, m_geometry.N3_31_13,
                         m_geometry.N3_22, m_Alpha, m_K, m_Lambda);
      auto newS3Action = static_cast<Gmpzf>(0);
      switch (move)
      {
        case move_tracker::move_type::TWO_THREE:
          // A (2,3) move adds a timelike edge
          // and a (2,2) simplex
          newS3Action =
              S3_bulk_action(m_geometry.N1_TL + 1, m_geometry.N3_31_13,
                             m_geometry.N3_22 + 1, m_Alpha, m_K, m_Lambda);
          break;
        case move_tracker::move_type::THREE_TWO:
          // A (3,2) move removes a timelike edge
          // and a (2,2) simplex
          newS3Action =
              S3_bulk_action(m_geometry.N1_TL - 1, m_geometry.N3_31_13,
                             m_geometry.N3_22 - 1, m_Alpha, m_K, m_Lambda);
          break;
        case move_tracker::move_type::TWO_SIX:
          // A (2,6) move adds 2 timelike edges and
          // 2 (1,3) and 2 (3,1) simplices
          newS3Action =
              S3_bulk_action(m_geometry.N1_TL + 2, m_geometry.N3_31_13 + 4,
                             m_geometry.N3_22, m_Alpha, m_K, m_Lambda);
          break;
        case move_tracker::move_type::SIX_TWO:
          // A (6,2) move removes 2 timelike edges and
          // 2 (1,3) and 2 (3,1) simplices
          newS3Action =
              S3_bulk_action(m_geometry.N1_TL - 2, m_geometry.N3_31_13,
                             m_geometry.N3_22 - 4, m_Alpha, m_K, m_Lambda);
          break;
        case move_tracker::move_type::FOUR_FOUR:
          // A (4,4) move changes nothing with respect to the action,
          // and e^0==1
#ifndef NDEBUG
          spdlog::trace("A2 is 1\n");
#endif
          return static_cast<long double>(1);
        default: break;
      }
 
      auto exponent        = currentS3Action - newS3Action;
      auto exponent_double = utilities::Gmpzf_to_double(exponent);
 
      // if exponent > 0 then e^exponent >=1 so according to Metropolis
      // algorithm return A2=1
      if (exponent >= 0) { return static_cast<long double>(1); }
 
      // Set precision for initialization and assignment functions
      mpfr_set_default_prec(PRECISION);
 
      // Initialize for MPFR
      mpfr_t r_1;
      mpfr_t a_2;
      mpfr_inits2(PRECISION, a_2, nullptr);  // NOLINT
 
      // Set input parameters and constants to mpfr_t equivalents
      mpfr_init_set_ld(r_1, exponent_double,  // NOLINT
                       MPFR_RNDD);            // r1 = exponent
 
      // e^exponent
      mpfr_exp(a_2, r_1, MPFR_RNDD);  // NOLINT
 
      // Convert mpfr_t total to Gmpzf result by using Gmpzf(double d)
      auto result = mpfr_get_ld(a_2, MPFR_RNDD);  // NOLINT
 
      // Free memory
      mpfr_clears(r_1, a_2, nullptr);  // NOLINT
 
#ifndef NDEBUG
      spdlog::trace("A2 is {}\n", result);
#endif
 
      return result;
    }
  }  // CalculateA2()
 
  auto try_move(move_tracker::move_type const move) -> bool
  {
    // Record the proposed move
    ++m_proposed_moves[as_integer(move)];
 
    // Calculate probability
    auto a_1 = CalculateA1(move);
 
    // Make move if random number < probability
    auto a_2 = CalculateA2<3>(move);
 
    if (auto const trial_value = utilities::generate_probability();
        trial_value <= a_1 * a_2)
    {
#ifndef NDEBUG
      spdlog::debug("{} called.\n", __PRETTY_FUNCTION__);
      spdlog::trace("Trying move.\n");
      spdlog::trace("Move type = {}\n", as_integer(move));
      spdlog::trace("Trial_value = {}\n", trial_value);
      spdlog::trace("A1 = {}\n", a_1);
      spdlog::trace("A2 = {}\n", a_2);
      spdlog::trace("A1*A2 = {}\n", a_1 * a_2);
      spdlog::trace("{}\n", trial_value <= a_1 * a_2 ? "Move accepted."
                                                     : "Move rejected.");
#endif
      // Accept the move
      ++m_accepted_moves[as_integer(move)];
      return true;
    }
 
    // Reject the move
    ++m_rejected_moves[as_integer(move)];
    return false;
 
  }  // try_move()
 
  [[nodiscard]] auto initialize(ManifoldType t_manifold)
      -> std::optional<MoveCommand<ManifoldType>>
  try
  {
    MoveCommand command(t_manifold);
    fmt::print("Making initial moves ...\n");
 
    // Make a move of each type
    for (auto move = 0;
         move < move_tracker::moves_per_dimension(ManifoldType::dimension);
         ++move)
    {
#ifndef NDEBUG
      spdlog::trace("Making move {} ...\n", move);
#endif
      command.enqueue(move_tracker::as_move(move));
      ++m_proposed_moves[move];
      ++m_accepted_moves[move];
    }
 
    // Execute the initial moves
    command.execute();
    command.print_successful();
    command.print_errors();
 
    // Update attempted, succeeded, and failed moves from MoveCommand
    m_attempted_moves += command.get_attempted();
    m_succeeded_moves += command.get_succeeded();
    m_failed_moves += command.get_failed();
 
    // Reset the move counters
    command.reset_counters();
 
    // print initial results
    auto initial_results = command.get_results();
    initial_results.print();
    initial_results.print_details();
 
    return command;
  }
  catch (std::system_error const& SystemError)
  {
    spdlog::debug("Metropolis initialization failed with {} ... exiting.\n",
                  SystemError.what());
    spdlog::trace("{}\n", SystemError.code().message());
    return std::nullopt;
  }
 
  auto operator()(ManifoldType const& t_manifold) -> ManifoldType
  {
#ifndef NDEBUG
    spdlog::debug("{} called.\n", __PRETTY_FUNCTION__);
#endif
 
    fmt::print(
        "Starting Metropolis-Hastings algorithm in {}+1 dimensions ...\n",
        ManifoldType::dimension - 1);
 
    auto command = initialize(t_manifold).value_or(MoveCommand(t_manifold));
 
    fmt::print("Making random moves ...\n");
 
    // Loop through m_passes
    for (auto pass_number = 1; pass_number <= m_passes; ++pass_number)
    {
      fmt::print("=== Pass {} ===\n", pass_number);
      auto total_simplices_this_pass = command.get_const_results().N3();
      // Attempt a random move per simplex
      for (auto move_attempt = 0; move_attempt < total_simplices_this_pass;
           ++move_attempt)
      {
        // Pick a move to attempt
 
        if (auto move = move_tracker::generate_random_move_3(); try_move(move))
        {
          command.enqueue(move);
        }
      }  // Ends loop through CurrentTotalSimplices
 
      // Do the moves
      command.execute();
 
      // Update attempted and failed moves
      this->m_attempted_moves += command.get_attempted();
      this->m_succeeded_moves += command.get_succeeded();
      this->m_failed_moves += command.get_failed();
 
      // Do stuff on checkpoint
      if (pass_number % m_checkpoint == 0)
      {
        fmt::print("=== Pass {} ===\n", pass_number);
        fmt::print("Writing to file.\n");
        print_results();
        auto result = command.get_results();
        utilities::write_file(result);
      }
    }  // Ends loop through m_passes
 
    // output results
    fmt::print("=== Run results ===\n");
    print_results();
    return command.get_results();
  }  // operator()
 
  void print_results()
  {
    if (ManifoldType::dimension == 3)
    {
      fmt::print("=== Move Results ===\n");
      fmt::print(
          "There were {} proposed moves with {} accepted moves and {} rejected "
          "moves.\n",
          m_proposed_moves.total(), m_accepted_moves.total(),
          m_rejected_moves.total());
      fmt::print(
          "There were {} attempted moves with {} successful moves and {} "
          "failed moves.\n",
          m_attempted_moves.total(), m_succeeded_moves.total(),
          m_failed_moves.total());
      fmt::print(
          "(2,3) moves: {} proposed ({} accepted and {} rejected) with {} "
          "attempted ({} successful and {} failed).\n",
          m_proposed_moves.two_three_moves(),
          m_accepted_moves.two_three_moves(),
          m_rejected_moves.two_three_moves(),
          m_attempted_moves.two_three_moves(),
          m_succeeded_moves.two_three_moves(),
          m_failed_moves.two_three_moves());
 
      fmt::print(
          "(3,2) moves: {} proposed ({} accepted and {} rejected) with {} "
          "attempted ({} successful and {} failed).\n",
          m_proposed_moves.three_two_moves(),
          m_accepted_moves.three_two_moves(),
          m_rejected_moves.three_two_moves(),
          m_attempted_moves.three_two_moves(),
          m_succeeded_moves.three_two_moves(),
          m_failed_moves.three_two_moves());
 
      fmt::print(
          "(2,6) moves: {} proposed ({} accepted and {} rejected) with {} "
          "attempted ({} successful and {} failed).\n",
          m_proposed_moves.two_six_moves(), m_accepted_moves.two_six_moves(),
          m_rejected_moves.two_six_moves(), m_attempted_moves.two_six_moves(),
          m_succeeded_moves.two_six_moves(), m_failed_moves.two_six_moves());
 
      fmt::print(
          "(6,2) moves: {} proposed ({} accepted and {} rejected) with {} "
          "attempted ({} successful and {} failed).\n",
          m_proposed_moves.six_two_moves(), m_accepted_moves.six_two_moves(),
          m_rejected_moves.six_two_moves(), m_attempted_moves.six_two_moves(),
          m_succeeded_moves.six_two_moves(), m_failed_moves.six_two_moves());
 
      fmt::print(
          "(4,4) moves: {} proposed ({} accepted and {} rejected) with {} "
          "attempted ({} successful and {} failed).\n",
          m_proposed_moves.four_four_moves(),
          m_accepted_moves.four_four_moves(),
          m_rejected_moves.four_four_moves(),
          m_attempted_moves.four_four_moves(),
          m_succeeded_moves.four_four_moves(),
          m_failed_moves.four_four_moves());
    }
  }  // print_results
};  // Metropolis
 
using Metropolis_3 =
    MoveStrategy<Strategies::METROPOLIS, manifolds::Manifold_3>;
using Metropolis_4 =
    MoveStrategy<Strategies::METROPOLIS, manifolds::Manifold_4>;
 
#endif  // INCLUDE_METROPOLIS_HPP_