---
title: "Code Associated with Explainable Machine Learning for Functional Data"
subtitle: "Code in Supplement"
output: 
  pdf_document:
    toc: true
header-includes: 
  - \renewcommand{\thefigure}{S\arabic{figure}}
---

```{r setup, message = FALSE}
# R markdown options
knitr::opts_chunk$set(
   message = FALSE,
   fig.align = "center",
   dpi = 600,
   warning = FALSE,
   error = FALSE
)
knitr::opts_knit$set(eval.after = "fig.cap")

# Using v0.1.7 of veesa R package. Uncomment line below to install from CRAN:
#install_version("veesa", version = "0.1.7")

# Load packages
library(cowplot)
library(dplyr)
library(forcats)
library(ggplot2)
library(latex2exp)
library(purrr)
library(randomForest)
library(reticulate)
library(stringr)
library(tidyr)
library(veesa)
library(wesanderson)

# Specify the conda environment to use
use_condaenv(condaenv = "veesa", required = T)

# Specify colors for groups
col_2groups = wes_palettes$Royal1[2:1]

# Specify colors for PC direction plots
col_pl1 = "#784D8C"
col_pl2 = "#A289AE"
col_pl3 = "#CAC0D2"
col_mi1 = "#EA9B44"
col_mi2 = "#EBBC88"
col_mi3 = "#F3DABC"
col_pcdir_2sd = c(col_pl2, col_pl1, "black", col_mi1, col_mi2)
col_pcdir_3sd = c(col_pl3, col_pl2, col_pl1, "black", col_mi1, col_mi2, col_mi3)
col_5groups <-
  wes_palette(name = "Royal1", n = 5, type = "continuous")[c(1,2,3,5,4)]

# Specify inkjet colors
col_inkjet = c("#0EB4AD", "#8f1883", "#d9c004")

# Specify some high/low colors for contrast throughout
col_low = wes_palette(name = "Royal1", n = 30, type = "continuous")[4]
col_high = wes_palette(name = "Royal1", n = 30, type = "continuous")[29]

# Specify the font for all plots
ff = "Helvetica"

# Specify a file path
# Assumes the working directory is the location of this script (code/)
fp = "../"

# Specify results folders
fp_res = paste0(fp, "results/")
fp_res_sp = paste0(fp_res, "shifted-peaks/")
fp_res_hct = paste0(fp_res, "hct/")
fp_res_inkjet = paste0(fp_res, "inkjet/")
```

# Comparing VEESA Pipeline to Cross-Sectional Approach with the Shifted Peaks Data {#shifted-peaks}

```{r shifted-peaks-data}
# Get shifted peaks data from veesa R package
sp = veesa::shifted_peaks

# Load the sim data
sp_data = sp$data |> mutate(group = factor(group))

# Create long versions of training/testing data sets
sp_train_long = sp_data |> filter(data == "Training")
sp_test_long = sp_data |> filter(data == "Testing")

# Create wide versions of training and testing data
sp_train_wide <- 
  sp_train_long |>
  select(-t) |>
  mutate(index = paste0("t", index)) |>
  pivot_wider(names_from = index, values_from = y)
sp_test_wide <- 
  sp_test_long |>
  select(-t) |>
  mutate(index = paste0("t", index)) |>
  pivot_wider(names_from = index, values_from = y)

# Determine number of observations in training/testing data
sp_train_nobs = dim(sp_train_wide)[1]
sp_test_nobs = dim(sp_test_wide)[1]

# Determine the time points where values are observed
sp_times = sort(unique(sp_train_long$t))

# Load ESA version of the train/test data
sp_train_esa_fp = paste0(fp_res_sp, "shifted-peaks-train-esa.rds")
sp_train_esa = readRDS(sp_train_esa_fp)
sp_test_esa_fp = paste0(fp_res_sp, "shifted-peaks-test-esa.rds")
sp_test_esa = readRDS(sp_test_esa_fp)
  
# Prepare data for model
sp_train_model_data <- 
  sp_train_esa$fpca_res$coef |> 
  data.frame() |> 
  mutate(group = factor(sp_train_wide$group))
sp_test_model_data <- 
  t(sp_test_esa$coef) |> 
  data.frame() |> 
  mutate(group = factor(sp_test_wide$group))

# Specify the number of replications used to compute PFI with the sim data
sp_K = 10
```

```{r shifted-peaks-true-means}
# True functional means
sp_true_means <- sp$true_means
```

```{r shifted-peaks-cross-sectional}
# Prepare cross-sectional data
sp_data_cs_train <-
  sp_train_wide |>
  select(-id, -data) |> 
  mutate(group = factor(group))
sp_data_cs_test <-
  sp_test_wide |> 
  select(-id, -data) |> 
  mutate(group = factor(group))

# Cross-sectional approach
sp_rf_cs_fp = paste0(fp_res_sp, "shifted-peaks-rf-cross-sectional.rds")
if (file.exists(sp_rf_cs_fp)) {
  sp_rf_cs = readRDS(sp_rf_cs_fp)
} else {
  set.seed(20210301)
  sp_rf_cs <-
    randomForest::randomForest(
      group ~ .,
      sp_data_cs_train,
      importance = TRUE
    )
  saveRDS(sp_rf_cs, sp_rf_cs_fp)
}

# Compute the model predictions on train/test data
sp_pred_cs_train = predict(sp_rf_cs, sp_data_cs_train)
sp_pred_cs_test = predict(sp_rf_cs, sp_data_cs_test)

# Compute the accuracy on train/test data
sp_acc_cs_train = sum(sp_pred_cs_train == sp_train_wide$group) / sp_train_nobs
sp_acc_cs_test = sum(sp_pred_cs_test == sp_test_wide$group) / sp_test_nobs

# Compute PFI using accuracy for the cross-sectional approach
sp_pfi_acc_cs_fp <-
  paste0(fp_res_sp, "shifted-peaks-pfi-acc-cross-sectional.rds")
if (file.exists(sp_pfi_acc_cs_fp)) {
  sp_pfi_acc_cs = readRDS(sp_pfi_acc_cs_fp)
} else {
  set.seed(20210921)
  sp_pfi_acc_cs <- 
    compute_pfi(
      x = sp_data_cs_test[, -1],
      y = sp_data_cs_test$group,
      f = sp_rf_cs,
      K = sp_K, 
      metric = "accuracy"
    )
  saveRDS(sp_pfi_acc_cs, sp_pfi_acc_cs_fp)
}

# Compute PFI using log-loss for the cross-sectional approach
sp_pfi_ll_cs_fp = paste0(fp_res_sp, "shifted-peaks-pfi-ll-cross-sectional.rds")
if (file.exists(sp_pfi_ll_cs_fp)) {
  sp_pfi_ll_cs = readRDS(sp_pfi_ll_cs_fp)
} else {
  set.seed(20210921)
  sp_pfi_ll_cs <- 
    compute_pfi(
      x = sp_data_cs_test[, -1],
      y = sp_data_cs_test$group,
      f = sp_rf_cs,
      K = sp_K, 
      metric = "logloss"
    )
  saveRDS(sp_pfi_ll_cs, sp_pfi_ll_cs_fp)
}
```

```{r shifted-peaks-jfpca}
# Random forest on jfPCA with aligned data
sp_rf_jfpca_fp = paste0(fp_res_sp, "shifted-peaks-rf-jfpca.rds")
sp_rf_jfpca = readRDS(sp_rf_jfpca_fp)

# Compute the model predictions on the test data
sp_train_pred_jfpca = predict(sp_rf_jfpca, sp_train_model_data |> select(-group))
sp_test_pred_jfpca = predict(sp_rf_jfpca, sp_test_model_data |> select(-group))

# Compute the accuracy on the training data
sp_train_acc_jfpca <-
  sum(sp_train_pred_jfpca == sp_train_model_data$group) / sp_train_nobs
sp_test_acc_jfpca <-
  sum(sp_test_pred_jfpca == sp_test_model_data$group) / sp_test_nobs

# Compute PFI using log-loss with the VEESA pipeline
sp_pfi_ll_jfpca_fp = paste0(fp_res_sp, "shifted-peaks-pfi-ll-jfpca.rds")
if (file.exists(sp_pfi_ll_jfpca_fp)) {
  sp_pfi_ll_jfpca = readRDS(sp_pfi_ll_jfpca_fp)
} else {
  set.seed(20210921)
  sp_pfi_ll_jfpca <-
    compute_pfi(
      x = sp_test_model_data |> select(-group),
      y = sp_test_model_data$group,
      f = sp_rf_jfpca,
      K = sp_K,
      metric = "logloss"
    )
  saveRDS(sp_pfi_ll_jfpca, sp_pfi_ll_jfpca_fp)
}
```

```{r shifted-peaks-check}
sp_pfi_acc_cs_sd <-
  data.frame(sp_pfi_acc_cs$pfi_single_reps) |> 
  summarise_all(.funs = sd) |> 
  as.numeric() |> 
  unique()

if (sp_pfi_acc_cs_sd > 0)
  stop("Some times have PFI standard deviation > 0.")
```

```{r shifted-peaks-pfi-cutoff}
sp_veesa_pfi_cutoff = 0.005
```

```{r figS1, fig.width = 6, out.width = '4.25in', fig.cap = "Spearman correlations between all pairs of cross-sectional variables from the shifted peaks data."}

# Specify the figure size (for determining font size)
fS1_fs = 14

# Prepare the correlations for the plot
corr_data <-
  sp_train_wide |>
  select(-id, -group, -data) |>
  cor(method = "spearman") |>
  as.data.frame() |>
  mutate(index_row = 1:n()) |>
  reshape2::melt(id = "index_row") |>
  rename(index_col = variable) |>
  mutate(index_row = as.numeric(str_remove(index_row, "t")),
         index_col = as.numeric(str_remove(index_col, "t"))) |>
  left_join(sp_data |>
              select(index, t) |>
              distinct(),
            by = c("index_row" = "index")) |>
  rename(t_row = t) |>
  left_join(sp_data |>
              select(index, t) |>
              distinct(),
            by = c("index_col" = "index")) |>
  rename(t_col = t) |>
  select(index_col, t_col, index_row, t_row, value)

# Plot the correlation matrix
corr_data |>
  ggplot(aes(x = t_col, y = t_row)) +
  geom_tile(aes(fill = value, color = value)) +
  coord_equal() +
  guides(color = guide_colourbar(barheight = 4, barwidth = 0.5)) +
  scale_color_gradient2(
    low = col_low,
    high = col_high,
    mid = "white",
    limits = c(-1, 1)
  ) +
  scale_fill_gradient2(
    low = col_low,
    high = col_high,
    mid = "white",
    limits = c(-1, 1)
  ) +
  guides(color = guide_colourbar(barheight = 6, barwidth = 1)) +
  labs(
    x = TeX('$t$'),
    y = TeX('$t$'),
    fill = "Correlation",
    color = "Correlation"
  ) +
  theme_bw(base_family = ff, base_size = fS1_fs) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.ticks = element_blank(),
    aspect.ratio = 1
  )
```

```{r figS2, fig.width = 15, fig.height = 10, out.width = '6.5in', fig.cap = "(Top Left) The solid lines represent the true group means. The dots and error bars represent the cross-sectional group means plus and minus one cross-sectional standard deviation. (Top Right) PFI values computed on the test data using the metric of accuracy for a random forest trained using the cross-sectional approach. (Bottom Left) PFI values computed on the test data using the metric of log-loss for a random forest trained using the cross-sectional approach. (Bottom Right) PFI values computed on the test data using the metric of log-loss for a random forest trained using the VEESA pipeline."}

# Specify the figure size (for determining font size)
fS2_fs = 16

# Specify limits associated with the region with high 
# importance from the cross-sectional approach
f1_ll = -0.7
f1_ul = 3.5

# Specify an aspect ratio for all plots relating to the sim data
sp_ar = 0.5

# Plot of cross-sectional means and standard deviations by group
plot_sp_cs_means <-
  sp_data |>
  group_by(t, group) |>
  summarise(mean_y = mean(y),
            sd_y = sd(y),
            .groups = "drop") |>
  ggplot(aes(x = t, y = mean_y)) +
  geom_errorbar(aes(
    ymin = mean_y - sd_y,
    ymax = mean_y + sd_y,
    color = group
  ), alpha = 0.6) +
  geom_point(aes(color = group)) +
  geom_line(
    data = sp_true_means,
    mapping = aes(
      x = t,
      y = mean_true,
      group = group,
      color = group
    ),
    size = 1.5,
    alpha = 0.5
  ) +
  labs(
    title = "Cross-sectional means plus/minus one cross-sectional \nstandard deviation",
    y = "y",
    x = TeX('$t$'),
    y = TeX('$y_{g,i}(t)$'),
    color = "Group",
    fill = "Group"
  ) +
  scale_color_manual(values = col_2groups) +
  scale_fill_manual(values = col_2groups) +
  theme_bw(base_family = ff, base_size = fS2_fs) +
  theme(legend.position = "bottom", aspect.ratio = sp_ar) +
  ylim(-0.15, 1.4)

# Plot the PFI from the cross-sectional approach (using accuracy as the metric)
plot_cs_pfi_acc <-
  data.frame(
    time = sort(unique(sp_data$t)),
    imp = sp_pfi_acc_cs$pfi
  ) |>
  ggplot(aes(x = time, y = imp)) +
  geom_point() +
  geom_segment(aes(xend = time, yend = 0)) +
  labs(
    title = "Feature importance for cross-sectional model \n(computed using log-loss)",
    x = TeX('$t$'),
    y = "PFI"
  ) +
  theme_bw(base_family = ff, base_size = fS2_fs) + 
  theme(aspect.ratio = sp_ar)

# Plot the PFI from the cross-sectional approach (using log-loss as the metric)
plot_cs_pfi_ll <-
  data.frame(
    time = sort(unique(sp_data$t)),
    imp = sp_pfi_ll_cs$pfi
  ) |>
  ggplot(aes(x = time, y = imp)) +
  geom_point() +
  geom_segment(aes(xend = time, yend = 0)) +
  labs(
    title = "Feature importance for cross-sectional model \n(computed using log-loss)",
    x = TeX('$t$'),
    y = "PFI"
  ) +
  theme_bw(base_family = ff, base_size = fS2_fs) + 
  theme(aspect.ratio = sp_ar)

# Plot the PFI from the VEESA pipeline (using log-loss as the metric)
plot_veesa_pfi_ll <-
  data.frame(fpc = 1:length(sp_pfi_ll_jfpca$pfi),
             imp = sp_pfi_ll_jfpca$pfi) |>
  ggplot(aes(x = fpc, y = imp)) +
  geom_point() +
  geom_segment(aes(xend = fpc, yend = 0)) +
  labs(title = "Feature importance for VEESA pipeline model \n(computed using log-loss)",
       x = TeX('jfPC'),
       y = "PFI") +
  theme_bw(base_family = ff, base_size = fS2_fs) + 
  theme(aspect.ratio = sp_ar)# + 
  #ylim(0, 0.15)

# Extract the legend
legend <- 
  get_plot_component(
    plot_sp_cs_means, 
    'guide-box-bottom', 
    return_all = TRUE
  )

# Join the plots
plot_grid(
  plot_grid(
    plot_sp_cs_means + 
      theme(legend.position = "none"),
    plot_cs_pfi_acc,
    plot_cs_pfi_ll,
    plot_veesa_pfi_ll,
    ncol = 2
  ),
  legend,
  ncol = 1,
  rel_heights = c(0.95, 0.05)
)
```

```{r figS3, fig.width = 16, fig.height = 14, out.width = "5.5in", fig.cap = paste0("Principal direction plots of the VEESA pipeline jfPCs with log-loss based PFI values greater than ", sp_veesa_pfi_cutoff, " (excluding jfPC 1, which is shown in the main text).")}

# Specify the figure size (for determining font size)
fS3_fs = 22

sp_veesa_pc_top <- 
  data.frame(pfi = sp_pfi_ll_jfpca$pfi) |>
  mutate(pc = 1:n()) |>
  filter(pfi > sp_veesa_pfi_cutoff) |>
  arrange(desc(pfi)) |>
  pull(pc)

plot_pc_directions(
  fpc = sp_veesa_pc_top,
  fdasrvf = sp_train_esa$fpca_res,
  fpca_method = "jfpca",
  times = -sp_times,
  alpha = 0.9, 
  nrow = 2, 
  linesizes = c(0.75, 1, 1.25, 1, 0.75)
) +
  theme_bw(base_family = ff, base_size = fS3_fs) +
  theme(
    strip.background = element_rect(color = "white", fill = "white"),
    legend.key.width = unit(1,"cm"),
    aspect.ratio = 0.4,
    legend.position = "bottom",
    axis.title = element_text(size = fS3_fs), 
    axis.text = element_text(size = fS3_fs), 
    title = element_text(size = fS3_fs),
    strip.text = element_text(size = fS3_fs),
    legend.text = element_text(size = fS3_fs),
    legend.title = element_text(size = fS3_fs)
  ) +
  scale_color_manual(values = col_pcdir_2sd) +
  guides(color = guide_legend(override.aes = list(linewidth = 1.25))) +
  labs(x = TeX('$t$'), y = TeX('$y_g(t)$'))
```

# H-CT Material Classification Additional Analyses {#hct}

```{r hct-load-res}
# Load subset of training data
hct_train_sub = py_load_object(paste0(fp, "data/hct-sub-train.pkl"))

# Load principal components from training data
hct_train_vfpca_pcs <- 
  py_load_object(paste0(fp, "results/hct/hct-train-vfpca-pc-dirs-sparam20.pkl"))

hct_train_jfpca_pcs <- 
  py_load_object(
    paste0(fp, "results/hct/hct-train-jfpca-centered-pc-dirs-sparam20.pkl")
  )

# Load latent values from training data
hct_train_vfpca_latent <- py_load_object(
  paste0(fp, "results/hct/hct-train-vfpca-latent-sparam20.pkl")
)
hct_train_jfpca_latent <- py_load_object(
  paste0(fp, "results/hct/hct-train-jfpca-latent-sparam20.pkl")
)

# Load PFI results
hct_jfpca_nn_pfi <- py_load_object(paste0(
  fp, "results/hct/hct-test-pfi-jfpca-sparam20.pkl"
))
hct_vfpca_nn_pfi <- py_load_object(paste0(
  fp, "results/hct/hct-test-pfi-vfpca-sparam20.pkl"
))

# Vector of smoothing parameter values
sparams = c("01", "05", "10", "15", "20", "25")

# Function for loading predictions and accuracies
load_results <- function(train_test, result, stage, sparam) {
  fp = paste0(
    fp, "results/hct/hct-", train_test, "-", 
    result, "-", stage, "-sparam", sparam, ".pkl"
  )
  py_load_object(fp)
}

# Load predictions and accuracies
hct_pred_and_metrics_train <- 
  set_names(c("jfpca", "vfpca", "hfpca")) |>
  map(.f = function(x) set_names(sparams) |>
  map(
    .f = load_results, 
    stage = x, 
    train_test = "train", 
    result = "pred-and-metrics")
  )

hct_pred_and_metrics_test <- 
  set_names(c("jfpca", "vfpca", "hfpca")) |>
  map(.f = function(x) set_names(sparams) |>
  map(
    .f = load_results, 
    stage = x, 
    train_test = "test", 
    result = "pred-and-metrics")
)

hct_pred_and_metrics_cs_pre_smooth <- 
  set_names(c("train", "test")) |>
  map(
    .f = load_results, 
    stage = "cs-pre-smooth",
    result = "pred-and-metrics", 
    sparam = "00"
  )

hct_pred_and_metrics_cs_pre_align <- 
  set_names(c("train", "test")) |>
  map(.f = function(x) set_names(sparams) |>
  map(
    .f = load_results, 
    stage = "cs-pre-align", 
    train_test = x, 
    result = "pred-and-metrics")
)

hct_pred_and_metrics_cs_post_align <- 
  set_names(c("train", "test")) |>
  map(.f = function(x) set_names(sparams) |>
  map(
    .f = load_results, 
    stage = "cs-post-align", 
    train_test = x, 
    result = "pred-and-metrics")
  )
```

```{r hct-cv-metrics}

# Prepare metrics for plotting 
hct_pred_and_metrics_train_df <-
  map_df(
    .x = hct_pred_and_metrics_train,
    .f = function(x)
      map_df(x, "acc"),
    .id = "pca"
  ) |>
  mutate(data = "training")

hct_pred_and_metrics_test_df <-
  map_df(
    .x = hct_pred_and_metrics_test,
    .f = function(x)
      map_df(x, "acc"),
    .id = "pca"
  ) |> 
  mutate(data = "testing")

hct_pred_and_metrics_cs_pre_smooth_df <-
  map_df(
    .x = hct_pred_and_metrics_cs_pre_smooth,
    .f = "acc"
  ) |> 
  mutate(pca = "CS pre-smooth") |>
  pivot_longer(cols = -pca, names_to = "data", values_to = "00") |>
  mutate(data = ifelse(data == "train", "training", "testing"))

hct_pred_and_metrics_cs_pre_align_df <-
  map_df(
    .x = hct_pred_and_metrics_cs_pre_align,
    .f = function(x)
      map_df(x, "acc"),
    .id = "data"
  ) |> 
  mutate(data = ifelse(data == "train", "training", "testing")) |>
  mutate(pca = "CS pre-align")

hct_pred_and_metrics_cs_post_align_df <-
  map_df(
    .x = hct_pred_and_metrics_cs_post_align,
    .f = function(x)
      map_df(x, "acc"),
    .id = "data"
  ) |> 
  mutate(data = ifelse(data == "train", "training", "testing")) |>
  mutate(pca = "CS post-align")
```

```{r hct-pfi-pcs}
# Function for identifying a certain number of top jfPCs
get_top_vars <- function(pfi, nvars) {
  data.frame(pfi_mean = pfi$importances_mean) |>
    mutate(x = 1:n()) |>
    arrange(desc(pfi_mean)) |>
    slice(1:nvars) |>
    pull(x)
}

# Specify number of top fPCs to consider 
hct_n_top_pcs = 9

# Identify the top fPCs
hct_top_jfpcs = get_top_vars(pfi = hct_jfpca_nn_pfi, nvars = hct_n_top_pcs)
hct_top_vfpcs = get_top_vars(pfi = hct_vfpca_nn_pfi, nvars = hct_n_top_pcs)

# Join PFI for jfPCA and vfPCA
hct_pfi <-
  bind_rows(
    data.frame(hct_jfpca_nn_pfi$importances) |>
      mutate(fpca = "Joint fPCA", fpc = 1:n()),
    data.frame(hct_vfpca_nn_pfi$importances) |>
      mutate(fpca = "Vertical fPCA", fpc = 1:n())
  ) |>
  pivot_longer(
    cols = X1:X5,
    values_to = "fi",
    names_to = "rep"
  ) |>
  mutate(rep = str_remove(rep, "X"))
```

## Smoothing Selection and Cross-Sectional Comparison {#hct-smoothing}

```{r figS4, fig.width = 13, fig.height = 7, out.width = '5.5in', fig.cap = "Model accuracies from neural networks applied using the VEESA pipeline and the cross-sectional approach. "}

# Specify the figure size (for determining font size)
fS4_fs = 20

# Create plot
bind_rows(
  hct_pred_and_metrics_train_df,
  hct_pred_and_metrics_test_df,
  hct_pred_and_metrics_cs_pre_smooth_df,
  hct_pred_and_metrics_cs_pre_align_df,
  hct_pred_and_metrics_cs_post_align_df
) |>
  pivot_longer(cols = -c(data, pca),
               names_to = "sparam",
               values_to = "acc") |>
  mutate(sparam = as.numeric(sparam)) |>
  mutate(pca = forcats::fct_relevel(
    pca,
    c(
      "CS pre-smooth",
      "CS pre-align",
      "CS post-align",
      "jfpca",
      "vfpca",
      "hfpca"
    )
  )) |>
  ggplot(aes(
    x = sparam,
    y = acc,
    color = pca,
    linetype = data,
    shape = data,
    group = factor(data):factor(pca)
  )) +
  geom_point(size = 3) +
  geom_line(size = 1) +
  theme_bw(base_family = ff, base_size = fS4_fs) +
  theme(legend.key.width = unit(2, "cm")) +
  gretchenalbrecht::scale_color_gretchenalbrecht(palette = "last_rays") +
  scale_linetype_manual(values = c("dashed", "solid")) +
  labs(
    x = "Number of Box-Filter Runs",
    y = "Accuracy",
    color = "Method",
    linetype = "Dataset",
    shape = "Dataset"
  )
```

## Analysis Times {#hct-times}

```{r hct-times}
hct_time_files = list.files(paste0(fp, "results/hct/"), pattern = "time")

hct_times_raw <-
  map(
    .x = hct_time_files, 
    .f = function(file) {
      data.frame(time = py_load_object(paste0(fp, "results/hct/", file))) |>
        mutate(file = file)
    }
  ) |> list_rbind() |>
    select(file, time)

hct_times <-
  hct_times_raw |>
  mutate(data = ifelse(
    str_detect(file, "both"),
    "train_and_test",
    ifelse(
      str_detect(file, "train"),
      "train", 
    ifelse(
      str_detect(file, "noisy-test"), 
      "noisy-test", 
      "test"
    )))) |> 
  mutate(
    file = str_remove(file, ".pkl"),
    file = str_remove(file, "hct-"),
    file = str_remove(file, "time-"),
    file = str_remove(file, "both-"),
    file = str_remove(file, "test-"),
    file = str_remove(file, "train-")
  ) |> 
  mutate(sparam = file) |> 
  mutate(
    sparam = str_remove(sparam, "split"),
    sparam = str_remove(sparam, "aligned-"),
    sparam = str_remove(sparam, "jfpca-"),
    sparam = str_remove(sparam, "vfpca-"),
    sparam = str_remove(sparam, "hfpca-"),
    sparam = str_remove(sparam, "pred-and-metrics-"),
    sparam = str_remove(sparam, "smoothed-"),
    sparam = str_remove(sparam, "pfi-"),
    sparam = str_remove(sparam, "model-"),
    sparam = str_remove(sparam, "sparam")
  ) |>
  mutate(sparam = as.numeric(sparam)) |> 
  mutate(
    pipeline_step = file,
    pipeline_step = str_remove(pipeline_step, "-sparam"),
    pipeline_step = str_remove(pipeline_step, "01"),
    pipeline_step = str_remove(pipeline_step, "05"),
    pipeline_step = str_remove(pipeline_step, "10"),
    pipeline_step = str_remove(pipeline_step, "15"),
    pipeline_step = str_remove(pipeline_step, "20"),
    pipeline_step = str_remove(pipeline_step, "25")
  ) |> 
  mutate(
    pipeline_step_clean = pipeline_step,
    pipeline_step_clean = str_replace(pipeline_step_clean, "aligned-hfpca", "Alignment and PCA"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "aligned-jfpca", "Alignment and PCA"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "aligned-vfpca", "Alignment and PCA"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "aligned", "Alignment"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "model-hfpca", "Model"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "model-vfpca", "Model"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "model-jfpca", "Model"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "pred-and-metrics-hfpca", "Predictions and Metrics"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "pred-and-metrics-jfpca", "Predictions and Metrics"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "pred-and-metrics-vfpca", "Predictions and Metrics"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "pfi-hfpca", "PFI"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "pfi-vfpca", "PFI"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "pfi-jfpca", "PFI"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "smoothed", "Smoothing"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "split", "Train/Test Split"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "hfpca", "PCA"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "jfpca", "PCA"),
    pipeline_step_clean = str_replace(pipeline_step_clean, "vfpca", "PCA")
  ) |>
  mutate(
    pipeline_step_clean = fct_relevel(
      pipeline_step_clean,
      "Train/Test Split",
      "Smoothing",
      "Alignment",
      "PCA",
      "Alignment and PCA",
      "Model",
      "Predictions and Metrics",
      "PFI"
    )
  ) |>
  mutate(
    pca = pipeline_step,
    pca = str_remove(pca, "split"),
    pca = str_remove(pca, "aligned-"),
    pca = str_remove(pca, "pfi-"),
    pca = str_remove(pca, "pred-and-metrics-"),
    pca = str_remove(pca, "smoothed"),
    pca = str_remove(pca, "model-"),
    pca = str_remove(pca, "aligned"),
    pca = ifelse(pca=="", "other-step", pca)
  ) |>
  mutate(hours = round(time / 60 / 60, 2), mins = round(time / 60, 2)) |> 
  select(
    data,
    pipeline_step,
    pipeline_step_clean,
    pca,
    sparam,
    time,
    hours,
    mins,
    file
  )

hct_times_ave <-
  hct_times |>
  summarise(time = mean(time), .by = c(data, pipeline_step_clean, pca)) |> 
  mutate(hours = round(time / 60 / 60, 2)) |>
  mutate(mins = round(time / 60, 2))
```

```{r figS5, fig.height = 5, fig.width = 12, out.width = "5.5in", fig.cap = "Analysis step times for the VEESA pipeline analyses of the H-CT data."}

# Specify the figure font size
fS5_fs = 16

hct_times |>
  filter(
    sparam == 1,
    data != "noisy-test",
    pipeline_step_clean != "Train/Test Split"
  ) |>
  ggplot(aes(
    y =  fct_rev(pipeline_step_clean), 
    x = hours, 
    color = data, 
    shape = data
  )) + 
  geom_point(size = 4) +
  scale_color_manual(values = col_5groups[c(3,1)]) +
  scale_shape_manual(values = c(17,19)) +
  theme_bw(base_size = fS5_fs, base_family = ff) +
  theme(
    strip.background = element_blank(),
    axis.title = element_text(size = fS5_fs),
    axis.text = element_text(size = fS5_fs),
    title = element_text(size = fS5_fs),
    strip.text = element_text(size = fS5_fs),
    legend.text = element_text(size = fS5_fs),
    legend.title = element_text(size = fS5_fs)
  ) +
  labs(
    y = "Pipeline Step",
    x = "Hours",
    color = "Data",
    shape = "Data"
  )
```

```{r figS6, fig.height = 5, fig.width = 12, out.width = "5.5in", fig.cap = "Analysis step times for the VEESA pipeline analyses of the H-CT data (with the alignment of the training data removed)."}

hct_times |>
  filter(
    sparam == 1,
    data != "noisy-test", 
    pipeline_step_clean != "Alignment",
    pipeline_step_clean != "Train/Test Split"
  ) |>
  ggplot(aes(
    y =  fct_rev(pipeline_step_clean), 
    x = hours, 
    color = pca, 
    shape = pca
  )) +
  geom_point(size = 4) +
  scale_color_manual(values = col_5groups[c(2,1,4,5)]) +
  scale_shape_manual(values = c(16,17,19,18)) +
  theme_bw(base_size = fS5_fs, base_family = ff) +
  theme(
    strip.background = element_blank(),
    axis.title = element_text(size = fS5_fs),
    axis.text = element_text(size = fS5_fs),
    title = element_text(size = fS5_fs),
    strip.text = element_text(size = fS5_fs),
    legend.text = element_text(size = fS5_fs),
    legend.title = element_text(size = fS5_fs)
  ) +
  labs(
    y = "Pipeline Step",
    x = "Hours",
    color = "PCA",
    shape = "PCA"
  )
```

```{r figS7, fig.height = 5, fig.width = 12, out.width = "5.5in", fig.cap = "Analysis step times for the VEESA pipeline analyses of the H-CT data (with the alignment and PCA computations of the training and testing data removed)."}

hct_times |>
  filter(
    sparam == 1,
    data != "noisy-test", 
    pipeline_step_clean != "Alignment",
    pipeline_step_clean != "Train/Test Split",
    pipeline_step_clean != "PCA",
    pipeline_step_clean != "Alignment and PCA"
  ) |>
  ggplot(aes(
    y =  fct_rev(pipeline_step_clean), 
    x = hours, 
    color = pca, 
    shape = pca
  )) +
  geom_point(size = 4) +
  scale_color_manual(values = col_5groups[c(2,1,4,5)]) +
  scale_shape_manual(values = c(16,17,19,18)) +
  theme_bw(base_size = fS5_fs, base_family = ff) +
  theme(
    strip.background = element_blank(),
    axis.title = element_text(size = fS5_fs),
    axis.text = element_text(size = fS5_fs),
    title = element_text(size = fS5_fs),
    strip.text = element_text(size = fS5_fs),
    legend.text = element_text(size = fS5_fs),
    legend.title = element_text(size = fS5_fs)
  ) +
  labs(
    y = "Pipeline Step",
    x = "Hours",
    color = "PCA",
    shape = "PCA"
  )
```

## PFI Variability {#hct-pfi-var}

```{r figS8, fig.width = 12, fig.height = 6, out.width = '5in', fig.cap = "PFI replicate values for vfPCs associated with the H-CT material classification example.", warning = FALSE}

# Specify the figure font size
fS8_fs = 18

# Create plots of proportion of variation explained and PFI values
hct_pfi |>
  filter(fpca == "Vertical fPCA") |>
  ggplot(aes(x = fpc, y = fi, color = rep)) +
  geom_point(size = 2, alpha = 0.5) +
  facet_grid(fpca ~ ., switch = "y") +
  labs(x = "Principal Component", color = "PFI replicate") +
  theme_bw(base_family = ff, base_size = fS8_fs) +
  theme(
    strip.placement = "outside",
    strip.background = 
      element_rect(color = "white", fill = "white"),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = fS8_fs),
    axis.text = element_text(size = fS8_fs),
    title = element_text(size = fS8_fs),
    strip.text = element_text(size = fS8_fs),
    legend.position = "bottom"
  ) + 
  scale_color_manual(values = wes_palettes$Zissou1[5:1])
```

## vfPCA Principal Directions (Additional) {#hct-vfPCs}

```{r figS9, fig.width = 15, fig.height = 3.5, out.width = '6.5in', fig.cap = "Proportion of variation and PFI values associated with vfPCs. The vfPCs with the nine highest PFI values are labeled and colored.", warning = FALSE}

# Specify the figure font size
fS9_fs = 16

# Get vfPCA prop vars
hct_prop_var_vfpca <-
  data.frame(latent = hct_train_vfpca_latent) |>
  mutate(index = 1:n(),
         prop_var = (latent^2) / (sum(latent^2)))

# Prepare proportion of variability and PFI for vfPCA for plotting 
hct_pv_pfi <-
  hct_prop_var_vfpca |>
  mutate(fpca = "Joint fPCA", pfi = hct_vfpca_nn_pfi$importances_mean) |>
  select(-latent) |>
  rename(fpc = index) |>
  pivot_longer(cols = -c(fpc, fpca), names_to = "variable") |>
  mutate(
    variable = fct_relevel(variable, "prop_var", "pfi"),
    variable = fct_recode(
      variable,
      "Proportion of \nVariation" = "prop_var",
      "PFI" = "pfi"
    )
  )

# Prepare labels
hct_pfi_labels <-
  hct_pv_pfi |>
  filter(variable == "PFI") |>
  filter(fpc %in% hct_top_vfpcs) |>
  mutate(label_x = fpc, label_y = value+ 0.02, label_id = 1:n())

# Create plots of proportion of variation explained and PFI values
left_join(
  hct_pv_pfi,
  hct_pfi_labels, 
  by = join_by(fpc, fpca, variable, value)
) |>
  mutate(has_label = ifelse(is.na(label_x), F, T)) |>
  ggplot(aes(x = fpc, y = value)) +
  geom_segment(aes(xend = fpc, yend = 0, color = has_label), size = 0.75) +
  facet_wrap(. ~ variable, scales = "free_y", switch = "y") +
  geom_text(aes(x = label_x, y = label_y, label = fpc, color = has_label),
            size = fS9_fs / 3) +
  labs(x = "Principal Component") +
  scale_color_manual(values = c("black", col_5groups[1])) +
  theme_bw(base_family = ff, base_size = fS9_fs) +
  theme(
    strip.placement = "outside",
    strip.background =
      element_rect(color = "white", fill = "white"),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = fS9_fs),
    axis.text = element_text(size = fS9_fs),
    title = element_text(size = fS9_fs),
    strip.text = element_text(size = fS9_fs), 
    legend.position = "none"
  )
```

```{r figS10, fig.width = 16, fig.height = 9, out.width = "6.5in", fig.cap = "Principal directions from the nine vfPCs with the highest PFI from the H-CT material classification example."}

# Specify the figure font size
fS10_fs = 18.5

plot_pc_directions(
  fpcs = hct_top_vfpcs,
  fdasrvf = list(
    "f_pca" = hct_train_vfpca_pcs[, 2:6, ],
    "latent" = hct_train_vfpca_latent
  ),
  times = unique(hct_train_sub$frequency_norm),
  fpca_method = "vfpca",
  alpha = 0.9,
  nrow = 3,
  linesizes = rep(0.75, 5), 
  linetype = FALSE
) +
  theme_bw(base_size = fS10_fs, base_family = ff) +
  theme(
    legend.key.width = unit(1, "cm"),
    strip.background = element_rect(color = "white", fill = "white"),
    aspect.ratio = 0.4,
    legend.position = "bottom",
    axis.title = element_text(size = fS10_fs), 
    axis.text = element_text(size = fS10_fs), 
    title = element_text(size = fS10_fs),
    strip.text = element_text(size = fS10_fs),
    legend.text = element_text(size = fS10_fs),
    legend.title = element_text(size = fS10_fs)
  )  +
  scale_color_manual(values = col_pcdir_2sd) +
  guides(color = guide_legend(override.aes = list(linewidth = 1.25))) +
  labs(
    x = "Normalized Frequency",
    y = "Intensity",
    color = "Material",
    size = "Material",
    linetype = "Material"
  )
```

```{r figS11, fig.width = 16, fig.height = 9, out.width = "6.5in", fig.cap = "Principal differences from the nine vfPCs with the highest PFI from the H-CT material classification example."}

plot_pc_diffs(
  fpcs = hct_top_vfpcs,
  fdasrvf = list(
    "f_pca" = hct_train_vfpca_pcs[, 2:6, ],
    "latent" = hct_train_vfpca_latent
  ),
  times = unique(hct_train_sub$frequency_norm),
  fpca_method = "vfpca",
  alpha = 0.9,
  alpha_fill = 0.2,
  nrow = 3,
  linesizes = rep(0.75, 5), 
  linetype = FALSE
) +
  theme_bw(base_size = fS10_fs, base_family = ff) +
  theme(
    legend.key.width = unit(1, "cm"),
    strip.background = element_rect(color = "white", fill = "white"),
    aspect.ratio = 0.4,
    legend.position = "bottom",
    axis.title = element_text(size = fS10_fs), 
    axis.text = element_text(size = fS10_fs), 
    title = element_text(size = fS10_fs),
    strip.text = element_text(size = fS10_fs),
    legend.text = element_text(size = fS10_fs),
    legend.title = element_text(size = fS10_fs)
  )  +
  scale_color_manual(values = col_pcdir_2sd[c(1,2,4,5)]) +
  scale_fill_manual(values = col_pcdir_2sd[c(1,2,4,5)]) +
  guides(color = guide_legend(override.aes = list(linewidth = 1.25))) +
  labs(
    x = "Normalized Frequency",
    y = "Intensity",
    color = "Material",
    size = "Material",
    fill = "Material",
    linetype = "Material"
  )
```

## jfPCA Principal Directions {#hct-jfPCs}

```{r figS12, fig.width = 15, fig.height = 3.5, out.width = '6.5in', fig.cap = "Proportion of variation and PFI values associated with jfPCs. The jfPCs with the nine highest PFI values are labeled and colored.", warning = FALSE}

# Get jfPCA prop vars
hct_prop_var_jfpca <-
  data.frame(latent = hct_train_jfpca_latent) |>
  mutate(index = 1:n(),
         prop_var = (latent^2) / (sum(latent^2)))

# Prepare proportion of variability and PFI for jfPCA for plotting 
hct_pv_pfi <-
  hct_prop_var_jfpca |>
  mutate(fpca = "Joint fPCA", pfi = hct_jfpca_nn_pfi$importances_mean) |>
  select(-latent) |>
  rename(fpc = index) |>
  pivot_longer(cols = -c(fpc, fpca), names_to = "variable") |>
  mutate(
    variable = fct_relevel(variable, "prop_var", "pfi"),
    variable = fct_recode(
      variable,
      "Proportion of \nVariation" = "prop_var",
      "PFI" = "pfi"
    )
  )

# Prepare labels
hct_pfi_labels <-
  hct_pv_pfi |>
  filter(variable == "PFI") |>
  filter(fpc %in% hct_top_jfpcs) |>
  mutate(label_x = fpc, label_y = value, label_id = 1:n()) |>
  mutate(label_x = ifelse(label_id == 1, label_x - 3, label_x)) |>
  mutate(label_y = ifelse(label_id == 2, label_y + 0.02, label_y)) |>
  mutate(label_x = ifelse(label_id == 2, label_x - 2, label_x)) |>
  mutate(label_y = ifelse(label_id == 3, label_y + 0.01, label_y)) |>
  mutate(label_x = ifelse(label_id == 3, label_x + 3.5, label_x)) |>
  mutate(label_y = ifelse(label_id == 4, label_y + 0.01, label_y)) |>
  mutate(label_y = ifelse(label_id == 5, label_y + 0.01, label_y)) |>
  mutate(label_y = ifelse(label_id == 6, label_y + 0.01, label_y))

# Create plots of proportion of variation explained and PFI values
left_join(
  hct_pv_pfi,
  hct_pfi_labels, 
  by = join_by(fpc, fpca, variable, value)
) |>
  mutate(has_label = ifelse(is.na(label_x), F, T)) |>
  ggplot(aes(x = fpc, y = value)) +
  geom_segment(aes(xend = fpc, yend = 0, color = has_label), size = 0.75) +
  facet_wrap(. ~ variable, scales = "free_y", switch = "y") +
  geom_text(aes(x = label_x, y = label_y, label = fpc, color = has_label),
            size = fS9_fs / 3) +
  labs(x = "Principal Component") +
  scale_color_manual(values = c("black", col_5groups[1])) +
  theme_bw(base_family = ff, base_size = fS9_fs) +
  theme(
    strip.placement = "outside",
    strip.background =
      element_rect(color = "white", fill = "white"),
    axis.title.y = element_blank(),
    axis.title.x = element_text(size = fS9_fs),
    axis.text = element_text(size = fS9_fs),
    title = element_text(size = fS9_fs),
    strip.text = element_text(size = fS9_fs), 
    legend.position = "none"
  )
```

```{r figS13, fig.width = 16, fig.height = 9, out.width = "6.5in", fig.cap = "Principal directions from the nine jfPCs with the highest PFI from the H-CT material classification example."}

plot_pc_directions(
  fpcs = hct_top_jfpcs,
  fdasrvf = list(
    "f_pca" = hct_train_jfpca_pcs[, 2:6, ],
    "latent" = hct_train_jfpca_latent
  ),
  times = unique(hct_train_sub$frequency_norm),
  fpca_method = "jfpca",
  alpha = 0.9,
  nrow = 3,
  linesizes = rep(0.75, 5), 
  linetype = FALSE
) +
  theme_bw(base_size = fS10_fs, base_family = ff) +
  theme(
    legend.key.width = unit(1, "cm"),
    strip.background = element_rect(color = "white", fill = "white"),
    aspect.ratio = 0.4,
    legend.position = "bottom",
    axis.title = element_text(size = fS10_fs), 
    axis.text = element_text(size = fS10_fs), 
    title = element_text(size = fS10_fs),
    strip.text = element_text(size = fS10_fs),
    legend.text = element_text(size = fS10_fs),
    legend.title = element_text(size = fS10_fs)
  )  +
  scale_color_manual(values = col_pcdir_2sd) +
  guides(color = guide_legend(override.aes = list(linewidth = 1.25))) +
  labs(
    x = "Normalized Frequency",
    y = "Intensity",
    color = "Material",
    size = "Material",
    linetype = "Material"
  )
```

```{r figS14, fig.width = 16, fig.height = 9, out.width = "6.5in", fig.cap = "Principal differences from the nine jfPCs with the highest PFI from the H-CT material classification example."}

plot_pc_diffs(
  fpcs = hct_top_jfpcs,
  fdasrvf = list(
    "f_pca" = hct_train_jfpca_pcs[, 2:6, ],
    "latent" = hct_train_jfpca_latent
  ),
  times = unique(hct_train_sub$frequency_norm),
  fpca_method = "jfpca",
  alpha = 0.9,
  alpha_fill = 0.2,
  nrow = 3,
  linesizes = rep(0.75, 5), 
  linetype = FALSE
) +
  theme_bw(base_size = fS10_fs, base_family = ff) +
  theme(
    legend.key.width = unit(1, "cm"),
    strip.background = element_rect(color = "white", fill = "white"),
    aspect.ratio = 0.4,
    legend.position = "bottom",
    axis.title = element_text(size = fS10_fs), 
    axis.text = element_text(size = fS10_fs), 
    title = element_text(size = fS10_fs),
    strip.text = element_text(size = fS10_fs),
    legend.text = element_text(size = fS10_fs),
    legend.title = element_text(size = fS10_fs)
  )  +
  scale_color_manual(values = col_pcdir_2sd[c(1,2,4,5)]) +
  scale_fill_manual(values = col_pcdir_2sd[c(1,2,4,5)]) +
  guides(
    color = guide_legend(override.aes = list(linewidth = 1.25)),
    fill = guide_legend(override.aes = list(linewidth = 1.25))
  ) +
  labs(
    x = "Normalized Frequency",
    y = "Intensity",
    color = "Material",
    fill = "Material",
    size = "Material",
    linetype = "Material"
  )
```

# Inkjet Printer Additional Analyses {#inkjet}

```{r inkjet-load-data-and-res}
# Load inkjet data
inkjet = read.csv(paste0(fp, "data/inkjet-cleaned.csv"))

# Load cross-validation folds 
inkjet_folds = read.csv(paste0(fp, "data/inkjet-cv-folds.csv"))

# Load CV results
inkjet_cv_preds = readRDS(paste0(fp_res_inkjet, "inkjet-cv-preds.rds"))

# Load CV results from Buzzini paper
inkjet_paper_res_best <-
  read.csv(paste0(fp_res_inkjet, "inkjet-paper-res-best.csv"))

# Load cross-validation results
inkjet_cv_acc_summary <-
  read.csv(paste0(fp_res_inkjet, "inkjet-cv-accuracy-summary.csv"))

# Load PFI results
inkjet_pfi_best <- py_load_object(
  paste0(fp_res_inkjet, "inkjet-pfi-best.csv"), 
  convert = TRUE
)
inkjet_pfi_worst <- py_load_object(
  paste0(fp_res_inkjet, "inkjet-pfi-worst.pkl"), 
  convert = TRUE
)

# Load CV scenarios
inkjet_smooth_rep_fold_color_list <- 
  readRDS(paste0(fp_res_inkjet, "inkjet-cv-srfc.rds"))

# Best and worst results
inkjet_res_best_fp = paste0(fp_res_inkjet, "inkjet-cv-res-best.csv")
inkjet_res_worst_fp = paste0(fp_res_inkjet, "inkjet-cv-res-worst.csv")
inkjet_res_best = read.csv(inkjet_res_best_fp)
inkjet_res_worst = read.csv(inkjet_res_worst_fp)
```

```{r inkjet-improve-worst}
# Specify a file path
inkjet_cv_preds_pfi_informed_fp <-
  paste0(fp_res_inkjet, "inkjet-cv-preds-pfi-informed.rds")

# Train models with best performing PCs
if (!file.exists(inkjet_cv_preds_pfi_informed_fp)) {

  # Function for loading inkjet CV objects
  load_inkjet <- function(rep, test_fold, color, data_type, fp, smoothing) {
    py_load_object(
      filename = paste0(
        fp_res_inkjet,
        "inkjet-cv-s",
        smoothing,
        "-rep_",
        rep,
        "-fold_",
        test_fold,
        "-color_",
        color,
        "-",
        data_type,
        ".pkl"
      )
    )
  }
  
  # Function for implementing cross-validation for inkjet data
  get_inkjet_rf_cv_preds_pfi_informed <- function(s, rep_curr, test_fold, color) {
  
    # Load training data
    train <-
      load_inkjet(
        smoothing = s, 
        rep = rep_curr, 
        test_fold = test_fold, 
        color = color, 
        data_type = "jfpca-train", 
        fp = fp
      )
    
    # Load testing data
    test <-
      load_inkjet(
        smoothing = s, 
        rep = rep_curr, 
        test_fold = test_fold, 
        color = color, 
        data_type = "aligned-jfpca-test",
        fp = fp
      )
    
    # Extract training data printer and sample values
    train_printer_sample <-
      inkjet_folds |> 
      filter(rep == rep_curr, fold != test_fold)
      
    # Extract testing data printer and sample values
    test_printer_sample <-
      inkjet_folds |>
      filter(rep == rep_curr, fold == test_fold)

    # Create a dataframe with training data
    rf_data <- data.frame(
      printer = factor(train_printer_sample$printer), 
      train$coef[,c(1:25, 75:100)]
    )
    
    # Train random forest
    set.seed(20230630)
    rf = randomForest(printer ~ ., data = rf_data, ntree = 500)
  
    # Get test fold predictions
    preds = predict(rf, data.frame(test$coef[,c(1:25, 75:100)]))
    
    # Put predictions in a dataframe
    preds_df <-
      data.frame(
        color = color,
        s = s,
        rep = rep_curr,
        test_fold = test_fold,
        printer = test_printer_sample$printer,
        sample = test_printer_sample$sample, 
        pred = as.numeric(as.vector(preds))
      )
    
    # Return the results
    return(preds_df)
      
  }
  
  # Implement CV for all specified scenarios
  inkjet_cv_preds_pfi_informed <- 
    pmap_df(
      .l = inkjet_smooth_rep_fold_color_list,
      .f = get_inkjet_rf_cv_preds_pfi_informed
    )
  
  # Save PFI informed model predictions
  saveRDS(inkjet_cv_preds_pfi_informed, inkjet_cv_preds_pfi_informed_fp)
  
}

# Load results
inkjet_cv_preds_pfi_informed = readRDS(inkjet_cv_preds_pfi_informed_fp)

# Compute test fold accuracy
inkjet_cv_acc_pfi_informed <-
  inkjet_cv_preds_pfi_informed |>
  summarise(
    acc = sum(pred == printer) / length(printer),
    .by = c(s, rep, test_fold, color)
  )

# Compute summary values of CV accuracies
inkjet_cv_acc_pfi_informed_summary <-
  inkjet_cv_acc_pfi_informed |>
  summarise(
    acc_ave = mean(acc),
    acc_sd = sd(acc),
    acc_min = min(acc),
    acc_max = max(acc), 
    .by = c(s, color)
  ) |>
  mutate(color = factor(color)) |>
  mutate(color = fct_recode(
    color, 
    "Cyan" = "c", 
    "Magenta" = "m", 
    "Yellow" = "y"
  ))
```

```{r inkjet-pc-plots}
# Specify the figure font size
fS17_fs = 28

plot_inkjet_pcs <- function(pfi_res) {
  pc_order <-
    pfi_res$pfi_df |>
    arrange(desc(pfi)) |>
    pull(pc)
  veesa::plot_pc_directions(
    fpc = pc_order[1:5],
    fdasrvf = pfi_res$jfpca,
    fpca_method = "jfpca",
    times = -unique(inkjet$spectra),
    alpha = 0.9,
    nrow = 1,
    linesizes = rep(1.5, 7), 
    linetype = "solid"
  ) +
    scale_x_reverse() +
    theme_bw(base_size = fS17_fs, base_family = ff) +
    theme(
      legend.position = "bottom",
      strip.background = element_rect(color = "white", fill = "white"),
      axis.title = element_text(size = fS17_fs), 
      axis.text = element_text(size = fS17_fs), 
      title = element_text(size = fS17_fs),
      strip.text = element_text(size = fS17_fs),
      legend.text = element_text(size = fS17_fs),
      legend.title = element_text(size = fS17_fs)
    )  +
    scale_color_manual(values = col_pcdir_3sd) +
    guides(
      color = guide_legend(nrow = 1), 
      size = guide_legend(nrow = 1), 
      override.aes = list(alpha = 1, size = 2)
    ) +
    labs(
      x = TeX("Wavenumber (cm$^{-1})$"), 
      y = "Intensity", 
      title = paste(pfi_res$pfi_df$color[1])
    ) 
}

plot_inkjet_pc_diffs <- function(pfi_res) {
  pc_order <-
    pfi_res$pfi_df |>
    arrange(desc(pfi)) |>
    pull(pc)
  plot_pc_diffs(
    fpc = pc_order[1:5],
    fdasrvf = pfi_res$jfpca,
    fpca_method = "jfpca",
    times = -unique(inkjet$spectra),
    alpha = 0.9,
    nrow = 1,
    linesizes = rep(1.5, 6),
    mean_linesize = 1.5,
    linetype = "solid",
    alpha_fill = 0.1
  ) +
    scale_x_reverse() +
    theme_bw(base_size = fS17_fs, base_family = ff) +
    theme(
      legend.position = "bottom",
      strip.background = element_rect(color = "white", fill = "white"),
      axis.title = element_text(size = fS17_fs),
      axis.text = element_text(size = fS17_fs),
      title = element_text(size = fS17_fs),
      strip.text = element_text(size = fS17_fs),
      legend.text = element_text(size = fS17_fs),
      legend.title = element_text(size = fS17_fs)
    )  +
    scale_color_manual(values = col_pcdir_3sd[c(1:3,5:7)]) +
    scale_fill_manual(values = col_pcdir_3sd[c(1:3,5:7)]) +
    guides(
      color = guide_legend(nrow = 1),
      fill = guide_legend(nrow = 1),
      size = guide_legend(nrow = 1),
      override.aes = list(alpha = 1, size = 2)
    ) +
    labs(
      x = TeX("Wavenumber (cm$^{-1})$"),
      y = "Intensity",
      title = paste(pfi_res$pfi_df$color[1])
    )
}
```

## Cross-Validation Procedure Details {#inkjet-cv}

```{r figS15, fig.height = 9.5, fig.width = 22, out.width = "6.5in", fig.cap = "Inkjet Data Cross-Validation Average Accuracies. Triangles pointing up and down highlight the highest and lowest average cross-validation accuracies from the VEESA pipeline, respectively, for each color. Horizontal dashed lines represent best average cross-validation accuracies for each color from Buzzini et al. (2021)."}

# Specify the figure font size
fS15_fs = 26

inkjet_cvs_for_plot <-
  inkjet_cv_acc_summary |>
  mutate(s = factor(s)) |>
  mutate(ntrees = paste(ntrees, "trees")) |>
  mutate(ntrees = factor(
    ntrees, 
    levels = paste(c("50", "100", "250", "500", "1000"), "trees")
  ))

inkjet_res_best_for_plot <-
  inkjet_res_best |>
  mutate(ntrees = paste(ntrees, "trees")) |>
  mutate(ntrees = factor(
    ntrees, 
    levels = paste(c("50", "100", "250", "500", "1000"), "trees")
  ))

inkjet_res_worst_for_plot <-
  inkjet_res_worst |>
  mutate(ntrees = paste(ntrees, "trees")) |>
  mutate(ntrees = factor(
    ntrees, 
    levels = paste(c("50", "100", "250", "500", "1000"), "trees")
  ))

ggplot() +
  geom_hline(
    data = inkjet_paper_res_best,
    mapping = aes(yintercept = acc_ave),
    linewidth = 0.5,
    linetype = "dashed"
  ) +
  geom_line(
    data = inkjet_cvs_for_plot,
    mapping = aes(
      x = pcs,
      y = acc_ave,
      color = s
    )
  ) +
  geom_point(
    data = inkjet_cvs_for_plot,
    mapping = aes(
      x = pcs,
      y = acc_ave,
      color = s
    )
  ) +
  geom_point(
    data = inkjet_res_best_for_plot,
    mapping = aes(
        x = pcs,
        y = acc_ave
      ),
    color = "black",
    fill = "black",
    size = 5,
    shape = 24,
    alpha = 0.5
  ) +
  geom_point(
    data = inkjet_res_worst_for_plot,
    mapping = aes(
        x = pcs,
        y = acc_ave
      ),
    color = "black",
    fill = "black",
    size = 5,
    shape = 25,
    alpha = 0.5
  ) +
  facet_grid(color ~ ntrees) +
  scale_color_manual(
    values = wes_palette(
      name = "Zissou1", 
      n = 16,
      type = "continuous"
    )[c(1:3,8:9,12:15)]) +
  theme_bw(base_size = 20) +
  theme(
    strip.background = element_blank(),
    axis.title = element_text(size = fS15_fs),
    axis.text = element_text(size = fS15_fs),
    title = element_text(size = fS15_fs),
    strip.text = element_text(size = fS15_fs),
    legend.text = element_text(size = fS15_fs),
    legend.title = element_text(size = fS15_fs)
  ) +
  labs(
    x = "Number of PCs",
    y = "Accuracy",
    color = "Number of \nbox-filter \nruns"
  )
```

## Printer Classification Details {#inkjet-conf}

```{r figS16, fig.height = 6, fig.width = 24, out.width = "6.5in", fig.cap = "Confusion matrices for the inkjet printer random forests." }

# Specify the figure font size
fS16_fs = 26

inkjet_plot_conf_cyan <-
  inkjet_cv_preds |>
  filter(
    color == "Cyan",
    s == 20, 
    pcs == 40,
    ntrees == 1000
  ) |>
  count(printer, pred, name = "n_preds") |>
  ggplot(aes(x = printer, y = pred, fill = n_preds)) +
  geom_tile() +
  geom_label(aes(label = n_preds), size = 6) +
  scale_fill_gradient(
    low = "grey95", 
    high = col_inkjet[1]
  ) +
  theme_bw(base_size = fS16_fs, base_family = ff) +
  theme(
    aspect.ratio = 1,
    axis.title = element_text(size = fS16_fs),
    axis.text = element_text(size = fS16_fs),
    title = element_text(size = fS16_fs),
    strip.text = element_text(size = fS16_fs),
    legend.text = element_text(size = fS16_fs),
    legend.title = element_text(size = fS16_fs),
    legend.key.height = unit(1, "cm")
  ) +
  labs(
    x = "Printer (true)",
    y = "Printer (prediction)", 
    fill = "Count"
  ) + 
  scale_y_continuous(breaks = seq(1, 11, by = 1)) + 
  scale_x_continuous(breaks = seq(1, 11, by = 1))

inkjet_plot_conf_magenta <- 
  inkjet_cv_preds |>
  filter(
    color == "Magenta", 
    s == 5, 
    pcs == 40, 
    ntrees == 250
  ) |>
  count(printer, pred, name = "n_preds") |>
  ggplot(aes(x = printer, y = pred, fill = n_preds)) +
  geom_tile() +
  geom_label(aes(label = n_preds), size = 6) +
  scale_fill_gradient(
    low = "grey95", 
    high = col_inkjet[2]
  ) +
  theme_bw(base_size = fS16_fs, base_family = ff) +
  theme(
    aspect.ratio = 1,
    axis.title = element_text(size = fS16_fs), 
    axis.text = element_text(size = fS16_fs), 
    title = element_text(size = fS16_fs),
    strip.text = element_text(size = fS16_fs),
    legend.text = element_text(size = fS16_fs),
    legend.title = element_text(size = fS16_fs),
    legend.key.height = unit(1, "cm")
  ) + 
  labs(
    x = "Printer (true)",
    y = "Printer (prediction)",
    fill = "Count"
  ) +
  scale_y_continuous(breaks = seq(1, 11, by = 1)) + 
  scale_x_continuous(breaks = seq(1, 11, by = 1))

inkjet_plot_conf_yellow <- 
  inkjet_cv_preds |>
  filter(
    color == "Yellow", 
    s == 20, 
    pcs == 20, 
    ntrees == 1000
  ) |>
  count(printer, pred, name = "n_preds") |>
  ggplot(aes(x = printer, y = pred, fill = n_preds)) +
  geom_tile() +
  geom_label(aes(label = n_preds), size = 6) +
  scale_fill_gradient(
    low = "grey95", 
    high = col_inkjet[3]
  ) +
  theme_bw(base_size = fS16_fs, base_family = ff) +
  theme(
    aspect.ratio = 1,
    axis.title = element_text(size = fS16_fs), 
    axis.text = element_text(size = fS16_fs), 
    title = element_text(size = fS16_fs),
    strip.text = element_text(size = fS16_fs),
    legend.text = element_text(size = fS16_fs),
    legend.title = element_text(size = fS16_fs),
    legend.key.height = unit(1, "cm")
  ) + 
  labs(
    x = "Printer (true)",
    y = "Printer (prediction)",
    fill = "Count"
  ) + 
  scale_y_continuous(breaks = seq(1, 11, by = 1)) + 
  scale_x_continuous(breaks = seq(1, 11, by = 1))

plot_grid(
  inkjet_plot_conf_cyan,
  inkjet_plot_conf_magenta,
  inkjet_plot_conf_yellow,
  nrow = 1
)
```

## Additional Principal Directions {#inkjet-jfPCs}

```{r figS17, fig.height = 15, fig.width = 25, out.width = "6in", fig.cap = "The principal direction and principal differences plots for the five most important jfPCs from the worst model for cyan inkjet signatures."}

inkjet_cyan_worst_pc_dirs <-
  plot_inkjet_pcs(inkjet_pfi_worst[[1]]) +
  labs(title = "Cyan (worst directions)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 2)))

inkjet_pfi_worst_for_plot = inkjet_pfi_worst[[1]]
inkjet_pfi_worst_for_plot$jfpca$f_pca = -inkjet_pfi_worst_for_plot$jfpca$f_pca

inkjet_cyan_worst_pc_diffs <-
  plot_inkjet_pc_diffs(inkjet_pfi_worst_for_plot) +
  labs(title = "Cyan (worst differences)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 2)))

plot_grid(
  inkjet_cyan_worst_pc_dirs + theme(legend.position = "none"),
  inkjet_cyan_worst_pc_diffs + theme(legend.position = "none"),
  get_plot_component(
    inkjet_cyan_worst_pc_dirs + 
      theme(legend.key.width = unit(2, 'cm')
    ),
    'guide-box-bottom',
    return_all = TRUE
  ),
  ncol = 1,
  rel_heights = c(0.3, 0.3, 0.05)
)
```

```{r figS18, fig.height = 15, fig.width = 25, out.width = "6in", fig.cap = "The principal direction and principal differences plots for the five most important jfPCs from the best model for magenta inkjet signatures."}

inkjet_magenta_best_pc_dirs <-
  plot_inkjet_pcs(inkjet_pfi_best[[2]]) + 
  labs(title = "Magenta (best directions)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 1.5))) + 
  theme(legend.key.width = unit(3.5, "line"))

inkjet_pfi_best_for_plot_magenta = inkjet_pfi_best[[2]]
inkjet_pfi_best_for_plot_magenta$jfpca$f_pca <-
  inkjet_pfi_best_for_plot_magenta$jfpca$f_pca

inkjet_magenta_best_pc_diffs <-
  plot_inkjet_pc_diffs(inkjet_pfi_best_for_plot_magenta) +
  labs(title = "Magenta (best differences)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 2)))

plot_grid(
  inkjet_magenta_best_pc_dirs + theme(legend.position = "none"),
  inkjet_magenta_best_pc_diffs + theme(legend.position = "none"),
  get_plot_component(
    inkjet_magenta_best_pc_dirs,
    'guide-box-bottom', 
    return_all = TRUE
  ),
  ncol = 1,
  rel_heights = c(0.3, 0.3, 0.05)
)
```

```{r figS19, fig.height = 15, fig.width = 25, out.width = "6in", fig.cap = "The principal direction and principal differences plots for the five most important jfPCs from the worst model for magenta inkjet signatures."}

inkjet_magenta_worst_pc_dirs <-
  plot_inkjet_pcs(inkjet_pfi_worst[[2]]) + 
  labs(title = "Magenta (worst directions)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 1.5))) + 
  theme(legend.key.width = unit(3.5, "line"))

inkjet_pfi_worst_for_plot_magenta = inkjet_pfi_worst[[2]]
inkjet_pfi_worst_for_plot_magenta$jfpca$f_pca <-
  inkjet_pfi_worst_for_plot_magenta$jfpca$f_pca

inkjet_magenta_worst_pc_diffs <-
  plot_inkjet_pc_diffs(inkjet_pfi_worst_for_plot_magenta) +
  labs(title = "Magenta (worst differences)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 2)))

plot_grid(
  inkjet_magenta_worst_pc_dirs + theme(legend.position = "none"),
  inkjet_magenta_worst_pc_diffs + theme(legend.position = "none"),
  get_plot_component(
    inkjet_magenta_worst_pc_dirs,
    'guide-box-bottom', 
    return_all = TRUE
  ),
  ncol = 1,
  rel_heights = c(0.3, 0.3, 0.05)
)
```

```{r figS20, fig.height = 15, fig.width = 25, out.width = "6in", fig.cap = "The principal direction and principal differences plots for the five most important jfPCs from the best model for yellow inkjet signatures."}

inkjet_yellow_best_pc_dirs <-
  plot_inkjet_pcs(inkjet_pfi_best[[3]]) + 
  labs(title = "Yellow (best directions)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 1.5))) + 
  theme(legend.key.width = unit(3.5, "line"))

inkjet_pfi_best_for_plot_yellow = inkjet_pfi_best[[3]]
inkjet_pfi_best_for_plot_yellow$jfpca$f_pca <-
  inkjet_pfi_best_for_plot_yellow$jfpca$f_pca

inkjet_yellow_best_pc_diffs <-
  plot_inkjet_pc_diffs(inkjet_pfi_best_for_plot_yellow) +
  labs(title = "Yellow (best differences)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 2)))

plot_grid(
  inkjet_yellow_best_pc_dirs + theme(legend.position = "none"),
  inkjet_yellow_best_pc_diffs + theme(legend.position = "none"),
  get_plot_component(
    inkjet_yellow_best_pc_dirs,
    'guide-box-bottom', 
    return_all = TRUE
  ),
  ncol = 1,
  rel_heights = c(0.3, 0.3, 0.05)
)
```

```{r figS21, fig.height = 15, fig.width = 25, out.width = "6in", fig.cap = "The principal direction and principal differences plots for the five most important jfPCs from the worst model for yellow inkjet signatures."}

inkjet_yellow_worst_pc_dirs <-
  plot_inkjet_pcs(inkjet_pfi_worst[[3]]) + 
  labs(title = "Yellow (worst directions)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 1.5))) + 
  theme(legend.key.width = unit(3.5, "line"))

inkjet_pfi_worst_for_plot_yellow = inkjet_pfi_worst[[3]]
inkjet_pfi_worst_for_plot_yellow$jfpca$f_pca <-
  inkjet_pfi_worst_for_plot_yellow$jfpca$f_pca

inkjet_yellow_worst_pc_diffs <-
  plot_inkjet_pc_diffs(inkjet_pfi_worst_for_plot_yellow) +
  labs(title = "Yellow (worst differences)") + 
  guides(color = guide_legend(override.aes = list(linewidth = 2)))

plot_grid(
  inkjet_yellow_worst_pc_dirs + theme(legend.position = "none"),
  inkjet_yellow_worst_pc_diffs + theme(legend.position = "none"),
  get_plot_component(
    inkjet_yellow_worst_pc_dirs,
    'guide-box-bottom', 
    return_all = TRUE
  ),
  ncol = 1,
  rel_heights = c(0.3, 0.3, 0.05)
)
```

## Improving Models Using PFI {#inkjet-pfi-remove}

```{r figS22, fig.height = 8, fig.width = 20, out.width = "6.5in", fig.cap = "Cross-validation average accuracies. Black lines represent CV accuracies from models with PCs selected via PFI."}

# Specify the figure font size
fS22_fs = 24

# Plot results in comparison to other accuracies (500 random forest trees):
ggplot() +
  geom_hline(
    data = inkjet_paper_res_best,
    mapping = aes(yintercept = acc_ave),
    linewidth = 1,
    linetype = "dashed"
  ) +
  geom_line(
    data = inkjet_cv_acc_summary |> filter(ntrees == 500),
    aes(
      x = s,
      y = acc_ave,
      group = factor(pcs, levels = seq(10, 100, 10)),
      color = factor(pcs, levels = seq(10, 100, 10))
    ),
    linewidth = 0.75
  ) +
  geom_line(
    data = inkjet_cv_acc_pfi_informed_summary,
    mapping = aes(x = s, y = acc_ave),
    color = "black",
    linewidth = 1.5
  ) +
  facet_grid(. ~ color) +
  scale_color_manual(
    values = wes_palette(
      name = "Zissou1",
      n = 10,
      type = "continuous"
    )
  ) +
  theme_bw(base_size = fS22_fs, base_family = ff) +
  theme(
    strip.background = element_blank(),
    axis.title = element_text(size = fS22_fs),
    axis.text = element_text(size = fS22_fs),
    title = element_text(size = fS22_fs),
    strip.text = element_text(size = fS22_fs),
    legend.text = element_text(size = fS22_fs),
    legend.title = element_text(size = fS22_fs)
  ) +
  guides(color = guide_legend(override.aes = list(linewidth = 1.5))) +
  labs(
    x = "Smoothing Parameter",
    y = "Accuracy",
    color = "Number of PCs"
  )

```
