# ============================================================
# MONTE CARLO
#  linear, sigmoide 3P, cosseno
# ============================================================

library(MASS)


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

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

cos_model <- function(x, mesor, amplitude, acrophase, period = 12) {
  mesor + amplitude * cos(2 * pi * (x - acrophase) / period)
}

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))
}

fit_cosine_dias <- function(x, y1, y2, period = 12) {
  start1 <- list(mesor = mean(y1), amplitude = (max(y1) - min(y1)) / 2, acrophase = 0)
  start2 <- list(mesor = mean(y2), amplitude = (max(y2) - min(y2)) / 2, acrophase = 0)
  
  m1 <- try(nls(y1 ~ mesor + amplitude * cos(2 * pi * (x - acrophase) / period),
                start = start1, control = nls.control(maxiter = 500)), silent = TRUE)
  m2 <- try(nls(y2 ~ mesor + amplitude * cos(2 * pi * (x - acrophase) / period),
                start = start2, control = nls.control(maxiter = 500)), 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) cos_model(x, p1[1], p1[2], p1[3], period) - cos_model(x, p2[1], p2[2], p2[3], period)
  roots <- find_roots(diff_func, 0, 12)
  if (length(roots) == 0) return(NULL)
  dias_x <- roots[1]
  dias_y <- cos_model(dias_x, p1[1], p1[2], p1[3], period)
  
  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) cos_model(x, mp[1], mp[2], mp[3], period) - cos_model(x, mp[4], mp[5], mp[6], period)
    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 <- cos_model(dias_x, mp[1], mp[2], mp[3], period)
    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))
}


true_linear <- list(a1 = 0, b1 = 1, a2 = 10, b2 = -1)   # (5,5)
true_sigmoid <- list(top1 = 100, bot1 = 0, lec1 = -7,
                     top2 = 80,  bot2 = 20, lec2 = -6.5) # (~-7.45, 26.0)
true_cosine <- list(mesor1 = 20, amp1 = 5, phase1 = 2,
                    mesor2 = 25, amp2 = 3, phase2 = 8)    # (~0.29, 23.1)

# ----------  Monte Carlo
run_mc_estimator <- function(model = "linear", n = 10, sigma = 1,
                             n_sim = 1000, seed = 2024) {
  set.seed(seed)
  
 
  if (model == "linear") {
    f1 <- function(x) true_linear$a1 + true_linear$b1 * x
    f2 <- function(x) true_linear$a2 + true_linear$b2 * x
    true_x <- (true_linear$a2 - true_linear$a1) / (true_linear$b1 - true_linear$b2)
    true_y <- f1(true_x)
    x_range <- c(0, 10)
  } else if (model == "sigmoid") {
    f1 <- function(x) sigmoide3(x, true_sigmoid$top1, true_sigmoid$bot1, true_sigmoid$lec1)
    f2 <- function(x) sigmoide3(x, true_sigmoid$top2, true_sigmoid$bot2, true_sigmoid$lec2)
    true_x <- uniroot(function(x) f1(x) - f2(x), c(-9, -5))$root
    true_y <- f1(true_x)
    x_range <- c(-9, -5)
  } else if (model == "cosine") {
    f1 <- function(x) cos_model(x, true_cosine$mesor1, true_cosine$amp1, true_cosine$phase1)
    f2 <- function(x) cos_model(x, true_cosine$mesor2, true_cosine$amp2, true_cosine$phase2)
    true_x <- find_roots(function(x) f1(x) - f2(x), 0, 12)[1]
    true_y <- f1(true_x)
    x_range <- c(0, 12)
  }
  
  x_vals <- seq(x_range[1], x_range[2], length.out = n)
  
  est_x <- numeric(n_sim); est_y <- numeric(n_sim)
  se_x <- numeric(n_sim); se_y <- numeric(n_sim)
  cover_x <- logical(n_sim); cover_y <- logical(n_sim)
  conv <- logical(n_sim)
  
  for (i in 1:n_sim) {
    y1 <- f1(x_vals) + rnorm(n, 0, sigma)
    y2 <- f2(x_vals) + rnorm(n, 0, sigma)
    
    res <- switch(model,
                  linear  = fit_linear_dias(x_vals, y1, y2),
                  sigmoid = fit_sigmoid_dias(x_vals, y1, y2),
                  cosine  = fit_cosine_dias(x_vals, y1, y2))
    
    if (is.null(res)) next
    conv[i] <- TRUE
    est_x[i] <- res$mu[1]
    est_y[i] <- res$mu[2]
    se_x[i]  <- sqrt(res$Sigma[1,1])
    se_y[i]  <- sqrt(res$Sigma[2,2])
    
    ci_x <- est_x[i] + c(-1.96, 1.96) * se_x[i]
    ci_y <- est_y[i] + c(-1.96, 1.96) * se_y[i]
    cover_x[i] <- (true_x >= ci_x[1]) & (true_x <= ci_x[2])
    cover_y[i] <- (true_y >= ci_y[1]) & (true_y <= ci_y[2])
  }
  
  valid <- which(conv)
  data.frame(
    model = model,
    n = n,
    sigma = sigma,
    conv = sum(conv),
    bias_x = mean(est_x[valid]) - true_x,
    bias_y = mean(est_y[valid]) - true_y,
    rmse_x = sqrt(mean((est_x[valid] - true_x)^2)),
    rmse_y = sqrt(mean((est_y[valid] - true_y)^2)),
    cov_x = mean(cover_x[valid]),
    cov_y = mean(cover_y[valid])
  )
}


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

results <- do.call(rbind, lapply(1:nrow(scenarios), function(i) {
  cat(sprintf("Running %s n=%d sigma=%.1f\n",
              scenarios$model[i], scenarios$n[i], scenarios$sigma[i]))
  run_mc_estimator(scenarios$model[i], scenarios$n[i], scenarios$sigma[i])
}))

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