source('./crt.R')


#Simulation 1

num_iters<- 1000


run_exps <- function(num_individuals,num_steps,lag,effect_vector){

    vec1 <- NULL
    vec2 <- NULL
    vec3 <- NULL
    vec4 <- NULL
    x <- rep(x = NA, times = num_iters)
    pb <- txtProgressBar(0, length(x), style = 3)
    for (seed in 1:num_iters){
      setTxtProgressBar(pb, seed)
      x[seed] 
      Sys.sleep(time = 0.001)
      
      data_xyz <- synthetic(num_individuals,num_steps,effect_vector,seed,normalize=TRUE,include_x=TRUE) 

      data_list_1 <- data_swd(data_xyz,lag,sample_split=FALSE)
      p_list_1 <- crt_swd(data_list_1,type='p_values')
      result1 <- global_test(data_list_1, p_list_1, method='bonferroni')
      result2 <- global_test(data_list_1, p_list_1, method='zscore')

      
      data_list_2 <-  data_swd(data_xyz,lag,sample_split=TRUE)
      p_list_2 <- crt_swd(data_list_2,type='p_values')
      result3 <- global_test(data_list_2, p_list_2, method='fisher')
      result4 <- global_test(data_list_2, p_list_2, method='zscore')
      
      vec1 <- c(vec1,result1) 
      vec2 <- c(vec2,result2) 
      vec3 <- c(vec3,result3) 
      vec4 <- c(vec4,result4) 
    }
  return(list(vec1,vec2,vec3,vec4))
}





#Number of units (Figure 3 a&d)

sample_list <- c(100, 200,300,400,500)
num_steps <- 8
lag <- 2
effect_vector <- rep(0, num_steps)

for (i in 1:2){
  
  if (i==2){
    effect_vector[lag+1] = 0.03
  }
  mat1 <- matrix(0, num_iters,5)
  mat2 <- matrix(0, num_iters,5)
  mat3 <- matrix(0, num_iters,5)
  mat4 <- matrix(0, num_iters,5)
  for (j in 1:5){
    num_individuals <- sample_list[j]
    cat("\n","num_individuals=",num_individuals,"\n")
    out <- run_exps(num_individuals,num_steps,lag,effect_vector)
    mat1[,j] <- out[[1]]
    mat2[,j] <- out[[2]]
    mat3[,j] <- out[[3]]
    mat4[,j] <- out[[4]]
    cat("\n")
    cat("rejection_rate_1:",colMeans(mat1))
    cat(" rejection_rate_2:",colMeans(mat2))
    cat(" rejection_rate_3:",colMeans(mat3))
    cat(" rejection_rate_4:",colMeans(mat4))
  }
  save(mat1,mat2,mat3,mat4, file= paste0("results/sim_1_units_",i, ".RData"))
}



#Number of time steps (Figure 3 b&e)

num_individuals <- 300
steps_list <- c(4,6,8,10,12)
lag <- 2
effect_vector <- rep(0, num_steps)

for (i in 1:2){
  
  mat1 <- matrix(0, num_iters,5)
  mat2 <- matrix(0, num_iters,5)
  mat3 <- matrix(0, num_iters,5)
  mat4 <- matrix(0, num_iters,5)
  for (j in 1:5){
    num_steps<- steps_list[j]
    cat("\n","num_steps=",num_steps,"\n")
    if (i==2){
      effect_vector <- rep(0, num_steps)
      effect_vector[lag+1] = 0.03
    }
    out <- run_exps(num_individuals,num_steps,lag,effect_vector)
    mat1[,j] <- out[[1]]
    mat2[,j] <- out[[2]]
    mat3[,j] <- out[[3]]
    mat4[,j] <- out[[4]]
    cat("\n")
    cat("rejection_rate_1:",colMeans(mat1))
    cat(" rejection_rate_2:",colMeans(mat2))
    cat(" rejection_rate_3:",colMeans(mat3))
    cat(" rejection_rate_4:",colMeans(mat4))
  }
  save(mat1,mat2,mat3,mat4, file= paste0("results/sim_1_steps_",i, ".RData"))
}







#Time lags (Figure 3 c&f)

num_individuals <- 300
num_steps <- 8
lag_list <- 0:4
effect_vector <- rep(0, num_steps)

for (i in 1:2){
  
  mat1 <- matrix(0, num_iters,5)
  mat2 <- matrix(0, num_iters,5)
  mat3 <- matrix(0, num_iters,5)
  mat4 <- matrix(0, num_iters,5)
  for (j in 1:5){
    lag <- lag_list[j]
    
    if (i==1){
      effect_vector <- rep(0, num_steps)
    }else {
      effect_vector <- rep(0, num_steps)
      effect_vector[lag+1] = 0.03
    }
    cat("\n","lag=",lag,"\n")
    out <- run_exps(num_individuals,num_steps,lag,effect_vector)
    mat1[,j] <- out[[1]]
    mat2[,j] <- out[[2]]
    mat3[,j] <- out[[3]]
    mat4[,j] <- out[[4]]
    cat("\n")
    cat("rejection_rate_1:",colMeans(mat1))
    cat(" rejection_rate_2:",colMeans(mat2))
    cat(" rejection_rate_3:",colMeans(mat3))
    cat(" rejection_rate_4:",colMeans(mat4))
  }
  save(mat1,mat2,mat3,mat4, file= paste0("results/sim_1_lags_",i, ".RData"))
}




#Effect sizes (Figure 3g)

num_individuals <- 300
num_steps <- 8
lag <- 2
effect_vector <- rep(0, num_steps)
size_list<- c(0.01, 0.02, 0.03, 0.04,0.05)
i <-2
  
mat1 <- matrix(0, num_iters,5)
mat2 <- matrix(0, num_iters,5)
mat3 <- matrix(0, num_iters,5)
mat4 <- matrix(0, num_iters,5)
for (j in 1:5){
  if (i==2){
    effect_vector[lag+1] = size_list[j]
  }
  cat("\n","effect_size=",size_list[j],"\n")
  out <- run_exps(num_individuals,num_steps,lag,effect_vector)
  mat1[,j] <- out[[1]]
  mat2[,j] <- out[[2]]
  mat3[,j] <- out[[3]]
  mat4[,j] <- out[[4]]
  cat("\n")
  cat("rejection_rate_1:",colMeans(mat1))
  cat(" rejection_rate_2:",colMeans(mat2))
  cat(" rejection_rate_3:",colMeans(mat3))
  cat(" rejection_rate_4:",colMeans(mat4))
}
save(mat1,mat2,mat3,mat4, file= paste0("results/sim_1_sizes_",i, ".RData"))

  
  
