python/dpt/materials.py

1"""Explicit, CPU-safe ingestion of primary attenuation coefficients.23There are no bundled physical tables. A caller supplies the values, identifies4where they came from and records the rights under which they may be used.5Interpolation is confined to the supplied support; repeated energies encode the6two one-sided values at an absorption edge rather than an interval to smooth.7"""89from __future__ import annotations1011# Public Python boundaries validate runtime inputs; workspace scratch stays module-owned.12# pyright: reportPrivateUsage=false, reportUnnecessaryIsInstance=false13import math14import re15from bisect import bisect_left, bisect_right16from dataclasses import dataclass17from itertools import pairwise18from typing import Literal1920from dpt.contracts import ContractError, finite_scalar, finite_tuple212223@dataclass(frozen=True, slots=True)24class Provenance:25    """Identity of the actual input, distinct from a bibliographic citation.2627    ``sha256`` identifies the original source bytes, before unit conversion.28    Analytic fixtures must say so in ``description`` and cite their defining29    formula in ``source``; a real material name is not fixture provenance.30    """3132    source: str33    sha256: str34    rights: str35    description: str3637    def __post_init__(self) -> None:38        for name in ("source", "rights", "description"):39            value = getattr(self, name)40            if not isinstance(value, str) or not value.strip():41                raise ContractError(f"provenance {name} must be non-empty")42        if not isinstance(self.sha256, str) or re.fullmatch(r"[0-9a-f]{64}", self.sha256) is None:43            raise ContractError("provenance sha256 must contain 64 lowercase hexadecimal digits")444546def validate_energy_grid(energies_kev: tuple[float, ...]) -> None:47    """Validate distinct, positive keV nodes shared by spectral inputs and response."""48    finite_tuple(energies_kev, "energies_kev", positive=True)49    if any(a >= b for a, b in pairwise(energies_kev)):50        raise ContractError("energy nodes must be strictly increasing")515253def _log_ratio(high: float, low: float) -> float:54    """Retain adjacent large energies without overflowing a widely separated ratio."""55    relative = (high - low) / low56    return math.log1p(relative) if math.isfinite(relative) else math.log(high) - math.log(low)575859# region book:material-table-contract60@dataclass(frozen=True, slots=True)61class MaterialTable:62    """Total primary attenuation at one declared material composition/density.6364    Arrays are immutable host metadata, never an alternative production backend.65    Energies are keV. Duplicate adjacent energies encode below/above edge values,66    in that order; at most two entries may share an energy. The density attached67    to a mass table is in g/cm³. Linear coefficients already include density.68    """6970    name: str71    energies_kev: tuple[float, ...]72    coefficients: tuple[float, ...]73    units: Literal["mm^-1", "cm^-1", "cm^2/g"]74    provenance: Provenance75    density_g_cm3: float | None = None76    interpolation: Literal["linear", "log-log"] = "log-log"7778    def __post_init__(self) -> None:79        if not self.name.strip():80            raise ContractError("a material table needs an unambiguous material/composition name")81        if not isinstance(self.provenance, Provenance):82            raise ContractError("material-table provenance is required")83        finite_tuple(self.energies_kev, "energies_kev", positive=True)84        finite_tuple(self.coefficients, "coefficients")85        if len(self.energies_kev) != len(self.coefficients):86            raise ContractError("energies and coefficients must have equal lengths")87        if self.units not in ("mm^-1", "cm^-1", "cm^2/g"):88            raise ContractError("unsupported coefficient units")89        if self.interpolation not in ("linear", "log-log"):90            raise ContractError("interpolation must be linear or log-log")91        if self.interpolation == "log-log" and any(c <= 0 for c in self.coefficients):92            raise ContractError("log-log interpolation requires strictly positive coefficients")93        for index, energy in enumerate(self.energies_kev[1:], 1):94            if energy < self.energies_kev[index - 1]:95                raise ContractError("energy support must be non-decreasing")96            if index > 1 and energy == self.energies_kev[index - 2]:97                raise ContractError("an absorption edge has exactly two one-sided entries")98        if self.units == "cm^2/g":99            density = self.density_g_cm3100            density = finite_scalar(density, "reference density in g/cm³", minimum=0)101            if density <= 0:102                raise ContractError("reference density must be strictly positive")103        elif self.density_g_cm3 is not None:104            raise ContractError(105                "linear coefficients already include density; do not apply it twice"106            )107108    @property109    def linear_mm_inverse(self) -> tuple[float, ...]:110        """Convert source values once, without inventing a different mixture model."""111        scale = 1.0112        if self.units == "cm^-1":113            scale = 0.1114        elif self.units == "cm^2/g":115            assert self.density_g_cm3 is not None116            scale = self.density_g_cm3 / 10.0117        converted = tuple(c * scale for c in self.coefficients)118        if any(not math.isfinite(c) for c in converted):119            raise ContractError("unit conversion overflowed; coefficient table is unusable")120        if any(121            source > 0 and target == 0122            for source, target in zip(self.coefficients, converted, strict=True)123        ):124            raise ContractError("unit conversion underflowed a positive attenuation coefficient")125        return converted126127    def at_energies(128        self,129        energies_kev: tuple[float, ...],130        *,131        edge_side: Literal["below", "above"],132    ) -> tuple[float, ...]:133        """Interpolate in mm⁻¹, requiring the one-sided convention at exact edges."""134        if edge_side not in ("below", "above"):135            raise ContractError("edge_side must be explicitly below or above")136        finite_tuple(energies_kev, "query energies", positive=True)137        grid = self.energies_kev138        coefficients = self.linear_mm_inverse139        result: list[float] = []140        for energy in energies_kev:141            if energy < grid[0] or energy > grid[-1]:142                raise ContractError(f"{self.name}: {energy} keV lies outside supplied support")143            left, right = bisect_left(grid, energy), bisect_right(grid, energy)144            if left != right:145                result.append(coefficients[left if edge_side == "below" else right - 1])146                continue147            low, high = left - 1, left148            if self.interpolation == "log-log":149                fraction = _log_ratio(energy, grid[low]) / _log_ratio(grid[high], grid[low])150                value = math.exp(151                    (1.0 - fraction) * math.log(coefficients[low])152                    + fraction * math.log(coefficients[high])153                )154            else:155                fraction = (energy - grid[low]) / (grid[high] - grid[low])156                value = (1.0 - fraction) * coefficients[low] + fraction * coefficients[high]157            result.append(value)158        return tuple(result)159160161# endregion book:material-table-contract162163164@dataclass(frozen=True, slots=True)165class MaterialBasis:166    """Fixed-density, non-negative dimensionless basis fields integrated in mm.167168    ``mixing_assumption`` states whether fields are volume fractions, material169    indicators or relative concentrations. There is no implicit HU conversion.170    Spatially varying mass density must be included in the declared basis field.171    """172173    tables: tuple[MaterialTable, ...]174    mixing_assumption: str175176    def __post_init__(self) -> None:177        if not isinstance(self.tables, tuple) or not self.tables:178            raise ContractError("a material basis requires an immutable tuple of tables")179        if any(not isinstance(table, MaterialTable) for table in self.tables):180            raise ContractError("material basis entries must be MaterialTable instances")181        if len({table.name for table in self.tables}) != len(self.tables):182            raise ContractError("material basis names must be unique")183        if not self.mixing_assumption.strip():184            raise ContractError("the material-basis mixing assumption must be recorded")185186    def coefficients_at(187        self, energies_kev: tuple[float, ...], *, edge_side: Literal["below", "above"]188    ) -> tuple[float, ...]:189        """Return material-major (material, energy) coefficients for explicit upload."""190        return tuple(191            value192            for table in self.tables193            for value in table.at_energies(energies_kev, edge_side=edge_side)194        )195