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

% [out, var] = SFFS_reg(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 estimation (reg) 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 regression algorithm, scroring function, etc.
%
% In this implementation the mean absolute error is used to decide whether 
% to include the feature in the selected feature set. Tolerance parameter
% t specifies how many features are allowed to be added without improving
% upon the best regression 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: classes are continuous numbers
%
%       ------------------------------------------------------------------
%       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
%                               estimation (reg) of the testing data. This 
%                               is specified in a structure holding assoc.
%                               member variables settings for the cls
%                               Default: .alg  - 'cart'
%
%           score:              The scoring algorithm that is used to find
%                               the features tobe included or excluded in 
%                               / from the best feature set.
%                               Default: 'mae' - mean absolute error
%
%           t:                  The number of features that is allowed to
%                               be included into the selected feature set
%                               without improving obtained best overall 
%                               estimation (reg) 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           - Regression score
%                               out.MAE           - Mean absolute error
%                               out.RMSE          - Root mean sqr. error
%                               out.EER           - Estimation err. (rel.)
%                               out.RSS           - Residual sum of squares
%                               out.COD           - Coeff. of determination
%                               out.RHO           - Spearman's corr. c. 
%                               out.P             - Rho significance lvl. 
%                               out.SCR_std       - Regression score std
%                               out.MAE_std       - Mean abs. err. std
%                               out.RMSE_std      - Root mean s. err. std
%                               out.EER_std       - Estim. err. (rel.) std
%                               out.RSS_std       - Resid. sum of sq. std
%                               out.COD_std       - Coeff. of determ. std
%                               out.RHO_std       - Spearman's corr. std
%                               out.P_std         - Rho signif. lvl. 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         - Regression score
%                               var.mae           - Mean squared error
%                               var.rmse          - Root mean squared err.
%                               var.eer           - Estimation error
%                               var.rss           - Residual sum of squares
%                               var.cod           - Coeff. of determination
%                               var.rho           - Spearman's corr. coeff
%                               var.p             - Significance level 
%                               var.num_features  - Num. of features in S
%
%           mdl:                Structure, that holds trained regression
%                               models for the selected feature sets. The
%                               internal structure of mdl is specific for
%                               the selected regression 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 estimation (reg) algorithms and also its
%      associated member variables. The supported regression algorithms 
%      are listed bellow with the description of settings + its default 
%      values selected if the option of the classifier is not specified, 
%      or is specified incorrectly (corrections are displayed into the 
%      console)

%      Regression algorithms and settings:
%           a) Classification and Regression Trees algorithm
%               cls.splitcriterion = Split criterion for Cart
%                   {'mse'} 
%                   (default: 'mse')
%               cls.minleaf = Minimum number of leaf node observations
%                   {1 ... positive integer value}
%                   (default: 1
%               cls.minparent = Minimum number of branch node observations
%                   {10 ... positive integer value}
%                   (default: 10)
%               cls.prune = Pruning flag
%                   {'on' 'off'}
%                   (default: 'on')
%               cls.prunecriterion = Pruning criterion
%                   {'mse'}
%                   (default: 'mse)
%
%           b) Ordinal Regression algorithm
%               cls.model = Type of model to fit 
%                   {'nominal' 'ordinal' 'hierarchical'}
%                   (default: 'ordinal')
%               cls.link  = Link function 
%                   {'logit' 'probit' 'comploglog' 'loglog'}
%                   (default: 'logit')
%
%           c) Linear Regression algorithm
%               cls.modelspec = Model specification
%                   {'constant', 'linear', 'interactions', ...
%                    'purequadratic', 'quadratic', 'polyijk'} 
%                   (default: 'linear')
%               cls.intercept = Indicator for constant term
%                   {'true', 'false'}
%                   (default: 'true')
%               cls.robustopts = Indicator of robust fitting type
%                   {'on', 'off', string}
%                    string: {'andrews', 'bisquare', 'cauchy', ...
%                             'fair', 'huber', 'logistic',     ...
%                             'ols', 'talwar', 'welsch'}
%                   (default (robustopts): 'off')
%                   (default (string): 'bisquare')
%
%           d) SVM Regression algorithm
%               cls.kernel = Kernel function
%                   {'linear', 'gaussian', 'rbf', 'polynomial'} 
%                   (default: 'linear')
%
%           e) Gaussian Process Regression algorithm
%               cls.basis = Explicit basis in the GPR model
%                    {'constant', 'none', 'linear', 'pureQuadratic'} 
%                    (default: 'constant')
%
%           e) Ensemble Regression
%               cls.method = Ensemble-aggregation method
%                   {'LSBoost', 'Bag'} 
%                   (default: 'LSBoost')
%
%           f) Random Forests Regression
%               cls.num_trees = number of trees 
%                   {1 ... positive integer value}
%                   (default: 50)
%               cls.prior = prior probability for the classes
%                   {'Empirical', 'Uniform', []} 
%                   (default: [])
%
%           g) Gradient Boosted Trees Regression
%               cls.loss = Shrinkage factor (max 1.0)
%                   {'squaredloss', 'logloss', 'exploss'} 
%                   (default: 'squaredloss')
%               cls.shrinkageFactor = Shrinkage factor (max 1.0)
%                   {0. ... positive float value} 
%                   (default: 0.1)
%               cls.subsamplingFactor = Subsampling factor
%                   {0. ... positive float value} 
%                   (default: 0.2)
%               cls.maxTreeDepth = Maximum depth of the trees
%                   {uint32(int) ... uint32 of integer value}
%                   (default: uint32(2))
%               cls.randSeed = Random seed
%                   {uint32(int) ... uint32 of integer value}
%                   (default: uint32(1))
%               cls.numIters = Number of fitting iterations
%                   {uint32(int) ... uint32 of integer value}
%                   (default: uint32(500))

classifier.def_alg                    = 'cart';
classifier.pos_alg                    = {'cart', 'or', 'linr', 'svmr', ...
                                         'gpr', 'emr', 'rfr', 'gbrt'};
classifier.def_alg                    = 'cart';
classifier.cart.alg                   = 'cart';
classifier.cart.pos_splitcriterion    = 'mse';
classifier.cart.def_splitcriterion    = 'mse';
classifier.cart.pos_minleaf           = 1:1:10;
classifier.cart.def_minleaf           = 1;
classifier.cart.pos_minparent         = 10:1:100;
classifier.cart.def_minparent         = 10;
classifier.cart.pos_prune             = {'on' 'off'};
classifier.cart.def_prune             = 'on';
classifier.cart.pos_prunecriterion    = 'mse';
classifier.cart.def_prunecriterion    = 'mse';
classifier.or.alg                     = 'or';
classifier.or.pos_model               = {'nominal' 'ordinal'           ...
                                         'hierarchical'};
classifier.or.def_model               =  'ordinal';
classifier.or.pos_link                = {'logit' 'probit' 'comploglog' ...
                                         'loglog'};
classifier.or.def_link                = 'logit';
classifier.or.pos_model               = {'nominal' 'ordinal'           ...
                                         'hierarchical'};
classifier.or.def_model               = 'ordinal';
classifier.linr.alg                   = 'linr';
classifier.linr.pos_modelspec         = {'linear', 'constant',         ...
                                         'polyijk' 'interactions',     ...
                                         'purequadratic', 'quadratic'};
classifier.linr.def_modelspec         = 'linear';
classifier.linr.pos_intercept         = [0, 1];
classifier.linr.def_intercept         = 1;
classifier.linr.pos_robustopts        = {'on' 'off'};
classifier.linr.def_robustopts        = 'on';
classifier.svmr.alg                   = 'svmr';
classifier.svmr.pos_kernel            = {'linear', 'gaussian', 'rbf',  ...
                                         'polynomial'};
classifier.svmr.def_kernel            = 'linear';
classifier.gpr.alg                    = 'gpr';
classifier.gpr.pos_basis              = {'constant', 'none', 'linear', ...
                                         'pureQuadratic'};
classifier.gpr.def_basis              = 'constant';
classifier.emr.alg                    = 'emr';
classifier.emr.pos_method             = {'LSBoost', 'Bag'};
classifier.emr.def_method             = 'LSBoost';
classifier.rfr.alg                    = 'rfr';
classifier.rfr.pos_num_trees          = 1:1:50;
classifier.rfr.def_num_trees          = 10;
classifier.rfr.pos_prior              = {'Empirical', 'Uniform', []} ;
classifier.rfr.def_prior              = [];
classifier.gbtr.alg                   = 'gbtr';
classifier.gbtr.def_loss              = 'squaredloss';
classifier.gbtr.pos_loss              = {'squaredloss', 'logloss',     ...
                                         'exploss'};
classifier.gbtr.def_shrinkageFactor   = 0.1;
classifier.gbtr.pos_shrinkageFactor   = 0.01:0.01:0.1;
classifier.gbtr.def_subsamplingFactor = 0.2;
classifier.gbtr.pos_subsamplingFactor = 0.2:0.1:1.0;
classifier.gbtr.def_maxTreeDepth      = uint32(2);
classifier.gbtr.pos_maxTreeDepth      = uint32(2):uint32(1):uint32(10);
classifier.gbtr.def_randSeed          = uint32(1);
classifier.gbtr.pos_randSeed          = uint32(1):uint32(1):uint32(10);
classifier.gbtr.def_numIters          = uint32(500);
classifier.gbtr.pos_numIters          = uint32(1):uint32(1):uint32(1000);

%  [B] Define the regression 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: 'mae'
%
%      Scoring algorithms (used in regression (estimation) process):
%           - ae    (absolute error)     
%           - mae   (mean absolute error)
%           - se    (squared error)
%           - sle   (squared log error)
%           - mse   (mean squared error)
%           - msle  (mean squared log error)
%           - rmse  (root mean squared error)
%           - rmsle (root mean squared log error)
%           - gii   (gini index)
%           - ngii  (normalized gini index)
%           - eem   (estimation error (using maximum value))
%           - eer   (estimation error (using data range))

score.pos_alg = {{'ae'    '(absolute error)'};               ...
                 {'mae'   '(mean absolute error)'};          ...
                 {'se'    '(squared error)'};                ...
                 {'sle'   '(squared log error)'};            ...
                 {'mse'   '(mean squared error)'};           ...
                 {'msle'  '(mean squared log error)'};       ...
                 {'rmse'  '(root mean squared error)'};      ...
                 {'rmsle' '(root mean squared log error)'};  ...
                 {'rss'   '(residual sum of squares)'};      ...
                 {'cod'   '(coefficient of determination)'}; ...
                 {'gii'   '(gini index)'};                   ...
                 {'ngii'  '(normalized gini index)'};        ...
                 {'eem'   '(estimation error (max. val.))'}; ...
                 {'eer'   '(estimation error (range))'};};
score.def_alg =  'mae';

%% 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;
    
    % Set the temporary saver filename
    % timestamp = Timestamp();
    % sett.tmp_saver_name = ['tmp_' timestamp];
    
    % Set the residual analysis option
    sett.residual_analysis = false;
    
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
    
    % Set the temporary saver filename
    % if (~isfield(sett, 'tmp_saver_name'))
    %     timestamp = Timestamp();
    %     sett.tmp_saver_name = ['tmp_' timestamp];
    % end
    
    % Set the residual analysis option
    if (~isfield(sett, 'residual_analysis'))
        sett.residual_analysis = false;
    else    
        if (isempty(sett.mean_residuals))
            error('sett.mean_residuals not set: residual analysis is on');
        end
        if (isempty(sett.label_original))
            error('sett.label_original not set: residual analysis is on');
        end
    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: Regression (estim) algorithm used in the feature selection
%              DEFAULT: Cart: Classification and Regression Trees 
if (~DEF_CLASSIFIER)
    STR_CLASSIFIER = sett.classifier.alg;
    STR_CLASSIFIER_STATE = ' (DEFAULT)';
    sett.classifier = ...
            check_regression_sett(sett.classifier, classifier);
else

    % check the classifier (set the regression algorithm)
    sett.classifier.alg = check_regression_alg(sett.classifier.alg);
    if (~isempty(sett.classifier.alg))
         STR_CLASSIFIER = sett.classifier.alg;
         STR_CLASSIFIER_STATE = '';
         sett.classifier = ...
            check_regression_sett(sett.classifier, classifier);
    else
        disp(' '                                               );
        disp('Unknown classifier, set to DEFAULT: cart'        );
        disp('Known classifiers:'                              );
        disp('   - cart  (Classification and Regression Trees)');
        
        STR_CLASSIFIER = DEFAULT_CLASSIFIER;
        STR_CLASSIFIER_STATE = ' (AUTOMATIC)';
        sett.classifier = ...
            check_regression_sett(sett.classifier, classifier);
    end
end

%      Handle: The ccoring method used used in the feature selection
%              DEFAULT: 'mae' (mean absolute error)
if (~DEF_SCORE)
    STR_SCORE = sett.score;
    STR_SCORE_STATE = ' (DEFAULT)';
else

    % check the regression score (set the regression score)
    sett.score = check_regression_score(sett.score);
    if (~isempty(sett.score))
         STR_SCORE = sett.score;
         STR_SCORE_STATE = '';
    else
        
        disp(' '                                           );
        disp('Unknown scoring method, set to DEFAULT: mae' );
        disp('Known scoring methods:'                      );
        disp('   - ae    (absolute error)'                 );            
        disp('   - mae   (mean absolute error)'            );
        disp('   - se    (squared error)'                  );
        disp('   - sle   (squared log error)'              );
        disp('   - mse   (mean squared error)'             );
        disp('   - msle  (mean squared log error)'         );
        disp('   - rmse  (root mean squared error)'        );
        disp('   - rmsle (root mean squared log error)'    );
        disp('   - rss   (residual sum of squares)'        );
        disp('   - cod   (coefficient of determination)'   );
        disp('   - gii   (gini index)'                     );
        disp('   - ngii  (normalized gini index)'          );

        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)'       );
        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 regression models
out.S        = [];       % Empty vector of best performing features (idx)
out.SCR      = [];       % Empty vector of regression scores (optional)
out.MAE      = [];       % Empty vector of mean absolute error (MAE)
out.MSE      = [];       % Empty vector of mean squared error (MSE)
out.RMSE     = [];       % Empty vector of root mean squared error (RMSE)
out.EER      = [];       % Empty vector of relative estimation error (EER)
out.RSS      = [];       % Empty vector of residual sum of squares (RSS)
out.COD      = [];       % Empty vector of coeff. of determination (COD)
out.RHO      = [];       % Empty vector of Spearman's correlation coeff.
out.P        = [];       % Empty vector of significance level of corr.
out.SCR_std  = [];       % Empty vector of regression scores std (optional)
out.MAE_std  = [];       % Empty vector of mean absolute error std (MAE)
out.MSE_std  = [];       % Empty vector of mean squared error std (MSE)
out.RMSE_std = [];       % Empty vector of root mean squ. err. std (RMSE)
out.EER_std  = [];       % Empty vector of relative estima. err. std (EER)
out.RSS_std  = [];       % Empty vector of residual sum of sq. std (RSS)
out.COD_std  = [];       % Empty vector of coeff. of determ. std (COD)
out.RHO_std  = [];       % Empty vector of Spearman's corr. coeff. std
out.P_std    = [];       % Empty vector of significance level of corr. std

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.reg_method = STR_CLASSIFIER;
mdl.reg_score  = STR_SCORE;
mdl.CV_folds   = STR_K_FOLD;
mdl.num_CV     = STR_NUM_CROSS_VAL;
mdl.reg_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 = [];
mae      = [];
mse      = [];
rmse     = [];
eer      = [];
rss      = [];
cod      = [];
rho      = [];
p        = [];

score_std = [];
mae_std   = [];
mse_std   = [];
rmse_std  = [];
eer_std   = [];
rss_std   = [];
cod_std   = [];
rho_std   = [];
p_std     = [];

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

save_iter    = 1;
% tmp_dirname  = ['_' sett.tmp_saver_name];
% tmp_filename = [tmp_dirname filesep sett.tmp_saver_name];

% if (~exist(tmp_dirname, 'dir'))
%     mkdir(tmp_dirname);
% end

%% Check the residual analysis
if (sett.residual_analysis)
    label_orig  = sett.label_original;
    label_resid = label_train;
    mean_resid  = sett.mean_residuals;
else
    label_orig  = label_train;
    label_resid = label_train;
    mean_resid  = 0;
end

%% 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;
    mae(length(out.S) + 1)       = 0;
    mse(length(out.S) + 1)       = 0;
    rmse(length(out.S) + 1)      = 0;
    eer(length(out.S) + 1)       = 0;
    rss(length(out.S) + 1)       = 0;
    cod(length(out.S) + 1)       = 0;
    rho(length(out.S) + 1)       = 0;
    p(length(out.S) + 1)         = 0;
    score_std(length(out.S) + 1) = 0;
    mae_std(length(out.S) + 1)   = 0;
    mse_std(length(out.S) + 1)   = 0;
    rmse_std(length(out.S) + 1)  = 0;
    eer_std(length(out.S) + 1)   = 0;
    rss_std(length(out.S) + 1)   = 0;
    cod_std(length(out.S) + 1)   = 0;
    rho_std(length(out.S) + 1)   = 0;
    p_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 regression 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_resid(train_idx);
                    test_tbl  = tbl_train(test_idx, sel);
                    test_lbl  = label_orig(test_idx);
                    
                    [~, pred_nested, reg_model] = ...
                        perf_regression(train_tbl, train_lbl, ...
                            test_tbl, cls);

                    predictions{fold} = pred_nested(:) + mean_resid;
                    presented{fold}   = test_lbl(:);
                end
                
                % Update the true and predicted labels
                pred_labels(cvi, :) = predictions;
                true_labels(cvi, :) = presented;
            end
            
            % c) Compute overall estimation (reg) performance score that 
            %    is achieved so far (default: 'mae' estim. scoring alg.) 
            MAE_cv  = zeros(ncv, fld);
            MSE_cv  = zeros(ncv, fld);
            RMSE_cv = zeros(ncv, fld);
            EER_cv  = zeros(ncv, fld);
            RSS_cv  = zeros(ncv, fld);
            COD_cv  = zeros(ncv, fld);
            RHO_cv  = zeros(ncv, fld);
            P_cv    = zeros(ncv, fld);          
            SCR_cv  = zeros(ncv, fld);

            for cvi = 1:ncv
                for fold = 1:fld
                    
                    true_l = true_labels{cvi, fold};
                    pred_l = pred_labels{cvi, fold};
                    
                    nans  = (isnan(true_l) | isnan(pred_l));
                    infs  = (isinf(true_l) | isinf(pred_l));
                    null  = (isempty(true_l) | isempty(pred_l));
                    leave = ~(nans | infs | null);
                    
                    true_l = true_l(leave);
                    pred_l = pred_l(leave);
                    
                    % Calculate the regression scores
                    SCR_cv(cvi, fold) = ...
                        calc_regression_score(pred_l, true_l, sett.score);
                    MAE_cv(cvi, fold) = ...
                        calc_regression_score(pred_l, true_l, 'mae');
                    MSE_cv(cvi, fold) = ...
                        calc_regression_score(pred_l, true_l, 'mse');
                    RMSE_cv(cvi, fold) = ...
                        calc_regression_score(pred_l, true_l, 'rmse');
                    RSS_cv(cvi, fold) = ...
                        calc_regression_score(pred_l, true_l, 'rss');
                    COD_cv(cvi, fold) = ...
                        calc_regression_score(pred_l, true_l, 'cod');
                    EER_cv(cvi, fold) = ...
                        calc_regression_score(pred_l, true_l, 'eer', ...
                            range(label_train), 1);
                    [RHO_cv(cvi, fold), P_cv(cvi, fold)] = ...
                        corr(pred_l, true_l, 'type', 'Spearman');
                end
            end
            
            % Calculate the mean scores over all cross validations
            MAE_temp  = mean(MAE_cv(:));
            MSE_temp  = mean(MSE_cv(:));
            RMSE_temp = mean(RMSE_cv(:));
            EER_temp  = mean(EER_cv(:));
            RSS_temp  = mean(RSS_cv(:));
            COD_temp  = mean(COD_cv(:));
            RHO_temp  = mean(RHO_cv(:));
            P_temp    = mean(P_cv(:));
            SCR_temp  = mean(SCR_cv(:));
            
            MAE_std   = std(MAE_cv(:));
            MSE_std   = std(MSE_cv(:));
            RMSE_std  = std(RMSE_cv(:));
            EER_std   = std(EER_cv(:));
            RSS_std   = std(RSS_cv(:));
            COD_std   = std(COD_cv(:));
            RHO_std   = std(RHO_cv(:));
            P_std     = std(P_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;
                mae(length(out.S) + 1)      = MAE_temp;
                mse(length(out.S) + 1)      = MSE_temp;
                rmse(length(out.S) + 1)     = RMSE_temp;
                eer(length(out.S) + 1)      = EER_temp;
                rss(length(out.S) + 1)      = RSS_temp;
                cod(length(out.S) + 1)      = COD_temp;
                rho(length(out.S) + 1)      = RHO_temp;
                p(length(out.S) + 1)        = P_temp;
                
                score_std(length(out.S) + 1) = SCR_std;
                mae_std(length(out.S) + 1)   = MAE_std;
                mse_std(length(out.S) + 1)   = MSE_std;
                rmse_std(length(out.S) + 1)  = RMSE_std;
                eer_std(length(out.S) + 1)   = EER_std;
                rss_std(length(out.S) + 1)   = RSS_std;
                cod_std(length(out.S) + 1)   = COD_std;
                rho_std(length(out.S) + 1)   = RHO_std;
                p_std(length(out.S) + 1)     = P_std;
                
                iterator                   = +iterator + 1;
                var.num_features(iterator) = length(out.S);
                var.score(iterator)        = SCR_temp;
                var.mae(iterator)          = MAE_temp;
                var.mse(iterator)          = MSE_temp;
                var.rmse(iterator)         = RMSE_temp;
                var.eer(iterator)          = EER_temp;
                var.rss(iterator)          = RSS_temp;
                var.cod(iterator)          = COD_temp;
                var.rho(iterator)          = RHO_temp;
                var.p(iterator)            = P_temp;
                
                var.score_std(iterator) = SCR_std;
                var.mae_std(iterator)   = MAE_std;
                var.mse_std(iterator)   = MSE_std;
                var.rmse_std(iterator)  = RMSE_std;
                var.eer_std(iterator)   = EER_std;
                var.rss_std(iterator)   = RSS_std;
                var.cod_std(iterator)   = COD_std;
                var.rho_std(iterator)   = RHO_std;
                var.p_std(iterator)     = P_std;
                
                statstr = ['SFFS(' num2str(length(out.S))   ...
                    '): ' num2str(i1) ':'                   ...
                    ' SCR = ' num2str(SCR_temp)             ...
                    ' +- ' num2str(SCR_std)                 ...
                    ' MAE = ' num2str(MAE_temp)             ...
                    ' +- ' num2str(MAE_std)                 ...
                    ' MSE = ' num2str(MSE_temp)             ...
                    ' +- ' num2str(MSE_std)                 ...
                    ' RMSE = ' num2str(RMSE_temp)           ...
                    ' +- ' num2str(RMSE_std)                ...
                    ' EER = ' num2str(EER_temp)             ...
                    ' +- ' num2str(EER_std)                 ...
                    ' RSS = ' num2str(RSS_temp)             ...
                    ' +- ' num2str(RSS_std)                 ...
                    ' COD = ' num2str(COD_temp)             ...
                    ' +- ' num2str(COD_std)                 ...
                    ' RHO = ' num2str(RHO_temp)             ...
                    ' +- ' num2str(RHO_std)                 ...
                    ' P = ' num2str(P_temp)                 ...
                    ' +- ' num2str(P_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.MAE  = [out.MAE mae(length(out.S))];
    out.MSE  = [out.MSE mse(length(out.S))];
    out.RMSE = [out.RMSE rmse(length(out.S))];
    out.EER  = [out.EER eer(length(out.S))];
    out.RSS  = [out.RSS rss(length(out.S))];
    out.COD  = [out.COD cod(length(out.S))];
    out.RHO  = [out.RHO rho(length(out.S))];
    out.P    = [out.P p(length(out.S))];
    
    out.SCR_std  = [out.SCR_std score_std(length(out.S))];
    out.MAE_std  = [out.MAE_std mae_std(length(out.S))];
    out.MSE_std  = [out.MSE_std mse_std(length(out.S))];
    out.RMSE_std = [out.RMSE_std rmse_std(length(out.S))];
    out.EER_std  = [out.EER_std eer_std(length(out.S))];
    out.RSS_std  = [out.RSS_std rss_std(length(out.S))];
    out.COD_std  = [out.COD_std cod_std(length(out.S))];
    out.RHO_std  = [out.RHO_std rho_std(length(out.S))];
    out.P_std    = [out.P_std p_std(length(out.S))];
    
    mdl.reg_models{1, length(out.S)} = reg_model;

    % save([tmp_filename '_' num2str(save_iter) '.mat'], ...
    %     'out', 'var', 'mdl');    
    save_iter = save_iter + 1;
                
    % [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 regression 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 regression  
                    train_tbl = tbl_train(train_idx, sel);
                    train_lbl = label_resid(train_idx);
                    test_tbl  = tbl_train(test_idx, sel);
                    test_lbl  = label_orig(test_idx);

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

                    predictions{fold} = pred_nested(:) + mean_resid;
                    presented{fold}   = test_lbl(:);
                end

                % Update the true and predicted labels
                pred_labels(cvi, :) = predictions;
                true_labels(cvi, :) = presented;
            end
   
            % c) Compute overall estimation (reg) performance score that 
            %    is achieved so far (default: 'mae' estim. scoring alg.)          
            SCR_cv = zeros(ncv, fld);

            for cvi = 1:ncv
                for fold = 1:fld
                    true_l = true_labels{cvi, fold};
                    pred_l = pred_labels{cvi, fold};

                    % Calculate the regression score
                    SCR_cv(cvi, fold) = ...
                        calc_regression_score(pred_l, true_l, 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 regression score. This means that 
            %    after excluding the least useful feature, the algorithm
            %    produces higher precision, 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.MAE(idx)  = [];
            out.MSE(idx)  = [];
            out.RMSE(idx) = [];
            out.EER(idx)  = [];
            out.RSS(idx)  = [];
            out.COD(idx)  = [];
            out.RHO(idx)  = [];
            out.P(idx)    = [];
            
            out.SCR_std(idx)  = [];
            out.MAE_std(idx)  = [];
            out.MSE_std(idx)  = [];
            out.RMSE_std(idx) = [];
            out.EER_std(idx)  = [];
            out.RSS_std(idx)  = [];
            out.COD_std(idx)  = [];
            out.RHO_std(idx)  = [];
            out.P_std(idx)    = [];
            
            best_val(length(out.S)) = exclusion_score;
            mdl.reg_models(idx) = [];
            
            % save([tmp_filename '_' num2str(save_iter) '.mat'], ...
            %     'out', 'var', 'mdl');
            save_iter = save_iter + 1;
            
            % 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.mae(iterator)          = mae(length(out.S));
            var.mse(iterator)          = mse(length(out.S));
            var.rmse(iterator)         = rmse(length(out.S));
            var.eer(iterator)          = eer(length(out.S));
            var.rss(iterator)          = rss(length(out.S));
            var.cod(iterator)          = cod(length(out.S));
            var.rho(iterator)          = rho(length(out.S));
            var.p(iterator)            = p(length(out.S));
            
            var.score_std(iterator) = score_std(length(out.S));
            var.mae_std(iterator)   = mae_std(length(out.S));
            var.mse_std(iterator)   = mse_std(length(out.S));
            var.rmse_std(iterator)  = rmse_std(length(out.S));
            var.eer_std(iterator)   = eer_std(length(out.S));
            var.rss_std(iterator)   = rss_std(length(out.S));
            var.cod_std(iterator)   = cod_std(length(out.S));
            var.rho_std(iterator)   = rho_std(length(out.S));
            var.p_std(iterator)     = p_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
[~, min_idx] = min(out.SCR);

out.S    = out.S(1:min_idx);
out.SCR  = out.SCR(1:min_idx);
out.MAE  = out.MAE(1:min_idx);
out.MSE  = out.MSE(1:min_idx);
out.RMSE = out.RMSE(1:min_idx);
out.EER  = out.EER(1:min_idx);
out.RSS  = out.RSS(1:min_idx);
out.COD  = out.COD(1:min_idx);
out.RHO  = out.RHO(1:min_idx);
out.P    = out.P(1:min_idx);

out.SCR_std  = out.SCR_std(1:min_idx);
out.MAE_std  = out.MAE_std(1:min_idx);
out.MSE_std  = out.MSE_std(1:min_idx);
out.RMSE_std = out.RMSE_std(1:min_idx);
out.EER_std  = out.EER_std(1:min_idx);
out.RSS_std  = out.RSS_std(1:min_idx);
out.COD_std  = out.COD_std(1:min_idx);
out.RHO_std  = out.RHO_std(1:min_idx);
out.P_std    = out.P_std(1:min_idx);

mdl.reg_models = mdl.reg_models(1:min_idx);

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