In FlukeAndFeather/cetaceanbcg: Cetacean Heart Beats Detected by Ballistocardiogram

knitr::opts_chunk$set(
  echo = FALSE,
  cache = TRUE,
  message = FALSE,
  fig.height = 5,
  fig.width = 5,
  dpi = 300
)
library(cetaceanbcg)
library(tidyverse)

bw_elg <- bw180905_53_elg %>% 
  mutate(
    depth_min = map2_dbl(start, stop, function(t1, t2) {
      bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% min()
    }),
    depth_max = map2_dbl(start, stop, function(t1, t2) {
      bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% max()
    }),
    depth_med = map2_dbl(start, stop, function(t1, t2) {
      bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% 
        first() %>% 
        median()
    }),
    ygyr_min = map2_dbl(start, stop, function(t1, t2) {
      bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), ]$gw[,2] %>% min()
    }),
    ygyr_max = map2_dbl(start, stop, function(t1, t2) {
      bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), ]$gw[,2] %>% max()
    }),
    ygyr_med = map2_dbl(start, stop, function(t1, t2) {
      bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), ]$gw[,2] %>% 
        first() %>% 
        median()
    }),
  )

bw_10hz <- bw180905_53_10hz %>% 
  mutate(rid_left = approx(bw_elg$start, 
                           bw_elg$region_id, 
                           dt, 
                           "constant")$y,
         rid_right = approx(bw_elg$stop, 
                            bw_elg$region_id, 
                            dt, 
                            "constant", 
                            yleft = 0)$y + 1,
         region_id = ifelse(rid_left == rid_right, rid_left, NA),
         dive_id = split_dives(dt, depth,
                          surface = 2, min_depth = 10, min_dur = 5 * 60))

bw_400hz <- bw180905_53_400hz %>%
  mutate(
    across(surge:heave, 
           filter_acc, fs = 400, upper = 10.0, 
           .names = "{.col}_filt"),
    jerk = jerk(cbind(surge_filt, sway_filt, heave_filt), 
                fs = 400, p = 4, n = 2 * 400 + 1),
    jerk_se = shannon_entropy(jerk),
    jerk_smooth = tma(jerk_se, 2 * 400),
    # Annotate regions
    rid_left = approx(bw_elg$start, 
                      bw_elg$region_id, 
                      dt, 
                      "constant")$y,
    rid_right = approx(bw_elg$stop, 
                       bw_elg$region_id, 
                       dt, 
                       "constant", 
                       yleft = 0)$y + 1,
    region_id = ifelse(rid_left == rid_right, rid_left, NA),
    # Zero-out signal in non-valid regions (i.e. remove movement artifacts)
    jerk_smooth = ifelse(is.na(region_id), 0, jerk_smooth)
  )

fs_hz <- 400
minperiod <- 2
cue <- bw_400hz$jerk_smooth
peaks <- pracma::findpeaks(cue, minpeakdistance = fs_hz * minperiod)
proms <- peak_prominences(cue, peaks[, 2])

dist_max <- sqrt((max(peaks[, 1]) - peaks[, 1])^2 + (max(proms) - proms)^2)
dist_dens <- stats::density(dist_max)
dist_thr <- dist_dens$x[pracma::findpeaks(-dist_dens$y)[, 2]]
is_beat <- peaks[, 2][dist_max <= dist_thr]
beats <- logical(length(cue))
beats[is_beat] <- TRUE

peak_tbl <- tibble(dt = bw_400hz$dt[peaks[, 2]],
                   height = peaks[, 1],
                   prominence = proms,
                   is_beat = dist_max <= dist_thr)

regionsample <- filter(bw_400hz, region_id == 6) %>% 
  mutate(secs = as.numeric(dt - min(dt), unit = "secs"))
peaksample <- filter(peak_tbl, between(dt, 
                                       first(regionsample$dt), 
                                       last(regionsample$dt))) %>% 
  arrange(dt) %>% 
  mutate(peak_lbl = seq(n()),
         secs = as.numeric(dt - min(regionsample$dt)))

max_height <- max(peaks[, 1])
thr_curve <- tibble(
  t = seq(0, 2 * pi, length.out = 1000),
  height = max_height + cos(t) * dist_thr,
  prominence = max_height + sin(t) * dist_thr,
) %>% 
  filter(height <= max_height, prominence <= max_height)

pA <- ggplot(regionsample, aes(secs, jerk_smooth)) +
  geom_line() +
  geom_point(aes(y = height, color = is_beat), data = peaksample) +
  geom_text(aes(y = height, label = peak_lbl), 
            data = peaksample, 
            vjust = -0.3) +
  scale_color_manual(values = c("red", "blue")) +
  labs(x = "Seconds", y = "Y") +
  theme_classic(base_size = 10) +
  theme(legend.position = "FALSE")

pB <- ggplot(peaksample, aes(height, prominence, color = is_beat)) +
  geom_point() +
  ggrepel::geom_text_repel(aes(label = peak_lbl), force = 3) +
  scale_color_manual(values = c("red", "blue")) +
  scale_x_continuous(n.breaks = 4) +
  labs(x = "Peak height",
       y = "Peak prominence") +
  coord_fixed() +
  theme_classic(base_size = 10) +
  theme(legend.position = "FALSE")

pC <- ggplot(peak_tbl, aes(dist_max)) +
  geom_density() +
  geom_vline(xintercept = dist_thr, linetype = 2, color = "purple") +
  scale_x_continuous(n.breaks = 4) +
  labs(x = "Distance from highest peak") +
  theme_classic(base_size = 10)

pD <- ggplot(peak_tbl, aes(height, prominence)) +
  geom_point(aes(color = is_beat), shape = 21) +
  geom_line(data = thr_curve, linetype = 2, color = "purple") +
  scale_color_manual(values = c("red", "blue")) +
  labs(x = "Peak height",
       y = "Peak prominence") +
  coord_fixed() +
  theme_classic(base_size = 10) +
  theme(legend.position = "FALSE")

# Figure S1 requires manual adjustments. Save to .svg, edit in Inkscape, then
# export back to PNG and insert here.
S1 <- cowplot::plot_grid(pA, pB, pC, pD, ncol = 2, labels = "AUTO")
ggsave(here::here("analysis/figures/S1.svg"), S1)
cowplot::ggdraw() +
  cowplot::draw_image(here::here("analysis/figures/S1_edited.svg"))