function [yhat,B,I,K] =  CRDA0(Xt,X,y,covtype,K,q,prior,B,mu,is_centered,Ind)
% CRDA0 is the basic utility function needed by the CRDA function. For 
% more details of the input parameters of CRDA0, type doc CRDA. 
%
% [yhat,B,I,K] =  CRDA0(Xt,X,y,covtype,K,q,prior,B,mu,is_centered,Ind)
%
% Output Arguments
% ----------------
% yhat     vector of size N x 1 of integers from 1, ..., G. 
%          specifies  the group to which each column of Xt has been 
%          assigned
%
% B        matrix of size p X G
%          coefficient matrix calculated by the function (or the one given 
%          in the optional 8th input argument)
%
% I        indices of the length of the row vectors of B organized in 
%          descending order, where the length is determined by optional 
%          argument 'q' (if q is scalar, q >=1, then the length is
%          determined by L_q-norm)
% 
% K        the value of K found by CRDA_coefmat or the value given as 
%          optional input to the function)
%
% AUTHORS
% -------
% Esa Ollila and Muhammad Naveed Tabassum, Aalto University, October 2019.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p = size(X,1);
n  = histcounts(y);	% class labels counts
G = length(n);

if nargin<4 
    covtype = 'ell2';
end

if nargin < 5
   % Default = use mean value of the row 'lengths' of B
    K = [];
end

if nargin<6 
    q = inf;
end

if nargin<7
    prior = 1/G * ones(1,G); % equal prior by default
end


if nargin < 8
    % Default = we need need to calcualte the coefficient matrix
    B = [];
end

if nargin < 9
    % Default: calculated the class means (mu = []);
    mu = [];
end

if nargin < 10 || isempty(is_centered)
    % Default: tranining data X is not centered w.r.t. class means
    is_centered = false; 
end

if nargin < 11
    % Default: We do not have indices of the row lengths of B
    Ind = [];
end

if isempty(mu)
    % need to calculate the class means 
    I = eye(max(y));    
    Y = I(y, :);	% nxG indicator matrix => Yij = 1, i belong to group j
    mu = X*Y ./ n; 
end

% We need to center the training data before computing the covariance
% matrix estimator 
if ~is_centered 
    for ii=1:G
        X(:,y==ii) = X(:,y==ii) - repmat(mu(:,ii),1,n(ii)); 
    end         
end

if isempty(B) 
    B = CRDA_coefmat(X,y,covtype,mu);
end

[Best,K,I] = hard_threshold(B,K,q,Ind);    
%% Predict the labels of test data
const = -0.5*sum(mu.*Best) + log(prior); % 1xK vector of constant in LDA disc. fnc
dlda = Xt'*Best + const;% scores
[~,yhat] = max(dlda,[],2);




