From 303daf513a272aefe722643da68c83200b7bef13 Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 00:07:18 -0500 Subject: [PATCH 01/13] Fixed missing _scalar on bernoulli mgf/cgf function calls --- python/fastdist/distributions/bernoulli.py | 4 ++-- python/tests/test_bernoulli.py | 17 +++++++++++++++++ 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/python/fastdist/distributions/bernoulli.py b/python/fastdist/distributions/bernoulli.py index bd87815..25822bd 100644 --- a/python/fastdist/distributions/bernoulli.py +++ b/python/fastdist/distributions/bernoulli.py @@ -559,7 +559,7 @@ def _mgf_scalar(cls, t: Real, p: Real) -> Real: cls._validate_params(p=p) cls._validate_inputs(_input=t, input_name="t") - return _core.bernoulli_mgf(float(t), float(p)) + return _core.bernoulli_mgf_scalar(float(t), float(p)) @classmethod def _cgf_scalar(cls, t: Real, p: Real) -> Real: @@ -588,7 +588,7 @@ def _cgf_scalar(cls, t: Real, p: Real) -> Real: cls._validate_params(p=p) cls._validate_inputs(_input=t, input_name="t") - return _core.bernoulli_cgf(float(t), float(p)) + return _core.bernoulli_cgf_scalar(float(t), float(p)) # ------------------------------------------------------------------------------------------------------------------ # Batch Instance Methods diff --git a/python/tests/test_bernoulli.py b/python/tests/test_bernoulli.py index 914659f..f103a9e 100644 --- a/python/tests/test_bernoulli.py +++ b/python/tests/test_bernoulli.py @@ -51,6 +51,23 @@ def test_pmf_scalar_valid(self, k, p): val = Bernoulli._pmf_scalar(k, p) assert isinstance(val, float) + @pytest.mark.parametrize("k, p, expected", [ + (0, 0.3, 0.7), + (1, 0.3, 1.0), + ]) + def test_cdf_scalar_valid(self, k, p, expected): + assert Bernoulli._cdf_scalar(k, p) == pytest.approx(expected) + + @pytest.mark.parametrize("t, p", [(0.0, 0.3), (1.0, 0.3), (-1.0, 0.8)]) + def test_mgf_scalar_valid(self, t, p): + expected = 1.0 - p + p * np.exp(t) + assert Bernoulli._mgf_scalar(t, p) == pytest.approx(expected) + + @pytest.mark.parametrize("t, p", [(0.0, 0.3), (1.0, 0.3), (-1.0, 0.8)]) + def test_cgf_scalar_valid(self, t, p): + expected = np.log(1.0 - p + p * np.exp(t)) + assert Bernoulli._cgf_scalar(t, p) == pytest.approx(expected) + @pytest.mark.parametrize("method_name, args", [ ("_pmf_scalar", (2, -0.1)), ("_cdf_scalar", (0, 1.5)), From 2c000e3670cb7fd699d23d1a546f15c4b2bc9165 Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 00:11:22 -0500 Subject: [PATCH 02/13] Added pycaches to gitignore Added missing packages to __init__.py and added a test suite to verify this --- .gitignore | 5 ++++- python/fastdist/__init__.py | 4 ++-- python/tests/test_package_api.py | 34 ++++++++++++++++++++++++++++++++ 3 files changed, 40 insertions(+), 3 deletions(-) create mode 100644 python/tests/test_package_api.py diff --git a/.gitignore b/.gitignore index 0159c5e..1cf4af0 100644 --- a/.gitignore +++ b/.gitignore @@ -11,4 +11,7 @@ /.idea/ /.venv3.14/ /.venv3.13/ -/.venv3.12/ \ No newline at end of file +/.venv3.12/ +__pycache__/ +*.py[cod] +.pytest_cache/ \ No newline at end of file diff --git a/python/fastdist/__init__.py b/python/fastdist/__init__.py index c62c58c..a07d63c 100644 --- a/python/fastdist/__init__.py +++ b/python/fastdist/__init__.py @@ -6,7 +6,7 @@ Poisson, Uniform, Utils ) -__all__ = ["Bernoulli", "Beta", "Binomial", +__all__ = ["Bernoulli", "Beta", "Binomial", "ChiSquare", "DiscreteUniform", "Exponential", "Gamma", - "Geometric", "Normal", "NegativeBinomial", + "Geometric", "NegativeBinomial", "Normal", "Poisson", "Uniform", "Utils"] diff --git a/python/tests/test_package_api.py b/python/tests/test_package_api.py new file mode 100644 index 0000000..fe2990c --- /dev/null +++ b/python/tests/test_package_api.py @@ -0,0 +1,34 @@ +import importlib + +import pytest + +import fastdist +import fastdist.distributions as distributions + +EXPECTED = [ + "Bernoulli", "Beta", "Binomial", "ChiSquare", + "DiscreteUniform", "Exponential", "Gamma", + "Geometric", "NegativeBinomial", "Normal", + "Poisson", "Uniform", "Utils", +] + + +def test_top_level_all_matches_expected(): + assert sorted(fastdist.__all__) == sorted(EXPECTED) + + +def test_distributions_all_matches_top_level(): + assert sorted(distributions.__all__) == sorted(fastdist.__all__) + + +@pytest.mark.parametrize("name", EXPECTED) +def test_every_exported_name_is_importable(name): + assert hasattr(fastdist, name), f"{name} missing from fastdist" + assert hasattr(distributions, name), f"{name} missing from fastdist.distributions" + + +@pytest.mark.parametrize("name", EXPECTED) +def test_every_class_has_a_module(name): + # every exported class must come from a real submodule, not a stale alias + cls = getattr(fastdist, name) + assert importlib.import_module(cls.__module__) is not None \ No newline at end of file From bf4eac950b5e100de53aa278c34980dc495427c9 Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 00:12:50 -0500 Subject: [PATCH 03/13] Updated github actions files to run on PRs with both develop and master --- .github/workflows/clang-format.yml | 8 +++----- .github/workflows/python-distro.yml | 2 +- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index 7406052..2703721 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -2,12 +2,10 @@ name: Check Clang Format on: push: - branches: - - master - - develop + branches: [ master, develop ] pull_request: - branches: - - master + branches: [ master, develop ] + jobs: clang-format: diff --git a/.github/workflows/python-distro.yml b/.github/workflows/python-distro.yml index 0fa757f..0800b15 100644 --- a/.github/workflows/python-distro.yml +++ b/.github/workflows/python-distro.yml @@ -4,7 +4,7 @@ on: push: branches: [ master, develop ] pull_request: - branches: [ master ] + branches: [ master, develop ] jobs: build: From c766c4bf21571ffb199c89c7b0e5e37de32aaf5b Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 00:17:49 -0500 Subject: [PATCH 04/13] Updated gpu_capacity check to be more user-friendly if a gpu driver is not present, added testing to cover config cases --- python/fastdist/config.py | 94 +++++++++++++++++++++++++++++------ python/tests/test_config.py | 99 +++++++++++++++++++++++++++++++++++++ requirements.txt | 3 +- 3 files changed, 180 insertions(+), 16 deletions(-) create mode 100644 python/tests/test_config.py diff --git a/python/fastdist/config.py b/python/fastdist/config.py index f600f87..821ced6 100644 --- a/python/fastdist/config.py +++ b/python/fastdist/config.py @@ -4,9 +4,13 @@ import time import numpy as np -import pynvml from pathlib import Path +try: + import pynvml +except ImportError: + pynvml = None + # ---------------------------------------------------------------------------------------------------------------------- # Constants and Defaults # ---------------------------------------------------------------------------------------------------------------------- @@ -560,24 +564,84 @@ def set_cuda_threshold(func_name: str, value: int) -> None: CUDA_THRESHOLDS[fd_class][fd_func] = value -def validate_gpu_capacity(array_size: int, dtype_item_size: int): - try: - pynvml.nvmlInit() - handle = pynvml.nvmlDeviceGetHandleByIndex(0) - info = pynvml.nvmlDeviceGetMemoryInfo(handle) +_NVML_STATE = None # None = not yet attempted, True = initialized, False = unavailable + +def _nvml_ready() -> bool: + """ + Initialize NVML once per process. + + Returns + ------- + bool + True if NVML is initialized and usable, False if it is unavailable. + """ + + global _NVML_STATE + + if _NVML_STATE is None: + if pynvml is None: + _NVML_STATE = False + else: + try: + pynvml.nvmlInit() + _NVML_STATE = True + atexit.register(_nvml_teardown) + except pynvml.NVMLError: + _NVML_STATE = False - required = array_size * dtype_item_size * 2 # Factor of 2 for input and output arrays + return _NVML_STATE - if required > info.free: - raise MemoryError( - f"GPU Memory Overflow: Required {required / 1e6:.2f}MB, " - f"but only {info.free / 1e6:.2f}MB is free." - ) +def _nvml_teardown() -> None: + """Shut down NVML at interpreter exit. Never raises.""" + + global _NVML_STATE + + if _NVML_STATE: + try: + pynvml.nvmlShutdown() + except pynvml.NVMLError: + pass + _NVML_STATE = False + +def validate_gpu_capacity(array_size: int, dtype_item_size: int, device_index: int = 0) -> None: + """ + Verify the GPU has enough free memory for the requested operation. + + Parameters + ---------- + array_size : int + Number of elements in the input array. + dtype_item_size : int + Size in bytes of a single element. + device_index : int, optional + Index of the GPU to query (default 0). + + Raises + ------ + MemoryError + If the operation would require more memory than is currently free. + + Notes + ----- + Silently returns if NVML is unavailable; the check is advisory, not a hard gate. + """ + + if not _nvml_ready(): + return + + try: + handle = pynvml.nvmlDeviceGetHandleByIndex(device_index) + info = pynvml.nvmlDeviceGetMemoryInfo(handle) except pynvml.NVMLError: - pass - finally: - pynvml.nvmlShutdown() + return + + required = array_size * dtype_item_size * 2 # Factor of 2 for input and output arrays + if required > info.free: + raise MemoryError( + f"GPU Memory Overflow: Required {required / 1e6:.2f}MB, " + f"but only {info.free / 1e6:.2f}MB is free." + ) # ---------------------------------------------------------------------------------------------------------------------- # Module Initialization diff --git a/python/tests/test_config.py b/python/tests/test_config.py new file mode 100644 index 0000000..b918119 --- /dev/null +++ b/python/tests/test_config.py @@ -0,0 +1,99 @@ +import pytest + +from fastdist import config + + +class FakeNVMLError(Exception): + pass + + +class _Mem: + def __init__(self, free): + self.free = free + + +def _install(monkeypatch, fake): + monkeypatch.setattr(config, "pynvml", fake) + monkeypatch.setattr(config, "_NVML_STATE", None) + + +def test_no_driver_returns_silently(monkeypatch): + class FakeNVML: + NVMLError = FakeNVMLError + + @staticmethod + def nvmlInit(): + raise FakeNVMLError("library not found") + + _install(monkeypatch, FakeNVML) + config.validate_gpu_capacity(1_000_000, 8) # must not raise + + +def test_pynvml_not_installed_returns_silently(monkeypatch): + _install(monkeypatch, None) + config.validate_gpu_capacity(1_000_000, 8) # must not raise + + +def test_raises_when_insufficient_memory(monkeypatch): + class FakeNVML: + NVMLError = FakeNVMLError + + @staticmethod + def nvmlInit(): + return None + + @staticmethod + def nvmlDeviceGetHandleByIndex(index): + return object() + + @staticmethod + def nvmlDeviceGetMemoryInfo(handle): + return _Mem(free=1_000) + + _install(monkeypatch, FakeNVML) + with pytest.raises(MemoryError, match="GPU Memory Overflow"): + config.validate_gpu_capacity(1_000_000, 8) + + +def test_passes_when_memory_sufficient(monkeypatch): + class FakeNVML: + NVMLError = FakeNVMLError + + @staticmethod + def nvmlInit(): + return None + + @staticmethod + def nvmlDeviceGetHandleByIndex(index): + return object() + + @staticmethod + def nvmlDeviceGetMemoryInfo(handle): + return _Mem(free=10 ** 12) + + _install(monkeypatch, FakeNVML) + config.validate_gpu_capacity(1_000, 8) # must not raise + + +def test_init_happens_only_once(monkeypatch): + calls = {"n": 0} + + class FakeNVML: + NVMLError = FakeNVMLError + + @staticmethod + def nvmlInit(): + calls["n"] += 1 + + @staticmethod + def nvmlDeviceGetHandleByIndex(index): + return object() + + @staticmethod + def nvmlDeviceGetMemoryInfo(handle): + return _Mem(free=10 ** 12) + + _install(monkeypatch, FakeNVML) + for _ in range(5): + config.validate_gpu_capacity(1_000, 8) + assert calls["n"] == 1 \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 57e693c..48a7ca9 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,4 +5,5 @@ wheel pre-commit pytest numpy -pynvml \ No newline at end of file +pynvml +nvidia-ml-py \ No newline at end of file From 6f8cd97ad510c3e691a79f9c47b4c70e5cebd4e3 Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 00:37:36 -0500 Subject: [PATCH 05/13] Fixed building errors and cmake pybind compiling bugs --- CMakeLists.txt | 27 ++++++++++++----------- README.md | 2 +- python/fastdist/config.py | 1 + python/fastdist/distributions/__init__.py | 5 +++++ setup.py | 2 +- 5 files changed, 22 insertions(+), 15 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5ea00c3..89683e6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,26 +19,27 @@ endif () # ----------------------------- # python # ----------------------------- -find_package(Python 3.14 REQUIRED COMPONENTS Interpreter Development) +find_package(Python 3.12 REQUIRED COMPONENTS Interpreter Development.Module) +message(STATUS "Building against Python ${Python_VERSION} at ${Python_EXECUTABLE}") # ----------------------------- # pybind11 # ----------------------------- -# Setting up python-pybind11 connection set(PYBIND11_FINDPYTHON ON) -set(PYBIND11_PYTHON_EXECUTABLE ${Python_EXECUTABLE}) -set(PYBIND11_PYTHON_VERSION ${Python_VERSION_STRING}) -# Initializing pybind11 library -include(FetchContent) - -FetchContent_Declare( - pybind11 - GIT_REPOSITORY https://github.com/pybind/pybind11.git - GIT_TAG v2.12.1 -) +# Use the pybind11 installed in the active environment (already declared in +# pyproject.toml build-system.requires) so headers always match the interpreter. +if (NOT pybind11_DIR) + execute_process( + COMMAND "${Python_EXECUTABLE}" -m pybind11 --cmakedir + OUTPUT_VARIABLE pybind11_DIR + OUTPUT_STRIP_TRAILING_WHITESPACE + COMMAND_ERROR_IS_FATAL ANY + ) +endif () -FetchContent_MakeAvailable(pybind11) +find_package(pybind11 CONFIG REQUIRED) +message(STATUS "Using pybind11 ${pybind11_VERSION} from ${pybind11_DIR}") pybind11_add_module(_fastdist python/bindings/bindings.cpp diff --git a/README.md b/README.md index 25c9cb8..cec682a 100644 --- a/README.md +++ b/README.md @@ -33,7 +33,7 @@ git submodule update --init --recursive 2. From the **project root**, build the Python wheel: ```bash -python3 python/setup.py bdist_wheel +python3 setup.py bdist_wheel ``` **Important:** diff --git a/python/fastdist/config.py b/python/fastdist/config.py index 821ced6..c358244 100644 --- a/python/fastdist/config.py +++ b/python/fastdist/config.py @@ -1,3 +1,4 @@ +import atexit import copy import json import platform diff --git a/python/fastdist/distributions/__init__.py b/python/fastdist/distributions/__init__.py index 30fffaf..bee997d 100644 --- a/python/fastdist/distributions/__init__.py +++ b/python/fastdist/distributions/__init__.py @@ -12,3 +12,8 @@ from .poisson import Poisson from .uniform import Uniform from .utils import Utils + +__all__ = ["Bernoulli", "Beta", "Binomial", "ChiSquare", + "DiscreteUniform", "Exponential", "Gamma", + "Geometric", "NegativeBinomial", "Normal", + "Poisson", "Uniform", "Utils"] \ No newline at end of file diff --git a/setup.py b/setup.py index 1a28af6..69b2074 100644 --- a/setup.py +++ b/setup.py @@ -59,7 +59,7 @@ def build_extension(self, ext: CMakeExtension) -> None: cmake_args = [ f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY={extdir}{os.sep}", - f"-DPYTHON_EXECUTABLE={sys.executable}", + f"-DPython_EXECUTABLE={sys.executable}", f"-DCMAKE_BUILD_TYPE={cfg}", ] From 0d67f7f418af9ee992be7be9d770fea56aa3ff4b Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 01:20:22 -0500 Subject: [PATCH 06/13] Updated versioning to have a unified setting (in CMakeLists.txt, near the top). This should cascade from the top to the rest of the doc --- CMakeLists.txt | 13 ++++++++++++- README.md | 6 +++++- include/fastdist/version.h | 20 -------------------- include/fastdist/version.h.in | 25 +++++++++++++++++++++++++ python/bindings/bindings.cpp | 3 ++- python/fastdist/__init__.py | 12 ++++++++++-- python/tests/test_version.py | 27 +++++++++++++++++++++++++++ requirements.txt | 1 - setup.py | 19 ++++++++++++++++++- 9 files changed, 99 insertions(+), 27 deletions(-) delete mode 100644 include/fastdist/version.h create mode 100644 include/fastdist/version.h.in create mode 100644 python/tests/test_version.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 89683e6..8e839be 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,9 +1,19 @@ cmake_minimum_required(VERSION 3.20) -project(fastdist LANGUAGES CXX) +project(fastdist VERSION 0.1.0 LANGUAGES CXX) set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED ON) +# ----------------------------- +# Generated version header +# ----------------------------- +configure_file( + ${PROJECT_SOURCE_DIR}/include/fastdist/version.h.in + ${PROJECT_BINARY_DIR}/generated/fastdist/version.h + @ONLY +) +message(STATUS "fastdist version ${PROJECT_VERSION}") + # ----------------------------- # CUDA Backend # ----------------------------- @@ -195,6 +205,7 @@ target_include_directories(fastdist_core PUBLIC $ $ + $ ${CUDAToolkit_INCLUDE_DIRS} $ ) diff --git a/README.md b/README.md index cec682a..980d6b4 100644 --- a/README.md +++ b/README.md @@ -43,7 +43,7 @@ python3 setup.py bdist_wheel 3. Install the generated wheel: ```bash -pip install .\dist\fastdist-0.0.1-cpXXX-cpXXX-win_amd64.whl --force-reinstall +pip install .\dist\fastdist--cpXXX-cpXXX-win_amd64.whl --force-reinstall ``` --- @@ -92,6 +92,10 @@ The normal distribution is a good reference for creating new distributions. \ --- +## Updating the Version +Only update the version from the CMakeLists.txt at the line: +`project(fastdist VERSION x.y.z LANGUAGES CXX)` + ## Building Wheels for Multiple Python Versions (3.12–3.14) To generate wheels for all currently supported Python versions: diff --git a/include/fastdist/version.h b/include/fastdist/version.h deleted file mode 100644 index 48417ff..0000000 --- a/include/fastdist/version.h +++ /dev/null @@ -1,20 +0,0 @@ -// Version header for the FastDist library -#pragma once - -// Version denoted as major.minor.patch -#define FASTDIST_VERSION_MAJOR 0 -#define FASTDIST_VERSION_MINOR 1 -#define FASTDIST_VERSION_PATCH 0 - -#ifdef __cplusplus -extern "C" { -#endif - -// Returns the major, minor, and patch version numbers respectively -int fd_version_major(void); -int fd_version_minor(void); -int fd_version_patch(void); - -#ifdef __cplusplus -} -#endif diff --git a/include/fastdist/version.h.in b/include/fastdist/version.h.in new file mode 100644 index 0000000..1d0dd8a --- /dev/null +++ b/include/fastdist/version.h.in @@ -0,0 +1,25 @@ +// Version header for the FastDist library +// +// GENERATED FILE — do not edit the copy in your build directory. +// Edit include/fastdist/version.h.in instead, and change the version +// number only in the project() call in CMakeLists.txt. +#pragma once + +// Version denoted as major.minor.patch +#define FASTDIST_VERSION_MAJOR @PROJECT_VERSION_MAJOR@ +#define FASTDIST_VERSION_MINOR @PROJECT_VERSION_MINOR@ +#define FASTDIST_VERSION_PATCH @PROJECT_VERSION_PATCH@ +#define FASTDIST_VERSION_STRING "@PROJECT_VERSION@" + +#ifdef __cplusplus +extern "C" { +#endif + +// Returns the major, minor, and patch version numbers respectively +int fd_version_major(void); +int fd_version_minor(void); +int fd_version_patch(void); + +#ifdef __cplusplus +} +#endif \ No newline at end of file diff --git a/python/bindings/bindings.cpp b/python/bindings/bindings.cpp index 9455258..5a2890e 100644 --- a/python/bindings/bindings.cpp +++ b/python/bindings/bindings.cpp @@ -1,6 +1,7 @@ // CPP file to link all other bindings #include #include +#include namespace py = pybind11; @@ -36,5 +37,5 @@ PYBIND11_MODULE(_fastdist, m) { bind_chi_square(m); bind_utils(m); - m.attr("__version__") = "0.0.1"; + m.attr("__version__") = FASTDIST_VERSION_STRING; } diff --git a/python/fastdist/__init__.py b/python/fastdist/__init__.py index a07d63c..b943f97 100644 --- a/python/fastdist/__init__.py +++ b/python/fastdist/__init__.py @@ -1,4 +1,6 @@ # python/fastdist/__init__.py +from importlib.metadata import PackageNotFoundError, version as _pkg_version + from . import _fastdist from .distributions import ( Bernoulli, Beta, Binomial, ChiSquare, DiscreteUniform, @@ -6,7 +8,13 @@ Poisson, Uniform, Utils ) -__all__ = ["Bernoulli", "Beta", "Binomial", "ChiSquare", +try: + __version__ = _pkg_version("fastdist") +except PackageNotFoundError: # running from an uninstalled source checkout + __version__ = "0.0.0+unknown" + +__all__ = ["__version__", + "Bernoulli", "Beta", "Binomial", "ChiSquare", "DiscreteUniform", "Exponential", "Gamma", "Geometric", "NegativeBinomial", "Normal", - "Poisson", "Uniform", "Utils"] + "Poisson", "Uniform", "Utils"] \ No newline at end of file diff --git a/python/tests/test_version.py b/python/tests/test_version.py new file mode 100644 index 0000000..f899aa1 --- /dev/null +++ b/python/tests/test_version.py @@ -0,0 +1,27 @@ +import re +from pathlib import Path + +import fastdist +import fastdist._fastdist as core + +CMAKELISTS = Path(__file__).resolve().parents[2] / "CMakeLists.txt" + + +def _cmake_version() -> str: + """The single source of truth: the project() call in CMakeLists.txt.""" + match = re.search( + r"project\s*\(\s*fastdist\s+VERSION\s+(\d+\.\d+\.\d+)", + CMAKELISTS.read_text(encoding="utf-8"), + ) + assert match is not None, "no VERSION found in the project() call" + return match.group(1) + + +def test_compiled_module_matches_cmake(): + """The C++ constant must match what CMake declared.""" + assert core.__version__ == _cmake_version() + + +def test_python_package_matches_compiled_module(): + """The wheel metadata must match the C++ constant.""" + assert fastdist.__version__ == core.__version__ \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 48a7ca9..d648bed 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,5 +5,4 @@ wheel pre-commit pytest numpy -pynvml nvidia-ml-py \ No newline at end of file diff --git a/setup.py b/setup.py index 69b2074..5d49229 100644 --- a/setup.py +++ b/setup.py @@ -147,10 +147,27 @@ def build_extension(self, ext: CMakeExtension) -> None: ["cmake", "--build", ".", *build_args], cwd=build_temp, check=True ) +def _read_version() -> str: + """ + Read the project version from the CMakeLists.txt project() call. + + CMakeLists.txt is the single source of truth for the version; the C++ + header and the Python package metadata are both derived from it. + """ + cmakelists = (Path(__file__).parent / "CMakeLists.txt").read_text(encoding="utf-8") + match = re.search( + r"project\s*\(\s*fastdist\s+VERSION\s+(\d+\.\d+\.\d+)", cmakelists + ) + if match is None: + raise RuntimeError( + "Could not parse the version from the project() call in CMakeLists.txt" + ) + return match.group(1) + setup( name="fastdist", # pip install fastdist - version="0.0.1", + version=_read_version(), author="Emanuel McGrail and Zachery Pipes", author_email="geometrydashgodwave@gmail.com", description="Manny!", From 653282e02635dfb79aa99d76431806c721b16fb4 Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 01:26:53 -0500 Subject: [PATCH 07/13] Updated missing setup.py information and faulty api tests --- python/tests/test_package_api.py | 6 ++++-- setup.py | 24 ++++++++++++++++++++++-- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/python/tests/test_package_api.py b/python/tests/test_package_api.py index fe2990c..e84974d 100644 --- a/python/tests/test_package_api.py +++ b/python/tests/test_package_api.py @@ -14,11 +14,13 @@ def test_top_level_all_matches_expected(): - assert sorted(fastdist.__all__) == sorted(EXPECTED) + # the top level exports every distribution class plus __version__ + assert sorted(fastdist.__all__) == sorted(EXPECTED + ["__version__"]) def test_distributions_all_matches_top_level(): - assert sorted(distributions.__all__) == sorted(fastdist.__all__) + # the distributions subpackage exports the classes only + assert sorted(distributions.__all__) == sorted(EXPECTED) @pytest.mark.parametrize("name", EXPECTED) diff --git a/setup.py b/setup.py index 5d49229..02e5b1e 100644 --- a/setup.py +++ b/setup.py @@ -170,11 +170,31 @@ def _read_version() -> str: version=_read_version(), author="Emanuel McGrail and Zachery Pipes", author_email="geometrydashgodwave@gmail.com", - description="Manny!", + description="High-performance probability distributions and statistical functions with C++ and CUDA backends", + long_description=(Path(__file__).parent / "README.md").read_text(encoding="utf-8"), + long_description_content_type="text/markdown", + url="https://github.com/ghosteau/fastdist", + license="Apache-2.0", + license_files=["LICENSE"], package_dir={"": "python"}, packages=["fastdist", "fastdist.distributions"], ext_modules=[CMakeExtension("fastdist._fastdist")], # (.pyd file) cmdclass={"build_ext": CMakeBuild}, zip_safe=False, - python_requires=">=3.7" + python_requires=">=3.7", + install_requires=[ + "numpy>=1.21", + ], + extras_require={ + "gpu": ["nvidia-ml-py"], + }, + classifiers=[ + "Development Status :: 3 - Alpha", + "Intended Audience :: Science/Research", + "Programming Language :: C++", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", + "Topic :: Scientific/Engineering :: Mathematics", + ], ) From ac71924d9a59821e7c44292c4ae00cc4452c218b Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 10:35:14 -0500 Subject: [PATCH 08/13] Revamped python tests to not use mock and instead test against the actual distribution Identified known bugs and flagged them for fixing further down the line --- python/tests/conftest.py | 156 +++++++++ python/tests/test_beta.py | 319 +++++++++++++----- python/tests/test_binomial.py | 364 ++++++++++++++------ python/tests/test_chi_square.py | 301 +++++++++++------ python/tests/test_discrete_uniform.py | 329 ++++++++++++------- python/tests/test_exponential.py | 326 +++++++++++++----- python/tests/test_gamma.py | 402 +++++++++++++++-------- python/tests/test_geometric.py | 296 ++++++++++------- python/tests/test_negative_binomial.py | 397 +++++++++++++--------- python/tests/test_poisson.py | 337 +++++++++++++------ python/tests/test_uniform.py | 401 +++++++++++++++++------ python/tests/test_utils.py | 437 +++++++++++++++++++------ 12 files changed, 2852 insertions(+), 1213 deletions(-) create mode 100644 python/tests/conftest.py diff --git a/python/tests/conftest.py b/python/tests/conftest.py new file mode 100644 index 0000000..488fb73 --- /dev/null +++ b/python/tests/conftest.py @@ -0,0 +1,156 @@ +""" +Shared pytest configuration and reference implementations for the fastdist suite. + +The suite exercises the *real* compiled extension. Every numeric assertion is +made against an independently derived closed form or a reference implementation +written here, rather than against a value captured from a previous run, so a +regression in the C++ or CUDA backend surfaces as a test failure instead of +passing silently. +""" + +import math +import sys +from pathlib import Path + +# Fall back to the source tree when the package has not been pip-installed. +try: + import fastdist # noqa: F401 +except ImportError: # pragma: no cover - exercised only in uninstalled checkouts + sys.path.insert(0, str(Path(__file__).resolve().parents[1])) + + +# ---------------------------------------------------------------------------- +# Numeric tolerances +# ---------------------------------------------------------------------------- +# Closed-form algebra (means, variances, MGFs) is accurate to machine epsilon. +EXACT = {"rel": 1e-12, "abs": 1e-15} + +# Iterative routines (regularized incomplete gamma/beta) converge to EPS=1e-12 +# as configured in include/fastdist/config.h; leave headroom above that. +ITERATIVE = {"rel": 1e-10, "abs": 1e-12} + + +# ---------------------------------------------------------------------------- +# Reference implementations +# ---------------------------------------------------------------------------- + +def regularized_lower_gamma(a: float, x: float) -> float: + """ + P(a, x), the regularized lower incomplete gamma function. + + Reference implementation using the standard series expansion for x < a+1 + and the Lentz continued fraction for x >= a+1 (Numerical Recipes 6.2). + Used to validate the backend's Gamma and Chi-square CDFs. + """ + if a <= 0.0 or x < 0.0: + raise ValueError("a must be positive and x non-negative") + if x == 0.0: + return 0.0 + + gln = math.lgamma(a) + + if x < a + 1.0: + # Series representation + ap = a + total = 1.0 / a + delta = total + for _ in range(1000): + ap += 1.0 + delta *= x / ap + total += delta + if abs(delta) < abs(total) * 1e-16: + break + return total * math.exp(-x + a * math.log(x) - gln) + + # Continued fraction representation (modified Lentz) + tiny = 1e-300 + b = x + 1.0 - a + c = 1.0 / tiny + d = 1.0 / b + h = d + for i in range(1, 1000): + an = -i * (i - a) + b += 2.0 + d = an * d + b + if abs(d) < tiny: + d = tiny + c = b + an / c + if abs(c) < tiny: + c = tiny + d = 1.0 / d + delta = d * c + h *= delta + if abs(delta - 1.0) < 1e-16: + break + return 1.0 - math.exp(-x + a * math.log(x) - gln) * h + + +def regularized_incomplete_beta(a: float, b: float, x: float) -> float: + """ + I_x(a, b), the regularized incomplete beta function. + + Reference implementation using the continued fraction of Numerical Recipes + 6.4 with the standard symmetry transformation. Used to validate the + backend's Beta CDF. + """ + if x <= 0.0: + return 0.0 + if x >= 1.0: + return 1.0 + + front = math.exp( + math.lgamma(a + b) + - math.lgamma(a) + - math.lgamma(b) + + a * math.log(x) + + b * math.log(1.0 - x) + ) + + # Strict comparison: at exactly the swap boundary with a == b, a non-strict + # test would swap to the identical point and recurse forever. + if x > (a + 1.0) / (a + b + 2.0): + return 1.0 - regularized_incomplete_beta(b, a, 1.0 - x) + + tiny = 1e-300 + c = 1.0 + d = 1.0 - (a + b) * x / (a + 1.0) + if abs(d) < tiny: + d = tiny + d = 1.0 / d + h = d + + for m in range(1, 1000): + m2 = 2 * m + + numerator = m * (b - m) * x / ((a + m2 - 1.0) * (a + m2)) + d = 1.0 + numerator * d + if abs(d) < tiny: + d = tiny + c = 1.0 + numerator / c + if abs(c) < tiny: + c = tiny + d = 1.0 / d + h *= d * c + + numerator = -(a + m) * (a + b + m) * x / ((a + m2) * (a + m2 + 1.0)) + d = 1.0 + numerator * d + if abs(d) < tiny: + d = tiny + c = 1.0 + numerator / c + if abs(c) < tiny: + c = tiny + d = 1.0 / d + delta = d * c + h *= delta + + if abs(delta - 1.0) < 1e-15: + break + + return front * h / a + + +def pytest_configure(config): + config.addinivalue_line( + "markers", + "known_bug: marks a test that documents a confirmed defect in the backend", + ) diff --git a/python/tests/test_beta.py b/python/tests/test_beta.py index 43b8f9f..f5e7e38 100644 --- a/python/tests/test_beta.py +++ b/python/tests/test_beta.py @@ -1,102 +1,228 @@ -# tests/python_modules/test_beta.py -from unittest.mock import patch +"""Numeric tests for the Beta distribution against closed-form references.""" + +import math import pytest -import fastdist.distributions.beta as beta_module +from conftest import EXACT, ITERATIVE from fastdist.distributions.beta import Beta -@pytest.fixture -def mock_core(): - """Patch the internal C++ core for the duration of each test.""" - with patch.object(beta_module, "_core", create=True) as mock: - yield mock +# --------------------------------------------------------------------------- +# Closed-form references +# --------------------------------------------------------------------------- + +def beta_fn(a: float, b: float) -> float: + """B(a, b) = Gamma(a) Gamma(b) / Gamma(a + b)""" + return math.gamma(a) * math.gamma(b) / math.gamma(a + b) + + +def beta_pdf(x: float, a: float, b: float) -> float: + """f(x; a, b) = x^(a-1) (1-x)^(b-1) / B(a, b)""" + return x ** (a - 1) * (1 - x) ** (b - 1) / beta_fn(a, b) + + +def beta_mean(a: float, b: float) -> float: + return a / (a + b) + + +def beta_variance(a: float, b: float) -> float: + return a * b / ((a + b) ** 2 * (a + b + 1)) + + +PARAMS = [(2.0, 3.0), (1.0, 1.0), (0.5, 0.5), (5.0, 2.0), (3.0, 3.0)] # --------------------------------------------------------------------------- -# Constructor & representation +# Constructor, properties, representation # --------------------------------------------------------------------------- def test_init_valid_parameters(): - b = Beta(alpha=2.0, beta=3.0) - assert b.alpha == 2.0 - assert b.beta == 3.0 - - -@pytest.mark.parametrize("alpha,beta", [ - (-0.1, 1.0), - (0, 1.0), - (-1, 1.0), - (1.0, -0.1), - (1.0, 0), - (1.0, -1), - (-1, -1), -]) -def test_init_invalid_parameters_raises(alpha, beta): + dist = Beta(alpha=2.0, beta=3.0) + assert dist.alpha == 2.0 + assert dist.beta == 3.0 + + +@pytest.mark.parametrize("alpha, beta", [(0, 1.0), (-1, 1.0), (1.0, 0), (1.0, -1)]) +def test_init_rejects_non_positive_parameters(alpha, beta): with pytest.raises(ValueError, match="must be positive"): Beta(alpha=alpha, beta=beta) +@pytest.mark.parametrize("bad", ["2.0", None, object()]) +def test_init_rejects_non_real_parameters(bad): + with pytest.raises((TypeError, ValueError)): + Beta(alpha=bad, beta=1.0) + + def test_repr(): - b = Beta(alpha=2.0, beta=3.0) - assert repr(b) == "Beta(alpha=2.0, beta=3.0)" + assert repr(Beta(alpha=2.0, beta=3.0)) == "Beta(alpha=2.0, beta=3.0)" + + +def test_alpha_setter_updates_and_validates(): + dist = Beta(alpha=2.0, beta=3.0) + dist.alpha = 4.0 + assert dist.alpha == 4.0 + with pytest.raises(ValueError, match="alpha must be positive"): + dist.alpha = -1.0 + + +# KNOWN BUG: the beta setter calls _validate_params(beta=value), leaving alpha +# as None. Because the `if alpha <= 0` check sits outside the `if alpha is not +# None` guard (beta.py line 47), the comparison None <= 0 raises TypeError for +# *every* assignment, valid or not. The beta property is unusable. +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="beta setter raises TypeError for any value") +def test_beta_setter_updates_value(): + dist = Beta(alpha=2.0, beta=3.0) + dist.beta = 5.0 + assert dist.beta == 5.0 + + +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="beta setter raises TypeError before validating") +def test_beta_setter_rejects_non_positive(): + dist = Beta(alpha=2.0, beta=3.0) + with pytest.raises(ValueError, match="beta must be positive"): + dist.beta = 0 # --------------------------------------------------------------------------- -# Validation behavior +# PDF # --------------------------------------------------------------------------- -@pytest.mark.parametrize("alpha,beta", [ - (0.1, 0.1), - (1.0, 1.0), - (2.5, 3.5), - (100, 200), -]) -def test_validate_params_accepts_valid_values(alpha, beta): - Beta._validate_params(alpha=alpha, beta=beta) +@pytest.mark.parametrize("alpha, beta", PARAMS) +@pytest.mark.parametrize("x", [0.1, 0.25, 0.5, 0.75, 0.9]) +def test_pdf_matches_closed_form(x, alpha, beta): + assert Beta(alpha, beta).pdf_scalar(x) == pytest.approx( + beta_pdf(x, alpha, beta), **EXACT + ) -@pytest.mark.parametrize("alpha,beta", [ - (0, 1.0), - (-1, 1.0), - (1.0, 0), - (1.0, -1), -]) -def test_validate_params_rejects_invalid_values(alpha, beta): - with pytest.raises(ValueError, match="must be positive"): - Beta._validate_params(alpha=alpha, beta=beta) +@pytest.mark.parametrize("x", [0.1, 0.5, 0.9]) +def test_pdf_of_uniform_special_case_is_one(x): + """Beta(1, 1) is the uniform distribution on [0, 1].""" + assert Beta(1.0, 1.0).pdf_scalar(x) == pytest.approx(1.0, **EXACT) + + +@pytest.mark.parametrize("alpha, beta", PARAMS) +@pytest.mark.parametrize("x", [0.05, 0.4, 0.95]) +def test_pdf_is_non_negative(x, alpha, beta): + assert Beta(alpha, beta).pdf_scalar(x) >= 0.0 + + +@pytest.mark.parametrize("alpha, beta", [(2.0, 3.0), (5.0, 2.0)]) +def test_pdf_is_symmetric_under_parameter_swap(alpha, beta): + """f(x; a, b) == f(1-x; b, a)""" + assert Beta(alpha, beta).pdf_scalar(0.3) == pytest.approx( + Beta(beta, alpha).pdf_scalar(0.7), **EXACT + ) # --------------------------------------------------------------------------- -# Instance method delegation +# Moments # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, value", [ - ("pdf_scalar", "beta_pdf_scalar", 0.5), - ("cdf_scalar", "beta_cdf_scalar", 0.3), - ("mean", "beta_mean", 0.4), - ("variance", "beta_variance", 0.048), - ("stddev", "beta_stddev", 0.219), - ("sample", "beta_sample", 0.42), -]) -def test_instance_methods_delegate_to_core(mock_core, method_name, core_method_name, value): - getattr(mock_core, core_method_name).return_value = value - b = Beta(alpha=2.0, beta=3.0) - method = getattr(b, method_name) - if "scalar" in method_name: - result = method(0.5) - getattr(mock_core, core_method_name).assert_called_once_with(0.5, 2.0, 3.0) - else: - result = method() - getattr(mock_core, core_method_name).assert_called_once_with(2.0, 3.0) - assert result == value +@pytest.mark.parametrize("alpha, beta", PARAMS) +def test_moments_match_closed_forms(alpha, beta): + dist = Beta(alpha, beta) + assert dist.mean() == pytest.approx(beta_mean(alpha, beta), **EXACT) + assert dist.variance() == pytest.approx(beta_variance(alpha, beta), **EXACT) + assert dist.stddev() == pytest.approx( + math.sqrt(beta_variance(alpha, beta)), **EXACT + ) + + +@pytest.mark.parametrize("alpha, beta", PARAMS) +def test_mean_lies_inside_the_support(alpha, beta): + assert 0.0 < Beta(alpha, beta).mean() < 1.0 + + +# --------------------------------------------------------------------------- +# CDF +# +# KNOWN BUG: the regularized incomplete beta in src/math/beta.cpp is incorrect. +# It disagrees with numerical integration at every tested point, returns 0.6534 +# for Beta(1, 1) at x=0.5 where the exact answer is 0.5, and returns a negative +# value (-0.5958) for Beta(0.5, 0.5) at x=0.5, which is impossible for a CDF. +# These tests assert the correct behaviour and are marked strict-xfail so they +# turn into failures the moment the backend is fixed and the marker goes stale. +# --------------------------------------------------------------------------- + +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="beta_cdf_scalar is numerically incorrect") +@pytest.mark.parametrize("x", [0.1, 0.3, 0.5, 0.75, 0.9]) +def test_cdf_of_uniform_special_case_is_identity(x): + """Beta(1, 1) is uniform on [0, 1], so its CDF is F(x) = x exactly.""" + assert Beta(1.0, 1.0).cdf_scalar(x) == pytest.approx(x, **ITERATIVE) + + +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="beta_cdf_scalar is numerically incorrect") +def test_cdf_matches_analytic_integral(): + """For Beta(2, 3), F(0.5) = 0.6875 by direct integration of 12x(1-x)^2.""" + assert Beta(2.0, 3.0).cdf_scalar(0.5) == pytest.approx(0.6875, **ITERATIVE) + + +# The CDF is wrong everywhere, but it still happens to be bounded and monotonic +# for most parameters. It breaks both properties only for alpha = beta = 0.5, +# where it returns values as low as -0.4273. Only that case is xfailed, so the +# structural guarantees stay enforced for every other parameter pair. +CDF_PROPERTY_PARAMS = [ + (2.0, 3.0), + (1.0, 1.0), + pytest.param( + 0.5, 0.5, + marks=[ + pytest.mark.known_bug, + pytest.mark.xfail( + strict=True, + reason="beta_cdf_scalar returns negative values for alpha=beta=0.5", + ), + ], + ), + (5.0, 2.0), + (3.0, 3.0), +] + + +@pytest.mark.parametrize("alpha, beta", CDF_PROPERTY_PARAMS) +def test_cdf_is_bounded(alpha, beta): + dist = Beta(alpha, beta) + for x in (0.05, 0.25, 0.5, 0.75, 0.95): + assert 0.0 <= dist.cdf_scalar(x) <= 1.0 + + +@pytest.mark.parametrize("alpha, beta", CDF_PROPERTY_PARAMS) +def test_cdf_is_monotonic(alpha, beta): + dist = Beta(alpha, beta) + values = [dist.cdf_scalar(x) for x in (0.05, 0.2, 0.4, 0.6, 0.8, 0.95)] + assert values == sorted(values) + + +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="beta_cdf_scalar is numerically incorrect") +def test_cdf_is_symmetric_for_symmetric_parameters(): + """For a == b the distribution is symmetric about 0.5, so F(0.5) == 0.5.""" + assert Beta(3.0, 3.0).cdf_scalar(0.5) == pytest.approx(0.5, **ITERATIVE) # --------------------------------------------------------------------------- -# Classmethod validation +# Classmethods # --------------------------------------------------------------------------- +@pytest.mark.parametrize("alpha, beta", PARAMS) +@pytest.mark.parametrize("x", [0.25, 0.5]) +def test_classmethods_agree_with_instances(x, alpha, beta): + dist = Beta(alpha, beta) + assert Beta._pdf_scalar(x, alpha, beta) == pytest.approx( + dist.pdf_scalar(x), **EXACT + ) + assert Beta._mean(alpha, beta) == pytest.approx(dist.mean(), **EXACT) + assert Beta._variance(alpha, beta) == pytest.approx(dist.variance(), **EXACT) + assert Beta._stddev(alpha, beta) == pytest.approx(dist.stddev(), **EXACT) + + @pytest.mark.parametrize("method_name, args", [ ("_pdf_scalar", (0.5, -0.1, 1.0)), ("_pdf_scalar", (0.5, 1.0, -0.1)), @@ -104,36 +230,47 @@ def test_instance_methods_delegate_to_core(mock_core, method_name, core_method_n ("_pdf_scalar", (0.5, 1.0, 0)), ("_cdf_scalar", (0.5, -1, 1.0)), ("_cdf_scalar", (0.5, 1.0, -1)), - ("_mean", (-0.5, 1.0)), - ("_mean", (1.0, -0.5)), - ("_variance", (0, 1.0)), - ("_variance", (1.0, 0)), - ("_stddev", (-1, 1.0)), - ("_stddev", (1.0, -1)), - ("_sample", (-0.1, 1.0)), - ("_sample", (1.0, -0.1)), + ("_mean", (0, 1.0)), + ("_variance", (1.0, -2.0)), + ("_stddev", (-3.0, 1.0)), + ("_sample", (0, 1.0)), ]) def test_classmethods_reject_invalid_parameters(method_name, args): - method = getattr(Beta, method_name) - with pytest.raises(ValueError, match=r"must be positive"): - method.__func__(Beta, *args) + with pytest.raises(ValueError, match="must be positive"): + getattr(Beta, method_name)(*args) # --------------------------------------------------------------------------- -# Classmethod delegation with valid parameters +# Validation edge case +# +# KNOWN BUG: in Beta._validate_params the `if alpha <= 0` check sits outside the +# `if alpha is not None` guard (beta.py line 47), so validating only `beta` +# raises TypeError comparing None to int. # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, args, value", [ - ("_pdf_scalar", "beta_pdf_scalar", (0.5, 2.0, 3.0), 1.5), - ("_cdf_scalar", "beta_cdf_scalar", (0.5, 2.0, 3.0), 0.6875), - ("_mean", "beta_mean", (2.0, 3.0), 0.4), - ("_variance", "beta_variance", (2.0, 3.0), 0.048), - ("_stddev", "beta_stddev", (2.0, 3.0), 0.219), - ("_sample", "beta_sample", (2.0, 3.0), 0.42), -]) -def test_classmethods_delegate_to_core(mock_core, method_name, core_method_name, args, value): - getattr(mock_core, core_method_name).return_value = value - method = getattr(Beta, method_name) - result = method.__func__(Beta, *args) - getattr(mock_core, core_method_name).assert_called_once_with(*args) - assert result == value +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="alpha <= 0 check sits outside the None guard") +def test_validate_params_accepts_a_single_named_parameter(): + Beta._validate_params(beta=3.0) + + +# --------------------------------------------------------------------------- +# Sampling +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("alpha, beta", [(2.0, 3.0), (1.0, 1.0), (5.0, 2.0)]) +def test_sample_lies_within_the_unit_interval(alpha, beta): + dist = Beta(alpha, beta) + for _ in range(200): + value = dist.sample() + assert math.isfinite(value) + assert 0.0 <= value <= 1.0 + + +# --------------------------------------------------------------------------- +# Slots +# --------------------------------------------------------------------------- + +def test_slots_prevent_dynamic_attributes(): + with pytest.raises(AttributeError): + Beta(2.0, 3.0).extra = 123 diff --git a/python/tests/test_binomial.py b/python/tests/test_binomial.py index 4e940df..9daac6f 100644 --- a/python/tests/test_binomial.py +++ b/python/tests/test_binomial.py @@ -1,149 +1,305 @@ -# tests/python_modules/test_beta.py -from unittest.mock import patch +"""Numeric tests for the Binomial distribution against closed-form references.""" + +import math import pytest -import fastdist.distributions.beta as beta_module -from fastdist.distributions.beta import Beta +from conftest import EXACT, ITERATIVE +from fastdist.distributions.binomial import Binomial + + +# --------------------------------------------------------------------------- +# Closed-form references +# --------------------------------------------------------------------------- + +def binom_pmf(k: int, n: int, p: float) -> float: + """P(X = k) = C(n, k) p^k (1-p)^(n-k)""" + return math.comb(n, k) * p ** k * (1 - p) ** (n - k) + + +def binom_cdf(k: int, n: int, p: float) -> float: + return sum(binom_pmf(i, n, p) for i in range(0, k + 1)) + +def binom_mgf(t: float, n: int, p: float) -> float: + """M(t) = (1 - p + p e^t)^n""" + return (1 - p + p * math.exp(t)) ** n -@pytest.fixture -def mock_core(): - """Patch the internal C++ core for the duration of each test.""" - with patch.object(beta_module, "_core", create=True) as mock: - yield mock + +PARAMS = [(1, 0.5), (10, 0.3), (10, 0.7), (5, 0.5), (20, 0.1)] # --------------------------------------------------------------------------- -# Constructor & representation +# Constructor, properties, representation # --------------------------------------------------------------------------- def test_init_valid_parameters(): - b = Beta(alpha=2.0, beta=3.0) - assert b.alpha == 2.0 - assert b.beta == 3.0 - - -@pytest.mark.parametrize("alpha,beta", [ - (-0.1, 1.0), - (0, 1.0), - (-1, 1.0), - (1.0, -0.1), - (1.0, 0), - (1.0, -1), - (-1, -1), -]) -def test_init_invalid_parameters_raises(alpha, beta): - with pytest.raises(ValueError, match="must be positive"): - Beta(alpha=alpha, beta=beta) + dist = Binomial(n=10, p=0.3) + assert dist.n == 10 + assert dist.p == 0.3 + + +@pytest.mark.parametrize("p", [-0.1, 1.1, 2.0]) +def test_init_rejects_out_of_range_p(p): + with pytest.raises(ValueError, match=r"p must be in the interval \[0, 1\]"): + Binomial(n=10, p=p) + + +def test_init_rejects_negative_n(): + with pytest.raises(ValueError, match="n must be a non-negative integer"): + Binomial(n=-1, p=0.5) + + +@pytest.mark.parametrize("bad_n", [2.5, "10"]) +def test_init_rejects_non_integer_n(bad_n): + with pytest.raises(TypeError, match="n must be an integer"): + Binomial(n=bad_n, p=0.5) + + +@pytest.mark.parametrize("bad_p", ["0.5", object()]) +def test_init_rejects_non_real_p(bad_p): + with pytest.raises(TypeError, match="p must be a real number"): + Binomial(n=10, p=bad_p) + + +@pytest.mark.parametrize("kwargs", [{"n": None, "p": 0.5}, {"n": 10, "p": None}]) +def test_init_rejects_none_parameters(kwargs): + """ + None is rejected, but only incidentally. _validate_params skips its type + check when a parameter is None, so the failure surfaces from an unguarded + comparison or from float() rather than as the intended message. + """ + with pytest.raises(TypeError): + Binomial(**kwargs) def test_repr(): - b = Beta(alpha=2.0, beta=3.0) - assert repr(b) == "Beta(alpha=2.0, beta=3.0)" + assert repr(Binomial(n=10, p=0.3)) == "Binomial(n=10, p=0.3)" + + +def test_n_setter_updates_and_validates(): + dist = Binomial(n=10, p=0.3) + dist.n = 20 + assert dist.n == 20 + with pytest.raises(ValueError, match="n must be a non-negative integer"): + dist.n = -5 + + +# KNOWN BUG: the p setter calls _validate_params(p=value), leaving n as None. +# Because the `if n < 0` check sits outside the `if n is not None` guard +# (binomial.py line 47), the comparison None < 0 raises TypeError for *every* +# assignment, valid or not. The p property is unusable. This is the same defect +# pattern as Beta._validate_params. +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="p setter raises TypeError for any value") +def test_p_setter_updates_value(): + dist = Binomial(n=10, p=0.3) + dist.p = 0.6 + assert dist.p == 0.6 + + +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="p setter raises TypeError before validating") +def test_p_setter_rejects_out_of_range(): + dist = Binomial(n=10, p=0.3) + with pytest.raises(ValueError, match=r"p must be in the interval \[0, 1\]"): + dist.p = 1.5 # --------------------------------------------------------------------------- -# Validation behavior +# PMF # --------------------------------------------------------------------------- -@pytest.mark.parametrize("alpha,beta", [ - (0.1, 0.1), - (1.0, 1.0), - (2.5, 3.5), - (100, 200), -]) -def test_validate_params_accepts_valid_values(alpha, beta): - Beta._validate_params(alpha=alpha, beta=beta) +@pytest.mark.parametrize("n, p", PARAMS) +def test_pmf_matches_closed_form(n, p): + dist = Binomial(n, p) + for k in range(n + 1): + assert dist.pmf_scalar(k) == pytest.approx(binom_pmf(k, n, p), **EXACT) -@pytest.mark.parametrize("alpha,beta", [ - (0, 1.0), - (-1, 1.0), - (1.0, 0), - (1.0, -1), -]) -def test_validate_params_rejects_invalid_values(alpha, beta): - with pytest.raises(ValueError, match="must be positive"): - Beta._validate_params(alpha=alpha, beta=beta) +@pytest.mark.parametrize("n, p", PARAMS) +def test_pmf_sums_to_one_over_the_support(n, p): + dist = Binomial(n, p) + total = sum(dist.pmf_scalar(k) for k in range(n + 1)) + assert total == pytest.approx(1.0, **ITERATIVE) + + +@pytest.mark.parametrize("n, p", PARAMS) +def test_pmf_is_a_probability(n, p): + dist = Binomial(n, p) + for k in range(n + 1): + assert 0.0 <= dist.pmf_scalar(k) <= 1.0 + + +def test_pmf_is_symmetric_for_fair_coin(): + """For p = 0.5, P(X = k) == P(X = n - k).""" + dist = Binomial(10, 0.5) + for k in range(11): + assert dist.pmf_scalar(k) == pytest.approx(dist.pmf_scalar(10 - k), **EXACT) # --------------------------------------------------------------------------- -# Instance method delegation +# log-PMF # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, value", [ - ("pdf_scalar", "beta_pdf_scalar", 0.5), - ("cdf_scalar", "beta_cdf_scalar", 0.3), - ("mean", "beta_mean", 0.4), - ("variance", "beta_variance", 0.048), - ("stddev", "beta_stddev", 0.219), - ("sample", "beta_sample", 0.42), -]) -def test_instance_methods_delegate_to_core(mock_core, method_name, core_method_name, value): - getattr(mock_core, core_method_name).return_value = value - b = Beta(alpha=2.0, beta=3.0) - method = getattr(b, method_name) - if "scalar" in method_name: - result = method(0.5) - getattr(mock_core, core_method_name).assert_called_once_with(0.5, 2.0, 3.0) - else: - result = method() - getattr(mock_core, core_method_name).assert_called_once_with(2.0, 3.0) - assert result == value +@pytest.mark.parametrize("n, p", [(10, 0.3), (20, 0.1), (5, 0.5)]) +def test_logpmf_is_log_of_pmf(n, p): + dist = Binomial(n, p) + for k in range(n + 1): + assert dist.logpmf_scalar(k) == pytest.approx( + math.log(binom_pmf(k, n, p)), **ITERATIVE + ) + + +# --------------------------------------------------------------------------- +# CDF +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("n, p", PARAMS) +def test_cdf_matches_summed_pmf(n, p): + dist = Binomial(n, p) + for k in range(n + 1): + assert dist.cdf_scalar(k) == pytest.approx(binom_cdf(k, n, p), **ITERATIVE) + + +@pytest.mark.parametrize("n, p", PARAMS) +def test_cdf_is_monotonic_and_bounded(n, p): + dist = Binomial(n, p) + values = [dist.cdf_scalar(k) for k in range(n + 1)] + assert values == sorted(values) + assert all(0.0 <= v <= 1.0 for v in values) + + +@pytest.mark.parametrize("n, p", PARAMS) +def test_cdf_reaches_one_at_the_top_of_the_support(n, p): + assert Binomial(n, p).cdf_scalar(n) == pytest.approx(1.0, **ITERATIVE) + + +# --------------------------------------------------------------------------- +# Moments +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("n, p", PARAMS) +def test_moments_match_closed_forms(n, p): + dist = Binomial(n, p) + assert dist.mean() == pytest.approx(n * p, **EXACT) + assert dist.variance() == pytest.approx(n * p * (1 - p), **EXACT) + assert dist.stddev() == pytest.approx(math.sqrt(n * p * (1 - p)), **EXACT) + + +@pytest.mark.parametrize("n, p", PARAMS) +def test_mean_equals_pmf_weighted_sum(n, p): + """E[X] = sum k P(X = k), computed independently of the mean() routine.""" + dist = Binomial(n, p) + expectation = sum(k * dist.pmf_scalar(k) for k in range(n + 1)) + assert dist.mean() == pytest.approx(expectation, **ITERATIVE) # --------------------------------------------------------------------------- -# Classmethod validation +# MGF / CGF # --------------------------------------------------------------------------- +@pytest.mark.parametrize("n, p", PARAMS) +@pytest.mark.parametrize("t", [-1.0, -0.1, 0.0, 0.1, 0.5]) +def test_mgf_matches_closed_form(t, n, p): + assert Binomial(n, p).mgf_scalar(t) == pytest.approx(binom_mgf(t, n, p), **EXACT) + + +@pytest.mark.parametrize("n, p", [(10, 0.3), (5, 0.5), (20, 0.1)]) +@pytest.mark.parametrize("t", [-1.0, -0.1, 0.0, 0.1, 0.5]) +def test_cgf_matches_closed_form(t, n, p): + assert Binomial(n, p).cgf_scalar(t) == pytest.approx( + n * math.log(1 - p + p * math.exp(t)), **EXACT + ) + + +@pytest.mark.parametrize("n, p", PARAMS) +def test_mgf_at_zero_is_one(n, p): + dist = Binomial(n, p) + assert dist.mgf_scalar(0.0) == pytest.approx(1.0, **EXACT) + assert dist.cgf_scalar(0.0) == pytest.approx(0.0, abs=1e-12) + + +@pytest.mark.parametrize("n, p", [(10, 0.3), (5, 0.5)]) +def test_cgf_is_log_of_mgf(n, p): + dist = Binomial(n, p) + for t in (-0.4, 0.2, 0.7): + assert dist.cgf_scalar(t) == pytest.approx( + math.log(dist.mgf_scalar(t)), **EXACT + ) + + +# --------------------------------------------------------------------------- +# Classmethods +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("n, p", PARAMS) +def test_classmethods_agree_with_instances(n, p): + dist = Binomial(n, p) + assert Binomial._pmf_scalar(2, n, p) == pytest.approx(dist.pmf_scalar(2), **EXACT) + assert Binomial._cdf_scalar(2, n, p) == pytest.approx( + dist.cdf_scalar(2), **ITERATIVE + ) + assert Binomial._mean(n, p) == pytest.approx(dist.mean(), **EXACT) + assert Binomial._variance(n, p) == pytest.approx(dist.variance(), **EXACT) + assert Binomial._stddev(n, p) == pytest.approx(dist.stddev(), **EXACT) + assert Binomial._mgf_scalar(0.1, n, p) == pytest.approx( + dist.mgf_scalar(0.1), **EXACT + ) + assert Binomial._cgf_scalar(0.1, n, p) == pytest.approx( + dist.cgf_scalar(0.1), **EXACT + ) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pmf_scalar", (1, 10, -0.1)), + ("_pmf_scalar", (1, 10, 1.5)), + ("_cdf_scalar", (1, 10, -0.5)), + ("_mean", (10, 2.0)), + ("_variance", (10, -1.0)), + ("_stddev", (10, 1.5)), + ("_mgf_scalar", (0.1, 10, -0.2)), + ("_cgf_scalar", (0.1, 10, 1.2)), + ("_sample", (10, -0.3)), +]) +def test_classmethods_reject_out_of_range_p(method_name, args): + with pytest.raises(ValueError, match=r"p must be in the interval \[0, 1\]"): + getattr(Binomial, method_name)(*args) + + @pytest.mark.parametrize("method_name, args", [ - ("_pdf_scalar", (0.5, -0.1, 1.0)), - ("_pdf_scalar", (0.5, 1.0, -0.1)), - ("_pdf_scalar", (0.5, 0, 1.0)), - ("_pdf_scalar", (0.5, 1.0, 0)), - ("_cdf_scalar", (0.5, -1, 1.0)), - ("_cdf_scalar", (0.5, 1.0, -1)), - ("_mean", (-0.5, 1.0)), - ("_mean", (1.0, -0.5)), - ("_variance", (0, 1.0)), - ("_variance", (1.0, 0)), - ("_stddev", (-1, 1.0)), - ("_stddev", (1.0, -1)), - ("_sample", (-0.1, 1.0)), - ("_sample", (1.0, -0.1)), + ("_pmf_scalar", (1, -10, 0.5)), + ("_cdf_scalar", (1, -1, 0.5)), + ("_mean", (-5, 0.5)), + ("_sample", (-2, 0.5)), ]) -def test_classmethods_reject_invalid_parameters(method_name, args): - method = getattr(Beta, method_name) - with pytest.raises(ValueError, match=r"must be positive"): - method.__func__(Beta, *args) +def test_classmethods_reject_negative_n(method_name, args): + with pytest.raises(ValueError, match="n must be a non-negative integer"): + getattr(Binomial, method_name)(*args) # --------------------------------------------------------------------------- -# Classmethod delegation with valid parameters +# Sampling # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, args, value", [ - ("_pdf_scalar", "beta_pdf_scalar", (0.5, 2.0, 3.0), 1.5), - ("_cdf_scalar", "beta_cdf_scalar", (0.5, 2.0, 3.0), 0.6875), - ("_mean", "beta_mean", (2.0, 3.0), 0.4), - ("_variance", "beta_variance", (2.0, 3.0), 0.048), - ("_stddev", "beta_stddev", (2.0, 3.0), 0.219), - ("_sample", "beta_sample", (2.0, 3.0), 0.42), -]) -def test_classmethods_delegate_to_core(mock_core, method_name, core_method_name, args, value): - getattr(mock_core, core_method_name).return_value = value - method = getattr(Beta, method_name) - result = method.__func__(Beta, *args) - getattr(mock_core, core_method_name).assert_called_once_with(*args) - assert result == value +@pytest.mark.parametrize("n, p", [(10, 0.3), (5, 0.5), (20, 0.9)]) +def test_sample_lies_within_the_support(n, p): + dist = Binomial(n, p) + for _ in range(300): + value = dist.sample() + assert isinstance(value, int) + assert 0 <= value <= n + + +def test_sample_is_deterministic_at_the_boundaries(): + assert all(Binomial(7, 0.0).sample() == 0 for _ in range(50)) + assert all(Binomial(7, 1.0).sample() == 7 for _ in range(50)) # --------------------------------------------------------------------------- -# Slots behavior +# Slots # --------------------------------------------------------------------------- def test_slots_prevent_dynamic_attributes(): - b = Beta(alpha=2.0, beta=3.0) with pytest.raises(AttributeError): - b.extra = 123 + Binomial(10, 0.3).extra = 123 diff --git a/python/tests/test_chi_square.py b/python/tests/test_chi_square.py index 33d109e..fd61a29 100644 --- a/python/tests/test_chi_square.py +++ b/python/tests/test_chi_square.py @@ -1,41 +1,69 @@ -# tests/python_modules/test_chi_square.py -from unittest.mock import patch +"""Numeric tests for the Chi-square distribution against closed-form references.""" + +import math import pytest -import fastdist.distributions.chi_square as chi_square_module +from conftest import EXACT, ITERATIVE, regularized_lower_gamma from fastdist.distributions.chi_square import ChiSquare -@pytest.fixture -def mock_core(): - """Patch the internal C++ core for the duration of each test.""" - with patch.object(chi_square_module, "_core", create=True) as mock: - yield mock +# --------------------------------------------------------------------------- +# Closed-form references +# --------------------------------------------------------------------------- + +def chi2_pdf(x: float, k: float) -> float: + """f(x; k) = x^(k/2-1) e^(-x/2) / (2^(k/2) Gamma(k/2))""" + return x ** (k / 2 - 1) * math.exp(-x / 2) / (2 ** (k / 2) * math.gamma(k / 2)) + + +def chi2_mgf(t: float, k: float) -> float: + """M(t) = (1 - 2t)^(-k/2), valid for t < 1/2""" + return (1 - 2 * t) ** (-k / 2) + + +def chi2_cgf(t: float, k: float) -> float: + """K(t) = -(k/2) ln(1 - 2t)""" + return -(k / 2) * math.log(1 - 2 * t) + + +DEGREES = [1.0, 2.0, 3.0, 5.0, 10.0] # --------------------------------------------------------------------------- -# Constructor & representation +# Constructor, properties, representation # --------------------------------------------------------------------------- def test_init_valid_parameter(): - cs = ChiSquare(k=5.0) - assert cs.k == 5.0 + assert ChiSquare(k=5.0).k == 5.0 @pytest.mark.parametrize("k", [-0.1, 0, -1, -10]) -def test_init_invalid_k_raises(k): +def test_init_rejects_non_positive_k(k): with pytest.raises(ValueError, match="k must be positive"): ChiSquare(k=k) +@pytest.mark.parametrize("bad", ["5.0", None, object()]) +def test_init_rejects_non_real_k(bad): + with pytest.raises(TypeError, match="k must be a real number"): + ChiSquare(k=bad) + + def test_repr(): - cs = ChiSquare(k=3.0) - assert repr(cs) == "ChiSquare(k=3.0)" + assert repr(ChiSquare(k=3.0)) == "ChiSquare(k=3.0)" + + +def test_k_setter_updates_and_validates(): + dist = ChiSquare(k=5.0) + dist.k = 8.0 + assert dist.k == 8.0 + with pytest.raises(ValueError, match="k must be positive"): + dist.k = -1.0 # --------------------------------------------------------------------------- -# Validation behavior +# Validation # --------------------------------------------------------------------------- @pytest.mark.parametrize("k", [0.1, 1.0, 5.0, 100.0]) @@ -50,146 +78,203 @@ def test_validate_params_rejects_invalid_values(k): # --------------------------------------------------------------------------- -# Instance method delegation +# PDF # --------------------------------------------------------------------------- -def test_pdf_delegates_to_core(mock_core): - mock_core.chi_square_pdf_scalar.return_value = 0.154 +@pytest.mark.parametrize("k", DEGREES) +@pytest.mark.parametrize("x", [0.5, 1.0, 3.0, 7.5]) +def test_pdf_matches_closed_form(x, k): + assert ChiSquare(k).pdf(x) == pytest.approx(chi2_pdf(x, k), **EXACT) - cs = ChiSquare(k=5.0) - result = cs.pdf(3.0) - mock_core.chi_square_pdf_scalar.assert_called_once_with(3.0, 5.0) - assert result == 0.154 +@pytest.mark.parametrize("k", DEGREES) +@pytest.mark.parametrize("x", [0.1, 1.0, 5.0, 20.0]) +def test_pdf_is_non_negative(x, k): + assert ChiSquare(k).pdf(x) >= 0.0 -def test_cdf_delegates_to_core(mock_core): - mock_core.chi_square_cdf_scalar.return_value = 0.416 +@pytest.mark.parametrize("x", [0.5, 1.0, 4.0]) +def test_pdf_of_two_degrees_is_exponential(x): + """Chi-square with k=2 is Exponential with mean 2, so f(x) = e^(-x/2) / 2.""" + assert ChiSquare(2.0).pdf(x) == pytest.approx(math.exp(-x / 2) / 2, **EXACT) - cs = ChiSquare(k=5.0) - result = cs.cdf(3.0) - mock_core.chi_square_cdf_scalar.assert_called_once_with(3.0, 5.0) - assert result == 0.416 +# --------------------------------------------------------------------------- +# CDF +# --------------------------------------------------------------------------- + +# KNOWN BUG: chi_square_cdf_scalar delegates to the same regularized lower +# incomplete gamma as Gamma.cdf_scalar, whose continued-fraction branch is +# incorrect. ChiSquare(3.0).cdf(7.5) returns 1.000498004, a probability greater +# than one. The failing points below were determined empirically against +# conftest.regularized_lower_gamma; k = 2 is correct throughout because that +# case reduces to the exact exponential form. +CF_BUG = ( + "regularized lower incomplete gamma is incorrect in the " + "continued-fraction branch (x/2 >= k/2 + 1)" +) +CHI2_X = (0.5, 1.0, 3.0, 7.5, 20.0, 50.0) -def test_mean_delegates_to_core(mock_core): - mock_core.chi_square_mean.return_value = 5.0 +_CDF_KNOWN_BAD = { + (1.0, 3.0), (1.0, 7.5), (1.0, 20.0), + (3.0, 7.5), (3.0, 20.0), + (5.0, 7.5), (5.0, 20.0), (5.0, 50.0), + (10.0, 20.0), (10.0, 50.0), +} - cs = ChiSquare(k=5.0) - result = cs.mean() - mock_core.chi_square_mean.assert_called_once_with(5.0) - assert result == 5.0 +def _cdf_case(k, x): + marks = [] + if (k, x) in _CDF_KNOWN_BAD: + marks = [ + pytest.mark.known_bug, + pytest.mark.xfail(strict=True, reason=CF_BUG), + ] + return pytest.param(k, x, marks=marks) -def test_variance_delegates_to_core(mock_core): - mock_core.chi_square_variance.return_value = 10.0 +CDF_CASES = [_cdf_case(k, x) for k in DEGREES for x in CHI2_X] - cs = ChiSquare(k=5.0) - result = cs.variance() - mock_core.chi_square_variance.assert_called_once_with(5.0) - assert result == 10.0 +@pytest.mark.parametrize("k, x", CDF_CASES) +def test_cdf_matches_reference(k, x): + assert ChiSquare(k).cdf(x) == pytest.approx( + regularized_lower_gamma(k / 2.0, x / 2.0), **ITERATIVE + ) -def test_stddev_delegates_to_core(mock_core): - mock_core.chi_square_stddev.return_value = 3.162 +@pytest.mark.parametrize("x", [0.5, 1.0, 3.0, 7.5, 20.0]) +def test_cdf_exact_for_two_degrees_of_freedom(x): + """For k=2 the CDF has the exact closed form 1 - e^(-x/2).""" + assert ChiSquare(2.0).cdf(x) == pytest.approx(1 - math.exp(-x / 2), **ITERATIVE) - cs = ChiSquare(k=5.0) - result = cs.stddev() - mock_core.chi_square_stddev.assert_called_once_with(5.0) - assert result == 3.162 +_CDF_PROPERTY_DEGREES = [ + 1.0, + 2.0, + pytest.param(3.0, marks=[ + pytest.mark.known_bug, + pytest.mark.xfail(strict=True, reason=CF_BUG + "; CDF exceeds 1.0"), + ]), + 5.0, + 10.0, +] -def test_mgf_scalar_delegates_to_core(mock_core): - mock_core.chi_square_mgf_scalar.return_value = 1.789 +@pytest.mark.parametrize("k", _CDF_PROPERTY_DEGREES) +def test_cdf_is_bounded(k): + dist = ChiSquare(k) + for x in (0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 30.0): + assert 0.0 <= dist.cdf(x) <= 1.0 - cs = ChiSquare(k=5.0) - result = cs.mgf_scalar(0.1) - mock_core.chi_square_mgf_scalar.assert_called_once_with(0.1, 5.0) - assert result == 1.789 +@pytest.mark.parametrize("k", _CDF_PROPERTY_DEGREES) +def test_cdf_is_monotonic(k): + dist = ChiSquare(k) + values = [dist.cdf(x) for x in (0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 30.0)] + assert values == sorted(values) -def test_cgf_scalar_delegates_to_core(mock_core): - mock_core.chi_square_cgf_scalar.return_value = 0.581 +# --------------------------------------------------------------------------- +# Moments +# --------------------------------------------------------------------------- - cs = ChiSquare(k=5.0) - result = cs.cgf_scalar(0.1) +@pytest.mark.parametrize("k", [1.0, 2.0, 5.0, 12.5, 100.0]) +def test_moments_match_closed_forms(k): + dist = ChiSquare(k) + assert dist.mean() == pytest.approx(k, **EXACT) + assert dist.variance() == pytest.approx(2 * k, **EXACT) + assert dist.stddev() == pytest.approx(math.sqrt(2 * k), **EXACT) - mock_core.chi_square_cgf_scalar.assert_called_once_with(0.1, 5.0) - assert result == 0.581 +# --------------------------------------------------------------------------- +# MGF / CGF +# --------------------------------------------------------------------------- -def test_sample_delegates_to_core(mock_core): - mock_core.chi_square_sample.return_value = 4.823 +@pytest.mark.parametrize("k", [1.0, 2.0, 5.0]) +@pytest.mark.parametrize("t", [-0.5, -0.1, 0.0, 0.1, 0.2]) +def test_mgf_matches_closed_form(t, k): + assert ChiSquare(k).mgf_scalar(t) == pytest.approx(chi2_mgf(t, k), **EXACT) - cs = ChiSquare(k=5.0) - result = cs.sample() - mock_core.chi_square_sample.assert_called_once_with(5.0) - assert result == 4.823 +@pytest.mark.parametrize("k", [1.0, 2.0, 5.0]) +@pytest.mark.parametrize("t", [-0.5, -0.1, 0.0, 0.1, 0.2]) +def test_cgf_matches_closed_form(t, k): + assert ChiSquare(k).cgf_scalar(t) == pytest.approx(chi2_cgf(t, k), **EXACT) + + +@pytest.mark.parametrize("k", [1.0, 3.0, 8.0]) +def test_mgf_at_zero_is_one(k): + assert ChiSquare(k).mgf_scalar(0.0) == pytest.approx(1.0, **EXACT) + assert ChiSquare(k).cgf_scalar(0.0) == pytest.approx(0.0, abs=1e-12) + + +@pytest.mark.parametrize("k", [1.0, 5.0]) +@pytest.mark.parametrize("t", [-0.3, 0.1, 0.25]) +def test_cgf_is_log_of_mgf(t, k): + """K(t) = ln M(t) - catches a sign or factor error in either one.""" + dist = ChiSquare(k) + assert dist.cgf_scalar(t) == pytest.approx(math.log(dist.mgf_scalar(t)), **EXACT) # --------------------------------------------------------------------------- -# Validation enforcement in classmethods +# Classmethods # --------------------------------------------------------------------------- -@pytest.mark.parametrize( - "method, args", - [ - (ChiSquare._pdf_scalar, (1.0, -1)), - (ChiSquare._pdf_scalar, (1.0, 0)), - (ChiSquare._cdf_scalar, (1.0, -5.0)), - (ChiSquare._cdf_scalar, (1.0, 0)), - (ChiSquare._mean, (-1,)), - (ChiSquare._mean, (0,)), - (ChiSquare._variance, (-5.0,)), - (ChiSquare._variance, (0,)), - (ChiSquare._stddev, (-1,)), - (ChiSquare._stddev, (0,)), - (ChiSquare._mgf_scalar, (0.1, -1)), - (ChiSquare._mgf_scalar, (0.1, 0)), - (ChiSquare._cgf_scalar, (0.1, -5.0)), - (ChiSquare._cgf_scalar, (0.1, 0)), - (ChiSquare._sample, (-1,)), - (ChiSquare._sample, (0,)), - ], -) -def test_classmethods_reject_invalid_parameters(method, args): +@pytest.mark.parametrize("k", [1.0, 5.0]) +@pytest.mark.parametrize("x", [0.5, 3.0]) +def test_classmethods_agree_with_instances(x, k): + dist = ChiSquare(k) + assert ChiSquare._pdf_scalar(x, k) == pytest.approx(dist.pdf(x), **EXACT) + assert ChiSquare._cdf_scalar(x, k) == pytest.approx(dist.cdf(x), **ITERATIVE) + assert ChiSquare._mean(k) == pytest.approx(dist.mean(), **EXACT) + assert ChiSquare._variance(k) == pytest.approx(dist.variance(), **EXACT) + assert ChiSquare._stddev(k) == pytest.approx(dist.stddev(), **EXACT) + assert ChiSquare._mgf_scalar(0.1, k) == pytest.approx(dist.mgf_scalar(0.1), **EXACT) + assert ChiSquare._cgf_scalar(0.1, k) == pytest.approx(dist.cgf_scalar(0.1), **EXACT) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pdf_scalar", (1.0, -1)), + ("_pdf_scalar", (1.0, 0)), + ("_cdf_scalar", (1.0, -5.0)), + ("_cdf_scalar", (1.0, 0)), + ("_mean", (-1,)), + ("_mean", (0,)), + ("_variance", (-5.0,)), + ("_variance", (0,)), + ("_stddev", (-1,)), + ("_stddev", (0,)), + ("_mgf_scalar", (0.1, -1)), + ("_mgf_scalar", (0.1, 0)), + ("_cgf_scalar", (0.1, -5.0)), + ("_cgf_scalar", (0.1, 0)), + ("_sample", (-1,)), + ("_sample", (0,)), +]) +def test_classmethods_reject_invalid_parameters(method_name, args): with pytest.raises(ValueError, match=r"k must be positive"): - method(*args) + getattr(ChiSquare, method_name)(*args) # --------------------------------------------------------------------------- -# Classmethod delegation with valid parameters +# Sampling # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, args, value", [ - ("_pdf_scalar", "chi_square_pdf_scalar", (3.0, 5.0), 0.154), - ("_cdf_scalar", "chi_square_cdf_scalar", (3.0, 5.0), 0.416), - ("_mean", "chi_square_mean", (5.0,), 5.0), - ("_variance", "chi_square_variance", (5.0,), 10.0), - ("_stddev", "chi_square_stddev", (5.0,), 3.162), - ("_mgf_scalar", "chi_square_mgf_scalar", (0.1, 5.0), 1.789), - ("_cgf_scalar", "chi_square_cgf_scalar", (0.1, 5.0), 0.581), - ("_sample", "chi_square_sample", (5.0,), 4.823), -]) -def test_classmethods_delegate_to_core(mock_core, method_name, core_method_name, args, value): - getattr(mock_core, core_method_name).return_value = value - method = getattr(ChiSquare, method_name) - result = method(*args) - getattr(mock_core, core_method_name).assert_called_once_with(*args) - assert result == value +@pytest.mark.parametrize("k", [1.0, 5.0, 20.0]) +def test_sample_is_positive_and_finite(k): + dist = ChiSquare(k) + for _ in range(200): + value = dist.sample() + assert math.isfinite(value) + assert value > 0.0 # --------------------------------------------------------------------------- -# Slots behavior +# Slots # --------------------------------------------------------------------------- def test_slots_prevent_dynamic_attributes(): - cs = ChiSquare(k=5.0) with pytest.raises(AttributeError): - cs.extra = 123 + ChiSquare(k=5.0).extra = 123 diff --git a/python/tests/test_discrete_uniform.py b/python/tests/test_discrete_uniform.py index bc9d311..8181a6f 100644 --- a/python/tests/test_discrete_uniform.py +++ b/python/tests/test_discrete_uniform.py @@ -1,195 +1,274 @@ -# tests/python_modules/test_discrete_uniform.py -from unittest.mock import patch +""" +Numeric tests for the Discrete Uniform distribution against closed-form +references. -import pytest +The support is the integers a, a+1, ..., b inclusive, so there are n = b - a + 1 +equally likely outcomes. +""" -import fastdist.distributions.discrete_uniform as du_module -from fastdist.distributions.discrete_uniform import DiscreteUniform +import math +import pytest -@pytest.fixture -def mock_core(): - """ - Patch the internal C++ core for the duration of each test. - """ - with patch.object(du_module, "_core", autospec=True) as mock: - yield mock +from conftest import EXACT, ITERATIVE +from fastdist.distributions.discrete_uniform import DiscreteUniform # --------------------------------------------------------------------------- -# Constructor & representation +# Closed-form references # --------------------------------------------------------------------------- -def test_init_valid_parameters(): - du = DiscreteUniform(a=1, b=10) - assert du.a == 1 - assert du.b == 10 - - -@pytest.mark.parametrize( - "a,b", - [ - (5, 5), # equal - (6, 5), # a > b - (10, 1), # a > b - ] -) -def test_init_invalid_range_raises(a, b): - with pytest.raises(ValueError, match="a must be less than b"): - DiscreteUniform(a=a, b=b) +def du_support_size(a: int, b: int) -> int: + return b - a + 1 -def test_repr(): - du = DiscreteUniform(a=3, b=8) - assert repr(du) == "DiscreteUniform(a=3, b=8)" +def du_mean(a: int, b: int) -> float: + return (a + b) / 2.0 -# --------------------------------------------------------------------------- -# Validation behavior -# --------------------------------------------------------------------------- +def du_variance(a: int, b: int) -> float: + """Var = ((b - a + 1)^2 - 1) / 12""" + n = du_support_size(a, b) + return (n ** 2 - 1) / 12.0 + + +def du_mgf(t: float, a: int, b: int) -> float: + """M(t) = (1/n) sum_{k=a}^{b} e^(tk); equals 1 at t = 0.""" + n = du_support_size(a, b) + return sum(math.exp(t * k) for k in range(a, b + 1)) / n -def test_validate_params_accepts_valid_ranges(): - DiscreteUniform._validate_params(a=1, b=2) - DiscreteUniform._validate_params(a=-10, b=-5) + +PARAMS = [(1, 6), (0, 1), (-3, 3), (2, 10), (-5, -1)] # --------------------------------------------------------------------------- -# Instance method delegation +# Constructor, properties, representation # --------------------------------------------------------------------------- -def test_pmf_delegates_to_core(mock_core): - mock_core.discrete_uniform_pmf_scalar.return_value = 0.1 +def test_init_valid_parameters(): + dist = DiscreteUniform(a=1, b=6) + assert dist.a == 1 + assert dist.b == 6 + + +@pytest.mark.parametrize("a, b", [(5, 5), (6, 1), (0, -1)]) +def test_init_rejects_non_increasing_bounds(a, b): + with pytest.raises(ValueError, match="a must be less than b"): + DiscreteUniform(a=a, b=b) + + +@pytest.mark.parametrize("bad", [1.5, "1"]) +def test_init_rejects_non_integer_bounds(bad): + with pytest.raises(TypeError, match="a must be an integer"): + DiscreteUniform(a=bad, b=10) + with pytest.raises(TypeError, match="b must be an integer"): + DiscreteUniform(a=1, b=bad) + + +def test_repr(): + assert repr(DiscreteUniform(a=1, b=6)) == "DiscreteUniform(a=1, b=6)" - du = DiscreteUniform(a=1, b=5) - result = du.pmf(3) - mock_core.discrete_uniform_pmf_scalar.assert_called_once_with(3, 1, 5) - assert result == 0.1 +# KNOWN BUG: the b setter assigns `self.b = value` instead of `self._b = value` +# (discrete_uniform.py), so it re-enters itself and raises RecursionError for +# every assignment. The b property is unusable. +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="b setter recurses into itself (self.b = value)") +def test_b_setter_updates_value(): + dist = DiscreteUniform(a=1, b=6) + dist.b = 10 + assert dist.b == 10 -def test_cdf_delegates_to_core(mock_core): - mock_core.discrete_uniform_cdf_scalar.return_value = 0.7 +# KNOWN BUG: the a setter stores float(value) even though a is an integer +# parameter that __init__ stores via int(). Setting a therefore changes the +# attribute's type from int to float. +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="a setter stores float(value) instead of int") +def test_a_setter_preserves_integer_type(): + dist = DiscreteUniform(a=1, b=6) + dist.a = 2 + assert dist.a == 2 + assert isinstance(dist.a, int) - du = DiscreteUniform(a=1, b=5) - result = du.cdf(3) - mock_core.discrete_uniform_cdf_scalar.assert_called_once_with(3, 1, 5) - assert result == 0.7 +def test_a_setter_validates(): + dist = DiscreteUniform(a=1, b=6) + with pytest.raises(TypeError, match="a must be an integer"): + dist.a = 1.5 # --------------------------------------------------------------------------- -# Statistical properties +# PMF # --------------------------------------------------------------------------- -def test_mean_delegates_to_core(mock_core): - mock_core.discrete_uniform_mean.return_value = 2.5 +@pytest.mark.parametrize("a, b", PARAMS) +def test_pmf_is_uniform_over_the_support(a, b): + dist = DiscreteUniform(a, b) + expected = 1.0 / du_support_size(a, b) + for k in range(a, b + 1): + assert dist.pmf(k) == pytest.approx(expected, **EXACT) - du = DiscreteUniform(a=1, b=4) - result = du.mean() - mock_core.discrete_uniform_mean.assert_called_once_with(1, 4) - assert result == 2.5 +@pytest.mark.parametrize("a, b", PARAMS) +def test_pmf_sums_to_one(a, b): + dist = DiscreteUniform(a, b) + total = sum(dist.pmf(k) for k in range(a, b + 1)) + assert total == pytest.approx(1.0, **EXACT) -def test_variance_delegates_to_core(mock_core): - mock_core.discrete_uniform_variance.return_value = 1.25 +@pytest.mark.parametrize("a, b", PARAMS) +def test_pmf_is_zero_outside_the_support(a, b): + dist = DiscreteUniform(a, b) + assert dist.pmf(a - 1) == pytest.approx(0.0, abs=1e-15) + assert dist.pmf(b + 1) == pytest.approx(0.0, abs=1e-15) + + +# --------------------------------------------------------------------------- +# CDF +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("a, b", PARAMS) +def test_cdf_matches_closed_form(a, b): + """F(k) = (k - a + 1) / n on the support.""" + dist = DiscreteUniform(a, b) + n = du_support_size(a, b) + for k in range(a, b + 1): + assert dist.cdf(k) == pytest.approx((k - a + 1) / n, **EXACT) - du = DiscreteUniform(a=1, b=4) - result = du.variance() - mock_core.discrete_uniform_variance.assert_called_once_with(1, 4) - assert result == 1.25 +@pytest.mark.parametrize("a, b", PARAMS) +def test_cdf_matches_summed_pmf(a, b): + dist = DiscreteUniform(a, b) + for k in range(a, b + 1): + assert dist.cdf(k) == pytest.approx( + sum(dist.pmf(i) for i in range(a, k + 1)), **ITERATIVE + ) -def test_stddev_delegates_to_core(mock_core): - mock_core.discrete_uniform_stddev.return_value = 1.118 +@pytest.mark.parametrize("a, b", PARAMS) +def test_cdf_is_monotonic_and_bounded(a, b): + dist = DiscreteUniform(a, b) + values = [dist.cdf(k) for k in range(a - 2, b + 3)] + assert values == sorted(values) + assert all(0.0 <= v <= 1.0 for v in values) - du = DiscreteUniform(a=1, b=4) - result = du.stddev() - mock_core.discrete_uniform_stddev.assert_called_once_with(1, 4) - assert result == 1.118 +@pytest.mark.parametrize("a, b", PARAMS) +def test_cdf_saturates_at_the_bounds(a, b): + dist = DiscreteUniform(a, b) + assert dist.cdf(a - 1) == pytest.approx(0.0, abs=1e-15) + assert dist.cdf(b) == pytest.approx(1.0, **EXACT) + assert dist.cdf(b + 5) == pytest.approx(1.0, **EXACT) -def test_mgf_delegates_to_core(mock_core): - mock_core.discrete_uniform_mgf_scalar.return_value = 2.345 +# --------------------------------------------------------------------------- +# Moments +# --------------------------------------------------------------------------- - du = DiscreteUniform(a=1, b=5) - result = du.mgf(0.5) +@pytest.mark.parametrize("a, b", PARAMS) +def test_moments_match_closed_forms(a, b): + dist = DiscreteUniform(a, b) + assert dist.mean() == pytest.approx(du_mean(a, b), **EXACT) + assert dist.variance() == pytest.approx(du_variance(a, b), **EXACT) + assert dist.stddev() == pytest.approx(math.sqrt(du_variance(a, b)), **EXACT) - mock_core.discrete_uniform_mgf_scalar.assert_called_once_with(0.5, 1, 5) - assert result == 2.345 +@pytest.mark.parametrize("a, b", PARAMS) +def test_mean_equals_pmf_weighted_sum(a, b): + dist = DiscreteUniform(a, b) + expectation = sum(k * dist.pmf(k) for k in range(a, b + 1)) + assert dist.mean() == pytest.approx(expectation, **ITERATIVE) -def test_cgf_delegates_to_core(mock_core): - mock_core.discrete_uniform_cgf_scalar.return_value = 0.852 - du = DiscreteUniform(a=1, b=5) - result = du.cgf(0.5) +# --------------------------------------------------------------------------- +# MGF / CGF +# --------------------------------------------------------------------------- - mock_core.discrete_uniform_cgf_scalar.assert_called_once_with(0.5, 1, 5) - assert result == 0.852 +@pytest.mark.parametrize("a, b", PARAMS) +@pytest.mark.parametrize("t", [-0.5, -0.1, 0.1, 0.5]) +def test_mgf_matches_closed_form(t, a, b): + assert DiscreteUniform(a, b).mgf(t) == pytest.approx(du_mgf(t, a, b), **ITERATIVE) -def test_sample_delegates_to_core(mock_core): - mock_core.discrete_uniform_sample.return_value = 3 +@pytest.mark.parametrize("a, b", PARAMS) +def test_mgf_at_zero_is_one(a, b): + """The t = 0 case is a removable 0/0 singularity in the usual closed form.""" + dist = DiscreteUniform(a, b) + assert dist.mgf(0.0) == pytest.approx(1.0, **EXACT) + assert dist.cgf(0.0) == pytest.approx(0.0, abs=1e-12) - du = DiscreteUniform(a=1, b=5) - result = du.sample() - mock_core.discrete_uniform_sample.assert_called_once_with(1, 5) - assert result == 3 +@pytest.mark.parametrize("a, b", PARAMS) +@pytest.mark.parametrize("t", [-0.5, -0.1, 0.1, 0.5]) +def test_cgf_is_log_of_mgf(t, a, b): + dist = DiscreteUniform(a, b) + assert dist.cgf(t) == pytest.approx(math.log(dist.mgf(t)), **ITERATIVE) # --------------------------------------------------------------------------- -# Validation enforcement in classmethods +# Classmethods # --------------------------------------------------------------------------- -@pytest.mark.parametrize( - "method, args", - [ - (DiscreteUniform._pmf_scalar, (3, 5, 2)), - (DiscreteUniform._cdf_scalar, (3, 6, 1)), - (DiscreteUniform._mean, (10, 1)), - (DiscreteUniform._variance, (5, 4)), - (DiscreteUniform._stddev, (7, 3)), - (DiscreteUniform._mgf_scalar, (0.5, 5, 2)), - (DiscreteUniform._cgf_scalar, (0.5, 6, 1)), - (DiscreteUniform._sample, (10, 1)), - ] -) -def test_classmethods_reject_invalid_range(method, args): +@pytest.mark.parametrize("a, b", PARAMS) +def test_classmethods_agree_with_instances(a, b): + dist = DiscreteUniform(a, b) + assert DiscreteUniform._pmf_scalar(a, a, b) == pytest.approx( + dist.pmf(a), **EXACT + ) + assert DiscreteUniform._cdf_scalar(a, a, b) == pytest.approx( + dist.cdf(a), **EXACT + ) + assert DiscreteUniform._mean(a, b) == pytest.approx(dist.mean(), **EXACT) + assert DiscreteUniform._variance(a, b) == pytest.approx(dist.variance(), **EXACT) + assert DiscreteUniform._stddev(a, b) == pytest.approx(dist.stddev(), **EXACT) + assert DiscreteUniform._mgf_scalar(0.1, a, b) == pytest.approx( + dist.mgf(0.1), **ITERATIVE + ) + assert DiscreteUniform._cgf_scalar(0.1, a, b) == pytest.approx( + dist.cgf(0.1), **ITERATIVE + ) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pmf_scalar", (1, 6, 1)), + ("_cdf_scalar", (1, 6, 1)), + ("_mean", (6, 1)), + ("_variance", (5, 5)), + ("_stddev", (10, 2)), + ("_mgf_scalar", (0.1, 6, 1)), + ("_cgf_scalar", (0.1, 6, 1)), + ("_sample", (6, 1)), +]) +def test_classmethods_reject_non_increasing_bounds(method_name, args): with pytest.raises(ValueError, match="a must be less than b"): - method(*args) + getattr(DiscreteUniform, method_name)(*args) # --------------------------------------------------------------------------- -# Classmethod delegation with valid parameters +# Sampling # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, args, value", [ - ("_pmf_scalar", "discrete_uniform_pmf_scalar", (3, 1, 5), 0.25), - ("_cdf_scalar", "discrete_uniform_cdf_scalar", (3, 1, 5), 0.75), - ("_mean", "discrete_uniform_mean", (1, 5), 3.0), - ("_variance", "discrete_uniform_variance", (1, 5), 1.333), - ("_stddev", "discrete_uniform_stddev", (1, 5), 1.155), - ("_mgf_scalar", "discrete_uniform_mgf_scalar", (0.5, 1, 5), 2.345), - ("_cgf_scalar", "discrete_uniform_cgf_scalar", (0.5, 1, 5), 0.852), - ("_sample", "discrete_uniform_sample", (1, 5), 3), -]) -def test_classmethods_delegate_to_core(mock_core, method_name, core_method_name, args, value): - getattr(mock_core, core_method_name).return_value = value - method = getattr(DiscreteUniform, method_name) - result = method(*args) - getattr(mock_core, core_method_name).assert_called_once_with(*args) - assert result == value +@pytest.mark.parametrize("a, b", [(1, 6), (-3, 3), (0, 1)]) +def test_sample_lies_within_the_support(a, b): + dist = DiscreteUniform(a, b) + for _ in range(300): + value = dist.sample() + assert isinstance(value, int) + assert a <= value <= b + + +def test_sample_eventually_covers_the_whole_support(): + """With 2000 draws from a 6-point support, every value should appear.""" + dist = DiscreteUniform(1, 6) + seen = {dist.sample() for _ in range(2000)} + assert seen == {1, 2, 3, 4, 5, 6} # --------------------------------------------------------------------------- -# Slots behavior +# Slots # --------------------------------------------------------------------------- def test_slots_prevent_dynamic_attributes(): - du = DiscreteUniform(a=1, b=5) with pytest.raises(AttributeError): - du.extra = 123 + DiscreteUniform(1, 6).extra = 123 diff --git a/python/tests/test_exponential.py b/python/tests/test_exponential.py index dc9ad9a..cd0a883 100644 --- a/python/tests/test_exponential.py +++ b/python/tests/test_exponential.py @@ -1,137 +1,303 @@ -# tests/python_modules/test_exponential.py -from unittest.mock import patch +""" +Numeric tests for the Exponential distribution against closed-form references. +`lambda_` is the rate parameter: mean = 1/lambda, f(x) = lambda e^(-lambda x). +""" + +import math + +import numpy as np import pytest -import fastdist.distributions.exponential as exp_module +from conftest import EXACT, ITERATIVE from fastdist.distributions.exponential import Exponential -@pytest.fixture -def mock_core(): - with patch.object(exp_module, "_core", create=True) as mock: - yield mock +# --------------------------------------------------------------------------- +# Closed-form references +# --------------------------------------------------------------------------- + +def exp_pdf(x: float, lam: float) -> float: + return lam * math.exp(-lam * x) + + +def exp_cdf(x: float, lam: float) -> float: + return 1.0 - math.exp(-lam * x) + + +def exp_mgf(t: float, lam: float) -> float: + """M(t) = lambda / (lambda - t), valid for t < lambda""" + return lam / (lam - t) + + +RATES = [0.5, 1.0, 2.0, 5.0] +XS = [0.0, 0.25, 1.0, 3.0, 7.5] # --------------------------------------------------------------------------- -# Constructor & representation +# Constructor, properties, representation # --------------------------------------------------------------------------- def test_init_valid_parameter(): - e = Exponential(lambda_=0.5) - assert e.lambda_ == 0.5 + assert Exponential(lambda_=2.0).lambda_ == 2.0 -@pytest.mark.parametrize("lambda_", [0, -1, -0.5]) -def test_init_invalid_lambda_raises(lambda_): +@pytest.mark.parametrize("lam", [0.0, -0.1, -5.0]) +def test_init_rejects_non_positive_rate(lam): with pytest.raises(ValueError, match="lambda_ must be positive"): - Exponential(lambda_=lambda_) + Exponential(lambda_=lam) + + +@pytest.mark.parametrize("bad", ["2.0", None, object()]) +def test_init_rejects_non_real_rate(bad): + with pytest.raises(TypeError, match="lambda_ must be a real number"): + Exponential(lambda_=bad) def test_repr(): - e = Exponential(lambda_=0.7) - assert repr(e) == "Exponential(lambda_=0.7)" + assert repr(Exponential(lambda_=2.0)) == "Exponential(lambda_=2.0)" + + +def test_rate_setter_updates_and_validates(): + dist = Exponential(lambda_=2.0) + dist.lambda_ = 4.0 + assert dist.lambda_ == 4.0 + with pytest.raises(ValueError, match="lambda_ must be positive"): + dist.lambda_ = 0.0 # --------------------------------------------------------------------------- -# Validation behavior +# Scalar PDF / CDF # --------------------------------------------------------------------------- -@pytest.mark.parametrize("lambda_", [0.0001, 1, 10]) -def test_validate_params_accepts_valid_values(lambda_): - Exponential._validate_params(lambda_=lambda_) +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("x", XS) +def test_pdf_matches_closed_form(x, lam): + assert Exponential(lam).pdf(x) == pytest.approx(exp_pdf(x, lam), **EXACT) + + +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("x", XS) +def test_cdf_matches_closed_form(x, lam): + assert Exponential(lam).cdf(x) == pytest.approx(exp_cdf(x, lam), **EXACT) + + +@pytest.mark.parametrize("lam", RATES) +def test_pdf_at_zero_equals_the_rate(lam): + assert Exponential(lam).pdf(0.0) == pytest.approx(lam, **EXACT) + + +@pytest.mark.parametrize("lam", RATES) +def test_cdf_is_monotonic_and_bounded(lam): + dist = Exponential(lam) + values = [dist.cdf(x) for x in (0.0, 0.1, 0.5, 1.0, 4.0, 20.0)] + assert values == sorted(values) + assert all(0.0 <= v <= 1.0 for v in values) + + +@pytest.mark.parametrize("lam", RATES) +def test_cdf_starts_at_zero_and_saturates(lam): + dist = Exponential(lam) + assert dist.cdf(0.0) == pytest.approx(0.0, abs=1e-15) + assert dist.cdf(500.0 / lam) == pytest.approx(1.0, abs=1e-12) + + +@pytest.mark.parametrize("lam", RATES) +def test_memorylessness(lam): + """P(X > s + t | X > s) == P(X > t), the defining property.""" + dist = Exponential(lam) + s, t = 1.5, 2.5 + survival = lambda z: 1.0 - dist.cdf(z) + assert survival(s + t) / survival(s) == pytest.approx(survival(t), **ITERATIVE) # --------------------------------------------------------------------------- -# Class method delegation to core +# Array API # --------------------------------------------------------------------------- -def test_pdf_scalar_delegates_to_core(mock_core): - mock_core.exponential_pdf_scalar.return_value = 0.3679 - result = Exponential._pdf_scalar(1.0, 1.0) - mock_core.exponential_pdf_scalar.assert_called_once_with(1.0, 1.0) - assert result == 0.3679 +@pytest.mark.parametrize("lam", RATES) +def test_pdf_array_matches_scalar_evaluation(lam): + dist = Exponential(lam) + xs = [0.0, 0.25, 1.0, 3.0, 7.5] + result = dist.pdf(xs) + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + assert result.shape == (len(xs),) + np.testing.assert_allclose(result, [exp_pdf(x, lam) for x in xs], rtol=1e-12) + +@pytest.mark.parametrize("lam", RATES) +def test_cdf_array_matches_scalar_evaluation(lam): + dist = Exponential(lam) + xs = [0.0, 0.25, 1.0, 3.0, 7.5] + np.testing.assert_allclose( + dist.cdf(xs), [exp_cdf(x, lam) for x in xs], rtol=1e-12 + ) -def test_cdf_scalar_delegates_to_core(mock_core): - mock_core.exponential_cdf_scalar.return_value = 0.6321 - result = Exponential._cdf_scalar(1.0, 1.0) - mock_core.exponential_cdf_scalar.assert_called_once_with(1.0, 1.0) - assert result == 0.6321 +def test_array_accepts_numpy_input(): + dist = Exponential(2.0) + xs = np.array([0.5, 1.0, 2.0]) + np.testing.assert_allclose(dist.pdf(xs), [exp_pdf(x, 2.0) for x in xs], rtol=1e-12) -def test_mgf_scalar_delegates_to_core(mock_core): - mock_core.exponential_mgf_scalar.return_value = 2.5 - result = Exponential._mgf_scalar(0.5, 2.0) - mock_core.exponential_mgf_scalar.assert_called_once_with(0.5, 2.0) - assert result == 2.5 +def test_empty_array_returns_empty_array(): + result = Exponential(2.0).pdf([]) + assert isinstance(result, np.ndarray) + assert result.size == 0 -def test_cgf_scalar_delegates_to_core(mock_core): - mock_core.exponential_cgf_scalar.return_value = 0.916 - result = Exponential._cgf_scalar(0.5, 2.0) - mock_core.exponential_cgf_scalar.assert_called_once_with(0.5, 2.0) - assert result == 0.916 + +def test_step_size_offsets_each_element_by_its_index(): + """With step_size s, element i is evaluated at x[i] + s*i.""" + dist = Exponential(2.0) + result = dist.pdf([0.0, 0.0, 0.0], 0.5) + np.testing.assert_allclose( + result, [exp_pdf(0.0, 2.0), exp_pdf(0.5, 2.0), exp_pdf(1.0, 2.0)], rtol=1e-12 + ) + + +def test_step_size_zero_is_a_plain_evaluation(): + dist = Exponential(2.0) + np.testing.assert_allclose(dist.pdf([0.5, 1.0], 0), dist.pdf([0.5, 1.0]), rtol=1e-15) + + +def test_array_rejects_two_dimensional_input(): + with pytest.raises(ValueError, match="must be 1-dimensional"): + Exponential(2.0).pdf([[1.0, 2.0], [3.0, 4.0]]) + + +def test_array_rejects_non_numeric_input(): + with pytest.raises(TypeError, match="must be numeric"): + Exponential(2.0).pdf(["a", "b"]) + + +def test_rejects_none_input(): + with pytest.raises(TypeError, match="must not be None"): + Exponential(2.0).pdf(None) # --------------------------------------------------------------------------- -# Instance methods that rely on self.lambda_ +# Moments # --------------------------------------------------------------------------- -def test_instance_methods_delegation(mock_core): - mock_core.exponential_sample.return_value = 0.347 - mock_core.exponential_mean.return_value = 2.0 - mock_core.exponential_variance.return_value = 4.0 - mock_core.exponential_stddev.return_value = 2.0 +@pytest.mark.parametrize("lam", RATES) +def test_moments_match_closed_forms(lam): + dist = Exponential(lam) + assert dist.mean() == pytest.approx(1.0 / lam, **EXACT) + assert dist.variance() == pytest.approx(1.0 / lam ** 2, **EXACT) + assert dist.stddev() == pytest.approx(1.0 / lam, **EXACT) - e = Exponential(lambda_=2.0) - # sample - result = e.sample() - mock_core.exponential_sample.assert_called_once_with(2.0) - assert result == 0.347 +@pytest.mark.parametrize("lam, override", [(2.0, 4.0), (1.0, 0.5)]) +def test_moment_parameter_override(lam, override): + """mean/variance/stddev accept an explicit rate that overrides the instance.""" + dist = Exponential(lam) + assert dist.mean(override) == pytest.approx(1.0 / override, **EXACT) + assert dist.variance(override) == pytest.approx(1.0 / override ** 2, **EXACT) + assert dist.stddev(override) == pytest.approx(1.0 / override, **EXACT) + assert dist.lambda_ == lam # the override must not mutate the instance - # mean - result = e.mean() - mock_core.exponential_mean.assert_called_once_with(2.0) - assert result == 2.0 - # variance - result = e.variance() - mock_core.exponential_variance.assert_called_once_with(2.0) - assert result == 4.0 +@pytest.mark.parametrize("lam", RATES) +def test_moment_override_validates(lam): + with pytest.raises(ValueError, match="lambda_ must be positive"): + Exponential(lam).mean(-1.0) + + +# --------------------------------------------------------------------------- +# MGF / CGF +# --------------------------------------------------------------------------- - # stddev - result = e.stddev() - mock_core.exponential_stddev.assert_called_once_with(2.0) - assert result == 2.0 +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("t_frac", [-2.0, -0.5, 0.0, 0.25, 0.5]) +def test_mgf_matches_closed_form(t_frac, lam): + t = t_frac * lam # keep t strictly below lambda + assert Exponential(lam).mgf(t) == pytest.approx(exp_mgf(t, lam), **EXACT) + + +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("t_frac", [-2.0, -0.5, 0.0, 0.25, 0.5]) +def test_cgf_matches_closed_form(t_frac, lam): + t = t_frac * lam + assert Exponential(lam).cgf(t) == pytest.approx( + math.log(exp_mgf(t, lam)), **EXACT + ) + + +@pytest.mark.parametrize("lam", RATES) +def test_mgf_at_zero_is_one(lam): + dist = Exponential(lam) + assert dist.mgf(0.0) == pytest.approx(1.0, **EXACT) + assert dist.cgf(0.0) == pytest.approx(0.0, abs=1e-12) + + +@pytest.mark.parametrize("lam", [1.0, 2.0]) +def test_mgf_array_matches_scalar_evaluation(lam): + dist = Exponential(lam) + ts = [-1.0, -0.25, 0.0, 0.25 * lam] + np.testing.assert_allclose( + dist.mgf(ts), [exp_mgf(t, lam) for t in ts], rtol=1e-12 + ) # --------------------------------------------------------------------------- -# Reject invalid lambda_ in instance methods +# Classmethods # --------------------------------------------------------------------------- -@pytest.mark.parametrize( - "method,invalid_lambda", - [ - (Exponential.mean, 0), - (Exponential.variance, -0.5), - (Exponential.stddev, -2), - (Exponential.sample, 0), - (Exponential.sample, -1), - ] -) -def test_instance_methods_reject_invalid_lambda(method, invalid_lambda): - e = Exponential(lambda_=1.0) +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("x", [0.25, 1.0, 3.0]) +def test_classmethods_agree_with_instances(x, lam): + dist = Exponential(lam) + assert Exponential._pdf_scalar(x, lam) == pytest.approx(dist.pdf(x), **EXACT) + assert Exponential._cdf_scalar(x, lam) == pytest.approx(dist.cdf(x), **EXACT) + assert Exponential._mgf_scalar(0.1, lam) == pytest.approx(dist.mgf(0.1), **EXACT) + assert Exponential._cgf_scalar(0.1, lam) == pytest.approx(dist.cgf(0.1), **EXACT) + + +@pytest.mark.parametrize("lam", [1.0, 2.0]) +def test_cpu_batch_classmethods_match_scalars(lam): + xs = [0.0, 0.5, 1.0, 4.0] + np.testing.assert_allclose( + Exponential._pdf_cpu(xs, lam), [exp_pdf(x, lam) for x in xs], rtol=1e-12 + ) + np.testing.assert_allclose( + Exponential._cdf_cpu(xs, lam), [exp_cdf(x, lam) for x in xs], rtol=1e-12 + ) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pdf_scalar", (1.0, 0.0)), + ("_pdf_scalar", (1.0, -2.0)), + ("_cdf_scalar", (1.0, 0.0)), + ("_mgf_scalar", (0.1, -1.0)), + ("_cgf_scalar", (0.1, 0.0)), +]) +def test_classmethods_reject_non_positive_rate(method_name, args): with pytest.raises(ValueError, match="lambda_ must be positive"): - method(e, invalid_lambda) # pass invalid lambda_ explicitly + getattr(Exponential, method_name)(*args) + + +def test_is_cuda_available_returns_bool(): + assert isinstance(Exponential.is_cuda_available(), bool) + + +# --------------------------------------------------------------------------- +# Sampling +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("lam", RATES) +def test_sample_is_positive_and_finite(lam): + dist = Exponential(lam) + for _ in range(200): + value = dist.sample() + assert math.isfinite(value) + assert value >= 0.0 # --------------------------------------------------------------------------- -# Slots behavior +# Slots # --------------------------------------------------------------------------- def test_slots_prevent_dynamic_attributes(): - e = Exponential(lambda_=0.5) with pytest.raises(AttributeError): - e.extra = 123 + Exponential(2.0).extra = 123 diff --git a/python/tests/test_gamma.py b/python/tests/test_gamma.py index 4b7a3ee..9146036 100644 --- a/python/tests/test_gamma.py +++ b/python/tests/test_gamma.py @@ -1,216 +1,332 @@ -# tests/python_modules/test_gamma.py -from unittest.mock import patch +""" +Numeric tests for the Gamma distribution against closed-form references. + +`theta` is the *scale* parameter: mean = alpha * theta, variance = alpha * theta^2. +""" + +import math import pytest -import fastdist.distributions.gamma as gamma_module +from conftest import EXACT, ITERATIVE, regularized_lower_gamma from fastdist.distributions.gamma import Gamma -@pytest.fixture -def mock_core(): - """Patch the internal C++ core for the duration of each test.""" - with patch.object(gamma_module, "_core", create=True) as mock: - yield mock +# --------------------------------------------------------------------------- +# Closed-form references +# --------------------------------------------------------------------------- + +def gamma_pdf(x: float, alpha: float, theta: float) -> float: + """f(x; a, th) = x^(a-1) e^(-x/th) / (th^a Gamma(a))""" + return x ** (alpha - 1) * math.exp(-x / theta) / (theta ** alpha * math.gamma(alpha)) + + +def gamma_mgf(t: float, alpha: float, theta: float) -> float: + """M(t) = (1 - th t)^(-a), valid for t < 1/th""" + return (1 - theta * t) ** (-alpha) + + +def gamma_cgf(t: float, alpha: float, theta: float) -> float: + """K(t) = -a ln(1 - th t)""" + return -alpha * math.log(1 - theta * t) + + +PARAMS = [(1.0, 1.0), (2.0, 3.0), (5.0, 2.0), (0.5, 4.0), (3.0, 0.5)] # --------------------------------------------------------------------------- -# Constructor & representation +# Constructor, properties, representation # --------------------------------------------------------------------------- def test_init_valid_parameters(): - g = Gamma(alpha=2.0, theta=3.0) - assert g.alpha == 2.0 - assert g.theta == 3.0 - - -@pytest.mark.parametrize("alpha,theta", [ - (-0.1, 1.0), - (0, 1.0), - (-1, 1.0), - (1.0, -0.1), - (1.0, 0), - (1.0, -1), - (-1, -1), -]) -def test_init_invalid_parameters_raises(alpha, theta): + dist = Gamma(alpha=2.0, theta=3.0) + assert dist.alpha == 2.0 + assert dist.theta == 3.0 + + +@pytest.mark.parametrize("alpha, theta", [(0, 1.0), (-1, 1.0), (1.0, 0), (1.0, -2.0)]) +def test_init_rejects_non_positive_parameters(alpha, theta): with pytest.raises(ValueError, match="must be positive"): Gamma(alpha=alpha, theta=theta) -def test_repr(): - g = Gamma(alpha=2.0, theta=3.0) - assert repr(g) == "Gamma(alpha=2.0, theta=3.0)" +@pytest.mark.parametrize("bad", ["2.0", object()]) +def test_init_rejects_non_real_parameters(bad): + with pytest.raises(TypeError, match="must be a real number"): + Gamma(alpha=bad, theta=1.0) -# --------------------------------------------------------------------------- -# Validation behavior -# --------------------------------------------------------------------------- +def test_init_rejects_none_parameter(): + """ + None is rejected, but only incidentally: _validate_params skips its checks + when a parameter is None, so the failure surfaces later as + "float() argument must be a string or a real number" rather than the + intended "alpha must be a real number". + """ + with pytest.raises(TypeError): + Gamma(alpha=None, theta=1.0) -@pytest.mark.parametrize("alpha,theta", [ - (0.1, 0.1), - (1.0, 1.0), - (2.5, 3.5), - (100, 200), -]) -def test_validate_params_accepts_valid_values(alpha, theta): - Gamma._validate_params(alpha=alpha, theta=theta) +def test_repr(): + assert repr(Gamma(alpha=2.0, theta=3.0)) == "Gamma(alpha=2.0, theta=3.0)" -@pytest.mark.parametrize("alpha,theta", [ - (0, 1.0), - (-1, 1.0), - (1.0, 0), - (1.0, -1), -]) -def test_validate_params_rejects_invalid_values(alpha, theta): - with pytest.raises(ValueError, match="must be positive"): - Gamma._validate_params(alpha=alpha, theta=theta) + +def test_property_setters_update_and_validate(): + dist = Gamma(alpha=2.0, theta=3.0) + dist.alpha = 4.0 + dist.theta = 5.0 + assert dist.alpha == 4.0 + assert dist.theta == 5.0 + with pytest.raises(ValueError, match="alpha must be positive"): + dist.alpha = 0 + with pytest.raises(ValueError, match="theta must be positive"): + dist.theta = -1.0 # --------------------------------------------------------------------------- -# Instance method delegation +# PDF +# +# NOTE: the instance method is named `pmf_scalar` even though Gamma is a +# continuous distribution and the corresponding classmethod is `_pdf_scalar`. +# The tests follow the shipped API; the naming inconsistency is tracked +# separately. # --------------------------------------------------------------------------- -def test_pmf_scalar_delegates_to_core(mock_core): - mock_core.gamma_pdf_scalar.return_value = 0.234 +@pytest.mark.parametrize("alpha, theta", PARAMS) +@pytest.mark.parametrize("x", [0.5, 1.0, 3.0, 8.0]) +def test_pdf_matches_closed_form(x, alpha, theta): + assert Gamma(alpha, theta).pmf_scalar(x) == pytest.approx( + gamma_pdf(x, alpha, theta), **EXACT + ) - g = Gamma(alpha=2.0, theta=3.0) - result = g.pmf_scalar(5.0) - mock_core.gamma_pdf_scalar.assert_called_once_with(5.0, 2.0, 3.0) - assert result == 0.234 +@pytest.mark.parametrize("theta", [0.5, 1.0, 3.0]) +@pytest.mark.parametrize("x", [0.25, 1.0, 4.0]) +def test_pdf_with_unit_shape_is_exponential(x, theta): + """Gamma(1, th) is Exponential with rate 1/th.""" + assert Gamma(1.0, theta).pmf_scalar(x) == pytest.approx( + math.exp(-x / theta) / theta, **EXACT + ) -def test_cdf_scalar_delegates_to_core(mock_core): - mock_core.gamma_cdf_scalar.return_value = 0.456 +@pytest.mark.parametrize("alpha, theta", PARAMS) +@pytest.mark.parametrize("x", [0.1, 2.0, 15.0]) +def test_pdf_is_non_negative(x, alpha, theta): + assert Gamma(alpha, theta).pmf_scalar(x) >= 0.0 - g = Gamma(alpha=2.0, theta=3.0) - result = g.cdf_scalar(5.0) - mock_core.gamma_cdf_scalar.assert_called_once_with(5.0, 2.0, 3.0) - assert result == 0.456 +# --------------------------------------------------------------------------- +# CDF +# +# KNOWN BUG: the regularized lower incomplete gamma used by gamma_cdf_scalar is +# correct in its series branch (x/theta < alpha+1) but wrong in the continued- +# fraction branch. Absolute errors reach 0.26, and the CDF can exceed 1.0 -- +# Gamma(1.5, 1.0).cdf_scalar(2.5) returns 1.000498004. +# +# The failing points below were determined empirically by comparing against +# conftest.regularized_lower_gamma. They are marked strict-xfail so they become +# failures the moment the backend is fixed and the markers go stale. Note that +# alpha = 1 is correct throughout, so the failure set is not simply +# "everything in the continued-fraction branch". +# --------------------------------------------------------------------------- +CF_BUG = ( + "regularized lower incomplete gamma is incorrect in the " + "continued-fraction branch (x/theta >= alpha+1)" +) -def test_mean_delegates_to_core(mock_core): - mock_core.gamma_mean.return_value = 6.0 +GAMMA_SHAPES = [(0.5, 1.0), (1.0, 1.0), (1.5, 1.0), (2.0, 3.0), (3.0, 1.0), (5.0, 2.0)] +GAMMA_X = (0.2, 0.8, 1.5, 2.5, 4.0, 10.0, 20.0) - g = Gamma(alpha=2.0, theta=3.0) - result = g.mean() +_CDF_KNOWN_BAD = { + (0.5, 1.0, 1.5), (0.5, 1.0, 2.5), (0.5, 1.0, 4.0), + (0.5, 1.0, 10.0), (0.5, 1.0, 20.0), + (1.5, 1.0, 2.5), (1.5, 1.0, 4.0), (1.5, 1.0, 10.0), (1.5, 1.0, 20.0), + (2.0, 3.0, 10.0), (2.0, 3.0, 20.0), + (3.0, 1.0, 4.0), (3.0, 1.0, 10.0), (3.0, 1.0, 20.0), + (5.0, 2.0, 20.0), +} - mock_core.gamma_mean.assert_called_once_with(2.0, 3.0) - assert result == 6.0 +def _cdf_case(alpha, theta, x): + marks = [] + if (alpha, theta, x) in _CDF_KNOWN_BAD: + marks = [ + pytest.mark.known_bug, + pytest.mark.xfail(strict=True, reason=CF_BUG), + ] + return pytest.param(alpha, theta, x, marks=marks) -def test_variance_delegates_to_core(mock_core): - mock_core.gamma_variance.return_value = 18.0 - g = Gamma(alpha=2.0, theta=3.0) - result = g.variance() +CDF_CASES = [ + _cdf_case(alpha, theta, x) + for (alpha, theta) in GAMMA_SHAPES + for x in GAMMA_X +] - mock_core.gamma_variance.assert_called_once_with(2.0, 3.0) - assert result == 18.0 +@pytest.mark.parametrize("alpha, theta, x", CDF_CASES) +def test_cdf_matches_reference(alpha, theta, x): + assert Gamma(alpha, theta).cdf_scalar(x) == pytest.approx( + regularized_lower_gamma(alpha, x / theta), **ITERATIVE + ) -def test_stddev_delegates_to_core(mock_core): - mock_core.gamma_stddev.return_value = 4.243 - g = Gamma(alpha=2.0, theta=3.0) - result = g.stddev() +@pytest.mark.parametrize("theta", [0.5, 1.0, 3.0]) +@pytest.mark.parametrize("x", [0.25, 1.0, 2.0, 6.0]) +def test_cdf_with_unit_shape_is_exponential(x, theta): + """Gamma(1, th) has the exact CDF 1 - e^(-x/th).""" + assert Gamma(1.0, theta).cdf_scalar(x) == pytest.approx( + 1 - math.exp(-x / theta), **ITERATIVE + ) - mock_core.gamma_stddev.assert_called_once_with(2.0, 3.0) - assert result == 4.243 +_BOUNDED_SHAPES = [ + (0.5, 1.0), + (1.0, 1.0), + pytest.param(1.5, 1.0, marks=[ + pytest.mark.known_bug, + pytest.mark.xfail(strict=True, reason=CF_BUG + "; CDF exceeds 1.0"), + ]), + pytest.param(2.0, 3.0, marks=[ + pytest.mark.known_bug, + pytest.mark.xfail(strict=True, reason=CF_BUG + "; CDF exceeds 1.0"), + ]), + (3.0, 1.0), + (5.0, 2.0), +] + + +@pytest.mark.parametrize("alpha, theta", _BOUNDED_SHAPES) +def test_cdf_is_bounded(alpha, theta): + dist = Gamma(alpha, theta) + for x in (0.1, 0.5, 1.0, 3.0, 8.0, 25.0): + assert 0.0 <= dist.cdf_scalar(x) <= 1.0 + + +_MONOTONIC_SHAPES = [ + (0.5, 1.0), + (1.0, 1.0), + pytest.param(1.5, 1.0, marks=[ + pytest.mark.known_bug, + pytest.mark.xfail(strict=True, reason=CF_BUG + "; CDF is non-monotonic"), + ]), + (2.0, 3.0), + (3.0, 1.0), + (5.0, 2.0), +] -def test_mgf_scalar_delegates_to_core(mock_core): - mock_core.gamma_mgf_scalar.return_value = 1.789 - g = Gamma(alpha=2.0, theta=3.0) - result = g.mgf_scalar(0.1) +@pytest.mark.parametrize("alpha, theta", _MONOTONIC_SHAPES) +def test_cdf_is_monotonic(alpha, theta): + dist = Gamma(alpha, theta) + values = [dist.cdf_scalar(x) for x in (0.1, 0.5, 1.0, 3.0, 8.0, 25.0)] + assert values == sorted(values) - mock_core.gamma_mgf_scalar.assert_called_once_with(0.1, 2.0, 3.0) - assert result == 1.789 +# --------------------------------------------------------------------------- +# Moments +# --------------------------------------------------------------------------- -def test_cgf_scalar_delegates_to_core(mock_core): - mock_core.gamma_cgf_scalar.return_value = 0.581 +@pytest.mark.parametrize("alpha, theta", PARAMS) +def test_moments_match_closed_forms(alpha, theta): + dist = Gamma(alpha, theta) + assert dist.mean() == pytest.approx(alpha * theta, **EXACT) + assert dist.variance() == pytest.approx(alpha * theta ** 2, **EXACT) + assert dist.stddev() == pytest.approx(math.sqrt(alpha) * theta, **EXACT) - g = Gamma(alpha=2.0, theta=3.0) - result = g.cgf_scalar(0.1) - mock_core.gamma_cgf_scalar.assert_called_once_with(0.1, 2.0, 3.0) - assert result == 0.581 +# --------------------------------------------------------------------------- +# MGF / CGF +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("alpha, theta", [(2.0, 3.0), (5.0, 2.0), (1.0, 1.0)]) +@pytest.mark.parametrize("t_over_theta", [-1.0, -0.25, 0.0, 0.1, 0.25]) +def test_mgf_matches_closed_form(t_over_theta, alpha, theta): + t = t_over_theta / theta # keep t strictly inside the radius of convergence + assert Gamma(alpha, theta).mgf_scalar(t) == pytest.approx( + gamma_mgf(t, alpha, theta), **EXACT + ) + + +@pytest.mark.parametrize("alpha, theta", [(2.0, 3.0), (5.0, 2.0), (1.0, 1.0)]) +@pytest.mark.parametrize("t_over_theta", [-1.0, -0.25, 0.0, 0.1, 0.25]) +def test_cgf_matches_closed_form(t_over_theta, alpha, theta): + t = t_over_theta / theta + assert Gamma(alpha, theta).cgf_scalar(t) == pytest.approx( + gamma_cgf(t, alpha, theta), **EXACT + ) -def test_sample_delegates_to_core(mock_core): - mock_core.gamma_sample.return_value = 5.823 +@pytest.mark.parametrize("alpha, theta", PARAMS) +def test_mgf_at_zero_is_one(alpha, theta): + dist = Gamma(alpha, theta) + assert dist.mgf_scalar(0.0) == pytest.approx(1.0, **EXACT) + assert dist.cgf_scalar(0.0) == pytest.approx(0.0, abs=1e-12) - g = Gamma(alpha=2.0, theta=3.0) - result = g.sample() - mock_core.gamma_sample.assert_called_once_with(2.0, 3.0) - assert result == 5.823 +@pytest.mark.parametrize("alpha, theta", [(2.0, 3.0), (5.0, 2.0)]) +def test_cgf_is_log_of_mgf(alpha, theta): + dist = Gamma(alpha, theta) + for t in (-0.5 / theta, 0.1 / theta, 0.25 / theta): + assert dist.cgf_scalar(t) == pytest.approx( + math.log(dist.mgf_scalar(t)), **EXACT + ) # --------------------------------------------------------------------------- -# Validation enforcement in classmethods +# Classmethods # --------------------------------------------------------------------------- -@pytest.mark.parametrize( - "method, args", - [ - (Gamma._pdf_scalar, (1.0, -1, 1.0)), - (Gamma._pdf_scalar, (1.0, 0, 1.0)), - (Gamma._pdf_scalar, (1.0, 1.0, -1)), - (Gamma._pdf_scalar, (1.0, 1.0, 0)), - (Gamma._cdf_scalar, (1.0, -5.0, 1.0)), - (Gamma._cdf_scalar, (1.0, 1.0, -5.0)), - (Gamma._mean, (-1, 1.0)), - (Gamma._mean, (1.0, -1)), - (Gamma._variance, (0, 1.0)), - (Gamma._variance, (1.0, 0)), - (Gamma._stddev, (-1, 1.0)), - (Gamma._stddev, (1.0, -1)), - (Gamma._mgf_scalar, (0.1, -1, 1.0)), - (Gamma._mgf_scalar, (0.1, 1.0, -1)), - (Gamma._cgf_scalar, (0.1, -5.0, 1.0)), - (Gamma._cgf_scalar, (0.1, 1.0, -5.0)), - (Gamma._sample, (-1, 1.0)), - (Gamma._sample, (1.0, -1)), - ], -) -def test_classmethods_reject_invalid_parameters(method, args): - with pytest.raises(ValueError, match=r"must be positive"): - method(*args) +@pytest.mark.parametrize("alpha, theta", PARAMS) +@pytest.mark.parametrize("x", [0.5, 3.0]) +def test_classmethods_agree_with_instances(x, alpha, theta): + dist = Gamma(alpha, theta) + assert Gamma._pdf_scalar(x, alpha, theta) == pytest.approx( + dist.pmf_scalar(x), **EXACT + ) + assert Gamma._cdf_scalar(x, alpha, theta) == pytest.approx( + dist.cdf_scalar(x), **ITERATIVE + ) + assert Gamma._mean(alpha, theta) == pytest.approx(dist.mean(), **EXACT) + assert Gamma._variance(alpha, theta) == pytest.approx(dist.variance(), **EXACT) + assert Gamma._stddev(alpha, theta) == pytest.approx(dist.stddev(), **EXACT) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pdf_scalar", (1.0, 0, 1.0)), + ("_pdf_scalar", (1.0, 1.0, 0)), + ("_cdf_scalar", (1.0, -1.0, 1.0)), + ("_cdf_scalar", (1.0, 1.0, -1.0)), + ("_mean", (0, 1.0)), + ("_variance", (1.0, 0)), + ("_stddev", (-2.0, 1.0)), + ("_mgf_scalar", (0.1, 0, 1.0)), + ("_cgf_scalar", (0.1, 1.0, -3.0)), + ("_sample", (0, 1.0)), +]) +def test_classmethods_reject_invalid_parameters(method_name, args): + with pytest.raises(ValueError, match="must be positive"): + getattr(Gamma, method_name)(*args) # --------------------------------------------------------------------------- -# Classmethod delegation with valid parameters +# Sampling # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, args, value", [ - ("_pdf_scalar", "gamma_pdf_scalar", (5.0, 2.0, 3.0), 0.234), - ("_cdf_scalar", "gamma_cdf_scalar", (5.0, 2.0, 3.0), 0.456), - ("_mean", "gamma_mean", (2.0, 3.0), 6.0), - ("_variance", "gamma_variance", (2.0, 3.0), 18.0), - ("_stddev", "gamma_stddev", (2.0, 3.0), 4.243), - ("_mgf_scalar", "gamma_mgf_scalar", (0.1, 2.0, 3.0), 1.789), - ("_cgf_scalar", "gamma_cgf_scalar", (0.1, 2.0, 3.0), 0.581), - ("_sample", "gamma_sample", (2.0, 3.0), 5.823), -]) -def test_classmethods_delegate_to_core(mock_core, method_name, core_method_name, args, value): - getattr(mock_core, core_method_name).return_value = value - method = getattr(Gamma, method_name) - result = method(*args) - getattr(mock_core, core_method_name).assert_called_once_with(*args) - assert result == value +@pytest.mark.parametrize("alpha, theta", [(1.0, 1.0), (2.0, 3.0), (5.0, 2.0)]) +def test_sample_is_positive_and_finite(alpha, theta): + dist = Gamma(alpha, theta) + for _ in range(200): + value = dist.sample() + assert math.isfinite(value) + assert value > 0.0 # --------------------------------------------------------------------------- -# Slots behavior +# Slots # --------------------------------------------------------------------------- def test_slots_prevent_dynamic_attributes(): - g = Gamma(alpha=2.0, theta=3.0) with pytest.raises(AttributeError): - g.extra = 123 + Gamma(2.0, 3.0).extra = 123 diff --git a/python/tests/test_geometric.py b/python/tests/test_geometric.py index 066e7df..060f3b6 100644 --- a/python/tests/test_geometric.py +++ b/python/tests/test_geometric.py @@ -1,198 +1,248 @@ -# tests/python_modules/test_geometric.py -from unittest.mock import patch +""" +Numeric tests for the Geometric distribution against closed-form references. + +This implementation uses the "number of trials until the first success" +convention: the support is k = 1, 2, 3, ... and P(X = k) = (1-p)^(k-1) p, +so the mean is 1/p. +""" + +import math import pytest -import fastdist.distributions.geometric as geometric_module +from conftest import EXACT, ITERATIVE from fastdist.distributions.geometric import Geometric -@pytest.fixture -def mock_core(): - """Patch the internal C++ core for the duration of each test.""" - with patch.object(geometric_module, "_core", create=True) as mock: - yield mock +# --------------------------------------------------------------------------- +# Closed-form references +# --------------------------------------------------------------------------- + +def geom_pmf(k: int, p: float) -> float: + """P(X = k) = (1-p)^(k-1) p for k >= 1""" + return (1 - p) ** (k - 1) * p + + +def geom_cdf(k: int, p: float) -> float: + """F(k) = 1 - (1-p)^k for k >= 1""" + return 1.0 - (1 - p) ** k + + +def geom_mgf(t: float, p: float) -> float: + """M(t) = p e^t / (1 - (1-p) e^t), valid for t < -ln(1-p)""" + return p * math.exp(t) / (1 - (1 - p) * math.exp(t)) + + +PROBS = [0.1, 0.3, 0.5, 0.75, 1.0] # --------------------------------------------------------------------------- -# Constructor & representation +# Constructor, properties, representation # --------------------------------------------------------------------------- def test_init_valid_parameter(): - g = Geometric(p=0.5) - assert g.p == 0.5 + assert Geometric(p=0.3).p == 0.3 -@pytest.mark.parametrize("p", [0, -0.1, -1, 1.1, 2]) -def test_init_invalid_p_raises(p): +@pytest.mark.parametrize("p", [0.0, -0.1, 1.1, 2.0]) +def test_init_rejects_out_of_range_p(p): with pytest.raises(ValueError, match=r"p must be in the interval \(0, 1\]"): Geometric(p=p) -def test_repr(): - g = Geometric(p=0.3) - assert repr(g) == "Geometric(p=0.3)" +@pytest.mark.parametrize("bad", ["0.5", None, object()]) +def test_init_rejects_non_real_p(bad): + with pytest.raises(TypeError, match="p must be a real number"): + Geometric(p=bad) -# --------------------------------------------------------------------------- -# Validation behavior -# --------------------------------------------------------------------------- - -@pytest.mark.parametrize("p", [0.001, 0.5, 1.0]) -def test_validate_params_accepts_valid_values(p): - Geometric._validate_params(p=p) +def test_repr(): + assert repr(Geometric(p=0.3)) == "Geometric(p=0.3)" -@pytest.mark.parametrize("p", [0, -0.5, 1.1, 2]) -def test_validate_params_rejects_invalid_values(p): +def test_p_setter_updates_and_validates(): + dist = Geometric(p=0.3) + dist.p = 0.8 + assert dist.p == 0.8 with pytest.raises(ValueError, match=r"p must be in the interval \(0, 1\]"): - Geometric._validate_params(p=p) + dist.p = 0.0 # --------------------------------------------------------------------------- -# Instance method delegation +# PMF # --------------------------------------------------------------------------- -def test_pmf_scalar_delegates_to_core(mock_core): - mock_core.geometric_pmf_scalar.return_value = 0.125 +@pytest.mark.parametrize("p", PROBS) +@pytest.mark.parametrize("k", [1, 2, 3, 5, 10]) +def test_pmf_matches_closed_form(k, p): + assert Geometric(p).pmf_scalar(k) == pytest.approx(geom_pmf(k, p), **EXACT) - g = Geometric(p=0.5) - result = g.pmf_scalar(3) - mock_core.geometric_pmf_scalar.assert_called_once_with(3, 0.5) - assert result == 0.125 +@pytest.mark.parametrize("p", PROBS) +def test_pmf_is_zero_below_the_support(p): + """The support starts at k = 1, so k = 0 has zero mass.""" + assert Geometric(p).pmf_scalar(0) == pytest.approx(0.0, abs=1e-15) -def test_cdf_scalar_delegates_to_core(mock_core): - mock_core.geometric_cdf_scalar.return_value = 0.875 +@pytest.mark.parametrize("p", [0.1, 0.3, 0.5, 0.75]) +def test_pmf_sums_to_one(p): + dist = Geometric(p) + total = sum(dist.pmf_scalar(k) for k in range(1, 2000)) + assert total == pytest.approx(1.0, **ITERATIVE) - g = Geometric(p=0.5) - result = g.cdf_scalar(3) - mock_core.geometric_cdf_scalar.assert_called_once_with(3, 0.5) - assert result == 0.875 +@pytest.mark.parametrize("p", PROBS) +@pytest.mark.parametrize("k", [1, 4, 9]) +def test_pmf_is_a_probability(k, p): + assert 0.0 <= Geometric(p).pmf_scalar(k) <= 1.0 -def test_mean_delegates_to_core(mock_core): - mock_core.geometric_mean.return_value = 2.0 +@pytest.mark.parametrize("p", [0.1, 0.3, 0.5]) +def test_pmf_is_decreasing(p): + dist = Geometric(p) + values = [dist.pmf_scalar(k) for k in range(1, 12)] + assert values == sorted(values, reverse=True) - g = Geometric(p=0.5) - result = g.mean() - - mock_core.geometric_mean.assert_called_once_with(0.5) - assert result == 2.0 +# --------------------------------------------------------------------------- +# CDF +# --------------------------------------------------------------------------- -def test_variance_delegates_to_core(mock_core): - mock_core.geometric_variance.return_value = 2.0 +@pytest.mark.parametrize("p", PROBS) +@pytest.mark.parametrize("k", [1, 2, 3, 5, 10]) +def test_cdf_matches_closed_form(k, p): + assert Geometric(p).cdf_scalar(k) == pytest.approx(geom_cdf(k, p), **ITERATIVE) - g = Geometric(p=0.5) - result = g.variance() - mock_core.geometric_variance.assert_called_once_with(0.5) - assert result == 2.0 +@pytest.mark.parametrize("p", [0.1, 0.3, 0.5, 0.75]) +@pytest.mark.parametrize("k", [1, 3, 7]) +def test_cdf_matches_summed_pmf(k, p): + dist = Geometric(p) + assert dist.cdf_scalar(k) == pytest.approx( + sum(dist.pmf_scalar(i) for i in range(1, k + 1)), **ITERATIVE + ) -def test_stddev_delegates_to_core(mock_core): - mock_core.geometric_stddev.return_value = 1.414 +@pytest.mark.parametrize("p", PROBS) +def test_cdf_is_monotonic_and_bounded(p): + dist = Geometric(p) + values = [dist.cdf_scalar(k) for k in range(1, 20)] + assert values == sorted(values) + assert all(0.0 <= v <= 1.0 for v in values) - g = Geometric(p=0.5) - result = g.stddev() - mock_core.geometric_stddev.assert_called_once_with(0.5) - assert result == 1.414 +# --------------------------------------------------------------------------- +# Moments +# --------------------------------------------------------------------------- +@pytest.mark.parametrize("p", PROBS) +def test_moments_match_closed_forms(p): + dist = Geometric(p) + assert dist.mean() == pytest.approx(1.0 / p, **EXACT) + assert dist.variance() == pytest.approx((1 - p) / p ** 2, **EXACT) + assert dist.stddev() == pytest.approx(math.sqrt((1 - p) / p ** 2), **EXACT) -def test_mgf_scalar_delegates_to_core(mock_core): - mock_core.geometric_mgf_scalar.return_value = 1.648 - g = Geometric(p=0.5) - result = g.mgf_scalar(0.2) +@pytest.mark.parametrize("p", [0.3, 0.5, 0.75]) +def test_mean_equals_pmf_weighted_sum(p): + dist = Geometric(p) + expectation = sum(k * dist.pmf_scalar(k) for k in range(1, 5000)) + assert dist.mean() == pytest.approx(expectation, rel=1e-8) - mock_core.geometric_mgf_scalar.assert_called_once_with(0.2, 0.5) - assert result == 1.648 +# --------------------------------------------------------------------------- +# MGF / CGF +# --------------------------------------------------------------------------- -def test_cgf_scalar_delegates_to_core(mock_core): - mock_core.geometric_cgf_scalar.return_value = 0.499 +@pytest.mark.parametrize("p", [0.3, 0.5, 0.75]) +@pytest.mark.parametrize("t", [-1.0, -0.2, 0.0, 0.1]) +def test_mgf_matches_closed_form(t, p): + assert Geometric(p).mgf_scalar(t) == pytest.approx(geom_mgf(t, p), **EXACT) - g = Geometric(p=0.5) - result = g.cgf_scalar(0.2) - mock_core.geometric_cgf_scalar.assert_called_once_with(0.2, 0.5) - assert result == 0.499 +@pytest.mark.parametrize("p", [0.3, 0.5, 0.75]) +@pytest.mark.parametrize("t", [-1.0, -0.2, 0.0, 0.1]) +def test_cgf_matches_closed_form(t, p): + assert Geometric(p).cgf_scalar(t) == pytest.approx( + math.log(geom_mgf(t, p)), **EXACT + ) -def test_sample_delegates_to_core(mock_core): - mock_core.geometric_sample.return_value = 3 +@pytest.mark.parametrize("p", [0.3, 0.5, 0.75]) +def test_mgf_at_zero_is_one(p): + dist = Geometric(p) + assert dist.mgf_scalar(0.0) == pytest.approx(1.0, **EXACT) + assert dist.cgf_scalar(0.0) == pytest.approx(0.0, abs=1e-12) - g = Geometric(p=0.5) - result = g.sample() - mock_core.geometric_sample.assert_called_once_with(0.5) - assert result == 3 +@pytest.mark.parametrize("p", [0.3, 0.5]) +def test_cgf_is_log_of_mgf(p): + dist = Geometric(p) + for t in (-0.5, -0.1, 0.05): + assert dist.cgf_scalar(t) == pytest.approx( + math.log(dist.mgf_scalar(t)), **EXACT + ) # --------------------------------------------------------------------------- -# Validation enforcement in classmethods +# Classmethods # --------------------------------------------------------------------------- -@pytest.mark.parametrize( - "method, args", - [ - (Geometric._pmf_scalar, (3, 0)), - (Geometric._pmf_scalar, (3, -0.1)), - (Geometric._pmf_scalar, (3, 1.1)), - (Geometric._cdf_scalar, (3, 0)), - (Geometric._cdf_scalar, (3, -0.5)), - (Geometric._cdf_scalar, (3, 2)), - (Geometric._mean, (0,)), - (Geometric._mean, (-0.1,)), - (Geometric._mean, (1.5,)), - (Geometric._variance, (0,)), - (Geometric._variance, (-1,)), - (Geometric._stddev, (0,)), - (Geometric._stddev, (1.1,)), - (Geometric._mgf_scalar, (0.2, 0)), - (Geometric._mgf_scalar, (0.2, -0.1)), - (Geometric._cgf_scalar, (0.2, 0)), - (Geometric._cgf_scalar, (0.2, 1.5)), - (Geometric._sample, (0,)), - (Geometric._sample, (-0.5,)), - ], -) -def test_classmethods_reject_invalid_parameters(method, args): +@pytest.mark.parametrize("p", [0.1, 0.5, 0.75]) +def test_classmethods_agree_with_instances(p): + dist = Geometric(p) + assert Geometric._pmf_scalar(3, p) == pytest.approx(dist.pmf_scalar(3), **EXACT) + assert Geometric._cdf_scalar(3, p) == pytest.approx( + dist.cdf_scalar(3), **ITERATIVE + ) + assert Geometric._mean(p) == pytest.approx(dist.mean(), **EXACT) + assert Geometric._variance(p) == pytest.approx(dist.variance(), **EXACT) + assert Geometric._stddev(p) == pytest.approx(dist.stddev(), **EXACT) + assert Geometric._mgf_scalar(0.1, p) == pytest.approx( + dist.mgf_scalar(0.1), **EXACT + ) + assert Geometric._cgf_scalar(0.1, p) == pytest.approx( + dist.cgf_scalar(0.1), **EXACT + ) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pmf_scalar", (1, 0.0)), + ("_pmf_scalar", (1, 1.5)), + ("_cdf_scalar", (1, -0.2)), + ("_mean", (0.0,)), + ("_variance", (1.5,)), + ("_stddev", (-1.0,)), + ("_mgf_scalar", (0.1, 0.0)), + ("_cgf_scalar", (0.1, 2.0)), + ("_sample", (0.0,)), +]) +def test_classmethods_reject_out_of_range_p(method_name, args): with pytest.raises(ValueError, match=r"p must be in the interval \(0, 1\]"): - method(*args) + getattr(Geometric, method_name)(*args) # --------------------------------------------------------------------------- -# Classmethod delegation with valid parameters +# Sampling # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, args, value", [ - ("_pmf_scalar", "geometric_pmf_scalar", (3, 0.5), 0.125), - ("_cdf_scalar", "geometric_cdf_scalar", (3, 0.5), 0.875), - ("_mean", "geometric_mean", (0.5,), 2.0), - ("_variance", "geometric_variance", (0.5,), 2.0), - ("_stddev", "geometric_stddev", (0.5,), 1.414), - ("_mgf_scalar", "geometric_mgf_scalar", (0.2, 0.5), 1.648), - ("_cgf_scalar", "geometric_cgf_scalar", (0.2, 0.5), 0.499), - ("_sample", "geometric_sample", (0.5,), 3), -]) -def test_classmethods_delegate_to_core(mock_core, method_name, core_method_name, args, value): - getattr(mock_core, core_method_name).return_value = value - method = getattr(Geometric, method_name) - result = method(*args) - getattr(mock_core, core_method_name).assert_called_once_with(*args) - assert result == value +@pytest.mark.parametrize("p", [0.1, 0.5, 0.9]) +def test_sample_lies_within_the_support(p): + dist = Geometric(p) + for _ in range(300): + value = dist.sample() + assert isinstance(value, int) + assert value >= 1 + + +def test_sample_is_deterministic_for_certain_success(): + assert all(Geometric(1.0).sample() == 1 for _ in range(50)) # --------------------------------------------------------------------------- -# Slots behavior +# Slots # --------------------------------------------------------------------------- def test_slots_prevent_dynamic_attributes(): - g = Geometric(p=0.5) with pytest.raises(AttributeError): - g.extra = 123 \ No newline at end of file + Geometric(0.3).extra = 123 diff --git a/python/tests/test_negative_binomial.py b/python/tests/test_negative_binomial.py index 585f1b2..8279e34 100644 --- a/python/tests/test_negative_binomial.py +++ b/python/tests/test_negative_binomial.py @@ -1,219 +1,322 @@ -# tests/python_modules/test_negative_binomial.py -from unittest.mock import patch +""" +Numeric tests for the Negative Binomial distribution against closed-form +references. + +This implementation uses the "number of failures before the r-th success" +convention: the support is k = 0, 1, 2, ... with +P(X = k) = C(k + r - 1, k) p^r (1-p)^k, so the mean is r(1-p)/p. +""" + +import math import pytest -import fastdist.distributions.negative_binomial as nb_module +from conftest import EXACT, ITERATIVE from fastdist.distributions.negative_binomial import NegativeBinomial -@pytest.fixture -def mock_core(): - """Patch the internal C++ core for the duration of each test.""" - with patch.object(nb_module, "_core", create=True) as mock: - yield mock +# --------------------------------------------------------------------------- +# Closed-form references +# --------------------------------------------------------------------------- + +def nbinom_pmf(k: int, r: int, p: float) -> float: + """P(X = k) = C(k + r - 1, k) p^r (1-p)^k for k >= 0""" + return math.comb(k + r - 1, k) * p ** r * (1 - p) ** k + + +def nbinom_cdf(k: int, r: int, p: float) -> float: + return sum(nbinom_pmf(i, r, p) for i in range(0, k + 1)) + + +def nbinom_mgf(t: float, r: int, p: float) -> float: + """M(t) = (p / (1 - (1-p) e^t))^r, valid for t < -ln(1-p)""" + return (p / (1 - (1 - p) * math.exp(t))) ** r + + +PARAMS = [(1, 0.5), (3, 0.5), (2, 0.3), (5, 0.7), (4, 0.9)] # --------------------------------------------------------------------------- -# Constructor & representation +# Constructor, properties, representation # --------------------------------------------------------------------------- def test_init_valid_parameters(): - nb = NegativeBinomial(r=5, p=0.6) - assert nb.r == 5 - assert nb.p == 0.6 - - -@pytest.mark.parametrize("r,p", [ - (-1, 0.5), - (0, 0.5), - (-5, 0.5), - (5, -0.1), - (5, 1.1), - (5, 2), - (0, -0.5), -]) -def test_init_invalid_parameters_raises(r, p): - with pytest.raises(ValueError): - NegativeBinomial(r=r, p=p) + dist = NegativeBinomial(r=3, p=0.5) + assert dist.r == 3 + assert dist.p == 0.5 + + +@pytest.mark.parametrize("r", [0, -1, -10]) +def test_init_rejects_non_positive_r(r): + with pytest.raises(ValueError, match="r must be positive"): + NegativeBinomial(r=r, p=0.5) + + +@pytest.mark.parametrize("p", [-0.1, 1.1, 2.0]) +def test_init_rejects_out_of_range_p(p): + with pytest.raises(ValueError, match=r"p must be in \[0, 1\]"): + NegativeBinomial(r=3, p=p) + + +@pytest.mark.parametrize("bad_r", [2.5, "3"]) +def test_init_rejects_non_integer_r(bad_r): + with pytest.raises(TypeError, match="r must be an integer"): + NegativeBinomial(r=bad_r, p=0.5) + + +@pytest.mark.parametrize("bad_p", ["0.5", object()]) +def test_init_rejects_non_real_p(bad_p): + with pytest.raises(TypeError, match="p must be a real number"): + NegativeBinomial(r=3, p=bad_p) def test_repr(): - nb = NegativeBinomial(r=3, p=0.7) - assert repr(nb) == "NegativeBinomial(r=3, p=0.7)" + assert repr(NegativeBinomial(r=3, p=0.5)) == "NegativeBinomial(r=3, p=0.5)" + + +def test_property_setters_update_and_validate(): + dist = NegativeBinomial(r=3, p=0.5) + dist.r = 5 + dist.p = 0.7 + assert dist.r == 5 + assert dist.p == 0.7 + with pytest.raises(ValueError, match="r must be positive"): + dist.r = 0 + with pytest.raises(ValueError, match=r"p must be in \[0, 1\]"): + dist.p = 1.5 # --------------------------------------------------------------------------- -# Validation behavior +# PMF # --------------------------------------------------------------------------- -@pytest.mark.parametrize("r,p", [ - (1, 0.0), - (5, 0.5), - (10, 1.0), - (5, 0.3), -]) -def test_validate_params_accepts_valid_values(r, p): - NegativeBinomial._validate_params(r=r, p=p) +@pytest.mark.parametrize("r, p", PARAMS) +@pytest.mark.parametrize("k", [0, 1, 2, 3, 5, 10]) +def test_pmf_matches_closed_form(k, r, p): + assert NegativeBinomial(r, p).pmf_scalar(k) == pytest.approx( + nbinom_pmf(k, r, p), **EXACT + ) -@pytest.mark.parametrize("r,p", [ - (0, 0.5), - (-1, 0.5), - (5, -0.1), - (5, 1.1), -]) -def test_validate_params_rejects_invalid_values(r, p): - with pytest.raises(ValueError): - NegativeBinomial._validate_params(r=r, p=p) +# The summation stops at k = 160 rather than running to convergence because the +# PMF overflows to inf/nan beyond that point; see the overflow section below. +# Every parameter set here has negligible mass past k = 160 (< 1e-20). +SUMMATION_LIMIT = 160 -# --------------------------------------------------------------------------- -# Instance method delegation -# --------------------------------------------------------------------------- +@pytest.mark.parametrize("r, p", PARAMS) +def test_pmf_sums_to_one(r, p): + dist = NegativeBinomial(r, p) + total = sum(dist.pmf_scalar(k) for k in range(0, SUMMATION_LIMIT)) + assert total == pytest.approx(1.0, **ITERATIVE) -def test_pmf_scalar_delegates_to_core(mock_core): - mock_core.negative_binomial_pmf_scalar.return_value = 0.186 - nb = NegativeBinomial(r=5, p=0.6) - result = nb.pmf_scalar(3) +@pytest.mark.parametrize("r, p", PARAMS) +@pytest.mark.parametrize("k", [0, 2, 7]) +def test_pmf_is_a_probability(k, r, p): + assert 0.0 <= NegativeBinomial(r, p).pmf_scalar(k) <= 1.0 - mock_core.negative_binomial_pmf_scalar.assert_called_once_with(3, 5, 0.6) - assert result == 0.186 +def test_pmf_with_unit_r_is_geometric_on_failures(): + """NegativeBinomial(1, p) gives P(X = k) = p (1-p)^k.""" + dist = NegativeBinomial(1, 0.3) + for k in range(10): + assert dist.pmf_scalar(k) == pytest.approx(0.3 * 0.7 ** k, **EXACT) -def test_cdf_scalar_delegates_to_core(mock_core): - mock_core.negative_binomial_cdf_scalar.return_value = 0.663 - nb = NegativeBinomial(r=5, p=0.6) - result = nb.cdf_scalar(3) +# --------------------------------------------------------------------------- +# CDF +# --------------------------------------------------------------------------- - mock_core.negative_binomial_cdf_scalar.assert_called_once_with(3, 5, 0.6) - assert result == 0.663 +@pytest.mark.parametrize("r, p", PARAMS) +@pytest.mark.parametrize("k", [0, 1, 3, 6, 12]) +def test_cdf_matches_summed_pmf(k, r, p): + assert NegativeBinomial(r, p).cdf_scalar(k) == pytest.approx( + nbinom_cdf(k, r, p), **ITERATIVE + ) -def test_mean_delegates_to_core(mock_core): - mock_core.negative_binomial_mean.return_value = 3.333 +@pytest.mark.parametrize("r, p", PARAMS) +def test_cdf_is_monotonic_and_bounded(r, p): + dist = NegativeBinomial(r, p) + values = [dist.cdf_scalar(k) for k in range(0, 25)] + assert values == sorted(values) + assert all(0.0 <= v <= 1.0 for v in values) - nb = NegativeBinomial(r=5, p=0.6) - result = nb.mean() - mock_core.negative_binomial_mean.assert_called_once_with(5, 0.6) - assert result == 3.333 +@pytest.mark.parametrize("r, p", [(3, 0.5), (5, 0.7), (4, 0.9)]) +def test_cdf_approaches_one_in_the_tail(r, p): + assert NegativeBinomial(r, p).cdf_scalar(150) == pytest.approx(1.0, abs=1e-9) -def test_variance_delegates_to_core(mock_core): - mock_core.negative_binomial_variance.return_value = 5.556 +# --------------------------------------------------------------------------- +# Numeric range +# +# KNOWN BUG: the PMF evaluates C(k + r - 1, k) using raw factorials, so the +# intermediate (k + r - 1)! overflows a double once k + r - 1 > 170. The result +# is inf at k = 170 and nan from k = 200 onward, and cdf_scalar(200) is nan. +# +# The coefficient itself is small -- for r = 3, k = 200 it is C(202, 200) = +# 20301 and the true PMF is 1.58e-57, comfortably inside double range. The +# Binomial implementation handles n = 1000 without difficulty, so this is a +# defect in the negative binomial routine rather than an inherent limit. +# --------------------------------------------------------------------------- - nb = NegativeBinomial(r=5, p=0.6) - result = nb.variance() +OVERFLOW_BUG = "C(k + r - 1, k) computed via raw factorials; overflows for k >= 170" - mock_core.negative_binomial_variance.assert_called_once_with(5, 0.6) - assert result == 5.556 +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason=OVERFLOW_BUG) +@pytest.mark.parametrize("k", [170, 200, 500, 1000]) +def test_pmf_stays_finite_in_the_far_tail(k): + value = NegativeBinomial(3, 0.5).pmf_scalar(k) + assert math.isfinite(value) + assert value >= 0.0 -def test_stddev_delegates_to_core(mock_core): - mock_core.negative_binomial_stddev.return_value = 2.357 - nb = NegativeBinomial(r=5, p=0.6) - result = nb.stddev() +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason=OVERFLOW_BUG) +@pytest.mark.parametrize("k", [200, 500]) +def test_pmf_matches_closed_form_in_the_far_tail(k): + assert NegativeBinomial(3, 0.5).pmf_scalar(k) == pytest.approx( + nbinom_pmf(k, 3, 0.5), **ITERATIVE + ) - mock_core.negative_binomial_stddev.assert_called_once_with(5, 0.6) - assert result == 2.357 +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason=OVERFLOW_BUG) +@pytest.mark.parametrize("k", [200, 500]) +def test_cdf_stays_finite_in_the_far_tail(k): + assert NegativeBinomial(3, 0.5).cdf_scalar(k) == pytest.approx(1.0, abs=1e-9) -def test_mgf_scalar_delegates_to_core(mock_core): - mock_core.negative_binomial_mgf_scalar.return_value = 1.845 - nb = NegativeBinomial(r=5, p=0.6) - result = nb.mgf_scalar(0.1) +# --------------------------------------------------------------------------- +# Moments +# --------------------------------------------------------------------------- - mock_core.negative_binomial_mgf_scalar.assert_called_once_with(0.1, 5, 0.6) - assert result == 1.845 +@pytest.mark.parametrize("r, p", PARAMS) +def test_moments_match_closed_forms(r, p): + dist = NegativeBinomial(r, p) + assert dist.mean() == pytest.approx(r * (1 - p) / p, **EXACT) + assert dist.variance() == pytest.approx(r * (1 - p) / p ** 2, **EXACT) + assert dist.stddev() == pytest.approx( + math.sqrt(r * (1 - p) / p ** 2), **EXACT + ) -def test_cgf_scalar_delegates_to_core(mock_core): - mock_core.negative_binomial_cgf_scalar.return_value = 0.612 +@pytest.mark.parametrize("r, p", [(3, 0.5), (2, 0.3), (5, 0.7)]) +def test_mean_equals_pmf_weighted_sum(r, p): + dist = NegativeBinomial(r, p) + expectation = sum(k * dist.pmf_scalar(k) for k in range(0, SUMMATION_LIMIT)) + assert dist.mean() == pytest.approx(expectation, rel=1e-8) - nb = NegativeBinomial(r=5, p=0.6) - result = nb.cgf_scalar(0.1) - mock_core.negative_binomial_cgf_scalar.assert_called_once_with(0.1, 5, 0.6) - assert result == 0.612 +# --------------------------------------------------------------------------- +# MGF / CGF +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("r, p", [(3, 0.5), (2, 0.3), (5, 0.7)]) +@pytest.mark.parametrize("t", [-1.0, -0.2, 0.0, 0.1]) +def test_mgf_matches_closed_form(t, r, p): + assert NegativeBinomial(r, p).mgf_scalar(t) == pytest.approx( + nbinom_mgf(t, r, p), **EXACT + ) + + +@pytest.mark.parametrize("r, p", [(3, 0.5), (2, 0.3), (5, 0.7)]) +@pytest.mark.parametrize("t", [-1.0, -0.2, 0.0, 0.1]) +def test_cgf_matches_closed_form(t, r, p): + assert NegativeBinomial(r, p).cgf_scalar(t) == pytest.approx( + math.log(nbinom_mgf(t, r, p)), **EXACT + ) -def test_sample_delegates_to_core(mock_core): - mock_core.negative_binomial_sample.return_value = 4 +@pytest.mark.parametrize("r, p", PARAMS) +def test_mgf_at_zero_is_one(r, p): + dist = NegativeBinomial(r, p) + assert dist.mgf_scalar(0.0) == pytest.approx(1.0, **EXACT) + assert dist.cgf_scalar(0.0) == pytest.approx(0.0, abs=1e-12) - nb = NegativeBinomial(r=5, p=0.6) - result = nb.sample() - mock_core.negative_binomial_sample.assert_called_once_with(5, 0.6) - assert result == 4 +@pytest.mark.parametrize("r, p", [(3, 0.5), (5, 0.7)]) +def test_cgf_is_log_of_mgf(r, p): + dist = NegativeBinomial(r, p) + for t in (-0.5, -0.1, 0.05): + assert dist.cgf_scalar(t) == pytest.approx( + math.log(dist.mgf_scalar(t)), **EXACT + ) # --------------------------------------------------------------------------- -# Validation enforcement in classmethods +# Classmethods # --------------------------------------------------------------------------- -@pytest.mark.parametrize( - "method, args", - [ - (NegativeBinomial._pmf_scalar, (3, 0, 0.5)), - (NegativeBinomial._pmf_scalar, (3, -1, 0.5)), - (NegativeBinomial._pmf_scalar, (3, 5, -0.1)), - (NegativeBinomial._pmf_scalar, (3, 5, 1.1)), - (NegativeBinomial._cdf_scalar, (3, 0, 0.5)), - (NegativeBinomial._cdf_scalar, (3, 5, -0.5)), - (NegativeBinomial._cdf_scalar, (3, 5, 1.5)), - (NegativeBinomial._mean, (0, 0.5)), - (NegativeBinomial._mean, (-1, 0.5)), - (NegativeBinomial._mean, (5, -0.1)), - (NegativeBinomial._mean, (5, 1.1)), - (NegativeBinomial._variance, (0, 0.5)), - (NegativeBinomial._variance, (5, -0.5)), - (NegativeBinomial._stddev, (-1, 0.5)), - (NegativeBinomial._stddev, (5, 1.5)), - (NegativeBinomial._mgf_scalar, (0.1, 0, 0.5)), - (NegativeBinomial._mgf_scalar, (0.1, 5, -0.1)), - (NegativeBinomial._cgf_scalar, (0.1, -1, 0.5)), - (NegativeBinomial._cgf_scalar, (0.1, 5, 1.1)), - (NegativeBinomial._sample, (0, 0.5)), - (NegativeBinomial._sample, (5, -0.5)), - ], -) -def test_classmethods_reject_invalid_parameters(method, args): - with pytest.raises(ValueError): - method(*args) +@pytest.mark.parametrize("r, p", PARAMS) +def test_classmethods_agree_with_instances(r, p): + dist = NegativeBinomial(r, p) + assert NegativeBinomial._pmf_scalar(3, r, p) == pytest.approx( + dist.pmf_scalar(3), **EXACT + ) + assert NegativeBinomial._cdf_scalar(3, r, p) == pytest.approx( + dist.cdf_scalar(3), **ITERATIVE + ) + assert NegativeBinomial._mean(r, p) == pytest.approx(dist.mean(), **EXACT) + assert NegativeBinomial._variance(r, p) == pytest.approx( + dist.variance(), **EXACT + ) + assert NegativeBinomial._stddev(r, p) == pytest.approx(dist.stddev(), **EXACT) + assert NegativeBinomial._mgf_scalar(0.1, r, p) == pytest.approx( + dist.mgf_scalar(0.1), **EXACT + ) + assert NegativeBinomial._cgf_scalar(0.1, r, p) == pytest.approx( + dist.cgf_scalar(0.1), **EXACT + ) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pmf_scalar", (1, 0, 0.5)), + ("_cdf_scalar", (1, -2, 0.5)), + ("_mean", (0, 0.5)), + ("_variance", (-1, 0.5)), + ("_stddev", (0, 0.5)), + ("_mgf_scalar", (0.1, 0, 0.5)), + ("_cgf_scalar", (0.1, -3, 0.5)), + ("_sample", (0, 0.5)), +]) +def test_classmethods_reject_non_positive_r(method_name, args): + with pytest.raises(ValueError, match="r must be positive"): + getattr(NegativeBinomial, method_name)(*args) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pmf_scalar", (1, 3, -0.1)), + ("_cdf_scalar", (1, 3, 1.5)), + ("_mean", (3, 2.0)), + ("_variance", (3, -1.0)), + ("_sample", (3, 1.7)), +]) +def test_classmethods_reject_out_of_range_p(method_name, args): + with pytest.raises(ValueError, match=r"p must be in \[0, 1\]"): + getattr(NegativeBinomial, method_name)(*args) # --------------------------------------------------------------------------- -# Classmethod delegation with valid parameters +# Sampling # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, args, value", [ - ("_pmf_scalar", "negative_binomial_pmf_scalar", (3, 5, 0.6), 0.186), - ("_cdf_scalar", "negative_binomial_cdf_scalar", (3, 5, 0.6), 0.663), - ("_mean", "negative_binomial_mean", (5, 0.6), 3.333), - ("_variance", "negative_binomial_variance", (5, 0.6), 5.556), - ("_stddev", "negative_binomial_stddev", (5, 0.6), 2.357), - ("_mgf_scalar", "negative_binomial_mgf_scalar", (0.1, 5, 0.6), 1.845), - ("_cgf_scalar", "negative_binomial_cgf_scalar", (0.1, 5, 0.6), 0.612), - ("_sample", "negative_binomial_sample", (5, 0.6), 4), -]) -def test_classmethods_delegate_to_core(mock_core, method_name, core_method_name, args, value): - getattr(mock_core, core_method_name).return_value = value - method = getattr(NegativeBinomial, method_name) - result = method(*args) - getattr(mock_core, core_method_name).assert_called_once_with(*args) - assert result == value +@pytest.mark.parametrize("r, p", [(3, 0.5), (2, 0.3), (5, 0.9)]) +def test_sample_lies_within_the_support(r, p): + dist = NegativeBinomial(r, p) + for _ in range(300): + value = dist.sample() + assert isinstance(value, int) + assert value >= 0 # --------------------------------------------------------------------------- -# Slots behavior +# Slots # --------------------------------------------------------------------------- def test_slots_prevent_dynamic_attributes(): - nb = NegativeBinomial(r=5, p=0.6) with pytest.raises(AttributeError): - nb.extra = 123 + NegativeBinomial(3, 0.5).extra = 123 diff --git a/python/tests/test_poisson.py b/python/tests/test_poisson.py index 935daa1..306c480 100644 --- a/python/tests/test_poisson.py +++ b/python/tests/test_poisson.py @@ -1,183 +1,316 @@ -# tests/python_modules/test_poisson.py -from unittest.mock import patch +"""Numeric tests for the Poisson distribution against closed-form references.""" +import math + +import numpy as np import pytest -import fastdist.distributions.poisson as poisson_module +from conftest import EXACT, ITERATIVE from fastdist.distributions.poisson import Poisson -@pytest.fixture -def mock_core(): - """Patch the internal C++ core for the duration of each test.""" - with patch.object(poisson_module, "_core", create=True) as mock: - yield mock +# --------------------------------------------------------------------------- +# Closed-form references +# --------------------------------------------------------------------------- + +def poisson_pmf(k: int, lam: float) -> float: + """P(X = k) = lambda^k e^(-lambda) / k!""" + return lam ** k * math.exp(-lam) / math.factorial(k) + + +def poisson_cdf(k: int, lam: float) -> float: + return sum(poisson_pmf(i, lam) for i in range(0, k + 1)) + + +def poisson_mgf(t: float, lam: float) -> float: + """M(t) = exp(lambda (e^t - 1))""" + return math.exp(lam * (math.exp(t) - 1)) + + +RATES = [0.5, 1.0, 4.0, 10.0] # --------------------------------------------------------------------------- -# Constructor & representation +# Constructor, properties, representation # --------------------------------------------------------------------------- def test_init_valid_parameter(): - p = Poisson(lambda_=3.5) - assert p.lambda_ == 3.5 + assert Poisson(lambda_=4.0).lambda_ == 4.0 -@pytest.mark.parametrize("lambda_", [0, -1, -0.5]) -def test_init_invalid_lambda_raises(lambda_): +@pytest.mark.parametrize("lam", [0.0, -0.1, -5.0]) +def test_init_rejects_non_positive_rate(lam): with pytest.raises(ValueError, match="lambda_ must be positive"): - Poisson(lambda_=lambda_) + Poisson(lambda_=lam) + + +@pytest.mark.parametrize("bad", ["4.0", None, object()]) +def test_init_rejects_non_real_rate(bad): + with pytest.raises(TypeError, match="lambda_ must be a real number"): + Poisson(lambda_=bad) def test_repr(): - p = Poisson(lambda_=2.0) - assert repr(p) == "Poisson(lambda_=2.0)" + assert repr(Poisson(lambda_=4.0)) == "Poisson(lambda_=4.0)" + + +def test_rate_setter_updates_and_validates(): + dist = Poisson(lambda_=4.0) + dist.lambda_ = 7.0 + assert dist.lambda_ == 7.0 + with pytest.raises(ValueError, match="lambda_ must be positive"): + dist.lambda_ = 0.0 # --------------------------------------------------------------------------- -# Validation behavior +# PMF # --------------------------------------------------------------------------- -@pytest.mark.parametrize("lambda_", [0.0001, 1, 10]) -def test_validate_params_accepts_valid_values(lambda_): - Poisson._validate_params(lambda_=lambda_) +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("k", [0, 1, 2, 5, 10, 20]) +def test_pmf_matches_closed_form(k, lam): + assert Poisson(lam).pmf(k) == pytest.approx(poisson_pmf(k, lam), **EXACT) -@pytest.mark.parametrize("lambda_", [0, -1, -0.5]) -def test_validate_params_rejects_invalid_values(lambda_): - with pytest.raises(ValueError, match="lambda_ must be positive"): - Poisson._validate_params(lambda_=lambda_) +@pytest.mark.parametrize("lam", RATES) +def test_pmf_sums_to_one(lam): + dist = Poisson(lam) + total = sum(dist.pmf(k) for k in range(0, 200)) + assert total == pytest.approx(1.0, **ITERATIVE) + + +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("k", [0, 3, 12]) +def test_pmf_is_a_probability(k, lam): + assert 0.0 <= Poisson(lam).pmf(k) <= 1.0 + + +@pytest.mark.parametrize("lam", RATES) +def test_pmf_at_zero_is_exp_minus_lambda(lam): + assert Poisson(lam).pmf(0) == pytest.approx(math.exp(-lam), **EXACT) + + +@pytest.mark.parametrize("lam", RATES) +def test_pmf_recurrence(lam): + """P(k) = P(k-1) * lambda / k, independent of the factorial formulation.""" + dist = Poisson(lam) + for k in range(1, 25): + assert dist.pmf(k) == pytest.approx(dist.pmf(k - 1) * lam / k, **ITERATIVE) # --------------------------------------------------------------------------- -# Instance method delegation +# CDF # --------------------------------------------------------------------------- -def test_pmf_scalar_delegates_to_core(mock_core): - mock_core.poisson_pmf_scalar.return_value = 0.2 +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("k", [0, 1, 3, 8, 20]) +def test_cdf_matches_summed_pmf(k, lam): + assert Poisson(lam).cdf(k) == pytest.approx(poisson_cdf(k, lam), **ITERATIVE) - result = Poisson._pmf_scalar(3, 2.0) - mock_core.poisson_pmf_scalar.assert_called_once_with(3, 2.0) - assert result == 0.2 +@pytest.mark.parametrize("lam", RATES) +def test_cdf_is_monotonic_and_bounded(lam): + dist = Poisson(lam) + values = [dist.cdf(k) for k in range(0, 40)] + assert values == sorted(values) + assert all(0.0 <= v <= 1.0 for v in values) -def test_cdf_scalar_delegates_to_core(mock_core): - mock_core.poisson_cdf_scalar.return_value = 0.8 +@pytest.mark.parametrize("lam", RATES) +def test_cdf_approaches_one_in_the_tail(lam): + assert Poisson(lam).cdf(150) == pytest.approx(1.0, abs=1e-9) - result = Poisson._cdf_scalar(3, 2.0) - mock_core.poisson_cdf_scalar.assert_called_once_with(3, 2.0) - assert result == 0.8 +# --------------------------------------------------------------------------- +# Array API +# --------------------------------------------------------------------------- +@pytest.mark.parametrize("lam", RATES) +def test_pmf_array_matches_scalar_evaluation(lam): + dist = Poisson(lam) + ks = [0, 1, 2, 5, 10] + result = dist.pmf(ks) + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + np.testing.assert_allclose(result, [poisson_pmf(k, lam) for k in ks], rtol=1e-12) -def test_mean_delegates_to_core(mock_core): - mock_core.poisson_mean.return_value = 2.0 - p = Poisson(lambda_=2.0) - result = p.mean() +@pytest.mark.parametrize("lam", RATES) +def test_cdf_array_matches_scalar_evaluation(lam): + dist = Poisson(lam) + ks = [0, 1, 3, 8] + np.testing.assert_allclose( + dist.cdf(ks), [poisson_cdf(k, lam) for k in ks], rtol=1e-10 + ) - mock_core.poisson_mean.assert_called_once_with(2.0) - assert result == 2.0 +def test_array_accepts_numpy_input(): + dist = Poisson(4.0) + ks = np.array([0, 2, 5]) + np.testing.assert_allclose( + dist.pmf(ks), [poisson_pmf(int(k), 4.0) for k in ks], rtol=1e-12 + ) -def test_variance_delegates_to_core(mock_core): - mock_core.poisson_variance.return_value = 2.0 - p = Poisson(lambda_=2.0) - result = p.variance() +def test_empty_array_returns_empty_array(): + result = Poisson(4.0).pmf([]) + assert isinstance(result, np.ndarray) + assert result.size == 0 - mock_core.poisson_variance.assert_called_once_with(2.0) - assert result == 2.0 +def test_step_size_offsets_each_element_by_its_index(): + """With step_size s, element i is evaluated at x[i] + s*i.""" + dist = Poisson(4.0) + result = dist.pmf([0, 0, 0], 1) + np.testing.assert_allclose( + result, + [poisson_pmf(0, 4.0), poisson_pmf(1, 4.0), poisson_pmf(2, 4.0)], + rtol=1e-12, + ) -def test_stddev_delegates_to_core(mock_core): - mock_core.poisson_stddev.return_value = 1.414 - p = Poisson(lambda_=2.0) - result = p.stddev() +def test_step_size_must_be_an_integer(): + with pytest.raises(TypeError, match="step_size must be an integer"): + Poisson(4.0).pmf([0, 1, 2], 0.5) - mock_core.poisson_stddev.assert_called_once_with(2.0) - assert result == 1.414 +def test_array_rejects_two_dimensional_input(): + with pytest.raises(ValueError, match="must be 1-dimensional"): + Poisson(4.0).pmf([[1, 2], [3, 4]]) -def test_mgf_scalar_delegates_to_core(mock_core): - mock_core.poisson_mgf_scalar.return_value = 7.389 - result = Poisson._mgf_scalar(0.5, 2.0) +def test_array_rejects_non_numeric_input(): + with pytest.raises(TypeError, match="must be numeric"): + Poisson(4.0).pmf(["a", "b"]) - mock_core.poisson_mgf_scalar.assert_called_once_with(0.5, 2.0) - assert result == 7.389 +def test_rejects_none_input(): + with pytest.raises(TypeError, match="cannot be None"): + Poisson(4.0).pmf(None) -def test_cgf_scalar_delegates_to_core(mock_core): - mock_core.poisson_cgf_scalar.return_value = 2.297 - result = Poisson._cgf_scalar(0.5, 2.0) +# --------------------------------------------------------------------------- +# Moments +# --------------------------------------------------------------------------- - mock_core.poisson_cgf_scalar.assert_called_once_with(0.5, 2.0) - assert result == 2.297 +@pytest.mark.parametrize("lam", RATES) +def test_moments_match_closed_forms(lam): + dist = Poisson(lam) + assert dist.mean() == pytest.approx(lam, **EXACT) + assert dist.variance() == pytest.approx(lam, **EXACT) + assert dist.stddev() == pytest.approx(math.sqrt(lam), **EXACT) -def test_sample_delegates_to_core(mock_core): - mock_core.poisson_sample.return_value = 3 +@pytest.mark.parametrize("lam", RATES) +def test_mean_equals_pmf_weighted_sum(lam): + dist = Poisson(lam) + expectation = sum(k * dist.pmf(k) for k in range(0, 200)) + assert dist.mean() == pytest.approx(expectation, rel=1e-9) - p = Poisson(lambda_=2.0) - result = p.sample() - mock_core.poisson_sample.assert_called_once_with(2.0) - assert result == 3 +@pytest.mark.parametrize("lam, override", [(4.0, 7.0), (1.0, 0.5)]) +def test_moment_parameter_override(lam, override): + dist = Poisson(lam) + assert dist.mean(override) == pytest.approx(override, **EXACT) + assert dist.variance(override) == pytest.approx(override, **EXACT) + assert dist.stddev(override) == pytest.approx(math.sqrt(override), **EXACT) + assert dist.lambda_ == lam # the override must not mutate the instance + + +def test_moment_override_validates(): + with pytest.raises(ValueError, match="lambda_ must be positive"): + Poisson(4.0).mean(-1.0) # --------------------------------------------------------------------------- -# Validation enforcement in classmethods +# MGF / CGF # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, invalid_lambda", [ - ("mean", 0), - ("mean", -1), - ("variance", -0.5), - ("sample", -1), -]) -def test_instance_methods_reject_invalid_params(method_name, invalid_lambda): - # You can't even initialize the class with these values based on your code - with pytest.raises(ValueError, match="lambda_ must be positive"): - p = Poisson(lambda_=invalid_lambda) - method = getattr(p, method_name) - method() +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("t", [-1.0, -0.25, 0.0, 0.1, 0.5]) +def test_mgf_matches_closed_form(t, lam): + assert Poisson(lam).mgf(t) == pytest.approx(poisson_mgf(t, lam), **EXACT) + + +@pytest.mark.parametrize("lam", RATES) +@pytest.mark.parametrize("t", [-1.0, -0.25, 0.0, 0.1, 0.5]) +def test_cgf_matches_closed_form(t, lam): + assert Poisson(lam).cgf(t) == pytest.approx( + lam * (math.exp(t) - 1), **EXACT + ) + + +@pytest.mark.parametrize("lam", RATES) +def test_mgf_at_zero_is_one(lam): + dist = Poisson(lam) + assert dist.mgf(0.0) == pytest.approx(1.0, **EXACT) + assert dist.cgf(0.0) == pytest.approx(0.0, abs=1e-12) + + +@pytest.mark.parametrize("lam", [1.0, 4.0]) +def test_cgf_is_log_of_mgf(lam): + dist = Poisson(lam) + for t in (-0.5, 0.1, 0.4): + assert dist.cgf(t) == pytest.approx(math.log(dist.mgf(t)), **EXACT) # --------------------------------------------------------------------------- -# Classmethod delegation with valid parameters +# Classmethods # --------------------------------------------------------------------------- -@pytest.mark.parametrize("method_name, core_method_name, lambda_val, value", [ - ("mean", "poisson_mean", 2.0, 2.0), - ("variance", "poisson_variance", 2.0, 2.0), - ("stddev", "poisson_stddev", 2.0, 1.414), - ("sample", "poisson_sample", 2.0, 3), +@pytest.mark.parametrize("lam", RATES) +def test_classmethods_agree_with_instances(lam): + dist = Poisson(lam) + assert Poisson._pmf_scalar(3, lam) == pytest.approx(dist.pmf(3), **EXACT) + assert Poisson._cdf_scalar(3, lam) == pytest.approx(dist.cdf(3), **ITERATIVE) + assert Poisson._mgf_scalar(0.1, lam) == pytest.approx(dist.mgf(0.1), **EXACT) + assert Poisson._cgf_scalar(0.1, lam) == pytest.approx(dist.cgf(0.1), **EXACT) + + +@pytest.mark.parametrize("lam", [1.0, 4.0]) +def test_cpu_batch_classmethods_match_scalars(lam): + ks = [0, 1, 2, 5] + np.testing.assert_allclose( + Poisson._pmf_cpu(ks, lam), [poisson_pmf(k, lam) for k in ks], rtol=1e-12 + ) + np.testing.assert_allclose( + Poisson._cdf_cpu(ks, lam), [poisson_cdf(k, lam) for k in ks], rtol=1e-10 + ) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pmf_scalar", (1, 0.0)), + ("_pmf_scalar", (1, -2.0)), + ("_cdf_scalar", (1, 0.0)), + ("_mgf_scalar", (0.1, -1.0)), + ("_cgf_scalar", (0.1, 0.0)), ]) -def test_instance_methods_delegate_to_core(mock_core, method_name, core_method_name, lambda_val, value): - # 1. Setup the mock return - getattr(mock_core, core_method_name).return_value = value +def test_classmethods_reject_non_positive_rate(method_name, args): + with pytest.raises(ValueError, match="lambda_ must be positive"): + getattr(Poisson, method_name)(*args) + + +def test_is_cuda_available_returns_bool(): + assert isinstance(Poisson.is_cuda_available(), bool) - # 2. Instantiate the class (this provides the 'self' that was missing) - p = Poisson(lambda_=lambda_val) - # 3. Call the method on the instance - method = getattr(p, method_name) - result = method() # No args needed, it uses p.lambda_ +# --------------------------------------------------------------------------- +# Sampling +# --------------------------------------------------------------------------- - # 4. Assert - getattr(mock_core, core_method_name).assert_called_once_with(lambda_val) - assert result == value +@pytest.mark.parametrize("lam", RATES) +def test_sample_is_a_non_negative_integer(lam): + dist = Poisson(lam) + for _ in range(300): + value = dist.sample() + assert isinstance(value, int) + assert value >= 0 # --------------------------------------------------------------------------- -# Slots behavior +# Slots # --------------------------------------------------------------------------- def test_slots_prevent_dynamic_attributes(): - p = Poisson(lambda_=2.0) with pytest.raises(AttributeError): - p.extra = 123 + Poisson(4.0).extra = 123 diff --git a/python/tests/test_uniform.py b/python/tests/test_uniform.py index 0ecd719..10a0e3e 100644 --- a/python/tests/test_uniform.py +++ b/python/tests/test_uniform.py @@ -1,154 +1,365 @@ -from unittest.mock import patch +""" +Numeric tests for the continuous Uniform distribution against closed-form +references. +The support is the interval [a, b], so f(x) = 1/(b-a) inside it and 0 outside. +""" + +import math + +import numpy as np import pytest -# Changed: Imported from uniform module -import fastdist.distributions.uniform as uniform_module +from conftest import EXACT, ITERATIVE from fastdist.distributions.uniform import Uniform -@pytest.fixture -def mock_core(): - """ - Patch the internal C++ core for the duration of each test. - """ - # Changed: Patching the uniform module's core - with patch.object(uniform_module, "_core", autospec=True) as mock: - yield mock +# --------------------------------------------------------------------------- +# Closed-form references +# --------------------------------------------------------------------------- + +def uniform_pdf(x: float, a: float, b: float) -> float: + return 1.0 / (b - a) if a <= x <= b else 0.0 + + +def uniform_cdf(x: float, a: float, b: float) -> float: + if x < a: + return 0.0 + if x > b: + return 1.0 + return (x - a) / (b - a) + + +def uniform_mgf(t: float, a: float, b: float) -> float: + """M(t) = (e^(tb) - e^(ta)) / (t(b-a)); the t=0 singularity is removable.""" + if t == 0.0: + return 1.0 + return (math.exp(t * b) - math.exp(t * a)) / (t * (b - a)) + + +PARAMS = [(0.0, 1.0), (1.0, 3.0), (-2.0, 2.0), (-5.0, -1.0), (0.5, 10.0)] # --------------------------------------------------------------------------- -# Constructor & representation +# Constructor, properties, representation # --------------------------------------------------------------------------- def test_init_valid_parameters(): - # Changed: Uses 'a' and 'b' parameters - u = Uniform(a=0.0, b=1.0) - assert u.a == 0.0 - assert u.b == 1.0 - - -# Changed: Tests for invalid 'a' and 'b' relation (a >= b) -@pytest.mark.parametrize( - "a, b", - [ - (1.0, 1.0), # a == b - (2.0, 1.0), # a > b - (0.0, 0.0), - (-1.0, -1.0), - ], -) -def test_init_invalid_params_raises(a, b): + dist = Uniform(a=1.0, b=3.0) + assert dist.a == 1.0 + assert dist.b == 3.0 + + +@pytest.mark.parametrize("a, b", [(1.0, 1.0), (3.0, 1.0), (0.0, -1.0)]) +def test_init_rejects_non_increasing_bounds(a, b): with pytest.raises(ValueError, match="a must be less than b"): Uniform(a=a, b=b) +@pytest.mark.parametrize("bad", ["1.0", object()]) +def test_init_rejects_non_real_bounds(bad): + with pytest.raises(TypeError, match="must be a real number"): + Uniform(a=bad, b=3.0) + + +def test_init_rejects_none_bound(): + """ + None is rejected, but only incidentally: _validate_params skips its checks + when a bound is None, so the failure surfaces from float() rather than as + the intended "a must be a real number". + """ + with pytest.raises(TypeError): + Uniform(a=None, b=3.0) + + def test_repr(): - # Changed: Uses 'a' and 'b' in the repr string - u = Uniform(a=1.5, b=2.5) - assert repr(u) == "Uniform(a=1.5, b=2.5)" + assert repr(Uniform(a=1.0, b=3.0)) == "Uniform(a=1.0, b=3.0)" + + +def test_property_setters_update_values(): + dist = Uniform(a=1.0, b=3.0) + dist.a = 0.0 + dist.b = 5.0 + assert dist.a == 0.0 + assert dist.b == 5.0 + + +# KNOWN BUG: each setter validates only the bound being assigned, never the +# a < b relationship against the other one. Uniform(1.0, 3.0) can be driven to +# a = 10.0, b = -10.0 -- a state the constructor rejects outright. Once there, +# pdf, cdf, mean, variance and sample all return nan rather than raising, so the +# corruption propagates silently. +CROSS_BOUND_BUG = "setters do not re-validate a < b against the other bound" + + +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason=CROSS_BOUND_BUG) +def test_a_setter_rejects_value_above_b(): + dist = Uniform(a=1.0, b=3.0) + with pytest.raises(ValueError, match="a must be less than b"): + dist.a = 10.0 + + +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason=CROSS_BOUND_BUG) +def test_b_setter_rejects_value_below_a(): + dist = Uniform(a=1.0, b=3.0) + with pytest.raises(ValueError, match="a must be less than b"): + dist.b = -10.0 + + +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason=CROSS_BOUND_BUG + "; results become nan") +def test_instance_stays_usable_after_setter_assignments(): + dist = Uniform(a=1.0, b=3.0) + try: + dist.a = 10.0 + except ValueError: + return # rejecting the assignment is the correct behaviour + assert math.isfinite(dist.mean()) # --------------------------------------------------------------------------- -# Validation behavior +# Scalar PDF / CDF # --------------------------------------------------------------------------- -# Changed: Tests for a >= b -@pytest.mark.parametrize( - "a, b", - [ - (1.0, 1.0), - (2.0, 1.0), - ], -) -def test_validate_params_rejects_a_greater_equal_b(a, b): - with pytest.raises(ValueError): - Uniform._validate_params(a=a, b=b) +@pytest.mark.parametrize("a, b", PARAMS) +def test_pdf_is_constant_inside_the_support(a, b): + dist = Uniform(a, b) + height = 1.0 / (b - a) + for frac in (0.0, 0.25, 0.5, 0.75, 1.0): + x = a + frac * (b - a) + assert dist.pdf(x) == pytest.approx(height, **EXACT) + +@pytest.mark.parametrize("a, b", PARAMS) +def test_pdf_is_zero_outside_the_support(a, b): + dist = Uniform(a, b) + assert dist.pdf(a - 1.0) == pytest.approx(0.0, abs=1e-15) + assert dist.pdf(b + 1.0) == pytest.approx(0.0, abs=1e-15) -# Changed: Tests for a < b -def test_validate_params_accepts_a_less_than_b(): - Uniform._validate_params(a=1.0, b=2.0) # should not raise + +@pytest.mark.parametrize("a, b", PARAMS) +def test_pdf_integrates_to_one(a, b): + """The density is constant, so the integral is just height * width.""" + dist = Uniform(a, b) + assert dist.pdf((a + b) / 2) * (b - a) == pytest.approx(1.0, **EXACT) + + +@pytest.mark.parametrize("a, b", PARAMS) +def test_cdf_matches_closed_form(a, b): + dist = Uniform(a, b) + for frac in (-0.5, 0.0, 0.25, 0.5, 0.75, 1.0, 1.5): + x = a + frac * (b - a) + assert dist.cdf(x) == pytest.approx(uniform_cdf(x, a, b), **EXACT) + + +@pytest.mark.parametrize("a, b", PARAMS) +def test_cdf_saturates_at_the_bounds(a, b): + dist = Uniform(a, b) + assert dist.cdf(a) == pytest.approx(0.0, abs=1e-15) + assert dist.cdf(b) == pytest.approx(1.0, **EXACT) + assert dist.cdf(a - 10.0) == pytest.approx(0.0, abs=1e-15) + assert dist.cdf(b + 10.0) == pytest.approx(1.0, **EXACT) + + +@pytest.mark.parametrize("a, b", PARAMS) +def test_cdf_is_monotonic_and_bounded(a, b): + dist = Uniform(a, b) + xs = [a + frac * (b - a) for frac in (-0.5, 0.0, 0.3, 0.6, 1.0, 1.5)] + values = [dist.cdf(x) for x in xs] + assert values == sorted(values) + assert all(0.0 <= v <= 1.0 for v in values) + + +@pytest.mark.parametrize("a, b", PARAMS) +def test_cdf_at_the_midpoint_is_one_half(a, b): + assert Uniform(a, b).cdf((a + b) / 2) == pytest.approx(0.5, **EXACT) # --------------------------------------------------------------------------- -# Instance method delegation +# Array API # --------------------------------------------------------------------------- -def test_pdf_scalar_delegates_to_core(mock_core): - # Changed: Mock the uniform specific core function - mock_core.uniform_pdf_scalar.return_value = 0.5 # Example for U(0, 2) at x=1 +@pytest.mark.parametrize("a, b", PARAMS) +def test_pdf_array_matches_scalar_evaluation(a, b): + dist = Uniform(a, b) + xs = [a - 1.0, a, (a + b) / 2, b, b + 1.0] + result = dist.pdf(xs) + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + np.testing.assert_allclose( + result, [uniform_pdf(x, a, b) for x in xs], rtol=1e-12, atol=1e-15 + ) + + +@pytest.mark.parametrize("a, b", PARAMS) +def test_cdf_array_matches_scalar_evaluation(a, b): + dist = Uniform(a, b) + xs = [a - 1.0, a, (a + b) / 2, b, b + 1.0] + np.testing.assert_allclose( + dist.cdf(xs), [uniform_cdf(x, a, b) for x in xs], rtol=1e-12, atol=1e-15 + ) + + +def test_array_accepts_numpy_input(): + dist = Uniform(1.0, 3.0) + xs = np.array([0.0, 2.0, 5.0]) + np.testing.assert_allclose(dist.cdf(xs), [0.0, 0.5, 1.0], rtol=1e-12, atol=1e-15) + + +def test_empty_array_returns_empty_array(): + result = Uniform(1.0, 3.0).pdf([]) + assert isinstance(result, np.ndarray) + assert result.size == 0 + + +def test_step_size_offsets_each_element_by_its_index(): + """With step_size s, element i is evaluated at x[i] + s*i.""" + dist = Uniform(0.0, 4.0) + result = dist.cdf([0.0, 0.0, 0.0], 1.0) + np.testing.assert_allclose(result, [0.0, 0.25, 0.5], rtol=1e-12, atol=1e-15) - u = Uniform(a=0.0, b=2.0) - result = u.pdf(1.0) - # Changed: Check for (x, a, b) call - mock_core.uniform_pdf_scalar.assert_called_once_with(1.0, 0.0, 2.0) - assert result == 0.5 +def test_array_rejects_two_dimensional_input(): + with pytest.raises(ValueError, match="must be 1-dimensional"): + Uniform(1.0, 3.0).pdf([[1.0, 2.0], [3.0, 4.0]]) -# Note: Uniform does not have logpdf_scalar +def test_array_rejects_non_numeric_input(): + with pytest.raises(TypeError, match="must be numeric"): + Uniform(1.0, 3.0).pdf(["a", "b"]) -def test_cdf_scalar_delegates_to_core(mock_core): - # Changed: Mock the uniform specific core function - mock_core.uniform_cdf_scalar.return_value = 0.5 - u = Uniform(a=0.0, b=2.0) - result = u.cdf(1.0) +def test_rejects_none_input(): + with pytest.raises(TypeError, match="must not be None"): + Uniform(1.0, 3.0).pdf(None) - # Changed: Check for (x, a, b) call - mock_core.uniform_cdf_scalar.assert_called_once_with(1.0, 0.0, 2.0) - assert result == 0.5 +# --------------------------------------------------------------------------- +# Moments +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("a, b", PARAMS) +def test_moments_match_closed_forms(a, b): + dist = Uniform(a, b) + assert dist.mean() == pytest.approx((a + b) / 2.0, **EXACT) + assert dist.variance() == pytest.approx((b - a) ** 2 / 12.0, **EXACT) + assert dist.stddev() == pytest.approx((b - a) / math.sqrt(12.0), **EXACT) + + +@pytest.mark.parametrize("a, b", [(0.0, 10.0), (-1.0, 1.0)]) +def test_moment_parameter_override(a, b): + dist = Uniform(1.0, 3.0) + assert dist.mean(a, b) == pytest.approx((a + b) / 2.0, **EXACT) + assert dist.variance(a, b) == pytest.approx((b - a) ** 2 / 12.0, **EXACT) + assert dist.stddev(a, b) == pytest.approx((b - a) / math.sqrt(12.0), **EXACT) + assert (dist.a, dist.b) == (1.0, 3.0) # override must not mutate the instance + + +def test_moment_override_validates(): + with pytest.raises(ValueError, match="a must be less than b"): + Uniform(1.0, 3.0).mean(5.0, 2.0) -# Note: Uniform does not have z_score # --------------------------------------------------------------------------- -# Statistical properties +# MGF / CGF # --------------------------------------------------------------------------- -def test_mean_delegates_to_core(mock_core): - # Changed: Mock the uniform specific core function - mock_core.uniform_mean.return_value = 1.5 +@pytest.mark.parametrize("a, b", PARAMS) +@pytest.mark.parametrize("t", [-1.0, -0.25, 0.25, 1.0]) +def test_mgf_matches_closed_form(t, a, b): + assert Uniform(a, b).mgf(t) == pytest.approx(uniform_mgf(t, a, b), **ITERATIVE) - u = Uniform(a=1.0, b=2.0) - result = u.mean() - # Changed: Check for (a, b) call - mock_core.uniform_mean.assert_called_once_with(1.0, 2.0) - assert result == 1.5 +@pytest.mark.parametrize("a, b", PARAMS) +def test_mgf_at_zero_is_one(a, b): + """t = 0 is a removable 0/0 singularity in the closed form.""" + dist = Uniform(a, b) + assert dist.mgf(0.0) == pytest.approx(1.0, **EXACT) + assert dist.cgf(0.0) == pytest.approx(0.0, abs=1e-12) -def test_variance_delegates_to_core(mock_core): - # Changed: Mock the uniform specific core function - mock_core.uniform_variance.return_value = 0.083333333 +@pytest.mark.parametrize("a, b", PARAMS) +@pytest.mark.parametrize("t", [-1.0, -0.25, 0.25, 1.0]) +def test_cgf_is_log_of_mgf(t, a, b): + dist = Uniform(a, b) + assert dist.cgf(t) == pytest.approx(math.log(dist.mgf(t)), **ITERATIVE) - u = Uniform(a=1.0, b=2.0) - result = u.variance() - # Changed: Check for (a, b) call - mock_core.uniform_variance.assert_called_once_with(1.0, 2.0) - assert result == 0.083333333 +def test_mgf_array_matches_scalar_evaluation(): + dist = Uniform(1.0, 3.0) + ts = [-1.0, -0.25, 0.0, 0.25, 1.0] + np.testing.assert_allclose( + dist.mgf(ts), [uniform_mgf(t, 1.0, 3.0) for t in ts], rtol=1e-10 + ) -def test_stddev_delegates_to_core(mock_core): - # Changed: Mock the uniform specific core function - mock_core.uniform_stddev.return_value = 0.288675135 +# --------------------------------------------------------------------------- +# Classmethods +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("a, b", PARAMS) +def test_classmethods_agree_with_instances(a, b): + dist = Uniform(a, b) + mid = (a + b) / 2 + assert Uniform._pdf_scalar(mid, a, b) == pytest.approx(dist.pdf(mid), **EXACT) + assert Uniform._cdf_scalar(mid, a, b) == pytest.approx(dist.cdf(mid), **EXACT) + assert Uniform._mgf_scalar(0.25, a, b) == pytest.approx( + dist.mgf(0.25), **ITERATIVE + ) + assert Uniform._cgf_scalar(0.25, a, b) == pytest.approx( + dist.cgf(0.25), **ITERATIVE + ) + + +@pytest.mark.parametrize("a, b", [(0.0, 1.0), (1.0, 3.0)]) +def test_cpu_batch_classmethods_match_scalars(a, b): + xs = [a - 0.5, a, (a + b) / 2, b, b + 0.5] + np.testing.assert_allclose( + Uniform._pdf_cpu(xs, a, b), + [uniform_pdf(x, a, b) for x in xs], + rtol=1e-12, atol=1e-15, + ) + np.testing.assert_allclose( + Uniform._cdf_cpu(xs, a, b), + [uniform_cdf(x, a, b) for x in xs], + rtol=1e-12, atol=1e-15, + ) + + +@pytest.mark.parametrize("method_name, args", [ + ("_pdf_scalar", (1.0, 3.0, 1.0)), + ("_cdf_scalar", (1.0, 5.0, 2.0)), + ("_mgf_scalar", (0.25, 3.0, 1.0)), + ("_cgf_scalar", (0.25, 3.0, 1.0)), +]) +def test_classmethods_reject_non_increasing_bounds(method_name, args): + with pytest.raises(ValueError, match="a must be less than b"): + getattr(Uniform, method_name)(*args) + - u = Uniform(a=1.0, b=2.0) - result = u.stddev() +def test_is_cuda_available_returns_bool(): + assert isinstance(Uniform.is_cuda_available(), bool) + + +# --------------------------------------------------------------------------- +# Sampling +# --------------------------------------------------------------------------- - # Changed: Check for (a, b) call - mock_core.uniform_stddev.assert_called_once_with(1.0, 2.0) - assert result == 0.288675135 +@pytest.mark.parametrize("a, b", PARAMS) +def test_sample_lies_within_the_support(a, b): + dist = Uniform(a, b) + for _ in range(300): + value = dist.sample() + assert math.isfinite(value) + assert a <= value <= b # --------------------------------------------------------------------------- -# Slots behavior +# Slots # --------------------------------------------------------------------------- def test_slots_prevent_dynamic_attributes(): - # Changed: Initialize Uniform - u = Uniform(a=0.0, b=1.0) with pytest.raises(AttributeError): - u.foo = 123 + Uniform(1.0, 3.0).extra = 123 diff --git a/python/tests/test_utils.py b/python/tests/test_utils.py index 39ba233..7f0bcff 100644 --- a/python/tests/test_utils.py +++ b/python/tests/test_utils.py @@ -1,148 +1,395 @@ -from unittest.mock import patch +"""Numeric tests for the Utils statistical helpers against closed-form references.""" +import math + +import numpy as np import pytest -# Imported from utils module -import fastdist.distributions.utils as utils_module +from conftest import EXACT, ITERATIVE from fastdist.distributions.utils import Utils -@pytest.fixture -def mock_core(): +# --------------------------------------------------------------------------- +# Chebyshev bound +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("variance, k", [(4.0, 2.0), (1.0, 1.0), (9.0, 3.0), (0.25, 0.5)]) +def test_chebyshev_bound_matches_closed_form(variance, k): + """P(|X - mu| >= k) <= sigma^2 / k^2""" + assert Utils.chebyshev_bound(variance, k) == pytest.approx( + variance / k ** 2, **EXACT + ) + + +def test_chebyshev_bound_decreases_with_k(): + values = [Utils.chebyshev_bound(4.0, k) for k in (1.0, 2.0, 4.0, 8.0)] + assert values == sorted(values, reverse=True) + + +# --------------------------------------------------------------------------- +# Bayes' rule +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("p_b_given_a, p_a, p_b", [ + (0.8, 0.25, 0.4), + (0.9, 0.1, 0.2), + (0.5, 0.5, 0.5), +]) +def test_bayes_rule_matches_closed_form(p_b_given_a, p_a, p_b): + """P(A|B) = P(B|A) P(A) / P(B)""" + assert Utils.bayes_rule(p_b_given_a, p_a, p_b) == pytest.approx( + p_b_given_a * p_a / p_b, **EXACT + ) + + +def test_bayes_rule_is_identity_when_evidence_matches_prior(): + """If P(B|A) == P(B), then A and B are independent and P(A|B) == P(A).""" + assert Utils.bayes_rule(0.4, 0.25, 0.4) == pytest.approx(0.25, **EXACT) + + +# --------------------------------------------------------------------------- +# Law of total probability +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("p_a, p_b_given_a", [ + ([0.3, 0.7], [0.2, 0.5]), + ([0.5, 0.5], [1.0, 0.0]), + ([0.2, 0.3, 0.5], [0.1, 0.4, 0.9]), +]) +def test_law_of_total_probability_matches_closed_form(p_a, p_b_given_a): + """P(B) = sum_i P(A_i) P(B|A_i)""" + expected = sum(a * b for a, b in zip(p_a, p_b_given_a)) + assert Utils.law_of_total_probability(p_a, p_b_given_a) == pytest.approx( + expected, **EXACT + ) + + +def test_law_of_total_probability_over_a_partition_is_bounded(): + result = Utils.law_of_total_probability([0.25, 0.25, 0.5], [0.9, 0.1, 0.4]) + assert 0.0 <= result <= 1.0 + + +# KNOWN BUG: the signature annotates both arguments as +# Union[Real, Sequence[Real]], but the implementation validates them as +# sequences and forwards them to a vector-only binding, so scalar inputs raise +# TypeError from pybind11. Either the annotation or the implementation is wrong. +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="scalar inputs are annotated but unsupported") +def test_law_of_total_probability_accepts_scalars(): + assert Utils.law_of_total_probability(0.3, 0.2) == pytest.approx(0.06, **EXACT) + + +# --------------------------------------------------------------------------- +# Sigmoid / logit +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("x", [-5.0, -1.0, 0.0, 1.0, 2.0, 5.0]) +def test_sigmoid_matches_closed_form(x): + assert Utils.sigmoid(x) == pytest.approx(1.0 / (1.0 + math.exp(-x)), **EXACT) + + +def test_sigmoid_at_zero_is_one_half(): + assert Utils.sigmoid(0.0) == pytest.approx(0.5, **EXACT) + + +@pytest.mark.parametrize("x", [-8.0, -2.0, 0.0, 2.0, 8.0]) +def test_sigmoid_is_symmetric(x): + """sigmoid(-x) == 1 - sigmoid(x)""" + assert Utils.sigmoid(-x) == pytest.approx(1.0 - Utils.sigmoid(x), **EXACT) + + +@pytest.mark.parametrize("x", [-10.0, -1.0, 0.0, 1.0, 10.0]) +def test_sigmoid_is_bounded(x): + assert 0.0 < Utils.sigmoid(x) < 1.0 + + +@pytest.mark.parametrize("p", [0.01, 0.25, 0.5, 0.75, 0.99]) +def test_logit_matches_closed_form(p): + assert Utils.logit(p) == pytest.approx(math.log(p / (1 - p)), **EXACT) + + +def test_logit_at_one_half_is_zero(): + assert Utils.logit(0.5) == pytest.approx(0.0, abs=1e-15) + + +@pytest.mark.parametrize("p", [0.05, 0.3, 0.5, 0.8, 0.95]) +def test_logit_inverts_sigmoid(p): + """sigmoid(logit(p)) == p""" + assert Utils.sigmoid(Utils.logit(p)) == pytest.approx(p, **ITERATIVE) + + +@pytest.mark.parametrize("x", [-3.0, -0.5, 0.0, 0.5, 3.0]) +def test_sigmoid_inverts_logit(x): + """logit(sigmoid(x)) == x""" + assert Utils.logit(Utils.sigmoid(x)) == pytest.approx(x, **ITERATIVE) + + +@pytest.mark.parametrize("p", [0.0, 1.0, -0.5, 1.5]) +def test_logit_returns_nan_outside_the_open_unit_interval(p): """ - Patch the internal C++ core for the duration of each test. + The backend signals a domain error by returning NaN rather than raising. + This pins the current contract; if it is changed to raise, update this test. """ - with patch.object(utils_module, "_core", autospec=True) as mock: - yield mock + assert math.isnan(Utils.logit(p)) + + +# KNOWN BUG: sigmoid is annotated Union[Real, Sequence[Real]] but its body calls +# float() on the validated input, so any sequence raises TypeError. The array +# path exists separately as sigmoid_cpu. +@pytest.mark.known_bug +@pytest.mark.xfail(strict=True, reason="sequence inputs are annotated but unsupported") +def test_sigmoid_accepts_a_sequence(): + result = Utils.sigmoid([0.0, 2.0]) + np.testing.assert_allclose(result, [0.5, 1.0 / (1.0 + math.exp(-2.0))], rtol=1e-12) + + +# --------------------------------------------------------------------------- +# Batch sigmoid / logit +# --------------------------------------------------------------------------- + +def test_sigmoid_cpu_matches_scalar_evaluation(): + xs = [-5.0, -1.0, 0.0, 1.0, 5.0] + result = Utils.sigmoid_cpu(xs) + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + np.testing.assert_allclose(result, [Utils.sigmoid(x) for x in xs], rtol=1e-12) + + +def test_logit_cpu_matches_scalar_evaluation(): + ps = [0.1, 0.25, 0.5, 0.75, 0.9] + result = Utils.logit_cpu(ps) + assert isinstance(result, np.ndarray) + np.testing.assert_allclose(result, [Utils.logit(p) for p in ps], rtol=1e-12) + + +def test_sigmoid_cpu_accepts_numpy_input(): + xs = np.array([-1.0, 0.0, 1.0]) + np.testing.assert_allclose( + Utils.sigmoid_cpu(xs), [Utils.sigmoid(float(x)) for x in xs], rtol=1e-12 + ) + + +def test_is_cuda_available_returns_bool(): + assert isinstance(Utils.is_cuda_available(), bool) # --------------------------------------------------------------------------- -# Core Utilities Delegation Tests +# Distances and similarity # --------------------------------------------------------------------------- -def test_chebyshev_bound_delegates_to_core(mock_core): - mock_core.chebyshev_bound.return_value = 0.75 - result = Utils.chebyshev_bound(variance=1.0, k=2.0) - mock_core.chebyshev_bound.assert_called_once_with(1.0, 2.0) - assert result == 0.75 +@pytest.mark.parametrize("x, y, expected", [ + ([0, 0], [3, 4], 5.0), + ([1, 2, 3], [1, 2, 3], 0.0), + ([0, 0, 0], [1, 1, 1], math.sqrt(3)), + ([-1, -1], [2, 3], 5.0), +]) +def test_euclidean_distance_matches_closed_form(x, y, expected): + assert Utils.euclidean_distance(x, y) == pytest.approx(expected, **EXACT) + + +@pytest.mark.parametrize("x, y, expected", [ + ([0, 0], [3, 4], 7.0), + ([1, 2, 3], [1, 2, 3], 0.0), + ([0, 0, 0], [1, 1, 1], 3.0), + ([-1, -1], [2, 3], 7.0), +]) +def test_manhattan_distance_matches_closed_form(x, y, expected): + assert Utils.manhattan_distance(x, y) == pytest.approx(expected, **EXACT) + + +@pytest.mark.parametrize("x, y", [ + ([1, 2, 3], [4, 5, 6]), + ([0, 1], [1, 0]), + ([-3, 7, 2], [5, -1, 4]), +]) +def test_euclidean_distance_matches_numpy(x, y): + assert Utils.euclidean_distance(x, y) == pytest.approx( + float(np.linalg.norm(np.array(x, float) - np.array(y, float))), **EXACT + ) + +@pytest.mark.parametrize("x, y", [([1, 2, 3], [4, 5, 6]), ([-3, 7, 2], [5, -1, 4])]) +def test_distances_are_symmetric(x, y): + assert Utils.euclidean_distance(x, y) == pytest.approx( + Utils.euclidean_distance(y, x), **EXACT + ) + assert Utils.manhattan_distance(x, y) == pytest.approx( + Utils.manhattan_distance(y, x), **EXACT + ) -def test_bayes_rule_delegates_to_core(mock_core): - mock_core.bayes_rule.return_value = 0.3 - result = Utils.bayes_rule(p_B_given_A=0.6, p_A=0.5, p_B=1.0) - mock_core.bayes_rule.assert_called_once_with(0.6, 0.5, 1.0) - assert result == 0.3 +@pytest.mark.parametrize("x, y", [([1, 2, 3], [4, 5, 6]), ([0, 1], [1, 0])]) +def test_euclidean_never_exceeds_manhattan(x, y): + """The L2 norm is bounded above by the L1 norm.""" + assert Utils.euclidean_distance(x, y) <= Utils.manhattan_distance(x, y) + 1e-12 -def test_sigmoid_delegates_to_core(mock_core): - mock_core.sigmoid.return_value = 0.5 - result = Utils.sigmoid(x=0.0) - mock_core.sigmoid.assert_called_once_with(0.0) - assert result == 0.5 +@pytest.mark.parametrize("x, y, expected", [ + ([1, 0], [0, 1], 0.0), + ([1, 2, 3], [1, 2, 3], 1.0), + ([1, 0], [-1, 0], -1.0), + ([1, 1], [2, 2], 1.0), +]) +def test_cosine_similarity_matches_closed_form(x, y, expected): + assert Utils.cosine_similarity(x, y) == pytest.approx(expected, **ITERATIVE) -def test_logit_delegates_to_core(mock_core): - mock_core.logit.return_value = 0.0 - result = Utils.logit(p=0.5) - mock_core.logit.assert_called_once_with(0.5) - assert result == 0.0 +@pytest.mark.parametrize("x, y", [([1, 2, 3], [4, 5, 6]), ([-3, 7, 2], [5, -1, 4])]) +def test_cosine_similarity_matches_numpy(x, y): + xv, yv = np.array(x, float), np.array(y, float) + expected = float(xv @ yv / (np.linalg.norm(xv) * np.linalg.norm(yv))) + assert Utils.cosine_similarity(x, y) == pytest.approx(expected, **ITERATIVE) -def test_euclidean_distance_delegates_to_core(mock_core): - mock_core.euclidean_distance.return_value = 5.0 - x = [3.0, 0.0] - y = [0.0, 4.0] +@pytest.mark.parametrize("x, y", [([1, 2, 3], [4, 5, 6]), ([0, 1], [1, 0])]) +def test_cosine_similarity_is_bounded(x, y): + assert -1.0 <= Utils.cosine_similarity(x, y) <= 1.0 - result = Utils.euclidean_distance(x, y) - called_args, _ = mock_core.euclidean_distance.call_args - import numpy as np - assert np.array_equal(called_args[0], np.array(x)) - assert np.array_equal(called_args[1], np.array(y)) +def test_cosine_similarity_is_scale_invariant(): + base = Utils.cosine_similarity([1, 2, 3], [4, 5, 6]) + scaled = Utils.cosine_similarity([10, 20, 30], [4, 5, 6]) + assert base == pytest.approx(scaled, **ITERATIVE) - assert result == 5.0 +def test_cosine_similarity_of_a_zero_vector_is_nan(): + """The zero vector has no direction; the backend signals this with NaN.""" + assert math.isnan(Utils.cosine_similarity([0, 0], [1, 1])) -def test_manhattan_distance_delegates_to_core(mock_core): - mock_core.manhattan_distance.return_value = 7.0 - x = [3.0, 0.0] - y = [0.0, 4.0] +@pytest.mark.parametrize("fn", [ + Utils.euclidean_distance, + Utils.manhattan_distance, + Utils.cosine_similarity, +]) +def test_distance_functions_reject_mismatched_lengths(fn): + with pytest.raises(ValueError, match="x and y must have the same length"): + fn([1, 2], [1, 2, 3]) - result = Utils.manhattan_distance(x, y) - called_args, _ = mock_core.manhattan_distance.call_args - import numpy as np - assert np.array_equal(called_args[0], np.array(x)) - assert np.array_equal(called_args[1], np.array(y)) +@pytest.mark.parametrize("fn", [ + Utils.euclidean_distance, + Utils.manhattan_distance, + Utils.cosine_similarity, +]) +def test_distance_functions_accept_numpy_input(fn): + assert fn(np.array([1.0, 2.0]), np.array([3.0, 4.0])) == pytest.approx( + fn([1.0, 2.0], [3.0, 4.0]), **EXACT + ) - assert result == 7.0 +# --------------------------------------------------------------------------- +# Summary statistics +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("mean, stddev", [(10.0, 2.0), (5.0, 5.0), (100.0, 1.0)]) +def test_coefficient_of_variation_matches_closed_form(mean, stddev): + """CV = sigma / mu""" + assert Utils.coefficient_of_variation(mean, stddev) == pytest.approx( + stddev / mean, **EXACT + ) -def test_coefficient_of_variation_delegates_to_core(mock_core): - mock_core.coefficient_of_variation.return_value = 0.2 - result = Utils.coefficient_of_variation(mean=10.0, stddev=2.0) - mock_core.coefficient_of_variation.assert_called_once_with(10.0, 2.0) - assert result == 0.2 +@pytest.mark.parametrize("mean_x, mean_y, e_xy", [ + (2.0, 3.0, 7.0), + (0.0, 0.0, 1.0), + (-1.0, 4.0, 2.0), +]) +def test_covariance_matches_closed_form(mean_x, mean_y, e_xy): + """Cov(X, Y) = E[XY] - E[X]E[Y]""" + assert Utils.covariance(mean_x, mean_y, e_xy) == pytest.approx( + e_xy - mean_x * mean_y, **EXACT + ) -def test_covariance_delegates_to_core(mock_core): - mock_core.covariance.return_value = 2.5 - result = Utils.covariance(mean_x=1.0, mean_y=2.0, E_xy=4.5) - mock_core.covariance.assert_called_once_with(1.0, 2.0, 4.5) - assert result == 2.5 + +def test_covariance_is_zero_for_independent_variables(): + """If E[XY] == E[X]E[Y] the covariance vanishes.""" + assert Utils.covariance(2.0, 3.0, 6.0) == pytest.approx(0.0, abs=1e-15) # --------------------------------------------------------------------------- -# Combinatorics and Special Functions Delegation Tests +# Combinatorics # --------------------------------------------------------------------------- -def test_choose_delegates_to_core(mock_core): - mock_core.choose.return_value = 10 - result = Utils.choose(n=5, k=2) - mock_core.choose.assert_called_once_with(5, 2) - assert result == 10 +@pytest.mark.parametrize("n, k", [(10, 3), (5, 0), (5, 5), (52, 5), (100, 50)]) +def test_choose_matches_reference(n, k): + assert Utils.choose(n, k) == pytest.approx(float(math.comb(n, k)), rel=1e-12) + + +@pytest.mark.parametrize("n, k", [(5, 7), (3, 10)]) +def test_choose_is_zero_when_k_exceeds_n(n, k): + assert Utils.choose(n, k) == pytest.approx(0.0, abs=1e-15) + + +@pytest.mark.parametrize("n, k", [(10, 3), (52, 5), (20, 10)]) +def test_choose_is_symmetric(n, k): + """C(n, k) == C(n, n-k)""" + assert Utils.choose(n, k) == pytest.approx(Utils.choose(n, n - k), rel=1e-12) + + +@pytest.mark.parametrize("n, k", [(10, 3), (5, 0), (5, 5), (20, 4)]) +def test_permutation_matches_reference(n, k): + assert Utils.permutation(n, k) == pytest.approx( + float(math.perm(n, k)), rel=1e-12 + ) -def test_permutation_delegates_to_core(mock_core): - mock_core.permutation.return_value = 60 - result = Utils.permutation(n=5, k=3) - mock_core.permutation.assert_called_once_with(5, 3) - assert result == 60 +@pytest.mark.parametrize("n, k", [(10, 3), (20, 4), (52, 5)]) +def test_permutation_equals_choose_times_factorial(n, k): + """P(n, k) == C(n, k) * k!""" + assert Utils.permutation(n, k) == pytest.approx( + Utils.choose(n, k) * math.factorial(k), rel=1e-10 + ) -def test_factorial_delegates_to_core(mock_core): - mock_core.factorial.return_value = 120 - result = Utils.factorial(n=5) - mock_core.factorial.assert_called_once_with(5) - assert result == 120 +@pytest.mark.parametrize("n", [0, 1, 5, 10, 20, 100, 170]) +def test_factorial_matches_reference(n): + assert Utils.factorial(n) == pytest.approx(float(math.factorial(n)), rel=1e-12) -def test_gamma_delegates_to_core(mock_core): - # Gamma(4) = 3! = 6 - mock_core.gamma.return_value = 6.0 - result = Utils.gamma(x=4.0) - mock_core.gamma.assert_called_once_with(4.0) - assert result == 6.0 +def test_factorial_overflows_beyond_the_double_range(): + """171! exceeds the maximum finite double, so inf is the correct result.""" + assert math.isinf(Utils.factorial(171)) -def test_log_gamma_delegates_to_core(mock_core): - # LogGamma(4) = log(6) ≈ 1.79176 - mock_core.log_gamma.return_value = 1.79176 - result = Utils.log_gamma(x=4.0) - mock_core.log_gamma.assert_called_once_with(4.0) - assert result == 1.79176 +@pytest.mark.parametrize("n, a, b", [(3, 1.0, 2.0), (5, 2.0, 3.0), (0, 4.0, 7.0)]) +def test_binomial_theorem_matches_closed_form(n, a, b): + """The binomial theorem expands to (a + b)^n.""" + assert Utils.binomial(n, a, b) == pytest.approx((a + b) ** n, rel=1e-10) # --------------------------------------------------------------------------- -# Slots behavior (Inherited from Normal/Uniform test structure) +# Special functions # --------------------------------------------------------------------------- -def test_utils_class_has_no_slots(): - # Since Utils uses only @classmethod, it doesn't need __slots__. - # A positive test ensures it does NOT have __slots__ defined (or if it did, it would pass). - # Since no __slots__ are defined, dynamic attributes should be allowed if instantiated, - # but the class itself is generally not meant for instantiation. - # The structure suggests testing the base case: no instance attributes to slot. - assert not hasattr(Utils, '__slots__') +@pytest.mark.parametrize("x", [0.5, 1.0, 2.0, 5.0, 7.5, 10.0]) +def test_gamma_matches_reference(x): + assert Utils.gamma(x) == pytest.approx(math.gamma(x), rel=1e-12) + + +@pytest.mark.parametrize("n", [1, 2, 3, 5, 8]) +def test_gamma_of_an_integer_is_a_factorial(n): + """Gamma(n) == (n-1)!""" + assert Utils.gamma(float(n)) == pytest.approx( + float(math.factorial(n - 1)), rel=1e-10 + ) + + +def test_gamma_of_one_half_is_root_pi(): + assert Utils.gamma(0.5) == pytest.approx(math.sqrt(math.pi), rel=1e-12) + + +@pytest.mark.parametrize("x", [0.5, 1.0, 2.5, 10.0, 100.0]) +def test_log_gamma_matches_reference(x): + assert Utils.log_gamma(x) == pytest.approx(math.lgamma(x), rel=1e-12) + + +@pytest.mark.parametrize("x", [0.5, 1.5, 4.0, 9.0]) +def test_log_gamma_is_log_of_gamma(x): + assert Utils.log_gamma(x) == pytest.approx(math.log(Utils.gamma(x)), **ITERATIVE) + + +@pytest.mark.parametrize("x", [1.5, 3.0, 7.25]) +def test_gamma_recurrence(x): + """Gamma(x + 1) == x * Gamma(x)""" + assert Utils.gamma(x + 1.0) == pytest.approx(x * Utils.gamma(x), rel=1e-10) From 924316d674238ced0ebc2f689240f5fc04909a68 Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 11:58:14 -0500 Subject: [PATCH 09/13] Fix CI build job broken by the pybind11 find_package migration The build job installed pybind11 via apt (pybind11-dev) but never pip installed it, while CMakeLists.txt probes for it with `python -m pybind11 --cmakedir` under COMMAND_ERROR_IS_FATAL ANY. That probe needs the pip package, so configure aborted. The link job was unaffected because it installs requirements.txt first. CMakeLists.txt: make the probe advisory rather than fatal. Capture the exit code, set pybind11_DIR only on success, and otherwise fall through to find_package so a system pybind11 (e.g. Debian pybind11-dev) still resolves. python-distro.yml: replace the apt pybind11-dev install with `pip install pybind11`, matching pyproject.toml build-system.requires and what the link job already does. Pin the configure step to the interpreter pip installed into with -DPython_EXECUTABLE="$(which python)" so CMake cannot bind to a different Python than the one being built for. Verified both paths locally with CMake 4.4.3: the pip path configures and generates cleanly, and an interpreter without pybind11 now reaches find_package instead of aborting at the probe. Co-Authored-By: Claude Opus 5 --- .github/workflows/python-distro.yml | 9 +++++++-- CMakeLists.txt | 12 ++++++++++-- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/.github/workflows/python-distro.yml b/.github/workflows/python-distro.yml index 0800b15..4238083 100644 --- a/.github/workflows/python-distro.yml +++ b/.github/workflows/python-distro.yml @@ -24,11 +24,16 @@ jobs: - name: Install system dependencies run: | sudo apt-get update - sudo apt-get install -y cmake build-essential pybind11-dev + sudo apt-get install -y cmake build-essential + + - name: Install Python build dependencies + run: | + python -m pip install --upgrade pip + pip install pybind11 - name: Configure CMake run: | - cmake -S . -B build + cmake -S . -B build -DPython_EXECUTABLE="$(which python)" - name: Build run: | diff --git a/CMakeLists.txt b/CMakeLists.txt index 8e839be..e65fa41 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,10 +42,18 @@ set(PYBIND11_FINDPYTHON ON) if (NOT pybind11_DIR) execute_process( COMMAND "${Python_EXECUTABLE}" -m pybind11 --cmakedir - OUTPUT_VARIABLE pybind11_DIR + OUTPUT_VARIABLE _pybind11_cmakedir OUTPUT_STRIP_TRAILING_WHITESPACE - COMMAND_ERROR_IS_FATAL ANY + RESULT_VARIABLE _pybind11_probe_result + ERROR_QUIET ) + if (_pybind11_probe_result EQUAL 0) + set(pybind11_DIR "${_pybind11_cmakedir}") + else () + # No pip pybind11 in this interpreter; fall back to a system install + # (e.g. Debian's pybind11-dev), which ships its own CMake config. + message(STATUS "pip pybind11 not found; falling back to a system installation") + endif () endif () find_package(pybind11 CONFIG REQUIRED) From d6d4c46ea7ecf5de8564f46588dc4fdadfe408dd Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Sun, 30 Aug 2026 11:58:34 -0500 Subject: [PATCH 10/13] Remove the unused pybind11 git submodule The submodule at python/pybind11 was never wired into the build. The original CMakeLists.txt resolved pybind11 through FetchContent, which downloaded its own copy of v2.12.1 at configure time; there was no add_subdirectory and no reference to the submodule path anywhere in CMakeLists.txt, setup.py, pyproject.toml or build_all.ps1. It has since been replaced by find_package against the pip package, so the submodule is dead weight in every configuration. Removing it reclaims 4.3 MB from the working tree and 15 MB of cached metadata under .git/modules, taking .git from roughly 16 MB to 940 KB. Also drops `submodules: recursive` from both CI checkout steps, since there is nothing left to fetch, and updates the README clone instructions and the release-notes line that described pybind11 as a submodule. Note for existing clones: run `git submodule deinit -f python/pybind11` after pulling if git leaves a stale empty directory behind. Fresh clones are unaffected. Co-Authored-By: Claude Opus 5 --- .github/workflows/python-distro.yml | 4 ---- .gitmodules | 3 --- README.md | 7 ++----- python/pybind11 | 1 - 4 files changed, 2 insertions(+), 13 deletions(-) delete mode 100644 .gitmodules delete mode 160000 python/pybind11 diff --git a/.github/workflows/python-distro.yml b/.github/workflows/python-distro.yml index 4238083..1dda1e3 100644 --- a/.github/workflows/python-distro.yml +++ b/.github/workflows/python-distro.yml @@ -13,8 +13,6 @@ jobs: steps: - name: Checkout repository uses: actions/checkout@v4 - with: - submodules: recursive - name: Set up Python 3.14 uses: actions/setup-python@v4 @@ -46,8 +44,6 @@ jobs: steps: - name: Checkout repository uses: actions/checkout@v4 - with: - submodules: recursive - name: Set up Python uses: actions/setup-python@v4 diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index 9b655d9..0000000 --- a/.gitmodules +++ /dev/null @@ -1,3 +0,0 @@ -[submodule "python/pybind11"] - path = python/pybind11 - url = https://github.com/pybind/pybind11.git diff --git a/README.md b/README.md index 980d6b4..940ee0f 100644 --- a/README.md +++ b/README.md @@ -13,11 +13,8 @@ ### Cloning the Repository -This project uses Git submodules. Clone the repository recursively: - ```bash -git clone --recurse-submodule https://github.com/ghosteau/fastdist.git -git submodule update --init --recursive +git clone https://github.com/ghosteau/fastdist.git ``` --- @@ -216,7 +213,7 @@ Testing and CI: Python Bindings: -- Pybind11 integrated as a submodule for modular C++/Python bindings +- Pybind11 resolved from the build environment for modular C++/Python bindings - Full Python support for all currently supported builds --- diff --git a/python/pybind11 b/python/pybind11 deleted file mode 160000 index 10f8708..0000000 --- a/python/pybind11 +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 10f8708979fb424163bc333eaa6f96f5f40a58cd From 5117ce892bb19fceef135cda00596eb00ae288db Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Mon, 31 Aug 2026 14:42:03 -0500 Subject: [PATCH 11/13] Changed config.py to only save the cuda limits if the user runs auto_tune(). Otherwise, it will use the default When a user runs auto_tune() it will write to CONFIG_FILE and return an error if it can't --- python/fastdist/config.py | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/python/fastdist/config.py b/python/fastdist/config.py index c358244..74e8e24 100644 --- a/python/fastdist/config.py +++ b/python/fastdist/config.py @@ -94,8 +94,6 @@ else: CONFIG_FILE = Path.home() / ".config" / "fastdist" / "config.json" -CONFIG_FILE.parent.mkdir(parents=True, exist_ok=True) - def _load_config(): """ @@ -117,21 +115,16 @@ def _load_config(): global CUDA_THRESHOLDS CUDA_THRESHOLDS = copy.deepcopy(_DEFAULT_CUDA_THRESHOLDS) - if CONFIG_FILE.exists(): - try: - data = json.loads(CONFIG_FILE.read_text()) - for key, value in data.items(): - if key in CUDA_THRESHOLDS and isinstance(value, dict): - CUDA_THRESHOLDS[key].update(value) - else: - CUDA_THRESHOLDS[key] = value - return - except json.JSONDecodeError: - pass - - # Corrupt or missing file, create default file - CONFIG_FILE.write_text(json.dumps(CUDA_THRESHOLDS, indent=4)) + try: + data = json.loads(CONFIG_FILE.read_text()) + except (OSError, json.JSONDecodeError): + return # absent, unreadable, or corrupt --> use defaults + for key, value in data.items(): + if key in CUDA_THRESHOLDS and isinstance(value, dict): + CUDA_THRESHOLDS[key].update(value) + else: + CUDA_THRESHOLDS[key] = value def _save_config(): """ @@ -144,6 +137,7 @@ def _save_config(): - Writes to the configuration file. """ + CONFIG_FILE.parent.mkdir(parents=True, exist_ok=True) CONFIG_FILE.write_text(json.dumps(CUDA_THRESHOLDS, indent=4)) From 1a9f8e76b3dc5503cf941ac6245817675ec0145a Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Tue, 1 Sep 2026 12:51:46 -0500 Subject: [PATCH 12/13] Fixed clang-format erroring out and added .cu/.cuh files to the formatter --- .pre-commit-config.yaml | 4 ++-- python/bindings/bindings.cpp | 2 +- src/cuda/bernoulli/cgf.cu | 15 +++++---------- src/cuda/bernoulli/mgf.cu | 15 +++++---------- src/cuda/bernoulli/pmf.cu | 15 +++++---------- src/cuda/exponential/cdf.cu | 16 ++++++---------- src/cuda/exponential/cgf.cu | 16 ++++++---------- src/cuda/exponential/mgf.cu | 16 ++++++---------- src/cuda/exponential/pdf.cu | 20 ++++++++------------ src/cuda/normal/logpdf.cu | 2 +- src/cuda/normal/mgf.cu | 2 +- src/cuda/poisson/cdf.cu | 13 ++++--------- src/cuda/poisson/cgf.cu | 13 ++++--------- src/cuda/poisson/mgf.cu | 13 ++++--------- src/cuda/poisson/pmf.cu | 13 ++++--------- src/cuda/uniform/cdf.cu | 17 ++++++----------- src/cuda/uniform/cgf.cu | 17 ++++++----------- src/cuda/uniform/mgf.cu | 17 ++++++----------- src/cuda/uniform/pdf.cu | 17 ++++++----------- src/cuda/utils/logit.cu | 3 ++- src/cuda/utils/sigmoid.cu | 3 ++- 21 files changed, 90 insertions(+), 159 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 48778f4..08483cc 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -3,5 +3,5 @@ repos: rev: v18.1.3 hooks: - id: clang-format - name: clang-format (C/C++) - files: \.(c|cc|cpp|cxx|h|hpp)$ + name: clang-format (C/C++/CUDA) + files: \.(c|cc|cpp|cxx|cu|cuh|h|hpp)$ diff --git a/python/bindings/bindings.cpp b/python/bindings/bindings.cpp index 5a2890e..136ff88 100644 --- a/python/bindings/bindings.cpp +++ b/python/bindings/bindings.cpp @@ -1,7 +1,7 @@ // CPP file to link all other bindings +#include #include #include -#include namespace py = pybind11; diff --git a/src/cuda/bernoulli/cgf.cu b/src/cuda/bernoulli/cgf.cu index fde234d..aac3176 100644 --- a/src/cuda/bernoulli/cgf.cu +++ b/src/cuda/bernoulli/cgf.cu @@ -5,13 +5,14 @@ #include #include #include -#include "fastdist/cuda/executor.cuh" #include "fastdist/cuda/bernoulli.cuh" +#include "fastdist/cuda/executor.cuh" #include "fastdist/math/constants.h" namespace fastdist::cuda::bernoulli { // CUDA kernel - __global__ void bernoulli_cgf_kernel(const double* t, double* output, const int n, const double p, const int stepSize, const int offset) { + __global__ void bernoulli_cgf_kernel(const double* t, double* output, const int n, const double p, + const int stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -31,13 +32,7 @@ namespace fastdist::cuda::bernoulli { void bernoulli_cgf_dispatcher(const double* t, double* output, const int n, const double p, const int stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - bernoulli_cgf_kernel, - t, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::COMPLEX_MATH, - p, - stepSize); + execute_cuda_kernel(bernoulli_cgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::COMPLEX_MATH, p, stepSize); } } // namespace fastdist::cuda::bernoulli diff --git a/src/cuda/bernoulli/mgf.cu b/src/cuda/bernoulli/mgf.cu index 06f32ca..d737db6 100644 --- a/src/cuda/bernoulli/mgf.cu +++ b/src/cuda/bernoulli/mgf.cu @@ -5,13 +5,14 @@ #include #include #include -#include "fastdist/cuda/executor.cuh" #include "fastdist/cuda/bernoulli.cuh" +#include "fastdist/cuda/executor.cuh" #include "fastdist/math/constants.h" namespace fastdist::cuda::bernoulli { // CUDA kernel - __global__ void bernoulli_mgf_kernel(const double* t, double* output, const int n, const double p, const int stepSize, const int offset) { + __global__ void bernoulli_mgf_kernel(const double* t, double* output, const int n, const double p, + const int stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -31,13 +32,7 @@ namespace fastdist::cuda::bernoulli { void bernoulli_mgf_dispatcher(const double* t, double* output, const int n, const double p, const int stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - bernoulli_mgf_kernel, - t, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::SIMPLE_MATH, - p, - stepSize); + execute_cuda_kernel(bernoulli_mgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::SIMPLE_MATH, p, stepSize); } } // namespace fastdist::cuda::bernoulli diff --git a/src/cuda/bernoulli/pmf.cu b/src/cuda/bernoulli/pmf.cu index 61f1457..48dcd64 100644 --- a/src/cuda/bernoulli/pmf.cu +++ b/src/cuda/bernoulli/pmf.cu @@ -5,13 +5,14 @@ #include #include #include -#include "fastdist/cuda/executor.cuh" #include "cuda/bernoulli.cuh" +#include "fastdist/cuda/executor.cuh" #include "fastdist/math/constants.h" namespace fastdist::cuda::bernoulli { // CUDA kernel - __global__ void bernoulli_pmf_kernel(const int* k, double* output, const int n, const double p, const int stepSize, const int offset) { + __global__ void bernoulli_pmf_kernel(const int* k, double* output, const int n, const double p, const int stepSize, + const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -37,13 +38,7 @@ namespace fastdist::cuda::bernoulli { void bernoulli_pmf_dispatcher(const int* k, double* output, const int n, const double p, const int stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - bernoulli_pmf_kernel, - k, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::SIMPLE_MATH, - p, - stepSize); + execute_cuda_kernel(bernoulli_pmf_kernel, k, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::SIMPLE_MATH, p, stepSize); } } // namespace fastdist::cuda::bernoulli diff --git a/src/cuda/exponential/cdf.cu b/src/cuda/exponential/cdf.cu index 86dc688..b453078 100644 --- a/src/cuda/exponential/cdf.cu +++ b/src/cuda/exponential/cdf.cu @@ -11,7 +11,8 @@ namespace fastdist::cuda::exponential { // CUDA kernel - __global__ void exponential_cdf_kernel(const double* x, double* output, const int n, const double lambda, const double stepSize, const int offset) { + __global__ void exponential_cdf_kernel(const double* x, double* output, const int n, const double lambda, + const double stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -33,16 +34,11 @@ namespace fastdist::cuda::exponential { } // Dispatcher - void exponential_cdf_dispatcher(const double* x, double* output, const int n, const double lambda, const double stepSize) { + void exponential_cdf_dispatcher(const double* x, double* output, const int n, const double lambda, + const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - exponential_cdf_kernel, - x, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::COMPLEX_MATH, - lambda, - stepSize); + execute_cuda_kernel(exponential_cdf_kernel, x, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::COMPLEX_MATH, lambda, stepSize); } } // namespace fastdist::cuda::exponential diff --git a/src/cuda/exponential/cgf.cu b/src/cuda/exponential/cgf.cu index f137cda..e79e8fb 100644 --- a/src/cuda/exponential/cgf.cu +++ b/src/cuda/exponential/cgf.cu @@ -11,7 +11,8 @@ namespace fastdist::cuda::exponential { // CUDA kernel - __global__ void exponential_cgf_kernel(const double* t, double* output, const int n, const double lambda, const double stepSize, const int offset) { + __global__ void exponential_cgf_kernel(const double* t, double* output, const int n, const double lambda, + const double stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -28,16 +29,11 @@ namespace fastdist::cuda::exponential { } // Dispatcher - void exponential_cgf_dispatcher(const double* t, double* output, const int n, const double lambda, const double stepSize) { + void exponential_cgf_dispatcher(const double* t, double* output, const int n, const double lambda, + const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - exponential_cgf_kernel, - t, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::COMPLEX_MATH, - lambda, - stepSize); + execute_cuda_kernel(exponential_cgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::COMPLEX_MATH, lambda, stepSize); } } // namespace fastdist::cuda::exponential diff --git a/src/cuda/exponential/mgf.cu b/src/cuda/exponential/mgf.cu index 8afff85..2c2b437 100644 --- a/src/cuda/exponential/mgf.cu +++ b/src/cuda/exponential/mgf.cu @@ -11,7 +11,8 @@ namespace fastdist::cuda::exponential { // CUDA kernel - __global__ void exponential_mgf_kernel(const double* t, double* output, const int n, const double lambda, const double stepSize, const int offset) { + __global__ void exponential_mgf_kernel(const double* t, double* output, const int n, const double lambda, + const double stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -28,16 +29,11 @@ namespace fastdist::cuda::exponential { } // Dispatcher - void exponential_mgf_dispatcher(const double* t, double* output, const int n, const double lambda, const double stepSize) { + void exponential_mgf_dispatcher(const double* t, double* output, const int n, const double lambda, + const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - exponential_mgf_kernel, - t, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::SIMPLE_MATH, - lambda, - stepSize); + execute_cuda_kernel(exponential_mgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::SIMPLE_MATH, lambda, stepSize); } } // namespace fastdist::cuda::exponential diff --git a/src/cuda/exponential/pdf.cu b/src/cuda/exponential/pdf.cu index ab311ae..5e1c304 100644 --- a/src/cuda/exponential/pdf.cu +++ b/src/cuda/exponential/pdf.cu @@ -5,13 +5,14 @@ #include #include #include -#include "fastdist/cuda/executor.cuh" #include "cuda/exponential.cuh" +#include "fastdist/cuda/executor.cuh" #include "fastdist/math/constants.h" namespace fastdist::cuda::exponential { - // CUDA kernel - __global__ void exponential_pdf_kernel(const double* x, double* output, const int n, const double lambda, const double stepSize, const int offset) { + // CUDA kernel + __global__ void exponential_pdf_kernel(const double* x, double* output, const int n, const double lambda, + const double stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -33,16 +34,11 @@ namespace fastdist::cuda::exponential { } // Dispatcher - void exponential_pdf_dispatcher(const double* x, double* output, const int n, const double lambda, const double stepSize) { + void exponential_pdf_dispatcher(const double* x, double* output, const int n, const double lambda, + const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - exponential_pdf_kernel, - x, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::COMPLEX_MATH, - lambda, - stepSize); + execute_cuda_kernel(exponential_pdf_kernel, x, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::COMPLEX_MATH, lambda, stepSize); } } // namespace fastdist::cuda::exponential diff --git a/src/cuda/normal/logpdf.cu b/src/cuda/normal/logpdf.cu index 90f8f2b..63583f2 100644 --- a/src/cuda/normal/logpdf.cu +++ b/src/cuda/normal/logpdf.cu @@ -35,7 +35,7 @@ namespace fastdist::cuda::normal { void normal_logpdf_dispatcher(const double* x, double* output, const int n, const double mu, const double sigma, const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel(normal_logpdf_kernel, x, output, ctx.dev_in, ctx.dev_out, n, + execute_cuda_kernel(normal_logpdf_kernel, x, output, ctx.dev_in, ctx.dev_out, n, StreamingThresholds::COMPLEX_MATH, mu, sigma, stepSize); } } // namespace fastdist::cuda::normal diff --git a/src/cuda/normal/mgf.cu b/src/cuda/normal/mgf.cu index 09a5359..5f9be6d 100644 --- a/src/cuda/normal/mgf.cu +++ b/src/cuda/normal/mgf.cu @@ -33,7 +33,7 @@ namespace fastdist::cuda::normal { void normal_mgf_dispatcher(const double* t, double* output, const int n, const double mu, const double sigma, const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel(normal_mgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, + execute_cuda_kernel(normal_mgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, StreamingThresholds::COMPLEX_MATH, mu, sigma, stepSize); } } // namespace fastdist::cuda::normal diff --git a/src/cuda/poisson/cdf.cu b/src/cuda/poisson/cdf.cu index 27551c5..85a8976 100644 --- a/src/cuda/poisson/cdf.cu +++ b/src/cuda/poisson/cdf.cu @@ -22,7 +22,8 @@ namespace fastdist::cuda::poisson { } // CUDA kernel - __global__ void poisson_cdf_kernel(const double* x, double* output, const int n, const double lambda, const int stepSize, const int offset) { + __global__ void poisson_cdf_kernel(const double* x, double* output, const int n, const double lambda, + const int stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -54,13 +55,7 @@ namespace fastdist::cuda::poisson { void poisson_cdf_dispatcher(const double* x, double* output, const int n, const double lambda, const int stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - poisson_cdf_kernel, - x, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::SIMPLE_MATH, - lambda, - stepSize); + execute_cuda_kernel(poisson_cdf_kernel, x, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::SIMPLE_MATH, lambda, stepSize); } } // namespace fastdist::cuda::poisson diff --git a/src/cuda/poisson/cgf.cu b/src/cuda/poisson/cgf.cu index 99b0002..e6cea32 100644 --- a/src/cuda/poisson/cgf.cu +++ b/src/cuda/poisson/cgf.cu @@ -11,7 +11,8 @@ namespace fastdist::cuda::poisson { // CUDA kernel - __global__ void poisson_cgf_kernel(const double* t, double* output, const int n, const double lambda, const int stepSize, const int offset) { + __global__ void poisson_cgf_kernel(const double* t, double* output, const int n, const double lambda, + const int stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -31,13 +32,7 @@ namespace fastdist::cuda::poisson { void poisson_cgf_dispatcher(const double* t, double* output, const int n, const double lambda, const int stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - poisson_cgf_kernel, - t, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::SIMPLE_MATH, - lambda, - stepSize); + execute_cuda_kernel(poisson_cgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::SIMPLE_MATH, lambda, stepSize); } } // namespace fastdist::cuda::poisson diff --git a/src/cuda/poisson/mgf.cu b/src/cuda/poisson/mgf.cu index 1efbec3..a5b3aa1 100644 --- a/src/cuda/poisson/mgf.cu +++ b/src/cuda/poisson/mgf.cu @@ -11,7 +11,8 @@ namespace fastdist::cuda::poisson { // CUDA kernel - __global__ void poisson_mgf_kernel(const double* t, double* output, const int n, const double lambda, const int stepSize, const int offset) { + __global__ void poisson_mgf_kernel(const double* t, double* output, const int n, const double lambda, + const int stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -31,13 +32,7 @@ namespace fastdist::cuda::poisson { void poisson_mgf_dispatcher(const double* t, double* output, const int n, const double lambda, const int stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - poisson_mgf_kernel, - t, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::COMPLEX_MATH, - lambda, - stepSize); + execute_cuda_kernel(poisson_mgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::COMPLEX_MATH, lambda, stepSize); } } // namespace fastdist::cuda::poisson diff --git a/src/cuda/poisson/pmf.cu b/src/cuda/poisson/pmf.cu index 5d78415..b0877ce 100644 --- a/src/cuda/poisson/pmf.cu +++ b/src/cuda/poisson/pmf.cu @@ -11,7 +11,8 @@ namespace fastdist::cuda::poisson { // CUDA kernel - __global__ void poisson_pmf_kernel(const double* x, double* output, const int n, const double lambda, const int stepSize, const int offset) { + __global__ void poisson_pmf_kernel(const double* x, double* output, const int n, const double lambda, + const int stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -37,13 +38,7 @@ namespace fastdist::cuda::poisson { void poisson_pmf_dispatcher(const double* x, double* output, const int n, const double lambda, const int stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - poisson_pmf_kernel, - x, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::COMPLEX_MATH, - lambda, - stepSize); + execute_cuda_kernel(poisson_pmf_kernel, x, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::COMPLEX_MATH, lambda, stepSize); } } // namespace fastdist::cuda::poisson diff --git a/src/cuda/uniform/cdf.cu b/src/cuda/uniform/cdf.cu index 0bd6936..8df84a4 100644 --- a/src/cuda/uniform/cdf.cu +++ b/src/cuda/uniform/cdf.cu @@ -11,7 +11,8 @@ namespace fastdist::cuda::uniform { // CUDA kernel - __global__ void uniform_cdf_kernel(const double* x, double* output, const int n, const double a, const double b, const double stepSize, const int offset) { + __global__ void uniform_cdf_kernel(const double* x, double* output, const int n, const double a, const double b, + const double stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -36,17 +37,11 @@ namespace fastdist::cuda::uniform { } // Dispatcher - void uniform_cdf_dispatcher(const double* x, double* output, const int n, const double a, const double b, const double stepSize) { + void uniform_cdf_dispatcher(const double* x, double* output, const int n, const double a, const double b, + const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - uniform_cdf_kernel, - x, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::SIMPLE_MATH, - a, - b, - stepSize); + execute_cuda_kernel(uniform_cdf_kernel, x, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::SIMPLE_MATH, a, b, stepSize); } } // namespace fastdist::cuda::uniform diff --git a/src/cuda/uniform/cgf.cu b/src/cuda/uniform/cgf.cu index d22e6cb..d5f8a4c 100644 --- a/src/cuda/uniform/cgf.cu +++ b/src/cuda/uniform/cgf.cu @@ -24,7 +24,8 @@ namespace fastdist::cuda::uniform { } // CUDA kernel - __global__ void uniform_cgf_kernel(const double* t, double* output, const int n, const double a, const double b, const double stepSize, const int offset) { + __global__ void uniform_cgf_kernel(const double* t, double* output, const int n, const double a, const double b, + const double stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -47,17 +48,11 @@ namespace fastdist::cuda::uniform { } // Dispatcher - void uniform_cgf_dispatcher(const double* t, double* output, const int n, const double a, const double b, const double stepSize) { + void uniform_cgf_dispatcher(const double* t, double* output, const int n, const double a, const double b, + const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - uniform_cgf_kernel, - t, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::SIMPLE_MATH, - a, - b, - stepSize); + execute_cuda_kernel(uniform_cgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::SIMPLE_MATH, a, b, stepSize); } } // namespace fastdist::cuda::uniform diff --git a/src/cuda/uniform/mgf.cu b/src/cuda/uniform/mgf.cu index cb6df94..b69d100 100644 --- a/src/cuda/uniform/mgf.cu +++ b/src/cuda/uniform/mgf.cu @@ -11,7 +11,8 @@ namespace fastdist::cuda::uniform { // CUDA kernel - __global__ void uniform_mgf_kernel(const double* t, double* output, const int n, const double a, const double b, const double stepSize, const int offset) { + __global__ void uniform_mgf_kernel(const double* t, double* output, const int n, const double a, const double b, + const double stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -33,17 +34,11 @@ namespace fastdist::cuda::uniform { } // Dispatcher - void uniform_mgf_dispatcher(const double* t, double* output, const int n, const double a, const double b, const double stepSize) { + void uniform_mgf_dispatcher(const double* t, double* output, const int n, const double a, const double b, + const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - uniform_mgf_kernel, - t, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::COMPLEX_MATH, - a, - b, - stepSize); + execute_cuda_kernel(uniform_mgf_kernel, t, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::COMPLEX_MATH, a, b, stepSize); } } // namespace fastdist::cuda::uniform diff --git a/src/cuda/uniform/pdf.cu b/src/cuda/uniform/pdf.cu index a2777bf..80979c5 100644 --- a/src/cuda/uniform/pdf.cu +++ b/src/cuda/uniform/pdf.cu @@ -11,7 +11,8 @@ namespace fastdist::cuda::uniform { // CUDA kernel - __global__ void uniform_pdf_kernel(const double* x, double* output, const int n, const double a, const double b, const double stepSize, const int offset) { + __global__ void uniform_pdf_kernel(const double* x, double* output, const int n, const double a, const double b, + const double stepSize, const int offset) { int idx = blockIdx.x * blockDim.x + threadIdx.x; int global_idx = idx + offset; @@ -34,17 +35,11 @@ namespace fastdist::cuda::uniform { } // Dispatcher - void uniform_pdf_dispatcher(const double* x, double* output, const int n, const double a, const double b, const double stepSize) { + void uniform_pdf_dispatcher(const double* x, double* output, const int n, const double a, const double b, + const double stepSize) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel( - uniform_pdf_kernel, - x, - output, ctx.dev_in, ctx.dev_out, - n, - StreamingThresholds::SIMPLE_MATH, - a, - b, - stepSize); + execute_cuda_kernel(uniform_pdf_kernel, x, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::SIMPLE_MATH, a, b, stepSize); } } // namespace fastdist::cuda::uniform diff --git a/src/cuda/utils/logit.cu b/src/cuda/utils/logit.cu index c59aeea..9a2dca6 100644 --- a/src/cuda/utils/logit.cu +++ b/src/cuda/utils/logit.cu @@ -24,6 +24,7 @@ namespace fastdist::cuda::utils { void logit_dispatcher(const double* p, double* output, const int n) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel(logit_kernel, p, output, ctx.dev_in, ctx.dev_out, n, StreamingThresholds::COMPLEX_MATH); + execute_cuda_kernel(logit_kernel, p, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::COMPLEX_MATH); } } // namespace fastdist::cuda::utils diff --git a/src/cuda/utils/sigmoid.cu b/src/cuda/utils/sigmoid.cu index 42b8855..efcf23f 100644 --- a/src/cuda/utils/sigmoid.cu +++ b/src/cuda/utils/sigmoid.cu @@ -31,6 +31,7 @@ namespace fastdist::cuda::utils { void sigmoid_dispatcher(const double* x, double* output, const int n) { DeviceContext& ctx = get_context(n); - execute_cuda_kernel(sigmoid_kernel, x, output, ctx.dev_in, ctx.dev_out, n, StreamingThresholds::COMPLEX_MATH); + execute_cuda_kernel(sigmoid_kernel, x, output, ctx.dev_in, ctx.dev_out, n, + StreamingThresholds::COMPLEX_MATH); } } // namespace fastdist::cuda::utils From 86f8389277e481cc94fef96ba53510ba053ce9d6 Mon Sep 17 00:00:00 2001 From: Zach Pipes Date: Tue, 1 Sep 2026 13:18:28 -0500 Subject: [PATCH 13/13] Fixed incorrect variable name passed in bernoulli.py/_pmf_cuda --- python/fastdist/distributions/bernoulli.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/fastdist/distributions/bernoulli.py b/python/fastdist/distributions/bernoulli.py index 25822bd..79dfc51 100644 --- a/python/fastdist/distributions/bernoulli.py +++ b/python/fastdist/distributions/bernoulli.py @@ -787,7 +787,7 @@ def _pmf_cuda(cls, k: Sequence[int], p: Real, step_size: int = 0) -> NDArray[np. validated_input = cls._validate_inputs(_input=k, input_name="k", step_size=step_size) config.validate_gpu_capacity(validated_input.size, 8) - return _core.bernoulli_pmf_cuda(x=validated_input, p=p, step_size=step_size) + return _core.bernoulli_pmf_cuda(k=validated_input, p=p, step_size=step_size) @classmethod def _cdf_cuda(cls, k: Sequence[int], p: Real, step_size: int = 0) -> NDArray[np.float64]: