clear
pwd

YX = readtable('isolet\isolet.txt');
YX = table2array(YX);
Y = YX(:,1);
[X,idxX] = Standardize(YX(:,2:end));

% Define number of repetitions
numRepeats = 50;
nnz_lasso_values = zeros(numRepeats,1); % Store selected features count
misclassification_rates = zeros(numRepeats,1); % Store misclassification rates
numFeatures = size(X, 2); % Number of features in X
B_best_matrix = zeros(numRepeats, numFeatures + 1); % Store B_best for each iteration (including intercept)
trainselct = zeros(size(YX,1),numRepeats); 
time = zeros(numRepeats,1);
% Repeat the process 50 times with a fixed seed for reproducibility
for i = 1:numRepeats
    tic
    rng(i+50); % Set seed to ensure reproducibility for each iteration

    % Stratified partitioning
    trainRatio = 0.8;
    cv = cvpartition(Y, 'HoldOut', 1 - trainRatio);
    trainIdx = training(cv);
    testIdx = test(cv);

    % Separate training and testing data
    Y_train = Y(trainIdx,:);
    Y_test = Y(testIdx,:);
    X_train = X(trainIdx,:);
    X_test = X(testIdx,:);
    trainselct(:,i) = trainIdx;
    % Perform LASSO logistic regression
    [B, FitInfo] = lassoglm(X_train, Y_train, 'binomial', 'CV', 5,'MaxIter', 1e4);

    % Select best model based on minimum cross-validation error
    bestLambda = FitInfo.LambdaMinDeviance;
    B_best = B(:, FitInfo.IndexMinDeviance);
    Intercept = FitInfo.Intercept(FitInfo.IndexMinDeviance);
    B_best = [Intercept; B_best];
    B_best_matrix(i, :) = B_best;

    % Predict probabilities
    %X_test = [ones(size(X_test, 1), 1) X_test]; % Add intercept column
    probabilities = glmval(B_best, X_test, 'logit');

    % Convert probabilities to binary predictions
    predicted_labels = probabilities > 0.5;

    % Compute misclassification rate and number of selected features
    nnz_lasso_values(i) = nnz(B_best(2:end)); % Number of selected features (excluding intercept)
    misclassification_rates(i) = mean(predicted_labels ~= Y_test);
    time(i)=toc;
    % Display progress
    fprintf('Iteration %d: nnz = %d, Misclassification Rate = %.4f\n', ...
            i, nnz_lasso_values(i), misclassification_rates(i));
end

% Combine results into a table
results_table = table((1:numRepeats)', nnz_lasso_values, misclassification_rates,time, ...
                      'VariableNames', {'Iteration', 'NNZ_Lasso', 'Misclassification_Rate','time'});

% Save results as a CSV file
csv_filename = 'isolet\lasso_results_iso.csv';
writetable(results_table, csv_filename);
csvwrite('isolet\B_best_matrix_lasso_iso.csv', B_best_matrix);

% summary statistics for non-zero coefficient, misclassification rates, running time 
median(nnz_lasso_values) 
mean(nnz_lasso_values)
std(nnz_lasso_values)/sqrt(50)
median(misclassification_rates)
mean(misclassification_rates)
std(misclassification_rates)/sqrt(50)
median(time)
mean(time)
std(time)/sqrt(50)