# Hyper Heuristics

> When the user wants to build a hyper-heuristic — a search method that selects or generates low-level heuristics instead of searching solutions directly — covering selection hyper-heuristics (heuristic selection plus move acceptance), low-level heuristic pool design, learning and reward schemes, and generation hyper-heuristics. Also use when the user mentions "hyper-heuristic," "operator selection," "low-level heuristics," "heuristic selection," "move acceptance," "choice function," or when many candidate operators exist and none dominates across instances. For adaptive destroy/repair weights, see large-neighborhood-search; for choosing the overall search paradigm, see metaheuristic-design-principles.

- Skill: `hajibabaie/hyper-heuristics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add hajibabaie/hyper-heuristics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/hajibabaie/hyper-heuristics/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: hajibabaie (https://skillmd.com/u/hajibabaie)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/hajibabaie/hyper-heuristics

---


# Hyper-Heuristics

You are an expert in hyper-heuristics for combinatorial optimization — search methods that operate on a *space of heuristics* rather than directly on the space of solutions. This skill covers selection hyper-heuristics (heuristic selection plus move acceptance), the design of low-level heuristic pools, online learning and credit assignment (reward schemes), and generation hyper-heuristics that assemble new heuristics from components. Use the framework below to take a user from "I have many candidate operators and no idea which to apply when" to a reproducible hyper-heuristic with measured operator usage, a calibrated acceptance criterion, and a defensible ablation against simpler baselines.

## Initial Assessment

Establish these facts before writing any hyper-heuristic code:

- **Why a hyper-heuristic at all.** If one well-understood metaheuristic with one strong neighborhood already works, a hyper-heuristic adds machinery without value. The case for a hyper-heuristic is: a pool of plausible operators with instance-dependent usefulness, heterogeneous instances, or a requirement for cross-domain reuse.
- **Pool inventory.** Which low-level heuristics already exist (moves, repair rules, construction rules)? A hyper-heuristic cannot fix a weak pool — it only arbitrates among the heuristics it is given. Aim for 4–12 heuristics with genuinely different behaviors.
- **Constructive or perturbative low-level heuristics.** Constructive heuristics extend a partial solution (pick the next packing rule, the next dispatching rule); perturbative heuristics modify a complete solution (move, swap, ruin-and-recreate). The loop structure differs; decide first.
- **Objective and scaling.** One scalar objective the engine can read through the domain barrier. If hard constraints are penalized, fix the penalty weights before any credit learning — rewards inherit the objective's scale.
- **Evaluation cost and budget.** Iterations available = time budget / (heuristic call + evaluation). Credit assignment needs hundreds of selections per heuristic to mean anything; with fewer than ~50 calls per heuristic, use uniform random selection.
- **Per-call cost asymmetry.** Does the pool mix microsecond moves with millisecond ruin-and-recreate heuristics? If yes, credit must be improvement *per unit time*, not per call, or cheap heuristics will be unfairly favored — and vice versa.
- **Online vs offline learning.** Online: learn during the run on this instance (selection rules below). Offline: tune selection/acceptance parameters or evolve heuristics on a training instance set beforehand. Most practical systems combine both.
- **No-op behavior.** Can a heuristic return the solution unchanged (e.g., no feasible swap found)? Decide how no-ops are credited (zero reward) and detected, or they silently distort the statistics.
- **Acceptance scale.** Move acceptance needs either a temperature (Metropolis) or a history length (late acceptance). Both must be set relative to the objective scale and the iteration budget.
- **Single-domain or cross-domain.** A one-problem project can let problem knowledge leak into selection. A cross-domain tool must keep the domain barrier strict: the engine sees only objective values and heuristic indices.
- **Baselines.** Always run (a) uniform random selection with the same pool and acceptance, and (b) the single best heuristic alone. The learning layer must beat both to justify itself.
- **Reproducibility.** One `np.random.default_rng(seed)` per run, seeds and parameters logged per run, usage statistics saved with results.

## Hyper-Heuristic Anatomy

### Heuristic space and the domain barrier

A hyper-heuristic searches over heuristics: "heuristics to choose heuristics" (Cowling, Kendall & Soubeiga 2001, "A Hyperheuristic Approach to Scheduling a Sales Summit"; the term itself appears in Denzinger, Fuchs & Fuchs 1997). The defining architectural idea is the **domain barrier**: the high-level strategy sees only (i) the indices of the low-level heuristics, (ii) the objective value returned after applying one, and (iii) bookkeeping such as elapsed time. It never sees the solution representation. Everything problem-specific lives below the barrier, inside the low-level heuristics and the evaluation function. This is what makes the high-level strategy reusable across domains, and it is enforced literally in the HyFlex benchmark framework (Ochoa et al. 2012, "HyFlex: A Benchmark Framework for Cross-domain Heuristic Search").

The standard classification (Burke et al. 2010, "A Classification of Hyper-heuristic Approaches"; surveyed in Burke et al. 2013, "Hyper-heuristics: A Survey of the State of the Art" and Drake et al. 2020, "Recent Advances in Selection Hyper-heuristics") has two axes:

| Axis | Options | Coverage in this skill |
|---|---|---|
| Nature of the heuristic space | **selection** (choose from a fixed pool) vs **generation** (build new heuristics from components) | selection in depth; generation as a working sketch |
| Nature of the low-level heuristics | **constructive** (extend partial solutions) vs **perturbative** (modify complete solutions) | perturbative engine + constructive bin-packing rules |
| Feedback | online learning, offline learning, no learning | online credit schemes; offline tuning/evolution |

### Selection hyper-heuristic = heuristic selection + move acceptance

A selection hyper-heuristic decomposes into two nearly independent components (Özcan, Bilgin & Korkmaz 2008, "A Comprehensive Analysis of Hyper-heuristics"): a **selection rule** that picks the next low-level heuristic, and a **move acceptance** rule that decides whether the modified solution replaces the incumbent. Both matter; in the CHeSC 2011 cross-domain competition, the winner AdapHH (Mısır et al. 2012) combined adaptive selection with an adaptive acceptance threshold, and ablations show acceptance often contributes more than selection.

**Credit assignment.** After applying heuristic $i$ to incumbent $x$ and obtaining $x'$, assign a reward and update a quality estimate $q_i$ by a recency-weighted average:

$$
r_t = \frac{\max\big(0,\; f(x) - f(x')\big)}{s_t} \;+\; \mathbf{1}\{f(x') < f(x^*)\}, \qquad
q_i \leftarrow (1-\lambda)\, q_i + \lambda\, r_t ,
$$

where $s_t$ is a running scale of $|f|$ (keeps rewards dimensionless) and the indicator adds a new-global-best bonus. This is non-stationary reinforcement learning in the sense of Nareyek (2003), "Choosing Search Heuristics by Non-Stationary Reinforcement Learning": the value of a heuristic changes as the search moves through the landscape, so old evidence must decay.

**Selection rules**, in increasing sophistication:

- *Uniform random* — the mandatory baseline; "simple random" in the HH literature.
- *$\varepsilon$-greedy* — exploit $\arg\max_i q_i$ with probability $1-\varepsilon$, explore uniformly otherwise.
- *Roulette (probability matching)* — select $i$ with probability $\propto q_i + q_{\min}$; smooth but slow to commit.
- *Upper confidence bound (UCB)* — treat selection as a multi-armed bandit (Fialho et al. 2010, "Analyzing Bandit-based Adaptive Operator Selection Mechanisms"):

$$
i_t = \arg\max_i \left( q_i + c \sqrt{\frac{\ln t}{n_i}} \right),
$$

with $n_i$ the usage count and $c$ the exploration strength. The decay $\lambda$ in $q_i$ supplies the non-stationarity handling that plain UCB lacks.

- *Choice function* — a learned score combining solo performance, pairwise (follow-up) performance, and time since last use (Cowling et al. 2001); implementation in Advanced Techniques.

**Move acceptance** rules, with their characters:

- *Improve-only* (`improve_only`): accept iff $f(x') \le f(x)$. Greedy; stalls in local optima but a clean ablation baseline.
- *Accept-all* (`accept_all`): the search becomes a random walk guided only by selection; useful with very strong, self-improving heuristics (e.g., a pool of ruin-and-recreate operators).
- *Metropolis*: accept worse moves with probability $e^{-(f(x')-f(x))/T}$, $T$ cooled geometrically — simulated-annealing acceptance.
- *Late acceptance* (LAHC; Burke & Bykov 2017, "The Late Acceptance Hill-Climbing Heuristic"): accept iff $f(x') \le f(x)$ **or** $f(x') \le f_{t-L}$, the incumbent value $L$ iterations ago, kept in a ring buffer. One parameter ($L$), no temperature calibration, strong empirical record — the recommended default.

### Decision guidance

- **Use a selection hyper-heuristic when** you have ≥4 heuristics whose relative usefulness varies by instance or by search phase, and an evaluation budget of at least a few thousand iterations.
- **Use the ALNS framing instead when** your heuristics decompose naturally into destroy and repair halves applied as a pair — that is a selection hyper-heuristic specialized to routing/scheduling; see large-neighborhood-search for segment-based weight updates.
- **Use a single tuned metaheuristic when** one neighborhood dominates (e.g., 2-opt on plain TSP) — see metaheuristic-design-principles for the selection logic.
- **Use a generation hyper-heuristic when** the deliverable is a *reusable fast rule* (a dispatching rule, a packing rule) to be applied online to many future instances, not one solution to one instance.
- **Complexity.** The high-level layer costs $O(k)$ per iteration for $k$ heuristics — negligible. Runtime is dominated by heuristic calls and objective evaluations; spend engineering effort below the barrier (delta evaluation inside low-level heuristics), not above it.

### Parameter guidance

| Parameter | Typical range | Increasing it buys | At the cost of |
|---|---|---|---|
| Pool size $k$ | 4–12 | Behavioral coverage | Credit split thinner, slower learning |
| $\varepsilon$ ($\varepsilon$-greedy) | 0.05–0.30 | Starvation protection, exploration | Budget spent on weak heuristics |
| UCB exploration $c$ | 0.3–2.0 | Re-testing of currently weak heuristics | Slower commitment to the best heuristic |
| Reward decay (memory) $1-\lambda$ | 0.90–0.99 | Stable estimates, long memory | Slow reaction to search-phase changes |
| New-best bonus | 0.5–2.0 (normalized units) | Bias toward record-breaking heuristics | Overrates single lucky events |
| LAHC list length $L$ | 50–5,000 | Diversification, better final quality | Slower convergence; scale $L$ with budget |
| Metropolis $T_0$ | 1–10% of $\|f(x_0)\|$ | Early exploration | Random-walk burn-in |
| Cooling rate | set so $T_{\text{end}}/T_0 \approx 10^{-3}$ over the budget | Final intensification | Long frozen tail |
| Segment length (batched updates) | 50–200 iterations | Stable per-segment scores | Slower adaptation |

Tune acceptance ($L$ or the temperature schedule) before tuning selection; it usually has the larger effect. For offline tuning of these parameters over an instance set, use the protocol in optuna-hyperparameter-tuning. For the design of the low-level moves themselves, reuse the operator catalogs in mutation-and-perturbation-operators and local-search-and-neighborhoods rather than re-inventing them.

## Reusable Selection Hyper-Heuristic Engine

The engine is problem-independent and enforces the domain barrier: it sees the solution only as an opaque object passed to `copy_solution`, the heuristics, and `objective`. Low-level heuristics have the uniform signature `(solution, rng) -> solution` and may mutate their argument (the engine always passes a copy).

```text
SELECTION-HH(x0, pool H = {h1..hk}, budget M)
  x <- x0;  x* <- x
  q_i <- 0, n_i <- 0  for all i                  // learned quality, usage count
  for t = 1 .. M:
      i  <- SELECT(q, n, t)                      // uniform | eps-greedy | roulette | UCB
      x' <- h_i(copy(x))                         // apply low-level heuristic
      r  <- max(0, f(x) - f(x')) / scale         // credit ...
             + bonus if f(x') < f(x*)            //   ... with new-best bonus
      q_i <- (1 - lambda) * q_i + lambda * r     // recency-weighted average
      n_i <- n_i + 1
      if ACCEPT(f(x'), f(x), memory):            // improve-only | accept-all | LAHC | Metropolis
          x <- x'
      if f(x') < f(x*):  x* <- x'
      update acceptance memory                   // LAHC ring buffer / temperature
  return x*, usage statistics
```

```python
"""Problem-independent single-point selection hyper-heuristic (minimization)."""
from dataclasses import dataclass, field
from typing import Any, Callable

import numpy as np

LowLevelHeuristic = Callable[[Any, np.random.Generator], Any]


@dataclass
class HHResult:
    """Outcome and bookkeeping of one hyper-heuristic run."""
    best: Any
    f_best: float
    usage: np.ndarray       # times each heuristic was selected
    accepted: np.ndarray    # times its candidate was accepted
    new_best: np.ndarray    # times it produced a new global best
    trace: list[float] = field(default_factory=list)  # best-so-far per iteration


class SelectionHyperHeuristic:
    """Heuristic selection + move acceptance over a pool of low-level heuristics.

    selection : 'uniform' | 'epsilon_greedy' | 'roulette' | 'ucb'
    acceptance: 'improve_only' | 'accept_all' | 'late_acceptance' | 'metropolis'
    """

    def __init__(
        self,
        heuristics: list[LowLevelHeuristic],
        objective: Callable[[Any], float],
        copy_solution: Callable[[Any], Any],
        selection: str = "ucb",
        acceptance: str = "late_acceptance",
        epsilon: float = 0.15,
        ucb_c: float = 0.5,
        reward_decay: float = 0.95,
        la_length: int = 100,
        t0: float | None = None,
        cooling: float = 0.999,
        seed: int = 0,
    ) -> None:
        self.heuristics = heuristics
        self.objective = objective
        self.copy_solution = copy_solution
        self.selection = selection
        self.acceptance = acceptance
        self.epsilon = epsilon
        self.ucb_c = ucb_c
        self.reward_decay = reward_decay
        self.la_length = la_length
        self.t0 = t0
        self.cooling = cooling
        self.rng = np.random.default_rng(seed)

    def _select(self, q: np.ndarray, n: np.ndarray, it: int) -> int:
        """Pick a heuristic index from quality estimates q and usage counts n."""
        k = q.size
        if self.selection == "uniform":
            return int(self.rng.integers(k))
        if self.selection == "epsilon_greedy":
            if self.rng.random() < self.epsilon:
                return int(self.rng.integers(k))
            return int(self.rng.choice(np.flatnonzero(q == q.max())))
        if self.selection == "roulette":
            w = q - q.min() + 1e-6
            return int(self.rng.choice(k, p=w / w.sum()))
        untried = np.flatnonzero(n == 0)  # selection == 'ucb'
        if untried.size:
            return int(untried[0])
        bound = q + self.ucb_c * np.sqrt(np.log(it + 1.0) / n)
        return int(self.rng.choice(np.flatnonzero(bound == bound.max())))

    def _accept(self, f_new: float, f_cur: float, f_late: float, t: float) -> bool:
        """Move acceptance: does the candidate replace the incumbent?"""
        if self.acceptance == "improve_only":
            return f_new <= f_cur
        if self.acceptance == "accept_all":
            return True
        if self.acceptance == "late_acceptance":
            return f_new <= f_cur or f_new <= f_late
        if f_new <= f_cur:  # acceptance == 'metropolis'
            return True
        return self.rng.random() < np.exp(-(f_new - f_cur) / max(t, 1e-12))

    def run(self, x0: Any, n_iters: int) -> HHResult:
        """Run the selection hyper-heuristic loop from initial solution x0."""
        k = len(self.heuristics)
        cur = self.copy_solution(x0)
        f_cur = self.objective(cur)
        best, f_best = self.copy_solution(cur), f_cur
        q = np.zeros(k)                      # recency-weighted reward per heuristic
        n = np.zeros(k)
        accepted = np.zeros(k, dtype=int)
        new_best = np.zeros(k, dtype=int)
        late = np.full(self.la_length, f_cur)
        t = self.t0 if self.t0 is not None else 0.05 * abs(f_cur) + 1e-9
        scale = abs(f_cur) + 1e-9            # running reward normalizer
        trace: list[float] = []
        for it in range(n_iters):
            i = self._select(q, n, it)
            cand = self.heuristics[i](self.copy_solution(cur), self.rng)
            f_cand = self.objective(cand)
            n[i] += 1
            reward = max(0.0, f_cur - f_cand) / scale
            if f_cand < f_best:
                reward += 1.0                # new-global-best bonus
            q[i] = self.reward_decay * q[i] + (1.0 - self.reward_decay) * reward
            if self._accept(f_cand, f_cur, float(late[it % self.la_length]), t):
                cur, f_cur = cand, f_cand
                accepted[i] += 1
            if f_cand < f_best:
                best, f_best = self.copy_solution(cand), f_cand
                new_best[i] += 1
            late[it % self.la_length] = f_cur
            scale = 0.99 * scale + 0.01 * (abs(f_cur) + 1e-9)
            t *= self.cooling
            trace.append(f_best)
        return HHResult(best, f_best, n.astype(int), accepted, new_best, trace)
```

Smoke test on a synthetic problem: minimize the quadratic pseudo-boolean function $f(x) = x^\top Q x$, $x \in \{0,1\}^{60}$, with a four-heuristic pool that mixes blind perturbations and one greedy move. The greedy move uses the closed-form flip delta $\Delta_j = (1-2x_j)(Q_{jj} + 2\sum_{i\ne j} Q_{ij}x_i)$.

```python
"""Smoke test: binary quadratic minimization with a four-heuristic pool."""
import numpy as np

rng0 = np.random.default_rng(7)
N = 60
Q = rng0.normal(0.0, 1.0, (N, N))
Q = (Q + Q.T) / 2.0


def f_bqp(x: np.ndarray) -> float:
    """Quadratic pseudo-boolean objective x^T Q x (minimize)."""
    return float(x @ Q @ x)


def flip_one(x: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """Flip one random bit."""
    x[int(rng.integers(x.size))] ^= 1
    return x


def flip_three(x: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """Flip three distinct random bits."""
    x[rng.choice(x.size, 3, replace=False)] ^= 1
    return x


def flip_block(x: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """Flip a contiguous block of five bits (strong perturbation)."""
    j = int(rng.integers(x.size - 5))
    x[j:j + 5] ^= 1
    return x


def best_single_flip(x: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """Greedy: flip the bit with the steepest objective decrease."""
    s = Q @ x - np.diag(Q) * x
    deltas = (1 - 2 * x) * (np.diag(Q) + 2 * s)
    x[int(np.argmin(deltas))] ^= 1
    return x


pool = [flip_one, flip_three, flip_block, best_single_flip]
x0 = rng0.integers(0, 2, N)
hh = SelectionHyperHeuristic(pool, f_bqp, np.copy, selection="epsilon_greedy",
                             acceptance="late_acceptance", la_length=80, seed=11)
res = hh.run(x0, n_iters=3000)
print(f"f(x0)={f_bqp(x0):.1f}  f_best={res.f_best:.1f}")
names = ["flip1", "flip3", "block5", "greedy"]
print("usage:", dict(zip(names, res.usage)))
print("new_best:", dict(zip(names, res.new_best)))
# Expected: f_best around -229 from a start near -33; usage concentrates on
# the greedy flip (~1,900 of 3,000 calls, ~35 of ~36 new bests), with
# flip_one second (~870 calls) as the cheap escape move.
```

Read the usage statistics before celebrating: a healthy $\varepsilon$-greedy run shows the intensifier exploited heavily while the exploration floor keeps every heuristic alive. Selection rules differ in how the learning shows up: on this same instance UCB finds the same $f^{best}$ but with *near-uniform usage* (750 calls each), because once rewards vanish its confidence radii equalize the counts — under UCB, read the `new_best` and `accepted` columns, not raw usage. If a blind perturbation dominates `new_best`, the reward scale or the acceptance criterion is broken.

## Worked Example 1: Exam Timetabling Selection Hyper-Heuristic

Timetabling is the classic hyper-heuristic domain (the field grew out of timetabling and rostering systems; see timetabling-and-rostering for full domain models and benchmark formats). The instance here: $n$ events, $T$ timeslots, and a symmetric conflict matrix $C$ where $C_{ij}$ counts students enrolled in both events $i$ and $j$. The solution is a slot vector $s \in \{0,\dots,T-1\}^n$. The cost combines hard conflicts (same slot) and a proximity soft cost — a linear simplification of the Carter cost $2^{4-d}$ (Carter, Laporte & Lee 1996, "Examination Timetabling: Algorithmic Strategies and Applications"):

$$
f(s) = W \sum_{i<j,\; s_i = s_j} C_{ij} \;+\; \sum_{i<j,\; 1 \le |s_i - s_j| \le 2} C_{ij}\,\big(3 - |s_i - s_j|\big), \qquad W = 1000 .
$$

```python
"""Exam timetabling: instance, vectorized cost, and a five-heuristic pool."""
import numpy as np


def make_timetabling_instance(n_events: int, n_slots: int, density: float,
                              seed: int) -> tuple[np.ndarray, int]:
    """Random symmetric conflict matrix; C[i, j] = students shared by events i, j."""
    rng = np.random.default_rng(seed)
    raw = np.where(rng.random((n_events, n_events)) < density,
                   rng.integers(1, 10, (n_events, n_events)), 0)
    C = np.triu(raw, 1)
    return C + C.T, n_slots


def timetable_cost(slots: np.ndarray, C: np.ndarray,
                   hard_weight: float = 1000.0) -> float:
    """Hard: conflicting events in one slot. Soft: proximity max(0, 3 - d)."""
    D = np.abs(slots[:, None] - slots[None, :])
    hard = (C * (D == 0)).sum() / 2.0
    soft = (C * np.maximum(0, 3 - D) * (D > 0)).sum() / 2.0
    return hard_weight * float(hard) + float(soft)


def make_timetabling_pool(C: np.ndarray, n_slots: int) -> list:
    """Five perturbative low-level heuristics closing over the instance."""
    n = C.shape[0]

    def event_penalty(slots: np.ndarray) -> np.ndarray:
        """Per-event share of the total cost (for targeting bad events)."""
        D = np.abs(slots[:, None] - slots[None, :])
        pen = C * ((D == 0) * 1000.0 + np.maximum(0, 3 - D) * (D > 0))
        return pen.sum(axis=1)

    def move_random(slots: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        """Move one random event to one random slot."""
        slots[int(rng.integers(n))] = int(rng.integers(n_slots))
        return slots

    def move_worst_greedy(slots: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        """Pick among the 5 worst events; place it in its cheapest slot."""
        e = int(rng.choice(np.argsort(event_penalty(slots))[-5:]))
        costs = np.empty(n_slots)
        for sl in range(n_slots):
            d = np.abs(slots - sl)
            costs[sl] = (C[e] * ((d == 0) * 1000.0
                                 + np.maximum(0, 3 - d) * (d > 0))).sum()
        slots[e] = int(rng.choice(np.flatnonzero(costs == costs.min())))
        return slots

    def swap_events(slots: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        """Exchange the slots of two random events."""
        a, b = rng.choice(n, 2, replace=False)
        slots[a], slots[b] = slots[b], slots[a]
        return slots

    def kempe_chain(slots: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        """Swap a Kempe chain between a random event's slot and a target slot."""
        e = int(rng.integers(n))
        a, b = int(slots[e]), int(rng.integers(n_slots))
        if a == b:
            return slots
        in_ab = (slots == a) | (slots == b)
        chain, frontier = {e}, [e]
        while frontier:
            v = frontier.pop()
            for u in np.flatnonzero((C[v] > 0) & in_ab):
                if int(u) not in chain:
                    chain.add(int(u))
                    frontier.append(int(u))
        idx = np.fromiter(chain, dtype=int)
        slots[idx] = np.where(slots[idx] == a, b, a)
        return slots

    def shuffle_slot(slots: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        """Scatter all events of one random slot (diversification)."""
        members = np.flatnonzero(slots == int(rng.integers(n_slots)))
        if members.size:
            slots[members] = rng.integers(n_slots, size=members.size)
        return slots

    return [move_random, move_worst_greedy, swap_events, kempe_chain, shuffle_slot]
```

The pool deliberately mixes intensifiers (`move_worst_greedy`, `kempe_chain` — the standard timetabling/coloring move that swaps a connected component between two slots and never splits a conflict pair across the wrong side) with diversifiers (`move_random`, `shuffle_slot`). That mix is what the selection layer is for: the right blend changes between the "repair hard conflicts" phase and the "polish proximity cost" phase, and the credit scheme tracks the change.

```python
"""Run the engine on a 50-event, 9-slot instance and report usage."""
import numpy as np

C, n_slots = make_timetabling_instance(n_events=50, n_slots=9,
                                       density=0.20, seed=3)
pool = make_timetabling_pool(C, n_slots)
names = ["move_random", "move_worst_greedy", "swap_events",
         "kempe_chain", "shuffle_slot"]

x0 = np.random.default_rng(1).integers(n_slots, size=50)
hh = SelectionHyperHeuristic(
    heuristics=pool,
    objective=lambda s: timetable_cost(s, C),
    copy_solution=np.copy,
    selection="epsilon_greedy",
    acceptance="late_acceptance",
    la_length=200,
    seed=1,
)
res = hh.run(x0, n_iters=6000)
D = np.abs(res.best[:, None] - res.best[None, :])
hard_left = int((C * (D == 0)).sum() / 2)
print(f"start cost={timetable_cost(x0, C):.0f}  final cost={res.f_best:.0f}  "
      f"hard conflicts left={hard_left}")
for nm, u, a, b in zip(names, res.usage, res.accepted, res.new_best):
    print(f"{nm:18s} used={u:5d} accepted={a:5d} new_best={b:3d}")
# Expected: hard conflicts reach 0 and the final cost is pure proximity
# penalty around 270 (start: ~158,600). Usage concentrates on
# move_worst_greedy (~2,600 calls, every one accepted, ~47 new bests) and
# kempe_chain (~1,700 calls); shuffle_slot stays a rare diversifier.
```

## Worked Example 2: Bin Packing with a Heuristic Pool

One-dimensional bin packing: items with sizes $w_i$, bins of capacity $c$, minimize the number of bins $B$. The raw bin count is a terrible search objective — it forms huge plateaus where almost no move changes the value. Use the Falkenauer fitness (Falkenauer 1996, "A Hybrid Grouping Genetic Algorithm for Bin Packing"), here written as a minimization cost:

$$
f(\text{assign}) = 1 - \frac{1}{B} \sum_{b=1}^{B} \left( \frac{\text{fill}_b}{c} \right)^2 ,
$$

which strictly prefers packings with a few very full bins plus one slack bin over evenly half-full bins, giving the search a gradient toward emptying bins. All pool heuristics below are feasibility-preserving, so no penalty term is needed.

```python
"""Bin packing: representation helpers, Falkenauer cost, pool, validator."""
import numpy as np


def compact(assign: np.ndarray) -> np.ndarray:
    """Relabel bins to 0..B-1, removing empty labels."""
    return np.unique(assign, return_inverse=True)[1]


def best_fit_insert(assign: np.ndarray, sizes: np.ndarray, cap: float,
                    items: np.ndarray) -> np.ndarray:
    """Reinsert `items` largest-first into the tightest feasible bin (best fit),
    opening a new bin when nothing fits. Mutates and returns `assign`."""
    assign[items] = -1
    placed = assign >= 0
    n_bins = int(assign[placed].max()) + 1 if placed.any() else 0
    fills = np.bincount(assign[placed], weights=sizes[placed],
                        minlength=n_bins).astype(float)
    for i in items[np.argsort(-sizes[items])]:
        room = cap - fills
        feas = np.flatnonzero(room >= sizes[i] - 1e-9)
        if feas.size:
            b = int(feas[np.argmin(room[feas])])
        else:
            b = fills.size
            fills = np.append(fills, 0.0)
        assign[i] = b
        fills[b] += sizes[i]
    return assign


def falkenauer_cost(assign: np.ndarray, sizes: np.ndarray, cap: float) -> float:
    """1 - mean((fill/cap)^2): plateau-free bin-minimization surrogate."""
    a = compact(assign)
    fills = np.bincount(a, weights=sizes)
    return 1.0 - float(np.mean((fills / cap) ** 2))


def validate_packing(assign: np.ndarray, sizes: np.ndarray, cap: float) -> int:
    """Independent feasibility check; returns the number of bins used."""
    fills = np.bincount(assign, weights=sizes)
    assert (fills <= cap + 1e-9).all(), "capacity violated"
    assert (np.bincount(assign) > 0).all(), "empty bin label present"
    return int(assign.max()) + 1


def make_binpacking_pool(sizes: np.ndarray, cap: float) -> list:
    """Five feasibility-preserving low-level heuristics over an assignment vector."""
    n = sizes.size

    def fills_of(assign: np.ndarray) -> np.ndarray:
        """Current load of every bin label."""
        return np.bincount(assign, weights=sizes,
                           minlength=int(assign.max()) + 1)

    def move_item(assign: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        """Move one random item to the tightest other bin with room."""
        i = int(rng.integers(n))
        room = cap - fills_of(assign)
        room[assign[i]] = -1.0
        feas = np.flatnonzero(room >= sizes[i])
        if feas.size:
            assign[i] = int(feas[np.argmin(room[feas])])
        return compact(assign)

    def swap_items(assign: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        """Swap two items between bins if both capacities allow it."""
        i, j = rng.choice(n, 2, replace=False)
        bi, bj = int(assign[i]), int(assign[j])
        if bi == bj:
            return assign
        f = fills_of(assign)
        if (f[bi] - sizes[i] + sizes[j] <= cap
                and f[bj] - sizes[j] + sizes[i] <= cap):
            assign[i], assign[j] = bj, bi
        return assign

    def empty_smallest(assign: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        """Empty the least-full bin; reinsert its items best-fit-decreasing."""
        items = np.flatnonzero(assign == int(np.argmin(fills_of(assign))))
        return compact(best_fit_insert(assign, sizes, cap, items))

    def empty_two_random(assign: np.ndarray, rng: np.random.Generator) -> np.ndarray:
        """Ruin-and-recreate: empty two random bins, reinsert best-fit-decreasing."""
        labels = np.unique(assign)
        if labels.size < 3:
            return assign
        chosen = rng.choice(labels, size=2, replace=False)
        items = np.flatnonzero(np.isin(assign, chosen))
        return compact(best_fit_insert(assign, sizes, cap, items))

    def largest_from_smallest(assign: np.ndarray,
                              rng: np.random.Generator) -> np.ndarray:
        """Relocate the largest item of the least-full bin (targets bin closure)."""
        f = fills_of(assign)
        b = int(np.argmin(f))
        items = np.flatnonzero(assign == b)
        i = int(items[np.argmax(sizes[items])])
        room = cap - f
        room[b] = -1.0
        feas = np.flatnonzero(room >= sizes[i])
        if feas.size:
            assign[i] = int(feas[np.argmin(room[feas])])
        return compact(assign)

    return [move_item, swap_items, empty_smallest,
            empty_two_random, largest_from_smallest]
```

The pool spans three granularities — single-item moves, pair swaps, and bin-level ruin-and-recreate — which is the bin-packing analogue of mixing small and large neighborhoods. `best_fit_insert` doubles as the constructor (called on an all-unassigned vector it performs best-fit-decreasing, BFD) and as the repair half of the ruin heuristics.

```python
"""Run: BFD start, hyper-heuristic improvement, independent validation."""
import numpy as np

rng = np.random.default_rng(42)
sizes = rng.integers(20, 70, size=60).astype(float)
cap = 100.0
lb = int(np.ceil(sizes.sum() / cap))           # L1 lower bound

x0 = best_fit_insert(np.full(60, -1), sizes, cap, np.arange(60))
bins0 = validate_packing(x0, sizes, cap)

hh = SelectionHyperHeuristic(
    heuristics=make_binpacking_pool(sizes, cap),
    objective=lambda a: falkenauer_cost(a, sizes, cap),
    copy_solution=np.copy,
    selection="epsilon_greedy",
    epsilon=0.2,
    acceptance="late_acceptance",
    la_length=300,
    seed=0,
)
res = hh.run(x0, n_iters=5000)
bins_final = validate_packing(res.best, sizes, cap)
names = ["move_item", "swap_items", "empty_smallest",
         "empty_two_random", "largest_from_smallest"]
print(f"L1 lower bound={lb}  BFD bins={bins0}  HH bins={bins_final}")
print("usage:", dict(zip(names, res.usage)))
# Expected: BFD opens 29 bins; the hyper-heuristic closes one and reaches the
# L1 bound of 28. swap_items dominates usage (~2,800 calls) by earning
# frequent small Falkenauer gains that skew fills; the bin-closing new-best
# events split between swap_items and ruin-and-recreate (empty_two_random).
```

If the run stalls at the BFD value, the usual cause is acceptance too greedy for a plateau-heavy problem — verify the Falkenauer cost is in use (not the bin count) and lengthen the LAHC list before touching the selection parameters.

## Advanced Techniques

### Choice-function selection

The choice function (Cowling et al. 2001) scores each heuristic by three terms: $f_1$ — decayed solo improvement; $f_2$ — decayed improvement when following the previously applied heuristic (captures productive pairs such as ruin→repair); $f_3$ — time since last use (a built-in diversification clock). The modified choice function of Drake, Özcan & Burke (2012, "An Improved Choice Function Heuristic Selection for Cross Domain Heuristic Search") replaces hand-tuned weights with a reactive $\phi$: jump to exploitation on any improvement, drift toward exploration otherwise.

```python
"""Choice-function heuristic selection (Cowling, Kendall & Soubeiga 2001),
with the reactive phi of Drake, Ozcan & Burke (2012)."""
import numpy as np


class ChoiceFunction:
    """Score_i = phi*(f1_i + f2[prev, i]) + (1 - phi)*(time since last use)."""

    def __init__(self, k: int, decay: float = 0.9) -> None:
        self.f1 = np.zeros(k)
        self.f2 = np.zeros((k, k))
        self.last_used = np.zeros(k)
        self.prev: int | None = None
        self.decay = decay
        self.phi = 0.5
        self.it = 0

    def select(self) -> int:
        """Greedy argmax of the choice function."""
        follow = self.f2[self.prev] if self.prev is not None else np.zeros_like(self.f1)
        score = (self.phi * (self.f1 + follow)
                 + (1.0 - self.phi) * (self.it - self.last_used))
        return int(np.argmax(score))

    def update(self, i: int, improvement: float) -> None:
        """Report f(x) - f(x') (positive = better) for the applied heuristic i."""
        self.it += 1
        self.last_used[i] = self.it
        self.f1[i] = improvement + self.decay * self.f1[i]
        if self.prev is not None:
            self.f2[self.prev, i] = improvement + self.decay * self.f2[self.prev, i]
        self.prev = i
        self.phi = 0.99 if improvement > 0 else max(0.01, self.phi - 0.01)


# Wiring into the engine loop: replace `i = self._select(q, n, it)` with
# `i = cf.select()` and the q-update with `cf.update(i, f_cur - f_cand)`.
```

The $f_2$ matrix is the cheapest way to learn heuristic *sequences*; if its entries dominate $f_1$ after a run, the pool has strong pair structure and a sequence-based method (below) is worth trying.

### Sequence-based selection

Single-step credit assumes the next-best heuristic is independent of what was just applied — often false (ruin is only good before repair; a greedy move is best right after a perturbation). The sequence-based selection hyper-heuristic of Kheiri & Keedwell (2015, "A Sequence-based Selection Hyper-heuristic Utilising a Hidden Markov Model") learns transition probabilities between heuristics plus an "end sequence / evaluate now" signal, effectively evolving short macros. The lightweight version: learn the $k \times k$ transition matrix exactly like the choice function's $f_2$, but *sample* the next heuristic from the softmax of the active row instead of taking an argmax, and only invoke acceptance at sequence ends.

### Acceptance engineering

Empirically, swapping the acceptance rule changes cross-domain performance more than swapping the selection rule (consistent with the ALNS finding of Santini et al. 2018 that acceptance criteria matter and simple ones do well — see large-neighborhood-search). Practical rules: LAHC's list length $L$ should scale with the iteration budget ($L \approx$ budget/100 is a sane start); Metropolis needs the warm-up calibration of $T_0$ from an accepted-gap target; accept-all only works when every heuristic embeds its own improvement step. For tight time budgets, record-to-record travel — accept iff $f(x') < f(x^*) + \theta$ with a shrinking $\theta$ — is a robust two-line alternative. Always re-tune acceptance after any change to the objective scaling (penalty weights included).

### Generation hyper-heuristics: evolving a packing rule

Generation hyper-heuristics search the space of heuristics themselves, classically with genetic programming over expression trees of problem features (Burke et al. 2009, "Exploring Hyper-heuristic Methodologies with Genetic Programming"; evolved packing and strip-packing rules that rediscover and beat first fit). The working sketch below keeps the idea and drops the GP machinery: parameterize an *index policy* — score every feasible bin by a weighted feature sum, pack into the argmax — and evolve the weights with a (1+1)-ES on a training set. The evolved object is a fast, reusable *online* heuristic, not a solution.

```text
GENERATION-HH (rule evolution)
  genome w  ->  rule "place item into feasible bin maximizing
                w0*fill + w1*leftover + w2*[perfect fit]"
  fitness(w) = total bins used by the rule over a TRAINING instance set
  search over w with (1+1)-ES (or CMA-ES / GP for richer rule spaces)
  validate the evolved rule on a held-out TEST instance set
```

```python
"""Generation hyper-heuristic sketch: evolve an online bin-packing index rule."""
import numpy as np


def pack_with_rule(w: np.ndarray, sizes: np.ndarray, cap: float) -> int:
    """Pack items in arrival order with the parameterized rule; return bins used."""
    fills: list[float] = []
    for s in sizes:
        best_b, best_score = -1, -np.inf
        for b, f in enumerate(fills):
            leftover = cap - f - s
            if leftover >= -1e-9:
                score = w[0] * f + w[1] * leftover + w[2] * (leftover < 1e-9)
                if score > best_score:
                    best_b, best_score = b, score
        if best_b < 0:
            fills.append(float(s))
        else:
            fills[best_b] += float(s)
    return len(fills)


def evolve_rule(train: list[np.ndarray], cap: float, iters: int = 150,
                seed: int = 0) -> np.ndarray:
    """(1+1)-ES over rule weights; fitness = total bins over training instances."""
    rng = np.random.default_rng(seed)
    w = rng.normal(0.0, 1.0, 3)
    f_w = sum(pack_with_rule(w, s, cap) for s in train)
    for _ in range(iters):
        cand = w + rng.normal(0.0, 0.3, 3)
        f_c = sum(pack_with_rule(cand, s, cap) for s in train)
        if f_c <= f_w:
            w, f_w = cand, f_c
    return w


rng = np.random.default_rng(5)
cap = 100.0
train = [rng.integers(20, 70, 60).astype(float) for _ in range(5)]
test = [rng.integers(20, 70, 60).astype(float) for _ in range(5)]
w_star = evolve_rule(train, cap, iters=150, seed=5)
best_fit = np.array([1.0, 0.0, 0.0])    # 'fullest feasible bin' = best fit
evolved_bins = sum(pack_with_rule(w_star, s, cap) for s in test)
bf_bins = sum(pack_with_rule(best_fit, s, cap) for s in test)
print(f"test bins: evolved={evolved_bins}  best-fit={bf_bins}")
# Expected: evolved=142, best-fit=142 on the held-out instances — the ES
# *rediscovers* best fi

…(truncated)