function [ranked,score] = fsrmrmr(X,y,varargin)
%FSRMRMR Importance of features (predictors) for regression using MRMR algorithm.
%   IDX=FSRMRMR(TBL,Y) ranks predictors for table TBL and response Y. The
%   ranks are learned using the Minimum Redundancy Maximum Relevance
%   algorithm.
%
%   TBL contains the predictor variables. Y can be any of the following:
%      1. An array of response values. Y can be a floating-point array,
%         single or double.
%      2. The name of a variable in TBL. This variable is used as the
%         response Y, and the remaining variables in TBL are used as
%         predictors.
%      3. A formula character vector such as 'y ~ x1 + x2 + x3' specifying
%         that the variable y is to be used as the response, and the other
%         variables in the formula are predictors. Any table variables not
%         listed in the formula are not used.
%
%   IDX=FSRMRMR(X,Y) is an alternative syntax that accepts X as an N-by-P
%   matrix of predictors with one row per observation and one column per
%   predictor. Y is the response and is a floating-point array of N
%   elements.
%
%   IDX is a 1-by-P vector for P predictors. IDX are indices of
%   columns in X ordered by importance, meaning IDX(1) is the index of
%   the most important predictor.  
%
%   [IDX,SCORES]=FSRMRMR(...) also returns predictor scores SCORES, a
%   1-by-P array for P predictors. SCORES have the same order as predictors
%   in the input data, meaning SCORES(1) is the score for the first
%   predictor. Large score indicates important predictor.
%
%   [...]=FSRMRMR(X,Y,'PARAM1',val1,'PARAM2',val2,...) specifies optional
%   parameter name/value pairs:
%       'CategoricalPredictors' - List of categorical predictors. Pass
%                        'CategoricalPredictors' as one of:
%                          * A numeric vector with indices between 1 and P,
%                            where P is the number of columns of X or
%                            variables in TBL.
%                          * A logical vector of length P, where a true
%                            entry means that the corresponding column of X
%                            or TBL is a categorical variable. 
%                          * 'all', meaning all predictors are categorical.
%                          * A string array or cell array of character
%                            vectors, where each element in the array is
%                            the name of a predictor variable. The names
%                            must match variable names in table TBL.
%                        Default: for a matrix input X, no categorical
%                        predictors; for a table TBL, predictors are
%                        treated as categorical if they are strings, cell
%                        arrays of character vectors, logical, or unordered
%                        of type 'categorical'.
%       'UseMissing'   - Logical value specifying whether missing values 
%                        in the predictors must be used or discarded. When
%                        'UseMissing' is set to true, FSRMRMR groups
%                        together missing values and uses them to compute
%                        mutual information between pairs of features and
%                        between each feature and the response variable.
%                        When 'UseMissing' is set to false, FSRMRMR omits
%                        missing values in each feature while computing the
%                        mutual information. Default: false
%       'Weights'      - Vector of observation weights, one weight per
%                        observation. FSRMRMR normalizes the weights to add
%                        up to one. Default: ones(size(X,1),1). For an
%                        input table TBL, the 'Weights' value can be the
%                        name of a variable in TBL.
%       'Verbose'      - Verbosity flag, a nonnegative integer:
%                         * 0  - FSRMRMR does not display any diagnostic
%                                messages (default).
%                         * 1  - FSRMRMR displays major diagnostic info.
%                         * >1 - FSRMRMR displays a lot of diagnostic info.
%
%   Example:
%       % Identify important predictors in the carsmall dataset:
%       load carsmall
%       X = [Cylinders Displacement Horsepower Weight];
%       [idx,scores] = fsrmrmr(X,MPG);
%       bar(scores(idx))
%       xlabel('Predictor rank')
%       ylabel('Predictor importance score')

%   Copyright 2020-2021 The MathWorks, Inc.

y = convertStringsToChars(y);

if nargin > 2
    [varargin{:}] = convertStringsToChars(varargin{:});
end

args = {'verbose' 'usemissing'};
defs = {        0        false};
[verbose,useNaNs,~,varargin] = internal.stats.parseArgs(args,defs,varargin{:});

validateattributes(verbose,{'numeric'},{'scalar' 'integer' 'nonnegative'},...
    'fsrmrmr','Verbose');
verbose = double(verbose);

useNaNs = internal.stats.parseOnOff(useNaNs,'UseMissing');

[X,y,w,dataSummary] = ...
  classreg.learning.regr.FullRegressionModel.prepareData(...
  X,y,varargin{:},'OrdinalIsCategorical',false);

% Accept only dense floating-point.
internal.stats.checkSupportedNumeric('X',X);

D = size(X,2);

% Get indices of categorical features
catpreds = dataSummary.CategoricalPredictors;
iscat = false(1,D);
iscat(catpreds) = true;

% Use at most 256 bins to bin continuous features. These will be lumped
% into larger bins by adaptive binning.
nbin = 256*ones(D,1);

tstart = tic;

% Bin response
ybinned = 1 + classreg.learning.treeutils.binPredictors(y,256);

% binPredictors returns -1 for missing values; these get transformed into
% zeros.
if ~any(iscat)
    Xbinned = 1 + classreg.learning.treeutils.binPredictors(X,nbin);
else
    Xbinned = zeros(size(X),'int32');
    Xbinned(:,~iscat) = 1 + classreg.learning.treeutils.binPredictors(...
        X(:,~iscat),nbin(~iscat));
    
    Xbinned(:,iscat) = classreg.learning.fsutils.indexCategoricals(X(:,iscat));
        
    if any(max(Xbinned(:,iscat),[],1) > intmax('int32')-1)
        error(message('stats:fscmrmr:TooManyCategories'))
    end
end

% Compute MI across X and between X and Y
[Ix,Iy,ranked,score,predsWithPositiveMI] = ...
    classreg.learning.fsutils.miXy(Xbinned,ybinned,w,iscat,false,useNaNs,verbose);

% Edge case: MI between all features and response is too small
if ~any(predsWithPositiveMI)
    return
end

telapsed = toc(tstart);

if verbose>0
    fprintf('Mutual information computed in %10.5e s.\n',telapsed)
end

tstart = tic;

% Convert logical indices to integer indices
predsWithPositiveMI = find(predsWithPositiveMI);

% Find ranks and MIQ ratios on the subset of features with positive MI
% values
[rankedPos,scorePos] = classreg.learning.fsutils.mrmrMIQsort(Iy,Ix);

telapsed = toc(tstart);

if verbose>0
    fprintf('Predictors are ranked in       %10.5e s.\n',telapsed)
end

% Index into ranks and scores for all features
Dgood = numel(predsWithPositiveMI);
ranked(1:Dgood)        = predsWithPositiveMI(rankedPos);
score(ranked(1:Dgood)) = scorePos;

end

% LocalWords:  usemissing
