function out = fsMRMR(X, Y, param)
% function out = fsMRMR(X, Y, param)
% 
% MIQ scheme according to MRMR
% X - the data, each row is an instance  
% Y - the class label in  form of 1 2 3 4 5 6 7 ...
% param.k - number of selected features
% param.pool - the number features will be considered in the second iteration 
% param.type = 1: MIQ, -1:MID 

if nargin < 3
    param.k = size(X,2); param.pool = 1000; param.type = 1;
end

if ~isfield(param, 'k')
    param.k = size(X,2);
end

if ~isfield(param, 'pool')
    param.pool = 1000;
end

if ~isfield(param, 'type')
    param.type = 1;
end

nF = size(X,2); fea = zeros(param.k,1);

if ~isinteger(X)
    t = weka.filters.supervised.attribute.Discretize();
    t.setOptions(wekaArgumentString({'-R','first-last'}));
    B = SY2MY(Y);
    cat = wekaCategoricalData(X, B);
    t.setInputFormat(cat);
    clear cat;
    cat = wekaCategoricalData(X, B);
    dat = weka.filters.Filter.useFilter(cat,t);
    clear cat;
    for i = 0:nF-1
        X(:,i+1) = dat.attributeToDoubleArray(i);
    end
end

% information gain
t = zeros(nF,1);
for i=1:nF 
   line = X(:,i);
   t(i) = mutualinfo(line, Y);
end

% add features provide lowest redundancy
[orMI, idxs] = sort(-t); fea(1) = idxs(1); KMAX = min(param.pool,nF); 
idxleft = idxs(2:KMAX); mi_array = zeros(max(idxleft), param.k);
curlastfea = 1;

for ft = 2:param.k
   ncand = length(idxleft); c_mi = zeros(ncand,1);
   for i=1:ncand
      mi_array(idxleft(i),curlastfea) = mutualinfo(X(:,fea(curlastfea)), X(:,idxleft(i)));
      c_mi(i) = mean(mi_array(idxleft(i), :)); 
   end
   
   switch (param.type)
       case 1
           [tmp, tmpidx] = max(orMI(idxleft) ./ (c_mi + 0.01));
       case -1
           [tmp, fea(ft)] = max(t_mi(1:ncand) - c_mi(1:ncand));
   end
   
   fea(ft) = idxleft(tmpidx); idxleft(tmpidx) = []; curlastfea = curlastfea + 1;
end

out.fList = fea; out.w = zeros(nF,1); out.W(fea) = param.k:-1:1; out.prf = -1;