pacman::p_load(tidyverse, ggthemes, gghighlight, ggrepel, patchwork, knitr, scales, colorspaces, ggh4x, RColorBrewer, ggdist)

df <- read_csv("data/respopagesex2024.csv")

# === STEP 1: Find the Maximum Age ===
max_age <- max(df$Age_Num, na.rm = TRUE)

# === STEP 2: Prepare the Data ===
plot_df <- df %>%
  mutate(
    AgeGroup = cut(
      Age_Num,
      breaks = if (max_age > 90) {
        c(seq(0, 85, by = 5), 90, max_age + 1)
      } else {
        c(seq(0, 85, by = 5), max_age + 1)
      },
      right = FALSE,
      labels = if (max_age > 90) {
        c(paste0(seq(0, 80, by = 5), "-", seq(4, 84, by = 5)), "90+")
      } else {
        c(paste0(seq(0, 80, by = 5), "-", seq(4, 84, by = 5)), "90+")
      }
    )
  ) %>%
  group_by(AgeGroup, Sex) %>%
  summarise(Population = sum(Population), .groups = "drop") %>%
  mutate(
    Population = ifelse(Sex == "Males", -Population, Population),
    PopulationAbs = abs(Population)
  )

# === STEP 3: Plot ===
ggplot(plot_df, aes(x = AgeGroup, y = Population, fill = Sex)) +
  geom_bar(stat = "identity", width = 0.90, color = "black") +
  geom_text(
    aes(label = format(PopulationAbs, big.mark = ",")),
    size = 4.5,
    color = "black",
    hjust = ifelse(plot_df$Sex == "Males", -0.05, 1.05),
    show.legend = FALSE
  ) +
  coord_flip() +
  labs(
    title = "Singapore's Population Pyramid (2024) by 5-Year Age Groups",
    x = "Age Group",
    y = "Population Size"
  ) +
  scale_y_continuous(
    labels = function(x) format(abs(x), big.mark = ",", scientific = FALSE),
    limits = c(-200000, 200000),
    breaks = seq(-200000, 200000, by = 50000),
    minor_breaks = seq(-200000, 200000, by = 25000)
  ) +
  scale_fill_manual(
    values = c("Males" = "#87CEFA", "Females" = "#FFB6C1"),
    breaks = c("Males", "Females")
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5),
    legend.position = "bottom",
    axis.text.y = element_text(hjust = 0),
    panel.grid.minor.x = element_line(color = "grey85", size = 0.3)
  )

# === STEP 1: Prepare Overall Population-weighted Data ===
overall_df <- df %>%
  filter(!is.na(Age_Num)) %>%
  group_by(Region, Age_Num) %>%
  summarise(Population = sum(Population), .groups = "drop") %>%
  uncount(weights = Population)

# === STEP 2: Compute Peak Age by Region ===
overall_peaks <- overall_df %>%
  group_by(Region) %>%
  summarise(
    d = list(density(Age_Num)),
    .groups = "drop"
  ) %>%
  rowwise() %>%
  mutate(
    peak_x = d$x[which.max(d$y)],
    label = paste0(Region, " - Peak Age: ", round(peak_x))
  ) %>%
  ungroup()

# === STEP 3: Combine Labels for Annotation ===
overall_annotation <- paste(overall_peaks$label, collapse = "\n")

# === STEP 4: Create Overall Density Plot ===
overall_plot <- ggplot(overall_df, aes(x = Age_Num, fill = Region)) +
  geom_density(alpha = 0.6, color = "black", linewidth = 0.7) +
  annotate("text", x = 5, y = Inf, vjust = 1.2,
           label = overall_annotation, hjust = 0, size = 3.5, fontface = "italic") +
  scale_x_continuous(breaks = seq(0, 90, 10)) +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "Age Distribution (2024): Overall",
    subtitle = "Density by Region with Y-axis scaled per chart",
    x = "Age", y = "Relative Density", fill = "Region"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
    plot.subtitle = element_text(hjust = 0.5, size = 16),
    axis.title.x = element_text(size = 16),
    axis.title.y = element_text(size = 16),
    axis.text.x  = element_text(size = 14),
    axis.text.y  = element_text(size = 14),
    panel.border = element_rect(color = "grey50", fill = NA),
    legend.position = "right"
  )

# === STEP 5: Render ===
overall_plot

# === STEP 1: Prepare Population-weighted Data ===
violin_df <- df %>%
  filter(!is.na(Region), !is.na(Age_Num)) %>%
  uncount(weights = Population) %>%
  mutate(Planning_Area = factor(Planning_Area))

# === STEP 2: Define Region Base Colors ===
region_base_colors <- c(
  "Central"     = "#66C2A5",  # teal
  "East"        = "#FC8D62",  # orange
  "North"       = "#8DA0CB",  # blue
  "North-East"  = "#E78AC3",  # pink
  "West"        = "#A6D854"   # green
)

# === STEP 3: Compute Median Age Helper ===
get_median_labels <- function(data, group_col) {
  data %>%
    group_by(across(all_of(group_col))) %>%
    summarise(median_age = median(Age_Num), .groups = "drop") %>%
    mutate(label = round(median_age, 1))
}

# === STEP 4: Region Function – Violin + Boxplot by Planning Area ===
make_violin_region <- function(region_name, base_color) {
  region_df <- violin_df %>% filter(Region == region_name)

  # Compute medians and sort Planning Areas in descending order
  pa_medians <- get_median_labels(region_df, "Planning_Area") %>%
    arrange(desc(median_age))

  # Reorder factors by descending median age
  region_df$Planning_Area <- factor(region_df$Planning_Area,
                                    levels = pa_medians$Planning_Area)
  pa_medians$Planning_Area <- factor(pa_medians$Planning_Area,
                                     levels = pa_medians$Planning_Area)

  # Reverse gradient (darkest color left)
  n_colors <- length(levels(region_df$Planning_Area))
  color_shades <- colorRampPalette(c(base_color, "white"))(n_colors)
  names(color_shades) <- levels(region_df$Planning_Area)

  ggplot(region_df, aes(x = Planning_Area, y = Age_Num, fill = Planning_Area)) +
    geom_violin(trim = FALSE, scale = "area", color = "black", alpha = 0.5) +
    geom_boxplot(width = 0.15, outlier.shape = NA, alpha = 0.8, color = "black") +
    ggrepel::geom_text_repel(
      data = pa_medians,
      aes(x = Planning_Area, y = median_age, label = label),
      inherit.aes = FALSE,
      size = 3.5, fontface = "italic",
      max.overlaps = Inf,
      min.segment.length = 0,
      nudge_y = 0.5
    ) +
    scale_fill_manual(values = color_shades) +
    labs(
      title = paste("Age Distribution by Planning Area: ", region_name),
      x = NULL,
      y ="Age"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
      axis.title.y = element_text(size = 16),
      axis.text.x  = element_text(angle = 90, hjust = 0.5, vjust = 0.5, size = 14),
      axis.text.y  = element_text(size = 14),
      panel.border = element_rect(color = "grey50", fill = NA),
      legend.position = "none"
    )
}

# === STEP 5: Region Violin Plots ===
central_violin     <- make_violin_region("Central", region_base_colors["Central"])
east_violin        <- make_violin_region("East", region_base_colors["East"])
north_violin       <- make_violin_region("North", region_base_colors["North"])
north_east_violin  <- make_violin_region("North-East", region_base_colors["North-East"])
west_violin        <- make_violin_region("West", region_base_colors["West"])

# === STEP 6: Assemble Vertical Layout ===
violin_boxplot_grid <- 
  central_violin /
  east_violin /
  north_violin /
  north_east_violin /
  west_violin +
  plot_layout(ncol = 1) +
  plot_annotation(
    title = "Violin + Boxplot of Age Distribution (2024)",
    subtitle = "Regional breakdown by Planning Areas sorted by Median Age (Descending)\nColor fades left to right by age rank",
    theme = theme(
      plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
      plot.subtitle = element_text(size = 11, hjust = 0.5)
    )
  )

# === STEP 7: Render ===
violin_boxplot_grid