---
title: "High-Dimensional Classification Boundaries"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{High-Dimensional Classification Boundaries}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse   = TRUE,
  comment    = "#>",
  fig.width  = 6,
  fig.height = 5
)
```

## The problem with more than two features

A boundary plot requires two axes. When a classifier is trained on more than two
features, there is no single 2D picture that faithfully represents the full
decision boundary. The boundary exists in a higher-dimensional space.

`classbound` provides two strategies for visualizing high-dimensional boundaries:

| Strategy | What it shows |
|---|---|
| **2D Slice** | Two selected features on the axes; all other features are fixed at reference values |
| **Projection** | Two axes that are linear combinations of all features; the regions are the model's predictions on that 2D plane through the data center |

Neither is universally better. Which to choose depends on your question.

## Setup

```{r setup_data, message=FALSE, warning=FALSE}
library(classbound)
library(palmerpenguins)

# Use three numeric features from palmerpenguins
penguins <- na.omit(palmerpenguins::penguins[
  ,
  c("species", "bill_length_mm", "bill_depth_mm", "flipper_length_mm")
])

# Fit a decision tree on all three features
m3 <- fit_model(
  penguins, species ~ bill_length_mm + bill_depth_mm + flipper_length_mm,
  rpart::rpart
)
m3
```

## 2D Slice

A 2D slice selects two features as the plot axes and holds all other features fixed
at their training-set median (for numeric features) or mode (for categorical features).

The slice shows how the classifier behaves at a particular cross-section of the
feature space. Two observations that appear in the same region may be separated
by the classifier in a dimension that is being held fixed.

```{r slice, message=FALSE, warning=FALSE}
# Visualize bill_length_mm vs bill_depth_mm
# flipper_length_mm is automatically fixed at its training-set median
m3_slice <- boundary_compute(
  m3,
  feature_range = list(bill_length_mm = c(30, 60), bill_depth_mm = c(10, 25)),
  resolution    = 60
)

plot_boundary(
  m3_slice,
  obs_data   = penguins,
  x_col      = "bill_length_mm",
  y_col      = "bill_depth_mm",
  true_label = "species"
)
```

The warning message tells you which features were automatically imputed, and at what
values. You can override this with the `reference` argument:

```{r slice_reference, message=FALSE, warning=FALSE}
# Fix flipper_length_mm at a specific value instead of the median
m3_slice2 <- boundary_compute(
  m3,
  feature_range = list(bill_length_mm = c(30, 60), bill_depth_mm = c(10, 25)),
  resolution    = 60,
  reference     = list(flipper_length_mm = 200)
)

plot_boundary(m3_slice2,
  obs_data = penguins,
  x_col = "bill_length_mm", y_col = "bill_depth_mm",
  true_label = "species"
)
```

The boundary shape changes as the reference value changes, because a different slice
through the 3D boundary is being rendered.

## Projection

A projection collapses all features into two new axes by multiplying the feature
matrix by a 2-column basis matrix. Rather than holding features fixed, it combines
them: each projected point summarizes information from all original features simultaneously.

The decision regions are computed differently from the points. The grid is built on the
plane spanned by the two basis columns, passing through the data center, and each grid
point is mapped back to the original feature space before prediction. The colored
regions are therefore a planar section (a zero-thickness slice) of the full decision
boundary along that plane, not a projection of it.

> **Note on Explorapp vs Scripts:** In the `explorapp()` interactive UI, Projection mode uses PCA automatically. In the programmatic API (like the script below), you provide the projection matrix yourself, allowing you to use PCA, `tourr` projections, or any other projection method of your choice.


The basis must be an orthonormal matrix (columns are unit vectors, perpendicular to
each other). PCA rotation matrices satisfy this by construction.

```{r projection, message=FALSE, warning=FALSE}
feat_cols <- c("bill_length_mm", "bill_depth_mm", "flipper_length_mm")

# Compute PCA on the three numeric features
pca <- prcomp(penguins[, feat_cols], scale. = TRUE)
basis <- pca$rotation[, 1:2] # first two principal components

# Project the training data manually to get axis ranges
x_mat <- scale(penguins[, feat_cols], center = pca$center, scale = pca$scale)
z_mat <- x_mat %*% basis

# Compute boundary in projected space, inverse-project for prediction
m3_proj <- boundary_compute(
  m3,
  feature_range = list(
    PC1 = range(z_mat[, 1]) + c(-0.5, 0.5),
    PC2 = range(z_mat[, 2]) + c(-0.5, 0.5)
  ),
  resolution = 60,
  projection = list(basis = basis, center = pca$center, scale = pca$scale)
)

plot_boundary(
  m3_proj,
  obs_data   = penguins,
  x_col      = "PC1",
  y_col      = "PC2",
  true_label = "species"
)
```

Points are forward-projected onto the plane for overlay. Points that are far from the
projection plane in the original 3D space are rendered with lower opacity (depth fading),
giving a visual indication of how faithfully each point's position is captured. A point
far from the plane can appear inside a region that does not match the model's prediction
for that point, because the region color describes the plane, not the point itself.

## When to use each strategy

**2D Slice** is the right choice when:
- You want to examine how a specific pair of features interacts
- Other dimensions are genuinely unimportant or you want to hold them constant
- You want to compare the boundary at different values of a third variable

**Projection** is the right choice when:
- All features contribute to the boundary
- You want a global overview of the multi-dimensional structure
- You are using `tourr` to animate through different projection angles

## High-dimensional data: data69_1

The package includes the UCI Waveform dataset (`data69_1`) with 21 numeric features
and 3 classes. This dataset is better suited than `palmerpenguins` when genuinely
high-dimensional boundary exploration is the goal.

```{r data69, eval=FALSE}
data(data69_1)
dim(data69_1) # 5000 x 22

# Fit on first 5 features for a manageable example
d <- data69_1[1:500, c("Y", "V1", "V2", "V3", "V4", "V5")]
d$Y <- as.factor(d$Y)
m_nd <- fit_model(d, Y ~ ., rpart::rpart)

# 2D slice: V1 vs V2, others at median
m_nd_slice <- boundary_compute(m_nd,
  feature_range = list(V1 = c(-3, 3), V2 = c(-3, 3)),
  resolution    = 50
)

plot_boundary(m_nd_slice,
  obs_data = d,
  x_col = "V1", y_col = "V2", true_label = "Y"
)
```

For the full high-dimensional tourr workflow using `data69_1`, see
`vignette("tourr-workflow")`.