diff --git a/.clang-format b/.clang-format new file mode 100644 index 00000000..06608729 --- /dev/null +++ b/.clang-format @@ -0,0 +1,9 @@ +# clang-format configuration for the CUDA/C++ extension (libs/ccc_cuda_ext). +# Microsoft base style already uses Allman braces and 4-space indentation, +# which matches the majority style of the existing sources, minimizing churn. +--- +BasedOnStyle: Microsoft +IndentWidth: 4 +ColumnLimit: 120 +PointerAlignment: Right +SortIncludes: CaseSensitive diff --git a/.clang-tidy b/.clang-tidy new file mode 100644 index 00000000..c26534d4 --- /dev/null +++ b/.clang-tidy @@ -0,0 +1,29 @@ +# clang-tidy configuration for the CUDA/C++ extension (libs/ccc_cuda_ext). +# +# Curated check list targeting the bug classes actually found in review +# (integer narrowing, unchecked returns). Run against the exported +# compile_commands.json, host-only pass: +# +# mamba run -n ccc-gpu-dev clang-tidy \ +# -p build \ +# --extra-arg=--cuda-host-only \ +# libs/ccc_cuda_ext/*.cu +# +# CI runs this report-only at first; promote to blocking once the baseline +# is clean. +Checks: > + -*, + clang-analyzer-*, + bugprone-*, + performance-*, + modernize-*, + readability-braces-around-statements, + readability-else-after-return, + cppcoreguidelines-init-variables, + misc-unused-parameters, + -modernize-use-trailing-return-type, + -modernize-avoid-c-arrays, + -bugprone-easily-swappable-parameters +WarningsAsErrors: '' +HeaderFilterRegex: 'libs/ccc_cuda_ext/.*\.cuh$' +FormatStyle: file diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs new file mode 100644 index 00000000..05ed6705 --- /dev/null +++ b/.git-blame-ignore-revs @@ -0,0 +1,11 @@ +# Revisions listed here are skipped by `git blame` so that large mechanical +# formatting commits do not obscure authorship of the surrounding code. +# +# Enable locally (one-time): +# git config blame.ignoreRevsFile .git-blame-ignore-revs +# +# GitHub honors this file automatically. +# +# Append the SHA of the isolated "apply ruff-format + clang-format" commit +# from the add-lint-tooling change below, e.g.: +# # style: apply ruff-format + clang-format (add-lint-tooling) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 00000000..f62e76d8 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,151 @@ +name: CI + +# Runs on every push to the default branch and on every pull request. +# +# NOTE: GPU tests (tests/gpu/**) are intentionally NOT run here — hosted GitHub +# runners have no CUDA device. They are a documented LOCAL command: +# mamba run -n ccc-gpu-dev python -m pytest tests/gpu -q -o addopts="" +# CI covers: (a) lint, (b) a no-device wheel/extension build smoke, and +# (c) the CPU-only test subset. +on: + push: + branches: [main, master] + pull_request: + +# Keep tool versions in sync with .pre-commit-config.yaml and the dev env. +env: + RUFF_VERSION: "0.15.1" + CLANG_FORMAT_VERSION: "22.1.8" + CUDA_IMAGE: "nvidia/cuda:12.5.1-devel-ubuntu22.04" + +concurrency: + group: ci-${{ github.ref }} + cancel-in-progress: true + +jobs: + # ---------------------------------------------------------------------- # + # (a) Lint: Python (ruff) + CUDA/C++ (clang-format). Pure CPU, no CUDA. + # ---------------------------------------------------------------------- # + lint: + name: Lint (ruff + clang-format) + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.11" + + - name: Install linters + run: | + python -m pip install --upgrade pip + python -m pip install "ruff==${RUFF_VERSION}" "clang-format==${CLANG_FORMAT_VERSION}" + + - name: ruff check + run: ruff check libs/ccc tests + + - name: ruff format --check + run: ruff format --check libs/ccc tests + + - name: clang-format --dry-run --Werror (CUDA/C++ extension) + run: | + clang-format --dry-run --Werror \ + libs/ccc_cuda_ext/*.cu libs/ccc_cuda_ext/*.cuh + + # ---------------------------------------------------------------------- # + # (b) Build smoke: compile the CUDA extension in a CUDA devel container. + # nvcc compiles for CMAKE_CUDA_ARCHITECTURES=75 without any GPU device + # present, so this catches compilation breakage on hosted runners. + # ---------------------------------------------------------------------- # + build: + name: Build smoke (no GPU device) + runs-on: ubuntu-latest + container: + image: nvidia/cuda:12.5.1-devel-ubuntu22.04 + steps: + - uses: actions/checkout@v4 + + - name: Install base tools + run: | + apt-get update + apt-get install -y --no-install-recommends git curl ca-certificates build-essential + + - name: Install Miniforge + run: | + curl -fsSL -o /tmp/miniforge.sh \ + https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh + bash /tmp/miniforge.sh -b -p /opt/conda + echo "/opt/conda/bin" >> "$GITHUB_PATH" + + - name: Build & install the package + env: + CUDACXX: /usr/local/cuda/bin/nvcc + CUDA_HOME: /usr/local/cuda + PATH: /opt/conda/bin:/usr/local/cuda/bin:/usr/bin:/bin + run: | + python -m pip install --upgrade pip + # Build isolation pulls the build backend (scikit-build-core, pybind11, + # cmake>=4, ninja) declared in pyproject [build-system]. + python -m pip install . -v + + - name: Import CPU path (smoke) + run: | + python - <<'PY' + import ccc + from ccc.coef import ccc as ccc_fn # triggers numba CPU compile at import + import numpy as np + print("CPU coefficient:", ccc_fn(np.random.rand(20), np.random.rand(20))) + print("import OK:", ccc.__name__) + PY + + # ---------------------------------------------------------------------- # + # (c) CPU-only tests: the CPU subset of the suite. GPU, network, and + # undeclared-dep tests are ignored (see BASELINE). The two known + # np.unicode_ test-only failures are deselected until restructure-tests + # (PR D) fixes the NumPy-2 alias. + # TODO(PR D): once pytest markers land, replace the --ignore/--deselect + # list with `-m "not gpu and not slow and not network"`. + # ---------------------------------------------------------------------- # + test-cpu: + name: CPU tests + runs-on: ubuntu-latest + container: + image: nvidia/cuda:12.5.1-devel-ubuntu22.04 + steps: + - uses: actions/checkout@v4 + + - name: Install base tools + run: | + apt-get update + apt-get install -y --no-install-recommends git curl ca-certificates build-essential + + - name: Install Miniforge + run: | + curl -fsSL -o /tmp/miniforge.sh \ + https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh + bash /tmp/miniforge.sh -b -p /opt/conda + echo "/opt/conda/bin" >> "$GITHUB_PATH" + + - name: Build & install the package (+ test deps) + env: + CUDACXX: /usr/local/cuda/bin/nvcc + CUDA_HOME: /usr/local/cuda + PATH: /opt/conda/bin:/usr/local/cuda/bin:/usr/bin:/bin + run: | + python -m pip install --upgrade pip + python -m pip install ".[test]" -v + + - name: Run CPU test subset + env: + PATH: /opt/conda/bin:/usr/local/cuda/bin:/usr/bin:/bin + run: | + python -m pytest tests/ \ + --ignore=tests/gpu \ + --ignore=tests/test_giant.py \ + --ignore=tests/test_methods.py \ + --ignore=tests/test_plots.py \ + --ignore=tests/test_corr.py \ + --deselect "tests/test_coef_pval.py::test_cm_numerical_and_categorical_features_perfect_relationship_pvalue" \ + --deselect "tests/test_coef_pval.py::test_cm_numerical_and_categorical_features_weakly_relationship_pvalue" \ + -q -ra -p no:cacheprovider -o addopts="" diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 32d67e6c..7841e739 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -10,19 +10,19 @@ repos: - id: check-toml - id: check-added-large-files - id: mixed-line-ending + # Python: ruff is the single source of truth for lint + format (black dropped). - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.8.5 # Ruff version. + rev: v0.15.1 # Ruff version (keep in sync with the dev env / CI). hooks: - id: ruff # Run the linter. types_or: [python, pyi] args: [--fix] - id: ruff-format # Run the formatter. types_or: [python, pyi] - - repo: https://github.com/psf/black - rev: 24.2.0 # Replace with desired version + # C++/CUDA: clang-format, scoped to the extension sources via file filter. + - repo: https://github.com/pre-commit/mirrors-clang-format + rev: v22.1.8 # clang-format version (matches the dev env / CI). hooks: - - id: black - types_or: [python, pyi] - - id: black-jupyter - types_or: [python, pyi, jupyter] - additional_dependencies: [".[jupyter]"] + - id: clang-format + types_or: [c++, c, cuda] + files: ^libs/ccc_cuda_ext/.*\.(cu|cuh|cpp|cc|h)$ diff --git a/CMakeLists.txt b/CMakeLists.txt index 3790c6a2..1f073ed6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,11 @@ cmake_minimum_required(VERSION 3.15...3.26) project(${SKBUILD_PROJECT_NAME} LANGUAGES CUDA CXX) +# Export compile_commands.json so clang-tidy (and editors/LSP) can analyze the +# CUDA/C++ sources. Run clang-tidy host-only, e.g.: +# clang-tidy -p build --extra-arg=--cuda-host-only libs/ccc_cuda_ext/*.cu +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + # Add this near the top of your file, after project() # Define the include directories for the whole project set(PROJECT_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/libs) diff --git a/libs/ccc/coef/__init__.py b/libs/ccc/coef/__init__.py index 29b5f439..e95faefb 100644 --- a/libs/ccc/coef/__init__.py +++ b/libs/ccc/coef/__init__.py @@ -1,7 +1,8 @@ +import numpy as np + from ccc.coef.impl import * # noqa: F403, F401 # Run CCC to initialize/compile its functions with numba from ccc.coef.impl import ccc -import numpy as np ccc(np.random.rand(10), np.random.rand(10)) diff --git a/libs/ccc/coef/impl.py b/libs/ccc/coef/impl.py index 24ed11e7..f2d92fbd 100644 --- a/libs/ccc/coef/impl.py +++ b/libs/ccc/coef/impl.py @@ -7,17 +7,16 @@ import os from collections.abc import Iterable from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor, as_completed -from typing import Union +import numpy as np from numba import njit from numba.typed import List +from numpy.typing import NDArray -import numpy as np from ccc.pytorch.core import unravel_index_2d from ccc.scipy.stats import rank from ccc.sklearn.metrics import adjusted_rand_index as ari from ccc.utils import DummyExecutor, chunker -from numpy.typing import NDArray @njit(cache=True, nogil=True) @@ -319,7 +318,7 @@ def get_coords_from_index(n_obj: int, idx: int) -> tuple[int]: def get_chunks( - iterable: Union[int, Iterable], n_threads: int, ratio: float = 1 + iterable: int | Iterable, n_threads: int, ratio: float = 1 ) -> Iterable[Iterable[int]]: """ It splits elements in an iterable in chunks according to the number of @@ -587,7 +586,7 @@ def get_n_workers(n_jobs: int | None) -> int: def ccc( x: NDArray, y: NDArray = None, - internal_n_clusters: Union[int, Iterable[int]] = None, + internal_n_clusters: int | Iterable[int] = None, return_parts: bool = False, n_chunks_threads_ratio: int = 1, n_jobs: int = 1, diff --git a/libs/ccc/coef/impl_gpu.py b/libs/ccc/coef/impl_gpu.py index 418578d8..9d720984 100644 --- a/libs/ccc/coef/impl_gpu.py +++ b/libs/ccc/coef/impl_gpu.py @@ -7,16 +7,15 @@ import os from collections.abc import Iterable from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor -from typing import Union import ccc_cuda_ext +import numpy as np from numba import njit from numba.typed import List +from numpy.typing import NDArray -import numpy as np from ccc.scipy.stats import rank from ccc.utils import chunker -from numpy.typing import NDArray @njit(cache=True, nogil=True) @@ -244,7 +243,7 @@ def get_feature_parts(params): def get_chunks( - iterable: Union[int, Iterable], n_threads: int, ratio: float = 1 + iterable: int | Iterable, n_threads: int, ratio: float = 1 ) -> Iterable[Iterable[int]]: """ It splits elements in an iterable in chunks according to the number of @@ -341,7 +340,7 @@ def get_n_workers(n_jobs: int | None) -> int: def ccc( x: NDArray, y: NDArray = None, - internal_n_clusters: Union[int, Iterable[int]] = None, + internal_n_clusters: int | Iterable[int] = None, return_parts: bool = False, n_chunks_threads_ratio: int = 1, n_jobs: int = 1, @@ -617,4 +616,4 @@ def ccc( return cm_values, max_parts, parts if pvalue_n_perms is not None and pvalue_n_perms > 0: return cm_values, cm_pvalues - return cm_values \ No newline at end of file + return cm_values diff --git a/libs/ccc/corr.py b/libs/ccc/corr.py index b4bd69cf..6c49a5fa 100644 --- a/libs/ccc/corr.py +++ b/libs/ccc/corr.py @@ -16,9 +16,8 @@ from __future__ import annotations -import pandas as pd - import numpy as np +import pandas as pd from sklearn.metrics import pairwise_distances @@ -60,9 +59,9 @@ def mic(data: pd.DataFrame, estimator="mic_approx", n_jobs=None) -> pd.DataFrame Compute the Maximal Correlation Coefficient (MIC). """ from minepy import pstats + from scipy.spatial.distance import squareform from ccc.methods import mic as mic_single - from scipy.spatial.distance import squareform if n_jobs is None: corr_mat = pstats( @@ -87,9 +86,10 @@ def ccc(data: pd.DataFrame, internal_n_clusters=None, n_jobs=1) -> pd.DataFrame: """ Compute the Clustermatch Correlation Coefficient (CCC). """ - from ccc.coef import ccc from scipy.spatial.distance import squareform + from ccc.coef import ccc + corr_mat = ccc( data.to_numpy(), internal_n_clusters=internal_n_clusters, @@ -110,9 +110,10 @@ def ccc_gpu(data: pd.DataFrame, internal_n_clusters=10, n_jobs=24) -> pd.DataFra """ Compute the Clustermatch Correlation Coefficient (CCC). """ - from ccc.coef.impl_gpu import ccc as ccc_gpu from scipy.spatial.distance import squareform + from ccc.coef.impl_gpu import ccc as ccc_gpu + corr_mat = ccc_gpu( data.to_numpy(), internal_n_clusters=internal_n_clusters, diff --git a/libs/ccc/giant.py b/libs/ccc/giant.py index ce4ff2c2..137c9734 100644 --- a/libs/ccc/giant.py +++ b/libs/ccc/giant.py @@ -2,12 +2,13 @@ Contains functions to interact with the REST API of HumanBase (GIANT networks): https://hb.flatironinstitute.org/ """ -from pathlib import Path -import tempfile + import json +import tempfile +from pathlib import Path -import requests import pandas as pd +import requests URL_GENE_INFO = "https://hb.flatironinstitute.org/api/genes/" URL_TISSUE_PREDICTION = "https://hb.flatironinstitute.org/api/integrations/relevant/" diff --git a/libs/ccc/log.py b/libs/ccc/log.py index d26434a7..b7307a1a 100644 --- a/libs/ccc/log.py +++ b/libs/ccc/log.py @@ -1,8 +1,10 @@ """ Provides logging functions. """ + import logging import logging.config + import yaml from ccc import conf @@ -10,7 +12,7 @@ def _get_logger_config(): """Reads the logging config file in YAML format.""" - with open(conf.GENERAL["LOG_CONFIG_FILE"], "r") as f: + with open(conf.GENERAL["LOG_CONFIG_FILE"]) as f: return yaml.safe_load(f.read()) diff --git a/libs/ccc/methods.py b/libs/ccc/methods.py index bb995640..9eb50f05 100644 --- a/libs/ccc/methods.py +++ b/libs/ccc/methods.py @@ -1,6 +1,7 @@ """ Contains other correlation methods. """ + import warnings from minepy.mine import MINE diff --git a/libs/ccc/plots.py b/libs/ccc/plots.py index 623f3448..73680512 100644 --- a/libs/ccc/plots.py +++ b/libs/ccc/plots.py @@ -39,13 +39,13 @@ from pathlib import Path -import numpy as np -from scipy import stats -import pandas as pd -from IPython.display import display import matplotlib as mpl import matplotlib.pyplot as plt +import numpy as np +import pandas as pd import seaborn as sns +from IPython.display import display +from scipy import stats from seaborn.distributions import _freedman_diaconis_bins from upsetplot import UpSet @@ -249,7 +249,7 @@ def jointplot( c = ccc(x_values, y_values) ax = grid.ax_joint - corr_vals = f"$r$ = {r:.2f}\n" f"$r_s$ = {rs:.2f}\n" f"$c$ = {c:.2f}" + corr_vals = f"$r$ = {r:.2f}\n$r_s$ = {rs:.2f}\n$c$ = {c:.2f}" bbox = dict(boxstyle="round", fc="white", ec="black", alpha=0.15) ax.text( 0.25, diff --git a/libs/ccc/pytorch/core.py b/libs/ccc/pytorch/core.py index a7776606..212b359d 100644 --- a/libs/ccc/pytorch/core.py +++ b/libs/ccc/pytorch/core.py @@ -78,6 +78,7 @@ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """ + import numpy as np from numba import njit diff --git a/libs/ccc/scipy/stats.py b/libs/ccc/scipy/stats.py index 69f65ac9..7e1ae916 100644 --- a/libs/ccc/scipy/stats.py +++ b/libs/ccc/scipy/stats.py @@ -57,6 +57,7 @@ (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """ + import numpy as np from numba import njit from numpy.typing import NDArray diff --git a/libs/ccc/sklearn/metrics.py b/libs/ccc/sklearn/metrics.py index efb1ca4a..690bab76 100644 --- a/libs/ccc/sklearn/metrics.py +++ b/libs/ccc/sklearn/metrics.py @@ -38,9 +38,8 @@ from __future__ import annotations -from numba import njit - import numpy as np +from numba import njit @njit(cache=True, nogil=True) diff --git a/libs/ccc/utils/utility_functions.py b/libs/ccc/utils/utility_functions.py index 7b5b2b36..f2ef1e8a 100644 --- a/libs/ccc/utils/utility_functions.py +++ b/libs/ccc/utils/utility_functions.py @@ -1,8 +1,9 @@ """ General utility functions. """ -import re + import hashlib +import re from pathlib import Path from subprocess import run @@ -122,13 +123,13 @@ def human_format(num): Formats numbers with a shortened style. Taken from: https://stackoverflow.com/a/45846841 """ - num = float("{:.3g}".format(num)) + num = float(f"{num:.3g}") magnitude = 0 while abs(num) >= 1000: magnitude += 1 num /= 1000.0 return "{}{}".format( - "{:f}".format(num).rstrip("0").rstrip("."), ["", "K", "M", "B", "T"][magnitude] + f"{num:f}".rstrip("0").rstrip("."), ["", "K", "M", "B", "T"][magnitude] ) diff --git a/libs/ccc_cuda_ext/binder.cu b/libs/ccc_cuda_ext/binder.cu index fe1e53c9..ae9501a8 100644 --- a/libs/ccc_cuda_ext/binder.cu +++ b/libs/ccc_cuda_ext/binder.cu @@ -1,8 +1,8 @@ #include #include -#include "metrics.cuh" #include "coef.cuh" +#include "metrics.cuh" namespace py = pybind11; @@ -11,10 +11,10 @@ using namespace pybind11::literals; PYBIND11_MODULE(ccc_cuda_ext, m) { m.doc() = "CUDA extension module for CCC"; - m.def("ari_int32", &ari, "CUDA version of Adjusted Rand Index (ARI) calculation", - "parts"_a, "n_features"_a, "n_parts"_a, "n_objs"_a, "batch_start"_a = 0, "batch_size"_a = 0); // _a is a shorter notation for named arguments - m.def("compute_coef", &compute_coef, - "CUDA version of CCC coefficient calculation", + m.def("ari_int32", &ari, "CUDA version of Adjusted Rand Index (ARI) calculation", "parts"_a, "n_features"_a, + "n_parts"_a, "n_objs"_a, "batch_start"_a = 0, + "batch_size"_a = 0); // _a is a shorter notation for named arguments + m.def("compute_coef", &compute_coef, "CUDA version of CCC coefficient calculation", "parts"_a, // numpy array of partitions "n_features"_a, // number of features "n_parts"_a, // number of partitions per feature diff --git a/libs/ccc_cuda_ext/coef.cu b/libs/ccc_cuda_ext/coef.cu index fa25134d..372b089b 100644 --- a/libs/ccc_cuda_ext/coef.cu +++ b/libs/ccc_cuda_ext/coef.cu @@ -1,33 +1,32 @@ -#include #include +#include +#include +#include #include -#include -#include -#include #include -#include #include -#include -#include +#include +#include +#include +#include +#include +#include #include -#include #include +#include #include #include -#include -#include -#include #include #include +#include #include "coef.cuh" -#include "metrics.cuh" #include "math.cuh" +#include "metrics.cuh" #include "utils.cuh" namespace py = pybind11; - /** * @brief CUDA kernel to find maximum ARI values and their corresponding partition pairs * @@ -43,10 +42,7 @@ namespace py = pybind11; * @param reduction_range Number of partition pairs to process per feature comparison */ template -__global__ void findMaxAriKernel(const T *aris, - uint8_t *max_parts, - T *cm_values, - const int n_partitions, +__global__ void findMaxAriKernel(const T *aris, uint8_t *max_parts, T *cm_values, const int n_partitions, const int reduction_range) { /* @@ -67,8 +63,8 @@ __global__ void findMaxAriKernel(const T *aris, */ typedef cub::KeyValuePair KeyValuePairT; KeyValuePairT thread_data; - thread_data.key = -1; // Initialize to invalid index - thread_data.value = -1.0f; // Initialize to very small value + thread_data.key = -1; // Initialize to invalid index + thread_data.value = -1.0f; // Initialize to very small value bool has_nan = false; /* @@ -92,7 +88,7 @@ __global__ void findMaxAriKernel(const T *aris, if (val > thread_data.value) { thread_data.value = val; - thread_data.key = i; // Store the local index + thread_data.key = i; // Store the local index } } @@ -108,13 +104,13 @@ __global__ void findMaxAriKernel(const T *aris, block_has_nan = false; } __syncthreads(); - + if (has_nan) { block_has_nan = true; } __syncthreads(); - + if (block_has_nan) { if (threadIdx.x == 0) @@ -135,10 +131,7 @@ __global__ void findMaxAriKernel(const T *aris, __shared__ typename BlockReduceT::TempStorage temp_storage; // Use standard ArgMax but rely on CUB's tie-breaking behavior - KeyValuePairT aggregate = BlockReduceT(temp_storage).Reduce( - thread_data, - cub::ArgMax() - ); + KeyValuePairT aggregate = BlockReduceT(temp_storage).Reduce(thread_data, cub::ArgMax()); /* * Result Writing @@ -194,21 +187,21 @@ __global__ void initRandomStates(curandState *states, const uint32_t n_states, c * @param n_perms Number of permutations * @param n_objects Number of objects to permute */ -__global__ void generatePermutations(curandState *states, uint32_t *perm_indices, - const uint32_t n_perms, const uint32_t n_objects) +__global__ void generatePermutations(curandState *states, uint32_t *perm_indices, const uint32_t n_perms, + const uint32_t n_objects) { uint32_t perm_idx = blockIdx.x * blockDim.x + threadIdx.x; if (perm_idx < n_perms) { curandState *state = &states[perm_idx]; uint32_t *perm = &perm_indices[perm_idx * n_objects]; - + // Initialize sequential indices for (uint32_t i = 0; i < n_objects; ++i) { perm[i] = i; } - + // Fisher-Yates shuffle for (uint32_t i = n_objects - 1; i > 0; --i) { @@ -240,21 +233,27 @@ __global__ void generatePermutations(curandState *states, uint32_t *perm_indices * @return ARI value */ template -__device__ R computePermutedARI(const T *part_i, const T *part_j, - const uint32_t *perm, const uint32_t n_objects, +__device__ R computePermutedARI(const T *part_i, const T *part_j, const uint32_t *perm, const uint32_t n_objects, const int k, int *scratch) { // Quick validation - if (n_objects == 0) return 0.0f; - if (k <= 0) return 0.0f; + if (n_objects == 0) + return 0.0f; + if (k <= 0) + return 0.0f; int *contingency = scratch; // k * k ints int *sum_rows = contingency + k * k; // k ints int *sum_cols = sum_rows + k; // k ints // Initialize scratch arrays - for (int i = 0; i < k * k; ++i) contingency[i] = 0; - for (int i = 0; i < k; ++i) { sum_rows[i] = 0; sum_cols[i] = 0; } + for (int i = 0; i < k * k; ++i) + contingency[i] = 0; + for (int i = 0; i < k; ++i) + { + sum_rows[i] = 0; + sum_cols[i] = 0; + } // Build contingency matrix and compute sum of squares in a single pass. long long sum_squares = 0; @@ -343,14 +342,14 @@ __device__ R computePermutedARI(const T *part_i, const T *part_j, * @param scratch_stride Ints per thread: k*k + 2*k */ template -__global__ void computePermutationCCC(const T *parts_i, const T *parts_j, - const uint32_t *perm_indices, R *perm_ccc_values, - const uint32_t perm_offset, const uint32_t perm_count, - const uint32_t n_partitions, const uint32_t n_objects, - const int k, int *scratch, const uint64_t scratch_stride) +__global__ void computePermutationCCC(const T *parts_i, const T *parts_j, const uint32_t *perm_indices, + R *perm_ccc_values, const uint32_t perm_offset, const uint32_t perm_count, + const uint32_t n_partitions, const uint32_t n_objects, const int k, int *scratch, + const uint64_t scratch_stride) { const uint32_t local_idx = blockIdx.x * blockDim.x + threadIdx.x; - if (local_idx >= perm_count) return; + if (local_idx >= perm_count) + return; const uint32_t perm_idx = perm_offset + local_idx; // Permutation for this thread. Scratch is indexed by the LOCAL id so the @@ -430,8 +429,8 @@ __global__ void computePermutationCCC(const T *parts_i, const T *parts_j, * @param n_perms Number of permutations per comparison */ template -__global__ void computePValues(const R *perm_ccc_values, const R *observed_ccc_values, - R *pvalues, const uint64_t n_comparisons, const uint32_t n_perms) +__global__ void computePValues(const R *perm_ccc_values, const R *observed_ccc_values, R *pvalues, + const uint64_t n_comparisons, const uint32_t n_perms) { const uint64_t comp_idx = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; if (comp_idx < n_comparisons) @@ -464,7 +463,8 @@ void check_feature_comp_bounds(const size_t n_features) { if (n_features > 1 && n_features > UINT64_MAX / (n_features - 1)) { - throw std::range_error("Feature comparison count would exceed maximum representable value: n_features too large"); + throw std::range_error( + "Feature comparison count would exceed maximum representable value: n_features too large"); } } @@ -479,12 +479,14 @@ void check_ari_count_bounds(const uint64_t n_feature_comp, const size_t n_partit { if (n_feature_comp > UINT64_MAX / n_partitions) { - throw std::range_error("ARI count would exceed maximum representable value: n_feature_comp * n_partitions too large"); + throw std::range_error( + "ARI count would exceed maximum representable value: n_feature_comp * n_partitions too large"); } const uint64_t temp = n_feature_comp * n_partitions; if (temp > UINT64_MAX / n_partitions) { - throw std::range_error("ARI count would exceed maximum representable value: n_feature_comp * n_partitions * n_partitions too large"); + throw std::range_error("ARI count would exceed maximum representable value: n_feature_comp * n_partitions * " + "n_partitions too large"); } } @@ -516,12 +518,8 @@ uint64_t calculate_total_aris(const uint64_t n_feature_comp, const size_t n_part } template -auto compute_coef(const py::array_t &parts, - const size_t n_features, - const size_t n_partitions, - const size_t n_objects, - const bool return_parts, - std::optional pvalue_n_perms) -> py::object +auto compute_coef(const py::array_t &parts, const size_t n_features, const size_t n_partitions, + const size_t n_objects, const bool return_parts, std::optional pvalue_n_perms) -> py::object { /* * Input validation (before any device work) @@ -533,24 +531,21 @@ auto compute_coef(const py::array_t &parts, py::buffer_info buffer = parts.request(); if (buffer.format != py::format_descriptor::format()) { - throw py::value_error( - std::string("Partitions array has an incompatible dtype: expected numpy format '") + - py::format_descriptor::format() + "', got '" + buffer.format + "'"); + throw py::value_error(std::string("Partitions array has an incompatible dtype: expected numpy format '") + + py::format_descriptor::format() + "', got '" + buffer.format + "'"); } - if (buffer.ndim != 3 || - buffer.shape[0] != static_cast(n_features) || + if (buffer.ndim != 3 || buffer.shape[0] != static_cast(n_features) || buffer.shape[1] != static_cast(n_partitions) || buffer.shape[2] != static_cast(n_objects)) { std::string got = buffer.ndim == 3 - ? ("(" + std::to_string(buffer.shape[0]) + ", " + - std::to_string(buffer.shape[1]) + ", " + - std::to_string(buffer.shape[2]) + ")") + ? ("(" + std::to_string(buffer.shape[0]) + ", " + std::to_string(buffer.shape[1]) + + ", " + std::to_string(buffer.shape[2]) + ")") : ("ndim=" + std::to_string(buffer.ndim)); throw py::value_error( "Partitions array shape mismatch: expected (n_features, n_partitions, n_objects) = (" + - std::to_string(n_features) + ", " + std::to_string(n_partitions) + ", " + - std::to_string(n_objects) + ") but got " + got); + std::to_string(n_features) + ", " + std::to_string(n_partitions) + ", " + std::to_string(n_objects) + + ") but got " + got); } } if (n_partitions == 0) @@ -566,15 +561,18 @@ auto compute_coef(const py::array_t &parts, } // Check for CCC_GPU_LOGGING environment variable to enable debug logging - const char* logging_env = std::getenv("CCC_GPU_LOGGING"); - if (logging_env != nullptr) { + const char *logging_env = std::getenv("CCC_GPU_LOGGING"); + if (logging_env != nullptr) + { spdlog::set_level(spdlog::level::debug); // Check CUDA info spdlog::debug("CUDA Device Info:"); print_cuda_device_info(); spdlog::debug("CUDA Memory Info:"); print_cuda_memory_info(); - } else { + } + else + { // Disable debug logging by default spdlog::set_level(spdlog::level::err); } @@ -665,8 +663,7 @@ auto compute_coef(const py::array_t &parts, */ for (uint64_t batch_start = 0; batch_start < n_aris; batch_start += batch_n_aris) { - spdlog::debug("Processing batch {} of {}", - (batch_start / batch_n_aris + 1), + spdlog::debug("Processing batch {} of {}", (batch_start / batch_n_aris + 1), (n_aris + batch_n_aris - 1) / batch_n_aris); spdlog::debug(" Start index: {}", batch_start); spdlog::debug(" Batch size: {}", batch_n_aris); @@ -678,22 +675,18 @@ auto compute_coef(const py::array_t &parts, spdlog::debug(" Current batch size: {}", current_batch_size); // Compute ARIs for this batch - const auto d_aris = ari_core_device( - parts, n_features, n_partitions, n_objects, batch_start, current_batch_size); + const auto d_aris = + ari_core_device(parts, n_features, n_partitions, n_objects, batch_start, current_batch_size); // Configure kernel launch parameters const int threadsPerBlock = 128; const int numBlocks = current_batch_size / (n_partitions * n_partitions); - spdlog::debug(" Launching reduction kernel with {} blocks, {} threads per block", - numBlocks, threadsPerBlock); + spdlog::debug(" Launching reduction kernel with {} blocks, {} threads per block", numBlocks, threadsPerBlock); // Launch kernel to find maximum values and their partition pairs findMaxAriKernel<<>>( - thrust::raw_pointer_cast(d_aris->data()), - thrust::raw_pointer_cast(d_max_parts.data()), - thrust::raw_pointer_cast(d_cm_values.data()), - n_partitions, - reduction_range); + thrust::raw_pointer_cast(d_aris->data()), thrust::raw_pointer_cast(d_max_parts.data()), + thrust::raw_pointer_cast(d_cm_values.data()), n_partitions, reduction_range); // Check for kernel errors cudaError_t kernelError = cudaGetLastError(); @@ -712,7 +705,8 @@ auto compute_coef(const py::array_t &parts, // Copy batch results back to host thrust::copy(d_cm_values.begin(), d_cm_values.begin() + current_batch_size / (n_partitions * n_partitions), batch_cm_values.begin()); - thrust::copy(d_max_parts.begin(), d_max_parts.begin() + (current_batch_size / (n_partitions * n_partitions)) * 2, + thrust::copy(d_max_parts.begin(), + d_max_parts.begin() + (current_batch_size / (n_partitions * n_partitions)) * 2, batch_max_parts.begin()); // Update main result arrays with batch results @@ -752,7 +746,8 @@ auto compute_coef(const py::array_t &parts, // removes the old MAX_CLUSTERS=16 cliff that silently returned ARI 0.0. int k_global = static_cast( thrust::reduce(d_parts.begin(), d_parts.end(), static_cast(-1), thrust::maximum()) + 1); - if (k_global < 1) k_global = 1; + if (k_global < 1) + k_global = 1; const uint64_t scratch_stride = static_cast(k_global) * k_global + 2ULL * k_global; // Per-permutation contingency/sum scratch. Its size scales as k^2 per @@ -765,9 +760,9 @@ auto compute_coef(const py::array_t &parts, const uint64_t scratch_elem_bytes = std::max(scratch_stride * sizeof(int), 1); const uint64_t scratch_budget = std::min(free_scratch / 4, static_cast(1) << 30); // cap at 1 GiB - uint32_t perm_batch = - static_cast(std::max(1, scratch_budget / scratch_elem_bytes)); - if (perm_batch > n_perms) perm_batch = n_perms; + uint32_t perm_batch = static_cast(std::max(1, scratch_budget / scratch_elem_bytes)); + if (perm_batch > n_perms) + perm_batch = n_perms; thrust::device_vector d_perm_scratch(static_cast(perm_batch) * scratch_stride); // Allocate device memory for p-value computation @@ -780,20 +775,14 @@ auto compute_coef(const py::array_t &parts, const uint32_t block_size = 256; const uint32_t grid_size_states = (n_perms + block_size - 1) / block_size; - initRandomStates<<>>( - thrust::raw_pointer_cast(d_rand_states.data()), - n_perms, - rand_seed - ); + initRandomStates<<>>(thrust::raw_pointer_cast(d_rand_states.data()), n_perms, + rand_seed); CUDA_CHECK_KERNEL("initRandomStates"); // Generate permutation indices - generatePermutations<<>>( - thrust::raw_pointer_cast(d_rand_states.data()), - thrust::raw_pointer_cast(d_perm_indices.data()), - n_perms, - n_objects - ); + generatePermutations<<>>(thrust::raw_pointer_cast(d_rand_states.data()), + thrust::raw_pointer_cast(d_perm_indices.data()), n_perms, + n_objects); CUDA_CHECK_KERNEL("generatePermutations"); /* @@ -809,8 +798,10 @@ auto compute_coef(const py::array_t &parts, const size_t perm_value_bytes = std::max(static_cast(n_perms) * sizeof(R), 1); const size_t budget = std::min(free_mem / 4, static_cast(512) * 1024 * 1024); uint64_t chunk_comps = static_cast(budget / perm_value_bytes); - if (chunk_comps < 1) chunk_comps = 1; - if (chunk_comps > n_feature_comp) chunk_comps = n_feature_comp; + if (chunk_comps < 1) + chunk_comps = 1; + if (chunk_comps > n_feature_comp) + chunk_comps = n_feature_comp; thrust::device_vector d_perm_ccc_values(static_cast(chunk_comps) * n_perms); @@ -830,39 +821,41 @@ auto compute_coef(const py::array_t &parts, uint64_t feat_i, feat_j; get_coords_from_index(static_cast(n_features), comp_idx, feat_i, feat_j); - R *out_slice = thrust::raw_pointer_cast(d_perm_ccc_values.data()) + - local * static_cast(n_perms); + R *out_slice = + thrust::raw_pointer_cast(d_perm_ccc_values.data()) + local * static_cast(n_perms); // Validate feature indices to prevent memory corruption. This // should not trigger for a valid comp_idx; fill the slice so the // stale reuse of the batch buffer cannot corrupt the p-value. if (feat_i >= n_features || feat_j >= n_features) { - spdlog::error("Invalid feature indices: feat_i={}, feat_j={}, n_features={}", - feat_i, feat_j, n_features); + spdlog::error("Invalid feature indices: feat_i={}, feat_j={}, n_features={}", feat_i, feat_j, + n_features); thrust::fill(d_perm_ccc_values.begin() + local * static_cast(n_perms), d_perm_ccc_values.begin() + (local + 1) * static_cast(n_perms), static_cast(0)); continue; } - const T* d_parts_i = thrust::raw_pointer_cast(d_parts.data()) + feat_i * n_partitions * n_objects; - const T* d_parts_j = thrust::raw_pointer_cast(d_parts.data()) + feat_j * n_partitions * n_objects; + const T *d_parts_i = thrust::raw_pointer_cast(d_parts.data()) + feat_i * n_partitions * n_objects; + const T *d_parts_j = thrust::raw_pointer_cast(d_parts.data()) + feat_j * n_partitions * n_objects; // Count valid partitions on host (small operation) uint32_t valid_count_i = 0, valid_count_j = 0; - const T* host_parts_i = parts.data() + feat_i * n_partitions * n_objects; - const T* host_parts_j = parts.data() + feat_j * n_partitions * n_objects; + const T *host_parts_i = parts.data() + feat_i * n_partitions * n_objects; + const T *host_parts_j = parts.data() + feat_j * n_partitions * n_objects; for (uint32_t p = 0; p < n_partitions; ++p) { - if (host_parts_i[p * n_objects] >= 0) valid_count_i++; - if (host_parts_j[p * n_objects] >= 0) valid_count_j++; + if (host_parts_i[p * n_objects] >= 0) + valid_count_i++; + if (host_parts_j[p * n_objects] >= 0) + valid_count_j++; } // Permute the feature that generated MORE valid partitions, // matching the CPU reference (impl.py compute_coef). - const T* d_parts_to_permute = (valid_count_i > valid_count_j) ? d_parts_i : d_parts_j; - const T* d_parts_fixed = (valid_count_i > valid_count_j) ? d_parts_j : d_parts_i; + const T *d_parts_to_permute = (valid_count_i > valid_count_j) ? d_parts_i : d_parts_j; + const T *d_parts_fixed = (valid_count_i > valid_count_j) ? d_parts_j : d_parts_i; // Process permutations in sub-batches bounded by the scratch budget. for (uint32_t p_off = 0; p_off < n_perms; p_off += perm_batch) @@ -870,18 +863,9 @@ auto compute_coef(const py::array_t &parts, const uint32_t p_cnt = std::min(perm_batch, n_perms - p_off); const uint32_t p_grid = (p_cnt + block_size - 1) / block_size; computePermutationCCC<<>>( - d_parts_fixed, - d_parts_to_permute, - thrust::raw_pointer_cast(d_perm_indices.data()), - out_slice, - p_off, - p_cnt, - n_partitions, - n_objects, - k_global, - thrust::raw_pointer_cast(d_perm_scratch.data()), - scratch_stride - ); + d_parts_fixed, d_parts_to_permute, thrust::raw_pointer_cast(d_perm_indices.data()), out_slice, + p_off, p_cnt, n_partitions, n_objects, k_global, + thrust::raw_pointer_cast(d_perm_scratch.data()), scratch_stride); CUDA_CHECK_KERNEL("computePermutationCCC"); } } @@ -891,10 +875,7 @@ auto compute_coef(const py::array_t &parts, computePValues<<>>( thrust::raw_pointer_cast(d_perm_ccc_values.data()), thrust::raw_pointer_cast(d_observed_ccc_values.data()) + chunk_start, - thrust::raw_pointer_cast(d_computed_pvalues.data()) + chunk_start, - chunk_len, - n_perms - ); + thrust::raw_pointer_cast(d_computed_pvalues.data()) + chunk_start, chunk_len, n_perms); CUDA_CHECK_KERNEL("computePValues"); } @@ -922,12 +903,10 @@ auto compute_coef(const py::array_t &parts, const auto cm_pvalues_py = pvalue_n_perms.has_value() ? py::object(py::array_t(cm_pvalues.size(), cm_pvalues.data())) : py::object(py::none()); - const auto max_parts_py = py::array_t(max_parts.size(), max_parts.data()).reshape({n_feature_comp, static_cast(2)}); + const auto max_parts_py = + py::array_t(max_parts.size(), max_parts.data()).reshape({n_feature_comp, static_cast(2)}); - return py::make_tuple( - cm_values_py, - cm_pvalues_py, - max_parts_py); + return py::make_tuple(cm_values_py, cm_pvalues_py, max_parts_py); } // Below is the explicit instantiation of the ari template function. @@ -937,8 +916,6 @@ auto compute_coef(const py::array_t &parts, // implementation of the template functions, we need to explicitly instantiate them here, so that they can be picked up // by the linker. template auto compute_coef(const py::array_t &parts, - const size_t n_features, - const size_t n_partitions, - const size_t n_objects, - const bool return_parts, - std::optional pvalue_n_perms) -> py::object; + const size_t n_features, const size_t n_partitions, const size_t n_objects, + const bool return_parts, std::optional pvalue_n_perms) + -> py::object; diff --git a/libs/ccc_cuda_ext/coef.cuh b/libs/ccc_cuda_ext/coef.cuh index 04e95731..987b43e1 100644 --- a/libs/ccc_cuda_ext/coef.cuh +++ b/libs/ccc_cuda_ext/coef.cuh @@ -1,14 +1,11 @@ #pragma once -#include #include +#include namespace py = pybind11; template -auto compute_coef(const py::array_t &parts, - const size_t n_features, - const size_t n_parts, - const size_t n_objs, - const bool return_parts, - std::optional pvalue_n_perms = std::nullopt) -> py::object; +auto compute_coef(const py::array_t &parts, const size_t n_features, const size_t n_parts, + const size_t n_objs, const bool return_parts, std::optional pvalue_n_perms = std::nullopt) + -> py::object; diff --git a/libs/ccc_cuda_ext/math.cuh b/libs/ccc_cuda_ext/math.cuh index 0ff1f117..db30e8fa 100644 --- a/libs/ccc_cuda_ext/math.cuh +++ b/libs/ccc_cuda_ext/math.cuh @@ -11,8 +11,8 @@ #pragma once -#include #include +#include /** * @brief Unravel a flat (linear) index into corresponding 2D indices @@ -33,8 +33,8 @@ * int row, col; * unravel_index(5, 3, &row, &col); // For a 3-column array, index 5 -> row=1, col=2 */ -__device__ __host__ inline void unravel_index(unsigned int flat_idx, unsigned int num_cols, - unsigned int &row, unsigned int &col) +__device__ __host__ inline void unravel_index(unsigned int flat_idx, unsigned int num_cols, unsigned int &row, + unsigned int &col) { // change int to uint32_t row = flat_idx / num_cols; // Compute row index @@ -65,49 +65,52 @@ __device__ __host__ inline void unravel_index(unsigned int flat_idx, unsigned in * [ 1 3 - 5 ] For idx=3, the function returns x=1, y=2 * [ 2 4 5 - ] representing the position in the original matrix */ -__device__ __host__ inline void get_coords_from_index(unsigned int n_obj, unsigned int idx, - unsigned int &x, unsigned int &y) +__device__ __host__ inline void get_coords_from_index(unsigned int n_obj, unsigned int idx, unsigned int &x, + unsigned int &y) { // Prevent overflow by using int64_t for intermediate calculations int64_t n_obj_64 = static_cast(n_obj); int64_t idx_64 = static_cast(idx); - + // Calculate 'b' based on the input n_obj int64_t b = 1 - 2 * n_obj_64; - + // Calculate discriminant using double precision to avoid overflow double discriminant = static_cast(b * b) - 8.0 * static_cast(idx_64); - + // Check for negative discriminant (invalid input) - if (discriminant < 0) { + if (discriminant < 0) + { x = 0; y = 0; return; } - + // Calculate 'x' using the quadratic formula part double x_float = (-static_cast(b) - sqrt(discriminant)) / 2.0; int64_t x_64 = static_cast(floor(x_float)); - + // Bounds checking for x - if (x_64 < 0 || x_64 > UINT32_MAX) { + if (x_64 < 0 || x_64 > UINT32_MAX) + { x = 0; y = 0; return; } - + x = static_cast(x_64); - + // Calculate 'y' based on 'x' and the index int64_t y_64 = idx_64 + x_64 * (b + x_64 + 2) / 2 + 1; - + // Bounds checking for y - if (y_64 < 0 || y_64 > UINT32_MAX) { + if (y_64 < 0 || y_64 > UINT32_MAX) + { x = 0; y = 0; return; } - + y = static_cast(y_64); } @@ -125,8 +128,7 @@ __device__ __host__ inline void get_coords_from_index(unsigned int n_obj, unsign * @param[out] x Reference to store the calculated row coordinate. * @param[out] y Reference to store the calculated column coordinate. */ -__device__ __host__ inline void get_coords_from_index(uint64_t n_obj, uint64_t idx, - uint64_t &x, uint64_t &y) +__device__ __host__ inline void get_coords_from_index(uint64_t n_obj, uint64_t idx, uint64_t &x, uint64_t &y) { const int64_t n_obj_64 = static_cast(n_obj); const int64_t idx_64 = static_cast(idx); @@ -134,8 +136,7 @@ __device__ __host__ inline void get_coords_from_index(uint64_t n_obj, uint64_t i const int64_t b = 1 - 2 * n_obj_64; // Use double precision for the discriminant to avoid intermediate overflow. - const double discriminant = - static_cast(b) * static_cast(b) - 8.0 * static_cast(idx_64); + const double discriminant = static_cast(b) * static_cast(b) - 8.0 * static_cast(idx_64); if (discriminant < 0) { diff --git a/libs/ccc_cuda_ext/metrics.cu b/libs/ccc_cuda_ext/metrics.cu index b800ae2c..a5d3b0a1 100644 --- a/libs/ccc_cuda_ext/metrics.cu +++ b/libs/ccc_cuda_ext/metrics.cu @@ -1,19 +1,19 @@ -#include #include +#include #include #include -#include #include -#include #include +#include +#include -#include -#include -#include -#include #include "metrics.cuh" #include "utils.cuh" +#include +#include +#include +#include namespace py = pybind11; @@ -125,7 +125,6 @@ __device__ void get_contingency_matrix_shared(T *part0, T *part1, int n_objs, in shared_cont_mat[i] = 0; } __syncthreads(); - #pragma unroll for (int i = tid; i < n_objs; i += n_block_threads) @@ -136,10 +135,11 @@ __device__ void get_contingency_matrix_shared(T *part0, T *part1, int n_objs, in const T col = part1[i]; // Bounds checking to ensure valid array access - if (row < 0 || row >= k || col < 0 || col >= k) { + if (row < 0 || row >= k || col < 0 || col >= k) + { continue; // Skip invalid indices } - + atomicAdd(&shared_cont_mat[row * k + col], 1); } __syncthreads(); @@ -175,18 +175,17 @@ __device__ void get_contingency_matrix_global(T *part0, T *part1, int n_objs, in const T col = part1[i]; // Add bounds checking - if (row < 0 || row >= k || col < 0 || col >= k) { + if (row < 0 || row >= k || col < 0 || col >= k) + { continue; // Skip invalid values } - + // Use atomic operations since we're writing to global memory atomicAdd(&global_cont_mat[row * k + col], 1); } __syncthreads(); - } - /** * @brief CUDA device function to compute the pair confusion matrix * @param[in] contingency Pointer to the contingency matrix @@ -196,13 +195,8 @@ __device__ void get_contingency_matrix_global(T *part0, T *part1, int n_objs, in * @param[in] k Number of clusters (assuming k is the max of clusters in part0 and part1) * @param[out] C Pointer to the output pair confusion matrix (2x2) */ -__device__ void get_pair_confusion_matrix( - const int *__restrict__ contingency, - int *sum_rows, - int *sum_cols, - const int n_objs, - const int k, - long long *C) +__device__ void get_pair_confusion_matrix(const int *__restrict__ contingency, int *sum_rows, int *sum_cols, + const int n_objs, const int k, long long *C) { // TODO: use block-level reduction @@ -267,7 +261,6 @@ __device__ void get_pair_confusion_matrix( C[2] = temp - sum_squares; // C[1,0] C[0] = (long long)n_objs * (long long)n_objs - C[1] - C[2] - sum_squares; // C[0,0] - } } @@ -285,17 +278,10 @@ __device__ void get_pair_confusion_matrix( * @param out Output array of ARIs */ template -__global__ void ari_kernel_global(T *parts, - const uint64_t n_aris, - const uint64_t n_features, - const uint64_t n_parts, - const uint64_t n_objs, - const uint64_t n_elems_per_feat, - const uint64_t n_part_mat_elems, - const uint32_t k, - const uint64_t batch_start, - int *global_cont_matrices, - float *out) +__global__ void ari_kernel_global(T *parts, const uint64_t n_aris, const uint64_t n_features, const uint64_t n_parts, + const uint64_t n_objs, const uint64_t n_elems_per_feat, + const uint64_t n_part_mat_elems, const uint32_t k, const uint64_t batch_start, + int *global_cont_matrices, float *out) { /* * Step 0: Compute global memory addresses for this block @@ -304,21 +290,21 @@ __global__ void ari_kernel_global(T *parts, const uint64_t cont_mat_size = k * k; const uint64_t sum_arrays_size = 2 * k; const uint64_t pair_confusion_size = 4; - + // Each block gets its own section of global memory // Calculate offsets to ensure proper alignment for long long const uint64_t ints_per_block = cont_mat_size + sum_arrays_size; - const uint64_t aligned_ints_per_block = ((ints_per_block + 1) / 2) * 2; // Align to 8-byte boundary - const uint64_t offset_per_block_bytes = aligned_ints_per_block * sizeof(int) + pair_confusion_size * sizeof(long long); + const uint64_t aligned_ints_per_block = ((ints_per_block + 1) / 2) * 2; // Align to 8-byte boundary + const uint64_t offset_per_block_bytes = + aligned_ints_per_block * sizeof(int) + pair_confusion_size * sizeof(long long); const uint64_t offset_per_block_ints = offset_per_block_bytes / sizeof(int); - + int *g_contingency = global_cont_matrices + block_id * offset_per_block_ints; int *g_sum_rows = g_contingency + cont_mat_size; int *g_sum_cols = g_sum_rows + k; // Ensure 8-byte alignment for long long int *aligned_start = g_contingency + aligned_ints_per_block; - long long *g_pair_confusion_matrix = (long long*)aligned_start; - + long long *g_pair_confusion_matrix = (long long *)aligned_start; /* * Step 1: Each thread unravels flat indices and loads the corresponding data @@ -329,10 +315,10 @@ __global__ void ari_kernel_global(T *parts, uint64_t i, j; get_coords_from_index(n_features, feature_comp_flat_idx, &i, &j); - + uint64_t m, n; unravel_index(part_pair_flat_idx, n_parts, &m, &n); - + T *t_data_part0 = parts + i * n_elems_per_feat + m * n_objs; T *t_data_part1 = parts + j * n_elems_per_feat + n * n_objs; @@ -372,8 +358,7 @@ __global__ void ari_kernel_global(T *parts, long long fn = g_pair_confusion_matrix[2]; long long tp = g_pair_confusion_matrix[3]; float ari = 0.0f; - - + if (fn == 0 && fp == 0) { ari = 1.0f; @@ -381,9 +366,9 @@ __global__ void ari_kernel_global(T *parts, else { double numerator = 2.0 * ((double)tp * (double)tn - (double)fn * (double)fp); - double denominator = ((double)tp + (double)fn) * ((double)fn + (double)tn) + ((double)tp + (double)fp) * ((double)fp + (double)tn); + double denominator = ((double)tp + (double)fn) * ((double)fn + (double)tn) + + ((double)tp + (double)fp) * ((double)fp + (double)tn); ari = (float)(numerator / denominator); - } out[blockIdx.x] = ari; } @@ -404,29 +389,22 @@ __global__ void ari_kernel_global(T *parts, */ // TODO: Parameterize the int type to allow using narrower int types for memory efficiency template -__global__ void ari_kernel(T *parts, - const uint64_t n_aris, - const uint64_t n_features, - const uint64_t n_parts, - const uint64_t n_objs, - const uint64_t n_elems_per_feat, - const uint64_t n_part_mat_elems, - const uint32_t k, - const uint64_t batch_start, - float *out) +__global__ void ari_kernel(T *parts, const uint64_t n_aris, const uint64_t n_features, const uint64_t n_parts, + const uint64_t n_objs, const uint64_t n_elems_per_feat, const uint64_t n_part_mat_elems, + const uint32_t k, const uint64_t batch_start, float *out) { /* * Step 0: Compute shared memory addresses */ extern __shared__ int shared_mem[]; - int *s_contingency = shared_mem; // k * k elements - int *s_sum_rows = s_contingency + (k * k); // k elements - int *s_sum_cols = s_sum_rows + k; // k elements + int *s_contingency = shared_mem; // k * k elements + int *s_sum_rows = s_contingency + (k * k); // k elements + int *s_sum_cols = s_sum_rows + k; // k elements // Align to 8-byte boundary for long long // Calculate offset in int units, then align to 8-byte boundary int offset_ints = k * k + 2 * k; - int aligned_offset_ints = ((offset_ints + 1) / 2) * 2; // Align to 8-byte (2 int) boundary - long long *s_pair_confusion_matrix = (long long*)(shared_mem + aligned_offset_ints); // 4 elements + int aligned_offset_ints = ((offset_ints + 1) / 2) * 2; // Align to 8-byte (2 int) boundary + long long *s_pair_confusion_matrix = (long long *)(shared_mem + aligned_offset_ints); // 4 elements /* * Step 1: Each thead, unravel flat indices and load the corresponding data into shared memory @@ -435,7 +413,8 @@ __global__ void ari_kernel(T *parts, const uint64_t ari_block_idx = blockIdx.x + batch_start; // obtain the corresponding parts and unique counts uint64_t feature_comp_flat_idx = ari_block_idx / n_part_mat_elems; // flat comparison pair index for two features - uint64_t part_pair_flat_idx = ari_block_idx % n_part_mat_elems; // flat comparison pair index for two partitions of one feature pair + uint64_t part_pair_flat_idx = + ari_block_idx % n_part_mat_elems; // flat comparison pair index for two partitions of one feature pair uint64_t i, j; // Unravel the feature indices @@ -450,9 +429,9 @@ __global__ void ari_kernel(T *parts, // Unravel the partition indices within the feature pair // For example, if n_parts = 3, n_part_mat_elems = n_parts * n_parts = 9 - // The partition indices of the pair being compared are (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2) - // i.e., the pairs being compared are part0-part1, part0-part2, part1-part0, part1-part1, part1-part2, part2-part0, part2-part1, part2-part2 - // The range of the flattened index is [0, n_part_mat_elems - 1] = [0, 8] + // The partition indices of the pair being compared are (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, + // 2) i.e., the pairs being compared are part0-part1, part0-part2, part1-part0, part1-part1, part1-part2, + // part2-part0, part2-part1, part2-part2 The range of the flattened index is [0, n_part_mat_elems - 1] = [0, 8] // Given the flat index, we compute the corresponding partition indices uint64_t m, n; unravel_index(part_pair_flat_idx, n_parts, &m, &n); @@ -489,7 +468,6 @@ __global__ void ari_kernel(T *parts, // int *s_contingency = shared_mem + 2 * n_objs; get_contingency_matrix_shared(t_data_part0, t_data_part1, n_objs, s_contingency, k); - /* * Step 3: Construct pair confusion matrix */ @@ -505,8 +483,7 @@ __global__ void ari_kernel(T *parts, long long fn = s_pair_confusion_matrix[2]; long long tp = s_pair_confusion_matrix[3]; float ari = 0.0f; - - + if (fn == 0 && fp == 0) { ari = 1.0f; @@ -514,9 +491,9 @@ __global__ void ari_kernel(T *parts, else { double numerator = 2.0 * ((double)tp * (double)tn - (double)fn * (double)fp); - double denominator = ((double)tp + (double)fn) * ((double)fn + (double)tn) + ((double)tp + (double)fp) * ((double)fp + (double)tn); + double denominator = ((double)tp + (double)fn) * ((double)fn + (double)tn) + + ((double)tp + (double)fp) * ((double)fp + (double)tn); ari = (float)(numerator / denominator); - } out[blockIdx.x] = ari; } @@ -539,31 +516,28 @@ __global__ void ari_kernel(T *parts, * @throws py::value_error on dtype, dimensionality, or shape mismatch */ template -T *process_input_array(const py::array_t &parts, - int64_t n_features = -1, int64_t n_parts = -1, int64_t n_objs = -1) +T *process_input_array(const py::array_t &parts, int64_t n_features = -1, int64_t n_parts = -1, + int64_t n_objs = -1) { py::buffer_info buffer = parts.request(); if (buffer.format != py::format_descriptor::format()) { - throw py::value_error( - std::string("Partitions array has an incompatible dtype: expected numpy format '") + - py::format_descriptor::format() + "' (e.g. int16), got '" + buffer.format + "'"); + throw py::value_error(std::string("Partitions array has an incompatible dtype: expected numpy format '") + + py::format_descriptor::format() + "' (e.g. int16), got '" + buffer.format + "'"); } if (buffer.ndim != 3) { - throw py::value_error( - "Partitions array must be 3-dimensional (n_features, n_parts, n_objs); got ndim=" + - std::to_string(buffer.ndim)); + throw py::value_error("Partitions array must be 3-dimensional (n_features, n_parts, n_objs); got ndim=" + + std::to_string(buffer.ndim)); } if (n_features >= 0 && n_parts >= 0 && n_objs >= 0) { if (buffer.shape[0] != n_features || buffer.shape[1] != n_parts || buffer.shape[2] != n_objs) { - throw py::value_error( - "Partitions array shape mismatch: expected (" + std::to_string(n_features) + ", " + - std::to_string(n_parts) + ", " + std::to_string(n_objs) + ") but got (" + - std::to_string(buffer.shape[0]) + ", " + std::to_string(buffer.shape[1]) + ", " + - std::to_string(buffer.shape[2]) + ")"); + throw py::value_error("Partitions array shape mismatch: expected (" + std::to_string(n_features) + ", " + + std::to_string(n_parts) + ", " + std::to_string(n_objs) + ") but got (" + + std::to_string(buffer.shape[0]) + ", " + std::to_string(buffer.shape[1]) + ", " + + std::to_string(buffer.shape[2]) + ")"); } } return static_cast(buffer.ptr); @@ -576,12 +550,8 @@ T *process_input_array(const py::array_t &parts, * @return std::unique_ptr to thrust device vector containing ARI values with type R */ template -auto ari_core_device(const T *parts, - const uint64_t n_features, - const uint64_t n_parts, - const uint64_t n_objs, - const uint64_t batch_start, - const uint64_t batch_size) -> std::unique_ptr> +auto ari_core_device(const T *parts, const uint64_t n_features, const uint64_t n_parts, const uint64_t n_objs, + const uint64_t batch_start, const uint64_t batch_size) -> std::unique_ptr> { /* * Show debugging and device information @@ -589,7 +559,10 @@ auto ari_core_device(const T *parts, // spdlog::debug("Max shared memory per block: {} bytes", get_max_shared_memory_per_block()); // Input validation - if (!parts || n_features == 0 || n_parts == 0 || n_objs == 0) { throw std::invalid_argument("Invalid input parameters"); } + if (!parts || n_features == 0 || n_parts == 0 || n_objs == 0) + { + throw std::invalid_argument("Invalid input parameters"); + } /* * Pre-computation @@ -599,7 +572,10 @@ auto ari_core_device(const T *parts, // Guard the unsigned subtraction below: check the bound BEFORE computing // (n_aris - batch_start), which would underflow to a huge value otherwise. - if (batch_start >= n_aris) { throw std::invalid_argument("Batch start index exceeds total number of ARIs"); } + if (batch_start >= n_aris) + { + throw std::invalid_argument("Batch start index exceeds total number of ARIs"); + } // Determine the actual batch size const auto actual_batch_size = batch_size == 0 ? n_aris : std::min(batch_size, n_aris - batch_start); @@ -629,8 +605,8 @@ auto ari_core_device(const T *parts, // Compute shared memory size // Align to 8-byte boundary for long long and add space for pair confusion matrix int offset_ints = k * k + 2 * k; - int aligned_offset_ints = ((offset_ints + 1) / 2) * 2; // Align to 8-byte (2 int) boundary - int s_mem_size = aligned_offset_ints * sz_int + 4 * sizeof(long long); // Total size + int aligned_offset_ints = ((offset_ints + 1) / 2) * 2; // Align to 8-byte (2 int) boundary + int s_mem_size = aligned_offset_ints * sz_int + 4 * sizeof(long long); // Total size // Check if shared memory size exceeds device limits auto [is_valid, message] = check_shared_memory_size(s_mem_size); @@ -647,46 +623,34 @@ auto ari_core_device(const T *parts, spdlog::debug("Memory before kernel launch: "); before_device_mem = print_cuda_memory_info(); - if (is_valid) { + if (is_valid) + { // Use shared memory kernel for smaller contingency matrices spdlog::debug("Using shared memory kernel for k={}", k); ari_kernel<<>>( - thrust::raw_pointer_cast(d_parts->data()), - actual_batch_size, - n_features, - n_parts, - n_objs, - n_parts * n_objs, - n_parts * n_parts, - k, - batch_start, - thrust::raw_pointer_cast(d_out->data())); + thrust::raw_pointer_cast(d_parts->data()), actual_batch_size, n_features, n_parts, n_objs, n_parts * n_objs, + n_parts * n_parts, k, batch_start, thrust::raw_pointer_cast(d_out->data())); CUDA_CHECK_KERNEL("ari_kernel"); - } else { + } + else + { // Use global memory kernel for larger contingency matrices spdlog::debug("Using global memory kernel for k={}", k); - + // Allocate global memory for contingency matrices // Each block needs: k*k (contingency) + 2*k (sum arrays) + 4 (pair confusion as long long) // Ensure 8-byte alignment for long long by using char vector - const size_t ints_per_block = k * k + 2 * k; // contingency + sum arrays - const size_t aligned_ints_per_block = ((ints_per_block + 1) / 2) * 2; // Align to 8-byte boundary + const size_t ints_per_block = k * k + 2 * k; // contingency + sum arrays + const size_t aligned_ints_per_block = ((ints_per_block + 1) / 2) * 2; // Align to 8-byte boundary const size_t mem_per_block = aligned_ints_per_block * sizeof(int) + 4 * sizeof(long long); const size_t total_global_mem = actual_batch_size * mem_per_block; - + auto d_global_cont_matrices = std::make_unique>(total_global_mem); - + ari_kernel_global<<>>( - thrust::raw_pointer_cast(d_parts->data()), - actual_batch_size, - n_features, - n_parts, - n_objs, - n_parts * n_objs, - n_parts * n_parts, - k, - batch_start, - reinterpret_cast(thrust::raw_pointer_cast(d_global_cont_matrices->data())), + thrust::raw_pointer_cast(d_parts->data()), actual_batch_size, n_features, n_parts, n_objs, n_parts * n_objs, + n_parts * n_parts, k, batch_start, + reinterpret_cast(thrust::raw_pointer_cast(d_global_cont_matrices->data())), thrust::raw_pointer_cast(d_out->data())); CUDA_CHECK_KERNEL("ari_kernel_global"); } @@ -710,16 +674,11 @@ auto ari_core_device(const T *parts, * @return std::unique_ptr to thrust device vector containing ARI values */ template -auto ari_core_device(const py::array_t &parts, - const size_t n_features, - const size_t n_parts, - const size_t n_objs, - const uint64_t batch_start, - const uint64_t batch_size) -> std::unique_ptr> +auto ari_core_device(const py::array_t &parts, const size_t n_features, const size_t n_parts, + const size_t n_objs, const uint64_t batch_start, const uint64_t batch_size) + -> std::unique_ptr> { - const auto parts_ptr = process_input_array(parts, - static_cast(n_features), - static_cast(n_parts), + const auto parts_ptr = process_input_array(parts, static_cast(n_features), static_cast(n_parts), static_cast(n_objs)); return ari_core_device(parts_ptr, n_features, n_parts, n_objs, batch_start, batch_size); } @@ -731,12 +690,8 @@ auto ari_core_device(const py::array_t &parts, * @return std::vector ARI values for each pair of partitions stored in host memory */ template -auto ari_core_host(const T *parts, - const size_t n_features, - const size_t n_parts, - const size_t n_objs, - const uint64_t batch_start, - const uint64_t batch_size) -> std::vector +auto ari_core_host(const T *parts, const size_t n_features, const size_t n_parts, const size_t n_objs, + const uint64_t batch_start, const uint64_t batch_size) -> std::vector { /* * Pre-computation @@ -747,7 +702,10 @@ auto ari_core_host(const T *parts, // Guard the unsigned subtraction below (see ari_core_device): validate the // bound before computing (n_aris - batch_start). - if (batch_start >= n_aris) { throw std::invalid_argument("Batch start index exceeds total number of ARIs"); } + if (batch_start >= n_aris) + { + throw std::invalid_argument("Batch start index exceeds total number of ARIs"); + } // Determine the actual batch size const auto actual_batch_size = batch_size == 0 ? n_aris : std::min(batch_size, n_aris - batch_start); @@ -784,16 +742,10 @@ auto ari_core_host(const T *parts, * @return std::vector All ARI values for each pair of partitions */ template -auto ari(const py::array_t &parts, - const size_t n_features, - const size_t n_parts, - const size_t n_objs, - const uint64_t batch_start, - const uint64_t batch_size) -> std::vector +auto ari(const py::array_t &parts, const size_t n_features, const size_t n_parts, + const size_t n_objs, const uint64_t batch_start, const uint64_t batch_size) -> std::vector { - const auto parts_ptr = process_input_array(parts, - static_cast(n_features), - static_cast(n_parts), + const auto parts_ptr = process_input_array(parts, static_cast(n_features), static_cast(n_parts), static_cast(n_objs)); return ari_core_host(parts_ptr, n_features, n_parts, n_objs, batch_start, batch_size); } @@ -806,28 +758,16 @@ auto ari(const py::array_t &parts, // by the linker. // Used for external python testing -template auto ari( - const py::array_t &parts, - const size_t n_features, - const size_t n_parts, - const size_t n_objs, - const uint64_t batch_start, - const uint64_t batch_size) -> std::vector; +template auto ari(const py::array_t &parts, const size_t n_features, const size_t n_parts, + const size_t n_objs, const uint64_t batch_start, const uint64_t batch_size) + -> std::vector; // Used for internal c++ testing -template auto ari_core_host( - const int *parts, - const size_t n_features, - const size_t n_parts, - const size_t n_objs, - const uint64_t batch_start, - const uint64_t batch_size) -> std::vector; +template auto ari_core_host(const int *parts, const size_t n_features, const size_t n_parts, const size_t n_objs, + const uint64_t batch_start, const uint64_t batch_size) -> std::vector; // Used in the coef API -template auto ari_core_device( - const py::array_t &parts, - const uint64_t n_features, - const uint64_t n_parts, - const uint64_t n_objs, - const uint64_t batch_start, - const uint64_t batch_size) -> std::unique_ptr>; +template auto ari_core_device(const py::array_t &parts, + const uint64_t n_features, const uint64_t n_parts, const uint64_t n_objs, + const uint64_t batch_start, const uint64_t batch_size) + -> std::unique_ptr>; diff --git a/libs/ccc_cuda_ext/metrics.cuh b/libs/ccc_cuda_ext/metrics.cuh index bbb9a664..131b21d9 100644 --- a/libs/ccc_cuda_ext/metrics.cuh +++ b/libs/ccc_cuda_ext/metrics.cuh @@ -1,10 +1,10 @@ #pragma once -#include -#include #include +#include #include #include +#include namespace py = pybind11; @@ -35,8 +35,7 @@ enum class PartPairValidity : int * fills an entire partition row with the same marker, inspecting index 0 is * sufficient and permutation of object order does not change the class. */ -template -__device__ __host__ inline PartPairValidity classify_partition_pair(T a0, T b0) +template __device__ __host__ inline PartPairValidity classify_partition_pair(T a0, T b0) { if (a0 == static_cast(-1) || b0 == static_cast(-1)) { @@ -51,27 +50,16 @@ __device__ __host__ inline PartPairValidity classify_partition_pair(T a0, T b0) // Used for external python testing template -auto ari(const py::array_t &parts, - const size_t n_features, - const size_t n_parts, - const size_t n_objs, - const uint64_t batch_start = 0, - const uint64_t batch_size = 0) -> std::vector; +auto ari(const py::array_t &parts, const size_t n_features, const size_t n_parts, + const size_t n_objs, const uint64_t batch_start = 0, const uint64_t batch_size = 0) -> std::vector; // Used for internal c++ testing template -auto ari_core_host(const T *parts, - const size_t n_features, - const size_t n_parts, - const size_t n_objs, - const uint64_t batch_start = 0, - const uint64_t batch_size = 0) -> std::vector; +auto ari_core_host(const T *parts, const size_t n_features, const size_t n_parts, const size_t n_objs, + const uint64_t batch_start = 0, const uint64_t batch_size = 0) -> std::vector; // Used in the coef API template -auto ari_core_device(const py::array_t &parts, - const uint64_t n_features, - const uint64_t n_parts, - const uint64_t n_objs, - const uint64_t batch_start = 0, - const uint64_t batch_size = 0) -> std::unique_ptr>; +auto ari_core_device(const py::array_t &parts, const uint64_t n_features, const uint64_t n_parts, + const uint64_t n_objs, const uint64_t batch_start = 0, const uint64_t batch_size = 0) + -> std::unique_ptr>; diff --git a/libs/ccc_cuda_ext/utils.cuh b/libs/ccc_cuda_ext/utils.cuh index 377ebeef..5cf3ecc5 100644 --- a/libs/ccc_cuda_ext/utils.cuh +++ b/libs/ccc_cuda_ext/utils.cuh @@ -18,15 +18,15 @@ #pragma once +#include #include -#include -#include -#include -#include -#include #include -#include +#include #include +#include +#include +#include +#include /** * @brief Internal function for CUDA error checking @@ -46,8 +46,8 @@ inline void gpu_assert(cudaError_t code, const char *file, int line, bool abort { if (code != cudaSuccess) { - const std::string msg = std::string("CUDA Error: ") + cudaGetErrorString(code) + - " at " + file + ":" + std::to_string(line); + const std::string msg = + std::string("CUDA Error: ") + cudaGetErrorString(code) + " at " + file + ":" + std::to_string(line); spdlog::error(msg); if (abort) { @@ -64,10 +64,10 @@ inline void gpu_assert(cudaError_t code, const char *file, int line, bool abort * * @param ans The CUDA operation to check (e.g., cudaMalloc, cudaMemcpy) */ -#define CUDA_CHECK(ans) \ - do \ - { \ - gpu_assert((ans), __FILE__, __LINE__, true); \ +#define CUDA_CHECK(ans) \ + do \ + { \ + gpu_assert((ans), __FILE__, __LINE__, true); \ } while (0) /** @@ -80,10 +80,10 @@ inline void gpu_assert(cudaError_t code, const char *file, int line, bool abort * @brief Optional error checking macro that only logs errors (never throws). * Use for cleanup/best-effort operations where failure can be tolerated. */ -#define CUDA_CHECK_OPTIONAL(ans) \ - do \ - { \ - gpu_assert((ans), __FILE__, __LINE__, false); \ +#define CUDA_CHECK_OPTIONAL(ans) \ + do \ + { \ + gpu_assert((ans), __FILE__, __LINE__, false); \ } while (0) /** @@ -97,29 +97,27 @@ inline void gpu_assert(cudaError_t code, const char *file, int line, bool abort * * @param kernel_name A human-readable name for the launched kernel. */ -#define CUDA_CHECK_KERNEL(kernel_name) \ - do \ - { \ - cudaError_t _launch_err = cudaGetLastError(); \ - if (_launch_err != cudaSuccess) \ - { \ - const std::string _msg = \ - std::string("CUDA kernel launch failed for '") + (kernel_name) + \ - "': " + cudaGetErrorString(_launch_err) + " at " + __FILE__ + \ - ":" + std::to_string(__LINE__); \ - spdlog::error(_msg); \ - throw std::runtime_error(_msg); \ - } \ - cudaError_t _sync_err = cudaDeviceSynchronize(); \ - if (_sync_err != cudaSuccess) \ - { \ - const std::string _msg = \ - std::string("CUDA kernel execution failed for '") + (kernel_name) + \ - "': " + cudaGetErrorString(_sync_err) + " at " + __FILE__ + \ - ":" + std::to_string(__LINE__); \ - spdlog::error(_msg); \ - throw std::runtime_error(_msg); \ - } \ +#define CUDA_CHECK_KERNEL(kernel_name) \ + do \ + { \ + cudaError_t _launch_err = cudaGetLastError(); \ + if (_launch_err != cudaSuccess) \ + { \ + const std::string _msg = std::string("CUDA kernel launch failed for '") + (kernel_name) + \ + "': " + cudaGetErrorString(_launch_err) + " at " + __FILE__ + ":" + \ + std::to_string(__LINE__); \ + spdlog::error(_msg); \ + throw std::runtime_error(_msg); \ + } \ + cudaError_t _sync_err = cudaDeviceSynchronize(); \ + if (_sync_err != cudaSuccess) \ + { \ + const std::string _msg = std::string("CUDA kernel execution failed for '") + (kernel_name) + \ + "': " + cudaGetErrorString(_sync_err) + " at " + __FILE__ + ":" + \ + std::to_string(__LINE__); \ + spdlog::error(_msg); \ + throw std::runtime_error(_msg); \ + } \ } while (0) /** @@ -183,21 +181,21 @@ inline std::tuple check_shared_memory_size(const size_t reque if (requested_size > max_shared_mem) { - return std::make_tuple(false, - std::string("Required shared memory (") + std::to_string(requested_size) + - " bytes) exceeds device limit (" + std::to_string(max_shared_mem) + " bytes)"); + return std::make_tuple(false, std::string("Required shared memory (") + std::to_string(requested_size) + + " bytes) exceeds device limit (" + std::to_string(max_shared_mem) + + " bytes)"); } // Optionally warn if close to limit (e.g., using more than 90%) if (requested_size > (max_shared_mem * 0.9)) { - return std::make_tuple(true, - std::string("Warning: Shared memory usage (") + std::to_string(requested_size) + - " bytes) is close to device limit (" + std::to_string(max_shared_mem) + " bytes)"); + return std::make_tuple(true, std::string("Warning: Shared memory usage (") + std::to_string(requested_size) + + " bytes) is close to device limit (" + std::to_string(max_shared_mem) + + " bytes)"); } return std::make_tuple(true, std::string("Required shared memory size: ") + std::to_string(requested_size) + - " bytes, within available: " + std::to_string(max_shared_mem) + " bytes"); + " bytes, within available: " + std::to_string(max_shared_mem) + " bytes"); } /** @@ -251,27 +249,19 @@ inline void print_cuda_device_info(int device_id = 0) cudaDriverGetVersion(&driverVersion); cudaRuntimeGetVersion(&runtimeVersion); - spdlog::debug(" CUDA Driver Version / Runtime Version {}.{} / {}.{}", - driverVersion / 1000, (driverVersion % 100) / 10, - runtimeVersion / 1000, (runtimeVersion % 100) / 10); - spdlog::debug(" CUDA Capability Major/Minor version number: {}.{}", - deviceProp.major, deviceProp.minor); + spdlog::debug(" CUDA Driver Version / Runtime Version {}.{} / {}.{}", driverVersion / 1000, + (driverVersion % 100) / 10, runtimeVersion / 1000, (runtimeVersion % 100) / 10); + spdlog::debug(" CUDA Capability Major/Minor version number: {}.{}", deviceProp.major, deviceProp.minor); spdlog::debug(" Total amount of global memory: {:.2f} GBytes ({} bytes)", - (float)deviceProp.totalGlobalMem / pow(1024.0, 3), - (unsigned long long)deviceProp.totalGlobalMem); + (float)deviceProp.totalGlobalMem / pow(1024.0, 3), (unsigned long long)deviceProp.totalGlobalMem); spdlog::debug(" GPU Clock rate: {:.0f} MHz ({:.2f} GHz)", - deviceProp.clockRate * 1e-3f, - deviceProp.clockRate * 1e-6f); - spdlog::debug(" Memory Clock rate: {:.0f} Mhz", - deviceProp.memoryClockRate * 1e-3f); - spdlog::debug(" Memory Bus Width: {}-bit", - deviceProp.memoryBusWidth); - spdlog::debug(" Shared Memory per Block: {:.2f} KB", - deviceProp.sharedMemPerBlock / 1024.0f); + deviceProp.clockRate * 1e-3f, deviceProp.clockRate * 1e-6f); + spdlog::debug(" Memory Clock rate: {:.0f} Mhz", deviceProp.memoryClockRate * 1e-3f); + spdlog::debug(" Memory Bus Width: {}-bit", deviceProp.memoryBusWidth); + spdlog::debug(" Shared Memory per Block: {:.2f} KB", deviceProp.sharedMemPerBlock / 1024.0f); spdlog::debug(" Shared Memory per Multiprocessor: {:.2f} KB", deviceProp.sharedMemPerMultiprocessor / 1024.0f); - spdlog::debug(" Number of Multiprocessors: {}", - deviceProp.multiProcessorCount); + spdlog::debug(" Number of Multiprocessors: {}", deviceProp.multiProcessorCount); } /** @@ -296,8 +286,7 @@ inline size_t print_cuda_memory_info(int device_id = 0) cudaMemGetInfo(&free_mem, &total_mem); // Helper function to format memory size with appropriate unit - auto format_memory = [](size_t bytes) -> std::string - { + auto format_memory = [](size_t bytes) -> std::string { const size_t KB = 1024; const size_t MB = KB * 1024; const size_t GB = MB * 1024; @@ -320,8 +309,7 @@ inline size_t print_cuda_memory_info(int device_id = 0) } }; - spdlog::debug("Free memory: {}, Total memory: {}", - format_memory(free_mem), format_memory(total_mem)); + spdlog::debug("Free memory: {}, Total memory: {}", format_memory(free_mem), format_memory(total_mem)); return free_mem; } diff --git a/pyproject.toml b/pyproject.toml index 58ab0de3..594bc2c6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -65,6 +65,96 @@ minversion = "6.0" addopts = "-ra -q" testpaths = ["tests"] +# --------------------------------------------------------------------------- # +# Ruff — single source of truth for Python lint + format. +# +# Scope: the installable package (`libs/ccc`) and the test suite (`tests`). +# `analysis/` and `docs/` are manuscript/notebook/tooling code intentionally +# left untouched and are excluded below. `force-exclude` makes the excludes +# apply even when pre-commit passes explicit file paths. +# --------------------------------------------------------------------------- # +[tool.ruff] +target-version = "py310" +line-length = 88 +force-exclude = true +extend-exclude = [ + "analysis", + "docs", + "build", + "_skbuild", + "*.ipynb", +] + +[tool.ruff.lint] +select = [ + "E", # pycodestyle errors + "F", # pyflakes + "W", # pycodestyle warnings + "B", # flake8-bugbear + "I", # isort + "ARG", # flake8-unused-arguments + "C4", # flake8-comprehensions + "EM", # flake8-errmsg + "PL", # pylint + "PT", # flake8-pytest-style + "PTH", # flake8-use-pathlib + "RET", # flake8-return + "SIM", # flake8-simplify + "UP", # pyupgrade + "NPY", # numpy-specific rules + "PD", # pandas-vet +] +# Intentionally NOT enabling isort's `required-imports` (no forced +# `from __future__ import annotations`) to avoid a semantic-adjacent mass edit. +ignore = [ + # --- Owned by the formatter / structural complexity (permanent) --- + "E501", # line length is handled by the formatter + "PLR0911", # too many return statements + "PLR0912", # too many branches + "PLR0913", # too many arguments — numeric kernels legitimately take many + "PLR0915", # too many statements + "PLR2004", # magic value comparison — pervasive in numeric code + + # --- Deferred to follow-up work (this change lands the tooling + + # mechanical format only; the fixes below are logic/behavior-adjacent + # and need manual review, so they are suppressed for now). --- + "NPY002", # legacy np.random.* — rewriting to Generator changes RNG streams/seeds + "PLC0415", # import-outside-top-level — many are intentional lazy imports (cupy/torch/heavy deps) + "EM101", # raw-string-in-exception — mechanical message extraction + "EM102", # f-string-in-exception — mechanical message extraction + "PT006", # pytest parametrize names as tuple vs string + "PT018", # composite assertion in tests + "PTH103", # os.makedirs -> Path.mkdir + "PTH118", # os.path.join -> Path + "PTH123", # open() -> Path.open() + "B905", # zip() without strict= — adding strict= is behavior-changing + "PLW1510", # subprocess.run without check= — adding check= is behavior-changing + "PLW2901", # redefined loop variable + "PLR1704", # redefined argument from local + "RET504", # unnecessary assignment before return + "C403", # unnecessary list comprehension (set) + "C404", # unnecessary list comprehension (dict) + "C408", # unnecessary dict/list/tuple call + "C409", # unnecessary literal within tuple call + "C417", # unnecessary map() usage + "SIM105", # suppressible exception -> contextlib.suppress + "SIM113", # use enumerate() + "SIM117", # combine nested with-statements + "ARG002", # unused method argument + "UP028", # yield-in-for-loop -> yield from +] + +[tool.ruff.lint.per-file-ignores] +"tests/**" = [ + "PLR2004", # magic values are expected in test assertions + "ARG", # pytest fixtures are passed but not always referenced + "B", # bugbear patterns are acceptable in tests + "PT011", # broad pytest.raises without match is acceptable in tests +] + +[tool.ruff.format] +docstring-code-format = true + [tool.cibuildwheel] build = "cp3{10,11,12,13,14}-manylinux_x86_64" manylinux-x86_64-image = "manylinux_2_28" diff --git a/tests/gpu/test_ari_gpu.py b/tests/gpu/test_ari_gpu.py index f607d193..ad24ddd3 100644 --- a/tests/gpu/test_ari_gpu.py +++ b/tests/gpu/test_ari_gpu.py @@ -1,8 +1,8 @@ import time -import pytest -import numpy as np -import ccc_cuda_ext +import ccc_cuda_ext +import numpy as np +import pytest from ccc.sklearn.metrics import ( adjusted_rand_index, ) diff --git a/tests/gpu/test_ccc_gpu.py b/tests/gpu/test_ccc_gpu.py index efb32513..20c78fb8 100644 --- a/tests/gpu/test_ccc_gpu.py +++ b/tests/gpu/test_ccc_gpu.py @@ -1,20 +1,21 @@ +import os import time -import pytest +from typing import Any + import numpy as np -from typing import Tuple, Optional, Dict, Any -import os import pandas as pd -from ccc.coef.impl_gpu import ccc as ccc_gpu +import pytest from ccc.coef.impl import ccc +from ccc.coef.impl_gpu import ccc as ccc_gpu from utils import clean_gpu_memory, generate_categorical_data def setup_logging( seed: int, - shape: Tuple[int, int], + shape: tuple[int, int], n_cpu_cores: int, generate_logs: bool, -) -> Optional[Dict[str, Any]]: +) -> dict[str, Any] | None: """Setup logging infrastructure if logging is enabled. Args: @@ -47,7 +48,7 @@ def setup_logging( return {"files": log_files, "log_file": open(log_files["log"], "w")} -def log_test_info(log_file, shape: Tuple[int, int], seed: int) -> None: +def log_test_info(log_file, shape: tuple[int, int], seed: int) -> None: """Log basic test information.""" print( f"\nTesting with {shape[0]} features, {shape[1]} samples, seed {seed}", @@ -69,9 +70,9 @@ def log_performance_metrics( def analyze_differences( c1: np.ndarray, c2: np.ndarray, - shape: Tuple[int, int], - log_file: Optional[Any] = None, -) -> Tuple[int, float, float, int]: + shape: tuple[int, int], + log_file: Any | None = None, +) -> tuple[int, float, float, int]: """Analyze differences between GPU and CPU results. Returns: @@ -101,7 +102,7 @@ def log_differences( c1: np.ndarray, c2: np.ndarray, not_close_indices: np.ndarray, - shape: Tuple[int, int], + shape: tuple[int, int], log_file: Any, ) -> None: """Log detailed information about differences between results.""" @@ -193,7 +194,7 @@ def log_statistics( @clean_gpu_memory def test_ccc_gpu_with_numerical_input( seed: int, - shape: Tuple[int, int], + shape: tuple[int, int], contain_singletons: bool, n_cpu_cores: int, max_not_close_percentage: float, @@ -277,14 +278,14 @@ def test_ccc_gpu_with_numerical_input( ((10, 20), 10, 2), ((20, 200), 50, 3), ((30, 300), 200, 4), - ((9, 10000), 500, 5) + ((9, 10000), 500, 5), ], ) @pytest.mark.parametrize("n_cpu_cores", [48]) @clean_gpu_memory def test_ccc_gpu_with_categorical_input( seed: int, - shape: Tuple[int, int], + shape: tuple[int, int], n_categories: int, str_length: int, n_cpu_cores: int, @@ -300,6 +301,7 @@ def test_ccc_gpu_with_categorical_input( gpu_df = pd.DataFrame(res_gpu.astype(np.float64)) pd.testing.assert_frame_equal(gpu_df, cpu_df, atol=1e-6, rtol=1e-6) + # @clean_gpu_memory # def test_ccc_gpu_with_mixed_input(): -# return \ No newline at end of file +# return diff --git a/tests/gpu/test_ccc_gpu_pvalue.py b/tests/gpu/test_ccc_gpu_pvalue.py index e3da0e79..0d7bc145 100644 --- a/tests/gpu/test_ccc_gpu_pvalue.py +++ b/tests/gpu/test_ccc_gpu_pvalue.py @@ -1,13 +1,8 @@ -import time -import pytest import numpy as np -from typing import Tuple, Optional, Dict, Any -import os -import pandas as pd -from sklearn.preprocessing import minmax_scale -from ccc.coef.impl_gpu import ccc as ccc_gpu +import pytest from ccc.coef.impl import ccc -from utils import clean_gpu_memory, generate_categorical_data +from ccc.coef.impl_gpu import ccc as ccc_gpu +from sklearn.preprocessing import minmax_scale # Original CCC test from tests/test_coef_pval.py @@ -220,60 +215,58 @@ def test_cm_single_argument_is_matrix(): assert pvalue[2] > 0.10 -@pytest.mark.parametrize( - "seed", [42, 123] -) +@pytest.mark.parametrize("seed", [42, 123]) @pytest.mark.parametrize( "shape, pvalue_n_perms", [ # Small cases with different permutation counts ((3, 50), 50), - ((5, 100), 100), + ((5, 100), 100), ((4, 80), 200), # Medium cases ((10, 200), 100), ((15, 150), 50), # ((100, 1000), 50), # Edge cases - ((2, 30), 20), # Minimum features + ((2, 30), 20), # Minimum features ], ) def test_ccc_gpu_vs_cpu_pvalue_comparison( seed: int, - shape: Tuple[int, int], + shape: tuple[int, int], pvalue_n_perms: int, ): """ Parameterized test to compare p-values generated by GPU and CPU implementations. - + This test verifies that: 1. GPU and CPU p-values are statistically similar 2. Both implementations handle various data shapes correctly 3. P-values are in valid range [0, 1] 4. The correlation between GPU and CPU p-values is high - + Args: seed: Random seed for reproducibility shape: Tuple of (n_features, n_samples) pvalue_n_perms: Number of permutations for p-value computation """ n_features, n_samples = shape - + # Generate reproducible test data np.random.seed(seed) - + # Create data with some linear relationships and some random data = np.random.rand(n_features, n_samples) - + # Add some linear relationships to make results more interesting if n_features >= 3: # Make feature 1 linearly related to feature 0 with noise data[1, :] = data[0, :] * 2.0 + np.random.normal(0, 0.1, n_samples) - + # Make feature 2 have a quadratic relationship with feature 0 plus noise if n_features >= 4: data[2, :] = data[0, :] ** 2 + np.random.normal(0, 0.1, n_samples) - + # Run GPU implementation try: gpu_result = ccc_gpu(data, pvalue_n_perms=pvalue_n_perms) @@ -281,15 +274,15 @@ def test_ccc_gpu_vs_cpu_pvalue_comparison( gpu_ccc, gpu_pvalues = gpu_result except Exception as e: pytest.fail(f"GPU implementation failed: {e}") - - # Run CPU implementation + + # Run CPU implementation try: cpu_result = ccc(data, pvalue_n_perms=pvalue_n_perms) assert len(cpu_result) == 2, "CPU should return (ccc_values, pvalues)" cpu_ccc, cpu_pvalues = cpu_result except Exception as e: pytest.fail(f"CPU implementation failed: {e}") - + # Convert scalar results to arrays for consistent handling if np.isscalar(gpu_ccc): gpu_ccc = np.array([gpu_ccc]) @@ -299,92 +292,115 @@ def test_ccc_gpu_vs_cpu_pvalue_comparison( gpu_pvalues = np.array([gpu_pvalues]) if np.isscalar(cpu_pvalues): cpu_pvalues = np.array([cpu_pvalues]) - + # Basic shape and type validations - assert gpu_ccc.shape == cpu_ccc.shape, f"CCC shapes don't match: GPU {gpu_ccc.shape} vs CPU {cpu_ccc.shape}" - assert gpu_pvalues.shape == cpu_pvalues.shape, f"P-value shapes don't match: GPU {gpu_pvalues.shape} vs CPU {cpu_pvalues.shape}" - + assert gpu_ccc.shape == cpu_ccc.shape, ( + f"CCC shapes don't match: GPU {gpu_ccc.shape} vs CPU {cpu_ccc.shape}" + ) + assert gpu_pvalues.shape == cpu_pvalues.shape, ( + f"P-value shapes don't match: GPU {gpu_pvalues.shape} vs CPU {cpu_pvalues.shape}" + ) + # Expected number of comparisons for symmetric matrix expected_comparisons = n_features * (n_features - 1) // 2 - assert len(gpu_ccc) == expected_comparisons, f"Expected {expected_comparisons} comparisons, got {len(gpu_ccc)}" - assert len(gpu_pvalues) == expected_comparisons, f"Expected {expected_comparisons} p-values, got {len(gpu_pvalues)}" - + assert len(gpu_ccc) == expected_comparisons, ( + f"Expected {expected_comparisons} comparisons, got {len(gpu_ccc)}" + ) + assert len(gpu_pvalues) == expected_comparisons, ( + f"Expected {expected_comparisons} p-values, got {len(gpu_pvalues)}" + ) + # Validate p-values are in [0, 1] range - assert np.all((gpu_pvalues >= 0) & (gpu_pvalues <= 1)), "GPU p-values should be in [0, 1]" - assert np.all((cpu_pvalues >= 0) & (cpu_pvalues <= 1)), "CPU p-values should be in [0, 1]" - + assert np.all((gpu_pvalues >= 0) & (gpu_pvalues <= 1)), ( + "GPU p-values should be in [0, 1]" + ) + assert np.all((cpu_pvalues >= 0) & (cpu_pvalues <= 1)), ( + "CPU p-values should be in [0, 1]" + ) + # Check CCC values are close (they should be nearly identical) np.testing.assert_allclose( - gpu_ccc, cpu_ccc, - rtol=1e-5, atol=1e-6, - err_msg="GPU and CPU CCC values should be nearly identical" + gpu_ccc, + cpu_ccc, + rtol=1e-5, + atol=1e-6, + err_msg="GPU and CPU CCC values should be nearly identical", ) - + # For p-values, we expect them to be close but not identical due to: # 1. Different random number generation # 2. Potential floating point differences # 3. Different permutation ordering - + # Check if any values are NaN (should be consistent) gpu_nan_mask = np.isnan(gpu_pvalues) cpu_nan_mask = np.isnan(cpu_pvalues) np.testing.assert_array_equal( - gpu_nan_mask, cpu_nan_mask, - err_msg="NaN patterns should be identical between GPU and CPU" + gpu_nan_mask, + cpu_nan_mask, + err_msg="NaN patterns should be identical between GPU and CPU", ) - + # For non-NaN values, check statistical properties valid_mask = ~(gpu_nan_mask | cpu_nan_mask) if np.any(valid_mask): gpu_valid = gpu_pvalues[valid_mask] cpu_valid = cpu_pvalues[valid_mask] - + # Check that most p-values are reasonably close # Use a more relaxed tolerance since p-values can vary due to randomness abs_diff = np.abs(gpu_valid - cpu_valid) max_diff = np.max(abs_diff) mean_diff = np.mean(abs_diff) - + # For small number of permutations, p-values can vary significantly # Use adaptive thresholds based on permutation count close_threshold = max(0.15, 1.0 / pvalue_n_perms) # At least 1/n_perms or 0.15 correlation_threshold = 0.5 if pvalue_n_perms < 50 else 0.7 - + close_mask = abs_diff < close_threshold close_percentage = np.mean(close_mask) - + # P-values should be reasonably correlated (only for multiple values) correlation = np.nan if len(gpu_valid) > 1: correlation = np.corrcoef(gpu_valid, cpu_valid)[0, 1] - assert correlation > correlation_threshold, f"P-value correlation {correlation:.3f} is too low (should be > {correlation_threshold})" - - print(f"P-value comparison stats:") + assert correlation > correlation_threshold, ( + f"P-value correlation {correlation:.3f} is too low (should be > {correlation_threshold})" + ) + + print("P-value comparison stats:") print(f" Shape: {shape}, Seed: {seed}, Permutations: {pvalue_n_perms}") print(f" Max difference: {max_diff:.4f}") print(f" Mean difference: {mean_diff:.4f}") if not np.isnan(correlation): print(f" Correlation: {correlation:.4f}") else: - print(f" Correlation: N/A (single value)") + print(" Correlation: N/A (single value)") print(f" Close (diff < {close_threshold:.3f}): {close_percentage:.1%}") - + # Adaptive thresholds: lower requirements for small permutation counts and small datasets # Note: For very small datasets with few permutations, p-values can vary significantly # due to the discrete nature of permutation tests - is_edge_case = (n_features == 2 and pvalue_n_perms <= 20) - + is_edge_case = n_features == 2 and pvalue_n_perms <= 20 + if is_edge_case: # Very lenient thresholds for edge cases - min_close_percentage = 0.0 # Allow any difference for single comparison with few permutations - max_diff_threshold = 1.0 # Allow full range of p-value differences + min_close_percentage = ( + 0.0 # Allow any difference for single comparison with few permutations + ) + max_diff_threshold = 1.0 # Allow full range of p-value differences else: min_close_percentage = 0.5 if pvalue_n_perms < 100 else 0.7 max_diff_threshold = 0.8 if pvalue_n_perms < 50 else 0.5 - + # At least min_close_percentage should be within threshold if min_close_percentage > 0: - assert close_percentage >= min_close_percentage, f"Only {close_percentage:.1%} of p-values are close (should be ≥{min_close_percentage:.0%})" - - # Maximum difference shouldn't be too extreme - assert max_diff < max_diff_threshold, f"Maximum p-value difference {max_diff:.4f} is too large (should be < {max_diff_threshold})" + assert close_percentage >= min_close_percentage, ( + f"Only {close_percentage:.1%} of p-values are close (should be ≥{min_close_percentage:.0%})" + ) + + # Maximum difference shouldn't be too extreme + assert max_diff < max_diff_threshold, ( + f"Maximum p-value difference {max_diff:.4f} is too large (should be < {max_diff_threshold})" + ) diff --git a/tests/gpu/test_ccc_gpu_return_parts.py b/tests/gpu/test_ccc_gpu_return_parts.py index 7c10d58c..69f0ec26 100644 --- a/tests/gpu/test_ccc_gpu_return_parts.py +++ b/tests/gpu/test_ccc_gpu_return_parts.py @@ -1,20 +1,21 @@ +import os import time -import pytest +from typing import Any + import numpy as np -from typing import Tuple, Optional, Dict, Any -import os import pandas as pd -from ccc.coef.impl_gpu import ccc as ccc_gpu +import pytest from ccc.coef.impl import ccc -from utils import clean_gpu_memory, generate_categorical_data +from ccc.coef.impl_gpu import ccc as ccc_gpu +from utils import clean_gpu_memory def setup_logging( seed: int, - shape: Tuple[int, int], + shape: tuple[int, int], n_cpu_cores: int, generate_logs: bool, -) -> Optional[Dict[str, Any]]: +) -> dict[str, Any] | None: """Setup logging infrastructure if logging is enabled. Args: @@ -47,7 +48,7 @@ def setup_logging( return {"files": log_files, "log_file": open(log_files["log"], "w")} -def log_test_info(log_file, shape: Tuple[int, int], seed: int) -> None: +def log_test_info(log_file, shape: tuple[int, int], seed: int) -> None: """Log basic test information.""" print( f"\nTesting with {shape[0]} features, {shape[1]} samples, seed {seed}", @@ -69,9 +70,9 @@ def log_performance_metrics( def analyze_differences( c1: np.ndarray, c2: np.ndarray, - shape: Tuple[int, int], - log_file: Optional[Any] = None, -) -> Tuple[int, float, float, int]: + shape: tuple[int, int], + log_file: Any | None = None, +) -> tuple[int, float, float, int]: """Analyze differences between GPU and CPU results. Returns: @@ -101,7 +102,7 @@ def log_differences( c1: np.ndarray, c2: np.ndarray, not_close_indices: np.ndarray, - shape: Tuple[int, int], + shape: tuple[int, int], log_file: Any, ) -> None: """Log detailed information about differences between results.""" @@ -182,7 +183,7 @@ def test_cm_return_parts_quadratic(): # - k == 3 for feature0 # - k == 2 for feature1 np.testing.assert_array_equal(max_parts, np.array([1, 0])) - + def test_cm_return_parts_linear(): # Prepare @@ -236,7 +237,7 @@ def test_cm_return_parts_linear(): @clean_gpu_memory def test_ccc_gpu_with_numerical_input( seed: int, - shape: Tuple[int, int], + shape: tuple[int, int], contain_singletons: bool, n_cpu_cores: int, generate_logs: bool, @@ -304,48 +305,57 @@ def test_ccc_gpu_with_numerical_input( # Check if return_parts is correctly implemented assert g_max_parts.shape == c_max_parts.shape assert g_parts.shape == c_parts.shape - + for i in range(len(g_parts)): - pd.testing.assert_frame_equal(pd.DataFrame(g_parts[i].astype(np.int16)), pd.DataFrame(c_parts[i]), check_exact=True) - + pd.testing.assert_frame_equal( + pd.DataFrame(g_parts[i].astype(np.int16)), + pd.DataFrame(c_parts[i]), + check_exact=True, + ) + # Validate max_parts: Both GPU and CPU choices should represent valid maxima # This handles tie-breaking differences between implementations for i in range(len(g_max_parts)): gpu_parts = g_max_parts[i] cpu_parts = c_max_parts[i] - + # If they match exactly, no need for further validation if np.array_equal(gpu_parts, cpu_parts): continue - + # For mismatches, verify both choices are valid maxima # We need to compute the ARI matrix to check this from ccc.coef.impl import cdist_parts_basic, get_coords_from_index - + # Get feature indices for this comparison feat_i, feat_j = get_coords_from_index(shape[0], i) - + # Compute ARI matrix for this feature pair ari_matrix = cdist_parts_basic(c_parts[feat_i], c_parts[feat_j]) max_ari = np.max(ari_matrix) - + # Check GPU choice gpu_ari = ari_matrix[gpu_parts[0], gpu_parts[1]] - # Check CPU choice + # Check CPU choice cpu_ari = ari_matrix[cpu_parts[0], cpu_parts[1]] - + # Both should be maximum values (within floating point tolerance) - assert np.abs(gpu_ari - max_ari) < 1e-8, f"GPU choice at comparison {i} is not maximum: {gpu_ari} vs {max_ari}" - assert np.abs(cpu_ari - max_ari) < 1e-8, f"CPU choice at comparison {i} is not maximum: {cpu_ari} vs {max_ari}" - + assert np.abs(gpu_ari - max_ari) < 1e-8, ( + f"GPU choice at comparison {i} is not maximum: {gpu_ari} vs {max_ari}" + ) + assert np.abs(cpu_ari - max_ari) < 1e-8, ( + f"CPU choice at comparison {i} is not maximum: {cpu_ari} vs {max_ari}" + ) + # Print info about the tie for debugging (optional) ties = np.where(np.abs(ari_matrix - max_ari) < 1e-8) num_ties = len(ties[0]) if num_ties > 1: - print(f"Tie detected at comparison {i} (features {feat_i} vs {feat_j}): {num_ties} positions with ARI={max_ari:.10f}") + print( + f"Tie detected at comparison {i} (features {feat_i} vs {feat_j}): {num_ties} positions with ARI={max_ari:.10f}" + ) print(f" GPU chose: {tuple(gpu_parts)}, CPU chose: {tuple(cpu_parts)}") - - + # @pytest.mark.parametrize( # "seed", [42] @@ -376,4 +386,3 @@ def test_ccc_gpu_with_numerical_input( # res_cpu = ccc(df, n_jobs=n_cpu_cores) # res_gpu = ccc_gpu(df) # assert np.allclose(res_cpu, res_gpu) - diff --git a/tests/gpu/test_coef_gpu.py b/tests/gpu/test_coef_gpu.py index 5e50a3c6..9e277fdb 100644 --- a/tests/gpu/test_coef_gpu.py +++ b/tests/gpu/test_coef_gpu.py @@ -1,6 +1,6 @@ +import ccc_cuda_ext import numpy as np import pandas as pd -import ccc_cuda_ext import pytest from ccc.coef.impl_gpu import ccc as ccc_gpu diff --git a/tests/gpu/test_cuda_correctness.py b/tests/gpu/test_cuda_correctness.py index 66d68cbe..a003398d 100644 --- a/tests/gpu/test_cuda_correctness.py +++ b/tests/gpu/test_cuda_correctness.py @@ -20,16 +20,14 @@ change's HANDOFF notes. """ +import ccc_cuda_ext import numpy as np import pandas as pd import pytest - -import ccc_cuda_ext -from ccc.coef.impl_gpu import ccc as ccc_gpu from ccc.coef.impl import ccc as ccc_cpu +from ccc.coef.impl_gpu import ccc as ccc_gpu from utils import clean_gpu_memory - # --------------------------------------------------------------------------- # 1. Input validation -> Python exception, no interpreter crash # --------------------------------------------------------------------------- diff --git a/tests/gpu/test_titanic_dataset.py b/tests/gpu/test_titanic_dataset.py index e2833800..ccca7571 100644 --- a/tests/gpu/test_titanic_dataset.py +++ b/tests/gpu/test_titanic_dataset.py @@ -1,9 +1,9 @@ -import pandas as pd import numpy as np +import pandas as pd import pytest -from scipy.spatial.distance import squareform from ccc.coef.impl import ccc from ccc.coef.impl_gpu import ccc as ccc_gpu +from scipy.spatial.distance import squareform from utils import clean_gpu_memory @@ -59,23 +59,30 @@ def test_ccc_gpu_with_titanic_dataset(titanic_data): # Verify results try: # Use numpy testing instead of pandas since these are numpy arrays - np.testing.assert_allclose(cpu_corrs, gpu_corrs, rtol=1e-5, atol=1e-6, - equal_nan=True, - err_msg="CPU and GPU correlation matrices should be nearly identical") + np.testing.assert_allclose( + cpu_corrs, + gpu_corrs, + rtol=1e-5, + atol=1e-6, + equal_nan=True, + err_msg="CPU and GPU correlation matrices should be nearly identical", + ) except AssertionError as e: # Print more detailed error information print("\nDetailed comparison:") print("Maximum absolute difference:", np.nanmax(np.abs(cpu_corrs - gpu_corrs))) print("Mean absolute difference:", np.nanmean(np.abs(cpu_corrs - gpu_corrs))) - + # Check NaN patterns cpu_nan_mask = np.isnan(cpu_corrs) gpu_nan_mask = np.isnan(gpu_corrs) print(f"CPU NaN count: {np.sum(cpu_nan_mask)}") print(f"GPU NaN count: {np.sum(gpu_nan_mask)}") - + if not np.array_equal(cpu_nan_mask, gpu_nan_mask): print("ERROR: NaN patterns differ between CPU and GPU!") - print("This indicates categorical data processing issues in GPU implementation") - + print( + "This indicates categorical data processing issues in GPU implementation" + ) + raise e diff --git a/tests/gpu/utils.py b/tests/gpu/utils.py index 8e8cc975..c933bcfe 100644 --- a/tests/gpu/utils.py +++ b/tests/gpu/utils.py @@ -1,7 +1,8 @@ import functools + import cupy as cp -import pandas as pd import numpy as np +import pandas as pd def clean_gpu_memory(func): @@ -12,11 +13,19 @@ def wrapper(*args, **kwargs): finally: mempool = cp.get_default_memory_pool() mempool.free_all_blocks() + return wrapper -def generate_categorical_data(n_features, n_objects, n_categories=3, categories=None, - str_length=None, random_state=None, feature_names=None): +def generate_categorical_data( + n_features, + n_objects, + n_categories=3, + categories=None, + str_length=None, + random_state=None, + feature_names=None, +): """ Generate random categorical data as a pandas DataFrame. @@ -61,10 +70,10 @@ def generate_categorical_data(n_features, n_objects, n_categories=3, categories= if categories is None: if str_length is not None: # Generate random string categories - letters = np.array(list('ABCDEFGHIJKLMNOPQRSTUVWXYZ')) + letters = np.array(list("ABCDEFGHIJKLMNOPQRSTUVWXYZ")) categories = [] for _ in range(n_categories): - cat = ''.join(np.random.choice(letters, size=str_length)) + cat = "".join(np.random.choice(letters, size=str_length)) categories.append(cat) else: categories = list(range(n_categories)) @@ -75,13 +84,15 @@ def generate_categorical_data(n_features, n_objects, n_categories=3, categories= random_indices = np.random.randint(0, n_categories, size=(n_objects, n_features)) # Convert indices to categories - categorical_data = np.array([[categories[idx] for idx in row] for row in random_indices]) + categorical_data = np.array( + [[categories[idx] for idx in row] for row in random_indices] + ) # Create feature names if not provided if feature_names is None: - feature_names = [f'feature_{i}' for i in range(n_features)] + feature_names = [f"feature_{i}" for i in range(n_features)] # Create DataFrame df = pd.DataFrame(categorical_data, columns=feature_names) - + return df diff --git a/tests/test_coef.py b/tests/test_coef.py index 1f2dc8bf..ccf14617 100644 --- a/tests/test_coef.py +++ b/tests/test_coef.py @@ -1,28 +1,26 @@ +import os +import time from concurrent.futures import ThreadPoolExecutor from random import shuffle from unittest.mock import patch -import time -import os import numpy as np import pandas as pd import pytest -from sklearn.preprocessing import minmax_scale -from sklearn.metrics import adjusted_rand_score as ari - from ccc.coef import ( ccc, - get_range_n_clusters, - run_quantile_clustering, - get_perc_from_k, - get_parts, - get_coords_from_index, cdist_parts_basic, cdist_parts_parallel, get_chunks, + get_coords_from_index, get_n_workers, + get_parts, + get_perc_from_k, + get_range_n_clusters, + run_quantile_clustering, ) - +from sklearn.metrics import adjusted_rand_score as ari +from sklearn.preprocessing import minmax_scale IN_GITHUB_ACTIONS = os.getenv("GITHUB_ACTIONS") == "true" @@ -796,6 +794,7 @@ def test_cm_values_equal_to_original_implementation(): # implementation (https://github.com/sinc-lab/clustermatch) plus some # patches (see tests/data/README.md about ccc data). from pathlib import Path + import pandas as pd # from pandas.testing import assert_frame_equal diff --git a/tests/test_coef_pval.py b/tests/test_coef_pval.py index 67a0c861..dabd6c00 100644 --- a/tests/test_coef_pval.py +++ b/tests/test_coef_pval.py @@ -4,9 +4,8 @@ import numpy as np import pandas as pd import pytest -from sklearn.preprocessing import minmax_scale - from ccc.coef import ccc +from sklearn.preprocessing import minmax_scale IN_GITHUB_ACTIONS = os.getenv("GITHUB_ACTIONS") == "true" @@ -276,9 +275,9 @@ def test_cm_large_n_objects_pvalue_permutations_is_parallelized(): # Validate that parallel execution is consistently faster for i in range(n_runs): - assert ( - multi_thread_times[i] < single_thread_times[i] - ), f"Run {i+1}: Multi-thread ({multi_thread_times[i]:.2f}s) not faster than single-thread ({single_thread_times[i]:.2f}s)" + assert multi_thread_times[i] < single_thread_times[i], ( + f"Run {i + 1}: Multi-thread ({multi_thread_times[i]:.2f}s) not faster than single-thread ({single_thread_times[i]:.2f}s)" + ) def test_cm_return_parts_quadratic_pvalue(): @@ -335,7 +334,7 @@ def test_cm_numerical_and_categorical_features_perfect_relationship_pvalue(): numerical_feature0_median = np.percentile(numerical_feature0, 50) # create a categorical variable perfectly correlated with the numerical one (this is actually an ordinal feature) - categorical_feature1 = np.full(numerical_feature0.shape[0], "", dtype=np.unicode_) + categorical_feature1 = np.full(numerical_feature0.shape[0], "", dtype=np.str_) categorical_feature1[numerical_feature0 < numerical_feature0_median] = "l" categorical_feature1[numerical_feature0 >= numerical_feature0_median] = "u" _unique_values = np.unique(categorical_feature1) @@ -380,7 +379,7 @@ def test_cm_numerical_and_categorical_features_weakly_relationship_pvalue(): numerical_feature0_perc = np.percentile(numerical_feature0, 2) # create a categorical variable strongly correlated with the numerical one - categorical_feature1 = np.full(numerical_feature0.shape[0], "", dtype=np.unicode_) + categorical_feature1 = np.full(numerical_feature0.shape[0], "", dtype=np.str_) categorical_feature1[numerical_feature0 < numerical_feature0_perc] = "l" categorical_feature1[numerical_feature0 >= numerical_feature0_perc] = "u" _unique_values = np.unique(categorical_feature1) diff --git a/tests/test_conf.py b/tests/test_conf.py index 9f4e8ac1..d3fb2621 100644 --- a/tests/test_conf.py +++ b/tests/test_conf.py @@ -1,9 +1,10 @@ """ Tests the conf.py module. """ + import os -import sys import runpy +import sys from unittest import mock import pytest @@ -18,9 +19,10 @@ def test_conf_module_load(): @mock.patch.dict(os.environ, {}, clear=True) def test_conf_entries(): - from ccc import conf import importlib + from ccc import conf + importlib.reload(conf) assert conf.ROOT_DIR is not None @@ -57,8 +59,9 @@ def test_conf_main(): reason="exporting variables is only supported in non-Windows platforms", ) def test_conf_export_variables(): - from pathlib import Path import subprocess + from pathlib import Path + from ccc import conf conf_filepath = Path(conf.__file__).resolve() @@ -103,9 +106,10 @@ def test_conf_export_variables(): @mock.patch.dict(os.environ, {"CM_MANUSCRIPT_DIR": "/tmp/some/dir"}) def test_conf_with_manuscript_dir(): - from ccc import conf import importlib + from ccc import conf + importlib.reload(conf) assert conf.MANUSCRIPT is not None @@ -116,9 +120,10 @@ def test_conf_with_manuscript_dir(): @mock.patch.dict(os.environ, {"CM_N_JOBS": ""}) def test_conf_cm_n_jobs_is_empty_string(): - from ccc import conf import importlib + from ccc import conf + importlib.reload(conf) assert conf.GENERAL is not None diff --git a/tests/test_corr.py b/tests/test_corr.py index 14b206dc..4743e4d6 100644 --- a/tests/test_corr.py +++ b/tests/test_corr.py @@ -1,9 +1,9 @@ """ Tests the corr.py module. """ + import numpy as np import pandas as pd - from ccc import corr @@ -250,6 +250,7 @@ def test_corr_clustermatch_outputs_same_as_original_clustermatch(): # implementation (https://github.com/sinc-lab/clustermatch) plus some # patches (see README.md in tests/data about ccc data). from pathlib import Path + from pandas.testing import assert_frame_equal input_data_dir = Path(__file__).parent / "data" diff --git a/tests/test_giant.py b/tests/test_giant.py index 2c41289e..8ceadb60 100644 --- a/tests/test_giant.py +++ b/tests/test_giant.py @@ -9,8 +9,7 @@ "Skipping REST test on GIANT in non-Linux systems", allow_module_level=True ) -from ccc.giant import gene_exists, predict_tissue, get_network - +from ccc.giant import gene_exists, get_network, predict_tissue # Gene mappings used in unit tests gene_mappings = pd.DataFrame( diff --git a/tests/test_methods.py b/tests/test_methods.py index c59e3cc5..b2c51451 100644 --- a/tests/test_methods.py +++ b/tests/test_methods.py @@ -1,5 +1,4 @@ import numpy as np - from ccc.methods import mic diff --git a/tests/test_plots.py b/tests/test_plots.py index 8cda7a5f..dfa66b8e 100644 --- a/tests/test_plots.py +++ b/tests/test_plots.py @@ -1,4 +1,5 @@ import sys + import pytest if sys.platform.startswith("win"): @@ -9,12 +10,11 @@ import numpy as np import pandas as pd - from ccc.plots import ( - plot_histogram, - plot_cumulative_histogram, - jointplot, MyUpSet, + jointplot, + plot_cumulative_histogram, + plot_histogram, ) diff --git a/tests/test_pytorch_core.py b/tests/test_pytorch_core.py index a74a92eb..4dfe5696 100644 --- a/tests/test_pytorch_core.py +++ b/tests/test_pytorch_core.py @@ -1,6 +1,5 @@ -import pytest import numpy as np - +import pytest from ccc.pytorch.core import unravel_index_2d diff --git a/tests/test_scipy_stats.py b/tests/test_scipy_stats.py index 34f059f4..0f45fe39 100644 --- a/tests/test_scipy_stats.py +++ b/tests/test_scipy_stats.py @@ -1,7 +1,6 @@ import numpy as np -from scipy import stats - from ccc.scipy.stats import rank +from scipy import stats def test_rank_no_duplicates(): diff --git a/tests/test_sklearn_metrics.py b/tests/test_sklearn_metrics.py index cd3e8a09..7d8acfcb 100644 --- a/tests/test_sklearn_metrics.py +++ b/tests/test_sklearn_metrics.py @@ -1,11 +1,10 @@ import numpy as np -from sklearn.metrics import adjusted_rand_score as sklearn_ari - from ccc.sklearn.metrics import ( adjusted_rand_index, get_contingency_matrix, get_pair_confusion_matrix, ) +from sklearn.metrics import adjusted_rand_score as sklearn_ari def test_get_contingency_matrix_k0_equal_k1(): diff --git a/tests/test_utils.py b/tests/test_utils.py index 131b9b89..c0950fa8 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,9 +1,10 @@ """ Tests the utility_functions.py module. """ + import sys -from unittest.mock import MagicMock from pathlib import Path +from unittest.mock import MagicMock import numpy as np import pandas as pd @@ -18,8 +19,8 @@ def reload_package(root_module): Taken and adapted from here: https://stackoverflow.com/a/2918951/3120414 """ - from importlib import import_module import types + from importlib import import_module package_name = root_module.__name__