function [metafeatures, Weights, GSorted, GSortedInd] = metafeatures(data,varargin)
%METAFEATURES - Attractor metagene learning.
% 
% M = metafeatures(X) implements attractor metagene learning, an iterative
% algorithm described in [1]. The algorithm creates attractor metafeatures
% that are weighted sums of the features in X. Weights are calculated by a
% distance like function based on mutual information. X is P-by-N numeric
% matrix. Columns of X correspond to different samples, like different
% patients or time points. Rows of X correspond to variables, such as
% measurements of gene expression. M, is a matrix of metafeatures
% calculated by the algorithm.
%
% M = metafeatures(T) takes T, a P-by-N numeric table.
% 
% [M, W] = metafeatures(...) returns W, a P-by-r matrix of the metafeature
% weights. r is the number of replicates performed by the algorithm as
% controlled by the 'Start' or 'Replicate' parameter-value-pairs described
% below. By default r is 1. If X is supplied then, M = W'*X, otherwise M =
% W'*T{:,:}
% 
% [M, W, GSorted] = metafeatures(X,G,...) returns Gsorted, a P-by-r cell
% string of variable names sorted by decreasing weight. The ith column of
% Gsorted lists the variable names in order of their contribution to ith
% metafeature. G is an optional input. If G is empty, or if G is not
% provided and Gsorted is requested, then metafeatures will automatically
% feature names as 'Var_i' associated with the ith row of X.
% 
% [M, W, GSorted] = metafeatures(T,...) returns Gsorted, a P-by-r cell
% array of variable names sorted by decreasing weight. The variable names
% are taken from the row names of T. The ith row of Gsorted lists the
% variable names in order of their contribution to ith metafeature.
% 
% [M, W, Gsorted, GsortedInd] = metafeatures(...) returns GsortedInd,
% indices such that GSorted = G(GsortedInd).
% 
% [...] = metafeatures(...,
% 'Parameter1','Value1','Parameter2','Value2',...) enables you to specify
% optional parameter value pairs that control elements of the algorithm.
% Valid parameter string are:
%
% 'Alpha' - A real number greater than 0. This parameter controls the
% non-linearity of the function used to calculated the weights. Concretely,
% it controls the number of metafeatures that exist for a dataset. As alpha
% increases the number of metafeatures will tend to increase. This is often
% the most important parameter to adjust in the analysis of a dataset.
% (default: 5)
%
% 'Start' - Initial weights are chosen between 0 and 1. 'Start' determines
%           the method for choosing the initial weights in the algorithm.
%           Options are:
%     
%           'random' - Initialize the weights to a vector of positive
%                      weights chosen uniformly at random and scaled such
%                      that they sum to 1. Choose a different initial
%                      weight vector for each replicate. (default)
%
%           'robust' - If X or T has N columns then run the algorithm N
%                      times. On the ith evaluation of the algorithm, the
%                      weights will be initialized to all zeros with the
%                      exception of the ith weight, which will be set to 1.
%                      This option is useful when attempting to find all
%                      (or nearly all) metafeatures of a dataset.
%
%           Matrix - N-by-r matrix of initial weights. The algorithm will
%                    be run r times. The weights in the ith run of the
%                    algorithm will be initialized to the ith column of the
%                    matrix.
%
% 'Replicates' - Number of times to repeat the algorithm. This option is
% only valid with the 'random' start option. (default: 1)
%
% 'ReturnUnique' - A logical value. If true, then only the unique
% metafeatures are returned. This is useful when the algorithm is evaluated
% more than once either by setting the replicates parameter, setting robust
% start, or specifying a Matrix start with more than 1 row, since the same
% metafeature will be discovered for many different initializations. A
% metafeature is considered unique if the correlation between it and all
% previously found metafeatures is less than 'UnqiueTolerance'. (default:
% true)
%
% 'UniqueTolerance' - A real number between 0 and 1 (default: .98)
%
% 'Options' - A structure specifying options for the iterative algorithm
%             used to find metafeatures. Fields of the structure are:
%
%     Display - Level of display to output. Choices are 'off' or 'iter'.
%     (default: off)
%
%     MaxIter - Maximum number of iterations allowed. (default: 100)
%
%     Tolerance - If M changes by less than Tolerance in an iteration, then
%     stop. (default: 1e-6)
%
%     Streams - A RandStream object. If you do not specify Streams,
%     metafeauture will use the default random stream.
%
%     UseParallel - Logical value indicating whether to perform
%     calculations in parallel if parallel pool is available. For problems
%     with large datasets relative to the available system memory,
%     parallelization may degrade performance. (default: false)
%
% See also: RELIEFF, SEQUENTIALFS, RANDFEATURES, RANKFEATURES

% References:
%
% [1] Cheng, Wei-Yi, Tai-Hsien Ou Yang, and Dimitris Anastassiou.
%     "Biomolecular events in cancer revealed by attractor metagenes." PLoS
%     computational biology 9.2 (2013): e1002920.
%

% Copyright MathWorks 2014 

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

if nargin < 1
    error(message('bioinfo:metafeatures:IncorrectNumberOfArguments'));
end

%% validate input, and convert it to a common form

% argument names and defaults
pnames = {'alpha', 'replicates', 'returnunique','uniquetolerance','options', 'start'};
dflts =  {      5,           [],           true,              .98,       [],'random'};

% determine if data is table or numeric
if ~(istable(data) || (ismatrix(data) && isnumeric(data)))
    error(message('bioinfo:metafeatures:InvalidData'));
end

%% Transpose into obs/features orientation

[numFeatures,numSamples] = size(data);

if istable(data)
    dataNames = data.Properties.RowNames;
    if isempty(dataNames)
        dataNames = fastUniqueNames(numFeatures);
    end
    data = data{:,:}'; % unpack and transpose
    
    if nargin>1
        % if PVPs supplied check to provide more informative error message
        if ~(ischar(varargin{1}) && ismember(lower(varargin{1}),pnames))
            error(message('bioinfo:metafeatures:TableDataBadPVP'))
        end
    end
else
    data = data'; % unpack and transpose
    
    % numeric case, check if data names were provided
    if nargin == 1
        % variable names unset
        dataNames = fastUniqueNames(numFeatures);
    elseif isempty(varargin{1})
        % variable names set to empty explicitly
        dataNames = fastUniqueNames(numFeatures);

        varargin(1) = []; % so varargin can be passed for further validation
    elseif  iscellstr(varargin{1})
        % variable names supplied
        if ~all(size(varargin{1}) == [size(data,2), 1]) % note, tranpose already occurred
            error(message('bioinfo:metafeatures:VariableNameSize'));
        end
        dataNames = varargin{1};
        varargin(1) = []; % so we can pass varargin to further validation
    else
        % not cell string, not empty, must be parameter name, check to give
        % more informative error message for this case
        if ~(ischar(varargin{1}) && any(strncmpi(varargin{1},pnames,length(varargin{1}))))
            error(message('bioinfo:metafeatures:NumericDataBadPVP'))
        end
        
        dataNames = fastUniqueNames(numFeatures);

    end
end

% Data is a 2d numeric array, dataNames is a cellstr of appropriate size.
% For illustration or debugging the following asserts will all pass if
% uncommented.
% assert(ismatrix(data))
% assert(isnumeric(data))
% assert(iscellstr(dataNames))
% assert(size(dataNames,1) == size(data,2));

% parse the parameter value pairs, then validate values

[alpha, replicates, returnUnique, uniqueTolerance, options, start] ...
    = internal.stats.parseArgs(pnames, dflts, varargin{:});

pvpParser = inputParser;

pvpParser.addRequired('alpha'          ,@validateAlpha)
pvpParser.addRequired('replicates'     ,@validateReplicates)
pvpParser.addRequired('returnUnique'   ,@validateReturnUnique)
pvpParser.addRequired('uniqueTolerance',@validateUniqueTolerance)
pvpParser.addRequired('options'        ,@validateOptions)
pvpParser.addRequired('start'          ,@(start) validateStart(start,numFeatures))

pvpParser.parse(alpha,replicates,returnUnique,uniqueTolerance,options,start);
pvpStruct = pvpParser.Results;

[display, maxIter, tol, useParallel, streams] = expandOptions(pvpStruct.options);

if useParallel
    poolSize = internal.stats.parallel.getParallelPoolSize('Guarded',true);
    if poolSize == 0
        useParallel = false;
        warning(message('bioinfo:metafeatures:PoolSizeZero'))
    end
end
% other parameter defaults not exposed to user
bins    = 6;
sporder = 3;

% Tranpose start into obs/features orientation if it's numeric
if isnumeric(pvpStruct.start)
    pvpStruct.start = pvpStruct.start';
    start = pvpStruct.start;
end

%%
if ~isempty(replicates) && (~strcmpi('random',pvpStruct.start))
    error(message('bioinfo:metafeatures:InvalidReplicatesWithStart'));
end
numStarts = [];
if ischar(pvpStruct.start) && strcmp(pvpStruct.start,{'robust'})
    numStarts = numFeatures;
end
if isnumeric(pvpStruct.start)
    numStarts = size(pvpStruct.start,1);
end
if isempty(numStarts) && ~isempty(pvpStruct.replicates)
    numStarts = pvpStruct.replicates;
elseif isempty(numStarts)
    numStarts = 1;
end

%% Caching self information

if display
    disp('Caching self information ...')
    selfInfoTiming = tic;
end

dataInfoMax = zeros(1,numFeatures);

if useParallel
    parfor ind = 1:numFeatures
        dataInfoMax(ind) = bioinfo.internal.splineMI(data(:,ind),data(:,ind),bins,sporder);
    end
else
    for ind = 1:numFeatures
        dataInfoMax(ind) = bioinfo.internal.splineMI(data(:,ind),data(:,ind),bins,sporder);
    end
end

if display
    disp(['... done. Took ' num2str(toc(selfInfoTiming)) ' seconds.']);
end

%% Caching entropy and binning information

if display
    disp('Caching entropy and binning information...')
    cachingTime = tic;
end

[~, miFunc] = bioinfo.internal.splineMI([],data,bins,sporder);

if display
    disp(['... done. Took ' num2str(toc(cachingTime)) ' seconds.']);
end

%% Run metagene attractor learning

metafeatures = zeros(numSamples,numStarts);
Weights   = zeros(numStarts,numFeatures);
GSorted     = cell(numStarts,numFeatures);
GSortedInd  = zeros(numStarts,numFeatures);

% setup dispfmt independent of if we're using it, to avoid an undefined
% variable error from the parfor.
if useParallel
    dispfmt = '%5d\t%5d\t  %1.2e\t    %4d\n';
    if display
        fprintf('\t\t\t\tnon-zero\n');
        fprintf('Worker \t iter\t      diff\t weights \n')
    end
else
    dispfmt = '%5d\t%5d\t  %1.2e\t    %4d\n';
    if display
        fprintf('\t\t\t\tnon-zero\n');
        fprintf('Found \t iter\t      diff\t weights \n')
    end
end

if useParallel
    parfor geneIter = 1:numStarts
        [metafeatures(:,geneIter), Weights(geneIter,:), GSorted(geneIter,:),GSortedInd(geneIter,:)] ...
            = coreAlgorithm(data,dataNames,geneIter,start,useParallel,streams,numFeatures,numSamples,tol,maxIter,alpha,display,dispfmt,miFunc,bins,sporder,dataInfoMax);
    end
else
    for geneIter = 1:numStarts
        [metafeatures(:,geneIter), Weights(geneIter,:), GSorted(geneIter,:),GSortedInd(geneIter,:)] ...
            = coreAlgorithm(data,dataNames,geneIter,start,useParallel,streams,numFeatures,numSamples,tol,maxIter,alpha,display,dispfmt,miFunc,bins,sporder,dataInfoMax);
    end
end

%% Data post processing

if pvpStruct.returnUnique
    [metafeatures, Weights, GSorted, GSortedInd] = ...
        getUniqueMetagene(metafeatures,Weights,GSorted,GSortedInd,pvpStruct.uniqueTolerance);
end

% transpose outputs back to features/obs orientation
metafeatures = metafeatures';
Weights   = Weights';
GSorted   = GSorted';
GSortedInd   = GSortedInd';

end

% core loop
function [thisMetagene, weights, dataNames, geneIndex] = coreAlgorithm(data,dataNames,geneIter,start,useParallel,streams,numFeatures,numSamples,tol,maxIter,alpha,display,dispfmt,miFunc,bins,sporder,dataInfoMax)
if useParallel
    t = getCurrentTask();
    workerID = t.ID;
end
% make start options have effect
if isnumeric(start)
    weights = start(geneIter,:);
elseif strcmpi('robust',start)
    weights = zeros(1,numFeatures);
    weights(geneIter) = 1;
else
    % random case
    if useParallel && iscell(streams)
        weights = rand(streams{labindex},1,numFeatures);
    else
        weights = rand(streams,1,numFeatures);
    end
end

weights = weights/sum(weights);

thisMetagene = sum(data.*(repmat(weights,numSamples,1)),2);

first = true;
iter  = 1;
oldMetagene = thisMetagene;
diff = norm(oldMetagene-thisMetagene);
block_inds = zeros(size(dataNames));
while   (diff>tol || first) && iter<=maxIter;
    
    oldMetagene = thisMetagene;
    
    metanorm = bioinfo.internal.splineMI(thisMetagene,thisMetagene,bins,sporder); % false, don't use parallel in this case
    
    weights = ...
        miFunc(thisMetagene)./max(dataInfoMax,metanorm);
    
    if first || diff>1e-3
        % only calculate disqualifications on basis of negative
        % correlation if we're not already somewhat converged
        block_inds = corr(thisMetagene, data) < 0;
    end
    weights(block_inds) = 0;
    
    weights = weights.^alpha;
    weights = weights/sum(weights);
    
    thisMetagene = sum(data.*(repmat(weights,numSamples,1)),2);
    
    diff = norm(oldMetagene-thisMetagene);
    
    if display
        if useParallel
            fprintf(dispfmt,workerID,iter,diff,sum(block_inds));
        else
            fprintf(dispfmt,geneIter,iter,diff,sum(block_inds));
        end
    end
    
    first = false;
    iter = iter+1;
end
if iter >= maxIter
    warning(message('bioinfo:metafeatures:EarlyTerm'));
end

[~, geneIndex] = sort(weights,2,'descend');
dataNames = dataNames(geneIndex); % in place reorder
end


% validation functions

function bool = validateAlpha(alpha)
bool = true;

if ~( isnumeric(alpha) && isscalar(alpha) && alpha > 0)
    error(message('bioinfo:metafeatures:InvalidAlpha'));
end
end

function bool = validateReplicates(reps)
bool = true;
if ~(isempty(reps) || (isnumeric(reps) && isscalar(reps) && reps>0 && floor(reps) == reps))
    error(message('bioinfo:metafeatures:InvalidReplicates'));
end
end

function bool = validateReturnUnique(returnUnique)
bool = true;
if ~( islogical(returnUnique) && isscalar(returnUnique))
    error(message('bioinfo:metafeatures:InvalidReturnUnique'));
end
end

function bool = validateUniqueTolerance(tol)
bool = true;
if ~( isnumeric(tol) && isscalar(tol) && tol>0 && tol<1)
    error(message('bioinfo:metafeatures:InvalidUniqueTolerance'));
end
end

function bool = validateOptions(opts)
bool = true;
if ~( isempty(opts) || isstruct(opts))
    error(message('bioinfo:metafeatures:InvalidOptions'));
end
end

function bool = validateStart(start,p)
bool = true;
if ~( isnumeric(start) || (ischar(start) && ismember(start,{'random','robust'})))
    error(message('bioinfo:metafeatures:InvalidStart'));
end

if isnumeric(start)
    if ~(ismatrix(start) && size(start,1) == p)
        error(message('bioinfo:metafeatures:InvalidNumericalStart'));
    end
end
end

function [display, maxIter, tol, useParallel,streams] = expandOptions(opts)
defaultOpts = struct('display','off','maxiter',100,'tolerance',1e-6,'useparallel',false,'streams',RandStream.getGlobalStream);
if isempty(opts)
    opts = defaultOpts;
end

% lower names of structure's fields
names    = cellfun(@lower,fieldnames(opts),'UniformOutput',false);
optsCell = struct2cell(opts);
opts     = cell2struct(optsCell,names);

namesDefault = fieldnames(defaultOpts);

for iter = 1:length(namesDefault)
    if ~isfield(opts,namesDefault{iter})
        % value wasn't supplied, set to default
        opts.(namesDefault{iter}) = defaultOpts.(namesDefault{iter});
    end
end

% all fields are defined, now unpack
if ischar(opts.display) && any(strcmpi(opts.display,{'iter','off'}))
    display = strcmpi(opts.display,'iter'); % value is 'iter' or 'off'
else
    error(message('bioinfo:metafeatures:InvalidDisplay'))
end
maxIter = opts.maxiter;
tol     = opts.tolerance;
useParallel = opts.useparallel;
streams = opts.streams;

% validate options
if ~isscalar(maxIter) || ~isnumeric(maxIter) || maxIter<1 || ~isfinite(maxIter)
    error(message('bioinfo:metafeatures:InvalidMaxIter'))
end

if ~isscalar(tol) || ~isnumeric(tol) || tol<=0
    error(message('bioinfo:metafeatures:InvalidTolerance'))
end

if ~islogical(useParallel)
    error(message('bioinfo:metafeatures:InvalidParallel'))
end

if ~isa(streams,'RandStream') || (iscell(streams) && all(cellfun(@(x) isa(x,'RandStream'),streams)))
    error(message('bioinfo:metafeatures:InvalidStreams'))
end

try
	pool = gcp('nocreate');
catch
   % gcp() will error if the Parallel Computing Toolbox is not
   %       installed, or if it is unable to check out a license	
   pool = []; % this will make it fail next when streams are given, otherwise we do not care
end

if iscell(streams) && (isempty(pool) || ~useParallel)
    error(message('bioinfo:metafeatures:CellStreamsNoPool'));
end

if iscell(streams) && (length(streams) ~= pool.NumWorkers) 
    % if the streams are a cell, we've not failed the previous check
    % so, there is a non-empty pool and the if statement is safe
    error(message('bioinfo:metafeatures:CellStreamsWrongLength'));
end
end

function [metaUnique, weightUnique, geneListUnique,geneIndUnique] = getUniqueMetagene(metafeatures,weights,gSorted,gSortedInd,tolerance)

[~, numCols] = size(metafeatures);

metaUnique     = zeros(size(metafeatures));
weightUnique   = zeros(size(weights));
geneListUnique = cell(size(gSorted));
geneIndUnique  = zeros(size(gSortedInd));

jter = 1;
iter = 1;

metaUnique(:,jter)     = metafeatures(:,iter);
weightUnique(jter,:)   = weights(iter,:);
geneListUnique(jter,:) = gSorted(iter,:);
geneIndUnique(jter,:)  = gSortedInd(iter,:);
jter = 2;

for iter = 2:numCols
    tmpCorr = corr(metafeatures(:,iter),metaUnique(:,1:(jter-1)));
    if all(tmpCorr < tolerance)
        metaUnique(:,jter)     = metafeatures(:,iter);
        weightUnique(jter,:)   = weights(iter,:);
        geneListUnique(jter,:) = gSorted(iter,:);
        geneIndUnique(jter,:)  = gSortedInd(iter,:);
        jter = jter + 1;
    end
    
end

metaUnique(:,jter:end)     = [];
weightUnique(jter:end,:)   = [];
geneListUnique(jter:end,:) = [];
geneIndUnique(jter:end,:) = [];

end

function uniqueNames = fastUniqueNames(numFeatures)
numcell = cellfun(@num2str,num2cell(1:numFeatures,1),'UniformOutput',false)';
uniqueNames = strcat('Var',numcell);
end



