use("shiny")
use("bslib")
use("mrgsolve")
use("dplyr")
use("ggplot2")
use("tidyr")
use("tibble")

mw_cyclophosphamide_g_mol <- 261.086

paper_regimen <- function(key) {
  switch(
    key,
    fixed_500 = list(
      label = "Fixed 500 mg",
      dose_mode = "mg",
      dose_value = 500,
      note = "Euro-Lupus style fixed-dose reference used in the paper's Figure 4 simulation grid."
    ),
    fixed_1000 = list(
      label = "Fixed 1000 mg",
      dose_mode = "mg",
      dose_value = 1000,
      note = "Upper fixed-dose reference used in the paper's Figure 4 simulation grid."
    ),
    nih_500 = list(
      label = "NIH 500 mg/m2",
      dose_mode = "mg_m2",
      dose_value = 500,
      note = "Lower NIH body-surface-area regimen used in Figure 4."
    ),
    nih_1000 = list(
      label = "NIH 1000 mg/m2",
      dose_mode = "mg_m2",
      dose_value = 1000,
      note = "Upper NIH body-surface-area regimen used in Figure 4."
    )
  )
}

dose_to_umol <- function(dose_mode, dose_value, bsa_m2) {
  total_mg <- if (identical(dose_mode, "mg_m2")) {
    dose_value * bsa_m2
  } else {
    dose_value
  }

  total_mg / mw_cyclophosphamide_g_mol * 1000
}

total_dose_mg <- function(dose_mode, dose_value, bsa_m2) {
  if (identical(dose_mode, "mg_m2")) {
    dose_value * bsa_m2
  } else {
    dose_value
  }
}

trap_auc <- function(time, conc) {
  if (length(time) < 2 || length(conc) < 2) {
    return(NA_real_)
  }

  dt <- diff(time)
  avg <- (head(conc, -1) + tail(conc, -1)) / 2
  sum(dt * avg)
}

estimate_terminal_half_life <- function(time, conc, start_h = 4, end_h = 24) {
  keep <- is.finite(conc) & conc > 0 & time >= start_h & time <= end_h
  if (sum(keep) < 3) {
    return(NA_real_)
  }

  fit <- stats::lm(log(conc[keep]) ~ time[keep])
  slope <- unname(stats::coef(fit)[2])

  if (!is.finite(slope) || slope >= 0) {
    return(NA_real_)
  }

  log(2) / abs(slope)
}

build_events <- function(dose_mode, dose_value, bsa_m2, interval_days, n_doses, infusion_h = 1) {
  dose_umol <- dose_to_umol(dose_mode, dose_value, bsa_m2)
  dose_times_h <- seq(0, by = interval_days * 24, length.out = n_doses)

  tibble::tibble(
    ID = 1,
    time = dose_times_h,
    amt = dose_umol,
    cmt = 1,
    evid = 1,
    rate = dose_umol / infusion_h
  )
}

model_code <- "
$PARAM @annotated
CL    : 13.3  : Parent cyclophosphamide clearance (L/h)
VC    : 59.82 : Parent central volume for male reference patient (L)
K12   : 8.97  : Parent central-to-peripheral distribution rate (1/h)
K21   : 3.49  : Parent peripheral-to-central distribution rate (1/h)
KE    : 4.3   : 4OH-cyclophosphamide elimination rate (1/h)
PMET  : 0.64  : Fraction of parent clearance forming 4OH metabolite
SEXF  : 1     : Female indicator (0 = male, 1 = female)
TAU   : 672   : Dosing interval (h)

$CMT @annotated
CENT    : Cyclophosphamide central amount (umol)
PERIPH  : Cyclophosphamide peripheral amount (umol)
METAB   : 4OH-cyclophosphamide amount (umol)

$MAIN
double vc_i = VC * pow(0.75, SEXF);
double cl_rate = CL / vc_i;

$ODE
double cp = CENT / vc_i;
dxdt_CENT = -(cl_rate + K12) * CENT + K21 * PERIPH;
dxdt_PERIPH = K12 * CENT - K21 * PERIPH;
dxdt_METAB = PMET * cl_rate * CENT - KE * METAB;

$TABLE
double VC_I = vc_i;
double CY = CENT / vc_i;
double OH4 = METAB / vc_i;

$CAPTURE @annotated
CY   : Cyclophosphamide concentration (uM)
OH4  : 4OH-cyclophosphamide concentration (uM)
VC_I : Individual central volume (L)
"

mod <- mcode("cyclophosphamide_glomerulonephritis_353", model_code, quiet = TRUE)

simulate_profile <- function(
  dose_mode,
  dose_value,
  sex,
  bsa_m2,
  interval_days,
  n_doses,
  infusion_h,
  horizon_h
) {
  events <- build_events(dose_mode, dose_value, bsa_m2, interval_days, n_doses, infusion_h)
  last_dose_time <- max(events$time)
  sim_end <- last_dose_time + horizon_h

  out <- mod %>%
    param(
      SEXF = if (identical(sex, "Female")) 1 else 0,
      TAU = interval_days * 24
    ) %>%
    data_set(events) %>%
    mrgsim(end = sim_end, delta = 0.05) %>%
    as.data.frame() |>
    mutate(
      time_since_last_dose = time - last_dose_time,
      cycle_day = time / 24,
      cy_uM = CY,
      oh4_uM = OH4
    )

  long <- out |>
    filter(time_since_last_dose >= 0, time_since_last_dose <= horizon_h) |>
    select(time_since_last_dose, cy_uM, oh4_uM) |>
    pivot_longer(
      cols = c(cy_uM, oh4_uM),
      names_to = "analyte",
      values_to = "conc_uM"
    ) |>
    mutate(
      analyte = recode(
        analyte,
        cy_uM = "Cyclophosphamide",
        oh4_uM = "4OH-Cyclophosphamide"
      )
    )

  list(
    sim = out,
    long = long,
    last_dose_time = last_dose_time
  )
}

compute_metrics <- function(sim_bundle, analyte_name) {
  analyte_col <- if (identical(analyte_name, "Cyclophosphamide")) "CY" else "OH4"

  sim <- sim_bundle$sim |>
    filter(time_since_last_dose >= 0, time_since_last_dose <= 24)

  conc <- sim[[analyte_col]]
  time_h <- sim$time_since_last_dose
  keep_unique <- !duplicated(time_h)

  list(
    cmax = max(conc, na.rm = TRUE),
    c24 = stats::approx(time_h[keep_unique], conc[keep_unique], xout = 24, rule = 2)$y,
    auc24 = trap_auc(time_h, conc),
    thalf = estimate_terminal_half_life(time_h, conc)
  )
}

build_exposure_summary <- function(sim_bundle) {
  analytes <- c("Cyclophosphamide", "4OH-Cyclophosphamide")

  dplyr::bind_rows(lapply(analytes, function(analyte) {
    met <- compute_metrics(sim_bundle, analyte)

    tibble::tibble(
      Analyte = analyte,
      Cmax_uM = round(met$cmax, 2),
      C24h_uM = round(met$c24, 3),
      AUC0_24h_uMh = round(met$auc24, 2),
      Half_life_h = round(met$thalf, 2)
    )
  }))
}

build_paper_grid <- function(sex, bsa_m2) {
  preset_keys <- c("fixed_500", "fixed_1000", "nih_500", "nih_1000")

  dplyr::bind_rows(lapply(preset_keys, function(key) {
    preset <- paper_regimen(key)
    sim <- simulate_profile(
      dose_mode = preset$dose_mode,
      dose_value = preset$dose_value,
      sex = sex,
      bsa_m2 = bsa_m2,
      interval_days = 28,
      n_doses = 1,
      infusion_h = 1,
      horizon_h = 48
    )

    cy <- compute_metrics(sim, "Cyclophosphamide")
    oh4 <- compute_metrics(sim, "4OH-Cyclophosphamide")

    tibble::tibble(
      Regimen = preset$label,
      Total_dose_mg = round(total_dose_mg(preset$dose_mode, preset$dose_value, bsa_m2), 1),
      CY_Cmax_uM = round(cy$cmax, 2),
      CY_AUC0_24h = round(cy$auc24, 2),
      OH4_Cmax_uM = round(oh4$cmax, 2),
      OH4_AUC0_24h = round(oh4$auc24, 2)
    )
  }))
}

app_theme <- bs_theme(
  version = 5,
  bootswatch = "flatly",
  primary = "#8b5cf6"
) |> bs_add_rules("
  .metric-card {
    background: #f8f9fa;
    border-radius: 10px;
    padding: 15px;
    margin: 5px;
    text-align: center;
    border: 1px solid #dee2e6;
    min-height: 110px;
  }
  .metric-value {
    font-size: 24px;
    font-weight: 700;
    color: #1f2937;
  }
  .metric-label {
    font-size: 12px;
    color: #6b7280;
    letter-spacing: 0.04em;
    text-transform: uppercase;
  }
  .metric-success .metric-value { color: #10b981; }
  .metric-warning .metric-value { color: #f59e0b; }
  .metric-primary .metric-value { color: #8b5cf6; }
  .metric-info .metric-value { color: #0ea5e9; }
  .ref-box {
    background: #f5f3ff;
    border-left: 4px solid #8b5cf6;
    padding: 12px 14px;
    border-radius: 6px;
    font-size: 13px;
    line-height: 1.45;
    margin-top: 10px;
  }
  .sidebar-section-note {
    color: #475569;
    font-size: 13px;
    line-height: 1.45;
  }
  .app-caption {
    color: #64748b;
    font-size: 13px;
    line-height: 1.5;
  }
")

ui <- page_sidebar(
  title = "Cyclophosphamide Glomerulonephritis PK Simulator",
  theme = app_theme,
  sidebar = sidebar(
    width = 340,
    accordion(
      accordion_panel(
        "Regimen",
        selectInput(
          "preset",
          "Paper regimen preset",
          choices = c(
            "Custom" = "custom",
            "Fixed 500 mg" = "fixed_500",
            "Fixed 1000 mg" = "fixed_1000",
            "NIH 500 mg/m2" = "nih_500",
            "NIH 1000 mg/m2" = "nih_1000"
          ),
          selected = "nih_500"
        ),
        conditionalPanel(
          "input.preset === 'custom'",
          radioButtons(
            "dose_mode",
            "Dose units",
            choices = c("Fixed dose (mg)" = "mg", "BSA-based (mg/m2)" = "mg_m2"),
            selected = "mg_m2",
            inline = TRUE
          ),
          numericInput("dose_value", "Dose amount", value = 500, min = 50, max = 5000, step = 50)
        ),
        sliderInput("infusion_h", "Infusion duration (h)", min = 0.25, max = 4, value = 1, step = 0.25),
        sliderInput("interval_days", "Dosing interval (days)", min = 7, max = 56, value = 28, step = 1),
        sliderInput("n_doses", "Number of doses", min = 1, max = 8, value = 1, step = 1),
        textOutput("preset_note"),
        div(class = "sidebar-section-note", textOutput("dose_summary"))
      ),
      accordion_panel(
        "Patient",
        radioButtons("sex", "Sex", choices = c("Female", "Male"), selected = "Female", inline = TRUE),
        sliderInput("bsa_m2", "Body surface area (m2)", min = 1.2, max = 2.8, value = 2.0, step = 0.05),
        div(
          class = "sidebar-section-note",
          "Only sex altered the final model in the paper. Weight, BSA, creatinine clearance, and serum albumin were tested but not retained as covariates."
        )
      ),
      accordion_panel(
        "Plot",
        checkboxGroupInput(
          "analytes",
          "Analytes to display",
          choices = c("Cyclophosphamide", "4OH-Cyclophosphamide"),
          selected = c("Cyclophosphamide", "4OH-Cyclophosphamide")
        ),
        selectInput(
          "metric_analyte",
          "Metric cards focus",
          choices = c("Cyclophosphamide", "4OH-Cyclophosphamide"),
          selected = "Cyclophosphamide"
        ),
        sliderInput("horizon_h", "Hours after last dose", min = 12, max = 72, value = 48, step = 4),
        checkboxInput("log_scale", "Use log concentration scale", value = FALSE)
      )
    ),
    div(
      class = "ref-box",
      tags$strong("Paper anchor"),
      tags$br(),
      "Iliopoulou VN, Charkoftaki G, Cooper JC, Dokoumetzidis A, Joy MS.",
      tags$br(),
      tags$em("Population pharmacokinetics of cyclophosphamide and 4-hydroxycyclophosphamide metabolite in patients with autoimmune glomerulonephritis."),
      tags$br(),
      "J Pharm Pharmacol. 2021;73(12):1683-1692. doi:10.1093/jpp/rgab135"
    ),
    div(
      class = "app-caption",
      "Built for PKPDBuilder with a 2-compartment parent model plus 1-compartment 4OH metabolite model. ",
      "Metabolite volume is fixed to the parent central volume exactly as described in the paper."
    ),
    div(
      class = "app-caption",
      tags$a(href = "https://pkpd-builder.vercel.app", target = "_blank", "PKPDBuilder")
    )
  ),
  layout_column_wrap(
    width = 1 / 4,
    fill = FALSE,
    div(
      class = "metric-card metric-success",
      div(class = "metric-value", textOutput("cmax")),
      div(class = "metric-label", "Cmax (uM)")
    ),
    div(
      class = "metric-card metric-warning",
      div(class = "metric-value", textOutput("c24")),
      div(class = "metric-label", "C24h (uM)")
    ),
    div(
      class = "metric-card metric-primary",
      div(class = "metric-value", textOutput("auc24")),
      div(class = "metric-label", "AUC0-24h (uM*h)")
    ),
    div(
      class = "metric-card metric-info",
      div(class = "metric-value", textOutput("thalf")),
      div(class = "metric-label", "t1/2 (h)")
    )
  ),
  card(
    full_screen = TRUE,
    card_header("Simulated concentration-time profile"),
    plotOutput("pk_plot", height = "500px")
  ),
  navset_card_underline(
    nav_panel("Exposure Summary", tableOutput("summary_table")),
    nav_panel("Paper Scenario Grid", tableOutput("paper_grid"))
  )
)

server <- function(input, output, session) {
  preset_config <- reactive({
    if (identical(input$preset, "custom")) {
      return(NULL)
    }

    paper_regimen(input$preset)
  })

  active_dose_mode <- reactive({
    cfg <- preset_config()
    if (is.null(cfg)) input$dose_mode else cfg$dose_mode
  })

  active_dose_value <- reactive({
    cfg <- preset_config()
    if (is.null(cfg)) input$dose_value else cfg$dose_value
  })

  output$preset_note <- renderText({
    cfg <- preset_config()
    if (is.null(cfg)) {
      "Custom regimen mode."
    } else {
      cfg$note
    }
  })

  output$dose_summary <- renderText({
    shiny::req(active_dose_mode(), active_dose_value(), input$bsa_m2)

    total_mg <- total_dose_mg(active_dose_mode(), active_dose_value(), input$bsa_m2)
    dose_label <- if (identical(active_dose_mode(), "mg_m2")) {
      paste0(round(active_dose_value(), 1), " mg/m2")
    } else {
      paste0(round(active_dose_value(), 1), " mg")
    }

    paste0(
      "Active dose: ", dose_label,
      "  |  Total delivered: ", round(total_mg, 1), " mg"
    )
  })

  sim_bundle <- reactive({
    shiny::req(input$sex, input$bsa_m2, input$interval_days, input$n_doses, input$infusion_h, input$horizon_h)

    simulate_profile(
      dose_mode = active_dose_mode(),
      dose_value = active_dose_value(),
      sex = input$sex,
      bsa_m2 = input$bsa_m2,
      interval_days = input$interval_days,
      n_doses = input$n_doses,
      infusion_h = input$infusion_h,
      horizon_h = input$horizon_h
    )
  })

  focus_metrics <- reactive({
    shiny::req(input$metric_analyte)
    compute_metrics(sim_bundle(), input$metric_analyte)
  })

  output$cmax <- renderText(sprintf("%.2f", focus_metrics()$cmax))
  output$c24 <- renderText(sprintf("%.3f", focus_metrics()$c24))
  output$auc24 <- renderText(sprintf("%.2f", focus_metrics()$auc24))
  output$thalf <- renderText({
    val <- focus_metrics()$thalf
    if (is.na(val)) "NA" else sprintf("%.2f", val)
  })

  output$pk_plot <- renderPlot({
    shiny::req(input$analytes)
    plot_df <- sim_bundle()$long |>
      filter(analyte %in% input$analytes)

    if (isTRUE(input$log_scale)) {
      plot_df <- plot_df |>
        filter(conc_uM > 0)
    }

    plot_obj <- ggplot(plot_df, aes(x = time_since_last_dose, y = conc_uM, color = analyte)) +
      geom_line(linewidth = 1.1) +
      scale_color_manual(
        values = c(
          "Cyclophosphamide" = "#2563eb",
          "4OH-Cyclophosphamide" = "#dc2626"
        )
      ) +
      labs(
        x = "Hours since last dose",
        y = if (isTRUE(input$log_scale)) "Concentration (uM, log scale)" else "Concentration (uM)",
        color = NULL,
        title = paste0("Post-dose profile for ", if (identical(input$sex, "Female")) "female" else "male", " virtual patient"),
        subtitle = paste0("BSA ", sprintf("%.2f", input$bsa_m2), " m2")
      ) +
      theme_minimal(base_size = 13) +
      theme(
        legend.position = "top",
        plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank()
      ) +
      NULL

    if (isTRUE(input$log_scale)) {
      plot_obj <- plot_obj + scale_y_log10()
    }

    plot_obj
  })

  output$summary_table <- renderTable({
    build_exposure_summary(sim_bundle())
  }, striped = TRUE, hover = TRUE, spacing = "m")

  output$paper_grid <- renderTable({
    build_paper_grid(input$sex, input$bsa_m2)
  }, striped = TRUE, hover = TRUE, spacing = "m")
}

shinyApp(ui, server)