---
title: "Vector Field Analysis and Visualization"
author: "Zaoqu Liu"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
    fig_width: 8
    fig_height: 6
vignette: >
  %\VignetteIndexEntry{Vector Field Analysis and Visualization}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  message = FALSE,
  warning = FALSE
)
```

## Introduction

One of the most powerful features of CellODE is its ability to infer and visualize **cellular velocity fields**. The vector field represents the instantaneous direction and magnitude of cellular state changes, providing insights into:

- Differentiation trajectories
- Cell fate decisions
- Developmental dynamics

This vignette demonstrates how to use CellODE's built-in functions for vector field analysis and visualization.

## Quick Start: Using CellODE Functions

```{r load-package}
library(CellODE)
library(ggplot2)
library(Seurat)
```

### Create Example Data

First, let's create a Seurat object with simulated trajectory data:

```{r create-data}
set.seed(42)
n_cells <- 500
n_genes <- 100

# Simulate trajectory
t_true <- seq(0, 1, length.out = n_cells)

# Create expression data (genes change along trajectory)
counts <- matrix(0, nrow = n_genes, ncol = n_cells)
for (i in 1:n_genes) {
    phase <- runif(1, 0, 2 * pi)
    amplitude <- runif(1, 3, 10)
    counts[i, ] <- rpois(n_cells, lambda = pmax(1, amplitude * (1 + sin(t_true * 2 * pi + phase))))
}
rownames(counts) <- paste0("Gene", seq_len(n_genes))
colnames(counts) <- paste0("Cell", seq_len(n_cells))

# Create Seurat object
srt <- CreateSeuratObject(counts = counts)

# Add UMAP embedding (curved trajectory)
umap1 <- t_true * 10 + rnorm(n_cells, sd = 0.5)
umap2 <- 4 * sin(t_true * 2 * pi) + rnorm(n_cells, sd = 0.3)

umap_embeddings <- matrix(c(umap1, umap2), ncol = 2)
rownames(umap_embeddings) <- colnames(srt)
colnames(umap_embeddings) <- c("UMAP_1", "UMAP_2")

srt[["umap"]] <- CreateDimReducObject(
    embeddings = umap_embeddings,
    key = "UMAP_",
    assay = "RNA"
)

# Add pseudotime
srt$pseudotime <- t_true + rnorm(n_cells, sd = 0.02)
srt$pseudotime <- pmax(0, pmin(1, srt$pseudotime))
```

### Simulate CellODE Model Output

In a real workflow, you would train a CellODE model. Here we simulate the output:

```{r simulate-model}
# Simulate latent space (normally from model output)
n_latent <- 5
zs <- matrix(rnorm(n_cells * n_latent), ncol = n_latent)

# Make latent space correlate with trajectory
zs[, 1] <- t_true * 2 + rnorm(n_cells, sd = 0.2)
zs[, 2] <- sin(t_true * pi) + rnorm(n_cells, sd = 0.1)
rownames(zs) <- colnames(srt)

# Simulate vector field in latent space
vf <- matrix(0, nrow = n_cells, ncol = n_latent)
vf[, 1] <- 1 + rnorm(n_cells, sd = 0.1)  # Constant velocity in z1
vf[, 2] <- cos(t_true * pi) * pi / 2 + rnorm(n_cells, sd = 0.1)  # Derivative of sin
rownames(vf) <- colnames(srt)

# Store in Seurat object (as CellODE does)
srt@misc$X_zs <- zs
srt@misc$X_VF <- vf
```

## Vector Field Computation with CellODE

### Step 1: Compute Cosine Similarity

Use `cosine_similarity()` to compute transition probabilities in latent space:

```{r cosine-sim}
# Compute cosine similarity matrix
T_mat <- cosine_similarity(
    zs = zs,
    vf = vf,
    n_neigh = 20
)

cat("Transition matrix dimensions:", dim(T_mat), "\n")
cat("Non-zero entries:", sum(T_mat != 0), "\n")
```

### Step 2: Compute Embedding Velocity

Use `vector_field_embedding()` to project velocities to UMAP space:

```{r embed-velocity}
# Get embedding coordinates
E <- Embeddings(srt, "umap")

# Compute velocity in embedding space
V_embed <- vector_field_embedding(
    T_mat = T_mat,
    E = E,
    scale = 10
)

cat("Embedding velocity dimensions:", dim(V_embed), "\n")
cat("Mean velocity magnitude:", mean(sqrt(rowSums(V_embed^2))), "\n")
```

### Step 3: Grid-Based Vector Field

Use `vector_field_embedding_grid()` for cleaner visualization:

```{r grid-vf}
# Compute grid-based vector field
grid_result <- vector_field_embedding_grid(
    E = E,
    V = V_embed,
    smooth = 0.5,
    stream = FALSE,
    density = 1.0
)

cat("Grid points:", nrow(grid_result$E_grid), "\n")
```

## Visualization with CellODE

### Plot Pseudotime with Direction Arrows

Use `plot_pseudotime()` to visualize pseudotime with trajectory arrows:
 
```{r plot-pseudotime, fig.width=9, fig.height=6}
# Basic pseudotime plot
p1 <- plot_pseudotime(
    srt, 
    time_key = "pseudotime",
    arrow_show = TRUE,
    arrow_density = 1.5,
    point_palette = "Spectral",
    title = "Pseudotime with Direction Arrows"
)
print(p1)
```

### Plot Vector Field - Grid Mode

Use `plot_vector_field()` with grid mode for clean visualization:

```{r plot-vf-grid, fig.width=9, fig.height=6}
# Vector field - grid mode
p2 <- plot_vector_field(
    srt,
    plot_type = "grid",
    color_by = "pseudotime",
    arrow_density = 0.8,
    arrow_color = "grey20",
    point_palette = "Spectral",
    title = "Vector Field (Grid Mode)"
)
print(p2)
```

### Plot Vector Field - Stream Mode

Stream mode provides colored arrows based on velocity magnitude:

```{r plot-vf-stream, fig.width=9, fig.height=6}
# Vector field - stream mode with colored arrows
p3 <- plot_vector_field(
    srt,
    plot_type = "stream",
    color_by = "pseudotime",
    arrow_density = 1.2,
    point_palette = "viridis",
    title = "Vector Field (Stream Mode)"
)
print(p3)
```

### Plot Vector Field - Raw Mode

Raw mode shows arrows at each cell position:

```{r plot-vf-raw, fig.width=9, fig.height=6}
# Vector field - raw mode (per-cell arrows)
p4 <- plot_vector_field(
    srt,
    plot_type = "raw",
    color_by = "pseudotime",
    arrow_density = 0.3,  # Subsample for clarity
    arrow_color = "navy",
    point_palette = "magma",
    title = "Vector Field (Raw Mode)"
)
print(p4)
```

## Customization Options

### Arrow Styling

```{r arrow-styles, fig.width=9, fig.height=6}
# Custom arrow styling
p_custom <- plot_vector_field(
    srt,
    plot_type = "grid",
    color_by = "pseudotime",
    arrow_type = "closed",  # Filled arrowheads
    arrow_color = "#E63946",  # Red arrows
    arrow_density = 0.6,
    arrow_length = 0.2,
    arrow_width = 0.8,
    arrow_angle = 30,
    point_palette = "RdYlBu",
    point_size = 0.8,
    title = "Custom Arrow Styling"
)
print(p_custom)
```

### Clean Publication Figure

```{r publication-fig, fig.width=9, fig.height=6}
# Clean figure for publication
p_pub <- plot_pseudotime(
    srt,
    time_key = "pseudotime",
    arrow_show = TRUE,
    arrow_density = 1.0,
    arrow_type = "open",
    arrow_color = "black",
    point_palette = "viridis",
    point_size = 1.2,
    show_axes = FALSE,
    title = NULL
)
print(p_pub)
```

## Mathematical Background

### From Latent Dynamics to Velocity

The Neural ODE in CellODE defines the velocity in latent space:

$$\mathbf{v}_\text{latent}(t) = \frac{d\mathbf{z}}{dt} = f_\theta(\mathbf{z}(t), t)$$

### Cosine Similarity

The `cosine_similarity()` function computes:

$$\text{cos\_sim}_{ij} = \frac{\mathbf{v}_i \cdot (\mathbf{z}_j - \mathbf{z}_i)}{||\mathbf{v}_i|| \cdot ||\mathbf{z}_j - \mathbf{z}_i||}$$

### Transition Probability

The transition probabilities are computed as:

$$T_{ij} = \text{sign}(\text{cos\_sim}_{ij}) \cdot (\exp(|\text{cos\_sim}_{ij}| \cdot s) - 1)$$

### Embedding Velocity

The `vector_field_embedding()` function computes:

$$\mathbf{v}_i^\text{embed} = \sum_j T_{ij} \cdot \frac{\mathbf{e}_j - \mathbf{e}_i}{||\mathbf{e}_j - \mathbf{e}_i||}$$

## Comparison with RNA Velocity

| Aspect | CellODE | RNA Velocity (scVelo) |
|--------|---------|----------------------|
| Input | Gene expression only | Spliced/unspliced counts |
| Model | Neural ODE | Steady-state/dynamic model |
| Time inference | Joint with dynamics | Separate |
| R package | Native R | Python (reticulate needed) |

## Best Practices

### 1. Neighbor Selection

```{r best-practice-neigh, eval=FALSE}
# More neighbors = smoother but less detailed
T_smooth <- cosine_similarity(zs, vf, n_neigh = 30)

# Fewer neighbors = more detail but noisier
T_detail <- cosine_similarity(zs, vf, n_neigh = 10)
```

### 2. Visualization Density

```{r best-practice-viz, eval=FALSE}
# Publication figure (fewer arrows)
plot_vector_field(srt, arrow_density = 0.3)

# Exploration (more arrows)
plot_vector_field(srt, arrow_density = 1.0)
```

### 3. Self-Transition

```{r best-practice-self, eval=FALSE}
# Include self-transition for smoother results
V_with_self <- vector_field_embedding(T_mat, E, self_transition = TRUE)
```

## Summary

CellODE provides a complete workflow for vector field analysis:

| Function | Purpose |
|----------|---------|
| `cosine_similarity()` | Compute transition probabilities |
| `vector_field_embedding()` | Project to embedding space |
| `vector_field_embedding_grid()` | Grid-based smoothing |
| `plot_vector_field()` | Visualize vector field |
| `plot_pseudotime()` | Visualize pseudotime with arrows |

## Session Info

```{r session}
sessionInfo()
```