diff --git a/docs-samples/data-science/gpu-accelerated-samples/README.md b/docs-samples/data-science/gpu-accelerated-samples/README.md new file mode 100644 index 00000000..2f15c9b7 --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/README.md @@ -0,0 +1,50 @@ +# GPU-Accelerated Samples + +Jupyter notebooks demonstrating GPU acceleration using NVIDIA RAPIDS libraries on Microsoft Fabric and Azure VMs with NVIDIA GPUs. + +> Based on examples from [NVIDIA RAPIDS](https://github.com/rapidsai) (Apache 2.0 License, Copyright NVIDIA Corporation). Modified and tested for Azure GPU VM environment. + +## Requirements + +- NVIDIA GPU (Tesla T4 or better) +- CUDA 12.x +- Python 3.12+ +- RAPIDS 25.x (cuDF, cuML, cuCIM) +- PyOD (`pip install pyod`) for anomaly detection demo +- Conda environment recommended + +## Notebooks + +| Notebook | Description | +|----------|-------------| +| [rapids-gpu-accelerated-demo.ipynb](rapids-gpu-accelerated-demo.ipynb) | End-to-end GPU vs CPU benchmark: DataFrames, strings, KMeans, Random Forest, text embeddings, and KNN search | +| [cudf-pandas-stock-analysis.ipynb](cudf-pandas-stock-analysis.ipynb) | Stock market data analysis using GPU-accelerated pandas (read, merge, resample, plot) | +| [cuml-scikit-learn-accelerator-demo.ipynb](cuml-scikit-learn-accelerator-demo.ipynb) | Drop-in GPU acceleration for scikit-learn (PCA, UMAP, KNN, HDBSCAN on activity recognition data) | +| [pyod-gpu-anomaly-detection-demo.ipynb](pyod-gpu-anomaly-detection-demo.ipynb) | GPU-accelerated anomaly detection with PyOD and cuml.accel (IForest, PCA, KNN, CBLOF, HDBSCAN) | + +### Notes + +- **cudf-pandas-stock-analysis.ipynb** uses `%load_ext cudf.pandas` which must be the very first executed statement before any `import pandas`. In Fabric, pandas may be pre-imported by the environment — if so, add `cudf.pandas` to the Fabric session pre-run script. +- **cuml-scikit-learn-accelerator-demo.ipynb** downloads the UCI HAR dataset to `/tmp/HAR_data/` on first run (requires internet access). +- **cudf-pandas-stock-analysis.ipynb** downloads stock price data from Yahoo Finance (requires internet access). +- **pyod-gpu-anomaly-detection-demo.ipynb** uses `cuml.accel` for GPU acceleration. Dataset is synthetic (generated in-notebook). Increase `N_TRAIN` to 100K+ for meaningful GPU speedups. + +## cuCIM Medical Imaging + +GPU-accelerated image processing for digital pathology and microscopy using [cuCIM](https://github.com/rapidsai/cucim). + +> Based on examples from [rapidsai/cucim](https://github.com/rapidsai/cucim) (Apache 2.0 License, Copyright NVIDIA Corporation). Modified for Azure GPU VM environment. + +| Notebook | Description | +|----------|-------------| +| [01-whole-slide-image-reading.ipynb](cucim_medical_imaging/01-whole-slide-image-reading.ipynb) | Reading and displaying multi-resolution whole-slide pathology images | +| [02-image-cache-performance.ipynb](cucim_medical_imaging/02-image-cache-performance.ipynb) | Tile caching strategies for faster repeated region access | +| [03-gabor-texture-classification.ipynb](cucim_medical_imaging/03-gabor-texture-classification.ipynb) | GPU vs CPU Gabor filter bank for texture classification | +| [04-random-walker-segmentation.ipynb](cucim_medical_imaging/04-random-walker-segmentation.ipynb) | GPU-accelerated random walker image segmentation | +| [05-vesselness-filter.ipynb](cucim_medical_imaging/05-vesselness-filter.ipynb) | GPU-accelerated vessel detection using Frangi/Sato filters | + +### Setup for cuCIM notebooks + +**Notebooks 01 and 02** auto-download a sample whole-slide image (~170MB from OpenSlide) on first run. No manual setup needed. + +**Notebooks 03, 04, and 05** are self-contained — they use synthetic data from scikit-image and do not require any additional files. diff --git a/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/01-whole-slide-image-reading.ipynb b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/01-whole-slide-image-reading.ipynb new file mode 100644 index 00000000..2a16a1b6 --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/01-whole-slide-image-reading.ipynb @@ -0,0 +1,337 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "eb420212", + "metadata": {}, + "source": [ + "# Basic Usage\n", + "\n", + "This notebook shows basic usage of cuCIM with a whole-slide pathology image." + ] + }, + { + "cell_type": "markdown", + "id": "39389384", + "metadata": {}, + "source": [ + "## Prerequisites" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "93157b57", + "metadata": {}, + "outputs": [], + "source": [ + "# !conda install -c conda-forge pillow" + ] + }, + { + "cell_type": "markdown", + "id": "6378c4a3", + "metadata": {}, + "source": [ + "## Read image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "auto-download", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import urllib.request\n", + "\n", + "os.makedirs(\"input\", exist_ok=True)\n", + "\n", + "IMAGE_PATH = \"input/image.tif\"\n", + "IMAGE_URL = \"https://openslide.cs.cmu.edu/download/openslide-testdata/Aperio/CMU-1.svs\"\n", + "\n", + "if not os.path.exists(IMAGE_PATH):\n", + " print(f\"Downloading sample whole-slide image (~170MB)...\")\n", + " urllib.request.urlretrieve(IMAGE_URL, IMAGE_PATH)\n", + " print(f\"Downloaded to {IMAGE_PATH}\")\n", + "else:\n", + " print(f\"Image already exists: {IMAGE_PATH}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c02bf93e", + "metadata": {}, + "outputs": [], + "source": [ + "from cucim import CuImage\n", + "\n", + "img = CuImage(\"input/image.tif\")\n", + "print(f\"Image loaded: {img.path}\")\n", + "print(f\"Dimensions (XY): {img.size('XY')}\")" + ] + }, + { + "cell_type": "markdown", + "id": "2b371170", + "metadata": {}, + "source": [ + "### See metadata" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ff52186", + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "\n", + "metadata = img.metadata\n", + "print(json.dumps(metadata, indent=2)[:2000])" + ] + }, + { + "cell_type": "markdown", + "id": "5ffeeeed", + "metadata": {}, + "source": [ + "### Read region\n", + "\n", + "Let's read the whole slide at the lowest resolution and save as a thumbnail." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec613ae8", + "metadata": {}, + "outputs": [], + "source": [ + "from PIL import Image\n", + "\n", + "resolutions = img.resolutions\n", + "level_dimensions = resolutions[\"level_dimensions\"]\n", + "level_count = resolutions[\"level_count\"]\n", + "\n", + "print(f\"Number of levels: {level_count}\")\n", + "print(f\"Level dimensions: {level_dimensions}\")\n", + "\n", + "# Read the whole slide at lowest resolution\n", + "region = img.read_region(\n", + " location=[0, 0], size=level_dimensions[level_count - 1], level=level_count - 1\n", + ")\n", + "\n", + "# Save and display\n", + "region.save(\"input/thumbnail.ppm\")\n", + "Image.open(\"input/thumbnail.ppm\")" + ] + }, + { + "cell_type": "markdown", + "id": "2f6da447", + "metadata": {}, + "source": [ + "Now let's read a 512x512 region at full resolution (level 0)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "daacf8cf", + "metadata": {}, + "outputs": [], + "source": [ + "from PIL import Image\n", + "\n", + "region = img.read_region([10000, 20000], [512, 512], 0)\n", + "region.save(\"input/test.ppm\")\n", + "Image.open(\"input/test.ppm\")" + ] + }, + { + "cell_type": "markdown", + "id": "6c367337", + "metadata": {}, + "source": [ + "### `__array_interface__` support\n", + "\n", + "A NumPy array has the `__array_interface__` property. Here is a simple example:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e09b53d6", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "np_arr = np.array([1, 2, 3])\n", + "print(np_arr.__array_interface__)" + ] + }, + { + "cell_type": "markdown", + "id": "f7a3f798", + "metadata": {}, + "source": [ + "As you can see from the above result, a NumPy array has `__array_interface__` property which describes the memory layout.\n", + "\n", + "**cuCIM also supports `__array_interface__`**, so you can directly convert a cuCIM region to a NumPy array and display it with PIL:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0804834e", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from PIL import Image\n", + "\n", + "region = img.read_region([10000, 20000], [512, 512], 0)\n", + "np_img_arr = np.asarray(region)\n", + "\n", + "print(f\"Array shape: {np_img_arr.shape}, dtype: {np_img_arr.dtype}\")\n", + "arr_iface = np_img_arr.__array_interface__\n", + "print(f\"__array_interface__ shape: {arr_iface['shape']}\")\n", + "\n", + "Image.fromarray(np_img_arr)" + ] + }, + { + "cell_type": "markdown", + "id": "0b702b4e", + "metadata": {}, + "source": [ + "### `__cuda_array_interface__` support\n", + "\n", + "A CuPy array has `__cuda_array_interface__` property. Here is a simple example:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0990d4d4", + "metadata": {}, + "outputs": [], + "source": [ + "import cupy as cp\n", + "\n", + "cp_arr = cp.array([1, 2, 3])\n", + "print(cp_arr.__cuda_array_interface__)" + ] + }, + { + "cell_type": "markdown", + "id": "a130d04f", + "metadata": {}, + "source": [ + "A Python object that has `__cuda_array_interface__` property is considered as a CUDA array-like object and can be converted to CuPy array through `cp.asarray(obj)` method.\n", + "\n", + "**`__cuda_array_interface__` is also supported in CuImage** if you specify `device='cuda'` in `read_region()`. The following code loads image data directly to GPU memory and visualizes it in the Jupyter notebook:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5a43d245", + "metadata": {}, + "outputs": [], + "source": [ + "import cupy as cp\n", + "from cucim import CuImage\n", + "from PIL import Image\n", + "\n", + "img = CuImage(\"input/image.tif\")\n", + "resolutions = img.resolutions\n", + "level_dimensions = resolutions[\"level_dimensions\"]\n", + "level_count = resolutions[\"level_count\"]\n", + "\n", + "# Read region directly to GPU memory\n", + "region_gpu = img.read_region(\n", + " [0, 0], level_dimensions[level_count - 1], level_count - 1, device=\"cuda\"\n", + ")\n", + "\n", + "print(f\"Device: {region_gpu.device}\")\n", + "cuda_iface = region_gpu.__cuda_array_interface__\n", + "print(f\"__cuda_array_interface__ shape: {cuda_iface['shape']}\")\n", + "\n", + "# Convert to CuPy array, then to CPU for display\n", + "cupy_arr = cp.asarray(region_gpu)\n", + "print(f\"CuPy array shape: {cupy_arr.shape}\")\n", + "Image.fromarray(cupy_arr.get())" + ] + }, + { + "cell_type": "markdown", + "id": "499f2e79", + "metadata": {}, + "source": [ + "### Associated images\n", + "\n", + "Some image formats such as Philips TIFF and Aperio SVS have associated images (Macro or Label images) in addition to the multi-resolution images.\n", + "\n", + "Let's check what associated images are available:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9714d55", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from PIL import Image\n", + "\n", + "print(f\"Associated images: {img.associated_images}\")\n", + "\n", + "for name in img.associated_images:\n", + " assoc_img = img.associated_image(name)\n", + " np_arr = np.asarray(assoc_img)\n", + " print(f\" '{name}': shape={np_arr.shape}\")\n", + "\n", + "# Display the macro image if available\n", + "if 'macro' in img.associated_images:\n", + " macro_image = img.associated_image('macro')\n", + " np_img_arr = np.asarray(macro_image)\n", + " print(f\"\\nDisplaying macro image: {np_img_arr.shape}\")\n", + " display(Image.fromarray(np_img_arr))\n", + "\n", + "# Check for label image\n", + "if 'label' in img.associated_images:\n", + " print(\"Label image exists\")\n", + "else:\n", + " print(\"\\nThere is no associated image named 'label'!\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/02-image-cache-performance.ipynb b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/02-image-cache-performance.ipynb new file mode 100644 index 00000000..0fdb6dbc --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/02-image-cache-performance.ipynb @@ -0,0 +1,546 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "bfef84ba", + "metadata": {}, + "source": [ + "# Using Cache (available since v21.06.00)\n", + "\n", + "## Need for Cache\n", + "\n", + "In many deep learning use cases, small image patches need to be extracted from the large image and they are fed into the neural network. \n", + "\n", + "If the patch size doesn't align with the underlying tile layout of TIFF image (e.g., AI model such as ResNet may accept a particular size of the image [e.g., 224x224] that is smaller than the underlying tile size [256x256]), redundant image loadings for a tile are needed (See the following two figures)\n", + "\n", + "![image](https://user-images.githubusercontent.com/1928522/118344267-333a4f00-b4e2-11eb-898c-8980c8725d32.png)\n", + "![image](https://user-images.githubusercontent.com/1928522/118344294-5238e100-b4e2-11eb-8f3a-4772ef055658.png)\n", + "\n", + "Which resulted in lower performance for unaligned cases as shown in our [GTC 2021 presentation](https://www.nvidia.com/en-us/gtc/catalog/?search=cuCIM)\n", + "\n", + "![image](https://user-images.githubusercontent.com/1928522/118344737-c07ea300-b4e4-11eb-9c95-15c2e5022274.png)\n", + "\n", + "\n", + "The proper use of cache improves the loading performance greatly, especially for **inference** use cases and when [accessing tiles sequentially (left to right, top to bottom) from one TIFF file](https://nbviewer.jupyter.org/github/rapidsai/cucim/blob/branch-21.06/notebooks/File-access_Experiments_on_TIFF.ipynb#1.-Accessing-tiles-sequentially-(left-to-right,-top-to-bottom)-from-one-TIFF-file).\n", + "\n", + "On the other hand, if the application [accesses partial tiles randomly from multiple TIFF files](https://nbviewer.jupyter.org/github/rapidsai/cucim/blob/branch-21.06/notebooks/File-access_Experiments_on_TIFF.ipynb#3.-Accessing-partial-tiles-randomly-from-multiple-TIFF-files) (this usually happens for **training** use cases), using a cache could be meaningless." + ] + }, + { + "cell_type": "markdown", + "id": "e952a222", + "metadata": {}, + "source": [ + "## Enabling cache\n", + "\n", + "Currently, cuCIM supports the following three strategies:\n", + "\n", + " - `nocache`\n", + " - `per_process`\n", + " - `shared_memory` (interprocess)\n", + "\n", + "\n", + "**1) `nocache`**\n", + "\n", + "No cache.\n", + "\n", + "By default, this cache strategy is used.\n", + "With this strategy, the behavior is the same as one before `v20.06.00`.\n", + "\n", + "**2) `per_process`**\n", + "\n", + "The cache memory is shared among threads.\n", + "\n", + "**3) `shared_memory`**\n", + "\n", + "The cache memory is shared among processes.\n", + "\n", + "### Getting cache setting\n", + "\n", + "`CuImage.cache()` would return an object that can control the current cache. The object has the following properties:\n", + "\n", + "- `type`: The type (strategy) name\n", + "- `memory_size`: The number of bytes used in the cache memory\n", + "- `memory_capacity`: The maximum number of bytes that can be allocated (used) in the cache memory\n", + "- `free_memory`: The number of bytes available in the cache memory\n", + "- `size`: The number of cache items used\n", + "- `capacity`: The maximum number of cache items that can be created\n", + "- `hit_count`: The cache hit count\n", + "- `miss_count`: The cache miss count\n", + "- `config`: A configuration dictionary that was used for configuring cache.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1845732", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "if not os.getcwd().endswith(\"cucim_medical_imaging\"):\n", + " os.chdir(os.path.join(os.getcwd(), \"cucim_medical_imaging\"))\n", + "print(f\"Working directory: {os.getcwd()}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac9aa319", + "metadata": {}, + "outputs": [], + "source": [ + "from cucim import CuImage\n", + "\n", + "cache = CuImage.cache()\n", + "\n", + "print(f\" type: {cache.type}({int(cache.type)})\")\n", + "print(f\"memory_size: {cache.memory_size}/{cache.memory_capacity}\")\n", + "print(f\"free_memory: {cache.free_memory}\")\n", + "print(f\" size: {cache.size}/{cache.capacity}\")\n", + "print(f\" hit_count: {cache.hit_count}\")\n", + "print(f\" miss_count: {cache.miss_count}\")\n", + "print(f\" config: {cache.config}\")" + ] + }, + { + "cell_type": "markdown", + "id": "f057a11a", + "metadata": {}, + "source": [ + "### Changing Cache Setting\n", + "\n", + "Cache configuration can be changed by adding parameters to `cache()` method.\n", + "\n", + "The following parameters are available:\n", + "\n", + "- `type`: The type (strategy) name. Default to 'no_cache'.\n", + "- `memory_capacity`: The maximum number of mebibytes (`MiB`, 2^20) that can be allocated (used) in the cache memory. Default to `1024`.\n", + "- `capacity`: The maximum number of cache items that can be created. Default to `5461` (= (\\ x 2^20) / (256x256x3)).\n", + "- `mutex_pool_capacity`: The mutex pool size. Default to `11117`.\n", + "- `list_padding`: The number of additional items used for the internal circular queue. Default to `10000`.\n", + "- `extra_shared_memory_size`: The size of additional memory allocation (in MiB) for shared_memory allocator in `shared_process` strategy. Default to `100`.\n", + "- `record_stat`: If the cache statistic should be recorded or not. Default to `False`.\n", + "\n", + "In most cases, `type`(required) and `memory_capacity` are used." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7d3090d", + "metadata": {}, + "outputs": [], + "source": [ + "from cucim import CuImage\n", + "\n", + "cache = CuImage.cache(\"per_process\", memory_capacity=2048)\n", + "print(f\" type: {cache.type}({int(cache.type)})\")\n", + "print(f\"memory_size: {cache.memory_size}/{cache.memory_capacity}\")\n", + "print(f\"free_memory: {cache.free_memory}\")\n", + "print(f\" size: {cache.size}/{cache.capacity}\")\n", + "print(f\" hit_count: {cache.hit_count}\")\n", + "print(f\" miss_count: {cache.miss_count}\")\n", + "print(f\" config: {cache.config}\")" + ] + }, + { + "cell_type": "markdown", + "id": "fbd45a9e", + "metadata": {}, + "source": [ + "## Choosing Proper Cache Memory Size\n", + "\n", + "It is important to select the appropriate cache memory size (capacity). Small cache memory size results in low cache hit rates. Conversely, if the cache memory size is too large, memory is wasted.\n", + "\n", + "For example, if the default tile size is 256x256 and the patch size to load is 224x224, the cache memory needs to be large enough to contain at least two rows of tiles in the image to avoid deleting the required cache entries while loading patches sequentially (left to right, top to bottom) from one TIFF file.\n", + "\n", + "![image](https://user-images.githubusercontent.com/1928522/120760720-4cbf2d00-c4c9-11eb-875b-b070203fd8e6.png)\n", + "\n", + "cuCIM provide a utility method (`cucim.clara.cache.preferred_memory_capacity()`) to calculate a preferred cache memory size for the given image (image size and tile size) and the patch size.\n", + "\n", + "Internal logic is available at \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bfb70aa4", + "metadata": {}, + "outputs": [], + "source": [ + "from cucim import CuImage\n", + "from cucim.clara.cache import preferred_memory_capacity\n", + "\n", + "img = CuImage(\"input/image.tif\")\n", + "\n", + "image_size = img.size(\"XY\") # same with `img.resolutions[\"level_dimensions\"][0]`\n", + "tile_size = img.resolutions[\"level_tile_sizes\"][0] # default: (256, 256)\n", + "patch_size = (1024, 1024) # default: (256, 256)\n", + "bytes_per_pixel = 3 # default: 3\n", + "\n", + "print(f\"image size: {image_size}\")\n", + "print(f\"tile size: {tile_size}\")\n", + "\n", + "# Below three statements are the same.\n", + "memory_capacity = preferred_memory_capacity(img, patch_size=patch_size)\n", + "memory_capacity2 = preferred_memory_capacity(\n", + " None, image_size, tile_size, patch_size, bytes_per_pixel\n", + ")\n", + "memory_capacity3 = preferred_memory_capacity(None, image_size, patch_size=patch_size)\n", + "\n", + "print(f\"memory_capacity : {memory_capacity} MiB\")\n", + "print(f\"memory_capacity2: {memory_capacity2} MiB\")\n", + "print(f\"memory_capacity3: {memory_capacity3} MiB\")\n", + "\n", + "cache = CuImage.cache(\n", + " \"per_process\", memory_capacity=memory_capacity\n", + ") # You can also manually set capacity` (e.g., `capacity=500`)\n", + "print(\"= Cache Info =\")\n", + "print(f\" type: {cache.type}({int(cache.type)})\")\n", + "print(f\"memory_size: {cache.memory_size}/{cache.memory_capacity}\")\n", + "print(f\" size: {cache.size}/{cache.capacity}\")" + ] + }, + { + "cell_type": "markdown", + "id": "5898b317", + "metadata": {}, + "source": [ + "### Reserve More Cache Memory\n", + "\n", + "If more cache memory capacity is needed in runtime, you can use `reserve()` method.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61801fe7", + "metadata": {}, + "outputs": [], + "source": [ + "from cucim import CuImage\n", + "from cucim.clara.cache import preferred_memory_capacity\n", + "\n", + "img = CuImage(\"input/image.tif\")\n", + "\n", + "memory_capacity = preferred_memory_capacity(img, patch_size=(256, 256))\n", + "new_memory_capacity = preferred_memory_capacity(img, patch_size=(512, 512))\n", + "\n", + "print(f\"memory_capacity : {memory_capacity} MiB\")\n", + "print(f\"new_memory_capacity: {new_memory_capacity} MiB\")\n", + "print()\n", + "\n", + "cache = CuImage.cache(\"per_process\", memory_capacity=memory_capacity)\n", + "print(\"= Cache Info =\")\n", + "print(f\" type: {cache.type}({int(cache.type)})\")\n", + "print(f\"memory_size: {cache.memory_size}/{cache.memory_capacity}\")\n", + "print(f\" size: {cache.size}/{cache.capacity}\")\n", + "print()\n", + "\n", + "cache.reserve(new_memory_capacity)\n", + "print(\"= Cache Info (update memory capacity) =\")\n", + "print(f\" type: {cache.type}({int(cache.type)})\")\n", + "print(f\"memory_size: {cache.memory_size}/{cache.memory_capacity}\")\n", + "print(f\" size: {cache.size}/{cache.capacity}\")\n", + "print()\n", + "\n", + "cache.reserve(memory_capacity, capacity=500)\n", + "print(\"= Cache Info (update memory capacity & capacity) =\")\n", + "print(f\" type: {cache.type}({int(cache.type)})\")\n", + "print(\n", + " f\"memory_size: {cache.memory_size}/{cache.memory_capacity} # smaller `memory_capacity` value does not change this\"\n", + ")\n", + "print(f\" size: {cache.size}/{cache.capacity}\")\n", + "print()\n", + "\n", + "cache = CuImage.cache(\"no_cache\")\n", + "print(\"= Cache Info (no cache) =\")\n", + "print(f\" type: {cache.type}({int(cache.type)})\")\n", + "print(f\"memory_size: {cache.memory_size}/{cache.memory_capacity}\")\n", + "print(f\" size: {cache.size}/{cache.capacity}\")" + ] + }, + { + "cell_type": "markdown", + "id": "f8ffe4cb", + "metadata": {}, + "source": [ + "## Profiling Cache Hit/Miss\n", + "\n", + "If you add an argument `record_stat=True` to `CuImage.cache()` method, cache statistics is recorded.\n", + "\n", + "Cache hit/miss count is accessible through `hit_count`/`miss_count` property of the cache object.\n", + "\n", + "You can get/set/unset the recording through `record()` method.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a75f92b", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91587c98", + "metadata": {}, + "outputs": [], + "source": [ + "from cucim import CuImage\n", + "from cucim.clara.cache import preferred_memory_capacity\n", + "\n", + "img = CuImage(\"input/image.tif\")\n", + "memory_capacity = preferred_memory_capacity(img, patch_size=(256, 256))\n", + "cache = CuImage.cache(\"per_process\", memory_capacity=memory_capacity, record_stat=True)\n", + "\n", + "img.read_region((0, 0), (100, 100))\n", + "print(f\"cache hit: {cache.hit_count}, cache miss: {cache.miss_count}\")\n", + "\n", + "region = img.read_region((0, 0), (100, 100))\n", + "print(f\"cache hit: {cache.hit_count}, cache miss: {cache.miss_count}\")\n", + "\n", + "region = img.read_region((0, 0), (100, 100))\n", + "print(f\"cache hit: {cache.hit_count}, cache miss: {cache.miss_count}\")\n", + "\n", + "print(f\"Is recorded: {cache.record()}\")\n", + "\n", + "cache.record(False)\n", + "print(f\"Is recorded: {cache.record()}\")\n", + "\n", + "region = img.read_region((0, 0), (100, 100))\n", + "print(f\"cache hit: {cache.hit_count}, cache miss: {cache.miss_count}\")\n", + "print()\n", + "\n", + "print(f\" type: {cache.type}({int(cache.type)})\")\n", + "print(f\"memory_size: {cache.memory_size}/{cache.memory_capacity}\")\n", + "print(f\"free_memory: {cache.free_memory}\")\n", + "print(f\" size: {cache.size}/{cache.capacity}\")\n", + "print()\n", + "\n", + "cache = CuImage.cache(\"no_cache\")\n", + "print(f\" type: {cache.type}({int(cache.type)})\")\n", + "print(f\"memory_size: {cache.memory_size}/{cache.memory_capacity}\")\n", + "print(f\"free_memory: {cache.free_memory}\")\n", + "print(f\" size: {cache.size}/{cache.capacity}\")" + ] + }, + { + "cell_type": "markdown", + "id": "d4b4b55e", + "metadata": {}, + "source": [ + "## Considerations in Multi-threading/processing Environment\n", + "\n", + "\n", + "### `per_process` strategy\n", + "\n", + "#### Cache memory\n", + "\n", + "If used in the multi-threading environment and each thread is reading the different part of the image sequentially, please consider increasing cache memory size than the size suggested by `cucim.clara.cache.preferred_memory_capacity()` to avoid dropping necessary cache items.\n", + "\n", + "If used in the multi-processing environment, the cache memory size allocated can be `(# of processes) x (cache memory capacity)`. \n", + "\n", + "Please be careful not to oversize the memory allocated by the cache.\n", + "\n", + "\n", + "#### Cache Statistics\n", + "\n", + "If used in the multi-processing environment (e.g, using `concurrent.futures.ProcessPoolExecutor()`), cache hit count (`hit_count`) and miss count (`miss_count`) wouldn't be recorded in the main process's cache object.\n", + "\n", + "\n", + "### `shared_memory` strategy\n", + "\n", + "In general, `shared_memory` strategy has more overhead than `per_process` strategy. However, it is recommended that you select this strategy if you want to use a fixed size of cache memory regardless of the number of processes.\n", + "\n", + "Note that, this strategy pre-allocates the cache memory in the shared memory and allocates more memory (as specified in `extra_shared_memory_size` parameter) than the requested cache memory size (capacity) for the memory allocator to handle memory segments.\n", + "\n", + "\n", + "#### Cache memory\n", + "\n", + "Since the cache memory would be shared by multiple threads/processes, you will need to set enough cache memory to avoid dropping necessary cache items.\n" + ] + }, + { + "cell_type": "markdown", + "id": "f5b247db", + "metadata": {}, + "source": [ + "## Setting Default Cache Configuration\n", + "\n", + "The configuration for cuCIM can be specified in `.cucim.json` file and user can set a default cache settings there.\n", + "\n", + "cuCIM finds `.cucim.json` file from the following order:\n", + "\n", + "1. The current folder\n", + "2. `$HOME/.cucim.json`\n", + "\n", + "The configuration for the cache can be specified like below.\n", + "\n", + "```jsonc\n", + "\n", + "{\n", + " // This is actually JSONC file so comments are available.\n", + " \"cache\": {\n", + " \"type\": \"nocache\",\n", + " \"memory_capacity\": 1024,\n", + " \"capacity\": 5461,\n", + " \"mutex_pool_capacity\": 11117,\n", + " \"list_padding\": 10000,\n", + " \"extra_shared_memory_size\": 100,\n", + " \"record_stat\": false\n", + " }\n", + "}\n", + "```\n", + "\n", + "You can write the current cache configuration into the file like below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1dda19e6", + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "from cucim import CuImage\n", + "\n", + "cache = CuImage.cache()\n", + "config_data = {\"cache\": cache.config}\n", + "json_text = json.dumps(config_data, indent=4)\n", + "print(json_text)\n", + "\n", + "# Save into the configuration file.\n", + "with open(\".cucim.json\", \"w\") as fp:\n", + " fp.write(json_text)" + ] + }, + { + "cell_type": "markdown", + "id": "60c2d934", + "metadata": {}, + "source": [ + "### Cache Mechanism Used in Other Libraries (OpenSlide and rasterio)\n", + "\n", + "Other libraries have the following strategies for the cache.\n", + "\n", + "- [OpenSlide](https://openslide.org/) \n", + " - 1024 x 1024 x 30 bytes (30MiB) per file handle for cache ==> 160 (RGB) or 120 (ARGB) 256x256 tiles\n", + " - Not configurable\n", + "- [rasterio](https://rasterio.readthedocs.io/en/latest/)\n", + " - 5% of available system memory per process by default (e.g., 32 GB of free memory => 1.6 GB of cache memory allocated).\n", + " - Configurable through [environment module](https://rasterio.readthedocs.io/en/latest/api/rasterio.env.html)\n" + ] + }, + { + "cell_type": "markdown", + "id": "5148543f", + "metadata": {}, + "source": [ + "## Results\n", + "\n", + "cuCIM has a similar performance gain with the aligned case when the patch and tile layout are not aligned.\n", + "\n", + "We compared performance against OpenSlide and rasterio.\n", + "\n", + "For the cache memory size(capacity) setting, we used a similar approach with rasterio (5% of available system memory).\n", + "\n", + "\n", + "### System Information\n", + "\n", + "- OS: Ubuntu 18.04\n", + "- CPU: [Intel(R) Core(TM) i7-7800X CPU @ 3.50GHz](https://www.cpubenchmark.net/cpu.php?cpu=Intel+Core+i7-7800X+%40+3.50GHz&id=3037)\n", + "- Memory: 64GB (G-Skill DDR4 2133 16GB X 4)\n", + "- Storage\n", + " - SATA SSD: [Samsung SSD 850 EVO 1TB](https://www.samsung.com/us/computing/memory-storage/solid-state-drives/ssd-850-evo-2-5-sata-iii-1tb-mz-75e1t0b-am/)\n", + " \n", + "### Experiment Setup\n", + "+ Use read_region() APIs to read all patches (256x256 size each) of a whole slide image (.tif) at the largest resolution level (92,344 x 81,017. Internal tile size is 256 x 256 with 95% JPEG compression quality level) on multithread/multiprocess environment.\n", + " - Original whole slide image (.svs : 1.6GB) was converted into .tif file (3.2GB) using OpenSlide & tifffile library in this experiment (image2.tif).\n", + " * Original image can be downloaded from here(https://drive.google.com/drive/u/0/folders/0B--ztKW0d17XYlBqOXppQmw0M2M , TUPAC-TR-488.svs)\n", + "+ Two different job configurations\n", + " - multithreading: spread workload into multiple threads\n", + " - multiprocessing: spread workload into multiple processes\n", + "+ Two different read configurations for each job configuration\n", + " - unaligned/nocache: (256x256)-patch-reads start from (1,1). e.g., read the region (1,1)-(257,257) then, read the region (257,1)-(513,257), ...\n", + " - aligned: (256x256)-patch-reads start from (0,0). OpenSlide's internal cache mechanism does not affect this case.\n", + "+ Took about 10 samples due to the time to conduct the experiment so there could have some variation in the results.\n", + "+ Note that this experiment doesn\u2019t isolate the effect of system cache (page cache) that we excluded its effect on C++ API benchmark[discard_cache] so IO time itself could be short for both libraries.\n", + "\n", + "### Aligned Case (`per_process`, JPEG-compressed TIFF file)\n", + "\n", + "![image](https://user-images.githubusercontent.com/1928522/120849255-ae63b380-c52a-11eb-80c3-8411990e6c25.png)\n", + "\n", + "\n", + "### Unaligned Case (`per_process`, JPEG-compressed TIFF file)\n", + "\n", + "![image](https://user-images.githubusercontent.com/1928522/120849176-92601200-c52a-11eb-9c38-92e55c3d413a.png)\n", + "\n", + "### Overall Performance of `per_process` Compared with `no_cache` for Unaligned Case\n", + "\n", + "![image](https://user-images.githubusercontent.com/1928522/118345306-494b0e00-b4e8-11eb-88ca-c835c70aa037.png)\n", + "\n", + "\n", + "The detailed data is available [here](https://docs.google.com/spreadsheets/d/1eAqs24p25p6iIzZdUlnWNlk_RsrdRfEkIOYB9Xgu67c/edit?usp=sharing).\n" + ] + }, + { + "cell_type": "markdown", + "id": "d1c665c5", + "metadata": {}, + "source": [ + "\n", + "## Room for Improvement\n", + "\n", + "### Using of a Memory Pool\n", + "\n", + "`per_process` strategy performs better than `shared_memory` strategy, and both strategies perform less than `nocache` strategy when underlying tiles and patches are aligned.\n", + "- `shared_memory` strategy does some additional operations compared with `per_process` strategy, and both strategies have some overhead using cache (such as memory allocation for cache item/indirect function calls)\n", + "\n", + "=> All three strategies (including `nocache`) can have benefited if we allocate CPU/GPU memory for tiles from a fixed-sized cache memory pool (using [RMM](https://docs.rapids.ai/api/rmm/stable/basics.html) and/or [PMR](https://en.cppreference.com/w/cpp/memory/synchronized_pool_resource)) instead of calling malloc() to allocate memory.\n", + "\n", + "### Supporting Generator (iterator)\n", + "\n", + "When patches to read in an image can be determined in advance (inference use case), we can load/prefetch entire compressed/decompressed image data to the memory and provide Python generator(iterator) to get a series of patches efficiently for inference use cases. \n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "bf312216", + "metadata": {}, + "source": [ + "## Appendix\n", + "\n", + "### Experiment Code\n", + "\n", + "Please see https://github.com/rapidsai/cucim/blob/branch-21.12/experiments/Using_Cache/benchmark.py to check out the code used for the experiment." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/03-gabor-texture-classification.ipynb b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/03-gabor-texture-classification.ipynb new file mode 100644 index 00000000..fa879830 --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/03-gabor-texture-classification.ipynb @@ -0,0 +1,191 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Gabor Texture Classification \u2014 GPU vs CPU\n", + "\n", + "Classify textures (brick, grass, gravel) using a Gabor filter bank. Compares GPU-accelerated cuCIM/CuPy against CPU-based scikit-image/SciPy and shows the speedup." + ], + "id": "title-cell" + }, + { + "cell_type": "code", + "execution_count": null, + "id": "accessible-promise", + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "import cupy as cp\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from skimage import data\n", + "\n", + "durations = {}\n", + "for use_gpu in (False, True):\n", + " if use_gpu:\n", + " from cupyx.scipy import ndimage as ndi\n", + " from cucim.skimage.util import img_as_float32\n", + " from cucim.skimage.filters import gabor_kernel\n", + "\n", + " xp = cp\n", + " asnumpy = cp.asnumpy\n", + " device_name = \"gpu\"\n", + " else:\n", + " from scipy import ndimage as ndi\n", + " from skimage.util import img_as_float32\n", + " from skimage.filters import gabor_kernel\n", + "\n", + " xp = np\n", + " asnumpy = np.asarray\n", + " device_name = \"cpu\"\n", + "\n", + " def compute_feats(image, kernels):\n", + " feats = xp.zeros((len(kernels), 2), dtype=np.double)\n", + " for k, kernel in enumerate(kernels):\n", + " filtered = ndi.convolve(image, kernel, mode=\"wrap\")\n", + " feats[k, 0] = filtered.mean()\n", + " feats[k, 1] = filtered.var()\n", + " return feats\n", + "\n", + " def match(feats, ref_feats):\n", + " min_error = np.inf\n", + " min_i = None\n", + " for i in range(ref_feats.shape[0]):\n", + " error = xp.sum((feats - ref_feats[i, :]) ** 2)\n", + " if error < min_error:\n", + " min_error = error\n", + " min_i = i\n", + " return min_i\n", + "\n", + " tstart = time.time()\n", + "\n", + " # prepare filter bank kernels\n", + " kernels = []\n", + " for theta in range(4):\n", + " theta = theta / 4.0 * np.pi\n", + " for sigma in (1, 3):\n", + " for frequency in (0.05, 0.25):\n", + " kernel = gabor_kernel(\n", + " frequency, theta=theta, sigma_x=sigma, sigma_y=sigma\n", + " )\n", + " kernels.append(kernel.real)\n", + "\n", + " # shrink = (slice(0, None, 3), slice(0, None, 3))\n", + " brick = img_as_float32(xp.asarray(data.brick())) # [shrink]\n", + " grass = img_as_float32(xp.asarray(data.grass())) # [shrink]\n", + " gravel = img_as_float32(xp.asarray(data.gravel())) # [shrink]\n", + " image_names = (\"brick\", \"grass\", \"gravel\")\n", + " images = (brick, grass, gravel)\n", + "\n", + " # prepare reference features\n", + " ref_feats = xp.zeros((3, len(kernels), 2), dtype=np.double)\n", + " ref_feats[0, :, :] = compute_feats(brick, kernels)\n", + " ref_feats[1, :, :] = compute_feats(grass, kernels)\n", + " ref_feats[2, :, :] = compute_feats(gravel, kernels)\n", + "\n", + " print(\"Rotated images matched against references using Gabor filter banks:\")\n", + "\n", + " print(\"original: brick, rotated: 30deg, match result: \", end=\"\")\n", + " feats = compute_feats(ndi.rotate(brick, angle=190, reshape=False), kernels)\n", + " print(image_names[match(feats, ref_feats)])\n", + "\n", + " print(\"original: brick, rotated: 70deg, match result: \", end=\"\")\n", + " feats = compute_feats(ndi.rotate(brick, angle=70, reshape=False), kernels)\n", + " print(image_names[match(feats, ref_feats)])\n", + "\n", + " print(\"original: grass, rotated: 145deg, match result: \", end=\"\")\n", + " feats = compute_feats(ndi.rotate(grass, angle=145, reshape=False), kernels)\n", + " print(image_names[match(feats, ref_feats)])\n", + "\n", + " def power(image, kernel):\n", + " # Normalize images for better comparison.\n", + " image = (image - image.mean()) / image.std()\n", + " return xp.sqrt(\n", + " ndi.convolve(image, kernel.real, mode=\"wrap\") ** 2\n", + " + ndi.convolve(image, kernel.imag, mode=\"wrap\") ** 2\n", + " )\n", + "\n", + " # Plot a selection of the filter bank kernels and their responses.\n", + " results = []\n", + " kernel_params = []\n", + " for theta in (0, 1):\n", + " theta = theta / 4.0 * np.pi\n", + " for frequency in (0.1, 0.4):\n", + " kernel = gabor_kernel(frequency, theta=theta)\n", + " params = \"theta=%d,\\nfrequency=%.2f\" % (theta * 180 / np.pi, frequency)\n", + " kernel_params.append(params)\n", + " # Save kernel and the power image for each image\n", + " results.append((kernel, xp.stack([power(img, kernel) for img in images])))\n", + "\n", + " dur = time.time() - tstart\n", + " print(f\"Duration {device_name} = {dur} s\")\n", + " durations[device_name] = dur\n", + "\n", + " fig, axes = plt.subplots(nrows=5, ncols=4, figsize=(12, 14))\n", + " plt.gray()\n", + "\n", + " fig.suptitle(\"Image responses for Gabor filter kernels\", fontsize=12)\n", + "\n", + " axes[0][0].axis(\"off\")\n", + "\n", + " # Plot original images\n", + " for label, img, ax in zip(image_names, images, axes[0][1:]):\n", + " ax.imshow(asnumpy(img))\n", + " ax.set_title(label, fontsize=9)\n", + " ax.axis(\"off\")\n", + "\n", + " for label, (kernel, powers), ax_row in zip(kernel_params, results, axes[1:]):\n", + " # Plot Gabor kernel\n", + " ax = ax_row[0]\n", + " ax.imshow(asnumpy(kernel.real))\n", + " ax.set_ylabel(label, fontsize=7)\n", + " ax.set_xticks([])\n", + " ax.set_yticks([])\n", + "\n", + " # Plot Gabor responses with the contrast normalized for each filter\n", + " vmin = float(powers.min())\n", + " vmax = float(powers.max())\n", + " for patch, ax in zip(powers, ax_row[1:]):\n", + " ax.imshow(asnumpy(patch), vmin=vmin, vmax=vmax)\n", + " ax.axis(\"off\")\n", + " plt.tight_layout()\n", + " plt.show()\n", + "\n", + "print(f\"GPU Acceleration = {durations['cpu'] / durations['gpu']:0.4f}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "durable-johnson", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/04-random-walker-segmentation.ipynb b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/04-random-walker-segmentation.ipynb new file mode 100644 index 00000000..5d1e7968 --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/04-random-walker-segmentation.ipynb @@ -0,0 +1,138 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Random Walker Segmentation \u2014 GPU vs CPU\n", + "\n", + "GPU-accelerated image segmentation using the Random Walker algorithm. Compares cuCIM on GPU against scikit-image on CPU for a 1500x1500 synthetic image." + ], + "id": "title-cell" + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dense-standard", + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "import cupy as cp\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "durations = {}\n", + "for use_gpu in (False, True):\n", + " length = 1500\n", + " blob_size_fraction = 0.025\n", + " if use_gpu:\n", + " from cucim import skimage\n", + " from cucim.skimage.exposure import rescale_intensity\n", + " from cucim.skimage.segmentation import random_walker\n", + "\n", + " try:\n", + " from cucim.skimage.data import binary_blobs\n", + "\n", + " blobs = binary_blobs(\n", + " length=length, rng=1, blob_size_fraction=blob_size_fraction\n", + " )\n", + " except ImportError:\n", + " from skimage.data import binary_blobs\n", + "\n", + " blobs = cp.asarray(\n", + " binary_blobs(\n", + " length=length, rng=1, blob_size_fraction=blob_size_fraction\n", + " )\n", + " )\n", + " asnumpy = cp.asnumpy\n", + " xp = cp\n", + " device_name = \"gpu\"\n", + " else:\n", + " import skimage\n", + " from skimage.data import binary_blobs\n", + " from skimage.exposure import rescale_intensity\n", + " from skimage.segmentation import random_walker\n", + "\n", + " blobs = binary_blobs(\n", + " length=length, rng=1, blob_size_fraction=blob_size_fraction\n", + " )\n", + " asnumpy = np.asarray\n", + " xp = np\n", + " device_name = \"cpu\"\n", + "\n", + " # Generate noisy synthetic data\n", + " data = skimage.img_as_float(blobs)\n", + " print(f\"data.shape = {data.shape}\")\n", + " sigma = 0.3\n", + " data += xp.random.normal(loc=0, scale=sigma, size=data.shape)\n", + " data = rescale_intensity(data, in_range=(-sigma, 1 + sigma), out_range=(-1, 1))\n", + " data = data.astype(np.float32, copy=False)\n", + "\n", + " print(f\"data.dtype={data.dtype}\")\n", + " # The range of the binary image spans over (-1, 1).\n", + " # We choose the hottest and the coldest pixels as markers.\n", + " markers = xp.zeros(data.shape, dtype=np.uint)\n", + " markers[data < -0.95] = 1\n", + " markers[data > 0.95] = 2\n", + "\n", + " tstart = time.time()\n", + " # Run random walker algorithm\n", + " labels = random_walker(data, markers, beta=5, mode=\"cg\", tol=1e-5)\n", + "\n", + " dur = time.time() - tstart\n", + " durations[device_name] = dur\n", + " print(f\"Duration {device_name} = {dur} s\")\n", + "\n", + " # Plot results\n", + " fig, (ax1, ax2, ax3) = plt.subplots(\n", + " 1, 3, figsize=(12, 4.8), sharex=True, sharey=True\n", + " )\n", + " ax1.imshow(asnumpy(data[:200, :200]), vmin=-0.5, vmax=1.5, cmap=\"gray\")\n", + " ax1.axis(\"off\")\n", + " ax1.set_title(\"Noisy data\")\n", + " ax2.imshow(asnumpy(markers[:200, :200]), cmap=\"magma\")\n", + " ax2.axis(\"off\")\n", + " ax2.set_title(\"Markers\")\n", + " ax3.imshow(asnumpy(labels[:200, :200]), cmap=\"gray\")\n", + " ax3.axis(\"off\")\n", + " ax3.set_title(\"Segmentation\")\n", + "\n", + " fig.tight_layout()\n", + " plt.show()\n", + "\n", + "print(f\"GPU Acceleration = {durations['cpu'] / durations['gpu']:0.4f}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "educational-paraguay", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/05-vesselness-filter.ipynb b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/05-vesselness-filter.ipynb new file mode 100644 index 00000000..d58a8a62 --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/cucim_medical_imaging/05-vesselness-filter.ipynb @@ -0,0 +1,189 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Vesselness Filter \u2014 GPU vs CPU\n", + "\n", + "Detect vessel-like structures using Frangi and Sato filters. Compares GPU-accelerated cuCIM against CPU-based scikit-image." + ], + "id": "title-cell" + }, + { + "cell_type": "code", + "execution_count": null, + "id": "hearing-gazette", + "metadata": {}, + "outputs": [], + "source": [ + "from time import time\n", + "\n", + "import numpy as np\n", + "\n", + "from skimage import data\n", + "from skimage import color\n", + "from skimage.filters import meijering, sato, frangi, hessian\n", + "import matplotlib.pyplot as plt\n", + "\n", + "\n", + "def identity(image, **kwargs):\n", + " \"\"\"Return the original image, ignoring any kwargs.\"\"\"\n", + " return image\n", + "\n", + "\n", + "retina = data.retina()[200:-200, 200:-200]\n", + "\n", + "image = color.rgb2gray(retina)\n", + "image = image.astype(np.float32)\n", + "# image = np.tile(image, (4, 4)) # tile to increase size to roughly (4000, 4000)\n", + "print(f\"image.shape = {image.shape}\")\n", + "\n", + "kwargs = {\"sigmas\": [2], \"mode\": \"reflect\"}\n", + "fig, axes = plt.subplots(2, 5, figsize=[16, 8])\n", + "cmap = plt.cm.gray\n", + "tstart = time()\n", + "for i, black_ridges in enumerate([1, 0]):\n", + " for j, func in enumerate([identity, meijering, sato, frangi, hessian]):\n", + " kwargs[\"black_ridges\"] = black_ridges\n", + "\n", + " result = func(image, **kwargs)\n", + " vmin, vmax = np.percentile(result, q=[1, 99.5])\n", + " axes[i, j].imshow(result, cmap=cmap, vmin=vmin, vmax=vmax, aspect=\"auto\")\n", + " if i == 0:\n", + " axes[i, j].set_title(\n", + " [\n", + " \"Original\\nimage\",\n", + " \"Meijering\\nneuriteness\",\n", + " \"Sato\\ntubeness\",\n", + " \"Frangi\\nvesselness\",\n", + " \"Hessian\\nvesselness\",\n", + " ][j]\n", + " )\n", + " if j == 0:\n", + " axes[i, j].set_ylabel(\"black_ridges = \" + str(bool(black_ridges)))\n", + " axes[i, j].set_xticks([])\n", + " axes[i, j].set_yticks([])\n", + "print(f\"duration = {time() - tstart} s\")\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "electric-macintosh", + "metadata": {}, + "outputs": [], + "source": [ + "from time import time\n", + "\n", + "import cupy as cp\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from skimage import data\n", + "\n", + "durations = {}\n", + "for use_gpu in (False, True):\n", + " if use_gpu:\n", + " from cucim.skimage import color\n", + " from cucim.skimage.filters import meijering, sato, frangi, hessian\n", + "\n", + " xp = cp\n", + " asnumpy = cp.asnumpy\n", + " device_name = \"gpu\"\n", + " else:\n", + " from skimage import color\n", + " from skimage.filters import meijering, sato, frangi, hessian\n", + "\n", + " xp = np\n", + " asnumpy = np.asarray\n", + " device_name = \"cpu\"\n", + "\n", + " def identity(image, **kwargs):\n", + " \"\"\"Return the original image, ignoring any kwargs.\"\"\"\n", + " return image\n", + "\n", + " retina = data.retina()[200:-200, 200:-200]\n", + "\n", + " # transfer image to the GPU\n", + " retina = xp.asarray(retina)\n", + "\n", + " image = color.rgb2gray(retina)\n", + " image = image.astype(np.float32)\n", + " # image = cp.tile(image, (4, 4)) # tile to increase size to roughly (4000, 4000)\n", + " print(f\"image.shape = {image.shape}\")\n", + "\n", + " cmap = plt.cm.gray\n", + "\n", + " kwargs = {\"sigmas\": [2], \"mode\": \"reflect\"}\n", + " fig, axes = plt.subplots(2, 5, figsize=[16, 8])\n", + "\n", + " tstart = time()\n", + " for i, black_ridges in enumerate([1, 0]):\n", + " for j, func in enumerate([identity, meijering, sato, frangi, hessian]):\n", + " kwargs[\"black_ridges\"] = black_ridges\n", + "\n", + " result = func(image, **kwargs)\n", + "\n", + " # transfer back to host for visualization with Matplotlib\n", + " result_cpu = asnumpy(result)\n", + " vmin, vmax = map(float, xp.percentile(result, q=[1, 99.5]))\n", + " axes[i, j].imshow(\n", + " result_cpu, cmap=cmap, vmin=vmin, vmax=vmax, aspect=\"auto\"\n", + " )\n", + " if i == 0:\n", + " axes[i, j].set_title(\n", + " [\n", + " \"Original\\nimage\",\n", + " \"Meijering\\nneuriteness\",\n", + " \"Sato\\ntubeness\",\n", + " \"Frangi\\nvesselness\",\n", + " \"Hessian\\nvesselness\",\n", + " ][j]\n", + " )\n", + " if j == 0:\n", + " axes[i, j].set_ylabel(\"black_ridges = \" + str(bool(black_ridges)))\n", + " axes[i, j].set_xticks([])\n", + " axes[i, j].set_yticks([])\n", + " dur = time() - tstart\n", + " print(f\"duration = {dur} s\")\n", + " durations[device_name] = dur\n", + " plt.tight_layout()\n", + " plt.show()\n", + "\n", + "print(f\"GPU Acceleration = {durations['cpu'] / durations['gpu']:0.4f}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "african-infrared", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/docs-samples/data-science/gpu-accelerated-samples/cudf-pandas-stock-analysis.ipynb b/docs-samples/data-science/gpu-accelerated-samples/cudf-pandas-stock-analysis.ipynb new file mode 100644 index 00000000..85a7b47b --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/cudf-pandas-stock-analysis.ipynb @@ -0,0 +1,784 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "kcF9ZWvjSybR" + }, + "source": [ + "# Introduction\n", + "\n", + "cuDF is a Python GPU DataFrame library (built on the Apache Arrow columnar memory format) for loading, joining, aggregating, filtering, and otherwise manipulating tabular data using a DataFrame style API in the style of pandas.\n", + "\n", + "cuDF includes a pandas accelerator mode (`cudf.pandas`), enabling you to accelerate your pandas workflows without requiring any code change.\n", + "\n", + "This notebook highlights the impact of GPU-acceleration for common operations and analytical questions analyzing a real-world dataset of stock prices.\n", + "\n", + "For a deeper introduction and insight into how things work under the hood, we encourage you to run the [10 Minutes to RAPIDS cuDF's pandas accelerator mode Colab notebook](https://nvda.ws/rapids-cudf) or visit https://rapids.ai/cudf-pandas/." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "SH_h6ci1Sx0u" + }, + "source": [ + "# \u26a0\ufe0f Verify your setup\n", + "\n", + "First, we'll verify that you are running with an NVIDIA GPU." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "Y2vPCtXcCvUR", + "outputId": "e8fbe51e-64c3-4c2c-820b-52490cc0ee4a" + }, + "outputs": [], + "source": [ + "!nvidia-smi # this should display information about available GPUs" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "KP0oc3PboQDv" + }, + "source": [ + "With our GPU-enabled Colab runtime active, we're ready to go. cuDF is available by default in the GPU-enabled runtime." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "lv0guwQAcpgX" + }, + "source": [ + "# Download the data\n", + "\n", + "The data we'll be working with is a subset of the [USA 514 Stocks Prices NASDAQ NYSE dataset](https://www.kaggle.com/datasets/olegshpagin/usa-stocks-prices-ohlcv) from Kaggle.\n", + "\n", + "We'll start by downloading the dataset from NVIDIA's Public Google Cloud Storage bucket to provide faster download speeds to Colab. This should take under 30 seconds." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "WmOguzNUcw4F", + "outputId": "dbe56095-403b-4527-8809-472c0561d403" + }, + "outputs": [], + "source": [ + "!if [ ! -f \"usa_stocks_30m.parquet\" ]; then curl https://storage.googleapis.com/rapidsai/colab-data/usa_stocks_30m.parquet -o usa_stocks_30m.parquet; else echo \"usa_stocks_30m.parquet found\"; fi" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Pq01z9FvJjxR" + }, + "source": [ + "# Analysis using Standard Pandas\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "> **Microsoft Fabric note:** `%load_ext cudf.pandas` must be the very first executed statement \u2014 before any `import pandas`. In Fabric, pandas may be pre-imported by the environment. If so, add `cudf.pandas` to the Fabric session pre-run script or restart the kernel and ensure this cell runs first." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext cudf.pandas" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "rZ2Ac34yIqe8" + }, + "outputs": [], + "source": [ + "import pandas as pd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "QrbZubaGJmv5", + "outputId": "fd5e8f27-cdd4-44c4-fe58-8bb9238bb561" + }, + "outputs": [], + "source": [ + "%time df = pd.read_parquet(\"usa_stocks_30m.parquet\")\n", + "df.info()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 206 + }, + "id": "_JGCNEbQfjOQ", + "outputId": "7bdaea5d-20c9-4439-e8cd-3d81f4e073be" + }, + "outputs": [], + "source": [ + "df.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "SXAyFbEQfVBX" + }, + "source": [ + "We've got about 36M rows and 7 columns, covering stock price and trading-related attributes (`open`, `close`, etc.) for various companies (`ticker`) at what looks like 30 minute intervals.\n", + "\n", + "Let's look at the more detailed summary statistics for the the numeric data." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "SHfPsuqtN09j" + }, + "source": [ + "To fit this analysis on Colab's free tier without running out of the 12GB of CPU memory, we'll only use the first 18 million rows of this dataset. If you'd like to process the full dataset, you can [sign up a Colab Pro account](https://colab.research.google.com/signup)!\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "IlKNHgyONzrk" + }, + "outputs": [], + "source": [ + "df = df.iloc[:18000000]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 336 + }, + "id": "DqqjcfcfJnvy", + "outputId": "2592a7f6-8777-4515-9737-4e0c48228dc8" + }, + "outputs": [], + "source": [ + "%time df.describe()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "iuatY3BFJzPO" + }, + "source": [ + "A little slow, but it gets the job done.\n", + "\n", + "Most of the price-related columns look well behaved with no crazy outliers based on the max values. Volume has a high variance, but that makes sense as some stocks will be traded **much** more often than others. We've also got a wide time range in this data, spanning 1998 to 2024.\n", + "\n", + "The time-range information feels important. This is per-stock data, so some stocks will have \"entered\" the dataset at different times. We should really group by each ticker for our analysis.\n", + "\n", + "To start, we can investigate the time periods for various stocks in the data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 522 + }, + "id": "b0tMz4nuJva6", + "outputId": "ff9c37b8-bbc8-4178-c1bc-ce5951e7df31" + }, + "outputs": [], + "source": [ + "%time df.groupby(\"ticker\").agg({\"datetime\": [\"min\", \"max\", \"count\"]})" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "JznDAW8tKOtg" + }, + "source": [ + "As expected, there's a pretty significant difference across stocks. Since this determines how many records exist for each ticker, we'll need to take that into account for future analysis.\n", + "\n", + "One way to do that is to use the time periods as part of any grouping or rolling window analysis we do on this data. For example, we can look at the minimum and maximum of each stock ticker at various time frequencies to understand trends." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 522 + }, + "id": "HZzNnYcfJ-nR", + "outputId": "67760552-e809-4e28-8ede-d317aa91480b" + }, + "outputs": [], + "source": [ + "%%time\n", + "df[[\"year\", \"week\", \"day\"]] = df.datetime.dt.isocalendar()\n", + "df.groupby([\"ticker\", \"year\", \"week\"]).agg({\"close\": [\"min\", \"max\"]})" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "aXEYOj-BLfXP" + }, + "source": [ + "Okay, that was a little slow, but it's manageable. We can easily aggregate the data up to weekly values and investigate further.\n", + "\n", + "But what if we want to ask slightly more complex questions?\n", + "\n", + "## More complex analysis\n", + "\n", + "\n", + "Let's say we want the daily rolling average of all of these values for each stock. With this, we could investigate how each 30 minute interval compares to the rolling average over the course of a day.\n", + "\n", + "In this **specific** demo dataset, we have exactly one record per 30 minutes, so we could use a fixed window size of 12 (12 periods per market day). But in practice, we often have messy data with inconsistent time frequencies -- not to mention some missing or duplicated data. Fixed time windows can potentially corrupt our analysis.\n", + "\n", + "Fortunately, solving this problem is actually pretty easy with pandas. We can use a fixed **time window** per ticker rather than fixed number of records.\n", + "\n", + "Unfortunately, it's **pretty** slow." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 241 + }, + "id": "42wbyIdaLh8w", + "outputId": "1931ebc7-dff4-45bb-fef7-815dddc8bde8" + }, + "outputs": [], + "source": [ + "%time result = df.set_index(\"datetime\").sort_index().groupby(\"ticker\").rolling(\"1D\").mean().reset_index()\n", + "result.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "7iuGSOSEgJZy" + }, + "source": [ + "**About** a 20 to 30 seconds just for a single question. That's not ideal. Especially for what comes next.\n", + "\n", + "Now that we have an average price per day, let's do this again to get the Simple Moving Averages (SMA) usually used in stock analysis. We'll first do the 50 Day SMA and then the 200 Day SMA." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 241 + }, + "id": "K1rexa2_gCTW", + "outputId": "2e265b6d-eba3-42ea-ff13-a6d2cbbdbd90" + }, + "outputs": [], + "source": [ + "%time fiftyDay = df.set_index(\"datetime\").sort_index().groupby(\"ticker\").rolling(\"50D\").mean().reset_index()\n", + "fiftyDay.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 241 + }, + "id": "-iyZercThFBb", + "outputId": "e80b6224-0d74-4d68-a396-0b0c3de04346" + }, + "outputs": [], + "source": [ + "%time twoHunDay = df.set_index(\"datetime\").sort_index().groupby(\"ticker\").rolling(\"200D\").mean().reset_index()\n", + "twoHunDay.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "9YQj8Kg1LlCs" + }, + "source": [ + "This took just about a minute and a half to get to just about where you can plot. \n", + "\n", + "Perfect time to switch to cudf.pandas. Let's do that and then run the **same code**." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "s4CC8LjKiJds" + }, + "source": [ + "# Analysis with cuDF Pandas\n", + "\n", + "Typically, you should load the `cudf.pandas` extension as the first step in your notebook, before importing any modules. Here, we explicitly restart the kernel to simulate that behavior." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Timing with cudf.pandas\n", + "\n", + "Since cudf.pandas was loaded at the top, all operations above were already GPU-accelerated.\n", + "The cells below re-run with explicit timing to show the speedup." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "rwqZKmGhiYXY" + }, + "source": [ + "We'll run the same code as above to get a feel what GPU-acceleration brings to pandas workflows." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "ci2waUePPhuo", + "outputId": "2a859dbe-9b29-4cab-a28a-0027536f0acd" + }, + "outputs": [], + "source": [ + "%time df = pd.read_parquet(\"usa_stocks_30m.parquet\")\n", + "df.info()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 206 + }, + "id": "BvQrz_ImiiLp", + "outputId": "a4d804db-31a5-47de-eb15-e188edebcf16" + }, + "outputs": [], + "source": [ + "df.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "BL5N9mYUil2u" + }, + "source": [ + "Let's look at the more detailed summary statistics for the the numeric data like we did before." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "vwzmxbzLjQ9B", + "outputId": "3a68612b-b905-45d7-9abc-9bc978cc7226" + }, + "outputs": [], + "source": [ + "df = df.iloc[:18000000]\n", + "df.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 336 + }, + "id": "JXl-zPKQijwp", + "outputId": "203f7abc-32d9-45b5-e7a1-e2649a9f4a92" + }, + "outputs": [], + "source": [ + "%time df.describe()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "szqehmcfi-BY" + }, + "source": [ + "First things first, we can see that results are the same, even though we're now using the GPU. Great.\n", + "\n", + "Second, that was much quicker, even with a little overhead from this being my first GPU-accelerated operation on the data.\n", + "\n", + "Let's do the groupby aggregations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 522 + }, + "id": "lh2h5u_FRv75", + "outputId": "b345e193-cf23-4f2c-9720-c87d9eafe60e" + }, + "outputs": [], + "source": [ + "%time df.groupby(\"ticker\").agg({\"datetime\": [\"min\", \"max\", \"count\"]})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 522 + }, + "id": "pSCoyXtoR2Kb", + "outputId": "d3d4c428-bcdd-49c7-9d6b-8824ca7e077b" + }, + "outputs": [], + "source": [ + "%%time\n", + "df[[\"year\", \"week\", \"day\"]] = df.datetime.dt.isocalendar()\n", + "df.groupby([\"ticker\", \"year\", \"week\"]).agg({\"close\": [\"min\", \"max\"]})" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "mQF_atfBSBTz" + }, + "source": [ + "20-50x faster with the same code. Nice!\n", + "\n", + "Let's do the long-running groupby rolling operation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 241 + }, + "id": "bbRTBWNYR8et", + "outputId": "909e4fc8-25d4-496a-b45d-4577115ef766" + }, + "outputs": [], + "source": [ + "%time result = df.set_index(\"datetime\").sort_index().groupby(\"ticker\").rolling(\"1D\").mean().reset_index()\n", + "result.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "iDG9Iz0gSJB4" + }, + "source": [ + "Nice! 25+ seconds down to about < 2 seconds.\n", + "\n", + "Let's finish off those SMAs. While we do that, let's profile the activity so you can see where the CPU and GPU\n", + "\n", + "We've only used pandas so far, and things have worked smoothly.\n", + "\n", + "But workflows often use multiple libraries, many of which are designed to accept pandas inputs. Fortunately, cudf.pandas works here, too." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 241 + }, + "id": "3PxHzWi1DFFJ", + "outputId": "540473bc-9b24-46c6-b4ac-35a0e1e37eeb" + }, + "outputs": [], + "source": [ + "%time fiftyDay = df.set_index(\"datetime\").sort_index().groupby(\"ticker\").rolling(\"50D\").mean().reset_index()\n", + "fiftyDay.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 241 + }, + "id": "it2UxeC_Pko0", + "outputId": "6a505b35-ac20-4fc2-d464-887a71149b4a" + }, + "outputs": [], + "source": [ + "%time twoHunDay = df.set_index(\"datetime\").sort_index().groupby(\"ticker\").rolling(\"200D\").mean().reset_index()\n", + "twoHunDay.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "ZpUGGkXMg8WJ" + }, + "source": [ + "Awesome, `cudf.pandas` got the same workflow done in about 10 seconds!\n", + "\n", + "It is also important to understand that this speed difference increases as you add all the cells and increase complexity. To do all a 200 Day SMA on all 36M cells will take **<12 seconds** on a the T4 using `cuDF.pandas` and well over a minute on CPU." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "RHWsDYuSaAmX" + }, + "source": [ + "## Using third-party libraries with cudf.pandas\n", + "\n", + "\n", + "When using cudf.pandas, you can pass pandas objects to third-party libraries just like you would when using regular Pandas.\n", + "\n", + "For example, we can pass this dataset to the visualization library plotnine and plot the yearly average closing price of the ticker `GOOG`. Let's join the results of our groupbys above, form the data so that `plotnine` can ingest it, and chart this data out!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "-0GIHLKSWqEf" + }, + "outputs": [], + "source": [ + "result = result.join(fiftyDay, rsuffix='_50')\n", + "result = result.join(twoHunDay, rsuffix='_200')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 429 + }, + "id": "GkUJTsI8XQzO", + "outputId": "65f9751a-e25f-4938-b3d2-0976e150039a" + }, + "outputs": [], + "source": [ + "result.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 429 + }, + "id": "vA_i-TY5Yqy3", + "outputId": "f0ee9ebc-feb1-479c-e5af-b96469584b95" + }, + "outputs": [], + "source": [ + "goog_closing_value = result.loc[result.ticker == \"GOOG\"]\n", + "goog_closing_value.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 206 + }, + "id": "ecmO09tFZOiE", + "outputId": "17df5cc5-b930-4e28-d97f-579d7fe13b4e" + }, + "outputs": [], + "source": [ + "goog_closing_value_p9 = goog_closing_value.melt(\n", + " id_vars=\"datetime\",\n", + " value_vars=[\"close\", \"close_50\", \"close_200\"],\n", + " var_name=\"SMA\",\n", + " value_name=\"price\"\n", + ").dropna()\n", + "\n", + "goog_closing_value_p9.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "YXxlTG4HbBTD" + }, + "source": [ + "Now let's bring in `plotnine` and visually analyze the stock and it's performance versus the 50 and 200 day moving averages!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "iDuSN_KmakTn" + }, + "outputs": [], + "source": [ + "from plotnine import *" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 588 + }, + "id": "EtUNRJsDYmF-", + "outputId": "4ed1ef9d-ecd9-44af-b0fa-b511602c2057" + }, + "outputs": [], + "source": [ + "# Gallery, lines\n", + "(\n", + " ggplot(goog_closing_value_p9, aes(x=\"datetime\", y=\"price\", color=\"SMA\"))\n", + " + geom_line()\n", + " # Styling\n", + " + scale_x_datetime(date_breaks=\"1 year\", date_labels=\"%Y\")\n", + " + theme_538()\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "J-rHjLf4YWNc" + }, + "source": [ + "# Summary" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "RikEB9oiYWDS" + }, + "source": [ + "With `cudf.pandas`, you can keep using pandas as your primary dataframe library if you enjoy using it but want faster performance. When things start to get a little slow, just load the `cudf.pandas` and run your existing code on a GPU!\n", + "\n", + "If you like Google Colab and want to get peak cudf.pandas performance to process even larger datasets, Google Colab's paid tier includes both L4 and A100 GPUs (in addition to the T4 GPU this demo notebook is using).\n", + "\n", + "To learn more about cudf.pandas, we encourage you to visit [rapids.ai/cudf-pandas](https://rapids.ai/cudf-pandas)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "h3O4hmE0YZoy" + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "accelerator": "GPU", + "colab": { + "gpuType": "T4", + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/docs-samples/data-science/gpu-accelerated-samples/cuml-scikit-learn-accelerator-demo.ipynb b/docs-samples/data-science/gpu-accelerated-samples/cuml-scikit-learn-accelerator-demo.ipynb new file mode 100644 index 00000000..eb7df329 --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/cuml-scikit-learn-accelerator-demo.ipynb @@ -0,0 +1,1101 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "ZhTWAPo7X5-I" + }, + "source": [ + "# Getting Started with cuML's accelerator mode (cuml.accel)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "BinKvOgMYOCp" + }, + "source": [ + "cuML is a Python GPU library for accelerating machine learning models using a scikit-learn-like API.\n", + "\n", + "cuML now has an accelerator mode (cuml.accel) which allows you to bring accelerated computing to existing workflows with zero code changes required. In addition to scikit-learn, cuml.accel also provides acceleration to algorithms found in umap-learn (UMAP) and hdbscan (HDBSCAN).\n", + "\n", + "This notebook is a brief introduction to cuml.accel." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "ylUzZjM-mRMH" + }, + "source": [ + "# \u26a0\ufe0f Verify your setup" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "e_AtDRMEZQd-" + }, + "source": [ + "First, we'll verfiy that we are running on an NVIDIA GPU:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext cuml.accel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "OKWtiMvJAS60", + "outputId": "3ad02328-d81d-4296-f2f9-9c97baf40114" + }, + "outputs": [], + "source": [ + "!nvidia-smi # this should display information about available GPUs" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "qC3fevZecnns" + }, + "source": [ + "With classical machine learning, there is a wide range of interesting problems we can explore. In this tutorial we'll examine 3 of the more popular use cases: classification, clustering, and dimensionality reduction." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "M37-8qsDa2Pe" + }, + "source": [ + "# Classification" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "vT5RNLwdce-O" + }, + "source": [ + "Let's load a dataset and see how we can use scikit-learn to classify that data. For this example we'll use the Coverage Type dataset, which contains a number of features that can be used to predict forest cover type, such as elevation, aspect, slope, and soil-type.\n", + "\n", + "More information on this dataset can be found at https://archive.ics.uci.edu/dataset/31/covertype." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "PKt2Lje5lYQw" + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "from sklearn.ensemble import RandomForestClassifier\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.metrics import classification_report, accuracy_score" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "rLHSPlLnv-1y" + }, + "outputs": [], + "source": [ + "url = \"https://archive.ics.uci.edu/ml/machine-learning-databases/covtype/covtype.data.gz\"\n", + "\n", + "# Column names for the dataset (from UCI Covertype description)\n", + "columns = ['Elevation', 'Aspect', 'Slope', 'Horizontal_Distance_To_Hydrology', 'Vertical_Distance_To_Hydrology',\n", + " 'Horizontal_Distance_To_Roadways', 'Hillshade_9am', 'Hillshade_Noon', 'Hillshade_3pm',\n", + " 'Horizontal_Distance_To_Fire_Points', 'Wilderness_Area1', 'Wilderness_Area2', 'Wilderness_Area3',\n", + " 'Wilderness_Area4', 'Soil_Type1', 'Soil_Type2', 'Soil_Type3', 'Soil_Type4', 'Soil_Type5', 'Soil_Type6',\n", + " 'Soil_Type7', 'Soil_Type8', 'Soil_Type9', 'Soil_Type10', 'Soil_Type11', 'Soil_Type12', 'Soil_Type13',\n", + " 'Soil_Type14', 'Soil_Type15', 'Soil_Type16', 'Soil_Type17', 'Soil_Type18', 'Soil_Type19', 'Soil_Type20',\n", + " 'Soil_Type21', 'Soil_Type22', 'Soil_Type23', 'Soil_Type24', 'Soil_Type25', 'Soil_Type26', 'Soil_Type27',\n", + " 'Soil_Type28', 'Soil_Type29', 'Soil_Type30', 'Soil_Type31', 'Soil_Type32', 'Soil_Type33', 'Soil_Type34',\n", + " 'Soil_Type35', 'Soil_Type36', 'Soil_Type37', 'Soil_Type38', 'Soil_Type39', 'Soil_Type40', 'Cover_Type']\n", + "\n", + "data = pd.read_csv(url, header=None)\n", + "data.columns=columns" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "53P_F5oHmh9F", + "outputId": "da81ae3e-9725-4a99-d11e-e65f57a99d73" + }, + "outputs": [], + "source": [ + "data.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "7Mz-yThWmlqg" + }, + "source": [ + "Next, we'll separate out the classification variable (Cover_Type) from the rest of the data. This is what we will aim to predict with our classification model. We can also split our dataset into training and test data using the scikit-learn train_test_split function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "heVSKrDLxdN3" + }, + "outputs": [], + "source": [ + "X, y = data.drop('Cover_Type', axis=1), data['Cover_Type']\n", + "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "OVl1oBxxm44k" + }, + "source": [ + "Now that we have our dataset split, we're ready to run a model. To start, we will just run the model using the sklearn library with a starting max depth of 5 and all of the features. Note that we can set n_jobs=-1 to utilize all available CPU cores for fitting the trees -- this will ensure we get the best performance possible on our system's CPU. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 115 + }, + "id": "T0Y2HUgykyLY", + "outputId": "5977bcbc-9c5e-45bd-9714-3d4814c22e72" + }, + "outputs": [], + "source": [ + "%%time\n", + "\n", + "clf = RandomForestClassifier(n_estimators=100, max_depth=5, max_features=1.0, n_jobs=-1)\n", + "clf.fit(X_train, y_train)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "_B7glDK8nWMQ" + }, + "source": [ + "In about 2 minutes, we were able to fit our tree model using scikit-learn. This is not bad! Let's use the model we just trained to predict coverage types in our test dataset and take a look at the accuracy of our model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "SLP-K7qynwWg", + "outputId": "d9239568-ca10-48e8-e6bb-4436d5e028a4" + }, + "outputs": [], + "source": [ + "y_pred = clf.predict(X_test)\n", + "accuracy_score(y_test, y_pred)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "UAFu5rDwn6VK" + }, + "source": [ + "We can also print out a full classification report to better understand how we predicted different Coverage_Type categories." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "AO3owg5PnwYW", + "outputId": "3dea9a3b-b08d-4206-952d-dcd30b7fa20d" + }, + "outputs": [], + "source": [ + "print(classification_report(y_test, y_pred))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "FRdxW-MpoJAt" + }, + "source": [ + "With scikit-learn, we built a model that was able to be trained in just a couple minutes. From the accuracy report, we can see that we predicted the correct class around 70% of the time, which is not bad but could certainly be improved.\n", + "\n", + "Often we want to run several different random forest models in order to optimize our hyperparameters. For example, we may want to increase the number of estimators, or modify the maximum depth of our tree. When running dozens or hundreds of different hyperparameter combinations, things start to become quite slow and iteration takes a lot longer.\n", + "\n", + "We provide some sample code utilizing GridSearchCV below to show what this process might look like. All of these combinations would take a LONG time to run if we spend 2 minutes fitting each model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "4qFknPIWpXcs" + }, + "outputs": [], + "source": [ + "\"\"\"\n", + "from sklearn.model_selection import GridSearchCV\n", + "\n", + "# Define the parameter grid to search over\n", + "param_grid = {\n", + " 'n_estimators': [50, 100, 200],\n", + " 'max_depth': [None, 10, 20, 30],\n", + " 'min_samples_split': [2, 5, 10],\n", + " 'min_samples_leaf': [1, 2, 4],\n", + " 'max_features': ['auto', 'sqrt', 'log2'],\n", + " 'bootstrap': [True, False]\n", + "}\n", + "\n", + "grid_search = GridSearchCV(estimator=clf, param_grid=param_grid, cv=5)\n", + "grid_search.fit(X_train, y_train)\n", + "\"\"\"" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "5-BtB1RdqKp0" + }, + "source": [ + "Now let's load cuml.accel and try running the same code again to see what kind of acceleration we can get." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# cuml.accel already loaded at top of notebook" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "TNLZVjtxqTmc" + }, + "source": [ + "After loading the IPython magic, we need to import the sklearn estimators we wish to use again." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "MVWHm7H9qS2T" + }, + "outputs": [], + "source": [ + "from sklearn.ensemble import RandomForestClassifier" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 115 + }, + "id": "8v0QxJXrViiR", + "outputId": "5510c2ad-5ea7-4cb7-c3b0-0527d836fc6b" + }, + "outputs": [], + "source": [ + "%%time\n", + "\n", + "clf = RandomForestClassifier(n_estimators=100, max_depth=5, max_features=1.0, n_jobs=-1)\n", + "clf.fit(X_train, y_train)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "fqJC2jtmqpva" + }, + "source": [ + "That was much faster! Using cuML we're able to train this random forest model in just seconds instead of minutes. One thing to note is that cuML's implementation of RandomForestClassifier doesn't utilize the `n_jobs` parameter like scikit-learn, but we still accept it which makes it easier to use this accelerator with zero code changes.\n", + "\n", + "Let's take a look at the same accuracy score and classification report to compare the model's performance." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "eVre6kav6iaS", + "outputId": "2c240c5f-b7c2-4eac-98b5-aee48bd4dcc7" + }, + "outputs": [], + "source": [ + "y_pred = clf.predict(X_test)\n", + "cr = classification_report(y_test, y_pred)\n", + "print(cr)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "8xkoz247VsIX" + }, + "source": [ + "Out of the box, the model performed about the same as the scikit-learn implementation. Because this model ran so much faster, we can quickly iterate on the hyperparameter configuration and find a model that performs better with excellent speedups." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 115 + }, + "id": "mfZamg7FVoPe", + "outputId": "19d44cd9-9863-495b-ac8d-067118fd9c22" + }, + "outputs": [], + "source": [ + "%%time\n", + "\n", + "clf = RandomForestClassifier(n_estimators=100, max_depth=30, max_features=1.0, n_jobs=-1)\n", + "clf.fit(X_train, y_train)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "E97LObEYVocu", + "outputId": "e560b1a6-1b12-4cb5-bc66-09b3530af692" + }, + "outputs": [], + "source": [ + "y_pred = clf.predict(X_test)\n", + "print(classification_report(y_test, y_pred))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "DmcwmqN4q4cv" + }, + "source": [ + "With a model that runs in just seconds, we can perform hyperparameter optimization using a method like the grid search shown above, and have results in just minutes instead of hours." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "9IXkjC6MLZoj" + }, + "source": [ + "# CPU Fallback\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "9zSIAHaKn7oa" + }, + "source": [ + "There are some algorithms and functionality from scikit-learn, UMAP, and HDBSCAN that are *not* implemented in cuML. For cases where the underlying functionality is not supported on GPU, the cuML accelerator will gracefully fall back and execute on the CPU instead.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "qtB_3HlFn-wx" + }, + "outputs": [], + "source": [ + "from sklearn.neighbors import KernelDensity\n", + "import numpy as np\n", + "\n", + "X = np.concatenate((np.random.normal(0, 1, 10000),\n", + " np.random.normal(5, 1, 10000)))[:, np.newaxis]\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 80 + }, + "id": "GYmEAsvXn-yv", + "outputId": "9e868e87-514c-4910-d19a-50cc4a67ad7d" + }, + "outputs": [], + "source": [ + "kde = KernelDensity(kernel='gaussian', bandwidth=0.5)\n", + "kde.fit(X)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "rBB16r0in-1F", + "outputId": "5a5cfef5-ea8e-4240-e866-4138b7f700b1" + }, + "outputs": [], + "source": [ + "print(kde.score_samples(X))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "OWcDT_9Cm6Mk" + }, + "source": [ + "Next, let's restart the kernel to unload the accelerator extension and observe the same performance comparisons on a few other algorithms." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# (cuml.accel already loaded at top)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "9Ql7UCNUW1AO" + }, + "source": [ + "We'll now take a look at a clustering example using HDBSCAN.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "pdjvwlZ-DaoW" + }, + "source": [ + "# Clustering" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "70sdcTIdfzkK" + }, + "source": [ + "Clustering is an important data science workflow because it helps uncover hidden patterns and structures within data without requiring labeled outcomes. In practice, with high dimensional data it can be difficult to discern whether the clusters we've chosen are good or not. One way to determine the quality of our clustering is with sklearn's [silhouette score](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_score.html#sklearn.metrics.silhouette_score), which we'll examine shortly.\n", + "\n", + "HDBSCAN is a popular density-based clustering algorithm that is highly flexible. We'll load a toy sklearn dataset to illustrate how HDBSCAN can be accelerated with cuml.accel." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "oIlRPHZ1DZ7q" + }, + "outputs": [], + "source": [ + "import hdbscan\n", + "import matplotlib.pyplot as plt\n", + "from sklearn.datasets import make_blobs\n", + "from sklearn.metrics import silhouette_score" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 449 + }, + "id": "RqXkNWRjDZ9n", + "outputId": "57435e6c-a303-4531-d64b-87ba14cce410" + }, + "outputs": [], + "source": [ + "N = 20000\n", + "K = 100\n", + "\n", + "X, y = make_blobs(\n", + " n_samples=N,\n", + " n_features=K,\n", + " centers=5,\n", + " cluster_std=[3,1,2,1.5,0.5],\n", + " random_state=42\n", + ")\n", + "\n", + "plt.scatter(X[:, 0], X[:, 1], c=y)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 184 + }, + "id": "536jWwBWg1Ou", + "outputId": "c589e021-e14e-4848-d8b5-2302c00d4ef7" + }, + "outputs": [], + "source": [ + "clusterer = hdbscan.HDBSCAN()\n", + "%time clusterer.fit(X)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "wrAzGgqrxHtY", + "outputId": "05493c84-d6cf-469a-e533-db36558c76ec" + }, + "outputs": [], + "source": [ + "print(silhouette_score(X, clusterer.labels_))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# cuml.accel already loaded at top of notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "vuwUVZeihJWh" + }, + "outputs": [], + "source": [ + "import hdbscan" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 115 + }, + "id": "BfDHJoiSDZ_9", + "outputId": "8e7f9de5-e545-4051-e1a7-698b734577cf" + }, + "outputs": [], + "source": [ + "clusterer = hdbscan.HDBSCAN()\n", + "%time clusterer.fit(X)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "B_tOMHd8y5ai" + }, + "source": [ + "Here we go from around 45 seconds to 1 second to fit the clustering model! This is a massive speed-up we got just from loading the `cuml.accel` extension." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "zZNpIYCqkk_8", + "outputId": "61451ac0-e21c-4b85-e718-0fbad3128bef" + }, + "outputs": [], + "source": [ + "print(silhouette_score(X, clusterer.labels_))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "91a7zkUIoE5d" + }, + "source": [ + "\n", + "It's important to note that on real-world datasets, the silhouette score produced by the GPU and CPU implementations of HDBSCAN will often have slight differences. The cuML implementation of HDBSCAN should provide equivalent results, but it is normal for the actual clusters to vary slightly when dealing with complex datasets." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "k6AvkVh7QMlw" + }, + "source": [ + "Lastly, let's take a look at how we can use cuml's accelerator mode for a third popular machine learning task -- dimensionality reduction. We'll restart the kernel to unload the extension yet again.\n", + "\n", + "Keep in mind that we don't normally need to restart the kernel when using `cuml.accel`, we just do it for the sake of showing the speed-ups in this demo. In practice, you'd just load the accelerator one time up front and be set." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# (cuml.accel already loaded at top)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "93qD18LqDiOj" + }, + "source": [ + "# Dimensionality Reduction" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "n0mwvzC0rdq4" + }, + "source": [ + "UMAP is a popular dimensionality reduction technique that is used for both data visualization and as preprocessing for downstream modeling due to its ability to balance preserving both local and global structure of high-dimensional data. To learn more about how it works, visit the [UMAP documentation](https://umap-learn.readthedocs.io/en/latest/).\n", + "\n", + "To explore how cuML can accelerate UMAP, let's load in another dataset from UCI. We'll use the Human Activity Recognition (HAR) dataset, which was created from recordings of 30 subjects performing activities of daily living (ADL) while carrying a waist-mounted smartphone with embedded inertial sensors." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "JOeW2dt9D2IZ", + "outputId": "3f687e7c-8652-4b8e-bfc3-43cd1b291531" + }, + "outputs": [], + "source": [ + "!wget https://archive.ics.uci.edu/ml/machine-learning-databases/00240/UCI%20HAR%20Dataset.zip -O /tmp/HAR_data.zip" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "TR78I9VuEGlu", + "outputId": "aad52a74-9c40-4be5-dc30-d6cda2b0f971" + }, + "outputs": [], + "source": [ + "!unzip /tmp/HAR_data.zip -d /tmp/HAR_data/" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import ssl\n", + "import urllib.request\n", + "import zipfile\n", + "import certifi\n", + "\n", + "HAR_DIR = \"/tmp/HAR_data/UCI HAR Dataset\"\n", + "HAR_URL = \"https://archive.ics.uci.edu/ml/machine-learning-databases/00240/UCI%20HAR%20Dataset.zip\"\n", + "\n", + "if not os.path.exists(HAR_DIR):\n", + " os.makedirs(\"/tmp/HAR_data\", exist_ok=True)\n", + " print(\"Downloading UCI HAR Dataset (~60MB)...\")\n", + " ctx = ssl.create_default_context(cafile=certifi.where())\n", + " zip_path = \"/tmp/HAR_data/har.zip\"\n", + " with urllib.request.urlopen(HAR_URL, context=ctx) as response:\n", + " with open(zip_path, 'wb') as f:\n", + " f.write(response.read())\n", + " with zipfile.ZipFile(zip_path, 'r') as z:\n", + " z.extractall(\"/tmp/HAR_data\")\n", + " os.remove(zip_path)\n", + " print(f\"Extracted to {HAR_DIR}\")\n", + "else:\n", + " print(f\"HAR dataset already exists: {HAR_DIR}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "QjGw0WDwCu2_" + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "X_train = pd.read_csv(\"/tmp/HAR_data/UCI HAR Dataset/train/X_train.txt\", sep=r\"\\s+\", header=None)\n", + "y_train = pd.read_csv(\"/tmp/HAR_data/UCI HAR Dataset/train/y_train.txt\", sep=r\"\\s+\", header=None)\n", + "X_test = pd.read_csv(\"/tmp/HAR_data/UCI HAR Dataset/test/X_test.txt\", sep=r\"\\s+\", header=None)\n", + "y_test = pd.read_csv(\"/tmp/HAR_data/UCI HAR Dataset/test/y_test.txt\", sep=r\"\\s+\", header=None)\n", + "labels = pd.read_csv(\"/tmp/HAR_data/UCI HAR Dataset/activity_labels.txt\", sep=r\"\\s+\", header=None)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "oesYc97pQ1vr", + "outputId": "e7dca260-5abb-4f27-d61e-e336955b5271" + }, + "outputs": [], + "source": [ + "X_train.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "F5Cle6mGOH7S" + }, + "source": [ + "Let's take a look at the activity labels to better understand the data we're working with. We can see that the sensors have grouped activities into 6 different classes." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 238 + }, + "id": "22YYMRvXOUUz", + "outputId": "b15af6d7-8c52-474e-f102-f573dd7f2d2e" + }, + "outputs": [], + "source": [ + "labels" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "8rA7SsEpCu5M" + }, + "outputs": [], + "source": [ + "from sklearn.preprocessing import StandardScaler\n", + "\n", + "# Scale the data before applying UMAP\n", + "scaler = StandardScaler()\n", + "X_train_scaled = scaler.fit_transform(X_train)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "_9ltJhG7tQt0" + }, + "source": [ + "Let's run UMAP with some basic parameters and explore a lower-dimensionality projection of this dataset." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "t4sR7x3sF_9a" + }, + "outputs": [], + "source": [ + "import umap\n", + "umap_model = umap.UMAP(n_neighbors=15, n_components=2, random_state=42, min_dist=0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "cGD49PFVlyrz", + "outputId": "81f5a19c-d2c6-498c-c673-facc3ced964d" + }, + "outputs": [], + "source": [ + "%%time\n", + "\n", + "# Fit UMAP model to the data\n", + "X_train_umap = umap_model.fit_transform(X_train_scaled)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "kwJC2woSl4U9" + }, + "source": [ + "It's often quite interesting to visualize the resulting projection of the embeddings created by UMAP. In this case, let's take a look at the now 2-dimensional dataset." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 718 + }, + "id": "PDRxaBEdF__X", + "outputId": "2ddd2b8b-8b1a-48e4-e06b-b5ffe0245f27" + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "# Plot the UMAP result\n", + "plt.figure(figsize=(10, 8))\n", + "plt.scatter(X_train_umap[:, 0], X_train_umap[:, 1], c=y_train.values.ravel(), cmap='Spectral', s=10)\n", + "plt.colorbar(label=\"Activity\")\n", + "plt.title(\"UMAP projection of the UCI HAR dataset\")\n", + "plt.xlabel(\"UMAP Component 1\")\n", + "plt.ylabel(\"UMAP Component 2\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "gKFBLuE7mFrs" + }, + "source": [ + "It's interesting to see how our different categories are grouped in relation to one another.\n", + "\n", + "We can look at the trustworthiness score to better understand how well the structure of the original dataset was preserved by our 2D projection" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "Say1EIdqMLoq", + "outputId": "3c3f1623-8ef2-4703-9dcb-1e2ba3e148d7" + }, + "outputs": [], + "source": [ + "from sklearn.manifold import trustworthiness\n", + "trustworthiness(X_train, X_train_umap, n_neighbors=15)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "1OzvVajcm2jX" + }, + "source": [ + "It looks like this projection is a great representation of our full dataset.\n", + "\n", + "Let's now run the same thing with the accelerator turned on." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# cuml.accel already loaded at top of notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "imdGmrWvnOLP", + "outputId": "0fbbd970-743f-4b82-8c39-87b064ff2ed8" + }, + "outputs": [], + "source": [ + "import umap\n", + "umap_model = umap.UMAP(n_neighbors=15, n_components=2, random_state=42, min_dist=0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "4yUnvHUOnONH", + "outputId": "e296feec-39e8-40f7-ab8c-9e703d62a64d" + }, + "outputs": [], + "source": [ + "%%time\n", + "\n", + "# Fit UMAP model to the data\n", + "X_train_umap = umap_model.fit_transform(X_train_scaled)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 718 + }, + "id": "K1qsC-Xmndgb", + "outputId": "439aba27-ab59-4f6d-9e09-882e1b1a48d8" + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "# Plot the UMAP result\n", + "plt.figure(figsize=(10, 8))\n", + "plt.scatter(X_train_umap[:, 0], X_train_umap[:, 1], c=y_train.values.ravel(), cmap='Spectral', s=10)\n", + "plt.colorbar(label=\"Activity\")\n", + "plt.title(\"UMAP projection of the UCI HAR dataset\")\n", + "plt.xlabel(\"UMAP Component 1\")\n", + "plt.ylabel(\"UMAP Component 2\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "AyFnouZInffh" + }, + "source": [ + "Note that while the projection here is not identical to the umap-learn plot, the quality of the results are equivalent. We can run the trustworthiness score again to compare and verify this claim." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "czIK4d6WndkA", + "outputId": "b45f3ae2-e39f-4ba3-d7f1-cb4001ab6f12" + }, + "outputs": [], + "source": [ + "from sklearn.manifold import trustworthiness\n", + "trustworthiness(X_train, X_train_umap, n_neighbors=15)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "nHVmnyGinOO3" + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "JEAEGvqHztTI" + }, + "source": [ + "For more information on getting started with `cuml.accel`, check out [RAPIDS.ai](https://rapids.ai/cuml-accel/) or the [cuML Docs](https://docs.rapids.ai/api/cuml/stable/)." + ] + } + ], + "metadata": { + "accelerator": "GPU", + "colab": { + "gpuType": "T4", + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/docs-samples/data-science/gpu-accelerated-samples/pyod-gpu-anomaly-detection-demo.ipynb b/docs-samples/data-science/gpu-accelerated-samples/pyod-gpu-anomaly-detection-demo.ipynb new file mode 100644 index 00000000..728fa4b0 --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/pyod-gpu-anomaly-detection-demo.ipynb @@ -0,0 +1,917 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "373f4311", + "metadata": {}, + "source": [ + "# PyOD with `cuml.accel` (GPU-backed scikit-learn)\n", + "\n", + "Many PyOD detectors wrap **scikit-learn** estimators (for example `IsolationForest`, `PCA`, `NearestNeighbors`, `KMeans`). The RAPIDS library provides **`cuml.accel`**, which swaps supported scikit-learn estimators for **cuML** implementations on the GPU where available, with automatic CPU fallback otherwise.\n", + "\n", + "**How it applies to PyOD**\n", + "\n", + "- PyOD imports scikit-learn inside each model module (for example `pyod.models.iforest` does `from sklearn.ensemble import IsolationForest` at module load time), and that name is looked up again from the module's own globals every time `.fit()`/`.__init__()` runs.\n", + "- Because of that, **enabling `cuml.accel` after PyOD's modules are already imported does nothing by itself** -- `cuml.accel` patches `sys.modules[\"sklearn.xxx\"]` in place, but the already-imported PyOD module still holds its own old, un-patched reference until you `importlib.reload()` it. Skipping this step is a *silent* failure: PyOD keeps working, just entirely on CPU, and a naive \"CPU time vs GPU time\" comparison will show ~1x speedup for every detector and looks like a benchmarking result rather than a bug. This notebook's own \"Enable `cuml.accel`...\" section below reloads every PyOD model module used above for exactly this reason -- do not skip it or re-order the cells.\n", + "- Alternatively, enable acceleration **before** any `sklearn` / PyOD import: `cuml.accel.install()`, the `CUML_ACCEL_ENABLED=1` environment variable, or `python -m cuml.accel your_script.py`. That avoids the reload step entirely, at the cost of not being able to show a CPU baseline first in the same session.\n", + "\n", + "For each of five detectors (IForest, PCA, KNN, CBLOF, HDBSCAN), we time **`fit`** and **`predict`** on held-out test data **separately** with normal PyOD (CPU sklearn), **then** (after `cuml.accel.install()` and reload) repeat the same workflow.\n", + "\n", + "**Requirements**\n", + "\n", + "- NVIDIA GPU, CUDA, and a RAPIDS/cuML environment ([RAPIDS install](https://docs.rapids.ai/install)).\n", + "- Compatible `pyod` and scikit-learn versions ([cuml.accel limitations](https://docs.rapids.ai/api/cuml/stable/cuml-accel/limitations/)).\n", + "\n", + "**Note:** Speedups depend on dataset size, GPU, and whether each estimator is fully accelerated or falls back to CPU. Small arrays can look slower on GPU because of transfer overhead. This demo uses 20K samples to keep runtime short; for meaningful GPU speedups, increase `N_TRAIN` to 100K+ (some algorithms like HDBSCAN may require significant CPU time at that scale).\n", + "\n", + "**Known CPU-only detectors (verified against cuml 26.06 on an NVIDIA GB10, not just read off the docs):**\n", + "\n", + "- `IForest` wraps `sklearn.ensemble.IsolationForest`, which is not in cuML's accelerated estimator list (only `RandomForestClassifier`/`RandomForestRegressor` are accelerated under `sklearn.ensemble`). Its \"GPU\" row is expected to match its CPU row.\n", + "- `HDBSCAN` (the PyOD detector) calls `sklearn.cluster.HDBSCAN.fit()` for the actual clustering, and a separate `sklearn.neighbors.NearestNeighbors` to score new points at predict time. Only the `NearestNeighbors` part is intercepted by `cuml.accel` -- `cuml.accel.profile()` never lists `HDBSCAN.fit` as a potentially-accelerated call at all, so the dominant cost (fitting) still runs on CPU. Its \"GPU\" fit time is expected to match its CPU fit time; only `predict` benefits.\n", + "\n", + "`PCA`, `KNN` (`NearestNeighbors`), and `CBLOF` (`KMeans`) *are* accelerated end to end (confirmed both via `cuml.accel`'s own GPU-call logging and by checking `cuml.accel.is_proxy()` on each detector's core fitted estimator) and show real speedups once the reload step below has run." + ] + }, + { + "cell_type": "markdown", + "id": "2a712762", + "metadata": {}, + "source": [ + "## Setup: synthetic data and timing helper" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "a3e85c38", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:26:48.411618Z", + "iopub.status.busy": "2026-07-30T01:26:48.411448Z", + "iopub.status.idle": "2026-07-30T01:26:49.140475Z", + "shell.execute_reply": "2026-07-30T01:26:49.139964Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "X_train: (20000, 32) float32 X_test: (5000, 32)\n" + ] + } + ], + "source": [ + "import time\n", + "import warnings\n", + "\n", + "import numpy as np\n", + "\n", + "warnings.filterwarnings(\"ignore\")\n", + "\n", + "from pyod.utils.data import generate_data\n", + "\n", + "RANDOM_STATE = 42\n", + "# N_TRAIN = 1_000_000\n", + "# N_TEST = 200_000\n", + "N_TRAIN = 20_000\n", + "N_TEST = 5_000\n", + "N_FEATURES = 32\n", + "CONTAMINATION = 0.1\n", + "\n", + "X_train, X_test, y_train, y_test = generate_data(\n", + " n_train=N_TRAIN,\n", + " n_test=N_TEST,\n", + " n_features=N_FEATURES,\n", + " contamination=CONTAMINATION,\n", + " random_state=RANDOM_STATE,\n", + " behaviour=\"new\",\n", + ")\n", + "X_train = np.asarray(X_train, dtype=np.float32)\n", + "X_test = np.asarray(X_test, dtype=np.float32)\n", + "\n", + "print(\"X_train:\", X_train.shape, X_train.dtype, \" X_test:\", X_test.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "1f8d5116", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:26:49.142159Z", + "iopub.status.busy": "2026-07-30T01:26:49.141975Z", + "iopub.status.idle": "2026-07-30T01:26:49.145005Z", + "shell.execute_reply": "2026-07-30T01:26:49.144503Z" + } + }, + "outputs": [], + "source": [ + "def _used_gpu_accel(detector) -> bool:\n", + " \"\"\"True if any fitted sub-estimator attribute is a cuml.accel GPU proxy.\n", + "\n", + " PyOD detectors stash their wrapped sklearn estimator under different\n", + " attribute names (detector_, neigh_, clustering_estimator_, ...), so this\n", + " checks all fitted attributes rather than hard-coding one name.\n", + " \"\"\"\n", + " from cuml.accel import is_proxy\n", + "\n", + " return any(is_proxy(v) for v in vars(detector).values())\n", + "\n", + "\n", + "def time_fit_and_predict(detector, X_tr, X_te, require_gpu=False):\n", + " \"\"\"Wall times (seconds) for fit on train, then predict on test.\n", + "\n", + " If require_gpu=True, asserts that fit() actually ran through a\n", + " cuml.accel GPU proxy. This exists to catch the case where cuml.accel\n", + " was enabled *after* PyOD's model modules were already imported --\n", + " PyOD then keeps calling the original CPU sklearn class silently, and\n", + " without this check that looks like a (bad) benchmark result instead\n", + " of a bug. See the \"How it applies to PyOD\" note above.\n", + " \"\"\"\n", + " t0 = time.perf_counter()\n", + " detector.fit(X_tr)\n", + " t_fit = time.perf_counter() - t0\n", + " if require_gpu:\n", + " assert _used_gpu_accel(detector), (\n", + " f\"{type(detector).__name__}.fit did not execute any cuml.accel \"\n", + " \"GPU calls -- did you run the 'Enable cuml.accel and reload PyOD \"\n", + " \"model modules' cells below *before* this one, in a session \"\n", + " \"where this detector's PyOD module had already been imported?\"\n", + " )\n", + " t0 = time.perf_counter()\n", + " detector.predict(X_te)\n", + " t_predict = time.perf_counter() - t0\n", + " return t_fit, t_predict\n", + "\n", + "\n", + "cpu_times = {}\n", + "gpu_times = {}" + ] + }, + { + "cell_type": "markdown", + "id": "9c491bbf", + "metadata": {}, + "source": [ + "---\n", + "## 1. Isolation Forest (`IForest`)\n", + "\n", + "Wraps `sklearn.ensemble.IsolationForest`." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "9d3a2423", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:26:49.146382Z", + "iopub.status.busy": "2026-07-30T01:26:49.146322Z", + "iopub.status.idle": "2026-07-30T01:26:49.854710Z", + "shell.execute_reply": "2026-07-30T01:26:49.854238Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "IForest CPU (sklearn): fit=0.5963s predict=0.0239s\n" + ] + } + ], + "source": [ + "from pyod.models.iforest import IForest\n", + "\n", + "clf = IForest(n_estimators=200, random_state=RANDOM_STATE, n_jobs=-1)\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test)\n", + "cpu_times[\"IForest\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"IForest CPU (sklearn): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "markdown", + "id": "1f0fcc25", + "metadata": {}, + "source": [ + "---\n", + "## 2. PCA reconstruction-based detector (`PCA`)\n", + "\n", + "Uses `sklearn.decomposition.PCA` internally." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "00a647ba", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:26:49.856201Z", + "iopub.status.busy": "2026-07-30T01:26:49.856132Z", + "iopub.status.idle": "2026-07-30T01:26:49.897000Z", + "shell.execute_reply": "2026-07-30T01:26:49.896601Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "PCA CPU (sklearn): fit=0.0358s predict=0.0008s\n" + ] + } + ], + "source": [ + "from pyod.models.pca import PCA\n", + "\n", + "clf = PCA(n_components=16, random_state=RANDOM_STATE)\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test)\n", + "cpu_times[\"PCA\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"PCA CPU (sklearn): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "markdown", + "id": "f6a724eb", + "metadata": {}, + "source": [ + "---\n", + "## 3. kNN distances (`KNN`)\n", + "\n", + "Uses `sklearn.neighbors.NearestNeighbors` (and related neighbor search)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "31d2bf24", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:26:49.898818Z", + "iopub.status.busy": "2026-07-30T01:26:49.898636Z", + "iopub.status.idle": "2026-07-30T01:26:50.237082Z", + "shell.execute_reply": "2026-07-30T01:26:50.236327Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "KNN CPU (sklearn): fit=0.2078s predict=0.1255s\n" + ] + } + ], + "source": [ + "from pyod.models.knn import KNN\n", + "\n", + "clf = KNN(n_neighbors=20, method=\"mean\", metric=\"euclidean\", n_jobs=-1)\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test)\n", + "cpu_times[\"KNN\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"KNN CPU (sklearn): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "markdown", + "id": "2f931858", + "metadata": {}, + "source": [ + "---\n", + "## 4. CBLOF with default `KMeans` (`CBLOF`)\n", + "\n", + "Uses `sklearn.cluster.KMeans` for clustering by default." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "6834adfc", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:26:50.238385Z", + "iopub.status.busy": "2026-07-30T01:26:50.238266Z", + "iopub.status.idle": "2026-07-30T01:26:51.041889Z", + "shell.execute_reply": "2026-07-30T01:26:51.041148Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CBLOF CPU (sklearn): fit=0.7807s predict=0.0014s\n" + ] + } + ], + "source": [ + "from pyod.models.cblof import CBLOF\n", + "\n", + "clf = CBLOF(n_clusters=16, random_state=RANDOM_STATE, n_jobs=-1)\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test)\n", + "cpu_times[\"CBLOF\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"CBLOF CPU (sklearn): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "markdown", + "id": "f612da94", + "metadata": {}, + "source": [] + }, + { + "cell_type": "markdown", + "id": "217eb60f", + "metadata": {}, + "source": [ + "---\n", + "## 5. HDBSCAN (`HDBSCAN`)\n", + "\n", + "Uses `sklearn.cluster.HDBSCAN` for clustering and `sklearn.neighbors.NearestNeighbors` for scoring new points. Requires a recent scikit-learn that includes `sklearn.cluster.HDBSCAN`." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "3be8c444", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:26:51.043693Z", + "iopub.status.busy": "2026-07-30T01:26:51.043544Z", + "iopub.status.idle": "2026-07-30T01:27:00.213952Z", + "shell.execute_reply": "2026-07-30T01:27:00.213644Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "HDBSCAN CPU (sklearn): fit=9.0811s predict=0.0864s\n" + ] + } + ], + "source": [ + "from pyod.models.hdbscan import HDBSCAN\n", + "\n", + "clf = HDBSCAN(\n", + " min_cluster_size=20,\n", + " min_samples=10,\n", + " metric=\"euclidean\",\n", + " n_jobs=-1,\n", + " contamination=CONTAMINATION,\n", + ")\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test)\n", + "cpu_times[\"HDBSCAN\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"HDBSCAN CPU (sklearn): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "markdown", + "id": "126a7d77", + "metadata": {}, + "source": [ + "---\n", + "## Enable `cuml.accel` and reload PyOD model modules\n", + "\n", + "Run the next **two** cells once (install, then reload), then run the GPU timing cells below.\n", + "\n", + "The reload step is not optional here: `pyod.models.iforest`, `.pca`, `.knn`, `.cblof`, and `.hdbscan` were already imported by the CPU section above, each with its own `from sklearn.xxx import Yyy` bound at that time. `cuml.accel.install()` patches `sklearn`'s module cache going forward, but it cannot retroactively fix names other modules already bound before it ran. `importlib.reload()` re-executes those `from sklearn... import ...` lines against the now-patched `sklearn` modules, which is what actually lets the GPU-backed proxies take effect for PyOD's wrapped estimators." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "b5e99a2a", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:27:00.215224Z", + "iopub.status.busy": "2026-07-30T01:27:00.215156Z", + "iopub.status.idle": "2026-07-30T01:27:01.049221Z", + "shell.execute_reply": "2026-07-30T01:27:01.048688Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] Enabled managed memory.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] Accelerator installed.\n" + ] + } + ], + "source": [ + "# %load_ext cuml.accel # equivalent IPython-magic form\n", + "import cuml\n", + "\n", + "cuml.accel.install(log_level=\"debug\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "b4715ede", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:27:01.051653Z", + "iopub.status.busy": "2026-07-30T01:27:01.051403Z", + "iopub.status.idle": "2026-07-30T01:27:01.055425Z", + "shell.execute_reply": "2026-07-30T01:27:01.054799Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reloaded PyOD model modules against the cuml.accel-patched sklearn.\n" + ] + } + ], + "source": [ + "import importlib\n", + "\n", + "import pyod.models.cblof\n", + "import pyod.models.hdbscan\n", + "import pyod.models.iforest\n", + "import pyod.models.knn\n", + "import pyod.models.pca\n", + "\n", + "for _mod in (\n", + " pyod.models.iforest,\n", + " pyod.models.pca,\n", + " pyod.models.knn,\n", + " pyod.models.cblof,\n", + " pyod.models.hdbscan,\n", + "):\n", + " importlib.reload(_mod)\n", + "\n", + "# Re-bind the detector classes in *this* notebook's namespace too, so `IForest`,\n", + "# `PCA`, etc. below point at the freshly reloaded module objects rather than the\n", + "# stale ones captured by the `from pyod.models... import ...` lines earlier.\n", + "from pyod.models.cblof import CBLOF\n", + "from pyod.models.hdbscan import HDBSCAN\n", + "from pyod.models.iforest import IForest\n", + "from pyod.models.knn import KNN\n", + "from pyod.models.pca import PCA\n", + "\n", + "print(\"Reloaded PyOD model modules against the cuml.accel-patched sklearn.\")" + ] + }, + { + "cell_type": "markdown", + "id": "38297f3e", + "metadata": {}, + "source": [ + "---\n", + "## Same five workflows with `cuml.accel` enabled\n", + "\n", + "Hyperparameters match the CPU section. Timings are **`fit`** on `X_train` and **`predict`** on `X_test` (binary labels; PyOD may call `decision_function` internally). `require_gpu=True` is passed for the three detectors verified (via `cuml.accel.profile()`, on cuml 26.06) to actually run their *core* estimator on GPU: `PCA`, `KNN`, `CBLOF`. `IForest` and `HDBSCAN` are left off the check -- see the notes at the top and by the `HDBSCAN` cell below for why." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "406c239e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:27:01.058188Z", + "iopub.status.busy": "2026-07-30T01:27:01.058068Z", + "iopub.status.idle": "2026-07-30T01:27:01.705077Z", + "shell.execute_reply": "2026-07-30T01:27:01.704513Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "IForest GPU (cuml.accel): fit=0.6202s predict=0.0240s\n" + ] + } + ], + "source": [ + "# No require_gpu=True here: sklearn.ensemble.IsolationForest is not on cuML's\n", + "# accelerated estimator list, so this legitimately always runs on CPU -- a\n", + "# fit time close to the CPU row above is expected, not a regression.\n", + "clf = IForest(n_estimators=200, random_state=RANDOM_STATE, n_jobs=-1)\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test)\n", + "gpu_times[\"IForest\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"IForest GPU (cuml.accel): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "adaada03", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:27:01.707105Z", + "iopub.status.busy": "2026-07-30T01:27:01.706988Z", + "iopub.status.idle": "2026-07-30T01:27:02.052021Z", + "shell.execute_reply": "2026-07-30T01:27:02.051425Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `PCA.fit` ran on GPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `PCA` fitted attributes synced to CPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "PCA GPU (cuml.accel): fit=0.3408s predict=0.0012s\n" + ] + } + ], + "source": [ + "clf = PCA(n_components=16, random_state=RANDOM_STATE)\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test, require_gpu=True)\n", + "gpu_times[\"PCA\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"PCA GPU (cuml.accel): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "e448a429", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:27:02.053429Z", + "iopub.status.busy": "2026-07-30T01:27:02.053360Z", + "iopub.status.idle": "2026-07-30T01:27:02.235236Z", + "shell.execute_reply": "2026-07-30T01:27:02.234425Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `NearestNeighbors.fit` ran on GPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `NearestNeighbors` fitted attributes synced to CPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `NearestNeighbors.kneighbors` ran on GPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `NearestNeighbors.kneighbors` ran on GPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "KNN GPU (cuml.accel): fit=0.1584s predict=0.0198s\n" + ] + } + ], + "source": [ + "clf = KNN(n_neighbors=20, method=\"mean\", metric=\"euclidean\", n_jobs=-1)\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test, require_gpu=True)\n", + "gpu_times[\"KNN\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"KNN GPU (cuml.accel): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "e9088266", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:27:02.236916Z", + "iopub.status.busy": "2026-07-30T01:27:02.236847Z", + "iopub.status.idle": "2026-07-30T01:27:02.409348Z", + "shell.execute_reply": "2026-07-30T01:27:02.408950Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `KMeans.fit` ran on GPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `KMeans` fitted attributes synced to CPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `KMeans.predict` ran on GPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CBLOF GPU (cuml.accel): fit=0.1666s predict=0.0026s\n" + ] + } + ], + "source": [ + "clf = CBLOF(n_clusters=16, random_state=RANDOM_STATE, n_jobs=-1)\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test, require_gpu=True)\n", + "gpu_times[\"CBLOF\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"CBLOF GPU (cuml.accel): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "39d0461f", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:27:02.411447Z", + "iopub.status.busy": "2026-07-30T01:27:02.411371Z", + "iopub.status.idle": "2026-07-30T01:27:12.455279Z", + "shell.execute_reply": "2026-07-30T01:27:12.454705Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `NearestNeighbors.fit` ran on GPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[cuml.accel] `NearestNeighbors.kneighbors` ran on GPU\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "HDBSCAN GPU (cuml.accel): fit=10.0339s predict=0.0058s\n" + ] + } + ], + "source": [ + "clf = HDBSCAN(\n", + " min_cluster_size=20,\n", + " min_samples=10,\n", + " metric=\"euclidean\",\n", + " n_jobs=-1,\n", + " contamination=CONTAMINATION,\n", + ")\n", + "# No require_gpu=True here: as of cuml 26.06, cuml.accel.profile() shows only\n", + "# the *auxiliary* NearestNeighbors call (used to score new points) is\n", + "# intercepted for this detector -- the dominant cost, the core\n", + "# sklearn.cluster.HDBSCAN.fit() clustering itself, never shows up in the\n", + "# profiler's list of potentially-accelerated calls and always runs on CPU.\n", + "# Asserting require_gpu=True here would still pass (vars(clf) also contains\n", + "# the accelerated .tree_ NearestNeighbors), which would hide that mismatch\n", + "# rather than catch it -- exactly the kind of misleading \"success\" this\n", + "# notebook is trying to avoid.\n", + "tf, tp = time_fit_and_predict(clf, X_train, X_test)\n", + "gpu_times[\"HDBSCAN\"] = {\"fit\": tf, \"predict\": tp}\n", + "print(f\"HDBSCAN GPU (cuml.accel): fit={tf:.4f}s predict={tp:.4f}s\")" + ] + }, + { + "cell_type": "markdown", + "id": "a119280d", + "metadata": {}, + "source": [ + "## Summary table" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "7412aae0", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-30T01:27:12.457281Z", + "iopub.status.busy": "2026-07-30T01:27:12.457216Z", + "iopub.status.idle": "2026-07-30T01:27:12.464904Z", + "shell.execute_reply": "2026-07-30T01:27:12.464530Z" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
detectorcpu_fit_scpu_predict_sgpu_fit_sgpu_predict_sfit_speedup_cpu/gpupredict_speedup_cpu/gpu
0IForest0.5963500.0238590.6201990.0240350.9615460.992673
1PCA0.0357930.0008470.3408330.0012460.1050180.679955
2KNN0.2078370.1255370.1583880.0197521.3121976.355562
3CBLOF0.7806930.0014480.1665510.0026274.6874250.551300
4HDBSCAN9.0811300.08640010.0338720.0058320.90504714.815217
\n", + "
" + ], + "text/plain": [ + " detector cpu_fit_s cpu_predict_s gpu_fit_s gpu_predict_s \\\n", + "0 IForest 0.596350 0.023859 0.620199 0.024035 \n", + "1 PCA 0.035793 0.000847 0.340833 0.001246 \n", + "2 KNN 0.207837 0.125537 0.158388 0.019752 \n", + "3 CBLOF 0.780693 0.001448 0.166551 0.002627 \n", + "4 HDBSCAN 9.081130 0.086400 10.033872 0.005832 \n", + "\n", + " fit_speedup_cpu/gpu predict_speedup_cpu/gpu \n", + "0 0.961546 0.992673 \n", + "1 0.105018 0.679955 \n", + "2 1.312197 6.355562 \n", + "3 4.687425 0.551300 \n", + "4 0.905047 14.815217 " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "rows = []\n", + "for name in (\"IForest\", \"PCA\", \"KNN\", \"CBLOF\", \"HDBSCAN\"):\n", + " cf = cpu_times[name][\"fit\"]\n", + " cp = cpu_times[name][\"predict\"]\n", + " gf = gpu_times[name][\"fit\"]\n", + " gp = gpu_times[name][\"predict\"]\n", + " sf = cf / gf if gf > 0 else float(\"nan\")\n", + " sp = cp / gp if gp > 0 else float(\"nan\")\n", + " rows.append((name, cf, cp, gf, gp, sf, sp))\n", + "\n", + "try:\n", + " import pandas as pd\n", + "\n", + " display(\n", + " pd.DataFrame(\n", + " rows,\n", + " columns=[\n", + " \"detector\",\n", + " \"cpu_fit_s\",\n", + " \"cpu_predict_s\",\n", + " \"gpu_fit_s\",\n", + " \"gpu_predict_s\",\n", + " \"fit_speedup_cpu/gpu\",\n", + " \"predict_speedup_cpu/gpu\",\n", + " ],\n", + " )\n", + " )\n", + "except Exception:\n", + " hdr = (\n", + " f\"{'detector':<10} {'cpu_fit':>10} {'cpu_pred':>10} \"\n", + " f\"{'gpu_fit':>10} {'gpu_pred':>10} {'fit_sp':>8} {'pred_sp':>8}\"\n", + " )\n", + " print(hdr)\n", + " for name, cf, cp, gf, gp, sf, sp in rows:\n", + " print(\n", + " f\"{name:<10} {cf:10.4f} {cp:10.4f} {gf:10.4f} {gp:10.4f} \"\n", + " f\"{sf:8.2f}x {sp:8.2f}x\"\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "a79edd00", + "metadata": {}, + "source": [ + "## Optional: other ways to enable acceleration\n", + "\n", + "- **Programmatic (before sklearn / pyod):** `cuml.accel.install()`\n", + "- **CLI:** `python -m cuml.accel your_script.py`\n", + "- **Environment:** `export CUML_ACCEL_ENABLED=1`\n", + "\n", + "See [logging and profiling](https://docs.rapids.ai/api/cuml/stable/cuml-accel/logging-and-profiling/) to see GPU vs CPU fallback per call." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs-samples/data-science/gpu-accelerated-samples/rapids-gpu-accelerated-demo.ipynb b/docs-samples/data-science/gpu-accelerated-samples/rapids-gpu-accelerated-demo.ipynb new file mode 100644 index 00000000..654bbcc1 --- /dev/null +++ b/docs-samples/data-science/gpu-accelerated-samples/rapids-gpu-accelerated-demo.ipynb @@ -0,0 +1,396 @@ +{ + "nbformat": 4, + "nbformat_minor": 5, + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.13.0" + } + }, + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "id": "cell-00", + "source": [ + "# RAPIDS GPU-Accelerated Data Science Demo\n", + "\n", + "End-to-end demonstration of GPU acceleration using NVIDIA RAPIDS for DataFrame operations, machine learning, text embeddings, and similarity search. Each section compares GPU vs CPU performance." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "id": "cell-01", + "source": [ + "## Setup" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "id": "cell-02", + "source": [ + "import time\n", + "import os\n", + "import numpy as np\n", + "import cupy as cp\n", + "import cudf\n", + "import pandas as pd\n", + "\n", + "print(\"=\" * 60)\n", + "print(\"HARDWARE\")\n", + "print(\"=\" * 60)\n", + "os.system(\"nvidia-smi --query-gpu=name,memory.total --format=csv,noheader\")\n", + "print(f\"\\ncuDF version: {cudf.__version__}\")\n", + "print(f\"CuPy version: {cp.__version__}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "id": "cell-03", + "source": [ + "## 1. DataFrame Operations: cuDF vs pandas\n", + "\n", + "Operations on 100M rows. cuDF keeps data in GPU memory and processes all rows in parallel." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "id": "cell-04", + "source": [ + "num_rows = 100_000_000\n", + "rng = np.random.default_rng(seed=42)\n", + "\n", + "print(f\"Creating DataFrame with {num_rows:,} rows...\")\n", + "pdf = pd.DataFrame({\n", + " \"numbers\": rng.standard_normal(num_rows).astype(np.float32),\n", + " \"category\": rng.integers(0, 1000, num_rows).astype(np.int32),\n", + "})\n", + "gdf = cudf.DataFrame(pdf)\n", + "print(f\"Done. Size in memory: ~{pdf.memory_usage(deep=True).sum()/1e9:.1f} GB\")\n", + "\n", + "# Warmup GPU (first operation includes initialization overhead)\n", + "_ = gdf[\"category\"].sum()\n", + "cp.cuda.Stream.null.synchronize()\n", + "\n", + "results = {}\n", + "\n", + "# value_counts\n", + "start = time.perf_counter()\n", + "_ = pdf[\"category\"].value_counts()\n", + "cpu_time = time.perf_counter() - start\n", + "\n", + "start = time.perf_counter()\n", + "_ = gdf[\"category\"].value_counts()\n", + "cp.cuda.Stream.null.synchronize()\n", + "gpu_time = time.perf_counter() - start\n", + "results[\"value_counts\"] = (cpu_time, gpu_time)\n", + "print(f\"value_counts: CPU {cpu_time:.3f}s GPU {gpu_time:.3f}s Speedup: {cpu_time/gpu_time:.0f}x\")\n", + "\n", + "# groupby mean\n", + "start = time.perf_counter()\n", + "_ = pdf.groupby(\"category\")[\"numbers\"].mean()\n", + "cpu_time = time.perf_counter() - start\n", + "\n", + "start = time.perf_counter()\n", + "_ = gdf.groupby(\"category\")[\"numbers\"].mean()\n", + "cp.cuda.Stream.null.synchronize()\n", + "gpu_time = time.perf_counter() - start\n", + "results[\"groupby_mean\"] = (cpu_time, gpu_time)\n", + "print(f\"groupby mean: CPU {cpu_time:.3f}s GPU {gpu_time:.3f}s Speedup: {cpu_time/gpu_time:.0f}x\")\n", + "\n", + "# sort\n", + "start = time.perf_counter()\n", + "_ = pdf.sort_values(\"numbers\")\n", + "cpu_time = time.perf_counter() - start\n", + "\n", + "start = time.perf_counter()\n", + "_ = gdf.sort_values(\"numbers\")\n", + "cp.cuda.Stream.null.synchronize()\n", + "gpu_time = time.perf_counter() - start\n", + "results[\"sort\"] = (cpu_time, gpu_time)\n", + "print(f\"sort: CPU {cpu_time:.3f}s GPU {gpu_time:.3f}s Speedup: {cpu_time/gpu_time:.0f}x\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "id": "cell-05", + "source": [ + "## 2. String Operations\n", + "\n", + "String processing is where GPU shows the most extreme speedups \u2014 each character operation on each row runs in parallel." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "id": "cell-06", + "source": [ + "import gc\n", + "del pdf, gdf\n", + "gc.collect()\n", + "\n", + "num_rows = 50_000_000\n", + "print(f\"Creating string Series with {num_rows:,} rows...\")\n", + "\n", + "pd_series = pd.Series(\n", + " rng.choice(\n", + " [\"hello world\", \"RAPIDS is fast\", \"gpu acceleration\", \"data science\", \"machine learning\"],\n", + " size=num_rows,\n", + " )\n", + ")\n", + "gd_series = cudf.Series(pd_series)\n", + "\n", + "# String contains\n", + "start = time.perf_counter()\n", + "_ = pd_series.str.contains(\"fast\")\n", + "cpu_time = time.perf_counter() - start\n", + "\n", + "start = time.perf_counter()\n", + "_ = gd_series.str.contains(\"fast\")\n", + "cp.cuda.Stream.null.synchronize()\n", + "gpu_time = time.perf_counter() - start\n", + "print(f\"str.contains: CPU {cpu_time:.3f}s GPU {gpu_time:.3f}s Speedup: {cpu_time/gpu_time:.0f}x\")\n", + "\n", + "# String upper\n", + "start = time.perf_counter()\n", + "_ = pd_series.str.upper()\n", + "cpu_time = time.perf_counter() - start\n", + "\n", + "start = time.perf_counter()\n", + "_ = gd_series.str.upper()\n", + "cp.cuda.Stream.null.synchronize()\n", + "gpu_time = time.perf_counter() - start\n", + "print(f\"str.upper: CPU {cpu_time:.3f}s GPU {gpu_time:.3f}s Speedup: {cpu_time/gpu_time:.0f}x\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "id": "cell-07", + "source": [ + "## 3. Machine Learning: KMeans Clustering\n", + "\n", + "2M points, 500 clusters \u2014 GPU computes all distances in parallel." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "id": "cell-08", + "source": [ + "del pd_series, gd_series\n", + "gc.collect()\n", + "\n", + "from cuml.cluster import KMeans as cuKMeans\n", + "from sklearn.cluster import KMeans as skKMeans\n", + "\n", + "n_samples = 2_000_000\n", + "n_features = 50\n", + "n_clusters = 500\n", + "\n", + "X_cpu = rng.standard_normal((n_samples, n_features)).astype(np.float32)\n", + "X_gpu = cp.asarray(X_cpu)\n", + "\n", + "print(f\"KMeans: {n_samples:,} samples, {n_features} features, {n_clusters} clusters\")\n", + "\n", + "start = time.perf_counter()\n", + "skKMeans(n_clusters=n_clusters, max_iter=10, n_init=1, random_state=42).fit(X_cpu)\n", + "cpu_time = time.perf_counter() - start\n", + "\n", + "start = time.perf_counter()\n", + "cuKMeans(n_clusters=n_clusters, max_iter=10, n_init=1, random_state=42).fit(X_gpu)\n", + "cp.cuda.Stream.null.synchronize()\n", + "gpu_time = time.perf_counter() - start\n", + "\n", + "print(f\"KMeans: CPU {cpu_time:.3f}s GPU {gpu_time:.3f}s Speedup: {cpu_time/gpu_time:.0f}x\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "id": "cell-09", + "source": [ + "## 4. Random Forest Classification\n", + "\n", + "100 trees built in parallel on GPU." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "id": "cell-10", + "source": [ + "del X_cpu, X_gpu\n", + "gc.collect()\n", + "\n", + "from cuml.ensemble import RandomForestClassifier as cuRF\n", + "from sklearn.ensemble import RandomForestClassifier as skRF\n", + "\n", + "n_samples = 500_000\n", + "n_features = 50\n", + "X_cpu = rng.standard_normal((n_samples, n_features)).astype(np.float32)\n", + "y_cpu = rng.integers(0, 5, n_samples).astype(np.int32)\n", + "X_gpu = cudf.DataFrame(X_cpu)\n", + "y_gpu = cudf.Series(y_cpu)\n", + "\n", + "print(f\"Random Forest: {n_samples:,} samples, {n_features} features, 100 trees\")\n", + "\n", + "start = time.perf_counter()\n", + "skRF(n_estimators=100, max_depth=12, random_state=42, n_jobs=-1).fit(X_cpu, y_cpu)\n", + "cpu_time = time.perf_counter() - start\n", + "\n", + "start = time.perf_counter()\n", + "cuRF(n_estimators=100, max_depth=12, random_state=42).fit(X_gpu, y_gpu)\n", + "cp.cuda.Stream.null.synchronize()\n", + "gpu_time = time.perf_counter() - start\n", + "\n", + "print(f\"Random Forest: CPU {cpu_time:.3f}s GPU {gpu_time:.3f}s Speedup: {cpu_time/gpu_time:.0f}x\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "id": "cell-11", + "source": [ + "## 5. Text Embedding & KNN Similarity Search\n", + "\n", + "Generate text embeddings on GPU using SentenceTransformers, then find nearest neighbors using cuML." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "id": "cell-12", + "source": [ + "del X_cpu, X_gpu, y_cpu, y_gpu\n", + "gc.collect()\n", + "\n", + "import warnings, logging\n", + "warnings.filterwarnings('ignore')\n", + "logging.disable(logging.WARNING)\n", + "os.environ['TOKENIZERS_PARALLELISM'] = 'false'\n", + "\n", + "from sentence_transformers import SentenceTransformer\n", + "from cuml.neighbors import NearestNeighbors as cuNearestNeighbors\n", + "\n", + "# Generate synthetic documents\n", + "import random\n", + "random.seed(42)\n", + "categories = [\"electronics\", \"food\", \"clothing\", \"books\", \"toys\"]\n", + "adjectives_good = [\"amazing\", \"excellent\", \"fantastic\", \"great\", \"wonderful\"]\n", + "adjectives_bad = [\"terrible\", \"awful\", \"poor\", \"disappointing\", \"broken\"]\n", + "templates = [\n", + " \"This {cat} product is {adj}. Highly recommend!\",\n", + " \"I bought this {cat} item and it was {adj}.\",\n", + " \"The {cat} quality is {adj}, will buy again.\",\n", + " \"{adj} {cat} purchase, exactly what I needed.\",\n", + "]\n", + "\n", + "documents = []\n", + "for _ in range(10000):\n", + " cat = random.choice(categories)\n", + " adj = random.choice(adjectives_good + adjectives_bad)\n", + " tmpl = random.choice(templates)\n", + " documents.append(tmpl.format(cat=cat, adj=adj))\n", + "\n", + "print(f\"Generated {len(documents):,} documents\")\n", + "print(f\"Sample: {documents[0]}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "id": "cell-13", + "source": [ + "# Generate embeddings on GPU\n", + "model = SentenceTransformer('all-MiniLM-L6-v2', device='cuda')\n", + "\n", + "start = time.perf_counter()\n", + "embeddings = model.encode(documents, batch_size=256, show_progress_bar=False, convert_to_numpy=True)\n", + "gpu_time = time.perf_counter() - start\n", + "print(f\"Embedding {len(documents):,} docs on GPU: {gpu_time:.2f}s ({len(documents)/gpu_time:.0f} docs/sec)\")\n", + "print(f\"Embedding shape: {embeddings.shape}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "code", + "metadata": {}, + "id": "cell-14", + "source": [ + "# KNN search on GPU\n", + "queries = [\n", + " \"delicious food with amazing taste\",\n", + " \"terrible electronics that broke quickly\",\n", + " \"comfortable sports clothing\",\n", + "]\n", + "\n", + "query_embeddings = model.encode(queries, convert_to_numpy=True)\n", + "k = 5\n", + "\n", + "# GPU KNN\n", + "embeddings_gpu = cp.asarray(embeddings)\n", + "query_gpu = cp.asarray(query_embeddings)\n", + "\n", + "start = time.perf_counter()\n", + "knn = cuNearestNeighbors(n_neighbors=k, metric='cosine')\n", + "knn.fit(embeddings_gpu)\n", + "distances, indices = knn.kneighbors(query_gpu)\n", + "gpu_time = time.perf_counter() - start\n", + "print(f\"GPU KNN search ({len(queries)} queries, {len(documents):,} corpus): {gpu_time:.4f}s\")\n", + "\n", + "# Show results\n", + "print(\"\\n\" + \"=\" * 60)\n", + "for i, query in enumerate(queries):\n", + " print(f\"\\nQuery: \\\"{query}\\\"\")\n", + " print(f\" Top {k} matches:\")\n", + " for j in range(k):\n", + " idx = int(indices[i][j])\n", + " dist = float(distances[i][j])\n", + " print(f\" {j+1}. [{dist:.3f}] {documents[idx]}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "id": "cell-15", + "source": [ + "## Summary\n", + "\n", + "| Category | Operation | Expected Speedup |\n", + "|----------|-----------|------------------|\n", + "| **DataFrame** | value_counts, groupby, sort | 30\u2013100x |\n", + "| **String** | contains, upper | 50\u2013150x |\n", + "| **ML** | KMeans (2M points) | 20\u201350x |\n", + "| **ML** | Random Forest (500K rows) | 5\u201315x |\n", + "| **Embedding** | SentenceTransformer encode | 10\u201330x vs CPU |\n", + "| **KNN Search** | cuML NearestNeighbors | 50\u2013200x |" + ] + } + ] +} \ No newline at end of file