#install.packages('R.matlab')
#install.packages("e1071")

rm(list=ls())
library(readxl) # read_excel
library(R.matlab) # readMat
library(stats) # p.adjust
library(e1071) # SVM


### 0. Data Loading and Preprocessing

# Set working directory to the folder containing the data
setwd("c:/Users/koritakova/Documents/Vyuka/Vicerozmerky a AKD/Analyza a klasifikace dat/2018-podzim/AKD_predn-07_cviceni/Data_klasifikace/")

# Load data from Excel sheet
data <- read_excel("gm_identifikatory.xls", sheet = "List1")
names <- data$id  # File names of the images
group <- data$group  # Group assignment (patients vs controls)

inxP <- which(group == 1)  # Indices for patients
inxC <- which(group == 0)  # Indices for controls

# Load brain mask (GM mask)
mask <- readMat("gm_mask_2D.mat")[[1]]
image(mask, main = "GM Mask") # Display the mask

# Load the first patient's image as an example
i <- 1
img <- readMat(paste0(names[inxP[i]], "_rez50.mat"))[[1]]
image(img, main = paste("Patient", i)) # Display the 2D image of the first patient's GM data

# Load all image data into a matrix
X <- matrix(0, nrow = length(inxP) + length(inxC), ncol = sum(mask == 1)) # Create matrix to hold the image data for all subjects (patients and controls)
for (i in 1:length(inxP)) {
  img <- readMat(paste0(names[inxP[i]], "_rez50.mat"))[[1]]
  X[i, ] <- img[mask == 1]
}
for (i in 1:length(inxC)) {
  img <- readMat(paste0(names[inxC[i]], "_rez50.mat"))[[1]]
  X[i + length(inxP), ] <- img[mask == 1]
}


### 1. Data Reduction and Classification: Resubstitution Validation

# Perform Mann-Whitney test for feature selection
pvals <- sapply(1:ncol(X), function(j) wilcox.test(X[group == 1, j], X[group == 0, j], exact=FALSE)$p.value)
sum(pvals<0.05) # 1457 significant pixels before correction

# Apply False Discovery Rate (FDR) correction
pvals_cor <- p.adjust(pvals, method = "fdr")
sum(pvals_cor<0.05) # 226 significant pixels after FDR correction
X_red <- X[, pvals_cor < 0.05] # Select features with corrected p-value < 0.05

# Display the significant pixels
significant_mask <- mask
significant_mask[mask == 1] <- 1 + (pvals_cor < 0.05)
image(significant_mask, main = "Significant Pixels")

# Train SVM classifier with the reduced features (X_red)
svm_model <- svm(X_red, as.factor(group))  # SVM training with C = 1
group_klasif <- predict(svm_model, X_red)  # SVM prediction
table(group_klasif,group) # Check accuracy
accuracy <- mean(group_klasif == group)
sensitivity <- mean(group_klasif[group == 1] == group[group == 1])
specificity <- mean(group_klasif[group == 0] == group[group == 0])
cat("Accuracy:", accuracy, "\n")
cat("Sensitivity:", sensitivity, "\n")
cat("Specificity:", specificity, "\n")
# 93.8%, 91.7%, 95.8%


# Train SVM model with BoxConstraint (C) parameter set to 100
svm_model_C100 <- svm(X_red, as.factor(group), cost = 100)
group_klasif_C100 <- predict(svm_model_C100, X_red)
accuracy_C100 <- mean(group_klasif_C100 == group)
sensitivity_C100 <- mean(group_klasif_C100[group == 1] == group[group == 1])
specificity_C100 <- mean(group_klasif_C100[group == 0] == group[group == 0])
cat("Accuracy with C=100:", accuracy_C100, "\n")
cat("Sensitivity with C=100:", sensitivity_C100, "\n")
cat("Specificity with C=100:", specificity_C100, "\n")
# 100%, 100%, 100%

# Train SVM model with BoxConstraint (C) parameter set to 0.1
svm_model_C0.1 <- svm(X_red, as.factor(group), cost = 0.1)
group_klasif_C0.1 <- predict(svm_model_C0.1, X_red)
accuracy_C0.1 <- mean(group_klasif_C0.1 == group)
sensitivity_C0.1 <- mean(group_klasif_C0.1[group == 1] == group[group == 1])
specificity_C0.1 <- mean(group_klasif_C0.1[group == 0] == group[group == 0])
cat("Accuracy with C=0.1:", accuracy_C0.1, "\n")
cat("Sensitivity with C=0.1:", sensitivity_C0.1, "\n")
cat("Specificity with C=0.1:", specificity_C0.1, "\n")
# 85.4%, 79.2%, 91.7%



### 2. Data Reduction and Classification: Leave-One-Out Cross-Validation (Incorrect Variant)

# Perform Mann-Whitney test for feature selection
pvals <- sapply(1:ncol(X), function(j) wilcox.test(X[group == 1, j], X[group == 0, j], exact=FALSE)$p.value)
sum(pvals<0.05) # 1457 significant pixels before correction

# Apply False Discovery Rate (FDR) correction
pvals_cor <- p.adjust(pvals, method = "fdr")
sum(pvals_cor<0.05) # 226 significant pixels before correction
X_red <- X[, pvals_cor < 0.05] # Select features with corrected p-value < 0.05

# Leave-One-Out Cross-Validation for feature selection and SVM
group_klasif2 <- numeric(length(group))

for (i in 1:length(group)) {
  cat("Processing subject", i, "\n")
  x_test_red <- X_red[i, ] # test data (1 subject)
  X_train_red <- X_red[-i, ] # train data
  group_train <- group[-i]
  
  svm_model_loo <- svm(X_train_red, as.factor(group_train)) # train model
  group_klasif2[i] <- predict(svm_model_loo, t(as.matrix(x_test_red)))
}

table(group_klasif2,group)
group_klasif2 <- as.numeric(group_klasif2)-1
table(group_klasif2,group)
accuracy_loo <- mean(group_klasif2 == group)
sensitivity_loo <- mean(group_klasif2[group == 1] == group[group == 1])
specificity_loo <- mean(group_klasif2[group == 0] == group[group == 0])
cat("Accuracy (LOO):", accuracy_loo, "\n")
cat("Sensitivity (LOO):", sensitivity_loo, "\n")
cat("Specificity (LOO):", specificity_loo, "\n")
# 79.2%; 79.2%; 79.2%


### 3. Data Reduction and Classification: Leave-One-Out Cross-Validation (Incorrect Variant)

# UKOL - vypocitejte spravnost, senzitivitu a specificitu pri pouziti
#        spravne krizove validace
#      - zamyslete se nad tim, jakym zpusobem by bylo mozne vizualizovat 
#        statisticky vyznamne pixely pri pouziti spravne krizove validace