featureselection.m

function [T_score,T_rank] = featureselection(T_final,ext)

% table variablenames
predictorNames = T_final.Properties.VariableNames(1:end-1);
Y_varname = T_final.Properties.VariableNames{end};
varnames = [{'FSmethod'},predictorNames];

% create a scorecount table
T_score = array2table(zeros(1,length(predictorNames)+1),"VariableNames",varnames);

% create a ranking table
T_rank = array2table(zeros(1,length(predictorNames)+1),"VariableNames",varnames);

% fsrftest
[idx,scores] = fsrftest(T_final,Y_varname);
normalized_scores = normalize(scores,"range");
score_append = [{'Ftest'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'Ftest'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% fsrmrmr
[idx,scores] = fsrmrmr(T_final,Y_varname);
normalized_scores = normalize(scores,"range");
score_append = [{'mRMR'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'mRMR'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% fsrnca* (set rng 0 to 10 and average)
a = 0:1:9;
X = T_final{:,1:end-1};
Y = T_final{:,end};
score_array = [];
for i = 1:length(a)
    rng(a(i))
    temp_mdl = fsrnca(X,Y,'Solver','sgd');
    temp_scores = temp_mdl.FeatureWeights; % Should be a column vector
    score_array = [score_array, reshape(temp_scores,[],1)];
end
scores = mean(score_array,2).';
[~,idx] = sort(scores,"descend");
normalized_scores = normalize(scores,"range");
score_append = [{'NCA'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'NCA'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% relieff
k = 1:10; % number of neighbors
X = T_final{:,1:end-1};
Y = T_final{:,end};
score_array = [];
for i = 1:length(k)
    [temp_idx,temp_scores] = relieff(X,Y,i);
    score_array = [score_array, reshape(temp_scores,[],1)];
end
scores = mean(score_array,2).';
[~,idx] = sort(scores,"descend");
normalized_scores = normalize(scores,"range");
score_append = [{'relieff'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'relieff'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% fitrgp (set rng 0 to 10 and average)
P = length(predictorNames);
sigmaL0 = sqrt(P)*ones(P,1); % Length scale for predictors
sigmaF0 = 1; % Signal standard deviation
sigmaN0 = 1; % Initial noise standard deviation
opts = statset('fitrgp');
opts.TolFun = 1e-2;
a = 0:1:9;
X = T_final{:,1:end-1};
Y = T_final{:,end};
score_array = [];
for i = 1:length(a)
    rng(a(i))
    temp_mdl = fitrgp(X,Y,'KernelFunction','ardsquaredexponential','Verbose',0, ...
        'Optimizer','lbfgs','OptimizerOptions',opts, ...
        'KernelParameters',[sigmaL0;sigmaF0],'Sigma',sigmaN0,'InitialStepSize',1);
    sigmaL = temp_mdl.KernelInformation.KernelParameters(1:end-1); % Learned length scales
    weights = exp(-sigmaL); % Predictor weights
    temp_scores = weights;
    score_array = [score_array, reshape(temp_scores,[],1)];
end
scores = mean(score_array,2).';
[~,idx] = sort(scores,"descend");
normalized_scores = normalize(scores,"range");
score_append = [{'GPR'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'GPR'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% fitrlinear (set rng 0 to 10 and average)
a = 0:1:9;
X = T_final{:,1:end-1};
Y = T_final{:,end};
% split into train-test
hpartition = cvpartition(height(X),"Holdout", 0.2);
idxTrain = training(hpartition);
XTrain = X(idxTrain,:);
YTrain = Y(idxTrain,:);
idxTest = test(hpartition);
XTest = X(idxTest,:);
YTest = Y(idxTest,:);
score_array = [];
Lambda = logspace(-5,10,50);
for i = 1:length(a)
    rng(a(i))
    temp_mdl = fitrlinear(XTrain,YTrain,'Lambda',Lambda,...
        'Learner','leastsquares','Solver','sparsa','Regularization','lasso','PredictorNames',predictorNames);
    mse = loss(temp_mdl,XTest,YTest);
    [~,idxMinMSE] = min(mse);
    temp_scores = temp_mdl.Beta(:,idxMinMSE); % the betas should be the scores
    score_array = [score_array, reshape(temp_scores,[],1)];
end
scores = mean(score_array,2).';
[~,idx] = sort(scores,"descend");
normalized_scores = normalize(scores,"range");
score_append = [{'LR'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'LR'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% lasso (set rng 0 to 10 and average)
a = 0:1:9;
X = T_final{:,1:end-1};
Y = T_final{:,end};
score_array = [];
for i = 1:length(a)
    rng(a(i))
    if height(X) > 10
        [temp_mdl,FitInfo] = lasso(X,Y,'CV',10,'PredictorNames',predictorNames);
    else
        [temp_mdl,FitInfo] = lasso(X,Y,'CV', 2,'PredictorNames',predictorNames);
    end
    idxLambdaMinMSE = FitInfo.IndexMinMSE;
    temp_scores = temp_mdl(:,idxLambdaMinMSE);
    score_array = [score_array, reshape(temp_scores,[],1)];
end
scores = mean(score_array,2).';
[~,idx] = sort(scores,"descend");
normalized_scores = normalize(scores,"range");
score_append = [{'lasso'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'lasso'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% SVM (set rng 0 to 10 and average)
a = 0:1:9;
X = T_final{:,1:end-1};
Y = T_final{:,end};
% split into train-test
hpartition = cvpartition(height(X),"Holdout", 0.2);
idxTrain = training(hpartition);
XTrain = X(idxTrain,:);
YTrain = Y(idxTrain,:);
idxTest = test(hpartition);
XTest = X(idxTest,:);
YTest = Y(idxTest,:);
score_array = [];
for i = 1:length(a)
    rng(a(i))
    temp_mdl = fitrsvm(XTrain,YTrain,'Standardize',true,'KernelFunction','linear');
    mse = loss(temp_mdl,XTest,YTest);
    [~,idxMinMSE] = min(mse);
    temp_scores = temp_mdl.Beta(:,idxMinMSE); % the betas should be the scores
    score_array = [score_array, reshape(temp_scores,[],1)];
end
scores = mean(score_array,2).';
[~,idx] = sort(scores,"descend");
normalized_scores = normalize(scores,"range");
score_append = [{'L-SVM'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'L-SVM'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% GAM (set rng 0 to 10 and average)
a = 0:1:9;
X = T_final{:,1:end-1};
Y = T_final{:,end};
% split into train-test
hpartition = cvpartition(height(X),"Holdout", 0.2);
idxTrain = training(hpartition);
XTrain = X(idxTrain,:);
YTrain = Y(idxTrain,:);
idxTest = test(hpartition);
XTest = X(idxTest,:);
YTest = Y(idxTest,:);
score_array = [];
query_idx = find(idxTest==1);
for i = 1:length(a)
    rng(a(i))
    temp_mdl = fitrgam(XTrain,YTrain,'Interactions','all','MaxPValue',0.05);
    q_score_array = [];
    for q = 1:height(query_idx)
        q_score = zeros(1, length(predictorNames));
        queryPoint = T_final{query_idx(q),1:end-1};
        results = lime(temp_mdl,'QueryPoint',queryPoint,'NumImportantPredictors',length(predictorNames), ...
            'SimpleModelType','tree');
        f = plot(results);
        b = findobj(f,'Type','bar');
        imp_q = b.YData;
        preds_q_text = f.CurrentAxes.YTickLabel;
        preds_qidx = str2double(extractAfter(preds_q_text,'x'));
        q_score(preds_qidx) = normalize(imp_q,"range");
        q_score_array = [q_score_array, reshape(q_score,[],1)];
        close % close the current figure
    end
    temp_scores = mean(q_score_array,2);
    score_array = [score_array, reshape(temp_scores,[],1)];
end
scores = mean(score_array,2).';
[~,idx] = sort(scores,"descend");
normalized_scores = normalize(scores,"range");
score_append = [{'GAM'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'GAM'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% RegressionBaggedEnsemble (set rng 0 to 10 and average)
a = 0:1:9;
X = T_final{:,1:end-1};
Y = T_final{:,end};
score_array = [];
for i = 1:length(a)
    rng(a(i))
    temp_mdl = fitrensemble(X,Y,'Method','Bag','NumLearningCycles',500);
    temp_scores = oobPermutedPredictorImportance(temp_mdl); % Should be a column vector
    score_array = [score_array, reshape(temp_scores,[],1)];
end
scores = mean(score_array,2).';
[~,idx] = sort(scores,"descend");
normalized_scores = normalize(scores,"range");
score_append = [{'RF'}, num2cell(normalized_scores)];
rank_array(idx) = 1:numel(idx);
rank_append = [{'RF'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% get rid of the first rows
T_score(1,:) = [];
T_rank(1,:) = [];

% overall score and rank
% overall score = mean(scores)
% overall rank: add all ranks, then map (sort the numbers and re-number them) the totals into 1:end
% omit methods with all zero scores
mean_row_score = mean(T_score{:,2:end},2); % mean scores of all methods
row_idx = find(mean_row_score > 0);

averge_score = mean(T_score{row_idx,2:end},1); % mean scores of all features
normalized_avgscores = normalize(averge_score,"range");
total_rank = sum(T_rank{row_idx,2:end},1); % total rank of all features
[sorted_rank,sorted_idx] = sort(total_rank);
rank_array(sorted_idx) = 1:numel(sorted_idx);
score_append = [{'Average'}, num2cell(normalized_avgscores)];
rank_append = [{'Average'}, num2cell(rank_array)];
T_score = [T_score; score_append];
T_rank = [T_rank; rank_append];

% savetables
filename = "Score table " + ext + ".xlsx";
writetable(T_score,filename)
filename = "Rank table " + ext + ".xlsx";
writetable(T_rank,filename)


% plot overall scores
figure("Units","normalized","OuterPosition",[0 0 1 1])
[~,idx] = sort(normalized_avgscores,"descend");
avg_rank_array(idx) = 1:numel(idx);
bar(normalized_avgscores(idx))
ylim([0 1])
xlabel("Predictor rank")
ylabel("Predictor importance score (Normalized)")
xticks(1:length(predictorNames))
xticklabels(strrep(predictorNames(idx),"_","\_"))
xtickangle(90)
grid on
titletext = "Arable sensor parameters overall score " + strrep(ext,"_"," ");
title(titletext)
% save figure
savefilename = titletext;
saveas(gcf,[savefilename + ".png"])
savefig(savefilename)

end