#rm(list=ls())

library(spls)
library(data.table)

train <- read.csv("trainselct_iso.csv",header=FALSE)  #load the indicators of training data to make sure the separate is the same as matlab code
YX <- fread("isolet.txt")
Y <-unlist(as.vector(YX[,1]))  
X <- YX[, -1]           # Exclude the first column (Y)


X <- scale(X)
YX <- cbind(Y,X)

# Set the number of repetitions and initialize storage
numRepeats <- 50
nnz_spls_values <- numeric(numRepeats)  # To store the number of selected features
misclassification_rates <- numeric(numRepeats)  # To store misclassification rates
B_best_matrix <- matrix(0, nrow = numRepeats, ncol = ncol(X) + 1)  # Matrix to store B_best coefficients
time <- numeric(numRepeats)

# Repeat the process 50 times
for (i in 1:numRepeats) {
  rst<- Sys.time()
  trainIDX <- which(train[,i]==1)
  Y_train <- Y[trainIDX]
  X_train <- X[trainIDX,]
  Y_test <- Y[-trainIDX]
  X_test <- X[-trainIDX,]
  
  # Fit SGPLS model using cross-validation
  
  # Fit cv.sgpls for each class
  cv_fit <- cv.sgpls(X_train,Y_train, fold=5, K=c(1:10), 
              eta=seq(0.1,0.95,0.05),plot.it=FALSE)  # 5-fold cross-validation
  
  # Get the best lambda (the one that minimizes cross-validation error)
  model.fit <- sgpls(X_train, Y_train, K=cv_fit$K.opt,eta=cv_fit$eta.opt)
  
  # Store the coefficients (including intercept) in the matrix
  B_best_matrix[i, ] <- model.fit$betahat
  
  # Make predictions on the test set
  pred <- predict.sgpls(model.fit, X_test,type="fit")
  
  # Compute misclassification rate and the number of non-zero coefficients
  nnz_spls_values[i] <- sum(model.fit$betahat[-1] != 0)   # Excluding intercept
  misclassification_rates[i] <- mean(pred != as.numeric(Y_test))  # Adjust factor levels

  rend<- Sys.time()
  time[i]<- rend-rst

  # Print progress
  cat(sprintf("Iteration %d: nnz = %d, Misclassification Rate = %.4f\n",
              i, nnz_spls_values[i], misclassification_rates[i]))
}

# Create a data frame with results (excluding B_best)
results_df <- data.frame(Iteration = 1:numRepeats, 
                         NNZ_SGPLS = nnz_spls_values, 
                         Misclassification_Rate = misclassification_rates,
                         Time = time)

# Save results to CSV
write.csv(results_df, "sgpls_results_R_iso.csv", row.names = FALSE)

# Save the B_best coefficients (if needed)
write.csv(B_best_matrix, "B_best_sgpls_R_iso.csv", row.names = FALSE)

# Print final summary
cat(sprintf("Results saved to sgpls_results_R_iso.csv and B_best_sgpls_R_iso.csv\n"))
cat(sprintf("Median NNZ: %.2f\n", median(nnz_spls_values)))
cat(sprintf("Median Misclassification Rate: %.4f\n", median(misclassification_rates)))
