Upload 11 files

__init__.py

@@ -0,0 +1,35 @@
+"""HPC-Quantize V2: factored Q2_K quantization in Python.
+
+The package keeps the full HPC decision stack (batched GGML seeds, WLS
+superblock solve, factored candidate field, sequential Sieve, 36-state
+Viterbi, spectral/DC shaping, local-safety floors) while replacing the
+original factored solver with a monotone, multi-start, fp16-exact version.
+"""
+
+from .config import HPCConfig
+from .factored import (
+    FactorResult,
+    factor_1d,
+    factor_candidates,
+    factor_only_solve,
+    factor_q2k_block,
+)
+from .pipeline import Q2KQuantizer, QuantizedTensor, TensorStats
+from .q2k import Q2KBlock, QK_K, dequantize_params, pack_codes, unpack_codes
+
+__all__ = [
+    "HPCConfig",
+    "FactorResult",
+    "Q2KBlock",
+    "Q2KQuantizer",
+    "QuantizedTensor",
+    "TensorStats",
+    "QK_K",
+    "dequantize_params",
+    "factor_1d",
+    "factor_candidates",
+    "factor_only_solve",
+    "factor_q2k_block",
+    "pack_codes",
+    "unpack_codes",
+]

__main__.py

@@ -0,0 +1,67 @@
+"""Command-line entry point for the V2 factored Q2_K quantizer.
+
+Examples::
+
+    python -m hpc_quantize_v2 --demo
+    python -m hpc_quantize_v2 --input tensor.npy --output tensor_q2k.npz
+"""
+
+from __future__ import annotations
+
+import argparse
+
+import numpy as np
+
+from .config import HPCConfig
+from .pipeline import Q2KQuantizer
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="HPC-Quantize V2 factored Q2_K")
+    parser.add_argument("--input", help=".npy file containing float weights")
+    parser.add_argument("--importance", help="optional .npy importance matrix")
+    parser.add_argument("--output", help="output .npz with scales/qs/dm")
+    parser.add_argument("--demo", action="store_true", help="quantize synthetic weights")
+    parser.add_argument("--blocks", type=int, default=64, help="demo block count")
+    parser.add_argument("--std", type=float, default=0.02, help="demo weight std")
+    parser.add_argument("--seed", type=int, default=1234, help="demo RNG seed")
+    parser.add_argument("--cand-d", type=int, default=None, help="candidate d count")
+    parser.add_argument("--cand-m", type=int, default=None, help="candidate min count")
+    parser.add_argument("--mse", action="store_true", help="disable HPC sequence search")
+    parser.add_argument("--verbose", action="store_true")
+    args = parser.parse_args()
+
+    if args.demo:
+        rng = np.random.default_rng(args.seed)
+        weights = rng.normal(0.0, args.std, size=(args.blocks, 256)).astype(np.float32)
+        importance = None
+    elif args.input:
+        weights = np.load(args.input)
+        importance = np.load(args.importance) if args.importance else None
+    else:
+        parser.error("pass --demo or --input")
+
+    cfg = HPCConfig(verbose=True)
+    if args.cand_d is not None:
+        cfg.n_cand_d = args.cand_d
+    if args.cand_m is not None:
+        cfg.n_cand_m = args.cand_m
+    if args.mse:
+        cfg.optimizer = "mse"
+
+    result = Q2KQuantizer(cfg).quantize(weights, importance=importance)
+    if args.output:
+        np.savez_compressed(
+            args.output,
+            scales=result.scales,
+            qs=result.qs,
+            dm=result.dm,
+            dequantized=result.dequantize(),
+        )
+        print(f"wrote {args.output}")
+    else:
+        _ = result.dequantize()
+
+
+if __name__ == "__main__":
+    main()

config.py

@@ -0,0 +1,92 @@
+"""Configuration for the HPC-Quantize V2 Q2_K engine.
+
+Every tunable that was a compile-time macro, environment variable or hidden
+global in the C engine is collected here as an explicit, validated field.
+The defaults are the C engine's production defaults unless marked
+``V2`` in the docstring, in which case they are deliberate improvements.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+
+@dataclass(slots=True)
+class HPCConfig:
+    # ---- block geometry -------------------------------------------------
+    qk_k: int = 256
+
+    # ---- candidate field ------------------------------------------------
+    # V2 improvement: the default is the 6x6 symbolic lattice (36 exact
+    # factored solves) instead of the 24x24 physical grid (576 solves).
+    # Set to 24/24 to reproduce the full C candidate field.
+    n_cand_d: int = 6
+    n_cand_m: int = 6
+    factor_pin: int = 3
+    factor_iters: int = 1
+    factor_standalone_iters: int = 16
+    factor_starts: int = 3
+    factor_ulp_radius: int = 1
+    factor_refine_1d: bool = True
+    factor_spectral_refine: bool = True
+    factor_exact_emit: bool = True
+
+    # ---- objective ------------------------------------------------------
+    use_spectral: bool = True
+    dc_lambda: float = 1.0
+    vw_lambda: float = 1.0
+    fold_depth: int = 1
+    fold_gamma: float = 2.0
+    dc_decay: float = 0.85
+    carry_cumulative: bool = True
+    sse_budget: float = 5.0e-4
+    greedy_override_ratio: float = 0.995
+
+    # ---- sieve ----------------------------------------------------------
+    same_penalty: float = 0.85
+    log_slack: float = 1.5
+    quhits_per_block: int = 2
+    graph_max_blocks: int = 8192
+
+    # ---- Viterbi --------------------------------------------------------
+    optimizer: str = "hpc"
+    viterbi_beta: float = 0.25
+    viterbi_alpha: float = 0.08
+    temperature_sample_blocks: int = 128
+
+    # ---- assembly -------------------------------------------------------
+    refine_passes: int = 3
+    ulp_radius: int = 4
+    subblock_states: int = 6
+    subblock_sieve: bool = True
+    use_factor_floor: bool = True
+    use_candidate_floor: bool = True
+    apply_dc_carry: bool = True
+
+    # ---- seed clamp -----------------------------------------------------
+    clamp_kurtosis_off: float = 6.0
+    clamp_kurtosis_soft: float = 20.0
+    clamp_sigma_soft: float = 6.0
+    clamp_sigma_hard: float = 4.0
+
+    # ---- reporting ------------------------------------------------------
+    verbose: bool = False
+
+    def validate(self) -> "HPCConfig":
+        if self.qk_k != 256:
+            raise ValueError("Q2_K requires qk_k == 256")
+        if not (1 <= self.n_cand_d <= 24 and 1 <= self.n_cand_m <= 24):
+            raise ValueError("candidate counts must lie in [1, 24]")
+        if self.optimizer not in ("hpc", "mse"):
+            raise ValueError("optimizer must be 'hpc' or 'mse'")
+        if not 0 <= self.factor_pin <= 3:
+            raise ValueError("factor_pin must be a 2-bit mask")
+        if self.fold_depth < 0:
+            raise ValueError("fold_depth must be non-negative")
+        if not 0.0 <= self.dc_decay <= 1.0:
+            raise ValueError("dc_decay must lie in [0, 1]")
+        if self.ulp_radius < 0 or self.factor_ulp_radius < 0:
+            raise ValueError("ULP radii must be non-negative")
+        if self.refine_passes < 1:
+            raise ValueError("refine_passes must be >= 1")
+        return self

factored.py

@@ -0,0 +1,547 @@
+"""Grid-free factored Q2_K solver (V2).
+
+The reference solver alternates four exact coordinate solves::
+
+    q       | (d, Ls, dmin, Lm)
+    s_j,t_j | q
+    d,Ls    | s
+    dmin,Lm | t
+
+V2 keeps that decomposition but fixes three weaknesses:
+
+1. ``factor_1d`` is no longer a fixed six-iteration Lloyd heuristic.  It
+   starts from several deterministic candidates and finishes with an exact
+   integer-code coordinate refinement, so its weighted error is monotone.
+2. The block solve is multi-start, accepts only objective-decreasing states,
+   and returns the best snapshot rather than the last iterate.
+3. The emitted fp16 ``(d, dmin)`` pair is polished in its own domain, and the
+   spectral/DC shaping stage is applied inside the solver.  The returned
+   objective therefore equals the objective of the bytes that are written.
+
+The result object is immutable and self-describing, replacing the C
+engine's many out-parameters.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from .config import HPCConfig
+from .fp16 import from_fp16_bits, fp16_ulp_neighbors, nearest_int, round_fp16, to_fp16_bits
+from .q2k import NS, QK_K, codes_from_params, derive_subscales
+from .spectral import objective, refine_codes_hadamard
+
+N_CAND_D = 24
+N_CAND_M = 24
+TOTAL_CANDIDATES = N_CAND_D * N_CAND_M
+
+MULTS_D = np.array(
+    [
+        0.780,
+        0.835,
+        0.880,
+        0.915,
+        0.943,
+        0.963,
+        0.978,
+        0.988,
+        0.994,
+        0.997,
+        0.999,
+        1.000,
+        1.002,
+        1.005,
+        1.011,
+        1.021,
+        1.035,
+        1.054,
+        1.080,
+        1.115,
+        1.160,
+        1.215,
+        1.275,
+        1.340,
+    ],
+    dtype=np.float32,
+)
+MULTS_M = np.array(
+    [
+        0.750,
+        0.800,
+        0.840,
+        0.870,
+        0.900,
+        0.920,
+        0.940,
+        0.955,
+        0.970,
+        0.985,
+        0.995,
+        1.000,
+        1.005,
+        1.015,
+        1.030,
+        1.045,
+        1.060,
+        1.080,
+        1.100,
+        1.130,
+        1.160,
+        1.200,
+        1.250,
+        1.300,
+    ],
+    dtype=np.float32,
+)
+
+CAND_TO_QUHIT = np.array(
+    [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5],
+    dtype=np.int64,
+)
+
+
+@dataclass(slots=True)
+class FactorResult:
+    d: np.float32
+    dmin: np.float32
+    d16: np.uint16
+    m16: np.uint16
+    ls: np.ndarray
+    lm: np.ndarray
+    codes: np.ndarray
+    sse: float
+    objective: float
+    iterations: int
+    starts: int
+
+    def subscales(self) -> np.ndarray:
+        return (self.ls & 0x0F) | ((self.lm & 0x0F) << 4)
+
+
+def _factor_1d_lloyd(v: np.ndarray, w: np.ndarray, nmax: int, iters: int):
+    vmax = float(v.max())
+    if vmax < 1.0e-30:
+        return 0.0, np.zeros(v.size, dtype=np.uint8)
+    best_d = vmax / nmax
+    best_codes = np.clip(nearest_int(v / best_d), 0, nmax)
+    best_obj = float(np.sum(w * (v - best_d * best_codes) ** 2))
+
+    starts = [vmax / nmax]
+    init_codes = np.clip(nearest_int(v / starts[0]), 0, nmax)
+    den0 = float(np.sum(w * init_codes * init_codes))
+    if den0 > 1.0e-30:
+        starts.append(float(np.sum(w * v * init_codes) / den0))
+
+    seen: set[float] = set()
+    for d0 in starts:
+        if not np.isfinite(d0) or d0 <= 0.0:
+            continue
+        key = float(np.float32(d0))
+        if key in seen:
+            continue
+        seen.add(key)
+        d = float(d0)
+        codes = np.clip(nearest_int(v / d), 0, nmax)
+        for _ in range(iters):
+            den = float(np.sum(w * codes * codes))
+            if den < 1.0e-30:
+                break
+            dn = float(np.sum(w * v * codes) / den)
+            if dn <= 0.0:
+                break
+            conv = abs(dn - d) <= 1.0e-7 * abs(d) + 1.0e-30
+            d = dn
+            codes = np.clip(nearest_int(v / d), 0, nmax)
+            if conv:
+                break
+        obj = float(np.sum(w * (v - d * codes) ** 2))
+        if obj < best_obj:
+            best_obj = obj
+            best_d = d
+            best_codes = codes.copy()
+    return best_d, best_codes
+
+
+def _refine_integer_codes(v, w, d, codes, nmax, sweeps=2):
+    v = np.asarray(v, dtype=np.float64)
+    w = np.asarray(w, dtype=np.float64)
+    codes = np.asarray(codes, dtype=np.int64).copy()
+    d = float(d)
+    if d <= 0.0:
+        return d, codes.astype(np.uint8)
+    obj = float(np.sum(w * (v - d * codes) ** 2))
+    options = np.arange(nmax + 1, dtype=np.int64)
+    for _ in range(sweeps):
+        improved = False
+        for i in range(v.size):
+            wi, vi, li = w[i], v[i], codes[i]
+            sum_lv = float(np.sum(w * v * codes))
+            sum_l2 = float(np.sum(w * codes * codes))
+            slv = sum_lv - wi * vi * li + wi * vi * options
+            sl2 = sum_l2 - wi * li * li + wi * options * options
+            valid = sl2 > 1.0e-30
+            if not np.any(valid):
+                continue
+            dd = np.where(valid, slv / np.where(valid, sl2, 1.0), 0.0)
+            trial = codes[None, :].repeat(options.size, axis=0)
+            trial[:, i] = options
+            e = v[None, :] - dd[:, None] * trial.astype(np.float64)
+            scores = np.sum(w[None, :] * e * e, axis=1)
+            scores = np.where(valid & (dd > 0.0), scores, np.inf)
+            k = int(np.argmin(scores))
+            if scores[k] < obj - 1.0e-18 and options[k] != li:
+                obj = float(scores[k])
+                d = float(dd[k])
+                codes[i] = options[k]
+                improved = True
+        if not improved:
+            break
+    return d, codes.astype(np.uint8)
+
+
+def factor_1d(v, w, nmax: int = 15, iters: int = 8, refine: bool = True):
+    """Weighted 1-D factorisation ``v ~= d * L`` with ``L in [0, nmax]``."""
+    v = np.maximum(np.asarray(v, dtype=np.float64), 0.0)
+    w = np.maximum(np.asarray(w, dtype=np.float64), 0.0)
+    d, codes = _factor_1d_lloyd(v, w, nmax, iters)
+    if refine and d > 0.0:
+        d, codes = _refine_integer_codes(v, w, d, codes, nmax)
+    return d, codes
+
+
+def _range_init(x, w, cfg: HPCConfig):
+    x = np.asarray(x, dtype=np.float64)
+    w = np.asarray(w, dtype=np.float64)
+    blocks = x.reshape(NS, 16)
+    wblocks = w.reshape(NS, 16)
+    s = np.zeros(NS)
+    t = np.zeros(NS)
+    codes = np.zeros(QK_K, dtype=np.uint8)
+    for j in range(NS):
+        sub = blocks[j]
+        xmin = min(float(sub.min()), 0.0)
+        xmax = float(sub.max())
+        span = xmax - xmin
+        sj = span / 3.0 if span > 1.0e-30 else 0.0
+        tj = -xmin
+        s[j] = sj
+        t[j] = tj
+        if sj > 1.0e-30:
+            codes[16 * j : 16 * (j + 1)] = np.clip(
+                nearest_int((sub + tj) / sj), 0, 3
+            ).astype(np.uint8)
+    sw = wblocks.sum(axis=1) if w is not None else np.full(NS, 16.0)
+    d, ls = factor_1d(s, sw, 15, cfg.factor_iters, cfg.factor_refine_1d)
+    m, lm = factor_1d(t, sw, 15, cfg.factor_iters, cfg.factor_refine_1d)
+    return d, m, ls, lm, codes
+
+
+def _seeded_init(x, scales, mins, d_seed, m_seed):
+    d = float(d_seed)
+    m = float(m_seed)
+    ls, lm = derive_subscales(scales, mins, d, m)
+    codes = codes_from_params(x, d, m, ls, lm)
+    return d, m, ls, lm, codes
+
+
+def _ls_init(x, w, cfg: HPCConfig):
+    x = np.asarray(x, dtype=np.float64)
+    w = np.asarray(w, dtype=np.float64)
+    blocks = x.reshape(NS, 16)
+    codes = np.zeros(QK_K, dtype=np.uint8)
+    for j in range(NS):
+        sub = blocks[j]
+        xmin = min(float(sub.min()), 0.0)
+        xmax = float(sub.max())
+        span = xmax - xmin
+        if span > 1.0e-30:
+            codes[16 * j : 16 * (j + 1)] = np.clip(
+                nearest_int((sub - xmin) * 3.0 / span), 0, 3
+            ).astype(np.uint8)
+    s, t = _sub_ls(x, w, codes)
+    sw = w.reshape(NS, 16).sum(axis=1)
+    d, ls = factor_1d(s, sw, 15, cfg.factor_iters, cfg.factor_refine_1d)
+    m, lm = factor_1d(t, sw, 15, cfg.factor_iters, cfg.factor_refine_1d)
+    codes = codes_from_params(x, d, m, ls, lm)
+    return d, m, ls, lm, codes
+
+
+def _sub_ls(x, w, codes):
+    x = np.asarray(x, dtype=np.float64)
+    w = np.asarray(w, dtype=np.float64)
+    q = np.asarray(codes, dtype=np.float64).reshape(NS, 16)
+    xb = x.reshape(NS, 16)
+    wb = w.reshape(NS, 16)
+    s = np.zeros(NS)
+    t = np.zeros(NS)
+    for j in range(NS):
+        ww = wb[j]
+        qq = q[j]
+        xx = xb[j]
+        Sqq = float(np.sum(ww * qq * qq))
+        Swq = float(np.sum(ww * qq))
+        Sw = float(np.sum(ww))
+        Sqx = float(np.sum(ww * qq * xx))
+        Sx = float(np.sum(ww * xx))
+        det = Sqq * Sw - Swq * Swq
+        if det > 1.0e-30:
+            sj = (Sqx * Sw - Swq * Sx) / det
+            tj = (sj * Swq - Sx) / Sw
+        else:
+            sj = 0.0
+            tj = 0.0
+        if sj < 0.0:
+            sj = 0.0
+        if tj < 0.0:
+            tj = 0.0
+            sj = Sqx / Sqq if Sqq > 1.0e-30 else 0.0
+            if sj < 0.0:
+                sj = 0.0
+        s[j] = sj
+        t[j] = tj
+    return s, t
+
+
+def _solve(
+    x,
+    w,
+    cfg: HPCConfig,
+    d,
+    m,
+    ls,
+    lm,
+    codes,
+    pin: int,
+    use_spectral: bool,
+    max_iter: int,
+    dc_carry: float,
+):
+    x = np.asarray(x, dtype=np.float32)
+    w = np.asarray(w, dtype=np.float32)
+    ls = np.asarray(ls, dtype=np.uint8).copy()
+    lm = np.asarray(lm, dtype=np.uint8).copy()
+    codes = np.asarray(codes, dtype=np.uint8).copy()
+    d = float(d)
+    m = float(m)
+
+    obj, _ = objective(x, w, d, m, ls, lm, codes, cfg, dc_carry)
+    best = (obj, d, m, ls.copy(), lm.copy(), codes.copy())
+    stall = 0
+    iters = 0
+    for it in range(max(1, max_iter)):
+        iters = it + 1
+        codes = codes_from_params(x, d, m, ls, lm)
+        s, t = _sub_ls(x, w, codes)
+        q = codes.astype(np.float64).reshape(NS, 16)
+        wb = w.reshape(NS, 16)
+        ws = np.sum(wb * q * q, axis=1)
+        wt = np.sum(wb, axis=1)
+
+        if pin & 1:
+            if d > 1.0e-30:
+                ls = np.clip(nearest_int(s / d), 0, 15).astype(np.uint8)
+            else:
+                ls = np.zeros(NS, dtype=np.uint8)
+        else:
+            d_new, ls_new = factor_1d(
+                s, ws, 15, cfg.factor_iters, cfg.factor_refine_1d
+            )
+            d, ls = float(d_new), ls_new
+
+        if pin & 2:
+            if m > 1.0e-30:
+                lm = np.clip(nearest_int(t / m), 0, 15).astype(np.uint8)
+            else:
+                lm = np.zeros(NS, dtype=np.uint8)
+        else:
+            m_new, lm_new = factor_1d(
+                t, wt, 15, cfg.factor_iters, cfg.factor_refine_1d
+            )
+            m, lm = float(m_new), lm_new
+
+        codes = codes_from_params(x, d, m, ls, lm)
+        cur, _ = objective(x, w, d, m, ls, lm, codes, cfg, dc_carry)
+        if cur < best[0] - 1.0e-15:
+            best = (cur, d, m, ls.copy(), lm.copy(), codes.copy())
+            stall = 0
+        else:
+            stall += 1
+        if stall >= 2:
+            break
+
+    obj, d, m, ls, lm, codes = best
+    return d, m, ls, lm, codes, iters
+
+
+def _emit_polish(x, w, cfg, d, m, ls, lm, codes, dc_carry):
+    d16 = np.uint16(to_fp16_bits(np.float32(d)))
+    m16 = np.uint16(to_fp16_bits(np.float32(m)))
+    d = float(from_fp16_bits(d16))
+    m = float(from_fp16_bits(m16))
+    codes = codes_from_params(x, d, m, ls, lm)
+    obj, sse = objective(x, w, d, m, ls, lm, codes, cfg, dc_carry)
+    best = (obj, sse, d, m, d16, m16, ls.copy(), lm.copy(), codes.copy())
+
+    if cfg.factor_exact_emit and cfg.factor_ulp_radius > 0:
+        d_choices = fp16_ulp_neighbors(d16, cfg.factor_ulp_radius)
+        m_choices = fp16_ulp_neighbors(m16, cfg.factor_ulp_radius)
+        for db in d_choices:
+            dd = float(from_fp16_bits(db))
+            for mb in m_choices:
+                mm = float(from_fp16_bits(mb))
+                c = codes_from_params(x, dd, mm, ls, lm)
+                o, s = objective(x, w, dd, mm, ls, lm, c, cfg, dc_carry)
+                if o < best[0] - 1.0e-15:
+                    best = (o, s, dd, mm, np.uint16(db), np.uint16(mb), ls.copy(), lm.copy(), c)
+    obj, sse, d, m, d16, m16, ls, lm, codes = best
+
+    if cfg.factor_spectral_refine:
+        refined = refine_codes_hadamard(
+            x, w, d, m, ls, lm, codes, cfg, dc_carry
+        )
+        o, s = objective(x, w, d, m, ls, lm, refined, cfg, dc_carry)
+        if o < obj - 1.0e-15:
+            obj, sse, codes = o, s, refined
+    return d, m, d16, m16, ls, lm, codes, sse, obj
+
+
+def factor_q2k_block(
+    x,
+    w=None,
+    cfg: HPCConfig | None = None,
+    init_s=None,
+    init_t=None,
+    d_seed: float = 0.0,
+    m_seed: float = 0.0,
+    pin: int = 0,
+    use_spectral: bool | None = None,
+    max_iter: int | None = None,
+    starts: int | None = None,
+    dc_carry: float = 0.0,
+) -> FactorResult:
+    cfg = (cfg or HPCConfig()).validate()
+    x = np.asarray(x, dtype=np.float32).reshape(QK_K)
+    w = (
+        np.ones(QK_K, dtype=np.float32)
+        if w is None
+        else np.asarray(w, dtype=np.float32).reshape(QK_K)
+    )
+    use_spectral = cfg.use_spectral if use_spectral is None else use_spectral
+    max_iter = cfg.factor_standalone_iters if max_iter is None else max_iter
+    starts = cfg.factor_starts if starts is None else starts
+
+    candidates = []
+    if init_s is not None and init_t is not None and d_seed > 1.0e-30 and m_seed > 1.0e-30:
+        d, m, ls, lm, codes = _seeded_init(x, init_s, init_t, d_seed, m_seed)
+        candidates.append((d, m, ls, lm, codes))
+    else:
+        d, m, ls, lm, codes = _range_init(x, w, cfg)
+        if d_seed > 1.0e-30:
+            d = float(d_seed)
+            if init_s is not None:
+                ls = np.clip(nearest_int(np.asarray(init_s) / d), 0, 15).astype(np.uint8)
+        if m_seed > 1.0e-30:
+            m = float(m_seed)
+            if init_t is not None:
+                lm = np.clip(nearest_int(np.asarray(init_t) / m), 0, 15).astype(np.uint8)
+        candidates.append((d, m, ls, lm, codes))
+        if starts > 1:
+            candidates.append(_ls_init(x, w, cfg))
+
+    best = None
+    total_iters = 0
+    for d, m, ls, lm, codes in candidates:
+        d, m, ls, lm, codes, it = _solve(
+            x, w, cfg, d, m, ls, lm, codes, pin, use_spectral, max_iter, dc_carry
+        )
+        total_iters += it
+        d, m, d16, m16, ls, lm, codes, sse, obj = _emit_polish(
+            x, w, cfg, d, m, ls, lm, codes, dc_carry
+        )
+        if best is None or obj < best.objective - 1.0e-15:
+            best = FactorResult(
+                d=np.float32(d),
+                dmin=np.float32(m),
+                d16=np.uint16(d16),
+                m16=np.uint16(m16),
+                ls=ls.astype(np.uint8),
+                lm=lm.astype(np.uint8),
+                codes=codes.astype(np.uint8),
+                sse=float(sse),
+                objective=float(obj),
+                iterations=total_iters,
+                starts=len(candidates),
+            )
+    return best
+
+
+def factor_candidates(
+    x,
+    w,
+    scales,
+    mins,
+    base_d,
+    base_m,
+    cfg: HPCConfig,
+    dc_carry: float = 0.0,
+) -> np.ndarray:
+    """Factored candidate field indexed exactly like the C 24x24 lattice."""
+    out = np.full(TOTAL_CANDIDATES, np.inf, dtype=np.float32)
+    nd = int(cfg.n_cand_d)
+    nm = int(cfg.n_cand_m)
+    for a in range(nd):
+        di = a if nd == N_CAND_D else (a * N_CAND_D) // nd
+        ds = float(round_fp16(np.float32(base_d) * MULTS_D[di]))
+        for b in range(nm):
+            mi = b if nm == N_CAND_M else (b * N_CAND_M) // nm
+            ms = float(round_fp16(np.float32(base_m) * MULTS_M[mi]))
+            res = factor_q2k_block(
+                x,
+                w,
+                cfg,
+                init_s=scales,
+                init_t=mins,
+                d_seed=ds,
+                m_seed=ms,
+                pin=cfg.factor_pin,
+                use_spectral=cfg.use_spectral,
+                max_iter=cfg.factor_iters,
+                starts=1,
+                dc_carry=dc_carry,
+            )
+            if nd == N_CAND_D and nm == N_CAND_M:
+                out[di * N_CAND_M + mi] = np.float32(res.objective)
+            else:
+                cb = di >> 2
+                fb = mi >> 2
+                for aa in range(4):
+                    for bb in range(4):
+                        idx = (cb * 4 + aa) * N_CAND_M + fb * 4 + bb
+                        if res.objective < out[idx]:
+                            out[idx] = np.float32(res.objective)
+    return out
+
+
+def factor_only_solve(
+    x, w, cfg: HPCConfig, dc_carry: float = 0.0
+) -> FactorResult:
+    return factor_q2k_block(
+        x,
+        w,
+        cfg,
+        pin=0,
+        use_spectral=cfg.use_spectral,
+        max_iter=cfg.factor_standalone_iters,
+        starts=cfg.factor_starts,
+        dc_carry=dc_carry,
+    )
+
+
+def candidate_pair(base_d, base_m, cidx: int) -> tuple[np.uint16, np.uint16]:
+    di = int(cidx) // N_CAND_M
+    mi = int(cidx) % N_CAND_M
+    d16 = np.uint16(to_fp16_bits(np.float32(base_d) * MULTS_D[di]))
+    m16 = np.uint16(to_fp16_bits(np.float32(base_m) * MULTS_M[mi]))
+    return d16, m16

fp16.py

@@ -0,0 +1,55 @@
+"""Bit-exact GGML numerical primitives.
+
+The C engine relies on two non-obvious primitives:
+
+* ``gguf_nearest_int`` implements round-to-nearest-even with the classic
+  ``float + 1.5*2**23`` mantissa trick rather than ``rint``.  The trick is
+  reproduced here so that every quantized code is bit-identical to the
+  reference implementation.
+* ``gguf_fp32_to_fp16`` / ``gguf_fp16_to_fp32`` are the emitted parameter
+  domain.  Every solver in V2 may optimise in float32, but it reports and
+  stores only fp16-rounded values.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+_MAGIC = np.float32(12582912.0)
+_F16_MAX_BITS = 0x7BFF
+
+
+def nearest_int(value) -> np.ndarray:
+    """GGML-compatible nearest-integer rounding (ties to even)."""
+    arr = np.asarray(value, dtype=np.float32)
+    shifted = arr + _MAGIC
+    bits = shifted.view(np.int32)
+    return ((bits & 0x007FFFFF) - 0x00400000).astype(np.int64)
+
+
+def to_fp16_bits(value) -> np.ndarray:
+    arr = np.asarray(value, dtype=np.float32)
+    return arr.astype(np.float16).view(np.uint16).astype(np.uint16)
+
+
+def from_fp16_bits(bits) -> np.ndarray:
+    arr = np.asarray(bits, dtype=np.uint16)
+    return arr.view(np.float16).astype(np.float32)
+
+
+def round_fp16(value) -> np.ndarray:
+    return from_fp16_bits(to_fp16_bits(value))
+
+
+def fp16_ulp_neighbors(bits, radius: int) -> np.ndarray:
+    """All non-negative finite fp16 bit patterns within ``radius`` ULPs."""
+    if radius < 0:
+        raise ValueError("radius must be non-negative")
+    b = int(bits)
+    lo = max(0, b - radius)
+    hi = min(_F16_MAX_BITS, b + radius)
+    return np.arange(lo, hi + 1, dtype=np.uint16)
+
+
+def epsilon_for(scale: np.ndarray, relative: float = 1.0e-7) -> np.ndarray:
+    return relative * np.abs(scale) + np.float32(1.0e-30)

pipeline.py

@@ -0,0 +1,559 @@
+"""End-to-end Q2_K pipeline preserving the HPC machinery.
+
+Pipeline stages (mirroring the C engine, with V2 fixes):
+
+1. kurtosis-gated outlier clamp and batched GGML seeds;
+2. vectorised WLS superblock solve;
+3. factored candidate field + standalone factored floor;
+4. sequential Sieve over coarse/fine symbolic states;
+5. 36-state Viterbi sequence optimisation;
+6. candidate-aware factor consistency pass;
+7. rolling DC residual carry;
+8. sub-block Sieve assembly, LS refit, fp16-ULP polish, spectral shaping;
+9. objective-consistent local-safety floors and DC nudge.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from math import ceil
+
+import numpy as np
+
+from .config import HPCConfig
+from .factored import (
+    CAND_TO_QUHIT,
+    N_CAND_M,
+    TOTAL_CANDIDATES,
+    FactorResult,
+    candidate_pair,
+    factor_candidates,
+    factor_only_solve,
+    factor_q2k_block,
+)
+from .fp16 import from_fp16_bits, nearest_int, round_fp16, to_fp16_bits
+from .q2k import (
+    NS,
+    QK_K,
+    Q2KBlock,
+    codes_from_params,
+    codes_from_params_batch,
+    dequantize_params,
+    dequantize_params_batch,
+    derive_subscales,
+    make_block,
+    pack_codes,
+    pack_scales,
+    unpack_codes,
+    weighted_sse,
+    weighted_sse_batch,
+)
+from .seeds import build_seeds, seed_clamp_value
+from .sieve import block_chain_edges, sieve_measure_chain, subblock_chain_edges
+from .spectral import dc_nudge_codes, objective, refine_codes_hadamard
+from .viterbi import robust_temperature, viterbi_select
+
+
+@dataclass(slots=True)
+class TensorStats:
+    n_blocks: int
+    total_sse: float
+    weighted_sse: float
+    rmse: float
+    sigma: float
+    rmse_over_sigma: float
+    temperature: float
+    factor_floor_used: int
+    candidate_floor_used: int
+
+
+@dataclass(slots=True)
+class QuantizedTensor:
+    blocks: list[Q2KBlock]
+    scales: np.ndarray
+    qs: np.ndarray
+    dm: np.ndarray
+    best_candidate: np.ndarray
+    factor_objectives: np.ndarray
+    config: HPCConfig
+    stats: TensorStats
+
+    def dequantize(self) -> np.ndarray:
+        return np.concatenate([b.dequantize() for b in self.blocks])
+
+
+def _phase2_wls(X, W, seeds, clamp_val):
+    n = X.shape[0]
+    clipped = np.clip(X, -clamp_val, clamp_val).astype(np.float32)
+    base_d = seeds["d"].astype(np.float64).copy()
+    base_m = seeds["dmin"].astype(np.float64).copy()
+    ls = seeds["ls"].copy()
+    lm = seeds["lm"].copy()
+    scales = seeds["scales"].astype(np.float64)
+    mins = seeds["mins"].astype(np.float64)
+    ww = np.ones_like(X, dtype=np.float64) if W is None else W.astype(np.float64)
+
+    for _ in range(5):
+        codes = codes_from_params_batch(clipped, base_d, base_m, ls, lm).astype(np.float64)
+        A = np.repeat(ls.astype(np.float64), 16, axis=1) * codes
+        B = np.repeat(lm.astype(np.float64), 16, axis=1)
+        Saa = np.sum(ww * A * A, axis=1)
+        Sab = np.sum(ww * A * B, axis=1)
+        Sbb = np.sum(ww * B * B, axis=1)
+        Sxa = np.sum(ww * clipped.astype(np.float64) * A, axis=1)
+        Sxb = np.sum(ww * clipped.astype(np.float64) * B, axis=1)
+        det = Saa * Sbb - Sab * Sab
+        valid = det > 1.0e-30
+        safe_det = np.where(valid, det, 1.0)
+        d_new = (Sbb * Sxa - Sab * Sxb) / safe_det
+        m_new = (Sab * Sxa - Saa * Sxb) / safe_det
+        upd_d = valid & (d_new > 0.0) & (d_new < 4.0 * (seeds["d"] + 1.0e-10))
+        upd_m = valid & (m_new > 0.0) & (m_new < 4.0 * (seeds["dmin"] + 1.0e-10))
+        base_d = np.where(upd_d, round_fp16(d_new.astype(np.float32)).astype(np.float64), base_d)
+        base_m = np.where(upd_m, round_fp16(m_new.astype(np.float32)).astype(np.float64), base_m)
+
+        safe_d = np.where(base_d > 1.0e-15, base_d, 1.0)
+        safe_m = np.where(base_m > 1.0e-15, base_m, 1.0)
+        ls = np.where(
+            base_d[:, None] > 1.0e-15,
+            np.clip(nearest_int(scales / safe_d[:, None]), 0, 15),
+            0,
+        ).astype(np.uint8)
+        lm = np.where(
+            base_m[:, None] > 1.0e-15,
+            np.clip(nearest_int(mins / safe_m[:, None]), 0, 15),
+            0,
+        ).astype(np.uint8)
+
+    return base_d.astype(np.float32), base_m.astype(np.float32)
+
+
+def _bin_geometry():
+    di = np.arange(TOTAL_CANDIDATES) // N_CAND_M
+    mi = np.arange(TOTAL_CANDIDATES) % N_CAND_M
+    return CAND_TO_QUHIT[di], CAND_TO_QUHIT[mi]
+
+
+def _sieve_marginals(candidate_errors, graph_blocks, stride, temperature, cfg):
+    qi_d, qi_m = _bin_geometry()
+    coarse = np.zeros((graph_blocks, 6), dtype=np.float64)
+    fine = np.zeros((graph_blocks, 6), dtype=np.float64)
+    n_blocks = candidate_errors.shape[0]
+    for gi in range(graph_blocks):
+        start = gi * stride
+        stop = min((gi + 1) * stride, n_blocks)
+        group = candidate_errors[start:stop].astype(np.float64)
+        finite = np.isfinite(group)
+        if np.any(finite):
+            agg = np.where(finite, group, np.nan).mean(axis=0)
+            valid = np.isfinite(agg)
+            mn = float(np.min(agg[valid])) if np.any(valid) else 0.0
+            coeff = np.where(
+                valid,
+                np.exp(-(agg - mn) / (2.0 * max(temperature, 1.0e-12))),
+                0.0,
+            )
+        else:
+            coeff = np.zeros(TOTAL_CANDIDATES, dtype=np.float64)
+        for v in range(6):
+            coarse[gi, v] = float(np.sum(coeff[qi_d == v]))
+            fine[gi, v] = float(np.sum(coeff[qi_m == v]))
+        cn = float(np.linalg.norm(coarse[gi]))
+        fn = float(np.linalg.norm(fine[gi]))
+        if cn > 1.0e-30:
+            coarse[gi] /= cn
+        if fn > 1.0e-30:
+            fine[gi] /= fn
+
+    amps = np.zeros((2 * graph_blocks, 6), dtype=np.float64)
+    amps[0::2] = coarse
+    amps[1::2] = fine
+    marg, _ = sieve_measure_chain(
+        amps,
+        block_chain_edges(graph_blocks),
+        same_penalty=cfg.same_penalty,
+        log_slack=cfg.log_slack,
+        deterministic=True,
+    )
+    return marg[0::2], marg[1::2]
+
+
+def _phase3_select(candidate_errors, cfg):
+    n_blocks = candidate_errors.shape[0]
+    if cfg.optimizer != "hpc" or n_blocks < 2:
+        best = np.argmin(candidate_errors, axis=1).astype(np.int64)
+        return best, 0.0, np.zeros(n_blocks, dtype=np.int64)
+
+    stride = max(1, ceil(n_blocks / cfg.graph_max_blocks))
+    graph_blocks = ceil(n_blocks / stride)
+    temperature = robust_temperature(candidate_errors, cfg)
+    coarse_marg, fine_marg = _sieve_marginals(
+        candidate_errors, graph_blocks, stride, temperature, cfg
+    )
+    best, path = viterbi_select(
+        candidate_errors,
+        coarse_marg,
+        fine_marg,
+        temperature,
+        stride,
+        graph_blocks,
+        cfg,
+    )
+    return best, temperature, path
+
+
+def _compute_dc_carry(X, seeds, best, base_d, base_m, cfg, row_width):
+    n = X.shape[0]
+    carry = np.zeros(n, dtype=np.float32)
+    blocks_per_row = row_width // QK_K if row_width and row_width % QK_K == 0 else 0
+    rolling = 0.0
+    for i in range(n):
+        if blocks_per_row > 0 and i % blocks_per_row == 0:
+            rolling = 0.0
+        d16, m16 = candidate_pair(base_d[i], base_m[i], best[i])
+        d = float(from_fp16_bits(d16))
+        m = float(from_fp16_bits(m16))
+        ls, lm = derive_subscales(seeds["scales"][i], seeds["mins"][i], d, m)
+        codes = codes_from_params(X[i], d, m, ls, lm)
+        xh = dequantize_params(d, m, ls, lm, codes)
+        carry[i] = np.float32(cfg.dc_decay * rolling)
+        dc_res = float(np.sum(X[i].astype(np.float64) - xh.astype(np.float64)))
+        rolling = carry[i] + dc_res if cfg.carry_cumulative else dc_res
+    return carry
+
+
+def _per_sub_top_states(x, w, d, m, k=6):
+    k = max(1, min(int(k), 6))
+    xb = np.asarray(x, dtype=np.float32).reshape(NS, 16)
+    if w is None:
+        wb = np.ones((NS, 16), dtype=np.float32)
+    else:
+        wb = np.asarray(w, dtype=np.float32).reshape(NS, 16)
+    state_ls = np.zeros((NS, k), dtype=np.uint8)
+    state_lm = np.zeros((NS, k), dtype=np.uint8)
+    state_err = np.full((NS, k), np.inf, dtype=np.float64)
+    ls_opts = np.arange(16, dtype=np.float32)[:, None]
+    lm_opts = np.arange(16, dtype=np.float32)[None, :]
+    d_sub = np.float32(d) * ls_opts
+    m_sub = np.float32(m) * lm_opts
+    safe = d_sub[:, :, None] >= np.float32(1.0e-15)
+    safe_d = np.where(safe, d_sub[:, :, None], np.float32(1.0))
+    for j in range(NS):
+        sub = xb[j]
+        wj = wb[j]
+        raw = (sub[None, None, :] + m_sub[:, :, None]) / safe_d
+        q = np.clip(nearest_int(raw), 0, 3).astype(np.float32)
+        q = np.where(safe, q, 0.0)
+        deq = d_sub[:, :, None] * q - m_sub[:, :, None]
+        err = np.sum(wj[None, None, :] * (sub[None, None, :] - deq) ** 2, axis=2)
+        flat = err.ravel()
+        order = np.argsort(flat, kind="stable")[:k]
+        state_err[j] = flat[order]
+        state_ls[j] = (order // 16).astype(np.uint8)
+        state_lm[j] = (order % 16).astype(np.uint8)
+    return state_ls, state_lm, state_err
+
+
+def _select_sub_states(state_ls, state_lm, state_err, cfg):
+    if not cfg.subblock_sieve:
+        return state_ls[:, 0].copy(), state_lm[:, 0].copy()
+    amps = np.zeros((NS, 6), dtype=np.float64)
+    for j in range(NS):
+        err0 = state_err[j, 0]
+        spread = float(state_err[j, 5] - state_err[j, 0])
+        temp = spread * 0.3 if spread > 1.0e-15 else 0.1
+        temp = max(temp, 1.0e-12)
+        amp = np.exp(-(state_err[j] - err0) / temp)
+        norm = float(np.linalg.norm(amp))
+        amps[j] = amp / norm if norm > 1.0e-30 else 1.0 / np.sqrt(6.0)
+    marg, _ = sieve_measure_chain(
+        amps,
+        subblock_chain_edges(NS),
+        same_penalty=cfg.same_penalty,
+        log_slack=cfg.log_slack,
+        deterministic=True,
+    )
+    pick = np.argmax(marg, axis=1)
+    ls = np.array([state_ls[j, pick[j]] for j in range(NS)], dtype=np.uint8)
+    lm = np.array([state_lm[j, pick[j]] for j in range(NS)], dtype=np.uint8)
+    return ls, lm
+
+
+def _ls_refit(x, w, d, m, ls, lm, codes, max_d, max_m):
+    x = np.asarray(x, dtype=np.float64)
+    w = np.ones(QK_K, dtype=np.float64) if w is None else np.asarray(w, dtype=np.float64)
+    A = np.repeat(np.asarray(ls, dtype=np.float64), 16) * np.asarray(codes, dtype=np.float64)
+    B = np.repeat(np.asarray(lm, dtype=np.float64), 16)
+    Saa = float(np.sum(w * A * A))
+    Sab = float(np.sum(w * A * B))
+    Sbb = float(np.sum(w * B * B))
+    Sxa = float(np.sum(w * x * A))
+    Sxb = float(np.sum(w * x * B))
+    det = Saa * Sbb - Sab * Sab
+    if det <= 1.0e-30:
+        return float(d), float(m)
+    d_new = (Sbb * Sxa - Sab * Sxb) / det
+    m_new = (Sab * Sxa - Saa * Sxb) / det
+    if d_new <= 0.0 or d_new >= 4.0 * (max_d + 1.0e-10):
+        d_new = d
+    if m_new <= 0.0 or m_new >= 4.0 * (max_m + 1.0e-10):
+        m_new = m
+    return float(round_fp16(np.float32(d_new))), float(round_fp16(np.float32(m_new)))
+
+
+def _ulp_polish(x, w, d, m, ls, lm, radius):
+    if radius <= 0:
+        codes = codes_from_params(x, d, m, ls, lm)
+        return float(d), float(m), codes, weighted_sse(x, w, d, m, ls, lm, codes)
+    d16 = np.uint16(to_fp16_bits(np.float32(d)))
+    m16 = np.uint16(to_fp16_bits(np.float32(m)))
+    best = None
+    for db in range(max(0, int(d16) - radius), min(0x7BFF, int(d16) + radius) + 1):
+        dd = float(from_fp16_bits(np.uint16(db)))
+        for mb in range(max(0, int(m16) - radius), min(0x7BFF, int(m16) + radius) + 1):
+            mm = float(from_fp16_bits(np.uint16(mb)))
+            c = codes_from_params(x, dd, mm, ls, lm)
+            s = weighted_sse(x, w, dd, mm, ls, lm, c)
+            if best is None or s < best[0] - 1.0e-15:
+                best = (s, dd, mm, c)
+    return best[1], best[2], best[3], best[0]
+
+
+def _per_sub_scan(x, w, d, m):
+    state_ls, state_lm, _ = _per_sub_top_states(x, w, d, m, k=1)
+    return state_ls[:, 0], state_lm[:, 0]
+
+
+def _assemble_block(x, w, seeds, base_d, base_m, cidx, dc_carry, cfg):
+    d16, m16 = candidate_pair(base_d, base_m, cidx)
+    d = float(from_fp16_bits(d16))
+    m = float(from_fp16_bits(m16))
+    ls, lm = derive_subscales(seeds["scales"], seeds["mins"], d, m)
+    max_d, max_m = d, m
+    prev = None
+    for _ in range(5):
+        s_ls, s_lm, s_err = _per_sub_top_states(x, w, d, m, cfg.subblock_states)
+        ls, lm = _select_sub_states(s_ls, s_lm, s_err, cfg)
+        codes = codes_from_params(x, d, m, ls, lm)
+        d, m = _ls_refit(x, w, d, m, ls, lm, codes, max_d, max_m)
+        cur = (int(to_fp16_bits(np.float32(d))), int(to_fp16_bits(np.float32(m))))
+        if cur == prev:
+            break
+        prev = cur
+
+    d, m, codes, _ = _ulp_polish(x, w, d, m, ls, lm, cfg.ulp_radius)
+    ls, lm = _per_sub_scan(x, w, d, m)
+    codes = codes_from_params(x, d, m, ls, lm)
+    codes = refine_codes_hadamard(x, w, d, m, ls, lm, codes, cfg, dc_carry)
+    obj, sse = objective(x, w, d, m, ls, lm, codes, cfg, dc_carry)
+    return d, m, ls, lm, codes, obj, sse
+
+
+def _candidate_config(seeds, base_d, base_m, cidx):
+    d16, m16 = candidate_pair(base_d, base_m, cidx)
+    d = float(from_fp16_bits(d16))
+    m = float(from_fp16_bits(m16))
+    ls, lm = derive_subscales(seeds["scales"], seeds["mins"], d, m)
+    return d, m, ls, lm
+
+
+def _apply_dc_nudge(X, W, blocks, cfg, row_width):
+    n = X.shape[0]
+    blocks_per_row = row_width // QK_K if row_width and row_width % QK_K == 0 else 0
+    rolling = 0.0
+    for i in range(n):
+        if blocks_per_row > 0 and i % blocks_per_row == 0:
+            rolling = 0.0
+        block = blocks[i]
+        carry = np.float32(cfg.dc_decay * rolling)
+        d = float(block.d_float)
+        m = float(block.dmin_float)
+        ls, lm = block.unpack_scales()
+        codes = unpack_codes(block.qs)
+        w = None if W is None else W[i]
+        codes = dc_nudge_codes(X[i], d, m, ls, lm, codes, w, float(carry), cfg)
+        blocks[i] = Q2KBlock(
+            scales=pack_scales(ls, lm),
+            qs=pack_codes(codes),
+            d=block.d,
+            dmin=block.dmin,
+        )
+        xh = blocks[i].dequantize()
+        dc_res = float(np.sum(X[i].astype(np.float64) - xh.astype(np.float64)))
+        rolling = float(carry) + dc_res if cfg.carry_cumulative else dc_res
+    return blocks
+
+
+class Q2KQuantizer:
+    def __init__(self, config: HPCConfig | None = None):
+        self.config = (config or HPCConfig()).validate()
+
+    def quantize(
+        self,
+        weights: np.ndarray,
+        importance: np.ndarray | None = None,
+        row_width: int | None = None,
+    ) -> QuantizedTensor:
+        cfg = self.config
+        x = np.asarray(weights, dtype=np.float32)
+        if x.size % QK_K != 0:
+            raise ValueError("weight count must be a multiple of 256")
+        if x.ndim == 2 and row_width is None and x.shape[1] % QK_K == 0:
+            row_width = int(x.shape[1])
+        x = x.ravel()
+        n_blocks = x.size // QK_K
+        X = x.reshape(n_blocks, QK_K)
+
+        if importance is not None:
+            imp = np.asarray(importance, dtype=np.float32).ravel()
+            if imp.size != x.size:
+                raise ValueError("importance must match the weight count")
+            W = imp.reshape(n_blocks, QK_K)
+        else:
+            W = None
+
+        seeds = build_seeds(X, W, cfg)
+        clamp_val = seed_clamp_value(x, cfg)
+        base_d, base_m = _phase2_wls(X, W, seeds, clamp_val)
+
+        candidate_errors = np.full((n_blocks, TOTAL_CANDIDATES), np.inf, dtype=np.float32)
+        standalone: list[FactorResult] = []
+        for i in range(n_blocks):
+            wi = None if W is None else W[i]
+            candidate_errors[i] = factor_candidates(
+                X[i],
+                wi,
+                seeds["scales"][i],
+                seeds["mins"][i],
+                base_d[i],
+                base_m[i],
+                cfg,
+            )
+            standalone.append(factor_only_solve(X[i], wi, cfg))
+
+        best, temperature, _ = _phase3_select(candidate_errors, cfg)
+        carry = _compute_dc_carry(X, seeds, best, base_d, base_m, cfg, row_width)
+
+        at_slot: list[FactorResult] = []
+        for i in range(n_blocks):
+            wi = None if W is None else W[i]
+            d16, m16 = candidate_pair(base_d[i], base_m[i], best[i])
+            at_slot.append(
+                factor_q2k_block(
+                    X[i],
+                    wi,
+                    cfg,
+                    init_s=seeds["scales"][i],
+                    init_t=seeds["mins"][i],
+                    d_seed=float(from_fp16_bits(d16)),
+                    m_seed=float(from_fp16_bits(m16)),
+                    pin=cfg.factor_pin,
+                    use_spectral=cfg.use_spectral,
+                    max_iter=cfg.factor_iters,
+                    starts=1,
+                    dc_carry=0.0,
+                )
+            )
+
+        blocks: list[Q2KBlock] = []
+        factor_floor_used = 0
+        candidate_floor_used = 0
+        for i in range(n_blocks):
+            wi = None if W is None else W[i]
+            d, m, ls, lm, codes, obj, sse = _assemble_block(
+                X[i], wi, seeds, base_d[i], base_m[i], best[i], carry[i], cfg
+            )
+
+            if cfg.use_candidate_floor:
+                cd, cm, cls, clm = _candidate_config(
+                    {k: v[i] for k, v in seeds.items() if k in ("scales", "mins")},
+                    base_d[i],
+                    base_m[i],
+                    best[i],
+                )
+                ccodes = codes_from_params(X[i], cd, cm, cls, clm)
+                cobj, csse = objective(
+                    X[i], wi, cd, cm, cls, clm, ccodes, cfg, float(carry[i])
+                )
+                if cobj < obj - 1.0e-15:
+                    d, m, ls, lm, codes, obj, sse = cd, cm, cls, clm, ccodes, cobj, csse
+                    candidate_floor_used += 1
+
+            if cfg.use_factor_floor:
+                fac = at_slot[i]
+                fobj, fsse = objective(
+                    X[i], wi, fac.d, fac.dmin, fac.ls, fac.lm, fac.codes, cfg, float(carry[i])
+                )
+                if fobj < obj - 1.0e-15:
+                    d = float(fac.d)
+                    m = float(fac.dmin)
+                    ls = fac.ls.copy()
+                    lm = fac.lm.copy()
+                    codes = fac.codes.copy()
+                    obj, sse = fobj, fsse
+                    factor_floor_used += 1
+
+            blocks.append(make_block(d, m, ls, lm, codes))
+
+        if cfg.apply_dc_carry and (cfg.dc_lambda > 0.0 or cfg.dc_decay > 0.0):
+            blocks = _apply_dc_nudge(X, W, blocks, cfg, row_width)
+
+        total_sse = 0.0
+        weighted = 0.0
+        for i, b in enumerate(blocks):
+            e = X[i].astype(np.float64) - b.dequantize().astype(np.float64)
+            total_sse += float(np.dot(e, e))
+            if W is None:
+                weighted += float(np.dot(e, e))
+            else:
+                weighted += float(np.dot(W[i].astype(np.float64), e * e))
+
+        sigma = float(np.sqrt(np.mean(x.astype(np.float64) ** 2)))
+        rmse = float(np.sqrt(total_sse / x.size))
+        stats = TensorStats(
+            n_blocks=n_blocks,
+            total_sse=total_sse,
+            weighted_sse=weighted,
+            rmse=rmse,
+            sigma=sigma,
+            rmse_over_sigma=rmse / sigma if sigma > 1.0e-15 else 0.0,
+            temperature=temperature,
+            factor_floor_used=factor_floor_used,
+            candidate_floor_used=candidate_floor_used,
+        )
+        scales, qs, dm = _stack_blocks(blocks)
+        factor_objectives = np.array([r.objective for r in standalone], dtype=np.float64)
+        result = QuantizedTensor(
+            blocks=blocks,
+            scales=scales,
+            qs=qs,
+            dm=dm,
+            best_candidate=best,
+            factor_objectives=factor_objectives,
+            config=cfg,
+            stats=stats,
+        )
+        if cfg.verbose:
+            _report(stats)
+        return result
+
+
+def _stack_blocks(blocks):
+    n = len(blocks)
+    scales = np.empty((n, NS), dtype=np.uint8)
+    qs = np.empty((n, QK_K // 4), dtype=np.uint8)
+    dm = np.empty((n, 2), dtype=np.uint16)
+    for i, b in enumerate(blocks):
+        scales[i] = b.scales
+        qs[i] = b.qs
+        dm[i, 0] = b.d
+        dm[i, 1] = b.dmin
+    return scales, qs, dm
+
+
+def _report(stats: TensorStats) -> None:
+    print(
+        f"Q2_K V2: blocks={stats.n_blocks} rmse={stats.rmse:.6e} "
+        f"rmse/sigma={stats.rmse_over_sigma:.5f} temp={stats.temperature:.3e} "
+        f"floors(cand={stats.candidate_floor_used}, factor={stats.factor_floor_used})"
+    )

q2k.py

@@ -0,0 +1,223 @@
+"""Bit-exact Q2_K block representation, packing and reconstruction.
+
+Layout (84 bytes per 256 weights)::
+
+    scales[16]  uint8   low nibble = Ls[j], high nibble = Lm[j]
+    qs[64]      uint8   four 2-bit codes per byte
+    d           fp16    super-block scale
+    dmin        fp16    super-block minimum scale
+
+The decoder reconstructs::
+
+    x_hat[16*j + k] = d * Ls[j] * L[16*j + k] - dmin * Lm[j]
+
+with ``L in {0,1,2,3}`` and ``Ls, Lm in {0..15}``.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from .fp16 import from_fp16_bits, nearest_int, to_fp16_bits
+
+QK_K = 256
+NS = QK_K // 16
+NQ = 4
+
+
+@dataclass(slots=True)
+class Q2KBlock:
+    scales: np.ndarray
+    qs: np.ndarray
+    d: np.uint16
+    dmin: np.uint16
+
+    @property
+    def d_float(self) -> np.float32:
+        return np.float32(from_fp16_bits(np.uint16(self.d)))
+
+    @property
+    def dmin_float(self) -> np.float32:
+        return np.float32(from_fp16_bits(np.uint16(self.dmin)))
+
+    def unpack_scales(self) -> tuple[np.ndarray, np.ndarray]:
+        scales = np.asarray(self.scales, dtype=np.uint8)
+        return scales & 0x0F, scales >> 4
+
+    def dequantize(self) -> np.ndarray:
+        ls, lm = self.unpack_scales()
+        return dequantize_params(
+            self.d_float, self.dmin_float, ls, lm, unpack_codes(self.qs)
+        )
+
+    def error(self, x: np.ndarray) -> float:
+        e = np.asarray(x, dtype=np.float64) - self.dequantize().astype(np.float64)
+        return float(np.dot(e, e))
+
+
+def unpack_codes(qs: np.ndarray) -> np.ndarray:
+    qs = np.asarray(qs, dtype=np.uint8)
+    codes = np.empty(QK_K, dtype=np.uint8)
+    for base in (0, 128):
+        q = qs[base // 4 : base // 4 + 32]
+        for shift_level in range(4):
+            lo = base + 32 * shift_level
+            codes[lo : lo + 32] = (q >> np.uint8(2 * shift_level)) & np.uint8(3)
+    return codes
+
+
+def pack_codes(codes: np.ndarray) -> np.ndarray:
+    codes = np.asarray(codes, dtype=np.uint8)
+    if codes.shape != (QK_K,):
+        raise ValueError("codes must have shape (256,)")
+    qs = np.empty(QK_K // 4, dtype=np.uint8)
+    for base in (0, 128):
+        q = (
+            codes[base : base + 32]
+            | (codes[base + 32 : base + 64] << np.uint8(2))
+            | (codes[base + 64 : base + 96] << np.uint8(4))
+            | (codes[base + 96 : base + 128] << np.uint8(6))
+        )
+        qs[base // 4 : base // 4 + 32] = q
+    return qs
+
+
+def pack_scales(ls: np.ndarray, lm: np.ndarray) -> np.ndarray:
+    return (np.asarray(ls, dtype=np.uint8) & 0x0F) | (
+        (np.asarray(lm, dtype=np.uint8) & 0x0F) << np.uint8(4)
+    )
+
+
+def dequantize_params(
+    d, dmin, ls: np.ndarray, lm: np.ndarray, codes: np.ndarray
+) -> np.ndarray:
+    ls = np.asarray(ls, dtype=np.float32)
+    lm = np.asarray(lm, dtype=np.float32)
+    codes = np.asarray(codes, dtype=np.float32)
+    d_sub = np.repeat(np.float32(d) * ls, 16)
+    m_sub = np.repeat(np.float32(dmin) * lm, 16)
+    return d_sub * codes - m_sub
+
+
+def codes_from_params(
+    x: np.ndarray, d, dmin, ls: np.ndarray, lm: np.ndarray
+) -> np.ndarray:
+    x = np.asarray(x, dtype=np.float32)
+    d_sub = np.repeat(np.float32(d) * np.asarray(ls, dtype=np.float32), 16)
+    m_sub = np.repeat(np.float32(dmin) * np.asarray(lm, dtype=np.float32), 16)
+    safe = d_sub >= np.float32(1.0e-15)
+    codes = np.zeros(QK_K, dtype=np.uint8)
+    if np.any(safe):
+        q = nearest_int((x[safe] + m_sub[safe]) / d_sub[safe])
+        codes[safe] = np.clip(q, 0, 3).astype(np.uint8)
+    return codes
+
+
+def derive_subscales(
+    scales: np.ndarray, mins: np.ndarray, d, dmin
+) -> tuple[np.ndarray, np.ndarray]:
+    scales = np.asarray(scales, dtype=np.float32)
+    mins = np.asarray(mins, dtype=np.float32)
+    if float(d) > 1.0e-15:
+        ls = np.clip(nearest_int(scales / np.float32(d)), 0, 15).astype(np.uint8)
+    else:
+        ls = np.zeros(NS, dtype=np.uint8)
+    if float(dmin) > 1.0e-15:
+        lm = np.clip(nearest_int(mins / np.float32(dmin)), 0, 15).astype(np.uint8)
+    else:
+        lm = np.zeros(NS, dtype=np.uint8)
+    return ls, lm
+
+
+def weighted_sse(
+    x: np.ndarray,
+    w: np.ndarray | None,
+    d,
+    dmin,
+    ls: np.ndarray,
+    lm: np.ndarray,
+    codes: np.ndarray,
+) -> float:
+    xh = dequantize_params(d, dmin, ls, lm, codes)
+    e = np.asarray(x, dtype=np.float64) - xh.astype(np.float64)
+    if w is None:
+        return float(np.dot(e, e))
+    return float(np.dot(np.asarray(w, dtype=np.float64), e * e))
+
+
+def make_block(d, dmin, ls, lm, codes) -> Q2KBlock:
+    return Q2KBlock(
+        scales=pack_scales(ls, lm),
+        qs=pack_codes(codes),
+        d=np.uint16(to_fp16_bits(np.float32(d))),
+        dmin=np.uint16(to_fp16_bits(np.float32(dmin))),
+    )
+
+
+def blocks_to_arrays(blocks: list[Q2KBlock]) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    n = len(blocks)
+    scales = np.empty((n, NS), dtype=np.uint8)
+    qs = np.empty((n, QK_K // 4), dtype=np.uint8)
+    dm = np.empty((n, 2), dtype=np.uint16)
+    for i, b in enumerate(blocks):
+        scales[i] = b.scales
+        qs[i] = b.qs
+        dm[i, 0] = b.d
+        dm[i, 1] = b.dmin
+    return scales, qs, dm
+
+
+def dequantize_blocks(blocks: list[Q2KBlock]) -> np.ndarray:
+    if not blocks:
+        return np.empty(0, dtype=np.float32)
+    return np.concatenate([b.dequantize() for b in blocks])
+
+
+def codes_from_params_batch(
+    x: np.ndarray, d: np.ndarray, dmin: np.ndarray, ls: np.ndarray, lm: np.ndarray
+) -> np.ndarray:
+    x = np.asarray(x, dtype=np.float32)
+    n = x.shape[0]
+    d_sub = np.repeat(
+        np.asarray(d, dtype=np.float32)[:, None] * np.asarray(ls, dtype=np.float32), 16, axis=1
+    )
+    m_sub = np.repeat(
+        np.asarray(dmin, dtype=np.float32)[:, None] * np.asarray(lm, dtype=np.float32), 16, axis=1
+    )
+    safe = d_sub >= np.float32(1.0e-15)
+    codes = np.zeros((n, QK_K), dtype=np.uint8)
+    if np.any(safe):
+        q = nearest_int((x[safe] + m_sub[safe]) / d_sub[safe])
+        codes[safe] = np.clip(q, 0, 3).astype(np.uint8)
+    return codes
+
+
+def dequantize_params_batch(
+    d: np.ndarray, dmin: np.ndarray, ls: np.ndarray, lm: np.ndarray, codes: np.ndarray
+) -> np.ndarray:
+    n = codes.shape[0]
+    d_sub = np.repeat(
+        np.asarray(d, dtype=np.float32)[:, None] * np.asarray(ls, dtype=np.float32), 16, axis=1
+    )
+    m_sub = np.repeat(
+        np.asarray(dmin, dtype=np.float32)[:, None] * np.asarray(lm, dtype=np.float32), 16, axis=1
+    )
+    return d_sub * np.asarray(codes, dtype=np.float32) - m_sub
+
+
+def weighted_sse_batch(
+    x: np.ndarray,
+    w: np.ndarray | None,
+    d: np.ndarray,
+    dmin: np.ndarray,
+    ls: np.ndarray,
+    lm: np.ndarray,
+    codes: np.ndarray,
+) -> np.ndarray:
+    xh = dequantize_params_batch(d, dmin, ls, lm, codes)
+    e = np.asarray(x, dtype=np.float64) - xh.astype(np.float64)
+    if w is None:
+        return np.sum(e * e, axis=1)
+    return np.sum(np.asarray(w, dtype=np.float64) * e * e, axis=1)

seeds.py

@@ -0,0 +1,190 @@
+"""Batched GGML seed quantizers and activation-aware seed weighting.
+
+The C engine performs ``make_qkx2_quants`` once per 16-weight sub-block and
+``make_qp_quants`` once per 16-sub-block superblock.  V2 batches both across
+every sub-block / superblock in the tensor with NumPy, which removes the
+per-block Python overhead while preserving the exact search rules.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .config import HPCConfig
+from .fp16 import nearest_int
+
+
+def make_qkx2_batch(
+    x: np.ndarray, w: np.ndarray, nmax: int = 3
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Batched ``make_qkx2_quants``.  Returns (scale, min_abs, codes)."""
+    x = np.asarray(x, dtype=np.float64)
+    w = np.asarray(w, dtype=np.float64)
+    rows, n = x.shape
+    xmin = x.min(axis=1)
+    xmax = x.max(axis=1)
+    xmin = np.minimum(xmin, 0.0)
+    span = xmax - xmin
+    degenerate = (xmax == xmin) | (span <= 1.0e-30)
+
+    sum_w = w.sum(axis=1)
+    sum_x = (w * x).sum(axis=1)
+
+    safe_span = np.where(span > 1.0e-30, span, 1.0)
+    iscale = nmax / safe_span
+    scale = 1.0 / iscale
+    codes = np.clip(nearest_int(iscale[:, None] * (x - xmin[:, None])), 0, nmax)
+    best_mad = np.sum(w * np.abs(scale[:, None] * codes + xmin[:, None] - x), axis=1)
+    best_scale = scale.copy()
+    best_min = xmin.copy()
+    best_codes = codes.copy()
+
+    for step in range(16):
+        try_iscale = (-0.5 + 0.1 * step + nmax) / safe_span
+        aux = np.clip(nearest_int(try_iscale[:, None] * (x - xmin[:, None])), 0, nmax)
+        sl = (w * aux).sum(axis=1)
+        sl2 = (w * aux * aux).sum(axis=1)
+        sxl = (w * aux * x).sum(axis=1)
+        det = sum_w * sl2 - sl * sl
+        valid = det > 1.0e-30
+        safe_det = np.where(valid, det, 1.0)
+        this_scale = (sum_w * sxl - sum_x * sl) / safe_det
+        this_min = (sl2 * sum_x - sl * sxl) / safe_det
+        positive_min = this_min > 0.0
+        this_min = np.where(positive_min, 0.0, this_min)
+        this_scale = np.where(
+            positive_min & (sl2 > 1.0e-30), sxl / np.where(sl2 > 1.0e-30, sl2, 1.0), this_scale
+        )
+        mad = np.sum(w * np.abs(this_scale[:, None] * aux + this_min[:, None] - x), axis=1)
+        take = valid & (mad < best_mad)
+        if np.any(take):
+            best_mad = np.where(take, mad, best_mad)
+            best_scale = np.where(take, this_scale, best_scale)
+            best_min = np.where(take, this_min, best_min)
+            best_codes[take] = aux[take]
+
+    best_scale = np.where(degenerate, 0.0, best_scale)
+    best_codes[degenerate] = 0
+    best_min = np.where(degenerate, xmin, best_min)
+    return best_scale, -best_min, best_codes.astype(np.uint8)
+
+
+def make_qp_batch(
+    values: np.ndarray, weights: np.ndarray, nmax: int = 15
+) -> tuple[np.ndarray, np.ndarray]:
+    """Batched ``make_qp_quants``.  Factors positive values into d * code."""
+    x = np.asarray(values, dtype=np.float64)
+    w = np.asarray(weights, dtype=np.float64)
+    rows, n = x.shape
+    xmax = x.max(axis=1)
+    dead = xmax < 1.0e-15
+    safe_max = np.where(dead, 1.0, xmax)
+
+    iscale = nmax / safe_max
+    codes = np.clip(nearest_int(iscale[:, None] * x), 0, nmax)
+    scale = 1.0 / iscale
+    best_mse = np.sum(w * (x - scale[:, None] * codes) ** 2, axis=1)
+    best_iscale = iscale.copy()
+
+    for step in range(-4, 5):
+        if step == 0:
+            continue
+        trial_iscale = (0.1 * step + nmax) / safe_max
+        trial_scale = 1.0 / trial_iscale
+        trial_codes = np.clip(nearest_int(trial_iscale[:, None] * x), 0, nmax)
+        mse = np.sum(w * (x - trial_scale[:, None] * trial_codes) ** 2, axis=1)
+        take = mse < best_mse
+        if np.any(take):
+            best_mse = np.where(take, mse, best_mse)
+            best_iscale = np.where(take, trial_iscale, best_iscale)
+
+    codes = np.clip(nearest_int(best_iscale[:, None] * x), 0, nmax)
+    sumlx = (w * x * codes).sum(axis=1)
+    suml2 = (w * codes * codes).sum(axis=1)
+
+    for _ in range(5):
+        changed = False
+        for i in range(n):
+            wi = w[:, i]
+            xi = x[:, i]
+            li = codes[:, i].astype(np.float64)
+            slx = sumlx - wi * xi * li
+            sl2 = suml2 - wi * li * li
+            valid = (slx > 0.0) & (sl2 > 0.0)
+            safe_slx = np.where(valid, slx, 1.0)
+            new_l = np.clip(nearest_int(xi * sl2 / safe_slx), 0, nmax)
+            trial = valid & (new_l != codes[:, i])
+            if not np.any(trial):
+                continue
+            slx2 = slx + wi * xi * new_l
+            sl2_2 = sl2 + wi * new_l * new_l
+            better = trial & (slx2 * slx2 * suml2 > sumlx * sumlx * sl2_2)
+            if np.any(better):
+                codes[better, i] = new_l[better].astype(np.uint8)
+                sumlx = np.where(better, slx2, sumlx)
+                suml2 = np.where(better, sl2_2, suml2)
+                changed = True
+        if not changed:
+            break
+
+    d = np.where(suml2 > 1.0e-30, sumlx / np.where(suml2 > 1.0e-30, suml2, 1.0), 0.0)
+    d = np.where(dead, 0.0, d)
+    codes[dead] = 0
+    return d, codes.astype(np.uint8)
+
+
+def seed_clamp_value(weights: np.ndarray, cfg: HPCConfig) -> float:
+    x = np.asarray(weights, dtype=np.float64).ravel()
+    mean_sq = float(np.mean(x * x))
+    if mean_sq <= 1.0e-30:
+        return float("inf")
+    mean_fourth = float(np.mean(x * x * x * x))
+    kurtosis = mean_fourth / (mean_sq * mean_sq)
+    sigma = float(np.sqrt(mean_sq))
+    if kurtosis <= cfg.clamp_kurtosis_off:
+        return float("inf")
+    if kurtosis <= cfg.clamp_kurtosis_soft:
+        return cfg.clamp_sigma_soft * sigma
+    return cfg.clamp_sigma_hard * sigma
+
+
+def seed_weights(
+    block: np.ndarray, importance: np.ndarray | None, sigma_sq: float, clamp_val: float
+) -> np.ndarray:
+    if importance is not None:
+        return np.asarray(importance, dtype=np.float64)
+    return np.sqrt(sigma_sq + np.minimum(block, clamp_val) ** 2)
+
+
+def build_seeds(
+    x: np.ndarray, importance: np.ndarray | None, cfg: HPCConfig
+) -> dict[str, np.ndarray]:
+    """Phase-1 seeds for all blocks: per-sub (scale,min) and shared (d,dmin)."""
+    n_blocks = x.shape[0]
+    blocks = np.asarray(x, dtype=np.float64).reshape(n_blocks, 16, 16)
+    if importance is None:
+        sigma_sq = np.mean(blocks * blocks, axis=(1, 2), keepdims=True)
+        weights = np.sqrt(sigma_sq + blocks * blocks)
+    else:
+        weights = np.asarray(importance, dtype=np.float64).reshape(n_blocks, 16, 16)
+
+    flat_x = blocks.reshape(-1, 16)
+    flat_w = weights.reshape(-1, 16)
+    scales, mins, codes = make_qkx2_batch(flat_x, flat_w, nmax=3)
+    scales = scales.reshape(n_blocks, 16)
+    mins = mins.reshape(n_blocks, 16)
+    codes = codes.reshape(n_blocks, 256)
+    sw = weights.sum(axis=2)
+    d, ls = make_qp_batch(scales, sw, nmax=15)
+    dmin, lm = make_qp_batch(mins, sw, nmax=15)
+    return {
+        "scales": scales.astype(np.float32),
+        "mins": mins.astype(np.float32),
+        "codes": codes.astype(np.uint8),
+        "sw": sw.astype(np.float32),
+        "d": d.astype(np.float32),
+        "dmin": dmin.astype(np.float32),
+        "ls": ls.astype(np.uint8),
+        "lm": lm.astype(np.uint8),
+        "weights": weights.astype(np.float32),
+    }

sieve.py

@@ -0,0 +1,106 @@
+"""Sequential six-state Sieve, ported from the HPC graph measurement.
+
+The C engine builds an :class:`HPCGraph` and runs a column-first sequential
+measurement with real back-action.  Only the following semantics are used by
+the Q2_K pipeline, so V2 implements them directly on a chain:
+
+* local Born probabilities ``|amplitude|^2``;
+* a multiplicative neighbour mix with ``same_penalty`` for equal bins;
+* a log-domain slack filter;
+* deterministic MAP collapse (the production setting);
+* back-action that biases still-unmeasured neighbours away from the outcome.
+
+This is mathematically identical to the production path while removing the
+graph allocation, edge bookkeeping and triality view machinery.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+
+def sieve_measure_chain(
+    amplitudes: np.ndarray,
+    edges: list[tuple[int, int]],
+    same_penalty: float = 0.85,
+    log_slack: float = 1.5,
+    deterministic: bool = True,
+    rng: np.random.Generator | None = None,
+) -> tuple[np.ndarray, np.ndarray]:
+    amps = np.array(amplitudes, dtype=np.float64, copy=True)
+    n_sites, n_states = amps.shape
+    if n_states != 6:
+        raise ValueError("sieve expects six symbolic states per site")
+    adj: list[list[int]] = [[] for _ in range(n_sites)]
+    for eid, (a, b) in enumerate(edges):
+        adj[a].append(eid)
+        adj[b].append(eid)
+    alive = np.ones(len(edges), dtype=bool)
+    marginals = np.zeros((n_sites, 6), dtype=np.float64)
+    measured = np.zeros(n_sites, dtype=np.int64)
+    generator = rng or np.random.default_rng(271828)
+
+    for site in range(n_sites - 1, -1, -1):
+        local = amps[site] * amps[site]
+        total = float(local.sum())
+        probs0 = local / total if total > 1.0e-30 else np.full(6, 1.0 / 6.0)
+        log_score = np.log(probs0 + 1.0e-30)
+
+        for eid in adj[site]:
+            if not alive[eid]:
+                continue
+            a, b = edges[eid]
+            partner = b if a == site else a
+            q = amps[partner] * amps[partner]
+            qt = float(q.sum())
+            q = q / qt if qt > 1.0e-30 else np.full(6, 1.0 / 6.0)
+            mix = np.empty(6)
+            for d in range(6):
+                mix[d] = sum(q[w] * (same_penalty if w == d else 1.0) for w in range(6))
+            log_score += np.log(mix + 1.0e-30)
+
+        best = float(log_score.max())
+        probs = np.where(best - log_score > log_slack, 0.0, np.exp(log_score - best))
+        total = float(probs.sum())
+        probs = probs / total if total > 1.0e-30 else np.full(6, 1.0 / 6.0)
+        marginals[site] = probs
+
+        if deterministic:
+            outcome = int(np.argmax(probs))
+        else:
+            outcome = int(np.searchsorted(np.cumsum(probs), generator.random(), side="right"))
+            outcome = min(outcome, 5)
+        measured[site] = outcome
+
+        amps[site] = 0.0
+        amps[site, outcome] = 1.0
+        for eid in adj[site]:
+            if not alive[eid]:
+                continue
+            a, b = edges[eid]
+            partner = b if a == site else a
+            amps[partner, outcome] *= same_penalty
+            norm = float(np.linalg.norm(amps[partner]))
+            if norm > 1.0e-30:
+                amps[partner] /= norm
+        for eid in adj[site]:
+            alive[eid] = False
+        adj[site] = []
+
+    return marginals, measured
+
+
+def block_chain_edges(n_blocks: int) -> list[tuple[int, int]]:
+    """The production two-site-per-block chain, including cross-block links."""
+    edges: list[tuple[int, int]] = []
+    for i in range(n_blocks):
+        s0, s1 = 2 * i, 2 * i + 1
+        edges.append((s0, s1))
+        if i + 1 < n_blocks:
+            edges.append((s0, 2 * (i + 1)))
+            edges.append((s1, 2 * (i + 1) + 1))
+    return edges
+
+
+def subblock_chain_edges(n_subblocks: int) -> list[tuple[int, int]]:
+    return [(i, i + 1) for i in range(n_subblocks - 1)]

spectral.py

@@ -0,0 +1,226 @@
+"""HPC spectral objective: folded error pyramid, DC/vesica penalty, shaping.
+
+The C engine scores a block with::
+
+    E = sum_i w_i e_i^2
+        + (dc_lambda / n) * (sum_i e_i)^2
+        + (vw_lambda / n) * sum_{pairs} (e_i + e_{i+n/2})^2
+
+and the fold tree generalises the vesica term to every periodic level.  V2
+keeps the exact same objective and adds a vectorised implementation of the
+D6/vesica greedy code descent.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .config import HPCConfig
+from .q2k import QK_K, dequantize_params, weighted_sse
+
+
+def fold_weights(n: int, cfg: HPCConfig) -> np.ndarray:
+    levels = int(np.log2(n))
+    if (1 << levels) != n:
+        raise ValueError("fold pyramid requires a power-of-two length")
+    weights = np.zeros(levels, dtype=np.float64)
+    gamma = 1.0
+    for k in range(1, levels + 1):
+        if k == levels:
+            weights[k - 1] = cfg.dc_lambda / n
+        elif k <= cfg.fold_depth:
+            weights[k - 1] = (cfg.vw_lambda / n) * gamma
+        gamma *= cfg.fold_gamma
+    return weights
+
+
+def fold_build(errors: np.ndarray) -> np.ndarray:
+    e = np.asarray(errors, dtype=np.float64)
+    n = e.size
+    py = np.empty(n - 1, dtype=np.float64)
+    src = e
+    offset = 0
+    m = n // 2
+    while m >= 1:
+        dst = src[:m] + src[m : 2 * m]
+        py[offset : offset + m] = dst
+        src = dst
+        offset += m
+        m //= 2
+    return py
+
+
+def fold_energy(pyramid: np.ndarray, weights: np.ndarray) -> float:
+    acc = 0.0
+    offset = 0
+    m = pyramid.size // 2 + 1
+    level = 0
+    while m >= 1:
+        if weights[level] != 0.0:
+            v = pyramid[offset : offset + m]
+            acc += float(weights[level] * np.dot(v, v))
+        offset += m
+        m //= 2
+        level += 1
+    return acc
+
+
+def spectral_penalty(
+    errors: np.ndarray, cfg: HPCConfig, dc_carry: float = 0.0
+) -> float:
+    if not cfg.use_spectral or (cfg.dc_lambda == 0.0 and cfg.vw_lambda == 0.0):
+        return 0.0
+    pyramid = fold_build(errors)
+    pyramid[-1] += dc_carry
+    return fold_energy(pyramid, fold_weights(errors.size, cfg))
+
+
+def objective(
+    x: np.ndarray,
+    w: np.ndarray | None,
+    d,
+    dmin,
+    ls: np.ndarray,
+    lm: np.ndarray,
+    codes: np.ndarray,
+    cfg: HPCConfig,
+    dc_carry: float = 0.0,
+) -> tuple[float, float]:
+    xh = dequantize_params(d, dmin, ls, lm, codes)
+    e = np.asarray(x, dtype=np.float64) - xh.astype(np.float64)
+    sse = float(np.dot(np.asarray(w, dtype=np.float64), e * e)) if w is not None else float(np.dot(e, e))
+    return sse + spectral_penalty(e, cfg, dc_carry), sse
+
+
+def refine_codes_hadamard(
+    x: np.ndarray,
+    w: np.ndarray | None,
+    d,
+    dmin,
+    ls: np.ndarray,
+    lm: np.ndarray,
+    codes: np.ndarray,
+    cfg: HPCConfig,
+    dc_carry: float = 0.0,
+) -> np.ndarray:
+    """Deterministic greedy descent of ``4*vesica^2 + (DC + carry)^2``.
+
+    Same acceptance rule as the C Phase-4 shaper, but the whole 256-element
+    search is evaluated as four vectorised candidate matrices per pass.
+    """
+    x = np.asarray(x, dtype=np.float64)
+    out = np.asarray(codes, dtype=np.uint8).copy()
+    ww = np.ones(QK_K, dtype=np.float64) if w is None else np.asarray(w, dtype=np.float64)
+    d_sub = np.repeat(np.float32(d) * np.asarray(ls, dtype=np.float32), 16).astype(np.float64)
+    m_sub = np.repeat(np.float32(dmin) * np.asarray(lm, dtype=np.float32), 16).astype(np.float64)
+    active = d_sub >= 1.0e-15
+    if not np.any(active):
+        return out
+
+    xh = d_sub * out.astype(np.float64) - m_sub
+    e_live = x - xh
+    sse_live = float(np.dot(ww, e_live * e_live))
+    sse_cap = sse_live * (1.0 + cfg.sse_budget)
+    v_live = e_live[: QK_K // 2] + e_live[QK_K // 2 :]
+    vesica = float(np.dot(v_live, v_live))
+    dc = float(np.sum(e_live)) + dc_carry
+    metric = 4.0 * vesica + dc * dc
+
+    idx = np.arange(QK_K)
+    pi = np.where(idx < QK_K // 2, idx, idx - QK_K // 2)
+    for _ in range(QK_K):
+        best_delta = 0.0
+        best_k = -1
+        best_q = 0
+        for q_try in range(4):
+            e_new = x - (d_sub * q_try - m_sub)
+            delta_e = e_new - e_live
+            sse_alt = sse_live + ww * (e_new * e_new - e_live * e_live)
+            mask = (sse_alt <= sse_cap) & active & (out != q_try)
+            if not np.any(mask):
+                continue
+            v_new = v_live[pi] + delta_e
+            vesica_alt = vesica - v_live[pi] * v_live[pi] + v_new * v_new
+            dc_alt = dc + delta_e
+            delta = metric - (4.0 * vesica_alt + dc_alt * dc_alt)
+            delta = np.where(mask, delta, -np.inf)
+            k = int(np.argmax(delta))
+            if delta[k] > best_delta:
+                best_delta = float(delta[k])
+                best_k = k
+                best_q = q_try
+        if best_k < 0:
+            break
+        old_e = e_live[best_k]
+        new_e = x[best_k] - (d_sub[best_k] * best_q - m_sub[best_k])
+        de = new_e - old_e
+        p = pi[best_k]
+        v_live[p] += de
+        vesica += v_live[p] * v_live[p] - (v_live[p] - de) * (v_live[p] - de)
+        dc += de
+        metric = 4.0 * vesica + dc * dc
+        sse_live += ww[best_k] * (new_e * new_e - old_e * old_e)
+        e_live[best_k] = new_e
+        out[best_k] = np.uint8(best_q)
+    return out
+
+
+def dc_nudge_codes(
+    x: np.ndarray,
+    d,
+    dmin,
+    ls: np.ndarray,
+    lm: np.ndarray,
+    codes: np.ndarray,
+    w: np.ndarray | None,
+    carry: float,
+    cfg: HPCConfig,
+    max_passes: int = 4,
+) -> np.ndarray:
+    """Spend a bounded SSE budget to cancel the rolling DC residual.
+
+    Equivalent in spirit to Phase 4.8 of the C engine: the reconstruction
+    target is always the true weight, only the codeword choice changes.
+    """
+    x = np.asarray(x, dtype=np.float64)
+    out = np.asarray(codes, dtype=np.uint8).copy()
+    ww = np.ones(QK_K, dtype=np.float64) if w is None else np.asarray(w, dtype=np.float64)
+    d_sub = np.repeat(np.float32(d) * np.asarray(ls, dtype=np.float32), 16).astype(np.float64)
+    m_sub = np.repeat(np.float32(dmin) * np.asarray(lm, dtype=np.float32), 16).astype(np.float64)
+    active = d_sub >= 1.0e-15
+    if not np.any(active):
+        return out
+    e = x - (d_sub * out.astype(np.float64) - m_sub)
+    sse = float(np.dot(ww, e * e))
+    cap = sse * (1.0 + cfg.sse_budget)
+    residual = float(np.sum(e)) + carry
+    for _ in range(max_passes):
+        best_improvement = 0.0
+        best_k = -1
+        best_q = 0
+        for q_try in range(4):
+            e_new = x - (d_sub * q_try - m_sub)
+            mask = active & (out != q_try)
+            sse_alt = sse + ww * (e_new * e_new - e * e)
+            mask &= sse_alt <= cap
+            if not np.any(mask):
+                continue
+            improvement = np.abs(residual) - np.abs(residual + (e_new - e))
+            improvement = np.where(mask, improvement, -np.inf)
+            k = int(np.argmax(improvement))
+            if improvement[k] > best_improvement:
+                best_improvement = float(improvement[k])
+                best_k = k
+                best_q = q_try
+        if best_k < 0:
+            break
+        old_e = e[best_k]
+        new_e = x[best_k] - (d_sub[best_k] * best_q - m_sub[best_k])
+        de = new_e - old_e
+        residual += de
+        sse += ww[best_k] * (new_e * new_e - old_e * old_e)
+        e[best_k] = new_e
+        out[best_k] = np.uint8(best_q)
+        if abs(residual) < 1.0e-12:
+            break
+    return out

viterbi.py

@@ -0,0 +1,164 @@
+"""36-state Viterbi selection over the Q2 symbolic lattice.
+
+This is a direct port of Phase 3 of the C engine.  Each graph block owns a
+``6 x 6`` coarse/fine lattice; candidate errors are reduced to a free-energy
+per bin, the Sieve marginals become a log-prior bonus, and an exact DP finds
+the minimum-cost state sequence with a Manhattan transition penalty.
+
+The forward pass is vectorised over the 36 states instead of the C engine's
+triple nested scalar loop.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .config import HPCConfig
+from .factored import CAND_TO_QUHIT, N_CAND_M, TOTAL_CANDIDATES
+
+N_STATES = 36
+
+
+def _median_upper(values: np.ndarray) -> float:
+    arr = np.asarray(values)
+    if arr.size == 0:
+        return 0.0
+    k = arr.size // 2
+    return float(np.partition(arr, k)[k])
+
+
+def robust_temperature(
+    candidate_errors: np.ndarray, cfg: HPCConfig
+) -> float:
+    n_blocks = candidate_errors.shape[0]
+    if n_blocks <= 0:
+        return 1.0e-10
+    step = max(1, n_blocks // cfg.temperature_sample_blocks)
+    medians = []
+    for b in range(0, n_blocks, step):
+        row = candidate_errors[b]
+        finite = row[np.isfinite(row)]
+        if finite.size == 0:
+            continue
+        medians.append(_median_upper(np.maximum(0.0, row - finite.min())))
+        if len(medians) >= cfg.temperature_sample_blocks:
+            break
+    t = _median_upper(np.asarray(medians))
+    return t if np.isfinite(t) and t > 1.0e-10 else 1.0e-10
+
+
+def _bin_geometry():
+    di = np.arange(TOTAL_CANDIDATES) // N_CAND_M
+    mi = np.arange(TOTAL_CANDIDATES) % N_CAND_M
+    qi_d = CAND_TO_QUHIT[di]
+    qi_m = CAND_TO_QUHIT[mi]
+    return qi_d, qi_m, qi_d * 6 + qi_m
+
+
+def viterbi_select(
+    candidate_errors: np.ndarray,
+    coarse_marg: np.ndarray,
+    fine_marg: np.ndarray,
+    temperature: float,
+    stride: int,
+    graph_blocks: int,
+    cfg: HPCConfig,
+) -> tuple[np.ndarray, np.ndarray]:
+    n_blocks = candidate_errors.shape[0]
+    qi_d, qi_m, bins = _bin_geometry()
+    bin_err = np.full((graph_blocks, N_STATES), np.inf, dtype=np.float64)
+    bin_cand = np.full((graph_blocks, N_STATES), -1, dtype=np.int64)
+    log_prior = np.zeros((graph_blocks, N_STATES), dtype=np.float64)
+    temp = max(float(temperature), 1.0e-12)
+
+    for gi in range(graph_blocks):
+        start = gi * stride
+        stop = min((gi + 1) * stride, n_blocks)
+        group = candidate_errors[start:stop].astype(np.float64)
+        flat = group.ravel()
+        finite = np.isfinite(flat)
+        group_min = np.full(N_STATES, np.inf)
+        if np.any(finite):
+            np.minimum.at(group_min, bins[finite], flat[finite])
+        counts = np.bincount(bins[finite], minlength=N_STATES)
+        sums = np.zeros(N_STATES, dtype=np.float64)
+        if np.any(finite):
+            shifted = flat[finite] - group_min[bins[finite]]
+            amps = np.exp(-shifted / (2.0 * temp))
+            sums = np.bincount(bins[finite], weights=amps, minlength=N_STATES)
+        valid = (counts > 0) & np.isfinite(group_min)
+        if np.any(valid):
+            mean_exp = np.maximum(sums[valid] / counts[valid], 1.0e-300)
+            bin_err[gi, valid] = group_min[valid] - 2.0 * temp * np.log(mean_exp)
+
+        rep = candidate_errors[start].astype(np.float64)
+        for s in range(N_STATES):
+            idx = np.where((bins == s) & np.isfinite(rep))[0]
+            if idx.size:
+                bin_cand[gi, s] = int(idx[int(np.argmin(rep[idx]))])
+
+        c_tot = float(coarse_marg[gi].sum())
+        f_tot = float(fine_marg[gi].sum())
+        pc = coarse_marg[gi] / c_tot if c_tot > 1.0e-30 else np.full(6, 1.0 / 6.0)
+        pf = fine_marg[gi] / f_tot if f_tot > 1.0e-30 else np.full(6, 1.0 / 6.0)
+        for s in range(N_STATES):
+            log_prior[gi, s] = np.log(pc[s // 6] + 1.0e-30) + np.log(pf[s % 6] + 1.0e-30)
+
+    finite_err = bin_err[np.isfinite(bin_err)]
+    scale_err = float(finite_err.mean()) if finite_err.size else 1.0e-10
+    scale_err = max(scale_err, 1.0e-20)
+
+    local = np.where(
+        np.isfinite(bin_err),
+        bin_err - cfg.viterbi_beta * scale_err * log_prior,
+        np.inf,
+    )
+
+    state_d = np.arange(N_STATES) // 6
+    state_m = np.arange(N_STATES) % 6
+    trans = np.empty((N_STATES, N_STATES), dtype=np.float64)
+    for sp in range(N_STATES):
+        trans[sp, :] = cfg.viterbi_alpha * scale_err * (
+            np.abs(state_d - state_d[sp]) + np.abs(state_m - state_m[sp])
+        )
+
+    dp = np.full((graph_blocks, N_STATES), np.inf, dtype=np.float64)
+    back = np.zeros((graph_blocks, N_STATES), dtype=np.int64)
+    dp[0] = local[0]
+    for gi in range(1, graph_blocks):
+        totals = dp[gi - 1][:, None] + trans
+        best = np.argmin(totals, axis=0)
+        best_val = totals[best, np.arange(N_STATES)]
+        dp[gi] = np.where(np.isfinite(best_val), best_val + local[gi], np.inf)
+        back[gi] = best
+
+    if not np.any(np.isfinite(dp[-1])):
+        path = np.zeros(graph_blocks, dtype=np.int64)
+    else:
+        path = np.zeros(graph_blocks, dtype=np.int64)
+        path[-1] = int(np.argmin(dp[-1]))
+        for gi in range(graph_blocks - 2, -1, -1):
+            path[gi] = back[gi + 1, path[gi + 1]]
+
+    best_candidate = np.zeros(n_blocks, dtype=np.int64)
+    for gi in range(graph_blocks):
+        s = int(path[gi])
+        qd, qm = s // 6, s % 6
+        rep_block = gi * stride
+        if bin_cand[gi, s] >= 0:
+            best_candidate[rep_block] = bin_cand[gi, s]
+        start = rep_block + 1
+        stop = min((gi + 1) * stride, n_blocks)
+        for b in range(start, stop):
+            mask = (qi_d == qd) & (qi_m == qm)
+            idx = np.where(mask)[0]
+            if idx.size:
+                best_candidate[b] = int(idx[int(np.argmin(candidate_errors[b, idx]))])
+
+    for b in range(n_blocks):
+        cur = candidate_errors[b, best_candidate[b]]
+        g = int(np.argmin(candidate_errors[b]))
+        if candidate_errors[b, g] < cur * cfg.greedy_override_ratio:
+            best_candidate[b] = g
+
+    return best_candidate, path