## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "##",
  fig.width = 6,
  fig.height = 4,
  out.width = "90%"
)

library(tidyverse)
library(viridisLite)

theme_set(theme_minimal() + theme(legend.position = "bottom"))

options(
  ggplot2.continuous.colour = "viridis",
  ggplot2.continuous.fill = "viridis",
  ggplot2.discrete.colour = "viridis_d",
  ggplot2.discrete.fill = "viridis_d"
)

library("tidyfun")
data(chf_df, package = "tidyfun")
data(dti_df, package = "tidyfun")

pal_5 <- viridis(7)[-(1:2)]
set.seed(1221)

## ----plot_chf-----------------------------------------------------------------
dti_df[1:10,] |>
  tf_ggplot(aes(tf = cca)) + geom_line(alpha = .3)

## -----------------------------------------------------------------------------
dti_df[1:3,] |>
  tf_ggplot(aes(tf = rcst)) + geom_line(alpha = .3) + geom_point(alpha= .3)

## -----------------------------------------------------------------------------
chf_df |>
  filter(id %in% 1:5) |>
  tf_ggplot(
    aes(tf = tf_smooth(activity, f = .05), # smoothed inputs for clearer viz
        color = gender)) +
  geom_line(alpha = 0.3)

## -----------------------------------------------------------------------------
chf_df |>
  filter((id %in% 1:10) & (day %in% c("Mon", "Sun"))) |>
  tf_ggplot(aes(tf = tf_smooth(activity, f = .05), color = gender)) +
  geom_line(alpha = 0.3, lwd = 1) +
  facet_grid(~day)

## ----dti-fig1-----------------------------------------------------------------
dti_df |>
  tf_ggplot(aes(tf = cca, col = case, alpha = 0.2 + 0.4 * (case == "control"))) +
  geom_line() + facet_wrap(~sex) +
  scale_alpha(guide = "none", range = c(0.2, 0.4))

## -----------------------------------------------------------------------------
chf_df |>
  group_by(gender, day) |>
  summarize(mean_act = mean(activity),
            .groups = "drop_last") |>
  mutate(smooth_mean = tfb(mean_act, verbose = FALSE)) |>
  filter(day %in% c("Mon", "Sun")) |>
  tf_ggplot(aes(color = gender)) +
  geom_line(aes(tf = smooth_mean), linewidth = 1.25) +
  geom_line(aes(tf = mean_act), alpha = 0.1) +
  geom_point(aes(tf = mean_act), alpha = 0.1, size = .1) +
  facet_grid(day~.)

## -----------------------------------------------------------------------------
dti_df |>
  group_by(sex, case) |>
  summarize(
    mean_cca = mean(tfb(cca, verbose = FALSE)), #pointwise mean function
    sd_cca = sd(tfb(cca, verbose = FALSE)), # pointwise sd function
    .groups = "drop_last"
  ) |>
  group_by(sex, case) |>
  mutate(
    upper_cca = mean_cca + 2 * sd_cca,
    lower_cca = mean_cca - 2 * sd_cca
  ) |>
  tf_ggplot() +
  geom_line(aes(tf = mean_cca, color = sex)) +
  geom_ribbon(aes(tf_ymin = lower_cca, tf_ymax = upper_cca, fill = sex), alpha = 0.3) +
  facet_grid(sex ~ case)

## -----------------------------------------------------------------------------
dti_df |> 
  tf_ggplot(aes(tf = cca, fill = case)) +
  geom_fboxplot(alpha = 0.35) +
  facet_grid(~ sex) + labs(title="MBD-based boxplot")

## -----------------------------------------------------------------------------
dti_df |>
  tf_ggplot(aes(tf = cca, colour = case)) +
  geom_fboxplot(depth = "FM", alpha = 0.3) +
  facet_grid(~ sex) + labs(title="FM-based boxplot")

## -----------------------------------------------------------------------------
dti_df |>
  tf_ggplot(aes(tf = cca, colour = case)) +
  geom_fboxplot(depth = "RPD", alpha = 0.3) +
  facet_grid(~ sex) + labs(title="RPD-based boxplot")

## -----------------------------------------------------------------------------
tf_ggplot(dti_df, aes(tf = rcst)) + geom_fboxplot()

## -----------------------------------------------------------------------------
tf_ggplot(dti_df, aes(tf = rcst)) + 
  geom_fboxplot(alpha = .5)
tf_ggplot(dti_df, aes(tf = rcst)) + 
  geom_fboxplot(alpha = .5, central = .2)
tf_ggplot(dti_df, aes(tf = rcst)) + 
  geom_fboxplot(alpha = .5, central = .2, outliers = FALSE)
tf_ggplot(dti_df, aes(tf = rcst)) +
  geom_fboxplot(orientation = "y", alpha = .3)

## -----------------------------------------------------------------------------
chf_df |>
  filter(day %in% c("Mon", "Sun")) |>
  gglasagna(activity)

## ----dti-fig2-----------------------------------------------------------------
dti_df |>
  gglasagna(
    tf = cca,
    order = tf_integrate(cca, definite = TRUE), #order by area under the curve
    arg = seq(0, 1, length.out = 101)
  ) +
  theme(axis.text.y = element_text(size = 6)) +
  facet_wrap(~ case:sex, ncol = 2, scales = "free")

## -----------------------------------------------------------------------------
canada <- data.frame(
  place = fda::CanadianWeather$place,
  region = fda::CanadianWeather$region,
  lat = fda::CanadianWeather$coordinates[, 1],
  lon = -fda::CanadianWeather$coordinates[, 2]
)

canada$temp <- tfd(t(fda::CanadianWeather$dailyAv[, , 1]), arg = 1:365)
canada$precipl10 <- tfd(t(fda::CanadianWeather$dailyAv[, , 3]), arg = 1:365) |>
  tf_smooth()

canada_map <-
  data.frame(maps::map("world", "Canada", plot = FALSE)[c("x", "y")])

## -----------------------------------------------------------------------------
ggplot(canada, aes(x = lon, y = lat)) +
  geom_capellini(aes(tf = precipl10),
    width = 4, height = 5, colour = "blue",
    line.linetype = 1
  ) +
  geom_capellini(aes(tf = temp),
    width = 4, height = 5, colour = "red",
    line.linetype = 1
  ) +
  geom_path(data = canada_map, aes(x = x, y = y), alpha = 0.1) +
  coord_quickmap()

## ----warning=FALSE------------------------------------------------------------
cca_fpc_tbl <- tibble(
  cca = dti_df$cca[1:30],
  cca_fpc = tfb_fpc(cca, pve = .8), 
  fpc_1 = map(coef(cca_fpc), 2) |> unlist(), # 1st PC loading
  fpc_2 = map(coef(cca_fpc), 3) |> unlist() # 2nd PC loading
) 
# rescale FPCs by sqrt of eigenvalues for visualization
cca_fpcs_1_2 <- 
  tf_basis(cca_fpc_tbl$cca_fpc, as_tfd = TRUE)[2:3] * 
    sqrt(attr(cca_fpc_tbl$cca_fpc, "score_variance")[1:2]) 
# scaled eigenfunctions look like this:
tibble(
   eigenfunction = cca_fpcs_1_2,
   FPC = factor(1:2)
) |> tf_ggplot() + 
  geom_line(aes(tf = eigenfunction, col = FPC)) + 
  geom_hline(yintercept = 0)

## ----warning=FALSE------------------------------------------------------------
ggplot(cca_fpc_tbl[1:40,], aes(x =  fpc_1, y = fpc_2)) +
  geom_point(size = .5, col = viridis(3)[2]) +
  geom_capellini(aes(tf =cca_fpc),width = .01, height = .01, line.linetype = 1) +
  labs(x = "FPC1 score", y = "FPC2 score")

## -----------------------------------------------------------------------------
cca <- dti_df$cca |>
  tfd(arg = seq(0, 1, length.out = 93), interpolate = TRUE)

layout(t(1:2))

plot(cca, type = "spaghetti")
lines(c(median(cca), mean = mean(cca)), col = viridis(3)[c(1, 3)])

plot(cca, type = "lasagna", col = viridis(50))

## ----ex-fig2------------------------------------------------------------------
cca_five <- cca[1:5]

cca_five |> plot(xlim = c(-0.15, 1), col = pal_5, lwd = 2)

text(
  x = -0.1, y = cca_five[, 0.07], labels = names(cca_five), col = pal_5, cex = 1
)

median(cca_five) |> lines(col = pal_5[3], lwd = 4)

## -----------------------------------------------------------------------------
pinch_reg <- tf::pinch |> tfb() |> #smooth before registration for better results
  tf_register() 
pinch_reg
summary(pinch_reg)
plot(pinch_reg)

## -----------------------------------------------------------------------------
layout(t(1:3))
plot(tf::pinch[1:5], col = pal_5, lwd = 2, points = FALSE)

plot(tf_inv_warps(pinch_reg)[1:5], col = pal_5, lwd = 2, points = FALSE)
abline(c(0, 1), col = "grey", lty = 2)

plot(tf_aligned(pinch_reg)[1:5], col = pal_5, lwd = 2)
lines(tf_template(pinch_reg), col = "black", lwd = 3, lty= 3)