# SENSIBILIDADE STEP SIZE h 
#  sigmoide 3P


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_sigmoid_dias_h <- function(x, y1, y2, h,
                               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])
  
  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(dias_x = dias_x,
              dias_y = dias_y,
              se_x = sqrt(Cov_XY[1,1]),
              se_y = sqrt(Cov_XY[2,2])))
}


set.seed(123)
n <- 10
x <- seq(-9, -5, length.out = n)
y1 <- sigmoide3(x, 100, 0, -7) + rnorm(n, 0, 2)
y2 <- sigmoide3(x, 80, 20, -6.5) + rnorm(n, 0, 2)

h_values <- c(1e-4, 1e-5, 1e-6, 1e-7, 1e-8)
results <- data.frame()

for (h in h_values) {
  res <- fit_sigmoid_dias_h(x, y1, y2, h)
  if (!is.null(res)) {
    results <- rbind(results, data.frame(
      h = h,
      dias_x = res$dias_x,
      dias_y = res$dias_y,
      se_x = res$se_x,
      se_y = res$se_y
    ))
  }
}

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