python/dpt/transport/model.py

1"""Host-only contracts for fixed-grid photon transport estimators.23Supplied linear coefficients in inverse millimetres define absorption plus either4isotropic elastic or free-electron Klein-Nishina/Compton scattering. The latter5changes energy and uses an explicit coefficient energy table; it does not include6bound-electron, polarisation or coherent-scattering corrections. The material7grid is piecewise constant, unlike the deterministic projector's trilinear field.8"""910from __future__ import annotations1112import math13from dataclasses import dataclass14from enum import IntEnum15from itertools import pairwise16from typing import Any, Literal, cast1718from dpt.contracts import ContractError, finite_scalar192021class TransportError(ContractError):22    """Compatibility subtype for transport physical, storage and execution contracts."""232425class IncompleteHistoryError(RuntimeError):26    """A resource guard stopped a history; the entire estimate is invalid."""272829class HistoryStatus(IntEnum):30    """Each launched history writes exactly one terminal status."""3132    ESCAPED = 033    ABSORBED = 134    EVENT_BUDGET = 235    CROSSING_BUDGET = 336    NUMERICAL_FAILURE = 437    ENERGY_SUPPORT = 538    ANGLE_BUDGET = 6394041def _positive_integer(value: int, name: str) -> None:42    if type(value) is not int or not 0 < value < 2**31:43        raise TransportError(f"{name} must be a positive signed 32-bit integer")444546def _physical(value: Any, name: str, *, positive: bool = False, nonnegative: bool = False) -> float:47    try:48        scalar = finite_scalar(value, name, minimum=0.0 if positive or nonnegative else None)49    except ContractError as error:50        raise TransportError(str(error)) from error51    if positive and scalar == 0:52        raise TransportError(f"{name} must be positive")53    return scalar545556def _freeze_fields(instance: Any, names: tuple[str, ...]) -> None:57    for name in names:58        object.__setattr__(instance, name, tuple(getattr(instance, name)))596061@dataclass(frozen=True, slots=True)62class MaterialGrid:63    """Axis-aligned cell faces; shape is (nz, ny, nx), device indices x fastest."""6465    origin_mm: tuple[float, float, float]66    spacing_mm: tuple[float, float, float]67    shape: tuple[int, int, int]68    materials: int6970    def __post_init__(self) -> None:71        _freeze_fields(self, ("origin_mm", "spacing_mm", "shape"))72        if len(self.origin_mm) != 3 or len(self.spacing_mm) != 3 or len(self.shape) != 3:73            raise TransportError("origin, spacing and shape must each have three components")74        _positive_integer(self.materials, "materials")75        for count, step, origin in zip(76            self.shape_xyz, self.spacing_mm, self.origin_mm, strict=True77        ):78            _positive_integer(count, "grid extent")79            _physical(step, "cell spacing", positive=True)80            _physical(origin, "cell origin")81            endpoint = origin + count * step82            if (83                not math.isfinite(endpoint)84                or origin + step == origin85                or endpoint - step == endpoint86            ):87                raise TransportError("cell faces must remain distinguishable in binary64")88        if math.prod(self.shape) >= 2**31:89            raise TransportError("the flattened material grid exceeds signed 32-bit indexing")9091    @property92    def shape_xyz(self) -> tuple[int, int, int]:93        """Geometric extents for device vectors; storage shape remains (z,y,x)."""94        return self.shape[2], self.shape[1], self.shape[0]9596    @property97    def upper_mm(self) -> tuple[float, float, float]:98        return cast(99            tuple[float, float, float],100            tuple(101                start + count * step102                for start, count, step in zip(103                    self.origin_mm, self.shape_xyz, self.spacing_mm, strict=True104                )105            ),106        )107108109@dataclass(frozen=True, slots=True)110class PlanarDetector:111    """An outward-facing z plane, pixel lower faces (x,y), and count score.112113    Pixel support is half-open. The detector must be above the entire material114    grid; an escaped straight ray can therefore hit it at most once. A recorded115    score includes both uncollided and scattered photons. Do not add a separate116    deterministic primary image to this estimate.117    """118119    lower_xy_mm: tuple[float, float]120    spacing_xy_mm: tuple[float, float]121    shape: tuple[int, int]  # width, height122    z_mm: float123124    def __post_init__(self) -> None:125        _freeze_fields(self, ("lower_xy_mm", "spacing_xy_mm", "shape"))126        _physical(self.z_mm, "detector z")127        if any(len(value) != 2 for value in (self.lower_xy_mm, self.spacing_xy_mm, self.shape)):128            raise TransportError("detector xy fields must have two components")129        for count, step, lower in zip(130            self.shape, self.spacing_xy_mm, self.lower_xy_mm, strict=True131        ):132            _positive_integer(count, "detector extent")133            _physical(step, "pixel spacing", positive=True)134            _physical(lower, "detector origin")135            end = lower + count * step136            if not math.isfinite(end):137                raise TransportError("detector endpoint must be finite")138            if lower + step == lower or end - step == end:139                raise TransportError("detector pixels must be distinguishable in binary64")140        if self.pixels >= 2**31:141            raise TransportError("the flattened detector exceeds signed 32-bit indexing")142143    @property144    def pixels(self) -> int:145        return math.prod(self.shape)146147148# region book:transport-model-contract149@dataclass(frozen=True, slots=True)150class TransportSpec:151    """A complete physical scope and finite launch-resource policy.152153    Source rays and nonnegative importance weights are supplied by the caller.154    Their sampling law must be fixed with respect to active material parameters;155    the library does not silently invent an emission distribution. Provenance156    identifies the user's coefficient source, not a fabricated material asset.157    """158159    grid: MaterialGrid160    detector: PlanarDetector161    energy_kev: float162    coefficient_provenance: str163    coefficient_energies_kev: tuple[float, ...] = ()164    scattering_law: Literal["isotropic-elastic", "free-electron-compton"] = "isotropic-elastic"165    scoring: Literal["photon-count", "energy-kev"] = "photon-count"166    max_angle_trials: int = 4096167    max_events: int = 4096168    max_crossings: int = 65536169    block_dim: Literal[64, 128, 256] = 128170    estimator: Literal["analogue", "continuous-absorption"] = "analogue"171172    def __post_init__(self) -> None:173        _freeze_fields(self, ("coefficient_energies_kev",))174        _physical(self.energy_kev, "source energy", positive=True)175        for node in self.coefficient_energies_kev:176            _physical(node, "coefficient energy", nonnegative=True)177        if type(self.coefficient_provenance) is not str or not self.coefficient_provenance.strip():178            raise TransportError("the supplied coefficients require a provenance identifier")179        if self.detector.z_mm <= self.grid.upper_mm[2]:180            raise TransportError("the detector plane must lie above the entire material grid")181        nodes = self.energy_nodes182        if any(right <= left for left, right in pairwise(nodes)):183            raise TransportError("coefficient energies must be strictly increasing")184        if not nodes[0] <= self.energy_kev <= nodes[-1]:185            raise TransportError("coefficient energies must contain the source energy")186        if self.scattering_law not in ("isotropic-elastic", "free-electron-compton"):187            raise TransportError("unsupported scattering law")188        if self.scattering_law == "free-electron-compton" and len(nodes) < 2:189            raise TransportError("Compton scattering needs energy-dependent coefficient tables")190        if self.scoring not in ("photon-count", "energy-kev"):191            raise TransportError("unsupported detector score")192        if self.estimator not in ("analogue", "continuous-absorption"):193            raise TransportError("unsupported transport estimator")194        _positive_integer(self.max_angle_trials, "max_angle_trials")195        _positive_integer(self.max_events, "max_events")196        _positive_integer(self.max_crossings, "max_crossings")197        if self.block_dim not in (64, 128, 256):198            raise TransportError("block_dim must be 64, 128 or 256")199200    @property201    def energy_nodes(self) -> tuple[float, ...]:202        """An omitted grid declares coefficients at the one elastic source energy."""203        return self.coefficient_energies_kev or (self.energy_kev,)204205206# endregion book:transport-model-contract207