diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 0000000..ab178bc --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,67 @@ +name: Build and Deploy Documentation + +on: + push: + branches: + - main + pull_request: + branches: + - main + workflow_dispatch: + +permissions: + contents: read + pages: write + id-token: write + +concurrency: + group: "pages" + cancel-in-progress: false + +jobs: + build: + runs-on: ubuntu-latest + + steps: + - name: Checkout repository + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.11' + + - name: Install dependencies + run: | + sudo apt-get update -qq + sudo apt-get install -y -qq pandoc + python -m pip install --upgrade pip + pip install numpy==2.3.5 numba==0.63.1 + pip install sphinx>=7.0.0 sphinx-rtd-theme>=2.0.0 nbsphinx>=0.9.0 ipython>=8.0.0 + + - name: Build documentation + run: | + cd docs + make html + touch build/html/.nojekyll + + - name: Upload artifact + uses: actions/upload-pages-artifact@v3 + with: + path: 'docs/build/html' + + deploy: + # Only deploy on push to main branch + if: github.event_name == 'push' && github.ref == 'refs/heads/main' + + environment: + name: github-pages + url: ${{ steps.deployment.outputs.page_url }} + + runs-on: ubuntu-latest + needs: build + + steps: + - name: Deploy to GitHub Pages + id: deployment + uses: actions/deploy-pages@v4 diff --git a/.gitignore b/.gitignore index 78fdef6..ff63761 100644 --- a/.gitignore +++ b/.gitignore @@ -26,6 +26,13 @@ venv/ ENV/ env/ +# Documentation builds +docs/build/ +docs/source/_build/ + +# Jupyter Notebook checkpoints +.ipynb_checkpoints + # IDE .vscode/ .idea/ diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d0c3cbf --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = source +BUILDDIR = build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 0000000..472cd75 --- /dev/null +++ b/docs/README.md @@ -0,0 +1,147 @@ +# Documentation Build Guide + +This guide explains how to build the CudAkima documentation locally. + +## Prerequisites + +You need to install the documentation dependencies: + +```bash +# Using pip +pip install -e ".[docs]" + +# Or install dependencies individually +pip install sphinx>=7.0.0 sphinx-rtd-theme>=2.0.0 nbsphinx>=0.9.0 ipython>=8.0.0 +``` + +## Building the Documentation + +To build the HTML documentation: + +```bash +cd docs +make html +``` + +The generated documentation will be in `docs/build/html/`. Open `docs/build/html/index.html` in your browser to view it. + +## Cleaning Build Files + +To clean the build directory: + +```bash +cd docs +make clean +``` + +## Other Output Formats + +Sphinx supports multiple output formats. Some useful ones: + +```bash +# PDF (requires LaTeX) +make latexpdf + +# Plain text +make text + +# ePub +make epub +``` + +## GitHub Pages Deployment + +The documentation is automatically built and deployed to GitHub Pages when changes are pushed to the `main` branch. The workflow is defined in `.github/workflows/docs.yml`. + +To enable GitHub Pages: + +1. Go to your repository settings on GitHub +2. Navigate to "Pages" in the left sidebar +3. Under "Source", select "GitHub Actions" +4. The documentation will be available at `https://.github.io//` + +## Documentation Structure + +``` +docs/ +├── source/ +│ ├── conf.py # Sphinx configuration +│ ├── index.rst # Main documentation page +│ ├── api.rst # API reference +│ ├── tutorial.rst # Tutorial page +│ ├── examples.rst # Examples page +│ ├── notebooks/ # Jupyter notebooks +│ │ └── tutorial.ipynb (symlink to examples/tutorial.ipynb) +│ ├── _static/ # Static files (CSS, images, etc.) +│ └── _templates/ # Custom templates +├── build/ # Generated documentation (git-ignored) +└── Makefile # Build commands +``` + +## Updating Documentation + +1. **Class/Function Documentation**: Update docstrings in the source code (`src/cudakima/`) +2. **Tutorial**: Edit `examples/tutorial.ipynb` +3. **Examples**: Edit `docs/source/examples.rst` +4. **Main Pages**: Edit `.rst` files in `docs/source/` + +After making changes, rebuild the documentation to see the updates: + +```bash +cd docs +make clean +make html +``` + +## Docstring Format + +CudAkima uses NumPy-style docstrings. Here's an example: + +```python +def my_function(x, y): + """ + Short description. + + Longer description if needed. + + Parameters + ---------- + x : array_like + Description of x + y : float + Description of y + + Returns + ------- + result : ndarray + Description of return value + + Examples + -------- + >>> my_function([1, 2, 3], 2.0) + array([2., 4., 6.]) + """ + pass +``` + +## Troubleshooting + +### Import Errors + +If you get import errors when building the documentation, make sure the package is installed: + +```bash +pip install -e . +``` + +### Notebook Execution Errors + +The notebooks are set to not execute during the build (`nbsphinx_execute = 'never'`). If you want to execute them, change this setting in `docs/source/conf.py`. + +### Missing Dependencies + +If you get warnings about missing dependencies, install them: + +```bash +pip install -e ".[docs]" +``` diff --git a/docs/source/.nojekyll b/docs/source/.nojekyll new file mode 100644 index 0000000..e69de29 diff --git a/docs/source/api.rst b/docs/source/api.rst new file mode 100644 index 0000000..f4a428e --- /dev/null +++ b/docs/source/api.rst @@ -0,0 +1,90 @@ +API Reference +============= + +This page documents the complete API for CudAkima. + +Main Classes +------------ + +.. currentmodule:: cudakima + +AkimaInterpolant1D +~~~~~~~~~~~~~~~~~~ + +.. autoclass:: cudakima.AkimaInterpolant1D + :members: + :special-members: __init__, __call__ + :show-inheritance: + +AkimaInterpolant1DMultiDim +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: cudakima.AkimaInterpolant1DMultiDim + :members: + :special-members: __init__, __call__ + :show-inheritance: + +AkimaInterpolant1DFlexible +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: cudakima.AkimaInterpolant1DFlexible + :members: + :special-members: __init__, __call__ + :show-inheritance: + +Kernel Functions +---------------- + +This section documents the low-level numba kernels used for GPU and CPU computation. + +GPU Kernels +~~~~~~~~~~~ + +.. currentmodule:: cudakima.kernels + +.. autofunction:: linearslope_gpu + +.. autofunction:: splineslope_gpu + +.. autofunction:: splineslope_gpu_optimized + +.. autofunction:: binary_search_gpu + +.. autofunction:: akima_spline_kernel_gpu + +.. autofunction:: akima_spline_kernel_optimized + +.. autofunction:: akima_linear_kernel + +.. autofunction:: precompute_slopes_kernel + +.. autofunction:: precompute_spline_slopes_kernel + +.. autofunction:: akima_linear_kernel_gpu_multidim + +.. autofunction:: akima_spline_kernel_gpu_multidim + +CPU Kernels +~~~~~~~~~~~ + +.. autofunction:: linearslope_cpu + +.. autofunction:: splineslope_cpu + +.. autofunction:: splineslope_cpu_optimized + +.. autofunction:: binary_search_cpu + +.. autofunction:: akima_spline_kernel_cpu + +.. autofunction:: akima_spline_kernel_cpu_optimized + +.. autofunction:: akima_linear_kernel_cpu + +.. autofunction:: precompute_all_linear_slopes_cpu + +.. autofunction:: precompute_all_spline_slopes_cpu + +.. autofunction:: akima_linear_kernel_cpu_multidim + +.. autofunction:: akima_spline_kernel_cpu_multidim diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 0000000..330db6b --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,88 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +import os +import sys +sys.path.insert(0, os.path.abspath('../../src')) + +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information + +project = 'CudAkima' +copyright = '2024, Alessandro Santini' +author = 'Alessandro Santini' +release = '0.2.2' +version = '0.2.2' + +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration + +extensions = [ + 'sphinx.ext.autodoc', + 'sphinx.ext.napoleon', + 'sphinx.ext.viewcode', + 'sphinx.ext.intersphinx', + 'sphinx.ext.mathjax', + 'nbsphinx', + 'sphinx_rtd_theme', +] + +templates_path = ['_templates'] +exclude_patterns = [] + +# Napoleon settings for Google/NumPy style docstrings +napoleon_google_docstring = True +napoleon_numpy_docstring = True +napoleon_include_init_with_doc = True +napoleon_include_private_with_doc = False +napoleon_include_special_with_doc = True +napoleon_use_admonition_for_examples = False +napoleon_use_admonition_for_notes = False +napoleon_use_admonition_for_references = False +napoleon_use_ivar = False +napoleon_use_param = True +napoleon_use_rtype = True +napoleon_preprocess_types = False +napoleon_type_aliases = None +napoleon_attr_annotations = True + +# Intersphinx mapping +intersphinx_mapping = { + 'python': ('https://docs.python.org/3', None), + 'numpy': ('https://numpy.org/doc/stable/', None), + 'numba': ('https://numba.readthedocs.io/en/stable/', None), +} + +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output + +html_theme = 'sphinx_rtd_theme' +html_static_path = ['_static'] +html_theme_options = { + 'logo_only': False, + 'display_version': True, + 'prev_next_buttons_location': 'bottom', + 'style_external_links': False, + 'collapse_navigation': False, + 'sticky_navigation': True, + 'navigation_depth': 4, + 'includehidden': True, + 'titles_only': False +} + +# -- Options for nbsphinx ---------------------------------------------------- +nbsphinx_execute = 'never' # Don't execute notebooks during build +nbsphinx_allow_errors = True # Continue building even if notebook has errors + +# -- Options for autodoc ----------------------------------------------------- +autodoc_member_order = 'bysource' +autodoc_typehints = 'description' +autodoc_default_options = { + 'members': True, + 'member-order': 'bysource', + 'special-members': '__init__, __call__', + 'undoc-members': True, + 'exclude-members': '__weakref__' +} diff --git a/docs/source/examples.rst b/docs/source/examples.rst new file mode 100644 index 0000000..cbf8952 --- /dev/null +++ b/docs/source/examples.rst @@ -0,0 +1,127 @@ +Examples +======== + +Additional Examples +------------------- + +For more examples and benchmarks, check out the `examples directory `_ +in the repository. + +Basic Interpolation +~~~~~~~~~~~~~~~~~~~ + +Single Array Interpolation +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: python + + from cudakima import AkimaInterpolant1D + import numpy as np + + # Simple single array interpolation + x = np.array([1, 2, 3, 4, 5]) + y = np.array([1, 4, 9, 16, 25]) + + interpolant = AkimaInterpolant1D() + x_new = np.linspace(1, 5, 50) + y_new = interpolant(x_new, x, y) + +Batch Interpolation with Different Lengths +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: python + + from cudakima import AkimaInterpolant1D + import numpy as np + + # Batch interpolation with NaN padding + x = np.array([ + [1, 2, 3, 4, 5, np.nan, np.nan], + [1, 2, 3, 4, 5, 6, np.nan], + [1, 2, 3, 4, 5, 6, 7] + ]) + y = np.array([ + [1, 4, 9, 16, 25, np.nan, np.nan], + [1, 4, 9, 16, 25, 36, np.nan], + [1, 4, 9, 16, 25, 36, 49] + ]) + + interpolant = AkimaInterpolant1D() + x_new = np.linspace(1, 5, 100) + y_new = interpolant(x_new, x, y) # Shape: (3, 100) + +Advanced Usage +~~~~~~~~~~~~~~ + +Linear Interpolation +^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: python + + from cudakima import AkimaInterpolant1D + + # Use linear interpolation instead of cubic + interpolant = AkimaInterpolant1D(order='linear') + y_new = interpolant(x_new, x, y) + +CPU-Only Mode +^^^^^^^^^^^^^ + +.. code-block:: python + + from cudakima import AkimaInterpolant1D + + # Force CPU execution (useful for debugging or when GPU is unavailable) + interpolant = AkimaInterpolant1D(use_gpu=False) + y_new = interpolant(x_new, x, y) + +Multidimensional x_new +^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: python + + from cudakima import AkimaInterpolant1DMultiDim + import numpy as np + + # Different interpolation points for each group + x = np.random.rand(2, 3, 10) # 2x3 batch, 10 points each + y = np.sin(x) + x_new = np.random.rand(2, 3, 50) # 2x3 batch, 50 interp points each + + interpolant = AkimaInterpolant1DMultiDim() + y_new = interpolant(x_new, x, y) # Shape: (2, 3, 50) + +Flexible Interpolator +^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: python + + from cudakima import AkimaInterpolant1DFlexible + + # Automatically chooses between standard and multidimensional modes + interpolant = AkimaInterpolant1DFlexible() + + # Standard mode (1D x_new) + y_new_standard = interpolant(x_new_1d, x, y) + + # Multidimensional mode (matching batch dims) + y_new_multidim = interpolant(x_new_multidim, x, y) + +Performance Tips +~~~~~~~~~~~~~~~~ + +1. **Presorting**: If your input data is already sorted, set ``sanitize=False`` to skip sorting: + + .. code-block:: python + + interpolant = AkimaInterpolant1D(sanitize=False) + +2. **Batch Size**: For optimal GPU performance, use larger batch sizes (more groups to interpolate). + +3. **Thread Configuration**: Adjust ``threadsperblock`` for your GPU (default is 64): + + .. code-block:: python + + interpolant = AkimaInterpolant1D(threadsperblock=128) + +4. **GPU Memory**: For very large datasets, consider processing in chunks to avoid out-of-memory errors. diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 0000000..2f4c73e --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,100 @@ +.. CudAkima documentation master file + +Welcome to CudAkima's Documentation +=================================== + +CudAkima is a Python package that provides a parallel, GPU-accelerated implementation +of Akima Splines. The package also includes CPU support for systems without CUDA/CuPy. + +**Key Features:** + +* GPU-accelerated parallel Akima spline interpolation +* CPU fallback for systems without GPU support +* Efficient batch interpolation of arrays with different lengths +* Support for both linear and cubic (Akima spline) interpolation +* Multidimensional interpolation support + +Quick Start +----------- + +Installation +~~~~~~~~~~~~ + +CudAkima requires ``numpy`` and ``numba``. For GPU support, you'll also need ``cupy``. + +.. code-block:: bash + + git clone https://github.com/asantini29/CudAkima.git + cd CudAkima + uv sync + +Basic Usage +~~~~~~~~~~~ + +.. code-block:: python + + from cudakima import AkimaInterpolant1D + import numpy as np + + # Create sample data + x = np.array([[1, 2, 3, 4, 5, np.nan, np.nan], + [1, 2, 3, 4, 5, 6, np.nan], + [1, 2, 3, 4, 5, 6, 7]]) + y = np.array([[1, 4, 9, 16, 25, np.nan, np.nan], + [1, 4, 9, 16, 25, 36, np.nan], + [1, 4, 9, 16, 25, 36, 49]]) + + # Create interpolator + interpolant = AkimaInterpolant1D() + + # Interpolate + x_new = np.linspace(1, 5, 100) + y_new = interpolant(x_new, x, y) + +About Akima Splines +------------------- + +`Akima Splines `_ are spline interpolants +that tend to show smoother behavior with respect to the widely used Cubic Splines. +Unlike cubic splines, Akima splines have discontinuous second derivatives, which can +be advantageous in certain applications. + +Why CudAkima? +~~~~~~~~~~~~~ + +While both ``scipy`` and ``cupy`` offer implementations of Akima splines, they only +support 1D x-arrays for interpolation. **CudAkima** enables fast, parallel interpolation +of batches of arrays with different lengths by padding shorter arrays with NaN values +and stacking them in multidimensional arrays. + +This makes CudAkima particularly suited for applications where: + +* The arrays to interpolate keep changing (e.g., parameter estimation) +* You need to interpolate many arrays in parallel +* Arrays in the batch have different lengths + +Performance +~~~~~~~~~~~ + +On CPU, CudAkima is approximately **3x faster** than a naive loop using scipy. +On GPU, CudAkima is approximately **20x faster** than using cupy in a loop. + +See the tutorial notebook in the examples directory for detailed benchmarks. + +Table of Contents +----------------- + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + api + tutorial + examples + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/source/notebooks/tutorial.ipynb b/docs/source/notebooks/tutorial.ipynb new file mode 120000 index 0000000..01fe282 --- /dev/null +++ b/docs/source/notebooks/tutorial.ipynb @@ -0,0 +1 @@ +/home/runner/work/CudAkima/CudAkima/examples/tutorial.ipynb \ No newline at end of file diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst new file mode 100644 index 0000000..874b4c7 --- /dev/null +++ b/docs/source/tutorial.rst @@ -0,0 +1,9 @@ +Tutorial +======== + +This tutorial demonstrates the usage of CudAkima with practical examples. + +.. toctree:: + :maxdepth: 1 + + notebooks/tutorial diff --git a/pyproject.toml b/pyproject.toml index 873f29c..8258dd7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,6 +13,14 @@ dependencies = [ "numba==0.63.1", ] +[project.optional-dependencies] +docs = [ + "sphinx>=7.0.0", + "sphinx-rtd-theme>=2.0.0", + "nbsphinx>=0.9.0", + "ipython>=8.0.0", +] + [build-system] requires = ["hatchling"] build-backend = "hatchling.build" diff --git a/src/cudakima/akima.py b/src/cudakima/akima.py index 1988577..2086165 100644 --- a/src/cudakima/akima.py +++ b/src/cudakima/akima.py @@ -48,25 +48,73 @@ class AkimaInterpolant1D(): """ GPU-accelerated parallel Akima Splines. - This class provides a parallel implementation of the Akima spline interpolation algorithm. It has a CPU version as well. - The interpolation can be performed simultaneously on array of different sizes. This is achieved by passing as inputs multidimensional arrays padded with NaN values. - The interpolation is performed along the last axis of the input arrays, which must have dimension equal to the one of the longest array in the batch. - - Example: - If you have a batch of 3 arrays with 5, 6 and 7 points respectively, you can interpolate them all at once by passing a 3D array with shape=(3, 7) and padding the arrays with NaN values.: - ``` - x = np.array([[1, 2, 3, 4, 5, np.nan, np.nan], [1, 2, 3, 4, 5, 6, np.nan], [1, 2, 3, 4, 5, 6, 7]]) - y = np.array([[1, 4, 9, 16, 25, np.nan, np.nan], [1, 4, 9, 16, 25, 36, np.nan], [1, 4, 9, 16, 25, 36, 49]]) - x_new = np.linspace(1, 5, 100) - interpolant = AkimaInterpolant1D() - y_new = interpolant(x_new, x, y) - ``` - - Parameters: - use_gpu (bool): If True, the interpolation is performed on the GPU if available. If False, the interpolation is performed on the CPU. - threadsperblock (int): The number of threads per block to use for the GPU implementation. This parameter is ignored if `use_gpu` is False. - sanitize (bool): If True, the input data are sorted in ascending order. If False, the input data must be already sorted. - verbose (bool): If True, print information about the interpolation process. Default is False. + This class provides a parallel implementation of the Akima spline interpolation + algorithm with both GPU and CPU support. The interpolation can be performed + simultaneously on arrays of different sizes by passing multidimensional arrays + padded with NaN values. The interpolation is performed along the last axis of + the input arrays. + + Parameters + ---------- + use_gpu : bool, optional + If True, use GPU acceleration if CUDA and CuPy are available. If False or + if GPU is unavailable, fall back to CPU implementation. Default is True. + threadsperblock : int, optional + Number of threads per block for GPU kernel execution. Ignored if using CPU. + Default is 64. + order : {'linear', 'cubic'}, optional + Interpolation order. 'cubic' uses Akima spline interpolation, 'linear' uses + linear interpolation. Default is 'cubic'. + sanitize : bool, optional + If True, sort input data in ascending order. If False, assumes input data + are already sorted. Set to False to improve performance when data is + pre-sorted. Default is False. + verbose : bool, optional + If True, print information about the interpolation process. Default is False. + + Attributes + ---------- + xp : module + Either numpy or cupy, depending on GPU availability + order : str + The interpolation order being used + sanitize : bool + Whether input sanitization is enabled + threadsperblock : int + Number of threads per GPU block + + Examples + -------- + Basic usage with batch interpolation: + + >>> import numpy as np + >>> from cudakima import AkimaInterpolant1D + >>> + >>> # Batch of 3 arrays with different lengths (5, 6, 7 points) + >>> x = np.array([[1, 2, 3, 4, 5, np.nan, np.nan], + ... [1, 2, 3, 4, 5, 6, np.nan], + ... [1, 2, 3, 4, 5, 6, 7]]) + >>> y = np.array([[1, 4, 9, 16, 25, np.nan, np.nan], + ... [1, 4, 9, 16, 25, 36, np.nan], + ... [1, 4, 9, 16, 25, 36, 49]]) + >>> + >>> # Create interpolator and interpolate + >>> interpolant = AkimaInterpolant1D() + >>> x_new = np.linspace(1, 5, 100) + >>> y_new = interpolant(x_new, x, y) # Shape: (3, 100) + + Using linear interpolation on CPU: + + >>> interpolant = AkimaInterpolant1D(use_gpu=False, order='linear') + >>> y_new = interpolant(x_new, x, y) + + Notes + ----- + - Requires at least 4 finite points for Akima spline interpolation due to + boundary conditions. Falls back to linear interpolation for fewer points. + - NaN values must be placed at the end of each array in the batch. + - All arrays in a batch must be padded to the same length (the length of the + longest array). """ def __init__(self, use_gpu=True, threadsperblock=64, order='cubic', sanitize=False, verbose=False): @@ -106,15 +154,26 @@ def sanitize(self, value): def sort_input(self, x, y): """ - Check that the input data have the right shape and sort them to ensure that all the arrays are in the correct order. - - Parameters: - x (ndarray): The x-values of the data points. Shape=(..., n). - y (ndarray): The y-values of the data points. Shape=(..., n). - - Returns: - x (ndarray): The sorted x values. - y (ndarray): The y values, sorted accordingly to `x`. + Sort input data in ascending order along the last axis. + + Parameters + ---------- + x : array_like + The x-values of the data points. Shape (..., n). + y : array_like + The y-values of the data points. Shape (..., n). + + Returns + ------- + x_sorted : ndarray + The sorted x values. + y_sorted : ndarray + The y values, sorted according to x. + + Notes + ----- + This method ensures that x values are in ascending order, which is required + for the interpolation algorithm to work correctly. """ x = self.xp.asarray(x) y = self.xp.asarray(y) @@ -129,17 +188,29 @@ def sort_input(self, x, y): def pass_input(self, x, y): - """ - Skip the sorting of the input data. This is useful to save computational time when the input data are already sorted. - - Parameters: - x (ndarray): The x-values of the data points. Shape=(..., n). - y (ndarray): The y-values of the data points. Shape=(..., n). - - Returns: - x (ndarray): The input `x` array. - y (ndarray): The input `y` array. + Pass input data through without sorting. + + This is used when sanitize=False to skip the sorting step for performance. + + Parameters + ---------- + x : array_like + The x-values of the data points. Shape (..., n). + y : array_like + The y-values of the data points. Shape (..., n). + + Returns + ------- + x : ndarray + The input x array (unchanged). + y : ndarray + The input y array (unchanged). + + Notes + ----- + When using this method, the user must ensure that input data are already + sorted in ascending order with NaN values at the end. """ return x, y @@ -156,18 +227,42 @@ def set_sanitize(self): def __call__(self, x_new, x, y, **kwargs): """ - Interpolates the values of `x_new` based on the given `x` and `y` data points. - - Parameters: - x_new (ndarray): The new x-values to interpolate. If `self.sanitize` is `False`, they must be sorted in ascending order. Shape=(n_f,). - x (ndarray): The x-values of the data points. If `self.sanitize` is `False`, they must be sorted in ascending order. Shape=(..., n). - If along some axis there are less than `n` points to interpolate, the remaining values must be NaN. - y (ndarray): The y-values of the data points. Shape=(..., n). - If along some axis there are less than `n` points to interpolate, the remaining values must be NaN. - **kwargs: Additional keyword arguments. Added for future flexibility. - - Returns: - ndarray: The interpolated values of `x_new`. Shape=(..., n_f). + Interpolate at new x-values. + + Parameters + ---------- + x_new : array_like, shape (n_f,) + New x-values at which to interpolate. If sanitize=False, must be sorted + in ascending order. + x : array_like, shape (..., n) + X-values of the data points. If sanitize=False, must be sorted in + ascending order along the last axis. Arrays with fewer than n points + should be padded with NaN values at the end. + y : array_like, shape (..., n) + Y-values of the data points. Shape must match x. Arrays with fewer + than n points should be padded with NaN values at the end. + **kwargs : dict, optional + Additional keyword arguments (reserved for future use). + + Returns + ------- + y_new : ndarray, shape (..., n_f) + Interpolated values at x_new positions. + + Examples + -------- + >>> import numpy as np + >>> x = np.array([1, 2, 3, 4, 5]) + >>> y = np.array([1, 4, 9, 16, 25]) + >>> interpolant = AkimaInterpolant1D() + >>> x_new = np.array([1.5, 2.5, 3.5]) + >>> y_new = interpolant(x_new, x, y) + + Notes + ----- + The method automatically handles batches of different-length arrays by + using NaN padding. The interpolation is performed in parallel across + all arrays in the batch. """ x, y = self.sanitize_input(x, y) @@ -192,7 +287,27 @@ def __call__(self, x_new, x, y, **kwargs): def linear_interpolate_gpu(self, x_new, x, y, nin, ngroups, nnans, result): """ - Optimized GPU interpolation with precomputed slopes and binary search + GPU implementation of linear interpolation with precomputed slopes. + + This method performs optimized linear interpolation on the GPU using + precomputed slopes and binary search for interval location. + + Parameters + ---------- + x_new : array_like + New x-values for interpolation + x : array_like + Flattened x data + y : array_like + Flattened y data + nin : int + Number of groups + ngroups : int + Maximum points per group + nnans : array_like + Number of NaN values per group + result : array_like + Output array for interpolated values """ # Allocate temporary arrays for precomputed slopes total_points = nin * ngroups @@ -221,7 +336,27 @@ def linear_interpolate_gpu(self, x_new, x, y, nin, ngroups, nnans, result): def linear_interpolate_cpu(self, x_new, x, y, nin, ngroups, nnans, result): """ - Optimized CPU interpolation with precomputed slopes and parallel execution + CPU implementation of linear interpolation with precomputed slopes. + + This method performs optimized linear interpolation on the CPU using + precomputed slopes, binary search, and parallel execution. + + Parameters + ---------- + x_new : array_like + New x-values for interpolation + x : array_like + Flattened x data + y : array_like + Flattened y data + nin : int + Number of groups + ngroups : int + Maximum points per group + nnans : array_like + Number of NaN values per group + result : array_like + Output array for interpolated values """ # Step 1: Precompute all linear slopes linear_slopes = precompute_all_linear_slopes_cpu(x, y, nin, ngroups, nnans) @@ -232,7 +367,27 @@ def linear_interpolate_cpu(self, x_new, x, y, nin, ngroups, nnans, result): def cubic_interpolate_gpu(self, x_new, x, y, nin, ngroups, nnans, result): """ - Optimized GPU interpolation with precomputed slopes and binary search + GPU implementation of Akima spline interpolation with precomputed slopes. + + This method performs optimized Akima spline interpolation on the GPU using + precomputed linear and spline slopes with binary search for interval location. + + Parameters + ---------- + x_new : array_like + New x-values for interpolation + x : array_like + Flattened x data + y : array_like + Flattened y data + nin : int + Number of groups + ngroups : int + Maximum points per group + nnans : array_like + Number of NaN values per group + result : array_like + Output array for interpolated values """ # Allocate temporary arrays for precomputed slopes total_points = nin * ngroups @@ -266,7 +421,27 @@ def cubic_interpolate_gpu(self, x_new, x, y, nin, ngroups, nnans, result): def cubic_interpolate_cpu(self, x_new, x, y, nin, ngroups, nnans, result): """ - Optimized CPU interpolation with precomputed slopes and parallel execution + CPU implementation of Akima spline interpolation with precomputed slopes. + + This method performs optimized Akima spline interpolation on the CPU using + precomputed linear and spline slopes with binary search and parallel execution. + + Parameters + ---------- + x_new : array_like + New x-values for interpolation + x : array_like + Flattened x data + y : array_like + Flattened y data + nin : int + Number of groups + ngroups : int + Maximum points per group + nnans : array_like + Number of NaN values per group + result : array_like + Output array for interpolated values """ # Step 1: Precompute all linear slopes linear_slopes = precompute_all_linear_slopes_cpu(x, y, nin, ngroups, nnans) @@ -280,44 +455,82 @@ def cubic_interpolate_cpu(self, x_new, x, y, nin, ngroups, nnans, result): class AkimaInterpolant1DMultiDim(AkimaInterpolant1D): """ - GPU-accelerated parallel Akima Splines with multidimensional x_new support. + Akima spline interpolator with multidimensional x_new support. - This class extends the base AkimaInterpolant1D to handle multidimensional x_new arrays, - where each group in the batch can have different interpolation points. + This class extends AkimaInterpolant1D to handle multidimensional x_new arrays, + where each group in the batch can have its own set of interpolation points. - Key difference from base class: - - x_new can have shape (..., n_f) matching the batch dimensions of x and y - - Each group gets its own set of interpolation points - - Output shape remains (..., n_f) as before + The key difference from the base class is that x_new can have shape (..., n_f) + matching the batch dimensions of x and y, allowing different interpolation + points for each group. + + Parameters + ---------- + use_gpu : bool, optional + If True, use GPU acceleration if available. Default is True. + threadsperblock : int, optional + Number of threads per block for GPU execution. Default is 64. + order : {'linear', 'cubic'}, optional + Interpolation order. Default is 'cubic'. + sanitize : bool, optional + If True, sort input data. Default is False. + verbose : bool, optional + If True, print information. Default is False. - Example: - If you have batch data with shape (M, L, K, n) and want different interpolation - points for each group: - ``` - x = np.random.rand(2, 3, 10) # 2x3 batch, 10 points each - y = np.sin(x) # same shape - x_new = np.random.rand(2, 3, 50) # 2x3 batch, 50 interp points each + Examples + -------- + With multidimensional x_new where each group has different interpolation points: - interpolant = AkimaInterpolant1DMultiDim() - y_new = interpolant(x_new, x, y) # shape: (2, 3, 50) - ``` + >>> import numpy as np + >>> from cudakima import AkimaInterpolant1DMultiDim + >>> + >>> # Batch data with shape (2, 3, 10) + >>> x = np.random.rand(2, 3, 10) + >>> y = np.sin(x) + >>> + >>> # Different interpolation points for each group (2, 3, 50) + >>> x_new = np.random.rand(2, 3, 50) + >>> + >>> interpolant = AkimaInterpolant1DMultiDim() + >>> y_new = interpolant(x_new, x, y) # Shape: (2, 3, 50) - Parameters are the same as AkimaInterpolant1D. + Notes + ----- + - x_new must have the same batch dimensions as x and y + - Each group gets its own set of interpolation points + - Inherits all methods and attributes from AkimaInterpolant1D + + See Also + -------- + AkimaInterpolant1D : Base class with 1D x_new + AkimaInterpolant1DFlexible : Automatically chooses between 1D and multidimensional modes """ def __call__(self, x_new, x, y, **kwargs): """ - Interpolates with multidimensional x_new arrays. - - Parameters: - x_new (ndarray): The new x-values to interpolate. Shape=(..., n_f) where ... - matches the batch dimensions of x and y. - x (ndarray): The x-values of the data points. Shape=(..., n). - y (ndarray): The y-values of the data points. Shape=(..., n). - **kwargs: Additional keyword arguments. - - Returns: - ndarray: The interpolated values. Shape=(..., n_f). + Interpolate with multidimensional x_new arrays. + + Parameters + ---------- + x_new : array_like, shape (..., n_f) + New x-values at which to interpolate. Batch dimensions (...) must + match those of x and y. + x : array_like, shape (..., n) + X-values of the data points. + y : array_like, shape (..., n) + Y-values of the data points. + **kwargs : dict, optional + Additional keyword arguments. + + Returns + ------- + y_new : ndarray, shape (..., n_f) + Interpolated values. + + Raises + ------ + ValueError + If x_new batch dimensions don't match x and y batch dimensions. """ # Validate input shapes if x.shape[:-1] != y.shape[:-1]: @@ -467,21 +680,56 @@ def cubic_interpolate_cpu_multidim(self, x_new, x, y, nin, ngroups, nnans, resul class AkimaInterpolant1DFlexible(AkimaInterpolant1D): """ - Flexible Akima interpolator that automatically detects x_new dimensionality. + Flexible Akima interpolator with automatic mode detection. + + This class automatically detects whether x_new is 1D (broadcast across all groups) + or multidimensional (different for each group) and dispatches to the appropriate + interpolation method. + + This provides a unified interface that handles both use cases seamlessly. + + Parameters + ---------- + use_gpu : bool, optional + If True, use GPU acceleration if available. Default is True. + threadsperblock : int, optional + Number of threads per block for GPU execution. Default is 64. + order : {'linear', 'cubic'}, optional + Interpolation order. Default is 'cubic'. + sanitize : bool, optional + If True, sort input data. Default is False. + verbose : bool, optional + If True, print information. Default is False. - This class automatically chooses between standard (1D x_new) and multidimensional - x_new interpolation based on the input array shapes. + Examples + -------- + Works seamlessly with both 1D and multidimensional x_new: - Usage: - ``` - interpolant = AkimaInterpolant1DFlexible() + >>> import numpy as np + >>> from cudakima import AkimaInterpolant1DFlexible + >>> + >>> interpolant = AkimaInterpolant1DFlexible() + >>> + >>> # Standard usage (1D x_new broadcast to all groups) + >>> x = np.random.rand(10, 20) + >>> y = np.sin(x) + >>> x_new_1d = np.linspace(0, 1, 100) + >>> y_new = interpolant(x_new_1d, x, y) # Shape: (10, 100) + >>> + >>> # Multidimensional usage (different x_new for each group) + >>> x_new_multi = np.random.rand(10, 100) + >>> y_new = interpolant(x_new_multi, x, y) # Shape: (10, 100) - # Standard usage (1D x_new for all groups) - y_new = interpolant(x_new_1d, x, y) + Notes + ----- + The class automatically determines the mode based on x_new shape: + - If x_new.shape[:-1] == x.shape[:-1], uses multidimensional mode + - Otherwise, uses standard mode (broadcasts x_new) - # Multidimensional usage (different x_new for each group) - y_new = interpolant(x_new_multidim, x, y) - ``` + See Also + -------- + AkimaInterpolant1D : Base class for standard interpolation + AkimaInterpolant1DMultiDim : Multidimensional interpolation """ @property @@ -499,7 +747,30 @@ def multidim_interpolator(self): def __call__(self, x_new, x, y, **kwargs): """ - Automatically dispatch to appropriate interpolation method based on x_new shape. + Automatically dispatch to appropriate interpolation method. + + Determines whether to use standard or multidimensional interpolation + based on the shape of x_new. + + Parameters + ---------- + x_new : array_like + New x-values for interpolation. Can be either: + - 1D array (n_f,) - broadcast to all groups + - Multidimensional (..., n_f) - different for each group + x : array_like, shape (..., n) + X-values of the data points. + y : array_like, shape (..., n) + Y-values of the data points. + **kwargs : dict, optional + Additional keyword arguments. + + Returns + ------- + y_new : ndarray + Interpolated values. Shape depends on x_new: + - If x_new is 1D: shape (..., n_f) + - If x_new is multidimensional: shape (..., n_f) """ x_new = self.xp.asarray(x_new) x = self.xp.asarray(x) diff --git a/src/cudakima/kernels.py b/src/cudakima/kernels.py index f5f27fb..f373be4 100644 --- a/src/cudakima/kernels.py +++ b/src/cudakima/kernels.py @@ -22,6 +22,26 @@ @cuda.jit(device=True) def linearslope_gpu(x, y, idx): + """ + Compute the linear slope between two consecutive points on GPU. + + This device function calculates the slope (derivative) between points + at indices idx and idx+1. + + Parameters + ---------- + x : array_like + X-coordinates of data points + y : array_like + Y-coordinates of data points + idx : int + Index of the first point in the pair + + Returns + ------- + float + Linear slope (dy/dx) between points idx and idx+1 + """ dx = x[idx + 1] - x[idx] dy = y[idx + 1] - y[idx] m = dy / dx @@ -29,6 +49,36 @@ def linearslope_gpu(x, y, idx): @cuda.jit(device=True) def splineslope_gpu(x, y, idx, start, stop): + """ + Compute the Akima spline slope at a given point on GPU. + + This device function calculates the spline slope using the Akima algorithm, + which provides smooth interpolation with special handling for boundary conditions. + Requires at least 4 points due to boundary conditions. + + Parameters + ---------- + x : array_like + X-coordinates of data points + y : array_like + Y-coordinates of data points + idx : int + Index at which to compute the spline slope + start : int + Starting index of the data segment + stop : int + Ending index of the data segment (exclusive) + + Returns + ------- + float + Akima spline slope at the specified point + + Notes + ----- + The Akima spline uses weighted averages of neighboring slopes with special + boundary conditions at the endpoints and near-endpoints. + """ #! with these boundary conditions I ALWAYS NEED AT LEAST FOUR POINTS @@ -64,6 +114,30 @@ def splineslope_gpu(x, y, idx, start, stop): @cuda.jit def akima_spline_kernel_gpu(x_new, x, y, n_in, ngroups, nnans, result): + """ + GPU kernel for Akima spline interpolation. + + This kernel performs parallel Akima spline interpolation across multiple + groups of data. Falls back to linear interpolation for groups with less + than 4 valid points. + + Parameters + ---------- + x_new : array_like, shape (n_f,) + New x-coordinates at which to interpolate + x : array_like, shape (n_in * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (n_in * ngroups,) + Flattened array of y-coordinates for all groups + n_in : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group (including NaNs) + nnans : array_like, shape (n_in,) + Number of NaN values at the end of each group + result : array_like, shape (n_in * n_f,) + Output array for interpolated values + """ start1 = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x increment1 = cuda.blockDim.x * cuda.gridDim.x @@ -129,6 +203,26 @@ def akima_spline_kernel_gpu(x_new, x, y, n_in, ngroups, nnans, result): @numba.jit(nopython=True) def linearslope_cpu(x, y, idx): + """ + Compute the linear slope between two consecutive points on CPU. + + This function calculates the slope (derivative) between points + at indices idx and idx+1. + + Parameters + ---------- + x : array_like + X-coordinates of data points + y : array_like + Y-coordinates of data points + idx : int + Index of the first point in the pair + + Returns + ------- + float + Linear slope (dy/dx) between points idx and idx+1 + """ dx = x[idx + 1] - x[idx] dy = y[idx + 1] - y[idx] m = dy / dx @@ -136,6 +230,36 @@ def linearslope_cpu(x, y, idx): @numba.jit(nopython=True) def splineslope_cpu(x, y, idx, start, stop): + """ + Compute the Akima spline slope at a given point on CPU. + + This function calculates the spline slope using the Akima algorithm, + which provides smooth interpolation with special handling for boundary conditions. + Requires at least 4 points due to boundary conditions. + + Parameters + ---------- + x : array_like + X-coordinates of data points + y : array_like + Y-coordinates of data points + idx : int + Index at which to compute the spline slope + start : int + Starting index of the data segment + stop : int + Ending index of the data segment (exclusive) + + Returns + ------- + float + Akima spline slope at the specified point + + Notes + ----- + The Akima spline uses weighted averages of neighboring slopes with special + boundary conditions at the endpoints and near-endpoints. + """ #! with these boundary conditions I ALWAYS NEED AT LEAST FOUR POINTS @@ -171,6 +295,30 @@ def splineslope_cpu(x, y, idx, start, stop): @numba.jit(nopython=True) def akima_spline_kernel_cpu(x_new, x, y, n_in, ngroups, nnans, result): + """ + CPU kernel for Akima spline interpolation. + + This function performs Akima spline interpolation across multiple groups + of data. Falls back to linear interpolation for groups with less than 4 + valid points. + + Parameters + ---------- + x_new : array_like, shape (n_f,) + New x-coordinates at which to interpolate + x : array_like, shape (n_in * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (n_in * ngroups,) + Flattened array of y-coordinates for all groups + n_in : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group (including NaNs) + nnans : array_like, shape (n_in,) + Number of NaN values at the end of each group + result : array_like, shape (n_in * n_f,) + Output array for interpolated values + """ for j in range(n_in): # find the first and last non-NaN values @@ -226,7 +374,28 @@ def akima_spline_kernel_cpu(x_new, x, y, n_in, ngroups, nnans, result): @cuda.jit(device=True) def binary_search_gpu(x, target, start, stop): - """Binary search for interval location - much faster than linear search""" + """ + Binary search for interval location on GPU. + + This device function performs binary search to find the interval [x[i], x[i+1]) + that contains the target value. Much faster than linear search for large arrays. + + Parameters + ---------- + x : array_like + Sorted array of x-coordinates + target : float + Value to search for + start : int + Starting index of the search range + stop : int + Ending index of the search range (exclusive) + + Returns + ------- + int + Index i such that x[i] <= target < x[i+1] + """ left = start right = stop - 1 @@ -245,7 +414,27 @@ def binary_search_gpu(x, target, start, stop): @cuda.jit def precompute_slopes_kernel(x, y, slopes, nin, ngroups, nnans): - """Precompute all linear slopes in parallel""" + """ + GPU kernel to precompute all linear slopes in parallel. + + This kernel computes linear slopes between consecutive points for all groups, + storing them for later use in interpolation kernels. + + Parameters + ---------- + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + slopes : array_like, shape (nin * ngroups,) + Output array for computed linear slopes + nin : int + Number of groups + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + """ idx = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x group_idx = cuda.blockIdx.y @@ -264,7 +453,29 @@ def precompute_slopes_kernel(x, y, slopes, nin, ngroups, nnans): @cuda.jit def precompute_spline_slopes_kernel(x, y, linear_slopes, spline_slopes, nin, ngroups, nnans): - """Precompute spline slopes in parallel""" + """ + GPU kernel to precompute Akima spline slopes in parallel. + + This kernel computes Akima spline slopes at all points using precomputed + linear slopes, storing them for later use in interpolation kernels. + + Parameters + ---------- + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + spline_slopes : array_like, shape (nin * ngroups,) + Output array for computed Akima spline slopes + nin : int + Number of groups + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + """ idx = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x group_idx = cuda.blockIdx.y @@ -282,7 +493,32 @@ def precompute_spline_slopes_kernel(x, y, linear_slopes, spline_slopes, nin, ngr @cuda.jit(device=True) def splineslope_gpu_optimized(x, y, linear_slopes, idx, start, stop): - """Optimized spline slope computation using precomputed linear slopes""" + """ + Optimized Akima spline slope computation using precomputed linear slopes on GPU. + + This device function computes Akima spline slopes more efficiently by reusing + precomputed linear slopes instead of recalculating them. + + Parameters + ---------- + x : array_like + X-coordinates of data points + y : array_like + Y-coordinates of data points + linear_slopes : array_like + Precomputed linear slopes between consecutive points + idx : int + Index at which to compute the spline slope + start : int + Starting index of the data segment + stop : int + Ending index of the data segment (exclusive) + + Returns + ------- + float + Akima spline slope at the specified point + """ # Boundary conditions if idx == start: return (3 * linear_slopes[idx] - linear_slopes[idx + 1]) / 2 @@ -320,7 +556,34 @@ def splineslope_gpu_optimized(x, y, linear_slopes, idx, start, stop): @cuda.jit def akima_spline_kernel_optimized(x_new, x, y, linear_slopes, spline_slopes, nin, ngroups, nnans, result): - """Optimized main interpolation kernel""" + """ + Optimized GPU kernel for Akima spline interpolation with precomputed slopes. + + This kernel performs fast parallel Akima spline interpolation using + precomputed linear and spline slopes. Uses binary search for efficient + interval location. + + Parameters + ---------- + x_new : array_like, shape (n_f,) + New x-coordinates at which to interpolate + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + spline_slopes : array_like, shape (nin * ngroups,) + Precomputed Akima spline slopes at all points + nin : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + result : array_like, shape (nin * n_f,) + Output array for interpolated values + """ i = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x # x_new index j = cuda.blockIdx.y # group index @@ -361,6 +624,31 @@ def akima_spline_kernel_optimized(x_new, x, y, linear_slopes, spline_slopes, @cuda.jit def akima_linear_kernel(x_new, x, y, linear_slopes,nin, ngroups, nnans, result): + """ + GPU kernel for linear interpolation with precomputed slopes. + + This kernel performs fast parallel linear interpolation using precomputed + slopes and binary search for interval location. + + Parameters + ---------- + x_new : array_like, shape (n_f,) + New x-coordinates at which to interpolate + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + nin : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + result : array_like, shape (nin * n_f,) + Output array for interpolated values + """ i = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x # x_new index j = cuda.blockIdx.y # group index @@ -385,7 +673,28 @@ def akima_linear_kernel(x_new, x, y, linear_slopes,nin, ngroups, nnans, result): @numba.jit(nopython=True) def binary_search_cpu(x, target, start, stop): - """Binary search for interval location - CPU version""" + """ + Binary search for interval location on CPU. + + This function performs binary search to find the interval [x[i], x[i+1]) + that contains the target value. Much faster than linear search for large arrays. + + Parameters + ---------- + x : array_like + Sorted array of x-coordinates + target : float + Value to search for + start : int + Starting index of the search range + stop : int + Ending index of the search range (exclusive) + + Returns + ------- + int + Index i such that x[i] <= target < x[i+1] + """ left = start right = stop - 1 @@ -404,7 +713,30 @@ def binary_search_cpu(x, target, start, stop): @numba.jit(nopython=True) def precompute_all_linear_slopes_cpu(x, y, nin, ngroups, nnans): - """Precompute all linear slopes for all groups""" + """ + Precompute all linear slopes for all groups on CPU. + + This function computes linear slopes between consecutive points for all groups, + storing them for later use in interpolation functions. + + Parameters + ---------- + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + nin : int + Number of groups + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + + Returns + ------- + ndarray, shape (nin * ngroups,) + Array of computed linear slopes + """ total_points = nin * ngroups linear_slopes = np.zeros(total_points) @@ -421,7 +753,32 @@ def precompute_all_linear_slopes_cpu(x, y, nin, ngroups, nnans): @numba.jit(nopython=True) def precompute_all_spline_slopes_cpu(x, y, linear_slopes, nin, ngroups, nnans): - """Precompute all spline slopes for all groups""" + """ + Precompute all Akima spline slopes for all groups on CPU. + + This function computes Akima spline slopes at all points using precomputed + linear slopes, storing them for later use in interpolation functions. + + Parameters + ---------- + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + nin : int + Number of groups + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + + Returns + ------- + ndarray, shape (nin * ngroups,) + Array of computed Akima spline slopes + """ total_points = nin * ngroups spline_slopes = np.zeros(total_points) @@ -439,7 +796,32 @@ def precompute_all_spline_slopes_cpu(x, y, linear_slopes, nin, ngroups, nnans): @numba.jit(nopython=True) def splineslope_cpu_optimized(x, y, linear_slopes, idx, start, stop): - """Optimized spline slope computation using precomputed linear slopes""" + """ + Optimized Akima spline slope computation using precomputed linear slopes on CPU. + + This function computes Akima spline slopes more efficiently by reusing + precomputed linear slopes instead of recalculating them. + + Parameters + ---------- + x : array_like + X-coordinates of data points + y : array_like + Y-coordinates of data points + linear_slopes : array_like + Precomputed linear slopes between consecutive points + idx : int + Index at which to compute the spline slope + start : int + Starting index of the data segment + stop : int + Ending index of the data segment (exclusive) + + Returns + ------- + float + Akima spline slope at the specified point + """ # Boundary conditions if idx == start: return (3 * linear_slopes[idx] - linear_slopes[idx + 1]) / 2 @@ -477,7 +859,34 @@ def splineslope_cpu_optimized(x, y, linear_slopes, idx, start, stop): @numba.jit(nopython=True, parallel=True) def akima_spline_kernel_cpu_optimized(x_new, x, y, linear_slopes, spline_slopes, nin, ngroups, nnans, result): - """Optimized CPU kernel with precomputed slopes and parallel execution""" + """ + Optimized CPU kernel for Akima spline interpolation with precomputed slopes. + + This function performs fast parallel Akima spline interpolation using + precomputed linear and spline slopes. Uses binary search for efficient + interval location and numba parallel execution. + + Parameters + ---------- + x_new : array_like, shape (n_f,) + New x-coordinates at which to interpolate + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + spline_slopes : array_like, shape (nin * ngroups,) + Precomputed Akima spline slopes at all points + nin : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + result : array_like, shape (nin * n_f,) + Output array for interpolated values + """ # Parallel loop over groups for j in numba.prange(nin): @@ -519,6 +928,31 @@ def akima_spline_kernel_cpu_optimized(x_new, x, y, linear_slopes, spline_slopes, @numba.jit(nopython=True, parallel=True) def akima_linear_kernel_cpu(x_new, x, y, linear_slopes,nin, ngroups, nnans, result): + """ + CPU kernel for linear interpolation with precomputed slopes. + + This function performs fast parallel linear interpolation using precomputed + slopes and binary search for interval location. + + Parameters + ---------- + x_new : array_like, shape (n_f,) + New x-coordinates at which to interpolate + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + nin : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + result : array_like, shape (nin * n_f,) + Output array for interpolated values + """ for j in range(nin): start = j * ngroups stop = (j + 1) * ngroups - nnans[j] @@ -544,7 +978,31 @@ def akima_linear_kernel_cpu(x_new, x, y, linear_slopes,nin, ngroups, nnans, resu @cuda.jit def akima_linear_kernel_gpu_multidim(x_new, x, y, linear_slopes, nin, ngroups, nnans, result): - """GPU kernel for multidimensional x_new arrays""" + """ + GPU kernel for linear interpolation with multidimensional x_new arrays. + + This kernel handles cases where each group has its own set of interpolation + points (x_new is 2D with shape (nin, n_f)). + + Parameters + ---------- + x_new : array_like, shape (nin, n_f) + New x-coordinates at which to interpolate, different for each group + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + nin : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + result : array_like, shape (nin * n_f,) + Output array for interpolated values + """ i = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x # x_new index within group j = cuda.blockIdx.y # group index @@ -579,7 +1037,31 @@ def akima_linear_kernel_gpu_multidim(x_new, x, y, linear_slopes, @numba.jit(nopython=True, parallel=True) def akima_linear_kernel_cpu_multidim(x_new, x, y, linear_slopes, nin, ngroups, nnans, result): - """CPU kernel for multidimensional x_new arrays""" + """ + CPU kernel for linear interpolation with multidimensional x_new arrays. + + This function handles cases where each group has its own set of interpolation + points (x_new is 2D with shape (nin, n_f)). + + Parameters + ---------- + x_new : array_like, shape (nin, n_f) + New x-coordinates at which to interpolate, different for each group + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + nin : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + result : array_like, shape (nin * n_f,) + Output array for interpolated values + """ # Parallel loop over groups for j in numba.prange(nin): @@ -608,7 +1090,33 @@ def akima_linear_kernel_cpu_multidim(x_new, x, y, linear_slopes, @cuda.jit def akima_spline_kernel_gpu_multidim(x_new, x, y, linear_slopes, spline_slopes, nin, ngroups, nnans, result): - """GPU kernel for multidimensional x_new arrays""" + """ + GPU kernel for Akima spline interpolation with multidimensional x_new arrays. + + This kernel handles cases where each group has its own set of interpolation + points (x_new is 2D with shape (nin, n_f)). + + Parameters + ---------- + x_new : array_like, shape (nin, n_f) + New x-coordinates at which to interpolate, different for each group + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + spline_slopes : array_like, shape (nin * ngroups,) + Precomputed Akima spline slopes at all points + nin : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + result : array_like, shape (nin * n_f,) + Output array for interpolated values + """ i = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x # x_new index within group j = cuda.blockIdx.y # group index @@ -658,7 +1166,33 @@ def akima_spline_kernel_gpu_multidim(x_new, x, y, linear_slopes, spline_slopes, @numba.jit(nopython=True, parallel=True) def akima_spline_kernel_cpu_multidim(x_new, x, y, linear_slopes, spline_slopes, nin, ngroups, nnans, result): - """CPU kernel for multidimensional x_new arrays""" + """ + CPU kernel for Akima spline interpolation with multidimensional x_new arrays. + + This function handles cases where each group has its own set of interpolation + points (x_new is 2D with shape (nin, n_f)). + + Parameters + ---------- + x_new : array_like, shape (nin, n_f) + New x-coordinates at which to interpolate, different for each group + x : array_like, shape (nin * ngroups,) + Flattened array of x-coordinates for all groups + y : array_like, shape (nin * ngroups,) + Flattened array of y-coordinates for all groups + linear_slopes : array_like, shape (nin * ngroups,) + Precomputed linear slopes between consecutive points + spline_slopes : array_like, shape (nin * ngroups,) + Precomputed Akima spline slopes at all points + nin : int + Number of groups to interpolate + ngroups : int + Maximum number of points per group + nnans : array_like, shape (nin,) + Number of NaN values at the end of each group + result : array_like, shape (nin * n_f,) + Output array for interpolated values + """ # Parallel loop over groups for j in numba.prange(nin):