===============
LIBRARY RULES
===============
From library maintainers:
- Prefer whole-tensor operations rather than element-by-element loops

### Install stablebear Package

Source: https://github.com/kthtda/stablebear/blob/main/README.md

Install the stablebear package using pip. Pre-built wheels with CUDA support are available.

```bash
pip install stablebear
```

--------------------------------

### Install StableBear Package

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/stablebear_intro_mnist_vis.ipynb

Installs the stablebear package. Use 'stablebear-cpu' for CPU-only or 'stablebear' if CUDA GPUs are available.

```python
# Install stablebear-cpu or, alternatively, stablebear if one or more CUDA GPUs are available
!pip install stablebear-cpu
```

--------------------------------

### Install stablebear using pip

Source: https://github.com/kthtda/stablebear/blob/main/docs/installing.md

Use this command to install the latest stable version of stablebear directly from PyPI.

```bash
pip install stablebear
```

--------------------------------

### Install stablebear from source using pip

Source: https://github.com/kthtda/stablebear/blob/main/docs/installing.md

Install stablebear from source using pip, which is useful if your CPU is older than x86-64-v3 or if you want a build tuned to your exact CPU.

```bash
pip install --no-binary=stablebear stablebear
```

--------------------------------

### Install CUDA Toolkit

Source: https://github.com/kthtda/stablebear/blob/main/docs/gpu.md

Install the CUDA runtime via pip if a system-wide CUDA installation is not present. This is required for GPU acceleration on supported platforms.

```default
pip install cuda-toolkit[cudart]
```

--------------------------------

### Visualize Point Cloud Data

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/pytorch_tda_classifier.ipynb

Visualizes one example point cloud from each class (disk and annulus) using matplotlib. Ensure matplotlib is installed.

```python
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(8, 3.5))
for ax, idx, title in zip(axes, [0, n_samples // 2], ["Disk (class 0)", "Annulus (class 1)"]):
    pts = pclouds[idx].to_numpy() if hasattr(pclouds[idx], 'to_numpy') else np.array(pclouds[idx])
    ax.scatter(pts[:, 0], pts[:, 1], s=10, alpha=0.7)
    ax.set_aspect("equal")
    ax.set_title(title)
plt.tight_layout()
plt.show()
```

--------------------------------

### Example: Combining PCF Generation and Reduction

Source: https://github.com/kthtda/stablebear/blob/main/docs/evaluation.md

Creates noisy sine and cosine functions, computes their means, and plots the results. This example demonstrates the practical application of PCF generation and reduction operations.

```python
import stablebear as sb
from stablebear.random import noisy_sin, noisy_cos
from stablebear.plotting import plot as plotpcf
import matplotlib.pyplot as plt


def plot_combining_example(sin_color="b", cos_color="r"):
    M = 10
    A = sb.zeros((2, M))

    A[0, :] = noisy_sin((M,), n_points=100)
    A[1, :] = noisy_cos((M,), n_points=15)

    fig, ax = plt.subplots(figsize=(6, 2))

    # Plot individual noisy functions
    plotpcf(A[0, :], ax=ax, color=sin_color, linewidth=0.5, alpha=0.4)
    plotpcf(A[1, :], ax=ax, color=cos_color, linewidth=0.5, alpha=0.4)

    # Compute and plot means
    Aavg = sb.mean(A, dim=1)
    plotpcf(Aavg[0], ax=ax, color=sin_color, linewidth=2, label="sin")
    plotpcf(Aavg[1], ax=ax, color=cos_color, linewidth=2, label="cos")

    ax.set_xlabel("t")
    ax.set_ylabel("f(t)")
    ax.legend()
    fig.tight_layout()
    return fig
```

--------------------------------

### Install Prerequisites for StableBear

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/stablebear_intro_mnist_vis.ipynb

Installs necessary Python packages for plotting, data handling, and barcode computation, which are required for the stablebear library.

```python
# For better plotting support in notebook
!pip install ipympl
!pip install tqdm
# For creation of sample dataset
!pip install scikit-learn
!pip install pandas
# For barcode computation
!pip install ripser
```

--------------------------------

### Complete Persistent Homology Example

Source: https://github.com/kthtda/stablebear/blob/main/docs/persistence.md

An end-to-end example demonstrating the creation of point clouds, computation of persistent homology, conversion to stable ranks, and visualization. Requires stablebear, numpy, and matplotlib.pyplot.

```default
import stablebear as sb
from stablebear import persistence
from stablebear.plotting import plot as plotpcf
import numpy as np
import matplotlib.pyplot as plt

shape = (10, 20)
pcloud_dim = 4

# Create and fill point cloud tensor
pclouds = sb.zeros(shape, dtype=sb.pcloud64)
for i in range(shape[0]):
    for j in range(shape[1]):
        n_points = np.random.randint(20, 100)
        pclouds[i, j] = np.random.randn(n_points, pcloud_dim)

# Compute persistent homology (H0 and H1)
bcs = persistence.compute_persistent_homology(pclouds, max_dim=1)
print(bcs.shape)   # (10, 20, 2)

# Convert to stable ranks
sranks = persistence.barcode_to_stable_rank(bcs)
print(sranks.shape) # (10, 20, 2)

# Plot H1 stable ranks for the first row of point clouds
plotpcf(sranks[0, :, 1])
plt.title('H1 stable ranks for pclouds[0, :, :]')
plt.show()

# Distance matrix between H1 stable ranks in the first row
D = sb.pdist(sranks[0, :, 1], verbose=False)
```

--------------------------------

### Build stablebear from a Git checkout

Source: https://github.com/kthtda/stablebear/blob/main/docs/installing.md

Clone the stablebear repository and install it from source. This allows for building specific branches or tags and fine-tuning the build for your CPU.

```bash
git clone https://github.com/kthtda/stablebear.git
cd stablebear

# optional: select a specific tagged version to build
git checkout tags/v0.4.1

pip install .
```

--------------------------------

### Full Mode Rectangular Matrix Example

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/cuda_engine.md

Shows block decomposition for a rectangular matrix in Full mode, where all blocks are computed. Illustrates block sizes based on row and column bands.

```text
Full mode (cdist), 2 row functions x 3 column functions, block side length 2:

         col 0-1   col 2
        +---------+------+
row 0-1 |    A    |   B  |
        +---------+------+
```

--------------------------------

### Using RandomGenerator with Poisson Distribution

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/deterministic_random.md

An example of a consumer function, `sample_poisson`, utilizing `parallel_walk` with a `RandomGenerator` to draw samples from a Poisson distribution.

```cpp
// Example: sample_poisson uses parallel_walk with a generator
sb::parallel_walk(out, gen,
  [&](const std::vector<size_t>& idx, auto& engine) {
    std::poisson_distribution<size_t> countDist(lambda);
    auto nPoints = countDist(engine);
    // ... fill point cloud using engine ...
  }, exec);
```

--------------------------------

### Plotting TDA Pipeline Example

Source: https://github.com/kthtda/stablebear/blob/main/docs/persistence.md

Demonstrates a full TDA pipeline: generating a point cloud, computing persistent homology, converting to stable ranks, and plotting the results. Requires matplotlib and persim.

```python
import numpy as np
import matplotlib.pyplot as plt
from stablebear import persistence
from stablebear.plotting import plot_barcode, plotpcf
import persim

def plot_tda_pipeline(h0_color="steelblue", h1_color="orangered"):
    # 1. Noisy circle (clear H1 topology)
    rng = np.random.RandomState(10)
    theta = rng.uniform(0, 2 * np.pi, 30)
    r = 1.0 + rng.normal(0, 0.15, 30)
    points = np.column_stack([r * np.cos(theta), r * np.sin(theta)]).astype(np.float64)

    # 2. Compute persistent homology
    bcs = persistence.compute_persistent_homology(points, max_dim=1, verbose=False)
    bc_h0, bc_h1 = bcs[0], bcs[1]

    # 3. Convert to stable rank
    sranks = persistence.barcode_to_stable_rank(bcs)

    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 3),
                                         gridspec_kw={"width_ratios": [1, 1, 1.2]})

    # Left: point cloud
    ax1.scatter(points[:, 0], points[:, 1], s=15, color="grey", edgecolors="black",
                linewidths=0.5)
    ax1.set_aspect("equal")
    ax1.set_title("Point cloud")

    # Middle: persistence diagram (via persim)
    import persim
    persim.plot_diagrams(
        [np.asarray(bc_h0), np.asarray(bc_h1)],
        ax=ax2, legend=True, show=False,
    )
    legend = ax2.get_legend()
    legend.get_frame().set_alpha(0)
    fg = ax2.xaxis.label.get_color()
    for text in legend.get_texts():
        text.set_color(fg)
    for line in ax2.get_lines():
        line.set_color(fg)
    ax2.set_title("Persistence diagram")

    # Right: stable rank
    plotpcf(sranks[0], ax=ax3, max_time=2, color=h0_color, linewidth=2,
            label="H0")
    plotpcf(sranks[1], ax=ax3, max_time=2, color=h1_color, linewidth=2,
            label="H1")
    ax3.set_xlabel("t")
    ax3.set_ylabel("rank")
    ax3.set_title("Stable rank")
    leg3 = ax3.legend(fontsize=8)
    leg3.get_frame().set_alpha(0)
    for text in leg3.get_texts():
        text.set_color(fg)

    fig.tight_layout(w_pad=1.5)
    return fig

```

--------------------------------

### Check GPU Availability and Count

Source: https://github.com/kthtda/stablebear/blob/main/docs/quickstart.md

Use `stablebear.gpu` functions to check if an NVIDIA GPU is available (`gpu.has_nvidia_gpu`) and to get the number of available GPUs (`gpu.nvidia_gpu_count`).

```python
from stablebear import gpu

gpu.has_nvidia_gpu()             # True/False
gpu.nvidia_gpu_count()           # number of available GPUs

```

--------------------------------

### Plot Accumulated Persistence Functions (APFs) Example

Source: https://github.com/kthtda/stablebear/blob/main/docs/persistence.md

This function generates a noisy circle point cloud, computes its persistent homology, converts barcodes to APFs, and visualizes the persistence barcode and the resulting APFs.

```python
def plot_apf_example(h0_color="steelblue", h1_color="orangered"):
    from stablebear import persistence
    from stablebear.plotting import plot_barcode

    # Noisy circle
    rng = np.random.RandomState(10)
    theta = rng.uniform(0, 2 * np.pi, 30)
    r = 1.0 + rng.normal(0, 0.15, 30)
    points = np.column_stack([r * np.cos(theta), r * np.sin(theta)]).astype(np.float64)

    bcs = persistence.compute_persistent_homology(points, max_dim=1, verbose=False)
    bc_h0, bc_h1 = bcs[0], bcs[1]
    apfs = persistence.barcode_to_accumulated_persistence(bcs, verbose=False)

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3))

    # Left: barcode
    y = plot_barcode(bc_h0, ax=ax1, color=h0_color, linewidth=2, label="H0")
    plot_barcode(bc_h1, ax=ax1, y_offset=y + 1, color=h1_color, linewidth=2,
                 label="H1")
    ax1.set_yticks([])
    ax1.set_xlabel("t")
    ax1.set_title("Persistence barcode")
    ax1.legend(fontsize=8)

    # Right: APF
    plotpcf(apfs[0], ax=ax2, max_time=2, color=h0_color, linewidth=2, label="H0")
    plotpcf(apfs[1], ax=ax2, max_time=2, color=h1_color, linewidth=2, label="H1")
    ax2.set_xlabel("m")
    ax2.set_ylabel("APF(m)")
    ax2.set_title("Accumulated persistence function")
    ax2.legend(fontsize=8)

    fig.tight_layout(w_pad=1.5)
    return fig
```

--------------------------------

### Plot PCF Arithmetic Results

Source: https://github.com/kthtda/stablebear/blob/main/docs/plotting.md

Visualize the result of arithmetic operations between PCFs. This example shows the addition of two PCFs and plots them alongside their individual components.

```python
def plot_arithmetic():
    f = sb.Pcf(np.array([[0, 1], [1, 3], [3, 1]], dtype=np.float32))
    g = sb.Pcf(np.array([[0, 2], [2, 0]], dtype=np.float32))

    fig, axes = plt.subplots(1, 3, figsize=(9, 2.5), sharex=True, sharey=True)
    for ax, pcf, title in [
        (axes[0], f, "f"),
        (axes[1], g, "g"),
        (axes[2], f + g, "f + g"),
    ]:
        plotpcf(pcf, ax=ax, max_time=5, linewidth=2)
        ax.set_title(title)
        ax.set_ylim(-0.3, 6)
    axes[0].set_ylabel("value")
    fig.tight_layout()
    return fig
```

--------------------------------

### Tensor with Symmetric Matrix dtype

Source: https://github.com/kthtda/stablebear/blob/main/docs/distances.md

Example of creating a tensor with a symmetric matrix dtype, specifically symmat64.

```default
T = sb.zeros((10,), dtype=sb.symmat64)
```

--------------------------------

### LowerTriangle Block Filtering Example

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/cuda_engine.md

Illustrates how blocks are filtered in LowerTriangle mode, skipping those entirely above the diagonal. Shows diagonal, below-diagonal, and skipped blocks.

```text
LowerTriangle, 9x9, blockSide=3:

         col 0-2   col 3-5   col 6-8
        +---------+---------+---------+
row 0-2 |    A    |    -    |    -    |
        +---------+---------+---------+
row 3-5 |    B    |    C    |    -    |
        +---------+---------+---------+
row 6-8 |    D    |    E    |    F    |
        +---------+---------+---------+
```

--------------------------------

### Plot Array Evaluation Example

Source: https://github.com/kthtda/stablebear/blob/main/docs/evaluation.md

Visualizes the evaluation of each PCF in the tensor X at multiple specified time points (1, 2, and 4). Each plot shows the PCF, vertical lines at evaluation times, and markers at the evaluated points.

```python
import matplotlib.pyplot as plt
from stablebear.plotting import plotpcf

def plot_tensor_eval_array():
    times = [1, 2, 4]
    fig, axes = plt.subplots(2, 2, figsize=(6, 4), sharex=True, sharey=True)

    for ax, pcf, label in [
        (axes[0, 0], f, "X[0,0] = f"),
        (axes[0, 1], g, "X[0,1] = g"),
        (axes[1, 0], 0.5 * g, "X[1,0] = 0.5g"),
        (axes[1, 1], f, "X[1,1] = f"),
    ]:
        plotpcf(pcf, ax=ax, max_time=5, linewidth=2)
        for t in times:
            val = pcf(t)
            ax.axvline(t, color="red", linestyle="--", linewidth=1, alpha=0.5)
            ax.plot(t, val, "ro", markersize=6, zorder=5)
        ax.set_title(label, fontsize=10)
        ax.set_ylim(-0.3, 5)

    fig.suptitle(f"X({times}) — shape (2, 2, 3)", fontsize=11)
    fig.tight_layout()
    return fig
```

--------------------------------

### Creating a Barcode Object

Source: https://github.com/kthtda/stablebear/blob/main/docs/persistence.md

Demonstrates how to instantiate the Barcode class from a NumPy array of (birth, death) pairs. Requires stablebear.persistence.

```default
from stablebear.persistence import Barcode

bc = Barcode(np.array([[0.0, 1.5],
                        [0.2, 3.0],
                        [0.5, 0.8]]))
```

--------------------------------

### Data Preparation for PyTorch DataLoader

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/pytorch_tda_classifier.ipynb

Prepares training and testing data and creates a PyTorch DataLoader for batching. Ensure X and y are already defined and split.

```python
# Train/test split
n_train = 160
X_train, X_test = X[:n_train], X[n_train:]
y_train, y_test = y[:n_train], y[n_train:]

train_loader = DataLoader(
    TensorDataset(X_train, y_train),
    batch_size=32, shuffle=True
)
```

--------------------------------

### Summary of Operation Configurations

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/cuda_engine.md

A table summarizing the configurations for different operations, including task class, output type, triangle skip mode, result writer, and block mode.

```text
| Operation   | Task class                    | Output type       | TriangleSkipMode        | ResultWriter                  | BlockMode       |
|-------------|-------------------------------|-------------------|-------------------------|-------------------------------|-----------------|
| `pdist`     | `CudaPairwiseIntegrationTask` | `DistanceMatrix`  | `LowerTriangleSkipDiag` | `DistanceMatrixResultWriter`  | `LowerTriangle` |
| `l2_kernel` | `CudaPairwiseIntegrationTask` | `SymmetricMatrix` | `LowerTriangle`         | `SymmetricMatrixResultWriter` | `LowerTriangle` |
| `cdist`     | `CudaCrossIntegrationTask`    | `Tensor`          | `None`                  | `DenseResultWriter`           | `Full`          |
```

--------------------------------

### Create a tensor of distance matrices

Source: https://github.com/kthtda/stablebear/blob/main/docs/distances.md

Distance matrices can be stored in tensors using the distmat32 or distmat64 dtypes. This example shows creating a 1-D tensor holding 10 distance matrices.

```python
import stablebear as sb

# A 1-D tensor holding 10 distance matrices
T = sb.zeros((10,), dtype=sb.distmat64)

# Assign a matrix into the tensor
m = sb.DistanceMatrix(5, dtype=sb.float64)
m[0, 1] = 3.14
T[0] = m
```

--------------------------------

### Constructing Non-Numeric Tensors (PCF)

Source: https://github.com/kthtda/stablebear/blob/main/docs/concepts.md

Shows how to construct Pcf and PcfTensor from lists of elements.

```default
# Non-numeric tensors from lists of elements
f = sb.Pcf([[0, 1.0], [1, 2.0]])
g = sb.Pcf([[0, 3.0], [2, 4.0]])
T = sb.PcfTensor([f, g])
```

--------------------------------

### Tensor Slicing and Copying

Source: https://github.com/kthtda/stablebear/blob/main/docs/distances.md

Demonstrates slicing, copying, flattening, and equality comparison for tensors, similar to other tensor types.

```default
sub = T[2:5]           # slice — shape (3,)
T2 = T.copy()          # independent deep copy
flat = T.flatten()     # shape (10,)
T == T2                # True
```

--------------------------------

### Plot Tensor Evaluation Example

Source: https://github.com/kthtda/stablebear/blob/main/docs/evaluation.md

Visualizes the scalar evaluation of each PCF within the tensor X at a specific time point (t=2). It highlights the evaluation time and the resulting value on each PCF plot.

```python
import matplotlib.pyplot as plt
from stablebear.plotting import plotpcf

def plot_tensor_eval_example():
    t_eval = 2
    fig, axes = plt.subplots(2, 2, figsize=(6, 4), sharex=True, sharey=True)

    for ax, pcf, label in [
        (axes[0, 0], f, "X[0,0] = f"),
        (axes[0, 1], g, "X[0,1] = g"),
        (axes[1, 0], 0.5 * g, "X[1,0] = 0.5g"),
        (axes[1, 1], f, "X[1,1] = f"),
    ]:
        plotpcf(pcf, ax=ax, max_time=5, linewidth=2)
        val = pcf(t_eval)
        ax.axvline(t_eval, color="red", linestyle="--", linewidth=1, alpha=0.7)
        ax.plot(t_eval, val, "ro", markersize=6, zorder=5)
        ax.set_title(label, fontsize=10)
        ax.set_ylim(-0.3, 5)

    fig.suptitle(f"X({t_eval}) = {X(t_eval).tolist()}", fontsize=11)
    fig.tight_layout()
    return fig
```

--------------------------------

### Rectangle Iteration Kernel Logic

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/cuda_engine.md

Illustrates the logic of the rectangle iteration kernel for simultaneously walking two PCFs. It shows how rectangles are formed and how operation-specific functions compute contributions.

```text
f: |---3---|---1---|---0---|
g: |----2----|---0----|

Rectangles: [0,1)x(3,2)  [1,2)x(1,2)  [2,3)x(1,0)  [3,4)x(0,0)
```

--------------------------------

### Operation-Specific Accumulation

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/cuda_engine.md

Provides examples of how operation-specific functions compute contributions within the rectangle iteration kernel for different distance metrics and inner products. Results are accumulated per CUDA thread.

```cpp
// L1 distance:
sum += (r - l) * |f_v - g_v|
// Lp distance:
sum += (r - l) * |f_v - g_v|^p,  finalize: sum^(1/p)
// L2 inner product:
sum += (r - l) * f_v * g_v
```

--------------------------------

### Compare Generators with Identical Seeds

Source: https://github.com/kthtda/stablebear/blob/main/docs/random.md

Demonstrates that two generators initialized with the same seed produce identical random outputs. This highlights the deterministic nature of the random generation.

```python
gen_a = sb.random.Generator(seed=123)
gen_b = sb.random.Generator(seed=123)

X = sb.random.noisy_sin((5, 10), generator=gen_a)
Y = sb.random.noisy_sin((5, 10), generator=gen_b)
# X and Y are identical
```

--------------------------------

### Create and Inspect a PCF

Source: https://github.com/kthtda/stablebear/blob/main/docs/api_base.md

Demonstrates how to create a PCF from a NumPy array and access its size. The PCF is defined by time-value pairs.

```python
import numpy as np
import stablebear as sb
f = sb.Pcf(np.array([[0.0, 1.0], [1.0, 2.0], [3.0, 0.0]], dtype=np.float32))
f.size
```

--------------------------------

### Tensor Indexing and Slicing

Source: https://github.com/kthtda/stablebear/blob/main/docs/concepts.md

Demonstrates basic tensor indexing and slicing operations, similar to NumPy. Returns a single element (Pcf) or a tensor view.

```default
X = sb.zeros((10, 5, 4))

# Single element -- returns a Pcf
f = X[3, 2, 1]

# Slicing -- returns a tensor (view)
row = X[3, :, :]        # shape (5, 4)
sub = X[2:8, 1:, 2]     # shape (6, 4)
```

--------------------------------

### Visualize Persistence Diagram and Stable Rank PCF

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/stablebear_intro_mnist_vis.ipynb

Plots the original persistence diagram alongside its corresponding stable rank PCF. Requires matplotlib and StableBear plotting utilities.

```python
from stablebear.plotting import plot as pcfplot

fig, axs = plt.subplots(1, 2, figsize=(6,3), constrained_layout=True)

plot_diagrams(diagrams, ax=axs[0], lifetime=True)
pcfplot(h1sr, ax=axs[1], color='tab:orange')
```

--------------------------------

### Create and Assign Point Clouds

Source: https://github.com/kthtda/stablebear/blob/main/docs/persistence.md

Initializes a tensor to hold point clouds and assigns random point clouds as NumPy arrays. Point clouds can have varying numbers of points and ambient dimensions.

```default
import stablebear as sb
import numpy as np

# A tensor that will hold 5 point clouds
pclouds = sb.zeros((5,), dtype=sb.pcloud64)

# Assign random point clouds (varying number of points, 3-dimensional)
for i in range(5):
    n_points = np.random.randint(20, 100)
    pclouds[i] = np.random.randn(n_points, 3)
```

--------------------------------

### Import Core Libraries

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/stablebear_intro_mnist_vis.ipynb

Imports essential Python libraries including numpy, matplotlib for plotting, ipympl for interactive widgets, and tqdm for progress bars.

```python
import numpy as np
import matplotlib.pyplot as plt
import ipympl
%matplotlib widget
from tqdm import trange, tqdm
```

--------------------------------

### Steps to Add New Function Type

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/cuda_engine.md

Outlines the four main steps required to add support for a new piecewise function type in the CUDA engine.

```text
1. **Define the GPU data layout** – create a point struct and a `CudaOffsetDataManager` type alias with a flattener function (analogous to `CudaPcfDataManager` + `init_pcf_data`).
2. **Write the kernel** – implement the device-side iteration function and integration kernel (analogous to `cuda_pcf_iterate_rectangles` and `cuda_pcf_block_integrate` in `cuda_pcf_kernel.cuh`).
3. **Create a BlockOp** – implement the `GpuStorage`/`init_gpu_storage`/`exec_block` interface (analogous to `PcfBlockOp`).
4. **Add task classes and factory functions** – create pairwise and cross integration task classes (analogous to `CudaPairwiseIntegrationTask` / `CudaCrossIntegrationTask`), and add factory functions in `cuda_matrix_integrate_api.hpp` / `cuda_matrix_integrate.cu`.
```

--------------------------------

### Create Higher-Dimensional Point Cloud Tensor

Source: https://github.com/kthtda/stablebear/blob/main/docs/persistence.md

Demonstrates the creation of a higher-dimensional tensor for storing point clouds, such as a 10x20 grid.

```default
# A 10 x 20 grid of point clouds
pclouds = sb.zeros((10, 20), dtype=sb.pcloud64)
```

--------------------------------

### stablebear.io.load

Source: https://github.com/kthtda/stablebear/blob/main/docs/api_base.md

Loads a tensor or object from a file in StableBear's binary format. The loaded item will retain its original type and data type.

```APIDOC
## stablebear.io.load(file)

### Description
Load a tensor or object from a file in stablebear’s binary format.
The returned item will have the same type and dtype as what was saved.

### Parameters
#### Path Parameters
None

#### Query Parameters
None

#### Request Body
None

### Parameters
* **file** (*str* *or* *file-like*) – A file path or an open file object in binary read mode.

### Returns
The loaded item.

### Return type
Tensor or [Pcf](#stablebear.functional.pcf.Pcf) or [Barcode](api_persistence.md#stablebear.persistence.barcode.Barcode) or DistanceMatrix or [SymmetricMatrix](#stablebear.symmetric_matrix.SymmetricMatrix)
```

--------------------------------

### Pinned Host Buffer Allocation

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/cuda_engine.md

Demonstrates the use of the PinnedHostBuffer class, which wraps cudaMallocHost and cudaFreeHost with RAII semantics for managing pinned host memory required for asynchronous memory transfers.

```c++
// PinnedHostBuffer class wraps cudaMallocHost/cudaFreeHost with RAII semantics.
// Required for truly asynchronous cudaMemcpyAsync.
```

--------------------------------

### Constructing Numeric Tensors

Source: https://github.com/kthtda/stablebear/blob/main/docs/concepts.md

Demonstrates creating FloatTensor and IntTensor from Python lists.

```default
# Numeric tensors
X = sb.FloatTensor([1.0, 2.0, 3.0])
Y = sb.IntTensor([[1, 2], [3, 4]])
```

--------------------------------

### Create a Stablebear Random Generator

Source: https://github.com/kthtda/stablebear/blob/main/docs/random.md

Initialize a Generator with a specific seed for reproducible random streams. This is useful for controlling random outputs in simulations or experiments.

```python
import stablebear as sb

gen = sb.random.Generator(seed=42)
```

--------------------------------

### Symmetric Matrix Storage and Access

Source: https://github.com/kthtda/stablebear/blob/main/docs/distances.md

Shows how to create and use a SymmetricMatrix, which stores only the lower triangle for efficiency. Symmetric access is demonstrated.

```default
from stablebear import SymmetricMatrix
from stablebear.typing import float32

m = SymmetricMatrix(100, dtype=float32)
m[3, 7] = 2.5
assert m[7, 3] == 2.5   # symmetric access
```

--------------------------------

### Create a Non-Deterministic Generator

Source: https://github.com/kthtda/stablebear/blob/main/docs/random.md

Initialize a Generator without providing a seed. This uses a non-deterministic seed sourced from the operating system, resulting in different outputs on each run.

```python
gen = sb.random.Generator()  # non-deterministic
```

--------------------------------

### Poisson Point Process Sampling with Custom Bounds and Seeding

Source: https://github.com/kthtda/stablebear/blob/main/docs/point_processes.md

Generates 3-D point clouds with a specified rate and custom bounding box. Uses a seeded generator for reproducible results.

```python
import stablebear as sb
gen = sb.random.Generator(seed=42)

X = sample_poisson(
    (10, 20),
    dim=3,
    rate=100.0,
    lo=[0.0, 0.0, -1.0],
    hi=[1.0, 2.0, 1.0],
    generator=gen,
)
```

--------------------------------

### Basic Poisson Point Process Sampling

Source: https://github.com/kthtda/stablebear/blob/main/docs/point_processes.md

Generates point clouds from a homogeneous spatial Poisson process. Supports specifying the number of clouds, dimension, and rate. Defaults to the unit cube sampling region.

```python
from stablebear.point_process import sample_poisson

# 100 point clouds in R^2, rate 50, in the unit square
X = sample_poisson((100,), dim=2, rate=50.0)
```

--------------------------------

### stablebear.system.force_cpu

Source: https://github.com/kthtda/stablebear/blob/main/docs/api_base.md

Forces all computations to execute on the CPU, overriding potential GPU usage.

```APIDOC
## stablebear.system.force_cpu(on: bool)

### Description
Set forced execution on CPU. By default, execution may happen on either CPU or GPU (if using a GPU-enabled build of stablebear).

### Parameters
#### Path Parameters
None

#### Query Parameters
None

#### Request Body
None

### Parameters
* **on** (*bool*) – If True, force execution on CPU for all operations. If False, execution may happen on either CPU or GPU (if using a GPU-enabled build of stablebear).
```

--------------------------------

### Converting Barcode to Stable Rank

Source: https://github.com/kthtda/stablebear/blob/main/docs/persistence.md

Shows how to convert a Barcode object to a stable rank (Pcf) using the barcode_to_stable_rank function. Requires stablebear.persistence.

```default
from stablebear.persistence import barcode_to_stable_rank

sr = barcode_to_stable_rank(bc)
# sr is a Pcf
```

--------------------------------

### PyTorch Feedforward Network Training

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/pytorch_tda_classifier.ipynb

Defines a simple feedforward neural network, sets up the optimizer and loss function, and runs a training loop. Prints training accuracy every 10 epochs.

```python
# A small feedforward network
model = nn.Sequential(
    nn.Linear(50, 32),
    nn.ReLU(),
    nn.Linear(32, 2),
)

optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
loss_fn = nn.CrossEntropyLoss()

# Training loop
for epoch in range(50):
    for X_batch, y_batch in train_loader:
        loss = loss_fn(model(X_batch), y_batch)
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
    if (epoch + 1) % 10 == 0:
        with torch.no_grad():
            train_acc = (model(X_train).argmax(1) == y_train).float().mean()
            print(f"Epoch {epoch+1:3d}  loss={loss.item():.4f}  train_acc={train_acc:.0%}")
```

--------------------------------

### Generate Point Cloud Data for Classification

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/pytorch_tda_classifier.ipynb

Samples points from a filled disk (class 0) and an annulus (class 1) to create point cloud data for classification tasks. Requires numpy and stablebear.

```python
import numpy as np
import stablebear as sb
from stablebear import persistence

rng = np.random.default_rng(42)
n_samples = 200
n_points = 60

def sample_disk(n, rng):
    """Uniform samples from a filled unit disk."""
    r = np.sqrt(rng.uniform(0, 1, n))
    theta = rng.uniform(0, 2 * np.pi, n)
    return np.column_stack([r * np.cos(theta), r * np.sin(theta)]).astype(np.float32)

def sample_annulus(n, rng):
    """Uniform samples from an annulus (r in [0.8, 1.2])."""
    r = np.sqrt(rng.uniform(0.8**2, 1.2**2, n))
    theta = rng.uniform(0, 2 * np.pi, n)
    return np.column_stack([r * np.cos(theta), r * np.sin(theta)]).astype(np.float32)

pclouds = sb.zeros((n_samples,), dtype=sb.pcloud32)
labels = np.zeros(n_samples, dtype=np.int64)

for i in range(n_samples):
    if i < n_samples // 2:
        pclouds[i] = sample_disk(n_points, rng)
        labels[i] = 0
    else:
        pclouds[i] = sample_annulus(n_points, rng)
        labels[i] = 1

print(f"Point clouds: {pclouds.shape}, Labels: {labels.shape}")
```

--------------------------------

### Create a Piecewise Constant Function (PCF)

Source: https://github.com/kthtda/stablebear/blob/main/docs/quickstart.md

Import stablebear and numpy to create a PCF from a NumPy array of (time, value) pairs. This PCF represents a function that changes value only at specific time points.

```python
import stablebear as sb
import numpy as np

# A PCF that equals 1 on [0,2), 3 on [2,5), and 0 on [5,7)
f = sb.Pcf(np.array([[0, 1],
                     [2, 3],
                     [5, 0]]))

```

--------------------------------

### Sequential Walk Function

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/deterministic_random.md

Demonstrates the sequential `walk` function for iterating through a tensor and providing a deterministically seeded engine to a callback for each element.

```cpp
// Sequential: visits every element in row-major order
walk(tensor, generator, [](const std::vector<size_t>& idx, EngineT& engine) {
  // engine is seeded deterministically from idx's flat position
});
```

--------------------------------

### Evaluating a PCF Tensor

Source: https://github.com/kthtda/stablebear/blob/main/docs/concepts.md

Shows how to evaluate all PCFs within a PCF tensor at given times.

```default
X = sb.zeros((3, 4), dtype=sb.pcf32)
# ... fill X with PCFs ...

X(1.5)        # shape (3, 4) -- one value per PCF
X([0, 1, 5])  # shape (3, 4, 3) -- each PCF evaluated at 3 times
```

--------------------------------

### stablebear.system.set_device_verbose

Source: https://github.com/kthtda/stablebear/blob/main/docs/api_base.md

Enable verbose device output. In this mode, when operations that may occur on GPU are invoked, a message is logged stating whether the operation will be performed on CPU or GPU.

```APIDOC
## stablebear.system.set_device_verbose(on: bool)

### Description
Enable verbose device output. In this mode, when operations that may occur on GPU are invoked, a message is logged stating whether the operation will be performed on CPU or GPU.

### Parameters
#### Path Parameters
- **on** (bool) - Required - Enable verbose device logging
```

--------------------------------

### Seed the Global Generator

Source: https://github.com/kthtda/stablebear/blob/main/docs/random.md

Use the global `sb.random.seed()` function to seed the default generator. This is convenient when explicit generator objects are not passed to random functions.

```python
sb.random.seed(42)
A = sb.random.noisy_sin((10, 20))

sb.random.seed(42)
B = sb.random.noisy_sin((10, 20))
# A and B are identical
```

--------------------------------

### Load MNIST Dataset

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/stablebear_intro_mnist_vis.ipynb

Fetches the MNIST dataset (version 1) containing 70,000 images of handwritten digits. Images are reshaped to 28x28 pixels.

```python
from sklearn.datasets import fetch_openml
Xin, yin = fetch_openml("mnist_784", version=1, return_X_y=True, as_frame=False)
Xin = Xin.reshape(Xin.shape[0], 28, 28)
print(Xin.shape)
```

--------------------------------

### Plot Poisson Samples

Source: https://github.com/kthtda/stablebear/blob/main/docs/point_processes.md

Generates and plots multiple samples from a Poisson point process. Useful for visualizing the distribution of points.

```python
def plot_poisson_samples(color="steelblue"):
    from stablebear.point_process import sample_poisson
    import stablebear as sb

    gen = sb.random.Generator(seed=42)
    X = sample_poisson((3,), dim=2, rate=50.0, generator=gen)

    fig, axes = plt.subplots(1, 3, figsize=(9, 3))
    for i, ax in enumerate(axes):
        pc = np.asarray(X[i])
        ax.scatter(pc[:, 0], pc[:, 1], s=8, color=color, alpha=0.7)
        ax.set_xlim(0, 1)
        ax.set_ylim(0, 1)
        ax.set_aspect("equal")
        ax.set_title(f"Sample {i + 1} ({pc.shape[0]} points)")
        ax.set_xlabel("$x_1$")
        ax.set_ylabel("$x_2$")
    fig.tight_layout()
    return fig
```

--------------------------------

### Display Sample Point Clouds

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/stablebear_intro_mnist_vis.ipynb

Visualizes the first six generated point clouds using matplotlib scatter plots. Each point cloud represents a digit from the filtered MNIST dataset.

```python
fig, axs = plt.subplots(1, 6, figsize=(10,3), tight_layout=True)
for i, ax in enumerate(axs):
    ax.scatter(Xcloud[i,:,0], Xcloud[i,:,1], s=2.0)
```

--------------------------------

### Limit GPU Count

Source: https://github.com/kthtda/stablebear/blob/main/docs/gpu.md

Specify the maximum number of GPUs stablebear should utilize. Useful when you have multiple GPUs but want to restrict usage.

```python
# Use at most 1 GPU
sb.system.limit_gpus(1)
```

--------------------------------

### Enable Verbose Device Logging

Source: https://github.com/kthtda/stablebear/blob/main/docs/gpu.md

Enable detailed logging to see which device (CPU or GPU) is used for each operation. This helps in understanding execution flow and debugging.

```python
sb.system.set_device_verbose(True)

# Now operations will print whether they run on CPU or GPU
D = sb.pdist(X)
# e.g. "Running on GPU"
```

--------------------------------

### Flat Index Computation in Parallel Walk

Source: https://github.com/kthtda/stablebear/blob/main/docs/internals/deterministic_random.md

Shows the internal implementation within `parallel_walk_impl` for reconstructing the multi-dimensional index from a flat index, ensuring consistency across threads.

```cpp
// In parallel_walk_impl:
size_t rem = flat;
for (ptrdiff_t i = ndim - 1; i >= 0; --i)
{
  idx[i] = rem % shape[i];
  rem /= shape[i];
}
```

--------------------------------

### Save and Load PCF Tensors

Source: https://github.com/kthtda/stablebear/blob/main/docs/quickstart.md

Persist PCF tensors to disk using `sb.io.save` and load them back into memory with `sb.io.load`.

```python
from stablebear.io import save, load

save(X, 'my_pcfs.sb')
X_loaded = load('my_pcfs.sb')

```

--------------------------------

### NumPy Equivalent for Multiple Boolean Masks

Source: https://github.com/kthtda/stablebear/blob/main/docs/indexing.md

Demonstrates the NumPy equivalent using `np.ix_` for achieving StableBear's outer indexing behavior with multiple boolean masks.

```default
# stablebear
X[row_mask, col_mask]

# NumPy equivalent
arr[np.ix_(row_mask, col_mask)]
```

--------------------------------

### Plotting Betti Curve Pipeline

Source: https://github.com/kthtda/stablebear/blob/main/docs/plotting.md

This function demonstrates the pipeline for generating Betti curves from persistent homology data. It visualizes the point cloud, persistence barcode, and the resulting Betti curves. Requires numpy, matplotlib, and stablebear.

```python
def plot_betti_pipeline(h0_color="steelblue", h1_color="orangered"):
    from stablebear import persistence
    from stablebear.plotting import plot_barcode

    # 1. Noisy circle
    rng = np.random.RandomState(10)
    theta = rng.uniform(0, 2 * np.pi, 30)
    r = 1.0 + rng.normal(0, 0.15, 30)
    points = np.column_stack([r * np.cos(theta), r * np.sin(theta)]).astype(np.float64)

    # 2. Compute persistent homology
    bcs = persistence.compute_persistent_homology(points, max_dim=1, verbose=False)
    bc_h0, bc_h1 = bcs[0], bcs[1]

    # 3. Convert to Betti curves
    bettis = persistence.barcode_to_betti_curve(bcs, verbose=False)

    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 3),
                                         gridspec_kw={"width_ratios": [1, 1, 1.2]})

    # Left: point cloud
    ax1.scatter(points[:, 0], points[:, 1], s=15, color="grey",
                edgecolors="black", linewidths=0.5)
    ax1.set_aspect("equal")
    ax1.set_title("Point cloud")

    # Middle: persistence barcode
    y = plot_barcode(bc_h0, ax=ax2, color=h0_color, linewidth=2, label="H0")
    plot_barcode(bc_h1, ax=ax2, y_offset=y + 1, color=h1_color, linewidth=2,
                 label="H1")
    ax2.set_yticks([])
    ax2.set_xlabel("t")
    ax2.set_title("Persistence barcode")
    ax2.legend(fontsize=8)

    # Right: Betti curves
    plotpcf(bettis[0], ax=ax3, max_time=2, color=h0_color, linewidth=2,
            label="H0")
    plotpcf(bettis[1], ax=ax3, max_time=2, color=h1_color, linewidth=2,
            label="H1")
    ax3.set_xlabel("t")
    ax3.set_ylabel("count")
    ax3.set_title("Betti curve")
    ax3.legend(fontsize=8)

    fig.tight_layout(w_pad=1.5)
    return fig
```

--------------------------------

### stablebear.gpu.detect_nvidia_gpus

Source: https://github.com/kthtda/stablebear/blob/main/docs/api_base.md

Detect NVIDIA GPUs present on the system. Uses OS-level tools (lspci, sysfs, PowerShell, system_profiler). Does not require CUDA or any NVIDIA drivers/libraries.

```APIDOC
## stablebear.gpu.detect_nvidia_gpus()

### Description
Detect NVIDIA GPUs present on the system. Uses OS-level tools (lspci, sysfs, PowerShell, system_profiler). Does not require CUDA or any NVIDIA drivers/libraries.

### Returns
#### Success Response
- **list[dict]** - A list of dicts, each with a "name" key describing the GPU. An empty list means no NVIDIA GPUs were found.
```

--------------------------------

### Initialize a StableBear Tensor

Source: https://github.com/kthtda/stablebear/blob/main/docs/tutorial_notebooks/stablebear_intro_mnist_vis.ipynb

Creates a multi-dimensional tensor filled with zeros using StableBear's tensor functionality. Useful for pre-allocating space for computations.

```python
import stablebear as sb

Z = sb.zeros((10, 3, 5))
print(Z.shape)
```

--------------------------------

### Boolean Masking for Tensor Selection

Source: https://github.com/kthtda/stablebear/blob/main/docs/concepts.md

Illustrates selecting tensor elements using boolean masks. Requires importing numpy.

```default
import numpy as np

X = sb.FloatTensor(np.arange(12, dtype=np.float32).reshape(3, 4))
mask = sb.BoolTensor(np.array([True, False, True, False]))

X[:, mask]   # shape (3, 2) -- select columns where mask is True
X[X > threshold]  # flat 1D -- all elements matching the condition
```

--------------------------------

### Load Tensor from File

Source: https://github.com/kthtda/stablebear/blob/main/docs/saving.md

Loads a tensor from a file saved in StableBear's binary format. The tensor is restored with its original type and dtype.

```python
X = sb.load('my_pcfs.sb')
```

--------------------------------

### Save Tensor to File

Source: https://github.com/kthtda/stablebear/blob/main/docs/saving.md

Saves a tensor to a file using StableBear's binary format. Ensure the tensor is defined before saving.

```python
import stablebear as sb
from stablebear.random import noisy_sin

X = noisy_sin((100,), n_points=50)
sb.save(X, 'my_pcfs.sb')
```

--------------------------------

### Query GPU Availability

Source: https://github.com/kthtda/stablebear/blob/main/docs/gpu.md

Check for the presence and count of NVIDIA GPUs, and retrieve detailed information about them. This helps in understanding the hardware capabilities for acceleration.

```python
from stablebear import gpu

gpu.has_nvidia_gpu()       # True/False
gpu.nvidia_gpu_count()     # Number of GPUs
gpu.detect_nvidia_gpus()   # Detailed GPU info
```

--------------------------------

### Generate Noisy Sine and Cosine PCFs

Source: https://github.com/kthtda/stablebear/blob/main/docs/random.md

Create tensors of piecewise constant functions approximating sine and cosine waves with added Gaussian noise. Use `n_points` to control the number of breakpoints.

```python
sines = sb.random.noisy_sin((200,), n_points=100)
cosines = sb.random.noisy_cos((10, 50), n_points=30)
```

--------------------------------

### Evaluating a Single PCF

Source: https://github.com/kthtda/stablebear/blob/main/docs/concepts.md

Demonstrates evaluating a single Pcf object at a specific time or multiple times.

```default
f = sb.Pcf([[0, 1], [2, 3], [5, 0]])

f(1.0)    # 1.0  (on the interval [0, 2))
f(3.5)    # 3.0  (on the interval [2, 5))
```

--------------------------------

### Creating Zero Tensors

Source: https://github.com/kthtda/stablebear/blob/main/docs/concepts.md

Create tensors of PCFs or scalars initialized to zero using sb.zeros(). Specify shape and dtype.

```python
import stablebear as sb

# A 1-D tensor of 100 "zero" PCFs (32-bit, the default)
X = sb.zeros((100,))

# A 2-D tensor of 64-bit PCFs
Y = sb.zeros((10, 50), dtype=sb.pcf64)

# A tensor of scalar floats
Z = sb.zeros((5, 5), dtype=sb.float32)
```

--------------------------------

### stablebear.gpu.has_nvidia_gpu

Source: https://github.com/kthtda/stablebear/blob/main/docs/api_base.md

Check whether the system has at least one NVIDIA GPU.

```APIDOC
## stablebear.gpu.has_nvidia_gpu()

### Description
Check whether the system has at least one NVIDIA GPU.

### Returns
#### Success Response
- **bool** - `True` if at least one NVIDIA GPU is detected.
```

--------------------------------

### Generate Random PCF Data

Source: https://github.com/kthtda/stablebear/blob/main/docs/quickstart.md

Use `noisy_sin` and `noisy_cos` from `stablebear.random` to generate collections of random PCFs for experimentation.

```python
from stablebear.random import noisy_sin, noisy_cos

sines   = noisy_sin((200,), n_points=100)   # 200 noisy sin functions
cosines = noisy_cos((10, 50), n_points=30)  # 10 x 50 noisy cosines

```