# Install from GitHub
install.packages("devtools")
devtools::install_github("cuihan1997/SensOverlap")


library(SensOverlap)

# Generate simulation data
gen.sample <- function(n) {
  samp <- data.frame()
  for (i in 1:n) {
    X <- runif(1, 0, 1)
    if (X <= 0.7) {
      U <- runif(1, 0, 1)
    } else {
      U <- runif(1, 0, 100)
    }
    Z_value <- rbinom(1, 1, 1/(1 + exp(X - 0.1*U)))
    Y_value <- 2*X + 3*U + 5*Z_value
    samp_add <- data.frame(Y = Y_value, Z = Z_value, X = X, U = U)
    samp <- rbind(samp, samp_add)
  }
  return(samp)
}

# Generate data
set.seed(316)
n <- 50
data_sim <- gen.sample(n)

Z_raw <- data_sim["Z"]
Y_raw <- data_sim["Y"]
X_raw <- data_sim["X"]

# Set sensitivity parameters
Lambda1 <- exp(1)
Lambda2 <-  exp(2)
delta1 <- 0.1
delta2 <- 0.2

# Define constraint and transformation functions
g_fun1 <- function(x) { x }
g_functions <- list(g_fun1)

s_fun1 <- function(x) { x }
s_functions <- list(s_fun1)

# Generate bootstrap weights
set.seed(316)
K <- as.numeric(rmultinom(1, dim(Z_raw)[1], prob = rep(1, dim(Z_raw)[1])))

# Create sensitivity analysis object
ow_model <- overlap.sen(
  Z_raw, X_raw, Y_raw, g_functions, 
  encoded_prop_score_information = list(get.fitted_prop_score, s_functions), 
  decode_method = decode.prop_score_information
)

# Compute sensitivity bounds
ow_model$ate.seperate.MILP(
    K, 
    c(Lambda1, Lambda2), 
    c(delta1, delta2), 
    M = 1e9, 
    integer_constraints_relaxed = TRUE
)



# Bootstrap Confidence Intervals
ow_model$boot_ci.seperate.MILP(
  c(Lambda1, Lambda2), 
  c(delta1, delta2), 
  alpha = 0.05, 
  B = 1000, 
  side = "LU", 
  base_seed = 2025
)