rm(list=ls())
library(shiny)
library(arrow)
library(dplyr)
library(stringr)
library(colorspace)
library(zoo)

develop <- FALSE
if (develop){
  results_directory <- "results/"
}else{
  results_directory <- "/data/running-calibration"
}

if (!exists("%||%")) {
  `%||%` <- function(x, y) {
    if (is.null(x) || identical(x, "")) y else x
  }
}

# --- Helper ---
adjust_transparency <- function(color, alpha=0.3) {
  adjustcolor(color, alpha.f=alpha)
}

safe_ylim <- function(values, probs = c(0.05, 0.95), fallback = c(-1, 1)) {
  vals <- as.numeric(values)
  vals <- vals[is.finite(vals)]
  if (length(vals) == 0) {
    return(fallback)
  }

  y <- as.numeric(quantile(vals, probs = probs, na.rm = TRUE))
  if (length(y) != 2 || !all(is.finite(y))) {
    return(fallback)
  }

  if (y[1] == y[2]) {
    pad <- ifelse(y[1] == 0, 1, abs(y[1]) * 0.1)
    y <- c(y[1] - pad, y[2] + pad)
  }
  y
}

expand_ylim <- function(ylim, values) {
  vals <- as.numeric(values)
  vals <- vals[is.finite(vals)]
  if (length(vals) == 0) {
    return(ylim)
  }

  y <- range(c(ylim, vals), na.rm = TRUE)
  if (y[1] == y[2]) {
    pad <- ifelse(y[1] == 0, 1, abs(y[1]) * 0.1)
    y <- c(y[1] - pad, y[2] + pad)
  }
  y
}

format_scientific_value <- function(x, digits = 6) {
  vals <- suppressWarnings(as.numeric(x))
  out <- rep(NA_character_, length(vals))
  keep <- is.finite(vals)
  out[keep] <- format(signif(vals[keep], digits), scientific = TRUE, trim = TRUE)
  out
}

format_fixed_value <- function(x, digits = 2) {
  vals <- suppressWarnings(as.numeric(x))
  out <- rep(NA_character_, length(vals))
  keep <- is.finite(vals)
  out[keep] <- formatC(vals[keep], format = "f", digits = digits)
  out
}

draw_calibration_uncertainty_bars <- function(calibration_uncertainty, plot_dates, sigma_level = 1, col = "#D55E00") {
  if (is.null(calibration_uncertainty) || nrow(calibration_uncertainty) == 0) {
    return(FALSE)
  }

  plot_dates <- as.Date(plot_dates)
  plot_dates <- plot_dates[!is.na(plot_dates)]
  if (length(plot_dates) == 0) {
    return(FALSE)
  }

  plot_min <- min(plot_dates)
  plot_max <- max(plot_dates)

  df <- calibration_uncertainty %>%
    mutate(
      cal_date = as.Date(cal_date),
      scref_uncert = abs(as.numeric(scref_uncert)) * as.numeric(sigma_level)
    ) %>%
    filter(!is.na(cal_date), is.finite(scref_uncert)) %>%
    arrange(cal_date)

  if (nrow(df) == 0) {
    return(FALSE)
  }

  valid_to <- c(df$cal_date[-1], plot_max)
  x0 <- as.Date(pmax(as.numeric(df$cal_date), as.numeric(plot_min)), origin = "1970-01-01")
  x1 <- as.Date(pmin(as.numeric(valid_to), as.numeric(plot_max)), origin = "1970-01-01")
  keep <- !is.na(x0) & !is.na(x1) & x1 >= x0
  if (!any(keep)) {
    return(FALSE)
  }

  segments(x0[keep], df$scref_uncert[keep], x1[keep], df$scref_uncert[keep],
           col = col, lwd = 3, lend = "butt")
  segments(x0[keep], -df$scref_uncert[keep], x1[keep], -df$scref_uncert[keep],
           col = col, lwd = 3, lend = "butt")
  TRUE
}

draw_synt_ref_periods <- function(calibration_uncertainty, plot_dates, col = "#D55E00") {
  if (is.null(calibration_uncertainty) || nrow(calibration_uncertainty) == 0) {
    return(FALSE)
  }

  plot_dates <- as.Date(plot_dates)
  plot_dates <- plot_dates[!is.na(plot_dates)]
  if (length(plot_dates) == 0) {
    return(FALSE)
  }

  df <- calibration_uncertainty %>%
    mutate(
      synt_ref_start = as.Date(synt_ref_start, tz = "UTC"),
      synt_ref_end = as.Date(synt_ref_end, tz = "UTC")
    ) %>%
    filter(!is.na(synt_ref_start), !is.na(synt_ref_end), synt_ref_end >= synt_ref_start)

  if (nrow(df) == 0) {
    return(FALSE)
  }

  x_min <- min(as.Date(plot_dates), na.rm = TRUE)
  x_max <- max(as.Date(plot_dates), na.rm = TRUE)
  plot_span_days <- max(as.numeric(x_max - x_min), 1)
  min_visible_days <- min(30, max(5, ceiling(plot_span_days * 0.005)))

  x0 <- df$synt_ref_start
  x1 <- pmax(df$synt_ref_end, df$synt_ref_start + 1)
  actual_width_days <- pmax(as.numeric(x1 - x0), 1)
  expand_days <- pmax(0, min_visible_days - actual_width_days) / 2
  x0 <- x0 - ceiling(expand_days)
  x1 <- x1 + floor(expand_days)
  keep <- is.finite(as.numeric(x0)) & is.finite(as.numeric(x1)) & x1 >= x_min & x0 <= x_max
  if (!any(keep)) {
    return(FALSE)
  }

  usr <- par("usr")
  rect(
    xleft = pmax(x0[keep], x_min),
    ybottom = usr[3],
    xright = pmin(x1[keep], x_max + 1),
    ytop = usr[4],
    col = adjust_transparency(col, alpha = 0.12),
    border = adjust_transparency(col, alpha = 0.45),
    lwd = 1.5
  )
  TRUE
}

plot_result <- function(x, result="intercept", title="", selected_day = NULL,
                        calibration_uncertainty = NULL,
                        show_calibration_uncertainty = FALSE,
                        sigma_level = 1) {
  windows <- length(x) + 1
  colors  <- qualitative_hcl(n = windows, h = c(-174, 123), c = 46, l = 29)

  y_vals <- unlist(lapply(x, function(df) df[[result]]), use.names = FALSE)
  plot_dates <- do.call(c, lapply(x, function(df) as.Date(df[["date"]])))
  calibration_values <- numeric(0)
  draw_uncertainty <- isTRUE(show_calibration_uncertainty) &&
    identical(result, "intercept") &&
    !is.null(calibration_uncertainty) &&
    nrow(calibration_uncertainty) > 0
  if (draw_uncertainty) {
    calibration_values <- c(
      abs(as.numeric(calibration_uncertainty$scref_uncert)) * as.numeric(sigma_level),
      -abs(as.numeric(calibration_uncertainty$scref_uncert)) * as.numeric(sigma_level)
    )
  }

  plot(as.Date(x[[length(x)]][["date"]]), x[[length(x)]][[result]],
       type = "n",
       ylim=expand_ylim(safe_ylim(y_vals), calibration_values),
       ylab=result,
       xlab="",
       main=title)
  grid()

  period_drawn <- FALSE
  if (draw_uncertainty) {
    period_drawn <- draw_synt_ref_periods(calibration_uncertainty, plot_dates)
  }

  uncertainty_drawn <- FALSE
  if (draw_uncertainty) {
    uncertainty_drawn <- draw_calibration_uncertainty_bars(calibration_uncertainty, plot_dates, sigma_level = sigma_level)
  }
  
  for (w in seq_along(x)) {
    points(x[[w]][[result]] ~ as.Date(x[[w]][["date"]]),
           col=adjust_transparency(colors[w], alpha=0.3),
           pch=16)
  }

  if (!is.null(selected_day)) {
    abline(v = as.Date(selected_day), col = "red", lwd = 2, lty = 2)
  }

  abline(h=0, lwd=2, lty=2)
  legend_labels <- names(x)
  legend_cols <- adjust_transparency(colors[seq_along(x)], alpha=0.4)
  legend_pch <- rep(16, length(legend_labels))
  legend_lty <- rep(NA, length(legend_labels))
  legend_lwd <- rep(NA, length(legend_labels))
  legend("topleft", legend_labels,
         col=legend_cols,
         bty="n", pch=legend_pch, lty=legend_lty, lwd=legend_lwd)

  if (uncertainty_drawn || period_drawn) {
    extra_labels <- character(0)
    extra_cols <- character(0)
    extra_lty <- numeric(0)
    extra_lwd <- numeric(0)
    extra_fill <- character(0)

    if (uncertainty_drawn) {
      extra_labels <- c(extra_labels, sprintf("Initial cal. uncertainty (%s-sigma)", as.character(sigma_level)))
      extra_cols <- c(extra_cols, "#D55E00")
      extra_lty <- c(extra_lty, 1)
      extra_lwd <- c(extra_lwd, 3)
      extra_fill <- c(extra_fill, NA_character_)
    }
    if (period_drawn) {
      extra_labels <- c(extra_labels, "Synthetic ref period")
      extra_cols <- c(extra_cols, NA_character_)
      extra_lty <- c(extra_lty, NA_real_)
      extra_lwd <- c(extra_lwd, NA_real_)
      extra_fill <- c(extra_fill, adjust_transparency("#D55E00", alpha = 0.12))
    }

    legend(
      "topright",
      legend = extra_labels,
      col = extra_cols,
      lty = extra_lty,
      lwd = extra_lwd,
      fill = extra_fill,
      border = NA,
      bty = "n"
    )
  }
}

plot_daily_mle_fit <- function(fit_data, selected_day = NULL, title = "Daily MLE Fit") {
  if (is.null(fit_data) || nrow(fit_data) == 0) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No MLE fit diagnostics found for this selection")
    return(invisible(NULL))
  }

  if ("date" %in% names(fit_data)) {
    fit_data$date <- as.Date(fit_data$date)
    if (!is.null(selected_day)) {
      fit_data <- subset(fit_data, date == as.Date(selected_day))
    } else {
      fit_data <- subset(fit_data, date == max(date, na.rm = TRUE))
    }
  } else if (!is.null(selected_day)) {
    fit_data$date <- rep(as.Date(selected_day), nrow(fit_data))
  }

  raw <- subset(fit_data, point_type == "raw")
  binned <- subset(fit_data, point_type == "binned")

  if (nrow(raw) == 0 && nrow(binned) == 0) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No points available for this day")
    return(invisible(NULL))
  }

  x_vals <- c(raw$amf, binned$amf)
  y_vals <- c(raw$sc, binned$sc, binned$fitted_sc)

  x_rng <- range(x_vals, na.rm = TRUE)
  y_rng <- range(y_vals, na.rm = TRUE)

  if (!all(is.finite(x_rng))) x_rng <- c(-1, 1)
  if (!all(is.finite(y_rng))) y_rng <- c(-1, 1)
  if (x_rng[1] == x_rng[2]) x_rng <- x_rng + c(-1, 1) * 0.1
  if (y_rng[1] == y_rng[2]) y_rng <- y_rng + c(-1, 1) * 0.1

  plot(raw$amf, raw$sc,
       pch = 16,
       cex = 0.8,
       col = adjust_transparency("black", 0.2),
       xlim = x_rng,
       ylim = y_rng,
       xlab = "AMF",
       ylab = "Slant Column [mol/m^2]",
       main = title)
  grid()

  fit_label <- "Linear fit"
  if (nrow(binned) > 0) {
    points(binned$amf, binned$sc, pch = 16, col = "blue")
    if (any(is.finite(binned$fitted_sc))) {
      binned <- binned[order(binned$amf), ]
      lines(binned$amf, binned$fitted_sc, col = "red", lwd = 2)

      fit_points <- binned[
        is.finite(binned$amf) & is.finite(binned$fitted_sc),
        ,
        drop = FALSE
      ]
      if (nrow(fit_points) >= 2 && length(unique(fit_points$amf)) >= 2) {
        fit_intercept <- coef(lm(fitted_sc ~ amf, data = fit_points))[1]
        if (is.finite(fit_intercept)) {
          fit_label <- sprintf("Linear fit (intercept = %.2e mol/m^2)", fit_intercept)
        }
      }
    } else {
      fit_points <- binned[
        is.finite(binned$amf) & is.finite(binned$sc),
        ,
        drop = FALSE
      ]
      if (nrow(fit_points) >= 2 && length(unique(fit_points$amf)) >= 2) {
        fit_intercept <- coef(lm(sc ~ amf, data = fit_points))[1]
        if (is.finite(fit_intercept)) {
          fit_label <- sprintf("Linear fit (intercept = %.2e mol/m^2)", fit_intercept)
        }
      }
    }
  }

  legend("topleft",
         legend = c("Raw points", "Binned percentile", fit_label),
         col = c(adjust_transparency("black", 0.4), "blue", "red"),
         pch = c(16, 16, NA),
         lty = c(NA, NA, 1),
         bty = "n")
}

plot_daily_columns <- function(day_data, selected_day = NULL, title = "Daily Column Data") {
  if (is.null(day_data) || nrow(day_data) == 0) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No daily column data found")
    return(invisible(NULL))
  }

  day_data$timestamp <- as.POSIXct(day_data$timestamp, tz = "UTC")
  vc_cols <- grep("^vc", names(day_data), value = TRUE)

  if (length(vc_cols) == 0) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No vc* columns in daily data")
    return(invisible(NULL))
  }

  if ("solar_zenith" %in% names(day_data)) {
    day_data$solar_zenith <- suppressWarnings(as.numeric(day_data$solar_zenith))
  }
  x_vals <- if ("solar_zenith" %in% names(day_data) && any(is.finite(day_data$solar_zenith))) {
    day_data$solar_zenith
  } else {
    day_data$timestamp
  }
  x_label <- if (is.numeric(x_vals)) "Solar Zenith [deg]" else "Time (UTC)"
  order_idx <- if (is.numeric(x_vals)) order(x_vals) else order(day_data$timestamp)
  day_data <- day_data[order_idx, ]
  x_vals <- x_vals[order_idx]
  y_vals <- unlist(day_data[, vc_cols, drop = FALSE], use.names = FALSE)
  y_vals <- as.numeric(y_vals[is.finite(y_vals)])

  if (length(y_vals) == 0) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "Daily column data is all NA")
    return(invisible(NULL))
  }

  y_rng <- range(y_vals, na.rm = TRUE)
  if (y_rng[1] == y_rng[2]) {
    pad <- ifelse(y_rng[1] == 0, 1, abs(y_rng[1]) * 0.1)
    y_rng <- c(y_rng[1] - pad, y_rng[2] + pad)
  }

  colors <- qualitative_hcl(n = length(vc_cols), h = c(-174, 123), c = 46, l = 40)

  plot(
    x_vals,
    day_data[[vc_cols[1]]],
    type = "n",
    ylim = y_rng,
    xlab = x_label,
    ylab = "Column amount [mol/m^2]",
    main = title
  )
  grid()
  abline(h = 0, lty = 2, col = "lightgrey", lwd = 2)

  for (i in seq_along(vc_cols)) {
    col_data <- day_data[[vc_cols[i]]]
    lines(x_vals, col_data, type = "o", pch = 16, col = colors[i])
  }

  legend("topleft", legend = vc_cols, pch = 16, lty = 1, bty = "n", col = colors)
}

pick_best_vertical_column <- function(day_data, vc_candidates = c("vc", "vc_tot", "vc_raw")) {
  if (is.null(day_data) || nrow(day_data) == 0) {
    return(NULL)
  }

  vc_cols <- vc_candidates[vc_candidates %in% names(day_data)]
  if (length(vc_cols) == 0 || !("amf" %in% names(day_data))) {
    return(NULL)
  }

  best_vc <- NULL
  best_n <- -1L
  amf_vals <- suppressWarnings(as.numeric(day_data$amf))
  for (vc_col in vc_cols) {
    vc_vals <- suppressWarnings(as.numeric(day_data[[vc_col]]))
    n_ok <- sum(is.finite(vc_vals) & is.finite(amf_vals) & amf_vals != 0)
    if (n_ok > best_n) {
      best_n <- n_ok
      best_vc <- vc_col
    }
  }

  if (!is.character(best_vc) || length(best_vc) != 1 || !nzchar(best_vc)) {
    return(NULL)
  }

  best_vc
}

prepare_daily_difference_data <- function(day_data, intercept_value = NA_real_) {
  if (is.null(day_data) || nrow(day_data) == 0) {
    return(NULL)
  }

  best_vc <- pick_best_vertical_column(day_data)
  if (is.null(best_vc) || !("vc_trop_int" %in% names(day_data))) {
    return(NULL)
  }

  timestamp <- if ("timestamp" %in% names(day_data)) {
    as.POSIXct(day_data$timestamp, tz = "UTC")
  } else {
    rep(as.POSIXct(NA, tz = "UTC"), nrow(day_data))
  }
  solar_zenith <- if ("solar_zenith" %in% names(day_data)) {
    suppressWarnings(as.numeric(day_data$solar_zenith))
  } else {
    rep(NA_real_, nrow(day_data))
  }

  original_vc <- suppressWarnings(as.numeric(day_data[[best_vc]]))
  trop_int_vc <- suppressWarnings(as.numeric(day_data$vc_trop_int))
  amf_vals <- suppressWarnings(as.numeric(day_data$amf))

  adjustment <- rep(NA_real_, length(amf_vals))
  if (is.finite(intercept_value)) {
    adjustment <- -as.numeric(intercept_value) / amf_vals
  }
  adjustment[!is.finite(adjustment)] <- NA_real_

  adjusted_vc <- original_vc + adjustment

  data.frame(
    timestamp = timestamp,
    solar_zenith = solar_zenith,
    original_vc = original_vc,
    adjusted_vc = adjusted_vc,
    trop_int_vc = trop_int_vc,
    raw_difference = original_vc - trop_int_vc,
    adjusted_difference = adjusted_vc - trop_int_vc,
    stringsAsFactors = FALSE,
    row.names = NULL
  ) %>%
    mutate(vc_col = best_vc)
}

get_daily_x_data <- function(df, prefer_solar_zenith = TRUE) {
  if (is.null(df) || nrow(df) == 0) {
    return(list(
      x = numeric(0),
      xlab = "Time (UTC)",
      order_idx = integer(0),
      use_numeric_axis = FALSE
    ))
  }

  if (isTRUE(prefer_solar_zenith) && "solar_zenith" %in% names(df)) {
    x_vals <- suppressWarnings(as.numeric(df$solar_zenith))
    if (any(is.finite(x_vals))) {
      return(list(
        x = x_vals,
        xlab = "Solar Zenith [deg]",
        order_idx = order(x_vals),
        use_numeric_axis = TRUE
      ))
    }
  }

  x_vals <- if ("timestamp" %in% names(df)) as.POSIXct(df$timestamp, tz = "UTC") else seq_len(nrow(df))
  list(
    x = x_vals,
    xlab = "Time (UTC)",
    order_idx = order(x_vals),
    use_numeric_axis = FALSE
  )
}

plot_daily_adjusted_column <- function(day_data, intercept_value, selected_day = NULL, title = "Adjusted Daily Vertical Column") {
  if (is.null(day_data) || nrow(day_data) == 0) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No daily column data found")
    return(invisible(NULL))
  }

  day_data$timestamp <- as.POSIXct(day_data$timestamp, tz = "UTC")

  if (!("amf" %in% names(day_data))) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No AMF column available for adjustment")
    return(invisible(NULL))
  }

  if (!is.finite(intercept_value)) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No finite intercept available for this day")
    return(invisible(NULL))
  }

  prepared <- prepare_daily_difference_data(day_data, intercept_value = intercept_value)
  if (is.null(prepared)) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No suitable vc/vc_trop_int data found in daily data")
    return(invisible(NULL))
  }

  best_vc <- unique(prepared$vc_col)[1]
  x_info <- get_daily_x_data(prepared, prefer_solar_zenith = TRUE)
  prepared <- prepared[x_info$order_idx, , drop = FALSE]
  keep <- is.finite(prepared$original_vc) &
    is.finite(prepared$adjusted_vc) &
    is.finite(as.numeric(if (x_info$use_numeric_axis) prepared$solar_zenith else prepared$timestamp))
  if (!any(keep)) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "Adjusted vertical column is all NA")
    return(invisible(NULL))
  }

  x <- if (x_info$use_numeric_axis) prepared$solar_zenith[keep] else prepared$timestamp[keep]
  y_original <- prepared$original_vc[keep]
  y_adjusted <- prepared$adjusted_vc[keep]
  y_trop_int <- prepared$trop_int_vc[keep]
  y_rng <- range(c(y_original, y_adjusted, y_trop_int), na.rm = TRUE)
  if (y_rng[1] == y_rng[2]) {
    pad <- ifelse(y_rng[1] == 0, 1, abs(y_rng[1]) * 0.1)
    y_rng <- c(y_rng[1] - pad, y_rng[2] + pad)
  }

  plot(
    x,
    y_original,
    type = "o",
    pch = 16,
    col = "grey35",
    ylim = y_rng,
    xlab = x_info$xlab,
    ylab = "Vertical Column [mol/m^2]",
    main = title
  )
  grid()
  abline(h = 0, lty = 2, col = "grey50", lwd = 2)
  lines(x, y_adjusted, type = "o", pch = 16, col = "firebrick", lwd = 2)
  if (any(is.finite(y_trop_int))) {
    lines(x, y_trop_int, type = "o", pch = 16, col = "steelblue", lwd = 2)
  }

  legend_labels <- c(
    sprintf("Original (%s)", best_vc),
    sprintf("Adjusted (%s - intercept/amf)", best_vc)
  )
  legend_cols <- c("grey35", "firebrick")
  legend_pch <- c(16, 16)
  legend_lty <- c(1, 1)
  if (any(is.finite(y_trop_int))) {
    legend_labels <- c(legend_labels, "Interpolated vc_trop_int")
    legend_cols <- c(legend_cols, "steelblue")
    legend_pch <- c(legend_pch, 16)
    legend_lty <- c(legend_lty, 1)
  }
  legend_labels <- c(legend_labels, sprintf("Intercept = %.3e mol/m^2", intercept_value))
  legend_cols <- c(legend_cols, NA)
  legend_pch <- c(legend_pch, NA)
  legend_lty <- c(legend_lty, NA)
  legend(
    "topleft",
    legend = legend_labels,
    col = legend_cols,
    pch = legend_pch,
    lty = legend_lty,
    bty = "n"
  )
}

plot_daily_difference <- function(day_data, intercept_value, selected_day = NULL, title = "Daily VC Difference") {
  if (is.null(day_data) || nrow(day_data) == 0) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No daily column data found")
    return(invisible(NULL))
  }

  prepared <- prepare_daily_difference_data(day_data, intercept_value = intercept_value)
  if (is.null(prepared)) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No suitable vc/vc_trop_int data found in daily data")
    return(invisible(NULL))
  }

  x_info <- get_daily_x_data(prepared, prefer_solar_zenith = TRUE)
  prepared <- prepared[x_info$order_idx, , drop = FALSE]
  keep <- is.finite(as.numeric(if (x_info$use_numeric_axis) prepared$solar_zenith else prepared$timestamp)) &
    is.finite(prepared$raw_difference)
  if (!any(keep)) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "Daily difference data is all NA")
    return(invisible(NULL))
  }

  x <- if (x_info$use_numeric_axis) prepared$solar_zenith[keep] else prepared$timestamp[keep]
  y_raw <- prepared$raw_difference[keep]
  best_vc <- unique(prepared$vc_col)[1]

  y_rng <- range(y_raw, na.rm = TRUE)
  if (y_rng[1] == y_rng[2]) {
    pad <- ifelse(y_rng[1] == 0, 1, abs(y_rng[1]) * 0.1)
    y_rng <- c(y_rng[1] - pad, y_rng[2] + pad)
  }

  plot(
    x,
    y_raw,
    type = "o",
    pch = 16,
    col = "grey35",
    ylim = y_rng,
    xlab = x_info$xlab,
    ylab = "Difference [mol/m^2]",
    main = title
  )
  grid()
  abline(h = 0, lty = 2, col = "grey50", lwd = 2)

  legend_labels <- sprintf("Raw (%s - vc_trop_int)", best_vc)
  legend_cols <- "grey35"
  legend_pch <- 16
  legend_lty <- 1

  if (is.finite(intercept_value)) {
    legend_labels <- c(legend_labels, sprintf("Intercept = %.3e mol/m^2", intercept_value))
    legend_cols <- c(legend_cols, NA)
    legend_pch <- c(legend_pch, NA)
    legend_lty <- c(legend_lty, NA)
  }

  legend("topleft", legend = legend_labels, col = legend_cols, pch = legend_pch, lty = legend_lty, bty = "n")
}

compute_daily_difference_series <- function(summary_df, day_dir = NULL) {
  if (is.null(summary_df) || nrow(summary_df) == 0 || !("date" %in% names(summary_df))) {
    return(NULL)
  }

  out <- summary_df %>%
    transmute(
      date = as.Date(date),
      raw_difference = if ("median_diff" %in% names(summary_df)) as.numeric(median_diff) else NA_real_
    ) %>%
    distinct(date, .keep_all = TRUE) %>%
    arrange(date)

  if (is.null(day_dir) || length(day_dir) == 0 || is.na(day_dir[1]) || !dir.exists(day_dir[1])) {
    return(out %>% select(date, raw_difference))
  }

  day_data <- tryCatch(
    open_dataset(day_dir[1], format = "parquet") %>%
      select(any_of(c("date", "timestamp", "solar_zenith", "vc", "vc_tot", "vc_raw", "vc_trop_int", "amf"))) %>%
      collect(),
    error = function(e) NULL
  )
  if (is.null(day_data) || nrow(day_data) == 0 || !("date" %in% names(day_data))) {
    return(out %>% select(date, raw_difference))
  }

  day_data$date <- as.Date(day_data$date)

  raw_rows <- lapply(out$date[!is.na(out$date)], function(day_value) {
    day_subset <- day_data %>% filter(as.Date(date) == as.Date(day_value))
    if (nrow(day_subset) == 0) {
      return(NULL)
    }

    prepared <- prepare_daily_difference_data(day_subset)
    if (is.null(prepared)) {
      return(data.frame(date = as.Date(day_value), raw_difference = NA_real_))
    }

    raw_vals <- prepared$raw_difference
    raw_vals <- raw_vals[is.finite(raw_vals)]
    data.frame(
      date = as.Date(day_value),
      raw_difference = if (length(raw_vals) > 0) median(raw_vals) else NA_real_
    )
  })

  raw_df <- bind_rows(raw_rows)
  if (nrow(raw_df) == 0 || !all(c("date", "raw_difference") %in% names(raw_df))) {
    return(out %>% select(date, raw_difference))
  }

  out %>%
    left_join(raw_df, by = "date", suffix = c("", "_computed")) %>%
    mutate(raw_difference = coalesce(raw_difference, raw_difference_computed)) %>%
    select(date, raw_difference)
}

plot_daily_difference_series <- function(series_data, selected_day = NULL, title = "Daily Difference Time Series") {
  if (is.null(series_data) || nrow(series_data) == 0) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "No daily difference series available")
    return(invisible(NULL))
  }

  series_data <- series_data %>%
    mutate(date = as.Date(date)) %>%
    arrange(date)

  y_vals <- series_data$raw_difference
  y_vals <- y_vals[is.finite(y_vals)]
  if (length(y_vals) == 0) {
    plot.new()
    title(main = title)
    text(0.5, 0.5, "Daily difference series is all NA")
    return(invisible(NULL))
  }

  y_rng <- range(y_vals, na.rm = TRUE)
  if (y_rng[1] == y_rng[2]) {
    pad <- ifelse(y_rng[1] == 0, 1, abs(y_rng[1]) * 0.1)
    y_rng <- c(y_rng[1] - pad, y_rng[2] + pad)
  }

  plot(
    series_data$date,
    series_data$raw_difference,
    type = "o",
    pch = 16,
    col = "grey35",
    ylim = y_rng,
    xlab = "Date",
    ylab = "Median difference [mol/m^2]",
    main = title
  )
  grid()
  abline(h = 0, lty = 2, col = "grey50", lwd = 2)

  legend_labels <- "Raw median (vc - vc_trop_int)"
  legend_cols <- "grey35"
  legend_pch <- 16
  legend_lty <- 1

  if (!is.null(selected_day) && !is.na(as.Date(selected_day))) {
    abline(v = as.Date(selected_day), col = "red", lwd = 2, lty = 2)
    legend_labels <- c(legend_labels, "Selected day")
    legend_cols <- c(legend_cols, "red")
    legend_pch <- c(legend_pch, NA)
    legend_lty <- c(legend_lty, 2)
  }

  legend("topleft", legend = legend_labels, col = legend_cols, pch = legend_pch, lty = legend_lty, bty = "n")
}

read_day_partition <- function(base_dir, day_value) {
  if (is.null(base_dir) || length(base_dir) == 0 || is.na(base_dir[1]) || !nzchar(base_dir[1])) {
    return(NULL)
  }

  day_chr <- as.character(as.Date(day_value))
  part_dir <- file.path(base_dir[1], paste0("date=", day_chr))
  if (!dir.exists(part_dir)) {
    return(NULL)
  }

  files <- list.files(part_dir, pattern = "\\.parquet$", full.names = TRUE)
  if (length(files) == 0) {
    return(NULL)
  }

  out <- bind_rows(lapply(files, read_parquet))
  if ("date" %in% names(out)) {
    out$date <- as.Date(out$date)
  } else {
    out$date <- rep(as.Date(day_chr), nrow(out))
  }
  if ("timestamp" %in% names(out)) {
    out$timestamp <- as.POSIXct(out$timestamp, tz = "UTC")
  }
  out
}

pick_single_choice <- function(current, choices) {
  choices <- choices[!is.na(choices)]
  if (length(choices) == 0) {
    return(character(0))
  }
  if (!is.null(current) && length(current) == 1 && current %in% choices) {
    return(current)
  }
  choices[[1]]
}

pick_optional_single_choice <- function(current, choices) {
  choices <- choices[!is.na(choices)]
  if (!is.null(current) && length(current) == 1 && current %in% choices) {
    return(current)
  }
  character(0)
}

pick_multi_choice <- function(current, choices) {
  choices <- choices[!is.na(choices)]
  if (length(choices) == 0) {
    return(character(0))
  }
  selected <- intersect(as.character(current %||% character(0)), as.character(choices))
  if (length(selected) > 0) {
    return(selected)
  }
  as.character(choices)
}

pick_preferred_window <- function(choices, preferred = "window_30") {
  choices <- as.character(choices %||% character(0))
  if (length(choices) == 0) {
    return(character(0))
  }
  if (preferred %in% choices) {
    return(preferred)
  }
  choices[[1]]
}

build_dataset_id <- function(pan_id, spec_id, location, blickp_version) {
  sprintf("pan%s_spec%s_%s_%s", pan_id, spec_id, location, blickp_version)
}

filter_result_metadata <- function(meta,
                                   dataset_id = NULL,
                                   cf_vector_dir = NULL,
                                   code = NULL,
                                    window_length = NULL) {
  if (is.null(meta) || nrow(meta) == 0) {
    return(meta)
  }

  out <- meta
  if (!is.null(dataset_id) && length(dataset_id) == 1 && nzchar(dataset_id)) {
    out <- out %>% filter(dataset_id == !!dataset_id)
  }
  if (!is.null(cf_vector_dir) && length(cf_vector_dir) == 1 && nzchar(cf_vector_dir)) {
    out <- out %>% filter(cf_vector_dir == !!cf_vector_dir)
  }
  if (!is.null(code) && length(code) == 1 && nzchar(code)) {
    out <- out %>% filter(code == !!code)
  }
  if (!is.null(window_length) && length(window_length) > 0) {
    window_values <- suppressWarnings(as.numeric(window_length))
    window_values <- window_values[is.finite(window_values)]
    if (length(window_values) > 0) {
      out <- out %>% filter(window_length %in% window_values)
    }
  }

  out
}

read_partitioned_metadata <- function(result_dir, dataset_name, dir_col) {
  files <- list.files(result_dir, pattern = "\\.parquet$", recursive = TRUE, full.names = TRUE)
  if (length(files) == 0) {
    return(NULL)
  }

  sep_pattern <- "[/\\\\]"
  nested_pattern <- paste0(
    "pan(\\d+)_spec(\\d+)_([^_/\\\\]+)_([^/\\\\]+)", sep_pattern,
    "cf_vector=([^/\\\\]+)", sep_pattern,
    "code=([^/\\\\]+)", sep_pattern,
    "window_length=(\\d+)", sep_pattern,
    dataset_name, sep_pattern,
    "date=([^/\\\\]+)", sep_pattern,
    "[^/\\\\]+\\.parquet$"
  )
  rooted_pattern <- paste0(
    "cf_vector=([^/\\\\]+)", sep_pattern,
    "pan(\\d+)_spec(\\d+)_([^_/\\\\]+)_([^/\\\\]+)", sep_pattern,
    "code=([^/\\\\]+)", sep_pattern,
    "window_length=(\\d+)", sep_pattern,
    dataset_name, sep_pattern,
    "date=([^/\\\\]+)", sep_pattern,
    "[^/\\\\]+\\.parquet$"
  )
  legacy_pattern <- paste0(
    "pan(\\d+)_spec(\\d+)_([^_/\\\\]+)_([^/\\\\]+)", sep_pattern,
    "code=([^/\\\\]+)", sep_pattern,
    "window_length=(\\d+)", sep_pattern,
    dataset_name, sep_pattern,
    "date=([^/\\\\]+)", sep_pattern,
    "[^/\\\\]+\\.parquet$"
  )

  parsed_rows <- lapply(files, function(file) {
    nested <- str_match(file, nested_pattern)
    if (!is.na(nested[1, 1])) {
      return(data.frame(
        dataset_dir = dirname(dirname(file)),
        cf_vector_dir = nested[1, 6],
        pan_id = nested[1, 2],
        spec_id = nested[1, 3],
        location = nested[1, 4],
        blickp_version = nested[1, 5],
        code = nested[1, 7],
        window_length = as.numeric(nested[1, 8]),
        date = as.Date(nested[1, 9]),
        stringsAsFactors = FALSE
      ))
    }

    rooted <- str_match(file, rooted_pattern)
    if (!is.na(rooted[1, 1])) {
      return(data.frame(
        dataset_dir = dirname(dirname(file)),
        cf_vector_dir = rooted[1, 2],
        pan_id = rooted[1, 3],
        spec_id = rooted[1, 4],
        location = rooted[1, 5],
        blickp_version = rooted[1, 6],
        code = rooted[1, 7],
        window_length = as.numeric(rooted[1, 8]),
        date = as.Date(rooted[1, 9]),
        stringsAsFactors = FALSE
      ))
    }

    legacy <- str_match(file, legacy_pattern)
    if (!is.na(legacy[1, 1])) {
      return(data.frame(
        dataset_dir = dirname(dirname(file)),
        cf_vector_dir = "(none)",
        pan_id = legacy[1, 2],
        spec_id = legacy[1, 3],
        location = legacy[1, 4],
        blickp_version = legacy[1, 5],
        code = legacy[1, 6],
        window_length = as.numeric(legacy[1, 7]),
        date = as.Date(legacy[1, 8]),
        stringsAsFactors = FALSE
      ))
    }

    NULL
  })

  parsed_rows <- Filter(Negate(is.null), parsed_rows)
  if (length(parsed_rows) == 0) {
    return(NULL)
  }

  bind_rows(parsed_rows) %>%
    mutate(
      cf_vector_dir = ifelse(is.na(cf_vector_dir) | !nzchar(cf_vector_dir), "(none)", cf_vector_dir),
      dataset_id = build_dataset_id(pan_id, spec_id, location, blickp_version)
    ) %>%
    rename(!!dir_col := dataset_dir)
}

read_daily_summary_metadata <- function(result_dir) {
  summary_dirs <- list.dirs(result_dir, recursive = TRUE, full.names = TRUE)
  if (length(summary_dirs) == 0) {
    return(list(windows = NULL))
  }

  summary_dirs <- summary_dirs[grepl("daily_summary$", summary_dirs)]
  if (length(summary_dirs) == 0) {
    return(list(windows = NULL))
  }

  sep_pattern <- "[/\\\\]"
  nested_pattern <- paste0(
    "pan(\\d+)_spec(\\d+)_([^_/\\\\]+)_([^/\\\\]+)", sep_pattern,
    "cf_vector=([^/\\\\]+)", sep_pattern,
    "code=([^/\\\\]+)", sep_pattern,
    "window_length=(\\d+)", sep_pattern,
    "daily_summary$"
  )
  rooted_pattern <- paste0(
    "cf_vector=([^/\\\\]+)", sep_pattern,
    "pan(\\d+)_spec(\\d+)_([^_/\\\\]+)_([^/\\\\]+)", sep_pattern,
    "code=([^/\\\\]+)", sep_pattern,
    "window_length=(\\d+)", sep_pattern,
    "daily_summary$"
  )
  legacy_pattern <- paste0(
    "pan(\\d+)_spec(\\d+)_([^_/\\\\]+)_([^/\\\\]+)", sep_pattern,
    "code=([^/\\\\]+)", sep_pattern,
    "window_length=(\\d+)", sep_pattern,
    "daily_summary$"
  )

  parsed_rows <- lapply(summary_dirs, function(dir_path) {
    nested <- str_match(dir_path, nested_pattern)
    if (!is.na(nested[1, 1])) {
      return(data.frame(
        summary_dir = dir_path,
        cf_vector_dir = nested[1, 6],
        pan_id = nested[1, 2],
        spec_id = nested[1, 3],
        location = nested[1, 4],
        blickp_version = nested[1, 5],
        code = nested[1, 7],
        window_length = as.numeric(nested[1, 8]),
        stringsAsFactors = FALSE
      ))
    }

    rooted <- str_match(dir_path, rooted_pattern)
    if (!is.na(rooted[1, 1])) {
      return(data.frame(
        summary_dir = dir_path,
        cf_vector_dir = rooted[1, 2],
        pan_id = rooted[1, 3],
        spec_id = rooted[1, 4],
        location = rooted[1, 5],
        blickp_version = rooted[1, 6],
        code = rooted[1, 7],
        window_length = as.numeric(rooted[1, 8]),
        stringsAsFactors = FALSE
      ))
    }

    legacy <- str_match(dir_path, legacy_pattern)
    if (!is.na(legacy[1, 1])) {
      return(data.frame(
        summary_dir = dir_path,
        cf_vector_dir = "(none)",
        pan_id = legacy[1, 2],
        spec_id = legacy[1, 3],
        location = legacy[1, 4],
        blickp_version = legacy[1, 5],
        code = legacy[1, 6],
        window_length = as.numeric(legacy[1, 7]),
        stringsAsFactors = FALSE
      ))
    }

    NULL
  })

  rows <- bind_rows(Filter(Negate(is.null), parsed_rows))
  if (is.null(rows) || nrow(rows) == 0) {
    return(list(windows = NULL))
  }

  windows <- rows %>%
    mutate(
      cf_vector_dir = ifelse(is.na(cf_vector_dir) | !nzchar(cf_vector_dir), "(none)", cf_vector_dir),
      dataset_id = build_dataset_id(pan_id, spec_id, location, blickp_version)
    ) %>%
    distinct(summary_dir, dataset_id, cf_vector_dir, pan_id, spec_id, location, blickp_version, code, window_length)

  list(windows = windows)
}

read_daily_fit_metadata <- function(result_dir) {
  rows <- read_partitioned_metadata(result_dir, "daily_mle_fit", "fit_dir")
  if (is.null(rows) || nrow(rows) == 0) return(NULL)

  rows %>%
    distinct(fit_dir, dataset_id, cf_vector_dir, pan_id, spec_id, location, blickp_version, code, window_length, date)
}

read_daily_columns_metadata <- function(result_dir) {
  rows <- read_partitioned_metadata(result_dir, "daily_columns", "day_dir")
  if (is.null(rows) || nrow(rows) == 0) return(NULL)

  rows %>%
    distinct(day_dir, dataset_id, cf_vector_dir, pan_id, spec_id, location, blickp_version, code, window_length, date)
}

derive_related_dir <- function(summary_dir, target_name) {
  if (is.null(summary_dir) || length(summary_dir) == 0) {
    return(character(0))
  }
  sub("daily_summary$", target_name, summary_dir)
}


get_cf_gas_slant <- function(instrument,
                             cfVector,
                             gas,
                             spectrometer = 1,
                             endpoint = "https://hetzner.pandonia-global-network.org/query",
                             field = "Slant columns in synthetic reference spectrum for U340") {
  requireNamespace("httr", quietly = TRUE)
  requireNamespace("jsonlite", quietly = TRUE)
  requireNamespace("zoo", quietly = TRUE)
  
  query <- sprintf(
    '{
      cfVector(instrument:%s, spectrometer:%s, cfVector:"%s") {
        calibrationFiles {
          content(fields:"%s")
          cfDate
          syntRefU340
        }
      }
    }',
    instrument, spectrometer, cfVector, field
  )
  
  response <- httr::POST(
    url = endpoint,
    body = list(query = query),
    encode = "json"
  )
  
  httr::stop_for_status(response)
  
  parsed <- jsonlite::fromJSON(
    httr::content(response, as = "text", encoding = "UTF-8"),
    simplifyVector = FALSE
  )
  
  if (!is.null(parsed$errors)) {
    stop("GraphQL error: ", jsonlite::toJSON(parsed$errors, auto_unbox = TRUE))
  }
  
  cf_vectors <- parsed$data$cfVector
  
  if (length(cf_vectors) == 0) {
    stop("No cfVector result returned.")
  }
  
  number_pattern <- "[+-]?(?:\\d+\\.?\\d*|\\.\\d+)(?:[eE][+-]?\\d+)?"
  
  gas_pattern <- paste0(
    "(^|\\s)",
    gas,
    ":(",
    number_pattern,
    ")\\((",
    number_pattern,
    ")\\),(",
    number_pattern,
    ")\\((",
    number_pattern,
    ")\\)K,(OD\\d+)"
  )
  
  parse_synt_ref_u340 <- function(x) {
    empty <- list(
      synt_ref_start = as.POSIXct(NA_real_, origin = "1970-01-01", tz = "UTC"),
      synt_ref_end = as.POSIXct(NA_real_, origin = "1970-01-01", tz = "UTC"),
      synt_ref_location = NA_character_
    )
    
    if (is.null(x) || length(x) == 0 || is.na(x) || x == "") {
      return(empty)
    }
    
    pattern <- paste0(
      "^",
      "([0-9]{8}T[0-9]{6}Z)",
      "\\s+to\\s+",
      "([0-9]{8}T[0-9]{6}Z)",
      "(?:\\s+at\\s+(.*))?",
      "$"
    )
    
    m <- regmatches(x, regexec(pattern, x, perl = TRUE))[[1]]
    
    if (length(m) == 0) {
      return(empty)
    }
    
    list(
      synt_ref_start = as.POSIXct(m[2], format = "%Y%m%dT%H%M%SZ", tz = "UTC"),
      synt_ref_end = as.POSIXct(m[3], format = "%Y%m%dT%H%M%SZ", tz = "UTC"),
      synt_ref_location = ifelse(length(m) >= 4 && nzchar(m[4]), m[4], NA_character_)
    )
  }
  
  out <- list()
  k <- 1
  
  for (i in seq_along(cf_vectors)) {
    calibration_files <- cf_vectors[[i]]$calibrationFiles
    
    if (length(calibration_files) == 0) {
      next
    }
    
    for (j in seq_along(calibration_files)) {
      cf_file <- calibration_files[[j]]
      
      content_string <- cf_file$content[[field]]
      cf_date_raw <- cf_file$cfDate
      synt_ref_raw <- cf_file$syntRefU340
      
      if (is.null(content_string) || is.null(cf_date_raw)) {
        next
      }
      
      match <- regmatches(
        content_string,
        regexec(gas_pattern, content_string, perl = TRUE)
      )[[1]]
      
      if (length(match) == 0) {
        values <- data.frame(
          slant_column = NA_real_,
          slant_column_uncertainty = NA_real_,
          effective_temperature = NA_real_,
          effective_temperature_uncertainty = NA_real_
        )
      } else {
        values <- data.frame(
          slant_column = as.numeric(match[3]),
          slant_column_uncertainty = as.numeric(match[4]),
          effective_temperature = as.numeric(match[5]),
          effective_temperature_uncertainty = as.numeric(match[6])
        )
      }
      
      synt_ref <- parse_synt_ref_u340(synt_ref_raw)
      
      out[[k]] <- data.frame(
        cfDate = as.Date(cf_date_raw, format = "%Y%m%d"),
        values,
        synt_ref_start = synt_ref$synt_ref_start,
        synt_ref_end = synt_ref$synt_ref_end,
        synt_ref_duration_hours = as.numeric(
          difftime(
            synt_ref$synt_ref_end,
            synt_ref$synt_ref_start,
            units = "hours"
          )
        ),
        synt_ref_location = synt_ref$synt_ref_location,
        stringsAsFactors = FALSE
      )
      
      k <- k + 1
    }
  }
  
  if (length(out) == 0) {
    stop("No calibration files with usable cfDate/content found.")
  }
  
  result_df <- do.call(rbind, out)
  result_df <- result_df[order(result_df$cfDate), ]
  rownames(result_df) <- NULL
  
  value_cols <- c(
    "slant_column",
    "slant_column_uncertainty",
    "effective_temperature",
    "effective_temperature_uncertainty"
  )
  
  metadata_cols <- c(
    "cfDate",
    "synt_ref_start",
    "synt_ref_end",
    "synt_ref_duration_hours",
    "synt_ref_location"
  )
  
  values_zoo <- zoo::zoo(
    result_df[, value_cols],
    order.by = result_df$cfDate
  )
  
  metadata <- result_df[, metadata_cols]
  rownames(metadata) <- NULL
  
  list(
    values = values_zoo,
    metadata = metadata
  )
}

#---- UI ----
ui <- fluidPage(
  titlePanel("Running Calibration Viewer"),
  sidebarLayout(
    sidebarPanel(
      width = 2,
      textInput("dir", "Results Directory", value = results_directory),
      selectInput("dataset_id", "Dataset", choices = character(0)),
      selectInput("cf_vector", "CF Vector", choices = character(0)),
      selectInput("code", "R Code", choices = character(0)),
      selectInput("windows", "Window Length(s)", choices = character(0), multiple = TRUE),
      selectInput("parameter", "Parameter", choices = c("intercept", "slope","ndata","groups","ndays","amf_dependency")),
      conditionalPanel(
        condition = "input.parameter == 'intercept'",
        wellPanel(
          checkboxInput("show_calibration_uncertainty", "Show initial calibration uncertainty", value = FALSE),
          conditionalPanel(
            condition = "input.show_calibration_uncertainty",
            numericInput("uncertainty_sigma", "Sigma level", value = 2, min = 0.1, step = 0.5)
          )
        )
      ),
      numericInput("min_ndays", "Min ndays (summary filter)", value = 0, min = 0, step = 1),
      numericInput("min_groups", "Min groups (summary filter)", value = 0, min = 0, step = 1),
      actionButton("load", "Load & Plot"),
      uiOutput("loading_status"),
      selectInput("fit_window", "Window Length for Daily MLE Fit", choices = NULL),
      selectInput("bp_window", "Window Length for Breakpoint Analysis", choices = NULL),
      actionButton("analyze_bp", "Analyze Breakpoints")
    ),
    mainPanel(
      width = 10,
      plotOutput("resultPlot"),#, height = "650px"),
      h4("Daily Difference Time Series"),
      plotOutput("dailyDifferenceSeriesPlot"),
      fluidRow(
        column(12,
        uiOutput("fit_day_slider"))
      ),
      fluidRow(
        column(
          width = 3,
          h4("Daily MLE Fit"),
          plotOutput("dailyFitPlot")
        ),
        column(
          width = 3,
          h4("Daily Columns"),
          plotOutput("dailyColumnsPlot")
        ),
        column(
          width = 3,
          h4("Adjusted Daily Vertical Column"),
          plotOutput("dailyAdjustedColumnPlot")
        ),
        column(
          width = 3,
          h4("Daily Difference to vc_trop_int"),
          plotOutput("dailyDifferencePlot")
        )
      ),
      hr(),
      h4("Breakpoint Analysis"),
      plotOutput("bpPlot"),#, height = "400px")
      hr(),
      uiOutput("calibration_uncertainty_panel")
      )
  )
)

#---- SERVER ----
server <- function(input, output, session) {
  status_message <- reactiveVal("Select a dataset and click 'Load & Plot'.")

  output$loading_status <- renderUI({
    tags$div(
      style = "margin-top: 10px; margin-bottom: 12px; color: #666; font-size: 0.9em;",
      status_message()
    )
  })
  
  # Scan directory and extract metadata
  metadata <- reactive({
    req(input$dir)
    status_message("Scanning available datasets...")
    withProgress(message = "Scanning results", detail = "Scanning summary windows...", value = 0, {
      incProgress(0.4)
      out <- read_daily_summary_metadata(input$dir)
      incProgress(0.6)
      out
    })
  })

  all_windows <- reactive({
    metadata()$windows
  })

  dataset_choices <- reactive({
    meta <- all_windows()
    req(meta)
    sort(unique(meta$dataset_id))
  })

  observe({
    choices <- tryCatch(dataset_choices(), error = function(e) character(0))
    if (length(choices) == 0) {
      status_message("No datasets found in the selected results directory.")
    } else if (!isTruthy(input$dataset_id) || input$load < 1) {
      status_message(sprintf("%d dataset(s) found. Select one and click 'Load & Plot'.", length(choices)))
    }
  })

  cf_vector_choices <- reactive({
    req(isTruthy(input$dataset_id))
    meta <- filter_result_metadata(all_windows(), dataset_id = input$dataset_id)
    req(meta)
    sort(unique(meta$cf_vector_dir))
  })

  code_choices <- reactive({
    req(isTruthy(input$dataset_id), isTruthy(input$cf_vector))
    meta <- filter_result_metadata(
      all_windows(),
      dataset_id = input$dataset_id,
      cf_vector_dir = input$cf_vector
    )
    req(meta)
    sort(unique(meta$code))
  })

  window_choices <- reactive({
    req(isTruthy(input$dataset_id), isTruthy(input$cf_vector), isTruthy(input$code))
    meta <- filter_result_metadata(
      all_windows(),
      dataset_id = input$dataset_id,
      cf_vector_dir = input$cf_vector,
      code = input$code
    )
    req(meta)
    sort(unique(meta$window_length))
  })

  observe({
    choices <- dataset_choices()
    updateSelectInput(
      session,
      "dataset_id",
      choices = choices,
      selected = pick_optional_single_choice(input$dataset_id, choices)
    )
  })

  observe({
    choices <- tryCatch(cf_vector_choices(), error = function(e) character(0))
    updateSelectInput(
      session,
      "cf_vector",
      choices = choices,
      selected = pick_optional_single_choice(input$cf_vector, choices)
    )
  })

  observe({
    choices <- tryCatch(code_choices(), error = function(e) character(0))
    updateSelectInput(
      session,
      "code",
      choices = choices,
      selected = pick_optional_single_choice(input$code, choices)
    )
  })

  observe({
    choices <- tryCatch(window_choices(), error = function(e) character(0))
    updateSelectInput(
      session,
      "windows",
      choices = choices,
      selected = pick_multi_choice(input$windows, choices)
    )
  })

  loaded_selection <- eventReactive(input$load, {
    req(
      isTruthy(input$dataset_id),
      isTruthy(input$cf_vector),
      isTruthy(input$code),
      length(input$windows) > 0
    )
    status_message("Loading selected dataset...")
    list(
      dataset_id = input$dataset_id,
      cf_vector = input$cf_vector,
      code = input$code,
      windows = as.character(input$windows)
    )
  }, ignoreInit = TRUE)

  selected_windows <- reactive({
    selection <- loaded_selection()
    req(selection)
    meta <- filter_result_metadata(
      all_windows(),
      dataset_id = selection$dataset_id,
      cf_vector_dir = selection$cf_vector,
      code = selection$code
    )
    req(meta)
    filter_result_metadata(meta, window_length = selection$windows)
  })

  fit_window_length <- reactive({
    req(input$fit_window)
    as.numeric(sub("window_", "", input$fit_window))
  })

  fit_dirs <- reactive({
    req(selected_windows(), fit_window_length())
    meta <- selected_windows() %>%
      filter(window_length == fit_window_length())
    if (nrow(meta) == 0 || !"summary_dir" %in% names(meta)) {
      return(character(0))
    }
    dirs <- unique(derive_related_dir(meta$summary_dir, "daily_mle_fit"))
    dirs[dir.exists(dirs)]
  })

  day_dirs <- reactive({
    req(selected_windows(), fit_window_length())
    meta <- selected_windows() %>%
      filter(window_length == fit_window_length())
    if (nrow(meta) == 0 || !"summary_dir" %in% names(meta)) {
      return(character(0))
    }
    dirs <- unique(derive_related_dir(meta$summary_dir, "daily_columns"))
    dirs[dir.exists(dirs)]
  })

  loaded_series_days <- reactive({
    dset <- dataset()
    req(dset)

    selected_name <- input$fit_window
    if (isTruthy(selected_name) && selected_name %in% names(dset)) {
      dates <- as.Date(dset[[selected_name]]$date)
    } else {
      dates <- do.call(c, lapply(dset, function(df) as.Date(df$date)))
    }

    dates <- sort(unique(dates[!is.na(dates)]))
    dates
  })

  output$fit_day_slider <- renderUI({
    days <- loaded_series_days()

    if (length(days) == 0) {
      return(helpText("Load a dataset to choose an MLE fit day."))
    }

    selected_day <- isolate(input$fit_day)
    if (is.null(selected_day) || !as.Date(selected_day) %in% days) {
      if (is.null(selected_day)) {
        selected_day <- max(days)
      } else {
        selected_day <- days[which.min(abs(as.numeric(days - as.Date(selected_day))))]
      }
    }

    tagList(
      fluidRow(
        column(
          width = 6,
          actionButton("fit_day_prev", "Previous day")
        ),
        column(
          width = 6,
          actionButton("fit_day_next", "Next day")
        )
      ),
      sliderInput(
        "fit_day",
        "MLE Fit Day",
        min = min(days),
        max = max(days),
        value = selected_day,
        step = 1,
        timeFormat = "%Y-%m-%d",
        animate = animationOptions(interval = 800, loop = FALSE)
      )
    )
  })

  move_fit_day <- function(offset) {
    days <- loaded_series_days()
    if (length(days) == 0 || is.null(input$fit_day)) {
      return(invisible(NULL))
    }

    current_day <- as.Date(input$fit_day)
    current_idx <- which(days == current_day)
    if (length(current_idx) == 0) {
      current_idx <- which.min(abs(as.numeric(days - current_day)))
    }

    target_idx <- min(max(current_idx[1] + offset, 1), length(days))
    updateSliderInput(session, "fit_day", value = days[target_idx])
    invisible(NULL)
  }

  observeEvent(input$fit_day_prev, {
    move_fit_day(-1)
  }, ignoreInit = TRUE)

  observeEvent(input$fit_day_next, {
    move_fit_day(1)
  }, ignoreInit = TRUE)

  #---- Load parquet datasets according to filters ----
  dataset <- reactive({
    req(selected_windows())
    meta <- selected_windows()
    if (nrow(meta) == 0) return(NULL)

    status_message(sprintf("Loading dataset '%s'...", loaded_selection()$dataset_id))
    showNotification("Loading selected dataset...", id = "dataset_load", duration = NULL, closeButton = FALSE)
    on.exit(removeNotification("dataset_load"), add = TRUE)

    result_list <- withProgress(message = "Loading dataset", detail = "Reading window summaries...", value = 0, {
      result_list <- list()
      n_meta <- nrow(meta)
      for (i in seq_len(n_meta)) {
        incProgress(
          amount = 1 / max(n_meta, 1),
          detail = sprintf("Loading window %s (%d of %d)", meta$window_length[i], i, n_meta)
        )
        df <- open_dataset(meta$summary_dir[i], format = "parquet") %>%
          collect()
        min_ndays <- if (is.null(input$min_ndays) || is.na(input$min_ndays)) 0 else as.numeric(input$min_ndays)
        min_groups <- if (is.null(input$min_groups) || is.na(input$min_groups)) 0 else as.numeric(input$min_groups)
        if ("ndays" %in% names(df) && is.finite(min_ndays) && min_ndays > 0) {
          df <- df %>%
            filter(is.finite(as.numeric(ndays)) & as.numeric(ndays) >= min_ndays)
        }
        if ("groups" %in% names(df) && is.finite(min_groups) && min_groups > 0) {
          df <- df %>%
            filter(is.finite(as.numeric(groups)) & as.numeric(groups) >= min_groups)
        }
        df <- df %>% arrange(as.Date(date))
        result_list[[paste0("window_", meta$window_length[i])]] <- df
      }
      result_list
    })

    status_message(sprintf(
      "Loaded %s / %s / %s with %d window(s).",
      loaded_selection()$dataset_id,
      loaded_selection()$cf_vector,
      loaded_selection()$code,
      length(result_list)
    ))
    result_list
  })
  
  
  # --- Update breakpoint window choices after dataset is loaded ---
  observeEvent(dataset(), {
    req(dataset())
    default_window <- pick_preferred_window(names(dataset()), preferred = "window_30")
    updateSelectInput(session, "bp_window",
                      choices = names(dataset()),
            selected = default_window)
    updateSelectInput(session, "fit_window",
                      choices = names(dataset()),
                      selected = default_window)
  })

  selected_calibration_context <- reactive({
    selection <- loaded_selection()
    req(selection)
    dset <- dataset()
    req(dset)

    available <- Filter(function(df) !is.null(df) && nrow(df) > 0, dset)
    if (length(available) == 0) {
      return(NULL)
    }

    row <- available[[1]][1, , drop = FALSE]
    pan_id <- as.character(row$pan_id[1] %||% NA_character_)
    spectrometer <- suppressWarnings(as.integer(row$spectrometer[1] %||% 1L))
    cf_vector <- as.character(row$cf_vector[1] %||% NA_character_)
    gas <- as.character(row$species[1] %||% NA_character_)

    if (!nzchar(pan_id) || is.na(pan_id) || !nzchar(cf_vector) || is.na(cf_vector) || !nzchar(gas) || is.na(gas)) {
      return(NULL)
    }

    list(
      pan_id = pan_id,
      spectrometer = if (is.finite(spectrometer)) spectrometer else 1L,
      cf_vector = cf_vector,
      gas = gas
    )
  })

  calibration_uncertainty_response <- reactive({
    if (!isTRUE(input$show_calibration_uncertainty) || !identical(input$parameter, "intercept")) {
      return(NULL)
    }

    ctx <- selected_calibration_context()
    if (is.null(ctx)) {
      return(NULL)
    }

    showNotification("Loading calibration uncertainty...", id = "cf_uncertainty_load", duration = NULL, closeButton = FALSE)
    on.exit(removeNotification("cf_uncertainty_load"), add = TRUE)
    cf_data <- withProgress(message = "Loading calibration uncertainty", detail = "Reading calibration factors...", value = 0, {
      incProgress(0.3)
      out <- tryCatch(
        get_cf_gas_slant(
          ctx$pan_id,
          ctx$cf_vector,
          ctx$gas,
          spectrometer = ctx$spectrometer
        ),
        error = function(e) NULL
      )
      incProgress(0.7)
      out
    })

    if (is.null(cf_data) || is.null(cf_data$values) || is.null(cf_data$metadata)) {
      return(NULL)
    }

    values_df <- as.data.frame(cf_data$values)
    if (nrow(values_df) == 0) {
      return(NULL)
    }

    values_df$cfDate <- as.Date(zoo::index(cf_data$values))
    metadata_df <- as.data.frame(cf_data$metadata)
    metadata_df$cfDate <- as.Date(metadata_df$cfDate)

    metadata_df %>%
      left_join(values_df, by = "cfDate") %>%
      transmute(
        cal_date = cfDate,
        slant_column = as.numeric(slant_column),
        scref_uncert = as.numeric(slant_column_uncertainty),
        effective_temperature = as.numeric(effective_temperature),
        effective_temperature_uncertainty = as.numeric(effective_temperature_uncertainty),
        synt_ref_start = as.POSIXct(synt_ref_start, tz = "UTC"),
        synt_ref_end = as.POSIXct(synt_ref_end, tz = "UTC"),
        synt_ref_duration_hours = as.numeric(synt_ref_duration_hours),
        synt_ref_location = as.character(synt_ref_location)
      ) %>%
      filter(!is.na(cal_date)) %>%
      arrange(cal_date)
  })

  calibration_uncertainty <- reactive({
    response <- calibration_uncertainty_response()
    if (is.null(response) || nrow(response) == 0) {
      return(NULL)
    }
    response %>%
      filter(is.finite(scref_uncert), scref_uncert >= 0) %>%
      distinct(cal_date, .keep_all = TRUE)
  })

  output$calibration_uncertainty_panel <- renderUI({
    if (!isTRUE(input$show_calibration_uncertainty) || !identical(input$parameter, "intercept")) {
      return(NULL)
    }

    tagList(
      h4("Queried Calibration Files"),
      tableOutput("calibrationUncertaintyTable")
    )
  })

  output$calibrationUncertaintyTable <- renderTable({
    response <- calibration_uncertainty_response()
    if (is.null(response) || nrow(response) == 0) {
      return(NULL)
    }

    response %>%
      transmute(
        cal_date = as.character(cal_date),
        slant_column = format_scientific_value(slant_column, digits = 6),
        slant_column_uncertainty = format_scientific_value(scref_uncert, digits = 6),
        effective_temperature = format_fixed_value(effective_temperature, digits = 2),
        effective_temperature_uncertainty = format_fixed_value(effective_temperature_uncertainty, digits = 2),
        synt_ref_start = as.character(synt_ref_start),
        synt_ref_end = as.character(synt_ref_end),
        synt_ref_duration_hours = format_fixed_value(synt_ref_duration_hours, digits = 2),
        synt_ref_location = synt_ref_location
      )
  }, striped = TRUE, bordered = TRUE, spacing = "xs")

  uncertainty_sigma_level <- reactive({
    sigma_level <- suppressWarnings(as.numeric(input$uncertainty_sigma))
    if (!is.finite(sigma_level) || sigma_level <= 0) {
      return(1)
    }
    sigma_level
  })
  
  
  # --- Plot ---
  output$resultPlot <- renderPlot({
    req(dataset(),input$fit_day)
    selection <- loaded_selection()
    req(selection)
    plot_title <- paste(selection$dataset_id, selection$cf_vector, selection$code)
    withProgress(message = "Rendering overview", detail = "Drawing selected parameter...", value = 0, {
      incProgress(0.4)
      plot_result(
        dataset(),
        result = input$parameter,
        title = plot_title,
        calibration_uncertainty = calibration_uncertainty(),
        show_calibration_uncertainty = isTRUE(input$show_calibration_uncertainty),
        sigma_level = uncertainty_sigma_level()
      )
      incProgress(0.6)
      abline(v=input$fit_day,col=2,lwd=3,lty=2)
    })
  })

  daily_summary_row <- reactive({
    req(dataset(), input$fit_window, input$fit_day)
    dset <- dataset()[[input$fit_window]]
    if (is.null(dset) || nrow(dset) == 0) return(NULL)
    day_df <- dset %>%
      filter(as.Date(date) == as.Date(input$fit_day))
    if (nrow(day_df) == 0) return(NULL)
    idx <- 1
    if ("ndata" %in% names(day_df)) {
      nvals <- suppressWarnings(as.numeric(day_df$ndata))
      if (any(is.finite(nvals))) {
        idx <- which.max(ifelse(is.finite(nvals), nvals, -Inf))
      }
    }
    day_df[idx, , drop = FALSE]
  })



  output$dailyFitPlot <- renderPlot({
    req(input$fit_day)
    dirs <- fit_dirs()
    if (length(dirs) == 0) {
      plot.new()
      title(main = "Daily MLE Fit")
      text(0.5, 0.5, "No daily MLE fit data found")
      return(invisible(NULL))
    }
    fit_data <- withProgress(message = "Loading daily fit", detail = "Reading daily MLE fit data...", value = 0, {
      incProgress(0.4)
      out <- read_day_partition(dirs[1], input$fit_day)
      incProgress(0.6)
      out
    })
    if (is.null(fit_data) || nrow(fit_data) == 0) {
      plot.new()
      title(main = "Daily MLE Fit")
      text(0.5, 0.5, "No daily MLE fit data found")
      return(invisible(NULL))
    }

    plot_daily_mle_fit(
      fit_data,
      selected_day = input$fit_day,
      title = paste("Window", fit_window_length())
    )
  })

  day_plot_data <- reactive({
    req(input$fit_day)
    dirs <- day_dirs()
    if (length(dirs) == 0) return(NULL)
    out <- withProgress(message = "Loading daily columns", detail = "Reading day-level column data...", value = 0, {
      incProgress(0.4)
      out <- read_day_partition(dirs[1], input$fit_day)
      incProgress(0.6)
      out
    })
    if (is.null(out) || nrow(out) == 0) {
      return(NULL)
    }
    if ("timestamp" %in% names(out)) {
      out$timestamp <- as.POSIXct(out$timestamp, tz = "UTC")
    }
    out
  })

  output$dailyColumnsPlot <- renderPlot({
    day_data <- day_plot_data()
    plot_daily_columns(
      day_data,
      selected_day = input$fit_day,
      title = paste("Window", fit_window_length(), "-", as.character(as.Date(input$fit_day)))
    )
  })

  output$dailyAdjustedColumnPlot <- renderPlot({
    day_data <- day_plot_data()
    row <- daily_summary_row()
    intercept_value <- NA_real_
    if (!is.null(row) && "intercept" %in% names(row)) {
      intercept_value <- suppressWarnings(as.numeric(row$intercept[1]))
    }
    plot_daily_adjusted_column(
      day_data,
      intercept_value = intercept_value,
      selected_day = input$fit_day,
      title = paste("Adjusted VC", fit_window_length(), "-", as.character(as.Date(input$fit_day)))
    )
  })

  daily_difference_series <- reactive({
    req(dataset(), input$fit_window)
    summary_df <- dataset()[[input$fit_window]]
    if (is.null(summary_df) || nrow(summary_df) == 0) {
      return(NULL)
    }

    dirs <- day_dirs()
    withProgress(message = "Loading difference time series", detail = "Preparing daily raw/adjusted difference...", value = 0, {
      incProgress(0.3)
      out <- compute_daily_difference_series(summary_df, day_dir = dirs)
      incProgress(0.7)
      out
    })
  })

  output$dailyDifferenceSeriesPlot <- renderPlot({
    plot_daily_difference_series(
      daily_difference_series(),
      selected_day = input$fit_day,
      title = paste("Window", fit_window_length(), "- median daily difference")
    )
  })

  output$dailyDifferencePlot <- renderPlot({
    day_data <- day_plot_data()
    row <- daily_summary_row()
    intercept_value <- NA_real_
    if (!is.null(row) && "intercept" %in% names(row)) {
      intercept_value <- suppressWarnings(as.numeric(row$intercept[1]))
    }
    plot_daily_difference(
      day_data,
      intercept_value = intercept_value,
      selected_day = input$fit_day,
      title = paste("Difference to vc_trop_int", fit_window_length(), "-", as.character(as.Date(input$fit_day)))
    )
  })
  
  # --- Breakpoint analysis and plotting ---
  bp_data <- eventReactive(input$analyze_bp, {
    req(dataset(), input$bp_window, input$parameter)
    status_message(sprintf("Running breakpoint analysis for %s...", input$bp_window))
    showNotification("Running breakpoint analysis...", id = "bp_load", duration = NULL, closeButton = FALSE)
    on.exit(removeNotification("bp_load"), add = TRUE)

    df <- withProgress(message = "Breakpoint analysis", detail = "Preparing selected series...", value = 0, {
      incProgress(0.3)
      out <- dataset()[[input$bp_window]]
      incProgress(0.2)
      out
    })
    req(!is.null(df), nrow(df) > 0)

    df <- df %>% arrange(as.Date(date))
    y_all <- suppressWarnings(as.numeric(df[[input$parameter]]))
    x_all <- as.Date(df$date)
    keep <- is.finite(y_all) & !is.na(x_all)
    y_vec <- y_all[keep]
    x_vec <- x_all[keep]
    if (length(y_vec) < 8) {
      return(list(error = "Need at least 8 finite points for breakpoint analysis"))
    }
    y <- zoo(y_vec, x_vec)

    n <- length(y)
    window_days <- suppressWarnings(as.integer(sub("window_", "", input$bp_window)))
    if (is.na(window_days) || window_days < 1L) {
      return(list(error = "Could not derive window length for breakpoint minimum distance"))
    }
    h_raw <- window_days / n
    h_min <- (2 + 1e-8) / n
    h_use <- max(h_raw, h_min)

    # breakpoints() requires h < 0.5 for constant fit with one regressor.
    if (h_use >= 0.5) {
      return(list(
        y = y,
        n_breaks = 0L,
        break_dates = as.Date(character(0)),
        h_use = h_use,
        window_days = window_days,
        note = "Series too short for breaks under this minimum distance",
        error = NULL
      ))
    }

    bp_obj <- withProgress(message = "Breakpoint analysis", detail = "Fitting breakpoint model...", value = 0, {
      incProgress(0.2)
      out <- tryCatch(
        pandoniaR::makebreakpoints(
          x = y,
          fit = "constant",
          plot = FALSE,
          ylab = input$parameter,
          h = h_use
        ),
        error = function(e) {
          e
        }
      )
      incProgress(0.3)
      out
    })
    if (inherits(bp_obj, "error")) {
      return(list(error = paste("Breakpoint fitting failed:", conditionMessage(bp_obj))))
    }

    bp_idx <- bp_obj$breakpoints
    if (length(bp_idx) == 0 || all(is.na(bp_idx))) {
      bp_idx <- integer(0)
    } else {
      bp_idx <- bp_idx[is.finite(bp_idx)]
    }
    break_dates <- if (length(bp_idx) > 0) index(y)[bp_idx] else as.Date(character(0))
    n_breaks <- as.integer(length(bp_idx))
    segment_edges <- c(1L, bp_idx + 1L, length(y) + 1L)
    segment_starts <- index(y)[segment_edges[-length(segment_edges)]]
    segment_ends <- index(y)[segment_edges[-1L] - 1L]
    segment_means <- vapply(seq_len(length(segment_starts)), function(i) {
      seg_idx <- segment_edges[i]:(segment_edges[i + 1L] - 1L)
      mean(as.numeric(y[seg_idx]), na.rm = TRUE)
    }, numeric(1))

    status_message(sprintf("Breakpoint analysis ready for %s with %d break(s).", input$bp_window, n_breaks))
    list(
      y = y,
      break_dates = break_dates,
      n_breaks = n_breaks,
      segment_starts = segment_starts,
      segment_ends = segment_ends,
      segment_means = segment_means,
      h_use = h_use,
      window_days = window_days,
      note = NULL,
      error = NULL
    )
  }, ignoreInit = TRUE)

  output$bpPlot <- renderPlot({
    bp <- bp_data()
    req(bp)

    if (!is.null(bp$error)) {
      plot.new()
      title(main = "Breakpoints")
      text(0.5, 0.5, bp$error)
      return(invisible(NULL))
    }

    bp_uncertainty <- calibration_uncertainty()
    show_bp_uncertainty <- isTRUE(input$show_calibration_uncertainty) &&
      identical(input$parameter, "intercept") &&
      !is.null(bp_uncertainty) &&
      nrow(bp_uncertainty) > 0
    sigma_level <- uncertainty_sigma_level()
    bp_ylim <- NULL
    if (show_bp_uncertainty) {
      bp_ylim <- expand_ylim(
        range(as.numeric(bp$y), na.rm = TRUE),
        c(
          abs(as.numeric(bp_uncertainty$scref_uncert)) * sigma_level,
          -abs(as.numeric(bp_uncertainty$scref_uncert)) * sigma_level
        )
      )
    }

    plot(
      bp$y,
      type = "b",
      pch = 16,
      col = adjust_transparency("black", 0.4),
      ylim = bp_ylim,
      ylab = input$parameter,
      xlab = "Time [UTC]",
      main = sprintf(
        "Breakpoints (%s, h=%.4f, min distance=%d days, breaks=%d)",
        input$bp_window, bp$h_use, bp$window_days, bp$n_breaks
      )
    )
    grid()
    abline(h=0,lty=1,col="lightgrey",lwd=2)

    if (!is.null(bp$segment_starts) && length(bp$segment_starts) > 0) {
      for (i in seq_along(bp$segment_starts)) {
        segments(
          x0 = bp$segment_starts[i],
          y0 = bp$segment_means[i],
          x1 = bp$segment_ends[i],
          y1 = bp$segment_means[i],
          col = "blue",
          lwd = 3
        )
      }
    }

    if (!is.null(bp$break_dates) && length(bp$break_dates) > 0) {
      for (break_date in bp$break_dates) {
        abline(v = break_date, col = "blue", lwd = 2, lty = 3)
      }
    }

    if (show_bp_uncertainty) {
      period_drawn <- draw_synt_ref_periods(bp_uncertainty, zoo::index(bp$y))
      uncertainty_drawn <- draw_calibration_uncertainty_bars(bp_uncertainty, zoo::index(bp$y), sigma_level = sigma_level)
      if (uncertainty_drawn || period_drawn) {
        legend(
          "topright",
          legend = c(
            if (uncertainty_drawn) sprintf("Initial cal. uncertainty (%s-sigma)", as.character(sigma_level)),
            if (period_drawn) "Synthetic ref period"
          ),
          col = c(
            if (uncertainty_drawn) "#D55E00",
            if (period_drawn) NA_character_
          ),
          lty = c(
            if (uncertainty_drawn) 1,
            if (period_drawn) NA_real_
          ),
          lwd = c(
            if (uncertainty_drawn) 3,
            if (period_drawn) NA_real_
          ),
          fill = c(
            if (uncertainty_drawn) NA_character_,
            if (period_drawn) adjust_transparency("#D55E00", alpha = 0.12)
          ),
          border = NA,
          bty = "n"
        )
      }
    }

    if (!is.null(bp$note)) {
      mtext(bp$note, side = 3, line = 0.2, col = "grey35", cex = 0.9)
    }
  })
}

shinyApp(ui, server)