DataDrivenTDAFinalPresentation.qmd

---
title: "The Shape of Digits"
subtitle: "A Bayesian Topological Data Analytic Approach to Classification of Handwritten Digits"
authors:
  - name: Thomas Reinke
    affiliation: 
      - name: Baylor University
        department: Statistical Science
        # city: Waco
        # state: TX
        # country: US
        url: https://www.baylor.edu
    # email: thomas_reinke1@baylor.edu
  - name: Theophilus A. Bediako
    affiliation: 
      - name: Baylor University
        department: Statistical Science
        # city: Waco
        # state: TX
        # country: US
        url: https://www.baylor.edu
  - name: Daniel Lim
    affiliation: 
      - name: Baylor University
        department: Statistical Science
        # city: Waco
        # state: TX
        # country: US
        url: https://www.baylor.edu
date: today
date-format: "MMMM D, YYYY"
format: 
  revealjs:
    theme: 
      - quarto-assets/baylor-theme.scss
    smaller: false
    scrollable: false
    show-slide-number: all
    toc: false
    toc-depth: 1
    preview-links: true
    slide-number: c/t
    multiplex: false
    embed-resources: true
    auto-animate: true
    #footer: "Thomas Reinke"
bibliography: references.bibtex
lightbox:
  match: auto
  effect: fade
  desc-position: bottom
  loop: true
logo: "quarto-assets/baylor.png"
license: "CC BY-NC"
copyright: 
  holder: Thomas Reinke
  year: 2025
editor: 
  markdown: 
    wrap: 72
fig-width: 15
---

```{r, setup}
#| include: false
#| message: false
library(quarto)
library(knitr)
library(tidyverse)
library(conflicted)
library(janitor)
library(ggtda)
# library(TDAvis)
library(patchwork)
library(gganimate)
library(ggforce)
library(simplextree) 
library(gifski)
library(magick)  
library(ripserr)
library(reshape2)
# remotes::install_github("maroulaslab/BayesTDA") Use this if package ‘BayesTDA’ is not available for this version of R
library(BayesTDA)
library(TDAstats)
library(mvtnorm)
library(kableExtra)
library(plotly)
library(DiagrammeR)
library(transport)
library(TDA)
library(RColorBrewer)
library(Rtsne)
library(keras)
library(furrr)
library(yardstick)
library(caret)
library(imager)
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflicts_prefer(ggtda::geom_simplicial_complex)
conflicted::conflicts_prefer(plotly::layout)
conflicts_prefer(magrittr::set_names)
knitr::opts_chunk$set(
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  cache = FALSE,
  echo = FALSE,
  tidy.opts = list(width.cutoff = 100),
  tidy = FALSE,
  fig.align = "center"
)
ggplot2::theme_set(ggplot2::theme_minimal())
ggplot2::theme_update(panel.grid.minor = ggplot2::element_blank())

#------------------------------------------------------------#
```

::: {.content-hidden}
$$
{{< include quarto-assets/_macros.tex >}}
$$
:::


```{r load_data}
#--------------------------------------------------------
#----------Load & Preprocess Data------------------------
#--------------------------------------------------------
mnist <- readRDS(file = "mnist_dataset")

train <- mnist$train
test <- mnist$test

train_images <- train$images
train_labels <- as.factor(train$labels)  

test_images <- test$images  
test_labels <- as.factor(test$labels)

train_images <- train_images / 255
test_images <- test_images / 255

train_images_list <- lapply(1:nrow(train_images), function(i) {
  matrix(train_images[i, ], nrow = 28) |> t()
})

test_images_list <- lapply(1:nrow(test_images), function(i) {
  matrix(test_images[i, ], nrow = 28) |> t()
})

plot_digit <- \(image_list = train_images_list, image_index = NULL, image_df = NULL, melted = FALSE){
  if(!melted){
    image_df <- melt(image_list[image_index])
    colnames(image_df) <- c("y", "x", "value")
  }
  ggplot(image_df, aes(x = x, y = y, fill = value)) + 
    geom_raster() +
  scale_fill_gradient(low = "white", high = "black") +
  scale_y_reverse() + 
  coord_equal() +
  theme_void() + 
  theme(legend.position = "none") 
}

# plot_digit(image_index = 8)
# paste0("Label: ", train$labels[8])

binarize_images <- function(images_list, threshold = 0.5) {
  lapply(images_list, function(mat) {
    ifelse(mat < threshold, 0, 1)
  })
}

train_images_binarized <- binarize_images(train_images_list)
test_images_binarized <- binarize_images(test_images_list)

# plot_digit(image_index = 8) + plot_digit(train_images_binarized, image_index = 8)

#------------------------------------------------------------#
```

# Contents

1. [MNIST EDA](#sec-EDA)
1. [Tradiotional ML](#sec-MLMethod)
1. [Proposed Methodology](#sec-PropMethod)
1. [Analysis](#sec-Analysis)
1. [TDA + ML](#sec-TDAML)
1. [Results/Future Work](#sec-Results)
1. [References](#sec-References)


# Exploratory Data Analysis {#sec-EDA}


## Distribution of training labels

```{r dist_labels, fig.dim=c(10,6)}
#--------------------------------------------------------
#----------Distribution of training labels---------------
#--------------------------------------------------------
#first run the qmd file
train_df <- train$images |> as.data.frame()
train_df$labels <- train$labels |> as.factor()


ggplot(train_df, aes(labels, fill = labels))+
  geom_bar()+
  labs(fill = "digit")

```

:::{.notes}
- Around 6000 digits in each class
- No class imbalance
:::

## Pixel Intensity

```{r pixel_intensity, fig.dim=c(10,6)}
#----------------------------------------------------------
#-----------Pixel intensity representation-----------------
#----------------------------------------------------------
# train_df_long <- train_df |>
#   pivot_longer(-labels, names_to  = "covariate", values_to = "intensity")
# saveRDS(train_df_long, "train_df_long.rds")
train_df_long <- readRDS("train_df_long.rds")

ggplot(train_df_long) +
  geom_histogram(
    aes(intensity),
    bins = 30,
    fill = "#2C7FB8",
    color = "white",
    alpha = 0.8
  ) 
# ggplot(train_df_long) +
#   geom_histogram(
#     aes(intensity),
#     bins = 30,
#     fill = "#2C7FB8",
#     color = "white",
#     alpha = 0.8
#   ) +
#   scale_y_log10() +
#   labs(y = "Count (log scale)")

```

:::{.notes}
- Most pixel values are zero, but out of the 784 pixels, 703 have an intensity > 0 at least once
:::

## Training Data tSNE Visualization

```{r tnse, fig.dim=c(10,6)}
#------------------------------------------------------------
#----------------distribution of mnist training--------------
#------------------------------------------------------------
#normalized pixel values and apply t-SNE
# tsne_results <- Rtsne(
#   train_images,
#   dim = 2,
#   perplexity = 30,
#   max_iter = 1000
# )
# df_tsne <- tibble(Dim1 = tsne_results$Y[, 1],
#                   Dim2 = tsne_results$Y[, 2],
#                   digit = train_df$labels)
#
# saveRDS(df_tsne, "df_tsne")
df_tsne <- readRDS("df_tsne")
ggplot(df_tsne, aes(Dim1, Dim2, color = digit)) +
  geom_point(alpha = .5)
```

:::{.notes}
- We adopt t-distributed Stochastic Neighbor Embedding(t-SNE) to represent the data in 2D.See [T-SNE Exploration by TusVasMit](https://rpubs.com/TusVasMit/T-SNEExploration) for more details.
- tSNE is a technique for dimension reduction that is well-suited for visualizing high-dimensional data in a lower-dimensional space
:::

# Traditional ML

## Neural networks

:::{.fragment}
Feedforward neural network with structure:
:::

:::{.incremental}
- Input layer: Consists of neurons that receives the input data. Each neuron in the input layer represents a feature of the input data
- Hidden layer: One(or more) hidden layers placed between the input and output layers, responsible for capturing complex patterns
- Output layer: Final output of the network; Number of neurons represents the number of digits
:::

:::{.notes}
- We only use one hidden layer in our case
:::


## NN with regularization

:::{.incremental}
- Depending on model, number of network weights > size of training data
    - This leads to overfitting
- We considered two approaches to overfitting:
  - Dropout learning: Like RF, randomly removes fraction of units in layer during model fitting
  - Regularization: Impose penalties on parameters like lasso, ridge, etc.
:::


## Specific NN models considered

:::{.fragment}
- NN with dropout regularization
:::
:::{.fragment}
- NN with ridge regularization
:::
:::{.fragment}
- NN with lasso regularization 
:::

:::{.notes}
- NN w/ dropout is a standard approach to prevent overfitting by randomly dropping out a fraction of neurons during training
- NN w/ ridge regularization adds a penalty term to the loss function that is proportional to the square of the weights, helping to prevent overfitting by discouraging large weights
- NN w/ lasso regularization adds a penalty term to the loss function that is proportional to the absolute value of the weights, promoting sparsity in the model by driving some weights to zero
:::

## Multinomial logistic regression

:::{.incremental}
- Multinomial logistic regression equivalently represented by NN with no hidden layers
- Output layer with softmax
    - $f_m(X) = Pr(Y = m | X) = \frac{e^{Z_m}}{\sum_\limits{l \in K}e^{Z_l}}$
:::

:::{.notes}
- $m$ is the class label
- $Z_m$ is the output of the model for class $m$
- $X$ is the input data
- $K$ is the set of all class labels
- $l$ is a class label in $K$
:::

## NN Fitting

:::{.incremental}
- Train the network for 30 epochs with a batch size of 128
    - SGD updates weights for each batch
- Images are presented in batches of 128 and SGD updates weights after each batch
- Each epoch processes all 60,000 training images
<!-- - After 30 epochs, the network’s error (loss) is visualized on the left graph -->
<!-- - Accuracy shown on the right graph -->
- Classification correct if largest output value matches target label
:::

:::{.notes}
- Epoch: One complete pass through the entire training dataset 
- SGD is an optimization algorithm used to minimize the loss function by updating the model's weights based on the gradients of the loss with respect to the weights
- This process iteratively adjusts the weights to improve the model's performance on the training data
:::

# Proposed Methodology {#sec-PropMethod}

## TDA Workflow 

```{tikz}
%| echo: false

\usetikzlibrary{
    positioning, 
    arrows.meta, 
    shapes.geometric, 
    fit, 
    calc
}

\begin{tikzpicture}[
    % Adjusted node distances for better spacing
    node distance = 1.2cm and 2cm,
    every node/.style={
        draw, 
        thick, 
        rounded corners, 
        align=center, 
        minimum height=1.3cm,
        font=\sffamily
    },
    data/.style={fill=green!20, text width=3cm},
    prior/.style={fill=yellow!30, text width=4cm},
    posterior/.style={fill=blue!20, text width=4cm},
    result/.style={fill=red!20, text width=3.5cm},
    process/.style={text width=4cm},
    arrow/.style={->, >=Stealth, thick},
    connector/.style={draw=none, font=\sffamily\Huge},
    % A dedicated style for labels on arrows (edges)
    edge_label/.style={draw=none, midway, fill=none, font=\sffamily}
]

% == Column 1 & 2: Data and PD Calculation ==
% Position nodes in the first two columns
\node[data] (train) {Train Data \\ (60,000 images)};
\node[process, right=of train] (calc_pd_train) {Calculate Train PDs \\ (for dim0 \& dim1)};

% Increased vertical distance for a clearer separation of train/test paths
\node[data, below=3.75cm of train] (test) {Test Data \\ (10,000 images)};
\node[process, right=of test] (calc_pd_test) {Calculate Test PDs \\ (for dim0 \& dim1)};

% == Column 3: Bayesian Model Training ==
% Position this block relative to the training data processing nodes
\node[process, right=of calc_pd_train] (likelihoods) {Likelihood Surfaces from Train PDs \\ (for digits 0-9)};
%\node[connector, right=of likelihoods] (update_op) {$\otimes$};
\node[connector, right=of likelihoods] (update_op) {$\odot$};
\node[prior, right=of update_op] (priors) {Uninformative Priors \\ (for digits 0-9)};
\node[posterior, below=of update_op] (posteriors) {Posterior Surfaces \\ (for digits 0-9)};

% Bounding box for the Bayesian update process
\node[draw, dashed, inner sep=0.4cm, fit=(priors) (likelihoods) (update_op) (posteriors), label={[font=\sffamily\bfseries]above:Bayesian Update}] (model_box) {};

% == Column 4: Classification ==
% Position the classification node vertically centered between its inputs for a balanced look
\node[process, below=of posteriors] (calc_dist) {Calculate Distances to all Posteriors \\ Distance = $d_{0} + d_{1}$};
\node[result, below=of calc_dist] (classify) {Classify as \\ argmin(Distance)};

% == Arrows ==
% Connect nodes with clearer, non-overlapping paths
\draw[arrow] (train) -- (calc_pd_train);
\draw[arrow] (test) -- (calc_pd_test);

% Bayesian model flow
\draw[arrow] (calc_pd_train) -- (likelihoods);
\draw[arrow] (priors) |- (posteriors);
\draw[arrow] (likelihoods) |- (posteriors);

% Classification flow
% Use |- routing to different anchors (north west and south west) to keep lines clean
\draw[arrow] (posteriors) -- (calc_dist);
\draw[arrow] (calc_pd_test) -- (calc_dist);

% Arrow with a nicely placed label for the distance formula
\draw[arrow] (calc_dist) -- (classify);
    %node[edge_label, right=0.2cm] {Distance = \\ $(1-\lambda)d_{0} + \lambda d_{1}$};
    %node[edge_label, right=0.2cm] {Distance = $d_{0} + d_{1}$};
    
\end{tikzpicture}
```

Based on framework by [@Maroulas2020-sp]

:::{.notes}
- images 28x28
- 60,000 train PDs made
- one uniform prior for all 10 digits which is a uniform distribution over the persistence diagram space $\pi(\theta) = \sum_\limits{j = 1}^J \omega_j \mathcal{N}(\theta | \mu_j, \Sigma_j)$
    - J = 1, $\omega = 1$
    - mean vector $[.5, .5]^\top$
    - covariance is $\Sigma\_1 = 1 \cdot I\_2 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$
- 20 posteriors, 10 for dim0, 10 for dim1
- 10k test PDs
- 100k distances for dim0, 100k for dim1
:::

# Analysis {#sec-Analysis}

## ML Model Performance

```{r mlmodelperf, fig.dim=c(10,6)}
# --- Load Data ---
nn_dropout_conf_mat <- readRDS("nn_dropout_conf_mat")
nn_ridge_conf_mat <- readRDS("nn_ridge_conf_mat")
nn_lasso_conf_mat <- readRDS("nn_lasso_conf_mat")
mlogit_conf_mat <- readRDS("mlogit_conf_mat")

# --- Prepare data frames for all plots ---
df_dropout <- as.data.frame.matrix(nn_dropout_conf_mat$table) %>%
  mutate(true_label = fct_rev(rownames(.))) |>
  pivot_longer(cols = -true_label, names_to = "predicted_label", values_to = "count")

df_ridge <- as.data.frame.matrix(nn_ridge_conf_mat$table) %>%
  mutate(true_label = fct_rev(rownames(.))) |>
  pivot_longer(cols = -true_label, names_to = "predicted_label", values_to = "count")

df_lasso <- as.data.frame.matrix(nn_lasso_conf_mat$table) %>%
  mutate(true_label = fct_rev(rownames(.))) |>
  pivot_longer(cols = -true_label, names_to = "predicted_label", values_to = "count")

df_mlogit <- as.data.frame.matrix(mlogit_conf_mat$table) %>%
  mutate(true_label = fct_rev(rownames(.))) |>
  pivot_longer(cols = -true_label, names_to = "predicted_label", values_to = "count")

# --- Determine the global scale limits for the fill color ---
global_limits <- range(
  bind_rows(df_dropout, df_ridge, df_lasso, df_mlogit) |> pull(count)
)

# --- Define a plotting function to reduce code repetition ---
create_heatmap <- function(data, title, limits) {
  data |>
    mutate(text_color = if_else(count < (0.6 * max(limits)), "black", "white")) |>
    ggplot(aes(x = predicted_label, y = true_label, fill = count)) +
    geom_tile(color = "gray50", linewidth = 0.5) +
    geom_text(aes(label = count, color = text_color), size = 2) +
    scale_fill_viridis_c(
      option = "magma",
      direction = -1,
      limits = limits # Apply global limits here
    ) +
    scale_color_manual(values = c("black" = "black", "white" = "white"), guide = "none") +
    coord_fixed() +
    theme_minimal(base_size = 12) +
    labs(
      x = "Predicted Label",
      y = "True Label",
      fill = "Count",
      title = title
    ) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
      axis.text.y = element_text(face = "bold"),
      panel.grid = element_blank()
    )
}

# --- Create each plot using the function ---
p_dropout <- create_heatmap(df_dropout, "Neural Net (Dropout)", global_limits)
p_ridge <- create_heatmap(df_ridge, "Neural Net (Ridge)", global_limits)
p_lasso <- create_heatmap(df_lasso, "Neural Net (Lasso)", global_limits)
p_mlogit <- create_heatmap(df_mlogit, "Multinomial Logistic", global_limits)


# --- Combine Plots with Patchwork, collecting guides ---
(p_dropout + p_ridge) / (p_lasso + p_mlogit) +
  plot_annotation(
    theme = theme(plot.title = element_text(hjust = 0.5, size = 18, face = "bold"))
  ) &
  theme(legend.position = "none")

```

## Proposed Method Analysis

```{r post_digit_plot0, fig.dim=c(10,6)}
posterior_list_dim0 <- readRDS("btda/posterior_list_dim0.rds")
posterior_dim0_df <- posterior_list_dim0 |>
  set_names(0:9) |>
  list_rbind(names_to = "digit")

ggplot(posterior_dim0_df, aes(x = birth, y = persistence, fill = intensity)) +
  geom_raster() +
  facet_wrap(~digit, ncol = 5) +
  scale_fill_viridis_c(option = "magma") +
  labs(
    title = "Posterior Densities for Digit Components (Dimension 0)",
    x = "Birth",
    y = "Persistence",
    fill = "Intensity"
  ) +
  theme_minimal() +
  theme(strip.text = element_text(size = 12, face = "bold"))
```

## Proposed Method Analysis

```{r post_digit_plot1, fig.dim=c(10,6)}
posterior_list_dim1 <- readRDS("btda/posterior_list_dim1.rds")
posterior_dim1_df <- posterior_list_dim1 |>
  set_names(0:9) |>
  list_rbind(names_to = "digit")

ggplot(posterior_dim1_df, aes(x = birth, y = persistence, fill = intensity)) +
  geom_raster() +
  facet_wrap(~digit, ncol = 5) +
  scale_fill_viridis_c(option = "magma") +
  labs(
    title = "Posterior Densities for Digit Loops (Dimension 1)",
    x = "Birth",
    y = "Persistence",
    fill = "Intensity"
  ) +
  theme_minimal() +
  theme(strip.text = element_text(size = 12, face = "bold"))
```

## Proposed Method Analysis

```{r tda01confmatrix, fig.dim=c(10,6)}
full_tda_01_results <- readRDS("btda/full_tda_01_results.rds")

conf_matrix <- table(
  true_label = full_tda_01_results$true_label,
  predicted_label = full_tda_01_results$predicted_label
)

as.data.frame.matrix(conf_matrix) %>%
  mutate(true_label = fct_rev(rownames(.))) |>
  pivot_longer(
    cols = -true_label,
    names_to = "predicted_label",
    values_to = "count"
  ) |>
  
  mutate(text_color = if_else(count < 450, "black", "white")) |>
  ggplot(aes(x = predicted_label, y = true_label, fill = count)) +
  geom_tile(color = "gray50", linewidth = 0.5) +
  
  geom_text(aes(label = count, color = text_color), size = 3.5) +
  scale_fill_viridis_c(option = "magma", direction = -1) +
  
  scale_color_manual(values = c("black" = "black", "white" = "white"), guide = "none") +
  coord_fixed() +
  theme_minimal(base_size = 12) +
  labs(
    x = "Predicted Label",
    y = "True Label",
    fill = "Count"
  ) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
    axis.text.y = element_text(face = "bold"),
    panel.grid = element_blank()
  )

# accuracy <- mean(full_results$true_label == full_results$predicted_label, na.rm = TRUE)

```

## Proposed Method Analysis

```{r tda01metrics}
multi_metrics <- metric_set(
  yardstick::accuracy, 
  yardstick::precision, 
  yardstick::recall, 
  yardstick::f_meas
)

results_for_metrics <- full_tda_01_results |>
  select(true_label, predicted_label) |>
  mutate(
    true_label = factor(true_label, levels = 0:9),
    predicted_label = factor(predicted_label, levels = 0:9)
  ) 

proposed_accuracy <- multi_metrics(results_for_metrics, truth = true_label, estimate = predicted_label) |> 
  janitor::clean_names() |>
  filter(metric == "accuracy") |> pull(estimate)

#run the Mnist ML method code first, and take results from there
# ml_results <- tibble(
#   method = c("multinomial", "dropout nn", "ridge nn", "lasso no"),
#   accuracy = c(mlogit_acc, nn_dropout_accu, nn_ridge_accu, nn_lasso_accu)
# )
# saveRDS(ml_results, "ml_results")

ml_results <- readRDS("ml_results") #still a tibble so easy to edit
ml_results <- ml_results %>%
  mutate(method = str_replace(method, "lasso no", "lasso nn"))
ml_results |>
  add_row(method = "proposed", accuracy = proposed_accuracy) |>
  kable(digits = 4)
```

# TDA + ML {#sec-TDAML}

Extension of [@DBLP:journals/corr/abs-1910-08345]

## Filtering {.smaller}

:::{.incremental}
- Grayscale: Uses the image's original pixel intensities directly.  Pixels are included in the complex as their intensity value passes a growing threshold
- Height: Assigns a value to each pixel based on its projection onto a chosen direction vector, essentially measuring its "height" from a specific angle
- Radial: Assigns a value to each pixel based on its distance from a chosen center point
- Dilation: Assigns each pixel a value corresponding to its shortest distance to a foreground (value=1) pixel
    - This has the effect of "growing" or "dilating" the digit
- Erosion: Inverse of dilation
    - Dilation of the inverted image 
    - This "shrinks" or "erodes" the digit
- Density: Assigns each pixel a value based on the number of foreground neighbors within a given radius
    - Combination of dilation & erosion
:::

:::{.notes}
- height: 8 heights, each corner and edge
- radial: 9 chosen
- dilation, erosion, density: multiple radii chosen, 2,4,6
:::

## Filtering 

```{r filtrations, fig.dim=c(10,6)}
image_index <- which(mnist$train$labels == 8)[100]
image_matrix <- matrix(mnist$train$images[image_index, ], nrow = 28, byrow = TRUE)

binary_image <- ifelse(image_matrix > 0.4, 1, 0)
img_c <- as.cimg(binary_image)
img_c_inv <- as.cimg(1 - binary_image)

coords <- as.matrix(expand.grid(y = 1:28, x = 1:28))

height_filt_td <- matrix(coords %*% c(0, 1), nrow = 28, byrow = TRUE) * binary_image
height_filt_lr <- matrix(coords %*% c(1, 0), nrow = 28, byrow = TRUE) * binary_image

radial_filt_c <- matrix(sqrt(rowSums(sweep(coords, 2, c(14, 14), "-")^2)), nrow = 28, byrow = TRUE) * binary_image
radial_filt_b <- matrix(sqrt(rowSums(sweep(coords, 2, c(14, 28), "-")^2)), nrow = 28, byrow = TRUE) * binary_image

dilation_filt <- as.matrix(distance_transform(img_c, 1))
erosion_filt <- as.matrix(distance_transform(img_c_inv, 1))
density_filt <- as.matrix(boxblur(img_c, 3))

plot_filtration <- function(data_matrix, title) {
  as.data.frame(data_matrix) |>
    mutate(y = row_number()) |>
    pivot_longer(-y, names_to = "x", values_to = "value") |>
    mutate(x = as.integer(gsub("V", "", x))) |>
    ggplot(aes(x = x, y = y, fill = value)) +
    geom_raster() +
    scale_y_reverse() +
    coord_equal() +
    scale_fill_viridis_c(option = "magma") +
    theme_void() +
    labs(title = title) +
    theme(plot.title = element_text(hjust = 0.5, face = "bold"), legend.position = "none")
}

p_orig <- plot_filtration(image_matrix, "Original Grayscale")
p_bin <- plot_filtration(binary_image, "Binarized")
p_height_td <- plot_filtration(height_filt_td, "Height (Top-Down)")
p_height_lr <- plot_filtration(height_filt_lr, "Height (Left-Right)")
p_radial_c <- plot_filtration(radial_filt_c, "Radial (Center)")
p_radial_b <- plot_filtration(radial_filt_b, "Radial (Bottom)")
p_dilation <- plot_filtration(dilation_filt, "Dilation")
p_erosion <- plot_filtration(erosion_filt, "Erosion")
p_density <- plot_filtration(density_filt, "Density (r=3)")

(p_orig + p_bin + p_dilation) /
(p_height_td + p_height_lr + p_erosion) /
(p_radial_c + p_radial_b + p_density)
```

## Analysis: TDA+ML Performance

```{r loadallfeatures}
all_features_train_data <- readRDS("rdss/all_features_train_data.rds")
all_features_test_data <- readRDS("rdss/all_features_test_data.rds")
all_features_rf_model <- readRDS("rdss/all_features_rf_model.rds")
all_features_svm_model<- readRDS("rdss/all_features_svm_model.rds")
```

```{r modelcomparison, fig.dim=c(10,6)}
model_comparison <- resamples(list(RandomForest = all_features_rf_model, SVM = all_features_svm_model))
# dotplot(model_comparison,
#   scales = list(x = list(rot = 45)))
tidy_comparison <- model_comparison$values |>
  pivot_longer(
    cols = -Resample,
    names_to = "key",
    values_to = "value"
  ) |>
  separate(key, into = c("Model", "Metric"), sep = "~")

ggplot(tidy_comparison, aes(x = value, y = Model, color = Model)) +
  geom_point(size = 4, alpha = 0.2) +
  scale_color_brewer(palette = "Set2") +
  scale_x_continuous(limits = c(NA, 1.0)) +
  facet_wrap(~ Metric, ncol = 4) +
  labs(
    x = "Metric Value",
    y = NULL,
    color = "Model Type"
  ) +
  theme_light(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 18, margin = margin(b = 5)),
    plot.subtitle = element_text(size = 13, margin = margin(b = 15)),
    strip.text = element_text(face = "bold", size = 10, color = "white"),
    strip.background = element_rect(fill = "#525252", color = "white"),
    legend.position = "bottom",
    panel.grid.major.y = element_blank(),
    panel.grid.minor.x = element_blank()
  ) +
  guides(color = guide_legend(override.aes = list(alpha = 1)))
```

:::{.notes}
* **Accuracy**: The proportion of all predictions that were correct. It's a great starting point but can be misleading for imbalanced datasets.
    * *Formula*: $(TP + TN) / (TP + TN + FP + FN)$

* **Balanced Accuracy**: The average of Sensitivity and Specificity. It provides a more robust measure when classes are imbalanced by giving equal importance to both positive and negative classes.
    * *Formula*: $(Sensitivity + Specificity) / 2$

* **Sensitivity (Recall or True Positive Rate)**: Answers: "Of all the actual positive cases, how many did we correctly identify?" High sensitivity is crucial when you cannot afford to miss a positive case (e.g., disease screening).
    * *Formula*: $TP / (TP + FN)$

* **Specificity (True Negative Rate)**: Answers: "Of all the actual negative cases, how many did we correctly identify?" High specificity is important when a false positive is very costly.
    * *Formula*: $TN / (TN + FP)$

* **Precision (Positive Predictive Value)**: Answers: "Of all the cases we predicted as positive, how many were actually positive?" High precision matters when the cost of acting on a false positive is high.
    * *Formula*: $TP / (TP + FP)$

* **Negative Predictive Value (NPV)**: Answers: "Of all the cases we predicted as negative, how many were actually negative?" It's the counterpart to Precision.
    * *Formula*: $TN / (TN + FN)$

* **F1 Score**: The harmonic mean of Precision and Sensitivity. It provides a single score that balances the concerns of both metrics, making it useful when both false positives and false negatives are important to minimize.
    * *Formula*: $2 \times (Precision \times Sensitivity) / (Precision + Sensitivity)$

* **Kappa (Cohen's Kappa)**: Measures the agreement between the model's predictions and the actual labels, corrected for the probability of agreement occurring by chance. A score of 1 indicates perfect agreement, 0 indicates agreement equivalent to random chance, and negative values indicate agreement worse than random.

* **Detection Rate**: The proportion of the entire dataset that consists of true positives. It's less commonly used than accuracy but gives a sense of how often the positive case is found in the population.
    * *Formula*: $TP / (TP + TN + FP + FN)$
:::

## Analysis: TDA+ML Confusion Matrix

```{r allfeaturesheatmap, fig.dim=c(10,6)}
test_indices <- 1:10000
predictions <- predict(all_features_svm_model, newdata = all_features_test_data)
results <- tibble(
  true_label = factor(mnist$test$labels[test_indices]),
  predicted_label = predictions
)

conf_matrix_tda_full <- table(
  true_label = results$true_label,
  predicted_label = results$predicted_label
)

as.data.frame.matrix(conf_matrix_tda_full) %>%
  mutate(true_label = fct_rev(rownames(.))) |>
  pivot_longer(
    cols = -true_label,
    names_to = "predicted_label",
    values_to = "count"
  ) |>
  
  mutate(text_color = if_else(count < 450, "black", "white")) |>
  ggplot(aes(x = predicted_label, y = true_label, fill = count)) +
  geom_tile(color = "gray50", linewidth = 0.5) +
  
  geom_text(aes(label = count, color = text_color), size = 3.5, fontface = "bold") +
  scale_fill_viridis_c(option = "magma", direction = -1) +
  
  scale_color_manual(values = c("black" = "black", "white" = "white"), guide = "none") +
  coord_fixed() +
  theme_minimal(base_size = 12) +
  labs(
    x = "Predicted Label",
    y = "True Label",
    fill = "Count"
  ) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
    axis.text.y = element_text(face = "bold"),
    panel.grid = element_blank()
  )

# final_accuracy <- mean(results$true_label == results$predicted_label, na.rm = TRUE)
# print(paste("Final Test Set Accuracy (Random Forest):", scales::percent(final_accuracy, accuracy = 0.1)))
```

:::{.notes}
- accuracy on 22 features rbf-svm: 93.1%
:::

## Analysis: TDA+ML Variable Importance

```{r varimportance, fig.dim=c(10,6)}
rf_model_importance <- varImp(all_features_rf_model, scale = TRUE)
svm_model_importance <- filterVarImp(
  x = all_features_train_data[, -which(names(all_features_train_data) == "label")],
  y = all_features_train_data$label
)

rf_plot <- rf_model_importance$importance |>
  as.data.frame() %>%
  mutate(variable = rownames(.)) |>
  pivot_longer(
    cols = -variable,
    names_to = "class",
    values_to = "importance"
  ) |>
  mutate(
    importance = (importance - min(importance)) / (max(importance) - min(importance))
  ) |>
  mutate(variable = fct_reorder(variable, importance, .fun = max, .desc = TRUE)) |>
  ggplot(aes(x = class, y = variable, fill = importance)) +
  geom_tile(color = "white", linewidth = 0.4) +
  scale_fill_viridis_c(
    option = "magma", 
    direction = -1,
    limits = c(0, 1)
  ) +
  labs(
    x = "Class",
    y = "Variable",
    fill = "Importance",
    title = "Random Forest Variable Importance"
  ) +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

svm_plot <- as.data.frame(svm_model_importance) %>%
  mutate(variable = rownames(.)) |>
  pivot_longer(
    cols = -variable,
    names_to = "class",
    values_to = "importance"
  ) |>
  mutate(
    importance = (importance - min(importance)) / (max(importance) - min(importance))
  ) |>
  mutate(
    variable = fct_reorder(variable, importance, .fun = max, .desc = TRUE),
    class = gsub("X", "", class)
    ) |>
  ggplot(aes(x = class, y = variable, fill = importance)) +
  geom_tile(color = "white", linewidth = 0.4) +
  scale_fill_viridis_c(option = "magma", direction = -1) +
  labs(
    x = "Class",
    y = NULL, 
    fill = "Importance",
    title = "SVM Variable Importance"
  ) +
  theme_minimal(base_size = 12) +
  theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      axis.text.y = element_blank() 
      )

rf_plot + svm_plot + 
  plot_layout(guides = 'collect')

```

:::{.notes}
- Persistent Entropy summarizes complexity of diagram
    - $PE(D) = - \sum_\limits{i=1}^n \frac{l_i}{L(D)}\log\left(\frac{l_i}{L(D)}\right)$
    - $D = \{(b_i, d_i)\}^n_{i=1}$ persistence diagram(barcode)
    - $l_i = d_i - b_i$
    - $L(D) = \sum_{i=1}^n l_i$
- Amplitude is distance to distance to empty diagram
    - Waserstein $L_p$ nrom of the vector of feature lifetimes
    - Bottleneck: max liftime among all features in diagram
    - Kernel based - transform to betti curve, persistence landscape, heat kernel, then take $L_1$ or $L_2$ norm
:::

# Future Work {#sec-Results}

## Future Work

:::{.incremental}
- Proposed Method
    - Use various distance metrics: entropy & amplitude
    - $\lambda$ to weight dimensional distance
    - Model different complex/filter as bayes PPP
- TDA + ML
    - Use other ML methods
    - Use finer tuning on paramters
    - Expand TDA features
    - Feature selection
:::

# References {#sec-References}

## References

::: {#refs .smaller}
:::