 clearvars -except stock_RFE_num e RFE_num file_name f;

load YOUR_DATA.mat

RFE_num= ; %Define the recursive feature elimination proportion
K_fold = 10; % 10-fold cross-vaidation
w = [];%save weight
deci_value=[];%save decision values
predicted_label_all=[];
label_nCV= [];
CrossValIndices_out = crossvalind('Kfold', label,K_fold);
h = waitbar(0,'Please wait..');
bag_train_data= [];
bag_train_label= [];
bag_perc = 1;%how many percentages of the data are going to be incluede in the reampled samples, 0-1
bag_num = 1000; %define number of resample
for i = 1:K_fold
    bag = 1;
    model_predicted_label=[];
    model_deci=[];
    %seperate train and test data
    new_DATA = data_all;
    new_label  = label;
    train_label = new_label(CrossValIndices_out~=i,:);train_data = new_DATA(CrossValIndices_out~=i,:);
    test_label = new_label(CrossValIndices_out==i,:);test_data = new_DATA(CrossValIndices_out==i,:);
    waitbar(i/K_fold,h,[num2str(i),'/',num2str(K_fold)]) %waitbar
    %data regularization
    [train_data,PS] = mapminmax(train_data',0,1);
    test_data = mapminmax('apply',test_data',PS);
    train_data = train_data';
    test_data = test_data';
    %%%feature selection
    step = 1;% how many features are eliminated in one step
    ftRank = SVMRFE(train_label,train_data, step,'-t 0');
    IX = ftRank(1:ceil(length(ftRank)*RFE_num)); %??????1%
    while bag <bag_num+1  
        %resample data
        length(train_label);
        %bag_index_non=truerand(length(find(train_label==-1))*bag_perc,min(find(train_label==-1)),max(find(train_label==-1)));
        %bag_index_ap=truerand(length(find(train_label==1))*bag_perc,min(find(train_label==1)),max(find(train_label==1)));
        bag_index_non = randsample([min(find(train_label==-1)):max(find(train_label==-1))],length(find(train_label==-1))*bag_perc,true);
        bag_index_ap = randsample([min(find(train_label==1)):max(find(train_label==1))],length(find(train_label==1))*bag_perc,true);
        
        
        bag_train_data = [train_data(bag_index_non,:);train_data(bag_index_ap,:)];
        bag_train_label = [train_label(bag_index_non,:);train_label(bag_index_ap,:)];
    
        % achive best C from inner "inner loop"
        [bestacc,bestc] = SVMcgForClass_NoDisplay_linear_nc(bag_train_label,bag_train_data(:,IX),-10,10,5,0.1);
        cmd = ['-t 0 ', ' -c ',num2str(bestc),' -w1 2 -w-1 1'];
        %train the model
        model = svmtrain(bag_train_label,bag_train_data(:,IX), cmd);
        model_w(bag,IX) = model.SVs'*model.sv_coef;
        
        
        [predicted_label, accuracy, deci] = svmpredict(test_label,test_data(:,IX),model);
        model_predicted_label(:,bag)=predicted_label;
        
        
        model_deci(:,bag)=deci;

    %     clear  test_data  train_data test_label train_label model IX k
        bag = bag +1;
    end
    fold_w=mean(model_w);
    w(i,IX)  = fold_w(1,IX); %????????
    fold_predicted_label = mode(model_predicted_label,2);
    fold_acc = sum(test_label == fold_predicted_label)*100/length(test_label);
    acc(i,1) = fold_acc;
    fold_deci = mean(model_deci,2);
    deci_value = [deci_value;fold_deci];
    label_nCV = [label_nCV; test_label];
end

w_msk = double(sum(w~=0,1)==size(w,1));
w = mean(w,1).*w_msk;% ??????????????????????????????????????????????????????
acc_final = mean(acc);% ????????
disp(['accuracy - ',num2str(acc_final)]);

% ROC
[X,Y,T,AUC] = perfcurve(label_nCV,deci_value,1);
figure;plot(X,Y);hold on;plot(X,X,'-');
xlabel('False positive rate'); ylabel('True positive rate');

for i=1:length(X)
    Cut_off(i,1) = (1-X(i))*Y(i);
end
[~,maxind] = max(Cut_off);
Specificity = 1-X(maxind);% ?????????????
Sensitivty = Y(maxind);% ??????????????
disp(['Specificity= ', num2str(Specificity)]);
disp(['Sensitivty= ', num2str(Sensitivty)]);
