library(shiny)
library(bslib)
library(ggplot2)
library(dplyr)
library(scales)
library(moments)

plot_theme <- theme_bw() +
  theme(
    axis.text   = element_text(size = 13),
    axis.title  = element_text(size = 14),
    legend.text = element_text(size = 10)
  )

ui <- page_navbar(
  title = "Uncertainty - Sampling Distributions",
  theme = bs_theme(bootswatch = "flatly", font_scale = 0.9),
  nav_panel(
    title = "Population",
    layout_sidebar(
      sidebar = sidebar(
        title = "Population Parameters",
        numericInput("pop_mean", "Mean:", value = 3.6, step = 0.1),
        numericInput("pop_sd", "Standard Deviation:", value = 0.5, min = 0, step = 0.1),
        numericInput("pop_skew", "Skew:", value = 0.3, step = 0.1),
        actionButton("pop_run", "Run")
      ),
      layout_columns(
        card(
          card_header("Flood Frequency Summary Table"),
          card_body(tableOutput("pop_table")),
          full_screen = TRUE
        ),
        card(
          card_header("Population LP3 Frequency Curve"),
          card_body(plotOutput("pop_freq_curve_plot")),
          full_screen = TRUE
        ),
        card(
          card_header("PDF of Flow"),
          card_body(plotOutput("prob_flow_plot")),
          full_screen = TRUE
        ),
        card(
          card_header("PDF of Log Flow"),
          card_body(plotOutput("prob_log_flow_plot")),
          full_screen = TRUE
        ),
        col_widths = c(6, 6, 6, 6)
      )
    )
  ),
  nav_panel(
    title = "Single Experiment",
    layout_sidebar(
      sidebar = sidebar(
        title = "Population Parameters",
        numericInput("sing_mean", "Mean:", value = 3.6, step = 0.1),
        numericInput("sing_sd", "Standard Deviation:", value = 0.5, min = 0, step = 0.1),
        numericInput("sing_skew", "Skew:", value = 0.3, step = 0.1),
        hr(),
        numericInput("sing_sample1", "Sample 1 Size (10-1,000):", value = 10, min = 10, max = 1000),
        numericInput("sing_sample2", "Sample 2 Size (10-1,000):", value = 50, min = 10, max = 1000),
        actionButton("sing_run", "Run")
      ),
      layout_columns(
        card(
          card_header("Sample Curve Parameters"),
          card_body(
            layout_column_wrap(
              width = 1/2,
              card(card_header("Sample 1"), tableOutput("sing_samp1_table")),
              card(card_header("Sample 2"), tableOutput("sing_samp2_table"))
            )
          ),
          full_screen = TRUE
        ),
        col_widths = 12
      ),
      layout_columns(
        card(
          card_header("PDF of Log Flow"),
          card_body(plotOutput("single_prob_log_flow_plot")),
          full_screen = TRUE
        ),
        card(
          card_header("Population and Sample LP3 Frequency Curves"),
          card_body(plotOutput("single_freq_curve_plot")),
          full_screen = TRUE
        ),
        col_widths = c(6, 6)
      )
    )
  ),
  nav_panel(
    title = "Multi Experiment",
    layout_sidebar(
      sidebar = sidebar(
        title = "Population Parameters",
        numericInput("mult_mean", "Mean:", value = 3.6, step = 0.1),
        numericInput("mult_sd", "Standard Deviation:", value = 0.5, min = 0, step = 0.1),
        numericInput("mult_skew", "Skew:", value = 0.3, step = 0.1),
        hr(),
        numericInput("mult_sample1", "Sample 1 Size (10-1,000):", value = 10, min = 10, max = 1000),
        numericInput("mult_sample2", "Sample 2 Size (10-1,000):", value = 50, min = 10, max = 1000),
        numericInput("mult_exp", "Number of Experiments: (10-1,000):", value = 1000, min = 10, max = 1000),
        actionButton("mult_run", "Run")
      ),
      layout_columns(
        card(
          card_header("Sample Curve Parameters"),
          card_body(
            radioButtons(
              "mult_table_view", label = NULL,
              choices  = c("Average Across Experiments" = "avg",
                           "All Experiments"            = "all"),
              selected = "avg", inline = TRUE
            ),
            conditionalPanel(
              condition = "input.mult_table_view == 'avg'",
              layout_column_wrap(
                width = 1/2,
                card(card_header("Sample 1"), tableOutput("mult_samp1_table")),
                card(card_header("Sample 2"), tableOutput("mult_samp2_table"))
              )
            ),
            conditionalPanel(
              condition = "input.mult_table_view == 'all'",
              layout_column_wrap(
                width = 1/2,
                card(
                  card_header("Sample 1"),
                  div(style = "overflow-y: auto; max-height: 380px;",
                      tableOutput("mult_samp1_all_table"))
                ),
                card(
                  card_header("Sample 2"),
                  div(style = "overflow-y: auto; max-height: 380px;",
                      tableOutput("mult_samp2_all_table"))
                )
              )
            )
          ),
          full_screen = TRUE
        ),
        col_widths = 12
      ),
      layout_columns(
        navset_card_underline(
          nav_panel(
            title = "Sampling Distribution of Parameters",
            card_body(
              layout_column_wrap(
                width = 1/3,
                card(
                  card_header("Mean"),
                  plotOutput("mult_mean_plot"),
                  uiOutput("mult_mean_notes"),
                  full_screen = TRUE
                ),
                card(
                  card_header("Standard Deviation"),
                  plotOutput("mult_sd_plot"),
                  uiOutput("mult_sd_notes"),
                  full_screen = TRUE
                ),
                card(
                  card_header("Skew"),
                  plotOutput("mult_skew_plot"),
                  uiOutput("mult_skew_notes"),
                  full_screen = TRUE
                )
              )
            )
          ),
          nav_panel(
            title = "Sampling Distribution of Q",
            card_body(
              layout_column_wrap(
                width = 1/2,
                selectizeInput("mult_flow_quantile", "Flow Quantile (Yrs):",
                               c(1.01, 1.11, 1.33, 2, 4, 10, 20, 50, 100, 200, 500), 100),
                div(style = "margin-top: 32px;",
                    actionButton("mult_quantile_run", "Generate Plots",
                                 class = "btn-secondary w-100"))
              ),
              h5(),
              layout_column_wrap(
                width = 1/3,
                card(card_header("Log Q"),  plotOutput("mult_logq_plot"), full_screen = TRUE),
                card(card_header("p(Q)"),   plotOutput("mult_pq_plot"),   full_screen = TRUE),
                card(card_header("Zp(Q)"),  plotOutput("mult_zpq_plot"),  full_screen = TRUE)
              )
            )
          ),
          nav_panel(
            title = "Sampled Curves",
            card_body(
              layout_column_wrap(
                width = 1/2,
                card(card_header("Sample 1"), plotOutput("mult_freq_samp1_plot"), full_screen = TRUE),
                card(card_header("Sample 2"), plotOutput("mult_freq_samp2_plot"), full_screen = TRUE)
              )
            )
          )
        )
      )
    )
  )
)

server <- function(input, output) {
  
  ##############################################################################
  ### population
  list_df_pop <- eventReactive(input$pop_run, {
    
    pop_mean <- input$pop_mean
    pop_sd   <- input$pop_sd
    pop_skew <- input$pop_skew
    
    calculate_stats <- function(exc_prob, mean_val, sd_val, skew_val) {
      data.frame(exc_prob = exc_prob) %>%
        mutate(
          Zp   = qnorm(1 - exc_prob),
          kp_g = 2 / skew_val * (((Zp - skew_val / 6) * skew_val / 6 + 1) ^ 3 - 1),
          flow = 10 ^ (mean_val + sd_val * kp_g)
        )
    }
    
    exc_prob_small <- c(0.99, 0.9, 0.75, 0.5, 0.25, 0.10, 0.05, 0.02, 0.01, 0.005, 0.002)
    df_pop_small <- calculate_stats(exc_prob_small, pop_mean, pop_sd, pop_skew)
    
    exc_prob_large <- c(0.9995, 0.999, 0.998, seq(0.995, 0.005, by = -0.005), 0.002, 0.001)
    df_pop_large <- calculate_stats(exc_prob_large, pop_mean, pop_sd, pop_skew) %>%
      mutate(
        log_flow               = log10(flow),
        flow_prob_per_unit     = -(lead(exc_prob) - exc_prob) / (lead(flow)     - flow),
        log_flow_prob_per_unit = -(lead(exc_prob) - exc_prob) / (lead(log_flow) - log_flow)
      )
    
    list(df_pop_small = df_pop_small, df_pop_large = df_pop_large)
  })
  
  output$pop_table <- renderTable({ list_df_pop()$df_pop_small })
  
  output$pop_freq_curve_plot <- renderPlot({
    df <- list_df_pop()$df_pop_small
    ggplot(df, aes(x = Zp, y = flow)) +
      geom_line(color = "darkgreen", linewidth = 1) +
      scale_y_log10(breaks = c(100, 1000, 10000, 100000, 1000000), labels = label_comma()) +
      scale_x_continuous(breaks = df$Zp, labels = as.character(df$exc_prob * 100)) +
      annotation_logticks(sides = "l") +
      labs(x = "Percent Chance Exceedance", y = "Annual Peak Streamflow (cfs)") +
      plot_theme
  })
  
  output$prob_flow_plot <- renderPlot({
    df <- list_df_pop()$df_pop_large
    ggplot(df, aes(x = flow, y = flow_prob_per_unit)) +
      geom_line(color = "darkgreen", linewidth = 1) +
      scale_x_log10(breaks = c(100, 1000, 10000, 100000), labels = label_comma()) +
      labs(x = "Flow (cfs)", y = "Prob per Unit") +
      plot_theme
  })
  
  output$prob_log_flow_plot <- renderPlot({
    df <- list_df_pop()$df_pop_large
    ggplot(df, aes(x = log_flow, y = log_flow_prob_per_unit)) +
      geom_line(color = "darkgreen", linewidth = 1) +
      labs(x = "Log Flow (cfs)", y = "Prob per Unit") +
      plot_theme
  })
  
  ##############################################################################
  ### single experiment
  list_single <- eventReactive(input$sing_run, {
    
    sing_mean <- input$sing_mean
    sing_sd   <- input$sing_sd
    sing_skew <- input$sing_skew
    
    calculate_sample_stats <- function(rand, mean_val, sd_val, skew_val) {
      data.frame(rand = rand) %>%
        mutate(
          Kg       = 2 / skew_val * (((qnorm(rand) - skew_val / 6) * skew_val / 6 + 1) ^ 3 - 1),
          log_flow = mean_val + sd_val * Kg,
          flow     = 10 ^ log_flow
        )
    }
    
    sing_sample1_size <- input$sing_sample1
    df_sample1 <- calculate_sample_stats(runif(sing_sample1_size), sing_mean, sing_sd, sing_skew) %>%
      mutate(prob_per_unit = 0.1)
    df_sample1_summary <- df_sample1 %>%
      mutate(mean = mean(log_flow), sd = sd(log_flow), skew = skewness(log_flow)) %>%
      distinct(mean, sd, skew)
    
    sing_sample2_size <- input$sing_sample2
    df_sample2 <- calculate_sample_stats(runif(sing_sample2_size), sing_mean, sing_sd, sing_skew) %>%
      mutate(prob_per_unit = 0.2)
    df_sample2_summary <- df_sample2 %>%
      mutate(mean = mean(log_flow), sd = sd(log_flow), skew = skewness(log_flow)) %>%
      distinct(mean, sd, skew)
    
    calculate_stats <- function(exc_prob, mean_val, sd_val, skew_val) {
      data.frame(exc_prob = exc_prob) %>%
        mutate(
          Zp   = qnorm(1 - exc_prob),
          kp_g = 2 / skew_val * (((Zp - skew_val / 6) * skew_val / 6 + 1) ^ 3 - 1),
          flow = 10 ^ (mean_val + sd_val * kp_g)
        )
    }
    
    exc_prob <- c(0.99, 0.95, seq(0.90, 0.10, by = -0.1), 0.05, 0.02, 0.01, 0.005, 0.002)
    
    df_sample1_lp3 <- calculate_stats(exc_prob, df_sample1_summary$mean,
                                      df_sample1_summary$sd, df_sample1_summary$skew)
    df_sample2_lp3 <- calculate_stats(exc_prob, df_sample2_summary$mean,
                                      df_sample2_summary$sd, df_sample2_summary$skew)
    
    exc_prob_large <- c(0.9995, 0.999, 0.998, seq(0.995, 0.005, by = -0.005), 0.002, 0.001)
    df_single_actual_lp3 <- calculate_stats(exc_prob_large, sing_mean, sing_sd, sing_skew) %>%
      mutate(
        log_flow               = log10(flow),
        flow_prob_per_unit     = -(lead(exc_prob) - exc_prob) / (lead(flow)     - flow),
        log_flow_prob_per_unit = -(lead(exc_prob) - exc_prob) / (lead(log_flow) - log_flow)
      )
    
    exc_prob_small <- c(0.99, 0.9, 0.75, 0.5, 0.25, 0.10, 0.05, 0.02, 0.01, 0.005, 0.002)
    df_single_actual_lp3_small <- calculate_stats(exc_prob_small, sing_mean, sing_sd, sing_skew)
    
    calculate_plot_pos <- function(n, df) {
      flow <- sort(df$flow, decreasing = TRUE)
      data.frame(sample_size_id = 1:n) %>%
        mutate(col2 = sample_size_id / (n + 1), Zp = qnorm(1 - col2), flow = flow[sample_size_id])
    }
    
    list(
      sing_mean                  = sing_mean,
      sing_sd                    = sing_sd,
      sing_skew                  = sing_skew,
      sing_sample1_size          = sing_sample1_size,
      sing_sample2_size          = sing_sample2_size,
      df_sample1                 = df_sample1,
      df_sample2                 = df_sample2,
      df_sample1_summary         = df_sample1_summary,
      df_sample2_summary         = df_sample2_summary,
      df_sample1_lp3             = df_sample1_lp3,
      df_sample2_lp3             = df_sample2_lp3,
      df_single_actual_lp3       = df_single_actual_lp3,
      df_single_actual_lp3_small = df_single_actual_lp3_small,
      df_sample1_plot_pos        = calculate_plot_pos(sing_sample1_size, df_sample1),
      df_sample2_plot_pos        = calculate_plot_pos(sing_sample2_size, df_sample2)
    )
  })
  
  output$sing_samp1_table <- renderTable({ list_single()$df_sample1_summary })
  output$sing_samp2_table <- renderTable({ list_single()$df_sample2_summary })
  
  output$single_prob_log_flow_plot <- renderPlot({
    
    df_single_actual_lp3 <- list_single()$df_single_actual_lp3
    df_sample1           <- list_single()$df_sample1
    df_sample2           <- list_single()$df_sample2
    sing_mean            <- list_single()$sing_mean
    sing_sample1_size    <- list_single()$sing_sample1_size
    sing_sample2_size    <- list_single()$sing_sample2_size
    mean_sample1         <- list_single()$df_sample1_summary$mean
    mean_sample2         <- list_single()$df_sample2_summary$mean
    
    leg_samp1_pt <- paste0("sample N=", sing_sample1_size)
    leg_samp2_pt <- paste0("sample N=", sing_sample2_size)
    leg_samp1_ln <- paste0("n = ", sing_sample1_size)
    leg_samp2_ln <- paste0("n = ", sing_sample2_size)
    
    manual_values <- setNames(
      c("darkgreen", "blue", "blue", "red", "red"),
      c("mean", leg_samp1_ln, leg_samp1_pt, leg_samp2_ln, leg_samp2_pt)
    )
    
    ggplot(df_single_actual_lp3, aes(x = log_flow, y = log_flow_prob_per_unit)) +
      geom_line(color = "darkgreen", linewidth = 1) +
      geom_point(data = df_sample1,
                 aes(x = log_flow, y = prob_per_unit, color = leg_samp1_pt),
                 size = 2.5, shape = 18) +
      geom_point(data = df_sample2,
                 aes(x = log_flow, y = prob_per_unit, color = leg_samp2_pt),
                 size = 2.5, shape = 18) +
      geom_vline(aes(xintercept = sing_mean,    color = "mean"),
                 linewidth = 1, key_glyph = "path") +
      geom_vline(aes(xintercept = mean_sample1, color = leg_samp1_ln),
                 linetype = "dashed", linewidth = 1, key_glyph = "path") +
      geom_vline(aes(xintercept = mean_sample2, color = leg_samp2_ln),
                 linetype = "dashed", linewidth = 1, key_glyph = "path") +
      scale_color_manual(
        name   = "Legend",
        values = manual_values,
        breaks = c("mean", leg_samp1_ln, leg_samp1_pt, leg_samp2_ln, leg_samp2_pt)
      ) +
      guides(color = guide_legend(override.aes = list(
        shape    = c(NA, NA, 18, NA, 18),
        linetype = c("solid", "dashed", "blank", "dashed", "blank")
      ))) +
      labs(x = "Log Flow (cfs)", y = "Prob per Unit") +
      plot_theme +
      theme(legend.key.width = unit(1.5, "cm"))
  })
  
  make_param_note <- function(s1_vals, s2_vals, n1, n2) {
    fmt <- function(vals) {
      paste0("Mean = ",  round(mean(vals),         2),
             ",   Stdev = ", round(sd(vals),        2),
             ",   Skew = ",  round(skewness(vals),  2))
    }
    tags$div(
      style = "margin-top: 2px; font-size: 0.88em; padding: 2px 8px;",
      tags$p(style = "margin: 0;", tags$strong(paste0("Sample 1:  ")), fmt(s1_vals)),
      tags$p(style = "margin: 0;", tags$strong(paste0("Sample 2:  ")), fmt(s2_vals))
    )
  }
  
  output$mult_mean_notes <- renderUI({
    req(list_mult())
    make_param_note(
      s1_vals = list_mult()$df_exp_samp1_dis$mean_log_flow,
      s2_vals = list_mult()$df_exp_samp2_dis$mean_log_flow,
      n1 = list_mult()$mult_sample1_size,
      n2 = list_mult()$mult_sample2_size
    )
  })
  
  output$mult_sd_notes <- renderUI({
    req(list_mult())
    make_param_note(
      s1_vals = list_mult()$df_exp_samp1_dis$std_log_flow,
      s2_vals = list_mult()$df_exp_samp2_dis$std_log_flow,
      n1 = list_mult()$mult_sample1_size,
      n2 = list_mult()$mult_sample2_size
    )
  })
  
  output$mult_skew_notes <- renderUI({
    req(list_mult())
    make_param_note(
      s1_vals = list_mult()$df_exp_samp1_dis$skew_log_flow,
      s2_vals = list_mult()$df_exp_samp2_dis$skew_log_flow,
      n1 = list_mult()$mult_sample1_size,
      n2 = list_mult()$mult_sample2_size
    )
  })
  
  output$single_freq_curve_plot <- renderPlot({
    
    df_single_actual_lp3_small <- list_single()$df_single_actual_lp3_small
    df_sample1_lp3             <- list_single()$df_sample1_lp3
    df_sample2_lp3             <- list_single()$df_sample2_lp3
    df_sample1_plot_pos        <- list_single()$df_sample1_plot_pos
    df_sample2_plot_pos        <- list_single()$df_sample2_plot_pos
    sing_sample1_size          <- list_single()$sing_sample1_size
    sing_sample2_size          <- list_single()$sing_sample2_size
    
    leg_samp1_pt <- paste0("sample N=", sing_sample1_size)
    leg_samp2_pt <- paste0("sample N=", sing_sample2_size)
    leg_samp1_ln <- paste0("n = ", sing_sample1_size)
    leg_samp2_ln <- paste0("n = ", sing_sample2_size)
    
    manual_values <- setNames(
      c("darkgreen", "blue", "blue", "red", "red"),
      c("population", leg_samp1_ln, leg_samp1_pt, leg_samp2_ln, leg_samp2_pt)
    )
    
    ggplot(df_single_actual_lp3_small, aes(x = Zp, y = flow)) +
      geom_line(aes(color = "population"), linewidth = 1) +
      scale_y_log10(breaks = c(100, 1000, 10000, 100000, 1000000), labels = label_comma()) +
      scale_x_continuous(breaks = df_single_actual_lp3_small$Zp,
                         labels = as.character(df_single_actual_lp3_small$exc_prob * 100)) +
      geom_line(data = df_sample1_lp3,
                aes(x = Zp, y = flow, color = leg_samp1_ln), linewidth = 1) +
      geom_line(data = df_sample2_lp3,
                aes(x = Zp, y = flow, color = leg_samp2_ln), linewidth = 1) +
      geom_point(data = df_sample1_plot_pos,
                 aes(x = Zp, y = flow, color = leg_samp1_pt), shape = 18) +
      geom_point(data = df_sample2_plot_pos,
                 aes(x = Zp, y = flow, color = leg_samp2_pt), shape = 18) +
      scale_color_manual(
        name   = "Legend",
        values = manual_values,
        breaks = c("population", leg_samp1_ln, leg_samp1_pt, leg_samp2_ln, leg_samp2_pt)
      ) +
      guides(color = guide_legend(override.aes = list(
        shape    = c(NA, NA, 18, NA, 18),
        linetype = c("solid", "solid", "blank", "solid", "blank")
      ))) +
      annotation_logticks(sides = "l") +
      labs(x = "Percent Chance Exceedance", y = "Annual Peak Streamflow (cfs)") +
      plot_theme
  })
  
  ##############################################################################
  ### multiple experiment
  
  q_plot_data  <- reactiveVal(NULL)
  pq_plot_data <- reactiveVal(NULL)
  
  list_mult <- eventReactive(input$mult_run, {
    
    q_plot_data(NULL)
    pq_plot_data(NULL)
    
    mult_mean <- input$mult_mean
    mult_sd   <- input$mult_sd
    mult_skew <- input$mult_skew
    exp_no    <- input$mult_exp
    
    calculate_stats <- function(rand, mean_val, sd_val, skew_val) {
      data.frame(rand = rand) %>%
        mutate(
          log_flow      = mean_val + sd_val * 2 / skew_val *
            (((qnorm(rand) - skew_val / 6) * skew_val / 6 + 1) ^ 3 - 1),
          mean_log_flow = mean(log_flow),
          std_log_flow  = sd(log_flow),
          skew_log_flow = skewness(log_flow)
        )
    }
    
    calculate_lp3 <- function(exc_prob, mean_val, sd_val, skew_val) {
      Zp   <- qnorm(1 - exc_prob)
      kp_g <- 2 / skew_val * (((Zp - skew_val / 6) * skew_val / 6 + 1) ^ 3 - 1)
      data.frame(exc_prob = exc_prob, Zp = Zp, flow = 10 ^ (mean_val + sd_val * kp_g))
    }
    
    compute_exp_curves <- function(df_dis, max_curves = 250) {
      exc_prob_pts <- c(0.99, 0.9, 0.75, 0.5, 0.25, 0.10, 0.05, 0.02, 0.01, 0.005, 0.002)
      n      <- min(nrow(df_dis), max_curves)
      df_sub <- df_dis[seq_len(n), ]
      do.call(rbind, lapply(seq_len(n), function(i) {
        d        <- calculate_lp3(exc_prob_pts,
                                  df_sub$mean_log_flow[i],
                                  df_sub$std_log_flow[i],
                                  df_sub$skew_log_flow[i])
        d$exp_no <- df_sub$exp_no[i]
        d
      }))
    }
    
    run_experiments <- function(sample_size, exp_no) {
      list_df <- vector("list", exp_no)
      for (i in seq_len(exp_no)) {
        df        <- calculate_stats(runif(sample_size), mult_mean, mult_sd, mult_skew)
        df$exp_no <- i
        list_df[[i]] <- df
      }
      bind_rows(list_df)
    }
    
    mult_sample1_size <- input$mult_sample1
    mult_sample2_size <- input$mult_sample2
    
    df_exp_samp1 <- run_experiments(mult_sample1_size, exp_no)
    df_exp_samp2 <- run_experiments(mult_sample2_size, exp_no)
    
    summarise_dis <- function(df, n) {
      df %>%
        distinct(mean_log_flow, std_log_flow, skew_log_flow, exp_no) %>%
        mutate(sample_size = n)
    }
    
    summarise_avg <- function(df) {
      df %>%
        mutate(
          mean_mean_log_flow = mean(mean_log_flow),
          mean_std_log_flow  = mean(std_log_flow),
          mean_skew_log_flow = mean(skew_log_flow)
        ) %>%
        distinct(mean_mean_log_flow, mean_std_log_flow, mean_skew_log_flow)
    }
    
    df_exp_samp1_dis <- summarise_dis(df_exp_samp1, mult_sample1_size)
    df_exp_samp2_dis <- summarise_dis(df_exp_samp2, mult_sample2_size)
    
    exc_prob_pts        <- c(0.99, 0.9, 0.75, 0.5, 0.25, 0.10, 0.05, 0.02, 0.01, 0.005, 0.002)
    df_pop_curve        <- calculate_lp3(exc_prob_pts, mult_mean, mult_sd, mult_skew)
    df_samp1_exp_curves <- compute_exp_curves(df_exp_samp1_dis, max_curves = 250)
    df_samp2_exp_curves <- compute_exp_curves(df_exp_samp2_dis, max_curves = 250)
    
    list(
      df_exp_samp1             = df_exp_samp1,
      df_exp_samp2             = df_exp_samp2,
      df_exp_samp1_dis         = df_exp_samp1_dis,
      df_exp_samp2_dis         = df_exp_samp2_dis,
      df_exp_dis               = bind_rows(df_exp_samp1_dis, df_exp_samp2_dis) %>%
        mutate(sample_size = as.factor(sample_size)),
      df_exp_samp1_avg_summary = summarise_avg(df_exp_samp1_dis),
      df_exp_samp2_avg_summary = summarise_avg(df_exp_samp2_dis),
      mult_sample1_size        = mult_sample1_size,
      mult_sample2_size        = mult_sample2_size,
      # Population parameters stored so observeEvent can use the frozen values
      mult_mean               = mult_mean,
      mult_sd                 = mult_sd,
      mult_skew               = mult_skew,
      df_pop_curve            = df_pop_curve,
      df_samp1_exp_curves     = df_samp1_exp_curves,
      df_samp2_exp_curves     = df_samp2_exp_curves
    )
  })
  
  # ── Average tables ───────────────────────────────────────────────────────────
  output$mult_samp1_table <- renderTable({
    list_mult()$df_exp_samp1_avg_summary %>%
      rename("Mean of Log Flow"  = mean_mean_log_flow,
             "Stdev of Log Flow" = mean_std_log_flow,
             "Skew Coefficient"  = mean_skew_log_flow)
  })
  
  output$mult_samp2_table <- renderTable({
    list_mult()$df_exp_samp2_avg_summary %>%
      rename("Mean of Log Flow"  = mean_mean_log_flow,
             "Stdev of Log Flow" = mean_std_log_flow,
             "Skew Coefficient"  = mean_skew_log_flow)
  })
  
  # ── Per-experiment tables ────────────────────────────────────────────────────
  output$mult_samp1_all_table <- renderTable({
    list_mult()$df_exp_samp1_dis %>%
      arrange(exp_no) %>%
      select(exp_no, mean_log_flow, std_log_flow, skew_log_flow) %>%
      rename("Experiment"        = exp_no,
             "Mean of Log Flow"  = mean_log_flow,
             "Stdev of Log Flow" = std_log_flow,
             "Skew Coefficient"  = skew_log_flow)
  })
  
  output$mult_samp2_all_table <- renderTable({
    list_mult()$df_exp_samp2_dis %>%
      arrange(exp_no) %>%
      select(exp_no, mean_log_flow, std_log_flow, skew_log_flow) %>%
      rename("Experiment"        = exp_no,
             "Mean of Log Flow"  = mean_log_flow,
             "Stdev of Log Flow" = std_log_flow,
             "Skew Coefficient"  = skew_log_flow)
  })
  
  # ── Histogram plots ──────────────────────────────────────────────────────────
  output$mult_mean_plot <- renderPlot({
    ggplot(list_mult()$df_exp_dis) +
      geom_histogram(aes(x = mean_log_flow, fill = sample_size), position = "dodge") +
      labs(x = "Mean of Log Flow") +
      plot_theme
  })
  
  output$mult_sd_plot <- renderPlot({
    ggplot(list_mult()$df_exp_dis) +
      geom_histogram(aes(x = std_log_flow, fill = sample_size), position = "dodge") +
      labs(x = "Stdev of Log Flow") +
      plot_theme
  })
  
  output$mult_skew_plot <- renderPlot({
    ggplot(list_mult()$df_exp_dis) +
      geom_histogram(aes(x = skew_log_flow, fill = sample_size), position = "dodge") +
      labs(x = "Skew Coefficient") +
      plot_theme
  })
  
  # ── Sampled Curves plots ─────────────────────────────────────────────────────
  make_sampled_curves_plot <- function(df_pop_curve, df_exp_curves, sample_size) {
    n_curves <- n_distinct(df_exp_curves$exp_no)
    ggplot() +
      geom_line(data = df_exp_curves,
                aes(x = Zp, y = flow, group = exp_no),
                color = "grey60", linewidth = 0.3, alpha = 0.7) +
      geom_line(data = df_pop_curve,
                aes(x = Zp, y = flow),
                color = "darkgreen", linewidth = 1) +
      scale_y_log10(breaks = c(100, 1000, 10000, 100000, 1000000), labels = label_comma()) +
      scale_x_continuous(breaks = df_pop_curve$Zp,
                         labels = as.character(df_pop_curve$exc_prob * 100)) +
      annotation_logticks(sides = "l") +
      labs(x        = "Percent Chance Exceedance",
           y        = "Annual Peak Streamflow (cfs)",
           subtitle = paste0("Curves drawn: ", n_curves,
                             "  \u2022  Sample size (n): ", sample_size)) +
      plot_theme +
      theme(plot.subtitle = element_text(size = 12))
  }
  
  output$mult_freq_samp1_plot <- renderPlot({
    make_sampled_curves_plot(list_mult()$df_pop_curve,
                             list_mult()$df_samp1_exp_curves,
                             list_mult()$mult_sample1_size)
  })
  
  output$mult_freq_samp2_plot <- renderPlot({
    make_sampled_curves_plot(list_mult()$df_pop_curve,
                             list_mult()$df_samp2_exp_curves,
                             list_mult()$mult_sample2_size)
  })
  
  ##############################################################################
  ### Sampling Distribution of Q
  
  observeEvent(input$mult_quantile_run, {
    
    req(list_mult()$df_exp_dis)
    
    mult_flow_quantile <- as.numeric(input$mult_flow_quantile)
    flow_quantile_fix  <- 1 / mult_flow_quantile   # exceedance probability
    
    # Population parameters frozen at the time "Run" was clicked
    mult_mean  <- list_mult()$mult_mean
    mult_sd    <- list_mult()$mult_sd
    mult_skew  <- list_mult()$mult_skew
    
    # ── True (population) log Q at the chosen return period ──────────────────
    Zp_true    <- qnorm(1 - flow_quantile_fix)
    kp_g_true  <- 2 / mult_skew *
      (((Zp_true - mult_skew / 6) * mult_skew / 6 + 1) ^ 3 - 1)
    true_log_q <- mult_mean + mult_sd * kp_g_true
    
    # ── Per-experiment estimates ──────────────────────────────────────────────
    # safe_cbrt handles negative radicands that arise with negative skew
    safe_cbrt <- function(x) sign(x) * abs(x) ^ (1/3)
    
    df_q <- list_mult()$df_exp_dis %>%
      mutate(
        # Estimated log Q: apply each experiment's fitted LP3 at the chosen AEP
        log_q = mean_log_flow + std_log_flow * 2 / skew_log_flow *
          (((qnorm(1 - flow_quantile_fix) - skew_log_flow / 6) * skew_log_flow / 6 + 1) ^ 3 - 1),
        
        # Invert each experiment's LP3 to find what probability it assigns
        # to the TRUE population Q value
        #   Step 1: kp_g from the true log Q and this experiment's parameters
        kp_g_inv = (true_log_q - mean_log_flow) / std_log_flow,
        #   Step 2: analytically invert kp_g -> Zp (Wilson-Hilferty inverse)
        inner    = kp_g_inv * skew_log_flow / 2 + 1,
        Zp_q     = (safe_cbrt(inner) - 1) * (6 / skew_log_flow) + skew_log_flow / 6,
        #   Step 3: convert Zp to exceedance probability
        p_q      = 1 - pnorm(Zp_q)
      ) %>%
      select(-kp_g_inv, -inner)   # drop working columns
    
    # Store everything the three render functions need
    q_plot_data(list(
      df             = df_q,
      true_log_q     = true_log_q,
      true_p         = flow_quantile_fix,
      Zp_true        = Zp_true,
      quantile_label = mult_flow_quantile
    ))
  })
  
  # ── Log Q histogram ──────────────────────────────────────────────────────────
  output$mult_logq_plot <- renderPlot({
    req(q_plot_data())
    qd <- q_plot_data()
    
    ggplot(qd$df) +
      geom_histogram(aes(x = log_q, fill = sample_size), position = "dodge") +
      geom_vline(xintercept = qd$true_log_q,
                 color = "darkgreen", linewidth = 1, linetype = "dashed") +
      labs(
        x        = paste0("log Q (", qd$quantile_label, "-yr)"),
        subtitle = paste0("True log Q", qd$quantile_label,
                          " (population): ", round(qd$true_log_q, 4))
      ) +
      plot_theme +
      theme(plot.subtitle = element_text(size = 12))
  })
  
  # ── p(Q) histogram ───────────────────────────────────────────────────────────
  # Each bar is the exceedance probability that an experiment's fitted LP3
  # assigns to the true population Q value.  True answer = 1 / return period.
  output$mult_pq_plot <- renderPlot({
    req(q_plot_data())
    qd <- q_plot_data()
    
    ggplot(qd$df) +
      geom_histogram(aes(x = p_q, fill = sample_size), position = "dodge") +
      geom_vline(xintercept = qd$true_p,
                 color = "darkgreen", linewidth = 1, linetype = "dashed") +
      labs(
        x        = paste0("p(Q", qd$quantile_label, ")"),
        subtitle = paste0("True p(Q", qd$quantile_label, "): ", qd$true_p)
      ) +
      plot_theme +
      theme(plot.subtitle = element_text(size = 12))
  })
  
  # ── Zp(Q) histogram ──────────────────────────────────────────────────────────
  # Same as p(Q) but expressed as the standard normal deviate Zp.
  # True answer = qnorm(1 - 1/return period).
  output$mult_zpq_plot <- renderPlot({
    req(q_plot_data())
    qd <- q_plot_data()
    
    ggplot(qd$df) +
      geom_histogram(aes(x = Zp_q, fill = sample_size), position = "dodge") +
      geom_vline(xintercept = qd$Zp_true,
                 color = "darkgreen", linewidth = 1, linetype = "dashed") +
      labs(
        x        = paste0("Zp(Q", qd$quantile_label, ")"),
        subtitle = paste0("True Zp(Q", qd$quantile_label,
                          "): ", round(qd$Zp_true, 4))
      ) +
      plot_theme +
      theme(plot.subtitle = element_text(size = 12))
  })
}

shinyApp(ui = ui, server = server)