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_
 
#include <cmath>
#include <cstdint>
#include <expected>
#include <optional>
#include <stdexcept>
#include <utility>
 
// CDT headers
#include "Ergodic_moves_3.hpp"
#include "Move_run.hpp"
#include "Move_strategy.hpp"
#include "Random.hpp"
#include "S3Action.hpp"
#include "Utilities.hpp"
 
namespace cdt
{
  template <typename ManifoldType>
    requires(ManifoldType::dimension == 3)
  class MoveStrategy<MoveStrategyKind::METROPOLIS, ManifoldType>
  {
    using Counter        = move_tracker::MoveTracker;
    using CommandResults = detail::MoveCommandResults<ManifoldType>;
 
    struct RunStatistics
    {
      Geometry<ManifoldType::dimension> geometry;

      std::uint64_t transition_trace{14695981039346656037ULL};

      std::uint64_t transition_count{};

      Counter proposed;

      Counter accepted;

      Counter rejected;
    };
 
    using PassResult = detail::MovePassResult<ManifoldType, RunStatistics>;

    s3_action::PhysicalParameters m_parameters;

    MoveRunCadence m_cadence;

    bool m_write_files{true};

    cdt::Random m_generator{
        cdt::Random{}.split(cdt::random_streams::transitions)};

    utilities::Reproducibility_metadata m_reproducibility;

    CommandResults m_command_results;

    RunStatistics m_run_statistics;

    Int_precision m_checkpoint_events{};

    Int_precision m_completed_passes{};

    bool                      m_resume_pending{};
 
    [[nodiscard]] static auto to_counts(Counter const& counter)
        -> utilities::Move_statistics::Counts
    {
      utilities::Move_statistics::Counts counts{};
      for (std::size_t index = 0; index < counts.size(); ++index)
      {
        counts[index] = counter[static_cast<gsl::index>(index)];
      }
      return counts;
    }
 
    [[nodiscard]] static auto from_counts(
        utilities::Move_statistics::Counts const& counts) -> Counter
    {
      Counter counter;
      for (std::size_t index = 0; index < counts.size(); ++index)
      {
        counter[static_cast<gsl::index>(index)] = counts[index];
      }
      return counter;
    }
 
    void restore_statistics(utilities::Reproducibility_metadata const& metadata)
    {
      if (!metadata.move_statistics || !metadata.transition_trace ||
          !metadata.transition_count)
      {
        throw std::invalid_argument(
            "Checkpoint resume requires cumulative transition statistics.");
      }
      auto const& saved                 = *metadata.move_statistics;
      m_command_results.attempted       = from_counts(saved.attempted);
      m_command_results.succeeded       = from_counts(saved.succeeded);
      m_command_results.failed          = from_counts(saved.failed);
      m_run_statistics.proposed         = from_counts(saved.proposed);
      m_run_statistics.accepted         = from_counts(saved.accepted);
      m_run_statistics.rejected         = from_counts(saved.rejected);
      m_run_statistics.transition_trace = *metadata.transition_trace;
      m_run_statistics.transition_count = *metadata.transition_count;
    }
 
    static void record_transition(
        RunStatistics& statistics, move_tracker::MoveType const move,
        ergodic_moves::MoveOutcome const outcome) noexcept
    {
      auto const append = [&statistics](std::uint8_t const value) {
        statistics.transition_trace ^= value;
        statistics.transition_trace *= 1099511628211ULL;
      };
      append(static_cast<std::uint8_t>(move));
      append(static_cast<std::uint8_t>(outcome));
      ++statistics.transition_count;
    }
 
   public:
    MoveStrategy() = delete;

    [[maybe_unused]] MoveStrategy(long double const alpha, long double const k,
                                  long double const   lambda,
                                  Int_precision const passes,
                                  Int_precision const checkpoint,
                                  bool const          write_files = true)
        : MoveStrategy{alpha,
                       k,
                       lambda,
                       passes,
                       checkpoint,
                       write_files,
                       cdt::Random{}.split(cdt::random_streams::transitions)}
    {}

    [[maybe_unused]] MoveStrategy(
        long double const alpha, long double const k, long double const lambda,
        Int_precision const passes, Int_precision const checkpoint,
        bool const write_files, cdt::Random random,
        std::optional<utilities::Reproducibility_metadata> reproducibility =
            std::nullopt,
        Int_precision const completed_passes = 0)
        : m_parameters{s3_action::make_physical_parameters(alpha, k, lambda)}
        , m_cadence{detail::parse_move_run_cadence(passes, checkpoint,
                                                   "Metropolis")}
        , m_write_files{write_files}
        , m_generator{std::move(random)}
        , m_reproducibility{reproducibility.value_or(
              utilities::Reproducibility_metadata{
                  .seed                = m_generator.seed(),
                  .alpha               = alpha,
                  .k                   = k,
                  .lambda              = lambda,
                  .configured_passes   = passes,
                  .checkpoint_interval = checkpoint})}
        , m_completed_passes{completed_passes}
    {
      if (m_completed_passes < 0)
      {
        throw std::invalid_argument("Completed passes cannot be negative.");
      }
      m_reproducibility.alpha  = m_parameters.alpha();
      m_reproducibility.k      = m_parameters.k();
      m_reproducibility.lambda = m_parameters.lambda();
      auto const total_passes  = static_cast<std::int64_t>(m_completed_passes) +
                                 static_cast<std::int64_t>(m_cadence.passes());
      if (!std::in_range<Int_precision>(total_passes))
      {
        throw std::out_of_range(
            "Total pass count exceeds the supported range.");
      }
      if (m_completed_passes > 0 || m_reproducibility.transition_random_state)
      {
        if (!m_reproducibility.configured_passes ||
            *m_reproducibility.configured_passes != total_passes ||
            !m_reproducibility.transition_random_state)
        {
          throw std::invalid_argument(
              "Checkpoint resume state does not match its pass range.");
        }
        if (m_reproducibility.seed != m_generator.seed() ||
            m_reproducibility.transition_stream != m_generator.stream() ||
            *m_reproducibility.transition_random_state !=
                m_generator.serialized_state())
        {
          throw std::invalid_argument(
              "Checkpoint resume generator does not match its recorded random state.");
        }
        restore_statistics(m_reproducibility);
        m_resume_pending = true;
      }
      m_reproducibility.seed              = m_generator.seed();
      m_reproducibility.transition_stream = m_generator.stream();
      m_reproducibility.configured_passes =
          static_cast<Int_precision>(total_passes);
      m_reproducibility.checkpoint_interval = m_cadence.checkpoint();
      m_reproducibility.transition_random_state.reset();
#ifndef NDEBUG
      spdlog::debug("{} called.\n", CDT_PRETTY_FUNCTION);
#endif
    }

    MoveStrategy(long double const alpha, long double const k,
                 long double const lambda, Int_precision const passes,
                 Int_precision const checkpoint, bool const write_files,
                 cdt::RandomSeed const seed)
        : MoveStrategy{
              alpha,
              k,
              lambda,
              passes,
              checkpoint,
              write_files,
              cdt::Random{seed, cdt::random_streams::transitions}
    }
    {}

    [[nodiscard]] auto alpha() const noexcept { return m_parameters.alpha(); }

    [[nodiscard]] auto k() const noexcept { return m_parameters.k(); }

    [[nodiscard]] auto lambda() const noexcept { return m_parameters.lambda(); }

    [[nodiscard]] auto passes() const noexcept { return m_cadence.passes(); }

    [[nodiscard]] auto checkpoint() const noexcept
    { return m_cadence.checkpoint(); }

    [[nodiscard]] auto checkpoint_events() const noexcept
    { return m_checkpoint_events; }

    [[nodiscard]] auto writes_files() const noexcept { return m_write_files; }

    [[nodiscard]] auto seed() const noexcept { return m_generator.seed(); }

    [[nodiscard]] auto stream() const noexcept { return m_generator.stream(); }

    [[nodiscard]] auto transition_trace() const noexcept
    { return m_run_statistics.transition_trace; }

    [[nodiscard]] auto transition_count() const noexcept
    { return m_run_statistics.transition_count; }

    [[nodiscard]] auto reproducibility_metadata(
        ManifoldType const& manifold, utilities::ArtifactKind const artifact,
        Int_precision const completed_passes) const
        -> utilities::Reproducibility_metadata
    {
      return make_reproducibility_metadata(manifold, artifact, completed_passes,
                                           m_command_results, m_run_statistics);
    }

    [[nodiscard]] auto proposed() const noexcept -> Counter const&
    { return m_run_statistics.proposed; }

    [[nodiscard]] auto accepted() const noexcept -> Counter const&
    { return m_run_statistics.accepted; }

    [[nodiscard]] auto rejected() const noexcept -> Counter const&
    { return m_run_statistics.rejected; }

    [[nodiscard]] auto attempted() const noexcept -> Counter const&
    { return m_command_results.attempted; }

    [[nodiscard]] auto succeeded() const noexcept -> Counter const&
    { return m_command_results.succeeded; }

    [[nodiscard]] auto failed() const noexcept -> Counter const&
    { return m_command_results.failed; }

    [[nodiscard]] auto geometry() const noexcept
        -> Geometry<ManifoldType::dimension> const&
    { return m_run_statistics.geometry; }

    [[nodiscard]] static constexpr auto reverse_move(
        move_tracker::MoveType const move) noexcept
        -> std::optional<move_tracker::MoveType>
    {
      using enum move_tracker::MoveType;
      switch (move)
      {
        case TWO_THREE: return THREE_TWO;
        case THREE_TWO: return TWO_THREE;
        case TWO_SIX: return SIX_TWO;
        case SIX_TWO: return TWO_SIX;
        case FOUR_FOUR: return FOUR_FOUR;
      }
      return std::nullopt;
    }

    [[nodiscard]] static constexpr auto proposal_site_count(
        Geometry<ManifoldType::dimension> const& geometry,
        move_tracker::MoveType const             move) noexcept -> Int_precision
    {
      using enum move_tracker::MoveType;
      switch (move)
      {
        case TWO_THREE: return geometry.N3_22;
        case THREE_TWO: return geometry.N1_TL;
        case TWO_SIX: return geometry.N3_13;
        case SIX_TWO: return geometry.N0;
        case FOUR_FOUR: return geometry.N1_SL;
      }
      return 0;
    }

    [[nodiscard]] static auto proposal_probability(
        Geometry<ManifoldType::dimension> const& geometry,
        move_tracker::MoveType const             move) -> mpfr_values::Value
    {
      auto const sites = proposal_site_count(geometry, move);
      if (sites <= 0) { return mpfr_values::zero(); }
      auto const move_count =
          mpfr_values::from_integer(move_tracker::NUMBER_OF_3D_MOVES);
      auto const site_count  = mpfr_values::from_integer(sites);
      auto const denominator = mpfr_values::multiply(move_count, site_count);
      return mpfr_values::divide(mpfr_values::from_integer(1), denominator);
    }

    [[nodiscard]] static auto hastings_ratio(
        Geometry<ManifoldType::dimension> const& current,
        Geometry<ManifoldType::dimension> const& proposed,
        move_tracker::MoveType const             move) -> mpfr_values::Value
    {
      auto const forward      = proposal_probability(current, move);
      auto const reverse_type = reverse_move(move);
      if (!reverse_type)
      {
        throw std::invalid_argument{"Cannot reverse an unknown move type."};
      }
      auto const reverse = proposal_probability(proposed, *reverse_type);
      if (mpfr_zero_p(forward.fr()) != 0 || mpfr_zero_p(reverse.fr()) != 0)
      {
        throw std::logic_error(
            "A successful reversible proposal must have nonzero forward and reverse probabilities.");
      }
      return mpfr_values::divide(reverse, forward);
    }

    [[nodiscard]] auto action_ratio(
        Geometry<ManifoldType::dimension> const& current,
        Geometry<ManifoldType::dimension> const& proposed) const
        -> mpfr_values::Value
    {
      auto const current_action = s3_action::s3_bulk_action(
          current.N1_TL, current.N3_31_13, current.N3_22, m_parameters);
      auto const proposed_action = s3_action::s3_bulk_action(
          proposed.N1_TL, proposed.N3_31_13, proposed.N3_22, m_parameters);
      return mpfr_values::exponential(
          mpfr_values::subtract(current_action, proposed_action));
    }

    [[nodiscard]] auto acceptance_probability(
        Geometry<ManifoldType::dimension> const& current,
        Geometry<ManifoldType::dimension> const& proposed,
        move_tracker::MoveType const move) const -> mpfr_values::Value
    {
      auto const ratio =
          mpfr_values::multiply(hastings_ratio(current, proposed, move),
                                action_ratio(current, proposed));
      auto const one = mpfr_values::from_integer(1);
      return mpfr_cmp(ratio.fr(), one.fr()) < 0 ? ratio : one;
    }
 
   private:
    [[nodiscard]] auto propose_candidate(ManifoldType const&          current,
                                         move_tracker::MoveType const move)
        -> ergodic_moves::MoveResult<ManifoldType>
    {
      using enum move_tracker::MoveType;
      switch (move)
      {
        case TWO_THREE:
          return ergodic_moves::propose_23_move(current, m_generator);
        case THREE_TWO:
          return ergodic_moves::propose_32_move(current, m_generator);
        case TWO_SIX:
          return ergodic_moves::propose_26_move(current, m_generator);
        case SIX_TWO:
          return ergodic_moves::propose_62_move(current, m_generator);
        case FOUR_FOUR:
          return ergodic_moves::propose_44_move(current, m_generator);
      }
      return std::unexpected{
          ergodic_moves::MoveError{
                                   .category       = ergodic_moves::MoveFailure::UNKNOWN_MOVE,
                                   .requested_move = move}
      };
    }
 
    [[nodiscard]] auto make_reproducibility_metadata(
        ManifoldType const& manifold, utilities::ArtifactKind const artifact,
        Int_precision const   completed_passes,
        CommandResults const& command_results,
        RunStatistics const&  statistics) const
        -> utilities::Reproducibility_metadata
    {
      auto metadata             = m_reproducibility;
      metadata.artifact         = artifact;
      metadata.completed_passes = completed_passes;
      metadata.transition_trace = statistics.transition_trace;
      metadata.transition_count = statistics.transition_count;
      metadata.move_statistics  = utilities::Move_statistics{
          .proposed  = to_counts(statistics.proposed),
          .accepted  = to_counts(statistics.accepted),
          .rejected  = to_counts(statistics.rejected),
          .attempted = to_counts(command_results.attempted),
          .succeeded = to_counts(command_results.succeeded),
          .failed    = to_counts(command_results.failed)};
      if (artifact == utilities::ArtifactKind::CHECKPOINT &&
          metadata.max_threads && *metadata.max_threads > 0)
      {
        metadata.transition_random_state =
            std::make_shared<std::string const>(m_generator.serialized_state());
      }
      else
      {
        metadata.transition_random_state.reset();
      }
      utilities::update_reproducibility_state(metadata, manifold);
      if (metadata.desired_simplices == 0)
      {
        metadata.desired_simplices = manifold.N3();
      }
      if (metadata.desired_timeslices == 0)
      {
        metadata.desired_timeslices = manifold.max_time();
      }
      return metadata;
    }
 
    auto resolve_transition(ManifoldType&                current,
                            CommandResults&              command_results,
                            RunStatistics&               statistics,
                            move_tracker::MoveType const move,
                            long double const            trial_value)
        -> ergodic_moves::MoveOutcome
    {
      if (!std::isfinite(trial_value) || trial_value < 0.0L ||
          trial_value > 1.0L)
      {
        throw std::invalid_argument("MH trial value must lie in [0, 1].");
      }
 
      if (!reverse_move(move))
      {
        return ergodic_moves::outcome_from(ergodic_moves::MoveError{
            .category       = ergodic_moves::MoveFailure::UNKNOWN_MOVE,
            .requested_move = move});
      }
 
      statistics.geometry = current.geometry();
      ++statistics.proposed[move];
      ++command_results.attempted[move];
 
      auto candidate = propose_candidate(current, move);
      if (!candidate)
      {
        ++command_results.failed[move];
        ++statistics.rejected[move];
        auto const outcome = ergodic_moves::outcome_from(candidate.error());
        record_transition(statistics, move, outcome);
        return outcome;
      }
      if (!ergodic_moves::detail::check_move(current, *candidate, move))
      {
        ++command_results.failed[move];
        ++statistics.rejected[move];
        record_transition(statistics, move,
                          ergodic_moves::MoveOutcome::EXECUTION_FAILED);
        return ergodic_moves::MoveOutcome::EXECUTION_FAILED;
      }
 
      ++command_results.succeeded[move];
      auto const probability = acceptance_probability(
          statistics.geometry, candidate->geometry(), move);
      if (mpfr_cmp_ld(probability.fr(), trial_value) >= 0)
      {
        swap(*candidate, current);
        statistics.geometry = current.geometry();
        ++statistics.accepted[move];
        record_transition(statistics, move,
                          ergodic_moves::MoveOutcome::METROPOLIS_ACCEPTED);
        return ergodic_moves::MoveOutcome::METROPOLIS_ACCEPTED;
      }
 
      ++statistics.rejected[move];
      record_transition(statistics, move,
                        ergodic_moves::MoveOutcome::METROPOLIS_REJECTED);
      return ergodic_moves::MoveOutcome::METROPOLIS_REJECTED;
    }
 
    [[nodiscard]] auto sample_transition(ManifoldType&   current,
                                         CommandResults& command_results,
                                         RunStatistics&  statistics)
        -> ergodic_moves::MetropolisTransition
    {
      auto const move = move_tracker::generate_random_move_3(m_generator);
      auto const trial_value = utilities::generate_probability(m_generator);
      return {move, resolve_transition(current, command_results, statistics,
                                       move, trial_value)};
    }
 
    [[nodiscard]] auto execute_pass(ManifoldType        current,
                                    RunStatistics       statistics,
                                    Int_precision const attempts) -> PassResult
    {
      auto command_results = CommandResults{};
      for (auto move_attempt = Int_precision{0}; move_attempt < attempts;
           ++move_attempt)
      {
        static_cast<void>(
            sample_transition(current, command_results, statistics));
      }
      return {.manifold        = std::move(current),
              .command_results = std::move(command_results),
              .strategy_state  = std::move(statistics)};
    }
 
    static void print_results(CommandResults const& command_results,
                              RunStatistics const&  statistics)
    {
      fmt::print("=== Move Results ===\n");
      fmt::print(
          "There were {} proposed moves with {} accepted moves and {} rejected "
          "moves.\n",
          statistics.proposed.total(), statistics.accepted.total(),
          statistics.rejected.total());
      fmt::print(
          "There were {} candidate construction attempts with {} successful "
          "candidates and {} failed candidates.\n",
          command_results.attempted.total(), command_results.succeeded.total(),
          command_results.failed.total());
      fmt::print(
          "(2,3) moves: {} proposed ({} accepted and {} rejected); candidate "
          "construction: {} attempted ({} succeeded and {} failed).\n",
          statistics.proposed.two_three_moves(),
          statistics.accepted.two_three_moves(),
          statistics.rejected.two_three_moves(),
          command_results.attempted.two_three_moves(),
          command_results.succeeded.two_three_moves(),
          command_results.failed.two_three_moves());
 
      fmt::print(
          "(3,2) moves: {} proposed ({} accepted and {} rejected); candidate "
          "construction: {} attempted ({} succeeded and {} failed).\n",
          statistics.proposed.three_two_moves(),
          statistics.accepted.three_two_moves(),
          statistics.rejected.three_two_moves(),
          command_results.attempted.three_two_moves(),
          command_results.succeeded.three_two_moves(),
          command_results.failed.three_two_moves());
 
      fmt::print(
          "(2,6) moves: {} proposed ({} accepted and {} rejected); candidate "
          "construction: {} attempted ({} succeeded and {} failed).\n",
          statistics.proposed.two_six_moves(),
          statistics.accepted.two_six_moves(),
          statistics.rejected.two_six_moves(),
          command_results.attempted.two_six_moves(),
          command_results.succeeded.two_six_moves(),
          command_results.failed.two_six_moves());
 
      fmt::print(
          "(6,2) moves: {} proposed ({} accepted and {} rejected); candidate "
          "construction: {} attempted ({} succeeded and {} failed).\n",
          statistics.proposed.six_two_moves(),
          statistics.accepted.six_two_moves(),
          statistics.rejected.six_two_moves(),
          command_results.attempted.six_two_moves(),
          command_results.succeeded.six_two_moves(),
          command_results.failed.six_two_moves());
 
      fmt::print(
          "(4,4) moves: {} proposed ({} accepted and {} rejected); candidate "
          "construction: {} attempted ({} succeeded and {} failed).\n",
          statistics.proposed.four_four_moves(),
          statistics.accepted.four_four_moves(),
          statistics.rejected.four_four_moves(),
          command_results.attempted.four_four_moves(),
          command_results.succeeded.four_four_moves(),
          command_results.failed.four_four_moves());
    }
 
   public:
    [[nodiscard]] auto attempt_transition(ManifoldType&                current,
                                          move_tracker::MoveType const move,
                                          long double const trial_value)
        -> ergodic_moves::MetropolisTransition
    {
      return {move, resolve_transition(current, m_command_results,
                                       m_run_statistics, move, trial_value)};
    }

    [[nodiscard]] auto attempt_transition(ManifoldType& current)
        -> ergodic_moves::MetropolisTransition
    { return sample_transition(current, m_command_results, m_run_statistics); }

    void initialize(ManifoldType const& manifold)
    { m_run_statistics.geometry = manifold.geometry(); }

    [[nodiscard]] auto operator()(ManifoldType const& t_manifold)
        -> ManifoldType
    {
#ifndef NDEBUG
      spdlog::debug("{} called.\n", CDT_PRETTY_FUNCTION);
#endif
 
      auto initial_statistics =
          m_resume_pending ? m_run_statistics : RunStatistics{};
      initial_statistics.geometry = t_manifold.geometry();
      auto initial_command_results =
          m_resume_pending ? m_command_results : CommandResults{};
      auto const initial_completed_passes =
          m_resume_pending ? m_completed_passes : Int_precision{};
      auto result = detail::execute_move_run(
          t_manifold, std::move(initial_command_results),
          std::move(initial_statistics), initial_completed_passes, m_cadence,
          detail::MoveRunIdentity{.algorithm = "Metropolis-Hastings",
                                  .seed      = seed(),
                                  .stream    = stream()},
          m_write_files,
          [this](ManifoldType current, RunStatistics statistics,
                 Int_precision const attempts) {
            return execute_pass(std::move(current), std::move(statistics),
                                attempts);
          },
          [](ManifoldType const&, CommandResults const& command_results,
             RunStatistics const& statistics) {
            print_results(command_results, statistics);
          },
          [this](ManifoldType const&   current,
                 CommandResults const& command_results,
                 RunStatistics const&  statistics,
                 Int_precision const   pass_number) {
            utilities::write_file(
                current, make_reproducibility_metadata(
                             current, utilities::ArtifactKind::CHECKPOINT,
                             pass_number, command_results, statistics));
          });
 
      m_command_results   = std::move(result.command_results);
      m_run_statistics    = std::move(result.strategy_state);
      m_checkpoint_events = result.checkpoint_events;
      m_completed_passes  = initial_completed_passes + m_cadence.passes();
      m_resume_pending    = false;
      return std::move(result.manifold);
    }

    void print_results() const
    { print_results(m_command_results, m_run_statistics); }
  };  // Metropolis

  using Metropolis_3 =
      MoveStrategy<MoveStrategyKind::METROPOLIS, manifolds::Manifold_3>;
}  // namespace cdt
 
#endif  // INCLUDE_METROPOLIS_HPP_