# ============================================================
#  BOOTSTRAP PARAMÉTRICO
# ============================================================

library(MASS)

 
line <- function(x, intercept, slope) intercept + slope * x

sigmoide3 <- function(x, top, bottom, lec50) {
  bottom + (top - bottom) / (1 + 10^(lec50 - x))
}

find_roots <- function(f, x_min, x_max, n_points = 1000, tol = 1e-6) {
  x_grid <- seq(x_min, x_max, length.out = n_points)
  y_grid <- f(x_grid)
  roots <- numeric(0)
  for (i in 1:(length(x_grid) - 1)) {
    if (y_grid[i] * y_grid[i + 1] < 0) {
      root <- tryCatch(uniroot(f, c(x_grid[i], x_grid[i + 1]), tol = tol)$root,
                       error = function(e) NA)
      if (!is.na(root)) roots <- c(roots, root)
    }
  }
  if (length(roots) > 1) {
    roots <- roots[order(roots)]
    roots <- roots[c(TRUE, diff(roots) > tol * 10)]
  }
  return(roots)
}


fit_linear_dias <- function(x, y1, y2) {
  m1 <- lm(y1 ~ x)
  m2 <- lm(y2 ~ x)
  p1 <- coef(m1); p2 <- coef(m2)
  c1 <- vcov(m1); c2 <- vcov(m2)
  
  denom <- p1[2] - p2[2]
  if (abs(denom) < 1e-12) return(NULL)
  dias_x <- (p2[1] - p1[1]) / denom
  dias_y <- p1[1] + p1[2] * dias_x
  
  h <- 1e-6
  J <- matrix(0, 2, 4)
  mt <- c(p1, p2)
  for (j in 1:4) {
    mp <- mt; mp[j] <- mp[j] + h
    denom_pert <- mp[2] - mp[4]
    if (abs(denom_pert) < 1e-12) return(NULL)
    root_pert <- (mp[3] - mp[1]) / denom_pert
    J[1, j] <- (root_pert - dias_x) / h
  }
  for (j in 1:2) {
    mp <- mt; mp[j] <- mp[j] + h
    y_pert <- mp[1] + mp[2] * dias_x
    J[2, j] <- (y_pert - dias_y) / h
  }
  
  Cov_all <- matrix(0, 4, 4)
  Cov_all[1:2, 1:2] <- c1
  Cov_all[3:4, 3:4] <- c2
  Cov_XY <- J %*% Cov_all %*% t(J)
  
  return(list(mu = c(dias_x, dias_y), Sigma = Cov_XY))
}

fit_sigmoid_dias <- function(x, y1, y2,
                             start1 = list(top = 100, bottom = 0, lec50 = -7),
                             start2 = list(top = 80, bottom = 20, lec50 = -6.5)) {
  m1 <- try(nls(y1 ~ bottom + (top - bottom) / (1 + 10^(lec50 - x)),
                start = start1,
                control = nls.control(maxiter = 500, minFactor = 1e-10)),
            silent = TRUE)
  m2 <- try(nls(y2 ~ bottom + (top - bottom) / (1 + 10^(lec50 - x)),
                start = start2,
                control = nls.control(maxiter = 500, minFactor = 1e-10)),
            silent = TRUE)
  if (inherits(m1, "try-error") || inherits(m2, "try-error")) return(NULL)
  
  p1 <- coef(m1); p2 <- coef(m2)
  c1 <- vcov(m1); c2 <- vcov(m2)
  
  diff_func <- function(x) sigmoide3(x, p1[1], p1[2], p1[3]) - sigmoide3(x, p2[1], p2[2], p2[3])
  
  x_center <- mean(c(p1[3], p2[3]))
  x_range <- 4
  x_min <- min(c(p1[3], p2[3])) - x_range
  x_max <- max(c(p1[3], p2[3])) + x_range
  
  roots <- find_roots(diff_func, x_min, x_max)
  if (length(roots) == 0) return(NULL)
  dias_x <- roots[1]
  dias_y <- sigmoide3(dias_x, p1[1], p1[2], p1[3])
  
  h <- 1e-6
  J <- matrix(0, 2, 6)
  mt <- c(p1, p2)
  for (j in 1:6) {
    mp <- mt; mp[j] <- mp[j] + h
    diff_pert <- function(x) sigmoide3(x, mp[1], mp[2], mp[3]) - sigmoide3(x, mp[4], mp[5], mp[6])
    root_pert <- tryCatch(uniroot(diff_pert, c(dias_x - 0.5, dias_x + 0.5), tol = 1e-6)$root,
                          error = function(e) dias_x)
    J[1, j] <- (root_pert - dias_x) / h
  }
  for (j in 1:3) {
    mp <- mt; mp[j] <- mp[j] + h
    y_pert <- sigmoide3(dias_x, mp[1], mp[2], mp[3])
    J[2, j] <- (y_pert - dias_y) / h
  }
  
  Cov_all <- matrix(0, 6, 6)
  Cov_all[1:3, 1:3] <- c1
  Cov_all[4:6, 4:6] <- c2
  Cov_XY <- J %*% Cov_all %*% t(J)
  
  return(list(mu = c(dias_x, dias_y), Sigma = Cov_XY))
}


wald_W <- function(dias1, dias2) {
  delta <- dias1$mu - dias2$mu
  cov_combined <- dias1$Sigma + dias2$Sigma
  W <- as.numeric(t(delta) %*% solve(cov_combined) %*% delta)
  return(W)
}


bootstrap_pvalue <- function(x, y1_g1, y2_g1, y1_g2, y2_g2,
                             model = "linear", B = 300, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  
  if (model == "linear") {
    m1_g1 <- lm(y1_g1 ~ x); m2_g1 <- lm(y2_g1 ~ x)
    m1_g2 <- lm(y1_g2 ~ x); m2_g2 <- lm(y2_g2 ~ x)
    beta1 <- (coef(m1_g1) + coef(m1_g2)) / 2
    beta2 <- (coef(m2_g1) + coef(m2_g2)) / 2
    n <- length(x)
    sig1 <- sqrt((sum(residuals(m1_g1)^2) + sum(residuals(m1_g2)^2)) / (2*n - 4))
    sig2 <- sqrt((sum(residuals(m2_g1)^2) + sum(residuals(m2_g2)^2)) / (2*n - 4))
    d1_obs <- fit_linear_dias(x, y1_g1, y2_g1)
    d2_obs <- fit_linear_dias(x, y1_g2, y2_g2)
  } else {
    m1_g1 <- try(nls(y1_g1 ~ bottom + (top - bottom) / (1 + 10^(lec50 - x)),
                     start = list(top = 100, bottom = 0, lec50 = -7),
                     control = nls.control(maxiter = 500, minFactor = 1e-10)),
                 silent = TRUE)
    m2_g1 <- try(nls(y2_g1 ~ bottom + (top - bottom) / (1 + 10^(lec50 - x)),
                     start = list(top = 80, bottom = 20, lec50 = -6.5),
                     control = nls.control(maxiter = 500, minFactor = 1e-10)),
                 silent = TRUE)
    m1_g2 <- try(nls(y1_g2 ~ bottom + (top - bottom) / (1 + 10^(lec50 - x)),
                     start = list(top = 100, bottom = 0, lec50 = -7),
                     control = nls.control(maxiter = 500, minFactor = 1e-10)),
                 silent = TRUE)
    m2_g2 <- try(nls(y2_g2 ~ bottom + (top - bottom) / (1 + 10^(lec50 - x)),
                     start = list(top = 80, bottom = 20, lec50 = -6.5),
                     control = nls.control(maxiter = 500, minFactor = 1e-10)),
                 silent = TRUE)
    if (any(inherits(m1_g1, "try-error"), inherits(m2_g1, "try-error"),
            inherits(m1_g2, "try-error"), inherits(m2_g2, "try-error"))) return(NA)
    beta1 <- (coef(m1_g1) + coef(m1_g2)) / 2
    beta2 <- (coef(m2_g1) + coef(m2_g2)) / 2
    n <- length(x)
    sig1 <- sqrt((sum(residuals(m1_g1)^2) + sum(residuals(m1_g2)^2)) / (2*n - 6))
    sig2 <- sqrt((sum(residuals(m2_g1)^2) + sum(residuals(m2_g2)^2)) / (2*n - 6))
    d1_obs <- fit_sigmoid_dias(x, y1_g1, y2_g1)
    d2_obs <- fit_sigmoid_dias(x, y1_g2, y2_g2)
  }
  
  if (is.null(d1_obs) || is.null(d2_obs)) return(NA)
  W_obs <- wald_W(d1_obs, d2_obs)
  
  W_boot <- numeric(B)
  for (b in 1:B) {
    if (model == "linear") {
      y1_b1 <- beta1[1] + beta1[2] * x + rnorm(n, 0, sig1)
      y2_b1 <- beta2[1] + beta2[2] * x + rnorm(n, 0, sig2)
      y1_b2 <- beta1[1] + beta1[2] * x + rnorm(n, 0, sig1)
      y2_b2 <- beta2[1] + beta2[2] * x + rnorm(n, 0, sig2)
      db1 <- fit_linear_dias(x, y1_b1, y2_b1)
      db2 <- fit_linear_dias(x, y1_b2, y2_b2)
    } else {
      y1_b1 <- sigmoide3(x, beta1[1], beta1[2], beta1[3]) + rnorm(n, 0, sig1)
      y2_b1 <- sigmoide3(x, beta2[1], beta2[2], beta2[3]) + rnorm(n, 0, sig2)
      y1_b2 <- sigmoide3(x, beta1[1], beta1[2], beta1[3]) + rnorm(n, 0, sig1)
      y2_b2 <- sigmoide3(x, beta2[1], beta2[2], beta2[3]) + rnorm(n, 0, sig2)
      db1 <- fit_sigmoid_dias(x, y1_b1, y2_b1)
      db2 <- fit_sigmoid_dias(x, y1_b2, y2_b2)
    }
    if (is.null(db1) || is.null(db2)) next
    W_boot[b] <- wald_W(db1, db2)
  }
  W_boot <- W_boot[!is.na(W_boot)]
  if (length(W_boot) == 0) return(NA)
  mean(W_boot >= W_obs)
}

# ---------- Simulação
run_mc_bootstrap <- function(model = "linear", n = 10, sigma = 1,
                             H1 = FALSE, delta = 2,
                             n_mc = 300, B = 300) {
  rejections <- 0
  valid <- 0
  
  if (model == "linear") {
    a1 <- 0; b1 <- 1
    a2 <- 10; b2 <- -1
  } else {
    top1 <- 100; bot1 <- 0; lec1 <- -7
    top2 <- 80;  bot2 <- 20; lec2 <- -6.5
  }
  
  for (i in 1:n_mc) {
    if (model == "linear") {
      x <- seq(0, 10, length.out = n)
      if (H1) {
        y1_g1 <- a1 + b1 * x + rnorm(n, 0, sigma)
        y2_g1 <- (a2 + delta) + b2 * x + rnorm(n, 0, sigma)
        y1_g2 <- a1 + b1 * x + rnorm(n, 0, sigma)
        y2_g2 <- a2 + b2 * x + rnorm(n, 0, sigma)
      } else {
        y1_g1 <- a1 + b1 * x + rnorm(n, 0, sigma)
        y2_g1 <- a2 + b2 * x + rnorm(n, 0, sigma)
        y1_g2 <- a1 + b1 * x + rnorm(n, 0, sigma)
        y2_g2 <- a2 + b2 * x + rnorm(n, 0, sigma)
      }
    } else {
      x <- seq(-9, -5, length.out = n)
      if (H1) {
        y1_g1 <- sigmoide3(x, top1, bot1, lec1) + rnorm(n, 0, sigma)
        y2_g1 <- sigmoide3(x, top2, bot2, lec2 + delta) + rnorm(n, 0, sigma)
        y1_g2 <- sigmoide3(x, top1, bot1, lec1) + rnorm(n, 0, sigma)
        y2_g2 <- sigmoide3(x, top2, bot2, lec2) + rnorm(n, 0, sigma)
      } else {
        y1_g1 <- sigmoide3(x, top1, bot1, lec1) + rnorm(n, 0, sigma)
        y2_g1 <- sigmoide3(x, top2, bot2, lec2) + rnorm(n, 0, sigma)
        y1_g2 <- sigmoide3(x, top1, bot1, lec1) + rnorm(n, 0, sigma)
        y2_g2 <- sigmoide3(x, top2, bot2, lec2) + rnorm(n, 0, sigma)
      }
    }
    
    p <- bootstrap_pvalue(x, y1_g1, y2_g1, y1_g2, y2_g2,
                          model = model, B = B)
    if (is.na(p)) next
    valid <- valid + 1
    if (p < 0.05) rejections <- rejections + 1
  }
  rejections / valid
}


scenarios <- expand.grid(
  model = c("linear", "sigmoid"),
  n = c(6, 20),
  sigma = c(1, 5)
)

results <- data.frame()

for (s in 1:nrow(scenarios)) {
  model <- scenarios$model[s]
  n <- scenarios$n[s]
  sigma <- scenarios$sigma[s]
  
  cat(sprintf("Model=%s, n=%d, sigma=%.1f: H0...\n", model, n, sigma))
  h0 <- run_mc_bootstrap(model, n, sigma, H1 = FALSE)
  
  delta_mod <- if (model == "linear") 2 else 0.5
  cat(sprintf("Model=%s, n=%d, sigma=%.1f: H1 moderate...\n", model, n, sigma))
  h1_mod <- run_mc_bootstrap(model, n, sigma, H1 = TRUE, delta = delta_mod)
  
  delta_strong <- if (model == "linear") 5 else 1.0
  cat(sprintf("Model=%s, n=%d, sigma=%.1f: H1 strong...\n", model, n, sigma))
  h1_strong <- run_mc_bootstrap(model, n, sigma, H1 = TRUE, delta = delta_strong)
  
  results <- rbind(results, data.frame(
    model = model, n = n, sigma = sigma,
    H0_rate = h0,
    H1_moderate_power = h1_mod,
    H1_strong_power = h1_strong
  ))
}

print(results, row.names = FALSE)
write.csv(results, "bootstrap_validation_results.csv", row.names = FALSE)