install.packages(c("ggplot2", "dplyr", "tidyr", "ggrepel", "patchwork", "ggpattern"))

library(ggplot2)
library(dplyr)
library(tidyr)
library(ggrepel)
library(patchwork)
library(ggpattern)
library(grid)

set.seed(123)

T <- 96
t <- 1:T
q <- 115

panel_top <- "Observed and predicted load \n (in kW)"
panel_bottom <- "Residual load after intervention \n (in kW)"

L <- 82 +
  5 * sin(2 * pi * t / 96) +
  115 * exp(-((t - 18)^2) / 120) +
  88 * exp(-((t - 72)^2) / 95) +
  rnorm(T, 0, 1.8)

L_hat <- 78 +
  5 * sin(2 * pi * (t + 2) / 96) +
  59 * exp(-((t - 18)^2) / 120) +
  56 * exp(-((t - 45)^2) / 125) +
  14 * exp(-((t - 72)^2) / 120)

df <- data.frame(t, L, L_hat) %>%
  mutate(
    e_actual = pmax(0, L - q),
    e_pred   = pmax(0, L_hat - q),
    e_eff    = ifelse(L_hat > q, e_pred, 0),
    e_res    = pmax(0, L - e_eff - q),
    L_res    = ifelse(L_hat > q, L - pmax(0, L_hat - q), L)
  )

y_min <- 70
y_max <- q + 110

K_baseline <- max(df$L)
K_new      <- max(df$L_res)
Delta_K    <- K_baseline - K_new

baseline_capacity <- K_baseline
residual_capacity <- max(df$L_res)

C_grid <- 250
p_max <- (C_grid * Delta_K) / sum(df$e_pred)

plot_df <- df %>%
  select(t, Observed_load = L, Predicted_load = L_hat) %>%
  pivot_longer(-t, names_to = "series", values_to = "value")

get_exact_phase_bounds <- function(t, y, q) {
  active <- y > q
  r <- rle(active)
  ends <- cumsum(r$lengths)
  starts <- ends - r$lengths + 1
  idx <- which(r$values)
  
  bounds <- c()
  
  for (j in idx) {
    s <- starts[j]
    e <- ends[j]
    
    x_start <- if (s > 1) {
      t[s - 1] + (q - y[s - 1]) * (t[s] - t[s - 1]) / (y[s] - y[s - 1])
    } else {
      t[s]
    }
    
    x_end <- if (e < length(t)) {
      t[e] + (q - y[e]) * (t[e + 1] - t[e]) / (y[e + 1] - y[e])
    } else {
      t[e]
    }
    
    bounds <- c(bounds, x_start, x_end)
  }
  
  bounds
}

predicted_bounds <- get_exact_phase_bounds(df$t, df$L_hat, q)
observed_bounds  <- get_exact_phase_bounds(df$t, df$L, q)

observed_intervals <- data.frame(
  start = observed_bounds[seq(1, length(observed_bounds), by = 2)],
  end   = observed_bounds[seq(2, length(observed_bounds), by = 2)]
) %>%
  mutate(
    mid = (start + end) / 2,
    pred_at_mid = approx(df$t, df$L_hat, xout = mid)$y
  )

fn_bounds <- observed_intervals %>%
  filter(pred_at_mid <= q) %>%
  select(start, end) %>%
  as.matrix() %>%
  as.vector()

predicted_intervals <- matrix(predicted_bounds, ncol = 2, byrow = TRUE)
fn_intervals_plot <- matrix(fn_bounds, ncol = 2, byrow = TRUE)

tp_start <- predicted_intervals[1, 1]
tp_end   <- predicted_intervals[1, 2]
fp_start <- predicted_intervals[2, 1]
fp_end   <- predicted_intervals[2, 2]

tp_mid <- if (nrow(predicted_intervals) >= 1) mean(predicted_intervals[1, ]) else NA
fp_mid <- if (nrow(predicted_intervals) >= 2) mean(predicted_intervals[2, ]) else NA
fn_mid <- if (nrow(fn_intervals_plot) >= 1) mean(fn_intervals_plot[1, ]) else NA

area_df <- df %>%
  mutate(
    requested_top = ifelse(
      (t >= tp_start & t <= tp_end) |
        (t >= fp_start & t <= fp_end),
      L_hat,
      q
    ),
    residual_top = q
  )

area_df_bottom <- df %>%
  mutate(
    panel = panel_bottom,
    residual_top_bottom = ifelse(L_res > q, L_res, q)
  )

shift_df <- df %>%
  mutate(
    panel = panel_bottom,
    shift_top = L,
    shift_bottom = L_res,
    is_fp = (t >= fp_start & t <= fp_end)
  )

panel_df <- bind_rows(
  plot_df %>%
    mutate(panel = panel_top),
  df %>%
    transmute(
      t,
      series = "Residual_load",
      value = L_res,
      panel = panel_bottom
    )
)

panel_df$panel <- factor(
  panel_df$panel,
  levels = c(
    panel_top,
    panel_bottom
  )
)

vline_df <- data.frame(
  x = c(predicted_bounds, fn_bounds),
  type = c(
    rep("Predicted positive", length(predicted_bounds)),
    rep("False negative", length(fn_bounds))
  )
)

ann_bold <- data.frame(
  panel = panel_top,
  x = c(tp_mid - 0.5, fp_mid, fn_mid),
  y = c(q + 59, q + 62, q + 76),
  label = c("True positive:", "False positive:", "False negative:")
) %>% na.omit()

ann_text <- data.frame(
  panel = panel_top,
  x = c(tp_mid - 0.5, fp_mid, fn_mid),
  y = c(q + 42, q + 48, q + 66),
  label = c(
    "underprediction;\nresidual load > q\nremains",
    "shifting incentivized;\nno actual exceedance",
    "residual peak remains"
  )
) %>% na.omit()

ggplot() +
  
  geom_ribbon_pattern(
    data = area_df %>% mutate(panel = panel_top),
    aes(
      x = t,
      ymin = q,
      ymax = requested_top,
      fill = "Hatched"
    ),
    pattern = "stripe",
    pattern_angle = 45,
    pattern_density = 0.25,
    pattern_spacing = 0.035,
    pattern_colour = "grey30",
    pattern_fill = "grey30",
    colour = NA
  ) +
  
  geom_ribbon(
    data = area_df_bottom,
    aes(x = t, ymin = q, ymax = residual_top_bottom, fill = "grau"),
    alpha = 0.3
  ) +
  
  geom_ribbon(
    data = shift_df %>% filter(is_fp),
    aes(x = t, ymin = shift_bottom, ymax = shift_top, fill = "grau"),
    alpha = 0.3
  ) +
  
  geom_vline(
    data = vline_df,
    aes(xintercept = x),
    linewidth = 0.4,
    linetype = "dotted"
  ) +
  
  geom_segment(
    data = data.frame(panel = panel_bottom, x = tp_start - 1),
    aes(x = x, xend = x - 5, y = q + 10, yend = q + 25),
    linewidth = 0.4,
    linetype = "11"
  ) +
  geom_label(
    data = data.frame(
      panel = panel_bottom,
      x = tp_start - 17,
      y = q + 38,
      label = "missed exceedance:\nno shifting intervention"
    ),
    aes(x = x, y = y, label = label),
    hjust = 0,
    size = 3,
    lineheight = 0.85,
    fill = "white",
    linewidth = 0,
    label.padding = unit(0.15, "lines"),
    inherit.aes = FALSE
  ) +
  
  geom_segment(
    data = data.frame(panel = panel_bottom, x = tp_start + 10),
    aes(x = x, xend = x + 10, y = q + 15, yend = q + 80),
    linewidth = 0.4,
    linetype = "11"
  ) +
  geom_label(
    data = data.frame(
      panel = panel_bottom,
      x = tp_start + 17,
      y = q + 92,
      label = "underpredicted exceedance:\ninsufficient shifting intervention"
    ),
    aes(x = x, y = y, label = label),
    hjust = 0,
    size = 3,
    lineheight = 0.85,
    fill = "white",
    linewidth = 0,
    label.padding = unit(0.15, "lines"),
    inherit.aes = FALSE
  ) +
  
  geom_segment(
    data = data.frame(panel = panel_bottom, x = tp_start + 19),
    aes(x = x, xend = x + 5, y = q + 10, yend = q + 25),
    linewidth = 0.4,
    linetype = "11"
  ) +
  geom_label(
    data = data.frame(
      panel = panel_bottom,
      x = tp_start + 20,
      y = q + 38,
      label = "missed exceedance:\nno shifting intervention"
    ),
    aes(x = x, y = y, label = label),
    hjust = 0,
    size = 3,
    lineheight = 0.85,
    fill = "white",
    linewidth = 0,
    label.padding = unit(0.15, "lines"),
    inherit.aes = FALSE
  ) +
  
  geom_segment(
    data = data.frame(panel = panel_bottom, x = tp_start - 1),
    aes(x = x + 36, xend = x + 36, y = q - 40, yend = q - 10),
    linewidth = 0.4,
    linetype = "11"
  ) +
  geom_label(
    data = data.frame(
      panel = panel_bottom,
      x = tp_start + 28,
      y = q - 10,
      label = "unneeded shifting:\neconomic inefficiency"
    ),
    aes(x = x, y = y, label = label),
    hjust = 0,
    size = 3,
    lineheight = 0.85,
    fill = "white",
    linewidth = 0,
    label.padding = unit(0.15, "lines"),
    inherit.aes = FALSE
  ) +
  
  geom_segment(
    data = data.frame(panel = panel_bottom, x = tp_start + 62),
    aes(x = x, xend = x - 10, y = q + 10, yend = q + 25),
    linewidth = 0.4,
    linetype = "11"
  ) +
  geom_label(
    data = data.frame(
      panel = panel_bottom,
      x = tp_start + 40,
      y = q + 38,
      label = "missed exceedance:\nno shifting intervention"
    ),
    aes(x = x, y = y, label = label),
    hjust = 0,
    size = 3,
    lineheight = 0.85,
    fill = "white",
    linewidth = 0,
    label.padding = unit(0.15, "lines"),
    inherit.aes = FALSE
  ) +
  
  geom_line(
    data = panel_df,
    aes(x = t, y = value, linetype = series),
    linewidth = 0.65,
    alpha = 0.9
  ) +
  
  geom_hline(
    yintercept = q,
    linewidth = 0.7,
    linetype = "11"
  ) +
  
  geom_text(
    data = data.frame(
      panel = c(panel_top, panel_bottom),
      x = max(t),
      y = q + 5,
      label = "Threshold q"
    ),
    aes(x = x, y = y, label = label),
    hjust = 1,
    vjust = -0.2,
    size = 3,
    fontface = "italic",
    inherit.aes = FALSE
  ) +
  
  geom_hline(
    data = data.frame(panel = panel_top),
    aes(yintercept = baseline_capacity),
    linewidth = 0.7,
    linetype = "dotdash",
    color = "grey30"
  ) + 
  geom_text(
    data = data.frame(
      panel = panel_top,
      x = max(t),
      y = baseline_capacity + 5,
      label = paste0(
        "Pre-intervation capacity requirement = ",
        round(baseline_capacity, 0), " kW"
      )
    ),
    aes(x = x, y = y, label = label),
    hjust = 1,
    vjust = -0.2,
    size = 3,
    fontface = "italic",
    inherit.aes = FALSE
  ) +
  
  geom_hline(
    data = data.frame(panel = panel_bottom),
    aes(yintercept = residual_capacity),
    linewidth = 0.7,
    linetype = "dotdash",
    alpha = 0.9
  ) +
  geom_text(
    data = data.frame(
      panel = panel_bottom,
      x = max(t),
      y = residual_capacity + 5,
      label = paste0(
        "Post-intervention capacity requirement = ",
        round(residual_capacity, 0), " kW"
      )
    ),
    aes(x = x, y = y, label = label),
    hjust = 1,
    vjust = -0.2,
    size = 3,
    fontface = "italic",
    inherit.aes = FALSE
  ) +
  
  geom_text(
    data = ann_bold,
    aes(x = x, y = y, label = label),
    fontface = "bold",
    size = 3.0,
    inherit.aes = FALSE
  ) +
  
  geom_text(
    data = ann_text,
    aes(x = x, y = y, label = label),
    size = 3.0,
    lineheight = 0.8,
    inherit.aes = FALSE
  ) +
  
  facet_grid(panel ~ ., scales = "fixed", switch = "y") +
  
  scale_linetype_manual(
    name = NULL,
    values = c(
      "Observed_load" = "solid",
      "Predicted_load" = "dashed",
      "Residual_load" = "solid"
    ),
    breaks = c("Observed_load", "Predicted_load"),
    labels = c(
      "Observed load - pre-intervention (upper panel)\nand post-intervention (lower panel)",
      "Predicted load"
    )
  ) +
  
  scale_fill_manual(
    name = NULL,
    values = c(
      "Hatched" = "white",
      "grau" = "grey50"
    ),
    breaks = c(
      "Hatched",
      "grau"
    ),
    labels = c(
      "Predicted threshold exceedances",
      "Inefficiencies from prediction errors"
    )
  ) +
  
  labs(
    x = "Time interval",
    y = NULL,
    linetype = NULL,
    caption = ""
  ) +
  
  guides(
    linetype = guide_legend(
      order = 1,
      nrow = 1,
      byrow = TRUE,
      keywidth = unit(0.9, "cm")
    ),
    
    fill = guide_legend(
      order = 2,
      nrow = 1,
      byrow = TRUE,
      override.aes = list(
        pattern_spacing = 0.02,
        pattern_density = 0.3
      ),
      keywidth = unit(0.6, "cm"),
      keyheight = unit(0.45, "cm")
    )
  ) +
  
  coord_cartesian(xlim = c(-3, 96), ylim = c(y_min, y_max), clip = "off") +
  
  theme_classic(base_size = 11) +
  theme(
    legend.position = "bottom",
    legend.box = "vertical",
    legend.direction = "horizontal",
    legend.box.just = "center",
    legend.spacing.x = unit(0.5, "cm"),
    legend.spacing.y = unit(0.15, "cm"),
    legend.key.width = unit(0.8, "cm"),
    legend.text = element_text(size = 8),
    legend.margin = margin(t = 2, b = 2),
    
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.y.left = element_text(
      angle = 90,
      face = "bold",
      margin = margin(r = 8)
    ),
    strip.text.y.right = element_blank(),
    panel.spacing.y = unit(0.9, "lines"),
    plot.title = element_text(face = "bold", size = 12, margin = margin(b = 16)),
    plot.subtitle = element_text(size = 10, margin = margin(b = 16)),
    plot.caption = element_text(size = 8, hjust = 0),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 9),
    plot.margin = margin(10, 10, 55, 10)
  )