# Fits the mechanism core's turnover half-lives and Pmax against the OPTIC
# trial's digitized proptosis-change-from-baseline curve
# (data/digitized/optic_trial_proptosis_change.csv) -- the only continuous,
# visit-by-visit endpoint in that trial that maps directly onto a model
# output (Proptosis_mm). See data/digitized/README.md and model-spec.md
# Section 8 for why CAS isn't fit the same way (no continuous mean-CAS
# table was reported, only a binary responder rate).
#
# Only 2 parameters are fit (effector_halflife, Pmax) against 8 points (2
# arms x 4 visits):
# - M50/A50 are FIXED at the structural default of 1 (not fit): an earlier
#   attempt to also fit M50 hit a hard identifiability wall -- since
#   Myo/Fat range roughly 0-1 in this model (normalized to 1 at
#   steady-state baseline drive), M50 >> 1 makes M/(M+M50) ~= M/M50 for the
#   entire trajectory, so Pmax and M50 are only jointly identified via
#   their ratio, never separately.
# - F_halflife is FIXED at 3 days (not fit): letting it float, it
#   repeatedly pushed to whatever the lower search bound was (1 day, then
#   even lower) -- the 8 visit-level data points can't distinguish "1 day"
#   from "instantaneous," so it isn't its own identifiable timescale here.
#   Once F is fast enough to just pass EffDrive through with negligible
#   delay relative to the ~3-week visit spacing, only the slower
#   effector_halflife (which the data DOES pin down cleanly, away from any
#   bound) matters for the visible trajectory shape.
#
# wM/wA are likewise not fit: with a single shared effector half-life for
# both the myofibroblast and adipocyte pools (M50 = A50), Myo(t) == Fat(t)
# identically, making the wM/wA split unidentifiable from proptosis alone
# (any split with wM+wA=1 gives the same sum). Left at 0.5/0.5.

source("../R/load_all.R", chdir = TRUE)

obs <- read.csv("../data/digitized/optic_trial_proptosis_change.csv")
obs <- obs[obs$week > 0, ]  # week 0 is 0 by construction, not informative

weeks <- c(6, 12, 18, 24)
days <- weeks * 7

F_HALFLIFE_FIXED <- 3

predict_changes <- function(par) {
  effector_halflife <- exp(par[1])
  Pmax <- exp(par[2])

  run <- function(drug_id) {
    sol <- ted_simulate(drug_id,
                         body_weight_kg = 70, trab_baseline = 1.0,
                         disease_duration_days = 190,
                         F_halflife = F_HALFLIFE_FIXED, effector_halflife = effector_halflife,
                         Pmax = Pmax, M50 = 1, A50 = 1,
                         sim_days = 170, dt = 1)
    baseline <- sol$Proptosis_mm[sol$time == 0]
    vapply(days, function(d) sol$Proptosis_mm[which.min(abs(sol$time - d))] - baseline, numeric(1))
  }

  data.frame(week = weeks,
             teprotumumab = run("teprotumumab"),
             placebo = run("none"))
}

objective <- function(par) {
  pred <- predict_changes(par)
  obs_tepro <- obs$mean_change_mm[obs$arm == "teprotumumab"]
  obs_pbo   <- obs$mean_change_mm[obs$arm == "placebo"]
  sum((pred$teprotumumab - obs_tepro)^2) + sum((pred$placebo - obs_pbo)^2)
}

start <- log(c(effector_halflife = 60, Pmax = 3))
lower <- log(c(effector_halflife = 5, Pmax = 0.5))
upper <- log(c(effector_halflife = 150, Pmax = 30))
fit <- optim(start, objective, method = "L-BFGS-B", lower = lower, upper = upper,
             control = list(maxit = 2000, factr = 1e7))

cat("convergence:", fit$convergence, " objective (sum sq. mm^2):", fit$value, "\n")
fitted <- exp(fit$par)
names(fitted) <- c("effector_halflife", "Pmax")
print(fitted)
cat("at bounds? effector:", fitted[["effector_halflife"]] %in% c(5, 150),
    " Pmax:", fitted[["Pmax"]] %in% c(0.5, 30), "\n")

pred <- predict_changes(fit$par)
cat("\nFitted vs observed proptosis change (mm):\n")
long_pred <- rbind(
  data.frame(week = pred$week, arm = "teprotumumab", predicted = pred$teprotumumab),
  data.frame(week = pred$week, arm = "placebo", predicted = pred$placebo)
)
comparison <- merge(long_pred, obs[, c("week", "arm", "mean_change_mm")], by = c("week", "arm"))
comparison <- comparison[order(comparison$arm, comparison$week), ]
print(comparison, row.names = FALSE)

rmse <- sqrt(mean((comparison$predicted - comparison$mean_change_mm)^2))
cat("\nRMSE:", round(rmse, 3), "mm\n")

saveRDS(fitted, "calibrated_mechanism_core_params.rds")