
revenue_multi <- function(r, u, p_set,p_a, c){
  if (length(p_set)==1){
    return(r[p_set[1]] * u[p_set[1]] - (1-u[p_set[1]]) * p_a * c)
  }else{
    current <- p_set[1]
    p_set <- p_set[-1]
    return(r[current] * u[current] - (1-u[current])*(p_a * c - (1-p_a) * revenue_multi(r,u,p_set,p_a,c)))
  }
}

optimal_sequence <- function(r, u, q, c){
  p_a <- 1 - q
  value <- (r*u-c*p_a*(1-u))/(max(1-q*(1-u),0.001))
  incre_index <- sort(value,index.return=TRUE, decreasing=TRUE)$ix
  optimal <- c()
  for (i in (1:length(incre_index))){
    if (value[incre_index[i]]>0){
      optimal <- c(optimal, incre_index[i])
    }
  }
  return(optimal)
}

simu_choice <- function(r,u, sq, p_a){
  for (i in (1:length(sq))){
    if (runif(1)<u[sq[i]]){
      return(c(i))
    }else{
      if (runif(1)<p_a){
        return(c(-1,i))
      }
    }
  }
  return(c(0))
}

read_data <- function(){
  history = read.csv(file="./data/history_1665.csv", header=TRUE, sep=",")
  df = history
  # Ensure df is a data.frame
  df <- as.data.frame(df)
  
  # Ensure numeric types
  df$clk <- as.numeric(df$clk)
  df$nonclk <- as.numeric(df$nonclk)
  df$price <- as.numeric(df$price)
  
  
  # Step 1: Aggregate CTR and price by campaign_id
  aggregated_df <- df %>%
    group_by(campaign_id) %>%
    summarise(
      total_clk = sum(clk),
      CTR = total_clk / (total_clk + sum(nonclk)),
      price = first(price)
    )
  
  # Step 2: Take top 100 campaigns by total clicks
  top_product <- aggregated_df %>%
    arrange(desc(total_clk)) %>%
    slice_head(n = 100) %>%
    dplyr::select(campaign_id, CTR, price)
}


initialize_env <- function(){
  n_product <<- 100
  c_cost <<- 6
  current_curve <<-0
  average_regret <<- 0
  T_total <<- 5000
  epoch_num <<- c(1:T_total)
  data_plot <<- c()
  average_estimation <<- 0
  data_plot <<- c()
  epsilon <<- 0.05
}

set.seed(1)
initialize_env()
top_product = read_data()
prod_feature <- top_product$price
r <- rep(1,n_product)
u <- rep(1,n_product)


alpha <- c(-4.8265064,0.1221894 ,0,  0.5344133  ,-0.1386934 ,-0.2627397 ,
           -0.3177119 ,-0.3058237 ,-0.2336395  ,0,  0.5232575 ,  0.3031273  ,0)
beta = c(-2.7655468143,-0.0001917805)



data_plot <- c()


current_curve <-0
average_regret <- 0
average_estimation <- 0
while (current_curve<curve_num){
  current_curve <- current_curve + 1
  #print(current_curve)
  data_abandon <-  c()
  abandon_feature <- c()
  data_prod <- c()
  choice_feature <- c()
  alpha_hat <- rep(1,13)
  beta_hat <-  c(1,1)
  uhat <- rep(1,n_product)
  qhat <- 1
  tilde_M <- diag(13)
  M_matrix <- diag(2)
  
  regret_total <- 0
  regret <- c()
  t <- 0
  explore_num <- 0
  while (t<T_total){
    t <- t + 1
    feature_cus <- c(1) 
    ogender <- sample(0:1,1)
    if (ogender==0){
      feature_cus <- c(feature_cus,1,0)
    }else{
      feature_cus <- c(feature_cus,0,1)
    }
    oage <- sample(0:6,1)
    feature_cus <- c(feature_cus,0,0,0,0,0,0,0)
    feature_cus[oage+4] <- feature_cus[oage+4] + 1
    oshopping <- sample(1:3,1)
    feature_cus <- c(feature_cus,0,0,0)
    feature_cus[oshopping + 10] <- feature_cus[oshopping + 10] + 1
    theta <- sum(alpha_hat*c(feature_cus))
    qhat <- exp(theta)/(1+exp(theta))
    theta <- sum(alpha*feature_cus)
    q <- exp(theta)/(1+exp(theta))
    for (i in 1:n_product){
      theta <- sum(beta*c(1,prod_feature[i]))
      u[i] <- exp(theta)/(1 + exp(theta))
      feature_prod <- c(1,prod_feature[i])
      theta <- sum(beta_hat * feature_prod)
      if (theta>10){
        uhat[i]=1
      }else{
        uhat[i] <- exp(theta)/(1+exp(theta))
      }
      
    }
    p_a = q
    q= 1-q
    optimal_sq <- optimal_sequence(r,u,q,c_cost)
    if (length(optimal_sq)>0){
      revenue_op <- revenue_multi(r, u, optimal_sq,p_a, c_cost)
    }else{
      revenue_op <- 0
    }
    
    sq <- optimal_sequence(r,uhat,qhat,c_cost)
    
    if (runif(1) < epsilon || is.null(sq)) {
      eligible_indices <- setdiff(1:n_product, sq)
      if (length(sq) < n_product){
        selected_index <- sample(eligible_indices, 1)
        sq <- c(selected_index,sq)
      }
    }
    revenue_sq <- revenue_multi(r, u, sq,p_a, c_cost)
    regret_total <- regret_total + revenue_op - revenue_sq
    regret <- c(regret, regret_total)
    choice <- simu_choice(r,u, sq, p_a)
    
    if (length(choice)==1){
      if (choice==0){
        for (i in 1:length(sq)){
          feature_prod <- c(1,prod_feature[sq[i]])
          data_abandon <- c(data_abandon,1)
          data_prod <- c(data_prod,0)
          abandon_feature <- rbind(abandon_feature,feature_cus)
          choice_feature <- rbind(choice_feature,feature_prod)
          tilde_M <- tilde_M + feature_cus %*% t(feature_cus)
          M_matrix <- M_matrix + feature_prod %*% t(feature_prod)
        }
      }else{
        for (i in (1:choice)){
          if (i<choice){
            data_abandon <- c(data_abandon,1)
            data_prod <- c(data_prod,0)
            abandon_feature <- rbind(abandon_feature,feature_cus)
            choice_feature <- rbind(choice_feature,c(1,prod_feature[sq[i]]))
          }
          if (i==choice){
            data_prod <- c(data_prod,1)
            choice_feature <- rbind(choice_feature,c(1,prod_feature[sq[i]]))
          }
          feature_prod <- c(1,prod_feature[sq[i]])
          M_matrix <- M_matrix + feature_prod %*% t(feature_prod)
          tilde_M <- tilde_M + feature_cus %*% t(feature_cus)
        }
      }
    }else{
      level <- choice[2]
      if (level>1){
        for (i in (1:(level-1))){
          feature_prod <- c(1,prod_feature[sq[i]])
          data_abandon <- c(data_abandon,1)
          data_prod <- c(data_prod,0)
          abandon_feature <- rbind(abandon_feature,feature_cus)
          choice_feature <- rbind(choice_feature,feature_prod)
          tilde_M <- tilde_M + feature_cus %*% t(feature_cus)
          M_matrix <- M_matrix + feature_prod %*% t(feature_prod)
        }
      }
      data_abandon <- c(data_abandon,0)
      data_prod <- c(data_prod,0)
      abandon_feature <- rbind(abandon_feature,feature_cus)
      tilde_M <- tilde_M + feature_cus %*% t(feature_cus)
      feature_prod <- c(1,prod_feature[sq[level]])
      choice_feature <- rbind(choice_feature,feature_prod)
      M_matrix <- M_matrix + feature_prod %*% t(feature_prod)
    }
    
    
    my_data <- as.data.frame(cbind(data_prod, choice_feature))
    colnames(my_data) <- c("pick","f1","f2")
    logit <- glm(formula= pick~ f2,
                 data = my_data,
                 family = binomial)
    if (nrow(my_data)>20){
      for (j in 1:2){
        beta_hat[j] <- logit$coefficients[j]
      }
    }
    
    my_data <- as.data.frame(cbind(data_abandon,abandon_feature))
    if (nrow(my_data)>500){
      colnames(my_data) <- c("pick","f1","f2","f3","f4","f5",
                             "f6","f7","f8","f9","f10",
                             "f11","f12","f13")
      logit <- glm(formula= pick~ -1 + f1+f2+f4+f5+f6+f7+f8+f9+f11+f12,
                   data = my_data,
                   family = binomial)
      for (j in 1:2){
        alpha_hat[j] <- logit$coefficients[j]
      }
      for (j in 4:9){
        alpha_hat[j] <- logit$coefficients[j-1]
      }
      for (j in 11:12){
        alpha_hat[j] <- logit$coefficients[j-2]
      }
      alpha_hat[3]=0
      alpha_hat[10] = 0
      alpha_hat[13] = 0
    }
  }
  
  average_regret <- average_regret + regret_total
  average_estimation <- average_estimation + sum(abs(uhat-u))+abs(qhat-q)
  epoch_num <- c(1:T_total)
  dd <- cbind(epoch_num,regret)
  epoch_num <- rep(current_curve,T_total)
  dd <- cbind(dd,epoch_num)
  data_plot <- rbind(data_plot,dd)
  
  
}
