---
title: "Backend provenance and CUDA diagnostics"
author: "Yaoxiang Li"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Backend provenance and CUDA diagnostics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

Use `cudacellr` to run single-cell PCA and exact nearest-neighbour search on an
NVIDIA GPU, then inspect the results through ordinary R objects. This guide
shows how to prepare CUDA, run the workflow, and check which stages used the
GPU. Normalization and highly variable feature (HVG) selection remain sparse
CPU operations, so the complete workflow is deliberately hybrid.

## Prepare CUDA

The native backend is included in `cudaverse`; you do not need torch or
LibTorch. A Windows or Linux GPU machine needs an NVIDIA driver, cuBLAS 12,
and cuSOLVER 11. Follow the
[platform setup guide](https://cudaverse.github.io/cudaverse/articles/gpu-setup.html)
to install and verify these libraries. Current CUDA execution is unavailable
on macOS.

```{r install, eval=FALSE}
# install.packages("pak")
pak::pak("cudaverse/cudacellr")
library(cudacellr)

health <- cudaverse::cuda_diagnostics()
health$summary
health$next_steps
health$selected_backend  # "native" for the lightweight CUDA path
cudaverse::cuda_select_device("cuda")
```

Run the CUDA chunks below on your configured GPU machine. They are displayed
without execution when this article is built so installing or checking the R
package does not require a GPU. The small CPU reference at the end is executed
during the article build.

## Prepare a count matrix

`cudacellr` expects features in rows and cells in columns. This small named
matrix makes the example reproducible; replace it with your own count matrix
for analysis.

```{r counts}
library(cudacellr)
library(Matrix)

set.seed(42)
counts <- matrix(
  rpois(60 * 30, lambda = 2),
  nrow = 60,
  ncol = 30,
  dimnames = list(
    paste0("gene_", seq_len(60)),
    paste0("cell_", seq_len(30))
  )
)
counts[1, ] <- counts[1, ] + 1
counts <- Matrix(counts, sparse = TRUE)
```

## Normalize, select features, and run PCA and kNN

One call reuses each preprocessing stage and returns normalized expression,
feature statistics, the PCA model, and exact neighbours:

```{r cuda-workflow, eval=FALSE}
fit <- cudacell_workflow(
  counts,
  n_hvg = 20,
  n_components = 5,
  k = 5,
  batch_size = 8,
  device = "cuda"
)

fit
head(fit$pca$x)
head(fit$neighbors$index)
head(fit$neighbors$distance)

identical(rownames(fit$pca$x), colnames(counts))
identical(rownames(fit$neighbors$index), colnames(counts))
```

You can also call stages independently. For example, pass the PCA model
directly to neighbour search:

```{r cuda-stages, eval=FALSE}
pca <- cuda_cell_pca(
  counts,
  n_hvg = 20,
  n_components = 5,
  device = "cuda"
)
neighbors <- cuda_cell_neighbors(
  pca,
  k = 5,
  batch_size = 8,
  device = "cuda"
)
```

Native PCA returns ordinary R score and loading matrices while retaining a
device-side score cache. Passing the unmodified PCA result directly to kNN
allows reuse of that cache without uploading the scores again. Distance
blocks and deterministic top-k selection run on CUDA; the final neighbour
indices and distances are returned as R matrices. The public PCA model is
also materialized on the host, so this is not a promise that only the final
neighbour output ever crosses the device boundary.

## Read the provenance table

`cuda_provenance()` returns the ordered compute stages that produced an
object, including inherited preprocessing:

```{r cuda-provenance, eval=FALSE}
workflow_provenance <- cuda_provenance(fit)
workflow_provenance
attr(workflow_provenance, "compute_device")
```

| Column | Meaning |
|---|---|
| `requested_device` | An explicit or automatic request, a fixed CPU stage, or an inherited choice. |
| `device` | Where that stage computed: `"cpu"` or `"cuda"`. |
| `backend` | The implementation, such as `Matrix`, `base`, `stats`, `native`, or the optional `torch` backend. |
| `selection_reason` | Why that device was used, including fixed CPU algorithms, input transfers, or CUDA availability. |
| `fallback` | Whether an `"auto"` request used CPU because CUDA was unavailable. |
| `output_device` | Where the stage result was materialized; a CUDA calculation can return an R object on CPU. |

The aggregate `compute_device` is `"cpu"` when all recorded stages use CPU,
`"cuda"` when all use CUDA, and `"hybrid"` when they contain both. A successful
`cudacell_workflow(..., device = "cuda")` is hybrid because normalization and
HVG selection run on CPU. Returning a result to R does not itself mean the
calculation fell back to CPU.

## Native stage boundaries

These semantics describe the native backend in `cudaverse` 0.4.1:

| Function | Stage | Device/backend semantics |
|---|---|---|
| `cuda_normalize_counts()` | `normalization` | CPU `Matrix` sparse operations. |
| `cuda_hvg()` | `hvg` | CPU `Matrix` feature statistics and base-R ranking. |
| `cuda_cell_pca()` | `normalization`, `hvg` | CPU sparse preprocessing; selected features are materialized as a dense input matrix. |
| `cuda_cell_pca()` | `pca_preprocessing`, `pca_decomposition` | CUDA centring/scaling and decomposition with `native`; the public PCA model is returned to R. |
| `cuda_cell_pca()` | `pca_scores_resident` | A retained CUDA score cache supports subsequent native operations. |
| `cuda_cell_neighbors()` | inherited source stages | Preserves PCA or other input provenance. |
| `cuda_cell_neighbors()` | `knn_distance`, `knn_neighbor_selection` | CUDA distance blocks and deterministic top-k with `native`; final neighbour matrices are returned to R. |
| `cudacell_workflow()` | all stages above | Hybrid CPU preprocessing and CUDA PCA/kNN. |

The optional torch compatibility backend can have different boundaries:
device-side neighbour selection depends on its stable-sort support, and
older torch APIs may use a recorded CPU selection path. Neither torch nor
LibTorch is needed for the native workflow shown here. Inspect the actual
provenance table rather than assuming that every CUDA backend is identical.

An upstream `cudasparse` input contributes prefixed source stages. Single-cell
normalization still needs a CPU `dgCMatrix`; materializing a CUDA sparse input
therefore records an `input_materialization` transfer.

## Device requests and failures

- `"cuda"` is strict: an unavailable runtime raises
  `cudaverse_cuda_unavailable` before single-cell preprocessing starts.
- `"auto"` uses CUDA when eligible; otherwise the requested numerical stages
  run on CPU and record fallback with a reason.
- `"cpu"` selects CPU computation for a reference run or a machine without
  CUDA.
- Fixed CPU stages such as normalization have `requested_device =
  "fixed-cpu"` and `fallback = FALSE`. They are part of the workflow, not a
  failed GPU request.
- A CUDA execution error after selection is reported; the workflow does not
  retry silently on CPU.

After changing drivers or CUDA library paths, restart R and run
`cudaverse::cuda_diagnostics()` again. Its `summary` and `next_steps` fields
help distinguish missing runtime libraries from a failed device self-test.

## Plan memory and workload size

Normalization and HVG calculations retain sparse `Matrix` storage.
`cuda_cell_pca()` then selects `n_hvg` rows and materializes a dense
cell-by-feature matrix. One double-precision copy of this selected matrix
uses roughly `8 * cells * n_hvg` bytes; centring, scaling, decomposition,
and host/device copies require additional memory. Lower `n_hvg` when dense
input memory is the constraint: the input is allocated before decomposition,
so reducing `n_components` alone does not make that matrix smaller.

The PCA model retains host scores/loadings and, for native CUDA, cached device
scores. Remove unused models and call `gc()` to allow their storage to be
released. `cudaverse::cuda_memory_info("cuda")` can help inspect device memory;
whole-device usage can include other applications.

Exact kNN remains quadratic in time but uses bounded distance blocks rather
than retaining the full cell-by-cell matrix. The double distance block uses
approximately `8 * min(batch_size, cells) * cells` bytes, plus backend
workspace. Returned distances use about `8 * cells * k` bytes plus an integer
index matrix. Lower `batch_size` to reduce peak memory; it does not change
the exact neighbour definition or its original-row-index tie breaking.

This tiny tutorial demonstrates the API, not a GPU speed advantage. Small
inputs may be dominated by startup and transfer costs. Benchmark your complete
workflow, including CPU normalization/HVG and transfers, before making
performance claims about single-cell analysis.

## A small CPU reference

The same API can produce a reference result without CUDA:

```{r cpu-reference}
reference <- cudacell_workflow(
  counts,
  n_hvg = 20,
  n_components = 5,
  k = 5,
  batch_size = 8,
  device = "cpu"
)

reference
identical(rownames(reference$pca$x), colnames(counts))
identical(rownames(reference$neighbors$index), colnames(counts))
cuda_provenance(reference)
```

For a GPU comparison, rerun the CUDA workflow above on the same input and
compare numerical results within an appropriate tolerance. PCA vectors may
differ in sign, so compare reconstruction or subspaces instead of raw loading
signs. Check neighbour indices and distances, including distance ties, and
retain the runtime diagnostics and provenance with the results.

Package installation and ordinary checks do not require CUDA. Hardware
validation must actually execute the CUDA calls, verify the selected backend
and recorded stages, and check numerical agreement; a skipped GPU example or
job is not evidence of a successful hardware test.