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_gpocre_values = zeros(numRepeats,1); % Store selected features count
misclassification_gpocre_rates = zeros(numRepeats,1); % Store misclassification rates
numFeatures_gpocre = size(X, 2); % Number of features in X
B_best_matrix = zeros(numRepeats, numFeatures_gpocre); % Store B_best for each iteration (including intercept)
trainselct = zeros(size(YX,1),numRepeats);
index = zeros(numRepeats,1);
time = zeros(numRepeats,1)
%X = log2(X); % Log transform
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);
    % Create logical indices for training and testing sets
    trainIdx = training(cv); % Logical index for training set
    testIdx = test(cv); % Logical index for test set
    Y_train = Y(trainIdx);
    Y_test = Y(testIdx);
    X_train = X(trainIdx,:);
    X_test = X(testIdx,:);
    %[X_train,idxX_train] = Standardize(YX_train(:,5:end));    % Standardizing all covariates
    trainselct(:,i) = trainIdx;
    % Settings:
    option.delta = 0.01;        % 0.1,0.05,0.02, or 0.01
    option.ptype = 'ebtz';
    nCVFolds = 5;

    tp = CVGPOCRE(Y_train,X_train,'binomial','logit',nCVFolds,option);
    disp(['  gPOCRE: Optimal Tuning Parameter = ' num2str(tp)]);
    tmpRes = gPOCRE(Y_train,X_train,'binomial','logit',tp,option);
    option.W = CalcOptW(tmpRes.mu+X_train*tmpRes.beta,'binomial','logit');
    [retRes,tmpStats] = gPOCRE(Y_train,X_train,'binomial','logit',tp,option);


    testYX = [Y_test,X_test];
    tmpEta = X_test*retRes.beta+retRes.mu;
    tmpEY = ginv(tmpEta,'logit');
    yPred = round(tmpEY);
    mr = sum(abs(Y_test-yPred))/size(X_test,1);
    nnz_gpocre_values(i) = nnz(retRes.beta); % Number of selected features (excluding intercept)
    misclassification_gpocre_rates(i) = mr;
    B_best_matrix(i, :) = retRes.beta;
    % Display progress
    fprintf('Iteration %d: nnz = %d, Misclassification Rate = %.4f\n', ...
            i, nnz_gpocre_values(i), misclassification_gpocre_rates(i));
    index(i) = i;
    time(i)= toc;
end
% Combine results into a table
results_table_gpocre = table(index, nnz_gpocre_values, misclassification_gpocre_rates, time,...
                      'VariableNames', {'Iteration', 'NNZ_gpocre', 'Misclassification_Rate','time'});

% Save results as a CSV file
csv_filename = 'isolet\gpocre_results_iso.csv';
writetable(results_table_gpocre, csv_filename);
csvwrite('isolet\B_best_matrix_gpocre_iso.csv', B_best_matrix);
csvwrite('isolet\trainselct_gpocre_iso.csv', trainselct);

median(nnz_gpocre_values)
mean(nnz_gpocre_values)
std(nnz_gpocre_values)/sqrt(50)
median(misclassification_gpocre_rates)
mean(misclassification_gpocre_rates)
std(misclassification_gpocre_rates)/sqrt(50)
median(time)
mean(time)
std(time)/sqrt(50)