clear
pwd

YX = load('breast\breast.txt');
Y_multi = categorical(YX(:,4));
[X,idxX] = Standardize(YX(:,5:end));

% Define number of repetitions
numRepeats = 55;
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(2*numRepeats, numFeatures_gpocre); % Store B_best for each iteration (including intercept)
trainselct = zeros(size(YX,1),numRepeats);
index = zeros(numRepeats,1);
%X = log2(X); % Log transform
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);
    % 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 = YX(trainIdx,1:3)
    Y_test = YX(testIdx,1:3);
    X_train = X(trainIdx,:);
    X_test = X(testIdx,:);
    %[X_train,idxX_train] = Standardize(YX_train(:,5:end));    % Standardizing all covariates

    % Settings:
    option.delta = 0.01;        % 0.1,0.05,0.02, or 0.01
    %option.covidx = [15,16];    % covariates of no interest
    option.ptype = 'ebtz';
    nCVFolds = 5;

    tp = CVGPOCRE(Y_train,X_train,'multibc','bclogit',nCVFolds,option);
    disp(['  gPOCRE: Optimal Tuning Parameter = ' num2str(tp)]);
    tmpRes = gPOCRE(Y_train,X_train,'multibc','bclogit',tp,option);
    option.W = CalcOptW(tmpRes.mu+X_train*tmpRes.beta,'multibc','bclogit');
    [retRes,tmpStats] = gPOCRE(Y_train,X_train,'multibc','bclogit',tp,option);
    if retRes.nCmp==0
        continue
    else


    retRes.pvalue = tmpStats.pvalue;
    Y_bc_test = Y_test(:,2:end);
    testYX = [Y_bc_test,X_test];
    tmpEta = X_test*retRes.beta+retRes.mu;
    tmpEY = arrayfun(@(i)(ginv(tmpEta(i,:),'bclogit')),...
                         [1:size(tmpEta,1)]','UniformOutput',false);
    tmpEY = cell2mat(tmpEY);
    tmpEY = [1-sum(tmpEY,2) tmpEY];
    yPred = bsxfun(@eq,tmpEY,max(tmpEY,[],2));
    mr = sum(sum(abs(Y_test-yPred)))/(2*size(Y_test,1));
    nnz_gpocre_values(i) = nnz(retRes.beta)/(size(Y_test,2)-1); % Number of selected features (excluding intercept)
    misclassification_gpocre_rates(i) = mr;
    B_best_matrix(2*i-1:2*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
    end
end

nnz_gpocre_values = nnz_gpocre_values([2:47,49:52]);   % iteration 1 and 48 don't have any results
misclassification_gpocre_rates = misclassification_gpocre_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_gpocre = table(index, nnz_gpocre_values, misclassification_gpocre_rates, ...
                      'VariableNames', {'Iteration', 'NNZ_gpocre', 'Misclassification_Rate'});

% Save results as a CSV file
csv_filename = 'breast\gpocre_results_br.csv';
writetable(results_table_gpocre, csv_filename);
csvwrite('breast\B_best_matrix_gpocre_br.csv', B_best_matrix);
csvwrite('breast\trainselct_gpocre.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)