python/dpt/validation/transmission.py

1"""Independent exact-input Decimal oracles for the CUDA transmission contract.23These scalar routines deliberately have no Warp dependency and no production4cutover. Inputs are quantised to binary32, reconstructed exactly by integer ratio,5and evaluated at the requested Decimal precision (100 digits by default).6"""78import math9import struct10from collections.abc import Callable, Iterable11from dataclasses import dataclass12from decimal import Decimal, localcontext131415def binary32(value: float) -> float:16    """Round a Python number to binary32; finite overflow is an explicit error."""17    try:18        result = struct.unpack("!f", struct.pack("!f", value))[0]19    except OverflowError as exc:20        raise ValueError("input exceeds finite binary32 range") from exc21    if not math.isfinite(result):22        raise ValueError("input must be finite binary32")23    return result242526def _from_bits(bits: int) -> float:27    return struct.unpack("!f", struct.pack("!I", bits))[0]282930def _bits(value: float) -> int:31    return struct.unpack("!I", struct.pack("!f", value))[0]323334def exact_input(value: float, *, nonnegative: bool = False) -> Decimal:35    """Recover the exact stored binary32 value, preserving signed zero."""36    stored = binary32(value)37    if nonnegative and stored < 0:38        raise ValueError("input must be non-negative")39    numerator, denominator = stored.as_integer_ratio()40    # A binary32 denominator is at most 2**149: 200 digits make this exact.41    with localcontext() as context:42        context.prec = 20043        result = Decimal(numerator) / Decimal(denominator)44    return result.copy_negate() if stored == 0 and math.copysign(1, stored) < 0 else result454647MIN_SUBNORMAL = Decimal.from_float(_from_bits(1))48MIN_NORMAL = Decimal.from_float(_from_bits(0x00800000))49MAX_FINITE = Decimal.from_float(_from_bits(0x7F7FFFFF))50with localcontext() as _context:51    _context.prec = 20052    U32 = Decimal(2) ** -2453    U64 = Decimal(2) ** -5354    OVERFLOW_MIDPOINT = MAX_FINITE + Decimal(2) ** 103555657def round_binary32(value: Decimal) -> float:58    """Round Decimal directly to nearest binary32, ties to even.5960    Binary search compares exact neighbours; converting through binary64 could61    incorrectly move a Decimal value onto a binary32 midpoint. Overflow follows62    round-to-nearest semantics and returns a signed infinity.63    """64    if not value.is_finite():65        raise ValueError("reference value must be finite")66    negative = value.is_signed()67    magnitude = value.copy_abs()68    if magnitude >= OVERFLOW_MIDPOINT:69        return -math.inf if negative else math.inf70    low, high = 0, 0x7F7FFFFF71    while low < high:72        middle = (low + high + 1) // 273        if Decimal.from_float(_from_bits(middle)) <= magnitude:74            low = middle75        else:76            high = middle - 177    chosen = low78    if low < 0x7F7FFFFF:79        # Exact endpoints need up to 150 decimal places, independent of caller context.80        with localcontext() as context:81            context.prec = max(200, len(value.as_tuple().digits) + 160)82            midpoint = (83                Decimal.from_float(_from_bits(low)) + Decimal.from_float(_from_bits(low + 1))84            ) / 285        if magnitude > midpoint or (magnitude == midpoint and low % 2):86            chosen += 187    return _from_bits(chosen | (0x80000000 if negative else 0))888990def _exponential(depth: Decimal) -> tuple[Decimal, Decimal]:91    """Return exp(-L), or zero with an explicit absolute bound for enormous L.9293    For L >= 1024, exp(-L) <= exp(-1024) < 2**-1477. Even two94    maximum binary32 factors and 2**31 contributions give < 2**-1190,95    far below half the smallest binary32 subnormal (2**-150). This96    bound covers every exponential weighted term in the declared contract.97    """98    if depth >= 1024:99        return Decimal(0), Decimal(2) ** -1477100    return (-depth).exp(), Decimal(0)101102103@dataclass(frozen=True)104class ForwardReference:105    """Mathematical values and the exponential tail bound, before storage rounding."""106107    T: Decimal108    counts: Decimal109    log_T: Decimal  # noqa: N815 - match the documented mathematical output name.110    removed: Decimal111    tail_bound: Decimal112113    def rounded(self) -> dict[str, float]:114        return {115            name: round_binary32(getattr(self, name))116            for name in ("T", "counts", "log_T", "removed")117        }118119120def forward_reference(depth: float, n0: float = 1.0, *, precision: int = 100) -> ForwardReference:121    """Evaluate four outputs from exact stored inputs, independently of device code."""122    optical_depth, beam = exact_input(depth, nonnegative=True), exact_input(n0, nonnegative=True)123    with localcontext() as context:124        context.prec = precision125        exponential, bound = _exponential(optical_depth)126        return ForwardReference(127            exponential,128            (beam * exponential).copy_abs(),129            optical_depth.copy_negate(),130            (Decimal(1) - exponential).copy_abs(),131            bound,132        )133134135def inverse_reference(delta: float, *, precision: int = 100) -> Decimal:136    """Reference -ln(1-delta); the supplied decrement is the input representation."""137    decrement = exact_input(delta, nonnegative=True)138    if decrement >= 1:139        raise ValueError("decrement must be less than one")140    with localcontext() as context:141        context.prec = precision142        return (-(Decimal(1) - decrement).ln()).copy_abs()143144145@dataclass(frozen=True)146class VJPReference:147    """Analytic weighted derivatives and the frozen per-pixel absolute budgets."""148149    grad_L: Decimal  # noqa: N815 - L denotes optical depth throughout the API.150    grad_n0: Decimal151    budget_L: Decimal  # noqa: N815 - pair the budget with grad_L.152    budget_n0: Decimal153154155def vjp_reference(156    depth: float,157    n0: float = 1.0,158    *,159    seed_T: float = 0.0,  # noqa: N803 - mirror the production cotangent name.160    seed_counts: float = 0.0,161    seed_log_T: float = 0.0,  # noqa: N803 - mirror the production cotangent name.162    seed_removed: float = 0.0,163    precision: int = 100,164) -> VJPReference:165    """Differentiate each mathematical output independently and contract its seed."""166    optical_depth, beam = exact_input(depth, nonnegative=True), exact_input(n0, nonnegative=True)167    seeds = [exact_input(value) for value in (seed_T, seed_counts, seed_log_T, seed_removed)]168    with localcontext() as context:169        context.prec = precision170        exponential, _ = _exponential(optical_depth)171        derivatives = (-exponential, -beam * exponential, Decimal(-1), exponential)172        terms = [seed * derivative for seed, derivative in zip(seeds, derivatives, strict=True)]173        beam_term = seeds[1] * exponential174        return VJPReference(175            sum(terms, Decimal(0)),176            beam_term,177            32 * U32 * sum((abs(term) for term in terms), Decimal(0)) + 2 * MIN_SUBNORMAL,178            32 * U32 * abs(beam_term) + 2 * MIN_SUBNORMAL,179        )180181182def inverse_vjp_reference(183    delta: float, seed: float, *, precision: int = 100184) -> tuple[Decimal, Decimal]:185    """Return the independent inverse derivative and magnitude-scaled error budget."""186    inverse_reference(delta, precision=precision)  # Validate the physical domain.187    decrement, weight = exact_input(delta), exact_input(seed)188    with localcontext() as context:189        context.prec = precision190        gradient = weight / (1 - decrement)191        return gradient, 32 * U32 * abs(gradient) + 2 * MIN_SUBNORMAL192193194def scalar_beam_reference(195    depths: Iterable[float],196    seeds: Iterable[float],197    *,198    addition_depth: int,199    precision: int = 100,200) -> tuple[Decimal, Decimal]:201    """Sum beam contributions with a bound using the actual reduction addition depth.202203    The supplied depth counts additions along the longest path, including within204    blocks. Kernel-launch count is not an acceptable substitute.205    """206    if addition_depth < 0:207        raise ValueError("addition_depth must be non-negative")208    with localcontext() as context:209        context.prec = precision210        references = [211            vjp_reference(depth, seed_counts=seed, precision=precision)212            for depth, seed in zip(depths, seeds, strict=True)213        ]214        product = addition_depth * U64215        if product >= 1:216            raise ValueError("addition depth exceeds the gamma-bound domain")217        total = sum((record.grad_n0 for record in references), Decimal(0))218        magnitude = sum((abs(record.grad_n0) for record in references), Decimal(0))219        element_budget = sum((record.budget_n0 for record in references), Decimal(0))220        final_rounding = U32 * abs(total) + MIN_SUBNORMAL / 2221        return total, element_budget + product / (1 - product) * magnitude + final_rounding222223224def ulp_distance(actual: float, expected: float) -> int:225    """Count representable binary32 steps; opposite signed zeros have distance zero."""226227    def ordered(value: float) -> int:228        bits = _bits(binary32(value))229        return 0x80000000 - (bits & 0x7FFFFFFF) if bits >> 31 else 0x80000000 + bits230231    return abs(ordered(actual) - ordered(expected))232233234@dataclass(frozen=True)235class ValueErrorRecord:236    absolute_error: Decimal237    ulps: int238    allowed_ulps: int239    regime: str240    passed: bool241242243def value_error(actual: float, expected: Decimal, *, counts: bool = False) -> ValueErrorRecord:244    """Apply frozen forward budgets, including one-ULP subnormal/zero checks."""245    target = round_binary32(expected)246    actual_exact = exact_input(actual)247    distance = ulp_distance(actual, target)248    subnormal = abs(expected) < MIN_NORMAL249    allowed = 1 if subnormal else (8 if counts else 4)250    with localcontext() as context:251        context.prec = 200252        absolute = abs(actual_exact - expected)253    return ValueErrorRecord(254        absolute, distance, allowed, "subnormal" if subnormal else "normal", distance <= allowed255    )256257258@dataclass(frozen=True)259class AnalyticCase:260    name: str261    optical_depth: Decimal262    definition: str263264265# region book:transmission-reference-cases266def analytic_cases() -> tuple[AnalyticCase, ...]:267    """Closed-form path cases, without claiming to validate a path integrator.268269    All coefficients and lengths below are exact dyadic mathematical stress inputs.270    The mm and cm expressions describe the same optical depth in different units.271    """272    coefficient, distance = Decimal("0.125"), Decimal(4)273    first_length, second_length = Decimal(1), Decimal(3)274    intercept, slope = Decimal("0.125"), Decimal("0.0625")275    return (276        AnalyticCase("empty", Decimal(0), "empty path"),277        AnalyticCase("zero-attenuation", Decimal(0), "mu = 0"),278        AnalyticCase("homogeneous", coefficient * distance, "mu*d"),279        AnalyticCase(280            "split-homogeneous",281            coefficient * first_length + coefficient * second_length,282            "mu*d1 + mu*d2",283        ),284        AnalyticCase("layered", Decimal("0.25") * 2 + Decimal("0.5") * 3, "mu1*d1 + mu2*d2"),285        AnalyticCase("millimetres", Decimal("0.125") * 4, "0.125/mm * 4 mm"),286        AnalyticCase("centimetres", Decimal("1.25") * Decimal("0.4"), "1.25/cm * 0.4 cm"),287        AnalyticCase(288            "added-segment",289            coefficient * distance + Decimal("0.25"),290            "mu*d + non-negative optical depth",291        ),292        AnalyticCase(293            "linear-coefficient",294            intercept * distance + slope * distance * distance / 2,295            "integral_0^d (a+b*s) ds",296        ),297    )298299300# endregion book:transmission-reference-cases301302303@dataclass(frozen=True)304class DirectionalRecord:305    step: float306    derivative: float307    absolute_error: float308    stencil: str309310311# region book:transmission-directional-check312def directional_check(313    evaluate: Callable[[float], float],314    anchor: float,315    analytic_derivative: float,316    *,317    exponents: Iterable[int] = range(2, 25),318) -> tuple[DirectionalRecord, ...]:319    """Sweep a real operator callback over admissible, distinct binary32 inputs.320321    Interior anchors use centred differences; zero uses a second-order one-sided322    stencil. No convergence claim is inferred from the smallest step: callers retain323    all records to expose truncation, agreement and storage-rounding regimes.324    """325    centre = binary32(anchor)326    if centre < 0:327        raise ValueError("anchor must be non-negative")328    records: list[DirectionalRecord] = []329    for exponent in exponents:330        step = 2.0**-exponent331        right = binary32(centre + step)332        if right == centre:333            continue334        if centre == 0:335            second = binary32(2 * step)336            if second <= right:337                continue338            derivative = (-3 * evaluate(centre) + 4 * evaluate(right) - evaluate(second)) / (339                2 * right340            )341            stencil = "one-sided-second-order"342        else:343            if centre - step < 0:344                continue345            left = binary32(centre - step)346            if left == centre or centre - left != right - centre:347                continue348            derivative = (evaluate(right) - evaluate(left)) / (right - left)349            stencil = "centred"350        records.append(351            DirectionalRecord(352                right - centre, derivative, abs(derivative - analytic_derivative), stencil353            )354        )355    return tuple(records)356357358# endregion book:transmission-directional-check359