library(keras)
library(tensorflow)
library(fst)
library(data.table)
library(progress)
library(depstats)
library(foreach)

sample_sizes <- c(30, 50, 100, 200, 300, 400)

image_model <- load_model_hdf5("../1_Models/image.h5")
score_model <- load_model_hdf5("../1_Models/score.h5")
combined_model <- load_model_tf("../1_Models/score_image.h5")


image_thresholds <- readRDS('../2_Threshold/threshold/image_thresholds.RDS')
score_thresholds <- readRDS('../2_Threshold/threshold/score_thresholds.RDS')
combined_thresholds <- readRDS('../2_Threshold/threshold/combined_thresholds.RDS')
other_thresholds <- readRDS('../2_Threshold/threshold/other_thresholds.RDS')

# —————————————————————————————————————————————
# define exactly the same slicing indices as in your first script
thr_cols  <- 19
pix_count <- 25*25*1       # = 625
start_img <- thr_cols + 2  # skip the 19 thr cols + the “n” col
end_img   <- start_img + pix_count - 1
score_dim <- as.integer(score_model$input_shape[[2]])
# —————————————————————————————————————————————

model_predict <- function(model, which, d) {
  if (which == "image") {
    # pull columns 21:645 and reshape into n × 25 × 25 × 1
    x_img <- array_reshape(
      as.matrix(d[, start_img:end_img]),
      c(nrow(d), 25L, 25L, 1L)
    )
    preds <- model$predict(x_img)

  } else if (which == "score") {
    # pull columns 1:20
    x_score <- as.matrix(d[, 1:score_dim])
    preds   <- model$predict(x_score)

  } else if (which == "combined") {
    # image branch
    x_img <- array_reshape(
      as.matrix(d[, start_img:end_img]),
      c(nrow(d), 25L, 25L, 1L)
    )
    # MLP branch on the first 19 thr columns
    x_comb_score <- as.matrix(d[, 1:thr_cols])
    preds         <- model$predict(list(x_img, x_comb_score))

  } else {
    stop("Unknown model type: ", which)
  }

  as.numeric(preds)
}


##########
#dependent
# "Linear", "Diamond", "Triangle", "Crescent", "Points", "Expo-nential",
# "Circles", "Cross", "Wedge", "Cubic", "W-shape", "Parabola", "Two-parabola",
# "Sine", "Doppler", "Heavy-sine", "Heart", "Spiral", "Taegeuk", "Sam-taegeuk"

pb <- progress_bar$new(total = length(sample_sizes) * 3 * 20)

powers <- foreach(n = sample_sizes, .combine = 'rbind') %do% {
  df <- read.fst(sprintf("../3_Testing/dependent_tests/%s/deptest.fst", n))

  foreach(k = 1:20, .combine = 'rbind') %do% {

    d <- df[(1 + 1000 * (k - 1)):(1000 * k),]

    image_power <- mean(model_predict(image_model, 'image', d) > image_thresholds[sample_size == n, thres])
    pb$tick()
    score_power <- mean(model_predict(score_model, 'score', d) > score_thresholds[sample_size == n, thres])
    pb$tick()
    combined_power <- mean(model_predict(combined_model, 'combined', d) > combined_thresholds[sample_size == n, thres])
    pb$tick()

    otherpower <- foreach(i = 1:19, .combine = 'cbind') %do% {
      mean(d[,i] >= other_thresholds[sample_size == n, get(paste0('score',i))])
    }

    c(combined_power, score_power, image_power, otherpower)
  }
}

powers <- data.table(powers)

colnames(powers) <- c('combined', 'score', 'image', paste0('score', 1:19))

powers[, `:=`(sample_size = rep(sample_sizes, each = 20),
              scenario = rep(1:20, 6))]

powers %>% saveRDS('dep_powers.RDS')

dep_powers <- powers

#data written into the package
#usethis::use_data(dep_powers, overwrite=TRUE)

##########
#additional
# 'Correlated Laplace A', 'Correlated Laplace B', 'Ishigami style A', 'Ishigami style  B', 'Tree ring A', 'Tree ring B', 'Escalating variance A', 'Escalating variance B',
# 'Infinity A', 'Infinity B', 'Pi A', 'Pi B'

pb <- progress_bar$new(total = length(sample_sizes) * 3 * 12)

powers <- foreach(n = sample_sizes, .combine = 'rbind') %do% {
  df <- cbind(read.fst(sprintf("../3_Testing/additional_tests/%s/add.fst", n)), n,
              read.fst(sprintf("../3_Testing/additional_tests/%s/add_image.fst", n)))

  foreach(k = 1:12, .combine = 'rbind') %do% {

    d <- df[(1 + 1000 * (k - 1)):(1000 * k),]

    image_power <- mean(model_predict(image_model, 'image', d) > image_thresholds[sample_size == n, thres])
    pb$tick()
    score_power <- mean(model_predict(score_model, 'score', d) > score_thresholds[sample_size == n, thres])
    pb$tick()
    combined_power <- mean(model_predict(combined_model, 'combined', d) > combined_thresholds[sample_size == n, thres])
    pb$tick()

    otherpower <- foreach(i = 1:19, .combine = 'cbind') %do% {
      mean(d[,i] >= other_thresholds[sample_size == n, get(paste0('score',i))])
    }

    c(combined_power, score_power, image_power, otherpower)
  }
}

powers <- data.table(powers)

colnames(powers) <- c('combined', 'score', 'image', paste0('score', 1:19))

powers[, `:=`(sample_size = rep(sample_sizes, each = 12),
              scenario = rep(1:12, 6))]

powers %>% saveRDS('add_powers.RDS')

add_powers <- powers

#data written into the package
#usethis::use_data(add_powers, overwrite=TRUE)

##########
#increasing noise
#Circles Cross Sine Spiral
#4 noise levels

pb <- progress_bar$new(total = length(sample_sizes) * 3 * 16)

powers <- foreach(n = sample_sizes, .combine = 'rbind') %do% {
  df <- cbind(read.fst(sprintf("../3_Testing/noise_tests/%s/noise.fst", n)), n,
              read.fst(sprintf("../3_Testing/noise_tests/%s/noise_image.fst", n)))

  foreach(k = 1:16, .combine = 'rbind') %do% {

    d <- df[(1 + 1000 * (k - 1)):(1000 * k),]

    image_power <- mean(model_predict(image_model, 'image', d) > image_thresholds[sample_size == n, thres])
    pb$tick()
    score_power <- mean(model_predict(score_model, 'score', d) > score_thresholds[sample_size == n, thres])
    pb$tick()
    combined_power <- mean(model_predict(combined_model, 'combined', d) > combined_thresholds[sample_size == n, thres])
    pb$tick()

    otherpower <- foreach(i = 1:19, .combine = 'cbind') %do% {
      mean(d[,i] >= other_thresholds[sample_size == n, get(paste0('score',i))])
    }

    c(combined_power, score_power, image_power, otherpower)
  }
}

powers <- data.table(powers)

colnames(powers) <- c('combined', 'score', 'image', paste0('score', 1:19))

powers[, `:=`(sample_size = rep(sample_sizes, each = 16),
              scenario = rep(1:16, 6))]

powers %>% saveRDS('increasingnoise_powers.RDS')

increasingnoise_powers <- powers

#data written into the package
#usethis::use_data(increasingnoise_powers, overwrite=TRUE)