function [out, var, mdl] = SFFS_cls(tbl_train, label_train, sett) 

% [out, var] = SFFS_cls(tbl_train, label_train, sett) 
%
% This function performs sequential floating forward selection using the
% input variables: (tbl_train label_train). The columns of tbl_train does
% correspond to the features and the rows corresponds to the observations
%
% Starting from an empty feature set one feature at a time is added based
% on which feature gives the best classification performance, when used
% together with the previously selected features. The implementation does
% the cross-validation process to select the best feature set using input
% settings of the classification algorithm, scroring function, etc.
%
% In this implementation the tradeoff between sensitivity and specificity
% is used to decide whether to include the feature in the selected feature
% set. A tolerance parameter t specifies how many features are allowed to
% be added without improving upon the best classification score obtained
% during the search. The maximum acceptable size of feature set can be set
% by parameter N. Increasing the size of the feature set stops when either
% of the conditions (related to t and N) is reached.
%
%       ------------------------------------------------------------------
%       Required inputs:
% 
%           tbl_train:          Input feature matrix. Columns: features; 
%                               rows: observations. The matrix is numeric
%                               and it represent the training and testing
%                               set of data for SFFS process.
%
%           label_train:        Input labels that corresponds to classes
%                               of the observations from tbl_train. This
%                               algorithm requires labels to be numeric,
%                               where: 1 - positive hypothesis (disease)
%                                      0 - negative hypothesis (healthy)
%
%       ------------------------------------------------------------------
%       Optional inputs:
%
%           num_cross_val.:     The number of cross validation runs that 
%                               the algorithm perform in every search for 
%                               the best features The score is consequen.
%                               averaged in order to obtain better estim.  
%                               of the model accuracy.
%                               Default: 100 cross-validation runs
%
%           k_fold:             The number of cross-validation folds in 
%                               every cross-validation run. 
%                               Default: 10-fold cross-validation
%
%           classifier:         The classifier to use in oreder to do the
%                               classification of the testing data. This
%                               is specified in a structure holding assoc.
%                               member variables settings for the cls
%                               Default: .alg   - 'tree'
%
%           score:              The scoring algorithm that is used to find
%                               the features tobe included or excluded in 
%                               / from the best feature set.
%                               Default: 'tss' - Tradeoff between SEN/SPE
%
%           t:                  The number of features that is allowed to
%                               be included into the selected feature set
%                               without improving obtained best overall 
%                               classification score previously seen with
%                               smaller feature sets.
%                               Default: 3.
%
%           N:                  The maximum number of features, that is to 
%                               be selected as the best possible feat. set
%                               Default: size(tbl_train, 2)
%
%       ------------------------------------------------------------------
%       Outputs:
%
%           out:                Structure holding the output variables:
%                               out.S             - feature indices
%                               out.SCR           - classification score
%                               out.ACC           - classification acc.
%                               out.SEN           - classification sen.  
%                               out.SPE           - classification spe.
%                               out.SCR_std       - classif. score std
%                               out.ACC_std       - classif. acc. std
%                               out.SEN_std       - classif. sen. std
%                               out.SPE_std       - classif. spe. std
%
%           var:                Structure, that holds the change of output
%                               variables during the function runtime. The
%                               score, acc, sen etc. are stored in the var
%                               var.score         - classification score
%                               var.accuracy      - classification acc.
%                               var.sensitivity   - classification sen.
%                               var.specificity   - classification spe.
%                               var.num_features  - num. of features in S
%
%           mdl:                Structure holding trained classification           
%                               models for the selected feature sets. The
%                               internal structure of mdl is specific for
%                               the selected classification technique
%
% Implemented according to:
% A. W. Whitney: "A direct method of nonparametric measurement selection",
% IEEE Transactions on Computers, vol. C-20, no. 9, pp. 1100-1103,
% September 1971.
%
%
%
% --
% ing. Zoltán Galáž
% xgalaz00@stud.feec.vutbr.cz       
% 
% Department of Telecommunications
% Faculty of Electrical Engineering and Communication
% Brno University of Technology

%% Define the optional variable settings
%  [A] Define all the possible classification algorithms and also its
%      associated member variables. The supported classifiers are listed
%      bellow with the description of settings + default values selected
%      if the option of the classifier is not specified, or is specified
%      incorrectly (corrections are displayed into the console)

%      Classification algorithms and settings:
%           a) Support Vector Machines classifier
%              cls.kernel = kernel used in SVM algorithm
%                   {'linear' 'quadratic' 'polynomial' 'rbf' 'mlp'}
%                   (default: 'linear')
%
%           b) Naive Bayes classifier
%               cls.distrb = distribution used in Bayes algorithm
%                   {'normal' 'kernel' 'mvmn' 'mn'} 
%                   (default: 'normal')
%               cls.prior  = prior probability for the classes
%                   {'empirical' 'uniform'} 
%                   (default: 'uniform')
% 
%           c) Linear Discriminant classifier
%               cls.prior = prior probability for the classes
%                   {'empirical' 'uniform'} 
%                   (default: 'uniform')
%               cls.discrm = discrimination type used in the classifier
%                   {'linear' 'quadratic' 'diagLinear' 'diagQuadratic'
%                    'pseudoLinear' 'pseudoQuadratic'} 
%                   (default: 'linear')
% 
%           d) K-Nearest Neighbors classifier
%               cls.k = the number of nearest neighbors 
%                   (default: 1)
%               cls.dist = distance calculated in knn classifier
%                    {'euclidean' 'seuclidean' 'cityblock' 'chebychev'
%                     'minkowski' 'mahalanobis' 'cosine' 'correlation'
%                     'spearman' 'hamming' 'jaccard'} 
%                    (default: 'euclidean')
% 
%           e) Random Forests classifier
%               cls.split = split criterion for ClassificationTree alg
%                   {'gdi' 'twoing' 'deviance'} 
%                   (default: 'gdi')
%               cls.prior = prior probability for the classes
%                   {'empirical', 'uniform'}
%                   (default: 'uniform')
% 
%           f) Gaussian Mixture Models classifier
%               cls.mixt = number of mixtures of GMM classifier
%                   (default: 1)

classifier.def_alg          = 'tree';
classifier.pos_alg          = {'tree' 'svm' 'bayes' 'lda' 'knn' 'gmm'};
classifier.svm.alg          = 'svm';
classifier.svm.pos_kernel   = {'linear' 'quadratic' 'polynomial' ...
                               'rbf' 'mlp'};
classifier.svm.def_kernel   = 'linear';
classifier.bayes.alg        = 'bayes';
classifier.bayes.pos_distrb = {'normal' 'kernel' 'mvmn' 'mn'};
classifier.bayes.def_distrb = 'normal';
classifier.bayes.pos_prior  = {'empirical' 'uniform'};
classifier.bayes.def_prior  = 'uniform';
classifier.lda.alg          = 'lda';
classifier.lda.pos_prior    = {'empirical' 'uniform'};
classifier.lda.def_prior    = 'uniform';
classifier.lda.pos_discrm   = {'linear' 'quadratic' 'diagLinear' ...
                               'diagQuadratic' 'pseudoLinear'    ...
                               'pseudoQuadratic'};
classifier.lda.def_discrm   = 'linear';
classifier.knn.alg          = 'knn';
classifier.knn.pos_k        = 1:size(tbl_train, 2);
classifier.knn.def_k        = 1;
classifier.knn.pos_dist     = {'euclidean' 'seuclidean' 'cityblock'  ...
                               'chebychev' 'minkowski' 'mahalanobis' ...
                               'cosine' 'correlation' 'spearman'     ...
                               'hamming' 'jaccard'};
classifier.knn.def_dist     = 'euclidean';
classifier.tree.alg         = 'tree';
classifier.tree.pos_split   = {'gdi' 'twoing' 'deviance'};
classifier.tree.def_split   = 'gdi';
classifier.tree.pos_prior   = {'empirical' 'uniform'};
classifier.tree.def_prior   = 'uniform';
classifier.gmm.alg          = 'gmm';
classifier.gmm.pos_mixt     = 1:size(tbl_train, 2);
classifier.gmm.def_mixt     = 1;

%  [B] Define the classification scoring algorithm supported. The list of
%      all algorithms is listed bellow. Default scoring algorithm is used
%      if input structure does not contain information about algorithm to
%      use or if it is incorrectly inserted. Default algorithm: 'tss'
%
%      Scoring algorithms (used in binary classification):
%           - acc (classification accuracy)            
%           - f1  (classification f1 score)
%           - inf (classification informedness)
%           - mar (classification markedness)
%           - pre (classification prevalence)
%           - sen (classification sensitivity)
%           - spe (classification sensitivity)
%           - dor (diagnostic odds ratio)
%           - fdr (false discovery rate)
%           - fnr (false negative rate)
%           - for (false omission rate)
%           - fpr (false positive rate)
%           - mcc (Mathews correlation coeff.)
%           - nlr (negative likelihood ratio)
%           - npv (negative predictive value)
%           - plr (positive likelihood ratio)
%           - ppv (positive predictive value)
%           - tss (tradeoff between SEN and SPE)
%           - yjs (Youden's J statistics)

score.pos_alg = {{'acc' '(classification accuracy)'};         ...
                 {'f1'  '(classification f1 score)'};         ...
                 {'inf' '(classification informedness)'};     ...
                 {'mar' '(classification markedness)'};       ...
                 {'pre' '(classification prevalence)'};       ...
                 {'sen' '(classification sensitivity)'};      ...
                 {'spe' '(classification specificity)'};      ...
                 {'dor' '(diagnostic odds ratio)'};           ...
                 {'fdr' '(false discovery rate)'};            ...
                 {'fnr' '(false negative rate)'};             ...
                 {'for' '(false omission rate)'};             ...
                 {'fpr' '(false positive rate)'};             ...
                 {'mcc' '(Mathews correlation coeff.)'};      ...
                 {'nlr' '(negative likelihood ratio)'};       ...
                 {'npv' '(negative predictive value)'};       ...
                 {'plr' '(positive likelihood ratio)'};       ...
                 {'ppv' '(positive predictive value)'};       ...
                 {'yjs' '(Youden''s J statistics)'};          ...
                 {'tss' '(tradeoff between SEN and SPE)'}};
score.def_alg =  'mcc';

%% Set the default/helper variables
%  [A] Helper variables that allow optional variables check and a console 
%      output. Those variables contain the helper variables (set / unset)
%      and default variables that hold the default values for unspecified
%      variables (empty, not specified or out of range). The algorithm is
%      implemented to apply the predefined value

DEF_NUM_CROSS_VAL  = true; % Helper var.: number of C.V. rounds specified
DEF_K_FOLD         = true; % Helper var.: k (number) C.V. folds specified
DEF_CLASSIFIER     = true; % Helper var.: classifier algorithm specified
DEF_SCORE          = true; % Helper var.: cls scoring algorithm specified
DEF_MAX_NUM_FEAT   = true; % Helper var.: max. num. of feat. is specified
DEF_ALLOWED_FEAT   = true; % Helper var.: num. look ahead feat. specified 
DEF_PLOT_VAR       = true; % Helper var.: plotting the var. specified

%  [B] Default settings for the SFFS algorithm. It includes the setting of
%      the number of cross-validation runs and k-folds of cross-validation
%      then the classifier algorithm, scoring function, number of features
%      to look ahead without closing the algorithm and the maximum number
%      of features to be inside the selected feature set

DEFAULT_NUM_CROSS_VAL = 100;                % DEFAULT: num. of C.V. rounds
DEFAULT_K_FOLD        = 10;                 % DEFAULT: nun. of C.V. folds
DEFAULT_CLASSIFIER    = classifier.def_alg; % DEFAULT: classif. algorithm
DEFAULT_SCORE         = score.def_alg;      % DEFAULT: scoring algorithm
DEFAULT_PLOT_VAR      = false;              % DEFAULT: plot var. changes
DEFAULT_DISP_STAT     = true;               % DEFAULT: display properties
DEFAULT_ALLOWED_FEAT  = ...
    floor(size(tbl_train, 2)/2) - 1;        % DEFAULT: look ahead features
DEFAULT_MAX_NUM_FEAT  = ...                 % DEFAULT: max. num. features
            size(tbl_train, 2);

%% Paths and variables
if ((nargin < 3) || (isempty(sett)))  
    
    % Set the number of cross-validations
    sett.num_cross_validations = DEFAULT_NUM_CROSS_VAL;
    DEF_NUM_CROSS_VAL = false;
    
    % Set the k-folds of cross-validation 
    sett.k_fold = DEFAULT_K_FOLD;
    DEF_K_FOLD = false;
    
    % Set the classification algorithm
    sett.classifier = struct();
    sett.classifier.alg = DEFAULT_CLASSIFIER;
    DEF_CLASSIFIER  = false;
    
    % Set the scoring function of classification
    sett.score = DEFAULT_SCORE;
    DEF_SCORE = false;
    
    % Set the maximum number of added (ahead) features
    sett.t = DEFAULT_ALLOWED_FEAT;
    DEF_ALLOWED_FEAT = false;
    
    % Set the maximum number of selected features
    sett.N = DEFAULT_MAX_NUM_FEAT;
    DEF_MAX_NUM_FEAT = false;
    
    % Set the switch for the variable changes plotting 
    sett.plot_var = DEFAULT_PLOT_VAR;
    DEF_PLOT_VAR = false;
    
    % Set the switch for the properties displaying 
    sett.disp_stat = DEFAULT_DISP_STAT;
    
else
    
    % Set the number of cross-validations
    if (~isfield(sett, 'num_cross_validations'))
        sett.num_cross_validations = DEFAULT_NUM_CROSS_VAL;
        DEF_NUM_CROSS_VAL = false; 
    end
    
    % Set the k-folds of cross-validation 
    if (~isfield(sett, 'k_fold'))
        sett.k_fold = DEFAULT_K_FOLD;
        DEF_K_FOLD = false;
    end
    
    % Set the classification algorithm 
    if (~isfield(sett, 'classifier'))
        sett.classifier = struct();
        sett.classifier.alg = DEFAULT_CLASSIFIER;
        DEF_CLASSIFIER  = false;
    end
    
    % Set the scoring function of classification
    if (~isfield(sett, 'score'))
        sett.score = DEFAULT_SCORE;
        DEF_SCORE = false;
    end
    
    % Set the maximum number of added (ahead) features
    if (~isfield(sett, 't'))
        sett.t = DEFAULT_ALLOWED_FEAT;
        DEF_ALLOWED_FEAT = false;
    end
    
    % Set the maximum number of selected features
    if (~isfield(sett, 'N'))
        sett.N = DEFAULT_MAX_NUM_FEAT;
        DEF_MAX_NUM_FEAT = false;
    end
    
    % Set the switch for the variable changes plotting
    if (~isfield(sett, 'plot_var'))
        sett.plot_var = DEFAULT_PLOT_VAR;
        DEF_PLOT_VAR = false;
    end
    
    % Set the switch for the properties displaying 
    if (~isfield(sett, 'disp_stat'))
        sett.disp_stat = DEFAULT_DISP_STAT;
    end
end

sett.score = lower(sett.score);
sett.classifier.alg = lower(sett.classifier.alg);

%% Test the input data   
%  [A] Check if the number of observations in 'tbl_train' is equal to the
%      number of observations in 'label_train'. The algorithm demands the 
%      same number of observations (rows) = e.g. patients (PD/HC, etc.)

if (size(tbl_train, 1) ~= size(label_train, 1))
    error('The number of observations in tbl_train ~= label_train');
end

%  [B] Optional variables handling routines to prepare the console output
%
%      Handle: The number of cross-validation runs
%              DEFAULT: 100 cross-validation runs
if (~DEF_NUM_CROSS_VAL)
    STR_NUM_CROSS_VAL = num2str(sett.num_cross_validations);
    STR_NUM_CROSS_VAL_STATE = ' (DEFAULT)';
elseif ((DEF_NUM_CROSS_VAL) && (sett.num_cross_validations <= 0))
    sett.num_cross_validations = DEFAULT_NUM_CROSS_VAL;
    STR_NUM_CROSS_VAL = num2str(sett.num_cross_validations);
    STR_NUM_CROSS_VAL_STATE = ' (AUTOMATIC)';
else
    STR_NUM_CROSS_VAL = num2str(sett.num_cross_validations);
    STR_NUM_CROSS_VAL_STATE = '';
end

%      Handle: The number of cross-validation folds (k)
%              DEFAULT: 10 cross-validation folds
if (~DEF_K_FOLD)
    STR_K_FOLD = [num2str(sett.k_fold) '-fold'];
    STR_K_FOLD_STATE = ' (DEFAULT)';
elseif ((DEF_K_FOLD) && ((sett.k_fold <= 0) || ...
        (sett.k_fold > floor(size(tbl_train, 1)/2))))
    sett.k_fold = DEFAULT_K_FOLD;
    STR_K_FOLD = [num2str(sett.k_fold) '-fold'];
    STR_K_FOLD_STATE = ' (AUTOMATIC)';
elseif ((DEF_K_FOLD) && (sett.k_fold == 1))    
    STR_K_FOLD = 'L-O-OUT';
    STR_K_FOLD_STATE = '';

    if (sett.num_cross_validations ~= 1)
        sett.num_cross_validations = 1;
        STR_NUM_CROSS_VAL = num2str(sett.num_cross_validations);
        STR_NUM_CROSS_VAL_STATE = ' (AUTOMATIC)';
    end
else       
    STR_K_FOLD = [num2str(sett.k_fold) '-fold'];
    STR_K_FOLD_STATE = '';
end

%      Handle: Classification algorithm used in the feature selection
%              DEFAULT: ClassificatioTree: (Random Forest) classifier
if (~DEF_CLASSIFIER)
    STR_CLASSIFIER = sett.classifier.alg;
    STR_CLASSIFIER_STATE = ' (DEFAULT)';
    sett.classifier = ...
            check_classification_sett(sett.classifier, classifier);
else

    % check the classifier (set the classification algorithm)
    sett.classifier.alg = check_classification_alg(sett.classifier.alg);
    if (~isempty(sett.classifier.alg))
         STR_CLASSIFIER = sett.classifier.alg;
         STR_CLASSIFIER_STATE = '';
         sett.classifier = ...
            check_classification_sett(sett.classifier, classifier);
    else
        disp(' '                                        );
        disp('Unknown classifier, set to DEFAULT: tree' );
        disp('Known classifiers:'                       );
        disp('   - tree  (Classification Tree)'         );
        disp('   - svm   (Support Vector Machines)'     );
        disp('   - knn   (K-Nearest Neighbors)'         );
        disp('   - gmm   (Gaussian Mixture Models)'     );
        disp('   - lda   (Linear discriminant Analysis)');
        disp('   - bayes (Naive Bayes Classifier)'      );
        
        STR_CLASSIFIER = DEFAULT_CLASSIFIER;
        STR_CLASSIFIER_STATE = ' (AUTOMATIC)';
        sett.classifier = ...
            check_classification_sett(sett.classifier, classifier);
    end
end

%      Handle: The ccoring method used used in the feature selection
%              DEFAULT: 'tss' (tradeoff between SEN and SPE)
if (~DEF_SCORE)
    STR_SCORE = sett.score;
    STR_SCORE_STATE = ' (DEFAULT)';
else

    % check the classification score (set the classification score)
    sett.score = check_classification_score(sett.score);
    if (~isempty(sett.score))
         STR_SCORE = sett.score;
         STR_SCORE_STATE = '';
    else
        disp(' '                                           );
        disp('Unknown scoring method, set to DEFAULT: tss' );
        disp('Known scoring methods:'                      );
        disp('   - acc (classification accuracy)'          );            
        disp('   - f1  (classification f1 score)'          );
        disp('   - inf (classification informedness)'      );
        disp('   - mar (classification markedness)'        );
        disp('   - pre (classification prevalence)'        );
        disp('   - sen (classification sensitivity)'       );
        disp('   - spe (classification sensitivity)'       );
        disp('   - dor (diagnostic odds ratio)'            );
        disp('   - fdr (false discovery rate)'             );
        disp('   - fnr (false negative rate)'              );
        disp('   - for (false omission rate)'              );
        disp('   - fpr (false positive rate)'              );
        disp('   - mcc (Mathews correlation coeff.)'       );
        disp('   - nlr (negative likelihood ratio)'        );
        disp('   - npv (negative predictive value)'        );
        disp('   - plr (positive likelihood ratio)'        );
        disp('   - ppv (positive predictive value)'        );
        disp('   - tss (tradeoff between SEN and SPE)'     );

        STR_SCORE = DEFAULT_SCORE;
        STR_SCORE_STATE = ' (AUTOMATIC)';
    end
end

%      Handle: Maximum number of features to look ahead
%              BOUNDS:  <0, floor(size(tbl_train, 2)/2) - 1) € Z 
%              DEFAULT: 3 allowed features (local max. possibility)
if (~DEF_ALLOWED_FEAT)
    STR_ALLOWED_FEAT = num2str(sett.t);
    STR_ALLOWED_FEAT_STATE = ' (DEFAULT)';
else
    if (sett.t >= floor(size(tbl_train, 2)/2) || (sett.t < 0))
        pref_value = floor(size(tbl_train, 2)/2) - 1;
        pref_str   = num2str(pref_value);
        disp(' '                                                     );
        disp(['Improper number of added features, set to: ' pref_str]);
        disp('Allowed values <0, floor(size(tbl_train, 2))/2 - 1)'   );
        disp('WARNING: values € Z'                                   );

        sett.t = DEFAULT_ALLOWED_FEAT;
        STR_ALLOWED_FEAT = num2str(sett.t);
        STR_ALLOWED_FEAT_STATE = ' (AUTOMATIC)';
    else
        STR_ALLOWED_FEAT = num2str(sett.t);
        STR_ALLOWED_FEAT_STATE = '';
    end
end

%      Handle: Maximum number of features selected into output feat. set
%              BOUNDS:  <1, size(tbl_train, 2)) € Z 
%              DEFAULT: max. num. of features == the number of features
if (~DEF_MAX_NUM_FEAT)
    STR_MAX_NUM_FEAT = num2str(sett.N);
    STR_MAX_NUM_FEAT_STATE = ' (DEFAULT)';
else
    if (sett.N >= size(tbl_train, 2) || (sett.N <= 0))
        str = 'size(tbl_train, 2)';
        disp(' '                                                    );
        disp(['Improper max. number features, set to DEFAULT: ' str]);
        disp('Allowed values <1, size(tbl_train, 2))'               );
        disp('WARNING: values € Z'                                  );

        sett.N = DEFAULT_MAX_NUM_FEAT;
        STR_MAX_NUM_FEAT = num2str(sett.N);
        STR_MAX_NUM_FEAT_STATE = ' (AUTOMATIC)';
    else
        STR_MAX_NUM_FEAT = num2str(sett.N);
        STR_MAX_NUM_FEAT_STATE = '';
    end
end

%      Handle: The switch to plot the variable changes in time
%              DEFAULT: turned off (false)
if (~DEF_PLOT_VAR)
    STR_PLOT_VAR = num2str(sett.plot_var);
    STR_PLOT_VAR_STATE = ' (DEFAULT)';
else
    STR_PLOT_VAR = num2str(sett.plot_var);
    STR_PLOT_VAR_STATE = '';
end

%  [C] Console output of the input variables and warnings (if occurred)
if (sett.disp_stat)
    fprintf('\n'                                                    );
    fprintf('INPUT VARIABLES\n'                                     );
    fprintf('training matrix:\t\t(OK): specified\n'                 );
    fprintf('training labels:\t\t(OK): specified\n'                 );
    fprintf('\n'                                                    );
    fprintf('OPTIONAL VARIABLES\n'                                  );
    fprintf('Cross validations:\t\t(OK): %s%s\n',  ...
        STR_NUM_CROSS_VAL, STR_NUM_CROSS_VAL_STATE                  );
    fprintf('K-fold       \t\t\t(OK): %s%s\n',     ...
        STR_K_FOLD, STR_K_FOLD_STATE                                );
    fprintf('Classifier:\t\t\t\t(OK): %s%s\n',     ...
        STR_CLASSIFIER, STR_CLASSIFIER_STATE                        );
    fprintf('Scoring method:\t\t\t(OK): %s%s\n',   ...
        STR_SCORE, STR_SCORE_STATE                                  );
    fprintf('Added features:\t\t\t(OK): %s%s\n',   ...
        STR_ALLOWED_FEAT, STR_ALLOWED_FEAT_STATE                    );
    fprintf('Max. num. features:\t\t(OK): %s%s\n', ...
        STR_MAX_NUM_FEAT, STR_MAX_NUM_FEAT_STATE                    );
    fprintf('Plot var. progress:\t\t(OK): %s%s\n', ...
        STR_PLOT_VAR, STR_PLOT_VAR_STATE                            );
    fprintf('\n'                                                    );
end

%% Set the variables
out         = struct(); % Empty struct of output variables (features, etc)
var         = struct(); % Empty struct of output variable changes in time
mdl         = struct(); % Empty struct of trained classification models
out.S       = [];       % Empty vector of best performing features (idx)
out.SCR     = [];       % Empty vector of classification scores (optional)
out.ACC     = [];       % Empty vector of classification accuracy (ACC)
out.SEN     = [];       % Empty vector of sensitivity of clasification
out.SPE     = [];       % Empty vector of specificity of clasification
out.SCR_std = [];       % Empty vector of classif. scores std (optional)
out.ACC_std = [];       % Empty vector of classif. accuracy std (ACC)
out.SEN_std = [];       % Empty vector of classif. sensitivity std (SEN)
out.SPE_std = [];       % Empty vector of classif. specificity std (SPE)

total_best_val = -inf;  % Total best value of SCORE in the SFFS func. call       
iterator = 0;           % Helper variable to iterate over variable changes
loop = 1 + sett.t;      % Control variable (counter of loops in while)

mdl.cls_method = STR_CLASSIFIER;
mdl.cls_score  = STR_SCORE;
mdl.CV_folds   = STR_K_FOLD;
mdl.num_CV     = STR_NUM_CROSS_VAL;
mdl.cls_models = cell(1, sett.N);

%% Set the possible classes (check classes)
ncv = sett.num_cross_validations;
cls = sett.classifier;

if (sett.k_fold == 1)
    fld = size(tbl_train, 1);
else
    fld = sett.k_fold;
end

best_val    = [];
accuracy    = [];
sensitivity = [];
specificity = [];

score_std       = [];
accuracy_std    = [];
sensitivity_std = [];
specificity_std = [];

pred_labels = cell(ncv, fld);
true_labels = cell(ncv, fld);

%% Perform Sequential Floating Forward Selection
while (loop > 0)

    % [A] EVALUATE THE UNION OF THE CURRENTLY SELECTED FEATURES WITH EACH 
    %     YET UNSELECTED FEATURE 
    bestidx                            = 0;
    best_val(length(out.S) + 1)        = -inf;
    accuracy(length(out.S) + 1)        = 0;  
    sensitivity(length(out.S) + 1)     = 0;
    specificity(length(out.S) + 1)     = 0;
    score_std(length(out.S) + 1)       = 0;
    accuracy_std(length(out.S) + 1)    = 0;  
    sensitivity_std(length(out.S) + 1) = 0;
    specificity_std(length(out.S) + 1) = 0;
    
    for i1 = 1:size(tbl_train, 2)
        if ~ismember(i1, out.S)
            
            % a) Select the best features so far plus a random feature and
            %    set the index vector as the union of S (feat. selected so
            %    far) and the actual iteration number (idx), if feature is
            %    not a member of S yet
            sel = [out.S i1];

            % b) Perform the cross validations (number of validations is
            %    specified by the settings of the algorithm). In each of
            %    the validation steps data is randomly shuffled and next
            %    set (training/testing) data's selected (stratification)
            for cvi = 1:ncv
                
                % Prepare cross validation data and labels for the
                % cross validation. From this set the k-folds will 
                % be executed. It also randomizes the data
                if (sett.k_fold == 1)
                    cv = cvpartition(numel(label_train), 'LeaveOut');
                else
                    cv = cvpartition(numel(label_train), ...
                        'KFold', sett.k_fold);
                end
                
                % Create temporary vector (predictions/presennted)
                predictions = cell(1, fld);
                presented   = cell(1, fld);
                
                % Iterate over all folds of the cross validation
                % and do the training and testing process of the
                % chosen rclassification algorithm
                for fold = 1:fld

                    % Set the indices for the cross-validation
                    train_idx = cv.training(fold);
                    test_idx  = cv.test(fold);

                    % Perform the estimation   
                    train_tbl = tbl_train(train_idx, sel);
                    train_lbl = label_train(train_idx);
                    test_tbl  = tbl_train(test_idx, sel);
                    test_lbl  = label_train(test_idx);

                    [~, pred_nested, cls_model] = ...
                        perf_classification(train_tbl, train_lbl, ...
                            test_tbl, cls);

                    predictions{fold} = pred_nested(:);
                    presented{fold}   = test_lbl(:);
                end
                
                % Update the true and predicted labels
                pred_labels(cvi, :) = predictions;
                true_labels(cvi, :) = presented;
            end
            
            % c) Compute overall classification performance score that 
            %    is achieved so far (default: 'tss' class. scoring alg.) 
            ACC_cv = zeros(ncv, fld);
            SEN_cv = zeros(ncv, fld);
            SPE_cv = zeros(ncv, fld);            
            SCR_cv = zeros(ncv, fld);

            for cvi = 1:ncv
                for fold = 1:fld
                    
                    true_D = (true_labels{cvi, fold} == 1);
                    detc_D = (pred_labels{cvi, fold} == 1);
                    true_H = (true_labels{cvi, fold} == 0);
                    detc_H = (pred_labels{cvi, fold} == 0);

                    % Calculate the values of the confusion matrix
                    TP = sum(true_D & detc_D);
                    TN = sum(true_H & detc_H);
                    FP = sum(true_H & detc_D);
                    FN = sum(true_D & detc_H);

                    % Calculate the classification scores
                    ACC_cv(cvi, fold) = calc_classifiation_score( ...
                        TP, FP, FN, TN, 'ACC', 1);
                    SEN_cv(cvi, fold) = calc_classifiation_score( ...
                        TP, FP, FN, TN, 'SEN', 1);
                    SPE_cv(cvi, fold) = calc_classifiation_score( ...
                        TP, FP, FN, TN, 'SPE', 1);
                    SCR_cv(cvi, fold) = calc_classifiation_score( ...
                        TP, FP, FN, TN, sett.score);
                end 
            end
            
            % Calculate the mean scores over all cross validations
            ACC_temp = mean(ACC_cv(:));
            SEN_temp = mean(SEN_cv(:));
            SPE_temp = mean(SPE_cv(:));
            SCR_temp = mean(SCR_cv(:));
            
            ACC_std  = std(ACC_cv(:));
            SEN_std  = std(SEN_cv(:));
            SPE_std  = std(SPE_cv(:));
            SCR_std  = std(SCR_cv(:));
            
            % d) Check whether obtained score is better than the old one 
            %    if yes set the currently included feature as it is best
            %    best candidate and update the best score
            if (SCR_temp > best_val(length(out.S) + 1))
                bestidx = i1;
                best_val(length(out.S) + 1)    = SCR_temp;
                accuracy(length(out.S) + 1)    = ACC_temp;
                sensitivity(length(out.S) + 1) = SEN_temp;
                specificity(length(out.S) + 1) = SPE_temp;
                
                score_std(length(out.S) + 1)       = SCR_std;
                accuracy_std(length(out.S) + 1)    = ACC_std;
                sensitivity_std(length(out.S) + 1) = SEN_std;
                specificity_std(length(out.S) + 1) = SPE_std;
                
                iterator                   = +iterator + 1;
                var.num_features(iterator) = length(out.S);
                var.score(iterator)        = SCR_temp;
                var.accuracy(iterator)     = ACC_temp;
                var.sensitivity(iterator)  = SEN_temp;
                var.specificity(iterator)  = SPE_temp;
                
                var.score_std(iterator)       = SCR_std;
                var.accuracy_std(iterator)    = ACC_std;
                var.sensitivity_std(iterator) = SEN_std;
                var.specificity_std(iterator) = SPE_std;                
                                
                statstr = ['SFFS(' num2str(length(out.S))   ...
                    '): ' num2str(i1) ':'                   ...
                    ' SCR = ' num2str(SCR_temp)             ...
                    ' +- ' num2str(SCR_std)                 ...
                    ' ACC = ' num2str(ACC_temp)             ...
                    ' +- ' num2str(ACC_std)                 ...
                    ' SEN = ' num2str(SEN_temp)             ...
                    ' +- ' num2str(SEN_std)                 ...
                    ' SPE = ' num2str(SPE_temp)             ...
                    ' +- ' num2str(SPE_std)];
                if (SCR_temp >= total_best_val)
                    statstr = [statstr ' *'];
                end     
                
                disp(statstr);
            end
        end
    end

    % [B] INCLUDE THE BEST FEATURE IN THE SET OF SELECTED FEATURES
    out.S   = [out.S bestidx];
    out.SCR = [out.SCR best_val(length(out.S))];
    out.ACC = [out.ACC accuracy(length(out.S))];
    out.SEN = [out.SEN sensitivity(length(out.S))];
    out.SPE = [out.SPE specificity(length(out.S))];
    
    out.SCR_std = [out.SCR_std score_std(length(out.S))];
    out.ACC_std = [out.ACC_std accuracy_std(length(out.S))];
    out.SEN_std = [out.SEN_std sensitivity_std(length(out.S))];
    out.SPE_std = [out.SPE_std specificity_std(length(out.S))];
    
    mdl.cls_models{1, length(out.S)} = cls_model;
    
    % [C] "FLOATING" EXCLUSION
    %     Exclude thus far selected features one at a time and find the 
    %     least useful one. Repeat this step as long as it's evaluation 
    %     score stay better than with the previous best score, by using
    %     feature set of the same size 
    first_exclusion = true;
    while (length(out.S) > 2)
        
        % (1) Find the least useful feature to be excluded from the set
        %     of the best festures selected so far by incrementing over
        %     the set and temporary excluding one feature at the time.
        exclusion_score = -inf; 
        excidx = 0;
        for i1 = 1:length(out.S)
            
            % a) Select one feature from the selected feature set S at 
            %    the time, to exclude it from the selected features to  
            %    find the least useful one ("floating" exclusion)
            S_x = setdiff(out.S, out.S(i1));
            
            % b) Perform the cross validations (number of validations is
            %    specified by the settings of the algorithm). In each of
            %    the validation steps data is randomly shuffled and next
            %    set (training/testing) data's selected (stratification)
            for cvi = 1:ncv
                    
                % Prepare cross validation data and labels for the
                % cross validation. From this set the k-folds will 
                % be executed. It also randomizes the data        
                if (sett.k_fold == 1)
                    cv = cvpartition(numel(label_train), 'LeaveOut');
                else
                    cv = cvpartition(numel(label_train), ...
                        'KFold', sett.k_fold);
                end
                
                % Create temporary vector (predictions/presennted)
                predictions = cell(1, fld);
                presented   = cell(1, fld);

                % Iterate over all folds of the cross validation
                % and do the training and testing process of the
                % chosen classification algorithm
                for fold = 1:sett.k_fold
                    
                    % Set the indices for the cross-validation
                    train_idx = cv.training(fold);
                    test_idx  = cv.test(fold);
                    
                    % Perform the classification
                    train_tbl = tbl_train(train_idx, sel);
                    train_lbl = label_train(train_idx);
                    test_tbl  = tbl_train(test_idx, sel);
                    test_lbl  = label_train(test_idx);

                    [~, pred_nested] = ...
                        perf_classification(train_tbl, train_lbl, ...
                            test_tbl, cls);
                    
                    predictions{fold} = pred_nested(:);
                    presented{fold}   = test_lbl(:);
                end
                
                % Update the true and predicted labels
                pred_labels(cvi, :) = predictions;
                true_labels(cvi, :) = presented;
            end
            
            % c) Compute overall classification performance score that 
            %    is achieved so far (default: 'tss' class. scoring alg.)                     
            SCR_cv = zeros(ncv, fld);
            
            for cvi = 1:ncv
                for fold = 1:fld
                    
                    true_D = (true_labels{cvi, fold} == 1);
                    detc_D = (pred_labels{cvi, fold} == 1);
                    true_H = (true_labels{cvi, fold} == 0);
                    detc_H = (pred_labels{cvi, fold} == 0);
                    
                    % Calculate the values of the confusion matrix
                    TP = sum(true_D & detc_D);
                    TN = sum(true_H & detc_H);
                    FP = sum(true_H & detc_D);
                    FN = sum(true_D & detc_H);
                    
                    % Calculate the classification scores
                    SCR_cv(cvi, fold) = calc_classifiation_score( ...
                        TP, FP, FN, TN, sett.score);
                end 
            end
            
            % Calculate the mean scores over all cross validations 
            SCR_temp = mean(SCR_cv(:)); 
            
            % d) look for the best feature to exclude from the selected
            %    feature set S. The condition is to find the feature to
            %    produce the lowest classification score (ACC). This
            %    means that after excluding the least useful feature,
            %    the classification produces higher accuracy, than if 
            %    other features are excluded
            if (SCR_temp > exclusion_score)
                exclusion_score = SCR_temp;
                excidx = out.S(i1);
            end
        end
        
        % (2) If the least useful feature is the most recently added
        %     feature, or no gain can be obtained by exclusion don't
        %     exclude anything and continue with forward selection
        if (first_exclusion && excidx == bestidx) || ...
                (exclusion_score <= best_val(length(out.S) - 1))
            break;
            
        % (3) If exclusion leads to better score, than what has been 
        %     previously obtained with feature set of the same size 
        %     exclude the least useful feature
        else
            idx          = find(out.S == excidx);
            out.S(idx)   = [];
            out.SCR(idx) = [];
            out.ACC(idx) = [];
            out.SEN(idx) = [];
            out.SPE(idx) = [];
            
            out.SCR_std(idx) = [];
            out.ACC_std(idx) = [];
            out.SEN_std(idx) = [];
            out.SPE_std(idx) = [];
            
            mdl.cls_models(idx) = [];
            
            best_val(length(out.S)) = exclusion_score;
            
            % reset counter after feature exclusion
            loop = 1 + sett.t; 
            
            iterator                   = +iterator + 1;
            var.num_features(iterator) = length(out.S);
            var.score(iterator)        = best_val(length(out.S));
            var.accuracy(iterator)     = accuracy(length(out.S));
            var.sensitivity(iterator)  = sensitivity(length(out.S));
            var.specificity(iterator)  = specificity(length(out.S));
            
            var.score_std(iterator)       = score_std(length(out.S));
            var.accuracy_std(iterator)    = accuracy_std(length(out.S));
            var.sensitivity_std(iterator) = sensitivity_std(length(out.S));
            var.specificity_std(iterator) = specificity_std(length(out.S));
            
            statstr = ['Exclusion (' num2str(length(out.S)) '-1): ' ...
                num2str(excidx) ': ' num2str(exclusion_score)];
            if (exclusion_score >= total_best_val)
                statstr = [statstr ' *'];
            end
            disp(statstr);
        end
        
        first_exclusion = false;
    end
    
    % [D] STOPPING CONDITIONS AND IMPROVEMENT CHECK
    %     At this phase, this algorithm controls if the length of the
    %     feature set (S) is not larger or equal to specified maximum 
    %     number of features. If the size wasn't specified it was set
    %     to be equal to the number of columns (features) in training
    %     and testing tables. The algorithm also checks, if there was
    %     any improvement obtained in the last iteration.
    
    % (1) Quit if the specified maximum number of features is reached
    if (length(out.S) >= sett.N)
        break;
    end
    
    % (2) Check, if any improvement was made after the last iteration
    if (best_val(length(out.S)) <= total_best_val)
        
        % a) The feature set (S) after inclusion and exclusion didn't 
        %    improve the overall best score => decrease the counter 
        loop = loop - 1;
    else
        
        % b) An improvement was obtained by the new feature set => do
        %    the reset of the counter to its starting value
        loop = 1 + sett.t; 
        total_best_val = best_val(length(out.S));
    end    
end

%% Back off to the best performing feature set
[~, max_idx] = max(out.SCR);

out.S   = out.S(1:max_idx);
out.SCR = out.SCR(1:max_idx);
out.ACC = out.ACC(1:max_idx);
out.SEN = out.SEN(1:max_idx);
out.SPE = out.SPE(1:max_idx);

out.SCR_std = out.SCR_std(1:max_idx);
out.ACC_std = out.ACC_std(1:max_idx);
out.SEN_std = out.SEN_std(1:max_idx);
out.SPE_std = out.SPE_std(1:max_idx);

mdl.cls_models = mdl.cls_models(1:max_idx);

if (sett.plot_var)
    plot_sffs_cls_progress(var);
end