diff --git a/.gitignore b/.gitignore index e4713a3..f56e0cd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,24 +1,12 @@ -# Temporary # -############# -.cache -.direnv +__pycache__/ +/.cache/ +/build-install/ +/ulp-report.html +/ulp-report.js +/ulp-report.json -# Compiled source # -################### -*.a -*.dll -*.exe -*.o -*.o.d -*.py[ocd] -*.so -*.mod - -# Logs # -######## -*.log - -# Patches # -########### -*.patch -*.diff +# meson subprojects: track the wrap manifests and the packagefile overlays +# they graft onto upstream checkouts; ignore the checkouts themselves. +/subprojects/* +!/subprojects/*.wrap +!/subprojects/packagefiles/ diff --git a/.spin/cmds.py b/.spin/cmds.py index e918884..18492e6 100644 --- a/.spin/cmds.py +++ b/.spin/cmds.py @@ -1,12 +1,192 @@ +import os import pathlib import sys + import click +import spin.cmds.meson as _meson curdir = pathlib.Path(__file__).parent rootdir = curdir.parent toolsdir = rootdir / "tools" sys.path.insert(0, str(toolsdir)) +PACKAGE = "numpy_sr" + +# Registered in `python/_conftest/hooks.py`; both are skipped unless opted into. +MARKERS = ("slow", "exhaustive") + + +def _markexpr(markers): + """`--marker slow,exhaustive` -> pytest's `-m "slow or exhaustive"`. + + The conftest only arms its skip-the-marked gate while pytest's mark + expression is empty, so `all` hands over a tautology: it disarms the gate + without narrowing the selection.""" + names = list( + dict.fromkeys(n for spec in markers for n in spec.split(",") if n.strip()) + ) + choices = (*MARKERS, "all") + unknown = [n for n in names if n not in choices] + if unknown: + raise click.BadParameter( + f"unknown marker(s): {', '.join(unknown)}; pick from {', '.join(choices)}" + ) + if not names: + return None + if "all" in names: + return f"{MARKERS[0]} or not {MARKERS[0]}" + return " or ".join(names) + + +def _is_path(arg): + return not arg.startswith("-") and ( + arg.endswith(".py") or "::" in arg or os.sep in arg + ) + + +def _as_path(arg, build_dir): + """`numpy_sr.tests.test_trig[::node]` -> its installed file. A `--pyargs` + module selection loads conftest.py too late to register --target & co, + so a selection has to reach pytest as a path. Returns None if unresolved.""" + if not (arg == PACKAGE or arg.startswith(f"{PACKAGE}.")): + return None + module, sep, node = arg.partition("::") + base = pathlib.Path(_meson._get_site_packages(build_dir), *module.split(".")) + for candidate in (base.with_suffix(".py"), base): + if candidate.exists(): + return f"{candidate}{sep}{node}" + return None + + +def _selection(pytest_args, tests, build_dir): + """Resolve test selections to paths and re-add the `--pyargs numpy_sr` that + upstream drops as soon as any pytest argument is given -- without it pytest + collects the whole of site-packages.""" + args = tuple(_as_path(a, build_dir) or a for a in pytest_args) + resolved = _as_path(tests, build_dir) if tests else None + if resolved: + args, tests = (*args, resolved), None + if tests or any(_is_path(a) for a in args): + return args, tests + return ("--pyargs", PACKAGE, *args), tests + + +@click.command( + help=f"""๐Ÿ”ง Run tests + + Wraps `spin.cmds.meson.test`, exposing the suite's own pytest options + directly and keeping collection rooted at `--pyargs {PACKAGE}`: + + \b + spin test --target avx2,sse4 + spin test --accuracy high --seed 12345 + spin test --collect-worst + spin test --marker slow + + Plain pytest options still go after `--`: + + \b + spin test --target sse2 -- -k trig -x + """ +) +@click.argument("pytest_args", nargs=-1) +@click.option( + "--target", + metavar="NAMES", + help="Restrict ISA targets; comma-separated (default: every supported " + "target except references).", +) +@click.option( + "--accuracy", + metavar="NAMES", + help="Restrict accuracy profiles; comma-separated (default: high,low).", +) +@click.option( + "--seed", + metavar="INT", + help="Root seed for every random stream (default: drawn fresh, echoed at " + "the end of the run).", +) +@click.option( + "--marker", + "-m", + "markers", + metavar="NAMES", + multiple=True, + help="Opt into marker-gated tests, comma-separated or repeated: " + f"{', '.join(MARKERS)}, or `all` to add them to the regular run. Naming " + "one runs only its tests (default: neither, both are skipped).", +) +@click.option( + "--collect-worst", + is_flag=True, + default=False, + help="Fold measured inputs into the worst-case book instead of only " + "replaying it as a gate.", +) +@click.option( + "-j", + "n_jobs", + metavar="N_JOBS", + default="1", + help="Number of xdist workers; `auto` for all cores.", +) +@click.option( + "--tests", + "-t", + metavar="TESTS", + help=f"Module, class or function to run, e.g. {PACKAGE}.tests.test_trig.", +) +@click.option("--verbose", "-v", is_flag=True, default=False) +@_meson.build_option +@_meson.build_dir_option +@click.pass_context +def test( + ctx, + *, + pytest_args, + target, + accuracy, + seed, + markers, + collect_worst, + n_jobs, + tests, + verbose, + build, + build_dir, +): + suite_args = [] + for name, value in ( + ("--target", target), + ("--accuracy", accuracy), + ("--seed", seed), + ): + if value is not None: + suite_args += [name, value] + if collect_worst: + suite_args.append("--collect-worst") + + markexpr = _markexpr(markers) + if markexpr: + if any(a == "-m" or a.startswith("-m") for a in pytest_args): + raise click.UsageError( + "--marker clashes with the `-m` handed to pytest after `--`; " + "keep just one of them." + ) + suite_args += ["-m", markexpr] + + pytest_args, tests = _selection(tuple(suite_args) + pytest_args, tests, build_dir) + ctx.invoke( + _meson.test, + pytest_args=pytest_args, + n_jobs=n_jobs, + tests=tests, + verbose=verbose, + build=build, + build_dir=build_dir, + ) + @click.command(help="Generate sollya c++/python based files") @click.option("-f", "--force", is_flag=True, help="Force regenerate all files") @@ -14,3 +194,72 @@ def sollya(*, force): import sollya # type: ignore[import] sollya.main(force=force) + + +# Markers in tools/ulp_report/index.html -> what replaces them in the bundle. +# Tabulator stays a CDN ' +_SIDECAR = '' +_VIEWER = '' + + +def _replace_once(html, marker, repl): + if html.count(marker) != 1: + raise SystemExit(f"index.html: expected exactly one {marker!r}") + return html.replace(marker, repl) + + +def _bundle(json_text): + viewdir = toolsdir / "ulp_report" + html = (viewdir / "index.html").read_text() + # "\n{(viewdir / 'viewer.css').read_text()}" + ) + html = _replace_once( + html, + _DATA_SLOT, + f'', + ) + html = _replace_once(html, _SIDECAR, "") + return _replace_once( + html, _VIEWER, f"" + ) + + +@click.command(help="Bundle the ULP JSON report + viewer into one self-contained HTML") +@click.argument("json_path", required=False, type=click.Path(path_type=pathlib.Path)) +@click.option( + "-o", + "--output", + type=click.Path(path_type=pathlib.Path), + default=pathlib.Path("ulp-report.html"), + show_default=True, + help="Output HTML path", +) +@click.option( + "--open/--no-open", + "show", + default=True, + show_default=True, + help="Open the result in the default browser (--no-open for CI)", +) +def ulp_report(json_path, output, show): + import webbrowser + + json_path = json_path or rootdir / "ulp-report.json" + if not json_path.exists(): + raise SystemExit( + f"{json_path}: not found. Run the test suite first, or pass a report path." + ) + json_text = json_path.read_text() + output.write_text(_bundle(json_text)) + print(f"{output} ({output.stat().st_size:,} bytes)") + # Sidecar for the live index.html over file://, where fetch() is CORS-blocked. + sidecar = json_path.with_suffix(".js") + sidecar.write_text(f"window.NPSR_ULP_DATA = {json_text};\n") + print(f"{sidecar} ({sidecar.stat().st_size:,} bytes)") + if show: + webbrowser.open(output.resolve().as_uri()) diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..cc59917 --- /dev/null +++ b/meson.build @@ -0,0 +1,21 @@ + +project( + 'npsr', + 'cpp', + version: '0.1.0', + meson_version: '>= 1.3.0', + license: 'BSD-3-Clause', + default_options: [ + 'cpp_std=c++17', + 'buildtype=release', + ], +) + +npsr_inc = include_directories('.') + +# Builds the `numpy_sr` Python package: one foreach_target TU per operation +# (python/_numpy_sr/.cpp) linked into a single pybind11 extension that +# exposes a submodule per runnable Highway target. +if get_option('tests').disable_auto_if(meson.is_subproject()).allowed() + subdir('python') +endif diff --git a/meson_options.txt b/meson_options.txt new file mode 100644 index 0000000..07887e7 --- /dev/null +++ b/meson_options.txt @@ -0,0 +1,2 @@ + +option('tests', type: 'feature', value: 'auto', description: 'Build per-operation shared libraries for the pytest harness') diff --git a/npsr/lut-inl.h b/npsr/lut-inl.h index b9e1cdf..97905ee 100644 --- a/npsr/lut-inl.h +++ b/npsr/lut-inl.h @@ -5,6 +5,7 @@ #define NPSR_LUT_INL_H_ #endif +#include #include #include "npsr/hwy.h" diff --git a/npsr/precise.h b/npsr/precise.h index ec4e617..93ee071 100644 --- a/npsr/precise.h +++ b/npsr/precise.h @@ -68,6 +68,13 @@ class FPExceptions { #else static constexpr auto kUnderflow = 0; #endif +#ifdef FE_INEXACT + static constexpr auto kInexact = FE_INEXACT; +#else + static constexpr auto kInexact = 0; +#endif + // kInexact is deliberately excluded: nearly every rounded operation raises + // it, so it carries no error signal for Raise()/Test() callers. static constexpr auto kAll = kInvalid | kDivByZero | kOverflow | kUnderflow; void Raise(int errors) noexcept { mask_ |= errors; } diff --git a/pyproject.toml b/pyproject.toml index 1338374..4e3da23 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,10 +7,60 @@ version = "1.0.0.dev0" description = "NumPy SIMD Routines" requires-python = ">=3.12" +[project.optional-dependencies] +# `spin test -j ` needs xdist; the plugin merges the workers' results itself. +test = ["pytest", "pytest-xdist"] + +[tool.ruff] +# Python style + lint: run `ruff check` and `ruff format`. 88 cols matches the +# width the current sources are written at. Build outputs, vendored subprojects +# and throwaway *_backup.py scratch files are excluded. +line-length = 88 +extend-exclude = ["build", "build-install", "subprojects", "*_backup.py"] + +[tool.ruff.lint] +# Bugs (F/B), style (E/W), imports (I), modernization (UP), simplify (SIM), +# comprehensions (C4), and ruff's own checks (RUF). +select = ["E", "W", "F", "I", "B", "C4", "UP", "SIM", "RUF"] + +[tool.pyright] +# Static type check for the package (run `pyright`); analyses the source tree +# without a build. The package installs as `numpy_sr` (meson subdir) but the +# source folder is `python/`, so absolute self-imports (`import numpy_sr`) are +# resolved against the built package under build-install/. That path carries the +# interpreter version, so bump the list when the Python minor changes. +include = ["python"] +exclude = [ + "**/__pycache__", + "**/*_backup.py", + "build", + "build-install", + "subprojects", +] +extraPaths = [ + "build-install/usr/lib/python3.12/site-packages", + "build-install/usr/lib/python3.13/site-packages", +] +pythonVersion = "3.12" +# The _numpy_sr extension is compiled-only and resolved via its .pyi stub; +# silence the "stub found without a source module" note that follows. +reportMissingModuleSource = "none" + [tool.spin] package = "numpy_sr" [tool.spin.commands] +"Build" = [ + "spin.cmds.meson.build", +] "Codegen" = [ ".spin/cmds.py:sollya", ] +"Testing" = [ + ".spin/cmds.py:test", + ".spin/cmds.py:ulp_report", +] +"Environments" = [ + "spin.cmds.meson.python", + "spin.cmds.meson.ipython", +] diff --git a/python/__init__.py b/python/__init__.py new file mode 100644 index 0000000..9bd0e62 --- /dev/null +++ b/python/__init__.py @@ -0,0 +1,21 @@ +from __future__ import annotations + +from . import argred, ieee754, testing +from ._numpy_sr import Accuracy, targets, ulp_error, ulp_error32 + +# Context manager over the flags-only `_numpy_sr.FPEnv`, which it shadows here. +from .fpenv import FPEnv + +mpfr = targets["mpfr"] + +__all__ = [ + "Accuracy", + "FPEnv", + "argred", + "ieee754", + "mpfr", + "targets", + "testing", + "ulp_error", + "ulp_error32", +] diff --git a/python/_conftest/__init__.py b/python/_conftest/__init__.py new file mode 100644 index 0000000..4ee9868 --- /dev/null +++ b/python/_conftest/__init__.py @@ -0,0 +1 @@ +"""Internal helpers for `numpy_sr/conftest.py`, which re-exports them all.""" diff --git a/python/_conftest/book.py b/python/_conftest/book.py new file mode 100644 index 0000000..e842295 --- /dev/null +++ b/python/_conftest/book.py @@ -0,0 +1,274 @@ +"""The frozen worst-case book: `WorstBook`, folded during a run.""" + +import pathlib +import warnings +from typing import ClassVar + +import numpy as np + + +class WorstBook: + """A capped CSV per (op, accuracy, dtype, fma, source) under tests/worst/.""" + + # dtype -> (bit-pattern uint, float, rows kept per source file). Headerless + # `,`: col0 is decimal, the only exact form np.loadtxt parses in C. + # f32 keeps as many rows as f64: only test_exhaustive books it, and a sweep of + # every finite input deserves the wider cap that sampling did not. + _SPEC: ClassVar[dict[str, tuple[type, type, int]]] = { + "float32": (np.uint32, np.float32, 10_000), + "float64": (np.uint64, np.float64, 10_000), + } + _EPS = 1e-5 # stored errors round to 6 decimals; reruns must stay idempotent + + def __init__(self, root: pathlib.Path) -> None: + self.root = root + # key -> (uint bits ascending, |ulp| err) of the live top cases + self._book: dict[str, tuple[np.ndarray, np.ndarray]] = {} + self._base_max: dict[str, float] = {} # per-key max ULP as committed + self._dirty: set[str] = set() + self._touched: set[str] = set() # keys this process has already measured + self.discoveries: dict[str, dict] = {} + self.restamped: dict[str, tuple[float, float]] = {} + + @staticmethod + def key( + operation: str, accuracy: str, dtype: str, have_fma: bool, source: str + ) -> str: + # `parse` splits on "_", so `operation` must not contain one. + fma = "fma" if have_fma else "nonfma" + return f"{operation}_{accuracy.lower()}_{dtype}_{fma}_{source}" + + @staticmethod + def parse(key: str) -> tuple[str, str, str, bool, str]: + """`key()` inverted: `(operation, accuracy, dtype, have_fma, source)`.""" + operation, accuracy, dtype, fma, source = key.split("_", 4) + return operation, accuracy, dtype, fma == "fma", source + + def all_keys(self) -> list[str]: + """Every book on disk, sorted; the gate parametrizes over these.""" + if not self.root.is_dir(): + return [] + return sorted( + f"{op_dir.name}_{path.stem}" + for op_dir in self.root.iterdir() + if op_dir.is_dir() + for path in op_dir.glob("*.csv") + ) + + @staticmethod + def source(label: str, *drop: str) -> str: + """`file.py::test_simple[SSE4-Low-_basic-float32-sin]` -> `test_simple[_basic]`. + + Each parametrization books separately, or one cap is shared by unrelated + data sets. `drop` removes the id segments the key already carries โ€” the + target most of all, since every non-FMA target would file the same rows. + """ + base, sep, bracket = label.split("::")[-1].partition("[") + dropped = {d.lower() for d in drop} + kept = ( + [s for s in bracket.rstrip("]").split("-") if s.lower() not in dropped] + if sep + else [] + ) + return f"{base}[{'-'.join(kept)}]" if kept else base + + def _spec(self, key: str) -> tuple: + """`(uint, float, cap)` of `____`.""" + return self._SPEC[self.parse(key)[2]] + + def _path(self, key: str) -> pathlib.Path: + op, rest = key.split("_", 1) # the op names its dir + return self.root / op / f"{rest}.csv" + + def _row_dtype(self, key: str) -> np.dtype: + return np.dtype([("bits", self._spec(key)[0]), ("err", "f8")]) + + @staticmethod + def _top(bits: np.ndarray, errs: np.ndarray, cap: int) -> tuple: + """The `cap` worst rows, bits ascending (the `_lookup` invariant).""" + if bits.size > cap: + keep = np.argpartition(-errs, cap)[:cap] + bits, errs = bits[keep], errs[keep] + order = np.argsort(bits) + return bits[order], errs[order] + + def _load(self, key: str) -> tuple[np.ndarray, np.ndarray]: + if key not in self._book: + u, _, cap = self._spec(key) + path = self._path(key) + if path.exists(): + with warnings.catch_warnings(): # an empty file is legitimate + warnings.simplefilter("ignore") + rec = np.loadtxt( + path, dtype=self._row_dtype(key), delimiter=",", ndmin=1 + ) + pair = self._top(rec["bits"], rec["err"], cap) + else: + pair = (np.array([], u), np.array([], np.float64)) + self._book[key] = pair + self._base_max[key] = float(pair[1].max(initial=0.0)) + return self._book[key] + + @staticmethod + def _member(bits: np.ndarray, other: np.ndarray) -> np.ndarray: + """Mask of `bits` that appear in the ascending `other`.""" + if other.size == 0: + return np.zeros(bits.shape, bool) + pos = np.clip(np.searchsorted(other, bits), 0, other.size - 1) + return other[pos] == bits + + @classmethod + def _lookup(cls, bits0: np.ndarray, err0: np.ndarray, query: np.ndarray): + """`err0` for each `query` bit in the ascending `bits0`, else -inf.""" + if bits0.size == 0: + return np.full(query.shape, -np.inf) + pos = np.clip(np.searchsorted(bits0, query), 0, bits0.size - 1) + return np.where(bits0[pos] == query, err0[pos], -np.inf) + + @staticmethod + def _dedupe(bits: np.ndarray, errs: np.ndarray) -> tuple: + """Worst error per input, bits ascending.""" + order = np.argsort(bits, kind="stable") + bits, first = np.unique(bits[order], return_index=True) + return bits, np.maximum.reduceat(errs[order], first) + + def _stale(self, key: str) -> bool: + """True once: `key` still holds disk errors, not this run's measurements.""" + first = key not in self._touched + self._touched.add(key) + return first + + def fold(self, key: str, x: np.ndarray, err: np.ndarray) -> None: + """Merge measured (x, |ulp|) into `key`'s book; `x` must be its dtype.""" + x = np.asarray(x) + err = np.asarray(err, dtype=np.float64) + keep = np.isfinite(x) & np.isfinite(err) + if not keep.any(): + return + u, _, cap = self._spec(key) + nb, ne = self._dedupe(np.ascontiguousarray(x[keep]).view(u), err[keep]) + # bound the merge: this batch can contribute at most its own worst `cap` + nb, ne = self._top(nb, ne, cap) + + b0, e0 = self._load(key) + # A live measurement supersedes the stored error for the same input: a + # kernel fix must be able to lower it. Later folds in one run take the + # max instead, so several targets sharing a key keep the worst. + drop = self._member(b0, nb) if self._stale(key) else np.zeros(b0.shape, bool) + bits, errs = self._dedupe( + np.concatenate([b0[~drop], nb]), np.concatenate([e0[~drop], ne]) + ) + bits, errs = self._top(bits, errs, cap) + + improved = errs > self._lookup(b0, e0, bits) + self._EPS + n_new = int(improved.sum()) + self._book[key] = (bits, errs) + # `or` short-circuits: allclose never sees a shape mismatch + if not np.array_equal(bits, b0) or not np.allclose( + errs, e0, rtol=0, atol=self._EPS + ): + self._dirty.add(key) + if n_new: + w = int(np.argmax(np.where(improved, errs, -np.inf))) + d = self.discoveries.setdefault( + key, + {"n_new": 0, "base_max": self._base_max[key], "err": -np.inf}, + ) + d["n_new"] += n_new + if errs[w] > d["err"]: + d["err"], d["bits"] = float(errs[w]), int(bits[w]) + + def refresh(self, key: str, x: np.ndarray, err: np.ndarray) -> None: + """Re-stamp `key`'s stored errors from a live measurement of its inputs. + + Keyed, never slice-wide: two sources may hold the same input, and each + book must record what its own replay measured. Without this the errors + only ever ratchet up: a kernel fix would leave every book at its + historical high-water mark, and `_top` would evict genuinely hard + inputs in favour of ones that are no longer hard. + """ + u, _, _ = self._spec(key) + keep = np.isfinite(x) & np.isfinite(err) + qb, qe = self._dedupe( + np.ascontiguousarray(np.asarray(x)[keep]).view(u), + np.asarray(err, dtype=np.float64)[keep], + ) + bits, errs = self._load(key) + measured = self._lookup(qb, qe, bits) + hit = measured > -np.inf + # first touch this run supersedes disk; later targets keep the worst + new = np.where( + hit if self._stale(key) else hit & (measured > errs), measured, errs + ) + if np.allclose(new, errs, rtol=0, atol=self._EPS): + return + self._book[key] = (bits, new) + self._dirty.add(key) + was = self._base_max[key] + self._base_max[key] = now = float(new.max(initial=0.0)) + self.restamped[key] = (self.restamped.get(key, (was, now))[0], now) + + def replay(self, key: str) -> np.ndarray: + """`key`'s stored inputs as its own dtype, bits ascending.""" + _, f, _ = self._spec(key) + return self._load(key)[0].view(f) + + def flush(self) -> list[str]: + """Rewrite every changed book (worst first); return the keys written.""" + written = sorted(self._dirty) + for key in written: + bits, errs = self._book[key] + order = np.argsort(-errs, kind="stable") + rec = np.empty(order.size, dtype=self._row_dtype(key)) + rec["bits"], rec["err"] = bits[order], errs[order] + path = self._path(key) + path.parent.mkdir(parents=True, exist_ok=True) + np.savetxt(path, rec, fmt="%d,%.6f") + self._dirty.clear() + return written + + def deltas(self) -> dict[str, tuple[bytes, bytes]]: + """Raw (bits, err) buffers of every book this process changed.""" + return { + key: (self._book[key][0].tobytes(), self._book[key][1].tobytes()) + for key in self._dirty + } + + def absorb( + self, + deltas: dict[str, tuple[bytes, bytes]], + restamped: dict[str, tuple[float, float]] | None = None, + ) -> None: + """Fold another process's `deltas()` in; discoveries recount from disk.""" + for key, (was, now) in (restamped or {}).items(): + self._load(key) # seeds _base_max from disk before it is superseded + prev = self.restamped.get(key) + self.restamped[key] = ( + (was, now) if prev is None else (prev[0], max(prev[1], now)) + ) + # a discovery on the controller ranks against the re-measured max + self._base_max[key] = self.restamped[key][1] + for key, (bits, errs) in deltas.items(): + u, f, _ = self._spec(key) + self.fold( + key, + np.frombuffer(bits, dtype=u).view(f), + np.frombuffer(errs, dtype=np.float64), + ) + + def report_lines(self) -> list[str]: + """One line per book this run re-stamped or added a case to.""" + out = [ + f"{key}: re-measured, max {was:.4f} -> {now:.4f} ULP" + for key, (was, now) in sorted(self.restamped.items()) + ] + for key, d in sorted(self.discoveries.items()): + u, f, _ = self._spec(key) + x = float(np.array(d["bits"], u).view(f)) + base, now = d["base_max"], float(self._book[key][1].max(initial=0.0)) + # book max before -> after; the new case need not be the book's worst + out.append( + f"{key}: +{d['n_new']} case(s), max {base:.4f} -> {now:.4f} ULP " + f"(worst new {d['err']:.4f} at x={x.hex()})" + ) + return out diff --git a/python/_conftest/fixtures.py b/python/_conftest/fixtures.py new file mode 100644 index 0000000..2a34159 --- /dev/null +++ b/python/_conftest/fixtures.py @@ -0,0 +1,149 @@ +"""Fixtures and per-target/-accuracy parametrization.""" + +import pathlib +import zlib +from collections.abc import ValuesView + +import numpy as np +import numpy_sr as sr +import pytest + +from .book import WorstBook +from .options import selected_accuracies, selected_targets +from .stash import BOOK_KEY, REPORT_KEY +from .targets import ACCURACIES, TARGETS +from .xdist import resolve_seed + + +@pytest.fixture(scope="session") +def targets(request: pytest.FixtureRequest) -> ValuesView: + """Every selected target at once, for tests that compare backends.""" + return selected_targets(request.config).values() + + +@pytest.fixture(scope="session") +def accuracies(request: pytest.FixtureRequest) -> list: + """Every selected accuracy at once, for tests that sweep them internally.""" + return selected_accuracies(request.config) + + +@pytest.fixture(params=list(TARGETS.values()), ids=list(TARGETS)) +def target(request: pytest.FixtureRequest): + """The backend under test; CLI selection deselects the rest at collection.""" + return request.param + + +@pytest.fixture( + params=list(ACCURACIES.values()), ids=[a.name for a in ACCURACIES.values()] +) +def accuracy(request: pytest.FixtureRequest): + """The `sr.Accuracy` under test; unselected ones deselected at collection.""" + return request.param + + +@pytest.fixture +def rng(request: pytest.FixtureRequest) -> np.random.Generator: + """Per-test PRNG; the salt drops the target id: all targets, same data.""" + node = request.node + base, sep, bracket = node.name.partition("[") + segments = bracket.rstrip("]").split("-") if sep else [] + if "target" in request.fixturenames: + target = request.getfixturevalue("target") + target_id = target.__name__.rsplit(".", 1)[-1] + segments = [s for s in segments if s != target_id] + parts = [node.path.name, base, *segments] + seed = resolve_seed(request.config) + return np.random.default_rng([seed, zlib.crc32("|".join(parts).encode())]) + + +@pytest.fixture +def worst_book(request: pytest.FixtureRequest) -> "WorstBook": + """The session's frozen worst-case book, for tests that replay it.""" + return request.config.stash[BOOK_KEY] + + +@pytest.fixture +def report_ulp(request: pytest.FixtureRequest): + """`assert_max_ulp` bound to the session's report and worst-case book.""" + nodeid = request.node.nodeid + records = request.config.stash[REPORT_KEY] + book = request.config.stash[BOOK_KEY] + collect = request.config.getoption("--collect-worst") + path, _, rest = nodeid.partition("::") + label = f"{pathlib.Path(path).name}::{rest}" if rest else path + + # operation/target/accuracy default to these; only what the test pulled in + global_kw = { + a: request.getfixturevalue(a) if a in request.fixturenames else None + for a in ("operation", "target", "accuracy") + } + + def sink(record: dict) -> None: + records.append(record) + + def check( + x, + computed, + ref, + maxulp, + *, + store_worst=None, + refresh_key=None, + operation=None, + target=None, + accuracy=None, + ): + # Explicit rather than **kw: a stray keyword must fail, not vanish. + __tracebackhide__ = True + given = {"operation": operation, "target": target, "accuracy": accuracy} + kw = {} + for k, v in given.items(): + kw[k] = v = global_kw[k] if v is None else v + if v is None: + raise TypeError( + f"{nodeid}: no `{k}` fixture in scope โ€” " + f"pass {k}=... to report_ulp()" + ) + + target = kw["target"] + target_name = getattr(target, "__name__", str(target)).rsplit(".", 1)[-1] + is_reference = getattr(target, "IS_REFERENCE", False) + a = np.asarray(computed) + # default: f64 only. Every f32 input is covered by test_exhaustive, which + # opts back in with store_worst=True; sampling f32 elsewhere burns the cap. + store = a.dtype == np.float64 if store_worst is None else store_worst + # Every write to the book goes through here, so `is_reference` is enforced + # once: only a non-reference target's error belongs in the book. + if collect and not is_reference and (store or refresh_key is not None): + xb = np.broadcast_to(np.asarray(x).astype(a.dtype, copy=False), a.shape) + err = np.abs(sr.testing.ulp_error(a, ref)) + xr, er = xb.ravel(), np.asarray(err).ravel() + if store: + axes = (kw["operation"], kw["accuracy"].name, str(a.dtype), target_name) + book.fold( + WorstBook.key( + *axes[:3], target.HAVE_FMA, WorstBook.source(label, *axes) + ), + xr, + er, + ) + if refresh_key is not None: + # a replay measures every stored input: stamp the book with what + # this kernel actually does, up or down + book.refresh(refresh_key, xr, er) + + meta = { + "x": x, + "label": label, + "target": target_name, + "accuracy": kw["accuracy"].name, + "operation": kw["operation"], + } + # A reference target is not the oracle (MPFR is): recorded, not asserted. + if is_reference: + record = sr.testing.measure(computed, ref, maxulp, **meta) + sink(record) + return record + return sr.testing.assert_max_ulp(computed, ref, maxulp, sink=sink, **meta) + + return check diff --git a/python/_conftest/hooks.py b/python/_conftest/hooks.py new file mode 100644 index 0000000..834d3bb --- /dev/null +++ b/python/_conftest/hooks.py @@ -0,0 +1,65 @@ +"""Session setup and collection hooks: stash, markers, selection, ordering.""" + +import pathlib + +import pytest + +from .book import WorstBook +from .options import selected_accuracies, selected_targets +from .stash import BOOK_KEY, REPORT_KEY +from .xdist import warn_unmerged + + +def _book_dir(config: pytest.Config) -> pathlib.Path: + """The committed `python/tests/worst`, not the build-install copy.""" + # `spin test` roots pytest in site-packages; walk up to the source ancestor. + start = pathlib.Path(config.rootpath) + for p in (start, *start.parents): + if (p / "python" / "tests").is_dir(): + return p / "python" / "tests" / "worst" + return start / "python" / "tests" / "worst" + + +def pytest_configure(config: pytest.Config) -> None: + config.stash[REPORT_KEY] = [] + config.stash[BOOK_KEY] = WorstBook(_book_dir(config)) + warn_unmerged(config) + config.addinivalue_line( + "markers", + "exhaustive: full f32 sweep, minutes per op; run with -m exhaustive", + ) + config.addinivalue_line( + "markers", + "slow: heavy worst-case generation, seconds-to-minutes; run with -m slow", + ) + + +def pytest_collection_modifyitems( + config: pytest.Config, items: list[pytest.Item] +) -> None: + # The fixtures parametrize over the whole registry; deselect (not skip, which + # the full reference set would flood) whatever the CLI did not select. + sel_targets = set(selected_targets(config).values()) + sel_accuracies = set(selected_accuracies(config)) + kept, dropped = [], [] + for item in items: + callspec = getattr(item, "callspec", None) + params = callspec.params if callspec else {} + keep = ("target" not in params or params["target"] in sel_targets) and ( + "accuracy" not in params or params["accuracy"] in sel_accuracies + ) + (kept if keep else dropped).append(item) + if dropped: + config.hook.pytest_deselected(items=dropped) + items[:] = kept + + # The frozen gate guards the whole suite: hoist it whatever file it lives in. + items.sort(key=lambda it: getattr(it, "originalname", it.name) != "test_worst_book") + if config.option.markexpr: + return + for item in items: + for mark in ["exhaustive", "slow"]: + if mark in item.keywords: + item.add_marker( + pytest.mark.skip(reason=f"{mark}; opt in with -m {mark}") + ) diff --git a/python/_conftest/options.py b/python/_conftest/options.py new file mode 100644 index 0000000..0fc3cfa --- /dev/null +++ b/python/_conftest/options.py @@ -0,0 +1,74 @@ +"""The `--target`/`--accuracy`/`--seed`/`--collect-worst` selection axes.""" + +import argparse +from collections.abc import Callable +from typing import cast + +import numpy as np +import pytest + +from .targets import ACCURACIES, DEFAULT_ACCURACIES, DEFAULT_TARGETS, TARGETS + + +def _csv_choice(choices: dict, label: str) -> Callable[[str], list]: + """An argparse `type` splitting on commas and validating each name.""" + + def parse(arg: str) -> list: + out = [] + for raw in arg.split(","): + key = raw.strip().lower() + if not key: + continue + if key not in choices: + raise argparse.ArgumentTypeError( + f"unknown {label} {raw.strip()!r}; choose from {', '.join(choices)}" + ) + out.append(choices[key]) + return out + + return parse + + +def pytest_addoption(parser: pytest.Parser) -> None: + parser.addoption( + "--target", + default=DEFAULT_TARGETS, + metavar="NAMES", + type=_csv_choice({n.lower(): n for n in TARGETS}, "target"), + help="restrict targets; comma-separated " + f"(default: {', '.join(DEFAULT_TARGETS)})", + ) + parser.addoption( + "--accuracy", + default=DEFAULT_ACCURACIES, + metavar="NAMES", + type=_csv_choice(ACCURACIES, "accuracy"), + help="restrict accuracies; comma-separated " + f"(default: {', '.join(a.name.lower() for a in DEFAULT_ACCURACIES)})", + ) + parser.addoption( + "--seed", + default=int(np.random.SeedSequence().generate_state(1)[0]), + metavar="INT", + type=lambda s: int(s, 0), + help="root seed for every random stream " + "(default: drawn fresh, echoed at the end of the run)", + ) + parser.addoption( + "--collect-worst", + action="store_true", + default=False, + help="fold measured inputs into the worst-case book " + "(off by default: the book is replayed as a gate but never rewritten)", + ) + + +def selected_targets(config: pytest.Config) -> dict: + # the dict drops repeats from e.g. --target avx2,avx2 + selected = cast(list, config.getoption("--target")) + return {n: TARGETS[n] for n in selected} + + +def selected_accuracies(config: pytest.Config) -> list: + selected = cast(list, config.getoption("--accuracy")) + return list(dict.fromkeys(selected)) diff --git a/python/_conftest/report.py b/python/_conftest/report.py new file mode 100644 index 0000000..8349710 --- /dev/null +++ b/python/_conftest/report.py @@ -0,0 +1,214 @@ +"""Run reporting and the ULP JSON payload for the HTML viewer.""" + +import json +import math +import time +from functools import partial + +import pytest + +from .options import selected_accuracies, selected_targets +from .stash import BOOK_KEY, REPORT_KEY +from .targets import TARGETS +from .xdist import resolve_seed, worker_output + +_OFFENDER_LIMIT = 10 +_SENTINEL = 1e300 # finite stand-in for non-finite sort keys +_TOKENS = {"nan": math.nan, "inf": math.inf, "-inf": -math.inf} + + +def pytest_report_header(config: pytest.Config) -> list: + """The selection echoed at the top and bottom of every run.""" + targets = [ + f"{name} (reference)" if getattr(target, "IS_REFERENCE", False) else name + for name, target in selected_targets(config).items() + ] + accuracies = [a.name for a in selected_accuracies(config)] + return [ + f"targets: {', '.join(targets)}", + f"accuracies: {', '.join(accuracies)}", + f"seed: {resolve_seed(config)}", + ] + + +def pytest_terminal_summary( + terminalreporter, exitstatus: int, config: pytest.Config +) -> None: + terminalreporter.write_sep("-", "selection") + for line in pytest_report_header(config): + terminalreporter.write_line(line) + + book = config.stash.get(BOOK_KEY, None) + lines = book.report_lines() if book is not None else [] + if lines: + terminalreporter.write_sep("-", "worst-case discoveries") + for line in lines: + terminalreporter.write_line(line) + + +def pytest_sessionfinish(session: pytest.Session, exitstatus: int) -> None: + config = session.config + book = config.stash.get(BOOK_KEY, None) + records = config.stash.get(REPORT_KEY, []) + config.stash[REPORT_KEY] = [] # drained; guard against a double dump + + out = worker_output(config) + if out is not None: + # The controller owns every write; ship the results and touch no files. + out["npsr_records"] = records + out["npsr_book"] = book.deltas() if book is not None else {} + out["npsr_restamped"] = book.restamped if book is not None else {} + return + + if book is not None: + book.flush() # persist any book that gained a worse case + + if records: + head = {"seed": resolve_seed(config), "created": time.time()} + text = json.dumps(head | _serialize_records(records), indent=2) + (config.rootpath / "ulp-report.json").write_text(text) + + +def _pnum(v: float | str) -> float: + """Back to a float for computation; `testing._num` tokenizes non-finite.""" + return _TOKENS[v] if isinstance(v, str) else float(v) + + +def _mag(v: float | str) -> float: + """`|v|` as a rank; non-finite (nan included) outranks any finite error.""" + f = _pnum(v) + return abs(f) if math.isfinite(f) else math.inf + + +def _fmt(v: float | str, fmt) -> str: + """`fmt` applied, except to non-finite values: those keep their token.""" + f = _pnum(v) + if math.isnan(f): + return "nan" + if math.isinf(f): + return "inf" if f > 0 else "-inf" + return fmt(f) + + +_fmt_ulp = partial(_fmt, fmt="{:.4g}".format) # ULP magnitude, ~4 sig figs +_fmt_val = partial(_fmt, fmt=repr) # operand: shortest round-tripping decimal + + +def _finite(f: float) -> float: + """A JSON-safe sort key: non-finite collapses to a large finite sentinel.""" + return f if math.isfinite(f) else math.copysign(_SENTINEL, f) + + +def _ratio(max_error: float | str, maxulp: float) -> float: + """Bar length + worst-first rank; non-finite error ranks above any bound.""" + m = _mag(max_error) + if math.isinf(m): + return math.inf + return 0.0 if math.isinf(maxulp) else m / maxulp + + +def _test_name(label: str) -> str: + """`file.py::test_luck[High-...]` -> `test_luck` (the group carries params).""" + return label.split("::")[-1].split("[")[0] + + +def _val(v: dict) -> dict: + """An offender operand as text plus its raw IEEE bits.""" + return {"text": _fmt_val(v["value"]), "bits": v["bits"]} + + +def _merge_offenders(recs: list) -> list: + """Worst offenders across every test in the group, highest |error| first.""" + flat = [(r, o) for r in recs for o in r["worst"]] + flat.sort(key=lambda t: _mag(t[1]["error"]), reverse=True) + return [ + { + "test": _test_name(r["label"]), + "x": _val(o["x"]), + "actual": _val(o["actual"]), + "expected": _val(o["expected"]), + "residual": f"{o['residual']:+.4f}", + "error": _fmt_ulp(o["error"]), + "over": math.isinf(m := _mag(o["error"])) or m > r["maxulp"], + } + for r, o in flat[:_OFFENDER_LIMIT] + ] + + +def _row(target, operation, dtype, accuracy, recs, have_fma) -> dict: + """One aggregate table row folding every test that hit this combination.""" + n = sum(r["n"] for r in recs) + n_over = sum(r["n_over"] for r in recs) + mean = sum(_pnum(r["mean_error"]) * r["n"] for r in recs) / n if n else 0.0 + lo, hi = min(r["maxulp"] for r in recs), max(r["maxulp"] for r in recs) + tests = sorted( + recs, key=lambda r: _ratio(r["max_error"], r["maxulp"]), reverse=True + ) + worst = tests[0] + ratio = _ratio(worst["max_error"], worst["maxulp"]) + return { + "operation": operation, + "target": target, + "fma": "FMA" if have_fma else "non-FMA", + "dtype": dtype, + "accuracy": accuracy, + "status": "pass" if all(r["passed"] for r in recs) else "fail", + "tests": len(recs), + # en dash is a deliberate range separator in the rendered report + "bound": _fmt_ulp(lo) if lo == hi else f"{_fmt_ulp(lo)}โ€“{_fmt_ulp(hi)}", # noqa: RUF001 + "boundSort": _finite(hi), + "maxErr": _fmt_ulp(worst["max_error"]), + "maxErrNum": _finite(_pnum(worst["max_error"])), + "barPct": min(1.0, ratio) * 100 if math.isfinite(ratio) else 100.0, + "barOver": ratio > 1, + "mean": _fmt_ulp(mean), + "meanNum": mean, + "over": f"{n_over}/{n}", + "nOver": n_over, + "n": n, + "search": " ".join(r["label"] for r in recs).lower(), + "perTest": [ + { + "label": r["label"], + "bound": _fmt_ulp(r["maxulp"]), + "maxErr": _fmt_ulp(r["max_error"]), + "mean": _fmt_ulp(r["mean_error"]), + "over": f"{r['n_over']}/{r['n']}", + "passed": r["passed"], + } + for r in tests + ], + "offenders": _merge_offenders(recs), + } + + +def _serialize_records(records) -> dict: + # Rows are grouped and ordered by these four dims. + dims = ("target", "operation", "dtype", "accuracy") + groups: dict = {} + for r in records: + groups.setdefault(tuple(r[d] for d in dims), []).append(r) + # first-appearance order, which for each dim matches record order + order = {d: list(dict.fromkeys(k[i] for k in groups)) for i, d in enumerate(dims)} + + # TARGETS (not sr.targets) so the numpy reference mock resolves too. + fma_of = {t: TARGETS[t].HAVE_FMA for t in order["target"]} + rows = [ + _row(target, op, dtype, acc, recs, fma_of[target]) + for (target, op, dtype, acc), recs in groups.items() + ] + + def present(field, opts): + seen = {r[field] for r in rows} + return [o for o in opts if o in seen] + + # Facet lists in display order; the viewer just renders and toggles them. + filter_by = { + "operation": order["operation"], + "target": order["target"], + "fma": present("fma", ["FMA", "non-FMA"]), + "dtype": order["dtype"], + "accuracy": order["accuracy"], + "status": present("status", ["pass", "fail"]), + } + return {"filter_by": filter_by, "records": rows} diff --git a/python/_conftest/stash.py b/python/_conftest/stash.py new file mode 100644 index 0000000..91c8b31 --- /dev/null +++ b/python/_conftest/stash.py @@ -0,0 +1,11 @@ +"""Session stash keys; one module so writer and readers share key identity.""" + +import pytest + +from .book import WorstBook + +# One record per ULP comparison (pass or fail), drained into the JSON report. +REPORT_KEY = pytest.StashKey[list]() + +# The frozen worst-case book (tests/worst/*.csv), folded during the run. +BOOK_KEY = pytest.StashKey[WorstBook]() diff --git a/python/_conftest/targets.py b/python/_conftest/targets.py new file mode 100644 index 0000000..33621c1 --- /dev/null +++ b/python/_conftest/targets.py @@ -0,0 +1,45 @@ +"""The selectable target and accuracy registries.""" + +from collections.abc import Callable + +import numpy as np +import numpy_sr as sr + + +class _TargetMock(type): + """A module as a reference target: swallows `accuracy`, ignores errstate.""" + + IS_REFERENCE = True + + def __new__(mcls, name: str, module, flags: dict) -> type: + return super().__new__(mcls, name, (), {"MODULE": module, **flags}) + + def __getattr__(cls, op: str) -> Callable: + fn = getattr(cls.MODULE, op) + + @np.errstate(all="ignore") + def call(*args, accuracy=None): + return fn(*args) + + return call + + +# `sr.targets` already holds the compiled references, so only numpy is manual. +# Excluded: what this CPU cannot run, and the MPFR oracle (returns tuples). +TARGETS: dict = { + name: target + for name, target in sr.targets.items() + if target.IS_SUPPORTED and not target.IS_ORACLE +} | { + "numpy": _TargetMock("numpy", np, {"HAVE_FLOAT64": True, "HAVE_FMA": True}), +} + +# References stay selectable via --target but out of the default run. +DEFAULT_TARGETS: list = [ + name + for name, target in TARGETS.items() + if not getattr(target, "IS_REFERENCE", False) +] + +ACCURACIES: dict = {a.name.lower(): a for a in sr.Accuracy} +DEFAULT_ACCURACIES: list = list(sr.Accuracy) diff --git a/python/_conftest/xdist.py b/python/_conftest/xdist.py new file mode 100644 index 0000000..cd7b0c0 --- /dev/null +++ b/python/_conftest/xdist.py @@ -0,0 +1,52 @@ +"""xdist glue: one seed per run, worker results merged back on the controller.""" + +from typing import cast + +import pytest + +from .stash import BOOK_KEY, REPORT_KEY + + +def worker_output(config: pytest.Config) -> dict | None: + """The channel back to the controller; `None` unless this is a worker.""" + return getattr(config, "workeroutput", None) + + +def resolve_seed(config: pytest.Config) -> int: + """The run's root seed: the controller's, else this process's own default.""" + wi = getattr(config, "workerinput", None) + if wi is not None and "npsr_seed" in wi: + return int(wi["npsr_seed"]) + return cast(int, config.getoption("--seed")) + + +def warn_unmerged(config: pytest.Config) -> None: + """A worker whose controller lacks this plugin: its results are dropped.""" + wi = getattr(config, "workerinput", None) + if wi is not None and "npsr_seed" not in wi: + config.issue_config_time_warning( + pytest.PytestConfigWarning( + "xdist controller is not running the numpy_sr plugin: workers " + "seed independently, and neither the ULP report nor the " + "worst-case book is kept" + ), + stacklevel=2, + ) + + +# optionalhook: the specs live in xdist, which need not be installed. +@pytest.hookimpl(optionalhook=True) +def pytest_configure_node(node) -> None: + """Controller -> worker: the one seed every random stream derives from.""" + node.workerinput["npsr_seed"] = resolve_seed(node.config) + + +@pytest.hookimpl(optionalhook=True) +def pytest_testnodedown(node, error) -> None: + """Absorb a finished worker's report records and book deltas.""" + out = getattr(node, "workeroutput", None) or {} # a crashed node sends none + config = node.config + config.stash[REPORT_KEY].extend(out.get("npsr_records", [])) + config.stash[BOOK_KEY].absorb( + out.get("npsr_book", {}), out.get("npsr_restamped", {}) + ) diff --git a/python/_numpy_sr.pyi b/python/_numpy_sr.pyi new file mode 100644 index 0000000..4eda64e --- /dev/null +++ b/python/_numpy_sr.pyi @@ -0,0 +1,51 @@ +"""Type stub for the compiled ``_numpy_sr`` extension.""" + +from __future__ import annotations + +import enum +from typing import Any, Protocol, overload + +import numpy as np +import numpy.typing as npt + +# Arguments are bound via SafeArg (no convert): dtype must match exactly, so ArrayLike +# would overpromise. Any layout is accepted; strided input is copied. +_F32 = npt.NDArray[np.float32] +_F64 = npt.NDArray[np.float64] +_Float = _F32 | _F64 +_RefResidual = tuple[npt.NDArray[Any], npt.NDArray[Any]] + +class Accuracy(enum.Enum): + """Accuracy profile selecting a kernel's error bound.""" + + High = 0 + Low = 1 + +class _Backend(Protocol): + """Flags every target submodule carries, kernel or reference alike. + + ``IS_SUPPORTED`` is False when the build included the target but this CPU + cannot run it -- such a submodule is bound *without* its ops. + ``HWY_TARGET`` is Highway's target bit (lower = better ISA, 0 when the + backend is not a Highway target).""" + + HAVE_FMA: bool + HAVE_FLOAT64: bool + IS_REFERENCE: bool + IS_ORACLE: bool + IS_SUPPORTED: bool + HWY_TARGET: int + +class _Target(_Backend, Protocol): + """One dispatched ISA target (e.g. ``numpy_sr.AVX2``).""" + def sin(self, x: _Float, accuracy: Accuracy = ...) -> npt.NDArray[Any]: ... + def cos(self, x: _Float, accuracy: Accuracy = ...) -> npt.NDArray[Any]: ... + +targets: dict[str, _Target] + +# Dispatched to the best supported target at call time. +@overload +def ulp_error(computed: _F32, ref: _F32, residual: _F64 | None = ...) -> _F64: ... +@overload +def ulp_error(computed: _F64, ref: _F64, residual: _F64 | None = ...) -> _F64: ... +def ulp_error32(computed: _F32, oracle: _F64) -> _F64: ... diff --git a/python/_numpy_sr/cos.cpp b/python/_numpy_sr/cos.cpp new file mode 100644 index 0000000..f90e602 --- /dev/null +++ b/python/_numpy_sr/cos.cpp @@ -0,0 +1,4 @@ +#define NPSR_OP Cos +#undef HWY_TARGET_INCLUDE +#define HWY_TARGET_INCLUDE "python/_numpy_sr/cos.cpp" +#include "python/_numpy_sr/op-inl.h" diff --git a/python/_numpy_sr/meson.build b/python/_numpy_sr/meson.build new file mode 100644 index 0000000..d0a5fd9 --- /dev/null +++ b/python/_numpy_sr/meson.build @@ -0,0 +1,70 @@ +hwy_proj = subproject( + 'highway', + default_options: [ + 'contrib=disabled', + 'examples=disabled', + # 'tests' alone still pulls googletest; 'test_standalone' cuts it too. + 'tests=disabled', + 'test_standalone=true', + ], +) +hwy_dep = hwy_proj.get_variable('hwy_dep') + +cpp = meson.get_compiler('cpp') +mpfr_dep = dependency('mpfr', required: false) +mpfr_rpath = '' +if mpfr_dep.found() + mpfr_rpath = mpfr_dep.get_variable(pkgconfig: 'libdir', default_value: '') +else + mpfr_dep = cpp.find_library('mpfr', has_headers: ['mpfr.h']) +endif + + +pybind11_dep = dependency('pybind11') + +# Intel SVML: an AVX-512 comparison reference. numpy/SVML ships x86-64 Linux +# assembly only; elsewhere the target binds no ops (NPSR_HAVE_SVML, target.h). +svml_enabled = ( + host_machine.system() == 'linux' and host_machine.cpu_family() == 'x86_64' +) +svml_link = [] +svml_cpp_args = ['-DNPSR_HAVE_SVML=0'] +if svml_enabled + svml_dep = subproject('svml').get_variable('svml_dep') + # Isolated in its own -mavx512f library so the module still loads at baseline + # ISA; these kernels only run when explicitly called (runtime-guarded). + svml_link = [static_library( + 'npsr_svml', + 'svml.cpp', + include_directories: npsr_inc, + # NPSR_HAVE_SVML must match module.cpp's view of TargetSVML::kIsEnabled. + cpp_args: ['-mavx512f', '-DNPSR_HAVE_SVML=1'], + gnu_symbol_visibility: 'hidden', + # hwy_dep: target.h reads HWY_AVX3 for the ISA gate (headers only). + dependencies: [pybind11_dep, hwy_dep, svml_dep], + pic: true, + )] + svml_cpp_args = ['-DNPSR_HAVE_SVML=1'] +endif + +# One foreach_target TU per op: SIMD variants build in parallel and touching +# one op recompiles only that op. +op_sources = [ + 'sin.cpp', + 'cos.cpp', +] + +py.extension_module( + '_numpy_sr', + ['module.cpp', 'mpfr.cpp', 'ulp.cpp'] + op_sources, + include_directories: npsr_inc, + cpp_args: svml_cpp_args, + dependencies: [pybind11_dep, hwy_dep, mpfr_dep], + link_with: svml_link, + install_rpath: mpfr_rpath, + # pybind11's namespace is hidden; matching it keeps our own glue types from + # out-scoping their members (PYBIND11_MODULE still exports PyInit_). + gnu_symbol_visibility: 'hidden', + subdir: 'numpy_sr', + install: true, +) diff --git a/python/_numpy_sr/module.cpp b/python/_numpy_sr/module.cpp new file mode 100644 index 0000000..3ddd1a4 --- /dev/null +++ b/python/_numpy_sr/module.cpp @@ -0,0 +1,113 @@ +#undef HWY_TARGET_INCLUDE +#define HWY_TARGET_INCLUDE "python/_numpy_sr/module.cpp" +#include "hwy/foreach_target.h" // IWYU pragma: keep + +#include + +#include "hwy/highway.h" +#include "npsr/precise.h" + +#include "python/_numpy_sr/target-inl.h" + +#if HWY_ONCE +#include + +#include +#include + +namespace npsr::py { +namespace { + +// A native enum.Enum, not pb::enum_/IntEnum: no int can stand in. Must be +// finalize()d before any def() defaulting to it. +template +void BindAccuracy(Module& m, std::tuple) { + pb::native_enum e(m, "Accuracy", "enum.Enum", + "Accuracy profile selecting a kernel's error " + "bound."); + (e.value(ACCS::kPyName, ACCS::kID), ...); + e.finalize(); +} + +// Flag primitives; scoping lives in numpy_sr.fpenv. Covers INEXACT, which the +// kernels leave unmanaged, so ERRORS masks it back out. +struct FPEnv { + static void Clear() noexcept { std::feclearexcept(FPExceptions::kAll); } + static int Test() noexcept { return std::fetestexcept(FPExceptions::kAll); } +}; + +// Bound without a constructor: a namespace of flag bits, not an object. +void BindFPEnv(Module& m) { + pb::class_ c(m, "FPEnv", + "Thread FP-exception flags: the FE_* bits plus the two " + "primitives acting on them."); + c.def_static("clear", &FPEnv::Clear, + "Clear all floating-point exception flags."); + c.def_static("test", &FPEnv::Test, + "Currently raised FP exception flags (FE_* bits)."); + c.attr("INVALID") = FPExceptions::kInvalid; + c.attr("DIVBYZERO") = FPExceptions::kDivByZero; + c.attr("OVERFLOW") = FPExceptions::kOverflow; + c.attr("UNDERFLOW") = FPExceptions::kUnderflow; + c.attr("INEXACT") = FPExceptions::kInexact; + c.attr("ERRORS") = FPExceptions::kAll; +} + +template +void BindOperation(Module& m) { + Binder binder(m, OP::kPyName); + T::template Bind(binder); +} + +template +void BindOperations(Module& m, std::tuple) { + (BindOperation(m), ...); +} + +template +void BindTarget(Module& m, pb::dict& d) { + Module sub = m.def_submodule(T::kPyName, ""); + bool supported = T::kIsEnabled; + if constexpr (T::kIsEnabled) { + if constexpr (T::kHighwayID != 0) { + supported = (::hwy::SupportedTargets() & T::kHighwayID) != 0; + } + if (supported) { + BindOperations(sub, typename T::Operations{}); + } + } + sub.attr("__name__") = T::kPyName; + sub.attr("HAVE_FMA") = T::kHaveFMA; + sub.attr("HAVE_FLOAT64") = T::kHaveF64; + sub.attr("IS_REFERENCE") = T::kIsReference; + sub.attr("IS_ORACLE") = T::kIsOracle; + sub.attr("IS_SUPPORTED") = supported; + // Highway's target bit: lower = better ISA, 0 for non-Highway. The dict is + // emit-order (alphabetic), so Python ranks targets with this. + sub.attr("HWY_TARGET") = T::kHighwayID; + d[T::kPyName] = sub; +} + +} // namespace +} // namespace npsr::py + +PYBIND11_MODULE(_numpy_sr, m) { + using namespace npsr::py; + m.doc() = "NumPy SIMD routines: per-target vectorized math kernels."; + + BindAccuracy(m, Accuracy::All{}); + BindFPEnv(m); + BindULP(m); + + pb::dict targets_dict; + BindTarget(m, targets_dict); + BindTarget(m, targets_dict); + +#define NPSR_BIND_HWY_TARGET(TARGET, NAMESPACE) \ + BindTarget(m, targets_dict); + HWY_VISIT_TARGETS(NPSR_BIND_HWY_TARGET) +#undef NPSR_BIND_HWY_TARGET + + m.attr("targets") = targets_dict; +} +#endif // HWY_ONCE diff --git a/python/_numpy_sr/mpfr.cpp b/python/_numpy_sr/mpfr.cpp new file mode 100644 index 0000000..ffad0f7 --- /dev/null +++ b/python/_numpy_sr/mpfr.cpp @@ -0,0 +1,97 @@ +// Correctly-rounded scalar oracle: (x) -> (ref, residual), where residual +// is the exact value's signed offset from ref in ULP of ref. +#include + +#include +#include +#include +#include +#include + +#include "python/_numpy_sr/target.h" + +namespace npsr::py { +namespace { + +using MpfrFn = int (*)(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t); + +template +void MpfrRefResidual(const T* src, T* ref_out, double* res_out, size_t n) { + constexpr bool kIsF32 = std::is_same_v; + constexpr mpfr_prec_t kPrec = kIsF32 ? 24 : 53; + const mpfr_exp_t emin = mpfr_get_emin(), emax = mpfr_get_emax(); + mpfr_set_emin(kIsF32 ? -148 : -1073); + mpfr_set_emax(kIsF32 ? 128 : 1024); + mpfr_t x, hi, ref, diff; + mpfr_init2(x, kPrec); + mpfr_init2(hi, kPrec + 64); // guard bits: hi is the exact value + mpfr_init2(ref, kPrec); + mpfr_init2(diff, kPrec + 64); + for (size_t i = 0; i < n; ++i) { + if constexpr (kIsF32) { + mpfr_set_flt(x, src[i], MPFR_RNDN); + } else { + mpfr_set_d(x, src[i], MPFR_RNDN); + } + Fn(hi, x, MPFR_RNDN); + const int t = mpfr_set(ref, hi, MPFR_RNDN); + mpfr_subnormalize(ref, t, MPFR_RNDN); + const T r = kIsF32 ? mpfr_get_flt(ref, MPFR_RNDN) + : static_cast(mpfr_get_d(ref, MPFR_RNDN)); + ref_out[i] = r; + if (!std::isfinite(r)) { + res_out[i] = 0.0; + continue; + } + const T a = std::fabs(r); + const double ulp = + a == T(0) + ? std::ldexp(1.0, kIsF32 ? -149 : -1074) + : static_cast( + std::nextafter(a, std::numeric_limits::infinity()) - a); + mpfr_sub(diff, ref, hi, MPFR_RNDN); + res_out[i] = mpfr_get_d(diff, MPFR_RNDN) / ulp; + } + mpfr_clear(x); + mpfr_clear(hi); + mpfr_clear(ref); + mpfr_clear(diff); + mpfr_set_emin(emin); + mpfr_set_emax(emax); +} + +template +pb::tuple MpfrRefUlp(Array x) { + const CArray xc = Contiguous(x); + const auto shape = xc.request().shape; + Array ref(shape); + Array res(shape); + const size_t n = static_cast(xc.size()); + { + pb::gil_scoped_release nogil; + MpfrRefResidual(xc.data(), ref.mutable_data(), res.mutable_data(), + n); + } + return pb::make_tuple(std::move(ref), std::move(res)); +} + +// Accuracy is taken for signature parity and ignored: always correctly rounded. +template +void BindOp(Binder& b) { + b([](Array x, AccuracyID) { return MpfrRefUlp(x); }); + b([](Array x, AccuracyID) { return MpfrRefUlp(x); }); +} + +} // namespace + +template <> +void TargetMPFR::Bind(Binder& b) { + BindOp(b); +} + +template <> +void TargetMPFR::Bind(Binder& b) { + BindOp(b); +} + +} // namespace npsr::py diff --git a/python/_numpy_sr/op-inl.h b/python/_numpy_sr/op-inl.h new file mode 100644 index 0000000..cbf3773 --- /dev/null +++ b/python/_numpy_sr/op-inl.h @@ -0,0 +1,18 @@ +// Body of a per-op foreach_target TU; no include guard, re-included once per +// target. The op .cpp defines NPSR_OP + HWY_TARGET_INCLUDE, then includes this. +#include "hwy/foreach_target.h" // IWYU pragma: keep +#include "hwy/highway.h" +#include "npsr/npsr.h" +#include "python/_numpy_sr/target-inl.h" + +HWY_BEFORE_NAMESPACE(); +namespace npsr::py::HWY_NAMESPACE { +namespace sr = ::npsr::HWY_NAMESPACE; + +template <> +void TargetHighway::Bind(Binder& b) { + BindKernel(b, [](auto& prec, auto v) { return sr::NPSR_OP(prec, v); }); +} + +} // namespace npsr::py::HWY_NAMESPACE +HWY_AFTER_NAMESPACE(); diff --git a/python/_numpy_sr/sin.cpp b/python/_numpy_sr/sin.cpp new file mode 100644 index 0000000..5c58e2d --- /dev/null +++ b/python/_numpy_sr/sin.cpp @@ -0,0 +1,4 @@ +#define NPSR_OP Sin +#undef HWY_TARGET_INCLUDE +#define HWY_TARGET_INCLUDE "python/_numpy_sr/sin.cpp" +#include "python/_numpy_sr/op-inl.h" diff --git a/python/_numpy_sr/svml.cpp b/python/_numpy_sr/svml.cpp new file mode 100644 index 0000000..7e52d2a --- /dev/null +++ b/python/_numpy_sr/svml.cpp @@ -0,0 +1,91 @@ +// Intel SVML as the `svml_avx512` reference target: a baseline, not an oracle. +// Built -mavx512f in its own library, behind TargetSVML's ISA gate. +#include + +#include + +#include "python/_numpy_sr/target.h" + +namespace npsr::py { +namespace { + +// Bare name = low accuracy, _ha = high. The vector ABI passes and returns in +// zmm0, so a plain C declaration and call is correct (as numpy does it). +extern "C" { +__m512d __svml_sin8(__m512d); +__m512d __svml_sin8_ha(__m512d); +__m512d __svml_cos8(__m512d); +__m512d __svml_cos8_ha(__m512d); +__m512 __svml_sinf16(__m512); +__m512 __svml_sinf16_ha(__m512); +__m512 __svml_cosf16(__m512); +__m512 __svml_cosf16_ha(__m512); +} + +using F64Fn = __m512d (*)(__m512d); +using F32Fn = __m512 (*)(__m512); + +// Full vectors, then a masked tail pass: padded lanes compute but never store. +void Run(F64Fn fn, const double* src, double* dst, size_t n) { + size_t i = 0; + for (; i + 8 <= n; i += 8) { + _mm512_storeu_pd(dst + i, fn(_mm512_loadu_pd(src + i))); + } + if (const size_t rem = n - i) { + const __mmask8 m = static_cast<__mmask8>((1u << rem) - 1u); + _mm512_mask_storeu_pd(dst + i, m, fn(_mm512_maskz_loadu_pd(m, src + i))); + } +} + +void Run(F32Fn fn, const float* src, float* dst, size_t n) { + size_t i = 0; + for (; i + 16 <= n; i += 16) { + _mm512_storeu_ps(dst + i, fn(_mm512_loadu_ps(src + i))); + } + if (const size_t rem = n - i) { + const __mmask16 m = static_cast<__mmask16>((1u << rem) - 1u); + _mm512_mask_storeu_ps(dst + i, m, fn(_mm512_maskz_loadu_ps(m, src + i))); + } +} + +// TargetSVML::Accuracies is {High, Low}; anything else is a Python-level error. +template +Array Apply(Fn ha, Fn la, Array x, AccuracyID acc_id) { + if (acc_id != Accuracy::High::kID && acc_id != Accuracy::Low::kID) { + throw pb::value_error("unsupported accuracy"); + } + const Fn fn = acc_id == Accuracy::High::kID ? ha : la; + const CArray xc = Contiguous(x); + Array out(xc.request().shape); + const T* src = xc.data(); + T* dst = out.mutable_data(); + const size_t n = static_cast(xc.size()); + { + pb::gil_scoped_release nogil; + Run(fn, src, dst, n); + } + return out; +} + +void BindOp(Binder& b, F32Fn s_ha, F32Fn s_la, F64Fn d_ha, F64Fn d_la) { + b([s_ha, s_la](Array x, AccuracyID acc_id) { + return Apply(s_ha, s_la, x, acc_id); + }); + b([d_ha, d_la](Array x, AccuracyID acc_id) { + return Apply(d_ha, d_la, x, acc_id); + }); +} + +} // namespace + +template <> +void TargetSVML::Bind(Binder& b) { + BindOp(b, __svml_sinf16_ha, __svml_sinf16, __svml_sin8_ha, __svml_sin8); +} + +template <> +void TargetSVML::Bind(Binder& b) { + BindOp(b, __svml_cosf16_ha, __svml_cosf16, __svml_cos8_ha, __svml_cos8); +} + +} // namespace npsr::py diff --git a/python/_numpy_sr/target-inl.h b/python/_numpy_sr/target-inl.h new file mode 100644 index 0000000..9e5a4ac --- /dev/null +++ b/python/_numpy_sr/target-inl.h @@ -0,0 +1,85 @@ +#if defined(NPSR_PY_TARGET_INL_H_) == defined(HWY_TARGET_TOGGLE) // NOLINT +#ifdef NPSR_PY_TARGET_INL_H_ +#undef NPSR_PY_TARGET_INL_H_ +#else +#define NPSR_PY_TARGET_INL_H_ +#endif + +#include + +#include "hwy/highway.h" +#include "python/_numpy_sr/target.h" + +HWY_BEFORE_NAMESPACE(); +namespace npsr::py::HWY_NAMESPACE { +namespace hn = ::hwy::HWY_NAMESPACE; + +struct TargetHighway : public Target { + using Operations = Operation::All; + using Accuracies = Accuracy::All; + // +2 skips the "N_" of HWY_NAMESPACE: TargetName(), but at compile time. + static constexpr const char* kPyName = HWY_STR(HWY_NAMESPACE) + 2; + static constexpr bool kHaveFMA = HWY_NATIVE_FMA != 0; + static constexpr bool kHaveF64 = HWY_HAVE_FLOAT64 != 0; + static constexpr bool kIsHighway = true; + static constexpr int64_t kHighwayID = HWY_TARGET; + + template + static void Bind(Binder& b); +}; + +// Masked LoadN/StoreN so full vectors and the tail share one op call site. +template +void Forward(Op op, const T* HWY_RESTRICT src, T* HWY_RESTRICT dst, + size_t len) { + Prec prec; + const hn::ScalableTag d; + const size_t N = hn::Lanes(d); + for (size_t i = 0; i < len; i += N) { + const size_t rem = len - i; + hn::StoreN(op(prec, hn::LoadN(d, src + i, rem)), d, dst + i, rem); + } +} + +template +Array Apply(std::tuple, Op op, Array x, AccuracyID acc_id) { + const CArray xc = Contiguous(x); + Array out(xc.request().shape); + const T* src = xc.data(); + T* dst = out.mutable_data(); + const size_t n = static_cast(xc.size()); + bool found = false; + { + pb::gil_scoped_release nogil; // released before the array_t return + const auto run = [&](auto acc) { + using Acc = decltype(acc); + if (Acc::kID == acc_id) { + found = true; + Forward(op, src, dst, n); + } + }; + (run(ACCS{}), ...); + } + if (!found) throw pb::value_error("unsupported accuracy"); + return out; +} + +// The f32/f64 overloads of one op; dtype is exact, so their order is free. +template +void BindKernel(Binder& b, Op op) { + const auto def = [&](auto zero) { + using T = decltype(zero); + b([op](Array x, AccuracyID acc_id) { + return Apply(TargetHighway::Accuracies{}, op, x, acc_id); + }); + }; + def(float{}); +#if HWY_HAVE_FLOAT64 + def(double{}); +#endif +} + +} // namespace npsr::py::HWY_NAMESPACE +HWY_AFTER_NAMESPACE(); + +#endif // NPSR_PY_TARGET_INL_H_ diff --git a/python/_numpy_sr/target.h b/python/_numpy_sr/target.h new file mode 100644 index 0000000..bdbf5b5 --- /dev/null +++ b/python/_numpy_sr/target.h @@ -0,0 +1,133 @@ +#ifndef NPSR_PY_TARGET_H_ +#define NPSR_PY_TARGET_H_ + +#include +#include + +#include +#include +#include + +#include "hwy/detect_targets.h" +#include "npsr/precise.h" + +// Set by meson; 0 keeps a TU that misses the flag from silently failing to +// resolve TargetSVML::kIsEnabled. +#ifndef NPSR_HAVE_SVML +#define NPSR_HAVE_SVML 0 +#endif + +namespace npsr::py { + +namespace pb = ::pybind11; + +struct Operation { + struct Sin { + static constexpr int kID = 0; + static constexpr const char* kPyName = "sin"; + }; + struct Cos { + static constexpr int kID = 1; + static constexpr const char* kPyName = "cos"; + }; + using All = std::tuple; +}; + +enum class AccuracyID { kHigh = 0, kLow = 1}; + +struct Accuracy { + struct High { + static constexpr AccuracyID kID = AccuracyID::kHigh; + static constexpr const char* kPyName = "High"; + using Prec = decltype(Precise{}); + }; + struct Low { + static constexpr AccuracyID kID = AccuracyID::kLow; + static constexpr const char* kPyName = "Low"; + using Prec = decltype(Precise{kLowAccuracy}); + }; + using All = std::tuple; +}; + +template +using Array = pb::array_t; + +template +using CArray = pb::array_t; + +// Same dtype, C-contiguous: copies strided or broadcast input, never casts. +template +CArray Contiguous(const Array& x) { + CArray c = CArray::ensure(x); + if (!c) throw pb::value_error("could not make the input C-contiguous"); + return c; +} + +// Bind every array argument through this: dropping forcecast blocks only +// *unsafe* casts, pybind's convert pass still promotes int/f32/list/scalar into +// the wider overload. Layout is left to Contiguous() instead. +inline pb::arg SafeArg(const char* name) { return pb::arg(name).noconvert(); } + +using Module = pb::module_; + +// The one def() site, so every op shares `(x, accuracy=Accuracy.High)`. +// Accuracy must be finalize()d before the default value is read. +class Binder { + public: + Binder(Module& m, const char* op_name) : mod_(m), op_name_(op_name) {} + + template + void operator()(F&& f) { + mod_.def(op_name_, std::forward(f), SafeArg("x"), + pb::arg("accuracy") = Accuracy::High::kID); + } + + private: + Module& mod_; + const char* op_name_; +}; + +struct Target { + static constexpr const char* kPyName = "?"; + static constexpr bool kHaveFMA = true; + static constexpr bool kHaveF64 = true; + static constexpr bool kIsReference = false; + static constexpr bool kIsOracle = false; + static constexpr bool kIsEnabled = true; + static constexpr bool kIsHighway = false; + // Non-zero = the HWY_* ISA bit this target needs; checked against + // hwy::SupportedTargets() before binding. + static constexpr int64_t kHighwayID = 0; +}; + +struct TargetSVML : public Target { + using Operations = Operation::All; + using Accuracies = std::tuple; + + static constexpr const char* kPyName = "svml_avx512"; + static constexpr bool kIsReference = true; + // x86-64 Linux only (see meson.build); AVX-512 kernels, hence the ISA gate. + static constexpr bool kIsEnabled = NPSR_HAVE_SVML != 0; + static constexpr int64_t kHighwayID = HWY_AVX3; + + template + static void Bind(Binder& b); +}; + +struct TargetMPFR : public Target { + using Operations = Operation::All; + using Accuracies = std::tuple; + + static constexpr const char* kPyName = "mpfr"; + static constexpr bool kIsOracle = true; + + template + static void Bind(Binder& b); +}; + +// Defined in ulp.cpp, dispatched to the best supported target at call time. +void BindULP(Module& m); + +} // namespace npsr::py + +#endif // NPSR_PY_TARGET_H_ diff --git a/python/_numpy_sr/ulp.cpp b/python/_numpy_sr/ulp.cpp new file mode 100644 index 0000000..0ce28fe --- /dev/null +++ b/python/_numpy_sr/ulp.cpp @@ -0,0 +1,243 @@ +// Fractional-ULP comparators: utilities, so they bind at module level and +// dispatch to the best supported target at call time. +#undef HWY_TARGET_INCLUDE +#define HWY_TARGET_INCLUDE "python/_numpy_sr/ulp.cpp" +#include "hwy/foreach_target.h" // IWYU pragma: keep + +#include // std::optional> residual + +#include +#include +#include +#include +#include + +#include "hwy/highway.h" +#include "python/_numpy_sr/target.h" + +HWY_BEFORE_NAMESPACE(); +namespace npsr::py::HWY_NAMESPACE { +namespace hn = ::hwy::HWY_NAMESPACE; + +template +V UlpError32(V c, V o) { + using namespace hn; + using D = DFromV; + using DU = RebindToUnsigned; + using VU = VFromD; + const D d; + const DU du; + const VU bias = Set(du, uint64_t{2069}); + const VU expmask = Set(du, uint64_t{0x7FF}); + const VU e = And(ShiftRight<52>(BitCast(du, o)), expmask); + const V inv = BitCast(d, ShiftLeft<52>(Sub(bias, e))); + return Mul(Abs(Sub(c, o)), inv); +} + +// f32 lanes widen to the double lane width; f64 already match. +template +hn::VFromD WidenTo(D d, V v) { + if constexpr (hwy::IsSame, hn::TFromD>()) { + return v; + } else { + return hn::PromoteTo(d, v); + } +} + +template +VD UlpError(V ct, V rt, VD res) { + using namespace hn; + using T = TFromV; + using TI = hwy::MakeSigned; + using DD = DFromV; + using DI = RebindToSigned
; + using DTI = Rebind; + using VTI = VFromD; + using MD = MFromD
; + const DD dd; + const DI di; + const DTI dti; + const VTI kmin = Set(dti, hwy::LimitsMin()); + const VD inf = Set(dd, std::numeric_limits::infinity()); + const VTI cbits = BitCast(dti, ct); + const VTI rbits = BitCast(dti, rt); + const VTI ckey = IfNegativeThenElse(cbits, Sub(kmin, cbits), cbits); + const VTI rkey = IfNegativeThenElse(rbits, Sub(kmin, rbits), rbits); + // int64 wraps like the numpy key subtraction; non-finites are fixed below. + const VFromD diff = Sub(WidenTo(di, ckey), WidenTo(di, rkey)); + VD err = Add(ConvertTo(dd, diff), res); + // Classify in f64: promotion preserves NaN-ness, Inf-ness and sign. + const VD c = WidenTo(dd, ct); + const VD r = WidenTo(dd, rt); + const MD c_nan = IsNaN(c), r_nan = IsNaN(r); + const MD c_inf = IsInf(c), r_inf = IsInf(r); + const MD sign_ne = RebindMask( + dd, Xor(Lt(BitCast(di, c), Zero(di)), Lt(BitCast(di, r), Zero(di)))); + const MD both_inf = And(c_inf, r_inf); + const MD match = Or(And(c_nan, r_nan), AndNot(sign_ne, both_inf)); + const MD mismatch = Or(Or(Xor(c_nan, r_nan), Xor(c_inf, r_inf)), + And(both_inf, sign_ne)); + err = IfThenZeroElse(match, err); + return IfThenElse(mismatch, inf, err); +} + +#if HWY_HAVE_FLOAT64 +// |computed - oracle| in f32-ULP units; 1/ulp32(oracle) = 2^(2069-E) built +// straight from the biased f64 exponent E, no frexp/divide. +void KernelUlpError32(const float* HWY_RESTRICT comp, + const double* HWY_RESTRICT oracle, + double* HWY_RESTRICT out, size_t n) { + const hn::ScalableTag dd; + const hn::Rebind df; + const size_t N = hn::Lanes(dd); + for (size_t i = 0; i < n; i += N) { + const size_t rem = n - i; + const auto c = hn::PromoteTo(dd, hn::LoadN(df, comp + i, rem)); + const auto o = hn::LoadN(dd, oracle + i, rem); + hn::StoreN(UlpError32(c, o), dd, out + i, rem); + } +} +// Signed fractional-ULP error: IEEE total-order key difference plus the +// nullable MPFR residual. Matching non-finites score 0, mismatches +inf. +template +void LoopUlpError(const T* HWY_RESTRICT comp, const T* HWY_RESTRICT ref, + const double* HWY_RESTRICT res, double* HWY_RESTRICT out, + size_t n) { + const hn::ScalableTag dd; + const hn::Rebind dt; + const size_t N = hn::Lanes(dd); + for (size_t i = 0; i < n; i += N) { + const size_t rem = n - i; + const auto ct = hn::LoadN(dt, comp + i, rem); + const auto rt = hn::LoadN(dt, ref + i, rem); + const auto rs = res ? hn::LoadN(dd, res + i, rem) : hn::Zero(dd); + hn::StoreN(UlpError(ct, rt, rs), dd, out + i, rem); + } +} +#else +// Bit-identical scalar fallbacks for targets without native float64 vectors. +void KernelUlpError32(const float* HWY_RESTRICT comp, + const double* HWY_RESTRICT oracle, + double* HWY_RESTRICT out, size_t n) { + for (size_t i = 0; i < n; ++i) { + const double o = oracle[i]; + const uint64_t e = (hwy::BitCastScalar(o) >> 52) & 0x7FF; + const double inv = hwy::BitCastScalar((uint64_t{2069} - e) << 52); + out[i] = std::fabs(static_cast(comp[i]) - o) * inv; + } +} + +template +void LoopUlpError(const T* HWY_RESTRICT comp, const T* HWY_RESTRICT ref, + const double* HWY_RESTRICT res, double* HWY_RESTRICT out, + size_t n) { + using TI = hwy::MakeSigned; + const auto key = [](T v) { + const TI bits = hwy::BitCastScalar(v); + return bits < 0 ? static_cast(hwy::LimitsMin() - bits) : bits; + }; + for (size_t i = 0; i < n; ++i) { + const T c = comp[i], r = ref[i]; + // unsigned subtraction wraps like the int64 vector/numpy key math + const int64_t diff = static_cast( + static_cast(key(c)) - static_cast(key(r))); + double err = static_cast(diff); + if (res) err += res[i]; + const bool c_nan = std::isnan(c), r_nan = std::isnan(r); + const bool c_inf = std::isinf(c), r_inf = std::isinf(r); + const bool sign_ne = + (hwy::BitCastScalar(c) < 0) != (hwy::BitCastScalar(r) < 0); + if ((c_nan && r_nan) || (c_inf && r_inf && !sign_ne)) err = 0.0; + if ((c_nan != r_nan) || (c_inf != r_inf) || (c_inf && r_inf && sign_ne)) + err = std::numeric_limits::infinity(); + out[i] = err; + } +} +#endif // HWY_HAVE_FLOAT64 + +// HWY_EXPORT needs one non-template symbol per dispatch table. +void KernelUlpErrorF32(const float* HWY_RESTRICT comp, + const float* HWY_RESTRICT ref, + const double* HWY_RESTRICT res, double* HWY_RESTRICT out, + size_t n) { + LoopUlpError(comp, ref, res, out, n); +} + +void KernelUlpErrorF64(const double* HWY_RESTRICT comp, + const double* HWY_RESTRICT ref, + const double* HWY_RESTRICT res, double* HWY_RESTRICT out, + size_t n) { + LoopUlpError(comp, ref, res, out, n); +} + +} // namespace npsr::py::HWY_NAMESPACE +HWY_AFTER_NAMESPACE(); + +#if HWY_ONCE +namespace npsr::py { +HWY_EXPORT(KernelUlpError32); +HWY_EXPORT(KernelUlpErrorF32); +HWY_EXPORT(KernelUlpErrorF64); + +namespace { +template +using UlpFn = void (*)(const T*, const T*, const double*, double*, size_t); + +template +Array UlpError(UlpFn fn, Array comp, Array ref, + std::optional> residual) { + if (comp.size() != ref.size()) + throw pb::value_error("computed/ref size mismatch"); + if (residual && residual->size() != ref.size()) + throw pb::value_error("residual size mismatch"); + const CArray c = Contiguous(comp), r = Contiguous(ref); + const std::optional> rs = + residual ? std::optional(Contiguous(*residual)) : std::nullopt; + Array out(c.request().shape); + const size_t n = static_cast(c.size()); + { + pb::gil_scoped_release nogil; + fn(c.data(), r.data(), rs ? rs->data() : nullptr, out.mutable_data(), n); + } + return out; +} +} // namespace + +void BindULP(Module& m) { + m.def( + "ulp_error32", + [](Array comp, Array oracle) { + if (comp.size() != oracle.size()) + throw pb::value_error("computed/oracle size mismatch"); + const CArray c = Contiguous(comp); + const CArray o = Contiguous(oracle); + Array out(o.request().shape); + const size_t n = static_cast(o.size()); + { + pb::gil_scoped_release nogil; + HWY_DYNAMIC_DISPATCH(KernelUlpError32) + (c.data(), o.data(), out.mutable_data(), n); + } + return out; + }, + SafeArg("computed"), SafeArg("oracle")); + + m.def( + "ulp_error", + [](Array comp, Array ref, + std::optional> residual) { + return UlpError(HWY_DYNAMIC_POINTER(KernelUlpErrorF32), comp, + ref, residual); + }, + SafeArg("computed"), SafeArg("ref"), SafeArg("residual") = pb::none()); + m.def( + "ulp_error", + [](Array comp, Array ref, + std::optional> residual) { + return UlpError(HWY_DYNAMIC_POINTER(KernelUlpErrorF64), comp, + ref, residual); + }, + SafeArg("computed"), SafeArg("ref"), SafeArg("residual") = pb::none()); +} +} // namespace npsr::py +#endif // HWY_ONCE diff --git a/python/argred.py b/python/argred.py new file mode 100644 index 0000000..907aa00 --- /dev/null +++ b/python/argred.py @@ -0,0 +1,193 @@ +"""Per binade, the float nearest a grid point: trig reduction worst cases.""" + +from __future__ import annotations + +import contextlib +import enum +import functools +import math +from collections.abc import Iterator +from fractions import Fraction +from typing import Any, cast + +import numpy as np +import numpy.typing as npt + +__all__ = ["Argred", "Grid"] + +# floor(pi * 2**1300). +_PI_1300 = int( + "3243f6a8885a308d313198a2e03707344a4093822299f31d0082efa98ec4e6c8" + "9452821e638d01377be5466cf34e90c6cc0ac29b7c97c50dd3f84d5b5b547091" + "79216d5d98979fb1bd1310ba698dfb5ac2ffd72dbd01adfb7b8e1afed6a267e9" + "6ba7c9045f12c7f9924a19947b3916cf70801f2e2858efc16636920d871574e6" + "9a458fea3f4933d7e0d95748f728eb658718bcd5882154aee7b54a41dc25a59b" + "59c30d", + 16, +) + +# 2**1300 / _PI_1300 == 1/pi, exactly representable as a rational. +_INV_PI = Fraction(1 << 1300, _PI_1300) +_INV_2PI = _INV_PI / 2 # 1/(2*pi) + + +class Grid(enum.Enum): + """Reduction grid (inv, offset): grid points are (k + offset)/inv.""" + + PI = (_INV_PI, Fraction(0)) # k*pi -- sin roots + PI_HALF = (_INV_PI, Fraction(1, 2)) # (k+1/2)*pi -- cos roots + TWO_PI = (_INV_2PI, Fraction(0)) # k*2*pi -- full period + + @property + def inv(self) -> Fraction: + return self.value[0] + + @property + def offset(self) -> Fraction: + return self.value[1] + + +def _frac(q: Fraction) -> Fraction: + """Fractional part in [0, 1).""" + return q - math.floor(q) + + +def _ceil(q: Fraction) -> int: + return -int(-q // 1) + + +def _min_frac(a: Fraction, b: Fraction, n: int) -> tuple[Fraction, int]: + """Exact min of {b - t*a mod 1} over 0 <= t < n; returns (d, t).""" + x = _frac(a) + y = 1 - x + d = _frac(b) + u = v = 1 + t = 0 + while x > 0: + if d >= x: + m = min(d // x, (n - 1 - t) // v) # every record reachable with step v + d -= m * x + t += m * v + if d >= x: + break + if v > n - 1 - t: + break + if y > x: # refine the above-step below x so the below-ladder can extend + k = _ceil((y - x) / x) + y -= k * x + u += k * v + # x - i*y stays a valid step only while positive + i = min(_ceil((x - d) / y), _ceil(x / y) - 1) + if i <= 0: + break + x -= i * y + v += i * u + return d, t + + +def _grid_distance(t: int, e: int, inv: Fraction, offset: Fraction, p: int) -> Fraction: + """Distance of x*inv to the grid Z + offset (x = binade e, mantissa t).""" + x = Fraction((1 << (p - 1)) + t) * Fraction(2) ** (e - p) + r = _frac(x * inv - offset) + return min(r, 1 - r) + + +def _binade_candidate( + e: int, inv: Fraction, offset: Fraction, above: bool, p: int +) -> int: + """Mantissa offset minimizing x*inv's distance below (or above) the grid.""" + a = _frac(Fraction(2) ** (e - p) * inv) + base = Fraction(2) ** (e - 1) * inv - offset + b = _frac(-base) + if above: + a, b = _frac(-a), _frac(-b) + _, t = _min_frac(a, b, 1 << (p - 1)) + return t + + +def _binade_has_point(e: int, inv: Fraction, offset: Fraction, emax: int) -> bool: + """Whether [2**(e-1), 2**e) holds a positive grid point (k + offset)/inv.""" + if not 1 <= e <= emax: + return False + k = _ceil(Fraction(2) ** (e - 1) * inv - offset) # first k with x >= 2**(e-1) + return k + offset < Fraction(2) ** e * inv + + +@functools.cache +def _champion_bits(e: int, grid: Grid, precision: int, emax: int) -> int: + """IEEE bits of the float in [2**(e-1), 2**e) closest to a `grid` point.""" + inv, offset = grid.inv, grid.offset + if not _binade_has_point(e, inv, offset, emax): + raise ValueError( + f"binade [{2 ** (e - 1)}, {2**e}) has no {grid.name} grid point" + ) + t = min( + ( + _binade_candidate(e, inv, offset, above, precision) + for above in (False, True) + ), + key=lambda t: _grid_distance(t, e, inv, offset, precision), + ) + # A true champion sits within ~2**-(p-1) of the grid, far under this slack. + assert 0 <= t < 1 << (precision - 1) + assert _grid_distance(t, e, inv, offset, precision) < Fraction(2) ** ( + 13 - precision + ) + return ((e + emax - 2) << (precision - 1)) | t + + +class Argred: + """Worst-case reduction inputs for one float type and grid.""" + + def __init__(self, dtype: npt.DTypeLike, grid: Grid = Grid.PI) -> None: + self.grid = grid + # DTypeLike widens to non-float dtypes; narrow to floating for finfo/.type. + self.dtype = cast("np.dtype[np.floating[Any]]", np.dtype(dtype)) + fi = np.finfo(self.dtype) + self.precision = int(fi.nmant) + 1 # incl. implicit bit + self._emax = int(fi.maxexp) + self._uint = np.dtype(f"uint{self.dtype.itemsize * 8}") + + def valid_binades(self) -> Iterator[int]: + """Binade exponents whose [2**(e-1), 2**e) contains a grid point.""" + inv, offset = self.grid.inv, self.grid.offset + return ( + e + for e in range(1, self._emax + 1) + if _binade_has_point(e, inv, offset, self._emax) + ) + + def worst_case(self, e: int) -> float: + """Closest float to a grid point in binade `e` (from `valid_binades`).""" + bits = _champion_bits(e, self.grid, self.precision, self._emax) + return float(np.array(bits, self._uint).view(self.dtype)) + + def grid_points(self, e: int, budget: int | None = None) -> list[float]: + """Nearest float to each grid point in binade `e`; `budget` samples evenly.""" + inv, offset = self.grid.inv, self.grid.offset + if not 1 <= e <= self._emax: + return [] + klo = _ceil(Fraction(2) ** (e - 1) * inv - offset) + khi = _ceil(Fraction(2) ** e * inv - offset) # exclusive + if budget is not None and khi - klo > budget: + span = khi - 1 - klo + ks: range | list[int] = ( + [klo] + if budget < 2 + else [klo + span * i // (budget - 1) for i in range(budget)] + ) + else: + ks = range(klo, khi) + pts = [] + to_dtype = self.dtype.type + with np.errstate(over="ignore"): + for k in ks: + x = (k + offset) / inv # exact grid position + if x <= 0: + continue + # overflow past this type's finite max: no window there, skip. + with contextlib.suppress(OverflowError): + xr = to_dtype(float(x)) + if np.isfinite(xr): + pts.append(float(xr)) + return pts diff --git a/python/conftest.py b/python/conftest.py new file mode 100644 index 0000000..7b2c688 --- /dev/null +++ b/python/conftest.py @@ -0,0 +1,45 @@ +"""Pytest plugin for numpy_sr; the re-exports below are the whole job.""" + +# Absolute imports: `spin test` runs pytest --pyargs from site-packages. + +from numpy_sr._conftest.fixtures import ( + accuracies, + accuracy, + report_ulp, + rng, + target, + targets, + worst_book, +) +from numpy_sr._conftest.hooks import ( + pytest_collection_modifyitems, + pytest_configure, +) +from numpy_sr._conftest.options import pytest_addoption +from numpy_sr._conftest.report import ( + pytest_report_header, + pytest_sessionfinish, + pytest_terminal_summary, +) +from numpy_sr._conftest.xdist import ( + pytest_configure_node, + pytest_testnodedown, +) + +__all__ = [ + "accuracies", + "accuracy", + "pytest_addoption", + "pytest_collection_modifyitems", + "pytest_configure", + "pytest_configure_node", + "pytest_report_header", + "pytest_sessionfinish", + "pytest_terminal_summary", + "pytest_testnodedown", + "report_ulp", + "rng", + "target", + "targets", + "worst_book", +] diff --git a/python/fpenv.py b/python/fpenv.py new file mode 100644 index 0000000..96baa2b --- /dev/null +++ b/python/fpenv.py @@ -0,0 +1,37 @@ +"""Scoped capture of the thread's floating-point exception flags.""" + +from __future__ import annotations + +from typing import Final + +# The stub omits the internal flags-only class; this one shadows it. +from ._numpy_sr import FPEnv as _flags # pyright: ignore[reportAttributeAccessIssue] + + +class FPEnv: + """Clear the FP exception flags on entry, capture on exit; none restored.""" + + INVALID: Final = _flags.INVALID + DIVBYZERO: Final = _flags.DIVBYZERO + OVERFLOW: Final = _flags.OVERFLOW + UNDERFLOW: Final = _flags.UNDERFLOW + INEXACT: Final = _flags.INEXACT + ERRORS: Final = _flags.ERRORS + + clear = staticmethod(_flags.clear) + test = staticmethod(_flags.test) + + def __init__(self) -> None: + self.raised = 0 + + @property + def errors(self) -> int: + return self.raised & self.ERRORS + + def __enter__(self) -> FPEnv: + self.raised = 0 + _flags.clear() + return self + + def __exit__(self, *exc: object) -> None: + self.raised = _flags.test() diff --git a/python/ieee754.py b/python/ieee754.py new file mode 100644 index 0000000..52b52db --- /dev/null +++ b/python/ieee754.py @@ -0,0 +1,220 @@ +"""Bit-level helpers for exhaustive IEEE-754 testing.""" + +from __future__ import annotations + +import functools +from collections.abc import Sequence +from typing import Any, cast + +import numpy as np +import numpy.typing as npt + +Float = np.floating[Any] +FloatArray = npt.NDArray[Float] +# A magnitude bound: numpy float scalar or plain float. +Bound = float | Float +Binades = Sequence[tuple[Bound, Bound]] + +_SIGNED = {np.dtype(np.float32): np.int32, np.dtype(np.float64): np.int64} + + +def _float_dtype(dtype: npt.DTypeLike) -> np.dtype[np.floating[Any]]: + """Narrow ``DTypeLike`` to a concrete float dtype for type-checkers.""" + return cast("np.dtype[np.floating[Any]]", np.dtype(dtype)) + + +def bits(a: np.ndarray) -> np.ndarray: + """Reinterpret floats as same-width unsigned ints for bit-exact compares.""" + a = np.asarray(a) + return a.view(np.dtype(f"uint{a.dtype.itemsize * 8}")) + + +def float_range( + dtype: npt.DTypeLike, + lo: float | Float, + hi: float | Float, +) -> FloatArray: + """Every representable value in [0 <= lo, hi], ascending and inclusive.""" + dt = _float_dtype(dtype) + uint = np.dtype(f"uint{dt.itemsize * 8}") + a = int(np.asarray(lo, dt).view(uint)) + b = int(np.asarray(hi, dt).view(uint)) + if a > b: + raise ValueError("need lo <= hi, and both non-negative") + return np.arange(a, b + 1, dtype=uint).view(dt) + + +@functools.cache +def binades( + dtype: npt.DTypeLike, + step: int = 1, + subnormals: bool = False, +) -> tuple[tuple[Float, Float], ...]: + """(lo, hi) pairs covering the finite positives, `step` binades per pair.""" + # Tuple, not list: the result is cached and handed to every caller. + dt = _float_dtype(dtype) + fi = np.finfo(dt) # raises on non-float dtypes + uint = np.dtype(f"uint{dt.itemsize * 8}") + bias = 1 - fi.minexp + + lows = [(e + bias) << fi.nmant for e in range(fi.minexp, fi.maxexp, step)] + if subnormals: + lows.insert(0, 0) # +0.0 + highs = [*lows[1:], (fi.maxexp + bias) << fi.nmant] # sentinel is +inf + highs = [b - 1 for b in highs] # one ULP below + + def as_float(bits: list[int]) -> FloatArray: + return np.array(bits, uint).view(dt) + + return tuple(zip(as_float(lows), as_float(highs), strict=True)) + + +def both_signs(x: np.ndarray) -> np.ndarray: + """Mirror to both signs and drop the non-finites bit-stepping can spawn.""" + x = np.concatenate([x, -x]) + return x[np.isfinite(x)] + + +def ulp_neighbors(x: np.ndarray, radius: int = 2) -> np.ndarray: + """Each finite positive x plus its +-radius ULP neighbors, by bit stepping.""" + i = np.ascontiguousarray(x).view(_SIGNED[x.dtype]) + off = np.arange(-radius, radius + 1, dtype=i.dtype) + return (i[:, None] + off).ravel().view(x.dtype) + + +def _uint(dtype: npt.DTypeLike) -> np.dtype: + return np.dtype(f"uint{np.dtype(dtype).itemsize * 8}") + + +def _neg(rng: np.random.Generator, n: int, sign: str) -> np.ndarray: + """`n`-long 0/1 is-negative mask per `sign` (positive/negative/mixed).""" + if sign == "mixed": + return rng.integers(0, 2, size=n) + if sign == "positive": + return np.zeros(n, dtype=np.int64) + if sign == "negative": + return np.ones(n, dtype=np.int64) + raise ValueError(f"sign must be 'positive', 'negative', or 'mixed', got {sign!r}") + + +def _signs( + rng: np.random.Generator, dtype: npt.DTypeLike, n: int, sign: str +) -> np.ndarray: + """`n` sign bits already shifted into the MSB, as uints, per `sign`.""" + u = _uint(dtype) + return _neg(rng, n, sign).astype(u) << (np.finfo(dtype).bits - 1) + + +def _normals( + rng: np.random.Generator, + dtype: npt.DTypeLike, + n: int, + binades: Binades, + sign: str, +) -> np.ndarray: + """`n` floats log-uniform over `binades`; `lo == 0` clamps to subnormal.""" + tiny = float(np.finfo(dtype).smallest_subnormal) + lows = np.maximum([float(lo) for lo, _ in binades], tiny) + highs = np.array([float(hi) for _, hi in binades]) + idx = rng.integers(0, len(binades), size=n) + e = rng.uniform(np.log2(lows[idx]), np.log2(highs[idx])) + s = np.where(_neg(rng, n, sign), -1.0, 1.0) + return (s * 2.0**e).astype(dtype) + + +def _subnormals( + rng: np.random.Generator, dtype: npt.DTypeLike, n: int, sign: str +) -> np.ndarray: + """Nonzero subnormals: random mantissa in [1, 2^nmant).""" + u = _uint(dtype) + m = rng.integers(1, 1 << np.finfo(dtype).nmant, size=n, dtype=u) + return (m | _signs(rng, dtype, n, sign)).view(dtype) + + +def _qnans( + rng: np.random.Generator, dtype: npt.DTypeLike, n: int, sign: str +) -> np.ndarray: + """Quiet NaNs: all-ones exponent, quiet bit set, random payload.""" + fi = np.finfo(dtype) + u = _uint(dtype) + quiet = u.type(1) << (fi.nmant - 1) + exp = ((u.type(1) << (fi.bits - 1 - fi.nmant)) - 1) << fi.nmant + payload = rng.integers(0, quiet, size=n, dtype=u) + return (exp | quiet | payload | _signs(rng, dtype, n, sign)).view(dtype) + + +def _snans( + rng: np.random.Generator, dtype: npt.DTypeLike, n: int, sign: str +) -> np.ndarray: + """Signaling NaNs: quiet bit clear, nonzero payload.""" + fi = np.finfo(dtype) + u = _uint(dtype) + quiet = u.type(1) << (fi.nmant - 1) + exp = ((u.type(1) << (fi.bits - 1 - fi.nmant)) - 1) << fi.nmant + payload = rng.integers(1, quiet, size=n, dtype=u) + return (exp | payload | _signs(rng, dtype, n, sign)).view(dtype) + + +def _infs( + rng: np.random.Generator, dtype: npt.DTypeLike, n: int, sign: str +) -> np.ndarray: + """+-inf, by bit assembly.""" + fi = np.finfo(dtype) + u = _uint(dtype) + exp = ((u.type(1) << (fi.bits - 1 - fi.nmant)) - 1) << fi.nmant + return (exp | _signs(rng, dtype, n, sign)).view(dtype) + + +def _zeros( + rng: np.random.Generator, dtype: npt.DTypeLike, n: int, sign: str +) -> np.ndarray: + """+-0.0, by bit assembly.""" + return _signs(rng, dtype, n, sign).view(dtype) + + +_SPECIALS = { + "subnormal": _subnormals, + "zero": _zeros, + "qnan": _qnans, + "snan": _snans, + "inf": _infs, +} + + +def rand( + rng: np.random.Generator, + dtype: npt.DTypeLike, + binades: Binades, + n: int, + *, + sign: str = "mixed", + subnormal: float = 0.0, + zero: float = 0.0, + qnan: float = 0.0, + snan: float = 0.0, + inf: float = 0.0, +) -> FloatArray: + """`n` shuffled values: log-uniform normals, each keyword its block share.""" + dt = _float_dtype(dtype) + fracs = { + "subnormal": subnormal, + "zero": zero, + "qnan": qnan, + "snan": snan, + "inf": inf, + } + parts = [] + used = 0 + for name, frac in fracs.items(): + if frac <= 0.0: + continue + k = min(round(frac * n), n - used) + if k <= 0: + continue + parts.append(_SPECIALS[name](rng, dt, k, sign)) + used += k + if used < n: + parts.append(_normals(rng, dt, n - used, binades, sign)) + out = np.concatenate(parts).astype(dt, copy=False) + rng.shuffle(out) + return out diff --git a/python/meson.build b/python/meson.build new file mode 100644 index 0000000..9005c4c --- /dev/null +++ b/python/meson.build @@ -0,0 +1,41 @@ +py = import('python').find_installation() +# install_sources flattens to basenames: tests/__init__.py would clobber the +# package __init__.py if shared with this call, so tests get their own. +py.install_sources( + [ + '__init__.py', + 'argred.py', + 'fpenv.py', + 'ieee754.py', + 'testing.py', + 'conftest.py', + '_numpy_sr.pyi', # type stub for the compiled extension (sibling of the .so) + ], + subdir: 'numpy_sr', +) +# conftest.py is a thin aggregator; its implementation lives here, one module +# per concern (install_sources flattens to basenames, so its own subdir). +py.install_sources( + [ + '_conftest/__init__.py', + '_conftest/book.py', + '_conftest/fixtures.py', + '_conftest/hooks.py', + '_conftest/options.py', + '_conftest/report.py', + '_conftest/stash.py', + '_conftest/targets.py', + '_conftest/xdist.py', + ], + subdir: 'numpy_sr/_conftest', +) +py.install_sources( + [ + 'tests/__init__.py', + 'tests/test_trig.py', + 'tests/test_worst.py', + ], + subdir: 'numpy_sr/tests', +) +subdir('_numpy_sr') + diff --git a/python/testing.py b/python/testing.py new file mode 100644 index 0000000..a955f50 --- /dev/null +++ b/python/testing.py @@ -0,0 +1,173 @@ +"""Fractional-ULP error vs the exact value; correct rounding scores ~0.5.""" + +from __future__ import annotations + +import math +from collections.abc import Callable +from typing import Any + +import numpy as np +import numpy_sr as sr + +Sink = Callable[[dict[str, Any]], None] + +# float dtype -> same-width signed / unsigned int, for bit reinterpretation. +_UNSIGNED = {np.dtype(np.float32): np.uint32, np.dtype(np.float64): np.uint64} + + +def ulp_error(computed: Any, ref: Any) -> np.ndarray: + """Signed fractional-ULP error; non-finite mismatches score +inf.""" + ref, residual = _unpack_ref(ref) + a = np.asarray(computed) + # sr.ulp_error is dtype-exact (SafeArg); only layout is its problem. + b = np.asarray(ref).astype(a.dtype, copy=False) + a, b = np.broadcast_arrays(a, b) + res = np.asarray(residual, dtype=np.float64) + r = np.broadcast_to(res, a.shape) if res.ndim or res != 0.0 else None + return sr.ulp_error(a, b, r) + + +def _unpack_ref(ref: Any) -> tuple[Any, Any]: + """Split `ref` into `(ref, residual)`; a bare reference gets 0.0.""" + if isinstance(ref, tuple) and len(ref) == 2: + return ref + return ref, 0.0 + + +def _num(v: Any) -> float | str: + """JSON-safe number: non-finite floats become tokens strict JSON accepts.""" + f = float(v) + if math.isnan(f): + return "nan" + if math.isinf(f): + return "inf" if f > 0 else "-inf" + return f + + +def _fmt(v: np.ndarray) -> dict[str, Any]: + """A float element as both its decimal value and exact IEEE bit pattern.""" + dt = v.dtype + bits = int(v.view(_UNSIGNED[dt])) + return {"value": _num(v), "bits": f"0x{bits:0{dt.itemsize * 2}x}"} + + +def _offenders( + x: np.ndarray, + computed: np.ndarray, + ref: np.ndarray, + residual: np.ndarray, + err: np.ndarray, + limit: int, +) -> list[dict[str, Any]]: + order = np.argsort(-np.abs(err), kind="stable")[:limit] + return [ + { + "x": _fmt(x[i]), + "actual": _fmt(computed[i]), + "expected": _fmt(ref[i]), + "residual": float(residual[i]), + "error": _num(err[i]), + } + for i in map(int, order) + ] + + +def measure( + computed: Any, + ref: Any, + maxulp: float, + *, + x: Any, + label: str = "", + limit: int = 16, + **meta: Any, +) -> dict[str, Any]: + """Compare and build a record without asserting; extra kwargs go in it.""" + ref, residual = _unpack_ref(ref) + a = np.asarray(computed) + b = np.asarray(ref).astype(a.dtype, copy=False) + a, b = np.atleast_1d(*np.broadcast_arrays(a, b)) + x = np.broadcast_to(np.asarray(x).astype(a.dtype, copy=False), a.shape) + res = np.broadcast_to(np.asarray(residual, dtype=np.float64), a.shape) + + err = ulp_error(a, (b, res)) + abserr = np.abs(err) + n = int(abserr.size) + finite = np.isfinite(abserr) + n_over = int(np.count_nonzero(abserr > maxulp)) # +inf mismatches count here + + record = { + "label": label, + "dtype": str(a.dtype), + "n": n, + "maxulp": float(maxulp), + "max_error": _num(abserr.max()) if n else 0.0, + "mean_error": float(abserr[finite].mean()) if finite.any() else 0.0, + "n_over": n_over, + "passed": n_over == 0, + "worst": _offenders(x, a, b, res, err, limit), + } + record.update(meta) + return record + + +def _exceeds(error: float | str, maxulp: float) -> bool: + # _num() returns a string only for a non-finite (always a violation). + return isinstance(error, str) or abs(error) > maxulp + + +def _format_failure(record: dict[str, Any]) -> str: + head = ( + f"{record['label']}: {record['n_over']}/{record['n']} over " + f"{record['maxulp']} ULP (max {record['max_error']}, {record['dtype']})" + ) + rows = [ + f" x={o['x']['value']} ({o['x']['bits']}) " + f"got={o['actual']['value']} ref={o['expected']['value']} " + f"residual={o['residual']:+.4f} err={o['error']}" + for o in record["worst"] + if _exceeds(o["error"], record["maxulp"]) + ] + return "\n".join([head, *rows]) + + +def assert_max_ulp( + computed: Any, + ref: Any, + maxulp: float, + *, + x: Any, + label: str = "", + limit: int = 16, + sink: Sink | None = None, + **meta: Any, +) -> dict[str, Any]: + """Assert every element scores <= `maxulp`; `sink` sees it either way.""" + __tracebackhide__ = True + record = measure( + computed, + ref, + maxulp, + x=x, + label=label, + limit=limit, + **meta, + ) + if sink is not None: + sink(record) + if not record["passed"]: + raise AssertionError(_format_failure(record)) + return record + + +def assert_bits(got: np.ndarray, want: np.ndarray, x: np.ndarray, msg: str): + """Fail unless `got` and `want` are bit-for-bit identical.""" + __tracebackhide__ = True + bad = sr.ieee754.bits(got) != sr.ieee754.bits(want) + if bad.any(): + i = int(np.argmax(bad)) + raise AssertionError( + f"{msg}: {int(bad.sum())}/{bad.size} bit mismatches; first at " + f"x={x[i]!r} (0x{int(sr.ieee754.bits(x)[i]):x}) " + f"got {got[i]!r} want {want[i]!r}" + ) diff --git a/python/tests/__init__.py b/python/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/python/tests/test_trig.py b/python/tests/test_trig.py new file mode 100644 index 0000000..c00fdfe --- /dev/null +++ b/python/tests/test_trig.py @@ -0,0 +1,245 @@ +"""sin/cos accuracy, IEEE specials, math identities and fenv discipline.""" + +import functools +import math +from collections.abc import Sequence +from itertools import product +from typing import Any, cast + +import numpy as np +import numpy_sr as sr +from numpy_sr import ieee754 as ie +from numpy_sr.tests.test_worst import register as _register_worst +from pytest import mark, param, skip + + +def _bound_low(target, dtype): + if dtype is np.float32: + return 2.0 + return 3.5 if target.HAVE_FMA else 4.0 + + +def _bound_high(target, dtype): + return 0.6 if dtype is np.float32 else 1.0 + + +_BOUNDS = { + sr.Accuracy.High: _bound_high, + sr.Accuracy.Low: _bound_low, +} + +_OPS = ["sin", "cos"] +_DTYPES = [np.float32, np.float64] +_MIN_BLOCK = 4096 + + +# Enroll sin/cos in the shared frozen worst-case gate (test_worst.py). +for _op in _OPS: + _register_worst(_op, lambda target, dtype, acc: _BOUNDS[acc](target, dtype)) + + +def _sweeps_binades(dtypes_jobs, subnormals=True): + """Binades dealt round-robin into `jobs` cases per dtype. + + Splitting serves two ends: the cases run in parallel under `pytest -n`, and + each books its own worst cases, widening the range the book covers. + """ + for dtype, jobs in dtypes_jobs: + bins = ie.binades(dtype, subnormals=subnormals) + for i in range(jobs): + yield param(dtype, bins[i::jobs], id=f"{dtype.__name__}-{i}") + + +def _report_combos(report_ulp, operation, dtype, x, ref, combos, store_worst=None): + """Score one block against every combo; `ref` is one shared MPFR pass.""" + for target, accuracy in combos: + computed = getattr(target, operation)(x, accuracy=accuracy) + bound = _BOUNDS[accuracy](target, dtype) + report_ulp( + x, + computed, + ref, + bound, + accuracy=accuracy, + target=target, + store_worst=store_worst, + ) + + +@functools.cache +def _basic(operation, dtype) -> tuple[np.ndarray, Any]: + """Named args over the quadrant seams, both signs, +-2 ULP neighbors.""" + fi, pi = np.finfo(dtype), math.pi + named = [0.0, pi / 6, pi / 4, pi / 3, pi / 2, 2 * pi / 3] + named += [pi, 3 * pi / 2, 2 * pi, fi.tiny, fi.smallest_subnormal] + x = ie.both_signs(ie.ulp_neighbors(np.array(named, dtype=dtype))) + return x, getattr(sr.mpfr, operation)(x) + + +@mark.parametrize("operation", _OPS) +@mark.parametrize("dtype", _DTYPES) +@mark.parametrize("data", [_basic]) +def test_simple(report_ulp, target, accuracy, operation, dtype, data): + op_fn = getattr(target, operation) + x, ref = data(operation, dtype) + computed = op_fn(x, accuracy=accuracy) + report_ulp(x, computed, ref, _BOUNDS[accuracy](target, dtype)) + + +def _seed_key(rng: np.random.Generator) -> tuple: + """A Generator hashes by identity; its seed is what the data depends on.""" + seed_seq = cast("np.random.SeedSequence", rng.bit_generator.seed_seq) + entropy = seed_seq.entropy + return tuple(entropy) if isinstance(entropy, Sequence) else (entropy,) + + +@functools.cache +def _rand(seed, operation, dtype, binades) -> tuple[np.ndarray, Any]: + """Cached on `_seed_key`, so every target shares one block and one MPFR pass.""" + x = ie.rand(np.random.default_rng(list(seed)), dtype, binades, 2**16) + return x, getattr(sr.mpfr, operation)(x) + + +@mark.parametrize("operation", _OPS) +@mark.parametrize( + "dtype,binades", list(_sweeps_binades([(np.float32, 2), (np.float64, 8)])) +) +def test_rand(report_ulp, rng, target, accuracy, operation, dtype, binades): + op_fn = getattr(target, operation) + x, ref = _rand(_seed_key(rng), operation, dtype, binades) + if not x.size: + skip(f"binades lie above the {accuracy.name} cap") + computed = op_fn(x, accuracy=accuracy) + report_ulp(x, computed, ref, _BOUNDS[accuracy](target, dtype)) + + +# Adaptive worst-ULP hunt +_HUNT_ROUNDS = 128 +_HUNT_POP = 2**16 +_HUNT_KEEP = 512 +_HUNT_RADIUS = 3 + + +def _hunt_worst(combos, operation, ref_op, rng, dtype, seeds): + """Hill-climb the ULP landscape: score vs MPFR, keep top-KEEP, re-probe.""" + binades = ie.binades(dtype) + fresh = ie.rand(rng, dtype, binades, _HUNT_POP, sign="positive") + pool = np.concatenate([np.asarray(seeds, dtype=dtype), fresh]) + best: list[tuple[np.ndarray, np.ndarray]] = [ + (np.empty(0, dtype), np.empty(0)) for _ in combos + ] + for _ in range(_HUNT_ROUNDS): + pool = np.unique(pool[np.isfinite(pool) & (pool > 0)]) + ref = ref_op(pool) # one MPFR pass, shared by every combo below + for i, (target, accuracy) in enumerate(combos): + op_fn = getattr(target, operation) + err = np.abs(sr.testing.ulp_error(op_fn(pool, accuracy=accuracy), ref)) + cand_x = np.concatenate([best[i][0], pool]) + cand_e = np.concatenate([best[i][1], err]) + keep = np.argsort(-cand_e, kind="stable")[:_HUNT_KEEP] + best[i] = cand_x[keep], cand_e[keep] + centers = np.unique(np.concatenate([bx for bx, _ in best])) + fresh = ie.rand(rng, dtype, binades, _HUNT_POP, sign="positive") + pool = np.concatenate([centers, ie.ulp_neighbors(centers, _HUNT_RADIUS), fresh]) + return np.unique(np.concatenate([bx for bx, _ in best])) + + +@mark.slow +@mark.parametrize("operation", _OPS) +@mark.parametrize("dtype", [np.float64]) +def test_worst_adaptive(report_ulp, rng, targets, accuracies, operation, dtype): + """`_hunt_worst` seeded from `argred`'s champions, +-2 ULP window at the end.""" + combos = list(product(targets, accuracies)) + ar = sr.argred.Argred( + dtype, sr.argred.Grid.PI if operation == "sin" else sr.argred.Grid.PI_HALF + ) + seeds = [ar.worst_case(e) for e in ar.valid_binades()] + + ref_op = getattr(sr.mpfr, operation) + champions = _hunt_worst(combos, operation, ref_op, rng, dtype, seeds) + + x = ie.both_signs(ie.ulp_neighbors(champions, radius=2)) + _report_combos(report_ulp, operation, dtype, x, ref_op(x), combos) + + +@mark.exhaustive +@mark.parametrize("operation", _OPS) +@mark.parametrize("dtype,binades", list(_sweeps_binades([(np.float32, 8)]))) +def test_exhaustive(report_ulp, targets, accuracies, operation, dtype, binades): + combos = list(product(targets, accuracies)) # reused across binades + for lo, hi in binades: + x = ie.float_range(dtype, lo, hi) + ref = getattr(sr.mpfr, operation)(x) + # the only f32 book worth keeping: this sweep sees every input there is + _report_combos(report_ulp, operation, dtype, x, ref, combos, store_worst=True) + del x, ref + + +@mark.parametrize("operation", _OPS) +@mark.parametrize("dtype", _DTYPES) +def test_parity(rng, target, accuracy, operation, dtype): + """sin odd / cos even: op(-x) must equal +-op(x) bit for bit.""" + op_fn = getattr(target, operation) + binades = ie.binades(dtype) + x = ie.rand(rng, dtype, binades, _MIN_BLOCK, sign="positive", subnormal=0.25) + pos = op_fn(x, accuracy=accuracy) + neg = op_fn(-x, accuracy=accuracy) + want = -pos if operation == "sin" else pos + sr.testing.assert_bits(neg, want, -x, f"{operation}(-x) parity") + + +@mark.parametrize("operation", _OPS) +@mark.parametrize("dtype", _DTYPES) +def test_tiny(rng, target, accuracy, operation, dtype): + """Below the round-to-1 branch, sin(x) == x and cos(x) == 1.0 bit-exactly.""" + # The stricter High round-to-1 threshold, |x| <= half ulp(pi/2). + cap = {np.float32: 2.0**-24, np.float64: 2.0**-53}[dtype] + tiny = np.finfo(dtype).smallest_subnormal + x = ie.rand(rng, dtype, [(tiny, cap)], _MIN_BLOCK, subnormal=0.5) + got = getattr(target, operation)(x, accuracy=accuracy) + want = x if operation == "sin" else np.ones_like(x) + sr.testing.assert_bits(got, want, x, f"{operation} tiny identity") + + +@mark.parametrize("operation", _OPS) +@mark.parametrize("dtype", _DTYPES) +def test_special(rng, target, accuracy, operation, dtype): + """Bit-exact IEEE specials scattered through a normal block.""" + specials = dict.fromkeys(("zero", "subnormal", "qnan", "snan", "inf"), 0.15) + x = ie.rand(rng, dtype, ie.binades(dtype), _MIN_BLOCK, **specials) + got = getattr(target, operation)(x, accuracy=accuracy) + + # sub-tiny lanes (+-0, subnormals): sin preserves, cos snaps to 1.0. + small = np.isfinite(x) & (np.abs(x) < np.finfo(dtype).tiny) + want = x[small] if operation == "sin" else np.ones(int(small.sum()), dtype) + sr.testing.assert_bits(got[small], want, x[small], f"{operation} sub-tiny") + + # inf / qNaN / sNaN all collapse to a quiet NaN. + nonfinite = ~np.isfinite(x) + assert np.isnan(got[nonfinite]).all(), f"{operation}: non-finite input -> non-NaN" + + +# `rand` configs; INVALID is raised solely for inf, so the expected outcome is +# read off the block rather than tabulated per case. +_FLAG_CASES = { + "normals": {}, + "zeros": {"zero": 1.0}, + "subnormals": {"subnormal": 1.0}, + "qnan": {"qnan": 1.0}, + "snan": {"snan": 1.0}, + "inf": {"inf": 1.0}, + "mixed": {"subnormal": 0.1, "qnan": 0.1, "snan": 0.1, "inf": 0.1}, +} + + +@mark.parametrize("operation", _OPS) +@mark.parametrize("dtype", _DTYPES) +@mark.parametrize("case", _FLAG_CASES) +def test_fp_errors(rng, target, accuracy, operation, dtype, case): + """fenv discipline: spurious flags suppressed, INVALID set only for inf.""" + if target.IS_REFERENCE: + skip("kernel fenv not applicable for reference targets") + x = ie.rand(rng, dtype, ie.binades(dtype), _MIN_BLOCK, **_FLAG_CASES[case]) + with sr.FPEnv() as f: + getattr(target, operation)(x, accuracy=accuracy) + assert f.errors == (sr.FPEnv.INVALID if np.isinf(x).any() else 0) diff --git a/python/tests/test_worst.py b/python/tests/test_worst.py new file mode 100644 index 0000000..6faed2c --- /dev/null +++ b/python/tests/test_worst.py @@ -0,0 +1,80 @@ +"""Reusable frozen worst-case gate; an op opts in via `register(op, bound)`.""" + +from collections.abc import Callable + +import numpy as np +import numpy_sr as sr +import pytest +from numpy_sr._conftest.book import WorstBook +from numpy_sr._conftest.stash import BOOK_KEY +from numpy_sr._conftest.targets import ACCURACIES, TARGETS + +# op name -> bound(target, dtype, accuracy) -> max |ULP|. Domains register here. +_OP_BOUNDS: dict[str, Callable] = {} + + +def register(operation: str, bound: Callable) -> None: + """Enroll `operation`; `bound(target, dtype, accuracy)` gives its max |ULP|.""" + _OP_BOUNDS[operation] = bound + + +def pytest_generate_tests(metafunc: pytest.Metafunc) -> None: + """One case per book file, against every target that book applies to. + + `target`/`accuracy` are parametrized here rather than taken from their + fixtures so the key's own axes pick the combinations; the CLI selection + still deselects the rest, which reads those same param names. + """ + if "worst_key" not in metafunc.fixturenames: + return + cases = [] + for key in metafunc.config.stash[BOOK_KEY].all_keys(): + operation, acc, dtype, have_fma, source = WorstBook.parse(key) + accuracy = ACCURACIES.get(acc) + if accuracy is None: # a book for an accuracy this build dropped + continue + for name, target in TARGETS.items(): + if have_fma != target.HAVE_FMA: + continue + if dtype == "float64" and not target.HAVE_FLOAT64: + continue + cases.append( + pytest.param( + target, + accuracy, + key, + id=f"{name}-{accuracy.name}-{operation}-{dtype}-{source}", + ) + ) + metafunc.parametrize("target,accuracy,worst_key", cases) + + +def test_worst_book(request, report_ulp, worst_book, target, accuracy, worst_key): + """Replay one frozen book; a violation aborts the session.""" + operation, _, dtype, _, _ = WorstBook.parse(worst_key) + bound = _OP_BOUNDS.get(operation) + if bound is None: + pytest.skip(f"{operation} has a book but no registered bound") + x = worst_book.replay(worst_key) + if not x.size: + pytest.skip("empty book") + + try: + report_ulp( + x, + getattr(target, operation)(x, accuracy=accuracy), + getattr(sr.mpfr, operation)(x), + bound(target, np.dtype(dtype).type, accuracy), + operation=operation, + target=target, + accuracy=accuracy, + # never folds new inputs; `refresh_key` re-stamps this book's errors + # under --collect-worst, and only for a non-reference target + store_worst=False, + refresh_key=worst_key, + ) + except AssertionError: + # shouldstop is the channel xdist propagates back from a worker; a + # pytest.exit() would abort one worker and leave the others running. + request.session.shouldstop = "frozen worst-case regression" + raise diff --git a/subprojects/highway.wrap b/subprojects/highway.wrap new file mode 100644 index 0000000..cd6fb7a --- /dev/null +++ b/subprojects/highway.wrap @@ -0,0 +1,8 @@ +[wrap-git] +url = https://github.com/google/highway.git +revision = master +depth = 1 + +[provide] +libhwy = hwy_dep +libhwy-test = hwy_test_dep diff --git a/subprojects/packagefiles/svml/meson.build b/subprojects/packagefiles/svml/meson.build new file mode 100644 index 0000000..01a6645 --- /dev/null +++ b/subprojects/packagefiles/svml/meson.build @@ -0,0 +1,26 @@ +project( + 'svml', + 'c', + license: 'BSD-3-Clause', + meson_version: '>=1.1', +) + +svml_srcs = files( + 'linux/avx512/svml_z0_cos_d_ha.s', + 'linux/avx512/svml_z0_cos_d_la.s', + 'linux/avx512/svml_z0_cos_s_ha.s', + 'linux/avx512/svml_z0_cos_s_la.s', + 'linux/avx512/svml_z0_sin_d_ha.s', + 'linux/avx512/svml_z0_sin_d_la.s', + 'linux/avx512/svml_z0_sin_s_ha.s', + 'linux/avx512/svml_z0_sin_s_la.s', +) + +svml_lib = static_library( + 'svml', + svml_srcs, + pic: true, +) + +svml_dep = declare_dependency(link_with: svml_lib) +meson.override_dependency('svml', svml_dep) diff --git a/subprojects/svml.wrap b/subprojects/svml.wrap new file mode 100644 index 0000000..a9a615a --- /dev/null +++ b/subprojects/svml.wrap @@ -0,0 +1,9 @@ +[wrap-git] +directory = svml +url = https://github.com/numpy/SVML.git +revision = main +depth = 1 +patch_directory = svml + +[provide] +svml = svml_dep diff --git a/tools/ulp_report/README.md b/tools/ulp_report/README.md new file mode 100644 index 0000000..b1a5de7 --- /dev/null +++ b/tools/ulp_report/README.md @@ -0,0 +1,33 @@ +# ULP report viewer + +Static viewer for the fractional-ULP JSON report the pytest suite writes +(`ulp-report.json` at the repo root). `spin ulp-report` bundles that JSON into a +single HTML file and also writes a `ulp-report.js` sidecar so `index.html` works +over `file://`. + +- **Local:** run the test suite, then `spin ulp-report` to refresh the sidecar + and open `index.html` directly โ€” no server. A report living elsewhere + (custom `$NPSR_ULP_JSON`, CI download) can be dragged onto the page instead. +- **CI:** `spin ulp-report` (optionally ` -o `) emits one HTML + artifact with the JSON inlined. + +The JSON is a render manifest: `python/_conftest/report.py::_serialize_records` does +every format, ULP ratio, mean and sort-key up front, so each entry in `records` +is already one table row of display strings (plus numeric aids for sorting). + +The page is just [Tabulator](https://tabulator.info) loaded from a CDN: +`viewer.js` maps `payload.records` to columns and lets the grid handle sort and +per-column filters โ€” no framework, no build, no vendored file. **The CDN means +the page needs network access to render** (both `index.html` and the bundled +artifact). + +- **Filters:** facet pills above the table (one group per categorical column) + drive `table.setFilter`, keeping the grid header clean; numeric columns sort + by their `*Num` field so `inf`/`1e+100` rank correctly. +- **Detail:** click a row to expand its per-test breakdown and worst offenders + (value + IEEE bits), built from the row's own `perTest`/`offenders` โ€” the + only custom markup, styled by `viewer.css`. `renderVertical: "basic"` renders + all rows (no virtual scroll) so print and Ctrl-F see everything. + +To pin or upgrade Tabulator, edit the two CDN ` + + + + + diff --git a/tools/ulp_report/viewer.css b/tools/ulp_report/viewer.css new file mode 100644 index 0000000..9846bb4 --- /dev/null +++ b/tools/ulp_report/viewer.css @@ -0,0 +1,152 @@ +/* Fully custom dark theme over Tabulator's structural (base) stylesheet โ€” the + palette and typography from the original viewer, recreated so nothing fights + a bundled theme. */ + +:root { + --page: #0d0d0d; + --surface: #1a1a19; + --ink: #ffffff; + --ink-2: #c3c2b7; + --muted: #898781; + --grid: #2c2c2a; + --border: rgba(255, 255, 255, 0.10); + --accent: #3987e5; /* magnitude bars */ + --good: #35c956; + --bad: #e0564f; + --good-wash: rgba(53, 201, 86, 0.12); + --bad-wash: rgba(224, 86, 79, 0.12); + --sans: system-ui, -apple-system, "Segoe UI", sans-serif; + --mono: ui-monospace, SFMono-Regular, "SF Mono", Menlo, Consolas, monospace; +} + +* { box-sizing: border-box; } + +body { + margin: 0 auto; + padding: 1.5rem; + max-width: 82rem; + background: var(--page); + color: var(--ink); + font: 14px/1.45 var(--sans); +} +h1 { font-size: 1.25rem; margin: 0 0 0.25rem; font-weight: 650; } +.meta { color: var(--ink-2); margin: 0 0 0.25rem; } +.hint { color: var(--muted); margin: 0 0 1rem; font-size: 0.85rem; } +#dropzone { margin: 1rem 0; color: var(--muted); } + +/* ---- facet pills (filtering; keeps the table header clean) ---- */ + +.filters { margin: 0 0 1rem; } +.facet { display: flex; align-items: baseline; gap: 0.5rem; margin: 0.3rem 0; } +.facet .name { + color: var(--muted); font-size: 0.72rem; text-transform: uppercase; + letter-spacing: 0.04em; min-width: 5.5rem; +} +.facet .pills { display: flex; flex-wrap: wrap; gap: 0.3rem; } +.pill { + border: 1px solid var(--grid); border-radius: 999px; + background: var(--surface); color: var(--muted); + font: inherit; font-size: 0.8rem; padding: 0.1rem 0.7rem; cursor: pointer; +} +.pill[aria-pressed="true"] { background: var(--accent); border-color: var(--accent); color: #fff; } + +/* ---- Tabulator: dark surface, our borders, single (browser) scrollbar ---- */ + +.tabulator { + background: var(--surface); + border: 1px solid var(--border); + border-radius: 10px; + color: var(--ink); + font: 14px/1.45 var(--sans); +} +.tabulator, .tabulator .tabulator-tableholder { height: auto !important; max-height: none !important; } +.tabulator .tabulator-tableholder { overflow: visible !important; } + +.tabulator .tabulator-header { + background: var(--surface); + border-bottom: 1px solid var(--grid); + color: var(--muted); +} +.tabulator .tabulator-header .tabulator-col { + background: var(--surface); + border-right: none; +} +.tabulator .tabulator-header .tabulator-col .tabulator-col-content { padding: 0.5rem 0.75rem; } +.tabulator .tabulator-header .tabulator-col .tabulator-col-title { + color: var(--muted); + font-size: 0.72rem; + text-transform: uppercase; + letter-spacing: 0.04em; + font-weight: 650; +} +.tabulator .tabulator-header .tabulator-col.tabulator-sortable.tabulator-col-sorter-element:hover, +.tabulator .tabulator-header .tabulator-col.tabulator-sortable:hover { + background-color: var(--page) !important; /* match the row hover, not the base gray */ +} +.tabulator .tabulator-col .tabulator-arrow { border-bottom-color: var(--muted) !important; } +.tabulator .tabulator-col.tabulator-sortable[aria-sort="ascending"] .tabulator-arrow, +.tabulator .tabulator-col.tabulator-sortable[aria-sort="descending"] .tabulator-arrow { border-bottom-color: var(--ink) !important; } + +.tabulator .tabulator-row, +.tabulator .tabulator-row.tabulator-row-even { background: var(--surface); } +.tabulator .tabulator-row { border-top: 1px solid var(--grid); } +.tabulator .tabulator-row:first-child { border-top: none; } +.tabulator .tabulator-row.tabulator-selectable:hover, +.tabulator .tabulator-row.tabulator-selected:hover, +.tabulator .tabulator-row:hover { background-color: var(--page) !important; } +.tabulator .tabulator-row .tabulator-cell { + color: var(--ink); + border: none; + padding: 0.5rem 0.75rem; +} +.tabulator-cell.mono { font-family: var(--mono); font-size: 0.9em; font-variant-numeric: tabular-nums; } + +/* ---- badges (fma / status) ---- */ + +.badge { + display: inline-block; + border-radius: 5px; + font-size: 0.75rem; + font-weight: 600; + padding: 0.05rem 0.45rem; +} +.badge.pass { color: var(--good); background: var(--good-wash); } +.badge.fail { color: var(--bad); background: var(--bad-wash); } +.badge.fma { color: var(--ink-2); border: 1px solid var(--border); font-weight: 400; } + +/* ---- max-error magnitude bar ---- */ + +.bar-wrap { display: flex; align-items: center; gap: 0.6rem; } +.bar-wrap .mono { + min-width: 3.4rem; text-align: right; + font-family: var(--mono); font-size: 0.9em; font-variant-numeric: tabular-nums; +} +.bar { + flex: 1; height: 6px; min-width: 30px; max-width: 4.2rem; + border-radius: 4px; background: var(--grid); overflow: hidden; +} +.bar > i { display: block; height: 100%; background: var(--accent); border-radius: 4px; } +.bar > i.over { background: var(--bad); } + +/* ---- expandable detail panel (value over bits, like the old viewer) ---- */ + +.detail { + flex: 0 0 100%; + padding: 0.55rem 0.85rem 0.8rem; + background: var(--page); + border-top: 1px solid var(--grid); + font-size: 0.85rem; + overflow-x: auto; +} +.detail .sub { width: 100%; border-collapse: collapse; } +.detail .sub th { + text-align: left; color: var(--muted); font-weight: 650; + text-transform: uppercase; font-size: 0.68rem; letter-spacing: 0.04em; +} +.detail .sub th.num, .detail .sub td.num { text-align: right; } +.detail .sub th, .detail .sub td { padding: 0.2rem 0.7rem; white-space: nowrap; vertical-align: top; } +.detail .sub td { font-family: var(--mono); color: var(--ink); } +.detail .v { display: block; } +.detail .bits { display: block; color: var(--muted); } +.detail .over td { color: var(--bad); } +.detail .sub-h { color: var(--muted); margin: 0.6rem 0 0.25rem; font-family: var(--mono); } diff --git a/tools/ulp_report/viewer.js b/tools/ulp_report/viewer.js new file mode 100644 index 0000000..022dedd --- /dev/null +++ b/tools/ulp_report/viewer.js @@ -0,0 +1,172 @@ +/* NPSR ULP report viewer. Takes the JSON and prints it with Tabulator (CDN). + * + * The payload is a render manifest (python/conftest.py::_serialize_records did + * every format, ratio, mean and sort-key): each record is one table row, with + * display strings for text and numeric fields for sorting. Filtering is by + * facet pills above the table (clean header, like the old viewer); a row click + * expands the per-test + worst-offender detail the manifest already carries. + * + * Data sources, strongest first (see index.html): + * 1. inline