"""Reusable search templates for AtCoder Heuristic Contest (AHC).

Python 3.10+ / PyPy compatible. The module has no third-party dependencies.

Beam search contract:
    expand(state) -> Iterable[tuple[action, next_state]]
    evaluate(state) -> numeric score (larger is better by default)

Simulated annealing contract (in-place variant):
    mutate(state, rng) -> move_token
        Mutates `state` in place and returns enough information to undo the move.
        Return None to indicate that no move was made.
    undo(state, move_token) -> None
        Restores the state after a rejected move.
    evaluate(state) -> numeric score
    candidate_score(state, move_token, current_score) -> numeric score (optional)
        Called after mutation; use it for O(1)/incremental score updates.

Both algorithms return their best solution found, not necessarily their last/current state.
"""

from __future__ import annotations

import copy
import math
import random
import time
from dataclasses import dataclass
from typing import Callable, Generic, Hashable, Iterable, Optional, TypeVar

S = TypeVar("S")  # State
A = TypeVar("A")  # Action
M = TypeVar("M")  # Reversible move token
K = TypeVar("K", bound=Hashable)


class Timer:
    """Small monotonic timer shared by the search routines."""

    __slots__ = ("_started_at", "time_limit")

    def __init__(self, time_limit: Optional[float] = None) -> None:
        if time_limit is not None and time_limit < 0:
            raise ValueError("time_limit must be >= 0 or None")
        self._started_at = time.perf_counter()
        self.time_limit = time_limit

    @property
    def elapsed(self) -> float:
        return time.perf_counter() - self._started_at

    @property
    def remaining(self) -> Optional[float]:
        if self.time_limit is None:
            return None
        return max(0.0, self.time_limit - self.elapsed)

    @property
    def expired(self) -> bool:
        return self.time_limit is not None and self.elapsed >= self.time_limit


def _is_better(a: float, b: float, maximize: bool) -> bool:
    """Return whether score a is strictly better than score b; NaN is worst."""
    a_nan = math.isnan(a)
    b_nan = math.isnan(b)
    if a_nan:
        return False
    if b_nan:
        return True
    return a > b if maximize else a < b


def _sort_key(score: float, maximize: bool) -> tuple[bool, float]:
    # NaNs always sort last. The second key is ascending for both objectives.
    return (math.isnan(score), -score if maximize else score)


@dataclass(slots=True)
class _BeamNode(Generic[S, A]):
    state: S
    score: float
    depth: int
    parent: Optional["_BeamNode[S, A]"]
    action: Optional[A]
    terminal: bool


@dataclass(slots=True)
class BeamResult(Generic[S, A]):
    """Result from :func:`beam_search`.

    `state` is the best terminal state found, if any. Otherwise, it is the
    highest-scoring state in the deepest beam reached. `actions` is None when
    keep_path=False; otherwise it is the action sequence from the root.
    """

    state: S
    score: float
    depth: int
    terminal: bool
    actions: Optional[list[A]]
    elapsed: float
    expanded: int
    generated: int
    reached_depth: int
    timed_out: bool


def beam_search(
    initial_state: S,
    expand: Callable[[S], Iterable[tuple[A, S]]],
    evaluate: Callable[[S], float],
    *,
    width: int = 100,
    max_depth: int = 100,
    is_terminal: Optional[Callable[[S], bool]] = None,
    state_key: Optional[Callable[[S], K]] = None,
    maximize: bool = True,
    time_limit: Optional[float] = None,
    keep_path: bool = True,
    time_check_interval: int = 64,
) -> BeamResult[S, A]:
    """Run a generic beam search.

    Args:
        initial_state: Root state. It is not copied; callers should not mutate it
            from outside while the search is running.
        expand: Given a state, yield `(action, next_state)` pairs. `next_state`
            must be independently usable and must not be a mutable alias of its
            parent or of another candidate unless that sharing is intentional.
        evaluate: Numeric score. Higher is better when maximize=True.
        width: Maximum number of non-terminal states retained at each depth.
        max_depth: Maximum number of transitions from the root.
        is_terminal: Optional predicate. Terminal states are saved separately
            and are never expanded further.
        state_key: Optional canonical key for deduplicating candidates *within
            the same depth*. States with the same key compete by score; the
            better one is retained. The key must be hashable.
        maximize: Set False for minimization.
        time_limit: Optional wall-clock budget in seconds.
        keep_path: Track parent links and return the action sequence. Disable it
            to reduce memory and per-candidate overhead.
        time_check_interval: Check the clock every this many generated children
            in addition to checking once per expanded state. Smaller is tighter
            on time but incurs slightly more overhead.

    Notes:
        * If any terminal state is found, the best terminal state is returned,
          even if a non-terminal state has a numerically better score.
        * If no terminal state is found, returns the best state in the deepest
          non-empty beam reached (not the best score across all depths).
        * If `max_depth == 0`, returns the initial state.
        * Ties are stable: earlier generated candidates win.
    """
    if width <= 0:
        raise ValueError("width must be >= 1")
    if max_depth < 0:
        raise ValueError("max_depth must be >= 0")
    if time_check_interval <= 0:
        raise ValueError("time_check_interval must be >= 1")

    timer = Timer(time_limit)
    terminal_test = is_terminal or (lambda _state: False)
    root_terminal = bool(terminal_test(initial_state))
    root = _BeamNode(
        state=initial_state,
        score=float(evaluate(initial_state)),
        depth=0,
        parent=None,
        action=None,
        terminal=root_terminal,
    )
    beam: list[_BeamNode[S, A]] = [root]
    best_terminal: Optional[_BeamNode[S, A]] = root if root_terminal else None
    expanded = 0
    generated = 0
    reached_depth = 0
    timed_out = timer.expired

    def make_result(node: _BeamNode[S, A]) -> BeamResult[S, A]:
        actions: Optional[list[A]]
        if not keep_path:
            actions = None
        else:
            rev: list[A] = []
            cur: Optional[_BeamNode[S, A]] = node
            while cur is not None and cur.parent is not None:
                # Every non-root node has exactly one action, possibly None if
                # the caller deliberately uses None as an action value.
                rev.append(cur.action)  # type: ignore[arg-type]
                cur = cur.parent
            rev.reverse()
            actions = rev
        return BeamResult(
            state=node.state,
            score=node.score,
            depth=node.depth,
            terminal=node.terminal,
            actions=actions,
            elapsed=timer.elapsed,
            expanded=expanded,
            generated=generated,
            reached_depth=reached_depth,
            timed_out=timed_out,
        )

    if root_terminal or max_depth == 0 or timed_out:
        return make_result(root)

    for depth in range(max_depth):
        candidates: list[_BeamNode[S, A]] = []
        seen: dict[K, int] = {} if state_key is not None else {}  # type: ignore[assignment]
        stop_generating = False

        for parent in beam:
            if timer.expired:
                timed_out = True
                stop_generating = True
                break
            expanded += 1
            for action, child_state in expand(parent.state):
                generated += 1
                child_score = float(evaluate(child_state))
                child_terminal = bool(terminal_test(child_state))
                child = _BeamNode(
                    state=child_state,
                    score=child_score,
                    depth=depth + 1,
                    parent=parent if keep_path else None,
                    action=action if keep_path else None,
                    terminal=child_terminal,
                )

                if state_key is None:
                    candidates.append(child)
                else:
                    key = state_key(child_state)
                    old_index = seen.get(key)
                    if old_index is None:
                        seen[key] = len(candidates)
                        candidates.append(child)
                    elif _is_better(child.score, candidates[old_index].score, maximize):
                        candidates[old_index] = child

                if generated % time_check_interval == 0 and timer.expired:
                    timed_out = True
                    stop_generating = True
                    break
            if stop_generating:
                break

        if not candidates:
            break

        # Keep the best terminal among candidates, across every depth.
        for child in candidates:
            if child.terminal and (
                best_terminal is None
                or _is_better(child.score, best_terminal.score, maximize)
            ):
                best_terminal = child

        # Terminal states are archived separately and not expanded further.
        non_terminal = [node for node in candidates if not node.terminal]
        non_terminal.sort(key=lambda node: _sort_key(node.score, maximize))
        next_beam = non_terminal[:width]
        if not next_beam:
            # We may still have found a terminal solution in this layer.
            break

        beam = next_beam
        reached_depth = depth + 1
        if stop_generating:
            break

    if best_terminal is not None:
        return make_result(best_terminal)
    # `beam` is the deepest non-empty beam reached. It contains one depth only.
    best_node = min(beam, key=lambda node: _sort_key(node.score, maximize))
    return make_result(best_node)


@dataclass(slots=True)
class AnnealingResult(Generic[S]):
    """Result from :func:`simulated_annealing`.

    `best_state` is an independent snapshot made with copy_state whenever a new
    best is discovered. `current_state` is the final mutable search state.
    """

    best_state: S
    best_score: float
    current_state: S
    current_score: float
    iterations: int
    accepted: int
    improved: int
    best_iteration: int
    elapsed: float
    timed_out: bool


def simulated_annealing(
    initial_state: S,
    mutate: Callable[[S, random.Random], Optional[M]],
    undo: Callable[[S, M], None],
    evaluate: Callable[[S], float],
    *,
    time_limit: Optional[float] = None,
    max_iterations: Optional[int] = None,
    t_start: float = 1.0,
    t_end: float = 1e-3,
    temperature_schedule: Optional[Callable[[float], float]] = None,
    maximize: bool = True,
    candidate_score: Optional[Callable[[S, M, float], float]] = None,
    copy_state: Callable[[S], S] = copy.deepcopy,
    seed: Optional[int] = None,
    time_check_interval: int = 32,
) -> AnnealingResult[S]:
    """Run simulated annealing using reversible in-place mutations.

    Args:
        initial_state: Initial solution. The algorithm first clones it using
            `copy_state` so the caller's object is not mutated.
        mutate: Mutate `state` in place and return an undo token. Return None
            only when no mutation was performed.
        undo: Undo a rejected mutation using its token. It must restore *all*
            state data modified by `mutate`.
        evaluate: Full score computation, used initially and as the default
            candidate evaluator.
        time_limit / max_iterations: At least one stopping limit is required;
            when both are supplied, whichever is reached first stops the search.
        t_start / t_end: Positive endpoints for the default geometric schedule.
        temperature_schedule: Optional function mapping progress in [0, 1] to
            temperature. Progress is based on elapsed time when time_limit is
            supplied, otherwise on iteration count.
        candidate_score: Optional incremental evaluator called after mutation
            as `candidate_score(state, move, current_score)`. When absent,
            `evaluate(state)` is used. Ensure its result matches full evaluation.
        copy_state: State-copy function. Use `list.copy`, a custom copy, etc.
            when deepcopy is too expensive. It is used only initially and when
            a new global best is found, not on every iteration.
        seed: Seed for reproducible random moves and acceptance decisions.
        time_check_interval: Time/temperature update frequency. Lower values
            enforce the deadline more tightly but add overhead.

    The move token must be non-None for a real mutation because None denotes a
    no-op. Scores can be maximized or minimized with `maximize`.
    """
    if time_limit is None and max_iterations is None:
        raise ValueError("specify time_limit, max_iterations, or both")
    if time_limit is not None and time_limit < 0:
        raise ValueError("time_limit must be >= 0")
    if max_iterations is not None and max_iterations < 0:
        raise ValueError("max_iterations must be >= 0")
    if t_start <= 0 or t_end <= 0:
        raise ValueError("t_start and t_end must be > 0")
    if time_check_interval <= 0:
        raise ValueError("time_check_interval must be >= 1")

    rng = random.Random(seed)
    timer = Timer(time_limit)
    state = copy_state(initial_state)
    current_score = float(evaluate(state))
    best_state = copy_state(state)
    best_score = current_score
    best_iteration = 0
    iterations = 0
    accepted = 0
    improved = 0
    timed_out = timer.expired
    temperature = t_start

    if max_iterations == 0 or timed_out:
        return AnnealingResult(
            best_state, best_score, state, current_score, iterations, accepted,
            improved, best_iteration, timer.elapsed, timed_out
        )

    while max_iterations is None or iterations < max_iterations:
        # Update time, progress, and temperature periodically to reduce overhead.
        if iterations % time_check_interval == 0:
            elapsed = timer.elapsed
            if time_limit is not None and elapsed >= time_limit:
                timed_out = True
                break
            if temperature_schedule is not None:
                if time_limit is not None:
                    progress = min(1.0, elapsed / time_limit) if time_limit > 0 else 1.0
                else:
                    assert max_iterations is not None
                    progress = min(1.0, iterations / max_iterations) if max_iterations > 0 else 1.0
                temperature = float(temperature_schedule(progress))
            else:
                if time_limit is not None:
                    progress = min(1.0, elapsed / time_limit) if time_limit > 0 else 1.0
                else:
                    assert max_iterations is not None
                    progress = min(1.0, iterations / max_iterations) if max_iterations > 0 else 1.0
                temperature = t_start * ((t_end / t_start) ** progress)

        move = mutate(state, rng)
        if move is not None:
            if candidate_score is None:
                new_score = float(evaluate(state))
            else:
                new_score = float(candidate_score(state, move, current_score))

            delta = new_score - current_score if maximize else current_score - new_score
            take = delta >= 0.0
            if not take and temperature > 0.0:
                # delta is negative here, so exp(delta / T) is in (0, 1).
                take = rng.random() < math.exp(delta / temperature)

            if take:
                current_score = new_score
                accepted += 1
                if _is_better(current_score, best_score, maximize):
                    best_state = copy_state(state)
                    best_score = current_score
                    best_iteration = iterations + 1
                    improved += 1
            else:
                undo(state, move)
        iterations += 1

    elapsed = timer.elapsed
    if time_limit is not None and elapsed >= time_limit:
        timed_out = True
    return AnnealingResult(
        best_state=best_state,
        best_score=best_score,
        current_state=state,
        current_score=current_score,
        iterations=iterations,
        accepted=accepted,
        improved=improved,
        best_iteration=best_iteration,
        elapsed=elapsed,
        timed_out=timed_out,
    )


__all__ = [
    "Timer",
    "BeamResult",
    "beam_search",
    "AnnealingResult",
    "simulated_annealing",
]

#pragma once

// Generic AHC search templates (C++17, header-only).
// Compile with -O2 or -O3. No external dependencies.
//
// Beam search callback contract:
//   expand(const State&, vector<pair<Action, State>>& out)
//   evaluate(const State&) -> numeric score
//   is_terminal(const State&) -> bool
//   state_key(const State&) -> hashable canonical key
//
// Simulated annealing callback contract:
//   mutate(State&, rng) -> Move   // mutates in-place and returns undo token
//   undo(State&, const Move&)     // fully restores rejected moves
//   evaluate(const State&) -> numeric score
//   candidate_score(const State&, const Move&, double current_score) -> score
//       Use evaluate(state) for the simple version, or implement delta scoring.

#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <optional>
#include <random>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>

namespace ahc {

class Timer {
 public:
  explicit Timer(double time_limit_sec = std::numeric_limits<double>::infinity())
      : start_(Clock::now()), limit_(time_limit_sec) {}

  double elapsed() const {
    return std::chrono::duration<double>(Clock::now() - start_).count();
  }

  double remaining() const {
    if (!std::isfinite(limit_)) return limit_;
    return std::max(0.0, limit_ - elapsed());
  }

  bool expired() const {
    return std::isfinite(limit_) && elapsed() >= limit_;
  }

 private:
  using Clock = std::chrono::steady_clock;
  Clock::time_point start_;
  double limit_;
};

struct BeamConfig {
  std::size_t width = 100;
  int max_depth = 100;
  double time_limit_sec = std::numeric_limits<double>::infinity();
  bool maximize = true;
  bool deduplicate = true;
  bool keep_path = true;
  std::size_t time_check_interval = 64;
};

template <class State, class Action>
struct BeamResult {
  State state;
  double score;
  int depth;
  bool terminal;
  std::vector<Action> actions;
  bool path_available;
  double elapsed_sec;
  std::size_t expanded;
  std::size_t generated;
  int reached_depth;
  bool timed_out;
};

namespace detail {
inline bool better(double a, double b, bool maximize) {
  const bool a_nan = std::isnan(a);
  const bool b_nan = std::isnan(b);
  if (a_nan) return false;
  if (b_nan) return true;
  return maximize ? (a > b) : (a < b);
}
}  // namespace detail

// state_key must return a key supported by std::hash and operator==.
// To disable deduplication at runtime, set config.deduplicate=false (the key
// type still needs to be hashable because it is part of this template's type).
template <class State, class Action, class Expand, class Evaluate,
          class IsTerminal, class StateKey>
BeamResult<State, Action> beam_search(
    State initial_state, Expand expand, Evaluate evaluate,
    IsTerminal is_terminal, StateKey state_key, BeamConfig config = {}) {
  using Key = std::decay_t<decltype(state_key(std::declval<const State&>()))>;

  if (config.width == 0) throw std::invalid_argument("BeamConfig.width must be >= 1");
  if (config.max_depth < 0) throw std::invalid_argument("BeamConfig.max_depth must be >= 0");
  if (config.time_limit_sec < 0.0) throw std::invalid_argument("time_limit_sec must be >= 0");
  if (config.time_check_interval == 0) config.time_check_interval = 1;

  struct Node {
    State state;
    double score;
    int depth;
    std::size_t parent;
    std::optional<Action> action;
    bool terminal;
  };
  struct Candidate {
    State state;
    double score;
    int depth;
    std::size_t parent;
    std::optional<Action> action;
    bool terminal;
  };

  Timer timer(config.time_limit_sec);
  std::vector<Node> arena;
  // Avoid an enormous up-front allocation when width/depth are large but the
  // time limit is short. The arena grows as selected nodes are retained.
  arena.reserve(std::min<std::size_t>(config.width, 4096) *
                std::min<std::size_t>(static_cast<std::size_t>(config.max_depth) + 1, 16) + 1);
  arena.push_back(Node{std::move(initial_state), 0.0, 0, 0, std::nullopt, false});
  arena[0].score = static_cast<double>(evaluate(arena[0].state));
  arena[0].terminal = static_cast<bool>(is_terminal(arena[0].state));

  auto extract_actions = [&](std::size_t node_id) {
    std::vector<Action> reversed;
    if (config.keep_path) {
      while (arena[node_id].depth > 0) {
        if (arena[node_id].action.has_value()) {
          reversed.push_back(*arena[node_id].action);
        }
        node_id = arena[node_id].parent;
      }
      std::reverse(reversed.begin(), reversed.end());
    }
    return reversed;
  };

  if (arena[0].terminal) {
    return BeamResult<State, Action>{arena[0].state, arena[0].score, 0, true,
        {}, config.keep_path, timer.elapsed(), 0, 0, 0, timer.expired()};
  }

  std::vector<std::size_t> beam{0};
  std::optional<State> best_terminal_state;
  std::vector<Action> best_terminal_actions;
  double best_terminal_score = std::numeric_limits<double>::quiet_NaN();
  int best_terminal_depth = -1;
  bool timed_out = timer.expired();
  std::size_t expanded = 0;
  std::size_t generated = 0;
  int reached_depth = 0;

  for (int depth = 0; depth < config.max_depth && !timed_out; ++depth) {
    std::vector<Candidate> candidates;
    std::unordered_map<Key, std::size_t> seen;
    std::vector<std::pair<Action, State>> successors;
    bool stop_generating = false;

    for (const std::size_t parent_id : beam) {
      if (timer.expired()) {
        timed_out = true;
        stop_generating = true;
        break;
      }
      ++expanded;
      successors.clear();
      expand(arena[parent_id].state, successors);

      for (auto& successor : successors) {
        ++generated;
        const double score = static_cast<double>(evaluate(successor.second));
        const bool terminal = static_cast<bool>(is_terminal(successor.second));
        std::optional<Action> action;
        if (config.keep_path) action.emplace(std::move(successor.first));
        Candidate candidate{std::move(successor.second), score, depth + 1,
                            parent_id, std::move(action), terminal};

        if (!config.deduplicate) {
          candidates.push_back(std::move(candidate));
        } else {
          Key key = state_key(candidate.state);
          auto it = seen.find(key);
          if (it == seen.end()) {
            seen.emplace(std::move(key), candidates.size());
            candidates.push_back(std::move(candidate));
          } else if (detail::better(candidate.score, candidates[it->second].score,
                                    config.maximize)) {
            candidates[it->second] = std::move(candidate);
          }
        }

        if (generated % config.time_check_interval == 0 && timer.expired()) {
          timed_out = true;
          stop_generating = true;
          break;
        }
      }
      if (stop_generating) break;
    }

    if (candidates.empty()) break;

    std::vector<std::size_t> non_terminal_indices;
    non_terminal_indices.reserve(candidates.size());
    for (std::size_t i = 0; i < candidates.size(); ++i) {
      Candidate& c = candidates[i];
      if (c.terminal) {
        if (!best_terminal_state.has_value() ||
            detail::better(c.score, best_terminal_score, config.maximize)) {
          best_terminal_score = c.score;
          best_terminal_depth = c.depth;
          best_terminal_state = c.state;
          best_terminal_actions = extract_actions(c.parent);
          if (config.keep_path && c.action.has_value()) {
            best_terminal_actions.push_back(*c.action);
          }
        }
      } else {
        non_terminal_indices.push_back(i);
      }
    }

    std::stable_sort(non_terminal_indices.begin(), non_terminal_indices.end(),
        [&](std::size_t ia, std::size_t ib) {
          const double a = candidates[ia].score;
          const double b = candidates[ib].score;
          const bool a_nan = std::isnan(a);
          const bool b_nan = std::isnan(b);
          if (a_nan != b_nan) return !a_nan;
          if (a_nan) return false;
          return config.maximize ? (a > b) : (a < b);
        });

    const std::size_t take = std::min(config.width, non_terminal_indices.size());
    std::vector<std::size_t> next_beam;
    next_beam.reserve(take);
    for (std::size_t i = 0; i < take; ++i) {
      Candidate& c = candidates[non_terminal_indices[i]];
      const std::size_t new_id = arena.size();
      arena.push_back(Node{std::move(c.state), c.score, c.depth,
          config.keep_path ? c.parent : 0, std::move(c.action), false});
      next_beam.push_back(new_id);
    }

    if (next_beam.empty()) break;
    beam = std::move(next_beam);
    reached_depth = depth + 1;
    if (stop_generating) break;
  }

  if (best_terminal_state.has_value()) {
    return BeamResult<State, Action>{std::move(*best_terminal_state), best_terminal_score,
        best_terminal_depth, true, std::move(best_terminal_actions), config.keep_path,
        timer.elapsed(), expanded, generated, reached_depth, timed_out};
  }

  // The beam is ordered from best to worst after each successful layer.
  const Node& best = arena[beam.front()];
  return BeamResult<State, Action>{best.state, best.score, best.depth, false,
      extract_actions(beam.front()), config.keep_path, timer.elapsed(), expanded,
      generated, reached_depth, timed_out};
}

struct AnnealConfig {
  double time_limit_sec = 1.8;
  std::size_t max_iterations = std::numeric_limits<std::size_t>::max();
  double t_start = 1.0;
  double t_end = 1e-3;
  bool maximize = true;
  std::uint64_t seed = 0;
  std::size_t time_check_interval = 32;
};

template <class State>
struct AnnealResult {
  State best_state;
  double best_score;
  State current_state;
  double current_score;
  std::size_t iterations;
  std::size_t accepted;
  std::size_t improved;
  std::size_t best_iteration;
  double elapsed_sec;
  bool timed_out;
};

// `mutate` always performs one valid mutation and returns an undo token.
// `candidate_score(state, move, current_score)` can call evaluate(state), or
// update the score incrementally in O(1) / O(size of affected region).
template <class State, class Mutate, class Undo, class Evaluate, class CandidateScore>
AnnealResult<State> simulated_annealing(
    State initial_state, Mutate mutate, Undo undo, Evaluate evaluate,
    CandidateScore candidate_score, AnnealConfig config = {}) {
  if (config.time_limit_sec < 0.0) throw std::invalid_argument("time_limit_sec must be >= 0");
  if (!std::isfinite(config.time_limit_sec) &&
      config.max_iterations == std::numeric_limits<std::size_t>::max()) {
    throw std::invalid_argument("provide a finite time limit or max_iterations");
  }
  if (config.t_start <= 0.0 || config.t_end <= 0.0) {
    throw std::invalid_argument("t_start and t_end must be > 0");
  }
  if (config.time_check_interval == 0) config.time_check_interval = 1;

  Timer timer(config.time_limit_sec);
  std::mt19937_64 rng(config.seed);
  std::uniform_real_distribution<double> uniform01(0.0, 1.0);
  State current = std::move(initial_state);
  double current_score = static_cast<double>(evaluate(current));
  State best = current;
  double best_score = current_score;
  std::size_t best_iteration = 0;
  std::size_t iterations = 0;
  std::size_t accepted = 0;
  std::size_t improved = 0;
  bool timed_out = timer.expired();
  double temperature = config.t_start;

  while (!timed_out && iterations < config.max_iterations) {
    if (iterations % config.time_check_interval == 0) {
      const double elapsed = timer.elapsed();
      if (std::isfinite(config.time_limit_sec) && elapsed >= config.time_limit_sec) {
        timed_out = true;
        break;
      }
      double progress;
      if (std::isfinite(config.time_limit_sec)) {
        progress = config.time_limit_sec > 0.0
            ? std::min(1.0, elapsed / config.time_limit_sec) : 1.0;
      } else if (config.max_iterations != std::numeric_limits<std::size_t>::max() &&
                 config.max_iterations > 0) {
        progress = std::min(1.0, static_cast<double>(iterations) /
                                      static_cast<double>(config.max_iterations));
      } else {
        progress = 0.0;
      }
      temperature = config.t_start *
          std::pow(config.t_end / config.t_start, progress);
    }

    auto move = mutate(current, rng);
    const double new_score = static_cast<double>(candidate_score(current, move, current_score));
    const double delta = config.maximize ? (new_score - current_score)
                                         : (current_score - new_score);
    bool take = delta >= 0.0;
    if (!take && temperature > 0.0) {
      take = uniform01(rng) < std::exp(delta / temperature);
    }

    if (take) {
      current_score = new_score;
      ++accepted;
      if (detail::better(current_score, best_score, config.maximize)) {
        best = current;
        best_score = current_score;
        best_iteration = iterations + 1;
        ++improved;
      }
    } else {
      undo(current, move);
    }
    ++iterations;
  }

  const double elapsed = timer.elapsed();
  if (std::isfinite(config.time_limit_sec) && elapsed >= config.time_limit_sec) {
    timed_out = true;
  }
  return AnnealResult<State>{std::move(best), best_score, std::move(current), current_score,
      iterations, accepted, improved, best_iteration, elapsed, timed_out};
}

}  // namespace ahc