clear
pwd
% Load data
YX = load('breast\breast.txt');
Y_multi = categorical(YX(:,4));
X = YX(:,5:end);
Y = YX(:,2:3);  % first class as baseline

% Standardization
[X, idxX] = Standardize(X); 

% Define number of repetitions
numRepeats = 55;
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(2*numRepeats, numFeatures + 1); % Store B_best for each iteration (including intercept)
trainselct = zeros(size(YX,1),numRepeats);
% Repeat the process 50 times with a fixed seed for reproducibility
for i = 1:numRepeats
    rng(i); % Set seed to ensure reproducibility for each iteration

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

    % Separate training and testing data
    Y_train = Y(trainIdx,:);
    Y_test = Y(testIdx,:);
    Y_multi_test = Y_multi(testIdx,:);
    X_train = X(trainIdx,:);
    X_test = X(testIdx,:);
    trainselct(:,i) = trainIdx;
    B = cell(2, 1);
    FitInfo = cell(2,1);
    bestLambda = cell(2,1);
    B_best = zeros(size(X_train,2),2);
    B_best_inter = zeros(size(X_train,2)+1,2);
    predictions = zeros(size(X_test, 1), 3);
    % Perform LASSO logistic regression
    for j = 1:2
        [B{j},FitInfo{j}] = lassoglm(X_train, Y_train(:,j), 'binomial', 'CV', 5);
        bestLambda{j} = FitInfo{j}.LambdaMinDeviance;
        B_best(:,j) = B{j}(:, FitInfo{j}.IndexMinDeviance);

        Intercept = FitInfo{j}.Intercept(FitInfo{j}.IndexMinDeviance);
        B_best_inter(:,j) = [Intercept; B_best(:,j)];
        predictions(:, j+1) = glmval(B_best_inter(:,j), X_test, 'logit');
    end
    predictions(:,1) = 1-predictions(:,2)-predictions(:,3);
    [~, ypred] = max(predictions, [], 2);
    mr = sum(Y_multi_test~=categorical(ypred))/size(Y_test,1);
    
    % Compute misclassification rate and number of selected features
    nnz_lasso_values(i) = nnz(abs(B_best(:,1))+abs(B_best(:,2))); % Number of selected features (excluding intercept)
    misclassification_rates(i) = mr;
    B_best_matrix(2*i-1:2*i, :) = B_best_inter';
    % Display progress
    fprintf('Iteration %d: nnz = %d, Misclassification Rate = %.4f\n', ...
            i, nnz_lasso_values(i), misclassification_rates(i));
end

nnz_lasso_values = nnz_lasso_values([2:47,49:52]);  % keep the same as gpocre
misclassification_rates = misclassification_rates([2:47,49:52]);
B_best_matrix = B_best_matrix([3:94,97:104],:);
index = index([2:47,49:52]);

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

% Save results as a CSV file
csv_filename = 'breast\lasso_results_br.csv';
writetable(results_table, csv_filename);
csvwrite('breast\B_best_matrix_lasso_br.csv', B_best_matrix);



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)