"""Place a picture on QR Scan Difficulty's scale from the five scanners' reads.

The page's scale comes from a two-parameter item-response fit of all eight scanners at seven frame sizes: each
(scanner, frame size) pair reads a code of difficulty b with chance 1 / (1 + exp(-a (theta - b))). Here the 35 curves
of the five scanners that run on Linux are frozen as fitted (model.json), and a picture's b is the fit's own estimate
for one code from its 35 reads: the most likely b under the fit's prior (normal, mean 0, sd 2.5). The log-posterior
has one peak (its slope falls as b grows), so bisection on the slope finds it; the fit's clipped Newton step can
bounce between two points when the readers are held fixed. A picture every pair read sits with the codes every
scanner read at every size (the left pile); one no pair read sits with the codes no scanner read (the right pile).
Otherwise d = 100 (b - lo) / (hi - lo), clipped to 0..100.
"""

import json
from pathlib import Path

import numpy as np

MODEL = json.loads(Path(__file__).with_name("model.json").read_text())
PILE_LO, PILE_HI = -6, 106


def map_b(y: np.ndarray, theta: np.ndarray, a: np.ndarray, prior_b: float) -> float:
    """The most likely b given reads y by readers (theta, a): where the log-posterior's slope, the fit's gradient,
    is zero."""
    def slope(b: float) -> float:
        p = 1.0 / (1.0 + np.exp(-np.clip(a * (theta - b), -30, 30)))
        return float(-(a * (y - p)).sum() - b / prior_b ** 2)

    lo, hi = -30.0, 30.0   # slope > 0 at lo and < 0 at hi for any reads
    for _ in range(80):
        mid = (lo + hi) / 2
        if slope(mid) > 0:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


def place(reads: dict[str, list[int]], frames: list[int]) -> dict:
    """{"d": difficulty, "pile": "lo" | "hi" | None} from {scanner: [read at each frame size]}."""
    order = [frames.index(f) for f in MODEL["frames"]]
    y = np.array([reads[k][i] for k in MODEL["decoders"] for i in order], dtype=float)
    if y.all():
        return {"d": PILE_LO, "pile": "lo"}
    if not y.any():
        return {"d": PILE_HI, "pile": "hi"}
    b = map_b(y, np.array(MODEL["theta"]).ravel(), np.array(MODEL["a"]).ravel(), MODEL["prior_b"])
    d = 100 * (b - MODEL["lo"]) / (MODEL["hi"] - MODEL["lo"])
    return {"d": round(min(100.0, max(0.0, d)), 1), "pile": None}