diff --git a/.github/workflows/fp-arena-ci.yml b/.github/workflows/fp-arena-ci.yml index 2e0bf90..54da94c 100644 --- a/.github/workflows/fp-arena-ci.yml +++ b/.github/workflows/fp-arena-ci.yml @@ -53,7 +53,7 @@ jobs: # FP-Arena without deps so it uses the DaCe just installed (any main), # rather than re-pulling pinned main from git. pip install --no-deps ./fp-arena - pip install pytest + pip install pytest scipy tqdm - name: Run FP-Arena tests run: | diff --git a/.gitignore b/.gitignore index 39c38a9..48ecb21 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,7 @@ out.sdfg .idea/ *.swp .DS_Store + + +#Experiment results database +*.db \ No newline at end of file diff --git a/corpus/__init__.py b/corpus/__init__.py new file mode 100644 index 0000000..07a7e98 --- /dev/null +++ b/corpus/__init__.py @@ -0,0 +1,2 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +"""A corpus of benchmark kernels written in the DaCe Python frontend.""" diff --git a/examples/heat3d.py b/corpus/heat3d.py similarity index 59% rename from examples/heat3d.py rename to corpus/heat3d.py index e478da0..94748dd 100644 --- a/examples/heat3d.py +++ b/corpus/heat3d.py @@ -1,20 +1,18 @@ # Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. -"""Run the heat3d stencil in reduced precision.""" +""" +The heat3d stencil kernel (Polybench), a 3D Jacobi-style heat-diffusion step. +""" import dace as dc -import numpy as np - -from fp_arena.transformations.change_and_propagate_fp_types import ( - change_and_propagate_fp_types, -) +# Default problem size: an N x N x N grid advanced over TSTEPS time steps. GRID_N, TSTEPS = 40, 20 N = dc.symbol("N", dtype=dc.int64) @dc.program -def kernel(TSTEPS: dc.int64, A: dc.float64[N, N, N], B: dc.float64[N, N, N]): +def heat3d_kernel(TSTEPS: dc.int64, A: dc.float64[N, N, N], B: dc.float64[N, N, N]): for t in range(1, TSTEPS): B[1:-1, 1:-1, 1:-1] = ( 0.125 * (A[2:, 1:-1, 1:-1] - 2.0 * A[1:-1, 1:-1, 1:-1] + A[:-2, 1:-1, 1:-1]) @@ -32,22 +30,3 @@ def kernel(TSTEPS: dc.int64, A: dc.float64[N, N, N], B: dc.float64[N, N, N]): * (B[1:-1, 1:-1, 2:] - 2.0 * B[1:-1, 1:-1, 1:-1] + B[1:-1, 1:-1, :-2]) + B[1:-1, 1:-1, 1:-1] ) - - -def main(): - sdfg = kernel.to_sdfg(simplify=True) - - change_and_propagate_fp_types( - sdfg, - {"A": dc.float16, "B": dc.float16}, - constant_type=dc.float16, - ) - - rng = np.random.default_rng(0) - A = rng.uniform(0, 100, (GRID_N,) * 3) - B = rng.uniform(0, 100, (GRID_N,) * 3) - sdfg(TSTEPS=TSTEPS, A=A, B=B, N=GRID_N) - - -if __name__ == "__main__": - main() diff --git a/examples/experiment_error.py b/examples/experiment_error.py new file mode 100644 index 0000000..eb61658 --- /dev/null +++ b/examples/experiment_error.py @@ -0,0 +1,55 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +"""Measure the heat3d stencil's error in reduced precision against an fp64 reference.""" + +from scipy import stats + +from corpus.heat3d import GRID_N, TSTEPS, heat3d_kernel +from fp_arena.experiment import ( + ErrorAnalysisConfig, + ExperimentConfig, + ResultStore, + run_error, +) + +N_SAMPLES = 3 + + +def main(): + experiment = ExperimentConfig( + name="heat3d", + program=heat3d_kernel.to_sdfg(simplify=True), + symbols={"N": GRID_N}, + scalar_args={"TSTEPS": TSTEPS}, + target="cpu", + inputs={ + "A": stats.uniform(0, 100), + "B": stats.uniform(0, 100), + }, + ) + + store = ResultStore("heat3d_results.db") + results = run_error( + ErrorAnalysisConfig( + experiment, + precisions=[{"A": "fp32", "B": "fp32"}, {"A": "fp16", "B": "fp16"}], + reference="fp64", + n_samples=N_SAMPLES, + ), + store=store, + ) + + for r in results: + label = " ".join(f"{k}={v}" for k, v in r.precision.items()) + for arr in sorted(r.errors): + e = r.errors[arr] + print( + f"{label:>16} {arr}: rel_mean {e.rel_mean:.3e} " + f"rel_max {e.rel_max:.3e} linf {e.linf:.3e}" + ) + + print(f"{len(store.query(experiment='heat3d', kind='error'))} rows in the database") + store.close() + + +if __name__ == "__main__": + main() diff --git a/examples/experiment_performance.py b/examples/experiment_performance.py new file mode 100644 index 0000000..6a80f62 --- /dev/null +++ b/examples/experiment_performance.py @@ -0,0 +1,54 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +"""Time the heat3d stencil across precision points with the experiment framework.""" + +import statistics + +from scipy import stats + +from corpus.heat3d import GRID_N, TSTEPS, heat3d_kernel +from fp_arena.experiment import ( + ExperimentConfig, + PerformanceAnalysisConfig, + ResultStore, + run_performance, +) + + +def main(): + experiment = ExperimentConfig( + name="heat3d", + program=heat3d_kernel.to_sdfg(simplify=True), + symbols={"N": GRID_N}, + scalar_args={"TSTEPS": TSTEPS}, + target="cpu", + inputs={ + "A": stats.uniform(0, 100), + "B": stats.uniform(0, 100), + }, + ) + + store = ResultStore("heat3d_results.db") + results = run_performance( + PerformanceAnalysisConfig( + experiment, + precisions=[ + {}, # the unmodified fp64 program + {"A": "fp32", "B": "fp32"}, + {"A": "fp16", "B": "fp16"}, + ], + n_warmup=1, + n_reps=5, + ), + store=store, + ) + + for r in results: + label = " ".join(f"{k}={v}" for k, v in r.precision.items()) or "fp64 baseline" + print( + f"{label:>16}: total {statistics.median(r.total_times):8.3f} ms " + f"kernel {statistics.median(r.kernel_times):8.3f} ms" + ) + + +if __name__ == "__main__": + main() diff --git a/examples/experiment_perturbation.py b/examples/experiment_perturbation.py new file mode 100644 index 0000000..968001c --- /dev/null +++ b/examples/experiment_perturbation.py @@ -0,0 +1,50 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +"""Measure the heat3d stencil's sensitivity to input noise with the experiment framework.""" + +from scipy import stats + +from corpus.heat3d import GRID_N, TSTEPS, heat3d_kernel +from fp_arena.experiment import ( + ExperimentConfig, + Noise, + PerturbationAnalysisConfig, + ResultStore, + run_perturbation, +) + + +def main(): + experiment = ExperimentConfig( + name="heat3d", + program=heat3d_kernel.to_sdfg(simplify=True), + symbols={"N": GRID_N}, + scalar_args={"TSTEPS": TSTEPS}, + target="cpu", + inputs={ + "A": stats.uniform(0, 100), + "B": stats.uniform(0, 100), + }, + ) + + noise = Noise(relative=1e-3, relative_dist=stats.uniform(-1.0, 2.0)) + store = ResultStore("heat3d_results.db") + results = run_perturbation( + PerturbationAnalysisConfig( + experiment, + noise={"A": noise, "B": noise}, + precisions=[{}], # the unmodified fp64 program + ), + store=store, + ) + + for r in results: + for arr in sorted(r.errors): + e = r.errors[arr] + print( + f"perturbed {r.perturbed} -> {arr}: rel_mean {e.rel_mean:.3e} " + f"rel_max {e.rel_max:.3e}" + ) + + +if __name__ == "__main__": + main() diff --git a/examples/reduced_precision.py b/examples/reduced_precision.py new file mode 100644 index 0000000..eb09b5d --- /dev/null +++ b/examples/reduced_precision.py @@ -0,0 +1,29 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +"""Run the heat3d stencil in reduced precision.""" + +import dace as dc +import numpy as np + +from corpus.heat3d import GRID_N, TSTEPS, heat3d_kernel +from fp_arena.transformations.change_and_propagate_fp_types import ( + change_and_propagate_fp_types, +) + + +def main(): + sdfg = heat3d_kernel.to_sdfg(simplify=True) + + change_and_propagate_fp_types( + sdfg, + {"A": dc.float16, "B": dc.float16}, + constant_type=dc.float16, + ) + + rng = np.random.default_rng(0) + A = rng.uniform(0, 100, (GRID_N,) * 3) + B = rng.uniform(0, 100, (GRID_N,) * 3) + sdfg(TSTEPS=TSTEPS, A=A, B=B, N=GRID_N) + + +if __name__ == "__main__": + main() diff --git a/fp_arena/experiment/__init__.py b/fp_arena/experiment/__init__.py new file mode 100644 index 0000000..244c0fb --- /dev/null +++ b/fp_arena/experiment/__init__.py @@ -0,0 +1,41 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +""" +FP-Arena experiment framework: precision performance and error analysis. +""" + +from fp_arena.experiment.config import ( + ErrorAnalysisConfig, + ExperimentConfig, + PerformanceAnalysisConfig, + PerturbationAnalysisConfig, + PrecisionMap, +) +from fp_arena.experiment.results import ( + ErrorResult, + ErrorStats, + PerfResult, + PerturbationResult, +) +from fp_arena.experiment.inputs import Noise +from fp_arena.experiment import registry +from fp_arena.experiment.runner import run_error, run_performance, run_perturbation +from fp_arena.experiment.store import ResultStore, StoredResult + +__all__ = [ + "ExperimentConfig", + "PerformanceAnalysisConfig", + "ErrorAnalysisConfig", + "PerturbationAnalysisConfig", + "PrecisionMap", + "PerfResult", + "ErrorResult", + "ErrorStats", + "PerturbationResult", + "Noise", + "run_performance", + "run_error", + "run_perturbation", + "registry", + "ResultStore", + "StoredResult", +] diff --git a/fp_arena/experiment/config.py b/fp_arena/experiment/config.py new file mode 100644 index 0000000..6cde3e9 --- /dev/null +++ b/fp_arena/experiment/config.py @@ -0,0 +1,87 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +""" +Experiment config. +""" + +from dataclasses import dataclass, field +from typing import Any, Dict, FrozenSet, List, Optional, Union + +import dace + +from fp_arena.experiment.inputs import DistributionLike, Noise + +#: ``{array_name: precision_key}`` -- pins a subset of arrays to a target format. +PrecisionMap = Dict[str, str] + + +@dataclass +class ExperimentConfig: + """ + The program under test plus everything needed to run it -- shared by all analyses. + + :param name: identifier used to group results in the database. + :param program: the ``dace.SDFG`` under test + :param inputs: per-array input distributions for the *read* arrays. + :param promotion_rules: promotion rules for resolving precision conflicts (e.g. ``{frozenset({fp16, fp32}): fp32}``); ``None`` uses ``DEFAULT_PROMOTION_RULES``. + :param symbols: values for the SDFG's free symbols (e.g. ``{"N": 100}``); + :param scalar_args: values for non-array scalar arguments. + :param target: ``"cpu"`` (default) or ``"gpu"``; + :param seed: base RNG seed for inputs and noise, shared across analyses. + :param gpu_block_size: GPU thread-block size ``[x, y, z]`` (x = contiguous dim) set on every GPU_Device map; ``None`` uses DaCe's default. + """ + + name: str + program: dace.SDFG + inputs: Dict[str, DistributionLike] = field(default_factory=dict) + promotion_rules: Optional[ + Dict[FrozenSet[dace.dtypes.typeclass], dace.dtypes.typeclass] + ] = None + symbols: Dict[str, int] = field(default_factory=dict) + scalar_args: Dict[str, Any] = field(default_factory=dict) + target: str = "cpu" + seed: int = 0 + gpu_block_size: Optional[List[int]] = None + + +@dataclass +class PerformanceAnalysisConfig: + """ + Measure wall-clock runtime across precision points. + """ + + experiment: ExperimentConfig + precisions: List[PrecisionMap] + noise: Dict[str, Noise] = field(default_factory=dict) + n_warmup: int = 1 + n_reps: int = 10 + + +@dataclass +class ErrorAnalysisConfig: + """ + Measure per-array error of each precision point against a high-precision + reference, aggregated over ``n_samples`` input realisations. + ``reference`` is a single key (e.g. ``"mpfr128"``, ``"fp64"``) or a per-array ``{name: key}`` map. + """ + + experiment: ExperimentConfig + precisions: List[PrecisionMap] + noise: Dict[str, Noise] = field(default_factory=dict) + reference: Union[str, Dict[str, str]] = "fp64" + n_samples: int = 1 + + +@dataclass +class PerturbationAnalysisConfig: + """ + Measure input sensitivity: perturb one input array at a time and compare + each written array against the clean run at the same precision point, + aggregated over ``n_samples`` input realisations. + ``noise`` names the inputs to perturb (each analysed separately); + ``precisions`` lists the points to analyse at (default: the unmodified program). + """ + + experiment: ExperimentConfig + noise: Dict[str, Noise] + precisions: List[PrecisionMap] = field(default_factory=lambda: [{}]) + n_samples: int = 1 diff --git a/fp_arena/experiment/inputs.py b/fp_arena/experiment/inputs.py new file mode 100644 index 0000000..9be5dd2 --- /dev/null +++ b/fp_arena/experiment/inputs.py @@ -0,0 +1,191 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +""" +What values go into a run: per-array input generators, optional additive noise, +and the materialisation of a run's call arguments. + +An input is specified in one of two ways: + +* **any SciPy distribution** (anything exposing ``rvs(size=, random_state=)``, + which covers the whole ``scipy.stats`` catalogue.) +* **an initialisation function** ``(shape, rng) -> ndarray``. + +As a default, any read array is sampled from a uniform distribution on [0, 1] without noise. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import ( + AbstractSet, + Any, + Callable, + Dict, + Optional, + Protocol, + Sequence, + TYPE_CHECKING, + Union, + runtime_checkable, +) + +import numpy as np +from scipy import stats as _stats + +import dace + +if TYPE_CHECKING: + from fp_arena.experiment.config import ExperimentConfig + +#: Default generator for unspecified read arrays (uniform on [0, 1]). +_DEFAULT_INPUT = _stats.uniform(0.0, 1.0) + + +@runtime_checkable +class SupportsRVS(Protocol): + """Structural type for SciPy distributions: a ``rvs`` sampling method.""" + + def rvs(self, *args, **kwargs) -> Any: ... + + +class InputDistribution: + """Internal base: produce one array realisation of a given shape/dtype.""" + + def sample( + self, shape: Sequence[int], dtype: np.dtype, rng: np.random.Generator + ) -> np.ndarray: + """Return an array of ``shape``/``dtype`` drawn from ``rng``.""" + raise NotImplementedError + + +@dataclass +class _ScipyDistribution(InputDistribution): + """Adapter around a SciPy distribution (frozen or unfrozen).""" + + dist: SupportsRVS + + def sample(self, shape, dtype, rng): + return np.asarray(self.dist.rvs(size=tuple(shape), random_state=rng)).astype( + dtype + ) + + +@dataclass +class _CallableDistribution(InputDistribution): + """Adapter around a user ``(shape, rng) -> ndarray`` initialisation function.""" + + fn: Callable[[Sequence[int], np.random.Generator], np.ndarray] + + def sample(self, shape, dtype, rng): + return np.asarray(self.fn(tuple(shape), rng)).astype(dtype) + + +DistributionLike = Union[ + SupportsRVS, Callable[[Sequence[int], np.random.Generator], np.ndarray] +] + + +def as_distribution(spec: DistributionLike) -> InputDistribution: + """Coerce a SciPy distribution or ``(shape, rng) -> ndarray`` fn into one.""" + if isinstance(spec, InputDistribution): + return spec + if isinstance(spec, SupportsRVS): + return _ScipyDistribution(spec) + if callable(spec): + return _CallableDistribution(spec) + raise TypeError( + f"Cannot interpret {spec!r} as an input generator; pass a scipy.stats " + f"distribution or a (shape, rng)->ndarray initialisation function" + ) + + +@dataclass +class Noise: + """ + Additive perturbation on an input array: ``x + relative*|x|*R + absolute*A``, + where ``R``/``A`` are draws from ``relative_dist``/``absolute_dist``. Attach + per array via an analysis config's ``noise`` field. + + Each term needs both its magnitude and its distribution; a half-specified + term (or a ``Noise`` that perturbs nothing) raises at construction. + """ + + relative: float = 0.0 + absolute: float = 0.0 + relative_dist: Optional[DistributionLike] = None + absolute_dist: Optional[DistributionLike] = None + + def __post_init__(self) -> None: + for term, mag, dist in ( + ("relative", self.relative, self.relative_dist), + ("absolute", self.absolute, self.absolute_dist), + ): + if (mag != 0.0) != (dist is not None): + raise ValueError( + f"Noise.{term} is half-specified: set both a nonzero " + f"{term} and {term}_dist, or neither" + ) + if self.relative == 0.0 and self.absolute == 0.0: + raise ValueError( + "Noise perturbs nothing: set a relative and/or absolute term" + ) + + def apply(self, arr: np.ndarray, rng: np.random.Generator) -> np.ndarray: + out = np.array(arr, copy=True) + if self.relative: + draw = as_distribution(self.relative_dist).sample(out.shape, out.dtype, rng) + out = out + self.relative * np.abs(out) * draw + if self.absolute: + draw = as_distribution(self.absolute_dist).sample(out.shape, out.dtype, rng) + out = out + self.absolute * draw + return out + + +def materialize_shape(shape, symbols: Dict[str, int]) -> tuple: + """Resolve a (possibly symbolic) shape to a tuple of ints using ``symbols``.""" + out = [] + for dim in shape: + if isinstance(dim, (int, np.integer)): + out.append(int(dim)) + continue + try: + out.append(int(dace.symbolic.evaluate(dim, symbols))) + except Exception as exc: + raise ValueError( + f"Cannot resolve dimension {dim!r} from symbols {symbols}; " + f"provide it in ExperimentConfig.symbols" + ) from exc + return tuple(out) + + +def make_call_args( + sdfg: dace.SDFG, + experiment: "ExperimentConfig", + rng: np.random.Generator, + noise: Optional[Dict[str, Noise]] = None, + reads: Optional[AbstractSet[str]] = None, +) -> Dict[str, Any]: + """ + Build the kwargs for one run: every non-transient array (read arrays sampled + from their distribution then optionally perturbed by ``noise``) plus the symbol and scalar values. + """ + noise = noise or {} + if reads is None: + reads, _ = sdfg.read_and_write_sets() + args: Dict[str, Any] = {} + for name, desc in sdfg.arrays.items(): + if desc.transient or not isinstance(desc, dace.data.Array): + continue + shape = materialize_shape(desc.shape, experiment.symbols) + np_dtype = desc.dtype.as_numpy_dtype() + if name in reads: + dist = as_distribution(experiment.inputs.get(name, _DEFAULT_INPUT)) + arr = dist.sample(shape, np_dtype, rng) + n = noise.get(name) + if n is not None: + arr = n.apply(arr, rng) + args[name] = arr + else: + args[name] = np.zeros(shape, np_dtype) + args.update(experiment.symbols) + args.update(experiment.scalar_args) + return args diff --git a/fp_arena/experiment/registry.py b/fp_arena/experiment/registry.py new file mode 100644 index 0000000..92143f0 --- /dev/null +++ b/fp_arena/experiment/registry.py @@ -0,0 +1,48 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +""" +Mapping between precision keys and DaCe typeclasses, and the promotion rules for them. +""" + +import re +from typing import Dict + +import dace + +import fp_arena # noqa: F401 + +_FIXED: Dict[str, dace.dtypes.typeclass] = { + "fp16": dace.float16, + "fp32": dace.float32, + "fp64": dace.float64, + "fp32sr": dace.float32sr, + "fp64sr": dace.float64sr, +} + +_MPFR_KEY = re.compile(r"^mpfr(\d+)$") + + +def to_typeclass(key: str) -> dace.dtypes.typeclass: + """Resolve a precision key (e.g. ``"fp32"``, ``"mpfr128"``) to its typeclass.""" + if key in _FIXED: + return _FIXED[key] + m = _MPFR_KEY.match(key) + if m: + return dace.mpfr(int(m.group(1))) + raise ValueError( + f"Unknown precision key {key!r}; known: {sorted(_FIXED)} or 'mpfr'" + ) + + +def is_mpfr(key: str) -> bool: + """Whether ``key`` names an MPFR precision (which requires linking libmpfr).""" + return _MPFR_KEY.match(key) is not None + + +def key_of(tc: dace.dtypes.typeclass) -> str: + """Inverse of :func:`to_typeclass`: the stable key for a typeclass.""" + if isinstance(tc, dace.mpfr): + return f"mpfr{tc.precision}" + for key, fixed in _FIXED.items(): + if tc == fixed: + return key + raise ValueError(f"No precision key registered for typeclass {tc!r}") diff --git a/fp_arena/experiment/results.py b/fp_arena/experiment/results.py new file mode 100644 index 0000000..38ef92b --- /dev/null +++ b/fp_arena/experiment/results.py @@ -0,0 +1,101 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +""" +Result records returned by the analysis drivers, one per precision point. Each +serialises to a plain dict (``to_dict``) for the results database. +""" + +from dataclasses import asdict, dataclass, field +from typing import Any, Dict, List + +from fp_arena.experiment.config import PrecisionMap + + +@dataclass +class ErrorStats: + """ + Error metrics of an array versus a reference, reduced over elements and samples. + + Calculations use the error vector `e` (array - reference) and the reference vector `r`. + + Metrics + ------- + abs_mean : Mean(|e|) + + rel_mean : Mean(|e| / |r|) + rel_max : Max(|e| / |r|) + + l1 : sum(|e|) + l2 : sqrt(sum(e**2)) + linf : max(|e|) + + l1_norm : sum(|e|) / sum(|r|) + l2_norm : sqrt(sum(e**2)) / sqrt(sum(r**2)) + linf_norm : max(|e|) / max(|r|) + snr : 10 * log10(sum(r**2) / sum(e**2)) + """ + + abs_mean: float + rel_mean: float + rel_max: float + l1: float + l2: float + linf: float + l1_norm: float + l2_norm: float + linf_norm: float + snr: float + + +@dataclass +class PerfResult: + """ + Per-repetition phase timings (milliseconds) for one precision point, with + ``total = h2d + cast_in + kernel + cast_out + d2h``. Only raw lists are kept; + derive medians/etc. from them. + """ + + precision: PrecisionMap + total_times: List[float] = field(default_factory=list) + h2d_times: List[float] = field( + default_factory=list + ) # host->device transfer (GPU only) + d2h_times: List[float] = field( + default_factory=list + ) # device->host transfer (GPU only) + cast_in_times: List[float] = field(default_factory=list) # input precision cast + cast_out_times: List[float] = field(default_factory=list) # output precision cast + kernel_times: List[float] = field(default_factory=list) # compute + seed: int = 0 + + def to_dict(self) -> Dict[str, Any]: + return asdict(self) + + +@dataclass +class ErrorResult: + """Per-array error for one precision point.""" + + precision: PrecisionMap + errors: Dict[str, ErrorStats] + n_samples: int + seed: int = 0 + + def to_dict(self) -> Dict[str, Any]: + return asdict(self) + + +@dataclass +class PerturbationResult: + """ + Per-array output deviation induced by perturbing the single input + ``perturbed``, versus the clean run at the same precision point. + """ + + precision: PrecisionMap + perturbed: str + errors: Dict[str, ErrorStats] + n_samples: int + seed: int = 0 + + def to_dict(self) -> Dict[str, Any]: + return asdict(self) diff --git a/fp_arena/experiment/retarget.py b/fp_arena/experiment/retarget.py new file mode 100644 index 0000000..a0cdd49 --- /dev/null +++ b/fp_arena/experiment/retarget.py @@ -0,0 +1,114 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +""" +The bridge between an :class:`ExperimentConfig` and DaCe: build a fresh SDFG and retarget its precision using :func:`change_and_propagate_fp_types`. +""" + +import copy +from typing import List, Optional + +import dace + +from fp_arena.transformations.change_and_propagate_fp_types import ( + change_and_propagate_fp_types, +) +from fp_arena.experiment import registry +from fp_arena.experiment.config import ExperimentConfig, PrecisionMap + + +def _is_fp(tc: dace.dtypes.typeclass) -> bool: + """:returns: whether ``tc`` is a registry float (native, SR, or MPFR).""" + try: + registry.key_of(tc) + return True + except ValueError: + return False + + +def fresh_sdfg(experiment: ExperimentConfig) -> dace.SDFG: + """ + A deep copy of the experiment's SDFG -- lowering mutates in place, so every precision point needs its own. + """ + if not isinstance(experiment.program, dace.SDFG): + raise TypeError( + f"program must be a dace.SDFG, got {type(experiment.program).__name__}" + ) + return copy.deepcopy(experiment.program) + + +def candidate_fp_arrays(sdfg: dace.SDFG) -> List[str]: + """:returns: non-transient floating-point array names (the boundary fp arrays).""" + return sorted( + name + for name, desc in sdfg.arrays.items() + if not desc.transient + and _is_fp(desc.dtype) + and isinstance(desc, dace.data.Array) + ) + + +def _validate_pins(sdfg: dace.SDFG, pin_map: PrecisionMap) -> None: + for name in pin_map: + if name not in sdfg.arrays: + raise ValueError(f"Pinned array {name!r} not found in SDFG {sdfg.name!r}") + if not _is_fp(sdfg.arrays[name].dtype): + raise ValueError(f"Pinned array {name!r} is not a floating-point array") + + +def _ensure_mpfr_linked() -> None: + """Add the MPFR library to DaCe's CPU link line.""" + libs = dace.Config.get("compiler", "cpu", "libs") or "" + if "mpfr" not in libs.split(): + dace.Config.append("compiler", "cpu", "libs", value=" mpfr") + + +def apply_precision( + sdfg: dace.SDFG, + pin_map: PrecisionMap, + promotion_rules, +) -> None: + """ + Retarget ``sdfg`` in place to the precisions in ``pin_map`` via :func:`change_and_propagate_fp_types` + """ + if not pin_map: + return + _validate_pins(sdfg, pin_map) + if any(registry.is_mpfr(key) for key in pin_map.values()): + _ensure_mpfr_linked() + typed = {name: registry.to_typeclass(key) for name, key in pin_map.items()} + change_and_propagate_fp_types(sdfg, typed, promotion_rules) + + +def apply_reference(sdfg: dace.SDFG, reference, promotion_rules) -> None: + """ + Lower ``sdfg`` in place to the reference precision (a single key for every fp array, or a per-array ``{name: key}`` map). + """ + if isinstance(reference, str): + ref_map = {name: reference for name in candidate_fp_arrays(sdfg)} + else: + ref_map = dict(reference) + apply_precision(sdfg, ref_map, promotion_rules) + + +def apply_target( + sdfg: dace.SDFG, + target: str, + gpu_block_size: Optional[List[int]] = None, +) -> None: + """ + Retarget ``sdfg`` in place for the execution target. + """ + if target == "cpu": + return + if target == "gpu": + # simplify=False keeps host<->device copies and the kernel in separate states so each timing phase is attributable + sdfg.apply_gpu_transformations(simplify=False) + if gpu_block_size is not None: + for state in sdfg.all_states(): + for node in state.nodes(): + if ( + isinstance(node, dace.nodes.MapEntry) + and node.map.schedule == dace.dtypes.ScheduleType.GPU_Device + ): + node.map.gpu_block_size = list(gpu_block_size) + return + raise ValueError(f"Unknown target {target!r}; expected 'cpu' or 'gpu'") diff --git a/fp_arena/experiment/runner.py b/fp_arena/experiment/runner.py new file mode 100644 index 0000000..85f4d1b --- /dev/null +++ b/fp_arena/experiment/runner.py @@ -0,0 +1,492 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +""" +Drivers that execute the experiment kinds. + +* :func:`run_performance` -- compile each precision point and time it. +* :func:`run_error` -- dual-execute each point against a high-precision reference (lockstep on identical inputs). +* :func:`run_perturbation` -- perturb one input at a time and compare against the clean run at the same precision point. +""" + +import math +from typing import Any, Dict, List, Optional, Tuple + +import numpy as np +from tqdm.auto import tqdm + +import dace + +from fp_arena.experiment.config import ( + ErrorAnalysisConfig, + PerformanceAnalysisConfig, + PerturbationAnalysisConfig, + PrecisionMap, +) +from fp_arena.experiment.results import ( + ErrorResult, + ErrorStats, + PerfResult, + PerturbationResult, +) +from fp_arena.experiment.inputs import make_call_args +from fp_arena.experiment.retarget import ( + apply_precision, + apply_reference, + apply_target, + fresh_sdfg, +) +from fp_arena.experiment.store import ResultStore + + +def _copy_args(args: Dict[str, Any]) -> Dict[str, Any]: + """Independent copy of the array arguments; scalars/symbols pass through.""" + return { + k: (np.array(v, copy=True) if isinstance(v, np.ndarray) else v) + for k, v in args.items() + } + + +def _reset_arrays(working: Dict[str, Any], source: Dict[str, Any]) -> None: + """Refill ``working``'s array buffers in place from ``source``.""" + for k, v in source.items(): + if isinstance(v, np.ndarray): + working[k][...] = v + + +def _sample_rngs(seed: int, n: int) -> List[np.random.Generator]: + """ + ``n`` independent, reproducible per-sample generators from one base seed. + """ + return [ + np.random.default_rng(child) for child in np.random.SeedSequence(seed).spawn(n) + ] + + +def _new_acc() -> Dict[str, float]: + return { + "abs_sum": 0.0, + "abs_cnt": 0, + "abs_max": 0.0, + "rel_sum": 0.0, + "rel_cnt": 0, + "rel_max": 0.0, + "sq_err_sum": 0.0, + "sq_ref_sum": 0.0, + "ref_abs_sum": 0.0, + "ref_abs_max": 0.0, + } + + +def _accumulate(acc: Dict[str, float], ref: np.ndarray, cand: np.ndarray) -> None: + """Fold one (reference, candidate) array pair into the running error stats.""" + r = np.asarray(ref, dtype=np.float64).ravel() + c = np.asarray(cand, dtype=np.float64).ravel() + diff = np.abs(c - r) + diff = np.where(np.isnan(diff), np.inf, diff) + denom = np.abs(r) + if diff.size: + acc["abs_sum"] += float(diff.sum()) + acc["abs_cnt"] += int(diff.size) + acc["abs_max"] = max(acc["abs_max"], float(diff.max())) + acc["sq_err_sum"] += float(np.square(diff).sum()) + acc["sq_ref_sum"] += float(np.square(r).sum()) + acc["ref_abs_sum"] += float(denom.sum()) + acc["ref_abs_max"] = max(acc["ref_abs_max"], float(denom.max())) + mask = denom > 0 + with np.errstate(invalid="ignore"): + rel = diff[mask] / denom[mask] + rel = np.where(np.isnan(rel), np.inf, rel) + if rel.size: + acc["rel_sum"] += float(rel.sum()) + acc["rel_cnt"] += int(rel.size) + acc["rel_max"] = max(acc["rel_max"], float(rel.max())) + + +def _ratio(num: float, den: float) -> float: + """``num/den`` with the zero-reference convention: ``0/0 -> 0``, ``x/0 -> inf``.""" + if den > 0.0: + return num / den + return 0.0 if num <= 0.0 else math.inf + + +def _finalize(acc: Dict[str, float]) -> ErrorStats: + abs_mean = acc["abs_sum"] / acc["abs_cnt"] if acc["abs_cnt"] else 0.0 + rel_mean = acc["rel_sum"] / acc["rel_cnt"] if acc["rel_cnt"] else 0.0 + err_power = acc["sq_err_sum"] + ref_power = acc["sq_ref_sum"] + if err_power <= 0.0: + snr = math.inf + elif ( + ref_power <= 0.0 or not math.isfinite(ref_power) or not math.isfinite(err_power) + ): + snr = -math.inf + else: + snr = 10.0 * math.log10(ref_power / err_power) + l1 = acc["abs_sum"] + l2 = math.sqrt(err_power) + linf = acc["abs_max"] + return ErrorStats( + abs_mean=abs_mean, + rel_mean=rel_mean, + rel_max=acc["rel_max"], + l1=l1, + l2=l2, + linf=linf, + l1_norm=_ratio(l1, acc["ref_abs_sum"]), + l2_norm=_ratio(l2, math.sqrt(ref_power)), + linf_norm=_ratio(linf, acc["ref_abs_max"]), + snr=snr, + ) + + +#: Device storage types; an edge crossing this boundary is a host<->device copy. +_GPU_STORAGE = ( + dace.dtypes.StorageType.GPU_Global, + dace.dtypes.StorageType.GPU_Shared, +) + + +def _transfer_direction(sdfg: dace.SDFG, state) -> Optional[str]: + """``"h2d"``/``"d2h"`` if ``state`` copies across the host/device boundary, else ``None``.""" + for e in state.edges(): + src, dst = e.src, e.dst + if isinstance(src, dace.nodes.AccessNode) and isinstance( + dst, dace.nodes.AccessNode + ): + src_dev = sdfg.arrays[src.data].storage in _GPU_STORAGE + dst_dev = sdfg.arrays[dst.data].storage in _GPU_STORAGE + if src_dev != dst_dev: + return "h2d" if dst_dev else "d2h" + return None + + +def _has_cast_map(state) -> bool: + """Whether ``state`` contains a precision-cast map.""" + return any( + isinstance(n, dace.nodes.MapEntry) and n.map.label.startswith("cast_map_") + for n in state.nodes() + ) + + +def _classify_state(sdfg: dace.SDFG, state) -> str: + """Assign ``state`` to a timing phase: h2d/d2h, cast_in/cast_out, or kernel.""" + direction = _transfer_direction(sdfg, state) + if direction is not None: + return direction + if _has_cast_map(state): + return "cast_out" if state.label.startswith("copy_out") else "cast_in" + return "kernel" + + +def _group_per_invocation( + samples: List[float], n_invocations: int, name: str +) -> List[float]: + """Sum a timer's per-execution samples into one value per invocation (loop bodies fire repeatedly).""" + if not samples: + return [] + if len(samples) % n_invocations != 0: + raise ValueError( + f"Timer {name!r} fired {len(samples)} times over {n_invocations} " + f"invocations; per-rep phase attribution requires a static " + f"per-invocation execution count" + ) + k = len(samples) // n_invocations + if k == 1: + return list(samples) + return [math.fsum(samples[i * k : (i + 1) * k]) for i in range(n_invocations)] + + +def _phase_series(report, name_pred, n_reps: int, n_invocations: int) -> List[float]: + """Per-rep milliseconds summed over matching timers, warmup invocations dropped.""" + out = [0.0] * n_reps + matched = False + if report is not None: + for names in report.durations.values(): + for name, tid_map in names.items(): + if not name_pred(name): + continue + for times_ms in tid_map.values(): + per_inv = _group_per_invocation(times_ms, n_invocations, name) + if not per_inv: + continue + for i, ms in enumerate(per_inv[-n_reps:]): + out[i] += ms + matched = True + return out if matched else [] + + +def _phase_breakdown( + sdfg: dace.SDFG, n_reps: int, n_invocations: int +) -> Dict[str, Any]: + """Reduce the latest report into per-rep phase series; ``total`` is their sum.""" + report = sdfg.get_latest_report() + # Map each state's report name -> phase category. + cat_of: Dict[str, str] = {} + for state in sdfg.all_states(): + cat_of.setdefault(f"State {state.label}", _classify_state(sdfg, state)) + + def series(category: str) -> List[float]: + return _phase_series( + report, lambda nm: cat_of.get(nm) == category, n_reps, n_invocations + ) + + phases = { + "h2d_times": series("h2d"), + "cast_in_times": series("cast_in"), + "kernel_times": series("kernel"), + "cast_out_times": series("cast_out"), + "d2h_times": series("d2h"), + } + nonempty = [p for p in phases.values() if p] + total = [math.fsum(col) for col in zip(*nonempty)] if nonempty else [] + return {"total_times": total, **phases} + + +def _fmt_pin(pin_map: PrecisionMap) -> str: + """Compact one-line summary of a precision point for progress display.""" + return " ".join(f"{k}={v}" for k, v in pin_map.items()) + + +def _output_arrays(sdfg: dace.SDFG) -> List[str]: + """Non-transient written array names -- the arrays error metrics are reduced over.""" + _, writes = sdfg.read_and_write_sets() + return sorted( + name + for name in writes + if name in sdfg.arrays + and isinstance(sdfg.arrays[name], dace.data.Array) + and not sdfg.arrays[name].transient + ) + + +#: Per input sample: (pristine call args, reference outputs of the written arrays). +ReferenceSamples = List[Tuple[Dict[str, Any], Dict[str, Any]]] + + +def _pin_tag(pin_map: PrecisionMap) -> str: + """Identifier-safe tag naming a precision point's build folder.""" + return "_".join(f"{k}_{v}" for k, v in sorted(pin_map.items())) or "baseline" + + +def _distinguish(sdfg: dace.SDFG, tag: str) -> None: + """Unique build folder per point.""" + safe = "".join(ch if ch.isalnum() else "_" for ch in tag) + sdfg.name = f"{sdfg.name}_{safe}" + + +def compile_reference(experiment, reference): + """Build and compile the high-precision reference SDFG once.""" + sdfg = fresh_sdfg(experiment) + apply_reference(sdfg, reference, experiment.promotion_rules) + apply_target(sdfg, experiment.target, gpu_block_size=experiment.gpu_block_size) + _distinguish(sdfg, "reference") + return sdfg.compile() + + +def run_reference( + experiment, reference, n_samples: int, seed: int, noise=None +) -> ReferenceSamples: + """ + Execute the reference once per input sample. + """ + ref_csdfg = compile_reference(experiment, reference) + program = experiment.program + reads, _ = program.read_and_write_sets() + outputs = _output_arrays(program) + samples: ReferenceSamples = [] + for rng in tqdm( + _sample_rngs(seed, n_samples), desc="reference", unit="smp", leave=False + ): + args = make_call_args(program, experiment, rng, noise, reads=reads) + ref_args = _copy_args(args) + ref_csdfg(**ref_args) + samples.append((args, {name: ref_args[name] for name in outputs})) + return samples + + +def measure_error( + experiment, + pin_map: PrecisionMap, + ref_samples: ReferenceSamples, + seed: int, +) -> ErrorResult: + """ + Execute one precision point on the reference's input samples and reduce the per-array error against the cached reference outputs. + """ + cand_sdfg = fresh_sdfg(experiment) + apply_precision(cand_sdfg, pin_map, experiment.promotion_rules) + apply_target(cand_sdfg, experiment.target, gpu_block_size=experiment.gpu_block_size) + _distinguish(cand_sdfg, _pin_tag(pin_map)) + cand_csdfg = cand_sdfg.compile() + + acc = {name: _new_acc() for name in _output_arrays(cand_sdfg)} + for args, ref_out in tqdm(ref_samples, desc="samples", unit="smp", leave=False): + cand_args = _copy_args(args) + cand_csdfg(**cand_args) + for name in acc: + _accumulate(acc[name], ref_out[name], cand_args[name]) + + errors = {name: _finalize(a) for name, a in acc.items()} + return ErrorResult( + precision=dict(pin_map), errors=errors, n_samples=len(ref_samples), seed=seed + ) + + +def run_performance( + cfg: PerformanceAnalysisConfig, store: Optional[ResultStore] = None +) -> List[PerfResult]: + """Time each precision point, optionally appending to ``store``.""" + results: List[PerfResult] = [] + target = cfg.experiment.target + # CPU: host std::chrono; GPU: CUDA events (on-device, not async-launch, time). + provider = ( + dace.InstrumentationType.GPU_Events + if target == "gpu" + else dace.InstrumentationType.Timer + ) + n_invocations = cfg.n_warmup + cfg.n_reps + + prev_each = dace.Config.get("instrumentation", "report_each_invocation") + dace.Config.set("instrumentation", "report_each_invocation", value=False) + prev_streams = dace.Config.get("compiler", "cuda", "max_concurrent_streams") + if target == "gpu": + dace.Config.set("compiler", "cuda", "max_concurrent_streams", value=-1) + + points = tqdm(cfg.precisions, desc="performance", unit="pt") + try: + for pin_map in points: + points.set_postfix_str(_fmt_pin(pin_map)) + sdfg = fresh_sdfg(cfg.experiment) + apply_precision(sdfg, pin_map, cfg.experiment.promotion_rules) + apply_target(sdfg, target, gpu_block_size=cfg.experiment.gpu_block_size) + _distinguish(sdfg, _pin_tag(pin_map)) + # Time every state; classified into a phase at readout. + for state in sdfg.all_states(): + state.instrument = provider + csdfg = sdfg.compile() + sdfg.clear_instrumentation_reports() + + rng = _sample_rngs(cfg.experiment.seed, 1)[0] + initial_args = make_call_args(sdfg, cfg.experiment, rng, cfg.noise) + args = _copy_args(initial_args) + runs = tqdm( + total=n_invocations, + desc="warmup", + unit="run", + leave=False, + ) + for _ in range(cfg.n_warmup): + _reset_arrays(args, initial_args) + csdfg(**args) + runs.update(1) + runs.set_description("reps") + for _ in range(cfg.n_reps): + _reset_arrays(args, initial_args) + csdfg(**args) + runs.update(1) + runs.close() + + csdfg.finalize() + + result = PerfResult( + precision=dict(pin_map), + seed=cfg.experiment.seed, + **_phase_breakdown(sdfg, cfg.n_reps, n_invocations), + ) + results.append(result) + if store is not None: + store.add( + cfg.experiment.name, + result, + symbols=cfg.experiment.symbols, + scalars=cfg.experiment.scalar_args, + ) + finally: + dace.Config.set("instrumentation", "report_each_invocation", value=prev_each) + dace.Config.set( + "compiler", "cuda", "max_concurrent_streams", value=prev_streams + ) + return results + + +def run_perturbation( + cfg: PerturbationAnalysisConfig, store: Optional[ResultStore] = None +) -> List[PerturbationResult]: + """ + Measure per-array output sensitivity to input noise: for each precision + point, execute clean and perturbed inputs (one noisy array at a time) on the + same compiled SDFG and reduce the output deviation over ``cfg.n_samples`` + input realisations. One result per (precision point, perturbed input). + """ + if not cfg.noise: + raise ValueError( + "PerturbationAnalysisConfig.noise must name at least one input array" + ) + exp = cfg.experiment + results: List[PerturbationResult] = [] + points = tqdm(cfg.precisions, desc="perturbation", unit="pt") + for pin_map in points: + points.set_postfix_str(_fmt_pin(pin_map)) + sdfg = fresh_sdfg(exp) + apply_precision(sdfg, pin_map, exp.promotion_rules) + apply_target(sdfg, exp.target, gpu_block_size=exp.gpu_block_size) + _distinguish(sdfg, _pin_tag(pin_map)) + csdfg = sdfg.compile() + + reads, _ = sdfg.read_and_write_sets() + for name in cfg.noise: + if name not in reads or name not in sdfg.arrays: + raise ValueError( + f"Perturbed array {name!r} is not a read input of SDFG {sdfg.name!r}" + ) + outputs = _output_arrays(sdfg) + acc = {pert: {out: _new_acc() for out in outputs} for pert in cfg.noise} + + for rng in tqdm( + _sample_rngs(exp.seed, cfg.n_samples), + desc="samples", + unit="smp", + leave=False, + ): + clean_args = make_call_args(sdfg, exp, rng, reads=reads) + base_args = _copy_args(clean_args) + csdfg(**base_args) + for pert_name, pert_noise in cfg.noise.items(): + pert_args = _copy_args(clean_args) + pert_args[pert_name] = pert_noise.apply(clean_args[pert_name], rng) + csdfg(**pert_args) + for out in outputs: + _accumulate(acc[pert_name][out], base_args[out], pert_args[out]) + + for pert_name in cfg.noise: + result = PerturbationResult( + precision=dict(pin_map), + perturbed=pert_name, + errors={out: _finalize(a) for out, a in acc[pert_name].items()}, + n_samples=cfg.n_samples, + seed=exp.seed, + ) + results.append(result) + if store is not None: + store.add( + exp.name, result, symbols=exp.symbols, scalars=exp.scalar_args + ) + return results + + +def run_error( + cfg: ErrorAnalysisConfig, store: Optional[ResultStore] = None +) -> List[ErrorResult]: + """Measure per-array error of each precision point, optionally appending to ``store`` database.""" + exp = cfg.experiment + ref_samples = run_reference( + exp, cfg.reference, cfg.n_samples, exp.seed, noise=cfg.noise + ) + results: List[ErrorResult] = [] + points = tqdm(cfg.precisions, desc="error", unit="pt") + for pin_map in points: + points.set_postfix_str(_fmt_pin(pin_map)) + result = measure_error(exp, pin_map, ref_samples, exp.seed) + results.append(result) + if store is not None: + store.add(exp.name, result, symbols=exp.symbols, scalars=exp.scalar_args) + return results diff --git a/fp_arena/experiment/store.py b/fp_arena/experiment/store.py new file mode 100644 index 0000000..a287c06 --- /dev/null +++ b/fp_arena/experiment/store.py @@ -0,0 +1,119 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. +""" +Append-only results database. +""" + +import json +import sqlite3 +from dataclasses import dataclass +from datetime import datetime, timezone +from typing import Any, Dict, List, Optional, Union + +from fp_arena.experiment.results import ErrorResult, PerfResult, PerturbationResult + +DEFAULT_DB_PATH = ".fp_arena_results.db" + +#: The ``kind`` column value for each result type. +_KIND_OF = { + PerfResult: "performance", + ErrorResult: "error", + PerturbationResult: "perturbation", +} + +_SCHEMA = """ +CREATE TABLE IF NOT EXISTS results ( + id INTEGER PRIMARY KEY AUTOINCREMENT, + kind TEXT NOT NULL, -- 'performance' | 'error' | 'perturbation' + experiment TEXT NOT NULL, -- ExperimentConfig.name + created_at TEXT NOT NULL, -- ISO-8601 UTC + precision TEXT NOT NULL, -- JSON: the precision map + symbols TEXT NOT NULL, -- JSON: the SDFG free-symbol values + scalars TEXT NOT NULL, -- JSON: the non-array scalar arguments + payload TEXT NOT NULL -- JSON: the full result record +); +CREATE INDEX IF NOT EXISTS idx_results_lookup ON results (experiment, kind); +""" + + +@dataclass +class StoredResult: + """One row from the database.""" + + id: int + kind: str + experiment: str + created_at: str + precision: Dict[str, Any] + symbols: Dict[str, Any] + scalars: Dict[str, Any] + payload: Dict[str, Any] + + +class ResultStore: + """A thin SQLite wrapper for appending and querying experiment results.""" + + def __init__(self, path: str = DEFAULT_DB_PATH): + self.path = path + self._con = sqlite3.connect(path, check_same_thread=False) + self._con.row_factory = sqlite3.Row + self._con.executescript(_SCHEMA) + self._con.commit() + + def close(self) -> None: + self._con.close() + + def add( + self, + experiment: str, + result: Union[PerfResult, ErrorResult, PerturbationResult], + symbols: Optional[Dict[str, Any]] = None, + scalars: Optional[Dict[str, Any]] = None, + ) -> int: + """Append one result record. :returns: the new row id.""" + kind = _KIND_OF.get(type(result)) + if kind is None: + raise TypeError(f"Cannot store a {type(result).__name__}") + with self._con: + cur = self._con.execute( + "INSERT INTO results (kind, experiment, created_at, precision, symbols, scalars, payload) " + "VALUES (?, ?, ?, ?, ?, ?, ?)", + ( + kind, + experiment, + datetime.now(timezone.utc).isoformat(), + json.dumps(result.precision), + json.dumps(symbols or {}, default=str), + json.dumps(scalars or {}, default=str), + json.dumps(result.to_dict()), + ), + ) + return int(cur.lastrowid) + + def query( + self, experiment: Optional[str] = None, kind: Optional[str] = None + ) -> List[StoredResult]: + """Fetch stored results, newest first, optionally filtered by experiment/kind.""" + clauses, params = [], [] + if experiment is not None: + clauses.append("experiment = ?") + params.append(experiment) + if kind is not None: + clauses.append("kind = ?") + params.append(kind) + where = (" WHERE " + " AND ".join(clauses)) if clauses else "" + rows = self._con.execute( + f"SELECT * FROM results{where} ORDER BY id DESC", params + ).fetchall() + return [ + StoredResult( + id=r["id"], + kind=r["kind"], + experiment=r["experiment"], + created_at=r["created_at"], + precision=json.loads(r["precision"]), + symbols=json.loads(r["symbols"]), + scalars=json.loads(r["scalars"]), + payload=json.loads(r["payload"]), + ) + for r in rows + ] diff --git a/fp_arena/transformations/change_and_propagate_fp_types.py b/fp_arena/transformations/change_and_propagate_fp_types.py index 7357e5a..0c20fe8 100644 --- a/fp_arena/transformations/change_and_propagate_fp_types.py +++ b/fp_arena/transformations/change_and_propagate_fp_types.py @@ -6,6 +6,7 @@ from dace.properties import CodeBlock from dace.sdfg import nodes, utils as sdfg_utils from dace.sdfg.state import AbstractControlFlowRegion, SDFGState +from tqdm.auto import tqdm from fp_arena.dtypes import float32sr, float64sr @@ -213,7 +214,7 @@ def _print_type_report( lines.append( f" {name:<{name_w}} {orig.to_string():<9} {final.to_string():<9}{marker}" ) - print("\n".join(lines)) + tqdm.write("\n".join(lines)) # Writes inferred types onto tasklet/map/library-node connectors. diff --git a/pyproject.toml b/pyproject.toml index 4bb188c..c561e3d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,6 +26,8 @@ keywords = ["dace", "floating-point", "stochastic-rounding", "reduced-precision" # that DaCe is importable; see the README. dependencies = [ "dace @ git+https://github.com/spcl/dace.git@extended", + "scipy", + "tqdm", ] [project.optional-dependencies] @@ -44,7 +46,7 @@ Issues = "https://github.com/spcl/FP-Arena/issues" [tool.setuptools.packages.find] where = ["."] -include = ["fp_arena*"] +include = ["fp_arena*", "corpus*"] [tool.setuptools.package-data] fp_arena = ["runtime/include/**/*.h"] diff --git a/tests/test_change_and_propagate_fp_types.py b/tests/test_change_and_propagate_fp_types.py index c01a4b8..cb3aedb 100644 --- a/tests/test_change_and_propagate_fp_types.py +++ b/tests/test_change_and_propagate_fp_types.py @@ -1,5 +1,6 @@ import dace import pytest +from corpus.heat3d import heat3d_kernel from dace.libraries.standard.nodes.reduce import Reduce from fp_arena.transformations.change_and_propagate_fp_types import ( DEFAULT_PROMOTION_RULES, @@ -638,35 +639,11 @@ def prog(A: dace.float64[n], B: dace.float64[n]): def test_heat3d_no_fp64_in_generated_code(): - """heat3d (examples/heat3d.py) lowered from fp64 to fp16: the generated + """heat3d (corpus/heat3d.py) lowered from fp64 to fp16: the generated C++ contains fp64 only at the preserved A/B interface.""" import re - N = dace.symbol("N", dtype=dace.int64) - - @dace.program - def heat3d(TSTEPS: dace.int64, A: dace.float64[N, N, N], B: dace.float64[N, N, N]): - for t in range(1, TSTEPS): - B[1:-1, 1:-1, 1:-1] = ( - 0.125 - * (A[2:, 1:-1, 1:-1] - 2.0 * A[1:-1, 1:-1, 1:-1] + A[:-2, 1:-1, 1:-1]) - + 0.125 - * (A[1:-1, 2:, 1:-1] - 2.0 * A[1:-1, 1:-1, 1:-1] + A[1:-1, :-2, 1:-1]) - + 0.125 - * (A[1:-1, 1:-1, 2:] - 2.0 * A[1:-1, 1:-1, 1:-1] + A[1:-1, 1:-1, :-2]) - + A[1:-1, 1:-1, 1:-1] - ) - A[1:-1, 1:-1, 1:-1] = ( - 0.125 - * (B[2:, 1:-1, 1:-1] - 2.0 * B[1:-1, 1:-1, 1:-1] + B[:-2, 1:-1, 1:-1]) - + 0.125 - * (B[1:-1, 2:, 1:-1] - 2.0 * B[1:-1, 1:-1, 1:-1] + B[1:-1, :-2, 1:-1]) - + 0.125 - * (B[1:-1, 1:-1, 2:] - 2.0 * B[1:-1, 1:-1, 1:-1] + B[1:-1, 1:-1, :-2]) - + B[1:-1, 1:-1, 1:-1] - ) - - sdfg = heat3d.to_sdfg(simplify=True) + sdfg = heat3d_kernel.to_sdfg(simplify=True) change_and_propagate_fp_types( sdfg, {"A": dace.float16, "B": dace.float16}, @@ -677,9 +654,9 @@ def heat3d(TSTEPS: dace.int64, A: dace.float64[N, N, N], B: dace.float64[N, N, N code_objects = sdfg.generate_code() # The computation is in fp16. - assert any( - "dace::float16" in co.clean_code for co in code_objects - ), "expected dace::float16 in the generated code" + assert any("dace::float16" in co.clean_code for co in code_objects), ( + "expected dace::float16 in the generated code" + ) # fp64 only on interface lines: A/B references and interface cast tasklets. allowed = re.compile(r"\b[AB]\b|\bdouble\s+_out;|=\s*double\(_in\)") @@ -692,6 +669,100 @@ def heat3d(TSTEPS: dace.int64, A: dace.float64[N, N, N], B: dace.float64[N, N, N assert not leaks, "fp64 leaked into the computation:\n" + "\n".join(leaks) +def test_end_to_end_demoted_written_array_runs(): + """A *written* array is demoted below the precision of the f64 computation that produces it. + + Chain A(f64) -> B -> C with B pinned to f32. The producer of B reads A (f64), + so its tasklet computes in f64 but must store into the f32 array B. + """ + import numpy as np + + n = 8 + sdfg = dace.SDFG("e2e_demoted_write") + sdfg.add_array("A", [n], dace.float64, transient=False) + sdfg.add_array("B", [n], dace.float64, transient=True) + sdfg.add_array("C", [n], dace.float64, transient=False) + + s1 = sdfg.add_state("s1") + s2 = sdfg.add_state("s2") + sdfg.add_edge(s1, s2, dace.InterstateEdge()) + + for st, src, dst in [(s1, "A", "B"), (s2, "B", "C")]: + me, mx = st.add_map("m", {"i": f"0:{n}"}) + t = st.add_tasklet( + "t", {"x"}, {"y"}, "y = x * 2.0;", language=dace.Language.CPP + ) + me.add_in_connector(f"IN_{src}") + me.add_out_connector(f"OUT_{src}") + mx.add_in_connector(f"IN_{dst}") + mx.add_out_connector(f"OUT_{dst}") + st.add_edge( + st.add_read(src), None, me, f"IN_{src}", dace.Memlet(f"{src}[0:{n}]") + ) + st.add_edge(me, f"OUT_{src}", t, "x", dace.Memlet(f"{src}[i]")) + st.add_edge(t, "y", mx, f"IN_{dst}", dace.Memlet(f"{dst}[i]")) + st.add_edge( + mx, f"OUT_{dst}", st.add_write(dst), None, dace.Memlet(f"{dst}[0:{n}]") + ) + + # Demote only the WRITTEN intermediate B; A stays f64, so the producer of B + # computes in f64 and must store into an f32 array (the heat3d crash pattern). + change_and_propagate_fp_types(sdfg, {"B": dace.float32}) + sdfg.validate() + + assert sdfg.arrays["B"].dtype == dace.float32, sdfg.arrays["B"].dtype + assert sdfg.arrays["A"].dtype == dace.float64 + + A = np.arange(1, n + 1, dtype=np.float64) + C = np.zeros(n, dtype=np.float64) + sdfg(A=A, C=C) + + # B = (f32)(A*2); C = (f32)(B*2). Small integers are exact in f32, so C == A*4. + b_ref = (A * 2.0).astype(np.float32) + c_ref = (b_ref.astype(np.float64) * 2.0).astype(np.float32).astype(np.float64) + np.testing.assert_allclose(C, c_ref, rtol=1e-6) + np.testing.assert_allclose(C, A * 4.0, rtol=1e-6) + + +def test_boundary_cast_inserted_for_fusion(): + """Regression: a fused map with mixed precision compiles via a map-boundary cast.""" + import numpy as np + from dace.sdfg import nodes + from dace.transformation.dataflow import MapFusion + + M = dace.symbol("M") + + @dace.program + def prog(A: dace.float64[M], B: dace.float64[M]): + B[:] = A * 2.0 + A[:] = B * 3.0 + + # Fuse the two statements through a transient holding B's value. + sdfg = prog.to_sdfg(simplify=True) + sdfg.apply_transformations_repeated(MapFusion) + + # fp16 transient written into fp32 B -> a cast must be inserted. + change_and_propagate_fp_types( + sdfg, {"A": dace.float16, "B": dace.float32}, DEFAULT_PROMOTION_RULES + ) + sdfg.validate() + casts = [ + n.label + for s in sdfg.all_states() + for n in s.nodes() + if isinstance(n, nodes.Tasklet) and "map_fusion_B_to_B" in n.label + ] + assert casts, "expected a cast on the fused transient's write into B" + + csdfg = sdfg.compile() + A = np.full(4, 3.0, dtype=np.float64) + B = np.zeros(4, dtype=np.float64) + csdfg(A=A, B=B, M=4) + # B = A*2 = 6, A = B*3 = 18 (exact in fp16/fp32). + np.testing.assert_allclose(B, 6.0) + np.testing.assert_allclose(A, 18.0) + + if __name__ == "__main__": test_transient_intermediate_propagates() test_all_nontransient_interface_preserved() @@ -715,4 +786,6 @@ def heat3d(TSTEPS: dace.int64, A: dace.float64[N, N, N], B: dace.float64[N, N, N test_constant_type_only_touches_python_tasklets() test_constant_type_end_to_end_float32_precision() test_heat3d_no_fp64_in_generated_code() + test_end_to_end_demoted_written_array_runs() + test_boundary_cast_inserted_for_fusion() print("All tests passed.") diff --git a/tests/test_experiment.py b/tests/test_experiment.py new file mode 100644 index 0000000..4da1b46 --- /dev/null +++ b/tests/test_experiment.py @@ -0,0 +1,503 @@ +# Copyright 2019-2026 ETH Zurich and the FP-Arena authors. All rights reserved. + +import dace +import numpy as np +import pytest +from scipy import stats +import fp_arena # noqa: F401 + +from fp_arena.experiment import ( + ErrorAnalysisConfig, + ExperimentConfig, + Noise, + PerformanceAnalysisConfig, + PerturbationAnalysisConfig, + ResultStore, + registry, + run_error, + run_performance, + run_perturbation, +) +from fp_arena.experiment.inputs import make_call_args +from fp_arena.experiment.retarget import apply_target, fresh_sdfg +from fp_arena.experiment.runner import ( + _accumulate, + _finalize, + _group_per_invocation, + _new_acc, +) + +N = dace.symbol("N") + + +@dace.program +def _axpy(a: dace.float64[N], b: dace.float64[N], c: dace.float64[N]): + for i in dace.map[0:N]: + c[i] = a[i] * b[i] + c[i] + + +_AXPY_SDFG = _axpy.to_sdfg(simplify=True) + + +def _exp(**kw): + base = dict( + name="axpy", + program=_AXPY_SDFG, + inputs={ + "a": stats.uniform(0.5, 1.0), + "b": stats.uniform(0.5, 1.0), + "c": stats.uniform(0.5, 1.0), + }, + symbols={"N": 64}, + ) + base.update(kw) + return ExperimentConfig(**base) + + +def _has_gpu() -> bool: + """Whether a runnable GPU device is present.""" + import shutil + import subprocess + + for smi, args in (("nvidia-smi", ["-L"]), ("rocm-smi", ["--showid"])): + if shutil.which(smi) is None: + continue + try: + out = subprocess.run( + [smi, *args], capture_output=True, text=True, timeout=15 + ) + except Exception: + continue + if out.returncode == 0 and out.stdout.strip(): + return True + return False + + +requires_gpu = pytest.mark.skipif(not _has_gpu(), reason="no GPU device available") + +#: Time-stepped axpy: the kernel runs T times per invocation (tests loop grouping). +T = dace.symbol("T") + + +@dace.program +def _axpy_loop(a: dace.float64[N], b: dace.float64[N], c: dace.float64[N]): + for _ in range(T): + for i in dace.map[0:N]: + c[i] = a[i] * b[i] + c[i] + + +_AXPY_LOOP_SDFG = _axpy_loop.to_sdfg(simplify=True) + + +def test_precision_roundtrip(): + for key in ("fp16", "fp32", "fp64", "fp32sr", "fp64sr", "mpfr128"): + assert registry.key_of(registry.to_typeclass(key)) == key + + +def test_store_roundtrip(): + db = ResultStore(":memory:") + res = run_error( + ErrorAnalysisConfig(_exp(), precisions=[{"a": "fp32"}], reference="fp64"), + store=db, + ) + rows = db.query(kind="error") + assert len(rows) == 1 + assert rows[0].precision == {"a": "fp32"} + assert rows[0].payload["errors"]["c"]["rel_mean"] == res[0].errors["c"].rel_mean + + +def test_error_decreases_with_precision(): + errs = run_error( + ErrorAnalysisConfig( + _exp(), + precisions=[{"a": "fp16"}, {"a": "fp32"}], + reference="fp64", + n_samples=2, + ), + ) + e16, e32 = errs[0].errors["c"], errs[1].errors["c"] + assert e32.rel_mean < e16.rel_mean + assert e32.rel_mean < 1e-5 + assert e16.linf >= 0 and e16.rel_max >= 0 + + +def test_error_zero_when_candidate_equals_reference(): + errs = run_error( + ErrorAnalysisConfig(_exp(), precisions=[{"a": "fp64"}], reference="fp64") + ) + assert errs[0].errors["c"].linf == 0.0 + + +def test_performance_runs(): + perfs = run_performance( + PerformanceAnalysisConfig( + _exp(), + precisions=[{"a": "fp32", "b": "fp32", "c": "fp32"}], + n_warmup=1, + n_reps=2, + noise={"a": Noise(absolute=0.1, absolute_dist=stats.norm(0.0, 1.0))}, + ), + ) + assert len(perfs) == 1 + assert len(perfs[0].total_times) == 2 + assert all(t > 0 for t in perfs[0].total_times) + + +def test_performance_phase_breakdown(): + perfs = run_performance( + PerformanceAnalysisConfig( + _exp(symbols={"N": 1 << 16}), + precisions=[{"a": "fp32", "b": "fp32", "c": "fp32"}, {}], + n_warmup=1, + n_reps=3, + ), + ) + cast, baseline = perfs + + assert len(cast.total_times) == 3 + assert len(cast.cast_in_times) == 3 + assert len(cast.cast_out_times) == 3 + assert len(cast.kernel_times) == 3 + assert all(t > 0 for t in cast.cast_in_times) + assert all(t >= 0 for t in cast.kernel_times) + # No transfers on CPU. + assert cast.h2d_times == [] + assert cast.d2h_times == [] + # total == sum of phases per rep. + for i in range(3): + expected = cast.cast_in_times[i] + cast.cast_out_times[i] + cast.kernel_times[i] + assert cast.total_times[i] == pytest.approx(expected) + + # No cast: total is kernel. + assert baseline.cast_in_times == [] + assert baseline.cast_out_times == [] + assert baseline.h2d_times == [] + assert baseline.d2h_times == [] + assert baseline.kernel_times == baseline.total_times + assert all(t > 0 for t in baseline.total_times) + + +@requires_gpu +def test_performance_gpu_runs(): + perfs = run_performance( + PerformanceAnalysisConfig( + _exp(symbols={"N": 1 << 16}, target="gpu"), + precisions=[{"a": "fp32", "b": "fp32", "c": "fp32"}], + n_warmup=1, + n_reps=3, + ), + ) + assert len(perfs) == 1 + assert len(perfs[0].total_times) == 3 + assert all(t > 0 for t in perfs[0].total_times) + + +@requires_gpu +def test_performance_gpu_phase_breakdown(): + perfs = run_performance( + PerformanceAnalysisConfig( + _exp(symbols={"N": 1 << 18}, target="gpu"), + precisions=[{"a": "fp32", "b": "fp32", "c": "fp32"}, {}], + n_warmup=1, + n_reps=3, + ), + ) + cast, baseline = perfs + + # Every phase is populated; transfers/cast/kernel are strictly positive. + for phase in ( + cast.h2d_times, + cast.cast_in_times, + cast.kernel_times, + cast.cast_out_times, + cast.d2h_times, + ): + assert len(phase) == 3 + assert all(t >= 0 for t in phase) + assert all(t > 0 for t in cast.h2d_times) + assert all(t > 0 for t in cast.d2h_times) + assert all(t > 0 for t in cast.cast_in_times) + assert all(t > 0 for t in cast.kernel_times) + + # total == sum of phases per rep. + for i in range(3): + expected = ( + cast.h2d_times[i] + + cast.cast_in_times[i] + + cast.kernel_times[i] + + cast.cast_out_times[i] + + cast.d2h_times[i] + ) + assert cast.total_times[i] == pytest.approx(expected) + + # No cast: transfers + kernel only. + assert baseline.cast_in_times == [] + assert baseline.cast_out_times == [] + assert all(t > 0 for t in baseline.h2d_times) + assert all(t > 0 for t in baseline.d2h_times) + assert all(t > 0 for t in baseline.kernel_times) + for i in range(3): + expected = ( + baseline.h2d_times[i] + baseline.kernel_times[i] + baseline.d2h_times[i] + ) + assert baseline.total_times[i] == pytest.approx(expected) + + +@requires_gpu +def test_performance_gpu_loop_kernel_grouped_per_invocation(): + # Kernel fires T times per invocation; must collapse to one value per rep. + perfs = run_performance( + PerformanceAnalysisConfig( + ExperimentConfig( + name="axpy_loop", + program=_AXPY_LOOP_SDFG, + inputs={ + "a": stats.uniform(0.5, 1.0), + "b": stats.uniform(0.5, 1.0), + "c": stats.uniform(0.5, 1.0), + }, + symbols={"N": 1 << 14, "T": 8}, + target="gpu", + ), + precisions=[{"a": "fp32", "b": "fp32", "c": "fp32"}], + n_warmup=1, + n_reps=3, + ), + ) + (r,) = perfs + assert len(r.kernel_times) == 3 + assert len(r.total_times) == 3 + assert all(t > 0 for t in r.kernel_times) + + +@requires_gpu +def test_error_gpu_zero_when_candidate_equals_reference(): + errs = run_error( + ErrorAnalysisConfig( + _exp(target="gpu"), precisions=[{"a": "fp64"}], reference="fp64" + ) + ) + assert errs[0].errors["c"].linf == 0.0 + + +def test_noise_half_specified_raises(): + # A term needs both its magnitude and its distribution; an inert Noise is an error. + with pytest.raises(ValueError, match="half-specified"): + Noise(absolute=1.0) + with pytest.raises(ValueError, match="half-specified"): + Noise(relative_dist=stats.norm(0.0, 1.0)) + with pytest.raises(ValueError, match="perturbs nothing"): + Noise() + + +def test_perturbation_error_scales_with_noise(): + def run(mag): + cfg = PerturbationAnalysisConfig( + _exp(), + noise={"a": Noise(relative=mag, relative_dist=stats.uniform(-1.0, 2.0))}, + n_samples=2, + ) + return run_perturbation(cfg)[0].errors["c"] + + small, big = run(1e-6), run(1e-2) + assert 0.0 < small.rel_mean < big.rel_mean + # axpy: c = a*b + c with positive inputs, so a relative perturbation of `a` + # bounded by mag moves c by at most mag relative. + assert big.rel_max <= 1e-2 * (1.0 + 1e-9) + + +def test_perturbation_one_input_at_a_time(): + noise = Noise(relative=1e-3, relative_dist=stats.uniform(-1.0, 2.0)) + res = run_perturbation( + PerturbationAnalysisConfig(_exp(), noise={"a": noise, "b": noise}) + ) + assert [r.perturbed for r in res] == ["a", "b"] + assert all(r.precision == {} for r in res) + assert all(r.errors["c"].linf > 0.0 for r in res) + + +def test_perturbation_store_roundtrip(): + db = ResultStore(":memory:") + res = run_perturbation( + PerturbationAnalysisConfig( + _exp(), + noise={"a": Noise(relative=1e-3, relative_dist=stats.norm(0.0, 1.0))}, + precisions=[{"a": "fp32"}], + ), + store=db, + ) + rows = db.query(kind="perturbation") + assert len(rows) == 1 + assert rows[0].precision == {"a": "fp32"} + assert rows[0].payload["perturbed"] == "a" + assert rows[0].payload["errors"]["c"]["rel_mean"] == res[0].errors["c"].rel_mean + + +def test_perturbation_unknown_input_raises(): + with pytest.raises(ValueError, match="not a read input"): + run_perturbation( + PerturbationAnalysisConfig( + _exp(), + noise={"z": Noise(absolute=1.0, absolute_dist=stats.norm(0.0, 1.0))}, + ) + ) + + +def test_perturbation_empty_noise_raises(): + with pytest.raises(ValueError, match="at least one input"): + run_perturbation(PerturbationAnalysisConfig(_exp(), noise={})) + + +def test_noise_perturbs_inputs(): + exp = _exp() + sdfg = fresh_sdfg(exp) + a_clean = make_call_args(sdfg, exp, np.random.default_rng(0))["a"] + a_noisy = make_call_args( + sdfg, + exp, + np.random.default_rng(0), + {"a": Noise(absolute=10.0, absolute_dist=stats.norm(0.0, 1.0))}, + )["a"] + assert not np.allclose(a_clean, a_noisy) + + +def test_noise_terms_take_distinct_distributions(): + exp = _exp() + sdfg = fresh_sdfg(exp) + a_clean = make_call_args(sdfg, exp, np.random.default_rng(0))["a"] + noisy = make_call_args( + sdfg, + exp, + np.random.default_rng(0), + {"a": Noise(relative=0.1, relative_dist=stats.uniform(-1.0, 2.0))}, + )["a"] + assert not np.allclose(a_clean, noisy) + assert np.all(np.abs(noisy - a_clean) <= 0.1 * np.abs(a_clean) + 1e-6) + + +def test_error_noise_is_per_analysis(): + errs = run_error( + ErrorAnalysisConfig( + _exp(), + precisions=[{"a": "fp32"}], + reference="fp64", + noise={"a": Noise(relative=1e-2, relative_dist=stats.norm(0.0, 1.0))}, + ), + ) + s = errs[0].errors["c"] + assert np.isfinite(s.rel_mean) and s.linf >= 0.0 + + +def test_is_mpfr(): + assert registry.is_mpfr("mpfr128") + assert not registry.is_mpfr("fp32") + + +def test_group_per_invocation(): + # k executions per invocation are summed into one value per invocation. + assert _group_per_invocation([1.0, 2.0, 3.0, 4.0], 2, "t") == [3.0, 7.0] + # One execution per invocation passes through. + assert _group_per_invocation([1.0, 2.0], 2, "t") == [1.0, 2.0] + assert _group_per_invocation([], 2, "t") == [] + # A data-dependent execution count cannot be attributed to reps: fail loudly. + with pytest.raises(ValueError, match="static per-invocation"): + _group_per_invocation([1.0] * 7, 11, "State s7") + + +def test_overflowing_reference_yields_inf_error_not_nan(): + acc = _new_acc() + _accumulate(acc, np.array([np.inf, 2.0]), np.array([1.0, 2.0])) + s = _finalize(acc) + assert s.rel_max == np.inf + assert not np.isnan(s.rel_mean) + + +def test_error_norms_match_hand_computed_values(): + # e = cand - ref = [0, -4]; r = [3, 4] -> ||r||_1 = 7, ||r||_2 = 5, max|r| = 4 + acc = _new_acc() + _accumulate(acc, np.array([3.0, 4.0]), np.array([3.0, 0.0])) + s = _finalize(acc) + + assert s.l1 == pytest.approx(4.0) + assert s.l2 == pytest.approx(4.0) + assert s.linf == pytest.approx(4.0) + + assert s.l1_norm == pytest.approx(4.0 / 7.0) + assert s.l2_norm == pytest.approx(0.8) + assert s.linf_norm == pytest.approx(1.0) + + assert s.snr == pytest.approx(10.0 * np.log10(25.0 / 16.0)) + # snr is the same quantity as the normalized L2 error, in decibels + assert s.snr == pytest.approx(-20.0 * np.log10(s.l2_norm)) + + +def test_error_norms_concatenate_across_accumulate_calls(): + # Folding two (ref, cand) pairs must equal one fold over their concatenation, + # matching the "all elements over all samples" reduction semantics. + ref1, cand1 = np.array([1.0, 2.0]), np.array([1.5, 2.0]) + ref2, cand2 = np.array([3.0, 4.0]), np.array([3.0, 5.0]) + + split = _new_acc() + _accumulate(split, ref1, cand1) + _accumulate(split, ref2, cand2) + whole = _new_acc() + _accumulate(whole, np.concatenate([ref1, ref2]), np.concatenate([cand1, cand2])) + + a, b = _finalize(split), _finalize(whole) + for f in ("l1", "l2", "linf", "l1_norm", "l2_norm", "linf_norm", "snr"): + assert getattr(a, f) == pytest.approx(getattr(b, f)) + + +def test_error_norms_zero_on_exact_match(): + acc = _new_acc() + _accumulate(acc, np.array([1.0, 2.0, 3.0]), np.array([1.0, 2.0, 3.0])) + s = _finalize(acc) + assert s.l1 == 0.0 and s.l2 == 0.0 and s.linf == 0.0 + assert s.l1_norm == 0.0 and s.l2_norm == 0.0 and s.linf_norm == 0.0 + assert s.snr == np.inf + + +def test_error_norms_infinite_when_reference_signal_is_zero(): + acc = _new_acc() + _accumulate(acc, np.zeros(2), np.array([0.0, 1.0])) + s = _finalize(acc) + assert s.l1_norm == np.inf + assert s.l2_norm == np.inf + assert s.linf_norm == np.inf + assert s.snr == -np.inf + + +def test_error_norms_non_finite_error_gives_minus_inf_snr(): + acc = _new_acc() + _accumulate(acc, np.array([1.0, 2.0]), np.array([np.nan, 2.0])) + s = _finalize(acc) + assert s.l2 == np.inf + assert s.snr == -np.inf + + +def test_error_metrics_persist_to_store(): + db = ResultStore(":memory:") + run_error( + ErrorAnalysisConfig(_exp(), precisions=[{"a": "fp16"}], reference="fp64"), + store=db, + ) + payload = db.query(kind="error")[0].payload["errors"]["c"] + for f in ("l1", "l2", "linf", "l1_norm", "l2_norm", "linf_norm", "snr"): + assert f in payload + + +def test_unknown_target_raises(): + with pytest.raises(ValueError, match="Unknown target"): + apply_target(fresh_sdfg(_exp()), "tpu") + + +def test_callable_init_function_covers_constant_and_fixed_data(): + exp = _exp( + inputs={ + "a": lambda shape, rng: np.full(shape, 2.0), + "b": stats.uniform(0.5, 1.0), + "c": stats.uniform(0.5, 1.0), + } + ) + args = make_call_args(fresh_sdfg(exp), exp, np.random.default_rng(0)) + assert np.allclose(args["a"], 2.0)