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..1dda1e3 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: @@ -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 @@ -24,11 +22,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: | @@ -41,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/.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/.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/.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/CMakeLists.txt b/CMakeLists.txt index 5ea00c3..e65fa41 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 # ----------------------------- @@ -19,26 +29,35 @@ 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_cmakedir + OUTPUT_STRIP_TRAILING_WHITESPACE + 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 () -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 @@ -194,6 +213,7 @@ target_include_directories(fastdist_core PUBLIC $ $ + $ ${CUDAToolkit_INCLUDE_DIRS} $ ) diff --git a/README.md b/README.md index 25c9cb8..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 ``` --- @@ -33,7 +30,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:** @@ -43,7 +40,7 @@ python3 python/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 +89,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: @@ -212,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/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..136ff88 100644 --- a/python/bindings/bindings.cpp +++ b/python/bindings/bindings.cpp @@ -1,4 +1,5 @@ // CPP file to link all other bindings +#include #include #include @@ -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 c62c58c..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", +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", "Normal", "NegativeBinomial", - "Poisson", "Uniform", "Utils"] + "Geometric", "NegativeBinomial", "Normal", + "Poisson", "Uniform", "Utils"] \ No newline at end of file diff --git a/python/fastdist/config.py b/python/fastdist/config.py index f600f87..74e8e24 100644 --- a/python/fastdist/config.py +++ b/python/fastdist/config.py @@ -1,12 +1,17 @@ +import atexit import copy import json import platform import time import numpy as np -import pynvml from pathlib import Path +try: + import pynvml +except ImportError: + pynvml = None + # ---------------------------------------------------------------------------------------------------------------------- # Constants and Defaults # ---------------------------------------------------------------------------------------------------------------------- @@ -89,8 +94,6 @@ else: CONFIG_FILE = Path.home() / ".config" / "fastdist" / "config.json" -CONFIG_FILE.parent.mkdir(parents=True, exist_ok=True) - def _load_config(): """ @@ -112,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(): """ @@ -139,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)) @@ -560,24 +559,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/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/python/fastdist/distributions/bernoulli.py b/python/fastdist/distributions/bernoulli.py index bd87815..79dfc51 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 @@ -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]: 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 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_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)), 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_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/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_package_api.py b/python/tests/test_package_api.py new file mode 100644 index 0000000..e84974d --- /dev/null +++ b/python/tests/test_package_api.py @@ -0,0 +1,36 @@ +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(): + # the top level exports every distribution class plus __version__ + assert sorted(fastdist.__all__) == sorted(EXPECTED + ["__version__"]) + + +def test_distributions_all_matches_top_level(): + # the distributions subpackage exports the classes only + assert sorted(distributions.__all__) == sorted(EXPECTED) + + +@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 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) 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 57e693c..d648bed 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,4 +5,4 @@ wheel pre-commit pytest numpy -pynvml \ No newline at end of file +nvidia-ml-py \ No newline at end of file diff --git a/setup.py b/setup.py index 1a28af6..02e5b1e 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}", ] @@ -147,17 +147,54 @@ 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!", + 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", + ], ) 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